A self-interaction-free local hybrid functional: accurate binding energies vis-à-vis accurate ionization potentials from Kohn-Sham eigenvalues

A self-interaction-free local hybrid functional: accurate binding energies
vis-à-vis accurate ionization potentials from Kohn-Sham eigenvalues

Tobias Schmidt Theoretical Physics IV, University of Bayreuth, 95440 Bayreuth, Germany    Eli Kraisler Department of Materials and Interfaces, Weizmann Institute of Science, Rehovoth 76100,Israel    Adi Makmal [ Department of Materials and Interfaces, Weizmann Institute of Science, Rehovoth 76100,Israel    Leeor Kronik Department of Materials and Interfaces, Weizmann Institute of Science, Rehovoth 76100,Israel    Stephan Kümmel Theoretical Physics IV, University of Bayreuth, 95440 Bayreuth, Germany
July 26, 2019
Abstract

We present and test a new approximation for the exchange-correlation (xc) energy of Kohn-Sham density functional theory. It combines exact exchange with a compatible non-local correlation functional. The functional is by construction free of one-electron self-interaction, respects constraints derived from uniform coordinate scaling, and has the correct asymptotic behavior of the xc energy density. It contains one parameter that is not determined ab initio. We investigate whether it is possible to construct a functional that yields accurate binding energies and affords other advantages, specifically Kohn-Sham eigenvalues that reliably reflect ionization potentials. Tests for a set of atoms and small molecules show that within our local-hybrid form accurate binding energies can be achieved by proper optimization of the free parameter in our functional, along with an improvement in dissociation energy curves and in Kohn-Sham eigenvalues. However, the correspondence of the latter to experimental ionization potentials is not yet satisfactory, and if we choose to optimize their prediction, a rather different value of the functional’s parameter is obtained. We put this finding in a larger context by discussing similar observations for other functionals and possible directions for further functional development that our findings suggest.

pacs:
31.15.ep, 31.15.eg, 31.10.+z, 71.15.Mb
thanks: These authors contributed equallythanks: These authors contributed equally

Present address: ]Institut für Theoretische Physik, Universität Innsbruck, Technikerstraße 25, A-6020 Innsbruck, Austria

I Introduction

Kohn-Sham (KS) density-functional theory (DFT)Hohenberg and Kohn (1964); Kohn and Sham (1965) has become one of the most frequently used theories for electronic structure calculations. It employs the electron ground-state density, , as the central quantity and accounts for all electronic interaction beyond the classical electrostatic (Hartree) repulsion, , via the exchange-correlation (xc) energy functional, Parr and Yang (1989); Dreizler and Gross (1995); Fiolhais et al. (2003). Even though the xc energy is typically the smallest component in the ground-state total energy, it governs binding properties, geometrical structures, and ionization processes Kurth and Perdew (2000); Burke (2012); Fiolhais et al. (2003). Thus, the quality of a DFT calculation depends decisively on the functional approximation put to task.

It has become popular to categorize density functional approximations according to the “Jacob’s ladder” scheme introduced in Ref. Perdew and Schmidt (2001). Typically, the accuracy of a density functional approximation (DFA) improves when more “ingredients” are allowed in the functional construction, at the price of increased complexity. The local spin-density approximation (LDA) Kohn and Sham (1965), which approximates based on the xc energy of the homogeneous electron gas Ceperley and Alder (1980); Perdew and Wang (1992); Perdew and Zunger (1981); Vosko et al. (1980), and even more so the semi-local generalized gradient approximations (GGAs) Perdew and Wang (1986); Becke (1988); Lee et al. (1988); Perdew (1991); Perdew et al. (1996a); Wu and Cohen (2006); Haas et al. (2011), which additionally take the density gradient into account, offer a favorable ratio of computational expense and accuracy Burke (2012). Hybrid functionals Becke (1993a, b); Perdew et al. (1996b); Stephens et al. (1994); Adamo and Barone (1999); Ernzerhof and Scuseria (1999) typically reach yet greater accuracy by combining a fixed percentage of Fock exchange

(1)

with (semi-)local exchange and correlation energy terms (Hartree atomic units are used throughout). Eq. (1) evaluated with the exact Kohn-Sham orbitals defines the exact Kohn-Sham exchange energy. Self-consistent Kohn-Sham calculations based on the energy of Eq. (1) use the optimized effective potential (OEP) equation (see Grabo et al. (1997); Engel and Dreizler (2011); Kümmel and Kronik (2008) and references therein). When we use the abbreviation EXX in the following, we always refer to this Kohn-Sham variant of exact exchange.

While the aforementioned functionals in many cases predict binding energies and bond-lengths reliably, semi-local DFAs and to some extent also hybrid functionals are less reliable for ionization processes, photoemission spectra and densities of states. Very early on it was realized that this problem is closely related to the (one-electron) self-interaction (SI) Perdew and Zunger (1981) error, i.e., to the fact that in exact DFT should vanish for any one-electron system, but does not do so for these DFAs. Due to the SI error and the fact that semi-local functionals “average over” the derivative discontinuity Perdew and Levy (1983); Tozer and Handy (1998), the Kohn-Sham eigenvalues of the above mentioned approximate functionals typically fulfill neither the exact condition that the highest occupied eigenvalue should match the first ionization potential (IP) Perdew et al. (1982); Levy et al. (1984); Almbladh and von Barth (1985); Perdew and Levy (1997), nor the approximate but for practical purposes equally important condition that upper valence eigenvalues are good approximations to higher IPs when they are calculated from accurate xc potentials Stowasser and Hoffmann (1999); Chong et al. (2002); Körzdörfer et al. (2009, 2010); Dauth et al. (2011); Bleiziffer et al. (2013).

Although the interpretation of occupied eigenvalues even with the exact xc potential is approximate (except for ), it is of great practical importance. For example, the band-structure interpretation of Kohn-Sham eigenvalues has had a great impact on solid-state physics and materials science Cohen and Chelikowsky (1988). In recent years the interpretation of eigenvalues has become particularly important in the field of molecular semiconductors and organic electronics. Efforts to understand, e.g., photoemission experiments, have revealed severe shortcomings of traditional DFAs that go considerably beyond a spurious global shift of the eigenvalue spectrum Dori et al. (2006); Mundt et al. (2006); Mundt and Kümmel (2007); Körzdörfer et al. (2009, 2010); Körzdörfer and Kümmel (2010); Dauth et al. (2011); Marom et al. (2008); Marom and Kronik (2009); Marom et al. (2010); Bisti et al. (2011). A similar problem is witnessed also in solid state systems Rinke et al. (2005); Fuchs et al. (2007); Fuchs and Bechstedt (2008); Rödl et al. (2009); Betzinger et al. (2012, 2013). We emphasize that these problems of interpretation arise already for the occupied eigenvalues, i.e., the issues are separate from the well known band-gap problem Perdew and Levy (1983); Sham and Schlüter (1983); Kümmel and Kronik (2008); Kronik et al. (2012) of Kohn-Sham theory. The KS EXX potential leads to band structures and eigenvalues that match experiments much better than the eigenvalues from (semi-)local approximations Bylander and Kleinman (1995, 1996); Grabo et al. (1997); Städele et al. (1999); Mundt et al. (2006); Engel and Schmid (2009).

A comparison to hybrid functionals is more involved, because already the occupied eigenvalue spectrum depends sensitively on whether one uses the hybrid functional in a KS or a generalized KS calculation Körzdörfer and Kümmel (2010). For well understood reasons Perdew and Levy (1983); Sham and Schlüter (1983), the differences between the KS and the generalized KS eigenvalues become yet larger for unoccupied eigenvalues (see, e.g., the review Kronik et al. (2012)).This article’s focus is on Kohn-Sham theory, therefore we do not discuss the comparison to hybrid functionals used in the generalized KS scheme in detail. We note, however, that in particular range-separated hybrid functionals used in the generalized KS approach can predict gaps and band structures quite accurately, as discussed, e.g., in Henderson et al. (2011); Kronik et al. (2012), but global hybrid functionals tend to yield a less reliable density of states for complex systems than self-interaction free Kohn-Sham potentials Körzdörfer and Kümmel (2010); Dauth et al. (2011).

Besides these practical benefits, EXX also appears as a natural component of DFAs because it may be considered attractive to treat as many energy components as possible exactly, and including EXX has shown to be beneficial for, e.g., describing ionization, dissociation and charge transfer processes  Kümmel and Kronik (2008). However, bare EXX is a very poor approximation for binding energies, (see, e.g., Refs. Engel et al. (2000),Engel and Dreizler (2011), and Fiolhais et al. (2003), chapter 2). Combining EXX with a (semi-)local correlation term in many situations leads to results of inferior quality compared to pure EXX or semi-local DFAs, because of an imbalance between the delocalized exchange hole and the localized correlation hole Gunnarsson and Lundqvist (1976); Perdew and Schmidt (2001); Prodan and Kohn (2005); Kümmel and Kronik (2008).

One promising approach for combining EXX with appropriate correlation in a balanced way is the local hybrid form Cruz et al. (1998); Jaramillo et al. (2003); Perdew et al. (2008)

(2)

Here, is the xc energy density per particle that yields the xc energy via . The quantities and denote exchange and correlation energy densities per particle, respectively, approximated with (semi-)local expressions, whereas represents the EXX energy density per particle deduced from Eq. (1). The function is the local mixing function (LMF). It is a functional of the density and a decisive part of the local hybrid concept.

Eq. (2) can be viewed as a generalization of the common (global) hybrids. Instead of a fixed amount of EXX, the local hybrid can describe different spatial regions of a system with varying combinations of EXX and semi-local xc, by means of (where ). For example, whereas one-electron regions are supposed to be well-described using EXX, regions of slowly varying density are expected to be captured appropriately by (semi-)local xc functionals. The idea of local hybrids can also be understood in terms of the adiabatic connection theorem Ernzerhof et al. (1997), because may offer further flexibility in an accurate construction of the coupling-constant-dependent xc energy Cruz et al. (1998), especially for small coupling constant values.

The local hybrid form was pioneered by Jaramillo et al. Jaramillo et al. (2003) with a focus on reducing the one-electron SI-error in single-orbital regions. Numerous further local hybrid constructions followed Arbuznikov et al. (2006); Arbuznikov and Kaupp (2007); Bahmann et al. (2007); Kaupp et al. (2007); Perdew et al. (2008); Arbuznikov et al. (2009); Haunschild et al. (2009); Haunschild and Scuseria (2010a, b); Theilacker et al. (2011). They proposed various LMFs with different one-electron-region indicators, suggested several (semi-)local exchange and correlation functionals to be used in the construction, and followed different procedures to satisfy known constraints and determine remaining free parameters.

In the present manuscript we propose a new local-hybrid approximation that combines full exact exchange with a compatible correlation functional. The development is guided by the philosophy of fulfilling known constraints Perdew et al. (2005): Our xc energy density per particle, , is one-electron SI-free, possesses the correct behaviour under uniform coordinate scaling, and has the right asymptotic behaviour at large distances. It includes one free parameter that is not determined uniquely from these constraints.

In difference to earlier work, our emphasis is not on improving further the accuracy of binding energies beyond the one that was achieved with global hybrids. Instead, we focus on whether it is possible to construct an approximation that yields binding energies of at least the same quality as established hybrids and at the same time affords other advantages, notably KS eigenvalues that approximate IPs reasonably well. We find that if we choose the parameter in our functional by optimizing the prediction of binding energies, the latter are obtained with an accuracy that is similar to the one reached with usual global hybrids. At the same time, we achieve a significant improvement in prediction of dissociation energy curves for selected systems. Improvement in prediction of the ionization energy via the highest occupied KS eigenvalue is also observed. It is especially large for alkali atoms. However, the quality of the ionization energy prediction is not yet satisfactory, and if we aim to optimize the prediction of the latter, a rather different value for the functional’s free parameter is obtained. We put this finding in a larger context by discussing similar observations for other functionals.

The paper is organized as follows: Section II is devoted to the description of the new local hybrid functional. Section III (and the Appendix) provide methodological and computational details. Sec. IV presents and discusses the results, and Sec. V offers conclusions and a summary.

Ii Construction of the functional

In the construction of our functional, we choose to concentrate on satisfying the following exact properties: (i) use the concept of full exact exchange, as defined by the correct uniform coordinate scaling Levy and Perdew (1985); Levy (1991) (see elaboration below); (ii) freedom from one-electron self-interaction Perdew and Zunger (1981); (iii) correct asymptotic behavior of the xc energy density per particle at  van Leeuwen and Baerends (1994); (iv) reproduction of the homogeneous electron gas limit. In addition, we wish to maintain an overall balanced non-locality of exchange and correlation Perdew et al. (2008).

Regarding property (i), under the uniform coordinate scaling the density transforms as , with its integral, , unchanged and the exchange scales as  Fiolhais et al. (2003); Kümmel and Kronik (2008), which implies . This scaling relation is fulfilled, e.g., by , the exchange energy density per particle in the local spin-density approximation (LSDA).

For the correlation functional, , no such simple scaling rule exists: the correlation scales as , where the superscript indicates a system with an electron-electron interaction that is reduced by a factor of Fiolhais et al. (2003). Additional scaling results for the correlation energy can be found in, e.g., Refs. Levy and Perdew (1985); Levy (1991). Here, we concentrate on the limiting case of high electron densities, i.e. , where the xc energy should be dominated by , Levy (1991)

(3)

A functional is said to use full exact exchange if it obeys Eq. (3) Perdew et al. (2008).

With this definition in mind, we return to Eq. (2). Using , we obtain

(4)

We now see that when scales in the high density limit as with , then it is clear that the first term on the RHS of Eq. (4) is a correlation contribution rather than an exchange term 111Note that there exists a stronger requirement on the correlation energy, namely (see Ref. Levy (1991), Eq.(12)), which here we do not strive to fulfill.. Assuming scales as with , the functional that fulfills this condition can therefore justly be viewed as a combination of EXX, namely, , and a compatible correlation term, .

The reduced density gradient Perdew et al. (1996a)

(5)

where is the Bohr radius, and is the spin polarization, is a natural ingredient to be used to construct a that aims at enforcing the uniform coordinate scaling, because in the high density limit . We make use of this quantity as described in detail below.

Property (ii) is reflected in the equation , where are one-spin-orbital densities, with denoting the -th KS-orbital in the spin-channel . One can attempt to realize such a one-spin-orbital condition by detecting regions of space in which the density is dominated by just one spin-orbital and making sure that full exact exchange and zero correlation is used there. Previous works have discussed the use of iso-orbital indicators Becke (1985); Dobson (1992); Perdew et al. (1999); Kümmel and Perdew (2003) for similar tasks. Here, we define a one-spin-orbital-region indicator by

(6)

where is the von Weizsäcker kinetic energy density and is the Kohn-Sham kinetic energy density.

For one-spin-orbital densities of ground-state character, , because and . For regions with slowly varying density, however, because tends to zero, whereas does not. In contrast to expressions suggested in the past Jaramillo et al. (2003), Eq. (6) does not classify a region of two spatially identical orbitals with opposite spins as a one-orbital region. It also avoids introducing Arbuznikov et al. (2006); Arbuznikov and Kaupp (2007); Bahmann et al. (2007); Kaupp et al. (2007); Perdew et al. (2008); Arbuznikov et al. (2009); Theilacker et al. (2011) any parameters in .

Despite our use of Eq. (6) and the frequent use of similar indicators in the past, we wish to point out two caveats before proceeding. First, it should be noted that formally there exists a difference between one-electron and one-spin-orbital regions. The former correspond to spatial regions in the interacting-electrons system where the probability density is such that one finds just one electron. The latter, however, correspond to spatial regions in the KS system dominated by a single KS spin-orbital Kurth et al. (1999). There is no guarantee that these two regions coincide, because, strictly speaking, the interacting system and the KS system have only the total electron density in common.

Our second caveat refers to the fact that orbital densities are typically not of ground-state character. Therefore the equivalence of and is not guaranteed for these over all space. It is reached, however, in the energetically relevant asymptotic region. We further note that it has recently been pointed out Hofmann and Kümmel (2012) that also the Perdew-Zunger SI correction Perdew and Zunger (1981) may have problems because of orbital densities not being ground-state densities. This may indicate that the question of how to associate orbitals with electrons for the purposes of eliminating self-interaction is a fundamental one, affecting all of the presently used concepts for self-interaction correction that we know of.

With the aim of fulfilling conditions (i) to (iv) we propose the following approximate form for our EXX-compatible correlation energy density per particle, :

(7)

In other words, we approximate the LMF function of Eq. (4) by

(8)

the semi-local exchange energy density per particle by its LSDA form Parr and Yang (1989) , and the semi-local correlation energy density per particle by

(9)

which is the LSDA correlation energy density per particle, multiplied by .

The proposed functional is one-electron SI-free, has the required asymptotic behavior for at , behaves correctly under uniform coordinate scaling, and reduces to the LSDA for regions of slowly varying density.

One-electron self-interaction is addressed via . When tends to , vanishes and the only remaining term is , which then cancels the Hartree repulsion. Note that the semi-local correlation part, which is the last term in Eq. (II), also vanishes for one-spin-orbital regions. This is assured by introducing the prefactor in front of . Otherwise, for one-orbital regions one would get the undesired, unbalanced combination of EXX and local correlation.

The correct uniform scaling is achieved due to the denominator in , which scales as , and cancels the -dependence of the exchange terms that multiply it. In addition, scales as (see Eq. (10) in Ref. Perdew and Wang (1992), Sec. II of Ref. Levy (1991)), which is slower than . Therefore, the limit in Eq. (3) is satisfied.

For slowly varying densities, and , which yields , reproducing the LSDA limit as required.

Finally, note that the proposed approaches the known exact limit at . Since EXX already has the right asymptotic decay of van Leeuwen and Baerends (1994), it suffices to verify that of Eq. (II) decays faster. Indeed, because the orbitals asymptotically tend to , where , the density is dominated by the highest occupied orbital, , and tends to . Because asymptotically and , one finds , which makes decay exponentially. Therefore, the correct asymptotic behavior at is achieved.

There remains one important point to be discussed. In Eq. (II) we are left with one undetermined parameter, . Unfortunately, we presently do not know of an ab initio constraint that would allow us to fix this parameter uniquely, although we do not rule out the possibility that future work may achieve this. The value of affects the amount of EXX that is used in a calculation and is therefore expected to have an influence in practical applications. One can therefore argue that not having determined from first principles is a disadvantage. However, with being a free parameter, the functional form contains some freedom which allows one to adjust it to specific many-electron systems. One can therefore argue that our yet undetermined is in line with the principle of reducing (but not eliminating) empiricism in DFT Burke (2012).

In this first study, the freedom of varying will be used deliberately to explore the properties of the proposed functional. We perform fitting of per system for a representative test set to observe how much its optimal value varies between the different systems, and whether a global fitting procedure, i.e. fitting for all systems combined, is at all justified. In particular, we wish to elucidate the question of whether good binding energies and good eigenvalues can be achieved with the suggested local hybrid functional form. As an aside we note that when is a fixed, system-independent parameter, the proposed functional is fully size-consistent and complications that are known to occur with system-specific adjustment procedures Karolewski et al. (2013) are avoided.

Iii Methods

The proposed functional was implemented and tested using the program package DARSECMakmal et al. (2009); Makmal (2010) an all-electron code, which allows for electronic structure calculations of single atoms or diatomic molecules on a real-space grid represented by prolate spheroidal coordinates. We therefore avoid possible uncertainties associated with the use of pseudopotentials or complicated basis sets in OEP calculations Kümmel and Kronik (2008) – an advantage for accurate functional testing.

DARSEC  allows the user to solve the KS equations self-consistently for density- as well as orbital-dependent functionals (ODFs), for example the proposed functional. For ODFs, the xc potential is constructed by using either the full optimized effective potential formalism (OEP) Kümmel and Kronik (2008); Grabo et al. (1997) via the S-iteration-method Kümmel and Perdew (2003a, b) or, with reduced computational effort, by employing the Krieger-Li-Iafrate (KLI) approximation Krieger et al. (1992). We note that other ways of defining approximations to the OEP exist della Sala and Görling (2001); Gritsenko and Baerends (2001); Ryabinkin et al. (2013). However, for pure exchange earlier works have shown that total energies and eigenvalues are obtained with very high accuracy in the KLI approximation Krieger et al. (1992); Grabo et al. (1997); della Sala and Görling (2001), and for our local hybrid we explicitly compare KLI results to full OEP results in Sec. IV.1 and find very good agreement.

In DARSEC , all computations were converged up to Ry in the total energy, , as well as in the highest occupied KS eigenvalue, , by appropriately choosing the parameters of the real-space grid and by iterating the self-consistent DFT cycle. For full OEP calculations, applying the S-iteration method to the KLI xc potential typically resulted in a reduction of the maximum value of the S-function Kümmel and Perdew (2003a) by a factor of 100. The spin and the axial angular momentum of the systems were taken as in experiment. Note that to this end, for some systems it was necessary to force the KS occupation numbers.

Numerical stability of self-consistent computations using ODFs, in the KLI- or OEP-scheme, mainly depends on the numerical realization of the functional derivative

(10)

which “conveys” the special character of the corresponding xc functional into the calculation of the xc potential. Because our functional approach results in a rather complicated function (see Appendix A, Eq. (52)), careful analytical restructuring was necessary in order to avoid diverging and unstable calculations. In particular, an explicit division by the KS orbitals or the electron density should be avoided, because their exponential decay Kreibich et al. (1998) leads to instabilities at outer grid points. A numerically stable was gained by such considerations, for example by replacing in Eq. (II) with the equivalent expression , or, in case division by the density cannot be avoided, by equally balancing density terms of the same power in numerator and denominator (for details see Appendix A).

All results using (semi-)local functionals (LSDA Perdew and Wang (1992), PBE Perdew et al. (1996a)) or the B3LYP hybrid functional Stephens et al. (1994) (evaluated within the generalized KS scheme Seidl et al. (1996)) were obtained with the Turbomole  program package TUR (), using the def2-QZVPP basis set. The pure EXX calculations were performed in DARSEC  by employing the functional derivative originating from Eq. (1) (as derived in Appendix A, Eq. (17)).

When evaluating a new functional, it is reasonable to concentrate on a class of relatively simple systems to keep computational costs low and to refrain from additional sources of error beyond the xc approximation, e.g., searching for an optimal geometry in systems with many degrees of freedom. However, the systems should not be too simple, so as to pose a significant challenge for the proposed functional. The class of systems has to be large enough, as success or failure for one particular system has very limited meaning. It should also be rich enough to try to represent other systems that are not included. Previous work Perdew et al. (1996a) has shown that a limited set of well selected small molecules can allow for meaningful exploration of a functional’s properties. For these reasons we focus on a set of 18 light diatomic molecules: H, LiH, Li, LiF, BeH, BH, BO, BF, CH, CN, CO, NH, N, NO, OH, O, FH, F, and their constituent atoms. The systems include single-, double-, and triple-bond molecules as well as atoms (no bonding).

Iv Results

iv.1 Comparison of KLI and OEP

While good agreement between the KLI and OEP scheme has been demonstrated before for ground-state energy calculations using EXX Li et al. (1993), the accuracy of the KLI approximation needs to be checked anew when it is applied to a previously untested functional. Table 1 provides this check for our functional. It compares the total energy and the highest occupied KS eigenvalue as obtained with the OEP and the KLI approximation for different values of the parameter (cf. Eq. (II)) for different systems, and lists the corresponding differences for EXX for comparison.

Table 1 shows that the requirement  Kümmel and Kronik (2008) is fulfilled independent of the value of employed. Unlike for the total energy, there is no theorem stating that the highest occupied KS eigenvalue found in the OEP scheme must be below its KLI counterpart. For example, for the C atom and the N molecule, we observe the opposite. Furthermore, because the suggested local hybrid with for spin-unpolarized systems () is exactly equivalent to the purely semi-local constituent functional, one would expect the KLI and OEP results to coincide. This is indeed fulfilled within numerical accuracy. A detailed listing of the total energies and eigenvalues of the highest occupied KS states obtained by the KLI approximation in comparison to full OEP can be found in Appendix B, Tables 4 and  5.

With increasing , a larger amount of EXX is employed and the functional gains more non-local character, leading to greater deviations between KLI and OEP results. Note that, within the considered -range, the deviations with our functional are consequently lower than those obtained for EXX. The last statement applies to both and . Furthermore, in agreement with Ref. Engel and Dreizler (2011) (p. 255), we observe an increasing difference between KLI and OEP results with growing number of electrons in the system.

To summarize, using the KLI approximation for our functional is as justified as it is for pure EXX. This observation is in agreement with the fact that EXX is the limiting case of the suggested functional for .

suggested functional EXX
system
C 0.0000 0.0002 0.0003 0.0004
-0.0005 0.0001 0.0003 0.0007
BH 0.0000 0.0002 0.0005 0.0006
0.0000 0.0003 0.0004 0.0010
Li 0.0000 0.0001 0.0002 0.0002
0.0000 0.0002 0.0005 0.0006
NH 0.0001 0.0005 0.0008 0.0011
0.0007 0.0013 0.0025 0.0055
N 0.0000 0.0009 0.0017 0.0023
0.0000 -0.0010 -0.0019 0.0018
Table 1: Comparison of total energy, , and highest occupied KS eigenvalue, , obtained with the suggested local hybrid functional and with pure EXX, within both the KLI and OEP schemes, as a function of (, ). All values are in Hartree.

iv.2 Fitting the parameter for each system

The proposed functional has one unknown parameter, . We aim to define a global value for , relying on fitting it such that for a group of selected systems, some predefined quantity is optimally predicted (possible choices are discussed in detail below). As a prerequisite, we obtain individual -values by optimizing the parameter for each of the systems separately. As a test for whether a global fitting procedure is meaningful, we verify that these individual -values are clustered within a reasonable numerical range.

In the following, we present two ways to fit . One possibility is fitting the dissociation energy: To find for the molecule AB, the total energies of the molecule and its constituent atoms have to be calculated, with the same . Then, the dissociation energy is fitted to its experimental value Lide (2011), , by varying . Alternatively, one can compute the total energy of the system for various values of and fit it to the experimental total energy. The latter is obtained for atoms as - the sum of all its experimental IPs, ; For molecules as . Unless explicitly stated otherwise, here and throughout molecular properties are calculated at their experimental bond lengths Lide (2011).

Table 2 presents optimized -values for various systems, obtained from both the -fitting and the -fitting procedures. The numerical uncertainty reported for the -values is due to the 1 mRy numerical accuracy in the total energy. The table confirms that the chosen numerical accuracy for the total energy is indeed sufficient. We note that in the -fitting there is a tendency for to increase with the electron number, which reflects a larger contribution of exact exchange. We attribute this to the fact that the energy of the core electrons (which is less important in -fitting) is more strongly dominated by exchange. For our purposes, the most important conclusion to be drawn from Table 2 is that for all systems examined in both approaches, optimal values for lie between 0 and 1, and are never larger than 2. This observation justifies our pursuit of a global value of .

System
H 0.552 0.002 0.537 0.012
LiH 0.642 0.005 0.556 0.004
Li 1.50 0.06 0.571 0.002
LiF 0.141 0.006 0.976 0.003
BeH 0.746 0.025 0.648 0.004
BH 0.590 0.010 0.685 0.004
BO 0.288 0.007 0.916 0.002
BF 0.578 0.027 0.943 0.002
CH 0.672 0.028 0.741 0.003
CN 0.146 0.005 0.908 0.002
CO 0.283 0.009 0.916 0.002
NH 0.667 0.027 0.811 0.004
N 0.107 0.009 0.908 0.003
NO 0.329 0.009 0.960 0.002
OH 1.20 0.07 0.942 0.004
O 0.472 0.009 1.004 0.002
FH 0.075 0.011 1.105 0.004
F 0.356 0.006 1.206 0.003
H any
Li 0.543 0.005
Be 0.644 0.005
B 0.698 0.003
C 0.757 0.002
N 0.848 0.005
O 0.925 0.004
F 1.067 0.003
Table 2: The parameter optimized for various systems, using the - and -fitting procedures.

iv.3 Determining a global value for the parameter

Following the conclusion that the parameter can indeed be fitted, we performed a series of calculations, obtaining the -dependent average relative errors

(11)

Here, can refer to the dissociation energy, , the total energy, , or the ionization potential evaluated via , the IP-theorem for the exact functional. The index runs over all the systems calculated 222The quantities and were obtained relying on all the molecules and atoms in the reference set (), while was obtained relying on the molecules only ()..

Figure 1: Average relative error of the dissociation energy, , as a function of the parameter (solid line). Relative errors for the LSDA (dashed), PBE (dash-dotted) and B3LYP (dotted) functionals are given for comparison. Pure EXX reaches an error of and exceeds the scale we chose here.

The functions and are plotted in Figs. 1 and 2, respectively, accompanied by the average relative errors for commonly used functionals: the LSDA, PBE, and B3LYP. As mentioned previously, the B3LYP results here and in the following were obtained in the generalized KS approach, which we, based on previous experience Kümmel and Kronik (2008), expect to yield total energies that are very similar to the ones from the KS approach for the systems studied here. For completeness, results obtained with pure EXX evaluated in the KLI approximation are also reported.

In both figures we observe clear minima for the proposed functional at the values of for and for , with minimal error values of and , respectively. These error values are close to those achieved with the B3LYP functional, and are significantly better than the PBE and LSDA results. Because optimizing and demands almost the same value for , a satisfying description of both properties is possible using a common parameter of . For this , the relative error in the dissociation energy is lowest for the BF molecule (0.7%) and highest for Li and F (14% and 17%, respectively). The relative error in the total energy is more evenly spread around 0.12%.

Figure 2: Average relative error of the total energy, , as a function of the parameter (blue solid line). Relative errors for the LSDA (dashed), PBE (dash-dotted) and B3LYP (dotted) functionals, as well as pure EXX(KLI) (purple solid line) are given for comparison.

The function shown in Fig. 3 exhibits a different behavior, reaching its minimum of 6 % at a higher value of . 333When calculating , the vertical experimental ionization potentials were used (see Ref. Lide (2011) and http://webbook.nist.gov) When evaluated at , the average relative error is . Although lower than 31% for B3LYP and 42% for both LSDA and PBE, such a deviation is rather significant. Therefore, Fig. 3 suggests that in order to reach good agreement between the experimental IP and , a larger amount of EXX is required.

Interestingly, when calculating NH and BO, we observed that the highest occupied state changes with varying the parameter from to at approximately for NH, and from to at for the BO molecule. Such systems could therefore be good candidates for checking the functional’s ability to predict physically meaningful orbitals in the sense of Ref. Dauth et al. (2011).

Figure 3: Average relative error of the IP predicted via the highest occupied KS eigenvalue, , as a function of the parameter (blue solid line). Relative errors for the LSDA (dashed), PBE (dash-dotted) and B3LYP (dotted) functionals, as well as pure EXX(KLI) (purple solid line) are given for comparison.

Last, we checked the previously made assumption that experimental bond lengths can be used, assuming they are not very different from those obtained by relaxation. To this end, all 18 molecules in the reference set were relaxed, and the obtained bond lengths were compared to the experimental values,  Lide (2011). It was found that for most systems lies within the computational error for and the difference is below 0.02 Bohr 444The numerical error in is governed by the accuracy of 1 mRy in the total energy, rather than by the convergence of the relaxation process., except F, where Bohr. 555Two exceptional cases are LiH and Li, which have an extremely shallow minimum. The uncertainty of 1 mRy in the total energy translates into a numerical uncertainty in the bond length of 0.16 Bohr and 0.26 Bohr, respectively. Therefore, the difference , being 0.04 Bohr and 0.22 Bohr, although large, has no actual meaning due to the large numerical uncertainty.

iv.4 Achievements of the suggested functional

In the following, we examine some of the proposed functional’s properties at the value of , which was determined in Sec.IV.3.

As the functional is one-electron SI-free (see Sec. II), it is important to investigate its behavior in systems that are known to suffer from a large self-interaction error when described by standard DFAs. First, for one-electron systems, like H, He, H, etc. the functional reduces analytically to the EXX functional, as can be seen from Eq. (II). Therefore, all the properties of these systems are obtained, by construction, exactly. This advantage is not shared by (semi-)local or most hybrid functionals.

In particular, Fig. 4 presents the dissociation curve of H, obtained with various functionals. It can be seen that, as expected, the curve obtained with the proposed functional perfectly agrees with the EXX curve, which provides the exact result in this case. In particular, our local hybrid does not exhibit a spurious maximum in the curve, which appears in conventional approximate functionals at bond lengths around 5-6 Bohr Bally and Sastry (1997); Ruzsinszky et al. (2005); Cohen et al. (2008); Livshits and Baer (2008); Nafziger and Wasserman (2013) and whose electrostatic origin has recently been discussed Dwyer and Tozer (2011). The dissociation of neutral H is a special challenge for most density functionals and is closely connected to the question of how static correlation is accounted for Cohen et al. (2012). Our local hybrid for H yields a binding curve that is qualitatively similar to the one obtained in pure exchange calculations, i.e., for large internuclear separation the lowest energy is obtained with a spin-polarized atomic density centered around each nucleus, which yields a total energy of 1 Hartree. Quantitatively, there are differences with respect to the EXX solution: The point from which on the spin-polarized solution has a lower energy than the spin-unpolarized one lies at about 3.6 Bohr with our local hybrid, and the minimum energy is -1.173 Hartree as compared to -1.134 Hartree obtained with EXX.

Figure 4: Dissociation curve of the H molecule, for the LSDA (squares), PBE (x’s), B3LYP (dots), EXX (solid line) and the suggested functional (circles).

Generally, delocalization in stretched molecular bonds is conceptually connected to the SI-error of DFAs Cohen et al. (2008). Therefore, reduction of this error marks a first step towards enhancing the description of dissociation processes and chemical reactions Ruzsinszky et al. (2006).

To examine this, the dissociation curve of the 3-electron molecule He is shown in Fig. 5. The curve achieved with the suggested functional is the closest to the reference result obtained with a highly accurate wavefunctions method (see Ref. Ruzsinszky et al. (2005) and references therein). Here again, only our local hybrid and the EXX curves do not possess the spurious maximum, which appears in the conventional approximations – LSDA, PBE and B3LYP around 4 Bohr. Unlike for the H system, here the proposed functional does not automatically reduce to the exact expression, and therefore the accurate prediction for He in Fig. 5 can be seen as a consequence of the strong reduction of SI-errors, in agreement with a previous study. Körzdörfer et al. (2008)

Figure 5: Dissociation curve of the He molecule, for the LSDA (squares), PBE (x’s), B3LYP (dots), EXX (rhombi), and the suggested functional with (circles), compared to CCSD(T) results from Ref. Ruzsinszky et al. (2005) (dashed line)

We further investigated how well corresponds to the experimental IP Lide (2011) for the atoms Li, Na, and K. These atoms can be considered as quasi-one-electron systems, consisting of electrons arranged in closed shells, which screen the charge of the nucleus, and one additional electron in the last open shell. Table 3 shows that the obtained from our functional evaluated with is closer to the experimental IP than the from LSDA, PBE, and B3LYP. Note that for these systems a remarkable improvement is achieved, as one would expect from their strong one-electron character.

IP
system Exp. LSDA PBE B3LYP suggested functional
Li 0.1163 0.1185 0.1311 0.1797
(-41 %) (-40 %) (-34 %) (-9 %)
Na 0.1131 0.1116 0.1251 0.1647
(-40 %) (-41 %) (-34 %) (-13 %)
K 0.0961 0.0930 0.1038 0.1334
(-40 %) (-42 %) (-35 %) (-16 %)
Table 3: The highest occupied eigenvalue compared to the experimental IP for Li, Na, and K, computed with four different functionals (LSDA, PBE, B3LYP and our suggested functional using ). The table contains the absolute numbers in Hartree as well as the relative error in %.

V Conclusions

In this article we presented the construction of a local hybrid functional that combines full exact exchange with compatible correlation. The functional respects the homogeneous electron gas limit and addresses the one-electron self-interaction error via a one-spin-orbital-region indicator. The qualitative improvement that is achieved with this construction is reflected in, e.g., dissociation energy curves for H and He that are much more realistic than the ones obtained from (semi-) local functionals and global hybrids. We investigated different conditions for fixing the undetermined parameter of the functional. When the parameter is fit to minimize binding energy errors or total energy errors, with respect to experiment, then our local hybrid reaches an accuracy that is better than LSDA or PBE and similar to the one afforded by the B3LYP global hybrid. Predicting the first ionization energy via the highest occupied eigenvalue is more accurate with our functional than with LSDA, PBE or B3LYP, but still not satisfactorily accurate.

When the parameter is fit such that should be as close as possible to the experimental first ionization potential, then the local hybrid functional achieves much smaller errors in this quantity than, e.g., B3LYP. However, the value obtained for the free parameter differs considerably from the one that was obtained by fitting to binding or total energies. As a result, prediction of these energies considerably differs from their experimental values. Therefore, our local hybrid does allow for reaching accurate binding energies or accurate highest eigenvalues, but not with the same functional parametrization.

Looking at this from a more general perspective, we note that many functionals can achieve good accuracy on one of the aforementioned properties or the other, but not on both properties at the same time.

A first example are global hybrids. With the usual 0.2 to 0.25 fraction of exact exchange they yield good binding energies, but highest eigenvalues that are considerably too small in magnitude. Increasing the aforementioned fraction to leads to improved gap prediction Sai et al. (2011); Imamura et al. (2011); Atalla et al. (2013). However, such a large fraction of exact exchange can compromise significantly the accuracy in thermochemical Adamo and Barone (1999); Ernzerhof and Scuseria (1999); Zhao and Truhlar (2008) or electronic structure properties. Marom et al. (2010); Körzdörfer and Kümmel (2010); Kronik et al. (2012)

A second example are range-separated hybrid functionals. When combined with a tuning procedure based on the IP theorem Stein et al. (2010); Refaely-Abramson et al. (2011, 2012); Kronik et al. (2012); Körzdörfer et al. (2011); Sini et al. (2011); Foster and Wong (2012); Phillips et al. (2012); Risko and Brédas (2013), they allow for obtaining eigenvalues that reflect ionization potentials very accurately by construction. However, the typical value of the range-separation parameter that is reached by such tuning is quite different from the one that is reached when atomization energy errors are minimized via the range-separation parameter Livshits and Baer (2007).

A third example is provided by the various self-interaction correction schemes. Different forms of self-interaction correction greatly improve the interpretability of the eigenvalues when compared to semi-local functionals Perdew and Zunger (1981); Körzdörfer et al. (2009); Pederson et al. (1985); Klüpfel et al. (2012); Ullrich et al. (2000); Körzdörfer et al. (2008); Vydrov and Scuseria (2005), but binding energies are not accurately predicted Vydrov and Scuseria (2005); Hofmann et al. (2012) unless the correction is “scaled down” Vydrov et al. (2006); Klüpfel et al. (2012).

This rather universal difficulty to achieve accurate eigenvalues and accurate binding energies at the same time may indicate that combining these two properties may require a type of physics that all present day functionals lack Verma and Bartlett (2012).

Recent work provides two new and interesting perspectives on this problem. On the one hand, it has been noted that even a semi-local functional can yield eigenvalues that are qualitatively similar to the ones obtained from bare EXX when the asymptotic properties of a GGA are carefully determined Armiento and Kümmel (2013). On the other hand, it has recently been shown that an ensemble perspective offers new and improved ways of interpreting eigenvalues and extracting information from semi-local functionals Kraisler and Kronik (2013). Exploring in particular this latter option, i.e., combining the ensemble approach with the present local hybrid functional, will be the topic of future work and may shed further light on the question of how to obtain accurate binding energies and Kohn-Sham eigenvalues from the same functional.

Acknowledgements.
S.K. gratefully acknowledges discussions with J. P. Perdew on local hybrids in general and on an early version of this functional in particular. We thank Baruch Feldman for fruitful discussions. Financial support by the DFG Graduiertenkolleg 1640, the European Research Council, the German-Israeli Foundation, and the Lise Meitner center for computational chemistry is gratefully acknowledged. E.K. is a recipient of the Levzion scholarship. T.S. acknowledges support from the Elite Network of Bavaria (“Macromolecular Science” program).

Appendix A Derivation of the correlation potential

In order to employ the OEP formalism Grabo et al. (1997); Kümmel and Kronik (2008) (or its KLI approximation Krieger et al. (1992)), one has to provide an analytical expression for the functional derivative of the explicitly orbital-dependent exchange-correlation energy, , with respect to the orbitals . For the functional proposed in the present contribution,

(12)

where is the exact exchange defined in Eq.(1), is the semi-local correlation energy, whose energy density per particle, , is given in Eq. (9), and equals

(13)

with the LMF function, , being defined in Eq. (8).

Due to the additive structure of Eq. (12), the functional derivative

(14)

can be split in three terms:

(15)

which are considered separately in the following.

a.1 The exact exchange contribution

The first term on the RHS of Eq. (15) can be computed directly from the exact-exchange expression (Eq. (1))

(16)

and simply reads

(17)

a.2 The semi-local correlation contribution

The self-interaction-free semi-local correlation energy contribution is defined by

(18)

where

(19)

and

(20)

For completeness, we list out all the quantities that are required to construct the function :

  1. kinetic energy density

    (21)
  2. Von Weizsäcker kinetic energy density

    (22)
  3. spin polarization

    .

Taking the functional derivative based on Eq. (18) results in two parts:

(23)

By denoting the constituent functions of the function by , and , chain rule arguments lead to the following expression:

(24)

Here we explicitly took into account the fact that depends on locally.

In order to obtain an analytical expression for , which can then be inserted into Eq. (23), one has to consider the three constituent functions separately:

(25)
(26)

(27)
(28)

(29)
(30)

Here the operator denotes a gradient relative to the coordinate and the quantity distinguishes between the two spin channels by

(31)

It is noted that the functional derivatives above have the presented form with respect to the occupied orbitals only; derivatives with respect to unoccupied orbitals equal zero.

The derived relations (25) - (30) now have to be inserted via Eq. (24) into Eq. (23). By further employing a chain rule argument for the second term on the RHS of Eq. (23)

one arrives at the final expression of the functional derivative of :

(33)

Equation (33) corresponds to the functional derivative the way it was implemented into the KLI/OEP-routine in the program package DARSEC.

a.3 The contribution

The correlation term is defined in Eq. (13). Let us denote

(34)

and recall that the LMF function equals

(35)

In addition to the quantities , and introduced above, the function additionally employs the so-called reduced density gradient Perdew et al. (1996a)

(36)

with

(37)

and

The exact relation

(38)

will be useful for later derivations.

Analogously to Eq. (23), the application of the functional derivative with respect to the KS-orbitals leads to two contributions:

(39)

Moreover, an analogous relation to Eq. (24) helps to rewrite the first part of this equation, only that now one has to consider also the functions and :

(40)

We evaluate each term separately and obtain:

(41)

(42)

(43)

(44)
(45)

(46)