346 Free energy calculations
that the sampling is of equal quality for each simulation. Note that this does not
necessarily mean that all simulations should be of equal length, as the relaxation of
the system along directions in configuration space orthogonal toqmight be longer
or shorter depending on the value ofq. The second assumption is that the biased
histogram in each umbrella windowH ̃k(q) is well estimated by simply applying the
biasing factor directly to the target distributionP(q), i.e.,
H ̃k(q)∝eβ(Ak−A^0 )e−βW(q,s(k))P(q), (8.8.16)which, again, will be approximately true if there is adequate sampling.Once these
assumptions are introduced into eqn. (8.8.15), the coefficients arefinally given by
Ck(q) =nkeβAke−βW(q,s(k))
∑n
j=1nje
βAje−βW(q,s(j)). (8.8.17)Therefore, the distribution becomes
P(q) =∑n
∑ k=1nkPk(q)
n
k=1nke
β(Ak−A 0 )e−βW(q,s(k)). (8.8.18)Although the WHAM procedure might seem straightforward, eqn. (8.8.18) only defines
P(q) implicitly because the free energy factors in eqn. (8.8.18) are directly related to
P(q) by
e−β(Ak−A^0 )=∫
dq P(q)e−βW(q,s(k)). (8.8.19)
Eqns. (8.8.18) and (8.8.19), therefore, constitute a set of self-consistent equations, and
the solution for the coefficients and the free energy factors mustbe iterated to self-
consistency. The iteration is usually started with an initial guess forthe free energy
factorsAk. Note that the WHAM procedure only yieldsAkup to an overall addi-
tive constantA 0. When applying the WHAM procedure, care must be taken that the
assumption of equal quality sampling in each umbrella window is approximately sat-
isfied. If this is not the case, the WHAM iteration can yield unphysicalresults which
might, for example, appear as holes in the final distribution. OnceP(q) is known, the
free energy profile is given by eqn. (8.6.5).
We note, finally, that K ̈astner and Thiel (2005) showed how to combine the um-
brella sampling and thermodynamic integration techniques. Their approach makes use
of the bias in eqn. (8.8.1) as a means of obtaining the free energy derivative dA/dqin
each umbrella window. The method assumes that if the bias is sufficiently strong to
keep the reaction coordinateq=f 1 (r) very close tos(k)in thekth window, then the
probability distributionP(q) can be well represented by a Gaussian distribution:
P(q)≈1
√
2 πσ^2 ke−(q− ̄qk)(^2) / 2 σk 2
. (8.8.20)
Here, ̄qkis the average value of the reaction coordinate in thekth window andσ^2 kis
the variance, both of which are computed via a molecular dynamics orMonte Carlo