Ihm 1988

Home Search Collections Journals About Contact us My IOPscience
Total energy calculations in solid state physics
This content has been downloaded from IOPscience. Please scroll down to see the full text.
1988 Rep. Prog. Phys. 51 105
(http://iopscience.iop.org/0034-4885/51/1/003)
View the table of contents for this issue, or go to the journal homepage for more
Download details:
IP Address: 128.122.253.228
This content was downloaded on 15/02/2015 at 18:06
Please note that terms and conditions apply.

Rep. Prog. Phys. 51 (1988) 105-142. Printed in the U K
Total energy calculations in solid state physics
J Ihmt
Bell Communications Research, Murray Hill, New Jersey 07974, USA
Abstract
Total energy calculations for the study of atomic and electronic structures of solids
are reviewed. A history of total energy calculations in solid state physics from the
emergence of quantum physics to the mid 1970s is briefly summarised. Important
developments in the last decade, the period during which computing capability has
grown explosively all over the world, are then described. Modern computational
method are discussed in detail, with emphasis on the tight-binding, quantum-chemical
cluster, and pseudopotential methods. Applications of the total energy method to
various properties of solids, including the cohesive energy, lattice constant, bulk
modulus, crystal structure and its phase transitions, phonons, surfaces, chemisorption,
interfaces, superlattices, and defects are described for various materials with particular
attention to semiconductors. Discussion of outstanding problems and concluding
remarks concerning prospects for future developments are given.
This review was received in its present form in July 1986.
t Present address: Department of Physics, Seoul National University, Seoul 151, Korea.
0034-4885/88/010105 +38$09.00 @ 1988 IOP Publishing Ltd 105

106 J Ihm
Contents
Page
1. Introduction 107
2. Formalism 109
2.1. Tight-binding method 109
2.2. Quantum-chemical cluster calculations 112
2.3. Pseudopotential method 115
2.4. Other theoretical tools 121
3. Applications 122
3.1. Structurally related properties of bulk materials 122
3.2. Phonons 124
3.3. Surface 125
3.4. Chemisorption, interface and superlattice 133
3.5. Defects and others 137
4. Discussion and conclusions 138
Acknowledgments 139
References 140
Total energy calculations 107
1. Introduction
The total energy calculation is not a new subject in the structural study of solids.
Although the total energy calculation of solids has received renewed attention because
of the growth in computing capability in recent years, the history of this field actually
dates back to the early 1930s. After quantum physics was firmly established in the late
1920s, mainly through the study of atomic physics, physicists tried to apply this new
physics to the solid state of materials. Perhaps the most fundamental question in solid
state physics is the origin of the cohesion of (gaseous) atoms to form a periodic
structure of solids. A pioneering quantitative study of the cohesive energy of alkali
metals (sodium and lithium) was done successfully by Wigner and Seitz (1933, 1934,
Seitz 1935). The band structure effect from the periodic potential was approximated
by a spherical ‘Wigner-Seitz’ cell. The ‘exhange hole’ effect originating from the
antisymmetric electron wavefunctions was correctly treated for the first time. Fuchs
(1935) followed the same approach to expand the work to the case of a noble metal,
copper. Some of the earlier attempts to understand the cohesion arid other structural
properties of metals were summarised by Brooks (1963).
For the case of molecular crystals (including crystallised noble gases) and ionic
crystals, a more empirical method has been very useful and well accepted. Here, one
assumes that the atoms or molecular units are classical particles localised at the
equilibrium position. The equilibrium position is determined by the balance between
the interatomic exclusion and the electric attraction. For molecular crystals like solid
argon, the origin of the attraction is the Van der Waals force due to the spontaneous
(dynamic) polarisation of the atomic or molecular units whereas for ionic crystals like
sodium chloride it is the Coulomb force between ions (in the presence of the electronic
background). The electron kinetic energy is a very small correction to this picture for
these materials. Since the Coulomb interaction is of long range, it is important to sum
the electrostatic energy of the lattice (Madelung energy) properly. The techniques for
the lattice energy summations have been developed by various authors (Ewald 1921,
Born and Misra 1940, Coldwell-Horsfall and Maradudin 1960). These results are
frequently used not only for ionic crystals but for metals and semiconductors as well.
A very serious problem exists, on the other hand, in the understanding of the
cohesion of the covalent (directionally bonded) semiconductors. None of the above
approximations are valid for covalent bonds being present in many semiconductors
and the quantitative understanding of the large cohesive energy in materials like silicon
remained an unsolved problem until recently. It seems that an actual band structure
calculation of the crystal is required to get reasonable results for the cohesive energy
and other structural properties of covalent materials. As a matter of fact, in order to
get accurate cohesive energies, one has to resort to full band structure calculations for
any kind of crystal. The importance of the correction of the electron eigenvalues due
to the periodic potential (the band structure effect) has been recognised from the
beginning of the quantum theory of solid state physics, and the cellular method
previously mentioned (Wigner and Seitz 1933) is in fact an approximate solution to
this complicated problem. Since this method was not satisfactory for most materials
other than alkali metals, more and more sophisticated methods were developed to
108 J Ihm
calculate the band structure of solids accurately. In fact, the development of the total
energy calculation is closely related to that of band structure computational methods.
(Of course, the total energy calculation also invokes other complicated problems
mentioned later and the development of the total energy calculation always falls behind
the corresponding development of the band structure calculation.) These are the
orthogonalised plane waves method (Herring 1940), the pseudopotential method (Phil-
lips and Kleinman 1959, Cohen and Heine 1970), the tight-binding or the linear
combination of atomic orbitals (LCAO) method (Slater and Koster 1954), the augmented
plane waves (APW) method (Slater 1937), the Green function ( K K R ) method (Korringa
1947, Kohn and Rostoker 1954), the k * p method (Kane 1966), and the linear muffin-tin
orbitals (LMTO) method (Andersen 1975).
Another essential feature in obtaining a correct cohesive energy of solids is the
adequate treatment of the exchange and the correlation energy. The exchange energy
originates from the Pauli exclusion principle between same-spin electrons and the
correlation energy results from many-body effects (mainly the Coulomb interaction
correction between opposite-spin electrons). They do not have an analogue in classical
physics. In fact, a significant part of the band structure energy is the exhange-correlation
energy, which is discussed separately here only for the convenience of presentation.
Important progress in this field was made when Hohenberg and Kohn (1964) proved
that the ground state of a system (atoms, molecules, or solids) is uniquely specified
by the electronic density distribution in space. From that theorem, a useful method
was extracted (Kohn and Sham 1965) to treat the exchange-correlation energy in terms
of a potential depending on the local density only. Slater and Johnson (1972) later
devised a statistical exchange method, the so-called X a method, with tremendous
reduction in computational effort. A separate correlation energy term has been investi-
gated intensively by other authors (Wigner 1934, Singwi et a1 1970, Hedin and Lundquist
1971, von Barth and Hedin 1972, Gunnarsson et a1 1974).
With the advance of the band structure computational techniques and of the
theoretical understanding of the exchange-correlation terms, more complicated
materials were being handled with improving accuracy by the mid 1970s. Total energy
calculations within the muffin-tin approximation (the A P W - X ~
method to be precise)
were reported by Averill (1972) on alkali metals, by Snow (1973) on copper, by Sabin
et a1 (1975) on neon, and by Janak et a1 (1975) and Janak and Williams (1976) on a
large number of metals including transition metals. This approach, however, is usually
confined to the study of close-packed materials for which the muffin-tin model of the
charge density (i.e. the spherical atomic charge approximation inside a sphere) is valid.
Generalising this method to include non-muffin-tin corrections seems to be very
complicated (Danese 1974). Usually the calculation of the cohesive energy becomes
involved by the need to treat accurately core energies which constitute a large part of
the total energy (e.g., more than 99.9% in tungsten) but introduce a negligible contribu-
tion to the cohesive energy. If the muffin-tin approximation is assumed, it is possible
to use a direct algebraic cancellation (Janak 1974) of these large energies. However,
when the charge density is not spherically approximated, the all-electron calculation
must be done with extreme accuracy. Despite this difficulty, self-consistent all-electron
calculations of the total energy without muffin-tin approximations have been performed
on a number of systems. These include the Hartree-Fock calculations on the atomic
hydrogen crystal (Harris and Monkhorst 1971) and diamond (Goroff and Kleinman
1970, Wepfer et a1 1974), density functional calculations on Tis2, LiF and BN (Zunger
and Freeman 1977a, b, 1978), and LCAO calculations on lithium (Ching and Callaway
1974). The pseudopotential method was also used in structure calculations of Si
(Weaire 1970) and several group IV and III-V semiconductors (Morita et a1 1972).
The pseudopotential theory was advanced and made applicable to a broader class of
materials by Moriarty (1974, 1977). The use of a model pseudopotential (Appelbaum
and Hamann 1973) was extended to the study of phonon modes of Si (Wendel and
Martin 1978). There was also a self-consistent pseudopotential calculation made by
means of a variational approximation of the valence band Wannier functions (Verges
and Tejedor 1979).
So far, the history of the total energy calculation in the structural study of solids
has been briefly reviewed up to the mid 1970s. (Works on rare-earth materials will
not be covered at all here.) This review is inevitably biased, although efforts have been
made to be as impartial as possible, due to the author’s limited knowledge of the field.
It is the author’s belief that this field received renewed attention by solid state physicists
in the late 1970s and progress has been tremendously accelerated since then because
of a rapid development in large-scale computers. In parallel with this development
interest in total energy calculations has shifted towards more complicated problems
such as phonons, structural changes of bulk materials, surfaces, interface, defects, etc.
Since the band structure energy is an essential part of the total energy of solids and
any approximate treatment of this term can lead to serious errors, it is nowadays a
common practice to compute the band structure (i.e. all the necessary eigenvalues of
the system) however co’mplicated the system may be. Therefore, one can imagine a
one-to-one correspondence between a scheme of a total energy calculation and that
of a band structure calculation. Such a correspondence does exist, but three methods
(the tight-binding, the quantum-chemical cluster, and the pseudopotential method)
seem to be used more widely than others for the purpose of the current total energy
calculations. There now exists a vast literature on total energy calculations of solids,
and it is impossible to cover them all in the present review. Some of the more popular
calculational techniques mentioned above are described here, with an emphasis on
the pseudopotential method in which the author’s expertise consists. For the description
of other techniques, references are provided for readers’ convenience. It is to be noted
that the treatment of the exchange-correlation energy has also improved remarkably
in recent years. However, since it is a complex subject, only brief remarks are given
in the final section.
The review is organised as follows. Theoretical methods are reviewed in 0 2,
important results obtained from these theoretical calculations are presented in 0 3 and
discussion and prospects for the future are given in 0 4.
2. Formalism
2.1. Tight-binding method

In 1978, a new scheme was proposed (Chadi 1978) to calculate the total energy of
solids based on tight-binding methods. Using this scheme, Chadi obtained the structure
of the GaAs ( 1 10) surface, giving excellent agreement with experiment. Since only a
relatively small amount of computation is involved, this method has proved useful to
study (and predict) new and complex structures where available information on the
structure is scarce. The input into the calculation is obtained from bulk electronic and
lattice dynamic properties. This approach employs a semiempirical model for the
variation of the total energy AE,,, with changes in atomic coordinates. In the bulk,
110 J Ihm
the elastic coefficients and phonon frequencies can be directly related to the variation
in E,,, with atomic displacements. Therefore, the accuracy of the model for AE,,, can
be easily tested. The formalism is derived below following the origional work by Chadi.
The total energy E,,, of an electron-ion system can be expressed as
Et,,= E,,+E,i+ Eii (2.1)

where E,, , Eeiand Eiidenote the electron-electron, electron-ion and ion-ion interaction
energies. It is useful to introduce the ‘band structure’ energy defined as
where the sum of single particle energies is taken over occupied states, i. In terms of
Ebs, equation (2.1) can be expressed as
= Ebs + ( Eii - (2.3)

making use of the fact that Ebs= Eei+ 2E,, (the factor 2 comes from the double counting
of the electron-electron Coulomb energy). The advantage of (2.3) over (2.1) is that
for two ions that are separated by a distance much larger than the Thomas-Fermi
screening length the combined ion plus screening electron system is nearly neutral and
( Eii- E,,) is close to zero. One would therefore expect that, to a good approximation,
this term can be described by a short-range force constant model.
To calculate E b s the Slater-Koster (1954) tight-binding method was used. We then
impose two conditions on E,,,, namely,
aE,,,/a V = 0 (2.4)
and
Va2E,,,/a V2 = B
where V denotes volume and B the bulk modulus. For a d-’ dependence of the
tight-binding parameters on nearest-neighbour distance d, the above equations can be
satisfied by taking
AEt,, = +U
with
U= C ( U,eij+U , e i )
i>j
where e,j = Adij/d,j is the fractional change in the distance d,, between two nearest-
neighbour atoms denoted by i and j . Once U , and U, are determined for a particular
material of interest (table 1) through this process, we can calculate, for example, the
frequency of phonon modes or determine the energy minimum configuration of surfaces
using the same parameters.
Table 1. Force constant parameters for Si, GaAs and ZnSe from Chadi (1978).
Si GaAs ZnSe
U, -23.19 -17.72 -14.15

U* 53.35 44.80 36.61
Without any additional parameters the above model for AE,,, also predicts the
correct angular-restoring forces, i.e. for the tight-binding parameters used in the
calculations AE,,, = AEbsfor atomic displacements which involve small (second-order)
variations in bond length. The elastic coefficients C,, - CI2and C44,and the zone
boundary transverse acoustic (TA) mode at X are related to angular restoring forces.
The agreement between the calculated and experimental values of these quantities
shown in table 2 is sufficiently close that no additional force constant parameters were
used. As discussed above the parameters U , and U, in equation (2.7) were determined
by considering the variation of E,,, with volume. To test the accuracy of the parameters
in predicting AE,,, for other situations the transverse optic (TO) mode at r has also
been calculated. This mode involves changes in both bond lengths and angles but
none in volume. The calculated and experimental values of the TO frequencies for Si,
GaAs and ZnSe shown in table 2 agree to within 4%. The model used in calculating
AE,,, is therefore sufficiently accurate for considering general types of atomic motion.
Table 2. The cubic lattice constant and bulk modulus B were fitted to the values shown
below. The values of some bulk elastic coefficients and phonon frequencies calculated
from E,,, are compared to experimental values (in parentheses). These results are repro-
duced from Chadi (1978).
Si GaAs ZnSe
a (A) 5.44 5.65 5.65

E (10" erg 9.78 7.48 5.95
c,,- clz(IO" erg ~ m - ~ ) 10.76 (10.18) 6.43 (6.49) 3.34 (3.22)
C , (10" erg ~ m - ~ ) 9.21 (7.96) 6.35 (5.92) 3.50 (4.41)
TA(X) (10" Hz) 6.66 (4.49) 3.12 (2.39) 2.38 (2.10)
TO(r) (10" Hz) 15.6 (15.6) 8.46 (8.1) 6.28 (6.21)
Since bond lengths are not constrained to remain fixed in the calculations of most
phonon modes or surface structures, the adequacy of the harmonic expansion in ( 2 . 7 )
has been further tested. For Si the coefficient U, of E ; in (2.7) obtained from the
results in the literature (Keating 1966, Wendel and Martin 1978) is found to be around
-166 eV. This indicates that the fractional change in bond lengths has to be as large
as 32% before the cubic term gives a contribution equal to that of the quadratic term,
hence the cubic term is negligible for practical purposes.
Application of the method to the determination of the energy minimum structure
of the GaAs (110) surface will be given in the next section. Because the method uses
empirical parameters, there is a limitation on the accuracy of the results. For instance,
the electronic structure obtained from this method is not self-consistent. Although the
tight-binding method can be done self-consistently (Ching and Callaway 1974), the
simplicity of the present energy minimisation method would be lost and the method
would become impractical if the self-consistency were to be implemented. A related
problem is that when there is charge transfer between atoms, it gives rise to an additional
Coulomb energy term contribution to AE,,, , Because of the non-self-consistency, there
is an uncertainty on the magnitude of this contribution. Therefore, it is advisable that
the method is not to be used when there is too much charge transfer (e.g., more than
a half electron) between atoms.
112 J Ihm
2.2 Quantum-chemical cluster calculations

Althogh the cluster calculation has been used widely to study molecules in quantum
chemistry, it was the Caltech group (Goddard et a1 1978) that applied this method
successfully to the study of solids, particularly to the surface and chemisorption
problems. Being originated from chemistry rather than solid state physics, this is not
one of the band structure calculational methods. Still, it. is capable of producing
excellent results for the total energy of solids. The interest of modern solid state physics
lies more in the system without three-dimensional periodicity than regular crystals and
the quantum-chemical cluster calculation is ideal for studying many such systems as
long as the system size is not too large (i.e. not more than tens of atoms). One can
mimic a non-periodic system such as a surface using a cluster of atoms and solve the
Schrodinger equation for the cluster. The calculational method used by the Caltech
group is a generalised valence bond (GVB) method. On top of that, the core is usually
replaced by a pseudopotential for heavy atoms to simplify the calculation as described
below. The GVB method can be best illustrated by considering the simplest case, HZ.
The Hartree-Fock (HF) method which is the basis for most quantum-chemical calcula-
tions will be described first for comparison with the GVB method. The presentation
published in the literature (Goddard et a1 1978) will be closely followed here.
The HF wavefunction for H2 is
a(4(1)4(2)a(l)P(2)) = 4(1)4(2)(aP - P a ) (2.8)

where a is the antisymmetriser or determinant operator, 4 is the solution of the
one-electron HF Hamiltonian, and a ( @ )is the spin-up (down) state. By symmetry, 4
is symmetric (gerade) and hence at large R (internuclear distance), 4 has the form
4=4/+4r (2.9)
where C$/ and 4r are localised on the left and right protons. Thus the spatial part of
the two-electron wavefunction is
Q H F= 4(1)4(2)= (416 + 474,) + (4141+ 66) (2.10)
= Qcov +@ion.
At large R there is an equal amount of ionic and covalent character, whereas the exact
wavefunction at R = a3 must be
= (4r4r+ 4r4r>(.P - P a ). (2.11)
As a result the H F wavefunction leads to a terrible description at large R, as indicated

in figure 1 (for H2 the error at R = 03 is 7.7 e v ! ) . The problem is that in the H F
wavefunction the motions of the two electrons in 4 are uncorrelated, resulting in too
high a probability of them being near each other. Generally, similar difficulties result
for processes involving bond formation and bond disruption. Hence, to study
chemisorption and the mechanism of oxidation, electron correlation or many-body
effects must be included.
An alternative form of H F theory is the unrestricted Hartree-Fock ( U H F ) wavefunc-
tion in which the up-spin and down-spin orbitals are allowed to be different and
non-orthogonal,
a ( C $ a ( l ) 4 b( 2 ) (1)P (2)) = C $ a + b a P - 4 b C $ a P a * (2.12)
‘)r
0.5 1.0 1.5 - n
r I ’ 1-6.0
I \ HF\ J4.0
t
I
f -0.90
8
:r
I
%
p -1.00
-+
E
I1 -1.10
-1.20
0
,
1.0 2.0
, ,
3.0 4.0
,
5.0
4.0
16.0
R [Bohr)
Figure 1. The energy of H, as a function of internuclear distance for various wavefunctions
from Goddard et al (1978).
If 4aand 4 b are different, this wavefunction describes a mixture of singlet and triplet,
whereas = 4 b leads to the closed-shell HF wavefunction of (2.8). As shown in figure
1, the UHF wavefunction is equivalent to the HF wavefunction until R = 2.6~2,= 1.4 A.
For larger R, the U H F wavefunction goes continuously to the correct limit at R =CO;
however, the bond energy is much too small and hence the U H F energy curve is not
quantitatively useful.
An alternative approach, called the valence bond method, starts with the exact
answer at R = C O , (2.11), and uses the same form for finite R,
qvB=(4/+r+ +r+/)(.~ (2.13)
where 4I and 4r are 1s orbitals of the H atom. Although exact for R = CO and
qualitatively correct for smaller R, this wavefunction is not quantitatively useful for
small R.
The generalised valence bond (GVB) wavefunction uses one orbital per electron, as
in VB, but solves self-consistently to obtain optimum orbitals as in HF. The resulting
wavefunction has the form
qGVB -(4/4r+4r4/)(43- P a )
- (2.14)
where the 4r and resemble H atom orbitals that have polarised, contracted, and
delocalised to some extent. For many-electron systems, it is convenient to transform
each G V B pair (equation (2.14)) into (orthogonal) natural orbitals
(2.15)
where S is the overlap of the GVB orbitals. With the GVB natural orbitals, equation
(2.14) becomes
qGvB - (Cg4€A1)4J2)- cu4u(1)4u(2))(aP- P a ) (2.16)
where cg and cu are suitable coefficients. Comparing with the HF wavefunction (2.8),
114 J Ihm
we see that one natural orbital, &, corresponds closely to the HF orbital, while the
other, 4u,is responsible for the electron correlation effects. Thus the G V B wavefunction
is the generalisation of the H F wavefunction in which electron correlation is incorpor-
ated and in which both the occupied and correlating orbitals are calculated self-
consistently. The form of (2.16) is easily generalised to allow more than two natural
orbitals per electron pair (in order to include additional electron correlation effects).
After determining the self-consistent GVB orbitals, a configuration interaction calcu-
lation is generally carried out
=1C I T l
TC’ (2.17)
I
in which all configurations TI formed by various occupations of the GVB orbitals are
allowed. This allows the inclusion of some additional electron correlation terms omitted
in the GVB description. The resulting wavefunction is referred to as GVB-CI.
Now we consider some simplifications of the calculation. It has long been recog-
nised that the core electrons (e.g. the He core of Li, the Ar core of Fe) do not change
significantly upon binding various atoms together and that considerable computational
simplicity could be obtained by replacing these core electrons with a pseudopotential
or effective potential. In order to illustrate some of the important considerations the
case of the Li atom will be discussed.
The ground state of Li has the configuration ( 1 ~ ) ~ ( 2and
s ) wavefunction
a((lCllsa.)(4lsP)(~2sa)). (2.18)
Since a annihilates any component of 1s contained in the 2s orbital, we can mix in
any amount of 1s character into the 2s orbital
42s= ( 4 2 s - A 4 * s ) / ( l + A 2 Y 2 (2.19)
to obtain a new wavefunction
u((~1s~)(~*~P)($2sa)) (2.20)
whose properties are all identical to those of (2.18). Because the pseudo-orbital $2s
is not orthogonal to 41s,the operator B of which 42sin (2.20) is an eigenfunction
contains new (non-local) terms. Here we will just write fi as
H =h + (2.21)
where QIs contains all the terms involving the 1s orbital. This approach is followed
in the calculations involving Si, Ga, and As below. The valence electrons are treated
explicitly so that only the Ne core
(lS)2(2S)2(2P)6 (2.22)
of Si and the Ni core
(ls)2(2s)2(2p)6(3s)2(3p)6(3d)10 (2.23)
of Ga and As are replaced by effective potentials. The approach used here is to first
carry out all electron ab initio calculations for several states of the atom and then to
find a core potential that leads to the same energies and shapes for the valence orbitals.
The procedure for As is illustrated below. We solve (all-electron) for the high-spin
states of the
[ Nil (4s)2(4P)3 (2.24)
(2.25)
(2.26)
configurations of As where [Nil is given in (2.23). The core potential is taken to be

of the form
v = vSfFs + Vdfgd + vf
+ vPfgp (2.27)
where Vsf, Vpf, Vdf, and Vf are local functions of r (distance from the As nucleus) and
g, is an angular momentum projection operator, projecting onto all components of
angular momentum 1. For s, p, and d states of the atom, the total potentials operating
on orbitals of the appropriate symmetry are
v,= V,f+ Vf
v, = VPf+ Vf (2.28)
vd = vdf+vf
while for all higher angular momenta the potential is Vf. As indicated above, the
valence 4s, 4p, and 4d orbitals are non-unique and the core character in these orbitals
was chosen so that (i) the valence-valence interactions are unchanged and (ii) the
long-range amplitude of the orbital is unchanged. With these conditions, a least-squares
procedure was used to obtain Vsf, Vpf, Vdf, and Vf so that the solution of the new
five-electron problem using V reproduces the ab initio orbitals and energies.
There are certain restrictions in this method. Since the calculation is limited to a
small number of atoms constituting the cluster, the band structure effect or any
long-range effect cannot be described accurately (or the error bound of these terms is
ha;d to estimate). However, this method seems to be superior in describing local
effects such as chemisorption problems as will be discussed in the next section.
2.3. Pseudopotential method

A method was developed by the Berkeley group (Ihm et a1 1979) to calculate the total
energy of solids within the pseudopotential formalism at about the same time as the
tight-binding energy minimisation method was developed by Chadi. This method,
when combined with a first-principles pseudopotential discussed later, turned out to
give structural energies of unprecedented accuracy, enabling us to study not only the
cohesive energy but such 'exotic' properties as the pressure-induced phase transitions.
The momentum-space formalism for the pseudopotential total energy published in the
literature (Ihm et a1 1979) will be presented below in detail with some modifications.
Following the conventional density functional formalism (Kohn and Sham 1965)
in a pseudopotential framework (Phillips and Kleinman 1959), the total crystal energy
(defined as the total energy difeerence between the solid and isolated cores) is given by
E,,, = T + V+
I
where the total kinetic energy, T, is
Exc(r ) d3r (2.29)
T = z + ? ( r ) ( - V 2 ) G i ( r ) d3r (2.30)
I
and the electrostatic potential energy, V, is

116 J Zhm
Rydberg atomic units are used here. $i( r ) is the (pseudo) wavefunction of the valence
electron where the index i denotes both the wavevector ki and the band index n and
runs over all occupied valence states. p ( r ) = E i $ ? ( r ) $ i ( r ) is the (pseudo) valence
electron density, 2 is the valence of the ion and the R, is the lattice vector to the
ionic site. Although we confine ourselves here to one kind of ion for notational
convenience, generalisation to many kinds of ions is straightforward. The first term
in (2.3 1) is the core-valence* interactton energy for angular-momentum-dependent
pseudopotentials (El U p s , lr() P I ,where PI is the projection operator on angular momen-
tum I ) . The second term is the (valence) electron-electron Coulomb energy and the
last term is the lattice (ion-ion) energy. The ion-ion interactions can be replaced
rigorously by point ion interactions as long as the ions are spherically symmetric and
non-overlapping. The last term in (2.29) is the density functional exchange-correlation
contribution to the total energy.
The one-electron Schrodinger equation derived variationally from (2.29) is
( -V2+ 1 Ups,,(r - R p ) k l+
P. I
2p0d3r‘+pxc(r ) ) m i ( r )= Ei$i(r)
Ir-r’I
(2.32)
where
Pxc(r) = a E x c ( r ) / J P ( r ) . (2.33)
The general form for E,, is not known. Employing the X a method (Slater and Johnson
1972), we obtain
3
p x c ( r )= - a - (3.ir2)”3(p(r))1’3 (2.34)
.ir
and
( p ( r ) ) 4 ’ 3= i p ( r ) p x c (r ) . (2.35)
For explicitness, (2.34) and (2.35) will be used to calculate the exchange correlation
in the present derivation. It should be straightforward, however, to modify the
expression in compliance with more general approximations for the exchange corre-
lation.
Thus far, we have not made use of the periodicity of the system. Now every quantity
will be expanded in terms of plane waves. With a smooth pseudopotential, the
plane-wave expansion converges rapidly. If we generalise it to the mixed basis formal-
ism, a few integrals (of the Gaussian type, for example) are to be added. The
momentum-space representations of the wavefunction, the charge density, the inter-
electronic Coulomb potential and the exchange-correlation potential are
G c ( r ) = I $ ( k , + G )ei(kl+c)‘r (2.36)
C
p(r)=Cp(~)e-’~‘~ (2.37)
C
V c o u l ( r=) 2p0d3r’= 1 VC,,,(C ) eic’r (2.38)

lr-r’l C
and
(2.39)
where the G are reciprocal lattice vectors. The exchange-correlation Fourier component
Pxc(G) is
j 2 1/3 1/3 G)
Pxc(G)=-93T) P ( . (2.40)
The Fourier component of the inter-electronic Coulomb potential, Vco,,( G), is

VCOdG) = 8TP(G)/G2 (2.41)
which is obtained by solving the Poisson equation
= --8Tp(r).
V2Vcou,(r) (2.42)
Note that Vco,,(G ) and pxc(G) have the dimension of energy. Let Cl denote the total
volume of the system and N denote the total number of atoms and Clat = Cl/ N. The
terms entering (2.30) can be reduced as
(2.44)
(2.45)
= $*(ki+G)$(ki+G') exp[-i(ki+G).(r+R,)]
i,w, 1, G,G'
x UPJ r ) A exp[i( k, + G') ( r + R,)] d3r

0
=Cl c
i,l,G,G'
+*(ki + G ) 4 ( k i + G') 1exp[i(G'- G).R,]/N
P
xL
Clat
[ exp[-i(ki+ G)*r]Ups,l(r)@l
exp[i(k,+ G').r] d3r
(2.46)
where S ( G'- G ) is the structure factor and
Ups,l,k, + G , k ,+ G' -
=- I exp[-i(ki+G)*r]Up,,,(r)filexp[i(ki+G).r]d3r
fiat
I
- 1 (21+1)4~ U,,,l(r)jl(~ki+G~r)j~(~ki+G'~r)r2dr.P,(cosy).
--
Rat
(2.47)
118 J Ihm
In the last line of (2.47) j , and PI are spherical Bessel functions and Legendre
polynomials respectively, with
( k i + G ) * ( k ,+ G ’ )
cos y = (2.48)
I k, + GI * I k, + G ’ J’
It is by virtue of the plane-wave expansion that we have such a simple analytic
expression for the non-local pseudopotentials.
In the local pseudopotential approximation, we have
C
hcL
+:(r)Ups(r-RP” d3r=C j U p s ( r - 4 ) d r ) d3r
= fl c S ( G )Up,(G)P(G).
G
(2.49)
It will be later shown that the multiple sum in (2.46) can be avoided in the final
expression for the total energy (equation (2.64)). For practical purposes, it is very
convenient to decompose the general non-local pseudopotentials into a purely local
part and non-local parts, namely,
c U p s , / ( r ) k Ups(r)+c( Ups,l(r) Ups(r))4
1
=
I
-
= Updr)-tC U ; s , l ( r ) 9 (2.50)
where r ) = Ups,/(
r ) - Up,(r ) , such that the purely local pseudopotential Ups(r )
-
takes care of the long-range interaction ( U,,( r ) - 2 z / r for large r ) thereby making
the non-local parts, U ; s , l ( r )of, short range. This procedure is always possible because
-
UPs,/(r) -2Z/r for any 1 (i.e. angular momentum dependence disappears for large
r ) . We still have infinite degrees of freedom for the small-r behaviour of Up,(r ) . Up,(r)
can be chosen such that both Ups(r)and U ; , , / ( r )be as smooth as possible. In the
local approximation, U;,,,( r ) = 0 identically. Now the Fourier transforms of U;,,/(r)
do not have a singularity at G = 0. Combining (2.46), (2.49) and (2.50),the core-valence
interaction energy becomes
4/41
J c L ~ ( r ) U p s , / ( r - R , ) ~ l + ,d3r
(r)
S ( G )Ups(G ) P( G )
= fl( G
+ C
I,I,G,G
cL*(k, + G)$’(ki
)
G ’ ) S ( G ’ -G )U ~ S , I , ~ , + G , ~ , + G (2.51)
Combining (2.36), (2.38), (2.39) and (2.50), the one-electron Schrodinger equation
(equation (2.32)) becomes
C [(kl + G ’ ) * ~ G+GVLGlcL(kl + G’>= ~ , $ ( k+, G )
G
(2.52)
where
VLG = V , , , , ( G f - G ) + ~ u , , ( G ‘ - G ) 4 S ( G ‘ - G ) U p , ( G ’ - G ) + C Uks,i,k,+c,k,+c
I
(2.53)
To simplify (2.54), we multiply on the left of (2.32) by +?(r), integrate over r and
sum over i, and substitute the result into (2.29),
(2.55)
In practice, some mathematical manipulations are necessary to calculate E,,, from
(2.55) because Vcoul(0),Up,(0) and
P#V
are infinite quantities. First, we solve (2.52) with Vcoul(0)and Ups(0)set equal to zero.
This corresponds to a constant shift of the potential. This arbitrary shift is compensated
for by the following procedure. For small 1 G I,
87
Ups(G ) = - -
T
z
+ cy, +(higher terms in G ) (2.56)
%,G2
that is,
C+O
), ’
U p , ( G ) + y --
GG
-
fiat
I( U p s ( r ) + y ) d3r. (2.57)
Note that (2.57) is no more than a formal expansion because G takes only discrete
values in periodic systems. Also notice that we do‘not have a term of order 1 / G in
the expansion in (2.56). If Ups(G) should have a term like l / G in the expansion near
G = 0 , then cy, defined by (2.57) would diverge. This means the total energy shown
later in (2.64) would diverge. One example of this behaviour is the Simons-Bloch
pseudopotential (Simons and Bloch 1973),
(2.58)
where the BI are constants depending on 1. cy1 diverges with this potential (unless an
artificial cut-off rmaxfor the integral is introduced) due to the long-range character of
-
the l / r 2 potential added to - 2 Z / r . (In fact, cy1 1 / G rather than a constant.) Even
120 J Ihm
with the ions to neutralise the whole system, the total average potential experienced
by the electrons is infinite in this pseudopotential. The Simons-Bloch pseudopotential
without cut-off rmax is inadequate for application to extended systems even though it
gives good results for localised systems like atoms or molecules. Physically, a 1 is a
measure of the repulsiveness of the pseudopotential which keeps the valence electrons
from collapsing into the core.
Another quantity relevant to the total energy is p, defined by
z
p ( G ) =-+pG2+(higher terms in G) (2.59)
Rat
for small G. From (2.41) we have
(2.60)
p is a rather complicated integral. Fortunately, p does not appear in the final

expression for the total energy because both the electron-electron interaction and the
ion-electron interaction contain p and they cancel out as will be shown in (2.63).
The lattice (ion-ion) energy per atom is usually expressed as
(2.61)
that is,
(2.62)
where the prime means that R , = 0 is excluded in the summation. Combining equations
(2.41), (2.56), (2.59) and (2.62), the energy per atom coming from three divergent terms
(with S(0) = 1 ) is:
In summary, we solve (2.52) with Vcoul(0)and U,,(O) set equal to zero. Then we
add ( a , Z + YEwald) to the total energy. The final expression for the total energy per
atom is
(2.64)
The first term on the right-hand side is the sum of the electron eigenvalues of the
occupied states, the second and the third terms correspond to the correction for
overcounting of the electron-electron interaction. The fourth term is the correction
coming from the ‘pseudo’ nature of the potential and the last term is the ion-ion
Coulomb energy. The non-local pseudopotential contribution (equation (2.46)) does
not explicitly appear here; it enters (2.64) indirectly through band structure
eigenvalues,
When first-principle pseudopotentials are employed in the total energy calculation
as described in the next section, this method seems to give more accurate results for
the structural properties of solids than any other method available. The present method
is widely used to obtain very detailed information on the energy and the valence charge
density. There are, however, two problems associated with this method. Firstly, since
this method solves the Schrodinger equation very accurately, the computing time is
enormous for a large unit cell and the calculation is impractical at the present if the
unit cell size is more than a few tens of atoms. This can be a serious restriction to the
study of a non-periodic system because the pseudopotential method uses a large
supercell geometry to mimic such a non-periodic system. Secondly, the local density
functional ( LDF) formalism is almost always used for the exchange-correlation potential
to reduce the complexity in the self-consistent pseudopotential calculation. The density
functional formalism is designed to describe the ground-state properties of solids. It
turns out that the excited states, especially the size of the gap, are poorly represented
with the LDF approximation as will be discussed in the final section. Therefore, the
first-principle pseudopotential methods with the LDF approximation should be used
with caution in studying conduction bands of solids.
2.4. Other theoretical tools

There are other theoretical methods of calculating the total energy which are not
discussed in detail here. The LMTO method (Andersen 1975) was mentioned in the
introduction and the same muffin-tin approximation was used with the KKR method
(Andersen and Jepsen 1977). The LAPW method (Posternak et a1 1980, Feibelman and
Hamann 1983) is actively used currently. A completely different approach, an exact
stochastic simulation of the Schrodinger equation for charged bosons and fermions,
was used to obtain the correlation energies to locate the transition to the lowest energy
phase at zero temperature (Ceperley and Alder 1980). Very impressive results were
obtained for the defect problem using the Green function method (Baraff and Schluter
1983). A first-principles LCAO method was developed and applied successfully to
diamond (Chelikowsky and Louie 1984).
Cross checks between different methods have been made as well. Adequacy of the
frozen core approximation in the pseudopotential method compared with the all-
electron calculation was discussed by von Barth and Gelatt (1980). Excellent agreement
between hard-core pseudopotential results and all-electron results was reported for
the atomic Si and the Si2 molecule within the LDF formalism (Harris and Jones 1978).
The same conclusion was reached for the bulk Si crystal by comparing pseudopotential
calculations and the all-electron LAPW (Hamann 1979) or the all-electron LMTO
122 J Zhm
(McMahan et a1 1981) calculations. Also, for the prototype case of the GaAs(ll0)
surface, the tight-binding (Chadi 1978), the quantum-chemical cluster (Barton et a1
1979), and the pseudopotential (Ihm and Joannopoulos 1981b, Pandey 1982) methods
all give essentially identical results for the energy minimum geometry.
3. Applications
3.1. Structurally related properties of bulk materials

As mentioned in the introduction, there have been tremendous efforts to understand
structurally related properties of bulk solids from the total energy point of view. These
properties include the cohesive energy, the equilibirum lattice constant, the bulk
modulus, the shear modulus, the Griineisen parameters, the ground-state structure,
and the critical pressure and the latent heat involved in the pressure-induced phase
transition between different structures, to mention a few. These theoretical efforts were
quite successful for most molecular crystals, ionic crystals, and simple metals but they
were less successful for transition metals and particularly for covalent semiconductors.
Recent development of the theoretical approach, however, has circumvented difficulties
in dealing with d electrons in the transition metals or directional bonds in important
semiconductors and now it is possible to study most elements in the periodic table
quite routinely using essentially the same computer program. In this section, we will
concentrate on the success of the new techniques. Among the three theoretical tech-
niques described in detail in the previous section, the tight-binding and the quantum-
chemical cluster method usually adjust parameters to fit bulk properties and use these
parameters to study more sophisticated systems such as a surface or a chemisorbed
system. Therefore, the pseudopotential method which does not use these parameters
is more appropriate to study bulk properties.
The momentum space formalism for the total energy with the pseudopotential
Hamiltonian derived in § 2 was first applied to bulk Si (Ihm and Cohen 1979, 1980).
The cohesive energy and the equilibrium lattice constant derived from this self-
consistent pseudopotential calculation agreed with experiment very well. Most notably,
among the BCC, FCC, HCP, white-tin and diamond structures, it was shown that the
diamond structure is the lowest in energy and the energy difference between different
structures is only a few tenths of an eV per Si atom. Although this result was very
impressive at that time, since the ionic pseudopotential used in the calculation (Schliiter
et a1 1975) was fitted to energy levels of the atom or the ion and did not necessarily
reproduce the all-electron valence wavefunctions correctly, the accuracy of the result
was still limited. For example, the calculated bulk modulus which is the second
derivative of the energy with respect to the volume change deviated from experiment
by 70%.
In the meantime, there was an important breakthrough in the development of the
pseudopotential theory. The norm-conserving pseudopotential was invented in 1979
(Hamann er a1 1979) which is capable of reproducing the all-electron wavefunctions
in the valence region exactly as well as the energy levels of the atom. As the amplitude
of the electron wavefunction in the valence region is correct, the net charge in the core
region must be correct as well when the total charge of an electron is normalised to
one, hence the name of the norm-conserving pseudopotential. Although the idea of
adjusting the pseudopotential to fit all-electron wavefunctions is not entirely unique
(Starkloff and Joannopoulos 1977, Zunger and Cohen 1979b, Kerker 1980), the norm-
conserving pseudopotential by Hamann-Schliiter-Chiang ( HSC) is most systematic in
derivation and most useful for practical applications and so is most widely used now.
The HSC pseudopotential is a truly first-principles pseudopotential in that one essentially
needs only the atomic number 2 to generate the pseudopotential.
By means of this powerful pseudopotential and the aforementioned momentum
space formalism, new results on the structural properties of solids have been produced
with unprecedented accuracy (Yin and Cohen 1980, 1982). They are summarised in
figure 2 and tables 3 and 4. In table 3, the comparison between the calculated and
measured properties of silicon in the diamond structure has been made. The agreement
is excellent for all the studied properties including the bulk modulus which is usually
hard to reproduce accurately. A more startling result is contained in figure 2 and table
4. The crystal energy for different structures is calculated and it is concluded from a
purely theoretical calculation that the silicon undergoes a pressure-induced phase
transition from the diamond to the white @)-tin structure at around 100 kbar, with a
Volume
Figure 2. The diamond, hexagonal diamond, p-tin, sc, HCP, B C C , and FCC structural
energies (in units of Ryd/atom) as a function of the atomic volume (normalised to the
measured free volume) for Si from Yin and Cohen (1982). The broken line is the common
tangent of the energy curves for the diamond and the p-tin structures, indicating the path
of the pressure-induced phase transition.
Table 3. Comparison of calculated and measured static properties of silicon, from Yin
and Cohen (1980).
Lattice Crystal Cohesive Bulk

constant energy energy mod u I u s
(A) (Ryd) (ev) (Mbar)
Calculation 5.451 -7.909 4.67 0.98

Experiment 5.429 -7.919 4.63 0.99
Difference (YO) 0.4% -0.1% 0.9% -1 %
124 J Ihm
Table 4. Comparison of the calculated and measureda volumes ( V t * pof ) the diamond and
@-tin phases at the transition pressure, their ratios ( Vf/ V f ) , and the transition pressure
from Yin and Cohen (1980). The volumes are normalised to the measured zero-pressure
volume.
v: VP VPl v:' P, (kbar)
Calculation 0.928 0.718 0.774 99

Experimenta 0.918 0.710 0.773 125
Deviation 1.1% 1.1% 0.1% -20%
a Jamieson (1963), Weinstein and Piermarini (1975).
volume reduction by 23%, in remarkable agreement with experiment ! When superstress

effects and the temperature dependence are considered, the agreement with experiment
in the critical pressure can be much closer in table 4. Such an amazing success indicates
the power of the norm-conserving potential and the momentum-space formalism and
it is also a triumph of the LDF formalism in describing the ground-state properties of
solids. More recently, superconductivity of the high-pressure phase Si has been
predicted with the first-principles pseudopotential calculation and confirmed by experi-
ment (Chang et a1 1985). Such a success, of course, is not limited to silicon. It
stimulated calculations for many different materials by different people and results
have been always impressive. Total energy calculations using this first-principles
pseudopotential method have been reported for other semiconductors such as Ge (Yin
and Cohen 1981a), Sn (Ihm and Cohen 1981), GaAs (Kunc and Martin 1981), GaAs
and AlAs (Ihm and Joannopoulos 1981a), Gap, Si, and C (Bachelet et a1 1981), Alp,
AlAs, Gap, and GaAs (Froyen and Cohen 1982); for simple metals such as A1 (Lam
and Cohen 1981), semimetals such as As (Needs et a1 1986), transition metals such
as M O and W (Zunger and Cohen 1979a), Zr, Nb, and MO (Ho et a1 1982); and for
insulators such as Se (Vanderbilt and Joannopoulos 1980) and C (Yin and Cohen
1981~).The list is still expanding quite steadily. There has also been a systematic
study of ionicity and the structural stability of solids using the total energy method
(Chelikowsky and Burdett 1986).
3.2. Phonons
The total energy method is useful to calculate the phonon frequency of solids as
demonstrated already in 0 2 when the tight-binding method was described. By compar-
ing the energy of an undistorted crystal and that of a distorted crystal corresponding
to a particular phonon mode, the phonon frequency is directly obtained. The Born-
Oppenheimer approximation (1927) is commonly used where the electronic system is
assumed to be in the ground state with respect to the instantaneous nuclear positions.
Wendel and Martin (1978) used a ionic model pseudopotential by Appelbaum and
Hamann (1973) to study To(r) and TA(X)phonon modes and the shear constant
Cll-C,*. The results were satisfactory except for the TA(X)mode which is more than
40% greater than the experimental value. The unusually low frequency experimentally
observed for the TA(X)mode has been well known and has led to the adiabatic shell
and bond charge models previously (Weber 1974, 1977). Using the tight-binding
method, Chadi (1978) also obtained good results for the elastic constants and the
phonon frequencies, but again the calculated TA(X)modes were too high as shown in
table 2. A different approach based on the dielectric screening was used to calculate
phonon dispersion curves of Si (Van Camp et a1 1979). The agreement for the TA(X)
mode was gratifying (20% too high). With the advent of the momentum space
formalism and the norm-conserving pseudopotential, the first-principles calculation of
the phonon frequency at special points has become straightforward. Since the first-
principles pseudopotential method solves the Schrodinger equation of the system
almost exactly in the valence region, no model building or ad hoc correction is required
to obtain the phonon frequency. The results from the first-principles pseudopotential
calculation for silicon (Yin and Cohen 1980) are presented in table 5. The mode
Griineisen parameter y is defined to be -d(lnf)/d(ln V). The agreement between the
calculation and experiment is excellent, without any problems related to the soft TA(X)
mode. The same method has also been applied to transition metals (Ho er a1 1982).
Here, phonon anomalies in MO, Nb, and BCC Zr have been studied in relation to
superconductivity. The results are in remarkable agreement with experiment, exhibiting
an anomalously soft phonon mode at the L ( f ,i, :) point. The Xerox group (Kunc and
Martin 1981) also reported good results for the phonon modes of GaAs using a
self-consistent pseudopotential method.
Table 5. Comparison of the calculated phonon frequencies (fCaJ and the mode-Griineisen
parameters (-yCa,J with experiment" (f,,,,, -ye.,,) from Yin and Cohen (1980). The devi-
ations from f,,,, are also presented. All frequencies are in units of 10" Hz.
LTo(r) TAW TdX) LOA(X)
fcak 15.16 4.45 13.48 12.16

fexpta 15.53 4.49 13.90 13.32
Percentage difference -2.4% -0.9% -3.0% - 1.3%
Ycalc 0.92 -1.50 1.34 0.92
Yexprb 0.98 -1.4 1.5 0.9
Measured at 296 K (Dolling 1963, Nilson and Nelin 1972).

yexptof LOA(X) is estimated from the thermal expansion coefficient (Weinstein and
Piermarini 1975).
Another way of obtaining the phonon frequency exists in the literature. By calculat-
ing the Hellmann-Feynman (Hellmann 1937, Feynman 1939) force on the nucleus,
one can get the phonon frequency directly. The Hellmann-Feynman theorem has been
rederived in case of the statistical exchange method (Slater 1972) and the pseudopoten-
tial method (Ihm et al 1979). Good results with essentially the same accuracy as the
total energy method have been obtained by means of the Hellmann-Feynman force
method (Ihm et al 1981). A bonus is that the anharmonicity can be studied much
more easily with the force method. The theory of the quantum mechanical force and
stress has been further developed by Nielsen and Martin (1985). The surface phonon,
on the other hand, has been studied for the Si (100) surface via empirical tight-binding
methods by Allan and Mele (1984).
3.3 Surface
The GaAs (1 10) surface, having been examined by various experimental techniques,
is perhaps the best understood semiconductor surface at present. Therefore it is a
126 J Ihm
convenient reference surface to test the validity of a theoretical technique. The tight-
binding, the quantum-chemical cluster, and the pseudopotential method all give similar
results for the stable structure of this surface in agreement with experiment. The
tight-binding calculation for GaAs and ZnSe surfaces (Chadi 1978) will be examined
here first.
The two-dimensional unit cell for the ( 1 10) surface of GaAs or ZnSe is shown in
figure 3. The changes in atomic coordinates (for the first three surface layers) from
unrelaxed, bulk terminated positions calculated for GaAs and ZnSe as well as several
LEED determined structures for GaAs are shown in table 6. The angle 0 formed by
nearest-neighbour cation-anion bonds at the surface with the bulk terminated plane
and AE,,, (in eV per surface atom) are also given in table 6. For GaAs the value of
-0.51 eV per surface atom for AE,,, is obtained. The absolute values of the normal
displacements of surface and subsurface atoms in GaAs are in good agreement with
the LEED results. For ZnSe LEED studies (Duke et a1 1977) indicate a surface structure
similar in many respects to that of GaAs but no definite structure for ZnSe has yet
Figure 3. The two-dimensional unit cell for the GaAs (110). The open circles denote
cations (Ga) and the filled circles denote anions (As). See table 6 for atomic displacements
from unrelaxed positions shown above.
Table 6. Atomic displacements in A from unrelaxed positions for the (110) surfaces of
GaAs and ZnSe reproduced from Chadi (1978). The superscripts a and c refer to anion
and cation respectively. The subscript 1 refers to the surface layer, 2 to the subsurface,
etc. The normal to the surface is in the z direction. AE,,, is in eV per surface atom. 6 is
the bond rotation angle at the surface.
GaAs ZnSe
Reference a Reference b Reference c Reference d Reference a

~~ ~~
0 27.3" 26.4" 34.8" 27.1" 25.6"

AEm -0.51 -0.39 -0.42 -0.46 -0.30
Ay;(AyY) 0.19 (0.35) - 0.34 (0.60) 0.32 (0.48) 0.18 (0.37)
Az;(Azf) 0.19 (-0.46) 0.20 (-0.50) 0.20 (-0.60) 0.10 (-0.55) 0.04 (-0.55)
Az;(AzS) -0.06 (0.07) -0.05 (0.05) - - -0.05 (0.12)
Az!(AzS) 0.02 (-0.04) 0.025 (-0.025) - - 0.04 (-0.04)
a Present work.
Mark et a1 (1978).
Lubinsky et al (1976).
Tong ef al (1978).
been determined. This calculation shows the main difference in the surface relaxations
of GaAs and ZnSe to be in the normal displacements of the surface anions; the Se
atoms do not move appreciably from their bulk terminated positions as compared to
As atoms. The bond length variations are calculated to be between -2.9 to 1.7% in
GaAs and -7 to 2.7% in ZnSe. The energy lowering arising from subsurface relaxations
was found to be very small (s0.02eV), typical of room temperature bulk phonon
frequencies. This suggests that low-temperature LEED measurements would be required
if subsurface relaxations are to be accurately determined.
The quantum-chemical cluster calculations using the generalised valence bond
method described in 0 2 will be examined next. In modelling this surface a finite
portion of the surface is sliced out of the solid and every Ga-As bond so cut is replaced
by a bond to an H atom. In this way the proper hybridisation is retained at each atom
and there are no spurious radical centres. Thus the minimal complex for studying
reconstruction of the GaAs surface consists of Ga,AslH,, where the Ga and As are
both surface atoms.
The geometric variations were carried out as follows.
(i) All second-layer atoms were fixed.
(ii) The surface Ga and the surface As were moved independently but with the
constraints that the bond length to the (virtual) second-layer atoms were fixed and
that there were no lateral displacements of this bond (that is, the surface atom remains
in the y z plane defined by this AsGa interlayer bond and the normal to the surface).
(iii) The virtual surface atoms being represented by H atoms are moved by the
same amount as the real surface atoms of item (ii).
With the above constraints there are two independent parameters remaining, one
for the Ga and one for the As. These will be written as SzGaand S z A S ,referring to the
projection of the diplacement upon the normal to the surface (positive Sz is away from
the surface); of course, S z Z O implies that there is also motion in the y direction.
Since one experimental probe of the surface atoms is the chemical shift in the 3d levels
of the Ga and As, we have included all electrons in the current calculations. The
optimum displacement is calculated as
SZG, = -0.44A SyG, = 0.48 A
8zAs = t-0.22 A SyA, = 0.37 A
with a total energy drop of 1.1 eV resulting from the reconstruction. The experimental
results are often cast in terms of an angle 8 corresponding to the angle between the
unreconstructed surface (110) plane and the line connecting a surface Ga to a surface
As, measured in the (11) plane perpendicular to the surface plane. The calculated 8
is 25.6", which is in good agreement with experimental results as well as with the
tight-binding calculations.
The close correspondence between calculated and experimental geometries indi-
cates that the surface reconstruction is dominated by rehybridisation or atomic valence
effects of orbitals localised at the surface. This validates the use of finite clusters for
calculation of reconstruction geometries. Indeed, the bond angles subtended at the
As are 90", 90", and 108" for an average of 96", close to the value of 92" found for free
ASH,. The bond angles at the Ga are 125", 125", and 108", leading to an average of
119", close to the value of 120" expected for free GaH,. Thus the surface atoms come
close to the valence angles expected for isolated trivalent molecules. The orbital
character and charge on the surface Ga and As atoms are different from that of the
bulk atoms; the valence distribution on these surface atoms is quite asymmetric. This
128 J Ihm
suggests that in dynamic LEED calculations, one should be cautious about using
sphericalIy symmetric, muffin-tin potentials for these surface atoms. However, the
good agreement between the calculated reconstruction and the experimental value
based on dynamic LEED calculations suggests that the usual approach does not lead
to errors affecting the geometry.
The core ionisation potentials shift as the surface distorts from the purely tetrahedral
position to the optimum reconstructed position (see table 7). The 3d level on As shifts
to higher binding energy by 0.24eV and the 3d level on Ga shifts to lower binding
energy by 0.19 eV. These shifts would be too small to resolve with current photoemission
experiments. The empty orbital localised above the Ga centre becomes more p-like
during reconstruction and moves up in energy (less bound) by 1.2 eV. Recent band
structure calculations (Chadi 1979b) have also shown that the relaxation moves empty
Ga-derived surface states up (and out of the band gap).
Table 7. Theoretical results for GaAs surface clusters from Barton et al (1979). Results
for oxidised surfaces in the following section are presented together.
Orbital energies (Hartree)a

Total
Twist angle energy" G a 3d As 3d Empty As lone
6 (ded (Hartrees) (average) (average) G a 4p pair
H,GaAs'
0.0 -4159.3761 -1.2441 -2.0998 -0.0244 -0.2936
14.3 -41 59.4079 -1.2382 -2.1063 0.0030 -0.3137
19.3 -4159.4135 - 1.2364 -2.1077 0.0109 -0.3192
27.0 -4159.4162 - 1.2367 -2.1086 0.0204 -0.3253
H,GaAsOb.'
0.0 -4234.1 600 -1.271 -2.196 -
14.3 -4234.1834 -1.2673 -2.2009 -0.0185
19.3 -4234.1850 -1.2675 -2.2018 -
27.0 -4234.1795 - 1.2677 -2.2005 -0.0032
a 1 Hartree=27.2116eV.
Energy of 0 atom in the same basis is O(3P) = 74.8006 Hartrees, and the As-0 bond length is 1.72 A.
Hartree-Fock wavefunctions.
The first-principles pseudopotential calculations (Ihm and Joannopoulos 1981b,

Pandey 1982) give essentially the same energy minimum structure and will not be
discussed further. Note that the energy reduction by relaxation is about 0.65 eV per
unit cell compared with 1.02 eV in the tight-binding method and 1.1 eV in the quantum-
chemical cluster calculation.
There are also interesting calculations for the Si (111) and Si (100) surfaces. The
relaxed Si (111) surface was first studied by Redondo et a1 (1976) using the quantum-
chemical cluster method and by Ihm and Cohen (1979, 1980) using the self-consistent
pseudopotential method. The reconstructed Si ( 111) (2 x 1) surface was studied by
Chadi (1978) using the tight-binding method within Haneman's buckling model
(Haneman 1961). After Pandey's 7-bonded chain model (Pandey 1981) was proposed
for the surface, the first-principles pseudopotential calculation was performed for this
reconstruction model, confirming the stability of the 7-bonded chain (Northrup and
Cohen 1982). The reconstructed Si (1 11) (7 x 7) surface is too difficult to use the total
energy method directly and will not be discussed here.
The total energy calculation for the Si (100) (2 x 1) surface was first done by means
of the tight-binding method (Chadi 1979a) and an asymmetric dimer model was
proposed for the (2 x 1) reconstruction based on the energy minimisation. Then the
self-consistent pseudopotential calculation was done by Ihm et a1 ( 1980a), providing
more detailed information on the asymmetric dimer model. A more accurate calculation
using the first-principles pseudopotential method followed (Yin and Cohen 198l b )
which confirmed the asymmetric dimer model with a modification in the dimer bond
length. The quantum-chemical cluster calculation, on the other hand, indicates that
the symmetric dimer is more stable (Redondo and Goddard 1982) than the asymmetric
dimer. A recent scanning tunnelling microscope picture (Tromp et a1 1985) shows
both symmetric and asymmetric dimers on the surface. It is not clear at this point
which of the two is the more stable configuration for the Si (100) (2 x 1) surface.
Another important problem associated with the Si (100) surface is a higher-order
reconstruction occasionally observed by means of LEED (Lander and Morrison 1962,
Poppendieck et a1 1978), the He atom diffraction (Cardillo and Becker 1978, 1980),
and the reflected high-energy electron diffraction experiment (Ichikawa and Ino 1979).
Based on the tight-binding total energy calculations combined with renormalisation
group techniques in statistical mechanics, it has been proposed that an order-disorder
phase transition from the c ( 4 x 2) structure to a disordered (2 x 1) structure exists at
around 250 K for both Si (100) and Ge (100) surfaces (Ihm et a1 1983b, c, 1985). This
occurs because an asymmetric dimer has two possible orientations in space and the
orientation of dimers on the surface is ordered with respect to each other at low
temperatures and becomes disordered at higher temperatures. Since this is a new area
of the application of the total energy method, more details (Ihm et al 1983a, b, c) will
be presented below.
All total energy calculations so far have been resiricted to zero-temperature
geometries because of the lack of a tractable scheme to incorporate the entropy
contribution at finite temperatures. However, if one were able to combine the energy
minimisation approach with the renormalisation-group technique which has had great
success in phase transition studies for two-dimensional systems, the behaviour of
semiconductor surface at finite temperatures could be predicted quantitatively from
first principles. With the use of a series of approximations described below, such a
scheme has been successfully developed and applied to the Si (100) surface, resolving
important questions regarding the structure of the Si (100) surface. For example, it is
shown that the (2 x 1) structure is not the ground state of the Si (100) surface and a
higher-order reconstruction should occur on the surface. Disappearance of these
higher-order spots near room temperature due to an order-disorder transition is
predicted. A phase diagram for the Si (100) surface is also obtained which can be
used simultaneously for systems belonging to the same universality class such as Ge
(100) and diamond (100).
Basic approximations involved are summarised as follows. (i) Use of the asymmetric
dimer model. Pairs of atoms at the surface relax by dimerising into one of the two
possible asymmetric configurations (figure 4 , inset). Different reconstructions result
from different arrangements of asymmetric dimers which are the building blocks of
the surface in this model. (ii) Use of the tight-binding energy minimisation method
by Chadi. The excited s orbital has been included in the basis set to improve the
accuracy of the calculation (Wang and Joannopoulos 1980, Vogl et a1 1983). (iii)
Mapping onto the spin-; Ising Hamiltonian. The two possible orientations of the
asymmetric dimer are represented by the two possible states of a spin, and the energy
130 Ihm
(2x1)
~14x2)
--.c c --c
c 4 c
c +
I
+ c +
Figure4. Reconstruction geometries of the Si (100) surface for ( a ) ‘(2x 1)’ and ( b )‘ c ( 2x 2)’
families reproduced from Ihm et al (1983a, b, c). Coupling constants for each family are
indicated. Side views of the asymmetric dimers are shown in the inset.
differences between different reconstructions are translated into a set of interaction

energies in an effective spin Hamiltonian. (iv) Use of the standard position-space
renormalisation group method (van Leeuwen 1975) to determine the phase diagrams.
The first and second approximations involve rather controversial issues. For
example, the atomic positions on the Si (100) surface are not precisely known.
Nevertheless, there is ample experimental support for the asymmetric dimer model
and it is the best choice available at present. Next, it is difficult to estimate the precision
of the energy differences between different reconstructions obtained from the tight-
binding calculations. It is known that the accuracy of the pseudopotential energy
calculation is better than 10 meV, but pseudopotential calculations for reconstructed
geometries greater than ( 2 x 1) or c(2 x 2 ) are intractable at the present. After extensive
tests and crosschecks with pseudopotential calculations for simpler geometries [ ( 2 x 1)
and c(2 x 2 ) structures], we find that the accuracy of the tight-binding energy calcula-
tions can be made comparable to that of the pseudopotential method when excited s
orbitals are properly included in the basis set. With this extra orbital, we reproduce
the energy difference between the symmetric and the asymmetric dimer (0.2 eV) to
within 10 meV and obtain the ( 2 x 1) structure lower in energy than the c(2 x 2 ) structure,
in agreement with both pseudopotential calculations and experimental observations.

(Without this improvement over Chadi’s original scheme, (2 x 1) would have a higher
energy than c(2 x 2).)
The third approximation is the most natural and simplest to make. A remaining
question relates to the number of neighbour interactions that must be included in the
effective spin Hamiltonian. In practice, this number is uniquely determined by the
number of different reconstruction geometries that need to be compared. Experi-
mentally, higher than quarter-order spots are never observed. Thus, the number of
possible reconstructions deserving analysis is small (4-5) and can be accommodated
by including up to the next-nearest-neighbour interactions in each dimension.
Regarding the fourth approximation, this class of spin systems has been studied
extensively and the description of the phase transitions using the renormalisation group
approach can be made as accurate as required.
Let us now turn to the specific problem at hand. Various diffraction experiments
at room temperature indicate that ( 2 x 1) is the basic reconstruction unit for the Si
(100) surface. Nevertheless, higher-order spots (up to quarter order) have been seen
and, occasionally, diffuse background as well as streaking is also observed. A family
of simple reconstruction geometries with a (2 x 1) backbone is shown in figure 4(a).
The spins indicate the two possible asymmetric dimer orientations. This set of structures
will henceforth be called the ‘2 x 1’ family.
The effective ‘spin’ Hamiltonian for this family may be written as
where all interactions up to twice the surface atom spacing are included and illustrated
in figure 4(a). The coupling U and F contribute equally to the ground-state energies
considered here and are taken to be zero. The remaining coupling constants are
extracted from the total energy differences of the individual structures. The total
energies at T = 0 are calculated with use of the aforementioned tight-binding approach
for infinite slabs of Si with twelve-layer thickness. A position-space renormalisation
group calculation is then performed, using a finite cluster of four cells, each containing
five sites. The corresponding flows are in the parameter space of V, H, 0,and F only.
There is another important family of reconstruction geometries based on a c(2 x 2)
backbone. This family is shown schematically in figure 4(b). The interaction para-
meters for the spin Hamiltonian of this ‘c(2 x 2)’ family are also illustrated in this
figure. These coupling constants are again extracted from total-energy calculations.
The coupling F’ is taken to be zero for simplicity.
The results of total energies and coupling constants for both families are given in
table 8. Several interesting features emerge. In the ‘2 x 1’ family the total energies for
the p ( 2 x 2) and c(4 x 2) structures are very similar. Since the error in the theoretical
calculations is of‘ the order of a few meV, it is impossible to tell which structure will
actually be realised at very low T. From a practical point of view, surface-preparation
conditions dictate the structure to be realised. A typical example is the cleaved Si
(111) surface which exhibits a (2 x 1) structure although the (7 x 7 ) structure is believed
to be slightly lower in energy. In the ‘c(2 x 2)’ family the energy difference between
structures is much larger. The corresponding phase transition in this family could be
predicted with certainty. However, since the c(2 x 2) structure is significantly higher
in energy than the p ( 2 x 2) or the c(4 x 2) structure in the (2 x 1) family, it is probably
not realised under normal conditions. Existence of steps in the surface further hampers
132 J Ihm
Table 8. Total energies for various reconstruction geometries of the Si(100) surface with
respect to the (2 x 1) surface and the coupling constants from Ihm et a [ (1983a, b, c).
‘(2 x 1)’ E (mev) ‘c(2 x 2)’ E (mev)
( 2 x 1) 0 c(2 x 2) 67
P(2 x 2) -36 P(2 x 2 ) 158
~(4x2) -3 1 ~ ( x 42 ) 135
( 4 x 1) 36 ( 4 x 1) 131
H 10 H’ 11
V -26 V’ -2
D 4 D’ 23
the realisation of the c(2 x 2) structure which usually has more dangling bonds across
the step and requires more energy to be created. If both (2 x 1) and c(2 x 2) families
happen to be created on the surface, the domain wall between these two structures is
going to be very stable. Breaking or displacing the wall costs a ‘huge’ amount of energy
(-1.7 eV per dimer). Thus each domain will behave independently, with different
phase transitions occurring in each.
The phase diagram for each family is shown in figures 5 ( a ) and 5 ( b ) . This is a
reduced parameter space obtained by dividing each coupling constant by kT (e.g.
d = D / kT, h H / kT, etc). The cross sections of constant d are shown to help visualise
the phase space. All phase boundaries shown in the figure are second order. As T
increases one moves toward the origin on a straight line in parameter space as indicated
in the figure.
For the ‘2 x 1’ family, this diagram shows that a second-order transition occurs
between layered antiferromagnetic p(2 x 2) and paramagnetic (disordered) phases at
roughly 250 K. As we have already mentioned, however, our theoretical calculations
Figure 5. Phase diagrams in reduced parameter space for ( a j ‘(2 x 1)’ and ( b ) ‘c(2 x 2)’
families from Ihm et a / (1983a, b, c). The straight lines passing through the origin corre-
spond to the Si(100) system as determined from the present calculations. P, paramagnetic;
F, ferromagnetic; AF, antiferromagnetic; LAF, layered antiferromagnetic.
are not accurate enough to distinguish between the p(2 x 2) and c(4 x 2) geometries.
Changes in energy of a few meV can turn the c(4 x 2) structure into the ground state
with a similar transition temperature to the paramagnetic phase. In any case, our
calculations predict that the pure 2 x 1 geometry is not the ground state and that
higher-order reconstructions of either ~ ( 2 x 2 or ) ~ ( 4 x 2 should
) be present at low
temperatures.
For the ' 4 2 x 2 ) ' family, the situation is much clearer. Figure 5 ( b ) predicts a
second-order transition between ferromagnetic c(2 x 2) and paramagnetic (disordered)
phases at -800 K. However, because the c(2 x 2) surface is significantly higher in
energy, the chance of observing the c(2 x 2) structure or the phase transition associated
with this structure would be small.
Since this work was published, there have been new developments on the subject.
Cardillo (1984) observed that the quarter-order spots of the Ge (100) surface corre-
sponding to the c(4 x 2) structure disappear much more rapidly than other spots at
elevated temperatures in his He atom diffraction experiments, suggesting an order-
disorder phase transition. Culbertson et a1 (1986) observed very clear c(4 x 2) quarter-
order spots at 100 K and the disappearance of them at 300 K on the Ge (100) surface.
More recently, Kevan and Stoffel (1984) and Kevan (1985) concluded, by monitoring
photoemission spectra in the energy gap, that there is an order-disorder transition on
the Ge (100) surface. Tabata et a1 (1985) reported the order-disorder phase transition
at 220 K from the c ( 4 x 2) to ( 2 x 1) structure on the Si (100) surface using the
temperature dependence of the LEED data. A Monte Carlo simulation (Saxena et a1
1985) using the parameters by Ihm et a1 (1983a, b, c) confirmed the accuracy of the
original renormalisation group calculation which predicted the critical temperature of
250 K.
There are also interesting total energy calculations on the reconstructed GaAs (100)
surfaces (Chadi et a1 1982, Ihm et a1 1983a) and the GaAs ( 1 1 1 ) (2 x 2) surface (Katnani
and Chadi 1985, Kaxiras et a1 1986) which will not be discussed in the present review.
3.4. Chemisorption, interface and superlattice

The chemisorption problem is perhaps most effectively coped with by the quantum-
chemical cluster calculation because of the very nature (the detailed local bonding
description) of the cluster method. The oxidation of the GaAs ( 1 10) surface (Barton
et a1 1979) is described as an example.
Calculations have been carried out for 0 bonded to As as a donor-acceptor complex
tabulated in table 7. Thus, upon oxidation, the GaAs relaxes from a twist angle 8 of
25.6" back to a twist of 18.6".
From the theoretical studies the As=O bond length of the donor-acceptor oxide
is expected to be 1.63 A (recall that the As-0 single bond is 1.80 A) oriented at 56"
with respect to the normal of the original unreconstructed surface. Preliminary results
from EXAFS studies (McMenamin et a1 1979) on the initial oxide layer on GaAs (110)
also indicate a short A s 0 bond length, -1.5 A. Extracting such bond lengths from
EXAFS requires phase shifts for scattering of the As x-ray electron by the neighbouring
0. Generally this phase shift is obtained from EXAFS studies of compounds with
known geometries. This is a problem for the As=O donor-acceptor bond since this
bond is much shorter than a normal As-0 bond and is quite different electronically.
Probably an effective experimental approach would be to determine the As-0 phase
134 J Zhm
shift from EXAFS studies on As203 and then use this phase shift for the three As-0
bonds to determine the As=O phase shift from EXAFS studies on As205.
It is also possible to compare core ionisation potential shifts with photoemission.
These calculations indicate that a donor-acceptor As=O bond shifts the As 3d orbitals
by 2.6 eV and the Ga 3d by 0.8 eV, both to higher binding energy. Thus, even though
0 is bonded only to As, it has a large effect upon Ga. The reasons is that bonding 0
to As (i) causes the As lone pair to rehybridise (become more p-like) and to delocalise
onto 0 to form a sigma bond; (ii) this rehybridisation affects the other three bonds to
As (less p) and causes the geometry to distort toward tetrahedral; (iii) the delocalisation
of the As lone pair toward 0 shifts the As core orbitals to deeper energy and polarises
the surface GaAs bonds, moving some charge from Ga to As and (iv) this latter effect
causes the Ga core orbitals also to shift to deeper energy. The experimental core shifts
are -2.9 and -1.0 eV for As and Ga, respectively (Brundle 1979, Goddard et a1 1976).
These values are in good agreement with the calculated values, supporting the sugges-
tion that the initial observed state of oxygenated GaAs (110) corresponds to the
donor-acceptor oxide.
Upon bonding of 0 to As of H,GaAs it is found that the empty orbital on the Ga
shifts down (higher binding energy) by 0.8 eV. Such a large shift downward is consistent
with the observed pinning of the Fermi energy upon oxidation of n-GaAs.
The interface problem has been studied as well with use of the total energy method.
A prototype problem is the Al-GaAs (1 10) interface formation at the initial stages of
A1 deposition. The determination of the A1 chemisorption site on the GaAs (110)
surface at the half-monolayer coverage (one A1 atom per GaAs surface unit cell) has
been a long-standing problem. However, Zunger (1981) proposed A1 cluster formation
with a very weak interaction between the adatom cluster and the substrate surface. A
complete energy minimisation calculation (Ihm and Joannopoulos 1981b, 1982) cover-
ing all the possible chemisorption sites plus the AI-A1 cluster geometry confirmed
Zunger’s proposal. Possible migration paths of A1 on the GaAs (110) surface to form
A1 clusters was also obtained. At elevated temperatures, an AI-Ga exchange interaction
was shown to be energetically most favourable, in agreement with a simple heat-of-
formation argument and also with the L E E D and the soft x-ray photoemission spectros-
copy (Duke et a1 1981). Details of the calculation are described below.
The geometries to be discussed here can be grouped into three classes: A1 chemisorp-
tion on the substrate, AI-A1 bonding, and Al-substrate exchange. In the first case,
one A1 atom exists per surface unit cell. Since the nearest-neighbour A1-A1 distance
(4.0A) is quite large, the major binding occurs only between the A1 atom and the
substrate. First-principles pseudopotential total energy calculations have been perfor-
med for nearly forty different arrangements of atoms, and the adsorption-energycontour
map for the entire surface is shown in figure 6. It is found that site cy gives the largest
chemisorption energy (2.3 eV) in this class of geometries. The substrate atoms tend
to recover the original bulk positions (in the ideal tetrahedral coordination) when A1
atoms come and saturate broken bonds on the surface. In fact, such an ‘unrelaxation’
of the substrate upon A1 adsorption lowers the energy in most cases studied.
*
The AI-As bond length is 2.38 0.03 A and the AI-Ga bond length is 2.49 0.03 A*
in this lowest energy chemisorption geometry. The A1 coordinate in the [ 1101 direction
is 0.65 A. The Al-substrate bond length turns out to be 2.45i0.1 A for other sites as
well. These calculated bond lengths and binding energies, being almost identical to
the corresponding bulk bond lengths (2.45 A) and binding energies (1.9 eV), are in
contrast to previous theoretical calculations (Barton et a1 1980, Zunger 1981) predicting
much longer and weaker bonding. Because the substrate atoms on the surface restore
the tetrahedral bonding configuration by A1 adsorption, the binding energy for the
chemisorbed A1 atom is expected to be comparable to the bulk AI-As binding energy
as obtained in the present study.
Two A1 atoms are then introduced per unit cell and the total energy is minimised
with respect to the adatom positions. In the lowest energy configurations, one A1 atom
occupies site /3 and the other occupies site y in figure 6. Therefore, A1 atoms saturate
all the broken bonds and form AI-A1 bonding chains. The corresponding binding
energy is 2.7 eV per AI atom. A moderately covalent bond between A1 atoms is formed
in this case. We also have calculated the cohesive energy of bulk A1 (to be regarded
as the large cluster limit) and obtained 3.6 eV in good agreement with the experimental
value of 3.4 eV per A1 atom.
Comparing binding energies calculated so far, it is concluded that any A1 chemisorp-
tion shown in figure 1 is unstable with respect to A1-A1 bond formation, and that a
large cluster is more stable than a single chain of A1 bonds. Because the ratio of the
lattice constant of GaAs to that of A1 is &' with only 1% error, it is likely that A1
clusters line u p in a preferential direction on the GaAs surface to minimise the strain
energy between them.
The A1 substrate exchange geometry is finally considered. Since the magnitude of
the heat of formation for AlAs is around 0.5 eV larger than GaAs per pair, it can be
immediately concluded that the A1 atoms eventually penetrate the GaAs surface and
replace the substrate G a atoms to make AI-As bonds in place of Ga-As bonds. Total
energy calculations for exchanged geometries have been done and compared with
corresponding geometries without exchange. The energy is lowered by 0.48 eV for the
- iiioi
Joannopoulos (1982). As, 0 ; Ga ..

Figure 6 . AI adsorption energy contour plot on the GaAs (110) surface from Ihm and
The substrate surface lies in the plane of the figure
and the adsorbed AI atom lies above the plane. One AI atom is assumed to exist per
surface unit cell. Important local extrema are labelled with Greek letters, a, 2.3 eV; /3,
1.9eV; y, 1.7eV; 6, 1.7eV; E , 1.SeV; 4', 1.SeV.
136 J Ihm
exchange of an A1 atom with the Ga atom in the first layer, and by 0.64eV for the
exchange in the second or third layer. These results are not surprising since the cohesive
energy difference between AlAs and GaAs was reproduced in a previous calculation
using the same method (Ihm and Joannopoulos 1981a).
Then the A1 atom migration on the surface at finite temperatures is studied. The
tunnelling probability is in general proportional to exp(- VB/ k T ) , where V , is the
potential barrier from one configuration to another. We are in a situation where
configurations with lower free energy are not easily accessible (i.e. they require activa-
tion energy), and it takes a finite time for the system to reach equilibrium. Under these
circumstances, describing how A1 atoms move around to arrive at lower energy
configurations is more meaningful than deriving a phase diagram which assumes a
thermal equilibrium. A non-trivial consequence immedia.tely conceivable is that at
very low temperature the A1 atoms do not move and remain chemisorbed. But at room
temperatures or higher, where experiments for this system are usually done, the cluster
formation or the exchange reaction does occur. Although the total energy calculation
is static in nature, it is capable of providing dynamical information in conjunction
with the Born-Oppenheimer approximation for the atomic motion as described below.
Let us consider the motion of A1 atoms on the GaAs surface. Going back to figure
6, two preferable channels for migration of A1 atoms can be identified. One such
channel is apparent between two zigzag chains of GaAs atoms (along the sites y - a - p -
a-y). This channel exists because the zincblende structure has an open space inside
it. Another channel exists above and along the zigzag chain (along the sites C - E - ~ - E - ~ ) .
Although this channel is higher in total energy than the first one, the activation barrier
along this channel is smaller. This is because the substrate atoms are located rather
compactly along the chain, and the chemical bonding between the adsorbed A1 and
the substrate atoms does not experience much change as the A1 atom proceeds along
the chain.
To help visualise the situation clearly, the same plot is drawn in a three-dimensional
view in figure 7. Two layers of the GaAs (110) surface system are also shown. This
figure illustrates how an A1 atom 'senses' the GaAs (110) surface. A1 atoms can migrate
through the deep valley as indicated by path 1. The potential barrier experienced by
an A1 atom travelling along this path is -0.5 eV. The second channel exists along the
plateau as indicated by path 2 . It provides a nice and smooth route for the surface
migration of A1 atoms and the barrier is within 0.3 eV. We note, however, that these
two channels are not completely separate. An A1 atom in channel 2 can slide down
and end up in channel 1. Nevertheless, it is obvious that the A1 atoms tend to move
favourably in a direction parallel to the chain (the [ i l O ] direction) rather than perpen-
dicular to it. This anisotropy in surface diffusion should be observable experimentally
unless the surface is significantly disordered.
The surface hopping diffusivity is expressed as ud2 exp(- VB/ k T ) where U is the
atomic vibration frequency and d is the lattice spacing. If the product of U and
exp( - VB/ k T ) is much greater than 1, the migration of A1 atoms is very active, and
the clusters which have lower free energy can be realised on the laboratory timescale.
The frequency v is typically -10'2s-', and clustering would be observable at room
temperature if VB is 0.6 eV or less. Since VB is smaller than 0.6 eV for both channels,
clustering should be observable immediately.
The exchange reaction, on the other hand, requires a large perturbation on the
surface. The activation energy should be as large as -1 eV even if due consideration
is given for possible surface imperfections which lower the activation barrier. Thus,
Figure 7. Three-dimensional plot of the total energy for an AI atom adsorbed on the GaAs
(110) surface from Ihm and Joannopoulos (1982). Ga, 0; As, 0 . Two favourable paths
for the surface migration of AI atoms are indicated as 1 and 2. It is clear in the figure that
channel 1 follows the valleys and channel 2 traces plateaus. A perspective view of atomic
positions for the first two surface layers is also illustrated.
the exchange reaction is active only at high temperatures as observed experimentally

(Duke et a1 1981).
Recently, the application of the total energy method has been extended to the
strained-layer superlattice system. Van der Walle and Martin (1985) have studied the
structural and electronic properties of pseudomorphic Si-Ge interfaces, where the
lattice spacing parallel to the interface is forced to be equal on both sides. They have
determined the minimum energy configurations and the relative line-up of Si and Ge
bands at the interface using the LDF first-principles pseudopotential method. This is
a new area of active research in solid state physics and more research results are
expected to be produced in the near future.
3.5. Defects and others

Total energy calculations have been performed for various kinds of defects in the last
few years with great success by many different groups. In this section, a few such
works will be described briefly.
When the silicon vacancy was proposed by the Bell Labs group (Baraff et al 1979)
to be a possible ‘Anderson negative U’ system (Anderson 1975) where the two
138 J Ihm
opposite-spin electrons are paired because of the net attractive interaction between
them, the theoretical scheme to estimate the total energy was still an empirical method
based on a specific force constant model (Keating 1966). Since then, there has been
remarkable progress in theoretical methods and in the understanding of defects on a
microscopic basis. The Bell Labs group (Baraff et a1 1981) developed a systematic
method to study defects (or impurities) using Green function and total energy calcula-
tions and predicted very weak G a - 0 bonds in the Gap: 0 system. Their scheme has
been further revised to incorporate the state of the art first-principles Green function
total energy method (Baraff and Schluter 1983), allowing for the study of isolated,
charged point defects.
The MIT group devised an ingenious first-principles pseudopotential method to
study a pair of oppositely charged defects using an idea of density matrices and applied
it to calculations of defects in selenium (Vanderbilt and Joannopoulos 1983), in contrast
to their previous work with a more empirical method (Vanderbilt and Joannopoulos
1979). A more sophisticated calculation of the silicon self-interstitial migration was
successfully done more recently (Bar-Yam and Joannopoulos 1984).
The IBM group has also been very active in studying the defect problem by means
of the total energy method. They proposed a microscopic theory of atomic diffusion
mechanisms in silicon (Car et a1 1984) or impurity-defect reactions and impurity
diffusion in silicon (Car et a1 1985) based on the self-consistent Green function total
energy calculations.
Defect calculations have been carried out with use of semiempirical methods (Masri
et al 1983) or unrestricted Hartree-Fock self-consistent-field methods (Sahoo et al
1983) in a finite cluster as well. A self-consistent, parameter-free approach based on
the Green function method was developed which allowed one to calculate the lattice
distortion around neutral and charged simple defects from the Hellmann-Feynman
theorem and applied to the Ga site single donor Ge in GaP (Scheffler et a1 1982).
There was a recent proposal that an exhange-correlation induced 'negative U' be
possible (Katayama-Yoshida and Zunger 1985). This was illustrated for an unrelaxed
interstitial Cr impurity in Si through a self-consistent local spin density Green function
calculation.
Dislocation has been studied by Northrup et a1 (1981) which will not be discussed
here. Ho et a1 (1980) have used the first-principles pseudopotential method to study
the effect of electromodulation of surface states on the electroreflectance of the Ag
(110) surface.
4. Discussion and conclusions
We have examined applications of the total energy methods to variety of fields in solid
state physics. In most cases, the results are impressive and the computational efforts
are rewarding. Looking back at the first work in the field by Wigner and Seitz (1933),
we recognise that the total energy calculation was based on the simple fact of statistical
physics that the equilibrium state of a system at the absolute zero temperature is the
state of the minimum total energy. Now the same total energy calculation has been
used to predict the order-disorder phase transition of the Si (100) surface at finite
temperatures and to describe such dynamic effects as the diffusion of the silicon
self-interstitials. Several different first-principles computational methods are available
nowadays and the claimed accuracy of the calculation is almost as good as the accuracy
in experiment in many cases.
It is appropriate on this occasion to point out a few outstanding problems in the

field. Each calculational method has its strong points and weak points. One of the
problems common to most methods is that as the number of inequivalent atoms in the
system increases, the computation of the total energy is quite time-consuming and
becomes intractable very soon even with the fastest computer available at the present.
The reduction of the problem to a few-body interaction, if successful, would simplify
the calculation dramatically and it would enable us to study such complicated systems
as amorphous materials. Such attempts are promising, but have not been quite success-
ful yet (Biswas and Hamann 1985).
The exchange-correlation potential is another ongoing problem of the research.
Some of the published works have been described in the introduction and there are
more expressions for the exchange-correlation potentials available now (Ceperley 1978,
Langreth and Perdew 1980, Langreth and Mehl 1981). There seems to be a limitation
to the description of the exchange-correlation effect using the one-electron effective
potential. A closely related problem is the accuracy of the density functional formalism
(Kohn and Sham 1965) widely used in total energy methods. In particular, the local
approximation to the true density functional formalism has been used in most practical
total energy calculations for computational convenience and the validity of the LDF
formalism has been questioned occasionally. It is a well known fact that the silicon
minimum gap obtained from the LDF approximation is 0.5eV compared with the
experimental values of -1.1 eV. This is not a totally unexpected result considering
that the density functional formalism is exact only for the description of the ground
state and there is no guarantee that it works for the excited state. There have been
successful attempts to calculate the gap correctly using methods beyond the simple
L.DF approximation (Wang and Pickett 1983, Sham and Schliiter 1983, Hybertsen and
Louie 1984). It is rather disappointing, though, to realise that the LDF approximation
is not capable of giving correct results for the valence band offset at the interface which
is genuinely a ground-state property of the system. This conclusion can be drawn from
a recent work (Hybertsen and Louie 1984) showing that the entire valence band shifts
by a different amount for different materials (e.g. between silicon and diamond) when
the exchange-correlation effect is accurately treated. This result casts doubt on band
offset calculations using the LDF approximations. The correct band offset calculations
beyond the LDF approximation for realistic systems seem to be too difficult at this point.
Even with these outstanding problems in the field, the future of the total energy
calculations is very bright. Physicists are in general not satisfied with hand waving
arguments or qualitative results and demanding more quantitative results from the
first-principles calculations these days. The modern total energy methods are supplying
such quantitative results directly comparable with experimental measurements. It
would be too pessimistic a viewpoint to worry that the physics would be lost in the
computer in the midst of the total energy calculations. Physical insights are always
required to set up the computational problem properly and extract physically meaning-
ful results from the calculations. Otherwise, a new area of physics which can be called
‘computational physics’ in contrast to theoretical physics and experimental physics,
to which the total energy methods probably belong, will not survive as a genuine
branch of physics in the future.
Acknowledgments
I thank Drs D J Chadi, A Redondo, J J Barton, W A Goddard I11 and T C McGill

for allowing me to quote their original works.
140 J Ihm
References
Allan D C and Mele E J 1984 Phys. Rev. Lett. 53 826

Andersen 0 K 1975 Phys. Rev. B 12 3060
Andersen 0 K and Jepsen 0 1977 Physica 91B 317
Anderson P W 1975 Phys. Rev. Lett. 34 953
Appelbaum J A and Hamann D R 1973 Phys. Rev. B 8 1777
Averill F W 1972 Phys. Rev. B 6 3637
Bachelet G B, Greenside H S, Baraff G A and Schluter M 1981 Phys. Rev. B 24 4745
Baraff G A, Kane E 0 and Schliiter M 1979 Phys. Rev. Lett. 43 956
-1981 Phys. Rev. Lett. 47 601
Baraff G A and Schliiter M 1983 Phys. Rev. B 28 2296
Barton J J, Goddard W A 111 and McGill T C 1979 J. Vac. Sci. Technol. 16 1178
Barton J J, Swarts C A, Goddard W A I11 and McGill T C 1980 J. Vac. Sci. Technol. 17 164, 869
Bar-Yam Y and Joannopoulos J D 1984 Phys. Rev. Lett. 52 1129
Biswas R and Hamann D R 1985 Phys. Rev. 'Lett. 55 2001
Born M and Misra R D 1940 Proc. Camb. Phil. Soc. 36 466
Born M and Oppenheimer J R 1927 A m . J. Phys. 84 451
Brooks H 1963 Trans. Metallurg. Soc. A I M E 221 546
Brundle C R 1979 J. Vac. Sci. Technol. 16 1186
Car R, Kelley P J, Oshiyama A and Pantelides S T 1984 Phys. Rev. Lett. 52 1814
-1985 Phys. Rev. Lett. 54 360
Cardillo M J 1984 private communication
Cardillo M J and Becker G E 1978 Phys. Rev. Lett. 40 1148
-1980 Phys. Rev. B 21 1497
Ceperley D 1978 Phys. Rev. B 18 3126
Ceperley D M and Alder B J 1980 Phys. Rev. Lett. 45 566
Chadi D J 1978 Phys. Rev. Lett. 41 1062
-1979a Phys. Rev. Lett. 43 43
- 1979b J. Vac. Sri. Technol. 16 1295
Chadi D J, Tanner C and Ihm J 1982 Sui$ Sci. 120 L425
Chang K J, Dacorogna M M, Cohen M L, Mignot J M, Choteau G and Martinez G 1985 Phys. Rev. Lett.
54 2375
Chelikowsky J R and Burdett J K 1986 Phys. Rev. Lett. 56 961
Chelikowsky J R and Louie S G 1984 Phys. Rev. B 29 3470
Ching W Y and Callaway J 1974 Phys. Rev. B 9 5115
Cohen M L and Heine V 1970 Solid State Phys. 24 37
Coldwell-Horsfall R A and Maradudin A A 1960 J. Math. Phys. 1 395
Culbertson R J, Kuk Y and Feldman L C 1986 Surf: Sci. 167 127
Danese J B 1974 J. Chem. Phys. 61 3071
Dolling G 1963 Inelastic Scattering of Neutrons in Solids and Liquids vol 2 (Vienna: IAEA) p 37
Duke C B, Lubinsky A R, Boon M, Cisneros G and Mark P 1977 J. Vac. Sci. Technol. 14 294
Duke C B, Paton A, Meyer R J, Brillson L J, Kahn A, Kanani D, Carelli J, Yeh J L, Margaritondo G and
Katnani A D 1981 Phys. Rev. Lett. 46 440
Ewald P P 1921 Ann. Phys., Lpz. 64 253
Feibelman P J and Hamann D R 1983 Phys. Rev. B 28 3092
Feynman R P 1939 Phys. Rev. 56 340
Froyen S and Cohen M L 1982 Solid State Commun. 43 447
Fuchs K 1935 Proc. R. Soc. A 151 585
Goddard W A 111, Barton J J, Redondo A and McGill T C 1978 J , Vac. Sci. Technof. 15 1274
Goddard W A 111, Redondo A and McGill T C 1976 Solid State Commun. 18 981
Goroff I and Kleinman L 1970 Phys. Rev. B 12574
Gunnarsson 0, Lundqvist B I and Wilkins J W 1974 Phys. Rev. B 10 1319
Hamann D R 1979 Phys. Rev. Lett. 42 662
Hamann D R, Schliiter M and Chiang C 1979 Phys. Rev. Lett. 43 1494
Haneman D 1961 Phys. Rev. 121 1093
Harris F E and Monkhorst H J 1971 Solid State Commun. 9 1449
Harris J and Jones R 0 1978 Phys. Rev. Lett. 41 191
Hedin L and Lundqvist B I 1971 J. Phys. C: Solid State Phys. 4 2064
Hellmann H 1937 Einfuhrung in die Quantenchemie (Leipzig: Deuticke)
Herring C 1940 Phys. Rev. 57 1169

Ho K M, Fu C L, Harmon B N, Weber W a n d Hamann D R 1982 Phys. Rev. Lett. 49 673
H o K M, Harmon B N and Liu S H 1980 Phys. Rev. Lett. 44 1531
Hohenberg P and Kohn W 1964 Phys. Rev. 136 B864
Hybertsen M S and Louie S G 1984 Phys. Rev. B 30 5777
Ichikawa T and Ino S 1979 Surf: Sci. 85 22
Ihm J, Chadi D J and Joannopoulos J D 1983a Phys. Rev. B 27 5119
Ihm J and Cohen M L 1979 Solid Stale Commun. 29 71 1
-1980 Phys. Rev. B 21 1527
_. 1981 Phys. Rev. B 23 1576
Ihm J, Cohen M L and Chadi D J 1980a Phys. Rev. B 21 4592
Ihm J and Joannopoulos J D 1981a Phys. Rev. B 24 4191
-1981b Phys. Rev. Lett. 47 679
-1982 Phys. Rev. B 26 4429
Ihm J, Lee D H, Joannopoulos J D and Berker A N 1983b J. Vac. Sci. Technol B 1 705
Ihm J, Lee D H, Joannopoulos J D and Xiong J J 1983c Phys. Rev. Left. 51 1872
-1985 Proc. 17th Int. Con$ Physics of Semiconductors ed D J Chadi and W A Harrison (New York:
Springer) p 51
Ihm J, Yin M T and Cohen M L 1981 Solid State Commun. 37 491
Ihm J, Zunger A and Cohen M L 1979 J. Phys. C: Solid State Phys. 12 4409
-1980b J. Phys. C: Solid State Phys. 13 3095
Jamieson J C 1963 Science 139 762, 845
Janak J F 1974 Phys. Rev. B 9 3985
Janak J F, Moruzzi V L and Williams A R 1975 Phys. Rev. B 12 1257
Janak J F and Williams A R 1974 Phys. Rev. B 9 2670
-1976 Phys. Rev. B 14 4199
Kane E 0 1966 Semiconductors and Semimetals vol 1, ed R K Willardson and A C Beer (New York:
Academic) p 75
Katayama-Yoshida H and Zunger A 1985 Phys. Rev. Lett. 55 1618
Katnani A D and Chadi D J 1985 Phys. Rev. B 31 2554
Kaxiras E, Bar-Yam Y , Joannopoulos J D and Pandey K C 1986 Phys. Rev. B 33 4406
Keating P N 1966 Phys. Rev. 145 637
Kerker G 1980 J. Phys. C: Solid State Phys. 13 L189
Kevan S D 1985 Phys. Rev. B 32 2344
Kevan S D and Stoffel N G 1984 Phys. Rev. Lett. 53 702
Kohn W and Rostoker N 1954 Phys. Rev. 94 1111
Kohn W a n d Sham L J 1965 Phys. Rev. 140 A1133
Korringa J 1947 Physica 13 392
Kunc K and Martin R M 1981 Phys. Rev. B 24 2311
Lam P K and Cohen M L 1981 Phys. Rev. B 24 4224
Lander J J and Morrison J 1962 J. Chem. Phys. 37 729
Langreth D C and Mehl M J 1981 Phys. Rev. Lett. 47 446
Langreth D C and Perdew J P 1980 Phys. Rev. B 21 5469
Lubinsky A R, Duke C B, Lee B W and Mark P 1976 Phys. Rev. Lett. 36 1058
Mark P, Kahn A, So E and Duke C B 1978 Bull. Am. Phys. Soc. 23 400
Masri P, Harper A H and Stoneham A M 1983 J. Phys. C: Solid State Phys. 16 L613
McMahan A K, Yin M T and Cohen M L 1981 Phys. Rev. B 24 7210
McMenamin J C , Bauer R S, Stohr J and Johansson L 1979 J. Vac. Sci. Technol. 16 1195
Moriarty J A 1974 Phys. Rev. B 10 3075
-1977 Phys. Rev. B 16 2537
Morita A, Soma T and Takeda T 1972 J. Phys. Soc. Japan 32 1972
Needs R J, Martin R M and Nielsen 0 H 1986 Phys. Rev. B 33 3778
Nielsen 0 H and Martin R M 1985 Phys. Rev. B 32 3780, 3792
Nilson G and Nelin G 1972 Phys. Rev. B 6 3777
Northrup J E and Cohen M L 1982 Phys. Rev. Lett. 49 1349
Northrup J E, Cohen M L, Chelikowsky J R, Spence J C H and Olsen A 1981 Phys. Rev. B 24 4623
Pandey K C 1981 Phys. Rev. Lett. 47 1913
-1982 Phys. Rev. Lett. 49 223
Phillips J C and Kleinman L 1959 Phys. Rev. 116 287
Poppendieck T D, Ngoc T C and Webb M B 1978 Surf: Sci. 75 287
142 J Zhm
Posternak M, Krakauer H, Freeman A J and Koeling D D 1980 Phys. Rev. B 21 5601

Redondo A and Goddard W A 111 1982 J. Vac. Sci. Technol. 21 344
Redondo A, Goddard W A 111, McGill T C and Surrat G T 1976 Solid State Commun. 20 733
Sabin J R, Worth J P and Trickey S B 1975 Phys. Rev. B 11 3658
Sahoo N, Mishra S K, Mishra K C, Coker A, Das T P, Mitra C K, Snyder L C and Glodeanu A 1983 Phys.
Rev. Left. 50 913
Saxena A, Gawlinski E T and Gunton J D 1985 Surf: Sci. 160 618
Schemer M, Vigneron J P and Bachelet G B 1982 Phys. Rev. Left. 49 1765
Schliiter M, Chelikowsky J R, Louie S G and Cohen M L 1975 Phys. Rev. B 12 4200
Seitz F 1935 Phys. Reo. 47 400
Sham L J and Schliiter M 1983 Phys. Rev. Lett. 51 1888
Simons G and Bloch A N 1973 Phys. Rev. B 7 2754
Singwi K S, Sjolander A, Tosi M P and Land R H 1970 Phys. Rev. B 1 1044
Slater J C 1937 Phys. Rev. 51 846
-1972 J. Chem. Phys. 57 6
Slater J C and Johnson K H 1972 Phys. Rev. B 5 844
Slater J C and Koster G F 1954 Phys. Rev. 94 1498
Snow E C 1973 Phys. Rev. B 8 5391
Starkloff T and Joannopoulos J D 1977 Phys. Reo. B 16 5212
Tabata T, Aruga T and Murata Y 1985 Techn. Rep. ISSP Ser. A 1595
Tong S Y, Lubinsky A R, Mrstik B J and Van Hove M A 1978 Phys. Rev. B 17 3303
Tromp R M, Hamers R J and Demuth J E 1985 Phys. Rev. Lett. 55 1303
Van Camp P E, Van Doren V E and Devreese J T 1979 Phys. Rev. Lett. 42 1224
Vanderbilt D and Joannopoulos J D 1979 Phys. Rev. Lett. 42 1012
-1980 Solid State Commun. 35 535
-1983 Phys. Rev. B 27 6302
Van der Walle C G and Martin R M 1985 J. Vac. Sci. Technol. B 3 1256
van Leeuwen J M 1975 Phys. Rev. Lett. 34 1056
Verges J A and Tejedor C 1979 Phys. Rev. B 20 4251
Vogl P, Hjalmerson H P and Dow D J 1983 J. Phys. Chem. Solids 44 364
von Barth U and Gelatt C D 1980 Phys. Rev. B 21 2222
von Barth U and Hedin L 1972 J. Phys. C: Solid State Phys. 5 1629
Wang C S and Pickett W E 1983 Phys. Rev. Lett. 51 597
Wang Y and Joannopoulos J D 1980 J. Vac. Sci. Technol. 17 997
Weaire D 1970 Phys. Status Solidi 42 767
Weber W 1974 Phys. Rev. Lett. 33 371
-1977 Phys. Rev. B 10 4789
Weinstein B A and Piermarini G J 1975 Phys. Rev. B 12 1172
Wendel H and Martin R M 1978 Phys. Reo. Lett. 40 950
Wepfer G G, Euwema R N, Surratt G T and Wilhite D L 1974 Phys. Reo. B 9 1670
Wigner E 1934 Phys. Rev. 46 1002
Wigner E and Seitz F 1933 Phys. Rev. 43 804
-1934 Phys. Rev. 46 509
Yin M T and Cohen M L 1980 Phys. Rev. Lett. 45 1004
-1981a Solid State Commun. 38 625
-1981b Phys. Rev. B 24 2303
-1981c Phys. Rev. B 24 6121
-1982 Phys. Rev. B 26 5675
Zunger A 1981 Phys. Rev. B 24 4372
Zunger A and Cohen M L 1979a Phys. Rev. B 19 568
-1979b, Phys. Rev. B 20 4082
Zunger A and Freeman A J 1977a Phys. Reo. B 16 906
-1977b Phys. Rev. B 16 2901
-1978 Phys. Rev. B 17 2030

Ihm 1988

Hochgeladen von

Copyright:

Verfügbare Formate

Ihm 1988

Hochgeladen von

Dokumentinformationen

Originaltitel

Copyright

Verfügbare Formate

Dieses Dokument teilen

Dokument teilen oder einbetten

Freigabeoptionen

Stufen Sie dieses Dokument als nützlich ein?

Sind diese Inhalte unangemessen?

Copyright:

Verfügbare Formate

Ihm 1988

Hochgeladen von

Copyright:

Verfügbare Formate

Home Search Collections Journals About Contact us My IOPscience

Total energy calculations in solid state physics

1988 Rep. Prog. Phys. 51 105

Please note that terms and conditions apply.

Total energy calculations in solid state physics

This review was received in its present form in July 1986.

0034-4885/88/010105 +38$09.00 @ 1988 IOP Publishing Ltd 105

2.1. Tight-binding method

Et,,= E,,+E,i+ Eii (2.1)

= Ebs + ( Eii - (2.3)

U, -23.19 -17.72 -14.15

a (A) 5.44 5.65 5.65

2.2 Quantum-chemical cluster calculations

a(4(1)4(2)a(l)P(2)) = 4(1)4(2)(aP - P a ) (2.8)

Q H F= 4(1)4(2)= (416 + 474,) + (4141+ 66) (2.10)

= (4r4r+ 4r4r>(.P - P a ). (2.11)

As a result the H F wavefunction leads to a terrible description at large R, as indicated

configurations of As where [Nil is given in (2.23). The core potential is taken to be

2.3. Pseudopotential method

and the electrostatic potential energy, V, is

V c o u l ( r=) 2p0d3r’= 1 VC,,,(C ) eic’r (2.38)

The Fourier component of the inter-electronic Coulomb potential, Vco,,( G), is

x UPJ r ) A exp[i( k, + G') ( r + R,)] d3r

where S ( G'- G ) is the structure factor and

for small G. From (2.41) we have

p is a rather complicated integral. Fortunately, p does not appear in the final

2.4. Other theoretical tools

3.1. Structurally related properties of bulk materials

Lattice Crystal Cohesive Bulk

Calculation 5.451 -7.909 4.67 0.98

v: VP VPl v:' P, (kbar)

Calculation 0.928 0.718 0.774 99

a Jamieson (1963), Weinstein and Piermarini (1975).

volume reduction by 23%, in remarkable agreement with experiment ! When superstress

LTo(r) TAW TdX) LOA(X)

fcak 15.16 4.45 13.48 12.16

Measured at 296 K (Dolling 1963, Nilson and Nelin 1972).

Reference a Reference b Reference c Reference d Reference a

0 27.3" 26.4" 34.8" 27.1" 25.6"

Orbital energies (Hartree)a

The first-principles pseudopotential calculations (Ihm and Joannopoulos 1981b,

differences between different reconstructions are translated into a set of interaction

in agreement with both pseudopotential calculations and experimental observations.

‘(2 x 1)’ E (mev) ‘c(2 x 2)’ E (mev)

3.4. Chemisorption, interface and superlattice

Joannopoulos (1982). As, 0 ; Ga ..

the exchange reaction is active only at high temperatures as observed experimentally

3.5. Defects and others

4. Discussion and conclusions

It is appropriate on this occasion to point out a few outstanding problems in the

I thank Drs D J Chadi, A Redondo, J J Barton, W A Goddard I11 and T C McGill

Allan D C and Mele E J 1984 Phys. Rev. Lett. 53 826

Herring C 1940 Phys. Rev. 57 1169

Posternak M, Krakauer H, Freeman A J and Koeling D D 1980 Phys. Rev. B 21 5601

Das könnte Ihnen auch gefallen