Exploring Avenues beyond Revised DSD Functionals: I. Range Separation, with xDSD as a Special Case

We have explored the use of range separation as a possible avenue for further improvement on our revDSD minimally empirical double hybrid functionals. Such ωDSD functionals encompass the XYG3 type of double hybrid (i.e., xDSD) as a special case for ω → 0. As in our previous studies, the large and chemically diverse GMTKN55 benchmark suite was used for evaluation. Especially when using the D4 rather than D3BJ dispersion model, xDSD has a slight performance advantage in WTMAD2. As in previous studies, PBEP86 is the winning combination for the semilocal parts. xDSDn-PBEP86-D4 marginally outperforms the previous “best in class” ωB97M(2) Berkeley double hybrid but without range separation and using fewer than half the number of empirical parameters. Range separation turns out to offer only marginal further improvements on GMTKN55 itself. While ωB97M(2) still yields better performance for small-molecule thermochemistry, this is compensated in WTMAD2 by the superior performance of the new functionals for conformer equilibria. Results for two external test sets with pronounced static correlation effects may indicate that range-separated double hybrids are more resilient to such effects.


INTRODUCTION
Kohn−Sham density functional theory (KS-DFT) 1 is presently by far the most widely used family of electronic structure methods. Its combination of reasonable accuracy and comparatively gentle computational cost scaling makes it an appealing choice for medium and large molecules; for small molecules, wavefunction ab initio (WFT) approaches still outperform it. 2,3 The accuracy of KS-DFT stands or falls with the exchange− correlation (XC) functional. Perdew 4 organized the plethora of available approaches into what he called a "Jacob's Ladder", arranged by the kinds of information employed in it: local density approximation (LDA) on the first rung, GGAs (generalized gradient approximations) on the second rung, meta-GGAs on the third rung (adding either the density Laplacian or the kinetic energy density), and hybrid functionals on the fourth rung (adding also the occupied orbital information). The fifth rung corresponds to the inclusion of virtual orbital information: the most widely used class of such methods are the so-called double hybrids (see refs 5−7 for reviews and most recently ref 8 by the present authors). As shown in refs, 8,9 their accuracy over the very large and diverse GMTKN55 (general main group thermochemistry, kinetics, and noncovalent interactions, with 55 problem sets) test suite 10 approaches that of WFT methods, yet the CPU cost increase over that of ordinary hybrid GGAs is actually quite modest if an RI (resolution of the identity 11,12 ) approximation is applied in the MP2 (second-order Møller−Plesset) part.
Generally speaking, there are two basic approaches available for double hybrids in the literature, which we shall denote gDH (after Grimme 13 ) and xDH (after the XYG3 functional 14,15 ) in the article. In gDH, an iterative KS calculation is carried out with a fraction (c X ′ ,HF ) of Hartree−Fock (HF) exchange and (1 − c X ′ ,HF ) of DFA (density functional approximation) exchange, plus the DFA correlation scaled by a coefficient c C,DFA . Next, using the converged orbitals from the KS step, a post-HF GLPT2 (second-order Gorling−Levy perturbation theory) 16 correlation energy term is evaluated on the basis of the KS orbitals and added in. (As with lower-rung DFT methods, a dispersion correction can optionally be added, though it generally needs a prefactor that is less than unity, since some dispersion is already captured in the GLPT2 term.) Double hybrids with nonlocal correlation terms other than PT2, such as the direct random phase approximation (dRPA, see ref 17 and references therein), are discussed by Kaĺlay and cow-orkers 18,19 as well as in the companion article 20 to the present work, while Manby and coworkers 21,22 very recently proposed a novel approach based on the Unsold approximation (UW12).
In contrast, for xDHs, the KS orbitals used for the evaluation of all energy terms at the final step are evaluated for a standard hybrid with the full DFA correlation (i.e., c C,DFA = 100%) and with c X as appropriate for a typical hybrid functional. It was argued 14,15 that such orbitals are more appropriate as a basis for GLPT2 than the damped-correlation orbitals in the gDH, though this argument has been refuted on empirical grounds by Goerigk and Grimme 23 and by Kesharwani et al. 24 = + + − where E N1e stands for the sum of nuclear repulsion and oneelectron energy terms; c X,HF and c C,XC are the fractions of exact exchange and semilocal correlation, respectively; E disp is the dispersion correction term (dependent upon parameters such as s 6 , s 8 , a 1 , a 2 , c ATM , and so on); and c 2ab and c 2ss are the two coefficients corresponding to the opposite-spin and same-spin GLPT2 correlation, respectively. The xDH version thereof, denoted xDSD, has been explored in ref 24 and found to offer only a minor advantage over the corresponding DSD. It must however be said that both the DSD and the xDSD functionals were originally parameterized and validated using quite modest training sets (for reasons of computational cost); furthermore, the weighting of the subsets is somewhat arbitrary, and experimentation on our part showed considerable dependence of the final parameters on the weights used there. In contrast, the much larger GMTKN55 dataset is not only over 10x larger but uses a more robust, unambiguously weighted performance metric in the guise of WTMAD2 (weighted mean absolute deviation, type 2, in which the weights of the subsets are corrected for the different energy scales of the reference data).
In ref 9, we were able to leverage GMTKN55 to obtain a family of more accurate revDSD functionals, with revDSD-PBEP86-D4 as the winner among them: just for the PBEP86 case, we also considered a single example of the xDH type and did find xrevDSD-PBEP86-D4 to be slightly more accurate still than revDSD-PBEP86-D4. One objective of the present article is to explore whether this is true more generally: specifically, we shall investigate xDSD-PBEPBE here; we also include xDSD-PBEPW91, xDSD-PBEB95, xDSD-BLYP, and xDSD-SCAN in our arsenal. Two types of dispersion correction, D3(BJ) 27,28 and more recent, flexible, and accurate D4, 29,30 will be considered, the latter with different many-body dispersion terms also. (We will also consider xDOD forms, in which same-spin GLPT2 is eliminated: this permits further acceleration for large systems through a Laplace transform algorithm. 31−36 ) The second objective is to investigate whether revDSD can be improved through introducing the range-separated HF exchange (RSH). In the long-distance limit, the exchange potential of global hybrids (GHs) behaves 37 like −c X /r 12 rather than the correct −1/r 12 term (r 12 being the interelectronic distance). Hence, Hirao and coworkers 38 proposed a scheme where the interelectronic repulsion operator 1/r 12 is partitioned into a short-range (SR) component to be treated by a (meta)GGA and a long-range (LR) component to be treated by the "exact" exchange, and to "cross-fade" from the SR to LR component using an error function (erf x) and the complementary error function (erfc x = 1 − erf x) of r 12 In this equation, ω represents the range separation parameter, which controls the transition between the LR and SR parts, α is the percentage of "exact" HF exchange in the short-range limit, and α + β is the corresponding percentage in the long-range limit. (Proper asymptotic behavior can be enforced through 39,40 α + β = 1/ε 0 , where the dielectric constant ε 0 = 1 in vacuoleading to β = 1 − αand ε 0 → ∞ for a perfect conductor.) ω can be determined empirically using a training set 37,41−45 or tuned non-empirically by minimizing the deviation from the conditions the exact KS functional must obey. 46,47 Following this approach, several empirical and non-empirical LC-DH functionals have been proposed, such as LC-PBE, 43 LC-ωPBE, 48 M11, 44 CAM-B3LYP, 37 ωΒ97, 45 ωB97X, 45 ωB97X-V, 49 ωB97M-V, 50 and many more.
We shall denote DSD-type double hybrids with RSH functionals ωDSD, where ω stands for the range separation parameter. Note that for ω = 0, ωDSD and ωDOD functionals reduce to the xDSD and xDOD forms, respectively, which ties our two objectives together.
The combination of range separation with GLPT2 for the correlation energy was first proposed by Ángyań and coworkers. 51 Chai and Head-Gordon 52 instead obtained orbitals from an RSH calculation and then evaluated the GLPT2 correlation on the basis of these orbitals, the final energy being a mix of the GGA exchange, HF exchange, GGA correlation, and GLPT2 correlation. Their most recent elaboration of this concept was the ωB97M(2) functional, 53 which for the GMTKN55 benchmark was found to have the lowest WTMAD2 of all functionals surveyed. 8 (To be fair, however, it has three times the number of empirical parameters of the next best performer, xrevDSD-PBEP86-D4. 8
After some experimentation, we settled on diet100 for the prescreening stage: based on this, we will decide which semilocal functionals to retain for in-depth investigation with the full GMTKN55 set. GMTKN55 10 is the updated and larger form of the Grimme group's previous GMTKN24 60 and GMTKN30 23 databases. This dataset consists of 55 types of chemical problems, which can be further categorized into five top-level subsets: thermochemistry of small-and medium-sized molecules, barrier heights, large molecule reactions, intermolecular interactions, and conformer energies. One full evaluation of the GMTKN55 needs 2459 single-point energy calculations (give or take a few duplicates) to generate 1499 unique energy differences. (Complete details of all 55 subsets and original references can be found in Table S1 in the Supporting Information.) Originally proposed by Goerigk et al., 10 WTMAD2 has been used as the primary metric for this work: where E i |Δ | is the mean absolute value of all the reference energies from i = 1 to 55, N i is the number of systems in each subset, and MAD i is the mean absolute difference between calculated and reference energies for each of the 55 subsets. Note that, from the statistical viewpoint, MAD (mean absolute deviation) is a more "robust" metric than rmsd (root-meansquare deviation) 61 as MAD is more resilient to a small number of large outliers than rmsd. For a normal distribution without a systematic error, rmsd ≈ 5MAD/4. 62 As one reviewer pointed out, the average absolute reaction energies for NBPRC and MB16-43 provided in the original GMTKN55 article 11 differ from the corresponding values calculated from the individual data supplied in the Supporting Information. If these corrected average absolute reaction energies were employed in the construction of eq 3, then their average, which appears in eq 3 as the overall scale factor, would be 57.76 rather than 56.84. However, as all previously published articles on GMTKN55 (such as refs 8,9,28,63−66 ) have used the original (56.84) coefficient, we are also retaining it for the sake of compatibility. It goes without saying that this will not affect the ranking between functionals; those who prefer WTMAD2 57.76 can simply multiply all WTMAD2 values by 1.0162.
Reference geometries were taken "as is" from ref 10 and not optimized further.
2.2. Electronic Structure Calculations. All the calculations were performed using the Q-CHEM 5.3 67 package (except ωB2GP-PLYP 58 and ωB2PLYP, 58 for which ORCA 4.2.1 68 has been used), running on the ChemFarm HPC cluster of the Weizmann Institute Faculty of Chemistry.
The Weigend−Ahlrichs 69 def2-QZVPP basis set was considered throughout with a few exceptions, such as the WATER27, RG18, IL16, G21EA, AHB21, BH76, and BH76RC subsets, where the diffuse-function-augmented def2-QZVPPD 70 basis set was used instead. However, for the computationally demanding C60ISO and UPU23 subsets, which have small weights in WTMAD2, the more economical def2-TZVPP 69 basis set was employed to curb the computational cost. The SG-3 71 integration grid was used across the board, except for the SCAN (strongly constrained and appropriately normed 72 meta-GGA type) variants, where due to SCAN's severe integration grid sensitivity, 73 an unpruned (150, 590) grid was employed. In the MP2-like step, the RI approximation was applied in conjunction with the def2-QZVPPD-RI fitting basis set. 74,75 For the ωB2GP-PLYP 58 and ωB2PLYP 58 functionals using ORCA, we have used the JK auxiliary basis set for Coulomb and exchange RI integrals (def2/JK). 76 In this project, while most of the calculations were completed using frozen inner-shell orbitals, we made two departures from this recipe to avoid unacceptably small orbital energy gaps between the highest frozen and lowest correlated orbitals. First, for the MB16-43, HEAVY28, HEAVYSB11, ALK8, CHB6, and ALKBDE10 subsets, we correlated the (n − 1)sp subvalence electrons of the metal and metalloid atoms. Second, for HAL59 and HEAVY28, the (n − 1)spd orbitals of the heavy p-block elements were kept unfrozen. Note that unlike the valence correlation consistent basis sets, the Weigend−Ahlrichs QZVPP basis set is multiple-zeta in the core as well and contains core−valence polarization functions (see Table 1 of ref 69). Semidalas and Martin (in the context of composite wavefunction calculations) considered 3 the impact of the core−valence correlation on GMTKN55 using correlation-consistent core−valence basis sets and found that its impact is on the order of 0.05 kcal/molwhich will be further reduced here as the PT2 correlation terms are scaled down in a double hybrid. Powell's BOBYQA (bound optimization BY quadratic approximation) derivative-free constrained optimizer 79 and a few scripts and Fortran programs developed in-house were used for optimizing the parameters.
Once a full GMTKN55 evaluation is finished with a fixed set of {c X,HF , c C,DFT , ω}, no further electronic structure calculation is needed to get an associated optimal set of (c−f); the latter set of parameters can be obtained in what amounts to an inner optimization loop, whereas c X,HF , c X,DFT , and ω (where applicable) can be minimized in an outer optimization loop. (We previously found in the revDSD article 9 that the coupling between (a) and (c,d) is too strong to permit placing c C,DFT in the inner loop and that for a fixed value of (a) convergence of the DFT correlation parameter, up to two decimal places can be achieved within two macroiterations.) The process is analogous to microiterations versus macroiterations in CASSCF algorithms (CI coefficients vs orbitals, see ref 80  The rate-determining step during the microiterations would normally be the evaluation of all the dispersion corrections for the entire GMTKN55 set, one after another, for a given combination of parameters. (When these were all done sequentially, the total wall clock time on our system was about 10−15 min for each microiteration, much of it due to the operating system overhead.) However, this step could be greatly accelerated by parceling out the individual D3BJ or D4 evaluations between all CPUs in a 40-core node. A minor I/O contention issue thus created was resolved by copying all required files onto a temporary RAM file system.
In the present case, the optimum value of the range separation parameter ω for a given fixed value of α is determined manually by interpolation. We repeated this process for six equally spaced α values, ranging from 0.57 to 0.72, to construct six different range-separated DSD (i.e., ωDSD) functionals.
From Figures 1 and S1 in the Supporting Information, it is clear that, for every c X,HF considered, the ωDSD n and ωDOD n -PBEP86-D3BJ variants benefited from range separation, whereas for the other exchange−correlation (XC) combinations, only c X,HF = 0.57 and to some extent c X,HF = 0.60 showed some advantage. Change of WTMAD2 (kcal/mol) for ωDSD X and ωDOD X -PBEP86-D3BJ with respect to range separation parameter ω (x-axis) for different c X,HF values. (Similar graphs for ωDSD X -PBEB95-D3BJ, ωDSD X -PBEPW91-D3BJ, and ωDSD X -PBEPBE-D3BJ and their ωDOD X versions can be found in Figure S1 in the Supporting Information.) Now, when we repeated the same experiment for other XC combinations, we found that none reached the accuracy of ωDSD n or ωDOD n -PBEP86-D3BJ in terms of WTMAD2 (see Table S3 in the Supporting Information). Therefore, we decided to proceed further with the range separation experiment only for PBEP86, where we considered the simultaneous variation of the short-range HF exchange (c X ) and range separation (ω) parametersusing full GMTKN55.
3.2. xDSD Functionals. In our previous study, 9 for the xrevDSD-PBEP86-D4 functional, we did not reoptimize the c X,HF parameter and instead took the earlier reported best value by Martin and coworkers 24 and optimized other linear parameters against GMTKN55. What if we also vary c X,HF ? (Note that a full set of GMTKN55 electronic structure calculations is necessary for each c X,HF value considered.) In the current study, we seek the c X,HF that minimizes WTMAD2.
We denote these new functionals xDSD n -PBEP86-Disp, where "n" stands for the percentage of HF exchange used. Following this notation, xDSD 69 -PBEP86-D4 is the same as the xrevDSD-PBEP86-D4 functional reported by us in ref 9.
xDSD n -PBEP86: We performed a full GMTKN55 evaluation for eight equally spaced c X,HF points, ranging from 0.57 to 0.78, followed by parameter optimization, taking both D3BJ and D4 dispersion corrections into account. Both with D3BJ and D4, c X,HF = 0.75 offered the lowest WTMAD2 instead of previously reported 24 c X,HF = 0.69 for xDSD-PBEP86which could be an artifact of optimizing c X,HF against a training set considerably smaller than that of GNTKN55.
With the D3BJ dispersion correction, the WTMAD2 for xDSD 75 -PBEP86-D3BJ is 2.144 kcal/molwhich is essentially identical to ωB97M(2) (WTMAD2 = 2.131 kcal/mol) but with fewer empirical parameters and (still) no range separation. (The latter is significant, considering that many codes are able to exploit density fitting in GHs but not rangeseparated hybrids.) Here, we should note that the ωB97M(2) functional was not trained against GMTKN55 but against a subset of the ca. 5000-point MGCDB84 (main group chemistry database 82 ); however, a fair amount of overlap exists between GMTKN55 and MGCDB84. Increased percentages of HF exchange in our optimized functional (i.e., going from c X,HF = 0.69 for xrevDSD-PBEP86-D3BJ to c X,HF = 0.75 for xDSD 75 -PBEP86-D3BJ) mainly benefited smallmolecule thermochemistry and intermolecular interactions (see Table S4 in the Supporting Information). Now, when we constrained c 2ss = 0 (i.e., the xDOD n -PBEP86-D3BJ functionals), c X,HF = 0.72 offered the lowest WTMAD2 (=2.231 kcal/mol); small-molecule thermochemistry suffered most of the deterioration resulted from applying this constraint.
Next, we repeated the same experiment for five different semilocal XC combinations, namely, PBEPBE, PBEB95, PBEPW91, SCAN, and BLYP (see Table 1 for WTMAD2 statistics and optimized parameters). Although, as expected, all xDSD n functionals surpassed the corresponding revDSD variants, none approached the accuracy of xDSD 75 -PBEP86-D3BJ; the only contender that came close was xDSD 77 -BLYP-D3BJ with 77% HF exchange. For the xDSD n -PBEPBE variants, one finds the "sweet spot" at c X,HF = 0.72, unlike the previously reported c X,HF = 0.68 in ref 24.
For the xDOD functionals (which permit us to use reducedscaling PT2 algorithms 31−36 as well as eliminate one empirical parameter), all except xDOD 69 -SCAN-D3BJ prefer a lesser percentage of exact exchange than the corresponding xDSD variants. The largest penalty for restricting c 2ss = 0 is paid by xDOD 74 -BLYP-D3BJthe WTMAD2 value drops from 2.254 kcal/mol to 2.564 kcal/mol. Upon further inspection, of all the 55 subsets, W4-11, TAUT15, and BSR36 are the three most affected ones. Similar to our previous observation for DSD-SCAN functionals, 9 xDSD 69 -SCAN-D3BJ sacrifices almost nothing when constraining c 2ss to be zero.
In our prior work, 9 for technical reasons, we adopted c ATM = s 6 , where c ATM is the prefactor for the Axilrod−Teller−Muto (ATM) 83,84 three-body correction term. If we allow c ATM as a variable, this somewhat reduces the WTMAD2 for xDSD n -PBEPBE-D4, xDSD n -PBEPW91-D4, and xDSDn-PBEB95-D4 and their xDOD variants. This leaves us with five adjustable dispersion parameters for the xDSD-D4 functionals. When we optimized all of them along with other parameters using BOBYQA, we noticed s 8 settling on values close to zero and c ATM on values close to one. Hence, if we constrain s 8 = 0 and c ATM = 1, the loss in accuracy is negligible, which permits eliminating two adjustable parameters. The only exceptions are Hence, going forward, we decided to freeze s 8 and c ATM throughout and optimize the remaining parameters (see Table  2). With D4 dispersion, the best performer of the xDSD family is again xDSD 75 -PBEP86-D4 (WTMAD2 = 2.119 kcal/mol), which now marginally outperforms Mardirossian and Head-Gordon's ωB97M (2)  (In response to a reviewer query, we have evaluated the impact of the recent revision 85 of D4, which corresponds to version 3 of the standalone dftd4 program, and found the difference for WTMAD2 to be a negligible 0.005 kcal/mol even for PBE0-D4, where s 6 = 1, unlike for the double hybrids.) Except for xDSD 75 -PBEP86-D4 and xDSD 77 -BLYP-D4, all other functionals prefer a very small fraction of the oppositespin MP2-like correlation. This is why for these functionals, we sacrifice very little by constraining it to zero, i.e., shifting from xDSD-D4 to xDOD-D4.
We should also mention the accidental similarity of xDSD 75 -PBEP86-D4 to Kaĺlay and coworkers' dRPA75 19 regarding the preferred percentage of exact exchange.
Can the performance of xDSDn-XC-D4 functionals be improved further by replacing the default three-body ATM term by the many-body MBD correction of Tkatchenko 86 and scaling it with a prefactor (now called c MBD rather than c ATM )? While we did find some improvement for one specific case, namely, xDSD 74 -PBEB95-D4MBD (the WTMAD2 drops from 2.403 to 2.288 kcal/mol), it appears that the molecules being considered here still are not large enough for higher-order MBD corrections to become statistically noteworthy and that our answer is hence inconclusive. (See Table S2 of the Supporting Information for some of our data.) In a recent study, Semidalas and Martin 3 have reported significant improvement for their composite methods by switching from the frozen core to the core−valence correlation and using the complete basis set (CBS) extrapolation from aug-ccpwCVTZ(-PP) and aug-cc-pwCVQZ(-PP) level calculations. Hence, we also checked whether further improvement of WTMAD2 statistics is possible by using a sufficiently large basis set and including the subvalence correlation in the MP2like part. Extrapolating from the core−valence aug-ccpwCVTZ(-PP) and aug-cc-pwCVQZ(-PP) energies for xDSD 75 -PBEP86-D3BJ using the L −3 formula for oppositespin and L −5 for same-spin MP2-like correlation, proposed by Halkier et al., 87 we found a change in WTMAD2 of up to three decimal places (0.00014 kcal/mol). We have therefore not explored this further for other double hybrids.
Similar to the previous section, here also, we checked how much performance we sacrificed by switching from ωDSD to ωDOD. With WTMAD = 2.204 kcal/mol, ωDOD 69 -PBEP86-D3BJ (ω = 0.10) appeared to be the best performer in this category. By and large, we gave up about 0.1 kcal/mol accuracy by the constraint c 2ss = 0; small-molecule thermochemistry is consistently the category the most affected by this restriction. Upon further inspection over all the 55 subsets, we found that BSR36 and TAUT15 are the main sources of this degradation. Both ωDSD 72 -PBEP86-D3BJ (ω = 0.13) and ωDOD 69 -PBEP86-D3BJ (ω = 0.10) only marginally outperform the corresponding ω = 0 variants xDSD 75 -PBEP86-D3BJ and xDOD 72 -PBEP86-D3BJ, respectively. Next, when instead of freezing a 2 at 5.5, we optimized it together with other parameters, we found almost no change in WTMAD2 statistics, neither for ωDSD nor for ωDOD.
Aiming for further improvement, we considered replacing the D3BJ term by D4 energy components. 30 Similar to what was mentioned earlier in this article, we found that imposing s 8 = 0 and c ATM = 1 caused only an insignificant increase of WTMAD2, although here optimal c ATM was somewhat further from unity. The lowest WTMAD2 we can get by shifting from D3BJ to D4 is 2.083 kcal/mol for ωDSD 72 -PBEP86-D4at the cost of eight adjustable parameters. The Journal of Physical Chemistry A pubs.acs.org/JPCA Article Now, we shift our focus from ωDSD-D4 to ωDOD-D4 functionals. The ωDOD 69 -PBEP86-D4 (ω = 0.10) functional is the best performer here with the WTMAD2 = 2.175 kcal/ mol. By having one parameter fewer (seven instead of eight), we sacrificed only 0.09 kcal/mol accuracy, and small-molecule thermochemistry is the reason behind this loss (see Table S4).
Similar to xDSD cases, here also, including the core−valence correlation for the MP2-like term or considering the MBD term beyond the three-body ATM correction did not help either.
We also considered eliminating the dispersion correction altogether. Similar to the global DH cases, 9 this approach significantly degrades the accuracy here too. The general trend shows the improvement in performance with increased HF exchange and the requirement of a higher ω value with respect to their ωDSD and ωDOD counterparts.
3.4. Benchmarks External to GMTKN55. The performances of the newly developed functionals were tested using four datasets, which are external to GMTKN55. These four test sets are MOBH35, originally proposed by Iron and Janes, 88 POLYPYR21, MPCONF196, 89 and CHAL336. 90 3.4.1. MOBH35. This database 88 comprises 35 reactions, including both early and late transition metal groups and 3d, 4d, and 5d transition metals. We extracted the best reported reference energies from the erratum 88 to the original 88  The Journal of Physical Chemistry A pubs.acs.org/JPCA Article MOBH35 article. The def2-QZVPP basis set was used with grids and auxiliary basis sets as described above in Section 2.2. Using a variety of multireference diagnostics, our group has recently found (E. Semidalas and J.M.L. Martin, unpublished) that reaction 9 exhibits severe static correlation in all three structures, which gets progressively worse from the reactant via the transition state to the product as the HOMO−LUMO gap narrows. Under these circumstances, as previously found for polypyrroles, 91 a large gap opens between canonical CCSD(T) and DLPNO-CCSD(T), and yet for the product, diagnostics are so large that one can legitimately question whether CCSD(T) itself is adequate for the problem. Hence, omitting this reaction from the MOBH35 dataset, we have recalculated MAD values for the remaining 34 reactions (see Figure 2).
In general, with the D3BJ dispersion correction, both rangeseparated and global DOD functionals perform better than their DSD counterparts. Shifting from D3BJ to D4 benefits the ωDSD (ωDOD) functionals across the board by 0.2−0.3 kcal/ mol. ωDSD 57 -PBEP86-D4 (ω = 0.22) achieves the lowest MAD of 1.0 kcal/mol, closely followed by the other range-separated DSD (DOD) functionals. Therefore, there is very little to choose among them.
Among the xDSD (xDOD) functionals with D3BJ, xDODs still do better than the xDSD variants. However, when we substitute the D4 dispersion correction, xDSDs are better performers than xDODs. The only exception is the PBEP86 XC combination, where xDOD 72 -PBEP86-D4 marginally outperforms xDSD 75 -PBEP86-D4 (see Figure 2). xDOD 72 -PBEP86-D4 offers the lowest MAD of 1.1 kcal/mol, close to what was found for ωDSD 57 -PBEP86-D4 (ω = 0.22). That being said, the other empirical range-separated double hybrids ωB2PLYP, ωB2GP-PLYP, and ωB97M(2) all have larger MAD values in the 1.6−1.7 kcal/mol range (see Figure 2). Finally, we can conclude that for the PBEP86 XC combination, shifting from the global to range-separated double hybrid is not so beneficial. Similar to what we found for the GMTKN55 test suite, considering higher-order MBD terms beyond the threebody ATM term has no perceptible benefit, though again, the systems under investigation may simply be too small. Now, the bimolecular reactions (i.e., reaction 17−20) could be problematic for a different reason, that is, because of  ) Therefore, if we also drop reactions 17−20 together with reaction 9 and recalculate MADs for the remaining 30 reactions, the MAD drops across the board. Also, while the MAD for ωB2PLYP and ωB2GP-PLYP remains elevated, that for ωB97M(2) now is in the same cohort as those of our best functionals (see Figure 3).

POLYPYR21
. This database 91 contains 21 unique structures of penta-, hexa-and heptaphyrins, which are [4n] πelectron expanded porphyrins that have generated considerable interest recently because of their potential application as molecular switches (see the introduction to ref 93 for a brief review). The structures are Huckel, Mobius, and figure-eight minima as well as the various transition states between them. Among them, the most troublesome are the Mobius rings, which exhibit a pronounced multireference character (for more details, see refs 91,93 ). CCSD(T)/CBS level reference energies were extracted from ref 93. We have used the def2-TZVP basis for all calculations here.
With the D3BJ dispersion correction, it appears that xDOD functionals perform noticeably better than their xDSD counterparts. In ref 93 for the problem at hand, as well as in the study by Iron and Janes 88 for MOBH35, the same trend was observed for DOD versus DSD functionals and ascribed to the greater resilience of spin-opposite-scaled GLPT2 to the static correlation. Now, if here we replace D3BJ by D4, the large difference between rmsd values of xDSD and xDOD functionals goes away (see Table 3). Finally, judging from the rmsd error statistics listed in Table 3, we observe that xDOD 74 -BLYP-D3BJ offers the lowest rmsd (1.64 kcal/mol) among the xDSD functionals.
In general, range-separated DSD double hybrids are better performers than the xDSD or xDOD variants (see Table 3). Switching from D3BJ to D4 dispersion correction deteriorates the performance for the ωDOD functional variants. Among all the xDSD (xDOD) and ωDSD (ωDOD) functionals tested, ωDOD 63 -PBEP86-D3BJ (ω = 0.16) offers the lowest rmsd = 0.45 kcal/mol, which, in fact, slightly outperforms the previously reported 93 top performer ωB97M(2) (0.63 kcal/ mol). However, in light of remaining uncertainties in the reference values, this difference should not be considered significant. Inspection of the Mobius structure data in isolation does reveal, across the board, that range-separated DHs cope with them much better than global double hybrids.  Together with all the new DHDFs we propose in the current study, we also report here the error statistics for our revDSD (revDOD) double hybrids, 9 Mardirossian and Head-Gordon's ωB97M(2), 53 ωB97M-V, 50 and ωB97X-V, 49 and Lars Goerigk and coworkers' ωB2GP-PLYP 58 and ωB2PLYP 58 functionals. All single-point energies were calculated using the def2-TZVPP 69 basis set throughout.
With the D3BJ dispersion correction, all new ωDSD and ωDOD functionals offer almost identical performance. Unlike what we saw for previous two external datasets, performancewise, there is practically no difference between revDSD and revDOD variants of a specific XC combination. Now, considering D4 instead of D3BJ benefits both rangeseparated and global DSD double hybrids across the board. As expected, small subsets are the least and large subsets are the most benefitted cases by this change. Specifically for large subsets, xDSD77-BLYP-D4 offers the lowest rmsd (0.49 kcal/ mol). However, for small subsets, the PBEB95 XC combination shows a particularly poor performance, even worse than that of lower-rung ωB97X-V. Between ωB97X-V and ωB97M-V, the former functional offers lower rmsd than the latter. Our new ωDSD (ωDOD)-PBEP86-D4 functionals are better performers than the combinatorially optimized, range-separated double hybrid, ωB97M (2).
For the PBEP86 XC combination, shifting from a global to a range-separated DSD-type functional does not offer any improvement in rmsd statistics. Finally, xDSD75-PBEP86-D4 is our best pick for the MPCONF196 database. However, we must acknowledge that most of the other global and rangeseparated DSD-D4 functionals are close contenders (see Table  4).
Similar to the previous two external datasets, and presumably for the same reasons, no benefit is seen from going beyond the three-body ATM term in D4.
3.4.4. CHAL336. While the present article was in peer review, Mehta et al. 90 proposed a comprehensive database of chalcogen bonding interactions, CHAL336. Consisting of molecules up to 49 atoms, CHAL336 contains 336 dimers, and a complete evaluation requires 1008 single-point energy calculations. These 336 dimer interaction energies can be subdivided into four categories: chalcogen−chalcogen, chalcogen−π, chalcogen−halogen, and chalcogen−nitrogen interactions. Mehta et al. have already assessed the performance of a large number of DFT methods (see Figure 9 of their article 90 ), and our revDSD-PBEP86-D3BJ 9 was among the best three performers (the differences between these are arguably a photo finish). Here, we evaluate the performance for CHAL336 of eight selected functionals, namely, xDSD 75 -PBEP86-D3BJ (D4), xDOD 72 -PBEP86-D3BJ (D4), ωDSD 72 -PBEP86-D3BJ (D4), and ωDOD 72 -PBEP86-D3BJ (D4). While we were at it, we also tested revDSD-PBEP86-D3BJ and revDSD-PBEP86-D4, the former to check the consistency with ref 90 and the latter as the D4 variant had not been included in ref 90. For the xDSD (xDOD) and ωDSD (ωDOD) functionals, we have used QCHEM 5.3, 67 whereas for the revDSD functional, we have used ORCA4 with the RIJCOSX (chain of spheres 100 ) approximation with the most accurate GRIDXS9. The same minimally augmented diffuse def2 basis set as in ref 90, ma-def2-QZVPP, 101 has been used across the board.
We found both xDSD75-PBEP86-D3BJ and xDOD72-PBEP86-D3BJ to perform slightly better than revDSD-PBEP86-D3BJ (see Figure 4 and Table S17 in the Supporting Information). However, for the subsets of systems where both canonical and DLPNO reference values were available in ref   Tables 1 and S4 in ref 90). Hence, the apparent improvement in the statistics of xDSD (xDOD) and ωDSD (ωDOD) functionals arguably does not rise above numerical noise.
From the ma-def2-QZVPP and ma-def2-TZVPP energies, two-point CBS extrapolation has been performed for each of the above functionals using the L −3 formula (where L is the cardinal number of the basis set), which works out in practice 7297. CBS extrapolation significantly improved the performance for both xDSD75-PBEP86-D4 (the MAD improves from 0.46 to 0.28 kcal/mol) and xDOD72-PBEP86-D4 (the MAD drops from 0.47 to 0.29 kcal/mol). Suffice to say that all functionals considered here are at least competitive with, and perhaps superior to, the best performers in the CHAL336 article, despite this dataset not having been involved in parameterization.

CONCLUSIONS
Aiming to improve our previous revDSD family functionals further, we have considered both range separation ωDSDs and xDSDsthe latter, while analogues of the XYG3 family of functionals, are also recovered as the ω = 0 limit of rangeseparated double hybrids. Concerning our first research objective: from an extensive survey, we can conclude the following: (a) xDSD-D3BJ functionals have a slight advantage over our prior revDSD family functionals, which can be further improved upon by replacing D3BJ with D4. (b) For D4, allowing c ATM , the prefactor for the three-body ATM term, to take on values different from one does not reduce WTMAD2 by an amount statistically significant enough that it would justify the introduction of the extra adjustable parameter. Replacing the ATM term by the many-body dispersion model of Tkatchenko and coworkers achieves no significant benefit, although the systems in GMTKN55 may simply not be large enough to rule this out. (c) For the xDSD n -PBEP86-D4 variants, c X,HF = 0.75 offers the lowest WTMAD2, unlike the previously reported c X,HF = 0.69 for DSD-PBEP86-D4. However, when we imposed the c 2ss = 0 constraint, WTMAD2 reaches a minimum at c X,HF = 0.72. (d) In terms of WTMAD2, xDSD 75 -PBEP86-D4 marginally outperforms ωB97M(2), 53 hitherto the "record holder" for the lowest WTMAD2, 8 but without its range separation and using just a half-dozen empirical parameters. In view of the uncertainty in the reference values, however, and the fact that xDSD 75 -PBEP86-D4 was trained on GMTKN55 itself rather than on a different albeit strongly overlapping set such as ωB97M (2), it is probably safer to say that the two functionals are competitive. Concerning the second research objective, applying range separation over the HF exchange part, we found the lowest WTMAD2 for c X,HF = 0.72 and ω = 0.13. With D3BJ, WTMAD2 is 2.108 kcal/mol, which can be lowered slightly further by substituting D4 (2.083 kcal/mol). Therefore, range separation helped us to improve the performance slightly beyond that of xDSD 75 -PBEP86-D3BJ(D4) and in turn a little further beyond that of ωB97M(2)again using just a halfdozen adjustable parameters. Although ωB97M(2) outperforms all the new ωDSD and ωDSD functionals for smallmolecule thermochemistry, this is outweighed in WTMAD2 by the superior performance of the new functionals for conformer equilibria.
All in all, however, the improvement for GMTKN55 from introducing range separation in DSD functionals is quite modest, somewhat surprisingly so. In some sense, this is convenient as GHs tend to be computationally more economical.
For some perspective beyond the comparison of WTMAD2 [where differences between, for example, ωB97M(2) and ωDSD 72 -PBEP86-D4 (ω = 0.13) may well be comparable to the residual uncertainty in the reference values], let us consider the performance of four representative functionals, namely, ωB97M(2), revDSD-PBEP86-D4, xDSD 75 -PBEP86-D4, and ωDSD 72 -PBEP86-D4 (ω = 0.13), for four external test sets not involved in the parameterization process. Two of these, metal− organic barrier heights (MOBH35) 88 and especially the isomer equilibria and interconversion barriers in polypyrroles (POLYPYR21), 91,93 put the functionals' performance in the presence of static correlation to the test. The two others are CHAL336, a very recently published 90 benchmark of chalcogen bonding interactions, and MPCONF196, 89 which features conformational energies of smaller peptides and of medium-sized macrocycles.
All four options perform very well for MOBH35 if pathologically multireference reaction 9 and bimolecular reactions 17−20 are removed (see above). For CHAL336, again, all four options perform excellently if the basis set extrapolation is performed to quench the effect of the BSSE; two of the options have smaller MAD values by about 0.05 kcal/mol, but in light of the residual uncertainty of about 0.15 kcal/mol rmsd in the reference data, this difference may be deemed insignificant.
For MPCONF196, ωB97M(2) still performs very well but less so than the other three options. This is consistent, actually, with the breakdown of WTMAD2 for GMTKN55 into five top-level subcategories: ωB97M(2) has an edge there for the small-molecule thermochemistry component, which is compensated by the better performance for the intramolecular interaction component.
This leaves POLYPYR21, where the two range-separated options are clearly superior and pronouncedly so for the Mobius structures (which have very pronounced static correlation 91,93 ). It may therefore be that range-separated double hybrids, be they ωB97M(2) or ωDSD, have an edge for these kinds of problems; moreover, as will be seen in the companion article (part II), the combination of range separation with post-PT2 corrections turns out to be more generally advantageous. 20 ■ ASSOCIATED CONTENT * sı Supporting Information breakdown of total WTMAD2 into five top-level subsets for all the xDSD and range-separated DSD functionals; MAD, MSD, and rmsd values as well as breakdown of total WTMAD2 by each subset for new functionals; MAD statistics of various DSD-DHs for the full CHAL336 and its four subsets; and QCHEM sample inputs for ωDSD (ωDOD)-D3BJ and ORCA sample inputs for xDSD (xDOD)-D3BJ functionals (PDF)