ChemShell tutorial: MD simulation using the QM/MM/GSBP method

MD simulation using the QM/MM/GSBP method

In this example of an MD simulation with the QM/MM/GSBP method, the system is prepared classically using the CHARMM programm as explained in the tutorial for QM/MM calculations of enzymes. Therefore, only those aspects of the system preparation that are specific to the application of the GSBP are explained in detail. Much of the theoretical background is mentioned only briefly. A detailed description is provided in J. Chem. Theory Comput. 2008, 9, 1600-1609 and J. Chem. Phys. 2001, 114, 2924-2937.

This chapter is divided in five parts:

System preparation with CHARMM
Definition of the inner and outer region
Definition of the dielectric interface
System truncation
MD simulation with QM/MM/GSBP

To run this example, a number of files are required, which are contained in gsbp.tgz.

1. System preparation with CHARMM

The system used for this example application is the chorismate mutase enzyme which was equilibrated and solvated previously with the CHARMM programm.

In the SMBP, the solvent molecules outside the inner region are represented by a polarisable dielectric continuum (PDC). Therefore, explicit water molecules must be deleted from the system before the QM/MM/SMBP calculations. Since CHARMM was used to prepare the system, also this step is performed with CHARMM. In this example, all water molecules with the oxygen atom more than 17 Å away atom number 5676 are assigend to the outer solvent region and deleted. This step is performed with the gsbp_prepare.inp script.

define outwater sele ( resname TIP3 ) .and. .not. .byres. ( type OH2 .and. ( ( bynu 5676 ) .around. 17.0)) end

dele atom sele outwater end

OPEN WRIT UNIT 41 CARD NAME @{jobname}.pdb
WRIT COOR pdb UNIT 41
close unit 41

open write unit 42 card name @{jobname}.psf
write  psf card unit 42
close unit 42

Before the MD, the full system is truncated to the inner region by means of a ChemShell procedure. This procedure, however, cannot update the .psf file which is necessary to define the force field in the DL_POLY module. Therefore, the system is initially truncated with CHARMM to create the .psf file of the truncated system. This step is performed with the gsbp_cut.inp script.

define outer sele all .and. .not. ( ( bynu 5676 ) .around. 19.0) end

dele atom sele outer end

OPEN WRIT UNIT 41 CARD NAME @{jobname}_in.pdb
WRIT COOR pdb UNIT 41
close unit 41

open write unit 42 card name @{jobname}_in.psf
write  psf card unit 42
close unit 42

2. Definition of the inner and outer region

The most important step of the QM/MM/GSBP setup is the definition of the inner and outer region of the macromolecule. All atoms in the outer region will be represented by the GSBP and only those atoms in the inner region are simulated explicitly. The selection of the inner region is closely related to the selection of the active atoms. All atoms that are needed to reproduce the flexibility of the system have to be incorporated into the active region. Similarly, all atoms that need to be simulated atomistically to reproduce the electrostatic potential have to be incorporated into the inner region. Typically, both active region and inner region have a spherical shape with a radius of 15 to 20 Å and enclose the active site of the enzyme. The electrostatic potential is very complex at the boundary separating inner and outer region , and therefore, not represented accurately by the GSBP. Hence, the explicit atoms in direct vicinity of the boundary should be frozen, i.e., the radius of the active region should be 2 to 3 Å smaller than that of the inner region.

	active region	inner region
selection criterion	cature flexibility around active site	capture electrostatic potential at the active site
shape	spherical	spherical
radius	15-20 Å	2-3 Å larger than the radius of the active region

In this example, all atoms within 19 Å of the C1 atoms of chorismate (atom number 5676) are assigned to the inner region, and all atoms within 16 Å of C1 belong to the active region. The do_selection.chm script performs all atom selections.

...

charmm_acquire_selection selec_var=in_at selection= \
"( ( bynu 5676 ) .around. 19.0)"


charmm_acquire_selection selec_var=act selection= \
" ( ( bynu 5676 ) .around. 16.0) "

...

set fil [open "inner_at" {RDWR CREAT TRUNC}]
puts $fil "set inner_atom_numbers [ list $in_at ]"
close $fil

...

set fil [open "active" {RDWR CREAT TRUNC}]
puts $fil "set active_atom_numbers [ list $act ]"
close $fil

...

3. Definition of the dielectric interface

Since the outer solvent molecules are represented by a PDC, the interface of macromolecule and solvent region has to be defined. In the GSBP, the dielectric interface is constructed as the envelope of the atomic van der Waals radii. A TCL list of van der Waals radii from the CHARMM22 force field can be obtained via the pb_types2vdwradii procedure which reads in a list of CHARMM atom types and returns a list of the respective van der Waals radii. Since all atoms in the inner region are modeled explicitly, the dielectric constant in the inner region is set to 1.0 corresponding to vacuum conditions.

...

set rad [ pb_types2vdwradii atom_types= $types ]

...

Moreover, one of the basic assumptions for the construction of the GSBP is the assumption that the dielectric interface is constant. It is indispensable to avoid any active explicit atom touching the dielectric interface with its van der Waals sphere. Thus, the dielectric cavity of the inner region is extended to ensure a smooth dielectric interface during the calculation. In this example, the extended dielectric cavity radius is set to 20.5 Å.

...

inner_radius= [ expr 20.5 / 0.5292 ] \

...

The code snippets used in the definition of the dielectric boundary are part of the input script for an MD simulation with the QM/MM/GSBP method (qmmmgsbp.chm) that is presented in the fifth part.

4. System truncation

To fully exploit the computational savings offered by the GSBP, it is necessary to truncate the system to the inner region. The gsbp.prepare procedure truncates the system and updates the TCL lists that define the atom types, charges, charge groups, etc.

gsbp.prepare coords= ${sys_name_id}.c  \
        groups = $groups \
        atom_charges= $charges \
        qm_region= $qmatoms \
        inner_region= $inner_atom_numbers \
        outer_region= $outer_atom_numbers \
        frozen= $frozen_atom_numbers \
        atom_types= $types \
        active= $active_atom_numbers \
        tot_coords= tcoords

In this example, the non-truncated system is saved in the tcoords fragment. The truncated system is stored in the ${sys_name_id}.c fragment object. The list of charges of the full system is always automatically saved in the variable tot_charges. The truncated lists of charges, atom types, etc. are stored in the lists with variable names prep_charges, prep_types, etc.

Again, the code snippets used to truncate the system are part of the MD input script (qmmmgsbp.chm).

5. MD simulation with QM/MM/GSBP

Now, after preparing the system in the previous steps, the QM/MM/GSBP method can be applied for MD simulation. The full input script is given in the file qmmmgsbp.chm. The GSBP is implemented in ChemShell only for application with QM/MM with mndo or pointcharge as qm_theory. In this example, the mndo program is used as qm_theory.

First, the TCL list defining the different region in the enzyme are sourced.

source ./qmatoms
source ./active
source ./frozen
source ./inner_at
source ./outer_at

Then, the centre of the inner region is defined,

set tmp [ get_atom_entry coords= ${sys_name_id}.c atom_number= 5676 ]
set xc [ lindex $tmp 1 ]
set yc [ lindex $tmp 2 ]
set zc [ lindex $tmp 3 ]
set centre [ list $xc $yc $zc ]

and the system is prepared for the GSBP (see above).

gsbp.prepare coords= ${sys_name_id}.c  \
	...

Finally, the dynamics driver is initialised with the GSBP by choosing gsbp= yes.