[PPT] - Adaptively Biased Molecular Dynamics Volodymyr Babin, Christopher PowerPoint Presentation

SLIDE 1

Adaptively Biased Molecular Dynamics

Volodymyr Babin, Christopher Roland and Celeste Sagui Department of Physics, NC State University, Raleigh, NC 27695-8202

SLIDE 2

1

The Talk Outline:

Problem statement
Metadynamics (+ Applications)
ABMD (+ Applications)

SLIDE 3

2

Problem Statement:

Collective Variable

σ : R3N → Q

Its probability distribution

p(ξ) =

δ [ξ −σ (r1,...,rN)]
Corresponding Free Energy (logdensity of ξ)

f(ξ) = −kBT ln p(ξ), ξ ∈ Q

SLIDE 4

3

The Free Energy f(ξ):

Describes the relative stability of

different states.

Provides useful insights into the

transitions between these states.

SLIDE 5

4

An Example: Ace-(Gly)2-Pro-(Gly)3-Nme

SLIDE 6

5

Long-Lived Conformations:

Globule β-turn

SLIDE 7

6

Collective Variable:

(radius of gyration) Rg =

∑

a

ma MΣ

ra −RΣ

2 RΣ = ∑

a

ma MΣ ra, MΣ = ∑

a

ma

SLIDE 8

7

Why the Rg ?

It provides a sensible description of the

conformations in terms of just one number.

SLIDE 9

8

The Conformations:

Rg ≈ 3.6 ˚ A Rg ≈ 4.4 ˚ A

SLIDE 10

9

The Probability Density:

3 3.5 4 4.5 5 Rg (˚ A)

T = 300 K Canonical Ensemble p(r1,...,rN) ∝ exp

− 1

kBT E(r1,...,rN)

SLIDE 11

10

Molecular Dynamics

200 400 600 800 1000 3 4 5 6

MD time (ns) Rg (˚ A)

T = 300 K

SLIDE 12

11

The Problems:

MD trajectory rarely jumps through the barriers (i.e.,

the MD is bad for sampling from the canonical dis- tribution; this can be “cured” by using, e.g., parallel tempering).

MD trajectory is trapped near the free energy minima

(canonical ensemble).

SLIDE 13

12

The Free Energy

3.5 4 4.5 5 5.5 6 6.5

30
25
20
15
10
5

Rg (˚ A) f(Rg) (kcal/mol)

T = 300 K

(kBT ≈ 0.6 kcal/mol)

≈ 5.5 kBT

SLIDE 14

13

“Classical” Remedies:

Better ways of sampling from the canonical

distribution (replica exchange).

Sampling from a biased distribution with the

bias that can be “undone” afterwards (umbrella sampling).

SLIDE 15

14

Non-Equilibrium Methods

Local Elevation (MD context)
T. HUBER, A. E. TORDA, AND W. F. VAN GUNSTEREN, Local elevation:

a method for improving the searching properties of molecular dynamics simulation, J. Comput. Aided. Mol. Des., 8 (1994), pp. 695–708.

Wang-Landau (MC context)
F. WANG AND D. P. LANDAU, Efficient, multiple-range random walk

algorithm to calculate the density of states, Phys. Rev. Lett., 86 (2001),

pp. 2050–2053.

SLIDE 16

15

Non-Equilibrium Methods

Adaptive Force Bias
E. DARVE AND A. POHORILLE, Calculating free energies using average

force, J. Chem. Phys., 115 (2001), pp. 9169–9183.

J. H´

ENIN AND C. CHIPOT, Overcoming free energy barriers using uncon-

strained molecular dynamics simulations., J. Chem. Phys., 121 (2004),

pp. 2904–2914.
Metadynamics
M. IANNUZZI, A. LAIO, AND M. PARRINELLO, Efficient exploration of

reactive potential energy surfaces using Car-Parrinello molecular dynamics, Phys. Rev. Lett., 90 (2003), pp. 238302–1.

SLIDE 17

16

Ingredients of a Non-Equilibrium Method:

Sampling Device (typically Molecular

Dynamics or Replica-Exchange Molec- ular Dynamics).

and

Non-stationary Biasing potential.

SLIDE 18

17

Evolving Biasing Potential:

t = 0 t = ∞ – biasing potential – “flattened” free-energy

SLIDE 19

18

Metadynamics References

A. LAIO AND M. PARRINELLO, Escaping free-energy minima, Proc.
Natl. Acad. Sci., 99 (2002), pp. 12562– 12566.
M. IANNUZZI, A. LAIO, AND M. PARRINELLO, Efficient exploration of

reactive potential energy surfaces using Car-Parrinello molecular dynamics, Phys. Rev. Lett., 90 (2003), pp. 238302–1.

SLIDE 20

19

Metadynamics Equations

M ¨ ξ +K

ξ −σ [r1,...,rN]
= − ∂

∂ξVh(ξ,t) ma¨ ra −K

ξ −σ [r1,...,rN]

∂ ∂ra σ [r1,...,rN] = Fa[r1,...,rN] ξ – additional dynamical variable harmonically coupled to the collective variable (σ[r1,...,rN]) Vh(ξ,t) – the “hills” potential

SLIDE 21

20 M ¨ ξ +K

ξ −σ [r1,...,rN]
= 0

If the dynamics of ξ is much slower than the dynamics of ra and the harmonic coupling (K) is strong enough, the motion

f ξ is driven by the free energy gradient:

δ (ξ −σ[r1,...,rN]) ≈ exp

− K

2 (ξ −σ[r1,...,rN])2 ∂ ∂ξ f(ξ) ∝ 1 T

t+T

t

dτ (ξ −σ[r1(τ),...,rN(τ)])

SLIDE 22

21

“Umbrella” Potential Vh(ξ,t)

t = 0 t = ∞ Vh(ξ,t) = ∑

τ<t

G

ξ −ξ(τ)
the metadynamics’ “hills” potential

is a sum of tiny bumps placed along the ξ(t) trajectory

SLIDE 23

22

As-Is Metadynamics is O(t2):

The number of terms (bumps) in Vh(ξ,t) (the

“hills” potential) at time t is proportional to t.

Vh(ξ,t) must be computed at every MD step.

SLIDE 24

23

Metadynamics / Applications

E. ASCIUTTO AND C. SAGUI, Exploring Intramolecular Reactions in

Complex Systems with Metadynamics: The Case of the Malonate An- ions, J. Phys. Chem. A, 109 (2005), pp. 7682–7687.

J. G. LEE, E. ASCIUTTO, V. BABIN, C. SAGUI, T. A. DARDEN AND
C. ROLAND, Deprotonation of Solvated Formic Acid: Car-Parrinello and

Metadynamics Simulations, J. Phys. Chem. B, 110 (2006) pp. 2325– 2331.

SLIDE 25

24

V. BABIN, C. ROLAND, T. A. DARDEN AND C. SAGUI, The free energy

landscape of small peptides as obtained from metadynamics with umbrella sampling corrections, J. Chem. Phys., 125 (2006), pp. 204909.

Metadynamics for AMBER (classical MD: significant entropy

contributions require long runs).

Fast implementation using kd-trees for the “hills” potential

(faster than na¨ ıve O(t2), but still not fast enough).

An equilibrium follow-up run to assess and improve the “raw”

metadynamics free energy.

SLIDE 26

25

The Hills Potential

W A Rc = 2W

Vh(ξ,t) =∑

n

An G

R(ξ
ξ (0)

n ,sn)

Wn
G(0)

R(ξ

ξ (0),s) =
∑

α

ξα −ξ (0)

α

sα 2 G(R) =

exp
−1

2R2

+P(R)exp

−1

2R2 c

, R < Rc

0, R ≥ Rc P(R) = 1 2R2

1+ 1

2R2

c − 1

4R2

− 1

2R2

c

1+ 1

4R2

c

−1

SLIDE 27

26

How to check the accuracy ?

SLIDE 28

27

Corrective Follow-Up:

1. Get the (equilibrium) biased probability density:

EB(r1,...,rN) = E(r1,...,rN)+Vh[σ(r1,...,rN)] pB(ξ) =

δ [ξ −σ(r1,...,rN)]
B
2. Use it to correct the free energy:

∆ f(ξ) = −kBT ln pB(ξ) f(ξ) = −Vh(ξ)+∆f(ξ)

SLIDE 29

28

Ace-(Gly)2-Pro-(Gly)3-Nme

(using Rg as the collective variable)

SLIDE 30

29

Metadynamics Trajectory:

50 100 150 3 4 5 6 7

MD time (ns) ξ(t) (˚ A)

SLIDE 31

30

“Raw” Metadynamics:

3 4 5 6 7

60
50
40
30
20
10

Rg (˚ A) −Vh (kcal/mol) 1 × 103 2 × 103 3 × 103 4 × 103 5 × 103

SLIDE 32

31

“Raw” Metadynamics Error:

3 4 5 6 7 2 3 4 5 6 7 8

Rg (˚ A) −kBT ln pB (kcal/mol)

SLIDE 33

32

“Raw” + “Correction”:

3 4 5 6 7

50
40
30
20
10

10 3.5 3.75 4 4.25 4.5

48
46
44
42

Rg (˚ A) f(Rg) (kcal/mol)

SLIDE 34

33

The “raw” error of ≈ 5kBT is unacceptable (it is comparable with the barrier height)! The equilibrium follow-up is therefore cru- cial to get accurate results.

SLIDE 35

34

Tri-Alanine

SLIDE 36

35

“Raw”:

14.5
13.5
12.5
11.5
10.5
9.5
8.5
7.5
6.5
5.5
4.5
3.5
2.5
1.5
0.5
140
60

+60 +140

140
60

+60 +140

ϕ ψ

140
60

+60 +140

ϕ

implicit explicit

SLIDE 37

36

“Raw” + “Correction”:

14.5
13.5
12.5
11.5
10.5
9.5
8.5
7.5
6.5
5.5
4.5
3.5
2.5
1.5
0.5
140
60

+60 +140

140
60

+60 +140

ϕ ψ

140
60

+60 +140

ϕ

implicit explicit

SLIDE 38

37

Can we do better ?

SLIDE 39

38

Vectors of Improvement:

Different biasing strategies (e.g., smoother

in both time and in Q biasing potential).

Better “sampling devices” (e.g., replica ex-

change).

SLIDE 40

39

V. BABIN, C. ROLAND, AND C. SAGUI, Adaptively biased molecular dynamics

for free energy calculations, J. Chem. Phys., 128 (2008), p. 134101.

ma d2ra dt2 = Fa − ∂ ∂ra U

t|σ (r1,...,rN)
∂U(t|ξ)

∂t = kBT τF G

ξ −σ (r1,...,rN)
T. LELI`

EVRE, M. ROUSSET, AND G. STOLTZ, Computation of free energy pro-

files with parallel adaptive dynamics, J. Chem. Phys., 126 (2007), p. 134111.

SLIDE 41

40

Discretization in Q

U(t|ξ) = ∑

m∈ZD

Um(t)B(ξ/∆ξ −m)

SLIDE 42

41

Discretization in Time

Um(t +∆t) = Um(t)+∆tkBT τF G

σ(t)/∆ξ −m
G(ξ) = 48

41   

1− ξ 2

4 2 , −2 ≤ ξ ≤ 2 0,

therwise

SLIDE 43

42

Advantages of the ABMD

It is fast: goes as O(t) since the numerical cost of

the U(t|ξ) does not depend on time t (it is O(1)).

It is memory efficient (if sparse arrays are used

for the Um(t)).

It is convenient (only two parameters: ∆ξ and τF).

SLIDE 44

43

Ace-(Gly)2-Pro-(Gly)3-Nme

(using Rg as the collective variable)

SLIDE 45

44

Reference Free Energies:

3.25 3.5 3.75 4 4.25 4.5 4.75 5 5.25

26
24
22
20
18
16

Rg (˚ A) f(Rg) (kcal/mol)

600 K 543 K 492 K 445 K 403 K 365 K 331 K 300 K

SLIDE 46

45

The RMS Free Energy Error:

ERMS =

1

b−a

b

a

dξ

f1(ξ)− f2(ξ)−∆

2 where ∆ = 1 b−a

b

a

dξ

f1(ξ)− f2(ξ)

SLIDE 47

46

ABMD vs Metadynamics:

50 100 150 200 250 300 1 2 3 4 5 6

MD time (ns) RMS Error (kcal/mol)

MTD ABMD τF = 90 ps 4∆ξ = 0.25˚ A

SLIDE 48

47

The slower, the better:

100 200 300 400 1 2 3 4 5

MD time (ns) RMS Error (kcal/mol) τF = 180 ps τF = 360 ps τF = 720 ps τF = 1440 ps

SLIDE 49

48

Multiple Walkers

∂U(t|ξ) ∂t = kBT τF ∑

α

G

ξ −σ
rα

1 ,...,rα N

(the sum runs over different MD trajectories)
P. RAITERI, A. LAIO, AND F. L. GERVASIO, C. MICHELETTI AND M. PAR-

RINELLO, Efficient Reconstruction of Complex Free Energy Landscapes by Mul-

tiple Walkers Metadynamics, J. Phys. Chem. B, 110 (2006), pp. 3533–3539.

SLIDE 50

49

Multiple Walkers for Fixed τF:

25 50 75 100 125 150 1 2 3 4 5

MD time (ns) RMS Error (kcal/mol) 8 walkers 4 walkers 2 walkers 1 walker

SLIDE 51

50

Replica Exchange – A Better “Sampling Device”

Y. SUGITA, A. KITAO AND Y. OKAMOTO, Multidimensional replica-exchange

method for free-energy calculations, J. Chem. Phys., 113 (2000), pp. 6042–6051.

SLIDE 52

51

N copies (replicas) in parallel.
Each replica may have different temperature

and/or collective variable.

Either stationary (τF = ∞) or evolving biasing

potential on per-replica basis.

SLIDE 53

52

Exchange Probability:

w(m|n) =

1,

∆ ≤ 0 exp(−∆), ∆ > 0 ∆= 1 kBTn − 1 kBTm

Em

p −En p

+

1 kBTm

Um(ξ n)−Um(ξ m)
− 1

kBTn

Un(ξ n)−Un(ξ m)

SLIDE 54

53

Parallel Tempering ABMD

Different temperatures.
Same collective variable(s) in all replicas.

SLIDE 55

54

Parallel Tempering ABMD

3 6 9 12 15 1 2 3 4 5

MD time (ns) RMS Error (kcal/mol) τF = 11.25 ps τF = 22.5 ps τF = 45 ps τF = 90 ps

600 K 543 K 492 K 445 K 403 K 365 K 331 K 300 K

SLIDE 56

55

Another Collective Variable:

0.5 1 1.5 2 2.5 3 0.5 1 1.5 2

x

1 − x6

1 − x12 NOH =

O,H 1−(rOH/r0)6 1−(rOH/r0)12

SLIDE 57

56

r0 = 2.5 ˚ A NOH ≈ 3

SLIDE 58

57

Hamiltonian Replica Exchange: 8 + 1

600 K 543 K 492 K 445 K 403 K 365 K 331 K 300 K

stationary U = U(Rg)

300 K τF < ∞ U = U(NOH, Rg)

NOH =

O,H 1−(rOH/r0)6 1−(rOH/r0)12

SLIDE 59

58

3 4 5 6 7 2 4 6 8

14
13
12
11
10
9
8
7
6
5
4
3
2
1

NOH Rg (˚ A) kcal/mol

SLIDE 60

59

Globular Conformations:

Rg ≈ 3.8 ˚ A, NOH ≈ 3.7 Rg ≈ 3.6 ˚ A, NOH ≈ 6.0

SLIDE 61

60

SLIDE 62

61

Exchange Probability:

w(m|n) =

1,

∆ ≤ 0 exp(−∆), ∆ > 0 ∆= 1 kBTn − 1 kBTm

Em

p −En p

+

1 kBTm

Um(ξ n)−Um(ξ m)
− 1

kBTn

Un(ξ n)−Un(ξ m)

SLIDE 63

62

Hamiltonian Replica Exchange: Explicit Solvent

All replicas at T = 300K.
10 replicas with U = U(rab) (distances

between the carbon atoms separated by at least two amino-acids).

11th replica with U = U(Rg)

SLIDE 64

63

The Free Energy:

3 4 5 6 7

18
15
12
9
6
3

"raw" "raw" + "correction"

Rg (˚ A) f(Rg) (kcal/mol)

SLIDE 65

64

Applications in Progress:

B ↔ Z DNA

(Vadzim Karpusenka and Mahmoud Moradi)

Left−Handed ↔ Right−Handed polyproline

(Mahmoud Moradi)

α ↔ 310 ↔ π helices by AAAAA(AAARA)3A

(Vadzim Karpusenka)

SLIDE 66

65

The ABMD is included in AMBER 10 (and freely available for AMBER 9)

SLIDE 67

66