1 Mixed Quantum-Classical simulation techniques
2 ---------------------------------------------
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.
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>`).
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>`.
34 Two approaches for describing the interactions between the QM and MM
35 subsystems are supported in this version:
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
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}},
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.
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
76 +E_{I+II}^{MM}-E_{I}^{MM},
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
88 To make use of the QM/MM functionality in |Gromacs|, one needs to:
90 #. introduce link atoms at the QM/MM boundary, if needed;
92 #. specify which atoms are to be treated at a QM level;
94 #. specify the QM level, basis set, type of QM/MM interface and so on.
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:
111 LA QMatom MMatom 1 0.65
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.
116 In addition, the bond itself is replaced by a constraint:
121 QMatom MMatom 2 0.153
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.
132 Specifying the QM atoms
133 ^^^^^^^^^^^^^^^^^^^^^^^
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:
146 Specifying the QM/MM simulation parameters
147 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
149 In the :ref:`mdp` file, the following parameters control a
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``.
163 | groups to be described at the QM level
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).
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.
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.
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.
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.
208 | The number of orbitals to be included in the active space when
209 doing a CASSCF computation.
212 | The number of electrons to be included in the active space when
213 doing a CASSCF computation.
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.
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.
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
241 What is needed as well is a transition state optimizer.