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 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

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