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 assemblies 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); additionally, 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 limitations 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 experiment. 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 essential. 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 ) represents 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 potentials 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 continuous, 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 modelled 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 bonding 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 − req2 (5a) + 1 2 X bend angles k θ ijk θijk − θeq2 (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 displacement 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 various strength parameters k and other constants therein. Actually, Eq. (5) is a considerable oversimplification. Molecular mechanics force-fields, aimed at accurately predicting 4