%Fig2_8b.m
%Essential Electron Transport for Device Physics
%tight-binding model in 1D
%random search for objective function
clear all; clc;
clf; %close all;

%plotting parameters
FS=12;  %label fontsize
LW=1;   %curve linewidth
MS=7.5; %marker size

% Change default axes fonts.
set(0,'DefaultAxesFontName', 'Times');
set(0,'DefaultAxesFontSize', FS);

% Change default text fonts.
set(0,'DefaultTextFontname', 'Times');
set(0,'DefaultTextFontSize', FS);

%parameters
A=3;            %total number of atoms on lattice 3 6 9
t=-1;           %hopping energy
alpha=3;        %long-range decay exponent 3
Gamma=0.3*abs(t);    %characteristic energy broadening

L=12;           %length of one-dimensional domain 4*A 10 20
%setup 
emax=2*abs(t);         %minimum energy in 1D band
emin=-2*abs(t);        %maximum energy in 1D band
Eobj=linspace(emin,emax,A)%evenly spaced eigenvalues as objective function
%Eobj=-2*t*cos(linspace(0,pi,A))%nearest neighbor tight-binding 
% note: rely on fact eigenenergies are in assending order   
%**************************************************************************
maxcount=100000;         %maximum number of iterations
x=L*rand(A,1);          %initialize atom positions x
H=zeros(A,A);           %initialize Hamiltonian H
J=zeros(maxcount,1);    %iniitialize cost function J
Jav=zeros(maxcount,1);  %initialize average cost function Jav
%**************************************************************************
%*******************         main loop         ****************************
%**************************************************************************
kavmax=100;             %average over 100 runs
for kav=1:kavmax
    for k=1:maxcount    %maxcount is bound on maximum number of iterations    
%**************************************************************************
% call EvalEig to find simulated eigenvalues, Esim
%**************************************************************************
Esim = EvalEig(x,A,t,alpha,L);
%*************************************************************************
%cost function is sum of differences in ordered eigenenergies squared
J(k) = sum((Esim'-Eobj).^2);%distance is abs(Esim'-Eobj) or (Esim'-Eobj).^2
%*************************************************************************
%*************************************************************************
    x=L*rand(A,1);
    end
%prepare cost function for convergence plot
for i=2:k
    if J(i)>J(i-1);
        J(i)=J(i-1);
    end
end
Jav=(Jav+J);      
end
ttl=['\rm Fig2.8b L=',num2str(L),', A=',num2str(A),', t=',num2str(t),...
     ', \alpha=',num2str(alpha),', J=',num2str(J(k)),...
     ', n_{count}=',num2str(k),', n_{av}=',num2str(kavmax)]; 
%*************************************************************************
%plot figure Jav from 10th interation to maxcount iteration
%*************************************************************************
figure(1)
loglog(Jav/kavmax);title(ttl)
grid on;hold on;
axis([1 maxcount 10^-7 10^4]);%
%axis([1 maxcount 0 Jav(10)]);%
xlabel('Iteration','Fontsize',FS);ylabel('Convergence','Fontsize',FS);
%*************************************************************************
% function EvalEig to find simulated eigenvalues, Esim
%*************************************************************************
function Esim = EvalEig(x,A,t,alpha,L)
Hltrig=zeros(A,A);%initialize lower-triangle and ensure H diagonal is zero
    for i=2:A
        for j=1:i-1
dx=abs(x(i)-x(j));% difference in position dx < L
Hltrig(i,j)=-t/(dx^alpha)-t/((L-dx)^alpha);%periodic boundary conditions in 1D
        end
    end   
%Hutrig=Hltrig', note upper-triangle is transpose of lower-triangle
H=Hltrig+Hltrig';%H with diagonals zero
Esim = eig(H);   %eigenergies note: rely on fact eigenenergies are in order   
end
