1549380323-Statistical Mechanics Theory and Molecular Simulation

198 Canonical ensemble

iL 2 =

∑N

i=1

Fi·

∂

∂pi

iLNHC=−

∑N

i=1

pη 1

Q 1

pi·

∂

∂pi

+

∑M

j=1

pηj

Qj

∂

∂ηj

+

M∑− 1

j=1

(

Gj−pηj

pηj+1

Qj+1

)

∂

∂pηj

+GM

∂

∂pηM

. (4.11.5)

Here, the thermostat “forces” are represented as

G 1 =

∑N

i=1

p^2 i

mi

−dNkT

Gj=

p^2 ηj− 1

Qj− 1

−kT. (4.11.6)

Note that the sumiL 1 +iL 2 in eqn. (4.11.4) constitutes a purely Hamiltonian subsys-
tem. The evolution of the full phase space vector

x = (r 1 ,...,rN,η 1 ,...,ηM,p 1 ,...,pN,pη 1 ,...,pηM) (4.11.7)

is given by the usual relation xt= exp(iLt)x 0 As was done in the Hamiltonian case,
we will employ the Trotter theorem to factorize the propagator exp(iL∆t) for a single
time step ∆t. Consider a particular factorization of the form

eiL∆t= eiLNHC∆t/^2 eiL^2 ∆t/^2 eiL^1 ∆teiL^2 ∆t/^2 eiLNHC∆t/^2 +O

(

∆t^3

)

. (4.11.8)

Note that the three operators in the middle are identical to those ineqn. (3.10.22). By
the analysis of Section 3.10, this factorization, on its own, would generate the velocity
Verlet algorithm. However, in eqn. (4.11.8), it is sandwiched betweenthe thermostat
propagators. This type of separation between the Hamiltonian andnon-Hamiltonian
parts of the Liouville operator is both intuitively appealing and, as will be seen below,
allows for easy implementation of both multiple time-scale (RESPA) schemes and
constraints.
The operatoriLNHCcontains many terms, so we still need to break down the
operator exp(iLNHC∆t/2) further. Experience has shown, unfortunately, that a simple
factorization of the operator based on the separate terms iniLNHCis insufficient to
achieve a robust integration scheme. The reason is that the thermostat forces in eqn.
(4.11.6) vary rapidly, thereby limiting the time step. To alleviate this problem, we can
apply the RESPA methodology of Section 3.11 to this part of the propagator. Once
again, experience shows that several hundred RESPA steps are needed to resolve the
thermostat part of the propagator accurately, so RESPA alone cannot easily handle the
rapidly varying thermostat forces. Consider, however, employinga higher-order (than
∆t^3 ) factorization together with RESPA to exp(iLNHC∆t/2). A judiciously chosen