High-level thermochemistry for the octasulfur ring: A converged coupled cluster perspective for a challenging second-row system

Sulfur clusters are challenging targets for high-level ab initio procedures. The heat of formation of the most common and energetically stable S8 allotrope (α-sulfur) has not been the subject of a high-level ab initio investigation. We apply the Weizmann-n computational thermochemistry protocols to the S8 sulfur cluster. We show that calculating the heat of formation with sub-chemical accuracy requires accurate treatment of post-CCSD (T), core-valence, scalar relativistic, and zero-point vibrational energy contributions. At the relativistic, allelectron CCSDT(Q)/CBS level of theory we obtain an enthalpy of formation at 0 K of ∆fH0 = 24.44 kcal mol–1, and at 298 K of ∆fH298 = 23.51 kcal mol–1. These values suggest that the experimental values from Gurvich (∆fH0 = 25.1 ± 0.5 kcal mol–1) and JANAF (∆fH0 = 24.95 ± 0.15 and ∆fH298 = 24.00 ± 0.15 kcal mol–1) represent overestimations and should be revised downward by 0.5–0.7 kcal mol–1. We also show that computationally economical composite ab initio protocols such as G4, G4(MP2), and CBS-QB3 are unable to achieve chemical accuracy relative to our best CCSDT(Q)/CBS heat of formation for S8.


Introduction
Sulfur has a large number of known allotropes and small sulfur clusters play a key role in the chemistry of next-generation lithium-sulfur batteries. [1,2,3,4,5,6,7] The chemistry of sulfur clusters has been the subject of many experimental [8,9,10,11,12,13,14] and theoretical [15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30] studies. The most energetically stable allotrope of sulfur is orthorhombic α-sulfur (hereinafter referred to as simply S 8 ) and perhaps the most fundamental thermochemical property of S 8 is its heat of formation. The experimental heat of formation of S 8 has been obtained by Guthrie et al. in 1954, [31] they arrived at ∆ f H 0 = 25.23 ± 0.05 kcal mol -1 and ∆ f H 298 = 24.35 kcal mol -1 based on the assumption that the saturated sulfur vapour in the experiments consisted of only S 6 and S 8 structures. In 1989, Gurvich et al. revised the value at 0 K slightly downwards to ∆ f H 0 = 25.1 ± 0.5 kcal mol -1 based on additional data available in the literature. [32] Despite the importance of the S 8 cluster, its heat of formation has not been the subject of a high-level theoretical investigation. With the advancement of computer hardware and theoretical methods over the past decade, it is now possible to calculate gas-phase thermochemical data with sub-kcal mol -1 (or even sub-kJ mol -1 ) accuracy using highly accurate composite ab initio methods (e.g., Gn, [33] ccCA, [34] HEAT, [35,36,37] Wn, [38,39,40] Wn-F12, [41,42] and others). [43,44,45] In particular, composite methods that approximate the FCI/CBS energy (i. e., full configuration interaction at the complete basis set limit) are capable of predicting atomization energies (or their cognates, molecular heats of formation) with sub-kJ mol -1 accuracy. [38,39,46,47] This level of accuracy is often better than that associated with many experimental thermochemical determinations, including many of the experimental values in the NIST Chemistry WebBook (i.e., most experimental values are associated with error bars larger than 1 kJ mol -1 , and often larger than 1 kcal mol -1 ). [48] Thus, highly accurate composite theories are instrumental in obtaining new thermochemical data with benchmark accuracy and in critical assessment of published thermochemical data.
In the present work we calculate the total atomization energy (TAE) of the S 8 cluster by means of W4lite theory. [38] This theory includes both secondary energetic contributions (e.g., core-valence, scalar relativistic, atomic spin-orbit, zero-point vibrational energy, and deviations from the Born-Oppenheimer approximation) and post-CCSD(T) contributions up to quasiperturbative quadruple excitations. We show that the T-(T), (Q), core-valence, scalar relativistic, atomic spin-orbit, and zero-point vibrational energy contributions to the TAE can exceed 1 kcal mol -1 and are needed for determining the TAE with sub-chemical accuracy.

Computational details
In order to obtain an accurate heat of formation for the orthorhombic S 8 ring in D 4d symmetry we carried high-level, ab initio calculations with the W2.2 (denoted by W2) and W4lite (denoted by W4L) thermochemical protocols. [38,40] Hereinafter, the cc-pV(n+d)Z and aug-cc-pV (n+d)Z correlation-consistent basis sets are denoted by VnZ and AVnZ, respectively, and the notation V{X,Y}Z indicates extrapolation from the VXZ and VYZ basis sets. [49,50,51] In the W4L procedure, the Hartree-Fock (HF) energy is extrapolated from the AV{5,6}Z basis sets, using the Karton-Martin HF extrapolation formula. [52] The valence CCSD correlation contribution is extrapolated from the same basis sets, where the singlet-and triplet-pair energies are extrapolated separately to the infinite basis set limit, using the E(L) = E ∞ + A/L α two-point extrapolation formula, with α = 3 and 5 for the singlet-and triplet-pair energies, respectively. [38,39,40] The (T) valence correlation component (i.e., CCSD(T) -CCSD) is extrapolated from the AVQZ and AV5Z basis sets using the above two-point extrapolation formula with α = 3. The CCSD(T) inner-shell contribution is extrapolated from the core-valence weighted correlation-consistent aug-cc-pwCV{T,Q}Z (APWCVnZ) basis sets of Peterson and Dunning. [53] The scalar relativistic contribution is obtained from second-order Douglas-Kroll-Hess [54,55] CCSD(T)/AVQZ-DK calculations. [56] The diagonal Born-Oppenheimer corrections (DBOCs) are calculated at the HF/AVTZ level of theory.
The CCSDT energy is calculated in conjunction with two truncated versions of the VTZ basis set. Namely, the spd part of the VTZ basis (denoted by VTZ(nof2d)) and the sp part of the VTZ basis set combined with the d function from the VDZ basis set denoted by VTZ(nof1d)). [40] These two basis sets are used in conjunction with the VDZ basis set to obtain the CCSDT-CCSD(T) component (denoted by T-(T)) using the two-point extrapolation formulas developed in [ref . 40], which do not depend on the highest angular momentum present in the basis set. The CCSDT(Q) energy is calculated in conjunction with the VDZ basis set and the (Q)/VDZ component is scaled by 1.1 as recommended in [refs. 38] and [40].
The S 8 geometry was optimized at the B3LYP, [57,58,59] B3LYP-D3BJ, [60] and CCSD(T) levels in conjunction with several valence and core-valence weighted correlation-consistent basis sets. All the large-scale CCSD(T) calculations were performed using the Molpro program suite, [61,62] the post-CSD(T) calculations were performed with the MRCC program, [63,64] the DBOC calculations were carried out with the CFOUR program, [65] and DFT geometry optimizations and economical composite ab initio calculations were performed with the Gaussian 16 program. [66]

Results and discussion
The principal conformer of the cyclic S 8 molecule is the crownshaped structure in D 4d symmetry. Before proceeding to a detailed discussion of the CCSDT(Q)/CBS results from W4L theory, it is instructive to examine the level of theory used for optimizing the reference geometry. Whilst the higher-level Wn protocols (e.g., W4L, W4, and above) [38,39] use a CCSD(T)/VQZ reference geometry, lower-level Wn theories (e.g., W1 and W2) prescribe the use of B3LYP/VTZ geometries. [40,41] In this context it is important to note that it was found that the CCSD(T)/VQZ level of theory reproduces CCSD(T)/AV{5,6}Z reference bond distances with a root-mean-square deviation of 0.0015 Å for a large set of 108 organic and inorganic molecules. [67] The bond distance and angle obtained at with the B3LYP and CCSD(T) levels are depicted in Table 1. Both the B3LYP/VTZ and B3LYP-D3BJ/VTZ levels of theory overestimate our best CCSD(T)/VQZ bond length, namely, by 0.017 and 0.012 Å, respectively. These results are consistent with recent observations that the B3LYP/VTZ level of theory tends to overestimate S-S bond distances relative to CCSD(T)/AV{5,6}Z reference values in small sulfur containing molecules such as S 2 , S 2 O, and S 2 H. [68] We can estimate the effect of using a B3LYP geometry rather than a CCSD(T)/VQZ geometry at the CCSD(T)/AV{T,Q}Z level of theory, i.e., by comparing the CCSD(T)/AV{T,Q}Z//B3LYP/VTZ and CCSD(T)/AV {T,Q}Z//CCSD(T)/VQZ energies. Using a B3LYP geometry rather than a CCSD(T)/VQZ geometry affects the CCSD(T)/AV{T,Q}Z energy by as much as 1.2 (B3LYP/VTZ) and 0.6 (B3LYP-D3BJ/VTZ) kcal mol -1 .
Moving to the CCSD(T) geometries, the CCSD(T)/VDZ level of theory significantly overestimates the CCSD(T)/VQZ bond length by 0.032 Å. This difference in geometry affects the CCSD(T)/AV{T,Q}Z energy by a staggering amount of 2.3 (!!) kcal mol -1 . This deviation in the S-S bond distance goes down to 0.011 Å at the CCSD(T)/VTZ level of theory. Nevertheless, the CCSD(T)/VTZ reference geometry still leads to a deviation in the CCSD(T)/AV{T,Q}Z energy by 0.5 kcal mol -1 relative to the CCSD(T)/VQZ geometry. The above results illustrate the problems associated with using DFT geometries or CCSD(T) geometries optimized with a double-ζ or even a triple-ζ basis set for a system containing many second-row atoms.
The above results raise the question of the effect of inner-shell correlation on the S 8 geometry. It has been shown that inclusion of innershell correlation in CCSD(T) geometry optimizations can affect bond distances of small second-row compounds by up to several milliangstroms. [39,53,69,70,71,72,73,74] The effect of these small geometric changes on the resulting total atomization energy (or heat of formation) has been examined in detail in [39] and [75]. It has been found that the effect of the reference geometry is largely confined to the CCSD(T) component and can affect the final total atomization energy by up to about 0.1 kcal mol -1 . Whilst we are unable to optimize the geometry of S 8 with all electrons correlated (apart from the 1 s deep-core orbitals) at the CCSD(T)/PWCVQZ level of theory, we can examine the geometry effect in conjunction with the PWCVTZ basis set. Using core-valence correlated CCSD(T)/PWCVDZ reference geometry reduces the S-S bond distance by 0.007 Å relative to the frozen-core CCSD (T)/VTZ bond distance (Table 1). This geometry change affects the CCSD (T)/AV{T,Q}Z energy by an appreciable amount of 0.34 kcal mol -1 . It is expected that the geometry effect of correlating the core-valence electrons will become smaller when moving to a quadruple-ζ quality basis set. Table 2 lists the component breakdown of the total atomization energy at the bottom of the well (TAE e ) for the D 4d symmetry S 8 molecule. The HF/AV{Q,5}Z (W2 theory) component amounts to 253.80 kcal mol -1 . However, at the HF/AV{5,6}Z level (W4L theory) it increases to 254.12 kcal mol -1 . Thus, the AV{Q,5}Z underestimates our best CBS value by 0.32 kcal mol -1 . We note that this underestimation is roughly twice as large as that observed for the S 4 cluster (in C 2v symmetry), namely, for S 4 an underestimation of 0.14 kcal mol -1 is observed. [46]  Theoretical geometries for the S 8 structure in D 4d symmetry, the bond distance (r) and angle (θ) are given in Å and degrees, respectively. Energy differences due to geometry effects (∆E geom ) are given in kcal mol -1 . The CCSD/AV{Q,5}Z correlation component (225.43 kcal mol -1 ) overestimates the CCSD/AV{5,6}Z correlation component (224.54 kcal mol -1 ) by as much as 0.89 kcal mol -1 . Again, it is interesting to note that this overestimation is about three times as large as that observed for the smaller S 4 molecule, namely for S 4 an overestimation of 0.27 kcal mol -1 is observed for the same basis set pairs. [46] The (T)/AV{T,Q}Z correlation component (38.18 kcal mol -1 ) underestimates the (T)/AV{Q,5}Z correlation component (38.60 kcal mol -1 ) by 0.42 kcal mol -1 . This underestimation is nearly twice as large as that observed for the S 4 molecule, namely for S 4 an underestimation of 0.23 kcal mol -1 is observed for the same basis set pairs. [46] The above results demonstrate that there is a fair amount of energy reshuffling between the HF, CCSD, and (T) components when moving from W2 to W4L theory. In particular, the HF/AV{Q,5}Z and (T)/AV{T, Q}Z components from W2 theory underestimate the HF/AV{5,6}Z and (T)/AV{Q,5}Z components from W4L theory by 0.32 and 0.42 kcal mol -1 , respectively (vide supra). On the other hand, the CCSD/AV{Q,5} Z correlation component from W2 theory overestimates the CCSD/AV {5,6}Z correlation component from W4L theory by as much as 0.89 kcal mol -1 . Thus, overall, despite the use of smaller basis sets in the CCSD(T) calculations in W2 theory, the difference between the CCSD(T)/CBS estimates from W2 and W4L theories amounts to just 0.15 kcal mol -1 .
Calculating the higher-order connected triples (T-(T)) component with the VDZ and VTZ(nof1d) basis sets, reduces the atomization energy by -0.84 kcal mol -1 . However, using the VDZ and VTZ(nof2d) basis sets reduces the atomization energy by a significantly larger amount of -2.99 kcal mol -1 . Thus, the VTZ(nof1d) basis set is clearly too small for obtaining useful estimation of the T-(T) component for this challenging system. Unfortunately, the CCSDT/VTZ calculation is not feasible for a system with eight sulfur atoms since it involves 32.6 × 10 9 amplitudes, compared to 12.7 × 10 9 and 5.1 × 10 9 amplitudes involved in the CCSDT/VTZ(nof2d) and CCSDT/VTZ(nof1d) calculations, respectively. To put these numbers in perspective, the CCSDT/VTZ(nof2d) calculation required 15 iterations to converge to a submicrohartree level and each iteration ran for about 30 h on a dual Intel Xeon machine with 40 cores and 1024 GB of RAM. Turning to the parenthetical connected quadruple excitations, the (Q)/VDZ contribution to the TAE e amounts to 2.89 kcal mol -1 and nearly perfectly cancels out with our best T-(T) component (-2.99 kcal mol -1 ). Overall, post-CCSD(T) contributions to the TAE e amount to merely 0.10 kcal/mol.
The core-valence correction amounts to 1.88 kcal mol -1 at the CCSD (T)/APWCV{T,Q}Z level of theory. Both the scalar relativistic and spinorbit contribution strongly destabilize the molecule and reduce the TAE e by -1.52 and -4.48 kcal mol -1 , respectively. The DBOC term is practically negligible and the so-called M-A correction amounts to less than 0.05 kcal mol -1 (Table 2) The above results demonstrate that, due to error cancelation between various components, the final heat of formation from W2 theory is in good agreement with that obtained from the computationally more intensive W4L theory. In particular, the difference between the two values amounts to just 0.22 kcal mol -1 . It is therefore of interest to examine the performance of more approximate composite ab initio procedures for the heat of formation of S 8 . In this context it should be pointed out that recent high-level ab initio studies found that the calculation of the total atomization energies of small sulfur clusters (e.g., S 3 and S 4 ) is a major challenge for high-level composite ab initio methods. [82,83] We consider the following composite ab initio procedures: G4, [84] G4(MP2), [85] G4(MP2)-6X, [86] G3, [87] G3(MP2), [88] G3B3, [89] G3(MP2)B3, [89] and CBS-QB3. [90,91] To enable a more meaningful analysis, apart from S 8 , we also include S 2 , S 3 , and S 4 in our error statistics. The TAE 0 reference values are taken from W4L theory for S 8 and from W4 theory for the smaller molecules. [46] Table 3 lists the individual deviations and mean absolute deviations (MADs) for the approximate composite ab initio procedures. Inspections of these results reveals that, consistent with results for other systems containing many second-row atoms, [92] the Gn(MP2)-type procedures outperform the Gn-type procedures. In particular, both G3(MP2) and G3(MP2)B3 attain MADs of 0.9 and 1.0 kcal mol -1 , respectively, and G4(MP2)-6X theory clocks in at a slightly higher MAD of 1.2 kcal mol -1 . It is noteworthy that the largest deviations for these three methods do not exceed 2.4 kcal mol -1 , and therefore they are recommended for the calculation of enthalpies of formation of larger sulfur compounds. The G4 method shows relatively poor performance (MAD = 1.9 kcal mol -1 ) at a higher computational cost and is therefore not recommended. The other composite ab initio procedures considered give MADs ranging between 2.3 (CBS-QB3) and 4.9 (G3) kcal mol -1 . It should be pointed out that MADs reported in Table 3 translate into much larger 95% confidence intervals for the Gn and CBS procedures. In particular, the conversion factors between MADs and 95% confidence intervals range between 2.5-3.5 depending on sample size and distribution [93]
Supplementary data: Optimized geometries for all the levels of theory considered in Table 1 (Table S1).

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Table 3
Deviations and mean absolute deviations (MADs) from W4 and W4L total atomization energies (TAE 0 ) values for popular Gn and CBS composite ab initio procedures (kcal mol -1 ). G3