Transcript Monte Carlo Methods - uni

Monte Carlo Methods
H. Rieger, Saarland University, Saarbrücken, Germany
Summerschool on Computational Statistical Physics, 4.-11.8.2010
NCCU Taipei, Taiwan
Monte Carlo Methods =
Stochastic evaluation of physical quantities using random numbers
Example: Pebble game on the beach of Monaco
= computation of  using random events
Kid‘s game:
„Direct sampling“
Adult‘s game:
„Markov chain sampling“
[after W. Krauth, Statistical Mechanics, Algorithms and Computation, Oxford Univ. Press]
Direct sampling
easy
„hard spheres in 2d“ – hard!
No direct sampling algorithm for hard spheres
(NOT random sequantial adsorption!)
Markov chain sampling
What to do,
when a stone from the lady
lands here?
1) Simply move on
2) Climb over the fence, and continue
until, by accident, she will reenter
the heliport
3) ….?
Move a pebble in a 3x3 lattice probabilistically
such that each site a is visited with the same probability p(a) = 1/9
w(ab) = ¼
ba
w(aa) = 0
w(aa) = 1/4
w(aa) = 1/2
Detailled Balance
b
a
The following condition for p(a) must hold:
c
Together with
This yields:
This condition is fulfilled when
etc.
This is called „detailled balance condition“
– and leads here to w(ab,c)=1/4, w(aa)=1/2
Rejections
Adult‘s pebbel game – solution: if a stone lands outside the heliport, stay there
put a stone on top of the present stone and continue,
i.e. reject move outside the square!

More or less large piles close to the boundary due to „rejections“
Master equation
Markov chain described by „Master equation“ (t = time step)
The time independent probability distribution p(a) is a solution of this equation,
if the transition probabilities w(ab) etc. fulfill
detailed balance:
For a given p(a) one can, for instance, choose
Which is the „Metropolis“ rule
Monte Carlo for Thermodynamic Equilibrium
a are configurations of a many particle system,
E(a) the energy of configuration a.
Thermodynamic equilibrium at temperature T is then described by the
Boltzman distribution
 = 1/kBT inverse temperature
is the normalization,
called partition function
Thus the Metropolis rule for setting up a Markov chain leading to
The Boltzmann distribution is
Hard spheres
a = (r1,r2, … , rN)
- configurations = coordinates of all N spheres in d dimension
All spheres have radius , they are not allowed to overlap,
otherwise all configurations have the same energy (no interactions):

H(a) =  if there is a pair (i,j) with |ri-rj|<2
H(a) = 0 otherwise
i.e.:
Define w(ab) in the following way:
In configuration a choose randomly a particle i
and displace it by a random vector  - this constitutes configuration b.
w(ab) = 1 if b is allowed (no overlaps),
w(ab) = 0 (reject) if displaced particle overlaps with some other particle
Hard spheres (2)
Tagged particle
Iteration: t
Tagged particle
t
Soft spheres / interacting particles
a = {(r1,p1),(r2,p2),…,(rN,pN)}
- configurations = coordinates and momenta of all N particles in d dimension
Energy:
Partition function
Peforming the
momentum integral
(Gaussian)
 Left with the
configuration integral I
Example:
LJ (Lennard-Jones)
L = box size
MC simulation for soft spheres: Metropolis
Choose randomly particle i, its position is ri
Define new position by ri‘=ri+,  a random displacement vector, [-,]3
All othe rparticle remain fixed.
Acceptance probability for the new postion:
if
otherwise
Repeat many times
Measurements:
Energy, specific heat, spatial correlation functions, structure function
Equlibration!
Note: Gives the same results
as molecular dynamics
Discrete systems: Ising spins
System of N interacting Ising spins Si  {+1,-1},
placed on the nodes of a d-dimensional lattice
a = (S1,S2,…,SN): spin configurations
Energy:
(i,j)
Jij = coupling strengths, e.g. Jij = J > 0 for all (i,j)  ferromagnet
h = external field strength
For instance 1d:
with periodic
poundary conditions
Quantities of interest / Measurements
Magnetization
Susceptibility
Average energy
Specific heat
How to compute:
where at are the configurations generated by the
Markov chain (the simulation) at time step t.
Ising spins: Metropolis update
H(S,S‘) = H(S‘)-H(S)
for
for
Procedure Ising Metropolis:
Initialize S = (S1,…,SN)
label Generate new configuration S‘
Calculate H= H(S,S‘)
if H  0
accept S‘ (i.e. S‘S)
else
generate random numer x[0,1]
if x<exp(-H)
accept S‘ (i.e. S‘S)
compute O(S)
goto label
Single spin flip Metropolis for 2d Ising
Procedure single spin flip
Input L, T, N=L*L
Define arrays: S[i], i=1,…,N, h[i], i=1,…,N, etc.
Initialize S[i], nxm[i], nxp[i],…., h[i]
step = 0
while (step<max_step)
choose random site i
calculate dE = 2*h[i]*S[i]
if ( dE <= 0 )
S[i]=-S[i]; update h[nxm[i]], h[nxp[i]],…
else
p = exp(-dE/T)
x = rand()
if ( x<p)
S[i]=-S[i]; update h[nxm[i]], h[nxp[i]],…
compute M(S), E(S), …
accumulate M, E, …
step++
Output m, e, …
Implementation issues
Periodic boundary conditions
Neighbor tables
e.g.:
if
if
Implementation issues (2)
With single spin flip E(S) and E(S‘) differ only by 4 terms in 2d (6 terms in 3d):
Flip spin i means Si‘ = -Si, all other spins fixed, i.e. Sj‘=Sj for all ji


Tabulate exponentials exp(-4), exp(-8)
to avoid transcendental functions in the innermost loop
Use array h[i] for local fields, if move (flip is rejected nothing to be done,
if flip accepted update Si and hnxm[i], hnxp[i], etc.
Study of phase transitions with MC
Ising model in d>1 has a 2nd order phase transition at h=0, T=Tc
Magnetization (order parameter):
h=0
2nd order
phase transition
as a function of T
at h=0!
Phase
diagram
1st order
phase transition
as a function
of h!
h0
h
T<Tc
m
Critical behavior
Singularities at Tc in the thermodynamic limit (N):
Magnetization:
Susceptibility:
Specific heat:
Correlation function:
Correlation length:
Scaling relations:
Finite size behavior
w. Periodic b.c.
Finite Size Scaling
FSS forms:
4th order cumulant:
Dimensionless (no L-dependent prefactor)
- Good for the localization of the critical point
Critical exponents of the d-dim. Ising model
Slowing down at the critical point
Quality of the MC estimats of therodynamic expectation values depends
on the number of uncorrelated configurations –
Need an estimate for the correlation time  of the Markov process!
Autocorrelationfunction
Schematically:
 for TTc
Configurations should decorrelate
faster than with single spin-flip!
Solution: Cluster Moves
Cluster Algorithms
Construction process of the clusters in the Wolff algorith:
Start from an initial + site, include other + sites with prbability p (left).
The whole cluster (gray) is then flipped
Here c1=10, c2=14
„A priori“ or construction probability:
p(a)
p(b)
Detailled balance condition: p(a) A(ab) w(ab) = p(b) A(ba) w(ba)

Wolff algorithm (cont.)
Once the cluster is constructed with given p, one gets c1 and c2 ,
with which one can compute the acceptance probability w(ab)
But with p = 1-e-2 the acceptance probability w(ab) becomes 1!
Thus with p=1-e-2 the constructed cluster is always flipped!
Remarkable speed up, no critical slowing down at the critical point!
Wolff cluster flipping for Ising
(1)
Randomly choose a site i
(2)
Draw bonds to all nearest neighbors j
with probability
(3)
If bonds have been drawn to any site j draw bonds to all
nearest neighbors k of j with probability
(4)
Repeat step (3) until no more bonds are created
(5)
Flip all spins in the cluster
(6)
Got to (1)
(Note
for S=S‘, and = 0 for SS‘, such that p=0 for SjSk)
Swendsen-Wang algorithm
Similar to Wolff, but
(1)
Draw bonds between ALL nearest neighbors
with probability
(2)
Identify connected clusters
(3)
Flip each individual cluster with probability 1/2