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 i<j i<j where p is the momentum,m the mass of the particles and u and u(3)are pair and three- body interaction potentials.H is an external part,which can include time dependent effects or external sources for a force.All simulated objects are defined within a model description.Often a precise knowledge of the interaction between atoms,molecules or sur- faces are not known and the model is constructed in order to describe the main features of some observables.Besides boundary conditions.which are imposed.it is the model which completely determines the system from the physical point of view.In classical simulations the objects are most often described by point-like centers which interact through pair-or multibody interaction potentials.In that way the highly complex description of electron dynamics is abandoned and an effective picture is adopted where the main features like the hard core of a particle,electric multipoles or internal degrees of freedom of a molecules are modeled by a set of parameters and(most often)analytical functions which depend on the mutual position of particles in the configuration.Since the parameters and functions give a complete information of the system's energy as well as the force acting on each particle through F=-VU,the combination of parameters and functions is also called a force field.Different types of force field were developed during the last ten years.Among them are e.g.MM313,MM414,Dreiding15,SHARP6,VALBON7,UFFI8,CFF9519 AMBER20 CHARMM21,OPLS22 and MMFF23.Typical examples for force field functions are summerized in Fig.2. There are major differences to be noticed for the potential forms.The first distinction is to be made between pair-and multibody potentials.In systems with no constraints,the interaction is most often described by pair potentials,which is simple to implement into a program.In the case where multibody potentials come into play,the counting of interaction partners becomes increasingly more complex and dramatically slows down the execution of the program.Only for the case where interaction partners are known in advance,e.g. in the case of torsional or bending motions of a molecule can the interaction be calculated efficiently by using neighbor lists or by an intelligent way of indexing the molecular sites. A second important difference between interactions is the spatial extent of the potential, classifying it into short and long range interactions.If the potential drops down to zero faster than r-d,where r is the separation between two particles and d the dimension of the problem,it is called short ranged,otherwise it is long ranged.This becomes clear by considering the integral drd oo:n≤d finite n>d (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<j u(ri , rj ) + X N i<j u (3)(ri , rj , rk) + . . . (2) where p is the momentum, m the mass of the particles and u and u (3) are pair and threebody interaction potentials. H1 is an external part, which can include time dependent effects or external sources for a force. All simulated objects are defined within a model description. Often a precise knowledge of the interaction between atoms, molecules or surfaces are not known and the model is constructed in order to describe the main features of some observables. Besides boundary conditions, which are imposed, it is the model which completely determines the system from the physical point of view. In classical simulations the objects are most often described by point-like centers which interact through pair- or multibody interaction potentials. In that way the highly complex description of electron dynamics is abandoned and an effective picture is adopted where the main features like the hard core of a particle, electric multipoles or internal degrees of freedom of a molecules are modeled by a set of parameters and (most often) analytical functions which depend on the mutual position of particles in the configuration. Since the parameters and functions give a complete information of the system’s energy as well as the force acting on each particle through F = −∇U, the combination of parameters and functions is also called a force field. Different types of force field were developed during the last ten years. Among them are e.g. MM313 , MM414 , Dreiding15 , SHARP16 , VALBON17 , UFF18 , CFF9519 , AMBER20 CHARMM21 , OPLS22 and MMFF23 . Typical examples for force field functions are summerized in Fig. 2. There are major differences to be noticed for the potential forms. The first distinction is to be made between pair- and multibody potentials. In systems with no constraints, the interaction is most often described by pair potentials, which is simple to implement into a program. In the case where multibody potentials come into play, the counting of interaction partners becomes increasingly more complex and dramatically slows down the execution of the program. Only for the case where interaction partners are known in advance, e.g. in the case of torsional or bending motions of a molecule can the interaction be calculated efficiently by using neighbor lists or by an intelligent way of indexing the molecular sites. A second important difference between interactionsis the spatial extent of the potential, classifying it into short and long range interactions. If the potential drops down to zero faster than r −d , where r is the separation between two particles and d the dimension of the problem, it is called short ranged, otherwise it is long ranged. This becomes clear by considering the integral I = Z drd r n = ∞ : n ≤ d finite : n > 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) i<j The long-range correction is thereby given as Uire=2πNpo dr r2g(r)u(r) (5) Re where po is the number density of the particles in the system and g(r)=p(r)/po is the radial distribution function.For computational reasons,g(r)is most often only calculated up to Re,so that in practice it is assumed that g(r)=1 forr>Re,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<j u(rij |rij < Rc) + Ulrc (4) The long-range correction is thereby given as Ulrc = 2πNρ0 Z ∞ Rc dr r 2 g(r)u(r) (5) where ρ0 is the number density of the particles in the system and g(r) = ρ(r)/ρ0 is the radial distribution function. For computational reasons, g(r) is most often only calculated up to Rc, so that in practice it is assumed that g(r) = 1 for r > 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