Monte Carlo Statistics for a 2D Ising Model with Periodic Boundary Conditions

(* Set the Size of the Lattice *)
rows=16; (* number of rows and columns - best set to be an even number *)

(* Generate "Interface" Initial Configuration *)
maxj=rows/2; (* the last row to be filled with +1's *)
nplus=rows*maxj; (* running count of +1's *)
A[0,0]=0;A[0,rows+1]=0;A[rows+1,0]=0;A[rows+1,rows+1]=0; (* the four corners *)
Do[A[i,j]=+1,{i,1,rows},{j,1,maxj}];
Do[A[i,j]=-1,{i,1,rows},{j,maxj+1,rows}];

(* Generate "Checkerboard" Initial Configuration *)
nplus=rows^2/2; (* ok if rows is an even number *)
A[1,1]=+1;
Do[A[i+1,1]=-A[i,1],{i,1,rows}];
Do[A[i,j+1]=-A[i,j],{i,1,rows},{j,1,rows}];

(* Impose Periodic Boundary Conditions, Calculate Averages *)
Do[{A[i,0]=A[i,rows],A[i,rows+1]=A[i,1]},{i,1,rows}];
Do[{A[0,j]=A[rows,j],A[rows+1,j]=A[1,j]},{j,1,rows}];
means[0]=0.; Do[means[0]=means[0]+A[i,j],{i,1,rows},{j,1,rows}]; means[0]=means[0]/rows^2;
corr[0]=0.; Do[corr[0]=corr[0]+A[i,j]*(A[i,j+1]+A[i+1,j]),{i,1,rows},{j,1,rows}];
corr[0]=corr[0]/(2*rows^2)-means[0]^2;

(* Generate and Displey Graphics from a Configuration *)
ipl=0; size=0.1;
Do[If[A[i,j]==1,
	{ipl=ipl+1,draw[ipl]=Graphics[Rectangle[{size*j,size*i},{size*(j+1),size*(i+1)}]]}],
{i,0,rows+1},{j,0,rows+1}];
Show[Array[draw,ipl],Axes->True,AspectRatio->1,PlotRange->{{0,size*(rows+2)},{0,size*(rows+2)}}]

(* Evolve the System and Calculate Average Spin *)
Tinv=0.60; ncount=10000;
(Do[{m=RandomInteger[{1,rows}],n=RandomInteger[{1,rows}], (* pick element to flip *)
S=-A[n,m],B=Tinv S (A[n-1,m]+A[n+1,m]+A[n,m-1]+A[n,m+1]) 2, (* B is Delta E *)
If[RandomReal[]<(E^B),{A[n,m]=S,nplus=nplus+S,If[n==1,A[rows+1,m]=S],
If[n==rows,A[0,m]=S],If[m==1,A[n,rows+1]=S],If[m==rows,A[n,0]=S]}], (* decide *)
means[ict]=0,
Do[means[ict]=means[ict]+A[i,j],{i,1,rows},{j,1,rows}],(* updated mean spin *)
means[ict]=((means[ict])/(rows^2)),Null}, {ict,1,ncount}];)

(* Display Evolution of Average Spin *)
ListPlot[Array[means,ncount],PlotRange->{-1,1}]

(* Calculate and Display Correlations of Spins One Site Apart *)
corr[ict]=0.;
	Do[corr[ict]=corr[ict]+A[i,j]*(A[i,j+1]+A[i+1,j]),{i,1,rows},{j,1,rows}];
	corr[ict]=corr[ict]/(2*rows^2)-means[ict]^2;

ListPlot[Array[corr,ncount],PlotRange->{-1,1}]