# Towards standard methods for benchmark quality ab initio thermochemistry — W1 and W2 theory

## I Introduction

Thermochemical data such as molecular heats of formation are among the most crucial quantitative chemical data. Thanks to great progress made in recent years in both methodology and computer technology, a broad range of empirical, semiempirical, density functional, and ab initio schemes now exist for this purpose (for a recent collection of reviews, see Ref.).

At present, only ab initio-based methods can claim ‘chemical accuracy’ ($`\pm `$1 kcal/mol). The most popular such schemes are undoubtedly the G2 and G3 theories of Pople and coworkers (which are based on a combination of additivity approximations and empirical corrections applied to relatively low-level calculations), followed by the CBS-Q and CBS/APNO methods which are intricate combinations of extrapolation and empirical correction schemes. With the exception of CBS/APNO (which allows for 0.5 kcal/mol accuracy, on average, but is restricted to first-row compounds), all these schemes allow for mean absolute errors of about 1 kcal/mol, although errors for some individual molecules (e.g. SO<sub>2</sub>, SiF<sub>4</sub>) can be much larger (e.g. about 8–12 kcal/mol for SiF<sub>4</sub> using G2 theory, and 4 kcal/mol using G3 theory).

In fact, many of the experimental data in the ”enlarged G2 set” employed in the parametrization of several of these methods (notably G3 theory and several of the more recent density functional methods) themselves carry experimental uncertainties larger than 1 kcal/mol.

The aim pursued in the present work is a more ambitious one than chemical accuracy. In light of the prevalent use of kJ/mol units in the leading thermochemical tables compendia (JANAF and CODATA), we shall arbitrarily define a mean absolute error of one such unit, i.e. 0.24 kcal/mol, as ‘calibration accuracy’ — with the additional constraint that no individual error be larger than the ‘chemical accuracy’ goal of 1 kcal/mol.

One of us has recently shown that this goal is achievable for small polyatomics using present technology. The approach followed employed explicit treatment of inner-shell correlation, coupled cluster calculations in augmented basis sets of $`spdfg`$ and $`spdfgh`$ quality, and extrapolation of the valence correlation contribution to the atomization energy using formulas based on the known asymptotic convergence behavior of pair correlation energies. In this manner, total atomization energies (TAE<sub>e</sub>) of about 15 first-row diatomics and polyatomics for which experimental data are known to about 0.1 kcal/mol could be determined to within 0.25 kcal/mol on average without any empirical parameters. (Upon introducing an empirical correction for A–N bonds, this could be improved to 0.13 kcal/mol, clearly within the target.) In fact, using this method, an experimental controversy concerning the heat of formation of gaseous boron — a quantity that enters any ab initio or semiempirical calculation of the heat of formation of any boron compound — could be resolved by a benchmark calculation of the total atomization energy of BF<sub>3</sub>.

Benchmark studies along similar lines by several other groups (e.g. those of Helgaker, Bauschlicher, Dunning) point in the same direction. Among those, Bauschlicher was the first to suggest that the inclusion of scalar relativistic corrections may in fact be essential for accurate results on second-row molecules.

High-accuracy results on second-row compounds can only be achieved in this manner — as has been shown repeatedly — if high-exponent $`d`$ and $`f`$ functions are added to the basis set. As shown by one of us, these ‘inner shell polarization functions’ address an SCF-level effect which bears little relationship to inner-shell correlation, and actually dwarfs the latter in importance (contributions as high as 10 kcal/mol having been reported).

All these approaches carry a dual disadvantage: their extravagant computational cost and their reliance on the quantum chemical expertise of the operator.

The target of the present study was to develop computational procedures that meet the following requirements:

* they should have mean absolute errors on the order of 0.25 kcal/mol or less, and problem molecules (if any) should be readily identifiable;
* the method should be applicable to at least first-and second-row molecules;
* it should be robust enough to be applicable in a fairly ‘black-box’ fashion by a nonspecialist;
* it should rely as little as possible (preferably not at all) on empirical parameters, empirical additivity corrections, or other ‘fudges’ derived from experimental data;
* relatedly, it should explicitly include all the physical effects that are liable to affect molecular binding energies of first-and second-row compounds, rather than rely upon absorbing them in empirical parametrization;
* last but not least, it should be sufficiently cost-effective that a molecule the size of, e.g., benzene should be treatable on a workstation computer.

In the course of this work, we will present two schemes which we shall denote W1 and W2 (for Weizmann-1 and Weizmann-2) theories. W2 theory yields about 0.2 kcal/mol (or better) accuracy for first-and second-row molecules with up to four heavy atoms, and involves no empirical parameters. W1 theory is applicable to larger systems (we shall present benzene and trans-butadiene as examples), yet still yields a mean absolute error of about 0.30 kcal/mol and includes only a single, molecule-independent, empirical parameter which moreover is derived from calculated rather than experimental results.

## II Computational details

Most electronic structure calculations reported in this work were carried out using MOLPRO 97.3 and MOLPRO 98.1 running on a Silicon Graphics (SGI) Octane workstation and on the SGI Origin 2000 of the Faculty of Chemistry. The full CCSDT (coupled cluster with all connected single, double, and triple substitutions) calculations were carried out using ACES II running on a DEC Alpha 500/500 workstation.

SCF and valence correlation calculations were carried out using correlation consistent polarized $`n`$-tuple zeta (cc-pV$`n`$Z, or V$`n`$Z for short) ($`n`$=D, T, Q, 5, 6) and augmented correlation consistent polarized $`n`$-tuple zeta (aug-cc-pV$`n`$Z, or AV$`n`$Z for short) ($`n`$=D, T, Q, 5, 6) basis sets of Dunning and coworkers. The maximum angular momentum parameter $`l`$, which occurs in the extrapolation formulas for the correlation energy, is identified throughout with the $`n`$ in V$`n`$Z and AV$`n`$Z. Except for the calculation of the electron affinity of hydrogen, regular V$`n`$Z basis sets were used throughout on hydrogen atoms.

Most valence correlation calculations were carried out using the CCSD (coupled cluster with all single and double substitutions) and CCSD(T) (i.e. CCSD followed by a quasiperturbative estimate of the effect of connected triple excitations) electron correlation methods. The CCSD(T) method is known to be very close to an exact solution within the given one-particle basis set if the wave function is dominated by dynamical correlation.

Where possible, imperfections in the treatment of connected triple excitations were estimated by comparing with full CCSDT calculations. The effect of connected quadruple and higher excitations were estimated by small basis set FCI (full configuration interaction) calculations — which represent exact solutions with a finite basis set.

Inner-shell correlation contributions were evaluated by taking the difference between valence-only and all-electron CCSD(T) calculations in special core-correlation basis sets. For first-row compounds, both Dunning’s ACVQZ (augmented correlation consistent core-valence quadruple zeta) basis set and the Martin-Taylor (MT) core correlation basis sets were considered; for second-row compounds only the MT basis sets. The latter are generated by completely decontracting a CV$`n`$Z or ACV$`n`$Z basis set, and adding one tight $`p`$ function, three high-exponent $`d`$ functions, two high-exponent $`f`$ functions, and (in the case of the MTv5z basis set) one high-exponent $`g`$ function to the basis set. The additional exponents were derived from the highest ones already present for the respective angular momenta, successively multiplied by 3.0. The smallest such basis set, MTvtz (based on VTZ) is also simply denoted MT.

Scalar relativistic corrections were calculated at the ACPF (averaged coupled pair functional) level as expectation values of the first-order Darwin and mass-velocity terms. An idea of the reliability of this approach is given by comparing a very recent relativistic (Douglas-Kroll) coupled cluster calculation of the relativistic contribution to TAE\[SiH<sub>4</sub>\], $``$0.67 kcal/mol, with the identical value of $``$0.67 kcal/mol obtained by means of the present approach. For GaCl, GaCl<sub>2</sub>, and GaCl<sub>3</sub> — where relativistic effects are an order of magnitude stronger than even in the second-row systems considered here — Bauschlicher found that differences between Douglas-Kroll calculations and the presently followed approach amounted to 0.12 kcal/mol or less on the binding energy.

Spin-orbit coupling constants were evaluated at the CASSCF-CI level using the $`spdf`$ part of the MTav5z basis set. (For a recent review of the methodology involved, see Ref..)

Density functional calculations for the purposes of obtaining certain reference geometries and zero-point energies were carried out using the Gaussian 98 package. Both the B3LYP (Becke 3-parameter-Lee-Yang-Parr) and B3PW91 (Becke 3-parameter-Perdew-Wang-1991) exchange-correlation functionals were considered.

Most geometry optimizations were carried out at either the CCSD(T)/VQZ+1 or the B3LYP/VTZ+1 (in some cases B3PW91/VTZ+1) levels of theory, where the notation V$`n`$Z+1 indicates the addition to all second-row atoms of a single high-exponent $`d`$-type ‘inner polarization function’ with an exponent equal to the highest $`d`$ exponent in the Dunning V5Z basis set. In the past this was found to recover the largest part of the effects of inner polarization on geometries and vibrational frequencies. (We note that for molecules consisting of first-row atoms only, the V$`n`$Z+1 basis sets are equivalent to regular V$`n`$Z basis sets.)

Past studies of the convergence behavior of the SCF energy have shown it to be very well described by a geometric extrapolation of the type first proposed by Feller, $`A+B/C^l`$. Clearly, for this purpose a succession of three SCF/AV$`n`$Z basis sets is required.

For the valence correlation CCSD and (T) energies, two extrapolation formulas were considered. The first, $`A+B/(l+1/2)^\alpha `$, was proposed by Martin — the philosophy being that using the extrapolation exponent as an adjustable parameter would enable inclusion of higher-order terms in the asymptotic expansion

$$A/(L+1)^3+B/(L+1)^4+C/(L+1)^5+\mathrm{}$$
(1)   
while the denominator shift of 1/2 was a compromise — for identification of the $`l`$ in cc-pV$`l`$Z with L — between hydrogen and nonhydrogen atoms. The second formula, simply $`A+B/l^3`$, was proposed by Helgaker and coworkers — where $`l`$ was identified with $`L1`$ throughout. Halkier et al. already noted that in terms of the extrapolated energy using $`A+B/(l+C)^D`$, the parameters C and D were very strongly coupled, and that it only made sense to vary one of them.

The combination of treatments for SCF, CCSD valence correlation, (T), imperfections in the T treatment, and connected quadruple and higher excitations is compactly denoted here by W\[p5;p4;p3;p2;p1\], in which p1 denotes the basis sets involved in the SCF extrapolation, p2 the basis sets involved in the CCSD extrapolation, p3 those in the (T) extrapolation (which may or may not be different from p2), p4 (if present) the basis sets used in correcting for imperfections in the treatment of connected triple excitations, and p5 (if present) those involved in evaluating the effect of connected quadruple and higher excitations. If any of the p’s consists of a single index, a simple additivity approximation is implied; two indices denote a two-parameter extrapolation of the type $`A+B/l^3`$, while three indices indicate a three-parameter extrapolation of the type $`A+B/(l+1/2)^\alpha `$ in the case of correlation contributions, and $`A+B/C^l`$ in the case of SCF contributions. For example, the level of theory used in the previous work of Martin and Taylor would be W\[TQ5;TQ5;TQ5\] in the present notation, while W\[D;Q;TQ5;TQ5;TQ5\] indicates W\[TQ5;TQ5;TQ5\]+CCSDT/AVQZ-CCSD(T)/AVQZ+FCI/AVDZ-CCSDT/AVDZ.

## III Atomic electron affinities as a ‘litmus test’

The electron affinities of the first-and second-row atoms have often been used as benchmarks for high-level electronic structure methods (see e.g. the introductions to Refs. for reviews). Because electron affinities involve a change in the number of electrons correlated in the system, they are very taxing tests for any electron correlation method; in addition, they involve a pronounced change in the spatial extent of the wave function, making them very demanding in terms of the basis set as well.

Until recently, three of the first-and second-row atomic electron affinities were imprecisely known experimentally (B, Al, and Si): this situation was changed very recently by high-precision measurements for recent experiments for B, Al, and Si.

The approach we have chosen here for the SCF and valence correlation components is summarized in our notation as W\[$`n`$,Q,56,56,Q56\] for the first-row atoms, and W\[$`n`$,Q,Q5,Q5,TQ5\] for the second-row atoms. The effect of inner-shell correlations was assessed at the CCSD(T)/MTav5z level, while Darwin and mass-velocity corrections were evaluated at the ACPF/MTav5z level. Finally, spin-orbit splittings were calculated at the CASSCF-CI level with the $`spdf`$ part of a MTav5z basis set. (For technical reasons, the $`h`$ functions were omitted in both the scalar relativistic and spin-orbit calculations, as were the $`g`$ functions in the latter.)

Our best computed results are compared with experiment in Table I, where results from recent calibration studies are also summarized.

Agreement between computed and observed values can be described without reservation as excellent: the mean absolute error amounts to 0.0009 eV. The fact that this accuracy is obtained systematically and across the board strongly suggest that the ’right result was obtained for the right reason’. Upon eliminating the corrections for imperfections in CCSD(T), i.e. restricting ourselves to W\[56,56,Q56\] for first-row atoms and W\[Q5,Q5,TQ5\] for second-row atoms, the mean absolute error increases by an order of magnitude to 0.009 eV, i.e. about 0.2 kcal/mol. As we shall see below, this is essentially the type of accuracy we can obtain for molecules without corrections for CCSD(T) imperfections.

The importance of Darwin and mass-velocity corrections increases, as expected, with increasing Z, and its contribution becomes quite nontrivial for atoms like Cl. It is therefore to be expected that, e.g., in polar second-row molecules like ClCN or SO<sub>2</sub> they will contribute substantially to TAE as well.

The importance of inner-shell correlation effects is actually largest for Al, because of the small gap between valence and sub-valence orbitals in the early second-row elements.

Table II compares the convergence behavior of the extrapolated valence correlation contributions as a function of the largest basis set used, both using the Martin three-term and Helgaker two-term formulas. While both formulas appear to give the same answer if the underlying basis sets are large enough, the two-term formula is by far the more stable towards reduction of the basis sets used in the extrapolation. Since the use of W\[Q56,Q56,Q56\] is hardly an option for molecules, the two-term formula appears to be the formula of choice.

Following the suggestion of a referee, we have considered (Table II) the performance of some other extrapolation formulas for the valence correlation energy. As a point of reference, we have taken an “experimental valence correlation contribution to EA”, which we derived by subtracting all computed contributions other than the valence correlation energy from the best experimental EAs. While some residual uncertainties may remain in some of the individual contributions, these should be reliable to 0.001 eV on average.

As seen in Table II, performance of the geometric series extrapolation $`A+B/C^n`$ is outright poor: in fact, for extrapolation from AV{D,T,Q}Z results the error is twice as large as that caused by not extrapolating at all. If AV{T,Q,5}Z basis sets are used, mean absolute error drops to 0.015 eV, which is still an order of magnitude larger than for the $`A+B/l^3`$ extrapolation, and only slightly better than not extrapolating at all. Finally, for AV{Q,5,6}Z basis sets, the error is three times smaller than complete omission of extrapolation, but three times larger than that of using any of the following formulas: $`A+B/l^3`$ , $`A+B/l^C`$ , or $`A+B/(l+1/2)^4+C/(l+1/2)^6`$ . All three of the latter yield a mean absolute error of about 0.001 eV, on about the same order of accuracy as the reference values. For the smallest basis set series AV{D,T,Q}Z, the mixed exponential-Gaussian extrapolation $`A+B/\mathrm{exp}(l1)+C/\mathrm{exp}((l1)^2)`$ represents a very substantial improvement over $`A+B/C^l`$, and actually exhibits slightly better performance than $`A+B/l^3`$. For the AV{T,Q,5}Z series which is of greatest interest here, the Halkier et al. $`A+B/l^3`$ formula by far outperforms the other formulas considered.

In short, it appears to be established that the two-term formula of Helgaker and coworkers is the extrapolation method of choice overall, with the Martin three-term formulas the second-best choice provided basis sets of AV{T,Q,5}Z quality are used. The mixed exponential-Gaussian formula performs slightly better than $`A+B/l^3`$ if only AV{D,T,Q}Z basis sets are used (see however Section VI.C below).

Computed spin-orbit contributions to the electron affinities are compared in Table III to values obtained from observed fine structures. While small deviations appear to persist, these may at least in part be due to higher-order spin-orbit effects which were neglected in the calculation rather than to deficiencies in the electronic structure treatment. At any rate, to the accuracy relevant for our purpose (establishing spin-orbit corrections to molecular binding energies) it appears to be immaterial whether the computed or the experimentally derived values are used.

Finally, the convergence of the SCF component is so rapid that it appears to be essentially irrelevant which extrapolation formula is used — the amount bridged by the extrapolation is on the order of 0.0001 eV.

## IV Results for molecules

Since application of electron correlation methods more elaborate than CCSD(T) would be well-nigh impossible for molecules of practical size, we have restricted ourselves to W\[Q5;Q5;TQ5\] and W\[TQ5;TQ5;TQ5\].

Inner-shell correlation contribution, as well as scalar relativistic corrections, were initially computed with the largest basis sets practicable – in most cases ACV5Z or MTavqz (see Table IV for details).

From a prerelease version of a re-evaluation of the experimental data in the G2/G3 set at the National Institute for Standards and Technology (NIST), we have selected 28 first-and second-row molecules which satisfy the following criteria: (a) the uncertainty in the experimental total atomization energy TAE is on the order of 0.25 kcal/mol or better; (b) the molecules are not known to exhibit severe multireference effects; (c) anharmonic vibrational zero-point energies are available from either experiment or high-level ab initio calculations (see footnotes to Table IV for details).

Geometries were optimized at the CCSD(T)/VQZ+1 level, and to all second-row atoms a complement of two tight d and one tight f function were added in every basis set to ensure saturation in inner-shell polarization effects. In all cases, the exponents were derived as even-tempered series $`\alpha \beta ^n`$ with $`\beta =3.0`$ and $`\alpha `$ the highest exponent already present for that angular momentum.

Computed (W\[Q5;Q5;TQ5\]) and observed results are compared in Table IV. The excellent agreement between theory and experiment is immediately apparent: in many cases, the computed results fall within the already quite narrow experimental error bars. Over the entire sample of molecules, the mean absolute error is 0.24 kcal/mol, with the largest errors being about 0.6 kcal/mol (O<sub>2</sub> and F<sub>2</sub>). Restricting our sample to first-row molecules only, we find a mean absolute error of 0.24 kcal/mol, which however gets reduced to 0.17 kcal/mol (maximum error 0.39 kcal/mol for N<sub>2</sub>) upon elimination of F<sub>2</sub>, NO, and O<sub>2</sub> as having known appreciable nondynamical correlation effects. Over the subset of second-row molecules in our sample MAE is 0.23 kcal/mol (maximum error 0.44 kcal/mol for H<sub>2</sub>S); upon elimination of H<sub>2</sub>S and SO<sub>2</sub> this is lowered to 0.20 kcal/mol.

It should be noted that these MAEs are comparable to those found by Martin and Taylor for a sample of first-row molecules, yet unlike their study no correction for N-containing bonds is required here.

The possibility that the errors in F<sub>2</sub>, NO, and O<sub>2</sub> are actually due to residual basis set incompleteness and/or that the excellent agreement with experiment for the other molecules is actually due to an error compensation involving deficiencies in the predicted basis set limit, was examined by carrying out W\[56;56;Q56\] calculations for H<sub>2</sub>O, F<sub>2</sub>, NO, O<sub>2</sub>, N<sub>2</sub>, HF, and CO. As seen in Table V, the predicted basis set limits do not differ materially from their W\[Q5;Q5;TQ5\] counterparts, strongly suggesting that the latter expression in fact does reach the basis set limit and that the residual errors are largely due to imperfections in the CCSD(T) method.

While molecules liable to exhibit such errors are readily identifiable from inspection of the largest coupled cluster amplitudes or evaluation of the $`𝒯_1`$ diagnostic, an even simpler criterion is apparently offered by the ratio TAE\[SCF\]/TAE\[SCF+val.corr.\]. In “well-behaved” molecules such as CH<sub>4</sub> and H<sub>2</sub>O, the SCF component makes up upwards of two-thirds of the binding energy, while in NO and in O<sub>2</sub> it makes up no more than a third and a fifth, respectively, of the total and F<sub>2</sub> is actually metastable towards dissociation at the SCF level. While for some molecules of this variety we actually obtain excellent results (e.g. ClF), this may be due to error compensation or to the binding energies being fairly small to begin with.

Further inspection of Table IV reveals that some of the ‘negligible’ contributions are in fact quite significant at the present precision level: for instance, Darwin and mass-velocity contributions in SO<sub>2</sub> amount to -0.71 kcal/mol (for SiF<sub>4</sub> a somewhat extravagant -1.88 kcal/mol was found), while atomic spin-orbit splitting in such compounds as Cl<sub>2</sub>, ClF, and SO<sub>2</sub> amounts to -1.68, -1.23, and -1.01 kcal/mol, respectively. Inner-shell correlation contributions of 2.36 (C<sub>2</sub>H<sub>4</sub>), 2.44 (C<sub>2</sub>H<sub>2</sub>), 1.68 (OCS), and 1.76 (ClCN) kcal/mol speak for themselves; interestingly (as noted previously), these effects on the whole do not seem to be more important in second-row than in first-row compounds.

Finally, we shall compare the performance of W\[TQ5;TQ5;TQ5\] and W\[Q5;Q5;TQ5\] (Table VI). In general, the results with the three-point valence correlation extrapolation are at best of the same quality as those with the two-point valence correlation extrapolation and in many cases agree less well with experiment. We therefor will use the two-point extrapolation exclusively henceforth.

## V W2 theory and its performance

Having established that our ’base level of theory’ can obtain the right results for the right reason, we shall now proceed to consider simplifications.

### A Inner-shell correlation

The use of the smaller MT basis set for the scalar relativistic contributions is found to have an effect of about 0.01 kcal/mol or less, with 0.02 kcal/mol being the largest individual cases. This approximation can therefore safely be made.

Using the same MT basis set for the core correlation contribution on average affects energetics by 0.03 kcal/mol, the largest individual effects being 0.07 kcal/mol for H2S, and 0.08 kcal/mol for OCS.

Even so, in fact, the core correlation calculations are quite CPU-time consuming, particularly for second-row compounds, due to the large number of electrons being correlated. Any further reduction would obviously be welcome — it should be noted that the MT basis set was developed not with efficiency, but with saturation (in the core-valence correlation energy) in mind. Further experimentation revealed that the tightest p, d, and f functions could safely be eliminated, but that further basis set reductions adversely affect the quality of the core correlation contribution computed. The reduced basis set shall be denoted as MTsmall, and in fact consists of a completely decontracted cc-pVTZ basis set with two tight $`d`$ and one tight $`f`$ functions added. Since this basis set only has about half the basis functions of the ACVQZ basis set per heavy atom, it represents a very significant CPU time savings (about 16 times) in a CCSD(T) calculation. The only molecule for which we see a substantial difference with the MT basis set is SO<sub>2</sub>, for which Bauschlicher and Ricca previously noted that the inner-shell correlation contribution is unusually sensitive to the basis set.

For the evaluation of the Darwin and mass-velocity corrections, differences with the larger MT basis set are less than 0.01 kcal/mol across the board.

A further reduction in CPU time for the core correlation contribution would have been achieved if MP2 or even CCSD calculations could be substituted for their CCSD(T) counterparts. However, as seen from Table VII, CCSD underestimates the CCSD(T) core correlation contributions for several molecules by as much as 50%. The behavior of MP2 is quite similar and the MP2–CCSD differences are substantially smaller than the (T) contribution, suggesting that it is the treatment of connected higher excitations that is the issue. Predictably, the largest (T) effects on the core correlation contribution occur in molecules where connected triple excitations are important for the valence binding energy as well, e.g. SO<sub>2</sub>, F<sub>2</sub>, Cl<sub>2</sub>, N<sub>2</sub>. Conversely, in CH<sub>3</sub> or CH<sub>4</sub>, which have quite small (T) contributions to the binding energy, CCSD does perform excellently for the core correlation contribution. In PH<sub>3</sub> and H<sub>2</sub>S, on the other hand, substantial errors in the core correlation are seen even as the (T) contribution to the valence correlation binding energy is quite small — it should be noted, however, that both the absolute inner-shell correlation energy and the (T) contribution to it are much more important in these second-row systems than in their first-row counterparts.

One may rightly wonder whether the inner-shell contributions are in fact converged at the CCSD(T) level. Unfortunately, if a more elaborate treatment is already impractical for the valence correlation, this would a fortiori be true for the inner-shell correlation. We did, however, carry out a CCSDT/MTsmall calculation on the N<sub>2</sub> molecule, which we chose as a representative case of failure of the CCSD approach for core correlation. The resulting CCSDT level core contribution, 0.87 kcal/mol, is only 0.05 kcal/mol larger than the CCSD(T) value of 0.82 kcal/mol, to be compared with 0.42 kcal/mol at the MP2 and 0.52 kcal/mol at the CCSD level. It cannot be ruled out a priori that connected quadruple and higher excitations might contribute to the inner-shell correlation energy. However, since apparently their importance for the valence correlation binding energy is very small (otherwise a treatment that completely ignores them would not yield the type of agreement with experiment found in this work), it seems unlikely that they would greatly contribute to the inner-shell correlation energy.

The “G3large” basis set used to evaluate, among other things, inner-shell correlation effects in G3 theory is still smaller than the MTsmall basis set, and its performance therefore is certainly of interest. Alas, in Table VII it is seen that in many cases it seriously overestimates the inner-shell correlation energy, almost certainly because of basis set superposition error which is apparently more of an issue for inner-shell correlation energies than for their valence counterparts. In G3 theory, the inner-shell correlation is evaluated at the MP2 level: hence the two errors cancel to a substantial extent.

### B Zero-point energy

Not in all cases is a complete anharmonic force field calculation feasible. We find in Table VIII that B3LYP/cc-pVTZ+1 harmonic zero-point energies scaled by 0.985 reproduce the rigorous anharmonic zero-point energies quite nicely. (The scaling factor is about halfway between what would be required for fundamentals, about 0.97, and harmonics, about 1.00.)

### C Separate extrapolation of CCSD and (T)

The (T) contribution makes up a relatively small part of the valence correlation energy, while its evaluation, in large basis sets and for systems with very many electrons, will dominate the CPU time. For instance, in a very recent study on SiF<sub>4</sub>, a CCSD(T) calculation with an AVQZ basis set on F and a VQZ+2d1f basis set on Si took 50h7 using MOLPRO on an SGI Octane workstation (768 MB of memory being allocated to the job), of which 41h30 were spent in the (T) step alone.

In addition, it was previously noted by Helgaker and coworkers that the (T) contribution appears to converge faster with the basis set than the CCSD correlation energy, for which reason they actually propose its separate evaluation in a smaller basis set. In the present case, we have considered extrapolating it from AVTZ+2d1f and AVQZ+2d1f results rather than AVQZ+2d1f and AV5Z+2d1f, respectively. In our adopted notation, this becomes W\[TQ;Q5;TQ5\].

As seen in Table VI, the difference in quality between W\[TQ;Q5;TQ5\] and W\[Q5;Q5;TQ5\] appears to be essentially negligible. This is an important conclusion, since it means that the largest basis set calculation to be carried out is only at the CCSD level, at a fraction of the cost of the full CCSD(T) counterpart — moreover it can be done using direct algorithms.

### D Protocol for W2 theory

The protocol obtained by introducing the succesful approximations given above will be denoted here as W2 (Weizmann-2) theory. Its steps consist of the following:

* geometry optimization at the CCSD(T)/VQZ+1 level, i.e. CCSD(T)/VQZ if only first-row atoms are present;
* zero-point energy obtained from a CCSD(T)/VTZ+1 anharmonic force field or, failing that, B3LYP/VTZ+1 frequencies scaled by 0.985 (vide infra);
* Carry out CCSD(T)/AVTZ+2d1f and CCSD(T)/AVQZ+2d1f single-point calculations;
* Carry out a CCSD/AV5Z+2d1f single-point calculation;
* the SCF component of TAE is extrapolated by $`A+B/C^l`$ from SCF/AVTZ+2d1f, SCF/AVQZ+2d1f, and SCF/AV5Z+2d1f results ($`l`$=3, 4, and 5, respectively);
* the CCSD valence correlation component is obtained from applying $`A+B/l^3`$ to CCSD/AVQZ+2d1f; and CCSD/AV5Z+2d1f valence correlation energies (l=4 and 5, respectively). It is immaterial whether this is done on total energies or on components to TAE;
* the (T) valence correlation component is obtained from applying $`A+B/l^3`$ results to CCSD(T)/AVTZ+2d1f and CCSD(T)/AVQZ+2d1f values for the (T) contribution. It is again immaterial whether this is done on total energies or on components to TAE;
* core correlation computed at CCSD(T)/MTsmall level
* scalar relativistic corrections (and, if necessary, spin-orbit splittings) computed at ACPF/MT level. To save CPU time, this can be combined into a single job with the previous step.

On a typical workstation at the time of writing (e.g. an SGI Octane with 1 GB of RAM and 2$`\times `$18 GB external disks) its applicability range would be about three heavy atoms in $`C_{2v}`$ symmetry, although the main limiting factor would be disk space and larger systems could be treated if a direct CCSD code were available.

We will illustrate the CPU time savings made in the W2 approach, compared to our most rigorous calculations, using two examples: a first-row diatomic (CO) and a second-row molecule (OCS). Using MOLPRO on an SGI Octane workstation, the most accurate calculations reported in this work (Table IV) required 21h36 for CO and no less than 362h12 for OCS. W2 theory yields essentially identical results at a cost of 1h12 (CO) or 13h42 (OCS) — a reduction by a factor of 20–30 which is typical of the other molecules.

## VI W1 theory and its performance

In an effort to obtain a method that is applicable to larger systems, we introduce a few further approximations. Relevant results can be found in Table VI.

### A Use of density functional reference geometries

B3LYP/cc-pVTZ+1 geometries are close enough to their CCSD(T)/VQZ+1 counterparts that their use does not cause a major effect on the final computed result. There is one notable exception to this rule for the molecules considered here: Cl<sub>2</sub>, for which the B3LYP/VTZ+1 bond distance of 2.0130 Å is quite different from its CCSD(T)/VQZ+1 counterpart, 1.9972 Å. (The experimental value is 1.987<sub>9</sub> Å.) While B3LYP and B3PW91 on the whole tend to produce essentially identical geometries and harmonic frequencies, it has been argued previously that the B3PW91 functional may be somewhat more reliable for systems with high electron density; and in fact, the B3PW91/VTZ+1 bond distance of 1.9912 Å is much closer to the CCSD(T)/VQZ+1 and experimental value.

Even so, the use of a B3LYP/VTZ+1 reference geometry still does not affect the computed $`D_e`$ by more than 0.14 kcal/mol. We conclude that CCSD(T) geometry optimizations, which will become fairly costly for larger molecules, can safely be replaced by B3LYP calculations, which can also serve for obtaining zero-point energies.

### B Further reduction of basis set sizes

The obvious first suggestion would be to also carry out the CCSD extrapolation using smaller basis sets, i.e. W\[TQ,TQ,DTQ\].

The effect for the SCF component of TAE is very small on condition that the extrapolation is carried out not on the individual total energies but rather on the computed SCF components of TAE themselves. Clearly error compensation occurs between the molecule and the constituent atoms.

The effect on the valence correlation component, unfortunately, is rather more significant. Over the ‘training set’, MAE rises to 0.37 kcal/mol even after SO<sub>2</sub> (which clearly is a pathological case here) has been eliminated. Aside from the latter, eliminating systems with mild nondynamical correlation does not lead to any significant reduction in MAE. Also noteworthy is that, on average, the binding energies appear to be somewhat overestimated: this is easily explained from the fact that basis sets like AVTZ+2d1f are not quite saturated in the radial correlation energy either, and that therefore the TAE\[AVQZ+2d1f\]$``$TAE\[AVTZ+2d1f\] gap will be an overestimate of the TAE\[L=4\]$``$TAE\[L=3\] gap.

### C Use of empirical extrapolation exponents

Truhlar considered the use of $`L`$-extrapolation formulas with empirical exponents, carried out from cc-pVDZ and cc-pVTZ calculations, as an inexpensive alternative to very large basis set calculations.

We will investigate here a variant of this suggestion adapted to the present framework. The valence correlation component to TAE will indeed be extrapolated using the formula $`A+B/l^\beta `$, in which $`\beta `$ is now an empirical parameter — we will denote this W\[Q5;Q5;TQ5\]$`\beta `$ and the like.

We then add in all the further corrections (core correlation, scalar relativistics, spin-orbit) that occur in W2 theory, and try to determine $`\beta `$ by minimizing MAE with respect to the experimental TAE values for our ‘training set’. Not surprisingly, for W\[Q5;Q5;TQ5\]$`\beta `$ this yields an optimum exponent ($`\beta `$=2.98) which differs insignificantly from the ‘ideal’ value of 3.0. Alternatively, $`\beta `$ could be optimized for the best possible overlap with the W\[56;56;Q56\] results: in fact, the same conclusion is obtained, namely that making $`\beta `$ an empirical parameter does not improve the quality of the results.

For W\[TQ;TQ;DTQ\] however, the situation is rather different. The optimum exponent $`\beta `$ is found to be 3.18 if optimized against the experimental TAE values, and 3.16 (insignificantly different) if optimized against the W\[Q5;Q5;TQ5\] results. In both cases, MAE drops to 0.30 kcal/mol, and on average no more overestimation occurs.

W\[TQ;TQ;DTQ\]3.18 represents a significant savings over W2 theory. Its time-determining step in molecules with many electrons will be the evaluation of the parenthetical triples in the AVQZ+2d1f basis set — their elimination would be most desirable.

A natural suggestion would then be W\[DT;TQ;DTQ\]$`\beta `$. Optimization of $`\beta `$ against the experimental TAE values yields $`\beta `$=3.26; the not greatly different $`\beta `$=3.22 is obtained by minimization of the deviation from W\[Q5;Q5;TQ5\] results for the training set. Since the latter does not explicitly depend on experimental results and therefore minor changes in the computational protocol do not require recalculation for the entire ‘training set’, we will opt for the latter alternative.

In either case, we obtain MAE=0.30 kcal/mol — for a calculation that requires not more than an AVTZ+2d1f basis set for the largest CCSD(T) calculation, and an AVQZ+2d1f basis set for the largest CCSD calculation. Again, the latter is amenable to a direct algorithm.

### D Protocol for W1 theory

We thus propose the following protocol for a computational level which we will call W1 (Weizmann-1) theory:

* geometry optimization at the B3LYP/VTZ+1 level (B3LYP/VTZ if only first-row atoms are present). Alternatively, the B3PW91 exchange-correlation functional may be preferable for some systems like Cl<sub>2</sub> — under normal circumstances, B3LYP/VTZ+1 and B3PW91/VTZ+1 should yield virtually identical geometries;
* zero-point energy obtained from B3LYP/VTZ+1 (or B3PW91/VTZ+1) harmonic frequencies scaled by 0.985;
* Carry out CCSD(T)/AVDZ+2d and CCSD(T)/AVTZ+2d1f single-point calculations;
* Carry out a CCSD/AVQZ+2d1f single-point calculation;
* the SCF component of TAE is extrapolated by $`A+B/C^l`$ from SCF/AVDZ+2d, SCF/AVTZ+2d1f, and SCF/AVQZ+2d1f components of TAE ($`l`$=2, 3, and 4, respectively)
* set $`\beta `$=3.22
* the CCSD valence correlation component is obtained from applying $`A+B/l^\beta `$ to CCSD/AVTZ+2d1f and CCSD/AVQZ+2d1f valence correlation energies (l=3 and 4, respectively). In both this and the next step, it is immaterial whether the extrapolation is carried out on components to the total energy or to TAE;
* the (T) valence correlation component is obtained from applying $`A+B/l^\beta `$ results to CCSD(T)/AVDZ+2d and CCSD(T)/AVTZ+2d1f values for the (T) contribution.
* core correlation contributions are obtained at the CCSD(T)/MTsmall level;
* scalar relativistic and, where necessary, spin-orbit coupling effects are treated at the ACPF/MTsmall level. As in W2 theory, this latter step can be combined in a single job with the previous step.

W1 theory can be applied to fairly large systems (see below). CPU times are dominated by the inner-shell correlation contribution (particularly for second-row compounds), which is reflected in the relatively small time reduction compared to W2 theory — e.g., from 1h12 to 24 for CO and from 13h42 to 8h48 for OCS. In addition — contrary to W2 theory — W1 theory exhibits a pronounced difference in performance between first-row and second-row compounds: for the species in Table VI, MAE is 0.26 kcal/mol for first-row, but 0.40 kcal/mol for second-row compounds. Since the CPU time gap between W1 and W2 theory is fairly narrow for second-row species, we conclude that for accurate work on second-row species — unless precluded by disk space or memory limitations — it may well be worth to ‘walk the extra mile’ and carry out a W2 rather than a W1 calculation. For first-row systems, on the contrary, W1 may well seem the more attractive of the two.

## VII Sample applications to larger systems

By way of illustration, we have carried out some W1 theory calculations on trans-1,3-butadiene and benzene. All relevant computed and observed results are summarized in Table IX. The example of benzene is representative and will be discussed here in detail — it should be mentioned that the calculation was carried out in its entirety on an SGI Octane workstation with 2x18 GB external SCSI disks.

The reference geometry was obtained at the B3LYP/cc-pVTZ level. The zero-point energy at that level, after scaling by 0.985, is found to be 62.04 kcal/mol.

The SCF component of TAE is predicted to be 1044.95 kcal/mol at the one-particle basis set limit, of which only 0.39 kcal/mol is covered by the geometric extrapolation. Of the CCSD valence correlation component of 291.07 kcal/mol, however, some 10.11 kcal/mol is covered by the extrapolation, which also accounts for 2.13 kcal/mol out of the 26.55 kcal/mol connected triple excitations contribution.

The inner-shell correlation contribution is quite sizable at 7.09 kcal/mol, although this number is not qualitatively different from that for three acetylenes or three ethylenes. Finally, Darwin and mass-velocity terms contribute a small but significant -0.96 kcal/mol, and atomic spin-orbit splitting another -0.51 kcal/mol. All adds up to 1367.95 kcal/mol at the bottom of the well, or 1305.92 kcal/mol at 0 K, which is in excellent agreement with the experimental value of 1306.1$`\pm `$0.12 kcal/mol from the NIST WebBook.

The CCSD/VQZ calculation took 10h10’ with MOLPRO on the Octane, the CCSD(T)/VTZ calculation 1h48’ on a single CPU on the Origin 2000. By far the most time-consuming part of the calculation was the inner-shell correlation contribution, at 67h46’, to which another 4h52’ should be added for the Darwin and mass-velocity contribution. We see similar trends in the results for trans-butadiene, which agree with experiment to virtually within the stated experimental uncertainty; for allene, we obtain a value intermediate between the two experimental values proposed in the WebBook.

We find for both molecules that the sum of core-correlation and relativistic contributions can be quite well estimated by additivity approximations. For instance, the core correlation and scalar relativistic contributions with the same basis set for C<sub>2</sub>H<sub>4</sub> are +2.360 and -0.330 kcal/mol, respectively, adding up to 2.030 kcal/mol. Assuming 2 and 3 times this ‘C=C bond equivalent’ for butadiene and benzene, respectively, yields estimated contributions of 4.06 (trans-butadiene), and 6.09 (benzene) kcal/mol, which agree excellently with the directly computed values of 4.02 and 6.12 kcal/mol, respectively. Considering that inner-shell correlation effects should be fairly local in character, such schemes should work quite well for larger organic systems where the valence calculation would still be feasible but the explicit inner-shell calculation would not be.

## VIII Conclusions

We have developed and presented two quasi-‘black box’ schemes for high-accuracy calculation of molecular atomization energies or, equivalently, molecular heats of formation, of first-and second-row compounds.

The less expensive scheme, W1 (Weizmann-1) theory, yields a mean absolute error of 0.30 kcal/mol and includes only a single, molecule-independent, empirical parameter. It requires no larger-scale calculations than CCSD/AVQZ+2d1f and CCSD(T)/AVTZ+2d1f (or, for nonpolar first-row compounds, CCSD/VQZ and CCSD(T)/VTZ). On workstation computers and using conventional coupled cluster algorithms, systems as large as benzene can be treated, while larger systems are feasible using direct coupled cluster methods.

The more expensive scheme, W2 (Weizmann-2) theory, contains no empirical parameters at all and yields a mean absolute error of 0.23 kcal/mol, which is lowered to 0.18 kcal/mol for molecules dominated by dynamical correlation. On workstation computers, molecules with up to three heavy atoms can be treated using conventional coupled cluster algorithms, while larger systems can still be treated using a direct CCSD code.

The inclusion of scalar relativistic (Darwin and mass-velocity) corrections is essential for good results in second-row compounds, particularly highly polar ones. Inclusion of inner-shell correlation contributions is absolutely essential: the basis set denoted as MTsmall (for Martin-Taylor small) appears to represent the best compromise between quality and computational expense. We do not recommend the use of lower-level electron correlation methods than CCSD(T) for the evaluation of the inner-shell contribution.

Among the several infinite-basis set extrapolation formulas for the correlation energy examined, the three-parameter $`A+B/(l+1/2)^\alpha `$ expression proposed by Martin and the $`A+B/l^3`$ expression proposed by Helgaker and coworkers yield the best results for sufficiently large basis sets, with the latter formula to be preferred on grounds of stability of the extrapolated results with the basis sets used. Geometric and mixed geometric-Gaussian extrapolation formulas are unsatisfactory when applied to the correlation energy, although they appear to be appropriate for the SCF component.

The main limiting factor for the quality of our calculations at this stage appears to be imperfections in the CCSD(T) method. This assertion is supported by the fact that the mean absolute error in the computed electron affinities of the atoms H, B–F and Al–Cl drops from 0.009 eV to 0.0009 eV if CCSDT and full CI corrections are included.

Extrapolation of the (T) contribution to the correlation energy can, at no loss in accuracy, be carried out using smaller basis sets than the CCSD contribution.

###### Acknowledgements.

JM is a Yigal Allon Fellow, the incumbent of the Helen and Milton A. Kimmelman Career Development Chair (Weizmann Institute), and an Honorary Research Associate (“Onderzoeksleider in eremandaat”) of the National Science Foundation of Belgium (NFWO/FNRS). GdO acknowledges the Feinberg Graduate School (Weizmann Institute) for a Postdoctoral Fellowship. This research was supported by the Minerva Foundation, Munich, Germany. The authors acknowledge enlightening discussions with (in alphabetical order) Drs. C. W. Bauschlicher Jr. (NASA Ames Research Center), Dr. Thom H. Dunning (PNNL), Prof. Trygve Helgaker (U. of Oslo, Norway), Dr. Frank Jensen (Odense U., Denmark), Dr. Timothy J. Lee (NASA Ames Research Center), and Prof. Peter R. Taylor (UCSD and National Partnership for Advanced Computing Infrastructure), and thank Dr. Peter Stern for technical assistance with various computer systems.