(* 1D Ising Model *)
(* Initialisation *)
J = 1.0; k=1.0;mu=1.0;initconfig=0;
L = 20;nmax=5*L;
T=25.1*(J/k);(*setting temperature in unit of J/k*)(*Reference temp is T=2.27;*)

(* Prepare initial config. *)
For[i=1,i≤L,i++,
If[initconfig==0,s[i]=1];
If[initconfig==1,s[i]=-1];
If[initconfig==2,s[i]=s[i]=2RandomInteger[]-1;];
];(*end for i*)
ts[0]=Table[s[i],{i,1,L}];
(* periodic BC*)
s[0]=s[L];s[L+1]=s[1];
(* end of periodic BC*)
TE=(-J/2.0)*Sum[s[i]*(s[i-1]+s[i+1]),{i,1,L}];
Print["initially, TE = ",TE];
(*visualisation*)
ArrayPlot[{ts[0]},Mesh->True,ColorRules->{1->White,-1->Black}];
(*end of visualisation*)


tnTE={};tnm={};

For[n=1,n≤nmax,n++,
Clear[iran,st];
iran=RandomInteger[{1,L}];
strial=-s[iran];
(* periodic BC*)
s[0]=s[L];s[L+1]=s[1];
(* end of periodic BC*)

(*TEtrial=TE-(-J*(s[iran]*s[iran-1]+s[iran]*s[iran+1]))+(-J*strial*(s[iran-1]+s[iran+1]));*)

TEtrial=TE+2J(s[iran]*s[iran-1]+s[iran]*s[iran+1]);

(*Print["n,TEtrial= ",{n,TEtrial}];*)

Eflip=TEtrial-TE;
Print["n, iran, Eflip, TE ",{n,iran, Eflip,TE}];

If[Eflip<0,s[iran]=strial,TE=TEtrial
,(*else*)
If[RandomReal[]<Exp[-Eflip/(k*T)],
s[iran]=strial;TE=TEtrial
];(*end if *)
];(*end if Eflip*)

(*Print["n, k*T , TE, TEtrial,RandomReal[],Exp[-Eflip/(k*T)]",{n,k*T,TE,TEtrial,RandomReal[],Exp[-Eflip/(k*T)]}];*)

(*update the spins*)
ts[n]=Table[s[i],{i,1,L}];
(*end of update the spins*)

tnTE=Append[tnTE,{n,TE}];
tnm=Append[tnm,{n,Sum[s[i],{i,1,L}]}/L];

];(*end for n*)

ListPlot[tnTE,PlotLabel->"TE vs n"]
ListPlot[tnm,PlotLabel->"m vs n"]

Manipulate[ArrayPlot[{ts[n]},Mesh->True,ColorRules->{1->White,-1->Black}],{n,0,nmax,1}]

For[n=1,n≤5,n++,
iran=RandomInteger[{1,10}];
Print[iran]
]