9 LECTURE Numerical Methodology
The preceding lecture focused on the description of a system comprising N particles, the intermolecular forces, the governing equations describing the dynamics of N particles, and the boundary conditions. In this lecture we will discuss some of the numerical methids used to solve the dynamics equations. The computation of the trajectories of N particles involves the solution of a system of 6N coupled, first-order ordinary differential equations. A variety of numerical methods are available to solve these equations. We will discuss some of these methods and other related issues relevant to MD simulations. These methods may be broadly classified into four groups:
1. 2. 3. 4. Runge-Kutta (R-K) Methods Leap Frog Methods Multi-Step and Predictor Corrector Methods Gear Predictor-Corrector Method
9.1 Runge-Kutta (R-K) Methods
R-K methods of various orders have been commonly used to solve ordinary differential equations (ODEs) involving initial value problems. In order to understand the fundamental concepts, let us start with the basic Euler method used to solve a first-order ODE:
over an interval x=a to x=b. As discussed in Lecture 8, the particle dynamics in a one-dimensional formulation is governed by two ODEs of the type given in Eq. (9.1).The integration of Eq. (9.1) yields
Here h=t represents the temporal step size and x the particle position. Various R-K methods differ essentially in the approximation used to evaluate the above integral. For the Euler method, f(x,t) is assumed to be constant and equal to its value f (xi,t) over the interval dt. This is equivalent to writing the Taylor series for x (t+h) as
The method is first-order accurate with the local truncation error being proportional to h2, or O(h2). It can easily be shown that while the local truncation error of O (h2), the global truncation error, which is accumulated over the entire computational interval (b-a) is O (h). The Euler method can also be considered as the first-order R-K method, since various R-K methods differ in how the slope is evaluated over the interval t and t + h. Each R-K method evaluates the slope at different points in the interval (h), and computes a weighted average that is used in Eq. (9.2). For example, a second-order R-K method yields
This method has a local truncation error (lte) of O (h3), and requires two computations of the function f(x,t) for each time step. Similarly a fourthorder R-K method, which is commonly used for solving a system of ODEs, has a lte of O (h5) and requires four computations of the function f(x,t) for each time step. Note that the particle velocity is calculated using the same approach, with f(x,t) now representing the force (or acceleration in terms of the non-dimensional variables). Further details can be found in Ref. 7.
R-K methods are self-starting, easy to program, and the step size can be automatically controlled to achieve the desired accuracy. However, these methods are less desirable for MD simulations due to two considerations. One is that higher-order R-K methods do not have any significant advantage over the low-order R-K methods. This is because the accuracy of trajectory calculations, due to the strong repulsive forces at short distances, puts an upper limit on the temporal step size (h). The second consideration is the number of evaluations of particle force or acceleration per time step. For example, a fourth-order R-K method involves four evaluations per time step. Consequently, when MD simulations involve a long time period, the R-K methods would require excessive amount of computational time. Then other methods such as predictor-corrector methods are preferred.
9.1.1 Various Types of Errors
Here we will discuss two major types of errors; truncation error and round-off error (roe). As indicated above, the Euler method, Eq. (9.3), has a local truncation error given by
A common terminology to quantify this error is to state that the local truncation error is of the order of h2 or O(h2). Also, it can easily be shown that the total truncation error accumulated over a finite number of steps is O(h). The corresponding local and total truncation errors for the second-order R-K method are O(h3) and O(h2). Thus the truncation error is determined by the time step (h= t) and the accuracy of the numerical algorithm used. The roe on the other hand depends on the number of significant figures used in the computation, the order by which the computations are performed, and the approximations used for calculating square roots, exponentials, logarithms etc. Since roe depends on the total number of computations or time steps in the simulation, it also depends on t. The global error is the sum of truncation error and round-off error, and, thus, clearly depends on t. Generally speaking as t is decreased, the te decreases, while the roe increases. Thus the choice of t should be such that it corresponds to the minimum total error, as shown in Fig. 9.1, taken from Haile [1].
Fig. 9.1 Global error ploted versus the non-dimensional time step. For this figure, the average global error is computed in terms of the nondimensional energy of the system and plotted versus the nondimensional time step. The average global error is defined as
where E* is the non-dimensional energy defined as E*=E/e, is the non-dimensional time, and M the number of time steps. The total error in energy is computed for a one-dimensional
collision between two Lennard-Jones atoms using the Gear's predictorcorrector method, which is discussed later in this lecture. Results are shown for simulations performed using both the single- and doubleprecision. As indicated in the figure, for t*5x10-3, the te makes the largest contribution to the global error, and thus, the single- and doubleprecision computations yield the same global error, i.e. the roe is of no consequence. However, for t*<10-3 global errors with single precision are dominated by round-off, but in double-precision roe remain unimportant for t* as small as 10-5. There are many ways to monitor the total error in the simulation. As indicated in the above example, a common procedure in MD simulations is to calculate the total error by monitoring the total energy as a function of time. For additional discussion on various types of errors, the reader is referred to Refs. [1, 7]. Note all the references pertaining to the material on MD simulations are listed in numerical order in Lecture 7.
9.2 Leap Frog Method
Writing the Taylor series for x (t+h) as
Here represent the particle velocity and acceleration respectively. The above expression can also be written as
where
is expressed as
The new particle position is calculated using Eqs. (9.6) and (9.7), while the particle velocity is calculated from
Here or the particle acceleration is obtained from the net force on the particle using the approach described in Lecture 8. Note that the particle position and velocity are being calculated at different times, which does not represent any problem. In describing initial conditions, is also ignored. Another the distinction between method which has been used more commonly compared to the leap frog method is the Verlet algorithm.
9.2.1 Verlet Algorithm
Writing the Taylor series for x (t-h) and combining it with Eq. (9.5) the yields Verlet formula
which is third-order accurate, and requires only one force evaluation per time step. The velocity is then obtained using
Higher-order algorithms are available to compute velocity, but generally are not used. The Verlet algorithm is a two-step method, and not selfstarting since it requires x(t-h). A special procedure is needed to obtain x at t=-h. The main advantages of this method are its efficiency, simplicity, relatively higher-order accuracy, and good stability characteristics.
9.3 Multi-Step and Predictor-Corrector Methods
These methods generally involve two steps for integration over one time step (h=t); the predictor step in which the particle position and velocity are predicted, followed by the corrector step in which the predicted values are corrected in order to improve the accuracy and/or the stability of the algorithm. Either of these steps can employ a self-starting method, such as Runge-Kutta methods, or a multi-step method, although the latter type has been used more commonly. In this section, we will first discuss the multi-step methods, followed by a discussion on the predictor-corrector algorithm. We will end the
section by discussing a specific predictor-corrector methodology that is commonly used in MD simulations. The starting point for deriving formulas for multi-step methods is the Taylor series for x (t+h)
we obtain a formula for a second-order multi-step method
Since this formula requires information at two time steps, t and t-h, the algorithm is not self-starting. Similarly, we can derive higher-order formulas by considering higher-order terms in the Taylor series, and using finite-difference approximations for these terms. For example, considering the third-order term in Eq. (9.11), and using approximation for and x(3)(t) as
yields a third-order multi-step algorithm, i.e.,
Note the superscript (3) with x is used to represent the third derivative of x. Again the algorithm is not self-starting, as it requires information
at three time steps. Multi-step algorithms with fourth- and higherorder accuracy can be obtained following a similar procedure. Similar formulas can be used for computing the particle velocity at the new step. The major advantage of multi-step methods is their efficiency since only one function evaluation is needed at each time step, although the storage requirement is increased. The above formulas provide explicit expressions for the unknown variable, and, therefore, are known as explicit methods. Another class of multi-step methods involves implicit expressions for the unknown. To derive such expressions, we write the Taylor series for x(t) as
and then rewrite it in terms of x(t+h)
Using the following finite-difference approximations for x(3)(t+h)
and
we obtain a third-order implicit formula for x(t+h
) Similarly we can obtain higher-order implicit formulas.
9.3.1 Predictor-Corrector Methods
The predictor-corrector methods consist of two steps; the predictor step and the correction step. The predictor step generally employs an explicit formula, while the correction step employs an implicit formula. For example, we may use the third-order explicit formula, Eq. (9.14), in the predictor step to calculate an intermediate value of particle position, i.e.,
A similar formula is used to calculate the intermediate value of particle velocity. These values are then corrected by using a thirdorder implicit formula, i.e., Eq. (9.17) as
represents the intermediate particle velocity where calculated in the predictor step. Again a similar formula is used to calculate the particle velocity in the corrector step. The particle acceleration, , required in the corrector step, is calculated using the particle position, , which has been obtained in the predictor step. Since many different explicit and implicit multi-step formulas can be used, a variety of predictor-corrector algorithms are available in the literature [7]. Here we will discuss a fifth-order predictor-corrector algorithm due to Gear [8], which has been commonly used in MD simulations.
9.4 Gear Predictor-Corrector Method
The algorithm consists of the following three steps. Predictor Step Taylor series is used to predict x(t+h) and higher derivatives
The higher order derivatives needed in the above equations are taken from the previous simulation. For starting the simulation, the...