Energetics of the AK13 Semi-Local Kohn-Sham Exchange Energy Functional
Abstract
The recent non-empirical semi-local exchange functional of Armiento and Kümmel, the AK13 [PRL 111, 036402 (2013)] incorporates a number of features reproduced by higher-order theory. The AK13 potential behaves analogously with the discontinuous jump associated with the derivative discontinuity at integer particle numbers. Recent works have established that AK13 gives a qualitatively improved orbital description compared to other semi-local methods, and reproduces a band structure closer to higher-order theory. However, its energies and energetics are inaccurate. The present work further investigates the deficiency in energetics. In addition to AK13 results, we find that applying the local-density approximation (LDA) non-self-consistently on the converged AK13 density gives very reasonable energetics with equilibrium lattice constants and bulk moduli well described across 14 systems. We also confirm that the attractive orbital features of AK13 are retained even after full structural relaxation. Hence, the deficient energetics cannot be a result of the AK13 orbitals having adversely affected the quality of the electron density compared to that of usual semi-local functionals; an improved orbital description and good energetics are not in opposition. We also prove that the non-self-consistent scheme is equivalent to using a single external-potential dependent functional in an otherwise consistent KS-DFT scheme. Furthermore, our results also demonstrate that, while an internally consistent KS functional is presently missing, non-self-consistent LDA on AK13 orbitals works as a practical non-empirical computational scheme to predict geometries, bulk moduli, while retaining the band structure features of AK13 at the computational cost of semi-local DFT.
I Introduction
Density-functional theory (DFT) [1], in its standard formulation of Kohn and Sham (KS) [2], is a very successful approach to the many-electron problem. It is generally accurate, and requires relatively little computational effort, compared to, e.g., wave-function based methods. While it is exact in principle, its key ingredient, the exchange-correlation (xc) energy functional is in practical calculations approximated. Many viable approximations have been proposed.
A very common class of functionals rely solely on semi-local information in the density, i.e., to some order , where is the electronic density. Semi-local DFT gives good results for structural properties. It is well established that there are a number of highly relevant exchange features that semi-local DFT omits, such as the derivative discontinuity (DD) [3] and the relative offset between potentials for well-separated subsystems [4]. These are closely related to the self-interaction error and lead, e.g., to over-delocalized KS orbitals, and which arguably constitutes a major deficiency in semi-local approximations. The missing features are to a large degree restored in higher-order methods such as hybrids [5], exact exchange (EXX) with the optimized effective potential (OEP) [6; 7], and many-body perturbation theory (in particular, the GW approximation [8]). These methods, however, increase the computational cost. This is in many situations undesirable, and hence it is a worthwhile goal of broad interest for applications across a large number of scientific fields to pursue a way to incorporate the exchange features into semi-local DFT.
Recently, one of us proposed a semi-local functional [9] (AK13) which achieves an orbital description more similar to that of higher order methods. It has been shown to improve the electronic structure, optical dielectric constants [10], and potential shape [12; 11], compared to other semi-local functionals. However, its energetics and total exchange energies are generally worse [9; 12]. In particular, bulk Al gives an equilibrium lattice constant larger than the exchange-only version of the functional by Perdew, Burke, and Ernzerhof [13] (PBE). In the present work we further investigate this deficiency in energetics. We find that the average error is about compared to the couple of per cent error for typical semi-local functionals.
The AK13 functional differs significantly from usual functionals on generalized gradient form (GGA) in that its refinement function is strongly divergent with the reduced density gradient, which heavily influences the resulting orbitals. With this background, one may ask if the deficiency in energetics stems from a misfeature directly in the AK13 electron density caused by qualitative difference in the orbitals, or, if the deficiency instead primarily stems from the energies assigned by the AK13 functional to those densities. This is a relevant question, because the former option means the deficiency would be inherently associated with the main beneficial property of the AK13 functional, making it very challenging to address. A main point of the present work is to show that this is not the case.
It is not clear how one can deduce the quality of the density for energetics from looking at the density itself. However, it has been a standard practice in the field to test new functionals by investigating their action on a known good electron density, e.g., a density from the local density approximation (LDA). Here, we reverse that test, and compare the results of inserting self-consistent densities from different functionals into the LDA functional. This method shows if reasonable energetics are retained for the AK13 density. As one may expect, we see that AK13 densities in the LDA functional does not visibly alter energies on the scale of typical total energies. However, the same turns out not to be true for the energetics (though reasonable energetics are acheived).
Hence, the primary contributions in this work to the field of semi-local methods with improved exchange features are: (i) we investigate more closely the deficiency in energetics of the AK13 functional, (ii) we answer in the affirmative that the qualitative changes to the AK13 orbitals have not adversely affected its density in a way that prevents reproducing reasonable energetics; (iii) we suggest a direction for how to proceed to alter the energy functional of AK13 without significantly changing its orbitals or density; and (iv) our results suggest that in lieu of a single consistent KS-DFT functional, post-corrected AK13 works in practice as a computational scheme that is predictive for crystal geometry while at the same time retains the band structure features of AK13 with, e.g., increased KS band gaps more similar to those of higher-order theories. We are only aware of prior works demonstrating the success of AK13 band structure for systems using experimental geometries.
The rest of the paper is organized as follows: Section II summarizes the derivation and properties of the AK13 exchange functional. In Sec. III we show how the scheme of using LDA on converged AK13 orbitals is equivalent to a single consistent external potential-dependent xc functional. In Sec. IV, we show the numerical results for atoms and solids. In Sec. IV we discuss our results and their implications for a future KS DFT functional. We summarize and conclude our work in Sec. VI.
Ii The AK13 exchange functional
We begin by briefly revisiting the derivation and motivation of the AK13 functional [9]. In the Perdew, Parr, Levy and Balduz ensamble extension of DFT [3], the total energy of a system with an electron density of electrons, with and moving in some external potential is
(1) |
where is the non-interacting kinetic energy, the classical Coulomb-interaction (Hartree) energy, and the exchange-correlation energy. The last term includes all remaining quantum many-body effects that we seek to approximate. A minimization of Eq. (1) results in the KS equations, given by (Hartree atomic units are used throughout)
(2) |
where are the KS orbitals, each corresponding to the eigenvalue , and the external potential from the atomic cores, the Hartree potential, and the xc potential, which is defined as , and which usually is separated into a sum of its parts, the exchange (x) and correlation (c).
The starting point of the derivation of the AK13 functional is the observation that the model potential of Becke and Johnson [14] (BJ), , reproduces a behaviour that stems from the derivative discontinuity, specifically it undergoes a constant but sudden shift when the particle number changes across integer particle numbers [15]. This is due to that, (i) the expression given by BJ takes an asymptotic value far outside a finite system that depends on the eigenvalue of highest occupied KS orbital , where (HOMO), and, (ii) BJ have stated that the potential may be manually shifted with the system-dependent constant to achieve an asymptotic value of zero. Various modifications of the BJ potential have been shown to give band structures closer to higher-order theory [16; 17], polarizabilities [15], as well as atomic and molecular properties [18; 19]. In contrast, the main feature of the AK13 functional is a potential that reproduces the same asymptotic behavior, but from a consistent energy-potential pair that avoids any theoretical or practical issue with model potentials [20].
A semi-local exchange-functional on generalized-gradient approximation (GGA) [21] form is defined as
(3) |
so that , and is a function of the reduced density gradient
(4) |
where is the Fermi wave vector. The choice corresponds to LDA exchange. In Ref. [9], the sought asymptotic behavior is shown to be fulfilled by
(5) |
where , and . This can be verified by inserting the asymptotic density outside a finite system [22], i.e.,
(6) |
where is a system-dependent constant, and is the occupation number of the HOMO orbital. The corresponding exchange potential
(7) |
which thus shifts discontinuously as a new orbital is occupied. Similar to the BJ model potential, one could propose a shift of the potential to zero alignment by adding a constant , i.e.,
(8) |
which is defined as
(9) |
It follows that the expected [3] asymptotic limit of , as is retrieved. An exact expression for in a finite system can be derived [9], and is given by
(10) |
where the plus sign applies for eigenvalues and . However, as it was argued in Ref. [9] there is no formal requirement in KS DFT to perform this shift. Without the shift, the energy functional and potential pair remains consistent. Leaving out the shift may seem sufficient, since it has no effect on the density, and only appears as a constant shift of equal magnitude of the eigenvalues. Nevertheless, if a functional could be constructed to give a shifted potential, that would most likely affect the energies. Similarly, in this work, we explore the consequences of a scheme that alters the energy without changing the orbitals. This idea is somewhat formalized in terms of a single density and potential-dependent functional in the next section.
Iii Non-self-consistent LDA energies from self-consistent AK13 orbitals
As we outlined in Secs. I and II, the aim of the present work is to investigate the energetics that result from a non-self-consistent evaluation of the LDA energy on orbitals obtained self-consistently with the AK13 functional. In practice, the implementation of this scheme is completely straightforward as an extra energy evaluation after electronic convergence in a regular AK13 calculation. In this section we discuss how this scheme can alternatively be seen as a single self-consistent functional, but which requires an external potential dependence.
Our starting point is a general, non-unique partition of the xc-functional
(11) |
consisting of semi-local (sl) and non-local (nl) parts. This relation defines the functional . It is allowed to be non-local in the sense that it may depend on the density throughout the entire system. Furthermore, we let , and, from now on, take , where is a suitable choice of a semi-local correlation-energy functional.
For a given external potential , define the density that results from a self-consistent solution of the KS DFT problem of only the sl part of the problem as , which via Eq. (4) also gives . It is now possible to make a definition of the contribution from to give exactly LDA energetics as
(12) |
and [23]. The definition in Eq. (12) is chosen to make all but the first term in Eq. (5) cancel out in the energy functional, but will by constraction not affect the orbitals. This definition is not a valid KS DFT functional due to its explicit dependence on the external potential from . However, such a dependence in the exchange-correlation functional is frequently seen in other schemes, and the present scheme is thus on equal fundamental footing with them. More importantly, Eq. (11) gives an explicit construction idea to explore for a future functional within KS DFT with the sought properties.
Iv Results on Atomic energies and the energetics of solids
Our first numerical results shows the energy change by using non-self-consitent LDA on converged AK13 orbitals for total atomic energies. We use a modified atomic DFT code, originating from an early work of Talman and Shadwick [7], to perform fully self-consistent spin-polarized DFT calculations on a one-dimensional grid of mesh points, . Figure 1 shows the total energy of a magnesium ion as a function of the electron occupation number, i.e, the system is successively filled with electrons, keeping the external potential of the atomic core fixed. In Figure 2 we show the total exchange energies for the same system. Observe that the results depicted in Fig. 2 have been obtained from exchange only calculations. In both cases, we find that self-consistent AK13 energies are consistently lower than those of both self-consistent LDA and exact exchange OEP (EXX-OEP). Furthermore, the energies calculated from LDA applied to the converged AK13 orbitals stay very close to those of the self-consistent LDA, i.e., for the LDA energy functional there is hardly any difference between using the LDA or AK13 orbitals.
Next we consider periodic solids. To this end we use Elk, a full-potential, linearized augmented plane wave (LAPW) [24] code. For each solid, we calculate the total energy at 7 different volumes in an approximate interval of around the equilibrium volume . The obtained values we fit to the Birch-Murnaghan equation of state (EOS) [25], from which we obtain the optimal lattice geometry and bulk modulus. Specifically, the lattice parameters of this work were calculated by choosing volume points in the vicinity of the equilibrium volume as obtained from self-consistent LDA. Since the optimal volume points for AK13 would be at a range of larger volumes than the points optimal for LDA, our calculated AK13 equilibrium lattice parameters and bulk moduli are not as accurate as the values reported for the other methods. However, the accuracy is still sufficient to illustrate the poor performance of AK13 for these quantities, and we thus find it unnecessary to optimize the volume points for greater accuracy. In all calculations we take to be the parameterized LDA functional of Perdew and Wang (PW) [23]. In all calculations we use -points, 50 empty bands, and a cutoff of , where is the muffin-tin radius. In Elk, the values of are set as default, which means that they are automatically rescaled when overlapping (however, the results were carefully checked so as to not be sensitive to minor changes in ).
Results for silicon in the diamond structure are shown in Figs. 3, and 4. Figure 3 shows how AK13 gives a much too large lattice minimum. When we use the orbitals of AK13 to calculate the non-self consistent LDA contribution, however, the result is a minimum much closer to self-consistent LDA. The shape of the energy vs. lattice constant curve in this case closely resembles that of self-consistent LDA. Furthermore, Fig. 4 shows two different sets of KS bands, each of which has been calculated at the equilibrium lattice constant when using non-self-consistent LDAxc on AK13 orbitals. This, together with the result shown in Fig. 3, illustrates that the main features of the AK13 orbitals are unaffected by structural relaxation.
To make sure Si is not a single fortuitous case, we repeat this analysis for 14 different solids and tabulate the equilibrium lattice constants and bulk moduli . The results are shown in Tables 1 and 2. For the lattice constants , we find that non-self consistent LDA on AK13 orbitals gives lattice constants that are slightly smaller, with a mean error (me) of about -0.2 Å, than those obtained from self-consistent LDA. The bulk moduli have also been substantially improved compared to self-consistent AK13. The calculated values of are reasonably close to those of self-consistent LDA. In Table 2, we have also included the result from the Perdew, Burke, and Ernzerhof (PBE) functional for comparison.
At the optimized geometries, we have calculated the KS gaps, . These are shown in Table 3. We emphasize that one cannot make a direct quantitative comparison between KS band gaps and experimental gaps. However, as has been argued previously [10], the significance of these results is that the present method retains the AK13 band structure that is more similar to that of higher-order theory.
V Discussion
The results in Table 1 demonstrate that the use of non-self-consistent LDA energies on converged AK13 orbitals results in reasonable values of equilibrium lattice constants (mean average error (mae) of 0.1 Å, mean average relative error (mare) of 0.02), although with a tendency to over-bind even somewhat more than self-consistent LDA. The over-binding seems to be a general trend among the studied systems, although somewhat surprisingly the reverse is true for Al. For the bulk moduli in Table 2, our approach gives an mae of 9.7 GPa (mare 0.3), which is a small difference, i.e., of the same magnitude that one can expect from using various implementations found in different computational codes.
The orbital features (i.e., features in the KS band structure and band gaps) appear to be more or less unaffected by the change in lattice constant due to structural relaxation using the LDA energy on converged AK13 orbitals as compared to AK13 on experimental geometries. While the quality of the DFT orbital description is a difficult topic, we have discussed in previous works [10] how we find band gaps and orbital features calculated with various methods to end up in one of two categories. Either (i) the orbital features have the typical LDA and PBE characteristics, e.g., with band gaps much smaller than the experimental ones and with a very reduced step structure in the KS potential for atoms; or (ii) the orbitals share features of higher-order methods (in particular, exact exchange DFT), i.e., they are more similar to actual quasiparticle states and, as a result, give expanded band gaps and optical properties more similar to experiments, etc. The orbitals of AK13 are in the latter category, both when calculated on experimental geometries and, as shown in this work, on relaxed geometries (when using LDA for the energies.)
The stated main focus of the present work is to establish if the inaccurate energetics of the AK13 functional is an inherent feature of the changed orbitals, as compared to other, general semi-local functional approximations. In light of our findings above, this appears to not be the case, since the LDA energy functional gives a physically acceptable energy minimum when applied to those orbitals. Furthermore, we have shown that the scheme of using AK13 to converge the orbitals and then LDA to calculate the energies is, in fact, equivalent to using a single consistent exchange-correlation functional, but, with an explicit external potential dependence, Eq. (12). From these results, we surmise that it is likely possible to find a density functional within KS DFT that is capable of incorporating both these features at the computational cost of semi-local DFT. We suggest that the external-potential dependent expression in Eq. (12) is a possible starting point for such a construct.
It is pertinent at this point to ask why the LDA energetics are slightly worsened when using orbitals we expect to be more accurate than the LDA orbitals (see Fig. 3 and Table 1). In the following we will speculate on a few possible reasons: (i) the use of non-self-consistent LDA on converged AK13 orbitals is meant as an approximation at what a fully self-consistent KS DFT scheme can achieve, and thus the lack of a fully self-consistent scheme may itself lead to worsened energetics; (ii) The performance of LDAxc relies on a cancellation of error between exchange and correlation [29; 30]. We have used LDA correlation when converging the AK13 orbitals, which is not ideal. A correlation functional better matched with AK13 exchange may influence the orbitals enough to recover LDA energetics or better; (iii) The AK13 functional may yield an improved orbital description in the sense of reproducing exchange features missing from other semi-local functionals, but it was never optimized for accurate energies or energetics. It may be that the AK13 form with its non-empirically derived values for and in Eq. (5) lead to other minor deficiencies in the orbitals that, in contrast to the main conclusion of this paper, inherently worsen energetics somewhat compared to LDA. However, if true, the conclusion would be that this effect is on a much smaller scale than what one is led to believe from the energetics of self-consistent AK13.
LDA | PBE | AK13 | This work | Expt.^{a} | |
Al | |||||
Si | |||||
-Sn | |||||
Ge | |||||
GaAs | |||||
InAs | |||||
Ca | |||||
Li | |||||
K | |||||
LiCl | |||||
AlN | |||||
TiN | |||||
AlAs | |||||
Sr | |||||
me | |||||
mae | |||||
mare |
LDA | PBE | AK13 | This work | Expt.^{a} | |
Al | |||||
Si | |||||
-Sn | |||||
Ge | |||||
GaAs | |||||
InAs | |||||
Ca | |||||
Li | |||||
K | |||||
LiCl | |||||
AlN | |||||
TiN | |||||
AlAs | |||||
Sr | |||||
me | |||||
mae | |||||
mare |
Vi Summary and Conclusions
In this work we have studied the deficiency of the energetics of the AK13 functional and, in doing so, taken a step towards a semi-local KS density functional that combines accurate energetics with the improved orbital description similar to that of higher-order theories. We have demonstrated that the description of exchange bonding in AK13 is generally highly unsatisfactory. We have shown that this deficiency cannot primarily be attributed to how the qualitatively different AK13 orbitals changes the electron density. Rather, one finds LDA-like energetics to be recovered when the self-consistent AK13 density is inserted into the LDA xc functional. We have also demonstrated that the changes in the AK13 band structure are robust under changes in geometry, and thus are retained when the structure undergoes structural relaxation based on the LDA energy.
The conclusion that the AK13 energetics cannot primarily be attributed to misfeatures in its self-consistent density opens for a hypothetical future functional that mostly retains the AK13 orbitals while improving the energetics. We take a first step towards such a functional by proving that the use of LDA xc on AK13 densities is equivalent to a single consistent xc functional that has an explicit dependence on the external potential. Future work is necessary to turn this approach into a functional fully within KS DFT. However, even without such a functional, our results suggest that LDA on self-consistent AK13 densities works in practice as a computational scheme that successfully combines improved band structure features with predictive lattice constants and bulk moduli. However, the tendency to over-bind is somewhat stronger even than ordinary self-consistent LDA.
Acknowledgements.
We thank Stephan Kümmel and Thilo Aschebrock for insightful discussions, and valuable input on an early version of the manuscript. R. A. gratefully acknowledges support from the Swedish Research Council (VR), Grant No. 621-2011-4249 as well as the Linnaeus Environment at Linköping on Nanoscale Functional Materials (LiLi-NFM) funded by VR.References
- Hohenberg and Kohn [1964] P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
- Kohn and Sham [1965a] W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965a).
- Perdew et al. [1982] J. P. Perdew, R. G. Parr, M. Levy, and J. L. Balduz, Phys. Rev. Lett. 49, 1691 (1982).
- Gritsenko and Baerends [1996] O. V. Gritsenko and E. J. Baerends, Phys. Rev. A 54, 1957 (1996).
- Becke [1993] A. D. Becke, J. Chem. Phys. 98, 5648 (1993).
- Sharp and Horton [1953] R. T. Sharp and G. K. Horton, Phys. Rev. 90, 317 (1953).
- Talman and Shadwick [1976] J. D. Talman and W. F. Shadwick, Phys. Rev. A 36, 36 (1976).
- Hedin [1965] L. Hedin, Phys. Rev. 139, A796 (1965).
- Armiento and Kümmel [2013] R. Armiento and S. Kümmel, Phys. Rev. Lett. 111, 036402 (2013).
- Vlcek et al. [2015] V. Vlcek, G. Steinle-Neumann, L. Leppert, R. Armiento, and S. Kümmel, Phys. Rev. B 91, 035107 (2015).
- Tran et al. [2015] F. Tran, P. Blaha, M. Betzinger, and S. Blügel, Phys. Rev. B 91, 165121 (2015).
- Cerqueira et al. [2014] T. F. T. Cerqueira, M. J. T. Oliveira, and M. A. L. Marques, J. Chem. Theory Comput. 10, 5625 (2014).
- Perdew et al. [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
- Becke and Johnson [2006] A. D. Becke and E. R. Johnson, J. Chem. Phys. 124, 221101 (2006).
- Armiento et al. [2008] R. Armiento, S. Kümmel, and T. Körzdörfer, Phys. Rev. B 77, 165106 (2008).
- Tran et al. [2007] F. Tran, P. Blaha, and K. Schwarz, J. Phys. Condens. Matter 19, 196208 (2007).
- Tran and Blaha [2009] F. Tran and P. Blaha, Phys. Rev. Lett. 102, 226401 (2009).
- Pittalis et al. [2010] S. Pittalis, E. Räsänen, and C. R. Proetto, Phys. Rev. B 81, 115108 (2010).
- Oliveira et al. [2010] M. J. T. Oliveira, E. Räsänen, S. Pittalis, and M. A. L. Marques, J. Chem. Theory Comput. 6, 3664 (2010).
- Karolewski et al. [2013] A. Karolewski, R. Armiento, and S. Kümmel, Phys. Rev. A 88, 052519 (2013).
- Perdew and Wang [1989] J. P. Perdew and Y. Wang, Phys. Rev. B 33, 8800 (1989).
- van Leeuwen and Baerends [1994] R. van Leeuwen and E. J. Baerends, Phys. Rev. A 49, 2421 (1994).
- Perdew and Wang [1992] J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992).
- Singh [1994] D. J. Singh, Planewaves, Pseudopotentials and the LAPW Method (Kluwer Academic Publishers, Boston, 1994).
- Fu and Ho [1983] C. L. Fu and K. M. Ho, Phys. Rev. B 28, 5480 (1983).
- Ogilvie and Wang [1992] J. F. Ogilvie and F. Y. H. Wang, J. Mol. Struct. 273, 277 (1992).
- Ogilvie and Wang [1993] J. F. Ogilvie and F. Y. H Wang, J. Mol. Struct. 291, 313 (1993).
- Vurgaftman et al. [2001] I. Vurgaftman, J. R. Meyer, and L. R. Ram-Mohan, J. Appl. Phys 89, 5815 (2001).
- Kurth et al. [1999] S. Kurth, J. P. Perdew, and P. Blaha, Int. J. Quan. Phys. 75, 889 (1999).
- Fiolhais et al. [2003] C. Fiolhais, F. Nogueira, and M. Marques (eds.), A Primer in Density Functional Theory(Springer, Berlin, 2003) .