John von Neumann Institute for Computing NIC Classical Molecular Dynamics Godehard Sutmann published in Quantum Simulations of Complex Many-Body Systems: From Theory to Algorithms,Lecture Notes, J.Grotendorst,D.Marx,A.Muramatsu(Eds.), John von Neumann Institute for Computing,Julich, NIC Series,Vol.10,1SBN3-00-009057-6,pp.211-254,2002. C2002 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/volume10
John von Neumann Institute for Computing Classical Molecular Dynamics Godehard Sutmann published in Quantum Simulations of Complex Many-Body Systems: From Theory to Algorithms, Lecture Notes, J. Grotendorst, D. Marx, A. Muramatsu (Eds.), John von Neumann Institute for Computing, Julich, ¨ NIC Series, Vol. 10, ISBN 3-00-009057-6, pp. 211-254, 2002. c 2002 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/volume10
Classical Molecular Dynamics Godehard Sutmann John von Neumann Institute for Computing Central Institute for Applied Mathematics Research Centre Julich,52425 Julich,Germany E-mail:g.sutmann@fz-juelich.de An introduction to classical molecular dynamics simulation is presented.In addition to some historical notes,an overview is given over particle models,integrators and different ensemble techniques.In the end,methods are presented for parallisation of short range interaction po- tentials.The efficiency and scalability of the algorithms on massively parallel computers is discussed with an extended version of Amdahl's law. 1 Introduction Computer simulation methods have become a very powerful tool to attack the many-body problem in statistical physics.physical chemistry and biophysics.Although the theoretical description of complex systems in the framework of statistical physics is rather well de- veloped and the experimental techniques for detailed microscopic information are rather sophisticated,it is often only possible to study specific aspects of those systems in great de- tail via the simulation.On the other hand,simulations need specific input parameters that characterize the system in question,and which come either from theoretical considerations or are provided by experimental data.Having characterized a physical system in terms of model parameters.simulations are often used both to solve theoretical models beyond certain approximations and to provide a hint to experimentalists for further investigations. In the case of big experimental facilities it is even often required to prove the potential outcome of an experiment by computer simulations.In that way one can say that the field of computer simulations has developed into a very important branch of science,which on the one hand helps theorists and experimentalists to go beyond their inherent limitations and on the other hand is a scientific field on its own. The traditional simulation methods for many-body systems can be divided into two classes of stochastic and deterministic simulations.which are largely covered by the Monte Carlo(MC)method and the molecular dynamics(MD)method,respectively.Monte Carlo simulations probe the configuration space by trial moves of particles.Within the so-called Metropolis algorithm,the energy change from step n to n+1 is used as a trigger to accept or reject the new configuration.Paths towards lower energy are always accepted,those to higher energy are accepted with a probability governed by Boltzmann statistics.In that way,properties of the system can be calculated by averaging over all Monte Carlo moves (where one move means that every degree of freedom is probed once on average). By contrast,MD methods are governed by the system's Hamiltonian and consequently Hamilton's equations of motion OH OH p=- qi= Opi (1) dgi 211
Classical Molecular Dynamics Godehard Sutmann John von Neumann Institute for Computing Central Institute for Applied Mathematics Research Centre Julich, ¨ 52425 Julich, ¨ Germany E-mail: g.sutmann@fz-juelich.de An introduction to classical molecular dynamics simulation is presented. In addition to some historical notes, an overview is given over particle models, integrators and different ensemble techniques. In the end, methods are presented for parallisation of short range interaction potentials. The efficiency and scalability of the algorithms on massively parallel computers is discussed with an extended version of Amdahl’s law. 1 Introduction Computer simulation methods have become a very powerful tool to attack the many-body problem in statistical physics, physical chemistry and biophysics. Although the theoretical description of complex systems in the framework of statistical physics is rather well developed and the experimental techniques for detailed microscopic information are rather sophisticated, it is often only possible to study specific aspects of those systems in great detail via the simulation. On the other hand, simulations need specific input parameters that characterize the system in question, and which come either from theoretical considerations or are provided by experimental data. Having characterized a physical system in terms of model parameters, simulations are often used both to solve theoretical models beyond certain approximations and to provide a hint to experimentalists for further investigations. In the case of big experimental facilities it is even often required to prove the potential outcome of an experiment by computer simulations. In that way one can say that the field of computer simulations has developed into a very important branch of science, which on the one hand helps theorists and experimentalists to go beyond their inherent limitations and on the other hand is a scientific field on its own. The traditional simulation methods for many-body systems can be divided into two classes of stochastic and deterministic simulations, which are largely covered by the Monte Carlo (MC) method and the molecular dynamics (MD) method, respectively. Monte Carlo simulations probe the configuration space by trial moves of particles. Within the so-called Metropolis algorithm, the energy change from step n to n + 1 is used as a trigger to accept or reject the new configuration. Paths towards lower energy are always accepted, those to higher energy are accepted with a probability governed by Boltzmann statistics. In that way, properties of the system can be calculated by averaging over all Monte Carlo moves (where one move means that every degree of freedom is probed once on average). By contrast, MD methods are governed by the system’s Hamiltonian and consequently Hamilton’s equations of motion p˙i = − ∂H ∂qi , q˙i = ∂H ∂pi (1) 211
are integrated to move particles to new positions and to get new velocities at these new positions.This is an advantage of MD simulations with respect to MC,since not only is the configuration space probed but the whole phase space which gives additional information about the dynamics of the system.Both methods are complementary in nature but they lead to the same averages of static quantities,given that the system under consideration is ergodic and the same statistical ensemble is used. Although there are different methods to obtain information about complex systems, particle simulations always require a model for the interaction between system constitu- ants.This model has to be tested against experimental results,i.e.it should reproduce or approximate experimental findings like distribution functions or phase diagrams,and theoretical constraints,i.e.it should obey certain fundamental or limiting laws like energy conservation. Concerning MD simulations the ingredients for a program are basically threefold: (i)As already mentioned,a model for the interaction between system constituants (atoms, molecules.surfaces etc.)is needed.Often,it is assumed that particles interact only pair- wise,which is exact e.g.for particles with fixed partial charges.This assumption greatly reduces the computational effort and the work to implement the model into the program. (ii)An integrator is needed,which propagates particle positions and velocities from time t to t ot.It is a finite difference scheme which moves trajectories discretely in time.The time step ot has properly to be chosen to guarantee stability of the integrator,i.e.there should be no drift in the system's energy (iii)A statistical ensemble has to be chosen,where thermodynamic quantities like pressure, temperature or the number of particles are controlled.The natural choice of an ensemble in MD simulations is the microcanonical ensemble (NVE),since the system's Hamiltonian without external potentials is a conserved quantity.Nevertheless,there are extensions to the Hamiltonian which also allow to simulate different statistical ensembles. These steps essentially define an MD simulation.Having this tool at hand,it is possible to obtain exact results within numerical precision.Results are only correct with respect to the model which enters into the simulation and they have to be tested against theoretical predictions and experimental findings.If the simulation results differ from the real system properties or are incompatible with solid theoretical manifestations,the model has to be refined.This procedure can be understood as an adaptive refinement which leads in the end to an approximation of a model of the real world at least for certain properties.The model itself may be constructed from plausible considerations,where parameters are cho- sen from neutron diffraction or NMR measurements.It may also result from first principle investigations,like quantum ab initio calculations.Although the electronic distribution of the particles is calculated very accurately,this type of model building contains also some approximations,since many-body interactions are mostly neglected(this would increase the parameter space in the model calculation enormously).However,it often provides a good starting point for a realistic model. An important issue of simulation studies is the accessible time-and length-scale cov- erable by microscopic simulations.Fig.1 shows a schematic representation for different types of simulations in a length-time-diagram.It is clear that the more detailed a simu- lation technique operates,the smaller is the accessibility of long times and large length scales.Therefore quantum simulations,where fast motions of electrons are taken into account,are located in the lower left corner of the diagram and typical length and time 212
are integrated to move particles to new positions and to get new velocities at these new positions. This is an advantage of MD simulations with respect to MC, since not only is the configuration space probed but the whole phase space which gives additional information about the dynamics of the system. Both methods are complementary in nature but they lead to the same averages of static quantities, given that the system under consideration is ergodic and the same statistical ensemble is used. Although there are different methods to obtain information about complex systems, particle simulations always require a model for the interaction between system constituants. This model has to be tested against experimental results, i.e. it should reproduce or approximate experimental findings like distribution functions or phase diagrams, and theoretical constraints, i.e. it should obey certain fundamental or limiting laws like energy conservation. Concerning MD simulations the ingredients for a program are basically threefold: (i) As already mentioned, a model for the interaction between system constituants (atoms, molecules, surfaces etc.) is needed. Often, it is assumed that particles interact only pairwise, which is exact e.g. for particles with fixed partial charges. This assumption greatly reduces the computational effort and the work to implement the model into the program. (ii) An integrator is needed, which propagates particle positions and velocities from time t to t + δt. It is a finite difference scheme which moves trajectories discretely in time. The time step δt has properly to be chosen to guarantee stability of the integrator, i.e. there should be no drift in the system’s energy. (iii) A statistical ensemble has to be chosen, where thermodynamic quantitieslike pressure, temperature or the number of particles are controlled. The natural choice of an ensemble in MD simulations is the microcanonical ensemble (NVE), since the system’s Hamiltonian without external potentials is a conserved quantity. Nevertheless, there are extensions to the Hamiltonian which also allow to simulate different statistical ensembles. These steps essentially define an MD simulation. Having this tool at hand, it is possible to obtain exact results within numerical precision. Results are only correct with respect to the model which enters into the simulation and they have to be tested against theoretical predictions and experimental findings. If the simulation results differ from the real system properties or are incompatible with solid theoretical manifestations, the model has to be refined. This procedure can be understood as an adaptive refinement which leads in the end to an approximation of a model of the real world at least for certain properties. The model itself may be constructed from plausible considerations, where parameters are chosen from neutron diffraction or NMR measurements. It may also result from first principle investigations, like quantum ab initio calculations. Although the electronic distribution of the particles is calculated very accurately, this type of model building contains also some approximations, since many-body interactions are mostly neglected (this would increase the parameter space in the model calculation enormously). However, it often provides a good starting point for a realistic model. An important issue of simulation studies is the accessible time- and length-scale coverable by microscopic simulations. Fig.1 shows a schematic representation for different types of simulations in a length-time-diagram. It is clear that the more detailed a simulation technique operates, the smaller is the accessibility of long times and large length scales. Therefore quantum simulations, where fast motions of electrons are taken into account, are located in the lower left corner of the diagram and typical length and time 212
10000 HD ● 1000 100 BD 10 MD QS fs ps ns us ms time Figure 1.Schematic comparison of time-and length-scales,accessible to different types of simulation techniques(quantum simulations(QM),molecular dynamics (MD),Brow- nian dynamics (BD)and hydrodynamics/fluid dynamics(HD)).The black dots mark the longest(≈1s)and the biggest(N>5×109,L≈0.4 m molecular dynamics simulations by Duan Kollman'and Roth2 respectively.) scales are of order of A and ps.Classical molecular dynamics approximates electronic distributions in a rather coarse-grained fashion by putting either fixed partial charges on interaction sites or by adding an approximate model for polarization effects.In both cases, the time scale of the system is not dominated by the motion of electrons,but the time of intermolecular collision events.rotational motions or intramolecular vibrations,which are orders of magnitude slower than those of electron motions.Consequently,the time step of integration is larger and trajectory lengths are of order ns and accessible lengths of order 10-100 A.If one considers tracer particles in a solvent medium,where one is not inter- ested in a detailed description of the solvent,one can apply Brownian dynamics,where the effect of the solvent is hidden in average quantities.Since collision times between tracer particles is very long,one may apply larger timesteps.Furthermore,since the solvent is not simulated explicitly,the lengthscales may be increased considerably.Finally,if one is in- terested not in a microscopic picture of the simulated system but in macroscopic quantities the concepts of hydrodynamics may be applied,where the system properties are hidden in effective numbers,e.g.density,viscosity,sound velocity. It is clear that the performance of particle dynamics simulations strongly depends on the computer facilities at hand.The first studies using MD simulation techniques were performed in 1957 by B.J.Alder and T.E.Wainright3 who simulated the phase transition 213
fs ps ns µs ms 1 10 100 1000 10000 QS MD BD QS MD BD HD length [ Å] time Figure 1. Schematic comparison of time- and length-scales, accessible to different types of simulation techniques (quantum simulations (QM), molecular dynamics (MD), Brownian dynamics (BD) and hydrodynamics/fluid dynamics (HD)). The black dots mark the longest (≈ 1 µs) and the biggest (N > 5 × 109 , L ≈ 0.4µm molecular dynamics simulations by Duan & Kollman1 and Roth2 respectively.) scales are of order of A˚ and ps. Classical molecular dynamics approximates electronic distributions in a rather coarse-grained fashion by putting either fixed partial charges on interaction sites or by adding an approximate model for polarization effects. In both cases, the time scale of the system is not dominated by the motion of electrons, but the time of intermolecular collision events, rotational motions or intramolecular vibrations, which are orders of magnitude slower than those of electron motions. Consequently, the time step of integration is larger and trajectory lengths are of order ns and accessible lengths of order 10 − 100 A. ˚ If one considers tracer particles in a solvent medium, where one is not interested in a detailed description of the solvent, one can apply Brownian dynamics, where the effect of the solvent is hidden in average quantities. Since collision times between tracer particles is very long, one may apply larger timesteps. Furthermore,since the solvent is not simulated explicitly, the lengthscales may be increased considerably. Finally, if one is interested not in a microscopic picture of the simulated system but in macroscopic quantities, the concepts of hydrodynamics may be applied, where the system properties are hidden in effective numbers, e.g. density, viscosity, sound velocity. It is clear that the performance of particle dynamics simulations strongly depends on the computer facilities at hand. The first studies using MD simulation techniques were performed in 1957 by B. J. Alder and T. E. Wainright3 who simulated the phase transition 213
of a system of hard spheres.The general method,however,was presented two years later" In this early simulation.which was run an an IBM-704,up to 500 particles could be simu- lated,for which 500 collisions per hour could be calculated.Taking into account 200000 collisions for a production run,these simulations lasted for more than two weeks.The propagation of hard spheres in a simulation is determined by the collision events between two particles.Therefore,the propagation is not based on an integration of the equations of motion.but rather the calculation of the time of the next collision.which results in a variable time step in the calculations The first MD simulation which was applied to atoms interacting via a continuous po- tential was performed by A.Rahman in 1964.In this case,a model system for Argon was simulated and not only binary collisions were taken into account but the interactions were modeled by a Lennard-Jones potential and the equations of motion were integrated with a finite difference scheme.This work may be considered as seminal for dynamical calcula- tions.It was the first work where an exact method (within numerical precision)was used to calculate dynamical quantities like autocorrelation functions and transport coefficients like the diffusion coefficient for a realistic system.Also more involved topics like the dy- namic van Hove function and non-Gaussian corrections to diffusion were evaluated.The calculations were performed for 864 particles on a CDC 3600,where the propagation of all particles for one time step took 45 s.The calculation of 50000 timesteps then took more than three weeks!a With the development of faster and bigger massively parallel architectures the accessi- ble time and length scales are increasing.In the case of classical MD simulations it was demonstrated by J.Roth in 1999 on the CRAY T3E-1200 in Juilich that it is possible to simulate more than 5 x 109 particles,corresponding to a length scale of several 1000 A. This was possible with the highly memory optimised MD program IMD5.2,which used the 512 nodes with 256 MB memory each,quite efficiently.However,the limits of such a demonstration became rather obvious,since for a usual production run of 10000 time steps a simulation time of a quarter of a year would be required(given that the whole machine is dedicated to one user).In another demonstration run Y.Duan and P.A.Kollman extended the time scale of an all atom MD simulation to 1 us,where they simulated the folding process of the subdomain HP-36 from the villin headpiece.1.The protein was modelled with a 596 interaction site model dissolved in a system of 3000 water molecules.Using a timestep of integration of 2 x 10-15s,the program was run for 5 x 108 steps.In order to perform this type of calculation,it was necessary to run the program several months on a CRAY T3D and CRAY T3E with 256 processors.It is clear that such kind of simulation is exceptional due to the large amount of computer resources needed,but is nonetheless a kind of milestone pointing to future simulation practices. Classical molecular dynamics methods are nowadays applied to a huge class of prob- lems,e.g.properties of liquids,defects in solids,fracture,surface properties,friction, molecular clusters,polyelectrolytes and biomolecules.Due to the large area of applica- bility,simulation codes for molecular dynamics were developed by many groups.On the internet homepage of the Collaborative Computational Project No.5(CCP5)?a lot of com- puter codes are assembled for condensed phase dynamics.During the last years several programs were designed for parallel computers.Among them,which are partly avail- "On a standard PC this calculation may be done within one hour nowadays! 214
of a system of hard spheres. The general method, however, was presented two years later4 . In this early simulation, which was run an an IBM-704, up to 500 particles could be simulated, for which 500 collisions per hour could be calculated. Taking into account 200000 collisions for a production run, these simulations lasted for more than two weeks. The propagation of hard spheres in a simulation is determined by the collision events between two particles. Therefore, the propagation is not based on an integration of the equations of motion, but rather the calculation of the time of the next collision, which results in a variable time step in the calculations. The first MD simulation which was applied to atoms interacting via a continuous potential was performed by A. Rahman in 1964. In this case, a model system for Argon was simulated and not only binary collisions were taken into account but the interactions were modeled by a Lennard-Jones potential and the equations of motion were integrated with a finite difference scheme. This work may be considered as seminal for dynamical calculations. It was the first work where an exact method (within numerical precision) was used to calculate dynamical quantities like autocorrelation functions and transport coefficients like the diffusion coefficient for a realistic system. Also more involved topics like the dynamic van Hove function and non-Gaussian corrections to diffusion were evaluated. The calculations were performed for 864 particles on a CDC 3600, where the propagation of all particles for one time step took ≈ 45 s. The calculation of 50000 timesteps then took more than three weeks! a With the development of faster and bigger massively parallel architectures the accessible time and length scales are increasing. In the case of classical MD simulations it was demonstrated by J. Roth in 1999 on the CRAY T3E-1200 in Julich ¨ that it is possible to simulate more than 5 × 109 particles, corresponding to a length scale of several 1000 A. ˚ This was possible with the highly memory optimised MD program IMD5, 2 , which used the 512 nodes with 256 MB memory each, quite efficiently. However, the limits of such a demonstration became rather obvious, since for a usual production run of 10000 time steps a simulation time of a quarter of a year would be required (given that the whole machine is dedicated to one user). In another demonstration run Y. Duan and P. A. Kollman extended the time scale of an all atom MD simulation to 1 µs, where they simulated the folding process of the subdomain HP-36 from the villin headpiece6, 1 . The protein was modelled with a 596 interaction site model dissolved in a system of 3000 water molecules. Using a timestep of integration of 2 × 10−15s, the program was run for 5 × 108 steps. In order to perform this type of calculation, it was necessary to run the program several months on a CRAY T3D and CRAY T3E with 256 processors. It is clear that such kind of simulation is exceptional due to the large amount of computer resources needed, but is nonetheless a kind of milestone pointing to future simulation practices. Classical molecular dynamics methods are nowadays applied to a huge class of problems, e.g. properties of liquids, defects in solids, fracture, surface properties, friction, molecular clusters, polyelectrolytes and biomolecules. Due to the large area of applicability, simulation codes for molecular dynamics were developed by many groups. On the internet homepage of the Collaborative Computational Project No.5 (CCP5)7 a lot of computer codes are assembled for condensed phase dynamics. During the last years several programs were designed for parallel computers. Among them, which are partly availaOn a standard PC this calculation may be done within one hour nowadays! 214
able free of charge,are,e.g.,Amber/Sander8,CHARMM,NAMD10,NWCHEM!1 and LAMMPS!2」 2 Models for Particle Interactions A system is completely determined through it's Hamiltoian =Hto +H1,where Ho is the internal part of the Hamiltonian,given as N 0=+∑4，)+u,,r)+ (2) =1 id (3) 215
able free of charge, are, e.g., Amber/Sander8 , CHARMM9 , NAMD10 , NWCHEM11 and LAMMPS12 . 2 Models for Particle Interactions A system is completely determined through it’s Hamiltoian H = H0 + H1, where H0 is the internal part of the Hamiltonian, given as H0 = X N i=1 p 2 i 2mi + X N i d (3) 215
i.e.a particles'potential energy gets contributions from all particles of the universe if n d,otherwise the interaction is bound to a certain region,which is often modeled by a spherical interaction range.The long range nature of the interaction becomes most important for potentials which only have potential parameters of the same sign,like the gravitational potential where no screening can occur.For Coulomb energies,where posi- tive and negative charges may compensate each other,long range effects may be of minor importance in some systems like molten salts.In the following two examples shall illus- trate the different treatment of short-and long range interactions. 2.1 Short Range Interactions Short range interactions offer the possibility to take into account only neigbored particles up to a certain distance for the calculation of interactions.In that way a cutoff radius is introduced beyond of which mutual interactions between particles are neglected.As an approximation one may introduce long range corrections to the potential in order to compensate for the neglect of explicit calculations.The whole short range potential may then be written as N U=>u(rijlrig Re)+Uire (4) iRe,which makes it possible for many types of potentials to calculate Uire analytically. Besides internal degrees of freedom of molecules,which may be modeled with short range interaction potentials(c.f.Fig.2),it is first of all the excluded volume of a particle which is of importance.A finite diameter of a particle may be represented by a steep repulsive potential acting at short distances.This is either described by an exponential function or an algebraic form,o r-",where n>9.Another source of short range interaction is the van der Waals interaction.For neutral particles these are the London forces arising from induced dipole interactions.Fluctuations of the electron distribution of a particle give rise to fluctuating dipole moments,which on average compensate to zero.But the instantaneous created dipoles induce also dipoles on neighbored particles which attract each other o r-6.Two common forms of the resulting interactions are the Buckingham potential uag(rij)=Aage-Boary Das r (6) and the Lennard-Jones potential (=)-(=)） (7) 216
i.e. a particles’ potential energy gets contributions from all particles of the universe if n ≤ d, otherwise the interaction is bound to a certain region, which is often modeled by a spherical interaction range. The long range nature of the interaction becomes most important for potentials which only have potential parameters of the same sign, like the gravitational potential where no screening can occur. For Coulomb energies, where positive and negative charges may compensate each other, long range effects may be of minor importance in some systems like molten salts. In the following two examples shall illustrate the different treatment of short- and long range interactions. 2.1 Short Range Interactions Short range interactions offer the possibility to take into account only neigbored particles up to a certain distance for the calculation of interactions. In that way a cutoff radius is introduced beyond of which mutual interactions between particles are neglected. As an approximation one may introduce long range corrections to the potential in order to compensate for the neglect of explicit calculations. The whole short range potential may then be written as U = X N i Rc, which makes it possible for many types of potentials to calculate Ulrc analytically. Besides internal degrees of freedom of molecules, which may be modeled with short range interaction potentials (c.f. Fig.2), it is first of all the excluded volume of a particle which is of importance. A finite diameter of a particle may be represented by a steep repulsive potential acting at short distances. This is either described by an exponential function or an algebraic form, ∝ r −n , where n ≥ 9. Another source of short range interaction is the van der Waals interaction. For neutral particles these are the London forces arising from induced dipole interactions. Fluctuations of the electron distribution of a particle give rise to fluctuating dipole moments, which on average compensate to zero. But the instantaneous created dipoles induce also dipoles on neighbored particles which attract each other ∝ r −6 . Two common forms of the resulting interactions are the Buckingham potential u B αβ(rij ) = Aαβe −Bαβrij − Dαβ r 6 ij (6) and the Lennard-Jones potential u LJ αβ (rij ) = 4αβ σαβ rij 12 − σαβ rij 6 ! (7) 216
6-9 van der Waals -e{e 6-12 van der Waals 6- Electrostatic 4)=99 Quadratic bond-stretching 6)上6-6 Morse bond-stretching u,g)=kh-e-明 Bond-bending 4e,)上®.a-喝 Improper dihedrals ea上naw-d Proper dihedrals )=k (+cos(n) Figure 2.Typical examples for potential terms as used in common force-fields. which are compared in Fig.3.In Eqs.6.7 the indices a,B indicate the species of the particles,i.e.there are parameters A,B,D in Eq.6 and e,o in Eq.7 for intra-species inter- actions (a =B)and cross species interactions (aB).For the Lennard-Jones potential the parameters have a simple physical interpretation:e is the minimum potential energy, located at r=21/6 and o is the diameter of the particle,since forr<o the potential becomes repulsive.Often the Lennard-Jones potential gives a reasonable approximation of a true potential.However,from exact quantum ab initio calculations an exponential type repulsive potential is often more appropriate.Especially for dense systems the too steep repulsive part often leeds to an overestimation of the pressure in the system.Since computationally the Lennard-Jones interaction is quite attractive the repulsive part is some- times replaced by a weaker repulsive term,like r-9.The Lennard-Jones potential has another advantage over the Buckingham potential,since there are combining rules for the parameters.A common choice are the Lorentz-Berelot combining rules 0aB=Jaa +aga 2 EaB VEaaEBB (8) This combining rule is,however,known to overestimate the well depth parameter.Two 217
Proper dihedrals Improper dihedrals Bond-bending Quadratic bond-stretching Morse bond-stretching 6-9 van der Waals 6-12 van der Waals Electrostatic ( ) σ − σ = ε 9 6 4 ij ij ij ij ij ij ij r r u r ( ) σ − σ = ε 12 6 4 ij ij ij ij ij ij ij r r u r ( ) ij i j ij ij r q q u r = bij ( ) ( ) 2 2 1 ij ij ij ij s uij r = k r − b ( ) ( ) ( ) 2 0 1 a r r ij ij ij u r k e − − = − ϑ ( ) ( ) 2 0 2 1 ijk ijk ijk ijk b ij u ϑ = k ϑ − ϑ i j k l ( ) ( ) 2 0 2 1 ξijkl = ijkl ξijkl − ξ id uij k i j k l ( ) ( ( )) 0 ϕijkl = ϕ 1+ cos ϕijkl − ϕ pd uij k n Figure 2. Typical examples for potential terms as used in common force-fields. which are compared in Fig.3. In Eqs.6,7 the indices α, β indicate the species of the particles, i.e. there are parameters A, B, D in Eq.6 and , σ in Eq.7 for intra-species interactions (α = β) and cross species interactions (α 6= β). For the Lennard-Jones potential the parameters have a simple physical interpretation: is the minimum potential energy, located at r = 2 1/6σ and σ is the diameter of the particle, since for r < σ the potential becomes repulsive. Often the Lennard-Jones potential gives a reasonable approximation of a true potential. However, from exact quantum ab initio calculations an exponential type repulsive potential is often more appropriate. Especially for dense systems the too steep repulsive part often leeds to an overestimation of the pressure in the system. Since computationally the Lennard-Jonesinteraction is quite attractive the repulsive part is sometimes replaced by a weaker repulsive term, like ∝ r −9 . The Lennard-Jones potential has another advantage over the Buckingham potential, since there are combining rules for the parameters. A common choice are the Lorentz-Berelot combining rules σαβ = σαα + σββ 2 , αβ = √ ααββ (8) This combining rule is, however, known to overestimate the well depth parameter. Two 217
Buckingham ------Lennard-Jones 12-6 …Lennard-Jones9-6 1 6 r [A] Figure 3.Comparison between a Buckingham-,Lennard-Jones(12-6)and Lennard-Jones (9-6)potential. other commonly known combining rules try to correct this effect,which are Kong24 rules 131 1/6 /12 1 12 e3033 0a6= 213 (9) EaQGGQE38038 (10) and the Waldman-Kagler25 rule 8。+B 1/6 CaaGoaE3303B 0a= (11) 2 Goa In a recent study26 of Ar-Kr and Ar-Ne mixtures,these combining rules were tested and it was found that the Kong rules give the best agreement between simulated and experimental pressure-density curves.An illustration of the different combining rules is shown in Fig.4 for the case of an Ar-Ne mixture. 2.2 Long Range Interactions In the case of long range potentials,like the Coulomb potential,interactions between all particles in the system must be taken into account,if treated without any approximation. 218
2 4 6 8 0 1 2 Buckingham Lennard-Jones 12-6 Lennard-Jones 9-6 U [arb. units] r [Å] Figure 3. Comparison between a Buckingham-, Lennard-Jones (12-6) and Lennard-Jones (9-6) potential. other commonly known combining rules try to correct this effect, which are Kong24 rules σαβ = 1 2 13 αασ 12 q αα αασ 6 ααββσ 6 ββ 1 + ββσ 12 ββ αασ 12 αα !1/13 13 1/6 (9) αβ = q αασ 6 ααββσ 6 ββ σ 6 αβ (10) and the Waldman-Kagler25 rule σαβ = σ 6 αα + σ 6 ββ 2 !1/6 , αβ = q αασ6 ααββσ 6 ββ σ 6 αβ (11) In a recent study26 of Ar-Kr and Ar-Ne mixtures, these combining rules were tested and it was found that the Kong rules give the best agreement between simulated and experimental pressure-density curves. An illustration of the different combining rules is shown in Fig.4 for the case of an Ar-Ne mixture. 2.2 Long Range Interactions In the case of long range potentials, like the Coulomb potential, interactions between all particles in the system must be taken into account, if treated without any approximation. 218
100 一Lorentz-Berthelot 60 ------Kong % …Waldman-Hagler 名 20 -20 -40 60 2 r [A] Figure 4.Resulting cross-terms of the Lennard-Jones potential for an Ar-Ne mixture. Shown is the effect of different combining rules (Eqs.8-11).Parameters used are Ar= 3.406A,eAr=119.4 K and oNe=2.75A,eNe=35.7K. This leads to an (N2)problem,which increases considerably the execution time of a pro- gram for larger systems.For systems with open boundary conditions(like liquid droplets). this method is straightforwardly implemented and reduces to a double sum over all pairs of particles.In the case when periodic boundary conditions are applied,not only the in- teractions with particles in the central cell but also with all periodic images must be taken into account and formally a lattice sum has to be evaluated N U= 1 2 (12) i,j=1 n ri -nLl where n is a lattice vector and means that for n =0 it is ij.It is,however,a well known problem that this type of lattice sum is conditionally convergent,i.e.the result depends on the sequence of evaluating the sum (see e.g.27).A method to overcome this limitation was invented by Ewald28.The idea is to introduce a convergence factor into the sum of Eq.12 which depends on a parameter s,evaluate the sum and put s0 in the end. A characterization of the convergence factors was given in Ref.29.30.A form which leads to the Ewald sum is an exponential en2 transforming Eq.12 into esn (13) i.j=1 n The evaluation and manipulation of this equation proceeds now by using the definition for the I-function and the Jacobi imaginary transform.A very instructive way of the derivation 219
2 4 6 8 -60 -40 -20 0 20 40 60 80 100 Lorentz-Berthelot Kong Waldman-Hagler εAr-Ne [K] r [Å] Figure 4. Resulting cross-terms of the Lennard-Jones potential for an Ar-Ne mixture. Shown is the effect of different combining rules (Eqs.8-11). Parameters used are σAr = 3.406 A˚ , Ar = 119.4 K and σNe = 2.75 A˚ , Ne = 35.7 K. This leads to an O(N2 ) problem, which increases considerably the execution time of a program for larger systems. For systems with open boundary conditions (like liquid droplets), this method is straightforwardly implemented and reduces to a double sum over all pairs of particles. In the case when periodic boundary conditions are applied, not only the interactions with particles in the central cell but also with all periodic images must be taken into account and formally a lattice sum has to be evaluated U = 1 2 X N i,j=1 X n 0 qiqj |rij − nL| (12) where n is a lattice vector and P n 0 means that for n = 0 it is i 6= j. It is, however, a well known problem that this type of lattice sum is conditionally convergent, i.e. the result depends on the sequence of evaluating the sum (see e.g.27). A method to overcome this limitation was invented by Ewald28 . The idea is to introduce a convergence factor into the sum of Eq.12 which depends on a parameter s, evaluate the sum and put s → 0 in the end. A characterization of the convergence factors was given in Ref.29, 30 . A form which leads to the Ewald sum is an exponential e −sn 2 , transforming Eq.12 into U(s) = 1 2 X N i,j=1 X n 0 qiqj |rij − nL| e −sn 2 (13) The evaluation and manipulation of this equation proceeds now by using the definition for the Γ-function and the Jacobi imaginary transform. A very instructive way of the derivation 219
