%Fig2_7.m
%Essential Electron Transport for Device Physics
%tight-binding model in 2D
%guided random walk towards objective function
clear all;
clf;

%plotting parameters
FS=12;  %label fontsize
LW=1;   %curve linewidth
MS=7.5; %marker size
npo=20; %number of printouts for iteration loop

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

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

tic;%initiate tic-toc timer

%switches
plotting=1; %option to generate figures during optimization
contplot=1; %option to generate contour plot at end
x_initial=1; %'1'=use random atom positions, '2'=use estimated atom positions from figure
            %using x_initial=2 requires that A1=12,A2=4,Lx=Ly=L=80,dx=0.1
x_obj=2; %'1'=use randomly generated atomic positions as objective;
         %'2'=use custom atomic configuration
         %'3'=use test configuration (to determine interatomic distance cutoff)
randtype=3; %type of randomization:
            %'1'=completely new random arrangement of all atoms
            %'2'=random walk, one atom selected at random
            %'3'=random walk, all atoms
maxcount=1e5;   %maximum number of iterations 1e4
Jth=0.0002;   %acceptance threshold for cost function 0.005

%miscellaneous options
periodic=1; %implement periodic boundary conditions ('1'=yes,'0'=no)
allints=0; %option to include sum of interactions between given atom and 
           %all other atoms in central cell and in every surrounding cell.
options=[periodic,allints];
levels=60; %number of contour lines (levels) for contour plot 20
conthresh=0.50; %threshold for contour plot 0.2

%parameters
A=8;        %total number of atoms 8
t=1;        %magnitude of hopping energy
alpha=2;    %long-range decay exponent 3 2 1
Gamma=0.2828*t;%characteristic energy broadening
EX=4;       %hopping strength multiplier (energy limits) 10
minE=-EX*t; %minimum energy 3
maxE= EX*t; %maximum energy -3
nE=1000;    %size of energy array in density of states plot
dminmaxE=(maxE-minE)/(nE); %prefactor density of states per unit energy interval
E=linspace(minE,maxE,nE);%energy range for density of states plot

L=40;       %length of two-dimensional square domain 40
Lx=L;       %length of 2D domain (along x-direction)
Ly=L;       %length of 2D domain (along y-direction)
dx=0.1;     %size of spatial increment (between lattice sites) 0.1

if x_initial==2 && isequal([A1,A2,Lx,Ly,dx],[12,4,80,80,0.1])==0
    error(['If x_initial==2, the following must be true: A1=12,A2=4,Lx=Ly=80,dx=0.1']);
else
    %store relevant parameters for parallel loop
    params=[randtype,A,Lx,Ly,t,alpha,dx];
end

%objective function setup
config_obj=sparse(Ly,Lx);
if x_obj==1
    %random starting configuration
    a=0;
    while a<A
        indx=randi(Lx,1);
        indy=randi(Ly,1);
        if config_obj(indy,indx)==0
            config_obj(indy,indx)=1;
            a=a+1;
        end
    end
elseif x_obj==2
    

    config_obj(5,5)=1;
    config_obj(5,15)=1;
    config_obj(5,25)=1;
    config_obj(5,35)=1;
    config_obj(20,5)=1;
    config_obj(20,15)=1;
    config_obj(20,25)=1;
    config_obj(20,35)=1;

%     config_obj(5,8)=1;
%     config_obj(5,16)=1;
%     config_obj(5,24)=1;
%     config_obj(5,32)=1;
%     config_obj(20,15)=1;
%     config_obj(20,23)=1;
%     config_obj(20,31)=1;
%     config_obj(20,39)=1;

elseif x_obj==3
    %test configuration    
    config_obj(1,1)=1;
    config_obj(1,11)=1;    
%     config_obj(1,Lx-1)=1;
%     config_obj(1,Lx)=1;
end
% full(config_obj)
Eobj = EvalEig2D(config_obj,params,options); %objective eigenenergies
%**************************************************************************
% maxcount=2;   %maximum number of iterations
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
NEobj=sum(Gamma2pi./((E-Eobj).^2+Gamma2));
NEobj=NEobj*dminmaxE;   %density of states per unit energy interval
    
Ecsw=1; %energy constraint switch
while Ecsw==1
%initial atom positions
config=sparse(Ly,Lx); %configuration
if x_initial==1 %random configuration
    a=0;
    while a<A
        indx=randi(Lx,1);
        indy=randi(Ly,1);
        if config(indy,indx)==0
            config(indy,indx)=1;
            a=a+1;
        end
    end
else %configuration estimated from figure    
    if A==16
        %inferred positions of atoms in Fig.1.2, pg. 4 of Optimal Device
        %Design (A.F.J. Levi and S. Haas), coordinate pair ordering: (y,x)
        
        %guessing atom positions (row,column), assuming 80x80 lattice
        ap=[8,18;... 
            6,36;...
            7,57;... 
            80,60;...
            15,13;...            
            16,22;...
            46,8;... 
            54,4;...
            51,36;...
            55,43;...
            43,47;...
            38,57;...
            45,72;...
            52,76;...
            64,32;...
            78,49]; 
        
        for i=1:16
            config(round(ap(i,1)),round(ap(i,2)))=1; 
        end
    end    
end
Esim = EvalEig2D(config,params,options);
if min(Esim)>=minE && max(Esim)<=maxE
    Ecsw=0;
end
end

J(1)=sum((Esim-Eobj).^2); %compute error
% Esim = -Esim; %for whatever reason, the hopping energy seems to be
                %flipped in textbook figure
NE=sum(Gamma2pi./((E-Esim).^2+Gamma2));%note factor 2 for spin
NE=NE*dminmaxE;   %density of states per unit energy interval
%**************************************************************************
%*******************         main loop         ****************************
%**************************************************************************
for k=1:maxcount %maxcount is bound on maximum number of iterations

if mod(k,round(maxcount/npo))==0
    percentdone=100*k/maxcount;
    disp([' Total iterations are ',num2str(percentdone),'% complete']);
end  
    
if J(k)<Jth %below threshold, minimum error achieved
    
    ncount=k;        %number of iterations
    tcomplete=toc;   %time to completion of calculation and plot
    break;
    
else %need to change atom positions    

Ecsw=1; %energy constraint switch
while Ecsw==1
    if randtype==1 %completely new random distribution

        ctemp=sparse(Ly,Lx); %configuration
        a=0;
        while a<A
            indx=randi(Lx,1);
            indy=randi(Ly,1);
            if ctemp(indy,indx)==0
                ctemp(indy,indx)=1;
                a=a+1;
            end
        end        
        
    elseif randtype==2 %random walk (1 random atom at a time)
                
        %randomly decide which direction for each atom to independently take
        ctemp=config; %create a temporary copy of the guess configuration to edit
        alocs=find(ctemp==1); %raw matrix coordinates of objective atoms
        x=ceil(alocs/Lx); %x-coordinates of objective atoms (columns)
        y=mod(alocs,Ly); %y-coordinates of objective atoms (rows)
        y(y==0)=Ly; %correct for atoms found at bottom (top) edge
        a2m=randi(A,1); %randomly pick an atom to move
        indx=x(a2m); %determine x-coordinate of randomly-selected atom
        indy=y(a2m); %determine y-coordinate of randomly-selected atom
        rs=0; %switch to accept random hop to adjacent site
        while rs==0
            indx_new=indx+randi(3,1)-2; %randomly move atom either left or right
            if indx_new==0
                indx_new=Lx; %send leftmost atom to right end (periodic BCs)
            elseif indx_new==Lx+1
                indx_new=1; %send rightmost atom to left end (periodicity)
            end
            indy_new=indy+randi(3,1)-2; %randomly move atom either up or down
            if indy_new==0
                indy_new=Ly; %send downmost atom to up end (periodic BCs)
            elseif indy_new==Ly+1
                indy_new=1; %send upmost atom to down end (periodicity)
            end
            ctemp(indy,indx)=0; %remove randomly selected atom from old site
            ctemp(indy_new,indx_new)=1; %add randomly selected atom to new site
            if sum(ctemp(:))==A %check to ensure a maximum of 1 atom per site            
                rs=1; %accept new position            
            end
        end
                
    elseif randtype==3 %random walk (all atoms at a time)
        
        %randomly decide which direction for each atom to independently take
        ctemp=config; %create a temporary copy of the guess configuration to edit
        alocs=find(ctemp==1); %raw matrix coordinates of objective atoms
        x=ceil(alocs/Lx); %x-coordinates of objective atoms (columns)
        y=mod(alocs,Ly); %y-coordinates of objective atoms (rows)
        y(y==0)=Ly; %correct for atoms found at bottom (top) edge
%         a2m=randi(A,1); %randomly pick an atom to move
        for a=1:A
            indx=x(a); %determine x-coordinate of randomly-selected atom
            indy=y(a); %determine y-coordinate of randomly-selected atom
            rs=0; %switch to accept random hop to adjacent site
            while rs==0
                indx_new=indx+randi(3,1)-2; %randomly move atom either left or right
                if indx_new==0
                    indx_new=Lx; %send leftmost atom to right end (periodic BCs)
                elseif indx_new==Lx+1
                    indx_new=1; %send rightmost atom to left end (periodicity)
                end
                indy_new=indy+randi(3,1)-2; %randomly move atom either up or down
                if indy_new==0
                    indy_new=Ly; %send downmost atom to up end (periodic BCs)
                elseif indy_new==Ly+1
                    indy_new=1; %send upmost atom to down end (periodicity)
                end
                ctemp(indy,indx)=0; %remove randomly selected atom from old site
                ctemp(indy_new,indx_new)=1; %add randomly selected atom to new site
                if sum(ctemp(:))==A %check to ensure a maximum of 1 atom per site            
                    rs=1; %accept new position            
                end
            end
        end        
    end
    
%**************************************************************************
% call EvalEig to find simulated eigenvalues, Esim
%**************************************************************************
    Esim = EvalEig2D(ctemp,params,options);
    if min(Esim)>=minE && max(Esim)<=maxE
        Ecsw=0;
    end
end    
%*************************************************************************
%cost function is sum of differences in ordered eigenenergies squared
    J(k+1)=sum((Esim-Eobj).^2); %compute cost for new configuration
%*************************************************************************    
    %decide whether to keep new configuration
    if J(k+1)<J(k)
        config=ctemp;
        
        NE=sum(Gamma2pi./((E-Esim).^2+Gamma2)); %note factor 2 for spin
        NE=NE*dminmaxE;   %density of states per unit energy interval
        
        ttl1={['L_{x}\timesL_{y}=',num2str(Lx),...
            '\times',num2str(Ly),', A=',num2str(A),', t=',num2str(t),...
            ', \Gamma=',num2str(Gamma),', \alpha=',num2str(alpha)];...
            ['J=',num2str(J(k+1)),', J_{th}=',num2str(Jth),...
            ', n_{count}=',num2str(k)]};
        ttl2={['obj(red), sim(blue), L_{x}\timesL_{y}=',num2str(Lx),...
            '\times',num2str(Ly),', A=',num2str(A),', t=',num2str(t),...
            ', \Gamma=',num2str(Gamma),', \alpha=',num2str(alpha)];...
            ['J=',num2str(J(k+1)),', J_{th}=',num2str(Jth),...
            ', n_{count}=',num2str(k)]};
               
        if plotting==1
            
        % plot figures
        figure(1); hold off;
        loglog(J);
        % axis([1,maxcount,0,10000*Jth]);
        axis([1,maxcount,0,max(J(2:end))]);
        xlabel('Iteration','Fontsize',FS);
        ylabel('Convergence','Fontsize',FS);
        title(ttl1); drawnow; hold off;
        
        figure(2);
        plot(E,NEobj,'r'); hold on; plot(E,NE);
        axis([minE,maxE,0,1.2*max(NEobj)]);
        title(ttl2)
        xlabel('Energy, E/t','Fontsize',FS);
        ylabel('Density of states, N(E)','Fontsize',FS);
        title(ttl2); drawnow; hold off;
        
        %plot atom positions on 2D lattice
        figure(3);
        spy(config_obj,'or',MS); hold on;
        spy(config,'*b',MS);
        ax=gca; ax.YDir='normal';
        axis([0,Lx,0,Ly]);
        xlabel('Position, x','Fontsize',FS);
        ylabel('Position, y','Fontsize',FS);
        title(ttl2); drawnow; hold off;
        
        end
        
    else
        J(k+1)=J(k);            
    end
    
end
end

disp(' ');
toc
disp(['Cost of final/best configuration: J_final = ',num2str(J(k))]);

clf;
% close all;

% plot figures
fig1=figure(1); hold on; grid on;
loglog(J,'-k','Linewidth',LW);
axis([1,maxcount,0,max(J)]);
xlabel('Iteration','Fontsize',FS);ylabel('Convergence','Fontsize',FS);
title(ttl1); drawnow; hold off;

fig2=figure(2); hold on; grid on;
plot(E,NEobj,'r','Linewidth',LW); 
plot(E,NE,'-b','Linewidth',LW);
axis([minE,maxE,0,1.2*max(NEobj)]);
xlabel('Energy, E/t_{hop}','Fontsize',FS);
ylabel('Density of states, N(E)','Fontsize',FS);
title(ttl2); drawnow; hold off;

%plot atom positions on 2D lattice
fig3=figure(3); hold on; grid on;
spy(config_obj,'or',MS); 
spy(config,'*b',MS);
ax=gca; ax.YDir='normal';
axis([0,Lx,0,Ly]);
xlabel('Position, x','Fontsize',FS);
ylabel('Position, y','Fontsize',FS);
title(ttl2); drawnow; hold off;
      
if contplot==1
%     clf;
    %extract positions (to compute distances) from configuration matrix
    alocs=find(config==1); %raw matrix coordinates of objective atoms
    x=ceil(alocs/Lx); %x-coordinates of objective atoms (columns)
    y=mod(alocs,Ly); %y-coordinates of objective atoms (rows)
    y(y==0)=Ly; %correct for atoms found at bottom (top) edge
    enland=sparse(Ly,Lx); %energy landscape in unit cell
    for j=1:Ly
        for i=1:Lx            
            stemp=0; %initialize sum of energy contributions from atoms in lattice
            for k=1:A                
                if periodic==0
                    deltax = i-x(k);
                    deltay = j-y(k);
                elseif periodic==1
                    deltax=[i-(x(k)-Lx),i-x(k),i-(x(k)+Lx),...
                            i-(x(k)-Lx),i-x(k),i-(x(k)+Lx),...
                            i-(x(k)-Lx),i-x(k),i-(x(k)+Lx)];
                    deltay=[j-(y(k)+Ly),j-(y(k)+Ly),j-(y(k)+Ly),...
                            j-y(k),j-y(k),j-y(k),...
                            j-(y(k)-Ly),j-(y(k)-Ly),j-(y(k)-Ly)];
                end
                dist=sqrt(deltax.^2+deltay.^2)*dx;%/sqrt(2); %compute array of distances
                temp=sum(-t./(dist.^alpha)); %compute potential energy contribution fro k-th atom
                stemp=stemp+temp;
            end
            enland(j,i)=-stemp;
        end
    end
    
    ind1=find(enland<Inf); %matrix indices of unoccupied sites
    ind2=find(enland==Inf); %matrix indices of occupied sites
    temp1=enland(ind1); %unoccupied subset
    max1=max(temp1(:)); %max value
    enland(ind2)=max1; %set values of all occupied sites to max surrounding value
    %contour line ranges:
    min_level=min(enland(:));
    rnge=max(enland(:))-min(enland(:)); %"range" function requires statistics toolbox
    max_level=min_level+conthresh*rnge;
    vec_level=full(linspace(min_level,max_level,levels));

    fig4=figure(4);
    contour(enland,vec_level);
    axis image; %axis square
    axis([0.5,Lx+0.5,0.5,Ly+0.5]);
    xlabel('Position, x','Fontsize',FS);
    ylabel('Position, y','Fontsize',FS);
    ax=gca; ax.YDir='normal'; title(ttl1);
    hold off;
end



%*************************************************************************
% function EvalEig2D to find simulated eigenvalues for 2D case, Esim
%*************************************************************************
function Esim = EvalEig2D(config,params,options)
%extract parameters from array
% randtype = params(1);
A        = params(2);
Lx       = params(3);
Ly       = params(4);
t        = params(5);
alpha    = params(6);
dx       = params(7);
%options
periodic = options(1); %whether to use 2D periodic boundary conditions
allints  = options(2); %use either NN or all neighboring cell interactions

H=zeros(A,A); %Hamiltonian matrix for random guess
alocs=find(config==1); %raw matrix coordinates of objective atoms
x=ceil(alocs/Lx); %x-coordinates of objective atoms (columns)
y=mod(alocs,Ly); %y-coordinates of objective atoms (rows)
y(y==0)=Ly; %correct for atoms found at bottom (top) edge

for i=1:A %loop through all atoms
    ja=1:A; ja(i)=[]; %other atom index array
    for j=1:length(ja) %for given ith atom, loop through all OTHER atoms (total: A-1)
        %implement periodic boundary conditions to find minimum distances
        %in 9 cells (8 surrounding cells + central cell)
        if periodic==0
            deltax=x(i)-x(ja(j));
            deltay=y(i)-y(ja(j));
        elseif periodic==1
            deltax=[x(i)-(x(ja(j))-Lx),x(i)-x(ja(j)),x(i)-(x(ja(j))+Lx),...
                    x(i)-(x(ja(j))-Lx),x(i)-x(ja(j)),x(i)-(x(ja(j))+Lx),...
                    x(i)-(x(ja(j))-Lx),x(i)-x(ja(j)),x(i)-(x(ja(j))+Lx)];
            deltay=[y(i)-(y(ja(j))+Ly),y(i)-(y(ja(j))+Ly),y(i)-(y(ja(j))+Ly),...
                    y(i)-y(ja(j)),y(i)-y(ja(j)),y(i)-y(ja(j)),...
                    y(i)-(y(ja(j))-Ly),y(i)-(y(ja(j))-Ly),y(i)-(y(ja(j))-Ly)];
        end
        dist=sqrt(deltax.^2+deltay.^2)*dx;%/sqrt(2); %compute array of distances
    %         sort_dist=sort(dist); dist=sort_dist(1:2); %select 4 shortest distances
        if allints==0 %only consider nearest-neighbors
            H(i,ja(j))=-t/(min(dist)^3);
        elseif allints==1 %include all nearest cell contributions
            H(i,ja(j))=sum(-t./(dist.^alpha));
        end
    end
end
Esim=eig(H); %extract eigenvalues
if numel(Esim)~=A
    error('Error! Number of non-zero eigenvalues does not equal A');
end
end

