On the piecewise convex or concave nature of ground state energy as a function of fractional number of electrons for approximate density functionals Free

Kohn-Sham density functional theory (KS-DFT) is perhaps the most widely used theoretical tool for ground state calculations of quantum systems.^1,2 Compared with wave-function theory, which solves an eigenvalue problem of the many-body Hamiltonian, and the wave-function which is defined on the 3N dimensional Hilbert space, with N being the number of electrons, density functional theory greatly reduces the problem to the minimization of a functional of the electron density, which is defined on the three dimensional Hilbert space, thus the computational complexity is tremendously simplified.

In doing this, one does not sacrifice the theoretical rigor for the integer-N systems. However, within such a formulation, it is critical to consider the extension to the fractional-N systems.^3–6 Within the N-particle Hilbert space and for a given external potential $v(𝒓)$⁠, the many-body ground state energy is explicitly dependent on the electron number N, which is restricted to the set of non-negative integers ℤ. As a consequence, the energy as a function of the number of electrons, $Ev$(N), is only defined on discrete integer points. In wavefunction theory, this is not a problem since a wavefunction associated with a non-integer number of electrons is not physically significant. In density functional theory, however, the electron densities have the freedom to integrate to a fractional number; furthermore, it has been realized in the past decades that the performance of approximate functionals are directly related to their behavior for fractional electron systems.^4–10 Thus it is most natural to consider the electron densities associated with an arbitrary non-negative electron number.

It is possible to extend to the Fock space to describe systems with a fractional electron number N, as established in the Perdew-Parr-Levy-Balduz (PPLB) formulation of the fractional N-systems:³ $Ev$(N) is a continuous function (a curve) of N. Moreover, for the exact functional, $Ev$(N) at N = M (M being an arbitrary non-negative integer) agrees with the wave-function theory in the Hilbert space, and in the interval [M, M + 1] is a straight line interpolation between $Ev$(M) and $Ev$(M + 1). Furthermore, it is assumed that $Ev$(N) is convex between integers, i.e., $Ev(M−1)+Ev(M+1)>2Ev(M)$⁠.^3,11 By this convexity assumption in conjunction with the PPLB condition^3,12–14 that $Ev$(N) is a series of straight line segments across integers, the convexity of the E vs N curve is still true, in the sense that $Ev(N1)+Ev(N2)≥2Ev(N1+N22)$ for any non-negative real numbers N₁ and N₂.

In practice, however, the exact functional is not known. With density functional approximations, the above mentioned convexity no longer holds true. In particular, by restricting to [M, M + 1], i.e., $Ev$(N) between neighboring integers, the PPLB condition is often violated by approximate functionals.^6–8,10 Here in KS-DFT, the calculation procedure for fractional M + n (⁠$0<n<1$⁠) systems can be modified from the integer calculations. Instead of imposing the idempotency condition of the KS reduced density matrix $ρs(𝒓,𝒓′)$⁠, we allow $ρs(𝒓,𝒓′)=∑i=1M+1niψi(𝒓)ψi(𝒓′)$⁠, where n_i = 1 for lower-lying orbitals (⁠$i=1,2,…,M$⁠) everywhere and n_i = n for the frontier orbital (i = M + 1), corresponding to an ensemble of determinants.

By implementing the fractional system calculations, it has been observed that semi-local density functionals display piecewise convex curves, while the Hartree Fock (HF) functional gives piecewise concave curves interpolating between integers.^5,7,15 Some of the consequences of such piecewise non-straight $Ev$(N) curves for real world integer systems are as follows. First, in the case of dissociating molecules, such as the dissociation of homo-nuclear diatomic cations, the energy of the symmetric dissociation product $Z0.5+⋯Z0.5+$ will be underestimated by semi-local functionals while overestimated by HF. In the case of heteronuclear diatomic dissociation, semi-local functionals may predict fictitious fractional charges on the atoms at large separation.^4,7–9 Second, in the process of electron ionization or affiliation for a system composed of multiple identical subsystems that are infinitely separated, the added or removed electron will tend to delocalize over as many subsystems as possible by semi-local density functionals, while localize at one subsystem by the HF functional. This is known as the delocalization error for semi-local density functionals and localization error for HF.⁵ They are the direct consequences of the E vs N curve; the behavior of the electron density follows the energetically most favorable, but physically incorrect pattern.

Furthermore, the non-straight E vs N curve on [M, M + 1] has implications on the energy gap errors. It has been proved that for finite and for bulk systems, in the KS scheme for explicit and continuous functionals of density, or the generalized Kohn-Sham (GKS) scheme for explicit and continuous functionals of density matrix, the chemical potential, or the energy derivative with respect to the number of electrons equals the frontier orbital energy in the corresponding KS or GKS, i.e., $∂E∂N=ϵf$⁠.¹⁵ For integer N systems, f refers to the highest occupied molecular orbital (HOMO) when the derivative is taken from the left side or the lowest unoccupied molecular orbital (LUMO) when the derivative is evaluated from the right side. It follows that $ϵHOMO=∂E∂N|N=M−$ and $ϵLUMO=∂E∂N|N=M+$⁠. For the exact functional which satisfies the PPLB condition, the frontier orbital energies should agree with the minus integer ionization potential and minus electron affinity, respectively, i.e., $ϵHOMO=−I≡E(M)−E(M−1)$ and $ϵLUMO=−A≡E(M+1)−E(M)$⁠. However, if the E vs N curve of semi-local density functionals is convex between neighboring integer points, the HOMO energy will overestimate −I while the LUMO energy will underestimate −A, i.e., $ϵHOMO>−I$ and $ϵLUMO<−A$⁠, and hence the band gap given by the difference of $ϵLUMO$ and $ϵHOMO$ will be underestimated. In contrast, piecewise concave functionals such as HF exhibit exactly opposite behaviors because of their concavity. For bulk systems, $Ev$(N) is linear in fractional numbers of electrons, even for approximate functionals and thus the discontinuity in chemical potentials, which is given by the energy gap in the KS orbitals for explicit and continuous functionals of density or by the energy gap in the GKS orbitals for explicit and continuous functionals of the density matrix, is the prediction of the fundamental bandgaps from the approximate functionals, which is typically underestimated.⁵ 

Although the piecewise concavity of HF and the piecewise convexity of semi-local functionals in terms of the E vs N curve are of vital significance, the conclusions are only made by numerical observations and have never been proved before. In this paper, we rigorously prove that HF is indeed piecewise concave, while the piecewise convexity of semi-local functionals is not generally guaranteed for all fractions. In Secs. II and III, we will first provide a mathematical proof for the HF functional and then perform analysis for the semi-local functionals using the example of local density approximation (LDA) with the exchange functional only to explain why the convexity cannot be proved. Finally we present some numerical examples which serve as a validation for our conclusion for the HF functional and also as counterexamples for the semi-local functionals.

For simplicity, in the following context, by concave (convex) we mean piecewise concave (convex) interpolating between integers. All the numerical calculations were performed using our in-house built QM4D program,¹⁶ with the basis set aug-cc-pVTZ throughout the calculations in the present study.

For simplicity of notations, in the following proof we suppress the spin variables. Given an external potential $v(𝒓)$⁠, we start by considering a system with a fractional number of electrons N = M + n₀ (M is an integer and $0<n0<1$⁠). Without the loss of generality, let us consider an arbitrary density functional in the GKS scheme, where the ground state energy is defined by the following functional minimization over admissible KS density matrix,

Here the admissible set follows the Aufbau principle, i.e.,

where we assume that the orbitals are real for simplicity. Eqs. (1) and (2) imply that we can write the energy as an explicit functional of the canonical orbitals $ψi$ and their occupation numbers n_i,

where the occupation numbers are given by

It follows that

subject to the orthonormality constraint. Note here that we have replaced the set {n_i} by the only variable n₀ in the set. Apparently, the minimizing orbitals are dependent on n₀, let us denote them as ${ψi(n0)}$⁠. The fact that $ψi(n0)≠ψi(n)$ for $n0≠n$ is usually referred to as orbital relaxation. Now we prove the following lemma regarding the effect of orbital relaxation on the energy curvature.

Lemma 1 Let $En0(n)=E[{ψi(n0)};n]$ be the M+n system energy evaluated using the minimizing orbitals of M + n₀ system, and let $Er(n)=E[{ψi(n)};n]$ be the M + n system energy with its relaxed orbitals. Assuming $En0(n)$ and E_r(n) are smooth functions of n with second derivatives, then

Proof 1By definition of the second derivative,

and similarly,

where F₀and F_rare the terms in the curly brackets of the second equalities of Eq.(7)and Eq.(8), respectively. Since$E[{ψi(n0)};n0±m]⩾E[{ψi(n0±m)};n0±m]$due to the orbital relaxation, we have

When taking the limit in Eqs.(7)and(8), we obtain the inequality of(6).

As a remark, the inequality of (6) suggests that the energy curve at fractional systems becomes more concave (curvature more negative) when the orbital relaxation effect is taken into account for any continuous functional in the KS or GKS schemes.

Another comment is that the effect of frozen orbitals on the first derivative of energy can be assessed similarly. However, instead of an inequality relation between $∂En0∂n$ and $∂Er∂n$⁠, we end up with an equality. This is because

Since ${ψi(𝒓)}$ are the minimizing orbitals of $E[ρs(n0)]$⁠, we have $δEδψi(𝒓)|n=n0=0$⁠, thus

This means that the orbital relaxation has no effect on the first derivative of energy, as shown previously.¹⁵ 

Now let us focus on the second derivatives computed with frozen orbitals. As a special case, we evaluate $∂2En0∂n2|n=n0$ for the HF functional. Since the linear terms in the HF functional have no contributions to the second derivatives, it suffices for us to evaluate only for the nonlinear terms, namely, the Coulomb and exchange functionals, whose expressions are given by the following:

and

where $ρs(𝒓,𝒓′)=∑iMψi(𝒓)ψi(𝒓′)+nψM+1(𝒓)ψM+1(𝒓′)$ and $ρ(𝒓)=ρs(𝒓,𝒓)$⁠. Under the frozen orbital assumption, it follows

and

Then it is obvious that

Therefore,

By Lemma 1, we thus have

for all $0<n0<1$⁠. The equality holds if and only if orbital relaxation has no effect. Eq. (18) implies that the HF E vs N curve is piecewise concave.

Here we remark that it has been recognized through numerical observations in the previous literature that orbital relaxation is responsible for making the HF E vs N curve concave.^7,17,18 Without orbital relaxation, it will display a linear curve. Now this is justified in our proof.

Semi-local density functionals differ from the HF functional in the exchange correlation part. Without considering orbital relaxation, the energy curvature at fractional systems is given by a similar expression as in Eq. (17), with the −K term replaced by $Exc$⁠. In particular,

where $f(𝒓)=|ψ¯M+1(𝒓)|2$ is the frontier orbital density from the self consistent field (SCF) calculation of the M + n₀ system. We note Eq. (19) agrees with the formula in Ref. 19 without considering the orbital relaxation, i.e., the energy curvature is computed using Eq. (7). Under this frozen orbital approximation, the concavity or convexity of the E vs N curve is dictated by the sign of the energy curvature in Eq. (19), which results from the competition between the positive Coulomb curvature (first term) and the negative exchange correlation curvature (second term).

In the following, we use exchange-only local density approximation (XLDA) as an example to demonstrate that the Coulomb curvature is not always the winner of the competition and hence $∂2En0∂n2|n=n0$ is not guaranteed to be non-negative. As a consequence, the E vs N curves for semi-local density functionals are not guaranteed to be convex under the frozen orbital assumption at least.

For XLDA, $Exc=−Cx∫[ρ(𝒓)]4/3d𝒓$⁠, where $Cx=32(34π)1/3$⁠, and Eq. (19) reads

If N = 0, i.e., the system contains less than one electron, $ρ(𝒓)=nf(𝒓)$⁠. It follows that

Note that the right hand side of Eq. (21) is not bounded from below and thus for small enough n₀, the energy curvature becomes negative.

When the orbital relaxation effect is taken into account, the above conclusion still holds. This is demonstrated by the SCF calculation of the E vs N curve and computing the numerical curvature for each point on the curve for various functionals through Eq. (8), where the limit definition is replaced by a finite difference calculation with m = 0.001. The SCF energy curvature is then plotted as a function of the number of electrons of H atom (we restrict $N∈[0,2]$⁠), as shown in Figure 1. Note that for both $N∈[0,1]$ and $N∈[1,2]$⁠, the number of electrons for the relevant spin is no greater than 1. For XLDA, as can be seen, for either $N∈[0,1]$ or $N∈[1,2]$⁠, there exists a turning point N₀, around 0.1 and 1.1, respectively, such that the energy curvatures are negative on its left side and become positive on its right side. Moreover, the negative curvature region (⁠$n∈[0,N0)$⁠) corresponds to adding a small fraction of electrons to the integer system. This is often overlooked in the literature, since the negative curvature in the small fraction region is not visible from the E(N) plot.

Such phenomena are observed for other more complicated functionals as well. For LDA implementing VWN5²⁰ as the correlation functional, when $0<N<1$⁠, i.e., the number of electrons is less than 1, the curvatures essentially agree with XLDA because of the fact that the correlation effect from the same spin has negligible contribution to the curvature; when $1<N<2$⁠, the LDA curvature becomes more positive than XLDA due to the correlation from the opposite spins, nevertheless there remains a negative curvature region near N = 1. For the generalized gradient approximation (GGA) functionals, such as BLYP,^21,22 the curvatures display a modest deviation from LDA, yet maintaining the same trend. The HF curvatures have been shown to be non-positive. This is verified in the curvature plot. When $0<N<1$⁠, the curvature is zero because HF is exact for systems with no more than one electron and hence displays a correct linear E vs N curve with zero curvature; when $1<N<2$⁠, HF is no longer exact and the curvatures become negative due to the orbital relaxation. Hybrid functionals, such as B3LYP,^22,23 have both semi-local and HF functional components, thus the curvature is expected to take values in between them. This is also manifested in Figure 1, where the B3LYP curvatures follow the trends of LDA and BLYP but are slightly more negative.

The analysis above is valid only for one-electron cases. For many-electron systems $(N⩾1)$⁠, however, we can modify our analysis as follows. To characterize the energy curvature in the right neighborhood of the integer point, i.e., $n→0+$⁠, let $ρ(𝒓)=ρ0(𝒓)+nf(𝒓)=ρ1(𝒓)−(1−n)f(𝒓)$⁠, where $ρ0(𝒓)=∑i=1M|ψi(𝒓)|2$⁠, $ρ1(𝒓)=∑i=1M+1|ψi(𝒓)|2$⁠, and $f(𝒓)=|ψM+1(𝒓)|2$⁠. In the frozen orbital approximation, $ψi(𝒓)(i=1,…,M)$ are the occupied orbitals of the integer N-system, while $ψM+1(𝒓)$ is its LUMO. It is obvious that

Let

and let

It follows that

and

Let us assume $A>B1$⁠, otherwise the E vs N curve in the right neighborhood of the integer-N system is completely concave. Let $ηLUMO=AB0$⁠, where the superscript LUMO highlights that the frontier orbital used in the evaluation is the LUMO of the integer-N system. It follows that if $ηLUMO⩾1$⁠, then the energy curvature is always positive and resulting in a convex E_v(N) curve. If $ηLUMO<1$⁠, since

there exists a turning point n₀ such that when $0<n<n0$⁠, the energy curve is concave due to the negative curvature; and when $n0<n<1$⁠, the curve becomes convex as a result of the positive curvature.

Similarly, to characterize the curvature in the left neighborhood of the integer-N system, we define $ηHOMO=AB0$⁠, where the corresponding A and B₀ are evaluated using the HOMO of the integer system. This extension is straightforward by considering the $n→0−$ limit and replacing LUMO by HOMO of the reference integer system. Details are omitted. Note that $ηLUMO$ and $ηHOMO$ are intrinsic properties of the reference integer-N system, and their values relative to 1 serve as a criterion for whether or not there is a concave region in the right (left) neighborhood of the integer point.

In Table I, we list $ηLUMO$ for a selected set of atoms and molecules (at their optimized structures) using XLDA calculations. Here we add a subscript “+” or “0” to distinguish the reference cation or neutral system. Noticing that for real systems, $ηLUMO$ is spin dependent, we also present the spin channel for the LUMO. Based on the numbers in Table I, it is obvious that $ηLUMO$ is not guaranteed to be greater than 1, and thus the E vs N curve could have concave regions on the right side of an integer point.

In Table II, we further present $ηHOMO$ for these systems with XLDA orbitals. Here the reference systems we use are neutral or anion systems, to characterize E_v(N) on the same [M, M + 1] interval as the LUMO analysis but from the $n→1−$ side. As can be seen, $ηHOMO$ are all greater than 1, suggesting that E_v(N) in the left neighborhood of the integer systems are completely convex. Moreover, by comparison between Tables I and II, we note that $ηLUMO$ are often smaller than $ηHOMO$⁠, suggesting that the curve of $n→0+$ is more likely to have concave regions than $n→1−$⁠. This is reasonable because the LUMO density is likely to spread more diffusively than the HOMO density, which gives rise to a larger B₀ through the integration of $[ρ0(𝒓)]−2/3[f(𝒓)]2$⁠, and hence a smaller η.

To assess the orbital relaxation effect on the energy curvatures, in Figure 3, we show the energy curvature as a function of the fractional number of electrons of the H₂O molecule obtained by SCF calculations. Here we focus on the small fraction regions, as illustrated in Figure 2. The curvatures are computed via Eq. (8) using the finite difference evaluation of the SCF energies with an interval m = 0.001. As can be seen, in Figs. 3(a), 3(b), and 3(d), the energy curvatures for XLDA are positive throughout, while in Fig. 3(c), there is a negative region in the right neighborhood of n = 0. This suggests that the piece of the E vs N curve going from the cation to the neutral molecule is completely convex, but going from the neutral to the anion will encounter a concave region. This agrees with the η analysis in Tables I and II neglecting the orbital relaxation effect, which validates the accuracy and usefulness of that simple judgment. Moreover, in Figure 3, the other semi-local and hybrid functionals follow the similar trend as XLDA, which implies that the correlation parts in these functionals indeed have minor effects. In contrast to the behaviors of these functionals, HF curvatures are non-positive throughout, as expected.

To further clarify the shape of the E vs N curve of the H₂O molecule for small fractions in the highlighted regions of Figure 2, in Figure 4 we present the relative LDA energies when adding or removing the small fractional electron (⁠$n→0+$ and $n→1−$⁠) through self consistent calculations. To pronounce the convexity or concavity behavior in such regions, instead of showing the absolute energy, we subtract a linear term and only plot the relative energy with respect to (i) the tangent line that passes through the energy point at n = 0 (or n = 1) and follows the right (or left) slope of its energy derivative (which is equal to the LUMO or HOMO energy) and in (ii) thelinear interpolation between the neighboring integer LDA energies. Here in (i) we focus on the concave/convex nature of the energy curve; while in (ii) we highlight the relationship between E_v(N) and the straight line segment connecting the integer energies. To enlarge the regions of interest, we further scale the relative energy by 10⁵ for (i) and 10³ for (ii) to visualize them in the same figure. In Figs. 4(a) and 4(b), the curves characterize the process of H₂O⁺→ H₂O; while in Figs. 4(c) and 4(d), they depict the process of H₂O → H₂O⁻. Since the relative LDA energies differ from the true LDA energies by a linear function of the fractional number n, the curvatures maintain the same up to the scaling factor. As manifested by the blue curves which highlight the curvatures, it is obvious that they are convex in Figs. 4(a), 4(b), and 4(d), while there is a transition from concavity to convexity in Fig. 4(c). Moreover, the turning point in Fig. 4(c) occurs around n = 0.01, which agrees with Figure 3(c), where at $n≈0.01$ the LDA curvature intersects with 0.

The deviation of LDA energy from the straight line interpolation between the integer energy points are shown in red in Figure 4. As manifested in all the subplots, the deviations are almost straight lines for small n and negative in sign, which suggests that the LDA energy derivative on the right (⁠$ϵLUMO$⁠) is underestimated, while the energy derivative on the left (⁠$ϵHOMO$⁠) is overestimated. This is true regardless of whether the E vs N curve starts with concavity or convexity.

It is worth noticing that the orbital energy is equal to the first derivative of the total energy with respect to the electron number, while the energy curvature is a quantity that relates to the second order energy derivative. Therefore they are independent properties of the function E(N). In terms of the relation of the E vs N curve in comparison with the reference straight line interpolation between integers, the formerdetermines whether the curve lies above or below the reference line for small fractions, while the latter determines whether the line will curve concavely or convexly. Furthermore, we highlight that the E vs N curve on [M, M+1] being convex is a sufficient but not necessary condition for the underestimation of the LUMO energy at N = M. In fact, we have found examples with non-convex E vs N curves, such as H₂O as demonstrated in Figures 3 and 4, and the atoms and molecules suggested in Table I with $ηLUMO<1$⁠, yet we have not found a counterexample where the LUMO energy, the chemical potential evaluated as an energy derivative with respect to the electron number, overestimates the minus electron affinity obtained by taking the difference of the integer energies. Similar argument applies for the HOMO energy. The conclusion that the energy gap is underestimated by semi-local functionals can still be valid even if the E vs N curve is not convex everywhere. This is because the fractional $Ev$(N) curve from semi-local functionals lies below the reference line is a weaker statement than the claim that $Ev$(N) is convex. And it is the former statement that better describes the situation.

So far our discussion has been focused on HF and semi-local functionals including global hybrids. Here as a remark, for other more sophisticated functionals, such as the random phase approximation (RPA), the self-consistent implementation is available but is not usually used. In the post-SCF implementation, it has been known that the particle-hole RPA (ph-RPA) suffers from larger delocalization error than LDA,²⁴ which implies a more positive curvature in their $Ev$(N) curve. The recently proposed particle-particle RPA (pp-RPA) eliminates most of the delocalization error in the post-SCF manner,²⁵ but their curvatures have not been investigated. This could be an interesting topic for future exploration.

In the literature, it has been inferred from the concavity or convexity of the E vs N curve of approximate density functionals that the former class of functionals suffers from the localization error while the latter suffers from the delocalization error.⁵ We note that this conclusion still holds, although whether a specific functional behaves convexly or concavely in the E vs N curve needs careful justification. In particular, indeed the HF functional suffers from the localization error due to its concavity throughout the E vs N curve. This means that for any line segment connecting a point on $Ev$(N) restricted to $N∈(M,M+1)$ and the integer point $Ev$(M) or $Ev$(M + 1), it should lie below the $Ev$(N) curve. As a consequence, suppose we add an electron to a system composed of multiple infinitely separated identical subsystems, HF predicts that the added electron density will localize to only one subsystem (suppose the curvature is not constantly zero). In contrast, if a semi-local functional displays a thoroughly convex E vs N curve, any line segments connecting an integer point and a fractional point should lie above the $Ev$(N) curve, and an added or removed electron density will delocalize to as many subsystems as possible. However, as shown in this paper, the E vs N curve from semi-local density functionals could display a complicated scenario, where there is a small but concave region at the beginning of the curve, corresponding to adding a small fraction of electrons to an integer system. Equivalently, not all line segments lie above the E vs N curve. In such cases, the electron density will neither localize to one nor delocalize to infinite number of subsystems, instead, it will delocalize to finite number of subsystems, and the extent of delocalization is determined by the energetically most favorable pattern. Note that functionals satisfying the linearity condition should not bias toward localization or delocalization, because they are degenerate in energy.

For example, in the case of adding one electron to M infinitely separated protons, this additional electron density will distribute onto 6 protons when $M>6$⁠. This is illustrated in Figure 5, where we plot the LDA energy of M infinitely separated protons with one electron, as a function of M. As can be seen, this function is not monotone decreasing with M (note that if it was, the LDA $Ev$(N) curve would be completely convex, proof is straightforward and omitted here) but has a minimum at M = 6 due to the presence of a concave region in the $Ev$(N) curve of a single H atom. Therefore, when $M⩽6$⁠, the added electron tends to delocalize its density to as many protons as possible (i.e., to M protons) in order to decrease its energy. However, when $M>6$⁠, further delocalization no longer displays energetic advantage. Therefore, the added electron will only delocalize its density over 6 of the protons, leaving the rest of the protons having nothing. Here we have assumed that the number of electrons near each proton can either take a fixed fractional value or zero, in order to reduce the degrees of freedom and simplify the problem. This has been verified by numerical calculations, where we set 15 protons on a circle (equally spaced by 1 m = 10¹⁰Å, where the large distance is to diminish the impact of Coulomb interactions between protons). When adding one electron to the system, by implementing gradient descent minimization, the electron would indeed evenly delocalize its density to 6 of the protons, leaving 9 bare protons, and the energy exactly matches the prediction in Figure 5. Note that the ways of distributing the electron is not unique, and there exist multiple local minimizers as well as degenerate global minimizers; what we find is one of the degenerate global minimizers.

In terms of the frontier orbital energies, the underestimation of the ionization energy I by $εHOMO$ and the overestimation of the electron affinity A by $εLUMO$ by semi-local functionals are direct consequences of the wrong first order derivative of the energy, although they are indirectly related to the energy curvature. The previous conclusions on the underestimation of ionization potential and overestimation of electron affinity, remain unchanged; the convexity of the E_$v$(N) curve for all fractions is not entirely necessary. Thus, instead of using the term convexity, a more rigorous description of the delocalization error in an approximate density functional is that its E_v(N) curve in the fractional region lies below the straight line segment across the neighboring integer points.

Support from the National Science Foundation (Grant No. CHE-13-62927) (C.L.) is gratefully appreciated. W.Y. appreciates the support as part of the Center for the Computational Design of Functional Layered Materials, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences under Award No. DE-SC0012575.