1549380323-Statistical Mechanics Theory and Molecular Simulation

250 Isobaric ensembles

×exp (iLǫ, 1 ∆t) exp (iL 1 ∆t)

×exp

(

iL 2

∆t

2

)

exp

(

iLǫ, 2

∆t

2

)

×exp

(

iLNHC−part

∆t

2

)

exp

(

iLNHC−baro

∆t

2

)

+O(∆t^3 ) (5.12.4)

(Tuckermanet al., 2006). In evaluating the action of this propagator, the Suzuki-
Yoshida decomposition developed in eqns. (4.11.16) and (4.11.17) is applied to the op-
erators exp(iLNHC−baro∆t/2) and exp(iLNHC−part∆t/2). The operators exp(iLǫ, 1 ∆t)
and exp(iLǫ, 2 ∆t/2) are simple translation operators. The operators exp(iL 1 ∆t) and
exp(iL 2 ∆t/2) are somewhat more complicated than their microcanonical or canoni-
cal ensemble counterparts due to the barostat coupling and needfurther explication.
The action of the operator exp(iL 1 ∆t) can be determined by solving the first-order
differential equation
r ̇i=vi+vǫri, (5.12.5)

keepingvi=pi/miandvǫ=pǫ/Wconstant with an arbitrary initial conditionri(0)
and then evaluating the solution att= ∆t. Note thatvimust not be confused with
the atomic velocityvi=r ̇i=vi+vǫri.vi=pi/mi, introduced here for notational
convenience to avoid having to writepi/miexplicitly everywhere. Solving eqn. (5.12.5)
yields the finite-difference expression

ri(∆t) =ri(0)evǫ∆t+ ∆tvievǫ∆t/^2

sinh(vǫ∆t/2)

vǫ∆t/ 2

. (5.12.6)

Similarly, the action of exp(iL 2 ∆t/2) can be determined by solving the differential
equation

v ̇i=

Fi

mi

−αvǫvi, (5.12.7)

keepingFiandvǫconstant with an arbitrary initial conditionvi(0) and then evaluating
the solution att= ∆t/2. This yields the evolution

vi(∆t/2) =vi(0)e−αvǫ∆t/^2 +

∆t

2 mi

Fie−αvǫ∆t/^4

sinh(αvǫ∆t/4)

αvǫ∆t/ 4

. (5.12.8)

In practice, the factor sinh(x)/xshould be evaluated by a power series for smallxto
avoid numerical instabilities.^3
Eqns. (5.12.4), (5.12.6) and (5.12.8), together with the Suzuki-Yoshida factorization
of the thermostat operators, completely define an integrator for eqns. (5.9.5). The
integrator can be easily coded using the direct translation technique.

(^3) The power series expansion of sinh(x)/xup to tenth order is
sinh(x)
x
≈

∑^5

n=0

a 2 nx^2 n, (5.12.9)

wherea 0 = 1,a 2 = 1/6,a 4 = 1/120,a 6 = 1/5040,a 8 = 1/362880,a 10 = 1/39916800.