1549380323-Statistical Mechanics Theory and Molecular Simulation

344 Free energy calculations

the equilibrium value of the harmonic restraining potential is chosen to be a set of
intermediate pointss(1),...,s(n)betweens(i)ands(f), so that the reaction coordinate
q =f 1 (r) is “driven” between the two endpoint values. Each molecular dynamics
or Monte Carlo simulation will then yield a biased probability distributionP ̃(s,s(k))
k= 1,...,nwiths(1) =s(i) ands(n)=s(f)about each points(k). From this set
of distribution functions, the true free energy profile across theentire range of the
reaction coordinate must be reconstructed.
When a biasing potential is used, reconstructing the full distribution requires an
unbiasing procedure. This parallels the use of a constraint in the bluemoon ensemble
approach, which also requires an unbiasing factor in ensemble averages. The technique
we will present in this section is known as theweighted histogram analysis method,
or WHAM (Kumaret al., 1992). WHAM begins by defining the biased probability
distribution generated from each molecular dynamics or Monte Carlosimulation

P ̃(q,s(k)) = eβAk

∫

dNre−βU(r)e−βW(f^1 (r),s

(k))

δ(f 1 (r)−q), (8.8.2)

whereAkis (apart from additive constants) the free energy associated with the biased
potential

e−βAk=

∫

dNre−βU(r)e−βW(f^1 (r),s

(k))

= e−βA^0

〈

e−βW(f^1 (r),s

(k))〉

. (8.8.3)

In eqn. (8.8.3), the average is taken with respect to the unbiased potentialU(r), and

e−βA^0 =

∫

dNre−βU(r) (8.8.4)

is just the unbiased configurational partition function. The factor exp(βAk) is, there-
fore, the correct normalization constant for eqn. (8.8.2).
Next, we define the unbiased probability distribution functionPk(q) corresponding
to the distributionP ̃(q,s(k)) as

Pk(q) = e−β(Ak−A^0 )eβW(q,s

(k)) ̃

P(q,s(k)). (8.8.5)

We now “glue” these distributions together to give the full probability distribution
functionP(q) by expressingP(q) as a linear combination of the distributionsPk(q):

P(q) =

∑n

k=1

Ck(q)Pk(q)

=

∑n

k=1

Ck(q)

[

e−β(Ak−A^0 )eβW(q,s

(k)) ̃

P(q,s(k))

]

, (8.8.6)

where{Ck(q)}is a set of coefficients that must be optimized to give the best repre-
sentation of the true distributionP(q). The coefficients must satisfy the constraint

∑n

k=1

Ck(q) = 1. (8.8.7)

In order to determine the coefficients, we seek to minimize the statistical error
in the distribution generated by the WHAM procedure. LetH ̃k(q) be the (biased)