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

Full Citation 
Material Information 

Title: 
Monte Carlo calculation of the BornOppenheimer potential between two helium atoms / 

Alternate Title: 
BornOppenheimer 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 
Subjects 

Subject: 
Wave functions  Measurement ( lcsh ) Helium ( lcsh ) Physics thesis Ph. D Dissertations, Academic  Physics  UF 

Genre: 
bibliography ( marcgt ) nonfiction ( marcgt ) 
Notes 

Thesis: 
Thesis (Ph. D.)University of Florida, 1980. 

Bibliography: 
Bibliography: leaves 100101. 

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 

Downloads 

Table of Contents 
Title Page
Page i
Acknowledgement
Page ii
Table of Contents
Page iii
Page iv
Abstract
Page v
Page vi
Page vii
Introduction
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
References
Page 100
Page 101
Biographical sketch
Page 102
Page 103
Page 104

Full Text 
MONTE CARLO CALCULATION OF THE
BORNOPPENHEIMER POTENTIAL BETWEEN
TWO HELIUM ATOMS
BY
REX EVERETT LOWTHER
A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL OF
THE UNIVERSITY OF FLORIDA
IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE
DEGREE OF DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
1980
ACKNOWLEDGMENTS
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
manuscript.
TABLE OF CONTENTS
Page
ACKNOWLEDGMENTS. ...... . . . . . ii
ABSTRACT . . . . . . . . . . .v
CHAPTER
I. INTRODUCTION. . . . . . . . 1
History. . . . . . . . 1
Synopsis . . . . . . . 2
II. BIASED SELECTION MONTE CARLO. . . . 8
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
III. FINDING THE BEST SET OF VARIATIONAL
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
V. DISCUSSION AND RESULTS. . . . ... 52
Differencing . . . . . .. 52
Curve Fit. . . . . . . . 55
Conclusion . . . . . . .. 60
CHECKS AND CORRECTIONS . . . . .
Is xspace Sampled Adequately? ... 62
Integration by Parts . . . .. 62
Comparison to the Hydrogen Molecule. 64
Difference between and the
Eigenenergy . . . . . . 65
BornOppenheimer Approximation . 66
Virial Theorem . . . . ... 68
APPENDIX THE COMPUTER CODE . . . . . .. 70
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
MONTE CARLO CALCULATION OF THE
BORNOPPENHEIMER POTENTIAL BETWEEN
TWO HELIUM ATOMS
By
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 BornOppenheimer 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 105 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 106 Ry. The extremely small
Monte Carlo standard deviation (a) of 3.0 x 106 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
herein.
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
189term Hylleraastype atomic wave function with molecular
terms containing dipoledipole, dipolequadrupole, and
further terms in the expansion of the interatomic potential
energy.
The BornOppenheimer 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.
CHAPTER I
INTRODUCTION
History
The importance and extremely small size of the well
depth has made the groundstate heliumhelium 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 dilutegas properties. Also in 1970, Bennewitz
et al.2 found a well depth of 10.4 K from total scattering
crosssection 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. Nuclearspin 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 BruchMcGee 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 heliumhelium potential
(depth = 9.1 K) by joining a repulsive energy term, which
1
worked well for small internuclear separations, with an
attractive dipoledipole interaction which they could
calculate at large distances. Margenau,7 also in 1931,
extended this formalism to include dipolequadrupole and
quadrupolequadrupole interactions. This lowered the curve
to a depth of 13.5 1.5 K, with the quadrupolequadrupole
term accounting for only 3% of the depth at the minimum.
Configurationinteraction (CI) calculations were carried out
in 19701972 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 (2002000
times the size of the well depth) part of the intraatomic
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.
Synopsis
At the beginning of this work, innovative biased
selection Monte Carlo techniques1517 (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
1820
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.81321 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
techniques.
2. The advantages of the MC method over the CI
method and other procedures.
We chose the ground state of the dihelium system because
it is a system to which a variety of techniques have been
114
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
2225
wave functions, and an appealing notion of how to con
67
struct molecular wave functions from them,7 the required
accuracy seemed attainable.
The variational method most often used81321 (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
N
'y(x) = cii(xp ) (ii)
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 (12)
H*ij = f j i (R)H j (k)dx ,
can be evaluated analytically. Since these O's are either
Gaussians or atomiclike 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 2002000 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 intraatomic
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 dipoledipole,
dipolequadrupole, quadrupolequadrupole terms, a term that
comes close to including all the terms in the above expan
sion, and a direct electronelectron interaction term for
electrons "on" differentnuclei. 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. (11), 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 (105106)
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 (13)
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
1820
ously.1820
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
7
of the BornOppenheimer potential the energy found was
7.10 0.30 Ry (or 11.20 0.47 K) at the NN 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).
CHAPTER II
BIASED SELECTION MONTE CARLO
General Methods
The biased selection method517 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 (21)
Nem
where
1 N f(xi)
I = (22)
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, (23)
0
so that
1
f f(x)dx = f h(z)dz j f[x()] dy. (24)
0 h[x(y)]
A Monte Carlo evaluation of this integral involves choosing
Yi at random, solving Eq. (23) for xi = x(Yi), and finding
w(xi) = h(xi) / h(z)dz. (25)
This is repeated N times and the value of the integral is
found from Eq. (22). 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. (24),
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.
(23) 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 3dimensional subspace of R.
2
And yet, since F( Irl) is 1dimensional, 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.
0.15
0.30
0.45
0.60
0.75
0.90
1.05
1.20
1.35
1.50
1.65
1.80
1.95
2.10
2.25
2 .25
2.40
2.55
2.70
2.85
3.00
3.15
3.30
3.45
3.60
3.75
3.90
4.05
4.20
4.J5
4.50
4.65
4.80
4.95
5.10
5.25
5.40
5.55
5.70
5.85
6.00
6.15
6.30
6.45
6.60
6.75
6.90
7.05
7.20
7.35
7.50
7.65
7.80
7.95
8.10
8.25
8.40
8.55
8.70
8.85
9.00
9.15
9.30
9.45
9.60
h (rm)
5. 1355382
5.13553E2
6.71723E2
7.09341E2
6.71994E2
5.98367E2
5.13377E2
4.30723E2
3.56702E2
2.93291E2
2.40218E2
1.96240E2
1.59859E2
1.29678E2
1.04546E2
8.3573723
6.60954E3
5. 16C83E3
3.97155E3
3.00809E3
2.99529E3
3.07932E3
3.11867E3
3.01786E3
2.87206E3
2.73679E3
2.61130E3
2.49487E3
2.38686E3
2.28665E3
2.19368E3
2. 10743E3
2.02741E3
1.95317E3
1.38430E3
1.82040E3
1.76112E3
1.70613E3
1.65510E3
1.60777B3
1.56385E3
1.52311E3
1.4853133
1.45025E3
1.41771E3
1.38753E3
1.35953E3
1.33355E3
1.30945E3
1.28709E3
1.26635E3
1.24710Z3
1.22925E3
1.21268E3
1. 19731E3
1.18306E3
1.16983B3
1. 15756E3
1.1461723
1.13561E3
1. 1258713
1.1167223
1.10829E3
1.10046E3
7.70329E3
1.5406622
2.4296222
3.46541E2
4.50142E2
5.45419E2
6.28799E2
6.9960722
.58664E2
8.07413E2
8.47426E2
9.80161E2
9.06868E2
9.28583E2
9.46150E2
9.60259E2
9. 71484E2
9.80312E2
9.87161E2
9.92396E2
9.96899E2
1.00145E1
1.00610E1
1.0107131
1.01512E1
1.01933E1
1.02334E1
1.02717E1
1.03083E1
1.03434E1
1.03770E1
1.04092E1
1.04402E1
1 .04701E1
1.04989E1
1.0526711
1.05535E1
1.05795E1
1.06047E1
1.06292E1
1.06530E1 .
1.06761E1
1.0698721
1.07207E1
1.0742221
1.07633E1
1.07839E1
1.08041E1
1.08239E1
1.08434E1
1.0862521
1.08814E1
1.08999E1
1.09183E1
1.09363E1
1.09542E1
1.0971811
1.09893E1
1.1006631
1.10237 1
1.1040661
1.10575E1
1.10741E1
1.10907E1
then defined as the linearly interpolated values between
the h Finally, an electron position is picked by the
recipe:
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
interpolation:
(a I )
m
r = rm +1/2
+ + [h m+1 m (a I )]
m m rm+l rm (26)
For the purpose of using more than one weight function,
one would like a proof of Eq. (24) for which no change of
variables is needed. Imagine a lattice of T different
values of R., each of which contribute equally to the
1
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 (27)
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. (27) explicitly. Taking the
limit T m so that all values of Ri are possible yields
1
Eq. (21). The positions xi may be picked with any relative
probability w(xi) as long as it is divided cut in Eq. (22).
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. (22) are replaced by
the overall probability of finding i.,
w(Ri) = plwl(Xi) + P2w2(Xi), (28)
we will have terms of total value
[Npll(xi)/T] [f(xi)/w(xi)] + [Np2w2 (i)/T] [f (5)/Ywxi)]
= Nf(x )/T (29)
which is exactly the value needed for Eq. (27).
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/ri] + 1/Ir  + g(?) (210)
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/ ra[2, ira < 1 (211)
0, ra > 1
w2(r) = c2/lrb rS < 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 (212)
We find then the average weight function given by
1 1 1 3
w(r) = + 2 + 1.5 x 10) (213)
which goes to infinity faster than f(r) for r close to a
or b. It can be seen from Eq. (22) that this eliminates
the singularities.
Suppose now that these weight functions are not aver
aged. We would then get, instead of Eq. (29), terms of
total value
w(ri) = [NPwl(ri)/T] [f (ri)/wl(ri)]
+ [Np2W2(ri)/T] [f(r )/w2(r )] + . . (214)
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 methodthe 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. (22). 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
1
to point a such that rial = 0.01. A new location picked
farther from the singularity, for example [ri+lal = 0.10,
will have a probability
P = w( i+l)/w(i) = (c/.10)2/(c/.01)2 = 0.01 (215)
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
1
method has no bearing on the position of the next point
xi+l'
Standard Deviation
The Monte Carlo standard deviation in the value found
in Eq. (22) can be found by squaring the equation
(216)
II = Z 6
i=l
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 = <(II ) > = < 6.6.> = < X 6.>
N N i=l j=l 1 3 j=l
N 2
E 6 (217)
j=l 3
The nondiagonal elements cancel over an ensemble average
because each MC point is totally uncorrelated with any
2
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(N1) w( N 1 w(
N ( f())2 1( 2
j=l N(N1)2 N(N)2 1
N [ f(x ) 2 f(i' 2
i w(x Ii1 wxTz (218)
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 (219)
f2 (K)di
Applying Eq. (22) to both integrals gives
N
Z (i.)HY (x.)/w ( .)
E i=l 1 (220)
N N 2
E (.i)/w(R i)
i=l
When summed over the same set of points xi, the numerator
and denominator are highly correlated. It can be seen from
Eq. (220) 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. (220)
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.
(218):
N
N N E x )w
2 2 i=l
o = E 6 =
N = 3=1 y2()/
i=l
N
Z E.Y (x.)/w.
i=l1 1 1
N 2
Z (Ti)/wi
i=l
2
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.
N
T2 (= .)/w.
i=l
2
2 2 (*j)/w.
N 2
SL2 (x)/wi
i=l
N
N (E.2 (.)/wj 
j=l
N 2
2 2 N 2 2
EN T2( )/w.)2/( Z T (i)/wi)
N i=l
(222)
(221)
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
eigenfunctions.
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 () (223)
This is the weight function most often used with the Markov
1820
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 xspace 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. (223)] into the regions where
T2 () is small.
Consider the onedimensional 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. (222) gives
o2 = [45(2120)2 + 45 20)20)2 + 5(3020)2
+ 5(1020) 2/(90+10)2
4
= (90+1000) x 104 = .109 (224)
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. (28)
w(x) = 2/3, x in R1 (225)
= 4 x in R2
20
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 _
J I
i
2 21202 1920 2 3020 2
a = [30(2 ) + 30() 2 + 20(
2/3 2/3 4
1020 2 1 1 2
+ 20( 4 ) ]/[60 2/ + 40 ]
4
= (135 +250)xl04 = 0.0385 (226)
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 dihelium 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)
2
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
xspace 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.]
i=l
e N 4 2 (227)
Si(xi)/w
i=l
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. (223) 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
HartreeFock (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
above.
4. With respect to any previously picked electron
r In this case the distance rkj = rkrjr was
chosen with constant probability up to a maximum
distance of 1.3 a These weight functions
(w5w8) are then proportional to l/rkj.
The above weight functions were normalized in the manner
of Eq. (25) and Eq. (212) such that (see Table 2)
fwi(r)dr = fwj(r)dr (228)
An overall normalization factor is unnecessary since the
factor fh(R)di appears in both numerator and denominator
in Eq. (220).
k
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) (229)
i=l
If not averaged further, the total weight function then
becomes
Table 2. Information relevant to weight functions wi(rk) used in Eq. (228) and Eq.
(229). ..M is the modified HartreeFock 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 911).
rk picked with normalization
i respect to: domain form for wi(rk) constants, ci
1 K. 1
1 nucleus A
rkA < 9.6 aB
2
gHFM(rkA) ,'
rkA > .3 aB
FM(.3)( ) rkA < .3 a
HFM kA
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
C4/rkM
2
C5/rkl
5 kl
C6/rk2
2
C7/r23
C8/rk4
8I/rAB
AB
I/1.3
I/1.3
I/1.3
I/1.3
k
Table 3. Probabilities Pi used in Eq. (229) to pick the k'th electron with the weight
functions w. (r).
k
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) (230)
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. (29)], 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 parallelspin 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) (231)
where the first two electrons picked are labeled spinup,
the last two, spindown. 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(232)
w(rl) w(r2) w(r3) w(r4), (232)
and lowered o by about 3% from the weight function in
Eq. (231).
CHAPTER III
FINDING THE BEST SET OF VARIATIONAL PARAMETERS
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)HY()d ] = 0, k = l,m. (31)
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. (12) are very precisely known. This
rules out terms that depend directly on interelectron
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
method:
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
mm
Suppose further that cl is capable of making the local
energy, HY(x)/N(I), at the point k. become
1
H (xi)
= + c a (32)
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 dispersionthe 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. (222).] 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
small.
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 aeven over a relatively small set of points. We
used 1000 points in {x}mi to find the 29 variational
min
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
min
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 (33)
E [T (xi)/w(xi)]
i=l
This, of course, introduces a new weight function used to
calculate a:
W(X ) Pw(x ). (34)
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 mdimensional 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
minimum.
2. The routine tends to get stuck in local minima
far above global minimum.
33
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.
CHAPTER IV
WAVE FUNCTIONS
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) = (1P12P34+P12P34)
xA (rllr3)(B (r2,r4)exp[U(rl,r3;r2,r4)/21
(41)
where Pij is the permutation operator between parallel
25
electrons i and j. The term PA(rlr3) is Schwartz's
189term Hylleraastype 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 dipoledipole term used by Slater.2627
Atomic Wave Function
22
Hylleraas,2 in 1929, obtained an upper bound of
5.80648 Ry for the helium atom by using six terms in the
form:
+ ,sm n ez's
PA(rlr2) = E c s u t (42)
a,m,n Z,m,n
where
s = (rlA+r2A)/aB'
t = (rlAr2A)/aB,
u = rl2/aB (43)
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
2325
bound with the terms used. Later authors,2325 with the
help of fast computers, tried many variations of this idea.
Higher accuracy came with the inclusion of a greater num
24
ber of terms in the expansion. Pekeris,24 using 1078 terms
in an expansion very similar to the above, and a 33term
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 189term, 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
9
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 higherby 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 189term 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.
1o00000000
1.902943036E+1
4.603636465+1*
1.037848041E 1
1.923960173E+2
2.8e53161671+2
6.698419033E+1
2.260694210E+2
1.215969591E+1
8.394132610E+1
5.3035970902E 1
9.116479191E+1
1.578020986E+2
6.162512424
6.0477680022+1
4.70399802511
1.724975592E1
1.429921553E*1
3.2384200308+1
9.943179167E2
2.029710542
5.036841074
5.94861411601
7.3150E526504
2.808997643E1
1.201623207
2.022013328E1
9.971180790E4
2.321560230E2
1.664063669E1
3.817256327E2
3.695103522E4
1.30310914725
8.616083210E4
1.015980123E2
3.060459698E3
5.101264671E5
5.047403976E6
5.676206449E1
1.244899486E*1
3.361217000
1.306100183E+1
2.770135115E+1
2.414376047E+2
2.913505181
7.999811368E+1
2.579764611E+1
1.305978447E+2
1.445389874
8.473931080E+1
3.004298892E*1
8.210878221E1
5.498596657E+1
3.120412013
8.879909237E2
2.540266224E0I
8.716991675E2
3.692937706E1
4.005407440
1.235055282E+1
8.737071941E2
3.836780819E3
7. 48460687E1
3.016330989
3.895799143E2
4.675354267E3
8.254948939E2
4.260428116E1
8.768071553E3
1.916123574E3
3.168170870E5
4.124028815E3
2.643277259E2
6.388872059E1
1.026802958
3.069398340
3.436179707
3.904835778E+1
5.531899045E+1
6.129375841
1.205385629E+2
8.255664738E1
9.848670095E+1
7.006484610E1
6.420011583E+1
6.738126236E*1
3.521437808E*+
1.059388848E+2
8.754168625
7.212391498E1
2.479212340E+1
2.511910999
1.858350881E1
7.476147452
1.721607346E1
3.245351907E1
4.364745962E3
1.459459700
7.635602081E2
1.452645159E1
5.237388480E3
1.693984877E1
1.566465228E2
3.288568040E2
2.154616496E3
1.789193139E5
8.832131764E3
1.276404132E3
7.933139599E4 2.994343269E3
2.634949065E4 3.016779448E4
1.208128981E5 6.751485425E6
3.601649402
7.632663324E1
1.469536530E*2
4.799940197
1.695078037E+1
8.110744231E41
1.724980398E1
7.982978682E+1
1.7706339232+2
1.677336023E.2
1.951982372
1.038851622E*2
1.252469963
3.366445734E11
2.452697604E+1
1.199648526
1.002670936E+I
5.029554086E*1
7.138739963
1.921603082E3
7.758764735
1.234795512
1.647217345E1
8.21523606814
1.599800986
3.677896122E1
7.594162605E2
9.668982706E4
1.951928160E1
6.109913795E2
1.793108591E2
4.103343655E4
3.132109060E7
1.066184506E2
4.347994712E3
1.714345730E3
6.119936532E5
1.161691551E7
2.018472944
6.4229609448+1
1.121025918E+2
2.399759883z+2
2.3748375568E1
4.52537247621
2.599999412E+2
1.2676170201+2
1.5943721663t2
2.4393856751+1
2.79959727681
8.683079074E*1
2.213914998
3.572976454E+1
5.679962058E+1
4.621098536E2
9.812668097
1.356885188E+1
1.066658882
2.038107176
1.646348439E+1
3.562579879
1.118357309z3
2.793244470E1
3.528033111
1.077487048
1.824394960E5
2.308777869E2
4.446743096E1
1.815735190E1
1.665615371E4
1651409751E6
8.678468784E4
2.49489490712
1.308412679E2
3.347126314E5
6.47428160497
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"
mm
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) (44)
The functions u(ri;r.) are given as
3
u(r ;r.) = E v (r i; ) + e(ri;rj) (45)
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. (44):
V(lr3;r2,r4) = v(rl;2) + v(rl;r4) + v(r3;r2)
+ v(r3;r4) (46)
where
v(ri;rj) = /RAB l/riB 1/rjA + 1/rij (47)
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) (48)
r
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) (49)
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 electronnuclear distances.
We were able to include similar
expanding the distances in Eq. (47)
about RAB:
of this term to vary
terms into u(ri;r.) by
in a Taylor series
21/2 1/2
1 2 2 ( 2 1/2 21 z .iA A
= [iA+YiA+(ABZiA R [12 +_]ri
iB B RAB RA
r2 1/2
1 =1 [1+2] ,
r R R 2
rjA RAB AB RAB
2 1/2
1 1 ( iA riA+rjB
rij R [1+2 R 2
SAB AB RAB
(410)
The nuclei are located on the zaxis 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
3
surviving term (proportional to RAB) is the dipoledipole
(DD) interaction:
Vl(ri;rj) = (xiAjB +iAyjB 2iAjB) fl(riAfl(rjB)
(411)
Again, the function fl was added to allow for greater
flexibility in the wave function. Similarly, the dipole
4
quadrupole (DQ) term (proportional to RAB) is given as
AD
2 2
V2 (ri;) = [riAZjB rjBZiA+ (ZiA ZjB)
(2xiAXjB + 2YiAjB 3ziAZjB)] 2(riA )f2(rjB)
(412)
The quadrupolequadrupole (QQ) term is
2 2 2
v3rA;rB) = {6.r r. r r AjB
( iA'rjB) = 6iArjB[riAjBriArjB 5iAZjB
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) (413)
The functions f are splines of the form
6 3
f(r) = E a.i(~ r)+ (414)
i=1
where
x, x>0
(x)+ = {0, x< (415)
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
Lim [riA(T )] 0 (416)
r iA A
Finally, to allow for close encounters of electrons not
already accounted for by the atomic wave functions, the
term
2
e(ri;rj) = e(r ) = a!( r. ) (417)
i i(417)
was added to u(ri;r ). Consider the region of i where
r.. is much smaller than all other distances. The dominant
12
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 (418)
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. (49), (414), and (417) which,
together with Schwartz's parameters, determine the wave functions.
S I 4.5 50 .6 6.6 I 7.5 I 9.0 15.0
DD 470892E2 1.55826E2 4.29590E4 7.4165E3 4.20371E2 6.1447282 6.43464E3
COEFS 4.15209E2 I 2.22784E2 1.07552E2 I 1.03047E2 6.39287E2 1.75601E1 I 2.07885E2
5.44734E2 I 7.87565E2 1.66395E1 I 9.0902012 7.21025B3 3.20039E1 1.24759E1
2.77328E1 3.2256981 4.37159E1 I 3.99165E1 3.01037E1 1.74799E2 5.01352E1
3.16033E1 ) 3.57847E1 4.20820E1 1 4.58833E1 4.27316E1 3.39207E1 5.60043E1
1.11122E1 I 1.25005E1 1.38560E1 I 1.62970E1 I 1.6196E1 I 1.6692981 1.96674E1
0D 4.09052E3 I 1.26733E3 2.98220E5 8.2964BE4 I 2.04350E4 7.58871E4 
COEFS 4.41566E3 I 2.02803E3 7.12454E4 7.1370534 I 1.88541E3 3.55547E3 
6.37596E3 7.9808483 1.27601E2 4.99981E3 1 2.60785E4 5.96739E3 
3.07022E2 3.14954E2 3.41041E2 2.18727E2 I 1.0597782 1.26135E3 
3.44949E2 3.45461E2 3.30412E2 I 2.51882E2 1 1.4425682 8.56932E3 
S1.20574E2 1.20107E2 I 1.09097E2 I 8.95646E3 1 5.40133E3 1 3.9646083 
OQ 1 2.48064E3 5.58420E4 1.24930E5 2.4690384 I 6.1789484 I 2.51864E4 
COEFs 2.14743E3 9.23865E4 2.4570484 3.87116E4 I 6.23788E4 7.197794 
1.79425E3 2.49016E3 4.70030E3 1.1790583 I 7.8279185 1.3112E3  
1.07968E2 1.16825E2 1.27553E2 7.67043E3 2.9010083 7.16191B5 
I 1.2540E2 1.33609E2 1.24277E2 9.61423E3 3.90153E3 1.39043E3 
4.38544E3 4.72081E3 4.11347E3 3.52716E3 1.45322E3 6.842494 
_I_ _I_____L I_ ____ __ _I I I
0D 5.61605E4 7.67956E5 1.5489016 I 1.47917E5 436035E4  
COEFS 274914E4  7.65038E5 5.93223E6 4.83552E5 1.0740313  
2.26788E4 2.78610E4 .434381E4 6.50308E5 3.38105E4  
1.31339E3 1.2741683 1.27742E3 6.84898E4 3.7789183  
1.50163E3 1.44960E3 1.27955E1 9.08712E4 5.5581283  
526785E4 5.11162E4 4.28494E4 3.39563E4 2.12206E3  
EE 2.81922 1 1.32947E 1 I 498646E1 3.01473 3.27098 1.99296 I 1.96270
COEFS 3.20660 8.09490E1 1.09985 2.39107 1.86162 1.13710 1. 13710
EE I 1.13759 2.60407EI 1 1.25888 9.34652E1 5.010131 I 8.26368E1 I 8.22202E1
KNOTS 9 1.11433 I 8.3781581 1 6.44535E1 I 1.11414 1 7.87467E1 I 1.20830 1 1.20830
Properties of the Interaction Terms
Consider the functions f normalized such that
33
u(r i;r) = 0(ri;r ) + RB(DD terms)fl(riA )f(r)jB
4
+ RAB(DQ terms)f2(riA)f2(rjB)
5
+ AB(QQ terms)f3(riA3(rjB (419)
Here the DD, DQ, and QQ 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,
S3 4
v(ri;rj) = R (DD terms) + R (DQ terms) + .
(420)
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
havior;
1. The contribution of the DD term (obtained by
adding the f0 curve to the fl curve), DQ 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. (420) 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 dipoledipole interaction can be obtained by
adding the fl curve to the f0 curve; similarly, the contribu
tion of the dipolequadrupole effects are found by adding the
f2 curve to the f curve, etc.
4
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
expected.
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 shapesbecoming
flatter with increasing order. It is estimated
that a was lowered by an extra 412% 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 eU/2 [Eq. (41)] account for the vast bulk of
the molecular behavior of the twoheliumatom 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.83x105 Ry
much closer to the standard deviation of the widely
separated atoms [a(RAB>15.0a ,N=6400) = 2.35x105 Ry].
Since the atomic coefficients [c,m,n in Eq. (42)] 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
(421)
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.35x105 Ry
(422)
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 functionthat 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
octupoleoctupole) 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 164term 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,
189term, 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)
10.0
8.0
6.0
4.0
2.0
10.0
15.0
5
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 164term atomic wave function.
The middle curve is for the same wave
functions but with Schwartz's 189term
function substituted for the 164term
wave function. The bottom curve is for
the wave functions optimized with the
139term 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]
(1Pl2P34 +12P34)A( 1,r3)B(r2,4) (423)
where the function U is completely symmetric:
U = Z u(r. ;r ) (424)
i 1j
We can see from Eq. (423) 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. (414)] to be
shortranged to prevent interference with the already ac
curate atomic wave functions; and, therefore, severely
limited the flexibility of the form shown in Eq. (423).
Before the v0 and v3 terms were included in (x),
several extended versions of the functional form were
tested:
1. The functions fl(r) and f2(r) in the v1 and v2
terms [Eq. (411) and Eq. (412)] were generalized
to f([x2 +y2+(za)2 /2) where E and a were al
lowed to vary with all the other variational
parameters.
2. The vl and v2 functions in the factor
eU/2 = exp{[v1v2f(r12)]/2} were expanded
to (lclv1/2c2v2/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. (417), two ee
functions were includedone for spinparallel
interactions, the other for spinantiparallel
interactionsand allowed to vary independently.
4. The factor of 2 multiplying ziAzjB in Eq. (411)
was allowed to vary.
5. A spline function s(zi) +s(z) was added to
u(ri;r") [Eq. (45)] 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 quadrupoleoctupole, 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. (411) and (412)],
o(RAB=5.6a ) was lowered by over a factor of 5; adding
the next two functions [v0 and v3 from Eqs. (49) and
(413)] into y(x) lowered this by not quite another factor
of 2.
CHAPTER V
DISCUSSION AND RESULTS
Differencing
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 (51)
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 (52)
f 2(R,x)d f2 [R' ,g(R)]dR
as smooth as possible. The function, g, is a product of
onedimensional 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) (53)
where
(2+R'R)zi/2, jzil < 1
g(ri) = g(zi) = z + (R'R)/2, z. 1
zi (R'R)/2, zi < 1 (54)
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
atoms:
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
E(Ri)E(Ri+)
41.19 0.68
7.33 0.34
2.92 0.22
2.09 0.13
1.34 0.07
0.65 0.07
E(Ri)E(Ri+2)
48.52 0.93
5.01 0.24
3.43 0.16
1.99 0.13
K(Ri)K(15)
643. 25.
216. 19.
61. 16.
19. 14.
17. 12.
15. 9.

2 2 2 2
o oA+M > A (55)
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
2
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. (220), 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 (56)
i=l XoMC(Ri, (AR)i)
where
AE(R (AR)i) = E(Ri+(AR)i) E(Ri) (57)
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 (R4.0)+(kiR)2 (58)
i=1
with {c } = +0.558029x105, 0.592108x106, 0.140675xl06
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 105 Ry.
R EBFL EMC
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
0
2
U)
iJJ
1 4
o
6
4 6 8 10 12 14
r(aB)
Figure 4. Upper bound to the BornOppenheimer 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. (218),
the Monte Carlo standard deviation in the fit can be found
by squaring the equation
N
ENEm = 6i (59)
i=l
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
<(ENE_) > =< 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.
5
Our outermost energy (0.1600.208)xl05 Ry at
R = 15.0 aB was in good agreement with the accepted
5
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.
Conclusion
The standard deviations of the curve in Fig. 4 and
Table 7 are all smaller than the total electronic energy
7
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 Rymore 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.40x105 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
5
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.93x105 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
accuracy.
CHAPTER VI
CHECKS AND CORRECTIONS
Is xspace 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. (13)] 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. (13). This is due to a small weight func
tion w. for these points appearing in the sums of Eq. (222)
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. (219)
yields an integral which becomes
62
N 2
E [(OT(S o (ic.l ) + VtY (x. ) ]/w.
i=l i
E= i N 21 1 (61)
z Y (x i)/w.
i=l
where o is the electronic gradient operator. As a check for
coding errors this expression was evaluated concurrently
with Eq. (220) over the same set of i.. To see the dif
1
ference in Eq. (61) from Eq. (220), consider the kinetic
energy term for the simple atomic helium wave function
2(rl+r2)
Y = e( (62)
For Eq. (61) this is
4(rl+r2)
OT.OT = 8e (63)
with nothing to cancel the integrable singularities in the
potential. The equivalent term for Eq. (220) is
2 4(rl+r2)
Y2o = 4(l/rl+ /r2 2)e ( (64)
with the 4/r terms canceling with potential terms making
the local energy relatively constant. This, of course,
makes Eq. (220) many times more precise using our methods
than Eq. (61). 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. (61) are in agreement with the results from Eq. (220)
(Table 6) but with standard deviations on the order of
0.019 Ry. Energy differences can also be calculated using
Eq. (61). This lowers the standard deviation to about
0.003 Rystill in agreement with the very accurate results
from Eq. (220).
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
EHH
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
(65)
31
where tAB is James and Coolidge's 11term 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 (66)
where
2 2 2 2 2 2 2 2
HABl2 V . (67)
AB2 = 1 2 rA rlB r2A 2B r12 AB
Note the absence of any crosspotential terms. The Monte
Carlo estimate for this energy was
80000
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
(68)
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
2
error is proportional to aM [Eq. (421)]. 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.33x105 Ry) vs.
2(N=1000)(10.0, 1.8, 0.9x08 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.18xl05 Ry.
BornOppenheimer 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) (69)
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
Hamiltonian:
12 2
in o +V +V +V ) T (X,R) = EYN(X)T
N e VNN Ne ee N X)e = N (X)(,),
(610)
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. (610) becomes
1 1 2 1 2
{ON *e" ONN MNONee +e ( N)
+ YN(0 1) + (VNN +VNe +Vee)Ne = ENe (611)
Since the nuclei are assumed to be fixed in space, the
Hamiltonian for the electrons can be defined as
2
H o +V +V +V (612)
e e NN Ne ee
making
1 2
H = +H (613)
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, ) (614)
where Ee is the electronic energy (including VNN). By
usingEqs. (612) and (614) in Eq. (611), and neglecting
the term in braces, Eq. (611) reduces to
(M2 +Ee)N = E\ (615)
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 BornOppenheimer 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 BornOppenheimer potential,
dividing by M makes this correction negligible.
Virial Theorem
Kinetic energy estimates corresponding to the first
term in Eq. (61) were summed concurrently with Eqs.
(220) and (222). The accuracy of these calculations was
enhanced by the differencing technique described in Ch. V
and by using the kinetic energy term from Eq. (61) rather
than from Eq. (220) which has singular terms. From these
values (Table 6), it is evident that our wave functions do
not satisfy the virial theorem. For a diatomic system in
the BornOppenheimer approximation, this is34
R = 2 (616)
AB SRA
where E is the BornOppenheimer 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) (617)
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. (616) 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.0xl09 Ry.
APPENDIX
THE COMPUTER CODE
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. (33).
3. Values of the weight function and the atomic
wave function and its derivatives are stored
along with these MC point positions.
C MAIN ROUTINE:
C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
C Monte Carlo summations corresponding to Egs(219), (220), and (61),
C and kinetic energies from Eq (61) are performed. Within this loop
C (statements 430012700); 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.
C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
IMPLICIT REAL*8(AH,OZ) 00000100
REAL*4 DNDMF 00000200
COMMON/LLL/DX,DZZ,VOL,RM,REM,HEEC,HMCON,HCON,DNN 00000300
COMMON/POL/CP(189) ,NSSC(189) ,NTSC(189) ,NUSC (189) 00000400
COMMON/RAN/IXX,IYY 00000500
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=ND1 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
CALL CONFIG(X,HINV,NNT) 00004600
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
ANOR(N)=HINV*TBNS*ANO 00007200
IF (ANO (N).L.1.D35) 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.D35) 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
CALL RSEED(IXX,IYY) 00009800
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=N1 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 +,E114) 00013100
30000 FORMAT(/6H DNN,5(10X,F1O.1)) 00013200
40000 FORMAT(1X,P5.1,5E20.10) 00013300
50000 FORMAT(6X,5(5X,E156)) 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
SUBROUTINE HAM (X,EK,E,PSI2,PAR,DNN) 00000100
C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
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 (219) and (61) (E).
C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
IMPLICIT REAL*8 (AH,OZ) 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
C<<<<<<<<<<
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 = V4.DO/R(I,J) 00003800
400 CONTINUE 00003900
DO 700 J = 2,4 00004000
Ji = J1 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
CALL POLY(X.R,XEE,REE,NNT,JU,JD,AN(I,J),DJAP,DDAN(I,J)) 00005800
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,3I) +R (2,7J) 00006400
IF (RSUM.GI.20.DO) GO TO 1200 00006500
CALL ALLFOL(X,n,XEEEEEE,JU,JD,VV,DV,DD,RP) 00006600
CALL APAHEP(XEE,BEE,JU,JD,VV,DV,DDV) 00006700
IF (DNN.GT.10.CO) GO TO 900 00006800
CALL DIPOLE(X,R,J[I,JD,VV,DV,DDV) 00006900
CALL QIADIP (X,R,J,JD, VV,D,DDV,R2) 00007000
IF (DNN.GT.7.DO) GO TO 900 00007100
CALL QUACUA(X,R,JU,JD, VV,DV,DDV,R2) 00007200
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
DADV = DADV+DAN(K,IE,JU,JD)*DV(K,IE) 00008300
1100 LAN (K,IZ,JU,JD) = EX*(DAN(K,IE,JU, JD).5DO*AN (JU,JD) *DV(K,IE)) 00C08400
PDAN(JU,JD) = EX*(DDAN(J.U,JD)DADV+AN(JU,JD)*(.25DO*DVDV.5CO*DDV)00008500
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.E34) A = 1.E34 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
SUBROUTINE ALLPOL(X,R,XEE,REE,JU,JD,U,DU,DDU,RP) 00000100
C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
C Input values are the electron positions and radial distances (X, R,
C XEE, and REE), RP used in Eq. (49), and Jo and JD which label the
C permutation from Eq. (41). The function U from Eq. (41) and its
C derivatives DU and DDU are initialized to zero then contributions
C from Eq. (49) are added. The parameter EPS, coefficients ALC, and
C knots ALK are from Table 5.
C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
IMPLICIT REAL*8(AH,OZ) 00C00200
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
COMMON/ALP/ALC(6),ALK(6),EPS,DNN 00000500
COMMCN/MOMN/MW (2,2) 00000600
MW(1,1) = JU 00000700
M (1,2) = JD 00000800
MW (2,1) = 3JU 00000900
MW(2,2) = 7JD 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=3N 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)**3XEE(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) **2lDO) /RBPE**33. DO* ((1((2,I) /Rp (2, 1)) **00o04loo
x 2 1 .DO) /RP (2, I)**3+ ( (B (1 ,J)/RP (1,J) ) **21. 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
CDU=DDU+ DDFF*ALL+2.00*DADF+FF*DDA 00004400
400 CONTINUE 00004500 00
BETURN 00004600 D
END 000047C0
SUBROUTINE QUADIP(X,R,JU,JD,U,DDDU,R2) 00000100
C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
C Dipolequadrupole effects, Eq.(412), are added to U,DU, and DDU,
C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
IMPLICIT RlEAL*8(AH.OZ) 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 = CCZZ 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 DQX(3,2,J)*CCX(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)+2DO*(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)+(CCP3.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
00
SUBBOUTINE CUAQUA(X,R,JU,JD,UU,DDU,E2) 00000100
C Quadrupcleguadrupole effects, Eq. (413), are added to U, DU, and DDU.
C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
IMPLICIT REAI*8(AH,0Z) 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*(CCRSUNM+5.DO*ZDIF2) +3.DO*ER*RR15.DO*(X(3,1,I)**2* 00002400
X a2(2,)+X(3,2,J)**2*2(1,I))+ZZ*(30.O*RSUM70.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*ZDIF26.DO* 00003100
X RSUM) +X(K,1,)*(6. DO*R2(2,J)+60.DO*ZZ30.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*ZDIF26.DO* 00003300
X RSUM) X(K,2,J)*(6.DO*R2(1,I)+60. DO*ZZ30.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 150DO*X(3,1,I)**2)70.DO*X(3,2,J)**324.DO*X(3,1,I)*R2 (2,J) 00004100
DQ(3,J)=48.DO*CC*ZDIF+X(3,1,I)*(30.DO+ZDIF2+24.DO*RSU.180.EO*ZZ00004200
S 150.DO*X(3,2,J)**2)70.DO*X(3,1,I)**324.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
SUBRCUTINE APAREP(XEE,REE,JU,JD,U,DU,DDU) 00000100
C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
C The direct electronelectron interactions, Eq.(417), are added to
C U, DU, and DDIJ.
C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
IMPLICIT REAL*8(AH,0Z) 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
CALL SPLINE (2,EEC,EEK,REE(I,J),F,DF,DDF,1) 00001100
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
SUBROUTINE POLY(X,R,XEE,REE,NNT,JU,JD,AJAP,DJAP,DDJAP) 00000100
C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
C Values for Schwartz's stu 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. (42) are determined ty
C NS, NT, and NU; the coefficients for these terms are given by CP.
C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
IMPLICIT REAL*8(AH,OZ) 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,3JU)+R(2,7JD) 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 = 3JI 00002100
JD = 7JD 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.D25) T =1.D25 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(I1)*SRS 00004900
DO 400 T = 2,9 00005000
400 UP(I) = UP(I1)*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*(RUPIRDNI) 00005700
D3 = SIOI* (2IT1+ZIT2) 00005800
D4 = TI*I*(2 I1ZIT2) 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* (PI1.DO) *SI2+DFLOAT (JP* (JP1) ) 00006700
X TI2+DFLOAT (KP* (KP1)) *UI2+I*D1+DFLOAT (JP)*D2+2D00*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.E35) B(NN) = 1.E35 00007400
DO 600 KX = 1,3 00007500
DP(KX,JU) = ((DS(NN) DT(NN))*X(KX,NN,JU)*RUPIDU(NN)**EE(KX,JU,JD)000076C0
X )/B(NN) 00007700
DP(KX;JD) = ((DS(NN)DT(NN)) *X(KX,IN,JD)*RDNI+DU(NN)*XEE(KX,JU,JD)00007800
X )/B(NN) 00007900
DUDP = DUDP+CP (KX,JO)*X(KX,NN,JU)*RUPI+DP(KX,JD)*X(KX,NN,JD)*RDNI 00008000
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
DDJAP = AJAP*(DDP+AK*(DUDP.5DO*DDU+.25DO*AK*DUDU)) 00008800
DO 800 NN = 1,2 00C08900
JU = 3JU 00009000
JD = 7JD 00009100
DO 800 KX = 1,3 00009200
DJAP(KX,JD) = AJAP*(.5DO*X(KX,NN,JD)/E(NN,JD)*AK+DP(KX,JD)) 000CC930
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
SUBROUTINE SPLINE(NK,C,AK,X,F,DE,DDF,NS) 00000100
C <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< <<<<<<<<<<<<<< <<< <<<<< <<<<<<
C Spline functions (F) for Eq. (414) and its derivatives (DF and EEF)
C are calculated. The coefficients and knot positions are given by
C C and AK.
C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
IMPLICIT REAL*8(AH,OZ) 00000200
DIMENSION C(NK),AK(NK) 00000300
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 = XAK(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
SUBROUTINE SETUP 00000100
C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
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.
C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
IMPLICIT REAL*8 (AH,OZ) 00000200
COMMON/LLL/DX,DZ,VOL,,RM,REM,HIEC, HCON, CON,DUN 00000300
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(I1) 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=QC(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 = 11 00003600
X(I)=.15DO*DFLOAT(I1) 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
SUBROUTINE CONFIG(X,WTAV,NNT) 00000100
C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
C Electron positions (X), oneelectron weight functions (H), and the
C final averaged weight function from Eq. (232) (1/WTAV) are determined.
C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
IMPLICIT REAL*8 (AH,OZ) 00000200
REAL*4 RNDMF 00000300
COMMON/LLL/DX,EZ,VOL,RM,REM,HEEC,HMCON,HCON,DNN 00000400
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(I1) 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(I1) 00002200
BEE = REMH*NDIF(1.) 00002300
AC = 2.D*O+NDMF(1.)1.D0 00002400
AS = DSQRT (1.DOAC*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*(DNNDZ) 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.DOAC*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

