Nonequilibrium molecular dynamics 521
iL 1 =
∑N
i=1
[
pi
mi
·
∂
∂ri
+γyi
∂
∂xi
]
iL 2 =
∑N
i=1
Fi·
∂
∂pi
iLmNHC=iLNHC−γ
∑N
i=1
pyi
∂
∂pxi
. (13.5.4)
Here,iLNHCis the Nos ́e–Hoover chain Liouville operator of eqn. (4.11.5). The classical
propagator exp(iL∆t) for a single time step is now factorized in a manner analogous
to eqn. (4.11.8):
eiL∆t= eiLmNHC∆t/^2
[
eiL^2 ∆t/^2 eiL^1 ∆teiL^2 ∆t/^2
]
eiLNHC∆t/^2 +O
(
∆t^3
)
. (13.5.5)
Consider, first, the three operators in the brackets. Ifγ= 0, these operators just
produce a step of velocity Verlet integration in eqns. (3.8.7) and (3.8.9). However, with
γ 6 = 0, the action of the operators is only slightly more complicated. Foreach particle,
the action of the operators can be deduced by solving the three coupled differential
equations
x ̇i=
pxi
mi
+γyi, y ̇i=
pyi
mi
, z ̇i=
pzi
mi
(13.5.6)
subject to initial conditionsxi(0),yi(0),zi(0) with constantpxi,pyi, andpzi. The
solution foryiandziis simply
yi(t) =yi(0) +
pyi
mi
t
zi(t) =zi(0) +
pzi
mi
t. (13.5.7)
Substituting eqn. (13.5.7) foryi(t) into eqn. (13.5.6) forxigives
x ̇i=
pxi
mi
+γ
[
yi(0) +
pyi
mit
]
xi(t) =xi(0) +
[
pxi
mi
+γyi(0)
]
t+
t^2
2
γ
pyi
mi
. (13.5.8)
The evolution forxican be understood by noting that the displacement ofxiin time
is elongated by an amount that depends on the positionyi(0), the shear rate, and
the momentumpyi. It is thisyi-dependent elongation that gives rise to the linear flow
profile. Settingt= ∆tand using the fact that exp(iL 2 ∆t/2) is a simple translation,
we can express the action of the operator in brackets in eqn. (13.5.5) in pseudocode as