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.

Biased Potential

Why do we need biased potential?

How to determine a biased potential?

• A basic but powerful method to improve the efficiency of such computations is to split the interval of computation along the reaction coordinate into subintervals, an approach termed stratification.
• In addition, in each window a biasing potential can be used in order to improve the sampling further. In this way, the energy barriers inside each window can be made smaller.
• If this biasing potential is a function of V(x), the new Hamiltonian for the system is H(x)+V(x).
  • In general the biasing potential needs to be guessed beforehand or can be gradually improved using an iterative refinement process.
    • Run a short simulation to estimate the free energy, A⁽⁰⁾(x).
    • Use this estimate to bias the system using V⁽⁰⁾(x)=-A⁽⁰⁾(x).
    • With this first bias, we improve the sampling and obtain a more accurate approximation of the free energy A⁽¹⁾(x).
    • The biasing potential V⁽ⁱ⁾(x) can be gradually improved in this fashion. Within each window, we obtain a relatively uniform sampling. This leads in general to small statistical errors.

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