clear all;
n =10; nnn = 20; nbin = 10;
T =1.0; JJ = -1.0; nmc = 1000;
for jk=1: nnn
st_dev_mg(jk)=0.0; end
rand('state', 1)
for mm=1: nnn % loop over temperature
mm
T = T + 0.1;
T
beta = 1/T; temper(mm)= T;
for kub=1:nbin % nbin loop to reduce variance
kub
%%%% 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
% 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
for kkk=1: nmc % loop over changes in spins
sum_en = 0.0; sum_mg = 0.0;
for j = 2: n+1 % loop over all mesh points
for i = 2: n+1 % loop over all mesh points
% Start: Metropolis Algorithm for spin flip -----
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 possiblereverting
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))
