%Fig2_8a.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 6 9
t=-1;       %hopping energy
alpha=3;    %long-range exponent 3
Gamma=0.3*abs(t);%characteristic energy broadening
minE=-4*abs(t);  %minimum energy 3
maxE=4*abs(t);   %maximum energy -3
nE=1000;    %size of energy array in density of states plot
E=linspace(minE,maxE,nE);%energy range for density of states plot

L=12;       %length of one-dimensional domain 12
%setup 
emax= 2;    %minimum energy in 1D band
emin=-2;    %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   
NEobj=zeros(1,nE);  %initialize objective density of states array
%**************************************************************************
maxcount=100000;   %maximum number of iterations
Jth=10^-6;  %acceptance threshold for cost function 0.1 0.03 0.01 0.001
H=zeros(A,A);       %initialize H
J=zeros(maxcount,1);%initialize J
Gamma2=(Gamma/2).^2;
Gamma2pi=2*Gamma/(2*pi);%note factor 2 for spin
%*************************************************************************
%calculate objective density of states over energy range of interest
for i=1:A
    NEobj=NEobj+(Gamma2pi./((E-Eobj(i)).^2+Gamma2));%note factor 2 for spin   
end
NEobj=NEobj*(maxE-minE)/(nE);%density of states per unit energy interval
%initial atom positions
x=L*rand(A,1);

%**************************************************************************
%*******************         main loop         ****************************
%**************************************************************************
for k=1:maxcount     %maxcount is bound on maximum number of iterations   
    for i=1:A
        for j=1:A
            if i==j
H(i,j)=0;%do not count i=j contributions, set diagonal to zero
            else
    dx=abs(x(i)-x(j));% difference in position dx < L
H(i,j)=-t/(dx^alpha)-t/((L-dx)^alpha);%periodic boundary conditions in 1D
%H(i,j)=-t/(dx^alpha);%non-periodic pairs
            end
        end
    end   
Esim= eig(H);% note: rely on fact eigen energies are in order   
%*************************************************************************
%cost function is sum of differences in ordered eigenenergies squared
    J(k)=sum((Esim'-Eobj).^2);
%*************************************************************************
%if J<Jth plot NE and quit
    if J(k)<Jth
%calculate density of states over energy range of interest
NE=zeros(1,nE);%initialize density of states array
for i=1:A
    NE=NE+(Gamma2pi./((E-Esim(i)).^2+Gamma2));%note factor 2 for spin    
end
NE=NE*(maxE-minE)/(nE);%density of states per unit energy interval
%int=sum(NE)%check integral
%prepare cost function for convergence plot
for i=2:k
    if J(i)>J(i-1);
        J(i)=J(i-1);
    end
end

ttl1=['\rm Fig2.8a L=',num2str(L),', A=',num2str(A),...
     ', t=',num2str(t),', \alpha=',num2str(alpha),', J=',num2str(J(k)),...
     ', J_{th}=',num2str(Jth),', n_{count}=',num2str(k)]; 
ttl2=['\rm Fig2.8a obj(red), L=',num2str(L),', A=',num2str(A),...
     ', t=',num2str(t),', \Gamma=',num2str(Gamma),...
     ', \alpha=',num2str(alpha),', J=',num2str(J(k)),...
     ', J_{th}=',num2str(Jth),', n_{count}=',num2str(k)]; 
%*************************************************************************
%plot figures
figure(1)
loglog(J);title(ttl1)
grid on;hold on
axis([1 maxcount 10^-7 10^4]);%
%axis([1 maxcount 0 10000*Jth]);
xlabel('Iteration','Fontsize',FS);ylabel('Convergence','Fontsize',FS);

figure(2)
plot(E,NEobj,'r');hold on;plot(E,NE);grid on;title(ttl2)
xlabel('Energy, E/t','Fontsize',FS);ylabel('Density of states, N(E)','Fontsize',FS);
hold off

x=sort(x) %sorted values of atom positions x
y=zeros(length(x),1);

figure(3)
plot(x,y,'or');grid on
axis([0 L -1 1]);
xlabel('Position, x','Fontsize',FS);title(ttl1);
  return
%*************************************************************************
%if k<maxcount then try new random values of position x    
    elseif k<maxcount
    x=L*rand(A,1);
    end
end


