EPM Theory Implementation
EPM Theory Implementation
EPM Theory Implementation
Dragica Vasileska
Professor
Arizona State University
Tempe, AZ 85287-5706, USA
1.
Introduction
The basis for discussing transport in semiconductors is the underlying electronic
band structure of the material arising from the solution of the many body Schrdinger
equation in the presence of the periodic potential of the lattice, which is discussed in a
host of solid state physics textbooks. The electronic solutions in the presence of the
periodic potential of the lattice are in the form of Bloch functions
n,k un k eik r
(1)
where k is the wavevector, and n labels the band index corresponding to different
solutions for a given wavevector. The cell-periodic function, u n k , has the periodicity of
the lattice and modulates the traveling wave solution associated with free electrons.
A brief look at the symmetry properties of the eigenfunctions would greatly
enhance the understanding of the evolution of the bandstructure. First, one starts by
looking at the energy eigenvalues of the individual atoms that constitute the
semiconductor crystal. All semiconductors have tetrahedral bonds that have sp
hybridization. However, the individual atoms have the outermost (valence) electrons in sand p-type orbitals. The symmetry (or geometric) properties of these orbitals are made
most clear by looking at their angular parts
s 1
x
3 sin cos
r
y
p y 3 sin sin
r
z
pz 3 cos
r
px
(2)
Lets denote these states by |S>, |X>, |Y> and |Z>. Once one puts the atoms in a crystal,
the valence electrons hybridize into sp3 orbitals that lead to tetrahedral bonding. The
crystal develops its own bandstructure with gaps and allowed bands. For semiconductors,
one is typically worried about the bandstructure of the conduction and the valence bands
only. It turns out that the states near the band-edges behave very much like the |S> and
the three p-type states that they had when they were individual atoms.
principles, i.e. without the need for empirical fitting parameters. In general, these
methods utilize a variational approach to calculate the ground state energy of a manybody system, where the system is defined at the atomic level. The original calculations
were performed on systems containing a few atoms. Today, calculations are performed
using approximately 1000 atoms but are computationally expensive, sometimes requiring
massively parallel computers.
In contrast to ab initio approaches, the second category consists of empirical
methods, such as the Orthogonalized Plane Wave (OPW) [ 2], tight-binding [3] (also
known as the Linear Combination of Atomic Orbitals (LCAO) method), the k p method
[4], and the local [5], or the non-local [6] empirical pseudopotential method (EPM). These
methods involve empirical parameters to fit experimental data such as the band-to-band
transitions at specific high-symmetry points derived from optical absorption experiments.
The appeal of these methods is that the electronic structure can be calculated by solving a
one-electron Schrdinger wave equation (SWE). Thus, empirical methods are
computationally less expensive than ab initio calculations and provide a relatively easy
means of generating the electronic band structure.
2.
(3)
where bk ,t are orthogonalization coefficients and k ,t are core wave functions. For Si14, the summation over t in Eq. (3) is a sum over the core states 1s2 2s2 2p6. Since the
crystal wave function is constructed to be orthogonal to the core wave functions, the
orthogonalization coefficients can be calculated, thus yielding the final expression
k k k ,t k k ,t .
(4)
p2
,
VC
2m
(5)
is applied to Eq. (4), where VC is the attractive core potential, and the following wave
equation results
p2
VC V R k E k ,
2m
(6)
VR
t
E Et
k ,t k k ,t .
k
(7)
Et in Eq. (7) represents the atomic energy eigenvalue, and the summation over t
represents a summation over the core states. The result given in Eq. (6) can be thought of
as wave equation for the pseudo-wave function, k , but the energy eigenvalue E
corresponds to the true energy of the crystal wave function k . Furthermore, as a result
of the orthogonalization procedure, the repulsive potential VR, which serves to cancel the
attractive potential VC, is introduced into the pseudo-wave function Hamiltonian. The
result is a smoothly varying pseudopotential VP = VC + VR. This result is known as the
Phillips-Kleinman cancellation theorem [9] which provides justification why the
electronic structure of strongly-bound valence electrons can be described using a nearlyfree electron model and weak potentials. To simplify the problem further, model
pseudopotenials are used in place of the actual pseudopotential. summarizes the various
models employed. Note that the 3D Fourier transforms (for bulk systems) of each of the
above-described model potentials are of the following general form
V (q) ~
Ze 2
cos qrc .
0q 2
(8)
This q-dependent pseudopotential is then used to calculate the energy band structure
along different crystallographic directions, using the procedure outlined in the following
section.
V ( r)
Ze
;
4 0 r
V (r )
rc
r rC
Ze 2
; r rC
4 0 rC
V ( r)
V (r )
Ze
; r rC
4 0 r
0;
r rC
V ( r)
Ze
; r rC
4 0 r
A;
r rC
V ( r)
Ze
1- exp r rC ; r rC
4 0 r
0;
r rC
V (r )
rc
V (r)
rc
rc
V P k E k ,
2m
(9)
(10)
To simplify further, the inner summation over can be expressed as the total potential, V0,
in the unit cell located at R. Eq. (10) then becomes
V P r V0 r R .
(11)
Because the crystal potential is periodic, the pseudopotential is also a periodic function
and can be expanded into a Fourier series over the reciprocal lattice to obtain
V P r V0 G e iG r ,
(12)
V0 G
r V0 r e iGr
(13)
(14)
where V1 and V2 are the atomic potentials of the cation and anion.
Figure 3. Diamond vs. zincblende lattice and definition of as the distance between two
atoms in a basis of the interpenetrating face centered cubic lattices. If the atoms in the
basis are the same one talks about diamond lattice. When the two atoms in the basis are
not the same then we have zincblende structure.
Substituting Eq. (14) into Eq. (13), and using the displacement property of Fourier
transforms, V0(r) can be recast as
V0 G e iG V1 (G ) e iG V2 (G )
(15)
,
sin
i G
V(G)
A
(16)
where the prefactors are referred to as the symmetric and antisymmetric structure factors.
The form factors above are treated as adjustable parameters that can be fit to
experimental data, hence the name empirical pseudopotential method. For diamondlattice materials, with two identical atoms per unit cell, the VA=0 and the structure factor
is simply cos G . For zinc-blende lattice, like the one in GaAs material system, VA0
and the structure factor is more complicated.
Now with the potential energy term specified, the next task is to recast the
Schrdinger equation in a matrix form. Recall that the solution to the Schrdinger wave
equation in a periodic lattice is a Bloch function, which is composed of a plane wave
component and a cell periodic part that has the periodicity of the lattice, i.e.
k r e ik r u k r e ik r U (G' )e iG' r .
G'
(17)
By expanding the cell periodic part uk(r) of the Bloch function appearing in Eq. (17) into
Fourier components, and substituting the pseudo-wave function k and potential V0 into
the Schrdinger wave equation, the following matrix equation results
10
2 k G 2
E U G V0 G G U G 0 .
2m
G
G'
(18)
The expression given in Eq. (18) is zero when each term in the sum is identically zero,
which implies the following condition
2 k G 2
E U G V0 G G U G 0 .
2m
G'
(19)
In this way, the band structure calculation is reduced to solving the eigenvalue problem
specified by Eq. (19) for the energy E. As obvious from Eq. (17), U G is the Fourier
component of the cell periodic part of the Bloch function. The number of reciprocal
lattice vectors used determines both the matrix size and calculation accuracy.
The eigenvalue problem of Eq. (19) can be written in the more familiar form
HU EU , where H is a matrix, U is a column vector representing the eigenvectors, and
E is the energy eigenvalue corresponding to its respective eigenvector. For the diamond
lattice, the diagonal matrix elements of H are then given by
H i, j
2
k Gi
2m
(20)
H i , j VS G i G j cos G i G j ,
(21)
for i j. Note that the term VS(0) is neglected in arriving at Eq. (20), because it will only
give a rigid shift in energy to the bands. The solution to the energy eigenvalues and
corresponding eigenvectors can then be found by diagonalizing matrix H.
11
12
difference between reciprocal lattice vectors. It can be shown that the square of the
difference between reciprocal lattice vectors will also form the set of integers previously
described. This means that VS is only needed at discrete points corresponding to nearestneighbor sites.
Therefore, its Fourier transform VS(q) is also a continuous function that is shown in
Figure 4. The points corresponding to the first three nearest neighbors are also indicated
on this figure.
q2=8
q2=11
q2=3
q G G'
Recall that the pseudopotential is only needed at a few discrete points along the
V(q) curve. The discrete points correspond to the q2-values that match the integer set
described previously.
13
For q2 = 4, the cosine term in Eq. (21) will always vanish. Furthermore, for
values of q2 greater than 11, V(q) quickly approaches zero. This comes from the fact that
the pseudopotential is a smoothly varying function, and only few plane waves are needed
to represent it. If a function is rapidly varying in space, then many more plane waves
14
[12].
simulations. The lattice constants specified for Si and Ge are 5.43 and 5.65,
respectively. This band-structure is obtained in the following manner.
Figure 5. Flow chart for the implementation of the Empirical Pseudopotential Method.
At the beginning of the simulation, the reciprocal lattice vectors G have to be generated
using the information given in Figure 3. The reciprocal lattice vectors given in Fig. 3 lead
to 89 G-vectors. These are labeled as follows: G1=(0,0,0), G2=(1,1,1), G3=(1,1,-1),
15
Figure 6. First Brillouin zone of FCC lattice that corresponds to the first Brillouin zone
for all diamond and Zinc-blende materials (C, Si, Ge, GaAs, InAs, CdTe, etc.). There are
8 hexagonal faces (normal to [111]) and 6 square faces (normal to [100]). The sides of
each hexagon and square are equal.
16
Si is an indirect band gap semiconductor. Its primary gap, i.e. minimum gap, is
calculated from the valence band maximum at the -point to the conduction band
minimum along the direction, 85% of the distance from to X. The band gap of Si,
using the parameters from Ref. [Error: Reference source not found], is calculated to be
EgSi = 1.08 eV, in agreement with experimental findings. Ge is also an indirect band gap
semiconductor. Its band gap is defined from the top of the valence band at to the
conduction band minimum at L. The band gap of Ge is calculated to be EgGe = 0.73 eV.
The direct gap, which is defined from the valence band maximum at to the conduction
band minimum at , is calculated to be 3.27 eV and 0.82 eV for Si and Ge, respectively.
Note that the curvature of the top valence band of Ge is larger than that of Si. This
corresponds to the fact that the effective hole mass of Si is larger than that of Ge. Note
that the inclusion of the spin-orbit interaction will lift the triple degeneracy of the bands
at the point, leaving doubly-degenerate heavy and light-hole bands and a split-off band
moved downward in energy by few 10's of meV (depending upon the material under
consideration).
17
Figure 7. Left panel: band structures of silicon. Right panel: band structure of
germanium.
3.
This code has been created by Ph.D. Santhosh Krishnan (Micron) who was a Ph.D.
student of Prof. Vasileska.
epm.m
clc;
clear;
load str_fac.txt;
load a0.txt;
Ry2eV=13.6056981d0;
hb=1.054D-34;
h=6.625D-34;
q=1.602d-19;
mo=9.11D-31;
t=hb*hb/2.0d0/mo/q;
18
Num_semi=14;
NN_neighbours=10;
N_pw=137;
N_bands=12;
N_ek=40;
ig2_mag(1:NN_neighbours)=[0, 3, 4, 8, 11, 12, 16, 19, 20, 24];
Lwork=2*N_pw-1;
sprintf('%s','This program calculates the bandstructure of
semiconductors based on the Empirical Pseudopotential Method (EPM)')
sprintf('%s','The following semiconductors are included')
sprintf('1 -> Si 2 -> Ge 3 -> Sn 4 -> GaP 5 -> GaAs 6 -> AlSb 7 -> InP')
sprintf('8 -> GaSb 9 -> InAs 10 -> InSb 11 -> ZnS 12 -> ZnSe 13 -> ZnTe
14 -> CdTe')
ichoose=input('Enter the semiconductor of your choice: ');
ao=a0(ichoose);
Vs3=str_fac(ichoose,1);
Vs8=str_fac(ichoose,2);
Vs11=str_fac(ichoose,3);
Va3=str_fac(ichoose,4);
Va8=str_fac(ichoose,5);
Va11=str_fac(ichoose,6);
ao=ao*1.0d-10;
Vs3=Vs3*Ry2eV;
Vs8=Vs8*Ry2eV;
Vs11=Vs11*Ry2eV;
Va3=Va3*Ry2eV;
Va8=Va8*Ry2eV;
Va11=Va11*Ry2eV;
icount=0;
for i=1:NN_neighbours
ig2=ig2_mag(i);
a=sqrt(double(ig2));
ia=floor(a);
ic=0;
for igx=-ia:ia
igx2=igx*igx;
for jgy=-ia:ia
jgy2=jgy*jgy;
igxy2=igx2+jgy2;
for kgz=-ia:ia
kgz2=kgz*kgz;
igxyz2=igxy2+kgz2;
if (igxyz2==ig2)
ic=ic+1;
igx_vec(ic+icount)=igx;
igy_vec(ic+icount)=jgy;
igz_vec(ic+icount)=kgz;
end
19
end
end
end
%call gen_gvecs(i,ia,ig2,ic)
icount=icount+ic;
sprintf('%d %d %d',ig2,ic,icount);
end
ik_curr=0;
sumx=0.0d0;
X0=0.5d0;
Y0=0.5d0;
Z0=0.5d0;
% print*
% print*,'L -> Gamma'
% print*
for ik=1:N_ek+1
ik+ik_curr
dkx=0.50d0*(1.0d0-double(ik-1)/double(N_ek));
dky=0.50d0*(1.0d0-double(ik-1)/double(N_ek));
dkz=0.50d0*(1.0d0-double(ik-1)/double(N_ek));
generate_Hamil;
sumx=sqrt((X0-dkx)*(X0-dkx)+(Y0-dky)*(Y0-dky)+(Z0-dkz)*(Z0-dkz));
k_vec(ik+ik_curr)=sumx;
end
ik_curr=ik;
% ! Gamma -> X
sumx0=sumx;
X0=0.0d0;
Y0=0.0d0;
Z0=0.0d0;
% print*
% print*,'Gamma -> X'
% print*
for ik=1:N_ek
%print*,ik
ik+ik_curr
dkx=double(ik)/double(N_ek);
dky=0.0d0;
dkz=0.0d0;
generate_Hamil;
sumx=sqrt((X0-dkx)*(X0-dkx)+(Y0-dky)*(Y0-dky)+(Z0-dkz)*(Z0-dkz));
k_vec(ik+ik_curr)=sumx+sumx0;
end
ik_curr=ik+ik_curr;
% ! X -> U
sumx1=sumx;
X0=1.0d0;
Y0=0.0d0;
20
Z0=0.0d0;
% print*
% print*,'X -> U'
% print*
for ik=1:N_ek
%print*,ik
ik+ik_curr
dkx=1.0d0;
dky=double(ik)/double(N_ek)*1.0d0/4.0d0;
dkz=double(ik)/double(N_ek)*1.0d0/4.0d0;
generate_Hamil;
sumx=sqrt((X0-dkx)*(X0-dkx)+(Y0-dky)*(Y0-dky)+(Z0-dkz)*(Z0-dkz));
k_vec(ik+ik_curr)=sumx+sumx1+sumx0;
end
ik_curr=ik+ik_curr;
% ! K -> Gamma
sumx2=sumx;
X0=0.75d0;
Y0=0.75d0;
Z0=0.0d0;
% print*
% print*,'K -> Gamma'
% print*
for ik=1:N_ek
%print*,ik
ik+ik_curr
dkx=0.75d0*(1.0d0-double(ik)/double(N_ek));
dky=0.75d0*(1.0d0-double(ik)/double(N_ek));
dkz=0.0d0;
generate_Hamil;
sumx=sqrt((X0-dkx)*(X0-dkx)+(Y0-dky)*(Y0-dky)+(Z0-dkz)*(Z0-dkz));
k_vec(ik+ik_curr)=sumx+sumx2++sumx1+sumx0;
end
plot(k_vec,ek);
clear all;
generate_Hamil.m
pi=4.0d0*atan(1.0d0);
k0=2.0d0*pi/ao;
H_EPM(1:N_pw,1:N_pw)=complex(0.0d0,0.0d0);
for i=1:N_pw
gx_i=double(igx_vec(i));
gy_i=double(igy_vec(i));
gz_i=double(igz_vec(i));
for j=1:N_pw
21
end
end
gx_j=double(igx_vec(j));
gy_j=double(igy_vec(j));
gz_j=double(igz_vec(j));
if (i==j)
gkx2=(gx_i+dkx)*(gx_i+dkx);
gky2=(gy_i+dky)*(gy_i+dky);
gkz2=(gz_i+dkz)*(gz_i+dkz);
Hreal= t*k0*k0*(gkx2+gky2+gkz2);
H_EPM(i,j)=complex(Hreal,0.0d0);
else
dgx=gx_i-gx_j;
dgy=gy_i-gy_j;
dgz=gz_i-gz_j;
idgx=round(dgx);
idgy=round(dgy);
idgz=round(dgz);
id2=idgx*idgx+idgy*idgy+idgz*idgz;
tau_fac=2.0d0*pi/8.0d0;
if ((id2==3))
V_s=Vs3;
V_a=Va3;
elseif ((id2==8))
V_s=Vs8;
V_a=Va8;
elseif ((id2==11))
V_s=Vs11;
V_a=Va11;
else
V_s=0.0d0;
V_a=0.0d0;
end
ct=cos((dgx+dgy+dgz)*tau_fac);
st=sin((dgx+dgy+dgz)*tau_fac);
Hreal= V_s*ct;
Himag= V_a*st;
H_EPM(i,j)=complex(Hreal,Himag);
end
E=sort(eig(H_EPM));
ek(ik_curr+ik,1:8)=E(1:8)';
semiconductors.txt
Si
Ge
Sn
GaP
GaAs
AlSb
InP
GaSb
InAs
InSb
ZnS
22
ZnSe
ZnTe
CdTe
a0.txt
5.43d0
5.66d0
6.49d0
5.44d0
5.64d0
6.13d0
5.86d0
6.12d0
6.04d0
6.48d0
5.41d0
5.65d0
6.07d0
6.41d0
str_fac.txt
-0.21d0
-0.23d0
-0.20d0
-0.22d0
-0.23d0
-0.21d0
-0.23d0
-0.22d0
-0.22d0
-0.20d0
-0.22d0
-0.23d0
-0.22d0
-0.20d0
0.04d0
0.01d0
0.00d0
0.03d0
0.01d0
0.02d0
0.01d0
0.00d0
0.00d0
0.00d0
0.03d0
0.01d0
0.00d0
0.00d0
0.08d0
0.06d0
0.04d0
0.07d0
0.06d0
0.06d0
0.06d0
0.05d0
0.05d0
0.04d0
0.07d0
0.06d0
0.05d0
0.04d0
0.00d0
0.00d0
0.00d0
0.12d0
0.07d0
0.06d0
0.07d0
0.06d0
0.08d0
0.06d0
0.24d0
0.18d0
0.13d0
0.15d0
0.00d0
0.00d0
0.00d0
0.07d0
0.05d0
0.04d0
0.05d0
0.05d0
0.05d0
0.05d0
0.14d0
0.12d0
0.10d0
0.09d0
0.00d0
0.00d0
0.00d0
0.02d0
0.01d0
0.02d0
0.01d0
0.01d0
0.03d0
0.01d0
0.04d0
0.03d0
0.01d0
0.04d0
plot_bs.m
load
load
load
load
load
load
load
load
load
load
load
load
k_vec
Si
Ge
Sn
GaP
GaAs
AlSb
InP
GaSb
InAs
InSb
ZnS
23
load ZnSe
load ZnTe
load CdTe
figure
Si=Si-Si(251,6);
plot(k_vec,Si)
axis([0 3.28 -5.5 6.5])
figure
Ge=Ge-Ge(251,6);
plot(k_vec,Ge)
axis([0 3.28 -5 7])
figure
Sn=Sn-Sn(251,6);
plot(k_vec,Sn)
axis([0 3.28 -4 6])
figure
GaP=GaP-GaP(251,6);
plot(k_vec,GaP)
axis([0 3.28 -4 7])
figure
GaAs=GaAs-GaAs(251,6);
plot(k_vec,GaAs)
axis([0 3.28 -4 7])
figure
AlSb=AlSb-AlSb(251,6);
plot(k_vec,AlSb)
axis([0 3.28 -4 7])
figure
InP=InP-InP(251,6);
plot(k_vec,InP)
axis([0 3.28 -4 7])
figure
GaSb=GaSb-GaSb(251,6);
plot(k_vec,GaSb)
axis([0 3.28 -3 6.5])
figure
InAs=InAs-InAs(251,6);
plot(k_vec,InAs)
axis([0 3.28 -4 7])
figure
InSb=InSb-InSb(251,6);
plot(k_vec,InSb)
axis([0 3.28 -3 6])
figure
ZnS=ZnS-ZnS(251,6);
plot(k_vec,ZnS)
24
References
25
[]
1999.
2
[]D. J. Chadi and M. L. Cohen, Phys. Stat. Sol. (b), 68 (1975) 405.
[]
[]
[]
10
[]
11
[]
12
[]
13
[]
14
[]
15
[]