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

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

1 The Talk Outline: • Problem statement • Metadynamics (+ Applications) • ABMD (+ Applications)

2 Problem Statement: • Collective Variable σ : R 3 N �→ Q • Its probability distribution � � p ( ξ ) = δ [ ξ − σ ( r 1 ,..., r N )] • Corresponding Free Energy (logdensity of ξ ) f ( ξ ) = − k B T ln p ( ξ ) , ξ ∈ Q

3 The Free Energy f ( ξ ) : • Describes the relative stability of different states. • Provides useful insights into the transitions between these states.

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

5 Long-Lived Conformations: Globule β -turn

6 Collective Variable: (radius of gyration) m a � � 2 ∑ � R g = r a − R Σ M Σ a m a R Σ = ∑ r a , M Σ = ∑ m a M Σ a a

7 Why the R g ? • It provides a sensible description of the conformations in terms of just one number.

8 The Conformations: R g ≈ 3 . 6 ˚ A R g ≈ 4 . 4 ˚ A

9 The Probability Density: T = 300 K Canonical Ensemble p ( r 1 ,..., r N ) ∝ − 1 � � k B T E ( r 1 ,..., r N ) exp 3 3.5 4 4.5 5 R g (˚ A)

10 Molecular Dynamics 6 5 T = 300 K A) R g (˚ 4 3 0 200 400 600 800 1000 MD time (ns)

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).

12 The Free Energy 0 -5 f ( R g ) (kcal/mol) T = 300 K -10 ( k B T ≈ 0 . 6 kcal/mol) -15 -20 ≈ 5 . 5 k B T -25 -30 3.5 4 4.5 5 5.5 6 6.5 R g (˚ A)

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).

14 Non-Equilibrium Methods • Local Elevation (MD context) T. H UBER , A. E. T ORDA , AND W. F. VAN G UNSTEREN , 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. W ANG AND D. P. L ANDAU , Efficient, multiple-range random walk algorithm to calculate the density of states , Phys. Rev. Lett., 86 (2001), pp. 2050–2053.

15 Non-Equilibrium Methods • Adaptive Force Bias E. D ARVE AND A. P OHORILLE , Calculating free energies using average force , J. Chem. Phys., 115 (2001), pp. 9169–9183. J. H´ ENIN AND C. C HIPOT , Overcoming free energy barriers using uncon- strained molecular dynamics simulations. , J. Chem. Phys., 121 (2004), pp. 2904–2914. • Metadynamics M. I ANNUZZI , A. L AIO , AND M. P ARRINELLO , Efficient exploration of reactive potential energy surfaces using Car-Parrinello molecular dynam- ics , Phys. Rev. Lett., 90 (2003), pp. 238302–1.

16 Ingredients of a Non-Equilibrium Method: • Sampling Device (typically Molecular Dynamics or Replica-Exchange Molec- ular Dynamics). and • Non-stationary Biasing potential.

17 Evolving Biasing Potential: t = 0 t = ∞ – biasing potential – “flattened” free-energy

18 Metadynamics References A. L AIO AND M. P ARRINELLO , Escaping free-energy minima , Proc. Natl. Acad. Sci., 99 (2002), pp. 12562– 12566. M. I ANNUZZI , A. L AIO , AND M. P ARRINELLO , Efficient exploration of reactive potential energy surfaces using Car-Parrinello molecular dynam- ics , Phys. Rev. Lett., 90 (2003), pp. 238302–1.

19 Metadynamics Equations = − ∂ � � M ¨ ξ + K ξ − σ [ r 1 ,..., r N ] ∂ξ V h ( ξ , t ) � ∂ � ξ − σ [ r 1 ,..., r N ] σ [ r 1 ,..., r N ] = F a [ r 1 ,..., r N ] m a ¨ r a − K ∂ r a ξ – additional dynamical variable harmonically coupled to the collective variable ( σ [ r 1 ,..., r N ] ) V h ( ξ , t ) – the “hills” potential

20 � � M ¨ ξ + K ξ − σ [ r 1 ,..., r N ] = 0 If the dynamics of ξ is much slower than the dynamics of r a and the harmonic coupling ( K ) is strong enough, the motion of ξ is driven by the free energy gradient: − K � 2 ( ξ − σ [ r 1 ,..., r N ]) 2 � δ ( ξ − σ [ r 1 ,..., r N ]) ≈ exp t + T ∂ ∂ξ f ( ξ ) ∝ 1 � d τ ( ξ − σ [ r 1 ( τ ) ,..., r N ( τ )]) T t

21 “Umbrella” Potential V h ( ξ , t ) t = 0 t = ∞ the metadynamics’ “hills” potential � � V h ( ξ , t ) = ∑ ξ − ξ ( τ ) G is a sum of tiny bumps placed along τ < t the ξ ( t ) trajectory

22 As-Is Metadynamics is O ( t 2 ) : • The number of terms (bumps) in V h ( ξ , t ) (the “hills” potential) at time t is proportional to t . • V h ( ξ , t ) must be computed at every MD step.

23 Metadynamics / Applications E. A SCIUTTO AND C. S AGUI , 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. L EE , E. A SCIUTTO , V. B ABIN , C. S AGUI , T. A. D ARDEN AND C. R OLAND , Deprotonation of Solvated Formic Acid: Car-Parrinello and Metadynamics Simulations , J. Phys. Chem. B, 110 (2006) pp. 2325– 2331.

24 V. B ABIN , C. R OLAND , T. A. D ARDEN AND C. S AGUI , 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 ıve O ( t 2 ) , but still not fast enough). (faster than na¨ • An equilibrium follow-up run to assess and improve the “raw” metadynamics free energy.

25 The Hills Potential W � A � � ξ ( 0 ) � �� V h ( ξ , t ) = ∑ R ( ξ n , s n ) G ( 0 ) A n G W n � n � � 2 ξ α − ξ ( 0 ) � � � � α � ξ ( 0 ) , s ) = R c = 2 W � ∑ R ( ξ � s α α � � 2 R 2 � � � − 1 − 1 2 R 2 exp + P ( R ) exp , R < R c c G ( R ) = 0 , R ≥ R c � � � � P ( R ) = 1 1 + 1 c − 1 − 1 1 + 1 2 R 2 2 R 2 4 R 2 2 R 2 4 R 2 − 1 c c

26 How to check the accuracy ?

27 Corrective Follow-Up: 1. Get the (equilibrium) biased probability density: E B ( r 1 ,..., r N ) = E ( r 1 ,..., r N )+ V h [ σ ( r 1 ,..., r N )] � � p B ( ξ ) = δ [ ξ − σ ( r 1 ,..., r N )] B 2. Use it to correct the free energy: ∆ f ( ξ ) = − k B T ln p B ( ξ ) f ( ξ ) = − V h ( ξ )+ ∆ f ( ξ )

28 Ace-(Gly) 2 -Pro-(Gly) 3 -Nme (using R g as the collective variable)

29 Metadynamics Trajectory: 7 6 A) ξ ( t ) (˚ 5 4 3 0 50 100 150 MD time (ns)

30 “Raw” Metadynamics: 0 -10 1 × 10 3 3 × 10 3 − V h (kcal/mol) -20 2 × 10 3 4 × 10 3 -30 5 × 10 3 -40 -50 -60 7 3 4 5 6 R g (˚ A)

31 “Raw” Metadynamics Error: 8 − k B T ln p B (kcal/mol) 7 6 5 4 3 2 7 3 4 5 6 R g (˚ A)

32 “Raw” + “Correction”: 10 -42 0 -44 f ( R g ) (kcal/mol) -10 -46 -20 -48 3.5 3.75 4 4.25 4.5 -30 -40 -50 7 3 4 5 6 R g (˚ A)

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

34 Tri-Alanine

35 “Raw”: -0.5 -1.5 -2.5 -3.5 -4.5 -5.5 -6.5 -7.5 -8.5 -9.5 -10.5 -11.5 -12.5 -13.5 -14.5 +140 +60 ψ -60 -140 ϕ ϕ -140 -60 +60 +140 -140 -60 +60 +140 implicit explicit

36 “Raw” + “Correction”: -0.5 -1.5 -2.5 -3.5 -4.5 -5.5 -6.5 -7.5 -8.5 -9.5 -10.5 -11.5 -12.5 -13.5 -14.5 +140 +60 ψ -60 -140 ϕ ϕ -140 -60 +60 +140 -140 -60 +60 +140 implicit explicit

37 Can we do better ?

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).

39 V. B ABIN , C. R OLAND , AND C. S AGUI , Adaptively biased molecular dynamics for free energy calculations , J. Chem. Phys., 128 (2008), p. 134101. d 2 r a d t 2 = F a − ∂ � � t | σ ( r 1 ,..., r N ) m a U ∂ r a ∂ U ( t | ξ ) = k B T � � ξ − σ ( r 1 ,..., r N ) G ∂ t τ F T. L ELI ` EVRE , M. R OUSSET , AND G. S TOLTZ , Computation of free energy pro- files with parallel adaptive dynamics , J. Chem. Phys., 126 (2007), p. 134111.

40 Discretization in Q U ( t | ξ ) = ∑ U m ( t ) B ( ξ / ∆ ξ − m ) m ∈ Z D

41 Discretization in Time U m ( t + ∆ t ) = U m ( t )+ ∆ tk B T � � σ ( t ) / ∆ ξ − m G τ F  � 2 � 1 − ξ 2 � G ( ξ ) = 48 , − 2 ≤ ξ ≤ 2  4 41 0 , otherwise 

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 U m ( t ) ). • It is convenient (only two parameters: ∆ ξ and τ F ).

43 Ace-(Gly) 2 -Pro-(Gly) 3 -Nme (using R g as the collective variable)