Improvement of functionals in density functional theory by the inverse Kohn--Sham method and density functional perturbation theory

TL;DR

This paper introduces a novel approach combining inverse Kohn-Sham and perturbation theory to enhance energy density functionals in density functional theory, achieving accurate reproduction of exchange-correlation functionals and ground-state properties.

Contribution

It presents a new method to improve energy density functionals by treating their difference as a perturbation, validated through benchmark calculations on noble-gas atoms.

Findings

Reproduced exchange-correlation functionals with high accuracy.

Achieved precise ground-state energies and densities.

Validated method on noble-gas atom systems.

Abstract

We propose a way to improve energy density functionals (EDFs) in the density functional theory based on the combination of the inverse Kohn--Sham method and the density functional perturbation theory. Difference between the known EDF and the exact one is treated as the first-order perturbation. As benchmark calculations, we reproduce the theoretical exchange and correlation functionals in the local density approximation. Systems of noble-gas atoms are used for benchmark calculations, and the ground-state energies and densities, as well as the functionals, are reproduced with good accuracies.

Tables

Table 1. Table 1: The coefficients λ ​ A n 𝜆 subscript 𝐴 𝑛 \lambda A_{n} and α n subscript 𝛼 𝑛 \alpha_{n} and the ground-state energies E gs n th superscript subscript 𝐸 gs n th E_{\scriptstyle{\text{{{{gs}}}}}}^{\scriptstyle{\text{{{{$n$th}}}}}} calculated in the n 𝑛 n th iteration shown in Eq.  ( 5n \@@italiccorr ) for the pair of atoms He He \mathrm{He} and Ne Ne \mathrm{Ne} . The Hartree functional given in Eq.  ( LABEL:eq:1stCase_tilde_EHxc \@@italiccorr ) is used for E ~ Hxc 0 th ​ [ ρ ] subscript superscript ~ 𝐸 0 th Hxc delimited-[] 𝜌 \tilde{E}^{\scriptstyle{\text{{{{$0$th}}}}}}_{\scriptstyle{\text{{{{Hxc}}}}}}\left[\rho\right] and the Hartree plus LDA exchange functional given in Eq.  ( LABEL:eq:1stCase_EHxc \@@italiccorr ) are used for E Hxc target ​ [ ρ ] subscript superscript 𝐸 target Hxc delimited-[] 𝜌 E^{\scriptstyle{\text{{{{target}}}}}}_{\scriptstyle{\text{{{{Hxc}}}}}}\left[\rho\right] . All units are in the Hartree atomic unit.

$n$	${\alpha }_n$	$\lambda A_n$	$E_{\text{gs}}^{n\text{th}}$ of $\mathrm{He}$	$E_{\text{gs}}^{n\text{th}}$ of $\mathrm{Ne}$
$0$			$-1.9517070$	$-116.99029$
$1$	$1.3199872$	$-0.7920448$	$-2.8010654$	$-127.95544$
$2$	$1.0125782$	$0.0459594$	$-2.7115372$	$-127.49446$
Target	$1.3333333$	$-0.7385588$	$-2.7237069$	$-127.49109$

Table 2. Table 2: Same as Table 1 but for the pair of atoms Xe Xe \mathrm{Xe} and Rn Rn \mathrm{Rn} .

$n$	${\alpha }_n$	$\lambda A_n$	$E_{\text{gs}}^{n\text{th}}$ of $\mathrm{Xe}$	$E_{\text{gs}}^{n\text{th}}$ of $\mathrm{Rn}$
$0$			$-7054.6485$	$-21479.344$
$1$	$1.3311445$	$-0.7558229$	$-7224.9365$	$-21852.010$
$2$	$1.0436323$	$0.0306234$	$-7223.0601$	$-21848.894$
Target	$1.3333333$	$-0.7385588$	$-7223.1853$	$-21848.954$

Table 3. Table 3: Same as Table  1 but the Hartree plus LDA exchange functional given in Eq.  ( LABEL:eq:2ndCase_tilde_EHxc \@@italiccorr ) is used for E ~ Hxc 0 th ​ [ ρ ] superscript subscript ~ 𝐸 Hxc 0 th delimited-[] 𝜌 \tilde{E}_{\scriptstyle{\text{{{{Hxc}}}}}}^{\scriptstyle{\text{{{{$0$th}}}}}}\left[\rho\right] and the Hartree plus LDA exchange-correlation functional given in Eq.  ( LABEL:eq:2ndCase_EHxc \@@italiccorr ) is used for E Hxc target ​ [ ρ ] subscript superscript 𝐸 target Hxc delimited-[] 𝜌 E^{\scriptstyle{\text{{{{target}}}}}}_{\scriptstyle{\text{{{{Hxc}}}}}}\left[\rho\right] .

$n$	${\alpha }_n$	$\lambda A_n$	$E_{\text{gs}}^{n\text{th}}$ of $\mathrm{He}$	$E_{\text{gs}}^{n\text{th}}$ of $\mathrm{Ne}$
$0$			$-2.7237069$	$-127.49109$
$1$	$1.1093168$	$-0.0697561$	$-2.8365158$	$-128.22944$
$2$	$0.6951714$	$0.0000544$	$-2.8362598$	$-128.22881$
Target	PZ81	PZ81	$-2.8343506$	$-128.22766$

Table 4. Table 4: Same as Table 3 but for the pair of atoms Xe Xe \mathrm{Xe} and Rn Rn \mathrm{Rn} .

$n$	${\alpha }_n$	$\lambda A_n$	$E_{\text{gs}}^{n\text{th}}$ of $\mathrm{Xe}$	$E_{\text{gs}}^{n\text{th}}$ of $\mathrm{Rn}$
$0$			$-7223.1853$	$-21848.954$
$1$	$1.0862074$	$-0.0737520$	$-7228.4020$	$-21857.981$
Target	PZ81	PZ81	$-7228.3628$	$-21857.954$

Table 5. Table 5: The coefficients λ ​ A 1 𝜆 subscript 𝐴 1 \lambda A_{1} and α 1 subscript 𝛼 1 \alpha_{1} calculated in the first iteration for all the pairs of noble-gas atoms. All units are in the Hartree atomic unit.

Pairs	Exchange ${\alpha }_1$	Exchange $\lambda A_1$	Correlation ${\alpha }_1$	Correlation $\lambda A_1$
$\mathrm{He}$ and $\mathrm{Ne}$	$1.3199872$	$-0.7920448$	$1.1093169$	$-0.0697561$
$\mathrm{He}$ and $\mathrm{Ar}$	$1.3209765$	$-0.7926638$	$1.1049002$	$-0.0694846$
$\mathrm{Ne}$ and $\mathrm{Ar}$	$1.3235352$	$-0.7841588$	$1.0916279$	$-0.0718811$
$\mathrm{He}$ and $\mathrm{Kr}$	$1.3227758$	$-0.7937863$	$1.1022591$	$-0.0693213$
$\mathrm{Ne}$ and $\mathrm{Kr}$	$1.3263436$	$-0.7779323$	$1.0941254$	$-0.0715841$
$\mathrm{Ar}$ and $\mathrm{Kr}$	$1.3290958$	$-0.7658732$	$1.0956532$	$-0.0711566$
$\mathrm{He}$ and $\mathrm{Xe}$	$1.3235844$	$-0.7942892$	$1.1001295$	$-0.0691892$
$\mathrm{Ne}$ and $\mathrm{Xe}$	$1.3270817$	$-0.7762984$	$1.0919878$	$-0.0718384$
$\mathrm{Ar}$ and $\mathrm{Xe}$	$1.3292187$	$-0.7654719$	$1.0921353$	$-0.0717899$
$\mathrm{Kr}$ and $\mathrm{Xe}$	$1.3294148$	$-0.7644846$	$1.0848635$	$-0.0742007$
$\mathrm{He}$ and $\mathrm{Rn}$	$1.3244450$	$-0.7948236$	$1.0980694$	$-0.0690610$
$\mathrm{Ne}$ and $\mathrm{Rn}$	$1.3279028$	$-0.7744818$	$1.0905537$	$-0.0720085$
$\mathrm{Ar}$ and $\mathrm{Rn}$	$1.3297748$	$-0.7636589$	$1.0902537$	$-0.0721279$
$\mathrm{Kr}$ and $\mathrm{Rn}$	$1.3303022$	$-0.7606336$	$1.0856559$	$-0.0739752$
$\mathrm{Xe}$ and $\mathrm{Rn}$	$1.3311445$	$-0.7558229$	$1.0862074$	$-0.0737520$
Target	$1.3333333$	$-0.7385588$	PZ81	PZ81

Equations

$$E_{\text{gs}}=T_0\left[{\rho }_{\text{gs}}\right]+\int V_{\text{ext}}\left(\mathbf{r}\right){\rho }_{\text{gs}}\left(\mathbf{r}\right)d\mathbf{r}+E_{\text{H}}\left[{\rho }_{\text{gs}}\right]+E_{\text{xc}}\left[{\rho }_{\text{gs}}\right],$$

$$\left[-\frac{{\hslash }^2}{2m}{\nabla }^2+V_{\text{KS}}\left(\mathbf{r}\right)\right]{\psi }_i\left(\mathbf{r}\right)={\epsilon }_i{\psi }_i\left(\mathbf{r}\right),$$

$$V_{\text{KS}}\left(\mathbf{r}\right)=V_{\text{ext}}\left(\mathbf{r}\right)+\frac{\delta E_{\text{Hxc}}\left[{\rho }_{\text{gs}}\right]}{\delta \rho \left(\mathbf{r}\right)}.$$

$$E_{\text{Hxc}}^{\text{exact}}\left[\rho \right]={\tilde{E}}_{\text{Hxc}}\left[\rho \right]+\lambda E_{\text{Hxc}}^{(1)}\left[\rho \right]+O\left({\lambda }^2\right),$$

$${\rho }_{\text{gs}}^{(1)}\left(\mathbf{r}\right)=\sum _{i=1}^N\left[{\psi }_i^{(1)\ast }\left(\mathbf{r}\right){\psi }_i\left(\mathbf{r}\right)+{\psi }_i^{\ast }\left(\mathbf{r}\right){\psi }_i^{(1)}\left(\mathbf{r}\right)\right]$$

$$\int {\tilde{\psi }}_i^{\ast }\left(\mathbf{r}\right){\psi }_i^{(1)}\left(\mathbf{r}\right)d\mathbf{r}=0.$$

$$\tilde{h}=-\frac{{\hslash }^2}{2m}{\nabla }^2+V_{\text{ext}}\left(\mathbf{r}\right)+\frac{\delta {\tilde{E}}_{\text{Hxc}}\left[{\tilde{\rho }}_{\text{gs}}\right]}{\delta \rho \left(\mathbf{r}\right)}$$

$$\tilde{h}{\tilde{\psi }}_i={\tilde{\epsilon }}_i{\tilde{\psi }}_i.$$

$$E_{\text{Hxc}}\left[\rho \right]={\tilde{E}}_{\text{Hxc}}\left[\rho \right]+\lambda E_{\text{Hxc}}^{(1)}\left[\rho \right].$$

$$E_{\text{Hxc}}^{(1)}\left[\rho \right]=A\int {\left[\rho \left(\mathbf{r}\right)\right]}^{\alpha }d\mathbf{r},$$

$$\lambda A\int \left\{{\left[{\tilde{\rho }}_{\text{gs}}\left(\mathbf{r}\right)\right]}^{\alpha }+\left(\alpha -1\right){\left[{\rho }_{\text{gs}}^{\text{exact}}\left(\mathbf{r}\right)\right]}^{\alpha }\right\}d\mathbf{r}=C\left[{\rho }_{\text{gs}}^{\text{exact}}\right].$$

$$E_{\text{Hxc}}^{\text{nth}}\left[\rho \right]={\tilde{E}}_{\text{Hxc}}^{\text{0th}}\left[\rho \right]+\sum _{k=1}^n\lambda E_{\text{Hxc}}^{\text{(1), kth}}$$

$$E_{\text{Hxc}}^{\text{nth}}\left[\rho \right]={\tilde{E}}_{\text{Hxc}}^{\text{0th}}\left[\rho \right]+\sum _{k=1}^n\lambda A_k\int {\left[\rho \left(\mathbf{r}\right)\right]}^{{\alpha }_k}d\mathbf{r},$$

Full text

††: J. Phys. B: At. Mol. Opt. Phys.

1 Introduction

Density functional theory (DFT) is one of the most successful approaches to the calculation of the ground-state properties of the quantum many-body problems, such as atoms, molecules, solids [1, 2, 3, 4], and nuclear systems [5, 6]. Since for the same level of accuracy the numerical cost of DFT is much less than those of other quantum many-body methods [7], such as the quantum Monte Carlo method [8, 9, 10, 11, 12], many-body perturbation theory [13, 14, 15], configuration interaction method [16, 17, 18], coupled-cluster method [19, 20, 21, 22, 23], and transcorrelated method [24, 25, 26]. Thus, DFT is applicable to larger-scale calculations [27, 28, 29]. Furthermore, in principle, the DFT gives the exact ground-state density $\rho_{\scriptstyle{\text{{{{gs}}}}}}$ and energy $E_{\scriptstyle{\text{{{{gs}}}}}}$:

[TABLE]

where $T_{0}$ is the Kohn–Sham (KS) kinetic energy, $V_{\scriptstyle{\text{{{{ext}}}}}}$ is the external field, and $E_{\scriptstyle{\text{{{{H}}}}}}\left[\rho\right]$ and $E_{\scriptstyle{\text{{{{xc}}}}}}\left[\rho\right]$ are the Hartree and exchange-correlation energy density functionals (EDFs), respectively [1, 2]. However, in practice, the accuracy of the DFT calculation depends on the accuracy of the approximations for $E_{\scriptstyle{\text{{{{xc}}}}}}\left[\rho\right]$, as it is unknown.

In electron systems, many approximations for $E_{\scriptstyle{\text{{{{xc}}}}}}\left[\rho\right]$ have been proposed from first principles, i.e., non-empirically. The widely used ones are the local density approximation (LDA) [30, 31, 32, 33] and generalized gradient approximation (GGA) [34, 35, 36, 37]. In the LDA, $E_{\scriptstyle{\text{{{{xc}}}}}}\left[\rho\right]$ is approximated as a functional of local density $\rho\left({\bm{r}}\right)$, whereas in the GGA, $E_{\scriptstyle{\text{{{{xc}}}}}}\left[\rho\right]$ is approximated as a functional of $\rho\left({\bm{r}}\right)$ and its gradient $\left|\nabla\rho\left({\bm{r}}\right)\right|$. Approximations beyond the GGA have also been developed, such as hybrid functionals, and they are being applied extensively [38], in particular, in quantum chemistry. Note that even in the GGA, $E_{\scriptstyle{\text{{{{xc}}}}}}\left[\rho\right]$ is constructed based on the LDA one. Thus, the exchange interaction is not fully included, and hence the physics of localized electrons (including materials containing $d$- and $f$-electrons) is not captured accurately with semi-local functionals. In order to avoid these problems, phenomenological methods, such as the $\mbox{LDA}+U$ [39, 40, 41], $\mbox{LDA}+\mbox{DMFT}$ [42], exact-exchange [43, 44] methods, and hybrid functionals [45, 46, 47, 48] have been proposed. Functionals with the long-range correction due to the van der Waals interaction [49] and semi-empirical functionals [50, 51] have also been discussed for a long time.

In contrast, in nuclear systems, the exact form of the interaction in the vacuum between nucleons is still under discussion [52, 53, 54, 55, 56, 57, 58]. Even if the exact form of the interaction in the vacuum were known, the nuclear interaction in the medium is different due to its highly non-perturbative property [59, 60]. Therefore, it is still difficult to derive the Hartree-exchange-correlation EDF $E_{\scriptstyle{\text{{{{Hxc}}}}}}\left[\rho\right]=E_{\scriptstyle{\text{{{{H}}}}}}\left[\rho\right]+E_{\scriptstyle{\text{{{{xc}}}}}}\left[\rho\right]$ from first principles, although the Hartree–Fock calculation from the interaction in the vacuum has been discussed for a long time [61, 62, 63, 64, 65]. Thus, $E_{\scriptstyle{\text{{{{Hxc}}}}}}\left[\rho\right]$ for the nuclear interaction is treated phenomenologically [66, 67, 68, 69], with fitting parameters to experimental data [5]. Since the fitting parameters are usually determined from the experimental data of several stable nuclei, different parameter sets can give totally different results for exotic nuclei [70]. Comparisons between parameter sets are still being discussed [71, 72].

Hence, the derivation or construction of accurate EDFs is one of the primary goals in DFT for both electron and nuclear systems, whose strategies are, however, under debate [73, 74]. Recently, a new microscopic way to derive EDFs based on the functional renormalization group was suggested [75, 76, 77, 78], while it is not ready for realistic systems yet.

As an alternative way to improve EDFs, the inverse approach of DFT, the so-called inverse Kohn–Sham (IKS) method, was proposed in Refs. [79, 80]. As mentioned in Ref. [81], the KS potential $V_{\scriptstyle{\text{{{{KS}}}}}}\left({\bm{r}}\right)$ is unique concerning the system, and $V_{\scriptstyle{\text{{{{KS}}}}}}\left({\bm{r}}\right)$ is calculated from the given ground-state density $\rho_{\scriptstyle{\text{{{{gs}}}}}}\left({\bm{r}}\right)$ in the IKS. The information provided by the IKS, such as the single-particle energies $\varepsilon_{i}$, is expected to be valuable for improving the accuracy of EDFs [82, 83, 84]. Although with such an expectation, the actual way for using IKS has not been pointed out explicitly, only the improvement of the numerical methods of IKS has been discussed [85]. Since the density $\rho\left({\bm{r}}\right)$ is usually the quantity which one interests in, once $\rho\left({\bm{r}}\right)$ of the system is known, tasks of the calculation for the system are almost achieved.

Nevertheless, improvement of the EDFs by using the IKS is promising since the EDF is, in principle, unique for all the electron systems, while the density of some systems, for instance, atoms and light molecules, can be determined from experiments or high-accuracy calculations, such as the coupled-cluster and the configuration interaction methods.

Moreover, in nuclear systems, it is known that the effective nuclear force, i.e., the nuclear force in medium, is substantially different from the bare nuclear force [59] and its details are still unknown. Therefore, the nuclear EDFs are derived from the fitting to the experimental data [5, 6]. In contrast, for the light or medium-heavy nuclei, the many-body methods beginning from the bare nuclear force, so-called ab initio methods, such as the no-core shell model [86, 87], the self-consistent Green’s function method [88], the lattice chiral effective field theory [89, 90, 91], the in-medium similarity renormalization group [92, 93, 94, 95], and the Brueckner Hartree–Fock method [61, 62, 63, 64, 65, 96, 97, 98], are available. These ab initio methods provide accurate density distributions for light nuclei, but meanwhile, their calculations for heavy nuclei are impossible even in the visible future. Therefore, how to make the best use of these ab initio methods for developing nuclear DFT is a worldwide hot topic in the nuclear community [98, 99, 100].

The density functional perturbation theory (DFPT) has been developed for decades. The DFPT was originally developed for the calculations of phonon and response properties from the first principle [101, 102, 103, 104]. Properties of solids are determined from phonons as well as electrons [105], but the original DFT is applicable only for the electronic structure. In the DFPT, the displacement of the external field $V_{\scriptstyle{\text{{{{ext}}}}}}$ caused by the displacement of the nuclei or external electric field is treated by a combination of the perturbation and linear response theories [104]. In principle, each order of the DFPT can be derived systematically [106], and so far, the DFPT up to the third order has been discussed [103]. The DFPT has also been applied to derive the EDF [107, 108].

In order to attack the open question about the practical use of the IKS, for the first time, a new strategy based on the combination of the DFPT and the IKS, the so-called IKS-DFPT, is proposed in this paper. The first-order DFPT, which is also called the Hellmann–Feynman theorem [109], is used in this paper. The known functional is improved by using the IKS-DFPT. The ground-state energy derived by two methods: One is a combination of the original and inverse Kohn–Sham scheme, and the other is the DFPT. In the latter way, the difference between the “exact” functional and the known functional is assumed to be small enough, and the difference is treated as the perturbation. Note that the DFPT is used as the perturbation for the functional itself instead of the external potential.

As benchmark calculations, we verify this method by reproducing both the LDA exchange functional [30] and the LDA correlation functional [32]. The iteration of IKS-DFPT is also discussed.

This paper is organized as follows: First, the theoretical framework of the IKS-DFPT is given in Sec. 2. Then, the benchmark calculations and discussion of their results are shown in Sec. 3. Finally, Sec. 4 is devoted to the conclusion and perspectives.

2 Theoretical Framework

2.1 Formalism of IKS-DFPT

In the DFT, the ground-state density $\rho_{\scriptstyle{\text{{{{gs}}}}}}\left({\bm{r}}\right)$ and energy $E_{\scriptstyle{\text{{{{gs}}}}}}$ of an $N$-particle system are obtained by solving the KS equations self-consistently

[TABLE]

where $m$ is the mass of particles, $\psi_{i}\left({\bm{r}}\right)$ and $\varepsilon_{i}$ are the single-particle orbitals and energies, respectively, and $\rho_{\scriptstyle{\text{{{{gs}}}}}}\left({\bm{r}}\right)=\sum_{i=1}^{N}\left|\psi_{i}\left({\bm{r}}\right)\right|^{2}$. Here, $V_{\scriptstyle{\text{{{{KS}}}}}}\left({\bm{r}}\right)$ is the KS effective potential defined as

[TABLE]

The IKS provides this KS effective potential $V_{\scriptstyle{\text{{{{KS}}}}}}$ for each system from the ground-state density $\rho_{\scriptstyle{\text{{{{gs}}}}}}$. In the novel method IKS-DFPT, the conventional Hartree-exchange-correlation functional $\tilde{E}_{\scriptstyle{\text{{{{Hxc}}}}}}$, such as the PZ81 [32] and PBE [110, 111] functionals, will be improved by using the IKS.

Here, $\tilde{E}_{\scriptstyle{\text{{{{Hxc}}}}}}\left[\rho\right]$ is assumed to be close enough to the exact one $E^{\scriptstyle{\text{{{{exact}}}}}}_{\scriptstyle{\text{{{{Hxc}}}}}}\left[\rho\right]$, since $\tilde{E}_{\scriptstyle{\text{{{{Hxc}}}}}}\left[\rho\right]$ is known to work well. Hence, the difference between $E_{\scriptstyle{\text{{{{Hxc}}}}}}^{\scriptstyle{\text{{{{exact}}}}}}$ and $\tilde{E}_{\scriptstyle{\text{{{{Hxc}}}}}}$ is treated as a perturbation. If the difference is not small enough to be treated as the perturbation, the final results would be unreasonable.

In this paper, the first-order perturbation theory is used for the treatment of the difference between $E_{\scriptstyle{\text{{{{Hxc}}}}}}^{\scriptstyle{\text{{{{exact}}}}}}$ and $\tilde{E}_{\scriptstyle{\text{{{{Hxc}}}}}}$ as

[TABLE]

with a small parameter $\lambda$. Then, the exact single-particle orbitals $\psi^{\scriptstyle{\text{{{{exact}}}}}}_{i}\left({\bm{r}}\right)$, ground-state density $\rho^{\scriptstyle{\text{{{{exact}}}}}}_{\scriptstyle{\text{{{{gs}}}}}}\left({\bm{r}}\right)$, and energy $E^{\scriptstyle{\text{{{{exact}}}}}}_{\scriptstyle{\text{{{{gs}}}}}}$ are also expanded perturbatively:

[TABLE]

where quantities shown with the tilde are given by $\tilde{E}_{\scriptstyle{\text{{{{Hxc}}}}}}$. According to these definitions,

[TABLE]

is hold. Also, the first-order perturbation term $\psi_{i}^{(1)}$ is assumed to be orthogonal to $\tilde{\psi}_{i}$;

[TABLE]

The perturbation is assumed not to affect the external field, i.e., $V^{\scriptstyle{\text{{{{exact}}}}}}_{\scriptstyle{\text{{{{ext}}}}}}\left({\bm{r}}\right)=\tilde{V}_{\scriptstyle{\text{{{{ext}}}}}}\left({\bm{r}}\right)$. Moreover, $\rho^{\scriptstyle{\text{{{{exact}}}}}}_{\scriptstyle{\text{{{{gs}}}}}}\left({\bm{r}}\right)$ is assumed to be given, and thus $\psi^{\scriptstyle{\text{{{{exact}}}}}}_{i}\left({\bm{r}}\right)$ are calculated from the IKS.

Under these assumptions, we calculate $E^{\scriptstyle{\text{{{{exact}}}}}}_{\scriptstyle{\text{{{{gs}}}}}}$ in two different ways. One way is based on the first-order DFPT, and the other way is based on the IKS and KS equation. In the former way, substitution of Eqs. (4@@italiccorr), (LABEL:eq:idea_psi@@italiccorr), and (LABEL:eq:idea_rho@@italiccorr) into Eq. (1@@italiccorr) gives

[TABLE]

\process@envbody

where the Taylor expansion for general functionals $F\left[\rho\right]$

[TABLE]

\process@envbody

is used. Here, the kinetic term $T_{0}$ satisfies

[TABLE]

\process@envbody

Combining Eqs. (5f@@italiccorr), (LABEL:eq:gs_1@@italiccorr), and (LABEL:eq:T0_exp@@italiccorr), we get

[TABLE]

\process@envbody

It should be noted that the single-particle orbitals $\tilde{\psi}_{i}$ are the eigenstates of the single-particle Hamiltonian for the non-perturbative system

[TABLE]

as

[TABLE]

In the latter way, Eq. (4@@italiccorr) and integration of the KS equation (2@@italiccorr) gives $E^{\scriptstyle{\text{{{{exact}}}}}}_{\scriptstyle{\text{{{{gs}}}}}}$:

[TABLE]

\process@envbody

where $\varepsilon_{i}^{\scriptstyle{\text{{{{exact}}}}}}$ are obtained from $\rho^{\scriptstyle{\text{{{{exact}}}}}}_{\scriptstyle{\text{{{{gs}}}}}}$ by using the IKS. By comparing these two expressions of the ground-state energy and neglecting $O\left(\lambda^{2}\right)$ term, the equation for $E^{(1)}_{\scriptstyle{\text{{{{Hxc}}}}}}\left[\rho\right]$ is obtained:

[TABLE]

\process@envbody

The right-hand side of this equation can be calculated from the known quantities and its value depends only on the exact ground-state density $\rho_{\scriptstyle{\text{{{{gs}}}}}}^{\scriptstyle{\text{{{{exact}}}}}}$ and the known functional $\tilde{E}_{\scriptstyle{\text{{{{Hxc}}}}}}$. Thus, hereafter the right-hand side of the equation is shown as $C\left[\rho\right]$.

Finally, solving Eq. (LABEL:eq:basic_eq@@italiccorr), the Hartree-exchange-correlation functional in the IKS-DFPT in the first-order, i.e., the IKS-DFPT1, is derived as

[TABLE]

Because Eq. (LABEL:eq:basic_eq@@italiccorr) is a functional equation, it is difficult to be solved directly. In this work, we introduce one of the simplest ansatze for $E_{\scriptstyle{\text{{{{Hxc}}}}}}^{(1)}\left[\rho\right]$ within the LDA,

[TABLE]

which has the same form as the LDA exchange functional (cf. Eq. (LABEL:eq:1stCase_EHxc@@italiccorr)). Here, the values of $A$ and $\alpha$ are to be determined, and then we get

[TABLE]

To determine $A$ and $\alpha$, two systems, Systems 1 and 2, are required. Here, $\rho_{1}$ and $\rho_{2}$ are the exact ground-state densities, and $\tilde{\rho}_{1}$ and $\tilde{\rho}_{2}$ are the ground-state densities calculated from $\tilde{E}_{\scriptstyle{\text{{{{Hxc}}}}}}\left[\rho\right]$ of the Systems 1 and 2, respectively 111 If there are $n$ parameters in the ansatz for $E^{(1)}_{\scriptstyle{\text{{{{Hxc}}}}}}$, densities for $n$ systems are required to solve Eq. (5l@@italiccorr)-like algrebraic equation.. Substituting $\rho_{i}$ and $\tilde{\rho}_{i}$ ($i=1,\mskip 3.0mu2$) for Eq. (LABEL:eq:basic_eq@@italiccorr), it leads to the two equations for $\lambda A$ and $\alpha$. In such a way, $\lambda A$ and $\alpha$ can be determined. Note that in principle the Hartree-exchange-correlation EDF is system independent, and thus any system can be used as Systems 1 and 2. In this paper, the noble-gas atoms are used for two systems due to simplicity of the spherical symmetry.

In summary, the flowchart of the IKS-DFPT method is shown as Fig. 1.

2.2 Iteration of IKS-DFPT

On the one hand, a functional derived by the IKS-DFPT1 $E_{\scriptstyle{\text{{{{Hxc}}}}}}^{\scriptstyle{\text{{{{1st}}}}}}=\tilde{E}_{\scriptstyle{\text{{{{Hxc}}}}}}+\lambda E_{\scriptstyle{\text{{{{Hxc}}}}}}^{(1)}$ is assumed to be improved from the original functional $\tilde{E}_{\scriptstyle{\text{{{{Hxc}}}}}}$. On the other hand, the functional $E_{\scriptstyle{\text{{{{Hxc}}}}}}^{\scriptstyle{\text{{{{1st}}}}}}$ may still be able to be improved more, if we repeat the same procedure. For that, the functional $E_{\scriptstyle{\text{{{{Hxc}}}}}}^{\scriptstyle{\text{{{{1st}}}}}}$ is treated as the conventional functional $\tilde{E}_{\scriptstyle{\text{{{{Hxc}}}}}}$ above and applied the IKS-DFPT1 again. The functional derived by the IKS-DFPT again, $E_{\scriptstyle{\text{{{{Hxc}}}}}}^{\scriptstyle{\text{{{{2nd}}}}}}=E_{\scriptstyle{\text{{{{Hxc}}}}}}^{\scriptstyle{\text{{{{1st}}}}}}+\lambda E_{\scriptstyle{\text{{{{Hxc}}}}}}^{(1)}$, is expected to be improved.

To reach the most improved functional in the IKS-DFPT1 with the ansatz (5k@@italiccorr), the IKS-DFPT is applied to the derived functional iteratively.

Formally, at the $n$th step of the iteration, the Hartree-exchange-correlation functional calculated in the IKS-DFPT1 is

[TABLE]

where E_{\scriptstyle{\text{{{{Hxc}}}}}}^{\scriptstyle{\text{{{{(1), kth}}}}}} is the derived term in $k$th step. In particular, under the ansatz (5k@@italiccorr), E_{\scriptstyle{\text{{{{Hxc}}}}}}^{\scriptstyle{\text{{{{nth}}}}}}\left[\rho\right] is defined as

[TABLE]

where \tilde{E}_{\scriptstyle{\text{{{{Hxc}}}}}}^{\scriptstyle{\text{{{{0th}}}}}}\left[\rho\right] is the original one $\tilde{E}_{\scriptstyle{\text{{{{Hxc}}}}}}\left[\rho\right]$ at the first step. We perform the iteration until convergence. This indicates that we cannot improve the EDF further by using this scheme and ansatz (5k@@italiccorr).

3 Benchmark Calculations and Discussions

As benchmark calculations, $\rho^{\scriptstyle{\text{{{{target}}}}}}_{\scriptstyle{\text{{{{gs}}}}}}\left({\bm{r}}\right)$ is calculated from the theoretical $E^{\scriptstyle{\text{{{{target}}}}}}_{\scriptstyle{\text{{{{Hxc}}}}}}\left[\rho\right]$, instead of experimental data, and we test whether $E^{\scriptstyle{\text{{{{target}}}}}}_{\scriptstyle{\text{{{{Hxc}}}}}}\left[\rho\right]$ is reproduced in this scheme. In this section, the superscript “target” is used instead of “exact” since the functional which should be reproduced is already somehow approximated. All the pairs of the isolated noble-gas atoms ($\mathrm{He}$, $\mathrm{Ne}$, $\mathrm{Ar}$, $\mathrm{Kr}$, $\mathrm{Xe}$, and $\mathrm{Rn}$) are used as Systems 1 and 2. In calculations, the Hartree atomic unit is used. The ADPACK code [112] is used for the DFT calculations of the isolated atoms. Hereafter, $\lambda A_{n}$ is denoted by $A_{n}$.

We analyze two cases: 1) \tilde{E}_{\scriptstyle{\text{{{{Hxc}}}}}}^{\scriptstyle{\text{{{{0th}}}}}}\left[\rho\right] is the Hartree functional, and $E^{\scriptstyle{\text{{{{target}}}}}}_{\scriptstyle{\text{{{{Hxc}}}}}}\left[\rho\right]$ is the Hartree plus LDA exchange functional [30]:

[TABLE]

2) \tilde{E}_{\scriptstyle{\text{{{{Hxc}}}}}}^{\scriptstyle{\text{{{{0th}}}}}}\left[\rho\right] is the Hartree plus LDA exchange functional, and $E^{\scriptstyle{\text{{{{target}}}}}}_{\scriptstyle{\text{{{{Hxc}}}}}}\left[\rho\right]$ is the Hartree plus LDA exchange-correlation functional, where the PZ81 [32] functional $E_{\scriptstyle{\text{{{{c}}}}}}^{\scriptstyle{\text{{{{PZ81}}}}}}\left[\rho\right]$ is used:

[TABLE]

In both cases, the external field $V^{\scriptstyle{\text{{{{target}}}}}}_{\scriptstyle{\text{{{{ext}}}}}}\left({\bm{r}}\right)=\tilde{V}_{\scriptstyle{\text{{{{ext}}}}}}\left({\bm{r}}\right)$ is the Coulomb interaction between the nucleus and electron.

First, let us consider the first case, i.e., the calculations from the Hartree approximation (LABEL:eq:1stCase_tilde_EHxc@@italiccorr) to the Hartree–Fock–Slater (LDA exchange) approximation (LABEL:eq:1stCase_EHxc@@italiccorr). In Tables 1 and 2, the coefficients $\lambda A_{n}$ and $\alpha_{n}$ and the ground-state energies E_{\scriptstyle{\text{{{{gs}}}}}}^{\scriptstyle{\text{{{{nth}}}}}} calculated in the $n$th iteration are shown for the pairs of atoms $\mathrm{He}$-$\mathrm{Ne}$ and $\mathrm{Xe}$-$\mathrm{Rn}$, respectively. It is found that $\lambda A_{1}$ and $\alpha_{1}$ are obtained within $7.2\mskip 3.0mu\%$ and $1.0\mskip 3.0mu\%$ errors in $\mathrm{He}$-$\mathrm{Ne}$, and within $2.3\mskip 3.0mu\%$ and $0.2\mskip 3.0mu\%$ errors in $\mathrm{Xe}$-$\mathrm{Rn}$, respectively, from their target values. The heavier atoms reproduce the coefficients better. The results of the other pairs are shown in the Appendix.

The exchange energy density calculated in the first iteration, \varepsilon^{\scriptstyle{\text{{{{1st}}}}}}_{\scriptstyle{\text{{{{x}}}}}}\left(r_{\scriptstyle{\text{{{{s}}}}}}\right), and the ratio to the target one, \varepsilon^{\scriptstyle{\text{{{{1st}}}}}}_{\scriptstyle{\text{{{{x}}}}}}\left(r_{\scriptstyle{\text{{{{s}}}}}}\right)/\varepsilon^{\scriptstyle{\text{{{{target}}}}}}_{\scriptstyle{\text{{{{x}}}}}}\left(r_{\scriptstyle{\text{{{{s}}}}}}\right), are shown as functions of $r_{\scriptstyle{\text{{{{s}}}}}}$ in Fig. 2 for the pairs of $\mathrm{He}$-$\mathrm{Ne}$ and $\mathrm{Xe}$-$\mathrm{Rn}$ with dashed and dot lines, respectively, while the target one is shown with a solid line. Here, the energy density $\varepsilon_{i}\left(\rho\right)$ and the Wigner-Seitz radius $r_{\scriptstyle{\text{{{{s}}}}}}$ are defined as $E_{i}\left[\rho\right]=\intop\nolimits\varepsilon_{i}\left(\rho\right)\mskip 3.0mu\rho\left({\bm{r}}\right)\mskip 3.0mud{\bm{r}}$ ($i=\mathrm{x}$, $\mathrm{c}$) and $r_{\scriptstyle{\text{{{{s}}}}}}=\left[3/\left(4\pi\rho\right)\right]^{1/3}$, respectively. The pair of $\mathrm{Xe}$-$\mathrm{Rn}$ reproduces the target functional within a few percents in the range of $0.01\mskip 3.0mu\mathrm{a.u.}\leq r_{\scriptstyle{\text{{{{s}}}}}}\leq 100\mskip 3.0mu\mathrm{a.u.}$, which is generally better than the pair of $\mathrm{He}$-$\mathrm{Ne}$. As comparing to $\mathrm{He}$-$\mathrm{Ne}$ and $\mathrm{Xe}$-$\mathrm{Rn}$ cases, since the polynomial form of the functional in Eq. (5k@@italiccorr) is more sensitive to the high-density region, better reproduction in the high-density region leads to better reproduction of the coefficients.

For the iterations, it is found in Tables 1 and 2 that the difference between E_{\scriptstyle{\text{{{{gs}}}}}}^{\scriptstyle{\text{{{{nth}}}}}} and the target $E^{\scriptstyle{\text{{{{target}}}}}}_{\scriptstyle{\text{{{{gs}}}}}}$ becomes smaller as the iteration proceeds further. The ground-state energies of $\mathrm{He}$, $\mathrm{Ne}$, $\mathrm{Xe}$, and $\mathrm{Rn}$ are finally reproduced within $0.4\mskip 3.0mu\%$, $0.003\mskip 3.0mu\%$, $0.002\mskip 3.0mu\%$, and $0.0003\mskip 3.0mu\%$ errors, respectively, comparing with $28\mskip 3.0mu\%$, $8\mskip 3.0mu\%$, $2\mskip 3.0mu\%$, and $2\mskip 3.0mu\%$ errors at the zeroth step, E_{\scriptstyle{\text{{{{gs}}}}}}^{\scriptstyle{\text{{{{0th}}}}}}. This indicates the iteration helps the improvement of the ground-state energy.

The Wigner-Seitz radii r^{\scriptstyle{\text{{{{nth}}}}}}_{\scriptstyle{\text{{{{s}}}}}} calculated in the zeroth, first, and second iterations and the target one $r_{\scriptstyle{\text{{{{s}}}}}}^{\scriptstyle{\text{{{{target}}}}}}$ for $\mathrm{Rn}$ are shown as functions of $r$ in Fig. 3 with dot-dashed, dashed, dot, and solid lines, respectively. The ratio of calculated Wigner-Seitz radius to the target one, r^{\scriptstyle{\text{{{{nth}}}}}}_{\scriptstyle{\text{{{{s}}}}}}/r^{\scriptstyle{\text{{{{target}}}}}}_{\scriptstyle{\text{{{{s}}}}}}, for each step is also shown in the insert of Fig. 3. It is found that the ground-state density at the first step is already much improved, and it is even further improved as the iteration proceeds further. This indicates the iteration also helps the improvement of the ground-state density.

Next, let us consider the second case, i.e., the calculations from the Hartree–Fock–Slater (LDA exchange) approximation (LABEL:eq:2ndCase_tilde_EHxc@@italiccorr) to the LDA exchange-correlation (LABEL:eq:2ndCase_EHxc@@italiccorr). In Tables 3 and 4, the coefficients $\lambda A_{n}$ and $\alpha_{n}$ and the ground-state energies E^{\scriptstyle{\text{{{{nth}}}}}}_{\scriptstyle{\text{{{{gs}}}}}} calculated in the $n$th iteration are shown for the pairs of atoms $\mathrm{He}$-$\mathrm{Ne}$ and $\mathrm{Xe}$-$\mathrm{Rn}$, respectively. It is found that the ground-state energies are already reproduced well at the first step. Here, for the pair of $\mathrm{Xe}$-$\mathrm{Rn}$, the convergence reaches. For the pair of $\mathrm{He}$-$\mathrm{Ne}$, the ground-state energies are further improved slightly as the iteration proceeds further. The ground-state energies of $\mathrm{He}$, $\mathrm{Ne}$, $\mathrm{Xe}$, and $\mathrm{Rn}$ are finally reproduced within $0.07\mskip 3.0mu\%$, $0.0009\mskip 3.0mu\%$, $0.0005\mskip 3.0mu\%$, and $0.0001\mskip 3.0mu\%$ error, respectively, comparing with $4\mskip 3.0mu\%$, $0.6\mskip 3.0mu\%$, $0.07\mskip 3.0mu\%$, and $0.04\mskip 3.0mu\%$ errors at the zeroth step, E_{\scriptstyle{\text{{{{gs}}}}}}^{\scriptstyle{\text{{{{0th}}}}}}. The results of the other pairs are shown in the Appendix.

In order to compare the calculated correlation functionals and the target one, the correlation energy density calculated in the first iteration \varepsilon^{\scriptstyle{\text{{{{1st}}}}}}_{\scriptstyle{\text{{{{c}}}}}}\left(r_{\scriptstyle{\text{{{{s}}}}}}\right) is shown as functions of $r_{\scriptstyle{\text{{{{s}}}}}}$ in Fig. 4 for the pairs of $\mathrm{He}$-$\mathrm{Ne}$ and $\mathrm{Xe}$-$\mathrm{Rn}$ with dashed and dot lines, respectively, while the target one is shown with a solid line. The non-polynomial PZ81 functional is reproduced better in the lower-density region from the pair of $\mathrm{He}$-$\mathrm{Ne}$, and in the higher-density region from the pair of $\mathrm{Xe}$-$\mathrm{Rn}$, since heavier atoms have higher-density region than lighter atoms.

The Wigner-Seitz radii r^{\scriptstyle{\text{{{{nth}}}}}}_{\scriptstyle{\text{{{{s}}}}}} calculated in the zeroth and first iterations and the target one $r_{\scriptstyle{\text{{{{s}}}}}}^{\scriptstyle{\text{{{{target}}}}}}$ for $\mathrm{Rn}$ are shown as functions of $r$ in Fig. 5 with dot-dashed, dashed, and solid lines, respectively. The ratio of calculated Wigner-Seitz radius to the target one, r^{\scriptstyle{\text{{{{nth}}}}}}_{\scriptstyle{\text{{{{s}}}}}}/r^{\scriptstyle{\text{{{{target}}}}}}_{\scriptstyle{\text{{{{s}}}}}}, for each step is also shown in the insert of Fig. 5. It is found that the ground-state density at the first step is already much improved.

Comparing with the above two cases, in the first case, we note that the difference between \tilde{E}_{\scriptstyle{\text{{{{Hxc}}}}}}^{\scriptstyle{\text{{{{0th}}}}}}\left[\rho\right] and $E_{\scriptstyle{\text{{{{Hxc}}}}}}^{\scriptstyle{\text{{{{target}}}}}}\left[\rho\right]$ is larger, and thus either more iteration steps are required, or we should consider the IKS with the second-order DFPT. Nevertheless, as the LDA exchange functional is polynomial, the energy density $\varepsilon_{\scriptstyle{\text{{{{x}}}}}}$ is reproduced well in wide-density range. In the second case, the difference between \tilde{E}_{\scriptstyle{\text{{{{Hxc}}}}}}^{\scriptstyle{\text{{{{0th}}}}}}\left[\rho\right] and $E_{\scriptstyle{\text{{{{Hxc}}}}}}^{\scriptstyle{\text{{{{target}}}}}}\left[\rho\right]$ is smaller, and thus the IKS-DFPT1 is enough to reproduce the ground-state energy and density. However, as the PZ81 functional is non-polynomial, the energy density $\varepsilon_{\scriptstyle{\text{{{{c}}}}}}$ is not reproduced well in wide-density range within the present polynomial ansatz. We should go beyond the polynomial ansatz in the future.

4 Conclusion and Perspectives

In summary, the way to improve conventional EDFs based on the combination of the IKS and the DFPT is proposed. As benchmark calculations, we test whether the LDA exchange and correlation functionals can be reproduced in this novel scheme IKS-DFPT1. It is found that with the present polynomial ansatz the polynomial functional can be well reproduced, while the non-polynomial one can be reproduced in the crucial density region. By improving the exchange and correlation functionals, the accuracy of the ground-state energies is improved by two to three orders of magnitude, and the accuracy of the ground-state densities is also improved one to two orders of magnitude. Therefore, the IKS-DFPT is promising to improve the conventional functionals.

In this paper, we mainly focus on the feasibility of this method by using a simple polynomial ansatz shown in Eq. (5k@@italiccorr) and the noble-gas atoms. To get more accurate functional, beyond polynomial ansatz of $E_{\scriptstyle{\text{{{{Hxc}}}}}}^{(1)}\left[\rho\right]$, including considering the density gradient, is one of the interesting topics. Also, in this paper, two combinations of two systems, $\mathrm{He}$-$\mathrm{Ne}$ and $\mathrm{Xe}$-$\mathrm{Rn}$, are used as a benchmark, while the results for the other pairs are shown in the Appendix. It is a future task to get optimized values across all possible pairs or among the essentially different systems. As another perspective, the second-order IKS-DFPT is interesting. It is also important to include the spin and isospin degrees of freedom for applications of spin-polarized electron systems or nuclear systems.

According to the Hohenberg–Kohn theorem, the functional is, in principle, system independent. In practice, due to approximations, such as the first-order perturbation theory and the ansatz (5k@@italiccorr), $\alpha_{n}$ and $\lambda A_{n}$ have slightly system dependence as shown in the Appendix. If these approximations are appropriate, the system dependence is quite small. Therefore, any system can be used for this purpose as long as the density is known, and once densities of several systems are known the EDFs can be improved by using this novel method IKS-DFPT for actual problems.

The authors appreciate Ryosuke Akashi, Gianluca Colò, Yixin Guo, Xavier Roca-Maza, Osamu Sugino, Shinji Tsuneyuki, and Dario Vretenar for stimulating discussions and valuable comments. T.N. and D.O. acknowledge the financial support from Computational Science Alliance, The University of Tokyo. T.N. and H.L. would like to thank the RIKEN iTHEMS program and the JSPS-NSFC Bilateral Program for Joint Research Project on Nuclear mass and life for unravelling mysteries of the $r$-process. T.N. acknowledges the JSPS Grant-in-Aid for JSPS Fellows under Grant No. 19J20543. H.L. acknowledges the JSPS Grant-in-Aid for Early-Career Scientists under Grant No. 18K13549.

Appendix A Results for All Pairs of Noble Gases

The coefficients $\lambda A_{1}$ and $\alpha_{1}$ calculated in the first iteration for all the pairs of noble-gas atoms are shown in Table 5.