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