Ising Model

parameters

width of the square lattice

size=10

temperature

Temp=2.5

plot spin configuration as black and white squares

reperesent the spin configuration as a Table that can be plotted using ArrayPlot,
with black squares for spin +1 and white squares for spin -1

spinconfig:=Table[(spin[i,j]+1)/2,{i,1,size},{j,1,size}]

energy difference from flipping the spin at site i, j

(periodic boundary conditions at the edges)

deltaU[i_,j_]:=Module[{up,down,left,right},
If[i==1,up=spin[size,j],up=spin[i-1,j]];
If[i==size,down=spin[1,j],down=spin[i+1,j]];
If [j==1,left=spin[i,size],left=spin[i,j-1]];
If[j==size,right=spin[i,1],right=spin[i,j+1]];
2*spin[i,j]*(up+down+right+left)]

update the spin at site i, j 
            using the Metropolis algorithm

calculate the energy difference Ediff  from flipping the spin
if flipping the spin lowers the energy, flip it
if flipping the spin raises the energy, flip it with probability given by the Boltzmann factor

metropolis[i_,j_]:=Module[{},
Ediff=deltaU[i,j];
If[Ediff<0,
    spin[i,j]=-spin[i,j],
    If[RandomReal[]<Exp[-Ediff/Temp],
        spin[i,j]=-spin[i,j],
        spin[i,j]=spin[i,j]]]];

initialize the spins randomly to +1 or -1

use a random number generator to initialize each spin to +1 or -1 
and plot the spin configuration

Do[spin[i,j]=Sign[2*RandomReal[]-1],{i,1,size},{j,1,size}];
Print[Table[spin[i,j],{i,1,size},{j,1,size}]];
ArrayPlot[spinconfig]

sweep randomly through the lattice, updating spins

the number of iterations in a sweep is 100 times the number of sites

iterations = 100*size^2

for each iteration in the sweep,
choose a site i, j  at random and flip the spin according to the Metropoois algorithm

sweep:=Do[
i=IntegerPart[RandomReal[]*size+1];j=IntegerPart[RandomReal[]*size+1];
metropolis[i,j],
{n,1,iterations}]

after sweeping through the lattice, 
plot the spin configuration

sweep
ArrayPlot[spinconfig]

sweep
ArrayPlot[spinconfig]

plot the spin configuration history for an entire sweep

record the spin configuration at each step of the next sweep

sweephistory=Table[(metropolis[IntegerPart[RandomReal[]*size+1],IntegerPart[RandomReal[]*size+1]]; spinconfig),{r,1,iterations}];

turn on Dynamic Updating Enabled under the Evaluation menu in order to Animate

Animate[ArrayPlot[sweephistory[[n]]],{n,1,iterations,1},AnimationRepetitions->1,AnimationRunning->False]