% The program computes conductance and local density of states for a
% graphene nanoribbon (GNR) attached to two clean semi-infinite GNR electrodes

% Self-energies due to GNR electrodes are computed numerically using
% Lopez-Sancho accelerated recursive method

% The retarded Green functions is computed recursively after Hamilotnian 
% reordering into block-tridiagonal form

% Authors: Farzad Mahfouzi and Branislav K. Nikolic

%% Variables:

clear all;
close all;

Nz=3;  % length of the ribon in terms of number of cells
N_y=4; % width of the ribon
t=1.0; % nearest neighbor hopping parameters in the central region
t_l=t; % hooping inside the leads
t_couple=t; % hopping between the lead and the sample
on_site_potential=0.0;

Ne=200; % number of energy points
emax=3.5; % upper limit of the energy interval
emin=-emax; % lower limit of the energy interval

%% Create Supercell Hamiltonian
N_y=2*N_y; % number of carbon atoms within the supercell
t_nnn=0.1*t; % next nearest neighbor hopping on honeycomb lattice (fitting DFT suggests 0.1*t)
[Tl1,h_l1,Tr1]=h_zigzag(on_site_potential,t_l,t_nnn,N_y);
[T1,h0,T2]=h_zigzag(on_site_potential,t,t_nnn,N_y);

%% Check the subband structure defined by the supercell Hamiltonian
figure(1),[emin,emax]=bstruct(T1,h0,T2,pi,100);
emin=emin-0.2*abs(emin);
emax=emax+0.2*abs(emax); 

%% Create the diagonal and upper/lower dagonal parts of H
H_c=zeros(size(h0,1),size(h0,1),Nz);
T_L=zeros(size(H_c));
T_R=zeros(size(H_c));
S=zeros(size(H_c));
onemat=zeros(size(H_c));
for m=1:size(H_c,3)
   onemat(:,:,m)=eye(size(H_c,1));
   H_c(:,:,m)=h0;
   T_L(:,:,m)=T1;
   T_R(:,:,m)=T2;
end

%% Create our sparse Hamiltonian using the following function
Hamil=blocktosparse(T_L,H_c,T_R,0);

%% Here we define the points that are going to be connected to leads
Nsc=size(H_c,1);

i=1:Nsc;
x(i)=i; % points that first lead will be connected

i1=Nsc*(size(H_c,3)-1)+i;
x(Nsc+i)=i1; % points that second lead will be connected

%% Now we get the following block diagonal matrices
[Tb1,Hb,Tb2,s,nfs]=sparsetoblock(Hamil,x);


%% Block-tridiagonal form for matrices 
S=zeros(size(H_c));
onemat=zeros(size(H_c));
for m=1:size(H_c,3)
   onemat(:,:,m)=eye(size(H_c,1));
end
Sb=zeros(size(Hb));
onematb=zeros(size(Hb));
for m=1:size(Hb,3)
   onematb(:,:,m)=eye(size(Hb,1));
end
%% Loop over energy
tic;
de=(emax-emin)/Ne;
epoint=0;
maxerror=-1;
for energy=emin:de:emax %loop over energy
   epoint=epoint+1
   ef(epoint)=energy;
   %% Obtain the self-energy matrix at each Fermi energy
   
   sigmal=Self(Tl1,h_l1,Tr1,energy);
   sigmar=Self(Tr1,h_l1,Tl1,energy);
   gammal=1i*(sigmal-sigmal');
   gammar=1i*(sigmar-sigmar');
   
   Sb(i,i,1)=sigmal;
   Sb(Nsc+i,Nsc+i,1)=sigmar;
   S(i,i,1)=sigmal;
   S(i,i,size(H_c,3))=sigmar;

   %% We can also add impurity to this system
   %Vimp=0.5;
   %Impurity=Vimp*reshape(onematb(:).*rand(numel(onematb),1),[size(onematb,1),size(onematb,2),size(onematb,3)]);
   
   %% Find block of the retarded Green function matrix needed for transission function 
   gb=energy*onematb-Hb-Sb; 
   Gr=invnn(-Tb1,gb,-Tb2,s);
   g_r=Gr(i,Nsc+i,1);
   
%  gb=energy*onematb-Hb-Sb-Impurity;
%  g=invb(-Tb1,gb,-Tb2,s);
%  g_r=g(i,Nsc+i);
   %% Find transmission function via the Caroli NEGF implementation of the Landauer-Buttiker formula    
   trans(epoint)=real(trace(gammal*g_r*gammal*g_r'));
   
   %% find total dos from the imaginary part of the retarded Green function
   dos(epoint)=-1/pi*imag(trace(sum(Gr,3)));

end %loop over energy
maxerror
toc;
figure(3), plot(ef,trans,'Marker','.');
grid;
title('Transmission of 2D conductor vs. Fermi energy');
xlabel('Fermi Energy');
ylabel('Conductance (2e^2/h)');
figure(4), plot(ef,dos,'Marker','.');
grid;
title('GNR Total DOS vs. Fermi energy');
xlabel('Fermi Energy');
ylabel('Density of States');

