Compiled functions

spinNRGnopbc=Compile[{{ng,_Integer,2}},Module[{sz},sz=Length[ng];-Sum[ng[[i,j]]*ng[[i+1,j]]+ng[[j,i]]*ng[[j,i+1]],{i,1,sz-1},{j,1,sz}]],{{sz,_Integer}}];
spinNRGpbc=Compile[{{ng,_Integer,2}},Module[{sz},sz=Length[ng];-Sum[ng[[i,j]]*ng[[i+1,j]]+ng[[j,i]]*ng[[j,i+1]],{i,1,sz-1},{j,1,sz}]
-Sum[ng[[1,i]]*ng[[sz,i]]+ng[[i,1]]*ng[[i,sz]],{i,1,sz}]
]
,{{sz,_Integer}}];

prepSTEPnopbc=Compile[{{ng,_Real,2}},Module[{sz,ix,iy,tn},sz=Length[ng];
ix=Floor[RandomReal[]*Length[ng]]+1;iy=Floor[RandomReal[]*Length[ng]]+1;
tn=If[ix>1,ng[[ix-1,iy]],0]+If[iy>1,ng[[ix,iy-1]],0]+If[ix<sz,ng[[ix+1,iy]],0]+If[iy<sz,ng[[ix,iy+1]],0];
{ix,iy,tn*ng[[ix,iy]]}
]
,{{sz,_Integer},{ix,_Integer},{iy,_Integer},{tn,_Integer}},CompilationTarget->"C"
];
multiSTEPnopbc=Compile[{{ng,_Integer,2},{temp,_Real,0},{b,_Real,0},{nsteps,_Integer,0}},
Module[{sz,ix,iy,tn,lg,dbeta},lg=ng;sz=Length[lg];
dbeta=2/temp;
Do[
ix=RandomInteger[{1,sz}];iy=RandomInteger[{1,sz}];
tn=If[ix>1,lg[[ix-1,iy]],0]+If[iy>1,lg[[ix,iy-1]],0]+If[ix<sz,lg[[ix+1,iy]],0]+If[iy<sz,lg[[ix,iy+1]],0];
If[RandomReal[{0,1}]<Exp[-dbeta*(tn+b)*lg[[ix,iy]]],lg[[ix,iy]]*=-1]
,{k,1,nsteps}];
lg
]
,{{ix,_Integer},{iy,_Integer},{tn,_Integer}},CompilationTarget->"C"];
prepSTEPpbc=Compile[{{ng,_Real,2}},Module[{sz,ix,iy,tn},sz=Length[ng];
ix=Floor[RandomReal[]*Length[ng]]+1;iy=Floor[RandomReal[]*Length[ng]]+1;
tn=If[ix>1,ng[[ix-1,iy]],ng[[sz,iy]]]+If[iy>1,ng[[ix,iy-1]],ng[[ix,sz]]]+If[ix<sz,ng[[ix+1,iy]],ng[[1,iy]]]+If[iy<sz,ng[[ix,iy+1]],ng[[ix,1]]];
{ix,iy,tn*ng[[ix,iy]]}
]
,{{sz,_Integer},{ix,_Integer},{iy,_Integer},{tn,_Integer}},CompilationTarget->"C"
];
multiSTEPpbc=Compile[{{ng,_Integer,2},{temp,_Real,0},{b,_Real,0},{nsteps,_Integer,0}},
Module[{sz,ix,iy,tn,lg,dbeta},lg=ng;sz=Length[lg];
dbeta=2/temp;
Do[
ix=RandomInteger[{1,sz}];iy=RandomInteger[{1,sz}];
tn=If[ix>1,lg[[ix-1,iy]],lg[[sz,iy]]]+If[iy>1,lg[[ix,iy-1]],lg[[ix,sz]]]+If[ix<sz,lg[[ix+1,iy]],lg[[1,iy]]]+If[iy<sz,lg[[ix,iy+1]],lg[[ix,1]]];
If[RandomReal[{0,1}]<Exp[-dbeta*(tn+b)*lg[[ix,iy]]],lg[[ix,iy]]*=-1]
,{k,1,nsteps}];
lg
]
,{{ix,_Integer},{iy,_Integer},{tn,_Integer}},CompilationTarget->"C"];

Mathematica wrappers

calcEnergy[st_,b_:0]:=(If[st["periodic"],(st["energy"]=(spinNRGpbc[st["grid"]]*1.0-b*Total[Flatten[st["grid"]]])/st["size"]^2),(st["energy"]=(spinNRGnopbc[st["grid"]]-b*Total[Flatten[st["grid"]]])*1.0/st["size"]^2)];
st["magnetization"]=Total[Flatten[st["grid"]]]*1.0/st["size"]^2;);
SetAttributes[calcEnergy,HoldFirst]

Module[{lst},
initState[size_,periodic_:True,mode_:"random"]:=(lst=<|"size"-> size,"grid"->
Switch[mode,"random",
(2*RandomInteger[{0,1},{size,size}]-1),
"up",ConstantArray[1,{size,size}],
"down",ConstantArray[-1,{size,size}],
"disk",Table[If[(i/size-0.5)^2+(j/size-0.5)^2<0.125,-1,1],{i,1,size},{j,1,size}],
"stripe",Table[If[(i/size-0.5)^2<0.125,-1,1],{i,1,size},{j,1,size}]
]
,"periodic"->periodic,"energy"->None,"magnetization"->None|>; calcEnergy[lst,0];lst)]

Module[{pstep},
mcStep[st_,temp_,b_,nsteps_]:=(If[st["periodic"],
st["grid"]=multiSTEPpbc[st["grid"],temp,b,nsteps],
st["grid"]=multiSTEPnopbc[st["grid"],temp,b,nsteps]];
calcEnergy[st,b];
);
SetAttributes[mcStep,HoldFirst]
];

state=initState[200,True,"random"];
sdata={{0,state["energy"],state["magnetization"]}}

Dynamic

helper (drawing) functions

Module[{revg},
isingDraw[st_,{hi_,hj_},b_:0]:=Graphics[
{EdgeForm[Black],
Table[{If[st["grid"][[i,j]]>0,RGBColor[1,0.5,0.5],RGBColor[0.5,0.5,1]],Rectangle[{i-1,j-1},{i,j}]},{i,1,st["size"]},{j,1,st["size"]}],
If[hi hj >0,{Opacity[0.5],Red,EdgeForm[Directive[Red,Thick]],Rectangle[{hi-1,hj-1},{hi,hj}],Opacity[1]},{}],
Arrowheads[0.05],Black,Table[{If[st["grid"][[i,j]]>0,Arrow[{{i-0.5,j-0.75},{i-0.5,j-0.25}}],Arrow[{{i-0.5,j-0.25},{i-0.5,j-0.75}}]]},{i,1,st["size"]},{j,1,st["size"]}],
Text[Style[ToString["E/site : "<>ToString[SetPrecision[st["energy"],4]]],Bold],Offset[{5,-5},{0,0}],{-1,1}],
Text[Style[ToString["μ/site : "<>ToString[SetPrecision[st["magnetization"],4]]],Bold],Offset[{5,-20},{0,0}],{-1,1}],
If[hi hj >0,revg=st;revg[[Key["grid"]]][[hi,hj]]*=-1;calcEnergy[revg,b];
Text[Style[ToString["ΔE : "<>ToString[Round[st["size"]^2(st["energy"]-revg["energy"])] ]],Bold,Large],Offset[{-5,-5},{st["size"],0}],{1,1}],{}
]
}]
]

Stepper

{i,j}={0,0}; de=0;
nsz=8;over={};sstate=initState[nsz];
stepper=Manipulate[Dynamic[Show[isingDraw[sstate,{i,j}],Graphics[over]]],
{{speri,True,"Periodic boundaries"},{True,False}},
{{sstyle,"random","Initialization style"},{"random","up","down","disk"}},
Button["reset",
sstate=initState[nsz,speri,sstyle],ImageSize->80],
{{temp,0.5,"Temperature"},0.001,5,Appearance->"Labeled"},
Button["step",
{i,j,de}=If[sstate["periodic"],prepSTEPpbc,prepSTEPnopbc][sstate["grid"]];
sstate[[Key["grid"] ]][[i,j]]*=-1;
calcEnergy[sstate];
Pause[1.5];
If[RandomReal[]<Exp[-de*2/temp],
over={Text[Style["ACCEPTED",Bold,Larger,Red],{nsz/2,nsz/2},{0,0},{1,1},Background->White]};
Pause[1.5];
,
over={Text[Style["REJECTED",Bold,Larger,Red],{nsz/2,nsz/2},{0,0},{1,1},Background->White]};
Pause[1.5];
sstate[[Key["grid"] ]][[i,j]]*=-1;calcEnergy[sstate];
];
{i,j}={0,0};over={};Pause[0.5];
,ImageSize->80,Method->"Queued"],
SaveDefinitions->True
]

CDFDeploy[NotebookDirectory[]<>"lab3_stepper.cdf",stepper]

Runner

(state=initState[100];
sdata={{0,state["energy"],state["magnetization"],state["energy"],state["magnetization"],state["energy"]^2,state["magnetization"]^2}};it=0;ii=1;
 buff=1000; fgo=False; fplot=True; gbuff={};iibuff=ii;)
runner=Manipulate[Dynamic[
If[fplot || sdyna,fplot=False;
iibuff=ii;
gbuff=Graphics[{Inset[
mpl=MatrixPlot[ state["grid"],FrameTicks->None,ImageSize->Full],Scaled[{0.25,0.5}],Center,Scaled[{0.5,1}]],
Inset[sp1=Show[ListPlot[{sdata[[All,{1,2}]],sdata[[All,{1,4}]]},Frame->True,Joined->True,Axes->False,PlotStyle->{Red,Directive[Red,Dashed]},FrameLabel->{"million steps","energy/site"},ImageSize->Full,AspectRatio->0.5],
Graphics[{
Text[Style["E="<>ToString[SetPrecision[state["energy"],4]],Bold],Scaled[{0.95,0.95}],{1,1},Background->White ],
Text[Style["<E>="<>ToString[SetPrecision[sdata[[-1,4]],4]],Bold],Scaled[{0.95,0.75}],{1,1},Background->White ],
Text[Style["(δ^2)E="<>ToString[SetPrecision[sdata[[-1,6]]-sdata[[-1,4]]^2,4]],Bold],Scaled[{0.95,0.55}],{1,1},Background->White ]
}],ImageSize->Full
],Scaled[{0.75,0.75}],Center,Scaled[{0.5,0.5}]],
Inset[sp2=Show[
ListPlot[{sdata[[All,{1,3}]],sdata[[All,{1,5}]]},Frame->True,Joined->True,Axes->False,PlotStyle->{Blue,Directive[Blue,Dashed]},
FrameLabel->{"million steps","magnetization/site"},ImageSize->Full,AspectRatio->0.5],
Graphics[{
Text[Style["M="<>ToString[SetPrecision[state["magnetization"],4]],Bold],Scaled[{0.95,0.95}],{1,1},Background->White],
Text[Style["<M>="<>ToString[SetPrecision[sdata[[-1,5]],4]],Bold],Scaled[{0.95,0.75}],{1,1},Background->White ],
Text[Style["(δ^2)M="<>ToString[SetPrecision[sdata[[-1,7]]-sdata[[-1,5]]^2,4]],Bold],Scaled[{0.95,0.55}],{1,1},Background->White ]
}],ImageSize->Full
],Scaled[{0.75,0.25}],Center,Scaled[{0.5,0.5}]]
},AspectRatio->0.5,ImageSize->Large]
,If[ii==iibuff,gbuff,Show[Graphics[{Inset[gbuff,Center,Center,Scaled[1]],Opacity[0.3],Gray,Rectangle[Scaled[{0,0}],Scaled[{1,1}]],
Red,Opacity[1],Text[Style["Computing.... "<>ToString[ii],Bold,Large],Scaled[{0.5,0.5}]]
}],PlotRange->All,AspectRatio->0.5,ImageSize->Large]
]
]
,
UpdateInterval->Infinity,TrackedSymbols:> {state, fplot, sdyna,ii}
],
{{ns,100,"Size"},4,500,5,Appearance->"Labeled"},
{{speri,True,"Periodic boundaries"},{True,False}},
{{sstyle,"random","Initialization style"},{"random","up","down","disk", "stripe"}},
{{sbuff,1000,"Buffer size"},100,10000,100,Appearance->"Labeled"},
Button["reset",
state=initState[ns,speri,sstyle];
sdata={{0,state["energy"],state["magnetization"],state["energy"],state["magnetization"],state["energy"]^2,state["magnetization"]^2}};buff=sbuff;it=0;ii=1;,ImageSize->80],
{{temp,0.5,"Temperature"},0.001,5,Appearance->"Labeled"},
{{field,0.0,"Magnetic Field"},-2,2,Appearance->"Labeled"},
{{stride,100000,"Number of MC steps"},{1,1000,100000,1000000,10000000}},
{{sdyna,True,"Dynamic update"},{True,False}},
Row[{Button["run",fgo=True;While[fgo,
mcStep[state,temp,field,stride];
If[ii<buff,
sdata=Append[sdata,{it,state["energy"],state["magnetization"],(sdata[[-1,4]]*ii+state["energy"])/(ii+1),(sdata[[-1,5]]*ii+state["magnetization"])/(ii+1),
(sdata[[-1,6]]*ii+state["energy"]^2)/(ii+1),(sdata[[-1,7]]*ii+state["magnetization"]^2)/(ii+1)
}],
sdata=Drop[Append[sdata,{it,state["energy"],state["magnetization"],
(sdata[[-1,4]]*buff+state["energy"]-sdata[[1,2]])/buff,(sdata[[-1,5]]*buff+state["magnetization"]-sdata[[1,3]])/buff,(sdata[[-1,6]]*buff+state["energy"]^2-sdata[[1,2]]^2)/buff,(sdata[[-1,7]]*buff+state["magnetization"]^2-sdata[[1,3]]^2)/buff
}],1]
]; (*switches to moving averages...*)
ii++;it+=stride*1.*^-6],ImageSize->80,Method->"Queued"],"  ",
Button["stop",fplot=True;fgo=False,ImageSize->80],"  ",
Button["reset avg",
sdata={{0,state["energy"],state["magnetization"],state["energy"],state["magnetization"],state["energy"]^2,state["magnetization"]^2}};buff=sbuff;it=0;ii=1;,ImageSize->80]
}],
Row[{
Button["save data",
Export[NotebookDirectory[]<>"/ising_image.png",mpl,"PNG"];
Export[NotebookDirectory[]<>"/ising_data.csv",Prepend[sdata[[;;,{1,2,3}]],{"step/1e6","energy/site","magnetization/site"}],"CSV"];
,ImageSize->80],"  ",
Button["print avg.",
MessageDialog["<E> = "<>ToString[CForm[sdata[[-1,4]]]] <>
"\nvar(E) = "<>ToString[CForm[sdata[[-1,6]]-sdata[[-1,4]]^2]] <>
"\n<M> = "<>ToString[CForm[sdata[[-1,5]]]] <>
"\nvar(M) = "<>ToString[CForm[sdata[[-1,7]]-sdata[[-1,5]]^2]] ,Selectable->True,Deployed->False
]
,ImageSize->80]}],
SaveDefinitions->True
]

CDFDeploy[NotebookDirectory[]<>"lab3_ising.cdf",runner]

mp=MatrixPlot[state["grid"],FrameTicks->None,ImageSize->Large]

Export["/tmp/ising.png",mp,"PNG"]

Graphics[{Inset[mpl,Scaled[{0.25,0.5}],Center,Scaled[{0.5,1}]],
Inset[sp1,Scaled[{0.75,0.75}],Center,Scaled[{0.5,0.5}]],
Inset[sp2,Scaled[{0.75,0.25}],Center,Scaled[{0.5,0.5}]]
 },AspectRatio->0.5]

GraphicsRow[{mpl,GraphicsColumn[{sp1,sp2},Frame->All,Spacings->0,ImagePadding->0,AspectRatio->Full,Method->{"ShrinkWrap"->True}]},Frame->All,Spacings->0,PlotRangePadding->0,ImagePadding->0,Method->{"ShrinkWrap"->True}]

Analyzer

RGReduceAVG[grid_]:=Table[
Sign[(
grid[[Mod[(i+j),Length[grid]]+1,Mod[(i-j),Length[grid]]+1 ]]+
grid[[Mod[(i+j)+1,Length[grid]]+1,Mod[(i-j),Length[grid]]+1 ]]+
grid[[Mod[(i+j)-1,Length[grid]]+1,Mod[(i-j),Length[grid]]+1 ]]+
grid[[Mod[(i+j),Length[grid]]+1,Mod[(i-j)+1,Length[grid]]+1 ]]+
grid[[Mod[(i+j),Length[grid]]+1,Mod[(i-j)-1,Length[grid]]+1 ]]
)]
,{i,0,Length[grid]-1},{j,0,Length[grid]-1}]
RGReduce[grid_]:=Table[
grid[[Mod[(i+j),Length[grid]]+1,Mod[(i-j),Length[grid]]+1 ]]
,{i,0,Length[grid]-1},{j,0,Length[grid]-1}]

compiled version

CRGReduce=Compile[{{grid,_Real,2}},
Table[grid[[Mod[(i+j),Length[grid]]+1,Mod[(i-j),Length[grid]]+1 ]]
,{i,0,Length[grid]-1},{j,0,Length[grid]-1}]]

RGPlot[grid_,acf_,log_:False]:=GraphicsRow[{
MatrixPlot[grid,FrameTicks->None,ImageSize->Medium],
actime=Max[Total[acf]*2+1,16];
Show[
(If[log,ListLogPlot,ListPlot])[acf,DataRange->{0,Length[acf]-1},Frame->True,Joined->True,Axes->False,PlotStyle->Blue,PlotRange->{{0,actime},All},
FrameLabel->{"Δi","<(Subscript[s, i])(Subscript[s, i + Δi])>"},ImageSize->Medium,AspectRatio->0.5]
]}]

Module[{fl},
doACF[l_,mean_:Automatic,var_:Automatic]:=(
fl=Fourier[l];
fl=fl*Conjugate[fl];
fl=Re[InverseFourier[fl]]/Sqrt[Length[fl]];
fl=(fl-If[mean===Automatic,Mean[l]^2,mean^2])/If[var===Automatic,Variance[l],var]
)]

(state=initState[256];fnrg=0;
facf=tacf=ConstantArray[0,Floor[Length[state["grid"]]/2]];ii=1; iibuff=1;iacf=0;
iacf=0;rngpanel=Show[GraphicsRow[{Graphics[{},ImageSize->Large],Graphics[{},ImageSize->Large]}]];
;
 fgo=False;fan=False;fplot=True; )
rngmanip=Manipulate[Dynamic[
If[fplot || sdyna,fplot=False;iibuff = ii;
gbuff=Column[
Join[
{RGPlot[state["grid"],tacf/Max[1,iacf] ,logscale ]},
{rngpanel}
]
],
If[ii==iibuff,gbuff,Show[Graphics[{Inset[gbuff,Center,Center,Scaled[1]],Opacity[0.3],Gray,Rectangle[Scaled[{0,0}],Scaled[{1,1}]],
Red,Opacity[1],Text[Style["Computing.... "<>ToString[ii],Bold,Large],Scaled[{0.5,0.5}]]
}]
]]
],
UpdateInterval->Infinity,TrackedSymbols:> {state, fplot, sdyna,logscale,ii}
]
,
{{ns,256,"Size"},Table[2^i,{i,4,9}],ControlType->SetterBar},
{{speri,True,"Periodic boundaries"},{True,False}},
{{sstyle,"random","Initialization style"},{"random","up","down","disk"}},
Button["reset",
state=initState[ns,speri,sstyle];
facf=tacf=ConstantArray[0,Floor[ns/2]];iacf=0;ii=1;fplot=True;
 fgo=False;iacf=0;rngpanel=Show[GraphicsRow[{Graphics[{},ImageSize->Large],Graphics[{},ImageSize->Large]}]];;
,ImageSize->80],
{{temp,2,"Temperature"},0.001,5,Appearance->"Labeled"},
{{field,0.0,"Magnetic Field"},-2,2,Appearance->"Labeled"},
{{stride,1000000,"Number of MC steps"},{1,100,10000,100000,1000000,1*^7}},
{{logscale,False,"Log y scale"},{True,False},Appearance->"Labeled"},
{{sdyna,True,"Dynamic update"},{True,False}},
Button["run",fgo=True;While[fgo,ii+=1; mcStep[state,temp,field,stride];]
,ImageSize->80,Method->"Queued"],
Button["stop",fgo=False;fplot=True,ImageSize->80],
Button["analyze",iacf++;
mg=Mean[Flatten[state["grid"]]];vg=Variance[Flatten[state["grid"]]];
facf=Take[(Mean[doACF[#,mg,vg]&/@state["grid"]]+Mean[doACF[#,mg,vg]&/@Transpose[state["grid"]]])/2,Floor[ns/2]];
tacf+=facf;fplot=True;
,ImageSize->80],
Button["clear",facf=tacf=ConstantArray[0,Floor[ns/2]];iacf=0;fplot=True;,ImageSize->80],
Delimiter,
{{nrg,0,"RNG Levels steps"},Table[i,{i,0,Log[2,ns]-2}],ControlType->SetterBar},
Panel[Column[{
Button["analyze frame",fgo=False;fan=True;fplot=True;
rngpanel=Column[{Row[
rp=state["grid"];rpacf=ConstantArray[0,nrg];
Table[
rp=CRGReduce[rp];
mm=Mean[Flatten[rp]];vv=Variance[Flatten[rp]];
facf=Take[(Mean[doACF[#,mm,vv]&/@rp]+Mean[doACF[#,mm,vv]&/@Transpose[rp]])/2,Floor[Length[rp]/2]];
rpacf[[i]]=facf;
Show[MatrixPlot[rp,ImageSize->8*72/nrg,FrameTicks->None],Graphics[{Text[Style[ToString[i],Bold,Large,RGBColor[i/nrg,0,1-i/nrg]],Offset[{5,5},{1,1}],{-1,-1},Background->White]}]]
,{i,1,nrg}]
],
cl=Floor[Total[Abs[tacf/Max[1,iacf]]]*2+5];
If[logscale,ListLogPlot,ListPlot][Evaluate[Join[{tacf/Max[1,iacf]},Table[Join[ri,ConstantArray[0,Length[tacf]-Length[ri]]],{ri,rpacf}]],DataRange->{0,Length[tacf]-1}]
,Frame->True,Joined->True,Axes->False,PlotStyle->Join[{Black},Table[RGBColor[i/nrg,0,1-i/nrg],{i,1,nrg}]],PlotRange->{{0,cl},All},
FrameLabel->{"Δi","<(Subscript[s, i])(Subscript[s, i + Δi])>"},ImageSize->Medium,AspectRatio->0.5]
}]
,ImageSize->120,Method->"Queued"],
Button["Save data",Export[NotebookDirectory[]<>"rng_data.csv",Transpose[rpacf]]]
}],Style["RNG ANALYSIS",Bold]],
SaveDefinitions->True
]

CDFDeploy[NotebookDirectory[]<>"lab3_rng.cdf",rngmanip]