Configurational-Bias Monte Carlo N interacting particles in volume V - PowerPoint PPT Presentation

Configurational-Bias Monte Carlo Thijs J.H. Vlugt [1] Random Sampling versus Metropolis Sampling (1) Configurational-Bias Monte Carlo N interacting particles in volume V at temperature T Thijs J.H. Vlugt Professor and Chair Engineering Thermodynamics Delft University of Technology Delft, The Netherlands t.j.h.vlugt@tudelft.nl January 7, 2019 January 9, 2020 • vector representing positions of all particles in the system: r N • total energy: U ( r N ) • statistical weight of configuration r N is exp[ � β U ( r N )] with β = 1 / ( k B T )

Configurational-Bias Monte Carlo Thijs J.H. Vlugt [2] Configurational-Bias Monte Carlo Thijs J.H. Vlugt [3] Random Sampling versus Metropolis Sampling (2) Random Sampling versus Metropolis Sampling (3) N interacting particles in volume V at temperature T Computing the ensemble average h · · · i of a certain quantity A ( r N ) pair interactions u ( r ij ) • Random Sampling of r N : P n i =1 A ( r N ⇥ � β U ( r N ⇤ i ) exp i ) h A i = lim P n ⇥ � β U ( r N ⇤ i =1 exp i ) n !1 Usually this leads to h A i = “ 0 ” / “ 0 ” = ??? • Metropolis sampling; generate n configurations r N with probability proportional ⇥ � β U ( r N ⇤ N � 1 N to exp i ) , therefore: U ( r N ) X X X = u ( r ij ) = u ( r ij ) i =1 j = i +1 i<j P n i =1 A ( r N i ) h A i = lim 1 Z d r N exp � β U ( r N ) n ⇥ ⇤ n !1 Q ( N, V, T ) = Λ 3 N N ! F ( N, V, T ) = � k B T ln Q ( N, V, T )

Configurational-Bias Monte Carlo Thijs J.H. Vlugt [4] Configurational-Bias Monte Carlo Thijs J.H. Vlugt [5] Simulation Technique (1) Simulation Technique (2) Bottoms up What is the ratio of red wine/white wine in the Netherlands?

Configurational-Bias Monte Carlo Thijs J.H. Vlugt [6] Configurational-Bias Monte Carlo Thijs J.H. Vlugt [7] Simulation Technique (3) Simulation Technique (4)

Configurational-Bias Monte Carlo Thijs J.H. Vlugt [8] Configurational-Bias Monte Carlo Thijs J.H. Vlugt [9] Simulation Technique (5) Simulation Technique (6) Bottoms up Liquor store Shoe shop

Configurational-Bias Monte Carlo Thijs J.H. Vlugt [10] Configurational-Bias Monte Carlo Thijs J.H. Vlugt [11] Metropolis Monte Carlo (1) Metropolis Monte Carlo (2) How to generate configurations r i with a probability proportional to Whatever our rule is to move from one state to the next, the equilibrium distribution N ( r i ) = exp[ � β U ( r i )] ??? should not be destroyed N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth. A.H. Teller and E. Teller, ”Equation of State Calculations by Fast Computing Machines,” J. Chem. Phys., 1953, 21, 1087-1092.

Configurational-Bias Monte Carlo Thijs J.H. Vlugt [12] Configurational-Bias Monte Carlo Thijs J.H. Vlugt [13] Move from the old state ( o ) to a new state ( n ) and back Detailed Balance (1) Requirement (balance): new 1 X X N ( o ) [ N ( n ) α ( n ! o )acc( n ! o )] [ α ( o ! n )acc( o ! n )] = new 5 n n new 2 Detailed balance : much stronger condition old N ( o ) α ( o ! n )acc( o ! n ) = N ( n ) α ( n ! o )acc( n ! o ) for every pair o , n new 4 new 3 new 1 new 5 leaving state o = entering state o new 2 X X N ( o ) [ α ( o ! n )acc( o ! n )] = [ N ( n ) α ( n ! o )acc( n ! o )] old n n N ( i ) : probability to be in state i (here: proportional to exp[ � β U ( r i )] ) α ( x ! y ) : probability to attempt move from state x to state y new 4 new 3 acc( x ! y ) : probability to accept move from state x to state y

Configurational-Bias Monte Carlo Thijs J.H. Vlugt [14] Configurational-Bias Monte Carlo Thijs J.H. Vlugt [15] Metropolis Acceptance Rule Detailed Balance (2) N ( o ) α ( o ! n )acc( o ! n ) = N ( n ) α ( n ! o )acc( n ! o ) General: acc( o ! n ) acc( n ! o ) = X • α ( x ! y ) ; probability to select move from x to y Choice made by Metropolis (note: infinite number of other possibilities) • acc( x ! y ) ; probability to accept move from x to y • often (but not always); α ( o ! n ) = α ( n ! o ) acc( o ! n ) = min(1 , X ) Therefore (note that ∆ U = U ( n ) � U ( o ) ): Note than min( a, b ) = a if a < b and b otherwise acc( n ! o ) = α ( n ! o ) exp[ � β U ( n )] acc( o ! n ) α ( o ! n ) exp[ � β U ( o )] = α ( n ! o ) α ( o ! n ) exp[ � β ∆ U ] • always accept when X � 1 In case that α ( o ! n ) = α ( n ! o ) • when X < 1 , generate uniformly distributed random number between 0 and 1 and accept or reject according to acc( o ! n ) acc( o ! n ) acc( n ! o ) = exp[ � β ∆ U ]

Configurational-Bias Monte Carlo Thijs J.H. Vlugt [16] Configurational-Bias Monte Carlo Thijs J.H. Vlugt [17] Monte Carlo Casino Smart Monte Carlo: α ( o ! n ) 6 = α ( n ! o ) Not a random displacement ∆ r uniformly from [ � δ , δ ] , but instead ∆ r = r (new) � r (old) = A ⇥ F + δ r F : force on particle A : constant δ r : taken from Gaussian distribution with width 2 A so P ( δ r ) ⇠ exp[ � ( δ r 2 ) / 4 A ] � ( r new � ( r old + A ⇥ F ( o ))) 2  � P ( r new ) ⇠ exp 4 A h � ( ∆ r � A ⇥ F ( o )) 2 i exp α ( o ! n ) 4 A α ( n ! o ) = h i � ( ∆ r + A ⇥ F ( n )) 2 exp 4 A

Configurational-Bias Monte Carlo Thijs J.H. Vlugt [18] Configurational-Bias Monte Carlo Thijs J.H. Vlugt [19] Chain Molecules Self-Avoiding Walk on a Cubic Lattice • 3D lattice; 6 lattice directions • only 1 monomer per lattice site (otherwise U = 1 ) • interactions only when | r ij | = 1 and | i � j | > 1

Configurational-Bias Monte Carlo Thijs J.H. Vlugt [20] Configurational-Bias Monte Carlo Thijs J.H. Vlugt [21] Simple Model for Protein Folding Number of Configurations without Overlap Random Chains: 20 by 20 interaction matrix ∆ ij P n i =1 R i exp[ � β U i ] h R i = lim P n i =1 exp[ � β U i ] n !1 YPDLTKWHAMEAGKIRFSVPDACLNGEGIRQVTLSN Fraction of chains without overlap decreases exponentially as a function of (E. Jarkova, T.J.H. Vlugt, N.K. Lee, J. Chem. Phys., 2005, 122, 114904) chainlength ( N ) total ( = 6 N � 1 ) N without overlap fraction no overlap 2 6 6 1 20 6 b) 6 7776 3534 0.454 5 8 279936 81390 0.290 15 4 10 10077696 1853886 0.183 z (b) 2 (dz) 10 3 12 362797056 41934150 0.115 prot1 13 2176782336 198842742 0.091 2 prot2 5 14 13060694016 943974510 0.072 prot3 1 prot4 15 78364164096 4468911678 0.057 0 0 0 2 4 6 8 10 0 2 4 6 8 10 16 470184984576 21175146054 0.045 Force (k B T/b) Force (k B T/b) 1 . 3 ⇥ 10 � 5 50 · · · · · ·

Configurational-Bias Monte Carlo Thijs J.H. Vlugt [22] Configurational-Bias Monte Carlo Thijs J.H. Vlugt [23] Rosenbluth Sampling (1) Rosenbluth Sampling (2) exp[ � β u ij ? ] P j ? = 1. Place first monomer at a random position P k j =1 exp[ � β u ij ] 2. For the next monomer ( i ), generate k trial directions ( j = 1 , 2 , · · · , k ) each with energy u ij 3. Select trial direction j ? with a probability exp[ � β u ij ? ] P j ? = P k j =1 exp[ � β u ij ] 4. Continue with step 2 until the complete chain is grown ( N monomers)

Configurational-Bias Monte Carlo Thijs J.H. Vlugt [24] Configurational-Bias Monte Carlo Thijs J.H. Vlugt [25] Rosenbluth Sampling (3) Rosenbluth Sampling (4) Probability to grow this chain ( N monomers, k trial directions) Probability to choose trial direction j ? for the i th monomer Q N i =1 exp[ � β u ij ? ( i ) ] = exp[ � β U chain ] P chain = exp[ � β u ij ? ] Q N P k W chain j =1 exp[ � β u ij ] P j ? = i =1 P k j =1 exp[ � β u ij ] Therefore, weightfactor for each chain i is the Rosenbluth factor W i : Probability to grow this chain ( N monomers, k trial directions) P n i =1 W i ⇥ R i h R i Boltzmann = lim P n N Q N i =1 W i i =1 exp[ � β u ij ? ( i ) ] = exp[ � β U chain ] n !1 Y P chain = P j ? ( i ) = Q N P k W chain j =1 exp[ � β u ij ] i =1 i =1 The unweighted distribution is called the Rosenbluth distribution: P n i =1 R i h R i Rosenbluth = lim n n !1 Of course: h R i Rosenbluth 6 = h R i Boltzmann

Configurational-Bias Monte Carlo Thijs J.H. Vlugt [26] Configurational-Bias Monte Carlo Thijs J.H. Vlugt [27] Intermezzo: Ensemble Averages at Di ff erent Temperatures Rosenbluth Distribution Di ff ers from Boltzmann Distribution Probability distribution for the end-to-end distance r Ensemble averages at β ? can (in principle) be computed from simulations at β : R d r N U ( r N ) exp[ � β U ( r N )] N = 10 N = 100 h U i � = d r N exp[ � β U ( r N )] R 0.4 0.1 d r N U ( r N ) exp[ � β ? U ( r N )] exp[( β ? � β ) ⇥ U ( r N )] R Rosenbluth distribution Boltzmann distribution = Rosenbluth distribution d r N exp[ � β ? U ( r N )] exp[( β ? � β ) ⇥ U ( r N )] R Boltzmann distribution 0.08 0.3 U ( r N ) exp[( β ? � β ) ⇥ U ( r N )] ⌦ ↵ � ? = 0.06 h exp[( β ? � β ) ⇥ U ( r N )] i � ? P(r) P(r) 0.2 ⌦ U ( r N ) exp[ ∆ β ⇥ U ( r N )] ↵ 0.04 � ? = h exp[ ∆ β ⇥ U ( r N )] i � ? 0.1 0.02 Useful or not??? 0 0 0 2 4 6 0 5 10 15 20 25 30 r r