Hybrid Ab-Initio / Force-Field Calculations

Introduction

The hybrid QM/MM implementation consists of a family of functions, hybrid.init hybrid.energy etc, following the general interface of energy/gradient evaluators of ChemShell. Usually the user will not invoke these functions directly, instead using a driver module such as dl-find to control a geometry optimisation, or dynamics to run a MD simulation as described in "Tasks" in this manual.

The details of the QM/MM calculation are controlled by the elements of a global Tcl array (named hybrid). There are three ways to establish this control information.

1. The argument theory=hybrid to any driver module will result in automatic setup based on one QM and one MM region. This option is recommended.
2. The hybrid_setup module may be used to set the array elements from input provided by a z-matrix style data stream, with extra fields added to the atom records to control the QM/MM calculation.
3. The array elements may be set up manually.

The setup of QM/MM calculations is generally more evolved than basic quantum chemical simulations. Examples of how to set up QM/MM simulations of an enzyme, a cluster of an inorganic bulk system, and a cluser representing a surface are avaiable in the tutorial.

Methodology

1. The total energy is taken to be the sum of the MM and QM energies. There is no subtraction as in (for example the IMOMO) scheme, and the MM energy is designed to exclude the internal energy of the MM part.
2. Link atoms are placed on the bond axis by the hybrid module, and are invisible to any controlling "driver" module (such as MD or optimisation). The link atom coordinates are thus not part of the molecule definition. The fact that the link atom positions are determined by the QM and MM atom positions allows the analytic derivative of the QM/MM energy to be corrected for forces acting on the link atoms.
3. The MM energy expression excludes terms at the junction where these are implicitly computed by the QM code as a result of the link atom constraints. e.g., changes in the angle formed by atoms M....Q....Q (M=MM atom, Q=QM atom) cause the movement of the link atom relative to the Q---Q direction, changing the QM energy. Thus an MM term is not needed for this.
4. The electrostatic coupling may be handled purely by the MM code (mechanical coupling) or by inclusion of the MM centres as point charges in the QM code. In the latter case a charge shift correction may be applied, in which charge present on the first layer of MM atoms (ie those connected to the QM region) is moved to the second layer. A dipole (constructed from two point charges) is placed on the recipient atom of this charge shift, so that the charge and dipole of the MM system are conserved after the shift. An alternative to the charge shift is the deletion of the charges on the first neutral charge group.
5. If shells are present on the QM fragment and an electrostatic coupling model is chosen, the QM electrostatic forces on the shells will be recovered, and a number of iterations of the sequence MM Energy, QM Energy will be performed to generate consistent shell positions and charge density. See the Polarisable MM page for further details.
6. QM/MM calculations using periodic boundary conditions are possible with some limitations. Mechanical embedding works. Electrostatic embedding is only available using a QM/MM cutoff. When using Ewald summation for the MM part, in this case, the MM electrostatics will not treat the QM part nor its periodic images.
7. The double link atom/Gaussian blur scheme suggested by Brooks is available but is in a developmental state and is not documented here. Please contact P. Sherwood if you wish to use it.
8. A QM/MM coupling based on the Direct Reaction Field is also under development but is incomplete, please ignore references to DRF.

Standard control arguments for theory=hybrid

Argument	Argument type	Mandatory	Default	To specify
mm_theory=	string	no	mm	module used to compute MM energy and forces
qm_theory=	string	no	mopac	module used to compute QM energy and forces
qm_region=	Tcl List	no	none	List of atoms in the QM part
coupling=	keyword	no	mechanical	how to compute QM/MM electrostatic interaction, see below.
coords=	Fragment tag	yes	none	initial configuration of the system
energy=	Matrix tag	yes	none	where to store the QM/MM energy (see note 1)
gradient=	Matrix tag	yes	none	where to store the QM/MM gradient
list_option=	Output keyword	no	medium	how much output to generate
groups=	Tcl List	no	none	Neutral group specification
cutoff=	real	no	none	Cutoff for QM/MM electrostatic interactions
conn=	Fragment tag	no	use coords=	Read in connectivity information from the specified object at initialisation.
calculate_connectivity=	Boolean	no	no	(Re-)calculate connectivity of the molecule at initialisation. Ignored if conn= is specified. Keywords scale and toler can be specified to pass to the connect command.

Note 1. If the hybrid module is running under control of another module (e.g. an optimiser) the user will not need to provide values for the energy, gradient and coords arguments as these arguments will be provided by the calling module.

mm_theory and qm_theory

The mm_theory argument must be set to either dl_poly, gulp, gromos, or charmm at present. For biological calculations dl_poly can be used with CHARMM and AMBER topology and force field parameters (see the dl_poly page for more details).

The qm_theory argument may be set to gamess, turbomole, mndo, molpro, orca, dmol3, fhi, qchem, gaussian, nwchem, dalton, mopac, or pointcharge, depending on the availability of QC codes. Interfaces to other QC codes may be added in a straightforward manner.

The function of these arguments is to allow, via the : syntax, arguments to be passed to the MM and the QM scripts.

Coupling models

Keyword	Coupling Scheme
mechanical	QM...MM electrostatic interactions handled by the MM code. No QM polarisation
shift (or polarised)	QM...MM electrostatic interactions are handled by the QM code with QM polarisation. Charge close to the link atoms is shifted away, and point dipoles added to compensate
L2	L2 scheme (Thiel classification) in which the charges of the first charge group are deleted from the QM/MM interaction.
rc, rcd	Redistributed charge (RC)/redistributed charge and dipole (RCD) schemes. (J. Phys. Chem. A 109, 3991 (2005)). The RCD scheme is similar to the shift method.

For electrostatic embedding calculations the usual choice is coupling=shift.

The Importance of Connectivity

1. Automatic addition of link atoms, using the qm_region argument. This will occur whenever there are connections between the QM and the MM region
2. The shifted electrostatic correction is applied according to bonds present between the first and second layers of MM atoms.

The calculation of connectivity information can be controlled using the conn and calculate_connectivity keywords.

Advanced theory=hybrid options

These methods are covered in more detail in the indicated pages.