Thesis
The electronic properties and optimized structures of
the alkali adsorbed Si(001) surface by using the rst
principles molecular dynamics

Kazuaki Kobayashi
Institute for Solid State Physics
University of Tokyo
December 25, 1990

Contents
1

Introduction

2

2

Method of calculation

2.1 Pseudopotential . . . . . . . . . . . . . . . . . . . . . . . .
2.1.1 Norm-conserving pseudopotential . . . . . . . . . .
2.1.2 Kleinman-Bylander separable form . . . . . . . . .
2.2 Total energy calculation . . . . . . . . . . . . . . . . . . .
2.3 Hellmann-Feynman force . . . . . . . . . . . . . . . . . . .
2.4 Car-Parrinello method . . . . . . . . . . . . . . . . . . . .
2.4.1 Steepest descent method . . . . . . . . . . . . . . .
2.4.2 Modi ed steepest descent method . . . . . . . . . .
2.5 Molecular dynamics for ionic part . . . . . . . . . . . . . .

8

8
8
11
14
18
20
22
24
27

3

Results and discussion

28

4Concluding remarks

40

3.1 Test calculations for bulk systems . . . . . . . . . . . . . .
3.2 Calculation of Si(001) surface . . . . . . . . . . . . . . . .
3.2.1 Clean surface . . . . . . . . . . . . . . . . . . . . .
3.2.2 Full coverage Na/Si(001)-221 surface . . . . . . . .
3.2.3 Half coverage Na/Si(001)-221 surface . . . . . . . .
3.2.4 61 coverage Na/Si(001)-221 surface . . . . . . . . .
3.2.5 Full coverage Li/Si(001)-221 surface . . . . . . . .
3.2.6 16 coverage Li/Si(001)-221 surface . . . . . . . . . .
3.3 Partial core correction . . . . . . . . . . . . . . . . . . . .

28
29
30
31
33
35
37
38
38

A

43

B

44

1

C

46

2

Chapter 1

Introduction
The structural determination is a very important rst step toward understanding the material properties from a microscopic view point. For perfect
crystal, very accurate structural information can be obtained by means of
the x-ray and neutron di raction experiments. On the other hand, thestructures with defects, surfaces and interfaces are generally very hard
to be determined experimentally. Therefore, it is highly desirable thattheories provide some useful complementary information in those cases.
However, it is necessary to consider a very large system in the theoreticalstructural determination.
In traditional and standard electronic structure calculations based on thelocal density approximation (LDA) in the density functional (DF) theory,
it is possible to obtain an accurate electronic properties (charge density,
dispersion curves, density of states, etc.). Several works have also demonstrated that the total energy di erence between di erent atomic structures
can be calculated fairly accurately by LDA [1, 2]. Therefore, it is in principle possible to determine the lowest energy structure by calculating thetotal energy for every possible structure. However, this brute force approach becomes practically impossible for systems with more than threedegrees of freedom.
It is necessary to calculate the forces to determine the structure efciently. In the standard electronic structure calculation by using thenorm-conserving pseudopotential and plane wave basis sets, it is easy to
calculate the Hellmann-Feynman forces [3, 4, 5] which are acting on theions. However, such an approach so far has generally used the standard
matrix diagonalization method in order to obtain the eigen-vectors andeigen-values. The matrix diagonalization method requires an N3 order2

task of calculation and an N2 order memory with N denoting the numberof plane waves. Therefore, the limitation of system size imposed by these
aspects is also severe in practice.
On the other hand, a standard molecular dynamics can be used to obtain
the optimized structure for a relatively large system. However, in thismethod, empirical or semi-empirical potentials are adopted whose accuracy
is not necessarily guaranteed. In addition, the electronic properties of the
optimized system are not given by this method.
Consequently, it is impossible to optimize the electronic and ionic structures simultaneously for a large system by using the traditional calculational methods. It has been hoped that a new method which solves above
problems would appear.
In 1985, a new method was developed by Car and Parrinello [6]. Usually,
this is called `Car-Parrinello method', but in this thesis, it will be called`The rst principles molecular dynamics (FPMD) method'. The conceptof this method is described as follows. In the density functional theory, thetotal energy of the system for a xed ionic con guration takes a minimum
for the true ground state. Therefore, solving the self-consistent solution of
the ground state electronic structure can be regarded as an optimization of
the total energy functional with regard to the electronic degrees of freedom.
Based on this idea, the FPMD method calculates the electronic and ionicparts simultaneously. The electronic part is treated in a similar way to theionic part by introducing an equation of motion like a Newton equation andits dynamics is totally ctitious. The wave function is solved dynamicallyby integrating the di erential equation given by a Lagrange equation of
motion and the process does not include matrix diagonalization. Therefore,
this FPMD method requires no more than an N2 order task of calculation
and an N order memory, i.e., this method is very fast and saves a memory
compared with the traditional one. Furthermore, it is possible to reduce
the task of calculation from N2 to NlogN( : numerical factor) as a resultof some technical considerations. It is possible to treat larger systems,
which cannot be imagined in the usual method.
In this thesis, we shall develop the program of the FPMD method andapply it to the Si(001) clean surface and the alkali adsorbed Si(001) surface systems. The surface is a very interesting system because it provides
varieties of interesting phenomena, such as catalysis, anticorrosion, silicide3formation at metal-semiconductor interfaces important in device applications, epitaxial growth, reconstruction and so on. Recently it has become
possible to perform very accurate experiments for clean surface systems because of the development of the ultra high vacuum technique. In addition,a new experimental technique `STM(Scanning Tunneling Microscope)' appeared in 1982 [7], which can clarify a local structure with single-atom
resolution. Therefore, the research of surface systems has been developed
remarkably both theoretically and experimentally.
Within the framework of the standard electronic structure calculation,
it is very hard to calculate the surface systems because these semi-in nite
systems do not have a periodicity. In order to solve this problem, a repeated slab geometry and supercell method are frequently adopted. The
supercell is a very large unit cell which includes several tens and sometimes
more than a hundred of atoms. In the repeated slab geometry, the unit
cell includes atomic layers and a vacuum region. Since this unit cell is verylarge, it is necessary to use a huge number of plane waves in this calculation. Therefore, the FPMD method is more favorable than the traditional
standard method for the study of surface problems.
We focus on the Si(001) surface system among a variety of surface systems, and study the clean and alkali adsorbed surfaces. A Si atom in theideal (001) surface has two dangling bonds as shown in Fig. 1.1(a). The
existence of dangling bonds unstabilizes the surface structure and a dimer
structure is formed by changing the atomic positions to reduce the number of dangling bonds. There are two kinds of dimer structures: one is asymmetric dimer as shown in Fig. 1.1(b) and the other is an asymmetricdimer structure. Figure 1.1(c) shows the p(222) structure, one of the possible arrangements of the asymmetric dimers. In STM experiments [8, 9],
the symmetric dimer structure is observed. On the contrary, in LEED
(Low-Energy Electron Di raction) [10] and ARUPS (Angle-Resolved Ultraviolet Photoelectron Spectroscopy) [11] experiments, the asymmetricdimer structure is observed.
On the other hand, in most of the theoretical electronic structure calculations, it is indicated that the asymmetric dimer structure is more stable
than the symmetric one [12, 13, 14, 15]. We also attempt to obtain the
optimized structure of the clean surface by allowing the ions to move automatically within the limitation of the (221) unit cell and calculate the4electronic structure of this system simultaneously. As a result, we ndthat the asymmetric dimer is more favorable in three di erent conditions(a cuto energy, k-points { [6.25Ry, 8k], [6.25Ry, 32k], [9.0Ry, 32k]). In
some other theoretical groups [13, 14, 15, 16, 17, 18, 19], they have obtained the results that the antiferro-type asymmetric dimer c(422) [15] orp(222) [17, 18, 19] structure is the most stable one by using a larger unit
cell [222, 422, etc.] than ours.
We shall extend our study to the alkali adsorbed Si(001)-221 and 223
surface systems. The kinds of adsorbed atoms are Li and Na. Several theoretical and experimental works on the alkali adsorbed Si(001) surface areavailable. Generally speaking, adsorption will be de nitely accompanied
by a structural relaxation of the substrate, which will modify the energetics of the adsorption. In this sense, theoretical works are not available so
far which can give reliable answers to the actual adsorption sites among
various possibilities. Therefore, the FPMD method is very suitable to thestudy of these cases.
Some fundamental aspects are controversial. For example, the saturation
coverage of Na or K has not been well de ned in the absolute scale: controversies are still going on about the single-layer model (like the Levine
model) [20, 21, 22, 23, 24] versus the double-layer model [25, 26, 27]. The
optimum adsorption sites have not been obtained yet. The alkali adsorbed
surface structure for the 223 (or 322) unit cell [23] is not known. It isalso interesting that the adsorption site of alkali adatoms depends strongly
on the coverage and the kind of alkali atoms. With the aim of obtaining aconsistent understanding of all these problems, we shall try in this thesis adetermination of optimum structures with full relaxation of ions and study
the dependences of the structures on the coverage and the kind of alkali
atoms.
We obtain the optimized structures and electronic properties for thesealkali adsorbed surface systems by using the FPMD method. The important results are summarized as follows. In the case of Na adsorption wetreat three coverage types: the rst one is a full coverage (two alkali atoms
per a 221 surface unit cell [=1] and it is called a `double-layer model'),the second one is a half coverage (one alkali atom per a 221 surface unit
cell [ = 0.5] and is called a `single-layer model'), and the third one is a
 = 61 coverage (one alkali atom per a 223 surface unit cell).
5For the full coverage Na case, the most stable structure is the combination of the hollow (HH) site along the dimer chain and trough (valley
bridge, T3) site on-top of the third layer Si atom. The optimized substrate
structure turns out to be a nearly symmetric dimer with a signi cant movement of the substrate atoms toward the ideal position. For the half coverage
Na case, the adsorption at the cave (T4) site is slightly more stable thanthe hollow (HH) one. In this case, the asymmetric dimer structure does
not change to the symmetric one unlike the full coverage case, but thedegree of asymmetry is reduced compared with the initial one of the cleansurface.
Next, we consider the Si(001)-223 surface system with  = 16 . In thiscase, the Na atom is adsorbed on the bridge (B2) site between the tworaised Si atoms of the neighboring Si-Si dimers. This adsorption site of Na
is consistent with the STM results for K and Cs [28, 29, 30].
In the case of Li adsorption, we treated two coverage types, one is afull coverage, the other is a  = 61 coverage. For the full coverage Licase, Li atoms are adsorbed at low symmetry positions. In this case, there
are two energetically almost degenerate structures where the two possible
adsorption sites along the trough are signi cantly o the valley bridge(T3) and the cave (T4) sites. Accordingly, the adsorption site of Li along
the dimer chain is also slightly o the hollow (HH) site. In the  = 16
coverage case, the result is almost similar to the Na adsorption case. This
is inconsistent with the STM result [28, 29, 30]. The following two reasonsare conceivable for the inconsistency: rst, the present coverage  = 61 ismuch higher than the experimental one (only a few percent); second, thenumber of the plane wave basis functions is not sucient to treat the deep
pseudopotential of Li. An improvement of the calculation in either of thetwo aspects requires too much of the computer resource and is left as a
future task.
In most of the present calculations, we assume that the charge density of
the core states is spatially separated from that of the valence states. This
allows us to use a simple treatment of the exchange-correlation potentialand energy, which works fairly well for Si. However, in the case of Na, theinteratomic distance is underestimated by this simpli cation. Therefore,
we implemented a partial core correction [31] for  = 1 and 0.5 in orderto improve the calculation. We nd that the calculated Na-Si bond length
6becomes closer to the experimental estimation [24].
The plan of the present thesis is as follows. In Chapter II, we will describe
the method of calculation about 1)pseudopotential, 2)total energy calculation by using the standard matrix diagonalization, 3)Hellmann-Feynman
forces, 4)the rst principles molecular dynamics method(Car-Parrinello
method), and 5)molecular dynamics method for the ionic system. The
results and discussion are presented in Chapter III. The results are shownfor 1)test calculations for bulk systems, 2)the clean Si(001)-221 surface
for symmetric and asymmetric dimer cases, and 3)alkali adsorption on theSi(001)-221, 223 surfaces. Chapter IV is devoted to concluding remarksof this thesis.

7Chapter 2

Method of calculation
2.1

Pseudopotential

2.1.1 Norm-conserving pseudopotential

The basic concept of pseudopotential is a frozen core approximation
which assumes that the electronic states of core electrons are insensitive tothe neighboring atomic con guration. The pseudopotential is constructedfrom the atomic potential so that the lowest eigenstate of the pseudopotential may correspond to the lowest valence state of the atomic potential. As
a consequence, the pseudopotential does not have a r1 -like part in the vicinity of r = 0 unlike an atomic potential. It is possible to use the plane wave
expansion in the electronic structure calculation because of the smoothed
potential in the core region.
We would like to make some comments on the use of the pseudopotential
and the plane wave expansion. As for the total energy per atom, ve
e ective digits are sucient in most cases in the pseudopotential approach
because of the absence of the core contribution, while at least eight e ective
digits are required in the all-electron calculation. This nice aspect of thepseudopotential approach is obscured to some extent by the plane wave
expansion because of very slow convergence in the total energy with respect
to the number of plane wave basis functions. However, we should notethat only the di erence in the total energy between di erent situations ismeaningful and that the convergence is much faster in the total energy
di erence than in the total energy itself. As force is essentially a totalenergy di erence, relatively good convergence may be expected in the force
calculation. Except such a convergence problem, forces acting on ions caneasily be calculated with the plane wave expansion by using the Hellmann8

Feynman theorem [3, 4, 5]
Traditionally, empirical pseudopotentials constructed by using some experimental data have been adopted [32, 33, 34, 35, 36]. Such empiricalpseudopotentials, however, have a weak point that the charge density does
not coincide with that of the real atom even outside the core region r > rc.
In 1979, Hamann, Schluter and Chiang proposed a new ab initio pseudopotential [37]. A slight modi cation was made and a complete tabulationof pseudopotentials was given by Bachelet, Hamann and Schluter (B-H-S)
in 1982 [38]. The following conditions are imposed in this new pseudopotential:
1. Real and pseudo valence eigen-values agree for a chosen atomic con guration.
2. In the region of r > rc, the amplitude of the pseudo-wave function
coincides with that of the real atomic-wave function.
3. The pseudo-wave function does not have a node, so the plane wave
basis can be used in the electronic structure calculation.
4. In the region of r < rc, the norm of the pseudo-wave function coincides
with that of the real wave function.
5. The rst order energy derivatives of the logarithmic derivative of thepseudo- and atomic-wave functions coincide with each other at rc.
Z r
1d @R(r) 12
2
0 2 r R(r) dE ( @r R(r) )jr=rc  0 c R3(r)R(r)r2dr;(2.1)
where R is a radial wave function. The Hartree atomic units (e2 =h = me = 1) are adopted throughout in this thesis. The cuto radius
rc is de ned as rc = rmax
cc with rmax is an outermost peak of the wave
function and cc is an adjustable parameter in the range of 1.5  2.0
as seen in Table I in the paper of B-H-S [38].
This B-H-S pseudopotential has a desirable feature that the charge density in the region of r > rc is suciently accurate. It has some weak points.
The plane wave expansion cannot be applied in practice to some materials
which have a strong core electron e ect. For example, B, C, N, O and mostof transition and noble metals have deep p or d potentials, which requires
a huge number of plane waves in the electronic structure calculation.
9

The following is the process of constructing the B-H-S pseudopotential.
Step1) The rst step pseudopotential V1(r) is made from an atomic
potential which is calculated by considering all electrons.
with

r!
r!
V1 (r) = Vatom(r)[1 0 fcut
] + Cfcut r ;
rc
c

(2.2)

r ! 0( rrc )
=e ;
fcut
(2.3)
rc
where Vatom(r) is an atomic potential, fcut(r) a cuto function, C an adjustable parameter, rc a core radius (cuto radius), and  = 3.5. The cuto

function should satisfy the conditions

r!
lim
f
= 1;
r!0 cut rc

and

(2.4)

r!
(2.5)
rlim
!1 fcut rc = 0:
The adjustable parameter C is determined so that the eigen-energies forthe pseudopotential coincide with those of valence states for the original
atomic potential.
Step 2) The rst step pseudo-wave function W1(r) is obtained by solving
the Schrodinger equation with V1(r).
Step 3) The second step pseudo-wave function is calculated to give thesame amplitude of the atomic-wave function in the region of r > rc.
W2 (r) =!

[W1(r) +  1 r(l+1)fcut r ];(2.6)
where and  are adjustable parameters which are determined to satisfy
the norm conservation.
Step 4) The second step pseudopotential is obtained by inverting theSchrodinger equation as follows,
V2 (r) = V1 (r) +rc

 1 r(l+1)fcut ( rrc )2W2(r)10
2 ( rrc )2 0 [2l + ( + 1)]( rrc )
+ 2 0 2V1(r) :
r2

2
(2.7)
Step 5) Vps(r) is obtained by subtracting the Hartree potential VH (r)and exchange-correlation potential Vxc(r) due to the valence electrons fromV2 (r),
Vps(r) = V2 (r) 0 VH (r) 0 Vxc(r):(2.8)
Step 6) The Vps(r) is divided into two parts, one is a local part Vlocal(r)land the other is a non-local part Vnon
0local (r) dependent on the angularmomentum l:
Vps (r) = Vlocal(r) +l
j l > Vnon
(2.9)
0local (r ) < l j:lThis Vps(r) is the nal form of the B-H-S pseudopotential. We shall use
this pseudopotential in the present FPMD. There are some calculations inwhich only the local part is used [39]. However the local part is generally
very deep, so that more plane waves are required. Further, the neglect of
l-dependent non-local parts destroys the norm conservation.
X

2.1.2 Kleinman-Bylander separable form

In this subsection, we focus on the non-local pseudopotential term
( ). Its matrix element islVnon
0local rXZli(~k+G~ )~rd~=e0i(~k+G~ )~rVnon
r0local (~r)elZ= X 4 (2l + 1)Pl(cos ~k+G;~ ~k+G~ ) 01 jl (j ~k + G~ j r)jl(j ~k + G~ 0 j r)ll
2
2 Vnon(2.10)
0local (r )r dr;
where Pl(cos ~k+G;~ ~k+G~ ) is a Legendre function, ~k+G;~ ~k+G~ is an angle between vectors ~k + G~ and ~k + G~ 0, and jl(x) is a spherical Bessel function.
The right hand side of eq. (2.10) is approximated by a separable form byusing the Kleinman and Bylander (K-B) method [40]. The detailed formulation will be given in Appendix A. Equation (2.10) is transformed into

V~kB+0G;~H~k0+SG~ 0

0

0

0

0

11

H 0S
V~kB+0G;~ ~k +G~ 0)

B
V~kK+0G;~ ~k +G~ 0

KB (G~ )S~KB (G~ 0 )
S~k;l
k;l
;
= 4(2l + 1)Pl(cos ~k+G;~ ~k+G~ 0 )WllX(2.11)

where l(r) is a pseudo-wave function and
Z 1l2~ j r)l (r)Vnon
SlKB (~k + G~ ) = 0 jl (j ~k + G
0local (r)r dr;
Z 1

l2
Wl = 0 3l (r)Vnon
0local (r)l (r )r dr:(2.13)P0 (cos ~k+G;~ ~k+G~ 0 ) = 1;

(2.14)In the s and p potentials, Legendre functions areand~k + G~ ) 1 (~k + G~ 0)(
P1 (cos ~k+G;:~ ~k+G~ ) = ~j k + G~ j j ~k + G~ 0 j
Therefore, the nal form is given by
NL (~NL (q~0)Vl=0q)Vl=0
K
0
B
V
j =
;
0

with

l=0

~q;q~0W0

KB
VlNL=0 (~q) = (4) 2 S~k;0 (~q);
1

and
with

(2.12)

j

V~qK;~qc0;q~B0;q~0c l=1

VlNL(~q; q~c)VlNL(
q~0 ; q~0 c)
=1
=1
=
;
W1

(2.15)(2.16)(2.17)(2.18)

1 qc
2
S~k;KB(~q);(2.19)
VlNL
(
~
q
;
q
)
=
(12

)
c
=1
q 1
where ~q = ~k + G~ , q =j ~q j, and qc is the c(= x; y; z ) component of ~q.

In the actual calculations, the following simpli cation is assumedl2
Vnon
0local (r ) = Vnon0local (r )12

for l  3:

(2.20)

This approximation may be allowed for Li, Na and Si for which mixing ofstates with l > 1 into the valence states is small. Then eq.(2.9) is rewrittenas
Vps (r) = Vlocal(r) +Xll
j l > Vnon
0local (r ) < l j

1
X

1
X

= Vlocal(r) + j
j l >< l j
non0local (r)
l=0
l=2
2
= fVlocal (r) + Vnon
0local (r)g
1
Xl
2
+ j l > fVnon(2.21)
0local (r ) 0 Vnon0local (r )g < l j;
l=0where use has been made of the relation P1l=2 j l >< l j= 1 0 P1l=0 j l >< l j.
2
Therefore, in our practical calculations, we take ( Vlocal (r) + Vnon
0local (r ) )2
as the local part of the pseudopotential and ( Vlocal (r) 0 Vnon0local (r) ) as
the non-local one. The results of test calculations for bulk systems of Li,
Na, Al and Si with these pseudopotentials will be mentioned in detail in
Chapter 3.
Other types of separable forms have also been proposed for the non-local
part of the pseudopotential. For example, corresponding to our VlNL, Hohl
et al. [41, 42, 43] proposed
Z
0;s
2
l=0:
j0 (j ~k + G~ j r)Vnon(2.22)
0local (r)Ti;~k (r )r dr;
whereTi;~k =

j0 (j ~k + G~ 0 j r)9i;~k (G~ 0 );(2.23)

1;s2
j1 (j ~k + G~ j r)Vnon
0local (r )Ui;~k (r )r dr;(2.24)

and
where

l=1:

Z

( ) <j +V 2

ll > Vnon0local r

X

G~ 0~k + G~ 0 )(Ui;~k =
j1 (j ~k + G~ 0 j r)9i;~k (G~ 0 ):(2.25)
0
~
~
j
k
+
G
j
~
G
These formulae require an integration of the radial direction in the real
space in every iteration, while such a radial part integration is performed
only at the initial setting in the K-B method.
X
0

13

2.2

Total energy calculation

We describe the total energy calculation for an electronic system byusing the norm-conserving pseudopotential and the plane wave basis setsbased on the local density approximation (LDA) [44, 45, 46] in the density
functional (DF) theory [47, 48].
At the rst step, we describe the DF theory and LDA. The DF theory[47] is one of the variational methods with the one-electron density n(~r) as
the variational parameter. The ground state energy for the electronic partis expressed asZ

E [n(~r)] = F [n(~r)] + Vext(~r)n(~r)d~r;(2.26)
where Vext(~r) is an external potential (a potential from nuclei) and F isa universal functional which does not depend on Vext(~r). The externalpotentials are the local and non-local pseudopotentials in this study. The
universal functional F is expressed as
Z Z n(~
1r)n(r~0 ) ~0
F [n(~r)] = T [n(~r)] +(2.27)
2 j ~r 0 r~0 j d~rdr + Exc[n(~r)];where T [n(~r)] is a kinetic energy in a non-interacting electron system, thesecond term is a classical Coulomb interaction and Exc[n(~r)] is an exchangecorrelation energy. The one-electron density n(~r) and the kinetic energy
are obtained by the following Kohn-Sham (K-S) equation [48]:[0 21 1 + Veff (~r)]9i(~r) = i9i(~r);(2.28)Zn(r~0) ~0 Exc [n(~r)]
Veff (~r) = Vext(~r) +dr +;(2.29)
n(~r)
j ~r 0 r~0 j
M ZX
T [n(~r)] =93i (~r)(0 21 1)9i(~r)d~r;(2.30)
i=1
andn(~r) =M
Xi

j 9i(~r) j2;

14(2.31)

where i speci es the one-electron eigenstate, 9i(~r) is its wave function and
M is a number of occupied states. The normalization constraint is imposedon the wave function as follows
Z

j 93i (~r) j2 d~r; = 1:
(2.32)As the functional of Exc[n(~r)] is not known in general, the LDA is commonly used which assumes the following expression.
Z

Exc [n(~r)] = xc (n(~r))n(~r)d~r;(2.33)
where xc(n(~r)) is called the exchange-correlation energy density and is afunction of the value of the one-electron density n(~r) at every ~r. As forthe actual functional form of xc(n(~r)), the results for the homogeneouselectron gas with the density n(~r) is utilized. Among many variations[44, 45, 49], the result by Ceperley and Alder [50] as parametrized by
Perdew and Zunger [46] is widely used recently. We adopt in this thesis
Wigner's interpolation formula [51].
We will express the total energy in the momentum space [52]. The wave
function is expanded in plane waves
9i;~k (r) = X 9i;~k (G~ )ei(~k+G~ )~r;~G

(2.34)

where i is a band index, ~k a wave vector in the rst Brillouin zone, G~ areciprocal lattice vector and 9i;~k (G~ ) a coecient of the expanded planewave. The K-S equation, eq. (2.28), in the momentum space is written as
X 1[ 2 (~k + G~ 0)2G;~ G~ + V~k+G;~ ~k+G~ ]9i;~k (G~ 0) = i;~k 9i;~k (G~ );(2.35)
~
G
with
0

0

0

~ 0 G~ 0) + Vxc (G~ 0 G~ 0 ) + Vlocal(G~ 0 G~ 0 )V~k+G;~ ~k+G~ 0 = VH (Gl~ ~k + G~ 0)P^l ;+ X Vnon0local (~k + G;l

(2.36)

where i;~k is the eigen-value and P^l is a projection operator for l.
The eigen-vectors and eigen-values are obtained by diagonalizing thematrix of the K-S equation, eq.(2.35). However, the matrix diagonalization15

method requires an N 3 order task of calculation and an N 2 order memoryin the case of N plane waves. It becomes virtually impossible to calculate
the electronic structure in the case of a large unit cell (supercell). For
example, in the case of double precision complex with N = 10000, thememory region amounts to an order of (N 2 2 16 =) 1.6Gbytes and thenumber of numerical operations an order of (N 3 =) 1012.
The charge density is needed for obtaining the Hartree and exchangecorrelation potentials or energies. As the task for obtaining the charge
density is very heavy in our calculation, we have to adopt an ecient
procedure of calculating the charge density. The charge densities in thereal and momentum spaces are(~r) =

occupied
X X

with
i;~k (~r) =j

and in G~ space

i;~k (~r);(2.37)

9i;~k (G~ )eiG~~ r j2;(2.38)

i

X

G~

~k

1 (r)e0iG~ 1~rd~r =(G~ ) =Z

with

occupied
X X
i

~k

~ );
i;~k (GZ1~r
i;~k (G) =i;~k (~r)e0iG~d~r= X 93i;~k (G~ 0 0 G~ )9i;~k (G~ 0);G~ 0(2.39)(2.40)

where the integration is over the unit cell volume .
We rst obtain 9i;~k (G~ ) as eigen-vectors of the eigen-value problem andthen calculate i;~k (G~ ), (G~ ), i;~k (~r) and (~r). If we adopt eq. (2.40) for
calculating i;~k (G~ ), preparation of the set f(G~ )g requires an N 2 ordertask of calculations. This can be reduced to an N log N order task by
calculating i;~k (~r) of eq.(2.38) by using inverse FFT and then i;~k (G~ ) again
by using FFT as
16

1  (~r)e0iG~~ rd~r:(2.41)
i;~k (G~ ) =i;~kThis procedure has another merit that (~r) is obtained as a byproduct of(G~ ). It should be noted that even if a plane wave expansion is adopted,
(~r) is a necessary quantity in the calculation of the exchange-correlationZ

potential and energy.
Finally, the expression for the total energy in our calculation is3(n01)(G~ ) (n) ~i;~k 0 4 (G)total =
2
~
G
i
~k
~
G6=0
3(n) ~
0 X Vxc3(n01)(G~ )(n)(G~ ) + 4 21 X  ~ (2G) (n)(G~ )
G
~
G
G~6=0
+ X Exc3(n)(G~ )(n)(G~ ) + 1Z + Ewald;(2.42)E (n)occupied
X X

X

~
G

where the superscript n implies the n-th step iteration, 1Z is the correction term of a local pseudopotential [53] and Ewald is the Ewald term
which is an electrostatic energy ofPoccupied
ions and the uniform part of the electron
density. In the above equation, imeans that a summation over theenergy bands at each k-point for sampling are up to the Fermi energy.
Although the last term of the total energy equation includes the diverging terms for the ionic and electronic parts, it can be properly evaluatedby the Ewald method. The expression for Ewald is given byX
1
4 X0 1 eiG~ (~s0~s )e0 G~ 0 1 
=
2
Z
Z
Ewald
2 s;s s s G~ j G~ j2
~+ X0 1 0 erf~( j l + ~s 0 ~s j) 0 2(  ) s;s ; (2.43)
j l + ~s 0 ~s j~lwhere ~l is a lattice vector, Zs the valence of the s-th ion, ~ s its basis
vector in the unit cell, a numerical factor for convergence in the EwaldP
summation and erf (x) an error function. In this equation, the prime ofmeans that it excludes the term of G~ = 0 and also the one of ~l = 0 whens = s0 .
The Fermi energy must be obtained before calculating the total energy
and charge density. The summation over k-points for sampling in the rst
0

0

j j2
4

0

1
2

0

0

17

1
2

0

Brillouin zone is usually performed by the tetrahedron method [54, 55] for
accurate calculations. However, in supercell calculations, the number ofsampling k-points is only a few in general so that the tetrahedron method
is not applicable. Therefore, in our calculation the Fermi energy f is
obtained by summing up the arti cially broadened energy spectrum: theenergy spectrum at each k-point is broadened with an adequate width. The
width is not necessary for non-metallic systems with a band gap. However,if the width is very small for metallic systems, the Fermi energy oscillates
in the process of calculation because of the possible interchange between
the occupied and unoccupied states.
After the eigen-vectors and eigen-values are obtained by using the matrix diagonalization, the Fermi energy is determined by counting the occupied states. The charge density is obtained from the eigen-vectors. Afterthese calculations, the total energy of the electronic part is calculated for agiven ionic con guration. The nal step is a convergence check and charge
mixing. The input and output charge densities are mixed for faster convergence. The calculation is continued until the di erence between input andoutput values (charge density, potential, total energy, etc) is smaller thana given criterion. Here, the ow chart of the above calculation is shown in
Fig. 2.1.
2.3

Hellmann-Feynman force

It is possible to calculate the forces [52] which act on the ions in the unit
cell on the basis of the Hellmann-Feynman theorem [3, 4, 5]. The force
acting on the s-th ion is divided into two contributions:@E
F~s = 0 total = F~se + F~sEwald:@~s(2.44)The electronic part F~se is given by
F~se = i

0 i

X~ s s~ )Ge~ 0iG~
3 (GVlocal (G~ )~G
occupied
X XX Xi~kl G;~ G~ 0(G~ 0 0 G~ )93i;~k (G~ )9i;~k (G~ 0)ei(G~ 0G~ )~s0l;s
~ ~k + G~ 0);
2 Vnon0local (~k + G;

18

(2.45)

and the ionic part F~sEwald originating from Ewald is given byX
4 X [ G~ sin [G~ (~ 0 ~ 0 )]e0 jG4~ j2 ]
F~sEwald = Zs Zs0
s
s
~ j2
s6=s0
~ 6=0 j G
G
1
X ~
1 2~
x 0 j~xj2 
xerfc( 2 j ~x j)
2++ (  ) j ~x j2 e
; (2.46)~x=~l+~s0~s0
j ~x j3~lwhere erfc(x) = 1 0 erf (x). In the above equations based on the plane
wave basis, it is interesting to note that several contributions originatingfrom the ~s dependence of the electron density (~r) cancel out and thatl
~ ~k + G~ 0) remain.
only the terms including Vlocal(G~ ) and Vnon
0local (~k + G;
In F~se, the second term in the right hand side of eq. (2.45) is a double
summation with respect to G~ and G~ 0. Therefore, the task of calculationof this term is very heavy. As is shown below, however, it is possible toreduce this task by using the K-B separable form of the non-local part ofthe pseudopotential [40]. In the second term of eq. (2.45), we consider thefollowing quantitye;NL
F~s;i;~k = 0i

X~ G~ 0G;(G~ 0 G~ 0)93i;~k (G~ )9i;~k (G~ 0)ei(G~ 0G~ )~s0l;s
~ ~k + G~ 0 ):2 Vnon0local (~k + G;
By using eqs. (2.17) and (2.19), we de neand(2.47)X~ s NL ~A~ s;i;~k = G~ 9i;~k (G~ )eiG~
Vs;l (k + G~ );(2.48)X~ s NL ~B~ s;i;~k = 93i;~k (G~ )eiG~
Vs;l (k + G~ );(2.49)~
G~
G

with which eq. (2.47) is reduced tos;NL~ ~~3 ~ 3F~e;i;~k = 0i (As Bs 0 As Bs ):(2.50)The task of the calculation for the non-local pseudopotential part in theforce is very heavy with the non-separable form. With the increase of the19number of atoms and plane waves, this task becomes heavier and heavier
so that it becomes dicult to perform the molecular dynamics for ions byusing these forces. Therefore, the reduction of the task from an N 2 orderto an N order by using the K-B separable form is crucial in practice.
2.4

Car-Parrinello method

So far, we have described the standard electronic structure calculation byusing the norm-conserving pseudopotential, the plane wave basis and thematrix diagonalization on the basis of the local density approximation in
the density functional theory. The problems in the matrix diagonalization
method are an N 3 order task of calculation and an N 2 order memory.
In 1985, Car and Parrinello developed a new and epochal method forthe electronic and ionic structure calculation [6]. In this method, the electronic and ionic parts of the system are solved simultaneously by using thedi erential equations like a Newton equation. It is possible to obtain the
ground state of the electronic structure and optimize the ionic structure
simultaneously. In addition, this method does not include matrix diagonalization and the task of the calculation and the memory are reduced. For
these reasons, it is advantageous to use this method in order to optimize
electronic and ionic structures in large systems.
In the optimization process for the ionic structure, the electronic stateshould be close to the Born-Oppenheimer surface of the ground state. Otherwise the accurate forces are not obtained and therefore, it is impossible
to perform accurate molecular dynamics for ionic systems. In the CarParrinello method, it is possible to keep the electronic state close to theBorn-Oppenheimer surface by adjusting some parameters so that the electronic degrees of freedom may follow the ionic motion quickly.
In the following, we describe the basic idea of the Car-Parrinello methodbrie y. We express the total energy of the system as E [f9i;~k g; fR~ I g]. f9i;~k gdenotes a set of trial Kohn-Sham orbitals occupied by electrons and fR~ I g
is a set of ionic positions. The ground state of the system for the optimized
ionic positions corresponds to the global minimum of E [f9i;~k g; fR~ I g] in theextended space spanned by f9i;~k g and fR~ I g. In this sense, the structural
optimization based on the DF electronic structure calculation corresponds
to the minimization of E [f9i;~k g; fR~ I g] with respect to f9i;~k g and fR~ I g.
20

Car and Parrinello adopted a molecular dynamics method for the minimization of E [f9i;~k g; fR~ I g] and introduced a ctitious equation of motion
for 9i;~k through the following Lagrangian.
1 M R~_ 2 0 E [f9 g; fR~ g]
I I
I
i;~k
I 2
i;~k
Z
+ X 3ij;~k ( 93i;~k (~r)9j;~k (~r)d~r 0 ij );(2.51)L =XZ

m j 9_ i;~k (~r) j2 d~r +Xi;j;~kwhere m is a ctitious mass, MI is the mass of an ion and the last term with
the Lagrange multiplier 3ij;~k originates from the orthonormality constraints
Z93i;~k (~r)9j;~k (~r)d~r = ij :With eq. (2.51), the equations of motion are given by
0

1

0

1

d B L C LA=@;
dt 9_ 3i;~k93i;~kandd @ @L A @L
= @R :
dt @ R_ I
I

Therefore i;~k (~r; t) = 0
m9
and

X
E
+
3 9 (~r; t);93i;~k (~r; t) j ij;~k j;~k

MI R~ I = 0 5R~ I E:

(2.52)(2.53)(2.54)(2.55)(2.56)

By using E9 = H 9, eq. (2.55) is reduced to i;~k (~r; t) = 0H 9i;~k (~r; t) + X 3ij;~k 9j;~k (~r; t);m9(2.57)jwhere H is the Kohn-Sham one-electron Hamiltonian. Equation (2.57) issolved numerically by the Verlet algorithm [56] combined with the Ryckaert
et al. method [57] for determine 3ij;~k . It should be noted here that the timedependence of the wave function does not have any real physical meaning.
321

It is completely a mathematical trick for solving the eigen-value problem
for the Kohn-Sham equation iteratively. In fact, if 9i;~k comes close to thesolution of the Kohn-Sham equation, 9i;~k becomes stationary and eq.(2.57)
is reduced to the K-S equation.
As for the optimization process for the structure of ions, the systemgenerally has some local minimum states. In the case of searching the
local or global minimum states of E [f9i;~k g; fR~ g], the rst step of the optimization process is to cooling down only the electronic state to the BornOppenheimer surface (ground state) for the initial xed ionic position.
This is shown as a process (1) in Fig. 2.2. After this step, the electronic
and ionic states are solved simultaneously as a process (2). Finally, theprocess of calculation will stop when the system comes suciently close to
the local or global minimum state as a process (3).
There are some methods which are di erent from the original CarParrinello method. The steepest descent (SD) and the modi ed steepest
descent (MSD) methods [58, 59, 60, 61] which originate from Williams
and Soler [62] are used in this study and will be described in detail in thefollowing subsections.
2.4.1 Steepest descent method

We adopted two methods, one is the steepest descent (SD) method which
will be described in this subsection, the other is the modi ed steepest
descent (MSD) method which will be described in the next subsection. In
the SD method, the eigen-vectors are solved dynamically by integrating
the rst order di erential equation.
~ t) = 0 1 j ~k + G~ j2 9i;~k (G;~ t) 9_ i;~k (G;2
0 XfVH (G~ 0 G~ 0) + Vxc(G~ 0 G~ 0) + Vlocal (G~ 0 G~ 0)g9i;~k (G~ 0; t)

0

G0

X

G0 ;l

l
~ ~k + G~ 0)9i;~k (G~ 0 ; t)Vnon
0local (~k + G;~ t);
+ X 3ij;~k 9j;~k (G;(2.58)jwhere  is a ctitious friction parameter. The simplest way of solving
eq. (2.58) is22

~ t + 1t) = 9i;~k (G;~ t) + 1t9_ i;~k (G;~ t);9i;~k (G;(2.59)
where 1t is a time step. The explicit expression for 3ij;~k , which is valid tothe rst order in 1t, is given byZ

3ij;~k (t) = 93j;~k (~r; t)H (t)9i;~k (~r; t)d~r:(2.60)The derivation of eq.(2.60) is given in Appendix B. In this thesis, we replace
3ij;~k with a diagonal matrix i;~k ij . This means that only the normalization
constraint for the wave functions is taken into account. The Gram-Schmidt
orthogonalization process is introduced after every time step in order to~ t + 1t)g. In order to save comguarantee the orthonormality for f9i;~k (G;putational time, we further introduce a following approximation, bearing
the following situation in mind. The initial stage of the calculation is fairlyecient in the sense that the system approaches the equilibrium con guration rather quickly. Then the force acting on each atom becomes small andthe time evolution of the system becomes slow. The calculation for thispart requires a lot of computer time to attain the convergence criterion.
The slow time evolution of the system in the latter part of the time steps
allows us to use the following expression for i;~k :i;~k (t + 1t) =Z

93i;~k (~r; t + 1t)H (t)9i;~k (~r; t + 1t)d~r;
which can be transformed into an ecient expression as
i;~k (t + 1t) =(2.61)~ t + 1t) < G~ j H (t) j G~ 0 > 9i;~k (G~ 0; t + 1t)
93i;~k (G;~ G~G;~ t)9i;~k (G;~ t + 1t) 0 93 ~ (G;~ t + 1t)9_ i;~k (G;~ t)= i;~k 0  Xf9_ 3i;~k (G;i;kG~3~ t)9_ i;~k (G;~ t)g + O(1t2);+9i;~k (G;(2.62)
where < G~ j means the bra of plane waves. This expression is derived
in Appendix C. Although the rst expression of i;~k (t + 1t) in eq. (2.62)
includes the double summation over G~ and G~ 0, the second approximate expression includes only the single summation over G~ . Therefore, an N 2 ordertask is reduced to an N order one. As all the expressions from eq. (2.59)~ t) is deto eq. (2.62) are approximations, the time dependence of 9i;~k (G;scribed only approximately. Nevertheless, the starting di erential equation
X

0

23

~ t) is an arti cial one and thereby it is not necessary to solve it
for 9i;~k (G;very accurately. We should emphasize that the nal stationary solution is
not a ected by the approximations introduced above. The ow chart of
the SD (or MSD) calculation is shown in Fig. 2.3.
2.4.2 Modi ed steepest descent method

In the MSD method, the di erential equation in eq. (2.58) is solved
analytically on the assumption that only the diagonal part of Hamiltonian,~ >, depends on time [58, 59, 60, 63]. Therefore, the allowed< G~ j H (t) j G
1t is the time interval during which 9i;~k (G~ 0) (G~ 0 6= G~ ) can be regarded
as a constant. As a result, it is possible to perform the calculation with alarger time step 1t than that of the SD case, so that the number of steps
in the MSD method for attaining a convergence is 21  15 of that in the SD
method.
The explicit expression for the MSD is derived as follows,~ t) = 0 [ 1 j ~k + G~ j2 +V (G~ = 0)
 9_ i;~k (G;2X l~ ~k + G~ ) 0 i;~k ]9i;~k (G;~ t)+ Vnon0local (~k + G;ll
~ ~k + G~ 0)]
0 X [V (G~ 0 G~ 0) + X Vnon
0local (~k + G;~)
G~ 0(6=Gl29i;~k (G~ 0; t);(2.63)
where V (G~ 0 G~ 0) = VH (G~ 0 G~ 0) + Vxc(G~ 0 G~ 0) + Vlocal (G~ 0 G~ 0), the rst
term is a diagonal part of Hamiltonian with the time dependence taken
into account and the second term is an o -diagonal part which is assumed
to be time independent. By de ning the following quantities W and R by
1  1 j ~k + G~ j2 +V (G~ = 0)
W =
 2

X l
~
~
~
~
+ Vnon0local (k + G; k + G) 0 i;~k ;(2.64)land
1 X [V (G~ 0 G~ 0) + X V l~ ~k + G~ 0)](~k + G;R =
 G~ 0 (6=G~ )l

24

non0local

29i;~k (G~ 0; t);eq. (2.63) becomes(2.65)~ t) = 0W 9i;~k (G;~ t) 0 R
9_ i;~k (G;~ t) + R ):= 0W (9i;~k (G;W
The solution of eq. (2.66) is(2.66)~ t + 1t) = e0W 1t9i;~k (G;~ t) + (e0W 1t 0 1) R :9i;~k (G;(2.67)W
In this treatment, the error originates from the assumption that R is timeindependent. From our experience, the time step 1t for the MSD calculation can be ve times as large as that of the SD case. If we use an evenlarger time step, the calculation does not converge.
In the original Car-Parrinello method, the sum of the kinetic energyof the electronic and ionic parts and the total energy of the electronicpart which includes the electrostatic term of ions (Ewald term) is constant exactly in the process of calculation if the Ryckaert method [57] fororthonormalization of eigen-vectors is adopted:

1 m Z j 9_ j2 d~r + X 1 M R~_ 2 + E [f9 g; fR~ g] = const: (2.68)
I I
tot
I
i;~k
i;~k
2
2
I
~
i;k
However, in the SD and MSD methods, such an energy conservation does
not hold. The terms in the left hand side in eqs. (2.58) and (2.63) havea friction parameter  , so that an energy dissipation should occur in theprocess of calculation.
This energy dissipation is explained in the following. Let us calculate
the following quantity:
d 1 X ~_ 2 dEtot
f M R g + dt :(2.69)
dt 2 I I I
The derivative of the total energy is
occupiedX

occupiedX X Z Etot d9i;~kEtot d93i;~kdEtot=(  9 dt + 93 dt )d~r
dt
i
~k
i;~k
i;~k

25

@Etot dR~ I~ I dt :I @R+(2.70)In this equation, the rst and second terms are reduced toEtot d9i;~k=(H 93i;~k )9_ i;~k ;(2.71)
9i;~k dt
Etot d93i;~k=(H 9i;~k )9_ 3i;~k ;(2.72)
3
9i;~k dt
and the last term is transformed into~I
X @Etot dRX
~ I R~_ I = 0 d ( 1 X MI j R~_ I j2):=0MR
(2.73)
I~ I dtdt 2 I
I @R
I
Therefore, the last term in eq.(2.73) cancels with the rst term in eq.(2.69).
The
remaining
terms
can
be
rewritten
by
using
the
equation 9_ i = 0H 9i +P
j 3ij 9j :XXi

< 9_ i j H j 9i > +X

i< 9_ i j H j 9i >3= f< 9_ i j Hij 0 3ij j 9j > + < 9_ i j 3ij j 9j >g + c:c:ij
= 02 X < 9_ i j 9_ i > + X 3ij f< 9_ i j 9j > + < 9i j 9_ j > g(;2.74)
X

i

ij

where the index ~k and the summation over it are implicitly included. Noticing that dtd < 9i j 9j >= 0 for the second term of eq. (2.74), we obtainoccupied=X Z
d < 1 X ~_ 2~
_ i;~k j2 d~r:
MI RI + Etot[f9i;~k g; fRI g]; = 02j9:
dt 2 Ii;~k8

9(2.75)

This relation suggests that the total energy of the ionic part is decreasing
monotonously.
In our calculation, the charge density is mixed in each iteration in orderto suppress the charge oscillation which seems to be enhanced in surface
systems. In our test calculations for some bulk systems, it is not necessary
to introduce this charge mixing. The e ect of this charge mixing on the26

energy balance problem discussed above is not suciently understood atpresent.
In our SD and MSD calculations, the initial wave function expanded inplane waves is given by using the matrix diagonalization. In this case, thecuto energy, which determines the number of plane waves, is set to beabout half of that in the following iterations in order to save memory andCPU-time. Of course, the nal results are not a ected by this initial wave
function.
2.5

Molecular dynamics for ionic part

The ionic motion is governed by the Hellmann-Feynman forces [3, 4, 5]
which were discussed in Section 2.3. The ionic coordinates in the unit cellare updated by the Verlet algorithm [56] as follows,
F~ (R~ )R~ (t + 1t) = 2R~ (t) 0 R~ (t 0 1t) + HF I (1t)2(2.76)I

I

I

MI

where R~ I (t + 1t) is the coordinate of I -th ion at time t + 1t, F~HF is the
Hellmann-Feynman force and 1t is a time step.
In our calculation, the ions are initially displaced from symmetry points
suciently in order to allow them to move toward a lower energy con guration. At rst, the ions are accelerated by the forces and temperature of
the ionic system increases. However, due to the relation given by eq.(2.75)
the sum of the total energy and kinetic energy of the ions starts to decrease
after some time. The temperature of the ionic part eventually decreases tozero Kelvin after taking some peaks. The process of calculation is terminated when the convergence criterion that forces acting on ions is less than
1:0 2 1003 Ry=a:u: is satis ed. A word of warning is necessary here. The
use of the SD and MSD method for the electronic part tends to quickly
reduce the kinetic energy of ions and therefore lead to a local minimum
instead of the global minimum. To obtain the global minimum con guration, it may be better to adopt the original Car-Parrinello method for theelectronic part.

27

Chapter 3

Results and discussion
3.1

Test calculations for bulk systems

In this section, we will describe the results of test calculations for bulk
systems in order to check the pseudopotentials and the program. The calculated systems are Li, Na, Al and Si for which some theoretical calculations
have already been performed with the same norm-conserving pseudopotential with the plane wave basis [64, 65, 66, 67]. We calculate the electronic
properties (the total energy and the dispersion curve) and obtain the equilibrium lattice constants and bulk moduli for these systems. The bulk
systems of Li and Na are metals and their crystal structures are bcc. Al
is also a metal and its crystal structure is fcc. Si is a semiconductor andtakes a diamond structure.
The conditions (the cuto energy Ecut and number of k-points for sampling Nk ) and results of the calculations are summarized in Table 3.1. The
results of Al and Si agree fairly well with the experimental data [68]. On
the other hand, the agreement is poor for Li and Na: the equilibrium lattice constants are underestimated and the bulk moduli are overestimatedaccordingly. It is known that a partial core correction [31] plays an important role in the case of alkali metals. As a matter of fact, the results are
improved signi cantly by considering the partial core correction for Na.
Most of the results presented in this thesis are those without the partial
core correction except some results presented in section 3.3.
Next, we attempt a test of the FPMD for the electronic and ionic parts.
The following simple model is set for a test calculation. We prepare a
cubic unit cell which includes two Al atoms. The Al atoms are initiallydisplaced along the [111] direction from the positions corresponding to the28bcc structure. We rst calculate only the electronic part to bring it to
the Born-Oppenheimer surface for the ground state corresponding to theinitial condition. After this procedure, the electronic and ionic parts are
solved simultaneously to optimize the structure. The displaced atoms start
to return to the equilibrium bcc positions in the process of calculation andthe nal equilibrium structure becomes the perfect bcc structure.
3.2

Calculation of Si(001) surface

The purpose of this thesis is the clari cation of the optimized structures of the alkali adsorbed Si(001) surfaces and their dependences on the
coverage and adsorbed atomic type.
The geometry of the possible adsorption sites for Na is schematicallyshown in Fig. 3.1. Here the following notations are used: HH and HBdenote the hollow site and the bridge site respectively along the hill of the
dimer chain, T3(T4) is the site along the trough on-top of the third(fourth)layer Si atom, B1 is the on-top site of the raised Si atom of an asymmetric
Si-Si dimer, and B2 denotes the bridge site between the two raised Si atomsof the neighboring Si-Si dimers. HH site is called a hollow or pedestal site,
HB is called a bridge site, T3 is called a valley bridge site, and T4 is calleda cave site by other research groups.
The starting substrate con guration is set to be the asymmetric dimer
for all systems, which is optimized separately for the clean surface. In eachcase, the alkali atoms are put slightly o (5%  12:5% of the side lengthof the unit cell) the symmetry positions HH, HB, T3, and T4 in order toallow them to move toward lower energy con gurations.
A repeated slab geometry and supercell method are adopted for surface
systems. A supercell is a large unit cell which includes several tens to ahundred of atoms. In a repeated slab geometry, the unit cell includes theatomic layers and vacuum region as schematically illustrated in Fig. 3.2.
The thickness of the vacuum region should be chosen in such a way thatthe interaction between the atomic slabs through the vacuum region canbe neglected and the number of the atomic layers should be chosen insuch a way that the interaction between two surfaces through bulk canbe neglected. In our calculation, the number of atomic (Si) layers is ten,and the thickness of the vacuum region is the same as that of the atomic
29

layers. The thickness of the vacuum region is sucient since the wave
functions for the occupied states decay into the vacuum exponentially andthe electron density near the center of the vacuum region is vanishingly
small. Figure 3.3 shows an example of the electron density for the fullcoverage Na adsorption case. On the other hand, the thickness of theatomic layers may not be sucient but is near the upper limit imposed by
the present computer capability.
The distributions of sampling k-points adopted in the following for the(221) unit cell are shown in Fig. 3.4. The cuto energy Ecut, the numberof plane waves Nw , the number of sampling k-points Nk , and some resultsof calculations are listed in Table 3.2. They include the total energies,
the heights of alkali atoms, the nearest interatomic distance between thealkali adatom and the substrate Si atom and the distance of the Si-Si dimer
atoms.
3.2.1 Clean surface

Some LEED experiments demonstrate that the clean surface of Si(001)exhibits a c(422) pattern which is generally believed to be an ordering
state of asymmetric dimers [10]. On the other hand, asymmetric dimers
can be clearly observed only near the step or some other defects by STM
measurements [8, 69]. In the clean surface region of Si(001) without stepor vacancies, STM observes a symmetric dimer. The symmetric dimers
observed by STM may be due to thermal oscillations of the asymmetricdimers, since it observes an average for time of the order of 1.021002
seconds [69].
We study whether the asymmetric dimer is stabilized or not with three
di erent conditions for the cuto energy and the k-space sampling and show
the results in Table 3.3. Although the energy di erence has not converged
yet with respect to the cuto energy Ecut, the asymmetric dimer has alower energy for all three cases, being consistent with other calculations[13, 14, 15, 16, 17, 18, 19]. The energy separation between the  and3 band gap states becomes wider by about 0.01 a.u. (0.27eV) when thesymmetric Si-Si dimer is transformed to the asymmetric one as seen in
Fig.3.5. The Si-Si dimer distance in our calculation is 2.30 A being shorter
than the experimental estimate of 2.45 A by 6%.
30

In Fig. 3.6, the electron density for the clean surface is shown. Its dimer
structure is asymmetric and the raised Si atom of the asymmetric dimeris slightly more populated by electrons than the lower one. This chargetransfer was rst discussed by Chadi [16] based on the empirical tightbinding calculation and also discussed by Zhu et al. [15] based on theconventional pseudopotential band calculation.
3.2.2 Full coverage Na/Si(001)-221 surface

The computational condition is such that Ecut = 6.25 Ry and Nk = 8
as seen in Table 3.2. As the system has an insulating gap for  = 1, thissmall Nk is expected to be sucient. We found two stable con gurations,the combinations HH-T3 and HH-T4. The former combination is lower
in energy by 28meV per dimer than the latter, being consistent with the
double layer model proposed for the K adsorption [26].
In contrast to the present result, Batra suggested that the T4 site isenergetically more stable than T3 [70]. He optimized the height of Na but
xed the substrate atomic con guration as the symmetric dimer model,
which was optimized separately for the clean surface. He also studied
the adsorption on the bulk terminated surface and concluded that suchan adsorption con guration would not be realized. We have con rmedthat the optimization of the Na position only gives us a result almost
identical to that of Batra. Recently, he found that the substrate relaxation
has a tendency of stabilizing the HH-T3 con guration to make it nearly
degenerate with HH-T4 [71].
The HB site of the Si-Si dimer is not a stable position for Na adsorption. The top views of the optimized atomic positions for the Na adsorption
with the full coverage in the HH-T3 and HH-T4 con gurations are shownin Fig. 3.7(a) and (b), respectively. Both for the HH-T3 and HH-T4 cases,a substantial relaxation of the Si substrate is induced (Fig.3.7(c)): the Si-Si
dimers are now symmetric and the dimer distance is 2.65 A (2.51 A) for
the HH-T3 (HH-T4) case. The details of geometrical parameters for the
optimized structures are listed in Tables 3.4 and 3.5. The total energy andtemperature of the ionic part in the process of calculation for the HH-T3
case are shown in Figs.3.8(a) and (b), respectively. In the rst hundred iteration steps, only the electronic part was cooled down to the ground state.
31

Subsequently, the electronic and ionic parts were solved simultaneously.
After iterations over 1500, the calculation attained a sucient convergenceand the optimized con guration was obtained. At rst, the temperatureof the ionic part increases rapidly and after several tens of steps, it begins
to decrease due to the friction e ect as described in Section 2.4.2. The
remaining small asymmetry in the positional parameters in the nal conguration may be some artifacts of our automatic optimization process.
As the atomic arrangement approaches the optimum one, the force acting
on each atom becomes vanishingly small (the criterion for stabilization is
1.0 2 1003 Ry/a.u.) so that atoms may not reach their exact optimum
positions. The disappearance of the asymmetry is a natural consequenceof the full occupancy of the 3 band-gap states, because the energy gain
caused by the lowering of the  band is cancelled out by the energy loss
caused by the counter movement of the 3 band as seen in Fig. 3.5.
Enta et al. [25] found that the surface with the saturation coverage of
K was insulating and suggested that the alkali atoms would be adsorbed
not only at the raised sites of the Levine model but also at the trough sites
between the dimer rows [26]. The saturation coverage was claimed to be
the full coverage  = 1.0. Recent experiments support the existence of two
di erent adsorption sites [27]. The Si(001)-221 surface with the saturation
coverage of Na was also shown to be insulating by Enta et al. [72]. However,some experiments claimed that the saturation coverage will not be  = 1
but  = 0.68 or  = 0.5 for the Na and K adsorption [22, 23, 73, 74]. A
realization of the adsorption pattern of the Levine model was claimed by
the LEED analysis [24].
These experiments were carried out at room temperature or at evenhigher temperature, whereas the present calculation is for zero Kelvin.
Therefore, a theoretical study on the e ect of the nite temperature isneeded to give a de nitive conclusion about the saturation coverage. Finally, we note that the shortest Na-Si distance in the HH-T3 con guration is 2.62 A (2.81 A) for Na at the HH (T3) site. This relatively short
Na-Si distance may imply that Na atoms are chemisorbed rather than physisorbed even for  = 1.
In Figs. 3.9(a) - (f), the contour maps of the electron density are shownfor the side and top views of the Na adsorbed Si(001) surface (HH-T3) for
the full coverage. The axes of x and y are shown in Fig. 3.1(a) by arrows,
32and the z axis is perpendicular to the Si(001) surface ([001] direction) as
shown in Fig. 3.1(b). The lled circles denote the positions of Na and Si
atoms projected on the cut plane. The Na atoms are indicated by arrows.
In the gure of the side view, the cut plane is shown in the attached small
gure. z(height) in Fig. 3.9(d) to (f) denotes a height of the cut plane forthe top view measured from the bottom of the unit cell. (With this origin
for z, the center of the substrate slab is located at z = 11.531 a.u.) The
heights of the adsorbed Na on the T3 and HH are 23.998, 25.329(a.u.),
respectively. In Fig. 3.9(a), the electron density in the dimer is symmetricand the covalent bond charge along the dimer, which is clearly seen in
Fig. 3.6 for the clean surface, does not exist. In Fig. 3.9(b), although thecut plane contains the Na atoms adsorbed at the HH site as indicatedby arrows, the electron density coming from the adsorbed Na cannot be
identi ed. Figures 3.9(c) and (d) show the contour maps of the di erenceelectron densities de ned as the total electron density of Na/Si(001) minus
the electron density of the Si substrate. In Fig. 3.9(c), it is found that thecovalent bond charge is formed between the Na adatom on the T3 site andthe second layer Si atoms. In Fig.3.9(d), which is a top view of the electron
density in the plane whose height is at the middle of the Na at HH andthe dimer Si, the formation of the bond charge between the Na and thedimer Si atoms is observed. Figures 3.9(e) and (f) are the contour map of
a top view and the bird's eye view for a plane at the height 30.329(a.u.)which is higher than that of the adsorbed Na on the HH site by 5.0(a.u.).
In the bird's eye view, the electron density coming from the adsorbed Naat the T3 site is slightly seen but that from the Na at the HH site cannot
be identi ed at all.
3.2.3 Half coverage Na/Si(001)-221 surface

Ecut = 9.0 Ry and Nk = 32 in this case. As the  3 band is half- lled for = 0.5, the system is metallic and a large Nk is required. The fundamental
question in the  = 0.5 case is whether the alkali atoms are adsorbed along

the hill of the dimer chain or along the trough between the dimer chains.
For Na adsorption, the pseudopotential calculation by Batra predicted thatthe HH site is more stable than T4 (T3) site by 0.12eV (0.22eV). On theother hand, a cluster calculation by Ling et al. [75] for K adsorption
33

predicted that the T4 site is more favored. In these calculations, however,the possibility of asymmetry in the Si-Si dimer was not taken into account.
The top views of the optimized atomic positions for the Na adsorption with the half coverage in the HH and T4 con gurations are shownin Fig. 3.10(a) and (b), respectively. As is shown in Figs. 3.10(c) and (d),the asymmetry of the Si-Si dimer remains for  = 0.5 in both cases of theadsorption along the hill and the trough of the dimer rows. There is still anenergy gain caused by the asymmetry-induced increase of the  0 3 band
separation as the 3 band is only half- lled. Even if the Na atoms are put
fairly close to the T3 sites, they move toward the T4 sites and stop there.
Along the hill of the dimer, the HH site is stable while the HB site is unstable. The present calculation shows that the T4 adsorption is lower thanthe HH adsorption by 12meV/dimer. Unfortunately, this energy di erenceis already close to the limit of the numerical accuracy. Analysis of experiments suggest that the preferential adsorption at T4 rather than HH may
not be plausible. It was shown by Ishida and Terakura [76] that the T4
site adsorption produces a much smaller lowering of the workfunction thanobserved. A recent LEED analysis concluded that the HH-site adsorption
is more favorable than the T4-site adsorption [24]. However, a caution may
be necessary to this LEED analysis. It was stated in the paper [24] thatthe double-layer model of the HH-T3 con guration is equally acceptable.
We calculate the binding energy of the Na adatom on the Si substrate inorder to obtain an insight to the saturation coverage. The binding energy
per one Na is de ned as follows. For the full coverage,
full
full
full
EBatra= 41 (ENa=Si
0 ESi0sub 0 ENa(3.1)
0lay );
full
full
full
0 41 [4ENa0atom 0 ENa
Ebind= EBatra(3.2)
0lay ];
and for the half coverage,
1 (E half 0 E
half
(3.3)
Si0sub 0 ENa0lay );
Na=Si
2
half
half
half
= EBatra0 21 [2ENa0atom 0 ENa
(3.4)
Ebind
0lay ];
full(half )is the binding energy adopted by Batra [71] between thewhere EBatrafull(half )Na layer and the Si substrate for the full (half) coverage, Ebind
ishalf
EBatra=34

thefulladsorptionenergy with respect to the assembly of isolated Na atoms,(
half )
ENa=Si is the total energy of the Na adsorbed Si surface, ESi0sub is thefull(half )total energy of the Si substrate, ENa
0lay is the total energy of isolated Nalayer corresponding to the full (half) coverage, and ENa0atom is the totalenergy of the Na (pseudo) atom. The factor 14 ( 12 ) takes account of thefact that the unit cell contains 4(2) Na atoms, half2(1) from each of the slab
surfaces, for the full (half) coverage. As for EBatra, the obtained binding
energy for the half coverage HH case is -2.21eV, being close to Batra's valueof -2.23eV [71]. The other quantities are shown in Table 3.6.
In this table, the binding energy for the full coverage (HH-T3) is larger
than that of the half coverage (HH or T4). And the binding energy forthe  = 61 is the largest. It should be noted that the condition (Ecut andNk ) for the full coverage is di erent from that for the half coverage andthat the binding energies may depend on this di erence. Nevertheless, an
important conclusion here is that the binding energy is still signi cant for
the full coverage.
In order to understand the real situation for the saturation coverage,
more detailed studies including temperature e ects are necessary. At thepresent stage, combining several experimental results as well as calculations, we speculate that the coexistence of the HH-site adsorption and theT4-(or T3-) site adsorption may be realized. The adatoms adsorbed along
the trough are expected to give only a small in uence on the workfunction,
LEED, and STM.
3.2.4

1
6

coverage Na/Si(001)-221 surface

Figure 3.11(a) is the STM image for low coverage K on the Si(001) surface[28, 29, 30], and the corresponding schematic diagram for the atomic arrangement is shown in Fig. 3.11(b). In Fig. 3.11(a), the bright large spots
correspond to adsorbed K and the neighboring zig-zag spots correspondto raised Si atoms along the alternating arrangement of the asymmetricdimers. This result suggests that the adsorption site for low coverage K
is none of those sites considered for  = 1 and 0:5. It is the B2 site in
our notation, i.e., the bridge site between the two raised Si atoms of theneighboring Si-Si dimers. It is also interesting to note that the adsorbed
K atoms seem to enhance the contrast of the STM image associated with35

the asymmetry of the neighboring Si-Si dimers [28, 29, 30]. We attempt
to analyze these aspects of the STM observations by using a supercell for = 61 .
Because of the large system size, the present calculation uses a relatively
small Ecut of 5.0 Ry corresponding to the number of plane waves of 3072.
Nk = 4 is also small. The present calculation is therefore a preliminary one.
The starting atomic con guration for the substrate is the 221-asymmetricdimer structure which was optimized for the clean surface. The Na atom is
initially put slightly o the HB or T3 site. The crosses in Fig.3.12(a) denote
two initial Na positions. As the optimization process proceeds startingfrom either of the two initial positions, the Na atom moves toward the interdimer bridge B2 site as denoted by a hatched sphere in Fig.3.12(a). This is
the position that the STM image [28, 29, 30] suggests as the adsorption sitefor K in the low coverage region. Figure 3.12(b) is a side view containing
Na and the Si atoms designated as A, A', B and B' in Fig. 3.12(a). In theabove calculation, a ferromagnetic ordering of the asymmetric dimers forthe substrate was assumed as a initial con guration and was not destroyedin the optimized con guration. However, as mentioned in Section 2.4.2, thekinetic energy of the ionic system dissipates quickly in the present MSD
approach so that the system may be easily trapped at a local minimum.
Therefore two additional initial arrangements of the asymmetric dimers are
considered. First, the reversal of the asymmetry for the dimer designatedas (C,C') in Fig. 3.12(a) lowers the total energy by 19.1meV per (223)
surface unit cell without signi cant change in the geometry in the vicinityof the adsorbed Na. This result is consistent with the zig-zag image along
the dimer row near the adsorbed K. However, for the initial con guration
shown in Fig. 3.13, the optimum Na position turns out to be the hollow
site as shown by a hatched sphere. However, the total energy of thiscon guration is higher than that of the previous one by 136meV per (223)
surface unit cell.
The present coverage  = 61 may not be low enough to make a directcomparison with the experiments and there may also be some di erence
between K and Na. Still it is very interesting that the present calculation
does indicate a dramatic change in the preferred adsorption site between
the high coverage and low coverage cases.
Here, we show the calculated electron density corresponding to the36con guration of Fig. 3.12 in order to compare with the results of STM[28, 29, 30]. The contour maps of the electron density in the side viewand top view are shown in Figs. 3.14(a)-(h). The cut plane of side viewis shown in the attached small gure. The Na atom is indicated by anarrow. In Fig.3.14(a), the electron density on the dimer is asymmetric andthat arising from the adsorbed Na is not identi ed as in the full coveragecase. In the side view of the electron density shown in Fig. 3.14(b), thearrows along the vertical axis denote the height of the cut plane for the
top views (Figs. 3.14(c)-(h)). In these gures for the top view which may
be compared with the STM results, the position of the adsorbed Na isaround the center of the gure as shown by the hatched circle. The heightof the cut plane in Fig. 3.14(c) is that of the raised Si atom of the asymmetric dimer. The height in Fig. 3.14(d) is that of the adsorbed Na atom.
In Figs. 3.14(e)-(h), the electron density of the adsorbed Na can be seenclearly. The position of the cut plane in Fig. 3.14(e) is higher than that of
the Na atom by 2.1 A. The bird's eye view corresponding to Fig. 3.14(f) is
shown in Fig. 3.14(g). The electron pro les given in Figs. 3.14(e)-(h) areconsistent with the observed STM image [28, 29, 30].
3.2.5 Full coverage Li/Si(001)-221 surface

As Li is chemically fairly di erent from other alkali metal elements,
it is interesting to see whether any qualitatively di erent aspect can befound between Li and Na as adsorbates on the Si surface. Even within thelimited size of the present unit cell, we have in fact found that Li behavesquite di erently from Na. The HB site is unstable for Li, too. However,unlike the case of Na, the Li atom along the trough is trapped at eitherof the two o -symmetry positions as shown in Fig. 3.15. Accordingly, theadsorption site of Li along the dimer chain is also slightly o the HH site.
The two con gurations are energetically almost degenerate with the energy
di erence of only 5.0meV. Nevertheless, the relaxation of the substrate ismore or less similar to that for the Na case. The positional parameters
corresponding to Fig. 3.15 are listed in Tables 3.7 and 3.8. Although themost stable atomic arrangement for Li at the saturation coverage should
be searched for with probably a much larger unit cell to allow more degrees
of freedom, the present calculation demonstrates very clearly a qualitative37

di erence between Li and Na as adsorbates on the Si(001) surface. It was
experimentally shown that Li behaves completely di erently from otheralkali atoms as adsorbates on the Si(001) surface in the high coverageregime (Silicide formation) [74].
coverage Li/Si(001)-221 surface

3.2.6

1
6

3.3

Partial core correction

The STM image for low coverage Li on the Si(001) surface [28, 29, 30]
and its schematic illustration are shown in Fig. 3.16. Unlike the case for K
shown in Fig. 3.11, the adsorption site for low coverage Li is the B1 site in
Fig. 3.1, i.e., the on-top site of the raised Si atom of an asymmetric Si-Si
dimer. The purpose of the study for this coverage of Li is the same as for the
Na  = 16 case. The condition of the calculation is also completely the sameas the Na case. The initial Li atoms are put slightly o the HB or T3 site
as in the case of Na in Fig. 3.12(a). As the optimization process proceedsfrom either of the two initial positions, the Li atom moves toward the interdimer bridge B2 site as in the Na case. The heights of adsorbed Li is lower
than that of adsorbed Na. In addition, the other substrate con guration
as the Na case of Fig. 3.13 is used and the optimized position is almost
the same (HH site) as in the Na case. (The height of Li is also lower thanthat of Na.) These results do not agree with the STM measurement for Li.
This inconsistency may arise from the fact that our surface unit cell is notlarge enough to simulate the experimental situation. Another problem isthat as Li has a relatively deep pseudopotential, the present cuto energy(5Ry) is de nitely insucient. More extended studies will be performed innear future.
In all of the above results, we noticed an inconvenient aspect in thetheoretically optimized structures : the interatomic distances in the surface region are signi cantly shorter than those estimated by experiments.
These shorter distances are caused by the following reason. The pseudopotential for the core part is obtained from the neutral pseudoatomic
potential by subtracting the Hartree and exchange-correlation potentials
for valence electrons as described in Section 2.1. The subtracted Hartree38

and exchange-correlation potentials are made of only valence electrons. Although this unscreening process is exact for the Hartree potential becauseof its linear dependence on the electron density, it is only an approximation for the exchange-correlation potential. This simpli ed process for theexchange-correlation potential is justi ed if the core electron density is spatially separated from the valence electron density. Such a situation is notnecessarily realized in alkali atoms and the approximate treatment for theexchange-correlation potential has a tendency of underestimating the bondlength. It is necessary to introduce a partial core correction [31] to correct
this underestimation.
We introduce this partial core correction for the adsorbed Na and thesubstrate Si atoms and study the di erence from previous results for thefull and half coverage cases. The optimized con gurations without thepartial core correction are adopted as the initial con gurations for thiscorrection. The obtained positional parameters and interatomic distances
are listed in Table 3.2.
For the full coverage cases (HH-T3 and HH-T4) as well as for the halfcoverage ones (HH and T4), the interatomic distances in the surface region
become signi cantly larger. The experimental values of the Na-Si distancefor the saturation coverage are 2.975 A [24] from LEED analysis and 2.8
A [77] by surface EXAFS. If we assume that the experimental saturation
coverage corresponds to  = 1.0, our new value of 3.13 A for the HH-T3 caseis slightly larger than the experimental values. A qualitatively new aspectis that the adsorbed atom on the T4 site moves to the T3 site and stops
there for the full coverage. Therefore, in the case of considering the partial
core correction, the relative stability of the optimized con guration of HHT3 with reference to that of HH-T4 is enhanced. For the half coverage, theoptimized con guration of T4 becomes signi cantly more stable than thatof HH by 100meV per dimer. The corresponding value without the partial
core correction is only 12meV per dimer. Nevertheless, we nd that thequalitative aspect of the most stable con guration is not a ected by thepartial core correction.

39

Chapter 4

Concluding remarks
We have performed the FPMD calculation for the clean and alkali adsorbed Si(001) surfaces. The results are summarized as follows.
(1) Clean Si(001)-221 surface
 The asymmetric dimer is more stable than the symmetric one.
(2) Na adsorbed Si(001) surfaces
1. Full coverage The con guration of HH-T3 is more favorable than that of HH-T4(consistent with the XPD(X-ray photoelectron di raction) result[26]).
 The HB site is unstable.
 The asymmetric dimer of substrate changes into the symmetric one.
2. Half coverage The con guration of T4 is more favorable than that of HH (inconsistent with the LEED result [24]).
 The HB site is unstable.
 The asymmetry in the Si-Si dimer becomes weaker than that of theclean substrate.
3. 61 coverage The optimized site of Na is on the bridge site (B2) between the tworaised Si atoms of the neighboring Si-Si dimers (consistent with the
STM result [28, 29, 30]).
(3) Li adsorbed Si(001) surfaces40

1. Full coverage The Li atoms are adsorbed on o -symmetry sites as shown in
Fig. 3.15.
2. 61 coverage The optimized site of Li is same as that of Na (inconsistent with
STM result [28, 29, 30]).
(4) Partial core correction
 The interatomic distances of Na-Si become much closer to experimental
values by considering the partial core correction [31].
We conclude from above results as follows. The stable con guration forthe half coverage is such that the Na atoms are adsorbed on the T4 sites
with forming linear chains of Na adatoms between dimer rows. Ling etal. [75] also arrived at the same conclusion for the K case, although they
did not take account of the asymmetry of the Si-Si dimer. Even for thefull coverage, the Na atoms can be chemisorbed on the Si substrate with asigni cant binding energy. The result seems to suggest that the saturation
coverage may be  = 1. For the low coverage ( = 61 ) case, the optimumsite of the adsorbed Na is the B2 which is di erent from the adsorption
sites for the higher coverage ( = 0:5 or 1:0) cases. The present FPMD
calculation has been successful in reproducing the di erence of the situationof adsorption for the low and high coverage cases but not in reproducing
the STM result for the low coverage of Li adsorption.
There remain controversies on the saturation coverage. Enta et al. [72]suggested that the saturation coverage for Na is  = 1:0. On the other
hand, Wei et al. [24] suggested that  = 0:5 is the saturation coverage.
Furthermore, there is also a possibility  = 0:68 as suggested by Glander
et al. [22, 23].
As mentioned above, our calculation supports Enta et al.'s [72] suggestion. However, we should note that the present calculations have covered
only very limited possibilities for the atomic arrangement. Further, we
assumed a priori the (2 2 1) unit cell for  = 0:5 and  = 1. Therefore,
we cannot rule out a possibility of incommensurate con guration in the
dimer direction: for example, four Na atoms on average per ve Si dimers.
41

Theoretical analysis of such a possibility requires very heavy calculations
and is left as a future task. Nevertheless, the present calculations strongly
suggest the necessity of modifying the Levine model. All calculations performed up to present (the present calculations both with and without thepartial core correction and the calculation by Ling et al. [75]) suggested
adsorption of Na and K along the trough. In most experiments, Na seems
to be adsorbed both at HH and T3 (or T4). This di erence is highly likely
to be due to nite temperature e ects.
The cuto energy, the number of k-points for sampling, the size of theunit cell, the number of plane waves, and the CPU-time of supercomputer(Hitac S820/80) for each calculation are listed in Table 4.1 in order toshow the performance of the present calculations. It is found that thecalculations for the half or 61 coverage case are already close to the practical
upper limit of the computer resource.
Some of the problems existing in the present calculations are listed below.
 The speed of the convergence for the SD calculation is slow for anelectronic part. Besides, it is rather dicult to obtain the global minimum of the structure for ions. A future test of some other algorithms
is necessary. (For example, the original Car-Parrinello method, theconjugate gradient method, etc [58, 59, 78])
 The temperature control by a constant temperature MD for an ionicpart should be introduced.
 A shape of the unit cell has been xed in the process of our FPMD
calculation. The constant pressure MD for an ionic part [79, 80] isrequired for a more strict test of the structural stability.
 The materials which we can study by pseudopotentials with plane wave
basis is fairly limited. For example, noble and transition metals, B,
C, O, N [with deep 2p potentials], and so on cannot be treated so
accurately. It is necessary to devise a new type of pseudopotentials[81, 82, 83, 84].
Finally, it should be noted that the FPMD method is useful for the
search and design of new materials and the present study is a rst step
toward its real application.
42

Appendix A
Separable form for non-local pseudopotentials

Equation (2.9) isVps(r) = Vlocal (r) +

(2:9)
l
l
In the second term of the right hand side of this eqaution, Vnon
0local (r) is
replaced by the Kleinman-Bylander form [85]
l
l
j
Vnon
j
l >< l j Vnon
l;KB
0
local
0local j ;Vnon0local (r) =(A.1)
l< l j Vnon
0local j l >l;KBwhere Vnon
0local (r) is a separable non-local pseudopotential and l (r ) is a
pseudo-wave function for an atomic eigen-state. Here, the pseudo-eigen
energies for the separable and non-separable forms are obtained as follows
Xlj l > Vnon
0local (r ) < l j :l;KBKB< l j Vnon
0local j l >= l ;

(A.2)

l
< l j Vnon
0local j l >= l ;where the KB
is equal to l .
l

(A.3)

But, the equality of these eigen-energies is not guaranteed in the generalcase with a wave function 9l;KBl< 9 j Vnon
0local j 9 >6=< 9 j Vnon0local j 9 > :(A.4)
The di erence in these eigen-energies can be reduced within a sucientaccuracy by choosing the local pseudopotential of eq. (2.9). At present,
there is no problem as for the accuracy of our calculation where the dpseudopotential is adopted for the local pseudopotential.
43

Appendix B
Lagrange multiplier in steepest descent method

The following three equations are our starting point:X 9_ i (t) = 0H (t)9i(t) + 3ij 9j (t);j9i(t + 1t) = 9i(t) + 1t9_ i(t);

(B.1)
(B.2)

Z

93i (t)9j (t)d~r = ij :(B.3)
The Lagrange multiplier 3ij is determined so as to satisfy the orthonormalization constraint at t + 1t
Z

93i (t + 1t)9j (t + 1t)d~r = ij :Inserting eq. (B.2) into the left hand side of eq. (B.4), we obtain
Z

(B.4)

f9i(t + 1t)g3f9j (t + 1t)gd~r
Z= Z f9i(t) + 1t9_ i(t)g3fZ9j (t) + 1t9_ j (t)gd~r= 93i (t)9Z j (t)d~r + 1t f9_ 3i (t)9j (t) + 93i (t)9_ j (t)gd~r+(1t)2 9_ 3i (t)9_ j (t)d~r:(B.5)
By neglecting the second order term of 1t, eq. (B.4) is reduced toZ
f9_ 3i (t)9j (t) + 93i (t)9_ j (t)gd~r = 0:
(B.6)
Substituting eq. (B.1) for 9_ i(t) into eq. (B.6), we obtain
44

Z

Z

(0H 9i(t) + X 3il 9l(t))39j (t)d~r + 93i (t)(0H 9j (t) + X
3jm9m(t))d~rml= 2(0 < 9i j H j 9j > +3ji) = 0:
(B.7)
Therefore, 3ji is given by
3ji =< 9i j H j 9j >:

45

(B.8)

Appendix C
Approximation to diagonal Lagrange multipliers

Here, the approximation in the calculation of i;~k in Section 2.4.1 isdescribed in detail. Equation (2.62) is written here again
i;~k (t + 1t) =~ t + 1t) < G~ j H (t) j G~ 0 > 9i;~k (G~ 0; t + 1t)
9i;~k (G;~ G~
G;~ t)9i;~k (G;~ t + 1t) 0 93 ~ (G;~ t + 1t)9_ i;~k (G;~ t)= i;~k 0  Xf9_ 3i;~k (G;i;kG~~ t)9_ i;~k (G;~ t)g:+93i;~k (G;(2:62)In order to derive this equation, we proceed in the following way:
X

~ G~ 0G;X

0~ t + 1t) < G~ j H (t) j G~ 0 > 9i;~k (G~ 0 ; t + 1t)
9i;~k (G;~ t) +  93 ~ (G;~ t)) < G~ j H (t) j G~ 0 > (9i;~k (G~ 0 ; t) + 9i;~k (G~ 0 ; t))= X (93i;~k (G;i;k~ G~ 0G;~ t) < G~ j H (t) 0 i;~k (t) j G~ 0 > 9i;~k (G~ 0; t)= X 93i;~k (G;~ G~ 0G;~ t) < G~ j H (t) 0 i;~k (t) j G~ 0 > 9i;~k (G~ 0 ; t)
+ X 93i;~k (G;~ G~ 0G;~ t) < G~ j H (t) 0 i;~k (t) j G~ 0 > 9i;~k (G~ 0 ; t)
+ X 93i;~k (G;+~ G~ 0G;

X

vecG;G~ 0~ t) < G~ j H (t) 0 i;~k (t) j G~ 0 > 9i;~k (G~ 0 ; t)93i;~k (G;~ t) +  93 ~ (G;~ t)) < G~ j i;~k (t) j G~ 0 > (9i;~k (G~ 0 ; t) + 9i;~k (G~ 0 ; t))+ X (93i;~k (G;
i;k
~ G~
G;
= i;~k (t)046

~ t) < G~ j H (t) 0 i;~k (t) j G~ 0 > 9i;~k (G~ 0; t)+ X 93i;~k (G;~ G~ 0G;+ X 93i;~k < G~ j H (t) 0 i;~k (t) j G~ 0 > 9i;~k (G~ 0; t)~ G~ 0G;~ t) < G~ j H (t) 0 i;~k (t) j G~ 0 > 9i;~k (G~ 0 ; t) (1t2 )+ X 93i;~k (G;~ G~G;= i;~k (t) 0  X 9_ 3i;~k (G~ 0; t)9i;~k (G~ 0; t)G~X
0  93i;~k (G~ 0; t)9_ i;~k (G~ 0; t) 0  X 9_ 3i;~k (G~ 0; t) 9i;~k (G~ 0; t) + (1t2)G~G~= i;~k (t)~ t)[9i;~k (G;~ t) +  9i;~k (G;~ t)]
0  X 9_ 3i;~k (G;~G
X 3~ t) + 93i;~k (G;~ t)]9_ i;~k (G;~ t)
0  [9i;~k (G;~G
X 3~ t)9_ i;~k (G;~ t):+  9i;~k (G;
0

0

0

0

~
G

(C.1)

The nal form is obtained as~ t)9i;~k (G;~ t + 1t)
f9_ 3i;~k (G;~G~ t + 1t)9_ i;~k (G;~ t) + 93 ~ (G;~ t)9_ i;~k (G;~ t)g; (2:62)
0 93i;~k (G;i;kwhere 9 means the change in the wave function after one time step including the Gram-Schmidt orthogonalization process. Schematically,i;~k (t + 1t) = i;~k (t) 0 X~ t + 1t) = 9i;~k (G;~ t) + 1t9_ i;~k (G;~ t);9i;~k (G;+Gram 0 Schmidt orthogonalization+~ t + 1t) = 9i;~k (G;~ t) +  9i;~k (G;~ t):9i;~k (G;We also made use of the relation(2:59)~ t + 1t)9i;~k (G;~ t + 1t) = 1:93i;~k (G;

(C.3)

X

~
G

47

(C.2)

Acknowledgements

I would like to express my sincere thanks to Professor K. Terakura forcontinual guidance and encouragement throughout the course of the workof this thesis. I wish to thank Dr. S. Blugel for many helpful discussions andsuggestions. I also wish to thank Mr. Y. Morikawa for helpful assistancein the calculations. The calculation by considering the partial core correction and the preparation of pseudopotentials with this correction werecarried out with his assistance. I am grateful to Dr. H. Ishida for givingme the starting version of the program for the electronic structure calculation. Several helpful discussions with Drs. T. Oguchi, T. Sasaki, Y. Ishii,
S. Takizawa, M. Ikeda and T. Yamazaki are gratefully acknowledged.
The numerical calculations were performed at the Computer Centers of
University of Tokyo, Institute for Molecular Science and Institute for Solid
State Physics.

48

Bibliography[1] Hans L. Skriver. Phys. Rev. B, 31(4):1909, 1985.
[2] N. E. Cristensen. Phys. Rev. B, 32(32):207, 1985.
[3] H. Hellmann. Einfuhrung in die Quanten Theorie. (Deuticke, Leipzig),
1937.
[4] R. P. Feynman. Phys. Rev., 56:340, 1939.
[5] J. C. Slater. J. Chem. Phys., 57:2389, 1972.
[6] R. Car and M. Parrinello. Phys. Rev. Lett., 55(22):2471, 1985.
[7] G. Binnig, H. Rohrer, Ch. Gerber and E. Weibel. Phys. Rev. Lett.,
50:120, 1983.
[8] R. M. Tromp, R. J. Hamers and J. E. Demuth. Phys. Rev. Lett.,
55:1303, 1985.
[9] R. J. Hamers, Ph. Avouris and F. Bozso. Phys. Rev. Lett., 59(18):2071,
1987.
[10] T. Tabata, T. Aruga and Y. Murata. Surf. Sci., 179:L63, 1987.
[11] Y. Enta, S. Suzuki and S. Kono. Phys. Rev. Lett., 65(21):2704, 1990.
[12] J. Ihm, Marvin L. Cohen and D. J. Chadi. Phys. Rev. B, 21(10):4592,
1980.
[13] M. T. Yin and Marvin L. Cohen. Phys. Rev. B, 24(4):2303, 1981.
[14] N. Roberts and R. J. Needs. J. Phys.: Condens. Matter, 1:3139, 1989.
[15] Zizhong Zhu, Nobuyuki Shima and Masaru Tsukada. Phys. Rev. B,
40(17):11868, 1990.
[16] D. J. Chadi. Phys. Rev. B, 43:43, 1979.
49

[17] N. Roberts and R. J. Needs. Surf. Sci., 236:112, 1990.
[18] Tsuyoshi Uda and Shigeo Ihara. Solid State Physics, 25(11):867, 1990.
[19] S. Ihara, S. L. Ho, T. Uda and M. Hirao. Phys. Rev. Lett., 65(15):1909,
1990.
[20] D. Levine. Surf. Sci., 34:90, 1973.
[21] T. Aruga, H. Tochihara and Y. Murata. Phys. Rev. Lett., 53:372, 1984.
[22] G. S. Glander and M. B. Webb. Surf. Sci., 222:64, 1989.
[23] G. S. Glander and M. B. Webb. Surf. Sci., 224:60, 1989.
[24] C. M. Wei, H. Huang, S. Y. Tong, G. S. Glander and M. B. Webb.
Phys. Rev. B, 42:11284, 1990.
[25] Y. Enta, T. Kinoshita, S. Suzuki and S. Kono. Phys. Rev. B, 36:9801,
1987.
[26] T. Abukawa and S. Kono. Phys. Rev. B, 37:9097, 1987.
[27] S. Tanaka, N. Takagi, N. Minami and M. Nishijima. Phys. Rev. B,
42:1868, 1990.
[28] Yukio Hasegawa, I. Kamiya, T. Hashizume, T. Sakurai, H. Tochihara,
M. Kubota and Y. Murata. Phys. Rev. B, 41(14):9688, 1990.
[29] T. Hashizume, I. Sumita, Y. Murata, S. Hyodo and T. Sakurai. to be
published in J. Vac. Sci. Technol.

[30] Toshio Sakurai and T. Hashizume. to be published.
[31] Steven G. Louie, S. Froyen and Marvin L. Cohen. Phys. Rev. B,
26:1738, 1982.
[32] B. J. Austin, V. Heine and L. J. Sham. Phys. Rev., 127(1):276, 1962.
[33] Volker Heine. Solid State Physics. Volume 24, New York, Academic
Press, 1970.
[34] Marvin L. Cohen and Volker Heine. Solid State Physics. Volume 24,
New York, Academic Press, 1970.
[35] Volker Heine and D. Weaire. Solid State Physics. Volume 24, New
York, Academic Press, 1970.
50

[36] Robert W. Shaw, Jr. and Walter A. Harrison. Phys. Rev., 163(3):604,
1967.
[37] D. R. Hamann, M. Schluter and C. Chiang. Phys. Rev. Lett.,
43(20):1494, 1979.
[38] G. B. Bachelet, D. R. Hamann and M. Schluter. Phys. Rev. B,
26(8):4199, 1982.
[39] J. Ihm and Marvin L. Cohen. Solid State Commun., 29(10):711, 1979.
[40] Leonard Kleinman and D. M. Bylander. Phys. Rev. Lett., 48(20):1425,
1982.
[41] D. Hohl, R. O. Jones, R. Car and M. Parrinello. Chem. Phys. Lett.,
139:540, 1987.
[42] D. Hohl, R. O. Jones, R. Car and M. Parrinello. J. Chem. Phys. ,89:6823, 1988.
[43] D. Hohl, R. O. Jones, R. Car and M. Parrinello. J. Am. Chem. Soc.,
111:826, 1989.
[44] O. Gunnerson, B. I. Lundqvist and S. Lundqvist. Solid State Commun., 11:149, 1972.
[45] O. Gunnerson and B. I. Lundqvist. Phys. Rev. B, 13:4274, 1976.
[46] John P. Perdew and Alex Zunger. Phys. Rev. B, 23:5028, 1981.
[47] P. Hohenberg and W. Kohn. Phys. Rev., 136:B864, 1964.
[48] W. Kohn and L. J. Sham. Phys. Rev., 140:A1133, 1965.
[49] L. V. Moruzzi, J. F. Janak and A. R. Williams. Pergamon, New York,
1978.
[50] D. M. Ceperley. Phys. Rev. B, 18:3126, 1978.
[51] E. P. Wigner. Phys. Rev., 46:1002, 1934.
[52] M. T. Yin and Marvin L. Cohen. Phys. Rev. B, 26(6):3259, 1982.
[53] J. Ihm, Alex Zunger and Marvin L. Cohen. J. Phys. C, 12(21):4409,
1979.
[54] G. Lehmann and M. Taut. phys. stat. sol., 54:469, 1972.
51

[55] J. Rath and A. J. Freeman. Phys. Rev. B, 11(6):2109, 1975.
[56] L. Verlet. Phys. Rev., 159:98, 1967.
[57] Jean-Paul Ryckaert, Giovanni Ciccotti and Herman J. C. Berendsen.
J. Comput. Phys., 23:327, 1977.
[58] I. Stich. PhD thesis, International School for Advanced Studies, 1987.
[59] I. Stich, R. Car, M. Parrinello and S. Baroni. Phys. Rev. B, 39(8):4997,
1989.
[60] M. C. Payne, J. D. Joannopoulos, D. C. Allan, M. P. Teter and David
H. Vanderbilt. Phys. Rev. Lett., 56(24):2656, 1986.
[61] M. C. Payne, M. Needels and J. D. Joannopoulos. Phys. Rev. B,
37(14):8138, 1988.
[62] A. R. Williams and J. Soler. Bull. Am. Phys. Soc., 32:562, 1987.
[63] M. C. Payne, D. C. Allan, M. P. Teter and J. D. Joannopoulos. to be
published.
[64] C. Woodward, B. I. Min, R. Benedek and J. Garner. Phys. Rev. B,
39(8):4853, 1989.
[65] Atsushi Oshiyama and Mineo Saito. J. Phys. Soc. Jpn., 56(6):2104,
1987.
[66] E. Holzschuh. Phys. Rev. B, 28(12):7346, 1983.
[67] Michel M. Dacorogna and Marvin L. Cohen. Phys. Rev. B, 34(8):4996,
1986.
[68] Charles Kittel. Introduction to Solid State Physis. John Wiley andSons, Inc., 1986.
[69] R. J. Hamers, R. M. Tromp and J. E. Demuth. Phys. Rev. B, 34:5343,
1986.
[70] Inder P. Batra. Phys. Rev. B, 39(6):3919, 1989.
[71] Inder P. Batra. to be published.
[72] Y. Enta, S. Suzuki and S. Kono. to be published in Surface Science,
1991.
52

[73] T. Aruga and Y. Murata. Phys. Rev. Lett., 53:372, 1984.
[74] Hiroshi Tochihara and Yoshitada Murata. Surf. Sci., 215:L323, 1989.
[75] Ye Ling, A. J. Freeman and B. Delley. Phys. Rev. B, 39(14):10144,
1989.
[76] H. Ishida and K. Terakura. Phys. Rev. B, 40(17), 1989.
[77] Patrick Soukiassian. to be published.
[78] Tamio Oguchi and Taizo Sasaki. Solid State Physics, 25(11):857, 1990.
[79] M. Parrinello and A. Rahman. Phys. Rev. Lett., 45:1196, 1980.
[80] M. Parrinello and A. Rahman. J. Appl. Phys., 52:7182, 1981.
[81] D. R. Hamann. Phys. Rev. B, 40(5):2980, 1989.
[82] David Vanderbilt. Phys. Rev. B, 32(12):8412, 1985.
[83] David Vanderbilt. Phys. Rev. B, 41(11):7892, 1990.
[84] E. L. Shirley, D. C. Allan, R. M. Martin and J. D. Joannopoulos.
Phys. Rev. B, 40:3652, 1989.
[85] Leonard Kleinman. Phys. Rev. B, 21(6):2630, 1980.

53

Figure Caption
 Fig.1.1. (a) A side view of the ideal clean Si(001) surface, (b) the clean
Si(001) surface with symmetric dimer, (c) the clean Si(001) surface
with p(222) con guration. These gures originate from
Introduction to Surface P hysics (in Japanese) by M. Tsukada (University of Tokyo Press, 1989).
 Fig. 2.1. A ow chart of the standard electronic structure calculation.
 Fig.2.2. A process of the FPMD calculation. The vertical line measuresa total energy. The horizontal line means a con guration of ions.
 Fig. 2.3. The ow chart of the SD or MSD calculation.
 Fig.3.1. Top (a) and side (b) views of the asymmetric dimer con guration of the Si(001)-221 surface. The lled and the open circles denote
the Si atoms in the topmost surface layer and those underneath, respectively. The possible adsorption sites for alkali atoms are denoted
as HB, HH, T3, T4, B1 and B2.
 Fig.3.2. A supercell including atomic layers and vacuum region for the
repeated slab model.
 Fig. 3.3. A bird's eye view of the electron density for the Na adsorbed
Si(001) surface for the full coverage. The cut plane is parallel to they-z plane and contains the Si-Si dimer. The axes x, y and z are shownin Fig. 3.1.
 Fig. 3.4. Positions of k-points for sampling in the whole surface Brillouin zone for Si(001)-221 unit cell. The open circles denote the kpoints. The numbers of k-points are 8 in (a) and 32 in (b).

54 Fig. 3.5. The energy dispersion curves of  and 3 bands for the symmetric and asymmetric dimer con gurations. Solid lines are for thesurface bands of the symmetric dimer and broken lines for those of
the asymmetric dimer. The regions of the bulk bands are indicated byhatching. The position of the Fermi energy for  = 0.0, 0.5 and 1.0 are
schematically shown.
 Fig.3.6. Contour map of the electron density for the clean surface withthe asymmetric dimer. The cut plane is parallel to the y-z plane andcontains the Si-Si dimer.
 Fig. 3.7. (a),(b) Top views of the optimized atomic positions for the
Na adsorption with the full coverage in the HH-T3 and HH-T4 congurations, with the hatched circles denoting Na. (c) A side viewcorresponding to the HH-T3 con guration. The dotted circles withand without cross denote the asymmetric dimer of the clean Si surface and the ideal bulk terminated Si surface, respectively. The circles
denoted by A and B are Na.
 Fig. 3.8. The total energy (a) and the temperature (b) curves for thecon guration of HH-T3. The vertical lines of (a) and (b) measure thetotal energy (Hartree) and the temperature (Kelvin), respectively. The
horizontal lines means a number of iteration steps in the process of thecalculation.
 Fig. 3.9(a). Contour map of the electron density for the Na adsorbed
Si(001) surface for the full coverage (HH-T3). The lled circles denotethe positions of the Na and Si atoms projected on the cut plane. The
arrows indicate Na atoms. This is a side view for the cut plane shownby the attached small gure.
 Fig. 3.9(b). Same as in Fig. 3.9(a).
 Fig.3.9(c). Same as in Fig.3.9(a). This is a di erence electron density.
 Fig.3.9(d). Same as in Fig.3.9(a). This is a top view for the di erenceelectron density. Z(height) is a height from the bottom of the unit cellas shown in Fig. 3.2.
 Fig. 3.9(e). Same as in Fig. 3.9(a). This is a top view. Z(height) is aheight from the bottom of the unit cell as shown in Fig. 3.2.
55

 Fig. 3.9(f). The bird's eye view corresponding to Fig. 3.9(e).
 Fig. 3.10. Top and Side views of the optimized atomic con guration
for the Na adsorption with  = 0.5. The hatched circles denote the Na
atoms and the others are the substrate Si atoms. The dotted circles
with and without cross denote the ideal bulk terminated Si surface
and the asymmetric dimer of the clean surface, respectively. (a) and(c) are for Na at the HH site and (b) and (d) for Na at the T4 site.
In the present calculation, case (b) is more stable than case (a) by
12meV/dimer.
 Fig. 3.11. (a) STM image for low coverage K on the Si(001) surface.
The bright large spots correspond to adsorbed K and the neighboring
zig-zag spots correspond to raised Si atoms. (b) The corresponding
schematic diagram for the atomic arrangement. The large double circle
denotes the adsorbed K atom and the lled circles denote the raised
Si atoms. (From reference [30]) Fig.3.12. (a) A top view of the 223 surface unit cell for the 61 coverageof Na. The crosses denote two initial positions for Na and the hatchedcircle denotes Na at the optimized adsorption site. The lled andopen circles denote the substrate Si atoms in the topmost surface layer
and the second layer, respectively. (b) The corresponding side view.
The hatched circle denotes the Na atom. The thin dotted circles withand without cross denote the ideal bulk terminated Si surface andthe asymmetric dimer of the clean surface, respectively. Appreciable
amount of displacements can also be seen in the second layer. The Si
atoms designated by A, A', B, B', C and C' correspond to those in (a).
 Fig. 3.13. A top view of the 223 surface unit cell for the 61 coverage of
Na. The initial con guration of the Si substrate is di erent from thatof Fig. 3.12(a). The cross denotes the initial position of Na and thehatched circle denotes Na at the optimized adsorption site. The lled
and open circles denote the substrate Si atoms in the topmost surface
layer and the second layer, respectively.

56 Fig. 3.14(a). Contour map of the electron density for the Na adsorbed
Si(001) surface for  = 16 with Na on the B2 site. The lled circles
denote the positions of the Na and Si atoms projected on the cut
plane. The arrow indicates the Na atom. This is a side view for thecut plane shown by the attached small gure.
 Fig. 3.14(b). Same as in Fig. 3.14(a). Each arrow along the verticalaxis indicates the height of the cut plane for the top view in each of
Figs. 3.14(c)-(h).
 Fig. 3.14(c). Same as in Fig. 3.14(a). This is a top view. Z(height) isa height from the bottom of the unit cell as shown in Fig. 3.2. The
hatched circle denotes the adsorbed Na atom.
 Fig. 3.14(d). Same as in Fig. 3.14(c).
 Fig. 3.14(e). Same as in Fig. 3.14(c).
 Fig. 3.14(f). Same as in Fig. 3.14(c).
 Fig. 3.14(g). The bird's eye view corresponding to Fig. 3.14(f).
 Fig. 3.14(h). Same as in Fig. 3.14(c).
 Fig.3.15. Top views of two possible optimized atomic positions (a) and(b) for the Li adsorption with the full coverage. The large lled andopen circles denote the Si atoms in the top-most and second substratelayers. The small hatched circles denote the adsorbed Li atoms. The
con gurations (a) and (b) are nearly degenerate with (a) having a lower
energy by only about 5 meV.
 Fig. 3.16. (a) STM image for low coverage Li on the Si(001) surface.
The bright large spots correspond to adsorbed Li and the neighboring
zig-zag spots correspond to raised Si atoms. (b) The corresponding
schematic diagram for the atomic arrangement. The double circle near
the center of the gure denotes the adsorbed Li atom and the lled
circles denote the raised Si atoms. (From reference [30])57