%Fig2_4a
%Essential Electron Transport for Device Physics
%TightBinding_with_boundary_potential.m
%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 7
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.5;          %Maximum value of end 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.01:vmax
v=real(j)*t;
%v=0;
A(1,1)=+v;
A(N,N)=+v;% -v
Hmatrix=Hmatrix1+A;	%add boundary terms
  
[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('Boundary hopping integral, v/|t|');
ylabel('Real (blue), imaginary (red) energy (E/|t|)');
grid on
ttl=['Fig2.4a tight binding with boundary hopping integral v. Bulk t=',num2str(t),...
    ', N_{atoms}=',num2str(N)];
title(ttl);

%**************************************************************************
%Plot density of states for v=0 and v=vmax
%**************************************************************************
npoints=200*N;                   %number of points in plot
emax=1.3*max(real(Energy));     %max energy in units of t
emin=1.3*min(real(Energy));     %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); %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;
grid on
ttl=['\rmFig2.4a 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)');

figure(3)
plot(diag(Hmatrix))
grid on
xlabel('Atom position, x/L'),ylabel('v');

figure(4)
plot(real(psi))
grid on
xlabel('Atom position, x/L'),ylabel('Re(\psi(x/L))');