Computational Techniques in Theoretical Physics
Section 12
Metropolis algorithm and application examples


 Key Point
Metropolis algorithm
An importance Markov chain Monte Carlo sampling for a specific distribution P(X) goes as follows:
  The starting state may be a unique state or may be generated at random with certain distribution P0.

The new state Xi+1 is based on the old state Xi. The new state is generated in such a way so that the probability that the system is in state Xi+1 is just W(Xi to Xi+1)


Requirement on W(Xi to Xi+1):

W(Xi to Xi+1) from Metropolis algorithm:
W(X -> X') = T(X -> X') min(1, P(X')/P(X))
where X =/= X', T(X -> X') = T(X' -> X)

The transition matrix T can be any distribution but must be symmetric.

Normally, it is chosen as:

                           const.    for X' inside a region around X
T(X -> X') = {
                           0            otherwise

The probability that the state does not change, W(X to X), is determined by the normalization condition:

sumX' W(X to X') = 1

Complete description of Metropolis algorithm:


Generate a sequence X0, X1, X2, ... , Xi, ... , with a probability distribution P(Xi).



Many physical quantities can be computed as expectation values.

Expectation values (or averages) can be estimated by:

< A(X) > = sumX A(X) P(X)
                ~ N-1 sum i = M+1 to M+N A(Xi)

It is important that we throw away the first M points in the average. We have to wait for the system approaching equilibrium before we use the points.

The reason for this is that the system is distributed according to Pk(X), and we want the distribution P(X). Pk(X) will converge to P(X) for sufficiently large k.

This convergence is usually fast. For efficiency consideration, usually we do not use every Xi even if we have reached equilibrium, since they are highly correlated. We calculate A(Xi) only every certain number of iterations or steps.


Examples of Metropolis Algorithm

Gaussian distribution
where xi is a uniformly distributed random number between 0 and 1; alpha is some constant, say, 0.3.
The newly proposed value x' is distributed uniformly inside an interval of length alpha, centered at the original value x.
exp{ - 1/2 (x'2 - x2)}
otherwise, the old value x is taken as the next point.
This generates a sequence of points, which should be distributed according to Gaussian with zero mean and unit variance.

Ising Model for Magnets

What is Ising model?

Ising model is the simplest model for ferromagnets.

Important features of a magnet:

The phase below Tc is called ferromagnetic phase and that above Tc is called paramagnetic phase. Tc is called critical temperature and there is a phase transition at this temperature.
Two-dimensional system as an example:

On a L x L square lattice, we put a spin sigmai at each lattice site i.

The spin here is a `classical' spin, not the quantum-mechanical operator. The spin can be thought as the z-component of the Pauli spin, taken only two values +1 and -1.

For each configuration denoted by:
{sigma} = (sigma1, sigma2, sigma3, ...)
which is the collection of all the spins, the system has a definite energy given by:

H(sigma) = - J sum <i, j> sigmai sigmaj
                    - B sum i sigmai

where J is called coupling constant; B is proportional to the external magnetic field.

The first summation runs over the nearest neighbors only, i.e., each bond is summed once.

Such a model is also referred as nearest neighbor Ising model.
The second summation is over all the lattice sites.

When J is positive, the model describes a ferromagnet; while J < 0 describes an anti-ferromagnet.

Let's consider just one term in the nearest neighbor interaction. The interaction energy of neighbors takes only two values, +J or -J. For ferromagnet, the energy is lower if the two spins are parallel. It costs energy if the spins are pointing in different directions.
Statistical-mechanical description of Ising model:
Due to thermal fluctuation, the system will have some probability distribution in states other than equilibrium one.

If the temperature T is a fixed parameter, we have the famous Boltzmann distribution (Canonical distribution):

P( { sigma } ) ~ exp [ -H( { sigma } ) / kB T ]
where k~ 1.38  x 10-23 Joule/Kelvin is called the Boltzmann constant.

Our goal is to calculate average quantities such as energy, magnetization, correlation functions, etc.

Definition of physical quantities:

u = 1/Ld < H >

   =  1/Ld sum { sigma } H( { sigma } )

m = 1/Ld < M > =
    =  1/Ld  sum { sigma }   | sum sigmai | P( { sigma } )
The summation is over all possible states.

Since each site can have two states (spin up and spin  down), we have 2Ld states for a system of linear size L in d dimensions.

c = du/dT

   = [ < H2 > - < H > 2 ] / kBT2 Ld

chi = dm/dB
               A ( < M2 > - < M >2 )   T < Tc
       =  {
               A < M2 >                          T > Tc
where A = 1/ kBT Ld , the angular brackets < >  denote average over the Boltzmann distribution and
g =  1/2 [ 3 - < M4 >/ < M2 >2 ]
which is useful in determining the critical temperature Tc.
Monte Carlo simulation of the Ising model:


The complete set of spins { sigma } is the abstract state X discussed in the previous sections.
The proposed state X' is the one with one spin at location i flipped. This amounts to take:
                                      1/Ld, if X and X' differ by 1 spin
             T( X to X') = {
                                      0,      otherwise