1549380323-Statistical Mechanics Theory and Molecular Simulation

Adiabatic free energy dynamics 323

Thus, if the region betweenλ= 0 andλ= 1 is transformed into a low-probability
region,λwill spend a very small fraction of its time there in a molecular dynamics
run, and sampling near the endpointsλ= 0 andλ= 1 will be enhanced.
To illustrate how we can achieve a barrier in the free energy profileA(λ), let us
examine a simple example. Consider taking two uncoupled harmonic oscillators and
usingλ-switching to “grow in” a bilinear coupling between them. Ifxandyrepresent
the coordinates of the two oscillators with massesmxandmyand frequenciesωxand
ωy, respectively, then the two potential energy functionsUAandUBtake the form

UA(x,y) =

1

2

mxωx^2 x^2 +

1

2

myω^2 yy^2

UB(x,y) =

1

2

mxωx^2 x^2 +

1

2

myω^2 yy^2 +κxy, (8.3.5)

whereκdetermines the strength of the coupling between the oscillators. For this
problem, the integration overxandycan be performed analytically (see Problem 8.1),
leading to the exact probability distribution function inλ

P(λ) =

C

√

mxω^2 xmyω^2 y(f(λ) +g(λ))^2 −κ^2 g^2 (λ)

, (8.3.6)

whereCis a constant, from which the free energy profileA(λ) =−kTlnP(λ) becomes

A(λ) =

kT

2

ln

[

mxωx^2 myω^2 y(f(λ) +g(λ))^2 −κ^2 g^2 (λ)

]

. (8.3.7)

The reader can easily verify that the free energy difference ∆A=A(1)−A(0) does
not depend on the choice off(λ) andg(λ). For concreteness, let us set the parameters
mx=my= 1,ωx= 1,ωy= 2,kT= 1, andκ= 1. First, consider switches of the form
f(λ) = (λ^2 −1)^2 andg(λ) = ((λ−1)^2 −1)^2. The solid line in Fig. 8.2(a) shows the
free energy profile obtained from eqn. (8.3.7). The free energy profile clearly contains
a barrier betweenλ= 0 andλ= 1. If, on the other hand, we choosef(λ) = (λ^2 −1)^4
andg(λ) = ((λ−1)^2 −1)^4 , the free energy profile appears as the dashed line in
Fig. 8.2(a), and we see that the profile exhibits a deep well. A well indicates a region
of high probability and suggests that in a molecular dynamics calculation,λwill spend
considerably more time in this irrelevant region than it will near the endpoints. In this
case, therefore, the quartic switches are preferable. For comparison, these two choices
forf(λ) (g(λ) is just the mirror image off(λ)) are shown in Fig. 8.2(b). It can be seen
that small differences in the shape of the switches lead to considerable differences in
the free energy profiles.
Suppose we now try to use a molecular dynamics calculation based on the Hamil-
tonian in eqn. (8.3.1) to generate a free energy profileA(λ) with a substantial barrier
betweenλ= 0 andλ= 1. We immediately encounter a problem: The probability for
λto cross this barrier becomes exponentially small! So have we actuallyaccomplished
anything by introducing the barrier? After all, what good is enhancing the sampling
in the endpoint regions if the barrier between them cannot be easily crossed? It seems
that we have simply traded the problem of inefficient sampling at the endpoints of the