Path sampling 595
q(t+ ∆t) =q(t) + ∆tv(t) +A(t)
v(t+ ∆t) =v(t) + ∆tf(q(t))
+
1
2
∆t^2 v(t)f′(q(t)) +σ
√
∆tξ(t)−∆tγv(t)−γA(t), (15.5.15)
where
A(t) =
1
2
∆t^2 (f(q(t))−γv(t)) +σ∆t^3 /^2
(
1
2
ξ(t) +
1
2
√
3
θ(t)
)
(15.5.16)
andξ(t) andθ(t) are Gaussian random variables sampled at timet. The appearance
ofv(t)f′(q(t)) in eqn. (15.5.15), which involves a force derivative, is inconvenient. This
term can be eliminated, however, by making use of eqn. (15.5.12). Substitutings=
t+∆tinto the expression forf(q(s)), we see that ∆tv(t)f′(q(t)) =f(q(t+∆t))−f(q(t)).
Using this fact in eqn. (15.5.15) gives the following solver for the Langevin equation:
q(t+ ∆t) =q(t) + ∆tv(t) +A(t)
v(t+ ∆t) =v(t) +
1
2
∆t[f(q(t+ ∆t)) +f(q(t))]
−∆tγv(t) +σ
√
∆tξ(t)−γA(t). (15.5.17)
Note that if we sett= 0, then the integrator can be recast using the convention of
eqns. (3.10.31) and (3.10.32):
q(∆t) =q(0) + ∆tv(0) +A(0)
v(∆t) =v(0) +
1
2
∆t[f(q(∆t)) +f(q(0))]
−∆tγv(0) +σ
√
∆tξ(0)−γA(0). (15.5.18)
The integrator in eqns. (15.5.17) and (15.5.18) reduce, as expected, to the velocity
Verlet integrator of eqns. (3.8.7) and (3.8.9) whenγ= 0 andσ= 0, which is the limit
of no bath coupling.
As an application of eqns. (15.5.17), we calculate the trajectory, phase space, and
position distribution functions of a harmonic oscillatorW(q) =μω^2 q^2 /2 withμ= 1,
ω= 1, andkT= 1 forγ= 0.5 andγ= 8; the results are shown in Fig. 15.3. The
Langevin equation is integrated for 10^8 steps with a time step of ∆t= 0.01. Fig. 15.3
shows how the trajectory changes betweenγ= 0.5 andγ= 8. Despite the different
values of the damping constant, the computed distribution functions agree with the
analytical distributions. Note that values ofγthat are too small or too large lead to
distortions in the probability distribution. In the present example, values ofγless than
10 −^3 or greater than 100 lead to such distortions.