Improving robustness of fifth-rung functionals with linearized ladder correlation

Abstract

Strong electron correlation effects due to multiconfigurational character remain a consistent challenge for modern density functional theory (DFT). Despite their prominent position atop Jacob's Ladder at rung 5, double-hybrid functionals are often even more sensitive to strong correlation than lower-rung functionals, casting doubt on their general usefulness for transition-metal complexes. Much of these problems can be traced to the second-order Møller–Plesset (MP2) correlation energy, which is the most common wavefunction-theoretic ingredient in double-hybrid functionals. Herein we introduce a novel, size-consistent double-hybrid functional that uses Coulomb-attenuated linearized hole–hole ladder coupled-cluster correlation as a robust alternative to MP2. The resultant ωB97X-L-V functional performs as well as other traditional double-hybrids on main-group chemistry, can smoothly dissociate covalent bonds, and performs well for reaction energies and bond dissociation energies of transition-metal complexes. Our work sets the stage for further improvements to ωB97X-L-V and future double-hybrid functionals that add many-body screening to the MP2 correlation energy.

---

Introduction

Density functional theory (DFT) remains the workhorse of modern quantum chemistry. While exact, in principle, practical DFT calculations invoke a wide variety of approximations for the exchange–correlation (XC) energy. These approximations are often grouped into classes, or “rungs” on the Jacob's Ladder¹ of DFT. Each rung is defined by the information about the one-electron density (or Kohn–Sham orbitals in generalized Kohn–Sham theory) that is used by the density functional in calculating the energy. By far, the most commonly used DFT approximations belong to rungs 2, 3, and 4, which constitute generalized gradient approximations (GGAs), meta-GGAs, and hybrid DFT, respectively. Sitting atop Jacob's Ladder at rung 5 are double-hybrid (DH) functionals,^2–4 which incorporate information about the density, its derivative(s), occupied Kohn–Sham orbitals via exact exchange, and virtual Kohn–Sham orbitals via wave function theory (WFT) correlation.

Inclusion of WFT correlation makes DH functionals more expensive but generally improves their accuracy relative to lower-rung approximations. The WFT of choice is often second-order Møller–Plesset perturbation theory (MP2),^5,6 as it is the most affordable option at [scr O, script letter O](N⁵), where N is the number of atomic orbital basis functions. However, the higher computational cost of DH functionals does not always ensure better accuracy, especially in systems where the underlying MP2 correlation energy breaks down. This was recently exemplified in a cautionary tale by Shee and co-workers regarding the significant errors encountered when applying DH functionals to the bond dissociation energies of transition-metal complexes.⁷ Despite these results, there are DH functionals that manage to describe transition-metal complexes well,^8–11 but many of the best performers use highly parameterized formulations of MP2 that empirically scale same-spin and opposite-spin correlation, sometimes even neglecting the former altogether.^12–16 This requisite empiricism suggests a deficiency in the MP2 correlation energy.

Another example of DH functionals under-performing for transition metal complexes comes from an appraisal of 240 density functionals on spin-state energetics and binding energies of Fe porphyrins in the POR21¹⁷ data set,¹⁸ where DH functionals exhibited the second-highest mean absolute deviations (MADs) of any rung on Jacob's Ladder, with only rung 1 (local density approximation) functionals performing worse. DFT (and DHs in particular) may perform poorly on POR21 because 3d metal complexes can be somewhat strongly correlated, and strong correlation remains a challenge for modern functionals. In lanthanides, for example, the multireference character of partially filled f-orbitals, high number of core electrons, and significant relativistic effects conspire to form a complex electronic structure that poses an exceptional challenge to current DFT methods.^19–22

While modern rung 1–4 functionals struggle with static correlation, MP2-based rung 5 functionals can perform even worse in comparison, defying the expectation set by their prominent position on Jacob's Ladder.^7,18 However, for manifestly single-reference transition-metal complexes like those in the MOR41⁹ and ROST61¹¹ data sets of closed- and open-shell metal–organic reaction energies, DH functionals do generally offer improvements. There is also evidence for DH functionals performing very well on MOBH28,²³ a refined subset of the original MOBH35¹⁰ set of metal–organic reaction barrier heights that removed several strongly-correlated systems. For these sets, many DH functionals even perform within 3d-transition-metal chemical accuracy,^9–11 defined as ±3 kcal mol⁻¹.²⁴ Even so, there remains a great deal of room to improve the performance of DFT on transition metal complexes that are not so carefully curated to exclude systems of (even moderate) multireference character.

Increasing the robustness of DFT approximations for strongly-correlated systems has been a hotbed of research.^25–32 Becke's foundational B05²⁵ and B13²⁶ functionals for strong correlation have inspired much of the work in the field.²⁷ For instance, the ratios between exchange energy densities used by local hybrid functionals that modulate the fraction of exact exchange in regions of space where strong correlation is pronounced have their foundation in KP16/B13.³⁰ Much of the work in this domain is well characterized by the concept of the “zero-sum game” in DFT, as introduced by Kaupp and co-workers.³⁰ This concept points out the apparent paradox wherein the addition of exact exchange to the XC functional reduces self-interaction error but evidently increases statistical errors when modeling strongly-correlated systems.

Whereas the aforementioned approaches attempt to modify the underlying density functional, there have been parallel efforts to move beyond MP2 correlation for the WFT component of DHs. These efforts span truly multireference approaches such as multiconfigurational pair DFT,^33–35 to interfaces with single-reference theories that are more adept at describing strongly correlated systems, like seniority-based coupled-cluster.^36–38 Of course, the random-phase approximation (RPA),^39–41 often in its direct (exchangeless) form has been explored as an alternative to MP2 correlation in the context of DFT.^42–51 Other work has explored the application of regularized MP2 methods^52–55 in place of canonical MP2,⁵⁶ finding that regularized MP2 does not improve the numerical accuracy of DHs for well-behaved main-group chemistry. There have also been important developments in the interface of DFT with particle–particle RPA (pp-RPA)^57,58 correlation through pairing matrix fluctuation.^59,60 Functionals paired with pp-RPA have been shown to be quite useful for modeling strongly-correlated systems and excited states of point defects.^61,62

Inspired by the latter, we construct a DH functional that attempts to get the most out of a single-reference WFT by using the hole–hole approximation to the linearized ladder coupled cluster doubles (LinLCCD(hh)) ansatz that one of us recently introduced.⁶³ This alternative to MP2 is a form of linearized pp-RPA that avoids catastrophic divergences in strongly-correlated systems while its linear form may be leveraged to reduce computational cost. To interface LinLCCD(hh) with DFT XC, we use adiabatic connection to construct a DH functional with only a single parameter connecting DFT to WFT correlation. We then re-optimize the parameters from the ωB97X functional for best performance with LinLCCD(hh).

We also consider a unique ansatz where we confine the LinLCCD(hh) contribution to the short-range part of the range-separated Coulomb potential, as we posit that the challenges posed to DFT in describing transition-metal complexes manifest locally near the metal center. This also permits us to reduce the cost of our functional even more in the future by means of aggressive integral screening akin to attenuated MP2.^64,65 Our results suggest that such a LinLCCD(hh)-based DH functional can perform about as well as other DH functionals on main-group chemistry. Our functional, which we call ωB97X-L-V, is also more robust than MP2-based DH functionals in cases where static correlation is pronounced, such as covalent-bond breaking.

Theory

Our functional follows the general design principles laid out by Kalai and Toulouse.⁶⁶ However, we interface our functional with a linearized coupled cluster wave function, and this choice imposes some subtleties that heretofore have been left undiscussed in the literature.

We begin with the universal density functional,⁶⁷

where T̂ is the kinetic energy operator, V̂ is the Coulomb operator, and the minimization is carried out over antisymmetric multideterminant wave functions Ψ to give fixed density ρ. We choose to range separate the Coulomb operator into a short-range (sr) and long-range (lr) component to give,where λ is an interaction strength parameter and E^sr,ω,λ_Hxc [ρ] is the (short-range) Hartree-exchange–correlation (Hxc) density functional contribution. This Hxc contribution can also be decomposed further as,

E^sr,ω,λ_Hxc [ρ] = E^sr,ω,λ_H [ρ] + E^sr,ω,λ_xc [ρ] (3)

The separation into long-range, V̂^lr,ω, and short-range V̂^sr,ω is not unique. In this work, we choose the error function approach to partition the Coulomb operator,

where the second and third terms correspond to long-range and short-range contributions, respectively. Note that due to its linearity in terms of the density, the short-range Hartree contribution from eqn (3) can be straightforwardly written as,

With the rudiments in place, we are now in a position to write the exact ground-state electronic energy of an N-electron system in the external field of the nuclei as described by v_ne(r),

where ρ ∈ N refers to the set of N-representable densities, Ψ → N refers to normalized N-electron multi-determinant wave functions, and ρ_Ψ denotes that our density must come from such a wave function. Minimizing the energy corresponds to solving the Schrödinger-like equation,

(T̂ + V̂^lr,ω + λV̂^sr,ω + V̂_ne + V̂^sr,ω,λ_Hxc [ρ_Ψ])|Ψ〉 = E|Ψ〉 (7)

Thus far, the treatment of the energy has been exact despite some introduced dependence on parameters ω and λ. We make our first approximation by imposing that our wave function is described by a single Slater determinant such that,

where the search is now over N-electron normalized single-determinant wave functions, Φ. Minimizing the single-determinant energy expression leads to the Kohn–Sham-like equation,⁶⁸

(T̂ + V̂_ne + V̂^lr,ω_Hx,HF [Φ] + λV̂^sr,ω_Hx,HF [Φ] + V̂^sr,ω,λ_Hxc [ρ_Φ])|Φ〉 = E₀|Φ〉 (9)

where V̂^lr,ω_Hx,HF [Φ] and V̂^sr,ω_Hx,HF [Φ] are the long-range and short-range Hartree–Fock potential operators, respectively.

From here, we make the specific choice to represent our wave function as |Ψ〉 = e^{[small tau, Greek, circumflex]}|Φ〉. This coupled cluster ansatz does not change any of the above argumentation surrounding the single-determinant approximation, but it gives us a means of incorporating electron correlation beyond Hartree–Fock in the parts of the potential that depend on the wave function. We will immediately specialize to the case of double-substitution clusters, taking [small tau, Greek, circumflex] = [small tau, Greek, circumflex]₂ where,

and â^†_p and â_p are creation and annihilation operators for orbital p, respectively.

This choice leads us to the modified Kohn–Sham equation,

(T̂ + V̂_ne + V̂^lr,ω_Hx,HF [Φ] + λV̂^sr,ω_Hx,HF [Φ] + V̂^sr,ω,λ_Hxc [ρ_Φ])e^{[small tau, Greek, circumflex]₂}|Φ〉 = E₀e^{[small tau, Greek, circumflex]₂}|Φ〉 (11)

Furthermore, we will choose a LinLCCD(hh) ansatz in which e^{[small tau, Greek, circumflex]₂} ≈ (1 + [small tau, Greek, circumflex]₂) and only hole–hole ladder diagrams are retained in the resulting amplitude equations.⁶³ This ansatz is more robust than MP2 in cases of strong correlation, and can be solved in [scr O, script letter O](n_occ⁴n_vir²) time.

Linearizing the wave function leads to the Kohn–Sham equation,

which can be solved in a Hermitian fashion by left-multiplying by 〈Φ|(1 + [small tau, Greek, circumflex]^†₂),⁶⁹then varying with respect to [small tau, Greek, circumflex]^†₂ to yield,

E = 〈Φ|(T̂ + V̂_ne + V̂^lr,ω + λV̂^sr,ω)(1 + [small tau, Greek, circumflex]₂)|Φ〉 + E^sr,ω,λ_Hxc [ρ_Φ] (14)

at stationarity. Evaluating this energy expression,

At this point we encounter an apparent difference in our short-range energy expression, which appears to only be linearly scaled by λ, while in ref. 66 they find that the short-range contribution is scaled by λ².

This apparent discrepancy is resolved on close inspection of the amplitude equations. To find the coupled-cluster amplitudes, we project onto the space of doubly-substituted determinants,

It is now clear that the amplitudes themselves carry a λ-dependence that results in factors of λ² in purely short-ranged parts of the correlation energy. This subtlety is unimportant in MP2-based DHs, but is essential in the non-linear equations of coupled-cluster theory, as the interaction-strength parameter λ scales the λ-dependent amplitudes in a non-linear fashion.

If we were to impose that our linear coupled-cluster ansatz contained only the so-called “driver” terms, our amplitude equations become equivalent to those of MP2,

Further neglecting any off-diagonal terms in the Kohn–Sham matrix elements f^q_p, we obtain,

where ε_p is the Kohn–Sham energy of orbital p. Under this MP2 approximation, we can show that the WFT-dependent part of our functional reduces to the correct form of the range-separated MP2 energy by inserting the MP2 amplitudes from eqn (18) into eqn (15) to give,where Δ^ab_ij = ε_a + ε_b − ε_i − ε_j. We re-write the above expression as,

E_c,MP2 = E^lr,ω_c,MP2 + λE^lr–sr,ω_c,MP2 + λ²E^sr,ω_c,MP2 (20)

where we have used the notation of ref. 66 for the long-range short-range cross-coupling term, E^lr–sr,ω_c,MP2. Therefore, by using our linearized coupled-cluster starting point, we are able to recover the exact form of the generalized DH functional of ref. 66 with an appropriate choice of amplitudes.

Now, we finalize the WFT portion of our functional by eliminating the long-range correlation contribution entirely. This has a host of benefits, including opening up the opportunity for aggressive screening of the integrals, as done in attenuated MP2.^64,65 This choice also allows us to learn about the importance of WFT correlation in the short range part of the Coulomb potential. Of course, this means that we will need to supplant our WFT correlation at long range with some functional that is capable of describing dispersion interactions.⁷⁰ For this, we choose the VV10 non-local correlation functional.⁷¹ By removing the long-range wave function correlation we also eliminate cross-coupling between the length scales, simplifying eqn (16) to,

which under the LinLCCD(hh) approximation⁶³ becomes,

This gives us the XC energy,

E_xc = E^lr,ω_x,HF + λE^sr,ω_x,HF + λE^sr,ω,λ_{c,LinLCCD(hh)} + E^sr,ω,λ_xc [ρ] (23)

Recall that the amplitudes are λ-dependent, so we have chosen to write the short-range LinLCCD(hh) energy as E^sr,ω,λ_{c,LinLCCD(hh)} to emphasize its own internal dependence on λ. The result is that the term λE^sr,ω,λ_{c,LinLCCD(hh)} is actually quadratic in λ.

The exchange part of the density functional is defined with the Kohn–Sham single determinant and is linear in λ,

E^sr,ω,λ_x [ρ] = 〈Φ|(1 − λ)V̂^sr,ω|Φ〉 − E^sr,ω,λ_H [ρ] = (1 − λ)E^sr,ω_x [ρ] (24)

While many approximations exist for the treatment of the short-ranged correlation energy,⁶⁶ the expression for E^sr,ω,λ_c [ρ] in the Coulomb (ω → 0) and high-density limits complements LinLCCD(hh) correlation in a few ways. Firstly, ladder correlation is particularly appropriate for describing the low-density limit of a uniform electron gas,^72,73 so it makes sense to rely on DFT correlation to describe the opposite extreme. For the second reason to make itself apparent, we write the correlation energy in the ω → 0 limit as,^74,75

E^sr,ω=0,λ_c = E_c [ρ] − λ²E_c [ρ_1/λ] (25)

where ρ_1/λ(r) = (1/λ²)ρ(r/λ) is the density with coordinates uniformly scaled by λ⁻¹. We will choose an approximation under which the density scaling is neglected,^76,77

E^sr,ω=0,λ_c ≈ (1 − λ²)E_c [ρ] (26)

Notably, the approximation in eqn (26) is sufficient to recover the exact high-density limit for an arbitrary value of ω.⁶⁶

We therefore combine our LinLCCD(hh) correlation with DFT correlation in such a way that they complement one another through a quadratic dependence on λ. Our choice of WFT correlation is most applicable in the limit of a low-density electron gas, which also nicely complements our choice of DFT correlation ansatz in eqn (26) that is most relevant in the high-density limit. Taken together, these approximations should capture some degree of strong correlation in the short-range regime, while DFT XC describes local and non-local dynamic correlation.

Lastly, we choose to modify the ωB97X⁷⁸ functional form to correctly incorporate WFT exchange and correlation, resulting in what we call ωB97X-L (L for linear-ladder correlation). The final ωB97X-L-V XC functional takes the form,

Computational details

We chose to fit the parameters for ωB97X-L-V in the Def2-ma-TZVPP basis set,^79,80 as this basis is the largest practical choice for LinLCCD(hh) correlation at the moment. The linear parameters within the density functional will account for some of the basis set incompleteness errors that are introduced by this choice. While basis set incompleteness might be somewhat mitigated due to our choice of Coulomb-attenuated operator for the WFT contribution,^64,81,82 our preliminary studies in Table S1 suggest investigation beyond this work is warranted and show a need for improving the computational efficiency of LinLCCD(hh).

We use a Euler-Maclaurin/Lebedev XC quadrature grid with 99 radial and 590 angular points to evaluate the XC potential. As ωB97X-L-V is a traditional DH functional, we first self-consistently converge the density without LinLCCD(hh) correlation to a root-mean square (RMS) error of 10⁻⁸, then we compute the LinLCCD(hh) correction using the converged Kohn–Sham orbitals. For all open-shell systems, we use restricted open-shell orbitals to converge the self-consistent field procedure, followed by a single unrestricted (Kohn–Sham) Fock build and semi-canonicalization routine prior to inputting the Kohn–Sham orbitals into an unrestricted coupled-cluster code. All non-covalent interactions calculations are done using the counterpoise correction of Boys and Bernardi.⁸³ All calculations reported in this work were carried out in a developer copy of the Q-Chem v6.3 software package.⁸⁴

We fit ωB97X-L-V to a large subset of the widely-used GMTKN55 dataset,⁸⁵ along with subsets of X40 (non-covalent interactions of halogenated systems),⁸⁶ A24 (non-covalent interactions of small molecules),⁸⁷ MOR41,⁹ ROST61,¹¹ and MOBH28.²³ We chose to fit on a subset of GMTKN55 that excluded the larger ISOL24,⁸⁸ C60ISO,⁸⁹ INV24,⁹⁰ HAL59,⁹¹ ADIM6,⁹² and UPU23⁹³ subsets for cost considerations. Our fitting set also omitted HEAVY28,⁹² four reactions from HEAVYSB11,⁸⁵ and any iodine-containing interactions in X40 in order to minimize errors imposed by effective core potentials in main-group chemistry.⁹⁴ We only included subsets of S22^95,96 and S66^97,98 as these were inherited from a smaller initial training set that we deemed insufficient. To represent transition metal thermochemistry, we fit to 27 reactions in MOR41, 24 reactions from ROST61, and 57 reaction energies along with forward/reverse barriers from MOBH28, again retaining the smallest systems. All of our abbreviated subset structures and benchmark values are reported in the SI. Overall, our fitting set has 1434 data points with 316 points reserved for testing, coming to a total of 1750 data points.

Our ωB97X-L-V functional contains several non-linear parameters including the range-separation parameter ω and the adiabatic connection parameter λ. The VV10 component also depends on parameters b and C which control short-range damping of the r⁻⁶ asymptote and the accuracy of the C₆ coefficients, respectively.⁹⁹ As the absolute errors do not change much after one least-squares optimization cycle,⁷⁸ we scanned over values of ω from 0a₀⁻¹ to 0.5a₀⁻¹, 0.2 ≤ λ ≤ 0.7, and from 6 ≤ b ≤ 10 for the VV10 b-parameter to find the lowest-error combination after one cycle. This procedure uncovered ω = 0.1a₀⁻¹, λ = 0.6, and b = 10 as the best combination. As for the VV10 damping parameter, we found that b = 10 performs within 0.1 kcal mol⁻¹ WTMAD-2 of b = 8 with results deteriorating for thermochemistry in GMTKN55 with b < 8. This result suggests that heuristically setting the b parameter in VV10 and relying on subsequent optimization of the linear parameters to compensate for errors – as done in past work on ωB97 functionals^100,101 – is indeed a suitable approach.

With the non-linear parameters determined, we employed an iterative linear least-squares fit for the remaining 15 linear parameters of the ωB97X-L functional. Before fitting, we fix the zero-order exchange coefficient c_xσ,0 = 1 − λ and the corresponding zero-order same-spin and opposite-spin correlation coefficients to c_σσ,0 = c_αβ,0 = 1 − λ² to satisfy the uniform electron gas limit. The former constraint was applied in the original definition of ωB97X,⁷⁸ but the constraints on the correlation functional are unique to this work, so we choose to call this ansatz ωB97X-L. This leaves us with 12 linear parameters to determine through least-squares fitting.

For the 12 fitted parameters, we used the following equation to guide our least-squares optimization,

W^1/2AΔx = W^1/2b (28)

implying that,

Δx = (A^TWA)⁻¹(A^TWb) (29)

where Δx = x_i−1 + x_i. The b vector is simply the signed error b_k = E_k − E_ωB97L-V, where E_k is the benchmark energy. Here, W is a diagonal matrix of training weights that were determined simply by,where the ratio on the right is that between the average root-mean square (RMS) energy for all points in the training set (〈ΔE_RMS〉) to the RMS energy of subset k (|ΔE_RMS|_k). The linear scaling factors w_k were set to 1 for all subsets except for ROST61 where w_k = 2, as the weight for ROST61 under our proposed scheme was only 0.7, which prohibited significant improvements upon fitting. Finally, the 1473 × 12 (n_points × n_coefs) A matrix is constructed with rows corresponding to the contribution of each coefficient (indexed across columns) to the total error, in a manner similar to ref. 102.

For the most flexible initial guess, we set all 12 of the linear parameters to 0 in our initial determination of the densities and Kohn–Sham orbitals. We then update the parameters using the scheme outlined above, followed by self-consistent re-optimization of the densities and Kohn–Sham orbitals, repeating this process until the parameters were well converged. Convergence was determined by a change in WTMAD-2 of less than 10 μE_H, leading to an RMS change in parameters of ∼10⁻², which is consistent with previous work.^78,103 This procedure required four cycles to converge.

Results & discussion

Throughout the coming analysis, we will cast our results against the backdrop of ref. 104, which documented the performance of 29 unique DH approximations on the GMTKN55 database. From this set, we focus on the semi-empirical class of functionals, as these were found to outperform non-empirical DHs and our ωB97X-L-V approach falls under the umbrella of the former. We will add to this set with a comparison against more recent DH functionals such as ωB97M(2),¹⁰⁰ ωDSD72-PBEP86-D4,¹⁰⁵ DH24,¹⁰⁶ and the non-empirical Pr²SCAN69-D4,¹⁰⁷ as their good performance on GMTKN55 sets the bar for high-accuracy DHs.

While this set is far from comprehensive, we compare against the 14 semi-empirical (and one non-empirical) DHs highlighted in Table 1. Table 1 also categorizes DHs by their use of additional empiricism through spin-component scaling of the WFT correlation. We emphasize that ωB97X-L-V is a semi-empirical DH that makes no use of spin-component scaling and that our objective here is to examine the performance of a DH based solely on short-range LinLCCD(hh) correlation without additional parameterization of the WFT component. Thus, semi-empirical DH functionals without spin-component scaling are the most valid comparison to ωB97X-L-V.

Table 1 DH functionals assessed in this work

Functional	Spin-component scaling	Ref.
a GMTKN55 data from ref. 104.b GMTKN55 data from ref. 116.c GMTKN55 data from this work.d GMTKN55 data from ref. 105.e GMTKN55 data from ref. 107.f GMTKN55 data from ref. 106.g Non-empirical functional.
B2PLYP^a	✗	5
mPW2PLYP^a	✗	108
B2GPPLYP^a	✗	109
B2NC-PLYP^a	✗	110
mPW2NC-PLYP^a	✗	110
B2PPW91^a	✗	111
ωB97M(2)^b	✗	100
ωB97X-L-V^c	✗	This work
mSD-PBEPBE^a	✓	13 and 112
PWPB95^a	✓	113
DSD-PBEP86^a	✓	114
DSD-PBEB95^a	✓	12
DSD-BLYP^a	✓	115
ωB97X-2^a	✓	101
ωDSD72-PBEP86-D4^d	✓	105
Pr^2SCAN69-D4^e^g	✓	107
DH24^f	✓	106

---

We will use the weighted mean absolute deviation (WTMAD-2) to compare the performance across functionals. The WTMAD-2 metric is defined as,⁸⁵

where is the averaged absolute energy in subset i, 56.84 kcal mol⁻¹ is the average over all 55 values, and N_i and MAD_i are the number of data points and mean absolute deviation for subset i, respectively. The WTMAD-2 should not be confused with the MAD as it is generally larger and reflects different information. Instead, the WTMAD-2 provides a statistically meaningful way of comparing results across many data sets for different density functionals. The WTMAD-2 can also be broken down into contributions from individual subsets, where the factor 56.84 kcal mol⁻¹ is retained, but the summations are truncated to include only the subsets of interest.

Optimized parameters for ωB97X-L-V

The parameters at each least-squares cycle for ωB97X-L-V are shown in Table 2. While the linear parameters are unconstrained, and therefore yield very little in physically-interpretable results, we can learn from the values of our other parameters. We note that ω = 0.1 a₀⁻¹ is quite satisfying, as we employed the approximation that ω → 0 to arrive at eqn (26), and the optimal ω is small. Our optimized λ = 0.6 is also close to the average coefficient of Hartree–Fock exchange used across double-hybrid functionals (λ = 0.64).¹⁰⁰ This value of λ also puts our percentage of WFT correlation at 36%, which again is close to the DH average of 32%.

Table 2 Optimized parameters for ωB97X-L-V

Parameter	Cycle 1	Cycle 2	Cycle 3	Final
Fixed parameters
b	10.0	10.0	10.0	10.0
C	0.01	0.01	0.01	0.01
ω	0.1a0^−1	0.1a0^−1	0.1a0^−1	0.1a0^−1
λ	0.6	0.6	0.6	0.6
cxσ,0	0.4	0.4	0.4	0.4
ccσσ,0	0.64	0.64	0.64	0.64
ccα β,0	0.64	0.64	0.64	0.64
Least-squares optimized parameters
cxσ,1	0.0	0.136	0.152	0.154
ccσσ,1	0.0	−1.399	−1.415	−1.417
ccαβ,1	0.0	−1.151	−1.194	−1.194
cxσ,2	0.0	3.910	3.901	3.884
ccσσ,2	0.0	4.585	4.697	4.716
ccαβ,2	0.0	−1.512	−1.364	−1.348
cxσ,3	0.0	−11.362	−11.357	−11.300
ccσσ,3	0.0	−2.703	−2.918	−2.956
ccαβ,3	0.0	12.659	11.964	11.869
cxσ,4	0.0	8.637	8.482	8.425
ccσσ,4	0.0	−0.985	−0.885	−0.861
ccαβ,4	0.0	−10.594	−9.673	−9.571

---

Regarding the quality of fit and potential for transferability, the average MAD across subsets in our fitting set was 1.4 kcal mol⁻¹ and the corresponding average MAD for our test set was 0.9 kcal mol⁻¹. We also computed modified WTMAD-2 values for the entirety of our fitting set (〈ΔE_RMS〉 = 66.158 kcal mol⁻¹) and for the test set (〈ΔE_RMS〉 = 24.208 kcal mol⁻¹), at 3.0 kcal mol⁻¹ and 1.6 kcal mol⁻¹, respectively. We remark that due to a large discrepancy in size and total energies for the reactions in our fitting and test sets, we needed to use different weights for the WTMAD-2 calculation to obtain meaningful results. The overall consistent performance between fitting and test sets suggests that our parameters may be reasonably transferable.

Performance on GMTKN55

Although we trained on a large portion of GMTKN55, we assess the performance of ωB97X-L-V on all 55 subsets to understand areas where our functional could be most improved. We emphasize that our GMTKN55 results should be interpreted with this in mind. To analyze our functional's performance more deeply, we break GMTKN55 into the usual five main categories: basic properties and reaction energies for small systems (BASIC), reaction energies for large systems and isomerisation reactions (LARGE + ISO), reaction barrier heights (BARRIERS), and inter/intra-molecular non-covalent interaction energies (INT/INTRA-NE).

The WTMAD-2 data for the BASIC category in Fig. 1 show that ωB97X-L-V performs about as well as other DH functionals. Overall, we report a WTMAD-2 for ωB97X-L-V of 2.1 kcal mol⁻¹ for the BASIC category. It is encouraging to see that ωB97X-L-V performs on par with the average WTMAD-2 for functionals that do not leverage spin-component scaling.

Fig. 1 WTMAD-2 across double-hybrid density functionals for basic properties and reaction energies in GMTKN55. Green bars indicate spin-component scaled functionals and orange bars are functionals without spin-component scaling.

Next, we consider the WTMAD-2 results for the LARGE + ISO category in Fig. 2. Here, ωB97X-L-V performs well with respect to functionals that omit spin-component scaling with a WTMAD-2 of 3.7 kcal mol⁻¹. We attribute the reasonable performance of ωB97X-L-V on LARGE + ISO to the fact that LinLCCD(hh) correlation outperforms MP2 correlation on large systems⁶³ for which MP2 diverges.¹¹⁷ Notably, all of the MP2-based DHs that perform best on the LARGE + ISO category make use of spin-component scaling. This tempers the divergent MP2 correlation contribution somewhat, allowing functionals like those in the DSD family to outperform simpler DHs that scale the MP2 correlation uniformly. Overall, the good performance of ωB97X-L-V without spin-component scaling is encouraging, suggesting there are gains to be made with such a treatment of the WFT contribution, and the results for MP2-based DHs suggest that adapting ωB97X-L-V with spin-component scaling may lead to further improvements.

Fig. 2 WTMAD-2 across double-hybrid density functionals for reaction energies of large systems and isomerization energies in GMTKN55. Green bars indicate spin-component scaled functionals and orange bars are functionals without spin-component scaling.

The results for the BARRIERS category are plotted in Fig. 3. Here, we report a WTMAD-2 of 2.5 kcal mol⁻¹ that once again performs similarly to other DHs that make no use of spin-component scaling. ωB97X-L-V is outdone by spin-component scaled functionals, suggesting again that future performance upgrades may be possible with a spin-component scaled adaptation.

Fig. 3 WTMAD-2 across double-hybrid density functionals for reaction barriers in GMTKN55. Green bars indicate spin-component scaled functionals and orange bars are functionals without spin-component scaling.

We conclude the analysis of GMTKN55 subsets with a discussion of the INTER-NE and INTRA-NE categories in Fig. 4 and 5, respectively. Overall, we find that non-covalent interactions (both inter- and intra-molecular) are a point where ωB97X-L-V could be improved the most. For inter- and intra-molecular interactions, ωB97X-L-V attains WTMAD-2 values of 3.9 kcal mol⁻¹ 3.8 kcal mol⁻¹, respectively. This result may have been expected, as ωB97X-L-V only incorporates WFT correlation at short-range, turning into more of a typical range-separated hybrid functional with non-local correlation in the long-range limit where non-covalent interactions are most relevant. This hypothesis is born out, as the rung-4 ωB97X-V functional has a WTMAD-2 of 3.0 kcal mol⁻¹ and 3.6 kcal mol⁻¹ for intra- and inter-molecular interactions, respectively.⁸⁵ Furthermore, the rung-5 ωB97X-2-D3 functional has a WTMAD-2 of 2.3 kcal mol⁻¹ for intra-molecular interactions and 5.3 kcal mol⁻¹ for inter-molecular interactions. Put together, ωB97X-L-V performs about as expected for non-covalent interactions. Overall, we cannot expect the benefits of the fifth-rung to significantly influence non-covalent interactions for ωB97X-L-V. This result motivates future parameterizations of ωB97X-L-V that perform better for non-covalent interactions.

Fig. 4 WTMAD-2 across double-hybrid density functionals for intermolecular interaction energies in GMTKN55. Green bars indicate spin-component scaled functionals and orange bars are functionals without spin-component scaling.

Fig. 5 WTMAD-2 across double-hybrid density functionals for intramolecular interaction energies in GMTKN55. Green bars indicate spin-component scaled functionals and orange bars are functionals without spin-component scaling.

Finally, we show the overall performance on the GMTKN55 data set in Fig. 6. Over the entire GMTKN55 data set, ωB97X-L-V performs reasonably well with a WTMAD-2 of 3.1 kcal mol⁻¹. This is very close to the best functional tested in ref. 104 (ωB97X-2-D3(BJ), WTMAD-2 = 3.0 kcal mol⁻¹), which is quite sensible as the underlying density functional is similar. Overall, we report that using attenuated LinLCCD(hh) correlation in a DH such as ωB97X-L-V can lead to a DH functional that is about as reliable for equilibrium main-group chemistry as MP2-based DHs.

Fig. 6 WTMAD-2 across double-hybrid density functionals for all of GMTKN55. Green bars indicate spin-component scaled functionals and orange bars are functionals without spin-component scaling.

Performance on transition-metal complexes

Our primary objective when we began designing ωB97X-L-V was not necessarily to improve upon the main-group chemistry represented in the GMTKN55 database, but to more accurately describe strongly-correlated systems. One domain in which conventional DHs can yield spurious results is in the description of transition metal complexes.⁷ To understand the robustness of ωB97X-L-V on such systems, we now evaluate our functional on MOR41, ROST61, and MOBH28 data sets.

In order to compare across functionals, we compute the WTMAD without energy scaling for these transition metal sets as done by Grotjahn and Kaupp,²⁹

We have chosen to compare ωB97X-L-V with PWPB95, ωB97M(2), and SOS-DH24 functionals, as the latter two are exemplary for GMTKN55 and PWPB95 is consistently the best for metal–organic reaction energies.^9,11 While data for MOR41 and ROST61 are abundantly available for many functionals, we needed to generate the PWPB95-D3(BJ) data for MOBH28. For this, we used the practical Def2-ma-TZVPP basis set.

The MAD and WTMAD results in Table 3 suggest that ωB97X-L-V is among these top performers on GMTKN55 for transition-metal complexes. Notably, 65% of MOR41 and <50% of ROST61 and MOBH28 were used to parameterize ωB97X-L-V, meaning that most of the reactions in this section come from outside our training set. Despite this, the MAD for ROST61 decreased by 0.5 kcal mol⁻¹ upon adding the test set data while the MAD for MOR41 and MOBH28 increased by no more than 0.1 kcal mol⁻¹, hinting that the ωB97X-L-V parameters may be reasonably transferable. Furthermore, transition metal complexes constituted only 7% of our training data, so we expect this to be a reasonable representation of the performance of ωB97X-L-V on such systems.

Table 3 Performance for representative functionals on transition-metal data sets (all results in kcal mol⁻¹)

Data set	MAD
	PWPB95-D3(BJ)	ωB97M(2)	SOS-DH24^a	ωB97X-L-V
a From ref. 106.b From ref. 9.c From ref. 11.d From ref. 10.e From ref. 118.
MOR41	1.9^b	2.8^d	2.2	2.7
ROST61	1.9^c	2.0^e	1.9	2.0
MOBH28	1.1	1.0^e	1.2	0.9
WTMAD	1.6	1.9	1.7	1.8

---

Across all three sets, ωB97X-L-V performs similarly to the highly-successful ωB97M(2) functional. Our comparison point of PWPB95-D3(BJ) was an excellent choice as its newly-reported MAD for MOBH28 of 1.1 kcal mol⁻¹ is lower than SOS-DH24, which itself yields the excellent MAD of 1.2 kcal mol⁻¹.¹⁰⁶ This low MAD of 1.1 kcal mol⁻¹ matches that of Pr²SCAN69,¹⁰⁷ another notably excellent functional for MOBH28. However, recall that PWPB95-D3(BJ) does not perform nearly as well for main-group chemistry as shown in Fig. 6. Finally, we note that one must beware that the differences in basis sets used to generate the data in Table 3 could potentially contaminate the direct comparison across approaches.

We now examine the performance of ωB97X-L-V on a set of experimental 3d transition-metal diatomic bond dissociation energies (BDEs) from ref. 119. We compare the performance of ωB97X-L-V against several of the best functionals on the main-group chemistry of GMTKN55. Top-performing spin-component scaled and non-spin-scaled functionals are included for additional contrast.

The 3d metal complexes from ref. 119 constitute a truly multireference testing ground for our ωB97X-L-V functional that lies entirely outside our training set. As an assessment of the multireference character, we report values for the T₁ metric for all 3d-metal complexes in the SI using the Def2-QZVPPD basis set. This screening suggests that 22 of the 43 metal complexes are at least moderately correlated (T₁ > 0.02), while another 19 complexes are strongly-correlated (T₁ > 0.05). While no single metric is perfect, the T₁ metric suggests that only two of these complexes are considered single-reference.

Our results in Fig. 7 suggest that ωB97X-L-V is indeed robust in cases of strong correlation, with a MAD of 3.8 kcal mol⁻¹. Notably, if we consider only the best performers on the GMTKN55 database, ωB97X-L-V is among the top two DH functionals for this set of 3d-metal BDEs with ωB97M(2) performing surprisingly well at a MAD of 3.1 kcal mol⁻¹. Importantly, ωB97X-L-V achieves a lower MAD than the related ωB97X-2-D3 functional. A source of error common to both of these ωB97X-based DH functionals is a systematic underestimation of the BDE of metal oxides, which appears to be remedied somewhat with DH meta-GGAs. However, the ωB97X-L-V max error is about half of that for ωB97X-2-D3, suggesting some improved robustness in multireference cases.

Fig. 7 Mean absolute deviation (MAD) of various functionals on 43 bond dissociation energies of 3d transition-metal diatomic species (top panel), and the associated maximum error for each functional (bottom panel). Green bars indicate spin-component scaled functionals and orange bars are functionals without spin-component scaling.

Owing to multiple instances of anomalously high errors that correspond almost entirely to the T₁ > 0.05 subset, the ωDSD-PBEP86-D4 and Pr²SCAN69-D4 functionals have MADs at least 2× larger than those of ωB97X-L-V. The largest errors across all functionals come from strongly-correlated metal oxides, but the magnitude of these errors is tempered significantly in ωB97X-L-V. Because ωB97X-2-D3 seems just as afflicted by these errors as other functionals, we conclude that the ωB97X-L-V robustness to static correlation leads to improved performance on strongly-correlated transition-metal diatomics.

While our parameterization of ωB97X-L-V could potentially be improved by resorting to spin-component scaling, its wholesale incorporation of attenuated-LinLCCD(hh) correlation rather than MP2 can teach us about essential elements of functional design. Despite its relatively average performance on main-group chemistry, the fact that ωB97X-L-V performs well on metal complexes implies that alternatives MP2 correlation in the short-range could be a beneficial avenue for future functional development.

Covalent bond breaking

To test ωB97X-L-V in a domain where other DHs fail catastrophically, we will examine bond dissociation energies (BDEs) for several molecules using the Def2-ma-TZVPP basis set. In particular, we will examine the BDEs of 19 closed-shell singlet systems curated from the multireference BDE subset of W4-11¹²⁰ alongside several others compiled from high-accuracy experimental data and CCSDTQ/CBS results (details in the SI). We compute the BDEs using restricted orbitals by the difference E (r = 10 000 Å) − E (r_eq), which is similar in spirit to the DISS10¹²¹ data set. Notably, much better results can be obtained by allowing for spin-symmetry breaking through the use of unrestricted reference orbitals, but the restricted case serves as a stringent test of methods in statically-correlated systems. This is because in covalent bond breaking, one encounters absolute near-degeneracy correlation,¹²² which is a potent case of static correlation that all but ensures the failure of MP2-based DHs. While we focus on asymptotic BDEs here, we also illustrate the performance of ωB97X-L-V for H₂, HF, and F₂ bond dissociation curves within Fig. S1 as well as the torsional profile of ethylene in Fig. S2.

Table 4 shows the BDE results for ωB97M(2), ωDSD72-PBEP86-D4, Pr²SCAN69-D4, ωB97X-2-D3, and ωB97X-L-V functionals. As expected, DH approximations built atop an MP2 ansatz diverge in over half of the systems. Notably, half of the divergent cases for MP2-based DH functionals correspond to asymmetrical bond cleavage, such as H⋯OF. Our ωB97X-L-V functional provides similar results to the other four DH functionals when none of the methods diverge. However, ωB97X-L-V provides finite results for every system tested, including the challenging case of the dinitrogen triple bond. Of course, for a quantitative picture of dinitrogen dissociation, one must include hextuple excitations, so despite the large error it is nonetheless noteworthy that ωB97X-L-V produces a finite result at all. Overall, ωB97X-L-V obtains a MAD of 53 kcal mol⁻¹ across all systems, and 33 kcal mol⁻¹ for BDEs of single bonds.

Table 4 Bond dissociation energies computed via E (r = 10 000 Å) − E (r_eq) (kcal mol⁻¹). Dashes indicate divergent results

System	Benchmark	Signed error
		ωB97M(2)	ωDSD72-PBEP86-D4	Pr^2SCAN69-D4	ωB97X-2-D3	ωB97X-L-V
a Computed at CCSD/aug-cc-pV5Z level.b Experimental data from ref. 124.c Computed with the Def2-QZVPPD basis set.d Computed at CCSDTQ/CBS level.e Experimental data from ref. 125.f Experimental data from ref. 126.g Experimental data from ref. 127.h Theoretical best estimate from ref. 128.i Experimental data from ref. 129, corrected with experimental zero-point energy from ref. 130.j Experimental data from ref. 131.k Theoretical best estimate from ref. 132.l Theoretical best estimate from ref. 120.
H⋯H	109.4^a	—	—	—	—	46.3
Be⋯Be	2.7^b	−3.6^c	−1.3^c	−0.2^c	−1.9^c	−0.4^c
B⋯B	52.7^d	—	—	—	—	16.4
C⋯C	142.5^e	62.2	37.1	45.3	31.6	49.7
N⋯N	228.4^f	—	—	—	—	167.8
O⋯O	93.5^d	115.5	120.7	118.6	117.5	118.9
F⋯F	37.6^g	—	—	—	—	16.8
Cl⋯Cl	56.5^h	—	—	—	—	30.0
OC⋯O	129.3^i	61.7	65.4	67.5	63.6	68.5
C⋯O	259.0^j	100.3	102.3	104.2	101.7	105.9
H⋯F	141.0^k	—	—	—	—	45.7
F⋯OOF	17.7^l	—	—	—	—	30.1
FO⋯OF	46.2^l	—	—	—	—	33.3
ClO⋯Cl	36.0^l	—	—	—	—	29.9
FO⋯F	40.7^l	—	—	—	—	23.9
H⋯OF	105.6^l	—	—	—	—	38.5
O2⋯O	26.6^l	89.5	94.9	93.9	91.0	92.9
S3⋯S	66.0^l	46.2	49.8	52.4	52.0	50.8
S2⋯ S2	25.9^l	30.9	35.3	37.4	35.8	36.0
MAD	—	—	—	—	—	52.7
MAD (single bonds)	—	—	—	—	—	33.1

---

Another interesting case is Be₂, which has an incredibly weak bond that is primarily due to static correlation.¹²³ Its weak, statically correlated character makes Be₂ a unique challenge for modern DFT and many functionals predict that Be₂ is simply unbound. Due to the high precision required, basis set errors were suppressed by using the Def2-QZVPPD functional. While the ωB97M(2) functional predicts that Be₂ is unbound, ωDSD72-PBEP86-D4 and ωB97X-2-D3 predict very small wells while Pr²SCAN69-D4 comes within 0.2 kcal mol⁻¹ of the correct well depth. ωB97X-L-V also provides a close estimate, within 0.4 kcal mol⁻¹ of the very weak 2.7 kcal mol⁻¹ BDE.

One might wonder how ωB97X-L-V improves BDEs computed in this way, given that its descriptor of static correlation (LinLCCD(hh)) is confined to the short-range part of the Coulomb operator. This is because the relevant correlation energy is not long-ranged (the Coulomb interaction between two infinitely-separated fragments for any size-consistent method is zero). For example, after converging an unrestricted Hartree–Fock calculation using the usual r⁻¹ Coulomb operator, Coulomb-attenuated CCSD with ω = 0.1a₀⁻¹ recovers 97% of the asymptotic CCSD correlation energy in F₂ because ω need only be set to recover correlation in the vicinity of each individual F atom. So, because ωB97X-L-V has the correct asymptotic exchange potential and LinLCCD(hh) short-range correlation, it is possible to accurately describe asymptotic BDEs without the catastrophic divergence of MP2-based functionals.

Conclusions

In this work, we introduced a new double-hybrid density functional called ωB97X-L-V that leverages attenuated LinLCCD(hh)⁶³ instead of the usual MP2 correlation. We justify the inclusion of one parameter to control the interface of DFT and LinLCCD(hh) correlation by means of adiabatic connection, but we adapt the parameters of ωB97X to best interface with LinLCCD(hh), making ωB97X-L-V a semi-empirical DH functional. The parameters were determined by fitting to a subset of the GMTKN55 database augmented with thermochemical data for transition-metal complexes.

To test ωB97X-L-V, we analyzed its performance on the whole GMTKN55 database, which revealed that it can perform about as well as MP2-based double-hybrid functionals for nearly every property (basic thermochemistry, large system reaction and isomerization energies, and barrier heights). We also found that ωB97X-L-V could be most improved for non-covalent interactions in main-group chemistry. This was expected, as ωB97X-L-V only leverages WFT correlation in the short-range part of the Coulomb potential. Over all of GMTKN55, ωB97X-L-V performed as well as most other double-hybrids that we examined, which is promising as it suggests that resorting to attenuated LinLCCD(hh) correlation instead of MP2 does not necessarily damage the results that are attainable for main-group chemistry.¹⁰⁴

For transition-metal reaction energies ωB97X-L-V appears to exceed the expectations set by its performance on main-group chemistry. This is especially true for transition-metal reaction barrier heights of MOBH28, for which ωB97X-L-V yields a MAD of 0.9 kcal mol⁻¹. The performance of ωB97X-L-V on transition-metal BDEs also remained similar to high-performing DHs with a notable reduction in errors relative to the similar ωB97X-2-D3 functional. Lastly, on main-group BDEs computed with restricted orbitals as E (r = 10 000 Å) − E (r_eq), we found that the use of LinLCCD(hh) correlation can provide finite results whereas MP2-based DH functionals diverge.

There remain many exciting avenues of future work on the ωB97X-L-V functional. Our group is currently testing a much faster algorithm for solving the LinLCCD(hh) amplitude equations (originally suggested in ref. 63), which leads to at least one order-of-magnitude speedups in our nascent implementation. We are also interested in parameterizations that leverage the complementary two error function (terfc) that is used in attenuated MP2 methods,^64,65 which can further improve the results of integral screening. Under the umbrella of accuracy considerations, we are currently pursuing spin-component scaled versions of LinLCCD(hh) to interface with a future parameterization of ωB97X-L-V. Overall, we view ωB97X-L-V as an intriguing scaffold upon which to build even more sophisticated density functional approximations that push the paradigm of double-hybrid DFT beyond MP2 correlation.

Author contributions

Kevin Carter-Fenk: conceptualization (lead), software (equal), supervision (lead), formal analysis (equal), investigation (equal), writing – original draft (equal), writing – review & editing (equal). Ella R. Ransford: software (equal), formal analysis (equal), investigation (equal), writing – original draft (equal), writing – review & editing (equal).

Conflicts of interest

There are no conflicts of interest to declare.

Data availability

All geometries and benchmark data for the abbreviated X40, S22, S66, A24, MOR41, ROST61, and MOBH28 sets used for fitting are provided (ZIP). Geometries for the ethylene torsion potential surface are provided (ZIP). This study was carried out using publicly available geometries and reference values from the GMTKN55, A24, X40, MOR41, ROST61, and MOBH28 data bases available at the corresponding references in the main text. All data required to reproduce the bar charts in the main text are reported (XLSX) in the supplementary information (SI) alongside all data for bond dissociation energies of 3d transition metal complexes. All error statistics and functional parameters are reported in the main text. Supplementary information: all geometries and reference data for abbreviated S22, S66, X40, A24, MOR41, ROST61, and MOBH28 subsets used for fitting. Potential energy surfaces for H₂, F₂, and HF bond dissociation and ethylene torsion. All data required to reproduce all bar charts presented in the text. T₁ metrics, zero-point vibrational corrections, spin–orbit corrections, and errors for all transition-metal bond dissociation energy data. Additional computational details. See DOI: https://doi.org/10.1039/d6cp00232c.

Acknowledgements

This work was supported by the University of Pittsburgh and the Center for Research Computing and Data RRID: SCR_022735, through the resources provided. Specifically, this work used the H2P cluster, which is supported by NSF award number OAC-2117681. E. R. R. acknowledges support from the Wass Undergraduate Research Fellowship. K. C.-F. acknowledges a fruitful discussion with Prof. Kieron Burke that helped guide our choices for our training set. K. C.-F. is also grateful for enthusiastic discussions with Dr Kimberly Daas and Dr Steven Crisostomo that sparked wonder and inspiration.

Notes and references

1. J. P. Perdew, MRS Bull., 2013, 38, 743–750.
2. J. C. Sancho-García and C. Adamo, Phys. Chem. Chem. Phys., 2013, 15, 14581–14594.
3. L. Goerigk and S. Grimme, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2014, 4, 576–600.
4. S. Ghosh, P. Verma, C. J. Cramer, L. Gagliardi and D. G. Truhlar, Chem. Rev., 2018, 118, 7249–7292.
5. S. Grimme, J. Chem. Phys., 2006, 124, 034108.
6. T. Benighaus, R. A. DiStasio, Jr., R. C. Lochan, J.-D. Chai and M. Head-Gordon, J. Phys. Chem. A, 2008, 112, 2702–2712.
7. J. Shee, M. Loipersberger, D. Hait, J. Lee and M. Head-Gordon, J. Chem. Phys., 2021, 154, 194109.
8. W. Jiang, M. L. Laury, M. Powell and A. K. Wilson, J. Chem. Theory Comput., 2012, 8, 4102–4111.
9. S. Dohm, A. Hansen, M. Steinmetz, S. Grimme and M. P. Checinski, J. Chem. Theory Comput., 2018, 14, 2596–2608.
10. M. A. Iron and T. Janes, J. Phys. Chem. A, 2019, 123, 3761–3781.
11. L. R. Maurer, M. Bursch, S. Grimme and A. Hansen, J. Chem. Theory Comput., 2021, 17, 6134–6151.
12. S. Kozuch and J. M. L. Martin, J. Comput. Chem., 2013, 34, 2327–2344.
13. L. M. Roch and K. K. Baldridge, Phys. Chem. Chem. Phys., 2017, 19, 26191–26200.
14. G. Santra, N. Sylvetsky and J. M. L. Martin, J. Phys. Chem. A, 2019, 123, 5129–5143.
15. J. M. L. Martin and G. Santra, Isr. J. Chem., 2020, 60, 787–804.
16. G. Santra, M. Bursch and L. Wittmann, J. Phys. Chem. A, 2025, 129, 7218–7228.
17. P. Morgante and R. Peverati, J. Comput. Chem., 2019, 40, 839–848.
18. P. Morgante and R. Peverati, Molecules, 2023, 28, 3487.
19. G. Schoendorff and A. K. Wilson, J. Chem. Phys., 2014, 140, 224314.
20. S. Grimmel, G. Schoendorff and A. K. Wilson, J. Chem. Theory Comput., 2016, 12, 1259–1266.
21. L. E. Aebersold, S. H. Yuwono, G. Schoendorff and A. K. Wilson, J. Chem. Theory Comput., 2017, 13, 2831–2839.
22. A. Rajabi, R. Grotjahn, D. Rappoport and F. Furche, Dalton Trans., 2024, 53, 410–417.
23. E. Semidalas and J. M. Martin, J. Chem. Theory Comput., 2022, 18, 883–898.
24. N. J. DeYonker, K. A. Peterson, G. Steyl, A. K. Wilson and T. R. Cundari, J. Phys. Chem. A, 2007, 111, 11269–11277.
25. A. D. Becke, J. Chem. Phys., 2005, 122, 064101.
26. A. D. Becke, J. Chem. Phys., 2013, 138, 074109.
27. J. Kong and E. Proynov, J. Chem. Theory Comput., 2015, 12, 133–143.
28. E. R. Johnson and A. D. Becke, J. Chem. Phys., 2017, 146, 211105.
29. R. Grotjahn and M. Kaupp, Isr. J. Chem., 2022, 63, e202200021.
30. M. Kaupp, A. Wodyński, A. V. Arbuznikov, S. Fürst and C. J. Schattenberg, Acc. Chem. Res., 2024, 57, 1815–1826.
31. A. Wodyński and M. Kaupp, J. Chem. Theory Comput., 2025, 21, 7419–7429.
32. A. M. M. Kai, T. T. Nguyen and R. Grotjahn, J. Chem. Theory Comput., 2025, 21, 11004–11020.
33. G. L. Manni, R. K. Carlson, S. Luo, D. Ma, J. Olsen, D. G. Truhlar and L. Gagliardi, J. Chem. Theory Comput., 2014, 10, 3669–3680.
34. L. Gagliardi, D. G. Truhlar, G. L. Manni, R. K. Carlson, C. E. Hoyer and J. L. Bao, ACS Cent. Sci., 2017, 50, 66–73.
35. P. Sharma, J. J. Bao, D. G. Truhlar and L. Gagliardi, Annu. Rev. Phys. Chem., 2021, 72, 541–564.
36. A. J. Garza, A. G. Sousa Alencar and G. E. Scuseria, J. Chem. Phys., 2015, 143, 244106.
37. A. J. Garza, I. W. Bulik, T. M. Henderson and G. E. Scuseria, Phys. Chem. Chem. Phys., 2015, 17, 22412–22422.
38. A. J. Garza, I. W. Bulik, A. G. S. Alencar, J. Sun, J. P. Perdew and G. E. Scuseria, Mol. Phys., 2016, 114, 997–1018.
39. G. E. Scuseria, T. M. Henderson and D. C. Sorensen, J. Chem. Phys., 2008, 129, 231101.
40. G. Jansen, R.-F. Liu and J. G. Ángyán, J. Chem. Phys., 2010, 133, 154106.
41. G. P. Chen, V. K. Voora, M. M. Agee, S. G. Balasubramani and F. Furche, Annu. Rev. Phys. Chem., 2017, 68, 421–445.
42. F. Furche, Phys. Rev. B: Condens. Matter Mater. Phys., 2001, 64, 195120.
43. M. Fuchs, Y.-M. Niquet, X. Gonze and K. Burke, J. Chem. Phys., 2005, 122, 094116.
44. B. G. Janesko, T. M. Henderson and G. E. Scuseria, J. Chem. Phys., 2009, 130, 081105.
45. J. Toulouse, W. Zhu, J. G. Ángyán and A. Savin, Phys. Rev. A: At., Mol., Opt. Phys., 2010, 82, 032502.
46. W. Zhu, J. Toulouse, A. Savin and J. G. Ángyán, J. Chem. Phys., 2010, 132, 244108.
47. P. Verma and R. J. Bartlett, J. Chem. Phys., 2012, 136, 044105.
48. H. Eshuis, J. E. Bates and F. Furche, Theor. Chem. Acc., 2012, 131, 1084.
49. A. Kreppel, D. Graf, H. Laqua and C. Ochsenfeld, J. Chem. Theory Comput., 2020, 16, 2985–2994.
50. D. Graf and A. J. W. Thom, J. Chem. Phys., 2023, 159, 174106.
51. C. Sheldon, J. Paier, D. Usvyat and J. Sauer, J. Chem. Theory Comput., 2024, 20, 2219–2227.
52. D. Stück and M. Head-Gordon, J. Chem. Phys., 2013, 139, 244109.
53. R. M. Razban, D. Stück and M. Head-Gordon, Mol. Phys., 2017, 115, 2102–2109.
54. J. Lee and M. Head-Gordon, J. Chem. Theory Comput., 2018, 14, 5203–5219.
55. J. Shee, M. Loipersberger, A. Rettig, J. Lee and M. Head-Gordon, J. Phys. Chem. Lett., 2021, 12, 12084–12097.
56. G. Santra and J. M. L. Martin, J. Phys. Chem. Lett., 2022, 13, 3499–3506.
57. D. Peng, S. N. Steinmann, H. Van Aggelen and W. Yang, J. Chem. Phys., 2013, 139, 104112.
58. G. E. Scuseria, T. M. Henderson and I. W. Bulik, J. Chem. Phys., 2013, 139, 104113.
59. H. Van Aggelen, Y. Yang and W. Yang, Phys. Rev. A: At., Mol., Opt. Phys., 2013, 88, 030501.
60. H. van Aggelen, Y. Yang and W. Yang, J. Chem. Phys., 2014, 140, 18A511.
61. J. Li, Y. Jin, J. Yu, W. Yang and T. Zhu, J. Chem. Theory Comput., 2024, 20, 7979–7989.
62. J. Li, Y. Jin, J. Yu, W. Yang and T. Zhu, J. Phys. Chem. Lett., 2024, 15, 2757–2764.
63. K. Carter-Fenk, J. Phys. Chem. A, 2025, 129, 7251–7260.
64. M. Goldey and M. Head-Gordon, J. Phys. Chem. Lett., 2012, 3, 3592–3598.
65. M. Goldey, A. Dutoi and M. Head-Gordon, Phys. Chem. Chem. Phys., 2013, 15, 15869–15875.
66. C. Kalai and J. Toulouse, J. Chem. Phys., 2018, 148, 164105.
67. M. Levy, Proc. Natl. Acad. Sci. U. S. A., 1979, 76, 6062–6065.
68. W. Kohn and L. J. Sham, Phys. Rev., 1965, 140, A1133–A1138.
69. A. G. Taube and R. J. Bartlett, J. Chem. Phys., 2009, 130, 144112.
70. S. Grimme, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2011, 1, 211–228.
71. O. A. Vydrov and T. Van Voorhis, J. Chem. Phys., 2010, 133, 244103.
72. J. J. Shepherd, T. M. Henderson and G. E. Scuseria, Phys. Rev. Lett., 2014, 112, 133002.
73. J. J. Shepherd, T. M. Henderson and G. E. Scuseria, J. Chem. Phys., 2014, 140, 124102.
74. M. Levy, Phys. Rev. A: At., Mol., Opt. Phys., 1991, 43, 4637–4646.
75. M. Levy, Density Functional Theory, Plenum Press, New York, 1995, pp. 11–32.
76. K. Sharkas, J. Toulouse and A. Savin, J. Chem. Phys., 2011, 134, 064113.
77. K. Sharkas, A. Savin, H. J. A. Jensen and J. Toulouse, J. Chem. Phys., 2012, 137, 044104.
78. J.-D. Chai and M. Head-Gordon, J. Chem. Phys., 2008, 128, 084106.
79. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
80. M. Gray and J. M. Herbert, J. Chem. Theory Comput., 2022, 18, 2308–2330.
81. M. Goldey and M. Head-Gordon, Chem. Phys. Lett., 2014, 608, 249–254.
82. Y. Huang, M. Goldey, M. Head-Gordon and G. J. O. Beran, J. Chem. Theory Comput., 2014, 10, 2054–2063.
83. S. F. Boys and F. Bernardi, Mol. Phys., 1970, 19, 553–566.
84. E. Epifanovsky, et al., J. Chem. Phys., 2021, 155, 084801.
85. L. Goerigk, A. Hansen, C. Bauer, S. Ehrlich, A. Najibi and S. Grimme, Phys. Chem. Chem. Phys., 2017, 19, 32184–32215.
86. J. Řezác, K. E. Riley and P. Hobza, J. Chem. Theory Comput., 2012, 8, 4285–4292.
87. J. Řezáč and P. Hobza, J. Chem. Theory Comput., 2013, 9, 2151–2155.
88. R. Huenerbein, B. Schirmer, J. Moellmann and S. Grimme, Phys. Chem. Chem. Phys., 2010, 12, 6940–6948.
89. R. Sure, A. Hansen, P. Schwerdtfeger and S. Grimme, Phys. Chem. Chem. Phys., 2017, 19, 14296–14305.
90. L. Goerigk and R. Sharma, Can. J. Chem., 2016, 94, 1133–1143.
91. S. Kozuch and J. M. L. Martin, J. Chem. Theory Comput., 2013, 9, 1918–1931.
92. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
93. H. Kruse, A. Mladek, K. Gkionis, A. Hansen, S. Grimme and J. Sponer, J. Chem. Theory Comput., 2015, 11, 4972–4991.
94. E. Rossomme, L. A. Cunha, W. Li, K. Chen, A. R. McIsaac, T. Head-Gordon and M. Head-Gordon, J. Chem. Theory Comput., 2023, 19, 2827–2841.
95. P. Jurečka, J. Šponer, J. Černý and P. Hobza, Phys. Chem. Chem. Phys., 2006, 8, 1985–1993.
96. T. Takatani, E. G. Hohenstein, M. Malagoli, M. S. Marshall and C. D. Sherrill, J. Chem. Phys., 2010, 132, 144104.
97. J. Řezác, K. E. Riley and P. Hobza, J. Chem. Theory Comput., 2011, 7, 2427–2438.
98. J. Řezác, K. E. Riley and P. Hobza, J. Chem. Theory Comput., 2011, 7, 3466–3470.
99. O. A. Vydrov and T. Van Voorhis, Phys. Rev. Lett., 2009, 103, 063004.
100. N. Mardirossian and M. Head-Gordon, J. Phys. Chem., 2018, 148, 241736.
101. J.-D. Chai and M. Head-Gordon, J. Chem. Phys., 2009, 131, 174105.
102. N. Mardirossian and M. Head-Gordon, J. Chem. Phys., 2016, 144, 214110.
103. N. Mardirossian and M. Head-Gordon, Phys. Chem. Chem. Phys., 2014, 16, 9904–9924.
104. N. Mehta, M. Casanova-Páez and L. Goerigk, Phys. Chem. Chem. Phys., 2018, 20, 23175–23194.
105. G. Santra, M. Cho and J. M. L. Martin, J. Phys. Chem. A, 2021, 125, 4614–4627.
106. A. D. Becke, J. Chem. Phys., 2024, 160, 204118.
107. L. Wittmann, H. Neugebauer, S. Grimme and M. Bursch, J. Chem. Phys., 2023, 159, 224103.
108. T. Schwabe and S. Grimme, Phys. Chem. Chem. Phys., 2006, 8, 4398–4401.
109. A. Karton, A. Tarnopolsky, J.-F. Lamère, G. C. Schatz and J. M. L. Martin, J. Phys. Chem. A, 2008, 112, 12868–12886.
110. F. Yu, J. Phys. Chem. A, 2014, 118, 3175–3182.
111. A. Mohajeri and M. Alipour, J. Chem. Phys., 2012, 136, 124111.
112. L. M. Roch and K. K. Baldridge, Phys. Chem. Chem. Phys., 2018, 20, 4606.
113. L. Goerigk and S. Grimme, J. Chem. Theory Comput., 2011, 7, 291–309.
114. S. Kozuch and J. M. L. Martin, Phys. Chem. Chem. Phys., 2011, 13, 20104–20107.
115. S. Kozuch, D. Gruzman and J. M. L. Martin, J. Phys. Chem. C, 2010, 114, 20801–20808.
116. A. D. Becke, G. Santra and J. M. L. Martin, J. Chem. Phys., 2023, 158, 151103.
117. B. Nguyen, G. P. Chen, M. M. Agee, A. M. Burow, M. Tang and F. Furche, J. Chem. Theory Comput., 2020, 16, 2258–2273.
118. J. Liang and M. Head-Gordon, J. Chem. Theory Comput., 2025, 21, 12601–12621.
119. Y. A. Aoto, A. P. d Batista, A. Köhn and A. G. S. de Oliveira-Filho, J. Chem. Theory Comput., 2017, 13, 5291–5316.
120. A. Karton, S. Daon and J. M. L. Martin, Chem. Phys. Lett., 2011, 510, 165–178.
121. A. Wodyński and M. Kaupp, J. Chem. Theory Comput., 2022, 18, 6111–6123.
122. J. W. Hollet and P. M. W. Gill, J. Chem. Phys., 2011, 134, 114111.
123. M. E. Khatib, G. L. Bendazzoli, S. Evangelisti, W. Helal, T. Leininger, L. Tenti and C. Angeli, J. Phys. Chem. A, 2014, 118, 6664–6673.
124. J. M. Merritt, V. E. Bondybey and M. C. Heaven, Science, 2009, 324, 1548–1551.
125. B. Visser, M. Beck, P. Bornhauser, G. Knopp, J. A. V. Bokhoven, P. Radi, C. Gourlaouen and R. Marquardt, Mol. Phys., 2018, 117, 1645–1652.
126. P. Wang, S. Gong, Y. Li and Y. Mo, J. Chem. Phys., 2024, 160, 014304.
127. P. Wang, S. Gong and Y. Mo, J. Phys. Chem. Lett., 2024, 15, 12594–12600.
128. J. Csontos and M. Kállay, J. Phys. Chem. A, 2011, 115, 7765–7772.
129. S. Gong, P. Wang and Y. Mo, J. Phys. Chem. Lett., 2024, 15, 10842–10848.
130. A. G. Császár and T. Furtenbacher, J. Phys. Chem. A, 2015, 119, 10229–10240.
131. S. Gong, P. Wang and Y. Mo, Phys. Rev. A: At., Mol., Opt. Phys., 2023, 108, 042802.
132. D. Huang, F. Wang, J. Yang, Q. Cao and M. Wan, Spectrochim. Acta, Part A, 2014, 128, 163–167.

---