High Precision Hy-CI Variational Calculations

for the Ground State of Neutral Helium and

Heliumlike ions

James S. Sims

National Institute of Standards and Technology

Gaithersburg, MD 20899-8911

Stanley A. Hagstrom

Department of Chemistry

Indiana University

Bloomington, Indiana 47405

Abstract

Hylleraas-Configuration Interaction (Hy-CI) method variational calculations

with up to 4648 expansion terms are reported for the ground 1S state of neutral

helium. Convergence arguments are presented to obtain estimates for the exact

nonrelativistic energy of this state. The nonrelativistic energy is calculated to be

-2.9037 2437 7034 1195 9829 99 a.u. Comparisons with other calculations and

an energy extrapolation give an estimated nonrelativistic energy of -2.9037 2437

7034 1195 9830(2) a.u., which agrees well with the best previous variational energy, -2.9037 2437 7034 1195 9829 55 a.u., of Korobov[1], obtained using the

universal (exponential) variational expansion method with complex exponents[2].

In addition to He, results are also included for the ground 1S states of H

−, Li+,

Be++, and B

+3

.

1
I. INTRODUCTION

In a review article on Computational Chemistry in 1996, Clementi and Corongiu [3]

stated that using an Hy-CI [4, 5] expansion to solve the dynamical correlation is nearly

impossible for more than 3 or 4 electrons. While that may have been true in 1996, its

validity today is being challenged by the availibility of cheap CPUs which can be connected in parallel to significantly (orders of magnitude) enhance both the CPU power

and the memory that can be brought to bear on the computational task. In this paper we

address some of the issues arising from an attempt to come up with a good technique

for obtaining very precise energies for few electron atomic systems using the Hy-CI

formalism. Not only do we address the issue of choice of terms in the wave function,

but we also address another fundamental mathematical issue arising in these types of

calculations. In any attempt to get very precise energies, large basis sets have to be

employed, which means that linear dependence in the basis set is never very far away.

To proceed to several thousand terms in a wave function, extended precision arithmetic

is needed to obviate the linear dependence problem, which in turn leads to high CPU

costs. This work employs a novel wave function, namely, a wave function consisting

of at most a single r12 raised to the first power combined with a conventional nonorthogonal CI basis. We believe that this technique can be extended to multielectron

systems [6], where the use of at most a single rij (to the first power) retains the power

of rij factors in the wave function without making the integral evaluation overly complicated. We use this technique to determine the nonrelativistic ground state energy of

neutral helium more precisely than the best previous calculation [1]

1

, we discuss how

we solved the extended precision problem, and we discuss the implications of this for

3 and 4 electron systems.

1Korobov[1] has done the best calculation to date on the ground state of neutral helium by considering

variational basis functions of the type exp(-αnr1 - βnr2 - γnr12) with complex exponents. See Thakkar and

Smith[7] and Frolov and Smith[2] for a general discussion of this method.

2
II. METHOD OF CALCULATION

Variational methods based on explicitly correlated wave functions (wave functions including rij terms) are known to give the most accurate upper bounds to energy states,

and the inclusion of rij terms in the wave function has become increasingly common,

at least for few electron atomic systems (N ≤ 4). For example, standard Hylleraas

(Hy) technique [8] calculations are perhaps best exemplified by the work of Drake and

collaborators [9, 10, 11] and employ factors of powers of rij in the wave function.

However, for four-electron atomic systems, there are already significant unresolved

integration problems when Hylleraas basis sets are employed conventionally [12], in

what may be referred to as an Hy-rij technique. In contrast, the Hy-CI method developed by us [4] and also independently by Woznicki  ́ [13] does not suffer from this

restriction. Woznicki  ́ and coworkers refer to the method as SCC (Superposition of Correlated Configurations) and have employed it for an accurate Li ground state[14], and

the approach has also been used to obtain accurate excited S states of Li [15]. For our

part, we have used this method to calculate energies, oscillator strengths, and polarizabilities of two-, three-, and four-electron systems (see [16] and references therein). In

these approaches, which we refer to as Hy-CI, the wave function is expanded as a linear

combination of configurations, each of which contains at most one rij to some power.

Using this type of wave function, the most cumbersome integrals that are needed for

atoms with an arbitrary number of electrons are dealt with already in the four-electron

problem. This point was noted by Clary and Handy [17] who demonstrated the utility of the Hy-CI method for N-electron atomic systems up to and including N = 10

by carrying out a calculation on atomic Ne. Clary [18, 19] and Clary and Handy [20]

also demonstated the utility of the technique for systems other than atoms by performing Hy-CI calculations on one-positron atomic systems and on many-electron diatomic

molecules. Nevertheless, in spite of all of the progress that has been made with Hy-rij

and Hy-CI calculations in recent years, there still are practical issues to be resolved.

Conventional Hy-rij calculations have proven to lead to better convergence estimates

than Hy-CI calculations, which for 3 and especially 4 or more electrons can best be

described as requiring educated guesswork to select the terms in the wave function,

making it very hard to obtain good estimates of upper bounds to the exact nonrelativistic energies. In this work we show how an Hy-CI calculation can do as well as

an Hy-rij calculation (indeed, we improve the Hy-rij energy estimate by two decimal

places), while retaining the ability to give good convergence estimates for the upper

bound to the exact nonrelativistic energy. In a future work we hope to make similar

progress for three and four-electron systems [6].

We also note that the best previous calculation to date on the He ground state is neither

Hy-rij nor Hy-CI, but one that employs a variational expansion in products of exponentials in the problem interparticle distances[1, 7, 2]. We made no real attempt in this

3
work to minimize the number of expansion terms (one would certainly not have this

luxury when dealing with larger atoms), resulting in an expansion that is approximately

twice as long as Korobov’s[1]. While the exponential products approach provides almost as good an energy and a more compact wave function for the He ground state

and is straightforward to implement, it is not clear how useful it will be for atoms with

more than two electrons.

4
III. WAVE FUNCTIONS

The Schrodinger  ̈ equation for a three-body system consisting of a nucleus of Ze, mass

M, and two electrons of charge -e and mass me has been treated thoroughly by Drake

[21]. In this work, we focus on the nonrelativistic energy, using a hamiltonian which,

for N electrons, is (in atomic units 2)

HNR =

X

N

i=1

(−

1

2

∇2

i −

Z

ri

) +

X

i<j

1

rij

(1)

Note that this hamiltonian does not include mass polarization, which is commonly

treated using first order perturbation theory [22].

For two electrons, the time-independent, nonrelativistic Schrodinger  ̈ equation becomes

HNRΨ(r1, r2) = EΨ(r1, r2) (2)

Since the Schrodinger  ̈ equation is not separable in the electron coordinates, basis sets

which incorporate the r12 interelectronic coordinate are most efficient. The method

we use to incorporate interelectronic coordinates in the calculations is our combined

Hylleraas - Configuration Interaction (Hy-CI) method [4]. The Hy-CI method wave

functions can be written in a form which is slightly different from conventional twoelectron Hy-rij calculations [23, 24, 25, 26]

Ψ(r1, r2) =

X

N

K=1

CKΦK(r1, r2), (3)

where the terms ΦK are specifically of the form

ΦK(r1, r2) = O(L

2

)(1 ± P12)r

i

1

r

j

2

r

k

12Y

m

l

(1)Y

m0

l

0 (2)e

−ξr1−ηr2

. (4)

P12 is the operator which permutes electrons 1 and 2, and the plus sign is for singlet

levels, the minus sign for triplet levels. In this work, the r12 power k is restricted to

be either 0 or 1, i.e., in contrast to conventional Hy-CI calculations the r12 power is

limited to the first power (a novel feature). Y

m

l

is a normalized spherical harmonic in

the Condon and Shortley phase condition [27].

Note that this form of the wave function differs from conventional two-electron Hy-rij

calculations in that k is ≤ 1 and instead of higher powers of rij there are higher spherical harmonics in the basis set, with the O(L

2

) idempotent orbital angular momentum

2The atomic unit of energy is chosen as

μe4

h ̄2 = 1 a.u. (of energy), where μ = memN /(me + mN )

5
projection operator [28] projecting out the proper symmetry. For two electrons, and

1S,3S(L = 0) symmetry, l0 = l. Furthermore, to within some constant factor,

O(L

2

)Y

0

l

(1)Y

0

l

(2) = Pl(cosθ12) (5)

where Pl(cosθ12) is a Legendre polynomial [5]. Therefore the wave function we actually used has

ΦK(r1, r2) = (1 ± P12)r

i

1

r

j

2

r

k

12Pl(cosθ12)e

−ξr1−ηr2

. (6)

Note that the Legendre polynomials in Eq. (6) can be expressed in terms of spherical

harmonics of (θ1, φ1) and (θ2, φ2) by the spherical harmonic addition theorem:[27].

Pl(cosθ12) =

4π

(2n + 1)

X

l

m=−l

Y

−m

l

(1)Y

m

l

(2) (7)

Substituting this expression in Eq. (6) and taking k = 0 gives a conventional CI configuration for S states. There is also a close relationship between our (Hy-CI) wave

function and Hy-rij wave functions. Since

r

2

12 = r

2

1 + r

2

2 − 2r1r2P1(cosθ12), (8)

r

k

12 for k>1, when expanded out, gives a polynomial in powers of Pl(cosθ12) times 1

or r12 depending on whether k is even or odd. The powers of Pl(cosθ12) can in turn be

expressed in terms of the Pl(cosθ12) polynomials, thus making explicit the connection

between Hy-rij and Hy-CI. Although other (radial) factors will be different in the two

methods, they are essentially equivalent methods, a fact we have used in coming up

with our final wave function.

Drake[29] has pointed out the need for “doubling” basis sets so there is a natural partition of the basis set into two distinct distance scales - one appropriate to the complex

correlated motion near the nucleus, and the other appropriate further out. Drake uses

just two sets of orbitals to accelerate convergence; we use two sets of basis functions

for each l, differing in the nonlinear parameters ξ and η. The first set has an orbital

with an orbital exponent ξ = η that makes it essentially a valence shell orbital. For the

second set ξ = η again, and the orbital has a large exponent. Valence (outer) orbitals

are optimized whereas the inner orbital set is optimized for s− and p− but not beyond since the contributions for l > 1 are not large and the energy is a slowly varying

function of the orbital exponents anyway.

To clarify this further, consider Table 1, where we list the basis set for our final 4648

term wave function. Each line in the table specifies both the {inner,inner} and {outer,outer}

set for each l quantum number ( = Lmax). In the table, only the minimum information

needed to specify the basis set is listed, namely, an exponent ζ (ξ = η = ζ), the l

quantum number for orbitals with that exponent, and the number of orbitals with that

exponent (i and j in Eq. (4) range from l to norbs + l - 1). This gives rise to norbs (norbs

+ 1) / 2 terms without r12, and the same number with r12. There are 380 terms for l=0

6
Table 1: Hy-CI results for 1

1S He

Lmax norbs ζ norbs’ ζ’ N Energy(E) in a.u.

0 19 2.20 19 25.0 760 -2.9034 9832 0585 8801 1976 93

1 19 3.05 19 40.5 1520 -2.9037 2426 8354 6045 9864 64

2 19 3.50 19 40.5 2280 -2.9037 2437 6954 6918 8919 88

3 18 3.90 18 40.5 2964 -2.9037 2437 7034 0541 4170 45

4 15 4.50 15 40.5 3444 -2.9037 2437 7034 1195 3899 21

5 14 5.20 14 40.5 3864 -2.9037 2437 7034 1195 9822 02

6 14 6.00 14 40.5 4284 -2.9037 2437 7034 1195 9829 62

7 13 6.50 13 40.5 4648 -2.9037 2437 7034 1195 9829 99

and ζ = 2.2, an {s

2}r12 set. There are also 380 terms for l=0 and ζ = 25.0, an {s’2}r12

set. Since the orbital exponent is the same for all terms in a set, we refer to terms of the

second one as {inner,inner} (large orbital exponent brings it in), and those of the first

set as {outer,outer}.

Table 1 also lists our final energy value for the He ground state. In the table we give

the energy as we progressively add in higher angular momentum basis sets. These are

what might be referred to as s-wave, p-wave, ..., out to and including j-wave sets times

both 1 and r12. In the table, Lmax = 0 refers to the l=0 {inner,inner} set plus the l=0

{outer,outer} set for a total of 380 + 380 = 760 terms. Similarly for succeeding rows

in the table ( Lmax = 1 includes l=0 and l=1, etc). In all of the tables, we use N to refer

to the number of terms in the wave function (see Eq. (3)).

In our final wave function, only {inner,inner} and {outer,outer} configurations were

included. {inner,outer} were tried but were found to be unimportant in addition to

greatly increasing the expansion size and degree of linear dependence in the configuration basis.

We have found the surprising result that one simply can’t get to our answer without the

{inner,inner} terms without a really massive increase in the expansion length. Many

authors, starting perhaps with Roothaan and Weiss[30], have emphasized that the wave

function should have a cusp-like behavior at r12 = 0 such that

(

1

Ψ

∂Ψ

∂r12

)r12=0 =

1

2

(9)

and attributed the slow CI convergence to the necessity of piling up higher s, p, d, ...

angular terms in attempting to represent the details of this behavior. No one has emphasized the need for {inner,inner} terms at higher l values. Apparently, the need for

these {inner,inner} terms at higher angular terms has not been evident in conventional

CI up to now, because conventional CI doesn’t describe the valence shell electron cusp

properly either and anything put in near the nucleus will try to contribute to the valence

7
shell cusp description (even an orbital with a large exponent will penetrate the valence

region with a large enough power of r). In any event, including both {inner,inner}

and {outer,outer} terms are key to the success of the current Hy-CI formalism. For

example, if one just drops the {inner,inner} terms from the expansion in the last line

of Table 1 and keeps all orbital exponents fixed, the result is E(2324 terms) = -2.9037

2437 7026 a.u., which is not as good as the Frankowski and Pekeris[31, 32] results

obtained in 1966.

8
IV. NONRELATIVISTIC ENERGIES -

RESULTS AND DISCUSSION

The wave function given by Eq. (3), Ψ(r1, r2), is a linear combination of configurations ΦK, where the coefficients CK are those which minimize the total energy, E,

given by

E =

< ΨHΨ >

< ΨΨ >

=

P

P

KL CKCLHKL

KL CKCLSKL

, (10)

where

HKL =< ΦK HΦL >; SKL =< ΦK ΦL > . (11)

The condition for the energy to be an extremum, δE = 0, is the well-known matrix

eigenvalue (secular) equation:

X

L

HKLCL =

X

L

SKLCL. (12)

Solving this equation is equivalent to solving the N-dimensional generalized eigenvalue

problem

HC = λSC (13)

where H and S have matrix elements HKL and SKL given by Eq. (11). We solve this

secular equation using both sequential and parallel inverse iteration algorithms, and

using real*16 and real*24 arithmetic, as discussed below.

We started out writing our own Microsoft3 ASM (MASM) quadruple precision (real*16,

∼ 32 digits) package simply because no Fortran 90 package for the PC had real*16 as

a native data type. So much of our early experience was with real*16 arithmetic. Then,

because of linear dependence problems (even at 32 decimal places), we extended the

real*16 package to real*24 extended precision (∼ 48 digits). Fortran 90, of course, has

user defined datatypes and operator extensions which make implementation of real*16

and real*24 floating point operations essentially automatic once the appropriate interface MODULEs are constructed (i.e., no need to explicitly call subroutines to carry

out the floating point operations). This is an important, very practical feature of the

calculations.

Except for some of the shorter expansions reported in Table 2, all results reported in

this paper were obtained using real*24 extended precision (∼ 48 digits) floating point

subroutines written in MASM including add, sub, mpy, div, exp, and log. There is a full

3The identification of any commercial product or trade name does not imply endorsement or recommendation by either the National Institute of Standards and Technology or Indiana University.

9
set of conversion routines between integers and single, double, and quad reals and our

extended precision format, as well as output routines in the extended precision format.

The parallel calculationsreported here were carried out at the National Institute of Standards and Technology on the instaNT.nist.gov NT Cluster, a 16 node cluster of Pentium

II systems running Microsoft Windows NT [Server 4.0] with 100-Mbit Fast Ethernet

for interconnection. The front-end machine is a 400MHz Pentium II with 780MB of

RAM and 19GB of disk storage. The back-end nodes are single-cpu 400MHz Pentium

IIs with 512MB of RAM and 8.4GB of local disk storage.

In the subsections that follow, we have endeavored to include all the data needed to

reproduce the results reported in Tables 1, 2, and 5, should future workers want to do

so.

10
A. Hy-rij results and comparison with earlier results

For a more detailed comparison with conventional Hy-rij methodology, we also undertook our own study of the Hy-rij method as applied to helium. Our results are

summarized in Table 2.

n

max

1

,nmax

2

, and n

max

12 are the maximum powers of r1, r2, and r12, respectively. N, the

number of terms, can best be explained with an example. For row 1 in Table 2, what we

did was to form all configurations from 4 outer s-orbitals with orbital exponent 2.0 and

then split the orbital exponent to 2.0, 2.2. Otherwise, one gets linear dependence very

quickly since 1s(2.0)1s(2.0), 1s(2.0)1s(2.2), and 1s(2.2)1s(2.2) are all very similar. The

way we did it, only 1s(2.0)1s(2.2) arises. And similarly for the {inner,inner} terms.

Hence one does not have 4 x 4 x 5 = 80 {outer, outer} terms but only [4(4+1)/2] x 5 =

50 terms, with a similar number of {inner, inner} terms, giving 100 terms in all.

These are real*16 (128 bit) calculations except for the last 5, which are real*24 (192

bit) calculations (and they are only a partial listing of our results). Our 2184 term result

agrees quite well with the 2114 term result of Drake[29] in Table 3, the best previous

Hy-rij calculation. As problem size increased, run times became a problem, especially

for the generalized eigenvalue problem step, so we decided to develop a portable parallel inverse iteration solver. Since the inverse iteration solver matrix representation

is a blocked one, we modified the secular equation step to generate H and S in the

appropriate block order. Then MPI[33] was used to run the same program on multiple

processors (on the same or different hosts) and give each host a block of the H (and

S) matrix, with no need to redistribute the matrices for the inverse iteration step. The

MPI code uses blocking sends and receives to do the equivalent of an MPI reduce and

then an MPI gather (on to the root process) in real*24 arithmetic. We have this parallel

inverse iteration solver running on the NIST NT cluster instaNT.nist.gov.

On a 3052 term calculation, run times were 41.51 seconds to construct the secular

equation and 6425.73 seconds to solve for one root using inverse interation. On 16

processors, solving the secular equation took only 555.66 seconds, i.e, a two hour

run on one processor became 10 minutes running on the cluster. We found that the

processing speed could be predicted, as a function of cluster size, by the simple scaling

law T = constant (s + (1 - s) /Nproc), where T is the runtime in seconds, constant

= 6419 in this case, and s is the inherently sequential part of the calculation. This

function is plotted in Fig. 1. We find that the sequential fraction s = 0.022, indicating

that the scaling is excellent and we could go to a larger number of processors (if we had

them). Finally, one last reason for going parallel is to circumvent the executable image

size limitation for sequential programs under both Windows and NT DOS modes. To

go beyond about 3300 terms in our calculations involves either writing intermediate

results out to disk or spreading the calculation (and memory requirements) across a

11
Table 2: Hy-rij results for 1

1S He (including orbital basis)

n

max

1 n

max

2 {outer,outer} {inner,inner} n

max

12 N Energy(a.u)

4 4 {2.0,2.2} {6.0,6.2} 4 100 -2.9037 2280 7177 558

6 6 {2.0,2.2} {6.0,6.5} 4 210 -2.9037 2437 3077 229

6 6 {2.0,2.2} {6.0,6.5} 6 294 -2.9037 2437 6814 373

6 6 {2.0,2.2} {6.0,6.5} 7 336 -2.9037 2437 6900 260

6 6 {2.0,2.2} {6.0,6.5} 8 378 -2.9037 2437 6947 786

7 7 {2.0,2.2} {6.0,6.5} 8 504 -2.9037 2437 7033 120

8 8 {2.0,2.2} {6.0,6.5} 8 648 -2.9037 2437 7033 898

8 8 {1.85,2.05} {8.6,8.8} 9 720 -2.9037 2437 7034 0967 2287

9 9 {1.95,2.15} {9.8,10.0} 9 900 -2.9037 2437 7034 1150 3

10 10 {1.95,2.15} {11.0,11,2} 9 1100 -2.9037 2437 7034 1184 93

10 10 {1.95,2.15} {11.0,11,2} 10 1210 -2.9037 2437 7034 1189 65

11 11 {2.10,2.30} {12.4,12.6} 10 1452 -2.9037 2437 7034 1194 42

11 11 {2.10,2.30} {12.4,12.6} 11 1584 -2.9037 2437 7034 1194 90

12 12 {2.20,2.40} {16.0,16.2} 11 1872 -2.9037 2437 7034 1195 8414

13 13 {2.20,2.40} {16.0,16.2} 11 2184 -2.9037 2437 7034 1195 9090

13 13 {2.20,2.40} {19.5,19.7} 12 2366 -2.9037 2437 7034 1195 9567

14 14 {2.20,2.40} {21.5,21.7} 12 2730 -2.9037 2437 7034 1195 9725 15

15 15 {2.20,2.40} {22.0,22.2} 12 3120 -2.9037 2437 7034 1195 9793 92

16 16 {2.20,2.40} {23.5,23.7} 12 3536 -2.9037 2437 7034 1195 9815 09

number of processors.

12
B. Extrapolation of Hy-CI results

The Hy-rij results gave us the insight we needed to do the actual Hy-CI runs discussed

previously. Table 3 compares our final energy value for the He ground state with previous results. The calculated value is a rigorous upper bound to the exact nonrelativistic

energy of this state. In Table 3, we have also included our extrapolation of the energies

given in Table 1, obtained using a formula from Pekeris and co-workers[24]

Eextrapolated = E1 +

(E1 − E0)(E2 − E1)

2E1 − E0 − E2

, (14)

where E0, E1, and E2 are the values at Lmax = 5, 6, and 7 respectively. This formula gives -2.9037 2437 7034 1195 9830 01 a.u for the extrapolated energy. Based

on the fact that extrapolations tend to undershoot the exact result, the sensitivity of E

to the non-linear parameters and the convergence in L, we estimate the exact nonrelativistic energy to be 2.9037 2437 7034 1195 9830 a.u. with an uncertainty of 1 or 2

in the last digit. This result is in complete agreement with the landmark calculations

of Korobov[1]. In addition to values for He, we have also included results for other

members of the He isoelectronic sequence, namely, H

−, Li+, Be++ and B+3

. For

these results, no exponent minimization was done, just a simple nuclear charge based

scaling of the wave function of Table 1. For example, the results for Li+ were obtained

by scaling the orbital exponents by a factor of 3.0 / 2.0 and for Be++ the orbital exponents were scaled by a factor of 4.0 / 2.0. For Li+ we varied the orbital exponents

extensively using shorter expansions and the result was that orbital exponent variations

are important in the 20th decimal place. Our results for H

−, Li+, and Be++ agree with

the results of Drake[21]

4

to the 16 decimal places he reports (he doesn’t have a value

for B

+3) and, except for H− 5

, are more precise by at least three decimal places. These

results were easy to obtain, which points out how flexible our expansion is.

We would like to point out again the close relationship between Hy-rij and Hy-CI. In

Table 4 we show this by tabulating El results for l = 0,1,2,... levels of truncation of the

wave function expansion used in Table 1 versus the equivalent Hy-rij truncation results

for the best wave function in Table 2

6

. The two energy columns are clearly related, the

minor differences being easily explained by the different functional forms, methods of

optimization of non-linear parameters, and the fact that these limits are not rigorously

defined (except for the first one).

4Drake’s value for Be++ contains a misprint. -13.1... should be -13.6...

5Convergence is slow for H−

6There is no entry on the Hy-rij side corresponding to l = 6 since 12 is the upper limit to nmax

12 in our

code.

13
Table 3: Comparison with previous explicitly correlated calculations for 1

1S He and

He-like ions

Author N Energy(a.u.)

Frankowski and Pekeris[31, 32] 246 -2.9037 2437 7032 6

Freund et.al.[34] 230 -2.9037 2437 7034 0

Thakkar and Koga[35] 308 -2.9037 2437 7034 1144

Baker et.al.[36] 476 -2.9037 2437 7034 1184

Goldman[37] 8066 -2.9037 2437 7034 1195 9382

Bur ̈ gers et.al.[38] 24497 -2.9037 2437 7034 1195 89

Drake[29] 2114 -2.9037 2437 7034 1195 9582

Korobov[1] 2200 -2.9037 2437 7034 1195 9829 55

Korobov[1] Extrapolated -2.9037 2437 7034 1195 9830 6(10)

This work 4648 -2.9037 2437 7034 1195 9829 99

Extrapolated ∞ -2.9037 2437 7034 1195 9830 01

Estimated -2.9037 2437 7034 1195 9830(2)

This work, H

− 4640 -0.5277 5101 6544 3771 2249

This work, Li+ 4284 -7.2799 1341 2669 3059 6489

This work, Be++ 4648 -13.6555 6623 8423 5867 0206

This work, B

+3 4648 -22.0309 7158 0242 7815 4163

Table 4: Comparison of Hy-rij and Hy-CI results for 1

1S He

l El(Hy-CI) in a.u. n

max

12 E(Hy-rij ) in a.u.

0 -2.9034 9832 0585 1 -2.9034 9825 7090

1 -2.9037 2426 8354 3 -2.9037 2426 8491

2 -2.9037 2437 6954 6918 5 -2.9037 2437 6955 3843

3 -2.9037 2437 7034 0541 7 -2.9037 2437 7034 0553

4 -2.9037 2437 7034 1195 3899 9 -2.9037 2437 7034 1195 4000

5 -2.9037 2437 7034 1195 9822 11 -2.9037 2437 7034 1195 9796

6 -2.9037 2437 7034 1195 9830 13

14
C. Configuration Interaction

As a test of our integral and eigenvalue routines and also to investigate the basis set linear dependence problem, we have carried out a number of conventional configuration

interaction calculations on He. In this connection, in Table 5 we report the results of a

4699 term expansion carried out using real*24 arithmetic. A typical symmetry adapted

configuration has the form

(1 + P12)r

i

1

e

−ξr1 r

j

2

e

−ηr2 Pl(cosθ12) (15)

where ξ = η for all configurations except the 49 ss terms in the third row. Nonlinear parameters were carefully chosen for terms up through ff. For this part of the

calculation there were considerable problems with linear dependence in real*16 mode,

forcing us to move to real*24 arithmetic. For higher l values we used an approximate

linear relationship between orbital exponent and l-value, with occasional checks on the

accuracy of the approximation. In any case, the energy contributions for the higher

l-values are very insensitive to the orbital exponents and depend much more on the

number of basis orbitals (powers of r) for a given l-value, even as the contribution for

an l-value decreases to less than microhartree values (e.g., l=16-18). The upper limit

on l of 18 was set by program and operating system limitations of the PC used to do

these particular calculations.

Also included in Table 5 for comparison purposes are the excellent recent results of

Jitrik and Bunge[39] for the He ground state. These results were obtained with much

smaller basis sets for each l-value but with essentially a separately optimized exponent

for each basis orbital. Our approach, however, has been to use a basis with as few

non-linear parameters as possible, that is, a basis very like that used in our Hy-rij and

Hy-CI calculations. Note that although one can get spectroscopic accuracy (1 cm−1

is approximately 5 microhartrees) for He using CI, microhartree accuracy is still not

available. We estimate that to get microhartree accuracy would require terms up to

l=30 and maybe up to 1200 more terms, although with a more efficient basis (more

non-linear parameters) the overall expansion length could be shortened considerably.

This and other options are currently being investigated[40]. Nevertheless, after some

70 years of trying, the best CI treatment for He is still accurate to only 5 decimal places!

15
Table 5: CI results for the 1S ground state of neutral He (a.u.)

l n ζ N El(a.u.) Previous work[39]

ss 0 19 4.1 190

ss 0 21 20.3 421

ss 0 7 4.1, 20.3 470 -2.8790 2875 65 2.8790 2862 7

pp 1 17 4.7

pp 1 21 20.9 854 -2.9005 1621 99 -2.9005 1577 4

dd 2 16 6.1

dd 2 21 22.5 1221 -2.9027 6680 53 -2.9027 6612 6

ff 3 16 7.0

ff 3 21 24.0 1588 -2.9033 2101 62 -2.9033 2011 0

gg 4 21 14.0 1819 -2.9035 1846 49 -2.9035 1734 1

hh 5 21 16.0 2050 -2.9036 0551 51 -2.9036 0419 6

ii 6 21 18.0 2281 -2.9036 4964 42 -2.9036 4818

jj 7 21 19.0 2512 -2.9036 7432 76 -2.9036 7269

kk 8 21 21.0 2743 -2.9036 8919 34 -2.9036 8741

ll 9 21 22.0 2974 -2.9036 9865 63 -2.9036 9675

mm 10 21 23.0 3184 -2.9037 0497 38 -2.9037 0296

nn 11 20 24.0 3415 -2.9037 0932 48 -2.9037 0720

oo 12 20 25.0 3625 -2.9037 1243 32 -2.9037 1019

pp 13 20 27.0 3835 -2.9037 1470 34 -2.9037 1237

qq 14 20 28.0 4045 -2.9037 1640 56 -2.9037 1396

rr 15 20 29.0 4255 -2.9037 1769 92 -2.9037 1515

ss 16 18 30.0 4426 -2.9037 1869 56

tt 17 17 31.0 4579 -2.9037 1948 04

uu 18 15 32.0 4699 -2.9037 2009 19

“Exact” -2.9037 2437

Error -0.0000 0428

16
V. CONCLUSIONS

The renewed interest in rij s is coming from the incredible accuracy of experiments

these days. To do as well theoretically requires explicitly correlated wave functions, in

general. Goldman’s [37] work is the exception for He, as he gets an energy of -2.9037

2437 7034 1195 938(50) a.u. (8066 term expansion). His estimate is in agreement with

our calculated result, -2.9037 2437 7034 1195 9829 99 a.u. However Goldman’s CI

is not the conventional kind of CI. Instead of powers of r1 and r2 he uses powers of

r< and r>, which is fine for He but leads to 3- and 4-electron integrals for Li and Be

(unlike conventional CI). Only Li has been done essentially as accurately as He, and

that only with a 6000 fold increase in CPU requirements7. And the integrals get messier

with conventional rij techniques. We have always thought that the idea behind the HyCI technique was to marry the relative ease of doing CI with the better convergence

obtained when rij factors are used. We have pointed out, for He, the close connection

between Hy-rij . It seems to us that the Hy-CI method selects the important term types

in a more natural manner. Also, the calculation, at least for He, is easier. If one restricts

the wave function to a single rij in each term (n≤1) then the most difficult integrals

are already dealt with at the N = 4 level and the calculation is greatly simplified. The

only problem has been figuring out a systematic technique for choosing basis orbitals

and configurations that leads to good energy estimates with a reasonable number of

expansion terms.

One of the results of this study has been a detailed comparison of Hy-rij and Hy-CI

techniques so that we are able to move on to calculations of Li and possibly Be atom

states to very good accuracy. As we move on to larger systems, computational times

will increase greatly, as well as memory requirements, so it is imperative to program

these routines using parallel programming techniques and MPI (for portability).

7Drake [21] says that only for He and Li have results of spectroscopic accuracy been obtained. He says

the demand on computer resources increases by a factor of about 6000 to reach spectroscopic accuracy for

Li over He.

17
Acknowledgements

We are deeply indebted to Carlos F. Bunge for providing us with the details of the Jitrik

and Bunge[39] He CI calculation, to Frederick W. King for bringing us up to date on rij

calculations on two, three, and four electron systems, to William Martin and Peter Mohr

for pointing out the Korobov work, and to Grzegorz Pestka for kindly informing us of

the activities of Professor W. Woznicki’  ́ s explicitly correlated wavefunction group. The

authors would also like to thank Charles Bouldin for Fig. 1, Barry J. Hershman and

Michael D. Indovina for putting together the NT cluster that was used in the parallel

runs reported in the paper, and Tom Kotowski of Microsoft for NT cluster support.

18
References

[1] V. I. Korobov. Phys. Rev. A, 61:064503, 2000.

[2] A. M. Frolov and Jr. V. H. Smith. J. Phys. B: At. Mol. Opt. Phys., 28:L449, 1995.

[3] E. Clementi and G. Corongiu. From Atoms to Large Chemical Systems with

Computational Chemistry. Revue de l’Institut Francais du Petrole, 51(1):19,

1996.

[4] J. S. Sims and S. A. Hagstrom. Phys. Rev. A, 4(3):908, 1971.

[5] J. S. Sims and S. A. Hagstrom. J. Chem. Phys., 55(10):4699, 1971.

[6] J. S. Sims and S. A. Hagstrom. in progress.

[7] A. J. Thakkar and Jr. V. H. Smith. Phys. Rev. A, 15(1):1, 1977.

[8] E. A. Hylleraas. Zeit. Phys., 54:347, 1929.

[9] G. F. W. Drake. In F. S. Levin and D. A. Micha, editors, Long-range Casimir

Forces: Theory and Recent Experiments on Atomic Systems, page 107. Plenum,

New York, 1993.

[10] G. F. W. Drake and Z.-C. Yan. Phys. Rev. A, 46:2378, 1992.

[11] G. F. W. Drake and Z.-C. Yan. Chem. Phys. Lett., 229:486, 1994.

[12] F. W. King. J. Chem. Phys., 99:3622, 1993.

[13] W. Woznicki.  ́ In A. Jucys, editor, Theory of electronic shells in atoms and

molecules, page 103. Mintis, Vilnius, 1971.

[14] G. Pestka and W. Woznicki.  ́ Chem. Phys. Lett., 255:281, 1996.

[15] J. Pipin and D. M. Bishop. Phys. Rev. A, 45(5):2736, 1992.

[16] J. S. Sims, D. R. Parmer, and J. M. Reese. J. Phys. B: Atom. Molec. Phys., 15:327,

1982.

[17] D. C. Clary and N. C. Handy. Phys. Rev. A, 14(5):1607, 1976.

[18] D. C. Clary. J. Phys. B: Atom. Molec. Phys., 9(17):3115, 1976.

[19] D. C. Clary. Mol. Phys., 34:793, 1977.

[20] D. C. Clary and N. C. Handy. Chem. Phys. Lett., 51:483, 1977.

[21] G. F. W. Drake. High Precision Calculations for Helium. In G. F. W. Drake,

editor, Atomic, Molecular, and Optical Physics Handbook, pages 154–171. API

press, New York, 1996.

19
[22] F. W. King. High-precision calculations for the ground and excited states of the

lithium atom. Advances in Atomic, Molecular, and Optical Physics, 40:57–112,

1999.

[23] C. L. Pekeris. Phys. Rev. A, 127:509, 1962.

[24] B. Schiff, H. Lifson, C. L. Pekeris, and P. Rabinowitz. Phys. Rev. A, 140:1104,

1965.

[25] Y. Accad, C. L. Pekeris, and B. Schiff. Phys. Rev. A, 4:516, 1971.

[26] M. T. Anderson and F. Weinhold. Phys. Rev. A, 9:118, 1974.

[27] E. U. Condon and G. H. Shortley. The Theory of Atomic Spectra. Cambridge U.

P., Cambridge, England, 1963.

[28] P. O. Lowdin.  ̈ Rev. Mod. Phys., 36:966, 1964.

[29] G. F. W. Drake. Physica Scripta, T83:83, 1999.

[30] C. C. Roothaan and A. W. Weiss. Revs. Modern Phys., 32:194, 1960.

[31] K. Frankowski and C. L. Pekeris. Phys. Rev., 146:46, 1966.

[32] K. Frankowski and C. L. Pekeris. Phys. Rev., 150:366(E), 1966.

[33] Message Passing Interface Forum. MPI: A message-passing interface standard.

Int. J. Supercomput. Appl., 8(3/4):159–416, 1994.

[34] D. E. Freund, B. D. Huxtable, and J. D. Morgan, III. Phys. Rev. A, 29:1984, 1984.

[35] A. J. Thakkar and T. Koga. Phys. Rev. A, 50:854, 1994.

[36] J. D. Baker, D. E. Freund, R. N. Hill, and J. D. Morgan, III. Phys. Rev. A, 41:1247,

1990.

[37] S. P. Goldman. Phys. Rev. A, 57(2):R677, 1998.

[38] A. Bur ̈ gers, D. Wintgen, and J.-M. Rost. J. Phys. B: At. Mol. Opt. Phys., 28:3163,

1995.

[39] O. Jitrik and C. F. Bunge. Phys. Rev. A, 56(4):2614, 1997.

[40] J. S. Sims, S. A. Hagstrom, and C. F. Bunge. in progress.

20
Figure 1: He Hy-CI Scaling with Cluster Size 21