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

Key Point

• Metropolis algorithm

• A useful method to generate an importance Markov chain in Monte Carlo sampling.
• Applications
• (1) generation of Gaussian distribution
• (2) Ising model for magnets

Metropolis algorithm

An importance Markov chain Monte Carlo sampling for a specific distribution P(X) goes as follows:

• An initial state (configuration): X0
• A sequence of states generated by a Markov chain through the transition probability W(Xi to Xi+1):
X0 => X1 => ... => Xi => ... => X
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):

• It satisfies the detailed balance condition.
• Easy to generate.
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).

Steps:

• Specify an initial state X0, set i = 0.
• Generate a new state X'i+1 according to T(Xi to X'i+1).
• Accept the new state with probability: min(1,P(X'i+1)/P(Xi)).
• That is, Xi+1 = X'i+1 if a random number between 0 and 1 is less than P(X'i+1)/P(Xi); retain the old state as (i+1)-th state, Xi+1 = Xi, otherwise.
• set i= i + 1, go to step 2.

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

• Starting with some arbitrary real value x

•
• we then propose a new value as
x' = x + alpha (xi - 1/2)
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.
This choice of x' is equivalent to take a particular functional form of T(x to x'). As long as T is symmetric, the exact form of T does not matter as far as generating correct distribution is concerned.

• The value x' will be accepted as the next point if a uniformly distributed random number between 0 and 1 is less than:
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:

• spontaneous magnetization at lower temperatures.
• the magnetization disappears at some higher temperature Tc.
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:

• Internal energy per spin

u = 1/Ld < H >

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

• Magnetization
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.

• Specific heat

c = du/dT

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

• Susceptibility
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

M = | sum sigmai |
is the absolute value of the total magnetization of a particular state.

• The fourth-order cumulant
g =  1/2 [ 3 - < M4 >/ < M2 >2 ]

which is useful in determining the critical temperature Tc.

Monte Carlo simulation of the Ising model:

Steps:

• Initialize sigmai at each site with arbitrary spin values, e.g., all spin up, or up or down with equal probability.
The complete set of spins { sigma } is the abstract state X discussed in the previous sections.

• Choose a site i at random and propose to flip the spin at that site, sigma'i = - sigmai.
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

• Calculate the energy increment
E'-E = H( { sigma' } ) - H( { sigma } )

=   2J sigmai sum neighbors of i sigmaj
+ 2B sigmai

• Accept the proposed state as the new state if a uniformly distributed random number (between 0 and 1) is less than exp[- (E'-E)/kB T]; retain the old state as the new state otherwise.
• go to step 2.

Homework