% Program to compute equilibrium electron density for a 1D nanowire modeled
% on a discrete lattice of points and with possible localized impurity potential
% within the wire

clear all %Clear all variables to avoid carrying over their values from previous calculations stored in the MATLAB workspace

%%
%Constants (all MKS, except energy which is in eV)
hbar=1.055e-34; %Planck constant
m=9.110e-31; %Electron mass
q=1.602e-19; %Electron charge
mu=0.25; %Chemical potential
kT=0.025; % 0.025 for Part (c),(e) and 0.0025 for Part (d),(f)

%%
%Lattice
Np=100; %Number of discrete lattice point along 1D box
a=2e-10; %Lattice spacing in meters; this means that we assume that the nanowire is 20 nm long
X=a*[1:1:Np]; %x coordinate in physical units
t0=(hbar^2)/(2*m*(a^2))/q; % off-diagonal matrix element of the discretized Hamiltonian 

%%
%Hamiltonian matrix
H(Np,Np)=0; %Set all elements of the Hamiltonian matrix of size Np x Np to zero
H=(2*t0*diag(ones(1,Np)))-(t0*diag(ones(1,Np-1),1))-(t0*diag(ones(1,Np-1),-1)); %MATLAB shortcut to fill the elements of discretized Hamiltonian matrix

%%
%Conventional programming (as employed in Fortran, C, C++) to fill the
%elements of the Hamiltonian matrix requires to loop over matrix elements
%of this tridiagonal sparse matrix (for non-sparse matrix two nested for
%loops are required, in general, to fill all elements of the matrix)
%H(1,1)=2*t0;
%H(Np,Np)=2*t0;
%for j=2:Np %loop over sites
%   H(j,j)=2*t0;
%   H(j,j-1)=-t0;
%   H(j-1,j)=-t0;
%end %loop over sites

%%
%H(1,Np)=-t0; H(Np,1)=-t0; %Uncomment this line to use periodic boundary conditions in panels (e), (f)
U=zeros(1,Np); 
%U(Np/2)=U(Np/2)+10.0; %Uncomment this line to introduce impurity potential for in panels (e), (f)

%%
%Diagonalize the Hamiltonina matrix and find the density matrix
[V,D]=eig(H+diag(U));
E=sort(diag(D)'); %Take diagonal of D -> make it a row vector -> sort eigenvalues within this vector from smallest to largest
D=diag(D)-mu; %Subtract chosen chemical potential value from each energy eigenvalue 
rho=1./(1+exp(D./kT)); %Compute the equilibrium density matrix in the basis of energy eigenstates 
rho=V*diag(rho)*V'; %Unitary transformation to get the density matrix in the original basis
rho=diag(rho)/a; %Electron density is given by the diagonal elements of rho; density is "per site", so divide by the lattice spacing to get it in physical units


%%
%Plotting the results:
figure(1);%Select Figure 1 to produce panel (a) in Fig. 4.3.1 of the textbook
h=plot(E); %Plot elements of the row vector E (energy eigenvalues) vs. their index 
grid;
hold on %Add another plot on the same figure
h1=plot(mu*ones(Np,1),'b'); %Plot chemical potential as a straight line
hold off
title('Energy eigenvalus and chemical potenial \mu'); %Use LaTeX command to get Greek letter mu for chemical potential
xlabel('Eigenvalue Number');
ylabel('Energy (eV)');
set(h,'linewidth',[2.0])


%%
figure(2); % Select Figure to produce panels (c)-(f) by redoing caculations above with different parameters
plot(X,rho)
grid;
title('n(x) vs. x for k_BT=0.025 eV and hard wall B.C.'); %Use LaTeX command to get subscripts
xlabel('X Coordinate (m)'); %Use LaTeX command to get superscript
ylabel('Electron density (m^{-1})');

%You can either reuse figure(2), by chaning its labels and title, or copy figure(2)
%with new labels, title, and number to plot the results for other parameters (such as B.C., temperature, impurity position and strength, ...)  
