High precision variational calculations for the Born-Oppenheimer energies

of the ground state of the hydrogen molecule

James S. Simsa

National Institute of Standards and Technology, Gaithersburg, Maryland 20899-8911

Stanley A. Hagstrom

Department of Chemistry, Indiana University, Bloomington, Indiana 47405

Received 23 November 2005; accepted 17 January 2006; published online 1 March 2006

Born-Oppenheimer approximation Hylleraas variational calculations with up to 7034 expansion

terms are reported for the 1

g

+ ground state of neutral hydrogen at various internuclear distances. The

nonrelativistic energy is calculated to be −1.174 475 714 2201hartree at R=1.4 bohr, which is

four orders of magnitude better than the best previous Hylleraas calculation, that of Wolniewicz

J. Chem. Phys. 103, 1792 1995. This result agrees well with the best previous variational energy,

−1.174 475 714 216 hartree, of Cencek personal communication, obtained using explicitly

correlated Gaussians ECGsCencek and Rychlewski, J. Chem. Phys. 98, 1252 1993; Cencek

et al., ibid. 95, 2572 1995; Rychlewski, Adv. Quantum Chem. 31, 173 1998. The uncertainty

in our result is also discussed. The nonrelativistic energy is calculated to be

−1.174 475 931 3991hartree at the equilibrium R=1.4011 bohr distance. This result also agrees

well with the best previous variational energy, −1.174 475 931 389 hartree, of Cencek and

Rychlewski Rychlewski, Handbook of Molecular Physics and Quantum Chemistry, edited by

S. Wilson Wiley, New York, 2003, Vol. 2, pp. 199–218; Rychlewski, Explicitly Correlated Wave

Functions in Chemistry and Physics Theory and Applications, edited by J. Rychlewski Kluwer

Academic, Dordrecht, 2003, pp. 91–147., obtained using ECGs. © 2006 American Institute of

Physics. DOI: 10.1063/1.2173250

I. INTRODUCTION

Variational methods based on explicitly correlated wave

functions are known to give the most accurate upper bounds

to energy states, and hence the inclusion of terms containing

the interelectronic distance rij in the wave function has become increasingly common, at least for few-electron atomic

systems N4so common, in fact, that a book dealing

entirely with explicitly correlated functions has recently been

produced1

. The milestone in the theory of the hydrogen

molecule, the simplest molecular system containing an electron pair bond, is the work of James and Coolidge.2 Following the work of Hylleraas3 on the helium atom, they employed factors of r12 in the hydrogen molecule wave function

the full bibliography on H2 calculations until 1960 is given

in Ref. 4 and established beyond a doubt the usefulness of

including the interelectronic distance explicitly in the wave

function. Kolos and co-worker5,6 and Wolniewicz7 generalized the approach of James and Coolidge to get a more accurate description of dissociation. Wave functions using their

“generalized James-Coolidge JC” wave functions are commonly referred to as Hylleraas Hyor more specifically

Kolos-Wolniewicz KWwave functions.1 In addition to Hy

wave function calculations, the Hylleraas-configuration interaction Hy-CItechnique developed by us8 and also independently by Woźnicki9

has been applied to diatomic molecules including H2by Clary and co-worker10–13 and

Clementi and co-workers.14–19 In Hy-CI calculations the

wave function is expanded as a linear combination of correlated configuration state functions CSFs, where the unique

part of each CSF is a product of orbitals ,,, etc.and at

most one rij raised to some power See the review article by

Rychlewski20 and Sec. 2.4 of Rychlewski1 for a discussion of

Hy-CI wave functions. However, the calculations of Clementi et al. for H2 did not achieve the accuracy of purely Hy

calculations, the best of which were the KW calculations of

Wolniewicz.7

In this work we extend the work of Kolos and

co-worker5,6 and Wolniewicz7 to calculate the energies of

1

g

+ states of H2 using expansions in confocal elliptical coordinates with explicit inclusion of interelectronic distance

coordinates up through r12

7 . We calculate Born-Oppenheimer

BOenergies for various internuclear distances in the range

of 0.4–6.0 bohr. We also determine the BO ground state

energy more precisely than the best previous calculation.1,21

The calculations reported are similar to those described

in the classic paper of Kolos and Roothaan5 see this paper

for algorithmic details of our calculationsbut go far beyond

those as a reflection of the improved capability of modern

computers. We also note that the best previous calculation to

date on the H2 ground state is neither Hy KWnor Hy-CI

but the one employing exponentially correlated Gaussian

ECGfunctions.1,20,22,23 While the results that have been

obtained with this technique have been impressive, this technique suffers from the same inability to represent the electron cusp behavior at rij=0 as a strictly orbital CI calcula- a

Electronic mail: jim.sims@nist.gov

THE JOURNAL OF CHEMICAL PHYSICS 124, 094101 2006

0021-9606/2006/1249/094101/7/$23.00 © 2006 American Institute of Physics 124, 094101-1

Downloaded 02 Mar 2006 to 129.6.80.55. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
tion. Many authors have emphasized that the wave function

should have a cusplike behavior24 at r12=0 such that

 1





r12



r12=0

= 1

2

. 1

The ECG wave function cannot properly represent this behavior. However, progress towards overcoming this defect

has recently been made by Pachucki and Komasa.25,26 They

add terms linear in ri and rij to an ECG calculation to produce what they call a linear ECG LECGbasis.

II. WAVE FUNCTIONS

Since the time-independent, nonrelativistic electronic

Schrödinger equation in the Born-Oppenheimer BOor the

so-called clamped nucleiapproximation is not separable in

the electron coordinates, basis sets which incorporate the r12

interelectronic coordinate are most efficient. The wave function we use for the H2 ground state is

r1,r2= 

K=1

N

CKKr1,r2. 2

The terms K have the form

Kr1,r2= 1 + P12r12

 	1

n


2

j

	1

 ̄

n


2

j

 ̄

e−	1+ ̄

	2

, 3

where 	 and 
 are confocal elliptic coordinates27 and n, j,  ̄

n,

j

 ̄, ,  ̄

, and  must be specified for each expansion term K.

For 1

g

+ states, j+ j

 ̄ must be even. Terms in the KW wave

function differ from ours in that they include exponential

factors e−
1+

 ̄


2in 
1 and 
2 providing proper asymptotic

behavior of the wave function at large internuclear distance

R recently the KW approach has been extended by Kubacki

and Komasa28 to allow multiple parameters, e.g., 132 in the

calculation reported in Table III. Setting = ̄

=0 in the KW

wave function limits the range of R for which high precision

energies can be obtained. For us that range turns out to be

0.4–6.0 bohr we could handle R
0.4 with the current basis

if we used normalized orbitals, but even then, for technical

reasons our existing code would not allow R to approach 0.

In practice for the ground state one can let = ̄

 in each term

K as James and Coolidge originally did. Thus our wave function has = ̄

 for each K, but to speed up convergence we

use two sets of terms K and K

, one set with = ̄

=1 and

the other with = ̄

=2 we refer to these as two alpha expansions.

A. Results and discussion

In Table I we give the results of a 7034 term calculation

for R=1.4 bohr. In the table we show in each row the results

of adding to the previous calculation terms with an rij power

given by  and containing all n and j powers of 	1 and 
1,

respectively, that satisfy n+ jk and all  ̄n and j

 ̄ powers of 	2

and 
2, respectively, satisfying  ̄n+ j

 ̄k. In addition for a 1

g

+

state j+ j

 ̄ must be even, so terms with odd j+ j ̄ are not included. The last term selection criterion is that n,  ̄

n, j, and j

 ̄

must be . This and the value of = ̄for each line completely specifies the wave function expansion in sufficient detail that the calculation can be repeated, in contrast to any of

the previous large scale calculations on H2. This is one reason why we made no attempt to further select the terms

being used it then becomes difficult to specify the final

wave function terms in print without giving a detailed list.

The first eight lines are for adding successively high

powers of rij with 1=1.320 75 and then terms corresponding to 2=6.320 75 are added. Energies for each expansion

length N are tabulated as well as the improvement over the

previous level.

We note that E for =7 is greater then that for =6,

whereas one might expect a monotone decrease as the power

of r12 is raised. However, the terms in our expansion are of

two types: those with even powers of r12, which are basically

complicatedCI-type terms, and those with odd powers of

r12, which treat electron-electron interactions in an essentially different fashion. The E contributions from both of

these term types are indeed monotone decreasing as expected

but there is nothing that requires the combination to be similarly monotone decreasing, especially since one has reached

the point where it is getting difficult to describe the electron

TABLE I. Calculations of the BO energy of the ground state of the H2 molecule at internuclear distance R

=1.4 bohr. See text for what is tabulated in each column.

 k   N N Energy EhartreeE

0 11 10 1.320 75 0 1491 −1.161 518 240 453 517 4

1 10 10 1.320 75 1131 2622 −1.174 437 033 907 773 7 −0.129 1810−1

2 9 10 1.320 75 790 3412 −1.174 475 267 311 513 6 −0.382 3310−4

3 8 10 1.320 75 535 3947 −1.174 475 706 160 320 7 −0.438 8410−6

4 7 10 1.320 75 346 4293 −1.174 475 713 814 860 6 −0.765 4510−8

5 7 10 1.320 75 346 4639 −1.174 475 714 083 184 5 −0.268 3210−9

6 6 10 1.320 75 214 4853 −1.174 475 714 100 185 3 −0.170 0010−10

7 6 10 1.320 75 214 5067 −1.174 475 714 145011 5 −0.448 2610−10

0 8 8 6.320 75 535 5602 −1.174 475 714 190 684 3 −0.456 7210−10

1 8 8 6.320 75 535 6137 −1.174 475 714 210 523 3 −0.198 3910−10

2 7 8 6.320 75 346 6483 −1.174 475 714 216 401 7 −0.587 8310−11

3 6 8 6.320 75 214 6697 −1.174 475 714 218 711 6 −0.230 9910−11

4 6 8 6.320 75 214 6911 −1.174 475 714 219 818 1 −0.110 6510−11

5 5 8 6.320 75 123 7034 −1.174 475 714 220 075 5 −0.257 410−12

094101-2 J. S. Sims and S. A. Hagstrom J. Chem. Phys. 124, 094101 2006

Downloaded 02 Mar 2006 to 129.6.80.55. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
distribution with a single nonlinear parameter. In this connection, we note that the reordered two alpha expansion results presented in Table II exhibit no such anomalous behavior.

Based on the rate of convergence that is observed

here, we conclude that the energy has converged to 12

decimal places and the 13th digit is 0 or 1, i.e.,

E=−1.174 475 714 2201hartree. Optimization of  at the

7034 term level changes things only in the 15th decimal

16th digit, i.e., the energy surface in  space is very flat.

Varying the expansion terms beyond 7034 systematically

gave energy improvements occurring only in the 13th and

14th decimal places, further evidence that our result has converged to 13 digits.

In Table II we again give results for the 7034 term calculation but with terms reordered to show the energy contribution of each power of rij . Energies for each expansion

length N are tabulated as well as the improvement over the

previous level.

Our final energy is also listed in Table III and compared

with previous results. The 2400 ECG result in Table III has

been communicated to us by W. Cencek to correct his calculation reported in Ref. 21. As communicated to us by

Cencek,29 he gets −1.174 475 714 216 hartree and then the

numerical stability of the wave function deteriorates so that

the last digit starts fluctuating. Which means that the

−1.174 475 714 223 hartree reported in Ref. 21 is not a true

last digit upper bound since the last digits suffer from

numerical instabilities. Cencek’s error estimate is

110−10 hartree. As stated previously, we believe that our

result has converged to 1.174 475 714 2201, two orders of

magnitude better than the error estimate of the best previous

calculations.21,29 In Table III we label our calculations with a

single  as one alpha and those with both an 1 and an 2 as

two alphas. The JC, KW, ECG, and Hy-CI designations in

this and following tables are those of Rychlewski.1

In Kolos and Roothaan’s original paper5 the highest

power of r12 used was 2. Terms with r12

2 introduce the

equivalent of 12CI-type terms. So only - type

terms were explicitly correlated in Kolos and Roothaan’s

original study. r12

3 factors will introduce the equivalent of

12r12 terms. In Table III it can be seen that using only

a single , we are able to compute the energy to eight digits

seven decimal placeswith r123. It is interesting that with

r125 r12

5 puts in 12r12 we are able to do better than

all previous Hy calculations employing many terms with

high powers of 	1 and 	2 as well as the accurate 1200 term

ECG calculation of Cencek and Kutzelnigg.30 The energy

calculated with the ECG wave function was obtained only

after a very time-consuming optimization process in which

five nonlinear parameters per basis function which means

6000 parameters for the 1200 term wave functionwere optimized. In contrast the optimization for our single  wave

function is based on the selection of the terms of the form of

Eq. 3and the optimization of only the single parameter .

By going up to r12

7 we effectively correlate up to

TABLE II. Calculations of the BO energy of the ground state of the H2 molecule at R=1.4 bohr expressed in

hartrees. Terms are organized by rij power ; 1=1.320 75 and 2=6.320 75.

 N N Energy EhartreeE

0 2026 −1.161 523 300 229 836 0

1 1656 3682 −1.174 437 336 488 608 6 −0.129 1410−1

2 1146 4828 −1.174 475 280 326 470 9 −0.379 4310−4

3 749 5577 −1.174 475 706 942 023 1 −0.426 6110−6

4 560 6137 −1.174 475 714 069 212 3 −0.712 7110−8

5 469 6606 −1.174 475 714 216 389 9 −0.147 1710−9

6 214 6820 −1.174 475 714 219 924 8 −0.353 4910−11

7 214 7034 −1.174 475 714 220 075 5 −0.150 710−12

TABLE III. Comparison with previous explicitly correlated calculations for the BO energy of the ground state

of the hydrogen molecule R=1.4 bohrexpressed in hartrees.

Technique Author N Energy Ehartree

Many alphas Kubacki and Komasa Ref. 2822a −1.174 474 77

JC Bishop and Cheung Ref. 45249 −1.174 475 65

KW Kolos includes r12

3 Ref. 6370 −1.174 475 685

One alpha This work includes r12

3 3947 −1.174 475 706 160

Two alphas This work includes r12

3 5577 −1.174 475 706 942

KW Wolniewicz includes r12

6 Ref. 7883 −1.174 475 713 565

ECG Cencek and Kutzelnigg Ref. 301200 −1.174 475 714 037

One alpha This work includes r12

5 4639 −1.174 475 714 083

One alpha This work includes r12

7 5067 −1.174 475 714 145

Two alphas This work includes r12

5 6606 −1.174 475 174 216

ECG Cencek Ref. 292400 −1.174 475 174 216

Two alphas This work includes r12

7 7034 −1.174 475 174 220

a

132 nonlinear parameters.

094101-3 Ground state of hydrogen molecule J. Chem. Phys. 124, 094101 2006

Downloaded 02 Mar 2006 to 129.6.80.55. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
12-type products. Table III shows that a single nonlinear parameter is adequate for a 0.1 nhartree level of accuracy, but is not as good as the 2400 term ECG wave function.

By adding a second  at the r12

5 level we were able to exceed

our single  r12

7 result as well as the 2400 ECG result. Going

up to r12

7 with 1 and r125 with 2 achieves an 0.001 nhartree

level of accuracy.

Table IV summarizes previous results at the equilibrium

bond distance R=1.4011 bohr. The improvement of this

wave function over the ECG 2400 term wave function of

Cencek and Rychlewski parallels the relationship of ECG to

Hy at R=1.4 bohr.

We find in this work that very good energies are obtained

with our wave functions up to and including R=6.0 bohrs.

The fact that we do so well for such large R values without a

 nonlinear parameter in the wave function is presumably

because we have numerous terms with high powers of 
1 and


2. In Table V we show that the comparison of single 

results versus ECG results holds not only for a single internuclear distance but for a whole range of R using the 5067

term single  wave function. We are basically assuming that

the expansion terms used for R=1.4 will be adequate for

other R values. We are also assuming that since 1 and 2 are

not tightly coupled, we can use the 1 values from Table VI

for a good single nonlinear parameter calculation. It is striking how well these results agree with the 1200 term ECG

results, which are better than the best previous Hy results, the

KW results of Wolniewicz.7

Drake31 has pointed out the practical need for “doubling” basis sets so there is a natural partition of the basis set

into two distinct distance scales: one appropriate to the complex correlated motion near the nucleus and the other appropriate for further out. Drake uses just two sets of orbitals to

accelerate convergence. In our He calculations we found that

the importance of the second orbital exponent came in

around the 13th digit. In our present calculations we found it

difficult to get the tenth digit 0.1 nhartree accuracywithout

taking into account the detailed description of the wave function near the protons by incorporating elliptical orbitals with

large orbital exponents to describe the charge distribution

near the nucleus. This effect did not enter in the He calculations until we wanted more than 13 decimal places. The

greater need for two scale parameters is probably to be expected in H2 where much of the electron correlation energy

is associated with the region far from the nuclei, unlike He

where the space close to the nucleus is energetically important for both electron-electron and electron-nucleus interactions.

In Table VI we give our final results 7034 term wave

functionsfor the BO energy of the ground state of the

H2 molecule at internuclear distances ranging from 0.4

to 6.0 bohrs using two nonlinear parameters, 1 and 2. In

Table VI we also tabulate the contribution of all the 2 terms

versus R the column labeled E21012. As previously

discussed, we estimate the error in the energy at the internuclear distance R=1.4 bohr to be no more than 0.001 nhartree.

We have not investigated this question for other R values. In

the table we give energies up to 16 digits for the record and

also for the benefit of anyone wanting to repeat a calculation.

The table is very interesting, showing clearly how as R increases the electrons interact less and less, particularly

around the nuclei. Each nucleus has an electron and the probability for both being around the same nucleus is small, as

one would expect. Being able to demonstrate this is a side

effect of our way of choosing the wave function. In a calculation with lots of nonlinear parameters one would not be

able to see this so clearly.

In this paper we have restricted our calculations to the

BO energies. Recently there has been much interest in a

direct nonadiabatic which is non-BOvariational approach

using correlated Gaussian basis sets. The best result using

this approach is that of Cafiero et al.32 This approach yields

energies of about ten digits accuracy. However, the feasibility of this approach to systems with more than two electrons

has not been demonstrated. And for the two-electron case,

we have been able to obtain results using traditional methods

that yield energies to about 13 digit accuracy.

III. Hy VS Hy-CI

As there is a close relationship between the Hy wave

functions we use and an essentially equivalent Hy-CI treatment, these Hy results can provide insight into what is

TABLE IV. BO energy for the ground state of the hydrogen molecule at R=1.4011 bohr.

Technique Author N Energy Ehartree

Hy-CI Frye et al. Refs. 17 and 18586 −1.174 474 67

KW Wolniewicz Ref. 7883 −1.174 475 930 742

ECG Rychlewski et al. Ref. 46700 −1.174 475 931 197

ECG Cencek and Rychlewski Ref. 212400 −1.174 475 931 389

Two alphas This work includes r12

7 7034 −1.174 475 931 399 84

TABLE V. BO energy in hartreesof the ground state of the hydrogen

molecule at selected internuclear distances expressed in bohrcalculated

using the 1200 ECG wave function Table 2.8 in Ref. 1and our 5067 term

single  wave function terms defined in Table I, lines 1–8. The 1 values

used here are those presented in Table VI.

R 1200 term ECG 5067 term one 

1.0 −1.124 539 719 5 −1.124 539 719 462 06

1.2 −1.164 935 243 3 −1.164 935 243 354 69

1.4 −1.174 475 714 0 −1.174 475 714 145 01

1.6 −1.168 583 373 3 −1.168 583 373 287 03

1.8 −1.155 068 737 5 −1.155 068 737 548 57

2.0 −1.138 132 957 0 −1.138 132 957 074 81

3.0 −1.057 326 268 7 −1.057 326 268 835 43

4.0 −1.016 390 252 8 −1.016 390 252 933 31

6.0 −1.000 835 707 6 −1.000 835 707 653 04

094101-4 J. S. Sims and S. A. Hagstrom J. Chem. Phys. 124, 094101 2006

Downloaded 02 Mar 2006 to 129.6.80.55. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
needed in Hy-CI calculations. This can be seen by considering what a 12r12

2 term expands out to using

r12

2 = r12 + r22 − 2r1r2 cos 12. 4

For high precision, Hy wave functions2,5–7 must have the r12

power k1, whereas Hy-CI wave functions need only have

r12 raised to at most the first power.

In elliptical coordinates,

r12

2 = 

R

2 

2

p − 2q cos1 − 2, 5

where

p = 	1

2 + 
12 + 	22 + 
22 −2−2	1
1	2
2,

6

q = 	1

2 − 11 − 
12

	2

2 − 11 − 
22

1/2.

Expanding Eq. 5, the first term is just a polynomial in

	1, 
1, 	2, and 
2. Hence multiplying it by 12just

produces another 12charge distribution with higher

powers of ni and ji

. The second term in Eq. 5is proportional to q cos1−2. But q cos1−2is a linear combination of 12products with M =0 but without the exponential factors. This gets multiplied by the 12part

which does have the exponential factors. The 12

term gets multiplied by the 12part of the 12r12

2

term, so the final result is a special linear combination of 

and  CI configurations. As long as the powers of 	 and 


in the CI part of Hy-CI are greater than or equal to powers of

	 and 
 in the r12

2 term in the Hy expansion for both  and

, one will get an Hy-CI answer just as good as an answer

with Hy. Similar results apply for r12

4 and r126 this same point

was made by Clary12.

Although other radialfactors will be different in the

two methods, they are essentially equivalent methods for

two-electron problems. We note that an equivalent Hy-CI

calculation would have orbitals up to and including  with

r12 raised to at most the first power. A really good Hy-CI for

H2 would start out with a CI calculation including all the

different term types ,,,, etc.and then add to

this each of the terms multiplied by r12 but in general with

different nonlinear parameters. In particular, we note that

where calculations have been done CI for H2 and Hy-CI for

Heorbital exponents tend to increase at least for  pp for

Heterms resulting in significantly improved convergence of

the expansion.

In this paper we go beyond both the molecular Hy-CI

calculations of Clementi et al. and the purely Hy calculations

to achieve the best results for both types of calculations.

Since the Hy and Hy-CI methods are computationally

equivalent for all practical purposes for two-electron molecular systems we have demonstrated this for the helium

atom33, a Hy-CI calculation should be able to achieve comparable or better energies than these results with shorter expansions.

A. Comparison with experiment

For a review of the status of our knowledge of the dissociation energy of H2 until the year 2000 theoretical and

experimental, see Stoicheff.34 Currently the best experimental value for the dissociation energy D0 for H2 is

36 118.06210cm−1.

35 The best theoretical value is

36 118.058 cm−1.

1 The previous best theoretical values were

36 118.069 cm−1 Ref. 36and 36 118.049 cm−1,

7

respectively. These theoretical values are obtained starting from

accurate values for the Born-Oppenheimer BOelectronic

energies for the ground state of the H2 molecule over a range

TABLE VI. BO energy of the ground state of the H2 molecule at internuclear distances R from 0.4 to 6.0 bohrs. N=7034. See text for Table I for

what is tabulated in each column.

R 1 2 Energy EhartreeE21012

0.4 0.507 5 0.982 86 −0.120 230 341 177 864 4 285.9

0.5 0.605 1.281 44 −0.526 638 758 742 317 2 152.1

0.6 0.687 5 1.65 −0.769 635 429 485 356 8 118.1

0.7 0.78 2.0 −0.922 027 461 527 015 2 93.8

0.8 0.861 25 2.445 −1.020 056 666 360 138 9 88.5

0.9 0.939 5 3.007 5 −1.083 643 239 958 508 7 86.6

1.0 1.017 5 3.679 4 −1.124 539 719 546 579 1 84.5

1.1 1.093 25 4.445 −1.150 057 367 738 288 5 83.0

1.2 1.16 5.197 5 −1.164 935 243 440 028 1 85.3

1.3 1.227 5 5.945 −1.172 347 149 037 780 0 86.1

1.4 1.320 75 6.320 75 −1.174 475 714 220 075 5 75.0

1.5 1.39 6.945 −1.172 855 079 578 144 7 74.3

1.6 1.43 7.351 25 −1.168 583 373 370 926 3 83.9

1.7 1.525 7.695 −1.162 458 726 897 808 8 72.6

1.8 1.625 7.82 −1.155 068 737 610 807 1 62.2

1.9 1.705 8.07 −1.146 850 697 028 721 0 58.2

2.0 1.775 8.32 −1.138 132 957 131 503 5 56.7

2.1 1.851 8.445 −1.129 163 836 099 972 1 53.8

2.2 1.94 8.445 75 −1.120 132 116 847 639 1 48.9

2.3 1.985 9.07 −1.111 181 765 202 644 8 51.2

2.4 2.05 9.445 −1.102 422 606 009 297 8 49.9

2.5 2.125 9.57 −1.093 938 129 953 599 8 47.0

2.6 2.225 9.32 −1.085 791 237 393 588 7 40.4

2.7 2.287 5 9.695 −1.078 028 484 181 047 9 40.2

2.8 2.35 9.82 −1.070 683 233 478 409 5 38.9

2.9 2.425 10.07 −1.063 778 008 802 791 6 36.3

3.0 2.5 10.195 −1.057 326 268 869 243 9 33.8

3.1 2.6 9.945 −1.051 333 772 264 451 6 29.9

3.2 2.65 10.32 −1.045 799 661 428 733 8 29.0

3.3 2.75 10.07 −1.040 717 365 347 598 5 25.7

3.4 2.825 10.21 −1.036 075 395 186 919 5 23.7

3.5 2.875 10.695 −1.031 858 084 851 223 0 22.5

3.6 2.95 10.757 5 −1.028 046 308 375 876 6 20.5

3.7 3.025 11.07 −1.024 618 188 407 147 2 18.6

3.8 3.115 10.82 −1.021 549 795 529 910 9 16.6

3.9 3.167 5 11.32 −1.018 815 827 692 849 8 15.3

4.0 3.23 11.601 25 −1.016 390 252 947 128 3 13.8

4.2 3.38 11.882 5 −1.012 359 959 679 918 9 11.0

4.4 3.525 12.15 −1.009 256 516 258 663 2 8.7

4.6 3.6 13.32 −1.006 895 223 820 121 1 6.9

4.8 3.75 13.32 −1.005 116 006 098 095 2 5.3

5.0 3.9 13.195 −1.003 785 658 581 988 9 4.0

5.2 3.95 14.75 −1.002 796 816 309 543 1 3.2

5.4 4.1 14.00 −1.002 065 057 208 235 3 2.4

5.6 4.15 14.82 −1.001 525 251 885 354 9 1.9

5.8 4.2 15.57 −1.001 127 880 851 321 4 1.5

6.0 4.3 15.00 −1.000 835 707 654 227 9 1.2

094101-5 Ground state of hydrogen molecule J. Chem. Phys. 124, 094101 2006

Downloaded 02 Mar 2006 to 129.6.80.55. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
of internuclear distances R. This defines a BO potential energy PEcurve. The results presented in Table VI are the

definitive values for these BO energies defining the PE

curve. These can be used with separate calculations of adiabatic, relativistic, and radiative QEDcorrections to obtain a

theoretical dissociation energy see, for example, Ref. 1.

Our results are good enough to pin down the theory to 0.001

wave number if one can calculate the other corrections accurately enough. As new values of adiabatic, relativistic, and

QED corrections are obtained, they can be combined with

our BO energies to obtain new theoretical values. So maybe

the experiment needs to be redone using 0.001 wave number

as the level of accuracy to shoot for, that is, an order of

magnitude reduction in the current experimental uncertainties.

ACKNOWLEDGMENTS

The authors are deeply indebted to Wojciech Cencek for

supplying the details of his H2 calculation including doing a

reoptimization to help us pin down the true nonrelativistic

energy at R=1.4 bohr; Peter Mohr, John Rumble, and Maria

Belen Ruiz for helpful comments on this work; and to Jacek

Komasa and Maria Belen Ruiz for getting the authors interested in an Hy calculation for H2. The authors would also

like to thank Aaron Fein, Denis Lehane, Carl Spangler, and

Michael Strawbridge NISTfor parallel systems support.

One of the authors J.S.S.would like to thank Judith Devaney Terrill for supporting this endeavor. This material is

based in part on work supported at the Indiana University by

the National Science Foundation under Grant No. EIA0116050.

APPENDIX: SPECIAL FEATURES

OF THE CALCULATION

For a discussion of integral techniques needed to solve

the H2 electronic secular equation using ellipsoidal coordinates, see Kolos and Roothaan5,37 and references therein and,

in particular, the definitive paper of Rüdenberg.38

A parallel eigensolver proved essential for obtaining results over a range of R from 0.4–6.0 bohrs. This involves

solving the N-dimensional generalized eigenvalue problem

HC = SC A1

by the inverse iteration method.39 We solve this secular equation using our own portable parallel inverse iteration

solver.40 The generation of the matrices H and S is of order

N2 while the solution of the secular equation Eq. A1is of

order N3

. Since for us N is large and the interprocessor communication is only of order N2

, the parallelization of the H

and S matrices is largely a matter of convenience since no

restructuring of the code is required, just a suitable selection

of columns by each process. We calculated the basic integrals on each process, since this takes up only a small fraction of the total time, but parallelized the H and S matrix

construction since that allows the total memory needed to be

spread across the processors and eliminates a need to communicate matrix elements between processors. The result

was that the whole process is almost “embarrassingly parallel” with near linear “speedup.” For example, for a 4190 term

wave function we achieved a factor of 30 speedup on 32

processors for the order N3 step running on the National

Institute of Standards and Technology NIST147 processor

cluster of Pentium the identification of any commercial

product or trade name does not imply endorsement or recommendation by either the National Institute of Standards

and Technology or Indiana University, Athlon, and Intel

processors running REDHAT LINUX.

Almost all results reported in this paper were obtained

using quadruple precision 30+ digitsfloating point subroutines written in FORTRAN 90 using Miller’s41 quadruple

precision package which we augmented to interface to other

data types where native quadruple precision exists and complete agreement with Miller’s package was always obtained.

Much of the present code was developed on a conventional

desktop personal computer PCrunning WINDOWS 98 and

using FORTRAN compilers that did not include quadruple precision as a native FORTRAN data type. It was here that

Miller’s package proved extremely useful. Recently Hida’s42

very efficient quadruple precision package coded in C

has become available. We have our own portable FORTRAN 90

version of this excellent package,40 but there was no need in

this work to switch from the Miller package which we have

been using from the very beginning.

In addition to Miller’s quadruple precision package, we

made use of Brent’s arbitrary multiple precision MPfloating point arithmetic package43,44 for the Rüdenberg 

function38 on which all integrals depend. The Rüdenberg 

function is computed in multiple precision using Brent’s MP

package interfaced to FORTRAN 90 using the user defined

data-type mechanism with associated interface block definitions for the MP arithmetic operations and intrinsic functions.

The Rüdenberg  function is

nn

Ml ,= nn

Ml ,

= − 1

M l − M!

l + M!

	

1



d	Ql

M		

2 − 1M/2


e−	 	

n

	

1

	

dxPl

Mxx2 − 1

M/2e− ̄x

x

 ̄

n

+ e− ̄

	 	

 ̄

n

	

1

	

dxPl

Mxx2 − 1

M/2e−x

xn

 ,

A2

where the Pl

M are Legendre polynomials and the QlM are

associated Legendre functions of the second kind see Ref.

38, footnote 4 for detailed definitions.

Rüdenberg gives complicated recursion schemes for

raising n, n, l, and M starting from 00

00. Upper limits on n

and n depend on the values of l and M required, which in

turn depend on the expansion basis being used maximum

powers of 	 and 
. In the present work values of l up to 40

and n and n up to 81 were used. The recursion relations are

unstable for all indices, but for l they are particularly pathological involving the loss of one to two decimal digits of

094101-6 J. S. Sims and S. A. Hagstrom J. Chem. Phys. 124, 094101 2006

Downloaded 02 Mar 2006 to 129.6.80.55. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
accuracy each time l is raised by 1. Raising M from 0 to its

maximum value of 7 in this workinvolves only minor loss

of accuracy through differencing. Conversion back to quadruple precision was then done and the rij

 integrals were

evaluated in quadruple precision. To address these differencing problems we systematically increased the number of

decimal digits used for only the  part of the calculation up

to a maximum of 160 decimal digits.

We note that Kolos and Wolniewicz6,7 worked in double

precision and evaluated the rij integrals numerically using an

efficient charge distribution scheme also introduced by

Rüdenberg. This numerical approach is not applicable in

this work due to the general need for higher precision values

for all the basic integrals.

The code development used a desktop PC Athlon,

0.5 Gbyte memory, WINDOWS 98until the memory capacity

of the system was exceeded, at which time we moved over to

larger systems with more memory and somewhat later, parallelized the code to run on cluster systems at Indiana University analysis and visualization of instrument-driven data

AVIDD and NIST Raritan. The parallelization proved

extremely helpful not only for speeding the calculation but

also by spreading the total memory needed across the nodes

of the cluster.

1

J. Rychlewski, in Explicitly Correlated Wave Functions in Chemistry and

Physics Theory and Applications, edited by J. Rychlewski Kluwar Academic, Dordrecht, 2003, pp. 91–147. References unpublished work of

Cencek and Rychlewski.

2H. M. James and A. S. Coolidge, J. Chem. Phys. 1, 825 1933. 3

E. A. Hylleraas, Z. Phys. 54, 347 1929. 4A. D. McLean, A. Weiss, and M. Yoshimine, Rev. Mod. Phys. 32, 211

1960. 5W. Kolos and C. C. J. Roothaan, Rev. Mod. Phys. 32, 219 1960. 6

W. Kolos, J. Chem. Phys. 101, 1330 1994. 7L. Wolniewicz, J. Chem. Phys. 103, 1792 1995. 8 J. S. Sims and S. A. Hagstrom, Phys. Rev. A 4, 908 1971. 9

W. Woźnicki, in Theory of Electronic Shells in Atoms and Molecules,

edited by A. Jucys Mintis, Vilnius, 1971, p. 103. 10D. C. Clary, J. Phys. B 9, 3115 1976. 11D. C. Clary and N. C. Handy, Chem. Phys. Lett. 51, 483 1977. 12D. C. Clary, Mol. Phys. 34, 793 1977.

13D. C. Clary and N. C. Handy, Phys. Rev. A 14, 1607 1976. 14D. Frye, G. C. Lie, S. J. Chakravorty, A. Preiskorn, and E. Clementi, in

Modern Techniques in Computational Chemistry, edited by E. Clementi

Escom, Leiden, 1989. 15A. Preiskorn, D. Frye, and E. Clementi, J. Chem. Phys. 94, 7204 1991. 16E. Clementi and G. Corongiu, Rev. Inst. Fr. Pet. 51, 19 1996. 17D. Frye, G. C. Lie, and E. Clementi, J. Chem. Phys. 91, 2366 1989. 18D. Frye, G. C. Lie, and E. Clementi, J. Chem. Phys. 91, 2369 1989. 19E. Clementi, G. Corongiu, D. Estrin, E. Hollauer, and O. G. Stradella, Int.

J. Quantum Chem. 45, 511 1993. 20 J. Rychlewski, Adv. Quantum Chem. 31, 173 1998. 21 J. Rychlewski, in Handbook of Molecular Physics and Quantum Chemistry, edited by S. Wilson Wiley, New York, 2003, Vol. 2, pp. 199–218.

References unpublished work of Cencek and Rychlewski. 22W. Cencek and J. Rychlewski, J. Chem. Phys. 98, 1252 1993. 23W. Cencek, J. Komasa, and J. Rychlewski, J. Chem. Phys. 95, 2572

1995. 24T. Kato, Commun. Pure Appl. Math. 10, 151 1957. 25K. Pachucki and J. Komasa, Chem. Phys. Lett. 389, 209 2004. 26K. Pachucki and J. Komasa, Phys. Rev. A 70, 022513 2004. 27F. E. Harris, J. Chem. Phys. 32, 3 1962. 28T. Kubacki and J. Komasa, Ann. Polish Chem. Soc. 3, 457 2004. 29W. Cencek private communication. 30W. Cencek and W. Kutzelnigg, J. Chem. Phys. 105, 5878 1996. 31G. W. F. Drake, Phys. Scr., T 83, 83 1999. 32M. Cafiero, S. Bubin, and L. Adamowitz, Phys. Chem. Chem. Phys. 5,

1491 2003. 33 J. S. Sims and S. A. Hagstrom, Int. J. Quantum Chem. 90, 1600 2002. 34B. P. Stoicheff, Can. J. Phys. 79, 165 2001. 35Y. P. Zhang, C. H. Cheng, J. T. Kim, J. Stanojevic, and E. E. Eyler, Phys.

Rev. Lett. 92, 203003 2004. 36W. Kolos and J. Rychlewski, J. Chem. Phys. 98, 3960 1993. 37W. Kolos and C. C. J. Roothaan, Rev. Mod. Phys. 32, 205 1960. 38K. Rüdenberg, J. Chem. Phys. 19, 1459 1951. 39B. N. Parlett, The Symmetric Eigenvalue Problem SIAM, Philadelphia,

1998, Chap. 4. 40 J. S. Sims and S. A. Hagstrom, unpublished. 41A. J. Miller, http://users.bigpond.net.au/amiller, accessed on 10 September 2003. 42Y. Hida, http://www.cs.berkeley.edu/yozo/, accessed on 15 September

2005. 43R. Brent, ACM Trans. Math. Softw. 4, 57 1978. 44R. Brent, http://web.comlab.ox.ac.uk/oucl/work/richard.brent/pub/

pub043.html, accessed on 25 January 2005. 45D. M. Bishop and L. M. Cheung, Phys. Rev. A 18, 1846 1978. 46 J. Rychlewski, W. Cencek, and J. Komasa, Chem. Phys. Lett. 229, 657

1994.

094101-7 Ground state of hydrogen molecule J. Chem. Phys. 124, 094101 2006

Downloaded 02 Mar 2006 to 129.6.80.55. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp