Local Density Approximation for the Short-Range Exchange Free Energy Functional.

Analytical expressions for the exchange free energy per particle of the uniform electron gas (UEG) associated with the short-range (SR) interelectronic interaction at the low- and high-temperature limits are examined, yielding an accurate analytical parametrization for the SR exchange free energy per particle of the UEG as a function of the uniform electron density, temperature, and range-separation parameter. This parametrization constitutes the local density approximation for the SR exchange free energy functional, which can be the first step toward finding generally accurate range-separated hybrid functionals in both finite-temperature density functional theory and thermally assisted-occupation density functional theory.

Introduction

Over the past decades, Kohn–Sham density functional theory (KS-DFT)[1,2] has been one of the most powerful quantum-mechanical methods for studying the ground-state properties of atoms, molecules, and bulk materials.[3,4] However, since the exact exchange–correlation (XC) energy functional Exc[ρ] in KS-DFT remains unknown, numerous efforts have been devoted to improving the accuracy of approximate XC energy functionals.

The uniform electron gas (UEG) is an important system from which the local density approximation (LDA) for the XC energy functional, ExcLDA, can be developed. However, a major shortcoming associated with the LDA is the self-interaction error (SIE),[5] which can be efficiently reduced by incorporating the Hartree–Fock (HF) exchange energy ExHF into the LDA XC energy functional ExcLDA. A global hybrid (GH) XC energy functional[6,7] may generally achieve this. However, to reduce the SIE problem more effectively, the range-separated hybrid (RSH) scheme can be adopted.[8−19] Particularly, in the long-range corrected hybrid (LCH) scheme, the Coulomb operator is partitioned into the long-range (LR) and short-range (SR) operators, commonly achieved with the use of the standard error function (erf) and the complementary error function (erfc)On the right-hand side of eq , the first term is the LR operator (i.e., the erf operator), the second term is the SR operator (i.e., the erfc operator), and ω is the range-separation parameter controlling the range of the SR operator (atomic units are used throughout the paper). The LCH XC energy functional can be commonly expressed as ExcLCH = ExLR-HF + ExSR-LDA + EcLDA, where ExLR-HF is the HF exchange energy associated with the LR operator, ExSR-LDA is the LDA exchange energy functional associated with the SR operator, and EcLDA is the LDA correlation energy functional (i.e., associated with the Coulomb operator). Accordingly, the SIE associated with an LCH XC energy functional can be greatly reduced, and the corresponding XC potential possesses the correct −1/r asymptote in the asymptotic regions of atomic and molecular systems. At ω = 0 (i.e., the SR interelectronic interaction reduces to the Coulomb interelectronic interaction), ExcLCH reduces to ExcLDA. Note also that in the LCH scheme, ExSR-LDA and EcLDA, which are the functionals based on the LDA (i.e., the simplest density functional approximation), may be replaced with those based on more sophisticated density functional approximations [e.g., the generalized-gradient approximations (GGAs)] for improved accuracy.

Finite-temperature density functional theory (FT-DFT) was first proposed by Mermin,[20] wherein the grand canonical ensemble is adopted to study the thermodynamic properties of a physical system at temperature T. To obtain the thermal equilibrium density in FT-DFT, one needs to solve the following self-consistent equationswhereis the effective potential, vext(r) is the external potential, and μxc(r) = δFxc[ρ]/δρ(r) is the XC potential (i.e., the functional derivative of the XC free energy functional Fxc[ρ]). The thermal equilibrium density is given bywhere μ is the chemical potential chosen to conserve the number of electrons and kB is the Boltzmann constant. At T = 0, FT-DFT reduces to KS-DFT and the XC free energy functional Fxc[ρ] reduces to the XC energy functional Exc[ρ].

Although FT-DFT has been proposed for several decades, only KS-DFT has received massive applications in solid-state physics and quantum chemistry. The main reason may be that the “thermal effect” is commonly regarded as an effect only on nuclei. According to the Born–Oppenheimer approximation, the Hamiltonian of electrons can be separated from that of nuclei. The coupling between the vibrational modes of nuclei and the motion of electrons can be properly described by some model Hamiltonian and treated as perturbation.[21] In many other cases, electrons are not in the extreme conditions and hence the thermal effect can simply be neglected. Therefore, the properties of atoms, molecules, and bulk materials are usually studied using KS-DFT (i.e., FT-DFT with T = 0). Although several efforts have been devoted to understanding the properties of thermal density functionals,[22−26] the applications of FT-DFT are only limited to hot plasmas and warm dense matter.[27−29] Besides, in conventional FT-DFT calculations, for simplicity, the XC free energy functional Fxc[ρ] has been commonly approximated by the “zero-temperature” XC energy functional Exc[ρ] evaluated with the thermal equilibrium density ρ(r) at temperature T.

In 1984, an analytical parametrization for the LDA XC free energy functional FxcLDA[ρ] was proposed by Perrot and Dharma-wardana,[30] based on the low-temperature limit for the exchange free energy derived by Horovitz and Thieberger.[31] In 2000, a classical mapping method[32] was adopted to study finite-temperature electron liquid.[33] Later, the classical representation of quantum systems at equilibrium was applied to study two special quantum systems.[34−37] Recently, the restricted path integral Monte Carlo (PIMC) approach was adopted to numerically evaluate the XC free energy per particle of the warm dense UEG.[38] The numerical data were soon adopted by Karasiev and co-workers for the development of an accurate parametrization for FxcLDA[ρ].[39] Very recently, the most advanced density matrix quantum Monte Carlo method was adopted to resolve the discrepancy between the configuration and restricted PIMC results.[40] Since these simulations are only limited to small systems, the finite-size errors were further corrected from the ab initio finite-N quantum Monte Carlo calculations and compared with the previous results where the significant difference was found.[41] In addition, a careful study on the free energy functional for a noninteracting electron system at temperature T was recently performed by Dufty and Trickey.[42]

Nevertheless, owing to the local approximation, the LDA XC free energy functional FxcLDA[ρ] should still suffer from the SIE problem. To reduce this, the aforementioned LCH scheme can be extended to FT-DFT, i.e., incorporating the HF exchange free energy FxHF into the LDA XC free energy functional FxcLDA. The resulting LCH XC free energy functional can be expressed as FxcLCH = FxLR-HF + FxSR-LDA + FcLDA, where FxLR-HF is the HF exchange free energy associated with the LR operator (e.g., the erf operator), FxSR-LDA is the LDA exchange free energy functional associated with the SR operator (e.g., the erfc operator), and FcLDA is the LDA correlation free energy functional (i.e., associated with the Coulomb operator). At ω = 0 (i.e., the SR interelectronic interaction reduces to the Coulomb interelectronic interaction), FxSR-LDA reduces to the LDA exchange free energy functional FxLDA and hence FxcLCH reduces to FxcLDA. Note that in FT-DFT, FxLR-HF is well-defined and FcLDA is readily available in the literature. Therefore, the focus of the present work is a reliably accurate parametrization for FxSR-LDA. For brevity, hereafter, we adopt “the SR LDA exchange free energy functional” (FxSR-LDA) for “the LDA exchange free energy functional associated with the SR operator” and adopt “the SR LDA exchange energy functional” (ExSR-LDA or FxSR-LDA with T = 0) for “the LDA exchange energy functional associated with the SR operator”.

On the other hand, in 2012, Chai developed thermally assisted-occupation density functional theory (TAO-DFT)[43] for the study of the ground-state properties of nanoscale systems with strong static correlation effects (which are “challenging systems” for traditional electronic structure methods).[44−51] In contrast to KS-DFT, an auxiliary system of N noninteracting electrons at some “fictitious temperature” is employed in TAO-DFT, wherein strong static correlation is explicitly approximated by the entropy contribution (e.g., see eq 26 of ref (43)). A self-consistent scheme for determining the fictitious temperature in TAO-DFT was recently proposed.[52] In 2017, the GH and RSH schemes in TAO-DFT were also developed for a wide range of applications.[53] Relative to local and semilocal density functionals in TAO-DFT, GH functionals in TAO-DFT were shown to possess reduced SIEs. However, to employ the RSH scheme in TAO-DFT, the SR (or LR) exchange free energy functional (at the fictitious temperature) should be further developed. Therefore, the development of SR LDA exchange free energy functional can also be the first step toward finding generally accurate RSH functionals in TAO-DFT, highlighting the value of the present work.

For the reasons given above, there is a strong need for developing the SR LDA exchange free energy functional for the RSH schemes in both FT-DFT and TAO-DFT. In Results and Discussion, we first examine analytical expressions for the SR exchange free energy per particle of the UEG at the low- and high-temperature limits. For the low-temperature limit, the first-order dependence on temperature is found to be absent, similar to that found for the exchange free energy per particle of the UEG. Moreover, the zero-temperature limit agrees exactly with that reported in the literature.[54] Based on these limits and our findings, we develop a reliably accurate analytical parametrization for the SR exchange free energy of the UEG as a function of the uniform electron density, temperature, and range-separation parameter, retaining the correct zero-temperature and high-temperature limits. Besides, the SR exchange potential of the UEG, obtained directly from the functional derivative of the parametrized SR exchange free energy of the UEG, is reliably accurate relative to the corresponding numerical value. For a general system, the SR LDA exchange free energy functional is developed. In the last section, we give our conclusions.

Results and Discussion

SR Exchange Free Energy Per Particle of the UEG

Consider a spin-unpolarized system containing N electrons associated with the SR interelectronic interaction (i.e., the erfc operator given by eq with the range-separation parameter ω) in a volume V at temperature T, with a positive background charge keeping the system neutral. Here, N and V are taken to infinity in the manner that keeps the electron density (ρ = N/V) finite. The SR exchange free energy of the UEG, FxSR, can be expressed in both momentum space and coordinate space[30,31]where n is given by the Fermi–Dirac distribution function (with β = 1/(kBT) and k = |k|)and G(x) (with x ≡ |x|) is defined asOn the basis of eq , one can replace the chemical potential μ with the uniform electron density ρ (i.e., the inversion of the equation below can be performed, as shown in refs (30) and (55))and numerically evaluate the SR exchange free energy per particle of the UEGas a function of the Fermi wave vector kF = (3π2ρ)1/3 and two dimensionless variables (t ≡ kBT/EF and λ ≡ ω/kF). Here, EF = kF2/2 is the Fermi energy. It is convenient to define the scaled SR exchange free energy per particle of the UEGwhich is a function of the dimensionless variables t and λ only. Here, μx(kF, 0) = −kF/π is the exchange potential of the UEG at the Fermi wave vector kF and t = 0.

In the following discussion, for the given values of ρ and ω, we examine the SR exchange free energy per particle of the UEG at the low-t and high-t limits. Therefore, the low-t and high-t limits discussed are equivalent to the low-temperature (low-T) and high-temperature (high-T) limits, respectively.

Low-Temperature Limit

Similar to the procedures of Horovitz and Thieberger,[31] we first express FxSR [given by eq ] as a function of z ≡ eβμwhere the lower limit of integration is determined by direct calculation of eq in the limit z → 0The derivative of FxSR with respect to z is given bywhere g(x) is defined asWe now proceed to evaluate FxSR [given by eq ] at the large-z limit, which corresponds to the low-t limit (e.g., see eq 5 of ref (55)). It can be shown thatwhere h(z) is a function of z only. Therefore, the integration of eq with respect to z from 0 to a constant z1 (see ref (31) for the reason that one can choose such a constant) yields at most a second-order temperature-dependent term, which can be expressed as δ/β2. Accordingly, one can separate the integration in eq with respect to z′ into two parts: (i) from 0 to z1 and (ii) from z1 to z. Using the generalized Sommerfeld’s lemma,[56] in the range of z > z1, one obtains the following low-temperature expansionwhere the variable is rescaled to . The integration in eq with respect to z′ from z1 to z can be transformed into integration with respect to k from to , which yields the SR exchange free energy per volume of the UEG at the low-temperature limithereandwhere the exponential integral Ei(x) is defined asThe first two terms of eq are ω-independent and are the same as the results of Horovitz and Thieberger.[31] The last two terms (i.e., those with I(k) and J(k)) are ω-dependent, which are the new results brought by the SR interelectronic interaction. From eq , we find that the first-order temperature-dependent term vanishes, which is similar to that found from the exchange free energy per volume of the UEG at the low-temperature limit.[31] Note also that in the last term of eq , J(k1)/β2, besides an analogue term with the logarithmic term of T2 ln T in ref (31) we arrive at an additional term of T2 Ei(−T), which originates from the SR interelectronic interaction. For an analysis of the divergent term, one may refer to ref (57). The dependence on k (i.e., (2μ)1/2) can be replaced by the uniform electron density ρ (see eq 12 of ref (31)), as T approaches zero.

In addition, as ω → 0 (i.e., the SR interelectronic interaction reduces to the Coulomb interelectronic interaction), eq correctly reduces to the exchange free energy per volume of the UEG at the low-temperature limitwhich was previously reported by Horovitz and Thieberger.[31] Note that the first-order temperature-dependent term in eq vanishes, which later guided the parametrization function for the LDA exchange free energy functional of Perrot and Dharma-wardana.[30]

Besides, on the basis of eq , as 1/β → 0 (i.e., as T or t reduces to zero), the scaled SR exchange free energy per particle of the UEG at the zero-temperature limit (first derived by Gill, Adamson, and Pople[54]) can be correctly obtained

High-Temperature Limit

Here, we examine the SR exchange free energy per particle of the UEG at the high-temperature limit. The Taylor series expansion of g(x) [see eq ] with respect to z ≡ eβμ yieldsAt the small-z limit, which corresponds to the high-t limit (e.g., see eq 5 of ref (55)), one can keep only the first term and obtainConsequently, the resulting FxSR can be expressed asUsing eq 2.3 of ref (30) to replace z with ρ at the high-temperature limit, the scaled SR exchange free energy per particle of the UEG at the high-temperature limit[58] can be obtainedwhere t = h stands for the high-temperature limit.

Note that as ω → 0 or λ → 0 (i.e., the SR interelectronic interaction reduces to the Coulomb interelectronic interaction), eq correctly reduces to the scaled exchange free energy per particle of the UEG at the high-temperature limit, i.e., 1/(3t).[30]

Parametrization for the SR Exchange Free Energy Per Particle of the UEG

In the work of Perrot and Dharma-wardana,[30] a fitting function was proposed to parametrize the numerical data of the exchange free energy per particle of the UEG. As the first-order temperature-dependent term vanishes at the low-temperature limit,[31] the temperature-dependent term in the fitting function of Perrot and Dharma-wardana starts from the t2 term (see eq 3.2 of ref (30)).

In the present work, to incorporate the correct zero-temperature limit f̃xSR(t = 0, λ) [see eq ] and high-temperature limit f̃xSR(t = h, λ) [see eq ], and also our findings that the first-order temperature-dependent term vanishes (i.e., the temperature-dependent term starts from the t2 term) at the low-temperature limit, we propose the following fitting functionto parametrize the numerical data of the scaled SR exchange free energy per particle of the UEG f̃xSR(t, λ) [given by eq ], where x(λ) (i = 1, 2, 3, and 4) is defined asand y(λ) (i = 1 and 2) is defined asThe fitting to the numerical data is performed in the range of 0 < t < 12 and 0 < λ < 20. The optimized parameters for x(λ) and y(λ) are shown in Tables and 2, respectively.

Table 1

Optimized Parameters for x(λ) (i = 1, 2, 3, and 4) [See Equation ]

i	1	2	3	4
ci0	0.3729	5.6674	0.2127	16.0023
ci1	0.0051	0.1777	0.0036	0.0894
ci2	0.0438	0.7474	0.0258	7.1526
ci3	0.3485	0.3471	0.2023	14.8795
ci4	2.6256	2.7513	1.9715	2.0447

Table 2

Optimized Parameters for y(λ) (i = 1 and 2) [See Equation ]

i	1	2
di0	0.2839	0.7331
di1	0.9912	0.9973
di2	–0.4510	–0.1511
di3	0.4941	0.6022
di4	1.0075	1.7642

Figure shows the surface plot for the parametrization of f̃xSR(t, λ) [given by eq ] and the scattered circles for the numerical data of f̃xSR(t, λ) [given by eq ]. Figure shows the fitting curves for the parametrization of f̃xSR(t, λ) [given by eq ] and the scattered circles for the numerical data of f̃xSR(t, λ) [given by eq ] at various λ values.

Figure 1

Scaled SR exchange free energy per particle of the UEG, f̃xSR(t, λ), as a function of t and λ. Surface: parametrization [given by eq ]. Circles: numerical data [given by eq ].

Figure 2

Scaled SR exchange free energy per particle of the UEG, f̃xSR(t, λ), as a function of t and λ. Lines: parametrization [given by eq ]. Circles: numerical data [given by eq ]. Here, magenta, red, yellow, cyan, green, and blue correspond to λ = 0.0, 0.2, 0.4, 0.6, 0.8, and 1.0, respectively.

To examine the accuracy of our parametrization, Figure shows the relative error of the parametrization of f̃xSR(t, λ) [given by eq ], where the relative error is defined as the absolute value of ((parametrization [given by eq ] – numerical data [given by eq ])/(numerical data [given by eq ])). Relative to the numerical data, our parametrization is reliably accurate. The relative error is vanishingly small in the low-t and high-t regions. However, in the intermediate-t region, the relative error is slightly larger, especially for the larger λ. The maximum relative error is 0.087 (i.e., the maximum percentage error = 8.7%) at t = 6 and λ = 20. Therefore, further investigation on the expression of f̃xSR(t, λ) at the large-λ limit can be essential for improved parametrization.

Figure 3

Relative error of the parametrization of the scaled SR exchange free energy per particle of the UEG, f̃xSR(t, λ), as a function of t and λ. Here, the relative error is defined as the absolute value of ((parametrization [given by eq ] – numerical data [given by eq ])/(numerical data [given by eq ])).

SR Exchange Potential of the UEG

In the work of Perrot and Dharma-wardana,[30] the exchange potential of the UEG was parametrized separately, which can, however, be inconsistent with the functional derivative of their parametrized exchange free energy functional of the UEG. It is worth mentioning that Karasiev and co-workers[55] recently proposed a more accurate parametrization (compared to the one from ref (30)) for the UEG exchange free energy functional with exchange potential calculated as the corresponding functional derivative.

For consistency, in the present work, the SR exchange potential of the UEG, μxSR, is obtained directly from the functional derivative of FxSRSubstituting fxSR = fxSR(kF, t, λ) = μx(kF, 0) f̃xSR(t, λ) = (−kF/π) f̃xSR(t, λ) [given by eq ] into eq , μxSR can be expressed asOn the basis of eq , it is convenient to define the scaled SR exchange potential of the UEGwhich is a function of the dimensionless variables t and λ only. As f̃xSR(t, λ) [given by eq ] is parametrized, μ̃xSR(t, λ) [given by eq ] can be evaluated analytically.

Figure shows the surface plot for μ̃xSR(t, λ) [given by eq ] and the scattered circles for the corresponding numerical data [given by differentiating eq ]. Figure shows the fitting curves for μ̃xSR(t, λ) [given by eq ] and the scattered circles for the corresponding numerical data [given by differentiating eq ].

Figure 4

Scaled SR exchange potential of the UEG, μ̃xSR(t, λ), as a function of t and λ. Surface: parametrization [given by eq ]. Circles: numerical data [given by differentiating eq ].

Figure 5

Scaled SR exchange potential of the UEG, μ̃xSR(t, λ), as a function of t and λ. Lines: parametrization [given by eq ]. Circles: numerical data [given by differentiating eq ]. Here, magenta, red, yellow, cyan, green, and blue correspond to λ = 0.0, 0.2, 0.4, 0.6, 0.8, and 1.0, respectively.

To assess the accuracy of our parametrization, Figure shows the relative error of the parametrization of μ̃xSR(t, λ) [given by eq ], where the relative error is defined as the absolute value of ((parametrization [given by eq ] – numerical data [given by differentiating eq ])/(numerical data [given by differentiating eq ])). A good agreement between our parametrization and the corresponding numerical data can be clearly seen from the figure. The relative error is vanishingly small in the low-t and high-t regions. Nonetheless, in the intermediate-t region, the relative error is slightly larger, especially for the larger λ. The maximum relative error is 0.108 (i.e., the maximum percentage error = 10.8%) at t = 7.7 and λ = 20. As mentioned previously, the accuracy of the parametrization may be further improved by investigating the expression of f̃xSR(t, λ) (and hence the corresponding expression of μ̃xSR(t, λ) given by eq ) at the large-λ limit.

Figure 6

Relative error of the parametrization of the scaled SR exchange potential of the UEG, μ̃xSR(t, λ), as a function of t and λ. Here, the relative error is defined as the absolute value of ((parametrization [given by eq ] – numerical data [given by differentiating eq ])/(numerical data [given by differentiating eq ])).

LDA for the SR Exchange Free Energy Functional

In the previous section, the SR exchange free energy per particle of the UEG fxSR(t, λ) has been discussed and parametrized, which can now be extended to a general system.

Consider a spin-unpolarized system containing N electrons associated with the SR interelectronic interaction (i.e., the erfc operator given by eq with the range-separation parameter ω) at temperature T, in the presence of an external potential vext(r). The LDA for the SR exchange free energy per particle can be obtained by replacing the uniform electron density ρ in eq with the local electron density ρ(r). Accordingly, kF, EF, t, λ, and μx(kF, 0) are replaced with kF(r) = [3π2ρ(r)]1/3, EF(r) = [kF(r)]2/2, t(r) = kBT/EF(r), λ(r) = ω/kF(r), and μx(kF(r), 0) = −kF(r)/π, respectively. Consequently, the SR LDA exchange free energy functional can be expressed aswhere Cx = −(3/π)1/3 and f̃xSR(t(r), λ(r)) [given by eq ] is the scaled SR exchange free energy per particle of the UEG at the local electron density ρ(r). The SR LDA exchange potential is given by the functional derivative of FxSR-LDA[ρ]Owing to the spin-scaling relation,[59] the extension of the SR LDA exchange free energy functional to a spin-polarized system (i.e., with the α-spin density ρα(r), β-spin density ρβ(r), temperature T, and range-separation parameter ω) is straightforwardwhere the spin-polarized functional FxSR-LDA[ρα, ρβ] [see eq ] can be conveniently expressed by the spin-unpolarized functional FxSR-LDA[ρ] [see eq ].

Conclusions

In summary, we have examined analytical expressions for the SR exchange free energy per particle of the UEG at the low- and high-temperature limits. The SR interelectronic interaction brings extra terms in the two limiting forms when compared with those for the Coulomb interelectronic interaction. At the low-temperature limit, the temperature-dependent term starts from the t2 term, which is similar to that found for the exchange free energy per particle of the UEG. An analytical fitting function has been proposed for the SR exchange free energy per particle of the UEG. Accordingly, the SR LDA exchange free energy functional for a general system has been developed, with which RSH functionals can be readily devised in both FT-DFT and TAO-DFT.

In the future, we plan to develop the SR exchange free energy functionals based on more sophisticated density functional approximations (e.g., GGAs) to further improve the accuracy of RSH functionals in both FT-DFT and TAO-DFT. Note that an accurate GGA XC free energy functional (i.e., associated with the Coulomb operator) has been recently developed.[60]