QM/MM Modelling with theory=hybrid
There are several ways to set up a QM/MM calculation in ChemShell, but the simplest approach is to use the automatic setup provided by the hybrid module, which sets up a QM/MM model based on one QM and one MM region. The following examples illustrate how it is used.
In these simple examples we will set up the force fields manually using ChemShell's native input format. In later tutorials we will show how more complex forcefields (e.g. CHARMM and AMBER) can be imported.
Mechanical embedding
The first example is pph3_mech.chm in samples/hybrid. This follows closely the PPh3 example in the DL_POLY tutorial. The only change here is that the computational method is chosen as hybrid and the qm_theory and mm_theory arguments are passed to indicate how the QM and MM calculations are performed.
The atoms in the QM region must also be specified using the keyword qm_region. In the example we specify a trivial one-atom QM region.
In this example mechanical embedding is used (coupling=mechanical), so there is no polarisation of the QM region by the MM atoms.
dl-find coords=c maxcycle=1000 theory=hybrid : { coupling=mechanical qm_region= { 1 } \
qm_theory=gamess \
mm_theory=dl_poly : mm_defs=ff }
If you do not have a copy of GAMESS-UK, try substituting the QM code you have been using with ChemShell into qm_theory. Note that most, but not all, of ChemShell's QM interfaces support QM/MM calculations. An error of the following kind indicates that the QM interface you have chosen may not currently support QM/MM:
For qm_theory = ... you need to provide a value (0/1) for the subtract_charge_charge_interactions argument
ChemShell supports QM/MM with at least GAMESS-UK, NWChem, ORCA, Dalton and LSDalton, MNDO (full version), Turbomole, MolPro, Gaussian, and DMol3. This is not a definitive list and support for other codes will be added over time.
Electrostatic embedding
In the second example, h2o_dimer_elec.chm, we consider two water molecules, one modelled by QM and one by an MM force field. In this example we use electrostatic embedding, in which the QM region is polarised by the MM environment. The MM atoms are represented in the QM hamiltonian as point charges taken from the partial charges defined in the MM force field.
Electrostatic embedding is specified by coupling=shift (we will see the significance of the term 'shift' in the third example).
First, we create a geometry as usual, this time a water dimer:
c_create coords=h2o_dimer.c {
coordinates angstrom
o 0.000000 0.000000 0.000000
h 0.000000 -0.751841 0.568201
h -0.000000 0.751841 0.568201
o -2.511474 0.000000 -1.450000
h -1.651960 0.000000 -1.063537
h -2.374291 0.000000 -2.382362
}
Then we need to define a DL_POLY force field which will be used to model one of the water molecules.
read_input ff.dat {
charge o -.796
charge h +.398
powers o o 12 159379 0 0. 0 0.
powers o h 12 44490 0 0. 0 0.
powers h h 12 9569 0 0. 0 0.
bond h o 748.6 0.9424
angle h o h 56.3 105.84
}
Now we are ready to calculate the QM/MM energy with electrostatic embedding. As before we specify theory=hybrid but this time with coupling=shift and a three-atom QM region (the first water molecule). The details of the point charges passed to the QM code to represent the MM environment are handled automatically by ChemShell.
energy energy=e coords=h2o_dimer.c \
theory= hybrid : { coupling=shift
qm_region= {1 2 3}
qm_theory=gamess : basis=3-21g
mm_theory=dl_poly : mm_defs=ff.dat
}
In the example we use GAMESS-UK as the QM driver but again you can try substituting in the QM code of your choice.
Link atoms and the charge shift scheme
In molecular QM/MM models we often face the problem that the QM and MM regions are connected by one or more chemical bonds, which require special handling. There are several good ways to account for these bonds, one of which is the link atom approach used by ChemShell. The link atom method terminates the QM region by the addition of hydrogen atoms along the QM-MM bond, thereby mimicking the effect of the broken bond. ChemShell automatically handles the adjustments to forces and force field terms that are required as a result.
In the region of bonds across the QM/MM boundary, some adjustment of the classical charge distribution near the QM region is also necessary to prevent over-polarisation of the QM region and link atoms. In ChemShell the standard approach is the charge shift scheme, which is automatically applied when coupling=shift is chosen. The charge shift scheme deletes charges on MM atoms bonded to QM atoms and redistributes the charge onto the neighbouring MM atoms. Extra charges are also added to compensate for the dipole created by this shift.
The example butanol_shift.chm is a calculation on butanol using electrostatic embedding with a C-C bond crossing the QM/MM boundary.
First, the geometry is specified as usual:
c_create coords=butanol.c {
coordinates angstrom
C 0.178280 0.355540 0.109010
H -0.009250 0.344290 1.187510
H 1.211650 0.046020 -0.076460
O 0.012900 1.685440 -0.363950
H 0.642900 2.241840 0.123500
C -0.800870 -0.552230 -0.621590
H -1.822920 -0.186120 -0.462760
H -0.623230 -0.483680 -1.702290
C -0.679920 -2.000860 -0.157410
H -0.856380 -2.062020 0.922780
H 0.336810 -2.367040 -0.341370
C -1.672420 -2.900920 -0.875880
H -1.502780 -2.887910 -1.957240
H -2.701650 -2.580310 -0.685340
H -1.568620 -3.933420 -0.527920
}
Next, we define the DL_POLY force field. Note there are four atom types in total as hydrogens bonded to carbon/oxygen are treated differently.
read_input butanol.ff {
declare CT 4
declare HC 1
declare OH 2 3
declare HO 1
subgroup CT CT_qm "qc carbon"
subgroup HC HC_qm "qc hydrogen"
bond CT CT 310.0 1.526
bond CT HC 340.0 1.090
bond CT OH 320.0 1.410
bond HO OH 553.0 0.960
angle CT CT CT 40.0 109.50
angle CT CT HC 50.0 109.50
angle CT CT OH 50.0 109.50
angle HC CT HC 35.0 109.50
angle CT OH HO 55.0 108.50
angle HC CT OH 35.0 109.50
angle HO OH HO 47.0 104.50
ptor * CT CT * 0.1555555 0.0 3 i-j-k-l
ptor * CT OH * 0.1666666 0.0 3 i-j-k-l
vdw HO HO 0.0000E+00 0.0000E+00
vdw HO HC 0.0000E+00 0.0000E+00
vdw HO OH 0.0000E+00 0.0000E+00
vdw HO CT 0.0000E+00 0.0000E+00
vdw HC HC 0.2173E+02 0.7516E+04
vdw HC OH 0.1253E+03 0.6828E+05
vdw HC CT 0.1269E+03 0.9717E+05
vdw OH OH 0.6997E+03 0.5818E+06
vdw OH CT 0.6931E+03 0.7915E+06
vdw CT CT 0.6756E+03 0.1043E+07
}
set forcefield butanol.ff
We then map the atoms to types and assign partial charges.
set types { CT_qm HC_qm HC_qm OH HO CT HC HC CT HC HC CT HC HC HC }
set charges { 0.0 0.15 0.15 -0.7 0.4 -0.18 0.09 0.09 -0.18 0.09 0.09 -0.27 0.09 0.09 0.09 }
These lists need to be passed to the hybrid module when the energy is calculated. The partial charges on the QM atoms will be automatically deleted by ChemShell once the QM region is defined.
Finally we perform a QM/MM energy and gradient calculation. The QM region is defined as the first five atoms. ChemShell uses the connectivity information generated when c_create was run to automatically identify the C-C bond broken, insert a link atom in the appropriate place, and shift the MM charges.
set args [ list coupling=shift \
atom_charges= $charges \
qm_region= { 1 - 5 } \
list_option=full \
mm_theory=dl_poly : [ list mm_defs=$forcefield \
atom_types= $types \
scale14 = { 0.5 0.5 } \
list_option=full \
debug=yes ] \
qm_theory=gamess : [ list hamiltonian=hf basisspec= {{6-31g *}} charge=0 ] ]
eandg coords=butanol.c energy=butanol.e gradient=butanol.g theory=hybrid : [ list $args ]
Note that atom_charges is passed as an argument to theory=hybrid, not as an argument to mm_theory=dl_poly. This is important as the hybrid module modifies the charges according to the charge shift scheme before passing them to the dl_poly module.
The scale14 option specifies how 1-4 electrostatic and VDW interactions are scaled. For illustration here they have both been scaled by 0.5, but other values will be necessary depending on the forcefield used.
The division of the molecule into QM and MM parts must be carefully considered in order to get sensible results. For more on the theory behind electrostatic embedding, link atoms and the charge shift scheme, please see the ChemShell paper (Sherwood et. al., J. Molec. Struc - THEOCHEM, v632 p1, 2003).
Advanced QM/MM modelling
This tutorial is intended to illustrate only the most basic QM/MM features of ChemShell. Setting up large-scale QM/MM calculations inevitably depends on the nature of the system you are trying to model. For example, ChemShell supports import of CHARMM and AMBER force fields for biological modelling, while the GULP interface supports polarisable MM (shell model) force fields for ionic solid state calculations. For further information please see the later QM/MM tutorials.