Integrator: Predictor-Corrector Algorithm

      The main algorithm of MD simulations is the integrator that calculates the trajectories with time using the equation of motion for interacting particles. There are a few types of numerical analysis methods for integrating Newton's equations, such as Verlet, leapfrog, and the predictor-corrector algorithm, which is used in the MD simulations discussed in study. The positions, velocities, and accelerations of particles are related with each other according to Eq. 1.
..........(1)

The predictor-corrector algorithm is based on a Taylor expansion.[1] If the position r, velocity v, acceleration a and time derivative of the acceleration (b) are known at time t, these quantities after time step Δt can be predicted as shown in the following equation:
..........(2a)
..........(2b)
..........(2c)
..........(2d)

If the Taylor expansions are truncated, so that only the terms shown explicitly in Eq. 2 are left, then the quantities can be called the predicted values rp, vp, ap, and bp. The force is computed by the gradient of potential at the predicted position rp, because the predicted values are not based on physics. The re-calculated acceleration with Eq. 1 is different from the predicted acceleration ap, and the difference between the two values is called the error signal, as shown in Eq. 3.
..........(3)

This term for the error signal is used to correct all the predicted quantities as follows:
..........(4a)
..........(4b)
..........(4c)
..........(4d)
All the corrected quantities are proportional to the error signal, and the proportional coefficients are determined to maximize the stability of the calculation. These corrected values are now better approximations of the true quantities, and are used to predict the quantities in the next iteration. The best choice for these coefficients depends on the order of both the differential equations and the Taylor series.[2] These coefficients are computed based on the order of the algorithm being used in the simulation. The algorithm used here is a third-order Nordsieck predictor-corrector algorithm, and the values are c0=1/6, c1=5/6, c2=1, and c3=1/3. Inherently, all the integrator algorithms generate errors due to the truncation of the Taylor expansion in addition to the round-off errors caused by the finite number of digits in the computer calculations. Therefore, the order of truncated-differential equations is chosen to balance accuracy and computational speed. In addition, the accuracy of the numerical integrator algorithms also depends on the time step size, which is typically on the order of fractions of femto-seconds (10-15 sec). Thus, the simulation as a whole is able to describe only short-time scale phenomena that last on the order of pico- (10-12) or nano-seconds (10-9 sec). However, these limitations have been abridged by the development of sophisticated computer and numerical algorithms.

References
[1] Frenkel, D. and B. Smit, Understanding Molecular Simulation. 1996, San Diego: Academic Press. 28.
[2] Gear, C.W., Numerical Initial Value Problems in Ordinary Differential Equations. 1971, Englewood Cliffs, NJ: Prentice-Hall.