understand a code to simulate ising model in 2D using montecarlo method

Hello everyone,

I found this code in Russian book , to simulate ising model in 2D using montecarlo method, but Franky I don't understand the 22 lines , alhtough it gives 2d square grid with a vector in each site, would anyone have an idea about this code ? I will be highly grateful for any remarks or feedback

Thanks in advance

with best

-------------------------- The code is :

clear all 
close all 
Nspin=49; % number of spins. Must be square of integer 
J=2; %  constant of interaction 
h=0; % external magnetic field 
Esi=10; % energy of system (microcanonic distribution) 
NTrial=1000; % number of trials 
[Es,Ed,SpM,A,S]=Ising2(Nspin,J,h,Esi,NTrial); 
% plot of energy of system 
% plot of energy of daemon 
[~,m]=size(Ed); % [~,m], that means that you just want the second output of your function, and do not care the first one. 
% is equal to Nspin*Ntrial and considered as time 
i=1:m; 
 figure(1);plot(i,Es(i)) 
 figure(2);plot(i,Ed(i)) 
%% 

% plot of spin coniguration 
m1=floor(m/3); 
m2=2*m1; 
m3=3*m1; 
    i=1:sqrt(Nspin); 
    j=1:sqrt(Nspin); 
            V(i,j)=0; 
            U(i,j)=0; 
            Z(i,j)=0; 
            figure 
            subplot(4,1,1);quiver3(Z,U,V,S(:,:,1)); colormap white; 
            title('Spins initial configuration') 
            subplot(4,1,2);quiver3(Z,U,V,S(:,:,m1)); colormap white 
            title(['Spins configuration for time ', num2str(m1)]) 
            subplot(4,1,3); quiver3(Z,U,V,S(:,:,m2)); colormap white 
            title(['Spins configuration for time ', num2str(m2)]) 
            subplot(1,1,1); quiver3(Z,U,V,S(:,:,m)); colormap white 
             title(['Spins configuration for time ', num2str(m)]) 

            mean(Es) % mean energy of system
----------------
Ising2.m 
function [Es,Ed,SpM,A,S] = Ising2(Nspin,J,h,Esi,NTrial)
% input parameters are described in square_grid1.m
% output parameters
% Es is instantaneous energy of system after each step of walk in each trial
%SpM,A,S are spin configuration
Ns=Nspin.^0.5; % number of spin for one side of square
s=ones(Ns,Ns); % initial spin configuration. All spin projections are 1
Esystem=-(J+h)*Nspin; % initial system energy
Edemon=4*J*floor((Esi-Esystem)/(4*J)); % initial energy of daemon.
% on each trial daemon energy is compared with energy of all spin knots of grid
% and energy of knot and daemon is changed not changed according some
% condition
Es(1)=Esystem; 
Ed(1)=Edemon;
S=s;
k=1;
for i=1:NTrial
  Accept=0;
  for j=1:Nspin
     % random choice of knot
     Ix=floor(Ns*rand(1)+1); 
     Iy=floor(Ns*rand(1)+1);
     % border conditions
     if Ix==1
        Left=Ns;
     else
        Left=Ix-1;
     end;
     if Ix==Ns
       Right=1;
     else;
        Right=Ix+1;
     end;
     if Iy==1
        Down=Ns;
     else;  
        Down=Iy-1;
     end;
     if Iy==Ns
        Up=1;
     else
        Up=Iy+1;
     end;
     % energy change
     de=2*s(Iy,Ix)*(-h+J*(s(Iy,Left)+s(Iy,Right)+s(Down,Ix)+s(Up,Ix)));
     if de<=Edemon % energy change accepted
        s(Iy,Ix)=-s(Iy,Ix);
        Accept=Accept+1;
        Edemon=Edemon-de;
        Esystem=Esystem+de;
     end;
     k=k+1;
     Es(k)=Esystem;
     Ed(k)=Edemon;
     A(k-1)=Accept;
     s1=sum(s);
     SpM(k)=sum(s1);
     S=cat(3,S,s);
  end;
end;
A=A/NTrial;

[Alan Stevens]

The heart of the calculation will lie in the function Ising2, which you haven't listed!

What aspect don't you understand? The Matlab syntax; Ising models, ...?

I have a two-dimensional Ising model MATAB program that I wrote many years ago that uses the Metropolis algorithm. I would be happy to let you have it, but it might not be any more meaningful to you than your Russian code. I can judge that better if you answer the question above. My code will be different from the Russian version as mine calculates the equilibrium configuration only, whereas the comments in the coding above hints at time dependence (I think!).

[Julia Dulie]

I modified the question , I added the code ising2.m in the end

Of course your code will be useful, as I am learning to do montecarlo using matlab with different approches, so I will be grateful if you share your code also

thank you very much Alan Stevens

[Alan Stevens]

Ok, here is my code. I've adjusted the values of J (the exchange energy), and L (L^2=Nspins) to match that of your Russian code (my original had J = 1 and L = 100) - of course, you can change them to whatever you like. This was written over 10 years ago, so there may well be more efficient ways of doing things these days!

.

% ISING.M
%
% Two-dimensional model of an Ising-spin ferromagnet.
% Calculates energies, magnetisations and heat capacities of
% equilibrium states.

rng(1,'twister');

% Basic data
L = 7;            % Number of spin sites on a side.
N = L^2;            % Total number of spin sites.
J = 2;              % Exchange constant +(-)ve for ferro(antiferro)magnet.
T = 5;              % Temperature (strictly, this is kT)
%T = 2/log(1+sqrt(2)); % Theoretical Tc (Curie temperature) for infinite system
% Monte-Carlo data
eqstart = 2000;          % Number of sweeps to establish equilibrium.
neq     = 1000;          % Number of equilibrium sweeps.
sweeps  = eqstart + neq; % Total number of sweeps.

% Initialisation
Spins  = -1+2*round(rand(L,L)); % Matrix of spins +/-1 at random.
M = zeros(neq,1);               % For recording equilibrium magnetisation.
E = zeros(neq,1);               % For recording equilibrium energies. 
list = (1:N)';                  % Vector of integers 1 to N in order.
tic
% Perform Monte-Carlo trials (sweeps)
for sweep = 1:sweeps
    
    % Randomise order of rows and columns
    n = ceil(N*rand(N,1));           % Random integers in the range 1 to N
    for i = 1:N                      % -----------------------------------
        tmp         = list(i);       %                               
        list(i)     = list(n(i));    %    Shuffle the integers
        list(n(i))  = tmp;           %                            
    end                              % -----------------------------------
    % Convert the list to lists of row and column indices
    [rowlist, collist] = ind2sub([L,L], list);
    
    % Flip spins using the Metropolis algorithm.
    % For all rows get the numbers of the row and neighbouring rows (using
    % cyclic boundary conditions).  Similarly for columns.
    % Sum the spins on the neighbour sites, calculate the energy change
    % needed to flip a spin, and use the Metropolis algorithm to decide
    % if the spin should be flipped.
    for counter = 1:N                      
        r   = rowlist(counter);
        lo  = r-1; if lo<1, lo = L; end     
        hi  = r+1; if hi>L, hi = 1; end
        c   = collist(counter);
        left  = c-1; if left<1,   left = L; end
        right = c+1; if right>L, right = 1; end                       
        % Sum neighbouring spins.
        SS = Spins(lo,c)+Spins(hi,c) ...
            +Spins(r,left)+Spins(r,right);
        % Calculate energy change if spin were to be flipped.
        DeltaE = 2*J*SS*Spins(r,c);
        % Metropolis algorithm.
        if DeltaE<0, Spins(r,c) = -Spins(r,c);
        elseif rand<exp(-DeltaE/T), Spins(r,c) = -Spins(r,c);
        end
    end % of counter
   
    % Start recording energies and magnetisations.
    % For all sites determine nearest neighbour sites. 
    % Sum the neighbour spins.  Accumulate the energies.
    if sweep>eqstart    
        % Calculate overall energy
        Energy = 0;
        for r = 1:L
            lo  = r-1; if lo<1, lo = L; end                 % nearest
            hi  = r+1; if hi>L, hi = 1; end                 % rows
            for c = 1:L
                left  = c-1; if left<1,   left = L; end     % nearest
                right = c+1; if right>L, right = 1; end     % columns
                % Sum neighbouring spins.
                SS = Spins(lo,c)+Spins(hi,c)+Spins(r,left)+Spins(r,right);
                % Calculate energy.
                Energy = Energy - 0.5*J*SS*Spins(r,c);
            end % of c loop
        end % of r loop
        % Calculate and store energy and magnetisation of the current state.
        t = sweep-eqstart;              % Equilibrium counter.
        E(t) = Energy/N;                % Total Energy (per spin)
        M(t) = sum(Spins(:))/N;         % Net magnetisation (per spin).
    end % of if sweep>eqstart
  
end % of sweeps
toc
% Calculate average equilibrium energy, heat capacity and magnetisation,
% all per spin.
Eave = mean(E);   
Cv = N*var(E)/T^2;
Mave = mean(M); 

varM = var(M);
varE = var(E);
% Magnetic susceptibility.
Chi = N*varM./T;
% Display results
disp(['T = ' num2str(T)])
disp(['E = ' num2str(Eave) ' sdE = ' num2str(sqrt(varE)) '    Cv = ' num2str(Cv)])
disp(['M = ' num2str(Mave) ' sdM = ' num2str(sqrt(varM)) '   Chi = ' num2str(Chi)])
% Plot final spin state
figure
surf(Spins),view(2),axis off
title(['Temperature = ' num2str(T)])

figure
subplot(2,1,1)
plot(E)
title(['Temperature = ' num2str(T)])
v = axis;
v(3) = -2; v(4) = 0;
axis(v);
ylabel('Energy')
subplot(2,1,2)
plot(M)
v = axis;
v(3) = -1; v(4) = 1;
axis(v)
xlabel('sweeps')
ylabel('Magnetisation')

[Julia Dulie]

thank you very much Alan Stevens , you are very nice,

even if it is old , it is a good starting point as tutorial

wish you have a nice day

[Alan Stevens]

The Ising model calculates the energy, E, of a number of spins on a lattice (2-dimensional here) using:

where J is the "exchange energy" and S is a spin either up or down (+1 or -1), and the summations are over nearest neighbours. The Russian code seems to use what is known as the Gibbs algorithm. Basically, it calculates the change in energy at a random lattice site that woud occur if the spin on that site changed from up to down (or vice versa). If the change in energy is negative then it accepts it and updates the overall energy. It repeats this at random sites a large number of times, eventually driving towards equilibrium (where the energy might fluctuate, but stays the same on average).

I noticed some errors in the Russian code. Here's a working version:

Esi=10; % energy of system (microcanonic distribution) 
NTrial=1000; % number of trials 
[Es,Ed,SpM,A,S]=Ising2(Nspin,J,h,Esi,NTrial); 
% plot of energy of system 
% plot of energy of daemon 
[~,m]=size(Ed); % [~,m], that means that you just want the second output of your function, and do not care the first one. 
% is equal to Nspin*Ntrial and considered as time 
i=1:m; 
 figure(1);plot(i,Es(i)) 
 figure(2);plot(i,Ed(i)) 
%% 
% plot of spin coniguration 
m1=floor(m/3); 
m2=2*m1; 
m3=3*m1; 
    i=1:sqrt(Nspin); 
    j=1:sqrt(Nspin); 
            V(i,j)=0; 
            U(i,j)=0; 
            Z(i,j)=0; 
            figure 
            subplot(4,1,1);quiver3(Z,U,V,S(:,:,1)); colormap white; 
            title('Spins initial configuration') 
            subplot(4,1,2);quiver3(Z,U,V,S(:,:,m1)); colormap white 
            title(['Spins configuration for time ', num2str(m1)]) 
            subplot(4,1,3); quiver3(Z,U,V,S(:,:,m2)); colormap white 
            title(['Spins configuration for time ', num2str(m2)]) 
            subplot(4,1,4); quiver3(Z,U,V,S(:,:,m)); colormap white 
             title(['Spins configuration for time ', num2str(m)]) 
            mean(Es) % mean energy of system
%----------------
%Ising2.m 
function [Es,Ed,SpM,A,S] = Ising2(Nspin,J,h,Esi,NTrial)
% input parameters are described in square_grid1.m
% output parameters
% Es is instantaneous energy of system after each step of walk in each trial
%SpM,A,S are spin configuration
Ns=Nspin^0.5; % number of spin for one side of square
s=ones(Ns,Ns); % initial spin configuration. All spin projections are 1
Esystem=-(J+h)*Nspin; % initial system energy
Edemon=4*J*floor((Esi-Esystem)/(4*J)); % initial energy of daemon.
% on each trial daemon energy is compared with energy of all spin knots of grid
% and energy of knot and daemon is changed not changed according some
% condition
Es = zeros(1,NTrial*Nspin); % Pre-allocate space
Ed = Es;                    % Ditto
A = Es;                     % Ditto
SpM = Es;                   % Ditto
Es(1)=Esystem; 
Ed(1)=Edemon;
S=s;
k=1;

for i=1:NTrial
  Accept=0;
  for j=1:Nspin
     % random choice of knot
     Ix=floor(Ns*rand(1)+1); 
     Iy=floor(Ns*rand(1)+1);
     % border conditions
     if Ix==1
        Left=Ns;
     else
        Left=Ix-1;
     end
     if Ix==Ns
       Right=1;
     else
        Right=Ix+1;
     end
     if Iy==1
        Down=Ns;
     else  
        Down=Iy-1;
     end
     if Iy==Ns
        Up=1;
     else
        Up=Iy+1;
     end
     % energy change
     de=2*s(Iy,Ix)*(-h+J*(s(Iy,Left)+s(Iy,Right)+s(Down,Ix)+s(Up,Ix)));
     if de<=Edemon % energy change accepted
        s(Iy,Ix)=-s(Iy,Ix);
        Accept=Accept+1;
        Edemon=Edemon-de;
        Esystem=Esystem+de;
     end
     k=k+1;
     Es(k)=Esystem;
     Ed(k)=Edemon;
     A(k-1)=Accept;
     s1=sum(s);
     SpM(k)=sum(s1);
     S=cat(3,S,s);
  end
end
A=A/NTrial;
end

[Julia Dulie]

The code is from this book, starting from page 502 , the code is mentioned explicity in page 521

https://drive.google.com/file/d/1nka7w4_RlwxvCqkeVyczcivPjFbHGQHM/view?usp=sharing

By the way , would you -if that is possible- reccomand some books (or any tutrial materials, such as youtube channel) to learn simulation of such systems using Matlab ? the level is between begginner up to above intermediate ...

Thank you very much Alan Steven, you are very nice , I highly appreciate your remarks, and help

wish you have very nice day, and stay safe

[Alan Stevens]

Here are a couple of YouTube videos links:

https://www.youtube.com/watch?v=OgO1gpXSUzU

https://www.youtube.com/watch?v=Xaymy3Blnq4&t=20s

However, if you just search monte carlo simulation on YouTube you will find several.