Local hybrid functionals are a more flexible class of density functional approximations, allowing for a position-dependent admixture of exact exchange. This additional flexibility, however, comes with a more involved mathematical form and a more complicated design. A common denominator for previously constructed local hybrid functionals is the usage of thermochemical benchmark data to construct these functionals. Herein, we design a local hybrid functional without relying on benchmark data. Instead, we construct it in a more ab initio manner, following the principles of modern meta-generalized gradient approximations and considering theoretical constraints. To achieve this, we make use of the density matrix expansion and a local mixing function based on an approximate correlation length. The accuracy of the developed density functional approximation is assessed for thermochemistry, excitation energies, polarizabilities, magnetizabilities, nuclear magnetic resonance (NMR) spin–spin coupling constants, NMR shieldings, and shifts, as well as EPR g-tensors and hyperfine coupling constants. Here, the new exchange functional shows a robust performance and is especially well suited for atomization energies, barrier heights, excitation energies, NMR coupling constants, and EPR properties, whereas it loses some ground for the NMR shifts. Therefore, the designed functional is a major step forward for functionals that have been designed from first principles.

The incorporation of exact exchange or Hartree–Fock (HF) exchange into semilocal density functional theory (DFT) was a milestone in computational chemistry.1–3 At the present time, DFT is probably the most widely applied computational method in chemistry and material science due to its favorable cost-accuracy ratio. Global hybrid functionals, such as B3LYP,4–6 PBE0,7,8 or TPSSh,9,10 use a fixed or static amount of HF exchange. Typically, global hybrid functionals use 5%–25% of HF exchange.3 While this leads to accurate results for many chemical properties, such a rigid admixture of exact exchange is not suited for charge-transfer excitations or the dissociation curve. Global hybrid functionals with the given amount of exact exchange still feature a self-interaction error. This is especially pronounced in one-electron regions. Therefore, range-separated hybrid (RSH) functionals, which typically use a significantly larger amount of exact exchange in the long-range region, were introduced.11–14 

Local hybrid functionals (LHFs) feature a position-dependent admixture of Hartree–Fock exchange via a local mixing function (LMF).15 Thus, the LMF increases the flexibility of the functional as the amount of HF exchange can be increased or decreased in certain regions and LHFs are a more general class of hybrid density functional approximations (DFAs). Compared to range-separated hybrid functionals, this functional form results in a smoother transition from spatial regions with a small amount of HF exchange to the regions with a large amount of HF exchange as it is possible to interpolate between the semilocal DFT and HF limits. However, LHFs introduce a gauge dependence of the exact-exchange energy density.16 Thus, a so-called calibration function16–18 (CF) is introduced in modern local hybrid functionals, such as LH14t-calPBE17 and LH20t.19 

The LMF is the key ingredient of LHFs, and different Ansätze were suggested.20,21 The first and still most prominent method to distinguish the electronic regions is the iso-orbital indicator (t-LMF).15 This allows us to identify one-electron regions; those accurate description requires an increased amount of exact exchange. Other approaches are based on the correlation length (z-LMF),22 reduced (spin) density gradients (s-LMF),23,24 the density overlap regions indicator (DORI-LMF),25,26 or the second-order Görling–Levy perturbation limit27 of the correlation (PSTS-LMF).28 We have recently carried out benchmark calculations of electric and magnetic response properties for LHFs based on different LMFs.29 This showed encouraging results for Johnson’s LHF22 based on the correlation length—especially for heavy elements. Compared to the PSTS-LMF, this LMF is also advantageous in numerically demanding cases and results in a smooth convergence behavior.30 However, only one functional based on this LMF was optimized,22 which consists of Becke’s 1988 exchange31 and correlation terms.32 Moreover, this functional does not make use of a calibration function, which is commonly included in the latest generation of LHFs based on the iso-orbital indicator. Therefore, we will first re-parameterize this functional with Becke’s 1995 correlation term33 and also add a CF.

In a second step, we will design a more sophisticated local hybrid functional from first principles. Modern (range-separated) hybrid and local hybrid functionals are typically optimized using (large) thermochemical test sets. While this may result in well-performing density functional approximations for the properties considered in the design, it may also come with a loss of generality and physical insight.34,35 In contrast, modern semilocal functionals36–39 are first designed to include (almost) all known theoretical constraints and then applied to large benchmark test sets. Herein, we will follow this route and adapt it to the form of local hybrid functionals.

In addition to electron correlation, an accurate description of the electronic structure throughout the Periodic Table of elements necessitates a proper treatment of relativistic effects.40–46 Herein, we will also generalize the two-component framework for open-shell systems developed in Ref. 47 to include the z-LMF22 and PSTS-LMF28 and higher-order derivatives to add the so-called calibration function of modern local hybrid functionals.19,21 It is shown that these functionals lead to substantially improved results for EPR properties, which are a problematic case for the first and second generation of t-LMF based functionals.47 

This paper is structured as follows. First, we will discuss the theoretical framework. In doing so, the optimization of a density functional approximation based on the correlation length from first principles is described. Then, the accuracy of the developed functionals is assessed for thermochemistry, excitations, polarizabilities, magnetizabilities, and NMR/EPR properties.

We will first review local hybrid functionals within an unrestricted Kohn–Sham (UKS) formalism in Subsection II A. Furthermore, the generalization to a two-component framework is described in Sec. II B. A local hybrid exchange functional is constructed from first principles in Subsection II C. In the supplementary material, we present an optimized LHF based on the correlation length similar to Johnson’s work.22 

In an UKS framework, the exchange–correlation (XC) energy is given by

EXCLHF=σ=α,βdr1aσ(r)eX,σDFT(r)+aσ(r)eX,σHF(r)+ECDFT(r),
(1)

where a is the LMF, eX,σDFT is the semilocal DFT exchange energy density, and eX,σHF denotes the exact-exchange or Hartree–Fock exchange energy density. ECDFT refers to the correlation energy. The LMF may depend on the spin quantities or total quantities. LMFs based on the latter include spin polarization and are called common LMFs.48 To carry out the integration in Eq. (1), Plessow and Weigend suggested to use a seminumerical scheme49 as this also computes the LMF, the exact-exchange energy contribution, and the semilocal DFT exchange energy density on a grid. The exact-exchange term is as follows:

eX,σHF(r)=12pqrsDpqσDrsσdrϕp*(r)ϕr(r)ϕq*(r)ϕs(r)|rr|,
(2)

where the integration with respect to r is evaluated. ϕp and Dpqσ denote the one-electron basis functions and the spin density matrix, respectively. The remaining integration with respect to r is performed numerically on a grid. Therefore, the position-dependent admixture of exact exchange is directly included in the DFT routines according to

EXLHF=σgwg1aσ(rg)eX,σHF(rg),
(3)

with g denoting a grid point. This LHF scheme was later applied self-consistently50 to various molecular properties, such as geometry gradients,29,51 excitation energies and polarizabilities,52–55 ionization potentials with the GW method,29,56 magnetizabilities,29 and nuclear magnetic resonance (NMR) parameters.29,30,57–60 The latter were also used to compute the magnetically induced current density and ring current strengths of aromatic systems.29 The reworked implementations of these properties were described in Refs. 29 and 55.

First, the LMF may be formulated based on the iso-orbital indicator.15 The respective t-LMF reads

aσ(r)=ctτσvWτσ=ct|ρσ|28ρστσ.
(4)

Here, the kinetic-energy density τ is compared to the von Weizsäcker approximation τvW.61 ρσ denotes the (spin) density. In the case of a common LMF, the total (kinetic-energy) densities are used. The prefactor ct in Eq. (4) is an optimized parameter.

Second, Johnson introduced a LMF using the correlation length z according to22 

aσ(r)=erf(czzσσ)
(5)

with the empirical parameter cz. Note that only the contribution of parallel spins is considered. Here, the correlation length is computed with the exchange potential UX and the exchange hole hX as

zσσ=cσσ|UX,σ|1+|UX,σ|1,
(6)
UX,σ=ds1s|hX,σ(r,s)|,
(7)

where cσσ are given by cαβ = 0.63 and cσσ = 0.88.32 Notably, this z-LMF explicitly depends on the underlying exchange functional.

A self-consistent treatment of spin–orbit coupling necessitates a generalized two-component Kohn–Sham formalism. This requires not only the spin excess density but also the complete spin magnetization vector, which is introduced according to

m=iφiσφi,
(8)

with the two-component spinor function φi and the vector σ containing the (2 × 2) Pauli spin matrices. The total electron density ρ and the non-collinear spin density ρs are given as62,63

ρ=iφiφi,
(9)
ρs=mm1/2.
(10)

In a basis set representation, the non-collinear exchange–correlation energy depends on the total density matrix M0 and the three spin vector density matrices Mi with i = x, y, z. These are defined as62,63

M0=Re(Pαα)+Re(Pββ),
(11)
Mx=Re(Pαβ)+Re(Pβα),
(12)
My=Im(Pαβ)Im(Pβα),
(13)
Mz=Re(Pαα)Re(Pββ).
(14)

The XC potential follows as

VμνXC[M0(r),M(r)]=δEXC[M0(r),M(r)]δM0,μν(r)+σδEXC[M0(r),M(r)]δMμν(r).
(15)

The generalization of a non-relativistic spin-density functional is done using the spin-up and spin-down densities,64,65

ρ=ρ+ρs/2,
(16)
ρ=ρρs/2.
(17)

Only considering M0 and Mz results in the UKS limit. Therefore, an unrestricted Kohn–Sham implementation of the XC potential and the XC energy can be extended straightforwardly.64 This involves the following major steps. First, the electron (spin) density and its derivatives, such as the gradient and the Laplacian or the kinetic-energy density, are evaluated at a given grid point. Second, the respective derivatives are multiplied with its spin density matrix contribution, and the sum of all three vector components is formed. The inverse of the total spin density is used as a prefactor. Third, the spin-up and spin-down contributions are constructed. Then, the exchange and correlation functional expressions can be evaluated similar to UKS. For further details, we refer to Refs. 64 and 66. Herein, we have implemented this scheme up to the Hessian of the electron density and the gradient of the kinetic-energy density to evaluate all calibration functions described in Ref. 18. A generalization of spin-dependent local mixing function is more involved, and many modern LHFs based on the t-LMF use a common local mixing function. Therefore, the implementation of the 2c t-LMF is still restricted to the common variant. Thus, LH07t-SVWN67 or LH14t-calPBE17 is not yet available. Moreover, current-dependent terms arise as spin–orbit coupling is closely related to magnetic induction, and these are neglected herein;68 however, their implementation will be presented elsewhere.

A common denominator for previously constructed local hybrid functionals is the usage of thermochemical benchmark data to construct these functionals. In this paper, we aim at constructing a local hybrid exchange functional without relying on benchmark data. Instead, we construct it in a more ab initio manner, following the principles of previous meta-GGA functionals.9,36,37

We define the local exchange part eX as usual using an enhancement factor FX,

EX=σ=α,βdrFX(ρσ,ρσ,τσ;r)eXunif(ρσ;r),
(18)

where the exchange energy density from the local (spin) density approximation is defined as

eXLSDA(ρσ;r)=ρσeXunif(ρσ;r)=3ρσ4π(3π2ρσ)1/3.
(19)

As in many recent density functional approximations (DFAs), we chose a semilocal enhancement factor, where FX is a function of the density ρσ, the gradient of the density ∇ ρσ, and the kinetic energy density τσ. As the exchange enhancement factor only depends on a single spin coordinate, it is dropped in the following equations. Particularly, we adopt the density matrix expansion (DME) approach of Tao and Mo,36 which has been shown to perform well for solids.69–71 Within the DME of Tao and Mo, the enhancement factor FXDME is given as

FXDME=1fX2+7RX9fX4.
(20)

The auxiliary quantities RX and fX are defined as

RX=1+59454yτ3(λX2λX+0.5)×ττunif|ρ|272ρ1τunif
(21)

and

fX=1+1070y27+βXy21/10
(22)

with τunif=(3/10)(3π2)2/3ρ5/3. λ is a real number between 0.5 and 1, describing the coordinate transformation. λ = 1 corresponds to the conventional exchange hole.36 Furthermore, y = (2λ − 1)2p is a scaled version of the reduced density gradient s depending on the coordinate transformation with p = s2. For a local hybrid functional, exact and local exchange are combined. Therefore, we choose λX = 1, corresponding to the untransformed exchange hole. The parameter βX will be determined later. As outlined by Tao and Mo, the DME is not exact for slowly varying densities. Hence, it is advisable to pair it with a fourth-order gradient correction.36 We use the same expressions as Tao and Mo for the correction in the slowly varying limit,

FXSC=1+10(10/81+50p/729)p+146q̃2/2025(73q̃/405)[3τvW/(5τ)](1τvW/τ)1/10
(23)

with q̃=(9/20)(α1)+2p/3 and α=ττvW/τunif.36 The interpolation between FXDME and FXSC is also unaltered from the Tao–Mo functional, yielding the final enhancement factor FX as

FX=wFXDME+(1w)FXSC
(24)

with the interpolation function w=[(τvW/τ)2+3(τvW/τ)3]/[1+(τvW/τ)3]2.

As missing piece to construct a local hybrid functional, a suitable mixing function is needed. The task of the mixing function is to augment the local exchange from the DME with exact exchange. We adopt an approach similar to the correlation length as suggested by Johnson.22 Using an approximated correlation length zDME,

zσσDME=(|Uσ|1+|Uσ|1),
(25)

we define U as

Uσ=cF[(1+ζ)ρσ]1/31fL2+7RL9fL4
(26)

with cF = 3/8 · 42/3(3/π)1/3 and ζ = (ρσρσ)/(ρσ + ρσ). Uσ is obtained by reversing the spin indices. The expressions of RL and fL are equivalent to those in Eqs. (21) and (22). The different subscripts hint at the parameters βL and λL possibly being different from those used in Eqs. (21) and (22). To map zσσDME to the interval {0, 1}, we define the local mixing function aDME as

aDME=1exp(cLzσσDME).
(27)

As stated in Ref. 47, it is advantageous to work with “common” mixing functions for 2c calculations using a single mixing function for both spin channels. Furthermore, common mixing functions have also shown to perform exceptionally well for predicting excited triplet states.72 Therefore, we exclusively use the common LMF for the designed functional in this work. Obtaining a common mixing function based on Eq. (27) is straightforward. Simply τσ is replaced by τσ + τσ and |∇ρσσ|/ρσ is replaced by |∇ρσσ + 2∇ρσσ + ∇ρσσ|/(ρσ + ρσ) during the evaluation of zσσDME.

What remains to be done is determining the parameters cL, λX/L, and βX/L for the exchange functional and the LMF. At first, a value of λX = 1.0 is set for the exchange enhancement factor in Eq. (20). This relates to the choice of parameterizing the untransformed (nonlocal) exchange hole. This is a natural choice as the exact exchange entering Eq. (1) is also calculated from the conventional hole. For the remaining parameters, two Ansätze can be chosen to achieve a set of parameters. In the first Ansatz, different parameters are chosen for βX/L and λX/L in the exchange enhancement factor and the LMF, leading to a total of five parameters to be optimized. In the second Ansatz, the conditions βX = βL and λX = λL are required, leading to only three parameters that need to be optimized.

Beginning with the first Ansatz, the parameters βL = 79.873 and λL = 0.6866 are set to the values obtained for the hydrogen atom by Tao and Mo.36 Therefore, the correlation length is measured for the transformed (localized) exchange hole. The localized exchange hole can be expected to be more appropriate for (semi-)local exchange,36 yielding a more reliable description of the switching parameter. What is left to be determined are the parameters βX from the DME, parameterizing the untransformed exchange hole, and cL, accounting for the inclusion of exact exchange through the LMF. The latter two parameters are obtained by fitting them to spin-unpolarized two-electron densities.73,74 As our aim is to design a functional without any prior exact knowledge of total exchange–correlation energies, even atomic ones, we fit the parameters to theoretical considerations. First, in the low density limit, which challenges the Lieb–Oxford bound more, the energy is constrained by eX=1.174eXLSDA.74 Second, in the high density limit, where exchange dominates, we expect the results to be closer to the mean-field Hartree–Fock solution, i.e., eX=1.16588eXLSDA.75 Both of these values can be derived analytically. To guide the optimization, we choose two two-electron systems as for those, Hartree–Fock exchange provides solutions close to the exact solution.73 For the low-density limit, He is used, while we chose Hg78+ as guide for the high-density limit. For both atoms, we evaluate the corresponding local spin density approximation (LSDA) exchange energies at the Hartree–Fock solution using saturated basis sets. For He and Hg78+, we find LSDA energies of −0.885 and −42.699 hartree. Both of these values are within 2 mhartree of those obtained for the exact density.73 Multiplying with the analytic prefactors yields −1.039 and −49.781 hartree, respectively. Finally, we optimize βX for a given value of cL, yielding possible pairs of solutions. We chose the pair with the highest value for cL, yielding cL = 0.18 and βX = 265.25.

In the second Ansatz, only two parameters need to be optimized. As one-electron densities are exact by construction of our exchange approximation, we drop the hydrogen atom as the reference system. For the three parameter model, we, therefore, only take into account the low- and high-density limits of spin-unpolarized two-electron systems. Repeating the same procedure as above yields βX = 262.25 and cL = 0.215. Indeed, βX is only a very weak function of the optimization procedure, and we find hardly any difference compared to the five parameter fit. cL = 0.215 hints at significantly more exact exchange being incorporated when a three parameter fit is desired. It shall be noted that while the Lieb–Oxford bound was taken into account during the evaluation of the parameters, the derived exchange approximation may still violate it for certain systems and basis sets as it is not strictly constrained by it. Furthermore, our functional certainly violates the conjectured local version of the Lieb–Oxford bound. The latter issue could be remedied by a gauge transformation.

For convenience, the parameters of the developed functional, which we will be refer to as Tao-Mo-Holzer-Franzke (TMHF) (five parameter fit) and TMHF-3P (three parameter fit), are summarized in Table I. We pair these local hybrid exchange functionals with the B95 correlation functional33 as the latter has shown to work well in conjunction with local hybrid exchange functionals.19 We note that the B95 correlation functional uses two parameters optimized for the He and Ne atoms.33 Attempts to employ the (modified) TPSS correlation9,36 failed as this yields very inaccurate results for thermochemical properties. The TPSS correlation is more suited for functionals with a small amount of exact exchange.10 Viable alternatives to B95 correlation may be obtained from the rung-3.5 Ansatz of Ramos and Janesko76 or from the recent non-dynamical correlation Ansatz of Becke.35 

TABLE I.

Parameters of the five parameter TMHF and three parameter TMHF-3P fits. cL, βL, and λL are employed for the LMF, while βX and λX are used for the exchange density. Note that all parameters were derived from physical constraints.

βLλLcLβXλX
TMHF 79.873 0.6866 0.180 265.25 1.0 
TMHF-3P βX λX 0.215 265.25 1.0 
βLλLcLβXλX
TMHF 79.873 0.6866 0.180 265.25 1.0 
TMHF-3P βX λX 0.215 265.25 1.0 

Figure 1 shows a comparison of the LMF defined by Johnson22 outlined in Eqs. (5) and (6) and our DME-based one outlined in Eqs. (27) and (25) for four diatomic molecules. In both approaches, the mixing fraction approaches zero near the nuclei. Furthermore, their behavior in the bonding region is similar, incorporating ∼20% exact exchange. Larger differences can be found in the tail regions of the density, where the correlation length defined in Eq. (5) grows approximately linear. Contrarily, zσσDME used in TMHF grows nearly exponentially, as outlined in Fig. 1(a).

FIG. 1.

(a) zσσ obtained from Eq. (25) (straight lines, DME) and Eq. (6) (dashed lines, Johnson). (b) Local mixing function a obtained from Eq. (27) (straight lines, DME) and Eq. (5) (dashed lines, Johnson). Properties calculated for four molecules at self-consistent TMHF/aug-cc-pVQZ orbitals and plotted along the internuclear axis.

FIG. 1.

(a) zσσ obtained from Eq. (25) (straight lines, DME) and Eq. (6) (dashed lines, Johnson). (b) Local mixing function a obtained from Eq. (27) (straight lines, DME) and Eq. (5) (dashed lines, Johnson). Properties calculated for four molecules at self-consistent TMHF/aug-cc-pVQZ orbitals and plotted along the internuclear axis.

Close modal

For another interesting example, the Li dimer at 5.0 and at 10.0 bohrs, the DME local mixing function is outlined in Fig. 2. It clearly shows the different behavior between bonded and stretched Li2. In the former case, the LMF has a local minimum at the center, while in the latter case, a local maximum is found. Comparably, at the bond center, the amount of exact exchange included is raised by roughly 25% in the stretched dimer. This effect will lead to significantly improved barrier heights while not degrading the overall performance for other properties. In Ref. 77, it was noted that due to the order-of-limits problem of the interpolation function, stretched Li2 will not converge with the Tao–Mo functional using the aug-cc-pVQZ basis set,78–80 which was confirmed by us. The root of this issue is a singularity located at the center between the two Li atoms, caused by the order-of-limits problem.77 TMHF, however, converges smoothly for the stretched Li2. This can be attributed to the large amount of exact exchange incorporated at the center of the stretched bond, as shown in Fig. 2, and the increased numerical rigor of our implementation. Still, the local exchange part of TMHF is plagued by the order-of-limits problem, but for stretched bonds, it is expected to be less harmful. Near the nucleus, where the LMF approaches zero, the singularity from the order-of-limits problem is, however, not countered. Therefore, no improvements of TMHF over the initial Tao–Mo functional (or other standard meta-GGA functionals) are expected for properties that are sensitive to the exchange contribution and the electron density in the vicinity of the nuclei.

FIG. 2.

DME local mixing function as a function of z at various distances between two Li atoms, calculated using the TMHF/aug-cc-pVQZ method. The distance between the two Li atoms is r = 5.0 bohrs (solid line) and r = 10.0 bohrs (dashed line).

FIG. 2.

DME local mixing function as a function of z at various distances between two Li atoms, calculated using the TMHF/aug-cc-pVQZ method. The distance between the two Li atoms is r = 5.0 bohrs (solid line) and r = 10.0 bohrs (dashed line).

Close modal

The outlined functionals were implemented into the TURBOMOLE program suite.81–84 The mathematical functions for the LMFs were generated using Maple scripts,85 while the individual exchange and correlation functional terms are computed with Libxc.86–88 For convenience, the Maple files of the TMHF local mixing function (separate spin-channels and common spin-channels) are part of the supplementary material and can be included in the routines for LHFs. Note that we exclusively use the common variant in this work.

The existing Kramers unrestricted two-component implementation47 was reworked for the efficiency and extended to include higher-order derivatives of the density, i.e., the Laplacian, Hessian, and the gradient of τ. This allows us to use the pig1, pig2, and tpig1 calibration function18 and general meta-GGAs for LHFs. Therefore, the LH20t functional19 is also now available in two-component open-shell calculations. We further added interfaces to Libxc86–88 to support (almost) all exchange and correlation functional ingredients. The revised thresholds of Ref. 30 are used for the numerical integration.89,90 We note that tighter thresholds and large grids are helpful for the SCF convergence and the numerical accuracy of the B95 correlation contribution. The efficiency was increased by using the routines developed in Refs. 29 and 55, i.e., the screening procedure and memory handling described therein are applied in the general two-component case.

We note that the reworked LHF gradient routines evaluate the high-angular momentum contributions similar to the integral routines developed for the finite nucleus model in relativistic all-electron calculations.91 Parallelization is available throughout with the OpenMP scheme.92,93

In this section, the performance of the newly developed five parameter TMHF and three parameter TMHF-3P functionals will be investigated. In addition, we will investigate a re-parameterized version of the z-LMF-based functional of Johnson,22 labeled LHJ-HF and LHJ-HFcal. Details about LHJ-HF(cal) can be found in the supplementary material, Sec. S1. We will limit the assessment of accuracy to thermochemistry, excitation energies, and EPR properties in the main text. Further studies on magnetizabilities, polarizabilities, NMR spin–spin coupling constants, and NMR shielding constants and shifts are presented in the supplementary material. The respective computational details are also given therein. To provide a comprehensive overview, we compare to a variety of density functionals commonly used at the present time and earlier designed local hybrid functionals that have been fitted to thermochemical data. This includes the Perdew–Burke–Ernzerhof generalized gradient approximation (PBE GGA);7 the meta-GGA TPSS;9 Tao–Mo;36 SCAN;37 the hybrid functionals PBE0,7,8 B3LYP,4–6 and TPSSh;9,10 the range-separated functionals CAM-B3LYP,14 LC-ωPBE,94 and ωB97X-D;95 and the local hybrid functionals LH12ct-SsirPW92,48 LH14t-calPBE,17 LH20t,19 LHJ14,22 and mPSTS-a1.28 

Atomization energies and barrier heights are important quality measures for DFAs. They yield a first overview of the general quality of functionals and are themselves commonly used when newly parameterized functionals are designed. Atomization energies were assessed for the W4-11 test set,96 and barrier heights were assessed for the BH76 test set.97–99 Both of those sets are subsets of the extensive “general main group thermochemistry, kinetics, and noncovalent interactions” (GMTKN) set.100 To yield values that can be directly compared to the GMTKN values, the def2-QZVP basis set101 was used throughout, in conjunction with a large integration grid (grid 4) for numerical integration.89,90 Results for other functionals are taken from Ref. 100.

For the excitation energies, we consider the benchmark test set of Suellen et al. with experimental reference results, which match third-order coupled cluster (CC3) values to within 0.05 eV or less.102 Structures were taken from Ref. 102. In line with this reference, we use the aug-cc-pVTZ basis set78–80 and large grids (grid size 4) for the numerical integration of the DFT parts.89,90 Tight SCF thresholds of 10−9 hartree for the energy and 10−7 a.u. for the root mean square of the density matrix are applied. Response equations are converged with a threshold of 10−5 a.u. for the norm of the residuum.103 The excitation energies are corrected by the B3LYP zero-point vibrational energies.104 Results with conventional functionals are taken from Ref. 102, while the results with all previously designed local hybrids are taken from Ref. 29.

EPR properties, such as the hyperfine coupling (HFC) constant and the g-tensor, are challenging properties for local hybrid functionals due to the high-density limit.47 We have recently benchmarked EPR properties in a self-consistent spin–orbit exact two-component (X2C) framework.105,106 Herein, we extent these studies to local hybrid functionals, and those general two-component implementation is described herein. To do so, we consider the 17 small transition-metal complexes of Ref. 107, namely, [MoNCl4]2−, [MoOF4], [MoOCl4], [MoOF5]2−, [MoOBr5]2−, [WOCl4], [WOF5]2−, [WOBr5]2−, [TcNF4], [TcNCl4], [TcNBr4], [ReNF4], [ReNCl4], [ReNBr4], [ReOBr4], [ReOF5], and [OsOF5]. We use the same computational parameters as in Refs. 105 and 106. In detail, the x2c-QZVPall-2c basis set108 is applied for all elements. Large grids (grid 5a) are chosen for the numerical integration.89,90,109 The conductor-like screening model (COSMO) is applied with the default parameters to compensate the negative charge.110,111 We use the X2C Hamiltonian in the diagonal local approximation to the unitary transformation (DLU).30,91,105,106,112–115 The restricted kinetic balance (RKB) condition116 is employed for the HFC, whereas both RKB and the restricted magnetic balance (RMB) conditions117 are employed for the g-tensor. An SCF threshold of 10−9 hartree is applied. In addition to the functionals used in Secs. III A and B and the Hartree–Fock method, we consider the following DFAs additionally to provide a more complete overview in this chapter as data are yet comparably rare: S-VWN,118–120 KT3,121 BP86,31,122 revTPSS,123,124 r2SCAN,38,39 BH&HLYP,5,31,125 PBE0, including 40% of HF exchange (PBE0-40HF),7,8,107 B97,126 B97-2,127 revTPSSh,123,124 TPSS0,10,128 r2SCANh,38,39,129 r2SCAN0,38,39,129 r2SCAN50,38,39,129 CAM-QPT-00,130 CAM-QTP-02,131 HSE06,132–134 LH12ct-SsifPW92,48 LH20t* (LH20t without calibration function),19 and mPSTS-noa2.28,29 Note that SCAN has been replaced by r2SCAN as HFC constants and g-tensors are considerably more sensitive to the integration grid. In the main text, results for the same functionals as above are shown. The complete results are listed in the supplementary material.

To assess the basic properties of the newly constructed functionals, their thermochemical properties are investigated. The W4-1196 and BH76 subsets97–99 from the GMTKN55 test set100 were chosen, outlining the principal capabilities of a density functional approximation (DFA) for the calculation of atomization energies (W4-11) and barrier heights (BH76). The results are displayed in Figs. 3 and 4.

FIG. 3.

Mean standard deviation (MSD), mean average deviation (MAD), and root mean square deviation (RMSD) for the atomization energies of the W4-11 test set. All values are in kcal/mol. PBE and LHJ14 are omitted due to errors of more than 10 kcal/mol.

FIG. 3.

Mean standard deviation (MSD), mean average deviation (MAD), and root mean square deviation (RMSD) for the atomization energies of the W4-11 test set. All values are in kcal/mol. PBE and LHJ14 are omitted due to errors of more than 10 kcal/mol.

Close modal
FIG. 4.

Mean standard deviation (MSD), mean average deviation (MAD), and root mean square deviation (RMSD) for the barrier heights of the BH76 test set. All values are in kcal/mol.

FIG. 4.

Mean standard deviation (MSD), mean average deviation (MAD), and root mean square deviation (RMSD) for the barrier heights of the BH76 test set. All values are in kcal/mol.

Close modal

TMHF is a major step forward for functionals that have been designed from first principles. With MAE/RMSD values of 2.78/4.60 kcal/mol for the W4-11 atomization energy subset and 2.80/3.36 kcal/mol for the barrier heights of the BH76 subset, it outperforms any other functionals that have been constructed from first principles. Popular thermochemically optimized functionals such as ωB97X-D exhibit similar errors for the W4-11 subset. Furthermore, previously reported thermochemically optimized local hybrid functionals, for example, LH20t, are not able to outperform TMHF. From a viewpoint of thermochemistry, TMHF and TMHF-3P perform as well as the best empirically parameterized functionals. While for barrier heights parameterized functionals still hold a slight edge, especially the slightly higher amount of exact exchange incorporated in TMHF-3P yields excellent results. The remaining differences in barrier heights are rather small and probably do not outweigh the loss of generality.

Table II summarizes the errors of thermochemical properties for a set of DFAs that has been derived from first principles. Here, a clear trend is observed following Jacob’s ladder. The most pronounced error reduction happens from GGA to meta-GGAs, followed by the step to global hybrids. An obvious exception is the meta-GGA SCAN, which is able to perform as well as global hybrids for atomization energies, still losing out for barrier heights. Barrier heights are more sensitive to the inclusion of exact exchange, generally preferring functionals that incorporate a higher amount of exact exchange. The non-empirical local hybrids presented in this work, TMHF and TMHF-3P, again lower the error bar significantly. The statistical errors in barrier heights are nearly halved. Atomization energies are also improved by more than 1 kcal/mol on average, yet yielding a less pronounced underbinding when compared to the global hybrids TPSSh and PBE0.

TABLE II.

Mean standard deviation (MSD), mean average deviation (MAD), and root mean square deviation (RMSD) for the atomization energies of the W4-11 test set and the barrier heights of the BH76 test set with functionals derived from first principles. All values are in kcal/mol.

W4-11BH76
MSDMADRMSDMSDMADRMSD
PBE 13.35 14.96 18.50 −9.11 9.15 10.39 
Tao–Mo 3.27 7.45 9.66 −8.21 8.24 9.24 
TPSS 3.27 5.11 6.65 −8.61 8.63 9.58 
TPSSh −1.62 4.41 6.12 −6.65 6.68 7.48 
SCAN −0.17 4.01 5.72 −7.36 7.66 8.37 
PBE0 −1.75 3.62 5.73 −3.17 4.62 5.90 
TMHF-3P −0.27 3.21 4.97 −1.83 2.09 2.59 
TMHF −0.73 2.78 4.60 −2.71 2.80 3.36 
W4-11BH76
MSDMADRMSDMSDMADRMSD
PBE 13.35 14.96 18.50 −9.11 9.15 10.39 
Tao–Mo 3.27 7.45 9.66 −8.21 8.24 9.24 
TPSS 3.27 5.11 6.65 −8.61 8.63 9.58 
TPSSh −1.62 4.41 6.12 −6.65 6.68 7.48 
SCAN −0.17 4.01 5.72 −7.36 7.66 8.37 
PBE0 −1.75 3.62 5.73 −3.17 4.62 5.90 
TMHF-3P −0.27 3.21 4.97 −1.83 2.09 2.59 
TMHF −0.73 2.78 4.60 −2.71 2.80 3.36 

The test set of Suellen et al., which is composed of 41 excitation energies with accurate experimental and approximate coupled cluster singles, doubles, and triples (CC3) reference values has become a popular way of benchmarking the capabilities to predict excitation energies.102 For these molecules, TMHF and TMHF-3P perform very well too, as shown in Fig. 5, with TMHF having the edge.

FIG. 5.

Mean standard deviation (MSD), mean average deviation (MAD), standard deviation (STD), root mean square deviation (RMSD), and maximum error (Max. Err.) for 41 excitation energies of small molecules. All values are in eV.

FIG. 5.

Mean standard deviation (MSD), mean average deviation (MAD), standard deviation (STD), root mean square deviation (RMSD), and maximum error (Max. Err.) for 41 excitation energies of small molecules. All values are in eV.

Close modal

Most obviously, the parameterized functionals, which performed well for thermochemistry, are no longer the top contenders. TMHF easily outperforms most other DFAs, although most statistical values are rather close for the top performers here. While performing very well overall, TMHF also features the by far lowest maximum error of all functional approximations that have been tested.102 A maximum error of 0.46 eV equals a reduction of 0.1–0.3 eV compared to other popular DFAs. Even compared to coupled-cluster singles and doubles (CCSD), which exhibits a maximum error of 0.43 eV, it remains competitive. Common to most local hybrid functionals, the average deviation of TMHF is again centered around 0 eV. This outlines the balanced interpolation between local and exact exchange of the newly constructed TMHF local hybrid. As outlined in the supplementary material, charge-transfer excitations are also well captured by TMHF and TMHF-3P, yielding errors comparable to range-separated hybrid functionals in these cases.

The results for the HFC constant and the g-tensor are shown in Figs. 6 and 7, respectively. Here, the statistical evaluation of the g-tensor is shown with the RMB condition. The previous functional generation with LH12ct-SsirPW92 and LHJ14 does not yield accurate results. These functionals exhibit larger errors than some of the conventional hybrid functionals. The mPSTS-a1 functional performs similar to established global hybrid functionals, such as B3LYP and TPSSh. However, it is outperformed by many range-separated hybrid functionals and global hybrid functionals with an increased amount of exact exchange such as TPSS0 and PBE0-40HF (see the supplementary material). Among the conventional hybrids, the latter shows the smallest errors with mean absolute percent-wise deviations (MAPDs) of 5%–6% for the HFC constant and also yield smaller errors for the g-tensor (MAPDs of 13%–18%).

FIG. 6.

Assessment of various density functional approximations for the HFC constant compared to the experimental findings for a set of 12 transition-metal complexes.107 MAPD and STD denote the mean absolute percent-wise deviation and its standard deviation. The data of the conventional functionals are taken from Ref. 105.

FIG. 6.

Assessment of various density functional approximations for the HFC constant compared to the experimental findings for a set of 12 transition-metal complexes.107 MAPD and STD denote the mean absolute percent-wise deviation and its standard deviation. The data of the conventional functionals are taken from Ref. 105.

Close modal
FIG. 7.

Assessment of various density functional approximations for the Δg-shift compared to the experimental findings for the set of 17 transition-metal complexes.107 MAPD and STD denote the mean absolute percent-wise deviation and its standard deviation. The data of the conventional functionals are taken from Ref. 106. [TcNCl4] and [ReNBr4] are neglected in the statistical evaluation.

FIG. 7.

Assessment of various density functional approximations for the Δg-shift compared to the experimental findings for the set of 17 transition-metal complexes.107 MAPD and STD denote the mean absolute percent-wise deviation and its standard deviation. The data of the conventional functionals are taken from Ref. 106. [TcNCl4] and [ReNBr4] are neglected in the statistical evaluation.

Close modal

TMHF performs remarkably well with an error of 4.77% for the HFC constant and 11.93% for the g-tensor. Thus, TMHF outperforms all other density functional approximations. TMHF-3P is again very similar to TMHF and only slightly inferior. This is also a remarkable improvement over the parent functional of Tao and Mo, which results in errors of 29.73% and 34.23%, respectively. The LH20t functional and the newly developed LHJ-HF functionals deliver similar results and are among the top performers with errors of about 7% for the HFC, while LH20t loses ground for the g-tensor and shows slightly larger errors than LH12ct-SsirPW92. For LHJ-HF, the improvements can be mainly attributed to the re-optimized parameter for the admixture of HF exchange as EPR properties are sensitive toward the amount of exact exchange.105,106

The impact of the calibration function on the EPR properties is negligible, whereas it may lead to an unfavorable SCF convergence behavior due to the increased numerical difficulties of the higher-order derivatives. This is especially important for two-component calculations, which require re-optimized thresholds for the numerical integration.30,109 Thus, we recommend functionals without the calibration function for two-component calculations of open-shell compounds.

Taking together, local hybrid functionals are generally able to deliver accurate results for EPR parameters. Overall, TMHF outperforms all other functionals. The LH20t and LHJ-HF families lead to similar errors for the HFC constant, and LHJ-HF shows improvements for the g-tensor. For both the HFC and the g-tensor, TMHF and LHJ-HF are the top performers. To illustrate the performance of local hybrid functionals for systems with more than one unpaired electron, we list the hyperfine coupling constants for the Lanthanide complex [TbPc2] in the supplementary material. For this example, any density functional approximation other than TMHF yields large errors. TMHF is, however, able to predict the hyperfine coupling constants for [TbPc2] with very good accuracy, stressing that the results shown in Fig. 6 are transferable.

We derived a new local hybrid exchange approximation, termed TMHF and TMHF-3P, from first principles. For their construction, we only take into account the low-density and high-density limits of the exchange energies of two-electron systems, and for TMHF, we additionally take the exact solution of the hydrogen atom. The derived functionals are, therefore, not fitted to any bound systems or reaction energies. Statistical errors of thermochemical properties reveal that, indeed, the TMHF exchange models are a significant step forward for density functional approximations designed from first principles. Being the next step on rung 4 of Jacob’s ladder,135 TMHF significantly outperforms all previously presented approximations constructed from first principles.

In the subsequent investigations of various properties, the assumption that density functional approximations from first principles are generally transferable could be verified. For various properties, such as the calculation of excited states, EPR parameters, or NMR coupling constants, TMHF is the best in class. It does provide not only significantly lower statistical errors but also far lower maximum errors as compared to other leading functionals. This leads to more reliable predictions across different molecules and properties, which may have very different needs. There are, of course, still points where future work is desperately needed. For example, a rather simple local correlation model is used, not being able to truly describe multi-configurational settings. A fully non-empirical correlation treatment could be derived or paired with the outlined exchange approach. Furthermore, core-related properties such as NMR shifts are not improved, given the lack of exact exchange near the nucleus from our model. Finally, van der Waals interactions still need to be included using a dispersion correction model as long range correlation can also not be modeled by our approach alone.

Despite the remaining deficiencies, we conclude that certainly TMHF is strikingly close to a one-for-all functional for the time being.

See the supplementary material for the re-parameterized version of Johnson’s local hybrid functional, individual results and statistical evaluation for magnetizabilities and polarizabilities, individual results and statistical evaluation for NMR spin–spin coupling constants, individual results and statistical evaluation for NMR shieldings and shifts, individual results for thermochemistry, individual results and statistical evaluation for excitation energies, results of TMHF and TMHF-3P for the QUEST#6 excitation energy test set of Ref. 136, individual results and statistical evaluation for hyperfine coupling constants and g-tensors, application to the isotropic EPR hyperfine coupling constant of [TbPc2], Maple files of TMHF to allow for an easy incorporation into quantum chemical programs, and Cartesian coordinates of optimized structures for NMR couplings and NMR shifts.

C.H. acknowledges funding from the Volkswagen Foundation. Y.J.F. acknowledges the financial support from TURBOMOLE GmbH.

The authors have no conflicts to disclose.

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

1.
G. I.
Csonka
,
J. P.
Perdew
, and
A.
Ruzsinszky
,
J. Chem. Theory Comput.
6
,
3688
(
2010
).
2.
A. D.
Becke
,
J. Chem. Phys.
140
,
18A301
(
2014
).
3.
N.
Mardirossian
and
M.
Head-Gordon
,
Mol. Phys.
115
,
2315
(
2017
).
4.
A. D.
Becke
,
J. Chem. Phys.
98
,
5648
(
1993
).
5.
C.
Lee
,
W.
Yang
, and
R. G.
Parr
,
Phys. Rev. B
37
,
785
(
1988
).
6.
P. J.
Stephens
,
F. J.
Devlin
,
C. F.
Chabalowski
, and
M. J.
Frisch
,
J. Phys. Chem.
98
,
11623
(
1994
).
7.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
77
,
3865
(
1996
).
8.
C.
Adamo
and
V.
Barone
,
J. Chem. Phys.
110
,
6158
(
1999
).
9.
J.
Tao
,
J. P.
Perdew
,
V. N.
Staroverov
, and
G. E.
Scuseria
,
Phys. Rev. Lett.
91
,
146401
(
2003
).
10.
V. N.
Staroverov
,
G. E.
Scuseria
,
J.
Tao
, and
J. P.
Perdew
,
J. Chem. Phys.
119
,
12129
(
2003
).
11.
P. M. W.
Gill
,
R. D.
Adamson
, and
J. A.
Pople
,
Mol. Phys.
88
,
1005
(
1996
).
12.
T.
Leininger
,
H.
Stoll
,
H.-J.
Werner
, and
A.
Savin
,
Chem. Phys. Lett.
275
,
151
(
1997
).
13.
Y.
Tawada
,
T.
Tsuneda
,
S.
Yanagisawa
,
T.
Yanai
, and
K.
Hirao
,
J. Chem. Phys.
120
,
8425
(
2004
).
14.
T.
Yanai
,
D. P.
Tew
, and
N. C.
Handy
,
Chem. Phys. Lett.
393
,
51
(
2004
).
15.
J.
Jaramillo
,
G. E.
Scuseria
, and
M.
Ernzerhof
,
J. Chem. Phys.
118
,
1068
(
2003
).
16.
J.
Tao
,
V. N.
Staroverov
,
G. E.
Scuseria
, and
J. P.
Perdew
,
Phys. Rev. A
77
,
012509
(
2008
).
17.
A. V.
Arbuznikov
and
M.
Kaupp
,
J. Chem. Phys.
141
,
204101
(
2014
).
18.
T. M.
Maier
,
M.
Haasler
,
A. V.
Arbuznikov
, and
M.
Kaupp
,
Phys. Chem. Chem. Phys.
18
,
21133
(
2016
).
19.
M.
Haasler
,
T. M.
Maier
,
R.
Grotjahn
,
S.
Gückel
,
A. V.
Arbuznikov
, and
M.
Kaupp
,
J. Chem. Theory Comput.
16
,
5645
(
2020
).
20.
T. M.
Henderson
,
B. G.
Janesko
, and
G. E.
Scuseria
,
J. Phys. Chem. A
112
,
12530
(
2008
).
21.
T. M.
Maier
,
A. V.
Arbuznikov
, and
M.
Kaupp
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
9
,
e1378
(
2018
).
22.
E. R.
Johnson
,
J. Chem. Phys.
141
,
124120
(
2014
).
23.
A. V.
Arbuznikov
and
M.
Kaupp
,
Chem. Phys. Lett.
440
,
160
(
2007
).
24.
T.
Schmidt
,
E.
Kraisler
,
A.
Makmal
,
L.
Kronik
, and
S.
Kümmel
,
J. Chem. Phys.
140
,
18A510
(
2014
).
25.
P.
de Silva
and
C.
Corminboeuf
,
J. Chem. Theory Comput.
10
,
3745
(
2014
).
26.
P.
de Silva
and
C.
Corminboeuf
,
J. Chem. Phys.
142
,
074112
(
2015
).
27.
A.
Görling
and
M.
Levy
,
Phys. Rev. A
50
,
196
(
1994
).
28.
J. P.
Perdew
,
V. N.
Staroverov
,
J.
Tao
, and
G. E.
Scuseria
,
Phys. Rev. A
78
,
052513
(
2008
).
29.
C.
Holzer
,
Y. J.
Franzke
, and
M.
Kehry
,
J. Chem. Theory Comput.
17
,
2928
(
2021
).
30.
Y. J.
Franzke
,
F.
Mack
, and
F.
Weigend
,
J. Chem. Theory Comput.
17
,
3974
(
2021
).
31.
A. D.
Becke
,
Phys. Rev. A
38
,
3098
(
1988
).
32.
A. D.
Becke
,
J. Chem. Phys.
88
,
1053
(
1988
).
33.
A. D.
Becke
,
J. Chem. Phys.
104
,
1040
(
1996
).
34.
M. G.
Medvedev
,
I. S.
Bushmarinov
,
J.
Sun
,
J. P.
Perdew
, and
K. A.
Lyssenko
,
Science
355
,
49
(
2017
).
35.
A. D.
Becke
,
J. Chem. Phys.
156
,
214101
(
2022
).
36.
J.
Tao
and
Y.
Mo
,
Phys. Rev. Lett.
117
,
073001
(
2016
).
37.
J.
Sun
,
A.
Ruzsinszky
, and
J. P.
Perdew
,
Phys. Rev. Lett.
115
,
036402
(
2015
).
38.
J. W.
Furness
,
A. D.
Kaplan
,
J.
Ning
,
J. P.
Perdew
, and
J.
Sun
,
J. Phys. Chem. Lett.
11
,
8208
(
2020
).
39.
J. W.
Furness
,
A. D.
Kaplan
,
J.
Ning
,
J. P.
Perdew
, and
J.
Sun
,
J. Phys. Chem. Lett.
11
,
9248
(
2020
).
40.
M.
Dolg
and
X.
Cao
,
Chem. Rev.
112
,
403
(
2012
).
41.
J.
Autschbach
,
J. Chem. Phys.
136
,
150902
(
2012
).
42.
T.
Saue
,
ChemPhysChem
12
,
3077
(
2011
).
43.
P.
Pyykkö
,
Annu. Rev. Phys. Chem.
63
,
45
(
2012
).
44.
K. G.
Dyall
and
K.
Fægri
, Jr.
,
Introduction to Relativistic Quantum Chemistry
(
Oxford University Press
,
New York
,
2007
).
45.
M.
Reiher
and
A.
Wolf
,
Relativistic Quantum Chemistry: The Fundamental Theory of Molecular Science
, 2nd ed. (
Wiley VCH
,
Weinheim, Germany
,
2015
).
46.
Handbook of Relativistic Quantum Chemistry
, edited by
W.
Liu
(
Springer
,
Berlin, Heidelberg
,
2017
).
47.
A.
Wodyński
and
M.
Kaupp
,
J. Chem. Theory Comput.
16
,
314
(
2020
).
48.
A. V.
Arbuznikov
and
M.
Kaupp
,
J. Chem. Phys.
136
,
014111
(
2012
).
49.
P.
Plessow
and
F.
Weigend
,
J. Comput. Chem.
33
,
810
(
2012
).
50.
H.
Bahmann
and
M.
Kaupp
,
J. Chem. Theory Comput.
11
,
1540
(
2015
).
51.
S.
Klawohn
,
H.
Bahmann
, and
M.
Kaupp
,
J. Chem. Theory Comput.
12
,
4254
(
2016
).
52.
T. M.
Maier
,
H.
Bahmann
, and
M.
Kaupp
,
J. Chem. Theory Comput.
11
,
4226
(
2015
).
53.
R.
Grotjahn
,
F.
Furche
, and
M.
Kaupp
,
J. Chem. Theory Comput.
15
,
5508
(
2019
).
54.
M.
Kehry
,
Y. J.
Franzke
,
C.
Holzer
, and
W.
Klopper
,
Mol. Phys.
118
,
e1755064
(
2020
).
55.
C.
Holzer
,
J. Chem. Phys.
153
,
184115
(
2020
).
56.
Y. J.
Franzke
,
C.
Holzer
, and
F.
Mack
,
J. Chem. Theory Comput.
18
,
1030
(
2022
).
57.
C. J.
Schattenberg
,
K.
Reiter
,
F.
Weigend
, and
M.
Kaupp
,
J. Chem. Theory Comput.
16
,
931
(
2020
).
58.
F.
Mack
,
C. J.
Schattenberg
,
M.
Kaupp
, and
F.
Weigend
,
J. Phys. Chem. A
124
,
8529
(
2020
).
59.
C. J.
Schattenberg
and
M.
Kaupp
,
J. Chem. Theory Comput.
17
,
1469
(
2021
).
60.
C. J.
Schattenberg
and
M.
Kaupp
,
J. Phys. Chem. A
125
,
2697
(
2021
).
61.
C. F.
von Weizsäcker
,
Z. Phys.
96
,
431
(
1935
).
62.
C.
Van Wüllen
,
J. Comput. Chem.
20
,
51
(
1999
).
63.
C.
Van Wüllen
,
J. Comput. Chem.
23
,
779
(
2002
).
64.
M. K.
Armbruster
,
F.
Weigend
,
C.
van Wüllen
, and
W.
Klopper
,
Phys. Chem. Chem. Phys.
10
,
1748
(
2008
).
65.
J.
Autschbach
,
Mol. Phys.
111
,
2544
(
2013
).
66.
A.
Baldes
and
F.
Weigend
,
Mol. Phys.
111
,
2617
(
2013
).
67.
H.
Bahmann
,
A.
Rodenberg
,
A. V.
Arbuznikov
, and
M.
Kaupp
,
J. Chem. Phys.
126
,
011103
(
2007
).
68.
A.
Pausch
and
C.
Holzer
,
J. Chem. Phys. Lett.
13
,
4335
(
2022
).
69.
Y.
Mo
,
R.
Car
,
V. N.
Staroverov
,
G. E.
Scuseria
, and
J.
Tao
,
Phys. Rev. B
95
,
035118
(
2017
).
70.
S.
Jana
,
A.
Patra
, and
P.
Samal
,
J. Chem. Phys.
149
,
044120
(
2018
).
71.
A.
Patra
,
S.
Jana
, and
P.
Samal
,
J. Phys. Chem. A
123
,
10582
(
2019
).
72.
R.
Grotjahn
and
M.
Kaupp
,
J. Chem. Phys.
155
,
124108
(
2021
).
73.
C. J.
Umrigar
and
X.
Gonze
,
Phys. Rev. A
50
,
3827
(
1994
).
74.
J. P.
Perdew
,
A.
Ruzsinszky
,
J.
Sun
, and
K.
Burke
,
J. Chem. Phys.
140
,
18A533
(
2014
).
75.
T.
Aschebrock
and
S.
Kümmel
,
Phys. Rev. Res.
1
,
033082
(
2019
).
76.
C.
Ramos
and
B. G.
Janesko
,
J. Chem. Phys.
153
,
164116
(
2020
).
77.
J. W.
Furness
,
N.
Sengupta
,
J.
Ning
,
A.
Ruzsinszky
, and
J.
Sun
,
J. Chem. Phys.
152
,
244112
(
2020
).
78.
T. H.
Dunning
,
J. Chem. Phys.
90
,
1007
(
1989
).
79.
R. A.
Kendall
,
T. H.
Dunning
, and
R. J.
Harrison
,
J. Chem. Phys.
96
,
6796
(
1992
).
80.
D. E.
Woon
and
T. H.
Dunning
,
J. Chem. Phys.
98
,
1358
(
1993
).
81.
Developers’ version of TURBOMOLE V7.6 2021, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989–2007, TURBOMOLE GmbH, since 2007; available from https://www.turbomole.org; retrieved January 12, 2022.
82.
R.
Ahlrichs
,
M.
Bär
,
M.
Häser
,
H.
Horn
, and
C.
Kölmel
,
Chem. Phys. Lett.
162
,
165
(
1989
).
83.
F.
Furche
,
R.
Ahlrichs
,
C.
Hättig
,
W.
Klopper
,
M.
Sierka
, and
F.
Weigend
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
4
,
91
(
2014
).
84.
S. G.
Balasubramani
,
G. P.
Chen
,
S.
Coriani
,
M.
Diedenhofen
,
M. S.
Frank
,
Y. J.
Franzke
,
F.
Furche
,
R.
Grotjahn
,
M. E.
Harding
,
C.
Hättig
,
A.
Hellweg
,
B.
Helmich-Paris
,
C.
Holzer
,
U.
Huniar
,
M.
Kaupp
,
A.
Marefat Khah
,
S.
Karbalaei Khani
,
T.
Müller
,
F.
Mack
,
B. D.
Nguyen
,
S. M.
Parker
,
E.
Perlt
,
D.
Rappoport
,
K.
Reiter
,
S.
Roy
,
M.
Rückert
,
G.
Schmitz
,
M.
Sierka
,
E.
Tapavicza
,
D. P.
Tew
,
C.
van Wüllen
,
V. K.
Voora
,
F.
Weigend
,
A.
Wodyński
, and
J. M.
Yu
,
J. Chem. Phys.
152
,
184107
(
2020
).
85.
Maple
,” version 2020, available from https://www.maplesoft.com/; retrieved December 6, 2020.
86.
Libxc
,” 5.1.7, available from https://www.tddft.org/programs/libxc/; retrieved July 30, 2021.
87.
M. A. L.
Marques
,
M. J. T.
Oliveira
, and
T.
Burnus
,
Comput. Phys. Commun.
183
,
2272
(
2012
).
88.
S.
Lehtola
,
C.
Steigemann
,
M. J. T.
Oliveira
, and
M. A. L.
Marques
,
SoftwareX
7
,
1
(
2018
).
89.
O.
Treutler
and
R.
Ahlrichs
,
J. Chem. Phys.
102
,
346
(
1995
).
90.
O.
Treutler
, “
Entwicklung und Anwendung von Dichtefunktionalmethoden
,” Dr. rer. nat. dissertation [
University of Karlsruhe (TH)
,
Germany
,
1995
].
91.
Y. J.
Franzke
,
N.
Middendorf
, and
F.
Weigend
,
J. Chem. Phys.
148
,
104110
(
2018
).
92.
OpenMP Architecture Review Boards
, “
OpenMP API shared-memory parallel programming
,” https://www.openmp.org; retrieved September 26, 2021.
93.
C.
Holzer
and
Y. J.
Franzke
, OpenMP version of ridft, rdgrad, and egrad with contributions to mpshift, dscf, and grad; improved OpenMP version of aoforce and escf, released with TURBOMOLE V7.4 and further improved in TURBOMOLE V7.5.
94.
O. A.
Vydrov
and
G. E.
Scuseria
,
J. Chem. Phys.
125
,
234109
(
2006
).
95.
J.-D.
Chai
and
M.
Head-Gordon
,
Phys. Chem. Chem. Phys.
10
,
6615
(
2008
).
96.
A.
Karton
,
S.
Daon
, and
J. M. L.
Martin
,
Chem. Phys. Lett.
510
,
165
(
2011
).
97.
Y.
Zhao
,
N.
González-García
, and
D. G.
Truhlar
,
J. Phys. Chem. A
109
,
2012
(
2005
).
98.
Y.
Zhao
,
B. J.
Lynch
, and
D. G.
Truhlar
,
Phys. Chem. Chem. Phys.
7
,
43
(
2005
).
99.
L.
Goerigk
and
S.
Grimme
,
J. Chem. Theory Comput.
6
,
107
(
2010
).
100.
L.
Goerigk
,
A.
Hansen
,
C.
Bauer
,
S.
Ehrlich
,
A.
Najibi
, and
S.
Grimme
,
Phys. Chem. Chem. Phys.
19
,
32184
(
2017
).
101.
F.
Weigend
and
R.
Ahlrichs
,
Phys. Chem. Chem. Phys.
7
,
3297
(
2005
).
102.
C.
Suellen
,
R. G.
Freitas
,
P.-F.
Loos
, and
D.
Jacquemin
,
J. Chem. Theory Comput.
15
,
4581
(
2019
).
103.
F.
Furche
,
B. T.
Krull
,
B. D.
Nguyen
, and
J.
Kwon
,
J. Chem. Phys.
144
,
174105
(
2016
).
104.
P.-F.
Loos
and
D.
Jacquemin
,
J. Chem. Theory Comput.
15
,
2481
(
2019
).
105.
Y. J.
Franzke
and
J. M.
Yu
,
J. Chem. Theory Comput.
18
,
323
(
2022
).
106.
Y. J.
Franzke
and
J. M.
Yu
,
J. Chem. Theory Comput.
18
,
2246
(
2022
).
107.
S.
Gohr
,
P.
Hrobárik
,
M.
Repiský
,
S.
Komorovský
,
K.
Ruud
, and
M.
Kaupp
,
J. Phys. Chem. A
119
,
12892
(
2015
).
108.
Y. J.
Franzke
,
L.
Spiske
,
P.
Pollak
, and
F.
Weigend
,
J. Chem. Theory Comput.
16
,
5658
(
2020
).
109.
Y. J.
Franzke
,
R.
Treß
,
T. M.
Pazdera
, and
F.
Weigend
,
Phys. Chem. Chem. Phys.
21
,
16658
(
2019
).
110.
A.
Klamt
and
G.
Schüürmann
,
J. Chem. Soc., Perkin Trans. 2
1993
,
799
.
111.
A.
Schäfer
,
A.
Klamt
,
D.
Sattel
,
J. C. W.
Lohrenz
, and
F.
Eckert
,
Phys. Chem. Chem. Phys.
2
,
2187
(
2000
).
112.
D.
Peng
and
M.
Reiher
,
J. Chem. Phys.
136
,
244108
(
2012
).
113.
D.
Peng
,
N.
Middendorf
,
F.
Weigend
, and
M.
Reiher
,
J. Chem. Phys.
138
,
184105
(
2013
).
114.
Y. J.
Franzke
and
F.
Weigend
,
J. Chem. Theory Comput.
15
,
1028
(
2019
).
115.
S.
Gillhuber
,
Y. J.
Franzke
, and
F.
Weigend
,
J. Phys. Chem. A
125
,
9707
(
2021
).
116.
R. E.
Stanton
and
S.
Havriliak
,
J. Chem. Phys.
81
,
1910
(
1984
).
117.
S.
Komorovský
,
M.
Repiský
,
O. L.
Malkina
,
V. G.
Malkin
,
I.
Malkin Ondík
, and
M.
Kaupp
,
J. Chem. Phys.
128
,
104101
(
2008
).
118.
P. A. M.
Dirac
,
Proc. R. Soc. London, Ser. A
123
,
714
(
1929
).
119.
J. C.
Slater
,
Phys. Rev.
81
,
385
(
1951
).
120.
S. H.
Vosko
,
L.
Wilk
, and
M.
Nusair
,
Can. J. Phys.
58
,
1200
(
1980
).
121.
T. W.
Keal
and
D. J.
Tozer
,
J. Chem. Phys.
121
,
5654
(
2004
).
122.
J. P.
Perdew
,
Phys. Rev. B
33
,
8822
(
1986
).
123.
J. P.
Perdew
,
A.
Ruzsinszky
,
G. I.
Csonka
,
L. A.
Constantin
, and
J.
Sun
,
Phys. Rev. Lett.
103
,
026403
(
2009
).
124.
J. P.
Perdew
,
A.
Ruzsinszky
,
G. I.
Csonka
,
L. A.
Constantin
, and
J.
Sun
,
Phys. Rev. Lett.
106
,
179902
(
2011
).
125.
A. D.
Becke
,
J. Chem. Phys.
98
,
1372
(
1993
).
126.
A. D.
Becke
,
J. Phys. Chem.
107
,
8554
(
1997
).
127.
P. J.
Wilson
,
T. J.
Bradley
, and
D. J.
Tozer
,
J. Chem. Phys.
115
,
9233
(
2001
).
128.
S.
Grimme
,
J. Phys. Chem. A
109
,
3067
(
2005
).
129.
M.
Bursch
,
H.
Neugebauer
,
S.
Ehlert
, and
S.
Grimme
,
J. Chem. Phys.
156
,
134105
(
2022
).
130.
Y.
Jin
and
R. J.
Bartlett
,
J. Chem. Phys.
145
,
034107
(
2016
).
131.
R. L. A.
Haiduke
and
R. J.
Bartlett
,
J. Chem. Phys.
149
,
131101
(
2018
).
132.
J.
Heyd
,
G. E.
Scuseria
, and
M.
Ernzerhof
,
J. Chem. Phys.
118
,
8207
(
2003
).
133.
J.
Heyd
,
G. E.
Scuseria
, and
M.
Ernzerhof
,
J. Phys. Chem.
124
,
219906
(
2006
).
134.
A. V.
Krukau
,
O. A.
Vydrov
,
A. F.
Izmaylov
, and
G. E.
Scuseria
,
J. Phys. Chem.
125
,
224106
(
2006
).
135.
J. P.
Perdew
and
K.
Schmidt
,
AIP Conf. Proc.
577
,
1
(
2001
).
136.
P.-F.
Loos
,
M.
Comin
,
X.
Blase
, and
D.
Jacquemin
,
J. Chem. Theory Comput.
17
,
3666
(
2021
).

Supplementary Material