 Title Page 
 Acknowledgement 
 Table of Contents 
 List of Tables 
 List of Figures 
 Abstract 
 Introduction 
 Exact Xa method and muffintin... 
 Muffintin charge density... 
 HellmannFeynman force calcula... 
 Dipole moment calculations 
 Conclusions 
 Appendix A: Derivation of the muffintin... 
 Appendix B: The multiple scattering... 
 Appendix C: HellmannFeynman theorem... 
 Appendix D: Proof of the independence... 
 Appendix E: Expansion theorems... 
 Appendix F 
 Appendix G: Evaluation of the sine... 
 Bibliography 
 Biographical sketch 

Full Citation 
Material Information 

Title: 
HellmannFeynman forces and dipole moments using multiplescattering Xa wavefunctions for diatomic molecules / 

Alternate Title: 
Multiplescattering Xa wavefunctions for diatomic molecules, HellmannFeynman forces and dipole moments using 

Physical Description: 
ix, 156 leaves : ill. ; 28 cm. 

Language: 
English 

Creator: 
Li, Choy Hing, 1946 

Publication Date: 
1975 

Copyright Date: 
1975 
Subjects 

Subject: 
Dipole moments ( lcsh ) Molecules ( lcsh ) Physics thesis Ph. D Dissertations, Academic  Physics  UF 

Genre: 
bibliography ( marcgt ) nonfiction ( marcgt ) 
Notes 

Thesis: 
ThesisUniversity of Florida. 

Bibliography: 
Bibliography: leaves 153155. 

General Note: 
Typescript. 

General Note: 
Vita. 

Statement of Responsibility: 
by Choy Hing Li. 
Record Information 

Bibliographic ID: 
UF00098149 

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  000171233 oclc  02954297 notis  AAT7658 

Downloads 

Table of Contents 
Title Page
Page i
Acknowledgement
Page ii
Table of Contents
Page iii
Page iv
List of Tables
Page v
List of Figures
Page vi
Page vii
Abstract
Page viii
Page ix
Introduction
Page 1
Page 2
Page 3
Exact Xa method and muffintin approximate Xa method
Page 4
Page 5
Page 6
Page 7
Page 8
Page 9
Page 10
Page 11
Page 12
Page 13
Page 14
Page 15
Page 16
Page 17
Page 18
Muffintin charge density approximation
Page 19
Page 20
Page 21
Page 22
Page 23
Page 24
Page 25
Page 26
Page 27
Page 28
Page 29
Page 30
Page 31
Page 32
Page 33
Page 34
Page 35
HellmannFeynman force calculations
Page 36
Page 37
Page 38
Page 39
Page 40
Page 41
Page 42
Page 43
Page 44
Page 45
Page 46
Page 47
Page 48
Page 49
Page 50
Page 51
Page 52
Page 53
Page 54
Page 55
Page 56
Page 57
Page 58
Page 59
Page 60
Page 61
Page 62
Page 63
Page 64
Page 65
Page 66
Page 67
Page 68
Page 69
Page 70
Page 71
Page 72
Page 73
Page 74
Dipole moment calculations
Page 75
Page 76
Page 77
Page 78
Page 79
Page 80
Page 81
Page 82
Page 83
Page 84
Page 85
Page 86
Page 87
Page 88
Page 89
Page 90
Page 91
Page 92
Page 93
Page 94
Page 95
Page 96
Page 97
Page 98
Page 99
Page 100
Page 101
Page 102
Page 103
Page 104
Page 105
Page 106
Page 107
Page 108
Page 109
Conclusions
Page 110
Page 111
Page 112
Page 113
Page 114
Page 115
Page 116
Page 117
Page 118
Page 119
Appendix A: Derivation of the muffintin Xa oneelectron equations
Page 120
Page 121
Page 122
Page 123
Page 124
Page 125
Page 126
Page 127
Page 128
Page 129
Appendix B: The multiple scattering method
Page 130
Page 131
Page 132
Page 133
Page 134
Appendix C: HellmannFeynman theorem and the Xa wavefunctions
Page 135
Page 136
Page 137
Page 138
Appendix D: Proof of the independence of dipole moment on the choice of the origin for a neutral molecule
Page 139
Page 140
Page 141
Appendix E: Expansion theorems for the Bessel functions and the Neumann functions
Page 142
Page 143
Page 144
Page 145
Appendix F
Page 146
Page 147
Page 148
Appendix G: Evaluation of the sine and cosine integrals
Page 149
Page 150
Page 151
Page 152
Bibliography
Page 153
Page 154
Page 155
Biographical sketch
Page 156
Page 157
Page 158
Page 159

Full Text 
HELLMANNFEYNMAN FORCES
AND DIPOLE MOMENTS USING
MULTIPLESCATTERING Xa WAVEFUNCTIONS
FOR DIATOMIC MOLECULES
By
Choy Hing Li
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
1975
ACKNOWLEDGEMENTS
The author would like to express his gratitude to
Professor John W. D. Connolly who suggested the topic of
this study and would like to thank him for his guidance,
instruction and criticisms during this research.
The author would also like to acknowledge the other
members of his advisory committee for their suggestions
and criticisms. Also, thanks are due to all my other
colleagues of the Quantum Theory Project at the University
of Florida. The author particularly extends his apprecia
tion to Dr. B. Danese for his many valuable discussions,
continuous encouragement and especially for providing the
author with his nonmuffintin total energy program.
Finally, the author would like to thank his wife
for her patience and encouragement.
This research was supported in part by the grant from
the National Science Foundation which is also appreciated.
TABLE OF CONTENTS
Page
ACKNOWLEDGEMENTS..........................................ii
LIST OF TABLES. .............................................v
LIST OF FIGURES............................................ vi
ABSTRACT .................................................vii
CHAPTER I INTRODUCTION................................ 1
CHAPTER II EXACT Xa METHOD AND MUFFINTIN
APPROXIMATE Xa METHOD ......................4
2.1 Introduction....................................... 4
2.2 Exact Xa OneElectron Schr6dinger Equations.......4
2.3 MuffinTin Approximate Xa OneElectron
Equations and the Solutions....................... 9
2.4 Relation Between the Exact Xa Method and
tne Muffintin Approximate Xa Method..............13
CHAPTER III MUFFINTIN CHARGE DENSITY APPRUXIMATION....19
3.1 Introduction.....................................19
3.2 MuffinTin Xa Total Energy Calculations..........20
3.3 Force Calculations...............................21
3.4 Dipole Moment Calculations.......................30
CHAPTER IV HELLMANNFEYNMAN FORCE CALCULATIONS........36
4.1 Introduction. ................. ................. 36
4.2 Prolate spheroidal Coordinates...................36
4.3 HellmannFeynman Theorem..........................39
4.4 Method of HellmannFeynman Force Calculations....41
4.5 Results and Discussion ...........................62
CHAPTER V DIPOLE MOMENT CALCULATIONS .................75
5.1 Introduction......................................75
5.2 Method of Dipole Moment Calculations.............75
5.3 Results and Discussion.............................99
CHAPTER VI CONCLUSIONS ............................... 110
APPENDIX A
APPENDIX B
APPENDIX C
APPENDIX D
APPENDIX E
APPENDIX F
APPENDIX G
Page
DERIVATION OF THE MUFFINTIN
Xa ONEELECTRON EQUATIONS...............120
THE MULTIPLESCATTERING METHOD..........130
HELLMANNFEYNMAN THEOREM AND THE Xa
WAVEFUNCTIONS........................... 135
PROOF OF THE INDEPENDENCE OF DIPOLE
MOMENT ON THE CHOICE OF THE ORIGIN
FOR A NEUTRAL MOLECULE..................139
EXPANSION THEOREMS FOR THE BESSEL
FUNCTIONS AND THE NEUMANN FUNCTIONS.....142
EXPANSION THEOREM FOR Yo(' )/r2 ........ 146
1 A A
EVALUATION OF THE SINE AND COSINE
INTEGRALS. ...............................149
BIBLIOGRAPHY ............................................153
BIOGRAPHICAL SKETCH .................................... 156
LIST OF TABLES
Table Page
4.1 Force Calculations for H 1( : I 2) 64
2 g g
4.2 Comparison of Ftotal and 3
4.3 ZDependence of HellmannFeynman Force
for H2 at R=1.4 a.u. 68
4.4 iValues in the MSXa Calculation for N2(1) 70
4.5 HellmannFeynman Forces for N2 at R=2 a.u.
Calculated with MSXa and HartreeFock
Wavefunctions 71
5.1 Dipole Moment Components for LiH(1 +)
at R=2.75, 2.836, 3.015, 3.25 a.u. 101
5.2 Dipole Moments and Their Derivatives of
LiH at R=3.015 a.u. 103
5.3 6 and eDependence of Dipole Moment for
LiH at R=3.015 a.u. 104
5.4 Dipole Moment for BH(1E+) at R=2.336 a.u. 107
5.5 Dipole Moment for CH( 2) at R=2.124 a.u. 108
LIST OF FIGURES
Figures Page
2.1 Flow Chart of the Conceptual Exact
Xa Method 7
2.2 Flow Chart of the Conceptual Exact
Xa Method and the MSXa Method 15
3.1 Partitioning of the Space of a Homonuclear
Diatomic Molecule in the MSXa Method 23
4.1 Representation of a Point in the Atomic
Coordinates 37
4.2 Relation of rA, R and 8A 49
4.3 Boundary of the Intersphere Region in Prolate
Spheroidal Coordinates with A and B as the Foci 49
4.4 Force Curves for H2 67
5.1 Schematic Representation of A Heteronuclear
Diatomic Molecule in the MSXa Calculation 76
5.2 Boundary of the Intersphere Region in
Prolate Spheroidal Coordinates with A and O
as the Foci 83
5.3 Boundary of the Intersphere Region in
Prolate Spheroidal Coordinates with O
and B as the Foci 89
5.4 Boundary of the Intersphere Region in
Prolate Spheroidal Coordinates with A and B
as the Foci 96
5.5 Dipole Moment v.s. Separation for LiH ( 1+) 102
6.1 Schematic Representation of a Triatomic
Molecule in the MSXa Method 113
6.2 Rotation of Atomic Coordinates to Align with
the internuclear Axis 115
Figures Page
D.1 Transformation of Coordinates 140
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
HELLMANNFEYNMAN FORCES
AND DIPOLE MOMENTS USING
MULTIPLESCATTERING Xa WAVEFUNCTIONS
FOR DIATOMIC MOLECULES
By
Choy Hing Li
December, 1975
Chairman: John W. D. Connolly
Major Department: Physics
The HellmannFeynman forces for H2 and N2, and the
dipole moments for LiH, BH and CH are calculated with the
MultipleScattering Xa wavefunctions. First, the muffintin
charge density approximation is used to compute these two
quantities. The results show that the internuclear forces
for homonuclear diatomic molecules are always repulsive, and
the dipole moments for LiH, BH and CH are far from the correct
values. Thus, it is necessary to use the nonmuffintin charge
density in such molecular property calculations. The main
difficulty in such nonmuffintin calculations is to evaluate
some threedimensional integrals over the intersphere region.
To achieve this, we express the integrands in Prolate Spher
oidal Coordinates and reduce the integrals into some onedi
mensional integrals. The results show that both the force cal
culations and the dipole moment calculations improve a lot
over the use of the muffintin charge density approximation.
For H2, the HellmannFeynman force curve is obtained and com
pared with the curve for the derivative of the total energy,
viii
thus checking whether the HellmannFeynman Theorem is valid
for the MSXa wavefunctions. For N2, the force components
contributed by each of the electronic states are calculated
at the experimental equilibrium separation and are compared
with the HartreeFock results. For LiH, BH and CH, the dipole
moments thus obtained are respectively 7.29 D, 1.24 D and
1.958 D at the experimental equilibrium separations. Also,
the derivative of the dipole moment for LiH at the equili
brium separation is obtained to be 1.27 D/a.u. We conclude
that the nonmuffintin correction is important in the force
and dipole moment calculations. Any further improvement has
to come from the correction to the MSXa wavefunctions.
CHAPTER I
INTRODUCTION
In the investigation of the electronic structure
and properties of atoms, molecules and solids, we often
encounter the problem of treating the exchangecorrelation
effect properly. The Xa theory suggested by Slater provides
us a good way to handle such a problem (see Slater(1972a)).
Together with the muffintin approximation, the Xa theory
can be used to investigate the molecular system in a prac
tical and accurate way. Such a method, called the Multiple
Scattering Xa method or the MSXa method, was first suggested
by Slater(1965) and then was implemented by Johnson(1966).
Many of the results obtained with the MSXa method
show that it gives better results in the description of
the oneelectron features than do the ab initio SCFLCAO
methods, and is significantly much faster in calculation
(see Connolly and Johnson(1973)). However, in many cases,
we find that the total energy as a function of molecular
conformation is inaccurate. The energy curves of many diatom
ic molecules, e.g., show that the systems are unbound. It
was pointed out by Connolly (see Connolly and Johnson(1973))
that such a discrepancy may be mainly due to the muffintin
approximation, not the Xa approximation, in the calculation
of the total energy. The current version for calculating
the total energy assumes the muffintin approximation for
the charge density. The result may be improved by taking
the charge density, instead of its muffintin part, gererated
from the converged wavefunctions in the calculation. This
was verified by Danese on the calculations of carbon and
neon molecules (see Danese(1973)).
If the discrepancy in the total energy calculation
is indeed due to the use of the muffintin part of the
charge density and not the Xa approximation, and the result
does improve by using the complete form of the charge densi
ty, then it will mean that the converged MSXa wavefunctions
are not too far from the true Xa wavefunctions and it may
be useful to employ the MSXa wavefunctions in other molecular
property calculations. In this dissertation, we shall inves
tigate the quality of these MSXa wavefunctions directly in
the dipole moment and force calculations of some diatomic
molecules.
It was proved by Slater (see Slater(1972b)) that the
HellmannFeynman Theorem is satisfied by the exact Xa wave
functions. That is, the quantum mechanical force acting on
a nuclear coordinate is equal to the electrostatic force
exerted by the other nuclei and the electronic charge cloud
which is formed by the exact Xa wavefunctions. Thus, it is
useful to check the quality of the MSXa wavefunctions by
testing if the theorem is also valid for these approximated
Xa wavefunctions. Also, as it is well known that the dipole
moment is more sensitive to the wavefunction than is the
energy, it is useful to check the quality of the MSXa wave
functions by doing some dipole moment calculations.
The difficulty in all these calculations is that
we have to evaluate some threedimensional multicenter
integrals over an awkward region. In this dissertation, we
solve this difficulty by a computational technique which
reduces the integrals into onedimensional integrals whose
integrands involve sine and cosine integrals. The main
point in this method is to employ the Prolate Spheroidal
Coordinates. It is hopeful that this method can be used in
some other molecular property calculations. We shall dis
cuss it in detail in Chapter IV and Chapter V.
In Chapter II, we shall give a brief discussion of
the exact Xa method and the approximate MSXa method. In
Chapter III, we shall discuss the effect of assuming the
muffintin charge density in the calculation of total energy,
force and dipole moment, thus showing the necessity of using
the nonmuffin charge density in the calculations. Chapter
IV will be devoted on the force calculations and Chapter V
will be devoted on the dipole moment calculations. Finally,
we shall conclude in Chapter VI.
CHAPTER II
EXACT Xa METHOD AND MUFFINTIN APPROXIMATE Xc METHOD
2.1 Introduction
The discussions in this chapter shall mainly concern
the basic Xa theory and also its muffintin approximate form
in the application to molecular systems. First, we shall
briefly discuss the basic Xa theory and its exact one
electron Schr6dinger equations, in which no specific restrict
ion is placed on the form of the charge density and the
potential. Secondly, we shall discuss the MSXa method, in
which the charge density and the potential are muffintin
averaged. The approximate Xa oneelectron Schr6dinger equa
tions and their solutions are discussed. The quality of
the solutions of the approximate Xa oneelectron equations
is one of our main concerns in this dissertation.
2.2 Exact Xa OneElectron Schrddinger Equations
The Xa method was suggested by Slater and has been
discussed in detail by many others (see Slater(1972a),
Slater and Johnson(1972)).
The Xa oneelectron Schrodinger equations for the
spinorbitals ui(r) can be derived by the variational princi
ple from an expression of spin up and spin down charge den
sity:
p() = Z niui()ui(7 )
i*
p() = niui()ui(r)
iP
where ni are the occupation numbers and the summation is over
the spin up and spin down orbitals respectively. Also we de
fine the total charge density as
p(r) = p (r) + p (r)
Then the Xa total energy is
f 2 p 1 2p(rl)p(r2) 1d
+ ; P (rl)drI + r dr drl
I J + t J 12
+ P(rl)Ux (r)dr + o(r5)Uxarl)dr
+ 2 2ZpZq(2.1)
Rpq
Here the nuclear coordinates are Rp and their atomic
numbers are Zp. Rydbergs are used as units of energy, atomic
units as units of distance. The first term is the kinetic
energy. The second term is the nuclearelectron coulomb
energy. The third term is the electronelectron coulomb
interaction term. The fourth and the fifth are the "exchange
correlation" energy terms of either spin. The last one is
the nuclearnuclear interaction term.
The above expression is exact except for the exchange
correlation term..Here we define it as:
Ut (?) = 9a((3/47)p (r))1/3 2.2)
S =/3(2.2)
and similarly for U (r).
xa
The Xa total energy is a function of the occupation
numbers n. and a functional of the spinorbitals u.. It can
1 1
be shown that if we vary the ui to minimize the total energy
with fixed occupation numbers ni, we can get the Xa one
electron Schrodinger equations as:
2
{V1 + Vc(Prl) + Vxa(P,r )} ui(r) = Ei ui(rl) (2.3)
where N 2Z p(2
V(p,l) = d + 2 r (2.4)
P'l r1 12
is the electrostatic potential at position r1 and
1/3
F 6c((3/47)p,(rl))/ for spin up ui
Vxa (P,rl) =
[ 1/3
6a((3/47)p (rl)) for spin down ui
(2.5)
is the exchange correlation potential at rl.
Since the potential in the Xa oneelectron equations
is a functional of the spinorbitals ui, the Xa method is a
selfconsistent method. We start by assuming a potential
Vc+Vx. Solve the differential equations and find the spin
Fig. 2.1 Flow Chart of the Conceptual Exact Xa Method
orbitals. From these spinorbitals we form the charge densi
ty. Then, by equations (2.4) and (2.5), we calculate the
potential V,+V,, .which is then compared with the starting
potential Vc+Vx,, If they are "almost" equal to each other,
the calculation is said to have converged. Otherwise we start
with V%+V, and repeat the cycle again until convergence is
achieved. Then we obtain the spinorbitals ui, and with
equation (2.1) we calculate the total energy . The
conceptual computational scheme is shown in figure (2.1).
There are some difficulties in this scheme, namely:
(1) There is a parameter a in equation (2.2)..In
atomic systems, the a value can be chosen in several ways,
all of which amount to almost the same total energy (see
Lindgren and Schwarz(1972)). One of them is to set the a
value such that the total energy calculated be equal to
the exact HartreeFock total energy. Another is to set the a
value such that the HartreeFock total energy calculated
with the Xa orbitals is a minimum. Still another is to set
the a value such that the virial theorem is satisfied.
However, for molecular systems, the space may not be as
homogeneous as that of the atomic systems, and we need
some other prescriptions to determine the a value.
(2) In general, the potential is a very structured
function. For atomic systems, the potential is spherically
symmetric and it is possible to separate the solutions in
spherical coordinates. However, for molecular systems, the
differential equations usually cannot be separated in any
coordinates, unless some approximations are made in the form
of the potential.
(3) Suppose the differential equations are solved
and the spinorbitals are obtained, we still have to calcu
late the potential by equations (2.4) and (2.5). and to
evaluate the total energy by equation (2.1). All of these
involve three and six dimensional integrals with compli
cated integrands.
Due to all these difficulties, the MSXa method is
suggested. With the muffintin approximation and the multi
ple scattering formalism, the differential equations can
be solved easily and the expressions for the total energy
and the potential become much more simple. We shall discuss
these in the next section.
2.3 MuffinTin Approximate Xa OneElectron Equations
and the Solutions
The MultipleScattering Xa method (MSXa) was suggest
ed by Slater in 1965 and its implementation was begun by
Johnson in 1966. There are three basic characteristics in
the MSXa method: the exchangecorrelation approximation, the
muffintin approximation for the charge density and the poten
tial, and the multiple scattering formalism.
First, we divide the coordinate space of a molecule
into three regions: (I) the contiguous spherical regions
surrounding the different nuclei, (II) an intersphere region
which is outside of the spheres of region(I) but inside
region (III), which is the exterior of a sphere enclosing
all the spheres in region (I). The muffintin approximation
is such that we approximate the potential and the charge
density in region (I) and region (III) by a spherical aver
age, and those in region (II) by a volume average, namely:
(r) = +ff(t)dQ if r, I or III
f(S) =
= f(r)dr if rE II
VII II
where f is the potential or the charge density function,
V is the volume of the intersphere region, and f(r)
is the averaged quantity.
With the muffintin approximation, it can be shown
that the difficulties mentioned in the previous section
can be solved, namely: the differential equations for the
spinorbitals are separable, the potential can be computed
easily, and the calculation of the total energy reduces to
the evaluation of onedimensional integrals.
Starting with the muffintin charge density p(r),
we have an expression for the approximate Xa total energy:
r C i%* 2 ()d +U + UT
= 2 ni ui(r)V )u(r)dr + UN +
where ui are the spinorbitals that minimize
UNN is the nuclearnuclear interaction energy,
UT is the sum of the nuclearelectron coulomb energy,
electronelectron coulomb energy and the exchangecorrelation
energy:
UT4 2Zp ) (rl)drl
pzI I 1ipI
+ CaJ (+(r )4/3 + ( /)4/3d
1 r r + 2( (
+ 2drlp(rl)Jdd 212
2 J r12
9 1/3
Ca = 9a(3/4T)1/3
With this expression for the total energy, it is
shown in Appendix A that the necessary and sufficient
condition for be a minimum is that the hi satisfy
the following differential equations:
2 + r e n, %%
V + V(p,r) + Vex() ui(r) = iui(r) (2.6)
where V (p, ) and Vex(p,$) are the muffintin average of
V (p,r) and Vex(p,r) as defined in equations (2.4) and (2.5).
We want to emphasize that the muffintin approximate
oneelectron equations (2.6) are derived from the fact that
the charge density is of muffintin form. There is only
one approximation (besides the Xa exchangecorrelation approx
imation) involved in this derivation, and that is the muffin
tin charge density approximation.
To solve the muffintin approximate Xa oneelectron
equations:
{ 2 + Vxa(p,)} u (r) = Eiui()
where we write V (p,r)+Ve () as V (p,r), we may use the
partial wave expansions of multiple scattering theory in
the different regions of the molecule. Within an atomic
sphere S (including the outer sphere which we call the 0th
sphere), the potential is spherically symmetric. Thus, we
may expand ui for state i into a linear combination of
products of a radial function and a spherical harmonic
(here we use the real spherical harmonics for computational
convenience):
"U m
ur)= CmR(rr, s)Ymr)
where C are expansion coefficients to be determined, and
the radial functions Rj(rg,Ei) have to satisfy the follow
ing differential equations:
1 d 2d (z+l) 
r dr (ridr V a(p,r)} R (rE i)
rB dr5 B dr5 3S B x
= EiR(rg,Ei)
In the intersphere region, the potential is a
constant, VII, and the differential equation for ui in this
region is:
2 + (V i)} u(r) = 0
which is just the equation of a free electron. Hence we may
express ui as a multicenter partial wave expansion:
r Z (ir)Y(r8) + A0 j(ir0)Y (ro)
ui() = if i >VII
Ai k (Kr )Y(r ) + A0 i ( rOYM r
paI Y, m ^U (kiro 0 z ro
if ci
where 1i (IViI i) '
n is the ordinary Neumann function,
j. is the ordinary Bessel function,
k. is the modified spherical Hankel function of the
first kind,
i. is the modified spherical Bessel function,
and Aim are the expansion coefficients to be determined.
With the requirements that the spinorbitals ui and
their derivatives have to be continuous on the boundaries of
the spheres, we obtain a secular equation from which the
eigenvalues Ei and the expansion coefficients { CZ' I and
i8
{ A m can be determined. Appendix B gives the details
of the procedure to determine the eigenvalues and the expan
sion coefficients.
2.4 Relation Between the Exact Xa Method
and the MuffinTin Approximate Xa Method
In this section, we shall discuss the relation
between the exact Xa spinorbitals ui and the spinorbitals
ui derived from the muffintin approximate oneelectron
equation (henceforth called MSXa spinorbitals). First of
Fig. 2,2 Flow Chart of the Conceptual Exact
Xa Method and the M6Xa Method.
 Exact Xa
0 MSXc
> NonMuffinTin Calculation
olve exac
Xc eq. to
obtain u 
I
Calc.
er Form p
Fp ),p converg..
'/
Generate
Vxa(
no
all, we shall look at the differences between the two self
consistent cycles for the conceptual exact Xa method and
the muffintin approximate Xa method (figure 2.2).
Both SCF cycles start with the superposition of
atomic potentials. In the exact Xa method we try to solve
the Schrodinger equations (2.3) to get ui directly from this
potential whereas in the MSXa method we solve the Schrodinger
equations (2.6) to get ui after the potential being muffin
tin averaged. Then the charge density p is formed from
either ui or ui. In the exact Xa method the potential is
formed from p while in the MSXa method the potential is
formed from the muffintin part of the charge density,
namely T. Then we repeat the cycle again until the conver
gence is achieved. Then the total energy is computed with
p or p as the charge density respectively in the exact Xa
method and the MSXa method.
As we have mentioned in the previous sections, the
basic difference between the exact Xa spinorbitals ui and
the MSXa spinorbitals u. is that they minimize two differ
ent total energy functionals, and . In
other words, we have two different hamiltonians: H, the
hamiltonian for the exact Xa method, and H, the hamiltonian
for the MSXa method. The ui and ui are the eigenfunctions
of H and H respectively.
Since we had employed the approximate hamiltonian
in our calculations, an important question should arise:
how is the quality of the MSXa spinorbtials ui, or how
good are they compared to the exact Xa spinorbitals ui?
The calculations of the nonmuffintin total energy by using
the MSXa spinorbitals ui done by Danese (see Danese(1973))
have given part of the answer to this question (we shall
give a more detailed discussion about this calculation in the
next chapter). In this dissertation, we want to look at the
quality of ui in two more aspects: the HellmannFeynman
force and the dipole moment calculations.
It was proved by Slater that the exact Xa spin
orbitals (with the restriction that the a parameter be a
constant) satisfy the HellmannFeynman Theorem. That is,
the electrostatic force calculated with ui should be equal
to the derivative of the quantum mechanical expectation
value of the hamiltonian H with respect to the coordinate of
the nucleus upon which the force is exerted. Now, instead
of ui, we have ui; and instead of the exact Xa total energy
, we have the nonmuffintin total energy .
We can check if the ui satisfy the HellmannFeynman Theorem
and test their quality.
Also, dipole moment calculations are performed on
some diatomic molecules by using ui as the wavefunctions. As
we know, the dipole moment is much more sensitive to the
wavefunction than the energy is. Therefore, the dipole moment
calculations should be a good check on the quality of the
spinorbitals too.
Since we do have some reasonable results for the
total energy by assuming the muffintin charge density p, we
wonder what would be the result if we assume the same kind
of approximation in the force and dipole moment calculations.
So, before we discuss the force and dipole moment calcula
tions by using the MSXa spinorbitals ui, we shall first
investigate these two kinds of calculations by assuming the
muffintin charge c(nsity p in the next chapter. We shall
see that the results are not satisfactory and hence it is
necessary to use the p formed from ui, instead of the crude
approximation p, in such calculations. The distinction be
tween these two different approaches is also shown in figure
(2.2).
CHAPTER III
MUFFINTIN CHARGE DENSITY APPROXIMATION
3.1 Introduction
In the last chapter, we discussed the exact Xa
oneelectron equation. Instead of solving this exact equa
tion, we solve the muffintin approximate Xa oneelectron
equation which is obtained by varying the Xa total energy
functional by assuming a muffintin charge density approxi
mation. Since we had assumed this approximation, we would
like to see how good it is. One way to check the validity
of the muffintin charge density approximation is to look
at some of the molecular property calculations using the
muffintin averaged charge density. In this chapter, we
shall discuss three such calculations: the total energy
calculations, the electrostatic force calculations and
the dipole moment calculations. With the results obtained,
we conclude that in most cases it is not satisfactory to
use the muffintin charge density in the molecular calcula
tions. Also, we shall discuss the nonmuffintin total
energy calculations done by Danese. With the success of
his calculations, it seems that we should use the charge
density generated by the spinorbitals ui, rather than the.
muffintin charge density, in the molecular property calcu
lations.
3.2 MuffinTin Xa Total Energy Calculations
As we have shown in the last chapter, we calculate
the Xa total energy by using the muffintin charge density.
In solid state systems, especially the closepacked system,
the use of the muffintin approximation is usually quite
successful. However, for molecular systems, it is not uni
formly good. Although the absolute values of the total ener
gy are quite close to the experimental results, we have some
bad results when we consider the dissociation energy and
the conformation of the systems. The calculation by Connolly
and Sabin(see Connolly and Sabin(1972)) on the water mole
cule led to a minimum total energy for a linear molecule
rather than for the bent form known to be correct experimen
tally. The work in many diatomic molecules indicates no
binding too.
It is noticed that the results are worse in loose
packed systems than in closepacked systems. It was suggest
ed that this is because the intersphere region is larger in
the previous case than in the latter case, and the muffin
tin effect is correspondingly more severe. This was verified
by the work done by Danese and Connolly, Danese performed
the calculations for the nonmuffintin Xa total energy of
carbon and neon molecule by using the MSXa spinorbitals
rather than the muffintin charge density (see Danese
(1973)). The result is dramatically improved for C2. One
finds binding very much closer to the experimental value
than in HartreeFock calculation, and somewhat closer than
the result of the configuration interaction calculation.
For Ne2, one also finds binding.
With the calculations on the total energy, one may
conclude that the approximation of charge density has quite
a large effect on the molecular properties. In the later
sections, we shall show such muffintin charge density effect
on the electrostatic force and the dipole moments.
3.3 Force Calculations
In the last section, we have shown that the use of
MSXa spinorbitals ?.causes a big improvement in the total
energy calculations. In this section, we shall look at ui
in a more elaborate way.
As we know, the error in the total energy calcula
tion is of second order in the error in the wavefunction.
We would like to look at some molecular properties which
are more sensitive to the wavefunction. The force exerted
on a nucleus by the electronic charge cloud and the other
nuclei is a good test of the quality of the wavefunctions.
The HellmannFeynman Theorem, which we shall discuss in
more detail in the next chapter, provides a way to calculate
this force. It was proved (see Slater(1972b)) that this
theorem holds for exact Xa spinorbitals u. (with the re
1
striction that the a parameter be a constant). However, it
is well known that for approximate wavefunctions the theorem
often gives rather unreliable results. Thus, it should be
a good test of the MSXa spinorbitals by employing them in
a calculation of the force exerted on a nucleus of a mole
cule by means of the HellmannFeynman Theorem.
We shall discuss such a calculation in Chapter IV.
Before that, we may wonder what the result would be if we
use the muffintin charge density in such a calculation,
just as in the total energy calculation. We shall show in
the following that for homonuclear diatomic molecules, the
force obtained with such a muffintin approximation would
be always be repulsive. This cannot be physically possible
because that would mean there is no binding for such a
system.
In the MSXa calculation of a homonuclear diatomic
molecule, the space is partitioned into three regions: two
atomic spheres surrounding the two nuclei A and B, each with
nuclear charge Z, an outer sphere enclosing the two atomic
spheres, and the intersphere region II, as shown in figure
(3.1).
By symmetry, the force on nucleus A is along the z
axis and is given by:
F = 2Z P() cosOA d3r 2 (3.1)
r A
A AB
Fig. 3.1 Partitioning of the Space of a Homonuclear
Diatomic Molecule in the MSXa Method.
where p is the charge density in atomic unit,
eA is the azimuthal angle measured with nucleus A at
the center,
rA is the distance from nucleus A,
RAB is the separation of the two nuclei,
FA is in Rydberg/a.u. units.
Now, if the muffintin charge density approximation
is assumed, the quantity p in equation (3.1) will be replaced
by p and the calculation for FA becomes much more simple.
First, we decompose the integral into four different parts,
one over each region:
Sp(r) cosOA 3
r dr = IA + IB + out +II
where I p(r ) coseA d3
where IA = d2 r
A JA rA
p(r) cosA 3
IB r 2 d3r
I p(r) coseA 3
out d r
Sp(r) coseA 3
III =  r d3r
1 i rA
We shall discuss each of the above integrals in the follow
ing:
Integral Over Sphere A: IA
Since p(r) is spherically symmetric inside sphere A,
we have:
Sp(r) cos9A 3
IA 2= r d3r
A
= T(r) rA drA cos sin O dSA d
where RA is the radius of the sphere A and is equal to one
half of the separation of the bond distance R.
Since the integral over the coordinate 0A is equal
to zero, we have IA = 0.
Integral Over Sphere B: IB
To evaluate the other integrals, we need to expand cosA
in terms of functions centered at nucleus B or at the cen
ter of the outer sphere. Such an expression can be derived
from an expansion theorem for Neumann functions (see Appen
ix F). The expression is as follows:
0
S9+1,(RAR B)
4t I(t+1,0;1,0;1,0)  r Yo:rB)
Y RB B Z BAB
Y(r RAB
 = (if rB < RAB) (3.2)
i Y (r~+)
4 I(10;11,0;,0) RAOZ1 Y O(RA YB +
(if r0 > RA0) (3.3)
where Ym are the real spherical harmonics,
I(. ,m ;Z ,m ;Z ,m ) are the Gaunt coefficients for
1 1 L 2 3 J
real spherical harmonics and are defined by:
Y I(r) Y 2(r) Y (r) d (3.4)
t i Z,2 9. j
and RAB, RAO are the distances between A.and B and
between A and 0 respectively.
Inside sphere B, we have rB < RAB and we can use
equation (3.2) for the evaluation of IB:
Sp(r) cos9A
I = r2 d3r
B B rA
= p(r) (47/3) Yl(r d3r
J rA
= (4i7/3)' p(rB)(4Ir)" I(t+1,0;1,0;, 0)
+ 
Y+1(RARB) o 3
rx +2 B Y(rB) d r
RAB
Y 01 AAB)
=(4/3) (4~) I(t+l,0;1,0;,0) eR+2
x RAB
R +2 O
x' (rB) rB drB Y(rB) dQB
0 r
r (4u)Bg) dQB ko) Q
I = (4Tr/3)2 (4i) I(1,0;1,0;0,0)
B
Since
and
Yl(RARB)
R'AB
x (4T) j P(rB) rB drg
I(1,0;1,0;0,0)=(1/4i)
00 ( (3/4)
Y (RA RB) = Y(0A = 1) = (3/4r)2cos =(3/4a)
1 2 2
therefore I P(r)4r2 dr = Q /R
B R2 B B B B AB
AB
where QB = IP(rB)4T rB drB is the electronic charge
inside sphere B.
Integral Over the Outer Sphere: Iout
In the outer sphere, ro > RAO, so we can use equation
cos0A
(3.3) for the expansion of " Then we have:
rA
S cosOA 3
Iout = Jo p() d3r
A
= (4T/3)2 P(r )4i I(1,0;1,0;,0)
Z1 0 Y( 0) 3
x RAO Yk1 RAR0) t+ d r
r0
= 4n(4/3) I(l,0;l,0;,0) RA0 Y RA 0
p(r0) 2 0
O 1 r0 dr0 Y(ro) d$0
0 1 )
Since Y(r ) d0 = (4 ) 60
and k j 0 in the summation, therefore Iout = O
Integral Over the Intersphere Region: III
In the intersphere region, the muffintin charge
density is a constant, so we can take p out of the integral:
If p( ) cosOA 3 
S= d r = p
I I r I
We further decompose the integral into three parts:
ScosA 3
J d r =
rI A
cosA 3
 d3r
rA
cose d3
6 A
Just as for IA and IB, wehave:
ScoseA
A A d3r = 0
A rA
cosA
d3r =
and 1
2
( 4r2 dr )
( B B
2
/RAB = 1TRAB/6
For the first term, we have:
ScosA 3
7 d3r
AtB*tII rA
= (4x/3) Y(rA) d3r +
SA
(4T/3)1 {(4T)
Azoo
JSr+2dr2
0 '
rY(A) d3r
, rat R A
y 0(AAO
I(Z+l,0;l,;Z,)ZA
KO+2
AO
YO(rO) .dQ
Z1 0 
+4,r I(zl,;lO;Z ,) RAO YZ1l(RARO)
P1
cosA d3
rAdr
rA
cos0A 3
 d3r
rA
29
.x 2   dr0 Y0(r0)d2O}
AO r0
= 4rRA0/3 = 2TRAB/3
Putting all of them together:
III= =II ( 2nRAB/3 TRAB/6 ) = p 7RAB/2
We can express II in terms of the total charge QI1
in the intersphere region and the volume VII;
3 3 3 3
Vi = (47/3) (RAB RA RB) = TRA
QII QII
Therefore I = V II RAB/2 = 2
II 2 RAB
Combining IA,IB,Iout' II and the nucleus force, we have:
A cosOA 3 2 2
FA = 2Z (r)  dr 2Z2/RAB
= 2Z (IA + I + It + I ) 2/RB
_B QII 2 2
=2Z (0 + + 0 2 ) 2Z /RAB
AB AB
2Z
= 2 (% + QII/2 z)
RAB
Since for a homonuclear diatomic molecule, we have:
QII = total 2 QB out
where Qout is the charge in the outer sphere.
out
Also, for a neutral homonuclear diatomic molecule,
Qtotal=2Z. Thus, we finally have:
2Z
F= 2 (Q + Qtota/2 Q Q /2 Z)
A R B total B out
2Z
AB
= AS <ut2 + Z Z)
Z out (3.5)
AB
Since Qout is always positive, we have the result
that the force exerted on the nucleus A is always negative,
or away from the other nucleus. Thus we prove that the force
is always repulsive for a homonuclear diatomic molecule if
we assume the muffintin charge density.
Due to such a result, we see that it is inadequate
to employ the muffintin charge density approximation in the
force calculations. As in the total energy calculations, we
hope that improvements may be achieved if we use the charge
density p generated by the MSXa spinorbitals, instead of
the muffintin charge density. Indeed, we have some improve
ments as we shall discuss in Chapter IV.
3.4 Dipole Moment Calculations
Another way to check the quality of the wavefunctions
is to look at the dipole moment calculations. As is well
known, the dipole moment is very sensitive to the wavefunc
tion employed in the calculations. We shall use the MSXa
spinorbitals for such an investigation. Before that, we
first derive an expression for the dipole moment calculated
from the muffintin charge density. We shall see that such a
charge density approximation is also inadequate, and we need
a more elaborate way to evaluate the dipole moment.
Consider a neutral diatomic molecule AB. We partition
the space as we do in the force calculation. Since we have a
heteronuclear diatomic molecule, sphere A and sphere B are
in general of different size.
The dipole moment is defined as:
= p(r)r d r + vN
where P(r) is the muffintin charge density,
UN is the nuclear contribution and is defined by:
UN = ZARA + ZBRB.
For a neutral molecule, it can be shown (Appendix D)
that the dipole moment is independent of the choice of the
origin. Thus, for computational convenience, we can choose
the origin at the center of the outer sphere. Furthermore, by
the axial symmetry of the system, we note that the only non
vanishing component is the zcomponent. Thus,
S= j (r)z0 d r + N
32
To evaluate the integral, we decompose it into four
parts again:
S(r)zO d3r
= p(r)z0d r + j (r)Zq d3r + () dr
+ J p(r)zo d3r
We shall consider them one by one:
33
(1) J' p(r)z0 dr
= J p(r) (rAcos0A RAO) d3r
p= (rA)rA drA 1 cosGA dA RAO I (r) d r
= 0 RAQA = RAQA
(2) p (r)zo d3r
= ] P(r) (rBgosOB + RB0) d3r
e 3
S (r )r drB cosoB dQB + RBO p(r) d3r
= 0 + RBOQB = RBOQB
(3) Iot P(?)z0 d3r
= lt o(r0)r0cosO0 d3r
= r p(r0)r dr0J cos00 d00 = 0
(4) J (r)z0 d3r
= PII 11Z d3r
= ( z0 d3r z0 d3r )
0 A ZO d ZO d
= { r" 3 dr0 Jcs0 d0 I (rAcsOA RAO) d3r
I (r cos + R ) d3r
J B B BO
= II ( 0 + RAOVA RBOVB )
where VA and VB are the volume of sphere A and sphere B
respectively.
Combining all the terms, we have:
I p(r)z0 d3r = RAQA + RBOQB + II(VARAO VBRBO
Since RAO = (RA + RB) RA = RB
RBO = (RA + RB) RB = RA
VA = 4r Ri/3
3
VB = 4r RB/3
PI1 = QII/VII = QII/ 4ARARB(RA+RB)
Therefore p(r)z0 d3r
= RBQA +
Also,
QII 4R 3 3RA
RAQB + 4RRARB(RA+RB) W (ARB RBRA)
= RBQA + RAQB + QII(RA RB)/3
PN = ZBRBO ZARAO = ZgRA ZAB
Thus, finally we have:
= { RBQA + RAQB + QII(RARB)/3} + ZBRA ZARB
= RB(QA + QII/3 ZA) RA(QB + QII/3 ZB)
(3.6)
We shall discuss the results computed with the above
equation in Chapter V. There we shall find that most of the
results are not satisfactory when compared with the experi
mental values or with the results obtained by other methods.
35
The reason is simple: the expression (3.5) depends entirely
on the crude approximation of muffintin charge density. We
hope that improvements may be achieved by using the charge
density generated by the MSXa spinorbitals instead of the
crude muffintin density. We shall devote Chapter V to such
calculations.
CHAPTER IV
HELLMANNFEYNMAN FORCE CALCULATIONS
4.1 Introduction
In the last chapter, we saw the necessity of using
the nonmuffintin charge density generated from the MSXa
spinorbitals ui in the force and dipole moment calculations.
To carry out these calculations, we need to evaluate a three
dimensional integral whose integrand is a multicenter func
tion. The basic difficulty of such calculations is the eval
uation of the integral over the intersphere region, which is
in general very awkward. In this chapter, we shall develop
an analytical way to perform such integrals for the case of
diatomic molecules, and we shall apply it to the force calcu
lations.
4.2 Prolate Spheroidal Coordinates
Before we go into the force and dipole moment calcu
lations, we first introduce the Prolate Spheroidal Coordinates
that are essential to our scheme of evaluating the integral
over the intersphere region.
Let A and B, with separation of R, be the two foci
of our coordinates. Also, let (rA, OA, iA) and (rB, 6B, (B)
Fig. 4.1 Representation of a Point in the Atomic Coordinates.
be the two spherical coordinates with centers located at A
and B respectively. OA and OB are both measured from the
major axis in the same direction, as shown in figure (4.1).
The Prolate Spheroidal Coordinates are d: fined .as:
5 = (rA + rB)/R
n = (rA + rB)/R
(4.1)
A = CA = dB
SA B
The reciprocal relations are:
rA = R( + n)/2
cosOA = (1 + +n)/( + n)
{A = 2 )(1 n2)}
+ n
rB = R( n)/2
and cosO = (1 5n)/(5 n)
= (2 1)(1 n2)
sin(B
(4.2)
Also, it can be shown that the volume element is:
d3r = (R/2)3(2 n2) d dn d4
(4.3)
If the integral is over the entire space, the inte
gration limits would be: 1 $ ( < <, 1 I q $1, and 0 2ir .
However, if we are evaluating an integral over a finite space
only, such as the integral over the intersphere region, the
integration limits have to be modified. We shall use the two
atomic centers or one atomic center together with the center
of the outer sphere as the two foci, depending on the form
of the integrands, in our treatment of the integrals.
4.3 HellmannFeynman Theorem
It was proved (see Slater (1972b)) that the exact Xa
spinorbitals, with the restriction that the a value be a
constant everywhere, rigorously obey the HellmannFeynman
9
Theorem. In other words, the force component
aXp
acting on a nuclear coordinate Xp, where is the exact
Xa total energy of the system as calculated quantummechani
cally from the exact Xa spinorbitals, is the force which
would be calculated by pure classical electrostatics acting
on this nuclear coordinate, resulting from the electronic
charge cloud, formed by the exact Xa spinorbitals, and all
other nuclei. We can express the theorem in the following
form:
9 ( 2Z ) d3r Z 2ZpZq
Xp x (r) ( ) d r  ( T p )
DX pJx p P @Xp RpRqI
(4.4)
th
where Zp is the atomic number of the p nucleus,
is the exact Xa total energy in Rydberg units,
Pxa is the electronic charge density formed by the
exact Xa spinorbitals,
R R are the nuclear coordinates.
In the above theorem, we want to emphasize that the
spinorbitals and the total energy are those derived from
the exact Xa method. As we have mentioned in the previous
chapters, the existing MSXa calculations are not exact Xa
calculations. They assume the muffintin approximation. With
such an approximation, we can ask the following two questions:
(1) Does the HellmannFeynman Theorem hold for the
approximate Xa spinorbitals ui? If not, then how severe is
the discrepancy? In other words, we want to see how large
is the difference between xa and Fp(pxa), where
Xp
is the Xa total energy obtained with MSXa spinorbitals ui
and Fp(pxa) is the electrostatic force exerted on the nuclear
coordinate Xp by the other nuclei and the electronic charge
cloud pxa formed from ui.
(2) How well does Fp(pxa) compare with the experimen
tal value of force acting on the nuclear coordinate Xp?
In order to answer the above questions, we made some
force calculations on homonuclear diatomic molecules. The
reasons that we pick the homonuclear diatomic molecules are
that the a values in these systems are constant throughout
the whole space which is the restriction of the Hellmann
Feynman Theorem, and that the computations are easier and
the results can be compared with the experimental values.
4.4 Method of HellmannFeynman Force Calculations
To investigate the HellmannFeynman Theorem, we have
to calculate two quantities, the total energy and the elec
trostatic force exerted by the electron cloud as expressed
by the integral on the left hand side of equation (4.4). The
calculation of the total energy can be done by using
Danese's program (see Danese (1973)), as we have discussed
in Section 2 of Chapter III. The main quantity that we wait
to evaluate in this discussion is the electrostatic force
Fe (xC)
2Zp
Fe 'a) = r) ( _7 ) d3r
3Xp IrRpl
Since p(r) = niui(r)ui(r) where ni is the occupa
tion number of the ith state, we can express Fe (xa) as:
i p I,2Zp
Fe(xa) = ni uI (i)i() ( TX ) d3r
p r p1
For a homonuclear diatomic molecule, the expression
for the zcomponent, which is the only nonvanishing part,
of the force exerted on nucleus A is reduced to:
S1* t > cosOA 3
Fe (xa) = 2 ZA ni i(r)ui(r) dr
i rA
To evaluate the above integral, we partition the
region of integration into four parts: one over sphere A,
one over sphere B, one over the outer sphere, and one over
the intersphere region:
S%* % cosGA 3r
ui() i(r) c d r = F + F + Fot + FII (4.5)
i U COSOAG 3
where F = ()ui() d3r j = A, B, out, II
We shall discuss each of the integrals in the follow
ing:
Integral Over Sphere A: FA
Inside sphere A, we have
a + iA +k m 9
ui(r) = C R (rA,,i) Y A(rA)
iA
where C.m are the Ccoefficients,
km
Ym(rA) are the real spherical harmonics centered at A.
Since we have used real functions as our basis, all
%* t
the Ccoefficients are real, and ui = ui. Therefore,
_i f cos6A 3
FA = ui(r)ui(r) OS d3r
A rA
= Cmi CiAm ( RyZ(rA)Ym (rA)
Ai f7, JA m 1 i
xR (rA)Y (rA) C }.d3r
r A
S iCA CiA fA R(rA)RZ(rA) r2 drA
R, miPr~rA
SmY r my (rA)COSEA dQA
x Y A 2 A A A
= CA C Am,(4t/3) I((4/2I2,m2;mlO )
IRR (rA)R (rA) drA
where I(L, ,m; 2,m2;i ,m3) are the Gaunt coefficients as
defined in equation (3.4).
Since I(t1,m ;Z2,m2;1,0) = 0 unless the arguments
satisfy the following conditions:
m = m = m
1 2
S+Z +1 is even
IL, I < 1 < Z + 2
therefore FA = C m C Am (4r/3)t
A M Yi 2 M
x I( ,m;2,m;1,0) R R(rA)R (rA) drA
By the axial symmetry of the diatomic molecule, the
Ccoefficients are nonzero only if m=mi for state i (hence
forth called m). Therefore, we have:
Fi = C iA (4i/3)1
CA 2 , m 2m
x I(,m;Z2,m;l,0) J RZ (rA)R 2(rA)drA
(4.6)
i
Integral Over Sphere B: FB
Inside sphere B, we have:
r C iB y(rg (B)
ui(r) = C R(rB) rB
So FB C Z C J2 R (r )Y (rB)R (r )Y (r)
B Is m B B 2 B B B
cos0A 3
A
2
By equation (3.2), we can expand cosCA/rA in terms
of functions centered at B:
0+
cosegA YI(rA)
r2 = (43/3) I rA
rA rA
= 47 (47/3)2 I( Z+1,0;1,0;l,0 )
0
S+1 (RARB)L Y0
R+2 B B
AB
Therefore F = iBC iB(4)(4/3)
S I(+l,0;l,0;l,0) Y (~1A+B)
d+2
AB
x R (rB)RL (rB)rg drB
X Y 1(rB)Y (rB )y(rB) dQB
j iB BB
F = C iB CiB (4T)(4iT/3) {I(+1,0;1,0;Z,0)
B m 2M (4(4/
x I(e,m;2,m;,0) Y+1(RARB)
AB
x R(rB)R2(rB)rB+2drB1 (4.7)
*/
Integral Over the Outer Sphere: Fout
Inside the outer sphere, we have:
ui (r) = C R (r0) Y(r
Thus Fut = C 00 {R (rO)Yi(r)R (r)Y (r)
out 1m 92m LA t xOOkOz
cosA 3
x 2 d3r
rA
By equation (3.3), we expand cosOA/r around the
center of the outer sphere 0. After integrating over QA'
we have:
Fi = 4(4w/3) C c0i 0 {I(,1,0;1,;0',o0)
out m '2m
A '.i, 1atPt"
x I(,,m;,m;,) RAO YIAr
R r dg (4.8
Ro r0) ro
i
Integral Over the Intersphere Region: FII
In the intersphere region, ui can be expressed as:
iA m + iBA
iA n (KrA)Y (rA) + A r n (KrB rm
2i1 Az1 A I 2 2 2
iO m .
mui() =
AA km (CrA)Y1 (rA) + TAiB k (rB)Yk (rB)
21 1 2 2 2
io )m
*+ A0 i (Kr0)Y (r0) if Ei < V
For diatomic molecules, the eigenvalues Ei for most
states (except the core.state) are greater than the constant
potential VII in the intersphere region. We shall look at
the positive energy casd (i>VI ) in more detail (the other
case can be treated in a similar way). For core states, the
electronic charge is concentrated around the atomic centers
and the contribution due to the intersphere region can be
neglected in the core states. Thus, we have the following
i
expression for FII:
F + cosOA 3
= ui(r)ui(r) dr
A
'iA AiA i 1 + AiB A iB i
= A m A2m IAA 1Z2m) + m Am BB ( 2m)
+ iO i A A Om I00( 2,m) + 2 A iA iA B i '
Rim ,m m) 2tAgm A2m IAB(i,2,m)
+ ,.i t.
iA iO i iB iO i
+ 2 AmA2mlAO (1, ,m) + 2 mA BO( ,,m)
(4.9)
i If m
where IA (' m) = in1 (KrA ) (rA)q 2(
cos0A d3r
A
rA
i m m
IAB ( 9M ,m) = An (
I i
COSOA d3
x  dr
A
i ((KZm) = Ii (Kr )Y mr ()
AO (i A r ,(A ) 22(r 0 2 0
xcosOA d3r
r2
rA
I 0(,)' 2m) 1 2i'n(rB)Y(rB) (Kr0Y(rO)
c osA 3
x d3r
r2
rA
i m m
I0O(1 2,m) = f 1(Kro)Y1(ro)j2(Kro)Y12(r)
cosOA d3r
x  dr
rA
To evaluate each of the above integrals, we employ
Prolate Spheroidal Coordinates as introduced in Section 2
at the beginning of this chapter. First of all, we want to
see what the limits of integration are for such coordinates.
For the intersphere region, we have the following
restrictions on E and n:
(1) rA + rB > R
(2) 1 rA rB I < R
(3) rA > RA = R/2
(4) rB > RB = R/2
or 5 = (rA + rB)/R > 1
or I n = (rA rB)/Rl 1
or + 1
S > 1
(5) From figure (4.2), we have:
2 2 2
r0 = rA + (R/2) 2rA(R/2)cos0A
2 2 2
= (R/2) (+n)2 + (R/2) 2(R/2)(+n)/2)(R/2)(+n)/(+n)
= (R/2)2 (E2 + n2 l)
2 2 2
= (R/2) ( + n 1)
Since r0 H R, we have 2 + (.5
Putting all of the restrictions together, we have
Fig. 4.2 Relation of rA, R and eA
TI
n=/5iZ
n= /5
Fig. 4.3 Boundary of the Intersphere Region in
Prolate Spheroidal Coordinates with A
and B as the Foci.
the transformed intersphere region as defined by: (figure
4.3)
f(S) $ C $ f(S) and 1 4 C J5
where r 1 if 1 5 S 2
f() =
(5 5 ) if 2 C VS 75
Before discussing the evaluation of the integrals,
we note that all the integrands, except for the factor of
the spherical harmonics, are 0independent, and, because of
the cylindrical symmetry of the intersphere region, all the
integrals can be reduced into twodimensional integrals
through the 4integration, such as:
m 2 3
n (KrA )Y(rA )n2(KrB) (rB)cOSA/rA d r
[ mm 23
= rJ (Kr )PA (cosS )nZ (Kr )P2(cosOB)cosOA/rAd r/d4
where P'j are the normalized Legendre polynomials.
To evaluate the integrals, we also introduce the fol
m
lowing expressions for P ri and j :
m = N in)m m m2
P (coseA) NZ(sinSA) ) (cos0A
51
i+1 1
n (z) = (1) z {Q(i+,z)cos(z+7) Q'(a+,z)sin(z+inw)}
jk(z) = z1 {Qo(+,z)sin(ziir) + Q'(+,z)cos(zjit)}
m m (2+1l)(m)!
where N = (1) { 2 (i+m)!
im =(_ (2i2v)!
v ( 2(tim2v)!(iv)!v!
Q(i+L,z) = (l)k(2+,2k)(2z)2k (4.10)
k=o
Q'(Z+1,z) = (1)k(k+J,2k+1)(2z)2k1
Q'(i+jz)= Zh) 2k1
Z = (l)k(t+,2k+l)(2z)
Ico
(+k)!
(R+i,k) = k(k)!
k!(9,k)!
and [ b] = the integer part of the real number b.
With the above expressions and some algebra, we have
i i i i i
the following results for IAA, IBB' IAB' IAO' IBO, and I00
as required by equation (4.9):
(1) IAA: ni(KrA)Y (rA)n (KrA)Y2(rA)cos5A/rA d r
This is a onecenter integral. With some straight
forward algebra, we have:
52
I ( 1 (l) 2 (R/2) d I 2 d,
AA = N 2 95 Ak IO S.0
{W1 m zm(_1)rCm(KR) 1+ 22m+12v+2r C R1+Z2m+l2v+2r
V1 V2 r r
(KR)s s (1 2 L +Z22m+12v+2rs
x {4Sgil2 (z+2+ 22v+2rs+4,0)
(2/KR) gl (I1+Z22m2v+2rs+3,))}} (4.11)
112
where v = v,+ v2
m = m!
Cr = (mr)!r!
Here the g1 functions are defined as:
g1 (2( ) = 5 7 (l)U+V
k I P 2 3 0
x A(u,v) (u v1) +n+u+v
n 1 2
where the Acoefficients are defined in terms of the (k+,k)
coefficients through the following relations:
S A(uv)(2z)n = QU(t +,z)QV(e +2,z)
A1 s 0)
53
and P(uv,1)(,0) are defined in terms of the sine and
1 2
cosine integrals as below:
P(u,v,1)(,U) = Icos(Z1+2)r (6u+v, 6u+v,2
KR(t _cost))
+ sin (Z +2 ) 6 }uv,l ( ) t dt
j u+vo1 f(P)) t
+ {sini( +R2)T (6u+v,26u+v,O
rrR(j ( ~ sin t
+ cos((Z +z2 ) 6 sin dt
1 2 u+v,f) t
+ {cos(I z )7T (6 ++6 )
1 2) ('u+v,0 u+v,2
+ sin(1 2) 6u+v,1(6vO6v1)
i m m 23
(2) IBB: 1(KrB )Y1(rB)n2 )Y (rB)COSOA/rA dr
(2) IBB: T1/^B ZI B Z2 B Z2 B A A
This is a twocenter integral. The result is as
following:
i = N N (1)z + Z+1(R/2) d "
BB I 2 ki 2o Ro p' d0
Xikm Pm m +2i2m2v+3r ,+,2m2v+2r
x{Vw 2 r () Cs
x (KR)S 2(1_2) s +22m2v+2rs {(1+2) g 2(s+)1
(S/KR) g92 2(s,)~) (4.12)
XZ XZ2
where v = v, + v,
and the g2functions are defined as:
t (1U+V
g 2(s,) = t Y, E (1)
US= v0o ,=o
x A(uv) p(uv,2)(s+2n+u+v+l,E)
n 2
i
where the Acoefficients are defined as in IAA, and
P1u v,2) (,C) are defined in terms of the sine and cosine
integrals as below:
p ,2 (e,) = .{cosj(i +z )T (6 6 2)
1i2 1 2 u+v,O u+v,2
rkR(Pflq))
+ sinji(1+,)T6 } os t dt
u+v,1 J KI()) (2KRCt) t
+{sinI(il+i2)n (6u+v,26u+v,)
+ cosj(ZJ+iZ)76 }t) sin t dt
u+v,1 ) (2KREt)2t
+{cosSi(,,)T (6u+v,o u+v,2)
+ sin(i1 2 ) (6 vo6 )6+v,1
fR() rf))
JKRE 1 dt
J' (tqI) (2
55
i 2
(3) IAB: (Kr )Y (rA (KrA )Y2(rA)OS0A/rA d3r
This is also a twocenter integral. The result is:
I N N (1) ,2(R/2) f d
AB 1 2 .2 av 0 rV
u+v+ 2v2 CCm 2u C m2v .1m m
=(1) C v V1 V2
2 C2m2v2 ) +i22v 2v22u+vsr v+r+s
s (KR)
x{2(1 2) g3 (. _2v,2u+vr+3,Z2v2s,E)
Sg3 ,L( 12v i2u+vr+2,i,2v2s,C)}} (4.13)
where the g3functions are defined as:
gi (Z,(I ) = (1)u+v+k,+k2 ( 1+1,2k,+u)
111X 2 W=o Vso *=O arO
((u,v,3)
x ( +,2k +v) P1u ,3) (+2k1+u, +2k +v,0
and P2U v3)(r,s,) are defined in terms of the sine and
cosine integrals as below:
(uv 3) ,
= {cos(KRC+2I2 )(6 +v,o+6u+v
2 u+v,0 u+v,2
rkr(j+kp)
+ sin(KRE+ 9 2 )6u+v'l(6vO6vl)} f cos t dt
K(JfCp) (2K Rt)rts
+{sin(KRC+1 )(6 u+v, 0+u+v,2
+ cos(KR + ) 6u+i,1(6u0o6ul1 r sin r
z Z2 sin t
+ cos(KRC+12 u+,v,( uO6ul)l (2KR t)rts dt
+
2 rU()VU U (2cR62t) t
+{cos(KRS++ 2
2 ) (6u+v,O u+v,2
+ sin(KRE2 w)6u+v,l (2KRt)rts dt
(4)1 i I(r r (KO)008r/r 2 d3r
(4) IA :z n(rA )Y(rA)i( 2 (ro)cosOA/rA dr
Here we have a twocenter integral. We want to ex
m +
press j 2(Kr)Y 2(rO) in terms of functions centered at A so
that we can reduce the integral into sums of onecenter inte
grals. Such an expansion is derived in Appendix E:
j2(KrO)Y2(ro) = (1) 24 7 i
x I(",m" ;!,m;,ii' ,m' ) (1)'
x j, (KRA)Y (RA l(KrA)Y (rA) (4.14)
for r0 in the intersphere region
With the above equation, we have:
AO A<. '.
0 [ m m M
X j ,(
coseA 3
X Co dr (4.15)
A
i
Thus, we reduce I to an infinite sum of onecenter
AO
integrals. Since j decreases rapidly with increasing 2',
we can expect that the sum would converge rapidly.
To evaluate the integral on the RHS of equation
(4.15), we employ the same technique as before and obtain
the following expression:
mI m 2 3
J (Kr)Y1 A) (K A) 2 (rA)coseA/rA d r
N i(KrA)YtN( r A) d({A)YZ(rAIsXA/
r + 1 L v r
Z ,P 2t 1\ Yt Y' =0D S=O
x(1)r m 2m Cm C c+22m2v+2r+l (KR) 2 )
) 2 r S (KR) (1i )
V 1 V 2 r S
XE1+z,2m+12v+2rs 2 4E g (s+3,) 4 (s+2,)}
Sz KR (s+2,)}}
(4.16)
where v = v + v2
and the g4functions are defined as:
SI ()uA (u, v) (u v,4P )(+2n+u+v,,)
Z 1 P 2 XAL 2
58
where P(u v,4)(,E) are defined in terms of the sine and
e i 2
cosine integrals as below:
S(,,4)(,) = {sini(L1Z2) (6 +6U )
1i2 u+v,O u+v,2
(KR(fwfP)
+ cos(Z,2) 6u+v,l(6 uul )} tdt
+{cosi(12)> ( ++v, ++v)
u+v,l ul uO KR[J(p) t
+{sin(1+ 2) (6U+v,26u+v,O)
+ cosi(+a). 6+u2 ) dt
i ( m m 23
(5) BO: (Kr )Y (rB) (Kr )Y2 (r)cosO A/rA d r
(5) .IBOB:2 0 d0rA
m
Expanding j 2(Kr )Y2(r0) in terms of functions cen
tered at at B, we have,similar to IA0, the following equation:
I1O = (l)24t 4 (1)k I(",m;,m;',0)
m + r m + .+
x j,(KRB)Y,(RB) n (KTr)Y (r )j,(KrB)Y,,(rB)
x cosOA/r d3r
59
The integral on the RHS of the above equation can be
evaluated as below:
2 3
11 (rB )Y (rB). ,(KrB )"(rB)cosEA/r d3r
= Nm ( ) +1(R/2) d( L
k1 2 I *, =0 1ro :O
+Z2m2v+3r 1m ,m m C1+2m2v+2r 2(KR)s
x(l) w Cr s
V1 V2 r s
x.(l2)s ,+122m2v+2rs {(i+2 )gS (s+l,) g (s,O)}
12 cKR Z2
where v = v + v
1 2
and the gsfunctions are defined as:
gkt (S,) = t Y_ (I) An (s+2n+u+v+l,C)
g (s, = Az ((unu (uvu) (uv,5) u
1L2 a2l V. ". n
where the P (u,,5)(t,) are defined in terms of sine and
cosine integrals as:
P v,5)( .,) = sin(1,2)l (6+v,O+6u+v,2)
+ cos(Z2)r 6 (66()}
+ cosi2( ,z)T 6u+v,1(6uO6ul
x cos t dt
(2KR(t) t
+{cosj(1,z2)T (6u+v,O+6 u+v,2)
u+v,l ul u0
+ sin(t,2)n 6)u+v,(6ul
x I sin t dt
.I()5jt)2KRSt) 2t
+{sin( 1+2 )T (6u+v,26u+v,0)
+ cosk +k )TT 6dt
+ cosi(+2) u 6u+,1} v(2KR t)2t
(6) Ii J (KrO)Y (r)Z(KrO)Y (Or)cosOA/r d3r
Expanding j (Kro)Y (r0)about A, we have:
i 0 0 Z ,Z2,m
I100 = E js(
where C ,2 ,m
ZI k3
= (l) ,+ (47)2
{ (l) I(P2I",(,m;,,m;R',0)
(m)
x I(,,m;k; ,m;R 3,0) I 1,
Z" I Z.4
and I(m
Z", 4
( i m )y 3
,(rA)Y (rA ) (KrA ) ( )COSA/rA d3r
4 4(
= j I(t",m;4,m;3,0)
0 2 3
3 (rA)cOSA/rA d r
, J"(CKrA)Ji (rA)
nIt
('+(
VIi'A
Sa~tIa
4b O~,4h~2rh
Thus, I00 is reduced to a sum of onecenter integrals
which, by employing the Prolate Spheroidal Coordinates,can
be further reduced as below:
0 2 3
j1(KrA )J (
2 L 30 2v+l r
= (2a3+1)2 .i (R/8) d{ L 30 Cr2v (KR)
x(1 2 r 32v+lr 2 KRf j(t)j0( dt
x(l 5 ) 5 {26 dt
(<f))) r+1
(I/KR) j j (Lit) dt}}
fR(cOf1))) tr
To conclude this section, we note that:
(i) We have reduced the threedimensional multicenter inte
grals into a sum of onedimensional integrals whose inte
grands involve the sine and cosine integrals. The sine and
cosine integrals can be evaluated quickly and accurately and
can be called directly from the IBM Scientific Subroutines
Package. We shall discuss such evaluations in Appendix G.
(ii) It looks as though the expressions for I I I,
AA BB' AB'
i i i
I, I and I are very formidable, yet in the actual cal
AO BO 00
culations they are much more simple. This is because the Aval
ues used in the basis of the MSXa spinorbitals are usually
small, and in most cases, do not exceed 3. Thus, for example,
if we use s and p atomic orbitals (. = 0,1) for the a state
(m = 0) in the MSXa calculation, then by equation (4.11),
62
there is only one term in the summation of v,, v2and r since
(im)/2 =0 .
(iii) By partitioning the region of integration into differ
ent parts, we can find the contributions to the total force
by:the electronic charge in each of the regions and for each
molecular state:
Total= Ni
= n(F + F + F + F ut) + FN
where i represents the molecular state i,
FA is the force exerted by the electronic charge
th i i
inside sphere A for the i state. Similarly for F F
B' II
and F
out'
FN is the force exerted by the other nucleus.
4.5 Results and Discussion
The HellmannFeynman force calculations have been
carried out for two different systems: the hydrogen molecule
and the nitrogen molecule. The hydrogen molecule is investi
gated because it is the simplest diatomic molecule and accu
rate results had been previously obtained (see Kolos and
Wolniewicz (1964)) so that we can make some comparison. The
nitrogen molecule is investigated because its electronic
configuration for the ground stateis a closed shell config
uration in which case the MSXa method should work better.
Also, it involves the tstates so that we can see how our
method works for a more complicated case.
Hydrogen molecule: H2
(A) The ground state for the hydrogen molecule H2
is z+ and the electronic configuration is log. Forces for
three different bond distances are computed in order to
plot a force curve and to obtain the equilibrium separation
Re and the force constant ke. First, for each bond distance,
a nonspinpolarized MSXa calculation is carried out by
using s and p orbitals in the atomic regions and s and d
.orbitals in the outer sphere. Then, the electrostatic force
is calculated, according to the method described in the pre
vious section, by using the converged wavefunctions. The re
sults are listed in table (4.1), and the curve is plotted in
figure (4.4).
(B) In order to check the validity of the Hellmann
Feynman Theorem of the MSXa spinorbitals, the nonmuffintin
total energy of H2 is calculated. From the values
a
of , we obtain the values of DR at the three
different separations by least squares fit. The results are
listed in table (4.2). Two such curves are plotted, one for
the wavefunction with basis: i=0,1 in the atomicsphere, i=0,2
in the outer sphere; another for the wavefunction with basis:
P=0 in the atomic sphere, Q=0,2 in the outer sphere. Also,
DE
to compare with the experimental result, a curve with 
against R where E is the total energy calculated by Kolos
and Wolniewicz (see Kolos and Wolniewicz (1964)) is plotted.
Such an H2 calculation is considered to be more accurate
TABLE 4.1
1+
Force calculations For H2 ( g
2FAg 2Fag 2Fg 2FCg
A B out II
2
: lg )
g
F Ftota
N total
R=1.4 0.264 0.271 0.026 0.442 1.020 0.068
R=1.44 0.262 0.266 0.U26 0.417 U.965 u.046
R=1.5 0.257 0.258 0.025 0.403 0.889 0.004
Basis: L=0,1 in the atomic spheres
=0,2 in the outer sphere.
a=0.77627
R in atomic unit; F in Rydberg/a.u.
TABLE 4.2
3
Comparison of Fto and R for H2
total aR
(a) (b) '(c)
Sa 3
total R dR
(d)
D
3R
R=1.4 0.068 0.024 0.035 0.004
R=1.44 0.046 0.004 0.012 0.032
R=1.5 0.004 0.046 0.024 0.068
Re 1.495 1.43 1.46 1.401
ke 0.72 U.71 0.60 0.73
(a) HellmannFeynman Force with basis: =0,1 in the atomic
region; L=0,2 in the outer sphere.
(b) Basis: =0,1 in the atomic region, =0,2 in the outer
sphere.
(c) Basis: =0 in the atomic region; =0,2 in the outer
sphere.
(d) Calculated with the KolosWolniewicz wavefuncqion
Force in Ry/a.u.; R in atomic units;ke in Hy/a.u.
Fig. 4.4 Force Curves for H2
1. Basis: =0,I in the H sphere
Z=u,2 In the outer sphere
2. Basis: i=0 in the H sphere
X=0,z in the outer sphere.
0.06 X Experimental
Q Exa() 1
aR
SDExa() 2
aR
0.04 Fx (x )
U.02
0
.0 __
0.04
O.U6
0.08
0.10
1.5 separation(a.u.)
1.3 1.4
TABLE 4.3
iD'ependence of HellmannFeynman Forces for H2 at R=1.4 a.u.
2FAg 2FBg zF 2F (
A B out II
FN Ftotal
gl 0.265 0.271 0.026 0.442 1.020 0.068
;2 0.268 0.274 0.026 0.442 1.020 0.062
#1: =0,1 in the atomic region; a=0,2 in the outer sphere
#2: =0,1,2 in the atomic region; =0,2 in the outer sphere.
Force in Ry/a/u.
than that experimentally observed. The equilibrium separa
tions Re and the force constants ke for each of the curves
are obtained. The results are listed in table (4.2) and the
curves are plotted in figure (4.4).
(C) To check the sensitivity of the force with re
spect to the basis functions, two different calculations, one
with s, p orbitals in the atomic region while the other with
with s, p, d orbitals (both with s, d functions in the outer
sphere), are carried at R=1.4 a.u. The results are in table
(4.3).
Nitrogen Molecule: N2
Bader,Henneker and Cadehave done some force calcula
tions by using HartreeFock wavefunctions as the basis (see
Bader, Henneker and Cade (1967)). In order to compare these
two methods, a force calculation for nitrogen molecule is
also done by using the MSXa wavefunctions.
Both the MSXa and the HartreeFock calculations are
carried out at experimental equilibrium separation Re=2.07
1+
a.u. The ground state for N2 is g and the electronic con
figuration is l1 la 2a 2 21 3a The spin polarized a value
g u g u u g
for nitrogen, c=0.74522 is used, and the values used for
the different molecular orbitals are listed as in table (4.4).
Here we treat the lag and lau as core states. The results of
the force calculations are as in table(4.5).
Discussion
When we compare Fxa(p), Fx (p) and the experimental
values for the force, we find that we have achieved a big
TABLE 4.4
ZVaiues in the MSXa Calculation for N2 ( Zg)
0,1,2
0,2
0,1,2
1,3
0,1,2
0,2
* o0 and lou are both treated as core states.
atomic
outer
TABLE 4.5
HellmannFeynman Forces for N2 at R=2 a.u.
Calculated with MSXa and HartreeFock Wavefunctions
2 2 2 2 4 2
Io la 2 2 0o ir 30
g u g u u g
fA 0. 1.513 0.501 1.278 0.462
fB 7. 2.227 0.9b2 2.317 1.685
f 0. U.03 0.844 0.160 0.788
out
fI 0. 3.316 U.0u3 3.274 0.452
fxa 7. 7.059 0.J95 7.029 1.811
fHarFock 7.858 9.387 1.621 8.512 0.525
Ftotal(Pxa) = 1.996
Ftotal(HartreeFock) = 0.263
* Total electronic force
Force in Ry/a.u.
improvement by using the nonmuffintin charge density to
calculate the force. The force constant k and the equili
e
brium separation for H2 calculated with this method involve
errors of 1.4% and 7% respectively, while the Fx,(P) indi
cates no binding at all.
Ex(>) we
However, when we compare Fx(p) and R we
find that there is some difference between them. This means
that the HellmannFeynman Theorem is not quite correct with
the approximate MSXa wavefunctions. Some reasons may be:
(1) Since the wavefunctions we used are hot exact
solutions to the Xa Schrodinger equations. It also assumes
the muffintin potential approximation and the potential is
generated only from a muffintin averaged charge density.
Thus, it is only an approximate Xa wavefunction and we cannot
claim that it should satisfy the HellmanFeynman Theorem very
well. Actually, due to the polarization effect, the potential
on both sides of the atomic spheres should be different. The
potential in the bonding region facing the other nucleus
should be less shallow than the potential in the outer re
ion. Thus, there should be more charge in the bonding region
than what is expected by using the muffintin averaged poten
tial. This means that the actual force should be more attrac
tive than the force we got by assuming the muffintin poten
tial. This polarization effect should be especially severe
for the core state, for which the electronic charge accumu
lates near the nucleus.
(2) Besides assuming the muffintin averaged poten
tial, we also use only a finite set of basis functions with
small i values. This may not be important in calculating
the total energy but it may be important in calculating the
force, because the force is much more sensitive to the wave
function than the total energy is. We further note that:
(i) In the force calculation for the hydrogen mole
cule ( table (4.2) and figure (4.4)), the values of Fxa(P)
by using the basis (H =0,1; out =0,2)are closer to the values
of Exa) calculated with the basis (Hp=0; 9out=0,2) than
3R.
to that obtained with the basis (tH=0,1; zout=0,2). Further
more, we have the .following relation:
S3Ex(P)
F x() ( =0,1; =0,2) < xa (i =O; R =0,2)
xci H out 'R I out
< 3Exa() ( ; =0,2) < F
aR H=0; out experimental
This can be expected since if the error in Fxa(P) (9H=0,1;
9out=0,2) is of first order of the errror in the wavefunction,
then the error in E x(p) (ZH=0,1; Zout=0.2) is of second
order and that the error in E (p) (. =0; out=0,2) should
xCi H out
be larger than that of Ex (P)(ZH=0,1; 9out=0,2) but smaller
than second order.
(ii) In table (4.3), we see that when we increase
the Zvalues, the net force becomes less repulsive and agrees
better with 3Exa(P) The improvement is very slight. However,
3R
if we use only s orbitals in the atomic spheres,then FA=0
and the net force becomes more repulsive. Thus, the p orbital
is very important in the force calculations, although the s
orbital is almost sufficient for the total energy calculation.
In the force calculation for the nitrogen molecule,
we see that the MSXa forces have the same characteristic as
the HartreeFock forces. Again, the MSXa forces are too re
pulsive. However, the relative values of forces contributed
by the different molecular states for both cases are compar
able, and the ou state in both calculations contributes re
pulsive force. Of course, we should not expect our result
could be the same as that of the HartreeFock calculation.
What we should expect is just a general characteristic of
these two kinds of forces.
In conclusion, the use of MSXa wavefunctions instead
of the muffintin averaged charge density to calculate the
force makes a big improvement. It gives a reasonable equili
brium separation and a good force constant. However, due to
the muffintin potential and the truncation of kvalues in
the basis set, the existing calculation of MSXa wavefunctions
still has to be modified to satisfy the HellmannFeynman
Theorem.
CHAPTER V
DIPOLE MOMENT CALCULATIONS
5.1 Introduction
In Chapter III,we have discussed how to compute the
dipole moment of a diatomic molecule by assuming the muffin
tin charge density. In this chapter, we shall discuss the
same calculation but using the nonmuffintin charge density,
i.e., the charge density formed from the MSXa wavefunctions.
We shall first discuss the computational procedure and then
give the results for some diatomic molecules. The results
show that the use of nonmuffintin charge density indeed
provides big improvements over the use of muffintin charge
density.
5.2 Method of Dipole Moment Calculations
To calculate the dipole moment of a diatomic mole
cule, we need to evaluate the integral p()z d3r contrib
uted by the electronic cloud. Again, in the MSXa calcula
tions, we have:
p(r ) = Z ni ui(r) ui(r)
As usual, we partition the coordinate space into three re
gions: the atomic region, the intersphere region, and the
II out
Fig. 5.1 Schematic Representation of A Heteronuclear
Diatomic Molecule in the MSXa Calculation.
outer sphere region (figure 5.1). As we have mentioned in
Chapter III, the value of the dipole moment is independent
of the choice of the origin for a neutral molecule, we can
choose the center of the outer sphere as the origin. Then
we have:
3 i i i
p()z dr = n (A B + out + II) (5.1)
i %* 3
where PA = J ui(r)ui(r)z0 d r
outi i *3
B = ui(r)ui (r)z0 d3r
11 = z ui(r)ui(r)zO d r
and the dipole moment u is: P = uN e
where PN = ZBRBO ZARAO is the nuclear contribution,
Pe = J p(r )z0 d3r is the electronic contribution,
ZA, ZB are respectively the atomic numbers of the two
nuclei A and B; p is in atomic units.
i
Integral Over Sphere A: uA
Inside sphere A, we have:
A, iA M
u. (r)= Ci R(r)Y(r )
1 Zm 2 A 9 A
and z0 = A RA0 = (4/3) rAY(rA) RAO
Thus, we have:
1A = ui(r ) ui(r) Z d r
= i CImA R (rAY m (rA)RZ(rA)Yr (rA)
x { (4T/3) rAY (rA) RAO d r
iA iAO
= C ciA c iC I (4r/3)2 I(Z1,m;92,m;l,0)
RA 3
x R (r)R(r)R )rA drA
A 2
RAO 62 RZ(rA)R2(rA)rA drA }
i
Integral Over Sphere B: up
Inside sphere B, we have
ai(?) = : Cis Re(rg) Y (?)
an = Z (r) Ym(rB)
and z0 = zB + RB0 = (4T/3) rB Yl(rB) + RB0
Thus, we have:
i = 3
B =J ui(r) u(r) d r
iB iB
= P CB CiB (47/3)2 I(t,1,m;2,,m;l,0)
= C 9m 92m
x R R(rB)R 2(rB)r drB + RB0o 6 [R (rB)R2(rB)rBdrB}
io
Integral Over the Outer Sphere: out
Inside the outer sphere, we have:
r iO
u.(r) = C R (r ) Y (ro)
1 Lm z 0 z 0
i %* r" 3
Thus = u.(r) u() z dr
out 1 1 0
= iCC. 1 C ou R (r )m(our)R (rOI (r
x (47/3) roY0(r0) d r
= C C m (47/3) I(k1,m;Z2,m;l,0)
ZIm Z2m
x R (ro)R (r0)ro dr0
Integral Over the Intersphere Region: ui
In the intersphere region, we can express ui(r) as
linear combinations of products of Neumann functions and
spherical harmonics as discussed in the force calculations.
In the same way as before, we obtain:
uI = (r) ui(r) z d3r
= Z lAmA 2m AA 1'2,m)' + AmAm IBB(1,2'm)
+ A 0A Om0( 1, 2,m) + 2 AAA AiB I 1 2,m)
Y 2m 00 il.m I2m AB
+2 AiA Ai0 Io( m) + 2 AiB Ai0
+2J, A m 2 m IAO(Il.2'm) + 2 A iBm A Im BO( 2'm)
m A 1 2 2 2m BO(t 1 2
(5.2)
where I AA(, 2 m) = T (rA )Y (r (KrA 2 A )z d 3r
AA 1 2 ti(A rA).(A Z2 A/Z2rA 0
m m ) 3
IBB I 1Z 2m) = (KrB B (r B 2 (KrB rB) zd r
IAB(" Z2,m) = n., (Kr r A )^^2 B9 31B
IAB(,IAm) = (KrA rA)2(KrB)Y(rB)Zd r
Sm 3
AO ') =2' l (KrA)Ym(r )j (,r )Ym (r )z d r
II(ll 2'm) 21A 2.1 A 2.2 0 9 2 0J0
IO( 2m) = I (Kr)Y ( )j (Kr )Ym ( )z d r
BO J I B i B 2 0 22 0 0
IB0 1' 2' j(Kr )Y B( 22 (Kr )Ymr)d
2m Z 1 (0 kr1 0 Z2 0 2.)2 0 0
To evaluate the above integrals, again we employ the
Prolate Spheroidal.Coordinates as we did in the force calcula
tions. However, since in this case we have a heteronuclear
diatomic case, so the limits of integration when transformed to
Prolate Spheroidal Coordinates would become more complex.
We shall discuss each of the integrals in the following:
() IAA 1(Kr A Yl rA 2 A Z2 A 0dr
This is a.twocenter integral. We shall express the
integrand in terms of Prolate Spheroidal Coordinates with A
and 0 as the two foci (figure(5.1)):
S= (rA + r )/RA
n = (rA r0)/RAO
S= A = 0O
The reciprocal relations are:
r = RAO ( + n) r0 = AO ( )
cosO = (l+5n)/(4+n) and coso = (l(n)/(5n)
A 0
(2_)( 2)}2 2
sinOA +n sinG = ( 1)(1n )
"+n o Sn
(5.3)
and the volume element is:
3 3 2 2
d r = (A RA)3 (E 2) dE dn d4
Since the two atomic spheres are in general of dif
ferent size, the geometry of the partition is a function of
two parameters, the bond distance R and the ratio of the
radii of any two spheres, in this case, we choose 6 = RB/R.
Then we can express the different quantities in terms of
these two parameters:
RA = (1)/R: RB = BR; Rout = R; RAO = BR;
RBO = (16)/R (5.4)
With the coordinate system as expressed in equation (5.3)
and also by the relations in (5.4), we have the following
restrictions for the intersphere region:
(1) rA > RA
implies ARAO(S+n) 3 RA or c+n > 2(1B)/B
(2) r0 < Rout
implies 2RAO(S+n) S Rout or En < 2/6
(3) rB RB
2 2 2
and r = rA + R 2r Rcose
B A A A
24 242
implies n (1 B)En + E (B +81) ) 0
(4)
(4) E 1
(5) l < n < 1
(i) B>1/2
(ii) 6<1/2
12 2
r
Fig. 5.2 Bounoary of the Intersphere Fegion in Prolate
Spheroidal Coordinates with A and O as the Foci.
Combining (1), (2), (3), (4) and (5), we have the
transformed region as defined by the following limits (figure
5.2):
1
0
2 3
1 2
1
if 2 2
if
if 6 <
where f (S) = +{(136)C {(1B)C + (B2+6)}1}
2/6 2 C if 5 < 2/81
and f (C) =
1
S 2/ if C > 2/81
Expanding the spherical harmonics and the Neumann
functions in terms of 5, q, p as we did in the force calcu
lation in Chapter IV, we have the following result after
some straightforward algebra:
A= N () 1 twi] ((fi~.)]
IA = k N (1) z+k2+1(RAO/2)4 dC .
2m2v+2r rs2
((_l)r Wim W2m m CC,+22m2+2r ( RAs2
vSe r1 2 R)
x 2(1_2 s E1+Z22m2v+2rs {2C(1+2 ) hi' (s+l,)
I1
g +1 h *1 (sR) + C/(KRA)2 h2 (1i}}
RAO 1 2 AO Z012
where v = v + v2
and the h'functions are defined as:
h' (s,) = Z Z 1)+ A ( )
U12 V:o t rO
x Q(u v,1)(s+2n+u+v,S)
~ 1~2
where the Acoefficients are also defined in Chapter IV and the
Qfunctions are defined in terms of sine and cosine integrals,
similarly to the corresponding terms in the force calcula
tions except for the integration limits:
Q(uv,1)(,') = ( cos(z + )r (6 6 )
Zi2 1 2 u+v,O u+v,2
I cost
+ sin(Z1+Z)n 2 u)vT1} 6R ,(1) r dt
+ {sini(1+L,2)7 (Su+v,26u+v,O)
+ cosi( ,+e,) 6 u+v,}1 ) intdt
jKg(+tfirp5) tL
+ {sin(1 2)r 6U+V,l (6vOvl )
)+ cs(2) (uv, u ,) ( )) dt
+ cosl(9LL2)r (6u+v,O u+v,2 1f(5tfp) t
kRllRJf,
(2) IA: J1 nTk(KrA)Y(A) J2(KrO)Y ((0)z0d3r
This is also a twocenter integrals, we can use the
same coordinate system as in the evaluation of IAA. The re
sult is:
tm m
IAA = N NM (1) (RAO/2)4(
1 2 1
z J )u+v m
n,.o a (so 1
SO v=o =
2m 1 22v
V (1)
2
m 2u ,m2v, 9.m2V2 12v 2u+vr+L 2v s
x C C C (KRAO) 1
u v r s AO
x (12 )2,m2vr+2v2s v+r+s {(1+2
x h 2 (Z 2v 2u+vr,922v2s, ) S/(KRAO)
x h 2( ,2v 2u+vrl, 22v,s,()}}
J, tCLa.)] t4')J
where h2 2(3,9,) = '
gso VgO kls
1 20 as voc
u+k.+k,
( (1) 1 +'(Z1+,2k,+u)(2+&,2k+v)
x Q^ Uv,2)(.3+2k,+u,Z,+2k2+v,S) }
and the Qfunctions are defined by:
Q(uv,2)(p,q,) = {sin(KRAAo (Z_2)I) (6 6 )
t2 AO u+v,0 u+v,2
+ cos(KRAOi(i12)r) u+v, }
x JK S cost
P dt
RoI~tf,ct) t(2KRA0,_t)q
+{cos(KRA0( Z1Z )) (u+v,2 6u+v,0)
+ sin(KRAO (ZL12)O) u+v,l)
Sjsint dt
J<^o( ,p) tP(2
+{sin(KRA0+( +k2)7) (6u+v,0+6 u+v,2)
+ cos(KRAO+*( +z 2 )T) 6 u,(6u 6 ul)
1
"4o(ftg1) tP(2KRAOt)q dt
(3) IBB: rliz(rB)Y(rB h(KrB)Y2rB)zO d3r
This is also a twocenter integral. We treat this
integral by considering the Prolate Spheroidal Coordinates
with 0 and B as the two foci:
S= (r0+ rB)/RBO = (r rB)/RBO = 0 = B .
With this coordinate system, we have the following restric
tion for the intersphere region:
(1) rB > RB
implies (RB0/2)(S n) > RB or n 5 28/(1 l )
(2) r0 < Rout
implies (RB0/2)( + n) ( Rout or 5 + n < 2/(1 8)
(3 rA > RA
2 2 2
and rA = rB + R + 2rgR cosOB
2 4 2 4 2
implies n + (1 2) n + + (1 y ) 0
Y Y
where y = 16 = RA/R
(4) 5 3 1
(5) 1 $ n ( 1
Combining all the above restrictions, we have the
transformed region as defined by the following (figure 5.3):
(1) Y<1/2 (B>1/2)
1 
n= +2= ;_=
F t r Ri i\
1 n=g()
aii) y> 1/2 (8<1/2)
nn= 2 =
_/ E 2/y~
1 2 13
/ ._^ ^
rln=gl (E)
I"
Fig. 5.3 Boundary of the Intersphere Region in
Prolate Spheroidal Coordinates with O
and B as the Foci.
2/Y 3. if y < (B> )
i >, =
1 if Y (B F )
2 2
C c 2/y + 1
g (C) n $ g2( )
2 2 (1_ 2))}
where g () = { (1y/2) + ((y)2 (Yy2))
g( 2/y + 2 if 5 < 2/y 1
g (5) =
2/y 5 if 5 > 2/y 1
Expanding the spherical harmonics and the Neumann
functions in , n, ), we have the following expression for
IBB:
BEB < < (1)N+2+1 (RBO/2)4 d"' I
ro , s o
(_ 1)z+2,2m+3r2v 2m 2m
V" 1 V2
xm C sZ+z22m+2r2v 2(Ro ) s2 E,+,22m+2r2vs
xr s 2(2BO
2
x (1s2 ) {2((1+E2) h ) 2(s+1,) 3 +1 h (s,
+ R/ B BO
+ f/(.RB. 2 h 1.2(s,1,W) }
where v = v, + v2
and the h3functions are defined as:
h (s, i) = L(1
2 uo ve a n=O
u+v (u,v) (u,v,3)
An (si+nu+v,C)
n1 2.ij
where the Acoefficients are defined as before and the Q
functions are defined as below:
= ( cos (Z,+k2)I (Su+v, u+v,2)
+ sin&(+Z +E2 uZ v sJt dt
u+v,1 t,
+ {sinj(tZ1+2) (6u+v,2 u+v,O)
+ cKRsoQOdl))
+ cos(Zi+k2)r au+vT jS }
k,( 2.'Oap)
sint
9dt
(4) IBO: f,1
+ icos(z ( U+v,0+6u+v,2
(*Reo~l ,tp) .dt
+ sin(1^ z) 6u+v, (6v6v)} t
1("KrB)Y'(rB)j2(OKro)Y (r0)z0 d3r
We evaluate this kind of integrals by using the same
coordinate system as defined in the treatment of IBB, i.e.,
with 0 and B as the two foci. The result after some straight
forward algebra is:
Q(u,v,3)( )
1i2

