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 quantummechanical methods for studying the ground-state properties of atoms, molecules, and bulk materials. 3,4 However, since the exact exchange−correlation (XC) energy functional E xc [ρ] 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, E xc LDA , 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 E x HF into the LDA XC energy functional E xc LDA . 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 longrange 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 1, 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 E xc LCH = E x LR-HF + E x SR-LDA + E c LDA , where E x LR-HF is the HF exchange energy associated with the LR operator, E x SR-LDA is the LDA exchange energy functional associated with the SR operator, and E c LDA 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), E xc LCH reduces to E xc LDA . Note also that in the LCH scheme, E x SR-LDA and E c LDA , 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 equations is the effective potential, v ext (r) is the external potential, and μ xc (r) = δF xc [ρ]/δρ(r) is the XC potential (i.e., the functional derivative of the XC free energy functional F xc [ρ]). The thermal equilibrium density is given by where μ is the chemical potential chosen to conserve the number of electrons and k B is the Boltzmann constant. At T = 0, FT-DFT reduces to KS-DFT and the XC free energy functional F xc [ρ] reduces to the XC energy functional E xc [ρ]. 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 F xc [ρ] has been commonly approximated by the "zerotemperature" XC energy functional E xc [ρ] evaluated with the thermal equilibrium density ρ(r) at temperature T.
In 1984, an analytical parametrization for the LDA XC free energy functional F xc LDA [ρ] 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 F xc LDA [ρ]. 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 finitesize 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 F xc LDA [ρ] 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 F x HF into the LDA XC free energy functional F xc LDA . The resulting LCH XC free energy functional can be expressed as F xc LCH = F x LR-HF + F x SR-LDA + F c LDA , where F x LR-HF is the HF exchange free energy associated with the LR operator (e.g., the erf operator), F x SR-LDA is the LDA exchange free energy functional associated with the SR operator (e.g., the erfc operator), and F c LDA 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), F x SR-LDA reduces to the LDA exchange free energy functional F x LDA and hence F xc LCH reduces to F xc LDA . Note that in FT-DFT, F x LR-HF is well-defined and F c LDA is readily available in the literature. Therefore, the focus of the present work is a reliably accurate parametrization for F x SR-LDA . For brevity, hereafter, we adopt "the SR LDA exchange free energy functional" (F x SR-LDA ) for "the LDA exchange free energy functional associated with the SR operator" and adopt "the SR LDA exchange energy functional" (E x SR-LDA or F x SR-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

ACS Omega
Article 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 1 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, F x SR , can be expressed in both momentum space and coordinate space 30,31 where n k is given by the Fermi−Dirac distribution function and On the basis of eq 5, 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 UEG as a function of the Fermi wave vector k F = (3π 2 ρ) 1/3 and two dimensionless variables (t ≡ k B T/E F and λ ≡ ω/k F ). Here, E F = k F 2 /2 is the Fermi energy. It is convenient to define the scaled SR exchange free energy per particle of the UEG which is a function of the dimensionless variables t and λ only.
Here, μ x (k F , 0) = −k F /π is the exchange potential of the UEG at the Fermi wave vector k F 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 F x SR [given by eq 5] as a function of z ≡ e βμ (11) where the lower limit of integration is determined by direct calculation of eq 5 in the limit z → 0 i k j j j j j j j y The derivative of F x SR with respect to z is given by where g(x) is defined as We now proceed to evaluate F x SR [given by eq 11] at the large-z limit, which corresponds to the low-t limit (e.g., see eq 5 of ref 55). It can be shown that where h(z) is a function of z only. Therefore, the integration of eq 13 with respect to z from 0 to a constant z 1 (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 11 with respect to z′ into two parts: (i) from 0 to z 1 and (ii) from z 1 to z. Using the generalized Sommerfeld's lemma, 56 in the range of z > z 1 , one obtains the following low-temperature expansion where the variable is rescaled to β = k z 2 ln / . The integration in eq 11 with respect to z′ from z 1 to z can be transformed into integration with respect to k from β = k z 2 ln / 1 1 to μ 2 , which yields the SR exchange free energy per volume of the UEG at the low-temperature limit

ACS Omega
and i k j j j j y where the exponential integral Ei(x) is defined as The first two terms of eq 17 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 17, 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 17, J(k 1 )/β 2 , besides an analogue term with the logarithmic term of T 2 ln T in ref 31 we arrive at an additional term of T 2 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 17 correctly reduces to the exchange free energy per volume of the UEG at the low-temperature limit Ä which was previously reported by Horovitz and Thieberger. 31 Note that the first-order temperature-dependent term in eq 21 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 17, 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 hightemperature limit. The Taylor series expansion of g(x) [see eq 14] with respect to z ≡ e βμ yields At 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 obtain Consequently, the resulting F x SR can be expressed as Using eq 2.3 of ref 30 to replace z with ρ at the hightemperature limit, the scaled SR exchange free energy per particle of the UEG at the high-temperature limit 58 can be obtained x SR F x F x SR 2 (26) where 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 26 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 Dharmawardana, 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 t 2 term (see eq 3.2 of ref 30).
In the present work, to incorporate the correct zerotemperature limit fx SR (t = 0, λ) [see eq 22] and hightemperature limit fx SR (t = h, λ) [see eq 26], and also our findings that the first-order temperature-dependent term vanishes (i.e., the temperature-dependent term starts from the t 2 term) at the low-temperature limit, we propose the following fitting function

ACS Omega
Article to parametrize the numerical data of the scaled SR exchange free energy per particle of the UEG fx SR (t, λ) [given by eq 10], where x i (λ) (i = 1, 2, 3, and 4) is defined as and y i (λ) (i = 1 and 2) is defined as The fitting to the numerical data is performed in the range of 0 < t < 12 and 0 < λ < 20. The optimized parameters for x i (λ) and y i (λ) are shown in Tables 1 and 2, respectively. Figure 1 shows the surface plot for the parametrization of fx SR (t, λ) [given by eq 27] and the scattered circles for the numerical data of fx SR (t, λ) [given by eq 10]. Figure 2 shows the fitting curves for the parametrization of fx SR (t, λ) [given by eq 27] and the scattered circles for the numerical data of fx SR (t, λ) [given by eq 10] at various λ values.
To examine the accuracy of our parametrization, Figure 3 shows the relative error of the parametrization of fx SR (t, λ) [given by eq 27], where the relative error is defined as the absolute value of ((parametrization [given by eq 27] − numerical data [given by eq 10])/(numerical data [given by eq 10])). 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,

ACS Omega
Article 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̃x SR (t, λ) at the large-λ limit can be essential for improved parametrization. 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 (30) Substituting x On the basis of eq 31, it is convenient to define the scaled SR exchange potential of the UEG x SR F x F x SR x SR x SR x SR (32) Figure 3. Relative error of the parametrization of the scaled SR exchange free energy per particle of the UEG, f̃x SR (t, λ), as a function of t and λ. Here, the relative error is defined as the absolute value of ((parametrization [given by eq 27] − numerical data [given by eq 10])/(numerical data [given by eq 10])).

ACS Omega
Article which is a function of the dimensionless variables t and λ only. As fx SR (t, λ) [given by eq 27] is parametrized, μx SR (t, λ) [given by eq 32] can be evaluated analytically. Figure 4 shows the surface plot for μx SR (t, λ) [given by eq 32] and the scattered circles for the corresponding numerical data [given by differentiating eq 10]. Figure 5 shows the fitting curves for μx SR (t, λ) [given by eq 32] and the scattered circles for the corresponding numerical data [given by differentiating eq 10].
To assess the accuracy of our parametrization, Figure 6 shows the relative error of the parametrization of μ̃x SR (t, λ) [given by eq 32], where the relative error is defined as the absolute value of ((parametrization [given by eq 32] − numerical data [given by differentiating eq 10])/(numerical data [given by differentiating eq 10])). 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 fx SR (t, λ) (and hence the corresponding expression of μx SR (t, λ) given by eq 32) at the large-λ limit.
LDA for the SR Exchange Free Energy Functional. In the previous section, the SR exchange free energy per particle of the UEG f x SR (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 1 with the range-separation parameter ω) at temperature T, in the presence of an external potential v ext (r). The LDA for the SR exchange free energy per particle can be obtained by replacing the uniform electron density ρ in eq 27 with the local electron density ρ(r). Accordingly, k F , E F , t, λ, and μ x (k F , 0) are replaced with k F (r) = [3π 2 ρ(r)] 1/3 , E F (r) = [k F (r)] 2 /2, t(r) = k B T/E F (r), λ(r) = ω/k F (r), and μ x (k F (r), 0) = −k F (r)/π, respectively. Consequently, the SR LDA exchange free energy functional can be expressed as r r r r r r r where C x = −(3/π) 1/3 and fx SR (t(r), λ(r)) [given by eq 27] is the scaled SR exchange free energy per particle of the UEG at the local electron density ρ(r  (34) 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 straightforward

■ CONCLUSIONS
In summary, we have examined analytical expressions for the SR exchange free energy per particle of the UEG at the lowand 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 t 2 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.  Here, the relative error is defined as the absolute value of ((parametrization [given by eq 32] − numerical data [given by differentiating eq 10])/(numerical data [given by differentiating eq 10])).

Article
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