of the Ewald sum may be found in Ref.29.30;a heuristic derivation is given Ref.31.For brevity,only the final form of the sum is given here qigjerfc(alrij -nL/L) krie-k2/4a3 2 i.j=1 n rii nL L3 Ureal N (14) erfc(an)+ e-nn2/a2 20 l n2 n≠0 The evaluation of the potential thus splits into four different terms,where the so called self- and surface-terms are constant and may be calculated in the beginning of a simulation. The first two sums,however,depend on the interparticle separations,which need to be evaluated in each time step.It is seen that the lattice sum is essentially split into a sum which is evaluated in real space and a sum over reciprocal space-vectors,k =2mn/L.The parameter o appears formally in the derivation as a result of an integral splitting but it has a very intuitive physical meaning.The first sum gives the potential of a set of point charges which are screened by an opposite charge of the same magnitude but with a Gaussian form factor where the width of the Gaussian is given by a.The second sum subtracts this screening charge,but the sum is evaluated in reciprocal space.Since erfc()=1-erf(z) decays as e-for large z,the first sum contains mainly short range contributions.On the other side,the second sum decays strongly for large k-vectors and thus contains mainly long range contributions.Most often,the parameter a is chosen in way to reduce the evaluation of the real space sum to the central simulation cell.Often,a spherical cutoff is then applied for this term,i.e.contributions of particle pairs,separated in a distance ri>Re are neglected.On the other hand,the reciprocal space sum is conventionally truncated after a maximum wave-vector kmar.All three parameters a,Re,kmaz may be chosen in an optimal way to balance the truncation error in each sum and the number of operations.This balancing even leads to the effect that the Ewald sum may be tuned32.33 to scale with O(N3/2)(for fast methods which have better scaling characteristics,see Ref.31). A detailed analysis of the individual errors occuring in the different sums was given in Ref.34.An alternative derivation of the Ewald sum starts directly by assuming a Gaussian form factor for the screening charge.This gives the opportunity to investigate also form factors,differing from a Gaussian.In these cases the convergence function is in general not known but it is assumed to exist.Different form factors were studied systematically in Ref.35. The present form of the Ewald sum gives an exact representation of the potential energy of point like charges in a system with periodic boundary conditions.Sometimes the charge distribution in a molecule is approximated by a point dipole or higher multipole moments. A more general form of the Ewald sum,taking into account arbitrary point multipoles was given in Ref.3.The case.where also electronic polarizabilities are considered is given in Ref.37 In certain systems,like in molten salts or electrolyte solutions,the interaction between charged species may approximated by a screened Coulomb potential,which has a Yukawa- 220
of the Ewald sum may be found in Ref.29, 30; a heuristic derivation is given Ref.31 . For brevity, only the final form of the sum is given here U = 1 2 ( X N i,j=1 X n 0 qiqj erfc(α|rij − nL|/L) |rij − nL| | {z } Ureal + 4πqiqj L3 X k 1 k 2 e ikrij e −k 2/4α 2 | {z } Ureciprocal + 1 L X n6=0 erfc(αn) |n| + e −π 2n 2 /α2 πn2 − 2α √ π X N i=1 q 2 i | {z } Uself + 4π L3 X N i=1 qi 2 | {z } Usurface ) (14) The evaluation of the potential thus splits into four different terms, where the so called selfand surface-terms are constant and may be calculated in the beginning of a simulation. The first two sums, however, depend on the interparticle separations rij , which need to be evaluated in each time step. It is seen that the lattice sum is essentially split into a sum which is evaluated in real space and a sum over reciprocal space-vectors, k = 2πn/L. The parameter α appears formally in the derivation as a result of an integral splitting but it has a very intuitive physical meaning. The first sum gives the potential of a set of point charges which are screened by an opposite charge of the same magnitude but with a Gaussian form factor where the width of the Gaussian is given by α. The second sum subtracts this screening charge, but the sum is evaluated in reciprocal space. Since erfc(x) = 1 − erf(x) decays as e −x 2 for large x, the first sum contains mainly short range contributions. On the other side, the second sum decays strongly for large k-vectors and thus contains mainly long range contributions. Most often, the parameter α is chosen in way to reduce the evaluation of the real space sum to the central simulation cell. Often, a spherical cutoff is then applied for this term, i.e. contributions of particle pairs, separated in a distance |rij | > Rc are neglected. On the other hand, the reciprocal space sum is conventionally truncated after a maximum wave-vector kmax. All three parameters α, Rc, kmax may be chosen in an optimal way to balance the truncation error in each sum and the number of operations. This balancing even leads to the effect that the Ewald sum may be tuned32, 33 to scale with O(N3/2 ) (for fast methods which have better scaling characteristics, see Ref.31). A detailed analysis of the individual errors occuring in the different sums was given in Ref.34 . An alternative derivation of the Ewald sum starts directly by assuming a Gaussian form factor for the screening charge. This gives the opportunity to investigate also form factors, differing from a Gaussian. In these cases the convergence function is in general not known but it is assumed to exist. Different form factors were studied systematically in Ref.35 . The present form of the Ewald sum gives an exact representation of the potential energy of point like charges in a system with periodic boundary conditions. Sometimes the charge distribution in a molecule is approximated by a point dipole or higher multipole moments. A more general form of the Ewald sum, taking into account arbitrary point multipoles was given in Ref.36 . The case, where also electronic polarizabilities are considered is given in Ref.37 . In certain systems, like in molten salts or electrolyte solutions, the interaction between charged species may approximated by a screened Coulomb potential, which has a Yukawa- 220
like form 1 U= e-klrygl (15) i,j=1 The parameter k is the inverse Debye length,which gives a measure of screening strength in the system.If k<1/L the potential is short ranged and usual cut-off methods may be used.Instead,if 1/L,or generally if u(r =L/2)is larger than the prescribed uncertainties in the energy,the minimum image convention in combination with truncation methods fails and the potential must be treated in a more rigorous way,which was pro- posed in Ref.38,where an extension of the Ewald sum for such Yukawa type potentials was developed. 3 The Integrator For a given potential model which characterizes the physical system,it is the integrator which is responsible for the accuracy of the simulation results.If the integrator would work without any error the simulation would provide exact model results within the errors occuring due to a finite number representation.However,any finite difference integrator is naturally an approximation for a system developing continuously in time.The require- ments for the integrator are therefore to be .accurate,in the sense that it approximates the true trajectory very well (this may be checked with simple model systems for which analytical solutions exist) stable,in the sense that it conserves energy and that small perurbations do not lead to instabilities robust,in the sense that it allows for large time steps in order to propagate the system efficiently through phase space In the following different types of integration schemes are presented.First,simple inte- grators based on Taylor expansions are shown.Later on integrators based on an operator splitting method are discussed which provide the possibility to introduce in an elegant way the integrations of motion on different time scales.Finally,attention is given to more complex situations where molecules with orientational degrees of freedom are considered. 3.1 Expansion Based Integrators The simplest and most straight forward way to construct an integrator is by expanding the positions and velocities in a Taylor series.The class of integrators which may be obtained in that way are called Verlet-Stormer integrators39.40.For a small enough time step 6t the expansion gives rt+60）=r0+v因t+5a倒+b()63+.… (16) vt+0=v0+a0+b0r+ce0r+ (17) 221
like form U = 1 2 X N i,j=1 qiqj e −κ|rij | |rij | (15) The parameter κ is the inverse Debye length, which gives a measure of screening strength in the system. If κ < 1/L the potential is short ranged and usual cut-off methods may be used. Instead, if κ > 1/L, or generally if u(r = L/2) is larger than the prescribed uncertainties in the energy, the minimum image convention in combination with truncation methods fails and the potential must be treated in a more rigorous way, which was proposed in Ref.38 , where an extension of the Ewald sum for such Yukawa type potentials was developed. 3 The Integrator For a given potential model which characterizes the physical system, it is the integrator which is responsible for the accuracy of the simulation results. If the integrator would work without any error the simulation would provide exact model results within the errors occuring due to a finite number representation. However, any finite difference integrator is naturally an approximation for a system developing continuously in time. The requirements for the integrator are therefore to be • accurate, in the sense that it approximates the true trajectory very well (this may be checked with simple model systems for which analytical solutions exist) • stable, in the sense that it conserves energy and that small perurbations do not lead to instabilities • robust, in the sense that it allows for large time steps in order to propagate the system efficiently through phase space In the following different types of integration schemes are presented. First, simple integrators based on Taylor expansions are shown. Later on integrators based on an operator splitting method are discussed which provide the possibility to introduce in an elegant way the integrations of motion on different time scales. Finally, attention is given to more complex situations where molecules with orientational degrees of freedom are considered. 3.1 Expansion Based Integrators The simplest and most straight forward way to construct an integrator is by expanding the positions and velocities in a Taylor series. The class of integrators which may be obtained in that way are called Verlet-Stormer ¨ integrators39, 40 . For a small enough time step δt the expansion gives r(t + δt) = r(t) + v(t) δt + 1 2 a(t)δt 2 + 1 6 b(t)δt 3 + . . . (16) v(t + δt) = v(t) + a(t) δt + 1 2 b(t)δt 2 + 1 6 c(t)δt 3 + . . . (17) 221
where a,b.c are the 2nd,3rd and 4th time derivative of the coordinates.In the same way the expansion may be performed for 6t--6t,which gives rt-0=r因-v因it+a国-b)it±… (18) 6 vt-油=v0-a)t+be)i-ac06t± (19) Adding up Eqs.16,18 and Eqs.17,19 gives for the new positions and velocities r(t+6t)=2r(t)-r(t-6t)+a(t)6t2+O(6t4) (20) v(t+6t)=2v(t)-v(t-6t)+b(t)6t2+O(6t4) (21) A method whose truncation varies as 6t(n+1)is called an n-th order method.Eqs.20,21 are therefore of 3rd order.The drawback of Eq.21 is,however,that it requires the 3rd deriva- tive of the coordinates with respect with to time which is not routinely calculated in MD simulations and thus introduces some additional computational and storage overhead.To overcome this drawback one can simply substract Eq.18 from Eq.16,giving the central difference scheme for the velocity v)=2应rt+6)-rt-)+0(6t) (22) This is,however,one order less in accuracy than Eq.21 and also provides velocities at timestep t,not at t +ot.Since this information i not required by Eq.20 to calculate accu- rately the positions,one may take Eq.22 as an estimate for the velocities from which the kinetic energy of the system is calculated. From the point of view of storage requirements,Egs.20,22 are not optimal,since infor- mation is required form positions not only at time t but also at time t-6t.An equivalent algorithm,which stores only information from one timestep is the so called velocity Verlet algorithm,whic reads rt+0=r)+v96t+2a()62 (23) v(t+8t)=v(t)+t(a()+a(t+8t)) (24) This scheme,however,requires the knowledge of the accelerations,a,at timestept+6t. One may therefore decompose Eq.24 into two steps.First calculate vt+t/2)=v()+2ta (25) then compute the actual forces on the particles at time t+ot and finish the velocity calcu- lation with vt+6t)=vt+t/2)+5a(t+)6t (26) At this point the kinetic energy may be calculated without a time delay of 6t,as it was in Eq.22.Several other schemes have been proposed in literature,such as the leap-frog41 222
where a, b, c are the 2nd, 3rd and 4th time derivative of the coordinates. In the same way the expansion may be performed for δt → −δt, which gives r(t − δt) = r(t) − v(t) δt + 1 2 a(t)δt 2 − 1 6 b(t)δt 3 ± . . . (18) v(t − δt) = v(t) − a(t) δt + 1 2 b(t)δt 2 − 1 6 c(t)δt 3 ± . . . (19) Adding up Eqs.16,18 and Eqs.17,19 gives for the new positions and velocities r(t + δt) = 2r(t) − r(t − δt) + a(t)δt 2 + O(δt 4 ) (20) v(t + δt) = 2v(t) − v(t − δt) + b(t)δt 2 + O(δt 4 ) (21) A method whose truncation varies as δt (n+1) is called an n-th order method. Eqs.20,21 are therefore of 3rd order.The drawback of Eq.21 is, however, that it requires the 3rd derivative of the coordinates with respect with to time which is not routinely calculated in MD simulations and thus introduces some additional computational and storage overhead. To overcome this drawback one can simply substract Eq.18 from Eq.16, giving the central difference scheme for the velocity v(t) = 1 2δt (r(t + δt) − r(t − δt)) + O(δt 3 ) (22) This is, however, one order less in accuracy than Eq.21 and also provides velocities at timestep t, not at t + δt. Since this information i not required by Eq.20 to calculate accurately the positions, one may take Eq.22 as an estimate for the velocities from which the kinetic energy of the system is calculated. From the point of view of storage requirements, Egs.20,22 are not optimal, since information is required form positions not only at time t but also at time t − δt. An equivalent algorithm, which stores only information from one timestep is the so called velocity Verlet algorithm, whic reads r(t + δt) = r(t) + v(t) δt + 1 2 a(t)δt 2 (23) v(t + δt) = v(t) + 1 2 δt(a(t) + a(t + δt)) (24) This scheme, however, requires the knowledge of the accelerations, a, at timestep t + δt. One may therefore decompose Eq.24 into two steps. First calculate v(t + δt/2) = v(t) + 1 2 δta(t) (25) then compute the actual forces on the particles at time t + δt and finish the velocity calculation with v(t + δt) = v(t + δt/2) + 1 2 a(t + δt)δt (26) At this point the kinetic energy may be calculated without a time delay of δt, as it was in Eq.22. Several other schemes have been proposed in literature, such as the leap-frog41 222
scheme or Beeman's42 algorithm.They all have the same accuracy and should produce identical trajectories in coordinate space 3.2 Operator Splitting Methods A more rigorous derivation,which in addition leads to the possibility of splitting the prop- agator of the phase space trajectory into several time scales,is based on the phase space description of a classical system.The time evolution of a point in the 6N dimensional phase space is given by the Liouville equation I(t)=eictr(0) (27) where I (q,p)is the 6N dimensional vector of generalized coordinates. q=q1,...,qN,and momenta,p=p1,...,pN.The Liouville operator,C,is defined as N iC={..,H}=】 (日+ Opj a t∂ai at Op. (28) =1 In order to construct a discrete timestep integrator,the Liouville operator is split into two parts,C=C1+C2,and a Trotter expansion3 is performed eiLot =ei(C1+Ca)ot (29) =eicibtp2eicaoteici8t/2+O(6t3) (30) The partial operators can be chosen to act only on positions and momenta.Assuming usual cartesian coordinates for a system of N free particles,this can be written as N iC1=∑Fp (31) N iC2=∑vm (32) j=1 Applying Eq.29 to the phase space vector I and using the property ea/arf(x)=f(+a) for any function f,where a is independent of x,gives v(t+6t/2)=v(t)+ Fi(t)6t mi 2 (33) r(t+t)=r(t)+v(t+6t/2)6t (34) vt+0）=vt+t/2+Et+的t (35) m12 which is the velocity Verlet algorithm,Eqs.23,25,26. bThis statement is derived from the point of view of accuracy.Since numerical operations are in general not associative a differnt implementation of an algorithm will have different round off errors and therefore the accu- mulation of the roundoff error will accumulate which will lead in practice to a deviation from the above statement. 223
scheme or Beeman’s 42 algorithm. They all have the same accuracy and should produce identical trajectories in coordinate space b 3.2 Operator Splitting Methods A more rigorous derivation, which in addition leads to the possibility of splitting the propagator of the phase space trajectory into several time scales, is based on the phase space description of a classical system. The time evolution of a point in the 6N dimensional phase space is given by the Liouville equation Γ(t) = e iLt Γ(0) (27) where Γ = (q, p) is the 6N dimensional vector of generalized coordinates, q = q1, . . . , qN , and momenta, p = p1, . . . , pN . The Liouville operator, L, is defined as iL = {. . . , H} = X N j=1 ∂qj ∂t ∂ ∂qj + ∂pj ∂t ∂ ∂pj (28) In order to construct a discrete timestep integrator, the Liouville operator is split into two parts, L = L1 + L2, and a Trotter expansion43 is performed e iLδt = e i(L1+L2)δt (29) = e iL1δt/2 e iL2δt e iL1δt/2 + O(δt 3 ) (30) The partial operators can be chosen to act only on positions and momenta. Assuming usual cartesian coordinates for a system of N free particles, this can be written as iL1 = X N j=1 Fj ∂ ∂pj (31) iL2 = X N j=1 vj ∂ ∂rj (32) Applying Eq.29 to the phase space vector Γ and using the property e a∂/∂xf(x) = f(x+a) for any function f, where a is independent of x, gives vi(t + δt/2) = v(t) + Fi(t) mi δt 2 (33) ri(t + δt) = ri(t) + vi(t + δt/2)δt (34) vi(t + δt) = vi(t + δt/2) + Fi(t + δt) mi δt 2 (35) which is the velocity Verlet algorithm, Eqs.23,25,26. bThis statement is derived from the point of view of accuracy. Since numerical operations are in general not associative a differnt implementation of an algorithm will have different round off errors and therefore the accumulation of the roundoff error will accumulate which will lead in practice to a deviation from the above statement. 223
In the same spirit,another algorithm may be derived by simply changing the definitions for C1C2 and C2C1.This gives the so called position Verlet algorithm r(t+t/2)=r:()+v(e)2 (36) vi(t+6t)=v(t)+F:(t+6t/2) (37) mi +)=(+6t/2)+(v()+v() (38) Here the forces are calculated at intermediate positions r;(t+6t/2).The equations of both the velocity Verlet and the position Verlet algorithms have the property of propagating velocities or positions on half time steps.Since both schemes decouple into an applied force term and a free flight term,the three steps are often called half-kick/drift/half kick for the velocity Verlet and correspondingly half-drift/kick/half-drift for the position Verlet algorithm. Both algorithms,the velocity and the position Verlet method,are examples for sym- plectic algorithms,which are characterized by a volume conservation in phase space. This is equivalent to the fact that the Jacobian matrix of a transform 'f(z,p)and p=g(r,p)satisfies （)(-0(）-(-6) (39) Any method which is based on the splitting of the Hamiltonian,is symplectic.This does not yet,however,guarantee that the method is also time reversible,which may be also be considered as a strong requirement for the integrator.This property is guaranteed by symmetric methods,which also provide a better numerical stability4 Methods,which try to enhance the accuracy by taking into account the particles'history(multi-step meth- ods)tend to be incompatible with symplecticness45.46,which makes symplectic schemes attractive from the point of view of data storage requirements.Another strong argument for symplectic schemes is the so called backward error ananlysis47-49.This means that the trajectory produced by a discrete integration scheme,may be expressed as the solution of a perturbed ordinary diffential equation whose rhs can formally be expressed as a power series in ot.It could be shown that the system,described by the ordinary differential equa- tion is Hamiltonian,if the integrator is symplecticso.51.In general,the power series in 6t diverges.However,if the series is truncated,the trajectory will differ only as O(otP)of the trajectory,generated by the symplectic integrator on timescales O(1/6t)52. 3.3 Multiple Time Step Methods It was already mentioned that the rigorous approach of the decomposition of the Liouville operator offers the opportunity for a decomposition of time scales in the system.Supposing that there are different time scales present in the system,e.g.fast intramolecular vibrations and slow domain motions of molecules,then the factorization of Eq.29 may be written in 224
In the same spirit, another algorithm may be derived by simply changing the definitions for L1 → L2 and L2 → L1. This gives the so called position Verlet algorithm ri(t + δt/2) = ri(t) + v(t) δt 2 (36) vi(t + δt) = v(t) + Fi(t + δt/2) mi (37) ri(t + δt) = ri(t + δt/2) + (v(t) + vi(t + δt)) δt 2 (38) Here the forces are calculated at intermediate positions ri(t+δt/2). The equations of both the velocity Verlet and the position Verlet algorithms have the property of propagating velocities or positions on half time steps. Since both schemes decouple into an applied force term and a free flight term, the three steps are often called half-kick/drift/half kick for the velocity Verlet and correspondingly half-drift/kick/half-drift for the position Verlet algorithm. Both algorithms, the velocity and the position Verlet method, are examples for symplectic algorithms, which are characterized by a volume conservation in phase space. This is equivalent to the fact that the Jacobian matrix of a transform x 0 = f(x, p) and p 0 = g(x, p) satisfies fx fp gx gp 0 I −I 0 fx fp gx gp = 0 I −I 0 (39) Any method which is based on the splitting of the Hamiltonian, is symplectic. This does not yet, however, guarantee that the method is also time reversible, which may be also be considered as a strong requirement for the integrator. This property is guaranteed by symmetric methods, which also provide a better numerical stability44 . Methods, which try to enhance the accuracy by taking into account the particles’ history (multi-step methods) tend to be incompatible with symplecticness45, 46 , which makes symplectic schemes attractive from the point of view of data storage requirements. Another strong argument for symplectic schemes is the so called backward error ananlysis47–49 . This means that the trajectory produced by a discrete integration scheme, may be expressed as the solution of a perturbed ordinary diffential equation whose rhs can formally be expressed as a power series in δt. It could be shown that the system, described by the ordinary differential equation is Hamiltonian, if the integrator is symplectic50, 51 . In general, the power series in δt diverges. However, if the series is truncated, the trajectory will differ only as O(δt p ) of the trajectory, generated by the symplectic integrator on timescales O(1/δt) 52 . 3.3 Multiple Time Step Methods It was already mentioned that the rigorous approach of the decomposition of the Liouville operator offersthe opportunity for a decomposition of time scales in the system. Supposing that there are different time scales present in the system, e.g. fast intramolecular vibrations and slow domain motions of molecules, then the factorization of Eq.29 may be written in 224