Supplemental Notebook for 2D Ising Model

Below is an implementation of the heat-bath algorithm discussed in class.  The code is commented so that it should be easy to follow and understand.  The code takes as inputs the size of the lattice Nspins, the temperature T, the Ising strength J, and the number of iterations Nev to update the system.  The function heatbath outputs the grid of up and down spins that can be plotted.  As the code runs, it displays the current spin system, with up spins blue and down spins yellow.  Periodic boundary conditions are imposed in this code.

heatbath[Nspins_,T_,J_,Nev_]:=Module[{spinstab,i,j,x,y,leftn,rightn,upn,downn,mval,rand,pup,plottab},

(*Initialize the lattice of spins to be all up (+1)*)
spinstab=Table[1,{i,1,Nspins},{j,1,Nspins}];
plottab=ArrayPlot[spinstab,ColorFunction->(If[#==1,Blue,Yellow]&)];

(*Run over all iterations*)
Monitor[
For[i=1,i≤Nev,i++,

(*Find a random lattice site*)
x=RandomInteger[{1,Nspins}];
y=RandomInteger[{1,Nspins}];

(*Determine the neighbors.  Need to consider if the chosen spin is at the boundary.  If so, we impose periodic boundary conditions.*)
leftn=If[x≠1,x-1,Nspins];
rightn=If[x≠Nspins,x+1,1];
upn=If[y≠1,y-1,Nspins];
downn=If[y≠Nspins,y+1,1];

(*Calculate the neighbors' spin sum*)
mval=spinstab[[leftn,y]]+spinstab[[rightn,y]]+spinstab[[x,upn]]+spinstab[[x,downn]];

(*Calculate the probability to have up (+1) spin*)
pup=((E^(J mval/T))/((E^(J mval/T))+(E^(-J mval/T))))//N;
rand=RandomReal[];
If[rand<pup,spinstab[[x,y]]=1,spinstab[[x,y]]=-1];

If[Mod[i,100000]==0,plottab=ArrayPlot[spinstab,ColorFunction->(If[#==1,Blue,Yellow]&),PlotLabel->i]];
],plottab];
Print[plottab];
Return[spinstab];
];

Here are some plots of the 2D Ising model on a grid of 200x200 with J=1, at various temperatures.

spins15=heatbath[200,1.5,1,10000000];

spins24=heatbath[200,2.4,1,10000000];

spins32=heatbath[200,3.2,1,10000000];

Below is a function that calculates the spin correlation for the 2D Ising model on a square lattice.  It takes arguments of the collection of spins to analyze (tab), the dimension of the spin lattice (Nspins) and the maximum distance to consider (lmax).

correlation[tab_,Nspins_,lmax_]:=Module[{i,x,y,d,rightn,leftn,upn,downn,spincorr,spinprod},
spincorr=Table[{i,0.0},{i,1,lmax}];
For[x=1,x≤Nspins,x++,
For[y=1,y≤Nspins,y++,
For[d=1,d≤Min[lmax,Nspins],d++,

leftn=If[x>d,x-d,Nspins+x-d];
rightn=If[x≤Nspins-d,x+d,x+d-Nspins];
upn=If[y>d,y-d,Nspins+y-d];
downn=If[y≤Nspins-d,y+d,y+d-Nspins];

spinprod=((tab[[x,y]])/4)*(tab[[leftn,y]]+tab[[rightn,y]]+tab[[x,upn]]+tab[[x,downn]]);

spincorr[[d,2]]+=(spinprod/(Nspins*Nspins));
];
];
];
Return[spincorr];
];

With this function, we can then calculate the spin correlations of the systems that we simulated above.  I’ve also plotted the spin resulting spin correlations as a function of the separation distance l.

corrtab32=correlation[spins32,200,20];corrtab24=correlation[spins24,200,20];corrtab15=correlation[spins15,200,20];

ListLinePlot[{corrtab15,corrtab24,corrtab32},PlotRange->All,Axes->False,Frame->True,PlotLabel->"Spin Correlation",PlotLabels->Placed[{"T = 1.5","T = 2.4","T = 3.2"},{Above}],FrameLabel->{"l","C(β,J,l)"}]

We can fit each of these expressions to the expected form of the spin correlation as an exponential decay and estimate the correlation length.

FindFit[corrtab15,a+b (E^(-x/ξ)),{a,b,ξ},x]
FindFit[corrtab24,a+b (E^(-x/ξ)),{a,b,ξ},x]
FindFit[corrtab32,a+b (E^(-x/ξ)),{a,b,ξ},x]

The correlation length is largest for T = 2.4, which is the closest to the critical temperature!