Molecular Simulation/Molecular Mechanical Force Fields

General Form of Potential Energy Function

The overall energy of a system is based on an enormous number of factors, and it is the sum of these individual attractive and repulsive forces between all the atoms in the system that results in a final "force field". The potential energy for a set of atomic coordinates (r) is the sum of bonded (ie intermolecular forces) and non-bonded (i.e., repulsion/dispersion and electrostatics) terms. This sum ultimately defines the molecular mechanical force field.

${\mathcal {V}}({\textrm {r}})=\sum _{bonds}{\frac {1}{2}}k_{bond}\left(r-r_{e}\right)^{2}+\sum _{angles}{\frac {1}{2}}k_{angle}\left(\theta -\theta _{e}\right)^{2}+\sum _{dihedrals}{k_{\phi }}({1+\cos(n\phi -\delta }))+\sum _{pairs:i,j}{4\epsilon _{ij}}\left[\left({\frac {\sigma _{ij}}{r_{ij}}}\right)^{12}-\left({\frac {\sigma _{ij}}{r_{ij}}}\right)^{6}\right]+{\frac {q_{1}q_{2}}{4\pi \epsilon _{0}r_{ij}}}$

Due to this complexity, calculating the potential energy for a set of molecules using a force field would require a complete set of known parameters, which also require several different methods of determination. These methods are outlined in this section:

Interaction	Term type	Expression	Parameters
Bonds	Bonded	${\frac {1}{2}}k_{bond}\left(r-r_{e}\right)^{2}$	spring constant (kbond), equilibrium length (re)
Angles	Bonded	${\frac {1}{2}}k_{angle}\left(\theta -\theta _{e}\right)^{2}$	spring constant (kangle), equilibrium angle (θe)
Dihedrals	Bonded	${k_{\phi }}({1+\cos(n\phi -\delta }))$	barrier height (kφ), periodicity (n), offset (δ)
Repulsion/Dispersion	Non-bonded	${4\epsilon _{ij}}\left[\left({\frac {\sigma _{ij}}{r_{ij}}}\right)^{12}-\left({\frac {\sigma _{ij}}{r_{ij}}}\right)^{6}\right]$	atomic radii (σ), well depth (ε)
Electrostatics	Non-bonded	${\frac {q_{1}q_{2}}{4\pi \epsilon _{0}r_{ij}}}$	partial atomic charges (q)

Parameterization

Empirical Parameters

For small molecules, parameters may be defined to fit experimental properties. Liquids will have different bulk physical properties for each set of non-bonded parameters. Each set of parameters in the force field will result in a different set of physical properties, the combinations of which can be tested until properly representing the properties of the liquid. An example would be determining the proper model for bulk water assuming no Lennard-Jones interactions for hydrogen atoms. Parameters would be tested together until the model properly represented the expected characteristics of bulk water.

Calculating Electrostatic Interactions

Starting with the non-bonded terms, we can look at electrostatics. Every intermolecular pair of atoms has this electrostatic interaction, and can be summed up in a form comparing all atoms in one molecule (A) to all atoms in another molecule (B). This quickly adds up to a massive sum even when considering systems of a few complex molecules:

${\mathcal {V}}(r_{ij})=\sum _{i}^{N_{A}}\sum _{j}^{N_{B}}{\frac {q_{1}q_{2}}{4\pi \epsilon _{0}r_{ij}}}$

In regards to partial charges, as there is no rigorous way to translate charge distribution to point charges, certain constraints are formed: sum of partial charges should equal formal charge of molecule (seen below), and equivalent atoms should carry the same charge. Charges can be chosen to match gas phase dipole moment, be treated as free parameters and adjusted for target properties, and to fit quantum mechanical electrostatics:

${\vec {\mu }}_{exptl}=\sum _{i}^{N_{atom}}q_{i}{\vec {r}}_{i}$

Charge fitting schemes fit partial charges to the quantum mechanical electrostatic potential (the potential energy a unit point charge has at position in space r). Potential can be calculated at a set of points around the molecule to generate electrostatic potential surface (ie: can map out the electronic density around the molecule and determine the localized area of charges):

Target Electrostatic Potential Data from Quantum Chemical Calculations

First term represents electron density, second term represents nuclear charges at positions R_i:

${\mathcal {V}}_{QM}(r)=\int {\frac {\rho (r')}{|r-r'|}}dr'+\sum _{i}^{N}{\frac {Z_{i}}{|r-R_{i}|}}$

Molecular Mechanical Partial Atomic Charges

Looks at partial atomic charges at positions R_i:

${\mathcal {V}}_{MM}(r)=\sum _{i}^{N}{\frac {q_{i}}{|r-R_{i}|}}$

In terms of electrostatic charges, point charges (q) must be defined as such that $V_{QM}(r)$ ≈ $V_{MM}(r)$ as closely as possible, thereby fitting both models.

Restraining terms can be introduced to give more practical models of molecules (known as Restrained Electronic Potential (RESP) fitting). Such restraints limits charges to moderate and chemically relevant values. Restraints are also placed to ensure chemically equivalent groups retain the same charges. These groups are given an atom type distinguished by element, hybridization and bonding partners, and parameters are defined for each type. For example, all aromatic carbons part of the same benzene ring should have the same parameters.

Generalized Amber Force Field (GAFF)

GAFF is a commonly used force field designed for reasonable parameters for simulations of small organic molecules. GAFF involves the use of the RESP method for assigning partial charges, and atom type defines the Lennard-Jones parameters and internal bonding terms that will be associated with the atom. Atom type definitions can be extremely specific (e.g., H on aliphatic C with 3 electron-withdrawing groups; inner sp2 C in conjugated ring systems).

Lennard-Jones Potential

A combination of the Pauli Repulsion and London dispersion attraction terms, the Lennard-Jones Potential is defined using two parameters: van der Waals radius (σ) and well depth (ε). These parameters are responsible for packing of molecules in liquids, and must compensate for errors in all other properties. van der Waals radius essentially sets the limit at which molecules may approach one another (potential energy > 0 when closer than σ). Well depth sets the energy well of the potential energy, occurring when distance between molecules is at R_min (ie radius related to minimum potential energy). First term represents Pauli repulsion, second term represents London dispersion attraction:

${\mathcal {V}}(r)={\frac {C_{12}}{r^{12}}}-{\frac {C_{6}}{r^{6}}}$ $={4\epsilon }\left[\left({\frac {\sigma }{r}}\right)^{12}-\left({\frac {\sigma }{r}}\right)^{6}\right]$ $={\epsilon }\left[\left({\frac {R_{min}}{r}}\right)^{12}-2\left({\frac {R_{min}}{r}}\right)^{6}\right]$

Bonding Interactions

Bond Stretch

Atoms within a molecule may vary their distances (stretch or compress), and this relationship can be described as a spring. This is why Hooke's Law can be applied to bond lengths to define the energy required to stretch/compress, according to spring constant (k_bond) and equilibrium length (r_e). Spring constant may be determined via the infrared vibrational frequency of the bond stretch (according to $\nu ={\frac {1}{2\pi }}\left({\frac {k}{\mu }}\right)^{\frac {1}{2}}$, while equilibrium bond lengths can be found through the use of rotational spectroscopy and crystallography.

${\mathcal {V}}_{bond}={\frac {1}{2}}k_{bond}\left(r-r_{e}\right)^{2}$

Angle Stretch

The potential energy of bending a bond angle can also be approximated by Hooke's law. Hooke's Law can again be applied to bond angles to define the energy required to stretch/compress, according to spring constant (k_angle) and equilibrium angle (θ_e). Again, the spring constant may be determined via the infrared vibrational frequency of the bond stretch, following the same formula, while the equilibrium angle can be found through the use of rotational spectroscopy and crystallography.

${\mathcal {V}}_{angle}={\frac {1}{2}}k_{angle}\left(\theta -\theta _{e}\right)^{2}$

Dihedral Parameters

Also known as torsional rotations, dihedral relationships are produced with 4 atoms. This relates to how the central bond is rotated and the overlap of the outer two atoms (designated as 1,4 interactions, with 2 and 3 being the central atoms). These overlaps may result in a higher or lower potential energy depending on confirmation (i.e. gauche, trans, periplanar, etc), and the relationship is dependent on barrier height (k_φ) (energy required to rotate the central axis), periodicity (n) and offset (δ). Determination of these parameters can be done using quantum chemistry or experimental isomer populations.

${\mathcal {V}}_{dihedral}={k_{\phi }}({1+\cos(n\phi -\delta }))$