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 following documentation assumes option 1.

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

The basic approach used by the hybrid module can be summarised as follows.
  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

The following arguments control the data flow and main program options for QM/MM calculations.

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

The connectivity of the molecule is used in the set up of the QM/MM calculation, and should be carefully checked. In particular it will control:
  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.

Argument Argument type Mandatory Default To specify
dipole_adjust= Boolean no no Make a bond dipole correction to the charges of the first layer of MM atoms. This allow you to replace a non-neutral group of atoms with a neutral H-terminated cluster. As programmed it will only work for aluminosilicates with q(Si)=1.2 and Q(O)=-0.6. See the zeolite tutorial for an example.
max_shell_cycles= integer no 5 (if shells present) Number of QM/shell relaxation cycles. See Polarisable MM.
mm_polcos= Boolean no no Charge-on-spring/drude oscillator model for mutual polarisation between QM and MM regions. See Polarisable MM.
mm_polboundary= Boolean no no Polarised-boundary RC(D) model for mutual polarisation between QM and MM regions. See Polarisable MM.
smbp= Boolean no no Whether to use the solvated macromolecule boundary potential (SMBP) to represent the outer macromolecule and solvent region. If the SMBP is applied further input has to be provided. See Boundary potentials.
gsbp= Boolean no no Wether to use the generalized solvent boundary potential (GSBP) to represent the outer macromolecule and solvent region. If the SMBP is applied QM/MM input has to be modified useing the gsbp.prepare procedure and further input has to be provided. See Boundary potentials.
npert= integer no 0 Number of perturbations in Free-Energy Perturbation calculations (FEP). 1 for forward perturbation and 2 for forward and backward perturbation. 0 for standard calculation (no FEP).
fep_linkmethod= integer no 4 Method for treating link atoms in FEP. Range: 1..4. See FEP calculations.
fep_qm_theory= string no qm_theory QM module used the perturbation in FEP calculations. Same syntax as qm_theory.
fep_denmat= string no see below Name of the file in which the density or ESP charge information is stored (has to be read by the QM code). See FEP calculations.
blur_width= Tcl List no none (Developmental option)
double_link_atom= Tcl List no none (Developmental option)




This manual was generated using htp, an HTML pre-processor Powered by htp