calculation, methods that use only one force evaluation per time step

are to be preferred. This rules out the well-known Runge“Kutta

methods, which moreover are also not symplectic and lead to erro-

neous behavior on longer time scales (Leimkuhler, 1999).

6.4 Solving the equations of motion 191

In the past practice of MD the Gear algorithm has been much used. The

Gear algorithm predicts positions and a number of derivatives based on a

Taylor expansion of previous values (how many depends on the order of

the algorithm); it then evaluates the accelerations from the forces at the

predicted position, and corrects the positions and derivatives on the basis

of the deviation between predicted and evaluated accelerations. There are

several variants and predictor“corrector algorithms of this kind have been

applied up to orders as high as eight. They are quite accurate for small time

steps but not very stable for larger time steps. When the forces are not very

precise, it does not help to use high orders. They are neither time-reversible

nor symplectic and have not survived the competition of the simpler, more

robust, reversible and symplectic Verlet or leap-frog algorithms.36 The latter

and its variants including multiple time-step versions, are derivable by the

Reference System Propagator Algorithms (RESPA) method of Tuckerman

et al. (1992) that uses simple operator algebra and is outlined below.

The original Verlet algorithm (Verlet, 1967) does not use the velocities,

and employs a simple discretization of the second derivative:

x(t ’ ”t) ’ 2x(t) + x(t + ”t)

x(t) ≈

¨ , (6.81)

(”t)2

leading to the predicted position (x stands for every particle coordinate;

f (t) = Fi (x(t))/mi is the corresponding force component, evaluated from

the positions at time t and “ for convenience “ divided by the mass)

x(t + ”t) = 2x(t) ’ x(t ’ ”t) + f (t)(”t)2 + O((”t)4 ). (6.82)

The velocity is found in retrospect from

v(t + ”t) ’ v(t ’ ”t)

+ O((”t)2 ),

v(t) = (6.83)

2”t

but plays no role in the evolution of the trajectory. It can be more accurately

estimated from37

v(t + ”t) ’ v(t ’ ”t) f (t ’ ”t) ’ f (t + ”t)

+ O((”t)3 ). (6.84)

v(t) = +

2”t 12

An equivalent scheme is the leap-frog algorithm, which uses positions at

integer time steps and velocities halfway in between time steps (Hockney

36 Gear algorithms (Gear, 1971) and their variants have been reviewed and evaluated for use in

MD including constraints by van Gunsteren and Berendsen (1977) and in relation to the Verlet

algorithm by Berendsen and van Gunsteren (1986).

37 Berendsen and van Gunsteren (1986).

192 Molecular dynamics

and Eastwood, 1988). Starting from v(t ’ 1 ”t) and x(t) the updates are

2

1 1

= v t ’ ”t + f (t)”t,

v t + ”t

2 2

1

x(t + ”t) = x(t) + v t + ”t ”t. (6.85)

2

It can be easily shown that this algorithm is equivalent to Verlet™s and

will generate the same trajectory, if the velocity v(t ’ 1 ”t) is started as

2

[x(t) ’ x(t ’ ”t)]/(”t). The velocity at integer time steps can be recovered

as the average of the two velocities at half time steps earlier and later, but

only to the O((”t)2 ) precision of (6.83).

In several applications, for example when velocity-dependent forces are

applied, it is desirable to know the velocity at the time the position is pre-

dicted, rather than a time step later. There are several algorithms, equiva-

lent to Verlet, that deliver equal-time velocities. One is Beeman™s algorithm

(Beeman, 1976):

2 1

f (t) ’ f (t ’ ”t) (”t)2 ,

x(t + ”t) = x(t) + v(t)”t +

3 6

1 5 1

v(t + ”t) = v(t) + f (t + ”t) + f (t) ’ f (t ’ ”t) ”t, (6.86)

3 6 6

but the most popular one is the velocity-Verlet algorithm (Swope et al.,

1982):

1

x(t + ”t) = x(t) + v(t)”t + f (t)(”t)2 ,

2

1

v(t + ”t) = v(t) + [f (t) + f (t + ”t)]”t. (6.87)

2

This algorithm needs the force at the new time step, but there is only one

force evaluation per step. Although it is not immediately obvious, all these

algorithms are equivalent (Berendsen and van Gunsteren, 1986).

The elegant operator technique considers the exponential Liouville oper-

ator to evolve the system in time. We start with (17.152) on page 493:

™

z = iLz, (6.88)

where z is the vector of generalized coordinates and conjugate momenta.

We apply (6.88) to the cartesian Newtonian equations (6.79) and introduce

the time-di¬erential operator D:

01

x

™ x x

= =D , (6.89)

ˆ

v

™ v v

f0

6.4 Solving the equations of motion 193

ˆ

where f is an operator acting on x that produces Fi (x)/mi , (i = 1, . . . 3N ).

The x, v vector has a length 6N for N particles. The solution is

x x x

(t) = etD (0) = U(t) (0). (6.90)

v v v

The exponential operator U = exp(tD) is time-reversible:

U(’t) = U’1 (t), (6.91)

and this solution is exact and symplectic.38 Unfortunately we cannot solve

(6.90) and we have to solve the evolution over small time steps by approxi-

mation. First split the operator D into two simple parts:

01 00

01

= + . (6.92)

ˆ ˆ

00

f0 f0

These two parts do not commute and we can use the Trotter“Suzuki expan-

sion (see Chapter 14, page 386)

et(A+B) ≈ e(t/2)A etB e(t/2)A , (6.93)

which can be further subdivided into higher products with higher-order ac-

curacy. Substituting the exponential operators A and B by

00

01

Uv (t) = t exp and Uf (t) = t exp , (6.94)

ˆ

00 f0

and using a ¬rst-order expansion of the exponential:

x x + vt

Uv (t) = , (6.95)

v v

x x

Uf (t) = , (6.96)

v v + ft

we can, for example, split U(”t) as follows: