%Fig2_5
%Essential Electron Transport for Device Physics
%Tight binding with PT
%N-atom linear chain tight binding nearest neighbor only
%t is matrix element coupling nearest neighbors
clear all;
clf;
%plotting parameters + fontsizes
FS      = 12;	%label fontsize 18
FSN     = 12;	%number fontsize 16
LW      = 1;	%linewidth

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

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

eye=complex(0.0,1.0);%square root minus one
s=char('y','k','r','g','b','m','c');	
t=1.0;              %t is matrix element coupling nearest neighbors
N=7;				%Number of atoms in chain
A=zeros(N,N);       %Initialize Hamiltonian matrix
d=zeros(1,N);       %Diagonal matrix element
offd=t*ones(1,N);	%off-diagonal matrix element
vmax=2.0;          %Maximum value of modulus of imaginary potential

t1=offd(1:N-1);							
Hmatrix=diag(t1,1);				%upper diagonal of Hamiltonian matrix							
Hmatrix=Hmatrix+diag(t1,-1);	%add lower diagonal of Hamiltonian matrix
t2=d(1:N);
Hmatrix1=Hmatrix+diag(t2,0);	%add diagonal of Hamiltonian matrix

xplot=ones(N,1);
for j=0:0.005:vmax
v=real(j);
A(1,1)=eye*v*t;    %put imaginary potential (source and drain)
A(N,N)=-eye*v*t;   %at end of chain
% A(2,2)=-eye*v; %boundary condition put source and drain at end of chain
% A(N-1,N-1)=eye*v;
% A(1,N)=t;	%upper right corner for periodic boundary conditions
% A(N,1)=t;	%lower left corner for periodic boundary conditions
Hmatrix=Hmatrix1+A;	%add boundary conditions
  
[psi,E]=eig(Hmatrix);
EnergyD=diag(E);
[Energy,I]=sort(EnergyD);			%sort eigenvalues

 if j==0
     E1=Energy;                     %save the eigenenergies with v=0
 end
 
figure(1)
plot(j*xplot,real(Energy),'.b');
hold on
plot(j*xplot,imag(Energy),'.r');
end
hold off
xlabel('v/|t|');
ylabel('Real (blue), imaginary (red) energy (E/|t|)');
ttl=['Fig2.5 tight binding with PT boundary term +/- i*v*t. Bulk t=',num2str(t),...
    ', N_{atoms}=',num2str(N)];
title(ttl);

%**************************************************************************
%Plot density of states for v=0 and v=vmax
%**************************************************************************
npoints=50*N;                   %number of points in plot
emax=1.3*max(real(E1));         %max energy in units of t
emin=1.3*min(real(E1));         %min energy in units of t
estep=(emax-emin)/(npoints-1);  %energy step

dos=zeros(1,npoints);   %Set dos  v=vmax vector to zero
dos1=zeros(1,npoints);  %Set dos1 v=0    vector to zero
w=emin:estep:emax;      %Energy array used in plot
gamma = 0.1;			%Energy broadening in units of t
gamma=gamma*abs(t);     %use units of t

	for ix=1:N
		Ek(ix)=real(Energy(ix)); %Eigenenergies with v=vmax
        E1(ix)=real(E1(ix));     %Eigenenergies with v=0
gamma22=gamma+imag(Energy(ix));const22=abs(gamma22)/pi;gamma22=(abs(gamma22)/2)^2;
gamma12=gamma+imag(E1(ix));const12=abs(gamma12)/pi;gamma12=(abs(gamma12)/2)^2;
        for j=1:npoints-1
dos(j) =dos(j) +const22/((w(j)-Ek(ix))^2+gamma22); %Lorentzian broadening
dos1(j)=dos1(j)+const12/((w(j)-E1(ix))^2+gamma12); %Lorentzian broadening
        end
	end
   
figure(2);
plot(w,dos,'r');
hold on;
plot(w,dos1,'b');
hold off;
ttl=['Nearest neighbor tight binding, t=',num2str(t),...
    ', \gamma=',num2str(gamma),'\times|t|, N_{atoms}=',num2str(N),...
    ', blue v=0, red v=',num2str(vmax)];
title(ttl);
xlabel('Energy, E/|t|'),ylabel('Density of states, N(E)');