%Fig2_3b.m
%Essential Electron Transport for Device Physics
%Calculation of 3D cubic lattice density of states (dos) from nearest 
%neigbor and next nearest nieghbor tight-binding dispersion relation

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);

natoms =100;	%number of lattice points (atoms) in each direction 100
npoints=300;	%number of points in plot 100
emax=7.5;		%max energy in units of t
emin=-7.5;		%min energy in units of t
if emax<=emin;error('Error: require emax > emin');end;
   estep=(emax-emin)/(npoints-1);

dos=zeros(1,npoints);	%set dos vector to zero
w=emin:estep:emax;		%energy 
t=-1;                   %set energy scale - for s-orbital, + for p-orbital
tprime=0.1*t;			%tprime
gamma = 0.1;			%Energy broadening in units of t
gamma=gamma*abs(t);     %use units of t
gamma2=(gamma/2)^2;
const=gamma/pi;

kx0=pi/natoms;			%kx constant
ky0=pi/natoms;			%ky constant
kz0=pi/natoms;			%kz constant
kx=0:kx0:kx0*natoms;
ky=0:ky0:ky0*natoms;
kz=0:kz0:kz0*natoms;
	for ix=1:natoms+1 
	for iy=1:natoms+1
	for iz=1:natoms+1
%Dispersion for the cubic lattice
        Ekprime=(t*cos(kx(ix)))+(t*cos(ky(iy)))+(t*cos(kz(iz)))+...
            (2*tprime*cos(kx(ix))*cos(ky(iy)))+...
            (2*tprime*cos(ky(iy))*cos(kz(iz)))+...
            (2*tprime*cos(kz(iz))*cos(kx(ix)));
        Ekprime=2.0*Ekprime;
        
for j=1:npoints
dos(j)=dos(j)+const/((w(j)-Ekprime)^2.0+gamma2); %Lorentzian broadening
end
	end
	end
	end
figure(1);
plot(w,dos,'b','LineWidth',LW);
ttl=['\rmFig2.3b, tight-binding density of states, t=',num2str(t),', t^\prime=',...
  num2str(tprime),', \Gamma=',num2str(gamma),'\times|t|, n_{atoms}=',num2str(natoms)];
title(ttl);
xlabel('Energy, E (t)'),ylabel('Density of states, N(E)');