Molecular Dynamics - Integration Algorithms

Introduction
Criteria of Good Integrators in Molecular Dynamics
Integrators in the Discover 2.9.5 and 95.0 Programs

Introduction

At its simplest, molecular dynamics solves Newton's familiar equation of motion:

Eq. 5-1:

where F_i is the force, m_i is the mass, and a_i is the acceleration of atom i.

The force on atom i can be computed directly from the derivative of the potential energy V with respect to the coordinates r_i:

Eq. 5-2:

Notice that classical equations of motion are deterministic. That is, once the initial coordinates and velocities are known, the coordinates and velocities at a later time can be determined. The coordinates and velocities for a complete dynamics run are called the trajectory.

A standard method of solving ordinary differential equation such as Eq. 5-2 numerically is the finite difference method. The general idea is as follows. Given the initial coordinates and velocities and other dynamic information at time t, the positions and velocities at time t + t are calculated. The timestep t depends on the integration method as well as the system itself. All the solution methods offered by the Discover program belong to this category.

Although the initial coordinates are determined in the input file or from a previous operation such as minimization, the initial velocities are randomly generated at the beginning of a dynamics run, according to the temperature and a random number seed. To repeat the exact same dynamics run, it is important to specify the random number seed used in the previous run, as found in the output file. This guarantees generating the same set of initial velocities in the current run. More details on the initial velocities are provided in the temperature section.

Criteria of Good Integrators in Molecular Dynamics

Molecular dynamics is usually applied to a large system. Energy evaluation is time consuming and the memory requirement is large. To generate the correct statistical ensembles, energy conservation is also important. Thus, the basic criteria for a good integrator for molecular simulations are as follows:

It should be fast, ideally requiring only one energy evaluation per timestep.
It should require little computer memory.
It should permit the use of a relatively long timestep.
It must show good conservation of energy.

Integrators in the Discover 2.9.5 and 95.0 Programs

Integrators provided in the Discover program are chosen according to those criteria. The Discover 2.9.5 program uses the Verlet leapfrog algorithm. In the Discover 95.0 program, you can choose one of several integrators: Verlet velocity, ABM4, and Runge-Kutta-4. The Verlet velocity method is the default, since it satisfies the above criteria best.

Verlet Leapfrog Integrator

Variants of the Verlet (1967) algorithm of integrating the equations of motion (Eq. 5-2) are perhaps the most widely used method in molecular dynamics. The Discover program uses the leapfrog version in release 2.9.5 and the velocity version for release 95.0. The advantages of Verlet algorithms is that it requires only one energy evaluation per step, requires only modest memory, and also allows a relatively large timestep to be used.

The Verlet leapfrog algorithm is as follows:

Given r(t), v(t -t/2), and a(t), which are (respectively) the position, velocity, and acceleration at times t, t -t/2, and t, compute:

where f(t + t) is evaluated from -dV/dr at r(t + t).
The Verlet leapfrog method has one major disadvantage: the positions and velocities calculated are half a timestep out of synchrony.

Verlet Velocity Integrator

To overcome the shortcoming of the Verlet leapfrog method, the Verlet velocity algorithm is implemented in release 95.0 of the Discover program.

The algorithm is as follows:
Given r(t), v(t), and a(t), which are (respectively) the position, velocity, and acceleration at time t, compute:

ABM4 Integrator

ABM4, which stands for Adams-Bashforth-Moulton fourth order, is a predictor and corrector method. It is a fourth-order method, meaning that the truncation error is to the fifth order of the timestep used. This method requires two energy evaluations per step and has to make use of the results of the previous three steps. It is thus not self starting--the first three steps are generated by the Runge-Kutta method. More memory has to be used, because previous information has to be stored.

The algorithm is as follows:

Let:

where subscripts (not shown above) 0, 1, -1, -2, and -3 indicate y at times t, t + t, t - t, t - 2t, and t - 3t.
The predictor is:

Now evaluate y'₁, using the predicted y₁, which involves one energy evaluation.
The corrector is:

Now evaluate y'₁, using the corrected y₁, which involves another energy evaluation.

Runge-Kutta-4 Integrator

Runge-Kutta-4 stands for the fourth-order Runge-Kutta method, which is one of the oldest numerical methods for solving ordinary differential equations. The method is self starting but requires four energy evaluations per step. From testing done at Biosym/MSI, the timestep has to be very small. This is thus not a very suitable integrator for molecular simulation. However, the method is very robust, meaning that it can deal with almost all kinds of equations, including stiff ones. This integrator is used to generate the trajectory for the first three steps for ABM4. Since this integrator should not be frequently used, the algorithm is not presented here. Details can be found in Press et al. 1986.

Main access page

Theory/Methodology access.

Dynamics access Dynamics - Background Dynamics - Choice of Timestep