Polarisable MM forcefields in QM/MM calculations
In the standard electrostatic embedding schemes (shift, L2, rc etc.), the QM region
is polarised by MM atoms via point charges in the QM Hamiltonian, but there is no polarisation of the MM region by the
QM atoms.
Three schemes are available in ChemShell to support polarisable MM forcefields: the shell model,
charge-on-spring/Drude oscillator, and polarised-boundary RC(D).
Shell model
The shell model of Dick and Overhauser allows ionic MM centres to be polarised by splitting them into two charged components: a
point charge core and a spherical shell of charge to model the valence electrons which contribute to the polarisability. The shell
positions can then be adjusted to simulate the effect of polarisation due to the QM atoms.
There is no specific option in theory=hybrid to select shell model polarisation; if there are shells in the fragment then it
will be selected automatically. Mutual self-consistency in electronic polarisation effects between the QM and MM regions
is achieved by a series of microiterations. The QM calculations in the embedding potential are followed by shell relaxations due to the
charge distribution in the QM region at the previous iteration. This is repeated until the shell positions converge. The option max_shell_cycles can be used to control the maximum number of shell relaxation iterations allowed.
Shell model polarisation is supported by the GULP MM interface and is usually used in the simulation of ionic materials. An example of solid state embedding using a shell model potential is provided in the ChemShell tutorial.
Shell model polarisation is tightly integrated into the ChemShell object system and there are a number of commands for
querying and manipulating shell properties - see the fragment commands page for a list of these.
In addition, the cluster generating procedure is fully compatible with fragments containing shells.
Charge-on-spring/Drude oscillator model
The polcos option in theory=hybrid indicates the mutual polarised QM/MM calculations by the GROMOS charge-on-spring(COS)(J. Chem. Theory Comput. 3, 1499 (2007))
and CHARMM Drude oscillator(J. Chem. Phys. 119, 3025 (2003)) models. In the two models,
a massless point charge is attached to the polarised MM heavy atom. Based on the induced dipoles and electric field,
the massless point charge position can be iteratively calculated.
Conceptually these methods are very similar to the
shell model also implemented in ChemShell, but the implementations are separate and the set up procedures are different.
In the polcos case the extra point charges are set up via arguments to the mm_polcos option in theory=hybrid
(fragments containing shells should not be used).
mm_polcos control arguments
| Argument |
Argument type |
Mandatory |
Default |
To specify |
| polcos_efield= |
string |
no |
cos |
Method to calculate the electric field. It has two options--drude(Drude oscillator model) and cos(COS model) |
| polcos_scale14= |
real |
no |
1.0 |
Scale factor of the electrostatic interaction |
| polcos_maxdx= |
real |
no |
0.000018 |
Maximum change(in a.u.) in the position of the massless charge |
| polcos_rmsdx= |
real |
no |
0.000012 |
Root mean square change(in a.u.) in the position of the massless charge |
| polcos_toler_energy= |
real |
no |
1.0e-6 |
Convergence criteria of the QM energy change |
| polcos_maxcycle= |
integer |
no |
10 |
Maximum number of outer iteration |
| polcos_inmaxcycle= |
integer |
no |
1000 |
Maximum number of inner iterations |
| polcos_atom_cosq= |
Tcl List
|
yes |
none |
Polarised atom ID, polarisability(in a.u.) and massless point charge(in a.u.) |
| polcos_rcutl= |
real |
no |
100000.0 |
Cut-off of long range interaction |
| polcos_rcutf= |
real |
no |
0.0 |
Cut-off radius of reaction-field calculation |
| polcos_epsrf= |
real |
no |
1.0 |
Epsilon of reaction-field calculation |
The COS model can only used with the GROMOS force field, and Drude oscillator model with the CHARMM force field.
Moreover, the two models should be applied with the electrostatic embedding schemes. polcos_scale14 must be equal to the value of the MM force field.
polcos_maxcycle indicates the iterative number of QM calculation, and polcos_inmaxcycle is the MM SCF interative number.
Polarisability and the COS point charges parameters of GROMOS force field can be directly used in the cos model after changing the unit of polarisability.
For the Drude oscillator model, the polarisability parameters can be taken from the CHARMM polarised force field after changing the unit.
The point charge parameters are calculated by q=sqrt(a_pol*k_d), where k_d is the force constant of the harmonic spring
from CHARMM polarised force field, and a_pol is polarisability parameters.
polcos_rcutl, polcos_rcutf and polcos_epsrf should be equal to the value of the corresponding keyword
controlling the GROMOS MM calculation. The three keywords are only used in the COS model calling GROMOS96 package.
The defaults of the three keywords mean that cut-off of long-range interaction and reaction-field model are not applied.
Thus, the rcrf and rcutl keywords to control the the GROMOS MM calculation should be set to large enough.
The epsrf keyword should be set to 1.0. In the Drude oscillator QM/MM calculation using the CHARMM parameters,
there is no cut-off of the long-range electrostatic interaction. For the consistency, the cut-off parameters in the CHARMM parameter file must be set to large enough.
However, one should be careful to set the cut-off parameters in the polarised QM/MM calculations using COS and Drude oscillator models.
mm_polcos example
set args [ list atom_types= $types \
atom_charges= $charges \
scale14= {1.0 0.0} \
exact_srf= yes \
use_charmm_psf=yes \
charmm_psf_file=wat1000.psf \
charmm_parameter_file=spc_water_param.inp \
charmm_mass_file=spc_water_topology.inp \
]
set polarg [ list polcos_scale14=1.0 \
polcos_maxcycle=20 \
polcos_toler_energy=1.0e-8 \
polcos_maxdx=2.0e-5 \
polcos_rmsdx=1.0e-5 \
polcos_efield=drude \
polcos_inmaxcyc=1000 \
polcos_atom_polcosq= $cospara \
]
set hybridargs [ list coupling=shift \
debug=off \
qm_region= $qm_atoms \
list_option=full \
mm_polcos=yes : [ list $polarg ] \
groups= $cgroups \
qm_theory= mndo : [ list hamiltonian= om3 \
charge=0 \
maxcyc=200 \
] \
mm_theory=dl_poly: $args \
]
eandg coords= wat1000.c energy= wat1000.e gradient=wat1000.g theory= hybrid : [ list $hybridargs ]
Polarised-boundary RC(D) model
The mm_polboundary option in theory=hybrid indicates the polarised-boundary RC and RCD treatments(J. Chem. Theory Comput. 3, 1378 (2007)).
In the treatment, the MM point charges on the QM/MM boundary-region atoms are adjusted based on the principles of the electronegativity equalization
until the charge distributions on the polarised MM atoms converge. It should only be used with the RC or RCD schemes and Gaussian03 package.
mm_polboundary control arguments
| Argument |
Argument type |
Mandatory |
Default |
To specify |
| polboundary_charge_method= |
string |
no |
qeqrg |
Method to calculate the charges on the polarised boundary-region MM atoms. It has three options: qeqrg(charge equalization method proposed by Rappe and Goddard), qeqbt(a modified version of the charge equalization method by Bakowies and Thiel), and eem(electronegativity equalization method of Mortier and coworkers). |
| polboundary_charge_region= |
Tcl List
|
no |
M2 and M3 |
MM atoms ID in the polarised-boundary region |
| polboundary_maxdq= |
real |
no |
0.005 |
Maximum change(in a.u.) of the partial atomic charges on the polarised-boundary MM atoms |
| polboundary_rmsdq= |
real |
no |
0.002 |
Root mean square change(in a.u.) of the partial atomic charges on the polarised-boundary MM atoms |
| polboundary_maxcycle= |
integer |
no |
10 |
Maximum number of iterations |
|
