Simple exchange hole models for long-range-corrected density functionals

TL;DR

This paper introduces simple, closed-form exchange hole models for long-range-corrected density functionals, improving the conversion of semilocal functionals into their short-range analogs with explicit formulas.

Contribution

It develops new exchange hole models using Hermite functions for better short-range functional approximations in density functional theory.

Findings

Models match the uniform electron gas limit

Energy densities are within 5% of each other

New models are non-oscillatory and simple to implement

Abstract

Density functionals with a range-separated treatment of the exchange energy are known to improve upon their semilocal forerunners and fixed-fraction hybrids. The conversion of a given semilocal functional into its short-range analog is not straightforward, however, and not even unique, because the latter has a higher information content that has to be recovered in some way. Simple models of the spherically-averaged exchange hole as an interpolation between the uniform electron gas limit and a few-term Hermite function are developed here for use with generalized-gradient approximations, so that the energy density of the error-function-weighted Coulomb interaction is given by explicit closed-form expressions in terms of elementary and error functions. For comparison, some new non-oscillatory models in the spirit of earlier works are also built and studied, their energy densities match rather closely (within less than 5%) but do lack the exact uniform electron gas limit.

Equations

$$E_{\mathrm{x}}=c_{\mathrm{x}}\int\rho^{4/3}({\rm\bf r})f\bigl{(}s({\rm\bf r})\bigr{)}\,\mathrm{d}^{3}{\rm\bf r},$$

$$c_{\mathrm{x}}=-\frac{3}{4}\sqrt[3]{\frac{3}{\pi }},$$

$$s({\rm\bf r})=\frac{\bigl{|}\nabla\rho({\rm\bf r})\bigr{|}}{2p({\rm\bf r})\rho({\rm\bf r})},$$

$$p(\mathbf{r})=\sqrt[3]{3{\pi }^2\rho (\mathbf{r})}.$$

$$E_{\mathrm{x}}=\frac{1}{2}\int \rho (\mathbf{r})\int \frac{{\rho }_{\mathrm{x}}\left(\mathbf{r},\mathbf{r}+\mathbf{R}\right)}{|\mathbf{R}|}{\mathrm{d}}^3\mathbf{R}{\mathrm{d}}^3\mathbf{r}$$

$$\rho_{\mathrm{x}}({\rm\bf r},{\rm\bf r}+{\rm\bf R})=\rho({\rm\bf r})\,q\bigl{(}p({\rm\bf r})\,|{\rm\bf R}|,s({\rm\bf r})\bigr{)}$$

$$E_{\mathrm{x}}^{\omega}={\textstyle\frac{1}{2}}\int\rho({\rm\bf r})\int\frac{\rho_{\mathrm{x}}\left({\rm\bf r},{\rm\bf r}+{\rm\bf R}\right)}{|{\rm\bf R}|}\operatorname{erfc}\bigl{(}\omega|{\rm\bf R}|\bigr{)}\,\mathrm{d}^{3}{\rm\bf R}\,\mathrm{d}^{3}{\rm\bf r}$$

$$E_{\mathrm{x}}^{\omega }=c_{\mathrm{x}}\int {\rho }^{4/3}(\mathbf{r})f\left(\frac{p(\mathbf{r})}{\omega },s(\mathbf{r})\right){\mathrm{d}}^3\mathbf{r}$$

$$f(\lambda ,s)=-\frac{8}{9}\underset{0}{\overset{\infty }{\int }}r\mathrm{erfc}\left(\frac{r}{\lambda }\right)q(r,s)\mathrm{d}r.$$

$$\frac{4}{3\pi }\underset{0}{\overset{\infty }{\int }}{r}^{2}q(r,s)\mathrm{d}r=-1,$$

$$\frac{8}{9}\underset{0}{\overset{\infty }{\int }}rq(r,s)\mathrm{d}r=-f(s),$$

$$\underset{r,s\to 0}{\lim }q(r,s)=-\frac{1}{2}+\left(\frac{1}{10}-\frac{1}{27}{s}^2\right)r^2+\dots$$

$$\tilde{q}(r)=-\frac{9}{2}{\left(\frac{\sin r-r\cos r}{r^3}\right)}^2$$

$$q(r,0)=\tilde{q}(r),$$

$$q(r,s)=\bigl{(}1-w(s)\bigr{)}\tilde{q}(r)+w(s)\,q_{2}\bigl{(}r,a(s),c(s)\bigr{)}$$

$$q_2(r,a,c)$$

$$c(s)=4\sqrt{\pi }{a}^5(s)-\frac{9}{2}{a}^4(s)\left(1+\frac{f(s)-1}{w(s)}\right)+\frac{1}{3}{a}^2(s),$$

$$w(s)=1-\exp \left(-\gamma s^2\right)$$

$$\gamma =\frac{6075{\alpha }^4\mu -20}{4320\sqrt{\pi }{\alpha }^5-6075{\alpha }^4+900{\alpha }^2-54}$$

$$\underset{s\to 0}{\lim }f(s)=1+\mu s^2+\dots$$

$$f(\lambda,s)=\bigl{(}1-w(s)\bigr{)}\tilde{u}(\lambda)+w(s)u_{2}\bigl{(}\lambda,a(s),c(s)\bigr{)}$$

$$\tilde{u}(\lambda )=1+\frac{2}{3{\lambda }^2}+\left(4-\frac{2}{{\lambda }^2}\right)\frac{1-\exp \left(-{\lambda }^2\right)}{3{\lambda }^2}-\frac{4\sqrt{\pi }\mathrm{erf}(\lambda )}{3\lambda }$$

$$\tilde{u}(\lambda )=\sum _{n=1}^{\infty }\frac{(-1{)}^{n-1}}{(n+2)!(n+\frac{1}{2})}{\lambda }^{2n},$$

$$u_2(\lambda ,a,c)$$

$$\begin{array}[]{lcl}v_{n}(t)&=&\upsilon_{n}\Bigl{(}\sqrt{1+t^{2}}\Bigr{)},\\ \upsilon_{0}(x)&=&1/\bigl{(}2x(1+x)\bigr{)},\\ \upsilon_{1}(x)&=&\bigl{(}1+2x\bigr{)}/\bigl{(}4x^{3}(1+x)^{2}\bigr{)},\\ \upsilon_{2}(x)&=&\bigl{(}3+9x+8x^{2}\bigr{)}/\bigl{(}8x^{5}(1+x)^{3}\bigr{)}.\end{array}$$

$$f(s)=1+\kappa -\kappa /\left(1+\mu s^2/\kappa \right),$$

$$f(s)=1+\kappa -\kappa \exp \left(-\mu s^2/\kappa \right),$$

$$\mu \approx 0.219514972757702295113021595185752214898986$$

$$\underset{s\to \infty }{\lim }c(s)=0,$$

$$24\sqrt{\pi }{A}^3-27(1+\kappa )A^2+2=0,$$

Full text

It is the generalized-gradient approximation Langreth and Mehl (1983); Becke (1986); Perdew and Wang (1986); Perdew et al. (1996) that paved the way for the density functional theory Hohenberg and Kohn (1964); Kohn and Sham (1965) into the mysterious kingdom of theoretical chemistry. Even more fruitful may seem to be the hybrids Becke (1993, 1993); Perdew et al. (1996) with a fixed fraction of exact exchange, they are widely used, but their “strange” asymptotic behavior of the effective potential is more than an æsthetic problem. Luckily, a wonderful solution Iikura et al. (2001) was found by splitting the two-electron interaction within the exchange energy into the short- and long-range parts and using a density-functional approximation for the former and the full exact exchange for the latter (the general idea has a longer history Gill et al. (1996); Leininger et al. (1997)). This was soon shown Tawada et al. (2004) to be even more helpful to the time-dependent Runge and Gross (1984); Petersilka et al. (1996) density functional theory where it greatly improves the calculated excited state properties and overcomes the failure for charge-transfer excitations Dreuw et al. (2003).

Given a well-tested semilocal density functional for exchange, it is not straightforward to get its short-range analog because the latter has a higher information content that cannot be recovered uniquely. The earliest studies Iikura et al. (2001); Tawada et al. (2004) took a somewhat simplistic shortcut that breaks the underlying sum rules, while a consistent construction should be based on an explicit model of the exchange hole — an entity deeply rooted in the adiabatic-connection approach Harris (1984). An elegant analytic model Ernzerhof and Perdew (1998) designed around a non-oscillatory (nodeless) approximation Perdew and Wang (1992) for the uniform electron gas was the first to be used Heyd et al. (2003) in this role, but the lack of closed-form expressions for the needed integrals led to a further work Henderson et al. (2008) where a similar but more computationally tractable nodeless function has been built and proved successful both in applications Rohrdanz et al. (2009); Weintraub et al. (2009) and as a starting point for more sophisticated developments Tao et al. (2017). Other models are known, those based on an atomic-like exchange hole Becke and Roussel (1989) are reported Modrzejewski et al. (2016) that satisfy fewer exact constraints, as well as oscillatory models Tao and Mo (2016); Patra et al. (2018) based on a density matrix expansion Negele and Vautherin (1972); Koehl et al. (1996); Tsuneda and Hirao (2000).

As we wanted to use the long-range corrected functionals for the good of chemistry, we could not blindly adopt any such model, we did not like the need for fitting a function to the numerical solution of a parametrized nonlinear equation Henderson et al. (2008), we were also slightly worried about the lack of the exact uniform electron gas limit by any nodeless exchange hole model. We have found new and simpler explicit solutions in closed form that should work no less well and are easy to deal with.

A generalized-gradient approximation for the exchange energy has a simple functional form

[TABLE]

[TABLE]

with all its wisdom condensed in the enhancement factor $f(s)$, a function of only one dimensionless variable

[TABLE]

[TABLE]

On the other hand, the exact exchange energy

[TABLE]

can be given in terms of the exchange hole $\rho_{\mathrm{x}}({\rm\bf r}_{1},{\rm\bf r}_{2})$ whose spherically-averaged part is only needed and is then approximated

[TABLE]

using the shape function $q(r,s)$ that holds more information than is otherwise hidden, by the integration, behind $f(s)$. If the shape function is known, the error-function-weighted short-range part of the exchange energy

[TABLE]

can be cast in the form

[TABLE]

with the new enhancement factor $f(\lambda,s)$ now being a function of two dimensionless variables ($\lambda$ is length-like),

[TABLE]

Finding a good shape function $q(r,s)$ given an enhancement factor $f(s)$ is the problem we want to solve here. In doing so, we should respect the sign of $q(r,s)\leq 0$, the normalization

[TABLE]

and the energy connection

[TABLE]

while the known on-top value and curvature Becke (1983)

[TABLE]

are very helpful to build a good overall shape. The uniform electron gas has an oscillatory function

[TABLE]

with a rather long tail of $-{\textstyle\frac{9}{4}}r^{-4}$, whereas finite band gap systems have it more localized and mostly smooth. What we have written up to here is the common knowledge Perdew and Wang (1992); Perdew et al. (1996); Ernzerhof and Perdew (1998) in the field, with all this in mind, we will now build and compare the new models of our own.

We want Eq. (9) to have the exact uniform electron gas limit at $s=0$, and the only way to meet this is when

[TABLE]

so our first model will be an interpolation

[TABLE]

between $\tilde{q}(r)$ and a three-term Hermite function

[TABLE]

It already follows Eq. (10), while from Eq. (11) we get

[TABLE]

and nothing seems to be more natural than

[TABLE]

with $\gamma$ set as

[TABLE]

to fulfill Eq. (12), here $\alpha\equiv a(0)$ and $\mu$ is from

[TABLE]

(understanding that $f(s)$ should always bee Becke (1988) an even function of $s$). After all this, we are left with the freedom to choose a good function $a(s)$ limited mainly by the sign of $q(r,s)\leq 0$.

The integral of Eq. (9) over the function of Eq. (15) has a simple closed-form expression

[TABLE]

with the known Gill et al. (1996) uniform electron gas function

[TABLE]

which for small $\lambda$ should be evaluated using (a few terms of) the series

[TABLE]

and the well-behaved functions

[TABLE]

[TABLE]

There are two kinds of enhancement factors: either bounded by a constant, $1\leq f(s)\leq 1+\kappa$, or unbounded $f(s)\to\infty$ as $s\to\infty$. We will deal first with those of the former kind, the simplest Becke (1986) and widely used Perdew et al. (1996)

[TABLE]

and another useful Hammer et al. (1999)

[TABLE]

both having only two parameters derivable from first principles, the gradient coefficient Ma and Brueckner (1968)

[TABLE]

(we have carefully computed the integrals Ma and Brueckner (1968) numerically to all digits given), and an estimated Lieb and Oxford (1981) $\kappa=0.804$ from the global lower bound Lieb (1979) on the exchange energy.

Our simplest $a(s)$ in Eqs. (15) and (17) is then a constant $a(s)=A$ whose value can be nailed down by setting

[TABLE]

so that $A$ is a root of the cubic equation

[TABLE]

for $\kappa=0.804$ we get

[TABLE]

and for $\mu$ of Eq. (28), from Eq. (19) with $\alpha=A$,

[TABLE]

This is our simplest model that can also work with other more flexible Adamo and Barone (2002) forms of $f(s)$ as long as they are bounded by a constant, it is straightforward to implement and it has, through Eq. (17), the input $f(s)$ as a multiplicative factor in the expression for $f(\lambda,s)$. Plots show that $c(s)>0$, monotonic for Eq. (27) but with a slight wave up and down for Eq. (26).

As a prototype of an unbounded $f(s)$, we take the most well-known and widely used Becke (1988)

[TABLE]

[TABLE]

where $\mu$ can be either adjusted Becke (1988) to fit some data or Campo (2016) the theoretical constant of Eq. (28), $\nu$ is fixed by the asymptotic behavior of the energy density, (and we must note that $\eta$ could as well have been an adjustable parameter — its value of Eq. (35) is nothing but arbitrary). Here, we should have an $a(s)$ that always grows with $s$, otherwise there would have been $c(s)<0$ and $q(r,s)>0$. To meet Eq. (29), it can be shown that the first two terms of

[TABLE]

would have been needed, and the third and higher terms would help $c(s)$ reach zero faster. We cannot take these first two terms exactly as written for $a(s)$, however, because there would be $c(s)<0$ for some small $s$, but the simplest

[TABLE]

already yields a working overall solution.

By the way, putting the bounded $f(s)$ of Eqs. (26) or (27) into Eq. (37) would also work and give us another $f(\lambda,s)$ of Eq. (21) that is clearly not the same as our first model with the constant $a(s)=A$, and when we plot the ratio of these $f(\lambda,s)$, we see that the one based on Eq. (37) is down to $15\%$ smaller for some $\lambda$ and $s$. This gives us a hint at their diversity and makes us think of how to narrow down the choice of $a(s)$. Besides Eq. (29), we can nail it down at the other end, $s=0$, by

[TABLE]

which makes $\alpha\equiv a(0)$ the root of the seventh-degree polynomial parametrized by $\mu$, that can be written as

[TABLE]

for $\alpha$ has to be solved for together with $\gamma$ of Eq. (19). For $\mu$ of Eq. (28) we get

[TABLE]

and this $\gamma$ is roughly $5$ times greater than that of Eq. (32), we think the greater $\gamma$ to be better because then the oscillatory $\tilde{q}(r)$ fades away more quickly in Eq. (15). It is easy to build monotonic interpolations $\alpha\leq a(s)\leq A$ for a bounded $f(s)$ to get a small $c(s)>0$: for Eq. (26)

[TABLE]

yields $c(s)<0.023$ for all $s$; whereas for Eq. (27)

[TABLE]

yields $c(s)<0.016$ for all $s$. Likewise, for Eq. (33)

[TABLE]

makes $c(s)<0.036$ for all $s$. This experience helps us get rid of $c(s)$ altogether, ending up with an even simpler model

[TABLE]

[TABLE]

that can be used in two ways: either by redefining

[TABLE]

for use with some $a(s)$ like in Eqs. (41), (43), and (Simple exchange hole models for long-range-corrected density functionals); or by holding true to Eq. (18) while fearlessly solving the cubic equation for $a(s)$ to get

[TABLE]

[TABLE]

This last idea is so strikingly simple that nothing is left to be shaved away with Ockham’s razor, and we like it the most.

Thus, given a $f(s)$ of any meaningful kind, we find its $\mu$ of Eq. (20), get $\alpha$ from Eq. (39), and $\gamma$ from Eq. (19), so we have $w(s)$ of Eq. (18), hence $a(s)$ of Eqs. (49) and (50), that yields us $f(\lambda,s)$ of Eq. (47) with Eqs. (22) and (25).

Here our tale would have had a happy end, but we feel that someone may call it a heresy to work with an oscillatory shape function having a thin but too long tail. In the spirit of the early works Perdew and Wang (1992); Ernzerhof and Perdew (1998); Henderson et al. (2008), we will now build some new and simple non-oscillatory models (of interest on their own) and compare the outcomes $f(\lambda,s)$ one-to-one to see only a small difference.

We begin with our two new amazingly beautiful non-oscillatory exchange hole models for the uniform electron gas that follow Eqs. (10), (11), (12) and have the $-{\textstyle\frac{9}{4}}r^{-4}$ tail from Eq. (13): the split-exponent version

[TABLE]

$\alpha$ and $\beta$ being roots of the polynomial system

[TABLE]

[TABLE]

and the shared-exponent version

[TABLE]

[TABLE]

[TABLE]

In both cases, all three functions $q(r)$, $rq(r)$, and $r^{2}q(r)$ have slim shapes without any shoulders for $r\geq 0$, whereas their forerunners Ernzerhof and Perdew (1998); Henderson et al. (2008), to become shoulderless, needed one more degree of freedom to be fixed by a sophisticated information-entropy-maximization Jaynes (1957) principle. We hope that our finding may help others in their future work.

From Eq. (Simple exchange hole models for long-range-corrected density functionals), we build

[TABLE]

with $a(0)=\alpha$ of Eq. (55), while $b(s)$ and $c(s)$ have to fulfill Eqs. (10) and (11),

[TABLE]

using the integrals over the first $r^{-4}$ term of Eq (57),

[TABLE]

[TABLE]

for $x\gg 1$, these should be computed as

[TABLE]

To follow the curvature of Eq. (12) at $s=0$, we need

[TABLE]

The first term in Eq. (57) should smoothly switch from having the $-{\textstyle\frac{9}{4}}r^{-4}$ tail to a short-range exponential behavior, to overcome the logarithmic singularity as $s\to 0$ in its integrals over $r$, such as in Eqs. (61) and (63), we multiply $a(s)$ by the healing function $h(\chi s)$,

[TABLE]

that has its value and all derivatives zero at $s=0$. For $\chi$, we can take either the greatest $\chi_{0}$ that still yields $b(s)\geq 0$, or the greatest $\chi_{1}$ that still yields a monotonic $b(s)$, by solving

[TABLE]

for $s_{n}$ and $\chi_{n}$, where $b^{(n)}(s;\chi)\equiv\partial^{n}b(s;\chi)/\partial s^{n}$; we can also set $\chi=0$ to see what happens. Given some $a(s)$, there is an explicit closed-form expression for Eq. (9),

[TABLE]

[TABLE]

[TABLE]

For $f(s)$ of Eq. (26), we take $a(s)$ from Eq. (41) with $\alpha$ from Eq. (56), $a^{\prime\prime}(0)=2a_{2}$ for $\gamma$ of Eq. (66), solve Eq. (68) to get

[TABLE]

and now we can compare $\bar{f}(\lambda,s)$ of Eq. (Simple exchange hole models for long-range-corrected density functionals) to our best $f(\lambda,s)$ of Eqs. (47) and (49). In Fig. 1 the two functions are plotted and we see how regular they are and how little they differ, even more impressive is the colorful family of curves in Fig. 2 for their ratio $\bar{f}(\lambda,s)/f(\lambda,s)$. This way, we get a measure of their similarity,

[TABLE]

In the uniform electron gas limit, the nodeless shape function of Eq. (Simple exchange hole models for long-range-corrected density functionals) yields the integral $\bar{f}(\lambda,0)$ that matches the exact $\tilde{u}(\lambda)$ of Eq. (22) to within $1.5\%$, while for $s>0$ the models are only a few times farther away from each other, being closest for $\chi=0$. Thus, $\chi$ plays no dramatic role, and it seems better to cut the tail in depth by $\exp\left(-\gamma s^{2}\right)$ than at length by $\exp\left(-h^{2}(\chi s)\,a^{2}(s)\,r^{2}\right)$, to enjoy a rewarding simplification of the equations. In this way, the function of Eq. (51) can also be used, and a cubic equation for $a(s)$ can then be set up, but we leave it out here to save space.

It is now clear that both kinds of shape functions — both the oscillatory of Eq. (46) and the non-oscillatory of Eq. (57) — would yield nearly the same integral output of Eq. (9) under the same constraints of Eqs. (10), (11), and (12). To our mind, the oscillatory function gives the best solution: we get an explicit closed-form expression for $f(\lambda,s)$ in terms of the given $f(s)$ using Eqs. (18), (49), (50), and (47); furthermore, it has the exact uniform electron gas limit. Nevertheless, our experience with the non-oscillatory functions was not in vain and these can be used in the further work on new functionals.

It might be time for a thorough benchmark of the new model on a wide set of molecules, but we put it off for now until we learn how to combine it with a dispersion-correction functional Dion et al. (2004); Vydrov and van Voorhis (2010).

We find our long-range corrected version of the PBE Perdew et al. (1996) functional with $\frac{1}{\omega}=3\;\mathrm{au}$ (an easy-to-remember whole number) to be already a good next step after its ${\textstyle\frac{1}{4}}$-fixed-fraction hybrid Perdew et al. (1996); Adamo and Barone (1999), and it can be used routinely in mechanistic studies of molecular structure and reactivity toward a full understanding of chemical kinetics.