Geometry optimisation using the QM/MM/SMBP method
In this example of an application of the QM/MM/SMBP 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 SMBP are explained in detail. Much of the theoretical background is mentioned only briefly. A detailed description is provided in J. Chem. Theory Comput. 2009, 5, 3114-3128.
This chapter is divided in four parts:
- System preparation with CHARMM
- Definition of the inner and outer region
- Definition of the dielectric interface
- Geometry optimisation with QM/MM/SMBP
To run this example, a number of files are required, which are contained in smbp.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 outside a radius of 19 Å around atom number 5676 are assigend to the outer region and deleted. This step is performed with the smbp_prep.inp script.
...
define outwater sele (resname TIP3 ) .and. .not. .byres. ( ( bynu 5676 ) .around. 19.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
2. Definition of the inner and outer region
The most important step of the QM/MM/SMBP setup is the definition of the inner and outer region of the macromolecule. All atoms in the outer region will be represented by the SMBP 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 SMBP. 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 SMBP, 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 SMBP 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 a geometry optimisation with the QM/MM/SMBP method (qmmmsmbp.chm) that is presented in the next part.
4. Geometry optimisation with QM/MM/SMBP
Now, after preparing the system in the previous steps, the QM/MM/SMBP method can be applied to a geometry optimisation. The full input script is given in the file qmmmsmbp.chm.
First, the TCL list defining the different region in the enzyme are sourced.
source ./qmatoms source ./active source ./frozen source ./inner_at source ./outer_atThen, 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 ]
Finally, the DL-FIND optimiser is called and the SMBP is invoked by the smbp= yes
keyword:
dl-find coords= ${sys_name_id}.c \
...
theory=hybrid : [ list \
coupling= $embedding_scheme \
...
smbp= yes : [ list \
inner_region= $inner_atom_numbers \
outer_region= $outer_atom_numbers \
centre= $centre \
inner_radius= [ expr 20.5 / 0.5292 ] \
vdw_radii= $rad \
grid_size= [ expr 1.25 / 0.5292 ] \
inner_grid_size= [ expr 0.6 / 0.5292 ] \
qm_charge= $qm_ch \
nrpc= 90 \
radius_extension= 10.0 \
spline_order= 3 \
eps_mol= 1.0 \
eps_solv= 80.0 ] \
qm_theory= $qm_theory : [ list \
...
] \
mm_theory= dl_poly : [ list \
...
] ]