Citation
Metallic hydrogen--a two-body approach to band theory and a quantum mechanical test of diophantine methods

Material Information

Title:
Metallic hydrogen--a two-body 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:
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 )
non-fiction ( marcgt )

Notes

Thesis:
Thesis (Ph. D.)--University of Florida, 1980.
Bibliography:
Includes bibliographical references (leaves 115-117).
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 non-profit 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 HYDROGEN--A TWO-BODY 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 TWO-BODY 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 Two-Electron 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 HYDROGEN--A TWO-BODY 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 2-electron wave func-

tions. Born-Oppenheimer 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 insulator-conductor
8B rs3
transition was found to occur between rs = 1.7aB and r = 1.9aB,

while the solid was still molecular. A molecular-atomic transition

was predicted, by extrapolation, to occur for rs= 1.1. Finally, a

significant broadening of the k-averaged Born-Oppenheimer 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 insulator-conductor transitions.

The second part of this work details the application of Diophan-

tine methods to the evaluation of multi-dimensional quantum-mechanical
1
expectation values. It shows that a convergence was found for a

substantial range of N, the number of multi-dimensional integration

points. For the same accuracy as biased selection, Diophantine

methods required one-fourth the number of points, for the remaining

range of N.































vi














INTRODUCTION

This work is most logically divided into two parts. Part I,

Metallic Hydrogen--A Two-Body Approach to Band Theory, details the

methods and results of a solid-state 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 multi-dimensional integrals that is faster than the Monte

Carlo convergence of 1/vT. Parts I and II are self-contained: 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 TWO-BODY 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 2-electron 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 k-averaged Born-Oppenheimer 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 molecular-atomic 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 atomic-hydrogen 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

insulator-conductor 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 atomic-hydrogen 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 H2-H2 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 molecular-atomic

transition density. This work casts doubts on this assumption near the

molecular-atomic transition. The H-H potentials were found to depend

on rs, which makes it highly likely that the H2-H2 potentials will

also. Furthermore, it is not clear how an H2-H2 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 Hartree-Fock [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

molecular-atomic transition pressure uncertain by a factor of about

two or three. This work includes the properties of a body-centered

cubic atomic-hydrogen 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 electron-electron distances explicitly.





6


The molecular-atomic transition pressure is usually calculated from

the slope of the common tangent to the atomic-solid energy-versus-

density and the molecular-solid energy-versus-density curves. This

pressure is also frequently assumed to be the insulator-conductor tran-

sition pressure. This work distinguishes the molecular-atomic tran-

sition, which at 00K did not occur for any density considered, from

the insulator-conductor transition. The molecular-atomic transition

takes place (at 00K) when there is no longer a minimum in the Born-

Oppenheimer potential for any nuclear-nuclear distance other than the

one corresponding to a body-centered cubic atomic lattice. Also, the

energy at the minimum of the k-averaged Born-Oppenheimer potential and

the energy at an internuclear separation corresponding the the body-

centered cubic lattice will agree at the rs where the molecular-atomic

transition occurs. The data of this work were extrapolated to find

this transition. The insulator-conductor 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 insulator-conductor transition in hydrogen [20].

However, a more recent experiment claims to have found one [21]. Also,

there are ongoing experiments to produce an insulator-conductor 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 N-particle 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 2-particle Hamiltonians for

the different cells, the N-particle wave function will be an anti-

symmetrized product of 2-electron wave functions for the cells.

In the Born-Oppenheimer 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 2-dimensional 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 B-l B B
A 2 A-2 A 2 A 2 2-A 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 2-dimensional 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 Ewald-type 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) = e-ik- 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 2-electron wave function instills a certain

amount of self-consistency in the wave function, without iterating.

Hartree-Fock type methods try to solve a 1-electron 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 self-consistency, 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

[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 Born-Oppenheimer approximation for an

"H2-molecule" 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-

ular-axis" 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 -
body-centered 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 k-dependence of i(r)

makes it possible to construct a N-electron wave function from a

product of V's with different k's, as will be further discussed later.

A 2-electron 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 6-dimensional 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 6-dimensional 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

Seitz--an ion, with a spherically symmetric potential, and a single

electron per unit cell--the above reduces to the Wigner-Seitz boundary

conditions:


uk(r) = constant for r > r

Sk = 0.
ar r=rs


Another important similarity exists between the Wigner-Seitz method

and this work. There is an implicit electron-electron correlation

that is assumed for the N-electron wave function. In the Wigner-Seitz

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 cell-by-cell 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

Wigner-Seitz 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 twenty-four 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 2-body form of uk(r) allows the effects of

2-body correlations to be calculated. The twenty-four 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
6-dimensional 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 twenty-four 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 parameters--in this case the grid has twenty-five

points in a 24-dimensional space--for each of which a is evaluated.

If vectors, from the origin to the points, are used to locate the

grid points, then twenty-four 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

twenty-five 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 Two-Electron 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 2-electron states.

For 1-electron 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
... (N-1)>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) ... (N-1)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 (N-1, N); and assuming that E~ < e 2 E t < E 3 etc., and
kN-I 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) ... (N-1 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 6-dimensional, here it consists of two 3-dimensional

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 insulator-conductor 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(2-x)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 s-orbital type of symmetry about the box's

center while PE has a p-orbital 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 s-orbital 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 p-orbital 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

many-body wave function from some 2-body wave functions, if Hartree-

Fock theory is correct. The difficulties can be most easily visualized

by trying to construct a 4-electron wave function from 2-electron wave

functions. The Hartree-Fock 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 1-electron states). Notice that if a
4-body state is constructed from $0 and 0I that three particles will

occupy the same 1-electron state and antisymmetrization will make this

4-body wave function identically zero. Given that 0 is occupied, the

next available state that a 4-body wave function can be constructed
from is 02 = :(1)1((2). A look at the Hartree-Fock 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 2-dimensional 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 2-electron wave functions. Therefore, if Hartree-Fock is

correct, the non-zero 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 2-electron states to form many-electron

states does not generate a "Pauli principle" as it does for 1-electron

states.

Now that the energy bands have been found the k-averaged potentials
2L3 >
can be calculated. This is accomplished by placing TSd electrons
in a volume, di, in k-space. The parabolic form of the bands means that

a k-sphere should be filled until it makes contact with the sides of the
first Brillouin zone. Then only that part of the k-sphere 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 Born-Oppenheimer 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 molecular-hydrogen solid is definitely a conductor at
this density. Note that even at this density the molecular solid is
favored over the atomic-hydrogen 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 Born-Oppenheimer 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 molecular-hydrogen 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 Born-Oppenheimer 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 molecular-hydrogen 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 Born-Oppenheimer 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 molecular-hydrogen 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 Born-Oppenheimer
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
4-1









1 2 3









Figure 9. Comparison of this work's Born-Oppenheimer 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 atomic-solid potential
1 2 3

r5 (aB)

Figure 9. Comparison of this work's Born-Oppenheimer 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 atomic-solid potential
for a body-centered cubic lattice. eM 7 this work's molecular-solid
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 insulator-conductor 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 k-space equal to the volume of a Brillouin zone has been

filled. When the procedure is completed, all N-electrons 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 Born-Oppenheimer

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 insulator-conductor

transition begins to appear. A significant broadening of the Born-

Oppenheimer potential takes place, the k-averaged 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 insulator-conductor transitions.





33







Conclusion

In this work it was found that the molecular-atomic and the

insulator-conductor transitions did not occur at the same density.

The data indicate that the insulator-conductor 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 molecular-atomic

transition occurs at rs = l.laB. Finally, the distinct widening of

the minimum of the k-averaged Born-Oppenheimer potential at rs = 1.97

is bound to have thermodynamical consequences beyond those normally

associated with insulator-conductor transitions.

In this work it was also successfully demonstrated that 2-elec-

tron wave functions could be used in band theory. This work also

opened the door to further possibilities. The Born-Oppenheimer

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. 21-41, 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. 21-41, 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

twenty-four 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 twenty-four 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 11-term James-

and-Coolidge 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
K-dimensions, 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[-i-2 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 14-digit 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
5-dimensional 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 nuclear-electron

separation--was estimated. This expectation value can be determined

by evaluating
c0 s u
L E27 2dsfdufdt(8su s2 + t2)e-2zs
0 0 0
Ss u
M = 2idsfdufdt e-Zu(s2 t2)[ + -] + 4usT}, (2)
0 00 0ss u s

and

CO u
Norm 2r2dsfdufdt u(s2 t)e-2zs
0 0 0

where T a [6]. Haselgrove's method, since it is designed to sample

the 3-dimensional 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)2z-1 2[1 u2t2+ 4 1 s)} (3)
0 zzn(l-s)]

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 difficult--the 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 multi-dimensional

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 11-parameter 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 6-dimensional 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 z-axis

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 6-dimensional xi, W(xi),

exactly corrects for the non-uniform 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 y-axis, E= the angle between the projection
of rc in the x-y plane and the x-axis, 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' = 2R5-l, 1 = 2TR6

w5(x2) = V/(2rMRE)

(rE 5 the electron-electron distance, ME e a constant). The function
g(r) was created by evaluating r2e-2r 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 8-dimensional 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 11-term James-and-Coolidge 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 nuclear-electron 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 particle-particle 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 6-dimensional integration

volume into 729 equal-volume hyper-cubes--i.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 hyper-cube 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 electron-electron 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 electron-electron 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 nuclear-electron 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 nuclear-electron correlations.





54








&
.4 -












/ A&
.2










U &



-3a -20 -a 0 a 2c 3a

Figure 3. N points were placed in a 6-dimensional 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 trial-and-error 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. In-the case of esti-

mating E0, the 5-decimal 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.172-183,
Vol. XX, Interscience Pub., 1966.

7. H. Niederreiter, "Diophantine Approximation and Its Applications."
pp. 132-194, 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, one-half to

one-third 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

2-electron 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 multi-dimensional

quantum-mechanical 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 non-uniformly 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 M-x 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-x.




1 Unless specifically indicated, all sums are from 1 to N.





63






I. ERROR ESTIMATION
The simple Monte Carlo estimate for a multi-dimensional 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- (xi-7)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) = loo-x-o0, 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 equal-volume hyper-cubes (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
multi-valued w(ix)'s used herein is discussed in Section III. The
case for single-valued 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)
T-N 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 N-o. 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 non-negative 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 multi-peaked 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 N-o
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


Sp-1 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 multi-valued 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)dx--leading
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 N-x, 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 = 1-P1. 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 Ns-1 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
4-1
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 two-dimensional 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 equal-area

sub-squares, Ta's, so that f(x,y) is a constant in each one of them.

Label the sub-squares 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) + N-f(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

non-zero 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 T9--from the given
pair of points--with 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 eleven-term James-and-Coolidge 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)) = 1Idh-RT) (20)

As an illustration, consider the calculation of the relative
probability for randomly selecting spherical coordinates:

r = R1M, p = 2R2-1, ( = 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 z-axis, and E
the angle between the projection of r in the x-y plane and the x-axis).
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 non-zero. Keeping the box

finite allows the random numbers to be generated easily without

introducing singularities. The two nuclei are fixed on the z-axis

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 nuclear-electron

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 electron-electron

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 =2R2-1, 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 Hartree-Fock nuclear-electronic correlation
function (in the example'(x) is the ground state hydrogen wave function).
The Hartree-Fock 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 + gj-l)(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 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 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 electron-electron 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 11-term James-and-
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, Hartree-Fock 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

non-uniform 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 TWO-BODY 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 TWO-BODY 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 Two-Electron 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 TWO-BODY 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 2-electron wave functions. Born-Oppenheimer 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 insulator-conductor transition was found to occur between r = 1.7an and r = 1.9an, s B s B' while the solid was still molecular. A molecular-atomic transition was predicted, by extrapolation, to occur for r = 1.1. Finally, a significant broadening of the k-averaged Born-Oppenheimer 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 insulator-conductor transitions. The second part of this work details the application of Diophantine methods to the evaluation of multi -dimensional quantum-mechanical expectation values. It shows that a tjconvergence was found for a substantial range of N, the number of multi-dimensional integration points. For the same accuracy as biased selection, Diophantine methods required one-fourth 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 Hydrogen--A Two-Body Approach to Band Theory, details the methods and results of a solid-state 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 multi-dimensional integrals that is faster than the Monte Carlo convergence of 1//R". Parts I and II are self-contained: 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 TWO-BODY 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 2-electron 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 l-7an and 1.9an (r for liquid hydrogen is 3.4an, -5-r^ ~ N"' '^ ^'^ ^^^ number of electrons contained in the crystal volume.

PAGE 9

L^). The k-averaged Born-Oppenheimer 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 molecular-atomic 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 atomic-hydrogen 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 insulator-conductor 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 atomic-hydrogen 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 molecular-atomic transition density. This work casts doubts on this assumption near the molecular-atomic transition. The H-H 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 H2-H2 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 Hartree-Fock [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 molecular-atomic transition pressure uncertain by a factor of about two or three. This work includes the properties of a body-centered cubic atomic-hydrogen 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 lectron-electron distances explicitly.

PAGE 12

The molecular-atomic transition pressure is usually calculated from the slope of the common tangent to the atomic-solid energy-versusdensity and the molecular-solid energy-versus-density curves. This pressure is also frequently assumed to be the insulator-conductor transition pressure. This work distinguishes the molecular-atomic transition, which at 0K did not occur for any density considered, from the insulator-conductor transition. The molecular-atomic transition takes place (at 0K) when there is no longer a minimum in the BornOppenheimer potential for any nuclear-nuclear distance other than the one corresponding to a body-centered cubic atomic lattice. Also, the energy at the minimum of the k-averaged Born-Oppenheimer potential and the energy at an internuclear separation corresponding the the bodycentered cubic lattice will agree at the r where the molecular-atomic transition occurs. The data of this work were extrapolated to find this transition. The insulator-conductor 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 insulator-conductor transition in hydrogen [20]. However, a more recent experiment claims to have found one [21]. Also, there are ongoing experiments to produce an insulator-conductor 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 N-particle 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 2-particle Hamiltonians for the different cells, the N-particle wave function will be an antisymmetrized product of 2-electron wave functions for the cells. In the Born-Oppenheimer 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 2-dimensional 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 2-dimensional 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 r-jg, 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 Ewald-type 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 2-electron wave function instills a certain amount of self-consistency in the wave function, without iterating. Hartree-Fock 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 self-consistency, 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 Born-Oppenheimer approximation for an "Hp-molecule" 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 "molecular-axis" 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 body-centered 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 k-dependence of $(r) makes it possible to construct a N-electron wave function from a product of $'s with different k's, as will be further discussed later. A 2-electron 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 6-dimensional 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 5-dimensional 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 Seitz--an ion, vn'th a spherically symmetric potential, and a single electron per unit cell--the above reduces to the Wigner-Seitz boundary conditions: u. (r) = constant for r ^r au,, S 9r = 0. Another important similarity exists between the Wigner-Seitz method and this work. There is an implicit electron-electron correlation that is assumed for the N-electron wave function. In the Wigner-Seitz 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 -by-cell 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 Wigner-Seitz 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 twenty-four 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(e-|rii + ^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 2-body form of ur(r) allows the effects of 2-body correlations to be calculated. The twenty-four 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 6-dimensional 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 twenty-four 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 parameters--in this case the grid has twenty-five points in a 24-dimensional space--for each of which a is evaluated. If vectors, from the origin to the points, are used to locate the grid points, then twenty-four 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 twenty-five 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 Two-Electron 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 2-electron states. For 1-electron wave functions and conventional band theory, the system's wave function is found by anti symmetrizing H' = -^ (4) ... *ir (N-l)(|)^ (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 (N-1, N); and assuming that e^ -^ ^ ^t t ^ ^t t ^ ^t t etc., and h-rn 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 6-dimensional here it consists of two 3-dimensional 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 insulator-conductor 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, ) iL-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(7-x)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 s-orbital type of symmetry about the box's center while P^ has a p-orbital 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 p-orbital 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 many-body wave function from some 2-body wave functions, if HartreeFock theory is correct. The difficulties can be most easily visualized by trying to construct a 4-electron wave function from 2-electron wave functions. The Hartree-Fock 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 4-body state is constructed from and *, that three particles will occupy the same 1 -electron state and antisymmetrization will make this 4-body wave function identically zero. Given that $q is occupied, the next available state that a 4-body wave function can be constructed from is ^^ = (j^-, (1 )(])-, (2) A look at the Hartree-Fock energies for ^sE-saijicK-

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 2-dimensional 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 e-p(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 3-8. 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 2-electron wave functions. Therefore, if Hartree-Fock is correct, the non-zero 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 2-electron states to form many-electron states does not generate a "Pauli principle" as it does for 1-electron states. Now that the energy bands have been found the k-averaged potentials can be calculated. This is accomplished by placing y^— ysdk electrons in a volume, dk, in k-space. The parabolic form of the bands means that a k-sphere should be filled until it makes contact with the sides of the first Brillouin zone. Then only that part of the k-sphere 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 Born-Oppenheimer 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 molecular-hydrogen solid is definitely a conductor at this density. Mote that even at this density the molecular solid is favored over the atomic-hydrogen solid.

PAGE 30

24 Figure 4. The tops, £2^0) and £](!<(,), and bottoms, £2^'^s^ ^"^ £-1(0) > of the first two bands, and the Born-Oppenheimer 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 e-j(i<(,) for an r^g that is less than the r^g at the minimum of |^, the molecular-hydrogen 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 Born-Oppenheimer potential, , versus r, energies calculated by Monte Carlo and their o's. not dip below s-rCkp) until r^o is larger than the r.o for the minimum of |^, the molecular-hydrogen 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 e-j (0) Figure 7. The tops, Eo^^) and e-](k ), and bottoms, ^n^^s of the first two bands and the Born-Oppenheimer 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 e-il^p) until an rnp well beyond the r.g for the minimum of , the molecular-hydrogen 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 Born-Oppenheimer 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 Born-Oppenheimer 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 atomic-solid potential for a body-centered cubic lattice, e.. = this work's molecular-solid 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 insulator-conductor 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 k-space equal to the volume of a Brillouin zone has been filled. When the procedure is completed, all N-electrons 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 Born-Oppenheimer 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 insulator-conductor transition begins to appear. A significant broadening of the BornOppenheimer potential takes place, the k-averaged 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 insulator-conductor transitions.

PAGE 39

33 Conclusion In this work it was found that the molecular-atomic and the insulator-conductor transitions did not occur at the same density. The data indicate that the insulator-conductor transition most likely occurs between r^ = 1.7ag and r = 1-9ag, 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 molecular-atomic transition occurs at r^ = l.lag. Finally, the distinct widening of the minimum of the k-averaged Born-Oppenheimer potential at r = 1.97 is bound to have thermodynamical consequences beyond those normally associated with insulator-conductor transitions. In this work it was also successfully demonstrated that 2-electron wave functions could be used in band theory. This work also opened the door to further possibilities. The Born-Oppenheimer 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. 21-41, 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. 21-41, 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 twenty-four 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 twenty-four 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 i-s, 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 11-term Jamesand-Coolidge 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 K-dimensions, 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 14-digit 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 5-dimensional 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. = nuclear-electron 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 3-dimensional 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 multi-dimensional 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 11-parameter 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 6-dimensional 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 z-axis 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 5-dimensional x. w(x.), exactly corrects for the non-uniform 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(x-j., 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 y-axis, (j) = the angle between the projection of r in the x-y plane and the X-axis, V = L^L^L^, \ is a constant) U^ = 2R2 1, (f)^ = 2-^R3 r is selected by requiring 1 rg-i 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 electron-electron 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 8-dimensional 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 James-and-Coolidge 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 nuclear-electron 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 particle-particle 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 6-dimensional integration volume into 729 equal-volume hyper-cubes — 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 hyper-cube 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 electron-electron 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 electron-electron 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 nuclear-electron 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 nuclear-electron correlations. i*3r%-^^^''<;:^i:~^^^'^

PAGE 60

54 ,4 .2 Figure 3. N points were placed in a 6-dimensional 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 -and-error 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 5-decimal 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. 172-183, Vol. XX, Interscience Pub., 1965. 7. H. Niederreiter, "Diophantine Approximation and Its Applications." pp. 132-194, 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, one-half to one-third 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 2-electron 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 multi-dimensional quantum-mechanical 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 non-uniformly 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 M-x^ 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 multi-dimensional integral is where the x.'s are chosen uniformly at random. As N-xo, 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 equal-volume hyper-cubes (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 multi-valued w(x.)'s used herein is discussed in Section III, The case for single-valued 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 T-H 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 non-negative 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

7-1 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-.*s-riLi.>-:'^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 multi-peaked 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 /Bl|y + P2 Pb2|V = 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^l-S* ^^ ..^^;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, p-jN 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 f-j^'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(A-5:^"-'^^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 1-1 ^^ 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 p-1 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^ia---AC^^..'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 multi-valued 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)dx--leading 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 g-j 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^ = l-p-,Since the probability of selecting g-j (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(Xp|g^) + P(g2)P(Xp|g2) = 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 two-dimensional 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 sub-squares, T, 's, so that f(x,y) is a constant in each one of them, a Label the sub-squares 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 non-zero 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 T9--from the given pair of points--with 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) 5-a 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 Cl-Tn3 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 CL-rfO CD to cn S_ ro >, c "•!•I— 4J •!— CTi n3 CQ -T+-> 1— D. M-U O 0) -r-= 3 •Icn-ic I— C tc O) •I•!Cl > ^ 4-J •!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-^r-co

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 -r-u 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 'S-ICTl =:r|cn •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 Q-M CiZ--^^ V?
PAGE 91

85 and I„ = sCN(8i)(l)4)(^ + ^;^) + N(8l)(l)(l){f) + N(8f)(.9878)(i){T-J^ +,f-g|55) + 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 eleven-term James-and-Coolidge 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)) = l|g^!. (20) As an illustration, consider the calculation of the relative probability for randomly selecting spherical coordinates: r = R^M, y = 2R2-I, 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 z-axis, and ^ e the angle between the projection of r in the x-y plane and the x-axis). 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) = R-jM), 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 non-zero. Keeping the box finite allows the random numbers to be generated easily without introducing singularities. The two nuclei are fixed on the z-axis 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 nuclear-electron 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 electron-electron 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 = R-iLx, 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 =2R2-1, = R32^ (26) (r E the electron's distance from the center). Using Eqs. (20) and (25) and letting h-j(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 Hartree-Fock nuclear-electronic correlation function (in the example ^(x) is the ground state hydrogen wave function). The Hartree-Fock 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 electron-electron 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 ll-term James-andCoolidge 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, Hartree-Fock 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 non-uniform 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), McGraw-Hill, 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, McGraw-Hill, 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, twenty-four 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 .f|6078 .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 .260237E-01+ -.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 -.960838E-05 .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 .855it2E-0li .262764 .1I|7525 .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 .829535E-04 -.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 -.200719E-04 -.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 .217I|5U -.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 .501903E-01* .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. -.000t|66l|51 -.00520523 .02203itlt ltll552E-0i+ .00271^881; .00031+0835 -.000190829 -.00011+781+3 379363E-0I+ -.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 .251657E-0I+ -.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+E-01+ .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 .112198E-ni+ .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 .006i|8935 -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 -.327607E-0I+ .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 .173906E-01* .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 .711075E-04 3. a 1.2 GK 1.71265 .000816345 .0103548 -.156562 -.011778 .0734655 -.195004E-04 -.0933257 .00199497 .00179311 -.0100391 3.5082 .0718827 .0112093 -.0389181 -.0135897 .00574354 .0430364 -.00941467 -.0477697 -.113945 -.740692 -.861249 .785136E-04 3. & 1.4 GK 1.87963 .00018834 .0144051 -.136536 -.0295948 .0531575 -.347293E-04 -.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 .77126ijE-05 -.00031^9501 .158105E-0I| .153ii52E-0l+ .176164 .00208872 -. 0001+11351 -.921+127 -.00207L>01 -.936073E-0t+ .317731E-06 -.204651 -.539558 -3.03406 3.41088 2.99772 .0350205 3. & 2.2 GK .583401 -.00309683 .000540415 .00580322 .500787E-04 .421803 -.0101896 -.00355381 21684 .000378205 -.00270008 .8170554-06 -.00112048 .784774E-04 00880374 -.00203501 -1.05312 .3750441E-06 -.301038 -1.04555 -4.51092 7.21316 2.48423 .104707 3. a 2.55 GK 2.1785 -.00432415 -.974541 -.000237863 .000471424 -.000303637 -.0529544 .721503E-05 .00197343 .000177692 -.286876E-05 1.2237 .0448264 -.000689399 -.160338 .0535872 .0031392 -.29754E-05 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 -.664804E-04 .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 -.433413E-04 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 .773192E-01+ -.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 .53239E-04 .000117607 -.000483835 -.00097995 .129331 .000156927 .00326745 -.160754 -.0294859 .0228599 -.505545 .000315916 -.179892E-05 .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. 21-41, 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. 172-183, 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, McGraw-Hill, 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. 21-4K Plenum, New York, 1977. Schneider, T., Helv. Phys. Acta 42, 956 (1969). Spiegel, M. R. "Probability and Statistics." (Schaum's Outline), McGraw-Hill, 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 UTF-8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchema-instance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd
INGEST IEID EBYSOOV5K_ROAE5X INGEST_TIME 2014-12-08T22:14:56Z PACKAGE AA00026482_00001
AGREEMENT_INFO ACCOUNT UF PROJECT UFDC
FILES