% Demonstrates computation of energy function for collection
% of spins defined on two-dimensional lattice.

% Number of spins per side ...
n = 5

% Coupling constant ...
Jij = 1

% Generate random values of spins, +/- 1 ...
s = rand(n);
for I =  1 : n
   for J = 1 : n
      if s(I,J) > 0.5 
         s(I,J) = 1;
      else
         s(I,J) = -1;
      end
   end
end
s

% For every spin arranged in the 2D lattice ...
for I = 1 : n
   for J = 1 : n
      % Compute energy with spin in original state ...
      E0 = 0;
      for K = 1 : n
         for L = 1 : n
            for M = 1 : n
               for N = 1 : n
                  E0 = E0 + Jij * s(K,L) * s(M,N);
               end
            end
         end
      end
      E0 = -E0
      % Compute energy with spin in flipped state ...
      s(I,J) = -s(I,J);
      E1 = 0;
      for K = 1 : n
         for L = 1 : n
            for M = 1 : n
               for N = 1 : n
                  E1 = E1 + Jij * s(K,L) * s(M,N);
               end
            end
         end
      end
      E1 = -E1
      % If energy increased, flip spin back to its original state ...
      if E1 > E0
         s(I,J) = -s(I,J);
      end
      s
   end
end
