docs/reference-manual/special/qmmm.rst

   1 Mixed Quantum-Classical simulation techniques
   2 ---------------------------------------------
   3
   4 In a molecular mechanics (MM) force field, the influence of electrons is
   5 expressed by empirical parameters that are assigned on the basis of
   6 experimental data, or on the basis of results from high-level quantum
   7 chemistry calculations. These are valid for the ground state of a given
   8 covalent structure, and the MM approximation is usually sufficiently
   9 accurate for ground-state processes in which the overall connectivity
  10 between the atoms in the system remains unchanged. However, for
  11 processes in which the connectivity does change, such as chemical
  12 reactions, or processes that involve multiple electronic states, such as
  13 photochemical conversions, electrons can no longer be ignored, and a
  14 quantum mechanical description is required for at least those parts of
  15 the system in which the reaction takes place.
  16
  17 One approach to the simulation of chemical reactions in solution, or in
  18 enzymes, is to use a combination of quantum mechanics (QM) and molecular
  19 mechanics (MM). The reacting parts of the system are treated quantum
  20 mechanically, with the remainder being modeled using the force field.
  21 The current version of |Gromacs| provides interfaces to several popular
  22 Quantum Chemistry packages (MOPAC :ref:`150 <refmopac>`,
  23 GAMESS-UK \ :ref:`151 <refgamess-uk>`, Gaussian \ :ref:`152 <refg03>` and
  24 CPMD \ :ref:`153 <refCar85a>`).
  25
  26 |Gromacs| interactions between the two subsystems are either handled as
  27 described by Field et al. :ref:`154 <refField90a>` or within
  28 the ONIOM approach by Morokuma and coworkers \ :ref:`155 <refMaseras96a>`,
  29 :ref:`156 <refSvensson96a>`.
  30
  31 Overview
  32 ^^^^^^^^
  33
  34 Two approaches for describing the interactions between the QM and MM
  35 subsystems are supported in this version:
  36
  37 #. **Electronic Embedding** The electrostatic interactions between the
  38    electrons of the QM region and the MM atoms and between the QM nuclei
  39    and the MM atoms are included in the Hamiltonian for the QM
  40    subsystem:
  41
  42    .. math::
  43
  44       H^{QM/MM} =
  45       H^{QM}_e-\sum_i^n\sum_J^M\frac{e^2Q_J}{4\pi\epsilon_0r_{iJ}}+\sum_A^N\sum_J^M\frac{e^2Z_AQ_J}{e\pi\epsilon_0R_{AJ}},
  46
  47    where :math:`n` and :math:`N` are the number of electrons and nuclei
  48    in the QM region, respectively, and :math:`M` is the number of
  49    charged MM atoms. The first term on the right hand side is the
  50    original electronic Hamiltonian of an isolated QM system. The first
  51    of the double sums is the total electrostatic interaction between the
  52    QM electrons and the MM atoms. The total electrostatic interaction of
  53    the QM nuclei with the MM atoms is given by the second double sum.
  54    Bonded interactions between QM and MM atoms are described at the MM
  55    level by the appropriate force-field terms. Chemical bonds that
  56    connect the two subsystems are capped by a hydrogen atom to complete
  57    the valence of the QM region. The force on this atom, which is
  58    present in the QM region only, is distributed over the two atoms of
  59    the bond. The cap atom is usually referred to as a link atom.
  60
  61 #. **ONIOM** In the ONIOM approach, the energy and gradients are first
  62    evaluated for the isolated QM subsystem at the desired level of *ab
  63    initio* theory. Subsequently, the energy and gradients of the total
  64    system, including the QM region, are computed using the molecular
  65    mechanics force field and added to the energy and gradients
  66    calculated for the isolated QM subsystem. Finally, in order to
  67    correct for counting the interactions inside the QM region twice, a
  68    molecular mechanics calculation is performed on the isolated QM
  69    subsystem and the energy and gradients are subtracted. This leads to
  70    the following expression for the total QM/MM energy (and gradients
  71    likewise):
  72
  73    .. math::
  74
  75       E_{tot} = E_{I}^{QM}
  76       +E_{I+II}^{MM}-E_{I}^{MM},
  77
  78    where the subscripts I and II refer to the QM and MM subsystems,
  79    respectively. The superscripts indicate at what level of theory the
  80    energies are computed. The ONIOM scheme has the advantage that it is
  81    not restricted to a two-layer QM/MM description, but can easily
  82    handle more than two layers, with each layer described at a different
  83    level of theory.
  84
  85 Usage
  86 ^^^^^
  87
  88 To make use of the QM/MM functionality in |Gromacs|, one needs to:
  89
  90 #. introduce link atoms at the QM/MM boundary, if needed;
  91
  92 #. specify which atoms are to be treated at a QM level;
  93
  94 #. specify the QM level, basis set, type of QM/MM interface and so on.
  95
  96 Adding link atoms
  97 ^^^^^^^^^^^^^^^^^
  98
  99 At the bond that connects the QM and MM subsystems, a link atoms is
 100 introduced. In |Gromacs| the link atom has special atomtype, called LA.
 101 This atomtype is treated as a hydrogen atom in the QM calculation, and
 102 as a virtual site in the force-field calculation. The link atoms, if
 103 any, are part of the system, but have no interaction with any other
 104 atom, except that the QM force working on it is distributed over the two
 105 atoms of the bond. In the topology, the link atom (LA), therefore, is
 106 defined as a virtual site atom:
 107
 108 ::
 109
 110     [ virtual_sites2 ]
 111     LA QMatom MMatom 1 0.65
 112
 113 See sec. :ref:`vsitetop` for more details on how virtual sites are
 114 treated. The link atom is replaced at every step of the simulation.
 115
 116 In addition, the bond itself is replaced by a constraint:
 117
 118 ::
 119
 120     [ constraints ]
 121     QMatom MMatom 2 0.153
 122
 123 **Note** that, because in our system the QM/MM bond is a carbon-carbon
 124 bond (0.153 nm), we use a constraint length of 0.153 nm, and dummy
 125 position of 0.65. The latter is the ratio between the ideal C-H bond
 126 length and the ideal C-C bond length. With this ratio, the link atom is
 127 always 0.1 nm away from the ``QMatom``, consistent with the carbon-hydrogen
 128 bond length. If the QM and MM subsystems are connected by a different
 129 kind of bond, a different constraint and a different dummy position,
 130 appropriate for that bond type, are required.
 131
 132 Specifying the QM atoms
 133 ^^^^^^^^^^^^^^^^^^^^^^^
 134
 135 Atoms that should be treated at a QM level of theory, including the link
 136 atoms, are added to the index file. In addition, the chemical bonds
 137 between the atoms in the QM region are to be defined as connect bonds
 138 (bond type 5) in the topology file:
 139
 140 ::
 141
 142     [ bonds ]
 143     QMatom1 QMatom2 5
 144     QMatom2 QMatom3 5
 145
 146 Specifying the QM/MM simulation parameters
 147 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 148
 149 In the :ref:`mdp` file, the following parameters control a
 150 QM/MM simulation.
 151
 152 ``QMMM = no``
 153     | If this is set to ``yes``, a QM/MM simulation is
 154       requested. Several groups of atoms can be described at different
 155       QM levels separately. These are specified in the QMMM-grps field
 156       separated by spaces. The level of *ab initio* theory at which the
 157       groups are described is specified by ``QMmethod`` and
 158       ``QMbasis`` Fields. Describing the groups at different
 159       levels of theory is only possible with the ONIOM QM/MM scheme,
 160       specified by ``QMMMscheme``.
 161
 162 ``QMMM-grps =``
 163     | groups to be described at the QM level
 164
 165 ``QMMMscheme = normal``
 166     | Options are ``normal`` and ``ONIOM``. This
 167       selects the QM/MM interface. ``normal`` implies that
 168       the QM subsystem is electronically embedded in the MM subsystem.
 169       There can only be one ``QMMM-grps`` that is modeled at
 170       the ``QMmethod`` and ``QMbasis`` level of
 171       * ab initio* theory. The rest of the system is described at the MM
 172       level. The QM and MM subsystems interact as follows: MM point
 173       charges are included in the QM one-electron Hamiltonian and all
 174       Lennard-Jones interactions are described at the MM level. If
 175       ``ONIOM`` is selected, the interaction between the
 176       subsystem is described using the ONIOM method by Morokuma and
 177       co-workers. There can be more than one QMMM-grps each modeled at a
 178       different level of QM theory (QMmethod and QMbasis).
 179
 180 ``QMmethod =``
 181     | Method used to compute the energy and gradients on the QM atoms.
 182       Available methods are AM1, PM3, RHF, UHF, DFT, B3LYP, MP2, CASSCF,
 183       MMVB and CPMD. For CASSCF, the number of electrons and orbitals
 184       included in the active space is specified by
 185       ``CASelectrons`` and ``CASorbitals``. For
 186       CPMD, the plane-wave cut-off is specified by the
 187       ``planewavecutoff`` keyword.
 188
 189 ``QMbasis =``
 190     | Gaussian basis set used to expand the electronic wave-function.
 191       Only Gaussian basis sets are currently available, i.e. STO-3G,
 192       3-21G, 3-21G\*, 3-21+G\*, 6-21G, 6-31G, 6-31G\*, 6-31+G\*, and
 193       6-311G. For CPMD, which uses plane wave expansion rather than
 194       atom-centered basis functions, the ``planewavecutoff``
 195       keyword controls the plane wave expansion.
 196
 197 ``QMcharge =``
 198     | The total charge in *e* of the ``QMMM-grps``. In case
 199       there are more than one ``QMMM-grps``, the total
 200       charge of each ONIOM layer needs to be specified separately.
 201
 202 ``QMmult =``
 203     | The multiplicity of the ``QMMM-grps``. In case there
 204       are more than one ``QMMM-grps``, the multiplicity of
 205       each ONIOM layer needs to be specified separately.
 206
 207 ``CASorbitals =``
 208     | The number of orbitals to be included in the active space when
 209       doing a CASSCF computation.
 210
 211 ``CASelectrons =``
 212     | The number of electrons to be included in the active space when
 213       doing a CASSCF computation.
 214
 215 ``SH = no``
 216     | If this is set to yes, a QM/MM MD simulation on the excited
 217       state-potential energy surface and enforce a diabatic hop to the
 218       ground-state when the system hits the conical intersection
 219       hyperline in the course the simulation. This option only works in
 220       combination with the CASSCF method.
 221
 222 Output
 223 ^^^^^^
 224
 225 The energies and gradients computed in the QM calculation are added to
 226 those computed by |Gromacs|. In the :ref:`edr` file there is a
 227 section for the total QM energy.
 228
 229 Future developments
 230 ^^^^^^^^^^^^^^^^^^^
 231
 232 Several features are currently under development to increase the
 233 accuracy of the QM/MM interface. One useful feature is the use of
 234 delocalized MM charges in the QM computations. The most important
 235 benefit of using such smeared-out charges is that the Coulombic
 236 potential has a finite value at interatomic distances. In the point
 237 charge representation, the partially-charged MM atoms close to the QM
 238 region tend to “over-polarize” the QM system, which leads to artifacts
 239 in the calculation.
 240
 241 What is needed as well is a transition state optimizer.