UWHAM and SWHAM

The unbinned weighted histogram analysis method (UWHAM) is an algorithm for multi-state free energy estimation and thermodynamic reweighting developed in the Levy Group by Emilio Gallicchio in collaboration with Zhiqiang Tan of Rutgers University. UWHAM is a binless extension of traditional WHAM.

Suppose \(M\) parallel (independent or coupled) simulations in the canonical ensemble are run at \(M\) \(\lambda\)-states. Let \(X_{i}^{(\alpha)}\) be the \(i^{th}\) microstate observed at the \(\alpha^{th}\) \(\lambda\)-state. The probability of observing \(X_{i}^{(\alpha)}\) at the \(\gamma^{th}\) \(\lambda\)-state is

\begin{equation} P_{\gamma}(X_{i}^{(\alpha)})\sim\frac{q_{\gamma}(\{x\}_{\alpha i})}{Z_{\gamma}} =\frac{\exp\{-\beta_{\gamma}E_{\gamma}(\{x\}_{\alpha i})\}}{Z_{\gamma}}\,, \end{equation}

where \(q_{\gamma}(\{x\}_{\alpha i})=\exp\{-\beta_{\gamma}E_{\gamma}(\{x\}_{\alpha i})\}\) is the unnormalized probability of observing \(X_{i}^{(\alpha)}\) during a simulation in the canonical ensemble at the \(\gamma^{th}\) \(\lambda\)-state; \(\{x\}_{\alpha i}\) is the coordinates of the microstate \(X_{i}^{(\alpha)}\); \(\beta_{\gamma}\) is the inverse temperature of the \(\gamma^{th}\) \(\lambda\)-state; \(E_{\gamma}(\{x\}_{\alpha i})\) is the potential energy of the microstate \(X_{i}^{(\alpha)}\) at the \(\gamma^{th}\) \(\lambda\)-state; and \(Z_{\gamma}\) is the partition function of the \(\gamma^{th}\) \(\lambda\)-state. The UWHAM equations are

\begin{align} \hat{Z}_{\alpha} & =\sum_{i=1}^{N}q_{\alpha}(u_{i})\hat{\Omega}(u_{i})\nonumber \\ \hat{\Omega}(u_{i}) & =\frac{1}{\sum_{\kappa=1}^{M}N_{\kappa}\hat{Z}_{\kappa}^{-1}q_{\kappa}(u_{i})}\,,\label{eq:uwham} \end{align}

where \(N=\sum_{\alpha=1}^{M}N_{\alpha}\) is the total number of observations; \(u_{i}\) is the reduced (energy) coordinate of the \(i^{th}\) observation and \(\Omega(u_{i})\) is its density of state.

There is a computational bottleneck in scaling up UWHAM. Suppose there are \(M\) \(\lambda\)-states and totally \(N\) observations to analyze. At minimum,numerical solution of the UWHAM equations requires evaluating M unnormalized density functions for each observation. The total number of function evaluations is of order \(nM^2\), where \(n=N/M\) is the average sample size per \(\lambda\)-state. These unnormalized density values need to be either computed during every iteration of the numerical solution or precomputed and stored in memory. Such a high computational cost presents a serious limitation on the use of UWHAM for large-scale simulations. To remove the computational bottleneck, Levy group and Zhiqiang Tan developed stochastic solver of UWHAM equations by using the protocols of Generalized-ensemble algorithms.

When applying UWHAM and its stochastic solvers, the base assumption is that the simulations at each λ-state are “approximately” equilibrated. However, this assumption might not always hold. To handle such situations, we developed an analysis tool called Stratified-UWHAM to compute free energy and expectations from a multi-state ensemble when the simulations at a subset of λ-states are far from being equilibrated.

The Stratified-UWHAM algorithm is based on coarse-graining. Firstly, we coarse-grained the whole phase space into macrostates, which correspond to energy basins separated by free energy barriers. Then all the \(\lambda\)-states in the system are divided into two groups \((S_{1}\:,S_{2})\). For the \(\lambda\)-states in the first group \(S_{1}\), the simulations are approximately equilibrated among all the macrostates. Namely, all the macrostates are well connected at each \(\lambda\)-state in group \(S_{1}\). For the \(\lambda\)-states in the second group \(S_{2}\), the simulations are only locally equilibrated. Namely, the coarse-grained set of macrostates forms a disconnected network of microstate clusters at each \(\lambda\)-state in group \(S_{2}\). We assume that the set of observations \(\{X_{i}^{(\alpha)}:i=1,2,\cdots N_{\alpha}\}\) observed at the \(\lambda\)-states in the \(S_{1}\) group is independent drawn from distribution

\begin{equation} P(X_{i}^{(\alpha)})\sim\frac{q_{\alpha}(\{x\}_{\alpha i})}{Z_{\alpha}}\,\,\,\textrm{for}\,\alpha\in S_{1}\,. \label{eq:dist_S1} \end{equation}

And the set of observations \(\{X_{i}^{(\alpha)}:X_{i}^{(\alpha)}\in R_{k},i=1,2,\cdots N_{\alpha}\}\) observed at the \(\lambda$\)-states in the \(S_{2}\) group is independent drawn from distribution

\begin{equation} P(X_{i}^{(\alpha)}|(X_{i}^{(\alpha)}\in R_{k}))\sim\frac{q_{\alpha k}(\{x\}_{\alpha i})}{Z_{\alpha k}}\,\,\,\textrm{for}\,\alpha\in S_{2}\,, \label{eq:dist_S2} \end{equation}

where \(q_{\alpha k}(\{x\}_{\alpha i})=q_{\alpha}(\{x\}_{\alpha i})\delta(\{x\}_{\alpha i}\in R_{k})\). Here \(\delta(\{x\}_{\alpha i})\in A\) is the indicator function for a macrostate cluster \(A\); \(Z_{\alpha}\) and \(Z_{\alpha k}\) are the partition functions. The likelihood function of the simulation data of this model is

\begin{equation} \prod_{\alpha\in S_{2}}\prod_{k}\prod_{i:X_{i}^{(\alpha)}\in R_{k}}\left[\frac{1}{Z_{\alpha k}}q_{\alpha}(u_{\alpha i}) \Omega(u_{\alpha i})\right]\times\prod_{\alpha\in S_{1}}\prod_{i=1}^{N_{\alpha}} \left[\frac{1}{Z_{\alpha}}q_{\alpha}(u_{\alpha i})\Omega(u_{\alpha i})\right]\,, \label{eq:SUWHAM_likelihood} \end{equation}

where \(u_{\alpha i}\) is the reduced (energy) coordinate of the microstate \(X_{i}^{(\alpha)}\).

The estimating equations from the maximization of Eq.(\ref{eq:SUWHAM_likelihood}) can be solved in the form of UWHAM equations Eq.(\ref{eq:uwham}) with an expanded set of \(\lambda\)-states. When solving the Stratified-UWHAM equations by solving the original UWHAM equations with an expanded set of \(\lambda\)-states, the total number of \(\lambda\)-states of the system increases from \(M_{1}+M_{2}\) to \(M_{1}+\sum_{\alpha=1}^{M_{2}}K_{\alpha}\), where \(K_{\alpha}\) is the total number of disconnected macrostate clusters at the \(\alpha^{th}\) \(\lambda\)-state. Apparently, it is more computationally expensive to solve the Stratified-UWHAM equations. We proposed a stochastic solver for the Stratified-UWHAM equations called Stratified RE-SWHAM, which is illustrated in the following figure. Stratified RE-SWHAM is a mutant of RE-WHAM. There are two differences between these two algorithms. First, when the observations observed at each \(\lambda\)-state are collected to serve as the database for that \(\lambda\)-state, each observation is tagged by the macrostate which it belongs to. The second difference is the Move procedure in the replica exchange cycle. For the \(\lambda\)-states in the \(S_{1}\) group, an observation in the database is randomly chosen with equal probability to associate with the replica at that \(\lambda\)-state. However, for the \(\lambda\)-states in the \(S_{2}\) group, an observation which is in the same connected macrostate cluster as the observation associated with the replica at that \(\lambda\)-state is chosen with equal probability. In Ref 4, we proved that the output os Stratified RE-SWHAM at a \(\lambda\)-state is the estimate of the equilibrium distribution of that \(\lambda\)-state according to the Stratified-UWHAM equations.

• UWHAM
• SWHAM
• Stratified RE-SWHAM

The python scripts used in these tutorials were tested by "Python 3.6.3 :: Anaconda, Inc.".

• 01) Umbrella Sampling: Double Well Potential
• 02) Umbrella Sampling: Alanine Dipeptide
• 03) Temperature Replcia Exchange: Alanine Dipeptide
• 04) Free Energe Perturbation: Solvation Free Energy
• 05) Hamitonian Replica Exchange: BEDAM and Binding
• 06) Temperature and Hamiltonian Replica Exchange
• 07) Two Binding Modes of the \(\beta\)-cyclodextrin Heptanoate Complex