Integration of the equations of motion#

Now that we have forces defined, we need to have a method for integrating the equations of motion. Let’s take a look of a few of them. Numerically, this is an initial value problem. There are several possible numerical integration methods, and all of the software packages have at least a few of them available. In general, there are several requirements for an algorithm, including

Requirements for numerical algorithms in MD

  • The algorithm should not violate physical laws.

  • Algorithms should always be checked that don’t violate conservation laws, for example energy, momentum (both linear and angular).

  • The algorithm must be time-reversible.

  • The algorithm should be symplectic (conservation of phase space volume).

  • The algorithm should be fast.

  • The algorithm should allow for a long time step.

  • Ideally, the algorithm should also be (relatively) easy to program.

Below we discuss several algorithms and illustrate some of the points listed above. The list of algorithms below is by no means exhaustive and many others, especially predictor-corrector type algorithms exist for MD.

../../_images/integration.svg

Fig. 9 Schematic showing integration. The system is propagate using a finite time step \(\Delta t\).#

Euler algorithm#

Let’s start the discussion from an algorithm that one should not use, the forward Euler algorithm. The algorithm is very easy to derive, one simply Taylor expands both the position and the velocity in the forward direction,

\[\begin{align*} \vec{r}_i(t + \Delta t) & = \vec{r}_i(t) + \frac{d\,\vec{r}_i(t)}{d\,t} \Delta t + \frac{1}{2}\frac{d^2\,\vec{r}_i(t)}{d\,t^2} (\Delta t)^2 \\ \vec{v}_i(t + \Delta t) & = \vec{v}_i(t) + \frac{d\,\vec{v}_i(t)}{d\,t} \Delta t + \frac{1}{2}\frac{d^2\,\vec{v}_i(t)}{d\,t^2} (\Delta t)^2 \end{align*}\]

Using \(\vec{F} = m \vec{a} = m d^2\vec{r}(t) /dt^2 \), and dropping the second order term for velocity, we can write the above equations also as

\[\begin{align*} %\left\{ \vec{r}_i(t + \Delta t) & = \vec{r}_i(t) + \vec{v}_i(t) \Delta t + \frac{\vec{F}_i(t)}{2m_i} (\Delta t)^2 \\ \vec{v}_i(t + \Delta t) & = \vec{v}_i(t) + \frac{\vec{F}_i(t)}{m_i} \Delta t % + \frac{1}{2}\frac{d^2\,\vec{v}_i(t)}{d\,t^2} (\Delta t)^2 %\right. \end{align*}\]

This also looks good, but let’s investigate one of the requirements, namely, that the algorithm is time-reversal. This is easy to do: We simply propagate the system forward in time by \(\Delta t\) and then propagate it backwards by \(-\Delta t\) and we should arrive to the starting point, that is, we should obtain

\[\begin{align*} r_\mathrm{b} (t_0) & = r_\mathrm{f} (t_0) \\ v_\mathrm{b} (t_0) & = v_\mathrm{f} (t_0). \\ \end{align*}\]

Let’s test the above for the position using the time \(t_0\) as the reference point. Let’s use one the dimensional case and drop the particle indices (\(i\)) for notational convenience. Then, using the above equations, we have

\[\begin{align*} % ok r_\mathrm{b} ( t_0 +\Delta t - \Delta t) & = r_\mathrm{f} (t_0 + \Delta t) + v(t_0 + \Delta t) (-\Delta t) + \frac{F(t_0 + \Delta t)}{2m} ( - \Delta t)^2 \\ r_\mathrm{b} (t_0) & = r_\mathrm{f} (t_0 + \Delta t) + v(t_0 + \Delta t) (-\Delta t) + \frac{F(t_0 + \Delta t)}{2m} ( - \Delta t)^2 \end{align*}\]

Substitute the expression from the above for \(r_\mathrm{f} (t_0 + \Delta t) \):

\[\begin{align*} r_\mathrm{b} (t_0) & = \left[ r_\mathrm{f}(t_0) + v(t_0) \Delta t + \frac{F(t_0)}{2m} (\Delta t)^2 \right] + v(t_0 + \Delta t) (-\Delta t) + \frac{F(t_0 + \Delta t)}{2m} ( - \Delta t)^2\\ % %r_\mathrm{fw} (t_0 + \Delta t) %+ v(t_0 + \Delta t) (-\Delta t) + \frac{F(t_0 + \Delta t)}{2m} ( - \Delta t)^2 \\ \end{align*}\]

and we do the same using \(\vec{v}_i(t + \Delta t)\),

\[\begin{align*} r_\mathrm{b} (t_0) & =r_\mathrm{f}(t_0) + \cancel{v(t_0) \Delta t} + \frac{F(t_0)}{2m} (\Delta t)^2 + \left[ \cancel{v(t_0)} + \frac{F(t_0)}{m} \Delta t \right] (-\Delta t) + \frac{F(t_0 + \Delta t)}{2m} ( - \Delta t)^2\\ r_\mathrm{b} (t_0) & =r_\mathrm{f}(t_0) + \frac{F(t_0)}{2m} (\Delta t)^2 - \frac{F(t_0)}{m} (\Delta t)^2 + \frac{F(t_0 + \Delta t)}{2m} ( - \Delta t)^2\\ r_\mathrm{b} (t_0) & =r_\mathrm{f}(t_0) + \frac{(\Delta t)^2}{2m} \left[ F(t_0 + \Delta t) - F(t_0) \right] \\ \end{align*}\]

This clearly violates the requirement that \(r_\mathrm{b} (t_0) = r_\mathrm{f} (t_0) \).

Properties of the forward Euler algorithm

  • It is not time-reversible.

  • The forward Euler algorithm doesn’t conserve energy.

  • Error (local) is \(\mathcal{O}(\Delta t^2)\)

  • Easy to derive.

  • Illustrates well some problems with algorithms.

  • Despite the problems, it is available in several software packages.

Verlet algorithm#

Unlike that forward Euler algorithm, the Verlet algorithm is well-behaved. The Verlet algorithm in MD is due to Loup Verlet, but the general algorithm goes back to 1791. The Verlet algorithm us not commonly used in MD simulations since in its basic forms, it propagates positions only. In it’s essence, the Verlet algorithm is central difference approximation of the second derivative.

The Verlet algorithm is obtained by addition of Taylor expanding the position \(\vec{r}(t)\) twice and then adding the expansions.

\[\begin{align*} \vec{r}(t + \Delta t) & = \vec{r}(t) + \frac{d\,\vec{r}(t)}{d\,t} \Delta t + \frac{1}{2}\frac{d^2\,\vec{r}(t)}{d\,t^2} (\Delta t)^2 + \frac{1}{3!}\frac{d^3\,\vec{r}(t)}{d\,t^3} (\Delta t)^3 + \mathcal{O} ((\Delta t)^4) \\ \vec{r}(t - \Delta t) & = \vec{r}(t) - \frac{d\,\vec{r}(t)}{d\,t} \Delta t + \frac{1}{2}\frac{d^2\,\vec{r}(t)}{d\,t^2} (\Delta t)^2 - \frac{1}{3!}\frac{d^3\,\vec{r}(t)}{d\,t^3} (\Delta t)^3 + \mathcal{O} ((\Delta t)^4) %\dot{\vec{r}}(t) \Delta t + \end{align*}\]

If we now add the two expansions, we obtain

\[\begin{align*} \vec{r}(t + \Delta t) + \vec{r}(t - \Delta t) & = 2 \vec{r}(t) + \frac{d^2\,\vec{r}(t)}{d\,t^2} (\Delta t)^2 + \mathcal{O} ((\Delta t)^4) \\ \vec{r}(t + \Delta t) & = 2 \vec{r}(t) - \vec{r}(t - \Delta t) + \frac{d^2\,\vec{r}(t)}{d\,t^2} (\Delta t)^2 + \mathcal{O} ((\Delta t)^4) \\ \vec{r}(t + \Delta t) & = 2 \vec{r}(t) - \vec{r}(t - \Delta t) + \frac{\vec{F}}{m} (\Delta t)^2 + \mathcal{O} ((\Delta t)^4) \end{align*}\]

The last line was obtained by noticing that \(\vec{F} = m \vec{a} = m d^2\vec{r}(t) /dt^2 \). The last term, \(\mathcal{O} ((\Delta t)^4) \) means that the error is of order \((\Delta t)^4\).

To resolve the velocities that are needed for kinetic energy, one simply uses the central difference formula

(8)#\[\begin{equation} \vec{v}_{i}(t)=\frac{\vec{r}_{i}(t+\Delta{t})-\vec{r}_{i}(t-\Delta{t})}{2\Delta{t}}. \end{equation}\]

Note that the accuracy of velocity is not \((\Delta t)^4\) but \((\Delta t)^2\) instead.

Properties of the Verlet algorithm

  • The Verlet algorithm has no explicit velocities, but velocities can be computed.

  • It is easy to code, it’s fast and doesn’t require much memory.

  • The symmetric expansion around \(\vec{r}(t)\) renders the algorithm time-reversible.

  • It is energy conserving. Note that several of the other conservation laws cannot be checked as the Verlet algorithm doesn’t have explicit velocities.

Leapfrog algorithm#

The leapfrog algorithm is based on the Verlet algorithm. It adds explicit velocities to the Verlet algorithm. The name leapfrog comes from the fact that in this algorithm, the positions and velocities are half time step off from each other, that is, they are leapfrogging each other. Leapfrog is included as an integration algorithm in most simulation packages. The algorithm is mathematically equivalent to the Verlet algorithm

(9)#\[\begin{eqnarray} \vec{v}_{i}\left(t+\frac{\Delta{t}}{2}\right) &=& \vec{v}_{i}\left(t-\frac{\Delta{t}}{2}\right)+ \vec{a}_{i}(t)\Delta{t} \nonumber\\ \vec{r}_{i}(t+\Delta{t}) &=& \vec{r}_{i}(t)+\vec{v}_{i}\left(t+\frac{\Delta{t}}{2}\right)\Delta{t} \nonumber\\ \end{eqnarray}\]

Velocities at the full time step are obtained analogously as what is done in the basic Verlet algorithm,

(10)#\[\begin{equation} \vec{v}_{i}=\frac{\vec{v}_{i}\left(t+\frac{\Delta{t}}{2}\right)+\vec{v}_{i}\left(t-\frac{\Delta{t}}{2}\right)}{2}. \end{equation}\]

Algorithm 2 (The leapfrog algorithm)

  1. Initialize: Need initial values (\(t=t_0\)) for \(\vec{r}_i\) and \(\vec{v}_i\).

    • Compute the forces (accelerations) using the interaction potentials and the initial values for the positions. That gives \( \vec{a}_i(t_0)\).

    • With the above we have \(\vec{r}_i(t=t_0)\), \(\vec{v}_i(t=t_0)\) and \(\vec{a}_i(t=t_0)\).

Properties of the leapfrog algorithm

  • Has explicit velocities, but the velocties and positions are not computed at the same step.

  • It is easy to code and it is fast.

  • The symmetric expansion around \(\vec{r}(t)\) renders the algorithm time-reversible.

  • It is energy conserving.

  • The error is of order \((\Delta t)^2\) for moment and \((\Delta t)^2\) for positions.

  • In the Verlet algorithm the velocity term is a step behind the position term, that is, it is not self-starting as it needs initialization for \(x(t- \Delta t)\).

Velocity-Verlet#

The velocity-Verlet algorithm is one of the most commonly used MD algorithms and it is available in essentially all software packages. It fulfills the requirements listed at the beginning of this section. As the name suggests, it is based on the Verlet algorithm and adds explicit velocities. It is also worth noticing that mathematically the velocity-Verlet algorithm is identical to the Verlet algorithm

Algorithm 3 (Velocity-Verlet)

  1. Initialize: Need initial values (\(t=t_0\)) for \(\vec{r}_i\) and \(\vec{v}_i\).

    • Compute the forces (accelerations) using the interaction potentials and the initial values for the positions. That gives \( \vec{a}_i(t_0)\).

    • With the above we have \(\vec{r}_i(t=t_0)\), \(\vec{v}_i(t=t_0)\) and \(\vec{a}_i(t=t_0)\).

  2. Loop over time. At each step:

    1. Compute half-step velocities: \(\vec{v}_i(t + \frac{\Delta t}{2}) = \vec{v}_i(t) + \frac{1}{2} \vec{a}_i(t) \Delta t\).

    2. Compute the (full step) positions: \( \vec{r}_i(t + \Delta t) = \vec{r}_i(t) + \vec{v}_i (t+\frac{\Delta t}{2}) \Delta t\).

    3. Compute the forces (accelerations) using the interaction potential. That gives \( \vec{a}_i(t + \Delta t)\).

    4. Update the velocities (full-step): \( \vec{v}_i(t + \Delta t) = \vec{v}_i (t+\frac{\Delta t}{2}) + \frac{1}{2} \vec{a}_i(t + \Delta t) \Delta t\).

Properties of the velocity-Verlet algorithm

  • It is time reversible and symplectic.

  • Very stable.

  • It is easy to code. At the end of the integration step, all the positions and velocities are at the same time.

  • Possibly the most used MD simulation algorithm.