INSTITUTE OF PHYSICS PUBLISHING JOURNAL OF PHYSICS B: ATOMIC, MOLECULAR AND OPTICAL PHYSICS

J. Phys. B: At. Mol. Opt. Phys. 37 (2004) 1519–1540 PII: S0953-4075(04)74803-8

Mathematical and computational science issues in

high precision Hylleraas-configuration interaction

variational calculations: I. Three-electron integrals

James S Sims1 and Stanley A Hagstrom2

1 National Institute of Standards and Technology, Gaithersburg, MD 20899-8911, USA

2 Department of Chemistry, Indiana University, Bloomington, IN 47405, USA

E-mail: jim.sims@nist.gov

Received 5 January 2004

Published 17 March 2004

Online at stacks.iop.org/JPhysB/37/1519 (DOI: 10.1088/0953-4075/37/7/012)

Abstract

The most difficult integral arising in Hylleraas-configuration interaction

(Hy-CI) calculations, the three-electron triangle integral, is discussed. We focus

on recursive techniques at both the double precision and quadruple precision

level of accuracy while trying to minimize the use of higher precision arithmetic.

Also, we investigate the use of series acceleration to overcome problems of slow

convergence of certain integrals defined by infinite series. We find that a direct

+ tail Levin u-transformation convergence acceleration overcomes problems

that arise when using other convergence acceleration techniques, and is the

best method for overcoming the slow convergence of the triangle integral. The

question of calibrating an acceleration method is also discussed, as well as

ways to improve our work.

1. Introduction

Frolov and Bailey have recently published [1] a scheme for evaluating few-body auxiliary

functions and four-body integrals with extremely high accuracy using a sophisticated arbitrary

multiple precision (MP) Fortran 90 package (MPFUN) developed by Bailey [2]. Some of

these integrals also occur in attempting to obtain very precise energies for few electron

atomic systems using the Hylleraas-configuration interaction (Hy-CI) formalism. Our Hy-CI

employs a novel wavefunction, namely, a wavefunction with terms consisting of at most a

single rij raised to the first power combined with a conventional non-orthogonal configuration

interaction (CI) basis3. We have recently used this technique to determine nonrelativistic

ground-state energies of helium and helium-like ions [3, 4]. We used this technique 30 years

ago to determine energy levels of Li [5], Li− [6] and Be [7, 8] much less accurately with

3 Which eliminates some of the most difficult integrals Frolov and Bailey discuss.

0953-4075/04/071519+22$30.00 © 2004 IOP Publishing Ltd Printed in the UK 1519
1520 J S Sims and S A Hagstrom

the supercomputers of that era. Here we discuss how to efficiently evaluate the only difficult

integrals arising when using this technique in the more general case of N (number of electrons)

3. We will be focusing on recursive techniques (section 3) at both the double precision (DP)

and quadruple precision (QP)4 level of accuracy while trying to minimize the use of higher

precision arithmetic. In this sense the focus of this work is substantially different from that of

Frolov and Bailey. Also, we investigate the use of series acceleration to overcome problems

of slow convergence of certain integrals defined by infinite series (section 4).

Many of the details of the integrals which arise in an Hy-CI calculation of atomic systems

have been discussed previously [9]. In this paper we discuss changes that we have made to

our integral methodology in our most recent work [10]. As discussed in [9], all integrals

can be expressed in terms of the standard Condon and Shortley coefficients (cks) [11, 9] and

auxiliary functions A, V ,W, and X. In this paper we update the treatment of V and W auxiliary

functions to reflect changes we have made to most effectively use modern day computers to

increase the size (number of terms) and accuracy of the calculations. We reserve treatment

of X auxiliary functions to a future paper, but we note that our X functions are easier than the

more general functions treated by Frolov and Bailey, a consequence of the restriction of the rij

terms to a single rij raised to the first power (in contrast to Hy-rij calculations where products

of rij terms in wavefunctions have to be dealt with).

Our emphasis on efficient methods for calculating these integrals arises because of how

we handle the H matrix assembly problem, which is still the hardest part of a high precision

atomic energy level calculation, especially for Be and beyond. In our energy calculations, we

essentially multipass the assembly of Hamiltonian matrix elements, ending up after pass one

with long lists of integrals needed, which we then sort by integral type, orbital type, nonlinear

parameters, and so on so that integrals depending on the same types of auxiliary functions are

effectively grouped together. To make the integral calculation efficient we need fast, stable

methods for calculating large arrays of auxiliary functions from which large blocks (tens of

thousands in some cases) of integrals can be calculated at one time. Recursive techniques for

calculating these arrays are the only viable methods found for implementing this approach.

2. Hy-CI triangle integral evaluation

The three-electron rij integrals which arise in our Hy-CI calculations are of the form

I =



1(r1)2(r2)3(r3)R dr1 dr2 dr3, (1)

where R can be any of the following rij products,

r12r13, r12r−1

13 , r12r−113 r23,

and the  are ‘charge distributions’ made up of products of unnormalized Slater-type orbitals

(STOs) φ(r)

φi(r) = rni−1 e−αir

Y mi

li (θ , φ), (2)

i(r) = φ∗

i (r)φi(r). (3)

In equation (2), Y m

l (θ , φ)is a normalized spherical harmonic in the Condon and Shortley phase

convention [11]. Of these rij integrals, the one arising when R = r12r−1

13 r23, the so-called

‘triangle integral’, stands out because its Legendre expansion is an infinite series requiring

special treatment to achieve high precision without prohibitive computational costs. The

4 When we refer to double, quadruple, or other precision it is with respect to a 32-bit word.
Mathematical and computational science issues in high precision Hy-CI calculations 1521

remaining two rij integrals, the ones when R = r12r13 and R = r12r−1

13 , lead to short, finite

expansions whose evaluation has been well covered in the literature and will not be discussed

further here5. For a very thorough treatment of three-electron integrals over s-type STOs6 (see

[1–25] of the paper by Pelzl and King [12] and [16–40] in the review article by King [13]).

Yan and Drake [14] treat the more general case of three-electron integrals over nonspherically

symmetric STOs.

The triangle integrals which have to be evaluated in our Hy-CI calculations are therefore

I = φ1(r1)φ2(r2)φ3(r3)|r12r−1

13 r23|φ

1(r1)φ2(r2)φ3(r3)

=



1(r1)2(r2)3(r3)r12r−1

13 r23 dr1 dr2 dr3. (4)

Expanding the charge distributions i(r) as linear combinations of STOs fi(ri)

fi(r) = rNi−1 e−wir

Y Mi

Li (θ , φ) (5)

with Ni = ni + n

i − 1, wi = αi + αi, Li = |li − l



i|, |li − l



i| + 2,...,li + l



i − 2, li + l



i, Mi =

m

i − mi, I becomes

I = δ(M1 + M2 + M3, 0) 

L1,L2,L3



3

i=1



(2Li + 1)

1

2 cLi

(l

i, mi;li, mi)



J (L1, L2, L3), (6)

where the cLi are standard Condon and Shortley coefficients (the so-called cks) [11, 9]. The

factor J (L1, L2, L3) is given by

J (L1, L2, L3) = 1

(4π )3/2



r12r−1

13 r23 

3

i=1

{fi(ri) dri}

= 1

(4π )3/2



r

N1−1

1 r

N2−1

2 r

N3−1

3

r12r23

r13

exp(−w1r1 − w2r2 − w3r3)

× Y M1

L1 (1)Y M2L2 (2)Y M3L3 (3) dr1 dr2 dr3. (7)

Inserting the Legendre expansions for rij and r−1

ij we get

J (L1, L2, L3) = ∞

q=0

q

+L1

ns=|q−L1|

min(q+



L3,ns+L2)

nu=max (|q−L3|,|ns−L2|)

A(L1, L2, L3, ns,q,nu)

× 

1

l=0



1

n=0

Bns,1

l Bnu,1 n R(N1, N2, N3, N12(l), N13(q), N23(n)), (8)

with N12(l) = ns + 1 − 2l,N13(q) = q + 1, N23(n) = nu + 1 − 2n. The B coefficients are

Bm,1

0 = 1/(2m + 3) and Bm,11 = −1/(2m − 1), while

A(L1, L2, L3, ns,q,nu) = 1

{(2ns + 1)(2q + 1)(2nu + 1)}1/2



mq

cq (ns, mq + M1;L1, M1)

× cns (nu, mq + M1 + M2;L2, M2)cnu (q, mq ;L3, M3) (9)

5 We have discussed how we do these integrals in a previous paper [9]. 6 The STOs we use are defined fully in [9]. An s-type STO has l = 0, a p-type STO has l = 1, a d-type STO has

l = 2, etc.
1522 J S Sims and S A Hagstrom

is a three-electron angular factor with mq running from max[−q, −(ns +M1), −(nu+M1+M2)]

to min(q, ns − M1, nu − M1 − M2). R is a special case of equation (32) of [9] and is

R(N1, N2, N3, N12, N13, N23)

=



r

N1+1

1 r

N2+1

2 r

N3+1

3

r

N12+1

12< r

N13−1

13< r

N23+1

23<

r

N12

12>r

N13

13>r

N23

23>

exp(−w1r1 − w2r2 − w3r3) dr1 dr2 dr3

= W (N1 + N12 + N13 + 1,N2 − N12 + N23 + 2,N3 − N13 − N23 + 1;w1,w2, w3)

+W (N1 + N12 + N13 + 1,N3 − N13 + N23 + 2,N2 + 1 − N12 − N23;w1,w3,w2)

+ W (N2 + N12 + N23 + 3, N1 − N12 + N13, N3 − N13 − N23 + 1; w2, w1, w3)

+ W (N2 + N12 + N23 + 3, N3 − N23 + N13, N1 − N12 − N13 + 1; w2, w3, w1)

+W(N3 + N13 + N23 + 1,N1 + N12 − N13 + 2,N2 + 1 − N12 − N23;w3,w1, w2)

+W (N3 + N13 + N23 + 1,N2 + N12 − N23 + 2,N1 − N12 − N13 + 1;w3,w2,w1),

(10)

where N12 ≡ N12(l), N13 ≡ N13(q), N23 ≡ N23(n), and the auxiliary function W, whose

efficient evaluation is the main focus of this paper, is defined in section 3 by equation (15).

Expanding equation (8) and substituting for the B coefficients, and then inserting the

resulting expression for J (L1, L2, L3) into equation (6) leads to our final expression for I

I = δ(M1 + M2 + M3, 0)

∞

q=0

A(q), (11)

where

A(q) = 

L1,L2,L3



3

i=1

{(2Li + 1)

1/2

cLi(l

i, mi;li, mi)}

×

q

+L1

ns=|q−L1|

min (q+



L3,ns+L2)

nu=max (|q−L3|,|ns−L2|)

A(L1, L2, L3, ns,q,nu)

×



R(N1, N2, N3, ns + 1, q + 1, nu + 1)

(2ns + 3)(2nu + 3)

− R(N1, N2, N3, ns + 1, q + 1, nu − 1)

(2ns + 3)(2nu − 1)

− R(N1, N2, N3, ns − 1, q + 1, nu + 1)

(2ns − 1)(2nu + 3)

+

R(N1, N2, N3, ns − 1, q + 1, nu − 1)

(2ns − 1)(2nu − 1)



. (12)

If all orbitals are s-type, L1 = L2 = L3 = 0 and the A(q) reduce to

A(q) = A(0, 0, 0,q,q,q)

×



R(N1, N2, N3, q + 1, q + 1, q + 1)

(2q + 3)(2q + 3) − R(N1, N2, N3, q + 1, q + 1, q − 1)(2q + 3)(2q − 1)

− R(N1, N2, N3, q − 1, q + 1, q + 1)

(2q − 1)(2q + 3)

+

R(N1, N2, N3, q − 1, q + 1, q − 1)

(2q − 1)(2q − 1)



.

(13)

If in addition all orbitals are 1s-type, our I integral reduces (except for a factor) to the integral

Z(w1, w2, w3) =



r12r−1

13 r23 exp(−w1r1 − w2r2 − w3r3) dr1 dr2 dr3 (14)
Mathematical and computational science issues in high precision Hy-CI calculations 1523

that Remiddi has evaluated in closed form [15]. See also Sims and Hagstrom [10, 16] and

Harris et al [17] for correction of some minor misprints in the original published formula

[15]. We have found Remiddi’s result extremely useful by providing exact answers with

which to compare both our truncated expansions and the results of various series convergence

acceleration attempts (see below).

Our I integral (equations (11) and (12)), of course, is a more general form than the Remiddi

integral with Nis > 1 and non-s STOs. While parametric differentiation with respect to the wis

in equation (14) would generate closed form expressions for integrals with 1s2s, 2s2s etc charge

distributions, thus far no convenient analytic expressions for the non-s charge distribution cases

have been derived (although Fromm and Hill [18] suggested an elegant method for deriving

such expressions). Even if such expressions existed actual implementation would be beyond

the scope of the present investigation.

The angular terms A(L1, L2, L3, ns,q,nu) in equation (12) present no real problems. For

accurate work they do need to be calculated to greater than QP, but they can be precomputed

in advance of the integral calculation, stored by orbital l quantum numbers and used multiple

times. The R in equation (12) are the R(N1, N2, N3, N12,q,N23) given by equation (10),

a sum over 6 Wfgh(αβγ ) terms. The Wfgh(αβγ ) terms will be considered in the following

section while the question of the convergence of the A(q) summation in equation (11) will be

taken up in section 4.

3. Computation of the A, V, and W auxiliary functions

We now turn our attention to the W auxiliary functions in equation (10) along with two other

auxiliary functions, A and V , used in the evaluation of the W functions. Here, f, g and h are

integers, and α, β, and γ are real and positive. We supplement the definitions of the functions

by giving sufficient conditions on the parameters f, g and h to ensure the convergence of

the integrals. These supplementary conditions can be derived easily by inspection of the

problematic integration regions which occur when the integration variables x, y, and z assume

values close to zero.

W (fgh; αβγ ) ≡ Wfgh(αβγ ) =

 ∞

0

xf e−αx dx

 ∞

x

yg e−βy dy

 ∞

y

zhe−γ z dz,

(f  0, f + g  −1, f + g + h  −2) (15)

V (mn; αβ) ≡ Vmn(αβ) =

 ∞

0

xm e−αx dx

 ∞

x

yn e−βy dy, (m  0, m + n  −1)

(16)

A(n; α) ≡ An(α) =

 ∞

0

xne−αx dx = n!

αn+1 , n  0. (17)

For the triangle integral there are two cases to consider, Wfgh(αβγ ) for f, g, h  0 and

f, g  0,h< 0. Only the h < 0 case presents any difficulties.

3.1. Recurrence relationships for Wfgh(αβγ ) for f, g, h  0

Recurrence relationships between the W auxiliary functions have been given by James and

Coolidge [19] and Ohrn and Nordling [20], and included in our earlier integral paper [9].  ̈
1524 J S Sims and S A Hagstrom

(19)

g

f

h

(20)

(22) start

(21)

(18)

Figure 1. Recursion scheme for calculating Wfgh(αβγ ) for f, g, h  0.

Wfgh(αβγ ) = γ −1[hWf,g,h−1(αβγ ) + Vf,g+h(α, β + γ )],

(f  0, f + g  −1, f + g + h  −1) (18)

Wfgh(αβγ ) = (β + γ )−1

[gWf,g−1,h(αβγ ) + hWf,g,h−1(αβγ ) + Vf +g,h(α + β, γ )],

(f  0, f + g  0, f + g + h  −1) (19)

Wfgh(αβγ ) = (α + β + γ )−1

[f Wf −1,g,h(αβγ ) + gWf,g−1,h(αβγ ) + hWf,g,h−1(αβγ )],

(f  1, f + g  0, f + g + h  −1) (20)

Wf −1,g,h(αβγ ) = f −1

[αWfgh(αβγ ) + Vf +g,h(α + β, γ )],

(f  1, f + g  0, f + g + h  −1) (21)

W000(αβγ ) = [(α + β + γ )(β + γ )γ )]

−1

. (22)

Equations (18)–(22) enable one to calculate the W in a stable way using equation (22) as

the starting point (see figure 1). Each leg of the path is labelled with the formula we use to

raise/lower the corresponding index.

For checking the accuracy of the recursions the explicit finite sum of Frolov and Smith

[21], equation (16), is useful:

Wfgh(αβγ ) = 

h

σ=0

Cσ

n Aσ (γ )Vf,g+h−σ (α + β + γ,β + γ ) (23)

where the Cσ

n are binomial coefficients, i.e., Cmn = n!

m!(n−m)!

.

3.2. Recurrence relationships for Vmn(αβ) for m, n  0

In the previous section, we also needed arrays of Vmn(αβ) for m, n  0. In equation (16),

integration by parts using e−βy leads to [20]

Vmn(αβ) = β−1

[nVm,n−1(αβ) + Am+n(α + β)], (m  0, m + n  0) (24)
Mathematical and computational science issues in high precision Hy-CI calculations 1525

which is useful for raising the second index n. Equation (24) is obviously stable for n  0.

For n = 0 we get

Vm0(αβ) = β−1

[Am(α + β)], (m  0). (25)

Similarly, integrating by parts using e−αx in equation (16) gives a recurrence relation for

lowering the m index:

Vm−1,n(αβ) = m−1

[αVmn(αβ) + Am+n(α + β)], (m  1, m + n  0). (26)

The above formulae lead to the following scheme for the efficient computation of Vmn(αβ)

when m and n are  0: use equation (25) to get Vm0(αβ) for m = 0,...,mmax, then use

equation (24) for mmax to get all Vmmax,n(αβ), followed by equation (26) to lower m from mmax

down to 0 for each value of n. All coefficients enter with positive sign, hence it is a stable

scheme.

For checking the accuracy of the Vmn(αβ) recursions the explicit finite sum form of Frolov

and Smith [21] can be used:

Vmn(αβ) = n

ν=0

Cν

nAm+n−ν (α + β)Aν (β) (27)

where the Cm

n are binomial coefficients. Equation (27) follows immediately from the

parametric derivative representation of Vmn(αβ), namely,

Vmn(αβ) =



− d

dα

	m 

− d

dβ

	n

V00(αβ)

=



− d

dα

	m 

− d

dβ

	n

A0(β)A0(α + β)

=



− d

dβ

	n

A0(β)Am(α + β) (28)

and the binomial theorem for differentiation. This calculation is obviously stable as there is

no loss of precision due to differencing anywhere in the formulae.

3.3. Recurrence relationships for Wfgh(αβγ ) for f, g  0,h< 0

The calculation of Wfgh(αβγ ) for f, g  0,h< 0 presents some serious difficulties over the

range of the f, g, h indices and α, β, γ arguments that arise in practical problems. The basic

problem is that some of the recursion relations that must be used involve negative signs and

are hence potentially unstable. There is no really good, stable recursion scheme starting with

W0,0,−1(αβγ ), especially if more than double precision accuracy is needed, although earlier

work along this line has been extensive [19, 20, 22–24, 9]. Our earlier scheme for evaluating

these integrals was adopted by others, for example, [25], but it needs to be revised to meet

present high precision requirements.

We have found, surprisingly, that if one starts out by calculating the integrals

Wfmax,g,h0 (αβγ ) for each g = 0,...,gmax (where fmax and gmax are the maximum values

needed for f and g, and h0 is suitably chosen), then potentially unstable recursions in fact turn

out to be stable (how we judge the accuracy of the W integrals will be explained towards the

end of this section), and there is no real difficulty in generating the entire array from ‘the top

down’, using only a relatively few starting W integral values. More to the point, suppose we

want to calculate the array W (0: fmax, 0: gmax, hmin: − 1).

7 For that, the following recursion

scheme has proven to be remarkably effective:

7 An array of values of Wf,g,h(α, β, γ ) integrals for a fixed α, β, γ and f, g, and h taking on the values

f = 0,...,fmax, g = 0,...,gmax, and h = hmin,..., −1.
1526 J S Sims and S A Hagstrom

(i) Calculate a starting function Wfmax,gmax,h0 (αβγ ), using the stable Larsson summation (see

below), where hmin < h0 < −1 and with h0 suitably chosen (see below).

(ii) Use equation (18) with f = fmax, g = gmax

Wfmax,gmax,h(αβγ ) = γ −1[hWfmax,gmax,h−1(αβγ ) + Vfmax,gmax+h(α, β + γ )],

(f  0, f + g  −1, f + g + h  −1) (29)

to raise h from h0 up to the value −1. In spite of the fact that h is negative, this recursion is

stable up to fairly large values of fmax(≈200–300), provided h0 is appropriately chosen.

We invert the above formula and use

Wfmax,gmax,h−1(αβγ ) = [γ Wfmax,gmax,h(αβγ ) − Vfmax,gmax+h(α, β + γ )]/h (30)

to lower h down to the value hmin. Remarkably, this recursion also is stable. The

final result after two steps, therefore, will be the vector (one-dimensional array)

W (fmax, gmax, hmin : −1).

(iii) Repeat (i) and (ii) for g = gmax − 1, gmax − 2,..., 1, 0 giving one the subarray (twodimensional array) W (fmax, 0 : gmax, hmin : −1). Note that the existing V integral values

may be used over again (which was the reason for starting with g = gmax).

(iv) Finally, complete the calculation by lowering the value of the f index down to its lowest

allowed value f = max(0, | h | −g −2) (recall that f  0, f + g  −1, and f + g + h 

−2 for the W auxiliary function to exist) for each g, h pair using the stable recursion

equation (21)

Wf −1,g,h(αβγ ) = f −1

[αWfgh(αβγ ) + Vf +g,h(α + β, γ )]. (31)

This recursion is actually hyperstable in that it is possible to start out with W values not of full

accuracy and after a few steps obtain fully accurate W values, provided of course that the V

are fully accurate.

This is the scheme we use in practical calculations, although it is not the most efficient

scheme since W functions are computed that never arise in actual triangle integrals. For

example, Wfmax,g,−1 will never arise in any triangle integral, although it is the starting

element for the generation of the vector W (0: fmax, g, −1). Furthermore, this is an integral

that can result in exponent overflow for some argument values. A more efficient scheme

would endeavour to avoid calculating such elements. For example, suppose we partition the

W (0 : fmax, 0 : gmax, hmin : −1) matrix using h1 approximately midway between hmin and −1.

Calculate W (fmax, g, hmin: h1) as outlined above, then lower f to get W (f1: fmax, g, h1), where

f1 is the maximum value of f for this value of the indices g and h1. Then raise the h index to get

the vector W (f1, g, h1 : −1) and finally lower f from fmax (case h h1) and f1 (case h  h1)

using equation (31). In this manner calculating the subarray W (f1 +1: fmax, g, h1 +1: −1)

(which is never needed) can be avoided. Further partitioning of the W array is obviously

possible with further savings. The limiting process here would be some sort of ‘diagonal’

recursion relating Wfgh to its neighbours Wf +1,g,h−1 and Wf −1,g,h+1 along with the associated

V functions. We have looked at this process but have been unable to find a good starting point,

since one will presumably have to recur in both directions along the diagonal.

The starting Wfmax,g,h0 (αβγ ), g = 0,...,gmax values in step (i) above can be calculated

using the stable (sum of positive terms) summation formula derived many years ago by Larsson

[24], equation (32), based on the work of Ohrn and Nordling [20],  ̈

Wfgh(αβγ ) = 

M

ν=1

αν−1 f !

(f + ν)!

Vf +g+ν,h(α + β, γ ) (32)
Mathematical and computational science issues in high precision Hy-CI calculations 1527

for each g in the range (0, gmax).

8 M depends on the values of f, g, h and the argument

(α + β)/(α + β + γ ) and should be chosen large enough for convergence to occur, yet not too

large. How does one choose h0? In practice the simple rule h0 = hmin/2 works well for values

of hmin not too negative. For hmin in the range (−200–80) the rule 3

8hmin works fairly well.

There is also a slight dependence on f and g that can usually be ignored.

The Larsson summation in equation (32) requires knowledge of the vector V (fmax +

1: mmax, h0; α + β, γ ), where mmax = fmax + gmax + M, and M depends on h0 and the value of

the parameter ratio s = (α + β)/(α + β + γ ) as well as, obviously, on the type of arithmetic

being used (double precision (DP), quadruple precision (QP) or multiple precision (MP)). One

first calculates V (mmax, h0; α +β, γ ) and then uses equation (26) to lower the first index down

to fmax + 1. For mid-range values of s and typical h0 values M < 75, but as s approaches

1.0, M can be as large as 300.

In practice M is set to a value based on s (and other factors) that allows the series to

converge in most cases. However, if the series still does not converge, a special coding option

exists to sum the remainder of the series, but this time using MP arithmetic (the MP calculations

are done using Richard Brent’s multiple precision package [26]) to handle exponent overflow

problems that might otherwise occur for large values of M).

The Vmmax,h0 (α + β, γ ) value required to calculate V (fmax + 1: mmax, h0; α + β, γ ) is

calculated using another Larsson summation formula [24], (equation (34)):

Vmn(αβ) = 

ν=1

αν−1m!

(m + ν)!

Am+n+ν (α + β)

= α−1

Am+n(α + β) 

m + n + 1

m + 1

s +

(m + n + 1)(m + n + 2)

(m + 1)(m + 2)

s2 + ···

= Am+n+1(α + β)

(m + 1) 2F1(1, m + n + 2; m + 2;s) (33)

where s = α/(α+β) and 2F1(a, b; c;s) is the usual hypergeometric function [27]. The second

form shows the behaviour of the series with respect to s, with convergence becoming very

slow as s nears 1 and slower still if the ratio (m + n)/m is close to 1 (which it may be since

mmax  |h0|). On the other hand, for small s the summation length is largely independent

of m and n. Any of the forms in equation (33) can be used depending on the circumstances.

For convenience, over most of the range of s we just sum the series brute force. However,

for s < 0.3 an expansion in −z/(1 − z) (see equation (15.3.4) in Abramowitz and Stegun

[27]) is shorter and just as accurate as the canonical expansion. For s > 0.9 and the ratio of

(m + n + 1)/(m + 1) close to 1, we use the closed form expression due to McKoy [22]:

Vmmax,−1(αβ) = −Ammax (α) 


ln(1 − s) +

m

max

ν=1

(sν /ν), s = α/(α + β) (34)

followed by downward recursion using the reverse of equation (24)

Vm,n−1(αβ) = n−1

[βVmn(αβ) − Am+n(α + β)], (m  0, m + n  0) (35)

to lower the second index down to h0. This recursion is in practice very stable up to large values

of m for all allowed values of the second index. There is mild differencing in equation (34),

which unfortunately increases with increasing m, so for this case we use extended precision

(QP or MP) to calculate just the Vm,−1(αβ) function and only if full machine precision is

required.

An expansion in (1 − s) for use as s nears 1.0 can be obtained from equation (15.3.11) in

Abramowitz and Stegun [27], namely

8 Note that for most of our work g is in the range (0–20).
1528 J S Sims and S A Hagstrom

2F1(1, m + n + 2; m + 2;s) = −(s − 1)−n−1s−m−1(m + 2)log(1 − s)

(−n)(m + n + 2)

− (1 − δ(n, −1)) 


m + 1

n + 1

−

n−2

k=0

(m + n + 2)k(1 − s)k

(n + 2)k



+ (1 − δ(m, 0)) 
(s − 1)

−n−1 (m + 2)

(m + n + 2)

× ∞

k=0

(1 − s)k(m + 1)k(−n)k(ψ(0)(k + 1) − ψ(0)(k + m + 1))

k!(k − n − 1)!



, (36)

with (k), (a)k and ψ(k) the usual gamma function, Pochhammer symbol, and logarithmic

derivative of the gamma function, respectively. This complicated (and costly to implement)

expression has several problems, including the fact that the first and third terms have opposite

signs and the first sum contains alternating signs (both of which can lead to significant

differencing). Only in the s-range (0.97–1.0) is this approach superior (and then not always)

to the McKoy scheme, described above, which we use. See appendix B for further details

regarding equation (36).

Following a suggestion by King [28, 29], we have also looked into evaluating the

hypergeometric functions using Taylor series expansions about points s = 0.8(0.002)1.0

over a grid of m, n value pairs. A minimum of 1 function value and derivatives up to order 6

when s < 0.9, and a maximum 1 function value plus derivatives up to order 15 for s > 0.9 are

required at the DP level of accuracy and approximately double this number when using QP.

This approach is completely stable and a factor of 5–20 times faster than any of the methods

mentioned above. The problem with this approach is obviously that the memory requirements

can be quite large, depending on the number of m, n pairs involved (for example, anywhere

from 2–32 Mb in our studies).

Finally we note that Frolov and Smith [21] and Drake and Yan [30] have both employed

the hypergeometic form in equation (33) in their studies of the V and W integrals (Frolov and

Smith’s A2 and A3 integrals). In this connection, we note that Drake and Yan give an elegant

backward recursion relation for evaluating 2F1(1, m + n + 2; m + 2;s):

2F1(1, m + n + 1; m + 1;s) = 1 + m + n + 1

m + 1 	

s2F1(1, m + n + 2; m + 2;s). (37)

This relation is extremely simple to code, completely stable, is self starting (use 1.0), and for

a given s the starting value of m depends essentially only on n. Unfortunately the number

of iterations (hence the operation count) turns out to be very similar to the number of terms

needed when using the defining series for 2F1(1, m + n + 2; m + 2;s).

Once one has V (mmax, h0; α + β, γ ) in hand, the m index is lowered from mmax down to

(−h0 − 1) using the stable recursion relation equation (26)

Vm−1,n(αβ) = m−1

[αVmn(αβ) + Am+n(α + β)], (m  1, m + n  0). (38)

This completes the calculation of the starting function Wfmax,g,h0 (αβγ ).

The V functions in equations (29) and (30) are calculated by a scheme analogous to that

used for the W, namely, calculate a starting function Vmmax,n0 with n0 carefully chosen, then

recur up/down on n to get V (mmax, nmin : −1) followed by downward recursion on m to its

minimum allowed value of (−n − 1). The choice of n0 turns out to be a bit more difficult

than the choice of h0 for the starting W function. However, we have been able to show that

for m in the range (10–400), s in the range (0.01–0.97), and (α + β) in the range (1.0–300.0),

we can find values of n0 for a starting function V (mmax, n0) which will allow us to generate
Mathematical and computational science issues in high precision Hy-CI calculations 1529

the V integrals needed in building the W-array from only one Larsson summation and stable

recursion schemes. Further details are available on request from the authors.

Our recursion scheme is best illustrated with an example. Suppose we want the array

W (0 : fmax, 0 : 20, hmin : −1), that is, with g in the range (0–20).

(i) Calculate W (fmax, 20, h0) using the Larsson [24] summation given by equation (32).

This is the hard part as the V (m, n) series is usually slowly converging (and there may

be overflow problems to avoid). But note that the V needed in this summation can be

calculated by just one call to the routine which computes a fixed V using the Larssontype summation described in equation (33) (or, in special cases, by McKoy’s formula,

equation (34)), followed by the stable recursion relation given by equation (35).

(ii) Raise the third index to get W (fmax, 20, h0 +1: −1) using equation (29). This is a stable

recursion. Lower the third index to get W (fmax, 20, hmin : h0 − 1) using equation (30).

This is also stable over the α, β, γ values we have tried. Repeat (i) and (ii) for another

value of g. The Larsson summation will be easier this time as the same set of V can be

used over again.

(iii) For each h in the range (hmin : −1) and for each g in the range (0–20) one can lower the

first index to get W (−h−1 : fmax, g, h) using the stable recursion formula equation (31).

So, for g in the range (0–20), which is typical for us, one only has to do 21 Larsson

summations (non-trivial) followed by the usual recursion relations.

The accuracy of the W integrals is checked by comparing W values (the entire array or a

specified subarray) calculated in QP with accurate 40 decimal digit MP values obtained using

converged Larsson summation and outputting a table giving the number of integrals calculated

versus binary bits of accuracy of the QP values. Ordinarily, only the W for small values of

f, g, h need be checked in practice since these depend on the W for larger values of the indices

through the various recursion equations. For example, we typically only check the subarray

W (0: 8, 0: 5, −6: − 1) or perhaps a section along the diagonal W (f1: f + 1 + 10, g, h1: − 1),

for small values of h1 and selected g values and where f1 = max(0, −1−g, −2−g−|h|). As

an additional test MP checks were made for the subarray W (fmax −5: fmax, g, hmin: 1) to test

the choice of h0 and subsequent recursions to raise/lower h. The V were similarly compared

with MP values and in much the same way. For the entire V and W arrays, 30 place agreement

with the 40 decimal place MP results were routine.

Finally, one additional detail of the W-array calculation is worth mentioning. Floating

point exponent overflow is an occasional problem whenever (α + β + γ ) is small and f and

h are large in magnitude. Although clever programming can handle problems of this sort, by

far the easiest thing to do in practice is just scale the nonlinear parameters of W using the fact

that

Wfgh(αβγ ) = sf +g+h+3Wfgh(sα, sβ, sγ ) (39)

where s is the scale factor, which is most conveniently taken to be a power of two depending

on the particular case (it being trival to multiply QP numbers in-line by powers of two).

Hopefully the scaled W functions will behave more sensibly than the unscaled W. The cost of

this scaling is, therefore, just one extra (trivial) multiplication for each W function calculated,

a modest cost for the computational stability obtained. We have found this to be an extremely

easy way to eliminate the exponent overflow problem.

4. Convergence acceleration

In tables 1 and 2 we have summarized the results obtained for the series expansion of a typical

triangle integral given by equations (11) and (12) using the procedures described in section 3,
1530 J S Sims and S A Hagstrom

Table 1. Levin u-transform acceleration applied to the series expansion of the integral

I = 1s(r1)1s(r2)1s(r3)|r12r−1

13 r23|1s(r1)1s(r2)1s(r3), where 1s(r) = exp(−0.9375r) and

1s(r) = exp(−3.6875r). In the notation of Remiddi this is the integral Z(1.875, 4.625, 1.875)

given in equation (14).

N I (N ) Levin (0, N)

0 0.26243 75873 41781 604 × 10−2

1 0.26503 67394 61564 128 × 10−2

2 0.26505 79473 16383 076 × 10−2

3 0.26505 91749 37799 976 × 10−2

4 0.26505 93291 00921 617 × 10−2

5 0.26505 93590 49782 063 × 10−2

6 0.26505 93667 68479 784 × 10−2

7 0.26505 93691 95071 254 × 10−2

8 0.26505 93700 79194 118 × 10−2 0.26505 93707 72107 87554 95682 74743 × 10−2

9 0.26505 93704 40200 502 × 10−2 0.26505 93707 72116 63440 45833 86022 × 10−2

10 0.26505 93706 01604 334 × 10−2 0.26505 93707 72116 21988 03288 33951 × 10−2

11 0.26505 93706 79306 139 × 10−2 0.26505 93707 72116 14894 71291 41561 × 10−2

12 0.26505 93707 19083 962 × 10−2 0.26505 93707 72116 15185 68254 40808 × 10−2

13 0.26505 93707 40532 094 × 10−2 0.26505 93707 72116 15250 20996 26353 × 10−2

14 0.26505 93707 52622 044 × 10−2 0.26505 93707 72116 15248 33212 49204 × 10−2

15 0.26505 93707 59703 982 × 10−2 0.26505 93707 72116 15247 74062 99741 × 10−2

16 0.26505 93707 63994 105 × 10−2 0.26505 93707 72116 15247 74993 62547 × 10−2

17 0.26505 93707 66671 112 × 10−2 0.26505 93707 72116 15247 75519 00254 × 10−2

18 0.26505 93707 68386 069 × 10−2 0.26505 93707 72116 15247 75517 63156 × 10−2

19 0.26505 93707 69510 860 × 10−2 0.26505 93707 72116 15247 75513 13851 × 10−2

20 0.26505 93707 70264 354 × 10−2 0.26505 93707 72116 15247 75513 08210 × 10−2

21 0.26505 93707 70778 867 × 10−2 0.26505 93707 72116 15247 75513 15893 × 10−2

22 0.26505 93707 71136 357 × 10−2 0.26505 93707 72116 15247 75512 97711 × 10−2

Procedure I

Remiddi ‘exact’ result 0.26505 93707 72116 15247 75513 12672 × 10−2

I (1) + Levin (2, 7) 0.26505 93707 71028 69329 30647 87266 × 10−2

I (1) + Levin (2, 12) 0.26505 93707 72116 15155 44538 30885 × 10−2

I (1) + Levin (2, 17) 0.26505 93707 72116 15247 75216 36081 × 10−2

I (1) + Levin (2, 22) 0.26505 93707 72116 15247 75513 12900 × 10−2

I (1) + Levin (2, 27) 0.26505 93707 72116 15247 75513 12679 × 10−2

I (2) + Levin (3, 8) 0.26505 93707 71546 90581 96662 95490 × 10−2

I (2) + Levin (3, 13) 0.26505 93707 72116 15108 59017 54022 × 10−2

I (2) + Levin (3, 18) 0.26505 93707 72116 15247 76056 40655 × 10−2

I (2) + Levin (3, 23) 0.26505 93707 72116 15247 75513 12081 × 10−2

I (2) + Levin (3, 28) 0.26505 93707 72116 15247 75513 12671 × 10−2

I (3) + Levin (4, 9) 0.26505 93707 71846 06434 29947 87907 × 10−2

I (3) + Levin (4, 14) 0.26505 93707 72116 14785 52033 98432 × 10−2

I (3) + Levin (4, 19) 0.26505 93707 72116 15247 76077 76597 × 10−2

I (3) + Levin (4, 24) 0.26505 93707 72116 15247 75513 13088 × 10−2

I (3) + Levin (4, 29) 0.26505 93707 72116 15247 75513 12672 × 10−2

I (4) + Levin (5, 10) 0.26505 93707 71986 85529 15619 33481 × 10−2

I (4) + Levin (5, 15) 0.26505 93707 72116 14441 75105 49427 × 10−2

I (4) + Levin (5, 20) 0.26505 93707 72116 15247 74384 49217 × 10−2

I (4) + Levin (5, 25) 0.26505 93707 72116 15247 75513 22971 × 10−2

I (4) + Levin (5, 30) 0.26505 93707 72116 15247 75513 12672 × 10−2
Mathematical and computational science issues in high precision Hy-CI calculations 1531

Table 2. Generalized (Hurwitz) zeta function acceleration applied to the series expansion of the

integral 1s(r1)1s(r2)1s(r3)|r12r−1

13 r23|1s(r1)1s(r2)1s(r3), where 1s(r) = exp(−0.9375r) and

1s(r) = exp(−3.6875r). In the notation of Remiddi this is the integral Z(1.875, 4.625, 1.875)

given in equation (14).

N I (N ) I (N ) + Z(8,N)

10 0.26505 93706 01604 334 × 10−2 0.26505 93707 72114 98348 41908 17171 × 10−2

12 0.26505 93707 19083 962 × 10−2 0.26505 93707 72116 13547 76363 87170 × 10−2

14 0.26505 93707 52622 044 × 10−2 0.26505 93707 72116 15173 03512 92726 × 10−2

16 0.26505 93707 63994 105 × 10−2 0.26505 93707 72116 15240 48596 20462 × 10−2

18 0.26505 93707 68386 069 × 10−2 0.26505 93707 72116 15246 67337 92125 × 10−2

20 0.26505 93707 70264 354 × 10−2 0.26505 93707 72116 15247 54794 42339 × 10−2

30 0.26505 93707 71994 834 × 10−2 0.26505 93707 72116 15247 75470 85112 × 10−2

40 0.26505 93707 72099 009 × 10−2 0.26505 93707 72116 15247 75512 55363 × 10−2

50 0.26505 93707 72112 431 × 10−2 0.26505 93707 72116 15247 75513 10544 × 10−2

60 0.26505 93707 72115 089 × 10−2 0.26505 93707 72116 15247 75513 12524 × 10−2

70 0.26505 93707 72115 785 × 10−2 0.26505 93707 72116 15247 75513 12656 × 10−2

80 0.26505 93707 72116 006 × 10−2 0.26505 93707 72116 15247 75513 12670 × 10−2

84 0.26505 93707 72116 048 × 10−2 0.26505 93707 72116 15247 75513 12671 × 10−2

88 0.26505 93707 72116 076 × 10−2 0.26505 93707 72116 15247 75513 12672 × 10−2

100 0.26505 93707 72116 121 × 10−2

Remiddi ‘exact’ result 0.26505 93707 72116 15247 75513 12672 × 10−2

along with results obtained for three series acceleration methods. The integral chosen involves

only 1s-orbitals so that direct comparison can be made with Remiddi’s analytic result [15] for

the all 1s orbital case. Nonlinear parameters used are typical of what one would use in an

actual Be atom energy calculation. Orbital notation follows our earlier Be work [7].

The triangle integral expansion given by equations (11) and (12) is well known to converge

slowly (see column 2 in both tables). Here I(N) is the direct sum value of I when the series

for I (given by equation (11)) is truncated at q = N. Convergence is logarithmic as shown

by Larsson [24]. As can be seen from table 2, column 2, direct summation of the series can

give full DP accuracy, but at the QP (30 decimal place) level direct summation is clearly not

an attractive option (we are only at DP accuracy with 100 terms in the direct sum9). Hence

the need for some sort of extrapolation/acceleration. A number of methods exist that attempt

to accelerate the convergence, the most promising of which for this work seemed to be the

Levin u-transform [32], generalized zeta function extrapolation [33, 34], and van Wijngaarden

[35] transformation followed by Levin u-transformation. The respective details for all three

of these methods are given in the appendix.

Pelzl and King [12] have applied the Levin u-transformation and the van Wijngaarden

method to the related integral:

Z(−1, −1, −1; w1, w2, w3) =



r−1

12 r−113 r−123 exp(−w1r1 − w2r2 − w3r3) dr1 dr2 dr3 (40)

which converges more slowly than the integral discussed here, obtaining good results at

the DP and DP+ level of accuracy. Drake and Yan [30] have used the generalized zeta

function extrapolation for the same integral with satisfactory results. Our results for the Levin

and zeta function acceleration methods are shown in tables 1 and 2. The van Wijngaarden

transformation is not a viable option since it involves ‘looking ahead’ on the series to calculate

9 Frolov [31] has reported that QP accuracy requires about 12 000 terms.
1532 J S Sims and S A Hagstrom

additional A(q) elements needed for the extrapolation and we have no way to efficiently

calculate these additional elements without calculating all the intervening elements.

In table 1, column 2, we list I (N ), N = 0,..., 22, with I (22) accurate only to 11

decimals. However, one should note the initial rapid convergence of the series, with I (3)

being accurate to five digits. Column 3 gives the results of a Levin u-transform acceleration

using these 23 partial sum values (alternatively the 23 A(q) values). Specifically, Levin (0, N)

denotes the result for the integral I when Levin u-transform acceleration is applied to the

sequence A(q), q = 0,...,N. Also shown is the Remiddi ‘exact’ result (accurate to 30

decimal digits) from which we see that the Levin acceleration gives 25 to 26 place agreement

with the exact answer. Although not shown in table 1, still larger Levin (0, N) transforms

(i.e., for N > 22) actually give progressively worse results, a result which is well understood

to be caused by differencing inherent in the Levin u-transform procedure. A similar result

was noted by Drake and Yan [30] in their work. However, comparison with I (22) indicates a

remarkable 14 digits of accuracy gained by using the Levin procedure.

It is possible to do better still if we partition the series in equation (11) as follows:

I = I (K) + 

N

q=K+1

A(q) (41)

and then apply the Levin u-transform acceleration to the second factor, or ‘tail’ of the series.

The results of this procedure, labelled

I (K) + Levin (K + 1, N ),

are given in table 1 also, for K = 1,..., 4 and selected values of N (up to 30). It is seen that

essentially full agreement with the Remiddi result is obtained for K = 2, 3, 4 and N = K +26,

while the K = 1 case differs from the Remiddi result by only 7 in the 30th digit. The reason

for the improved agreement is obvious. Since I (K), k  2 is accurate to 5 digits and the A(q)

have all been computed to full accuracy, the erroneous digits in the Levin(K + 1,N) term,

when added to I (K), do not enter into the 30 digits of the final result. Comparing I (30) from

table 2 with I (4) + Levin (5, 30) we see that 18 digits have been gained by the acceleration.

In table 2 we present the results of a generalized zeta function extrapolation/acceleration

for this same integral. Column 2 gives the direct sum results I (N ) for selected values of N

up to N = 100. Column 3 lists values of I (N ) + Z(8,N), where Z(8,N) is the generalized

zeta function estimate for the I-series ‘tail’ ∞

q=N+1 A(q). Essentially full QP accuracy

(30 decimals) is obtained for N  80. A Z(7,N) extrapolation would work equally well here

(see the appendix for details). Note that zeta function extrapolation does not appear to have

the sort of inherent differencing errors one encounters with the Levin u-transform method.

On the other hand, the expansion length N is approximately three times as long as for the tail

Levin acceleration for the same QP level of accuracy, leading one to prefer Levin in practical

calculations, with the zeta function extrapolation serving mainly as a check on the Levin result.

Whereas the ‘cost’ of calculating the W-array is spread over the cost of calculating maybe

thousands of triangle integrals, each integral undergoing ‘acceleration’ incurs a separate, nontrivial fixed cost for the acceleration. For example, a Levin u-transform of the sequence

A(0), . . . , A(N ) requires (5N + 2) floating point arithmetic operations. Thus, the I (3) +

Levin (4, 29) scheme requires 5 × 25 + 2 = 127 arithmetic operations. On the other

hand, the I (88) + Z(8, 88) zeta function acceleration, which gives the same result, requires

2(8+1)(8+2) = 180 arithmetic operations. A Z(7, 88) extrapolation would require 144

operations. Here again Levin u-transformation would seem to be the method of choice. Of

course, if one does not need full QP accuracy, one can lower the ‘sizes’ of the extrapolations

for both methods and reduce the per integral acceleration cost accordingly.
Mathematical and computational science issues in high precision Hy-CI calculations 1533

Table 3. Series acceleration applied to the expansion of the integral Z(N1, N2, N3, w1, w2, w3) = 
 r12r−1

13 r23r

N1−1

1 r

N2−1

2 r

N3−1

3 exp(−w1r1 − w2r2 − w3r3) dr1 dr2 dr3.

N1 N2 N3 ω1 ω2 ω3 Procedure I

1 1 2 1.875 1.875 1.875 I (3) + L(4, 29) 0.13082 09812 08397 73522 35202 82063

1 1 2 1.875 1.875 1.875 I (88) + Z(8, 88) 0.13082 09812 08397 73522 35202 82063

1 1 3 1.875 1.875 4.625 I (3) + L(4, 29) 0.33721 55183 97029 92662 01726 35923 × 10−2

1 1 3 1.875 1.875 4.625 I (88) + Z(8, 88) 0.33721 55183 97029 92662 01726 35922 × 10−2

3 3 3 1.875 1.875 4.625 I (3) + L(4, 29) 0.60013 12193 11404 67291 98499 11050 × 10−1

3 3 3 1.875 1.875 4.625 I (88) + Z(8, 88) 0.60013 12193 11404 67291 98499 11051 × 10−1

1 1 3 1.875 1.875 7.375 I (3) + L(4, 29) 0.34451 77030 77120 20124 05543 27182 × 10−3

1 1 3 1.875 1.875 7.375 I (88) + Z(8, 88) 0.34451 77030 77120 20124 05543 27181 × 10−3

1 2 3 1.875 1.875 7.375 I (3) + L(4, 29) 0.84633 71300 85042 97747 68064 59143 × 10−3

1 2 3 1.875 1.875 7.375 I (88) + Z(8, 88) 0.84633 71300 85042 97747 68064 59142 × 10−3

3 3 3 1.875 1.875 7.375 I (3) + L(4, 29) 0.57815 58605 15428 39144 73631 36233 × 10−2

3 3 3 1.875 1.875 7.375 I (88) + Z(8, 88) 0.57815 58605 15428 39144 73631 36232 × 10−2

3 3 5 1.875 1.875 7.375 I (3) + L(4, 29) 0.32093 93187 20061 10583 85837 28259 × 10−2

3 3 5 1.875 1.875 7.375 I (88) + Z(8, 88) 0.32093 93187 20061 10583 85837 28260 × 10−2

3 5 5 1.875 7.375 7.375 I (3) + L(4, 29) 0.55456 45336 69859 54897 01541 65994 × 10−6

3 5 5 1.875 7.375 7.375 I (88) + Z(8, 88) 0.55456 45336 69859 54897 01541 65998 × 10−6

3 5 7 1.875 7.375 7.375 I (3) + L(4, 29) 0.64154 30023 06312 85317 08759 36157 × 10−6

3 5 7 1.875 7.375 7.375 I (88) + Z(8, 88) 0.64154 30023 06312 85317 08759 36158 × 10−6

In table 3 we give values for a selection of triangle integrals over Ns orbitals computed

using QP arithmetic corresponding to the integral

Z(N1, N2, N3, w1, w2, w3)

=



r12r−1

13 r23r

N1−1

1 r

N2−1

2 r

N3−1

3 exp(−w1r1 − w2r2 − w3r3) dr1 dr2 dr3. (42)

This table shows that the I (3) + Levin (4, 29) procedure works for QP accuracy in the

calculations for Ns orbitals as well. In this table L(4, 29) ≡ Levin (4, 29). Generalized

zeta function extrapolation (I (88) + Z(8, 88)) is used as the ‘exact’ result for estimating

the accuracy of the direct + tail Levin u-transforms for these integrals, since only in the all

1s-orbitals case can direct comparison with Remiddi’s analytic result be made. Comparison

with the Remiddi ‘exact’ value in table 2 indicates that comparison with the I (8) + Z(8, 88)

zeta function expansion is a good indicator of accuracy for these integrals. We judge these

results to be accurate to essentially 30 digits. In practical calculations we use an I (3) + Levin

(4, 29) acceleration procedure for QP accuracy calculations for triangle integrals over Ns

orbitals.

For non-s orbitals, the situation is more complicated. For all s orbitals, the series A(q)

in equations (11) and (12) is monotone decreasing starting with A(0). Integrals of the

form (ss,ss, nln

l), l > 0 behave similarly10, since such integrals are easily shown to be

proportional to the corresponding integral with nln

l replaced by nsns. But for integrals of

the form (sl,ss,sl), l = p, d,f, the A(q) series initially increases in value with the maximum

occurring at A(1) for p orbitals, A(2) for d-orbitals, A(3) for f orbitals, etc. This implies

that logarithmic convergence will not obtain until further out in the series, meaning one must

take a larger direct sum I (K) before applying an appropriate tail Levin transform. That is,

the acceleration procedure is now a function of the type of integral. This is shown in table 4,

10 The notation for the integrals lists the orbitals in charge distribution form, orbitals for electrons 1, 2 and 3,

respectively, separated by commas.
1534 J S Sims and S A Hagstrom

Table 4. Direct + tail Levin u-transform applied to the series expansion of the integrals

I = φ1(r1)φ2(r2)φ3(r3)|r12r−1

13 r23|φ

1(r1)φ2(r2)φ3(r3), where the s and non-s STOs φi(r)

are defined in equation (2). ‘Charge distrib’ lists the orbitals used in the integrals in charge

distribution form, the product of orbitals for electrons 1, 2 and 3, respectively, separated by

commas. α = 3.6875 for orbitals labelled with 0.9375 otherwise. ωi is the sum of the α for

electron i.

ω1 ω2 ω3 Charge distrib I

1.875 1.875 4.625 (1s2 p0, 1s1s, 1s2 p

0) 0.10066 84878 43444 84864 63592 909 × 10−2

1.875 1.875 4.625 (1s2 p0, 1s1s, 1s4 p

0) 0.11833 16774 13366 58115 19439 408 × 10−2

1.875 1.875 1.875 (1s2 p0, 1s1s, 1s6 p

0) 0.25743 89694 78797 30692 11152 245 × 10−2

4.625 1.875 4.625 (1s3 d

0, 1s1s, 1s3 d0) 0.14346 47386 47689 32252 51073 608 × 10−4

4.625 1.875 4.625 (1s3 d

0, 1s1s, 1s4 d0) 0.14703 92116 35576 65591 18724 228 × 10−4

4.625 1.875 4.625 (1s4 d

0, 1s1s, 1s4 d0) 0.16231 07687 95971 07033 07366 828 × 10−4

4.625 1.875 4.625 (1s4f 

0 , 1s1s, 1s4 f0 ) 0.90510 72608 44242 71161 49757 6 × 10−5

4.625 1.875 4.625 (1s4 f

0 , 1s1s, 1s5 f0 ) 0.11374 35574 81047 80859 46671 6 × 10−4

4.625 1.875 4.625 (1s5 f

0 , 1s1s, 1s5 f0 ) 0.15286 81862 39593 07275 31049 5 × 10−4

1.875 1.875 4.625 (2 p02 p0, 1s1s, 2 p02 p

0) 0.98583 53341 16785 59122 05368 963 × 10−2

1.875 1.875 4.625 (2 p02p0, 1s1s, 2 p03 p

0) 0.11134 31134 51262 35929 64403 064 × 10−1

1.875 1.875 4.625 (2 p02 p0, 1s1s, 3 p03 p

0) 0.15035 21732 19263 87212 62852 604 × 10−1

1.875 1.875 4.625 (3 d03 d0, 1s1s, 3 d03 d

0) 0.11449 39119 42576 97455 82583 88

1.875 1.875 4.625 (3 d03 d0, 1s1s, 3 d04 d

0) 0.18406 70664 81700 38830 73834 80

1.875 1.875 4.625 (3 d03 d0, 1s1s, 4 d04 d

0) 0.33719 54635 42604 81407 30305 24

where we list our best estimates for the values of a number of triangle integrals involving non-s

type orbitals. It can be seen that the accuracy for integrals containing p orbitals is 28 digits, 27

digits for d orbitals, and 26 digits when f orbitals are used. We point out that (pp, ss, pp) type

integrals have an estimated accuracy similar to that of (sd, ss, sd) integrals. This follows from

the fact that pp charge distributions can be expanded in terms of ss and sd charge distributions,

meaning that (pp, ss, pp) type integrals can be expressed in terms of (ss, ss, ss) and (sd, ss, sd)

type integrals, from which one can infer the probable accuracy of any acceleration attempt.

Similarly, (dd, ss, dd) can be expressed in terms of (ss, ss, ss), (sd, ss, sd) and (sg, ss, sg) type

integrals.

Also, the accuracy of the non-s orbital integrals is not full QP accuracy. We have done

extensive experimentation but have been unable to achieve QP accuracy or to understand why

the accuracy of the series acceleration falls off as the orbital l value increases. Fortunately, in

actual calculations these integrals do not need to be known to such high accuracy because they

ordinarily enter into matrix elements that are less important in the energy calculation. Certainly,

for the Be atom the (ss, ss, ss) type triangle integrals are far and away the most important.

However we are still looking into the problem, especially with respect to calculating the A

factors, more accurate zeta expansions (which we use to judge the tail Levin extrapolations),

and perhaps using quad-double (QD) to get some ‘known’ values against which we can

calibrate acceleration schemes.

For DP accuracy, table 1 suggests that either I (1) + Levin (2, 12) or I (2) + Levin (3, 13)

accelerations will give full DP accuracy (when the A(q) values are computed to (full) DP

accuracy). And in table 2, I (12) + Z(8, 12) also gives full DP accuracy. That is, 14 terms

in the I expansion are required for DP accuracy using Levin acceleration while 13 terms are

needed for the zeta function acceleration. For expansions of this length, many of the problems

mentioned in section 3.3 no longer obtain, including exponent overflow or the need for QP

arithmetic. Since both acceleration methods perform similarly at the DP level, both are made
Mathematical and computational science issues in high precision Hy-CI calculations 1535

available in our codes as user selected options. We also note that Levin (0,10) or Levin

(0,11), which give essentially DP accuracy, are not attractive options as the A(q) values must

be calculated in QP arithmetic for the acceleration to work. Finally we note that the above

observations are made with data from s orbital calculations. The problems with non-s orbitals

affecting the accuracy of the QP calculation do not affect the DP calculation of non-s orbitals.

These can be handled the same way as the all s orbital cases, except that we may need to go

to I (3) + Levin (4, 21) to compensate for the fact that the series is not monotone decreasing.

5. Discussion

The techniques we have discussed, which include doing integrals by blocks, recursive

calculation of V and W, special techniques for slowly convergent regions, direct plus Levin

u convergence acceleration, and minimizing the number of MP calculations are all important

in any attempt to extend high precision calculations to Be and beyond (and probably for Li

too). The W integrals have been discussed extensively in the past [19, 20, 22–24, 9, 36, 30,

21, 1], but no one has discussed them in the context of both accuracy (30 digits) and speed.

Calculation of the W to QP accuracy is absolutely essential in order for the triangle integral

direct sum to be Levin extrapolated.

It is also important to know how ‘good’ our extrapolated results really are, especially

in the absence of any Remiddi style analytic results with which to calibrate an acceleration

method. Certainly, agreement between Levin and Zeta function methods is one way to proceed,

assuming there is a systematic pattern to the acceleration data. But we still have a problem in

the non-s orbital cases where we have not been able to get as good a zeta extrapolation value as

that given by tail Levin extrapolation. We are investigating larger zeta function extrapolations

to get calibration values for the non-s cases. It may not be possible to get these two very

different extrapolation schemes to agree (zeta and tail Levin), of course, in which case what

we really need are analytic Remiddi style formulae. Or we could use the brute force methods

of Frolov and Bailey [1] (Frolov and Bailey do not discuss the recursive evaluation of the vast

majority of the W needed, a major part of this paper), there being little difference between an

all s-orbital and an (sp, ss, sp) integral. We cannot be really certain of the non-s results beyond

what is suggested by the ‘pattern’ of convergence, but the important all s-orbital case seems

solved to QP accuracy. It is important to note that either acceleration scheme will give DP

accuracy with no doubt whatsoever. With full DP accuracy, one can go a long way at the DP

level. Specifically, we should have accurate enough integrals for a 10−9 uncertainty in the total

energy.

There is considerable room for improving the W-array calculation. Specifically,

(i) Minimizing the calculation of W that are never used but arise because of the recursion

used.

(ii) More use of mixed mode arithmetic, i.e., calculate early terms in the I-series in QP but

the higher terms in DP.

(iii) Use of additional data types (beyond DP and QP) to essentially eliminate the need for

MP. Note that the time consumed in MP code is negligible overall since only a very small

number of W are calculated using Larsson summation and the rest of the W are calculated

on average with one floating point multiply and 1 floating point add plus an integer divide

(fast in QP) (or they are calculated in QP in the case of an DP accuracy calculation).

Software support for a number of new data types for Fortran 90 has recently become

available, including double precision with exponent (DPE), quadruple precision with

exponent (QPE) (for handling the exponent overflow problem), quad–double precision
1536 J S Sims and S A Hagstrom

(QD), and quad–double with exponent (QDE). And triple–double (TD) and triple–double

with exponent (TDE) are under development for general purpose MP [37].

(iv) It is important to tune the DP accuracy option for maximum speed, since the bulk of the

time will be spent in this mode. The overhead associated with QP arithmetic is such that

‘tuning’ is not really feasible in practice.

The results of our convergence acceleration study can be summarized as follows: Levin

u-transforms have problems because of alternating signs, zeta expansions are too long, and

B-transforms, while they solve the sign problem, introduce another problem. The direct + tail

Levin extrapolation used in this work largely overcomes all these problems.

The methods developed in this paper should be relevant to the calculation of the very

difficult X functions arising in conventional Hy-rij calculations. And they should be relevant

to the acceleration of the infinite series that arise in three and four-electron systems. A Larsson

summation formula exists for the X function in terms of W functions, so that one can start a

recursion scheme essentially identical to the one used here for the W. Of course, one needs

lots of W but we have an efficient scheme to get them.

Finally, we point out that the three-electron triangle integrals discussed here have been

the real bottleneck to highly accurate Hy-CI calculations. Now that this bottleneck has been

removed, doing really accurate calculations on atoms with N  5 becomes a real possibility.

This is unlike the situation with Hy-rij , where there are considerable difficulties already at

N = 5 due to the large number of five-electron integrals [38].

Acknowledgments

We are deeply indebted to Alexei M Frolov, Frederick W King and Krzysztof Pachucki for

helpful correspondence on this work. Ettore Remiddi and Robert B Bohn have offered kind

support. Ulrich Jentschura, Daniel Lozier, Peter Mohr and Roldan Pozo have all made helpful

comments on this work. We wish to thank the referee who kindly pointed out the approach

given by equation (36) as well as the further transformation given in appendix B. One of us

(JSS) would also like to thank Judith Devaney for supporting this endeavour.

Appendix A. Convergence acceleration methods

A.1. Levin u-transform

A Levin u-transformation [32] is perhaps the most widely known technique for accelerating

logarithmically convergent series. Consider equation (11), written as

I = ∞

q=0

A(q). (A.1)

This becomes I = 

k uk with the Levin u-transformation given by

uk =

k

j=0 cj (k, A(j ))Sj

k

j=0 cj (k, A(j )) . (A.2)

In equation (A.2) Sj is a partial sum of I given by

Sj = 

j

w=0

A(w) (A.3)
Mathematical and computational science issues in high precision Hy-CI calculations 1537

and

cj (k, A(j )) = (−1)

jCk

j (j + 1)

k−2

A(j )−1 (A.4)

with Ck

j the binomial coefficient given by

Ck

j = j !

k!(k − j )!

. (A.5)

Equation (A.2) is Levin’s u-transformation. Note that the Levin u-transformation has

coefficients that alternate in sign [39] which is responsible for the cancellation problems

in the Levin scheme.

A.2. van Wijngaarden transform

If the series of interest

S∞ = ∞

w=0

A(w) (A.6)

can be converted to a series of the form

S∞ = ∞

w=0

(−1)

wB(w) (A.7)

followed by a Levin u-transform of the new series, then the alternating sign in the new

series will cancel the alternating sign in equation (A.4) and the Levin ‘differencing’ problem

essentially goes away.

Pelzl and King [12] were able to come up with an expression for the B(w) which allowed

them to apply convergence acceleration techniques in a very effective manner by eliminating

the alternate sign problem. Jentschura, Mohr, Soff and Weniger [40], however, point out that

the Pelzl–King transform is a rederivation of an older transform of van Wijngaarden [35],

where

B(w) = ∞

t=0

2tA(2t(w + 1) − 1). (A.8)

The Pelzl–King function, as well as all the functions Jentschura et al consider, has first to be

converted into an alternating series. This involves sampling of the complete series by directly

(rather than recursively) calculating the terms with large index, as is seen from equation (A.8).

Thus this transform cannot be used in our situation as this involves knowing terms in the

expansion we do not have and cannot get in any reasonable manner. The reason the van

Wijngaarden transform gives such good results, of course, is just this ‘lookahead’ feature, but

it is also the feature that prevents it from being useful to us.

A.3. Generalized zeta function

Suppose one has the series

I = ∞

s=0

A(s) = 

N

s=0

A(s) + ∞

s=N+1

A(s) (A.9)

I = I (N ) + ∞

s=N+1

a(s)

s6 , (A.10)
1538 J S Sims and S A Hagstrom

where the rate of convergence is geometric. Larsson [24] has shown that the triangle integral

has a rate of convergence that has this asymptotic form for A(s), s  N + 1. If now we assume

that the analytic functions a(s) can be written

a(s) = 

L

n=0

α(n)

sn , (A.11)

then equation (A.10) becomes

I = I (N ) +

∞

N+1

1

s6



L

n=0

α(n)

sn (A.12)

= I (N ) +



L

n=0

α(n) ∞

s=N+1

1

sn+6

(A.13)

= I (N ) +



L

n=0

α(n)ζ (n + 6, N + 1), (A.14)

where ζ (n + 6, N + 1) is the generalized (Hurwitz) zeta function, which in this case is just

ζ (n + 6) minus the first N terms of the series. Observe that ζ (n + 6, N + 1) ≈ 1

(N+1)n+6 . The

rate of convergence therefore remains geometric but this rate can now be as large as desired

by taking N large enough. To determine the constants α(n), n = 0, 1,...,L we use the last

(L + 1) values of A(s) in the series for I and assume that

A(s) = 

L

n=0

α(n)

sn+6 , s = N − L, N − L = 1,...,N. (A.15)

This is a linear system that can be solved for the α(n). Let Ms(n) = 1

sn+6 be the coefficient

matrix. Then

A = Mα α = M−1

A (A.16)

L should be picked first and foremost to get the best acceleration but also to minimize the

cost of calculating the α(n) constants. In the present work L = 8 seemed to be optimum.

Once a value of L has been decided on, the Hurwitz functions in equation (A.14) (whose

computation is non-trivial) can be precomputed and entered into the program at compile

time. And, of course, the (L + 1) by (L + 1) inverse matrix M−1 is a one time only

calculation. Given M−1 and the Hurwitz function values, the per integral cost of a generalized

zeta function acceleration is 2(L + 1)(L + 2) floating point operations, a not inconsiderable

expense.

Note added in proof. Equation (36) and appendix B (below) are very recent developments aided by a referee’s

comments.

Appendix B. Comments on equation (36)

Both of the summations in equation (36) are of hypergeometric type and accordingly can

be collapsed to simpler closed form using any of the modern computer algebra systems.

Specifically, the second definite sum in equation (36) is reducible with a little prodding by

either Mathematica or Maple to the following:

−m

k=1

sk

k

sm+1(−n) (B.1)
Mathematical and computational science issues in high precision Hy-CI calculations 1539

Inserting this expression into equation (36) and collecting terms gives the simple closed form

expression

2F1(1, m + n + 2; m + 2;s) = −(s − 1)−n−1s−m−1(m + 2)

(−n)(m + n + 2)




log(1 − s) +

m

k=1

sk

k



− m + 1

n + 1

−

n−2

k=0

(m + n + 2)k(1 − s)k

(n + 2)k

(B.2)

which is valid for n −1 if we adopt the convention that there is no sum when the upper

limit is less than the lower limit. The similarity with the McKoy formula, equation (34), is

obvious. We have been unable to reduce the remaining indefinite sum in equation (B.2) to

a simpler form, although a number of different interesting expressions can be generated, all

of which (unfortunately) contain essentially (−n − 1) terms. That is, no real simplification

is obtained. On the other hand, we have not been able to show that a suitable ‘closed form’

expression for this sum is not possible. A detailed numerical comparison between equation

(B.2) and the combination McKoy formula followed by recursion to lower n gives essentially

identical results (with respect to accuracy and floating point operation counts). Finally, we

point out that for large values of m and for s near 1.0 the summation form in equation (36)

may in fact be computationally more efficient than its closed form counterpart.

References

[1] Frolov A M and Bailey D H 2003 Highly accurate evaluation of the few-body auxiliary functions and four-body

integrals J. Phys. B: At. Mol. Opt. Phys. 36 1857–67

[2] Bailey D H 1995 ACM Trans. Math. Soft. 21 379

[3] Sims J S and Hagstrom S A 2002 High-precision Hy-CI variational calculations for the ground state of neutral

helium and heliumlike ions Int. J. Quantum Chem. 90 1600–9

[4] Sims J S et al 2002 Accelerating scientific discovery through computation and visualization II J. Res. Natl. Inst.

Stand. Technol. 107 223–45

[5] Sims J S and Hagstrom S A 1975 Combined configuration-interaction-Hylleraas studies of atomic states: III.

the four lowest 2S and four lowest 2P states of lithium Phys. Rev. A 11 418–420

[6] Sims J S, Hagstrom S A, Munch D and Bunge C F 1976 Nonrelativistic energy of the ground state of Li− and

the electron affinity of lithium Phys. Rev. A 13 560–2

[7] Sims J S and Hagstrom S A 1971 A combined configuration interaction–Hylleraas wavefunction study of the

ground state of the beryllium atom Phys. Rev. A 4 908–916

[8] Sims J S and Hagstrom S A 1975 Combined CI-HY studies of atomic states: II. Compact wavefunctions for the

Be ground state Int. J. Quantum Chem. 9 149–56

[9] Sims J S and Hagstrom S A 1971 One center rij integrals over Slater-type orbitals J. Chem. Phys. 55 4699–710

[10] Sims J S and Hagstrom S A 2003 Comment: analytic value of the atomic three-electron integral with Slater

wavefunctions (Phys. Rev. A 44 5492 (1991)) Phys. Rev. A 68 016501

[11] Condon E U and Shortley G H 1963 The Theory of Atomic Spectra (Cambridge: Cambridge University Press)

[12] Pelzl P J and King F W 1998 Convergence accelerator approach for the high-precision evaluation of threeelectron correlated integrals Phys. Rev. E 57 7268–73

[13] King F W 1999 High-precision calculations for the ground and excited states of the lithium atom Adv. At. Mol.

Opt. Phys. 40 57–112

[14] Yan Z-C and Drake G W F 1997 Computational methods for three-electron atomic systems in Hylleraas

coordinates J. Phys. B: At. Mol. Opt. Phys. 30 4723–50

[15] Remiddi E 1991 Analytic value of the atomic three-electron correlation integral with Slater wavefunctions Phys.

Rev. A 44 5492–502

[16] Sims J S and Hagstrom S A 2003 Analytic value of the atomic three-electron integral with Slater wave functions

(Phys. Rev. A 68 016501) Phys. Rev. 68 059903

[17] Harris F E, Frolov A M and Smith V H Jr 2004 Comment on ‘Analytic value of the atomic three-electron

correlation integral with Slater wave functions’ Phys. Rev. A at press

[18] Fromm D M and Hill R N 1987 Analytic evaluation of three-electron integrals Phys. Rev. A 36 1013–44
1540 J S Sims and S A Hagstrom

[19] James H M and Coolidge A S 1936 On the ground state of lithium Phys. Rev. 49 688–95

[20] Ohrn Y and Nordling J 1963 On the calculation of some atomic integrals containing functions of  ̈ r12, r13, and

r23 J. Chem. Phys. 39 1864–72

[21] Frolov A M and Smith V H Jr 1997 Exact finite series for the few-body auxiliary functions Int. J. Quantum

Chem. 63 269–78

[22] McKoy V 1965 Closed form of infinite series used in some atomic integrals containing r12, r13, and r23 J. Chem.

Phys. 42 2959–61

[23] Byron F W Jr and Joachain C J 1966 Phys. Rev. 146 1

[24] Larsson S 1968 Calculations on the 2S ground state of the lithium atom using wave functions of the Hylleraas

type Phys. Rev. 169 49–54

[25] Clary D C and Handy N C 1976 CI-Hylleraas variational calculation on the ground state of the neon atom Phys.

Rev. A 14 1607–13

[26] Brent R 1978 A Fortran multiple-precision arithmetic package ACM Trans. Math. Software (TOMS) 4

57–70

[27] Abramowitz M and Stegun A 1964 Handbook of Mathematical Functions With Formulas, Graphs, and

Mathematical Tables (Washington, DC: Government Printing Office)

[28] King F W 2004 private communication

[29] Pelzl P J, Smethells G J and King F W 2002 Improvements on the application of convergence accelerators for

the evaluation of some three-electron atomic integrals Phys. Rev. E 65 (036707) (see p 036707-7, section B)

[30] Drake G W F and Yan Z-C 1995 Asymptotic-expansion method for the evaluation of correlated three-electron

integrals Phys. Rev. A 52 3681–5

[31] Frolov A M 2003 private communication

[32] Levin D 1973 Development of nonlinear transformations for improving convergence of sequencesInt. J. Comput.

Math. B 3 371–88

[33] Flajolet P and Vardi I 1996 Zeta function expansion of classical constants, Online publication available at:

http://www.lix.polytechnique.fr/ilan/publications.html

[34] Vardi I 1991 Computational Recreations in Mathematica (Reading, MA: Addison-Wesley)

[35] van Wijngaarden A 1965 Cursus: Wetenschappelijk Rekenen B Process Analyse (Amsterdam: Stichting

Mathematisch Centrum) pp 55–60

[36] Berk A, Bhatia A K, Junker B R and Temkin A 1986 Phys. Rev. A 34 4591

[37] Hagstrom S A unpublished

[38] Busse G, Kleindienst H and L  ̈ uchow A 1998 Nonrelativistic energies for the Be atom: double-linked Hylleraas–  ̈

CI calculation Int. J. Quantum Chem. 66 241–7

[39] Fessler T and Ford W F 1983 HURRY: an acceleration algorithm for scalar sequences and series ACM Trans.

Math. Software(TOMS) 9 346–54

[40] Jentschura U D, Mohr P J, Soff G and Weniger E J 1999 Convergence acceleration via combined nonlinearcondensation transformations Comput. Phys. Commun. 116 28–54