1549380323-Statistical Mechanics Theory and Molecular Simulation

298 Monte Carlo

this consistency. We can achieve this by accompanying each “bad” molecular dynamics
move by an acceptance step based on the change ∆Hin the Hamiltonian. Since hybrid
Monte Carlo uses the full phase space, the acceptance probabilityis expressed in terms
of the change from one phase space point (r,p) to another (r′,p′), wherepandrare
sets ofNmomenta and coordinates, respectively, and similarly forp′andr′. Thus,
the acceptance probability is given by

A(r′,p′|r,p) = min

[

1 ,e−β{H(r

′,p′)−H(r,p)}

]

= min

[

1 ,e−β∆H

]

. (7.4.1)

Eqn. (7.4.1) ensures that the acceptance is based on the correctcanonical distribution
f(r,p)∝exp[−βH(r,p)].
There is a subtlety in eqn. (7.4.1). As with the other M(RT)^2 schemes we have
considered, eqn. (7.4.1) assumes that the probability distributionT(r′,p′|r,p) for trial
moves satisfies
T(r′,p′|r,p) =T(r,−p|r′,−p′), (7.4.2)

i.e., that the original configurationrcan be reached fromr′by simply reversing the
momenta and running the integrator backwards. In addition, detailed balance requires
that phase space volume be preserved. These conditions will only bemet if the molec-
ular dynamics move is carried out using a symplectic, reversible integration algorithm
such as the velocity Verlet integrator of eqns. (3.8.7) and (3.8.9). IfHis given by

H(r,p) =

∑N

i=1

p^2 i

2 mi

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

then the molecular dynamics move can be expressed in terms of the Trotter-factorized
classical propagator scheme of Section 3.10. IfFi=−∂U/∂riis the force on particle
i, then the Liouville operator is

iL=

∑N

i=1

pi

mi

·

∂

∂ri

+

∑N

i=1

Fi·

∂

∂pi

=iL 1 +iL 2 , (7.4.4)

and the trial molecular dynamics move can be expressed as
(
p′
r′

)

=

[

eiL^2 ∆t/^2 eiL^1 ∆teiL^2 ∆t/^2

]m(p

r

)

, (7.4.5)

which in general, allows the point (r′,p′) to be generated from the initial condition
(r,p) frommiterations of the velocity Verlet integrator using the time step ∆tbefore
eqn. (7.4.1) is applied.
When implementing hybrid Monte Carlo, the values ofm, ∆t, and the target
average acceptance probability need to be optimized. As with the M(RT)^2 algorithms
previously presented, it is difficult to provide guidelines for choosingmand ∆t, as
optimal values are strongly dependent on the system and efficiencyof the computer
program. It is important to note, however, that a single trial movegenerated via eqn.
(7.4.5) requiresmfull force evaluations, in comparison to the relatively inexpensive