% Program computes total DOS of a finite length 1D clean nanowire via Green functions
% Definition: DOS(E) = -1/pi * Trace [Im G(E)]
% Definition: G(E)=[E+eta-H]^{-1} where eta is small imaginary number required to defind the retarded GF

clear all

%%
%Parameters and constants for clean wire
Np=100; %Number of atoms long the wire
t=1.0; %Nearest-neighbor hopping as the unit of energy
eta=i*10^(-3); %Small imaginary number needed to define retarded G^r and advanced G^a Green functions; 
% if eta is too small (try it), it is impossible numerically do differentiate between G^r=inv((E+eta)*eye(Np)-H) and G^a=inv((E-eta)*eye(Np)-H)
% remember not to use "i" as variable other than i^2=-1 when working with complex numbers

%%
%Hamiltonian matrix for clean wire
H(Np,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));
H(1,Np)=-t; H(Np,1)=-t; %Use periodic boundary conditions

%%
%OPTIONAL: Add single impurity potential
%U=zeros(1,Np);
%U(round(Np/2))=3.0*t;
%H=H+diag(U);


%%
%Find DOS from the trace of imaginary G
j=1
for E=-5:0.01:5 %energy loop
    gr=inv((E+eta)*eye(Np)-H);
    dos(j)=-1/pi*trace((gr-gr')/(2*i)); %Im part of a matrix G is Im G = (G-G')/2i
    j=j+1
end %energy loop

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