• Applications

• (1) Ising model for magnets
• (2) Percolation simulation based on Ising model
• (3) Other applications: simple fluid and polymer

C-code:
 

The following is an implementation of the two-dimensional Ising model without the magnetic field (B=0).

This function performs one Monte Carlo step by trying to flip spins L² times. It picks a site at random and tries to flip it. The flip is actually performed if the energy change is negative or zero, or if the random number is less than the probability exp{-(E'-E)/kBT} when the energy change is positive.
 

• N is the total number of spins,
• Z is coordination number.
• Spin s[ ] and temperature are passed globally.

#include  <math.h>                              exponential function and
                                                                           random number need this
                                                                           macro definitions

#define L 16                                          lattice linear size
#define N (L*L)                                    total number of spins
#define Z 4                                            coordination number

int s[N];                                                 spin +1 or -1
double T = 2.26918531;                        dimensionless temperature

void monte_carlo()
{
  int i, j, k, e;                                          i is the center site
  int nn[Z];                                             the names of neighbors

  for(k = 0; k < N; ++k) {                       do one Monte Carlo step
    i = drand48() * (double) N;               pick the center site at random
    neighbor(i, nn);                                 find neighbors of site i
    for(e = 0, j = 0; j < Z; ++j)                 go over the neighbors
      e += s[nn[j]];                                   sum of the neighbor spins
      e *= 2 * s[i];                                    2 times the center spin
      if (e <= 0)                                        when energy change is less than
        s[i] = - s[i];                                   or equal to zero, spin is flipped
      else if (drand48() < exp(-e/T))      otherwise, it is flipped
        s[i] = - s[i];                                   with probability less than one
    }
 }

The neighbor() function is the same as in the previous section for the cluster labelling algorithm.
 
 
 

 

Swendsen-Wang multi-cluster algorithm
 

One starts with a spin configuration { sigma } and generates a percolation configuration based on the spin configuration by the rule described below. Then the old spin configuration is thrown away and a new spin configuration { sigma'} is generated based on the percolation configuration.

The rule is such that the detailed balance is satisfied. Or more generally, the transition leaves the equilibrium probability invariant.

The following is the algorithm for nearest neighbor Ising model without magnetic field (B=0). It can be generalized to other models and with magnetic field.

 
 

Swendsen-Wang Algorithm:
 

• Start with some arbitrary state {sigma }.

• Go through each nearest neighbor connection of the lattice, create a bond between the two neighboring sites i and j with probability 
               p = 1 - exp{-2J/(k_BT)}
  only if the spins are the same, sigma_i = sigma_j.

Never put a bond between the sites if the spin values are different. We are creating a bond percolation configuration with probability p on a subset of lattice sites where the spins are either all pointing up or down.

• Identify clusters as a set of sites connected by bonds, or isolated sites. Two sites are said to be in the same cluster if there is a connected path of bonds joining them. Every site has to belong to one of the clusters.

• After the clusters are found, each cluster is assigned a new Ising spin chosen with equal probability between +1 and -1. The old spin values now can be `forgotten'. And all the sites in a cluster now take the value of spin assigned to the cluster.

• Repeat from step 2 for the next iteration.

The performance of the algorithm in terms of correlation time in comparison with Metropolis algorithm is remarkable. The dynamic critical exponent z ~ 2 for Metropolis algorithm is almost independent of the dimensionality.

The Swendsen-Wang algorithm in one dimension gives z = 0 (tau approaches a constant as T -> 0). For the two-dimensional Ising model, the dynamic critical exponent z is less than 0.3 or possibly zero (but tau -> ln L). In three dimensions it is about 0.5. At and above four dimensions it is 1.
 
 
 

Wolff single-cluster algorithm

In the Swendsen-Wang algorithm, we generated many clusters and then flip these clusters. In Wolff algorithm one picks a site at random, and then generates one single cluster by growing a cluster from the seed.

The neighbors of the seed site will also belong to the cluster if the spins are in parallel and a random number is less than p = 1 - exp{-2J/(kB T)}.

That is, the neighboring site will be in the same cluster as the seed site with probability p if the spins have the same sign.

If the spins are different, neighboring site will never belong to the cluster. Neighbors of each new site in the cluster are tested for membership. This testing of membership is performed on pair of sites (forming a nearest neighbor bond) not more than once.

The recursive process will eventually terminate. The spins in the cluster are turned over with probability 1.

Single-cluster algorithm appears more efficient than the multi-cluster one.

The following is a fairly efficient way of implementing the Wolff algorithm. The monte_carlo_step(n) function does the Wolff single-cluster flip n times.

#include <math.h>           drand48() need this macro definitions

#define L 16                      lattice linear size
#define N (L*L)                total number of spins
#define Z 4                        coordination number
                                                    = 2 x d global variables
int s[N];                             spin +1 or -1
double p = 0.5;                   percolation probability
                                                    p=1-e^{-2J/(kBT)}

void monte_carlo_steps(int n)
{
   int i, k;
   for(k = 0; k < n; ++k)                 do the Wolff single-cluster
       {                                                         flip n times
     i = drand48() * (double) N;     pick a seed site at random
     flip(i, s[i]);                                flip a cluster starting from
                                                                  the seed
     }
  }

The next function is the core part which performs a Wolff single-cluster flip. This function is recursive. The array for spins s[ ], percolation probability p, and coordination number Z are passed globally.

The first argument i of the flip() function is the site to be flipped, the second argument s0 is the spin of the cluster before flipping. The function neighbor() is used again as in single-spin-flip Ising model or percolation program.

void flip(int i, int s0)
{
   int j, nn[Z];
   s[i] = - s0;                                     flip the spin immediately
   neighbor(i, nn);                             find nearest neighbor of i
   for(j = 0; j < Z; ++j){                      flip the neighbor if
   if(s0 == s[nn[j]] && drand48() < p)  spins are equal and
    flip(nn[j], s0);                               random number is smaller
                                                                      than p
  }

 

Metropolis algorithm and similar algorithms have wide applications in physics. Examples are quantum system and quantum field theory, fluids, macro-molecules and polymer systems. In this section, we'll describe briefly how to simulation fluid and lattice model of polymer systems.
 

Simple fluid
 

Let's imagine that we have a (two-dimensional) box of dimension L x L. Inside, we put certain number of particles. Depending on the density of the particles, we might get a fluid state or a solid state. We assume that the particles can be described as classical point particles. The particles interact with each other via a Lennard-Jones potential energy:

V(r) = 4e [ (sigma/r)12 - 2 (sigma/r)6 ]

where sigma ~ 3.5A (A = 10-10 m) and e/kb ~ 118K.
 

We are interested in the calculation of the thermodynamic quantities like internal energy, specific heat, equation of state (pressure as a function of temperature and volume), entropy, etc. A Monte Carlo sampling generates the particle configurations according to the Boltzmann distribution P ~ exp{-E/kB T}. This can be achieved by the Metropolis algorithm.

To perform a Monte Carlo move, one picks a particle at random and try to move it inside certain small region, say a square of side alpha with coordinates changing to:

x -> x + alpha (xi1 - 1/2)

y -> y + alpha (xi2 - 1/2)

where xi1 and xi2 are two independent uniformly distributed random numbers between 0 and 1.

We then calculate the change in energy of the system E'-E, due to this move.
 

• If E'-E < 0 , i.e., if the move would bring the system to a state of lower energy, we allow the move and put the particle in the new location.

• If E'-E > 0, we allow the move with probability 
  exp{ -(E'-E)/kBT}

i.e., we take a random number xi between 0 and 1, and if xi < exp{-(E'-E)/kBT}, we move the particle to the new position. Otherwise we return to its old position.

Note, however, whether we moved the particle or not, we always count it as one Monte Carlo step. Once we have an equilibrium distribution, it is possible to calculate the equation of state by virial theorem (relating pressure with certain average of force).

Simulation of polymer melt
 

Polymers are more complex than simple fluids. Again as before, we try to model the system as simple as possible. The lattice polymer models have been studied extensively in the literature.
 
 

To represent a polymer system, we take a square lattice (or cubic lattice in three dimensions) and occupy a set of N sites connected by bonds, forming a one-dimensional chain. The sites or bonds can not be occupied more than once. Each single occupied site is called a monomer, which is supposed to represent an atom or a group of atoms. Monomers of the same polymer or different polymers will have interactions. Usually, it is taken as nearest neighbor interactions.

There are a number of Monte Carlo simulation methods for polymers. Probably the easiest to implement is the reptation method, also called slithering snake algorithm. We try to move a polymer by erasing a monomer at one end of the chain and attaching it to the other end.

Again, we use the Metropolis method to decide whether we accept the move.
 

• If the total energy of the system decreases or stays the same, we allow the move.

• If the energy increases by an amount E'-E, we make the move only with a probability exp{-(E'-E)/kB T}.

Unsuccessful trial of move is also counted as a move (or Monte Carlo step). The polymer simulation algorithms can be used to study many properties of polymer systems.

Homework