clear all;
n =10; % number of atoms in a row
n_temp = 20; % number of steps for temperature loop
nbin = 5; % number of bins to average energyies etc
nmc = 500; % number of monte carlo steps for one bin
T = 0.9; % initial temperature
JJ = -1.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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rand('state', 1)
for mm=1: n_temp % loop over temperature
T = T + 0.1;
T
beta = 1/T; temper(mm)= T;
for kub=1:nbin % nbin loop to reduce variance
% kub
% 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 '
