Our objective is to construct the underlying unbiased free energy G(x) that is most consistent with the observed simulation data. For this purpose we aim to determine the unbiased equilibrium probability distribution {p_l} (l =1, …, M), representing the probability of finding the coordinate x in each bin when no biasing potential is applied. Then G(x) is given by where k_B is the Boltzmann constant and Δ_l is the width of bin l. Given the {p_l}, the expected probability distribution {p_il} in simulation i is also known:⁶ where c_il is determined by the biasing potential at the center of each bin: and f_i is a normalization factor to ensure that {p_il} sum to one: f_i is thus the reciprocal of the partition function of simulation i, Given the probabilities {p_il} for simulation i, the likelihood of having {n_il} counts in the respective bins obeys the multinomial distribution:⁶ The overall likelihood for jointly observing the counts in the S independent simulations is then given by Substituting Eqs. (5) and (9) into Eq. (10) and then taking the logarithm, we obtain⁶ in which the negative log-likelihood A includes all the terms containing the {p_l}: where M_l is the total count in the l -th bin, summed over all simulations: Finding the maximum likelihood estimates of {p_l} is then equivalent to the minimization of A. We note that general discussions of the maximum likelihood approach applied to biased sampling can be found in the statistical literature.^13–15

To obtain the minimum of A, we take the derivative of A with respect to each individual p_l, noting that in Eq. (12) f_i is a function of the {p_l} through Eq. (7): By demanding these derivatives to be zero, we obtain for each bin. Eqs. (7) and (15) jointly form the set of coupled nonlinear WHAM equations, which are traditionally solved by iteratively substituting {p_l} and {f_i} to obtain a new set of values until self-consistency is achieved. This direct iteration method is numerically inefficient. In the following we introduce techniques with superior convergence rate that directly aim to minimize A, as given in Eq. (12). In addition, for a given problem, A can be used to compare the accuracy of different results, with those closer to the exact solution having smaller A values.

To reduce the dimensionality of the optimization problem, we rewrite the minimization of A with respect to the {p_l} into an equivalent minimization with respect to the {f_i}. Since there are typically far fewer simulations than bins, S < M, this reformulation reduces the computational cost. Substitution of Eq. (15) into Eq. (12) results in a new function of the {f_i} alone: By taking the gradient with respect to the {f_i}, it is then straightforward to show that the minimum of this new function coincides with the solution of the WHAM equations: where we used Eq. (15) to rewrite the gradient in terms of the {p_l}. At the minimum of these derivatives are zero, and we thus recover Eq. (7). Solving the WHAM equations is equivalent to minimizing Eq. (16) with respect to the {f_i}, which is numerically more efficient than minimizing Eq. (12) with respect to the {p_l} in typical cases.

One may also work on the logarithm of {f_i} by a change of variables: According to Eq. (8), the {g_i} actually correspond to the total free energies of the system in the presence of the biasing potential i. Using the {g_i} as the free variables eliminates the possibility of having (unphysical) negative f_i values in the search of the minimum. The optimization function becomes: Furthermore, it can be shown¹³ that is a convex function with nonnegative second derivatives everywhere, and thus has a single minimum.¹³ The derivatives of with respect to the {g_i} are given by Because multiplying all {f_i} by a constant factor, or adding a constant to all {g_i}, does not alter the free energy profile G(x), we can set g₁ to zero without loss of generality, and minimize with respect to g₂, …, g_S.

In typical umbrella sampling setups, the histogram of a particular umbrella window only substantially overlaps with the histograms of neighboring windows. Consequently, the derivative for window i is predominantly determined by those values of {g_j} with |j − i| being small, and is hardly affected by distant windows. The derivative therefore mainly reflects the local consistency of the free energy across a small range of umbrella windows, but poorly indicates the real distance of g_i to its true value in the exact solution. To improve the performance of the optimization, one may use the incremental changes of the {g_i} over adjacent windows as free variables in the minimization: The function to be minimized and its derivatives with respect to the {Δg_i} (i=1, …, S −1) are then given by Our experimentation shows that the change of variables to {Δg_i} makes the convergence considerably faster, especially when the number of umbrella windows is large.

To estimate the statistical errors in the {f_i}, we first define the total free energy of the system in the presence of a harmonic biasing potential K(x − r)² / 2 centered at r: Using Eq. (8), we have with g_i defined in Eq. (18). Therefore the {f_i} or {g_i} are determined by the values of G_r(r) at the discrete points {r_i}.

Furthermore, taking the derivative of G_r(r), we have: This derivatives at the {r_i} are thus given by where F_i = K(r_i − x) is the restraining force arising from the harmonic potential w_i(x) in Eq. (3), and _i denotes the ensemble average in the presence of the biasing potential w_i(x). The i -th umbrella simulation actually represents such an ensemble, and thus can be used to obtain the ensemble averages: where and are the mean position of x and the mean restraining force in the simulation, respectively. If the biasing potentials are not harmonic, the estimators of the free energy gradient have to be modified appropriately (see Appendix B). We can integrate the gradients at the {r_i} using the trapezoidal rule:

Eqs. (24) and (28) show that the {f_i} or {g_i} can be obtained approximately by integrating the gradients of G_r(r), or the mean restraining forces from the simulations,¹⁸ analogous to thermodynamic integration. This method can also be conveniently applied in multidimensional cases, and therefore has been adopted in conjunction with the string method¹⁹ to obtain a free energy profile in a high-dimensional coordinate space. For the one-dimensional PMFs considered in this study, the method is expected to give similar results as WHAM does, if G_r(r) varies smoothly as a function of r between windows. Whereas the systematic errors arising from discretization of the chosen coordinate are expected to be smaller in WHAM than in the gradient-based method, the statistical errors in the latter can nonetheless be taken as reasonable estimates of the statistical errors in the WHAM result. We expect the quality of this approximation to decrease when the windows are spaced too far compared to features in G_r(r) .

Because the {f_i} or {g_i} can be expressed as functions of the , their statistical errors can also be calculated from the variances of the , which according to Eq. (27) are determined by the variances of the mean coordinates : If the {r_i} are equally spaced, i.e., Eq. (28) can be simplified into The resulting variance of the estimator for the free energy difference over repeated measurements is then given simply by The corresponding standard deviation, given by the square root of the variance, is typically taken as the statistical error of the free energy difference above. Here we assumed that the errors in the are uncorrelated. For some simulation protocols, the of different windows are correlated and one should then include their covariances in addition to the variances. In Hamiltonian replica exchange simulations²⁰ the {x_i} of different windows are statistically independent of each other in the joint probability distribution. For sufficiently long Hamiltonian replica exchange simulations that fully explore the combined phase space, the covariances of the are thus expected to vanish, and Eq. (32) approximately holds in such cases.

If the harmonic springs are relatively stiff (with large K), as is often the case in umbrella sampling simulations, p(x) would vary little within a narrow umbrella window. In such cases we have G_r(r) in Eq. (23) can then be simplified into Therefore under the stiff spring approximation,²¹ the {G_r(r_i)} directly represent the PMF, and thus also share the same statistical errors as the PMF.

The variance var can be estimated from the trajectory of x in simulation i. For a correlated trajectory, this variance depends on the autocorrelation function of the data. By block averaging, however, var can be estimated without explicitly computing the autocorrelation function.¹² In this method the time series x_iα of the x - coordinate in window i is divided into n blocks of sufficiently large size m, so that the averages of the n blocks are independent of each other, and then var is estimated from the variance of these averages:¹² For small n, this estimate tends to be noisy; for large n, the blocks are correlated. While n between 5 and 10 are typically used in practice, for accurate estimates it is advisable to plot Eq. (35) as a function of log(n) to identify a plateau¹² at small n and use that for the estimate of var.

From the variances of in each individual simulation, the uncertainty in the {g_i} (assuming g₁=0) can be obtained via Eqs. (24) and (32): This estimate is the central result of the second part of this paper, and clearly shows the accumulation of the statistical error over the intermediate umbrella windows. Under the stiff spring approximation discussed earlier, G_r(r_i) or {g_i} directly represent the PMF (Eq. 34), and therefore Eq. (36) also gives the error estimate for the PMF, thus leading to Eq. (1) except for minor differences arising from the first and last windows.

Consider two simulations in adjacent umbrella windows centered at r₁ and r₂. From the simulation trajectories, we obtain the probability distributions p₁(x) and p₂(x), and aim to test whether they are consistent with a common underlying unbiased probability distribution p(x). In principle, we may reconstruct p(x) from either p₁(x) or p₂(x), and then compare the two resulting distributions for consistency. In practice, however, the p(x) obtained from p_i (x) is only accurate in the vicinity of r_i, and will bear increasingly large statistical uncertainties at increasing distances from r_i. Indeed, in the WHAM equations p₁(x) and p₂(x) primarily contribute to the shape of p(x) around the region between r₁ and r₂, and the comparison should therefore also be focused on this range. For this purpose, we imagine a virtual umbrella window at the midpoint r* =(r₁ + r₂) / 2, with the biasing potential given by as in Eq. (3). The probability distribution of x in this virtual umbrella window can be estimated from either p₁(x) or p₂ (x): If p₁(x) and p₂ (x) indeed arise from a common unbiased probability distribution, (x) and (x) should be identical within statistical uncertainties. Furthermore, (x) and (x) are peaked in the region where p₁(x) and p₂ (x) have maximum overlap, and could be compared on the same footing. A significant disagreement of (x) and (x) would then indicate that the two simulations are probably sampling different free energy surfaces.

Several statistical methods can be applied to test (x) and (x). In this study we adopt a simple protocol based on the Kolmogorov-Smirnov test. The method uses the maximum difference between the two respective cumulative distribution functions to quantify the deviation between (x) and (x). For two histograms arising from a common probability distribution, the expected value of this maximum difference asymptotically approaches zero following the inverse square root of the sample size, as the latter approaches infinity. Therefore we define an inconsistency coefficient, θ_1,2, as an empirical measure for the discrepancy between (x) and (x): Note that N₁ and N₂ above should be the effective number of independent samples, which can be estimated from Eq. (37) for correlated data. An abnormally high θ_1,2 indicates inconsistency between the results of the two simulations. For any pair of adjacent umbrella windows i and i +1, we can similarly calculate an inconsistency coefficient θ_i,i+1. We note that other methods, such as Pearson's χ² test, can also be used to examine (x) and (x).

We first compare the performance of the superlinear (trust region and BFGS) algorithms and the traditional direct iteration method in solving the WHAM equations. In all methods the iterations start with given initial values of the probabilities {p_l} or the normalization factors {f_i}. Although it is common practice to assign a constant initial value to all {p_l} or {f_i}, it was shown that using more accurate estimates of the free energy as initial values can significantly speed up the convergence in the direct iteration method.²³ In fact, approximate {f_i} or {g_i} can be calculated by integrating the gradients, or the mean restraining forces (Eqs. 24, 31). Here we test each method, first for constant initial values, and then by using the coarse-grained profile determined from the mean forces. In our implementation of the direct iteration method, the convergence is deemed achieved when the relative change in {p_l} for any l by a WHAM iteration is smaller than a given threshold δ. We test each case with a larger (10⁻³) and a smaller (10⁻⁶) δ value.

As shown in Table 1, the trust region and BFGS numerical algorithms yield more accurate results than the traditional direct iteration method, indicated by smaller values of the target function A. In fact, the direct iteration method with a convergence threshold of δ = 10⁻³ leads to a rather inaccurate free energy profile (Fig. 1A, blue dashed curve), with an error of more than 2 k_BT, when starting with uniform initial values. When the gradient-based estimates of the {f_i} are used as initial values, the direct iteration method indeed converges faster, with the threshold of δ = 10⁻³ almost satisfied already at the start of the iterations. To achieve better accuracy with a more stringent threshold of δ = 10⁻⁶, however, a large number of iterations is still required, taking a computational time at least an order of magnitude longer than the trust region and BFGS methods. The example also demonstrates that apparently small variations between WHAM iterations, on the order of δ, do not necessarily mean that convergence at that level of accuracy has been achieved. At the end of the WHAM iterations in our test, the free energies change by less than 10⁻⁶ k_BT in a single iteration, although their absolute errors are still on the order of 0.01 k_BT, implying that to gain one more significant digit in the result, one needs to perform thousands of iterations. In contrast, the convergence rate of superlinear algorithms increases when approaching the final solution; in our case only ~500 more equivalent iterations are required in the trust region method when the termination tolerance is decreased from 10⁻⁶ to ~2×10⁻¹⁶ (the minimum relative change that can be represented by a double-precision floating-point number).

We thank Dr. Edina Rosta for helpful discussions regarding Eq. (38). This research was supported by the Intramural Research Programs of the NIDDK, NIH, and utilized the high-performance computational capabilities of the Biowulf Linux cluster at the National Institutes of Health, Bethesda, MD (http://biowulf.nih.gov).

In this appendix, the estimate of the statistical error in Eqs. 1, 32, and 36 is generalized to higher-dimensional umbrella sampling in a space spanned by N order parameters x_α, sampled in S simulations with added harmonic bias potentials for simulation i. From each simulation, we obtain an estimate of the N components of the free energy gradient from the mean restraining forces, . As in Eq. 28, we want to combine these gradients into an estimate of the relative free energies G_i=k_BTg_i of the restrained simulations. However, in multiple dimensions the numerically estimated free energy gradients can be integrated through multiple paths, the results of which may not match exactly, due to both statistical errors in the estimated gradients and discretization errors of the integration. A path-independent result can be obtained by minimizing the squared deviation of the free energy differences G_j −G_i between adjacent windows i and j from the corresponding averaged gradients, where the sum is over pairs (i,j) of simulations with nearby restraining centers and , such that the linear approximation in Eq. A1 is applicable. For restraining points on a rectangular grid, minimization of χ² results in a discrete Poisson equation with boundary condition G₁=0 (where we arbitrarily chose i=1 as reference of the free energies, without loss of generality). In general, setting the derivatives of χ² with respect to the G_i to zero results in a set of linear equations that can be written in matrix form asMG = BF, where M amounts to a discrete version of the Laplace operator, G is the vector of (unknown) free energies G_i, and B is a matrix that produces the appropriate linear combination of the restraining forces forming the vector F, consistent with Eq. A1. The formal solution then is G = M⁻¹BF. From the errors in the mean positions of the order parameters, we can again estimate the uncertainties in the corresponding restraining forces, var, and more generally estimate their covariances cov. The covariances of the estimated free energies are then obtained as linear combinations of the covariances in the restraining forces, where c_ikα = (M⁻¹B)_i,kα. Here we assumed that there are no correlations between the results of different simulations i. If we assume further that the restraining forces are uncorrelated, we obtain an estimate of the uncertainty in the free energies, In practice, it may be simpler to circumvent the matrix computations and instead calculate c_ikα directly by repeated minimization of the quadratic target function χ², in what amounts to a variational construction of the Green's function M⁻¹. Specifically, if all are set to zero in Eq. A1, except , then numerical minimization of the resulting χ² with respect to the G_{i_}(as above with the reference G₁=0) directly produces the required solutions c_ikγ = G_i. By repeated minimization of χ² for all k and γ, one can thus build up the entire set of coefficients entering the error formulas Eqs. A2 and A3. We note that for one-dimensional umbrella sampling, with the (i,j) sum in Eq. A1 restricted to nearest neighbors, one recovers the results of the main text, in particular Eqs. 32 and 36.

In this appendix we show that our analysis of the statistical errors can be similarly applied to umbrella sampling with more general biasing potentials w(x,r), such as those with nonuniform spring constants, or anharmonic potentials. In the general case we define the coarse-grained free energy whose derivative is We define The derivative of G_r(r) at r_i is then given by in which is obtained by averaging F(x,r_i) in simulation i with biasing potential w(x,r_i). Then the difference of G_r (r) between any two windows is again given approximately by Eq. (28), with the statistical uncertainty accordingly determined from the estimated variances of the .

For the harmonic biasing potentials w(x,r) = K(x − r)² / 2 discussed in the main text, we have , and thus coincides with the average restraining force of the biasing potential.