# The heat of atomization of sulfur trioxide, SO3 — a benchmark for computational thermochemistry

## I Introduction

Neither the sulfuric anhydride (SO<sub>3</sub>) molecule, nor its importance in atmospheric and industrial chemistry, require any introduction to the chemist.

SO<sub>3</sub> displays somewhat unusual bonding. While it is often cited as a ‘hypervalent molecule’ in undergraduate inorganic chemistry textbooks, quantitative theories of chemical bonding such as atoms-in-molecules unequivocally show (see Ref. for a lucid review and discussion) that there are no grounds for invoking violation of the octet rule in SO<sub>3</sub> (or, for that matter, most second-row molecules), and that bonding in SO<sub>3</sub> is best seen as a combination of moderately polar $`\sigma `$ bonds with highly polar $`p_{\pi ,S},p_{\pi ,O}`$ bonds.

Previous experience on BF<sub>3</sub> and SiF<sub>4</sub> suggests that in molecules with several strong and very polar bonds, basis set convergence will be particularly slow. In addition, in a recent calibration study on the anharmonic force field of SO<sub>3</sub> it was found that the molecule represented a fairly extreme example of a phenomenon noted previously for second-row molecules — namely the great sensitivity of the SCF part of computed properties to the presence of so-called ‘inner polarization functions’, i.e. high-exponent $`d`$ and $`f`$ functions.

Very recently, Martin and de Oliveira published a standard protocol known as W2 (Weizmann-2) theory that was able to predict total atomization energies of a fairly wide variety of molecules (including SO<sub>2</sub>, which is relevant for this work) to better than 0.23 kcal/mol on average (0.18 kcal/mol for molecules dominated by a single reference configuration). Application of this method to SO<sub>3</sub> requires a CCSD (coupled cluster with all single and double excitations) calculation with 529 basis functions in the $`C_{2v}`$ nondegenerate subgroup, which was well beyond our available computational resources, particularly in terms of disk space.

Very recently, however, Schütz et al. developed a general implementation of integral-direct correlated methods that made possible, inter alia, CCSD calculations on basis sets this size on workstation computers. Consequently, we carried out a benchmark calculation on the heat of atomization of SO<sub>3</sub>, which is reported in the present work.

Having obtained the benchmark ab initio value, we will assess the performance of some less computationally demanding schemes. This includes W1 theory, which is much more cost-effective than W2 theory but performs much less well for second-row than for first-row compounds. From an analysis of the SO<sub>3</sub> results, we will derive a minor modification (denoted W1 theory) which in effect largely removes this disadvantage.

## II Methods

Most electronic structure calculations were carried out using MOLPRO98.1 (with integral-direct code installed) running on a DEC Alpha 500/500 workstation at the Weizmann Institute of Science. Some additional calculations were carried out using GAUSSIAN 98 running on the same platform.

As in our previous work on SO<sub>2</sub>, the CCSD(T) electron correlation method, as implemented by Hampel et al., has been used throughout. The acronym stands for coupled cluster with all single and double substitutions augmented by a quasiperturbative account for triple excitations. From extensive studies (see  for a review) this method is known to yield correlation energies very close to the exact $`n`$-particle solution within the given basis set as long as the Hartree-Fock determinant is a reasonably good zero-order reference wave function. None of the usual indicators ($`𝒯_1`$ diagnostic, largest excitation amplitudes, or natural orbital occupancies of first few HOMOs and LUMOs) suggest a significant departure from the single-reference regime. (For the record, $`𝒯_1`$=0.018 for SO<sub>3</sub>.)

Valence correlation basis sets are built upon the augmented correlation-consistent polarized $`n`$-tuple zeta (aug-cc-pV$`n`$Z, or AV$`n`$Z for short) basis sets of Dunning and coworkers. In this work, we have considered AVDZ, AVTZ, AVQZ, and AV5Z basis sets, with maximum angular momenta $`l`$=2 ($`d`$), 3 ($`f`$), 4 ($`g`$), and 5 ($`h`$), respectively. The effect of inner polarization was accounted for by adding ‘tight’ (high-exponent) $`d`$ and $`f`$ functions with exponents that follow even-tempered series $`\alpha \beta ^n`$, with $`\alpha `$ the tightest exponent of that angular momentum in the underlying basis set and $`\beta `$=2.5. Such basis sets are denoted AV$`n`$Z+d, AV$`n`$Z+2d, and AV$`n`$Z+2d1f. The largest basis set considered in the present work, AV5Z+2d1f, corresponds to $`[8s7p7d5f3g2h]`$ on sulfur and $`[7s6p5d4f3g2h]`$ on oxygen (148 and 127 contracted basis functions, respectively), adding up to 529 basis functions for the entire molecule. The CCSD calculation in this basis set was carried out using the newly implemented direct algorithm; all other CCSD and CCSD(T) calculations were done conventionally.

The effect of inner-shell correlation was considered at the CCSD(T) level using two specialized core correlation basis sets, namely the Martin-Taylor (MT) basis set used in previous work on SO<sub>2</sub>, and the somewhat more compact MTsmall basis set that is used in the W2 protocol for this purpose. Correlation from the sulfur ($`1s`$) orbital was not considered, since this lies too deep to meaningfully interact with the valence orbitals. Scalar relativistic effects were computed as expectation values of the first-order Darwin and mass-velocity corrections for the ACPF (averaged coupled pair functional) wave function with the abovementioned core correlation basis sets. (All electrons were correlated in these calculations since relativistic effects are most important for the electrons closest to the nucleus.)

The CCSD(T)/VQZ+1 reference geometry used throughout this work, $`r_{SO}`$=1.42279 Å, was taken from the earlier spectroscopic work on SO<sub>3</sub>, as was the anharmonic zero-point energy of 7.794 kcal/mol.

## III Results and discussion

The most striking feature of the basis set convergence at the SCF level (Table 1) is certainly the great importance of inner polarization functions: augmenting the AVDZ basis set with two tight functions on S has an effect of no less than 40.5 kcal/mol! The same operation affects the AVTZ SCF binding energy by 15.7 kcal/mol, and even from AVQZ to AVQZ+2d the effect is still 8.6 kcal/mol, probably the largest such effect hitherto observed. In addition augmenting the basis set by a tight $`f`$ function has an effect of 1.1 kcal/mol from AVTZ+2d to AVTZ+2d1f, but only 0.16 kcal/mol from AVQZ+2d to AVQZ+2d1f. Presumably the effect from AV5Z+2d to AV5Z+2d1f will be next to negligible.

Not surprisingly, this translates into a substantial effect on the extrapolated SCF limit. A geometric extrapolation from the AV{D,T,Q}Z results would yield 153.64 kcal/mol as the SCF limit, 6.3 kcal/mol less than the AV{T,Q,5}Z+2d1f limit employed in W2 theory. The AV{D,T,Q}Z+2d limit, on the other hand, if fairly close to the latter at 159.7 kcal/mol. (Our best SCF limit is 159.90 kcal/mol, of which the extrapolation accounts for 0.15 kcal/mol.)

This type of variability is almost completely absent for the correlation energy, where AV$`n`$Z, AV$`n`$Z+2d and AV$`n`$Z+2d1f largely yield the same answers. Following the W2 protocol, the CCSD correlation energy is extrapolated using the $`A+B/l^3`$ extrapolation formula of Halkier et al. to CCSD/AV{Q,5}Z+2d1f energies (for which $`l`$={4,5}). (For a fairly comprehensive review of theoretical and empirical arguments in favor of this type of extrapolation, see Ref. and references therein.) We thus obtain 165.94 kcal/mol as our best estimate for the CCSD correlation contribution to TAE. It should be noted that the extrapolation accounts for 3.2 kcal/mol of this amount: basis set convergence is indeed quite slow. We note that the largest direct CCSD calculation took a solid two weeks of CPU time on the DEC Alpha — a conventional calculation would have required about 60 GB of temporary disk space, as well as a much higher I/O bandwidth if a reasonable wall time to CPU time ratio were to be attained.

As a general rule, the (T) contribution converges much more rapidly with basis set (besides being smaller to begin with) and therefore, we were able to dispense entirely with the CCSD(T)/AV5Z+2d1f calculation. From CCSD(T)/AV{T,Q}+2d1f results and the $`A+B/l^3`$ formula, we obtain a basis set limit for the (T) contribution of 20.17 kcal/mol, in which the extrapolation accounts for 0.57 kcal/mol. Together with the CCSD results, this adds up to a valence correlation contribution to TAE\[SO<sub>3</sub>\] of 186.11 kcal/mol, of which 3.75 kcal/mol is covered by extrapolations.

The inner-shell correlation contribution (Table 2) at the CCSD(T) level using the Martin-Taylor core-correlation basis set, was found to be 0.89 kcal/mol with the Martin-Taylor core correlation basis set, and 0.96 kcal/mol with the somewhat more compact MTsmall basis set used in W2 theory. Bauschlicher and Ricca found that basis set superposition error significantly affects the inner-shell correlation contribution in SO<sub>2</sub>. It was evaluated here using the site-site counterpoise method ; we thus found counterpoise-corrected core correlation contributions of 0.73 kcal/mol with the Martin-Taylor and 0.68 kcal/mol with the MTsmall basis sets.

Scalar relativistic effects were obtained as expectation values of the mass-velocity and Darwin operators for the ACPF (averaged coupled pair functional) wavefunction. Their effect on the computed TAE (with either core correlation basis set) is -1.71 kcal/mol, comparable to the -1.88 kcal/mol previously found for SiF<sub>4</sub>. Atomic spin-orbit splitting adds another -1.23 kcal/mol to the result. (These latter two terms together imply a relativistic contribution of -2.94 kcal/mol, or nearly 1% of the atomization energy.)

Finally, we obtain a W2 total atomization energy at the bottom of the well, TAE<sub>e</sub>, of 344.03 kcal/mol; using the BSSE-corrected inner shell correlation contribution, this value drops to 343.76 kcal/mol. In combination with the very accurate ZPVE=7.795 kcal/mol, we finally obtain, at absolute zero, TAE<sub>0</sub>=336.17 kcal/mol without, and 335.96 kcal/mol with, BSSE correction on the core correlation contribution. This latter value is in perfect agreement with the experimental TAE<sub>0</sub>=335.92$`\pm `$0.19 listed in the Gurvich compilation. We thus see once more the importance of including BSSE corrections for the inner-shell correlation part of TAE: it should be noted that while the inner-shell contribution to TAE is small, the S($`2s,2p`$);O($`1s`$) absolute correlation energy is comparable with the valence correlation energy in SO<sub>3</sub>. BSSE on the valence contribution is much less of an issue since the basis sets used for valence correlation are much more saturated to begin with, and furthermore the valence correlation energy is being extrapolated to the infinite-basis limit where it should vanish by definition.

The performance of more approximate computational thermochemistry schemes is of some interest here (Table 3). G1 theory is in error by no less than -11.4 kcal/mol, which goes down to -6.9 kcal/mol for G2 theory and -5.45 kcal/mol for G3 theory. (Only the latter includes spin-orbit splitting as part of the protocol: none of these methods consider scalar relativistic effects.) G2(MP2) performs relatively well as a result of error compensation (-2.4 kcal/mol). The CBS-Q scheme underestimates the true binding energy by only 1 kcal/mol, while CBS-QB3 is only 0.2 kcal/mol above experiment. It should be noted that neither CBS-Q nor CBS-QB3 include relativistic effects of any kind as part of the standard protocol; therefore the excellent performance put in by these methods is to a large extent thanks to error compensation. Finally, the W1 theory of Martin and de Oliveira — which yields a mean absolute error of about 0.3 kcal/mol for a wide variety of compounds — has an error in TAE<sub>0</sub>\[SO<sub>3</sub>\] of -1.13 kcal/mol. (W1 theory includes both scalar relativistic and spin-orbit contributions.)

The largest calculations involved in the W1 protocol are CCSD/AVQZ+2d1f and CCSD(T)/AVTZ+2d1f, which is still rather more demanding than the steps in any of the G$`n`$ or CBS methods. Hence this performance is rather disappointing — a failure of W1 theory was also noted for SO<sub>2</sub> in the original paper. Balance considerations may lead us to wonder whether an AVTZ+2d1f basis set is not rather top-heavy on inner polarization functions. Using the AV$`n`$Z+2d series favored by Bauschlicher and coworkers (e.g.) indeed reduces the discrepancy with experiment by 0.55 kcal/mol (of which 0.20 kcal/mol in the SCF part). The alternative sequence {AVDZ+2d,AVTZ+2d,AVQZ+2d1f} yields even better agreement with experiment (and the more rigorous calculations): in fact, the final value thus obtained falls within the experimental error bar. Particularly encouraging is the fact that the predicted SCF limit is now within 0.04 kcal/mol of our best estimate. Preliminary calculations on other second-row systems suggest that this procedure, which we will label W1 theory, may be preferable over standard W1 theory for second-row systems with strong inner shell polarization. (The two variants are equivalent for first-row compounds.)

As a test, we have taken three molecules for which W1 yields fairly large errors (CS, SO, and SO<sub>2</sub>) and repeated the calculation using W1 theory. Deviations from experiment drop from -0.92, -0.62, and -1.01 kcal/mol, respectively, to -0.56, -0.32, and -0.02 kcal/mol, respectively, which is not qualitatively different from the vastly more expensive W2 calculations which yielded deviations of -0.51, +0.02, and +0.23 kcal/mol for these molecules. We conclude that W1 theory indeed represents an improvement, and recommend it for future work on second-row systems instead of W1 theory.

## IV Conclusions

Benchmark ab initio calculations using direct coupled cluster methods predict the total atomization energy at 0 K of SO<sub>3</sub> to be 335.96 (observed 335.92$`\pm `$0.19) kcal/mol. The computed results includes extrapolation to the basis set limit (3.75 kcal/mol), relativistic effects (-2.94 kcal/mol), inner-shell correlation (0.68 kcal/mol after BSSE correction), and anharmonic zero-point energy (7.755 kcal/mol). Inner polarization functions make very large (40 kcal/mol with $`spd`$, 10 kcal/mol with $`spdfg`$ basis sets) contributions to the SCF part of the binding energy. The molecule presents an unusual hurdle for less computationally intensive theoretical thermochemistry methods and is proposed as a benchmark for them. A slight modification of W1 theory is proposed which appears to result in improved performance for second-row systems with strong inner-shell polarization effects.

###### Acknowledgements.

JM is a Yigal Allon Fellow, an Honorary Research Associate (“Onderzoeksleider in eremandaat”) of the National Science Foundation of Belgium (NFWO/FNRS), and the incumbent of the Helen and Milton A. Kimmelman Career Development Chair. He thanks Prof. Peter J. Knowles (Birmingham University, UK) for assistance with the installation of the direct coupled cluster code, and Dr. Charles W. Bauschlicher Jr. (NASA Ames Research Center, Moffett Field, CA) for critical reading of the manuscript prior to submission. This research was supported by the Minerva Foundation, Munich, Germany.