# Molecular mechanics

A force field is used to minimize the bond stretching energy of this ethane molecule.

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.

All-atomistic molecular mechanics methods have the following properties:

• Each atom is simulated as one particle
• Each particle is assigned a radius (typically the van der Waals radius), polarizability, and a constant net charge (generally derived from quantum calculations and/or experiment)
• Bonded interactions are treated as springs with an equilibrium distance equal to the experimental or calculated bond length

Variants on this theme are possible. For example, many simulations have historically used a united-atom representation in which each terminal methyl group or intermediate methylene unit was considered one particle, and large protein systems are commonly simulated using a bead model that assigns two to four particles per amino acid.

## Functional form

Molecular mechanics potential energy function with continuum solvent.

The following functional abstraction, termed an interatomic potential function or force field in chemistry, calculates the molecular system's potential energy (E) in a given conformation as a sum of individual energy terms.

${\displaystyle \ E=E_{\text{covalent}}+E_{\text{noncovalent}}\,}$

where the components of the covalent and noncovalent contributions are given by the following summations:

${\displaystyle \ E_{\text{covalent}}=E_{\text{bond}}+E_{\text{angle}}+E_{\text{dihedral}}}$

${\displaystyle \ E_{\text{noncovalent}}=E_{\text{electrostatic}}+E_{\text{van der Waals}}}$

The exact functional form of the potential function, or force field, depends on the particular simulation program being used. Generally the bond and angle terms are modeled as harmonic potentials centered around equilibrium bond-length values derived from experiment or theoretical calculations of electronic structure performed with software which does ab-initio type calculations such as Gaussian. For accurate reproduction of vibrational spectra, the Morse potential can be used instead, at computational cost. The dihedral or torsional terms typically have multiple minima and thus cannot be modeled as harmonic oscillators, though their specific functional form varies with the implementation. This class of terms may include improper dihedral terms, which function as correction factors for out-of-plane deviations (for example, they can be used to keep benzene rings planar, or correct geometry and chirality of tetrahedral atoms in a united-atom representation).

The non-bonded terms are much more computationally costly to calculate in full, since a typical atom is bonded to only a few of its neighbors, but interacts with every other atom in the molecule. Fortunately the van der Waals term falls off rapidly. It is typically modeled using a 6–12 Lennard-Jones potential, which means that attractive forces fall off with distance as r−6 and repulsive forces as r−12, where r represents the distance between two atoms. The repulsive part r−12 is however unphysical, because repulsion increases exponentially. Description of van der Waals forces by the Lennard-Jones 6–12 potential introduces inaccuracies, which become significant at short distances.[1] Generally a cutoff radius is used to speed up the calculation so that atom pairs which distances are greater than the cutoff have a van der Waals interaction energy of zero.

The electrostatic terms are notoriously difficult to calculate well because they do not fall off rapidly with distance, and long-range electrostatic interactions are often important features of the system under study (especially for proteins). The basic functional form is the Coulomb potential, which only falls off as r−1. A variety of methods are used to address this problem, the simplest being a cutoff radius similar to that used for the van der Waals terms. However, this introduces a sharp discontinuity between atoms inside and atoms outside the radius. Switching or scaling functions that modulate the apparent electrostatic energy are somewhat more accurate methods that multiply the calculated energy by a smoothly varying scaling factor from 0 to 1 at the outer and inner cutoff radii. Other more sophisticated but computationally intensive methods are particle mesh Ewald (PME) and the multipole algorithm.

In addition to the functional form of each energy term, a useful energy function must be assigned parameters for force constants, van der Waals multipliers, and other constant terms. These terms, together with the equilibrium bond, angle, and dihedral values, partial charge values, atomic masses and radii, and energy function definitions, are collectively termed a force field. Parameterization is typically done through agreement with experimental values and theoretical calculations results. Allingers force field in the last MM4 version calculate for hydrocarbons heats of formation with a rms error of 0.35 kcal/mol; vibrational spectra with a rms error of 24 cm-1 , rotational barriers with a rms error of 2.2 %, C-C bond lengths within 0.004 o and C-C-C angles within 1o [2]; later MM4 versions cover also compounds with heteroatoms , such as aliphatic amines [3]

Each force field is parameterized to be internally consistent, but the parameters are generally not transferable from one force field to another.

## Areas of application

The main use of molecular mechanics is in the field of molecular dynamics. This uses the force field to calculate the forces acting on each particle and a suitable integrator to model the dynamics of the particles and predict trajectories. Given enough sampling and subject to the ergodic hypothesis, molecular dynamics trajectories can be used to estimate thermodynamic parameters of a system or probe kinetic properties, such as reaction rates and mechanisms.

Another application of molecular mechanics is energy minimization, whereby the force field is used as an optimization criterion. This method uses an appropriate algorithm (e.g. steepest descent) to find the molecular structure of a local energy minimum. These minima correspond to stable conformers of the molecule (in the chosen force field) and molecular motion can be modelled as vibrations around and interconversions between these stable conformers. It is thus common to find local energy minimization methods combined with global energy optimization, to find the global energy minimum (and other low energy states). At finite temperature, the molecule spends most of its time in these low-lying states, which thus dominate the molecular properties. Global optimization can be accomplished using simulated annealing, the Metropolis algorithm and other Monte Carlo methods, or using different deterministic methods of discrete or continuous optimization. While the force field represents only the enthalpic component of free energy (and only this component is included during energy minimization), it is possible to include the entropic component through the use of additional methods, such as normal mode analysis.

Molecular mechanics potential energy functions have been used to calculate binding constants,[4][5][6][7][8] protein folding kinetics,[9] protonation equilibria,[10] active site coordinates,[6][11] and to design binding sites.[12]

## Environment and solvation

In molecular mechanics, several ways exist to define the environment surrounding a molecule or molecules of interest. A system can be simulated in vacuum (termed a gas-phase simulation) with no surrounding environment, but this is usually undesirable because it introduces artifacts in the molecular geometry, especially in charged molecules. Surface charges that would ordinarily interact with solvent molecules instead interact with each other, producing molecular conformations that are unlikely to be present in any other environment. The best way to solvate a system is to place explicit water molecules in the simulation box with the molecules of interest and treat the water molecules as interacting particles like those in the molecule. A variety of water models exist with increasing levels of complexity, representing water as a simple hard sphere (a united-atom model), as three separate particles with fixed bond angles, or even as four or five separate interaction centers to account for unpaired electrons on the oxygen atom. As water models grow more complex, related simulations grow more computationally intensive. A compromise method has been found in implicit solvation, which replaces the explicitly represented water molecules with a mathematical expression that reproduces the average behavior of water molecules (or other solvents such as lipids). This method is useful to prevent artifacts that arise from vacuum simulations and reproduces bulk solvent properties well, but cannot reproduce situations in which individual water molecules have interesting interactions with the molecules under study.

## Software packages

This is a limited list; many more packages are available.

## References

1. ^ Zgarbova M, et al. (2010). "Large-scale compensation of errors in pairwise-additive empirical force fields: comparison of AMBER intermolecular terms with rigorous DFT-SAPT calculations". Phys. Chem. Chem. Phys. 12 (35): 10476–10493. Bibcode:2010PCCP...1210476Z. doi:10.1039/C002656E. PMID 20603660.
2. ^ Allinger, N. L.; Chen, K.; Lii, J.-H. J. Comput. Chem. 1996, 17, 642 https://onlinelibrary.wiley.com/doi/abs/10.1002/%28SICI%291096-987X%28199604%2917%3A5/6%3C642%3A%3AAID-JCC6%3E3.0.CO%3B2-U
3. ^ Kuo‐Hsiang Chen ,Jenn‐Huei Lii, Yi Fan, Norman L. Allinger J. Comput. Chem. 2007, 28, 2391 https://onlinelibrary.wiley.com/doi/full/10.1002/jcc.20737
4. ^ Kuhn B, Kollman PA (October 2000). "Binding of a diverse set of ligands to avidin and streptavidin: an accurate quantitative prediction of their relative affinities by a combination of molecular mechanics and continuum solvent models". Journal of Medicinal Chemistry. 43 (20): 3786–91. doi:10.1021/jm000241h. PMID 11020294.
5. ^ Huo S, Massova I, Kollman PA (January 2002). "Computational alanine scanning of the 1:1 human growth hormone-receptor complex". J Comput Chem. 23 (1): 15–27. doi:10.1002/jcc.1153. PMID 11913381.
6. ^ a b Mobley DL, Graves AP, Chodera JD, McReynolds AC, Shoichet BK, Dill KA (August 2007). "Predicting absolute ligand binding free energies to a simple model site". J Mol Biol. 371 (4): 1118–34. doi:10.1016/j.jmb.2007.06.002. PMC 2104542. PMID 17599350.
7. ^ Wang J, Kang X, Kuntz ID, Kollman PA (April 2005). "Hierarchical database screenings for HIV-1 reverse transcriptase using a pharmacophore model, rigid docking, solvation docking, and MM-PB/SA". Journal of Medicinal Chemistry. 48 (7): 2432–44. doi:10.1021/jm049606e. PMID 15801834.
8. ^ Kollman PA, Massova I, Reyes C, et al. (December 2000). "Calculating structures and free energies of complex molecules: combining molecular mechanics and continuum models". Acc Chem Res. 33 (12): 889–97. CiteSeerX 10.1.1.469.844. doi:10.1021/ar000033j. PMID 11123888.
9. ^ Snow CD, Nguyen H, Pande VS, Gruebele M (November 2002). "Absolute comparison of simulated and experimental protein-folding dynamics". Nature. 420 (6911): 102–6. Bibcode:2002Natur.420..102S. doi:10.1038/nature01160. PMID 12422224.
10. ^ Barth P, Alber T, Harbury PB (March 2007). "Accurate, conformation-dependent predictions of solvent effects on protein ionization constants". Proc Natl Acad Sci USA. 104 (12): 4898–903. Bibcode:2007PNAS..104.4898B. doi:10.1073/pnas.0700188104. PMC 1829236. PMID 17360348.
11. ^ Chakrabarti R, Klibanov AM, Friesner RA (July 2005). "Computational prediction of native protein ligand-binding and enzyme active site sequences". Proc Natl Acad Sci USA. 102 (29): 10153–8. Bibcode:2005PNAS..10210153C. doi:10.1073/pnas.0504023102. PMC 1177389. PMID 15998733.
12. ^ Boas FE, Harbury PB (July 2008). "Design of Protein-Ligand Binding Based on the Molecular-Mechanics Energy Model". J Mol Biol. 380 (2): 415–24. doi:10.1016/j.jmb.2008.04.001. PMC 2569001. PMID 18514737.
13. ^ ACEMD - GPU MD
14. ^ Ascalaph
15. ^ COSMOS
16. ^ CytoSolve
17. ^ StruMM3D (STR3DI32)
18. ^ Zodiac Archived 2009-12-16 at the Wayback Machine