1549380323-Statistical Mechanics Theory and Molecular Simulation

258 Isobaric ensembles

“ROLL scalars.” (In the fully flexible cell case, these scalars are replaced by 3×3 matri-
ces.) Note the three operators exp(iLt∆t) exp(iL 2 ∆t/2) exp(iLǫ, 2 ∆t/2) also generate
the following half-step velocities:

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

∆t

2 mi

[

Fi(0) +

∑

k

λkF(c,ik)(0)

]

e−αvǫ∆t/^4

sinh(αvǫ∆t/4)

αvǫ∆t/ 4

≡Rvv(λ,0)v(NHC)i +

∆t

2 mi

RFv(λ,0)

[

Fi(0) +

∑

k

λkF(c,ik)(0)

]

, (5.13.5)

where we have introduce the ROLL scalarsRvv(λ,0) andRFv(λ,0).
The first half of the ROLL algorithm is derived by requiring that the coordinates
in eqn. (5.13.4) satisfy the constraint conditionsσk(r 1 (∆t),...,rN(∆t)) = 0. That is,
eqns. (5.13.4) are inserted into the conditionsσk(r 1 (∆t),...,rN(∆t)) = 0, which are
then solved for the Lagrange multipliersλ. Once the multipliers are determined, they
are substituted into eqns. (5.13.4), (5.13.5), and (5.13.1) to generate final coordinates,
half-step velocities, and the virial contribution to the pressure. Unfortunately, unlike
theNV EandNV Tcases, where the coordinates and velocities depend linearly on the
Lagrange multipliers, the highly nonlinear dependence of eqn. (5.13.4) onλcomplicates
the task of solving for the multipliers. To see how we can solve this problem, we begin
by lettingλ ̃k= (∆t^2 /2)λk. We now seed the ROLL algorithm with a guess{ ̃λ(1)k }for

the multipliers and write the exact multipliers asλ ̃k=λ ̃
(1)
k +δ
̃λ(1)
k. We also assume, at
first, that the ROLL scalars are independent of the multipliers. Thus, when this ansatz
for the multipliers is substituted into eqn. (5.13.4), the coordinatescan be expressed
as

ri(∆t) =r

(1)

i +

1

mi

RFx(λ,0)

∑

k

δλ ̃

(1)

k F

(k)

c,i(0), (5.13.6)

wherer(1)i contains everything except theδλ ̃(1)k -dependent term. Since we are ignoring

the dependence of the ROLL scalars on the multipliers,r(1)i has no dependence on

δ ̃λ(1)k. The constraint conditions now become

σl

(

r(1) 1 +

1

m 1

RFx(λ,0)

∑

k

δ ̃λ(1)k F(c,k 1 )(0),...,

r(1)N +

1

mN

RFx(λ,0)

∑

k

δλ ̃(1)k F(c,Nk)(0)

)

= 0. (5.13.7)

As we did in eqn. (3.9.13), we linearize these conditions using a first-order Taylor
expansion:

σl(r(1) 1 ,...,r(1)N) +

∑N

i=1

∑Nc

k=1

F(c,ik)(1)·

1

mi

RFx(λ,0)δ ̃λ(1)k F(c,ik)(0)≈ 0 , (5.13.8)

whereF(c,ik)(1) =∇iσk(r(1) 1 ,...,r(1)N) are the constraint forces evaluated at the positions

r(1)i. As noted in Section 3.9, we can either solve the full matrix equation ineqn.