Mathematical Physics - Volume II - Numerical Methods

Chapter 6. Introduction to Computational Mechanics of Discontinua

228

φ LJ = φ spl , φ ′ LJ = φ ′ spl and φ ′′

′′ spl . Note that, since the cubic spline is defined in r 2

LJ = φ

instead of r , rooting and division operations with r are avoided in the MD simulation, while differences in the appearance of the function are often practically imperceptible [27]. Physical quantities are, of course, somewhat influenced by this shortening of the range of potentials. Overall, the eventual assessment of these effects is a problem specific to each application (e.g., [16], [17]). 6.2.3 Integration Algorithms The core of the MD program is the integration algorithm necessary for the integration of the equations of motion (6.1) 1 . These algorithms are based on the finite difference method in which time is discretized within a finite network and the time step, δ t , represents the distance between successive points on the network. Knowing the position and some of its time derivatives in the current time ( t ), allows calculation of the same quantities in the next time ( t + δ t ) . In addition to being as accurate as possible, the integration algorithm should be both fast and of modest memory requirements, to allow the application of the longest possible time step, and to be simple for implementation. Several different methods are available for this purpose but only one will be discussed herein. The Verlet algorithm is still the most widely used, although more powerful techniques for integrating finite difference equations are available. The derivation of this algorithm is based on the development of the atomic position function in the Taylor series at an arbitrary time t

d 3 r ( t ) d t 3

r i ( t + δ t ) = r i ( t )+ v i ( t ) δ t + a i ( t ) δ t 2 +

δ t 3 + O ( δ t 4 )

that is

d 3 r ( t ) d t 3

r i ( t − δ t ) = r i ( t ) − v i ( t ) δ t + a i ( t ) δ t 2 −

δ t 3 + O ( δ t 4 ) .

Adding the previous two expressions yields

2 + O ( δ t 4 ) .

r i ( t + δ t )+ r i ( t − δ t ) = 2 r i ( t )+ 2 a i ( t ) δ t

The time-reversible equation for calculating the next position has the form

f i j ( r i , r j ) m i

2 ∑ j̸ = i

+ O ( δ t 4 ) .

r i ( t + δ t ) = 2 r i ( t ) − r i ( t − δ t )+ δ t

(6.9)

This method of integration is very compact and easy to implement. Since no dissipative forces act between the atoms, the dynamic system is conservative. Therefore, the force by which the atom j acts on the atom i is f i j = − ∇ i φ ( r i j ) , and its calculation is by far the most demanding part of MD simulations. The Verlet method seems to be the least time consuming and thus the most suitable for modest computing resources, which explains its popularity. The main disadvantages of the Verlet algorithm are clumsy handling of velocities v i ( t ) = r i ( t + δ t ) − r i ( t − δ t ) 2 δ t + O ( δ t 2 )