
Citation 
 Permanent Link:
 http://ufdc.ufl.edu/AA00026482/00001
Material Information
 Title:
 Metallic hydrogena twobody approach to band theory and a quantum mechanical test of diophantine methods
 Alternate title:
 Quantum mechanical test of diophantine methods
 Creator:
 Burdick, Stephen A ( Stephen Arthur ), 1952
 Publication Date:
 1980
 Language:
 English
 Physical Description:
 vi, 118 leaves : ill. ; 28 cm.
Subjects
 Subjects / Keywords:
 Coordinate systems ( jstor )
Electrons ( jstor ) Energy ( jstor ) Estimation bias ( jstor ) Estimation methods ( jstor ) Hydrogen ( jstor ) Integrands ( jstor ) Point estimators ( jstor ) Wave functions ( jstor ) Weighting functions ( jstor ) Diophantine analysis ( lcsh ) Dissertations, Academic  Physics  UF Hydrogen  Analysis ( lcsh ) Physics thesis Ph. D
 Genre:
 bibliography ( marcgt )
theses ( marcgt ) nonfiction ( marcgt )
Notes
 Thesis:
 Thesis (Ph. D.)University of Florida, 1980.
 Bibliography:
 Includes bibliographical references (leaves 115117).
 Additional Physical Form:
 Also available online.
 General Note:
 Typescript.
 General Note:
 Vita.
 Statement of Responsibility:
 by Stephen A. Burdick.
Record Information
 Source Institution:
 University of Florida
 Holding Location:
 University of Florida
 Rights Management:
 Copyright Stephen A. Burdick. Permission granted to the University of Florida to digitize, archive and distribute this item for nonprofit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
 Resource Identifier:
 07391544 ( OCLC )
0023447918 ( ALEPH )

Downloads 
This item has the following downloads:

Full Text 
METALLIC HYDROGENA TWOBODY APPROACH TO BAND THEORY
AND
A QUANTUM MECHANICAL TEST OF DIOPHANTINE METHODS
BY
STEPHEN A. BURDICK
A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL
OF THE UNIVERSITY OF FLORIDA IN
PARTIAL FULFILLMENT OF THE REQUIREMENTS
FOR THE DEGREE OF DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
1980
Copyright 1980
by
Stephen A. Burdick
ACKNOWLEDGEMENTS
The author wishes to thank the Kirkland Data Center for the
Division of Highway Safety and Motor Vehicles for the enormous amounts
of computer time this work required. The Northeast Regional Data
Center was extremely helpful in making the connections and transmissions
of data between Gainesville and the D.H.S.M.V. at Tallahassee. The
author is also indebted to Dr. Coldwell for many lively discussions,
his development of the biased selection method and coding of the
minimization routines, and the hours of time he had to put in to get
the required amounts of computer time. The author also had numerous
enlightening discussions with Dr. Broyles and Dr. Burdick. The many
hours of typing were put in by hiswife, Linda. Finally, the author
would like to express his sincerest appreciation to his entire family
for their continual support and understanding throughout this most
challenging period.
iii
TABLE OF CONTENTS
Page
ACKNOWLEDGEMENTS iii
ABSTRACT v
INTRODUCTION 1
PART I: METALLIC HYDROGEN
A TWOBODY APPROACH TO BAND THEORY 2
History 4
The Model Hamiltonian 7
The Wave Function 11
The Monte Carlo Error and Finding
the Eigenfunctions of Equation (4) 15
A TwoElectron Approach to Band Theory 17
Results 32
Conclusion 33
List of References 34
PART II: A QUANTUM MECHANICAL TEST OF DIOPHANTINE METHODS 36
Haselgrove's Method 38
Haselgrove's Method Used by Itself 42
Haselgrove's Method Used in Combination
with Biased Selection 47
Further Tests of the Blend of Haselgrove's
and the Biased Selection Methods 51
Monte Carlo Estimation of the Error 55
List of References 58
CONCLUSIONS OF DISSERTATION 59
APPENDIX A: THE BIASED SELECTION METHOD 61
APPENDIX B: SOME Cm's WHICH FORM SOLUTIONS TO
EQUATION (4) OF PART I 96
ALPHABETIZED BIBLIOGRAPHY 115
BIOGRAPHICAL SKETCH 118
iv
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
METALLIC HYDROGENA TWOBODY APPROACH TO BAND THEORY
AND
A QUANTUM MECHANICAL TEST OF DIOPHANTINE METHODS
By
Stephen A. Burdick
December 1980
Chairman: Arthur A. Broyles
Major Department: Physics
The first part of this work demonstrates the feasibility of cal
culating the properties of bulk hydrogen from a periodically repeated
cell of two "atoms" using biased selection. It further shows how
conventional band theory can be extended to use 2electron wave func
tions. BornOppenheimer potentials were generated for a few 's for
the first three bands and for rs = 1.23aB, 1.48aB, 1.97aB, 2.96aB and
10.84aB ( r 3 mean volume per electron). An insulatorconductor
8B rs3
transition was found to occur between rs = 1.7aB and r = 1.9aB,
while the solid was still molecular. A molecularatomic transition
was predicted, by extrapolation, to occur for rs= 1.1. Finally, a
significant broadening of the kaveraged BornOppenheimer potential
takes place at rs = 1.97aB, where the first signs of an insulator
v
conductor transition occur. (The energy is between 1.17Ry/atom and
1.15Ry/atom when the internuclear separation is between 1.25aB and
2.20aB.) This effect is bound to have thermodynamical effects beyond
those normally associated with insulatorconductor transitions.
The second part of this work details the application of Diophan
tine methods to the evaluation of multidimensional quantummechanical
1
expectation values. It shows that a convergence was found for a
substantial range of N, the number of multidimensional integration
points. For the same accuracy as biased selection, Diophantine
methods required onefourth the number of points, for the remaining
range of N.
vi
INTRODUCTION
This work is most logically divided into two parts. Part I,
Metallic HydrogenA TwoBody Approach to Band Theory, details the
methods and results of a solidstate type of calculation for bulk
hydrogen. Part II, A Quantum Mechanical Test of Diophantine Methods,
details an exploration into the possibility of getting a convergence
for general multidimensional integrals that is faster than the Monte
Carlo convergence of 1/vT. Parts I and II are selfcontained: the
equations and references are numbered as if the other part were not
present, and the understanding of one is not required for the under
standing of the other.
1
PART I:
METALLIC HYDROGEN
A TWOBODY APPROACH TO BAND THEORY
The fact that hydrogen becomes metallic under high pressure is
interesting for several reasons. Firstly, theoretical calculations
exist that have found the superconducting transition to be close to
room temperature [1], [2]. Secondly, the interior of Jupiter contains
hydrogen under high pressure (P > 30Mb) [3]. Thirdly, the conceptual
distinction between an insulator and a conductor is clear cut. If,
experimentally, hydrogen is placed under pressure with a small
voltage applied to attached electrodes, either a current will flow
and it is a conductor, or no current will flow and it is an insula
tor. If, theoretically, the first two energy bands are calculated,
either they overlap and hydrogen is a conductor, or they are separated
and it is an insulator at 00K.
In this work the usual method for calculating the properties of
a crystal was extended, from treating the electrons one at a time to
treating them two at a time. The explicit inclusion of electron
electron distances in the wave function allowed the bulk of the corre
lation energy to be included. An approximate wave function for the
system was constructed from 2electron wave functions in analogy with
conventional band theory. The overlapping of the first two energy
bands indicated that at 00K the molecular solid becomes a conductor
for an rs between 1.7aB and 1.9aB (rs for liquid hydrogen is = 3.4aB,
47r 3 = L3
r3 r= N is the number of electrons contained in the crystal volume,
2
3
L3). The kaveraged BornOppenheimer potential indicated that at 00K
the solid was still molecular for the smallest rs used, 1.23aB.
However, if the data are extrapolated to smaller rs, then the molecular
atomic transition appears to occur at rs = l.1aB. The energies found
here can be used to form the partition function, from which the prop
erties of the crystal can be calculated for temperatures less than
the ionization temperature.
4
History
J. D. Bernal was the first to suggest that insulators would be
come conductors under high pressures [4]. Prior to 1974, most calcu
lations were made assuming that the molecularatomic and insulator
conductor transitions coincided for solid hydrogen. Atomic hydrogen
has one electron per fundamental cell. Since each energy band can
accomodate two electrons per fundamental cell, the valence band is
only half full. This means that an atomichydrogen solid is neces
sarily a metal. Therefore, it is not clear whether assuming the two
transitions coincided was a matter of belief or merely computational
convenience. However, since 1974, theoretical evidence has been
mounting, which this work further supports, that the solid becomes a
conductor while it is still molecular [5], [6]. This makes the
insulatorconductor transition a continuous function of rs, second
order, because there is no abrupt change in symmetry, from molecular
to atomic, to create a discontinuity.
Wigner and Huntington were the first to calculate the properties
of an atomichydrogen solid [4]. They used the cellular, or Wigner
Seitz, method and found an energy minimum of 1.05 Ry/atom at an rs
of 1.63aB. Later calculations essentially reproduced these results.
The solid hydrogen calculations fall into three general categories.
There are calculations to determine the properties of the molecular
solid, generally by summing over predetermined H2H2 potentials [7],
sometimes by using a low density expansion such as Brueckner's method [8].
These methods commonly assume the potentials to be density independent
when they extrapolate the equation of state to the molecularatomic
transition density. This work casts doubts on this assumption near the
molecularatomic transition. The HH potentials were found to depend
on rs, which makes it highly likely that the H2H2 potentials will
also. Furthermore, it is not clear how an H2H2 potential should be
defined when the molecules begin to dissociate. Some of the molecular
solid calculations accurately reproduce the experimentally available
equation of state (0 < P < .025Mb) [7].
Other calculations are on the properties of an atomic solid; 0
function [9]; Wannier function [10]; Green's function [10], augmented
plane wave [11]; exact HartreeFock [12], [13], [14]; perturbation
theory [15], [16]; cellular [4], [17]; alternate molecular orbitals
[18]; and quantum statistics [19]. Ross and McMahan [11] and McMahan
[5] state that the uncertainty in the correlation energy is the out
standing problem for these methods (the energies agree to within
.04 Ry/atom while the uncertainty in the correlation energy is about
.13 Ry/atom). They further state that this uncertainty makes the
molecularatomic transition pressure uncertain by a factor of about
two or three. This work includes the properties of a bodycentered
cubic atomichydrogen solid as a subset.
Finally, there have been a few calculations since 1974 that, like
this work, use the same formalism to calculate the properties of both
the atomic and the molecular forms of solid hydrogen [5], [6]. The
potential and periodicity are treated somewhat more crudely here than
in other methods, but this work uses much better wave functions that
contain electronelectron distances explicitly.
6
The molecularatomic transition pressure is usually calculated from
the slope of the common tangent to the atomicsolid energyversus
density and the molecularsolid energyversusdensity curves. This
pressure is also frequently assumed to be the insulatorconductor tran
sition pressure. This work distinguishes the molecularatomic tran
sition, which at 00K did not occur for any density considered, from
the insulatorconductor transition. The molecularatomic transition
takes place (at 00K) when there is no longer a minimum in the Born
Oppenheimer potential for any nuclearnuclear distance other than the
one corresponding to a bodycentered cubic atomic lattice. Also, the
energy at the minimum of the kaveraged BornOppenheimer potential and
the energy at an internuclear separation corresponding the the body
centered cubic lattice will agree at the rs where the molecularatomic
transition occurs. The data of this work were extrapolated to find
this transition. The insulatorconductor transition takes place when
the bands overlap at the equilibrium separation of the nuclei.
As of 1977, Ruoff felt that there was not any convincing experi
mental evidence for an insulatorconductor transition in hydrogen [20].
However, a more recent experiment claims to have found one [21]. Also,
there are ongoing experiments to produce an insulatorconductor tran
sition in hydrogen [22].
7
The Model Hamiltonian
It is currently impossible to solve the N (~1023) particle
Schrbdinger equation for a general system. A standard method of
dealing with th.is problem is to model the real system with a peri
odically repeated unit cell of fewer particles. In this work, two
particles are placed in the unit cell and the Nparticle Hamiltonian,
HN, becomes
HN(rl ... rN) = H2( + =), (1)
where the coordinates of the ith electron are given by ri, r =
(rx', rly, r Iz' r2x' r2y r2z) and 1 (nlx' nly, n"lz n2x' n2y, n2z)a.
The n's are integers and "a" is the length of a side of the cubical
unit cell. Since HN is then a sum of 2particle Hamiltonians for
the different cells, the Nparticle wave function will be an anti
symmetrized product of 2electron wave functions for the cells.
In the BornOppenheimer approximation, H2(F) is written as
HV2 + + 2 2 + V'(F),
Al B2 AB 12 PBl PA2
where A and B refer to the nuclei, 1 and 2 refer to the electrons,
and V'(F) represents that part of the potential that is not explicitly
written down. The P's are the shortest distances between two par
ticles, which can be in adjacent cells, of a given kind. Figure 1
illustrates this for a periodically repeated 2dimensional unit cell.
8
I B B I B B B B
A 2 A 2 A 2 A 2A 2A 2
I B B I B Bl B B
A 2 A2 A 2 A 2 2A 2
B IB I B Bi B I B
A 2 A 2 A 2 A 2A 2A 2
(a) (b)
Figure 1. A periodically repeated 2dimensional unit cell is pictured.
The interparticle interactions that are considered in this work are
represented by lines that connect the particles. A configuration where
all interparticle distances, P's, are within the unit cell is pictured
in (a). In (b) some P's are between particles in the unit cell while
others are between a particle in the unit cell and a particle in an
adjacent cell. In both cases, the P's are between the closest particles
of a given type.
9
If lines were drawn to connect each particle in the unit cell with all
of those particles to which it is not already connected, then these
connections would represent the same interparticle interactions that
1 1
V'(r) does. A general group of terms would look like  + r
1 1
+ where the primes are for particles not explicitly
rlB, r121
included in H2(F). The 1 in the numerator is because the remain
ing half of the interaction will be included in H2(F) for the cell
which the primed particle occupies. Note that this group of terms is
approximately zero.
In this work, V'(r) is neglected. It could have been included
by using an Ewaldtype sum, but the desirability of including V'(r)
is questionable. For one thing, the imposed periodicity certainly
misrepresents nature. Thus,including V'(r) might actually lessen the
correspondence between the model and the real system. Secondly,
the real system V'(r) is bound to be small due to charge neutrality
and screening, especially for a conductor.
The periodicity of H2(f), i.e., H2() = H2(r + I), was used to
advantage in constructing solutions to
H2(r)O(r) = EO(r). (2)
If D(F) is written as
=(e) = eik'ruk( ) (3)
(this can always be done; e.g.,let uk(r) = eik r(F),then ek'ru(k) =
ik. ik. #
e e (r) = Q(r)),then the translational invariance of H2(r)
demands that uk(r) = uk(( + I). That is to say that uk(r) will be
periodic without approximation. However, uk(r) might only be approxi
10
mately periodic for the real system. The benefits of using (3) will
be discussed later.
Solving (2) for a 2electron wave function instills a certain
amount of selfconsistency in the wave function, without iterating.
HartreeFock type methods try to solve a 1electron equation where
V2 2 2 2 + V'(F) in (2) is replaced by the interaction
SP B2 PA2 P12
of electron 1 with an electron charge density, the coulomb and ex
change potentials. This makes it necessary to solve these equations
iteratively to achieve selfconsistency, since the electron charge
density depends on the 1electron wave functions that are being solved
for.
If (3) is substituted into (2) and V'(r) is approximated by zero,
the following equation is obtained
[72 2 2 (k. + k + 2 2
1 P 2 P 1 2i 2 2) + P
Al 82 AB 12
A2 + (k1 + k 2 E)]ui() = 0 (4)
PA2 PB1 1 2
This equation was solved in the BornOppenheimer approximation for an
"H2molecule" in a cubical unit cell of side "a". The boundary con
ditions were satisfied by restricting the class of functions used to
represent uk(r) (further discussion will be given later). The "molec
ularaxis" was along a body diagonal with its center at the cell's
center. This positioning of the "molecule" allows the crystal lattice
a I
to be varied from a simple cubic molecular lattice, rAB < 32 to a
Aa 
bodycentered cubic atomic lattice, rAB = 32. Two is the smallest
number of "atoms" in the unit cell which will allow correlational
effects to be calculated.
The Wave Function
The imposed periodicity allows an exact representation of the
solutions of (2) to be written as in (3). The kdependence of i(r)
makes it possible to construct a Nelectron wave function from a
product of V's with different k's, as will be further discussed later.
A 2electron wave function allows the major portion of the correlation
energy to be calculated. The relatively simple and yet accurate wave
function of James and Coolidge, which explicitly includes r12 terms,
was used as a basis for constructing uk(r) [23]. This choice meant
that 6dimensional integrals had to be evaluated to find the expec
tation value of H2(F). These integrals were evaluated using the Monte
Carlo method of Appendix A. The accuracy of the wave function allowed
the energy to be determined to within chemical accuracy, .01 Ry/molecule,
by using less than 500,000 functional evaluations. This was accomplished
by making Ho/D relatively constant, which makes the Monte Carlo error
very low.
The periodic boundary conditions can be incorporated in uk(r) as
follows:
U() = 7g(r + I), (5)
where g(r + i) is a 6dimensional generalization of
g(r) = g(r) for r < a
constant for r
= constant for r > a
12
This form of g(r) is not mandatory, but it greatly restricts the sum in
(5) to only include the unit cell. Explicit forms of g(F) will be
detailed shortly. Notice that in the case considered by Wigner and
Seitzan ion, with a spherically symmetric potential, and a single
electron per unit cellthe above reduces to the WignerSeitz boundary
conditions:
uk(r) = constant for r > r
Sk = 0.
ar r=rs
Another important similarity exists between the WignerSeitz method
and this work. There is an implicit electronelectron correlation
that is assumed for the Nelectron wave function. In the WignerSeitz
method, the assumed function is zero if any cell contains more than
one electron and i.s otherwise unity. In this work the assumed wave
function is zero if any unit cell contains more than two electrons
and is otherwise unity. (Notice that this function imposes charge
neutrality on the system on a cellbycell basis.) This function
could be made more flexible by allowing it to vary between zero and
one as the number of particles per unit cell varies. Lastly, the
WignerSeitz method does not consider interactions between particles
in different cells. But for one electron per unit cell, neither does
this work: the distance between particles in adjacent cells is greater
than .
2
13
A twentyfour parameter function, modeled after James and Coolidge's
function, is used to represent u (r):
uk(r) = B(xl)B(yl)B(zl)B(x2)B(y2)B(z2)SAlSA2SB1SB2PJC,
where PJC has two forms,
PG = c + c2(n2 + ) + c3( 1 + 2) + 2c4n 2 + 2c5g +
C+6(Elr E2) + C7nln2(E1 + E2) + c8(E12n + E2l) +
c9(C + E2) + 2ci0u2 + 2ci u,
and
PE = (cl + c4s 12 + 5u + c8U2)(l n2) + c2( 2n2) +
c3( E2 E2n1) + (c6 + cl0u)(e1nq e22) +
(c7 + cllu)(eq2 e2rl) + c9(In n'),
with
u = 2f(P12)
1 f(PAl) + f(PBI)
S= f(PA) f(P1)
f(x)= 12 x) 2(a 2 x) +  x) (a 2] + b}/rAB
and b = 12{ae 2(a' .2)3 + (a' .4)3], a' = Furthermore,
20
SAl = 1 Cm(Km rAl)+ + c24'
m=12
with Km (m 1)MIN(2, 5)
S10 2'
14
(x) = x for x > 0
= 0 for x > 0
(the same c 's and K 's were used for all of the S's). Finally,
m m
2r 4rr
1 x 1
B(x) = c2 + c2e + c23e
The form of f(x) insures that the PJC's satisfy the conditions of (5)
the knots are less than which drives the function to b and its
derivatives to zero for x > The same condition is fulfilled for
the S's, in the same way, except that now the constant is c24. Finally,
the B's are constructed so that they are explicitly periodic. The
reason for two PJC's is discussed in the section on band theory. As
previously mentioned, the 2body form of uk(r) allows the effects of
2body correlations to be calculated. The twentyfour variational
parameters, c 's, are tabulated for the calculated wave functions
in Appendix B.
15
The Monte Carlo Error and Finding
the Eigenfunctions of Equation (4)
The biased selection method (see Appendix A) was used to estimate
the expectation value of the Hamiltonian,
= IEmAm/A2 (6)
where Em [H m) m)/ m), Am = u (Rm)/w(m), and w(%m) is a
function that exactly compensates for the biased manner used to select
rm. Note that Em will vary unless u (r) is an exact solution of
(4). The estimate of the Monte Carlo error made in estimating is
S= [Z(EmA2)2 2EmA' + 2Am]2/A2A
where the above sums extend from m = 1 to m = M, M is the number of
6dimensional Monte Carlo points used. Note that if u&(F) were an
eigenfunction that a becomes [E2 Am 2E2ZAm + E2Am]2/1Am which is
zero. In other words, (6) would have zero Monte Carlo error. This
fact was employed to find the eigenfunctions of (4) by varying the
cm 's to minimize a. After minimizing, u(rF) is almost, but not quite,
an eigenfunction. Therefore,a will be proportional to M 2 with a
small constant of proportionality.
Equation (4) was solved by varying the twentyfour cm's in uk(F)
to minimize c. During the minimization, the in a had to be replaced
by a constant that was less than the eigenenergy of the state being
sought. Two techniques were used to vary the c 's: (1) a routine
16
adapted from Nelder and Mead's simplex algorithm [24], (2) a routine
that used analytically determined first derivatives, approximate
second derivatives, and a matrix inversion with smoothing [25].
The first method starts with a test grid that has one more point
than the number of parametersin this case the grid has twentyfive
points in a 24dimensional spacefor each of which a is evaluated.
If vectors, from the origin to the points, are used to locate the
grid points, then twentyfour of the vectors must be linearly indepen
dent. The routine then generally evaluates a for a point lying along
the vector extending from the point where the largest a was found
through the "center of gravity" of the remaining points. If the a
for this trial point is less than the second highest a (of the current
twentyfive a's), then it is kept and the highest is rejected.
The second routine heads in the direction for which a decreases
the most rapidly. Smoothing is used to make the distance moved towards
the predicted minimum only a percentage of the predicted distance.
This is done both to prevent the divergences frequently encountered
in gradient methods and to make the expansion of a in terms of
derivatives reasonably accurate. Additionally, smoothing allows the
unimportant coefficients to become zero without making the matrix
singular.
17
A TwoElectron Approach To Band Theory
By writing the wave function as e uk(F), full advantage can be
taken of the similarities between this method and conventional band
theory. The importance of k here, as for the i in conventional band
theory, lies in the fact that it is a good quantum number of the system
and can therefore be used to label the different 2electron states.
For 1electron wave functions and conventional band theory, the system's
wave function is found by antisymmetrizing TN = (1)= (2)c (3)c (4)
1 2 3 4
... (N1)>t (N). If it is assumed that e(il e) < (...< E(N )
kN N ) < 2) < <
then the ground state of the system will be given by antisymmetrizing
N = 1 (k) > (2)> (3)1 (4) ... (N1)l)N (N). Similar results hold
S1 k1 k2 2 N kN
2 2
for two electron wave functions: TN = t2(1, 2)k 3* (3, 4) ...
N_ N (N1, N); and assuming that E~ < e 2 E t < E 3 etc., and
kNI N 1 tI 1 kt2 k 2 2 2 3
that a mn can only accommodate two electrons,makes the ground state
km n
the antisymmetrization of TNG (1, 2) (3, 4) ... (N1 N).
2 2
Finally, the properties that im' e( ), and (a) and (b) have in
m m m m
conventional band theory will also be assumed to be the properties of
km (k) (a, b) (km has the components of km as its first three
m
components and as its next, and final, three components also). Thus,
the eigenvalues of (4), E(k ), can be used to find the energy bands.
These bands are continuous functions of k within each Brillouin zone,
m
18
with possible discontinuities at the zone boundaries. Remember that
even though km is 6dimensional, here it consists of two 3dimensional
t's that are exactly analogous to the k's in conventional band theory.
Additionally, a km labels unique electron states only as long as it
is contained within the first Brillouin zone. If k extends beyond
the first zone, then it can always be reduced back to the first zone
by adding a reciprocal lattice vector to it. Therefore, km's outside
the first zone merely relabel a state already labeled by a km within
the first zone.
For k = 0, the center of the first Brillouin zone, (4) has a
spectrum of eigenenergies, em s. These energies can be used to gener
ate the various energy bands, m (k),as k is varied over the first zone.
It is interesting to note that if k were varied from the center of the
first zone to the center of the second Brillouin zone that ui(r), if
it were sufficiently flexible, could vary continuously from the ground
state solution of (4) to either the first excited state or back to the
ground state solution. This was the reason two PJC's were used. Once
the ground state u (r) was found (using PG), then a method was needed
that would allow the first excited state to be found while guaranteeing
that a poor approximation to the ground state was not found. This was
accomplished by using a function, PE, that was orthogonal to the
ground state. Two states were required because the first two bands
were needed to find the insulatorconductor transition.
Since PG and PE are the two eigenstates of lowest energy for the
hydrogen molecule, they were used to generate the first two energy
bands. The bands generated by PG and PE' EG and EE, turned out to
be two approximately parallel parabolas of the form
19
EmE() = m(0) + dmkk. (7)
(This form is supported by conventional band theory calculations to the
desired accuracy [5], [9], [10], [15].) Apparently a band had been
missed that should have roughly connected E(0) with SG(ks),
k e (1, 0, 0, 1, 0, 0). This is what nearly free electron band
s a
theory would lead one to predict, and there are also conventional band
theory calculations that have this state. If the state existed, it
could be found by directly calculating the band gap, 6 at ks from PG
and seeing if it was less than E(ks) G(ks). When it was calculated,
6 turned out to be less than EE(0) EG(ks), and thus the missing state
had been found.
The band gap at k, was calculated by forcing B(xl)B(y1)B(zl)
ik .r ik r
B(x2)B(y2)B(z2), b(r), to be b = e e This creates the
equivalent of the familiar forms for the wave function at the zone
boundary in conventional band theory, cos(2x)x(x) and sin(~x)k(x).
a a
The uk (r) with the b+ would then be orthogonal to the uks (r) with
the b_. This assured that the b form would give s2(k ) and not a
poor approximation to l(ks). The existence of this new energy band
was further verified by varying k slightly away from ks and finding
energies lying on the new curve. The new energy band was taken to be
an inverted parabola. (i.e. d2 < 0 in (7)) connecting SE(0) and the
newly found E2(ks) = EG(ks)+ 5. The first and third bands were found
in a similar fashion: d3 was chosen to connect E(0) with eE(kc),
= (1, 1, 1, 1, 1, 1) and d3 = Ec) E() and dl was chosen
kc kc
to connect EG(0) with EG(kc), dl = G(c) G(0) A few additional
kc kc
20
e(k)'s were calculated and found to fall on the predicted curves.
These are the energy bands required to calculate the properties of
solid hydrogen, almost.
Before continuing, a brief discussion on why 3Ck) was found in
stead of e2(k) is in order. This problem was caused by symmetry
discrepancies between PG, and PE and b(r). The ground state wave
function, PG' has an sorbital type of symmetry about the box's
center while PE has a porbital type. Thus, for PE to be able to
give e2(k ), and not e3(kc), b(r) had to assume a form at kc that
would make PEb (r) have an sorbital type of symmetry about the box's
b E ikc r ik~ir
center, b(r) = e e Furthermore, P, b would have to
have a porbital type of symmetry about the box's center. Unfortunately,
b(r) was not flexible enough to have these forms. Figure 2 illustrates
these symmetry conflicts.
The eigenenergies do not necessarily give the solid's energy
bands. Close inspection reveals difficulties in constructing the
manybody wave function from some 2body wave functions, if Hartree
Fock theory is correct. The difficulties can be most easily visualized
by trying to construct a 4electron wave function from 2electron wave
functions. The HartreeFock representations of PG and PE are equivalent
to, respectively, 0 0 = (1)0(2) and 1 = 0(I) (2) %0(2)cI(1)
(where 00 and )1 are different 1electron states). Notice that if a
4body state is constructed from $0 and 0I that three particles will
occupy the same 1electron state and antisymmetrization will make this
4body wave function identically zero. Given that 0 is occupied, the
next available state that a 4body wave function can be constructed
from is 02 = :(1)1((2). A look at the HartreeFock energies for
21
SC ,S \C/S C C / S C
ds B s B sB
/ C C \S C C \ C.
C S C
C S CS C S / C ,
SC
>A scB >A s B ;A s B
/ C S C / C S C / C \S C
C C s\ / S \ C, C S CI
A s ; A s s A s s
/ 1 /
",,C,/ S ",,C,, S c C"S B/
\C/ S $ C.' C S //C S
/B AB B
/C\S / C 1 S 0 / C C A)
C S., C N, C S. C I, C S', C c
Figure 2. The nodal lines are drawn for the first excited state
wave function, dashed, and for cos(x ), c, and sin('x ), s, in a
a a
periodically repeated 2dimensional unit cell. Note that neither
the c's nor the s's fall on the dashed lines. This symmetry con
flict, due to the inflexibility of b(r), was the reason that the
eigenenergies homed in on E3(k), instead of e2(k). Note that the
extra nodes introduced by either c or s increase the kinetic energy
of the electron.
22
these states will show how E2(0) and e3(k) must be corrected: 0
has energy E00 = 2eO' D1 has energy E01 EO + 1, and D2 has energy
Ell = 2e1. Thus, to find the energies of the next occupiable states
from the unoccupiable states given by PE, the correct e2(0) is given
by 2cE(0) EG(O) and the correct E3(k) is given by 2EE(k) EG(k).
The bottom of the second band, e2(ks), did not need correcting because
both electrons had been placed in an excited state. The resulting
band energies are plotted in Figures 3 8.
As a possibly very important aside, it should be pointed out
that the antisymmetrization of 0O(1, 2), (3, 4) does not give zero for
general 2electron wave functions. Therefore, if HartreeFock is
correct, the nonzero part left after antisymmetrization must be a
bad approximation to the wave function. Thus,the energy of the anti
symmetrized 0o(1, 2)1 (3, 4) will be much larger than E00 + E1ll
However, conventional band theory is known to be incorrect for certain
systems [26], [27]. This breakdown could be due to the neglect of such
states as 00(1, 2)1I(3, 4) which might actually have energies such
as e00+ E01 At any rate, it is rather interesting that the anti
symmetrizing of a product of 2electron states to form manyelectron
states does not generate a "Pauli principle" as it does for 1electron
states.
Now that the energy bands have been found the kaveraged potentials
2L3 >
can be calculated. This is accomplished by placing TSd electrons
in a volume, di, in kspace. The parabolic form of the bands means that
a ksphere should be filled until it makes contact with the sides of the
first Brillouin zone. Then only that part of the ksphere inside the
first zone will continue to be filled until the energy l1(k) is equal
23
2/
0 1 1
S \
'1(0)
0 1 2
rAB (aB)
Figure 3. The top, El(kc), and bottom, El(O), of the first band, the
bottom, c2(ks), of the second band, and the BornOppenheimer potential,
k, versus internuclear separation are plotted, rs = 1.23. The
vertical lines represent the energies calculated by Monte Carlo and
their a's. Since E2(ks) < l(kc) for the entire range of likely
rAB's, the molecularhydrogen solid is definitely a conductor at
this density. Note that even at this density the molecular solid is
favored over the atomichydrogen solid.
24
\j
2 (0)
(0)
2 
t
0 1 2 3
rAB (aB)
Figure 4. The tops, E2(0) and El(kc), and bottoms, E2(ks) and E1(0),
of the first two bands, and the BornOppenheimer potential, k,
versus internuclear separation are plotted, rs = 1.48. The vertical
lines represent the energies calculated by Monte Carlo and their o's.
Since E2(ks) first dips below El(kc) for an rAB that is less than the
rAB at the minimum of k, the molecularhydrogen solid is a conductor
at this density also.
25
2
1 (0)
0
1 c
c\/"
0 1 2 3
rAB (aB)
Figure 5. The tops, E2(0) and e1(kc), and bottoms, e2(ks) and e6(0),
of the first two bands, and the BornOppenheimer potential, k,
versus rAB are plotted, rs = 1.97. The vertical lines represent the
energies calculated by Monte Carlo and their o's. Since e2(k ) does
not dip below el(kc) until rAB is larger than the rAB for the minimum
of k, the molecularhydrogen solid will not become a conductor until
the temperature is raised slightly above 00K.
26
II i i .. I i I i i
2 
1
0
4)
o
) 00
Si e(), r = 2.46
Sl(O), rs = 3.45
1(0), rs = 2.46
2 
0 1 2 3 4 5
rAB (aB)
Figure 6. The top, El(kc), and bottoms, E1(0), of the first band for
two densities, rs = 2.46 and rs = 3.45, versus internuclear separation
are plotted. The vertical lines represent the energies calculated by
Monte Carlo and their a's.
27
C I
\ 
S\
2(0)
= ^EI (kc)
SelO)
E 2(ks)
1____k
2
I I I I I I I I I I I
0 1 2 3 4 5 6
rAB (aB)
Figure 7. The tops, e2(0) and el(kc), and bottoms, E2(ks) and el(0),
of the first two bands and the BornOppenheimer potential, k,
versus internuclear separation are plotted, rs = 2.96. The vertical
lines represent the energies calculated by Monte Carlo and their o's.
Since E2(ks) does not dip below El(kc) until an rAB well beyond the
rAB for the minimum of k, the molecularhydrogen solid will not
be a conductor until a substantial temperature is reached.
28
2
E
O
o
1
2
2 0
I I I I I 1, I I
0 1 2 3 4 5
rAB (aB)
Figure 8. The first two Kolos and Wolniewicz BornOppenheimer
potentials (solid lines) are compared with those of this work.
The dots represent the energies calculated by Monte Carlo and
their o's. As expected, these potentials did not depend on k.
29
.5
E
0
41
1 2 3
Figure 9. Comparison of this work's BornOppenheimer potential versus
density, 3A and M, with the atomic potentials of Bellemans and
De Leener [19], triangles; Wigner and Huntington [4], and Carr [16],
squares; Calais 18], circles. A this work's atomicsolid potential
1 2 3
r5 (aB)
Figure 9. Comparison of this work's BornOppenheimer potential versus
density, EA and eM, with the atomic potentials of Bellemans and
DeLeener[19], triangles; Wigner and Huntington [4], and Carr [16],
squares; Calais [18], circles. 6 this work's atomicsolid potential
for a bodycentered cubic lattice. eM 7 this work's molecularsolid
potential for a simple cubic lattice. If EM is extrapolated to smaller
rs, it will intersect EA at rs = l.1aB; this signals the molecular
atomic transition.
30
2
12 3
0
o
1
2
1 2 3
rs (aB)
Figure 10. The band gap, 6 = e2(k ) El(kc), versus density is plotted.
The insulatorconductor transition occurs when 6 = 0. The arrows
indicate the extreme limits of this transition, rs = 1.35 and rs = 2.
The vertical lines represent the largest possible range of energies
that 6 could have at each density. The information for this figure
came from Figures 3, 4, 5, and 7.
31
to E2(ks). At this time the volume inside the first zone should continue
to be filled in the same manner while the volume in the second zone
should be filled by occupying those parts outside of the first but in
side of the second zone. (Each time a zone boundary is encountered,
the above procedure would be followed.) The process is continued until
a volume in kspace equal to the volume of a Brillouin zone has been
filled. When the procedure is completed, all Nelectrons of the solid
will have been placed in states.
32
Results
When (4) was solved for several k's at each rAB, several rAB s
at each rs, and several rs's for the first three energy bands, the
following results were obtained. The first two BornOppenheimer
curves of Kolos and Wolniewicz [28] were reproduced for H2 with an
rs of lO.8aB, see Figure 8. The top and bottom of the first energy
band were found to agree, as expected. As rs was decreased, the
bottom of the band decreased in energy while the top of the band
increased, for all three bands. From rs = 10.8aB to rs = 1.23aB,
the full range of densities tested, there were no signs of a molecular
atomic transition at 00K. However, if the data are extrapolated to
smaller rs, the transition would be predicted to occur at rs = 1.1,
see Figure 9. The first signs of conductivity, the overlapping of
the first two bands at the equilibrium rAB, appeared at rs = 1.97aB.
By the time rs was 1.23aB, there was no doubt that the solid was a
conductor at 00K, see Figure 10. As seen in Figure 5, a very inter
esting phenomenom occurs at rs = 1.97aB, where the insulatorconductor
transition begins to appear. A significant broadening of the Born
Oppenheimer potential takes place, the kaveraged energy is between
1.17 Ry/atom and 1.15 Ry/atom for rAB between 1.25aB and 2.20aB.
The fact that the nuclei at this, and only this, density are free to
move 27% of their maximum separation possible and only change energy
by 16% of the well depth, is bound to have thermodynamic consequences
beyond those normally associated with insulatorconductor transitions.
33
Conclusion
In this work it was found that the molecularatomic and the
insulatorconductor transitions did not occur at the same density.
The data indicate that the insulatorconductor transition most likely
occurs between rs = 1.7aB and rs = 1.9aB, with extreme limits of
rs = 1.35aBand rs = 2.0aB (c.f. experiment: rs = 1.3aB to 1.4ab
[21]). Extrapolation of the date indicates that the molecularatomic
transition occurs at rs = l.laB. Finally, the distinct widening of
the minimum of the kaveraged BornOppenheimer potential at rs = 1.97
is bound to have thermodynamical consequences beyond those normally
associated with insulatorconductor transitions.
In this work it was also successfully demonstrated that 2elec
tron wave functions could be used in band theory. This work also
opened the door to further possibilities. The BornOppenheimer
potentials can be used to determine the nuclear motion and quantum
solid effects. The data can be used to find the Slater sum, from
which the thermodynamics of bulk hydrogen can be determined. The
kind of "atoms" in the unit cell can be altered so that new systems
could be explored. Even the number of particles per unit cell could
be increased.
34
List Of References
1. N. W. Ashcroft, Phys. Rev. Lett. 21, 1748 (1968).
2. M. D. Whitmore, J. P. Carbotteand R. C. Shukla, Can. J. Phys.
57, 1185 (1979).
3. S. Mitton, "The Cambridge Encyclopaedia of Astronomy." Crown,
New York, 1977.
4. E. Wigner and H. B. Huntington, J. Chem. Phys. 3, 764 (1935).
5. A. K. McMahan, "High Pressure and Low Temperature Physics."
pp. 2141, Plenum, New York, 1977.
6. D. E. Ramaker, L. Kumar, and F. E. Harris, Phys. Rev. Lett. 34,
812 (1975).
7. R. D. Etters, R. Danilowicz, and W. England, Phys. Rev. A 12,
2199 (1975).
8. E. Ostgaard, Z. Physik 252, 95 (1972).
9. H. W. Myron, M. Y. Boon, and F. M. Mueller, Phys. Rev. B 118,
3810 (1978).
10. W. Andreoni, Phys. Rev. B 14, 4247 (1976).
11. M. Ross and A. K. McMahan, Phys. Rev. B 13, 5154 (1976).
12. J. Oddershede, L. Kumar, and H. J. Monkhorst, Int. J. Quantum
Chem. Symp. No. 8, 447 (1974).
13. H. J. Monkhorst and J. Oddershede, Phys. Rev. Lett. 30, 797 (1973).
14. F. E. Harris, L. Kumar, and H. J. Monkhorst, Phys. Rev. B 7, 2850
(1973).
15. T. Schneider, Helv. Phys. Acta. 42, 957 (1969).
16. W. J. Carr, Phys. Rev. 128, 120 (1962).
17. G. A. Neece, F. J. Rogers, and W. G. Hoover, J. Comput. Phys. 7,
621 (1971).
18. J.L. Calais, Arkiv Fysik 29, 255 (1965).
35
19. A. Bellemans and M. De Leener, Phys. Rev. Lett. 6, 603 (1961).
20. A. L. Ruoff, "High Pressure and Low Temperature Physics."
pp. 2141, Plenum, New York, 1977.
21. P. S. Hawke, T. J. Burgess, D. E. Duerre, J. G. Huebel, R. N. Keeler,
H. Klapper, and W. C. Wallace, Phys. Rev. Lett. 41, 994 (1978).
22. Experiments are being conducted at Cornell University and Los
Alamos Scientific Labroratory.
23. H. M. James and A. S. Coolidge, J. Chem. Phys. 1, 825 (1933).
24. J. Kolwalik and M. R. Osborne, "Methods for Unconstrained
Optimization Problems." American Elsevier, New York, 1968
25. R. L. Coldwell and M. A. Pokrant, J. Computational Phys. 20,
365 (1976); the minimization routine was developed by R. L.
Coldwell.
26. W. A. Harrison, "Solid State Theory." Dover, New York, 1979.
27. J. M. Ziman, "Principles of the Theory of Solids." Cambridge
Univ. Press, New York, 1972.
28. W. Kolos and L. Wolniewicz, J. Chem. Phys. 43, 2429 (1965).
36
PART II:
A QUANTUM MECHANICAL TEST OF DIOPHANTINE METHODS
The wave functions in Part I were calculated by varying the
twentyfour parameters to minimize the Monte Carlo error in estimating
the energy,using roughly 1,000 functional evaluations. However, to
reduce the Monte Carlo error in estimating the energy for the final
set of twentyfour parameters, 500,000 functional evaluations were
required. Thus, 499,000 functional evaluations are used not in im
proving the wave function but only in reducing the Monte Carlo error
of the energy estimate to acceptable levels. Obviously, it would
be a big help if the required number of functional evaluations could
be reduced. In this part of the thes4~ a method is presented which
can decrease the required number of functional evaluations by a
factor of four.
In particular, a numerical experiment was conducted to compare
the biased selection method [1] of evaluating the expectation value
of the Hamiltonian for a dihydrogen molecule using the 11term James
andCoolidge wave function, YlC, [2] with a blend of the biased
selection and Haselgrove's method [3]. In its essence, the blend
amounted to using Haselgrove's method of generating numbers in place
of the random number generator used by the biased selection method
(the latter method is a Monte Carlo scheme). The results were mixed.
The "blend" gave an error proportional to either or depending
on whether N was, respectively, less than or greater than 60,000
(N the number of functional evaluations). The "blend" also gave
an error that went from being larger than the biased selection error
(for N < 10,000) to being only half as large for N > 60,000. These
numbers should be compared with the 10,000 points used in finding an
optimum set of Haselgrove's parameters. It is reasonable to expect
that if 100,000 points were used in the optimization of Haselgrove's
parameters that the above results would be found with 60,000 replaced
1
by 600,000. At any rate, a faster rate of convergence than was
found for a large range of accessible N.
38
Haselgrove's Method
Haselgrove's method was investigated to see if it could be used
to give a more rapid rate of convergence in estimating the multi
dimensional integrals encountered in quantum mechanics than the Monte
Carlo methods in current use. His method promises an N (with i > 1)
convergence for integrands which are fourier expandable and whose
derivatives satisfy certain continuity conditions. The method uses
N
SN = f(o) + 2 Y cmf(2Lmall, ..., 2LmaKI) (1)
m=1
to estimate the integral (2L)K ... ff(x)dx, call it I, when f(x) is
an even function (K = the dimensionality, f(x) E f(xl, ..., xK),
dx = dx1 ... dxK). The cm's (which are functions of N, m, and i)
determine the rate of convergence, the value of i The a.'s are
irrational numbers, which must be found and which admit no solutions
K
to 1 + M M.a. = 0 when the M.'s are integers.
i=l K
For i equal to two the integrand W(x) E 7 (1 Ixi )2, where
K1 1 Ki=l
2 .../ W(x)dx = 3 gives the worst possible error, w (to within
1 1
a multiplicative constant that depends on the specific integral being
evaluated but not on N) of all the integrands satisfying Haselgrove's
conditions. The key to finding the ai's lies in minimizing E This
wasdone by selecting K random numbers (between 0 and 1) and varying
them by a simplex method [4] to minimize Tables I and II give
some a's and corresponding e 's as functions of N.
39
Table I. This table contains the aKM found by minimizing w over
Kdimensions, using either my method (second subscript M) of finding
{2ma } or Haselgrove's method (second subscript H).
a optimized over 1,000 points a optimized over 10,000 points
a5M a6M
.45260537885008 .43657287149021
.47561868232093 .59181442054365
.57787664217299 .050243216161086
.49493808300031 .84377789758502
.21929292303371 .38121353509683
.78500128467888
a optimized over 10,000 points
a8M
a5M
.73778770979932
.45241479419455 .083137623602196
.47565271662758 .84742517966315
.5782437781132 .88984806744176
.49509111688911 .80256318005657
.21929571968924 .27985939980508
.67436009334498
"5H .5303179412788
.957286852634
.86730302040274
.097247969465186
.31301711424498
.48486758446139
40
Table II. This table gives the absolute value of the difference,
multiplied by 10 between W and Haselgrove's estimate of W, using
my method of finding (2ma I (EWM) and Haselgrove's (EWH).
N GWH EWM
(a5M optimized over 1,000 points)
1000 9774 5570
3000 2338 1876
7000 887 596
15000 237 249
31000 125 118
59000 68 47
(a5M optimized over 10,000 points)
1000 7818 7404
3000 3173 2225
7000 1289 541
15000 299 191
31000 112 72
59000 42 35
(a5H optimized over 10,000 points)
1000 5980 5434
3000 1326 1179
7000 405 447
15000 143 223
31000 63 95
59000 38 43
(a from Haselgrove's paper)
1000 5933 5716
3000 1592 1808
7000 455 1086
15000 282 732
31000 178 462
59000 114 241
It can be seen that eWH and eWM have been substantially lowered from
those given by Haselgrove's original a. Also note that eWH and EWM
differ by less as N is increased, regardless of which method was used
to find {2mai}. Furthermore, the EW's for $5H and I5M, optimized over
the same number of points, differ by less as N is increased.
41
N+1 Iml 1
If c is equal to 'N+I or 2 then an i equal to 2 or 1
m (N+l 2N+l
should be found. The main reason an i equal to 1, instead of 2,
might be used is because Haselgrove's method places less stringent
demands on the integrand for smaller i's. In the experiment, both
i's were found to give the same convergence.
Eventually, 2mai will fall outside the integration region as m
is increased. Two distinct methods were used to reduce 2ma. to the
ma.
integration region: Haselgrove's used 2[i2 for each x. ([x] is the
signed distance to the nearest integer, e.g., [1.9] = .1 and [2.1] =
+.l), my way used the periodic equivalent of 2mai lying within the
interval from 1 to +1. (In either case, denote this reduced
number by {2ma i.) The a that is found by varying an initial a to
minimize w will depend on the method used to find {2mai.. If one
method is used to find {2mail for the minimization and the other
method is used to find {2mai} in estimating I, then reasonable,
although not optimal results, will be found. In the experiment,
both methods gave equivalent error estimates.
The ai's used in the experiment had a 14digit computer repre
sentation. Since the ai's are represented by a computer, they are
not irrational. This fact will be a problem when N is roughly 107 [5].
42
Haselgrove's Method Used by Itself
Haselgrove's method was used to estimate the known values of
some definite integrals. His a was used to estimate the value of the
1 1 x1x X3x4x5
5dimensional integral f ... f e 1 1 dx dx2dx34dx5, and his
0 0
results were reproduced exactly.
The method was then tested on a problem that required a change
of variables. The expectation value of the Hamiltonian for a helium
zs 27
atom using T = e", z = 7, s rI + r2, r nuclearelectron
separationwas estimated. This expectation value can be determined
by evaluating
c0 s u
L E27 2dsfdufdt(8su s2 + t2)e2zs
0 0 0
Ss u
M = 2idsfdufdt eZu(s2 t2)[ + ] + 4usT}, (2)
0 00 0ss u s
and
CO u
Norm 2r2dsfdufdt u(s2 t)e2zs
0 0 0
where T a [6]. Haselgrove's method, since it is designed to sample
the 3dimensional space uniformly, will sample the important region
of the integrand with a very low probability. In order to increase
the probability of sampling the important regions,a change of
variables, which amounts to importance sampling [5], can be employed
s' = Zn(l s), t' = u = u
u S
43
These transformations alter Eqs. (2) to
1 1 12z1
L = 2n2fds/dtfdu{u(l s) [8u 1 + u2t2] kn4(1 s)}
0 0 0
M = 22z2fdsdtdu{(l s)2z1 2[1 u2t2+ 4 1 s)} (3)
0 zzn(ls)]
12z
Norm = 22/fdsdtdu{u2(l x)2Z[ u2t2]n5s(1 s)}.
0
In Eqs. (2) and (3), M, L, and Norm are the expectation of V2, the
potential, and 1, respectively. Table III gives some results.
This problem has an interesting feature: the variable transfor
mation introduces singularities into the first derivative of the
integrand. Even with this violation of Haselgrove's conditions,
the results in Table III seem to give a better convergence than 1
The results with helium were good enough to warrant moving on
to a problem which is more difficultthe evaluation of the expec
tation value of the Hamiltonian for the dihydrogen molecule, using 'YJC
First, Haselgrove's method alone was used, with the results given
in Table IV. The lack of any signs of a 2 or even a N convergence
was not unexpected for the N used (the problem, in helium, of samp
ling the important region frequently enough is made more severe here
due to the increased dimensionality). Next,Haselgrove's and the
biased selection methods were combined with the hope of gaining the
advantages of both. Biased selection would be a way of sampling the
important regions more frequently, taking the place of the importance
sampling used for helium. By replacing the random number generator
with Haselgrove's method of generating points in a multidimensional
space, we hoped to gain the advantages of sampling the important
44
Table III. Energy bounds, /Hydr~dr2//T2dr1dr2, for a helium atom
were calculated using Haselgrove's & and two different methods.
N Haselgrove's Method 4EEW My Method 4EEW
1000 5.697346 .002644 5.700216 .003694
3000 5200 611 5.695513 570
5000 026 319 516 228
7000 189 182 350 137
9000 358 114 354 91
11000 338 68 29 68
13000 305 68 15 46
15000 278 68 23 46
17000 285 46 27 46
19000 298 46 15 23
For comparison the exact value of the integrand, which is not the
exact ground state energy, is 5.695313. All answers are in
Rydbergs. The error estimates, 4EeW, were calculated by assuming
that the error in each integral in Eq. (3) was E, and then by
taking small W limit.
45
Table IV. This table gives estimates of the expectation value of the
hydrogen molecule, using James and Coolidge's 11parameter wavefunction
and Haselgrove's method without the aid of biased selection (all
energies are in Rydbergs and the correct answer is 02.34609).
N H H
1000 2.31939 2.49532
3000 2.22855 2.22018
7000 2.22166 2.23877
15000 2.23085 2.24075
31000 2.23844 2.25364
59000 2.25539 2.23981
The second column gives the results when Haselgrove's method was used
to find {2mai}. The third column gives the results when my method
was used to find {2mai.. Because it can be difficult to sample the
important regions of a 6dimensional integrand, the poor results were
not unexpected. Compare these results with the results of Table V
where the help of biased selection was employed. The a8M from
Table I was used.
46
regions with more uniformity [5], [7]. In other words, an error
1 1
estimate which was proportional to rather than seemed possible.
N vI
47
Haselgrove's Method Used in Combination
with Biased Selection
>. 2 >2
The evaluation of E = f ic HY jdx/fTcdx was carried out over an
integration region which was a box of length Lx = 20a (Ly = Lz = 22a,
aB E Bohr radius) in the x (y, z) direction. This box is large enough
to insure the desired accuracy. The two nuclei are fixed on the zaxis
with the point halfway between them coincident with the center of the
box.
If the biased selection method were going to be used to estimate
E, then the following could be done. Select the coordinates of the
first electron in four ways and those of the second in five ways (four
ways being the same as for the first electron). The estimator for
E would be
iN
BN = N1if(,)/w(xi), (4)
where the relative probability of selecting a 6dimensional xi, W(xi),
exactly corrects for the nonuniform selection of the x.'s. If each of
1
the twenty ways (four ways for the first electron times five ways for
the second) of selecting an xi is used with equal probability then
W(xli, x2i) = {wl(xl) + 2(xli) w3(xli) + w4(xli)] x
W(2i) + w2(2i) + w3(2i) + w4(x2i) + 5(2i)]
These wi's, etc., are determined by the methods used to select the
.48
electrons' coordinates:
x = R1 L, y = R2Ly z = R3L
wl(x) = 1
(Ri random number uniformly selected between 0 and 1, a total of six
R.'s are needed for each x. in Eq. (4))
r = /ciMc, = 2R2 1, 2 = 2TR3
w2(x1) = V/(2 Mcrc)
(rc = the electron's distance from the center, p = the cosine of the
angle between rc and the yaxis, E= the angle between the projection
of rc in the xy plane and the xaxis, V = LxLLz' Mc is a constant)
Va = 2R2 1, = 2nR3
rai is selected by requiring
rai 8
f g(y)dy//g(y)dy = R
0 0
w3(x) = Vg(rAi)
8
4Tri/fg(y)dy
0
Vg(rB)
w4(xi) = V Bi
4Tr'ifg(y)dy
0
(a = A or B, i = 1 or 2, rai the distance between nucleus a and
electron i)
49
rE = 4ME' = 2R5l, 1 = 2TR6
w5(x2) = V/(2rMRE)
(rE 5 the electronelectron distance, ME e a constant). The function
g(r) was created by evaluating r2e2r at 64 points and connecting
these values with line segments.
The estimator that is used to estimate E0 is
20 N
I J JC(xim)HYJC(X m)/Wi(Xm)
c i=l m=l
N 20 N (5)
STY JC(xim)/ i(Xm)
i=l m=l
where w() 1wl(rl) + w2(r) + w3(rl) + w4(r)] x
{wl(r2) + w2(r2) + w3(r2) + w4(r2) + W5(r2)]
(ri electron's coordinates, x represents an rI, r2 pair). The
i in Eq. (5) signifies which of the twenty ways (four ways for the
first electron times five ways for the second) is to be used to gen
erate x. from {2m}, for a given m. (This same {2mO} is used to
generate a total of twenty im's.) This procedure is used (instead
of using an 8dimensional a and using the two extra a.'s to randomly
select which of the four (five) ways would be used to select the
first (second) electron's coordinates, as can be done in the usual
biased selection) to insure that im will be close to .in if {2ma
is close to {2na}. Note that for cN to be fairly compared with a
Monte Carlo estimate that the Monte Carlo estimate would have to use
20N functional evaluations.
50
Table V. Energy bounds (in Rydbergs) for an H2 molecule calculated
using the 11term JamesandCoolidge wave function and using a8M
from Table I to generate the "random numbers" that were used to
determine the electron positions.
N Energy Bound _4E_
1000 2.34684 .006982
3000 560 1736
7000 585 619
15000 596 241
31000 601 105
57000 607 57
The blend of Haselgrove's and the biased selection method was used.
The exact answer, to six digits, for the given wave function, is
2.34609 Ryd [2].
51
Further Tests of the Blend of HaselGrove's
and the Biased Selection Methods
One must be careful to insure that the method of generating the
xi s is not creating correlations between the xi's. A commonly used
random number generator, furnished by IBM with their FORTRAN compiler,
was found to have this defect [8]; the random number generator used
for the calculation herein did not. Therefore, Haselgrove's method was
used to calculate the nuclearelectron and, especially, the electron
electron correlation functions. These functions were then compared
with the correlation functions generated by using biased selection.
The functions found by the two methods agreed, indicating that the
configurations generated from a did not have any correlations that
were not due to particleparticle interactions. The correlation
functions can be compared in Figures 1 and 2.
Unwanted correlations were also tested for by seeing if the con
figurations generated by using Haselgrove's scheme, without the help
of biased selection, were uniformly distributed throughout the
volume. This was done by dividing the 6dimensional integration
volume into 729 equalvolume hypercubesi.e., three equally spaced
divisions were made in each dimensionand noting the distribution
of points within them. Figure 3 illustrates that the number of points
per hypercube was normally distributed and that this distribution has
less deviation than uniformly and randomly distributed points have, as
expected.
52
1/
.0 4
.02 /
0 1 2 3 4 5
ree (aB)
Figure 1. Comparison of the biased selection method (I) and the
blend of biased selection and Haselgrove's method (.) for calcu
lating the electronelectron correlation function. The biased se
lection method used 40,000 points and the 'blend' used the equivalent
of 200,000 Monte Carlo points. The excellent agreement demonstrates
that the points generated by the 'blend' did not introduce any non
physical electronelectron correlations.
53
.050 i \
/ \
.025 /
/ \1*
/ 5
0 1 2 3 4 5
rne (aB)
Figure 2. Comparison of the biased selection method (1) and the
blend of biased selection and Haselgrove's method (*) for calcu
lating the nuclearelectron correlation function. The biased se
lection method used 40,000 points and the 'blend' used the equivalent
of 200,000 Monte Carlo points. The excellent agreement demonstrates
that the points generated by the 'blend' did not introduce any non
physical nuclearelectron correlations.
54
&
.4 
/ A&
.2
U &
3a 20 a 0 a 2c 3a
Figure 3. N points were placed in a 6dimensional volume that was
divided into 729 equal cells, using Haselgrove's method. The number
of points falling in each cell was counted, and the mean and standard
deviation of these 729 numbers were calculated. The fraction of cells
containing a number of points which was a given number of standard
deviations away from the mean was plotted. The squares are the re
sults of a 10,000 point run: mean = 13.7, a = 2.53. The triangles
are the results of a 20,000 point run: mean = 27.4, a = 3.73. A
normal curve is drawn for comparison. Also for comparison, a purely
random set of 10,000 (20,000) points would be expected to have a
mean of 13.7 (27.4) and a of 3.70 (5.23).
55
Monte Carlo Estimation of the Error
Since w at best, can only determine an error estimate for Eqs.
(1), (4), and (5) to within an undetermined multiplicative constant,
which is independent of N, some trialanderror searching would be re
quired to find this constant. But there is a very real danger in this
type of an approach. If the constant, B, is large enough to ensure
that is larger than the difference between cN and the correct answer
for some NB, then the error might be overestimated for smaller N and,
at the same time, be underestimated for larger N. Inthe case of esti
mating E0, the 5decimal place accuracy quoted by James and Coolidge
was approximately achieved for N = 59,000. Therefore, the B that would
be correct here could not be tested by increasing N. In order to obtain
a reliable error estimate, a Monte Carlo procedure was used.
In place of Eq. (1), the following estimator
N
SN (i) E mf(2mal + Rli, .. 2maK + RKiI)
m=l
can be used without violating any of Haselgrove's conditions ( si is
a set of K random numbers selected uniformly from the interval between
0 and 1). By selecting M r.'s and estimate can be made of the average
and the standard deviation of the distribution of SN( i)'s. The
standard deviation can be used as an error estimate for any given SN,
SM M
i 1=N )2
i=l J=l
56
An estimate of the error in estimating e is
To prevent the overhead in finding E and 6 from going to waste, the
following refined estimates of I and its error can be used:
i=1
with the error in the above estimate given by
Figure 4 displays some results.
It is interesting to note that for N less than or on the order
of 10,000 the number of configurations used to optimize a that
the convergence appears to be ,but for N very much larger than 10,000
the convergence becomes This suggests the possibility of extending
the 1 convergence to larger N by increasing the number of points used
to optimize a.
57
8000
4000 "
2000 
1000 
400
200 
"1
100 
40
+ 20 
S20
x 10 III I
0 .1 .2 .4 1 2 4 10 20 40 100
N 104
Figure 4. Comparison of the error using biased selection (straight
line prediction from 40,000 points) and the error using the blend
of biased selection and Haselgrove's method (vertical lines). The
size of the vertical lines represents the possible error, 6, in esti
mating the error, c.
58
List of References
1. see Appendix A.
2. H. M. James and A. S. Coolidge, J. Chem. Phys. 1, 825 (1933).
3. C. B. Haselgrove, Math of Computations 15, 323 (1961).
4. J. Kolwalik and M. R. Osborne, "Methods for Unconstrained
Optimization Problems." American Elsevier, New York, 1968.
5. J. M. Hammersley and D. C. Handscomb, "Monte Carlo Methods."
Chapman and Hall, London, 1979.
6. B. L. Moiseiwitsch, "Variational Principles." pp.172183,
Vol. XX, Interscience Pub., 1966.
7. H. Niederreiter, "Diophantine Approximation and Its Applications."
pp. 132194, Academic Press, New York and London, 1973.
8. R. L. Coldwell, J. Computational Phys. 14, 223 (1974).
CONCLUSIONS OF DISSERTATION
The calculations of Part I required an enormous amount of computer
time, roughly 1200 hours of Amdahl 470 cpu time. While experience
would allow these calculations to be redone with, perhaps, onehalf to
onethird of the computer time, an enormous amount of computer time
would still be needed. This time was mainly used in the minimization
of a to find the wave functions. Less than 5% of the time was used
to find the energy of the final wave functions for the equilibrium
separation with a Monte Carlo error that was less than chemical
accuracy.
For the required amount of computer time to be reduced by more
than the factor of two to three, three things would have to be done,
separately or jointly. Firstly, the minimization would have to be
made faster. Frequently, more than 500 calls had to be made to the
simplex algorithm to find the best wavefunction. Dr. Coldwell is
making continual progress along this line. Secondly, the number of
Monte Carlo points required to evaluate the integrals accurately
enough to allow meaningful minimizations would have to be reduced.
Part II was a first step and suggested further steps to take along
this line. Thirdly, better wave functions would have to be used.
The form of b(P) was not flexible enough to allow the wave function
for the first energy band to have the same symmetry as the wave
function for the second band: b(r) should have been Tcge (G
reciprocal lattice translation vector). Furthermore, the form of
59
60
ui(r) became more inappropriate as rs 0, more and more Monte Carlo
points were required to minimize a effectively.
This work showed that it is possible to use Monte Carlo and
2electron wave functions to calculate the properties of bulk matter.
It also gave encouraging signs, in Part II, that the 1/Al convergence
of Monte Carlo can be beat in the evaluation of multidimensional
quantummechanical expectation values.
APPENDIX A
THE BIASED SELECTION METHOD
INTRODUCTION
This appendix presents in one place for the first time detailed
proofs, statements, and examples of the biased selection method of
evaluating integrals and expectation values. Biased selection is,
as the name suggests, a method of nonuniformly sampling an integrand,
f(x), in order to estimate ff(x)dx. It is an extension of importance
sampling, as given by Hammersley and Handscomb in reference [1], in
the sense that it uses sums instead of integrals in order to exactly
correct for the bias introduced by using a weight function, g(x).
In addition, biased selection allows one to use several simple weight
functions simultaneously, one for each of several important regions
of the integrand, in place of a single complicated weight function.
Biased selection is harder to code than the Metropolis scheme, [1],
[2], [3], but more flexible. Coldwell and Lowther used it in their
calculations on the dihelium system [4].
A weight function, g(x) > 0, that mimics f(x) can be used to
lower the standard deviation in the Monte Carlo estimate of ff(x)dx.
For example, let
y = /fg(z)dz/fg(z)dz
so that
61
62
ff(x)dx = ff(x(y)) dy
w(x(y))
(all integrals, from now on, are from 0 to 1, unless explicitly stated
otherwise) where
w(x) g(x)/fg(z)dz.
The estimate of the integral is then 1 f(x(yi)) where the y's are
w(x(yi))
picked uniformly at random. When the exact value of fg(x)dx is used
the above is importance sampling. It will be shown that biased
selection gives similar results but does not require the exact
evaluation of any integrals in order to exactly correct for the bias
introduced by the nonuniform selection of the x's. This appendix
will also show how several different w's can be averaged and used in
place ofa single w. This averaging allows the combination of several
simple weight functions, rather than one complicated weight function,
to sample the space adequately.
The methods detailed in this appendix often select one point
from M points and repeat this N times to get the integral estimates.
As Mx the biased selection method becomes importance sampling,
owing to the fact that certain sums become integrals. It must be
emphasized, and will be proved, that for M fixed, e.g., for M = 3
the method converges on the correct answer almost surely as Nx.
1 Unless specifically indicated, all sums are from 1 to N.
63
I. ERROR ESTIMATION
The simple Monte Carlo estimate for a multidimensional integral
is
IN = f(i) (1)
where the xi s are chosen uniformly at random. As N>, IN * f(x)dX.
Since the x.'s are random, sampling theory can be used to derive (for
N large [5], [6]) an error estimate for IN
'{ 1 1 (2)
a = (f(xi) IN)2)T~' (2)
and an error estimate for Na2
S= N2o4}7 / T, (3)
where = (f() IN) [1], [2]. For these errors to assume their
standard meanings (e.g. 'IN ff(x)dxl
time) then IN and a2 must be approximately normally distributed about
the values that they would take on if an infinite number of points
were sampled, in this case 6 = /2Na2[I], [7], [8]. For simple Monte
Carlo it is generally claimed that N should be greater than about
30 for the IN'S and greater than about 100 for the o2's to become
almost normally distributed. For both biased selection and importance
sampling, Eqs. (1), (2), and (3) are correct after the following
changes: replace f(t) by f(3)/w(2) and select the x's according
to w() (in importance sampling this can be done by selecting y uni
64
formly at random and then calculating x, see introduction, in the case
of biased selection the techniques of section III can be used). Biased
selection requires fewer points than simple Monte Carlo does, in
order to obtain accurate estimates of the integral and its error. It
reduces a by reducing its proportionality constant, not its inherent
1//T' dependence.
Restrictions must be placed on the integrand for a and 6 to be
used as error estimates and maintain their common meanings. In partic
ular, a and should exist and be bounded,and a should tend to zero
4
as N .
It is suprisingly easy to violate these conditions. Consider
I1
the integral /x 2dx, it exists but a and p are infinite. Numerically,
4
this shows up as a a that decreases as N increases,and then rises when
a random point is chosen which comes much closer to the origin than
all previous points; and as a 6 which becomes much larger than a as
1 2
N is increased. If the proper weight function, e.g. g(x) = x3, is
chosen, the integral and a can be interpreted in the usual way. In
this case, the estimator i.s
1 3 1
N (xi7)3xiT = 7 xiT
and a, /N,, and 6,/, behave properly.
Suppose that nothing more were known about the integrand than
that the largest contributions come from the origin. Observe what
happens if the weight function favors the origin too much. For example,
1 99 1 1 99
if g(x) = looxo0, the integral becomes /x y (100 x(y)1OO)ldy
which reduces to 100l"ydy. The problem at x = 0 has become a problem
at x = 1, y = 1, which, for good numerical results, still requires N
65
to be very large (this is now necessary in order to sample the region
near x = 1 adequately, due to the low probability of placing a point
there).
For a more useful example, consider evaluating the expectation
value of a quantum mechanical Hamiltonian, H. Let E0 be the
eigenvalue of H, and Tt(x) be an approximate (normalized) eigen
function. Then HYt(x) = E(x)Yt(x). Note that if Tt(x) were an exact
eigenfunction that E(x) would be constant. However, for our purposes,
suppose that Tt(x) is not an exact eigenfunction. We will assume
that it is very good for Yt(x) large and very bad for PT(x) small.
If we use importance sampling and let the sampling function be
g(x) = T2(x) + 6 then the integral
2(X(y)) +
S(f C(z) + 6)dz)
can be simplified and the estimator written as
= ( + 6) E(xi) EO (5)
1 + 65/y(x.)
t ,
It might appear that the bestsampling function would be YT(x) (i.e.,
6 = 0), and, in fact, the Metropolis scheme uses this weight function.
But, since most of the error in HYt(x) is assumed to be for small
values of Pt(x), say for Y2(x) less than a tenth of its maximum value
(2 M), then IE(x) EO will be the largest there. The error term
will then be dominated by an area that is unlikely to be sampled (because
Y2(x) is small). This means that a large N is needed to lower the
error in because most of the points contribute little to
the error while a very few make large contributions, since w(;i) is
very small for them. If 6 were kept finite (say 6 = .lw2M) instead
of being set to zero, then the "misses" would be lowered in size to
66
E(x) E0 instead of E(x) E0. This is, of course, exactly compen
1 + 6/ x)
sated for by the fact that there are now more of these "misses".
Notice, in this case, that the term contributing to the error in
is approximately halved. Each "miss" is cut by a factor of
4 in its contribution to a, but there are twice as many of them. So,
by sampling the "unimportant" regions more (i.e., no region is sampled
with a probability that is less than, in this case, .1T2 ),a can be
reduced. It is interesting to observe that a was reduced by using a
seemingly worse weight function, and that a large number of points
appear to be wasted because (E(x) EO)/(l + 6/(Y)) is (as T*0)
zero for the most part.
67
II. BIASED SELECTION
In biased selection, a positive definite weight function, g(x),
is used to distribute the xi's, with a probability PB(i) (or proba
bility density, as the case may be) which must be computable
(explicit schemes for doing this are detailed in Section III). The
bias is compensated for by dividing each f(xi) in Eq. (1) by the
relative probability of choosing xi,
w(x i) PB(x i)/PU(x) (6)
PU(x ) is the unbiased (uniform) probability of choosing xi. One
should note, as the discussion develops, that w(xi) can be a multi
valued function. This is because there can be more than one way to
select any given x. and each way can have a unique w(xi) associated
with it. The biased selection estimator is
1 1 fyx )
IN = f( i) (7)
To prove that Eq. (7) gives If(x)d? as N*o, the following can be done.
Divide the integration region into T equalvolume hypercubes (label
them Ta where a = 1, 2, ..., T) and select N x.'S, each with relative
1
probability w(x'), from it, N>>T [4]. Now, break up the sum into
two sums one over the Ta's and the second over the points within a
given Ta 
IT na
IN I f(I ).
a=1 i=l
68
To derive Eq. (8), two cases need to be considered. The case for the
multivalued w(ix)'s used herein is discussed in Section III. The
case for singlevalued w(xi)'s follows. As T becomes infinite, the
volume of a given Ta becomes vanishingly small. Therefore, the
continuity of f(x) and w(x) requires that j for all
w( i) w( x) for all i
and xj within a given Ta'1 Therefore, combining the na equivalent
terms within each Ta gives
T +
I __l n f(xa) (8)
TN T = W(a )
a=1 a
Since N>>T, the law of large numbers can be used to show that a PB(a,)
SN
almost surely, a.s., as No. Also note that PU(xa) = These two
facts and Eq. (6) imply n N(T)w(Xa); which when substituted into
Eq. (8) gives
IN 1 N( )w( xa) (almost surely)
a=1 w(Xa)
T
= f(xa)"
a= 1
This is the definition of the Riemann integral in the limit of infinite
IN T ff()d.
Note that if w(x) were equal to unity that the above proof reduces
to a proof for simple Monte Carlo.
1 The results remain valid if f()) and g(x) have a countably infinite
set of regions over which they are continuous, but the explanations
are too detailed for this paper's purpose.
69
The closer the weight function can be made to f(x), the more
constant f(x)/w(x) will be; and, consequently, the smaller a will
be (approaching zero for nonnegative f(x) and f(x)/w(x) absolutely
constant). While w(x) = If(x)//f(x)Idx does minimize a, to evaluate
w(x) one would have to know the integral that, by hypothesis, is not
known. It should be cautioned that a bad weight function can greatly
increase the number of points before Eqs. (2) and (3) accurately
indicate how far IN is from ff(x)dx. For example, if w concentrates
the points in "unimportant" regions, the integral effectively becomes
the integral of f(x.) over these "unimportant" regions.. The important
regions will not be sampled for small N because of their small
probability of being selected. Thus, the averages and error estimates
are only for a part of the total integral being sought. This results
in a low a which could lead to the erroneous conclusion that IN was
close to ff(x)dx. instead of the correct conclusion that IN was a
good estimate of the integral of f over the "unimportant" regions.
This misconception will not be revealed until N is large enough for
the "important" regions to finally be sampled. (Note: For the
important regions to eventually be sampled, in general, g(") must
never be zero.)
Let us assume that we have an integrand, f(x), with the locations
of some peaks known and others unknown. Then a measure of the number
of points making a significant contribution to the estimate of the
integral is
Neff f(i)2 f~i)2. (9)
w(xi) w(xi)
70
The number of "effective" points can range from Neff 1, when one
value of f(x.)/w(x.) is much larger than the others (and therefore
makes the major contribution to Eq. (7)), to Neff = N, when all N
values of f(xi)/w(xi) are equal. Since it has been assumed that we
do not know the location of all the peaks, it is undesirable to have
Neff = N. This is because such a large Neff can only be obtained by
concentrating the points where peaks of the integrand were thought to
lie. This concentration precludes any chance of selecting a point
from one of the unknown peaks. A good weight function to use is
g(x) = b(x) + d, where b(x) incorporates, as much as possible, all of
the known properties of the integrand and d is a constant. This way,
if there exists some x0 where f(x0) is large but was thought to be
small (i.e. b(x0) is small), then x0 still has a chance of being selected.
The selection of this x0 causes Neff to drop to approximately unity.
The statistics of this run become highly suspect. However, the precip
itous drop in Neff tells us that we didn't know all of the regions
in which f(x) was large, and therefore, allows us to select a better
b("), enabling us to eventually find the correct value of ff(x)dx.
Notice that it was the presumably wasted points that enabled us to
locate the unsuspected peak.
In practice we have found that when Neff < .2N it pays to improve
b(x), and that for N eff > .3N it pays to sample more of the regions
where b(1) is small (increase d).
The following example illustrates this. Suppose that f(x) has
two equal peaks (where f(x) = 106) that each extend over a length equal
to .001 (where the total integrand is from x = 0 to x = 1) and is one
elsewhere. Then consider a b(x) that equals f at, one peak and zero
71
elsewhere. If d = l,then information sampling would have roughly a
chance in 1000 of sampling where b(x) = 0,and then only a chance in
1000 of that point hitting the second peak. Therefore, for N = 104
there is a 99% chance of not sampling the second peak. This results
in an IN = If(x)dx and an Nef 104. If N is increased to 106.
N 2 eff
there is a 73% chance of hitting the second peak and getting IN
ff(x)dx and an Neff = 4. Finally, if d is increased to 5000,then
there is approximately an 83% chance of selecting a point where b(x) = 0.
Then for N = 104,roughly 8 points would hit the second peak giving
IN = ff(x)dx and Neff = 30. So, with fewer points than was required
for d/b<
has been gathered to alter b(x) so that it contains the second peak.
The Metropolis scheme for evaluating expectation values, e.g.
= /h(')f(x)dx/ff(x)dx, uses If(x)l to sample the integrand, but
does not give an estimate for ff(x)d as biased selection does. This
weight function gives the minimum error in only when h(1) = 1.
The Metropolis method is easier to code than biased selection, since
there is no need to calculate probabilities, biases, or weight functions,
and since it is very simple to generate xi+ from i.. But each xi+
'1+ 1 i+l
is strongly correlated with each xi. This means that a large number
of functional evaluations must be discarded before the remaining can
be used to estimate the integrand. However, all of the functional
evaluations in biased selection are statistically independent, giving
the biased selection method an advantage when N must be small. By
increasing the complexity of the Metropolis scheme it is possible to
use a different weight function, g(x). The integral, ig(x)dx, generally
is not evaluated in the Metropolis scheme. For this reason the Metro
72
polis scheme doesn't have the ability to use multiple weight functions
simultaneously. This gives biased selection an edge in evaluating the
integrals of multipeaked functions. Biased selection can use several
weight functions to generate a set of points which are not as concen
trated as those generated by Metropolis' method. This set of points
(and associated w() )'s) can be saved and used to evaluate the integrals
of a wider variety of integrands providing useful correlations in the
integral estimates and saving the expense of generating a new set of
x 's and w(xi)'s.
The class of functions which can be used for w(x) in Eq. (7) can
be greatly extended to include, not only relative probabilities, but
also weighted sums of relative probabilities. For example, a point,
x, can be selected by two different schemes one which has biased
and unbiased probabilities of selecting an x of PBl() and PU(
the second which has biased and unbiased probabilities of selecting
an x of PB2(X) and PU2(X). This can happen in cases where one scheme
rarely selects points from a region that the second scheme frequently
does. It will be shown that if the point is selected by one method
some fraction, pI, of the time and by the second method the remainder,
P2 = 1 p1, of the time, then
P (x ) P (x)
w(x ) p Bl(xp + p2 B2(p
P 1 \ 2 T
UIl p U2(Xp
= Pl(xp) + p2w2(xp) (10)
is the correct function to use in place of w(x) in Eq. (7). (The w(p )
in Eq. (10) is the relative probability of choosing Xp only if PUl(Xp) =
PU2(xp).) It should be emphasized that the w in Eq. (6) and the w in
73
Eq. (10) represent totally different ways of selecting points. In
Eq. (6) each of the points is selected using the same weight function
g(x). But in Eq. (10) each of the xi's is chosen either using gl
(in which case this point, xli, would be used to evaluate the contri
bution, f(xli)/w(xli), to the sum) or using g2 (giving f(x2i)/w(x2i))
but not by both.
The following shows that w(x) does indeed give the correct results.
First Eq. (7), with Eq. (10) used in place of Eq. (6), is broken into
two sums, one is over the points chosen by w1 and the second is over
the points chosen w2,
plN p2N
1 1 +
IN N= I f(xli)/w(xli) + f(x2i)/w(x2i)]
i=l N =
Breaking these into a sum over the T 's and a sum over f. ( the
a ia
fia s are the fi's within a given T ) and then using continuity to
carry out the sum over the f. 's gives
a.s. 1T
IN T a N [NplWl(ja)/T + Np2w2(a)/T]f(xa)/w(xa)
a=
1 +
T mc
What's more, this estimator is still correct if some of the T 's
cannot be reached by one or the other of the weight functions, as long
as each T is reachable by at least one of the weight functions. For
a
example, if Ta was unreachable by w, (i.e., w,(1a) = 0) then
f(Xla)(/wla) = 0, since there would not be any la's. Furthermore,
74
f(2a)/w(x2a) = (NP2w2(x2a)/T)f(2a)/(O + P2W2(2a)
= (N/T)f(x2a)
and everything follows through as before. All of these methods have
the same biases that simple Monte Carlo does: they are exact as No
and their errors are proportional to 1/vN. The advantage of these
methods is that they can greatly reduce the proportionality constant
for 1/vl.
75
III. DETAILED BIAS SELECTION SCHEMES
The following will give explicit details on how to choose an x.
with biased probability PB(Xi), once a positive definite g(x) has been
decided upon. We will not give recipes for finding the best g(x) for
a given problem, but the "physics" of the situation frequently suggests
which regions are the "important" ones. First, Ns i.'s are selected
uniformly at random, and a g(x.) is calculated for each of them. Next,
a series of partial sums is calculated
SO = 0
Si = Si + g(xi)
Ns ,
SN = g(ix
s i=l
A random number, R, is then chosen uniformly from the interval R = 0
to R = SNs. Then a "p" is determined such that
Sp1 R Sp
Finally, x is picked and the remaining N 1 x.'s are discarded.
p s 1
The biased probability of choosing xp is proportional to g( p)
PB(x) = g(/SN (11)
while the uniform probability of choosing i is
P ) N (12)
76
The relative probability of choosing x is found by substituting Eq. (11)
and Eq. (12) into Eq. (6) to get
Ss s
Note that,. since SN depends on which set of Ns x's is selected,
w(x ) is a multivalued function. But also notice that w(x ) is dis
tributed about 9(xp) in exactly the same fashion, and for exactly
:g()dx 1 ) d
the same reasons, as s
the same reasons, as N SN is distributed about fg(x)dxleading
s s
to the following interesting consequences. Consider what happens to
1 na
Sf(xi) as T~ (causing T 0) and as n (and therefore N) x~.
i=1 w(xi) a a
As in the derivation of Eq. (8) continuity of f(x) and g(x) demands
that
1 naf(), f ) na 1 Ns
a i g(x.)
iw(=T N g i=1 Ns j=1 J
n a)m N g ('% a
a n
But note that
w( a) . g( a) .
Ns fg(7)dx
Therefore, as T becomes infinite
1 f(x1.) a.s. f(xa)
w1  n a
N na a
i= a w
and as Nx, IN If(x)dx for all values of Ns. For small values of Ns
the size of a is partly due to the sampling variance of SNs and partly
due to the fact that g does not imitate f perfectly. Note that for
Ns 1, simple Monte Carlo is exactly recovered, while for N ,
importance sampling is obtained. The best value of Ns depends on how
closely g mimics f and on the relative costs of evaluating g and f.
77
Frequently, it is more convenient to use two (or more) independent
weight functions, gl and g2 [4].
An example of this is the calculation of the energy
integral for an H2+ molecule. The potential energy causes
the integrand to be singular for two separate instances 
when the separation of the electron from either of the two
nuclei goes to zero. Then one weight function could be
designed to follow the singularity due to the electron's
close approach to one nucleus, and the second weight
function could follow the singularity due to the electron's
close approach to the second nucleus. Neither could be
suitably used by itself, because the remaining singularity
would make a excessively large.
Then gl could be used (in the same way g was) with probability p, to
select a fraction, plN, of the points, and g2 could be used to select
the remaining p2N points, p2 = 1P1. Since the probability of selec
ting g1 (g2) is Pl (P2), and the probability of selecting x once
g1 (g2) has been selected is gl(X )/SN (92(Xp)/S2N ), the total
(s S
biased probability of choosing x is
p
PB(xp) = P(gl)P(x pgi) + P(g2)P(x g2)
=Plgl(Xp)/SNs + P292( )/S2Ns (14)
where P(x gl) is the probability that p is chosen given that g
was. The uniform probability of choosing x was ; therefore, the
relative probability of choosing x is
w(Xp) = plg1(Xp)/(SlN ) + P2(x)/(S2 ). (15)
Note that Eq. (15) is just a particular case of Eq. (10).
78
One problem with biasedly selecting an xp by the above methods
is that a can be increased by rare occurrences [9]. For example,
suppose a large number of the Ns x.'s happened to be selected uniformly
at random in some small, but "important", region, Ta. Then if X
happened to be within Ta, PB(x ) would be anomalously small because
a large number of the g(x.)'s would be the same size as g(x ). This
would, in turn, cause w(x ) to be small; which would make f(x )/w( p)
much larger than the other contributions, f(xi)/w(xi), to Eq. (7)
from Ta where more probable sets of x.'s were selected. The effect
a 1
of the anomalous behaviour of f(xi)/w(xi) within an "important" Ta
on a can be mediated by the following averaging process. Once an
x has been chosen by g(x) and a biased probability of selecting x ,
PY( p), has been calculated, then select a new set of Ns xi's by
keeping x and reselecting the other Ns1 xi.s. Then calculate the
)t s 1
new PB(Xp), call it P"(x ), assuming that x was also chosen in the
second case. The biased probability to use in Eq. (6) will be
+ 1 > .
PB(p) = ) + PpB ()
1 g 9(xp)
2 +
g(xp)
2 g( +. + g(P) + + g(Ns
g(x,) + ... +g(x) + ...+ g(x
The uniform probability of choosing x is still 1/Ns, and therefore,
the relative probability of choosing x is
p) B(Xp
W(xp) = PB(xp)/(N ).
S
79
Great care must be taken to perform the averages precisely as
indicated in order to be sure that the biased probability is correct.
It is only when PB is calculated correctly that w, as defined, corrects
41
for the biasedly selected xi's. For example, although it might reduce
the variance of PB to average the two denominators and then invert
them the resulting number would no longer be the biased probability
of selecting x It should be noted that this must be done for every
xp, and as many PB(x )'s can be averaged as desired. Furthermore,
f(x) is only calculated one time, for each P This offers an inex
pensive way (when f(x) is expensive to evaluate while g(x) is not)
to reduce the dispersion introduced by the above methods of calculating
PB(x).
80
IV. EXAMPLES WITH NUMBERS
This section provides an example that explicitly demonstrates
the subtleties of the biased selection methods, and how the biased
probabilities must be calculated and averaged in order to obtain the
correct results. Furthermore, it demonstrates that the results are
exact, even for finite Ns.
First, suppose the value of a twodimensional integral, ff(x,y)dxdy,
is desired. Assume for simplicity that the unit square, over which
the integral is being evaluated, can be divided into nine equalarea
subsquares, Ta's, so that f(x,y) is a constant in each one of them.
Label the subsquares as follows:
y
1
2 T7 T8 T9
T4 T5 T6
3
0 T1 T2 T3
1 2
3 3
and let the integrand be
f(x,y) = 0 (x,y) j T9
= 9 (x,y) e T9 (16)
Three different Monte Carlo integrations (simple, biased selection
without averaging, and biased selection with averaging) will be per
formed with this function to see how close each comes to the answer, 1.
This allows a comparison of the three methods with one another.
8T
First, simple Monte Carlo: each of the points, (x,y)'s, will
be selected uniformly at random from the unit square. Therefore,
each of the Ta's will be equally likely to be sampled, with a proba
bility of . If a large enough number of points, N, is used to
evaluate Eq. (7) then
N [ 1 1
IN = [Nf(TI) + Nf(T2) + ... + N (T9)]= 1
and
= Nf (Tl) + ... + Nf 2(T9)] = 9,
where, when x e T1, f(x) is written f(T1). Using Eq. (3) and the two
above equations give
81
Second, bias selection without averaging: two points will be
selected from which one point will be selected with a biased proba
bility determined by some g(x,y), as in Eq. (11). This point will
then be used to evaluate one of the N contributions to Eq. (6). Let
the weight function be
g(x,y) = .01 for (x,y) e T2, T4, T7, or T8
= .09 for (x,y) e Tl, T3, T5, or T6
= .81 for (x,y) e T9. (17)
If the points are selected uniformly at random from the unit square
there will be nine equally. probable cells they can be selected from.
This gives 9(9) or 81 possible ways of choosing two points from the
nine cells. Therefore, the probability of selecting two points from
82
1
a particular pair of cells is 8. Only if the pair of points is
chosen from one of the seventeen pairs of cells in Table I can a
nonzero contribution be made to Eq. (7). The contribution to Eq. (7)
from each given value of f( )/w(x) is
N x (Probability of selecting a pair with given w) x (Bias
probability of selecting the point from T9from the given
pair of pointswith the given w) x (f(x)/w(x)).
Using this and Table I yields
1 1 9 9
I N()()( ) + N(8)(. + N(8 )(.9878)( 9756 1
<(f/w)2>N = 5.2469.
Using Eq. (3) and the two above equations gives
1
= 4.2469 2
k N
The ratio of the error estimate for simple Monte Carlo to the error
estimate for biased selection without averaging is approximately .73.
For the above integrand and weight function biased selection without
averaging gives a 27% smaller error estimate than simple Monte Carlo.
Third, biased selection with averaging: after a point has been
selected from a pair as in the second method a new pair of points is
generated by keeping the previously selected point and selecting
the second point uniformly at random from the unit square. The two
relative probabilities obtained from the two pairs of points will then
be averaged, as in Eq. (15), to obtain w(x, y). Table II,when used in
conjunction with Table I, can be used to calculate IN and a
Table I. The necessary probabilities and relative probabilities for the evaluation of Eq. (7) (using
Eqs. (6), (13), (16), and (17)) are given.
17 pairs that [g(l)+g(2)], I Probability of Bias probabil Uniform probabil w for case where
can give non (2) is for point selecting the ity of selec ity of selec selected point
zero f(x,y) from 1st (2nd) given pair with ting point from ting point from was from T9
cell of pair the given w T9, given the T9, given the
pair and w pair and w
1 1
T9 T9 1.62 1 1 1
81 1
T9 T1
T9 T3
T9 T5 8 .81 1 .9 1
T9 T6 81 .90 2 .5
T1 T9
T3 T9
T5 T9
T6 T9
T9 T2
T9 T4
T9 T7 8 .81 78 1 .9878 1.
.T9 T8 82 8 . =.9878 5 1.9756
T9 T8 81 .82 2 .5
T2 T9
T4 T9
T7 T9
T8 T9
CO
Table II. The second relative probability of selecting the same point by a different path.
Possible regions Probability of w2 where one Average of w1 and w2, where w1 has the indicated
to rechoose 2nd selecting point p t is in T9 values
point is in T9
point from from given re and the other
gion with given is in given w2 + 1 2 + 1.8 w2 + 1.9756
w2 region 2 2 2
1
T9 9 1 1 1.4 1.4878
T1
T3 4 .81
T5 4 1 .8 =1.8 1.4 1.8 1.8878
T6 (.81+.09)
T2
T4 4 .81
T7 i 81 1.9756 1.4878 1.8878 1.9756
T8 y(.81+.01)
00
P.
85
i = N(81 1.(4 1.4878) + N(8)(1)()() +
4 _9 9 8 19
N( )(.9)( )(8 + .8878) + N(81)(.9)( )(
9981 9 1.) .878 .
N( )(.9878)( )( 7) + N(8 )(.9878)(1)( 9
81 91.8878 +1.9756~ U81 9 1.4878
= 1
and
< ( T)2>N = :5.1144.
Using Eq. (3) and the two above equations gives
4.1144 1
N
Y = ( N)2
The ratio of the error estimate for biased selection with averaging
to biased selection without averaging is approximately .98. For the
above integral and weight function biased selection with averaging
gives a 2% smaller estimate than biased selection without averaging.
Note that it also approximately doubles the biased selection overhead.
How useful this averaging is depends on the relative expense of eval
uating f and w and is therefore highly dependent on the particular
problem. The more expensive f is to evaluate, relative to w, the
smaller the improvement in c needs to be for the averaging to pay for
itself.
86
V. HYDROGEN MOLECULE
As a final example, this section contains a calculation of the
ground state energy, using the eleventerm JamesandCoolidge wave
function [10], for dihydrogen at a nuclear separation of 1.4 Bohr
radii (the potential minimum). Note that James and Coolidge did not
use Monte Carlo methods. Instead, they analytically evaluated the
more than 140 integrals involved. This example incorporates a number
of biased selection methods which can also, in principle, be used for
other atomic and molecular systems.
The following is useful for calculating the required relative
probabilities. The probability density function, p(x), of a random
variable, x, gives the probability of selecting a given x, x0, be
tween x0 and x0 + dx as
P(x0) = p(xo)dx. (18)
Frequently, the density function pl(h(x)), for a function, h(x), of
a random variable, x (whose probability function is known to be
Po(x)), is desired. In this case [7]
pl(h(x)) = p0(x)l x (19)
For the special case where x is a random variable, Ri, which is
selected uniformly at random from the interval R. = 0 to
Ri = 1, the probability density function is po(Ri) = 1,
therefore, the probability density function for some
87
function, h(R), of this random variable is
Pl(h(R)) = 1IdhRT) (20)
As an illustration, consider the calculation of the relative
probability for randomly selecting spherical coordinates:
r = R1M, p = 2R21, ( = 27R3 (21)
(Ri's = different random numbers between zero and one, M = the maximum
r, u the cosine of the angle between r and the zaxis, and E
the angle between the projection of r in the xy plane and the xaxis).
Selecting spherical coordinates randomly would have an advantage over
selecting cartesian coordinates randomly when, for example, the
integrand was basically a function of r (e.g., the r2 part of PU
cancels the 1/r singularity of the potential). Since the R.'s
are chosen uniformly at random and the spherical coordinates are
functions of R (e.g., r = h(R) = R1M), the probability density func
tions are, using Eqs. (20) and (21),
pl(r) = 1/M, p2(p) = 1/2, p3(G) = 1/27. (22)
From Eqs. (18) and (22) the biased probability of selecting these
coordinates is
PB(r, dr) l (23)
M 2 27 (23)
The uniform probability of selecting these coordinates is the ratio
of the infinitesimal volume around them to the total volume, V,
PU(r, p, () = r2drdpdq/V. (24)
88
Finally, the relative probability of randomly selecting spherical
coordinates is, using Eqs. (6), (23) and (24),
wsphere(x) E Wsphere(r, , ) = V/(4wMr2). (25)
The integration region will be a large box of length Lx (Ly or Lz)
in the x (y or z) direction. The box is large enough to contain all
of the region where Y(x) is essentially nonzero. Keeping the box
finite allows the random numbers to be generated easily without
introducing singularities. The two nuclei are fixed on the zaxis
with the point halfway between them coincident with the box's center.
The coordinates of the first electron will be selected in four ways
and those of the second in five (four ways being the same as for the
first electron). One way, which guarantees that the entire integration
region is sampled, is to select the electron's coordinates uniformly
at random from the box. A second way, which accounts for the rather
large probability of the electrons being there, is to select the
electron's coordinates with respect to the midpoint between the
nuclei. Third and fourth ways, which account for the nuclearelectron
correlation, are to select the electron's coordinates with respect to
one or the other of the two nuclei. The fifth way, for the second
electron, which accounts for the singularity due to the electronelectron
potentials, is to select the second electron's coordinates with respect
to the first electron's.
An electron's coordinates can be selected uniformly at random
from the box by letting its coordinates be
x = R1Lx, y = R2Ly, z = R3Lz.
89
Using Eqs. (6) and (20) gives the expected relative probability of
selecting the coordinates uniformly at random from the box of
wB(x) = 1.
An electron's coordinates can be selected with respect to the
box's center by letting
rc = .iMc, p =2R21, R = R327 (26)
(rc = the electron's distance from the center). Using Eqs. (20) and
(26) and letting h (Rc) = vicM gives the following probability density
function
Pl(r) = 2/ /Mc. (27)
With Eqs. (18),(22) and (27), the biased probability of selecting
rc, ~, and ( can be calculated; then, using Eqs. (6) and (24), the
relative probability of selecting the electron's coordinates with
respect to the center is calculated to be
wc(r, i, p) = V/(2TMrc). (28)
An electron's coordinates can be selected with respect to either
of the two nuclei by using a weight function, g(y). An rNl is
selected by requiring
m rN1
A(R1) R 1g(y)dy = / g(y)dy. (29)
0 0
The remaining coordinates, t and , are selected as in Eq. (26). The
probability density function for rNl can be calculated using Eq. (10),
if it is first realized that rNl is a function of A which is, in turn,
90
a function of R. Using Eqs. (19),(20) and (29) leads to
Pl(rN1) = p(A) IdrA(A)
=1dR dA
dA(R1) drNl(A)
S g(rN ) (30)
fg(y)dy
0
Therefore, using Eqs. (18), (26), and (30),
PB(rNI ) = g(rN)drNld du (31)
47T/g(y)dy
0
and the relative probability of selecting the electron's coordinates
with respect to the first and second nucleus, respectively, can be
calculated (using Eqs. (6), (24), and (31)) to be
WN1(rN1' ll)' g(rN1)V
4Trr~1g(y)dy
N10
and
wN2(rN2 2' 2) = g(rN2)V (32)
4rN2f/g(y)dy
In practice, a good weight function is g(x) = x22(x), where y2(x)
is, in general, the HartreeFock nuclearelectronic correlation
function (in the example'(x) is the ground state hydrogen wave function).
The HartreeFock functions can be found in tables such as reference [11].
Note that the x2 cancels the rl1 in Eq. (32) and that Y follows the
system's wave function quite well for large nuclear separations. To
greatly lessen the amount of computer time required to biasedly select
an rNl, x 22(x) is evaluated at some number of points, say K, and g(x)
is formed by connecting these K values by line segments. To determine
xi 1 i
an rN1, R1K (Qi =X g(x)dx = I (gj + gjl)(xj xj. ), and each xj
N1 K o j Jg d J 2
is one of the K used to evaluate xi2(x.)) is calculated, and then a
p is located such that
Q < RQK Q+
then
rN1 = gp/A + [g + 2(A)Q6]T/A
(A E 9p+l Sp and Q6 p+l Qp).
Xp+1 
If rNl is not one of the K x 's then g(rNl) in Eq. (32) is
g(r1) gp + rNl x (gp gp)
Xp+1 p
Since what happens as an electron approaches a nucleus closely is of
interest, g(r) in Eq. (32) is altered slightly at the originit is
made greater than zero, usually we take g(0) = g(l )to make it more
probable to place an electron close to the nucleus.
The second electron's coordinates can be selected with respect
to the first's as in Eq. (26), with Mc replaced by ME Therefore,
the relative probability of selecting the second electron's coordinates
with respect to the first's is
wE(rE, P, ,) = V/(27TMrE).
Note that rE, in wE(x), is exactly what is needed to account for the
singularity of the electronelectron potential.
92
If each of these methods is used with equal probability then
W(Xl, x2) = WB(Xl) + wc(X) + WNl(X) + wN2(X)] x
1 ) + w *
[WBg(X2) + wc(2() + ( WN2(X2) + WE(2)]. (33)
If desired, one more averaging can be done
w'(x~, x2) = w(x, x2) + w(x2, X)].
The numbers in Table III were calculated using the relative
probability in Eq. (33). The 10,000 point estimate required 14 seconds
of Amdahl 470 time, and gave a .002 Rydberg accuracy. The calculation
of serves as a check that VT was calculated properly, while the
calculation of serves as a check that the relative probabilities
were evaluated correctly. Note that being correct does not
guarantee that the relative probabilities are correct. This is because
HY. = EY. for all i, and therefore,
= E /w/f /wi = E
regardless of the w.'s used. It is interesting to note, by comparison,
that an unbiased Monte Carlo estimate for 10,000 points has an
Neff < 2, = 2.6941, = 7.3062, and = 2.4828. In addition,
the error estimates are unreliable instead of decreasing as 1/VT
they remained constant from N = 5,000 to N = 10,000. This is because
smaller and smaller interparticle separations are sampled as N is
made larger and larger. Even though these separations are sampled
the least frequently they make the largest contributions. These contri
butions are increasing enough, as N increases, so that they keep a from
dropping.
93
Table III. Results for the H2 molecule using the 11term Jamesand
Coolidge wave function for a nuclear separation of 1.4 Bohr radii.
Number of
Monte Carlo
points a a a
N I.
100 eff 21.2 26.0 21.2
T 2.277 .086 2.344 ; .083 2.369 .099
V 4.40 .28 4.16 .44 4.75 .39
H 2.354 .019 2.369 015 2.350 .033
I I ,
I I I
200 eff 43.6 45.0 41.5
T 2.333 .069 2.339 .062 2.273 .070
V 4.69 .54 3.94 .28 4.90 i .28
H 2.361 .017 2.3520' .0097 2.336 .020
400 eff 81.2 94.7 87.2
T 2.306 .049 2.340 .038 2.271 .046
V 5.02 .44 4.16 .21 4.85 .21
H 2.361 .013 2.33521 .0069 2.355 .015
10,000 eff 2222.0
T 2.3475 .0089
V 4.682 .055
H 2.3446 .0023
Three independent runs are tabulated. Neff (E iH'i/wi)2/(iHi/wi)2,
E , H the system's Hamiltonian. Note the anomalous accuracy
of H and that the estimated error, a, does not display 1//IT conver
gence until Neff = 43. Also, note that =  and = 2, in
accord with the virial theorem. All expectation values are in Rydbergs.
94
CONCLUSION
Systematic methods, such as trapezoidal, of evaluating multi
dimensional integrals have the disadvantage that they require an enor
mous number of points before they can adequately sample the integrand.
Biased selection allows one to reduce N in exchange for using his
knowledge of the integrand. For example, HartreeFock wave functions
can be used as weight functions to find the integrals of more exact
wave functions. In addition, biased selection allows one to adequately
sample the important regions where the electrons approach the nucleus,
or each other, closely. (These regions are important because small
1
errors in the wave function will be amplified by the dependence of
the potential. But since these regions occupy a relatively small
volume, they will not be likely to be sampled. This will in turn
cause the error estimate to be deceptively small.) The price paid
for this flexibility is that the points must be distributed in a
nonuniform manner and the correct weight function to compensate for
this must be calculated. A final important aspect of biased selection
is that it allows for the possibility of finding unexpectedly large
values of f/w. (One technique for doing this is to uniformly and
randomly sample the regions believed to unimportant, for dimension
ality < 12. For more dimensions other methods must be devised:
e.g. T2", where a<1, could be used.) These will undoubtedly ruin
a good many papers which seemingly needed just slightly smaller
deviations, but the discovery also enables one to find better weight
functions, and eventually, to find the correct value of the integral.

Full Text 
PAGE 1
METALLIC HYDROGENÂ— A TWOBODY APPROACH TO BAND THEORY AND A QUANTUM MECHANICAL TEST OF DIOPHANTINE METHODS BY STEPHEN A. BURDICK A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1980
PAGE 2
Copyright 1980 by Stephen A. Burdick
PAGE 3
ACKNOWLEDGEMENTS The author wishes to thank the Kirkland Data Center for the Division of Highway Safety and Motor Vehicles for the enormous amounts of computer time this work required. The Northeast Regional Data Center was extremely helpful in making the connections and transmissions of data between Gainesville and the D.H.S.M.V. at Tallahassee. The author is also indebted to Dr. Coldwell for many lively discussions, his development of the biased selection method and coding of the minimization routines, and the hours of time he had to put in to get the required amounts of computer time. The author also had numerous enlightening discussions with Dr. Broyles and Dr. Burdick. The many hours of typing were put in byhiswife, Linda. Finally, the author would like to express his sincerest appreciation to his entire family for their continual support and understanding throughout this most challenging period. n 1
PAGE 4
TABLE OF CONTENTS Page ACKNOWLEDGEMENTS iii ABSTRACT v INTRODUCTION 1 PART I: METALLIC HYDROGENA TWOBODY APPROACH TO BAND THEORY 2 History 4 The Model Hamiltonian 7 The Wave Function 11 The Monte Carlo Error and Finding the Eigenfunctions of Equation (4) 15 A TwoElectron Approach to Band Theory 17 Results 32 Conclusion 33 List of References 34 PART II: A QUANTUM MECHANICAL TEST OF DIOPHANTINE METHODS 36 Hasel grove's Method 38 Hasel grove's Method Used by Itself 42 Hasel grove's Method Used in Combination with Biased Selection 47 Further Tests of the Blend of Haselgrove's and the Biased Selection Methods 51 Monte Carlo Estimation of the Error 55 List of References 58 CONCLUSIONS OF DISSERTATION 59 APPENDIX A: THE BIASED SELECTION METHOD 61 APPENDIX B: SOME c^'s WHICH FORM SOLUTIONS TO EQUATION (4) OF PART I 96 ALPHABETIZED BIBLIOGRAPHY 115 BIOGRAPHICAL SKETCH 118 IV
PAGE 5
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 METALLIC HYDROGENÂ—A TWOBODY APPROACH TO BAND THEORY AND A QUANTUM MECHANICAL TEST OF DIOPHANTINE METHODS By Stephen A. Burdick December 1980 Chairman: Arthur A. Broyles Major Department: Physics The first part of this work demonstrates the feasibility of calculating the properties of bulk hydrogen from a periodically repeated cell of two "atoms" using biased selection. It further shows how conventional band theory can be extended to use 2electron wave functions. BornOppenheimer potentials were generated for a few t's for the first three bands and for r^ = 1.23ag, 1.48an, 1.97an, 2.96an and 10.84ag (J r^ e mean volume per electron). An insulatorconductor transition was found to occur between r = 1.7an and r = 1.9an, s B s B' while the solid was still molecular. A molecularatomic transition was predicted, by extrapolation, to occur for r = 1.1. Finally, a significant broadening of the kaveraged BornOppenheimer potential takes place at r = 1.97an, where the first signs of an insulator
PAGE 6
cofiductor transition occur. (The energy is between 1 .17Ry/atoin and 1.15Ry/atom when the internuclear separation is between 1.25an and 2,20ag.) This effect is bound to have thermodynamical effects beyond those normally associated with insulatorconductor transitions. The second part of this work details the application of Diophantine methods to the evaluation of multi dimensional quantummechanical expectation values. It shows that a tjconvergence was found for a substantial range of N, the number of multidimensional integration points. For the same accuracy as biased selection, Diophantine methods required onefourth the number of points, for the remaining range of N. VT
PAGE 7
INTRODUCTION This work is most logically divided into two parts. Part I, Metallic HydrogenA TwoBody Approach to Band Theory, details the methods and results of a solidstate type of calculation for bulk hydrogen. Part II, A Quantum Mechanical Test of Diophantine Methods, details an exploration into the possibility of getting a convergence for general multidimensional integrals that is faster than the Monte Carlo convergence of 1//R". Parts I and II are selfcontained: the equations and references are numbered as if the other part were not present, and the understanding of one is not required for the understanding of the other.
PAGE 8
PART I: METALLIC HYDROGENÂ— A TWOBODY APPROACH TO BAND THEORY The fact that hydrogen becomes metallic under high pressure is interesting for several reasons. Firstly, theoretical calculations exist that have found the superconducting transition to be close to room temperature [1], [2]. Secondly, the interior of Jupiter contains hydrogen under high pressure (P > 30Mb) [3]. Thirdly, the conceptual distinction between an insulator and a conductor is clear cut. If, experimentally, hydrogen is placed under pressure with a small voltage applied to attached electrodes, either a current will flow and it is a conductor, or no current will flow and it is an insulator. If, theoretically, the first two energy bands are calculated, either they overlap and hydrogen is a conductor, or they are separated and it is an insulator at 0K. In this work the usual method for calculating the properties of a crystal was extended, from treating the electrons one at a time to treating them two at a time. The explicit inclusion of electronelectron distances in the wave function allowed the bulk of the correlation energy to be included. An approximate wave function for the system was constructed from 2electron wave functions in analogy with conventional band theory. The overlapping of the first two energy bands indicated that at 0K the molecular solid becomes a conductor for an r between l7an and 1.9an (r for liquid hydrogen is 3.4an, 5r^ ~ N"' '^ ^'^ ^^^ number of electrons contained in the crystal volume.
PAGE 9
L^). The kaveraged BornOppenheimer potential indicated that at 0K the solid was still molecular for the smallest r used, 1.23aÂ„. However, if the data are extrapolated to smaller r then the molecularatomic transition appears to occur at r = l.lag. The energies found here can be used to form the partition function, from which the properties of the crystal can be calculated for temperatures less than the ionization temperature.
PAGE 10
History J. D. Bernal was the first to suggest that insulators would become conductors under high pressures [4]. Prior to 1974, most calculations were made assuming that the molecularatomic and insulatorconductor transitions coincided for solid hydrogen. Atomic hydrogen has one electron per fundamental cell. Since each energy band can accomodate two electrons per fundamental cell, the valence band is only half full. This means that an atomichydrogen solid is necessarily a metal. Therefore, it is not clear whether assuming the two transitions coincided was a matter of belief or merely computational convenience. However, since 1974, theoretical evidence has been mounting, which this work further supports, that the solid becomes a conductor while it is still molecular [5], [6]. This makes the insulatorconductor transition a continuous function of r second order, because there is no abrupt change in symmetry, from molecular to atomic, to create a discontinuity. Wigner and Huntington were the first to calculate the properties of an atomichydrogen solid [4]. They used the cellular, or WignerSeitz, method and found an energy minimum of 1.05 Ry/atom at an r of 1.63ag. Later calculations essentially reproduced these results. The solid hydrogen calculations fall into three general categories. There are calculations to determine the properties of the molecular solid, generally by summing over predetermined H^H^ potentials [7], sometimes by using a low density expansion such as Brueckner's method [8].
PAGE 11
These methods commonly assume the potentials to be density independent when they extrapolate the equation of state to the molecularatomic transition density. This work casts doubts on this assumption near the molecularatomic transition. The HH potentials were found to depend on r which makes it highly likely that the H^H^ potentials will also. Furthermore, it is not clear how an H2H2 potential should be defined when the molecules begin to dissociate. Some of the molecular solid calculations accurately reproduce the experimentally available equation of state (0 Â£ P ^ .025Mb) [7]. Other calculations are on the properties of an atomic solid; 0function [9]; Wannier function [10]; Green's function [10], augmented plane wave [11]; exact HartreeFock [12], [13], [14]; perturbation theory [15], [16]; cellular [4], [17]; alternate molecular orbitals [18]; and quantum statistics [19]. Ross and McMahan [11] and McMahan [5] state that the uncertainty in the correlation energy is the outstanding problem for these methods (the energies agree to within .04 Ry/atom while the uncertainty in the correlation energy is about .13 Ry/atom). They further state that this uncertainty makes the molecularatomic transition pressure uncertain by a factor of about two or three. This work includes the properties of a bodycentered cubic atomichydrogen solid as a subset. Finally, there have been a few calculations since 1974 that, like this work use the same formalism to calculate the properties of both the atomi c and the molecu lar forms of solid hydrogen [5], [5]. The potential and periodicity are treated somewhat more crudely here than in other methods, but this work uses much bette r wave functions that contain e lectronelectron distances explicitly.
PAGE 12
The molecularatomic transition pressure is usually calculated from the slope of the common tangent to the atomicsolid energyversusdensity and the molecularsolid energyversusdensity curves. This pressure is also frequently assumed to be the insulatorconductor transition pressure. This work distinguishes the molecularatomic transition, which at 0K did not occur for any density considered, from the insulatorconductor transition. The molecularatomic transition takes place (at 0K) when there is no longer a minimum in the BornOppenheimer potential for any nuclearnuclear distance other than the one corresponding to a bodycentered cubic atomic lattice. Also, the energy at the minimum of the kaveraged BornOppenheimer potential and the energy at an internuclear separation corresponding the the bodycentered cubic lattice will agree at the r where the molecularatomic transition occurs. The data of this work were extrapolated to find this transition. The insulatorconductor transition takes place when the bands overlap at the equilibrium separation of the nuclei. As of 1977, Ruoff felt that there was not any convincing experimental evidence for an insulatorconductor transition in hydrogen [20]. However, a more recent experiment claims to have found one [21]. Also, there are ongoing experiments to produce an insulatorconductor transition in hydrogen [22].
PAGE 13
The Model HamHtonlan It is currently impossible to solve the N (10^^) particle Schrbdinger equation for a general system. A standard method of dealing with this problem is to model the real system with a periodically repeated unit cell of fewer particles. In this work, two particles are placed in the unit cell and the Nparticle Hamiltonian, H^, becomes H^(r^, ... r^) = iH^ir + I), (1) I where the coordinates of the i electron are given by r. r = ^'^Ix' "ly Iz' ^2x' ^2y' ^2z) ^"^ ^ ^"ix' "ly' "iz' "2x' "2y' "2z)^The n's are integers and "a" is the length of a side of the cubical unit cell. Since Hj^ is then a sum of 2particle Hamiltonians for the different cells, the Nparticle wave function will be an antisymmetrized product of 2electron wave functions for the cells. In the BornOppenheimer approximation, Hr,(r) is written as H2(r) Â— V^ p^V^ p^+^+ p^/p^+ V(r), ^ ^Al ^ ^B2 '^AB ^12 ^Bl ^A2 where A and B refer to the nuclei, 1 and 2 refer to the electrons, and V'(r) represents that part of the potential that is not explicitly written down. The P's are the shortest distances between two particles, which can be in adjacent cells, of a given kind. Figure 1 illustrates this for a periodically repeated 2dimensional unit cell.
PAGE 14
1 B A 2 1 B A 2 1 B A 2 1 B A 2 [xi i B A 2 1 B A 2 1 B A 2 1 B A 2 B A 2 (a) (b) ^1 B A 2 B A 2 A B 2 B A" B A 2 I B A 2 B 'A 2 1 B A 2 Figure 1. A periodically repeated 2dimensional unit cell is pictured. The interparticle interactions that are considered in this work are represented by lines that connect the particles, A configuration where all interparticle distances, P's, are within the unit cell is pictured in (a). In (b) some P's are between particles in the unit cell while others are between a particle in the unit cell and a particle in an adjacent cell. In both cases, the P's are between the closest particles of a given type.
PAGE 15
If lines were drawn to connect each particle in the unit cell with all of those particles to which it is not already connected, then these connections would represent the same interparticle interactions that V'(r) does. A general group of terms would look like Â— + '^lA' ^ir + where the primes are for particles not explicitly rjg, r^2' included in H^(r). The 1 in the numerator is because the remaining half of the interaction will be included in H2(r) for the cell which the primed particle occupies. Note that this group of terms is approximately zero. In this work, V'(r) is neglected. It could have been included by using an Ewaldtype sum, but the desirability of including V'(r) is questionable. For one thing, the imposed periodicity certainly misrepresents nature. Thus, including V'(r) might actually lessen the correspondence between the model and the real system. Secondly, the real system V'(r) is bound to be small due to charge neutrality and screening, especially for a conductor. The periodicity of Wn{r) i.e., H2(r) = H2(r + I) > was used to advantage in constructing solutions to H2(r)*(r) = E$(r). (2) If $(r) is written as Hr) = e'^'^'^Uj^Cr) (3) (this can always be done; e.g., let Uj(r) = e~ *'^ (r), then e Uj^(r) = e e' *(r) = (r)), then the translational invariance of H2(r) demands that ur(r) = ur(r + i) That is to say that u^(r) will be periodic without approximation. However, ur(r) might only be approxi
PAGE 16
10 mately periodic for the real system. The benefits of using (3) will be discussed later. Solving (2) for a 2electron wave function instills a certain amount of selfconsistency in the wave function, without iterating. HartreeFock type methods try to solve a 1 electron equation where l\ p^ Â— ^ p^+ V'(r) in (2) is replaced by the interaction ^ ^B2 ^hl ^12 of electron 1 with an electron charge density, the coulomb and exchange potentials. This makes it necessary to solve these equations iteratively to achieve selfconsistency, since the electron charge density depends on the 1 electron wave functions that are being solved for. If (3) is substituted into (2) and V'(r) is approximated by zero, the following equation is obtained 1 F^l ^ Pg2 \ \ L c r^g K^2 "4" 4'^^^ '^^ t^ Â• t^ E)]u~,(r) = (4) This equation was solved in the BornOppenheimer approximation for an "Hpmolecule" in a cubical unit cell of side "a". The boundary conditions were satisfied by restricting the class of functions used to represent ur(r) (further discussion will be given later). The "molecularaxis" was along a body diagonal with its center at the cell's center. This positioning of the "molecule" allows the crystal lattice a 1. to be varied from a simple cubic molecular lattice, rÂ„n < o" 3^, to a a bodycentered cubic atomic lattice, r,g = ^ 3^. Two is the smallest number of "atoms" in the unit cell which will allow correlational effects to be calculated.
PAGE 17
n The Wave Function The imposed periodicity allows an exact representation of the solutions of (2) to be written as in (3). The kdependence of $(r) makes it possible to construct a Nelectron wave function from a product of $'s with different k's, as will be further discussed later. A 2electron wave function allows the major portion of the correlation energy to be calculated. The relatively simple and yet accurate wave function of James and Coolidge, which explicitly includes r^2 terms, was used as a basis for constructing ur(r) [23]. This choice meant i that 6dimensional integrals had to be evaluated to find the expec[ tation value of H2(r). These integrals were evaluated using the Monte Carlo method of Appendix A. The accuracy of the wave function allowed the energy to be determined to within chemical accuracy, .01 Ry/molecule, by using less than 500,000 functional evaluations. This was accomplished by making H$/c[) relatively constant, which makes the Monte Carlo error very low. The periodic boundary conditions can be incorporated in ur(r) as follows: Uj^(r) = Igir + I), (5) I where g(r + 1) is a 5dimensional generalization of g(r) = g(r) for r < j = constant for i^ 1 f
PAGE 18
12 This form of g(r) is not mandatory, but it greatly restricts the sum in (5) to only include the unit cell. Explicit forms of g(r) will be detailed shortly. Notice that in the case considered by Wigner and Seitzan ion, vn'th a spherically symmetric potential, and a single electron per unit cellthe above reduces to the WignerSeitz boundary conditions: u. (r) = constant for r ^r au,, S 9r = 0. Another important similarity exists between the WignerSeitz method and this work. There is an implicit electronelectron correlation that is assumed for the Nelectron wave function. In the WignerSeitz method, the assumed function is zero if any cell contains more than one electron and is otherwise unity. In this work the assumed wave function is zero if any unit cell contains more than two electrons and is otherwise unity. (Notice that this function imposes charge neutrality on the system on a cell bycell basis.) This function could be made more flexible by allowing it to vary between zero and one as the number of particles per unit cell varies. Lastly, the WignerSeitz method does not consider interactions between particles in different cells. But for one electron per unit cell, neither does this work: the distance between particles in adjacent cells is greater than I"
PAGE 19
13 A twentyfour parameter function, modeled after James and Coolidge's function, is used to represent ur(r): uj;(r) = B(x^)BCy^)B(z^),B(x2)B(y2)B(z2)S^^S^22g^Sg2Pjc where P ,p has two forms 2 2. Pg = c + C2(n^ + n2) + ^^^^^ + ^2^ ^ ^^4^1^2 "^ ^"^5^ "*" Cg(erii + ^211.2) + CJ^^^2^{e^ + Â£2) + Cg(Â£:^n2 + ^2^^ "*" Cg(Â£^ + Â£2) + 2cigU^ + 2c^^Ti,n2U, and P^ = (c^ + c^Â£^e2 + ^s^ ^ Cgu2)(Ti^ ri2) + C2(e^^ri^ Â£2^2^ "^ C3(Â£^^n2 Â£2'^i) + (Cg + c^Qu)(Â£^n e2'^2^ ^ (c^ + c^^u)(Â£^n2 Â£2^1^ "^ '^g^^i ^^2^' with u = 2f(P^2) ^1 = f(PAi^ f^Pei^ i f(x) = C^^^Ca' x) 2(a' .2 x) + Ca" .4 x)] + b}/r^g, and b = l^a 2(a'.2)' + (a'.4)^], a' = . Furthermore, 20 [ ^Al = ?/m^^^Al)+^^24' ^ m=12 with K^ = (n!l^)MIN(, 5) I
PAGE 20
14 (x) = X for X >^ = for X > (the same c 's and K 's were used for all of the S's). Finally, i^x 1Â— X B(x) = C21 + c^oS + Â€236 The form of f(x) insures that the Pip's satisfy the conditions of (5)the knots are less than pwhich drives the function to b and its derivatives to zero for x _> p. The same condition is fulfilled for the S's, in the same way, except that now the constant is 024. Finally, the B's are constructed so that they are explicitly periodic. The reason for two Pip's is discussed in the section on band theory. As previously mentioned, the 2body form of ur(r) allows the effects of 2body correlations to be calculated. The twentyfour variational parameters, c 's, are tabulated for the calculated wave functions in Appendix B.
PAGE 21
15 The Monte Carlo Error and Finding the Eigenfunctions of Equation (4) The biased selection method (see Appendix A) was used to estimate the expectation value of the Hamiltonian, IE aVIa' (6) ^ m m '^ m ^ 2 2 where E = [Hr,(r )'5(r_)]/$(r ), A = ur(r_)/w(r ), and w(r ) is a m '2^ m' ^ m'" ^ m m k m" ^ m' ^ m function that exactly compensates for the biased manner used to select r^. Note that E will vary unless ur(r) is an exact solution of m m ^ k^ (4). The estimate of the Monte Carlo error made in estimating is a = [^E a')2 2yE A' + ^yA']vyA' ^ m m ^ m m ^ m ^ m where the above sums extend from m = 1 to m = M, M is the number of 6dimensional Monte Carlo points used. Note that if ur:(r) were an eigenfunction that a becomes [E^TA^ 2E^yA'' + E^jA^lvyA^ which is zero. In other words, (6) would have zero Monte Carlo error. This fact was employed to find the eigenfunctions of (4) by varying the c 's to minimize a. After minimizing, uf;(r) is almost, but not quite, m K an eigenfunction. Therefore, a will be proportional to M~ ^ with a small constant of proportionality. Equation (4) was solved by varying the twentyfour c 's in ur;(r) 111 K to minimize a. During the minimization, the in a had to be replaced by a constant that was less than the eigenenergy of the state being sought. Two techniques were used to vary the c 's: (1) a routine
PAGE 22
16 adapted from Nelder and Mead's simplex algorithm [24], (2) a routine that used analytically determined first derivatives, approximate second derivatives, and a matrix inversion with smoothing [25], The first method starts with a test grid that has one more point than the number of parametersin this case the grid has twentyfive points in a 24dimensional spacefor each of which a is evaluated. If vectors, from the origin to the points, are used to locate the grid points, then twentyfour of the vectors must be linearly independent. The routine then generally evaluates a for a point lying along the vector extending from the point where the largest a was found through the "center of gravity" of the remaining points. If the a for this trial point is less than the second highest a (of the current twentyfive a's), then it is kept and the highest is rejected. The second routine heads in the direction for which a decreases the most rapidly. Smoothing is used to make the distance moved towards the predicted minimum only a percentage of the predicted distance. This is done both to prevent the divergences frequently encountered in gradient methods and to make the expansion of a in terms of derivatives reasonably accurate. Additionally, smoothing allows the unimportant coefficients to become zero without making the matrix singular.
PAGE 23
17 A TwoElectron Approach To Band Theory By writing the wave function as e ur(r), full advantage can be taken of the similarities between this method and conventional band theory. The importance of k here, as for the k in conventional band theory, lies in the fact that it is a good quantum number of the system and can therefore be used to label the different 2electron states. For 1electron wave functions and conventional band theory, the system's wave function is found by anti symmetrizing H' = ^ (4) ... *ir (Nl)()^ (N). If it is assumed that e(k,) < e{k^) < < e{tj then the ground state of the system will be given by anti symmetrizing "^N ^ *k ^^^*k ^^^*k ^2^*k ^^^ Â•Â•Â• *t (^'^)*lt (^)' Similar results hold 2 2 I for two electron wave functions: \ = $^ ^^ (1 2)^ Â•}> (3, 4) ^ t (N1, N); and assuming that e^ ^ ^ ^t t ^ ^t t ^ ^t t etc., and hrn n I ^12 2*^2 ^2^^3 [ that a 4>]t" I? can only accommodate two electrons, makes the ground state \ m n i the antisymmetrization of ^j^g *t it ^'' ^^^ l< ^^' ^^ Â• Â• ^ f ^^""' ^^[ 2 2 Finally, the properties that t e{t ), and <^^ (a) and (^^ (b) have in [ Â•^ "^ "^m ^m conventional band theory will also be assumed to be the properties of y i k e(k ), (br (a, b) (k has the components of k^ as its first three m ^ m k m '^ m m components and as its next, and final, three components also). Thus, the eigenvalues of (4), e(k ), can be used to find the energy bands. These bands are continuous functions of k within each Brillouin zone.
PAGE 24
18 with possible discontinuities at the zone boundaries. Remember that even though k is 6dimensional here it consists of two 3dimensional t's that are exactly analogous to the k's in conventional band theory. Additionally, a k labels unique electron states only as long as it is contained within the first Brillouin zone. If k_ extends beyond m the first zone, then it can always be reduced back to the first zone by adding a reciprocal lattice vector to it. Therefore, k's outside the first zone merely relabel a state already labeled by a k^ within the first zone. For k = 0, the center of the first Brillouin zone, (4) has a spectrum of eigenenergies, Â£ 's. These energies can be used to generate the various energy bands, e (k), as k is varied over the first zone. It is interesting to note that if k were varied from the center of the first zone to the center of the second Brillouin zonethat ur(r), if it were sufficiently flexible, could vary continuously from the ground state solution of (4) to either the first excited state or back to the ground state solution. This was the reason two Pjp's were used. Once the ground state ur(r) was found (using Pg), then a method was needed that would allow the first excited state to be found while guaranteeing that a poor approximation to the ground state was not found. This was accomplished by using a function, Pj, that was orthogonal to the ground state. Two states were required because the first two bands were needed to find the insulatorconductor transition. Since Pp and Pp are the two eigenstates of lowest energy for the hydrogen molecule, they were used to generate the first two energy bands. The bands generated by P^ and P^., e^ and e^, turned out to be two approximately parallel parabolas of the form
PAGE 25
19 Bj.k) = sJO) +d^k.k. (7) (This form is supported by conventional band theory calculations to the desired accuracy [5], [9], [10], [15].) Apparently a band had been missed that should have roughly connected Â£r(0) with ^^(kg). k^ = Â— (1, 0, 0, 1, 0, 0). This is what nearly free electron band s a theory would lead one to predict, and there are also conventional band theory calculations that have this state. If the state existed, it could be found by directly calculating the band gap, 6 at k from Pg and seeing if it was less than Â£r(ke) ~ ^g^'^s^* ^^^^ it was calculated, 5 turned out to be less than er(0) Â£p(k ), and thus the missing state had been found. The band gap at k was calculated by forcing B(x, )B(y, )B(z, ) iLr ik r B(x2)B(y2)B(z2), b(r), to be b^ = e e This creates the equivalent of the familiar forms for the wave function at the zone boundary in conventional band theory, cos(7x)c()(x) and sin(^x)(j)(x) a a The Uj(r) with the b would then be orthogonal to the ur; (r) with s ^ s the b_. This assured that the b form would give Sp^l^c) ="^^ '^'^^ ^ poor approximation to Â£,(k ). The existence of this new energy band was further verified by varying k slightly away from k and finding energies lying on the new curve. The new energy band was taken to be an inverted parabola, (i.e. d^ < in (7)) connecting Â£r(0) and the newly found e^{k ) = ep(k ) + 5. The first and third bands were found in a similar fashion: d^ was chosen to connect eAO) with Â£r(k ), L = 7 (1, 1, 1, 1, 1, 1) and do = ^ii^JLlii^ and d. was chosen kc Â• kc to connect Â£g(0) with Â£Q(k^) d^ = ^G^V_J_^G^ A few additional kr kr
PAGE 26
20 e(k)'s were calculated and found to fall on the predicted curves. These are the energy bands required to calculate the properties of solid hydrogen, almost. Before continuing, a brief discussion on why ^^oCk) was found instead of ep(k) is in order. This problem was caused by symmetry discrepancies between Pp, and P^ and h{r) The ground state wave function, Pp, has an sorbital type of symmetry about the box's center while P^ has a porbital type. Thus, for Pp to be able to give EpCk ), and not Â£3('< ), b(r) had to assume a form at k that would make P^b (r)^have an S;3orbital type of symmetry about the box's ik r "''"'^r**^ center, b(r) = e e Furthermore, Ppb would have to have a porbital type of symmetry about the box's center. Unfortunately, b(r) was not flexible enough to have these forms. Figure 2 illustrates these symmetry conflicts. The eigenenergies do not necessarily give the solid's energy bands. Close inspection reveals difficulties in constructing the manybody wave function from some 2body wave functions, if HartreeFock theory is correct. The difficulties can be most easily visualized by trying to construct a 4electron wave function from 2electron wave functions. The HartreeFock representations of Pp and Pr are equivalent to, respectively, $q = ^qO)^^{2) and $, = *q(1)cJ)^(2) *o(2)*i(l) (where ^^ and cj), are different 1 electron states). Notice that if a 4body state is constructed from and *, that three particles will occupy the same 1 electron state and antisymmetrization will make this 4body wave function identically zero. Given that $q is occupied, the next available state that a 4body wave function can be constructed from is ^^ = (j^, (1 )(]), (2) A look at the HartreeFock energies for ^sEsaijicK
PAGE 27
21 Sx C /'S C ^i i\ X^s B J $/ c ^X *^ ^ s'x c .^s c ; 5 Xs B ^ 5'' c ''^s c ; I A ''>'^ i ? s/c\; y 5x C ^ S C 5! ^ ^ D ; ^ s D s ^ X / C ^ C 11 5 'c^ s /C^. :; s'^x^^ s B ^' ^^"^ s B '' \^' s B ^> ;j,''''c ^^s c j,/''c^xS c ;5,'' c^'sS c ;; s c s^^c /;; c s n c^^s c Sx c ^^s s A s K ;; A s X s A s X s / X X <; r ,=y n '^;; p <::'' r <; r q^ r x<; $f c X c /^r C 3^ Q ^i C ^^ C /!; ;; A s V's A s^x''' A s V ;5 ;; c s^^ c ns c s^^ c ^^; c s,^ c xj; Figure 2. The nodal lines are drawn for the first excited state wave function, dashed, and for cos(x), c, and sin(x), s, in a a a periodically repeated 2dimensional unit cell. Note that neither the c's nor the's's fall on the dashed lines. This symmetry conflict, due to the inflexibility of b(r), was the reason that the eigenenergies homed in on e^i^) instead of S2(i<)'''lote that the extra nodes introduced by either c or s increase the kinetic energy of the electron.
PAGE 28
22 these states will show how ep(O) and sJk) must be corrected: $Â„ has energy Sqq = 2eg, $, has energy Sg, = e^ + e, and ^ ^^^ energy Â£,, = 2Â£, Thus, to find the energies of the next occupiable states from the unoccupiable states given by Pr, the correct Eo^'^) ""^ given by 2ec(0) EqCO) and the correct SoCk) is given by 2er(i<) Â£p(k). The bottom of the second band, e^^k), did not need correcting because both electrons had been placed in an excited state. The resulting band energies are plotted in Figures 38. As a possibly very important aside, it should be pointed out that the antisymmetrization of ^ni'^ 2)$, (3, 4) does not give zero for general 2electron wave functions. Therefore, if HartreeFock is correct, the nonzero part left after antisymmetrization must be a bad approximation to the wave function. Thus, the energy of the antisymmetrized '5q(1 2)$, (3, 4) will be much larger than e^q + e,,. However, conventional band theory is known to be incorrect for certain systems [26], [27]. This breakdown could be due to the neglect of such states as "^0(1, 2)$, (3, 4) which might actually have energies such as eQQ+ Eq^ At any rate, it is rather interesting that the antisymmetrizing of a product of 2electron states to form manyelectron states does not generate a "Pauli principle" as it does for 1electron states. Now that the energy bands have been found the kaveraged potentials can be calculated. This is accomplished by placing y^Â— ysdk electrons in a volume, dk, in kspace. The parabolic form of the bands means that a ksphere should be filled until it makes contact with the sides of the first Brillouin zone. Then only that part of the ksphere inside the first zone will continue to be filled until the energy e(k) is equal
PAGE 29
23 o M n3 01 Dl S_ CD o 1 1 2 Â•"ab ^^B^ Figure 3. The top, Â£^(k^), and bottom, e^(0), of the first band, the bottom, Â£2(1^5^' ^ ^^^ second band, and the BornOppenheimer potential, ^, versus internuclear separation are plotted, r = 1.23. The vertical lines represent the energies calculated by Monte Carlo and their a's. Since Â£2^'^$^ "^ '^l^'^c^ ^^ ^^^ entire range of likely r^jg's, the molecularhydrogen solid is definitely a conductor at this density. Mote that even at this density the molecular solid is favored over the atomichydrogen solid.
PAGE 30
24 Figure 4. The tops, Â£2^0) and Â£](!<(,), and bottoms, Â£2^'^s^ ^"^ Â£1(0) > of the first two bands, and the BornOppenheimer potential, <Â£>, K versus internuclear separation are plotted, r = 1.48. The vertical lines represent the energies calculated by Monte Carlo and their a's. Since e2^^^s'^ first dips below ej(i<(,) for an r^g that is less than the r^g at the minimum of ^, the molecularhydrogen solid is a conductor at this density also.
PAGE 31
25 2 1 o CO ?^ 1 2 Figure 5, The tops, ep(O) and Â£,(](), and bottoms, e^i^^) and ei(0), of the first two bands, and the BornOppenheimer potential, , versus r, energies calculated by Monte Carlo and their o's. not dip below srCkp) until r^o is larger than the r.o for the minimum of ^, the molecularhydrogen solid will not become a conductor until the temperature is raised slightly above 0K. ,n are plotted, r = 1.97. The vertical lines represent the Since e^{k ) does 2 s
PAGE 32
26 i 4) (13 (/) s~ a >, 1 T Â— r eAk), r = 2.45> ei(0), r^ = 3.45 Â£^(0), r^ = 2.46 J L ^AB ^^B^ Figure 6. The top, e^(k^), and bottoms, Â£^(0), of the first band for two densities, r^ = 2.46 and r = 3.45, versus internuclear separation are plotted. The vertical lines represent the energies calculated by Monte Carlo and their o's.
PAGE 33
27 2 1 o +> (13 s1 AB 3 (a. and ej (0) Figure 7. The tops, Eo^^) and e](k ), and bottoms, ^n^^s of the first two bands and the BornOppenheimer potential, , versus internuclear separation are plotted, r = 2.96. The vertical lines represent the energies calculated by Monte Carlo and their a's. Since eo(k ) does not dip below eil^p) until an rnp well beyond the r.g for the minimum of , the molecularhydrogen solid will not be a conductor until a substantial temperature is reached.
PAGE 34
28 o fd cn s> Figure 8. The first two Kolos and Wolniewicz BornOppenheimer potentials (solid lines) are compared with those of this work. The dots represent the energies calculated by Monte Carlo and their o's. As expected, these potentials did not depend on k.
PAGE 35
29 s o CO CD s0) o 1 ^ (^b) Figure 9. Comparison of this work's BornOppenheimer potential versus density, e. and e.., with the atomic potentials of Bellemans and De Leener [19], triangles; Wigner and Huntington [4], and Carr [16], squares; Calais [18], circles, e. = this work's atomicsolid potential for a bodycentered cubic lattice, e.. = this work's molecularsolid potential for a simple cubic lattice. If e^^ is extrapolated to smaller r it will inter; atomic transition, r^, it will intersect e^ at r^ = Lla^; this signals the molecular
PAGE 36
30 1 o +> (/) en s_ a Figure 10. The band gap, 6 = ^2'^^s^ ^1^'^c^' ^^^^'^^ density is plotted, The insulatorconductor transition occurs when 6 = 0. The arrov/s indicate the extreme limits of this transition, r =1.35 and r = 2. s s The vertical lines represent the largest possible range of energies that 6 could have at each density. The information for this figure came from Figures 3, 4, 5, and 7.
PAGE 37
31 to Â£o(k ). At this time the volume inside the first zone should continue c s to be filled in the same manner while the volume in the second zone should be filled by occupying those parts outside of the first but inside of the second zone. (Each time a zone boundary is encountered, the above procedure would be followed.) The process is continued until a volume in kspace equal to the volume of a Brillouin zone has been filled. When the procedure is completed, all Nelectrons of the solid will have been placed in states.
PAGE 38
32 Results When (4) was solved for several k's at each r,n, several ling's at each r and several r 's for the first three energy bands, the following results were obtained. The first two BornOppenheimer curves of Kolos and Wolniewicz [28] were reproduced for H^ with an r of 10, Sag, see Figure 8. The top and bottom of the first energy band were found to agree, as expected. As r was decreased, the bottom of the band decreased in energy while the top of the band increased, for all three bands. From r = lO.Sap to r^ = 1.23an the full range of densities tested, there were no signs of a molecularatomic transition at 0K. However, if the data are extrapolated to smaller r the transition would be predicted to occur at r = 1.1, see Figure 9. The first signs of conductivity, the overlapping of the first two bands at the equilibrium rn, appeared at r = 1.97an. By the time r was 1.23an, there was no doubt that the solid was a S D conductor at 0K, see Figure 10. As seen in Figure 5, a very interesting phenomenom occurs at r = 1.97an, where the insulatorconductor transition begins to appear. A significant broadening of the BornOppenheimer potential takes place, the kaveraged energy is between 1.17 Ry/atom and 1.15 Ry/atom for r.g between 1.25an and 2.20an. The fact that the nuclei at this, and only this, density are free to move 27% of their maximum separation possible and only change energy by 16% of the well depth, is bound to have thermodynamic consequences beyond those normally associated with insulatorconductor transitions.
PAGE 39
33 Conclusion In this work it was found that the molecularatomic and the insulatorconductor transitions did not occur at the same density. The data indicate that the insulatorconductor transition most likely occurs between r^ = 1.7ag and r = 19ag, with extreme limits of r_ = 1.35anand r^ = 2.0an (c.f. experiment: r = 1.3an to 1.4 J D s D s D ao [21]). Extrapolation of the date indicates that the molecularatomic transition occurs at r^ = l.lag. Finally, the distinct widening of the minimum of the kaveraged BornOppenheimer potential at r = 1.97 is bound to have thermodynamical consequences beyond those normally associated with insulatorconductor transitions. In this work it was also successfully demonstrated that 2electron wave functions could be used in band theory. This work also opened the door to further possibilities. The BornOppenheimer potentials can be used to determine the nuclear motion and quantum solid effects. The data can be used to find the Slater sum, from which the thermodynamics of bulk hydrogen can be determined. The kind of "atoms" in the unit cell can be altered so that new systems could be explored. Even the number of particles per unit cell could be increased.
PAGE 40
34 List Of References 1. N. W. Ashcroft, Phys. Rev. Lett. 21, 1748 (1968). 2. M. D. Whitmore, J. P. Carbotte, and R. C. Shukla, Can. J. Phys. 57, 1185 (1979). 3. S. Mitton, "The Cambridge Encyclopaedia of Astronomy." Crown, New York, 1977. 4. E. Wigner and H. B. Huntington, J. Chem. Phys, 3_> 764 (1935). 5. A. K. McMahan, "High Pressure and Low Temperature Physics." pp. 2141, Plenum, New York, 1977. 6. D. E. Ramaker, L. Kumar, and F. E. Harris, Phys. Rev. Lett. 34, 812 (1975). "~ 7. R. D. Etters, R. Danilowicz, and W. England, Phys. Rev. A 12, 2199 (1975). 8. E. Ostgaard, Z. Physik 252, 95 (1972). 9. H. W. Myron, M. Y. Boon, and F. M. Mueller, Phys. Rev. B 118, 3810 (1978). 10. W. Andreoni, Phys. Rev. B 14, 4247 (1976). 11. M. Ross and A. K. McMahan, Phys. Rev. B 13, 5154 (1976). 12. J, Oddershede, L. Kumar, and H. J. Monkhorst, Int. J. Quantum Chem. Symp. No. 8, 447 (1974). 13. H. J. Monkhorst and J. Oddershede, Phys. Rev. Lett. 30, 797 (1973). 14. F. E. Harris, L. Kumar, and H. J, Monkhorst, Phys. Rev. B 7, 2850 (1973). 15. T. Schneider, Helv. Phys. Acta. 42, 957 (1969). 16. W. J. Carr, Phys. Rev. 128, 120 (1962). 17. G. A. Neece, F. J. Rogers, and W. G. Hoover, J. Comput. Phys. 7_, 621 (1971). 18. J.L. Calais, Arkiv Fysik 29, 255 (1965).
PAGE 41
35 19. A. Bellemans and M. De Leaner, Phys. Rev. Lett. _6, 603 (1961). 20. A. L. Ruoff, "High Pressure and Low Temperature Physics." pp. 2141, Plenum, New York, 1977. 21. P. S. Hawke, T. J. Burgess, D. E. Duerre, J. G. Huebel R. N. Keeler. H. Klapper, and W. C. Wallace, Phys. Rev. Lett. 41. ^94 (1978). 22. Experiments are being conducted at Cornell University and Los Alamos Scientific Labroratory. 23. H. M. James and A. S. Coolidge, J. Chem. Phys. 1, 825 (1933). 24. J. Kolwalik and M. R. Osborne, "Methods for Unconstrained Optimization Problems." American Elsevier, New York, 1968 25. R. L. Coldwell and M. A. Pokrant, J. Computational Phys. 20, 365 (1976); the minimization routine was developed by R. L. Coldwell. 26. W. A. Harrison, "Solid State Theory." Dover, New York, 1979. 27. J. M. Ziman, "Principles of the Theory of Solids." Cambridge Univ. Press, New York, 1972. 28. W. Kolos and L. Wolniewicz, J. Chem. Phys. 43, 2429 (1965).
PAGE 42
36 PART II: A QUANTUM MECHANICAL TEST OF DIOPHANTINE METHODS The wave functions in Part I were calculated by varying the twentyfour parameters to minimize the Monte Carlo error in estimating the energy, using roughly 1,000 functional evaluations. However, to reduce the Monte Carlo error in estimating the energy for the final set of twentyfour parameters, 500,000 functional evaluations were required. Thus, 499,000 functional evaluations are used not in improving the wave function but only in reducing the Monte Carlo error of the energy estimate to acceptable levels. Obviously, it would be a big help if the required number of functional evaluations could be reduced. In this part of the t h e s is, a method is presented which can decrease the required number of functional evaluations by a factor of four. In particular, a numerical experiment was conducted to compare the biased selection method [1] of evaluating the expectation value of the Hamiltonian for a dihydrogen molecule using the 11term JamesandCoolidge wave function, ^Pjp, [2] with a blend of the biased selection and Haselgrove's method [3]. In its essence, the blend amounted to using Haselgrove's method of generating numbers in place of the random number generator used by the biased selection method (the latter method is a Monte Carlo scheme). The results were mixed. The "blend" gave an error proportional to either i;j or Â— dependina on whether N was, respectively, less than or greater than 50,000
PAGE 43
,37 (N = the number of functional evaluations). The "blend" also gave an error that went from being larger than the biased selection error (for N < 10,000) to being only half as large for N > 50,000. These numbers should be compared with the 10,000 points used in finding an optimum set of Hasel grove's parameters. It is reasonable to expect that if 100,000 points were used in the optimization of Haselgrove's parameters that the above results would be found with 50,000 replaced by 600,000. At any rate, a faster rate of convergence than j^ was found for a large range of accessible N.
PAGE 44
38 Hasel grove's Method Haselgrove's method was investigated to see if it could be used to give a more rapid rate of convergence in estimating the multidimensional integrals encountered in quantum mechanics than the Monte Carlo methods in current use. His method promises an N""" (with i ^ 1) convergence for integrands which are fourier expandable and whose derivatives satisfy certain continuity conditions. The method uses M S^ = f(o) + 2 y c fCl2LmaJ, ..., 2LmaJ) (1) to estimate the integral (2L)" /... /f(x)dx, call it I, when f(x) is an even function (K e the dimensionality, f(x) e f (x, ..., x,,) dx E dx, ... dX.). The c 's (which are functions of N, m, and i) determine the rate of convergence, the value of i The a.'s are irrational numbers, which must be found and which admit no solutions K to 1 + )^ M.a. when the M.'s are integers. i=l ^ ^ ^ K For i equal to two the integrand W(x) e tt (1 Ix.l)^, where _K 1 1 ^^^ .K 1 = 1 2 /.../ W(x)dx = 3 gives the worst possible error, e (to within 1 1 w a multiplicative constant that depends on the specific integral being evaluated but not on N) of all the integrands satisfying Haselgrove's conditions. The key to finding the a.'s lies in minimizing e This 1 ^ w was done by selecting K random numbers (between and 1) and varying them by a simplex method [4] to minimize e Tables I and II give w ^ some a's and corresponding e s as functions of M. ^ w
PAGE 45
39 Table I. This table contains the a^,^ found by minimizing e over Kdimensions, using either my method (second subscript M) of finding {2ma. } or Haselgrove's method (second subscript H). a optimized over 1,000 points ^5M .45260537885008 .47561868232093 .57787664217299 .49493808300031 .21929292303371 a optimized over 10,000 points "5M .45241479419455 .47565271662758 .5782437781132 .49509111688911 .21929571968924 "5H .957286852634 .86730302040274 .097247969465186 .31301711424498 .48486758446139 a optimized over 10,000 points "6M ,43657287149021 ,59181442054365 ,050243215161086 ,84377789758502 ,38121353509683 ,78500128467888 a 8M ,73778770979932 ,083137623602196 ,84742517965315 ,88984806744175 ,80256318005657 ,27985939980508 ,67436009334498 ,5303179412788
PAGE 46
40 Table II. This table gives the absolute value of the difference, multiplied by 10 between W and Haselgrove's estimate of W, using my method of finding i2ma.} (em^^) and Haselgrove's (euu). ^WH (ar,. optimized over 1,000 points) 1000 9774 5570 3000 2338 1876 7000 887 596 15000 237 249 31000 125 118 59000 58 47 {acM optimized over 10,000 points) 1000 7818 7404 3000 3173 2225 7000 1289 541 15000 299 191 31000 112 72 59000 42 35 (arn optimized over 10,000 points) 1000 5980 5434 3000 1326 1179 7000 405 447 15000 143 223 31000 63 95 59000 38 43 (a from Haselgrove's paper) 1000 5933 5716 3000 1592 1808 7000 455 1086 15000 282 732 31000 178 462 59000 114 241 It can be seen that e^n and e^,., have been substantially lowered from those given by Haselgrove's original a. Also note that Â£,,n and e^., differ by less as N is increased, regardless of which method was used to find {2ma}. Furthermore, the e,, 's for 'S.^^, and dn.,, optimized over the same number of points, differ by less as M is increased.
PAGE 47
41 N+l Inl 1 If c is equal to (\!^\\ or ^inrr then an i equal to 2 or 1 should be found. The main reason an i equal to 1, instead of 2, might be used is because Haselgrove's method places less stringent demands on the integrand for smaller i's. In the experiment, both i's were found to give the same convergence. Eventually, 2ma. will fall outside the integration region as m is increased. Two distinct methods vjere used to reduce 2ma. to the ma. integration region: Haselgrove's used 2[2^ for each x. ([x] is the signed distance to the nearest integer, e.g., [1.9] = .1 and [2.1] = +.1), my way used the periodic equivalent of 2ma. lying within the interval from 1 to +1. (In either case, denote this reduced number by {2ma.}.) The a that is found by varying an initial a to minimize e, will depend on the method used to find {2ma. }. If one method is used to find {2ma. } for the minimization and the other method is used to find {2ma.} in estimating I, then reasonable, although not optimal results, will be found. In the experiment, both methods gave equivalent error estimates. The a. 's used in the experiment had a 14digit computer representation. Since the a's are represented by a computer, they are not irrational. This fact will be a problem when N is roughly 10'^ [5],
PAGE 48
42 Haselqrove's Method Used by Itself Hasel grove's method was used to estimate the known values of some definite integrals. His a was used to estimate the value of the 5dimensional integral / ... / e ^ ^ "^ ^dxTdx^dx^dx^dXr, and his 12 3 4 5 results were reproduced exactly. The method was then tested on a problem that required a change of variables. The expectation value of the Hamiltonian for a helium zs 27 atom Â— using = e z = j^, s = r, + r2, r. = nuclearelectron separation was estimated. This expectation value can be determined by evaluating CO s u L = 27TVds/du/dt(8su s' + t')e"^^^ s u ^ 2 TA^ rA.. rAAÂ„ZSr,,/^2 4.2\rii; J. ^ M E 27TVds/du/dt e ^^{u(s' t')[^ + ] + 4us }, (2) ss u s and CO s u Norm = 27rVds/du/dt u(s' t')e"^^^ where "^^ = gj [5]. Haselgrove's method, since it is designed to sample the 3dimensional space uniformly, will sample the important region of the integrand with a yery low probability. In order to increase the probability of sampling the important regions, a change of variables, which amounts to importance sampling [5], can be employed s' = Â£n(l s), f = ^ u' = I
PAGE 49
43 These transformations alter Eqs. (2) to 111 L = 2TTVds/dt/du{u(l s)'^^"'[8u 1 + uH'] Â£nMl s)} M = Z^^zVdsdtduld s)^^"^u'[l uH^+r4T Â— rl^n'Cl s)} (3) Norm = 27TVdsdtdu{u2(l x)^^"^[l uH^J^n^d s)}. In Eqs. (2) and (3), M, L, and Norm are the expectation of V^ the potential, and 1, respectively. Table III gives some results. This problem has an interesting feature: the variable transformation introduces singularities into the first derivative of the integrand. Even with this violation of Haselgrove's conditions, the results in Table III seem to give a better convergence than ^. vH The results with helium were good enough to warrant moving on to a problem which is more difficult the evaluation of the expectation value of the Hamiltonian for the dihydrogen molecule, using ^jp. First, Haselgrove's method alone was used, with the results given in Table IV. The lack of any signs of a 772 or even a jt, convergence was not unexpected for the M used (the problem, in helium, of sampling the important region frequently enough is made more severe here due to the increased dimensionality). Next, Haselgrove's and the biased selection methods were combined with the hope of gaining the advantages of both. Biased selection would be a way of sampling the important regions more frequently, taking the place of the importance sampling used for helium. By replacing the random number generator with Haselgrove's method of generating points in a multidimensional space, we hoped to gain the advantages of sampling the important
PAGE 50
44 Table III. Energy bounds, /'i'H'i'dr,dr2//'i'^dr,dr2, for a helium atom were calculated using Haselgrove's S and two different methods. N Haselgrove's Method 5.597346 .002644 My 5 Method 700216 4E^W 1000 .003694 3000 5200 611 5 695513 570 5000 026 319 516 228 7000 189 182 350 137 9000 358 114 354 91 11000 338 68 29 68 13000 305 68 15 46 15000 278 68 23 46 17000 285 46 27 46 19000 298 46 15 23 For comparison the exact value of the integrand, which is not the exact ground state energy, is 5.695313. All answers are in Rydbergs. The error estimates, 4EÂ£,,, were calculated by assuming that the error in each integral in Eq. (3) was e<,, and then by taking small e,, limit.
PAGE 51
45 Table IV. This table gives estimates of the expectation value of the hydrogen molecule, using James and Coolidge's 11parameter wavefunction and Haselgrove's method without the aid of biased selection (all i energies are in Rydbergs and the correct answer is 02,34609). ,, ,, N H H 1000 2.31939 2.49532 3000 2.22855 2.22018 7000 2.22166 2.23877 15000 2.23085 2.24075 31000 2.23844 2.25354 59000 2.25539 2.23981 The second column gives the results when Haselgrove's method was used to find {2ma. }. The third column gives the results when my method was used to find {2ma^. }. Because it can be difficult to sample the important regions of a 6dimensional integrand, the poor results were not unexpected. Compare these results with the results of Table V where the help of biased selection was employed. The ar.y. from Table I was used. ^^
PAGE 52
46 regions with more uniformity [5], [7]. In other words, an error estimate which was proportional to rr rather than = seemed possible.
PAGE 53
47 Haselgrove's Method Used in Combination with Biased Selection The evaluation of E = /^ipH^,Â„dx//Y,pdx was carried out over an integration region which was a box of length L = 20aD5 (L = L = 22an, an = Bohr radius) in the x (y, z) direction. This box is large enough to insure the desired accuracy. The two nuclei are fixed on the zaxis with the point halfway between them coincident with the center of the box. If the biased selection method were going to be used to estimate E, then the following could be done. Select the coordinates of the first electron in four ways and those of the second in five ways (four ways being the same as for the first electron). The estimator for E would be BN=ll^f(J,)/Sa.), (4) where the relative probability of selecting a 5dimensional x. w(x.), exactly corrects for the nonuniform selection of the x.'s. If each of the twenty ways (four ways for the first electron times five ways for the second) of selecting an x. is used with equal probability then w(xj., x^.) = 4[w(x^^. ) + W2(x^^. ) + W3(x^.) + W4(x].)] x ^w^(x2^. ) + W2(x2^. ) + W3(x2.j) + yi/^{y^2^) + vt^ix^.U. These Wi's, etc., are determined by the methods used to select the
PAGE 54
.48 electrons' coordinates: = '^i 4' y = h^r ^ h^z w,(x^) = 1 (R. = random number uniformly selected between and 1 a total of six R. 's are needed for each x. ir 1 Eq. (4)) ^C = ^l\' ^ = 2R21, = Z.R3 W2(x^) = V/(2uMjr^) (r = the electron's distance from the center. u = the cosine of the angle between r and the yaxis, (j) = the angle between the projection of r in the xy plane and the Xaxis, V = L^L^L^, \ is a constant) U^ = 2R2 1, (f)^ = 2^R3 r is selected by requiring 1 rgi 8 / g(y)dy//g(y)dy = R, 1 1 Wo(x.) = Vg(r^.) 4TrrL/g(y)dy AIq w.(x.) = Vg(rg.) ^ ^ 8 4iTr2 /g(y)dy BIq 1 1 (a = A or B, i 1 or 2, r = an the distance betweer 1 nucleus a and electron i)
PAGE 55
49 W5(X2) = V/(27tM^R^) (rr = the electronelectron distance, Mp e a constant). The function g(r) was created by evaluating re at 64 points and connecting these values with line segments. The estimator that is used to estimate E^ is 20 N I I ^Jc(^m)H^Jc(^m)^(^m) Â•^N 20 N ^ ,1 I J. ^oct^J/^t^.) where w(x) = ^fJ^ir^) + W2(r^) + W3(r,) + W4(r^)] x 5{w^(r2) + W2(r2) + W3(r2) + W4(r2) + Wg(r2)] (r. = electron's coordinates, x represents an r, rp pair). The i in Eq. (5) signifies which of the twenty ways (four ways for the first electron times five ways for the second) is to be used to generate x. from {2ma}, for a given m. (This same {2ma} is used to generate a total of twenty x. 's.) This procedure is used (instead of using an 8dimensional a and using the two extra a.'s to randomly select which of the four (five) ways would be used to select the first (second) electron's coordinates, as can be done in the usual biased selection) to insure that x. will be close to x. if {2ma} is close to {2na}. Note that for c., to be fairly compared with a Monte Carlo estimate that the Monte Carlo estimate would have to use 20N functional evaluations.
PAGE 56
50 Table V. Energy bounds (in Rydbergs) for an H^ molecule calculated using the 11 term JamesandCoolidge wave function and using a^^ from Table I to generate the "random numbers" that were used to determine the electron positions. 4EÂ£:, N Energy Bound 2.34684 """W 1000 .006982 3000 560 1736 7000 585 619 15000 596 241 31000 501 105 57000 607 57 The blend of Hasel grove's and the biased selection method was used, The exact answer, to six digits, for the given wave function, is 2.34609 Ryd [2].
PAGE 57
51 Further Tests of the Blend of HaselGrove's and the Biased Selection Methods One must be careful to insure that the method of generating the x^'s is not creating correlations between the x. 's, A commonly used random number generator, furnished by IBM with their FORTRAN compiler, was found to have this defect [8]; the random number generator used for the calculation herein did not. Therefore, Haselgrove's method was used to calculate the nuclearelectron and, especially, the electronelectron correlation functions. These functions were then compared with the correlation functions generated by using biased selection. The functions found by the two methods agreed, indicating that the configurations generated from a did not have any correlations that were not due to particleparticle interactions. The correlation functions can be compared in Figures 1 and 2. Unwanted correlations were also tested for by seeing if the configurations generated by using Haselgrove's scheme, without the help of biased selection, were uniformly distributed throughout the volume. This was done by dividing the 6dimensional integration volume into 729 equalvolume hypercubes Â— 1 .e. three equally spaced divisions were made in each dimensionÂ— and noting the distribution of points within them. Figure 3 illustrates that the number of points per hypercube was normally distributed and that this distribution has less deviation than uniformly and randomly distributed points have, as expected.
PAGE 58
52 .04 ,02 ^ee ^^b) Figure 1. Comparison of the biased selection method () and the blend of biased selection and Haselgrove's method (Â•) for calculating the electronelectron correlation function. The biased selection method used 40,000 points and the 'blend' used the equivalent of 200,000 Monte Carlo points. The excellent agreement demonstrates that the points generated by the 'blend' did not introduce any nonphysical electronelectron correlations.
PAGE 59
53 ^ne ^^B^ Figure 2. Comparison of the biased selection method (]) and the blend of biased selection and Haselgrove's method (Â•) for calculating the nuclearelectron correlation function. The biased selection method used 40,000 points and the 'blend' used the equivalent of 200,000 Monte Carlo points. The excellent agreement demonstrates that the points generated by the 'blend' did not introduce any nonphysical nuclearelectron correlations. i*3r%^^^''<;:^i:~^^^'^
PAGE 60
54 ,4 .2 Figure 3. N points were placed in a 6dimensional volume that was divided into 729 equal cells, using Haselgrove's method. The number of points falling in each cell was counted, and the mean and standard deviation of these 729 numbers were calculated. The fraction of cells containing a number of points which was a given number of standard deviations av/ay from the mean was plotted. The squares are the results of a 10,000 point run.: mean = 13.7, a = 2.53. The triangles are the results of a 20,000 point run: mean = 27.4, a = 3.73. A normal curve is drawn for comparison. Also for comparison, a purely random set of 10,000 (20,000) points would be expected to have a mean of 13.7 (27.4) and a of 3.70 (5.23).
PAGE 61
55 Monte Carlo Estimation of the Error Since e at best, can only determine an error estimate for Eqs. w (1), (4), and (5) to within an undetermined multiplicative constant, which is independent of N, some trial anderror searching would be required to find this constant. But there is a very real danger in this type of an approach. If the constant, B, is large enough to ensure that rr is larger than the difference between c^, and the correct answer for some Nn, then the error might be overestimated for smaller N and, at the same time, be underestimated for larger N. Inthe case of estimating Ep, the 5decimal place accuracy quoted by James and Coolidge was approximately achieved for N = 59,000. Therefore, the B that would be correct here could not be tested by increasing N. In order to obtain a reliable error estimate, a Monte Carlo procedure was used. In place of Eq. (1), the following estimator N ^^(Ri) = I V(2ma^ + Ri^U Â•Â•Â• > I2ma^+ f^Ki ^ m=l .> can be used without violating any of Hasel grove's conditions (R. is a set of K random numbers selected uniformly from the interval between and 1). By selecting M ^.'s and estimate can be made of the average and the standard deviation of the distribution of S^(t. )'s. The standard deviation can be used as an error estimate for any given Sj^,
PAGE 62
56 An estimate of the error in estimating e is fi To prevent the overhead in finding Â£ and 6 from going to waste, the foil owing refined estimates of I and its error can be used: with the error in the above estimate given by Â£ Figure 4 displays some results. It is interesting to note that for N less than or on the order of 10,000 the number of configurations used to optimize a that the convergence appears to be rr, but for N yery much larger than 10,000 the convergence becomes p This suggests the possibility of extending /N the TT convergence to larger N by increasing the number of points used >Â• to optimize a.
PAGE 63
I" 57 8000 4000 2000 1000 400 200 100 40 o 20 10 Figure 4. Comparison of the error using biased selection (straight line prediction from 40,000 points) and the error using the blend of biased selection and Haselgrove's method (vertical lines). The size of the vertical lines represents the possible error, 6, in estimating the error, z.
PAGE 64
58 List of References 1 see Appendix A. 2. H.M.James and A. S. Coolidge, J. Chem. Phys. 1, 825 (1933). 3. C. B. Haselgrove, Math of Computations J5_, 323 (1961). 4. J, Kolwalik and M. R. Osborne, "Methods for Unconstrained Optimization Problems." American Elsevier, New York, 1968. 5. J. M. Hammersley and D. C. Handscomb, "Monte Carlo Methods." Chapman and Hall, London, 1979. 6. B. L. Moiseiwitsch, "Variational Principles." pp. 172183, Vol. XX, Interscience Pub., 1965. 7. H. Niederreiter, "Diophantine Approximation and Its Applications." pp. 132194, Academic Press, New York and London, 1973. 8. R. L. Coldwell, J. Computational Phys. M, 223 (1974). <;=^'iC,o's:J3:^j'5j'^f"
PAGE 65
CONCLUSIONS OF DISSERTATION The calculations of Part I required an enormous amount of computer time, roughly 1200 hours of Amdahl 470 cpu time. While experience would allow these calculations to be redone with, perhaps, onehalf to onethird of the computer time, an enormous amount of computer time would still be needed. This time was mainly used in the minimization of a to find the wave functions. Less than 57o of the time was used to find the energy of the final wave functions for the equilibrium separation with a Monte Carlo error that was less than chemical accuracy. For the required amount of computer time to be reduced by more than the factor of two to three, three things would have to be done, separately or jointly. Firstly, the minimization would have to be made faster. Frequently, more than 500 calls had to be made to the simplex algorithm to find the best wavefunction. Dr. Coldwell is making continual progress along this line. Secondly, the number of Monte Carlo points required to evaluate the integrals accurately enough to allow meaningful minimizations would have to be reduced. Part II was a first step and suggested further steps to take along this line. Thirdly, better wave functions would have to be used. The form of b(r) was not flexible enough to allow the wave function for the first energy band to have the same symmetry as the wave function for the second band: b(r) should have been IcpB (G e reciprocal lattice translation vector). Furthermore, the form of 59
PAGE 66
60 ur;(r) became more inappropriate as r^ 0, more and more Monte Carlo points were required to minimize o effectively. This work showed that it is possible to use Monte Carlo and 2electron wave functions to calculate the properties of bulk matter. It also gave encouraging signs, in Part II, that the l/i^N convergence of Monte Carlo can be beat in the evaluation of multidimensional quantummechanical expectation values.
PAGE 67
APPENDIX A THE BIASED SELECTION METHOD INTRODUCTION This appendix presents in one place for the first time detailed proofs, statements, and examples of the biased selection method of evaluating integrals and expectation values. Biased selection is, as the name suggests, a method of nonuniformly sampling an integrand, f(x), in order to estimate /f(x)dx. It is an extension of importance sampling, as given by Hammersley and Handscomb in reference [1], in the sense that it uses sums instead of integrals in order to exactly correct for the bias introduced by using a weight function, g(x). In addition, biased selection allows one to use several simple weight functions simultaneously, one for each of several important regions of the integrand, in place of a single complicated weight function. Biased selection is harder to code than the Metropolis scheme, [1], [2], [3], but more flexible. Coldwell and Lowther used it in their calculations on the dihelium system [4]. A weight function, g(x) > 0, that mimics f(x) can be used to lower the standard deviation in the Monte Carlo estimate of /f(x)dx. For example, let y = /^g(z)dz//g(z)dz so that 61
PAGE 68
62 (all integrals, from now on, are from to 1 unless explicitly stated otherwise) where w(x) = g(x)//g(z)dz. The estimate of the integral is then irry__^!j_ where the y.'s are ^"w(x(yi)) picked uniformly at random. When the exact value of /g(x)dx is used the above is importance sampling. It will be shown that biased selection gives similar results but does not require the exact evaluation of any integrals in order to exactly correct for the bias introduced by the nonuniform selection of the x's. This appendix will also show how several different w's can be averaged and used in place of a single w. This averaging allows the combination of several simple weight functions, rather than one complicated weight function, to sample the space adequately. The methods detailed in this appendix often select one point from M points and repeat this N times to get the integral estimates. As Mx^ the biased selection method becomes importance sampling, owing to the fact that certain sums become integrals. It must be emphasized, and will be proved, that for M fixed, e.g., for M = 3 the method converges on the correct answer almost surely as N^=. Unless specifically indicated, all sums are from 1 to
PAGE 69
63 I. ERROR ESTIMATION The simple Monte Carlo estimate for a multidimensional integral is where the x.'s are chosen uniformly at random. As Nxo, Ij^ ^ /f(x)dx. Since the x.'s are random, sampling theory can be used to derive (for N large [5], [6]) an error estimate for I^ a = 4l(f(x.) I^VA^, (2) and an error estimate for Na^ & =' [\i NVWv^, (3) 4 where u e Ij (f(x) L,)"* [1], [2]. For these errors to assume their Standard meanings (e.g. Ilj. ff(x)dt\
PAGE 70
64 formly at random and then calculating x, see introduction, in the case of biased selection the techniques of section III can be used). Biased selection requires fewer points than simple Monte Carlo does, in order to obtain accurate estimates of the integral and its error. It reduces a by reducing its proportionality constant, not its inherent 1//R" dependence. Restrictions must be placed on the integrand for a and 5 to be used as error estimates and maintain their common meanings. In particular, o and y should exist and be bounded, and a should tend to zero as N^. It is suprisingly easy to violate these conditions. Consider the integral /x'^dx, it exists but ct and y are infinite. Numerically, this shows up as a a that decreases as N increases, and then rises when a random point is chosen which comes much closer to the origin than all previous points; and as a 6 which becomes much larger than a as N is increased. If the proper weight function, e.g. g(x) = ^x 3 is chosen, the integral and a can be interpreted in the usual way. in this case, the estimator is 775" 7735 and o,v~r, and 6,/^, behave properly. Suppose that nothing more were known about the integrand than that the largest contributions come from the origin. Observe what happens if the weight function favors the origin too much. For example, if g(x) = y^o^"''^, the integral becomes //^^ (750^^^^"^^"'^^ which reduces to 100/y'*^dy. The problem at x = has become a problem at X = 1, y = 1, which, for good numerical results, still requires N :A ^CJ;,i/titij ,i^ J,.^i
PAGE 71
55 to be very large (this is now necessary in order to sample the region near x ~ 1 adequately, due to the low probability of placing a point there) For a more useful example, consider evaluating the expectation value of a quantum mechanical Hamiltonian, H. Let E^ be the eigenvalue of H, and ^.(x) be an approximate (normalized) eigenfunction. Then H (x) = E(x)^,(x). Note that if \(x) were an exact eigenfunction that E(x) would be constant. However, for our purposes, suppose that ^4(x) is not an exact eigenfunction. We will assume that it is very good for 'i'^(x) large and very bad for ^^(x) small. If we use importance sampling and let the sampling function be g(x) = "^Ux) + 5 then the integral = /[ (x(y))(H EjW,(x(y)) Lt!^^}^]lll_{'ldy (4) ^ ^ ^ ^ V(^J(z) + 6)dz^ can be simplified and the estimator written as = 1(1 + 6) l!iV_^. (5) ^ 1 + 6/r^(x.) It might appear that the best sampling function would be '^Ux) (i.e., 6=0), and, in fact, the Metropolis scheme uses this weight function. But, since most of the error in H,(x) is assumed to be for small values of 'l'^(x), say for ^^(x) less than a tenth of its maximum value (^j^^), then E(x) Eq will be the largest there. The error term will then be dominated by an area that is unlikely to be sampled (because ^^(x) is small). This means that a large N is needed to lower the error in because most of the points contribute little to the error while a very few make large contributions, since w(x.) is very small for them. If 6 were kept finite (say 6 = I'^^u) instead of being set to zero, then the "misses" would be lowered in size to
PAGE 72
66 ^^' instead of E(x) E^. This is, of course, exactly compen1 + 6/j(x) ^ sated for by the fact that there are now more of these "misses". Notice, in this case, that the term contributing to the error in is approximately halved. Each "miss" is cut by a factor of 4 in its contribution to a, but there are twice as many of them. So, by sampling the "unimportant" regions more (i.e., no region is sampled with a probability that is less than, in this case, .l^Â„),a can be reduced. It is interesting to observe that a was reduced by using a I seemingly worse weight function, and that a large number of points appear to be wasted because (E(x) Eg)/(1 + 5/('i'^)) is (as j^) i zero for the most part.
PAGE 73
67 II. BIASED SELECTION In biased selection, a positive definite weight function, g(x), is used to distribute the x. 's, with a probability Pd(X:j) (or probability density, as the case may be) which must be computable (explicit schemes for doing this are detailed in Section III), The bias is compensated for by dividing each f(x.) in Eq, (1) by the relative probability of choosing x. w(x.) E Pg(?.)/P^J(X.), (6) P (x.) is the unbiased (uniform) probability of choosing x.. One should note, as the discussion develops, that w(x.) can be a multivalued function. This is because there can be more than one way to select any given x. and each way can have a unique w(x.) associated with it. The biased selection estimator is 1 r f(x. ) /^N To prove that Eq. (7) gives /f(x)dx as N>co, the following can be done. Divide the integration region into T equalvolume hypercubes (label them T^ where a = 1, 2, ..., T) and select N x. 'S, each with relative probability w(x^. ), from it, NT [4]. Now, break up the sum into two sums one over the T 's and the second over the points within a a given T, a In=iI l'f(^)/w(t). '^ '^a=li=l ^ ^
PAGE 74
68 To derive Eq. (8), two cases need to be considered. The case for the multivalued w(x.)'s used herein is discussed in Section III, The case for singlevalued w(x.)'s follows. As T becomes infinite, the volume of a given T becomes vani shingly small. Therefore, the ^ > ^ ^ >. ffy^ ) fix) continuity of f(x) and w(x) requires that i J ^Â„^ ^in z w(x^. ) w(Xj) 1 and X. within a given T,,^ Therefore, combining the n. equivalent J a a terms within each T gives a n, '^IL (8) 'N TH N ^ "a r^v a=l w(x^} a Since NT, the law of large numbers can be used to show that _J. ^ Pd(x,) 1 N ^ almost surely, a.s,, as N^. Also note that PmCx^) = j. These two facts and Eq. (6) imply n^ Â— ^Â— ^ N(t)w(x,); which when substituted into a I a Eq. (8) gives hlS^nl ^(T)^(^a)^ (almost surely) a=l w(Xa) a=l This is the definition of the Riemann integral in the limit of infinite T Note that if w(x) were equal to unity that the above proof reduces to a proof for simple Monte Carlo. ^ The results remain valid if f(x) and g(x) have a countably infinite set of regions over which they are continuous, but the explanations are too detailed for this paper's purpose. '."^5^i;^^'St^s;^&:.^=5S^^' S' ^^j?v>MJ=^'Â£^it5i^ ssi"
PAGE 75
69 The closer the weight function can be made to f(x), the more constant f(x)/w(x) will be; and, consequently, the smaller a will be (approaching zero for nonnegative f(x) and f(x)/w(x) absolutely constant). While w(x) = f (x) //f (x) dx does minimize a, to evaluate w(x) one would have to know the integral that, by hypothesis, is not' known. It should be cautioned that a bad weight function can greatly increase the number of points before Eqs. (2) and (3) accurately indicate how far Ij^ is from /f(x)dx. For example, if w concentrates the points in "unimportant" regions, the integral effectively becomes the integral of f(x.) over these "unimportant" regions. The important regions will not be sampled for small N because of their small probability of being selected. Thus, the averages and error estimates are only for a part of the total integral being sought. This results in a low a which could lead to the erroneous conclusion that Ij^ was close to /f(x)dx instead of the correct conclusion that l^ was a good estimate of the integral of f over the "unimportant" regions. This misconception will not be revealed until N is large enough for the "important" regions to finally be sampled. (Note: For the important regions to eventually be sampled, in general, g(x) must never be zero.) Let us assume that we have an integrand, f(x), with the locations of some peaks known and others unknown. Then a measure of the number of points making a significant contribution to the estimate of the integral is %ff'~ii\'^\r/n!!^\ (9) w(Xt) w(xi) ,i;.ii;il.72Ei.Â£.
PAGE 76
70 The number of "effective'^ points can range from N ^^ = 1 when one '> \ ,^ value of f(x.)/w(x.) is much larger than the others (and therefore makes the major contribution to Eq. (7)), to N ^. = N, when all N values of f(x. )/w(x.) are equal. Since it has been assumed that we do not know the location of all the peaks, it is undesirable to have N^^ N. This is because such a large N^^ can only be obtained by concentrating the points where peaks of the integrand vyere thought to lie. This concentration precludes any chance of selecting a point from one of the unknown peaks. A good weight function to use is g(x) = b(x) + d, where b(x) incorporates, as much as possible, all of the known properties of the integrand and d is a constant. This way, if there exists some Xq where f(xQ) is large but was thought to be small (i.e. b(Xg) is small), then x^ still has a chance of being selected. The selection of this Xg causes N .^ to drop to approximately unity. The statistics of this run become highly suspect. However, the precipitous drop in N^^ tells us that we didn't know all of the regions in which f(x) was large, and therefore, allows us to select a better b(x), enabling us to eventually find the correct value of /f(x)dx. Notice that it was the presumably wasted points that enabled us to locate the unsuspected peak. In practice we have found that when N ^^ < .2N it pays to improve b(x), and that for M^^^ > .3N it pays to sample more of the regions where b(x) is small (increase d). The following example illustrates this. Suppose that f(x) has two equal peaks (where f(x} = 10^) that each extend over a length equal to .001 (where the total integrand is from x = to x = 1 ) and is one elsewhere. Then consider a b(x) that equals f atone peak and zero
PAGE 77
71 elsewhere. If d l,then information sampling would have roughly a chance in 1000 of sampling where b(x) = 0, and then only a chance in 1000 of that point hitting the second peak. Therefore, for N 10'* there is a 99% chance of not sampling the second peak. This results in an I^ ^ffMdx and an N .^ ~ 10\ If N is increased to 10^ there is a 73% chance of hitting the second peak and getting I^ /f(x)dx and an N .^ 4. Finally, if d is increased to 5000, then there is approximately an 83% chance of selecting a point where b(x) = 0. Then for M lO'*, roughly 8 points would hit the second peak giving Ip^ /f(x)dx and N^^^ 30. So, with fewer points than was required for d/bl, the estimates were improved. Furthermore, enough evidence has been gathered to alter b(x) so that it contains the second peak. The Metropolis scheme for evaluating expectation values, e.g. = /h(x)f(x)dx//f(x)dx, uses f(x) to sample the integrand, but does not give an estimate for /f(x)dx as biased selection does. This weight function gives the minimum error in only when h(x) = 1. The Metropolis method is easier to code than biased selection, since there is no need to calculate probabilities, biases, or weight functions, and since it is very simple to generate x.,, from x.. But each x.,^. 1 + 1 1 T + I is strongly correlated with each x^. This means that a large number of functional evaluations must be discarded before the remaining can be used to estimate the integrand. However, all of the functional evaluations in biased selection are statistically independent, giving the biased selection method an advantage when N must be small. By increasing the complexity of the Metropolis scheme it is possible to use a different weight function, g(x). The integral, /g(x)dx, generally is not evaluated in the Metropolis scheme. For this reason the Metro^5E^:^.:&dS;Jaiv^Â•:c^>^"^.:^ci^*: .j^;r,Â£5a5.*sriLi.>:'^Si
PAGE 78
72 polls scheme doesn't have the ability to use multiple weight functions simultaneously. This gives biased selection an edge in evaluating the integrals of multipeaked functions. Biased selection can use several weight functions to generate a set of points which are not as concentrated as those generated by Metropolis' method. This set of points (and associated w(x^. )'s) can be saved and used to evaluate the integrals of a wider variety of integrands providing useful correlations in the integral estimates and saving the expense of generating a new set of X. 's and w(x. ) 's. The class of functions which can be used for w(x) in Eq. (7) can be greatly extended to include, not only relative probabilities, but also weighted sums of relative probabilities. For example, a point, x, can be selected by two different schemes one which has biased and unbiased probabilities of selecting an x of Pg, (x) and Pmi(x), the second which has biased and unbiased probabilities of selecting an X of Pg2(x) 3"^ Py2(x)This can happen in cases where one scheme rarely selects points from a region that the second scheme frequently does. It will be shown that if the point is selected by one method some fraction, p^ of the time and by the second method the remainder, P2 = 1 P] of the time, then w(Xp) = p /Bly + P2 Pb2V = PlW^(Xp) + P2W2(Xp) (10) is the correct function to use in place of w(x) in Eq. (7). (The w(x ) in Eq. (10) is the relative probability of choosing x only if P..(x ) = Pu2(^p)') I* should be emphasized that the w in Eq. (6) and the w in Â•y.'Sl^ia'S^'^i^lS* ^^ ..^^;v5:^va:^3tS.^Â£<=t^^i&: ^"
PAGE 79
73 Eq. (10) represent totally different ways of selecting points. In Eq. (6) each of the points is selected using the same weight function g(x). But in Eq. (10) each of the x. 's is chosen either using g, (in which case this point, x, would be used to evaluate the contribution, f(x^^. )/w(x^^. ) to the sum) or using g^ (giving f (x^)/w(x2 ) ) but not by both. The following shows that w(x) does indeed give the correct results, First Eq. (7), with Eq. (10) used in place of Eq. (6), is broken into two sums, one is over the points chosen by w, and the second is over the points chosen vu, pjN P2N I^ =^.I^f(x^,.)/w(x^.) +.I^f(x2^.)/w(^2i)^Breaking these into a sum over the T 's and a sum over f. ( the a 1 a fj^'s are the f/s within a given T^) and then using continuity to carry out the sum over the f. 's gives hY^l HNp^w^(x^)/T + Np2W2(^^)/T]f(?^)/w(?^) aÂ— I aÂ— I What's more, this estimator is still correct if some of the T 's a cannot be reached by one or the other of the weight functions, as long as each T is reachable by at least one of the weight functions. For example, if T was unreachable by w, (i.e., w, (1 ) = 0) then a I I a f(^]g)Ai'(Xig) = 0, since there would not be any t. 's. Furthermore,
PAGE 80
74 f(x2^)/w(x2^) = (Np2W2(x2^)/T)f(x2^)/(0 + P2W2(x23)) = (N/T)f(x2g) and everything follows through as before. All of these methods have the same biases that simple Monte Carlo does: they are exact as N^' and their errors are proportional to 1/^. The advantage of these methods is that they can greatly reduce the proportionality constant for l/v. .^i55S3Siu;s52;:uW'ittTrT^i?'i(A5:^"'^^5y^^ rv=^^.ijc^.
PAGE 81
75 III. DETAILED BIAS SELECTION SCHEMES The following will give explicit details on how to choose an x. with biased probability PgCx), once a positive definite g(x) has been decided upon. We will not give recipes for finding the best g(x) for a given problem, but the "physics" of the situation frequently suggests which regions are the "important" ones. First, N x.'s are selected uniformly at random, and a g(x.) is calculated for each of them. Next, a series of partial sums is calculated Sq = S. = S. + g(x.) 1 11 ^^ 1 '^ i=l ^ A random number, R, is then chosen uniformly from the interval R = to R = Sj^ Then a "p" is determined such that s SÂ„ , < R < S p1 p Finally, x is picked and the remaining N 1 x.'s are discarded. The biased probability of choosing x is proportional to g(x ) PB(Xp) = g(Xp)/S^j^, (11) while the uniform probability of choosing x is %{\) l^. (12) ; j55^iaAC^^..'s^3^;s;;jS>& R.S2:
PAGE 82
76 The relative probability of choosing x is found by substituting Eq. (11) and Eq. (12) into Eq. (6) to get w(Xp) = g(Xp)/(i S^^). (13) Note that, .since S., depends on which set of N x.'s is selected, s w(x ) is a multivalued function. But also notice that w(x ) is distributed about "^ P^ in exactly the same fashion, and for exactly the same reasons, as jr S^ is distributed about /g(x)dxleading s s to the following interesting consequences. Consider what happens to 1 "^ f(xO 77 y Â— L_Li as T^ (causing T, > 0) and as n (and therefore N) x^. 1 Â• I /> \ a a 1=1 w(xi) ^ As in the derivation of Eq. (8) continuity of f(x) and g(x) demands that i=iwlf^ gT^i=i ^ j=i ^ But note that w(t)^^9^V Ng^ /g(^)dx Therefore, as T becomes infinite and as N^<, l^ ^ /f(x)dx for all values of N For small values of N the size of a is partly due to the sampling variance of S^ and partly s due to the fact that g does not imitate f perfectly. Note that for N = 1, simple Monte Carlo is exactly recovered, while for N ^<, ^ s importance sampling is obtained. The best value of N depends on how closely g mimics f and on the relative costs of evaluating g and f.
PAGE 83
77 Frequently, it is more convenient to use two (or more) independent weight functions, g, and g^ [4]. An example of this is the calculation of the energy integral for an ^2^ molecule. The potential energy causes the integrand to be singular for two separate instances when the separation of the electron from either of the two nuclei goes to zero. Then one weight function could be designed to follow the singularity due to the electron's close approach to one nucleus, and the second weight function could follow the singularity due to the electron's close approach to the second nucleus. Neither could be suitably used by itself, because the remaining singularity would make a excessively large. Then gj could be used (in the same way g was) with probability p, to select a fraction, p,N, of the points, and gÂ„ could be used to select the remaining p^N points, p^ = lp,Since the probability of selecting gj (g^) is p^ (p^), and the probability of selecting x once g^ (g2) has been selected is g](xQ)/S,^ ^92^S^^^2N ^' ^^^ ^^^'' biased probability of choosing x is PB(Xp) = P(gT)P(Xpg^) + P(g2)P(Xpg2) = Pigi(Xp)/S^. ^P2^2^V^^2H (14) ^ where P(x g^) is the probability that x is chosen given that g, was. The uniform probability of choosing x was jr ; therefore, the p N3 relative probability of choosing x is w(? ) = Pigi(x )/({s^^^ ) + P2g2(^p)/(^S ). (15) "^ s s ^ s s Note that Eq. (15) is just a particular case of Eq. (10).
PAGE 84
78 One problem with biasedly selecting an x by the above methods is that a can be increased by rare occurrences [9]. For example, suppose a large number of the M x.'s happened to be selected uniformly at random in some small, but "important", region, T Then if x happened to be within T Pn(xÂ„) would be anomalously small because a ts p a large number of the g(x.)'s would be the same size as g(x ). This r would, in turn, cause w(x ) to be small; which would make f(x )/w(x ) r P P much larger than the other contributions, f(x.)/w(x.)5 to Eq. (7) from T, where more probable sets of x.'s were selected. The effect a 1 of the anomalous behaviour of f(x)/w(x. ) within an "important" T II a on o can be mediated by the following averaging process. Once an x has been chosen by g(x) and a biased probability of selecting x Pr(x_), has been calculated, then select a new set of N x.'s by b p s 1 ^ keeping x and reselecting the other N 1 x.'s. Then calculate the new Pglx ), call it Pg(x ), assuming that x was also chosen in the second case. The biased probability to use in Eq. (6) will be 'B^y l^^B^^p) ^ PS^^p^^ g(xp) g(x.) + ... + g(x^) + ... + g(x^ ) g^) g(xl) + ... + g(Xp) + ... + g(XfJ ) The uniform probability of choosing x is still 1/N p sthe relative probability of choosing x is and therefore. w(Xp) = P3(Ip)/(l ) ,::z^'<^tii.i^.s^>crt
PAGE 85
79 Great care must be taken to perform the averages precisely as indicated in order to be sure that the biased probability is correct. It is only when ?Â„ is calculated correctly that w, as defined, corrects for the biasedly selected x.'s. For example, although it might reduce the variance of Pg to average the two denominators and then invert them the resulting number would no longer be the biased probability of selecting x It should be noted that this must be done for e'^ery X and as many Pp(x )'s can be averaged as desired. Furthermore, r D jJ f(x) is only calculated one time, for each x This offers an inexpensive way (when f(x) is expensive to evaluate while g(x) is not) to reduce the dispersion introduced by the above methods of calculating ^i^ ;^^Si*^S^iJ.SgOKJ ^ '7^ y.j^i'i?'?yAi5:jlnjÂ„irt*^^^''
PAGE 86
80 IV. EXAMPLES WITH NUMBERS This section provides an example that explicitly demonstrates the subtleties of the biased selection methods, and how the biased probabilities must be calculated and averaged in order to obtain the correct results. Furthermore, it demonstrates that the results are exact, even for finite M s First, suppose the value of a twodimensional integral, /f(x,y)dxdy, is desired. Assume for simplicity that the unit square, over which the integral is being evaluated, can be divided into nine equal area subsquares, T, 's, so that f(x,y) is a constant in each one of them, a Label the subsquares as follows: y 1 3 i. 3 T7 T8 T9 T4 T5 T6 Tl T2 T3 V ^ ^ 1 3 3 and let the integrand be f(x,y) = = 9 (x,y) i T9 (x,y) Â£ T9 (16) Three different Monte Carlo integrations (simple, biased selection without averaging, and biased selection with averaging) will be performed with this function to see how close each comes to the answer, 1. This allows a comparison of the three methods with one another.
PAGE 87
81 First, simple Monte Carlo: each of the points, (x,y)'s, will be selected uniformly at random from the unit square. Therefore, each of the T 's will be equally likely to be sampled, with a probaa "I bility of Q. If a large enough number of points, N, is used to evaluate Eq. (7) then In = M^4^(^^^ + 4^(T2) + ... + M^f(T9)]= 1 and ^j = M^'^Tl) + ... + N^f^Tg)] = 9, where, when x e Tl f(x) is written f(Tl). Using Eq. (3) and the two above equations give = (^)" Second, bias selection without averaging: two points will be selected from which one point will be selected with a biased probability determined by some g(x,y), as in Eq. (11). This point will then be used to evaluate one of the N contributions to Eq. (6). Let the weight function be g(x,y) = .01 for (x,y) e T2, T4, T7, or T8 = .09 for (x,y) e Tl T3, T5, or T6 n for (x,y) e T9. (17) = .01 If the points are selected uniformly at random from the unit square there will be nine equallyprobable cells they can be selected from. This gives 9(9) or 81 possible ways of choosing two points from the nine cells. Therefore, the probability of selecting two points from Vi^S5rfiiAt'E='^Â£HCii=(i#Â£s>^riin:'E;,i(j;\V'>
PAGE 88
82 a particular pair of cells is ^. Only if the pair of points is chosen from one of the seventeen pairs of cells in Table I can a nonzero contribution be made to Eq. (7). The contribution to Eq. (7) from each given value of f(x)/w(x) is N X (Probability of selecting a pair with given w) x (Bias probability of selecting the point from T9from the given pair of pointswith the given w) x (f(x)/w(x)). Using this and Table I yields I, = ^N(3i)(l)(i) N(3^)(.9)(^) N(3^)(.9878)(j73g)] = 1 <(f/w)'>^ = 5.2469. Using Eq. (3) and the tv^o above equations gives The ratio of the error estimate for simple Monte Carlo to the error estimate for biased selection without averaging is approximately .73. For the above integrand and weight function biased selection without averaging gives a 27% smaller error estimate than simple Monte Carlo. Third, biased selection with averaging: after a point has been selected from a pair as in the second method a new pair of points is generated by keeping the previously selected point and selecting the second point uniformly at random from the unit square. The two relative probabilities obtained from the two pairs of points will then be averaged, as in Eq. (15), to obtain w(x, y). Table II, when used in conjunction with Table I, can be used to calculate Ij,, and a
PAGE 89
83 +o s o +J 3 > O) +J o to JD rO JD O SD, Â•rÂ— +> ITS ^Â— O) 5a E c OJ n3 > Â• r00 cn cu OJ Â•f! SÂ•rna rÂ— Â•IÂ— ^^^ JQ ,*> fO ^Q Â— O Â•^^ sQ.0 c >, re sfO n to Â— s. c/) (X3 (U 1Â— u *_ S rO CT> CU Q.1Â— to fd "O E U O) O +> SSO MO O) 41Â— to O) m 3 to 3 I fO I o (J SOJ C SO ClTn3 O OO 4CD S_ r>, C 'iC 4J .,G~i (B ::3 T+j IÂ— oI E I o O SO) Â•r(1) I+jC JD rÂ— t" (O 0) M 2 JD CO E C o Â•> C Q. O CLrfO CD to cn S_ ro >, c "Â•!Â•IÂ— 4J Â•!Â— CTi n3 CQ T+> 1Â— D. MU O 0) r= 3 Â•Icnic IÂ— C tc O) Â•IÂ•!Cl > ^ 4J Â•!fO U E OT O IÂ— > OJ S_ tU Â•!JE Q. to cn+J ^al .' Â— K ^ vÂ— rO CVI o Q. Â— '<+M cr> (o t*+ to rO .Â— > ^ Â— E 1 Â— Â— Â— o f^ cncM SOJ 1 Â— iÂ—to 1 jz o ,Â— ^ +> E to QJ X s_ > Â•rÂ— 4rO C7) Q. o E s_ r. 03 OJ 1 Â— CJ N CO cn un OS rÂ— O CO cn CO CO IÂ— co CM o to to CO CO cn LO CO CO cn rÂ— CM CO CO CO CO CO cn I Â— coLntocncTicncTi cvj^r^cocncricTicn cn CTicTicncTii Â— roLnto cncnCTicriCM^rco
PAGE 90
84 s: (d o. c i<+c o Q. (U B to w (U en o w o J2 o O> M ? a c o o 0) t/) (U "3 T3 UD c LO Â•+3 r~ (0 CTi 00 CO Â• r. i. ^Â— co 0) <3CO CO f + CM Â• s CM "" '" "" A 5 CM 5 a c fC 1 3 ^CO o CO OJ + ^3=*Â• CD 01 CM IÂ— Â• Â• (C QJ CM ^^ rÂ— i3 5 0) 1Â— > (B =3: > r*^ CO cn CTi Â• Â• 1Â— S"" Â— cu OJ II 21 c Â£=.Â£::: O ru oi ^"Â—^ Â„,Â— ^ o > cri ^Â— O) 01 Â•rÂ— o o S_ T c o S e Â•r" TCO CO Â•r O Q1 C\J O C e 4C 1 OJ O 1QJ > O S>, Q. CD 4> c Â•rCn 0) c rÂ— E > 4) 1 Â— \CTl 'SICTl =:rcn Â•rÂ— '^ Â•^ Â•r~ ^ +> en s fO o J3 0) H c O ro o SO) s_ Â•1CM D. (/) MCD 2 00 c a O E Â•1CM crt cn 1 Â— CO Ln isD CM <* r~CO 01 Sc hIÂ— 1Â— 1Â— 1Â— 1Â— 1Â— 1Â— 1Â— co rÂ— o o o QM CiZ^^ V?
PAGE 91
85 and IÂ„ = sCN(8i)(l)4)(^ + ^;^) + N(8l)(l)(l){f) + N(8f)(.9878)(i){TJ^ +,fg55) + NCgfjC .9878) (i){^)] = 1 <^''>N = 5.1144. Using Eq. (3) and the two above equations gives The ratio of the error estimate for biased selection with averaging to biased selection without averaging is approximately ,98. For the above integral and weight function biased selection with averaging gives a 1% smaller estimate than biased selection without averaging. Note that it also approximately doubles the biased selection overhead. How useful this averaging is depends on the relative expense of evaluating f and w and is therefore highly dependent on the particular problem. The more expensive f is to evaluate, relative to w, the smaller the improvement in a needs to be for the averaging to pay for itself. ^T>fc"^*i*^JK'iVSi^"^^=i>f^3?i'lCii" wrt
PAGE 92
85 V. HYDROGEN MOLECULE As a final example, this section contains a calculation of the ground state energy, using the eleventerm JamesandCoolidge wave function [10], for dihydrogen at a nuclear separation of 1,4 Bohr radii (the potential minimum). Note that James and Coolidge did not use Monte Carlo methods. Instead, they analytically evaluated the more than 140 integrals involved. This example incorporates a number of biased selection methods which can also, in principle, be used for other atomic and molecular systems. The following is useful for calculating the required relative probabilities. The probability density function, p(x), of a random variable, x, gives the probability of selecting a given x, x^., between Xq and x^ + dx as P(Xq) = p(xQ)dx. (18) Frequently, the density function p(h(x)), for a function, h(x), of a random variable, x (whose probability function is known to be Pq(x)), is desired. In this case [7] P,(h(x))=pÂ„{x)g^[. (19) For the special case where x is a random variable, R., which is selected uniformly at random from the interval R. = to R^ = 1 the probability density function is pr)(R.) = 1, therefore, the probability density function for some
PAGE 93
87 function, h(R), of this random variable is Pl(h(R)) = lg^!. (20) As an illustration, consider the calculation of the relative probability for randomly selecting spherical coordinates: r = R^M, y = 2R2I, cf) = 27TR3 (21) (RVs = different random numbers between zero and one, M e the maximum r, u E the cosine of the angle between r and the zaxis, and ^ e the angle between the projection of r in the xy plane and the xaxis). Selecting spherical coordinates randomly would have an advantage over selecting cartesian coordinates randomly when, for example, the integrand was basically a function of r (e.g., the r^ part of P,. cancels the 1/r singularity of the potential). Since the R.'s are chosen uniformly at random and the spherical coordinates are functions of R (e.g., r = h(R) = RjM), the probability density functions are, using Eqs. (20) and (21), p^(r) = 1/M, P2(y) = 1/2, p^{^) = 1/27T. (22) From Eqs. (18) and (22) the biased probability of selecting these coordinates is P^Cr, Â„. *) = ^
PAGE 94
88 Finally, the relative probability of randomly selecting spherical coordinates is, using Eqs. (6), (23) and (24), The integration region will be a large box of length Lx (Ly or Lz) in the x (y or z) direction. The box is large enough to contain all of the region where 'i'(x) is essentially nonzero. Keeping the box finite allows the random numbers to be generated easily without introducing singularities. The two nuclei are fixed on the zaxis with the point halfway between them coincident with the box's center. The coordinates of the first electron will be selected in four ways and those of the second in five (four ways being the same as for the first electron). One way, which guarantees that the entire integration region is sampled, is to select the electron's coordinates uniformly at random from the box. A second way, which accounts for the rather large probability of the electrons being there, is to select the electron's coordinates with respect to the midpoint between the nuclei. Third and fourth ways, which account for the nuclearelectron correlation, are to select the electron's coordinates with respect to one or the other of the two nuclei. The fifth way, for the second electron, which accounts for the singularity due to the electronelectron potentials, is to select the second electron's coordinates with respect to the first electron's. An electron's coordinates can be selected uniformly at random from the box by letting its coordinates be x = RiLx, y = R2Ly, z = RoLz. rifSs.'xi;eCa^^a^K!;:^'
PAGE 95
89 Using Eqs. (5) and (20) gives the expected relative probability of selecting the coordinates uniformly at random from the box of Wg(x) = 1. An electron's coordinates can be selected with respect to the box's center by letting r^ = /R^M^, y =2R21, = R32^ (26) (r E the electron's distance from the center). Using Eqs. (20) and (25) and letting hj(R^) = ^i^^ gives the following probability density function p^(r^) = 2/R^/M^. (27) With Eqs. (18), (22) and (27), the biased probability of selecting r^, y, and cj) can be calculated; then, using Eqs. (5) and (24), the relative probability of selecting the electron's coordinates with respect to the center is calculated to be wjr^, y, (j)) = V/(2TTMjr^). (28) An electron's coordinates can be selected with respect to either of the two nuclei by using a weight function, g(y). An r^^ is selected by requiring ^Nl A(R ) E R,/g(y)dy = / g(y)dy. (29) 'O The remaining coordinates, y and (j>, are selected as in Eq. (26). The probability density function for r^^j can be calculated using Eq. (10), if it is first realized that r^^ is a function of A which is, in turn. ^^ i36^*i^tfiVr.
PAGE 96
90 a function of R. Using Eqs. (19), (20) and (29) leads to h^^Nl) P(A) dA^ (30) /g(y)dy Therefore, using Eqs. (18), (26), and (30), 47T/g(y)dy and the relative probability of selecting the electron's coordinates with respect to the first and second nucleus, respectively, can be calculated (using Eqs. (5), (24), and (31)) to be (31) jr'Ni' ^I'^i 4iTrJ,/g(y)dy INIq and w^l2^^N2' *2' ^2^ ~ g(r,2)V ^^r^2^g{y)dy (32) In practice, a good weight function is g(x) = x^^^(x), where ^(x) is, in general, the HartreeFock nuclearelectronic correlation function (in the example ^(x) is the ground state hydrogen wave function). The HartreeFock functions can be found in tables such as reference [11], Note that the x^ cancels the r^^ in Eq. (32) and that ^ follows the system's wave function quite well for large nuclear separations. To greatly lessen the amount of computer time required to biasedly select an r^^, x^>i'^(x) is evaluated at some number of points, say K, and g(x)
PAGE 97
91 is formed by connecting these K values by line segments. To determine an r^^, R^Q^ (Q. e ;g(x)dx = 1^9^ + 9^_^)U^ x._^), and each x. is one of the K used to evaluate x^'i'^(x.)) is calculated, and then a p is located such that ^p rK ^p+l then ^Nl = ^p 9p/" K ^ 2(A)Q,]VA (A = gp+1 gp and Q^ EQp^^ Qp). ^p+1 ^ If r^. is not one of the K x. 's then g(^Mi) in Eq. (32) is 9(^,1 ) = S ^ Nl ^ (gp^^ gp). ^p+1 ^p Since what happens as an electron approaches a nucleus closely is of interest, g(r) in Eq. (32) is altered slightly at the origin it is made greater than zero, usually we take g(0) = g(l ) Â— to make it more probable to place an electron close to the nucleus. The second electron's coordinates can be selected v/ith respect to the first's as in Eq. (26), with M replaced by Mr Therefore, the relative probability of selecting the second electron's coordinates with respect to the first's is WÂ£(r^, ^, y) = V/(2uM^r^). Note that rp, in Wp(x), is exactly what is needed to account for the singularity of the electronelectron potential.
PAGE 98
92 If each of these methods is used with equal probability then w(x^, x^) = l[wg(x,) + w^(x^) + w^^(x^) + Wj^2(^i)^ ^ 5[wg(x2) + w^(x2) + w^,(x2) + w,^2^^2^ ^ w^(x2)]. (33) If desired, one more averaging can be done w'(x], x^) = 2"^^^^!' ^2^ ^ ^^^2' ^r^' The numbers in Table III were calculated using the relative probability in Eq. (33). The 10,000 point estimate required 14 seconds of Amdahl 470 time, and gave a .002 Rydberg accuracy. The calculation of serves as a check that V was calculated properly, while the calculation of serves as a check that the relative probabilities were evaluated correctly. Note that being correct does not guarantee that the relative probabilities are correct. This is because Hy. E^. for all i, and therefore, = Eil^'./vj./l^/vi.) = E regardless of the w.'s used. It is interesting to note, by comparison, that an unbiased Monte Carlo estimate for 10,000 points has an Ngf^ < 2, = 2.6941, = 7.3062, and = 2.4828. In addition, the error estimates are unreliable instead of decreasing as l/ZN" they remained constant from N = 5,000 to N = 10,000. This is because smaller and smaller interparticle separations are sampled as N is made larger and larger. Even though these separations are sampled the least frequently they make the largest contributions. These contributions are increasing enough, as N increases, so that they keep a from dropping. se^EV (=r(T^j=~i^ '
PAGE 99
93 Table III. Results for the H^ molecule using the llterm JamesandCoolidge wave function for a nuclear separation of 1.4 Bohr radii. Number of Monte Carlo points o a 100 ^eff T V H 21.2 2.277 4.40 2.354 .086 .28 .019 26.0 2.344 4.16 2.369 .083 .44 .015 21.2 2.369 4.75 2.350 .099 .39 .033 200 ^eff T V H 43.6 2.333 4.59 2.361 .069 .54 .017 45.0 2.339 3.94 2.3520 .062 .28 .0097 41.5 2.273 4.90 2.336 .070 .28 .020 400 ^eff T V H 81.2 2.306 5.02 2.361 .049 .44 .013 94.7 2.340 4.16 2.3352 .038 .21 .0069 87.2 2.271 4.85 2.355 .046 .21 .015 10,000 ^'eff T V H 2222.0 2.3475 4.682 2.3446 .0089 .055 .0023 Three independent runs are tabulated. N^^ e (^v.H^./w. )^/I(^H*i'./w. )^, E , H E the system's Hamiltonian. Note the anomalous accuracy of H and that the estimated error, a, does not display l/i/FT convergence until N^. ~ 43. Also, note that =  and = 2, in accord with the virial theorem. All expectation values are in Rydbergs.
PAGE 100
94 CONCLUSION Systematic methods, such as trapezoidal, of evaluating multidimensional integrals have the disadvantage that they require an enormous number of points before they can adequately sample the integrand. Biased selection allows one to reduce N in exchange for using his knowledge of the integrand. For example, HartreeFock wave functions can be used as weight functions to find the integrals of more exact wave functions. In addition, biased selection allows one to adequately sample the important regions where the electrons approach the nucleus, or each other, closely. (These regions are important because small errors in the wave function will be amplified by the Â— dependence of the potential. But since these regions occupy a relatively small volume, they will not be likely to be sampled. This will in turn cause the error estimate to be deceptively small.) The price paid for this flexibility is that the points must be distributed in a nonuniform manner and the correct weight function to compensate for this must be calculated. A final important aspect of biased selection is that it allows for the possibility of finding unexpectedly large values of f/w. (One technique for doing this is to uniformly and randomly sample the regions believed to unimportant, for dimensionality < 12. For more dimensions other methods must be devised: e.g. Y^^', where a
PAGE 101
.95 LIST OF REFERENCES 1. J. M. Hammersley and D. C. Handscomb, "Monte Carlo' Methods," Chapman and Hall, London, 1979. 2. N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller and E. Teller, J. Chetn. Phys. 21, 1087 (1953). 3. W. W. Wood and F. R. Parker, J. Chem. Phys. 27, 720 (1957). 4. R. L. Coldwell and R. E. Lowther, Int. J. Quantum Chem. Symp 22, 329 (1978): R. E. Lowther and R. L! Coldwell, Phys. Rev. A 22, 14 (1980). 5. R. V. Hogg and A. T. Craig, "Introduction to Mathematical Statistics." Vol. 1, Macmillan, Mew York, 1979. 6. M. 6. Kendall and A. Steward, "The Advanced Theory of Statistics." Vol. 1, Hafner Pub. Co., New York, 1979. 7. M. R. Spiegel, "Probability and Statistics," (Schaum's Outline), McGrawHill, New York, 1975. 8. R. A. Waller, "Statistics: An Introduction to Numerical Reasoning." Holdenday, San Francisco, 1979. 9. R. L. Coldwell, Int. J. Quantum Chem. Symp. ]3_, 705 (1979). 10. L. Pauling and E. B. Wilson, "Introduction to Quantum Mechanics." p. 351, McGrawHill, New York and London, 1935. 11. E. Clementi, "Tables of Atomic Functions." IBM J. Res. Develop. (Suppl.) 9, 2 (1965).
PAGE 102
APPENDIX B SOME Ci's WHICH FORM SOLUTIONS TO EQUATION (4) OF PART I This appendix is a compilation of some sets of parameters, twentyfour c.'s, which form solutions to equation (4) in Part I. The headings to the sets of c.'s are laid out as follows: 22. & 1.0 G, where a = 22. ag, rÂ„n = 1.0, P.^ = Pg, and the lack of a "K" signifies that k = 0. If "GK" were used in place of "G" then k = g" C > U 1, 1, 1,1). Note that r = j^ The c.'s are from left to right in each row. 96
PAGE 103
22 G 97 5.1^5253 1.71587 1.50369 .1111^1 .218721 .0976309 .00575232 .01^^7098 .0352 .007157tt .111+696 .0135826 .00657276 1.2810U .0U697S7 66202^+ .512976 .0672783 .0251fi32 .0135535 .000883105 .0084121+7 .003720^+1+ .000189845 22 1,U G 3.39767 2.I+S038 1.16991 .87532 .49554 .426513 .258449 .112614 .0324457 .0232371 .0465648 .0667345 .00111758 .0185579 .0293377 .00739149 .00431471 .000247546 .0151876 .0123873 .00679898 1.52693 .188704 .067455 22 2.5 2.44178 2 .0734526 .00995919 .00831653 .00379001 80047 1.07003 .97485 .475162 .203358 .0317995 .00533546 .00699875 .0254261 .0122042 .00133731 .00175191 .00496228 .00687356 .00771743 .455529 .0635454 .0351453 00357605 2 2 h 5. G .86484 2.68932 1.35956 1.51112 .551501 .508982 .108947 .895791 .51926 .239661 1.82782 1.12354 .482346 .13721 .108638 .0467131 .000768374 .00751892 .0126527 .259243 .569531 .293583 .243716 0142958 22 a 9.79604 2.41877 .221294 2, .00294875 .472025 .0557871 .0290706 .0625061 .0655341 .00116609 .666173 .0281946 .0123794 1.45267 .0114822 , 24931 2.5539 .228159 .259165 .00257433 .00965766 .0175874 .0128155 0443525 22 (i 1.4 6.51765 2.47922 .195927 1.08057 1.29258 .691783 .50457 .250655 .103556 .159826 .180081 .244985 .0870172 .0104181 .0511503 .00527028 .00923478 .00193199 .0073204 .00610576 .0420515 5.80021 .203152 .0199858
PAGE 104
22 a 2.5 2.26643 3.71881 1.23873 98 1.25281 .0935228 1. 571^06 9.31001 5.07G62 3.00956 6.itt;8U9 2.80698 2.276023 2.559859 8.17701 1.71^733 152882 .0U17232 .050691^6 .10U6281 .16979t 10U8 9.99352 .315278 .ltl58JiU 22 & 5. 7.286i^2 1.57907 .246736 1.90526 1.45214 1.178 5.51162 5.62043 2.38639 3.68849 2.50488 4.54574 2.99831 .967757 .638429 0408065 .0420574 .0555465 .093742 .0566238 247029 2.25731 1.8555 1.25996 7 a .808634 .0405503 .0603024 1.37737 7. 23978 G 367838 1.10213 .0577388 .09284 .015956 .100718 4.5581 3.34419 1.1032 1.81124 .927221 5.55354 .848395 .071288 .0384779 296987 .0153895 .788383 4.45407 7 & 1.4 G .481469 .270501 .0183554 .147389 .0558767 .0210455 .0233409 .0115855 .0257772 .00495709 .00158293 .0682038 .0192743 .0099876 .0239055 .0119419 .0111319 .0397716 .156139 .133501 .31089 1.0741 .0541095 .0188915 7 & 2.1 G .732829 .488502 .00924068 .0348024 .518622 .0151572 .0553333 .0407695 .08142126 .0941925 .00598963 .00765814 .000380395 .0000107138 .000217986 .0006464 .000573316 .00268044 .00598939 .00419664 .00811611 .416217 .01156929 .00712904 7 Si 4. G .943901 .783469 .516351 .0475801 .206923 .0736034 .0115854 .257859 .0745104 .0505936 .114812 .515653 .144165 .0030547 .0463784 .0265587 .0129225 .0353934 .111192 .103978 .374468 1.84667 .132754 .00390177
PAGE 105
93 7 h 6. G .3Ii0903 .01621+07 .3771+01+ .356173 .0556194 .0t+85f+77 .0822796 .238177 .0986887 .00683235 .0001310J+1 .00797737 .0115565 .0297882 .007382t+7 .000180389 .077^+891+ .00108238 .0258003 .493195 3.33552 9.89799 .281797 .0378118 6. & .7 G .999872 .301+633 .816051 .0282304 .183647 .0117536 .0267538 .0019381 .0950404 .0140295 .00763626 4.07512 3.03542 .676836 .331668 1.23952 1.56327 .619112 5.90178 4.5099 5.56695 .882361 .120174 .0316187 1. G .993274 .336831 .0322635 .135719 .222058 .0278307 .0215007 .0205931 .0492358 .0190133 .00378005 .335957 .175089 .0811609 .16057 .0847402 .092241 .160525 .48705 .33228 .751135 .977751 .221827 .0192453 6. & 1.4 G .403898 .0425998 .00223583 .000507758 .12307 .0178114 .00533254 .00938574 .0375081 .0139143 .0196725 .00148256 .0570519 .000568257 .00562378 .0138714 .0310901 .093716 .243856 .143064 .143064 .647833 .0349958 .0191335 5. a 2.1 G .512238 .153144 .00874024 .0889581 .0744925 .0755689 .0266955 .00830042 .0420707 .0054518 .0127181 .0216736 .00173564 .000060526 .0015734 .0023727 .0012383 .00382614 .00598454 .00453195 .00972202 .864395 .084437 .000799673 6. & 2.8 1.15887 .499014 .803254 .424744 .154854 .01589995 .116603 .114598 .109916 .0370554 .000694444 .054303 .201445 .0990926 .00510171 .0202202 .01375 .0884194 .19433 .131791 .522135 .802187 .131343 .000880853
PAGE 106
6. & 3.5 G TOO .582931+ .0583701 .126138 .0883582 .f6078 .0689593 .123313 .035809 .161211+ .118069 .Oi+33726 .0180514 .03117U5 .00205878 .000530559 .00071^4579 .00180071 .0001813 .010178 .0132202 .0i+3!;i+58 .908775 .0928176 .0182599 6. a k.S G 2.21822 .175551+ 1.28856 .382912 .33901+1 .00152078 .t+19925 .00197797 .0532977 .0^+7046 .019Â£+078 2.t+2218 .0710333 .01+35891+ .12091+8 .00577958 .20023J+ .00930i+79 .0201+1+99 .551289 [+.01+1+88 i+. 69091+ .10051 .0122861 & GK .508922 .10828 .27076 .OOlt+Oi+2 .15582 .02881+5 .01+61+1+9 .03791 .001+9973 .0022803 .031719 .375981+ .950555 .061+52 .5238 .21+81+2 .83971+ .0099339 2.2307 1.5691 1.6572 .67207 .59126 .032787 6. a 1.1+ GK 5.25938 6.871+176 1.53923 .1+53272 1.5301+7 .0201+71+1+ 1.09572 1.11+1+09 .551053 .129076 .590799 .116889 .152821 .0366009 .0329876 .529631+ .138629 1.92551 6.73739 i+. 90517 7.31+396 t+. 1+5156 5.08985 .51301+3 6. 2.1 GK l+. 78711 1.90755 .219757 .0979517 .8051+52 .1+36595 .01+58215 .0608261+ .1+50621 .13636 .0803552 .00886701 .001511+91+ .001+15771+ .00277251 .00651+569 .0081+9572 .02681+62 .0252691 .0137081+ .0903088 .51+3297 .1+93302 .0223178 6. & 2.8 GK 2.73789 .230026 .13772 .31+581+1 .998111+ .1+82971+ .176575 .0851+19 .1+1+9929 .176539 .0375613 .11+9793 .1+1+0619 .0955906 .108536 .01+79309 .0893655 .1+11+912 .68761+1 .1+78835 1.91+372 7.79768 6.97091 .133261
PAGE 107
.101 6. ?i 3.5 GK 3.769111 .260786 .119386 .53721+8 .212816 .131+583 .0551912 .163327 .471025 .0588323 .108191 .6750U3 3.73598 7.12777 6.91+929 381572 .712753 95111+8 .123978 . 03921+37 .ltf709!+ .970255 .7356i+5 0992998 6. & V.5 6.1+6281 .11+8621 3.1+9655 GK .0321666 3.52096 Â• .593509 .551+595 .151+165 .895801+ 11+1926 1.31+993 .0383212 1.751+82 .0163387 .00290803 .286738 .015798 .118612 .51913 .111+581 i+. 02877 .321+98 .11279 00701+361 6. Si 1. E .35911+9 .386965 .0333688 .0595855 1.06056 .000158135 .012125 .0755991 .051+1+618 .01+21097 .000911+361 .215127 .820923 .0165551 .0132192 .03861+58 .0150526 .222921+ I.IOI+II+ .872019 1.72058 .137521 .01251+66 Â•1.20168 5. & 1.'^ E .93706 .00536058 .000633 .00070221+5 .22921+2 .0131+293 .0168112 .0131+81+2 .01+37538 .0200026 .031+9277 .0019821+7 .0655356 .260237E01+ .0133587 .0021+0967 .0201225 .0155215 .17771+6 .170ttl2 .191+97 .71+01+7 .063606 .011201+8 6. & 2.1 E .903011+ .109996 .0120655 .15721+2 .128985 .135232 .11061+9 .00651+293 .0763787 .0121+782 .02281+08 .0368515 .000123393 .960838E05 .001+81+971 .00322753 .0020571+3 .000117265 .00563137 .00588031 .0291+353 1.21+902 .1871+26 .0005351+53 5 '5 1+ 5 E 1.591+01 .166177 .81365 .369392 .31+1835 .00120815 .1+35381+ .00206602 .01+53201 .01+13191+ .0217937 2.26881+ .0727125 .031+9618 .116659 .001+1+597 .115035 .00733502 .0172796 .1+74062 3.55928 5.36281+ .13304 .0124572
PAGE 108
102 5. & 1. EK 1.02382 .3121^99 .00583193 .0825135 .828956 .0175545 .0287117 .0974467 .0403062 .00233446 .000946246 .0266974 .859021 .00579986 .355382 .0170031 .00373953 .157419 .918115 .734089 .505029 .465357 .338817 .00100792 1.4 EK 4.14915 .122743 .0180227 .00130893 .6558013 .201955 .0483587 .0248006 .289425 .0850231 .0588636 .00335191 .443182 .000178247 .056419 .0797879 .161631 .737539 .0721878 .0710252 2.21844 1.01373 .346956 .0453018 5. a .7 G 1.05627 .259298 .741654 .0139225 .330218 .00436738 .0348232 .00674219 .127492 .0226627 .00398958 7.3068 2.30589 1.5895 .339471 .752712 1.17604 .720449 5.53218 5.2129 10.9573 .910543 .233424 .0142728 5. & 1.4 G .351895 .256686 .111507 .095252 .0572138 .0549598 .0235298 .0127587 .00211684 .00811348 .000553647 .035532 .119704 .0289219 .0426376 .0316559 .00729109 .0758738 .191986 .123062 .241364 .538105 .138458 .00475326 2.8 G 1.01928 .509853 .781018 .298329 .272656 .0590349 .0914175 .0753137 .123082 .07433555 .0722998 .311707 .342193 .125315 .00954091 .0943528 .0196233 .155103 .171516 .0793897 .19793 .759801 .0901513 .0101635 5. & 4. .890985 .0419985 .642047 .035573 .484857 .356288 .222302 1.59593 .430792 0191797 .315474 .0385253 9.85889 1.85158 .0119251 1,04209 204482 0100374 .0910031 119823 .677723 .176957 2.33009 .05767
PAGE 109
5. a .7 Gl< 103 .880133 .0083i^526 .0207695 .000597221 .0222559 .000275672 .00209088 .000320598 .0078611^5 .000712109 .855it2E0li .262764 .1I7525 .12U2908 .01^7284 .0128878 .0900391 .0599751 .361025 .318289 .716716 1.22457 .462259 .0191604 5. & 1.4 GK 1.205 .336949 .0121651 .0932783 .0132759 .100572 .019481 .0504815 .00580886 .000614979 .829535E04 .0614424 .658208 0490575 .725265 .0012495 .01815 .496898 1.40315 .909643 1.60067 1.71658 1.7921 .025795 5. 2.1 GK 8.21846 .225419 .563867 .355504 .338896 .0219977 .273293 .00517749 .698519 .376634 .232002 .160929 .000593355 .000453545 .00212904 .000192353 .00172613 .00679883 .0076754 .00240923 .0277538 4.18932 1.42478 .000654507 5. & 4. GK .537897 .00730111 .283427 .0222135 .202403 .00208664 .385397 .0297808 .0135025 .0109745 .174762 .357908 .417361 .0422085 .779853 .00762719 .059207 .00106257 .0123419 .557983 2.85379 1.88094 .0177811 .0447337 4. ci 1,34137 .0295193 .0719105 .0243236 .00036063 .0498882 .00474802 .00101707 .00472455 .000190988 .0108448 2.26266 .0347953 .5594155 .735385 .800162 .588527 .46292 .00518431 1.0059 2.3849 2.01827 .248512 .00210553 & 1. .477832 .0200384 .155021 .0580049 .00251883 .0201761 .000688027 .00518321 .00565811 .200719E04 .0109771 3.43527 1.31028 .0977203 .312343 .00389155 .785948 .045583 .675351 .125076 7.9765 6.14652 1.16283 .0334169
PAGE 110
k. & 1.'4 104 .i+6it054 .0167439 .217I5U .00U57966 .110493 .0280303 .0430283 .0227969 .00195698 3.60013 1.42279 .19831 .14783 .922565 .180047 .808834 7.94649 1. 00001 .0492232 .0079566 .449941 .00775902 4. & 1.75 G 1.58386 .257572 .414904 .658333 .199658 .200211 .123333 .0435489 .0554079 .0293466 .244442 .0857579 1.45333 .34844 .080133 .0158724 .0959436 .0534119 .553248 .712648 1.4507 1.97703 .249875 .00709993 4. & 2.25 G 1.47644 .395068 .546909 .647157 .18311 .0164863 .0724931 .21139 .00561677 .0137878 .436996 .138887 1.51684 .342682 .112284 .0575647 .0545064 .0973971 .506628 .761794 1.38135 1.98354 .194017 .023575 4. & 2.8 G 1.43478 .443678 .574064 .626631 .0990152 .0690079 .0450037 .275788 .000362574 .0216476 .497174 .149369 1.56209 .31789 .105251 .0479535 .0574437 .0833989 .519053 .747503 1.39776 1.99109 .0637323 .00651043 4. & 3.45 G 1.51982 .138012 .622281 .4611906 .115573 .248374 .162853 .197609 .0873137 .0710092 .229392 .986832 .10661 .102881 .128541 .034485 .0357477 .124313 .00225508 .291915 .795946 .847298 .00740728 .00568203 4. & 7 GK 1.2186 .106932 .0668585 .0273917 .106244 .02666651 .0115813 .0533193 .0167068 .0104495 .00147221 4.20263 .052549 .979817 .330732 1.57531 .0255937 .136045 .000998754 1.2519 6.51349 .535701 .488841 .017845
PAGE 111
k. & 1. GK 105 1.11592 .15621 .0870161 .0351861 .10629!; .033635 .0211912 .07lf6079 .0325898 .0087i;lU6 .0020378 5.82896 .0765675 l.hhll .53707 2.3U23 .0303509 .167822 .0011li7i;6 1.70953 9.36861 .509093 .1+81717 .0097323 h. & l.Â£; GK .9191 , .0351*289 .00218033 131*083 .0120291 5. 57211* 0755982 .07138 .05t*i*314 0310222 .037941*6 1.28115 .06t*6937 .00735838 .1*57757 1.16516 .02971*32 .12915 .00113571 1,521*98 9.78121 .729531 .657791* .0338287 1*. & 1.75 GK 1.19517 .1*1*3559 .0330208 .153753 .356383 .107757 .01*18882 .01*53108 .162186 .0752388 .000836007 .511*238 .01*5803 .000550776 .1271 .000160521* .00172179 .01*00537 .1271*88 ,29822 1.31*006 2.13725 2.18828 .0565107 1*. a 2.25 GK .1*76618 .00203571* .0180382 .000817231* 0100676 .0700191 .00381*932 .0551531 .0263839 .0260537 .0832772 .00706662 00232655 .515817 .0137673 .0605093 .501903E01* .0022957 .0227172 .19251*9 1*. 1391*1 3.1*9581* .03681*1*7 1*. & 2.8 GK 1.3005 .00061*607 .231*951* .016261*9 .00198002 .0165335 .197211 .0520056 .21*851 .0171*01*5 .0297961* .0765833 i*. 06888 .0363056 .279387 .0239215 .000986573 .00932823 .136613 .52521*7 2.0173 I*. 95212 2.1*51*51* .0670788 1*. a 3.1*5 GK 1.1*1*876 .119917 .808179 .1*76972 .27387 .351*731* .318395 .31205 .08611*72 .1*621*1*5 .271526 .181311* .0271339 .106382 .321*276 .01*50396 .021731*9 .183537 .00201017 .273571+ 1.211*1*2 .3953 .0133977 .00613955
PAGE 112
106 it. & .7 1. .000t66l51 .00520523 .02203itlt ltll552E0i+ .00271^881; .00031+0835 .000190829 .00011+781+3 379363E0I+ .00308891+ 1. .li+7 .1321+56 2.71289 2.87831 2.5ttig .7G3603 .0205172 3.30906 6.88461 1.191+65 .15801+1 .00117802 t+. & l.k E 1.0l+95i+ .0738691 .0298693 .05i+0t+86 .000225375 .0001+93738 .251657E0I+ .012523 .00715601+ .0001+57317 .0021+6386 .1+981+56 .0316727 .11+2918 .11+203 .287652 .0397597 .0031+0031 .192288 .000889265 i+. 33326 7.73088 .00237376 .06101+1+1+ i+. & 1.75 E 6.02688 .89091+5 .00939682 .1+01+658 .0131+625 .0181+903 .00606713 .0212773 232101+E01+ .031+3392 .031069 .01+01281 .699581+ .01+6861+5 .115585 .0179792 .0021+9136 .1571+1+3 .031+2815 .0003591+1+7 1.05119 3.93391+ .108828 .000356911+ i+. & 2.25 7.85916 .531162 .025901+1+ 1+32213 209371 1+3156 .607039 .395779 .338201+ .021+8839 .166152 1.31903 .97761+1+ .157058 .295362 .01151+01 1.0655U 376396 .191+229 .75117 ,665527 7.15127 .059679 .0183839 i+. & 2.8 E .508639 .0982552 .0222123 .051+91+557 .0001+29258 .00193939 .01+09091+ .07681+09 .112198Eni+ .0011388 .00372063 .0218031+ .15691+ .0031+1+1+55 .0890257 .0173317 .0151559 .0601161+ .113709 .198512 .716933 2.1+181+8 .00771705 .00105531 i+. & 3.1+5 8.71381 .33111+9 , .0201789 .310502 .769017 302317 .6901+98 067951+5 1.05231+ 5.67635 1.27287 .975751 1.23709 53721+2 .130226 .18221+6 .1+31521+ .11+1376 .352753 1.1+0271 .0308805 .860379 0321+755 .009531+38
PAGE 113
107 3. & .3 G 4.5i+957 ~.0t^69112 .00571881^ .22997G .0010253f .0881^06 .0208195 .39252J .1871^ii9 .030/t^59 .006i8935 2.31205 .^126022 .22li352 .707107 .138503 .328661 1.92961 .7583393 3.02612 1^.06172 2.031:^9 .101+738 .000636008 3. & 1. G 7.31958 3.071+11 .91+805 .6581+3 1.51335 .658162 .651+536 .01+61+1+8 .1+50327 .1589536 .11+51+93 3.8731+5 .011+5591 .8711+71+ .1+66132 .518968 .110711+ .167211+ .576905 .0338673 1.15158 i+. 57591+ .291272 .001+79572 3. & 1.2 7.37126 .61+1+31+ .01+71+759 .1691+35 09256 1.00755 .551031 1.1211+1+ 716007 .071+7596 .53826 .155878 l+. 69158 .0255857 .575091 .15QG96 371+616 .297299 .535181+ .0813895 .555101+ I+. 52951 .590035 .0321+559 3. iS 1.1+ 3.18182 .000171+876 .00093291+9 .000902208 .531+191+ .0006251+1+5 .000271+685 .0021+3311 .0007261+96 .00630803 .000512782 1.65 1.26702 .00265873 .01+52621+ .115553 .0110213 .292955 .31+711+9 .939683 1.01366 .59381+1 .059868 .0071+0655 3. a 1.7 G 3.18153 .000600121+ .000697515 .000982331+ .6357 .327607E0I+ .000380021+ .00335713 .00161169 .000275971+ .000620208 1.65198 1.26789 .000601359 .01+1081+9 .11101+9 .00658103 .288831+ .31+31+89 .91+2622 1.01067 .5961+7 .00316225 .000355323 3. a 1.95 G .59359 .1551+29 .218206 .0579716 .0980123 .0661581+ .107275 ,117531 .0401+155 .0032351+5 .0187101+ .3891+07 3.2221+3 .35721+8 2.65393 2.17258 .335101+ .238991 .9011+65 3.61+053 6.0388 7.11+031 .1+85591+ .111+81+1
PAGE 114
108 3. & 2.2 G .2557It .00791^733 .000105811 .0111*902 .031661*3 .0020656 .00183107 .0055899 .000157365 .01*51*615 .0078063 21*. 261*1* .4931*01 .0255167 1.0921*1* .0680123 .0261596 .00529 .565882 .751697 1.051*19 2.33021* .00167775 .00058601*1* 3. & 2.55 G .525261* .237911* .196799 .0753161* .00121*258 .331*862 .00186811 .0223232 .0235775 .173906E01* .01*08157 .50025 2.79318 .0091871*1* 2.991*1*8 2.38021* .278783 .265502 .881*892 3.59059 5.9718 7.15325 .114133 01*781*15 3. & .8 5.377 .1*1*9175 .0122926 .122373 GK 00235305 .0159123 .0688899 .125051* .00325769 .00091721*1 .0291*095 .00177857 2.1861*8 .2581*57 .0205866 .0285267 01*90895 171355 .000132097 .196505 .339371* 1.05071 1.05027 .000392251 3. & 1. GK 1.37915 .00125133 .0128358 .175917 .0229163 .05631*18 .000513313 .058718 .00633293 .00167228 .0106791 2.50361 .0470755 .106002 .200307 .0177135 .03911*81 .0249471 . 00594399 .0392893 .09606 .569114 .744293 .711075E04 3. a 1.2 GK 1.71265 .000816345 .0103548 .156562 .011778 .0734655 .195004E04 .0933257 .00199497 .00179311 .0100391 3.5082 .0718827 .0112093 .0389181 .0135897 .00574354 .0430364 .00941467 .0477697 .113945 .740692 .861249 .785136E04 3. & 1.4 GK 1.87963 .00018834 .0144051 .136536 .0295948 .0531575 .347293E04 .076596 .00135121 .00299131 .00837353 3.82941 .0570866 .0506361 .00737841 .0121464 .0279308 .0332995 .0105837 .0513134 .13864 1.01433 1.090294 .000118932
PAGE 115
109 3. a 1.75 GK ,50kkk6 .00155579 .13191l .0002!i8318 .00120853 .0001591^25 .000579926 .77126ijE05 .00031^9501 .158105E0I .153ii52E0l+ .176164 .00208872 . 0001+11351 .921+127 .00207L>01 .936073E0t+ .317731E06 .204651 .539558 3.03406 3.41088 2.99772 .0350205 3. & 2.2 GK .583401 .00309683 .000540415 .00580322 .500787E04 .421803 .0101896 .00355381 21684 .000378205 .00270008 .817055406 .00112048 .784774E04 00880374 .00203501 1.05312 .3750441E06 .301038 1.04555 4.51092 7.21316 2.48423 .104707 3. a 2.55 GK 2.1785 .00432415 .974541 .000237863 .000471424 .000303637 .0529544 .721503E05 .00197343 .000177692 .286876E05 1.2237 .0448264 .000689399 .160338 .0535872 .0031392 .29754E05 1.25398 5.86435 6.09737 1.12854 6.34452 .202719 3. & .8 E 9.85374 .0903145 .425509 1.35953 .556749 .0017652 .0775017 .0448593 .0152553 .233185 .457196 1.00992 .84681 .043007 .137123 .0233745 .00923187 .052861 .0585851 .0653404 .253065 5.9881 590429 00959872 3. & 1.3 E .248764 .0121822 .00157445 .0115359 .101752 .010358 .664804E04 .000532329 .000513565 .0656477 .00110905 .0525031 .00286657 .0175151 .0259384 .0290269 .000900645 .049882 .0360132 .0309902 .445575 .952933 .0588593 .000558181 3. & 1.5 E 1.65403 .0424279 .062289 .0350777 .0369542 .704132 .00708947 .236553 .429861 .00907606 .0555782 .233732 .388013 .422118 .500828 .0289222 .0179963 .179528 .0584534 .049494 1.5994 4.00586 .130799 .0226223
PAGE 116
a 1.7 110 1.03661^ .019329 .0520827 .0929843 .33155^^ .0180385 .01668it9 . 015f+I+6l^2 .00792921^ .1161^72 .0761892 2.5171+8 .0225852 .0771914 .0318474 .684515 .0068679 .2072 .430328 .741016 1.4782 5.93666 .170697 .0197555 3. & 2.2 .505918 .0797304 .027905 259495 .0183721 .0391596 .00302436 .00244577 .0304951 .0391204 .433413E04 730746 .189141 .0185869 2.48587 .365944 .11239 .125445 1.10173 1.72249 3.89127 8.13134 .222414 .0315895 3. a 2.55 E 1.14117 .927488 .328315 .010948 .184447 1.37797 .00374476 .121085 .0691217 .100349 .215946 .0219576 .00484935 .000372127 3.06324 .0328806 .0480139 .306817 1.74555 1.97162 5.31645 7.57 .148788 .0336732 3. & .8 EK 9.00884 .0104799 .0283241 .263324 .0675639 .001448 .0167284 .167606 .508565 .0296039 .000528373 .525359 .977103 .0703406 .0955005 .150209 .23801 .0564694 .21159 .172021 .572865 7.17687 8.77046 .0355244 3. & 1. EK 3.88425 .872749 .236505 1.15792 1.19953 .1+1+^1+63 .161751 .0189474 1.45377 .234358 .095499 9.45776 1.11258 .253461 .735305 .123893 .0142571 .121973 .889028 .00592028 1.15349 9.20654 1.59249 .00247565 3. & 1.7 EK .438773 .0015238 .00507751 .000509462 .110447 .000636213 .00486389 .00411117 .000606257 .00319358 .000268826 1.90009 3.191919 .02482198 1.47752 .803197 .0286415 .976997 .252868 1.40383 .875955 4.85769 .00509545 .00309971
PAGE 117
3. & 2.2 m EK .767891 .35020ft .221753 .0637921; .00l;85957 .366187 .001+68831; .0236328 .lli;7i;l .00016502^; .01;83025 2.30821 2.58719 .1632i;2 2.7t;302 1.66928 1.556 .121151 1.37077 2.55i;67 2.49i;38 7.1921*4 .589523 .lt;3086 2.5 & .7 G 3.12731; .3i;3355 .0108685 1.039t;l .0130322 .00687696 .00509265 .0459078 .158252 .0051975 .309554 5.63389 .6056576 .0308444 3.30814 .0198433 .0749837 .400315 3.1478 .170127 5.21059 3.11545 .199584 .00347407 2.5 a G 2.31433 .354808 .0108066 1.1705 .020443 .00947228 .00933654 .0502323 .192271 .0172044 .510483 6.45437 2.42413 .0420169 1.76911 .0276455 .0768791 .46458 1.94692 .1656 2.71905 2.28952 .13177 .0034102 2.5 a 1,. 2 1.92097 .106647 .217264 .21453 3.5945 G .319094 .00712546 .148525 .0346882 6.211991 2.14314 .0314638 .370219 2, 638403 .0856298 .232141 .00210723 .0539591 3.70569 2.3795 .0876014 14221 .125532 00362669 2.5 & 1.4 G 1.5743 .219832 .048025 .550603 .0841845 .189597 .0958006 .0204073 .289116 .0204817 .320556 7.96244 3.25056 .139658 2.53559 .140543 .0023666 .188339 1.30422 .0807703 1.3336 1.7692 .0709753 .000227715 2.5 & 1.75 G 1.90362 .126961 .246556 1.80271 .00187333 .467704 .377433 .0623958 .34587 .0372722 1.20556 1.73606 .89766 .80353 .139785 .39685 .0243075 .419039 1.06141 .117238 .745717 .939662 .000413687 .0031294
PAGE 118
112 2.5 a 2, ,880395 .371339 .569392 .503315 .1103U1 .37862it . 021+3136 .00591095 .193692 .0133993 1.151+1!+ .515512 .713726 .237087 .600598 .985535 .0258201 .00197f+65 .3615f+6 .2161+52 .376331+ 1. 0. 0. 2.5 0! 2,15 732913 291+305 561+126 t+87061 .lf+595 .17531+5 .585525 .0288397 .263993 .0091+9506 .0151023 l.t+Ui+55 .181+851 .00169938 535112 .1+51972 0050575 518621+ .857519 17321+8 .291257 .883707 .00100188 .00091+575 2.5 a .7 GK 1.69269 .21+7301 .307191 .323369 .00821718 .00561+997 .0283391+ .00550127 .0023061+3 .773192E01+ .0777091 5.580711 .139181 .18251+1 .87069 .1+35138 .00052781+2 .1+25563 1.1+2891 .0555672 i+. 30958 1.57651 1.771+02 .000300383 2.5 & 1. GK 1.85835 .282381 .509886 .1+35532 .0118781 .0112552 .0522782 .001+77231+ .00091+8503 .000375387 .11+121+6 l+. 85359 .1921+72 .1+6251+7 .831582 .1+81+277 .00305192 1.33921 1.93031+ .0329698 3.1+5821+ .590357 .636391+ .000111688 2.5 & 1.25 GK 1.89625 .01+33578 .0181619 .01+9779 .008261+05 .0001+51+365 .075751+8 .00151659 .26323 .0091+5891+ .01+50531+ .335803 2.931+95 .631096 .618001 .0585302 .00175778 .892575 .981799 .0363507 .77771+5 .67315 .537619 .01051+08 2.5 & 1.1+ 2.27823 .228791+ .552801 GK .1251+13 .0188551 6.11+832 .1361+63 1.28001 .11676 .0195081 .231865 .001+1+0673 .305329 .071227 .051+9223 .2591+33 .00925521 .0929985 .325338 .1+97591+ 2.362875 1.5979 1.5S7I+5 .0222885
PAGE 119
2,5 & 1.75 113 GK 2.75524 .00290017 .00115782 .098551 .000738376 .000179532 .0073079 .00229562 .530769 .00113228 .0115516 .8061+3 .677491 .0224857 3.114214 .0047878 .000110941 .398021 2.07058 .0389272 4.27656 2.19704 1.39469 .030199 2.5 & 215 GK .370562 .00222894 .000673023 .0128396 .53239E04 .000117607 .000483835 .00097995 .129331 .000156927 .00326745 .160754 .0294859 .0228599 .505545 .000315916 .179892E05 .0612475 .818405 .0138183 1.19826 .308589 5.54238 .00452586 2.5 a 1. 5.1194 .0385783 .259048 .375664 .247376 3.48253 .0700439 .328775 .0290115 .169508 .033565 .017309 6.79051 1.40541 1.17974 .335617 .165585 2.5251 1.6515 261532 .942128 1.17803 .0583131 .00326765 2.5 a 1.1 5.12517 .213843 .0311153 .319889 .0139709 .301488 .371758 .198875 .0596792 .0325275 .209978 4.34689 5.15337 .843351 1.12143 .0580962 .0970094 2.57285 1.5605 .00231555 .826448 1.17927 .017953 .00372299 2.5 & 1.2 E 3.07485 .332756 .204362 .551465 .743215 .124694 .109272 3.85797 10.3012 175381 .327187 .320149 .0729359 .497307 5.09774 .0434727 .036299 2.78676 2.59686 .125148 .953589 .850428 .0232027 .00247386 2.5 a 1.4 1.24877 .168445 .141241 .332148 .0188992 .523531 .0914282 .0433416 .450358 .0390543 .0573283 1.69359 1.32811 .120077 4.09392 1.07539 .0580292 2.34624 3.26934 .447635 2.00775 2.20387 .0789998 .0149201
PAGE 120
2.5 1.75 114 9.30867 .597198 .187005 1.87098 .01^22725 .108313 .28355 .236104 .1+61+553 .20278 .0805101 2.1U811 .104295 1.01489 2.24482 .261669 .125258 7.5701 5.24055 .150421 1.20365 1.5141 .0690978 .00250879 2.5 & 2.15 3.92124 .968814 .0760087 .796784 .0287254 .172552 .426077 .232567 .376058 .0506417 .121489 1.94693 1.25445 .411459 1.72555 .243261 .611487 4.01324 3.52209 .509019 .75868 4.18482 .059053 .0587498
PAGE 121
ALPHABETIZED BIBLIOGRAPHY Andreoni, W.,.Phys. Rev. B. U, 4247 (1975). Ashcroft, N. W. Phys. Rev. Lett. 2]_, 1748 (1968). Bellemans, A., and M. De Leener, Phys. Rev. Lett. 6_, 603 (1961). Calais, J.L., Arkiv Fysik 29, 255 (1965). Carr, W. J., Phys. Rev. J28, 120 (1962). Clementi, E., "Tables of Atomic Functions." IBM J. Res. Develop. (Suppl.) 9, 2 (1965). Coldwell, R. L., Int. J. Quantum Chem. Symp. 13_, 705 (1979) Coldwell, R. L., and R. E. Lowther, Int. J. Quantum Chem, Symp. Jl, 329 (1978): R. E. Lowther and R. L. Coldwell, Phys. Rev. A 22, 14 (1980). Coldwell, R. L., and M. A. Pokrant, J. Computational Phys. 20, 365 (1976). Â— Coldwell, R. L., J. Computational Phys. 24, 223 (1974). Etters, R. D., R. Danilowicz, and W. England, Phys. Rev. A 12, 2199 (1975). Â— Hammersley, J. M., and D. C. Handscomb, "Monte Carlo Methods." Chapman and Hall, London, 1979. Harris, F. E., L. Kumar, and H. J. Monkhorst, Phys. Rev. B 7, 2850 (1973). Harrison, W. A., "Solid State Theory." Dover, 1979. Haselgrove, C. B., Math of Computations, 15_, 323 (1961). Hawke, P. S., T. J. Burgess, D. E. Duerre, J. G. Huebel R. N. Keeler, H. Klapper, and W. C. Wallace, Phys. Rev. Lett. 41_, 994 (1978). Hogg, R. V., and A. T. Craig, "Introduction to Mathematical Statistics," Vol. 1, Macmillan, Mew York, 1979. James, H. M. and A. S. Coolidge, J. Chem. Phys. 1, 825 (1933). Kendall, M. G., and A. Steward, "The Advanced Theory of Statistics." Vol. 1, Hafner Pub. Co., New York, 1969. 115
PAGE 122
116 Kolos, W., and L. Wolniewicz, J. Chem. Phys. 43_, 2429 (1965). Kolwalik, J., and M. R. Osborne, "Methods for Unconstrained Optimization Problems,." American Elsevier, New York, 1968. McMahan, A. K. "High Pressure and Low Temperature Physics." p. 2141, Plenum, New York, 1977. Metropolis, N., A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, J. Chem. Phys. 21, 1087 (1953). Mitton, S., "The Cambridge Encyclopaedia of Astronomy." Crown, New York, 1977. Moiseiwitsch, B. L., "Variational Principles." pp. 172183, Vol. XX, Interscience Pub., 1966. Monkhorst, H. J., and J. Oddershede, Phys. Rev. Lett. 30, 797 (1973) Myron, H. W., M. Y. Boon, and F. M. Mueller, Phys. Rev. B 118 3810 (1978). Neece, G. A., F. J. Rogers, and W. 6. Hoover, J. Comput. Phys. 7_, 621 (1971). Niederreiter, H., "Diophantine Approximation and Its Applications." Academic Press, New York and London, 1973. Oddershede, J., L. Kumar, and H. J. Monkhorst, Int. J. Quantum Chem. Symp. No. 8, 447 (1974). Ostgaard, E., Z. Physik 252, 95 (1972). Pauling, L., and E. B. Wilson, "Introduction to Quantum Mechanics." p. 351, McGrawHill, New York and London, 1935. Ramaker, D. E., L. Kumar, and F. E. Harris, Phys. Rev. Lett. 34_s 812 (1975). Ross, M., and A. K. McMahan, Phys. Rev. B U, 5154 (1976). Ruoff, A. L. "High Pressure and Low Temoerature Physics." pp. 214K Plenum, New York, 1977. Schneider, T., Helv. Phys. Acta 42, 956 (1969). Spiegel, M. R. "Probability and Statistics." (Schaum's Outline), McGrawHill, New York, 1975. Waller, R. A., "Statistics: An Introduction to Numerical Reasoning.' Holdenday, San Francisco, 1979. Whitmore, M. D., J. P. Carbotte, and R. C. Shukla, Can. J. Phys. 57 1185 (1979).
PAGE 123
Wigner, E., and H. B. Huntington, J. Chem Phys. 3_, 755 (1935). Wood, W. W., and F. R. Parker, J. Chem. Phys. 27_, 720 (1957). Ziman, J. M., "Principles of the Theory of Solids." Cambridge Univ. Press, 1972. 117
PAGE 124
BIOGRAPHICAL SKETCH Stephen A. Burdick was born in San Antonio, Texas, on July 27, 1952. He attended and completed his undergraduate work at the Georgia Institute of Technology on a tennis scholarship. After receiving his B.S. degree from Georgia Tech, he continued his studies at the University of Florida, from which he received his Ph.D. in December 1980. He and his wife, Linda, have one daughter, Leslie. 118
PAGE 125
I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. / ^^ Arthur A. Broyles, Chairman/ Professor of Physics I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. F. Eugene Dunnam Professor of Physics I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. Kwan Chen Professor of Astronomy
PAGE 126
This dissertation was submitted to the Graduate Faculty of the Department of Physics in the College of Liberal Arts and Sciences and to the Graduate Council, and was accepted as partial fulfillment of the requirements for the degree of Doctor of Philosophy. December 1980 Dean for Graduate Studies and Research
xml version 1.0 encoding UTF8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchemainstance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd
INGEST IEID EBYSOOV5K_ROAE5X INGEST_TIME 20141208T22:14:56Z PACKAGE AA00026482_00001
AGREEMENT_INFO ACCOUNT UF PROJECT UFDC
FILES

