We construct an orbital-free non-empirical meta-generalized gradient approximation (GGA) functional, which depends explicitly on density through the density overlap regions indicator [P. de Silva and C. Corminboeuf, J. Chem. Theory Comput. **10**, 3745 (2014)]. The functional does not depend on either the kinetic energy density or the density Laplacian; therefore, it opens a new class of meta-GGA functionals. By construction, our meta-GGA yields exact exchange and correlation energy for the hydrogen atom and recovers the second order gradient expansion for exchange in the slowly varying limit. We show that for molecular systems, overall performance is better than non-empirical GGAs. For atomization energies, performance is on par with revTPSS, without any dependence on Kohn-Sham orbitals.

In Kohn-Sham spin-density functional theory,^{1} in order to evaluate the ground-state energy of a system of electrons moving in an external potential, it is necessary to approximate the exchange-correlation energy functional *E _{xc}*[

*n*

_{↑},

*n*

_{↓}]. The simplest approximation, as well as the starting point for many others, is the Local Density Approximation (LDA), which depends only on the values of spin-densities at each point in space. LDA constitutes the first rung of the “Jacob’s ladder”

^{2}of density functional approximations, in which each higher rung requires additional ingredients beyond spin-densities. The second rung, Generalized Gradient Approximation (GGA), takes into account density gradients, thereby fulfilling more known properties of the exact functional. While LDA and GGAs are designed primarily to respect the paradigm densities of condensed matter physics, i.e., the uniform and slowly varying electron gases, the third rung of functionals (meta-GGAs (MGGAs)) also satisfy some of the exact constraints for the hydrogen atom, the paradigm density of quantum chemistry. Meta-GGAs achieve this by taking advantage of the exact kinetic energy density, in addition to spin-densities and their gradients, giving the following form:

where *n* is the total electron density,

is kinetic energy density, and *ϕ*_{σi}(**r**) are Kohn-Sham orbitals.

Inclusion of the kinetic energy density provides more flexibility in the functional form, which led to the development of accurate semi-empirical^{3–6} and non-empirical^{7–14} meta-GGAs. Due to explicit orbital dependence, most meta-GGAs are only implicit functionals of the density. Nevertheless, functionals that make use of density Laplacian instead of *τ* have also been proposed.^{3,8} The latter are explicit functionals of the density and, in principle, give direct access to exchange correlation potential, without the need of resorting to the optimized effective potential (OEP) method. In this letter, we propose a new class of meta-GGAs that are explicit density functionals. Instead of the kinetic energy density or the electron density Laplacian, we use an inhomogeneity parameter *θ*(**r**), which is the main ingredient of a recently proposed density-dependent bonding descriptor, the density overlap regions indicator (DORI).^{15} The function *θ*(**r**) was designed to reveal regions of space where atomic or molecular densities overlap significantly and its formula reads

where the square of a vector denotes inner product and the right hand side of Eq. (3) is written in terms of the local wave vector **k**(**r**) = ∇*n*(**r**)/*n*(**r**).^{16–18}

For an exponential density, as in a hydrogen-like atom, *θ*(**r**) = 0. For a multi-center and multi-electron finite system, *θ*(**r**) → 0 at the nuclear cusps and density tail, while *θ*(**r**) → ∞ at bond critical points. Therefore, *θ* can be treated as a local measure of deviations from single-exponential behavior and is a modification of previously introduced single-exponential decay detector (SEDD),^{18–20} which uses the uniform electron gas as a reference system and succeeds in resolving semi-quantitatively the atomic shell structure in real space.^{18} Although it may seem that *θ*(**r**) is indeterminate for the uniform electron gas, its value depends on how the uniform limit is approached. In general, this limit will depend on the relation between the first and second derivatives of the density and will vanish under certain assumptions (see Ref. 21 for details). The properties of *θ*(**r**) allow distinguishment between rapidly and slowly varying densities, making this function a suitable ingredient for density functionals. For instance, *θ*(**r**) has been recently used to construct orbital-free mixing functions of local hybrid functionals.^{21} The resulting functionals uniformly perform better than respective global hybrids. Here, we show *θ*(**r**) can be used in semi-local functionals to build in exact constraints relevant to the density of a hydrogen atom. Since calculation of *θ*(**r**) requires second density derivatives, the resulting functional sits on the third rung of the Jacob’s ladder.

One of the commonly used *τ*-dependent variables utilised in the construction of meta-GGA functionals is the iso-orbital detector,

where *τ _{W}*(

**r**) = |∇

*n*|

^{2}/8

*n*is the von Weizsäcker kinetic energy density and

*τ*(

**r**) =

*τ*

_{↑}+

*τ*

_{↓}is the exact kinetic energy density. Wherever electron density is dominated by a single Kohn-Sham orbital, like at nuclei and tails,

*g*(

**r**) tends to 1. On the contrary, in bonding regions, where the density gradient is small, the iso-orbital indicator is close to 0. The same asymptotic behavior can be achieved with a

*θ*-based function by a suitable mapping to the [0–1] range. To this end, we use the following transformation:

Fig. 1 compares $ g \u0303 ( r ) $ and *g*(**r**) for the N_{2} molecule. Both functions go to 1 at nuclei and infinity, but $ g \u0303 ( r ) $ approaches this limit much faster. This steeper behavior of $ g \u0303 ( r ) $ is also reflected in the bonding region where both indicators tend to 0 at the density critical point. Additionally, $ g \u0303 ( r ) $ exhibits a sharp peak in the intermediate regions between cores and tails. This is due to the stationary point of the local wave vector $ \u2207 n ( r ) n ( r ) $, which does not decay strictly monotonically with distance from the nucleus. The ability of $ g \u0303 ( r ) $ to distinguish between regions with different density characters is a promising prerequisite in the development of meta-GGA-level explicit density functionals.

We construct a *θ*-based meta-GGA (*θ*-MGGA) exchange energy functional as an interpolation between different parametrizations of the PBE^{22–24} functional. Consequently, we propose the following form of the exchange energy functional, where, thanks to spin scaling relations,^{25} it is sufficient to consider only the spin-unpolarized version,

where $ \u03f5 x UEG ( n ) =\u2212 3 4 \pi ( 3 \pi 2 n ) 1 / 3 $ is the exchange energy per electron of the uniform electron gas, *p* = |∇*n*|/[4(3*π*^{2})^{2/3}*n*^{8/3}] = *s*^{2} is the squared reduced density gradient, and *θ* is defined in Eq. (3).The enhancement factor *F _{x}* has the same form as in PBE to ensure the Lieb-Oxford bound

where *κ* = 0.804 and *x*(*p*, *θ*) will be defined later. In PBE-type GGAs, *x* = *μp*, where *μ* is a parameter fixed to satisfy a chosen exact constraint. In the original PBE,^{22}^{,}*μ* = *μ _{PBE}* = 0.219 51 was chosen to recover the local spin density approximation (LSDA) linear response by canceling the gradient contribution in the second order expansion of the correlation functional. Such a choice violates the known second order expansion of exchange, which requires

*μ*=

*μ*= 10/81 = 0.123 46. This constraint is recovered instead in the PBEsol

_{GE}^{23}functional, which results in improved accuracy for solids. In the PBEmol

^{24}parametrization,

*μ*=

*μ*= 0.275 83 was chosen to recover the exact exchange energy of the hydrogen atom, i.e.,

_{H}*E*= 0.3125 hartree, the value needed to completely cancel spurious Coulomb self-repulsion.

_{x}The PBEsol and PBEmol parametrizations correspond to two different density regimes, namely, that of slowly varying electron gas and rapidly varying density of the hydrogen atom. Accurate description of the latter case, at the GGA level, requires more than a two times greater gradient contribution. While PBEsol is accurate for solids, PBEmol gives improved results for molecular properties. Nevertheless, even in molecules, the density of bonding regions is closer to the slowly varying limit, as the reduced density gradient is small. Since *θ* can differentiate between the two paradigm densities, it can be used to switch between different parametrizations to correctly account for the local character of the density. To this end, we propose the following formula for *x*(*p*, *θ*) in Eq. (7):

where the mixing function $ g \u0303 ( \theta ) = ( 1 + \theta 2 ) \u2212 1 $ is the same as Eq. (5). Eqs. (6)–(8) define an exchange functional that simultaneously recovers the correct second order gradient expansion and is exact for the hydrogen atom. These two constraints cannot possibly be satisfied at the GGA level and, until now, were not, to the best of our knowledge, exactly satisfied by any explicit density functional.

To visualize how $ g \u0303 ( \theta ) $ enables switching between different parameters, we plot PBEmol, PBEsol, and *θ*-MGGA exchange enhancement factors *F _{x}* (Eq. (7)) for the N

_{2}molecule (Fig. 2). As expected, the

*θ*-MGGA and PBEmol enhancement factors are indistinguishable around the nuclei and far from the molecule. In the bonding region,

*θ*-MGGA closely follows the PBEsol enhancement factor. This is better seen in Fig. 3, which shows the same data focused on the relevant region. The sharp peak in $ g \u0303 ( \theta ) $ (Fig. 1) is reflected by a kink in the enhancement factor. Otherwise, the transition between two density regimes is smooth.

For the correlation functional, we propose a modified version of TPSS correlation,^{7} which is a fairly sophisticated self-interaction correction to the PBE correlation functional. To detect one-electron regions, TPSS uses spin-polarization *ζ*(**r**) = (*n*_{↑} − *n*_{↓})/(*n*_{↑} + *n*_{↓}) and the iso-orbital indicator *g*(**r**). In our *τ*-independent meta-GGA for correlation, *g*(**r**) is replaced by $ g \u0303 \theta ( r ) = ( 1 + \theta ( r ) 2 ) \u2212 1 $ wherever it appears.

PBE-type exchange-correlation functionals recover the LSDA linear response if the following relation between parameters *μ* in the exchange part and *β* (Eq. (4) of Ref. 22) in the correlation part is fulfilled,

As TPSS correlation is built on top of PBE, we retain this property for the limit *θ* → ∞ by replacing the parameter *β* with an interpolation analogous to Eq. (8),

where *β _{H}* = 0.083 84 and

*β*= 0.0375 are fixed to satisfy Eq. (9) for parameters

_{GE}*μ*and

_{H}*μ*, respectively.

_{GE}Contrary to TPSS, the *θ*-based correlation functional *E _{c}*[

*n*] is not exactly self-correlation-free for any one-electron system; however, the form of $ g \u0303 ( r ) $ ensures that the correlation energy is zero for a hydrogen-like atom. For multi-center systems, the self-correlation correction does not switch on in high density overlap regions and

*E*[

_{c}*n*] tends to pure PBE-type correlation with

*β*=

*β*= 0.0375. The corresponding exchange functional tends to PBEsol exchange (

_{GE}*μ*=

*μ*= 10/81); therefore, the exchange and correlation functionals are expected to be well-balanced in terms of error cancellation, as the overall self-interaction correction to PBE becomes effective in the same regions of density.

_{GE}We validate the performance of our *θ*-MGGA functional on a selected set of chemical databases. Since *θ*(**r**) should distinguish exponential regions of the density, we chose to use a slater-type basis set for benchmarking to circumvent possible irregularities due to the more common Gaussian-type basis sets. To this end, we implemented the functional in a locally modified version of the ADF^{26} package and selected the AUG/ATZ2P basis set. All the energy computations are performed post-SCF with PBE orbitals. We verified that differences between self-consistent and non-self-consistent energies are insignificant for other non-empirical meta-GGAs and assume that the same will hold for the *θ*-MGGA. For benchmarking, we chose several subsets from the Minnesota^{27} and GMTKN30^{28} databases. Selections were made to include the main properties of interest for molecular systems, atomization energies, ionization potentials, electron affinities, reaction barrier heights, and covalent and non-covalent bonding energies. With the exception of non-covalent interaction energies data set, we deliberately omitted sets for which dispersion interactions may be substantial for the property of interest. If a subset for a given property was present in both Minnesota and GMTKN30 databases, we chose the one that contained more elements. The data sets included are: 109 main group atomization energies (MGAE109), 9 single-reference metal bond energies (SRMBE9), 10 multi-reference bond energies (MRBE10), 36 adiabatic ionization potentials (G21IP), 25 adiabatic electron affinities (G21EA), 76 barrier heights (BH76), and 31 non-covalent complexation energies (NCCE31). SRMBE9 is a subset of the SRMBE13 set, where compounds containing 5th row elements have been removed due to missing basis functions in the AUG/ATZ2P basis set. Additionally, we included two sets representing cases problematic for density functional approximation: 11 self-interaction error sensitive energies (SIE11) from GMTKN30 and 9 difficult cases (DC9) from the Minnesota database. Table I shows mean absolute errors (MAEs) and mean signed errors (MSEs) for our functional and selected other non-empirical GGA and meta-GGA functionals (PBE, PBEsol, PBmol, TPSS, and revTPSS) for comparison.

. | PBE . | PBEsol . | PBEmol . | TPSS . | revTPSS . | θ-MGGA
. | ||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|

Data set . | MAE . | MSE . | MAE . | MSE . | MAE . | MSE . | MAE . | MSE . | MAE . | MSE . | MAE . | MSE . |

MGAE109 | 13.24 | 11.91 | 36.08 | 36.01 | 7.21 | 0.34 | 4.32 | 1.92 | 3.97 | 0.89 | 3.98 | −0.50 |

SRMBE9 | 7.57 | −6.35 | 10.63 | −10.35 | 6.59 | −4.47 | 6.92 | −5.39 | 7.05 | −5.44 | 6.89 | −3.81 |

MRBE10 | 13.24 | 11.91 | 20.99 | −20.38 | 9.78 | −4.56 | 9.23 | −0.62 | 9.96 | −0.36 | 6.86 | −1.20 |

G21IP | 3.76 | 0.22 | 3.66 | −0.99 | 4.01 | 0.74 | 3.79 | −0.54 | 3.80 | −1.23 | 3.57 | 0.21 |

G21EA | 3.48 | 3.05 | 3.26 | 2.52 | 3.64 | 3.22 | 2.21 | 0.46 | 2.13 | −0.24 | 3.55 | 3.03 |

BH76 | 8.90 | −8.85 | 11.22 | −11.16 | 7.81 | −7.78 | 8.18 | −8.15 | 7.81 | −7.81 | 7.84 | −7.8 |

NCCE31 | 1.17 | 0.32 | 1.69 | 1.04 | 1.12 | 0.06 | 1.17 | −0.26 | 1.14 | −0.17 | 1.31 | −0.60 |

SIE11 | 11.96 | 10.88 | 13.39 | 13.37 | 11.36 | 9.76 | 10.66 | 9.85 | 10.63 | 9.85 | 11.39 | 9.55 |

DC9 | 38.14 | 12.44 | 80.32 | 67.14 | 26.46 | −12.83 | 17.15 | −13.56 | 20.08 | −12.08 | 34.01 | −32.07 |

. | PBE . | PBEsol . | PBEmol . | TPSS . | revTPSS . | θ-MGGA
. | ||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|

Data set . | MAE . | MSE . | MAE . | MSE . | MAE . | MSE . | MAE . | MSE . | MAE . | MSE . | MAE . | MSE . |

MGAE109 | 13.24 | 11.91 | 36.08 | 36.01 | 7.21 | 0.34 | 4.32 | 1.92 | 3.97 | 0.89 | 3.98 | −0.50 |

SRMBE9 | 7.57 | −6.35 | 10.63 | −10.35 | 6.59 | −4.47 | 6.92 | −5.39 | 7.05 | −5.44 | 6.89 | −3.81 |

MRBE10 | 13.24 | 11.91 | 20.99 | −20.38 | 9.78 | −4.56 | 9.23 | −0.62 | 9.96 | −0.36 | 6.86 | −1.20 |

G21IP | 3.76 | 0.22 | 3.66 | −0.99 | 4.01 | 0.74 | 3.79 | −0.54 | 3.80 | −1.23 | 3.57 | 0.21 |

G21EA | 3.48 | 3.05 | 3.26 | 2.52 | 3.64 | 3.22 | 2.21 | 0.46 | 2.13 | −0.24 | 3.55 | 3.03 |

BH76 | 8.90 | −8.85 | 11.22 | −11.16 | 7.81 | −7.78 | 8.18 | −8.15 | 7.81 | −7.81 | 7.84 | −7.8 |

NCCE31 | 1.17 | 0.32 | 1.69 | 1.04 | 1.12 | 0.06 | 1.17 | −0.26 | 1.14 | −0.17 | 1.31 | −0.60 |

SIE11 | 11.96 | 10.88 | 13.39 | 13.37 | 11.36 | 9.76 | 10.66 | 9.85 | 10.63 | 9.85 | 11.39 | 9.55 |

DC9 | 38.14 | 12.44 | 80.32 | 67.14 | 26.46 | −12.83 | 17.15 | −13.56 | 20.08 | −12.08 | 34.01 | −32.07 |

Atomization energies are perhaps the most sensitive probe of a functionals’ accuracy, for which meta-GGAs proved to be significantly better than GGAs. *θ*-MGGA has an absolute error of 3.98 kcal/mol, meaning that it is on par with revTPSS, which has a far more complicated exchange enhancement factor. Our functional clearly outperforms all considered GGAs. For bonding energies, ionization potentials, and reaction barriers, meta-GGAs do not show a systematic or significant improvement. The *θ*-MGGA performs comparably to other non-empirical functionals, with errors for some properties being a bit smaller and a bit larger for others. However, the simple (generalized PBE-like) exchange enhancement factor of *θ*-MGGA appears to be inferior to TPSS/revTPSS exchange for electron affinities and problematic cases. The latter functionals improve over PBE, while *θ*-MGGA gives similar results in terms of absolute errors. On the other hand, *θ*-MGGA has the smallest error for multi-reference bond energies. To complement the performance comparison given in Table I, we list the set of constraints that each of the functionals satisfies (Table II). Constraints satisfied by all of these functionals are not included for brevity.

Constraint . | PBE . | PBEsol . | PBEmol . | TPSS . | revTPSS . | θ-MGGA
. |
---|---|---|---|---|---|---|

LSDA linear response | Yes | No | Yes | No | Yes^{a} | Yes |

2nd order gradient expansion for E _{c} | Yes | No | No | Yes | Yes^{b} | No |

2nd order gradient expansion for E _{x} | No | Yes | No | Yes | Yes | Yes |

4th order gradient expansion for E _{x} | No | No | No | Yes | Yes | No |

E = − 0.3125 hartree for the hydrogen atom _{x} | No | No | Yes | Yes | Yes | Yes |

E = 0 for one-electron systems _{c} | No | No | No | Yes | Yes | Yes^{c} |

Constraint . | PBE . | PBEsol . | PBEmol . | TPSS . | revTPSS . | θ-MGGA
. |
---|---|---|---|---|---|---|

LSDA linear response | Yes | No | Yes | No | Yes^{a} | Yes |

2nd order gradient expansion for E _{c} | Yes | No | No | Yes | Yes^{b} | No |

2nd order gradient expansion for E _{x} | No | Yes | No | Yes | Yes | Yes |

4th order gradient expansion for E _{x} | No | No | No | Yes | Yes | No |

E = − 0.3125 hartree for the hydrogen atom _{x} | No | No | Yes | Yes | Yes | Yes |

E = 0 for one-electron systems _{c} | No | No | No | Yes | Yes | Yes^{c} |

^{a}

For *r _{s}* → ∞.

^{b}

For *r _{s}* → 0.

^{c}

Only for hydrogen-like atoms.

In this letter, we introduce a new class of meta-GGA functionals that depend on the inhomogeneity parameter, *θ*(**r**) (Eq. (3)). This enables construction of a non-empirical explicit density functional, which simultaneously satisfies some exact constraints for the slowly varying electron gas and the hydrogen atom. While meta-GGA functionals satisfying these constraints were previously proposed, they depend on Kohn-Sham orbitals through the kinetic energy density, therefore, are only implicit density functionals. Our new functional takes advantage of *θ*(**r**) to interpolate between different parametrizations of PBE and to replace dependence on the iso-orbital indicator in TPSS correlation. For molecular properties, which are the sole item considered in this work, it generally performs better than PBE-like GGAs and often matches TPSS/revTPSS results.

Funding from EPFL and from Swiss NSF Grant Nos. “200021_137529” and “200021_156001” is gratefully acknowledged.