Do Double Hybrid Functionals Benefit from Regularization in the PT2 term? Observations from an Extensive Benchmark
that MP2 regularization might improve the performance of double-hybrid density functionals. Using the very large and chemically diverse GMTKN55 benchmark, we find that κ-regularization is indeed beneficial at lower percentages of Hartree-Fock exchange, especially if spin-component scaling is not applied (such as in B2GP-PLYP or ωB97M( 2)). This benefit dwindles for DSD and DOD functionals, and vanishes entirely in the ca. 70% HF exchange region optimal for them.
Introduction:<!>Computational Details:<!>Results and Discussion:<!>A B C D<!>Conclusions:
<p>Double-hybrid density functional (DHDF) theory (for reviews see [1][2][3][4] ) represents a special case of fifth-rung functionals on 'Jacob's Ladder' 5 (the fifth rung is where dependence on unoccupied orbitals enters). As such, DHDF theory resides on the seamline between DFT and wavefunction approaches. DHDFs are the most accurate DFT methods available to date for main group energetics, [6][7][8] transition metal catalysis, [9][10][11][12] electronic excitation spectroscopy, [13][14][15][16][17][18][19][20][21] external magnetic [22][23][24][25] and electric field 26,27 induced properties, and vibrational frequencies. 28 Its reliance on a second-order perturbation theory term for part of the correlation energy creates an Achilles' Heel for molecules with small band gaps, due to the presence of orbital energy differences in the denominator. Using the spin-orbital notation, the MP2-like Görling-Levy 29 nonlocal correlation term has the following form:</p><p>where i and j are the occupied orbital, a and b are virtual orbital indices, and the energy denominator</p><p>One remedy that has been proposed in the past, in the context of single-or multireference MP2, has been DCPT (degeneracy-corrected perturbation theory 30 ). Another is regularization of the above expression to remove the singularity for ∆ $% &' → 0. Stück and Head-Gordon 31 proposed a simple level-shift regularizer (δ) in the context of ROOMP2 [32][33][34] (restricted orbital-optimized second-order perturbation theory). They found that this works well for single-bond breaking, but the required level shifts for multiple bonds were so large that they disrupted thermochemistry results. 35 Lee and Head-Gordon 36 have recently proposed two energy-gap-dependent regularizers: σ and κ. Although these forms were developed initially with OOMP2 in mind, Head-Gordon and coworkers 37,38 have shown that even without orbital optimization, σ-and κ-MP2 can achieve better accuracy than ordinary (unregularized) MP2 for both main-group and transition-metal thermochemistry, barrier heights, and non-covalent interactions. They also found that both sets of regularizers performed comparably: hence we will focus here only on κ-regularized MP2 correlation, which has the following expression:</p><p>In which κ is a fixed regularization parameter. In the large-κ or large-∆ $% &' limits, the regularization factor approaches unity, while in the small-κ or small-∆ $% &' limits, the corresponding term in the energy summation approaches zero, i.e., the Hartree-Fock energy is recovered. Ref. 37 states that "We are optimistic that this study could pave the way for future development of double-hybrid density functionals based on nonlocal correlation expressions that are more appropriate than conventional MP2 for large dispersion-bound systems and organometallic bonding, yet still free of self-correlation errors. κMP2, σ 2 MP2, and σMP2 are promising candidates in this regard."</p><p>While such DHDFs as our own minimally empirical revDSD, 37 revDOD, 37 and xDSD 39 families or the more heavily parametrized ωB97M(2) range-separated DHDF 40 reach accuracies approaching those of composite wavefunction approaches (see Ref. 41 for a head-to-head comparison), there might still be room for further improvement, particularly in terms of resilience for systems with small band gaps, and hence significant type A static correlation (a.k.a., absolute near-degeneracy correlation 42 ).</p><p>In the present work, we will attempt to confirm or refute the above conjecture from Ref. 37 That is: to find out whether or not using spin-component-scaled κ-GLPT2 instead of unregularized (i.e., conventional) PT2 correlation can further improve the performance of DSD functionals. For the κ regularized DSD functionals, the final energy has the following expression: where EN1e stands for the sum of nuclear repulsion and one-electron energy terms; EX,HF represents the exact exchange and cX,HF the corresponding coefficient; EX,DFT is the exact exchange energy component from the semilocal generalized gradient approximation (GGA) and cX,DFT is the corresponding parameter. EC,DFT represents the semilocal GGA correlation component, and cC,DFT is the coefficient for that energy part. E2ab,κ-PT2 and E2ss,κ-PT2 are the opposite-spin and same-spin κ-GLPT2 correlation energies, and their respective linear coefficients are c2ab and c2ss. Finally, Edisp is a dispersion correction such as D3(BJ) [43][44][45] with any associated adjustable parameters: in the present work (as in Ref. 46 ), the nonlinear damping-function shape parameters are fixed at a1=0, a2=5.5, and the higher-order coefficient is fixed at s8=0 (as we have found 46 to be appropriate for double hybrids).</p><!><p>We have used the GMTKN55 benchmark suite 7 (General Main-Group Thermochemistry, Kinetics, and Noncovalent Interactions) throughout. It comprises 55 types of chemical problems, which can be further divided into five major subcategories: thermochemistry of small and medium-sized molecules, barrier heights, large-molecule reactions, intermolecular interactions, and intramolecular interactions (or conformer energies). WTMAD2 (weighted mean absolute error) as defined in Eq. (2) of Ref. 7 has been used as our primary metric:</p><p>where |∆𝐸| FFFFFF . is the mean absolute value of all the reference energies from 𝑖 = 1 to 55, 𝑁 . is the number of systems in each subset, 𝑀𝐴𝐷 . is the mean absolute difference between calculated and reference energies for each of the 55 subsets. For the details of all 55 subsets with proper references, see Table 1 in ref. 7 All electronic structure calculations were performed using the QCHEM 5.4 47 package on the ChemFarm HPC cluster in the Faculty of Chemistry at the Weizmann Institute of Science. The Weigend−Ahlrichs def2-QZVPP 48 basis set was used for all subsets except seven, where we used the diffuse-function augmented variant def2-QZVPPD 49 instead: the rare gas clusters RG18 and the six anion-containing subsets WATER27, IL16, G21EA, BH76, BH76RC, and AHB21. For the C60ISO and UPU23 subsets (which have comparatively small weights in WTMAD2), we settled for the def2-TZVPP basis set to reduce computational cost. For the remaining 46 subsets, the def2-QZVPP 48 basis set was used. The corresponding standard RI-MP2 50 and Coulomb-exchange fitting (RI-JK) 51 basis sets were employed throughout to reduce the computational cost further. The large pruned integration grid SG-3 was used across the board. 52 Generally, inner-shell orbitals were frozen, but in subsets where such orbitals come close enough to the valence shell to qualify as 'honorary valence orbitals', 53 the same frozen-core settings were used as in Refs. 46,54,55 Similar to the unregularized (x)DSD functionals, one fully optimized κ-DSD also has six adjustable linear parameters: cX,HF, cC,DFT, c2ab, c2ss(=c2aa+c2bb), and for the D3(BJ) dispersion correction, one prefactor s6 and one parameter a2 for the damping function (like in refs 56,57 we constrain a1=s8=0). On top of that, our new κ(x)DSD functionals also contain the PT2 regularization parameter κ.</p><p>Powell's BOBYQA 58 (Bound Optimization BY Quadratic Approximation) derivative-free constrained optimizer together with inhouse written scripts and Fortran programs were used to optimize all parameters.</p><!><p>In order to find our bearings, we first explore the effect of using κ-GLPT2 energies in the revDSD-PBEP86-D3BJ and xDSD75-PBEP86-D3BJ double hybrids using the "Diet-100" statistical reduction 59 version of GMTKN55. The unregularized same-spin and opposite-spin PT2 energies were replaced by the corresponding κ-PT2 terms, evaluated for a fixed κ value, then all four linear parameters were reoptimized to obtain the regularized double hybrids, κDSD-PBEP86-D3BJ and κxDSD75-PBEP86-D3BJ. Including κ = ∞ (i.e., unregularized GLPT2), we calculated WTMAD2 for seventeen κ values ranging from 0.9 to 10. For both κDSD-PBEP86-D3BJ and κxDSD75-PBEP86-D3BJ, WTMAD2 decreases with increasing κ, and beyond κ = 4.0, WTMAD2 approaches WTMAD2κ=∞ (see Figure S1 in the Supporting Information). Interestingly, for κDSD-PBEP86-D3BJ, we observe a very shallow dip near κ = 4.0, but WTMAD2κ=4.0 is only 0.03 kcal/mol lower than WTMAD2κ=∞.</p><p>Having initially scanned κ-space for one gDH and one xDH, we then used full GMTKN55 for the final parametrization of the κ-regularized DSD functionals. For this purpose, we selected nine κ values (i.e., κ = 1.1, 1.45, 2.0, 2.2, 2.5, 3.0, 3.5, 4.0, and ∞) and three different exchangecorrelation (XC) combinations: PBE-P86, PBE-PBE, and B88-LYP. Among the three regularized DSD (i.e., κDSD-PBEP86-D3BJ, κDSD-PBEPBE-D3BJ, and κDSD-BLYP-D3BJ) and one regularized xDSD (i.e., κxDSD75-PBBEP86-D3BJ) functionals, only κDSD-BLYP-D3BJ marginally benefits from the PT2 regularization. We obtained the lowest WTMAD2 (2.34 kcal/mol) for κDSD-BLYP-D3BJ at κ = 3.0, which is just 0.09 kcal/mol lower than the WTMAD2 of the unregularized counterpart, revDSD-BLYP-D3BJ, see Figure 1 and Table S1 in the Supporting Information. (By way of perspective, in the Supporting Information of Ref. 55 we applied the Bayesian Information Criterion 60,61 to see what reduction in WTMAD2 could be considered "decisive", and found 3.5%, or (in this context) about 0.09 kcal/mol.) Now, partitioning each WTMAD2 into five major subcategories, we found that small-molecule thermochemistry, barrier heights, and intramolecular interactions do not benefit at all from PT2 regularization and κ = ∞ always offers the best performance (see Table S1 in the Supporting Information). However, for the large-molecule reactions, all four κDSD functionals benefit from regularization, most prominently for κDSD-BLYP-D3BJ. Finally, for the intermolecular interactions, κ-regularization improved performance for κxDSD75-PBEP86-D3BJ, κDSD-PBEP86-D3BJ, and κDSD-PBEPBE-D3BJ. We obtained the best results at κ=2.5 while using κxDSD75-PBEP86-D3BJ, and at κ=1.45 for κDSD-PBEP86-D3BJ and κDSD-PBEPBE-D3BJ. However, the performance of those three functionals deteriorates with increasing κ value. The unregularized variant wins the race for this subset when κDSD-BLYP-D3BJ is considered (see Figure S2 and Table S1 in the Supporting Information).</p><!><p>What if we only allow one degree of freedom between c2ab and c2ss but self-consistently reoptimize the remaining three parameters (i.e., cC,DFT, c2ab, and s6)? We considered two possibilities: first, imposing the condition c2ab=c2ss, like in simple double hybrids like B2PLYP, 62 B2GPPLYP, 63 PBE0-2, 64 etc. as well as the combinatorially optimized ωB97M(2). 40 The nomenclature, 'DSD', is no longer appropriate in this case and will be replaced by κDH-XC-D3BJ in the present manuscript. The second option is to exclude same-spin correlation entirely (c2ss=0), a.k.a., the DOD forms. These are of interest because they are amenable to reduced-scaling opposite-spin MP2 techniques like Laplace transform MP2 of Häser and Almlöf, 65 which scales as Ο(N 4 ) with system size, or the tensor hypercontraction approach of Song and Martinez, 66 which scales as Ο(N 3 ). As expected, constraining c2ab=c2ss always increases WTMAD2 over the κDSD forms. Both for κDH-PBEPBE-D3BJ and κDH-PBEP86-D3BJ, we found the lowest WTMAD2 at κmin=2.5, corresponding to 0.11 and 0.10 kcal/mol reductions, respectively, compared to their unregularized forms. For κDH-BLYP-D3BJ, the WTMAD2 gap between κmin (κ value for which we obtain the minimum WTMAD2) and κ=∞ is ~0.08 kcal/mol. Interestingly, for the κDOD functionals, GLPT2 regularization does more harm than good, and the unregularized forms always offer the best performance (see Figure 1).</p><p>The use of the more modern D4 67,68 dispersion correction instead of D3BJ for the κ(x)DSD functionals does not affect any trends with respect to the regularizer. For κDSD-PBEPBE and κDSD-PBEP86, the WTMAD2 gaps between the D3BJ and D4 corrected forms increase gradually with rising κ value. However, switching from D3BJ to D4 has no significant effect on the performance of the κxDSD75-PBEP86-D4 and κDSD-BLYP-D4 functionals (see Table S2 in the Supporting Information). The same story is repeated for the κDOD and κxDOD functionals with D4 dispersion correction (see Table S3 in the Supporting Information).</p><p>Thus far, for the regularized DSD functionals, we have used the parameter for exact exchange same as for their unregularized forms (i.e., cX,HF = 0.69, 0.68, 0.71, and 0.75 for revDSD-PBEP86-D3BJ, 46 revDSD-PBE-D3BJ, 46 revDSD-BLYP-D3BJ, 46 and xDSD75-PBEP86-D3BJ, 39 respectively). Our next objective is to check whether or not the cX,HF Hartree-Fock exchange prefactors of the unregularized forms are still optimal for the new regularized functionals? To answer this question, we have considered seven cX,HF ranging from 0.5 to 1.0, and the κDSDX-PBEP86-D3BJ functional together with its three other variants: κDHX-PBEP86-D3BJ, κDODX-PBEP86-D3BJ, and κnoDispSDX-PBEP86 (X represents the percentage of HF exchange used. Only one joint set of electronic structure calculations is required for all variants).</p><p>It turns out that regularization in κDSD-PBEP86-D3BJ becomes gradually more beneficial as cX,HF is decreased, with κmin decreasing concomitantly. For example, we obtain the lowest WTMAD2 near κ=1.67 when cX,HF=0.55, but when cX,HF=0.5 is used, the κmin drops to 1.45 (see Figure 2). Now, splitting each WTMAD2 into five major subcategories, we found that regularization does more harm than good across the board for small-molecule thermochemistry. However, the extent of performance deterioration with respect to the κ values becomes less prominent as cX,HF in κDSD-PBEP86-D3BJ is reduced. Using a small κ value can severely harm the performance for barrier heights at elevated percentages of HF exchange, but much less so at lower percentages. At κ=1.1, the κDSD50-PBEP86-D3BJ marginally outperforms the unregularized variant. For large-molecule reactions, κDSD-PBEP86-D3BJ is a better choice across the board compared to the revDSD-PBEP86-DBJ functional, and κmin decreases gradually with the increase of cX,HF. The κDSD69-PBEP86-D3BJ and κDSD66-PBEP86-D3BJ functionals offer the best performance near κ=3.0 and 2.0, respectively. With κ=1.45, the ΔWTMAD2 (i.e., the contribution from a particular subset for total WTMAD2) for κDSD55-PBEP86-D3BJ is very close to what we obtained for the unregularized functional form. For intramolecular interactions, the trends are largely the same as what we obtained for the small-molecule thermochemistry subsets. Finally, for intermolecular interactions, regularized forms always outperform the unregularized alternatives, and κmin decreases with increasing of cX,HF. For this subset κDSD69-PBEP86-D3BJ and κDSD63-PBEP86-D3BJ are the two best picks at κ=1.45 and 1.33, respectively. (For optimized parameters, total WTMAD2 for full GMTKN55, and its decomposition into five major subsets see Table S4 in the Supporting Information.)</p><p>As we found no material change from D3BJ to D4 in the previous section, we have decided not to explore that avenue for different cX,HF values.</p><p>Now, what happens if we impose c2ab=c2ss , i.e., a simple double hybrid rather than a DSD or DOD form? (As always, parameters are re-optimized self-consistently.) The results can be found in Figure S3 and Table S5 in the Supporting Information. Even at cX,HF=0.75 we find a shallow WTMAD2 reduction (0.04 kcal/mol) at κmin; as cX,HF is decreased, this 'well' is deepened until it reaches 0.36 kcal/mol at cX,HF=0.50. (κmin decreases in tandem with decreasing cX,HF.) Among the five major subsets, at low cX,HF the WTMAD2 component from noncovalent interactions (NCI) decreases as κmin decreases, while this is detrimental to small-molecule thermochemistry and (at κmin =1.1 to 1.45) for barrier heights as well: the former tendency grows weaker, and the latter stronger, as cX,HF is increased. We note that the ωB97M(2) combinatorially optimized range separated double hybrid of Mardirossian and Head-Gordon 40 has c2ab=c2ss and might hence benefit. (The way spin-component-scaled MP2 behaves differently from standard MP2 has been rationalized to some degree as approximate higher-order effects. 69,70 ) Interestingly, the behavior seen for the κDODX-PBEP86-D3BJ functionals (i.e., where c2ss=0) is fairly similar, with lower percentages of HF exchange significant benefiting from PT2 regularization. For example, the WTMAD2 'well' for κDOD50-PBEP86-D3BJ is 0.39 kcal/mol, but for κDOD63-PBEP86-D3BJ, this shrinks to just 0.04 kcal/mol. Except for small molecule thermochemistry, we found the same trend as κDSDX-PBEP86-D3BJ for the remaining subsets. At small cX,HF, regularization seems slightly beneficial for the small molecule thermochemistry subsets (see Figure S3 and Table S6 in the Supporting Information).</p><p>For the sake of completeness: for dispersion-uncorrected κnoDispSDX-PBEP86 functionals, we likewise found that at elevated fractions of HF exchange, κmin approaches infinity (see Figure S3 and Table S7 in the Supporting Information). In this case, obviously the balance is tipped more strongly toward high percentages of PT2 correlation (and of HF exchange, as the two are well known 63 to be linearly related) as long-range dispersion is not covered by anything else.  Thus far, we have used the same percentage of HF and semilocal exchange for orbital generation and final energy calculation. What happens when we use PBE0P86 (i.e., 0.25HFx + 0.75PBEx + 1.0P86c) orbitals in the final energy evaluation of κDSD69-PBEP86-D3BJ? These new DSD functionals are found not to benefit from the PT2 regularization (see Table S8 and Figure S4 in the Supporting Information). With WTMAD2=2.28 kcal/mol, the unregularized version, revDSD-PBEP86-D3BJ@PBE0P86, slightly outperforms the original revDSD-PBEP86-D3BJ (WTMAD2=2.38 kcal/mol): nearly all of that gain comes from RSE43 (radical-stabilization energies), due to substantially reduced spin contamination thanks to the smaller cX,HF in the orbitals.. If we constrain c2ab=c2ss, the WTMAD2 gap between κmin and κ=∞ is 0.13 kcal/molmarginally larger than what we obtained for κDH69-PBEP86-D3BJ (0.10 kcal/mol). κmin also increases from 2.0 (for κDH69-PBEP86-D3BJ) to 3.0. For the κDOD variants of these new functionals, we saw the same trend as the κDSD forms.</p><!><p>From an extensive survey of regularized DHDFs using the large and chemically diverse GMTKN55, we can conclude the following:</p><p>• The benefits of PT2 regularization for intermolecular interactions and large molecule reactions are negated by the losses for small molecule thermochemistry, barrier heights, and conformer energies. Hence, κ-GLPT2 correlation causes no significant reduction in WTMAD2 compared to the unregularized revDSD-PBEP86-D3BJ, revDSD-PBE-D3BJ, and xDSD75-PBEP86-D3BJ functionals. However, the significantly better performance of κDSD-BLYP-D3BJ at κ=2.2 for large-molecule reactions has enough impact on WTMAD2 that overall, it marginally outperforms revDSD-BLYP-D3BJ. Replacing D3BJ by D4 dispersion does not affect those trends. • If we eliminate spin-component scaling (i.e., c2ab=c2ss), the WTMAD2 gap between the κmin and κ=∞ (i.e., unregularized) forms of κ(x)DH-XC-D3BJ are more significant than that we obtained for the κ(x)DSD functionals. In contrast, for the κ(x)DOD forms, unregularized functionals always perform better. • Regularization of the GLPT2 terms in double-hybrids is most helpful at lower percentages (e.g., 50%) of HF exchange. At higher percentages of HF exchange, the benefits for intermolecular interactions and large-molecule reactions are outweighed by the deterioration in the remaining three subsets. At lower percentages of HF exchange, the benefits are heightened and the deterioration is mitigated, hence an overall beneficial effect.</p><p>In special cases where DHs with a small fraction of HF exchange might be more resilient (e.g., systems with strong static correlation or prone to severe spin contamination error), κregularized double hybrids will offer advantages over their unregularized counterparts.</p>