**From:** Chris Chipot (*chipot_at_ks.uiuc.edu*)

**Date:** Sat Jan 24 2009 - 16:28:37 CST

**Next message:**Falgun Shah: "RMSD and merging trajectories"**Previous message:**Chris Harrison: "Re: CUDA version status"**Messages sorted by:**[ date ] [ thread ] [ subject ] [ author ] [ attachment ]

Jorgen,

the issue of stratification in free-energy calculations along

a reaction coordinate, or model thereof, is recurrently raised

on this list. Intuitively, either in ABF or in umbrella sampling

simulations, one is tempted, for efficiency purposes, to break

down the reaction path into a series of sequential windows. Beyond

hand-waving and heuristic arguments, one can, however, rationalize

this choice through simple considerations. Let us imagine a reaction

path of length L. Convergence of the free energy over the entire

range L is achieved after t0. Let us now consider that the reaction

path is broken down in N windows (non-overlapping with ABF, on

account of the continuity of the average force) of lengths L1,...,LN,

for which convergence is attained after t1',...,tN'. We want to show

that t0 > Σi ti'. For all intents and purposes, let us drop the drift

part of the ABF equation and only retain the Brownian motion of the

latter, dXt = dBt. Let us further assume that the probability

distribution of configurations, Ψ(t,x), follows:

|Ψ(t,x) - Ψ∞(x)| ~ exp[-λ(L) t]

If we compare the free-energy profile at some value of the reaction

coordinate obtained from a single simulation and from a stratification

scheme, the difference is equal to some tolerance, ε. Typically:

exp[-λ(L1) t1'] = ε

which we can rewrite as:

λ(L1) t1' = ε' = -ln ε, or t1' = ε' / λ(L1)

More generally, Σi ti' = ε' Σi 1 / λ(Li) and t0 = ε' / λ(Σi Li).

From the above, it is rather obvious that if L increases, λ should

decrease accordingly. For the sake of arguments, let us now choose

a law where λ varies with 1 / L² and consider the Fokker-Planck

equation for pure diffusion (i.e. in the absence of drift), ∂tΨ = ΔΨ.

When t tends towards ∞, Ψ is expected to be constant. We can write

the probability distribution of configurations at time t as:

Ψ(t,x) = Ψ = Σk Ψk(x) Ak(t)

where Ψk(x) is a solution of ΔΨk = -λk Ψk. From the latter, it

follows that:

Σk dAk/dt Ψk = -Σk Ak λk Ψk

which are independent vectors. From dAk/dt = -Ak λk, one can infer

the solutions:

Ak(t) = Ak(t=0) exp(-λk t)

which are valid over the interval [0;L]. As an example, let us

consider Ψk(x) = sin(kπ/L)x. The second derivative of Ψk(x) is

-(kπ/L)² sin(kπ/L)x and λk = (kπ/L)². For k=1, the eigenvalue is

(π/L)², which is in line with our law for λ in 1 / L². Over the

series of sequential windows, (Σi Li)² = (L1 + L2 + ...)² +

cross-terms > L1² + L2². It follows that (Σi Li)² is always greater

than (Σi Li²). The time required for convergence in a single run,

t0 = ε' (Σi Li)², and Σi ti = ε' Σi Li². As a result, t0 > Σi ti'.

Chris

_______________________________________________________________________

Chris Chipot, Ph.D.

on leave from Nancy Université, CNRS

Theoretical and Computational Biophysics Group

Beckman Institute

University of Illinois at Urbana-Champaign

405 North Mathews Phone: (217) 244-5711

Urbana, Illinois 61801 Fax: (217) 244-6078

E-mail: chipot_at_ks.uiuc.edu

Christophe.Chipot_at_edam.uhp-nancy.fr

Web: http://www.ks.uiuc.edu

http://www.edam.uhp-nancy.fr

To sin by silence when we should protest makes cowards out of men

Ella Wheeler Wilcox

_______________________________________________________________________

**Next message:**Falgun Shah: "RMSD and merging trajectories"**Previous message:**Chris Harrison: "Re: CUDA version status"**Messages sorted by:**[ date ] [ thread ] [ subject ] [ author ] [ attachment ]

*
This archive was generated by hypermail 2.1.6
: Wed Feb 29 2012 - 15:50:25 CST
*