QM/MM

The hybrid QM/MM (quantum mechanics/molecular mechanics) approach is a molecular simulation method that combines the strengths of ab initio QM calculations (accuracy) and MM (speed) approaches, thus allowing for the study of chemical processes in solution and in proteins. The QM/MM approach was introduced in the 1976 paper of Warshel and Levitt.^[1] They, along with Martin Karplus, won the 2013 Nobel Prize in Chemistry for "the development of multiscale models for complex chemical systems".^[2][3]

An important advantage of QM/MM methods is their efficiency. The cost of doing classical molecular mechanics (MM) simulations in the most straightforward case scales as O(N²), where N is the number of atoms in the system. This is mainly due to electrostatic interactions term (every particle interacts with everything else). However, use of cutoff radius, periodic pair-list updates and more recently the variations of the particle mesh Ewald (PME) method has reduced this to between O(N) to O(N²). In other words, if a system with twice as many atoms is simulated then it would take between twice to four times as much computing power. On the other hand, the simplest ab initio calculations formally scale as O(N³) or worse (restricted Hartree–Fock calculations have been suggested to scale ~O(N^2.7)). Here in the ab initio calculations, N stands for the number of basis functions rather than the number of atoms. Each atom has at least as many basis functions as is the number of electrons. To overcome the limitation, a small part of the system that is of major interest is treated quantum-mechanically (for instance, the active site of an enzyme) and the remaining system is treated classically.^[4][5]

The energy of the combined system may be calculated in two different ways. The simplest is referred to as the 'subtractive scheme' which was proposed by Maseras and Morokuma in 1995. In the subtractive scheme the energy of the entire system is calculated using a molecular mechanics force field, then the energy of the QM system is added (calculated using a QM method), finally the MM energy of the QM system is subtracted.

${\displaystyle E=E^{QM}(QM)+E^{MM}(QM+MM)-E^{MM}(QM)}$

In this equation ${\displaystyle E^{MM}(QM)}$ would refer to the energy of the QM region as calculated using molecular mechanics. In this scheme, the interaction between the two regions will only be considered at a MM level of theory.

In practice, a more widely used approach is a more accurate, additive method. The equation for this consists of three terms: ${\displaystyle {\begin{aligned}E(QM/MM)&=\sum _{I'=1}^{MM}\left[\int \operatorname {d} \mathbf {r} {q_{I'}\rho (\mathbf {r} ) \over \left\vert \mathbf {R} _{I'}-\mathbf {r} \right\vert }+\sum _{I=1}^{QM}{\frac {q_{I'}q_{I}}{\left\vert \mathbf {R} _{I'}-\mathbf {R} _{I}\right\vert }}\right]\\&+\sum _{\text{non-bonded pairs}}\left({\frac {A_{II'}}{R_{II'}^{12}}}-{\frac {B_{II'}}{R_{II'}^{6}}}\right)+\sum _{\text{bonds}}k_{r}(R_{II'}-r_{0})^{2}\\&+\sum _{\text{angles}}k_{\theta }(\theta -\theta _{0})^{2}+\sum _{\text{torsions}}\sum _{n}k_{\phi ,n}[\cos(n\phi +\delta _{n})+1]\end{aligned}}}$

The index ${\displaystyle I}$ labels the nuclei in the QM region whereas ${\displaystyle I'}$ labels the MM nuclei. The first two terms represent the interaction between the total charge density (due to electrons and cores) in the QM region and classical charges of the MM region. The third term accounts for dispersion interactions across the QM/MM boundary. Any covalent bond-stretching potentials that cross the boundary are accounted for by the fourth term. The final two terms account for the energy across the boundary that arises from bending covalent bonds and torsional potentials. At least one of the atoms in the angles ${\displaystyle \theta }$ or ${\displaystyle \phi }$ will be a QM atom with the others being MM atoms.^{[6]: 422–3}

Evaluating the charge-charge term in the QM/MM interaction equation given previously can be very computationally expensive (consider the number of evaluations required a system with 10⁶ grid points for the electron density of the QM system and 10⁴ MM atoms). A method by which this issue can be mitigated is to construct three concentric spheres around the QM region and evaluate which one of these spheres the MM atoms lie within. If the MM atoms reside within the innermost sphere their interactions with the QM system are treated as per the equation for ${\displaystyle E(QM/MM)}$. The MM charges that lie within the second sphere (but not the first) interact with the QM region by giving the QM nuclei constructed charges. These charges are determined by the RESP approach in an attempt to mimic electron density. Using this approach the changing charges on the QM nuclei during the course of a simulation are accounted for.

In the third outermost region the classical charges interact with the multipole moments of the quantum charge distribution. By calculating charge-charge interactions by using successively more approximate methods it is possible to obtain a very significant reduction in computational cost whilst not suffering a significant loss in accuracy.^{[6]: 423–4}

Electrostatic interactions between the QM and MM region may be considered at different levels of sophistication. These methods can be classified as either mechanical embedding, electrostatic embedding or polarized embedding.

Mechanical embedding treats the electrostatic interactions at the MM level, though simpler than the other methods, certain issues may occur, in part due to the extra difficulty in assigning appropriate MM properties such as atom centered point charges to the QM region. The QM region being simulated is the site of the reaction thus it is likely that during the course of the reaction the charge distribution will change resulting in a high level of error if a single set of MM electrostatic parameters is used to describe it. Another problem is the fact that mechanical embedding will not consider the effects of electrostatic interactions with the MM system on the electronic structure of the QM system.^[7]

Electrostatic embedding does not require the MM electrostatic parameters for the QM. This is due to it considering the effects of the electrostatic interactions by including certain one electron terms in the QM regions Hamiltonian. This means that polarization of the QM system by the electrostatic interactions with the MM system will now be accounted for. Though an improvement on the mechanical embedding scheme it comes at the cost of increased complexity hence requiring more computational effort. Another issue is it neglects the effects of the QM system on the MM system whereas in reality both systems would polarize each other until an equilibrium is met.

In order to construct the required one electron terms for the MM region it is possible to utilize the partial charges described by the MM calculation. This is the most popular method for constructing the QM Hamiltonian however it may not be suitable for all systems.^[7]

Whereas electrostatic embedding accounts for the polarisation of the QM system by the MM system, neglecting the polarization of the MM system by the QM system, polarized embedding accounts for both the polarization of the MM system by the QM. These models allow for flexible MM charges and fall into two categories. In the first category, the MM region is polarized by the QM electric field but then does not act back on the QM system. In the second category are fully self-consistent formulations which allow for mutual polarization between the QM and the MM systems. Polarized embedding schemes have scarcely been applied to bio-molecular simulations and have essentially been restricted to explicit solvation models where the solute will be treated as a QM system and the solvent a polarizable force field.^[7]

Even though QM/MM methods are often very efficient, they are still rather tricky to handle. A researcher has to limit the regions (atomic sites) which are simulated by QM, however methods have been developed that allow particles to move between the QM and MM region.^[8] Moving the limitation borders can both affect the results and the time computing the results. The way the QM and MM systems are coupled can differ substantially depending on the arrangement of particles in the system and their deviations from equilibrium positions in time. Usually, limits are set at carbon-carbon bonds and avoided in regions that are associated with charged groups, since such an electronically variant limit can influence the quality of the model.^[9]

Directly connected atoms, where one is described by QM and the other by MM are referred to as Junction atoms. Having the boundary between the QM region and MM region pass through a covalent bond may prove problematic however this is sometimes unavoidable. When it does occur it is important that the bond of the QM atom be capped in order to prevent the appearance of bond cleavage in the QM system.^[9]

In systems where the QM/MM boundary cuts a bond three issues must be dealt with. First, the dangling bond of the QM system must be capped, this is because it is undesirable to truncate the QM system (treating the bond as if it has been cleaved will yield very unrealistic calculations). The second issue relates to polarisation, more specifically for electrostatic or polarized embedding it is important to ensure that the proximity of the MM charges near the boundary does not cause over-polarisation of the QM density. The final issue is the bonding MM terms must be carefully selected in order to prevent double counting of interactions when looking at bonds across the boundary.^[9]

Overall the goal is to obtain a good description of QM-MM interactions at the boundary between the QM and the MM system and there are three schemes by which this can be achieved.^[9]

Link atom schemes introduce an additional atomic centre (usually a hydrogen atom). This atom is not part of the real system. It is covalently bonded to the atom being described by quantum mechanics which serves to saturate its valency (by replacing the bond that has been broken).^[9]

In boundary atom schemes, the MM atom which is bonded across the boundary to a QM atom is replaced with a special boundary atom which appears in both the QM and the MM calculation. In the MM calculation, it simply behaves as an MM atom but in the QM system it mimics the electronic character of the MM atom bounded across the boundary to the QM atom.^[9]

These schemes place hybrid orbitals at the boundary and keep some of them frozen. These orbitals cap the QM region and replace the cut bond.^[9]

BuRNN (Buffer Region Neural Network) approach was developed as an alternative to QM/MM methods. Its focus is to reduce artifacts that are created in between QM and MM region by introducing buffer region between them. Buffer region experiences full electronic polarization by the QM region and together with QM region is described by NN (neural network) trained on QM calculations. The substitution of QM calculations for NN also speeds up overall simulation. BuRNN was introduced in the 2022 paper of Lier, Poliak, Marquetand, Westermayr, and Oostenbrink.^[10]

• ONIOM: "our own n-layered integrated molecular orbital and molecular mechanics"
• List of quantum chemistry and solid state physics software
• List of software for molecular mechanics modeling