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.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.

$<exp(-\beta \Delta U)>_{0}={\frac {\int \,exp(-\beta \Delta U)exp(-\beta \Delta U_{0})dx}{\int \,exp(-\beta \Delta U_{0})dx}}$

$<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 $<...>_{w}$ denotes the weighted ensemble average and the probability distribution function can be shown as:

$\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: $w=exp(-\beta V)$, where $V$ means the bias potential.

$<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) $\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 $\xi _{0}$)

$\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:

$\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)

$\Delta A(\xi )=-k_{B}Tln{\frac {\rho _{w}(\xi )w(\xi _{0})}{\rho _{w}(\xi _{0})w(\xi )}}$

$\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.

$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
• Generalized by Kumar et al.: Weighted Histogram Analysis Method

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 a practical way to do histogram analysis and calculate free energy under discretized language.

Applications

• US combined with FEP

• Potential of Mean Force

Examples

Butane

• Protocol
  • 18 independent simulations
    • 500ps
    • Restraint spring constant = 0.02 kcal/mol-deg
  • WHAM
    • 90 bins (40/bin)
    • Enforced periodcity

• Histograms from Individual Trajectories

• Histograms of Combined Trajectories

• Butane PMF

Ion Channel