3 The MD Algorithm Solving Newton's equations of motion does not immediately suggest activity at the cutting edge of research.The molecular dynamics algorithm in most common use today may even have been known to Newton35.Nonetheless,the last decade has seen a rapid develop- ment in our understanding of numerical algorithms;a forthcoming review36 and a book37 summarize the present state of the field. Continuing to discuss,for simplicity,a system composed of atoms with coordinates rN(r2...N)and potential energy u().we introduce the atomic momenta pN-(p P2,...PN),in terms of which the kinetic energy may be written K(pN)- mThen the energy.or hamiltonian.may be written as a sum of kinetic and potential terms H=K+u.Write the classical equations of motion as ri=pi/mi and pi=fi (7) This is a system of coupled ordinary differential equations.Many methods exist to perform step-by-step numerical integration of them.Characteristics of these equations are:(a)they are 'stiff',i.e.there may be short and long timescales,and the algorithm must cope with both;(b)calculating the forces is expensive,typically involving a sum over pairs of atoms, and should be performed as infrequently as possible.Also we must bear in mind that the advancement of the coordinates fulfils two functions:(i)accurate calculation of dynamical properties,especially over times as long as typical correlation times Ta of properties a of interest (we shall define this later);(ii)accurately staying on the constant-energy hypersur- face,for much longer times Trun>Ta,in order to sample the correct ensemble. To ensure rapid sampling of phase space,we wish to make the timestep as large as possible consistent with these requirements.For these reasons,simulation algorithms have tended to be of low order (i.e.they do not involve storing high derivatives of positions, velocities etc.):this allows the time step to be increased as much as possible without jeop- ardizing energy conservation.It is unrealistic to expect the numerical method to accurately follow the true trajectory for very long times Tn.The 'ergodic'and'mixing'proper- ties of classical trajectories,i.e.the fact that nearby trajectories diverge from each other exponentially quickly,make this impossible to achieve. All these observations tend to favour the Verlet algorithm in one form or another,and we look closely at this in the following section.For historical reasons only,we mention the more general class of predictor-corrector methods which have been optimized for classical mechanical equations38.39;further details are available elsewhere34.40.41. 3.1 The Verlet Algorithm There are various,essentially equivalent,versions of the Verlet algorithm,including the original method7.42 and a'leapfrog'form43.Here we concentrate on the 'velocity Verlet' algorithm44,which may be written Pi(t+t)=pi(t)+tfi(t) (8a) ri(t+8t)=ri(t)+8tpi(t+t)/mi (8b) pP:(t+6t)=p:(t+6t)+6tf,(t+t) (8c) 73 The MD Algorithm Solving Newton’s equations of motion does not immediately suggest activity at the cutting edge of research. The molecular dynamics algorithm in most common use today may even have been known to Newton35 . Nonetheless, the last decade has seen a rapid develop￾ment in our understanding of numerical algorithms; a forthcoming review 36 and a book37 summarize the present state of the field. Continuing to discuss, for simplicity, a system composed of atoms with coordinates r N = (r1, r2, . . . rN ) and potential energy U(r N ), we introduce the atomic momenta p N = (p1 , p2 , . . . pN ), in terms of which the kinetic energy may be written K(p N ) = PN i=1  pi   2 /2mi . Then the energy, or hamiltonian, may be written as a sum of kinetic and potential terms H = K + U. Write the classical equations of motion as r˙ i = pi/mi and p˙ i = fi (7) This is a system of coupled ordinary differential equations. Many methods exist to perform step-by-step numerical integration of them. Characteristics of these equations are: (a) they are ‘stiff’, i.e. there may be short and long timescales, and the algorithm must cope with both; (b) calculating the forces is expensive, typically involving a sum over pairs of atoms, and should be performed as infrequently as possible. Also we must bear in mind that the advancement of the coordinates fulfils two functions: (i) accurate calculation of dynamical properties, especially over times as long as typical correlation times τa of properties a of interest (we shall define this later); (ii) accurately staying on the constant-energy hypersur￾face, for much longer times τrun  τa, in order to sample the correct ensemble. To ensure rapid sampling of phase space, we wish to make the timestep as large as possible consistent with these requirements. For these reasons, simulation algorithms have tended to be of low order (i.e. they do not involve storing high derivatives of positions, velocities etc.): this allows the time step to be increased as much as possible without jeop￾ardizing energy conservation. It is unrealistic to expect the numerical method to accurately follow the true trajectory for very long times τrun. The ‘ergodic’ and ‘mixing’ proper￾ties of classical trajectories, i.e. the fact that nearby trajectories diverge from each other exponentially quickly, make this impossible to achieve. All these observations tend to favour the Verlet algorithm in one form or another, and we look closely at this in the following section. For historical reasons only, we mention the more general class of predictor-corrector methods which have been optimized for classical mechanical equations38, 39; further details are available elsewhere34, 40, 41 . 3.1 The Verlet Algorithm There are various, essentially equivalent, versions of the Verlet algorithm, including the original method7, 42 and a ‘leapfrog’ form43 . Here we concentrate on the ‘velocity Verlet’ algorithm44 , which may be written pi (t + 1 2 δt) = pi (t) + 1 2 δtfi (t) (8a) ri(t + δt) = ri(t) + δtpi (t + 1 2 δt)/mi (8b) pi (t + δt) = pi (t + 1 2 δt) + 1 2 δtfi (t + δt) (8c) 7