Title Page
 Table of Contents
 Biased selection Monte Carlo
 Finding the best set of variational...
 Wave functions
 Discussion and results
 Checks and corrections
 Appendix: The computer code
 Biographical sketch

Title: Monte Carlo calculation of the Born-Oppenheimer potential between two helium atoms /
Full Citation
Permanent Link: http://ufdc.ufl.edu/UF00099516/00001
 Material Information
Title: Monte Carlo calculation of the Born-Oppenheimer potential between two helium atoms /
Alternate Title: Born-Oppenheimer potential between two helium atoms
Physical Description: vii, 102 leaves : ; 28 cm.
Language: English
Creator: Lowther, Rex Everett, 1951-
Publication Date: 1980
Copyright Date: 1980
Subject: Wave functions -- Measurement   ( lcsh )
Helium   ( lcsh )
Physics thesis Ph. D
Dissertations, Academic -- Physics -- UF
Genre: bibliography   ( marcgt )
non-fiction   ( marcgt )
Thesis: Thesis (Ph. D.)--University of Florida, 1980.
Bibliography: Bibliography: leaves 100-101.
General Note: Typescript.
General Note: Vita.
Statement of Responsibility: by Rex Everett Lowther.
 Record Information
Bibliographic ID: UF00099516
Volume ID: VID00001
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: alephbibnum - 000100186
oclc - 07316125
notis - AAL5647


This item has the following downloads:

PDF ( 3 MBs ) ( PDF )

Table of Contents
    Title Page
        Page i
        Page ii
    Table of Contents
        Page iii
        Page iv
        Page v
        Page vi
        Page vii
        Page 1
        Page 2
        Page 3
        Page 4
        Page 5
        Page 6
        Page 7
    Biased selection Monte Carlo
        Page 8
        Page 9
        Page 10
        Page 11
        Page 12
        Page 13
        Page 14
        Page 15
        Page 16
        Page 17
        Page 18
        Page 19
        Page 20
        Page 21
        Page 22
        Page 23
        Page 24
        Page 25
        Page 26
        Page 27
    Finding the best set of variational parameters
        Page 28
        Page 29
        Page 30
        Page 31
        Page 32
        Page 33
    Wave functions
        Page 34
        Page 35
        Page 36
        Page 37
        Page 38
        Page 39
        Page 40
        Page 41
        Page 42
        Page 43
        Page 44
        Page 45
        Page 46
        Page 47
        Page 48
        Page 49
        Page 50
        Page 51
    Discussion and results
        Page 52
        Page 53
        Page 54
        Page 55
        Page 56
        Page 57
        Page 58
        Page 59
        Page 60
        Page 61
    Checks and corrections
        Page 62
        Page 63
        Page 64
        Page 65
        Page 66
        Page 67
        Page 68
        Page 69
    Appendix: The computer code
        Page 70
        Page 71
        Page 72
        Page 73
        Page 74
        Page 75
        Page 76
        Page 77
        Page 78
        Page 79
        Page 80
        Page 81
        Page 82
        Page 83
        Page 84
        Page 85
        Page 86
        Page 87
        Page 88
        Page 89
        Page 90
        Page 91
        Page 92
        Page 93
        Page 94
        Page 95
        Page 96
        Page 97
        Page 98
        Page 99
        Page 100
        Page 101
    Biographical sketch
        Page 102
        Page 103
        Page 104
Full Text








I would like to express my deep gratitude to Dr. R. L.

Coldwell who developed the techniques that made this work

possible, who worked closely with me in all stages of this

work, and who was a constant source of guidance and en-

couragement. I would like to extend my appreciation to

Professor A. A. Broyles for his continuous support and

encouragement and whose comments and criticisms were ex-

tremely helpful in this work. In addition, I would like to

thank both Professor A. A. Broyles and R. L. Coldwell for

their careful corrections of the first draft of this dis-

sertation. I would like to thank Ron Schoenau and the

NERDC computing center for making available the hours of

CPU time necessary for the completion of this work. I

would also like to thank the NERDC operations staff (many

of whom I have become friends with) for service far beyond

the call of duty. Finally, I thank my wife, Joni, for the

use of her electrons during the preparation of this



ACKNOWLEDGMENTS. ...... . . . . . ii

ABSTRACT . . . . . . . . . . .v


I. INTRODUCTION. . . . . . . . 1

History. . . . . . . . 1
Synopsis . . . . . . . 2


General Methods. .. . ..... 8
Removal of Singularities . . .. 12
A Monte Carlo Method not Used Herein 14
Standard Deviation . . . . . 16
Calculating the Energy . . . .. 17
Weight Functions . . . . .. 23

PARAMETERS . . . . . . .. 28

Criterion Used . . . . ... 28
Selection of{X}min . . .. . 31
Optimizing Routine . . . . . 32

IV. WAVE FUNCTIONS . . . . . ... 34

General Form . . . . . ... 34
Atomic Wave Function . . . .. 34
Interaction Terms . . ....... 37
Properties of the Interaction Terms. 43
Discarded Forms for the Trial Wave
Function . . . . . . .. 49


Differencing . . . . . .. 52
Curve Fit. . . . . . . . 55
Conclusion . . . . . . .. 60


Is x-space Sampled Adequately? ... 62
Integration by Parts . . . .. 62
Comparison to the Hydrogen Molecule. 64
Difference between and the
Eigenenergy . . . . . . 65
Born-Oppenheimer Approximation . 66
Virial Theorem . . . . ... 68


REFERENCES . . . . . . . . ... . 100

BIOGRAPHICAL SKETCH. . . . . . . . .. 102

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



Rex Everett Lowther

December, 1980

Chairman: Dr. Arthur A. Broyles
Major Department: Physics

The results of the calculations of extremely accurate

wave functions for the ground state of two helium atoms,

including energies obtained from these wave functions, are

presented herein. These energies provide a variational

upper bound to the Born-Oppenheimer potential curve for

this system. The necessary expectation values were calcu-

lated by biased Monte Carlo techniques at seven internuclear

distances. The energy obtained from the trial wave function

at the potential minimum is -11.6149685 0.0000030 Ry

giving a well depth of -7.10 0.30 x 10-5 Ry at the nuclear

separation distance of 5.6 Bohr radii (aB). It is estimated

that this energy is above the energy of the exact wave func-

tion by no more than 1.8 x 10-6 Ry. The extremely small

Monte Carlo standard deviation (a) of 3.0 x 10-6 Ry was made

possible through a combination of the three factors:

1. Evaluation of the integrands for many (over 106)

Monte Carlo points. For the seven internuclear

distances this took a total of about 50 hours of

CPU time on an Amdahl 470/V6.

2. Monte Carlo methods (which allowed for analytic

removal of all singularities) for finding good

weight function.

3. The extremely accurate wave functions reported


These wave functions, in fact, were found by minimizing,

rather than the energy (), the standard deviation in this

energy (o) which is zero for a perfect wave function. This

enabled us to optimize the set of values for the 29 varia-

tional parameters by using very few Monte Carlo points and,

therefore, made this step financially feasible.

Monte Carlo evaluation of the integrals allows total

freedom to choose a natural and concise expansion for the

wave functions. The wave functions used combine Schwartz's

189-term Hylleraas-type atomic wave function with molecular

terms containing dipole-dipole, dipole-quadrupole, and

further terms in the expansion of the interatomic potential


The Born-Oppenheimer potential curve found in this work

is in rough agreement with the experimental results of

Burgmans, Farrar, and Lee (BFL). The greatest departure is

at the nuclear separation distance of 5.6 aB where the

potential found is 1.3o below the BFL result of -6.70 Ry.

Therefore the upper bound found herein should be considered

to be in agreement with the BFL potential curve with just a

hint that the exact curve is deeper than the BFL curve.



The importance and extremely small size of the well

depth has made the ground-state helium-helium pair potential

the subject of much theoretical and experimental attention

with estimates of this depth ranging from -9.1 to -13.5 K.

The recent experimental curves begin with Bruch and McGee,l

who in 1970 fitted a pair potential with a well depth of

-10.75 K to dilute-gas properties. Also in 1970, Bennewitz

et al.2 found a well depth of -10.4 K from total scattering

cross-section measurements. Differential scattering cross-

section measurements were made in 1973 by Farrar and Lee,3

who found a well depth of -11.0 0.2 K. Burgmans, Farrar,

and Lee4 (hereafter BFL) in 1976 revised this experiment,

obtaining a depth of -10.57 K. Nuclear-spin relaxation in

dilute gases has also been recognized as a sensitive means

of studying intermolecular forces. Chapman in 1975 per-

formed measurements of this on dilute helium gas, finding

a potential of the Bruch-McGee form but with a deeper well

depth of -11.5 K.

The theoretical work begins in 1931 with a paper by

Slater and Kirkwood,6 who found a helium-helium potential

(depth = 9.1 K) by joining a repulsive energy term, which


worked well for small internuclear separations, with an

attractive dipole-dipole interaction which they could

calculate at large distances. Margenau,7 also in 1931,

extended this formalism to include dipole-quadrupole and

quadrupole-quadrupole interactions. This lowered the curve

to a depth of -13.5 1.5 K, with the quadrupole-quadrupole

term accounting for only 3% of the depth at the minimum.

Configuration-interaction (CI) calculations were carried out

in 1970-1972 by McLaughlin and Schifer8 (-12.0 K),

Bertoncini and Wahl9 (-12.0 K), and others. Attempts to

correct or account for the neglect of a large (200-2000

times the size of the well depth) part of the intra-atomic

correlation energy missed by these calculations were done

by Liu and McLeanl0 (-11.0 K), Bertoncini and Wahl11

(-10.8 K), Dacre12 (-10.54 K), and Burton13 (-10.55 K). A

good discussion of the problem this correlation energy

poses to CI calculations is given in Ref. 10. A calculation

(1976) using perturbation theory was carried out by

Chalasinski and Jeziorski,14 who found a lower "bound" of

-13.4 K. When they approximately corrected for intra-

atomic correlation effects, an upper "bound" of -10.7 K

was obtained.


At the beginning of this work, innovative biased

selection Monte Carlo techniques15-17 (developed by R. L.

Coldwell in this department) were available for solving the

Schrodinger equation. Coldwell and I felt that these meth-

ods are often more precise by orders of magnitude than other
presently used Monte Carlo (MC) methods. 0 We also felt

that these techniques could compete favorably with other

existing procedures of evaluating molecular properties such

as the configuration integral (CI) method.8-1321 We

therefore set out to find (successfully) a molecular system

that, when solved, would best show:

1. The accuracy that is attainable with these new


2. The advantages of the MC method over the CI

method and other procedures.

We chose the ground state of the di-helium system because

it is a system to which a variety of techniques have been
applied by many groups and because the extremely small

well depth required a very high accuracy (1 part in 106).

Furthermore, with the existence of extremely accurate atomic
wave functions, and an appealing notion of how to con-
struct molecular wave functions from them,7 the required

accuracy seemed attainable.

The variational method most often used8-1321 (CI

methods) to solve the Schrodinger equation for electronic

systems involves the diagonalization of large N x N matrices

for wave functions of the form

'y(x) = cii(xp ) (i-i)

to find the set of N coefficients {ci} that minimizes the

energy. One problem with this method is that the choice of

gi's is limited to a set for which the integrals,

Nij = fi (i )j (x)dx (1-2)

H*ij = f j i (R)H j (k)dx ,

can be evaluated analytically. Since these O's are either

Gaussians or atomic-like orbitals, the expansion of the wave

function for systems of more than one nucleus is not a

natural one. Thus it becomes necessary to make N so large

that the matrix manipulations become difficult. For the di-

helium system, the accuracy requirements are so great that

the energies found by this method are higher than the true

energies by amounts 200-2000 times the size of the well

depth. The depth is then determined by differencing the

energies from that of the infinitely separated atoms based

on the assumption that errors in the energy cancel.

If the integrals are evaluated by Monte Carlo, the

aforementioned restrictions no longer apply. This allowed

us to choose an atomic wave function25 that included essen-

tially all (to within 5.0 x 109 Ry) of the intra-atomic

correlation energy. Since differencing from the infinite

nuclear separation distance is therefore unnecessary, the

energies calculated remain variational upper bounds. These

atomic wave functions were then appropriately combined with

a natural expansion for the molecular terms containing 29

variational parameters. These terms include dipole-dipole,

dipole-quadrupole, quadrupole-quadrupole terms, a term that

comes close to including all the terms in the above expan-

sion, and a direct electron-electron interaction term for

electrons "on" different-nuclei. The exact form for the

wave functions and all the constants necessary to duplicate

these wave functions are reported herein.

Since the parameters do not appear linearly in our

form for the wave functions, as in Eq. (1-1), they had to

be optimized by a method other than by a single matrix

inversion. Furthermore, due to the Monte Carlo uncertainty

of the energy, optimization by minimizing the energy is

quite difficult. However, as the perfect wave function is

approached, not only does the energy minimize, but the

dispersion in this energy becomes zero. Minimizing the

Monte Carlo standard deviation, a, of the energy (which is

very similar to this dispersion) over a limited number of

MC points is relatively easy because it is a sum of squares.

This allowed us to use a much smaller set of MC points (1000)

to find the wave functions than were necessary (10-5-10-6)

to find the energies of these wave functions. The simplex

method26 was then used to find the parameter set with the

lowest a. The simplex method is slow and sometimes unre-

liable at finding the absolute minimum; but, since there

were only 29 parameters to find, it still consumed much less

computer time than is needed to get the same result by

diagonalizing the large matrices of the CI method.

Monte Carlo calculations obey the relation

a = c//N (1-3)

for large N where N is the number of evaluations of the

integrand (this ranged from 377,000 to 1159,000 in the

final energy estimates) and c is a proportionality con-

stant. Trying to make this standard deviation small by

increasing N is feasible only to a certain point, however.

Beyond this point the effort is better spent at making the

proportionality constant, c, small. We made c small by a

combination of the accurate wave functions mentioned above

and the biased selection techniques used.

The biased selection method chooses each MC point with

a probability proportional to a weight function, then divides

by this weight function to correct (rigorously) for the bias

in picking it. To produce a small c, the weight function

is adjusted to sample most heavily the regions where the

integrand and variations in the integrand are large. This

work uses recent improvements in the MC method that can

choose MC points totally uncorrelated from the last point

and with more flexible weight functions than used previ-

Total energies (representing variational upper bounds

to the true energies) and kinetic energies were found for

the wave functions at the internuclear separation distances

of 4.5, 5.0, 5.6, 6.6, 7.5, 9.0, 15.0 a At the minimum


of the Born-Oppenheimer potential the energy found was

-7.10 0.30 Ry (or -11.20 0.47 K) at the N-N distance

of 5.6 aB. This result is 1.3 standard deviations below

the experimental result of -6.70 Ry obtained by Burgmans,

Farrar, and Lee4 (BFL).


General Methods

The biased selection method5-17 is used here to esti-

mate multidimensional integrals. The essence of the method

is that the multidimensional points xi are selected with

relative probability w(5i) and then corrected for this bias

so that the integrals are given as

I = f f(*)di = lim IN (2-1)


1 N f(xi)
I = (2-2)
N N i 1 w(R.i)

In its simplest form this is importance sampling.27

In one dimension it amounts to introducing a positive defi-

nite function h(x) and the variable change:

y = fx h(z)dz/ fh(z)dz, (2-3)

so that

f f(x)dx = f h(z)dz j f[x()] dy. (2-4)
0 h[x(y)]

A Monte Carlo evaluation of this integral involves choosing

Yi at random, solving Eq. (2-3) for xi = x(Yi), and finding

w(xi) = h(xi) / h(z)dz. (2-5)

This is repeated N times and the value of the integral is

found from Eq. (2-2). The probability of selecting a

given point xi is proportional to h(xi). If h(x) is large

where f(x) is large and small where it is small, points

will be concentrated in the important region of f. In ad-

dition, f(xi)/w(xi) will tend to be constant which will

reduce the Monte Carlo error. Notice that if, in Eq. (2-4),

h[x(y)] was made equal to f[x(y)], the dispersion in the

integrand would disappear. This can never be practical,

though, since that leaves the integral over h(z) in Eq.

(2-3) identical to the original integral. Nevertheless, it

is often possible to choose a weight function for which

fh(z)dz is easy to evaluate and which matches the function

f(x) closely so that the standard deviation is decreased by

several orders of magnitude.

A good example is the radial weight function used to

pick an electron position r with respect to a nucleus.

This weight function, which is similar to the square of the

Hartree Fock atomic wave function (HF( rl), closely fits

S2(R) in this region of the 3-dimensional subspace of R.
And yet, since F( Irl) is 1-dimensional, h(r) and its

integral are easy to calculate. We do this by evaluating

our best choice for the weight function at 65 values of

r = Ir, and setting up a table (Table 1) of rm, h(rm)
2 r
= rh(rm), and Im = Jm h(r)dr. The function h(r) is
m m 0

Table 1. The values of rm, h(rm), and Im used
to pick electron positions with re-
spect to the nuclei.

2 .25

h (rm)

5. 135538-2
5. 16C83E-3
2. 10743E-3
1. 19731E-3
1. 15756E-3
1. 125871-3

9. 71484E-2
1 .04701E-1
1.06530E-1 .
1.10237 --1

then defined as the linearly interpolated values between

the h Finally, an electron position is picked by the


1. Choose a random number, o, between zero and 165'

2. Determine m such that I< << Im+I

3. The radial distance is determined from the


(a I )
r = rm +1/2

+ + [h m+1 m (a I )]
m m rm+l rm (2-6)

For the purpose of using more than one weight function,

one would like a proof of Eq. (2-4) for which no change of

variables is needed. Imagine a lattice of T different

values of R., each of which contribute equally to the
integral and which has the probability w(xi) of being

selected. It must be shown that

ST N f(5.)
T il i Nm il (2-7)
i"' '-" i=l

The prime indicates a sum over each of the T points, while

the unprimed sum contains the same point many times. In a

set of points large compared to T, the point xi will be

found Nw(2.)/T times since, by assumption, its relative

probability of being selected is w(xi). Grouping the same

values of xi together, there are Nw(xi)/T terms of value

f(2.)/w(i ). The relative probability cancels to leave a

total value for each point xi of Nf(xi)/T. This is exactly

the value needed to give Eq. (2-7) explicitly. Taking the

limit T m so that all values of Ri are possible yields
Eq. (2-1). The positions xi may be picked with any relative

probability w(xi) as long as it is divided cut in Eq. (2-2).

One last step is needed to arrive at the MC method as

used here. Suppose that there are two independent methods

for choosing points and that method 1 is used with prob-

ability p, and method 2 with probability p2 = 1 pl. The

N points will contain Nplwl( i)/T terms picked with method

1. If both wl( i) and w2(qi) in Eq. (2-2) are replaced by

the overall probability of finding i.,

w(Ri) = plwl(Xi) + P2w2(Xi), (2-8)

we will have terms of total value

[Npll(xi)/T] [f(xi)/w(xi)] + [Np2w2 (i)/T] [f (5)/Ywxi)]

= Nf(x )/T (2-9)

which is exactly the value needed for Eq. (2-7).

Removal of Singularities

This last result allows us to choose points with dif-

ferent probability distributions and average over the proba-

bilities of getting final points. This made it possible to

find the weight function needed to produce small standard

deviations and to remove the singularities from H'(x).

For example, suppose we want the three dimensional

integral in the region rl < 10 of the function

f(r) = 1/|r-i] + 1/Ir -| + g(?) (2-10)

where g(r) has no singularities. The singularities would

normally introduce large dispersions owing to the fact that

points are eventually picked arbitrarily close to a and b

which requires an enormous number of points to average

down. If, however, we introduce the weight functions

w1(r) = Cl/ r-a[2, ir-a < 1 (2-11)

0, |r-a > 1

w2(r) = c2/lr-b |r-S < 1

0, |r- | > 1

w3(r) = c3, |rl < 10

0, Ir > 10

with the probabilities Pl = p2 = p3 = 1/3, we are then

required to pick these w's which choose r with equal

probability. (To prevent ever picking any points outside

the region of integration, assume \a, Il < 9.0.) The

constants c. are determined by the condition

fw(r)dr = 1 (2-12)

We find then the average weight function given by

1 1 1 3
w(r) = + 2 + 1.5 x 10-) (2-13)

which goes to infinity faster than f(r) for r close to a

or b. It can be seen from Eq. (2-2) that this eliminates

the singularities.

Suppose now that these weight functions are not aver-

aged. We would then get, instead of Eq. (2-9), terms of

total value

w(ri) = [NPwl(ri)/T] [f (ri)/wl(ri)]

+ [Np2W2(ri)/T] [f(r )/w2(r )] + . . (2-14)

However, it is still possible to pick points arbitrarily

close to a or b with the weight function w3(ri). These
3 J-
terms, since they are divided by w3(ri) rather than w(ri),

still have singularities. Terms like this with large

values of f(x) but small values of w(x) (producing high

standard deviations) are eliminated by forming a w(xi)

averaged over all possible ways of finding xi. By taking

advantage of symmetries in f(x), we were even able to

average over some values of x that produced values identical

to f(x ).

A Monte Carlo Method not Used Herein

In the biased selection method explained above, each

MC point xi is picked totally independently from all other

points. There is another method--the Markov chain

method --which picks the next point xi+l randomly within

a small box centered on the last point x.. This move is

then accepted or rejected according to a weight function

w(x) :

If w( i+l) > w(2i), the move is accepted.

If w(j i+) < w(.i), the move is accepted with

probability w( i+l)/w(i.).

When the move is rejected, Xi+l is set equal to i It can

be shown that this weight function corresponds exactly to

the weight function used in Eq. (2-2). This proof depends

upon the ergodicity condition that any position 2 can

eventually be reached from any other position.

By choosing w(x) equal to the weight function used in

the biased selection method, it is possible for the Markov

chain to achieve the same results. The singularities in

w(2), however, can cause troubles for the Markov chain.

Suppose, as in the example above, that r. is picked close
to point a such that |ri-al = 0.01. A new location picked

farther from the singularity, for example [ri+l-al = 0.10,

will have a probability

P = w( i+l)/w(i) = (c/.10)2/(c/.01)2 = 0.01 (2-15)

of being accepted. This shows that by repeatedly rejecting

points farther away from I than ri the Markov chain can

become "stuck" in these regions near the singularities.

Remedying this problem is possible, but introduces com-

plexities. Most have avoided it by leaving out the

singular terms in w. The singularities in f must then be

handled by other methods. The problem is avoided altogether

in this work because an R. chosen by the biased selection
method has no bearing on the position of the next point


Standard Deviation

The Monte Carlo standard deviation in the value found

in Eq. (2-2) can be found by squaring the equation

I-I = Z 6

and taking an ensemble average (where 6. is the error

introduced into IN by omitting the i'th point):

2 N N N 2
aN = <(I-I ) > = < 6.6.> = < X 6.>
N N i=l j=l 1 3 j=l

N 2
E 6 (2-17)
j=l 3

The nondiagonal elements cancel over an ensemble average

because each MC point is totally uncorrelated with any
previously picked point. The quantity 26., because it is
j 3
a sum of squares, is very close to its ensemble average.

In practice we found that, for 10 runs with N = 2000, this

quantity varied by around 2% and by 1% for N = 10,000. The

standard deviation in I is then

2 N 2
j=1 j

N N f(Ri) i N f(i) f() 2

= I N(N-1) w( N- 1 w(

N ( f())2 1( 2
j=l N-(N-1)2 N(N-)2 1

N [ f(x ) 2 f(i' 2
i w(x Ii1 w-xTz (2-18)

for large N. The dependence of o2 on N is readily seen
as a 1//.

Calculating the Energy

The energy of a trial wave function P (c) is

E = = f()H(R)dR (2-19)
f2 (K)di

Applying Eq. (2-2) to both integrals gives
Z (i.)HY (x.)/w ( .)
E i=l 1 (2-20)
N N 2
E (.i)/w(R i)
When summed over the same set of points xi, the numerator

and denominator are highly correlated. It can be seen from
Eq. (2-20) that this eliminates the part of the dispersion

resulting from variations in 2(2)/w(2), and leaves only

contributions from variations in HY(x)/N(2) [weighted by

2 (2)/w(R)]. Since Hy ()/ (X) becomes constant for the

correct wave function, the dispersion of EN in Eq. (2-20)

will be much smaller than the dispersions of either the

numerator or denominator. This standard deviation in EN

can be calculated by the same procedure leading to Eq.

N N E x )w
2 2 i=l
o = E 6 =
N = 3=1 y2()/

Z E.Y (x.)/w.
i=l1 1 1
N 2
Z (Ti)/wi


T 2 / i1

where Ei = HY(/i)/(5i.) is the local energy, and wi = w(Ri).

Expanding the second term to first order gives

2 N E.T2 (x)/w.
N = E N 2
=1 Z (x.)/w.
i=1 1

N 2
Z EI ( x)/w.

T2 (= .)/w.

2 2 (*j)/w.
N 2
SL2 (x)/wi

N (E.2 (.)/wj -

N 2
2 2 N 2 2
EN T2( )/w.)2/( Z T (i)/wi)
N i=l



It can be seen from this expression that when Y(x) is an

eigenfunction of H then the standard deviation is zero.

This was, in fact, the quantity minimized (see Ch. III) to

find wave functions as close as possible to the correct


Since the quantity we now want to calculate is a ratio

of two integrals, the criterion for choosing the optimum

w(x) is slightly different. In fact, it is now impossible

to eliminate the MC dispersion by choosing a weight function

that perfectly matches the integrand (since there are two

integrands). However, from this criterion alone, one can

construct a pretty good weight function by coming close to

the relation

w(x) = Y2 () (2-23)

This is the weight function most often used with the Markov
chain method. Based on the assumption that variations

of HY(R)/T(*) do not change much from region to region,

this is also the best w(x).

We have, however, found repeatedly that these variations

become greater as 2 (R) becomes small. This is because

it is always most efficient to optimize wave functions

with greatest emphasis in the regions where Y2(R) is

large. This is especially true since there is a much

greater volume in x-space where Y2 () is small than where

it is large. Consequently, it is harder to fit 2(R) to

this vast region with a limited set of variational

parameters. Therefore, we found that it is profitable to

choose a weight function that selects more points [than

the weight function in Eq. (2-23)] into the regions where

T2 () is small.

Consider the one-dimensional example over a region

where 2(x) is constant and HT(x)/g(x) is given by Fig. 1.

Out of 100 MC points, a constant weight function obeying

the relation w(x) = 2(x) = 1 will pick on the average 90

points in region 1 and 10 points in region 2. There will

be 45 points each for which HW(x)/T(x) equals 21 and 19,

and 5 points each for the values 10 and 30. Substituting

this into Eq. (2-22) gives

o2 = [45(21-20)2 + 45 20)20)2 + 5(30-20)2

+ 5(10-20) 2/(90+10)2

= (90+1000) x 104 = .109 (2-24)

Notice that over 90% of the contribution to o2 is coming

from region 2 (R2). If we now choose a weight function

that quadruples the number of points in region 2 at the

expense of Rl, we have from Eq. (2-8)

w(x) = 2/3, x in R1 (2-25)

= 4 x in R2


HW (x)/T (x)

0.9 1.0

Figure 1. Values of Hi(x)/Y(x) for an example one-
dimensional trial wave function. Positions
between 0.0 and 0.9 are in region and in
region 2 if between 0.9 and 1.0.

I / / T I / j I _



2 21-202 19-20 2 30-20 2
a = [30(2 ) + 30() 2 + 20(
2/3 2/3 4

10-20 2 1 1 2
+ 20( 4 ) ]/[60 2/ + 40 ]

= (135 +250)xl0-4 = 0.0385 (2-26)

The contribution to a2 from R1 is only slightly increased

while the contribution from R2 is cut by 75%. By throwing

extra MC points into the region where the oscillations in

Hy(x)/y(x) are large, we have found a weight function sig-

nificantly better than w(x) = 2(x).

The weight function we used for the di-helium system

was a compromise between these two considerations. About

1/4 of our MC points were picked in the region where 2 (S)

is large; the rest of the points were picked in regions

where [although y2(x) is small] oscillations in HY(i)/T(R)
were expected to be large. Although the contribution to 2

from the region where 2(() is large is higher than if

w(R) = 2 (2), the overall standard deviation is reduced

significantly. We also gain assurance that all regions of

x-space are sampled adequately.

A good measure of the number of points that are chosen

in the regions where 2(() is large is given by the quantity

we call the effective N:
N 2
[ T (. )/w.]
e N 4 2 (2-27)

If the points were chosen totally at random, one large

value would dominate both sums resulting in Ne = 1. If

the weight function in Eq. (2-23) was used, all points

would be equal in both sums giving Ne = N. As implied

above, we found that N was about 1/4 of the total number

of points N.

Weight Functions

In this work, an electron position rk in the point 5

had probabilities of being picked in the different ways:

1. With respect to nucleus A (or B). In this case

the distance rkA = rkA was chosen with a con-

stant probability for rkA < 0.30 aB [introducing

a 1/rkA into wI (or w2)] and for larger distances
k^A 1 2
with probability r2kA times the square of the

Hartree-Fock (HF) atomic wave function [introducing

the square of the HF atomic wave function into w1

(or w2)]. A slight alteration was made in the HF

wave functions to throw extra points (about 5%)

into regions where rkA [and variations in

Hy(x)/Y(x)] is large.

2. With respect to a point M midway between the

nuclei. In this case the distance rkM was chosen

with probability rkM (introducing a 1/rkM into

w3). The maximum rkM that could be picked with

this weight function (w3) was one half the inter-

nuclear separation, RAB.

3. With respect to a large box (40 x 40 x 60 aB)

centered on the nuclei. In this case the positions

rk were selected randomly within the box. The

weight function w4 is then just a constant de-

termined by the normalization condition described


4. With respect to any previously picked electron

r In this case the distance rkj = |rk-rjr was

chosen with constant probability up to a maximum

distance of 1.3 a These weight functions

(w5-w8) are then proportional to l/rkj.

The above weight functions were normalized in the manner

of Eq. (2-5) and Eq. (2-12) such that (see Table 2)

fwi(r)dr = fwj(r)dr (2-28)

An overall normalization factor is unnecessary since the

factor fh(R)di appears in both numerator and denominator

in Eq. (2-20).
The probabilities pi of picking an electron with the

weight functions w (rk) (Table 3) depend on the order in

which this electron is picked. The average weight function

for the k'th picked electron is then determined by

8 k
k k k
w(rk) Piwi(rk) (2-29)

If not averaged further, the total weight function then


Table 2. Information relevant to weight functions wi(rk) used in Eq. (2-28) and Eq.
(2-29). ..M is the modified Hartree-Fock wave function. The constant I is
M 9.6 2
defined by I = fJ wl(r)r dr and is calculated numerically within the computer
code (as explained on pages 9-11).

rk picked with normalization
i respect to: domain form for wi(rk) constants, ci
1 K. 1

1 nucleus A

rkA < 9.6 aB

gHFM(rkA) ,'

rkA > .3 aB

FM(.3)( ) rkA < .3 a

2 nucleus B

rkB < 9.6 a

similar to w above

4 1I
(40) (40) (60)

(xkl, 3Yk < 20 a B

Izkl < 30 aB

4 midpoint
between A and B

5 electron 1

6 electron 2

7 electron 3

8 electron 4

rkM < RAB/2

rk < 1.3 aB

rk2 < 1.3 aB

rk3 < 1.3 aB

rk4 < 1.3 aB

3 box


5 kl









Table 3. Probabilities Pi used in Eq. (2-29) to pick the k'th electron with the weight
functions w. (r).

Pi 1st nucleus 2nd nucleus box midpoint previously picked electrons

1st electron 85/110 5/110 3/110 17/110 0

2nd electron 5/107 85/107 3/107 7/107 7/107

3rd electron 85/114 5/114 3/114 7/114 (2x7)/114

4th electron 5/121 85/121 3/121 7/112 (3x7)/121

1 + 2 3 4 *
w() = w(r ) w(r2) w(r3) w(r4) (2-30)

However, it is not only possible to average over all weight

functions leading to the chosen value of x [see proof lead-

ing to Eq. (2-9)], but it is also correct to average over

weight functions w(x') with different positions for which

the integrands are identical such that f(x') = f(X). And,

since our wave function is antisymmetric with respect to

interchange of parallel-spin electrons, the values

P(i)HY(i) and 2 (R) are invariant with respect to these

permutations. Thus it is possible to average w(k) over

these permutations:

1 -) 2 3 4w
w(x) = ( + P2 +P34 +p1234) lw(rl) w(r2) w(r3) (r4) (2-31)

where the first two electrons picked are labeled spin-up,

the last two, spin-down. Since there are no external fields

present, the integrands also remain invariant when all elec-

tron spins are flipped. This final averaging led to the

weight function actually used,

w(x) = ( +P13P24) (1 +P2 +P34 +P12P34)

Sr 2 3 4(2-32)
w(rl) w(r2) w(r3) w(r4), (2-32)

and lowered o by about 3% from the weight function in

Eq. (2-31).


Criterion Used

After deciding on a form for the trial wave function

(Ch. IV), the optimum set of m variational parameters {c}

must be found. When all coefficients appear linearly in

T(R), as in the CI method,21 a single matrix inversion is

sufficient to solve the set of m simultaneous equations

[3 (R)H-Y()d ] = 0, k = l,m. (3-1)
3ck f2 (x)dx

This produces the lowest energy possible with the form used

and finds the best set of coefficients {c} in one step. The

terms in '(R), however, are limited to a form for which the

matrix elements in Eq. (1-2) are very precisely known. This

rules out terms that depend directly on inter-electron

distances and introduces orthogonality constraints.

Here the linear form is dropped in favor of a more

natural expansion of Y(i) with 29 variational parameters.

Our wave functions, then, are such that all integrals must

be evaluated by the Monte Carlo method and such that the

parameters {c} must be found by other than a single matrix

inversion. Also, since the MC method can only find an

estimate to the integral, computational difficulties arise

when minimizing the energy.

Suppose we were to minimize the energy by the following


1. Select a group of MC points {x}min.

2. Use a searching routine to find the parameter

set {c) which gives the lowest MC estimate to the

energy over {} min
Suppose further that cl is capable of making the local

energy, HY(x)/N(I), at the point k. become

H (xi)
= + c a (3-2)
T(x ) 1

(where a is a positive constant) and with no effect on any

other point in {X} min (Since {x} m was restricted to
min min
1000 MC points, they are spread sparsely enough to make

this a realistic situation.) The MC estimate to the energy

can be made arbitrarily small by making cl + -m, even

though we know that the correct value for cl is zero. The

true energy gets no lower, of course, because the local

energy for unsampled regions of x increases as the few

sampled regions go low. By increasing the number of points

in {X}min so that a small variation in any parameter affects

the local energy at many of these points, it is still pos-

sible to optimize T(x) by minimizing the energy. But,

since several hundred sets of variational parameters must

be tested by the optimizing routine, it is too expensive to

let {R}min become this large.

Fortunately, it is not necessary to optimize the trial

wave functions by minimizing the energy. Since all eigen-

functions of the Hamiltonian have zero dispersion in the

local energy, we were able to minimize a quantity closely

related to this dispersion--the MC standard deviation.

[The standard deviation actually used (described more fully

in Ch. V) is the uncertainty in the difference of two en-

ergies for different nuclear separation distances and is

very similar to that defined in Eq. (2-22).] A similar

local energy method was used by Frost et al.28 but without

weight functions. This forced them to waste constants

making H(/)/y(2) = in regions of k where 2(R) is


By minimizing a rather than in the above example,

the constant cl is set to the correct value of zero. This

demonstrates that a smaller number of MC points in {} min

is sufficient. Since 2 is a sum of squares, there is no

way to lower this quantity by letting a small number of

points dominate over the rest. Consequently, the points

must be fitted more equally. This greatly lowers the pos-

sible values of {cI for which '(x) has an artificially low

value of a. The minimization routine is forced to find a

wave function close to the eigenfunction in order to get a

low a--even over a relatively small set of points. We

used 1000 points in {x}mi to find the 29 variational

parameters in y(R). The obvious price we pay is that the

parameter set {c} that minimizes a is slightly different

from the {c} that gives the lowest possible upper bound to

the energy. The wave functions were so accurate, however,

that these energies were above the true energies by a very

small amount (estimated in Ch. VI).

Selection of {fmi

Since we are limited to about 1000 points in {}mi
when optimizing the parameters in (x)), it is important to

select many of these in the region where 2 (R) is large.

It is also important to choose some MC points in regions

where HN(x)/Y(x) X . Since this is more frequent when

Y2 () is small, the selection of points into the region

where 12 () is large, should not be too efficient. This

was carried out by selecting each point in {x}min out of

a group of 2 with the probability

['22 )/w(Xp )]1/2
p =p p
2 2 1/2 (3-3)
E [T (xi)/w(xi)]

This, of course, introduces a new weight function used to

calculate a:

W(X ) Pw(x ). (3-4)

The selected points were then stored, along with the weight

functions, for use by the minimization routine in testing

the various sets of variational parameters (c}.

Optimizing Routine

The simplex method26 was used to find the set of m

parameters {c} that minimized a. This routine first guesses

at m+l sets of parameters, "points," then systematically

moves them around in the m-dimensional parameter space

until the minimum a is found. The scheme takes the "high-

est point" ({c} with highest value of a) and reflects it in

the parameter space about the midpoint of all the other

"points." If this point is found to be lower than all the

other points, then another point further in this direction

is tested; if it is higher than the second highest, then

a point closer to the midpoint is selected. The lowest of

the tested points then replaces the highest previous point

and the process is repeated until a minimum is reached.

If none of the points tested above are smaller than the

second highest point, then this routine multiplies the

distance of each point from the midpoint by a random num-

ber (between -b to +b where b varies from 8 to 1/8) before

the next iteration.

As with all nonlinear optimization routines, simplex

has the two related problems:

1. It is slow. About 100 evaluations of a per

variational parameter are necessary to find the


2. The routine tends to get stuck in local minima

far above global minimum.


Several attempts29 to speed the convergence of simplex have

only increased the second problem. It is important to keep

the routine somewhat insensitive to fluctuations in a so

that a wide region of the parameter space is sampled to

prevent it from naively converging onto a local minimum.


General Form

The wave function we used combined extremely accurate

atomic wave functions with terms accounting for interac-

tions between the two atoms:

(rl,r2,r3,r4) = (1-P12-P34+P12P34)

xA (rllr3)(B (r2,r4)exp[-U(rl,r3;r2,r4)/21

where Pij is the permutation operator between parallel
electrons i and j. The term PA(rlr3) is Schwartz's

189-term Hylleraas-type atomic wave function (Table 4)

for electrons 1 and 3 on nucleus A. The function U ac-

counts for interactions between the two atoms by including

terms similar to the dipole-dipole term used by Slater.2627

Atomic Wave Function
Hylleraas,2 in 1929, obtained an upper bound of

-5.80648 Ry for the helium atom by using six terms in the


+ ,sm n -ez's
PA(rlr2) = E c s u t -(4-2)
a,m,n Z,m,n


s = (rlA+r2A)/aB'

t = (rlA-r2A)/aB,

u = rl2/aB (4-3)

The term riA = 1i- iRA is the distance between electron

i and nucleus A; z' is a constant (variational parameter);

and, to make !A(rl;r3) symmetric with respect to inter-

change of the electrons, t is found in even powers only.

A matrix inversion is then performed to find the set of

coefficients c ,m,n that give the lowest possible upper
bound with the terms used. Later authors,23-25 with the

help of fast computers, tried many variations of this idea.

Higher accuracy came with the inclusion of a greater num-
ber of terms in the expansion. Pekeris,24 using 1078 terms

in an expansion very similar to the above, and a 33-term

recursion relation, found an essentially exact answer

(-5.807448750 Ry). After this, the game became one of

getting the same answer with fewer coefficients.

We used Schwartz's25 189-term, z' = 1.75 wave func-

tion (Table 4) which provides an upper bound to the energy

of -5.80744875232 Ry. This was estimated to be above the
true energy by 2.0x10 Ry.

Although Schwartz found extreme accuracy in the ex-

pectation value for the energy, we found fluctuations in

the local energy to be much higher--by about a factor of

0.5 x106. Since a is proportional to fluctuations in the

Table 4. Parameters cZ,m n reading from left to right in
Schwartz's 189-term atomic wave function. The
"integers" (Z,m,n) are in the order: (0,0,0);
(1,0,0) (0,1,0); (3/2,0,0) (1/2,1,0); (2,0,0),
(1,1,0), (0,2,0), (0,0,2); (5/2,0,0), (3/2,1,0),
(1/2,2,0), (1/2,0,2); etc. Note s includes half-
integer powers while t includes even powers only.

-1.037848041E 1


5.3035970902E 1










7. 48460687E-1










7.933139599E-4 -2.994343269E-3
2.634949065E-4 -3.016779448E-4
1.208128981E-5 -6.751485425E-6

















local energy, and since the contribution to these fluctua-

tions from the atomic wave function is sizable (37% of a2

at the well minimum), it was profitable to use the most

accurate atomic wave function available (Schwartz's 189-

term i). These 189 coefficients were held fixed during the

minimization of a Because they were not varied, and

because iA (and its first and second derivatives) were

easy to store along with x and w(*), it was necessary to

evaluate these only once over {x}min"

Interaction Terms

The function U, which accounts for interactions between

the two atoms, may be broken into a sum of pair interac-

tions with electrons from opposite nuclei:

U(rlr ;r2,r4) = u(rl;r2) + u(rl;r4) + u(r3;r2)

+ u(r3;r4) (4-4)

The functions u(ri;r.) are given as

u(r ;r.) = E v (r i; ) + e(ri;rj) (4-5)
i ] v=O

where the v terms (defined next paragraph) are different

orders in the expansion of the interaction potential energy

between the two atoms. This interaction potential energy

can also be put into the form of Eq. (4-4):

V(lr3;r2,r4) = v(rl;2) + v(rl;r4) + v(r3;r2)

+ v(r3;r4) (4-6)


v(ri;rj) = /RAB l/riB 1/rjA + 1/rij (4-7)

Since wave functions generally have the highest amplitude

in the regions of low potential energy, we wanted somehow

to incorporate into the wave function the ability to

respond flexibly to this potential. It would have been

difficult, however, to put v(ri;r ) directly into u(ri;rj)

due to the delta functions in this term:

V 2 = -4 6(r) (4-8)

It was easy, though, to include into u(ri;rj) a func-

tion very similar to v(r;r):

v [2 -1/2 2 -1/2 2 -1/2
v0i;) = [(RAB+) (r iB+a) (rjA+t)

2 -1/2
+ (r j+a)) ]fo(riA )f(rjB) (4-9)

It can be seen that for a = 0 the expression in brackets is

equal to v(ri;r ). The variable parameter a, given in

Table 5, was added to eliminate the singularities mentioned

above and had little effect otherwise. The function f0(riA)

was included to allow the importance

with the electron-nuclear distances.

We were able to include similar

expanding the distances in Eq. (4-7)

about RAB:

of this term to vary

terms into u(ri;r.) by

in a Taylor series

2-1/2 -1/2
1 2 2 ( 2 -1/2 21 z .iA A
= [iA+YiA+(AB-ZiA R [1-2 +_]ri

r2 -1/2
1 =1 [1+2] ,
r R R 2

2 -1/2
1 1 ( iA riA+rjB
rij R [1+2 R 2


The nuclei are located on the z-axis and XiA, YiA' iA are

the x, y, z components of riA, etc. Terms multiplied by

equal powers of RAB were then grouped together. The first
surviving term (proportional to RAB) is the dipole-dipole

(D-D) interaction:

Vl(ri;rj) = (xiAjB +iAyjB 2iAjB) fl(riAfl(rjB)


Again, the function fl was added to allow for greater

flexibility in the wave function. Similarly, the dipole-
quadrupole (D-Q) term (proportional to RAB) is given as

2 2
V2 (ri;) = [riAZjB rjBZiA+ (ZiA- ZjB)

(2xiAXjB + 2YiAjB- 3ziAZjB)] 2(riA )f2(rjB)

The quadrupole-quadrupole (Q-Q) term is

2 2 2
v3rA;rB) = {6.r r. r r AjB
( iA'rjB) = 6iArjB[riAjB-riA-rjB 5iA-ZjB

2 2 2 2
15(ziArjB + jBriA

2 2 2 2
+ ZiAZjB[30(riA+jB)- 70(ZiA+jB + 105ziAjB

2 2
+ 3riArjBf3(riA)f3(rjB) (4-13)

The functions f are splines of the form

6 3
f(r) = E a.i(~- r)+ (4-14)


x, x>0
(x)+ = {0, x< (4-15)

The functions f therefore become zero for large values of

their arguments and thereby tend to make the v terms a

function of the interactions of the electrons "on" nucleus

A with those "on" nucleus B. Since the v functions appear

exponentially in Y(i), this feature is necessary for the

wave function to obey the boundary condition that, for

all i:

Lim [riA(T )] 0 (4-16)
r iA A

Finally, to allow for close encounters of electrons not

already accounted for by the atomic wave functions, the


e(ri;rj) = e(r ) = a!( -r. ) (4-17)
i i(4-17)

was added to u(ri;r ). Consider the region of i where

r.. is much smaller than all other distances. The dominant
term in the local energy, HY(:)/T(1), then becomes

HT(()/T(x) [(2/rij 7V- 2)exp(-e (r ij)/2)]exp[+e(r ij)/2]
ij 1 j ij ij

= 2/rij + 2e' (rij)/rij + e"(r i)/2 (4-18)

The derivative e'(0), therefore, was set equal to -1 to

cancel the singularities and to make the local energy more

constant. The coefficients ai, a!, and knots Ei, given in

Table 5, are parameters with respect to which the trial

wave functions were optimized (as described in Ch. III).

The knots .i were kept fixed at 2.0, 4.0, 6.0, 7.0, 8.0,

and 9.0 aB.

Table 5. Variable parameters for Eqs. (4-9), (4-14), and (4-17) which,
together with Schwartz's parameters, determine the wave functions.

S I 4.5 5-0 .6 6.6 I 7.5 I 9.0 15.0

D-D -4-70892E-2 -1.55826E-2 -4.29590E-4 7.-4165E-3 -4.20371E-2 -6.144728-2 -6.43464E-3
COEFS 4.15209E-2 I 2.22784E-2 -1.07552E-2 I 1.03047E-2 6.39287E-2 1.75601E-1 I 2.07885E-2
5.44734E-2 I 7.87565E-2 1.66395E-1 I 9.090201-2 -7.21025B-3 -3.20039E-1 1.24759E-1
-2.77328E-1 -3.225698-1 -4.37159E-1 I -3.99165E-1 -3.01037E-1 1.74799E-2 -5.01352E-1
3.16033E-1 ) 3.57847E-1 4.20820E-1 1 4.58833E-1 4.27316E-1 3.39207E-1 5.60043E-1
-1.11122E-1 I -1.25005E-1 -1.38560E-1 I -1.62970E-1 I -1.6196E-1 I -1.669298-1 -1.96674E-1

0-D -4.09052E-3 I -1.26733E-3 -2.98220E-5 8.2964BE-4 I 2.04350E-4 7.58871E-4 -
COEFS 4.41566E-3 I 2.02803E-3 -7.12454E-4 7.137053-4 I 1.88541E-3 -3.55547E-3 -----
6.37596E-3 7.980848-3 1.27601E-2 4.99981E-3 1 2.60785E-4 5.96739E-3 -
-3.07022E-2 -3.14954E-2 -3.41041E-2 -2.18727E-2 I -1.059778-2 1.26135E-3 ------
3.44949E-2 3.45461E-2 3.30412E-2 I 2.51882E-2 1 1.442568-2 -8.56932E-3 ----
S-1.20574E-2 -1.20107E-2 I -1.09097E-2 I -8.95646E-3 1 -5.40133E-3 1 3.964608-3 ----

O-Q 1 -2.48064E-3 -5.58420E-4 -1.24930E-5 2.469038-4 I -6.178948-4 I 2.51864E-4 ------
COEFs 2.14743E-3 9.23865E-4 -2.457048-4 3.87116E-4 I 6.23788E-4 -7.19779-4 -----
1.79425E-3 2.49016E-3 4.70030E-3 1.179058-3 I 7.827918-5 1.3112E-3 -- ----
-1.07968E-2 -1.16825E-2 -1.27553E-2 -7.67043E-3 -2.901008-3 -7.16191B-5 ------
I 1.2540E-2 1.33609E-2 1.24277E-2 9.61423E-3 3.90153E-3 -1.39043E-3 ------
-4.38544E-3 -4.72081E-3 -4.11347E-3 -3.52716E-3 -1.45322E-3 6.84249-4 ----------
_I_ _I_____L I_ ____ __ _I I I
0-D -5.61605E-4 -7.67956E-5 -1.548901-6 I -1.47917E-5 -4-36035E-4 ---------- -------
COEFS 274914E-4 | 7.65038E-5 -5.93223E-6 -4.83552E-5 1.074031-3 --------- -----
2.26788E-4 2.78610E-4 .434381E-4 -6.50308E-5 -3.38105E-4 ---------- -----
-1.31339E-3 -1.274168-3 -1.27742E-3 6.84898E-4 -3.778918-3 ---------- ----
1.50163E-3 1.44960E-3 1.27955E-1 -9.08712E-4 5.558128-3 ---------- -----
-5-26785E-4 -5.11162E-4 -4.28494E-4 3.39563E-4 -2.12206E-3 ------ ------

E-E -2.81922 1 1.32947E 1 I 498646E-1 -3.01473 -3.27098 -1.99296 I -1.96270
COEFS 3.20660 -8.09490E-1 -1.09985 2.39107 1.86162 1.13710 1. 13710

E-E I 1.13759 2.60407E-I 1 1.25888 9.34652E-1 5.01013-1 I 8.26368E-1 I 8.22202E-1
KNOTS 9 1.11433 I 8.378158-1 1 6.44535E-1 I 1.11414 1 7.87467E-1 I 1.20830 1 1.20830

Properties of the Interaction Terms

Consider the functions f normalized such that

u(r i;r) = 0(ri;r ) + RB(D-D terms)fl(riA )f(r)jB

+ RAB(D-Q terms)f2(riA)f2(rjB)

+ AB(Q-Q terms)f3(riA3(rjB (4-19)

Here the D-D, D-Q, and Q-Q terms are obtained directly

from the expansion of the interaction potential without

being multiplied by any constants. By comparing this to

the similar expansion of the interaction potential,

S-3 -4
v(ri;rj) = R (D-D terms) + R (D-Q terms) + .


a first guess for the functions fl, f2, f3 would have them

all roughly equal and roughly constant for small r.

Figure 2 shows several departures from this simple be-


1. The contribution of the D-D term (obtained by

adding the f0 curve to the fl curve), D-Q term,

etc. to the ratio u(ri;r.)/v(ri;r.) decreases

with increasing order. From the sole considera-

tion that u(ri;rj) lowers the energy by increasing

T(X) in the regions where the potential is low,

these ratios would be equal. The functions

1 2 3 4 5 6 7 8 aB

Figure 2. Relative values for the normalized functions fv(riA)f,(r.)
defined in Eq. (4-20) as a function of riA for the nuclear
separation distance of 5.6 aB. The radial distance rjB is
set here at the typical value of 1.0 ag. The total contribu-
tion of the dipole-dipole interaction can be obtained by
adding the fl curve to the f0 curve; similarly, the contribu-
tion of the dipole-quadrupole effects are found by adding the
f2 curve to the f curve, etc.

u(ri;r ), however, also have the effect of increas-

ing the kinetic energy. Since the kinetic energy

is proportional to the second derivative of the

wave function, smoother terms introduced into

u(ri;rj) will be more efficient at lowering the

total energy. This helps to explain why the more

oscillatory higher order terms in the expansion

of v(ri;rj) have less importance than first


2. The functions f decrease significantly from

r = 3 aB to r = 0 aB. Again, from the lone as-

sumption that u(ri;r.) is proportional to

v(ri;r ), these curves would be roughly flat.

[Of course they must become zero for large r to

satisfy the boundary condition at r = -, Eq. (4-

16).] This decrease in the electron's response

to molecular effects from the opposite atom is

apparently due to the increased influence of the

atomic effects from the nearer nucleus. Increased

shielding as r 0 from the other electron "on"

this nucleus may also be an influence in the de-

creased importance of these terms.

3. The functions f have different shapes--becoming

flatter with increasing order. It is estimated

that a was lowered by an extra 4-12% by allowing

these shapes to vary independently. This may be

related to the trend that, as the order of the

term increases, so does the value of r for which

that term has maximum effect.

Beyond the radial distance of r = 4 aB, where the wave

function is small, the exact shape of these functions has

very little importance.

The molecular terms described above and included in

the factor e-U/2 [Eq. (4-1)] account for the vast bulk of

the molecular behavior of the two-helium-atom system. Since

the Monte Carlo standard deviation, a, is a measure of the

accuracy of the wave function, the degree to which this is

true can be measured by comparing this quantity over a

fixed set of MC points. At the separation distance of

RAB = 5.6 aB, without the inclusion of any terms in U, the

standard deviation over 6400 MC points is (RAB=5.6aB,N=6400)

= 40.2x10- Ry. With the inclusion of the molecular terms,

this dispersion drops by over a factor of 10 to 3.83x10-5 Ry

--much closer to the standard deviation of the widely

separated atoms [a(RAB>15.0a ,N=6400) = 2.35x10-5 Ry].

Since the atomic coefficients [c,m,n in Eq. (4-2)] are

held fixed, and since the atomic and molecular contribu-

tions to variations in the local energy can be considered

to be roughly independent, the total standard deviation can

be written:

2 2 2
a (RAB,N) = A(N) + M(RABN)
A131 A M A


2 2
where CA and oM are the atomic and molecular contributions

to a2 and

CAtomic(N=6400) = o(RAB>15.aB,N=6400) = 2.35x10-5 Ry


for all internuclear separation distances. From this rela-

tion, we can see that a can only increase as RAB decreases

due to the increasing molecular effects. This is basically

caused by the form used for the wave function--that of two

atoms weakly interacting through multiple forces. There

will, in fact, be some point, as RAB gets smaller, where

the wave function used fails to accurately describe the

actual eigenfunction. This is due partly to the breakdown

in the accuracy of the basic form for Y(x) and partly due

to the increasing importance of higher order terms (such as

octupole-octupole) not included in U. As shown in Fig. 3,

this begins at about RAB = 5.0 aB (where the potential

"barrier" starts up).

For a long time we used Schwartz's 164-term atomic wave

function. But, with the inclusion of v0 and v3 into the
2 2
wave function, aM became small enough, compared to a,2

that it became advantageous to switch to his most accurate,

189-term, atomic wave function. Merely switching the atomic

wave function did improve the total wave function, but a

further improvement came by reoptimizing the variable

(molecular) parameters with the new atomic wave function

O(RAB, N= 6400)








Figure 3. Standard deviations in units of 10 Ry
for some trial wave functions evaluate
over 6400 Monte Carlo points. Top curve
is for wave functions optimized with
Schwartz's 164-term atomic wave function.
The middle curve is for the same wave
functions but with Schwartz's 189-term
function substituted for the 164-term
wave function. The bottom curve is for
the wave functions optimized with the
139-term atomic wave function.

_ r

(Fig. 3). It is apparent that reducing the "noise" (OA)

from the atomic effects made it easier for the optimization

routine to reduce a .

Discarded Forms for the Trial Wave Function

The first wave function we constructed has the form,16

T(R) = exp[-U(rlr2,r3,r4)/2]

(1-Pl2-P34 +12P34)A( 1,r3)B(r2,4) (4-23)

where the function U is completely symmetric:

U = Z u(r. ;r ) (4-24)
i 1j

We can see from Eq. (4-23) that the terms with the atomic

wave function 1A(r1;r3), for example, are multiplied by the

factor exp[-u( r;3)/2 u(r3;rl)/2]. But, since electrons

1 and 3 are "on" the same atom in these terms, the factor

above accounts for atomic rather than molecular effects.

This forced us to make the functions f [Eq. (4-14)] to be

short-ranged to prevent interference with the already ac-

curate atomic wave functions; and, therefore, severely

limited the flexibility of the form shown in Eq. (4-23).

Before the v0 and v3 terms were included in (x),

several extended versions of the functional form were


1. The functions fl(r) and f2(r) in the v1 and v2

terms [Eq. (4-11) and Eq. (4-12)] were generalized

to f([x2 +y2+(z-a)2 /2) where E and a were al-

lowed to vary with all the other variational


2. The vl and v2 functions in the factor

e-U/2 = exp{-[v1-v2-f(r12)]/2} were expanded

to (l-clv1/2-c2v2/2 + quadratic and cubic

terms)exp[-f(rl2)/2] and the coefficients c

allowed to vary. Incidentally, with all the

coefficients set equal to 1, the two forms

above were equivalent in the energy to 9 decimal

places. This is a good measure of the weakness

of these Van der Waals forces on the atoms.

3. Instead of using one function for all electron-

electron interactions,-as in Eq. (4-17), two e-e

functions were included--one for spin-parallel

interactions, the other for spin-antiparallel

interactions--and allowed to vary independently.

4. The factor of 2 multiplying ziAzjB in Eq. (4-11)

was allowed to vary.

5. A spline function s(zi) +s(z) was added to

u(ri;r") [Eq. (4-5)] to allow the electrons to

increase (or decrease) the probability of being

in the region between the nuclei.

None of these attempted improvements gave (x) the flexi-

bility needed at the time to lower a which, of course, was

finally provided by the inclusion of the v0 and v3 func-

tions. With much of the contribution to a eliminated by

the v0 and v3 terms, it might be argued that the improve-

ments above could more effectively reduce a further. But,

in each case, these generalizations to T(x) failed to

lower a at all. So, even with the inclusion of the v0

and v3 terms, it would be surprising if they had any

effect now.

The most likely way that (k) might be improved would

be to carry the expansion of the interaction potential en-

ergy to higher order. I believe that with the present

wave function, however, this would be more trouble than it

is worth. With the inclusion of the v0, v1, v2, v3 functions

into u, even quadrupole-octupole, v4, effects are accounted

for. Including the next function, v4, into u(ri;rj) would

free v0 to adjust for v5 effects; but, due to the very

large number of terms in v4, this would be very difficult

to code and time consuming to evaluate on a computer. An

estimate of the improvement made by the evaluation of fur-

ther terms can be made by extrapolation. By including into

the wave function the first two terms in the expansion of

v(r ;r ) [vI and v2 from Eqs. (4-11) and (4-12)],

o(RAB=5.6a ) was lowered by over a factor of 5; adding

the next two functions [v0 and v3 from Eqs. (4-9) and

(4-13)] into y(x) lowered this by not quite another factor

of 2.



For the internuclear separation distances RAB > 5.6 aB,

Fig. 3 shows that the standard deviations are not much

greater than that due to the dispersion in the atomic wave

functions. This implies that taking the difference in the

energies at two different nuclear separation distances

(RAB = R,R') will subtract out a good portion of the atomic
effects; and, as long as both energies are evaluated over

a very similar set of MC points, this will yield standard

deviations smaller than for either individual energy. The

key to actually doing this is the simple relationship

J f[g(x)]g' (x)dx = f_ f(y)dy (5-1)

which means that the free function g can be used to make

AE = (f (R,x )HY(R,)dR f- [R',g(x)]HY[R',g(*)]di (5-2)
f 2(R,x)d f2 [R' ,g(R)]dR

as smooth as possible. The function, g, is a product of

one-dimensional functions that serve to make all the spread-

ing and contracting occur in the relatively unimportant

region between the nuclei:

g(x) = g(rl)g(r2g(3)9(r3 4) (5-3)


(2+R'-R)zi/2, jzil < 1

g(ri) = g(zi) = z + (R'-R)/2, z. 1

zi (R'-R)/2, zi < -1 (5-4)

and zi is the distance of the i'th electron from the plane

bisecting the nuclei. For Izil 1, this transformation

makes the distances of g(ri) and r. to the nearest nucleus

exactly equal. Note that g(z)-z looks like a ramp. Results

from the final MC sums that determined the energies, these

energy differences, and the corresponding standard devia-

tions are reported in Table 6. One sum over 782,000 MC

points was taken with R = 5.6 aB and transformed to the

distances of R' = 6.6, 7.5, 9.0, and 15.0 ag; another was

taken over 377,000 points with R = 4.5 aB and transformed

to R' = 5.0 and 5.6 aB.

This differencing technique was also used during the

optimization of the wave functions. Since the form of the

wave function is that of two atoms weakly interacting

through molecular terms which contain all the variable

parameters, optimizing these parameters can yield a standard

deviation no less than that for the two widely separated


Table 6. Energies, energy difference s, and kinetic energies with standard
deviations in units of 10 Ry.

Ri E(Ri)

4.5 41.35 1.29

5.0 0.16 0.82

5.6 -7.16 0.33

6.6 -4.24 0.29

7.5 -2.15 0.25

9.0 -0.81 0.22

15.0 -0.16 0.21


41.19 0.68

7.33 0.34

-2.92 0.22

-2.09 0.13

-1.34 0.07

-0.65 0.07


48.52 0.93

-5.01 0.24

-3.43 0.16

-1.99 0.13


643. 25.

216. 19.

61. 16.

19. 14.

17. 12.

15. 9.


2 2 2 2
o oA+M > A (5-5)

Since optimizing the set of parameters can do nothing to
2 2
lower A, this contribution to a is an undesirable "noise"

that the minimization routine must sieve through. Most of
GA can be eliminated, however, by minimizing the standard

deviation, od' of (rather than the energy) the difference

of two energies. It can be seen from Fig. 3 that there is

an increasing (as RAB decreases) contribution to 02 that

the variable parameters cannot get rid of. Consequently,

this, too, is extraneous noise and should be subtracted

out along with aAtomic. The wave functions at each inter-

nuclear separation were, therefore, optimized by minimizing

ad with respect to the next wave function. For example,

T(RAB=4.5a ,x) was found by minimizing ad(R=4.5a ,R'=5.Oag),

T(RAB=6.6agB,), by minimizing ad(R=6.6a ,R'=7.5aB), etc.

When tested over a set of MC points independent from {x}min'

wave functions found by minimizing ad were found to have

significantly lower ad's and lower o's than wave functions

found by minimizing a.

Curve Fit

Values for the energies from Eq. (2-20), the energy

differences mentioned above, and the corresponding standard

deviations (Table 6) were used in a curve fit for minimizing

2 7 EFit(Ri)-EM(Ri) 2
X i= MC(Ri)

5 AEFit(Ri' (AR)i)- AEMC(Ri (AR) i) 2
+ l t (5-6)
i=l XoMC(Ri, (AR)i)


AE(R (AR)i) = E(Ri+(AR)i) -E(Ri) (5-7)

and Ma (Ri,(AR)) is the standard deviation of

AEMC(R ,(AR)i). The factors .i and .i made it possible to

test for the relative importance of these terms. The best

final result had error bars about 3% less than would have

been found with Xi = i. = 1. This yields

52 2
EFit(R) = EBFL(R) = c c (R-4.0)+(ki-R)2 (5-8)

with {c } = +0.558029x105, -0.592108x106, -0.140675xl0-6
-7 -9 4
-0.349129x10 -0.709027x10 Ry/aB and {k.} = 5.6, 6.6,

7.5, 9.0, 15.0 aB. The term EBFL(R) is the experimental

curves of Burgmans, Farrar, and Lee4 which was used because

it goes to the accepted30 limits at small (R < 4.0 a ) and

large (R > 10.0 aB) internuclear distances. The energies

EFit(R) and EBFL(R) for 5.0 aB < R < 9.5 aB are given in
Table 7 and Fig. 4.

Table 7. Energies with standard deviations from the curve
fit along with the experimental results of
Burgmans, Farrar, and Lee in units of 10-5 Ry.

5.0 0.36 0.26 + 0.52
5.1 -2.35 -2.51 0.48
5.2 -4.22 -4.44 0.43
5.3 -5.44 -5.73 0.38
5.4 -6.18 -6.53 0.34
5.5 -6.57 -6.96 0.31
5.6 -6.70 -7.10 0.30
5.7 -6.65 -7.04 0.29
5.8 -6.48 -6.87 0.28
5.9 -6.23 -6.61 0.27
6.0 -5.93 -6.30 0.26
6.1 -5.60 -5.94 0.25
6.2 -5.25 -5.57 0.24
6.3 -4.90 -5.20 0.23
6.4 -4.54 -4.82 0.22
6.5 -4.20 -4.46 0.21
6.6 -3.88 -4.13 0.19
6.7 -3.59 -3.82 0.18
6.8 -3.31 -3.53 0.17
6.9 -3.06 -3.27 0.16
7.0 -2.82 -3.02 0.15
7.1 -2.60 -2.79 0.14
7.2 -2.40 -2.57 0.14
7.3 -2.22 -2.38 0.13
7.4 -2.05 -2.20 0.13
7.5 -1.89 -2.04 0.12
7.6 -1.75 -1.89 0.12
7.7 -1.62 -1.75 0.12
7.8 -1.50 -1.63 0.11
7.9 -1.39 -1.51 0.10
8.0 -1.29 -1.41 0.10
8.1 -1.20 -1.31 0.09
8.2 -1.12 -1.22 0.09
8.3 -1.04 -1.13 0.08
8.4 -0.97 -1.05 0.08
8.5 -0.90 -0.98 0.08
8.6 -0.83 -0.91 0.07
8.7 -0.78 -0.84 0.07
8.8 -0.72 -0.79 0.07
8.9 -0.67 -0.74 0.07
9.0 -0.63 -0.69 0.07
9.1 -0.58 -0.65 0.07
9.2 -0.55 -0.61 0.07
9.3 -0.51 -0.57 0.07
9.4 -0.48 -0.54 0.07
9.5 -0.45 -0.51 0.07




1 -4


4 6 8 10 12 14

Figure 4. Upper bound to the Born-Oppenheimer potential (hatched
curve). The curve is two of our standard deviations
wide. Solid line is the experimental curve of Burgmans,
Farrar, and Lee.

In a manner similar to the derivation of Eq. (2-18),

the Monte Carlo standard deviation in the fit can be found

by squaring the equation

EN-Em = 6i (5-9)

and taking an ensemble average. The quantity E_ is the

true (N = m) energy for the trial wave function and 6i is

the error introduced into EN by omitting the i'th point.

This yields the formula for the standard deviation of the

energy at a nuclear separation distance R on the curve:

2 N N N 2 N 2
<(EN-E_) > =< Z <6 .> =< 6> .
i=l j=1 i=l i=l1

Our estimate for this quantity (Table 7) was found by

breaking our total run into 1159 partial runs, each rep-

resenting a point i, and redoing the curve fit successively

leaving each of these out.
Our outermost energy (-0.1600.208)xl05 Ry at

R = 15.0 aB was in good agreement with the accepted
asymptotic results (-0.027x10- Ry, as found in BFL and

references therein). Since the other theories are more

accurate in this region, our energy was made equal to this

value at this point. Through differencing, this had the

effect of slightly raising the rest of the curve. This

also reduced the standard deviation in the fit to about

equal to the standard deviation in the difference between

these energies from the energy at R = 15.0 a Consequently,

for the internuclear separations for which E(R) was evalu-

ated directly, the curve fit (Table 7) gives a better es-

timate to these energies than the direct calculations

(Table 6) since more information went into it.

Between these seven calculated energies there was

error introduced by the looping of the fitting function.

It was found that variation of the form of the fitting

function and of the knot locations ki produced variations

in E(R) between these points on the order of 0.3 standard

deviations. The directly calculated points, however, were

independent of these variations (to within 0.la). The

present fit was used because it looked smooth and it had

the qualitative features expected from the input.


The standard deviations of the curve in Fig. 4 and

Table 7 are all smaller than the total electronic energy
by factors less than 5.0x107. We have therefore demon-

strated the ability to produce extremely accurate energies

using Monte Carlo techniques. Although the accuracy

achieved here was largely due to Schwartz's extremely pre-

cise atomic wave function, chemical accuracy for most

systems is around 0.1 Ry--more than 103 times the standard

deviations reported herein.

As befits the term Monte Carlo, these calculations are

inherently risky and expensive since the best trial wave

function is found without knowing the energy bound that it

predicts. The upper bound calculated here is lower than

the BFL4 result in the region RAB > 5.1 aB with a difference

of 0.40x10-5 Ry (1.33 standard deviations) at the potential

minimum. Since this is at the 82% confidence level, this

curve should be considered in agreement with the BFL

result with just a hint that the true curve is deeper than

the BFL curve. In the region 4.5 aB < RAB < 5.0 a how-

ever, our result is above the BFL curve with a difference
of 1.07x105 Ry (0.83 standard deviations) at RAB = 4.5 a .

In fact, the calculated energy difference from Table 6

between E(RAB=4.5aB) and E(RAB=5.6aB) is 48.52 0.93x10-5 Ry

or 1.66 standard deviations above the corresponding BFL

value. This difference is probably due to the relative

inaccuracy of our wave functions in the region of small

RAB. And yet, Fig. 4 shows that this difference is unim-

portant in determining the curve: it is so steep in this

area that a difference in the energy produces an extremely

small change in the position of the potential barrier. It

is, therefore, very hard for the scattering experiments

(since they are sensitive to the overall shape of the

potential) to measure the energy in this region to this



Is x-space Sampled Adequately?

As long as the standard deviation, a, is accurate,

there is a definite bound on how far off the energy can be.

This standard deviation may, however, be artificially small

if a region of the position space x is inadequately sampled.

One way to be sure that this is not the case is to test the

convergence property [Eq. (1-3)] that a should have for a

large number, N, of MC points. If new points are picked in

previously inadequately sampled regions, then the standard

deviation for a longer run will be larger than that pre-

dicted from Eq. (1-3). This is due to a small weight func-

tion w. for these points appearing in the sums of Eq. (2-22)

To check this, the standard deviations of various energies

for ten partial runs with N = 2000 and N = 10,000 were com-

pared to the average of ten partial runs with N = 50,000.

From the runs with N = 2000 and N = 10,000, the values

predicted for a(N=50000) were equal to u(N=50000) to

within uncertainties of 2% and 1%.

Integration by Parts

An integration by parts of the integral in Eq. (2-19)

yields an integral which becomes


N 2
E [(OT(S o- (ic.l ) + VtY (x. ) ]/w.
i=l i
E= i N 21 1 (6-1)
z Y (x i)/w.

where o is the electronic gradient operator. As a check for

coding errors this expression was evaluated concurrently

with Eq. (2-20) over the same set of i.. To see the dif-
ference in Eq. (6-1) from Eq. (2-20), consider the kinetic

energy term for the simple atomic helium wave function

Y = e-( (6-2)

For Eq. (6-1) this is

OT.OT = 8e (6-3)

with nothing to cancel the integrable singularities in the

potential. The equivalent term for Eq. (2-20) is

2 -4(rl+r2)
-Y2o = 4(l/rl+ /r2 2)e ( (6-4)

with the 4/r terms canceling with potential terms making

the local energy relatively constant. This, of course,

makes Eq. (2-20) many times more precise using our methods

than Eq. (6-1). Since these two equations are different

and also weight differently the different regions of space,

agreement between them is a good test for errors and for

adequate sampling of the space of x. These values from

Eq. (6-1) are in agreement with the results from Eq. (2-20)

(Table 6) but with standard deviations on the order of

0.019 Ry. Energy differences can also be calculated using

Eq. (6-1). This lowers the standard deviation to about

0.003 Ry--still in agreement with the very accurate results

from Eq. (2-20).

Comparison to the Hydrogen Molecule

As a check on the weight functions and differencing

technique, our first 80,000 MC points at R = 6.6 aB were

scaled to the system of two infinitely separated H2

molecules. The relevant expression is


1 flAB(rlr2)2CD(r,4)HAB(l 1,2)CD (rr4)dld 2d 3dr4
2 2 2 d
AB(rlr2 )CD(rg,r4)drldr2dr3dr4
where tAB is James and Coolidge's 11-term trial H2

wave function for electrons 1 and 2 on nuclei A and B,

and H is the sum of the Hamiltonians for the two inde-

pendent H2 molecules:

H = HAB12 = HCD34 (6-6)


2 2 2 2 2 2 2 2
HABl2 V -. (6-7)
AB2 = 1 2 rA rlB r2A 2B r12 AB

Note the absence of any cross-potential terms. The Monte

Carlo estimate for this energy was

w w(+P23+P24) AB (li'2iCDr3i'r 4i (AB12

E +HCD34) AB(rli'r2i) CD (3i'4i
800001 2 2
2 = (1+P23+P24)9AB rlir2i CD r3ir4i
i=l 1


where the ri are electron positions originally picked with

respect to the 2He system, then scaled to the 2H2 system

using the differencing technique in Ch. V; and wi is the

weight function used to pick these original points multiplied

by a scaling factor. The permutations were used to include

all possible combinations with the same weight function and

had the effect of reducing the standard deviation by smooth-

ing the integral and by increasing the number of evalua-

tions. At the nuclear separations (RAB = RCD) of 1.2, 1.4,

1.5, and 1.7 aB the energies were -4.416 0.013,

-4.674 0.013, -4.612 0.014, and -4.350 0.010 eV

compared with James and Coolidge's -4.41, -4.68, -4.63,

and -4.35 eV.

Difference between and the Eigenenergy

Since our trial wave functions are not exact eigen-

functions, the variational bounds would be above the correct

eigenfunctions if the calculations were made over an in-

finite number of Monte Carlo evaluations. An estimate of

the size of this effect can be made by assuming that this
error is proportional to aM [Eq. (4-21)]. During our search

for the best wave function [lowest a (N=1000)], energies

were found for two of the earlier forms of the wave function

at RAB = 5.6 aB. A straight line fitted through the plot
of E (-5.07 0.60, -5.95 0.41, -.716 0.33x10-5 Ry) vs.

2(N=1000)(10.0, 1.8, 0.9x0-8 Ry2) for these wave functions

predicts that a perfect wave function [a (N=1000) = 0]

would have an energy lower than our present wave function

by at most 0.18xl0-5 Ry.

Born-Oppenheimer Approximation

As described in Ref. 32, we have made the Born-

Oppenheimer approximation that, from a classical viewpoint,

the motion of the electrons is the same as it would be if

the nuclei were held fixed in space. This same approxima-

tion is stated quantum mechanically by the relation,

T(X,x) = TN(X)e (X,x) (6-9)

where YN(X) and Te() are the nuclear and electronic wave

functions. In order to test the validity of this assump-

tion, this wave function can be operated on by the exact


-12 2
-in -o +V +V +V ) T (X,R) = EYN(X)T
-N- e VNN Ne ee N X)e = N (X)(,),


where oN and 0e are the nuclear and electronic gradient

operators and M is the mass of the nuclei. After expand-

ing the derivatives Eq. (6-10) becomes

-1 1 2 -1 2
{ON *e" ONN MNONee +e ( N)

+ YN(-0 1) + (VNN +VNe +Vee)Ne = ENe (6-11)

Since the nuclei are assumed to be fixed in space, the

Hamiltonian for the electrons can be defined as

H -o +V +V +V (6-12)
e e NN Ne ee


1 2
H = +H (6-13)
MN e

(Many authors separate H such that VNN is not part of He.)

The function e(X,x), therefore, is defined by the equation

H e~ (X,x) = E (X, ) (6-14)

where Ee is the electronic energy (including VNN). By

usingEqs. (6-12) and (6-14) in Eq. (6-11), and neglecting

the term in braces, Eq. (6-11) reduces to

(M2 +Ee)N = E\ (6-15)
MN e N N

This shows that, if the approximation is valid, the elec-

tronic energy Ee behaves as the potential energy between

the nuclei. This is the Born-Oppenheimer potential cal-

culated here.

It is not necessary, with the methods used here, to

make the above approximations; but it is convenient since

the corrections are negligible. These corrections have

been examined by Laue33 who finds that the leading term

which varies as RAB is of the form of an expectation value

of the interaction potential between the atoms divided by

the nuclear masses, M. Since this expectation value is of

the order of the dip in the Born-Oppenheimer potential,

dividing by M makes this correction negligible.

Virial Theorem

Kinetic energy estimates corresponding to the first

term in Eq. (6-1) were summed concurrently with Eqs.

(2-20) and (2-22). The accuracy of these calculations was

enhanced by the differencing technique described in Ch. V

and by using the kinetic energy term from Eq. (6-1) rather

than from Eq. (2-20) which has singular terms. From these

values (Table 6), it is evident that our wave functions do

not satisfy the virial theorem. For a di-atomic system in

the Born-Oppenheimer approximation, this is34

R = -2 (6-16)

where E is the Born-Oppenheimer potential, and and

are the values of the kinetic and potential energies from

the electronic wave function. By introducing a variational

scaling parameter into the wave function (Ti() from Tables

4 and 5), and minimizing the energy with respect to this

new parameter produces a new wave function,

X(x) = W(kx) (6-17)

which has a lower energy and satisfies the virial theorem.

At the potential minimum of RAB = 5.6 aB, where the term on

the left of Eq. (6-16) is zero, the values for the kinetic

and total energies from Table 6 were used to calculate a

value for k of 0.999977. This gives X(X) an energy lower

than that for T(R) by 6.0xl0-9 Ry.


The following computer code was used to perform the

Monte Carlo sums necessary to determine the energies and

kinetic energies from the wave functions already found.

The routine used to find the wave functions is very similar

to this code but with the following exceptions:

1. The simplex2 optimization routine is included

to minimize the standard deviation of the trial

wave functions.

2. A subset of MC points,{R}min is reselected from

a larger set as in Eq. (3-3).

3. Values of the weight function and the atomic

wave function and its derivatives are stored

along with these MC point positions.

C Monte Carlo summations corresponding to Egs-(2-19), (2-20), and (6-1),
C and kinetic energies from Eq (6-1) are performed. Within this loop
C (statements 4300-12700); electron positions and weight functions are
C obtained from CONFIG, transformed to the ND internuclear separation
C distances (DN), then fed into HAM where the wave function and its
C derivatives are calculated. At intervals of NTDT iterations, partial
C sums are printed and random number' table (RSEED) is changed.
IMPLICIT REAL*8(A-H,O-Z) 00000100
REAL*4 DNDMF 00000200
COMMON/POL/CP(189) ,NSSC(189) ,NTSC(189) ,NUSC (189) 00000400
DIMENSION ELPSM(5,5),EUPSM(5,5) ,EMPSH(5,5),PSM(5,5) ,PS(5) ,EPS(5) 00000600
X,ESPS(5),DN(5),XP(3,2,4),PAR(41),E(3) ,ES(5),ANOR(5),PAnP(41,5) 00000700
X ,EAV (5) ,VAR (5),EDIF(5),SIGDIF(5),X(3,2,4) 00000800
DIMENSION ET(5),EPT(5),ESPT(5),ESPF(5) ,EP (5),EU(5),PSK (5), 00000900
XPFK(5),PFKS(5) 00001000
EATA DN/5.6D0,6.600,7.5D0,9.DO,15.DO/ 00001100
IX=454314877 00001200
IY=1971698501 00001300
CALL RSEED(IX,IY) 00001400
DNN=5.6DO 00001500
CALL SETUP 00001600
ND=3 00001700
NDM=ND-1 00001800
BD=1.DO 00001900
BEAD (5,90000) NSSC,NTSC,NUSC 00002000
READ (5,10000) CP 00002100
DO 200 J=1,ND 00002200
DO 100 I=1,ND 00002300
ELPSH (I,J)=0.DO 00002400
EUPS((I,J) =0.DO 00002500
EMPSM (1,J) =0.DO 00002600

100 PSM(I,J) =.DO 00002700
PS(J)=0. DO 00002800
EPS(J)=O.DO 00002900
EPT(J)=0.DO 00003000
ESPT(J)=O.DO 00003100
ESPF(J)=0.DO 00003200
EPF(J)=O.DO 00003300
PFK(J)=O.DO 00003400
PSK(J)=O.DO 00003500
PFKS(J)=O.DO 00003600
READ (5,10000) (PABP(I,J),I=1,41) 00003700
WRITE (6,60000) (PAEP(I,J),I=1,41) 00003800
200 CONTINUE 00003900
NT=8OO 00004000
NTDT=100 00004100
N T=0 00004200
C<<<<<<<<<<<<<<<<<<<<<<< DO 1100 I=1,NT 00004400
NNT=NUT+1 00004500
DO 800 N=1,ND 00004700
00 300 K=1,41 00004800
300 PAR(K)=2ARP(K,N) 00004900
DO 400 J=1,4 00005000
DO 400 L=1,2 00005100
Do 400 K=1,3 00005200
400 XP(K,L,J)=X(K,L,J) 00005300
C<<<<<<<<<<<<<<<<<<<<<<< TBNS=1.D10 00005500
DO 500 J=1,4 00005600
ZM=X(3,1,J)-.5DO*DNN 00005700
IF (ZM.Ll.-RD) XP(3,2,J)=X (3,1,J) -DN (N) 00005800
IF (ZI.GT.RD) XP(3,1,J)=X(3,2,J)+DN(N) 00005900
IF (DABS(ZM).GT.RD) GO TO 500 00006000
DZ=1.D0+.5DO*(DN () -DNN)/RD 00006100
ZMP=DZ*ZM 00006200

TRNS=TRNS*DZ 00006300
XP (3, 1,J)=ZHP+.5)DO*DN(N) 00006400
XP(3,2,J)=ZMP-.5DO*DN(N) 00006500
500 CONTINUE 00006600
C<<<<<<<<<<<<<<<<<<<<<<< CALL HAM(XP,EK,E,ANO,PAR,DN(N)) 00006800
ES(N)=E(1) 00006900
ET(N)=E(2) 00007000
EU(N) =E(3) 00007100
IF (ANO (N).L.1.D-35) GO TO 700 00007300
EPT (N) =EPT (N) +ET (U) *ANOR (N) 00007400
ESPT (N)=ESPT (N)+ES(N)*ES(N)*ANOB(N) 00007500
EPF(N)=EPF (N) +ET(N) *ANOR(N) *ANOB(N) 00007600
ESPF(N)=ESPF(N) +ET(N)*ET(N) *ANOR(N) *ANOR(N) 00007700
PS(N)=PS (N)+ANOR(N) 000C7800
EPS( S (=EPS (N)ES(N)*ANOR(N) 00007900
PSK(N)=PSK(N)+EU(N)*ANOR(N) 00C008000
PFK (N)=PFK (N)+EU(N) *ANOR (N)*ANOR(N) 00008100
PFKS(N) =PFKS(N)+EU (N)*EU(N)*ANOH(N)*ANOR(N) 00008200
DO 600 M=1,N 00008300
IF (ANOR(M).L-.I.D-35) GO TO 600 00008400
AMN=A'NO(M) *ANOrI(N) 00008500
ELPSM (N, N) =ELPSM (M,N) +ES (M)*ANMN 00008600
EUPSM (M,N)=EUP (3, N) +ES (N) *ANMN 00008700
EilPSn (I, N) =EMPS ( ES(M, N) + ES () *ES(N)*ANN 00008800
PSM(M,N)=PSM(M,N)+A!'IN 00008900
600 CONTINUE 00009000
700 CONTINUE 00009100
800 CONTINUE 00009200
IF (NNT.NE.NTDT) GO TO 1100 00009300
NNT=0 00009400
WRITE (6,80000) PS 00009500
WRITE (6,80000) ESPT 00009600
WRITE (6,70000) IXX,IYY 00009700

C<<<<<<<<<<<<<<<<<<<<<<< DO 900 N=1,NC 00010000
EFFN=PS (N)*PS (1) /PSM(N,N) 00010100
EAV(N)=EPS (N)/PS(N) 00010200
VAR(N)=EMPSM (N,N) -2. DO*ELPS (N, N) *EAV (N) +PS (N,N) *EAV (N) *EAV(N) 00010300
SIG=DSQRT(VAB (N))/PS(N) 00010400
E2=EPT(N)/PS (N) 00010500
VAR2=ESPF(N)-2.DO*EPF(N)*E2+PSH(N,N)*E2E2 00010600
SIG2=DSQRT(VAR2)/PS (N) 00010700
EK=PSK(N)/PS(N) 00010800
VARK=PFKS(N)-2.DO*PFK(N) *EK+PSM(N,N)*EK*EK 00010900
SIGK=DSQRT(VARK) /PS (N) 00011000
WRITE (6,20000) DN(N),EAV(N),SIG,EFFN,E2,SIG2,EK,SIGK 00011100
900 CONTINUE 00011200
WRITE (6,30000) (DN(J),J=1,NDM) 00011300
DO 1000 N=1,ND 00011400
DO 1000 M=1,ND 00011500
IF (M.GE.N) GO TO 1000 00011600
EDIF (:) =EAV () -EAV (N) 00011700
VARDIF=VAR (N) *PS (M) *PS (1) +VAR (M) *PS (N) *PS (N)-2. DO*PS (M) *PS(N)* 00011800
X (EMPSM(M,N) +EAV (M *EAV (N) *PS (M,N) -ELPSI (M, N) *EAV (N) -EUPSM (a,00011900
X N)*EAV(M)) 00012000
SIGDIF(1)=DSQRT(VARDIF)/(PS(M)*PS(N)) 00012100
NM=N-1 00012200
IF (M.LT.NM) rO TO 1000 00012300
WRITE (6,40000) DN(N), (EDIF(L) ,L=1,NM) 00012400
WRITE (6,50000) (SIGDIF(L),L=1,NM) 00012,500
1000 CONTINUE 00012600
1100 CONTINUE 00012700
STOP 00012800
10000 FORMAT(4E20.10) 00012900
20000 FORBAT(3H E(,F4.1,2H)=,E19.10,3H +-,E14.7,E16.7 00013000
X ,2K,E14.7,3H +-,E11.4,E16.7,3H +-,E11-4) 00013100
30000 FORMAT(/6H DNN,5(10X,F1O.1)) 00013200
40000 FORMAT(1X,P5.1,5E20.10) 00013300
50000 FORMAT(6X,5(5X,E15-6)) 00013400

60000 FOREAT(IX,10E13.5) 00013500
70000 FORMAT(7H IX,IY=,2I14) 00013600
80000 FORMAT (SE20.12) 00013700
90000 FORMAT(7(25I3/),14I3) 00013800
END 00013900

C Input parameters are the electron positions (X), the variable
C parameters from Table 5 (PAR), and the nuclear separation distance
C (DNN). Returned values are the square of the wave function (PSI2),
C the locil kinetic energy and local energies corresponding to Eqs.
C (2-19) and (6-1) (E).
IMPLICIT REAL*8 (A-H,O-Z) 00000200
DIMENSION X(3,2,4) ,R(2,4),VNE(2,4) ,E(3),XEE(3,4,4),REE(4,4) 00000300
X,AN(2,4) ,DDAN(2,4) ,DAN (3,4,2,4) DA (3,4),DU(3,4),DJAP(3,4) 00000400
X ,R2 (2,4) ,F (2,4) 00000500
COMMOU/ALP/ALC(6),ALK(6) ,EPS,DDD 00000600
COMMON/DIP/DIPC(6) ,DIPK(6) 00000700
COMMON/QDP/QDC(6) ,QDK(6) 00000800
COMMON/QCQ/QQC(6),QQK(6) 00000900
COaMMO[/APA/EEC(3), EEK(3) 00001000
DIMENSION DQU(3,4),PAR (41),ZC(4),ZK(4), DV(3,4),AKNOT(6) 00001100
CATA AKNCT/2.DO,4.DO,6.DO,7.DO,8.DO,9.DO/ 00001200
CO 100 I = 1,6 00001300
ALC(I)=PAR(I) 00001400
DIPC(I) = PAB(I+6) 00001500
QDC(I)=PAR(I+12) 00001600
QQC(I)=PAR (1+18) 00001700
ALK(I)=AKNOT(I) 00001800
DIPK(I)=AKNOT(I) 00001900
CDK(I)=AKNOT(I) 00002000
QQK(I)=AKNOT(I) 00002100
100 CONTINUE 00002200
DO 200 I = 1,2 00002300
EEC(1) = PAR (1+24) 00002400
200 EEK(I) = PAB(I+26) 00002500
EPS=PAR(41) 00002600
tDD=DNN 00002700
V = 8.DO/DNN 00002800

DO 400 J = 1,4 00003200
DO 300 I = 1,2 00003300
R2(I,J) =X(1,I,,J)* *2+X(2,I,) **2+X(3,I,J)**2 00003400
R(I,J)=DSQRT (E2 (T,J)) 00003500
RP(I,J)=CSQR1 (R2(I,J) #EPS) 00003600
IF (R(I,J).LT..003DO) R(I,J) = .002DO 00003700
300 V = V-4.DO/R(I,J) 00003800
400 CONTINUE 00003900
DO 700 J = 2,4 00004000
Ji = J-1 00004100
DO 700 I = 1,JH 00004200
DO 500 K = 1,3 00004300
XEE(K,I,J) = X(K,1,J)-X(K,1,I) 00004400
500 XEE(K,J,I) =-XEE(K,I,J) 00004500
REE(I,J) = DSQT (XEE(1,I,J)**2+XEE(2,I,J)**2+XEE(3,I,J)**2) 00004600
IF (REE(I,J).LT..003DO) BEE(I,J) = .002D0 00004700
REE(J,I) = REE(I,J) 00004800
VEE = 2.DO/BEE(I,J) 00004900
600 V = V+VEE 00005000
700 CONTINUE 00005100
DO 1200 T = 1,2 00005400
DO 1200 J = 3,4 00005500
JU = I 00005600
JD = J 00005700
JU=I 00005900
JD=J 00006000
DO 800 L = 1,4 00006100
DO 800 K = 1,3 00006200
B00 DAN(K,L,I,J) = DJAP(K,L) 00006300
RSUM = R (1,1)+R (1,J) +R(2,3-I) +R (2,7-J) 00006400
IF (RSUM.GI.20.DO) GO TO 1200 00006500
IF (DNN.GT.10.CO) GO TO 900 00006800

CALL QIADIP (X,R,J,JD, VV,D,DDV,R2) 00007000
IF (DNN.GT.7.DO) GO TO 900 00007100
900 CONTINUE 00007300
DVCV = O.DO 00007400
CADV = O.DO 00007500
EX = 0.00 00007600
IF (VV.GT.150.) GO TO 1000 00007700
EX = DEXP(-.5DO*VV) 00007800
1000 CONTINUE 00007900
DO 1100 IE = 1,4 00008000
DO 1100 K = 1,3 00008100
DVDV = DVDV+DV(K,IE)*DV(K,IE) 00008200
1100 LAN (K,IZ,JU,JD) = EX*(DAN(K,IE,JU, JD)-.5DO*AN (JU,JD) *DV(K,IE)) 00C08400
X ) 00008600
AN(JU,JD) = EX*AN(JO,JD) 00CC8700 m
1200 CONTINUE 00008800
A = AN(1,3) -AN (2,3) -AN (1,4)+AN (2,4) 000089C0
IF (DABS(A).LT.1.E-34) A = 1.E-34 00009000
AI = 1.DO/A 00009100
DO 1300 I = 1,4 00009200
DO 1300 K = 1,3 00009300
1300 EA(K,I) = (DAN(K,I,1,3)-DAN(K,I,2,3)-DAN(K,I,1,4)+DAN(K,I,2,4))*AI00009400
EDA = (DEAN(1,3)-DDAN(2,3)-DDAN(1,4)+DDAN(2,4))*AI 00009500
C<<<<<<<<<<<<<<< LADA = O.DO 00C09900
DO 1400 I = 1,4 00010000
DO 1400 K = 1,3 00010100
1400 DADA = DA(K,I)*DA(K,I)+DADA 00010200
PSI2 = A*A 00010300
E(1) = -IEA+V 00010400
E(2) = DADAV 00010500
E(3)=DADA 00010600
BETURN 00010700

C Input values are the electron positions and radial distances (X, R,
C XEE, and REE), RP used in Eq. (4-9), and Jo and JD which label the
C permutation from Eq. (4-1). The function U from Eq. (4-1) and its
C derivatives DU and DDU are initialized to zero then contributions
C from Eq. (4-9) are added. The parameter EPS, coefficients ALC, and
C knots ALK are from Table 5.
DIMENSION X(3,2,4) ,R(2,4),DU(3,4) ,F(2,4 ) ,DF(2,4),DDF,(2,4) ,RP(2,4) 00000300
X,XEE(3,4,4) ,REE(4,4) 00000400
COMMCN/MOMN/MW (2,2) 00000600
MW(1,1) = JU 00000700
M (1,2) = JD 00000800
MW (2,1) = 3-JU 00000900
MW(2,2) = 7-JD 00001000
U = 0.DO 00001100
DDU = O.DO 00001200
DO 100 I = 1,4 00001300
DO 100 K = 1,3 00001400
100 DU(K,I) = 0.DO 00001500
DO 200 IT = 1,2 00001600
CO 200 N = 1,2 00001700
I =MW(N,II) 00001800
M=3-N 00001900
CALL SPLINE(6,ALC,ALK,R(N,I),F(N,I),DF(N,I),DDF(N,I),1) 00002000
DF(N,I)=DF(N,I)/R(N,I) 00002100
200 CONTINUE 00002200
DO 400 IA = 1,2 00002300
DO 400 IB = 1,2 00002400
I =M9(1,IA) 00002500
J =MW(2,IB) 00002600
RPE=DSQRT(REE(I,J)**2+EPS) 00002700
FF = F(1,I)*F(2,J) 00002800

ALL=1 .DO/DSQBT (DNN**24.EPS) -1. DO /RP (1 J) -1. D0/p(2,,I) + 1. DO/RPE 00002900
ADOF=0.00 00003000
00 300 K=1,3 00003100
DA1X (K,2,1)/HP(2,I)**3+XEE(K,I,J)/RPE**3 00003200
CA2=X(K,1,J)/OP(1,J)**3-XEE(K,I,J)/8PE**3 00003300
DF1=F (2,J)*DF(1,I)*X(K,1, I) 00003400
DF2=F(1.1)*DF(2,J)*X(K,2,J) 00003500
EAEF= EAI+EA1*LF +DA2*DF2 00003600
DU (R,I)= DA1+FF+ALL*DF1 +DU (K,I) 00003700
DU (K,)= DA2*FFf AT. L*F2 +0 D(K, J) 00003800
300 CONTB1OE 00003900
U=U+ALL*FF 00004000
EPA=6. DO*( (REE (I, J) /UP E) **2-lDO) /RBPE**3-3. DO* ((1((2,I) /Rp (2, 1)) **00o04loo
x 2- 1 .DO) /RP (2, I)**3+ ( (B (1 ,J)/RP (1,J) ) **2-1. DO) /EP (1,J) **3) 00004200
DDFF=F (2,J) (DDF (1 I) +2. DO *DF ( 1,I)) +F(1, I)* (DF(2,J) +2. DO*DF (2,J) ) 00 O4300
400 CONTINUE 00004500 00
BETURN 00004600 D
END 000047C0

C Dipole-quadrupole effects, Eq.(4-12), are added to U,DU, and DDU,
IMPLICIT RlEAL*8(A-H.O-Z) 00000200
COMMCN/MilO[W/W (2,2) 00000300
DIMENSION X(3,2,4) ,R(2,4),R2(2,4),DU(3,4),F(2,4),DF(2,4) ,DF(2,4) 00000400
COMMO N/QDP/DC(6) ,QDK(6) 00000500
DO 100 II = 1,2 00000600
DO 100 N = 1,2 00000700
I =MW(N,II) 00000800
CALL SPLINE(6,QDC,QDK,R(N,I),F(N,I),DF(N,I),DDF(N,I),1) 00000900
DF(N,I) = DF(N,I)/R (N,I) 00001000
100 CONTINUE 00001100
DO 300 IA = 1,2 00001200
DO 300 IB = 1,2 00001300
I =MW(1,IA) 00001400 c
J =MW(2,IB) 00001500
ZDIF = X(3,1,I)-X(3,2,J) 00001600
ZZ = X(3,1,I)*X(3,2,J) 00001700
CC = 2.DO* (X(1,1,I)*X(1,2,J)+X(2,1,I)*X(2,2,J))-3.DO*ZZ 00001800
CCP = CC-ZZ 00001900
DQ = R2 (1,)*X(3,2,J)-R2(2,J)*X (3,1,I)ZDIF*CC 00002000
FF = F(1,I)*F(2,J) 00002100
U = U+DQ*FF 00002200
DDU = DDU+F(2,J)*(DDF(1,I)*DQ+2.DO*DF(1,I)*(2.DO*DQ+X (3,1,I)*CC+ 00002300
X X(3,2,J)*R2(1,I)))+F(1,I)* (DDF(2,J)*DQ2.DO*DF (2,J)*(2.DO* 00002400
X DQ-X(3,2,J)*CC-X(3,1,I) *R2 (2,J))) 00002500
DO 200 K = 1,2 00002600
DUI(K,I) = DU(K,I)+2.DO*(X(3,2,J)*X(K,1,I)+ZDIF*X(K,2,J))*FF+DQ* 00002700
X F(2,J)*DF(1,I)*X(K, ,1I) 00002800
200 DO(K,J) = DU(K,J)+2-DO*(-X(3,1,I)*X(K.2,J)+ZDIF*X(K,1,I))*FF+DQ* 00C02900
X F(1,I)*DF(2,J)*X(K,2,J) 00003000
DU(3,1) = DU(3,I) (CCP+3.DO*X(3,2,J)*X(3,2,J)-R2(2,J))*FF+DQ*F(2, 00003100
X J) *F (1,I)*X(3,1, ) 00003200
DU(3,J) = Di[(3,J)+(-CCP-3.DO*X(3,1,I)*X(3,1,I) +R2(1,I))*FF+ C*F(1 ,0003300

X I)*DF (2,J)*X(3,2, ) 00003400
300 CONTINDE 00003500
RETURN 00003600
END 00003700



C Quadrupcle-guadrupole effects, Eq. (4-13), are added to U, DU, and DDU.
IMPLICIT REAI*8(A-H,0-Z) 00000200
COMnMON/MCNM /MW (2,2) 00000300
DIMENSION X(3,2,4) (2,4) 82(2,4 ),DU(3,4),DA(3,4),F(2,4),DF(2,4), 00000400
X DDF(2,4) ,E (3,4) 00000500
COMMON/CCQ/QQC(6),QQK(6) 00000600
DO 100 11=1,2 00000700
DO 100 N=1,2 00C00800
I=MW(N,II) 00000900
CALL SPLINE(6,QQC,QQK,R(N,I),F(N,I),D (N,I)(N ,I) F(N,I),1) 00001000
DF(N,I)=DF(N,I)/R(N,I) 00001100
100 CONTINUE 00001200
DO 300 IA=1,2 00001300
DO 300 IB=1,2 00001400 mo
I=MW ( 1,IA) 00001500 i
J=MW(2,IE) 00001600
RSUM=R (1,I)* (1,I) +R(2,J)*RB(2,J) 00001700
ZSUM=X(3,1,I)*X(3,1,I)+X(3,2,J)*X(3,2,J) 00001800
RB=R(1,I)*R(2,J) 00001900
ZZ=X(3,1,I)*X(3,2,J) 00002000
CC=X(1,1,I)*X(1,2, J)+X(2, 1,I)*X (2,2,J) +ZZ 00002100
ZDIF=X(3,1,I)-X(3,2,J) 00002200
ZDIF2=ZDIF*ZDIF 00002300
QQ=6.DO*CC*(CC-RSUNM+5.DO*ZDIF2) +3.DO*ER*RR-15.DO*(X(3,1,I)**2* 00002400
X a2(2,)+X(3,2,J)**2*2(1,I))+ZZ*(30.O*RSUM-70.DO*ZSUiM+ 00002500
X 105.DO*ZZ) 00002600
FF=F(1,I)*F(2,J) 00002700
U=0*QQ*FF 00002800
CQDF=O.DO 00002900
DO 200 K=1,2 00003000
Q(K,I)=12. CO*CC*(X(K,2,J)-X(K,1,I))+X(K,2,J) (30. DO*ZDIF2-6.DO* 00003100
X RSUM) +X(K,1,)*(6. DO*R2(2,J)+60.DO*ZZ-30.DO*X(3,2,J)**2) 00003200
DQ(K,J)=12.DO*CC*(X(K,1,I) -X(K,2,))+X(K,1,I)*(30.DO*ZDIF2-6.DO* 00003300

X RSUM) X(K,2,J)*(6.DO*R2(1,I)+60. DO*ZZ-30.DO*X(3.1,I)**2) 00003400
DU(K,I)= DQ (K,I)*FF+QQ*DF(1.I)*F(2,J)*X(K,1,I) +D[U(K,I) 00003500
DU(K,J)= DQ(K,J)*FF+QQ*DF(2,J)*F(1,I)*X(K,2,J) +DU(K,J) 00003600
DQDF=DQDF*DQ(K,I)*F(2,J)*DF(1,I)*X(K,1,I)+DQ(K,J)*F(1,I)*DF(2,J)* 00C03700
X X(K,2,J) 00003800
200 CONTINUE 00003900
DQ(3,I)=48.DO*CC*ZDIF+X(3,2,J)* (30.DO*ZDIF2+24.DO*RSUM+180.DO*Z2- 00004000
X 150-DO*X(3,1,I)**2)-70.DO*X(3,2,J)**3-24.DO*X(3,1,I)*R2 (2,J) 00004100
S 150.DO*X(3,2,J)**2)-70.DO*X(3,1,I)**3-24.DO*X(3,2,J) *B2 (1,I) 00004300
DQDF=DQDF+DQ(3,I)*F(2,J)*DF(1,I)*X(3.1,I)+DQ(3,J) *F (1,I *DF (2,J) 00004400
X X(3,2,J) 00004500
DU(3,I)= DQ (3,I)*FF+QQ*DF(I,I)*F(2,J)*X(3,1,I) +DU(3,I) 00004600
DU(3,J)= DQ(3,J)*FF+QQ*DF(2,J)*F(1,I)*X(3,2,J) +DU(3,J) 00004700
DDU= 2 DO*DQDF*QQ* (F (2,J)*(DDF(1,I)+2.DO*DF(1,I)) +F(1,I)*(E F(2,J)00004800
X +2.DO*DF(2,J))) +DDU 00004900
300 CONTINUE 00005000 C
BETURN 00005100
END 00005200

C The direct electron-electron interactions, Eq.(4-17), are added to
C U, DU, and DDIJ.
IMPLICIT REAL*8(A-H,0-Z) 00000200
COMMON/MCaW/MW (2,2) 00000300
DIMENSION XEE(3,4,4),REE(4,4),DU(3,4) 00000400
COMMCN/APA/EEC(3),EEK(3) 00000500
DO 200 IA = 1,2 00000600
DO 200 IB = 1,2 00000700
I=MW(1,IA) 00000800
J=MH(2,IB) 00000900
IF (REE(I,J).GT.3.DO) GO TO 200 00001000
DF = CF/REE(I,J) 00001200
U = UtF 00001300 oo
DDU = DDU+2.DO*DDF+4.DO*DF 00001400 0i
DO 100 K = 1,3 00001500
DU(K,I) = DU(K,I)-DF*XEE(K,I,J) 00001600
100 DU(K,J) = DU(K,J)+DF*XEE(K,I,J) 00001700
200 CONTINUE 00001800
RETURN 00001900
END 00002000

C Values for Schwartz's s-t-u expansion for the atomic wave function
C (AJAP) and its derivatives (DJAP and DDJAP) are calculated for both
C atoms. The values of 1, m, and n from Eq. (4-2) are determined ty
C NS, NT, and NU; the coefficients for these terms are given by CP.
IMPLICIT REAL*8(A-H,O-Z) 00000200
DIMENSION R(2,4),X(3,2,4),XEE(3,4,4),REE(4,4),DB(2),DDB(2), 00000300
XDS(2) ,DT(2),DII(2),DP(3,4),B(2),DJAP(3,4) 00000400
DIMENSION SP(18),UP(9) ,TP(5) 00000500
COMcON/PCL/CP(189) ,NS(189),NT(189),NU(189) 00000600
AK = 3.500 00000700
DDJAP = 0.00 00000800
AJAP = 0.00 00000900
DO 100 I 1,4 00001000
DO 100 K = 1,3 00001100 co
100 DJAP(K,I) = O.CO 00001200
RSUM = (1,JU)+R(1,JD)+R(2,3-JU)+R(2,7-JD) 00001300
IF (RSUM.G1.20.DO) GO TO 900 00001400
DUDP = 0.DO 00001500
DDU = O.DO 00001600
AK = 3.5DO 00001700
EX = DEXP(-AK*RSUM*.5DO) 00001800
DO 700 NN = 1,2 00001900
IF (NN.EQ.1) GO TO 200 00C02000
JU = 3-JI 00002100
JD = 7-JD 00002200
200 CONTINUE 00002300
E(NN) = 1.DO 00002400
CEB(Nl) = O.00 00002500
DS(NN) 0.DO 00002600
DT(NN) = 0.D0 00002700
DU(NN) = 0.CO 00002800
U = REE(JU,JD) 00002900
S = R(NN,JU)+ (NN,JD) 00003000

SRS = DSOR (S) 00003100
1 = R(NN,JU)-E(NN,JD) 00003200
IF (DABS(T).LI.1.D-25) T =1.D-25 00003300
SI = I.DO/S 00003400
TI = 1.DO/I 00003500
UI = 1.DO/U 00003600
RUPI = 1.DO/F (NN,JO) 00003700
EDNI = 1.DO/R(NN,JD) 00003800
SI2 = ST*SI 00003900
T2 = T*T 00004000
112 = TI*TI 00004100
012 = 01*UT 00004200
ZIT1 = (R(NN,JU)*R(NN,JU)-B(NN,JD)*R(NN,JD) +U*U)*EUPI*UI 00004300
ZIT2 = (R(NN,JD) (NN,JD)-R(NN,JU)*R(NN,JU)+U*U)*RDNI*UI 00004400
SP(1) = 1.-D 00004500
UP(1) = 1.DO 00004600
TP(1) = 1.DO 00004700
DO 300 I = 2,18 00C04800
300 SP(I) = SP(I-1)*SRS 00004900
DO 400 T = 2,9 00005000
400 UP(I) = UP(I-1)*U 00005100
TP(2) = 12 00005200
TP(3) = T2*T2 00005300
TP(4) = TP(3)*T2 00005400
TP(5) = TP(4)*T2 00005500
01 = SI* (RUPI+RDNI) 00005600
D2 = TI*(RUPI-RDNI) 00005700
D3 = SIOI* (2IT1+ZIT2) 00005800
D4 = TI*I*(2 I-1-ZIT2) 00005900
DO 500 NFOLY = 2,189 00006000
IP = NS(NPOLY) 00006100
PI = DFLCAT(IP)*.5DO 00006200
JP = NT(NPCLY) 00006300
KP = NU(NPOLY) 00006400
ROOT = SP(IP+1)*TP(JP/2+1) *P (KP+1)*CP(NPOLY) 00006500
B(NN) = B(NN)+ROOT 00006600

CDB(NN) = DCB (NN) +ROOT*(2. DO* (PI* (PI-1.DO) *SI2+DFLOAT (JP* (JP-1) ) 00006700
X TI2+DFLOAT (KP* (KP-1)) *UI2+I*D1+DFLOAT (JP)*D2+2-D00*DLOAT(KP) 00006800
X *OI2) +DFLCAT(KP)*PI+D3+DFLOAT(JP*KP)*D4) 00C06900
DS(NN) = DS(NS)+PI*SI*ROOT 00007000
DT(NN) = DT(NN) +DFLOAT(JP)*TI*BCOT 00007100
DU(NN) = DO (NN)+DFLOAT (KP) *UI2*EOOT 00007200
500 CONTINUE 00007300
IF (DABS (B (NN)).LT.1.E-35) B(NN) = 1.E-35 00007400
DO 600 KX = 1,3 00007500
X )/B(NN) 00007700
X )/B(NN) 00007900
600 CONTINUE 00008100
EDU=DDU2. DO*(aUPI+BDNI) 000C8200
700 CONTINUE 00008300
DUDU = 4.DO 00008400
P = B(1)*B (2) 00008500
DDP = DDB(1)/B(1)+DDB(2)/B(2) 00C08600
AJAP = EX*P 00008700
DO 800 NN = 1,2 00C08900
JU = 3-JU 00009000
JD = 7-JD 00009100
DO 800 KX = 1,3 00009200
800 EJAP(KX,JU) = AJAP*(-.5DO*X(KX,NN,JU)/R(NN,JU)*AK+DF(KX,JU)) 00009400
900 CONTINUE 00009500
BETURN 00009600
END 00009700

C <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< <<<<<<<<<<<<<< <<< <<<<< <<<<<<
C Spline functions (F) for Eq. (4-14) and its derivatives (DF and EEF)
C are calculated. The coefficients and knot positions are given by
C C and AK.
IMPLICIT REAL*8(A-H,O-Z) 00000200
F = O.DO 00000400
DF = 0.DO 00000500
DDF = O.DO 00000600
DO 100 I = 1,NK OOCO0700
IF (X.GT.AK(I).AND.NS.EQ.1) GO TO 100 00000800
IF (X.LT.AK(I) .AND.NS.EQ.-1) GO TO 100 00000900
Bi = X-AK(I) 0000100C
EDB = C(I) *EM 00001100
EB = DD*RM 00001200
F = F+DB*RM 00001300
DF = DF+EB 00001400
DDF = DDF+DDB 00001500
100 CONTINUE 00001600
DF = 3.DO*DF 00001700
DDF = 6.DO*DDF 00001800
RETUD 00001900
END 00002000

C This routine generates values for X, PB, P, and XI (Table 1) used to
C select electron positions with respect to the nuclei. The Hartree-
C Fock parameters (CHF and EHF) are found in Ref. 35. Some values from
C Table 2 are also defined.
IMPLICIT REAL*8 (A-H,O-Z) 00000200
DIMENSION PB(65) ,EHF(5) ,CHF(5) 00000400
COIMON/AUX/X(65),XI (65),P(65) ,DEPI(65) 00000500
DIMENSION C(3) 00000600
DATA CHF/.78500,.2028DO,.03693DO,-.00293DO,.00325DO/ C0000700
DATA EHF/1.43DO,2.441DO,4.10DO,6.484D0,.7978D0/ 00000800
DO 100 I = 1,5 000C0900
100 CHF(I) = DSQBT((EHF(I)*2.DO)**3/2.DO)*.282095*CHF(I) 00001000
WRITE (6,10000) 00001100
X(1) = O.D0 00001200 o
XI(1)= 0.DO 00001300
DO 300 I = 2,65 00001700
Q=.15DO*DFLOAT(I-1) 00001800
PB(I) = 0.00 00001900
DO 200 J = 1,5 00002000
200 PB(I) = PB(I)+CHF(J)*DEXP(-Q*EHF(J)) 00002100
PB(I) = PB(I)**2 00002200
P(I) = PE(I)*Q**2 00002300
P(I)=DMAX1 (P(I) ,.001D) 00002400
C(1)=3.DO 00002500
C(2)=.5DO 00002600
C(3)=3.5DO 00002700
Cl 1 =C (1)*DSQRT (C(2)) 00002800
CM=Q-C(3) 00002900
IF (QM.G7.0.DO) QM=DSQRT(QM) 00003000
P(I)=P(I)*(1.DO+C111*DEXP(-C(2)*QM*QM)) 00C03100
300 CONTINUE 00003200
P(2) = P(3) 00003300

P(1) = P(2) 00003400
DO 400 I = 2,65 00003500
In = 1-1 00003600
X(I)=.15DO*DFLOAT(I-1) 00003700
XI(I)=.075DO* (P(I) +P(IM)) +XI (IM) 00003800
DPI (I) = (P(I) -P (IM)) 15DO 00003900
WRITE (6,20000) X(I),PB(I) ,P(I) ,XI(I) 00004000
400 CONTINDE 00004100
DZ=60.DO 00004200
DX=40.DO 00004300
VOL=DX*DX*DZ 00004400
BM=.5DO*DNN 00004500
REM=1.3DO 00004600
HEEC=XI(65)/EEM 00004700
HMCON=16.DO*XI(65) /DNN**2 00004800
HCON=XI(65)*12.56637062DO/VOL 00004900
BETUR I 00005000
10000 FORMAT(3H WT FUN,11X, 11X,18X,2HPB,19X,1HP,18X,211XI) 00005100
20000 FORMAT(4E20.6) 00005200
END 00005300

C Electron positions (X), one-electron weight functions (H), and the
C final averaged weight function from Eq. (2-32) (1/WTAV) are determined.
IMPLICIT REAL*8 (A-H,O-Z) 00000200
REAL*4 RNDMF 00000300
DIMENSION X(3,2,4) ,H (8,4) HAV (4) ,CON (4) 00000500
HEC = 7.DO 00000600
DO 800 I = 1,4 00000700
NCN = 5 00000800
IF (I.EQ.1.OR.I.EQ.3) NON = 85 00000900
II = I 00001000
C<<<<<<<<<<<<<<<<<<<<<<<<< IP = 100.DO+HEC*DFLOAT(I-1) 00001400
IF (I.EQ.1) TP = 110.DO 00001500
100 CONTINUE 00001600
IW = RNDMF(1.)*TP 00001700
IF (IW.LI.90) GO TO 200 00001800
IF (IW.L'.93) GO TO 300 00001900
IF (IW.LT.100) GO TO 400 00002000
IF (I.EQ.1) GO TO 400 00002100
IE = 1.DO+RNDMF(1.)*DFLOAT(I-1) 00002200
BEE = REMH*NDIF(1.) 00002300
AC = 2.D*O+NDMF(1.)-1.D0 00002400
AS = DSQRT (1.DO-AC*AC) 00002500
PHI = 6.28318530830*aNDMF(1.) 00002600
X(1,1,I) = X(1,1,IE)+REE*AS*DCOS(PHII) 00002700
X(2,1,I) = X(2,1,IE) +REE*AS*DSIN(PHI) 00002800
X(3,1,I) = X(3,1,IE)+REE*AC 00002900
IF (DABS(X(1,1,I)) .GT..5DO*DX) GO TO 100 00003000
IF (DABS(X(2,1,I)).GT..5DO*DX) GO TO 100 00003100
IF (X(3,1,1).LI.RM-.5DO*DZ.OR.X (3,1,I).GT.RM+. 5DO*DZ) GO TO 100 00003200
GO TO 500 00003300
200 CONTINUE 000034C00

CALL XFIND(X,II) 00003500
IF (IW.GE.NCN) X(3,1,I) = X(3,1,I) DNN 00003600
GO TO 500 00003700
300 CONTINUE 00003800
X(1,1,I) = (BNDMF(1.)-.5DO)*DX 00003900
X(2,1,I) = (FNDMF(1.I-.5DO)*DX 00004000
X(3,1,I) = RNDMF(1.)*DZ+.5DO*(DNN-DZ) 00004100
GO TO 500 00004200
400 CONTINUE 00004300
FAN = RNDMF(1.) 00004400
RW = RMDSQB' (RAN) 00004500
AC = 2.DO*RNDMF(1.)-1.DO 00004600
AS = DSQRT(1.DO-AC*AC) 00004700
PHI = 6.2831e5308DO*RlDMF(1.) 00004800
X(1,1,I) = RW*AS*DCOS(PHI) 00004900
X(2,1,I) = RW*AS*DSIN(PHI) 00005000
X(3,1,I) = RW*AC+Ri 00005100
500 CONTINUE 00005200
C<<<<<<<<<<<<<<<<<<<<<<<<< C(<<<<<<<<<<<<<< X(1,2,I) = x(1,1,I) 00005700
X(2,2,I) = X(2,1,I) 00005800
X(3,2,I) = X(3,1,I)-DNN 00005900
DO 600 K = 4,8 00006000
600 H(K,I) = O.DO 000061CO
R1 = DSQBT (X(1,1,I)**2+X(2,1,I) **2+X(3,1,I)**2) 00006200
R2 = DSQBT (X(1,2,I)**2+X(2,2,I)**2+X(3,2,I)**2) 00006300
CALL FINDH(B1,TNI1) 00006400
CALL FINDH(R2,HN2) 00006500
H(1,I) = HN1 00006600
H(2,I) = HN2 00006700
RP = DSQRT(4.DO*(X(1,1,I)**2+X(2,1,I)**2)+(X(3,1,I)+X(3,2,I))**2) 00006800
IF (RP.LT.ENN) H(4,I) = HMCON/RP 00006900
H(3,I) = HCCN 00007000
CON(I) = 3.DO*H(3,1)+7.DO*H(4,I) 00007100
DO 700 IE = 1,I 00007200

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

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