PHYSICS 509, Assignment 4

Montecarlo Simulation of the 2-d Ising Model

Announcement: There will be no lecture on Tuesday April 11th Friday April 14th and Tuesday April 18th. Classes will resume on Friday April 21st

This problem set will explore the statistical mechanics of the two dimensional Ising model defined by:

H = −J ∑_{<i,j>} S_{i}S_{j}  (1)

The integers i and j range from 1 to N = L^{2} and describe a two dimensional lattice as described in class. We will assume periodic boundary condition whose meaning and implementation ws discussed in class. The sum in equation( ??) is over pairs of nearest neighbors on the square lattice. Each pair appears on the sum only once.

Equation ?? can be written more explicitly as

H = −J/2 ∑_{i} ∑_{δ=±x,±y} S_{i}S_{i+δ}  (2)

We will set J=1 and measure the temperature T in units of J. β = 1/T is the inverse temperature. In this units this model has (in the thermodynamical limit N → ∞ ) a phase transition when βJ = .5ln(1 + \sqrt{2}) ≈ .44069.

The goal is to evaluate the sum over spin configurations [S] as usual denotes a whole spin configuration.

< O >= \frac{∑_{S} O[S]exp−βH[S]}{∑_{S} exp−βH[S]}  (3)

using the Montecarlo simulation method.

An essential component of the simulation is the random number generator. It is instructive to write one, to see really how it works.

1) Program the simplest random number generator. An algorithm that generates a sequence of random numbers is:

real function randu(ibm) integer mult, modulo,ibm real rmodulo mult=1277 modulo= 2**17 rmodulo=modulo ibm=ibm*mult ibm=mod(ibm,modulo) randu=ibm/rmodulo end

This algorithm is infamous because it was widely used as a standard in older machines.

2) This random number generator can be used to generate a sequence of random numbers in (0,1) x_{1}, x_{2}, x_{3}, x_{4}, x_{5}, ...........

Test this random number generator by ploting pairs of consecutive random numbers, (x_{1}, x_{2}), (x_{3}, x_{4}), ...... Plot something of the order of 2000-3000 pairs, can you discern the correlations?

3) Now repeat step 2 with ‘ran3’, the random number generator provided by the numerical recipes, and compare it with the results obtained in 2.

Now we can begin to explore the Ising model proper.

4) Write a Montecarlo simulation program for the 2d ising model on an L x L lattice using the metropolis algorithm and periodic boundary conditions.

The various steps were discussed in class. It is customary to set up a look-up table of exp(βh), with the finite and small number of values that the internal field h can take, to avoid computing the exponential function each time.

5) Run your program first on a system of very modest size L=4 (16 sites). Measure the magnetization and the energy as a function of temperature for β in the interval range from .1 to 1.0. As discuss in class to estimate the magnetization measure 1/N| ∑_{i} s_{i}|.

6) for the size L=4 one can evaluate ?? exactly using exact enumeration using the Gray code (discussed in class and provided with the numerical recipes), instead of the Montecarlo procedure. Pick a single temperature, since the only objective of this step is to make sure that your code really works.

Explore how the statistical error is reduced as you increse the length of the Montecarlo run, by repeating item 5 with different lenths of your montecarlo run.(Start with short runs) Watch the results of the simulations approach the exact enumeration result when the runs get longer.

7) Now that you have a rough idea of what the system does as a function of temperature and you are sure that your code really works, its time to estimate the errors in the measurements performed in 5. Estimate the errors in your measurements and the correlation time. As discussed in class this can be accomplish by studying how the variance depends on the size of the ensemble of measurements used to calculate the variance.

8) One gain insights into the statistical mechanics of this problem by obtaining a histogram of the variableM = ∑_{i} s_{i} The range of this variable is [−L^{2}, L^{2}] and it can take only even values. Compute the histogram of M at β = .3 and β = .7, β = .55, and β = .1. Give a physical interpretation of the results.

9)Phase transitions occur only in the thermodynamic limit, so it is important to understand the dependence of the results on the size L. Repeat step 5, with a size L=10 (100 sites)and compare the results with those obtained in step 5.

DIRECTIONS FOR FURTHER STUDY. THIS PART IS OPTIONAL. If you went as far as doing part 8) you may be curious what would the

probability distribution of the M = ∑_{i} s_{i} looks like if you are near the critical point. (That is when the correlation length is larger than the sample size). That cannot be studied with the simple single flip algorithm, because of the dramatic increase of the autocorrelation time. However the CLUSTER algorithm that we described in the lecture is fairly easy to program (see Werner Krauth’s lectures, I put a link to it in my Web Page). and would allow you to see what the histogram of the magnetizations looks like near the transition.