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 222where 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 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(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 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 − δ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 calcu￾lation 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