1549380323-Statistical Mechanics Theory and Molecular Simulation

Time evolution operator 113

In the same way, the third operator, which involves another derivative with respect to
momentum, yields:

exp

(

∆tF(x)

∂

∂p

)





x+ ∆tmp

p+∆ 2 tF

(

x+ ∆tmp

)





=





x+∆mt

(

p+∆ 2 tF(x)

)

p+∆ 2 tF(x) +∆ 2 tF

(

x+∆mt

(

p+∆ 2 tF(x)

))



. (3.10.30)

Using the fact thatv=p/m, the final positionx(∆t) can be written as (replacing the
“(0)” label on the initial conditions):

x(∆t) =x(0) + ∆tv(0) +

∆t^2

2 m

F(x(0)). (3.10.31)

Eqn. (3.10.31) is equivalent to a second-order Taylor expansion ofx(∆t) up to second
order in ∆tand is also the position update part of the velocity Verlet algorithm. From
eqn. (3.10.31), the momentum update step can be written compactly (usingv=p/m
and replacing the “(0)” label) as

v(∆t) =v(0) +

∆t

2 m

[F(x(0)) +F(x(∆t))], (3.10.32)

which is the velocity update part of the velocity Verlet algorithm.
The above analysis demonstrates how we can obtain the velocity Verlet algorithm
via the powerful formalism provided by the Trotter factorization scheme. Moreover,
it is now manifestly clear that the velocity Verlet algorithm constitutes a symplectic,
unitary, time-reversible solver that preserves the important symmetries of classical
mechanics (see also Problem 1.5 in Chapter 1). Since eqns. (3.10.31) and (3.10.32)
together constitute a symplectic algorithm, each term of the product implied by eqn.
(3.10.20) is symplectic. Finally, in the limit ∆t→0 andP → ∞, exact classical
mechanics is recovered, which allows us to conclude that time evolution under Hamil-
ton’s equations is, indeed, symplectic, as claimed in Section 1.6. While it may seem
as though we arrived at eqns. (3.10.31) and (3.10.32) using an overlycomplicated
formalism, the power of the Liouville operator approach will be readilyapparent as
we encounter increasingly complex numerical integration problems inthe upcoming
chapters.
Extending the analysis toN-particle systems inddimensions is straightforward.
For the standard Hamiltonian in Cartesian coordinates

H=

∑N

i=1

p^2 i

2 mi

+U(r 1 ,...,rN), (3.10.33)

the Liouville operator is given by