Subsections

• Verlet Algorithm
• Velocity Verlet MD algorithm

Numerical Integration. Finite Difference methods

If we have a system of N particles, then equations (1) represent a system of 3N second-order differential equations. It can be reduced to the system of 6N first-order differential equations by introducing additional variable for the first derivative of coordinate. This leads to (2).

Both of these systems can be solved numerically using a finite difference technique. The general idea of the finite difference method is the following: given the positions, velocities, and other possible dynamic information at the time t, we attempt to obtain the position, velocities, etc. at a later time $t+\delta t$ with a sufficient degree of accuracy. Trajectory ($x(t), \dot{x}(t),$ etc.) is obtained by the step-by-step evaluation of these quantities. The choice of time step $\delta t$ depends upon the numerical integration method used and on the required accuracy. There exist many numerical methods to solve the ordinary differential equations. Unfortunately, not all of them satisfy the specific requirements of the molecular dynamics simulations:

• Large number of degrees of freedom requires the method to be easy on memory.
• The method should permit a long time step $\delta t$
• It should satisfy the known conservation laws for energy and momentum and be time-reversible.
• Expensive force calculations practically exclude high order and predictor-corrector methods.
• The methods should be simple in form and easy to program.

Verlet Algorithm

This algorithm, which exists in various versions, is the most popular due to its simplicity and efficiency.

Let's expand the coordinates at times $t+\delta t$ and at $t - \delta t$:

$$\begin{aligned} \mathbf{r}_i(t+\delta t) =&\mathbf{r}_i(t) +... ...ta t^3 \,\dddot{\mathbf{r}}_i(t) + O(\delta t^4) \end{aligned}$$ (3)

By subtracting and adding them, we obtain the main equations of the Verlet algorithm:  

$$\begin{aligned} \mathbf{r}_i(t+\delta t) = & 2 \mathbf{r}_i(... ...f{r}_i(t-\delta t) }{2\,\delta t} + O(\delta t^2) \end{aligned}$$ (4)

Its main strength is that it allows to calculate coordinates with the accuracy of $O(\delta t^4)$, storing only three arrays: two coordinates and accelerations, and calculating forces only once per time step. This scheme is symmetric with respect to time ($\mathbf{r}_i(t+\delta t)$ and $\mathbf{r}_i(t-\delta t)$ enter the equations (4) symmetrically), and hence the algorithm is time reversible. Velocities are not present explicitly in this version. Although they can be calculated, but with a second order accuracy. This is inconvenient, as this introduce some additional ``error'' in the energy or some other properties that depend on it. Another potential drawback is that because the small, of the order of $O(\delta t^2)$, number is added to the difference in the $O(\delta t^0)$ values, the rounding errors may be needlessly introduced. Better estimate for the velocity can be obtained by the following modification, also known as the ``velocity Verlet'' algorithm, which stores the positions, velocities and accelerations at the same time t:  

$$\mathbf{r}_i(t+\delta t) = \mathbf{r}_i(t) + \delta t \,\dot{\mathbf{r}}_i(t) + (1/2) \delta t^2\, \ddot{\mathbf{r}}_i(t)$$ (5)

 

$$\dot{\mathbf{r}_i}(t+\delta t) = \dot{\mathbf{r}_i}(t) + (1... ...dot{\mathbf{r}}_i(t) + \ddot{\mathbf{r}}_i(t+\delta t) \right ]$$ (6)

Instead of storing new and old accelerations, the velocity advancement is split in two steps to save memory:

1.

Before the forces (accelerations) are updated  

$$\dot{\mathbf{r}_i}(t+(1/2)\delta t) = \dot{\mathbf{r}_i}(t) + (1/2)\delta t \, \ddot{\mathbf{r}}_i(t)$$ (7)

and

2.

after the coordinates are updated and new forces (accelerations) $\ddot{\mathbf{r}}_i(t+\delta t)$are available:  

$$\dot{\mathbf{r}_i}(t+\delta t) = \dot{\mathbf{r}_i}(t+(1/2)\delta t) + (1/2)\delta t \, \ddot{\mathbf{r}}_i(t+\delta t)$$ (8)

Velocity Verlet MD algorithm

1.

prepare the initial configuration $\mathbf{r}_i(t)$ and velocities $\dot{\mathbf{r}}_i(t)$

2.

calculate energy and forces (accelerations) $\ddot{\mathbf{r}}_i(t)$, for this configuration

3.

calculate new positions $\mathbf{r}_i(t+\delta t)$ according to (5)

4.

advance velocities by equation (7).

5.

calculate potential energy and forces (accelerations) at the next step $\ddot{\mathbf{r}}_i(t+\delta t)$

6.

complete velocities update (equation (8))

7.

calculate new kinetic and total energies