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

J. Phys. B: At. Mol. Opt. Phys. 40 (2007) 1575–1587 doi:10.1088/0953-4075/40/8/009

Mathematical and computational science issues in

high precision Hylleraas–configuration interaction

variational calculations: II. Kinetic energy and

electron–nucleus interaction 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

Received 22 January 2007, in final form 13 March 2007

Published 30 March 2007

Online at stacks.iop.org/JPhysB/40/1575

Abstract

Three-electron kinetic energy and electron–nucleus interaction integrals arising

in Hylleraas–configuration interaction (Hy–CI) calculations are discussed. We

show that the electron–nucleus interaction operators introduce no new integral

complications and that the kinetic energy operators lead to integrals which

have much in common with electron interaction integrals treated in the first

paper of this series. We discuss the efficient evaluation of the kinetic energy

integrals and give selected values of the integrals to 30 digit accuracy. The

one new W auxiliary function which arises in treating the kinetic energy

operator is discussed in appendix B. We also give a few details on accurate

evaluations of the Condon and Shortley cks which arise in Hy–CI calculations in

appendix A.

1. Introduction

In this paper, we continue the discussion, begun in the first paper of this series [1], of the

three-electron integrals arising in high precision energy level calculations on few electron

atoms and ions using the Hylleraas–configuration interaction (Hy–CI) formalism. The Hy–

CI method uses a wavefunction consisting of a linear combination of configuration state

functions multiplied by a single rij raised to some power. The basic idea is to combine the

relative computational ease and orbital picture of CI with the better convergence obtained by

introducing interelectronic coordinates rij directly into the wavefunction. We have recently

used this technique to determine extremely accurate nonrelativistic ground-state energies of

helium-like ions [2, 3]. As pointed out in the first paper of this series, the major difficulty

in Hy–CI calculations as far as integrals are concerned comes from the electron interaction

operators r−1

ij in the Hamiltonian leading to two-, three- and four-electron integrals. We

0953-4075/07/081575+13$30.00 © 2007 IOP Publishing Ltd Printed in the UK 1575
1576 J S Sims and S A Hagstrom

postpone discussion of four-electron integrals to a future paper [4] and continue the threeelectron discussion with a brief discussion of electron–nucleus interaction integrals and a

detailed treatment of the remaining kinetic energy integrals. For a very thorough treatment of

three-electron integrals over spherically symmetrical Slater-type orbitals (STOs), see [1–25]

of the paper by Pelzl and King [5] and [16–40] in the review article by King [6]. Yan and

Drake [7] treat the more general case of three-electron integrals over nonspherically symmetric

STOs. Ruiz [8] also discusses kinetic energy in Hy–CI. She transforms the kinetic energy

part of the Hamiltonian into the mutually independent distance coordinates ri and rij and the

polar angles θi and φi. Our treatment requires no transformations and leads to integrals which

are special cases of integrals which we have already treated [9, 1]. The special nature of

the Hy–CI expansion terms leads to much simpler integral formulae than the more general

relations of Yan and Drake [7], a consequence of the restriction of at most a single rij operator

in each Hy–CI term (in contrast to Hy-rij calculations where products of rij operators have to

be dealt with).

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

have been discussed previously [9]. As discussed in [9], all two- and three-electron integrals

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

auxiliary functions A, V and W [1, 9]. The auxiliary functions have been extensively discussed

(again see the references in the paper by Pelzl and King [5] and the review article by King [6]),

often under different names [11, 12]. Here we extend our paper I treatment of W-auxiliary

functions to include the case of Wfgh(αβγ ) for g < 0, f, h  0, a case which arises in treating

the kinetic energy integrals but does not arise in other three-electron interaction integrals. We

also briefly discuss the three-electron integrals involving the one-electron electron–nucleus

interaction operators (which present no difficulties). In this paper, we focus on the remaining

untreated three-electron integrals, those involving the one-electron kinetic energy operators.

Finally we give a few details on accurate evaluations of the Condon and Shortley cks in

appendix A.

2. Kinetic energy integrals

The nonrelativistic Hamiltonian H is, in atomic units (au),3

HNR = 

N

i=1

Hi +



i<j

r−1

ij , (1)

where Hi = Ti + Vi is a one-electron operator for electron i consisting of a kinetic energy part

Ti = −1/2∇2

i and a electron–nucleus interaction part Vi = −Z/ri.

In Hy–CI calculations, the major integral complications arise from the electron repulsion

operators r−1

ij which result in two-, three- and four-electron integrals that we have already

discussed [9, 1, 4]. The other operators in the Hamiltonian are electron–nucleus interaction

operators introducing no new integral complications and kinetic energy operators which lead

to integrals which have much in common with the electron interaction integrals, as will be

shown in this paper.

The secret to efficient calculation of the kinetic energy integrals is to avoid as much as

possible differentiating rij terms arising on the right-hand side of the kinetic energy matrix

elements. The method for doing this uses the Hermitian character of the kinetic energy

operator. There are three cases, namely,

3 The atomic unit of energy is chosen as μe4

h ̄ 2 = 1 au (of energy), where μ = memN /(me + mN ).
Math and computational science issues in high precision Hy–CI calculations 1577

Case φ1(r1)φ2(r2)|T1r12|φ

1(r1)φ2(r2)=φ1(r1)φ2(r2)r12|T1|φ1(r1)φ2(r2).

Case φ1(r1)φ2(r2)r12|T1 + T2|r12φ

1(r1)φ2(r2)=φ1(r1)φ2(r2)|φ1(r1)φ2(r2) + 1/2{φ1

(r1)φ2(r2)r2

12|T1 + T2|φ1(r1)φ2(r2) + φ1(r1)φ2(r2)r212|T1 + T2|φ1(r1)φ2(r2)} where we

have used the transformation of Kolos and Roothaan [13].

Case φ1(r1)φ2(r2)φ3(r3)r12|T1 + T2 + T3|r13φ

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

r12T1r13|φ

1(r1)φ2(r2)φ3(r3)+φ1(r1)φ2(r2)φ3(r3)|r12r13|T2φ1(r1)φ2(r2)φ3(r3)+φ1(r1)

φ

2(r2)φ3(r3)|r13r12|T3φ1(r1)φ2(r2)φ3(r3).

The first two cases involve only products of elementary one-electron integrals. In the

third case, the connection with the electron interaction integrals φ1(r1)φ2(r2)φ3(r3)|r12r13|

φ

1(r1)φ2(r2)φ3(r3) is apparent for all except the first term. We show below that this

connection also extends to this more complicated term even after the differentiations involved

with T1.

3. IT1 = φ1(r1)φ2(r2)φ3(r3)|r12T1r13|φ

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

The φ(r) and φ

(r) are unnormalized, nonorthogonal Slater-type orbitals (STOs) (here, and in

equations which follow, we drop parameter dependences where that would make the formulae

excessively cluttered)

φ(r) = φ(r, θ , φ) = rn−1 e−αrY m

l (θ , φ),

φ

(r) = φ(r, θ , φ) = rn

−1 e−αr

Y m

l (θ , φ),

(2)

where the radial factor depends on the orbital exponent α and the radial quantum number

n. The l and m quantum numbers define the order and degree of the orthonormal spherical

harmonics Y m

l defined in the Condon and Shortley phase convention [10].4

We use a multipole expansion for ∇2

i φ1(r1)r13 and do the resulting integral in such a way

that it is a special case of integrals we have already treated [1]. In the following l

2 is the

squared resultant one-electron orbital angular momentum operator l

2 = l2

x + l

2

y + l

2

z which,

when operating on a φi(ri), has the property that

l

2

i φi(ri) = li(li + 1)φi(ri). (3)

Writing the one-electron operator T1 (equation (1)) as

T1 = −1

2

∂2

∂r2

1

+

2

r1

∂

∂r1



+

l

2(1)

2r2

1

= T 0

1 +

l

2(1)

2r2

1

, (4)

IT1 can be expressed as the sum of two integrals, I 

T1 and I T1

, given by

I 

T1 = φ1(r1)φ2(r2)φ3(r3)|r12T 0

1 r13|φ1(r1)φ2(r2)φ3(r3) (5)

and

I 

T1 = φ1(r1)φ2(r2)φ3(r3)|r12

l

2(1)

2r2

1

r13|φ

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

4 The STOs we use are defined fully in [9]. An s-type STO has l = 0, a p STO has l = 1, a d STO has l = 2, etc. A

p0 STO has m = 0.
1578 J S Sims and S A Hagstrom

3.1. I 

T1 = φ1(r1)φ2(r2)φ3(r3)|r12T 0

1 r13|φ

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

Consider

T 0

1 φ1(r1)r13 =



−1

2

∂2

∂r2

1

+

2

r1

∂

∂r1

 φ

1(r1)r13. (7)

To carry out the derivatives, we use the expansion for r13 derived in [9],

r13 = ∞

n=0

Bn(r13, 1)Pn(cos θ13)

= ∞

n=0

Pn(cos θ13)



1

m=0

Bn,1 m

rn+2−2m

13<

rn+1−2m

13>

, (8)

with

Bn,1

0 = 1

2n + 3, Bn,1

1 = − 1

2n − 1

. (9)

This is just a particular case of an expansion given by Perkins [14], which is a rewriting of the

well-known expansion of Sack [15]. Here r13< and r13> refer to the lesser and greater of r1

and r3, respectively. The term T 0

1 φ

1(r1)r13 becomes, after taking the derivatives,

T 0

1 φ1(r1)r13 = ∞

n=0

Pn(cos θ13)



−1

2

∂2

∂r2

1

+

2

r1

∂

∂r1

 φ

1(r1)



1

m=0

Bn,1 m

rn+2−2m

13<

rn+1−2m

13>

= φ

1(r1)

∞

n=0

Pn(cos θ13)



1

m=0

Bn,1 m

×



rn+2−2m

1

rn+1−2m

3



−(n

1 + n + 2 − 2m)(n1 + n + 1 − 2m)

2r2

1

+

α

1(n1 + n + 2 − 2m)

r1

− α2

1

2



r1<r3

+ rn+2−2m

3

rn+1−2m

1



−(n

1 − n − 1+2m)(n1 − n − 2+2m)

2r2

1

+

α

1(n1 − n − 1+2m)

r1

− α2

1

2



r3<r1

	

. (10)

Comparing equation (10) with the expansion of r13 in equation (8) makes it clear that

· · · |r12T 0

1 r13| · · · leads to the same kinds of integrals as does · · · |r12r13| · · ·. Indeed, our

integral could be expressed as a linear combination of appropriate · · · |r12r13| · · · integrals

but this is not the approach usually taken. Now introduce ‘charge distributions’

i(ri) = φ∗

i (ri)φi(ri) (11)

which we expand in terms of STO-like functions fi(ri)

fi(ri) = r

Ni−1

i e−wiri

Y Mi

Li (θi, φi); (12)

the formula is

i(ri) = 

Li

(2Li + 1)

4π

1/2

cLi(l

i, mi;li, mi)fi(ri), (13)
Math and computational science issues in high precision Hy–CI calculations 1579

where Ni = ni + n

i − 1, wi = αi + αi, and Li satisfies the triangular inequality

|li − l



i| Li li + l



i with li + l



i +Li = 2g (g integral) and with the further restriction

Li  |Mi|, Mi = m

i − mi. The cLi are standard Condon and Shortley coefficients (the

so-called cks) which we discuss in appendix A.

Using the spherical harmonic addition theorem [10] (see the Addenda)

Pn(cos θ13) = 4π

2n + 1

n

mn=−n

Y mn∗ n (θ1, φ1)Y mn

n (θ3, φ3), (14)

expanding r12 using equation (8), and expanding the products of STOs using equation (13) we

get

I 

T1 = 

L1,L2,L3




3

i=1

{(2Li + 1)

1/2

cLi(l

i, mi;li, mi)}J 2(L1, L2, L3), (15)

where

J 

2(L1, L2, L3) = 1

(4π ) 3

2



r12T 0

1 r13 


3

i=1

{fi(ri) dri}

= A

2(L1, L2, L3)



1

l=0



1

m=0

BL2,1

l BL3,1 m

× R(N1, N2, N3, N12(l), N13(m); w1, w2, w3) (16)

with N12(l) = L2 + 1 − 2l,N13(m) = L3 + 1 − 2m, Ni = ni + n

i − 1, wi = αi + αi, and

A

2(L1, L2, L3) = δ(M1 + M2 + M3, 0)(−1)

M2 cL3 (L2, −M2;L1, M1)

(2L2 + 1)(2L3 + 1)

1

2

. (17)

A

2(L1, L2, L3) will be 0 unless the triangular inequality |L2 − L1| L3 L1 + L2,

with L1 + L2 + L3 = 2g (g integral), is satisfied and the sum of the Mi equals 0.

R(N1, N2, N3, N12(l), N13(m); w1, w2, w3)) is given by

R(N1, N2, N3, N12(l), N13(m); w1, w2, w3) =

 r

N12(l)+1

12<

r

N12(l)

12>

×



r

N13(m)+1

1

r

N13(m)

3



−(n

1 + N13(m) + 1)(n1 + N13(m))

2r2

1

+

α

1(n1 + N13(m) + 1)

r1

− α2

1

2



r1<r3

+ r

N13(m)+1

3

r

N13(m)

1



−(n

1 − N13(m))(n1 − N13(m) − 1)

2r2

1

+

α

1(n1 − N13(m))

r1

− α2

1

2



r3<r1

	

× e−w1r1−w2r2−w3r3 r

N1+1

1 r

N2+1

2 r

N3+1

3 dr1 dr2 dr3. (18)

Let

R1(N1, N2, N3, N12(l), N13(m), p; w1, w2, w3)

= W (N1 +3+ N12(l) + N13(m) + p, N2 + 1 − N12(l), N3 + 1 − N13(m); w1, w2, w3)

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

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

(19)
1580 J S Sims and S A Hagstrom

and

R2(N1, N2, N3, N12(l), N13(m), p; w1, w2, w3)

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

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

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

(20)

where the W integrals are defined by

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

 ∞

0

xf e−αx dx

 ∞

x

yg e−βy dy

 ∞

y

zh e−γ z dz,

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

The W auxiliary integrals are discussed in [1]. Then

R(N1, N2, N3, N12(l), N13(m); w1, w2, w3)

= −1/2(n

1 + N13(m) + 1)(n1 + N13(m))R1(N1, N2, N3, N12(l), N13(m),−2; w1, w2, w3)

− 1/2(n

1−N13(m))(n1−N13(m)−1)R2(N1, N2, N3, N12(l), N13(m),−2; w1, w2, w3)

+ [α

1(n1 + N13(m) + 1)]R1(N1, N2, N3, N12(l), N13(m), −1; w1, w2, w3)

+ [α

1(n1 − N13(m))]R2(N1, N2, N3, N12(l), N13(m), −1, w1, w2, w3)

− α2

1

2 {R1(N1, N2, N3, N12(l), N13(m), 0; w1, w2, w3)

+ R2(N1, N2, N3, N12(l), N13(m), 0; w1, w2, w3)}. (22)

Substituting equation (16) into equation (15) and rearranging terms, I 

T1 becomes

I 

T1 = δ(M1 + M2 + M3, 0)(−1)

M2



L1,L2,L3

A2(L1, L2, L3)

×



1

l=0



1

m=0

BL2,1

l BL3,1 m R

(N1, N2, N3, N12(l), N13(m); w1, w2, w3), (23)

where

A2(L1, L2, L3) =

(2L1 + 1)

(2L2 + 1)

1

2

cL3 (L2, −M2;L1, M1)




3

i=1

cLi(l

i, mi;li, mi)}, (24)

and |L2 − L1| L3 L1 + L2, with L1 + L2 + L3 = 2g (g integral).

3.2. I 

T1 = φ1(r1)φ2(r2)φ3(r3)|r12

l

2(1)

2r2

1

r13|φ

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

Note that we cannot immediately expand the orbital product for electron 1 since we must first

operate with l

2(1) on φ

1(r1)r13. Expanding the charge distributions for electrons 2 and 3, I 

T1

becomes

I 

T1 = 

L2,L3




3

i=2



(2Li + 1)

1

2 cLi

(l

i, mi;li, mi)




J 

2 (L2, L3), (25)

where J 

2 (L2, L3) is given by

J 

2 (L2, L3) = 1

4π



φ∗

1 (r1)r12

l

2(1)

2r2

1

r13φ

1(r1)




3

i=2

{fi(ri) dri} dr1. (26)
Math and computational science issues in high precision Hy–CI calculations 1581

To carry out the operation of l

2(1) on φ

1(r1)r13, we again use the expansion for r13

(equation (8)), so that

l

2

(1)φ

1(r1)r13 = ∞

n=0

Bn(r13, 1)l

2

(1)Pn(cos θ13)φ

1(r1). (27)

Using the spherical harmonic addition theorem (equation (14)) and the formula for expanding

a product of spherical harmonics (one of which is complex conjugated),

Y m∗ l (θ , φ)Y m

l (θ , φ) = 

L

2L + 1

4π

1

2

cL(l, m;l, m)Y M

L (θ , φ), (28)

we get

Pn(cos θ13)φ

1(r1) = {4π}

1

2

2n + 1r

n

1−1

1 e−α

1r1

n

mn=−n

Y mn

n (θ3, φ3)

×

L

{2L + 1}

1

2 cL(l

1, m1; n, mn)Y ML (θ1, φ1), (29)

where L satisfies the triangular inequality |l1 − n| L l1 + n, with l1 + n + L = 2g

(g integral) and with the further restriction L  |M|, M = m

1 − mn. But

l

2

Y M

L (θ , φ) = L(L + 1)Y ML (θ , φ), (30)

so

l

2

(1)φ

1(r1)r13 = ∞

n=0

Bn(r13, 1)l

2

(1)Pn(cos θ13)φ

1(r1)

= {4π}

1

2 r

n

1−1

1 e−α

1r1

∞

n=0

n

mn=−n

1

(2n + 1)

Bn(r13, 1)Y mn

n (θ3, φ3)

×

L

L(L + 1){2L + 1}

1

2 cL(l

1, m1; n, mn)Y ML (θ1, φ1). (31)

Substituting this expression into equation (26), expanding r12, and integrating over angular

coordinates leads to

J 

2 (L2, L3) = A(L2, L3)



1

l=0



1

m=0

BL2,1

l BL3,1 m R(N1 − 2, N2, N3, N12(l),

N13(m); w1, w2, w3), (32)

where N12(l) = L2 + 1 − 2l and N13(m) = L3 + 1 − 2m as before and

A

2(L2, L3) = δ(M1 + M2 + M3, 0)(−1)

M3



L

L(L + 1)(2L + 1)

2(2L2 + 1)(2L3 + 1)

× cL(l

1, m1;L3, −M3)cL(l1, m1;L2, M2). (33)

L needs to satisfy all of the following conditions:

(i) |l



1 − L3| L l



1 + L3, l1 + L3 + L = 2g (g integral)

(ii) |l1 − L2| L l1 + L2, l1 + L2 + L = 2g (g integral)

(iii) L  |m

1 + M3|, L  |m1 − M2|,L> 0
1582 J S Sims and S A Hagstrom

and the range of L from the first condition must overlap the range of L from the second

condition. For most L2, L3 cases arising in practice, the range of L will be quite limited.

R(N1 − 2, N2, N3, N12(l), N13(m); w1, w2, w3) is given by

R(N1 − 2, N2, N3, N12(l), N13(m); w1, w2, w3) =



r

N1−1

1 r

N2+1

2 r

N3+1

3

r

N12+1

12< r

N13+1

13<

r

N12

12>r

N13

13>

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

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

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

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

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

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

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

R is a special case of a more general formula for three-electron interaction radial integrals

R(N1, N2, N3, N12, N13, N23, s,t, u; w1, w2, w3) discussed elsewhere [9] which handles the

operator rs

12rt13/ru23 and reduces to equation (34) for s = t = 1, N23 = 0 and u = 0. Note that

in R(N1 − 2, N2, N3, N12(l), N13(m); w1, w2, w3) the power that r1 is raised to is reduced

by 2 due to the 1

2r2

1

factor; i.e., the radial integral is obviously the same sort of term as arises

in · · · |r12r13| · · · except for N1 being replaced by N1 − 2. This reduction of N1 by 2 has an

effect on the Wfgh(αβγ ) function evaluation, namely the case g < 0, f, h  0 can arise now,

whereas it cannot in the other three-electron interaction integrals treated in the first paper in

this series [1]. In appendix B, we provide a suitable recursion scheme for this new W-function

case.

Our final expression for I 

T1 is thus

I 

T1 = δ(M1 + M2 + M3, 0)(−1)

M3



L2,L3

A

T1

(L2, L3)

×



1

l=0



1

m=0

BL2,1

l BL3,1 m R(N1 − 2, N2, N3, N12(l), N13(m); w1, w2, w3) (35)

where

A

T1

(L2, L3) =

3

i=2{cLi(l

i, mi;li, mi)}

2{(2L2 + 1)(2L3 + 1)}

1

2

×

L

L(L + 1)(2L + 1)cL(l

1, m1;L3, −M3)cL(l1, m1;L2, M2) (36)

with the conditions on L given above.

IT1 is the sum of two integrals, I 

T1 and I T1

, treated above. In table 1 we give values

for the kinetic energy operator T1 operating on r13 and both spherically and non-spherically

symmetric STOs computed using quadruple precision (QP) arithmetic (when we refer to

double, quadruple, or other precision it is with respect to a 32 bit word). In this connection,

software support for a number of new data types for Fortran 90 has recently become available,

including double precision with exponent (DPE) and quadruple precision with exponent (QPE)

(for handling the exponent overflow problem) as well as quad-double precision (QD) and quaddouble with exponent (QDE). As a check on the accuracy of the results tabulated in table 1

we computed all integrals in QD and found that the QP results agree with the QD results with

an error of at most 1–2 in the last digit reported.
Math and computational science issues in high precision Hy–CI calculations 1583

Table 1. Selected values of the integral IT1 = φ1(r1)φ2(r2)φ3(r3)|r12T1r13|φ

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

computed in quadruple precision, where the s and non-s STOs φi(r) are defined in equation (2).

‘Charge distribution’ lists the orbitals used in the integrals in charge distribution form (e.g.,

equation (11)), the product of orbitals for electrons 1, 2 and 3, respectively, separated by commas.

α = 1.4 for orbitals labelled with , 2.86 otherwise. ωi is the sum of the α’s for electron i. T1

operates on both r13 and φ

1(r1).

ω1 ω2 ω3 Charge distribution IT1

5.72 4.26 4.26 (1s1s,1s2s,2s1s) 0.13126 72235 19561 47025 67125 92522 × 10−4

5.72 5.72 2.8 (1s1s,1s1s,2s2s) 0.64207 60509 60589 25934 71629 40533 × 10−4

2.8 5.72 2.8 (2p

02p0,1s1s,2s2s) 0.57720 73518 58580 32536 07477 41306 × 10−3

5.72 4.26 2.8 (1s1s,1s2p

0,2s2p0) 0.38140 71178 87570 42824 91516 47944 × 10−5

5.72 5.72 4.26 (1s2p0,1s1s,1s4p

0) −0.14502 00857 39528 49418 94063 13172 × 10−6

5.72 5.72 4.26 (1s2p0,1s1s,1s6p

0) −0.32844 14491 49562 97251 56025 55653 × 10−6

4.26 5.72 4.26 (1s3d

0,1s1s,1s3d0 ) 0.79465 21564 41320 36387 60727 96091 × 10−7

4.26 5.72 4.26 (1s3d

0,1s1s,1s4d0 ) 0.12075 23527 76977 23170 64933 06356 × 10−6

4.26 5.72 4.26 (1s3d

0,2p03p0,1s4d0) 0.67304 67088 08933 49971 96261 05823 × 10−7

5.72 5.72 4.26 (2p02p0,1s1s,2s2s) 0.29871 77066 84990 94807 38892 10487 × 10−5

5.72 5.72 4.26 (2p02p0,1s1s,2p02p

0) 0.29519 96062 34340 69487 85811 15133 × 10−5

5.72 5.72 4.26 (2p02p0,2p03p0,2p03p

0) 0.19213 21210 92616 69396 65303 73606 × 10−5

5.72 5.72 4.26 (3p03p0,1s1s,3p03p

0) 0.40535 29068 76436 24808 08465 07749 × 10−5

5.72 5.72 4.26 (3d03d0,1s1s,3d03d

0) 0.87231 87107 67592 04016 01662 18286 × 10−5

5.72 5.72 4.26 (3d03d0,1s1s,3d04d

0) 0.15815 57944 02140 20530 85237 26766 × 10−4

5.72 5.72 4.26 (3d03d0,1s1s,4d04d

0) 0.32596 11204 87135 13502 55123 63283 × 10−4

5.72 5.72 4.26 (3d03d0,2p03p0,4d04d

0) 0.13472 45687 15817 35201 25998 75627 × 10−4

5.72 5.72 4.26 (3d03d0,3d04d0,4d04d

0) 0.20760 58133 97401 23636 32696 31433 × 10−4

5.72 5.72 4.26 (2p12p1,1s1s,2p12p

1) 0.29783 81815 72328 38477 50621 86648 × 10−5

5.72 5.72 4.26 (3d23d2,1s1s,3d23d

2) 0.87337 99889 07121 70685 75959 02697 × 10−5

5.72 5.72 4.26 (2p12p1,2p12p1,2p12p

1) 0.14390 09188 74908 52676 79501 50271 × 10−5

5.72 5.72 4.26 (3d23d2,3p13p1,3d23d

2) 0.41359 57078 58469 80859 64162 38576 × 10−5

5.72 5.72 4.26 (2p02p0,2p02p−1,2p02p

1) −0.75104 03404 78207 96905 89118 93078 × 10−9

5.72 5.72 4.26 (3d03d0,3p03p−1,3d04d

1) −0.11247 28271 12214 02357 71049 25523 × 10−8

4. Nuclear attraction integrals

The electron–nucleus interaction part of the one-electron operator for electron i is Vi = −Z/ri.

The three-electron integrals involving the one-electron electron–nucleus interaction can be

treated like electron interaction integrals with a modified ‘charge distribution’; i.e., a separate

option in a subroutine for · · · |r12r13| · · · can handle the · · · |r12( 1

r1 + 1r2 + 1r3

)r13| · · ·.

Alternatively, in the kinetic energy formulae the kinetic energy operator T1 can be converted

into a corresponding H1 operator incorporating electron–nucleus interaction, resulting in

formulae for the one-electron Hamiltonian operator H1 which reduce to the kinetic energy

formulae we have given in the special case of Z = 0; i.e., T1 = H1 (Z = 0). In

equation (18) one simply adds −Z to the numerators of both r−1

1 terms to define a

R

Z(N1, N2, N3, N12(l), N13(m); w1, w2, w3)). That is

R

Z(N1, N2, N3, N12(l), N13(m); w1, w2, w3)

= −1/2(n

1 + N13(m) + 1)(n1 + N13(m))R1(N1, N2, N3, N12(l), N13(m),−2; w1, w2, w3)

− 1/2(n

1−N13(m))(n1−N13(m)−1)R2(N1, N2, N3, N12(l), N13(m),−2; w1, w2, w3)

+ [α

1(n1 + N13(m) + 1) − Z]R1(N1, N2, N3, N12(l), N13(m), −1; w1, w2, w3)
1584 J S Sims and S A Hagstrom

+ [α

1(n1 − N13(m)) − Z]R2(N1, N2, N3, N12(l), N13(m), −1, w1, w2, w3)

− α2

1

2 {R1(N1, N2, N3, N12(l), N13(m), 0; w1, w2, w3)

+ R2(N1, N2, N3, N12(l), N13(m), 0; w1, w2, w3)}. (37)

Using this R

Z instead of R in equation (23) yields an expression for I 

H1 which can be added

to I 

T1 to form IH1 .

5. Discussion

In this paper, we have completed the discussion, begun in the first paper of this series [1], of

the three-electron integrals arising in high precision energy level calculations on few electron

atoms and ions using the Hylleraas–configuration interaction (Hy–CI) formalism. We provide

a brief discussion of electron–nucleus interaction integrals and a detailed treatment of the

remaining untreated three-electron integrals, those involving the one-electron kinetic energy

operators.

The formulae we have given can be easily generalized to handle the integrals arising in

a strictly Hylleraas calculation. We have deliberately chosen not to do so, since the resulting

formulae are unnecessarily complicated for coding purposes when doing Hy–CI (due to the

special nature of Hy–CI expansion terms). With the present formulation, the three-electron

kinetic energy integrals are computed as special cases of three-electron electron repulsion

integrals. We typically compute all three-electron integrals in large blocks, but for checking

purposes we also have codes which allow us to calculate the needed auxiliary integrals entirely

using Larsson sums. This slows down the calculation but verifies the correctness of the faster

techniques.

Acknowledgments

Maria Belen Ruiz checked many of the integrals in the table and pointed out the mistake in

equation (23) of the first paper in this series. William Mitchell and Peter Mohr have made

helpful comments on this work. One of us (JSS) would also like to thank Judith (Devaney)

Terrill for supporting this endeavour.

Appendix A. Condon and Shortley cks

A Condon and Shortley coefficient [10] is defined by

ck

(l, m;l



, m) =

 2

2k + 1  π

0

(km − m)(lm)(lm)sin θ dθ. (A.1)

This can be expressed equivalently in terms of the spherical harmonics Y m

l (θ , φ) as

ck

(l, m;l



, m) =

 4π

2k + 1 

Y m−m

k (θ , φ)Y m∗ l (θ , φ)Y m

l (θ , φ)sin θ dθ dφ. (A.2)

There is a substantial body of work on the evaluation of ck coefficients (actually on the

closely related Wigner 3j -symbols) at the double precision level of accuracy but very little

in the way of numerical results at the quadruple precision level of accuracy. Our approach

is to put the ck integral into the ‘Gaunt’ form5, since the Gaunt coefficient is equivalent to

5 [10, p 52, equation (18) and p 176, equation (11).]
Math and computational science issues in high precision Hy–CI calculations 1585

the ck coefficient except for a phase factor and is an efficient approach in terms of number of

arithmetic operations. As the l,l arguments increase; however, there can be a problem with

exponent and/or mantissa overflow. This has long been known, and several recursive methods

that circumvent the problem have been developed [16–18]. This work has been done at the

DP level and has been particularly effective for large values (>100) of l and l



. However, these

methods are not especially efficient and in this work we are not interested in large values of

the ck arguments.

The advantage of the Gaunt formulation, using ordinary factorials, is that the calculation

can be kept ‘exact’ for l,l 7 using DP arithmetic and for l,l 14 using QP arithmetic (with

attention to detail in the numerics). These limits are set by mantissa overflow. A significant

improvement in these l,l limits is possible if one uses a prime factor representation (PFR)

for the factorials. In this case, exact results are obtained for l,l 15 using DP and l,l 

30 using QP. The reason for this is that a large common factor can be trivially factored out of

each term in the Gaunt sum, thus delaying the point at which mantissa overflow occurs in the

summation. Once mantissa overflow occurs the substantial differencing inherent in the Gaunt

summation becomes a problem, even in the PFR case. For example, in the PFR case, for

l,l = 14 four digits are lost in taking the sum, while for l,l = 30 over 9 digits are lost.

Emphasis in this work has been on accuracy and not execution speed. However, since the

exact Gaunt formulation is a minimum operation one, the execution times turn out to be quite

acceptable, namely, approximately 0.25–0.30 μs per coefficient using DP and approximately

3 μs per coefficient using QP (these times are for a generic PC (1.90 GHz Pentium III)).

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

As explained in section 3.2, the case g < 0, f, h  0 can arise in the three-electron kinetic

energy integrals, whereas it cannot in the other three-electron interaction integrals treated in

the first paper in this series [1]. A suitable recursion scheme for this new W-function case

involves the V -function defined by

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

 ∞

0

xm e−αx dx

 ∞

x

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

The W-function is defined by

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

 ∞

0

xf e−αx dx

 ∞

x

yg e−βy dy

 ∞

y

zh e−γ z dz,

(f  0, f + g  −1, f + g + h  −2). (B.2)

For h = 0, i.e, in the f, g plane, the W’s are related to the V ’s in a very simple way, namely,

Wfg0(αβγ ) = γ −1

Vfg(α, β + γ ), (f  0, f + g  −1). (B.3)

Then we can use

Wfgh(αβγ ) = γ −1

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

(f  0, f + g  −1, f + g + h  −1, h  1) (B.4)

to raise h in a stable fashion. Note that we need Vf,g+h(α, β + γ ) for f = 0, fmax and

g + h = gmin, gmax + hmax. The accurate calculation of this array has been discussed in the first

paper in this series. The recurrence scheme is stable (the V ’s are all  0 and the recurrence

relation has no subtractions) and the integral accuracy will be as accurate as the calculation of

the V integrals.
1586 J S Sims and S A Hagstrom

Addenda

In equation (23) in the first paper in this series [1] we incorrectly stated the explicit finite sum

of Frolov and Smith [19] for the Wfgh functions for the special case of f, g, h  0. The

correct formula is

Wfgh(αβγ ) = 

h

σ=0

Cσ

h Aσ (γ )Vf,g+h−σ (α, β + γ ). (B.5)

In equation (7) in [2] we gave the spherical harmonic addition theorem to be

Pl(cos θ12) = 4π

(2n + 1)



l

m=−l

Y −m

l (1)Y m

l (2). (B.6)

That should have been

Pl(cos θ12) = 4π

(2l + 1)



l

m=−l

(−1)

mY −m

l (1)Y m

l (2). (B.7)

This can also be written as

Pl(cos θ12) = 4π

(2l + 1)



l

m=−l

Y m∗ l (1)Y m

l (2) (B.8)

since the Condon–Shortley phase condition is [10]

Y m∗ l (θ , φ) = (−1)

mY −m

l (θ , φ). (B.9)

References

[1] Sims J S and Hagstrom S A 2004 Mathematical and computational science issues in high precision Hylleraas–

configuration interaction variational calculations: I. Three-electron integrals J. Phys. B: At. Mol. Opt.

Phys. 37 1519–40

[2] 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

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

Stand. Technol. 107 223–45

[4] Sims J S and Hagstrom S A 2007 Mathematical and computational science issues in high precision Hylleraas–

configuration interaction variational calculations: III. Four-electron integrals J. Phys. B: At. Mol. Opt. Phys.

to be submitted

[5] 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

[6] 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

[7] 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

[8] Ruiz M B 2005 Hylleraas method for many-electron atoms: I. The Hamiltonian Int. J. Quantum Chem.

101 246–60

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

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

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

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

[12] Frolov A M 2004 Highly accurate calculation of the auxiliary functions of the fourth order and five-body

integrals J. Phys. B: At. Mol. Opt. Phys. 37 2103–18

[13] Kolos W and Roothaan C C J 1960 Accurate electronic wave functions for the H2 molecule Rev. Mod. Phys.

32 219–32

[14] Perkins J F 1968 Atomic integrals containing rλ

23r

μ

31rν

12 J. Chem. Phys. 48 1985–88
Math and computational science issues in high precision Hy–CI calculations 1587

[15] Sack R A 1964 Generalization of Laplace’s expansion to arbitrary powers and functions of the distance between

two points J. Math. Phys. 5 245–51

[16] Luscombe J H and Luban M 1998 Simplified recursive algorithm for Wigner 3j - and 6j -symbols Phys. Rev.

E 57 7274–7

[17] Schulten K and Gordon R G 1975 Exact recursive evaluation of 3j - and 6j -coefficients for quantum-mechanical

coupling of angular momenta J. Math. Phys. 16 1961–70

[18] Weniger E J and Steinborn E O 1982 Programs for the coupling of spherical harmonics Comput. Phys.

Commun. 25 149–57

[19] 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