Potential energy functions

Evaluating the force is the most computationally demanding part of molecular dynamics. The force is the negative gradient of a scalar potential energy function,

$$\vec{F}(\vec{r}) = -\nabla U(\vec{r}),$$ (1)

and, for systems of biomolecules, this potential function involves the summing,

$$U(\vec{r}) = \sum U_{\text{bonded}}(\vec{r}) + \sum U_{\text{nonbonded}}(\vec{r}),$$ (2)

over a large number of bonded and nonbonded terms. The bonded potential terms involve 2-, 3-, and 4-body interactions of covalently bonded atoms, with $O(N)$ terms in the summation. The nonbonded potential terms involve interactions between all pairs of atoms (usually excluding pairs of atoms already involved in a bonded term), with $O(N^2)$ terms in the summation, although fast evaluation techniques are used to compute good approximations to their contribution to the potential with $O(N)$ or $O(N \log N)$ computational cost.

Bonded potential energy terms

The bonded potential terms involve 2-, 3-, and 4-body interactions of covalently bonded atoms.

The 2-body spring bond potential describes the harmonic vibrational motion between an $(i,j)$-pair of covalently bonded atoms,

$$U_{\text{bond}} = k (r_{ij} - r_0)^2,$$ (3)

where $r_{ij} = \Vert\vec{r}_j - \vec{r}_i\Vert$ gives the distance between the atoms, $r_0$ is the equilibrium distance, and $k$ is the spring constant.

The 3-body angular bond potential describes the angular vibrational motion occurring between an $(i,j,k)$-triple of covalently bonded atoms,

$$U_{\text{angle}} = k_{\theta} (\theta - \theta_0)^2 + k_{\text{ub}} (r_{ik} - r_{\text{ub}})^2,$$ (4)

where, in the first term, $\theta$ is the angle in radians between vectors $\vec{r}_{ij} = \vec{r}_j - \vec{r}_i$ and $\vec{r}_{kj} = \vec{r}_j - \vec{r}_k$, $\theta_0$ is the equilibrium angle, and $k_{\theta}$ is the angle constant. The second term is the Urey-Bradley term used to describe a (noncovalent) spring between the outer $i$ and $k$ atoms, active when constant $k_{\text{ub}} \neq 0$, where, like the spring bond, $r_{ik} = \Vert\vec{r}_k - \vec{r}_i\Vert$ gives the distance between the pair of atoms and $r_{\text{ub}}$ is the equilibrium distance.

The 4-body torsion angle (also known as dihedral angle) potential describes the angular spring between the planes formed by the first three and last three atoms of a consecutively bonded $(i,j,k,l)$-quadruple of atoms,

$$U_{\text{tors}} = \begin{cases}k (1 + \cos(n \psi + \phi)) & \text{if $n > 0$,} \\ k (\psi - \phi)^2 & \text{if $n = 0$,} \end{cases}$$ (5)

where $\psi$ is the angle in radians between the $(i,j,k)$-plane and the $(j,k,l)$-plane. The integer constant $n$ is nonnegative and indicates the periodicity. For $n > 0$, $\phi$ is the phase shift angle and $k$ is the multiplicative constant. For $n = 0$, $\phi$ acts as an equilibrium angle and the units of $k$ change to potential$/$rad$^2$. A given $(i,j,k,l)$-quadruple of atoms might contribute multiple terms to the potential, each with its own parameterization. The use of multiple terms for a torsion angle allows for complex angular variation of the potential, effectively a truncated Fourier series.

Nonbonded potential energy terms

The nonbonded potential terms involve interactions between all $(i,j)$-pairs of atoms, usually excluding pairs of atoms already involved in a bonded term. Even using a fast evaluation methods the cost of computing the nonbonded potentials dominates the work required for each time step of an MD simulation.

The Lennard-Jones potential accounts for the weak dipole attraction between distant atoms and the hard core repulsion as atoms become close,

$$U_{\text{LJ}} = (-E_{\text{min}}) \left[ \left( \frac{R_{\text{mi... ...ij}} \right)^{12} - 2 \left( \frac{R_{\text{min}}}{r_{ij}} \right)^{6} \right],$$ (6)

where $r_{ij} = \Vert\vec{r}_j - \vec{r}_i\Vert$ gives the distance between the pair of atoms. The parameter $E_{\text{min}} = U_{\text{LJ}}(R_{\text{min}})$ is the minimum of the potential term ( $E_{\text{min}} < 0$, which means that $-E_{\text{min}}$ is the well-depth). The Lennard-Jones potential approaches 0 rapidly as $r_{ij}$ increases, so it is usually truncated (smoothly shifted) to 0 past a cutoff radius, requiring $O(N)$ computational cost.

The electrostatic potential is repulsive for atomic charges with the same sign and attractive for atomic charges with opposite signs,

$$U_{\text{elec}} = \epsilon_{14} \frac{C q_i q_j}{\epsilon_0 r_{ij}},$$ (7)

where $r_{ij} = \Vert\vec{r}_j - \vec{r}_i\Vert$ gives the distance between the pair of atoms, and $q_i$ and $q_j$ are the charges on the respective atoms. Coulomb's constant $C$ and the dielectric constant $\epsilon_0$ are fixed for all electrostatic interactions. The parameter $\epsilon_{14}$ is a unitless scaling factor whose value is 1, except for a modified 1-4 interaction, where the pair of atoms is separated by a sequence of three covalent bonds (so that the atoms might also be involved in a torsion angle interaction), in which case $\epsilon_{14} = \varepsilon$, for a fixed constant $0 \leq \varepsilon \leq 1$. Although the electrostatic potential may be computed with a cutoff like the Lennard-Jones potential, the $1/r$ potential approaches 0 much more slowly than the $1/r^6$ potential, so neglecting the long range electrostatic terms can degrade qualitative results, especially for highly charged systems. There are other fast evaluation methods that approximate the contribution to the long range electrostatic terms that require $O(N)$ or $O(N \log N)$ computational cost, depending on the method.