Extrapolation to the Complete Basis Set Limit without Counterpoise

Dec 3, 2009 - realistic curve for the helium dimer by extrapolation to the complete basis set limit. ... for the BSSE by extrapolating the raw energie...
0 downloads 0 Views 506KB Size
Download PDF
J. Phys. Chem. A 2010, 114, 8505–8516

8505

Extrapolation to the Complete Basis Set Limit without Counterpoise. The Pair Potential of Helium Revisited† A. J. C. Varandas* Departamento de Quı´mica, UniVersidade de Coimbra, 3004-535 Coimbra, Portugal ReceiVed: September 12, 2009; ReVised Manuscript ReceiVed: October 26, 2009

Calculations using the full configuration interaction and coupled cluster methods are combined to yield a realistic curve for the helium dimer by extrapolation to the complete basis set limit. Extrapolations of the full configuration interaction energies alone are also done by reassignment of the basis hierarchic numbers. In no case has the common procedure of correcting the energies for the basis set superposition error by counterpoise been employed. The curves so generated agree well with the best available estimates, in one case within 1 cK (centikelvin) at the minimum. This corroborates our previous finding that such an approach provides an alternative to counterpoise that is more economical and general, besides yielding a more accurate curve. It may also conceivably justify that extrapolation to the complete basis set limit is, per se, all that is necessary to get a reliable potential function. 1. Introduction An accurate knowledge of the potential energy function (curve or hypersurface) that regulates atomic or molecular interactions is key in many fields of science and technology. Naturally, the accuracy requirements can vary with the system and property under study. For the helium dimer, such a requirement can be extreme by going as far down as a millikelvin, whereas for studying chemical reactions, it often suffices to attain 1 kcal mol-1 ) 503.219 K. With quantum chemistry at the heart of any rigorous calculation of such potentials, it comes therefore as no surprise that different methods are chosen according to the system and magnitude of the interaction involved. For example, molecular orbital (MO) methods are commonly utilized to calculate the potential functions of reactive systems, whereas perturbation theory can be more adequate in treating weak interactions. Although this was the conventional wisdom up to the 1990s, significant methodological advances and progress in computer hardware and software have conferred to MO methods a status of general recommendation. Because longrange forces can play a key role in reactive systems and even control the rates of chemical reactions, it seems therefore time to investigate the recurrent issue of whether MO methods can be used to treat with high accuracy, efficiency, and costeffectiveness such regions of the molecular potential. For this, the barely bound helium dimer can serve as an excellent prototype. In particular, we will examine to what extent the ubiquitous basis set superposition error (BSSE) that has been plaguing ab initio MO calculations since their inception can be eschewed without counterpoise1 (CP) since this is cumbersome and often unaffordable. Due to the weakness of van der Waals (vdW) interactions, the chosen electronic structure method should offer as much as possible a balanced description of both the complex and separated fragments. Stated concisely, the method should be size-consistent, with the Hartree-Fock (HF), full configuration interaction (FCI), and coupled cluster (CC) methods2,3 being examples that satisfy such a requirement. Conversely, truncation Part of the ″Klaus Ruedenberg Festschrift″. * E-mail: [email protected].

†

of the single-reference configuration interaction (CI) expansion such as to include only single and double (SD) electronic excitations is known to fail for size-consistency.2,3 The same holds for the popular MRCI (multireference CI) approach including single and double electronic excitations (if HF orbitals are used, the single-reference averaged coupled pair functional method, ACPF,4 which may be regarded as a modification of MRCI, is size-extensive for a system composed of equal subsystems). However, as noted above, an even more acute difficulty in the study of vdW interactions is the BSSE that manifests as a spurious attraction even in the absence of interaction. Stemming from the necessity of using finite sets of basis functions, BSSE has been studied a lot over the years3,5-11 (the list is by no means exhaustive), being commonly accounted for by the popular CP1 approach that consists of calculating the monomer energies with the full dimer basis. Interestingly, BSSE is general in the sense that it affects the weak interaction (in general, any) regions of any potential function while appearing unaffected by whether the method satisfies size-consistency or not.12 Regarding the noble gas dimers, past approaches have almost exclusively used atomic basis sets that are heavily augmented by diffuse orbitals,13,14 with additional functions in the bond region often being used also (e.g., see ref 15). The result is that such basis sets can generate substantial BSSE, which can then been dealt with by applying CP. Alternative approaches16-20 to reduce BSSE have also been developed. Although a full elimination of BSSE may be elusive, the question remains of how far one can go without explicitly accounting for it, an issue that will be analyzed in detail in the present work. Recently,12 we have suggested that one may largely correct for the BSSE by extrapolating the raw energies without CP (NCP) to the complete basis set (CBS) limit. Although such a work focused on the helium dimer, the same conclusion has been reassured21,22 in studies of interactions between a helium atom and fullerene molecules. A similar conclusion from thermochemistry calculations23 also came to our attention. Although we have further shown12 that CBS extrapolation of CP and NCP energies leads to the same extrapolated values, at one stage or another, results from CP have been invoked to

10.1021/jp908835v  2010 American Chemical Society Published on Web 12/03/2009

8506

J. Phys. Chem. A, Vol. 114, No. 33, 2010

work out the final result. More recently, Bytautas and Ruedenberg24 have determined an accurate potential energy curve for the neon dimer by CBS extrapolation of the energies obtained with nonaugmented basis sets and NCP. The analysis of their curve has furthermore yielded the correct dispersion coefficients. They did, however, report difficulties in using NCP energies when large basis sets were employed. This finding is somehow to be expected from the rather large BSSE that is met with such basis, especially for small cardinal numbers. Difficulties were also encountered in earlier studies on hydrogen-bonded complexes,6 where it was noted that CP and NCP HF and correlation energies converged unsystematically. Such a behavior has been attributed to the fact that the BSSE and the error from the incomplete description of the electronic Coulomb cusp were both present. This would suggest that an accurate description of weak interactions may be feasible only at the expense of a huge computational effort due to CP (for an heteronuclear pair, it amounts to six molecular calculations per geometry). However, to our knowledge, there is not a unique way of treating BSSE via CP in the case of many-body (fragment) systems due to many-body effects.25-28 For example, it is unclear whether two-body interactions in many-body systems should be computed in the basis set of the whole system or that of the dimer alone. Because (n g 3)-body effects are known to be considerably smaller than two-body ones, an even higher accuracy would be required for their calculation. Moreover, the Wigner-Witmer rules may allow more than one set of electronic states for the fragments, thus posing another difficulty29 to surmount especially for non-vdW interactions. Thus, although recognizing the positive aspects of CP in CBS extrapolation,12 it is convenient to set in perspective the type of errors that can be encountered when CP is left aside, as this is possibly the rule rather than the exception in the field of reaction dynamics. Although it has been found30 that CP may not warrant by itself a systematic improvement of the calculations, we shall not dwell on the view6,8 that CBS can only perform accurately together with CP. We will rather insist that CBS extrapolation should, by definition, render it unnecessary (at least in principle). Of course, one may assume that not all correlation-consistent basis sets will allow for accurate CBS results, but this is a weaker statement than the need to use both treatments. To investigate this issue, we will select here both the cost-effective correlation-consistent polarized valence (cc-pVXZ or VXZ) basis set of Dunning and his double diffusely augmented version (daug-cc-pVXZ or d-AVXZ) and show that rather accurate results can be obtained even when the calculated energies are not corrected with CP. Although calculations with two basis sets will be utilized, their cost can be smaller than that required by CP. Because He2 is only barely bound, the accurate computation of its potential energy curve provides an enormous challenge in electronic structure theory. Recent calculations have utilized the explicitly correlated multireference-averaged coupled pair functional (r12-MR-ACPF) method31 and MRCI32 and CC methods including corrections for both the FCI and the basis set limit.33 Additionally, refs 32 and 33 included corrections for the basis set limit taken from CCSD(T)-R12 calculations.34,35 Calculations using the quantum Monte Carlo36,37 (QMC) method and symmetry-adapted perturbation theory38 (SAPT) have also been performed. Moreover, Klopper39 reported an accurate CBS extrapolation of the He2 potential using CC methods that included single and double excitations (CCSD40) as well as perturbative corrections for connected triple excitations [CCSD(T)41]. Most of these works directly31,42,43 or indirectly (by

Varandas correcting for basis set unsaturation32,33) made use of explicitly correlated methods. Calculations15,44 have also been performed using supermolecular Gaussian geminal theory. In these, the bulk of the interaction energy was estimated using the Gaussian geminal implementation of CC theory with double excitations (CCD), while the smaller contribution due to single, triple, and quadruple excitations has been estimated using CCSD(T) and FCI methods with large atomic d-AVXZ basis sets and cardinal numbers up to X ) 6. To our knowledge, most such accurate MO calculations corrected for BSSE via CP. The accuracy that has been achieved is impressive, with the best estimates of the binding energy differing by small fractions of a kelvin and the most accurate interaction energy at a distance of R ) 5.6 a0 being estimated to be15 -11.0037 ( 0.0031 K. Because there is little room for numerical improvement, the focus of the present research will be on the estimation of how much BSSE remains after CBS extrapolation with NCP raw energies is done. As it will be shown, it resumes to a few percent of the well depth for the very weakly bound helium dimer, conceivably justifying the claimed accuracy of many MO calculations where such a technique has not been employed. The paper is organized as follows. In section 2, we describe the ab initio methodology that has been utilized to obtain the helium pair potential, while the results are reported and discussed in section 3. The conclusions are in section 4. For convenience, the interaction energies will be given either in microhartree or kelvin units (1 µEh ) 0.315775 K). 2. Methodology Let the theoretical methods be denoted as M ) Ma, Mb, ... and, by increasing size, let the basis sets be denoted as B ) B1, B2, and so forth. For example, B1 may be the VXZ basis set, B2 the d-AVXZ one, and so forth, while M stands for CI, CC, and so on. By adding the difference between the CBS/Ma/ B2 and CBS/Ma/B1 energies to the CBS/Mb/B1 one, a reliable estimate of the CBS/Mb/B2 energy will hopefully be obtained. Ma and Mb may in principle be any correlated methods, but if Mb is chosen to be the FCI method, then an accurate estimate of the FCI energy assumes the form FCI HF HF EMfFCI ∼ ECBS/B1 + ECBS/M - ECBS/M + b/B2 a/B1 dc dc B2 B1 - ECBS/M + δM - δM ECBS/M b/B2 a/B1 bfFCI afFCI

(1)

HF FCI where ECBS/B1 and ECBS/B1 are CBS extrapolations of converged estimates of the energies indicated in superscript with basis set B1 (correspondingly for B2). As usual, the superscript dc stands for the (dynamical) correlation energy; this and the HF energy B1 and add together to give the total energy. Similarly, δMfFCI B2 are corrections needed to go at each geometry from δMfFCI methods Ma and Mb to FCI with the indicated basis set. These too should represent CBS extrapolations (see later). If the FCI energy is also split into the HF and correlation components, eq 1 may then be written as

HF dc dc EMfFCI ∼ ECBS/M + ECBS/FCI/B1 + (ECBS/M b/B2 b/B2 dc B2 B1 ) + (δM - δM ) ECBS/M a/B1 bfFCI afFCI

(2)

with the difference on the right-hand side being denoted as B2/B1 B2 B1 B1 ) δMfFCI - δMfFCI ∆EMfFCI . If the correction δM were afFCI known exactly, this term added to the corresponding correlation dc energy, Edc CBS/Ma/B1, would cancel ECBS/FCI/B1. Removing the latter may, however, not be advantageous from the practical point of view (see discussion below on the BSSE). Moreover, we may

The Pair Potential of Helium Revisited

J. Phys. Chem. A, Vol. 114, No. 33, 2010 8507

Figure 1. Hartree-Fock and correlation (as obtained from the corresponding FCI calculations) energy differences calculated with the d-AVXZ and VX basis sets. The lines are mere spline-fits to guide the eye and hence illustrate better the rather different rates of convergence of the HF and correlation energies with X. After 4.5 Å, the various lines flatten out and are hardly distinguishable within the scale of the plot. B2 B1 later wish to ignore the difference ∆EB2/B1 MfFCI ) δMfFCI - δMfFCI. Otherwise, it would give

HF dc B2 EMfFCI ∼ ECBS/M + ECBS/M + δM b/B2 b/B2 bfFCI

(3)

an equation that has been successfully utilized in ref 12 both with CP and without it, with Mb being there the CCSD(T) method. Note from eqs 2 and 3 that, if the basis B1 and B2 are the same (the cardinal numbers are then also identical), one FCI to simply presupposes that B2 is sufficiently large for ECBS/B2 be approximately the exact FCI energy. If the two basis sets differ, then one further presupposes from eq 2 that the amount of recovered correlation energy is essentially independent of the basis. We will call the above approach, which follows the spirit of the popular focal point analysis,45 method I. In the present work, Mb in eq 2 will be the CCSD(T) method, while the basis sets are B1 ≡ VXZ and B2 ≡ d-AVXZ. Method II refers to direct extrapolation of FCI energies calculated with the largest affordable basis set. Thus, if eq 3 is utilized with Mb ≡ FCI, it results in the mere addition of the extrapolated HF and FCI correlations to basis set B2 since B2 ≡ 0. Such a brute force scheme is hardly affordable δFCIfFCI for X g 5, even for a four-electron system such as the title one. A further remark to note is that method I can partly account for the BSSE in the dynamical correlation. To see this, assume for simplicity that Ma ≡ Mb. Let now the wave function associated with basis set B2, Ψ2, be written as Ψ2 ) Ψ1 + φ where φ is a small correction, which should be a good representation for basis functions of the correlation-consistent family. The energy difference (E2 - E1) will then be ∆E ) 2〈Ψ1|H|φ〉 +〈φ|H|φ〉, where H ) HA + HB is the total Hamiltonian for the A and B subsystems. Thus, to first-order in φ, the effect of BSSE in ∆E will be proportional to the overlap of the basis set enhancement in Ψ1 (due to BSSE) with φ. This should be small if Ψ1 has a small BSSE. Indeed, this appears to be the case when B1 is the VXZ basis. It is also apparent from Figure 1, which shows that the HF and correlation energy differences display a more regular decreasing pattern in absolute value with X. Of course, CBS extrapolation of the HF energy with CP should be relatively inexpensive and might possibly enhance accuracy,6,8,12 but it will not be pursued in the present work.

2.1. Ab Initio Calculations. All energies have been calculated with the Molpro suite of electronic structure programs.46 Since the FCI energies can be extremely time-consuming, we give a sample of the calculated points in Table 1. All others at the HF and CCSD(T) levels have either been given12 or are relatively easy to calculate and hence need not be tabulated. Almost all calculated geometries are indicated by the symbols in Figure 2. The largest bond distance considered here is R ) 891 Å, although a bond distance of R ) 23 Å gives a total energy that differs from that beyond the last reported digit, thus a negligibly small difference for the purpose of the present analysis. 2.2. CBS Extrapolation Schemes. CBS extrapolation techniques have been discussed in detail elsewhere.3 Suffice it to say that, for total energies, the convergence of the correlation part is significantly slower than that of the HF part, and hence, one should treat them separately during such a procedure. Additionally, the formalism used for the CBS extrapolation of the correlation energy finds solid theoretical support47 (see also refs 3 and 48, and references therein), while the one utilized for the HF energy relies more on physically intuitive expectations and the accuracy of the CBS results that are obtained49-54 when judged against benchmark values. Thus, only a brief survey of the formalisms employed in the present work will be presented. The uniform singlet and triplet pair extrapolation (USTE48) scheme (see also refs 55 and 56) for the CBS extrapolation of the correlation energy is based on the three-parameter protocol

EXcor ) E∞cor +

A3 (X + R)

3

+

A5 (X + R)5

(4)

with A5 defined by

A5 ) Ao5 + cAm 3

(5)

where E∞cor, A5(0) ) A5o, A3, c, and R are parameters. By fixing R, A5o, c, and m from other criteria, eq 4 can be transformed into an effective two-parameter (E∞cor, A3) rule,48 with the parameters determined from a fit to the raw dynamical correlation energies. Using the USTE model, we have shown57 that both the full correlation in systems studied by the popular singlereference Møller-Plesset (MP2) and coupled cluster [CCD and CCSD] methods as well as its dynamical part in MRCI(Q) calculations48 or even correlation energies obtained by correlation energy extrapolation via intrinsic scaling58 can be accurately extrapolated to the CBS limit. From the dynamical correlation of 24 calibrating systems studied48 by the MRCI(Q) method, the optimum values of the “universal-like” parameters were found to be A5o ) 0.003769 and c ) -1.1784771 Eh-5/4, with m ) 5/4 and R ) -3/8. Corresponding parameters have been given for other methods with R ) -3/8, namely, A5o ) 0.1660699, c ) -1.4222512 Eh-1, and m ) 1 for the CC-type methods. Naturally, one may expect these coefficients to vary with the method used for the electronic structure calculations and the basis set. Hopefully, such a dependence is not significant for methods and basis sets that belong to related families, and we have explored the extendibility (“universality”) by showing that they lead to accurate results even for systems that have not been part of the calibrating set.57,59 For example, we have shown12 that the coefficients (A5, A3) determined from the CI(Q)/dAVXZ, CCSD/d-AVXZ, and CCSD(T)/d-AVXZ calculations for the title system are all well described by eq 5. Because there may be difficulties close to the origin (small values of A3 and A5), we have suggested a convenient variant of the method that

8508

J. Phys. Chem. A, Vol. 114, No. 33, 2010

Varandas

TABLE 1: Sample of Calculated Energies for the Helium Pair Interaction at the FCI Level of Theorya VXZ R/Reb

R/Å

X)D

0.50 1.485 -5.75489395 0.75 2.228 -5.77456448 0.90 2.673 -5.77515762 0.98 2.911 -5.77519251 1.00 2.970 -5.77519466 1.02 3.029 -5.77519567 1.04 3.089 -5.77519594 1.10 3.267 -5.77519481 1.20 3.326 -5.77519424 1.25 3.712 -5.77519158 1.50 4.455 -5.77519022 3.00 8.910 -5.77518967 4.00 11.880 -5.77518966 8.00 23.760 -5.77518966

d-AVXZ

T

Q

5

6

D

T

Q

5

-5.78079210 -5.79976614 -5.80040274 -5.80045846 -5.80046414 -5.80046802 -5.80047059 -5.80047341 -5.80047342 -5.80047029 -5.80046607 -5.80046436 -5.80046434 -5.80046434

-5.78560950 -5.80418315 -5.80477610 -5.80482386 -5.80482841 -5.80483141 -5.80483327 -5.80483475 -5.80483446 -5.80483016 -5.80482451 -5.80482180 -5.80482176 -5.80482176

-5.78731815 -5.80571604 -5.80627262 -5.80631359 -5.80631697 -5.80631898 -5.80631999 -5.80631970 -5.80631903 -5.80631351 -5.80630713 -5.80630382 -5.80630378 -5.80630377

-5.78798222 -5.80631287 -5.80684606 -5.80688111 -5.80688354 -5.80688474 -5.80688506 -5.80688318 -5.80688211 -5.80687505 -5.80686800 -5.80686429 -5.80686425 -5.80686424

-5.75892844 -5.77862239 -5.77921158 -5.77924615 -5.77924688 -5.77924604 -5.77924415 -5.77923564 -5.77923259 -5.77921676 -5.77920178 -5.77918935 -5.77918874 -5.77918872

-5.78215801 -5.80071610 -5.80123348 -5.80125992 -5.80126089 -5.80126085 -5.80126013 -5.80125583 -5.80125408 -5.80124300 -5.80122812 -5.80121661 -5.80121631 -5.80121626

-5.78619214 -5.80457463 -5.80508794 -5.80511392 -5.80511432 -5.80511360 -5.80511213 -5.80510566 -5.80510331 -5.80509064 -5.80508016 -5.80507334 -5.80507323 -5.80507321

-5.78757859 -5.80590749 -5.80641585 -5.80644137 -5.80644186 -5.80644129 -5.80644000 -5.80643418 -5.80643205 -5.80642023 -5.80640948 -5.80640399 -5.80640390 -5.80640389

a Distances in angstroms; energies in hartree. calculations have been performed.

b

Ratios (relative to the experimental equilibrium distance of Re ) 2.97 Å) for which the

reliability when CBS extrapolations from small cardinal numbers are performed, as in method II. The procedure begins by first recognizing that an extrapolation scheme is best discriminated as performing well or poorly if the extrapolation assumes a linear form. To achieve this, we rewrite eqs 4 and 5 as

EXcor ) E∞cor +

A3

(6)

X˜3

which is equivalent to define the following new “cardinal number” (denoted heretofore as the hierarchic number since it is generally a real positive number rather than an integer)

(

X˜ ) (X + R) 1 + Figure 2. Calculated raw interaction potentials for the helium dimer with VXZ (open symbols) and d-AVX (solid symbols) basis sets. The solid lines represent HFACE fits to the points and serve to guide the eye by connecting them. Legend as that in Figure 1. Shown free of symbols are the CBS-extrapolated curves (CBS/FCI/VXZ, dashed; method I, solid; method II, dash-dot). The inset shows the corresponding HF energies and curves.

overcomes such difficulties, so-called GUSTE (generalized USTE). Since an accuracy as high as one possibly can get is a clear demand for the title system, most extrapolations will be done here by choosing (E∞cor, A3) such as to reproduce the calculated raw energies for the highest affordable cardinal numbers. Yet, to maintain the scheme tractable for other common interactions, we limit such values to the pair (5, 6) for the smallest basis (VXZ). It turns out that USTE and GUSTE give essentially the same results under such conditions, and hence, we will employ only the USTE scheme in the calculations reported here. Since no set of parameters has been given for FCI calculations, one may wonder which set to use. Clearly, the difference is expected to be very small since both the CI and CC sets should perform rather accurately for the title system. This is particularly true when extrapolating from the highest affordable pair of cardinal numbers, the so-called cardinal number pair. We will therefore utilize it to CBS extrapolate the FCI energies in method I. However, FCI calculations become prohibitively expensive when using the d-AVXZ basis with X g 5. Thus, we suggest next a route that may help in the endeavor of enhancing

A5 /A3 (X + R)2

)

-1/3

(7)

where A5 and A3 are obtained as described above. Recalling now that the d-AVXZ basis differs from VXZ in that the former is augmented with a double diffusely set of basis functions, one expects any reassignment of hierarchic numbers that is optimal for one will also be adequate for the other. We have chosen here to make this reassignment from a fit to the FCI energies with X ) T, Q, 5, and 6 to eqs 6 and 7 since this gives reassigned hierarchic numbers (denoted as h) at the equilibrium geometry of 2.97 Å that differ typically by ,1% for all such cardinal numbers. As illustrated in Figure 3, the new hierarchic numbers are chosen such that they fall on the straight line so obtained. Such h numbers can be found by solving numerically the root of a simple nonlinear equation, yielding h ) 2.067731, 2.999982, 3.995798, 5.039294, and 5.941800), in clear correspondence to D, T, and so on. Thus, they show deviations of 3.4, -0.0, -0.1, 0.1, and -0.1% relative to the original cardinal numbers. As seen, except for X ) D, such deviations are quite small. However, it is visible from Figure 3 that extrapolations from at least the two lowest cardinal numbers can be problematic. Note also that the VDZ basis appears to be better calibrated (saturated) than the VTZ and VQZ when the V5Z and V6Z ones are selected for redefining the hierarchy, and therefore, there is not a good criterion for excluding the former from the set used for the CBS extrapolation procedure. We emphasize that h is found from the requirement that it falls on the straight line defined from a least-squares fit to the FCI/VXZ correlation energies for X ) T, Q, 5, and 6. Of course, slightly different assignments would occur if the reference were chosen differently. For a similar reason, the hierarchic numbers will depend

The Pair Potential of Helium Revisited

J. Phys. Chem. A, Vol. 114, No. 33, 2010 8509

δE∞cor

)

3 cor 3 Ecor 2 X2 - E1 X1

X31 - X32

+

cor 3 3 Ecor 2 (δ2 + X2) - E1 (δ1 + X1)

(δ2 + X2)3 - (δ1 + X1)3

Figure 3. Set up of a new hierarchy number from FCI calculations with the VXZ basis. Shown in the inset is a plot showing fits to various pairs of cardinal numbers (note that the curves fitted using cardinal numbers miss the T and/or Q points). See text.

on the set of coefficients used to parametrize USTE. For example, if the set of CC coefficients is used instead to rehierarchize the basis from the same FCI energies, one obtains at the above geometry h ) 2.169354, 3.002737, 3.963168, 5.043094, and 6.070349, with errors of 8.5, 0.1, -0.1, 0.1, and 1.2% relative to the cardinal numbers. For an easier notation, if no confusion arises, we will refer to them by the small capitals of the Greek numerical prefixes for the closest integers, h ) d, t, q, p, h (if an ambiguity arises such as hexa and hepta, we must resort to the numbers themselves, define primed letters, etc.). The observed trend is rather similar to what has been found above, with the DZ basis set corroborating common wisdom by not following the pattern of all others. Conversely, by definition, all h numbers follow a straight line pattern. The CBS extrapolations described later for the FCI/d-AVXZ correlation energies in method II will then consist of fits to an expression similar to eq 6 but reading

EXcor ) E˜∞cor +

A˜3 h˜3

(9)

Note that δE∞cor can be surprisingly large. For example, if (X1, X2) ) (T, Q) with d-AVXZ correlation energies of (-0.07817805, -0.08182562) at the equilibrium geometry of He2 and δ1 ) δ2 ) 0.01, the error on the extrapolated correlation energy is -3.6 K, just about one-third of the magnitude of the well depth itself. However, if the errors are of opposite sign, (δ1, δ2) ) (-0.01, 0.01), then δE∞cor ) -26 K, which is over twice the well depth. Indeed, even with (δ1, δ2) ) (0, 0.01), the error is in absolute value (-14.7K) more than the well depth. For a given error in X2 (say δ2 ) 0), the nearly linear variation with the error in X1 can be understood by fixing δ2 ) 0, yielding

δE∞cor

) ) -

cor 2 2 3 δ1(Ecor 1 - E2 )(δ1 + 3X1δ1 + 3X1)X2

(X31 - X32)[(δ1 + X1)3 - X32] 3(Ecor 1

2 3 - Ecor 2 )X1X2

(X31

-

X32)2

δ1 + O(δ21)

(10) where the second term on the right-hand side represents the first-order term in a Taylor series expansion. Clearly, for a pair (X1, X2) with associated energies (E1cor, E2cor), δE∞cor varies approximately linearly with δ1 for small values of the latter. Of course, such errors are in the absolute correlation energy, being hopefully smaller in the interaction potential when Ecor ∞ (R ) ∞) is subtracted. As could have been anticipated, the error decreases with increasing pair (X1, X2) as this produces a larger denominator in eq 9 and hence a smaller unsigned value of δEcor ∞ . A similar analysis could have been done with USTE, although it had to be performed numerically due to the dependence of A5 on A3. However, the general trend should be similar as the term in X-5 is a mere (yet relevant) correction. Thus, the results recommend the hierarchy of the basis to be set at a single geometry such as to enhance error cancellation. This will be the approach followed here by choosing the reassignment at R

(8)

where E˜∞cor and A˜3 are parameters to be determined in general from a least-squares fit to the calculated energies. The results for R ) 2.97 Å of USTE (d-p) are compared with those that originated from the traditional USTE (D-5) and (Q, 5) fits in Figure 4. The important observation is that all calculated energies still fall nicely over a straight line, which contrasts with the fits obtained by the traditional procedure (shown by the dashed line). It is now pertinent to ask whether the reassignment of the basis to new hierarchic numbers should be done pointwise (i.e., at each geometry). Although the answer would seem obvious and affirmative, this is not quite so due to the extreme weakness of the bonding. This can be demonstrated analytically. Consider the popular CBS rule3,6,39 based on the inverse cubic term alone, that is, eq 6 with X˜3 as is or replaced by the cubed cardinal number itself (X3). Let us consider the latter. Assume now that cor the extrapolation is from the pair (Ecor 1 , E2 ) for cardinal numbers (X1, X2) and that δ1 and δ2 are errors due to the hierarchy assignment of the basis. The error in the CBS-extrapolated energy will then be

Figure 4. Illustration of a CBS extrapolation of FCI/d-AVXZ energies using cardinal numbers (dashed lines) and the new hierarchical setting (solid lines). The open circles for X ) T-6 are essentially coincident with the solid dots. A similar (or even slightly better) performance of the h hierarchical number is observed for other geometries of the dimer.

8510

J. Phys. Chem. A, Vol. 114, No. 33, 2010

Varandas

) 2.97 Å. For further consistency, CBS extrapolations will be done by fitting (d-p) energies rather than the pair (p, q) as numerical noise is enhanced in this case [the denominator of eq 10 is smaller]. Of course, errors may still arise through the cardinal numbers themselves. Extrapolation schemes are also available for the HF and CAS uncorrelated energies. One is the three-parameter exponential rule,49-51 EXCAS ) E∞CAS + A exp(-bX), where E∞CAS, A, and b are determined from a fit to the raw energies, typically for the three or more highest affordable cardinal numbers. In ref 48, we have suggested a variant of this rule whereby the value ECAS ∞ is obtained from a fit to X ) T, Q, 5, and 6 energies followed by averaging with the raw HF energy for the highest fitted value of X (typically 6). Another CBS scheme is52,53,60

EXHF ) A + B(X + 1) exp(-9√X)

(11)

which has been amply tested by Karton and Martin and will be here referred to as KM. Since the HF energies predicted from our protocol show good agreement with the two-point KM one for the (5,6) case, we will adopt here their own. A final remark to note is that CBS extrapolation can be combined with correlation scaling61,62 to obtain significant computer time savings.59,62 To limit errors to a minimum, such a hybrid (correlation scaling/extrapolation) procedure will not be followed in the present work. 3. Results and Discussion We now turn to the analysis of the CBS-extrapolated potentials from the present work and a comparison with others available in the literature. Before going into the details, we should note that a comparison with the correlation energies extrapolated from the traditional X-3 law using the (5, 6) cardinal numbers has been made elsewhere,12 where it has been noted that both yield comparable results for extrapolations based on the (5, 6) correlation pair. This is not surprising since USTE has been primarily targeted for extrapolations from low hierarchic numbers. Since not much more needs to be said about the CBS extrapolation schemes, we now address the correction term B2/B1 B2 B1 ) δMfFCI - δMfFCI . In our previous paper on the title ∆EMfFCI dimer,12 we have followed Klopper39 by noting that MRCI calculations in the IO301 basis,32 FCI calculations in the AVQZ and d-AVTZ basis,33 MR-ACPF calculations in the d-AV6Z basis,63 and the Gaussian geminal implementation of CCSD theory including effects of triple and quadruple excitations44 (and employing the conventional orbital approach and very large augmented correlation-consistent basis sets extended by sets of bond functions) indicate that at R ) 5.6 a0, the step from CCSD(T) to FCI adds δE(T)fFCI ) -(323 ( 5) mK to the well depth. In fact, by using even larger basis sets and extrapolations, this correction has been slightly reduced to15 δE(T)fFCI ) -(318 ( 3) mK, a value that has been utilized in ref 12. For method I of the present work, we needed such corrections as a function of bond distance and hence calculated them with both the VXZ and d-AVXZ basis. As Figure 5 shows for R ) 2.7 Å, there is no conceivable simple relation with X that warrants a safe CBS extrapolation. This can be reinforced by examining plots for other bond distances. A similar observation comes from refs 44 and 64, where the authors have obtained estimates of δE(T)fFCI using various Dunning-type basis set sequences restricted to X e 5 (except the AVXZ sequence which used up to X e 6) and some interaction-energy-optimized basis sets. They have also found that extrapolations were not reliable due to a significant spread of the extrapolated values. For R ) 5.6

Figure 5. Variation of the energy correction δE(T)fFCI with X for the VXZ and d-AVXZ basis sets. The error bars in solid show the estimated uncertainties (the dashed one is the more ample one that embraces the values for the T-6 energies; see text.)

a0, the values computed44,64 with the largest basis sets used ranged from -321 to -328 mK, from where they have made VXZ ) (-323 ( 5) mK to the above recommendation of δE(T)fFCI cover both ranges. From their more recent analysis with X e 6 d-AVXZ seemed to increase results, the computed values of δE(T)fFCI monotonically with X (becoming less negative), a trend also corroborated by our calculations. In fact, for R ) 5.6 a0, their computed value in the d-AV6Z basis (326 orbitals) was -324mK. They were then led to conclude that the CBS limit should be still smaller in magnitude, which they estimated in the range between -318 and -321 mK. Because an unsigned linear growing (larger negative values) with X is physically unsound and the results for other geometries do not warrant a monotonous decreasing behavior with X, the following strategy . First, we have accepted has been used to estimate ∆EVXZ/d-AVXZ (T)fFCI as an estimate for the VXZ basis the average of the corrections obtained with the three largest cardinal numbers (Q-6). The VXZ 〉 so calculated are given in Table 2, average values 〈δE(T)fFCI VXZ ) (-178 ( 59) mK. with the value at 5.7 Å being δE(T)fFCI For the d-AVXZ basis, we base our corrections on the values obtained with X ) 5 but slightly correct it to take into account the new d-AV6Z result of ref 15. Then, we note that our d-AV5Z ) -332 mK, with prediction at 2.97 Å ) 5.613 a0 is δE(T)fFCI an unsigned value slightly larger (-351mK) at 2.94 Å. From these two results, we have estimated by linear interpolation a value of 338 mK at 5.6 a0, which is in fair agreement with the result reported by Cencek et al.44 for the d-AV5Z basis. By further recognizing that the values for this basis allow CBS extrapolation, a conjecture is made that the difference VXZ/d-AVXZ follows a similar pattern (this implies an expected ∆E(T)fFCI VXZ 〉 determined similar variation of the average value 〈δE(T)fFCI above, although carrying the associated uncertainties, for example, (59 mK at R ) 2.97 Å). The situation is unlikely to change by carrying out CP as the BSSE should be similar at the CCSD(T) and FCI levels, as indeed suggests the significant spread encountered in the calculated data.44,64 Accordingly, the VXZ/d-AVXZ has been downscaled by a factor above value of ∆E(T)fFCI d-AVV5Z value into of 0.941 ()318/338) such as to bring the δE(T)fFCI agreement with the corresponding CBS (5, 6)-extrapolated values for the same basis. In summary, the major source of error in our calculations arises from the corrections δE(T)fFCI, particularly with the VXZ basis set. A very conservative estimate of the error would then be (59 mK, which embraces the X )

The Pair Potential of Helium Revisited

J. Phys. Chem. A, Vol. 114, No. 33, 2010 8511

TABLE 2: Some CBS extrapolated energies by methods I and II for the helium pair interaction potentiala R/Reb

R

HFc

cor(I)

VXZ δE(T)fFCI

0.50 0.75 0.90 0.98 1.00 1.02 1.04 1.10 1.20 1.25 1.50 3.00 4.00 8.00

1.485 2.228 2.673 2.911 2.970 3.029 3.089 3.267 3.326 3.712 4.455 8.910 11.880 23.760

6584.882 262.116 35.076 11.801 8.982 6.820 5.187 2.262 1.721 0.290 0.003 -0.002 0.000 0.000

-651.446 -105.050 -37.953 -22.655 -19.996 -17.662 -15.637 -10.943 -9.756 -4.820 -1.497 -0.019 -0.004 0.000

-6.123 ( 0.303 -1.116 ( 0.236 -0.359 ( 0.102 -0.204 ( 0.064 -0.178 ( 0.059 -0.157 ( 0.052 -0.138 ( 0.048 -0.096 ( 0.034 -0.085 ( 0.032 -0.044 ( 0.019 -0.016 ( 0.013 -0.004 ( 0.005 -0.004 ( 0.008 -0.005 ( 0.007

d-AVXZ δE(T)fFCI

d

-6.619 -1.566 -0.609 -0.376 -0.332 -0.297 -0.262 -0.189 -0.171 -0.088 -0.035 -0.010 -0.006 -0.006

VXZ/d-AVXZ ∆E(T)fFCI

e

-0.496 -0.450 -0.250 -0.172 -0.154 -0.140 -0.124 -0.093 -0.086 -0.044 -0.019 -0.006 -0.002 -0.001

cor(II)

total(I)f

total(II)

-662.196 -106.160 -38.392 -23.014 -20.308 -17.949 -15.879 -11.088 -9.862 -4.739 -1.373 -0.013 -0.002 -0.000

5933.436 157.066 -2.877 -10.853 -11.014 -10.843 -10.449 -8.681 -8.035 -4.529 -1.495 -0.020 -0.004 0.000

5922.674 155.956 -3.316 -11.212 -11.326 -11.129 -10.691 -8.826 -8.141 -4.448 -1.371 -0.015 -0.002 -0.000

a Distances in angstroms; energies in kelvin. b See Table 1. c CBS/HF/d-AVXZ energy, applicable to both methods I and II. The negative value at 8.910 Å is a remnant of the strong BSSE with the d-AVXZ basis and possible numerical truncation errors. d Estimate of error available VXZ/d-AVXZ d-AVXZ VXZ ) 〈δE(T)fFCI 〉 - 〈δE(T)fFCI 〉 before scaling by the factor of 318/338. See the text. only for the equilibrium geometry. e Value of ∆E(T)fFCI f Differences of (1 mK relative to the sum of the tabulated parts may be obtained in method I due to truncation of the numbers at the third decimal figure. See also the text.

V5Z T-6 estimates of δE(T)fFCI at R ) 2.97 Å, although a more realistic choice would be half or so of such a value (see later). Near equilibrium, this amounts to an uncertainty of (30 mK, which turns out to be just slightly larger than the difference d-AV5Z at R ) 5.6 a0 and the lowest between our value of δE(T)fFCI CBS estimate (318 mK) reported in ref 15. Such a similarity of values, and the observation that the decreasing trend (larger negative values) is far from being observed at other geometries with the VXZ basis, suggests that the above average value for VXZ d-AVXZ once affected by the error in δE(T)fFCI is possibly a δE(T)fFCI fair measure of the former energy difference. Thus, the original uncertainty has been cut down by 30% to leave still a large margin of confidence. Of course, the large error bar still remaining is the price to pay for the aim at obtaining an accurate energy with a modest basis sets such as VXZ, particularly when dealing with one of the weakest bonds. Another source of error is due to having used the default arithmetic of Molpro, which may add up to (3 mK, perfectly negligible at the equilibrium region when compared with the previous one. Additionally, there may be a systematic error due to having not gone beyond the (5, 6) cardinal number pair. This should manifest predominantly in the correlation component of the energy and will be tentatively examined later. Regarding the variances quoted in the previous paragraph, they compound quadratically15 to an error of (30 mK or so near Re. As above, a reduction of 30% has been considered to estimate the error at other geometries. As Table 4 and Figure 6 show (the analytic curve in this figure will be described later), the CBS extrapolations of method I predict a value close to the best available estimate at the minimum. Of course, at this level of accuracy, one should care also about relativistic and even QED corrections, not to mention nonadiabatic effects, but this is outside of the scope of the present work. We now address method II, that is, the direct CBS extrapolation of FCI correlation energies calculated with the expensive d-AVXZ basis set (a point with the largest affordable 5Z basis takes about two weeks in a Quad core 2.4 MHz processor). Figure 4 shows that the method using hierarchic numbers calibrated from FCI calculations with the VXZ basis set works quite well, giving confidence to the CBS-extrapolated values based on the d-AVXZ basis. If the CBS extrapolated d-AVXZ correlation energies are added to the CBS/HF/d-AVXZ (X ) D-6) ones, the interaction energy falls short (more negative) of the target value at the equilibrium geometry of 2.97 Å by

TABLE 3: Dispersion Coefficientsa Ib coefficient

with ∆c

without ∆c

C6 C8 C10 C11 C12 C13 C14 C15 C16

1.4646 14.112 178.13 -76.7e 3093 -3806.0e 72016 -171000.0e 2276994

1.4541 14.011 176.85 -76.7e 3071 -3806.0e 71499 -171000.0e 2260636

IIb

otherd

1.4872 14.330 180.89 -76.7e 3141 -3806.0e 73128 -171000.0e 2312168

1.460977837725 14.11785737 183.691 075 -76.7e 3372 -3806.0 85340 -171000.0 2860000

a In atomic units. For coefficient Cn, such units are Ehan0. b This work, calculated with R0 ) 4.353922 a0. c ∆ stands for the VX/Zd-AVXZ ∆E(T)fFCI correction; see the text. d The errors of the consecutive Cn are65 0.000000000002, 0.00000002, 0.000001, 0.04, 1, 1, 50, 1000, and 10000. e Taken from ref 65.

334 mK. This is visible from Figure 2, where the curves obtained here are compared among themselves and with the accurate one of Jeziorska et al.65 Interestingly, the deviations at the geometries calculated in the present work do not exceed 50 cK or so. Given the still appreciable differences, it seems unjustified discussing other corrections. Although the deviation from the commonly accepted well depth in method II amounts to only ∼3% of the well depth, this is a somewhat frustrating result in the sense that a much better estimate has been obtained with the more-cost-effective method I. However, it is interesting that good results have been obtained despite the fact that the extrapolated HF and correlation energies were contaminated by large BSSEs. Thus, although the FCI/d-AVXZ calculations constituted a major part of computer costs in the present work, such calculations by themselves led to a significantly less accurate potential than the one obtained by method I. Unfortunately, the authors of ref 15 did not report the absolute values of the FCI energies calculated with the d-AVXZ basis, which might allow an improvement of our estimate of the correlation energy. Indeed, such calculations are unaffordable and would probably not add much to the conclusions extracted here. However, the convergence properties in ref 12 at the CCSD(T) level leave us the optimistic view that a closer agreement would be obtained at the CBS/FCI(p,h) NCP level. The message to take is that the performance of the d-AVXZ basis is relatively modest, a fact already noticed in

8512

J. Phys. Chem. A, Vol. 114, No. 33, 2010

Varandas

TABLE 4: Interaction Energies (in Kelvin) for the Helium Dimer at Three Bond Distances method

R ) 4.0 a0

5.6 a0

7.0 a0

Korona et al. Jeziorska et al.b Komasac Komasad Andersone Andersonf van Mourik and Dunningg van de Bovenkamp and van Duijnveldth Gdanitzi Klopperj Cencek et al.k Patkowski et al.l Jeziorska et al.m Varandasn Varandaso this workp

291.64 ( 0.87

-11.059 ( 0.03 -11.000 ( 0.011 -10.798 -10.981 -10.98 ( 0.02 -10.998 ( 0.005 -11.004 ( 0.03 -10.99 ( 0.02 -10.980 ( 0.004 -10.99 ( 0.02 -11.009 ( 0.008 -11.0037 ( 0.0031 -11.006 ( 0.004 -11.003 ( 0.03 -11.003 ( 0.03 -11.015(-11.005) ( 0.03

-4.629 ( 0.03

a

292.784

293.496 292.72 ( 0.2 292.75 ( 0.01 292.6 ( 0.3 292.54 ( 0.04 292.58 292.20 ( 0.4 292.31 ( 0.4 292.4 ( 0.3

-4.583

-4.620 ( 0.002 -4.619 ( 0.007 -4.62 -4.649 ( 0.02 -4.637 ( 0.02 -4.594 ( 0.02

a From symmetry-adapted perturbation theory.38 b From symmetry-adapted perturbation theory,65 cited in ref 15. c Variational upper bounds using exponentially correlated Gaussian functions.42 d Variational upper bound.43 e From quantum Monte Carlo36 calculations. f From quantum Monte Carlo37 calculations. g Extrapolated CCSD(T) values obtained by combining earlier results with nonextrapolated CCSDT and FCI energies.33 h MRCI results from ref 32. i From r12-MR-ACPF calculations.31 j Extrapolated results from CCSD(T) and FCI calculations in orbital bases.39 k From supermolecular Gaussian geminal calculations.44 l From supermolecular Gaussian geminal calculations.15 m From the analytical fit in ref 65, with well depth calculated at the minimum (2.967 Å). n From ref 12 by direct extrapolation of the total interaction energy and the corrections to FCI from ref 15; δE(T)fFCI ) -(1.883 ( 0.011) K at 4a0; δE(T)fFCI ) -(0.08084 ( 0.00031) K at 7a0. o Result obtained12 from separate extrapolations of the NCP and CP raw CCSD(T) energies with the GUSTE (5, 6) method.48 Results corrected for FCI as in footnote n. p From method I of the present work. See the text for the best estimate (in parentheses), corrections at other distances, and also for method II.

attractive dispersion forces. Any appearance of a well at the HF level is mostly due to BSSE. Note further that the size of the latter can be up to an order of magnitude larger than the depth of the vdW well in He2, in agreement with similar findings3,24 for other rare gas dimers. Figure 8 compares the calculated correlation energy with the dispersion series expansion

Vcor(R) ) -

∑

Cnχn(R)R-n

(12)

n)6,8,10-16

where χn(R) are dispersion damping functions, and the series has here been truncated at C16. In a similar analysis for Ne2, Bytautas and Ruedenberg24 have taken χn(R) as modified Tang-Toennies forms,66 whereas we utilize our own damping functions67 Figure 6. A comparison of the analytic potential energy curves of the present work for the helium dimer with the one of ref 65. This is indistinguishable from the curves of methods I and II in the main plot, with the differences [in centikelvin (cK)] shown in the inset; the open circles (triangles) refer now the difference at the calculated points with (without) consideration of the ∆EVX/Zd-AVXZ correction (this is abbreviated (T)fFCI as ∆). Note that the curves refer only to the fit itself and do not contain the corrections mentioned at the end of section 3.

ref 15, where the authors estimated the uncertainty at R ) 5.6 a0 to be as large as 64 mK. The use of midbond functions may partly alleviate the problem,15 but likely at an increased cost. Of course, a test of such basis sets is outside of the scope of the present work. As in most ab initio studies on the noble gas dimers, a major goal here has been to predict as accurately as possible the binding energy and the equilibrium distance. We may now aim at knowing whether the supermolecular approach can yield a correlation potential curve that decays in the long-range region as the dispersion expansion. This is crucial to obtain an accurate curve over the full configuration space. Note that the induction energy is rather small, and hence, the minimum is dictated mostly from a balance between the HF repulsion and the

χn(R) ) [1 - exp(-AnR/F - BnR2 /F2)]n

(13)

where An and Bn in eq 13 are auxiliary functions68,69 defined by

An ) R0n-R1

(14)

Bn ) β0 exp(-β1n)

(15)

and R0, β0, R1, and β1 are universal dimensionless coefficients for all isotropic interactions; R0 ) 16.36606, R1 ) 0.70172, β0 ) 17.19338, and β1 ) 0.09574. Note that the scaling parameter F is defined by F ) 5.5 + 1.25R0, where R0 ) 2(〈rX2 〉1/2 + 〈rY2 〉1/2) and 〈rX2 〉 is the expectation value of the squared radii for the outermost electrons in atom X. For helium, this has the value of70 R0 ) 4〈rX2 〉1/2 ) 4.353922 a0. The only fitting parameter is C6, with all other dispersion coefficients predicted from relations available in the literature. For n ) 8 and 10, we use69,71

Cn ) κnRa(n-2)/2 0

(16)

with κ8 ) 1, κ10 ) 1.13, and a ) 1.31. The others have been estimated with66,69,72-75

The Pair Potential of Helium Revisited

J. Phys. Chem. A, Vol. 114, No. 33, 2010 8513

( )

Cn+4 ) κn+4Cn-2

Cn+2 Cn

2

(17)

where75 κ12 ) 1.028 and κ14 ) 0.975. Since we are not aware of any other value, a fair assumption is κ16 ) (κ12 + κ14)/2. All points beyond the equilibrium geometry (with a weight of 108) carried an equal weight of 104 in the least-squares fitting procedure, while the others carried a unit weight. The fit included the coefficients C11, C13, and C15 as reported in ref 65. The advantage of such a fitting approach is three-fold. First, the well-established linear dependence of the coefficients makes it impossible to determine in a unique manner the various terms in eq 13. In fact, it has been observed24 that even the two coefficients C8 and C10 cannot be determined unambiguously, if the values of the long-range potential have possible errors of up to 0.01 µEh, which is the type of accuracy that one often gets from ab initio calculations. Moreover, extensive trade-offs are expected to occur in polyatomic systems, where terms in R-7 and R-9 can be part of the expansion. Second, there has been extensive backup work66,69,71-75 asserting the reliability of the above semiempirical rules. Third, it is a one-parameter fit, although slightly more elaborate but also problematic, variants (with k16 or a, or so on in the optimization list) may be open to experimentation. It may be argued that C6 dominates at asymptotic distances and hence should not be fixed as far inside of the potential as the vdW minimum. This is not the case since all coefficients are interrelated and the dispersion suitably damped to account for charge-overlap effects. In fact, this may be another advantage of the present fitting scheme as it may be calibrated from energies for relatively short distances. Indeed, this has been our choice to avoid the relatively large uncertainties in the data at large values of R (see the inset plot of Figure 7). A further remark to note is that removing C11, C13, and C15 from the fit leads to a value of C6 that falls within 1% of the one given in Table 3. However, such coefficients may play a subtle role in getting the balance of the various terms right and hence

Figure 8. Calculated raw correlation energies showing that the d-AVXZ basis sets for X ) D, T produce a rather strong BSSE (the same holds for the HF energies as shown in Figure 2). The same conclusion is even more clear from the inset plot. In this, the intercepts at the origin indicate the predicted values of C6. Note that the CBS-extrapolated correlation (dispersion) curves, shown free of symbols, are indistinguishable (method I, solid; method II, dash-dot).

should be kept whenever available if the expansion is truncated as high up as in the case here reported. As shown in Figures 7 and 8, the agreement with the accurate dispersion interaction of Jeziorska et al.65 is quite good over the entire configuration space. Quantitatively, eq 12 predicts at R ) 2.97 Å for method I (II) a value of -19.996 (-20.308) K, which compares with the value of -20.080 K calculated from the curve of ref 65. Also notable is the agreement obtained for the various dispersion coefficients, with the unsigned errors obtained from a fit to the energies of method I (similarly for II) varying for increasing even values of n as