NONMUFFINTIN CORRECTIONS TO THE TOTAL
MOLECULAR ENERGY IN THE MULTIPLESCATTERING Xa
METHOD OF MOLECULAR CALCULATION
By
John Bryan Danese
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
1973
1III111IIIII1II II1I1I1UII1III11111111111 1111111111111111111111
ACKNOWLEDGEMENTS
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
acknowledged.
iii
TABLE OF CONTENTS
Page
. . ii
*.. V1
. ... vii
. .i
.. .. .1
.. .. 1
.. .. 3
.. .. 6
. . .
. .. .9
. .. .0
. . .15
ACKNOWLEDGEMENTS. ........
LIST OF TABLES* ********
LIST OF FIGURES ................
ABSTRACT* ........
CHAPTER I
INTRODUCTION. ..........
1.1 Molecular Methods ..........
1.2 The MuffinTin (MT) Approximation ..
1.3 The NonMuffinTin (NMT) Correction
to the MuffinTin (MT) Xce Total
Energy. . . . . . . .
THE MOLECULAR Xa TOTAL ENERGY
AND ONEELECTRON EQUATIONS. ...
CHAPTER II
2.1 Introduction. ............
2.2 Xa Theory ..............
2.3 Derivation of the Xa OneElectron
Equations .. .... .. ... .
CHAPTER III
THE MT APPROXIMATION TO THE Xa
TOTAL ENERGY AND ONEELECTRON
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 OneElectron Orbitals and
the MultipleScattering (MS)
Formal ism . . . . . . .
3.6 Explicit Expressions for the Components
of the NMT Correction a.....
. . .26
. . .26
. . .31
. . .38
. . .46
..49
. .
CHAPTER IV
THE RELATIONSHIP BETWEEN THE MT Xci
TOTAL ENERGY AND THE Xar TOTAL
ENERGY .... ....
. . 8
BIBLIOGRAPHY. ...............
BIOGRAPHICAL SKETCH ............
Page
4.1 Introduction. ...............
4.2 The Explicit Relationship Between
and . . . . . . .
4.3 Firs aOrder Co~rections to the Eigenvalues.
4.4 The SCF Cycle ...............
IER V NUMERICAL INTEGRATION AND IMPORTANCE
SAMPLING. . . . . .. . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
..68
..69
..78
..82
..88
..88
..91
..97
. 109
. 114
. 116
. 116
. 122
. 124
. 151
.172
. 179
. 185
. 192
. 199
CHAPT
5.1
5.2
5.3
5.4
5.5
Introduction. ..............
Importance Sampling in OneDimension. ..
Importance Sampling in nDimensions ...
Some Specific Sampling Functions. ....
A Note on the Haselgrove and Conroy
Methods of Numerical Integration. ...
CHAPTER VI RESULTS AND CONCLUSIONS .......
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. ..
APPENDIX A FUNCTIONAL DERIVATIVES. .......
APPENDIX B THE MULTIPLE'SCATTERING FORMALISM ..
APPENDIX C THE NMT DENSITY ...........
APPENDIX D EXPLICIT EXPRESSIONS FOR AaI
aV IN,IN and AVoIN c.....
APPENDIX E THE STATISTICAL EVALUATION OF DEFINITE
INTEGRALS .. .. .. .. . .. .
APPENDIX F THE 1/rl2 SINGULARITY ........
APPENDIX G THE BOND SAMPLING FUNCTION. .....
. . .
. .
.. 216
..224
.. 229
.. 235
. .240
LIST OF TABLES
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
LIST OF FIGURES
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 XZ 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 XZ Quadrant of
the C2 Molecule. 133
6.6 Contour Plot of AICL(r) in the XZ Quadrant
of the C2 Molecule. 138
6.7 Line Plot of AICL(;) for the C2 Molecule. 139
6.8 Schematic Division of the XZ Quadrant of the
C2 Molecule into Sampling Regions. 140
6.9 Contour Plot of alX(r) in the XZ Quadrant
of the C2 Molecule. 143
Vll
Figure
Page
6.10 Line Plot of AI ( )' for the C2 Molecule. 144
6.11 Contour Plot of AI (r) in the XZ 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 AtomicInner Vector Relationships Used
in the Derivation of nV a,IN and aV IN,IN, 201
D.2 AtomicAtomic Vector Relationships Used
in the Derivation of AV a,IN and AV IN,IN. 204
c c
D.3 AtomicOuter 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
viii
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
NONMUFFINTIN CORRECTIONS TO THE TOTAL
MOLECULAR ENERGY IN THE MULTIPLESCATTERING Xcl
METHOD OF MOLECULAR CALCULATION
By
John Bryan Danese
August, 1973
Chairman: J.C. Slater
Cochairman: John R. Sabin
Major Department: Physics
An investigation of the muffintin (]MT) approximation
to the total electron density in the multiplescattering
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 nonmuffintin (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 sixdimensional
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 moleculesC2 and
N2which it is argued should provide definitive tests for
the method of improving the MT XcL total energy with the
NMT correction. The multiplescattering 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 .
CHAPTER I
INTRODUCTION
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 everincreasing
size, weight, and complexity.
There are several categories or methods of molecular
calculation which can be broadly classified as ab initio,
semiempirical, and multiplescattering 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 "muffintin"
(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 LCAOMOSCF
(linear combination of atomic orbitalsmolecular orbitalself
consistent field) HartreeFock 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 semiquantitative.
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 LCAOMYOSCF HartreeFock 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 "muffintin" (MT) approximation to the
total electron density and potential, and the multiple
scattering expansion of the molecular orbitals. The APW
(augmentedpla newave) method and the KKR (KohnKorringa
Rostoker) methodsof solid state (see Slater (1963b)) employ
an exactly analogous muffintin approximation. Hence similar
conclusions may be drawn for these methods from the study
of the muffintin 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 multiplescattering
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 MuffinTin (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 oneelectron Xa( potential
and the total electron density are spherically averaged.
Between the atomic spheres and within the outer sphere in
the socalled 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 formulaa
functional of p which we will designate produces the
MT Xa total energy formulaa functional of the MT density
p which we will call . All of the integrals in
may be reduced to onedimensional 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 oneelectron 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 selfconsistent solutions of the
MT Xa oneelectron 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 nonmuffintin (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 semiquantitative
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
good.
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 nonmuffintin 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
tion.
__
1.3 The NonMuffinTin (NMT) Correction to the
M~uffinTin (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 partsthe 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
tion6 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 nonexistent (unbound) would be a good test of
this method.
The NMT correction A takes the form of three
and sixdimensional integrals which must be evaluated numeri
cally. In Chapter V we develop an importance sampling
scheme which combined with the M~onteCarlo 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 moleculesC2 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
itself.
CHAPTER II
THE MOLECULAR Xa TOTAL ENERGY AND ONEELECTRON EQUATIONS
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 muffintin 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 muffintin 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 oneelectron selfconsistent 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 oneelectron dif
ferential equations which a set of orbitals must satisfy in
order to make a minimum.
Xc9
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 oneelectron selfconsistent 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 spinpolarized 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 oneelectron 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 HartreeFock 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( rl1RaI
1 ii2p(rl)p(r2 1 + 1 +t + +
+ drldr2 + pkrl)U a(rlr
'12
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 nuclearelectron 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
" exchangecorre lation terms (see Slater (1972a)) The
last term is the classical nuclear interaction term. The
exchangecorrelation 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 nonspin 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
12
equilibrium position in the ground state exists in a closed
shell configuration, the results of spinpolarized and non
spinpolarized 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 nonspinpolarized assumption does not
preclude the possibility of fractional occupation numbers
or nonclosed 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,
and
UXaf = UXaf = UXcL (2.4)
where
+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 rlRai
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 nonspinpolarized assumption and the resulting slightly
simpler above form for will be used in the rest of
this study.
The exchangecorrelation 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
twothirds 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 HartreeFock 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
molecule.
Although there are similarities the Xa method is dif
ferent from HartreeFock 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 HartreeFock 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 nonintegral
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 OneElectron 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 nonspin
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(rr1 )*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),
6Ut(r)
= VTp~rl) + Ut(r)dr'dr+(Complex Conjugate)
II6p(r') VTUt~rr)
1 (2.14)
1
Using expression (2.13) above
6V p( ) N 2Z 6p(? )
T a .i
rl2r
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 rRa I rr2
+ 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.
(2.20)
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 (22 + 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
1/2
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
Irr 21
[V + +
a=1 I rR 
= ~Z C n1/2U.r.
k j ikhk~ kj(r
(2.28)
We may write the right hand side of therms of a matrix ele
ment
C~eik1/2 ~~~1/2
k j ikk
(2.29)
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 .
(2.33)
___~_ __~__________1^_1
1/2
Inserting this in Eq. (2.28) and cancelling n. em
N2Z 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 rRoi  jrr2 I
(2.35)
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 oneelectron self
consistentfield (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
treeFock equations except that in HartreeFock 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 socalled transitionn state" in which
the orbital occupation number is reduced by onehalf. (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.
(1972)).
Second, the solutions to the Eqs.(2.36) are selfcon
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 threedimen
sional functions which are not necessarily separable in a
spherical or any other coordinate system. The only
_~~ ~__ __~_I__~_ __
25
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 threedimensions for a molecule or solid. It
may approach a muffintin form but in reality of course
it does not have this simple form.
__ I
CHAPTER III
THE MT APPROXIMIATION TO THE Xa TOTAL ENERGY
AND ONEELECTRON EQUATIONS
3.1 Introduction
In this chapter we will discuss the muffintin (MT)
approximation to the Xa total energy for a molecule and
the MT oneelectron 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 nonmuffintin (NMT), a, part. From
the energy functional we will derive the oneelec
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),
27
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
oneelectron 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 threedimensional
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 nonmuffintin (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)
IN IN
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 atomiclike 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
region:
Afnrc~dr = Ifr): f(r)]dr
14r f Ir')sin9'd6'dg' sinededfr2dr.
(3.5)
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 atomiclike 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 nonspin polarized case)
+ 2 + + ac
+ 1 II2p(rl)p(r2 +t +
*drl + drlr
'12
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)
Nuclearelectron, Ee
Net
E = pr )d ; (3.9b)
ENel +C +t 12 tP~l 1
a=1 1 a
Coulomb, Ec
2p(rl)p(r2)
E j dr dr ; (3.9c)
c 2 1 2
'12
Exchange, Ex
1 +P +f +
Ex 2 prl)UXa rl drl ; (3.9d)
NuclearNuclear, 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 nuclearnuclear 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
(3.12)
We may define a nuclearelectron 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__~ __ ___ _
34
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
to
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
rl2
1 ( p 2P(rl) (r2 + + 2prlp(r
=2 dr2drl +l 2prlPr dr2drl
rl2 rl2
1 t 2A Apr
+ III2pr)Pr dr2drl (3.16)
'12
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 onehalf 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
rl2
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)
rl2
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
rrdr
'12
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 atomiclike
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 onedimensional integrals. Thus the first term in
Eq. (3.20) is in effect a twodimensional integral, the
second term in effect a fourdimensional integral, and the
last term is a sixdimensional 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)
where
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 +;
(3.22)
where the first term is that retained in the M~T approximation
to
~ 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)
where
~ ~ + 2 +i j fI1 +2 t
= n U (rl lU (rl)dr prVNe(rl drl
(3.23)
1 Ij2bp(rl1:)np(r2) t 4/3 4/3 +
+ 2 dr 2drl + Ca i(p(rl W(rl )drl
12
(3.24)
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
and
b = Ap~r )nVNe (rl)drl + Ap(rl)Abc(p, rl~drl
3.4 The Derivation of the MT Xcr OneElectron 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 oneelectron 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 nuclearnuclear term does not contribute
anything. Hence
6~ = n. Ur)[ V U.(r )]dr + 6VT (3.25)
Xa 1 1
where
V~~~~~~ = ~ 7er)r (rl Vc(pIrl)dr l
 I I
+ Ca p(r ) dr a,.. (3.26)
4/3
As indicated previously the functions p, V~e, Vc, and p are
radial functions in the atomiclike 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 atomiclike regions and the
volume integration in the intersphere region:
S1 + 4/3 2
ar=0
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
a=0
40
+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)
(3.29)
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(rlr)Nelrdl
(IN Ne,IN) =Ne,IN (
.30)
.31)
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 )
12
= 6(r r')dr drl
'12
(3.33)
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
1
6p(r')
= ~ ~ ~  6( 'V(pr)r.
(3.34)
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
regions
6 1 jl 2p(rl)p(r2) + + B P''1)frao
2 dr2drl1
6p(r') rl2
1 V NVcI for the
intersphere
region
(3.35)
The last term of 9T to consider is the exchange term. For
an atomic region it is
6j p~4/3r 2 41/3 2
4xC 6r) rdr x C rl) 6rr)lr
4 1/3
=4r C p(r')
3 a
(3.36)
and for the intersphere region it is
6 4/3 4 1/3
CUV p= V Cp .
aININ IN 3 a IN
6p
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') /]
ar=0
6p(r') 2 +t
x ~, rt' dr'dr
6U.(r)
~* ~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 atomiclike
and in the intersphere regions. In an atomiclike 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 interatomic region:
__
IN
6U.(2)
1 n.U C')U.(f')d?'
V
IN IN i
6
~*i r
1 ~ i(~
V
IN
1 1,
V
(3.39)
Substituting the expressions ofEqs*(3.38) and (3.39) in
Eq. (3.37) we have
 N
6V = 1
a=0
~* 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 .
(3.40)
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.
(3.41)
In this expression it is understood that p, VNe, Vc, and
1/3
p are radial functions in the atomiclike 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 muffintin 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 muffintin 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)),
a
Xa
E. = (3.45)
an.
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 OneElectron Orbitals and
the MultipleScattering (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~ellknow: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 multiplescattering
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 multiplescattering 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)
1,m
where the Cl are expansion coefficients to be determined,
Rm
the Y~m are real spherical harmonics, and the atomiclike
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
(3.47)
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
(3.49)
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 muffintin approximation
to the density, and the multiple scattering formalism as
outlined above together constitute the multiplescattering
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
by
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 multiplescattering 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 termsAV~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 moleculeatomic, outer, and
intersphereis a tediomibut straightforward 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 nuclearelectron terms and these need no derivation.
The other terms of nV (p, r) have no counterpart in the
nuclearelectron terms and the derivation of these terms is
relegated to Appendix D.
The NMT Nuclearelectron TermsaV (e
The explicit definition of VNe( ) is given, once again
as
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 muffintin
average in an atomic or the outer sphere region is an angular
averagean average performed in the region of interest.
The muffintin 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
Ne
electron in the ath atomic region and the
Bth nucleus located in the 8th atomic
region.
Vo,a1
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 muffintin average for each one of these
potentials.
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 
2Z
a (3.50)
ra
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)
'Ne
For Va, we have
Ne
2Z,
V",58() = ,(3 53
Ne + + .5
r RI(
53
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 
by
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
(3.56)
(3.57)
Hence
and from Eq. (3.55) and Eq. (3.56)
2Z
'Ne
Ru
(3.58)
IYem(r)dG = 4n6(R)6(m)Y .
a, r +B 1 t dG ,
'Ne 
4x Ir Rcf 
I IO
(a)
( b)
Figure 3.2.
Vector Relationships Used for the Derivation
of (a) V ,and (b) VO 4
INe ~e'
I _
Ke~nce
Br 2Zt 2ZB
ide **
Icr Ra Rcr
(3.59)
Then,within an atomic sphere a the total potential AVaT(
Ne
is given by
N
dV~e(r) =t 2Z( ,1 _
B=1f a R g
1 ) .
Ra
(3.60)
For Vroea(r) the point r is in the outer sphere
region
'Ne + + 
r R
O C1O
(3.61)
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 .
(3.62)
+ +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
I
2o4a I 1t dGt (3.64)
or
2Z
vyo(.a a (3.65)
~Ne
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 threedimensional
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
2Z
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
ro
for rI ro>l jRao
(3.68)
and
1 4i
Sro Relo toO
(3.69)
Hence, the volume integration of VIa inside the outer sphere
is
bo
rI2dr +4nj r dro }
R
CLo
R
ao
2Z di 2Z { 4r (
r r R R
o toL ao 0
(3.70)
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
theother type case where we integrate over; the other (N1)
atomic spheres B. Within the oth sphere r R=r and
4nb2
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)
Then
1 + 41T 2 g
2Z. II dr =2Z rd
R aB1 aB 0
=2Z (3.73)
S3R
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
ac=1
IN +t +
Then aVNe(rlN) at any point rI in the intersphere region
is given by
N 2Z
IN + e IN
AV (r)=1V ( (3.76)
N ea=1 IIN acl
The NMT Coulomb TermsaV (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 moleculeatomic, 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 nuclearelectron 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
B=1
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
Cr,
in the other (N1) f3 atomic spheres Va,
in the outer sphere Va'o
cI
in the intersphere region VN
Using the same kind of notation and definitions we may
divide the potential VI in the intersphere region into
parts
VIN ~ IN,a + IN,o + IN,IN (.9
cr=1
Similarly for Vo the potential in the outer sphere region
Vo= V~ +V~ VN.(3.80)
c c c c
al=1
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
two
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
2Q~
N
B=1
8=1
N
B=1
AVa ( o) = T
8=1
20 ( 1 1 )
r Ra  Rol
(3. 82b)
0 (r')
I ouT dr'
r'
b
o
VC1o(r) 2
VfCL o (r ) =C Vo (o
CD C
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
3
b
8 1
Ru
N
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
(3.82d)
N N
a=1 al=1
2Q,
r R
IN oa
2
N N b
INa+ 8n e( o
INcl IN
V 2
al=1 IN a=1
2
R
oa
6
2
b
a
2
3
b
Rcr
N
B=1
2
b
2QI 1 471[ o
r I Rea VI 2
2
R
oa )
6
2
b
ar
2
N N
A V I~c(rI )= C
a=1 a=1
3
b0
Ra
N
1
B=1
(3.82e)
CO
a (r')
ouT dr
r'
cN IN 
IN'o +C
IN, o
V c r )
c IN ~
c ININ
(3.82f)
3
b
a
It )
r 
IN oai
N
2
a=1
5
b
o
5
IN IN _
c C IN 4.rIN C o rIN
"gl"' (r N) = 4v~ (b2 4w
3IN
2
b
a
5
2
R
oa1
3
N
Sb3
a
a=1
3
b
B1)
Ra
(3.82g)
2
b
3 o
b [
a 2
2
b
a
2
2
R
oa
6
N
1 v
B=1
8xr N
 I
3V
ILN ar=3
 V~MIN,IN~r)
AV (INr ) =VININ~r N)
N N
~ cO', o L
a=1 cr=1
N  N 2
Vc r o)= Z
a=1 a=1
203
ii
0 C
Q,
N N
A~~tr)= Q( 1 1 )
r R ro
cr=1 ar=1 o oa
(3.82h)
~o b
a (r')
ouT
dr'
r'
r
o
a (r')dr' + 2
ouT
o,o + o0o
Ve (ro)= Vc' (ro)
avPOI~o)o= )=0
(3.82i)
4vrb N
VOIN~ + o
7 (r)= 2p 
3r
o ac=1
4nb N
V,~ (r )= 2pIN
3r
o a=1
3r R
4rb 2pV
a IN IN
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
1~
___
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 nuclearelectron 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 nuclearelectron and
Coulomb terms).
CHAPTER IV
THE RELATIONSHIP BETWEEN THE MT Xa TOTAL ENERGY
AND THE Xa TOTAL ENERGY
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 muffintin 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 muffintin
approximation to the total density in the total energy eval
uation. This is accomplished by using the set {U } to
evaluate the nonmuffintin correction A to the muffin
tin total energy . We shall show that this approximation
to is analagous to performing RayleighSchroedinger
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 muffintin 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
fintin approximations are performed, and how the approxi
mate method for calculating fits into the present
muffintin scheme for calculating
Only the development of an efficient multidimensional
multicenter numerical sampling and integration scheme has
made possible the practical evaluation of the nonmuffin
tin correction A since it involves up to sixdimen
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 muffintin and nonmuffintin parts
= + a
(3.11)
where
1 + + + +4/
+~~~~~~ p~lV(,r)r a rl) drl
2Z Z
+ r (3.23)
and
n = bp(r ) nVNE(r )dr + bp(rll) AVclp, rl~drl
1 II2ap(rl ) p(r2 + +4/
+ 2dr2drl C. I[p(rl~4
'12
+4/3
p(r ) ]drl (3.24)
The definition of ap has given as the difference between
the density formed from the orbitals {U.} and the muffintin
density (see Eqs. (3.1) and (3.3)) formed from the same
orbitals:
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)
where
+ a 2 4 + 1/3
V (p, r) =C +2 dr + C pr)
xa ~ ~ ** r * .* 3 a
a=1 2
(2.35)
In Chapter III we showed that the set of orbitals {U } which
make a minimum are solutions to
xa
2 t ++ +
[V +V (p, r) ]U. (r = ..(),(.4
Xa 1u 1 r) 134
where
 + + 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 sixdimensional
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_ muffintin
multiplescattering 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 casesC2
being an examplethe 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 nonmuffintin
correction a to the muffintin 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 nonmuffintin correction a and
xa
therefore indicate whether the exclusion of this term is a
major source of error in current MSXa calculations.
This approximationusing the {U.}to calculate
73
is quite analogous to first order RayleighSchroedinger
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)
and
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
given
C = cO + El (4.5)
and
E = .
(4.6)
_~__ ___ ___~___ _II__ _~^
In analogy, in the present case we have instead of an
expectation value a functional which we may separate
xar
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 :
xa
6 = 0 (4.8)
Using the set {U.} we may calculate n which we see to be
analogous to the first order correction in RayleighSchroedinger
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 muffintin 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 orbitalsa
wellknown 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
_ _
75
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 HellmanFeynman 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 singleexcitation 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 oneelectron 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 selfconsistently obtained minimize the energy
functional defined by this oneelectron 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
functional.)
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 selfconsistently 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 oneelectron 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
offdiagonal 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 muffintin.approxi
mation gives a rather bad result for the binding energy.
Muffintin 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 muffintin energy of C2 to three figures is
150 Rydbergs). We estimate then that the muffintin 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 muffintin 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 muffintin 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 nonmuffintin 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 Ne2a molecule qualitatively different from
C2has corroborated this trend for the nonmuffintin
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 nonmuffintin correction A to
the muffintin 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 muffintin single
electron Xal equations. (This of course may be done for
transition state eigenvalues as well as any other con
figuration with integral or nonintegral occupation numbers.)
If we compare the exact Xar singleelectron equations
of Eq. (2.36) with the muffintin singleelectron 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 selfevident. However the second term needs further
elucidation. Beginning with the potential VC(p,r) we may
split it into a muffintin and a nonmuffintin 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)
+ (~)13p r)1/) (4.14)
80
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 muffintin and non
muffintin 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~__ ___
81
SPit r:)V (CpI r dr = p (r )VC( pr dr (4.17)
Ii2i (r )ap(r') dr'dI?
rrr'
+ 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 muffintin 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 nonmuffintin correction to the
muffintin energy no calculations to date have been per
formed for these eig~envalue corrections.
4.4 The SCF Cycle
Selfconsistency is implied in the muffintin Xa
equations Eqs.(3.44) and in the exact Xa equations Eqs.(2.36)
since in both the XaL potentialVXa(p) or VXa(p)depends on
the orbitals through respectively the muffintin 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 muffintin 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 muffintin Xcc calculation.
To begin the iteration cycle to reach selfconsistency
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 muffintin potential is formed and from this potential via
the multiplescattering (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
spacethey are not themselves muffintin 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 muffintin
case is similar except that as indicated we have two ad
ditional stepsof applying the muffintin average to the
density and to the potential. As indicated in Figure 4.1
we form the energy from the muffintin density p.
In actual MSXa( calculations orbitals, density and the
"unmuffintinned" potential are never formed. Instead the
muffintin density is formed directly from the A and C coef
ficients and the muffintin 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 nonmuffintin 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
form
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.
CHAPTER V
NUMERICAL INTEGRATION AND IMPORTANCE SAMPLING
5.1 Introduction
The NMT correction A to the MT Xa energy
xar
as derived in Chapter III (Eq. (3.24)) and discussed
xa
in Chapter IV is in the form of threeand sixdimensional
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 onedimensional function
f(x)
F = f (x) dx (5.1)
