PHYSICAL REVIEW A 83, 032518 (2011)

Hylleraas-configuration-interaction study of the 1S ground state of neutral beryllium

James S. Sims

National Institute of Standards and Technology, Gaithersburg, Maryland 20878-9957, USA

Stanley A. Hagstrom

Departments of Chemistry and Computer Science, Indiana University, Bloomington, Indiana 47405, USA

(Received 14 December 2010; published 31 March 2011)

Hylleraas-configuration-interaction (Hy-CI) method variational calculations are reported for the 1S ground state

of neutral beryllium. The best nonrelativistic energy obtained was −14.667 356 4 hartree, which is estimated to

be accurate to a tenth of a microhartree.

DOI: 10.1103/PhysRevA.83.032518 PACS number(s): 31.15.ve, 31.15.−p

I. INTRODUCTION

Beryllium, with its four electrons and strong mixing of the

1s2 2s2 and 1s2 2p2 configurations, has long been the subject

of investigation. The first reasonably accurate nonrelativistic

calculations on the ground state were the configurationinteraction (CI) studies of Watson [1] and Weiss [2]. Szasz

and Byrne [3] were the first to explicitly correlate the wave

function by adding intrashell rij coordinates in a Hylleraas

(Hy)-type wave function. Gentner and Burke [4] improved on

Szasz and Byrne’s work by including intershellrij coordinates.

Sims and Hagstrom [5] combined the Hy and CI methods in a

Hylleraas-configuration-interaction (Hy-CI) [5] calculation in

1971 with a considerable improvement in the energy. Bunge

did his first Be calculation in 1968 [6] and his CI calculations

for Be improved steadily over the years, culminating in a

benchmark calculation of over 2.6 million terms in 2010 [7]. A

multiconfiguration Hartree-Fock (MCHF) correlation study of

Be was presented by Froese Fischer and Saxena in 1974 [8] and

then was further refined by Froese Fischer in 1993 [9], leading

finally to a calculation of over 650 000 terms in 2010 [10].

While significant progress was being made in Hy treatments

of three-electron systems, integral problems with the conventional Hy expansions [11–13] brought calculations to a halt

as far as four-electron systems were concerned until finally, in

1998, Busse  ̈ et al. [14] overcame serious integral problems and

made major improvements through the use of Hy-type doubly

linked terms in the wave-function expansions. This was the

first major Hy-type calculation of near microhartree accuracy.

In 1995 the exponentially correlated Gaussian (ECG) method

was introduced by Komasa et al. [15]. The relative simplicity

of the resultant integrals led to a sequence of impressively

accurate calculations, culminating in the work of Adamowitz

and co-workers [16,17] which stands as the best work to date

(our work is second only to theirs). Table I summarizes Be

ground-state variational results listed in order of increasing

accuracy. In this paper, we present preliminary results from a

very large Hy-CI calculation (over 40 000 symmetry adapted

expansion functions) for the Be ground state.

II. METHOD OF CALCULATION

For four electrons, the Hy-CI wave function used is

Be = 

K

CKK, (1)

where

K = 



r

νK

ij 

4

s=1

{φKs (rs)}K



= OasOL,MLOS,MS



r

νK

ij 

4

s=1

{φKs (rs)}K



(2)

denotes the Kth antisymmetrized spin and angular momentum

projected configuration state function (CSF). OL,ML and

OS,MS are idempotent orbital and spin angular momentum

projection operators of the Lowdin type [  ̈ 24] for a state of

total quantum numbers L,ML,S,MS [Russell-Saunders (LS)

coupling is assumed]. In practice it is sufficient to take

νK equal to 0 or 1, with νK = 0 the CI case. K is a

primitive spin product function for the term K and φKs (rs)

represents the sth basis orbital in the Kth term. The basis

orbitals are taken to be un-normalized Slater-type orbitals.

Oas is the idempotent antisymmetry projection operator. For

four-electron singlet states there exist two linearly independent

primitive spin functions 1 = αβαβ and 2 = ααββ. It is

possible to converge on the exact wave function using only

the 1 product. Similar observations have been made by

Larsson [25] and by Sims and Hagstrom [26] for Li, and by

Sims and Hagstrom [5] for Be. Cencek and Rychlewski [27]

have given the general proof that only one primitive spin

function is needed to ensure convergence of eigenvalues to

the exact root of the Hamiltonian. The Appendix contains

a discussion of how to handle the antisymmetrization, spin,

and angular momentum projections involved in computing the

Hamiltonian and overlap matrix elements.

The CK coefficients in Eq. (1) are found by solving the

generalized eigenvalue problem HC = λSC, where HKL =

K|H|Land SKL = K|L, using the familiar inverse

iteration method. Root λ will be an upper bound to the exact

energy. The nonrelativistic Hamiltonian H used in this work is

H = 

4

i=1

Hi +

4

i<j

r−1

ij , (3)

where Hi = Ti + Vi is a one-electron operator (electron i)

consisting of a kinetic-energy part Ti = −1/2∇2

i and a nuclear

attraction part Vi = −Z/ri. Quadruple precision was used

throughout the calculations. As in our recent work [26,28–30],

MPI [31] was used to parallelize the code. The number of

1050-2947/2011/83(3)/032518(7) ©2011 American Physical Society 032518-1
JAMES S. SIMS AND STANLEY A. HAGSTROM PHYSICAL REVIEW A 83, 032518 (2011)

TABLE I. Comparison of theoretical Be ground-state nonrelativistic energies (in hartree). N is the number of terms in the expansion. All

energies are variational except those labeled “Estimated exact.”

Technique Author(s) N Energy E (hartree)

Hy Szasz and Byrne (1967) [3] 28 −14.6565

CI Watson (1960) [1] 37 −14.657 40

Hy Gentner and Burke (1968) [4] 25 −14.6579

CI Weiss (1961) [2] 55 −14.660 90

CI Bunge (1968) [6] 180 −14.664 19

MCHF Froese Fischer and Saxena (1974) [8] 52 −14.665 87

Hy-CI Sims and Hagstrom (1971) [5] 107 −14.666 547

CI Bunge (1976) [18] 650 −14.666 902

MCHF Froese Fischer (1993) [9] 3381 −14.667 113

MCHF Froese Fischer (2010) [10] 652 683 −14.667 114 52

CI Jitrik and Bunge (1997) [19] 212 710 −14.667 275 57

CI Bunge (2010) [7] 2 614 689 −14.667 347 30

Hy-CI this work, s, p rij basis 20 330 −14.667 350 33

Hy Busse  ̈ et al. (1998) [14] 5306 −14.667 354 7

ECG Komasa et al. (1995) [15] 1200 −14.667 355 021

ECG Komasa (2002) [20] 3700 −14.667 355 627

ECG Pachucki and Komasa (2006) [21] 4600 −14.667 355 748

Hy-CI this work, s, p, d rij basis 40 784 −14.667 356 359

Hy-CI this work, s, p, d, f rij basis 41 871 −14.667 356 411

ECG Stanke et al. (2007) [16] 4600 −14.667 356 458

ECG Stanke et al. (2009) [17] 10 000 −14.667 356 486

Estimated exact Bunge (2010) [7] −14.667 355(1)

Estimated exact Pachucki and Komasa (2004) [22] −14.667 355 7(1)

Estimated exact Davidson et al. (1991) [23] −14.667 36

Estimated exact Komasa et al. (1995) [15] −14.667 360(2)

processes used varied from 12 to 48 with 32 being typical.

This is not large by current standards and scaling to 128–256

processes with the existing code would be feasible on a routine

basis provided one can be assured of adequate resources

per process. The parallel environment for this work was the

National Institute of Standards and Technology (NIST)’s 394

processor Linux cluster.

For Be the generation of the matrices H and S takes much

longer than the solution of the secular equation due to the very

large number of four-electron integrals. For large expansion

lengths N the ratio of matrix element build to eigenvalue

solve phases of the calculation was typically about 20 to

1 using an efficient solver (inverse iteration) but relatively

inefficient integral and matrix element build packages. Inverse

iteration scales as expected but matrix element build scales

rather poorly (we do not really know why). Four-electron

integrals totally dominate the integral calculation, 98 percent

of the total integral time being typical. Improvements to the

code will come in this area. Details of how we distribute the

workload, allocate memory, manage mass storage, and deal

with serial aspects of the algorithm will be described in a

future publication.

III. RESULTS

The CSFs given by Eq. (2) can be written

K =  [FK(r1,r2,r3,r4)K] (4)

in terms of spatial and spin functions FK(r1,r2,r3,r4) and K.

Only one spin function K = 1 = αβαβ is used in this work

and the spatial part of a CSF is given by a particular choice of

Hartree orbital product and rij factor:

FK(r1,r2,r3,r4) = r

νK

ij 

4

s=1

{φKs(rs)}, (5)

where νK is either 0 or 1.

Table II lists our best energy results (column 4) for various

expansion lengths Ntot shown in column 3. Column 2 gives

the number of CSFs added for the block shown in column 1.

The energy improvements for each CSF block type are shown

in column 5. Column 1 lists the basis orbitals that are used

to generate the CSFs for each block type in the order electron

1 (α spin), electron 2 (β spin), electron 3 (α spin), electron 4

(β spin). For example, in the first line 1:8sK for the first electron

means the basis orbitals are 1sK through 8sK orbitals (where

K indicates an orbital exponent appropriate for a K-shell

electron). Products of four orbitals are built up by taking one

pair of orbitals from the K-shell set and the other pair from the

L-shell set. All the orbital promotions are within the shell, not

between shells. This leaves out a substantial number of CSFs

that turn out to be of no importance in this work. The choice of

terms is highly regular, there having been no serious attempt

to eliminate energetically unimportant terms. The number of

CSFs in a block can be computed from the listed basis orbitals

and the condition that Nrsum, the sum of the powers of r for

each Hartree product (HP), must be less than or equal to 16

(i.e., Nrsum  16). For example, consider 2:8pKp 2:8pKp 1:8sL

1:8sL r12. There are (7 × 8)/2 = 28 unique pairs of orbitals

for electrons 1 and 2, and (8 × 9)/2 = 36 unique pairs of

032518-2
HYLLERAAS-CONFIGURATION-INTERACTION STUDY OF ... PHYSICAL REVIEW A 83, 032518 (2011)

orbitals for electrons 3 and 4. Since the K and L shells have

no orbital in common, there will be 28 × 36 = 1008 different

CSF terms for this block. Applying the condition that the sum

of the powers of r for these orbitals has to be 16, the number

of terms is reduced to 630.

The leading CSF term types in Table II are what one

would expect from adding explicit correlation to a typical

Be CI expansion. For some of the remaining term types,

the Legendre expansion of typical Hy rij product terms

helped in determining possible Hy-CI term types to try [26].

For example, expansion of s1s2s3s4r2

23r234 (only even power

products need to be considered) gives rise to a spdp term,

in addition to ssss, spsp, sspp, and spps terms. Similarly,

r2

12r223 gives rise to pdps and r213r223 to ppds. These are not

exactly obvious important CI promotions. Of course, other

orbital products arise from the appropriate Hy terms, e.g., ppss

from r2

12, spsp from r224, and so on. From the possible r2

ij r2kl

products, there are six linearly independent orbital products

used in this work, namely, sppd, ppsd, sdpp, pspd, dpsp,

and pdsp.

Column 6 in Table II lists ω, the sum of the powers of the rij

products for the Hy equivalent of the Hy-CI term. Except for

the sK1 blocks toward the end of Table II, the blocks are ordered

by increasing ω. Except for ppdd, ddpp, and pppp, all CSF

types give rise to only one 1S state obtained from projecting

the orbital product (l10, l20, l30, l40). For ppdd, ddpp, and

pppp, three 1S states are obtained by projection of (0, 0, 0, 0),

(1, 1, −1, −1) and (1, −1, 1, −1), where only the m quantum

numbers are shown. For ddpp and pppp, all three projections

were used. For ppdd only the first two cases proved to be

energetically important. CSF blocks were tested individually

for importance essentially by trial and error, keeping in mind

that the contribution depends on the order in which blocks

are added. All blocks which lowered the energy by more than

5 nanohartree were kept.

Final orbital exponent optimization was done at the N =

19 076 expansion level for s and p orbitals and at the N =

27 468 level for d orbitals. f orbital exponents were taken

over from an earlier optimization. Due to limitations in the

integral codes, it was not possible to give every orbital type

its own orbital exponent. In a practical sense this was not a

problem, since the dependence of the final energy surface on

the orbital exponents turned out to be extremely flat for the

5-nanohartree cutoff used. Careful minimization of the ζ s was

not a particular problem for 0.1 microhartree accuracy.

IV. DISCUSSION

Table I summarizes earlier Be ground-state variational

results listed in order of increasing accuracy and includes

three important limiting energy values from the present work

inserted to facilitate comparison with the earlier work. Note

that 6-microhartree accuracy is obtained with just an s,p rij

basis, which is better than the best CI result of Bunge [7]

using over 2.6 million CSFs with l up to 30. A 40 784 term

s,p,d rij basis wave function (all the blocks in Table II

except for the three blocks containing f orbitals) results in

better than microhartree accuracy (0.2 microhartree). Finally

adding f orbitals gives an s,p,d,f rij energy of E(41 871) =

−14.667 356 411 hartree. Use of an s,p,d,f rij basis is

consistent with what we know from our earlier Li calculations

[26]. However, it should be noted that the ECG result of Stanke

et al. [17] is fully 75 nanohartree below our current best. From

the various different estimates of the exact energy it should

be clear that there is no doubt about the exact energy at the

microhartree level, but beyond that it gets problematic.

In Table II are listed results for the Be wave function

expansions at various levels of truncation. The most distinctive

feature of the table is the obvious importance of CSF blocks

with r12 and r34 reflecting the pronounced Be shell structure.

There are a few important blocks incorporating intershell

correlation (via r13- and r14-type terms), but there are also

intershell blocks with anomalously small contributions, for

which there is no ready explanation. There are 18 blocks

that contribute at the microhartree level out of a total of 71

given in the table. Many more block types (approximately 78)

were tested but not included in the final wave function due to

the 5-nanohartree cutoff. The spdf {1,r12,r34,r13,r14} blocks

are an example. Many blocks were not tried; many of the

ones tried and dropped because of the 5-nanohartree cutoff

could contribute at the 1-nanohartree level. Some blocks which

lowered the energy by less than 5 nanohartree were kept

because they function as “sanity” checks on the calculation.

One example is the 4:8fKf 4:8fKf 1:7sL1:7sLr34 block which

is involved in our attempt to include an r12r34-type correlation,

as discussed below. Another example is that the last two

blocks in the table are part of the sK1 set of expansion

terms, which introduces additional correlation into the K

shell. After we arrived at the 41 871 term wave function,

further tests were done to check the choice of blocks in the

wave function. For example, blocks containing terms in which

K- and L-shell orbitals are not both doubly occupied and

containing either r23 or r24 were added to the 41 871 term

wave function, but none of these blocks added more than 2

nanohartree. As another example, our algorithm for building

HPs omits HPs with {msK,m

sK}{nsL,nsL} m m,n n.

These terms were also tested but not included in the final wave

function because they led to only very small improvements;

they were not important for Li [26] either.

The Hy calculation by Busse  ̈ et al. [14] in Table I is the calculation closest in spirit and results to the current calculation. A

direct comparison is difficult, however, because of differences

in choice of spin function as well as, of course, the presence of

products of rij factors raised to odd powers (so-called odd-odd

terms, e.g., rij rik and rij rkl), which appear in Hy expansions

but are not explicitly represented in the corresponding Hy-CI

expansions. For Li [26] it has been shown that odd power

linked product terms (rij products with one index in common)

are unimportant at the nanohartree level, and by inference

probably not a problem at the four-electron level, although

there are more such terms in the Be case, and one might expect

their contribution to be somewhat greater. The unlinked rij rkl

(no indices in common) term types, which first occur in the

four-electron case, are expected to be of major importance,

however. Support for this can be inferred from Busse  ̈ et al.’s

[14] ω = 2 results. There is a big 2000-microhartree improvement due to odd power term types of which r12r34 would

be expected to be the most important. In Hy-CI this effect

can hopefully be effectively represented using a superposition

of normal Hy-CI term types. In particular we suggest that

032518-3
JAMES S. SIMS AND STANLEY A. HAGSTROM PHYSICAL REVIEW A 83, 032518 (2011)

TABLE II. Hy-CI results for the Be ground-state energy (in hartree). In the table, N is the number of CSF terms added with each CSF

block, Ntot is the cumulative number of terms, and ω is the sum of the rij powers in an equivalent Hy expansion. All terms are r-sum filtered;

Nrsum = 16.

Terms addeda,b N Ntot E(Ntot) (hartree) −E (microhartree) ω

1:8sK 1:8sK 1:8sL 1:8sL 896 896 −14.591 724 231 0

1:8sK 1:8sK 1:8sL 1:8sL r12 896 1792 −14.659 886 529 68 162.298 0 + 1

1:8sK 1:8sK 1:8sL 1:8sL r34 896 2688 −14.662 973 074 3040.778 0 + 1

1:8sK 1:8sK 1:8sL 1:8sL r13 896 3584 −14.662 930 798 3.491 0 + 1

1:8sK 1:8sK 1:8sL 1:8sL r14 896 4480 −14.662 934 922 4.124 0 + 1

2:8pKp 2:8pKp 1:8sL 1:8sL 630 5110 −14.664 420 851 1485.930 2

2:8pKp 2:8pKp 1:8sL 1:8sL r12 630 5740 −14.664 504 277 83.426 2 + 1

2:8pKp 2:8pKp 1:8sL 1:8sL r34 630 6370 −14.666 606 396 2102.119 2 + 1

2:8pKp 2:8pKp 1:8sL 1:8sL r13 630 7000 −14.666 816 908 210.512 2 + 1

2:8pKp 2:8pKp 1:8sL 1:8sL r14 630 7630 −14.666 947 723 130.816 2 + 1

1:8sK 1:8sK 2:8pLp 2:8pLp 630 8260 −14.667 138 747 191.023 2

1:8sK 1:8sK 2:8pLp 2:8pLp r12 630 8890 −14.667 289 253 150.506 2 + 1

1:8sK 1:8sK 2:8pLp 2:8pLp r34 630 9520 −14.667 289 795 0.543 2 + 1

1:8sK 1:8sK 2:8pLp 2:8pLp r13 630 10 150 −14.667 306 512 16.717 2 + 1

1:8sK 1:8sK 2:8pLp 2:8pLp r14 630 10 780 −14.667 314 875 8.363 2 + 1

1:6sK 2:6pKp 1:6sL 2:6pLp 865 11 645 −14.667 314 971 0.096 2

1:6sK 2:6pKp 1:6sL 2:6pLp r12 865 12 510 −14.667 316 045 1.074 2 + 1

1:6sK 2:6pKp 1:6sL 2:6pLp r34 865 13 375 −14.667 316 303 0.258 2 + 1

1:6sK 2:6pKp 1:6sL 2:6pLp r13 865 14 240 −14.667 316 355 0.051 2 + 1

1:6sK 2:6pKp 2:6pLp 1:6sL 865 15 105 −14.667 316 402 0.048 2

1:6sK 2:6pKp 2:6pLp 1:6sL r12 865 15 970 −14.667 316 503 0.100 2 + 1

1:6sK 2:6pKp 2:6pLp 1:6sL r14 865 16 835 −14.667 316 618 0.115 2 + 1

2:7pKp 2:7pKp 2:7pLp 2:7pLp 987 17 822 −14.667 348 007 31.389 4

2:6pKp 2:6pKp 2:6pLp 2:6pLp r12 627 18 449 −14.667 350 243 2.236 4 + 1

2:6pKp 2:6pKp 2:6pLp 2:6pLp r34 627 19 076 −14.667 350 287 0.044 4 + 1

2:6pKp 2:6pKp 2:6pLp 2:6pLp r13 627 19 703 −14.667 350 326 0.040 4 + 1

2:6pKp 2:6pKp 2:6pLp 2:6pLp r14 627 20 330 −14.667 350 336 0.010 4 + 1

1:6sK 2:5pKp 2:5pLp 3:6dLd 379 20 709 −14.667 352 968 2.632 4

1:6sK 2:6pKp 2:5pLp 3:6dLd r12 465 21 174 −14.667 353 413 0.444 4 + 1

1:6sK 2:5pKp 2:5pLp 3:6dLd r34 379 21 553 −14.667 353 496 0.084 4 + 1

1:6sK 2:5pKp 2:5pLp 3:6dLd r13 379 21 932 −14.667 353 568 0.071 4 + 1

1:6sK 2:5pKp 2:5pLp 3:6dLd r14 379 22 311 −14.667 353 644 0.076 4 + 1

2:6pKp 2:6pKp 1:7sL 3:6dLd 375 22 686 −14.667 353 689 0.045 4

2:6pKp 2:6pKp 1:7sL 3:6dLd r12 375 23 061 −14.667 353 803 0.114 4 + 1

2:6pKp 2:6pKp 1:7sL 3:6dLd r14 375 23 436 −14.667 353 822 0.019 4 + 1

2:6pKp 2:6pKp 1:7sL 3:6dLd r13 375 23 811 −14.667 353 832 0.010 4 + 1

2:6pKp 2:6pKp 1:7sL 3:6dLd r34 375 24 186 −14.667 353 837 0.005 4 + 1

1:6sK 3:7dKd 2:7pLp 2:7pLp 506 24 692 −14.667 353 992 0.155 4 + 1

1:6sK 3:7dKd 2:7pLp 2:7pLp r12 506 25 198 −14.667 353 999 0.007 4 + 1

1:6sK 3:7dKd 2:7pLp 2:7pLp r13 506 25 704 −14.667 354 005 0.006 4 + 1

1:6sK 3:7dKd 2:7pLp 2:7pLp r34 506 26210 −14.667 354 042 0.037 4 + 1

2:6pKp 1:6sK 2:7pLp 3:7dKd 776 26 986 −14.667 354 049 0.007 4 + 1

2:6pKp 1:6sK 2:7pLp 3:7dKd r14 776 27 762 −14.667 354 154 0.104 4 + 1

2:6pKp 1:6sK 2:7pLp 3:7dKd r13 776 28 538 −14.667 354 169 0.015 4 + 1

3:7dKd 2:6pkp 2:7sL 2:7pLp 702 29 240 −14.667 354 173 0.005 4 + 1

2:6pKp 3:7dKd 1:7sL 2:7pLp r12 851 30 091 −14.667 354 180 0.007 4 + 1

2:6pKp 3:7dKd 1:7sL 2:7pLp r14 851 30 942 −14.667 354 187 0.007 4 + 1

2:6pKp 3:7dKd 1:7sL 2:7pLp 851 31 793 −14.667 354 194 0.007 4 + 1

1:6sK 1:6sK 3:7dLd 3:7dLd 266 32 059 −14.667 354 437 0.243 4 + 1

1:6sK 1:6sK 3:7dLd 3:7dLd r12 266 32 325 −14.667 355 886 1.449 4 + 1

1:6sK 1:6sK 3:7dLd 3:7dLd r34 266 32 591 −14.667 355 913 0.027 4 + 1

1:6sK 1:6sK 3:7dLd 3:7dLd r13 266 32 857 −14.667 355 941 0.028 4 + 1

1:6sK 1:6sK 3:7dLd 3:7dLd r14 266 33 123 −14.667 355 956 0.015 4 + 1

3:8dKd 3:8dKd 1:7sL 1:7sL 378 33 501 −14.667 356 065 0.110 4 + 1

3:8dKd 3:8dKd 1:7sL 1:7sL r34 378 33 879 −14.667 356 107 0.042 4 + 1

3:8dKd 3:8dKd 1:7sL 1:7sL r12 378 34 257 −14.667 356 122 0.015 4 + 1

3:8dKd 3:8dKd 1:7sL 1:7sL r13 378 34 635 −14.667 356 130 0.008 4 + 1

032518-4
HYLLERAAS-CONFIGURATION-INTERACTION STUDY OF ... PHYSICAL REVIEW A 83, 032518 (2011)

TABLE II. Continued.

Terms addeda,b N Ntot E(Ntot) (hartree) −E (microhartree) ω

3:8dKd 3:8dKd 1:7sL 1:7sL r14 378 35 013 −14.667 356 135 0.005 4 + 1

1:6sK 3:7dKd 3:7dLd 1:7sL 851 35 864 −14.667 356 140 0.005 4

1:6sK 3:7dKd 3:7dLd 1:7sL r14 851 36 715 −14.667 356 146 0.006 4 + 1

3:8dKd 3:8dKd 2:7pLp 2:7pLp 732 37 447 −14.667 356 179 0.034 6

3:8dKd 3:8dKd 2:7pLp 2:7pLp r34 732 38 179 −14.667 356 302 0.123 6 + 1

3:8dKd 3:8dKd 2:7pLp 2:7pLp r12 732 38 911 −14.667 356 308 0.006 6 + 1

2:7pKp 2:7pKp 3:8dLd 3:8dLd 488 39 399 −14.667 356 320 0.012 6

2:7pKp 2:7pKp 3:8dLd 3:8dLd r12 488 39 887 −14.667 356 324 0.004 6 + 1

1:6sK 1:6sK 4:8fLf 4:8fLf r12 210 40 097 −14.667 356 371 0.047 6 + 1

4:8fKf 4:8fKf 1:6sL 1:6sL r34 232 40 329 −14.667 356 374 0.003 6 + 1

1:6sK 4:8fKf 4:8fLf 1:6sL r34 645 40 974 −14.667 356 378 0.004 6

1:6sK1 1:6sK1 2:6sL 2:6sL 299 41 273 −14.667 356 404 0.026 0

1:6sK1 1:6sK1 2:6sL 2:6sL r12 299 41 572 −14.667 356 407 0.003 0 + 1

1:6sK1 1:6sK1 2:6sL 2:6sL r34 299 41 871 −14.667 356 411 0.004 0 + 1

a

K-shell orbital exponents are K = 3.3, K1 = 10.0, Kp = 4.65, Kd = 4.65. Kf = 4.5 bL-shell orbital exponents are L = 1.6, and Lp = 2.0, Ld = 2.6, Lf = 3.0.

the forms (ss + pp + dd +···)K(ss + pp + dd +···)L

r34 + (ss + pp + dd +···)Kr12(ss + pp + dd +···)L will

accomplish this. In this representation one electron pair is

being correlated primarily by the rij factor while the other

pair is correlated by a CI pair expansion of the sort used

for He [28], which is unfortunately slowly converging in the

CI case, but may be better behaved in Be. In Table III are

gathered together results from Table II, with additional test

results using (ss)Kr12(gg)L and (ff )K(pp)Lr34 term types

added specifically to explore this point, suggesting that this is

indeed the case. For example, the series (ss)Kr12(pp + dd +

ff )L converges nicely as does (pp + dd + ff )K(ss)Lr34.

(pp)Kr12(pp + dd + ff )L and (pp + dd + ff )K(pp)Lr34

contribute similarly. There is one anomalous case for which

we have no explanation. Note also that the numbers in Table III

depend on the block order used in Table II, an ordering that was

not explicitly designed to show the r12r34 convergence effect.

Compared to the slow, cusp-connected convergences in typical

CI calculations, this is unusually fast convergence, suggesting

that this correlation type can be accurately represented within

the Hy-CI model.

However, further research into this point is warranted. In

this connection, an Hy treatment of (hopefully) nanohartree

accuracy both with and without unlinked cluster terms, analogous to the way the effect of linked cluster terms like r12r13 was

examined in Li [26], would be most useful, especially since the

issue of unlinked cluster terms is qualitatively different from

the issue of linked cluster terms. Generalization of Hy-CI to

include unlinked rij rkl terms would certainly settle the issue,

but at the cost of greatly complicating the calculation of matrix

elements.

A 5-nanohartree threshold was used to determine which

term types to include in the final wave function. Blocks

which contribute less than this threshold value were dropped.

Tests to date suggest that (20–25) nanohartree could probably

be picked up using a 1-nanohartree cutoff leaving about

50 nanohartree to reach or surpass the Stanke et al. [17] value,

but to do substantially better, assuming that r12r34 is not a

problem, will involve addressing issues of (i) greater flexibility

in the atomic orbital basis (the number and choice of nonlinear

parameters is very important at the nanohartree level for

Li [26], and there is no reason to believe this will not be the case

for Be as well); (ii) more careful optimization of the nonlinear

parameters; (iii) focusing on the most important (core) CSF

types involving s and p orbitals; and (iv) better CSF filtering

to reduce the expansion lengths. Extensive experimentation

TABLE III. Convergence of r12 r34 as represented by Hy-CI.

Block Energy E (nanohartree) Block Energy E (nanohartree)

(ss)K r12(ss)L 68 162 298 (ss)K (ss)Lr34 3 040 778

(ss)K r12(pp)L 151 506 (pp)K (ss)Lr34 2 102 114

(ss)K r12(dd)L 1449 (dd)K (ss)Lr34 42

(ss)K r12(ff )L 47 (ff )K (ss)Lr34 3

(ss)K r12(gg)L 0

(pp)K r12(ss)L 83 426 (ss)K (pp)Lr34 543

(pp)K r12(pp)L 2236 (pp)K (pp)Lr34 44

(pp)K r12(dd)L 4 (dd)K (pp)Lr34 123

(ff )K (pp)Lr34 3

(dd)K r12(ss)L 15 (ss)K (dd)Lr34 27

(dd)K r12(pp)L 6 (pp)K (dd)Lr34 4

032518-5
JAMES S. SIMS AND STANLEY A. HAGSTROM PHYSICAL REVIEW A 83, 032518 (2011)

will be needed to find the best combination of CSFs and

orbital parameters upon which to base a much larger so-called

“full, complete expansion.” It should be noted that the present

calculations are very computationally expensive with the H

and S matrix element build step taking an order of magnitude

longer than the eigenvalue solve step. This is due to the

complete dominance of the four-electron integral computation

time (by a factor of 100 or more) over the time required for

all other integrals. Fortunately, substantial improvements are

possible in both the four-electron integral package and in

the eigenvalue routines currently being used. Ideally these

improved codes will lead to more efficient computational

access to not only the ground state of Be, but to excited states

(both S and non-S, singlet and triplet) as well as Be-like ions.

V. CONCLUSION

The ability of Hy-CI calculations to achieve 0.1-

microhartree accuracy for beryllium is shown in this work.

The convergence of r12r34 term types has been investigated,

and Hy-CI appears to accurately represent this term type,

although more research is needed on this point. Given the

close relationship between Hy-CI and Hy calculations when

the Hy expansion terms contain at most a single odd power

of rij , it should be possible using Hy expansions to achieve

comparable accuracy for beryllium without using the odd-odd

power products of rij , which leads to severe integral problems

in Hy calculations for four or more electrons. Thus even though

the present results are preliminary, they significantly advance

prospects for accurate Hy-CI calculations on four or more

electron atoms1.For Hy-CI calculations to go beyond the 0.1-

microhartree accuracy achieved in this work, further research

is needed on r12r34 representation, more flexible atomic orbital

basis sets, and better CSF filtering techniques to control

expansion lengths.

ACKNOWLEDGMENTS

We are indebted to Frederick W. King for keeping us up

to date on rij calculations on two-, three-, and four-electron

systems, and to Carlos Bunge for doing the same with respect

to accurate atomic CI calculations (and for carrying out the

benchmark CI calculations in Ref. [7] at our request). We also

thank Arne Luchow, Krzysztof Pachucki, Ludwik Adamowicz,  ̈

and Maria Belen Ruiz for correspondence concerning their

calculations, and Carl Spangler, Michael Strawbridge, and

Denis Lehane (NIST) for parallel system support. One of

us (J.S.S.) would like to thank Judith Devaney Terrill for

supporting this endeavor.

APPENDIX: PROJECTION REDUCTION

In the generalized eigenvalue problem HC = λSC, matrix

element HKL is

HKL = K|H|L= 

P

K



H

P

L



, (A1)

1An extension of Hy-CI to the five-electron boron case is underway

[32].

where P

K denotes the Kth primitive (unprojected) function

r

νK

ij

4

s=1{φKs (rs)}K. SKL is similar, with H replaced by the

unit operator 1.

Now notice that

N!Oas = 

p

(−1)pP

=





u

(−1)uPα

u

 

ν

(−1)νPβ

ν



+

γ

(−1)γ Pαβ

γ

= AαAβ + Aαβ , (A2)

where N is the number of electrons, Pα refers to those

permutations which permute only α spins, Pβ refers to those

permutations which permute only β spins, and Pαβ refers to

only those permutations which permute an α spin and a β spin.

If we rewrite the bra and ket functions in Eq. (A1) in terms

of a single spin function, the HKL matrix element will end

up involving only AαAβ since permutations between α and β

give zero results when integrating over spin.

In our case K = L = 1 = α(1)β(2)α(3)β(4) =

αβαβ, so we can write

P

K = OL,MLOasOS,MS gKfK1, (A3)

where gK = r

νK

ij and fK is the Hartree product 4

s=1φKs (rs).

Using the quantum-mechanical “turnover rule” [33], the

commutativity of  and H, and the idempotency condition

† = , HKL reduces to

HKL = 

P

K



H

P

L



. (A4)

Projecting on 1 with the Lowdin spin projection operator [  ̈ 24]

OS,MS (S = 0,MS = 0), we get

OS,MS1 = OS,MSαβαβ

= 1

3αβαβ − 16 (ααββ + αββα + βααβ + ββαα)

+ 1

3βαβα, (A5)

where the subscripts (not shown) on the α and β are always in

the order 1,2,3,4.

This can be written

OS,MS1 = 	 1

3 − 16

	

Pσ

23 + Pσ34 + Pσ12 + Pσ14


+ 1

3Pσ

12P34
σ

1,

(A6)

where Pσ is used to denote permutations over spin coordinates

(Pr will be used later to denote permutations over spatial

coordinates). Using the identity

Oas = (−1)pOas(Pr

Pσ )

−1

, (A7)

HKL becomes

HKL = 

fKgK1|H|OL,MLOas

× 	 1

3 + 16

	

Pr

23 + Pr34 + Pr12 + Pr14


+ 1

3Pr

12P34
r

× gLfL1



. (A8)

032518-6
HYLLERAAS-CONFIGURATION-INTERACTION STUDY OF ... PHYSICAL REVIEW A 83, 032518 (2011)

Now we can integrate over spin, effectively removing spin

from the matrix element:

HKL = gKfK|H|AαAβOL,ML Br(gLfL). (A9)

We have used the fact that OL,ML commutes with Oas in the

above equation. Here the “reduced” antisymmetrizer AαAβ ,

AαAβ = (1 − P13)(1 − P24), (A10)

and

B = 1

3 + 16 (P23 + P34 + P12 + P14) + 13P12P34 (A11)

operate on the spatial coordinates. r has been dropped from

the terms in B since only spatial coordinates remain.

Using Eq. (A11), we obtain for our final expression for HKL

HKL = gKfK|H|(1 − P13)(1 − P24)OL,ML B(gLfL). (A12)

In Eq. (A12) one can apply OL,ML either before or after

applying B. We first apply the permutations B to gLfL,

then we project on the resulting terms with OL,ML . The

OL,ML projection is applied only on the orbital products since

gL = r

νL

ij commutes with OL,ML [5,34].

Depending on the structure of gLfL, (1 − P13)(1 − P24)

B(gLfL) may further simplify, and this should be checked

before application of OL,ML . In practice we routinely swap

the bras and kets if doing so will reduce the number of terms

on the right-hand side in Eq. (A12) or otherwise simplify the

process of assembling matrix elements.

[1] R. E. Watson, Phys. Rev. 119, 170 (1960).

[2] A. W. Weiss, Phys. Rev. 122, 1826 (1961).

[3] L. Szasz and J. Byrne, Phys. Rev. 158, 34 (1967).

[4] R. F. Gentner and E. A. Burke, Phys. Rev. 176, 63 (1968).

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

[6] C. Bunge, Phys. Rev. 168, 92 (1968).

[7] C. F. Bunge, Theor. Chem. Acc. 126, 139 (2010).

[8] C. Froese Fischer and K. M. S. Saxena, Phys. Rev. A 9, 1498

(1974).

[9] C. Froese Fischer, J. Phys. B 26, 855 (1993).

[10] C. Froese Fischer, J. Phys. B 43, 074017 (2010).

[11] M. Puchalski and K. Pachucki, Phys. Rev. A 73, 022503 (2006).

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

[13] F. E. Harris, A. M. Frolov, and V. H. Smith Jr., J. Chem. Phys.

119, 8823 (2003).

[14] G. Busse, H. Kleindienst, and A. L  ̈ uchow,  ̈ Int. J. Quantum Chem.

66, 241 (1998).

[15] J. Komasa, W. Cencek, and J. Rychlewski, Phys. Rev. A 52,

4500 (1995).

[16] M. Stanke, D. Kedziera, S. Bubin, and L. Adamowicz, Phys.

Rev. Lett. 99, 043001 (2007).

[17] M. Stanke, J. Komasa, S. Bubin, and L. Adamowicz, Phys. Rev.

A 80, 022514 (2009).

[18] C. Bunge, Phys. Rev. A 14, 1965 (1976).

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

[20] J. Komasa, Chem. Phys. Lett. 363, 307 (2002).

[21] K. Pachucki and J. Komasa, Phys. Rev. A 73, 052502 (2006).

[22] K. Pachucki and J. Komasa, Phys. Rev. Lett. 92, 213001 (2004).

[23] E. R. Davidson, S. A. Hagstrom, S. J. Chakravorty, V. M. Umar,

and C. Froese Fischer, Phys. Rev. A 44, 7071 (1991).

[24] P. O. Lowdin,  ̈ Rev. Mod. Phys. 36, 966 (1964).

[25] S. Larsson, Phys. Rev. 169, 49 (1968).

[26] J. S. Sims and S. A. Hagstrom, Phys. Rev. A 80, 052507

(2009).

[27] W. Cencek and J. Rychlewski, J. Chem. Phys. 98, 1252 (1993).

[28] J. S. Sims and S. A. Hagstrom, Int. J. Quantum Chem. 90, 1600

(2002).

[29] J. S. Sims et al., J. Res. Natl. Inst. Stand. Technol. 107, 223

(2002).

[30] J. S. Sims and S. A. Hagstrom, J. Chem. Phys. 124, 094101

(2006).

[31] Message Passing Interface Forum, Int. J. Supercomput. Appl. 8,

159 (1994).

[32] M. B. Ruiz, Habilitation thesis, University of Erlangen,

Erlangen, Germany, 2010.

[33] F. L. Pilar, Elementary Quantum Chemistry, 2nd ed.

(McGraw-Hill, New York, 1980).

[34] J. S. Sims, Ph.D. dissertation, Bloomington, Indiana, 1969.

032518-7