The curious case of DMSO: A CCSD(T)/CBS(aQ56+d) benchmark and DFT study

Dimethyl sulfoxide [DMSO; (CH₃)₂SO] is a common polar aprotic solvent and adopts a pyramidal geometry (left in Fig. 1), whereas methyl methanesulfenate (CH₃SOCH₃), a structural isomer of DMSO, adopts a bent geometry (right in Fig. 1). Prior work has examined DMSO and its structural isomer, as well as the cases where X, Y = O or S (where X ≠ Y) and R, R′ = H or CH₃. The corresponding relative energetics from previous theoretical investigations show a wide range of relative energies and sometimes even qualitative discrepancies (select data are given in Table I). All the previously reported relative energies can be found in Tables S69–S74 of the supplementary material.

DMSO (14), relative to the bent methyl methanesulfenate (16) isomer, was determined to be higher in energy by as much as 16.6 kcal mol⁻¹ or lower in energy by as much as 8.7 kcal mol⁻¹ with post-SCF methods and larger basis sets.^1–5 A significant basis set dependence was reported for these systems where smaller basis sets were found to be unreliable as they provided qualitatively incorrect relative energies.^3,4 For example, the reported second-order Møller–Plesset perturbation theory (MP2) relative energies (ΔE) for 14 (relative to 16) were 4.9 kcal mol⁻¹ with the aug-cc-pVDZ basis set and −6.6 kcal mol⁻¹ with the larger aug-cc-pVQZ basis set. CCSD(T)/6-311G(d) energetics also indicate DMSO (14) to be higher in energy than 16. According to MP2 and B3LYP complete basis set (CBS) limit extrapolations, DMSO (14) is lower in energy by 7.9 kcal mol⁻¹ and 3.5 kcal mol⁻¹, respectively, indicating a non-trivial method dependence on the quantitative result.⁴ For the analogous systems of HSOH (9), H₂SO (7), and H₂OS (8), one study implements basis sets as large as cc-pV(6+d)Z with the coupled-cluster with single, double, and perturbative triple excitations [CCSD(T)] method but did not report the relative electronic energies of these isomers.⁶ However, the reported relative enthalpies (ΔH) demonstrate that the pyramidal isomers (7 and 8) are higher in energy than the bent conformer (9) and a diverse set of computational studies support this conclusion.^1,6–12 Interestingly, an investigation of the chlorinated variants of 7 and 9 (i.e., Cl₂SO and ClSOCl, respectively) showed the pyramidal isomer to be lower in energy by 47.0 kcal mol⁻¹ according to CCSD(T)-F12/aug-cc-pV(T+d)Z computations.¹³ This demonstrates a significant effect of the substituent, R and R′, on the qualitative energy ordering of these structural isomers.

The computed relative energetics of RR′SO and RSOR′ (where R = CH₃ and R′ = H) indicate the pyramidal conformer (10) to be 21.5 kcal mol⁻¹ higher in energy than the bent conformer (13) according to the MP2/6-311+G** level of theory.¹⁴ However, when R = H and R′ = CH₃ (giving 10 and 12), a significant basis set dependence on the relative energies is observed.¹⁵ The range of reported relative energies of 10 (relative to 12) using double-ζ to quadruple-ζ basis sets spans from +17.7 to −0.6 kcal mol⁻¹ for this system (Table I).

The similar pyramidal and bent isomers of disulfanes [RS(S)R and RSSR, respectively] have also been explored with computational approaches.^9,12,16–20 The pyramidal structure is conclusively higher in energy. When R = H, the pyramidal form (17) is 27–36 kcal mol⁻¹ higher in energy than the bent form (18), whereas when R = CH₃, the pyramidal form (21) is about 20 kcal mol⁻¹ higher in energy than the bent form (22). There is also recent interest in sulfenic acids (RSOH) and their sulfenyl and sulfinyl tautomers as well as persulfide adducts (RSSH) due to their role in biological systems.^22–23 Recent computational studies examining thioperoxides, disulfides, and thiosulfonates in biologically relevant reactions^24,25 have modeled a few of these simple molecules but did not report the relative energetics. Others have examined the effect of varying the substituents (R = H, Li, Na, F, Cl, Br, and I) on peroxides and disulfides in both the pyramidal and bent forms and determined the bent conformer to always be lower in energy except in the case of the molecular formula FSSF.²⁶ 

Interest in these types of systems has naturally extended beyond the isolated molecule state and is moving toward characterizing the bulk phase. Recent work has characterized clusters of DMSO (14),^27,28 HOSH (9),^19,29 CH₃OSH (12),²⁹ and HSSH (18).^19,30 The interactions of DMSO with non-polar liquids²⁷ and ionic liquids^31–34 have also been examined.

While an abundance of prior work exists on these types of systems, none have uniformly characterized the relative energetics of these structural isomers with state-of-the-art methodologies in one consistent study. For example, explicitly correlated F12 methods are designed to accelerate the convergence to the CBS limit.^35–38 In addition, flexible basis sets, particularly the correlation consistent basis sets that include tight-d functions for second row main group elements, have been shown to dramatically improve computational results including atomization energies, binding energies, and dissociation energies.^39,40

The aim of this work is to systematically characterize the relative energetics of a series of small, bent molecules of the forms ROOR′, RSSR′, and ROSR′ (where R, R′ = H, CH₃) relative to their pyramidal conformers at the CCSD(T) CBS limit in conjunction with explicitly correlated F12 methods and a wide range of density functional theory (DFT) methods. These structures are given in Fig. 2 where the substituents (R, R′) are organized by columns and the chalcogen–chalcogen bonding framework (referred to here as the core with the following considered: peroxides, OO-core; thioperoxides, OS-core; and disulfides, SS-core) is organized by rows. The peroxides are included in this dataset in addition to the disulfides to serve as a bracket to the thioperoxides, which are difficult to properly characterize as clearly demonstrated by prior work. Reliable MP2/CBS and CCSD(T)/CBS energetics will be obtained given the basis set dependence identified in a couple of these systems. The pathological behavior of OS-core structures is examined via a systematic mixing of correlation consistent basis sets. Furthermore, an analysis of the performance of 22 density functional approximations used in conjunction with 12 different basis sets will be presented. Finally, the performance of various composite methodologies (ccCA, ccCA-TM, and G3B3) is evaluated. From this work, one will be able to see how basis set dependence, higher order correlation effects, and different substituents affect the relative energetics of these systems.

Full geometry optimizations and corresponding harmonic vibrational frequency computations were performed on each of the compounds in Fig. 2 in the reported point groups using the following 22 density functional theory (DFT) methods as implemented in Gaussian 09:⁴¹ APF-D,⁴² B3LYP,^43,44 B3P86,^44,45 B3PW91,^44,46 B97-D,⁴⁷ M05,⁴⁸ M05-2X,^49,50 M06,^49,51 M06-2X,^49,51 M06-HF,^49,52 M11,⁵³ MN12SX,⁵⁴ mPW1LYP,⁵⁵ mPW1PBE,⁵⁵ mPW1PW91,⁵⁵ mPW3PBE,⁵⁵ N12SX,⁵⁴ PBE0,⁵⁶ PBEPBE,⁵⁷ TPSSh,^58,59 TPSSTPSS,⁶⁰ and ωB97X-D.⁶¹ 

In addition, the 22 geometries in Fig. 2 were fully characterized in their reported point groups with two ab initio methods, second-order Møller–Plesset perturbation theory (MP2)^62–67 and the coupled cluster method that includes all single and double substitutions, as well as a perturbative treatment of the connected triple excitations [CCSD(T)].^68,69 MP2 optimizations and corresponding unscaled harmonic vibrational frequency computations were performed with the Gaussian 09 software package, whereas the corresponding CCSD(T) computations were performed with the CFOUR⁷⁰ software package.

For DFT and MP2 computations, the following 12 correlation consistent basis sets of double-, triple-, and quadruple-ζ quality were used: cc-pVXZ,⁷¹ aug-cc-pVXZ,^72,73 heavy-aug-cc-pVXZ,^72,73 and aug-cc-pV(X+d)Z⁷⁴ (where X = D, T, Q). The heavy-aug-cc-pVXZ basis set incorporates diffuse functions for the heavy (i.e., non-hydrogen) atoms (i.e., cc-pVXZ for H and aug-cc-pVXZ for all other atoms). The aug-cc-pV(X+d)Z basis sets contain a set of tight-d functions for second-row atoms. The basis sets cc-pVXZ, aug-cc-pVXZ, heavy-aug-cc-pVXZ, and aug-cc-pV(X+d)Z are abbreviated as XZ, aXZ, haXZ, and a(X+d)Z, respectively, where X is the cardinal number of the basis set. For CCSD(T) optimizations, the following basis sets were used: aug-cc-pVXZ, aug-cc-pV(X+d)Z, and aug-cc-pCVXZ⁷⁵ (where X = D, T), and CCSD(T) frequency calculations were only performed on geometries obtained with the aug-cc-pCVXZ basis sets (where X = D, T). The aug-cc-pCVXZ basis set, denoted as aCVXZ, includes core functions for non-hydrogen atoms to recover core–core and core–valence correlation effects.

Single-point energy computations were performed on the CCSD(T)/aCVTZ optimized geometries. The canonical MP2 and CCSD(T) methods and the explicitly correlated MP2-F12⁷⁶ and CCSD(T)-F12³⁵ methods were used in conjunction with the aXZ and aCVXZ series of basis sets (where X = D, T, Q, 5, or sometimes 6). In addition, for sulfur-containing systems, the a(X+d)Z, and aCV(X+d)Z series of basis sets (where X = D, T, Q, and sometimes 5 or 6 were used). The explicitly correlated F12 methods were used in conjunction with the VXZ-F12⁷⁷ series of basis sets (where X = D, T, Q). The aXZ, a(X+d)Z, and VXZ-F12 single-point energy computations were performed with Molpro 2012.1⁷⁸ and Molpro 2015.1.⁷⁹ The aCVXZ single-point energy computations were run using CFOUR and Molpro 2015.1, while the aCV(X+d)Z single-point energy computations were performed with Gaussian 16.⁸⁰ All single-point energies performed in Molpro utilized basis sets as implemented in Molpro. The aCV(X+d)Z basis sets were taken from the Psi4⁸¹ basis set library and can be found in the supplementary material (Tables S87–S90). All explicitly correlated computations include the default density fitting (DF) and resolution of the identity (RI) basis sets implemented in Molpro 2012.1 and Molpro 2015.1. The canonical and explicitly correlated MP2 and CCSD(T) energies were used to estimate the electronic energy at the complete basis set (CBS) limit via established extrapolation procedures. The three-parameter exponential function by Feller is used to extrapolate to the Hartree–Fock (HF) CBS limit,⁸² 

Equation (1) can then be reduced to an algebraic form that includes the three largest correlation consistent basis sets used in this study,⁸³ 

To describe the convergence of the correlation energy (E_corr) to the CBS limit, the two-parameter cubic function of Helgaker et al.⁸⁴ is used,

Equation (3) can also be written in an algebraic form for the two largest correlation consistent basis sets used in this study,⁸³ 

To obtain the CBS limit of the MP2 and CCSD(T) electronic energies, $EHFCBS$ is added to the corresponding $EcorrCBS$⁠.

For purposes of basis set testing, additional single-point energy computations were performed for the H₂O₂, H₂S₂, and H₂OS structures on the CCSD(T)/aCVTZ optimized geometries using CCSD(T) in conjunction with various combinations of the aXZ basis set for oxygen and sulfur and the a(X+d)Z basis sets for sulfur (where X = D, T, Q, 5). For example, for the H₂OS isomers, four sets of computations were performed: (1) The a5Z basis set was used for sulfur and the basis set for oxygen was varied with the aXZ basis sets (where X = D, T, Q), (2) the a(5+d)Z basis set was used for sulfur and the basis set for oxygen was varied with the aXZ basis set (where X = D, T, Q), and (3) the a5Z basis set was used for oxygen and the basis set for sulfur was varied with either the aXZ or (4) a(X+d)Z basis sets (where X = D, T, Q). For these computations, the aDZ basis set was used for hydrogen and symmetry was disabled (because of the use of mixed basis sets). Central and terminal labels for pyramidal structures are defined in Fig. 3 to denote the choice of basis sets for these atoms.

The 22 DFT methods used in this investigation were used to determine additional single-point energies on the CCSD(T)/aCVTZ optimized structures in conjunction with the aQZ, a(Q+d)Z, def2-QZVPD,⁸⁵ and aug-pc-3^86–88 basis sets. Five DFT methods (i.e., B3LYP, ωB97X-D, mPW1LYP, M06, MN12SX) were used to determine additional single-point energies on the CCSD(T)/aCVTZ optimized structures with an OS-core with the following basis sets: 6-311++G(2d,p),^89–92 6-311++G(2df,2dp),^92,93 6-311++G(3df,3pd),⁹¹ def2-TZVPD,⁸⁵ def2-TZVPPD,⁸⁵ def2-QZVPPD,⁸⁵ aTZ, a(T+d)Z, aug-pc-2,^86–88 aug-pcseg-2,⁸⁶ and aug-pcseg-3.⁸⁶ 

The correlation consistent composite approach (ccCA)^94,95 that is recommended for systems containing main group elements was also applied. For the ccCA energetics, the following single-point energy computations were performed at the HF⁹⁶/aug-cc-pVXZ (X = D, T, Q), MP2/aug-cc-pVXZ (X = D, T, Q), MP2(FC1)/aug-cc-pCVTZ, CCSD(T)/cc-pVTZ, and MP2/cc-pVTZ-DK^71,97 levels of theory on the B3LYP/cc-pVTZ optimized geometries. The energy described by ccCA is represented by the following equation:

Given the use of MP2 to describe core–core and core–valence effects in ccCA, the correlation consistent composite approach for transition metals (ccCA-TM)^94,98,99 was also considered as it uses CCSD(T) rather than MP2 to characterize these effects. For the ccCA-TM energetics, the following single-point computations were performed at the HF/aug-cc-pVXZ-DK^72,97 (X = D, T, Q), MP2/aug-cc-pVXZ-DK (X = D, T, Q), MP2/cc-pVTZ-DK, CCSD(T)/cc-pVTZ-DK (X = D, T), CCSD(T)/aug-cc-pwCVDZ-DK,¹⁰⁰ CCSD(T,FC1)/aug-cc-pVDZ-DK, and CCSD(T,FC1)/aug-cc-pwCVDZ-DK levels of theory on the B3LYP/cc-pVTZ optimized geometries. The energy described by ccCA-TM is determined by the following equation:

E_ref(ccCA) in Eqs. (5) and (6) represents the energy at the HF CBS limit added to the correlation energy at the MP2 CBS(aDTQ) limit. In a study by Halkier et al.,¹⁰¹ it was observed that the b parameter in Eq. (1) was roughly constant in their test of main group species and found that a value of 1.63 worked well for two-point extrapolations with smaller correlation consistent basis sets. The HF CBS limit can be computed with a two-point extrapolation of HF energies with the aTZ and aQZ basis sets with the following equation:

Equation (7) can also be written in an algebraic form as follows:¹⁰² 

$ΔECC$ is a correction that implements CCSD(T) to account for higher order dynamic correlation effects. $ΔECV$ represents a basis set correction that is used to describe the core–core and core–valence electron interactions. $ΔEZPE$ uses harmonic vibrational frequencies with a scaling factor of 0.989 and is the result of zero-point energy and thermal corrections at 298.15 K. $ΔESR$ is the scalar relativistic correction, it is only in the ccCA equation because all the terms in ccCA-TM include the effects due to scalar relativity. $ΔESO$ is the atomic spin–orbit coupling correction. For all ccCA and ccCA-TM computations involving sulfur-containing systems, the basis sets including the tight-d functions [i.e., aCV(X+d)Z-DK] were used instead of aCVXZ-DK based on the recommendations for compounds containing second-row elements from Wilson and co-workers.⁹² 

The G3B3¹⁰³ composite methodology was also applied on the B3LYP/cc-6-31G(d) optimized geometries where the following single-point computations were performed QCISD(T,FC)¹⁰⁴/6-31G(d),^105,106 MP4^107,108(FC)/6-31+G(d),^91,92 MP4(FC)/6-31G(2df,p),⁹¹ and MP2(full)/G3large.^109,110

Spherical harmonic basis functions (5d, 7f, etc.) were used rather than their Cartesian counterparts (6d, 10f, etc.) for all computations. The magnitudes of the components of the residual Cartesian gradients of the optimized structures were less than $2.14×10−4Eha0−1$⁠. The frozen core approximation was invoked for all computations (1s, 2s, 2p-like orbitals on oxygen and sulfur and 1s-like orbitals on carbon) with the XZ, aXZ, haXZ, and a(X+d)Z basis sets, while all electrons were correlated (i.e., full) in computations implementing the aCVXZ, aCV(X+d)Z, and aCV(X+d)Z-DK basis sets. The aCV(X+d)Z-DK basis sets were constructed from the cc-pV(X+d)Z-DK basis sets and the aug-cc-pCVXZ-DK basis sets. These basis sets can be found in the supplementary material (Tables S66–S68), and the changes are highlighted in red. Additionally, all MP2 data are given in the supplementary material (Tables S2, S3, S31, and S44). We note that the MP2 convergence behavior is similar to CCSD(T) results. The MP2 CBS limits deviate from the CCSD(T) CBS limits from around 0.5 to 5 kcal mol⁻¹. It is suggested that the T1 diagnostic should be less than or equal to 0.02,¹¹¹ the D1 diagnostic should be less than or equal to 0.05,¹¹² and the T1/D1 ratio should be less than or equal to 0.40.¹¹³ All the T1 and D1 diagnostic values as well as the T1/D1 ratios in this study (given in Tables S75–S86 of the supplementary material) are within or reasonably close to the recommended values.

All reported relative energies in this study are defined as the energy of a pyramidal isomer relative to an appropriate bent structural isomer. The CCSD(T) and corresponding explicitly correlated F12 relative energies with a variety of basis sets for select CCSD(T)/aCVTZ optimized geometries are plotted in Fig. 4 to illustrate the convergence behavior of these systems to the CCSD(T) CBS limit. The arrangement of the nine plots in Fig. 4 correlates with the layout of the structures given in Fig. 2, and the y-axes across all nine plots are consistently sized.

The CCSD(T) relative energies for all of the CCSD(T)/aCVTZ optimized geometries are given in Table II along with the corresponding CCSD(T) CBS limits. The bold numbers in Table II represent the best computed relative energy for each isomer pair and will serve as the benchmark values in this study. We note that the best CBS limits are determined using a(Q+d)Z, a(5+d)Z, and a(6+d)Z results [e.g., CBS(aQ56+d)] except in the case of isomers with an OO-core where CBS[CV](aTQ5) limits are used. CBS(aQ56) extrapolations were not performed for all the OO-core structures as the relative energies were converged at the aTQ5 level of theory as demonstrated with 1/2 where the a5Z and a6Z relative energies (in kcal mol⁻¹) agreed to two decimal places. Finally, all CCSD(T)-F12 relative energies are given in Table III.

The relative energies for these systems decrease as hydrogen substituents are systematically substituted with methyl groups (left to right in Fig. 4). Looking at the differences in the best CBS limits (Table II), the relative energies of compounds with an OO-core (row 1) decrease by 2.44 kcal mol⁻¹ as one hydrogen is substituted (3/4) and by 3.44 kcal mol⁻¹ as both hydrogens are substituted (5/6). Compounds with an SS-core (row 3) show larger decreases in the relative energies of 6.29 kcal mol⁻¹ and 7.22 kcal mol⁻¹ as one hydrogen is substituted (19/20) and both hydrogens are substituted (21/22), respectively.

The compounds with an OS-core (row 2) show a more dramatic decrease in relative energies than their OO- or SS-core counterparts (rows 1 and 3, respectively). The relative energy of 10 and 12 (Fig. 4) and 14 and 16 (Fig. 4) qualitatively changes when characterized with larger basis sets with the pyramidal conformers being lower in energy than their bent counterparts.

Compounds with an OO-core (the first row of Fig. 4) show rapid convergence to the CCSD(T) CBS[CV](aTQ5) limit. The deviations from the CBS[CV](aTQ5) limit with the smallest basis set in this series (i.e., aDZ) are 1.06 (1/2), 0.97 (3/4), and 1.29 (5/6) kcal mol⁻¹. Compounds with an SS-core (the third row of Fig. 4) show fast, smooth convergence to the CCSD(T) CBS(aQ56+d) limit. The deviations from the CBS(aQ56+d) limit with an aDZ basis set are 4.07 (17/18), 3.89 (19/20), and 3.56 (21/22) kcal mol⁻¹. A less than 1 kcal mol⁻¹ deviation is observed with an aQZ basis set, and these deviations are 0.99 (17/18), 0.98 (19/20), and 0.94 (21/22) kcal mol⁻¹.

Surprisingly, compounds with an OS-core (the second row in Fig. 4) exhibit the pathological behavior of extremely slow convergence of the relative energy to the CCSD(T) CBS(aQ56+d) limit. These systems exhibit large relative energy deviations as large as 5.86 kcal mol⁻¹, even at the CCSD(T)/aTZ level of theory. The deviations from the CBS(aQ56+d) limit with an aDZ basis set are 12.52 (7/9), 13.92 (10/12), and 13.96 (14/16) kcal mol⁻¹, an order of magnitude larger when compared to compounds with an OO- or SS-core. A less than 1 kcal mol⁻¹ deviation from the benchmark CBS limit is not observed until the a5Z basis set is used, necessitating the need for a Q56 extrapolation to obtain converged relative energies for these systems. The deviations from the CBS(aQ56+d) limit with an a5Z basis set are 0.80 (7/9), 0.82 (10/12), and 0.80 (14/16) kcal mol⁻¹. Additionally, the relative energy for 10/12 qualitatively changes from positive to negative with the a6Z basis set.

The inclusion of tight-d functions in the basis set for sulfur accelerates the convergence of the relative energies to the CCSD(T) CBS(aQ56+d) limit for the sulfur-containing structures. This acceleration is more dramatic in the OS-core compounds. The deviation from the CBS limit with an a(D+d)Z basis set is approximately halved when compared to the aDZ basis set. For example, compounds with an SS-core show deviations from the CCSD(T) CBS limit of 4.07 (17/18), 3.89 (19/20), and 3.56 (21/22) kcal mol⁻¹ with an aDZ basis set and only 2.27 (17/18), 2.05 (19/20), and 1.75 (21/22) kcal mol⁻¹ with an a(D+d)Z basis set. The deviations drop below 1 kcal mol⁻¹ with an a(T+d)Z basis set and are 0.91 (17/18), 0.93 (19/20), and 0.94 (21/22) kcal mol⁻¹.

For compounds with an OS-core, the deviations from the CCSD(T) CBS limit are 12.52 (7/9), 13.92 (10/12), and 13.96 (14/16) kcal mol⁻¹ with an aDZ basis set and only 5.87 (7/9), 7.22 (10/12), and 7.30 (14/16) kcal mol⁻¹ with an a(D+d)Z basis set. The deviations drop to 1 kcal mol⁻¹ with an a(Q+d)Z basis set and are 0.94 (7/9), 1.03 (10/12), and 1.00 (14/16) kcal mol⁻¹. The inclusion of tight-d functions leads to deviations of about 1 kcal mol⁻¹ from the CCSD(T) CBS limit with a basis set that is one cardinal number smaller than the corresponding standard basis set.

The inclusion of core basis functions to recover core–core and core–valence correlation effects in a fully correlated computation shows trivial differences on the relative energetics of compounds with an OO-core with respect to the standard valence basis sets (i.e., aXZ). For compounds with an SS- or OS-core, these basis sets mildly accelerate the convergence to the CBS limit with respect to the standard valence basis sets; however, computed relative energies show a similar convergence behavior seen with the a(X+d)Z basis sets. In fact, a(D+d)Z results are dramatically improved over their aCVDZ counterparts.

Similar relative energetics between the aCVXZ and a(X+d)Z basis sets appear when X = T, Q, and 5. For all SS- and OS-core compounds, the differences between the aCVXZ and a(X+d)Z basis sets (where X = T, Q, 5) were determined. The average absolute value of the 33 differences (11 isomer pairs and 3 basis sets) is 0.19 kcal mol⁻¹. These results demonstrate that core–core and core–valence effects on the relative energies are minimal in these sulfur-containing systems at an increased computational cost.

The inclusion of both core and tight-d functions in a fully correlated computation for sulfur-containing compounds gives similar relative energies as those obtained with a(X+d)Z basis sets with only minor improvement. For example, the deviation from the CCSD(T) CBS limit for the relative energies of 7 and 9 with an a(D+d)Z basis set is 5.87 kcal mol⁻¹ and with an aCV(D+d)Z basis set is 5.77 kcal mol⁻¹. For all SS- and OS-core compounds, the differences between the a(X+d)Z and aCV(X+d)Z basis sets (where X = D, T, Q) were determined. The average absolute value of the 33 differences (11 isomer pairs and 3 basis sets) is 0.25 kcal mol⁻¹.

The aCV(X+d)Z basis sets provide only marginally better results when compared to the a(X+d)Z basis sets overall, but at significantly increased computational cost. Even though the aCV(D+d)Z basis set outperforms the a(D+d)Z basis set in the computed relative energies with respect to the CCSD(T) CBS limit, the results are still unreliable with deviations from the corresponding benchmark values exceeding 7 kcal mol⁻¹ as seen with DMSO (14) and 16. These accuracy improvements quickly diminish as the cardinal number of the basis set is increased with respect to the a(X+d)Z basis sets used in computations invoking the frozen core approximation.

The relative energies calculated at the CCSD(T) CBS limits are given at the bottom of Table II. For OO- and SS-core compounds, the CCSD(T)/CBS relative energies are usually 0.1 kcal mol⁻¹ or less from the benchmark values (given in bold) although they can grow as large as −0.50 kcal mol⁻¹ as seen for 17/18 with CBS(aTQ5). All CCSD(T)/CBS results incorporating the tight-d function [i.e., a(X+d)Z] for sulfur are within 0.07 kcal mol⁻¹ of the benchmark values.

The CCSD(T) CBS limits for OS-core compounds illustrate the extent of the challenging nature of characterizing the relative energies of these systems with standard aXZ basis sets. For example, the CBS(aDTQ) relative energy of 7 and 9 is 1.86 kcal mol⁻¹ higher than the CBS(aQ56+d) benchmark result, whereas the CBS(aTQ5) is −1.86 kcal mol⁻¹ below the benchmark value (left in Fig. 5). The relative energies tabulated in Table II demonstrate interesting pathological behavior. For the relative energies of isomer pairs 7/9, 8/9, 10/12, and 14/16, the deviations of the CCSD(T)/CBS(aDTQ) either over- or under-estimate the CCSD(T)/CBS(aQ56+d) limit by +1.86, −1.01, +1.68, and +1.59 (kcal mol⁻¹). For this same set of compounds, the CCSD(T)/CBS(aTQ5) results either under- or over-estimate the CCSD(T)/CBS(aQ56+d) limit by −1.86, +0.42, −1.67, and −1.48 (kcal mol⁻¹). The CBS(aQ56) limit lies within 0.14 kcal mol⁻¹ of the CBS(aQ56+d) benchmark. Similar deviations are seen in 10/12, 10/13, and 14/16 with these basis sets. All CBS extrapolations incorporating tight-d functions or core functions in the basis set show smaller deviations (less than 1 kcal mol⁻¹) from the benchmark. The CBS[CV] limits deviate from the benchmark values as much as 0.5 kcal mol⁻¹ as seen with CBS[CV](aDTQ) for 7 and 9 (left in Fig. 5) but show better agreement with the benchmark for DMSO (14) and 16 (right in Fig. 5).

It is known that explicitly correlated methods accelerate the convergence to the CBS limit. An F12 computation with a particular basis set, XZ, has been demonstrated to give results that have similar quality to a standard method with an augmented basis set two cardinal numbers larger [i.e., (X+2)Z].¹¹⁴ For compounds with an OO-core, the explicitly correlated CCSD(T)-F12 method with an aDZ basis set gives a deviation of about 0.5 kcal mol⁻¹ from the CCSD(T) CBS limit and gives very similar results relative to the canonical CCSD(T) method as the size of the basis set increases. For compounds with an SS- or OS-core, the CCSD(T)-F12 method significantly accelerates the convergence to the CBS limit. The average absolute deviations (AADs) from the CCSD(T) CBS limit for all the sulfur-containing compounds are quite small and are 0.71 kcal mol⁻¹ with an aDZ basis set and 0.46 kcal mol⁻¹ with an a(D+d)Z basis set. These AADs decrease to 0.27 kcal mol⁻¹ with an aTZ basis set and 0.16 kcal mol⁻¹ with an a(T+d)Z basis set.

The VXZ-F12 basis sets perform similarly to the aXZ basis sets for compounds with an OO- and SS-core. In fact, the maximum deviation from the corresponding CCSD(T) CBS limit for OO- and SS-structures did not exceed 0.8 kcal mol⁻¹ for the VDZ-F12 basis set. For the OS-core compounds, this maximum deviation grew to 1.70 kcal mol⁻¹ as seen with 8/9.

Overall, the explicitly correlated CCSD(T)-F12 method brings significant improvement to the relative energies of these systems. When coupled with the aDZ or a(D+d)Z double-zeta basis sets, the CCSD(T)-F12 method gave relative energies that were within 1.00 kcal mol⁻¹ of the benchmark CBS limit. VDZ-F12 basis sets gave deviations as large as 1.70 as seen with 8 and 9.

The relative energies of OS-core systems are challenging to characterize relative to the OO- or SS-core systems as they exhibit the pathological behavior of extremely slow convergence to the CCSD(T) CBS(aQ56+d) limit, deviations that are an order of magnitude larger than the OO- or SS- counterparts, and a significant basis set dependence. To gain insight into why the OS-core systems behave pathologically, additional CCSD(T) single-point energy computations were performed on the CCSD(T)/aCVTZ optimized geometries using various combinations of the aXZ and a(X+d)Z basis sets for the H₂O₂, H₂S₂, and H₂OS structures (because the substituents H and CH₃ demonstrate a similar convergence behavior). These single-point computations apply an a5Z basis set to oxygen or sulfur or an a(5+d)Z basis set to sulfur and a series of aXZ basis sets (X = D, T, Q, 5) to the second heavy atom (an aDZ basis set is used for H throughout this series). The absolute and relative energies are given in Tables S92–S97 of the supplementary material. The deviations of these energies from the best CBS limit (Table II) for the OO-, SS-, and OS-core isomers are given in Tables IV–VI, respectively. With the OO-core, no pathological behavior is observed. For all mixed basis set combinations, the OO-core isomers give relative energies that deviate by no more than 1.74 kcal mol⁻¹ (Table IV) when compared to the CCSD(T) CBS[CV](aTQ5) limit.

The SS-core isomers reveal that the choice of basis set is critical to characterizing these relative energies; the basis set on the central sulfur appears to be most sensitive to this basis set choice (Table V). When the basis set on the terminal sulfur is held at a5Z and the basis set for the central sulfur is increased from aDZ to aQZ, the relative energy deviation from the CCSD(T) CBS(aQ56+d) limit decreases from 3.64 to 1.50 kcal mol⁻¹. When the central sulfur has a5Z and the terminal sulfur is increased from aDZ to aQZ, the maximum deviation is 1.77 kcal mol⁻¹ and the minimum deviation is 0.53 kcal mol⁻¹. When tight-d functions are added to both sulfurs, first keeping the terminal sulfur at a(5+d)Z and increasing the basis set for the central sulfur from a(D+d)Z to a(Q+d)Z, the deviation decreases to a maximum of 1.20 kcal mol⁻¹. In the opposite case, when keeping the central sulfur at a(5+d)Z and increasing the basis set for the terminal sulfur from a(D+d)Z to a(Q+d)Z, the deviation decreases to a maximum of 0.99 kcal mol⁻¹. In all cases, the results demonstrate that the basis set on the central sulfur is more important to obtaining the smallest relative energy deviation from the CCSD(T) CBS(aQ56+d) limit.

The choice of basis set is even more critical for characterizing the relative energies of the OS-core isomers (Table VI). The first set of cases to consider is when no tight-d functions are utilized in the basis sets. For 7/9 or 8/9, when the basis set on oxygen is essentially saturated (a5Z) and the basis set on sulfur is varied from aDZ to aQZ, the deviations range from 6.88 to 1.73 kcal mol⁻¹. When the basis set on sulfur is essentially saturated (a5Z) and the basis set on oxygen is varied from aDZ to aQZ, the results are vastly improved (1.83–0.14 kcal mol⁻¹). The addition of tight-d functions to sulfur improves the results, although not in a satisfying way. Using a(D+d)Z on sulfur and a5Z on oxygen yields fortuitously close relative energy differences of 0.51 (7/9) and 0.84 (8/9) kcal mol⁻¹. However, in this case, increasing the basis set on sulfur does not improve the results. The alternative combination, saturating sulfur with a(5+d)Z and using aDZ on oxygen, yields improvement, but less satisfying results: relative energy differences of 2.21 (7/9) and 1.95 (8/9) kcal mol⁻¹. Increasing the basis set on oxygen does improve the results.

The average absolute deviation (AAD) of the relative energies from the CCSD(T) CBS limit of the isomer pairs from fully characterized geometries determined with 22 density functionals and 12 basis sets (categorized by an OO-, OS-, or SS-core) is given in Fig. 5. All DFT relative energies are given in the supplementary material (Tables S31–S44). For all the isomer pairs in this study, the smallest basis set considered (i.e., DZ) gives unreliable relative energies across all 22 density functionals with deviations from the CCSD(T) CBS limit that are generally 6 kcal mol⁻¹ or larger.

For compounds with an OO-core (left in Fig. 6), the AADs generally decrease as the size of the basis set is increased from double-ζ quality to quadruple-ζ quality. The inclusion of diffuse functions on heavy atoms (i.e., haXZ) dramatically improves the deviations in general. The inclusion of diffuse functions on all the atoms (i.e., aXZ) also improves the deviations in general but shows marginal improvement with respect to haXZ results. The smallest basis set that gives a deviation of about 1 kcal mol⁻¹ from the CCSD(T) CBS limit for most density functional approximations considered is the aDZ basis set.

Compounds with an SS-core (right in Fig. 6) have AADs that behave similarly to those seen in the compounds with an OO-core with respect to the haXZ and aXZ series of basis sets. The inclusion of tight-d functions in the basis set for sulfur have the most dramatic improvement, in general, on the AADs of the relative energies. These results indicate the importance of tight-d functions on the relative energies for these sulfur-containing systems. In fact, 17 of the 22 methods with an a(D+d)Z basis set give AADs that are less than 1 kcal mol⁻¹ from the CCSD(T) CBS limit.

Compounds with an OS-core (middle in Fig. 6) are challenging to accurately characterize using density functional approximations. Basis sets of double-ζ quality consistently give unreliable relative energies. The corresponding AADs for these systems do generally decrease as the size of the basis set is increased from double-ζ quality to quadruple-ζ quality as seen with the OO- and SS-core compounds. The inclusion of tight-d functions in the basis set for sulfur decreases the AADs relative to the other family of basis sets used. However, the resulting AADs are still quite large overall. All the AADs for DFT/a(D+d)Z results on OS-core compounds never fall below 1 kcal mol⁻¹. In fact, only three AADs fall below 1 kcal mol⁻¹ with an a(T+d)Z basis set and only four AADs fall below 1 kcal mol⁻¹ with an a(Q+d)Z basis set. These results indicate the challenging nature of accurately characterizing the relative energies of these sulfoxide systems.

Single-point energies of the CCSD(T)/aCVTZ optimized geometries were computed with the 22 DFT methods used in this investigation in conjunction with the aQZ, def2-QZVPD, and aug-pc-3 basis sets on the compounds with an OO-core. The a(Q+d)Z basis set was also used in addition to the three basis sets listed above for compounds with an SS-core. On average, these basis sets give deviations of about 1 kcal mol⁻¹ from the CCSD(T) CBS limit for these systems. All the relative energies are available in the supplementary material (Tables S47–S56).

For the challenging OS-core compounds, the 5 DFT methods that gave an AAD of about 1 kcal mol⁻¹ from the CCSD(T) CBS limit (i.e., B3LYP, ωB97X–D, mPW1LYP, M06, and MN12SX) were used in additional single-point energy computations on the CCSD(T)/aCVTZ optimized geometries with the following basis sets: 6-311++G(2d,p), 6-311++G(2df,2pd), 6-311++G(3df,3pd), aTZ, a(T+d)Z, aQZ, a(Q+d)Z, apc-2, apc-seg-2, apc-3, apc-seg-3, def2-TZVPD, def2-TZVPPD, def2-QZVPD, and def2-QZVPPD. The average absolute deviations (AADs) and maximum absolute deviations (MADs) for each level of theory are given in Table VII. The B3LYP and ωB97X–D methods give the best AADs overall, while ωB97X–D and MN12SX yielded the lowest MADs overall. At minimum, when using density functional approximations, robust triple-ζ basis sets are required to characterize the relative energies of these systems with minimal error such as augmented correlation consistent basis sets with tight-d functions. Interestingly, none of the triple-ζ quality polarization consistent basis sets (i.e., apc-2 and apc-seg-2) yield AADs that are less than 1 kcal mol⁻¹, whereas the def2 basis sets consistently lie below this threshold.

The ccCA, ccCA-TM, and G3B3 composite methodologies were applied to all 22 compounds given in Fig. 2, and the relative energies of each isomer pair are given in Table VIII. The deviations of these relative energies from the best CCSD(T) CBS limit are given in Table VIII. These composite methodologies consistently give deviations that are less than 1 kcal mol⁻¹ from the CCSD(T) CBS limit for compounds with an OO- or SS-core, and we note that G3B3 gives the smallest deviations on average (Table VIII).

For compounds with an OS-core, these composite methodologies give deviations that are less than about 2 kcal mol⁻¹ from the CCSD(T) CBS limit. None of these composite methods consistently perform better than the other. It could be expected that ccCA-TM would perform better than ccCA since ccCA-TM incorporates CCSD(T) rather than MP2 to describe the core valence effects. However, ccCA-TM uses CCSD(T) with an aCVDZ-DK basis set, while ccCA uses MP2 with an aCVTZ basis set to describe the core valence effects. This paper has demonstrated that a basis set of double-ζ quality is not of sufficient size to accurately describe the energetics of these OS-core compounds, which could explain why ccCA performs better than ccCA-TM.

The curious nature and pathological behavior of the relative energy of DMSO with respect to its bent structural isomer reveal the challenging nature of characterizing the relative energies of these small systems. This study systematically characterizes the relative energetics of small, bent molecules of the forms ROOR′, RSSR′, and ROSR′ (where R, R′ = H, CH₃) relative to their pyramidal structural isomer forms in order to provide a computational benchmark to determine which method and basis set to implement for these challenging systems. The substituents H and CH₃ have little effect on the convergence behaviors of these systems.

This study thoroughly demonstrates the pathological nature of the OS-core compounds and highlights the serious basis set dependence for sulfur (the behavior is slightly worse when the sulfur is in the central position in the pyramidal structure). This behavior can also be seen in the SS-core compounds; however, the effects are tempered.

Due to the significant basis set dependence, basis sets as large as sextuple-ζ are necessary to obtain reliable relative energies at the CCSD(T) CBS(aQ56+d) limit for OS-core systems. The inclusion of tight-d basis functions for sulfur provided by the a(X+d)Z basis sets significantly improves the corresponding relative energies with respect to the benchmark CBS limit. The inclusion of core functions for non-hydrogen atoms provided by the aCVXZ basis sets to recover core–core and core–valence correlation effects in fully correlated computations improves the computed relative energies with respect to computations invoking the frozen core approximation with the standard aXZ series of basis sets. However, these fully correlated approaches offer little advantage to using a(X+d)Z basis sets, given the higher computational cost. Including both core and tight-d functions provided by aCV(X+d)Z basis sets yielded trivial improvement over the a(X+d)Z series of basis sets. The augmented basis sets with only the inclusion of tight-d functions are recommended as they have the most significant impact on the accuracy of the relative energetics for the compounds with an SS- or OS-core. The results from the explicitly correlated CCSD(T)-F12 method, even with aDZ, closely agree with the CCSD(T) CBS limit results (consistently within 1 kcal mol⁻¹). To our delight, the explicitly correlated methods with aDZ basis sets outperform the canonical methods with aQZ basis sets. Given these results, CCSD(T)-F12 is recommended, surprisingly, with a double-ζ quality basis set for characterizing the energies of these systems.

Due to the qualitative and quantitative discrepancies, DFT energetics are generally unreliable. However, the following DFT methods with the a(T+d)Z basis sets are recommended for compounds with an OS-core: B3LYP, ωB97X-D, mPW1LYP, M06, and MN12SX.

All composite methodologies gave deviations that were within 1 kcal mol⁻¹ of the CCSD(T) CBS limit for compounds with an OO- and SS-core and deviations that were within 2 kcal mol⁻¹ for compounds with an OS-core. ccCA is recommended for compounds with an OO- and SS-core, and G3B3 is recommended for main group compounds with an OS-core.

See the supplementary material for all relative and absolute energies; Cartesian coordinates for all optimized geometries; the aCV(X+d)Z-DK constructed basis sets for sulfur; the aCV(X+d)Z basis sets from the Psi4 basis set library; common names for all 22 compounds in this study; the number of contracted basis functions for 21 basis sets used for 14, 15, and 16; the T1 and D1 diagnostics; and the T1/D1 ratios. Additionally, interactive CCSD(T) extrapolations and DFT heat maps are given in HTML format.

The authors acknowledge the Mississippi Center for Supercomputing Research (MCSR) for providing access to their computational resources. This work was funded, in part, by the National Science Foundation (Grant Nos. CHE 1800201 and OIA 1539035). The authors would also like to thank Professor Gregory Tschumper (University of Mississippi) for providing access to his computational resources for some of the CCSD(T)/a6Z and CCSD(T)/a(6+d)Z single-point energy computations, Professor Kirk Peterson (University of Washington) for his insights into the CBS nomenclature, and Professor Steven Gwaltney for stimulating conversations. L.N.O. thanks Dr. Robert Lamb for answering a million questions. C.E.W. thanks Professor Robert Pascal (Tulane University) for an off-hand comment that led to this work.

The data that support the findings of this study are available within the article and its supplementary material.