% Simulated annealing approach to optimize the arrangement
% of particles on a linear chain where each site is coupled (with
% a random strength) to every other site. An "on-site" term is also 
% present for each lattice site.
%
% MCC Summer School 2005
% Gus Hart Northern Arizona Univesity
% mail: gus.hart "at" nau.edu

clc; close
clear

load Precomputed_EandH  % Load Hamiltonian and precomputed energies

N = length(H);          % Number of sites in the 1D lattice
BFmin = min(energies);  % Minimum found in brute-force search 
oldguess = ceil(rand(1,N)-.5);  % Random guess to start from
k = 0;          % Iteration counter
kT = 2;         % Starting emperature
dt = 1-.005;    % Temperature cooling factor (should be <1)
max_its = 2000; % Maximum number of iterations to run

energy = compute_energy(H,oldguess);  % Compute starting energy

while energy > BFmin  % Continue searching till we find the minimum
    rnd = rand(1);
    newguess = oldguess;  % Save old guess
    if rnd < .5 % try switching two sites
        site1=ceil(rand(1)*N);  % Select two sites 
        site2=ceil(rand(1)*N);  % to swap
        while oldguess(site1)==oldguess(site2)  % Keep selecting sites 
            site2=ceil(rand(1)*N);              % till we've selected both an 
        end                                     % empty and occupied site             
        temp = oldguess(site1);
        oldguess(site1) = oldguess(site2);
        oldguess(site2) = temp;
    else % Try adding or removing a particle
        site = ceil(rand(1)*N);
        oldguess(site) = mod(oldguess(site)+1,2);  % make 1 --> 0 or 0 --> 1
    end
    deltaE = compute_energy(H,oldguess) - energy;
    if deltaE < 0
        newguess = oldguess;
    elseif rand(1) < exp(-deltaE/kT)
        newguess = oldguess;
    else
        oldguess = newguess;
    end
    k = k + 1;
    if k > max_its
        break
    end
    energy = compute_energy(H,newguess);
    MCe(k) = energy;
    kT = kT*dt;
end
plot(MCe,'.-')
fprintf('Global min: %15.7f\n Found min: %15.7f\n',min(energies),energy)
if (min(energies)<energy)
    disp('Not converged')
else
    disp('Converged')
end

plot(MCe,'.-')