% Ising_lattice_MC.m
% This MATLAB program performs Monte Carlo simulation on an
% 2-D Ising lattice.
% Input:
% MCOpts - a data structure with members
%    .N = # of sites in each direction
%    .mu = site magnetization
%    .H = external magnetic field
%    .J = coupling constant
%    .kb_T = thermal energy
%    .S_init = initial site spins (as 2D array)
%    .Nequil = # equilibration MC steps
%    .Nsamples = # total number of samples
%    .freq_sample = # MC steps between samples
%    .make_plots = if non-zero, plot trajectories and final state
%
% Output:
%--------
% m_PROPS - data structure with magnetization information
%     .avg = average net magnetization
%     .var = variance of net magnetization
%     .order = 0 if no order, +1 if all + spins, -1 if all - spins
%     .traj = trajectory of net magnetization
% ener_PROPS - data structure with energy information
%     .avg = average energy
%     .var = variance of energy
%     .traj = trajectory of energy
% S = final state spins
%
% K. Beers. MIT ChE. 11/7/03
% v 11/14/03.A. corrected bug in magnetization definition
% v 12/23/03.A. added optional display of progress to screen
function [mag_PROPS,ener_PROPS,S,iflag_main] = Ising_lattice_MC(MCOpts);
iflag_main = 0;

% Extract simulation parameters from MCOpts and use defaults
% display
try
    display = MCOpts.display;
catch
    display = 0;
end
% Nequil
try
    Nequil = MCOpts.Nequil;
catch
    Nequil = 1e5;
end
if(display)
    Nequil,
end
% Nsamples
try
    Nsamples = MCOpts.Nsamples;
catch
    Nsamples = 1e5;
end
if(display)
    Nsamples,
end
% freq_sample
try
    freq_sample = MCOpts.freq_sample;
catch
    freq_sample = 10;
end
if(display)
    freq_sample,
end
% N
try
    N = MCOpts.N;
catch
    N = 100;
end
if(display)
    N,
end
% mu
try
    mu = MCOpts.mu;
catch
    mu = 1;
end
if(display)
    mu,
end
% H
try
    H = MCOpts.H;
catch
    H = 0;
end
if(display)
    H,
end
% J
try
    J = MCOpts.J;
catch
    J = 1;
end
if(display)
    J,
end
% kb_T
try
    kb_T = MCOpts.kb_T;
catch
    kb_T = 1;
end
if(display)
    kb_T,
end
% make_plots
try
    make_plots = MCOpts.make_plots;
catch
    make_plots = 0;
end
if(display)
    make_plots,
end

% Next, set initial spin states (use random dist. if no initial state input)
try 
    S = MCOpts.S_init;
catch
    S = rand(N,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
end

% Next, compute energy and magnetization of initial state
[ener,mag] = Ising_lattice_calc_ener(S,N,mu,H,J);


% ---------------------------
% Perform equilibration stage
if(display)
    disp(' '); disp('Beginning equilibration ...');
end
for iMC = 1:Nequil
    
    % every so often, recompute the total energy to avoid creeping errors
    if(~mod(iMC,N))
        [ener,mag] = Ising_lattice_calc_ener(S,N,mu,H,J);
    end
    if(display)
        if(~mod(iMC,10000))
            disp(['     Passing iMC = ', int2str(iMC)]);
        end
    end
    
    % generate new state by randomly changing one spin
    % and compute delta in system energy
    [i_change,j_change,delta_ener,delta_mag] = propose_MC_move(S,N,mu,H,J);
    
    % compute acceptance probability and see whether to accept move
    prob_ratio = exp(-delta_ener/kb_T);
    alpha = min(1, prob_ratio);
    rand_var = rand;  % random variable uniformly distributed on [0,1]
    if(rand_var <= alpha)  % accept move
        S(i_change,j_change) = -S(i_change,j_change);
        ener = ener + delta_ener;
        mag = mag + delta_mag;
    end
    
end


% --------------------
% Begin sampling stage
if(display)
    disp('Beginning sampling ...');
end

% allocate space for magnetization trajectory
mag_PROPS.traj = zeros(Nsamples,1);
ener_PROPS.traj = zeros(Nsamples,1);

for isample = 1:Nsamples
    
    % compute new energy to avoid creeping error
    % in energy and magnetization
    [ener,mag] = Ising_lattice_calc_ener(S,N,mu,H,J);

    if(display)
        if(~mod(isample,1000))
            disp(['     Passing isample = ', int2str(isample)]);
        end
    end

    
    
    % perform some number of MC steps between samples
    for iMC = 1:freq_sample
        % generate new state by randomly changing one spin
        % and compute delta in system energy
        [i_change,j_change,delta_ener,delta_mag] = ...
            propose_MC_move(S,N,mu,H,J);
        
        % compute acceptance probability and see whether to accept move
        prob_ratio = exp(-delta_ener/kb_T);
        alpha = min(1, prob_ratio);
        rand_var = rand;  % random variable uniformly distributed on [0,1]
        if(rand_var <= alpha)  % accept move
            S(i_change,j_change) = -S(i_change,j_change);
            ener = ener + delta_ener;
            mag = mag + delta_mag;
        end
    end
    % record current result in trajectory
    ener_PROPS.traj(isample) = ener;
    mag_PROPS.traj(isample) = mag;
    
end

% ----------------------------------------------------------
% compute statistical properties of magnetization and energy
mag_PROPS.avg = mean(mag_PROPS.traj);
mag_PROPS.var = var(mag_PROPS.traj);
mag_PROPS.order = mag_PROPS.avg / (N^2) / mu;
ener_PROPS.avg = mean(ener_PROPS.traj);
ener_PROPS.var = var(ener_PROPS.traj);

% ---------------------------------------------------------
% if desired, make plots of final state and the trajectory
% information of the magnetization and energy
if(make_plots)    
    
    % final state space
    figure;
    spy(S + ones(size(S)));
    xlabel('lattice x-site #i');  ylabel('lattice y-site #j');
    title_phrase = ['Ising lattice + spins: ', ...
            'N = ', int2str(N), ...
            ', \mu = ', num2str(mu), ...
            ', H = ', num2str(H), ...
            ', J = ', num2str(J), ...
            ', k_bT = ', num2str(kb_T)];
    title(title_phrase);
   
    % magnetization properties
    figure;
    subplot(2,1,1);
    plot(mag_PROPS.traj);
    title_phrase = ['Ising lattice magnetization : ', ...
            'N = ', int2str(N), ...
            ', \mu = ', num2str(mu), ...
            ', H = ', num2str(H), ...
            ', J = ', num2str(J), ...
            ', k_bT = ', num2str(kb_T)];
    xlabel('MC sample');  ylabel('m');
    title(title_phrase);
    subplot(2,1,2);
    hist(mag_PROPS.traj/Nsamples,20);
    xlabel('m');  ylabel('Pr(m)');
    title(['avg. = ', num2str(mag_PROPS.avg),', variance = ', num2str(mag_PROPS.var), ...
            ', order param = ', num2str(mag_PROPS.order)]);
    
    % energy properties
    figure;
    subplot(2,1,1);
    plot(ener_PROPS.traj);
    title_phrase = ['Ising lattice energy : ', ...
            'N = ', int2str(N), ...
            ', \mu = ', num2str(mu), ...
            ', H = ', num2str(H), ...
            ', J = ', num2str(J), ...
            ', k_bT = ', num2str(kb_T)];
    xlabel('MC sample');  ylabel('energy');
    title(title_phrase);
    subplot(2,1,2);
    hist(ener_PROPS.traj/Nsamples,20);
    xlabel('E');  ylabel('Pr(E)');
    title(['avg. = ', num2str(ener_PROPS.avg),', variance = ', num2str(ener_PROPS.var)]);

end   


iflag_main = 1;
return;


% ===============================================
% This routine computes the total system energy
% and net magnetization of the 2-D Ising lattice
function [ener,mag] = Ising_lattice_calc_ener(S,N,mu,H,J);

% First, compute energy due to coupling with external
% field
mag = mu*sum(sum(S));  % net magnetization
ener_ext_H = H*mag;

% Next, compute energy associated with coupling
% between each state
sum1 = 0;
for i=1:N
    for j=1:N
        % (i-1,j)
        m = i-1;
        if(m==0)  % enforce periodic BC
            m = N;
        end
        n = j;
        sum1 = sum1 + S(i,j)*S(m,n);
        % (i+1,j)
        m = i+1;
        if(m==(N+1))
            m = 1;
        end
        n = j;
        sum1 = sum1 + S(i,j)*S(m,n);
        % (i,j-1)
        m = i;
        n = j-1;
        if(n==0)
            n = N;
        end
        sum1 = sum1 + S(i,j)*S(m,n);
        % (i,j+1)
        m = i;
        n = j+1;
        if(n==(N+1))
            n = 1;
        end
        sum1 = sum1 + S(i,j)*S(m,n);
    end
end
ener_J_couple = -J/2*sum1;

% get total system energy
ener = ener_ext_H + ener_J_couple;

return;


% ==================================
% This routine selects at random a spin to change and
% computes the energy of doing so.
function [i_change,j_change,delta_ener,delta_mag] = ...
    propose_MC_move(S,N,mu,H,J);

i_change = round(  1 + rand*(N-1));
j_change = round(  1 + rand*(N-1));

% compute current and new energies of this site interacting
% with its neighbors
sum1 = 0;
sum1_new = 0;

% (i_change-1,j_change)
m = i_change-1;
if(m==0)  % enforce periodic BC
    m = N;
end
n = j_change;
sum1 = sum1 + S(i_change,j_change)*S(m,n);
sum1_new = sum1_new - S(i_change,j_change)*S(m,n);
% (i_change+1,j_change)
m = i_change+1;
if(m==(N+1))
    m = 1;
end
n = j_change;
sum1 = sum1 + S(i_change,j_change)*S(m,n);
sum1_new = sum1_new - S(i_change,j_change)*S(m,n);
% (i_change,j_change-1)
m = i_change;
n = j_change-1;
if(n==0)
    n = N;
end
sum1 = sum1 + S(i_change,j_change)*S(m,n);
sum1_new = sum1_new - S(i_change,j_change)*S(m,n);
% (i_change,j_change+1)
m = i_change;
n = j_change+1;
if(n==(N+1))
    n = 1;
end
sum1 = sum1 + S(i_change,j_change)*S(m,n);
sum1_new = sum1_new - S(i_change,j_change)*S(m,n);

ener = -J*sum1;
ener_new = -J*sum1_new;
delta_ener = ener_new - ener;

% Than, add change in energy in external field and get
% total magnetization
delta_mag = -2*mu*S(i_change,j_change);
delta_ener = delta_ener + H*delta_mag;


return;