One-electron self-interaction and the asymptotics of the Kohn–Sham potential: an impaired relation

Abstract

One-electron self-interaction and an incorrect asymptotic behavior of the Kohn–Sham exchange–correlation potential are among the most prominent limitations of many present-day density functionals. However, a one-electron self-interaction-free energy does not necessarily lead to the correct long-range potential. This is shown here explicitly for local hybrid functionals. Furthermore, carefully studying the ratio of the von Weizsäcker kinetic energy density to the (positive) Kohn–Sham kinetic energy density, τ_W/τ, reveals that this ratio, which frequently serves as an iso-orbital indicator and is used to eliminate one-electron self-interaction effects in meta-generalized-gradient approximations and local hybrid functionals, can fail to approach its expected value in the vicinity of orbital nodal planes. This perspective article suggests that the nature and consequences of one-electron self-interaction and some of the strategies for its correction need to be reconsidered.

---

1 Density functional approximations and their Kohn–Sham potentials

During the past few decades, Kohn–Sham density-functional theory (DFT)^1,2 has evolved into a standard tool for electronic structure calculations of atoms, molecules and solids. The decisive quantity of DFT is the exchange–correlation (xc) energy functional, E_xc, which contains all electronic interactions beyond the classical electrostatic Hartree contribution, E_H. E_xc in practice has to be approximated, and the approximation used governs the accuracy of a DFT calculation.^3,4 It is one of the puzzles of DFT that explicit density functionals such as the generalized gradient approximations (GGAs) can predict binding energies and bond lengths of complex many-electron systems reliably, but make substantial errors in describing simple one-electron systems. The underlying problem is well-known as the one-electron “self-interaction problem”:⁵ for the exact functional, E_xc + E_H will vanish for any one-electron ground-state density because one electron does not interact with itself – but most approximate functionals yield a spurious finite value for this case. Following ref. 5 a functional is considered to be one-electron self-interaction free if it fulfills the condition

E_H[n_iσ] + E_xc[n_iσ] = 0, (1)

where n_iσ = |φ_iσ(r)|² designates a single spin-orbital density.

Self-interaction plays a decisive (although not the only) role in the (un)reliability of density functional theory calculations, and its consequences are particularly pronounced, e.g., in questions of orbital localization,^5–8 ionization processes,^9–12 charge transfer,^13–15 and the interpretability of eigenvalues and orbitals, e.g., as photoemission observables.^16–23

Many of these observables can also be directly related to properties of the Kohn–Sham exchange–correlation potential, which is defined as the functional derivative of the xc energy with respect to the ground-state density n(r), i.e., . It is generally expected that there is a close relation between freedom from self-interaction and xc potential features. The field-counteracting term that is important for obtaining correct response properties is one such feature.^24,25 Another example, and probably the most prominent one, is the long-range asymptotic behavior of the xc potential,^26,27

(Hartree units are used here and throughout.) In this perspective article we focus exclusively on Kohn–Sham theory, i.e., on a local multiplicative xc potential, as opposed to orbital-specific (non-multiplicative) potentials that arise in generalized Kohn–Sham theory²⁸ and are used in the standard application of hybrid functionals. In the Kohn–Sham approach, the local xc potential models the interaction of one particle with all others and it therefore appears intuitively plausible that a functional that is not self-interaction-free cannot show the correct −1/r long-range asymptotic behavior: as one particle of a finite, overall electrically neutral system ventures out to infinity, it will “feel” the hole of charge 1 that it left behind in the total charge. This gives rise to the −1/r potential asymptotics (see, e.g., ref. 4, p. 242 for a more detailed argument along these lines). However, a particle that spuriously self-interacts will “feel” itself, and thus not the proper hole. Consequently, the potential will not have the proper long-range decay.

The correct asymptotics of the xc potential has proven to be important for a variety of physical quantities. It plays a prominent role for obtaining stable anions in DFT, it leads to a Rydberg series in the Kohn–Sham eigenvalues and generally to unoccupied eigenvalues of improved interpretability, and as a consequence allows for improved accuracy in the prediction of various response properties.^29–33 The correct asymptotic behavior is also important for the ionization potential (IP) theorem,^26,27,34,35 which states that the negative of the highest occupied Kohn–Sham eigenvalue −ε_ho should correspond to the vertical IP, and for developing functionals that allow for approximately predicting IPs from ground-state eigenvalues.^36,37

There have been fruitful attempts to incorporate the correct behavior in the limit |r| → ∞ directly into the xc potential,^38–41 leading to improvements in the description of some of the aforementioned properties. However, since directly designed potential expressions are typically not functional derivatives of any energy functional, the use of such “potential only” approximations is necessarily limited, as discussed, e.g., in detail in ref. 42–44.

A functional that combines freedom from self-interaction and the correct asymptotics of the potential is exact exchange (EXX), being defined as the Fock integral evaluated using Kohn–Sham orbitals φ_iσ(r), where i labels orbitals in spin channel σ:

Here, N_σ is the number of electrons with spin σ. Treating exchange exactly with a local Kohn–Sham potential leads to a significant improvement in the quality of Kohn–Sham eigenvalues when compared to (semi-)local functionals.^20,45,46 EXX also tends to increase Kohn–Sham gaps,^8,47–50 leads to a desired particle number discontinuity in static⁵¹ and time-dependent⁵² situations, and improves the description of charge transfer,^24,25 dissociation⁵³ and ionization processes.⁵²

However, using bare EXX is known for its rather poor description of binding energies and structural properties (see, e.g., ref. 54 and 55, and ref. 4, chapter 2). Adding a (semi-)local correlation term to EXX hardly improves the situation and typically leads to results that are inferior to the ones from (semi-)local functionals. The reason for this failure is the well-known incompatibility of the fully non-local Fock exchange with a purely (semi-)local correlation term.⁵⁶

A class of approximations which has been designed to remedy this incompatibility is that of local hybrid functionals,^57,58 sometimes also called hyper-GGAs.⁵⁶ Whereas global hybrid functionals^59–64 mix a constant, fixed fraction of Fock exchange with (semi-)local exchange and correlation, local hybrids replace the fixed fraction by a density dependent local mixing function (LMF). Both types of hybrids originate from the concept of the coupling-constant integration, i.e., the adiabatic connection scheme.^59,65 Global hybrids are successful in modeling the coupling-constant averaged, integrated energy. Local hybrids can go one step further and aim to model the coupling-constant curve itself⁶⁶ instead of just the integral. Thus, in contrast to the global hybrid functionals which are used in practical applications of DFT and combine GGA components with a fixed fraction of exact exchange, local hybrids can incorporate full exact exchange and can be fully one-electron self-interaction-free.

An early local hybrid with reduced one-electron self-interaction error showed promising results for dissociation curves and reaction barriers, but its accuracy for binding energies was limited.⁵⁸ A self-consistent implementation of a local hybrid functional was given in ref. 67, and over the years several local hybrids were constructed, using different LMFs and (semi-)local exchange and correlation functionals,^68–76 striving to reach greater accuracy by refining the position-dependent mixing of nonlocal and local components. Many of these functionals rely on the concept of an iso-orbital indicator, i.e., a functional that allows one to distinguish regions of space in which the density is dominated by one orbital shape from regions of space where several orbitals of different shapes contribute to the density. The most prominent iso-orbital indicator, which goes back to a long tradition of using kinetic energy densities in density functional construction,^77–79 is the ratio of the von Weizsäcker kinetic energy density τ_W to the positive (as opposed to other possible definitions, see, e.g., ref. 80) Kohn–Sham kinetic energy density τ, discussed in detail below.

By using full EXX and an iso-orbital indicator, local hybrids aim at being one electron self-interaction-free and producing a Kohn–Sham potential with the proper long-range asymptotic decay. They are a paradigm class of functionals designed for simultaneously curing both of these two prominent problems of (semi-)local density functionals. In the following, we therefore use the example of a local hybrid functional to shed light on the relation between a functional's self-interaction and its potential asymptotics, as well as the properties of the τ_W/τ indicator. We argue that quite generally a one-electron self-interaction-free energy does not guarantee the correct long-range potential, and that τ_W/τ loses its indicator ability in the vicinity of nodal planes of the highest-occupied molecular orbital (HOMO).

2 Correlation compatible with exact exchange: the local hybrid approach

The xc energy functional can be written as

E_xc[n] = ∫n(r)e_xc([n];r)d³r, (4)

with e_xc([n];r) denoting the xc energy density per particle. The definition of e_xc(r) is not unique and subject to a gauge-dependence.⁷² Yet, for local hybrid functionals it has become common to define this energy in the form

e^lh_xc(r) = e^ex_x(r) + f(r)(e^sl_x(r) − e^ex_x(r)) + e^sl_c(r). (5)

Here, e^ex_x(r) marks the exchange energy density per particle corresponding to the EXX energy of eqn (3). This non-local term is mixed with (semi-)local exchange and correlation energy densities e^sl_x(r) and e^sl_c(r), respectively. The position dependent mixing ratio f(r), which is itself a density functional, marks the LMF.

Often, the LMF is designed in a way that aims at eliminating the one-electron self-interaction error of eqn (1) that is inherent in most density functionals. An established method for reducing self-interaction effects is to detect regions of space where a single Kohn–Sham orbital shape dominates the density (“iso-orbital regions”), and then enforce eqn (1) in these regions. One of the most popular^{58,67,70,71,74,76,80–83} indicator functions for detecting iso-orbital regions is

where τ_W(r) = |∇n(r)|²/(8n(r)) denotes the von Weizsäcker kinetic energy density and is the positive Kohn–Sham kinetic energy density. In iso-orbital regions, τ(r) → τ_W(r) and therefore g(r) → 1. In the case of a slowly varying density, τ_W(r) → 0 and, since τ(r) remains finite, g(r) → 0. This indicator function is typically a decisive ingredient in the LMF, f(r), of local hybrids. With its help one can construct f(r) such that eqn (5) reduces to correct limiting cases, e.g., e^sl_x(r) + e^sl_c(r) for slowly varying densities, and e^ex_x(r) for single orbital regions. The latter case additionally requires that e^sl_c(r) vanishes in single-orbital regions, a condition that we discuss below.

In the asymptotic limit, |r| → ∞, the xc energy density for a finite system should be dominated by e^ex_x(r). When e^sl_c(r) vanishes sufficiently fast in the asymptotic region (a condition that is usually fulfilled), then

is the requirement that one aims at, because it leads to the correct asymptotic limit of the xc energy density per particle

(Note the difference to the asymptotic limit of the xc potential, see ref. 38).

Since for a finite system each Kohn–Sham orbital decays exponentially with an exponent set by its eigenvalue,⁸⁴ the density is asymptotically dominated by the HOMO density, i.e., it becomes of iso-orbital character. Therefore, g(r) can be used in the construction of the LMF to realize eqn (8).

Considerations of the type discussed above are inherent to many density functional constructions. As a particular example for a local hybrid functional we here use a recently proposed, physically motivated LMF,⁸³ which reads

The function g(r) in the numerator is multiplied by the squared spin polarization ζ(r) = (n_↑(r) − n_↓(r))/(n_↑(r) + n_↓(r)), which lets the LMF not only identify iso-orbital regions, but also correctly distinguish between true one-orbital regions, and regions with two identical spin-orbitals. The function g(r) is used in such a way that f_t(r) vanishes for one-orbital regions, as required. The use of the reduced density gradient

where a₀ is the Bohr radius and Φ(ζ(r)) = ½((1 + ζ)^2/3 + (1 − ζ)^2/3), in the denominator of f_t(r), ensures the correct behavior of E_xc under uniform coordinate scaling r → γr.^85,86 The density transforms as n_γ(r) = γ³n(γr) and as a consequence eqn (9) uses full exact exchange in the sense of ref. 72

The function t²(r) is multiplied by a parameter c that we cannot determine, at least presently, from fundamental constraints. It allows for adjustments in the functional ansatz. In the case of slowly varying densities, f_t(r) → 1 and eqn (5) reduces to its purely (semi-)local components. As an aside we note that this LMF comprises the one of ref. 58 as the special case when c = 0 and ζ(r) = 1 ∀ r. We denote this case by f₀(r), i.e., .

For the semi-local exchange we use the LSDA,⁸⁷i.e., e^sl_x(r) = e^LSDA_x(r), whereas . The additional multiplication with the numerator of eqn (9) consistently reduces eqn (5) to pure EXX in the one-spin-orbital case, where e^LSDA_c(r) alone does not vanish.

The general questions that we discuss in this perspective article, i.e., whether there is a relation between self-interaction and the xc potential asymptotics and how far the iso-orbital indicator τ_W/τ can be used to enforce freedom from self-interaction, can be scrutinzed with the local hybrid of eqn (9) as an instructive example.

3 The Kohn–Sham exchange–correlation potential of local hybrid functionals

In order to implement local hybrids self-consistently within the Kohn–Sham scheme, one has to find the local multiplicative xc potential corresponding to the energy of eqn (4) and (5). The fact that local hybrids use EXX and typically also τ(r) makes them explicitly orbital-dependent. Therefore, the local xc potential must be obtained from the optimized effective potential (OEP) equation (see, e.g., ref. 45, 55 and 88–91). The computational effort can be reduced significantly by employing the approximation of Krieger, Li and Iafrate (KLI).⁹² For the local hybrid of eqn (9) it has been shown that the total energy E_tot and the highest occupied Kohn–Sham eigenvalue ε_ho obtained using the KLI approximation agree quite well with the ones from the full OEP.⁸³ Furthermore, it is a general finding⁸⁴ that the KLI approximation does not affect the potential asymptotics to leading order. In the actual calculations presented in the following we therefore always use the KLI approximation.

In the OEP (and KLI) scheme the chain rule for functional derivatives⁴⁵ relates the derivative with respect to the density to the derivatives with respect to the orbitals,

From the structure of the OEP equation it further follows that to first order

i.e., the functional derivative with respect to the HOMO in general determines the potential asymptotics.⁸⁴ Therefore, investigating the HOMO functional derivative is the key to determining the asymptotic behavior of an orbital dependent functional's xc potential. When one takes the functional derivative (with respect to the orbital) of a local hybrid one obtains three terms, corresponding to the three addends in eqn (5):

u^lh_iσ(r) = u^exx_iσ(r) + u^c-nl_iσ(r) + u^c-sl_iσ(r) (14)

Evaluating the asymptotical behavior of each of these three terms for the highest occupied orbital allows one to predict the potential asymptotics.

The first term can be derived directly from eqn (3) and reads

This term evaluated for the HOMO indeed provides the correct asymptotic behavior⁴⁵

The third term u^c-sl_iσ(r), on the other hand, does not contribute to the asymptotics of eqn (16) as it decays exponentially due to its purely (semi-)local nature.

Evaluating the second term on the right-hand side of eqn (14) requires careful consideration. Intuitively, one might expect that an asymptotically vanishing LMF will surpress any asymptotic contribution of this term to the potential. In the following we check this expectation. Details of the underlying calculation for both LMFs used in this work, i.e., f_t(r) and f₀(r), can be found in ref. 83 and in Appendix B, eqn (27), respectively.

By defining P(r) = n(r)e^sl_x(r) and Q(r) = n(r)e^ex_x(r) one can write

The first two terms consist of (semi-)local components and thus vanish exponentially. Evaluating the third term on the other hand is not as trivial as it contains the non-local quantity Q(r) as well the functional derivative of the LMF with respect to the corresponding Kohn–Sham orbital. For the LMFs addressed in this perspective we find that this term does not contribute to the asymptotics either (see ref. 83 and Appendix B for details).

Thus, only the fourth term in eqn (17) is relevant in the asymptotic limit and therefore

The first term in this equation equals −u^exx_iσ(r) of eqn (15), locally multiplied by f(r)/2. Due to eqn (7) it vanishes faster than the leading term of u^lh_iσ, which is given in eqn (16).

The second term, however, is of a different structure, as it evaluates the LMF under the integral. By considering the HOMO level, its asymptotic limit is

This corresponds to a Hartree-like potential caused by the spin-orbital density of the HOMO averaged over all space, with the LMF as a weighting function. Thus, this term gives a finite contribution in the asymptotic limit despite eqn (7).

Now, when adding the asymptotically significant components, eqn (16) and (19), for the evaluation of eqn (13), we arrive at

Here, the parameter γ_σ denotes the reduced slope of the potential asymptotics, which can numerically be extracted from a self-consistent Kohn–Sham calculation via

Eqn (21) is a central result of this work, as it demonstrates that a local hybrid of the form of eqn (5) does not lead to the exact asymptotic behavior of the xc potential. Eqn (20) holds for all f(r) that vanish in the asymptotic limit and for which the third term of eqn (17) does not contribute to the asymptotics of the functional derivative u^c-nl_iσ(r), i.e., under very general conditions. Further details of the calculation, specifically regarding the question of the xc potential asymptotics in different spin-channels, are given in Appendix B.

The LMF is limited between 0 ≤ f(r) ≤ 1 and therefore the asymptote is bound between ½ < γ_σ ≤ 1. Consequently, the exact value γ_σ = 1 can only be reached by setting f(r) = 0 ∀ r, which corresponds to the trivial case of using EXX, “as is” or combined with a purely (semi-)local correlation functional.

A different extreme case, f(r) = 1 ∀ r, does not, as one could naïvely believe due to eqn (21), lead to γ_σ = ½. Here, we have to take the neglected first term of eqn (18) into account again, and from this we see that γ_σ actually vanishes. This is to be expected, since in this case the local hybrid reduces to a purely (semi-)local functional.

Fig. 1 shows a numerical verification of the above analytical considerations (see Appendix A for numerical details). It depicts the xc (KLI) potential corresponding to the local hybrid of eqn (9) in comparison with the asymptotic decay according to eqn (20) and (21) for the carbon atom. An additonal curve indicates the exact −1/r decay, which is clearly not reached. The xc potential, instead of decaying with γ_↑ = 1 as one would intuitively expect,⁶⁹ approaches the predicition of eqn (21) (γ_↑ = 0.716) quite rapidly.

Fig. 1 The xc potential v_xc↑(r) for the C atom along the x-axis (see Appendix A for definition), computed using f_t(r) with c = 0.5. Also displayed is the asymptotic curve according to eqn (20) with γ_↑(c = 0.5) = 0.716 and the correct asymptotic −1/r. The inset shows the potential plotted along the complete axis.

Fig. 2 shows the xc energy density e_xc(r) for the same system in comparison to its correct asymptotic of −1/(2r). Clearly the xc energy density shows the correct asymptotic, cf.eqn (8). We thus see that while the behavior of e_xc can directly be controlled via the LMF in eqn (5), the process of finding the local xc potential via functional differentiation leads to non-local evaluations of the LMF that decisively impact the potential's asymptotics.

Fig. 2 xc energy density per particle e_xc(r) for the C atom along the x-axis (see Appendix A for definition), computed using f_t(r) with c = 0.5. Also displayed is the corresponding asymptotic slope of −1/(2r). The inset shows the energy density plotted along the complete axis.

A physically meaningful quantity closely related to the asymptotics of the xc potential is the highest occupied eigenvalue ε_ho. Table 1 shows −ε_ho compared to the experimental IP for the carbon and fluorine atoms for different functionals, together with the corresponding value of γ_σ of the xc potential from eqn (21).

Table 1 Comparison of the highest occupied Kohn–Sham eigenvalue −ε_ho to the experimental vertical IP⁹³ for the C (ε_4↑) and F atoms (ε_4↓) using different functionals. All values are in hartrees

System	Functional	γ σ ho	−εho	Exp. IP
C	LSDA	—	0.2249	0.4138
	f t (r)(c = 0)	0.6098	0.2740
	f t (r)(c = 0.5)	0.7162	0.3067
	f t (r)(c = 1.0)	0.7678	0.3302
	f t (r)(c = 2.5)	0.8441	0.3688
	f t (r)(c = 5.0)	0.8966	0.3970
	f 0(r)	0.8309	0.3530
	EXX	1.0000	0.4378
F	LSDA	—	0.3808	0.6403
	f t (r)(c = 0)	0.5055	0.3810
	f t (r)(c = 0.5)	0.6665	0.4724
	f t (r)(c = 1.0)	0.7390	0.5269
	f t (r)(c = 2.5)	0.8365	0.6060
	f t (r)(c = 5.0)	0.8971	0.6570
	f 0(r)	0.7927	0.5798
	EXX	1.0000	0.6779

---

The LSDA, as generally known, significantly underestimates the IP due to the wrong potential asymptotics and the inherent self-interaction error. Using pure EXX with the correct asymptotic decay and no self-interaction error leads to a much better prediction of the IP. When employing a local hybrid with the LMF f_t(r), the explicit dependence on the parameter c becomes evident: with growing c, the asymptotic value γ_σ grows and the description of the IP improves. Fig. 3 sheds further light on the situation. It shows potentials of local hybrids which are all based on eqn (9) but use different values of c. Growing values of c increase the amount of EXX and lead to an overall deeper potential. This explains that the eigenvalues become more negative.

Fig. 3 Asymptotics of the xc potential v_xc↑(r) for the C atom along the x-axis, computed with pure EXX and local hybrids using f_t(r) from eqn (9) with parameters c = 0.5 and c = 5.0. Also displayed are the complete potentials in the inset.

We thus see that while all of the local hybrids used here can (so far, see caveat in the next section) be thought of as being one-electron self-interaction free, they show different potential asymptotics and their highest occupied eigenvalues predict the IP with significantly different reliability. The relation between freedom from self-interaction, potential asymptotics, and physical interpretability of the highest occupied eigenvalue as the negative IP is therefore much less clear than intuitively believed. This observation also calls for taking a closer look at the iso-orbital indicator g(r) that is used in enforcing freedom from self-interaction. This is the topic of the next section.

4 The implications of orbital nodal planes

As explained in the preceding sections, many local hybrids and other functionals such as meta-GGAs rely on the function g(r) tending to 1 to detect regions of space in which a single orbital shape dominates the density, and then, e.g., correct for self-interaction in such regions. However, a first caveat that one has to take note of is that g(r) → 1 holds for one-particle densities of ground-state character. This is a possibly far reaching restriction for the use of g(r) because electron orbital densities typically have nodes, i.e., are not of ground-state character. As a specific, illustrative example, consider an atom where the HOMO has an azimuthal quantum number m i.e., is expressed in spherical coordinates as φ_ho(r) = R(r,θ)e^imϕ. In the region where the density is HOMO-dominated, τ_W(r) is therefore given approximately by ½|∇R(r,θ)|². However, in the same region τ(r) is given approximately by ½|∇φ_ho(r)|² = ½(|∇R(r,θ)|² + m²R²(r,θ)/(r sin θ)²). Thus, if m = 0, then τ(r) = τ_W(r), but for m ≠ 0 this is no longer the case. Instead, τ(r) and τ_W only approach each other asymptotically, as the m²-dependent-term of τ(r) decays to zero with large r.

One may counter-argue that this restriction is not so severe because in density functional construction, the condition τ(r) → τ_W(r) is mostly used to detect those regions of space in a finite system which are far from all nuclei, and where the density decays nodelessly. However, we here show that even in such regions the condition τ(r) → τ_W(r) can be violated. This leads to a second caveat about the reliability of the g(r) indicator. It is rooted in the existence of orbital densities that have nodal planes or nodal axes. Fig. 4 illustrates this case. It shows g(r) evaluated in the (xz)-plane for the carbon atom (see Appendix A for numerical details, including grid setup). The density here was obtained using f_t(r)(c = 0.5), but the density features relevant here are not sensitive to functional details. The important observation is that g(r) approaches 1 in the asymptotic limit in every direction – except for in the vicinity of x = 0.

Fig. 4 The function on the numerical grid for the C atom.

The first step towards an understanding of this finding is to note that the z-axis is a nodal axis, being the intersection of the nodal planes of the two HOMOs of carbon, which are degenerate and of p-orbital character.

The consequences of the existence of nodal planes can be studied analytically. To this end we look at a schematic density that is dominated by the HOMO φ_ho(r), but also take the next lower lying orbital φ_ho−1(r) into account, i.e. n(r) ∼ |φ_ho(r)|² + |φ_ho−1(r)|². With this ansatz one finds

and

τ ∼ ½(|∇φ_ho|² + |∇φ_ho−1|²). (23)

These two terms combined and evaluated on or close to a nodal plane (denoted by ), where φ_ho → 0, yield

Even though φ_ho vanishes on the nodal plane, its gradient still yields a finite value and keeps the function g(r) from approaching 1.

Fig. 4 shows that the deviation from 1 has a noticeable spatial extension of a few a.u. This raises the question of how well the use of the iso-orbital indicator g(r) leads to freedom from self-interaction, as in some regions that so far have been considered as iso-orbital ones, e.g., all space far from the system's center, self-interaction effects may not be eliminated fully when the indicator aberrates due to the presence of a nodal plane or axis. A different interpretation of Fig. 4 would be to reconsider one's expectation of where iso-orbital regions are, or what they are. The traditional point of view has been that all space far from a finite system's center is of iso-orbital nature. Fig. 4 and eqn (24) may be interpreted to show that this is not the case when the HOMO has a nodal plane/axis extending to infinity. From this perspective one might say that g(r) does exactly what it is supposed to be doing, i.e., it indicates that the nodal plane region is not of iso-orbital character. Yet, also from this perspective Fig. 4 reveals a surprising finding, namely that even infinitely far from a finite system's center, the density may not be of iso-orbital character.

The nodal plane observation also forces us to take a yet closer look at the central topic of this perspective, the potential asymptotics. Nodal planes can influence the asymptotics of a local hybrid's xc potential in two ways. First, it has been argued that all orbital-dependent functionals show non-vanishing asymptotic constants in their xc potential along nodal planes of the highest occupied Kohn–Sham orbital that extend to infinity. This was first discussed in ref. 94 and 95 for the case of pure EXX, and the occurring shift was determined to be

C_σ = v̄_xcMσσ − ū_xcMσσ, (25)

with v̄_xciσ = ∫φ_iσ*(r)v_xcσ(r)φ_iσ(r)d³r and ū_xciσ = ∫φ_iσ*(r)u_xciσ(r)φ_iσ(r)d³r. The index M_σ denotes the highest lying Kohn–Sham orbital that does not show a vanishing spin-orbital density along the nodal plane of the HOMO. Since eqn (25) follows from the KLI (OEP) equation without referring to a specific functional, non-vanishing asymptotic constants on nodal planes of the HOMO are expected on rather general grounds.

Second, the fact that g(r) → 1 is not guaranteed on a nodal plane can also affect the potential. For the sake of clarity, we again discuss this effect for the specific example of local hybrids. When the LMF tends to zero on the nodal plane, i.e., and eqn (7) is obeyed, then the non-vanishing constant of eqn (25) is the only effect. An example for this case is the LMF f_t(r) with a finite value of the parameter c. It is depicted in Fig. 5 for the C atom density in the (xz)-plane, and one sees that there are no asymptotic features. This is because the reduced density gradient in the denominator causes f_t(r) to vanish in the asymptotic limit, regardless of the occurrence of a nodal plane. The potential decays like −γ_σ/r in all directions, but along the z-axis a non-vanishing constant

appears. This is shown in Fig. 6 for f_t(r)(c = 0.5) and the F atom. One can clearly see how v_xc↑(r) decays with γ_↑ = 0.6650, but, instead of zero, approaches a constant of C_↑ = 0.0244, in agreement with eqn (25).

Fig. 5 The LMF , evaluated with c = 0.5, on the numerical grid for the C atom.

Fig. 6 Asymptotics of the xc potential v_xc↑(r) for the F atom along the (projected) z-axis (denoted z*, see Appendix A for definition), computed using f_t(r) with the parameter c = 0.5.

A different situation occurs when , i.e., the behavior of the indicator function along a nodal plane/axis of the HOMO prevents the LMF from reaching its intended limit. This happens, e.g., for f₀(r) or f_t(r)(c = 0) and is depicted in Fig. 7, again for the C atom. The occurrence of a nodal axis here very clearly affects the LMF. Since in this case eqn (7) is violated in the direction of the z-axis, the previous derivations cannot be used to predict the potential's asymptotic behavior. However, we have numerically checked the xc potential's behavior. On the nodal axis it neither tends to −1/r, nor to , but rather tends to some other value. Thus, the nodal axis in this case has a very noticeable influence on the potential asymptotics, which is hard to predict a priori.

Fig. 7 The LMF on the numerical grid for the C atom.

5 Conclusions

With local hybrid functionals serving as an explicit example we have argued that freedom from self-interaction in the sense of eqn (1) does not necessarily lead to the expected −1/r decay of the local Kohn–Sham xc potential. We have further argued that the ratio of the von Weizsäcker kinetic energy density to the positive Kohn–Sham kinetic energy density, which is frequently used in functional construction for indicating iso-orbital regions and eliminating self-interaction effects in these, may not serve its intended purpose because it is very sensitive to excited state features such as orbital nodal planes that are present in Kohn–Sham orbitals that construct ground-state densities of many-electron systems.

These findings have immediate and somewhat discomforting consequences for the local hybrid approach. For a large class of functionals one has to accept that the correct long-range xc potential simply cannot be obtained. This observation plays a role in explaining why it is very hard to construct a local hybrid that yields good binding energetics and physically meaningful eigenvalues with the same functional form and set of parameters.⁸³ However, the impaired relation between self-interaction and the xc potential's asymptotics, and also the impact of nodal planes, stand in a context that is much larger than the local hybrid one. The iso-orbital indicator g(r) has been used in many functionals, not only local hybrids. Nodal planes are known to impact the exact exchange potential in surprising ways.^94,95 They have appeared here as a prominent feature in kinetic energy ratios, and we expect⁹⁶ that they play a much larger role in the exchange potential than has been realized so far. The observation that a one-electron self-interaction-free energy can go together with a potential that does not fall of like −1/r is not only a feature of local hybrids, but has also been reported for a “scaled down” version of the Perdew–Zunger self-interaction correction.⁹⁷ One may therefore wonder whether semi-local indicator functionals are in some sense incompatible with the fully non-local self-interaction correction that is achieved by EXX or full Perdew–Zunger-type correction approaches. It has also been pointed out recently⁹⁸ that eqn (1) itself, which is the basis of the present definition of one-electron self-interaction, leads to questions when evaluated for orbital densities, because E_xc[n] is intended to be used with ground state densities, whereas orbital densities are excited state densities. Further conceptual questions about eqn (1) relate to its inherent identification of orbitals with electrons and its unitary variance.^5,99,100 The success of self-interaction corrections schemes that rely on eqn (1) tells us that the equation is meaningful. However, the sum of the insights into its limitations that emerged over the years suggests that there is more to the question of self-interaction in density functional theory.

While the above considerations point out areas that require further thought and work, one should also note that there have been developments in DFT that shine a bright light into the future. The concept of many-electron self-interaction^101,102 is not as straightforward to use as eqn (1), but it avoids the conceptual questions that are associated with this equation. Range-separated hybrids yield the correct asymptotic potential and have proven to be a very successful concept, without being self-interaction-free.^{23,36,103–117} There have been successful functional constructions that can be seen as combinations of the local hybrid and the range-separation idea.^118,119 Ensemble corrections¹²⁰ allow to extract information from functionals in an unexpected way, and can, e.g., further improve IP prediction. Finally, it has recently been shown¹²¹ that a new type of a generalized gradient approximation can show features that were so far thought of as being associated only with exact exchange, such as step structures and surprising nodal plane features,⁹⁶ and understanding potentials in terms of xc charges has provided new insights.^122–124 Therefore, the battle against DFT's old foe, the self-interaction error, and its surprisingly independent side-kick, the wrong potential fall-off, is far from being lost.

Appendix A: numerical details

We used the all-electron code DARSEC¹²⁵ for all calculations presented in this perspective. This code exploits the rotational symmetry of diatomic molecules along the interatomic axis z, treating the azimuthal angle ϕ analytically and thus effectively reducing the problem of solving the Kohn–Sham equations in two dimensions. The equations are represented on a real-space grid of prolate-spheroidal coordinates. In such a coordinate system, the nuclear position(s) coincide with the focal point(s) of the grid located at z = ±R/2, with R being the bond length of the diatomic molecule. This is the case also for calculations of single atoms: the position of the nucleus is not equivalent to the origin of the coordinate system, but it is located at z = −R/2, where R was set to 0.5 a.u. e.g., the C atom in our plots is centered at r_C = (x,z) = (0, −0.25). The x-axis is defined as perpendicular to the z-axis, crossing the latter at z = 0, i.e., at a point being equidistant from the focal points of the grid (see ref. 125 for details).

In order to avoid numerical instabilities due to singularities in the Laplacian, the grid was chosen such that it does not include the actual z-axis, i.e. the interatomic axis. As a consequence, in this direction all quantities can only be plotted along a projected z*-axis, which takes into account all grid points that are closest to the actual z-axis. Since the discrepancy between the projected and the real z-axis decreases with increasing number of grid points, we made sure that the difference between z and z* is small by choosing sufficiently dense and large grids.

Appendix B: the asymptotic decay of the exchange–correlation potential in detail

In the following, we present considerations about the asymptotics of the xc potential in the spin channel that carries the global HOMO (σ_ho), as compared to the other spin channel ([small sigma, Greek, macron]_ho). Section 3 used the condition that f(r) needs to vanish at a sufficient rate in the derivation of eqn (20). In the present work, we investigated two possibilities for the decay of the LMF.

First, f_t(r) for a finite value of the parameter c vanishes exponentially because . In this case, all individual terms in each functional derivative, u^c-nl_iσ (see eqn (17)), vanish exponentially in the asymptotic limit as well, except for the second term in eqn (18). Eventually, this remaining term is responsible for the reduced asymptotic decay of eqn (20) due to the non-local evaluation of f_t(r). Consequently, the xc potential in both spin channels decays with −γ_σ/r.

However, a different picture emerges when evaluating . This function decays much more slowly than f_t(r) with finite c, as Fig. 8 shows for the carbon atom. Consequently, not all terms in the functional derivative originating from f₀(r) vanish individually and more detailed investigations are necessary.

Fig. 8 Comparison of f₀(r) and f_t(r) with c = 0.5 for the C atom along the x-axis.

Defining K(r) = n(r)(e^sl_x(r) − e^ex_x(r)), the functional derivative in this case reads

with . Therefore, both reach the same absolute value in the asymptotic limit, but show opposite signs. Now, we have to distinguish between the spin channels: If one looks at in the spin channel that has the global HOMO, i.e. n(r) ∼ |φ_Nσhoσho(r)|², then one can see from eqn (27) that the fourth and fifth terms are equivalent in the asymptotic limit except for the sign. Therefore, they cancel each other and, since the first and second terms decay fast enough, only the third term remains, leading to the limit of −γ_σho/r. In the other spin channel however, the fourth and fifth terms do not cancel anymore, since the density is still dominated by φ_Nσhoσho(r), whereas the fourth term features φ_{N[small sigma, Greek, macron]ho}[small sigma, Greek, macron]_ho(r). Therefore, in the other spin channel yet another asymptotic limit is obtained, again strictly following from the evaluation of the functional derivative.

This feature can be corrected by using a spin-polarized ansatz with an indicator function that is a spin-polarized LMF of the form , with τ_Wσ(r) and τ_σ(r) being the kinetic energy spin densities. In this case, a functional derivative that does not feature the total density n(r) follows and therefore the aforementioned effect does not occur. However, since for the spin channel σ_ho all derivations made are valid independently of the form of the LMF and since this spin channel features the physical meanigful quantity −ε_ho, it suffices for this work to consider the more simple LMFs instead of their spin-polarized counterparts.

Acknowledgements

We acknowledge financial support from the German–Israeli Science foundation. T.S. acknowledges support from the Elite Network of Bavaria (“Macromolecular Science” program). S.K. acknowledges support from Deutsche Forschungsgemeinschaft. E.K. is a recipient of the Levzion scholarship. L.K. acknowledges support from the Lise Meitner Minerva Center for Computational Chemistry.

References

1. P. Hohenberg and W. Kohn, Phys. Rev. Sect. B, 1964, 136, 864–871.
2. W. Kohn and L. Sham, Phys. Rev. Sect. A, 1965, 140, 1133.
3. S. Kurth and J. P. Perdew, Int. J. Quantum Chem., 2000, 77, 814–818.
4. A Primer in Density Functional Theory, ed. C. Fiolhais, F. Nogueira and M. A. Marques, Springer, 2003, vol. 620.
5. J. P. Perdew and A. Zunger, Phys. Rev. B: Condens. Matter Mater. Phys., 1981, 23, 5048–5079.
6. W. M. Temmerman, Z. Szotek and H. Winter, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 11533.
7. P. Strange, A. Svane, W. M. Temmerman, Z. Szotek and H. Winter, Nature, 1999, 399, 756.
8. E. Engel and R. Schmid, Phys. Rev. Lett., 2009, 103, 036404.
9. X.-M. Tong and S.-I. Chu, Phys. Rev. A: At., Mol., Opt. Phys., 1997, 55, 3406.
10. S.-I. Chu, J. Chem. Phys., 2005, 123, 062207.
11. C. A. Ullrich, P.-G. Reinhard and E. Suraud, Phys. Rev. A: At., Mol., Opt. Phys., 2000, 62, 053202.
12. D. A. Telnov, T. Heslar and S.-I. Chu, Chem. Phys., 2011, 391, 88.
13. T. Körzdörfer, M. Mundt and S. Kümmel, Phys. Rev. Lett., 2008, 100, 133004.
14. D. Hofmann, T. Körzdörfer and S. Kümmel, Phys. Rev. Lett., 2012, 108, 146401.
15. D. Hofmann and S. Kümmel, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 86, 201109(R).
16. P. Duffy, D. Chong, M. E. Casida and D. R. Salahub, Phys. Rev. A: At., Mol., Opt. Phys., 1994, 50, 4707.
17. D. P. Chong, O. V. Gritsenko and E. J. Baerends, J. Chem. Phys., 2002, 116, 1760.
18. A. Pohl, P.-G. Reinhard and E. Suraud, Phys. Rev. A: At., Mol., Opt. Phys., 2004, 70, 023202.
19. T. Körzdörfer, S. Kümmel, N. Marom and L. Kronik, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 201205(R).
20. T. Körzdörfer and S. Kümmel, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 155206.
21. M. Dauth, T. Körzdörfer, S. Kümmel, J. Ziroff, M. Wiessner, A. Schöll, F. Reinert, M. Arita and K. Shimada, Phys. Rev. Lett., 2011, 107, 193002.
22. P. Klüpfel, P. M. Dinh, P.-G. Reinhard and E. Suraud, Phys. Rev. A: At., Mol., Opt. Phys., 2013, 88, 052501.
23. L. Kronik and S. Kümmel, in First Principle Approaches to Spectroscopic Properties of Complex Materials, ed. C. D. Valentin, S. Botti and M. Cococcioni, Springer, Berlin, to be published.
24. S. J. A. van Gisbergen, P. R. T. Schipper, O. V. Gritsenko, E. J. Baerends, J. G. Snijders, B. Champagne and B. Kirtman, Phys. Rev. Lett., 1999, 83, 694.
25. S. Kümmel, L. Kronik and J. P. Perdew, Phys. Rev. Lett., 2004, 93, 213002.
26. M. Levy, J. P. Perdew and V. Sahni, Phys. Rev. A: At., Mol., Opt. Phys., 1984, 30, 2745–2748.
27. C.-O. Almbladh and U. von Barth, Phys. Rev. B: Condens. Matter Mater. Phys., 1985, 31, 3231–3244.
28. A. Seidl, A. Görling, P. Vogl, J. A. Majewski and M. Levy, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 53, 3764.
29. D. J. Tozer and N. C. Handy, J. Chem. Phys., 1998, 109, 10180.
30. D. J. Tozer, R. D. Amos, N. C. Handy, B. O. Ross and L. Serrano-Andres, Mol. Phys., 1999, 97, 859.
31. M. E. Casida and D. R. Salahub, J. Chem. Phys., 2000, 113, 8918.
32. M. A. L. Marques, A. Castro and A. Rubio, J. Chem. Phys., 2001, 115, 3006.
33. A. D. Teale and D. J. Tozer, Chem. Phys. Lett., 2004, 59, 383.
34. J. P. Perdew, R. G. Parr, M. Levy and J. L. Balduz, Phys. Rev. Lett., 1982, 49, 1691–1694.
35. J. P. Perdew and M. Levy, Phys. Rev. B: Condens. Matter Mater. Phys., 1997, 56, 16021–16028.
36. L. Kronik, T. Stein, S. Refaely-Abramson and R. Baer, J. Chem. Theory Comput., 2012, 8, 1515.
37. P. Verma and R. J. Bartlett, J. Chem. Phys., 2012, 136, 044105.
38. R. van Leeuwen and E. J. Baerends, Phys. Rev. A: At., Mol., Opt. Phys., 1994, 49, 2421–2431.
39. D. J. Tozer and N. Handy, J. Chem. Phys., 1998, 109, 10180.
40. A. D. Becke and E. R. Johnson, J. Chem. Phys., 2006, 124, 221101.
41. W. Cencek and K. Szalewicz, J. Chem. Phys., 2013, 139, 024104.
42. A. P. Gaiduk, S. K. Chulkov and V. N. Staroverov, J. Chem. Theory Comput., 2009, 5, 699.
43. A. Karolewski, R. Armiento and S. Kümmel, J. Chem. Theory Comput., 2009, 5, 712.
44. A. Karolewski and S. Kümmel, Phys. Rev. A: At., Mol., Opt. Phys., 2013, 88, 052519.
45. T. Grabo, T. Kreibich and E. K. U. Gross, Mol. Eng., 1997, 7, 27.
46. F. della Sala and A. Görling, J. Chem. Phys., 2001, 115, 5718.
47. M. Städele, M. Moukara, J. A. Majewski, P. Vogl and A. Görling, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 10031.
48. A. Makmal, R. Armiento, E. Engel, L. Kronik and S. Kümmel, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 80, 161204.
49. E. Engel, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 80, 161205.
50. M. Betzinger, C. Friedrich and S. Blügel, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 88, 075130.
51. J. B. Krieger, Y. Li and G. J. Iafrate, Phys. Rev. A: At., Mol., Opt. Phys., 1992, 45, 101.
52. M. Mundt and S. Kümmel, Phys. Rev. Lett., 2005, 95, 203004.
53. A. Makmal, S. Kümmel and L. Kronik, Phys. Rev. A: At., Mol., Opt. Phys., 2011, 83, 062512.
54. E. Engel, A. Höck and R. Dreizler, Phys. Rev. A: At., Mol., Opt. Phys., 2000, 62, 042502.
55. E. Engel and R. Dreizler, Density Functional Theory: An Advanced Course, Springer, 2011.
56. J. P. Perdew and K. Schmidt, Density Functional Theory and Its Application to Materials, Melville, NY, 2001.
57. F. G. Cruz, K.-C. Lam and K. Burke, J. Phys. Chem. A, 1998, 102, 4911–4917.
58. J. Jaramillo, G. E. Scuseria and M. Ernzerhof, J. Chem. Phys., 2003, 118, 1068.
59. A. D. Becke, J. Chem. Phys., 1993, 98, 1372.
60. A. D. Becke, J. Chem. Phys., 1993, 98, 5648.
61. J. P. Perdew, M. Ernzerhof and K. Burke, J. Chem. Phys., 1996, 105, 9982.
62. P. J. Stephens, F. J. Devlin, C. F. Chabalowski and M. J. Frisch, J. Phys. Chem., 1994, 98, 11623–11627.
63. C. Adamo and V. Barone, J. Chem. Phys., 1999, 110, 6158.
64. M. Ernzerhof and G. E. Scuseria, J. Chem. Phys., 1999, 110, 5029–5036.
65. M. Ernzerhof, J. P. Perdew and K. Burke, Int. J. Quantum Chem., 1997, 64, 285.
66. P. Mori-Sánchez, A. J. Cohen and W. Yang, J. Chem. Phys., 2006, 124, 091102.
67. A. V. Arbuznikov, M. Kaupp and H. Bahmann, J. Chem. Phys., 2006, 124, 204102.
68. B. G. Janesko and G. E. Scuseria, J. Chem. Phys., 2007, 127, 164117.
69. A. V. Arbuznikov and M. Kaupp, Chem. Phys. Lett., 2007, 440, 160–168.
70. H. Bahmann, A. Rodenberg, A. V. Arbuznikov and M. Kaupp, J. Chem. Phys., 2007, 126, 011103.
71. M. Kaupp, H. Bahmann and A. V. Arbuznikov, J. Chem. Phys., 2007, 127, 194102.
72. J. P. Perdew, V. N. Staroverov, J. Tao and G. E. Scuseria, Phys. Rev. A: At., Mol., Opt. Phys., 2008, 78, 052513.
73. R. Haunschild, B. G. Janesko and G. E. Scuseria, J. Chem. Phys., 2009, 131, 154112.
74. A. V. Arbuznikov, H. Bahmann and M. Kaupp, J. Phys. Chem. A, 2009, 113, 11898–11906.
75. R. Haunschild and G. E. Scuseria, J. Chem. Phys., 2010, 133, 134116.
76. K. Theilacker, A. V. Arbuznikov, H. Bahmann and M. Kaupp, J. Phys. Chem. A, 2011, 115, 8990–8996.
77. A. D. Becke, Int. J. Quantum Chem., 1985, 27, 585.
78. J. F. Dobson, J. Phys.: Condens. Matter, 1992, 4, 7877.
79. A. D. Becke, J. Chem. Phys., 1998, 109, 2092.
80. S. Kümmel and J. P. Perdew, Mol. Phys., 2003, 101, 1663.
81. J. Perdew, S. Kurth, A. Zupan and P. Blaha, Phys. Rev. Lett., 1999, 82, 2544.
82. J. M. Tao, J. P. Perdew, V. N. Staroverov and G. E. Scuseria, Phys. Rev. Lett., 2003, 91, 146401.
83. T. Schmidt, E. Kraisler, A. Makmal, L. Kronik and S. Kümmel, J. Chem. Phys., 2014, 140 DOI:10.1063/1.4865942.
84. T. Kreibich, S. Kurth, T. Grabo and E. K. U. Gross, Adv. Quantum Chem., 1998, 33, 31.
85. M. Levy and J. P. Perdew, Phys. Rev. A: At., Mol., Opt. Phys., 1985, 32, 2010.
86. M. Levy, Phys. Rev. A: At., Mol., Opt. Phys., 1991, 43, 4637.
87. R. G. Parr and W. Yang, Density-Functional Theory of Atoms and Molecules, Oxford University Press, New York, 1989.
88. R. T. Sharp and G. K. Horton, Phys. Rev., 1953, 90, 317.
89. J. D. Talman and W. F. Shadwick, Phys. Rev. A: At., Mol., Opt. Phys., 1976, 14, 36.
90. V. Sahni, J. Gruenebaum and J. P. Perdew, Phys. Rev. B: Condens. Matter Mater. Phys., 1982, 26, 4371–4377.
91. S. Kümmel and L. Kronik, Rev. Mod. Phys., 2008, 80, 3–60.
92. J. Krieger, Y. Li and G. Iafrate, Phys. Rev. A: At., Mol., Opt. Phys., 1992, 46, 5453–5458.
93. CRC Handbook of Chemistry and Physics, ed. D. R. Lide, CRC, London, 92nd edn, 2011.
94. F. della Sala and A. Görling, Phys. Rev. Lett., 2002, 89, 33003.
95. S. Kümmel and J. P. Perdew, Phys. Rev. B: Condens. Matter Mater. Phys., 2003, 68, 035103.
96. T. Aschebrock, T. Schmidt and S. Kümmel, unpublished results.
97. O. A. Vydrov, G. E. Scuseria, J. P. Perdew, A. Ruzsinszky and G. I. Csonka, J. Chem. Phys., 2006, 124, 094108.
98. D. Hofmann and S. Kümmel, J. Chem. Phys., 2012, 137, 064117.
99. T. Körzdörfer, M. Mundt and S. Kümmel, J. Chem. Phys., 2008, 129, 014110.
100. D. Hofmann, S. Klüpfel, P. Klüpfel and S. Kümmel, Phys. Rev. A: At., Mol., Opt. Phys., 2012, 85, 062514.
101. P. Mori-Sánchez, A. J. Cohen and W. Yang, J. Chem. Phys., 2006, 125, 201102.
102. A. Ruzsinszky, J. P. Perdew, G. I. Csonka, O. A. Vydrov and G. E. Scuseria, J. Chem. Phys., 2006, 125, 194112.
103. A. Savin and H.-J. Flad, Int. J. Quantum Chem., 1995, 56, 327–332.
104. T. Leininger, H. Stoll, H.-J. Werner and A. Savin, Chem. Phys. Lett., 1997, 275, 151–160.
105. T. Yanai, D. P. Tew and N. C. Handy, Chem. Phys. Lett., 2004, 393, 51.
106. J.-W. Song, T. Hirosawa, T. Tsuneda and K. Hirao, J. Chem. Phys., 2007, 126, 154105.
107. E. Livshits and R. Baer, Phys. Chem. Chem. Phys., 2007, 9, 2932.
108. J.-D. Chai and M. Head-Gordon, J. Chem. Phys., 2008, 128, 084106.
109. M. A. Rohrdanz, K. M. Martins and J. M. Herbert, J. Chem. Phys., 2009, 130, 054112.
110. U. Salzner and R. Baer, J. Chem. Phys., 2009, 131, 231101.
111. T. Stein, H. Eisenberg, L. Kronik and R. Baer, Phys. Rev. Lett., 2010, 105, 266802.
112. S. Refaely-Abramson, R. Baer and L. Kronik, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 075144.
113. T. Körzdörfer, J. S. Sears, C. Sutton and J.-L. Brédas, J. Chem. Phys., 2011, 135, 204107.
114. A. Karolewski, T. Stein, R. Baer and S. Kümmel, J. Chem. Phys., 2011, 134, 151101.
115. M. A. Rohrdanz and J. M. Herbert, J. Chem. Phys., 2008, 129, 034107.
116. L. Pandey, C. Doiron, J. S. Sears and J.-L. Brédas, Phys. Chem. Chem. Phys., 2012, 14, 14243–14248.
117. S. Refaely-Abramson, S. Sharifzadeh, N. Govind, J. Autschbach, J. B. Neaton, R. Baer and L. Kronik, Phys. Rev. Lett., 2012, 109, 226405.
118. B. G. Janesko, A. V. Krukau and G. E. Scuseria, J. Chem. Phys., 2008, 129, 124110.
119. R. Haunschild and G. E. Scuseria, J. Chem. Phys., 2010, 132, 224106.
120. E. Kraisler and L. Kronik, Phys. Rev. Lett., 2013, 110, 126403.
121. R. Armiento and S. Kümmel, Phys. Rev. Lett., 2013, 111, 036402.
122. S. V. Kohut and V. N. Staroverov, J. Chem. Phys., 2013, 139, 164117.
123. N. I. Gidopoulos and N. N. Lathiotakis, J. Chem. Phys., 2012, 136, 224109.
124. X. Andrade and A. Aspuru-Guzik, Phys. Rev. Lett., 2011, 107, 183002.
125. A. Makmal, S. Kümmel and L. Kronik, J. Chem. Theory Comput., 2009, 5, 1731–1740.

---