clear all;
n =3; nnn = 1; nbin = 200;
T =4.; JJ = -1.0; nmc = 1000; beta = 1/T;
rand('state', 0)
%----------- generate random spins for the first time ----
for ip= 2: n+1
for jp= 2: n+1
rr = 2.*rand - 1.0 ;
s(ip,jp) = sign(rr); %s(ip,jp) = 1;
end
end
s
%------------- fix the boundary conditions as periodic ------
for ja=2: n+1
s(1,ja) = s(n+1,ja); s(n+2,ja) = s(2, ja); end
for ia=2: n+1
s(ia,1)= s(ia,n+1); s(ia,n+2) = s(ia,2); end
% ------------------------------------------------------------
%i = 2; j = 2;
for j = 2: n+1 % loop over all mesh points
for i = 2: n+1 % loop over all mesh points
en_initial = JJ*s(i,j)*(s(i+1,j)+s(i-1,j)+s(i,j+1)+s(i,j-1));
s_old = s(i,j); % save this for possible reverting
if(s(i,j)== +1) snew = -1; end
if(s(i,j)== -1) snew = +1; end
s(i,j) = snew; %(spin flip has been done)
% refresh the boundary values for the flip
for jj=2: n+1
s(1,jj) = s(n+1,jj); s(n+2,jj) = s(2, jj); end
for ii=2: n+1
s(ii,1) = s(ii,n+1); s(ii,n+2) = s(ii,2); end
en_final = JJ*s(i,j)*(s(i+1,j)+s(i-1,j)+s(i,j+1)+s(i,j-1));
delta_e = en_final - en_initial;
if(delta_e<=0) s(i,j) = snew; % keep the flip
'spin is flipped; as delta_e<=0 '
end
if(delta_e>0)
rr = rand;
if(rr<= exp(-beta*delta_e))
s(i,j)= snew; % again keep the flip
'spin is flipped when delta_e>0'
elseif(rr>exp(-beta*delta_e))
s(i,j) = s_old; % remove the flip
'spin is not flipped'
end
end
% refresh the boundary values after the flip -------------
for jk=2: n+1
s(1,jk) = s(n+1,jk); s(n+2,jk) = s(2, jk); end
for ik=2: n+1
s(ik,1) = s(ik,n+1); s(ik,n+2) = s(ik,2); end
s
