Transcript matut 5431
An Introduction to Monte Carlo Methods in Statistical Physics Kristen A. Fichthorn The Pennsylvania State University University Park, PA 16803 Monte Carlo Methods: A New Way to Solve Integrals (in the 1950’s) “Hit or Miss” Method: What is p? Algorithm: C 1 B y •Generate uniform, random x and y between 0 and 1 •Calculate the distance d from the origin d x2 y 2 •If d ≤ 1, thit = thit + 1 •Repeat for ttot trials 0 x 1 A 4 x Area Under Curve CA p Area of Square OABC 4t hit t tot Monte Carlo Sample Mean Integration To Solve: x2 F dx f ( x) x1 We Write: Then: f ( x) ( x) F dx (x) x1 x2 f ( ) F ( ) trials When on Each Trial We Randomly Choose from Monte Carlo Sample Mean Integration: Uniform Sampling to Estimate p To Estimate F π 2 x 2 1 dx (1 x 2 1 / 2 ) x1 0 Using a Uniform Distribution 1 ( x) x2 x1 , x1 x x2 Generate ttot Uniform, Random Numbers p 1 t 2( x 2 x1) t tot t tot t 1 2 1 / 2 , 0 t 1 Monte Carlo Sample Mean Integration in Statistical Physics: Uniform Sampling Z NVT dr exp U (r ) L Quadrature e.g., with N=100 Molecules 3N=300 Coordinates 10 Points per Coordinate to Span (-L/2,L/2) 10300 Integration Points!!!! L L Uniform Sample Mean Integration •Generate 300 uniform random coordinates in (-L/2,L/2) •Calculate U V N t tot Z NVT exp U (t ) •Repeat ttot times… t tot t 1 Problems with Uniform Sampling… Z NVT N t tot V t tot exp U (t ) t 1 Too Many Configurations Where exp U (r ) 0 Especially for a Dense Fluid!! L L L What is the Average Depth of the Nile? Integration Using… t max d d t wt t 1 t max wt t 1 1 , if in the Nile w(t ) 0 , else Quadrature or Uniform Sampling vs. Importance Sampling Adapted from Frenkel and Smit, “Understanding Molecular Simulation”, Academic Press (2002). Importance Sampling for Ensemble Averages If We Want to Estimate an Ensemble Average Efficiently… A A NVT NVT NVT A A trials trials A NVT dr NVT (r ) A(r ) exp U (r ) NVT (r ) Z NVT We Just Need to Sample It With NVT !! Importance Sampling: Monte Carlo as a Solution to the Master Equation dP(r , t ) p (r r ' ) P(r , t ) p (r ' r ) P(r ' , t ) dt r' r' P (r , t ) : Probability to be at State r at Time t : Transition Probability per Unit Time p (r r ' ) from r r to r' r' The Detailed Balance Criterion After a Long Time, the System Reaches Equilibrium dP(r , ) 0 p ( r r ' ) P ( r , ) p ( r ' r ) P ( r ' , ) dt x x At Equilibrium, We Have: P(r , ) NVT (r ); P(r ' , ) NVT (r ' ) This Will Occur if the Transition Probabilities p Satisfy Detailed Balance p (r ' r ) NVT (r ' ) p (r r ' ) NVT (r ) p (r r ' ) NVT (r ' ) exp U (r ' ) U (r ) p (r ' r ) NVT (r ) Metropolis Monte Carlo Let p Take the Form: p (r r ' ) (r r ' ) acc(r r ' ) = Probability to Choose a Particular Move acc = Probability to Accept the Move NVT (r ' ) , if NVT (r ' ) NVT (r ) α(r r ') Use: NVT (r ) p (r r ' ) α(r r ') , if NVT (r ' ) NVT (r ) (r r ' ) (r ' r ) 1/ N With: N. Metropolis et al. J. Chem. Phys. 21, 1087 (1953). Metropolis Monte Carlo Use: 1 N exp (U 'U ) , if U U ' p (r r ' ) 1 , if U U ' N Detailed Balance is Satisfied: p (r r ' ) exp (U 'U ) p (r ' r ) Metropolis MC Algorithm Initialize the Positions Atot 0; t tot 0 Calculate the Ensemble Average A A Select a Particle at Random, U ( r ) U Calculate the Energy Atot Atot A(r ) t tot t tot 1 Give the Particle a Random Displacement, Calculate the New Energy U (r ' ) U ' tot t tot Yes No Finished ? Accept the Move with 1 p (r r ' ) min exp U 'U Periodic Boundary Conditions L d L If d>L/2 then d=L-d It’s Like Doing a Simulation on a Torus! Nearest-Neighbor, Square Lattice Gas A B Interactions eAA 0.0 -1.0kT eBB 0.0 0.0 eAB -1.0kT 0.0 When Is Enough Enough? Run it Long... 100 200 Energy 300 …and Longer! 400 500 600 700 800 0 100000 200000 Trials 300000 400000 When Is Enough Enough? 0.3 0.1 0.4 0.2 Energy Energy 0.5 0.6 0.7 0.3 0.4 0.5 0.8 0.6 0.9 0.7 1 0 2500 5000 7500 10000 Trials 12500 15000 Run it Big… 0 100000 300000 …and Bigger! t tot Estimate the Error 200000 Trials 1/ 2 xt x Error in x t 1 t tot t tot 1 2 400000 When Is Enough Enough? Make a Picture! When Is Enough Enough? Try Different Initial Conditions! Phase Behavior in Two-Dimensional Rod and Disk Systems TMV and spheres E. coli Bottom-up assembly of spheres Nature 393, 349 (1998). Electronic circuits Use MC Simulation to Understand the Phase Behavior of Hard Rod and Disk Systems Lamellar Nematic Smectic Miscible Nematic Isotropic Miscible Isotropic Hard Systems: It’s All About Entropy A = U – TS Hard Core Interactions U = 0 if particles do not overlap U = ∞ if particles do overlap Maximize Entropy to Minimize Helmholtz Free Energy Depletion Zones Overlap Volume Ordering Can Increase Entropy! Metropolis Monte Carlo Old Configuration U old (r ) 0 Perform Move at Random Pold new exp U new U old exp 0 New Configuration Ouch! A Lot of Infeasible Trials! U new (r ' ) Small Moves or… Configurational Bias Monte Carlo Old Rosenbluth & Rosenbluth, J. Chem. Phys. 23, 356 (1955). k W(old) exp βU or(b j ) j 1 Move Center of Mass Randomly Generate k-1 New Orientations bj New Final k W (new) exp U or (b j ) j 1 Select a New Configuration with exp U or (bn ) P(bn ) W (new) Accept the New Configuration with W (new) Pold new min 1, W (old ) Configurational Bias Monte Carlo and Detailed Balance Recall we Have of the Form: p p (o n) (o n) acc(o n) exp U (n) P (bn ) (o n) W ( n) exp U (o) ( n o) W (o ) The Probability of Choosing a Move: The Acceptance Ratio: Detailed Balance acc(o n) W (n) acc(n o) W (o) p (o n) exp U (n) U (o) p (n o) Properties of Interest Nematic Order Parameter Q 1 N 2u (i)u (i) N i 1 Radial Distribution Function A g r 2 N 2 r ri r j N N i 1 j 1 i j Orientational Correlation Functions g 2 r cos2 0 r g 4 r cos4 0 r Snapshots 800 rods ρ = 3.5 L-2 1257 rods ρ = 5.5 L-2 Snapshots 6213 rods ρ = 6.75 L-2 8055 rods ρ = 8.75 L-2 Energy Accelerating Monte Carlo Sampling x How Can We Overcome the High Free-Energy Barriers to Sample Everything? Accelerating Monte Carlo Sampling: Parallel Tempering E. Marinari and and G. Parisi, Europhys. Lett. 19, 451 (1992). System N at TN … System 3 at T3 System 2 at T2 System 1 at T1 TN >…>T3 >T2 >T1 Metropolis Monte Carlo Trials Within Each System Swaps Between Systems i and j 1 P(o n) min exp ( )( U U ) i j j i Parallel Tempering in a Model Potential x 2 1 (1 sin( 2px)) 2 x 1.25 2 (1 sin( 2px)) 1.25 x 0.25 U ( x) 3 (1 sin( 2px)) 0.25 x 0.75 4 (1 sin( 2px)) 0.75 x 1.75 1.75 x 2 5 (1 sin( 2px)) x2 90% Move Attempts within Systems 10% Move Attempts are Swaps Adapted from: F. Falcioni and M. Deem, J. Chem. Phys. 110, 1754 (1999). System 3 at kT3=5.0 System 2 at kT2=0.5 System 1 at kT1=0.05 Good Sources on Monte Carlo: D. Frenkel and B. Smit, “Understanding Molecular Simulation”, 2nd Ed., Academic Press (2002). M. Allen and D. J. Tildesley, “Computer Simulation of Liquids”, Oxford (1987).