QM/MM

Last updated

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]

Contents

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 scales as O(N2), 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(N2). 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(N3) or worse (restricted Hartree–Fock calculations have been suggested to scale ~O(N2.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]

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, then the energy of the QM system is added (calculated using a QM method), finally the MM energy of the QM system is subtracted.

In this equation 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:

The index labels the nuclei in the QM region whereas 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 or will be a QM atom with the others being MM atoms. [6] :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 106 grid points for the electron density of the QM system and 104 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 . 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

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. [7]

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. 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]

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. [7]

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. [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]

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. [9]

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. [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]

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. [9]

Localized-orbital schemes

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]

See also

Related Research Articles

<span class="mw-page-title-main">Dipole</span> Electromagnetic phenomenon

In physics, a dipole is an electromagnetic phenomenon which occurs in two ways:

An intermolecular force (IMF) is the force that mediates interaction between molecules, including the electromagnetic forces of attraction or repulsion which act between atoms and other types of neighbouring particles, e.g. atoms or ions. Intermolecular forces are weak relative to intramolecular forces – the forces which hold a molecule together. For example, the covalent bond, involving sharing electron pairs between atoms, is much stronger than the forces present between neighboring molecules. Both sets of forces are essential parts of force fields frequently used in molecular mechanics.

In physics, screening is the damping of electric fields caused by the presence of mobile charge carriers. It is an important part of the behavior of charge-carrying fluids, such as ionized gases, electrolytes, and charge carriers in electronic conductors . In a fluid, with a given permittivity ε, composed of electrically charged constituent particles, each pair of particles interact through the Coulomb force as

<span class="mw-page-title-main">Molecular dynamics</span> Computer simulations to discover and understand chemical properties

Molecular dynamics (MD) is a computer simulation method for analyzing the physical movements of atoms and molecules. The atoms and molecules are allowed to interact for a fixed period of time, giving a view of the dynamic "evolution" of the system. In the most common version, the trajectories of atoms and molecules are determined by numerically solving Newton's equations of motion for a system of interacting particles, where forces between the particles and their potential energies are often calculated using interatomic potentials or molecular mechanical force fields. The method is applied mostly in chemical physics, materials science, and biophysics.

<span class="mw-page-title-main">AMBER</span>

Assisted Model Building with Energy Refinement (AMBER) is a family of force fields for molecular dynamics of biomolecules originally developed by Peter Kollman's group at the University of California, San Francisco.

Chemistry at Harvard Macromolecular Mechanics (CHARMM) is the name of a widely used set of force fields for molecular dynamics, and the name for the molecular dynamics simulation and analysis computer software package associated with them. The CHARMM Development Project involves a worldwide network of developers working with Martin Karplus and his group at Harvard to develop and maintain the CHARMM program. Licenses for this software are available, for a fee, to people and groups working in academia.

<span class="mw-page-title-main">Molecular mechanics</span> Use of classical mechanics to model molecular systems

Molecular mechanics uses classical mechanics to model molecular systems. The Born–Oppenheimer approximation is assumed valid and the potential energy of all systems is calculated as a function of the nuclear coordinates using force fields. Molecular mechanics can be used to study molecule systems ranging in size and complexity from small to large biological systems or material assemblies with many thousands to millions of atoms.

<span class="mw-page-title-main">Molecular modelling</span> Discovering chemical properties by physical simulations

Molecular modelling encompasses all methods, theoretical and computational, used to model or mimic the behaviour of molecules. The methods are used in the fields of computational chemistry, drug design, computational biology and materials science to study molecular systems ranging from small chemical systems to large biological molecules and material assemblies. The simplest calculations can be performed by hand, but inevitably computers are required to perform molecular modelling of any reasonably sized system. The common feature of molecular modelling methods is the atomistic level description of the molecular systems. This may include treating atoms as the smallest individual unit, or explicitly modelling protons and neutrons with its quarks, anti-quarks and gluons and electrons with its photons.

Mulliken charges arise from the Mulliken population analysis and provide a means of estimating partial atomic charges from calculations carried out by the methods of computational chemistry, particularly those based on the linear combination of atomic orbitals molecular orbital method, and are routinely used as variables in linear regression (QSAR) procedures. The method was developed by Robert S. Mulliken, after whom the method is named. If the coefficients of the basis functions in the molecular orbital are Cμi for the μ'th basis function in the i'th molecular orbital, the density matrix terms are:

<span class="mw-page-title-main">Force field (chemistry)</span> Concept on molecular modeling

In the context of chemistry and molecular modelling, a force field is a computational method that is used to estimate the forces between atoms within molecules and also between molecules. More precisely, the force field refers to the functional form and parameter sets used to calculate the potential energy of a system of atoms or coarse-grained particles in molecular mechanics, molecular dynamics, or Monte Carlo simulations. The parameters for a chosen energy function may be derived from experiments in physics and chemistry, calculations in quantum mechanics, or both. Force fields are interatomic potentials and utilize the same concept as force fields in classical physics, with the difference that the force field parameters in chemistry describe the energy landscape, from which the acting forces on every particle are derived as a gradient of the potential energy with respect to the particle coordinates.

Ewald summation, named after Paul Peter Ewald, is a method for computing long-range interactions in periodic systems. It was first developed as the method for calculating the electrostatic energies of ionic crystals, and is now commonly used for calculating long-range interactions in computational chemistry. Ewald summation is a special case of the Poisson summation formula, replacing the summation of interaction energies in real space with an equivalent summation in Fourier space. In this method, the long-range interaction is divided into two parts: a short-range contribution, and a long-range contribution which does not have a singularity. The short-range contribution is calculated in real space, whereas the long-range contribution is calculated using a Fourier transform. The advantage of this method is the rapid convergence of the energy compared with that of a direct summation. This means that the method has high accuracy and reasonable speed when computing long-range interactions, and it is thus the de facto standard method for calculating long-range interactions in periodic systems. The method requires charge neutrality of the molecular system to accurately calculate the total Coulombic interaction. A study of the truncation errors introduced in the energy and force calculations of disordered point-charge systems is provided by Kolafa and Perram.

Car–Parrinello molecular dynamics or CPMD refers to either a method used in molecular dynamics or the computational chemistry software package used to implement this method.

Implicit solvation is a method to represent solvent as a continuous medium instead of individual “explicit” solvent molecules, most often used in molecular dynamics simulations and in other applications of molecular mechanics. The method is often applied to estimate free energy of solute-solvent interactions in structural and chemical processes, such as folding or conformational transitions of proteins, DNA, RNA, and polysaccharides, association of biological macromolecules with ligands, or transport of drugs across biological membranes.

<span class="mw-page-title-main">Water model</span> Aspect of computational chemistry

In computational chemistry, a water model is used to simulate and thermodynamically calculate water clusters, liquid water, and aqueous solutions with explicit solvent. The models are determined from quantum mechanics, molecular mechanics, experimental results, and these combinations. To imitate a specific nature of molecules, many types of models have been developed. In general, these can be classified by the following three points; (i) the number of interaction points called site, (ii) whether the model is rigid or flexible, (iii) whether the model includes polarization effects.

Drude particles are model oscillators used to simulate the effects of electronic polarizability in the context of a classical molecular mechanics force field. They are inspired by the Drude model of mobile electrons and are used in the computational study of proteins, nucleic acids, and other biomolecules.

Localized molecular orbitals are molecular orbitals which are concentrated in a limited spatial region of a molecule, such as a specific bond or lone pair on a specific atom. They can be used to relate molecular orbital calculations to simple bonding theories, and also to speed up post-Hartree–Fock electronic structure calculations by taking advantage of the local nature of electron correlation. Localized orbitals in systems with periodic boundary conditions are known as Wannier functions.

Free energy perturbation (FEP) is a method based on statistical mechanics that is used in computational chemistry for computing free energy differences from molecular dynamics or Metropolis Monte Carlo simulations.

Biology Monte Carlo methods (BioMOCA) have been developed at the University of Illinois at Urbana-Champaign to simulate ion transport in an electrolyte environment through ion channels or nano-pores embedded in membranes. It is a 3-D particle-based Monte Carlo simulator for analyzing and studying the ion transport problem in ion channel systems or similar nanopores in wet/biological environments. The system simulated consists of a protein forming an ion channel (or an artificial nanopores like a Carbon Nano Tube, CNT), with a membrane (i.e. lipid bilayer) that separates two ion baths on either side. BioMOCA is based on two methodologies, namely the Boltzmann transport Monte Carlo (BTMC) and particle-particle-particle-mesh (P3M). The first one uses Monte Carlo method to solve the Boltzmann equation, while the later splits the electrostatic forces into short-range and long-range components.

<span class="mw-page-title-main">Interatomic potential</span> Functions for calculating potential energy

Interatomic potentials are mathematical functions to calculate the potential energy of a system of atoms with given positions in space. Interatomic potentials are widely used as the physical basis of molecular mechanics and molecular dynamics simulations in computational chemistry, computational physics and computational materials science to explain and predict materials properties. Examples of quantitative properties and qualitative phenomena that are explored with interatomic potentials include lattice parameters, surface energies, interfacial energies, adsorption, cohesion, thermal expansion, and elastic and plastic material behavior, as well as chemical reactions.

In computational chemistry, a solvent model is a computational method that accounts for the behavior of solvated condensed phases. Solvent models enable simulations and thermodynamic calculations applicable to reactions and processes which take place in solution. These include biological, chemical and environmental processes. Such calculations can lead to new predictions about the physical processes occurring by improved understanding.

References

  1. Warshel A, Levitt M (May 1976). "Theoretical studies of enzymic reactions: dielectric, electrostatic and steric stabilization of the carbonium ion in the reaction of lysozyme". Journal of Molecular Biology. 103 (2): 227–49. doi:10.1016/0022-2836(76)90311-9. PMID   985660.
  2. "The Nobel Prize in Chemistry 2013" (PDF) (Press release). Royal Swedish Academy of Sciences. October 9, 2013. Retrieved October 9, 2013.
  3. Chang K (October 9, 2013). "3 Researchers Win Nobel Prize in Chemistry". New York Times . Retrieved October 9, 2013.
  4. Brunk E, Rothlisberger U (June 2015). "Mixed Quantum Mechanical/Molecular Mechanical Molecular Dynamics Simulations of Biological Systems in Ground and Electronically Excited States". Chemical Reviews. 115 (12): 6217–63. doi:10.1021/cr500628b. PMID   25880693.
  5. Morzan UN, Alonso de Armiño DJ, Foglia NO, Ramírez F, González Lebrero MC, Scherlis DA, Estrin DA (April 2018). "Spectroscopy in Complex Environments from QM-MM Simulations". Chemical Reviews. 118 (7): 4071–4113. doi:10.1021/acs.chemrev.8b00026. PMID   29561145.
  6. 1 2 Allen MP, Tildesley DJ (August 2017). Computer Simulation of Liquids (2nd ed.). Oxford: Oxford University Press. ISBN   9780198803201.
  7. 1 2 3 Lin H, Truhlar DG (February 2007). "QM/MM: what have we learned, where are we, and where do we go from here?". Theoretical Chemistry Accounts. 117 (2): 185–199. doi: 10.1007/s00214-006-0143-z .
  8. Kerdcharoen T, Liedl KR, Rode BM (1996). "A QM/MM simulation method applied to the solution of Li+ in liquid ammonia". Chemical Physics. 211 (1): 313–323. Bibcode:1996CP....211..313K. doi:10.1016/0301-0104(96)00152-8.
  9. 1 2 3 4 5 6 7 Senn HM, Thiel W (2009). "QM/MM methods for biomolecular systems". Angewandte Chemie. 48 (7): 1198–229. doi:10.1002/anie.200802019. PMID   19173328.