10.4 Histogram Reweighting

To further generalize our discussion of free-energy methods, it is convenient to introduce the concept of the “order parameter” \(z\) which is computed by a mapping function \(\theta (\rb ^N)\), which is generally just any function of configuration (for example, \(\Delta \mathscr {U} = \mathscr {U}_1-\mathscr {U}_0\)). For a system obeying the potential \(\mathscr {U}_0\), we can define the order parameter probability density as \begin{equation} p_0 = \frac {\ds \int \rb ^N e^{-\beta \mathscr {U}_0}\delta \left [z-\theta (\rb ^N)\right ]}{\ds Z_0} \end{equation} where we are using the shorthand \begin{equation} Z_0 = \int \rb ^N e^{-\beta \mathscr {U}_0} \end{equation} to represent (now the unscaled) configurational integral of the Boltzmann factor. The “Landau” free energy is then expressed as \begin{equation} \label {eq:Landau} F_{\rm Landau}(z)=-k_BT\ln p_0(z) \end{equation} We’d generally like to be able to compute \(p_0\) or, alternatively, \(F_{\rm Landau}\), but as we have already seen, it is often impossible to perform adequate sampling over an entire range of order parameter. Suppose we restrict sampling to a region in order parameter space using a “bias potential” \(W_i(\theta (\rb ^N))\). Under the action of this bias, we can compute directly a biased probability density: \begin{equation} \label {eq:hw:biasedp} p_i(z) = \frac {\ds \int \rb ^N e^{-\beta [\mathscr {U}_0+W_i(z)]}\delta \left [z-\theta (\rb ^N)\right ]}{\ds Z_i} \end{equation} where we use the shorthand \begin{equation} Z_i = \int \rb ^N e^{-\beta [\mathscr {U}_0+W_i(\theta (\rb ^N))]} \end{equation} Note that in Eq. 269, we have made use of the fact that the Dirac delta function only permits \(\theta (\rb ^N)\) to equal \(z\), so we do not need to express the argument of \(W_i\) as \(\theta (\rb ^N)\) in the numerator. This is important, because it allows us to express the unbiased probability density: \begin{equation} p_0(z) = \exp (+\beta W_i(z))\frac {Z_i}{Z_0}p_i(z) \end{equation}

Let’s now consider that we compute \(p_i\) by histogramming into bins of constant width \(\Delta z\). Let’s call the histogram \(H_i\) and the total number of hits in the histogram \(M_i\). Then one way to write the relation between the approximate probability density and the histogram is \begin{equation} p_i\Delta z = \frac {H_i}{M_i} \end{equation}

Now we imagine we can use several different \(W_i\)’s together, each of which focus sampling on a particular region of order-parameter space, and we propose that the unbiased probability density can be constructed by a superposition of the biased probabilities: \begin{equation}\label {eq:hw:super} \begin {aligned} p_0(z)& =\sum _{i=1}^n a_i(z)\exp (+\beta W_i(z))\frac {Z_i}{Z_0}p_i(z) \\ \Rightarrow p_0(z)\Delta z& = \sum _{i=1}^n a_i(z)\exp (+\beta W_i(z))\frac {Z_i}{Z_0}\frac {H_i(z)}{M_i} \end {aligned} \end{equation} subject to \begin{equation} \sum _{i=1}^n a_i(z) = 1 \ \ \mbox {for all $z$} \end{equation} In Eq. 273, we have supposed that \(p_0\Delta z\) is a normalized histogram with the same bin size \(\Delta z\) as all \(H_i\).

Via minimization of the squared fluctuations of \(p_0(z)\) (since it will fluctuate in a simulation), one can arrive at an expression for \(a_i\) (math not shown):

where the normalization condition is met by \begin{equation} \alpha (z) = \frac {1}{\ds \sum _{i=1}^n \exp [-\beta W_i(z)]M_i\frac {Z_0}{Z_i}} \end{equation}

This gives a reweighted histogram of \begin{equation} p_0(z)\Delta z = \frac {\ds \sum _{i=1}^n H_i}{\ds \sum _{i=1}^n \exp [-\beta W_i(z)]M_i\frac {Z_0}{Z_i}} \end{equation}

Now, let’s define \(F_i\) via \begin{equation}\label {eq:wham:F_cont} \begin {aligned} \exp \beta F_i & = -\frac {Z_i}{Z_0} \\ & = \frac {\ds \int d\rb ^N e^{-\beta (\mathscr {U}_0+W_i)}}{\ds Z_0} \\ & = \frac {\ds \int dz \int \rb ^N e^{-\beta (\mathscr {U}_0+W_i)}\delta (\theta (\rb ^N)-z)}{\ds Z_0} \\ & = \int dz\ e^{-\beta W_i(z)}\frac {\ds \int \rb ^N e^{-\beta \mathscr {U}_0}\delta \left [z-\theta (\rb ^N)\right ]}{Z_0} \\ & = \int dz\ p_0(z) e^{-\beta W_i(z)} \\ \Rightarrow F_i & = -\frac {1}{\beta }\ln \int dz\ p_0(z) e^{-\beta W_i(z)} \end {aligned} \end{equation}

This allows us to write the unbiased histogram as \begin{equation} \label {eq:wham:p0} p_0(z)\Delta z = \frac {\ds \sum _{i=1}^n H_i}{\ds \sum _{j=1}^n \exp [-\beta (W_j(z)-F_j)]M_j} \end{equation} Given that we have a finite bin width, the integration in Eq. 278 can be approximated as \begin{equation} \label {eq:wham:F} F_j = -\frac {1}{\beta }\ln \sum _{k=0}^{n_{\rm bins}-1} p_0(z_k)\exp [-\beta W_j(z_k)]\Delta z \end{equation} where \(z_k = z_{\rm min}+(k+\frac {1}{2})\Delta z\), with \(z_{\rm min}\) the lower bound of the histogram domain and \(n_{\rm bins}\) is the number of bins.

Eq. 279 and 280 are called (by some) the WHAM equations (for Weighted Histogram Analysis Method [45]). They can be solved iteratively to yield an unbiased probability provided with a set of biased histograms \(H_i\) obtained from running MC simulations on \(\mathscr {U}_0+W_i\). The standard WHAM approach is

• Set \(F_j\) = 0 for all \(j\).
• Compute the first estimate of \(p_0\) using Eq. 279 from the \(H_i\)’s and \(M_i\)’s.
• Compute a new estimate for the \(F_j\)’s using Eq. 280.
• Repeat the last two steps until the \(F_j\)’s stop changing within some tolerance.

As an example, consider a single degree of freedom \(x\) (our “particle”) moving under Brownian dynamics (BD) on a quartic two-well potential: \begin{equation}\label {eq:bd:v} \begin {aligned} \dot {x} = -\frac {1}{\gamma }\frac {dV}{dx} + \sqrt {2k_BT/\gamma }\eta ,\ \ \mbox {where} \\ V(x) = a(x^2-1)^2 + bx^2 + cx \end {aligned} \end{equation} Here, \(\gamma \) is the BD friction, \(T\) is the temperature, and \(\eta \) is a random variate centered on zero with unit variance. A symmetrical two-well potential with a barrier at \(x=0\) of about 13 can be had by using \(a\) = 0.02, \(b\) = -1, and \(c\) = 0. Since there is only one degree of freedom, it can easily be shown that \(F(x) = V(x)\). So, if we run BD, sampling \(x\) and then histogramming it into \(p_0\), then we should see that \(F=-k_BT\ln p_0\) and \(V\) agree. Computing \(F\) from a directly histogrammed variable is sometimes called “Boltzmann inversion”.

Fig. 42 shows the calculation of \(F\) using this approach, with BD on \(V\) using the code bd-w.c. We show that sampling with a reduced temperature of 10 requires about 2 billion BD timesteps to reproduce \(V\). If we reduce the temperature to 1, the large barrier between the two wells is not surmounted and \(V\) is therefore not correctly reproduced.

Fig. 43 shows that \(V\) can be reproduced at \(T\) of 1.0 using WHAM using 21 windows and harmonic bias potentials \(W_i\) spaced uniformly along \(z\), each with a \(k\) of 20. That is, each bias potential is of the form \begin{equation} W_i(z) = \frac {1}{2}k(z-z_{0,i})^2 \end{equation} where \(z_{0,i}\)’s are uniformly spaced between -10 and 10. (Actually, -10 and 10 are the domain boundaries; this domain of size 20 is divided into 21 windows, and each \(z_0\) is in the center of its window. So, \(z_00 = -10+0.5(20/21) = -9.5238\), etc. An odd number of windows is a good idea so that we can better resolve the tippy top of the barrier.)

Here, only two million BD steps were run per window, for a total of 42 million steps. This means this WHAM dataset required more than twenty times less BD sampling at a temperature of 1.0 than the brute force sampling approach could at a temperature of 10.0.

As a second example of WHAM, consider MD simulation of butane molecule in vacuum at \(T\) = 310 K. One order parameter we may consider here is the C\(_1\)-C\(_4\) distance. When the molecule is in trans, the C\(_1\)-C\(_4\) distance is about 4.5 Å, while when it is in gauche, the C\(_1\)-C\(_4\) distance is about 3.2 Å. We expect there to be a small free-enegy barrier between these two states when characterized using the C\(_1\)-C\(_4\) distance. Fig. 44 shows results of both WHAM and long MD simulations, showing indeed that the trans and gauche states are separated by a small barrier of about 2 kcal/mol at 310 K, easily surmountable in MD. WHAM reconstruction of \(F\) from 20 biased histograms generated from MD restrained using harmonic window potentials with spring constants of 100 kcal/mol-Å\(^2\) shows perfect agreement with the Boltzmann-inverted result.

Butane is an almost trivially simple case; the C1-C4 distance is relatively easily sampled in standard MD at 310 K in few million time-steps. This is of course evident because we can generate a Boltzmann-inverted free energy directly from the histogram of this distance generated by MD. But it is important to note that the order parameter here, the C1-C4 distance, has a small amount of degeneracy: every value of this parameter except for its minimum and maximum has two major realizations determined by positive and negative senses of the dihedral angle. (There are of course some minor realizations due to angle stretching.) That means that each window’s simulation should sample both the positive and negative regions of the order parameter space, and we must take care that the window potential is not so strong that these dihedral angle fluctuations cannot occur; if so, we would undersample the gauche states for sure. Fig. 45 shows traces of the C1-C4 distance for each window and corresponding traces of the dihedral angle; clearly for this value of \(k\), each window is indeed able to sample both positive and negative values of the dihedral angle.