EPM Theory Implementation

Download as doc, pdf, or txt
Download as doc, pdf, or txt
You are on page 1of 27

Empirical Pseudopotential Method:

Theory and Implementation Details

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.

Figure 1: The typical bandstructure of semiconductors. For direct-gap semiconductors,


the conduction band state at k=0 is s-like. The valence band states are linear
combinations of p-like orbitals. For indirect-gap semiconductors on the other hand, even
the conduction band minima states have some amount of p-like nature mixed into the slike state.
Electronic band structure calculation methods can be grouped into two general
categories [1]. The first category consists of ab initio methods, such as Hartree-Fock or
Density Functional Theory (DFT), which calculate the electronic structure from first
3

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.

The Empirical Pseudopotential Method


The concept of pseudopotentials was introduced by Fermi [ 7] to study high-lying

atomic states. Afterwards, Hellman proposed that pseudopotentials be used for


calculating the energy levels of the alkali metals [ 8]. The wide spread usage of
pseudopotentials did not occur until the late 1950s, when activity in the area of

condensed matter physics began to accelerate. The main advantage of using


pseudopotentials is that only valence electrons have to be considered. The core electrons
are treated as if they are frozen in an atomic-like configuration. As a result, the valence
electrons are thought to move in a weak one-electron potential.
The pseudopotential method is based on the orthogonalized plane wave (OPW)
method due to Herring [Error: Reference source not found]. In this method, the crystal
wavefuntion k is constructed to be orthogonal to the core states. This is accomplished
by expanding k as a smooth part of symmetrized combinations of Bloch functions k
, augmented with a linear combination of core states. This is expressed as
k k bk ,t k ,t ,

(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)

To obtain a wave equation for k , the Hamiltonian operator


H

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)

where VR represents a short-range, non-Hermitian repulsion potential, of the form

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)

(a) Constant effective potential in the core region:

Ze
;
4 0 r

V (r )

rc

r rC

Ze 2
; r rC
4 0 rC

V ( r)

(b) Empty core model:

V (r )

Ze
; r rC
4 0 r
0;
r rC
V ( r)

(c) Model potential due to Heine and Abarenkov:

Ze
; r rC
4 0 r
A;
r rC

(d) Lin and Kleinman model potentials:

V ( r)

Ze
1- exp r rC ; r rC
4 0 r
0;
r rC

V (r )

rc

V (r)

rc

rc

Figure 2. Various model potentials.

2.1 Description of the Empirical Pseudopotential Method


Recall from the previous section that the Phillips-Kleinman cancellation theorem
provides a means for the energy band problem to be simplified into a one-electron-like
problem. For this purpose, Eq. (6) can be re-written as
p2

V P k E k ,
2m

(9)

where VP is the smoothly varying crystal pseudopotential. In general, VP is a linear


combination of atomic potentials, Va, which can be expressed as summation over lattice
translation vectors R and atomic basis vectors to arrive at the following expression
V P r Va r R .
R

(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)

where the expansion coefficient is given by

V0 G

r V0 r e iGr

(13)

and is the volume of the unit cell.


To apply this formalism to the zincblende lattice, it is convenient to choose a twoatom basis centered at the origin (R = 0). If the atomic basis vectors are given by 1 = =
- 2, where , the atomic basis vector, is defined in terms of the lattice constant a0 as =
a0(1/8,1/8,1/8), V0(r) can be expressed as
V0 r V1 r V2 r ,

(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)

Writing the Fourier coefficients of the atomic potentials in terms of symmetric


(VS(G)=V1+V2)) and antisymmetric (VA(G)=V1-V2)) form factors, V0(G) is given by
V0 G cos G V(G)
S

,
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)

for i = j, and the off-diagonal matrix elements of H are given by

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

2.2 Implementation of the Empirical Pseudopotential Method for Si and Ge


For a typical semiconductor system, 137 plane waves are sufficient, each
corresponding to vectors in the reciprocal lattice, to expand the pseudopotential. The
reciprocal lattice of a face-centered cubic (FCC), i.e. diamond or zinc-blende structure, is
a body-centered cubic (BCC) structure. So these reciprocal lattice vectors correspond to
a body centered cubic lattice. Reciprocal lattice vectors up to and including the 10 thnearest neighbor from the origin are usually considered which results in 137 plane waves
for the zinc-blende structure. Smaller number of nearest neighbors (in other words plane
waves) gives less accurate results.

Figure 3. Reciprocal lattice vectors up to 8 th nearest neighbor that need to be entered in


the code to get first order approximation of the bandstructure. (1,1,1) has eight
permutations since there are eight possible G-vectors: (1,1,1), (1,1,-1), (1,-1,1), (-1,1,1),
(1,-1,-1), (-1,-1,1), (-1,1,-1), (-1,-1,-1).
The square of the distance from the origin to each equivalent set of reciprocal
lattice sites is an integer in the set |G2| = 0, 3, 4, 8, 11, 12, where |G2| is expressed in
units of (2/ao)2. Note that the argument of the pseudopotential term VS in Eq. (21) is the

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.

The pseudopotential, on the other hand, is a continuous quantity.

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

Figure 4. Fourier transform of the pseudopotential. (Note that

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

Table 1: Local Pseudopotential Form Factors.

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

would be required. Another advantage of the empirical pseudopotential method is that


only three parameters are needed to describe the band structure of non-polar materials.
Using the form factors listed in Table 1, where the Si form factors are taken from
[10] and the Ge form factors are taken from [ 11], the band structures for Si and Ge are
plotted in Figure

[12].

Note that spin-orbit interaction is not included in these

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

G4=(1,-1,1), G5=(-1,1,1), G6=(1,-1,-1), G7=(-1,-1,1), G8=(-1,1,-1), G9=(-1,-1,-1)


G89=(-3,-3,-1). Once the G-vectors are being generated, one is ready with the calculation
of the bandstructure. In this procedure, illustrated in Figure 5, first one defines a k-vector
based on the information along which high-symmetry line we want to plot the
bandstructure. The high symmetry points and the corresponding high symmetry lines are
given in Figure 6. Then, the elements of H ij are calculated where i=1,2, 89 and j=1,2,
89 according to the formulas (17-18). After the matrix H is constructed, the eigenvalue
solver in MATLAB is called using eigen(H). 89 eigenvalues will be the solution and
these are sorted in ascending order and the first 8 (lowest 8) are kept and stored. Then, if
all high-symmetry directions have been visited the process is terminated and one plots the
eigenvalues as a function of k. If not, k is incremented along the same high symmetry
direction and the process is repeated until complete. The high symmetry points for zincblende materials are given in Figure 6.

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.

In summary, the local empirical pseudopotential method described in this section


is rather good for an accurate description of the optical gaps. However, as noted by
Chelikowsky and Cohen [13], when these local calculations are extended to yield the
valence-band electronic density of states, the results obtained are far from satisfactory.
The reason for this discrepancy arises from the omission of the low cores in the
derivation of the pseudopotential in the previous section. This, as previously noted,
allowed the usage of a simple plane wave basis. To correct for the errors introduced, an
energy-dependent non-local correction term is added to the local atomic potential. This
increases the number of parameters needed but leads to better convergence and more
exact band-structure results [14,15].

3.

Source Code of the Local Empirical Pseudopotential Method

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

axis([0 3.28 -3 10])


figure
ZnSe=ZnSe-ZnSe(251,6);
plot(k_vec,ZnSe)
axis([0 3.28 -3 9])
figure
ZnTe=ZnTe-ZnTe(251,6);
plot(k_vec,ZnTe)
axis([0 3.28 -3 9])
figure
CdTe=CdTe-CdTe(251,6);
plot(k_vec,CdTe)
axis([0 3.28 -3 8])

References

25

[]

P. Y. Yu and M. Cardona, Fundamentals of Semiconductors, Springer-Verlag, Berlin,

1999.
2

[]C. Herring, Phys. Rev., 57 (1940) 1169.

[]D. J. Chadi and M. L. Cohen, Phys. Stat. Sol. (b), 68 (1975) 405.

[]J. Luttinger and W. Kohn, Phys. Rev., 97 (1955) 869.

[]M. L. Cohen and T. K. Bergstresser, Phys. Rev., 141 (1966) 789.

[]J. R. Chelikowsky and M. L. Cohen, Phys Rev. B, 14 (1976) 556.

[]

E. Fermi, Nuovo Cimento, 11 (1934) 157.

[]

H. J. Hellman, J. Chem. Phys., 3 (1935) 61.

[]

J. C. Phillips and L. Kleinman, Phys. Rev., 116 (1959) 287.

10

[]

J. R. Chelikowsky and M. L. Cohen, Phys Rev. B, 10 (1974) 12.

11

[]

L. R. Saravia and D. Brust, Phys. Rev., 176 (1968) 915.

12

[]

S. Gonzalez, Masters Thesis, Arizona State University, 2001.

13

[]

J. R. Chelikowsky and M. L. Cohen, Phys. Rev. B, 10 (1974) 5059.

14

[]

K. C. Padney and J. C. Phillips, Phys. Rev. B, 9 (1974) 1552.

15

[]

D. Brust, Phys. Rev. B, 4 (1971) 3497.

You might also like