《统计热力学》课程教学资源（参考文献）Allen. M.P.Introduction to MD

John von Neumann Institute for Computing NIC Introduction to Molecular Dynamics Simulation Michael P.Allen published in Computational Soft Matter:From Synthetic Polymers to Proteins, Lecture Notes, Norbert Attig,Kurt Binder,Helmut Grubmuller,Kurt Kremer(Eds.), John von Neumann Institute for Computing,Julich, N1 C Series,VolL.23,1SBN3-00-012641-4,pp.1-28,2004. C 2004 by John von Neumann Institute for Computing Permission to make digital or hard copies of portions of this work for personal or classroom use is granted provided that the copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page.To copy otherwise requires prior specific permission by the publisher mentioned above. http://www.fz-juelich.de/nic-series/volume23

John von Neumann Institute for Computing Introduction to Molecular Dynamics Simulation Michael P. Allen published in Computational Soft Matter: From Synthetic Polymers to Proteins, Lecture Notes, Norbert Attig, Kurt Binder, Helmut Grubmuller ¨ , Kurt Kremer (Eds.), John von Neumann Institute for Computing, Julich, ¨ NIC Series, Vol. 23, ISBN 3-00-012641-4, pp. 1-28, 2004. c 2004 by John von Neumann Institute for Computing Permission to make digital or hard copies of portions of this work for personal or classroom use is granted provided that the copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. To copy otherwise requires prior specific permission by the publisher mentioned above. http://www.fz-juelich.de/nic-series/volume23

Introduction to Molecular Dynamics Simulation Michael P.Allen Centre for Scientific Computing and Department of Physics, University of Warwick,Coventry CV4 7AL,United Kingdom E-mail:m.p.allen@warwick.ac.uk In this chapter a summary is given of the key ingredients necessary to carry out a molecular dynamics simulation,with particular emphasis on macromolecular systems.We discuss the form of the intermolecular potential for molecules composed of atoms,and of non-spherical sub-units,giving examples of how to compute the forces and torques.We also describe some of the MD algorithms in current use.Finally.we briefly refer to the factors that influence the size of systems,and length of runs,that are needed to calculate statistical properties. 1 The Aims of Molecular Dynamics We carry out computer simulations in the hope of understanding the properties of assem- blies of molecules in terms of their structure and the microscopic interactions between them.This serves as a complement to conventional experiments,enabling us to learn something new,something that cannot be found out in other ways.The two main families of simulation technique are molecular dynamics(MD)and Monte Carlo(MC);addition- ally,there is a whole range of hybrid techniques which combine features from both.In this lecture we shall concentrate on MD.The obvious advantage of MD over MC is that it gives a route to dynamical properties of the system:transport coefficients,time-dependent responses to perturbations,rheological properties and spectra. Computer simulations act as a bridge(see Fig.1)between microscopic length and time scales and the macroscopic world of the laboratory:we provide a guess at the interactions between molecules,and obtain 'exact'predictions of bulk properties.The predictions are 'exact'in the sense that they can be made as accurate as we like,subject to the limita- tions imposed by our computer budget.At the same time,the hidden detail behind bulk measurements can be revealed.An example is the link between the diffusion coefficient and velocity autocorrelation function(the former easy to measure experimentally,the latter much harder).Simulations act as a bridge in another sense:between theory and experi- ment.We may test a theory by conducting a simulation using the same model.We may test the model by comparing with experimental results.We may also carry out simulations on the computer that are difficult or impossible in the laboratory (for example,working at extremes of temperature or pressure). Ultimately we may want to make direct comparisons with experimental measurements made on specific materials,in which case a good model of molecular interactions is essen- tial.The aim of so-called ab initio molecular dynamics is to reduce the amount of fitting and guesswork in this process to a minimum.On the other hand,we may be interested in phenomena of a rather generic nature,or we may simply want to discriminate between good and bad theories.When it comes to aims of this kind,it is not necessary to have a perfectly realistic molecular model;one that contains the essential physics may be quite suitable

Introduction to Molecular Dynamics Simulation Michael P. Allen Centre for Scientific Computing and Department of Physics, University of Warwick, Coventry CV4 7AL, United Kingdom E-mail: m.p.allen@warwick.ac.uk In this chapter a summary is given of the key ingredients necessary to carry out a molecular dynamics simulation, with particular emphasis on macromolecular systems. We discuss the form of the intermolecular potential for molecules composed of atoms, and of non-spherical sub-units, giving examples of how to compute the forces and torques. We also describe some of the MD algorithms in current use. Finally, we briefly refer to the factors that influence the size of systems, and length of runs, that are needed to calculate statistical properties. 1 The Aims of Molecular Dynamics We carry out computer simulations in the hope of understanding the properties of assem￾blies of molecules in terms of their structure and the microscopic interactions between them. This serves as a complement to conventional experiments, enabling us to learn something new, something that cannot be found out in other ways. The two main families of simulation technique are molecular dynamics (MD) and Monte Carlo (MC); addition￾ally, there is a whole range of hybrid techniques which combine features from both. In this lecture we shall concentrate on MD. The obvious advantage of MD over MC is that it gives a route to dynamical properties of the system: transport coefficients, time-dependent responses to perturbations, rheological properties and spectra. Computer simulations act as a bridge (see Fig. 1) between microscopic length and time scales and the macroscopic world of the laboratory: we provide a guess at the interactions between molecules, and obtain ‘exact’ predictions of bulk properties. The predictions are ‘exact’ in the sense that they can be made as accurate as we like, subject to the limita￾tions imposed by our computer budget. At the same time, the hidden detail behind bulk measurements can be revealed. An example is the link between the diffusion coefficient and velocity autocorrelation function (the former easy to measure experimentally, the latter much harder). Simulations act as a bridge in another sense: between theory and experi￾ment. We may test a theory by conducting a simulation using the same model. We may test the model by comparing with experimental results. We may also carry out simulations on the computer that are difficult or impossible in the laboratory (for example, working at extremes of temperature or pressure). Ultimately we may want to make direct comparisons with experimental measurements made on specific materials, in which case a good model of molecular interactions is essen￾tial. The aim of so-called ab initio molecular dynamics is to reduce the amount of fitting and guesswork in this process to a minimum. On the other hand, we may be interested in phenomena of a rather generic nature, or we may simply want to discriminate between good and bad theories. When it comes to aims of this kind, it is not necessary to have a perfectly realistic molecular model; one that contains the essential physics may be quite suitable. 1

Experimental Results Intermolecular Complex Fluid potential Test model (real system) Structure v(r g(r Simulation Make model Results Complex Fluid Test theory (model system) c(t) Dynamics Theoretical Predictions Figure 1.Simulations as a bridge between (a)microscopic and macroscopic;(b)theory and experiment. 2 Molecular Interactions Molecular dynamics simulation consists of the numerical,step-by-step,solution of the classical equations of motion,which for a simple atomic system may be written miri=fi fi=-Ori (1) For this purpose we need to be able to calculate the forces f;acting on the atoms,and these are usually derived from a potential energy u(rN),where rN =(r1,r2,...rN)repre- sents the complete set of 3N atomic coordinates.In this section we focus on this function (rN),restricting ourselves to an atomic description for simplicity.(In simulating soft condensed matter systems,we sometimes wish to consider non-spherical rigid units which have rotational degrees of freedom:rotational equations of motion and interaction poten- tials will be considered in section 5). 2.1 Non-bonded Interactions The part of the potential energy Uhon-bonded representing non-bonded interactions between atoms is traditionally split into 1-body,2-body,3-body...terms: a-bondod(r）=∑a(r)+∑(r,r)+ (2) i j>i The u(r)term represents an externally applied potential field or the effects of the container walls;it is usually dropped for fully periodic simulations of bulk systems.Also,it is usual to concentrate on the pair potential v(ri,rj)=v(r)and neglect three-body (and higher order)interactions.There is an extensive literature on the way these potentials are determined experimentally,or modelled theoretically1-4. In some simulations of complex fluids,it is sufficient to use the simplest models that faithfully represent the essential physics.In this chapter we shall concentrate on continu- ous,differentiable pair-potentials(although discontinuous potentials such as hard spheres 2

Intermolecular T r t c(t) P Structure Dynamics Phases s l g g(r) r v(r) potential Complex Fluid (real system) Complex Fluid (model system) Make model Simulation Experimental Results Results Predictions Theoretical Test theory Test model Figure 1. Simulations as a bridge between (a) microscopic and macroscopic; (b) theory and experiment. 2 Molecular Interactions Molecular dynamics simulation consists of the numerical, step-by-step, solution of the classical equations of motion, which for a simple atomic system may be written mir¨i = fi fi = − ∂ ∂ri U (1) For this purpose we need to be able to calculate the forces f i acting on the atoms, and these are usually derived from a potential energy U(r N ), where r N = (r1, r2, . . . rN ) repre￾sents the complete set of 3N atomic coordinates. In this section we focus on this function U(r N ), restricting ourselves to an atomic description for simplicity. (In simulating soft condensed matter systems, we sometimes wish to consider non-spherical rigid units which have rotational degrees of freedom: rotational equations of motion and interaction poten￾tials will be considered in section 5). 2.1 Non-bonded Interactions The part of the potential energy Unon-bonded representing non-bonded interactions between atoms is traditionally split into 1-body, 2-body, 3-body . . . terms: Unon-bonded(r N ) = X i u(ri) + X i X j>i v(ri , rj ) + . . . . (2) The u(r) term represents an externally applied potential field or the effects of the container walls; it is usually dropped for fully periodic simulations of bulk systems. Also, it is usual to concentrate on the pair potential v(ri , rj ) = v(rij ) and neglect three-body (and higher order) interactions. There is an extensive literature on the way these potentials are determined experimentally, or modelled theoretically1–4 . In some simulations of complex fluids, it is sufficient to use the simplest models that faithfully represent the essential physics. In this chapter we shall concentrate on continu￾ous, differentiable pair-potentials (although discontinuous potentials such as hard spheres 2

Lennard-Jones 4斯6 WCA 0 -1 0 0.5 1.5 2 2.5 3 r/6 Figure 2.Lennard-Jones pair potential showing.Also shown is the WCA shifted repulsive part of the potential. and spheroids have also played a role,see e.g.5).The Lennard-Jones potential is the most commonly used form: )=[)）“-() (3) with two parameters:o,the diameter,and s,the well depth.This potential was used, for instance,in the earliest studies of the properties of liquid argon.7 and is illustrated in Fig.2.For applications in which attractive interactions are of less concern than the excluded volume effects which dictate molecular packing,the potential may be truncated at the position of its minimum,and shifted upwards to give what is usually termed the WCA models.If electrostatic charges are present,we add the appropriate Coulomb potentials Coulomb(r）= QiQ2 (4) 4Teor where Q1,Q2 are the charges and eo is the permittivity of free space.The correct handling of long-range forces in a simulation is an essential aspect of polyelectrolyte simulations, which will be the subject of the later chapter of Holm. 2.2 Bonding Potentials For molecular systems,we simply build the molecules out of site-site potentials of the form of Eg.(3)or similar.Typically,a single-molecule quantum-chemical calculation may be used to estimate the electron density throughout the molecule,which may then be mod- elled by a distribution of partial charges via Eq.(4),or more accurately by a distribution of electrostatic multipoles4.10.For molecules we must also consider the intramolecular bond- ing interactions.The simplest molecular model will include terms of the following kind:

0 0.5 1 1.5 2 2.5 3 r / σ -2 -1 0 1 2 3 4 v / ε Lennard-Jones -4r-6 4r-12 WCA Figure 2. Lennard-Jones pair potential showing the r−12 and r−6 contributions. Also shown is the WCA shifted repulsive part of the potential. and spheroids have also played a role, see e.g. 5). The Lennard-Jones potential is the most commonly used form: v LJ(r) = 4ε σ r 12 − σ r 6  . (3) with two parameters: σ, the diameter, and ε, the well depth. This potential was used, for instance, in the earliest studies of the properties of liquid argon6, 7 and is illustrated in Fig. 2. For applications in which attractive interactions are of less concern than the excluded volume effects which dictate molecular packing, the potential may be truncated at the position of its minimum, and shifted upwardsto give what is usually termed the WCA model8 . If electrostatic charges are present, we add the appropriate Coulomb potentials v Coulomb(r) = Q1Q2 4π0r , (4) where Q1, Q2 are the charges and 0 is the permittivity of free space. The correct handling of long-range forces in a simulation is an essential aspect of polyelectrolyte simulations, which will be the subject of the later chapter of Holm9 . 2.2 Bonding Potentials For molecular systems, we simply build the molecules out of site-site potentials of the form of Eq. (3) or similar. Typically, a single-molecule quantum-chemical calculation may be used to estimate the electron density throughout the molecule, which may then be mod￾elled by a distribution of partial charges via Eq. (4), or more accurately by a distribution of electrostatic multipoles4, 10 . For molecules we must also consider the intramolecular bond￾ing interactions. The simplest molecular model will include terms of the following kind: 3

中1234 234 3 Figure 3.Geometry of a simple chain molecule,illustrating the definition of interatomic distance r23.bend angle 0234.and torsion angle 1234. Uintramolecular= ∑5w-r网 (5a) bonds +互∑号9-)2 (5b) bend angles +号∑∑g量1+smo,u-n》 (5c) mem The geometry is illustrated in Fig.3.The "bonds"will typically involve the separation r;=r:-ri between adjacent pairs of atoms in a molecular framework,and we assume in Eq.(5a)a harmonic form with specified equilibrium separation,although this is not the only possibility.The"bend angles"k are between successive bond vectors such as ri-rj and rj-rk,and therefore involve three atom coordinates: cosk==()12()2() where=r/r.Usually this bending term is taken to be quadratic in the angular dis- placement from the equilibrium value,as in Eq.(5b),although periodic functions are also used.The"torsion angles"are defined in terms of three connected bonds,hence four atomic coordinates: cos ijkl=-元k·元jkl,where nijk=T)XTjk,njkl=Tjk X Tkl, and n=n/n,the unit normal to the plane defined by each pair of bonds.Usually the torsional potential involves an expansion in periodic functions of order m =1,2,.... Eq.(5c). A simulation package force-field will specify the precise form of Eq.(5),and the var- ious strength parameters k and other constants therein.Actually,Eq.(5)is a consider- able oversimplification.Molecular mechanics force-fields,aimed at accurately predicting

1 3 φ 1234 234 4 2 θ 23 r Figure 3. Geometry of a simple chain molecule, illustrating the definition of interatomic distance r23, bend angle θ234, and torsion angle φ1234. Uintramolecular = 1 2 X bonds k r ij ￾ rij − req
2 (5a) + 1 2 X bend angles k θ ijk ￾ θijk − θeq
2 (5b) + 1 2 X torsion angles X m k φ, m ijkl ￾ 1 + cos(mφijkl − γm)
(5c) The geometry is illustrated in Fig. 3. The “bonds” will typically involve the separation rij = |ri −rj | between adjacent pairs of atoms in a molecular framework, and we assume in Eq. (5a) a harmonic form with specified equilibrium separation, although this is not the only possibility. The “bend angles” θijk are between successive bond vectors such as ri − rj and rj − rk, and therefore involve three atom coordinates: cos θijk = rˆij · rˆjk = ￾ rij · rij
−1/2￾ rjk · rjk
−1/2￾ rij · rjk
where rˆ = r/r. Usually this bending term is taken to be quadratic in the angular dis￾placement from the equilibrium value, as in Eq. (5b), although periodic functions are also used. The “torsion angles” φijkl are defined in terms of three connected bonds, hence four atomic coordinates: cos φijkl = −nˆijk · nˆjkl , where nijk = rij × rjk , njkl = rjk × rkl , and nˆ = n/n, the unit normal to the plane defined by each pair of bonds. Usually the torsional potential involves an expansion in periodic functions of order m = 1, 2, . . ., Eq. (5c). A simulation package force-field will specify the precise form of Eq. (5), and the var￾ious strength parameters k and other constants therein. Actually, Eq. (5) is a consider￾able oversimplification. Molecular mechanics force-fields, aimed at accurately predicting 4

structures and properties,will include many cross-terms (e.g.stretch-bend):MM311-13 and MM414-16 are examples.Quantum mechanical calculations may give a guide to the"best" molecular force-field;also comparison of simulation results with thermophysical proper- ties and vibration frequencies is invaluable in force-field development and refinement.A separate family of force fields,such as AMBER7.18,CHARMM19 and OPLS20 are geared more to larger molecules(proteins,polymers)in condensed phases;their functional form is simpler,closer to that of Eq.(5),and their parameters are typically determined by quan- tum chemical calculations combined with thermophysical and phase coexistence data.This field is too broad to be reviewed here;several molecular modelling texts21-23(albeit target- ted at biological applications)should be consulted by the interested reader.The modelling 200 150 FENE+WCA 50 FENE ··WCA --harmonic 0 0.5 1.5 r/σ Figure 4.The FENE+WCA potential,with its separate FENE (attractive)and WCA (repulsive)components, between bonded atoms in a coarse-grained polymer chain.Also shown is the equivalent harmonic potential. Unlike the harmonic spring,the FENE potential cannot be extended beyond a specified limit,here Ro=1.5. For more details see Ref.24. of long chain molecules will be of particular interest to us,especially as an illustration of the scope for progressively simplifying and"coarse-graining"the potential model.Var- ious explicit-atom potentials have been devised for the n-alkanes25.More approximate potentials have also been constructed26-28 in which the CH2 and CHa units are represented by single "united atoms".These potentials are typically less accurate and less transfer- able than the explicit-atom potentials,but significantly less expensive;comparisons have been made between the two approaches29.For more complicated molecules this approach may need to be modified.In the liquid crystal field,for instance,a compromise has been suggested30:use the united-atom approach for hydrocarbon chains,but model phenyl ring hydrogens explicitly. In polymer simulations,there is frequently a need to economize further and coarse- grain the interactions more dramatically:significant progress has been made in recent years in approaching this problem systematically31.32.Finally,the most fundamental prop- erties,such as the entanglement length in a polymer melt33,may be investigated using a 5

structures and properties, will include many cross-terms (e.g. stretch-bend): MM311–13 and MM414–16 are examples. Quantum mechanical calculations may give a guide to the “best” molecular force-field; also comparison of simulation results with thermophysical proper￾ties and vibration frequencies is invaluable in force-field development and refinement. A separate family of force fields, such as AMBER17, 18 , CHARMM19 and OPLS20 are geared more to larger molecules (proteins, polymers) in condensed phases; their functional form is simpler, closer to that of Eq. (5), and their parameters are typically determined by quan￾tum chemical calculations combined with thermophysical and phase coexistence data. This field is too broad to be reviewed here; several molecular modelling texts21–23 (albeit target￾ted at biological applications) should be consulted by the interested reader. The modelling 0 0.5 1 1.5 r / σ 0 50 100 150 200 v / ε FENE+WCA FENE WCA harmonic Figure 4. The FENE+WCA potential, with its separate FENE (attractive) and WCA (repulsive) components, between bonded atoms in a coarse-grained polymer chain. Also shown is the equivalent harmonic potential. Unlike the harmonic spring, the FENE potential cannot be extended beyond a specified limit, here R0 = 1.5σ. For more details see Ref. 24. of long chain molecules will be of particular interest to us, especially as an illustration of the scope for progressively simplifying and “coarse-graining” the potential model. Var￾ious explicit-atom potentials have been devised for the n-alkanes25 . More approximate potentials have also been constructed26–28 in which the CH2 and CH3 units are represented by single “united atoms”. These potentials are typically less accurate and less transfer￾able than the explicit-atom potentials, but significantly less expensive; comparisons have been made between the two approaches29 . For more complicated molecules this approach may need to be modified. In the liquid crystal field, for instance, a compromise has been suggested30: use the united-atom approach for hydrocarbon chains, but model phenyl ring hydrogens explicitly. In polymer simulations, there is frequently a need to economize further and coarse￾grain the interactions more dramatically: significant progress has been made in recent years in approaching this problem systematically31, 32 . Finally, the most fundamental prop￾erties, such as the entanglement length in a polymer melt33 , may be investigated using a 5

simple chain of pseudo-atoms or beads(modelled using the WCA potential of Fig.2,and each representing several monomers),joined by an attractive finitely-extensible non-linear elastic (FENE)potential24 which is illustrated in Fig.4. FENE(r）= -kRo In(1-(r/Ro)2)r<Ro (6) T≥Ro The key feature of this potential is that it cannot be extended beyond r=Ro,ensuring(for suitable choices of the parameters k and Ro)that polymer chains cannot move through one another. 2.3 Force Calculation Having specified the potential energy function (rN).the next step is to calculate the atomic forces =-0uy For site-site potentials this is a simple exercise.For the intramolecular part of the potential, it is a little more involved,but still a relatively straightforward application of the chain rule.Examples of how to do it are given in appendix C of Ref.34.As a simple illustration, consider one of the bending potential terms for the polymer of Fig.3.supposing that it may be written v=-kc0s0234=-k(r23·r23)-1/2(r34·r34)-1/2(r23·r34) This will contribute to the forces on all three atoms.To calculate these,we need: 0 -(r23·r34）=T34 0r -(r23T34)=T23-T34 0 O -(r23·T34)=-T23 gr23=2r2 分r23·r23)=-2r23 mr23r2a)=0 元4r3=0 8r94·r34)=2r3别 [r94r3a)=-2r3 and hence cos0234=r23r3 T23·T3 0r T34- T -T23 0 ∂r c0s0234=r23r3 T23:T34r23 T T23T34r31+T23-T34 T34 COs0234=23 T34 T23·T34 0r4 2T34-T23 T34 A similar approach applied to the torsional potential gives the forces on all four involved atoms. 6

simple chain of pseudo-atoms or beads (modelled using the WCA potential of Fig. 2, and each representing several monomers), joined by an attractive finitely-extensible non-linear elastic (FENE) potential24 which is illustrated in Fig. 4. v FENE(r) = ( − 1 2 kR2 0 ln￾ 1 − (r/R0) 2
r < R0 ∞ r ≥ R0 (6) The key feature of this potential is that it cannot be extended beyond r = R0, ensuring (for suitable choices of the parameters k and R0) that polymer chains cannot move through one another. 2.3 Force Calculation Having specified the potential energy function U(r N ), the next step is to calculate the atomic forces fi = − ∂ ∂ri U(r N ) For site-site potentials this is a simple exercise. For the intramolecular part of the potential, it is a little more involved, but still a relatively straightforward application of the chain rule. Examples of how to do it are given in appendix C of Ref. 34. As a simple illustration, consider one of the bending potential terms for the polymer of Fig. 3, supposing that it may be written v = −k cos θ234 = −k(r23 · r23) −1/2 (r34 · r34) −1/2 (r23 · r34) This will contribute to the forces on all three atoms. To calculate these, we need: ∂ ∂r2 (r23 · r34) = r34 ∂ ∂r3 (r23 · r34) = r23 − r34 ∂ ∂r4 (r23 · r34) = −r23 ∂ ∂r2 (r23 · r23) = 2r23 ∂ ∂r3 (r23 · r23) = −2r23 ∂ ∂r4 (r23 · r23) = 0 ∂ ∂r2 (r34 · r34) = 0 ∂ ∂r3 (r34 · r34) = 2r34 ∂ ∂r4 (r34 · r34) = −2r34 and hence ∂ ∂r2 cos θ234 = r −1 23 r −1 34  r34 − r23 · r34 r 2 23 r23 ∂ ∂r3 cos θ234 = r −1 23 r −1 34  r23 · r34 r 2 23 r23 − r23 · r34 r 2 34 r34 + r23 − r34 ∂ ∂r4 cos θ234 = r −1 23 r −1 34  r23 · r34 r 2 34 r34 − r23 A similar approach applied to the torsional potential gives the forces on all four involved atoms. 6

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) 7

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 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

After step (8b).a force evaluation is carried out,to give fi(t+ot)for step (8c).This scheme advances the coordinates and momenta over a timestep 6t.A piece of pseudo-code illustrates how this works: do step =1,nstep pp 0.5*dt*f r =r dt*p/m f force(r) p =p +0.5*dt*f enddo As we shall see shortly there is an interesting theoretical derivation of this version of the algorithm.Important features of the Verlet algorithm are:(a)it is exactly time reversible; (b)it is symplectic (to be discussed shortly);(c)it is low order in time,hence permitting long timesteps:(d)it requires just one(expensive)force evaluation per step:(e)it is easy to program. 3.2 Constraints It is quite common practice in classical computer simulations not to attempt to represent intramolecular bonds by terms in the potential energy function,because these bonds have very high vibration frequencies (and arguably should be treated in a quantum mechanical way rather than in the classical approximation).Instead,the bonds are treated as being constrained to have fixed length.In classical mechanics,constraints are introduced through the Lagrangian5 or Hamiltonian46 formalisms.Given an algebraic relation between two atomic coordinates,for example a fixed bond length b between atoms 1 and 2,one may write a constraint equation,plus an equation for the time derivative of the constraint X(r1,T2)=(r1-T2)·(r1-T2)-b2=0 (9a) X(T1,T2)=2(v1-v2)·(r1-T2)=0. (9b) In the Lagrangian formulation,the constraint forces acting on the atoms will enter thus: miri=fi+Agi where A is the undetermined multiplier and 0y=-2(r1-r2) g1=-r1 _0X=2(r1-r2) 92=-0r2 It is easy to derive an exact expression for the multiplier A from the above equations; if several constraints are imposed,a system of equations(one per constraint)is obtained. However,this exact solution is not what we want:in practice,since the equations of motion are only solved approximately.in discrete time steps,the constraints will be increasingly violated as the simulation proceeds.The breakthrough in this area came with the pro- posal to determine the constraint forces in such a way that the constraints are satisfied exactly at the end of each time step47-49.For the original Verlet algorithm,this scheme is called SHAKE.The appropriate version of this scheme for the velocity Verlet algorithm is called RATTLES0.Formally,we wish to solve the following scheme,in which we combine 8

After step (8b), a force evaluation is carried out, to give f i (t + δt) for step (8c). This scheme advances the coordinates and momenta over a timestep δt. A piece of pseudo-code illustrates how this works: do step = 1, nstep p = p + 0.5*dt*f r = r + dt*p/m f = force(r) p = p + 0.5*dt*f enddo As we shall see shortly there is an interesting theoretical derivation of this version of the algorithm. Important features of the Verlet algorithm are: (a) it is exactly time reversible; (b) it is symplectic (to be discussed shortly); (c) it is low order in time, hence permitting long timesteps; (d) it requires just one (expensive) force evaluation per step; (e) it is easy to program. 3.2 Constraints It is quite common practice in classical computer simulations not to attempt to represent intramolecular bonds by terms in the potential energy function, because these bonds have very high vibration frequencies (and arguably should be treated in a quantum mechanical way rather than in the classical approximation). Instead, the bonds are treated as being constrained to have fixed length. In classical mechanics, constraints are introduced through the Lagrangian45 or Hamiltonian46 formalisms. Given an algebraic relation between two atomic coordinates, for example a fixed bond length b between atoms 1 and 2, one may write a constraint equation, plus an equation for the time derivative of the constraint χ(r1, r2) = (r1 − r2) · (r1 − r2) − b 2 = 0 (9a) χ˙(r1, r2) = 2(v1 − v2) · (r1 − r2) = 0 . (9b) In the Lagrangian formulation, the constraint forces acting on the atoms will enter thus: mir¨i = fi + Λgi where Λ is the undetermined multiplier and g1 = − ∂χ ∂r1 = −2(r1 − r2) g2 = − ∂χ ∂r2 = 2(r1 − r2) It is easy to derive an exact expression for the multiplier Λ from the above equations; if several constraints are imposed, a system of equations (one per constraint) is obtained. However, this exact solution is not what we want: in practice, since the equations of motion are only solved approximately, in discrete time steps, the constraints will be increasingly violated as the simulation proceeds. The breakthrough in this area came with the pro￾posal to determine the constraint forces in such a way that the constraints are satisfied exactly at the end of each time step47–49 . For the original Verlet algorithm, this scheme is called SHAKE. The appropriate version of this scheme for the velocity Verlet algorithm is called RATTLE50 . Formally, we wish to solve the following scheme, in which we combine 8

(r1,r2)into r,(p,p2)into p,etc.for simplicity: p(t+8t)=p(t)+tf(t)+xg(t) r(t+6t)=r(t)8tp(t+t)/m choosing A such that:0=x(r(t +6t)) (10a) p(t+t)=p(t+号6t)+6tf(t+6t)+μg(t+6t) choosing u such that:0=(r(t+6t),p(t +6t)) (10b) Step(10a)may be implemented by defining unconstrained variables (t+26t)=p(t)+26tf(t),(t+6t)=r(t)+6tpt+2t)/m then solving the nonlinear equation for A X(t+6t)=x((t+t)+λ6tg(t)/m)=0 and substituting back p(t+t)=p(t+t)+xg(t);r(t+6t)=(t+6t)+6tAg(t)/m Step(10b)may be handled by defining p(t+t)=p(t+号t)+号tf(t+6t) solving the equation for the second Lagrange multiplier u X(t+6t)=(r(t+6t),(t+t)+g(t+6t)）=0 (which is actually linear,since x(r,p)=-g(r).p/m)and substituting back p(t+6t)=(6t)+ug(t+t) In pseudo-code this scheme may be written do step 1,nstep pp (dt/2)*f r =r dt*p/m lambda_g shake(r) pp lambda g r=r dt*lambda_g/m f force(r) p =p (dt/2)*f mu_g rattle(r,p) p =p mug enddo The routine called shake here calculates the constraint forces Ag;necessary to ensure that the end-of-step positions ri satisfy Eq.(9a).For a system of many constraints,this calculation is usually performed in an iterative fashion,so as to satisfy each constraint in turn until convergence;the original SHAKE algorithm was framed in this way.These constraint forces are incorporated into both the end-of-step positions and the mid-step mo- menta.The routine called rattle calculates a new set of constraint forces ug;to ensure that the end-of-step momenta satisfy Eq.(9b).This also may be carried out iteratively

(r1, r2) into r, (p1 , p2 ) into p, etc. for simplicity: p(t + 1 2 δt) = p(t) + 1 2 δtf(t) + λg(t) r(t + δt) = r(t) + δtp(t + 1 2 δt)/m choosing λ such that: 0 = χ(r(t + δt)) (10a) p(t + δt) = p(t + 1 2 δt) + 1 2 δtf(t + δt) + µg(t + δt) choosing µ such that: 0 = χ˙(r(t + δt), p(t + δt)) (10b) Step (10a) may be implemented by defining unconstrained variables p¯(t + 1 2 δt) = p(t) + 1 2 δtf(t), r¯(t + δt) = r(t) + δtp¯(t + 1 2 δt)/m then solving the nonlinear equation for λ χ(t + δt) = χ ￾ r¯(t + δt) + λδtg(t)/m
= 0 and substituting back p(t + 1 2 δt) = p¯(t + 1 2 δt) + λg(t), r(t + δt) = r¯(t + δt) + δtλg(t)/m Step (10b) may be handled by defining p¯(t + δt) = p(t + 1 2 δt) + 1 2 δtf(t + δt) solving the equation for the second Lagrange multiplier µ χ˙(t + δt) = χ˙ ￾ r(t + δt), p¯(t + δt) + µg(t + δt)
= 0 (which is actually linear, since χ˙(r, p) = −g(r) · p/m) and substituting back p(t + δt) = p¯(δt) + µg(t + δt) In pseudo-code this scheme may be written do step = 1, nstep p = p + (dt/2)*f r = r + dt*p/m lambda_g = shake(r) p = p + lambda_g r = r + dt*lambda_g/m f = force(r) p = p + (dt/2)*f mu_g = rattle(r,p) p = p + mu_g enddo The routine called shake here calculates the constraint forces λgi necessary to ensure that the end-of-step positions ri satisfy Eq. (9a). For a system of many constraints, this calculation is usually performed in an iterative fashion, so as to satisfy each constraint in turn until convergence; the original SHAKE algorithm was framed in this way. These constraint forces are incorporated into both the end-of-step positions and the mid-step mo￾menta. The routine called rattle calculates a new set of constraint forces µgi to ensure that the end-of-step momenta satisfy Eq. (9b). This also may be carried out iteratively. 9