Title: Non-muffin-tin correction to the total molecular energy in the multiple-scattering Xa method of molecular calculation
Full Citation
Permanent Link: http://ufdc.ufl.edu/UF00098370/00001
 Material Information
Title: Non-muffin-tin correction to the total molecular energy in the multiple-scattering Xa method of molecular calculation
Physical Description: xi, 240 leaves. : illus. ; 28 cm.
Language: English
Creator: Danese, John Bryan, 1944-
Publication Date: 1973
Copyright Date: 1973
Subject: Quantum theory   ( lcsh )
Molecular dynamics   ( lcsh )
Physics thesis Ph. D
Dissertations, Academic -- Physics -- UF
Genre: bibliography   ( marcgt )
non-fiction   ( marcgt )
Bibliography: Bibliography: leaves 235-239.
General Note: Typescript.
General Note: Thesis -- University of FLorida.
General Note: Vita.
 Record Information
Bibliographic ID: UF00098370
Volume ID: VID00001
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: alephbibnum - 000582456
oclc - 14114258
notis - ADB0831


This item has the following downloads:

nonmuffintincorr00dane ( PDF )

Full Text



John Bryan Danese




1III111IIIII1II II1I1I1UII1III11111111111 1111111111111111111111


I would like to acknowledge Professor J. W. D.

Connolly who suggested the topic for this study and would

like to thank him for his many suggestions, criticisms,

and continuous encouragement during this research.

I would like to acknowledge the members of my

advisory committee, particularly Dr. J. R. Sabin, for their

interest and criticisms. I would like to extend special

thanks to Dr. S. B. Trickey for his discussions of this

work and communication of pertinent information. Thanks

are due to my other colleagues of the Quantum Theory Project

at the University of Florida for their piquant but stimulating

discussions and criticisms. I wish to extend special

appreciation to Dr. N. H. F. Beebe for general programming

assistance but most especially for providing me with the

line plotter program which proved so useful in the analysis

of my results (this program is originally due to Dr. H.

McIntosh of the University of Mexico). I would also like

to express my appreciation to the members of the North

East Regional Data Center for their cooperation and

assistance. Finally, I would like to thank my wife for

her patience and encouragement.

The financial support from the National Science

Foundation for part of this research is gratefully





. . ii

*.. V1

. ... vii

. .i

.. .. .1

.. .. 1
.. .. 3

.. .. 6

. . .

. .. .9
. .. .0

. . .15


LIST OF TABLES* ********

LIST OF FIGURES ................

ABSTRACT* ........


INTRODUCTION. ..........

1.1 Molecular Methods ..........
1.2 The Muffin-Tin (MT) Approximation ..
1.3 The Non-Muffin-Tin (NMT) Correction
to the Muffin-Tin (MT) Xce Total

Energy. . . . . . . .



2.1 Introduction. ............
2.2 Xa Theory ..............
2.3 Derivation of the Xa One-Electron
Equations .. .... .. ... .


EQUATIONS ............

. .. ..26

3.1 Introduction. .............
3.2 The MT Approximation. .........
3.3 The Separation of the Xcl Total Energy
into MT and NMT Parts ........
3.4 The Derivation of the MT Xcr One-
Electron Equations. .........
3.5 The MT Xa One-Electron Orbitals and
the Multiple-Scattering (MS)
Formal ism . . . . . . .
3.6 Explicit Expressions for the Components
of the NMT Correction a.....

. . .26
. . .26

. . .31

. . .38

. . .46


. .


ENERGY .... ....

. . 8

BIBLIOGRAPHY. ...............



4.1 Introduction. ...............
4.2 The Explicit Relationship Between
and . . . . . . .
4.3 Firs aOrder Co~rections to the Eigenvalues.
4.4 The SCF Cycle ...............

SAMPLING. . . . . .. . .














. 109

. 114

. 116

. 116
. 122

. 124

. 151

. 179

. 185

. 192

. 199



Introduction. ..............
Importance Sampling in One-Dimension. ..
Importance Sampling in n-Dimensions ...
Some Specific Sampling Functions. ....
A Note on the Haselgrove and Conroy
Methods of Numerical Integration. ...


6.1 Introduction. ..............
6.2 The C2 and Ne2 Molecules. ....
6.3 Features of the C2 Molecule and
A at an Internuclear Separation
of 2.48 Bohr. . .. .. .. .. .
6.4 The Resul ~s for C2 and Ne2 aa'~nd Prncipal
Conclusions ..............
6.5 Further Discussion and Future Work. ..




aV IN,IN and AVoIN c.....

INTEGRALS .. .. .. .. . .. .



. . .

. .

.. 216


.. 229

.. 235

. .240


Table Page

6.1 The NMT Correction for C2 and Its Regional
Parts as a Function of Internuclear Separation 153

6.2 The NMT Binding Energy for C2asaunto
of Internuclear Separation 154

6.3 Spectroscopic Constants for C2 159

6.4 The NMT Correction for Ne2 and Its Regional
Parts at the Nuclear Separation 6.0 (Bohr) 161

6.5 The NMT Binding Energy of Ne at the
Internuclear Separation 6.0 Bohr) 162

6.6 Comparisons of the NMT Binding Energy at the
Internuclear Separation 6.0 (Bohr) with
Other Sources 164


Figure Page

3.1 Schematic Representation of the C2 Mlcl
in the MT Approximation. 28

3.2 Vector Relationships Used for the Derivation
of (a) Va,8and (b) Vora. 54
Ne Ne
4.1 Flow Diagram for Conceptual Exact Xa Calcu-
lation and a Conceptual MT Xa Calculation. 84

4.2 Flow Diagram for an Actual MT Xa Calculation
and the Major Steps in the NMT Calculation. 87

5.1 Contour Plot of AlI(r) in the X-Z Quadrant
of the C2 Molecule. 102

5.2 Schematic Plot (a) and Contour Plot (b) of
AI(r) 110

6.1 Line Plot of the ag Valence Orbital of the C2
Molecule. 92127

6.2 Line Plot of the oz Valence Orbital of the C2
Molecule. 128

6.3 Line Plot of the Total Electron Density of
the C2 Molecule. 130

6.4 Line Plot of Ap(r) for the C2 Molecule. 132
6.5 Contour Plot of Ap(r) in the X-Z Quadrant of
the C2 Molecule. 133

6.6 Contour Plot of AICL(r) in the X-Z Quadrant
of the C2 Molecule. 138

6.7 Line Plot of AICL(;) for the C2 Molecule. 139
6.8 Schematic Division of the X-Z Quadrant of the
C2 Molecule into Sampling Regions. 140
6.9 Contour Plot of alX(r) in the X-Z Quadrant
of the C2 Molecule. 143




6.10 Line Plot of AI ( )' for the C2 Molecule. 144

6.11 Contour Plot of AI (r) in the X-Z Quadrant
of the C2 Molecule. 147

6.12 Line Plot of 6I(r) for the C2 Molecule. 148

6.13 MT, NMT, and CI Potential Curves for the
1C + Ground State of C 157

D.1 Atomic-Inner Vector Relationships Used
in the Derivation of nV a,IN and aV IN,IN, 201

D.2 Atomic-Atomic Vector Relationships Used
in the Derivation of AV a,IN and AV IN,IN. 204
c c
D.3 Atomic-Outer Vector Relationships Used in
the Derivation of aV oIN, 214

G.1 Volume Integral Relationships Used in the
Derivations of the Bond Sampling Function. 230


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



John Bryan Danese

August, 1973

Chairman: J.C. Slater
Co-chairman: John R. Sabin
Major Department: Physics

An investigation of the muffin-tin (]MT) approximation

to the total electron density in the multiple-scattering

Xac method of molecular calculations is carried out. A

method is developed to improve upon the shortcomings of the

MT approximation in the calculation of the Xa total energy

by including the non-muffin-tin (NMT) correction to the

density in the Xar total energy calculation. This results

in a small NMT correction term to the total MT Xa energy

which takes the form of numerical three- and six-dimensional

integrals. It is shown that the method of including the NMT

correction terms is analogous to the first order perturbation

calculation and it is inferred from this and the smallness

of the correction that the resulting total energy is quite

close to the exact Xal energy. An extensive numerical samp-

ling and integration scheme is developed with which to inte-

grate the NMT correction term whose integrands are rather

complex functions.

Calculations are performed on two molecules--C2 and

N2-which it is argued should provide definitive tests for
the method of improving the MT XcL total energy with the
NMT correction. The multiple-scattering Xcr method in the

MT approximation predicts that the potential curves for
both molecules in their ground states are repulsive. In

reality both molecules in their ground state bind with C2
binding with 0.468 Rydbergs and Ne2 with 0.00023 Rydbergs.

(The extremely small binding of Ne2 is attributed to van der
Waals attraction.)

It is found that in the case of both molecules with

the NMT correction included a significant qualitative and

quantatative improvement is made in the total energy. For

C2 three points on the potential curve near the equilibrium
internuclear separation are calculated with the NMT correc-

tion and a thorough analysis of the correction and its com-

ponent parts is made. For C2 a binding potential curve is
predicted which compares quite favorably with a CI calcula-
tion on C2 and which predicts spectroscopic constants in

substantial agreement with experiments. For Ne2 only one
point on the potential curve is calculated very near the
experimental equilibrium internuclear separation. With the

NMT correction binding is definitely predicted for Ne2 but
due to the difficulty of this particular calculation the

error in the binding is such that the precise amount of
binding cannot be specified.

In light of these results for these two dissimilar

and difficult test case molecules we argue finally that not

only does the NMT correction give results very close to

the exact Xcc results but also that these results provide

additional evidence for the accuracy of the Xar approxima-

tion itself .



1.1 Molecular Methods

The understanding of the chemical properties of

materials from first principles has been the quest of quantum

theory since its inception in the 1920's. In this quest

quantum theory has provided the crucial frame of understanding

and the calculational tools for many other fields of endeavor

in physics, chemistry, and more recently biology. Within the

field of quantum theory the main direction taken has been

toward the calculation, prediction, and understanding of the

properties of molecules (and of solids) of ever--increasing

size, weight, and complexity.

There are several categories or methods of molecular

calculation which can be broadly classified as ab initio,

semi-empirical, and multiple-scattering Xa(MSXa). This

dissertation is concerned with a particular aspect of the

MSXa method. Specifically, the purpose of this dissertation

is to investigate the effect of the so called "muffin-tin"

(MT) approximation to the density in the total energy

calculation in the MISXaL method and to provide remedies for

the shortcomings of this approximation.

The MSXu method is a relatively new method first

suggested by Slater (1965) and its implementation begun by

Johnson (1966J. The ab initio methods such as the LCAO-MO-SCF

(linear combination of atomic orbitals-molecular orbital-self-

consistent field) Hartree-Fock method are capable of usefully

accurate calculations on smaller molecules such. as water

(H20) or ammonia (NH3) (see Slater (1963a)) but due to the
computational times involved cannot be realistically applied

to larger molecules such as benzene (C6H6). The semi-

empirical methods (see Pople and Beveridge (1970)) have

been applied to systems such as benzene and even somewhat

larger systems but their results are usually semi-quantitative.

Further, they are parameter dependent (sometimes involving

many parameters) and the selection of these parameters is

something more of an art than a science. The MSXa method

is a promising method since it is capable of an accuracy

in ionization, transition, and total energies quite comparable

to that of the ab initio LCAO-MYO-SCF Hartree-Fock method

but is about two orders of magnitude faster (see Slater and

Johnson (1972)). This of course enhances the usefulness

of the method making possible accurate calculations on larger

systems such as sulfur hexaflouride (SF6) (see Connolly and
Johnson (1971)).

There are three distinguishing features of the MSXa

method. It makes use of the Xa local density exchange

approximation, the "muffin-tin" (MT) approximation to the

total electron density and potential, and the multiple-

scattering expansion of the molecular orbitals. The APW

(augmented-pla ne-wave) method and the KKR (Kohn-Korringa-

Rostoker) methodsof solid state (see Slater (1963b)) employ

an exactly analogous muffin-tin approximation. Hence similar

conclusions may be drawn for these methods from the study

of the muffin-tin approximation in the MSX~method. The

MSXa method is a single particle molecular orbital self-

consistent field method. The details of the method as

well as the relevant references will be given in Chapter II

of this dissertation. Two particularly comprehensive

recent reviews on the subjects respectively of the Xa method

and the MSXa theory are those of Slater (1972a) and Johnson

(1973). The Xa approximation and the multiple-scattering

formalism have been the subject of much study and even

though reasonably well understood are still under development.

However, the state of knowledge of the effects of the muffin-

tin approximation in the MSXa theory and the developemnt

of schemes for overcoming this approximation are neither

as advanced nor as well developed as the former two.

1.2 The Muffin-Tin (MT) Approximation

In the MT approximation the coordinate space of

a molecule is schematically divided into contiguous spherical

regions surrounding the nuclei and a large outer sphere

surrounding the whole molecule. Within the atomic spheres

and outside the outer sphere the one-electron Xa( potential

and the total electron density are spherically averaged.

Between the atomic spheres and within the outer sphere in

the so-called intersphere or interatomic region the potential

and density are volume averaged to a constant.

The motivation behind the MT approximation of the

density and potential is simple. Using the MT form of

the total electron density p in the Xa total energy formula--a

functional of p which we will designate --produces the
MT Xa total energy formula--a functional of the MT density

p which we will call . All of the integrals in may be reduced to one-dimensional numerical integrals

whereas those of involving p are three and six
dimensional numerical integrals. From the MT Xa energy

functional we variationally derive (see Chapter III)
the MT Xar one-electron equations. The MT potential in this

equation is related variationally to the MT approximation

to the density which defines the functional . Hence
the MT approximation of the Xa potential follows from that

of the density and there is really only one MT approximation--

that of the density. The self-consistent solutions of the

MT Xa one-electron equations are the set of orbitals then

which minimize the MT Xax energy functional. In making the

MT approximation to the Xa total energy we leave out the

effect of the non-muffin-tin (NMT) components of the density

upon the total Xa energy. We call the difference between

the Xa total energy and MT total energy the NMT correction

to the MT total energy.

Although it has been found that the ionization and

transition energies calculated by the MSERi method do not

seem to be particularly sensitive to the MT approximation

(see Connolly et al. (1971) ) it has not been found that the

MT Xa total energy calculations have been uniformly good.

For example, fairly good symnmetric mode potential curves

for the ground states of Li2 (see Johnson et al.(1972)) and

CHq (see Danese (1972)) have been calculated in the MT

approximation. Two solid state APW calculations employing

the MT approximation by Averill (1971b)and Hattox (1972) on

the compressibilities as well as other volume dependent

properties of cesium and vanadium produced semi-quantitative

results. However a conformational study of the water

molecule by Connolly and Sabin (1972) using the MT approxima-

tion predicted the result that the water molecule was linear.

These and other evidence have led to the conclusion that the

MT approximation to the total energy is not always particularly


This dissertation has two purposes. One purpose is

to investigate the effect of the MT approximation to the

density upon the Xa total energy through an analysis of

the non-muffin-tin corrections to the MT Xa total energy.

The second purpose is to develop a general and practical

method of evaluating the NMT correction to the MT Xac

energy and thereby to improve upon the MT Xac energy calcula-


1.3 The Non-Muffin-Tin (NMT) Correction to the
M~uffin-Tin (NiT) Xc: Total Energy

In Chapter III w!e howJ hpw for a given set of

assumed orbitals the Xa total energy may be separated

into two parts--the MdT Xa total energy part and a

second part which is the NMt~T part of and which we

designate as h In Chapter IV we examine the relation-

ship between and via A. We show that

using the set of orbitals which minimizes to calculate

is analogous to a first order perturbation calcula-
tion--6 being the perturbation to Hence the

energy calculated in this manner has only a second order

dependence upon the difference between the exact Xa orbitals

and the MT Xa orbitals. Because of this and because the

NMT correction (estimated) in most cases appears to be

quite small we argue that our approximation to the exact

Xa energy should be fairly good. However, even though the

correction may be small it can be quite important if one

is interested in a small effect such as the binding of a

molecule. Hence to be able to calculate binding energies

accurately with this correction in cases in which the

binding predicted by the MT approximation is either very

bad or non-existent (unbound) would be a good test of

this method.

The NMT correction A takes the form of three-

and six-dimensional integrals which must be evaluated numeri-

cally. In Chapter V we develop an importance sampling

scheme which combined with the M~onte-Carlo method of

numerical integration results in an admittedly complex but

nevertheless workable integration scheme whereby the numerical

integrals of n may be performed in very reasonable
calculational times.

In Chapter VI we give the results and conclusions.

We chose as our test cases two relatively dissimilar but

simple molecules--C2 and Ne2. However, they are difficult
cases. Both are actually bound in their ground state but

neither is bound in the MT Xa approximation. In fact

for the ground state of C2 the MT Xa approximation predicts

a distinctly repulsive potential curve whereas in reality

in its ground state C2 is bound by 0.468 Rydbergs. This

represents a rather drastic failure of the MT Xamethod.

The Ne2 molecule is bound by the extremely small amount

of 0.00023 Rydbergs. This binding is usually described

as van der Waals binding. It might appear questionable

as to whether the exact Xa method let alone the approximate

Xa method with NMT correction would produce binding for

this molecule. Our findings are however that the NMT

correction results in binding for both cases. For C2

it predicts a potential curve and spectroscopic constants

for the ground state which compare quite favorable with a

CI calculation by Fougere and Nesbet (1966) and with

experiment. For Ne2 definite binding is predictedlalthough
due to statistical as well as other errors the amount of

binding can only be determined approximately. From these

______ ___~_~I~_~____I_ ____~___~ I___

results and further analysis in Chapter VI we conclude that

the NMT correction h provides a significant qualita-

tive and quantitative improvement to the MIT Xa total energy.

Further, due to the significant correlation of the results

with experiment we feel these results represent additional

evidence for the accuracy of the exact Xar approximation




2.1 Introduction

The discussion in this chapter will be concerned with

the basic Xar theory. This chapter is given as background

for and an introduction to the next two chapters in which

will be discussed first the muffin-tin approximation of

the XaL total energy and then the relation between this

energy and the exact Xa total energy. In this chapter we

wish to emphasize that in exact Xa theory no restrictions

are placed on the spatial form of the electron density p
which is used in the evaluation of the Xa total energy

. In the evaluation of the muffin-tin Xar total energy,

, a particular averaged form of the density, p, is
used. Since the form for the energy in the two cases is

different it will be shown that the differential equations

which the one-electron self-consistent orbitals must satisfy

will also be different. Besides giving some of the general

characteristics of the Xa theory in this chapter we will

derive by the calculus of variations the one-electron dif-

ferential equations which a set of orbitals must satisfy in

order to make a minimum.


2.2 Xax Theory

The Xa method as applied to atoms, molecules and solids

is covered in detail in a forthcoming book by Slater (1973).

Two review articles covering the subject have been re-

cently published, one by Slater (1972a) and the other by

Slater and Johnson (1972). The former contains a compre-

hensive bibliography of the Xa theory and its applications

since 1965.

The Xa method is a one-electron self-consistent field

(SCF) method and as such it is an approximate method and

cannot give an exact treatment of many body problems. For

example, in general it cannot give a complete description

of multiplets, although spin-polarized Xa ~calculations can

give partial information about multiple separations (see
Slater (1968)).

The Xa total energy for a molecule is expressed as a

functional of the one-electron orbitals U.. We may define

the spin up and spin down charge densities (ptf and pt) as

ptr = C n.U iU.; p+ = n.UU.v; p = pt. + pt. (2.1)

The summation for each spin is carried out over the occupied

spin orbitals for that particular spin. The occupation

numbers n. may have any value between zero and one (a fea-
ture different from Hartree-Fock in which the ni must be
either zero or one). The formula for total energy of a

molecular system in the Xa method is (in Rydberg units)

+ 12 Ui+ +dl- i +C2 +tPr~
= n U? r U(l]dlp(l r
i Q a1( rl1-RaI

1 ii2p(rl)p(r2 1 + 1 +t + +
+ drldr2 + pkrl)U a(rlr

N 2Z Z
+ -pl(r.)U c~(r )drl + -rB .c (2.2)
2 X 1 1 2 R

In this formula the first term is the kinetic energy term

and the. summation in this term is over all states regardless

of spin. The second term is the nuclear-electron Coulomb

interaction term and the summation in this term is over all

the nuclei a. The third term is the classical Coulomb

interaction term and the fourth and fifth terms are the

" exchange-corre lation terms (see Slater (1972a)) The

last term is the classical nuclear interaction term. The

exchange-correlation potential UXaf is defined explicitly as

+ 1/3
UXcl(rl) = -9t[(3/4x)pt(r )] (2.3)

and similarly for UXat'

For simplicity in the present treatment we will assume

that the orbitals for spin up and spin down electrons are

identical. This results in a so called non-spin polarized

version of the formalism. For a system like the ground

state of the C2 molecule, the main system of interest in

the present work, in which the number of spin up and spin

down electrons is equal and which near the molecular


equilibrium position in the ground state exists in a closed

shell configuration, the results of spin-polarized and non-

spin-polarized calculations are identical (i.e., they pro-

duce the same orbitals). Everything discussed here and in

the following sections may be generalized for the spin-

polarized case. The non-spin-polarized assumption does not

preclude the possibility of fractional occupation numbers

or non-closed shells. We may still think of the occupation

numbers as varying from zero to one for a given orbital

but are equal for both spins. If one does assume that the

orbitals for spin up and spin down electrons are the same

for every state and the number (fractional or not) of spin

up and spin down electrons is the same, then

pf= pf= p/2,


UXaf = UXaf = UXcL (2.4)


+t + 1/3
UXa(rl) = -9t[(3/8n)p(r )] (2.5)

The Xa total energy will then be given by

N -2Z + Prld
EXQ> C i (1,- 1 ]a- 1 +1 1
i a=1| rl-Rai

1 j2p(r )p(r2 ) t lr
1 1r 2 +1 + +
+ rl drld2 Z p(rl)UXa(rl)drl
N 2Z Z
+ 1 6 (2.6)
2c, cR


The non-spin-polarized assumption and the resulting slightly

simpler above form for will be used in the rest of
this study.

The exchange-correlation potential is linearly depen-

dent on a parameter the "a" parameter. At the present

point of development of the Xa theory this parameter is
chosen to be different for each atom varying from slightly

less than one (0.98) for hydrogen to slightly greater than

two-thirds for the heavier elements (see Schwarz (1972)).

The means of determining a for free atoms has been discussed

in the literature and is still under development (see

Lindgren and Schwarz (1972)). It consists presently of

adjusting the value of a until the Xa energy of the free

atom equals the Hartree-Fock energy for the atom.

A further point must be discussed. As indicated the

method of determining the a parameter is well defined for

free atoms. However, when atoms are bound together into a

molecule a prescription must be used to fix the value of

a used in different parts of the molecule. A very simple

and sensible way of doing this is to divide the molecule

into regions and use the value of a for the free atom in

a spherical region surrounding the atom or nucleus in the

molecule. In the region left over between the spherical

atomic regions use a weighted value of a between the atomic

a values used. If the molecule is homonuclear, as for the.

cases C2 and Ne2 which are treated in this study, this is not

a problem since the same value of a is used throughout the

Although there are similarities the Xa method is dif-

ferent from Hartree-Fock in at least two important ways.

First, the Xa theory as applied to molecules in the MSXar

method treats the separation of molecules into their atomic

constituents in a manner which leads to more physically

meaningful separation limits (see Slater and Johnson (1972)).

Thus for the H2 molecule (to take the simplest possible
case) the MSXa method "dissociates" the molecule into two

neutral hydrogen atoms whereas Hartree-Fock theory predicts

a separation limit which is a mixture of ionized and neutral

states. (This point will be further discussed in Chapter

VI where the NMT calculations for the potential curve of

C2 are discussed.) Second, the eigenvalues in the Xar method
have a different relation to the total energy in that the

eigenvalue E. for a state i is related to the Xa total energy
as a partial derivative with respect to the occupation

number n. for the state i (see Slater et al. (1969))

8. = Xa (2.7)

This relation leads to a different way of calculating transi-

tion and ionization energies, a subject to be briefly dis-

cussed later in this chapter. It also ensures that the

method satisfies Fermi statistics. That is, the system

has its lowest total energy when the levels of the lowest

eigenvalues are filled with n =1 up to the Fermi energy EF
and all the levels above the Fermi energy are empty. Rela-

tion (2.7) also allows for the possibility of non-integral

occupation numbers at the Fermi energy (see Slater et al.

(1969) and Slater (1972a)). The point emphasized here is

that the Xar theory is not just an approximation to Hartree-

Fock theory but stands as a separate and in some ways

superior theory (see also Slater and Johnson (1972)).

2.3 Derivation of the Xa One-Electron Equations

The Xae total energy as expressed in Eq. (2.6) is a

functional of the orbital functions U.. To find the ground

state energy of the system represented by we must find

the set of orbitals {U.} such that ' is stationary with

respect to arbitrary variations of the U The variation

of must take place subject to the restraints imposed
by the orthonormality conditions of the orbitals. We can
handle the restraint conditions by introducing a Lagrange

multiplier, Xi for each condition between two orbitals
i and j. Note that we are working under the non-spin-

polarized assumption and under this assumption the spin up
and spin down orbitals are identical. I am also assuming

for now that the number of spin up and spin down electrons

for each orbital is identical, although perhaps fractional.

Since the occupation number n. for an orbital i may be

fractional it is convenient to introduce a factor n. /n.1/
with each Lagrange multiplier X... (It is clear that t~he

facor n may be absorbed into the Lagrange multi-

plier X.. but this leads to awkward expressions later.)

In order to find the minimum of with respect to

arbitrary variation of the {U.} we demand

GCE~L +C 1i/2n 1/2X + +t +t-
6{ Xa j 1 l 1- j

Since the U.'s are required to be orthonormal we take the

A..s to be hermitean. (The summations over i and j are

over all states, spin up and spin down.) If is a

minimum, it is a minimum for the variation of any orbital

U hence we vary Eq. (2.8) with respect to an orbital U .
We first consider the variation of Eq. (2.6) using the func-

tional derivative notation

6 + +
6 =U 6U i(r)dr+(Complex Conjugate).(2.9)
Xa 6Utr) X

Here 6Ut(r) is an arbitrary variation of the function U"(r)
1 1
over its whole range of definition. Consider the kinetic-

energy part of in Eq. (2.6) first

6 +t 2 +t ++ -+
n~ U(rl) lU (rl)]6U i(r)drldr+(Complex Conjugate)

ii+~ + +t 2 +
= n.Ut~)6(rr )-V U.(r )]dr dr

+i~i n.t1 )[V(r-r1 )*6U.(r)]dr [ dr. (2.10)

See Appendix A.

Integrating over rl and using the Hermitean property of
the kinetic energy operator

n. Url[VU()]Urdldr +(Complex Conjugate)
6Ut(r) 1

= n. 6Ut~r)[-V2U (r~)]dr + (Complex Conjugate). (2.11)

Taking the variation of the Lagrange multiplier part in

Eq. (2.8) gives

"1/2 j1/2 + + *'U
6 ( 1 n. Cn. 1. ter)U.( )dr )

-~~~~- 3 1

=~~ n 1/2 1/ 6.(r)U.(r~dr
+~ (CmlxCojgt).(.2

sv =~jNC -2Za+ + +~d 1 2p~ )r22 +t +d
o=1 r -R 12

+ Ca cl ~~l,4/3drl}, (2.13)

where Ca = -(9/2)a(3/8T)13 Each of the terms in V i

functional of the electron density p, hence we may write

6V =V (p(rl))6Ut(r)dr + (Complex Conjugate),

= VTp~rl) + Ut(r)dr'dr+(Complex Conjugate)
II6p(r') VTUt~rr)
1 (2.14)


Using expression (2.13) above

6V p( ) N -2Z 6p(? )
T -a .i

6p(rl) +
*~~d dr. 2.5
612p( ')1

It should be noted that the second term occurs without the

factor of 1/2 which it had in expression (2.13) above.

This is so since the variables of the integral are dummy

variables and the variation actually gives two terms which

are identical. Replacing 6p(rl)/6p(r') by 6(rlr)an

integrating over rl we have that

6V (p(rl)) N -2Za 2p(r2 d
6p(f') a=1 i'-al -2 2

+ Cap( ')1/3.(16

Still working on expression (2.14) consider the term

6p(r' )
=U(r nU (2) 6 -).(.7

__ _~ _ _I~__~__~_______ ________~_~

Substituting expressions (2.16) and (2.17) into expression

(2.13) we have that

sv(r N -2Z r 2p(r )
T 1= 1 +R,I "t +r 2
=1r -Rr -rI1

4 **1/3 +t 9r *t 3 *
+ -Cap(r ) U()(rrdrr
3 1v

+ (Complex Conjugate). (2.18)

Performing the integration over r' we are left with

Nt -2Z 2p(r )
6V= .U*(r) 2 dr
ac=1 r-Ra I r-r2

+ Cap( r)13Ui ()d r + (Complex Conjugate). (2.19)

If the first terms in expressions (.1,(2.12) and (2.19)

are zero then so are their complex conjugates. Hence we

drop the complex conjugate terms and Eq. (2.8) becomes

121/2 -+ + +
0 = 6{ + C n. n. X1.. IU (r ')U.(r)dr },

"i"1 v I)E[V N -2Z 2p(r

4 + 1/3 1/2 +1/2 + +
+-Cap(r) ]ni U.(r) + / n. X U.(r)}dr.

Since this expression must be zero for arbitrary variation

of the orbital U.(r) we may take the entire expression inside
the curly bracket to be zero, hence:

2 N 2a (2-2 + 4 1/ 12
[-V+ +dr +-Cclp( r) 1/]n, U.2(r)
ai 2

=-% n.1/2.U. (r) (2.21)

We obtain the same type equation for each U .
In its present form Eq. (2.21) is a set of coupled

equations due to the mixing by the Lagrange multipliers on

the right hand side of the equation. The Xi matrix is
Hermitean and hence may be diagonalized uniquely by some

unitrary transformation. This transformation results in a

transformation of the orbitals which may be defined such

that the density remains invariant. If the density remains

invariant then so does the operator on the left hand side

of Eq. (2.21). This is easily shown. Note that the orbi-

tals U. in Eq. (2.24) occur with a coefficient n. we
must take this into account. We postulate a transformation

on the orbitals U. involving a unitary transformation C
such that (using matrix notation)

1/2 1/2 ~.2
n U' = n U. (.2

In this expression, C is an nxn unitary matrix, gl/ is an
nxn diagonal matrix whose (i,i) element is ni and U'
and U are the orbital vectors. Note that the same matrix

nl/ occurs on both sides of Eq. (.2,adE.(.2

serves as the defining relation for the new orbitals U .
Another way of writing Eq. (2.22) is

.1/2U = .n1/2U for all i. (2.23)

Since C is unitary

~+~ ~~+
C C CC =I. (2.24)

We defined p (the diagonal part of the first order

density matrix) earlier as

p" = i i .UiU = (n.i1/2Ut)n1/U) (2.25)
1 1

which we may write as

p = (nl/2 )+ nl/20). (2.26)

Now consider B where

( = (nl/2') (nl/2')

=(Cnl/2U) iSnl/2U)

=(nl/20 C C(nl/2g,

1/2~ + 1/2
=(n U) (n U)

( = p. (2.27)

This means that the unitary transformation postulated does

not alter the form of the operator on the left hand side of

Eq. (2.21).

Hence, taking each of the equations for the Uk, mul-

tiplying by Cik' and summing over k on both sides of the

_ ~

equation leads us to

dr + -Calp(r) ]n. .r
+ 2 3 i 1il
Ir-r 21

[-V + +
a=1 I r-R |

= ~Z C n1/2U.r.
k j ikhk~ kj(r


We may write the right hand side of therms of a matrix ele-


C~eik1/2 ~~~1/2
k j ikk


where X is the matrix of the A..'s. But inserting unity,

~~1/2~ ~~~+~~1/2 ~.0
(Chn U) = (CAC Cn U)..(230

We define a new matrix h' by

-' -+ (2.31)

and remembering that

we get

~~1/2~ ~1/2 (.2
(C~n U)i =h (n U')..(232

However, since C is an arbitrary unitary transformation,

we may pick C so as to diagonalize the ), mratrix

~~1/2~~ 6. 1/2 U "1/2UI
(C~n U) = A .. U A . U .


___~_ __~__________1^_1

Inserting this in Eq. (2.28) and cancelling n. em
N-2Z r2p(r2 )
2 a 4 *1/
[-V +=CI -. 1 + -r dr +-t Cap(r) ]U (r)
+ + 2 3 i

= -A!. !(r),(2.34)

and it is customary to denote -X.. as E..

For simplicity we may define a new potential VXe(p,r)

(where the p indicates the dependence of VXa upon p)

N -22 12p(r )
+ a 2 +) 4 + 1/3
Vt (pr)= dr2 + Cclp(r) .
a= 1 r-Roi | jr-r2 I


Then Eq. (2.34) may be written

[-V2 +x(.)u() Via 'r)}U (r) = EU (r). (2.36)

The potential VXa(p,r) is usually referred to as the "Xa
potential" and the Eqs. (2.36) are the one-electron self-
consistent-field (SCF) equations for the set of orbitals

{U } which make a minimum. Several things should be
noted about these equations.

First, the quantity Ei in Eqs.(2.36) is at this point

only a Lagrange multiplier. The same is true of the Har-

tree-Fock equations except that in Hartree-Fock theory the

E. are identified with ionization energies in an approxi-
mate way through the Koopman Theorem. In Xa theory the

eigenvalues Ei are related to the total energy as a

__ __~~___~_1_____1~11

partial derivative (see Slater (1972a) and Slater.and

Johnson (1972)),

Eix = Xa (2.37)

This relation causes the Xa method to obey Fermi statistics

and also leads to a different method of calculating ioni-

zation energies. In Xa theory the correspondence is made

between the ionization energy and the orbital "energy"

"i calculated for a so-called transitionn state" in which

the orbital occupation number is reduced by one-half. (See

Slater (1972a) for the derivation of this concept. Further

discussion here leads far afield from the main purpose of

this dissertation.) This is a convenient method and yields

comparatively good results where it has been tested against

rigorous methods which take the difference of total energies

of an atom or molecule and its positive ion as the ioniza-

tion energy (see Slater and Wood (1971) and Connolly et al.


Second, the solutions to the Eqs.(2.36) are self-con-

sistent solutions since the potential VXa depends upon the

orbitals Ui through the density p. The SCF cycle as it is

applied will be discussed in greater detail in the last

section of the next chapter.

Third, the solutions to Eqs.(2.36) are three-dimen-

sional functions which are not necessarily separable in a

spherical or any other coordinate system. The only

_~~ ~__ __~_I__~_ __


restrictions put on the U. so far is that they be the same

for spin up and spin down electrons (a restriction we may

relax) and that they satisfy Eqs.(2.36) above. The Xa

potential VXa (r) in Eqs. (2.36) may be a very structured
function in three-dimensions for a molecule or solid. It

may approach a muffin-tin form but in reality of course

it does not have this simple form.

__ I



3.1 Introduction

In this chapter we will discuss the muffin-tin (MT)

approximation to the Xa total energy for a molecule and

the MT one-electron equations which may be derived varia-

tionally from MT Xac energy functional. In this approxi-

mation one assumes a MT form for the electron density in

the evaluation of the energy. First in the chapter we de-

fine the MT form in terms of a spatial averaging process

and derive some of its properties. Then we will show how

the exact Xal total energy may be separated into a

MT, , and a non-muffin-tin (NMT), a, part. From

the energy functional we will derive the one-elec-

tron SCF equations which must be satisfied by a set of

orbitals {U.}in order to make ' a minimum.

3.2 The MT Approximation

In the MT approximation, the space of a molecule is

schematically divided into three types of regions.

Spherical regions surround the atomic nuclei (region I),


the whole molecule is surrounded by a large outer sphere

(region IIII, and the region left over within the outer

sphere and between the atomic spheres (region II) is called
the interatomic or intersphere region (Johnson (1966) and

Johnson (1967)). Thus a molecule like C2 may be represented

as in Figure 3.1 where the points labeled "C" locate the

carbon nuclei and the point labeled "O" the center of the

outer sphere.

It is argued on intuitive physical ground that the

potential V~c and density p for a molecule must be very close
to spherical in regions I and III and close to a constant

in region II. This is what is assumed in the MT approxima-

tion. It will be shown that under this assumption the Xa

total energy evaluation reduces to a simple form and the

one-electron equations for the orbitals reduce to a separable

central field form in regions I and III and to the form of

an electron in a constant potential in region II. Actually

it will be shown that both of these results follow simply

from the assumption of a MT density p in the evaluation of

the total energy.

The spherical average f(r) of any three-dimensional

function f(2) is defined in regions I or III as

T(2) 1 f() sineded$.(31

Thus f(r) is only a function of the radial coordinate r.

The non-muffin-tin (NMT) part of f(r), Af(r), is defined


Figure 3.1. Schematic Representation of the C2 Molecule
in the MiT Approximation.

_ __ I

Af Cr1 = f r) ICr) 3.2)

In the intersphere region the MT part of a function f(r) is

defined as

-t) 1 f (r)dr (3.3)

where the integral is only over the intersphere region and

VIN is the volume of the intersphere region. It should be
noted that TY(r) for the. intersphere region is a constant

(although for notational convenience a dependence on r

is included). The NMT part af(r) is defined in the inter-

sphere region as

af(r) = f(r) f(r) (3.4)

just as above for the atomic-like regions.
From its definition the NMT part of the function

af(r) in an atomic or the intersphere region integrates

to zero. Consider the integral of A~f(r) in an atomic


Afnrc~dr = Ifr): f(r)]dr

14r f Ir')sin9'd6'dg' sinededfr2dr.

The integral over the angles 6 and Q contributes a factor

of 47 and we have

Afb r) dr = fIr 2sinedsdadr f rc~r2sineded~dr

Saf(r)dr = (3.6)

Since this is true for the atomic-like as well as the

intersphere region, the integral of Af(2) over the whole

molecule is zero. Similarly, the integral of the product

of a MT and NMT function is also zero. Let g(r) be the MT

part of a function g(r) and Af(r) the NMT part of a function
f(r). Then

A f(r g( r dr = (3.7)

and vice versa. By these relations the total energy evalua-

tion may, in the MT approximation,be reduced to a much

simpler form.
To clarify the meaning of the averages we will

take an example. (See Figure 3.1) Consider a function

1/rl where rl is the radial distance from center I at
position R1 measured from some arbitrary origin. Inside
the sphere surrounding center 1 the MT average of this

function is trivially 1/rl. Inside the sphere surrounding
center 2 at position R2 the MT average is performed as an

angular average about center 2. In this case, the average

of 1/rl about center 2 is 1/1 R1 R2 ; that is, a constant
throughout the region about Center 2.

3.3 The Separation of the Xa Total Energyv
into MT and NMT Parts

The Xa total energy was given in the previous
chapter (Eq. (2.6) as (in the non-spin polarized case)

-+ 2 -+ + ac

-+ 1 II2p(rl)p(r2 +t -+
*drl + drlr

22 2
1 + +r- +Ua(l)r- 1 a B
+ r) r)r+. (2.6)
2 1 a1 1 2

This is obtained using some set of orbitals {U.} (not

necessarily the set which minimize )and density p
formed from these orbitals. For convenience we divide this

expression into parts

= EK + Ee + Ec + Ex + E:: (3.8)

where parts are identified:

__ __ __.______ ~____I_ I I __

Kinetic, EK

E = U (r )[- U.( )]r ;(3.9a)

Nuclear-electron, Ee


E = pr )d ; (3.9b)
ENel +C +t 12 tP~l 1
a=1 1 a

Coulomb, Ec

E j--- dr- dr- ; (3.9c)
c 2 1 2

Exchange, Ex

1 +P +f -+
Ex 2 prl)UXa rl drl ; (3.9d)

Nuclear-Nuclear, EN

1 a B
E (3.9e)

We will now treat these terms one by one in the MT

approximation to the density.

In the MT approximation to the density, P in Eq. (2.6)

is replaced by P. Here we will instead break p into its MT

and NMIT parts and consider the separation of Eq. (2.6) into

its MIT and NMT parts:

_~ ___~~_ ___~_X____I __

p (r) = (r) + 6p(r) ,(3.10)

and similarly

= + a (3.11)
XaXa Xa

where is the result obtained from Eq. (2.6) where p

is replaced by and a is the difference between

and The kinetic EK and nuclear-nuclear ENN terms
will not be affected at all. Thus we must consider the nuclear-

electron ENe, Coulomb Ec, and exchange Ex terms.

Considering ENe first we may write Eq. (3.9b) as

N -2Z-
ENe= +C at2 (P(r ) + np(r ))drl
a=1 1r Ra

We may define a nuclear-electron potential VNe as

N -2Z
V e(r ) = ( 3.13)
a=1 I 1 ac

Just as for p, we may divide this potential in any region

of the molecule and over the whole molecule into MT, V9 ,

and NMT, AVe parts

V e(r ) =Ve (rl~ )aVe (r ). (3.14)

II_____ ___I__~ __ ___ _


As shown in Section (3.2) the integral over the whole molecule

of a MT with a NMIT function is zero. Hence Eq. (3.12) reduces


ENe = V $e(rl)prrl)drl + A VNe~r)A l~r (315

This expression is exact. (The explicit forms for Ap,

AVNe, and any other such quantities will be given in the
next chapter and Appendices C and D.) In the MT approx~ima-

tion to only the first term in Eq. (3.15) is retained.

Next we consider the Coulomb term Ec. Substituting

for p we get

1 Ij2 (W (rl) + Ap (rl))( (r2) + ap (r2) + +
E dr drl

1 ( p 2P(rl) (r2 -+- + 2prlp(r
=2 dr2drl +l 2prlPr dr2drl
rl2 rl2

1 -t 2A- Apr
+ III2pr)Pr dr2drl (3.16)

There are twJo cross terms from the multiplication; but,

since they are identical by interchange of dummy variables, the

second integral has a factor of one instead of one-half as

originally. We may rearrange the first and second terms
into a similar form

1 c + (1 2p(r2 3 pr
Ecp~l)dr2dl+ pr) dr2dr
rl2 rl2

1 I ( 2Ap(rl)ap(r2 *
+ dr2drl .(.7

We define a new potential Vc(p, rl) (the p is included to
indicate that Vc is calculated using the MT density)

Vc p'rl) = pr dr2 .(3.18)

This potential itself may be separated into MT and NMT parts,

Vc(p,rl) and CV (p,rl) respectively

Vc(Prrl)= Vc(p'rl) + AVc(p~rl) 3.19

Inserting these into expression (3.17) and taking integrals
of MT with NMT functions to be zero we are left with

Ec 1 p~r)Vc(p,rl)drl + (r)c(p, rl)drl

1 26p(r )Ap(r2

Again, this expression for Ec is exact but only the first
term is retained in the MT approximation to . Note

that V (p, r ) is a radial function in an atomic-like
region and a constant in the intersphere region. Further

Vic(p,rl) may be evaluated at any point in the molecu e in
terms of one-dimensional integrals. Thus the first term in

Eq. (3.20) is in effect a two-dimensional integral, the
second term in effect a four-dimensional integral, and the

last term is a six-dimensional integral.

Continuing, we have the exchange term of the total

energy Ex which we handle in a different manner. We write

Eq. (3.9d) as

Ex =~ C p(r ) /dr (3.21)


Ca~ = -(9/2)ar(3/8r)1/ (3.22)

+ 4/3
Adding and subtracting Ca i rl) to the integrand of
Eq. (3.21) we have

-+ 4/3 + 4/3 4/3 +;


where the first term is that retained in the M~T approximation


~ I

Now we have completed the separation, for a given

set of orbitals {U.} and density p, of into its parts

and 6
= + a (3.11)


~ ~ + 2 +i j fI1- +2 t
= n U (rl lU (rl)dr prVNe(rl drl


1 Ij2bp(rl1:)np(r2) -t 4/3 4/3 -+
+ 2 dr 2drl + Ca i(p(rl W(rl )drl


Again, we emphasize that no approximations have been made

in Eq. (3.11)--only a separation.

1~- + + 4/3 -+
+~~~ ~~~ 2 pIlVIr~dl+C ~l drl

2Z Z


b = Ap~r )nVNe (rl)drl + Ap(rl)Abc(p, rl~drl

3.4 The Derivation of the MT Xcr One-Electron Equations

In the MT approximation for the Xar method the orbitals

U. from which p is formed are determined as those orbitals

which make expression of Eq. (3.23) a minimum. We

designate these orbitals as "U." to distinguish them from

the "U." which in Chapter II were the orbitals which made

a minimum. The two sets of orbitals although perhaps
nearly identical are not the same. We will derive the one-

electron equations for the set {0.} in the manner of Sec.

(2.2);that is using functional derivatives. For simpli-

city in the following derivation we will leave out the

Lagrange multipliers and the complex conjugate terms

since they are handled in an exactly analogous manner as

they were in Section (2.2). We wish to concentrate here

on the difference between p and p as well as on the form

that the potential will take in the one-electron differen-

tial equations.

In forming the variation of ,' the form of the
kinetic energy term will be the same as for the~ exact

Xcr method and the nuclear-nuclear term does not contribute

anything. Hence

6~ = n. Ur)[- V U.(r )]dr + 6VT (3.25)
Xa 1 1


V~~~~~~ = ~ 7er)r (rl Vc(pIrl)dr l

-- I I

+ Ca p(r ) dr a,.. (3.26)

As indicated previously the functions p, V~e, Vc, and p are

radial functions in the atomic-like regions (regions I and

III) and constants in the intersphere region. Hence we may

divide VT into a sum over the various regions performing
the angular integration in the atomic-like regions and the

volume integration in the intersphere region:

S1 + 4/3 2


1 -- 4/3
VIN PIN Ne,IN 2 IN c,IN (P INP 3.7

where a=0 corresponds to the outer sphere region, VI is

the volume of the intersphere region, and the subscript

"IN" indicate volume averages over the intersphere region.

The variation of VT may be written as

69 =1 7 6U(r)dr .(3.28)
6U. (r)

Since VT depends on U. Only through 7we may apply the chain
rule of differentiation generalized to functional derivatives

to find 6V Hence

N ** 6 (1-
67; = 6, Ui-*(r)- 4x12lN~l pr)cpr


+e4/3 2 2 +
+~(" Ca (l rdrl r' dr' dr
6U. (r)

~* 6 1 -4/3
+i 6U.(rT) ---V6 I[p NVJe + p V I+ Cp ]

6p -+
6U. (r)


Taking the terms one by one we have first the nuclear-

electron term in an atomic region

6 -2 6 rl -2
1 N 11 1 ,Ne 1r~ 1 1l

sp~r= r1jB~6 P(rl-r)Nelrdl

(IN Ne,IN) =Ne,IN (



In taking the functional derivative of the Coulomb

part of VT it is easier to return to the earlier form of
this term which is given

4x p1lVcprlrdrl

SI2p(r )p(r2 )t-
drdr (3.32)
2 r1 2 1l

Then the functional derivative is

6 1 jl2p(rl)p(r2 +

-~d dr d
2 r 2 1

2p(r )

= 6(r- r')dr drl


since rl and r2 are dummy variables.

V (p,r ) we have

By the definition of

2p rl'p'r2 + +
'12 dr drl

6 1

= ~ ~ ~ -- 6( 'V(pr)r.


12 2 (C

- I I

Doing the angular integration if in an atomic region and

multiplying and dividing by 4n, and doing the volume integra-

tion and multiplying and dividing by VI if in the intersphere

region we get

6 1 jl 2p(rl)p(r2) -+ + B P''1)frao
2 dr2drl1
6p(r') rl2
1 V NVcI for the


The last term of 9T to consider is the exchange term. For

an atomic region it is

6j p~4/3r 2 4-1/3 2
4xC 6r) rdr x C rl) 6r-r)lr

4 -1/3
=4r- C p(r')
3 a


and for the intersphere region it is

6 -4/3 4 -1/3
CUV p= V Cp .

Finally, collecting terms for the functional derivative

of VT in Eq. (3.29) we have


N~* + -41/
6VT =E 6U. (:)47i[Vi~ (r') + V (p, r') + -C (r') /]

6p(r') 2 +t
x ~, rt' dr'dr

~* ~4 -1/3
+ 6 .() [ + C) -Cp ]
I 6Ur)VNVIN NeI Vc,IN p 3 ar IN

x - dr (3.37)

Now, we must consider the term 65/6U.(r)in the atomic-like

and in the intersphere regions. In an atomic-like region

67(') 1 ~* ~ +

GU (r) U.(r) 4x .

1 1 1 i(-

1 n.6r' r)U.(r )sin6'd6'd (3.38)

Note that we may combine the angular integration over the

primed coordinates with the radial integration over the

primied coordinate in Eq. (3.37). Now consider the functional

derivative of p in the inter-atomic region:




1 n.U C')U.(f')d?'


~*i r

1 ~ i(~

1 1,


Substituting the expressions ofEqs*(3.38) and (3.39) in

Eq. (3.37) we have

- N
6V = 1

~* 4 r1/3]
ii6U.(2 )[V (r') + V (P r) + C Wr'

x n.6(r' r)U.( )dr 'dr

~*+- 4 1/3]
1ivr [Ne,IN Vc PIN 3 aC IN

x n.U.( r)dr .


We may now recombine the integration over all the regions

-- I

and perform the primed integration in Eq;. (3.40) to get

67T = n.62() [ (2 + V (p, r) + 4 C a Pr)1/3U.(r dr.


In this expression it is understood that p, VNe, Vc, and
p are radial functions in the atomic-like regions and

constants in the intersphere region. We may conbine the

potentials above into one potential VX,;

-t + 4 -+ 1/3
V (p, r) V (r) + V (p, r) C p (r) (3.42)
Xa e 3 a

Finally we combine Eq. (3.41) and with Eq. (3.25). In

doing so we include the diagonalized Lagrange multiplier

.i which was left out of Eq. (3.25) but understood. (Note

that we use 7. here in contrast to the E. of Eq. (2.36) ):

~*2 + -t t -++
0 = 6U.*(1) [-V + 9 (pr) E.] U.(r) dr .(3.43)

Since dU. is arbitrary it implies that the rest of the

integrand in Eq. (3.4 3) is zero:

[-V +V xa(p, r)]U.(r) =i~ e.r) (3.44)

The solutions to this equation, U (r), occupied from that of

_~ ~ _^ ~_______ II____ __ _

the lowest energy up to the Fermi level so as to be con-

sistent with the total number of electrons in the mole-

cule are the set of orbitals (U.} which make as

given in Eq. (3.23) a minimum. It should be noted that the

potential of Eq. (3.44) is of muffin-tin form and hence

the solutions to Eq. (3.44) are separable in spherical

coordinates. Both of these facts follow from the assump-

tion of a muffin-tin form for the density. The E. eigen-

values have the same property as the c. for the exact Xar

method; namely (see Eq. (2.37)),

E. = (3.45)

The orbitals U. are SCF orbitals as is implicit in Eq.

(3.44). The details of the SCF cycle will be discussed

at the end of the next chapter (Section 4.4).

3.5 The MT Xa One-Electron Orbitals and
the Multiple-Scattering (MS) Formalism

The MT form of the Xcl theory is that used for compu-

tation since the equations for the U. are separable in

spherical coordinates and hence readily solvable as one-

dimensional second order differential equations. Histori-

cally this fact was noted first; that is, that if the poten-

tial were
solved using w~ell-know:n techn~iques. As shown above in Section

(3.4) this result follows simply from the assumption of

a MT density which transformstle energy function from

to . In the following paragraphs of this
section a brief discussion of the multiple-scattering

formalism will be given. For a more complete exposition

consult Appendix B, Johnson (1966), and Johnson (1967).

With the separable form of the differential equations

for the U. (Eq. 13.44)) we may use the partial wave expansions
of multiple-scattering theory in the three regions of the

molecule. Within an atomic sphere a or the outer sphere

(a = 0) we may expand the orbital function Ui for state
for eigenvalue ,..

aq = ClanR (ru E)Y m(r ) (3.46)

where the Cl are expansion coefficients to be determined,
the Y~m are real spherical harmonics, and the atomic-like

RR satisfy the central field differential equation

R(I1 + 1)
[- 2 (r~c d )+ Z + V ,(p, r ) E.]R (r ,'E ) = 0
r dr dr r


In the intersphere region the potential VXa(p, r)

is averaged to a constant VI in the MT approximation and the
differential equation to be solved is

IV~ + e. V JU.r) = 0. (3.48)

The solution to this equation may be expressed as a

multicenter partial wave expansion which for c. < ViI is
given by

U. (r) =C i Aak (1) (Kir )Y (r )C + Alot Kr)Y(

ar=1 R,m R,m


where the Atm are expansion coefficeints to be determined,

the Y~m are real spherical harmonics, ra = r~ ital (where Ba
is the location of the aith nucleus), and i. (c. 1/2
1 1 VIN)
The k(1) and it are respectively the modified spherical

Hankel function of the first kind and the modified spherical

Bessel function. For i > VIN the expansion for the orbital
in the intersphere region is done respectively in terms of

the ordinary Neumann and ordinary Bessel functions.

To find the eigenvalues and the A and C expansion

coefficients one requires the orbitals U. and their first
derivatives to be continuous across the sphere boundaries.

This leads to a secular equation from which the eigenvalues

Ei and the A and C coefficients can be determined. For the
details refer to Appendix B.

The Xa approximation, the muffin-tin approximation
to the density, and the multiple scattering formalism as

outlined above together constitute the multiple-scattering

Xa (MSXa) method of molecular calculation as it is cur-

rently implemented.

3.6 Explicit Expressions for the Components
of the NMT Correction a

The expression for the NMT correction A is given

n = Anpcr )AVNe(rl~drl + ApP(rl)AVclpIrlldrl

1 I 26p(rl)Ap(r2+
+1 dr drl

+4/3 +t 4/3 -+
+ C (~rl -p(rl) )drl. (3.24)

Of the terms contained in the integrands, the meaning of

Ap(rl) is clear and its explicit expression in terms of
the partial wave expansions of the multiple-scattering theory

is given in Appendix C. The third term whose integrand

contains the second order product Ap(rl)ap(r2) may be argued
on intuitive grounds to be small compared to the other

terms. However, we have actually performed a calculation

of this term and found it to be quite small. The result

of this check calculation will be given in Chapter VI and
this term will not be considered further here. (The manner

of performing the six dimensional integration for this

term, particularly the method used to handle the 1/r1

singularity, is given in Appendix F). The meaning of

the integrand of the exchange term is transparent presenting

no conceptual or calculational difficulties. The remaining

two terms--AV~e(r) and aVc(p, r)--are more complicated and
will be derived in detail.

Deriving the expressions for these two terms in the

three types of regions of a molecule--atomic, outer, and

intersphere--is a tediomibut straight-forward exercise in

electrostatics. Of the two, aVNe(') is the simpler and we
will give the derivation for its explicit regional expressions

in this section. The Coulomb term nV (p, r) is more involved

and we will give in this section the explicit expressions

for this term in each part of the molecule. It will be

shown that some of these terms are exactly analogous to

the nuclear-electron terms and these need no derivation.

The other terms of nV (p, r) have no counterpart in the
nuclear-electron terms and the derivation of these terms is

relegated to Appendix D.

The NMT Nuclear-electron Terms--aV (e

The explicit definition of VNe( ) is given, once again

N -2Z
Ve (r) = .(3.13)
Ne~ r-+ I

The function nV ,(r) is defined by

aV ( ) = Vr ( ) Ve ( ) .(3.14)

Hence our task is to derive expressions for V e( ) in the

three regions of the molecule. Repeating,the muffin-tin

average in an atomic or the outer sphere region is an angular

average--an average performed in the region of interest.

The muffin-tin average in the intersphere region is a

volume average over the intersphere region.

We may divide VNe(r) into four types of potentials

depending upon where the point r is located and where the
nuclei are located: VaL~aL VaC,B O~a an ICa hs
Ne 'Ne 'Ne Ne
are defined:

VaNa The potential of interaction between an
Ne electron in the oth atomic region and the
ath nucleus located in that atomic region.

Va,8' : The potential of interaction between an
electron in the ath atomic region and the
Bth nucleus located in the 8th atomic

Ne :The potential of interaction between an
electron in the outer sphere region and
the ath nucleus.

VIN, a: The potential of interaction between an
Ne electron in the intersphere region and
the ath nucleus.

We need to find the muffin-tin average for each one of these


In the following the nuclei are located by vectors


R[(a=1,..., N) measured from an arbitrary origin "O". The

point r in the potential V~Je(r) may be located in any of

the three regions: 1) In an atomic region "a" it is

designated rn where rc is the vector from the nucleus al to

the point r. 2) In the outer sphere region "o" it is

designated r where r is the vector from the center of the
o o
outer sphere to the point r. 3) In the intersphere region

"IN" it is designated rI (sometimes ro for convenience)

where rI is measured from the center of the outer sphere.

Beginning with V~.C we have

Ne -


a (3.50)

Since this is already a spherical function in the ath

region the angular average makes no difference and

Ti ra~r) =7 Vaa(r) (3.51)
'Ne 'Ne

and hence

6V~(r) = (3.52)

For Va, we have

V",58() = ,(3 53
Ne -+ + .5
r -RI(


where r is in the oth region. We may use another set of

vectors which are illustrated in Figure 3.2(a):

r R = r C1+ R -

=r -R .(3.54)

R ~i where r', C < R ( is given

The expansion for 1/lr -


03 m=+R r
1 4r ar Y m(r)Y m(IR ) .(3.55)
+-t R+1 E ma
Sa RaB I =Om=-R R R,

In this expression Ra is of course fixed and r ctis

variable. The angular average is performed in the ath

region and in performing this average we note that




and from Eq. (3.55) and Eq. (3.56)



IYem(r)dG = 4n6(R)6(m)Y .

-a, r -+B 1 t dG ,
'Ne -
4x Ir Rcf |



( b)

Figure 3.2.

Vector Relationships Used for the Derivation
of (a) V ,and (b) VO 4
INe ~e'

I _


Br -2Zt -2ZB
ide **
Icr Ra Rcr


Then,within an atomic sphere a the total potential AVaT(

is given by

dV~e(r) =t -2Z( ,1 _
B=1f -a R g

1 ) .


For Vroea(r) the point r is in the outer sphere


'Ne -+ + -
|r -R


where the vectors r and R are defined in Figure 3.2(b)
o to

+ + +
r = r- R
o o

-+ -++
and R r= R R .


-+ +t + +
Again we expand 1/ r -F R where ) ro > R

(a m=+R R
1 4n to *
1 1Y m(IR )Y (r ).(3.63)
Ir =0 m=-R r
Io- B to o9+

Taking the spherical average in the outer sphere region we get


-2o4a I 1t dGt (3.64)


vyo(.a a- (3.65)

Kence in the outer sphere region AVle~r isgienb

avo ( )= = -2Z, ( 1 1) (3.66)
Ir -Ra r0
a=1 o t

The derivation of VN (r) in the intersphere region

is a bit more involved due to the odd shape of the inter-

sphere region. We must perform a complete three-dimensional

integration over this region and this is accomplished by

first integrating the potential over the whole region inside

the outer sphere and then subtracting the integral of the

potential over each atomic region. The potential VIN'"a~

is defined

VINq a + (367
VNel (' *6

where r is in the intersphere region. First we will integrate

this function over the volume inside the outer sphere. To

do this we use the same set of coordinates as for the case

Vo,a aoeol o 2)myb eso rae hnRI
Ne abv nyn Iol ma tosso raerta
We get (skipping the by now familiar expansion step):

- 4r

for rI ro>l jRao



1 4-i

Sro -Relo toO


Hence, the volume integration of VIa inside the outer sphere


rI2dr +4nj r dro }

-2Z di -2Z { 4r (
r -r R |R
o toL ao 0


where bo is the radius of the outer sphere. Continuing

2 2
-22 arb R
-23 t dr = 2 37
r Rao 2 6

o3 o

In integrating 1/lr R ( over the atomic regions
o tco
there are two cases to consider: one in which we.-integrate


I ~1 a
ro t ol

for (r3 I < R3 c

over the ath sphere where the oth nucleus is located and

the-other type case where we integrate over; the other (N-1)

atomic spheres B. Within the oth sphere r -R=r and

1 -t a
-2Z a d r =-2Z (3.72)

Within the EBth sphere we decompose the vector r as

before in Eq. (3.54) (See Figure 3.2 (a))

r rB + Rg R

=r R *' (3.72)


1 -+ 41T 2 g
-2Z. II dr =-2Z rd
R aB1 aB 0

=-2Z (3.73)

where bg is the radius of the Bth sphere. Hence with the
result of this last equation, Eq. (3.71), and Eq. (3.72)
-IN +t
we have that VNe(r) in the intersphere region is given by

2 2 2 3
N b R b N b
-IN 4r o -2o, a iO i 1 C
VNe IN 3
V 2 6 2R
I:N =1 B=1 RuB

where VI is the volume of the intersphere region and is

given by

V _47 ( b3 b3 ). (3.75)
IN o ai

IN +t +
Then aVNe(rlN) at any point rI in the intersphere region

is given by

N -2Z
IN + e -IN
AV (r)=1-V ( (3.76)
N ea=1 IIN acl

The NMT Coulomb Terms--aV (p, r)

The expressions for the function AV (p, r) in

different regions of the molecule are somewhat more complicated

than those for aVe (r) The basic integral for V (p, r)

is given by

V (pr) = dr.(3.77)

The point r may be in any of the three types of regions of


the molecule--atomic, outer, or intersphere. In each case

of course the integration over r' is over all space. For

a given point r we may divide the total potential at r

into a contribution from the region in which 2 is located

and contributions from the other individual regions of the

molecule. To find aV (p, rZ) at each point r we must first

derive V (p, r) in the region in which r is located for

the contribution of the integral from each region of the

molecule. The derivation proceeds in the same manner as

above for aV ,(r) with the basic difference that instead of

a potential produced by nuclear point charges we now deal

with a potential produced by spherical radial charge

distributions in the atomic and outer sphere regions and by

a constant charge throughout the intersphere region. Just

as we divided the nuclear-electron term into four types

we divide the Coulomb term into ten types of contribu-

tions. ~(For a reference covering many aspects of this

section see J. W7. D. Connolly & J. R. Sabin (1972)).

In the orth atomic sphere the potential Va may be

divided into parts as follows

V" = V"'a + V",$ + Va,o + V",IN (.8
c c c c c

where the terms on the right hand side correspond to contri-

butions to the potential at a point in the ath sphere arising

from MT chage distributions

_ _

--in the ath sphere itself V"'a
--in the other (N-1) f3 atomic spheres Va,

--in the outer sphere Va'o
--in the intersphere region VN

Using the same kind of notation and definitions we may

divide the potential VI in the intersphere region into


VIN ~ IN,a + IN,o + IN,IN (.9

Similarly for Vo the potential in the outer sphere region

Vo= V~ +V~ VN.(3.80)
c c c c

In the following we will give the explicit expression

for each of these ten types of potentials and also the MT

form for each, AV- of course being the difference of the

two. At the end of the list we will indicate the parallel

nature of some of these MT averages to those of the nuclear-

electron case which were derived above. For the rest, the

derivation of the MT averages will be given in Appendix D.

Before giving the list of expressions we give here the

definitions (some repeated) of the important symbols which

are used:

b_ -- radius of atomic sphere (a=0 for outer sphere)

go -- total charge contained in sphere a

R --distance between two atomic sphere centers
aS a and B

R --distance of atomic center a from outer
sphere center

PN- average charge density in intersphere region

VIN- volume of intersphere region
0 () =4v27 r) radial charge density in ath

sphere. (3.81)

Since there are ten types of potentials they will begiven in

the order of potential, MT potential, and NMT part of

potential without comments or other divisions:

r b -
a to~ (r')
Vaa+ _2t r)r 2 a dr'
c a a ,
r r

ca a cf aLCC

gya"(2 )= 0 (3.82a)

N N 2Q2
C Vca-+ Ir R







AVa ( o) = T

20 ( 1 1 )
r -Ra | Rol

(3. 82b)

0 (r')
I ouT dr'

VC1o(r) 2

VfCL o (r ) =C Vo (o

AvaL'o (S ) = 0
c a

(3. 82c)

N b
a, (IC) qlIN + 2 2 1 -+ +taO 2 1 2 2 C
c a I o a to a 3 R I
B=1 a a

8 1



aINb + 2 12
V (r ) = 4 xp (~b b R.
c a INo at 3 ro

a,"'IN *;, 1 .* 2 2 1 2
c I 3 a o o 3 a

23 1 1
b ( -
8=1 +a RaBI as



a=1 al=1

r -R
IN oa

N N b
INa+ 8n e( o
V 2
al=1 IN a=1







2QI 1 471[ o
r I- Rea VI 2

oa )


A V I~c(rI- )= C
a=1 a=1






a (r')
ouT dr

cN IN -

IN'o +C

--IN, o
V c r )

c IN ~



It -)
r -
IN oai




c C IN 4.rIN C o rIN

"gl"' (r N) = 4v~ (b2 4w






3 o
b [
a 2



1 v


8xr N
-- -I
ILN ar=3

- V~MIN,IN~r)

AV (INr ) =VININ~r N)

~ cO', o L
a=1 cr=1

N -- -N 2
Vc r o)= Z
a=1 a=1


0 C


A~~tr)= Q( 1 1 )
r -R ro
cr=1 ar=1 o oa


~o b

a (r')


a (r')dr' + 2

o,o -+ o0o
Ve (ro)= Vc' (ro)

avPOI~o)o= )=0


4vrb N
VOIN~ + o
7 (r)= 2p -
o ac=1

4nb N
V,~ (r )= 2pIN
o a=1

3|r R

4rb 2pV
3r r
o o

av (r )~i=-2p;, ~- ( L 1) (3.82j)
3 Ir -R r
a=1 o oa O' o

In examining these terms it will be noted that the

expressions for -7' (an V,0~ a,
c c c

-IN a --o,a
, V and V
c c



-a,a -0,8B IN,a
exactly parallel the expressions for VNe VNe Ne'

and Wo~a The difference is that -Za is replaced in

the latter case by Qa in the former. Of the remaining terms
AVIN,o and 6Va'o are zero and are relatively simple expressions
c c
to derive. This leaves three NMT potentials whose MT average

is relatively complicated and which have no counterpart

in the nuclear-electron potentials derived above: AVa, IN

AVIN,IN and AVoIN These three potentials, their MT
c c
average,and their NMT parts are derived in detail in

Appendix D.

This completes the survey of the component terms

which form the NMT correction A to the MT energy
In actual calculations they are not necessarily used in their

present form and many are combined together for convenience

(especially the counterparts of the nuclear-electron and

Coulomb terms).



4.1 Introduction

In Chapter II we derived the differential equation

which must be satisfied by a set {U.} in order to minimize

the Xar total energy . In Chapter III we showed how,

for a given set of orbitals {U }, could be separated

exactly into a MT part and a NMT part A; we also
derived the muffin-tin differential equation which must be

satisfied by a set {U.} in order to make ' a minimum.
In this chapter we shall examine the relationship between

and KEXa> more closely in order to elucidate (which
is our main purpose) more fully the effect of the muffin-tin

approximation to the total density in the total energy eval-

uation. This is accomplished by using the set {U } to

evaluate the non-muffin-tin correction A to the muffin-

tin total energy . We shall show that this approximation

to is analagous to performing Rayleigh-Schroedinger
first order perturbation theory in ordinary quantum mechanics

and hence our result for is correct to second order

in the differences between the set {U } and the true Xa orbi-

tals (U.}. We will discuss the extension of this method to
excited states and will also show the analogous second and

~ ~

higher order corrections to do not exist in our pres-

ent method. We shall indicate how a first order approxi-

mation similar to that implemented in the calculation of the

approximate total energy may also be used to calculate

corrections to the muffin-tin eigenvalues, although no

calculations of these corrections have been performed in

the present study.

At the end of the chapter we discuss the SCF cycle,

at which points in the hypothetically exact SCF cycle muf-

fin-tin approximations are performed, and how the approxi-

mate method for calculating fits into the present

muffin-tin scheme for calculating
Only the development of an efficient multi-dimensional

multi-center numerical sampling and integration scheme has

made possible the practical evaluation of the non-muffin-

tin correction A since it involves up to six-dimen-

sional numerical integration. This sampling and integra-

tion scheme is discussed in detail in Chapter V.

4.2 The Explicit Relationship Between and

In Chapter IIIwe showed that could, for a given

set of orbitals {U.} and density p formed from the set {U.},

be broken into muffin-tin and non-muffin-tin parts

= + a



1 + + + +4/
+~~~~~~ p~lV(,r)r a rl) drl

2Z Z
+ r (3.23)


n = bp(r ) nVNE(r )dr + bp(rll) AVclp, rl~drl-

1 II2ap(rl ) p(r2 + +4/
+ 2dr2drl C. I[p(rl~4

-p(r ) ]drl (3.24)

The definition of ap has given as the difference between

the density formed from the orbitals {U.} and the muffin-tin
density (see Eqs. (3.1) and (3.3)) formed from the same


ap = p P (3.10)

The terms OVNE and a'Vc6) are defined respectively in Eqs.
(3.14) and (3.19).

_ _

In Chapter II we showed that the set of orbitals

{U.) which make a minimum are solutions to:

12 Xt

[-V2 + V ,(p, r)] U.(') =E.U. (r), (2.36)


+ a 2 4 -+ 1/3
V (p, r) =C +2 dr + -C pr)
xa ~ ~ ** r * .* 3 a
a=1 2

In Chapter III we showed that the set of orbitals {U } which

make a minimum are solutions to

2 --t ++ -+
[-V +V (p, r) ]U. (r = ..(),(.4
Xa 1u 1 r) 134


- + + 4 +t 1/3
V (pr) Ve (r) + V (p, r) ~)(.0
VxaP Ne ccPr (3 a0

with the definitions of Ve and V (p) given in Section 3.3.

Ideally we would like to solve Eq. (2.36) for the

set (U } and use these orbitals to form . However, in

practice this has not been done up to the present for two

reasons: 1) we have not had a practical means for solving

Eq. (2.36) for the (U.} and 2) the integration in Eq.

(2.2) where is defined involve up to six-dimensional

numerical integrals. Concerning the first reason as pointed

out in Section 2.3 we have at present no practical means

available for determining the solutions to Eq. (2.36). This

is of course not to say that the solutions can not be

obtained. Concerning the second reason, we will show in

thenext chapter how a practical scheme may be devised to

perform the numerical integration involved in a direct

evaluation of or a.
xa xa
In the present implementation of the X_ muffin-tin

multiple-scattering method (MSXa for short) we solve for

the set {U.} and the quantity . Hence, although we

can not obtain the (U 1 to minimize exactly, we can

easily obtain the {U.} which minimize ; and, use
them to calculate 6 and therefore in a first
xar xa
order approximation. As indicated before, for some molecules--

CH4 and Li2 for example (Johnson, et al. (1972))--the total

energies calculated from with the {U } are quite

good and it has been found that in even the worst cases--C2
being an example--the relative error (estimated) in the

total energy is always small. This being the case, we might

expect that using the set {U.} to calculate the non-muffin-tin

correction a to the muffin-tin energy should give

improved results; that is, results somewhat closer to the

true Xa energy. At any rate this should give a good estimate
of the size of the non-muffin-tin correction a and
therefore indicate whether the exclusion of this term is a

major source of error in current MSXa calculations.

This approximation--using the {U.}to calculate


--is quite analogous to first order Rayleigh-Schroedinger
perturbation theory. In the usual quantum mechanical treatment

(see, for example, Merzbacher (1965)) one has a Hamiltonian

H which may be conveniently separated into a zero order

Hamiltonian Ho and a perturbation V:

H = H + V (4.1)

Further one has a solution for the ground state Yo which

minimizes the expectation valve of Ho

'o = (4.2)


6 = (4.3)

The first order correction q1 to the zero order energy Eo
is the diagonal matrix element of the perturbation V

El = (4*4)

Hence the total energy E in the first order approximation is


C = cO + El (4.5)


E = .


_~__ ___ ___~___ _II__ _~^

In analogy, in the present case we have instead of an

expectation value a functional which we may separate
into a zero order part and a perturbation a:
xa xa

= + A .(4.7)
xa xci xa

We have a set of orbitals {U.} which minimizes the

functional :

6 = 0 (4.8)

Using the set {U.} we may calculate n which we see to be

analogous to the first order correction in Rayleigh-Schroedinger

perturbation theory. Hence we may calculate with the

set {U.)

+ n (4.9)

which is seen to be analogous to Eq. (4.6). If we take

the difference between each hypothetical exact X_ orbital

U. and the corresponding muffin-tin XQ orbital U. to be 6U.

then the energy of Eq. (4.9) is correct up to (SU )2 effects

for all the orbitals. This is just to say that the total

energy is not very sensitive to details of the orbitals--a

well-known fact. This of course does not imply that the indi-

vidual parts of the total energy in this first order

approximation will be accurate to the same order as the total

_ _


energy. In fact, since we are using the {U } to calculate

the kinetic energy we may expect that the first order ap-

proximation will not satisfy the Virial Theorem although

an exact Xa calculation for a molecule with the same value

of a throughout does satisfy this theorem.(See Slater (1972b)

for the proof of this result and further discussion~.) We

will return to this point again in the discussion of the

results for C2 in Chapter VI. A more accurate way of test-

ing the quality of the orbitals would be to calculate the

forces on the nuclei using the Hellman-Feynman Theorem

since the forces depend on integrals in which the density

occurs linearly. However, our main concern in this study

is the total energy.

The total Xa energy for a single-excitation excited

state is calculated, self consistently (see Section 4.4),

with an electron removed from one of the occupied orbitals

of the ground state to one of the unoccupied higher energy

orbitals of the molecule (see Slater et al. (1969), Slater

and Wood (1971),r and Slater (1972a)). The total energy of

the excited state is a functional, just as for the ground

state, of the one-electron density matrix; for an excited

state it is defined with one electron removed from a ground

state orbital to one of the higher energy orbitals. The

orbitals self-consistently obtained minimize the energy

functional defined by this one-electron density matrix.

Double excitations are treated similarly. (Considering

for a moment the orbitals obtained, as an example, for a

singly excited state, it should be noted that these orbitals

will be orthogonal to each other but not necessarily orthog-

onal to the orbitals obtained for the ground state Xa energy


If one were to draw an analogy with a configuration

interaction (CI) calculation it would be seen that within

the Xae (that is, pl/3) approximation the XaL total energy

for the ground or an excited state corresponds to a diagonal

energy matrix element obtained in the CI expansion. In

contrasting the diagonal energy elements of a CI expansion

to the Xa energy for an excited state or the ground state

a particular qualitative difference exists in that in the
Xa method the energy is calculated self-consistently whereas

in a CI calculation the energy matrix elements are basis set

dependent since a limited basis is always used. (This self-

consistency allows the Xa theory to satisfy the Hellman-

Feynman and Virial Theorems. See Chapter VI for further com-
ment on this subject.)

In the present method described in this chapter we

calculate (both MT and NMT parts) with the set of

orbitals {Ui) which minimize . Since the NMT correction

A is a functional of the one-electron density, the NMT
correction for an excited state is calculated in exactly the

same way as for the ground state. There is no difference.
Since we have shown our present method to be equiv-

alent to a first order perturbation calculation the question

naturally arises as to whether the method may be extended to

higher orders. The answer is no. No, because to extend the

method to higher orders would require the calculation of

off-diagonal matrix elements of the "Hamiltonian" which,

as we have indicated above,.do not exist in nor are con-

sistent with the basic Xa theory. From this we conclude

that any fundamental improvements to be made in the cal-

culation of exact Xac energies beyond our present approxi-

mate method must come from improvement in the methods of

orbital calculation (see Chapter VI Section 6.5 for some

further discussion).

For our main test calculations we picked the C2

molecule since for this molecule the muffin-tin.approxi-

mation gives a rather bad result for the binding energy.

Muffin-tin calculations give the result that C2 at its

experimental equilibrium separation is unbound by 0.90

Rydbergs whereas in reality it is bound by 0.47 Rydbergs

(the total muffin-tin energy of C2 to three figures is

150 Rydbergs). We estimate then that the muffin-tin re-

sult is inaccurate by about 1.2 Rydbergs. We may expect

the C2 molecule to be a bad case for two reasons. First,
it is a homonuclear diatomic molecule which means that

75% of the volume inside the outer sphere is .included in

the constant potential intersphere region (see Figure 3.1).

The muffin-tin approximation is more severe in the inter-

sphere region since at the equilibrium separation 2 3 of

the 12 electrons of the C2 molecule are in that region
where the density and potential are both averaged to a

constant. Secondly, the C2 molecule has a oU and nU

valence states. The nu state is a highly structured state

with 42 percent of its density the intersphere region.

However, the muffin-tin approximation averages away this

fine structure and its effect on the total energy. Hence

we feel that C2 is a good test case which should stretch
the first order approximation to the limit of its validity.

To anticipate the results given later, the calculation of

the non-muffin-tin correction A for C2 has resulted
in a significant qualitative and quantitative improvement

in the potential curve for C2. In addition a point cal-
culation at the approximate experimental equilibrium sep-

aration for Ne2--a molecule qualitatively different from

C2--has corroborated this trend for the non-muffin-tin
correction. Hence on the basis of the relative insensi-

tivity of the functional to changes in the orbitals
and on the basis of the significant improvement of the re-

sults we feel that the non-muffin-tin correction A to

the muffin-tin Xa total energy represents a signifi-
cant and practical improvement to current MSXa calculational

methods. Further discussion with complete results will be

given in Chapter VI.

4.3 First Order Corrections to the Eigenvalues

We may also calculate a first order perturbation

correction to the eigenvalues of the muffin-tin single-

electron Xal equations. (This of course may be done for

transition state eigenvalues as well as any other con-

figuration with integral or non-integral occupation numbers.)

If we compare the exact Xar single-electron equations

of Eq. (2.36) with the muffin-tin single-electron Xa equa-

tions of Eq. (3.44) we see that the two differ in the po-

tential. We may take the difference between the two poten-

tials VX,(p) and VXa(p) to be AVXa ')

+ + + 4 + 1/3 +t 1/3
AV~a(p,r) = VNe(r) + AVC (prr) + 3 Ca(p(r) -p (r). (4.10)

The first term AVNe(r) is defined in Eq. (3.14) and the last
term is self-evident. However the second term needs further

elucidation. Beginning with the potential VC(p,r) we may

split it into a muffin-tin and a non-muffin-tin part--

VC(p~r)= VC (p)r + AVC,~ ) (4.11)

Separating the density into p and Ap in each part--

+ + + -+
VC(pI r) VC(p~r) + VC(Ap r) + VC (p,r) +VC (3pr).(4.12)

Substracting VC(p,r) we get--

VC(p r) VC (p r) VC(pr) + VC (p r) +VC(Aplr).(4.13)

Hence, we rewrite Eq. (4.10) as--

CVXa (P,r) = AVNe(r) + VC (Ap,r) + AVC (P,r) + bV(AP,r)

+ (~)13-p r)1/) (4.14)


We may consider nVXa(p,r) to be the perturbation to VXa(p,r)

in Eq. (3.4 4) Then the first order perturbation correction

Aci to the eigenvalue of Eq. (3.44)which we call E. is--

Ac = U.r)V p~)U. (r)dr (4.15)

or since no derivatives are involved

AE. =Ii p.r)CIVa (p r)dr.

Substituting Eq. (4.14) we have--

AC A ()AVe~~d + ()V(Apir)dr (4.16)

+ Api ()V(prr)d r +i p(r)V(Ap ,r)dr

+4 + 1/3 1/3 -+
+ p (r) Co (p (r) (r) ) dr.

We have made use of the properties of muffin-tin and non-

muffin-tin functions in reducing Eq. (4.14) to Eq. (4.16).

The last term and first term may not be simplified further.

The integrand of the second term involves the func-

tion VC(ap,r) and since this is an inconvenient function to

evaluate we rearrange this term as follows:

_~_1_1~__ ___


SPit r:)V (CpI r dr = p (r )VC( pr dr (4.17)

Ii2i (r )ap(r') dr'dI?


+ -f +

I '- '*= Ap(r)iVC 1,r:)dr
Pi(r)VC (Ap,r)dr = A~)V ~~r

The last function IVC~pi,r) is easier to calculate than the

original function VC(ap,r). Hence we rewrite Eq. (4.16) as--

Ae. =Ap.()av (rdr + ap(r)av, (p ,rdr (4.18)

+ Api (r) ~VC (plr)dr + Apyqr)6nVC(Apir)dr

4 +t + 1/3 +1/
3 a I ir P

It should be observed that this correction ae. to

the eigenvalue .i has a very similar form to the correction

6 to the muffin-tin energy and it can be easily
calculated in the same way. The effect of such first order

corrections to the eigenvalues is expected to be significant

for molecules such as benzene (C6H ) which have relatively

large intersphere regions. Since the present work is concerned

~____~_____~____I~ __ __ ~__ __

with the effect of the non-muffin-tin correction to the

muffin-tin energy no calculations to date have been per-

formed for these eig~envalue corrections.

4.4 The SCF Cycle

Self-consistency is implied in the muffin-tin Xa

equations Eqs.(3.44) and in the exact Xa equations Eqs.(2.36)

since in both the XaL potential--VXa(p) or VXa(p)--depends on

the orbitals through respectively the muffin-tin or exact

Xcr electron density. In Figure 4.1 we give a conceptual

flow chart for the principal steps in the SCF cycle for both

the exact and the muffin-tin Xa calculation. The dashed

lines in the diagram represent the steps taken in an exact

Xa calculation and the solid lines correspond to the con-

ceptual steps taken in a muffin-tin Xcc calculation.

To begin the iteration cycle to reach self-consistency

it is the usual practice to assume for the molecular Xa po-

tential a superposition of calculated Xal atomic potentials.

It is not necessary to use such potentials but they prove to

be a good starting point. From these superposed potentials

a muffin-tin potential is formed and from this potential via

the multiple-scattering (MS) equations the eigenvalues and

A and C expansion coefficients (see Sec. (3.5) and Appendix

B) are determined. From the A and C coefficients we form

orbitals {Ui }' orUi }which it should be emphasized are con-
tinuous and have continuous first derivatives throughout all

space--they are not themselves muffin-tin functions and hence

Figure 4.1. The dotted lines in this figure represent the
flow path for an exact Xal calculation and the
solid lines represent the conceptual steps
taken in a MT Xa calculation.

_ _

Figure 4.1. Flow Diagram for Conceptual Exact Xa Calcula-
tion and a Conceptual FIT Xa Calculation.

neither is the density formed from these orbitals. For the

exact Xa case having formed the density we then form the

new potential VXal(p) and the energy. After this we
test the convergence criterion which is usually related to

the new potential or density and if the convergence criterion

is met we are finished. If not then we return to solve the

MS equations again and continue iterating in this manner
until the convergence criterion is met. The muffin-tin

case is similar except that as indicated we have two ad-

ditional steps--of applying the muffin-tin average to the

density and to the potential. As indicated in Figure 4.1

we form the energy from the muffin-tin density p.
In actual MSXa( calculations orbitals, density and the

"un-muffin-tinned" potential are never formed. Instead the

muffin-tin density is formed directly from the A and C coef-

ficients and the muffin-tin potential VXa(p) is formed directly
from (See Appendices C and D). This calculational cycle

is shown in Figure 4.2 by the solid lines. In the same figure

the dashed lines indicate the steps performed in the calcula-

tion of A, the non-muffin-tin first order correction to

the functional . In this branch of the diagram we form

p and ap from the A anc C coefficients and from these we

calculate A and A

Figure 4.2. The dotted lines in this figure represent
the flow path for the calculation of A
and the solid lines represent the flow pX~
for an actual MT Xa calculation.

I_ _

Figure 4.2. Flow Diagram for an Actual M~T Xa Calculation
and the M~ajor Steps in the NMT Calculation.



5.1 Introduction

The NMT correction A to the MT Xa energy
as derived in Chapter III (Eq. (3.24)) and discussed
in Chapter IV is in the form of three-and six-dimensional

integrals. These integrals must be evaluated numerically

and there are two basic types of methods which might be

used: the Monte Carlo method or the Haselgrove and Conroy

methods. Although the Haselgrove and Conroy methods of

numerical integration have enjoyed some degree of success in

applications to quantum theoretical problems (see Haselgrove

(1961), Conroy (1967), and Ellis and Painter (1970)) they are

not adapted to integrating discontinuous functions as are

the integrals of A Hence we have used the Monte

Carlo method in these calculations. For a discussion of the

theoretical basis for the Monte Carlo evaluation of definite

integrals and statistical error analysis see Appendix E

and references contained therein.

For a definite integral F of a one-dimensional function


F = f (x) dx (5.1)

University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - - mvs