Improvement of Functionals in Density Functional Theory by the Inverse Kohn-Sham Method and Density Functional Perturbation Theory

Improvement of Functionals in Density Functional Theory by the Inverse Kohn-Sham Method and Density Functional Perturbation Theory

Daisuke Ohashi (\CJKfamilymin大橋大介) Department of Physics, Graduate School of Science, The University of Tokyo, Tokyo 113-0033, Japan RIKEN Nishina Center, Wako 351-0198, Japan    Tomoya Naito (\CJKfamilymin内藤智也) Department of Physics, Graduate School of Science, The University of Tokyo, Tokyo 113-0033, Japan RIKEN Nishina Center, Wako 351-0198, Japan    Haozhao Liang (\CJKfamilygbsn梁豪兆) haozhao.liang@riken.jp RIKEN Nishina Center, Wako 351-0198, Japan Department of Physics, Graduate School of Science, The University of Tokyo, Tokyo 113-0033, Japan
July 20, 2019
Abstract

We propose the way to improve energy density functionals in the density functional theory based on the combination of the inverse Kohn-Sham method and the density functional perturbation theory. 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.

preprint: RIKEN-QHP-403preprint: RIKEN-iTHEMS-Report-18{CJK*}

UTF8

The density functional theory (DFT) is one of the most successful approaches to calculate the ground-state properties of the quantum many-body problems, such as atoms, molecules, solids Hohenberg and Kohn (1964); Kohn and Sham (1965); Kohn (1999); Jones (2015), and nuclear systems Bender et al. (2003); Nakatsukasa et al. (2016). Since for the same level of accuracy, the numerical cost of DFT is much less than those of other quantum many-body methods Engel and Dreizler (2011), the DFT is applicable to larger-scale calculations Iwata et al. (2010); Soler et al. (2002); Ozaki (2006). Furthermore, in principle, the DFT gives the exact ground-state density and energy :

(1)

where is the Kohn-Sham (KS) kinetic energy, is the external field, and and are the Hartree and exchange-correlation energy density functionals (EDFs), respectively Hohenberg and Kohn (1964); Kohn and Sham (1965). However, in practice, the accuracy of the DFT calculation depends on the accuracy of the approximations for , as it is unknown.

In electron systems, many approximations for have been proposed from the first principle, i.e., non-empirically. The widely used ones are the local density approximation (LDA) Dirac (1930); Vosko et al. (1980); Perdew and Zunger (1981); Perdew and Wang (1992) and generalized gradient approximation (GGA) Becke (1988); Perdew et al. (1992, 1996a, 2008). In the LDA, is approximated as a functional of local density , whereas in the GGA, is approximated as a functional of and its gradient . Approximations beyond the GGA have also been discussed Perdew and Schmidt (2001). Even in the GGA, is constructed based on the LDA one. Thus, the exchange interaction is not fully included, and hence localized systems like -electron systems are sometimes unfavorable. In order to avoid these problems, phenomenological methods, such as the Anisimov et al. (1991); Liechtenstein et al. (1995); Petukhov et al. (2003), Anisimov et al. (1997), exact-exchange Görling (1996); Städele et al. (1997) methods, and hybrid functionals Becke (1993); Perdew et al. (1996b); Adamo and Barone (1998); Heyd et al. (2003) have been proposed. Functionals with the long-range correction due to the van der Waals interaction Iikura et al. (2001) and semi-empirical functionals Chai and Head-Gordon (2008); Peverati and Truhlar (2012) 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 Wiringa et al. (1995); Ishii et al. (2007); Aoki et al. (2012); Hatsuda and Kunihiro (1994); Holt et al. (2016); Fujita and Miyazawa (1957); Pieper et al. (2001). Even if the exact form of the interaction in the vacuum was known, the nuclear interaction in the medium is different due to its highly non-perturbative property Brueckner et al. (1954); Goldstone (1957). Therefore, it is still difficult to derive the Hartree-exchange-correlation EDF from the first principle, although the Hartree-Fock calculation from the interaction in the vacuum has been discussed for a long time Coester et al. (1970); Day (1967); Shen et al. (2016, 2017, 2018). Thus, for the nuclear interaction is treated phenomenologically Vautherin and Brink (1972); Dechargé and Gogny (1980); Long et al. (2006); Meng et al. (2006), with fitting parameters to the experimental data Bender et al. (2003). Since the fitting parameters are usually determined from the experimental data of several stable nuclei, different parameter sets can give totally different results for the exotic nuclei Kruppa et al. (2000). Comparisons between parameter sets are still being discussed Yoshida et al. (1998); Roca-Maza and Paar (2018).

Hence, the derivation of accurate EDFs is one of the most important topics for the fundamental DFT both in electron and nuclear systems. Recently, a new microscopic way to derive EDFs based on the functional renormalization group was suggested Liang et al. (2018); Yokota et al. (2018a, b); Yokota and Naito (2018), while it is not ready for realistic systems yet. Phenomenologically or empirically, even though the ways to derive EDFs have been discussed, there are still big debates on this topic, see e.g. Medvedev et al. (2017); Kepp (2017).

As an alternative way to improve EDFs, the inverse approach of DFT, the so-called inverse Kohn-Sham (IKS) method, was proposed Wang and Parr (1993); Zhao and Parr (1993). In the IKS, the KS potential is calculated from the given ground-state density , so that the information such as the single-particle energies is expected to be valuable for improving the accuracy of EDFs Almbladh and Pedroza (1984); Pedroza (1986); Umrigar and Gonze (1994). 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 Jensen and Wasserman (2018).

In order to attack this open question, a new strategy based on the density functional perturbation theory (DFPT), the so-called IKS-DFPT, will be shown in this paper. The DFPT was originally developed for the calculations of phonon and response properties from the first principle Baroni et al. (1987); Gonze (1995a); Gonze and Vigneron (1989); Baroni et al. (2001). This is because properties of solids are determined from phonons as well as electrons Ashcroft et al. (2016), but the original DFT is applicable only for the electron structure. In the DFPT, displacement of the external field caused by displacement of the nuclei or external electric field is treated by combination of the perturbation and linear response theories Baroni et al. (2001). In principle, each order of the DFPT can be handled Gonze (1995b), and so far the DFPT up to the third order has been discussed Gonze and Vigneron (1989).

In this Letter, for the first time, the way to improve EDFs from the IKS, the IKS-DFPT, is proposed by the combination with the DFPT. In the IKS-DFPT, the DFPT is used as the perturbation for the instead of , and the perturbed term is the deviation of the conventional EDF to the exact EDF . As benchmark calculations, we verify this method by reproducing both the LDA exchange functional Dirac (1930) and the LDA correlation functional Perdew and Zunger (1981). The iteration of IKS-DFPT is also discussed.

In the DFT, the ground-state density and energy of an -particle system are obtained by solving the self-consistent KS equation

(2)

where is the mass of particles, and are the single-particle orbitals and energies, respectively, and . Here, is the KS effective potential defined as

(3)

So as to improve a conventional Hartree-exchange-correlation functional, , from the IKS, is assumed to be close enough to the exact one, . We consider the effect of the difference between them in the first-order perturbation theory, as

(4)

with a small parameter . Then, the exact orbitals and ground-state density and energy are also expanded perturbatively:

(5a)
(5b)
(5c)

where quantities shown with the tilde are given by . The perturbation is assumed not to affect the external field, i.e., . Moreover, is assumed to be given, and thus are calculated from the IKS.

Under these assumptions, we calculate 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), (5a), and (5b) into Eq. (1) gives . In the latter way, Eq. (4) and integration of the KS equation (2) gives :

(6)

where are obtained from by using the IKS. By comparing these two expressions of the ground-state energy and neglecting term, the equation for is obtained:

(7)

Finally, the Hartree-exchange-correlation functional in the IKS-DFPT in the first-order, i.e., the IKS-DFPT1, is derived as

(8)

Because Eq. (8) is a functional equation, it is difficult to be solved directly. In this work, we assume

(9)

with the values of and to be determined, and then we get

(10)

To determine and , two systems, Systems 1 and 2, are required. Here, and are the exact ground-state densities, and and are the ground-state densities calculated from of the Systems 1 and 2, respectively. It leads to the two equations for and . In such a way, and can be determined.

Moreover, we also consider the iteration of the IKS-DFPT1. At the -th step of the iteration, the Hartree-exchange-correlation functional calculated in the IKS-DFPT1, in Eq. (8), is defined as

(11)

and in Eq. (4) is now , where is the original one at the first step. We perform the iteration until convergence.

As benchmark calculations, is calculated from the theoretical , instead of experimental data, and we test whether is reproduced in this scheme. All the pairs of the isolated noble-gas atoms (, , , , , and ) are used as Systems 1 and 2. In calculations, the Hartree atomic unit is used. The ADPACK code Ozaki et al. (2011) is used for the DFT calculations of the isolated atoms. Hereafter, is denoted by .

We calculate two cases: 1)  is the Hartree functional, and is the Hartree plus LDA exchange functional Dirac (1930):

(12a)
(12b)

2)  is the Hartree plus LDA exchange functional, and is the Hartree plus LDA exchange-correlation functional, where the PZ81 Perdew and Zunger (1981) functional is used:

(13a)
(13b)

In both cases, the external field is the Coulomb interaction between the nucleus and electron.

First, let us consider the first case, i.e., the calculations with Eqs. (12a) and (12b). In Tables 1 and 2, the coefficients and and the ground-state energies calculated in -th iteration are shown for the pairs of atoms - and -, respectively. It is found that and are obtained within and errors in -, and within and errors in -, respectively, from their exact values. The heavier atoms reproduce the coefficients better. The results of the other pairs are shown in the Supplemental Material.

The exchange energy density calculated in the first iteration, , and the ratio to the exact one, , are shown as a function of in Fig. 1 for the pairs of - and - with dashed and dot lines, respectively, while the exact one is shown with a solid line. Here, the energy density and the Wigner-Seitz radius are defined as (, ) and , respectively. The pair of - reproduces the exact functional within a few percents in the range of , which is generally better than the pair of -. Since the polynomial form of the functional in Eq. (9) 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 and the exact becomes smaller as the iteration proceeds further. The ground-state energies of , , , and are finally reproduced within , , , and errors, respectively, comparing with , , , and errors at the zeroth step, . This indicates the iteration helps the improvement of the ground-state energy.

The Wigner-Seitz radii calculated in the zeroth, first, and second iterations and the exact one for are shown as a functions of in Fig. 2 with dot-dashed, dashed, dot, and solid lines, respectively. The ratio of calculated Wigner-Seitz radius to the exact one, , for each step is also shown in the insert of Fig. 2. 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.

of of
Exact
Table 1: The coefficients and and the ground-state energies calculated in -th iteration for the pair of atoms and . The Hartree functional given in Eq. (12a) is used for and the Hartree plus LDA exchange functional given in Eq. (12b) are used for . All units are in the Hartree atomic unit.
of of
Exact
Table 2: Same as Table 1 but for the pair of atoms and .
Figure 1: Energy density for the LDA exchange functional calculated in the first iteration as a function of . Results for the pairs of - and - are shown with dashed and dot lines, respectively. The exact one is shown with a solid line. Ratios of are shown in the insert.
Figure 2: Wigner-Seitz radii as a function of for . Results calculated in the zeroth, first, and second iterations are shown with dot-dashed, dashed, and dot lines, respectively. The exact one is shown as a solid line. Ratios of are shown in the insert.

Next, let us consider the second case, i.e., the calculations with Eqs. (13a) and (13b). In Tables 3 and 4, the coefficients and and the ground-state energies calculated in -th iteration are shown for the pairs of atoms - and -, respectively. It is found that the ground-state energies are already reproduced well at the first step. Here, for the pair of -, the convergence reaches. For the pair of -, the ground-state energies are further improved slightly as the iteration proceeds further. The ground-state energies of , , , and are finally reproduced within , , , and error, respectively, comparing with , , , and errors at the zeroth step, . The results of the other pairs are shown in the Supplemental Material.

In order to compare the calculated correlation functionals and the exact one, the correlation energy density calculated in the first iteration is shown as a function of in Fig. 3 for the pairs of - and - with dashed and dot lines, respectively, while the exact one is shown with a solid line. The non-polynomial PZ81 functional is reproduced better in the lower-density region from the pair of -, and in the higher-density region from the pair of -, since heavier atoms have higher-density region than lighter atoms.

The Wigner-Seitz radii calculated in the zeroth and first iterations and the exact one for are shown as functions of in Fig. 4 with dot-dashed, dashed, and solid lines, respectively. The ratio of calculated Wigner-Seitz radius to the exact one, , for each step is also shown in the insert of Fig. 4. 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, the difference between and is larger, and thus 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 is reproduced well in wide-density range. In the second case, the difference between and 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 is not reproduced well in wide-density range within the present polynomial ansatz. We should go beyond the polynomial ansatz in the future.

of of
Exact PZ81 PZ81
Table 3: Same as Table 1 but the Hartree plus LDA exchange functional given in Eq. (13a) is used for and the Hartree plus LDA exchange-correlation functional given in Eq. (13b) is used for .
of of
Exact PZ81 PZ81
Table 4: Same as Table 3 but for the pair of atoms and .
Figure 3: Same as Fig. 1 but for .
Figure 4: Wigner-Seitz radii as a function of for . Results calculated in the zeroth and first iterations are shown with dot-dashed and dashed lines, respectively. The exact one is shown as a solid line. Ratios of are shown in the insert.

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. As perspectives, the second-order IKS-DFPT and beyond polynomial ansatz are interesting topics. It is also important to include the spin and isospin degrees of freedom for applications of spin-polarized electron systems or nuclear systems.

Acknowledgements.
The authors appreciate Ryosuke Akashi, Gianluca Colò, Xavier Roca-Maza, Osamu Sugino, Shinji Tsuneyuki, and Dario Vretenar for stimulating discussions and valuable comments. DO and TN acknowledge the financial support from Computational Science Alliance, The University of Tokyo. TN and HL 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. HL acknowledges the JSPS Grant-in-Aid for Early-Career Scientists under Grant No. 18K13549.

References

Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
340699
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description