The Monte-Carlo 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 Monte-Carlo 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 (so-called Markov chain of configurations). This wandering generates states, which are the most probable from energy-entropy point of view (importance sampling).

The Metropolis algorithm

In the Metropolis algorithm we try to turn over a single spin direction with transition probability
    W12 = exp[(E1-E2)/T]     if   E1 < E2
    W12 = 1 if   E1 > E2
where E1, E2 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 W12. To get equilibrium final spins distribution it is enough that Wij satisfy the detailed balance conditions
    W12 exp(-E1/T) = W21 exp(-E2/T) or
    W12 / W21 = exp[(E1-E2)/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) Et, Mt values are calculated
    <M>M = 1/T ∑ t=1,T Mt ,     <E>M = 1/T ∑ t=1,T Et .
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/T2 (<E2> - <E>2) ,
    χ = limH→ 0 dM/dT = n/T (<M2> - <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 Tc) or press "Enter" to set a new T value from the text field.

High non-equilibrium spin flip dynamics

You see below two features in the Metropolis algorithm iterations. First of all, as since transition probability W1->2 = exp[(E1-E2)/T] -> 1 for high temperature T >> (E1 - E2) ~ J , therefore spins oscillate in big solid clusters. Second, you see growth of dendrite structures in the left-bottom 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:cond-mat/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