MatLab Code for Band Structure of Aluminium

Tuesday, September 28, 2010

Band structure of Al using pseudopotential.
MatLab code for plot band structure of Al



clc
clear
close all

U0 = -31.3; % in eV
d = 0.35 ;% in Angstroms
Rc = 0.943; % in Angstroms
a = 4.05; %in Angstrom
c = 3.8; %in eV*(Angstrom^2) h^2/2m
N=50; % Number of k values in first brillouin zone
bands=6; % Number of bands that need to plot



% Al has fcc structure. reciprocal lattice has bcc structure. Structure factor for fcc will gives the h,k,l values

t=3;
r = 1;
for h = -t:1:t
for k = -t:1:t
for l = -t:1:t

stfact = 1 + exp(i*pi*(h+k)) + exp(i*pi*(h+l)) + exp(i*pi*(l+k)); %structure factor

if real(stfact) == 4
%now defining an array of g-vectors for h,k,l values
g(r,1) = h;
g(r,2) = k;
g(r,3) = l;
r = r + 1;
end

end
end
end



%calculate K and U================================================
for ri=1:length(g)
for rj=1:length(g) %calculate distance between two rows of g
K(ri,rj)=sqrt((g(ri,1)-g(rj,1)).^2+(g(ri,2)-g(rj,2)).^2+(g(ri,3)-g(rj,3)).^2)*2*pi/a;
%calculate pseudopotential ( Marder equation 10.57)
if K(ri,rj)==0 %if K is zero
U(ri,rj)=U0*exp(-Rc/d)*(Rc/d+1);
else
U(ri,rj)=U0*exp(-Rc/d)*(sin(Rc*K(ri,rj))+K(ri,rj)*d*cos(K(ri,rj)*Rc))/(d*K(ri,rj)*(d^2*K(ri,rj)^2+1));
end

end
end
% end of calculation of K and U=====================================



%Calculate Pannel for Gamma to X************************************
B=2*pi/a; % Gamma=(0,0,0) X=(1,0,0) L=1/2(1,1,1) with units of 2*pi/a.
k=B/N;
kx=0:k:B; %array of k


for z=1:N+1
for ri=1:length(g)
for rj=1:length(g)
if K(ri,rj) == 0
E(ri,rj) = c*(((z-1)*k-B*g(ri,1))^2+(B*g(ri,2))^2+(B*g(ri,3))^2) + U(ri,rj);
else
E(ri,rj) = U(ri,rj);
end
end
end
%Ex=energy eigen values for Gamma to X
Ex(z,1:length(g)) = eig(E);
end





%Calculate Pannel for L to Gamma++++++++++++++++++++++++++++++++++++++
Bl=pi/a; % Gamma=(0,0,0) L=1/2(1,1,1) with units of 2*pi/a.
k=Bl/N;
kl=0:k:Bl;


for z=1:N+1
for ri=1:length(g)
for rj=1:length(g)
if K(ri,rj) == 0
E(ri,rj) = c*((Bl-(z-1)*k-B*g(ri,1))^2+(Bl-(z-1)*k-B*g(ri,2))^2 +(Bl-(z-1)*k-B*g(ri,3))^2)+ U(ri,rj);
else
E(ri,rj) = U(ri,rj);
end
end
end
%Ex=energy eigen values for L to Gamma
El(z,1:length(g)) = eig(E);
end


1 comment

Interesting article.
How about calculating fluid mechanism?

December 31, 2010 at 10:55 PM

Post a Comment