Reweighting Techniques for Monte Carlo and Molecular Dynamics - PowerPoint PPT Presentation

The 13th KIAS Protein Folding Winter School High 1 Resort, Korea January 19-24, 2014 Reweighting Techniques for Monte Carlo and Molecular Dynamics Simulations I Yuko OKAMOTO （岡本 祐幸） Department of Physics and Structural Biology Research Center Graduate School of Science and Center for Computational Science Graduate School of Engineering and Information Technology Center NAGOYA UNIVERSITY （名古屋大学） e-mail: okamoto{a}phys.nagoya-u.ac.jp URL: http://www.tb.phys.nagoya-u.ac.jp/

Papers 1. A.M. Ferrenberg and R.H. Swendsen, “New Monte Carlo Technique for Studying Phase Transitions,” Physical Review Letters 61 , 2635-2638 (1988). --- Single-Histogram Reweighting Techniques 2. S. Kumar, D. Bouzida, R.H. Swendsen, P.A. Kollman, and J.M. Rosenberg, “The Weighted Histogram Analysis Method for Free-Energy Calculations on Biomolecules. I. The Method,” Journal of Computational Chemistry 13 , 1011-1021 (1992). --- Multiple-Histogram Reweighting Techniques (WHAM) 3. M.R. Shirts and J.D. Chodera, “Statistically Optimal Analysis of Samples from Multiple Equilibrium States,” Journal of Chemical Physics 129 , 124105 (10 pages) (2008). --- A Variant of WHAM (Multistate Bennett Acceptance Ratio Estimator: MBAR)

Contents 1. Metropolis Monte Carlo Method 2. Simulated Annealing and Related Methods 3. Single-Histogram Reweighting Techniques 4. Generalized-Ensemble Algorithms I (MUCA) 5. Multiple-Histogram Reweighting Techniques (WHAM) 6. Multistate Bennett Acceptance Ratio Estimator (MBAR) 7. Generalized-Ensemble Algorithms II (REM)

1. Metropolis Monte Carlo Method

Microcanonical Ensemble ミクロカノニカルアンサンブル Isolated System ： E tot = const 孤立系： E tot = 一定

Canonical Ensemble カノニカルアンサンブル System in Heat Bath (Exchange Energy w/ Heat Bath) 熱浴中の系（熱浴とエネルギーをやりとり）： T = const = 一定

Canonical Ensemble at n(E) Temperature T P B (E) = n(E)W B (E) E Density of States W B (E) = exp(-  E ) E Canonical Probability Distribution E Boltzmann Factor cano

MONTE CARLO State x    x , , ,  , ; , , ,  , q q q q p p p p 1 2 3 1 2 3 N N Generate states one after another.            1 2 3 v v 1        x x x  x x x x j k Suppose at the  -th step the state was x j and the candidate for the (  +1)-th step is x k . Suppose also that the transition probability w satisfies the following: Markov Chain                 v   v 1  1 2 v 1   w x x x x ; x , x ,  , x w x x j k j k

By definition, the transition probability w satisfies             v 1  v  P x P x w x x k j j k j where P (  ) ( x ) is the probability distribution of state x at the  -th step. The equilibrium probability distribution P eq ( x ) should satisfy          P x P x w x x eq k eq j j k j A sufficient condition for this to be satisfied is to impose a detailed balance condition:            P x w x x P x w x x eq j j k eq k k j

One solution to this detailed balance condition is ( Metropolis method ): N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller & E. Teller, J. Chem. Phys. 21 , 1087 (1953).

In Canonical Ensemble, where the equilibrium probability distribution is ∝ the Boltzmann weight factor, we have

Suppose at the  -th step the state was x j and the candidate for the (  +1)-th step is x k .

Example of MC simulations: Lennard-Jones fluid   12 6             E 4         LJ   r r      i j   ij ij r   q q where ij i j   x  , ,  , q q q 1 2 N where Trial Move: q old → q new    q q q i i i     q i          x y z random numbers         0.5 0.5 x y z

Met-Enkephalin: Global Minimum Structure in Gas Phase Tyr-Gly-Gly-Phe-Met (N = 5) E = -12 kcal/mol potential energy: ECEPP/2 degrees of freedom: dihedral angles

Canonical 1000K 30 20 10 E 0 -10 0 50000 100000 150000 200000 MC Sweeps

Canonical 600K 30 20 10 E 0 -10 0 50000 100000 150000 200000 MC Sweeps

Canonical 300K 30 20 10 E 0 -10 0 50000 100000 150000 200000 MC Sweeps

Canonical 50K 30 20 10 E 0 -10 0 50000 100000 150000 200000 MC Sweeps

Cf. Molecular Dynamics Method

MOLECULAR DYNAMICS Newton’s equations of motion: Microcanonical Ensemble  E     m q f i i  q i Nose’s method: Canonical Ensemble at temperature T     E s s         m m m q q f q  i i i i  s s  q i    2  s       2  Qs s m 3 Nk T Q s q    i B    i S. Nose, Mol. Phys. 52 , 255 (1984); J. Chem. Phys. 81 , 511 (1984).

2. Simulated Annealing and Related Methods

徐冷法 (Simulated Annealing) S. Kirkpatrick, C. Gelatt, Jr. & M. Vecchi, Science 220 , 671 (1983). Reproduce a Crystal-Making Process on a Computer SA-2

Application of Simulated Annealing to Systems of Biopolymers • H. Kawai, T. Kikuchi & Y.O., Protein Eng . 3 , 85 (1989). See also: • S. Wilson, W. Cui, J. Moskovitz & K. Schmidt, Tetrahedron Lett. 29 , 4373 (1988). • C. Wilson & S. Doniach, Proteins 6 , 193 (1989). • A. Brunger, J. Mol. Biol. 203 , 803 (1988). • M. Nilges, G. Clore & A. Gronenborn, FEBS Lett. 229 , 317 (1988). For a review see: • Y.O., Recent Res. Devel. In Pure & Applied Chem. 2 , 1 (1998).

Simulated Annealing （徐冷法） P B (E) = n(E)W B (E) P B (E) High T E E P B (E) Low T E E min SA-2

Simulated Annealing: T = 1000 K  50 K Canonical 1000K Canonical 50K 30 20 10 E 0 -10 0 50000 100000 150000 200000 MC Sweeps

C-peptide-1 C-Peptide (N = 13): X-Ray C-Peptide of Ribonuclease A Amino-Acid Sequence: Lys-Glu-Thr-Ala-Ala-Ala-Lys-Phe-Glu-Arg- Gln-His-Met

M. Masuya & Y.O., unpublished. Cp-first

M. Masuya & Y.O., unpublished. Cp-last

C-Peptide (N = 13): X-Ray Level 0 Solvation (Gas Phase) RMSD = 1.9 A Level 1 Solvation (distant- Level 2 Solvation (SASA) dependent dielectric) RMSD = 0.8 A RMSD = 1.4 A M. Masuya & Y.O., unpublished.

Problems of Simulated Annealing (SA) Because we lower the temperature during the simulation, we are always breaking thermal equilibrium (and detailed balance conditions). Hence, thermodynamic quantities obtained from SA are not reliable.