222 Molecular dynamics simulations
powers oftand equating equal powers ofton the left and right hand sides of the
equality [30]. Applying this formula withA=hT ̃andB=hU ̃to increasing orders
of commutators, we find
exp(hJ∇H)=exp(hT ̃)exp(hU ̃)+O(h^2 ) (8.64a)
exp(hJ∇H)=exp(hT ̃/ 2 )exp(hU ̃)exp(hT ̃/ 2 )+O(h^3 ) (8.64b)
etc.,but the extra terms are often tedious to find. AsT ̃andU ̃appear in the exponent,
these expressions do not seem very useful. However, as it follows directly from
Eq. (8.60)that applyingT ̃andU ̃more than once gives zero, we have simply
exp(ahT ̃)= 1 +ahT ̃ (8.65)and similarly for exp(bhU ̃). Therefore, the first order integrator is
p(t+h)=p(t)−h{∂U[x(t)]/∂x}; (8.66a)
x(t+h)=x(t)+h{∂T[p(t+h)]/∂p} (8.66b)which is recognised as the Verlet algorithm (although with a less accurate definition
of the momentum).
The second order integrator is given by
p(t+h/ 2 )=p(t)−h{∂U ̃[x(t)]/∂x}/2; (8.67a)
x(t+h)=x(t)+h{∂T ̃[p(t+h/ 2 )]/∂p}; (8.67b)
p(t+h)=p(t+h/ 2 )−h{∂U ̃[x(t+h)]/∂x}/2. (8.67c)
Applying this algorithm successively, the first and third step can be merged into one,
and we obtain precisely the Verlet algorithm in leap-frog form with a third order
error in the time steph. This error seems puzzling since we know that the Verlet
algorithm gives positions with an error of orderh^4 and momenta with an error of
orderh^2. The solution to this paradox lies in the interpretation of the variablep.
If at timet,p(t)is the continuous time derivative of the continuum solutionx(t),
the above algorithm gives usx(t+h)andp(t+h)both with errorh^3. If however
p(t)is defined as[x(t+h)−x(t−h)]/( 2 h), the algorithm is equivalent to the
velocity-Verlet algorithm and hence gives the positionsx(t+h)with an error of
orderh^4 andp(t+h)is according to its definition given with ah^2 error. The way
in which initial conditions are given defines which case we are in.
Finding higher order algorithms is nontrivial as we do not know the form of
the higher order expansion terms of the operators exp(hT ̃)and exp(hU ̃). However,
Yoshida [21] proposed writing the fourth order integrator in the following form:
S 2 (αh)S 2 (βh)S 2 (αh) (8.68)