% Program computes local and total DOS of a 1D nanowire attached to macroscopic reservoir(s)
% Definition: LDOS(E) = -1/pi * [Im G (E)]_ii | DOS(E) = -1/pi * Trace [Im G(E)]
% Definition: G(E)=[E-H-Sigma]^{-1} where Sigma is the self-energy of the reservoir(s)

clear all

%%
%Parameters and constants for modeling of clean wire
Np=100; %Number of atoms long the wire
t=1.0; %Nearest-neighbor hopping as the unit of energy
t_lead=t; %Hopping within the semi-infinite lead modeling macroscopic reservoir of fermions
t_couple=t; %Hopping between the wire and the macroscopic reservoir
site=Np/2; %LDOS will be computed in the middle of the wire

%%
%Hamiltonian matrix for clean wire
H(1:Np,1:Np)=0; %Set all elements of the Hamiltonian matrix of size Np x Np to zero
H=-(t*diag(ones(1,Np-1),1))-(t*diag(ones(1,Np-1),-1)); %This assumes hard wall boundary conditions

left_edge_potential=H(1,1); %store on-site potential on the left edge of the wire
right_edge_potential=H(Np,Np); %store on-site potential on the right edge of the wire
%%
%OPTIONAL: Add single impurity potential
%U=zeros(1,Np);
%U(round(Np/2))=3.0*t;
%H=H+diag(U);


%%
%Find LDOS from the diagonal elements of Im[G]
j=1 % count the number of energy points
for E=-5:0.01:5 %energy loop
    
    %Create self-energy of a single reservoir (modeled as semi-infinite lead) 
    rad=4*t_lead^2-E^2;
    if (rad>0.0)
        self_energy=(E-i*sqrt(rad));
    else
        self_energy=(E-sign(E)*sqrt(-rad));
    end
    self_energy=self_energy*t_couple^2/(2*t_lead^2); % final expression for the lead self-energy 
     
    H(1,1)=left_edge_potential + self_energy; %attach reservoir self-energy to the left edge of the wire
    %H(Np,Np)=right_edge_potential + self_energy; %attach reservoir self-energy to the right edge of the wire
    
    gr=inv(E*eye(Np)-H); % compute retarded Green function (self-energy added above knows about infinitesimal eta)
    im_g=-1/pi*((gr-gr')/(2*i)); %Im part of a matrix G is Im G = (G-G')/2*i
    ldos(j)=im_g(site,site); %LDOS at a chosen site
    dos(j)=-1/pi*trace((gr-gr')/(2*i)); %Total DOS=sum_over_sites LDOS
    j=j+1
    
end %energy loop

%%
%Plot LDOS and DOS versus Energy using the same grid as in the for loop above
E=-5:0.01:5; %Energy values at which dos will be plotted
subplot(2,1,1); % Subdivide the figure into 2 quadrants
plot(E,ldos);
grid;
title('Density of States for 1D nanowire composed of 100 atoms');
xlabel('Energy (t)');
ylabel('LDOS at site in the middle');
subplot(2,1,2); 
plot(E,dos);
grid;
xlabel('Energy (t)');
ylabel('Total DOS');

