AN LCAO LOCAL DENSITY FUNCTIONAL APPROACH TO

SURFACE ELECTRONIC STRUCTURE CALCULATIONS

BY

JOHN WALLACE MINTMIRE

A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL

OF THE UNIVERSITY OF FLORIDA

IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE

DEGREE OF DOCTOR OF PHILOSOPHY

UNIVERSITY OF FLORIDA

1980

This dissertation is dedicated to those members of my family who

wondered when I would ever graduate.

ACKNOWLEDGMENTS

I express my deepest appreciation to Prof. J. R. Sabin for the

support and guidance which has made this dissertation possible. I would

also like to express my gratitude to all the members of the Quantum

Theory Project for the camaraderie which has made my period of graduate

study a pleasant and productive one, with a special recognition of Prof.

P.-O. Lbwdin for the assistance he has provided in the development of my

professional background. I am deeply indebted to Prof. S. B. Trickey,

Dr. B. 1. Dunlap, and Dr. J. P. Worth for their invaluable advice and

help in the course of my research. The editorial assistance of us.

Cynthia Kinney is most sincerely appreciated.

The financial support of the Northeast Regional Data Center and the

National Science Foundation is gratefully acknowledged. I also acknowledge

the use of computational facilities provided by the Department of Highway

Safety and Motor Vehicles of the State of Florida.

TABLE OF CONTENTS

PAGE

ACKNOWLEDGMENTS .................................................... iii

ABSTRACT ........................................................... vi

CHAPTER

I. INTRODUCTION .............................................. I

1-1. General Remarks................................... 1

1-2. Historical Background............................. 3

1-3. Organization of Dissertation ...................... 8

II. MATHEMATICAL FORMALISM AND THEORY......................... 9

2-1. Symmetry Properties of Thin Films................. 9

2-2. The Xa Method.......... ..... .......... ............ 17

2-3. LCAO Methods for Extended Systems................. 21

2-4. Fitting Methods for Thin Films.................... 24

2-5. Self-Consistent Solution of the Secular

Equation.......................................... 34

Ill. COMPUTATIONAL METHODS ..................................... 38

3-1. Choice of Basis Sets.............................. 38

3-2. Overlap and Kinetic Energy Integrals.............. 39

3-3. Fitting Function Integrals........................ 41

3-4. Numerical Integration Techniques.................. 44

IV. RESULTS ................................................... 47

4-1. Atomic Hydrogen Monolayer......................... 47

4-2. Atomic Beryllium Monolayer ........................ 65

V. SUMMARY ................................................... 78

APPENDICES

A. HERMITE-GAUSSIAN FUNCTIONS AND INTEGRALS.................. 82

B. MULTIPLE EXPANSIONS OF COULOMB INTEGRALS................. 92

C. RECIPROCAL LATTICE EXPANSIONS......................... 102

D. NUMERICAL INTEGRATION GRID............................113

BIBLIOGRAPHY ................................................... 120

BIOGRAPHICAL SKETCH............................................ 124

Abstract of Dissertation Presented to the Graduate Council

of the University of Florida in Partial Fulfillment of the

Requirements for the Degree of Doctor of Philosophy

AN LCAO LOCAL DENSITY FUNCTIONAL APPROACH TO

SURFACE ELECTRONIC STRUCTURE CALCULATIONS

By

JOHN WALLACE MINTMIRE

DECEMBER 1980

Chairman: John R. Sabin

Major Department: Physics

A new method for calculating the electronic structure of thin films

is presented using the local density functional formalism in the linear

combination of atomic orbitals (LCAO) approximation. The self-consistent

solution of the resulting secular equations is made tractable through the

use of fitting procedures for approximating the charge density (for the

Coulomb potential) and the cube root of the charge density (for the ex-

change potential) as linear combinations of fitting functions, allowing

the construction of a more efficient (in terms of computational effort)

procedure than would otherwise be possible. Hermite-Gaussian functions

comprise the basis sets, with different sets of Hermite-Gaussian functions

used for the orbital basis set, for the charge density fitting functions,

and for the functions used to fit the cube root of the charge density.

Results are presented for two systems, the atomic hydrogen monolayer and

the atomic beryllium monolayer, which demonstrate the feasibility and

applicability of this approach. Techniques are discussed for improve-

ment of the computational procedure in order to enable calculations to

be made on larger and more complex systems than those presented in this

dissertation.

CHAPTER I

INTRODUCTION

1-1. General Remarks

The study of the electronic properties of solid surfaces is becoming

increasingly more important as techniques dependent on surface properties

(such as catalysis techniques and ultra-thin-film devices) occupy more

extensive roles in modern technology. Experimental interest in these

systems has expanded greatly in recent years, in large measure due to

improved vacuum techniques over the past two decades. The theoretical

study of the electronic properties of surfaces thus faces a challenge to

develop new methods for understanding the behavior of these systems.

Many experimental methods for studying surfaces have been developed

or improved in the last decade using a multitude of probe techniques.

Electron beams impinging upon the surface of interest are used in Auger

Electron Spectroscopy (AES), Low Energy Electron Diffraction (LEED), and

the Scanning Electron Microscope (SEM) to yield information about the

electronic structure and nuclear geometry of the surface, while Electron

Spectroscopy for Chemical Analysis (ESCA or XPS) and Ultraviolet Photon

Spectroscopy (UPS) use photons (X-ray and UV respectively) as a probe to

provide additional information about the electronic structure of surfaces

(for a historical review of these methods, see Lee, 1977).

In this work we propose a new computational approach to calculate

the approximate electronic two-dimensional band structure and total energy

of a slab having finite thickness in one direction and two-dimensional

periodicity in the plane normal to the direction of finite thickness.

The slab system will thus possess both a surface and the periodicity in

directions parallel to the surface characteristic of the bulk solid.

Thus by increasing the thickness of the slab, the properties of the semi-

infinite solid may be approached as a limit.

Methods will be described in this work for solving the secular

equation resulting from a Linear Combination of Atomic Orbitals (LCAO)

tight-binding approach to the Xa formalism (for a review see Connolly,

1976). The computational difficulties arising from the evaluation of

matrix elements involving pl/3(r) (where p is the electronic charge den-

sity) are treated by fitting p1/3 to a linear combination of two-

dimensionally periodic fitting functions in a manner equivalent to that

of Dunlap (Dunlap, Connolly, and Sabin, 1979a and 1979b) for molecular

systems. In order to reduce the number of computationally expensive

integrals which must be evaluated, the electronic charge density p is also

fit to a linear combination of a different set of two-dimensionally per-

iodic fitting functions. A mathematical formalism is derived for fitting

p based on the procedure introduced by Dunlap (Dunlap, Connolly, and

Sabin, 1979a) and Mintmire (1979) for molecular systems which finds the

variational extremum of an approximate electron repulsion energy func-

tional of the fit to p. This approximate electron repulsion energy

functional is bounded by the electron repulsion energy of p with itself

and is equal to the electron repulsion energy of p with itself to second

order in the difference between the exact density p and the fitted density.

The formalism derived by Dunlap and Mintmire in the above mentioned

references for the molecular case unfortunately cannot be used directly

for extended systems due to the long range nature of the Coulomb inter-

action resulting in certain integrals tending toward infinity in the

limit as a finite molecular system increases in size to approach the two-

dimensionally periodic slab. Thus changes were introduced into the

original formalism to yield a scheme equivalent to the original molecular

density scheme (Dunlap, Connolly, and Sabin, 1979a), but with proper

limiting behavior.

This new computational scheme is tested on two systems: the atomic

hydrogen monolayer and the atomic beryllium monolayer. The results for

the hydrogen monolayer demonstrate the feasibility of the method and

illustrate some of the effects of changes in basis sets and in the number

of points considered in the surface Brillouin zone. The atomic beryllium

results further demonstrate the feasibility of the method and selected

results are compared with experimentally determined results such as the

work function and Auger spectra of the beryllium surface.

1-2. Historical Background

Early calculations of the electronic structure of surfaces (excluding

cluster calculation approaches) normally fell into one or two categories:

1) treating the system as an electron gas and neglecting the

electron-nuclear lattice interactions (by assuming a homogeneous

positive charge density), such as the results by Lang (1973)

using a jellium model with a step potential to simulate the

potential at the surface; or

2) band structure calculations taking into account the interactions

between electrons and the nuclear lattice, although these

calculations almost universally were not self-consistent and

frequently were semi-empirical methods using experimental data

or bulk band structure data.

We wish to review briefly the history of calculations using the band

structure approach (with self-consistent procedures) over the last several

years in order to see how our proposed method compares with previous

work. One of the earliest calculations to report self-consistent calcu-

lations for the electronic structure of a thin film system was by

Alldredge and Kleinman (1972). These two workers calculated the band

structure of a 13 atom thick layer lithium system, taking into account

the effects of the nuclear lattice, using a non-self-consistent plane

wave solution which was later extended to a self-consistent scheme for

calculations on the same lithium system (Alldredge and Kleinman, 1974).

The lack of translational symmetry in the direction normal to the surface

was treated by forcing all orbitals to vanish at planes located some

large distance from the surface (Alldredge and Kleinman used planes lo-

cated three times the layer separation away from the exterior nuclear

site.), equivalent to invoking an impenetrable barrier for electrons

away from the surface. Core states were not treated directly in this

method, but core effects were included in a pseudopotential model.

Applebaum and Hamann (1972) introduced a model of the semi-infinite

surface by requiring all orbitals to match the bulk orbitals as boundary

conditions on one side of a thin film system. Model potentials including

core effects were used to generate orbitals in the Laue representation

(von Laue, 1931), where a Fourier expansion in the parallel coordinates

is combined with the coordinate representation in the normal direction,

using a numerical grid for the coordinate representation component of

all evaluated functions. This approach was used primarily to study re-

laxation effects on the (111) surface (Appelbaum and Hamann, 1973) and

the (100) surface (Appelbaum, Baraff, and Hamann, 1975a and 1975b) of

silicon.

Cooper (1973) reported calculations on the copper (100) monolayer

using a muffin-tin version of the KKR method adapted to two-dimensionally

periodic systems. The lack of translational symmetry in the surface

normal direction was treated with an impenetrable barrier in a similar

manner to that of Alldredge and Kleinman (1972). Various refinements to

this basic approach were introduced by Kohn (1975) and Kar and Soven

(1975) who simultaneously reported muffin-tin KKR approaches that did not

require the use of an impenetrable barrier located away from the surface.

Krakauer and Cooper (1977) further adapted this approach into a local

combination of muffin-tin orbitals (LCMTO) approach incorporating non-

muffin-tin corrections to the potential.

Gay, Smith, and Arlinghaus (1977) introduced a self-consistent band

structure method in an LCAO approach using Gaussian basis functions.

Their approach utilized an effective potential constructed from two terms:

1) the initial potential generated by overlapping atomic potentials;

and

2) differences between the effective potential at each iteration

in the self-consistent process and the initial potential cal-

culated as a Fourier expansion.

This approach has been very successful in describing the band structure

of transition metals, leading to computed results for a nine layer model

of the copper (100) surface (Smith, Gay, and Arlinghaus, 1980) and a

nine layer model of the nickel (100) surface (Arlinghaus, Gay and Smith,

1980). This method unfortunately appears not to admit a straightforward

method for evaluating the total energy of the system due to the nature

of the Fourier transform techniques.

Wang and Freeman (1978) introduced another LCAO approach for thin

films based on the Discrete Variation Method (DVM) of Ellis and Painter

(1971) where the orbitals are expressed in a numerical basis set. The

charge density used for generating the secular matrix of each iteration

in this self-consistent procedure is constructed by fitting the current

true charge density to a charge density resulting from superimposed over-

lapping spherically symmetric atomic charge densities. This method has

been used to calculate the band structure of one, three, and five layer

models of the nickel (001) surface (Wang and Freeman, 1979) and a nine

layer model of the nickel (001) surface (Wang and Freeman, 1980). No

total energies have been reported using this method, although no inherent

difficulties appear to prevent this approach from yielding total energies

in a straightforward fashion.

In addition to the method introduced by them mentioned previously,

an alternative self-consistent LCAO procedure using Gaussian basis func-

tions has been introduced by Appelbaum and Hamann (1978), with calculated

results for the copper (100) monolayer. Their approach is similar to that

derived in this work in the use of fitting techniques for the charge den-

sity p and for p1/3. Appelbaum and Hamann fit the charge density with a

least-squares fitting procedure to a periodic sum of spherically symmetric

Gaussians using "floating" Gaussians at sites other than the nuclear

centers in order to describe asymetric components of the charge density.

Using the same basis of spherically symmetric Gaussians as used for

fitting the charge density, Appelbaum and Hamann then fit the sum of the

local density functional exchange potential and the total Coulomb poten-

tial minus a screened nuclear attraction term, which is treated outside

the fitting procedure. This approach has been used to calculate the

band structure of the eleven layer model of the titanium (0001) surface

(Feibelman, Appelbaum, and Hamann, 1979) and the same titanium system

with a chemisorbed hydrogen monolayer (Feibelman and Hamann, 1980) and

with a chemisorbed nitrogen monolayer (Feibelman and Himpsel, 1980).

Recent calculations (Feibelman, Hamann, and Himpsel, 1980) report cohesive

energy curves for the eleven layer model of the titanium (0001) surface

with a chemisorbed hydrogen monolayer at various distances from the ex-

terior layer of titanium. Unfortunately complete details as to how the

total energy terms were evaluated in this formalism have not been pub-

lished at this time.

The formalism described in this work is similar to the techniques

discussed by Feibelman, Appelbaum, and Hamann (1979); the primary dif-

ferences arise from the fitting techniques utilized. In our approach

a basis set of Gaussian functions (not restricted to spherically symmetric

Gaussians) is used to fit the pl/3 portion of the local density functional

exchange potential by using least-squares fitting methods. The charge

density is then fitted using the variational methods mentioned in Section

1-1, which results in the evaluated Coulombic energy portion of the total

energy agreeing with the Coulombic energy resulting from the non-fitted

density to second order in the error in the charge density introduced by

the fitting procedure. We feel this is a more rational scheme than a

least-squares fitting method for fitting the charge density when the

fitted charge density will be used to construct the Coulombic potential

due to the electron charge and to evaluate the Coulombic contribution

to the total energy.

1-3. Organization of Dissertation

There are two fundamental topics to cover in this dissertation.

The first topic is a basic description of the electronic band structure

approach. The mathematical formalism of our approach is presented in

Chapter II, followed by a brief description of the computational tech-

niques used to implement the formalism given in Chapter Ill. The results

are the second fundamental topic to be discussed. Computational results

for the atomic hydrogen monolayer and the atomic beryllium are presented

and discussed in Chapter IV, with a summary of our method and results

in Chapter V.

The appendices at the end of this dissertation discuss in extended

detail several topics relating to the computational procedure. Appendix

A defines the Hermite-Gaussian functions used in the work and integrals

involving these functions are reduced to forms which may be evaluated

using algorithms described in this appendix and later appendices. The

multiple expansion techniques used in evaluating Coulomb integrals are

discussed in Appendix B, while equations necessary for the evaluation of

integrals using reciprocal lattice expansions are derived in Appendix C.

Finally, Appendix D describes the three-dimensional grid used for the

numerical integration algorithms which are part of the computational

fitting procedure.

CHAPTER II

MATHEMATICAL FORMALISM AND THEORY

2-1. Symmetry Properties of Thin Films

As mentioned in Chapter I, the physical systems of actual interest

were studied using an ideal periodic thin film model. Given an ideal

three-dimensional periodic lattice of nuclei, the thin film may be vis-

ualized as that portion of the three-dimensional lattice enclosed between

two parallel planes of finite separation as illustrated in Figure 2-1.

A more formal definition of the systems under consideration is expressed

by the following two conditions on the positions of the nuclei in the

nuclear lattice:

1) There exist two independent vectors R and R such that a

translation of coordinates by R = mR + nR where m and n are

integers, results in an equivalent lattice;

2) and there exist two planar surfaces, parallel to the plane of

periodicity containing RI and R with finite separation such

that all of the nuclear sites are contained within the region

enclosed by the two parallel planes.

The first condition requires two-dimensional periodicity while the second

condition requires that the thin film indeed has finite thickness.

The set of all possible translations R,

R = mR + nR (m and n integers) (2-1)

1 -2

Figure 2-1. Schematic representation of a side view of a two-dimension-

ally periodic slab with three layers.

defines a two-dimensional geometric lattice of points, the planar Bravais

lattice. There are five possible types of planar Bravais lattices, which

have different point group symmetries associated with them in addition to

the translational symmetry, as illustrated in Figure 2-2. An elementary

treatment of these systems is given in Blakemore (1974).

Once a planar Bravais lattice is chosen, the treatment of symmetry

is similar to that of the three-dimensional crystalline lattice as des-

cribed by Tinkham (1964) or by Slater (1972a). A unit cell may be con-

structed as shown in Figure 2-3 in a manner analogous to that of the

three-dimensional unit cell. Since the thin film is actually a three-

dimensional system while the planar Bravais lattice is only a two-

dimensional geometrical construct, the thin film unit cell is a parallel-

opiped which is finite in directions parallel to the plane of periodicity

but infinite in extent in the directions normal to the plane of period-

icity. The arrangement of nuclear centers inside each unit cell is

unrestricted except for the previous constraint that the thin film have

finite thickness. Therefore the complete arrangement of nuclei need not

have the full point group symmetry of the Bravais lattice.

Thus the symmetry group of the total system (or rather the plane

group, since we are currently considering the set of two-dimensional

symmorphic symmetry operations) is a subgroup of the plane group of the

Bravais lattice. Using an approach similar to that of Koster (1957) the

plane group consists of symmetry operators of the form {P R} which cor-

respond to a point group coordinate transformation followed by a trans-

lation of coordinates:

{PIR} r = Pr + R

(2-2)

0 T (,-od -.tro. u..

A Th"tee-old lro- o als

0 Four- fold torsi*2

O Sn-told nNKM 0

OblIque Littice lallbl. a 90

b7

S-----J

Rectangug., Labtc. Ia|ibI. a QO

0r

r I

Ctfrtcred Rectantulu a cosffa/2

I M-roymnroiry pfw

+ Omhogol mi pfnoc

* Mi.r- pl.,o. -.ty .5*

./i

0--

ii

The svmrtry 1s tO wsme

to ihY mnim m svmmi n .c

of any otlque lrsco.

Figure 2-2. The five two-dimensional lattices and their symmetries.

(From Solid State Physics by J. S. Blakemore. Copyright (c)

1969 by W. B. Saunders Company. Reprinted by permission of

Holt, Rinehart, and Winston.)

Figure 2-3. Illustration of parallelopiped unit cell for the two-dimen-

sionally periodic monolayer. Note that the unit cell extends

to infinite distance in both directions parallel to the z

axis.

The operator {(PIR can be seen to have the inverse operator {PIR}-

such that:

{PIR} = {P -P1 R} (2-3)

Given these definitions the operation of a symmetry operator upon a

function f(r) may be defined:

{PIR} f(r) = f({PIR} r) (2-4)

Bloch's theorem applies to the two-dimensionally periodic thin film

in a manner analogous to that for the three-dimensional crystalline case.

From group theory (Tinkham, 1964) one can show that for an operator M

A A

which commutes with all the translation operators R (where R denotes the

subgroup of translation operators of the plane group given by {EIR}, if

E denotes the identity point group operator) the eigenfunctions of M can

be chosen as members of the irreducible representations of the translation

group. These irreducible representations may be labelled with wave

vectors k that lie parallel to the plane of periodicity such that if

the function D(r; k ) belongs to the k irreducible representation,

then:

R i(r; k,,) = (r R; k,1) = exp(ik /.R) O(r; k11) (2-5)

In a manner similar to the three-dimensional crystalline case, the

set of wave vectors k is not a set of unique labels for each individual

irreducible representation. For example, if the vector c is normal to

both primitive lattice vectors R and R (e.g., c = R x R ), then two

S -reciproca2 as:

reciprocal lattice vectors may be defined as:

R 2 2( ( x c) (2-6)

-1 -2 x C2

and

2 = 2r (R x c) (2-7)

-2 R2. (R x c)T 1

These reciprocal lattice vectors have the property that:

R. -K. = 27r 6.. (2-8)

-I -J IJ

Thus if we define a general reciprocal lattice vector K such that:

K = hK + kK (2-9)

-- -1 -2

then any two wave vectors k and k' related by:

k' = k + K (2-10)

label the same irreducible representation since:

(r + R; k') = exp[i(k/ + K)*R] D(r; k') (2-11)

= exp[ik *R] P(r; kl )

Thus we may restrict our choice of wave vectors to the two-dimensional

central Brillouin zone of wave vectors k closer to the origin than to

any other point described by an arbitrary reciprocal lattice vector K / 0.

With the Brillouin zone defined we may restate Bloch's theorem

in a more straightforward form for use in solving finite secular matrix

problems. The matrix elements of an operator M, which commutes with all

translation operators R, will be zero for matrix elements between functions

belonging to different irreducible representations labelled by k and k'/.

Recalling that the symmetry operators are unitary, or:

= (2-12)

We next consider the matrix element

= <(r; k/) R- R M .(r; kl)> (2-13)

=

= exp[i(k' k l ).R]

Since R may represent any lattice point and by the hypothesis we are

given that k' k K, then the matrix element must equal zero:

<.(r; kl/)I M j(r; k //)> = 0 (2-14)

A star of the wave vector k may be defined as the set of wave

vectors generated by the operation of each symmetry operator in the point

symmetry subgroup of the plane group of the thin film system on the wave

vector k Since the wave vectors k are merely labels for the irre-

ducible representations of the translation group, it is understood that

if a symmetry operation yields a wave vector k" outside the central

Brillouin zone then the wave vector k' = k" + K equivalent to the wave

vector k" is used as the member of the star. Thus it can be shown that

the matrix elements of an operator M which commutes with the point group

symmetry operators P (where P denotes (PIO}) have the property that:

< .(r; k//) M .(r; = k )> = < .(r; k' p- P M .(r; k //)> (2-15)

=

I- -I/ j -II

We notice that the new function P(i.(r; k1/) belongs to the irreducible

representation labelled by Pk since if R is a general translation

operator:

RP 4.(r; ki ) = R (P- r; kl) (2-16)

= .(P [r R]; k //)

= exp[-ik1, P R] 4 i(P r; k,,)

= exp[-i(Pk )-R] P P.(r; kI)

From group theory we recall that the eigenvectors of an operator M

which commute with the operators of the symmetry group will be degenerate

if a wave vector star possesses more than one unique member, since each

irreducible representation k in the wave vector star will contain one

degenerate eigenfunction. Equation 2-16 implies that this relationship

is also true for the finite secular matrix problem if the basis sets are

such that if the basis set for irreducible k contains a function

<.(r; k ), then the basis set for the Pk irreducible representation

must contain the function P (r; k ).

2-2. The Xa Method

The Xa method is an approximate method for calculating the electronic

structure of atoms, molecules, and extended systems using an approach

originated by Slater (1951). The Born-Oppenheimer approximation (1927)

is used to simplify the problem of finding the eigenvalues of the elec-

tronic Hamiltonian:

[H(R) E(R)]M(r.) = 0 (2-17)

where (using Hartree atomic units throughout):

192 1 m mn } (218)

H(R) = f.--- Z R-R } 2- I8)

1 2 r. j> r.. _r-R Rnm n -Rn-

i -i j> Ij m- m>n-m -n

The coordinates R refer to nuclear coordinates and the coordinates r.

-m -I

refer to electronic coordinates.

The Xa method originated as a means of approximating the Hartree-

Fock model by replacing the non-local exchange terms with a local func-

tional of the electron density (the origins of the Xa method are discussed

at length in Slater (1974)). A different interpretation of the local

exchange potential may be made using an approach based on the local density

functional formalism of Kohn and Sham (1965). This local density func-

tional formalism is based on the theorem by Hohenberg and Kohn (1964).

The Hohenberg-Kohn theorem states that for a Hamiltonian of the general

form

H = E{-- V + v(r ) + Z } (2-19)

2 r. ( r..

i -I j>I Ij

there exists an energy functional G[p] such that if p(r) is the electron

density resulting from Y(r.), then

3nr

G[p] = fd nr '(r.) [H Ev(r.)] Y(r.) (2-20)

Furthermore if p0(r) corresponds to the ground state 'Yo(r.) of the Ham-

iltonian H, then the total energy functional:

E[p] = G[p] + fd3r p(r) v(r) (2-21)

has a variational minimum about the electron density po(r).

If instead of considering just the electron density p(r) but considering

the first order density matrix p(r; r') which may without loss of gener-

ality be assumed to be of the form:

p(r; r') = E n. .(r') < (r) (2-22)

where the functions 4. form an orthonormal basis and the occupation num-

bers have the property 0 n. n 1, then we may define a new energy func-

tional E xcp]:

E [p] = G[p] fd3r fd3 r' 6(r r')[- -V2] p(r; r')

xc 2 r

p(r) p(r')

-1 fd3r fd r' -- -- (2-23)

IL -

Ir r'

where p(r) = p(r; r). The second term on the right-hand side of Equation

2-23 corresponds to the expectation value of the kinetic energy of the

first order density matrix p(r; r') and the third term corresponds to the

Coulomb interaction of the charge density p(r) with itself. Thus we see

that the energy functional E [p] will correspond to the exchange and

xc

correlation energy components of the total energy.

The local density functional formalism is based on the approximation

of the exchange correlation energy functional, which has an unknown

analytic form, by a prescribed functional of the form:

E xcp] = fd3r {pt(r) V xc(pt) + p+(r) V (pf)} (2-24)

where V (r) is a function of the electron density of spin up for pt and

of the electron density of spin down for p+. The potential Vx (r) is

most properly denoted as the exchange-correlation potential. However,

for brevity and historical continuity we shall refer to V (r) henceforth

as the exchange potential.

In the Xa model the exchange potential V is chosen to be of the

XC

form introduced by Slater (1951) with the addition of an adjustable

parameter a with a value between 2/3 and I (Slater and Wood, 1971):

1/3

V (pt) 9 1/ (2-25)

For convenience we shall henceforth refer to an exchange potential operator

V which when operating on p(r) has the form:

XCa

Vxap = V X(pt)p+ + V (p4)pt (2-26)

This set of approximations yields an energy functional of the form:

S3 V p(r) p(r')

E[p] = fd3r fd r' {6(r r') p(r; r') + r

2 r Ir r'

Z Z Z

+ fd r {V mp(r) + E mn (2-27)

m m>n m n

where the total energy includes the nuclear repulsion term introduced by

the Born-Oppenheimer approximation.

The orbitals i.(r) are chosen by assuming as an approximate form the

linear combination of a finite set of basis functions:

i.(r) = E c i.(r) (2-28)

m

and the energy functional E[p] is minimized using the variational principle

subject to the constraint that the orbitals 4. are orthonormalized. This

procedure yields the familiar secular equation:

Z H c = Z S c e. (2-29)

mn ni mn ni i

n n

where c. are the orbital eigenvalues and

Smn = fd3r m(r) n(r) (2-30)

Hmn = fd3r (r) H (r) n (r) (2-31)

The one-electron Hamiltonian Heff(r) is given by:

S2 Z p(r') 4

H eff(r) = -V2 + E + fd3r' +- V (p) (2-32)

m r R - 3 X

where p corresponds to pt- if and n correspond to spin up, and p

o m n o

corresponds to p+ if im and yn correspond to spin down.

2-3. LCAO Methods for Extended Systems

For infinitely extended systems the limit of the total energy of a

finite section of the infinite system as the size of the finite section

is systematically increased is not a useful definition of the total energy

of the periodic system. It is easily seen that the total energy defined

in this fashion will in most cases increase linearly with the cross-

sectional area of the finite section (for two-dimensionally periodic

systems as defined in Section 2-1) and thus will not have a convergent

limit. The total energy per unit cell is a far more useful quantity that

is usually convergent in the above mentioned limiting process for systems

that possess zero total charge. In the Xa model the functional yielding

the energy per unit cell may be expressed:

E[p] = f d r fd3r' 6(r r') V2 p(r; r')

2 r -

1 3 P (r) p(r') 3 Z

+ fdr fd3r' d3r (r) E m

2 r' p R R

R m -m -

ZZ

+ d3r V (r) + E I mn (2-33)

x -h ' | K| (2-33)

r xa 2 nR R R

R m,n -m -n -

where the primed summation over the indices m and n indicates that the

terms corresponding to m = n when R = 0 are to be omitted. The inte-

grations over the region S denote an integration over one three-dimensional

unit cell (as defined in Section 2-1).

For computational convenience we assume Born-von Karman (1912)

periodic boundary conditions, such that all eigenfunctions of periodic

operators are periodic over translations of the general form M R + N R ,

1 -2

where M and N are integers and R and R are the primitive lattice vectors

1 -2

for the system of interest. Periodic boundary conditions reduce the

problem involving an infinite number of electrons to the problem of a

finite number of electrons contained in M by N unit cells, such that the

orbitals obey periodic boundary conditions. This yields a procedure

which is effectively a cluster calculation on an M by N unit cell system

with an external environment outside being a periodically repeated version

of the internal system rather than empty space.

Let us approximate the orbitals 4.(r) by a linear combination of a

finite number of orbitals m (r). The spin-up one-electron effective

Hamiltonian for this procedure may be expressed,

p(r') Z 1/3

Ht (r) = -- v2+ fd3 r' 7 EE m 3a 3p]

- R m m -

(2-34)

with the spin down effective Hamiltonian having an analogous form. From

equation 2-14 we see that the logical choice of basis functions will be

Bloch functions m (r; k,/), or now labelling functions with the wave

vector k//:

p(r; r') = n.(k/ ) i.(r; k d) .i(r; k ) (2-35)

i k i

4.(r; r ) = E cim ) (r; k ) (2-36)

m

We see that the Born-von Karman periodic boundary conditions require that

the orthonormalization condition be expressed:

fS1d r i(r; k' ) .(r; ,/) = 6 6

k // i k, (2-37)

where the integration region Q' corresponds to the region defined by the

M by N unit cells. The periodic boundary conditions also restrict the

choice of wave vectors to a discrete set, since

P(r + MR + NR2; ki) = exp[ik -(MR1 + NR )] i.(r; k/ ) (2-38)

= i.(r; k /)

This implies

k M + n K (m and n integers) (2-29)

are the only allowed choices of k,,, in addition to the condition that we

only consider wave vectors in the central Brillouin zone.

The Bloch basis functions p.(r; k//) may be constructed from localized

functions U.(r):

ki(r; 11) = E exp[ik ,.R] U.(r R) (2-40)

V R /

This yields a secular equation of the form

E mn(k c (k) = S (k ) c (k //) (2-41)

mn (/ n -n/ ni -//'

n n

where

S (k E) = exp[ik//*R] fd3r U (r+R) U (r) (2-42)

mn R/-- m--- n -

R

n (k ) = E exp[i1/.R] fd3r U(r+R) Hff(r) U (r) (2-43)

mn -/ R- m nf--

2-4. Fitting Methods for Thin Films

The effective Hamiltonian in Equation 2-41 contains a term propor-

tional to p ; one cannot usually express integrals involving this term

in an analytic closed form. One way of alleviating this difficulty for

molecular systems is the approximation of the functions p (r) by a linear

combination of fitting functions G m(r), henceforth referred to as the

exchange fit,

/3(r) E gm G (r) (2-44)

m

m

This approach to solving the problems involving the analytic behavior of

1/3

p was originally suggested by Sambe and Felton (1975) with extensions

to the formalism and computational method later introduced by Dunlap,

Connolly, and Sabin (1977) and Mintmire (1979). In these previous works

the total charge density p(r) was also approximated using a linear com-

bination of a set of charge fitting functions F (r), this procedure being

henceforth referred to as the charge fit,

p(r) = fm F (r) (2-45)

m

The rationale behind approximating p(r) is to reduce the number of analytic

integrals which must be computed, thus achieving a savings in computa-

tional effort.

For extended systems the fitting functions m (r) and Fm (r) are com-

posed of periodic sums of localized functions G (r) and F (r), (where a

m m

localized function has the property that lim r3 G (r) = 0), that may be

r+-> m

expressed:

m (r) = E G (r R) (2-46)

R- m m- -

R

m (r) = F m(r R) (2-47)

R

The exchange fitting coefficients gm are chosen using a least squares

fitting procedure in which the set of conditions

L [ f d3r{pl/3(r) g (r)}2 = 0 (2-48)

i m

yield the matrix equation

EG g = x (2-49)

mn n m

n

where

G = E fd3r G (r) G (r R) (2-50)

mn m n -

R

and

xm = dr p/3 r) Gm (r) (2-51)

=fd3 r p/3 (r) G m (r)

The matrix elements G may be evaluated in analytic form if we choose

mn

the functions Gm (r) to be some analytic function such as Hermite-Gaussian

functions. The integrals x however, require numerical integration for

most cases where only the charge density is known as a linear combination

of analytic functions. This procedure does allow a reduction of the

total number of numerical integration required compared to a process

not using fitting by an approximate factor of N /Nx, where N is the number

of orbital basis functions and N is the number of exchange fitting func-

tions G (r).

m --

The total charge density p(r) is also approximated using a linear

combination of analytic functions such as Hermite-Gaussian functions.

The rationale for fitting the charge density as well as p1/3 is that the

fitting process will reduce the total number of computationally time

consuming integrals necessary for computing the Coulomb contributions to

the total energy and the secular matrix by an approximate factor of N2 /N

where N is the number of charge fitting functions F (r). Since the value

of N is typically the same order of magnitude as N, approximating the

density may produce a substantial saving in computational time.

Before explaining the procedure for choosing the coefficients fm'

let us digress for a moment to the problem of approximating charge den-

sities in the general case using the methods introduced for molecular

claculations by Dunlap (Dunlap, Connolly, and Sabin, 1979a) and by

Mintmire (1979). Consider a prototypical molecular charge density p(r)

-3

such that p decreases more rapidly than r in the limit of large r.

Then define the electrostatic interaction U of a charge density with it-

self to be U = I [p1p], where

[plp] = fd3r fd3r' ( p2(r) (2-52)

If the specified charge density p is approximated by another density p

which is dependent upon several parameters which may be varied independ-

ently, and if

p(r) = p(r) + Ap(r) (2-53)

then

[pip] = 2[plp] [P1p] + [Ap|Ap] (2-54)

Let us define an approximate electrostatic interaction energy U in terms

of the approximate density p and the exact charge density p:

U = [p ] 1 [ i] (2-55)

Since the difference AU = U = [Ap|Ap] is greater than or equal to

zero for any choice of p, a reasonable criterion for choosing p is the

minimization of AU. This constraint will be satisfied if the variation

of AU with respect to allowed variations of p is zero:

6AU = 6U = 0 (2-56)

leading to the relationship expressed in terms of 6p:

6U = [p|6p] [pI6p] = 0 (2-57)

We see that if a variation of the form 6p = p de is possible within

the range of variations allowed to p and de arbitrary, then the additional

relationship

{[pl ] [p10]} de = 0

is valid. This equation implies the equality

[p1p] = [I]

which shows that for the choice of p satisfying Equation 2-57, then

U = [PP] =

2 2

(2-58)

(2-59)

(2-60)

For p expressed as a linear combination of fitting functions Fm (r):

p(r) = E f F (r)

-- m m

m

(2-61)

we may restate Equation 2-57 in terms of the quantities f and F (r):

m m

6U = E {[pIp-] [ p- I-]} 6f,

m m m

(2-62)

E {[PIF ] E f [F IF 1m 6f = 0

m n

Since the variations 6f are arbitrary and independent, each term in the

summation over the index m must equal zero. Rearranging terms yields:

E[F Fn fn = [F MP]

n

E B f = a

mn n n

n

where

(2-63)

(2-64)

B = [F IF ]

mn m n

am = [mP

(2-65)

With the proper choice of fitting functions F (r) and form of the charge

density p, these integrals may be evaluated analytically.

Extension of these methods to extended systems requires a few mod-

ifications. The fitting functions are now the periodic fitting functions

F (r) formed from the periodic sum of the localized functions F (r) as

m m

stated in Equation 2-47. The charge interaction expression [pl|P2] must

also be redefined in order to avoid dealing with infinite terms. The

most convenient definition of this term, if p, and p2 are periodic extended

charge densities, is the electrostatic interaction per unit cell:

= p (r) p2(r')

[P1jP2] = fd3r fd3r' (2-66)

The charge interaction term [pl|p2] is still non-negative for p = p2

and p1 having no singularities. In addition, this expression is still

symmetric with respect to the order of pl and p2 (both of which are

periodic), since

P 2 = f fd3r' p (r) p 2(r')

pr r r' (2-67)

d3r d3 P p (r) p 2(r' R)

= Z fd r fd r'

R Ir r' + RI

Sdr f d3r (r + R) p2(r')

= E fL d r f d d-r'

R Ir+ R r'l

3r 3 P 2(r') p I(r)

= f drf d = [p21p1]

r r'-

This expression for the charge interaction is also convenient because

it satisfies two conditions:

1) for m(r) a periodic function formed from the periodic sum of

of localized functions F (r) as in Equation 2-47, then

fd3r F (r) F (r R)

[F jF ] = E fd r fd 3r m n (2-68)

m n R L 'l(2-68)

which is analytically integrable in closed form for F (r) being

a Hermite-Gaussian function;

2) [pIIP2] is finite if and only if one or both of the charge

distributions pl and p2 satisfy the charge neutrality condition:

Jfd r p(r) = 0 (2-69)

The first condition is evident upon inspection. The second condition

follows from consideration of a limiting procedure which begins with a

finite assembly of identical unit cells. The initial assembly is con-

structed with no embedded voids. The limiting process then consists of

arbitrary enlargement of the finite assembly, with the infinitely periodic

(in two-dimensions) slab the ultimate result. Any such finite assembly

has a periodic charge density which satisfies the explicit constraint on

asymptotic behavior required by the definition of the molecular charge

density interaction [p 1 2] of Equation 2-52. It then follows from

Equation 2-52 that the charge interaction for such a finite assembly with

the charge density belonging to a single unit cell is

3 d3 1p (r) p 2(r')

[PlP2] = f d r fd r' (2-70)

R r r' + Rr

where the primed summation indicates a finite number of terms. In order

to carry through a physically reasonable limiting process the single

unit cell is conceived as being the most nearly central unit cell of the

finite assembly. Similarly the origin of coordinates is associated, for

convenience, with the central unit cell. Then for each term of Equation

2-70 save for the one corresponding to R = 0 the contribution may be

exactly as a multiple expansion (see Appendix B for a more detailed

discussion). If the charge neutrality condition is not satisfied then

the leading term in the multiple expansion will be a monopole-monopole

-l

interaction term proportional to R For a two-dimensional lattice of

vectors R an infinite sum of monopole-monopole interactions is not con-

vergent. Therefore in the limit of an arbitrarily extended slab, the

charge neutrality condition is necessary if the charge interaction is to

converge to finite limit. Similarly if the charge neutrality condition

is satisfied, the leading term in the multiple expansion for given R

-2

is a monopole-dipole term proportional to R2. It is shown in Appendix

B, however, that this multiple term has zero total contribution to the

charge interaction term for two-dimensional lattices with inversion sym-

metry. Thus the leading order term with a non-zero contribution to the

charge interaction term will be a term proportional to R 3, which leads

to a convergent sum in Equation 2-70 in the infinite limit under consid-

eration.

Thus the second condition leads immediately to the observation that

Equations 2-63 through 2-65 as derived are not usable in a direct manner,

since the density involved is purely electronic and manifestly not neutral.

However, Equation 2-65 is valid for the extended system so long as charge

densities neutral by cells are employed. For any finite system Equation

2-63 may be rewritten trivially as

[F p E f F ] = 0 (2-71)

m n n

n

irrespective of charge neutrality. In the limit described above, however,

the charge interaction expression becomes that of Equation 2-67 and the

necessary and sufficient conditions on charge neutrality come into play

with the result that if at least one of the F is not charge neutral,

then the right-hand factor in Equation 2-71 must be charge neutral,

where

fd3r p(r) = E f fd3r F (r) = Z f fd3r F (r) (2-72)

n n

or expressed differently

Sfm n = 1 (2-73)

n

where

nm = {fd3r F (r)} / {f>d3r p(r)} (2-74)

We can produce a matrix equation equivalent to Equation 2-64 if we intro-

duce a nuclear lattice charge pN to each electronic density, thereby

giving the requisite cell-by-cell neutrality, to the interaction term in

Equation 2-64 and then subtract out terms contained in Equation 2-71 and

analyze the identity

[Fm nmPNIP [Fm nmPN' fn(Fn PN) (2-75)

n

= [Fmp E fn F {l E fn}[Fm mPNPN] nm[PN - fnFn]

n n n

The first term of the right hand side of Equation 2-75 is zero according

to Equation 2-71. Since Equation 2-73 states that {l E fnn } = 0 and

n n

in addition we see that [Fm nmPN PN, is finite, then the second term of

in ad it o we see

the right-hand side of Equation 2-75 must equal zero. The term

[p Np Z f F ] must also be finite (again by charge neutrality) for

n n

thin films; therefore Equation 2-75 reduces to:

{E f [Fm nmPNn GN Fm nmN PN] = nm (2-76)

n

where

A = 1[NP E F n (2-77)

n

is an undetermined constant independent of the index m in Equation 2-76.

Restating Equation 2-76 in matrix notation yields:

E B f = a + An (2-78)

mn n m m

n

where

Bmn = [Fm nmIFn nn N (2-79)

am = [Fm nm N (2-80)

This set of equations is similar to the equations resulting from a

least-squares fitting procedure which includes a linear constraint equation.

Thus we may solve for A by rearranging Equation 2-78:

-1 -1

f {B }mn a = {B } n (2-81)

m mn n mn n

n n

Em f E n {B -} a =A3 nm {B I r28

nm m m mnn m mn n (2-82)

m m n mn

resulting in the value of

X = {1 Z E nm (B-}mn a } / {E E n {B } n } (2-83)

n m mn n

mn mn

The solution for the coefficients f is then given by

-1

f = E {B mn (an + n ) (2-84)

n

We see that the approximate electrostatic interaction energy defined in

Equation 2-55 now has the form:

~ 1 -l

U = E f a - E f {B } f (2-85)

n n 2 n mn n

n mn

= [pip] [pip] [pIP + I [pp 1 (2-86)

which is the electrostatic energy of the combined electronic and nuclear

charges with only the electron repulsion term being replaced by an approx-

imate interaction.

2-5. Self-Consistent Solution of the Secular Equation

The orbital coefficients cni (k //) are obtained by solving the secular

equation in Equation 2-41. Since the matrix elements H (k /) are de-

mn 2//

pendent upon the coefficients cni (k ) for which we seek a solution,

Equation 2-41 is solved using the familiar iterative procedure known as

the self-consistent field (SCF) method, described by Slater (1963) for

atomic and molecular systems. In the SCF method one starts with an

initial choice of c .(k ) to generate the matrix elements H (k ). The

n i -// mn -//

secular equation is then solved to yield a new set of values for the

coefficients cni (k ) which is used to begin a new iteration. This process

is repeated until some predetermined criterion of convergence for the

iterative process is satisfied.

Due to the introduction of the fitting procedures, the secular

equation defined in Equations 2-41 through 2-43 will not use the effective

Hamiltonian as it stands in Equation 2-34. Instead the effective Ham-

iltonian will be of the form (for the spin up effective Hamiltonian):

1 2 3 m a

H+ (r) = -- V + E f fd3 r E -

eff 2 m r r a R-

m Ra -a -

3a[ ]/3 E g G (r) (2-87)

,T m m

m

Thus given an initial choice of cni (k //) and n.(k//), the logical

order of each iteration will be:

1) Compute the fitting coefficients f and g

m m

2) Use the fitting coefficients to generate the secular matrix

defined in Equation 2-87;

3) Solve the secular Equation 2-41 to yield a new set of coefficients

c (k//) and occupation numbers n.(k,).

ni -// i -"//

Since this procedure differs from those of typical LCAO calculations in

the use of a fitting formalism, we present briefly a discussion of the

necessary equations for the above-mentioned procedure for each iteration,

beginning with an initial set of coefficients c .(k /) and occupation

numbers n. (k//).

First the charge density p of Equation 2-35 is used with the fitting

scheme of Equations 2-78 through 2-84 to compute the charge fitting

coefficients f The B matrix elements of Equation 2-79 are evaluated

m mn

directly while the values for a are constructed:

m

a = {E Z n.(k ) E c..(k//)c (k ) V (k )} C (2-88)

m i 1ji -// ki -//i jkm-/ m

ik/ j k

where

Vjkmk//) = R R exp[ik -R] fd3r U r + R) U (r)

Vjkmk/ R RI k k

F (r' RI) Z

x {fd r' n m Z R a M} (2-89)

and

C = {[fd3r F (r R) ZE E ab} (2-90)

m m r a Ri_ | T]lm R (9 R

R a --a a b -a -b -

The primed summation indicates that nuclear self-interaction terms for

which R is the zero vector and indices a and b denote identical nuclear

sites are to be omitted.

Next the cube root of the density is fit using Equation 2-49. The

G matrix elements are evaluated in the form given in Equation 2-50.

mn

The x integrals are evaluated using a numerical integration procedure.

For each numerical grid point r. the value of the density p(ri.) is generated

from the precalculated values of Bloch functions P.(r.; k/) and the

exchange functions G (r.):

P(r.) = E Z n (k ) I E c. (k, ) .(r.; k ) |2 (2-91)

j I -

which is then used to compute pl/3(r.). The numerical cube root of the

density is used then with the values of the exchange fitting functions

G (r.) to evaluate x integral numerically:

m --i m

x = Z w. pl/3(r.) G (r.) (2-92)

where the weights w. and the grid points r. are appropriately chosen

I --I

grid points for a three-dimensional numerical integration over the unit

cell volume Q.

The fitting coefficients fm and gm are then used to evaluate the

elements of the secular matrix defined by Equation 2-43. Using the

expression for the effective Hamiltonian of Equation 2-87, the matrix

elements for the spin up effective Hamiltonian have the form:

H. (k ) = T..(k//) + E f V (k )-3a[ ]11/3 X. (k ) (2-93)

,j I j m ijmml m ijm

m m

where

T..(k, ) = E exp[ik R_ fd U (r + R){- V2} IJ (r) (2-94)

IR // R // 1 2-

and

X.. (k//) = Z exp[ik// R] E fd3r U*(r + R) U.(r) G (r R') (2-95)

ijm / R R m

This form of the secular matrix is diagonalized and new occupation numbers

n.(k/) are computed, using the eigenvalues of the secular equation to

find the lowest eigenstates of the effective Hamiltonian to occupy. At

this point the agreement between the current orbital coefficients c ni (k/)

and those present at the beginning of the current iteration is tested for

convergence of the iterative sequence. If the convergence criteria are

not satisfied, a new iteration begins.

CHAPTER III

COMPUTATIONAL METHODS

3-1. Choice of Basis Sets

In any atomic, molecular, or solid state calculation of electronic

structure using a linear combination of a finite set of basis functions

as approximate orbitals, the particular choice of functions is a critical

part of the computational procedure. The first decision in choosing

basis functions that must be made is the type of functions to be used.

This work is based on the use of Hermite-Gaussian functions, which are

described in Appendix A. These functions are specific linear combinations

of the more familiar Cartesian Gaussians of the form x ymznexp(-ar2),

where 1, m, and n are nonnegative integers and the exponent coefficient

a (commonly denoted the "orbital exponent") is constant for all Cartesian

Gaussians which are components of the Hermite-Gaussian. That is, all

Hermite-Gaussians may be expressed as P(x,y,z)exp(-ar 2), where P(x,y,z)

is a polynomial of finite order of the coordinate variables x, y, and z.

The total number of Hermite-Gaussian functions in the basis and the

choice of the orbital exponents a. is also a critical choice. Fortu-

nately there is an extensive literature of results of atomic and molec-

ular calculations using Cartesian Gaussian functions that may be consulted

in making these decisions so that Cartesian basis sets may be converted

to Hermite-Gaussian basis sets. The conversion for most systems from

Cartesian Gaussian basis sets to Hermite-Gaussian basis sets is usually

trivial, since s and p type Cartesian Gaussian functions correspond to

equivalent Hermite-Gaussian basis functions. Cartesian d type functions

cannot be expressed as a single Hermite-Gaussian function; instead they

must be expressed as a sum of two Hermite-Gaussian functions. This is

not a serious problem, since we may include both Hermite-Gaussian func-

tions and contract them to yield the Cartesian Gaussian. Thus we may

use Cartesian Gaussian basis sets expressed as contracted Hermite-

Gaussians if necessary. The charge fitting and exchange fitting basis

sets are then chosen using the method described by Dunlap, Connolly, and

Sabin (1979a) for atomic and molecular calculations.

3-2. Overlap and Kinetic Energy Integrals

The overlap matrix elements S.. (k ) between normalized Hermite-

Ij -!-

Gaussian basis functions will be of the form:

S. (k,/) = (N.N.)-1/2 E exp[ik / R]

ij -/ i j R -- -

fd3r f(i.,a.,A.-R; r) f(n.,a ,A.; r) (3-1)

where N. and N. are normalization constants for the localized Hermite-

I J

Gaussians. The kinetic energy matrix elements T .(k. ) are evaluated at

the same time as the overlap matrix elements since they are of the same

general form as the overlap matrix elements such that:

T..(k ) = (N.N )-/2 E exp[ik / R]

d3

fd r f(n~ +26 ,a. ,A.-R; r) f(n.,a.,A.; r) (3-2)

s=i s i j J -

where the integer sets G refer to the sets (1,0,0), (0,1,0), and (0,0,1)

for values of s equal to 1, 2, and 3 respectively. For values of a and

b such that (a+b)R2 is greater than 1 (the value of I is rather arbi-

max

trary and may be adjusted in order to speed computational time as more

computational benchmarks are accumulated) these integrals will be eval-

uated using direct lattice sums of integrals of the form given in Equa-

tion A-5 in Appendix A. For values of (a+b)R2 (where R is the

max max

greater of R and R the primitive lattice vectors) less than or equal

1 -2

to 1, the overlap matrix elements and kinetic energy matrix elements are

evaluated using reciprocal lattice space techniques described in Appendix

C using expansions of the form given by Equation C-8.

L6wdin orthonormalization (L6wdin, 1956) is used to transform the

matrix elements of the secular matrix to an orthonormal basis so that the

S matrix in the secular equation will be the unit matrix. First the

unitary matrix U is found which diagonalizes the overlap matrix:

Ut S U = d (3-3)

The overlap matrix may then be transformed to the unit matrix by use of

the matrix T:

T S T = 1 (3-4)

where

T = U d-1/2 (3-5)

The secular equation of Equation 2-41 then transforms according to

Tt H T (T-1 C) = (T-1 C) c (3-6)

Approximate linear dependencies (which arise from the limits to precision

inherent in any computational procedure) in the basis may be handled by

excluding eigenvectors of the overlap matrix which have eigenvalues less

than some specified tolerance (for the purposes of this work the tolerance

was chosen to be 10-5).

3-3. Fitting Function Integrals

The integrals B.., C., Vijm. (k ), and X.jm (k/ ) currently are eval-

ij j ijm -// ijm -//

uated using direct lattice expansion techniques as discussed in Appendix

A. The matrix elements B.., defined in Equation 2-79, may be reduced to

the form:

n. 5/2 a.a.

B.. = E{(-1) J 2T R(n.+n., J A.-R-A.)

R a.aia.+a.)12 a.+a.

-- Ii I I J

n. n.

Z'2rZ a -[ R(n.,a ,A.+R-R ) + -- R(n.,a.,A.-R-R )]

Sa a. j -j -a a. i I -i ---a

Za Zb

SiJ ab l__b 3

using Equations A-40 and A-46. The C. integrals, defined in Equation

J

2-90, may be reduced to

2rZ Z Z

C. = E { Z a R(n.,a.,A.-R-R ) n. Z Z b (3-8)

J R a a a b -a -b

Equations A-41 and A-47 may be used to convert the matrix elements

Vijm (k ) to the form:

ijm -"//

V.. (N -1/2 E g

V. j (k /) = (N.N.)-1/2 Z exp[ik_ /-R] E g(K,n.,n.,a.,a.)

m J// R K I

a.a. n 5/2

f( ,A.-R-A.; 0) 1 {(-1) m 2-

a j R' (a.+a.)a (a.+a.+a )/2

(a.+a.)a a .(A.-R)+a.A.

R(n.+n.+n -K, -- A -R')

i j m a.+a.+a a.+a. -m -

I j m I j

27rn a.i (A .-R)+a.A.

-E Z -m R(n.+n.-K, a.+a., -- J-J R -R')}

a a.+a. i j i j a.+a. -a -

a i j i j

(3-9)

For the direct lattice expansions displayed in Equations 3-7 and 3-8 and

the inner lattice sum over R' in Equation 3-9, multiple expansion tech-

niques as described in Appendix B are used to increase the computational

speed of evaluating these integrals.

The exchange fitting integrals X. (k//), defined in Equation 2-95,

are also evaluated in a double lattice sum having the form:

-1/2

X.-. (k ) = (N.N.) E exp[ik_ R] E g(K,ni,n.,a ,a )

jm -// R j

a.a. n

(K, A.-R-A.; 0) E T( )3/2 m

Sa. +a -- -j R a.+a .+a

(a.+a.)a a.(A.-R)-a.A.

f(n.+n.+n -, i--- A -R'; 0)

I j m a.+a.+a a.+a. -m -

a j m I j

(3-10)

All of the above integrals could also be evaluated using reciprocal

lattice expansions as discussed in Appendix C. Since Vi. (k ) and

ijm -//

X. (k//) each contain two lattice sums as expressed in Equations 3-9

ljm -I/

and 3-10, we may choose to convert one or both of the direct lattice

summations to reciprocal lattice expansions. An approach using double

reciprocal lattice expansions has been tested for the case of all spher-

ically symmetric Gaussian functions for the Vijm. (k) and Xijm (k)

ijm -// ijm -//

integrals (which are the computationally most time consuming integrals

to evaluate). The results from these tests indicate that the large number

of sine and cosine function evaluations required make the double lattice

reciprocal lattice expansion algorithm at best equal in computational

time to the direct lattice expansion algorithm for integrals involving

diffuse Gaussian functions, which is the case requiring the most computer

time using the direct lattice expansion algorithm. For Gaussian functions

which were not extremely diffuse, the double reciprocal lattice expansion

algorithms required considerably more computational time than algorithms

using the double direct lattice expansions.

Another approach would be to convert only the inner lattice sums

over R' in Equations 3-9 and 3-10 to reciprocal lattice expansions. The

formulas necessary for incorporating this approach into computational

algorithms are derived in Appendix C. This approach has not yet been

incorporated into the computer code, leaving open one avenue for increasing

the computational speed of the computer code for calculating the thin

layer electronic structure.

3-4. Numerical Integration Techniques

In order to fit the exchange potential [or more precisely, the

function p (r)] to a linear combination of fitting functions m (r)

we must numerically integrate the expression:

Xm = f d r p (r) G (r) (3-11)

This leads to the approximate form of x :

m

x = E (r ) (r.) (3-12)

m I m -1

The integration points r. and weights w., which are chosen for numerical

integration over the parallelopiped primitive unit cell defined in

Section 2-1, are discussed in Appendix D.

The values of the primitive orbitals (r.; k ) and the values of

the periodic fitting functions G (r.) must be evaluated at each integra-

m --

tion point. For a primitive Bloch orbital p(r.; k /) constructed from

localized Hermite-Gaussian functions:

(r.; k//) = N-1/2 E exp[ik -R] f(n,a,A+R; r.) (3-13)

R

we may evaluate P(r ; k ) either by using the above direct lattice

expansion or by using the reciprocal lattice expansion:

(r.; k ) 1 an H [a /2 (z-A )]

-i -// aQ n3 z

nl+n2 n1 n2

E (i) (K+k ) (K +k//y)

exp[-(K+k/) /24a] exp[i(K+k)-((r.-A)] (3-14)

discussed in Appendix C. Evaluation of the periodic fitting function

G (r.), defined in Equation 2-46, is equivalent in form to that of a

Bloch function with a wave vector k,, equal to the zero vector. Thus

we may evaluate the fitting function with either a direct lattice expan-

sion or a reciprocal lattice expansion:

G(r.) = E f(n,a,A.r.) (3-15)

R

n3/2 1/2 n +n2

=- a H [a (z-A )] Z (-) 1 2

aQ n3 z K

n1 n2 2

K K exp[-K /4a] exp[iK-(r.-A)]

x y -1 -

The matrix elements G defined in Equation 2-50, may be evaluated

either analytically or numerically. We have chosen to evaluate G

mn

numerically using the same integration grid as used for the evaluation

of x This particular choice of the same numerical integration grid is

advantageous since the use of G in Equation 2-49 evaluated as

G = w. G (r.) G (r.) (3-16)

mn i m(ri n i

i

is equivalent to a least-squares fitting procedure over a weighted discrete

set of points. That is to say, the use of G as it is expressed in

mn

Equation 3-16 in the solution of Equation 2-49 yields a set of coeffi-

cients for the exchange fitting functions gm which minimize the quantity:

1/3 ~2

S w. {p/3 (r.) Z gm G (r.)} (3-17)

i m

We find that this particular means of evaluating G yields a method

mn

which requires fewer integration points (and thus less computational

time) for the same precision in the exchange fitting procedure than if

we evaluate G analytically. This is not difficult to understand when

one considers that the process which uses the numerically evaluated

values of G will solve the least-squares fitting procedure over a dis-

mn

create set of points with an error of the order of the precision of the

evaluated values of .(r.; k//) and G (r.), which are typically of the

J --I -/-

_5 8

order 10 to 10 The use of an analytically evaluated G matrix

mn

yields a process which will have an error of the order of the precision

of the numerical integration, which is typically of the order 10-1 to

10 As long as the number of grid points exceeds the number of fitting

functions Gm (r) such that the discrete least-squares fitting procedure

is not capable of yielding solutions with severe oscillations between

the discrete set of grid points, it appears that the use of a numerically

integrated value of G is preferable.

CHAPTER IV

RESULTS

4-1. Atomic Hydrogen Monolayer

The atomic hydrogen monolayer was chosen as the initial test case

for the approach outlined in Chapters 2 and 3 for calculating the elec-

tronic structure of thin films because atomic hydrogen possesses only

one electron and may be treated adequately using only a small number of

basis functions. Calculations have been performed on two different lat-

tice structures: a square lattice and a hexagonal lattice. Preliminary

calculations were performed (Mintmire and Sabin, 1980b) using a minimal

basis set consisting of three s-type Gaussian functions taken from van

Diujneveldt (1971) for the orbital basis. The basis set used for fitting

the charge density contains three s-type Gaussian functions with exponents

equal to twice the exponents used in the orbital basis. The basis set

used for fitting the exchange potential similarly contained three s-type

Gaussian functions with exponents equal to one-third the exponents used

in the charge fitting basis. This set will henceforth be referred to as

Basis 1. Table 4-1 presents the orbital exponents used in the orbital

basis set as well as the orbital exponents for the Gaussian functions used

to fit the charge density and the exchange potential.

The value of 2/3 was used for a in all calculations in this work,

except where explicitly noted. Other values of a could of course be used,

such as the avt or aHF values for the atomic systems given by Schwarz (1972).

Table 4-1. Basis sets for the atomic hydrogen monolayer.

Type of basis

Orbital

Charge fitting

Exchange fitting

Function

Type

s

s

s

x

s

2+ 2

x +y

2 2

x +y

Basis 1

0.151398

0.681440

4.50180

0.302796

1.362888

9.00360

0.100932

0.454296

3.00120

Exponents

Basis 2

0.151398

0.681440

4.50180

1.000

1.000

1.000

0.302796

1.362888

9.00360

2.000

0.100932

0.454296

3.00120

1.666667

aThe nomenclature for the function type is of the form xnyn2z n3 for

non-spherically symmetric Hermite-Gaussians. The superscripts nl, n2,

and n3 refer to the set of integers n = (nj, n2, n3) which define the

Hermite-Gaussian function and do not represent the Cartesian Gaussian

function one would normally associate with this nomenclature. Expres-

sions of the form "x2+y2" are to be interpreted as the sum of two

Hermite-Gaussian functions of type x2 and y2.

The formalism described in this work assumes that the value of a is con-

stant over all regions of space as opposed to the muffin-tin multiple-

scattering Xa method described by Johnson (1973), where a may have dif-

ferent values in different regions of space. Since one of the ultimate

goals of this method is to describe the electronic structure of systems

containing more than one type of atomic center, we feel that it is simpler

at this time to avoid the complicated question of what value of a is

appropriate for a particular system and use the value of a resulting from

the original derivation of the local density functional formalism approach

of Kohn and Sham (1965) and Gaspar (1954).

The values of the wave vector k were chosen such that all orbitals

were periodic over a translation equal to four times that of any primitive

lattice translation, which is equivalent to a Born-von Karman region

equal to a four by four array of contiguous primitive unit cells. Thus

this choice yields sixteen evenly spaced non-equivalent (by translation

symmetry) wave vectors in the central Brillouin zone. For the square

lattice, point group symmetry reduces the number of non-equivalent wave

vectors to seven. The Brillouin zone of the square lattice is depicted

in Figure 4-1, with the irreducible region of the Brillouin zone defined

by the triangle (rXM). Figure 4-2 illustrates the location of the seven

non-equivalent wave vectors used in the calculations. The hexagonal lat-

tice has only four non-equivalent wave vectors in the Brillouin zone.

The hexagonal Brillouin zone and the four non-equivalent wave vectors are

illustrated in Figures 4-3 and 4-4. The cohesive energies per particle

were calculated by taking the difference of the total energy at a given

lattice constant and the energy of atomic hydrogen calculated using the

molecular LCAO-Xa program (Mintmire, 1979). These results for the

M

Figure 4-1. Central Brillouin zone for the square lattice.

X M

Figure 4-2. Seven representative points in the irreducible triangle of

the central Brillouin zone for the square lattice.

M T' K

Figure 4-3. Central Brillouin zone for the hexagonal lattice.

M K

Figure 4-4. Four representative points of the irreducible triangle of

the central Brillouin zone for the hexagonal lattice.

cohesive energies per particle versus the lattice constant are depicted

in Figure 4-5. Initially these results were interpreted to imply that

the square lattice structure possessed a lower equilibrium cohesive

energy per particle (-0.035 Hartrees versus -0.005 Hartrees) and thus

a more stable equilibrium geometry relative to the hexagonal structure.

The effects of the choice of a were tested by first using the molec-

ular LCAO-Xa program to find the value of a which yields a total energy of

atomic hydrogen equal to -0.5 Hartrees using Basis 1. Further calculations

were performed using this value of a (equal to 0.781508) and it was found

that the cohesive energy curve underwent only minimal changes, with the

equilibrium cohesive energy per particle of the square structure changing

by less than 0.001 Hartrees with respect to the equilibrium cohesive energy

per particle of 0.035 Hartrees for the square lattice calculations with

a equal to 2/3.

More extensive calculations were performed using a (3s,1p) basis set

for the orbital basis set where p-type functions with exponents equal to

1.0 were added to Basis I as polarization functions. This basis set

(henceforth referred to as Basis 2) is presented in Table 4-1 along with

Basis 1. In addition to the slight improvement in the basis set, a finer

grid of wave vectors was used in the Brillouin zone, increasing the number

of non-equivalent (by translation symmetry) points in the central Bril-

louin to 256. This increases the number of non-equivalent points, when

point group symmetry is considered, in the Brillouin zone from 4 points

to 30 points for the hexagonal lattice and from 7 points to 45 points for

the square lattice. Figures 4-6 and 4-7 display these two sets of points

in their irreducible regions of the Brillouin zone for the square and

hexagonal lattices respectively.

1.5

1.0

0.5

- 0.5" ,

-s-a- - a -

1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2

ITEN: tlntosh Plat ,ten Me co 150D9180 01l. .48

Figure 4-5. Cohesive energy per particle versus lattice constant for

square (E[]) and hexagonal (A) lattices.

X M

000000

O000OO

(00000

000

00

S0

Figure 4-6. Forty-five representative points of the irreducible triangle

of the Brillouin zone for the square lattice.

M K

O

0O

Figure 4-7. Thirty representative points of the irreducible triangle of

the Brillouin zone for the hexagonal lattice.

The resulting cohesive energies per particle and virial ratios are

presented as a function of the lattice constant for the square and hex-

agonal lattices in Tables 4-2 and 4-3 respectively. A graphical display

of the cohesive energies per particle is made in Figure 4-8. The cohe-

sive energies are calculated as the total energy minus the energy of the

atomic monolayer at a lattice separation of 100 a.u. As can be easily

seen, these results differ markedly from the results of our preliminary

calculations. Our earlier inference that the square lattice is the ener-

getically preferred geometry at equilibrium is obviously refuted by these

more recent results. In fact by plotting the cohesive energy per particle

versus the square root of the area of the cross section in the periodic

plane of the parallelopiped unit cell defined in Section 2-1 (thus giving

the cohesive energy per particle as a function of the density of the

lattice sites) as displayed in Figure 4-9, we see that the results for

the two different geometries are almost indistinguishable. Most of the

difference between the two sets of results may be attributed to the in-

crease in the number of Brillouin zone points, since the improvement in

the basis from a (3s) basis to a (3s,lp) basis is only a modest improve-

ment compared to the increase of the total number of points in the Bril-

louin zone from 16 to 256.

It is to be noted that both Table 4-2 and Table 4-3 contain two sets

of results for the lattice constant of 5.0 a.u. Using the spin-polarized

version of the computer code at this particular lattice separation (and

only at this particular choice) leads to two stable (in the self-consistent

iterative procedure) solutions that correspond to a non-spin-polarized

(NSP) and completely spin-polarized (SP) solutions. Which result is ob-

tained depends on the initial choice of fitting coefficients to generate

Table 4-2. Calculated results for the square atomic hydrogen monolayer.

Lattice

Separation

(a.u.)

2.25

2.50

2.65

2.75

3.00

5.00a

5.00b

100.00

Total Energy

(Hartrees)

-0. 4883

-0.4921

-0.4918

-0.4909

-0. 4868

-0.4327

-0. 4521

-0. 4530

-V/T

1.867

1.995

2.062

2.103

2.188

2.187

1.934

1.972

Binding Energy

(Hartrees) (eV)

-0.0354

-0.0391

-0.0389

-0.0380

-0.0338

+0.0202

+0.0009

-0.963

-1.065

-1.057

-1.033

-0.920

+0.550

+0.024

anon-spin-polarized solution

spin-polarized solution

Table 4-3. Calculated results for the hexagonal atomic hydrogen monolayer.

Lattice

Separation

(a.u.)

2.25

2.50

2.65

2.75

3.00

5.00a

5.00b

100.00

Total Energy

(Hartrees)

-0.4823

-0.4917

-0.4924

-0.4922

-0.4892

-0.4378

-0.4503

-0.4530

-V/T

1.784

1.919

1.990

2.033

2.131

2.229

1.923

1.971

Binding Energy

(Hartrees) (eV)

-0.0293

-0.0382

-0.0394

-0.0392

-0.0363

+0.0151

+0.0027

-0.797

-1.040

-1.073

-1.066

-0.987

+0.412

+0.072

anon-spin-polarized solution

bspin-polarized solution

2.2 2.4 2.6 2.8

lciif;8? E-,,,_ ,=am

3.0 3.2 3.4

MEN, rtinsh Pic, slzte- realm II5/IWj 01l.Q.4,7

Figure 4-8. Cohesive energy per particle versus lattice constant for

square (]) and hexagonal (A) lattices.

DU I I

-0.70 F

-1.00

-I.10

0

2

-1.20 a

2.0

L' -: :. '- ,-_' .T L;r

1

-0.60 -

-0.70

s-0.9 -

-1.10 -,

-1.20 I I

1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2

IhEN: ttrntosh Plat Sltep mtiexlC ISO15,'S 01.A.47

Figure 4-9. Cohesive energy per particle versus modified lattice constant

for square (EI) and hexagonal (A) lattices.

the effective potential, since using fitting coefficients corresponding

to an SP system of overlapping atomic potentials yields the SP result,

while starting with coefficients from the NSP state at 3.0 a.u. yields

the NSP result. For the NSP results in general, it was necessary to

rerun the calculation in the NSP mode of the computer code in order to

confirm this assertion, since technical details involving choosing occu-

pation numbers at wave vectors near the Fermi surface prevented the SP

version from converging on a totally NSP result. The NSP results cal-

culated in both the NSP and SP modes of the computer code differed by

less than 10-5 Hartrees in the total energy. The NSP and SP results at

a lattice constant equal to 5.0 a.u. evidently indicate that the bands

corresponding to the different spins collapse at a lattice constant

slightly less than 5.0 a.u. This contention was tested by calculating

the cohesive energies per particle for both the NSP and SP modes at

lattice constants equal to 6.0 a.u., 8.0 a.u., and 12.0 a.u. Figure

4-10 illustrates the results for these calculations and demonstrates the

apparent crossing of the NSP state and the SP state at a lattice constant

between 4.0 a.u. and 5.0 a.u. A similar phenomenon has been reported

for the crystalline vanadium system (Hattox, Conklin, Slater, and Trickey,

1973) which exhibits the same collapse of bands of different spins to a

doubly degenerate band as the lattice constant is decreased from the

atomic limit.

64

2.0

1.5 -

NSP

/

0, /

-0.5 -

-1.5 I I

0 2 4 6 8 10 12 14

Figure 4-10, Spin-polarized and non-spin-polarized solution cohesive

energies per particle versus lattice constant for the square

lattice. Dashed portion of spin-polarized curve is

extrapolated from the solid portion of the line.

4-2. Atomic Beryllium Monolayer

Beryllium is an interesting system to study using the methods des-

cribed in this work for several reasons. Recent experimental measurements

of the work function of the (0001) surface of beryllium by Green and

Bauer (1978) indicate that beryllium has a work function of 5.1 eV, a

value much higher than previous measurements had indicated. Until the

results of Green and Bauer the recommended value of the work function

according to the standard reference (Fomenko, 1966) was only 3.92 eV.

Although Green and Bauer attribute this discrepancy in the work function

to possible oxidation of the beryllium surface in the earlier work, these

differences in experimental results pose a question worthy of theoretical

investigation. The interpretation of certain peaks in the Auger spectra

of beryllium (Suleman and Pattinson, 1971) also poses questions appropriate

for theoretical investigation using surface electronic structure computa-

tional methods.

In addition to the questions resulting from experimental study, beryl-

lium is an interesting prototype system for studying the surface of metals.

With only four electrons per atom, beryllium should require relatively

modest amounts of computational effort for an adequate treatment of the

surface. Also beryllium is only weakly bound as a dimer, with a dissoci-

ation energy which is experimentally estimated to be about 0.7 eV (Gaydon,

1968). The weak binding of the dimer is also indicated by theoretical

calculations such as the configuration interaction calculations by Bender

and Davidson (1967), which yield a strictly repulsive ground state poten-

tial energy curve, and the self-consistent electron pair/coupled electron

pair approximation (SCEP/CEPA) calculations by Dykstra, Schaefer, and

Meyer (1977) which leads them to estimate the calculated dissociation

energy to be about 0.03 eV at an equilibrium separation of about 8.5 bohr.

Hartree-Fock calculations by Bauschlicher, Liskow, Bender, and Schaefer

(1975) and by Jordan and Simons (1977) indicate the smallest cluster that

exhibits appreciable binding (0.46 eV per atom) contains four beryllium

atoms. This indicates that the cohesive energy of the solid (and pre-

sumably the monolayer) may be primarily due to the delocalization of the

electron charge over a collection of sites. The full effect of this de-

localization will be investigated by comparing cohesive energies of the

infinitely extended monolayer using the methods described in this work

with the binding energies calculated by Bauschlicher (1976) on beryllium

clusters and with experimental results for the real surface.

We have performed calculations on the atomic beryllium monolayer

using the (6s, 2p) basis given in Table 4-4. The 6 s-type functions are

taken from van Diujneveldt's (1971) 6s basis for beryllium. The 2 p-type

functions are the same as used by Bauschlicher (1976) in his cluster cal-

culations where the p-type functions were used with van Duijneveldt's 9s

basis for beryllium. The p-type functions are introduced as polarization

functions, since no p orbitals will be occupied in the ground state of the

beryllium atom. However since beryllium contains the greatest number of

electrons of the elements that have only s orbitals occupied, we may antic-

ipate that the p functions will have a greater contribution than equivalent

calculations for hydrogen, helium, or lithium. Thus questions of the bal-

ance between the s functions and p functions should be considered. In

his study of basis set effects, van Diujneveldt (1971) states that s/p

ratios greater than 5/2 (for 2 p-type functions) lead to unbalanced

behavior in Hartree-Fock calculations on the nitrogen dimer. Since

67

Table 4-4. Basis sets for the beryllium monolayer.

Orbital basis

Type Exponent

s 0.067376

s 0.198210

s 1.767558

s 6.819528

s 30.827597

s 204.906144

x 0.509

x 0.118

y 0.509

y 0.118

z 0.509

z 0.118

Charge fitting

basis

Type Exponent

s

S

s

s

s

s

s

s

x2+y2

x2+y2

z2

z2

0.134752

0.396420

3.535116

13.639056

61.655194

409.81229

1 .018

0.236

1.018

0.236

1.018

0.236

Exchange fitting

basis

Type Exponent

s

s

s

s

s

s

s

s

x2+y2

x2+y2

z2

z2

0.0449173

0.132140

1.178372

4.546352

20.551731

136.604096

0.33933

0.079667

0.33933

0.079667

0.33933

0.079667

beryllium has no occupied p orbitals in the atomic ground state (while

nitrogen does) we feel that our 6s/2p ratio is acceptable, especially in

the light of the results by Bauschlicher (1976) that a 9s/2p ratio yields

reasonable results. Preliminary calculations on the dimer were performed

with the above basis using the molecular LCAO-Xa program (Mintmire, 1979)

and a plot of the binding energy versus the nuclear separation is presented

in Figure 4-11. These results are consistent with the previously men-

tioned results of Dykstra, Schaefer, and Meyer and demonstrate the weak

binding of the beryllium dimer as well as provide a check on the balance

of the (6s, 2p) basis set.

Since beryllium possesses a hexagonal close packed structure in

the solid, the calculations were performed on a hexagonal lattice of beryl-

lium atoms. This structure is equivalent to a single layer from the (0001)

surface of the crystalline solid. The beryllium lattice possesses a Bril-

louin zone with the same symmetry as that of the hexagonal atomic hydrogen

monolayer previously displayed in Figure 4-3. We used the same array of

30 non-equivalent points in the Brillouin zone for the beryllium as for

the hexagonal atomic hydrogen monolayer calculations.

Table 4-5 presents some results of our calculations giving the co-

hesive energies and virial ratios of the beryllium monolayer as a function

of the lattice spacing. These results and Figure 4-12 imply that the

beryllium monolayer has an equilibrium lattice separation of 4.1 bohrs and

an equilibrium cohesive energy of -2.64 eV per particle. Since the virial

theorem is valid for the Xa method (Slater, 1972b), the virial ratio -V/T

should equal 2 at the minimum of the cohesive energy curve where the de-

rivative of the cohesive energy with respect to the lattice spacing equals

zero. We see in Table 4-5 that the location of the cohesive energy

Table 4-5. Calculated results for the atomic beryllium monolayer.

Lattice

Separation

(a.u.)

3.80

3.90

4.00

4.10

4.25

4.321

4.50

5.00

6.00

100.00

Total Energy

(Hartrees)

-14.293628

-14.296508

- 14.298019

-14.298351

-14.297454

-14.296885

-14.293275

-14.275816

-14.234844

-14.201443

-V/T

1.99370

1.99553

1.99880

2.00149

2.00586

2.00566

2.00837

2.01493

2.01263

1.99901

Binding

(Hartrees)

-0.092185

-0.095065

-0.096576

-0.096908

-0.096011

-0.095442

-0.091832

-0.074373

-0.033401

Energy

(eV)

-2.5084

-2.5867

-2.6278

-2.6369

-2.6125

-2.5970

-2.4987

-2.0237

-0.9088

4.0 5.0 E.0 7.0

N:. Mits Plot S~Ler rexi~ 2VC/O,13

02.0S.01

Figure 4-11. Binding energy versus internuclear separation for the

beryllium dimer.

MUMUM M-t'EtIEJI

0O

/

/

4.5 5.0 5.5

I.EN, ,tlntaoh PIot S.ptn "~'tm 20/I09,2 21.42.12

Figure 4-12. Cohesive energy per particle versus lattice constant for

the hexagonal beryllium monolayer.

-0.5 -

r-a

E8_

-2

-2.5 1-

*3.0 '

3.5

4.0

I I I I I

minimum (predicted by interpolation of the virial ratios to yield a pre-

dicted lattice constant of 4.05 a.u. and an estimated cohesive energy per

particle at that lattice constant of -2.63 eV) is in close agreement with

that predicted by the actual cohesive energy curve. These equilibrium

results differ slightly from the bulk equilibrium separation of 4.321

bohrs and equilibrium cohesive energy of -3.32 eV per particle (Kittel,

1976), but not more than could be attributed to the difference between

monolayer and bulk solid beryllium. Bauschlicher (1976) has performed

calculations using the Hartree-Fock method on various size beryllium clus-

ters with the same internuclear separation as the bulk solid. For the

various monolayer clusters, Bauschlicher reports binding energies of +0.22

eV per particle for a three atom cluster, +0.07 eV per particle for a six

atom cluster, -0.16 eV per particle for a seven atom cluster, and -0.63 eV

per particle for a fourteen atom cluster. The largest cluster Bauschlicher

considered was a twenty-two atom cluster in two layers which yielded a

A

computed binding energy of -0.89 eV per particle. These cluster results

apparently indicate a trend in binding energies as the cluster size in-

creases which should be consistent with our equilibrium cohesive energy

of -2.64 eV per particle, although one may infer that much larger clusters

are necessary for a proper treatment of the cohesive energy of the beryl-

lium monolayer. This is a point in favor of our earlier contention that

a primary cause of the binding energy of the beryllium monolayer is the

delocalization of the electronic charge over many nuclear sites.

Band energies were interpolated for 4096 points in the Brillouin zone

using a scheme adapted from the method proposed by Monkhorst and Pack

(1976) for three-dimensional Brillouin zone interpolations. This method

uses a discrete Fourier transform over the 256 evenly spaced points in the

complete Brillouin zone for which the band energies are computed. Band

energies at intermediate points were then evaluated using the resulting

Fourier expansion. Figure 4-13 displays the valence bands and lower empty

bands computed for the beryllium monolayer at a lattice spacing of 4.10 a.u.

The shape of the bands resemble qualitatively the bands reported by Loucks

and Cutler (1964) for crystalline beryllium, if we compare our results

with the solid energy bands for wave vectors lying in planes parallel to

the plane of periodicity of the monolayer. One may assume that the dis-

crepancy between our bands and those of Loucks and Cutler is due to the

differences between the monolayer and the crystal, and that the addition

of more layers would lead to multilayer bands more in agreement with those

of the crystalline beryllium solid. Further work is currently progressing

in the direction of calculations involving more than one layer of beryllium.

We may estimate the work function of beryllium from the negative of

the Fermi energy of our calculations (Kittel, 1976). We find that our

Fermi energy for the monolayer of -3.80 eV agrees quite well with the

recommended value of Fomenko (1966). However the discussion by Trickey

and Worth (1977) indicates that while the choice of a equal to 2/3 will

yield reasonable cohesive energies and band widths, one expects that this

value of a will lead to underestimating band gaps and the Fermi energy.

Thus the fact that we anticipate our Fermi energy to underestimate the

work function, and that increasing the number of layers will probably

change the Fermi energy in a manner not easily predictable, our estimated

work function of 3.80 eV is not necessarily in conflict with the reported

experimental work function of 5.10 eV (Green and Bauer, 1978).

The density of states for the beryllium monolayer was computed using

the interpolated set of band energies and is presented in Figure 4-14.

74

20 -

20

-20 -------------------------------------------

r E M T' K T F

IM|, rtinlant PIcL SELp ",exico 09/1OI 22.7.44

Figure 4-13. Band energies for the hexagonal beryllium monolayer. The

core band located at -100.98 eV is omitted.

-20 -15 -10 -5 a S 1 15 20

iMEN: rtintm~ Plot aptem i1e'xto 09/1 MG8 22.18.41

Figure 4-14. Total density of states for the hexagonal beryllium mono-

layer. The solid line shown in the inset is the total

density of states computed for crystalline beryllium by

Loucks and Cutler (1964).

MULLIM IONTLOM

The results are also consistent and in qualitative agreement with the

density of states computed for the crystalline beryllium solid (Loucks

and Cutler, 1964). The relatively flat region in the density of states

at energies around -10. eV in Figure 4-14 is understandable when one no-

tices that the bands in Figure 4-13 are nearly free electron like in char-

acter. One can demonstrate that the density of states for a two-dimen-

sional electron gas will be a step function. Thus the flat region of our

density of states is due to the nearly free electron like behavior of the

bottom of the lowest band in Figure 4-12.

Musket and Fortner (1971) have reported Auger Electron Spectra (AES)

results for beryllium indicating two peaks at 92 eV and 104 eV. Later

work by Suleman and Pattinson (1971) indicates that the AES peak at 92 eV

was evidently due to impurities on the beryllium surface, although they

agreed with the earlier conclusion of Musket and Fortner that the peak at

104 eV was due to a (2p, 2p) KVV transition at the solid surface. The

energy of this transition will equal the energy difference of an electron

in the 2p band (c 2p) and the energy of the hole it fills in the core K

shell (K ) minus the energy required to excite an electron from some dif-

ferent location in the 2p band (cp ). We may thus find an estimate for

the peak of the above mentioned transition by taking the difference

AE = 2p + e'p E K To find appropriate bounds for our estimate (since

2p 2p K

we do not anticipate monolayer results more accurate than 2 or 3 eV in the

light of our previous discussion of the effects of adding more layers to

the beryllium system) we may consider c2p and e2p jointly equal either to

the bottom of the second band in Figure 4-13, which corresponds to the

2pz band, or to the Fermi energy. Since we anticipate that the band energy

for the core band will not be very accurate, we will use the reported

results of -111 eV plus the work function (Siegbahn, et al., 1967).

Depending upon whether we use the experimental work function or our esti-

mated work function, we find a range of 105 to 109 eV for our estimate

of the (2p, 2p) KVV Auger peak. Thus we find one more source of agree-

ment of our monolayer results with those of the bulk solid, indicating

that even these simple monolayer calculations provide useful information

for describing the surface of the bulk solid.

These calculations were performed on an Amdahl V7 computer and re-

quired approximately 20 minutes of CPU time for each value of the lattice

constant. Preliminary calculations have been attempted for the beryllium

double layer with relaxed precision tolerances in order to reduce the CPU

time required. Unfortunately these preliminary calculations required

approximately 4 hours of CPU time per lattice constant for the double

layer and yielded unsatisfactory results due to precision problems involv-

ing the more diffuse basis functions. It is felt that the inclusion of

routines to incorporate the mixed lattice expansion algorithm for evalu-

ating the Vi. (k//) and X.. (k //) integrals described in Chapter III will

1Jk // I Jk -"//

both reduce the total required computational time and increase the pre-

cision of the evaluated integrals involving diffuse functions. Steps are

currently being taken to incorporate these techniques into the computer

code in order to perform calculations on the beryllium double layer and

other more complex two-dimensionally periodic systems.

CHAPTER V

SUMMARY

The primary purpose of this work has been to develop and test a new

approach for calculating the electronic structure of thin films, with the

ultimate objective of studying the electronic properties of solid sur-

faces. The calculations presented in this dissertation hopefully have

indicated some of the possible applications of our approach, although a

complete understanding of the strengths and weaknesses of the approach

outlined in this dissertation will require much more extensive testing.

The basic approach discussed herein has relied on the following assump-

tions:

1) Properties characteristic of the surface can be reliably

studied using a thin layer model;

2) the coupling of nuclear motion and electronic motion may be

neglected so that electronic properties can be calculated

reliably within the Born-Oppenheimer approximation; and

3) the local density functional (Xa) formalism yields results

which are useful in understanding the electronic properties

of surfaces.

The first assumption is intuitively reasonable given sufficiently

many layers in the film, and its use is supported by work on d-band

metals which indicates that the surface effects exhibited by films 15 to

to 20 layers thick accurately simulate the effects of the bulk surface

(Cooper, 1977). The Born-Oppenheimer approximation is a common one in

both molecular and solid state calculations, and barring evidence to the

contrary for a particular case, we may assume this approximation to be

reasonable. The third assumption has been shown empirically to yield

results in reasonable agreement with experiment in a multitude of cases

(for a review of applications of the Xa formalism see Slater, 1974, and

Connolly, 1976). The final test of all three assumptions, of course, is

that of comparison of results for a series of different systems using

the above model with the experimental results for those systems.

In addition to the previously mentioned assumptions used in estab-

lishing the mathematical model, various approximations have been made in

finding a computational solution for the model. These approximations,

which are capable of refinement to any limit of precision (or at least to

the limit of precision of the computational equipment) with a correspond-

ing increase in computational time, fall in four basic categories:

1) the representation of orbitals with a finite basis (the LCAO

approximation);

2) the general computational approximation of evaluating quantities

to within some predetermined tolerance;

3) the procedure for fitting the charge density and the cube root

of the charge density; and

4) the approximate treatment of the Brillouin zone introduced through

the use of periodic boundary conditions.

Of these approximations, the first two are well known and no detailed dis-

cussion of these points will be made. The effects of finite basis sets

are discussed by Schaefer (1972, pp. 56-83) while the effects of precision

tolerances are discussed by Clementi and Mehl (1971). No attempts have

been made to optimize basis sets, although bases optimized for atomic

hydrogen and beryllium in other works have been used in our reported

calculations. Appropriate tolerances were found so that the total ener-

gies were calculated to an estimated precision of 10-5 Hartrees. The

effects of using different numbers of points in the Brillouin zone has

been examined briefly for the hydrogen monolayer, with the conclusion

that an insufficient number of points can yield radically different co-

hesive energies from results using an adequate number of Brillouin zone

points. Examining the difference of band energies between adjacent points

in the Brillouin zone for the beryllium calculations using 256 points in

the Brillouin zone yields an upper bound to the error in the total energy

of about 0.05 Hartrees for this particular source of error. No testing

has been made of the errors introduced by the charge density fitting

procedures, although we anticipate that errors in the total energy due to

this source of error will be less than 0.001 Hartrees if the fitting

scheme for thin films behaves similarly to the fitting scheme for the

molecular case (Dunlap, Connolly, and Sabin, 1979a).

The physically reasonable results produced for the hydrogen mono-

layer indicate that there are no inherent flaws in the mathematical for-

malism. The beryllium monolayer results further indicate the feasibility

of our approach as well as demonstrating the applicability to experimen-

tally interesting systems. The comparison of our theoretically estimated

values for the equilibrium lattice spacing, equilibrium cohesive energy

per particle, and work function with the experimentally determined quan-

tities for the beryllium system and with other calculated results indicate

that even such a simple model as the monolayer yields useful results.

The primary difficulty with our approach at this time is the amount of

computer time required to produce adequately precise results for systems

larger than the atomic beryllium monolayer. Current work to incorporate

the reciprocal lattice expansion techniques described in Appendix C as

well as structural changes in the computer code to optimize the order of

operation should lead to a substantial reduction in required computer

time.

In summary we may state that this dissertation hopefully represents

the beginning of an extended study of the electronic properties of sur-

faces. Needless to say, the methods presented herein are not a finished

product, but rather the first stage of development of what we believe to

be a line of study that can be extremely fruitful in terms of understanding

the basic properties of solid surfaces.

APPENDIX A

HERMITE-GAUSSIAN FUNCTIONS AND INTEGRALS

A-l. Definition of Hermite-Gaussian Functions

The basis functions used in this work are functions which are com-

posed of specific combinations of Cartesian Gaussians (which are of the

I m n 2-a2

form x y z exp(-ar ), where 1, m, and n are non-negative integers).

These functions may be referred to as Hermite-Gaussian functions, which

is the nomenclature introduced by Zivkovic and Maksic (1968) in their

discussion of the use of these functions for molecular calculations. A

Hermite-Gaussian function f(^, a, A; r) is defined:

f(n, a, A; r) = n exp(-air Al2) (A-1)

where n refers to a set of three non-negative integers (nl, n2, n3) which

define the operation an:

n = a an2 (A-2)

A- aAnl 3An2 aAn3

x y z

This notation is consistent with that of Zivkovic and Maksic, although

other workers (Golebiewski and tMrozek, 1973; McMurchie and Davidson,

1978) have introduced slightly different notations. The term "Hermite-

Gaussian functions" derives from the fact that given the form of the

Hermite polynomials:

H (x) = (-1)n exp(x2) 3n exp(-x2)

n x

(A-3)

the Hermite-Gaussian function f(^, a, A; r) may be expressed:

f(n, a, A; r) = an/2 H n[a /2(x-Ax)] Hn2[a /2(y-Ay)]

H n3[al/2(z-Az)] exp(-air A 2)

n 3 -

(A-4)

where n denotes the sum of the components of n (i.e., n = nI + n2 + n3)

These functions are useful in computation since they enable us to per-

form the calculation of integrals and other numerical data, as pointed

out by Zivkovic and Maksic (1968), in three steps:

1) for any integration or summation which does not involve A and

is linear in f(6, a, A; r), bring the differentiation operator

outside the integration or summation;

2) calculate the appropriate integral or sum for spherically sym-

metric Gaussians in closed form; and

3) apply the differentiation operator to the analytic form of the

resultant expression.

As an example we may compute the overlap integral between two

Hermite-Gaussian functions centered on two different centers A and B.

d3r f(A,a,A; r) f(n',b,B; r)

= n 3n fd3r exp[-alr-A 2 b|r-B 21

= n ;n' [r/(a+b)]3/2 exp[- ab A BI2/(a+b)]

SA BB /(

= (-1)n' [7r/(a+b)]3/2 f(6+6', ab/(a+b), A-B; 0)

(A-5)

84

A-2. Products of Two Hermite-Gaussian Functions

The overlap integrals of Equation A-5 are a simple example of inte-

grals which involve the product of two Hermite-Gaussian functions. Such

integrals may be evaluated with less complexity by decomposing the prod-

uct of two Hermite-Gaussian functions into two factors. Consider the

product of two Hermite-Gaussian functions:

f(6, a, A; r) f(6', b, B; r)

aB= n exp[-a

= B

=A AB exp[

Defining the quantitities P and Q:

P = A B

Q = (aA + bB) / (a + b)

Ir A 2] exp[-blr- B 2]

(A-6)

ab 2 aA + bB]2

ab A-B 2] exp[-(a+b) r a -+ b- ]

+b (a+b)

(A-7)

(A-8)

we may express the product of two Hermite-Gaussian functions in terms

of functions g(P) and h(r Q):

f(n, a, A; r) f(n', b, B; r) = 3A 3n g(P) h(r Q) (A-:

where

g(P) = exp -ab IPl2]

(A-10)

h(r Q) = exp[ -(a+b)ir Q 2]

(A- ll)

Repeated taking of derivatives yields:

g(P) h(r-Q) =

nl n2

E Ek

kj=0 k2=0

n3 nI' n2' n3'

E E E E

k3=0 kI'=0 k2'=0 k3'=0

k3 A B )

x (nk) (n

k, k2

an-k an'-k'

A B-

h(r-Q) }

(A-12)

We may rearrange the sums over the indices k. and k.', using the new

I I

indices K. and K.', where

I I

K. = k. + k.'

I I I

K.' = k.'

I I

This rearrangement yields the expression:

(A-13)

(A-14)

g(P) h(r-Q) =

nl+nl' n2+n2'

E E

KI=O K2=0

min[nl ,Kl]

x E

K]'=max[0,Kl-nl]

min[n2' ,K2]

K2 '=max[0,K2-n2]

min[n3' ,K3]

K3 '=max [0,K3-n3]

A A

K' g(P n-K+K'

B A

Using the relationships

ag = ag = _ag

aP. aA. aB.

I I I

^ ^n

an an

A B

A a

n 3+n 3'

E

K3=0

n'-K'

B

h(r-Q) }

(A-15)

(A-16)

n3 nl ') n2 )

k3 kj' k2'

( n3 > ni n2: (n31

K3-K3') K1 Kg2 K3'

x 1 n1 ) 2 ,)

(I-KI K2-K2.

86

and

ah a + b ah a + b ah

3Q. a aA. b aB. (A-7)

I -I

we may rewrite Equation A-15:

aA aB g(P) h(r-Q) = E g(K, n, n', a, b) {p g(P)}

K

x n+n'-K h(r-Q)} (A-18)

Q

The function g(K, n, n', a, b) is defined:

^ ^ ^ 3 n-K+K' n'-K'

g(K, n n', a, b) = H { E a nn'-K Kn K.) n( .) (A-19)

i=1 K.' (a+b)n+n'-K

where the indices K.' have as their respective ranges:

max[O, K. n.] < K.' < min[n.' K.] (A-20)

Furthermore we may use the relationships:

aK g(P) = f(K, ab/(a+b), P; 0) (A-21)

and

A A A

n+n'-K

aQ h(r-Q) = f(n+n'-K, a+b, Q; r) (A-22)

to yield the expression

f(n, a, A; r) f(n', b, B; r)

= E g(K, n, n', a, b) f(K, ab/(a+b), P; 0)

K

f(n+n'-K, a+b, Q; r)

(A-23)

Thus we see that Equation A-23 may be used to reduce integrals involving

the product of two Hermite-Gaussian functions of the same coordinate r

to a sum of integrals involving a single Hermite-Gaussian function of the

coordinate r.

A-3. Nuclear Attraction Integrals

The electrostatic interaction integral of a charge distribution

described by a Hermite-Gaussian function f(n, a, A; r) with a nuclear

center of charge ZC located at the position vector C may be denoted:

I3 1

[f(n,a,A;r) Z] = Z fd3r f(n,a,A;r) 1 (A-24)

C C Ir CLI

We may use Boys' (1950) relation for spherically symmetric Gaussian func-

tions:

[f(0,a,A;f) I ZC] = X Fo(aT2) (A-25)

to evaluate this integral, where A is the constant:

A = 27f ZC / a (A-26)

and F0(s) is the incomplete gamma function:

F0(s) = J du exp[-su 2] (A-27)

0

The vector T is defined as

T = A C (A-28)

Thus the electrostatic interaction integral in Equation A-34 may be

expressed:

[f(n, a, A; r) I Z = aA Fo(aT2) (A-29)

nA

McMurchie and Davidson (1978) have introduced the function R(n, a, r):

^ 2

R(n, a, r) = an FO(ar) (A-30)

where we have employed slightly different notation from that of McMurchie

and Davidson. This function may be expressed in the expansion:

R(n, a, r) = ^ an-k 2n-2k (-1)k xn1-2k n2-2k2 zn3~2k3 Fnk(ar2

k

n n2! n3! (A-31)

^ -- _I _2 n^ (A-31)

k1! (n1-2kl)! k2! (n2-2k2)! k3! (n3-2k3)!

where the incomplete gamma functions F (s) may be defined:

F (s) = 1 du u2n exp[-su2I (A-32)

n 0

and the interval over which the indices k. range is:

0 5 2k. < n. (A-33)

I I

The function R(n, a, r) may be evaluated using recursion relationships

if we define the quantity T

n,1 ,m,j

Tnl m = (-j)n+1+m 1 du {H [ux/a] H [uy/a] H [uz/a]

n,1,m j 0 n 1 m

u n+]+m+2j e 2j} (A-32 D

(A-34)

This quantity satisfies the following relationships:

T F.(ar2)

and

T ,n2 3,0 = R(n, a, r)

nl,n2,n3,0

(A-35)

(A-36)

If we use the method described by McMurchie and Davidson (1978) to eval-

uate the incomplete gamma functions F (s), the functions R(n, a, r) may

be evaluated using the recursion relationships:

T ,n3+1 = -2a {zT ,n2,n,j+l + n3Tn1 n2,n

nl ,n2,n3+l ,j n1 ,n2,n3,j+l n1 ,n2,n3-1 ,j+l

nnl,n2+l,n3,j

Tnl+1,n2,n3,j

and

S-2a {yTni,n2,n3,j+1 + n2Tn2-1 ,n3,j+l }

= -2a {xTn ,n2,n3,j+1 + nT -1,n2,n3,j+l }

(A-37)

(A-38)

(A-39)

Thus we see that the electrostatic interaction integral of a Hermite-

Gaussian function with a nuclear center may be expressed:

[f(n,a,A;r) I ZC] = (27r/a) ZC R(n, a, A C) (A-40)

Using the relationship expressed in Equation A-23 we may also express

the electrostatic interaction integral of the product of two Hermite-

Gaussian functions with a nuclear center:

2 A A A A

[f(n,a,A;r) f(n',b,B;r)|Z ] = Z E g(K,n,n',a,b)

-- (a + b) -- ^

K

f(K,ab/(a+b),A-B;O) R(n+n'-K,a+b,Q-C) (A-41)

where Q is the vector defined in Equation A-8.

A-4. Electron Repulsion Integrals

The electrostatic interaction integral of two Hermite-Gaussian func-

tions with each other may be expressed:

[f(n,a,A;r) I f(n",c,C;r')] = fd r fd3r' -

\L~- I

f(n,a,A;r) f(n",c,C;r') (A-42)

We may again use the relationships derived by Boys (1950) for spherically

symmetric Gaussians:

[f(O,a,A;r) I f(O,c,C;r')] = \ F0( acT2/[a+c]) (A-43)

where X is now the new constant

A = 27n5/2 / {ac (a + c)1/2} (A-44)

and the vector T is defined:

T = A C (A-45)

The methods used in Section A-3 for evaluating the nuclear interaction

integrals may be extended trivially to yield:

A I )n" 5/2 1/2

[f(n, a, A; r) f(n", c, C; r')] = (-1) {2 /2/[ac (a+c)1/2

R(n+n", ac/(a+c), A C) (A-46)

Similarly Equation A-23 may be used to evaluate the electrostatic inter-

action integral of two Hermite-Gaussian functions with a single Hermite-

Gaussian function:

[f(n, a, A; r)f(n', b, B; r) I f(n", c, C; r')]

= (-1)n" {25/2/[(a+b)c (a+b+c)1/2]} E g(K, n, n, a, b)

K

f(K, ab/(a+b), A-B; 0) R(n+n'+n"-K, (a+b)c Q-C) (A-47)

a+b+c

This procedure could be extended straightforwardly to the case of two

Hermite-Gaussian functions interacting with the product of two other

Hermite-Gaussian functions. This has not been done since no integrals

of this form are necessary for the current work.

APPENDIX B

MULTIPLE EXPANSIONS OF COULOMB INTEGRALS

B-1. Multipole Expansions of Electrostatic Interactions

Consider the sum of electrostatic interaction terms of a charge den-

sity p (r) centered about the origin and a charge density p2(r) centered

about a lattice site R of a two-dimensionally periodic Bravais lattice:

U = E fd3r fd3r2 p1(r ) p2(r-2) 1 (B-1)

R L2 -

where the coordinates'r and r2 are illustrated in Figure B-1. These

-1 --2

charge densities are assumed to be localized about their respective cen-

ters, such that p(r) goes to zero faster than any polynomial of |jL in

the asymptotic limit as Ifr tends toward infinity. This condition may

be more restrictive than necessary, but it is easily satisfied by the

Gaussian functions used in this work as well as by exponential functions.

If this expansion is summed in order of increasing |RI, the convergence

of the Coulomb terms is usually slow so that many terms must be included

in a truncated expansion so that the sum may be computed within a spec-

ified precision.

Piela and Delhalle (1978) have found that the use of a multiple

expansion as an approximation for all terms having |R[ greater than some

radius R' for one-dimensionally periodic systems such as polymers is an

effective method for speeding the computational procedure. The multi-

pole expansion considered in this work uses the bipolar expansion of the