A local hybrid exchange functional approximation 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.³ 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).¹⁵ 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.¹⁶ Thus, a so-called calibration function^16–18 (CF) is introduced in modern local hybrid functionals, such as LH14t-calPBE¹⁷ and LH20t.¹⁹ 

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).¹⁵ 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),²² 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 limit²⁷ of the correlation (PSTS-LMF).²⁸ We have recently carried out benchmark calculations of electric and magnetic response properties for LHFs based on different LMFs.²⁹ This showed encouraging results for Johnson’s LHF²² 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.³⁰ However, only one functional based on this LMF was optimized,²² which consists of Becke’s 1988 exchange³¹ and correlation terms.³² 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 term³³ 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 functionals^36–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-LMF²² and PSTS-LMF²⁸ 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.⁴⁷ 

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.²² 

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

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.⁴⁸ To carry out the integration in Eq. (1), Plessow and Weigend suggested to use a seminumerical scheme⁴⁹ 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:

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

with g denoting a grid point. This LHF scheme was later applied self-consistently⁵⁰ 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,²⁹ 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.²⁹ 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.¹⁵ The respective t-LMF reads

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

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

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

where c_σσ′ are given by c_αβ = 0.63 and c_σσ = 0.88.³² 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

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 as^62,63

In a basis set representation, the non-collinear exchange–correlation energy depends on the total density matrix M₀ and the three spin vector density matrices M_i with i = x, y, z. These are defined as^62,63

The XC potential follows as

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

Only considering M₀ and M_z results in the UKS limit. Therefore, an unrestricted Kohn–Sham implementation of the XC potential and the XC energy can be extended straightforwardly.⁶⁴ 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-SVWN⁶⁷ or LH14t-calPBE¹⁷ is not yet available. Moreover, current-dependent terms arise as spin–orbit coupling is closely related to magnetic induction, and these are neglected herein;⁶⁸ 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 e_X as usual using an enhancement factor F_X,

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

As in many recent density functional approximations (DFAs), we chose a semilocal enhancement factor, where F_X 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,³⁶ 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

The auxiliary quantities R_X and f_X are defined as

and

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.³⁶ Furthermore, y = (2λ − 1)²p is a scaled version of the reduced density gradient s depending on the coordinate transformation with p = s². 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.³⁶ We use the same expressions as Tao and Mo for the correction in the slowly varying limit,

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

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.²² Using an approximated correlation length z^DME,

we define U as

with c_F = 3/8 · 4^2/3(3/π)^1/3 and ζ = (ρ_σ − ρ_σ′)/(ρ_σ + ρ_σ′). U_σ′ is obtained by reversing the spin indices. The expressions of R_L and f_L 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 a^DME as

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.⁷² 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 c_L, λ_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.³⁶ 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,³⁶ 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 c_L, 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$⁠.⁷⁴ 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$⁠.⁷⁵ 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.⁷³ For the low-density limit, He is used, while we chose Hg⁷⁸⁺ 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 Hg⁷⁸⁺, 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.⁷³ Multiplying with the analytic prefactors yields −1.039 and −49.781 hartree, respectively. Finally, we optimize β_X for a given value of c_L, yielding possible pairs of solutions. We chose the pair with the highest value for c_L, yielding c_L = 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 c_L = 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. c_L = 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 functional³³ as the latter has shown to work well in conjunction with local hybrid exchange functionals.¹⁹ We note that the B95 correlation functional uses two parameters optimized for the He and Ne atoms.³³ Attempts to employ the (modified) TPSS correlation^9,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.¹⁰ Viable alternatives to B95 correlation may be obtained from the rung-3.5 Ansatz of Ramos and Janesko⁷⁶ or from the recent non-dynamical correlation Ansatz of Becke.³⁵ 

Figure 1 shows a comparison of the LMF defined by Johnson²² 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).

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 Li₂. 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 Li₂ 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.⁷⁷ TMHF, however, converges smoothly for the stretched Li₂. 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.

The outlined functionals were implemented into the TURBOMOLE program suite.^81–84 The mathematical functions for the LMFs were generated using Maple scripts,⁸⁵ 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 implementation⁴⁷ 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 function¹⁸ and general meta-GGAs for LHFs. Therefore, the LH20t functional¹⁹ is also now available in two-component open-shell calculations. We further added interfaces to Libxc^86–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.⁹¹ 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,²² 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);⁷ the meta-GGA TPSS;⁹ Tao–Mo;³⁶ SCAN;³⁷ the hybrid functionals PBE0,^7,8 B3LYP,^4–6 and TPSSh;^9,10 the range-separated functionals CAM-B3LYP,¹⁴ LC-ωPBE,⁹⁴ and ωB97X-D;⁹⁵ and the local hybrid functionals LH12ct-SsirPW92,⁴⁸ LH14t-calPBE,¹⁷ LH20t,¹⁹ LHJ14,²² and mPSTS-a1.²⁸ 

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,⁹⁶ 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.¹⁰⁰ To yield values that can be directly compared to the GMTKN values, the def2-QZVP basis set¹⁰¹ 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.¹⁰² Structures were taken from Ref. 102. In line with this reference, we use the aug-cc-pVTZ basis set^78–80 and large grids (grid size 4) for the numerical integration of the DFT parts.^89,90 Tight SCF thresholds of 10⁻⁹ hartree for the energy and 10⁻⁷ a.u. for the root mean square of the density matrix are applied. Response equations are converged with a threshold of 10⁻⁵ a.u. for the norm of the residuum.¹⁰³ The excitation energies are corrected by the B3LYP zero-point vibrational energies.¹⁰⁴ 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.⁴⁷ 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, [MoNCl₄]²⁻, [MoOF₄]⁻, [MoOCl₄]⁻, [MoOF₅]²⁻, [MoOBr₅]²⁻, [WOCl₄]⁻, [WOF₅]²⁻, [WOBr₅]²⁻, [TcNF₄]⁻, [TcNCl₄]⁻, [TcNBr₄]⁻, [ReNF₄]⁻, [ReNCl₄]⁻, [ReNBr₄]⁻, [ReOBr₄], [ReOF₅]⁻, and [OsOF₅]. We use the same computational parameters as in Refs. 105 and 106. In detail, the x2c-QZVPall-2c basis set¹⁰⁸ 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) condition¹¹⁶ is employed for the HFC, whereas both RKB and the restricted magnetic balance (RMB) conditions¹¹⁷ are employed for the g-tensor. An SCF threshold of 10⁻⁹ 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,¹²¹ BP86,^31,122 revTPSS,^123,124 r²SCAN,^38,39 BH&HLYP,^5,31,125 PBE0, including 40% of HF exchange (PBE0-40HF),^7,8,107 B97,¹²⁶ B97-2,¹²⁷ revTPSSh,^123,124 TPSS0,^10,128 r²SCANh,^38,39,129 r²SCAN0,^38,39,129 r²SCAN50,^38,39,129 CAM-QPT-00,¹³⁰ CAM-QTP-02,¹³¹ HSE06,^132–134 LH12ct-SsifPW92,⁴⁸ LH20t* (LH20t without calibration function),¹⁹ and mPSTS-noa2.^28,29 Note that SCAN has been replaced by r²SCAN 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-11⁹⁶ and BH76 subsets^97–99 from the GMTKN55 test set¹⁰⁰ 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.

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.

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.¹⁰² For these molecules, TMHF and TMHF-3P perform very well too, as shown in Fig. 5, with TMHF having the edge.

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.¹⁰² 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%).

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 [TbPc₂]⁻ 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 [TbPc₂]⁻ 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,¹³⁵ 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 [TbPc₂]⁻, 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.