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. 223scheme 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 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 Γ(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 accu￾mulation of the roundoff error will accumulate which will lead in practice to a deviation from the above statement. 223