% 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

%%
%Static impurity introduces on-site scalar potential
U=zeros(1,Np);
simpsite=50;
%U(simpsite)=U(simpsite)+10.0; %Uncomment this line to introduce scalar impurity potential for in panels (e), (f)
H=H+diag(U); %Add potential on the diagonal of Hamiltonian matrix

%%
%Expand Hamiltonian into orbital x spin space
H=kron(H,eye(2)); %Kronecker product of Hamiltonian matrix and unit 2x2 matrix in spin space

%% 
%Magnetic impurity described by unit classical vector (S_x,S_y,S_z) introduces spin-dependent 
%potentail Jsd*(S_x,S_y,S_z)*(sigma_x,sigma_y,sigma_z)
Jsd=-2.0; %strength of spin-dependent potential
Sx=0; Sy=1; Sz=0;
sigmax=[0 1; 1 0];
sigmay=[0 -i; i 0];
sigmaz=[1 0; 0 -1];
Umag=Jsd*(Sx*sigmax+Sy*sigmay+Sz*sigmaz);
mimpsite=50; %lattice site hosting magnetic impurity
%H(2*mimpsite-1:2*mimpsite,2*mimpsite-1:2*mimpsite)=H(2*mimpsite-1:2*mimpsite,2*mimpsite-1:2*mimpsite)+Umag;

%%
%Diagonalize Hamiltonian matrix
[V,D]=eig(H);
E=sort(diag(D)'); %Take diagonal of D -> make it a row vector -> sort eigenvalues within this vector from smallest to largest

%%
%Find equilibrium density matrix
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

%%
%Compute particle and spin density 
rho=rho/a; %density is "per site", so divide by the lattice spacing to get it in physical units
for site=1:Np
    n(site)=trace(rho(2*site-1:2*site,2*site-1:2*site)*eye(2)); %particle density at site
    spinx(site)=trace(rho(2*site-1:2*site,2*site-1:2*site)*sigmax); %spin density S_x at site
    spiny(site)=trace(rho(2*site-1:2*site,2*site-1:2*site)*sigmay); %spin density S_y at site
    spinz(site)=trace(rho(2*site-1:2*site,2*site-1:2*site)*sigmaz); %spin density S_z at site
end

%%
%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 Eigenvalues and Chemical Potential \mu'); %Use LaTeX command to get Greek letter mu for chemical potential
xlabel('Eigenvalue No.');
ylabel('Energy (eV)');
set(h,'linewidth',[2.0])

figure(2); %Select Figure 2 to plot particle and spin density as a grid of 4 plots 
subplot(2,2,1);
plot(X,n)
grid;
title('n vs. x');
xlabel('X (m)'); 
ylabel('Electron density (m^{-1})'); %Use LaTeX command to get superscript

subplot(2,2,2);
plot(X,spinx)
grid;
title('S_x vs. x'); %Use LaTeX command to get subscripts
xlabel('X (m)');
ylabel('Spin density (m^{-1})');

subplot(2,2,3);
plot(X,spiny)
grid;
title('S_y vs. x'); 
xlabel('X (m)'); 
ylabel('Spin density (m^{-1})');

subplot(2,2,4);
plot(X,spinz)
grid;
title('S_z vs. x'); 
xlabel('X (m)');
ylabel('Spin density (m^{-1})');