View source for QM/MM

{{Short description|Molecular simulation method}}
{{cs1 config|name-list-style=vanc|display-authors=6}}
The hybrid '''QM/MM''' ([[quantum mechanics]]/[[molecular mechanics]]) approach is a molecular [[simulation]] method that combines the strengths of [[Ab initio quantum chemistry methods|''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 [[Arieh Warshel|Warshel]] and [[Michael Levitt (biophysicist)|Levitt]].<ref>{{cite journal | vauthors = Warshel A, Levitt M | title = Theoretical studies of enzymic reactions: dielectric, electrostatic and steric stabilization of the carbonium ion in the reaction of lysozyme | journal = Journal of Molecular Biology | volume = 103 | issue = 2 | pages = 227–49 | date = May 1976 | pmid = 985660 | doi = 10.1016/0022-2836(76)90311-9 }}</ref>  They, along with [[Martin Karplus]], won the 2013 [[Nobel Prize in Chemistry]] for "the development of multiscale models for complex chemical systems".<ref>{{cite press release
 | title = The Nobel Prize in Chemistry 2013
 | publisher = Royal Swedish Academy of Sciences
 | date = October 9, 2013
 | url = https://www.nobelprize.org/nobel_prizes/chemistry/laureates/2013/press.pdf
 | access-date = October 9, 2013}}</ref><ref name="NYT-20131009">{{cite news | vauthors = Chang K |title=3 Researchers Win Nobel Prize in Chemistry |url=https://www.nytimes.com/2013/10/10/science/three-researchers-win-nobel-prize-in-chemistry.html |date=October 9, 2013 |work=[[New York Times]] |access-date=October 9, 2013 }}</ref>

== Efficiency ==
An important advantage of QM/MM methods is their efficiency. The cost of doing classical molecular mechanics (MM) simulations in the most straightforward case [[Big O notation|scales]] as ''O''(''N''<sup>2</sup>), 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''<sup>2</sup>). 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''<sup>3</sup>) or worse (restricted [[Hartree–Fock method|Hartree–Fock calculations]] have been suggested to scale ~''O''(''N''<sup>2.7</sup>)). Here in the ''ab initio'' calculations, ''N'' stands for the number of [[Basis set (chemistry)|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.<ref>{{cite journal | vauthors = Brunk E, Rothlisberger U | title = Mixed Quantum Mechanical/Molecular Mechanical Molecular Dynamics Simulations of Biological Systems in Ground and Electronically Excited States | journal = Chemical Reviews | volume = 115 | issue = 12 | pages = 6217–63 | date = June 2015 | pmid = 25880693 | doi = 10.1021/cr500628b | url = http://infoscience.epfl.ch/record/207721 }}</ref><ref>{{cite journal | vauthors = Morzan UN, Alonso de Armiño DJ, Foglia NO, Ramírez F, González Lebrero MC, Scherlis DA, Estrin DA | title = Spectroscopy in Complex Environments from QM-MM Simulations | journal = Chemical Reviews | volume = 118 | issue = 7 | pages = 4071–4113 | date = April 2018 | pmid = 29561145 | doi = 10.1021/acs.chemrev.8b00026 | hdl = 11336/89583 | hdl-access = free }}</ref>

== Calculating the energy of the combined system ==
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 (chemistry)|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.

<math>E = E^{QM}(QM)+ E^{MM}(QM+MM)-E^{MM}(QM)</math>

In this equation <math>E^{MM}(QM)</math> 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:
<math>
\begin{aligned}
E(QM/MM) & = \sum^{MM}_{I'=1}\left[\int \operatorname{d}\mathbf{r}{q_{I'} \rho(\mathbf{r}) \over \left \vert \mathbf{R}_{I'}-\mathbf{r}  \right \vert }+\sum^{QM}_{I=1}\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^{12}_{II'}} - \frac{B_{II'}}{R^{6}_{II'}}\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}</math>

The index <math>I </math> labels the nuclei in the QM region whereas <math>I' </math> 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 [[London dispersion force|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 <math>\theta </math> or <math>\phi </math> will be a QM atom with the others being MM atoms.<ref name = "Allen_2017" />{{rp|422–3}}

=== Reducing the computational cost of calculating QM-MM interactions ===
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<sup>6</sup> grid points for the electron density of the QM system and 10<sup>4</sup> MM atoms). A method by which this issue can be mitigated is to construct three [[Concentric objects|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 <math>E(QM/MM) </math>. 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.<ref name = "Allen_2017">{{cite book |title=Computer Simulation of Liquids | vauthors = Allen MP, Tildesley DJ | date= August 2017 |publisher=Oxford University Press|isbn=9780198803201|edition= 2nd|location=Oxford }}</ref>{{rp|423–4}}

== The electrostatic QM-MM interaction ==
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 ===
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.<ref name="Lin_2017">{{Cite journal| vauthors = Lin H, Truhlar DG |title=QM/MM: what have we learned, where are we, and where do we go from here? |journal=Theoretical Chemistry Accounts | date = February 2007 | volume=117 |issue=2 |pages=185–199 |doi=10.1007/s00214-006-0143-z |doi-access=free }}</ref>

=== Electronic embedding ===
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 (quantum mechanics)|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.<ref name="Lin_2017" />

=== Polarized embedding ===

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.<ref name="Lin_2017" />

== Problems involved with QM/MM ==
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.<ref name="Kerdcharoen1996">{{cite journal | vauthors = Kerdcharoen T, Liedl KR, Rode BM | title = A QM/MM simulation method applied to the solution of Li+ in liquid ammonia | journal = Chemical Physics | volume = 211 | issue = 1 | pages = 313–323 | date = 1996 | doi = 10.1016/0301-0104(96)00152-8 | bibcode = 1996CP....211..313K }}</ref> 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.<ref name="Senn_2009">{{cite journal | vauthors = Senn HM, Thiel W | title = QM/MM methods for biomolecular systems | journal = Angewandte Chemie | volume = 48 | issue = 7 | pages = 1198–229 | date = 2009 | pmid =  19173328| doi = 10.1002/anie.200802019 }}</ref>

== Covalent bonds across the QM-MM boundary ==
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.<ref name="Senn_2009" />

=== Boundary schemes ===
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.<ref name="Senn_2009" />

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.<ref name="Senn_2009" />

==== Link atom schemes ====
Link atom schemes introduce an additional atomic centre (usually a hydrogen atom). This atom is not part of the real system. It is [[Covalent bond|covalently bonded]] to the atom being described by quantum mechanics which serves to saturate its valency (by replacing the bond that has been broken).<ref name="Senn_2009" />

==== Boundary atom schemes ====
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.<ref name="Senn_2009" />

==== Localized-orbital schemes ====
These schemes place [[Orbital hybridisation|hybrid orbitals]] at the boundary and keep some of them frozen. These orbitals cap the QM region and replace the cut bond.<ref name="Senn_2009" />

== BuRNN ==
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.<ref>{{cite journal | vauthors = Lier B, Poliak P, Marquetand P, Westermayr J, Oostenbrink C | title = BuRNN: Buffer Region Neural Network Approach for Polarizable-Embedding Neural Network/Molecular Mechanics Simulations | journal = The Journal of Physical Chemistry Letters | volume = 13 | issue = 17 | pages = 3812–3818 | date = May 2022 | pmid = 35467875 | pmc = 9082612 | doi = 10.1021/acs.jpclett.2c00654 }}</ref>

== See also ==
* [[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]]

== References ==
{{Reflist}}

{{DEFAULTSORT:Qm Mm}}
[[Category:Computational chemistry]]
[[Category:Molecular dynamics]]