Ammm:Umbrella sampling and WHAM
Reaction Coordinate
- A reaction coordinate is an abstract one-dimensional coordinate which represents progress along a reaction pathway.
- Real coordinate system: bond length, bond angle, torsion
- Non geometric parameters: bond order, Hydrogen bonds, RMSD
- Reaction coordinates are often plotted against free energy to demonstrate in some schematic form the potential energy profile associated to the reaction.
Umbrella Sampling (US)
Umbrella sampling is a method in computational physics and chemistry, which can enhance the sampling of different systems when it is hard to realize ergodicity due to the energy landscape. It involves the importance sampling in statistical mechanics and was first recommended by Torrie and Valleau in 1977[1].In many cases, the configuration space may have multiple high-energy barriers which are hard to conquer in the limited simulation time with common molecular dynamics (MD). Therefore, it is of high probability that the trajectory is trapped in the local minimum (local state) and leaves other inaccessible important states unsampled. In order to remove the impact of the energy barriers (activation case), the researchers can chosose the umbrella sampling (US) to bridge the gap between different configuration states by adding an extra bias potential. This bias potential should cancel out the influence of these energy barriers so that the energy landscape will become more smooth, making it easier to reach other states.
Why do we need biased potential?
Features of method
- Pre-determined collective variables are required to describe the configuration space related to reaction.
- The target of this method is to handle activation problems in sampling, where high energy barriers can be smoothened by adding extra bias potential. (Convert activation into diffusion)
- The extra bias potential causes the loss of temporal properties.
- This method supports parallel MD simulations with different bias potentials at the same time so that computational cost is reduced.
- This method demands post-processing to remove the effects of bias potential.
Weighted ensemble
The original ensemble can be modified by using a weight function. The weight function should be positive here.
- Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle <exp(-\beta \Delta U)>_{0}={\frac {\int \,exp(-\beta \Delta U)exp(-\beta \Delta U_{0})dx}{\int \,exp(-\beta \Delta U_{0})dx}}}
- Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle <exp(-\beta \Delta U)>_{0}={\frac {\int \,exp(-\beta \Delta U)ww^{-1}exp(-\beta \Delta U_{0})dx}{\int \,exp(-\beta \Delta U_{0})ww^{-1}dx}}={\frac {<w^{-1}exp(-\beta \Delta U)>_{w}}{<w^{-1}>_{w}}}}
In the equation above, the Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle <...>_{w}} denotes the weighted ensemble average and the probability distribution function can be shown as:
- Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle \rho _{w}(x)={\frac {wexp(-\beta \Delta U_{0})}{\int \,wexp(-\beta \Delta U_{0})dx}}}
To reflect the connection between bias potential and the weight function, we can represent the weight function as an exponential form: Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle w=exp(-\beta V)} , where means the bias potential.
- Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle <exp(-\beta \Delta U)>_{0}={\frac {<exp(\beta V)exp(-\beta \Delta U)>_{w}}{<exp(\beta V)>_{w}}}}
Free energy calculation
After defining a collective variable (CV) Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle \xi } , we can calculate the free energy at any point along the direction of CV by using the probability distribution function. (The reference point is chosen at Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle \xi _{0}} )
- Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle \Delta A(\xi )=-k_{B}Tln{\frac {\rho _{0}(\xi )}{\rho _{0}(\xi _{0})}}}
When the bias potential is added, the free energy change under the weighted ensemble can be given by:
- Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle \Delta A_{w}(\xi )=-k_{B}Tln{\frac {\rho _{w}(\xi )}{\rho _{w}(\xi _{0})}}}
In this case, the original free energy change should be recovered by combining these two equations above: (the denominators have been offset in the ratio)
- Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle \Delta A(\xi )=\Delta A_{w}(\xi )+k_{B}T(lnw(\xi )-lnw(\xi _{0}))=\Delta A_{w}(\xi )-[V(\xi )-V(\xi _{0})]}
The best bias potential can help remove all the possible barriers on the energy landscape, which means the free energy change under this bias potential should be 0 everywhere. So, it is obvious to know that the best bias potential must counterbalance the true free energy.
- Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle V(\xi )-V(\xi _{0})=-\Delta A(\xi )}
Challenges of Umbrella Sampling
- Pick up the best bias potential along the CVs to realize the sufficiently flat energy landscape
- Ensure complete scanning of all the configuration space along the CVs under bias potential
- Perform discretization of the probability distribution along CVs
WHAM
Historty
- Initially developed in 1989 by Ferrenberg and Swenden: Ferrenberg-Swenden reweighting [2]
- Generalized by Kumar et al.: Weighted Histogram Analysis Method [3]
Basic Idea
The Weighted Histogram Analysis Method (WHAM) is a practical technique for the calculation of potentials of mean force (PMFs) from a group of umbrella sampling simulations. It can be applied to analyze the energy data based on not only real collective variables (have relations to the atom coordinates) but also virtual variables (temperature or interactions in the system). This technique can give out a systematic way to improve the choice of the bias potential from a set of different bias potentials so that it yields uniform sampling along the collective variables. Meanwhile, the multiple simulations can start from different points along collective variables, then the full sampling can be realized. Also, the underlying theory of this method provides us with a practical way to do histogram analysis and calculate free energy under discretized language.
Discretization of configuration space along collective variable
Assume we have run Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle R} simulations at the same time, with each using a different or same bias potential and temperature. Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle }
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle U_{k}=U_{0}+V_{k}(\xi),k=1,...,R}
The whole configuration space can be divided into Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle L} bins. The middle bins share the same width Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle d} while the beginning and end bins have the width Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle d/2} . After the division, the configuration space becomes a histogram graph, and the energy value for each bin Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle i} in simulation Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle k} should be given by the following equation: (energy for middle bins is set to be the value of the middle point of this bin, while the energy of the beginning bin is set to be the value at the starting point and the energy of the end bin is set to be the value at the endpoint)
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle V_{k,i}= \begin{cases} V_{k}(\xi_{i}), & (i=1,L+1) \\ V_{k}(\frac{\xi_{i}+\xi_{i+1}} {2}), & (i=2,...,L) \end{cases}}
Then, the bias potential energy of the simulation Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle k} along collective variable can be expressed in discretized language:
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle V_{k}(\xi)=\sum_{i=1}^L V_{k,i} \chi_{[\xi_{i},\xi_{i+1}]}(\xi)}
Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \chi_{[\xi_{i},\xi_{i+1}]}} is characteristic function: if Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \xi} is in the range Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle [\xi_{i},\xi_{i+1}]} then it takes value Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle 1} , otherwise it is Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle 0} .
From the operation of dividing the whole configuration space into discrete bins, it is obvious that more bins indicate higher resolution, but more snapshots from the simulations are required. To make the expression more clear, we directly use Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle V_{k,i}(\xi)} to denote the function Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle V_{k,i} \chi_{[\xi_{i},\xi_{i+1}]}(\xi)} in the following.
State density & probability distribution function
The most important property in WHAM is state density (or energy state density) Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Omega} . In an NVT ensemble, the generalized state density can be shown as:
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Omega(N,V,U)=\frac{1} {h^{3N} N!} \iint\limits_{V^{3N}} \delta (H(\mathbf{q},\mathbf{p})-U)\, d\mathbf{q}\,d\mathbf{p}}
then Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Omega(N,V,U)dU} represents the volume of microstates in Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle [U-dU/2,U+dU/2]} . The partition function can be illustrated by using state density:
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle Q(N,V,T)=\int\limits exp(-\beta U) \Omega(N,V,U)\, dU}
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle Q(N,V,T)=\frac{1} {h^{3N} N!} \iint\limits_{V^{3N}} exp(-\beta U(\mathbf{q},\mathbf{p}))\, d\mathbf{q}\,d\mathbf{p}=\frac{1} {\Lambda^{3N} N!} \int\limits exp(-\beta U(\mathbf{q}))\, d\mathbf{q}=\frac{1} {\Lambda^{3N} N!} Z(N,V,T)}
Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle Z(N,V,T)} is the configurational integral. From these equations above, the probability distribution function can be expressed as:
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \rho(U;N,V,T)=\frac{exp(-\beta U) \Omega(N,V,U)} {Z(N,V,T)}}
In WHAM, the probability and state density can be calculated by using statistical data: Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle N_{k}(\xi)} is the number of counts at Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \xi} in simulation Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle k} , Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle n_{k}} is the total number of snapshots in simulation Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle k} .
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \rho_{k}(U(\xi),T_{k})=\frac{N_{k}} {n_{k}}=\frac{exp(-\beta_{k} U) \Omega_{k}(U)} {Z_{k}}}
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Omega_{k}(U)=a_{k} \frac{N_{k}} {n_{k}} exp(-\beta_{k} U)}
Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle a_{k}} denotes the proportionality constant unique for k-th simulation (related to configurational integral), which can be modified further to make the following equations more pretty: Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle a_{k}=exp(-\beta_{k} f_{k})} .
Derivation of WHAM[4]
When bias potential is shown explicitly, the state density becomes:
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Omega_{k}(\xi)=\frac{N_{k}} {n_{k}} exp(\beta_{k} U_{0}(\xi)+\sum_{i=1}^L \beta_{k} V_{k,i}(\xi)-\beta_{k} f_{k})}
Because the original energy distribution is the same for all the simulations and will be reduced in the following deduction, then we can remove that in state density.
- Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle \Omega _{k}(\xi )={\frac {N_{k}}{n_{k}}}exp(\sum _{i=1}^{L}\beta _{k}V_{k,i}(\xi )-\beta _{k}f_{k})}
The state density with the best bias potential can be given by a linear combination of state density functions from all Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle R} simulations, with each one being assigned a weight . The sum of all the weights should be .
- Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle \Omega (\xi )=\sum _{k=1}^{R}\omega _{k}\Omega _{k},\sum _{k=1}^{R}\omega _{k}=1}
To get the best state density, we need to minimize the mean square deviation (MSD) Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle \delta ^{2}\Omega }
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \delta^2 \Omega=\sum_{k=1}^R \omega_{k}^2 \delta^2 \Omega_{k}}
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \delta^2 \Omega_{k}=\frac{1} {n_{k}^2} exp(\sum_{i=1}^L \beta_{k} V_{k,i}-\beta_{k} f_{k}) \delta^2 N_{k}}
The MSD of Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle N_{k}} can be given by multiplication of a correlation factor and its average[5]: (the correlation factor is roughly equal if the division and running time are the same for all simulations, so it can be canceled)
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \delta^2 N_{k}(\xi)=g_{k} \bar {N_{k}}(\xi)}
An estimate can be made to obtain the average Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \bar N_{k}} :
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Omega(\xi)=\frac{\bar {N_{k}}(\xi) exp(\sum_{i=1}^L \beta_{k} V_{k,i}(\xi)-\beta_{k} f_{k})} {n_{k}}}
which can be compared with Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Omega_{k}(\xi)=\frac{ N_{k}(\xi) exp(\sum_{i=1}^L \beta_{k} V_{k,i}(\xi)-\beta_{k} f_{k})} {n_{k}}} .
Then: Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \bar {N_{k}}(\xi)=n_{k} \Omega(\xi) exp(-\sum_{i=1}^L \beta_{k} V_{k,i}(\xi)+\beta_{k} f_{k})} .
The minimization can be finished by using the method of Lagrangian multipliers:
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \frac{\partial (\delta^2 \Omega)} {\partial \omega_{k}}=\frac{\partial} {\partial \omega_{k}} (\delta^2 \Omega+\lambda (\sum_{i} \omega_{i}-1))}
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle 2\omega_{k} \delta^2 \Omega_{k}+\lambda=0(k=1,...,R)}
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \omega_{1} \delta^2 \Omega_{1}=\omega_{2} \delta^2 \Omega_{2}=...=\omega_{R} \delta^2 \Omega_{R}=C}
Because Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \sum{i} \omega_{i}=1} , then Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle C=\frac{1} {\sum_{k} \frac{1} {\delta^2 \Omega_{k}}}} .
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \omega_{k}=\frac{\delta^2 \Omega_{k}} {\sum_{k} \frac{1} {\delta^2 \Omega_{k}}}=\frac{g_{k}^{-1} n_{k} exp(-\sum_{i=1}^L \beta_{k} V_{k,i}(\xi)+\beta_{k} f_{k})} {\sum_{m=1}^R g_{m}^{-1} n_{m} exp(-\sum_{i=1}^L \beta_{m} V_{m,i}(\xi)+\beta_{m} f_{m})}} .
Then:
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Omega(\xi)=\sum_{k=1}^R \omega_{k} \Omega_{k}(\xi)=\frac{\sum_{k=1}^R g_{k}^{-1} N_{k}(\xi)} {\sum_{m=1}^R g_{m}^{-1} n_{m} exp(-\sum_{i=1}^L \beta_{m} V_{m,i}(\xi)+\beta_{m} f_{m})}}
In bin Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle i} :
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Omega_{,i}=\frac{\sum_{k=1}^R g_{k}^{-1} N_{k,i}} {\sum_{m=1}^R g_{m}^{-1} n_{m} exp(-\sum_{i=1}^L \beta_{m} V_{m,i}+\beta_{m} f_{m})}}
The un-normalized density at bin Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle i} in simulation Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle K} can be given by:
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \tilde \rho_{K,i}=\frac{exp(-\beta_{K} V_{K,i}) \sum_{k=1}^R g_{k}^{-1} N_{k,i}} {\sum_{m=1}^R g_{m}^{-1} n_{m} exp(-\sum_{i=1}^L \beta_{m} V_{m,i}+\beta_{m} f_{m})}\ \ (*)}
Because Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle f_{K}} is related to configurational integral:
Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle exp(-\beta_{K} f_{K})=\sum_{i=1}^L \tilde \rho_{K,i}\ \ (**)}
The most common case is that all the temperatures are the same, and the correlation factors are cancelled out. Then the equations Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle (*)} and Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle (**)} can be simplified into:
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \tilde \rho_{K,i}=\frac{exp(-\beta V_{K,i}) \sum_{k=1}^R N_{k,i}} {\sum_{m=1}^R n_{m} exp(-\sum_{i=1}^L \beta V_{m,i}+\beta f_{m})}\ \ (*)}
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle exp(-\beta f_{K})=\sum_{i=1}^L \tilde \rho_{K,i}\ \ (**)}
Self-consistent algorithm
In equations Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle (*)} and Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle (**)} , the unknown values are Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle f_{K}} and Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \tilde \rho_{K,i}} . We can work out their values by using the self-consistent method.
At first, we can assume Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle f_{K}=0(K=1,...,R)} and put them into Equation Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle (*)} to obtain Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \tilde \rho_{K,i}(K=1,..,R;i=1,...,L)} .
Next step: Put all Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \tilde \rho_{K,i}} into Equation Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle (**)} to get the values of Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle f_{K}} .
Repeat the two steps above until the deviations of all these values are smaller than pre-set tolerance. The true free energy change between any two bins in simulation Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle K} without bias potential can be calculated by:
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Delta A_{K,j \to i}=-\beta_{K}^{-1} ln(\frac{\rho_{K,i}} {\rho_{K,j}})}
Procedure of umbrella sampling
- Define collective variables appropriate to describe the reaction space
- Divide the whole reaction space into bins with proper resolution
- Generate different configurations at each bin by using steered MD (or directly do pulling during different simulations)
- Set bias potential for each configuration (commonly used: harmonic restrain)
- Run multiple MDs with different bias potentials and record the trajectories & energy
- Put data into the WHAM processing program (use the module from Gromacs, …)
Examples
Butane[6]

- Protocol
- 18 independent simulations
- 500ps
- Restraint spring constant = 0.02 kcal/mol-deg
- WHAM
- 90 bins (40/bin)
- Enforced periodcity
- 18 independent simulations
- Histograms from Individual Trajectories
- Histograms of Combined Trajectories
- Butane PMF
Ion Channel [7]

Tutorial of US & WHAM
The Gromacs provides us with a good tutorial to show the way to perform Umbrella Sampling & WHAM post-processing: http://www.mdtutorials.com/gmx/umbrella/ [8]
References
- ↑ Torrie, G. M.; Valleau, J. P. (1977). "Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling". Journal of Computational Physics. 23 (2): 187–199. doi:10.1016/0021-9991(77)90121-8
- ↑ Ferrenberg; Swenden, Optimized Monte Carlo data Analysis, Phys.Rev.Lett. 1989, 63, 1195-1198
- ↑ Kumar S, Bouzida D, Swendsen RH, Kollmann PA, Rosenberg JM: The weighted histogramfckLRanalysis method for free-energy calculations on biomolecules. I. The method. J Comp ChemfckLR1992, 13: 1011-1021.
- ↑ Roux B: Extension to the weighted histogram analysis method: combining umbrella sampling with free energy calculations. Comp Phys Commu 2001, 135: 40-57.
- ↑ Ferrenberg, A. M., and Swendsen, R. H. (1989) Optimized Monte Carlo Data Analysis. Phys. Rev. Lett 63, 1195–1198.
- ↑ membrane.urmc.rochester.edu/wham/wham_talk.pdf
- ↑ Zhifeng Jing, Joshua A. Rackers , Lawrence R. Pratt, Chengwen Liu, Susan B. Rempe and Pengyu Ren. Thermodynamics of ion binding and occupancy in potassium channels. Chem. Sci., 2021, 12, 8920-8930
- ↑ Justin A. Lemkul and David R. Bevan. Assessing the Stability of Alzheimer's Amyloid Protofibrils Using Molecular Dynamics. J. Phys. Chem. B 2010, 114, 4, 1652–1660