Sample code 8.1. A numerical code to solve the implicit function 

for <s>in 1D Ising Model.

mu=1;(* in unit of a Bohr magneton *)
z=4; (*# of nn in 2-D *)
J=+1;B=0;k=1;T=0.1;
P = mu*B/(k*T);Q=J*z/(k*T);
f[s_]:=s-Tanh[P+Q*s];
fprime[s_]:=1-Q*Sech[P+ Q*s]^2;
tol=1.0*10^(-4);
Deltaroot=1.0*10^6;
s=10.0;(*initial guess of root*)
While[Deltaroot>tol,
sold=s;
s=s-f[s]/fprime[s];
Deltaroot=Abs[s-sold];
Print[Deltaroot,"  ",s];
];(*end while*)
root=s;
Print["the root is ",root," Deltaroot=",Deltaroot];
Print["f[sroot]=",f[root]];

Sample code 8.2. Now the code can be modified to find <s> as  a fucntion of kT.

mu=1;(* in unit of a Bohr magneton *)
z=4; (*# of nn in 1-D *)
J=+1;B=0.0;k=1;(*T=30;*)
P := mu*B/(k*T);Q:=J*z/(k*T);
f[s_]:=s-Tanh[P+Q*s];
fprime[s_]:=1-Q*Sech[P+ Q*s]^2;
tol=1.0*10^(-4);

NT=500;
Tlast=10; Tini=0; DeltaT=N[(Tlast-Tini)/NT];
epsilon=10^(-6);
tTroot={};
For[iT=0,iT≤NT,iT++,
T=iT*DeltaT+epsilon;

s=10.0;(*initial guess of root*)
Deltaroot=1.0*10^6;

While[Deltaroot>tol,
sold=s;
s=s-(f[s]/fprime[s]);
Deltaroot=Abs[s-sold];
(*Print[Deltaroot,"  ",s];*)
];(*end while*)
root=s;
tTroot=Append[tTroot,{k*T,root}];
(*Print["T, root :",T, " ",root];*)
];(*end for iT*)
ListPlot[tTroot,PlotRange->{0.0,1},PlotLabel->"<s> vs kT",AxesLabel->{"kT","<s>"}]

Sample code 8.3. MC code for Ising model to simulate the microscopic sampling of the system at a fixed temperature. Total energy and magnetisation as function of MCStep are calculated.
No averaaging is performed.

L = 10; (* number of spins along a grid line *)
k = 1; J = 1; mu = 1.; 

MCStep = L^2; 
nMClast = 50*MCStep; 
T = 2.2*(J/k);(*setting temperature in unit of J/k *)
(*Reference temp is T = 2.27;*)
 
Initial = "Random"; (*chose initial condition among  "Allup, Alldown, Random"*)

If[Initial == "Random",
 For[ix = 1, ix <= L, ix++, 
     For[iy = 1, iy <= L, iy++,
 s0 = 2*RandomInteger[] - 1; s[ix, iy] = s0; 
];]; 
 ];(*end if "Random"*)


If[Initial == "Alldown",
 s0 = -1;
 For[ix = 1, ix <= L, ix++, 
     For[iy = 1, iy <= L, iy++, 
s[ix, iy] = s0; 
]; ]; 
]; (*end if Alldown*)

If[Initial == "Allup", 
s0 = Plus[1]; 
For[ix = 1, ix <= L, ix++, 
     For[iy = 1, iy <= L, iy++,
 s[ix, iy] = s0;
 ]; ]; 
];(*end if Allup*) 


(*Assigning periodic boundary conditions*)
For[ix = 1, ix <= L, ix++,
 s[ix, L + 1] = s[ix, 1];  s[ix, 0] = s[ix, L]; 
];(*end for ix*)
For[iy = 1, iy <= L, iy++, 
s[0, iy] = s[L, iy]; s[L + 1, iy] = s[0, iy]; 
]; (*end for iy*)

(*summing up all spins to calculate total magnetisation*)
M = mu*Sum[s[ix, iy], {ix, 1, L}, {iy, 1, L}]; 

m = N[M/L^2]; (*calculating magnetisation per spin*)

tnMCm = {{0, m}}; (*storing initial m at timestep n=0*)

(*summing up the total spin-spin interaction energy. Be noted of the factor 1/2 *)
TE =(1/2) Sum[(-J)*(s[ix, iy]*s[ix, iy - 1] + 
      s[ix, iy]*s[ix, iy + 1] + s[ix, iy]*s[ix - 1, iy] + 
      s[ix, iy]*s[ix + 1, iy]), {ix, 1, L}, {iy, 1, L}]; 

tnMCTE[0] = {0, TE}; (*storing the initial TE at timestep n=0*)


tixiys[0] = Table[s[ix, iy], {ix, 1, L}, {iy, 1, L}]; (*storing the configuration of the spins at n=0*)

ArrayPlot[tixiys[0], Mesh -> True, ColorRules -> 
    {1 -> White, -1 -> Black}]; (*ArrayPlot the spin configuration at n=0*)


(*start the MC time-stepping*)
For[nMC = 1, nMC <= nMClast, nMC++, 
   ixRan = Random[Integer, {1, L}]; 
    iyRan = Random[Integer, {1, L}]; 
strial = -s[ixRan, iyRan]; 
    
(*imposing periodic boundary condition*)
For[iy = 1, iy <= L, iy++, 
s[L + 1, iy] = s[1, iy]; s[0, iy] = s[L, iy]; 
]; (*end for iy*)

For[ix = 1, ix <= L, ix++, 
     s[ix, L + 1] = s[ix, 1]; s[ix, 0] = s[ix, L]; 
]; (*end for ix*)
    (*end of imposing periodic boundary condition*)



TEtrial = TE - (-J)*(s[ixRan, iyRan]*s[ixRan, iyRan - 1] + 
        s[ixRan, iyRan]*s[ixRan, iyRan + 1] + s[ixRan, iyRan]*
         s[ixRan - 1, iyRan] + s[ixRan, iyRan]*
         s[ixRan + 1, iyRan]) + 
      (-J)*(strial*s[ixRan, iyRan - 1] + 
        strial*s[ixRan, iyRan + 1] + strial*s[ixRan - 1, iyRan] + 
        strial*s[ixRan + 1, iyRan]); 

If[TEtrial < TE, 
     s[ixRan, iyRan] = strial; TE = TEtrial; 
, (*else if *)
     {DeltaE = TEtrial - TE; q = RandomReal[]; 
       Boltzmanfactor = Exp[-DeltaE/(k*T)]; 
       If[Boltzmanfactor > q, s[ixRan, iyRan] = strial; 
         TE = TEtrial; , {}]; }];(*end if*)
 
    M = mu*Sum[s[ix, iy], {ix, 1, L}, {iy, 1, L}]; 
m = N[M/L^2]; 
    
tnMCm = Append[tnMCm, {nMC, m}]; tnMCTE[nMC] = {nMC, TE};
    tixiys[nMC] = Table[s[ix, iy], {ix, 1, L}, {iy, 1, L}]; ]; 

(* Arrayplot the current spin configuration at timestep n *)
Manipulate[ArrayPlot[tixiys[nMC], Mesh -> True, 
   ColorRules -> {1 -> White, -1 -> Black}], {nMC, 0, nMClast, 1}]

ListPlot[Table[tnMCTE[nMC], {nMC, 0, nMClast}], 
  PlotLabel -> "Total Energy vs. MCStep", 
PlotRange->All,
  PlotStyle -> {AxesLabel -> {"MCStep", "Total Energy"}, 
    PointSize -> 0.001}]

ListPlot[tnMCm, PlotRange -> {-1.01, 1.01}, 
  PlotLabel -> "Total Magnetisation vs. MCStep", 
  PlotStyle -> {AxesLabel -> {"MCStep", "Total Magnetisation"}, 
    PointSize -> 0.001}]

Sample code 8.4. Next, we modify the code 8.3 so that it can calculate the average energy, average magentisation and specific heat capacity at a fixed T. Start sampling of data for averaging only after a threshold time step, nMCthr.

nMCthr=100*MCStep;aveTE=0;sumTE=0;tnMCaveTE={};avem=0;summ=0;tnMCavem={};tnMCC={};tnMCaveTEsq={};sumTEsq=0;

L = 4;
k = 1; J = 1; mu = 1.; 
MCStep = L^2; 
nMClast = nMCthr+200*MCStep; 
T = 3.5*(J/k); 
Initial = "Allup"; 

If[Initial == "Random", For[ix = 1, ix <= L, ix++, 
     For[iy = 1, iy <= L, iy++, s0 = 2*RandomInteger[] - 1; 
        s[ix, iy] = s0; ]; ]; ]; 

If[Initial == "Alldown", s0 = -1; For[ix = 1, ix <= L, ix++, 
     For[iy = 1, iy <= L, iy++, s[ix, iy] = s0; ]; ]; ]; 

If[Initial == "Allup", s0 = Plus[1]; For[ix = 1, ix <= L, ix++, 
     For[iy = 1, iy <= L, iy++, s[ix, iy] = s0; ]; ]; ]; 


For[ix = 1, ix <= L, ix++, s[ix, L + 1] = s[ix, 1]; 
    s[ix, 0] = s[ix, L]; ]; 
For[iy = 1, iy <= L, iy++, s[0, iy] = s[L, iy]; 
    s[L + 1, iy] = s[0, iy]; ]; 


M = mu*Sum[s[ix, iy], {ix, 1, L}, {iy, 1, L}]; 
m = N[M/L^2]; 
tnMCm = {{0, m}}; 
tnMCM = {{0, M}}; 
TE = (1/2)Sum[(-J)*(s[ix, iy]*s[ix, iy - 1] + 
      s[ix, iy]*s[ix, iy + 1] + s[ix, iy]*s[ix - 1, iy] + 
      s[ix, iy]*s[ix + 1, iy]), {ix, 1, L}, {iy, 1, L}]; 
tnMCTE[0] = {0, TE}; 
tixiys[0] = Table[s[ix, iy], {ix, 1, L}, {iy, 1, L}]; 
ArrayPlot[tixiys[0], Mesh -> True, ColorRules -> 
    {1 -> White, -1 -> Black}]; 
TE; 
tnMCTE[0]



(*Begin MC Stepping*)
For[nMC = 1, nMC <= nMClast, nMC++, 
   ixRan = Random[Integer, {1, L}]; 
    iyRan = Random[Integer, {1, L}]; strial = -s[ixRan, iyRan]; 
    
(*boundary condition*)
For[iy = 1, iy <= L, iy++, s[L + 1, iy] = s[1, iy]; 
      s[0, iy] = s[L, iy]; 
]; (*end for iy*)

For[ix = 1, ix <= L, ix++, 
     s[ix, L + 1] = s[ix, 1]; s[ix, 0] = s[ix, L]; 
]; (*end for ix*)
(*end of boundary condition*)

   TEtrial = TE - (-J)*(s[ixRan, iyRan]*s[ixRan, iyRan - 1] + 
        s[ixRan, iyRan]*s[ixRan, iyRan + 1] + s[ixRan, iyRan]*
         s[ixRan - 1, iyRan] + s[ixRan, iyRan]*
         s[ixRan + 1, iyRan]) + 
      (-J)*(strial*s[ixRan, iyRan - 1] + 
        strial*s[ixRan, iyRan + 1] + strial*s[ixRan - 1, iyRan] + 
        strial*s[ixRan + 1, iyRan]); 

(* test the trial state if it is to be accepted *)
If[TEtrial < TE, 
     s[ixRan, iyRan] = strial; TE = TEtrial; (*trial state is accepted*)
M=M+2*mu*strial;
, (*else*)(*trial state is to be further tested *)
     {DeltaE = TEtrial - TE; q = RandomReal[]; 
       Boltzmanfactor = Exp[-DeltaE/(k*T)]; 
       If[Boltzmanfactor > q, s[ixRan, iyRan] = strial; 
         TE = TEtrial; 
M=M+2*mu*strial;, {}
](*end if Boltzmanfactor > q,*);
 }];(*end if TEtrial < TE*)
(* end of test the trial state if it is to be accepted *)
(*
M = mu*Sum[s[ix, iy], {ix, 1, L}, {iy, 1, L}];*)
m = N[M/L^2]; 
(*tnMCM = Append[tnMCM, {nMC, M}]; *)
tnMCm = Append[tnMCm, {nMC, m}]; 
(*tnMCM = Append[tnMCM, {nMC, M}]; *)
tnMCTE[nMC] = {nMC, TE}; 
tixiys[nMC] = Table[s[ix, iy], {ix, 1, L}, {iy, 1, L}]; 

(*Averaging*)

If[nMC>nMCthr,
sumTE=sumTE+TE;
aveTE=sumTE/(nMC-nMCthr);
tnMCaveTE = Append[tnMCaveTE, {nMC, aveTE}]; 


sumTEsq=sumTEsq+TE^2;
aveTEsq=sumTEsq/(nMC-nMCthr);
tnMCaveTEsq = Append[tnMCaveTEsq, {nMC, aveTEsq}]; 

varE=aveTEsq-aveTE^2;
Cspec=varE/(k*T^2);
tnMCC = Append[tnMCC, {nMC,Cspec}]; 



summ=summ+m;
avem=summ/(nMC-nMCthr);
tnMCavem = Append[tnMCavem, {nMC, avem}]; 

];(* end if *)
(*End of Averaging*)

]; (*end for nMC *)


Manipulate[ArrayPlot[tixiys[nMC], Mesh -> True, 
   ColorRules -> {1 -> White, -1 -> Black}], {nMC, 0, nMClast, 1}];

ListPlot[Table[tnMCTE[nMC], {nMC, 0, nMClast}], PlotRange->All,
  PlotLabel -> {"TE vs. MCStep at kT =", k*T}, 
  PlotStyle -> {AxesLabel -> {"MCStep", "TE"}, 
    PointSize -> 0.001}]

ListPlot[tnMCaveTE, PlotRange -> All (*{-1.01, 1.01}*), 
  PlotLabel -> {"aveTE vs. MCStep at kT = ", k*T}, 
  PlotStyle -> {AxesLabel -> {"MCStep", "aveTE"}, 
    PointSize -> 0.001}]


ListPlot[tnMCaveTEsq, PlotRange -> All, 
  PlotLabel -> {"average TEsq vs. MCStep at kT = ",k*T}, 
  PlotStyle -> {AxesLabel -> {"MCStep", "TEsq"}, 
    PointSize -> 0.001}]


ListPlot[tnMCC, PlotRange -> All, 
  PlotLabel -> {"C vs. MCStep at kT = ",k*T}, 
  PlotStyle -> {AxesLabel -> {"MCStep", "TEsq"}, 
    PointSize -> 0.001}]


ListPlot[tnMCm, PlotRange -> All, 
  PlotLabel -> {"m vs. MCStep at kT = ",k*T}, 
  PlotStyle -> {AxesLabel -> {"MCStep", "Total Magnetisation"}, 
    PointSize -> 0.001}]

ListPlot[tnMCavem, PlotRange -> All, 
  PlotLabel -> {"average m vs. MCStep at kT = ",k*T}, 
  PlotStyle -> {AxesLabel -> {"MCStep", "Total Magnetisation"}, 
    PointSize -> 0.001}]

Print["kT =", k*T," Cspec="," ",Cspec, " avem="," ",avem];

Sample code 8.5. Sample code 8.4 is modified to calculate specific heat, magnetisation m as a function of temperature T. Basically, to do so we simply put the sample code 8.4 into a loop to vary for temperature.

nMCthr=100*MCStep;

nave=100;(*#, in unit of MSCstep, of step for sample averaging *)

tnMCaveTE={};tnMCavem={};tnMCC={};
tkTavem={};tkTCspec={};tnMCaveTEsq={};
tkTaveTE={};


Tini=0.001;Tfinal=5;nTStep=25;
DeltaT=(Tfinal-Tini)/nTStep;

L = 10;(*Grid size*)
k = 1; J = 1; mu = 1.; 
MCStep = L^2; 
nMClast = nMCthr+nave*MCStep; (*nMCthr MCStep for initial thermalisation, *)

(*T = 3.5*(J/k); *)
Initial = "Alldown"; 

If[Initial == "Random", For[ix = 1, ix <= L, ix++, 
     For[iy = 1, iy <= L, iy++, s0 = 2*RandomInteger[] - 1; 
        s[ix, iy] = s0; ]; ]; ]; 
If[Initial == "Alldown", s0 = -1; For[ix = 1, ix <= L, ix++, 
     For[iy = 1, iy <= L, iy++, s[ix, iy] = s0; ]; ]; ]; 
If[Initial == "Allup", s0 = Plus[1]; For[ix = 1, ix <= L, ix++, 
     For[iy = 1, iy <= L, iy++, s[ix, iy] = s0; ]; ]; ]; 

For[ix = 1, ix <= L, ix++, s[ix, L + 1] = s[ix, 1]; 
    s[ix, 0] = s[ix, L]; ]; 
For[iy = 1, iy <= L, iy++, s[0, iy] = s[L, iy]; 
    s[L + 1, iy] = s[0, iy]; ]; 

M = mu*Sum[s[ix, iy], {ix, 1, L}, {iy, 1, L}]; 
m = N[M/L^2]; 
tnMCm = {{0, m}}; 
tnMCM = {{0, M}}; 
TE = (1/2)Sum[(-J)*(s[ix, iy]*s[ix, iy - 1] + 
      s[ix, iy]*s[ix, iy + 1] + s[ix, iy]*s[ix - 1, iy] + 
      s[ix, iy]*s[ix + 1, iy]), {ix, 1, L}, {iy, 1, L}]; 
tnMCTE[0] = {0, TE}; 
tixiys[0] = Table[s[ix, iy], {ix, 1, L}, {iy, 1, L}]; 
ArrayPlot[tixiys[0], Mesh -> True, ColorRules -> 
    {1 -> White, -1 -> Black}]; 
TE; 
tnMCTE[0] 

(*Begin stepping temperature*)
For[iT=0,iT≤nTStep,iT++,
T=Tini+iT*DeltaT;

sumTEsq=0;aveTE=0;sumTE=0;
summ=0;avem=0;

M = mu*Sum[s[ix, iy], {ix, 1, L}, {iy, 1, L}]; 
m = N[M/L^2]; 
TE = (1/2)Sum[(-J)*(s[ix, iy]*s[ix, iy - 1] + 
      s[ix, iy]*s[ix, iy + 1] + s[ix, iy]*s[ix - 1, iy] + 
      s[ix, iy]*s[ix + 1, iy]), {ix, 1, L}, {iy, 1, L}]; 

(*Begin MC Stepping*)


If[iT==0,nMClast = nMCthr+nave*MCStep];
If[iT>0,nMClast = nMCthr+nave*MCStep];(*after iT≥1, thermalisation has been achieved hence take less # of step for thermalisation purpose *)

For[nMC = 1, nMC <= nMClast, nMC++, 
   ixRan = Random[Integer, {1, L}]; 
    iyRan = Random[Integer, {1, L}]; strial = -s[ixRan, iyRan]; 
    
(*boundary condition*)
For[iy = 1, iy <= L, iy++, s[L + 1, iy] = s[1, iy]; 
      s[0, iy] = s[L, iy]; 
]; (*end for iy*)

For[ix = 1, ix <= L, ix++, 
     s[ix, L + 1] = s[ix, 1]; s[ix, 0] = s[ix, L]; 
]; (*end for ix*)
(*end of boundary condition*)

(*Mnew=M;*)
    TEtrial = TE - (-J)*(s[ixRan, iyRan]*s[ixRan, iyRan - 1] + 
        s[ixRan, iyRan]*s[ixRan, iyRan + 1] + s[ixRan, iyRan]*
         s[ixRan - 1, iyRan] + s[ixRan, iyRan]*
         s[ixRan + 1, iyRan]) + 
      (-J)*(strial*s[ixRan, iyRan - 1] + 
        strial*s[ixRan, iyRan + 1] + strial*s[ixRan - 1, iyRan] + 
        strial*s[ixRan + 1, iyRan]); 

(* test the trial state if it is to be accepted *)
If[TEtrial < TE, 
     s[ixRan, iyRan] = strial; TE = TEtrial; (*trial state is accepted*)
M=M+2*mu*strial;
, (*else*)(*trial state is to be further tested *)
     {DeltaE = TEtrial - TE; q = RandomReal[]; 
       Boltzmanfactor = Exp[-DeltaE/(k*T)]; 
       If[Boltzmanfactor > q, s[ixRan, iyRan] = strial; 
         TE = TEtrial; 
M=M+2*mu*strial;, {}
](*end if Boltzmanfactor > q,*);
 }];(*end if TEtrial < TE*)
(* end of test the trial state if it is to be accepted *)
(*
M = mu*Sum[s[ix, iy], {ix, 1, L}, {iy, 1, L}];*)
m = N[M/L^2]; 
(*tnMCM = Append[tnMCM, {nMC, M}]; *)
tnMCm := Append[tnMCm, {nMC, m}]; 
(*tnMCM = Append[tnMCM, {nMC, M}]; *)
tnMCTE[nMC] := {nMC, TE}; 
tixiys[nMC] := Table[s[ix, iy], {ix, 1, L}, {iy, 1, L}]; 

(*Averaging*)

If[nMC>nMCthr,
sumTE=sumTE+TE;
aveTE=sumTE/(nMC-nMCthr);
tnMCaveTE := Append[tnMCaveTE, {nMC, aveTE}]; 


sumTEsq=sumTEsq+TE^2;
aveTEsq=sumTEsq/(nMC-nMCthr);
tnMCaveTEsq := Append[tnMCaveTEsq, {nMC, aveTEsq}]; 

varE=aveTEsq-aveTE^2;
Cspec=varE/(k*T^2);
tnMCC := Append[tnMCC, {nMC,Cspec}]; 




summ=summ+m;
avem=summ/(nMC-nMCthr);
tnMCavem := Append[tnMCavem, {nMC, avem}]; 

];(* end if *)
(*End of Averaging*)
]; (*end for nMC *)
Print["kT =", k*T," Cspec="," ",Cspec, " avem="," ",avem, " aveTE/L^2 ", N[aveTE/(L^2)]];
tkTavem=Append[tkTavem,{k*T,avem}];
(*Print["iT=",iT, " tkTavem",tkTavem];*)
tkTCspec=Append[tkTCspec,{k*T,Cspec}];
tkTaveTE=Append[tkTaveTE,{k*T,aveTE/(L^2)}];
(*Print["iT=",iT, " tkTCspec",tkTCspec];*)

];(*end for iT*)

ListPlot[tkTavem, PlotRange -> All,Mesh->True,(*Joined->True,*)
  PlotLabel -> {"avem vs. kT "}, 
  PlotStyle -> {AxesLabel -> {"kT", "avem"}, 
    PointSize -> 0.01}]

ListPlot[tkTCspec, PlotRange -> All,Mesh->True,(*Joined->True,*)
  PlotLabel -> {"Cspec vs. kT "}, 
  PlotStyle -> {AxesLabel -> {"kT", "Cspec"}, 
    PointSize -> 0.01}]

Sample code 8.6 : Now we include the magnetic field into the MC code to calculate magnetisation as a function of magnetic field at fixed temperature.The magnetic field comes into the picture through the definition of the total energy TE = TE (spin - spin) + TE (magnetic). Below T_c it is show that found that magnetisation m undergoes a first order phase transition as the external magnetic field swips through the value distinctly different from zero.

nMCthr=10*MCStep;

nave=20;(*#, in unit of MSCstep, of step for sample averaging *)

tnMCaveTE={};tnMCavem={};tnMCC={};
tBavem={};tBCspec={};tnMCaveTEsq={};
tBaveTE={};


Bini=-5;Bfinal=5;nBStep=40;
DeltaB=(Bfinal-Bini)/nBStep;
T = 5; (*default T=1*)

L = 10;(*Grid size*)
k = 1; J = 1;(* J = 0; *)
mu = 1.; 
MCStep = L^2; 
nMClast = nMCthr+nave*MCStep; (*nMCthr MCStep for initial thermalisation, *)


Initial = "Random"; 
If[Initial == "Random", For[ix = 1, ix <= L, ix++, 
     For[iy = 1, iy <= L, iy++, s0 = 2*RandomInteger[] - 1; 
        s[ix, iy] = s0; ]; ]; ]; 
If[Initial == "Alldown", s0 = -1; For[ix = 1, ix <= L, ix++, 
     For[iy = 1, iy <= L, iy++, s[ix, iy] = s0; ]; ]; ]; 
If[Initial == "Allup", s0 = Plus[1]; For[ix = 1, ix <= L, ix++, 
     For[iy = 1, iy <= L, iy++, s[ix, iy] = s0; ]; ]; ]; 
For[ix = 1, ix <= L, ix++, s[ix, L + 1] = s[ix, 1]; 
    s[ix, 0] = s[ix, L]; ]; 
For[iy = 1, iy <= L, iy++, s[0, iy] = s[L, iy]; 
    s[L + 1, iy] = s[0, iy]; ]; 
M = mu*Sum[s[ix, iy], {ix, 1, L}, {iy, 1, L}]; 
m = N[M/L^2]; 
tnMCm = {{0, m}}; 
tnMCM = {{0, M}}; 
(*tnMCMnew = {{0, M}}; *)
TE =((1/2) Sum[(-J)*(s[ix, iy]*s[ix, iy - 1] + 
      s[ix, iy]*s[ix, iy + 1] + s[ix, iy]*s[ix - 1, iy] + 
      s[ix, iy]*s[ix + 1, iy]), {ix, 1, L}, {iy, 1, L}])-Bini*M; 

tnMCTE[0] = {0, TE}; 
tixiys[0] = Table[s[ix, iy], {ix, 1, L}, {iy, 1, L}]; 
ArrayPlot[tixiys[0], Mesh -> True, ColorRules -> 
    {1 -> White, -1 -> Black}]; 
TE; 
tnMCTE[0]; 

(*Begin stepping temperature*)
For[iB=0,iB≤nBStep,iB++,
B=Bini+iB*DeltaB;

sumTEsq=0;aveTE=0;sumTE=0;
summ=0;avem=0;

M = mu*Sum[s[ix, iy], {ix, 1, L}, {iy, 1, L}]; 
m = N[M/L^2]; 
TE =( (1/2)Sum[(-J)*(s[ix, iy]*s[ix, iy - 1] + 
      s[ix, iy]*s[ix, iy + 1] + s[ix, iy]*s[ix - 1, iy] + 
      s[ix, iy]*s[ix + 1, iy]), {ix, 1, L}, {iy, 1, L}]) -mu*B*M; 

(*Begin MC Stepping*)


If[iB==0,nMClast = nMCthr+nave*MCStep];
If[iB>0,nMClast = nMCthr+nave*MCStep];(*after iT≥1, thermalisation has been achieved hence take less # of step for thermalisation purpose *)

For[nMC = 1, nMC <= nMClast, nMC++, 
   ixRan = Random[Integer, {1, L}]; 
    iyRan = Random[Integer, {1, L}]; strial = -s[ixRan, iyRan]; 
    
(*boundary condition*)
For[iy = 1, iy <= L, iy++, s[L + 1, iy] = s[1, iy]; 
      s[0, iy] = s[L, iy]; 
]; (*end for iy*)

For[ix = 1, ix <= L, ix++, 
     s[ix, L + 1] = s[ix, 1]; s[ix, 0] = s[ix, L]; 
]; (*end for ix*)
(*end of boundary condition*)

(*Mnew=M;*)
    TEtrial = TE - (-J)*(s[ixRan, iyRan]*s[ixRan, iyRan - 1] + 
        s[ixRan, iyRan]*s[ixRan, iyRan + 1] + s[ixRan, iyRan]*
         s[ixRan - 1, iyRan] + s[ixRan, iyRan]*
         s[ixRan + 1, iyRan]) + 
      (-J)*(strial*s[ixRan, iyRan - 1] + 
        strial*s[ixRan, iyRan + 1] + strial*s[ixRan - 1, iyRan] + 
        strial*s[ixRan + 1, iyRan])+2mu*B*s[ixRan,iyRan]; 

(* test the trial state if it is to be accepted *)
If[TEtrial < TE, 
     s[ixRan, iyRan] = strial; TE = TEtrial; (*trial state is accepted*)
M=M+2*mu*strial;
, (*else*)(*trial state is to be further tested *)
     {DeltaE = TEtrial - TE; q = RandomReal[]; 
       Boltzmanfactor = Exp[-DeltaE/(k*T)]; 
       If[Boltzmanfactor > q, s[ixRan, iyRan] = strial; 
         TE = TEtrial; 
M=M+2*mu*strial;, {}
](*end if Boltzmanfactor > q,*);
 }];(*end if TEtrial < TE*)
(* end of test the trial state if it is to be accepted *)
(*
M = mu*Sum[s[ix, iy], {ix, 1, L}, {iy, 1, L}];*)
m = N[M/L^2]; 
(*tnMCM = Append[tnMCM, {nMC, M}]; *)
tnMCm := Append[tnMCm, {nMC, m}]; 
(*tnMCM = Append[tnMCM, {nMC, M}]; *)
tnMCTE[nMC] := {nMC, TE}; 
tixiys[nMC] := Table[s[ix, iy], {ix, 1, L}, {iy, 1, L}]; 

(*Averaging*)

If[nMC>nMCthr,
(*sumTE=sumTE+TE;
aveTE=sumTE/(nMC-nMCthr);
tnMCaveTE = Append[tnMCaveTE, {nMC, aveTE}]; *)
(*
sumTEsq=sumTEsq+TE^2;
aveTEsq=sumTEsq/(nMC-nMCthr);
tnMCaveTEsq := Append[tnMCaveTEsq, {nMC, aveTEsq}]; 
*)

(*varE=aveTEsq-aveTE^2;
Cspec=varE/(k*T^2);
tnMCC := Append[tnMCC, {nMC,Cspec}]; 
*)

summ=summ+m;
avem=summ/(nMC-nMCthr);
tnMCavem := Append[tnMCavem, {nMC, avem}]; 
];(* end if *)
(*End of Averaging*)
]; (*end for nMC *)
(*Print["B =", B," Cspec="," ",Cspec, " avem="," ",avem, " aveTE/L^2 ", N[aveTE/(L^2)]];*)
tBavem=Append[tBavem,{B,avem}];
(*Print["iT=",iT, " tkTavem",tkTavem];*)
(*tBCspec=Append[tBCspec,{B,Cspec}];
tBaveTE=Append[tBaveTE,{B,aveTE/(L^2)}];*)
Print["iB=",iB, " avem  ",avem];

];(*end for iT*)

ListPlot[tBavem, PlotRange -> All,Mesh->True,
  PlotLabel -> {"avem vs. B "}, 
  PlotStyle -> {AxesLabel -> {"B", "avem"}, 
    PointSize -> 0.01}]