The MonteCarlo method
As we know, it is impossible to calculate statistical averages directly.
One could generate spins configurations at random and approximate the real
thermal averages by MonteCarlo averages. The next problem here is, most
randomly chosen states will make a negligible contribution to the sum.
Random walks are used to take into account only important spins configurations
by the introduction of a fictitious dynamics (socalled Markov chain of
configurations). This wandering generates states, which are the most probable
from energyentropy point of view (importance sampling).
The Metropolis algorithm
In the Metropolis algorithm we try to turn over a single spin direction with
transition probability
W_{12} =
exp[(E_{1}E_{2})/T]

if E_{1} < E_{2}

W_{12} = 1

if E_{1} > E_{2}

where E_{1}, E_{2} are energies of the old and new
configurations (see details in the Gould and Tobochnik book).
Statistical averages may then be computed as simple arithmetic means.
The Metropolis algorithm can by applied easy e.g. to the XY model simulation.
There are many possible choices for W_{12}.
To get equilibrium final spins distribution it is enough that
W_{ij} satisfy the detailed balance conditions
W_{12} exp(E_{1}/T) =
W_{21} exp(E_{2}/T) or
W_{12} / W_{21} =
exp[(E_{1}E_{2})/T].
This Java applet is based on the Gould and Tobochnik book.
A successive loop over all sites (referred as one Metropolis sweep or iteration)
is used. Then "momentary" (in Markov's time) E_{t}, M_{t}
values are calculated
<M>_{M} = 1/T
∑ _{t=1,T} M_{t} ,
<E>_{M} = 1/T
∑ _{t=1,T} E_{t} .
Susceptibility χ and specific heat C
of the system can be found by the thermal fluctuations of magnetization
M and energy E
C = d<E>/dT = n/T^{2}
(<E^{2}>  <E>^{2}) ,
χ = lim_{H→ 0} dM/dT
= n/T (<M^{2}>  <M>^{2}) .
Controls Click mouse into lattice to get a new spins configuration.
Drag by mouse "thermometer" to the right (the black bar corresponds to
T_{c}) or press "Enter" to set a new T value from
the text field.
High nonequilibrium spin flip dynamics
You see below two features in the Metropolis algorithm iterations.
First of all, as since transition probability W_{1>2} =
exp[(E_{1}E_{2})/T] > 1 for high
temperature T >> (E_{1}  E_{2}) ~ J , therefore
spins oscillate in big solid clusters. Second, you see growth of dendrite
structures in the leftbottom part of the pictures,
because all spins are flipped by lines from this corner
(to avoid this feature one could choose flipped spins at random).
The thermostat algorithm
In the thermostat algorithm we get a single spin (all the rest are
fixed) into contact with big thermostat at temperature T. Then we
get for transition probabilities (the Glauber formula)
P_{+} = exp(E_{+} /T) /
[exp(E_{+} /T) + exp(E_{} /T)] =
exp((E_{+}E_{} )/T) /
[exp((E_{+}E_{} )/T) + 1] ,
P_{} = exp(E_{} /T) /
[exp(E_{+} /T) + exp(E_{} /T)] = 1  P_{+} .
The Glauber transition probabilities satisfy the detailed balance conditions
too. For high temperature P_{+} , P_{} → 1/2
and it leads to fast and smooth "cluster melting" (see the Ising applet with
the thermostat algorithm below).
But for low temperature (e.g. T = 1) and random initial
configurations you can see formation of metastable clusters. Then they are
decreasing slowly and one single cluster will appear. In the Metropolis method
"freezing" is faster, but sometimes we get two clusters with a big metastable
domain wall.
Due to ambiguity of choice of transition probabilities this dynamics
seems me a bit artificial, but it influences on spins relaxation greatly.
G.T.Barkema, M.E.J.Newman New Monte Carlo algorithms for classical
spin systems arXiv:condmat/9703179
See also Kenji Harada's
Java
applet demonstrating Metropolis, Swendsen and Wang, Wolff algorithms
for Ising model.
Contents
Previous: 2D Ising model
Next: Relaxation time
updated 27 Dec 2001