Mixed Quantum-Classical simulation techniques

In a molecular mechanics (MM) force field, the influence of electrons is expressed by empirical parameters that are assigned on the basis of experimental data, or on the basis of results from high-level quantum chemistry calculations. These are valid for the ground state of a given covalent structure, and the MM approximation is usually sufficiently accurate for ground-state processes in which the overall connectivity between the atoms in the system remains unchanged. However, for processes in which the connectivity does change, such as chemical reactions, or processes that involve multiple electronic states, such as photochemical conversions, electrons can no longer be ignored, and a quantum mechanical description is required for at least those parts of the system in which the reaction takes place.

One approach to the simulation of chemical reactions in solution, or in enzymes, is to use a combination of quantum mechanics (QM) and molecular mechanics (MM). The reacting parts of the system are treated quantum mechanically, with the remainder being modeled using the force field. The current version of GROMACS provides interfaces to several popular Quantum Chemistry packages (MOPAC 150, GAMESS-UK 151, Gaussian 152 and CPMD 153).

GROMACS interactions between the two subsystems are either handled as described by Field et al. 154 or within the ONIOM approach by Morokuma and coworkers 155, 156.

Overview

Two approaches for describing the interactions between the QM and MM subsystems are supported in this version:

  1. Electronic Embedding The electrostatic interactions between the electrons of the QM region and the MM atoms and between the QM nuclei and the MM atoms are included in the Hamiltonian for the QM subsystem:

    \[H^{QM/MM} = 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}},\]

    where \(n\) and \(N\) are the number of electrons and nuclei in the QM region, respectively, and \(M\) is the number of charged MM atoms. The first term on the right hand side is the original electronic Hamiltonian of an isolated QM system. The first of the double sums is the total electrostatic interaction between the QM electrons and the MM atoms. The total electrostatic interaction of the QM nuclei with the MM atoms is given by the second double sum. Bonded interactions between QM and MM atoms are described at the MM level by the appropriate force-field terms. Chemical bonds that connect the two subsystems are capped by a hydrogen atom to complete the valence of the QM region. The force on this atom, which is present in the QM region only, is distributed over the two atoms of the bond. The cap atom is usually referred to as a link atom.

  2. ONIOM In the ONIOM approach, the energy and gradients are first evaluated for the isolated QM subsystem at the desired level of ab initio theory. Subsequently, the energy and gradients of the total system, including the QM region, are computed using the molecular mechanics force field and added to the energy and gradients calculated for the isolated QM subsystem. Finally, in order to correct for counting the interactions inside the QM region twice, a molecular mechanics calculation is performed on the isolated QM subsystem and the energy and gradients are subtracted. This leads to the following expression for the total QM/MM energy (and gradients likewise):

    \[E_{tot} = E_{I}^{QM} +E_{I+II}^{MM}-E_{I}^{MM},\]

    where the subscripts I and II refer to the QM and MM subsystems, respectively. The superscripts indicate at what level of theory the energies are computed. The ONIOM scheme has the advantage that it is not restricted to a two-layer QM/MM description, but can easily handle more than two layers, with each layer described at a different level of theory.

Usage

To make use of the QM/MM functionality in GROMACS, one needs to:

  1. introduce link atoms at the QM/MM boundary, if needed;
  2. specify which atoms are to be treated at a QM level;
  3. specify the QM level, basis set, type of QM/MM interface and so on.

Specifying the QM atoms

Atoms that should be treated at a QM level of theory, including the link atoms, are added to the index file. In addition, the chemical bonds between the atoms in the QM region are to be defined as connect bonds (bond type 5) in the topology file:

[ bonds ]
QMatom1 QMatom2 5
QMatom2 QMatom3 5

Specifying the QM/MM simulation parameters

In the mdp file, the following parameters control a QM/MM simulation.

QMMM = no
If this is set to yes, a QM/MM simulation is requested. Several groups of atoms can be described at different QM levels separately. These are specified in the QMMM-grps field separated by spaces. The level of ab initio theory at which the groups are described is specified by QMmethod and QMbasis Fields. Describing the groups at different levels of theory is only possible with the ONIOM QM/MM scheme, specified by QMMMscheme.
QMMM-grps =
groups to be described at the QM level
QMMMscheme = normal
Options are normal and ONIOM. This selects the QM/MM interface. normal implies that the QM subsystem is electronically embedded in the MM subsystem. There can only be one QMMM-grps that is modeled at the QMmethod and QMbasis level of * ab initio* theory. The rest of the system is described at the MM level. The QM and MM subsystems interact as follows: MM point charges are included in the QM one-electron Hamiltonian and all Lennard-Jones interactions are described at the MM level. If ONIOM is selected, the interaction between the subsystem is described using the ONIOM method by Morokuma and co-workers. There can be more than one QMMM-grps each modeled at a different level of QM theory (QMmethod and QMbasis).
QMmethod =
Method used to compute the energy and gradients on the QM atoms. Available methods are AM1, PM3, RHF, UHF, DFT, B3LYP, MP2, CASSCF, MMVB and CPMD. For CASSCF, the number of electrons and orbitals included in the active space is specified by CASelectrons and CASorbitals. For CPMD, the plane-wave cut-off is specified by the planewavecutoff keyword.
QMbasis =
Gaussian basis set used to expand the electronic wave-function. Only Gaussian basis sets are currently available, i.e. STO-3G, 3-21G, 3-21G*, 3-21+G*, 6-21G, 6-31G, 6-31G*, 6-31+G*, and 6-311G. For CPMD, which uses plane wave expansion rather than atom-centered basis functions, the planewavecutoff keyword controls the plane wave expansion.
QMcharge =
The total charge in e of the QMMM-grps. In case there are more than one QMMM-grps, the total charge of each ONIOM layer needs to be specified separately.
QMmult =
The multiplicity of the QMMM-grps. In case there are more than one QMMM-grps, the multiplicity of each ONIOM layer needs to be specified separately.
CASorbitals =
The number of orbitals to be included in the active space when doing a CASSCF computation.
CASelectrons =
The number of electrons to be included in the active space when doing a CASSCF computation.
SH = no
If this is set to yes, a QM/MM MD simulation on the excited state-potential energy surface and enforce a diabatic hop to the ground-state when the system hits the conical intersection hyperline in the course the simulation. This option only works in combination with the CASSCF method.

Output

The energies and gradients computed in the QM calculation are added to those computed by GROMACS. In the edr file there is a section for the total QM energy.

Future developments

Several features are currently under development to increase the accuracy of the QM/MM interface. One useful feature is the use of delocalized MM charges in the QM computations. The most important benefit of using such smeared-out charges is that the Coulombic potential has a finite value at interatomic distances. In the point charge representation, the partially-charged MM atoms close to the QM region tend to “over-polarize” the QM system, which leads to artifacts in the calculation.

What is needed as well is a transition state optimizer.