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.

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

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_{w}(\xi)=-k_{B} T ln \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 (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(\xi)=-k_{B} T ln \frac {\rho_{w}(\xi) w(\xi_{0})} {\rho_{w}(\xi_{0}) w(\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 \Delta A(\xi)=\Delta A_{w}(\xi)+k_{B} T (ln w(\xi)-ln w(\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 (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(\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

Umbrella Sampling

Overview

• As in the quasi-static method, the aim of the umbrella sampling simulation is to predict the macroscopic properties of a material from its microscopic state and energy barriers by determining the energy landscape. The method uses molecular dynamic simulations to determine the probability for the system to be in a given conformation. Then the energy landscape is determined by inverting the Boltzmann distribution.
• Umbrella sampling attempts to overcome the sampling problem by modifying the Hamiltonian so that the unfavorable states are sampled sufficiently. The modification can be written as a perturbation.

Derivation

File:Ergodicity.gif

• Add a biasing potential to the total potential energy:

     U(R,q) → U(R,q) + V(q)

     V(q) = V^(k)(q) = k(q-q^(k))²

   where: R : a vector of all coordinates
             q : reaction coordinate

• The unbiased probability distribution P(q₀)

Both sides are clearly equal since we have done nothing more than adding and subtracting V(q) from the exponent and multiplying and dividing by the same integral.

• left ratio

Here we replaced an average done with a weight of exp[-U/T] by an average with our new weight which is exp[-(U+V)/T] that can be tuned to get optimal sampling.

• right ratio

The average in the denominator is (again) performed with the adjusted potential.

• In sum, we have

The denominator is not straightforward to calculate. This is since we sample configuration with a weight of exp[−(U+V)/T]=exp[−U/T]exp[−V/T] while the entity we attempt to average is exp[V/T] . The weight is exponentially small in V while the function we attempt to average is exponentially large in the same variable. This guarantees that most of the sampled points will be a miss, most of the times either the function or the weight will be essentially zero. It is therefore highly desirable to find a way to avoid the calculation of the denominator.

Recover Unbiased Properties

• we can compute the average <δ(q−q₀)>_u in more than one way using two different biasing potentials.

Note that the right hand side of the equation includes only terms that are independent of q₀. So miraculously all the q₀ dependence on the left hand side must cancel out to yield a constant. If we use now the same potentials to consider free energy difference of set of points that are still sampled adequately. We can write as an example the same equation for another point say q₁ in which the right hand side of the equation will remain the same constant.

            

• Free Energy
  

          ßF(q₀)=-ßV(q₀)-ln(<δ(q−q₀)>)_U+V

   Since we only care about the free energy difference for different values of q, we can simply shift the values we got from one of the windows to make the F continuous.

          

• How to stitch the pieces from different windows to obtain the free energy along the reaction coordinate?

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 contribution of each run to a reweighting estimate should be weighed based on the magnitude of errors in their histograms. That is, runs that have greater overlap with the reweighting conditions should contribute more to the estimation of property averages.

Derivation

• Theory

• Derivation of the WHAM expressions

• Back to the WHAM expression...

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