%Fig2_2b.m
%Essential Electron Transport for Device Physics
%Calculation of 2D square 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
npoints=300;			%number of points in plot
emax=6.5;				%max energy in units of t
emin=-6.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
kx=0:kx0:kx0*natoms;
ky=0:ky0:ky0*natoms;
Ekprime=zeros(natoms,natoms);
	for ix=1:natoms+1
	for iy=1:natoms+1
% dispersion for the square lattice
        Ekprime(ix,iy)=(2.0*t*cos(kx(ix)))+(2.0*t*cos(ky(iy)))+...
(2.0*tprime*cos(kx(ix))*cos(ky(iy)))+(2.0*tprime*cos(ky(iy))*cos(kx(ix)));
        for j=1:npoints
			dos(j)=dos(j)+const/(((w(j)-Ekprime(ix,iy))^2.)+gamma2); %Lorentzian broadening
        end
	end
	end
figure(1);
plot(w,dos,'b','LineWidth',LW);
ttl=['\rmFig2.2b, 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)');

figure(2);
surf(Ekprime);
ttl=['\rmFig2.2b, tight-binding dispersion, \Gamma =',num2str(gamma),...
    '\times|t|, n_{atoms}=',num2str(natoms)];
shading interp;
title(ttl);
xlabel('Wave vector, k_x (\Deltax \pi/L)');
ylabel('Wave vector, k_y (\Deltax \pi/L)');
zlabel('E_k');
