From: Chris Chipot (chipot_at_ks.uiuc.edu)
Date: Sat Jan 24 2009 - 16:28:37 CST
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
_______________________________________________________________________
This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:52:17 CST