1 Introduction

Density functional theory (DFT) is, at present, the method of choice for the description of the ground-state (GS) properties of quantum systems with relatively low computational complexity. A practical and effective approach for electronic structure calculations is offered by Kohn–Sham (KS) method. The great success of KS-DFT is due to the fact that simple density functional approximations (DFAs) perform remarkably well for a wide range of problems in chemistry and physics, particularly for the prediction of the structure and thermodynamic properties of molecules and solids. DFT owes its power to an efficient formalism founded on the single-determinant state function that should in principle be able to yield the exact energy and density. All the complexity of this theory is hidden in one term––the exchange–correlation functional. This term holds ‘the key to the success or failure of DFT’ [1]. The main problem is that the explicit form of the exact functional is unknown, and only certain desired properties which should be satisfied by approximate functionals can be formulated [2, 3]. Failures in practical calculations are not the breakdown of the theory itself, but are only due to deficiencies in commonly used approximations [1]. Common approximations for the exchange–correlation functional have been found to give big errors for the linearity condition of fractional charges (electron numbers), leading to delocalization error, and for the constancy condition of fractional spins (the energy of the ensemble constructed from the components of the spin doublet is independent of the spin number), leading to static correlation error [1, 4, 5]. Violation of these conditions underlies the major failures of all known approximate functionals to describe the energy gaps and related properties in strongly correlated systems [6]. The linearity condition of the energy versus the number of electrons N follows from the zero-temperature limit of the temperature extension of the DFT [7–9] or, alternatively, from the pure states for a collection of replicas at zero temperature [2, 10]. Standard local or semi-local DFAs such as local density approximation (LDA), generalized gradient approximation (GGA) and meta-GGA, result in a very poor description, giving a smooth, convex, almost parabolic, interpolation between the values at integers. Some corrections applied to these functionals straighten the curve for fractional charges but significantly worsen the description at integers. Hartree–Fock (HF) results have the opposite behavior, with a concave interpolation between the integers. However, the hybrid functionals which combine these two ingredients are still quite convex [11, 12]. There have been two strategies to tackle this problem. The first one is the scaling correction method [13–15] or the restoration scheme [16, 17] applied to existing functionals. The second is the construction of new functionals which satisfy the energy linearity condition [11, 18, 19] or the constancy condition of the frontier orbital energy (interpreted as the ionization energy) [20–22]. The “dilemma” of how to define the exchange–correlation potential for fractional N was discussed in early paper [23], while the papers [24, 25] represent the most recent investigations related to the fractional-N problem for density functionals. Recently, an orbital functional based on the particle–particle random phase approximation to many-body theory has been shown to exhibit minimal errors for fractional charges and fractional spins and have many promising features [26, 27].

In this article, we investigate four properties related to fractional number of electrons to test the quality of the exchange–correlation functional. As illustrative example, the results for DFT method with five exchange–correlation functionals (SVWN5, PBE, B3LYP, CAM-B3LYP, rCAM-B3LYP) and the Hartree–Fock method are presented.

2 Theoretical background

In the exact DFT [28], the total energy of the system is a functional of the electron density and the external potential v(r),

$$\begin{gathered} E\left[ {\rho ,v} \right] = T\left[ \rho \right] + E_{\text{ee}} \left[ \rho \right] + \int {v\left( {\mathbf{r}} \right)\rho \left( {\mathbf{r}} \right)d{\mathbf{r}}} \\ = T_{s} \left[ \rho \right] + E_{\text{es}} \left[ \rho \right] + E_{\text{xc}} \left[ \rho \right] + \int {v\left( {\mathbf{r}} \right)\rho \left( {\mathbf{r}} \right)d{\mathbf{r}}} \\ \end{gathered}$$
(1)

where T[ρ] and E ee[ρ] are the kinetic energy and the electron–electron interaction energy functionals, respectively. Following the Kohn–Sham (KS) scheme, the sum of these two functionals can be represented as the sum of the non-interacting reference KS kinetic energy, T s [ρ], the electrostatic energy, E es[ρ], and the exchange–correlation energy, E xc[ρ]. The GS energy, E[N, v], is the minimum E[N, v] = Min ρ → N E[ρ, v] with the constraint \(N = \smallint \rho \left( {\mathbf{r}} \right)d{\mathbf{r}}.\) The minimizer is the GS density. The main result of the fractional-number extension of the DFT at zero temperature is that the GS energy is continuous and piecewise linear function of the electron number (charge) [8]. Namely, for the system with fractional charge N = J ± δ, where J is an integer and δ ∊ [0, 1], the GS energy is

$$E_{ \pm \delta }^{{}} \left[ {J,v} \right] \equiv E\left[ {J \pm \delta ,v} \right] = \delta E\left[ {J \pm 1,v} \right] + \left( {1 - \delta } \right)E\left[ {J,v} \right],$$
(2)

and the GS density is a linear mixture of the GS densities of the J- and the (J ± 1)-electron number systems,

$$\rho_{ \pm \delta }^{{}} \left( {{\mathbf{r}};J,v} \right) \equiv \rho \left( {{\mathbf{r}};J \pm \delta ,v} \right) = \delta \,\rho \left( {{\mathbf{r}};J \pm 1,v} \right) + \left( {1 - \delta } \right)\rho \left( {{\mathbf{r}};J,v} \right).$$
(3)

The total energy is an implicit functional of the density via the dependence of the molecular orbitals and their occupancies on the density. The term representing the coupling between the electrons and the external potential is an exact, explicit functional of the electron density. The kinetic energy term, expressed as an explicit functional of spinorbitals and occupation numbers, is

$$\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{T}_{s} \left[ {\left\{ {n_{p} } \right\},\left\{ {\chi_{p} } \right\}} \right] = \sum\limits_{p} {n_{p} \int {\chi_{p}^{*} \left( {\mathbf{r}} \right)\left( { - \frac{1}{2}\nabla_{{\mathbf{r}}}^{2} } \right)\chi_{p}^{{}} \left( {\mathbf{r}} \right)d{\mathbf{r}}} } .$$
(4)

E es represents the classical electrostatic energy

$$E_{\text{es}} \left[ \rho \right] = \frac{1}{2}\int {\int {\frac{{\rho \left( {\mathbf{r}} \right)\rho \left( {{\mathbf{r^{\prime}}}} \right)}}{{\left| {{\mathbf{r}} - {\mathbf{r^{\prime}}}} \right|}}d{\mathbf{r}}d{\mathbf{r^{\prime}}}} } .$$
(5)

Since such E es[ρ] is not linear in ρ(r), it immediately follows that its nonlinearity has to be compensated by the appropriate exchange–correlation term E xc[ρ], or the standard energy expression corrected by some additional terms [29, 30]. Alternatively, both E es and E x should be expressed in the ensemble terms [31]. In general, if the electrostatic energy for the (J ± δ) system is defined as a linear combination of the electrostatic energies for the ‘end’ points, the form of the electrostatic energy as a functional of density is unknown [32]. Although no explicit form is available for the exact E xc[ρ], much is known about the “proper” way in which approximations for this energy term should be constructed. One of the main challenges for DFT is to keep, as its cornerstone, some element of simplicity [4]. Unfortunately, a simple density functional adopted from the solid-state physics, the LDA, does not perform well in many areas of chemistry. The introduction of the density gradient into the form of the exchange–correlation functional (the GGA or meta-GGA), allows to receive satisfactory results in the chemical applications. The next major advances came with the inclusion of a fraction of the HF exchange in the E xc functional (hybrid functionals). One of the recent developments in functionals is due to the range separation. The idea is to separate the electron–electron interaction into two parts, the long-range one and the short-range one, and then to treat these parts with different functionals (see Ref. [4] for review).

In practice, the direct extension (to fractional N) of DFT is applied, for which the ground-state energy is evaluated as

$$E^{\text{DFA}} \left[ {N,v} \right] \equiv \mathop {\text{Min}}\limits_{{\left\{ {n_{p} } \right\},\left\{ {\chi_{p} } \right\}}} \left\{ {\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{E} \left[ {\left\{ {n_{p} } \right\},\left\{ {\chi_{p} } \right\};v} \right]\left| {\sum\limits_{p} {n_{p} } = N,n_{p} \in \left[ {0,1} \right],\left\langle {{\chi_{p} }} \mathrel{\left | {\vphantom {{\chi_{p} } {\chi_{q} }}} \right. \kern-0pt} {{\chi_{q} }} \right\rangle = \delta_{pq} } \right.} \right\}$$
(6)

where \(\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{E} \left[ {\left\{ {n_{p} } \right\},\left\{ {\chi_{p} } \right\}} \right]\) is the energy functional defined in terms of commonly used approximations

$$\begin{gathered} \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{E} \left[ {\left\{ {n_{p} } \right\},\left\{ {\chi_{p} } \right\},v} \right] = \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{T}_{s} \left[ {\left\{ {n_{p} } \right\},\left\{ {\chi_{p} } \right\}} \right] + \int {v\left( {\mathbf{r}} \right)\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\rho } \left( {{\mathbf{r}};\left\{ {n_{p} } \right\},\left\{ {\chi_{p} } \right\}} \right)d{\mathbf{r}}} \\ \; + E_{\text{es}} \left[ {\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\rho } \left[ {\left\{ {n_{p} } \right\},\left\{ {\chi_{p} } \right\}} \right]} \right] + b_{\text{xc}} E_{\text{xc}} \left[ {\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\rho } \left[ {\left\{ {n_{p} } \right\},\left\{ {\chi_{p} } \right\}} \right]} \right] \\ \; - \frac{1}{2}\sum\limits_{pq} {n_{p} n_{q} \left\{ {c_{\text{x}} \left\langle {{\chi_{p} \chi_{q} }} \mathrel{\left | {\vphantom {{\chi_{p} \chi_{q} } {\chi_{q} \chi_{p} }}} \right. \kern-0pt} {{\chi_{q} \chi_{p} }} \right\rangle + b_{\text{x}}^{\text{LR,HF}} \left\langle {{\chi_{p} \chi_{q} }} \mathrel{\left | {\vphantom {{\chi_{p} \chi_{q} } {\chi_{q} \chi_{p} }}} \right. \kern-0pt} {{\chi_{q} \chi_{p} }} \right\rangle^{\mu } } \right\}} \\ \end{gathered}$$
(7)

with a set of orthonormal spinorbitals, {χ p } and their occupation numbers, {n p }. Here, the non-integer electron number N(fractional charge) is realized indirectly by the constraint on the spinorbital occupancies, {n p }. The density of the system (the non-interacting KS density) is

$$\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\rho } \left( {{\mathbf{r}};\left\{ {n_{p} } \right\},\left\{ {\chi_{p} } \right\},v} \right) = \sum\limits_{p} {n_{p} \left| {\chi_{p}^{{}} \left( {\mathbf{r}} \right)} \right|^{2} } ,$$
(8)

and can integrate to any nonnegative real number N = ∑ p n p . The GS density ρ(r; N, v) is obtained from Eq. (8) with {χ p } and {n p } being the minimizers in Eq. (6). \(E_{\text{xc}} [ {\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\rho } } ]\) represents the part of the exchange–correlation energy which depends locally on the density (and on the density derivatives). If b xc = 0, there is no correlation part and all exchange is represented by the HF-like exchange. The parameters c x and \(b_{\text{x}}^{{{\text{LR}},{\text{HF}}}}\) determine the portion of the HF exchange in the hybrid and the long-range correction functional, respectively. The parameter μ defines separation of regions. The spinorbitals {χ p } and orbital energies {ɛ p } are the eigensolutions of the one-electron Hamiltonian with the effective potential v eff(x). The v eff(x) is a local (multiplicative) potential in the original KS method or a non-local potential in a generalized KS method (with hybrid functionals) and in the HF method.

As a consequence of the linearity condition, Eq. (2), the chemical potential (the first derivative of the GS energy with respect to N) is constant between the integers

$$\mu \left[ {J \pm \delta ,\nu } \right] \equiv \left. {\left( {\frac{{\partial E\left[ {N,\nu } \right]}}{\partial N}} \right)_{\nu } } \right|_{N = J \pm \delta } = \pm \left( {E\left[ {J \pm 1,\nu } \right] - E\left[ {J,\nu } \right]} \right),\;\;\delta \in \left( {0,1} \right) ,$$
(9)

and it is related to the frontier (f) eigenvalue [8, 12, 33–35]

$$\mu \left[ {J \pm \delta ,v} \right] = \varepsilon_{f}^{\text{DFA}} \left[ {J \pm \delta ,v} \right] .$$
(10)

Now, we can formulate the first two indicators for testing the quality of functionals: the relative deviation from the global linearity condition (GLC), Eq. (2), as

$$\Delta_{ \pm \delta }^{\text{GLC}} = \frac{{E_{ \pm \delta }^{\text{DFA}} }}{{\delta E_{ \pm 1}^{\text{DFA}} + \left( {1 - \delta } \right)E_{0}^{\text{DFA}} }} - 1,$$
(11)

and the relative deviation from the global constancy condition (GCC), Eqs. (9) and (10), as

$$\Delta_{ \pm \delta }^{\text{GCC}} \equiv \frac{{\varepsilon_{f}^{\text{DFA}} \left[ {J \pm \delta } \right]}}{{ \pm \left( {E_{ \pm 1}^{\text{DFA}} - E_{0}^{\text{DFA}} } \right)}} - 1.$$
(12)

Here, \(E_{ \pm \delta }^{\text{DFA}} = E^{\text{DFA}} \left[ {J \pm \delta , v} \right]\).

Besides the linearity condition for the total energy and the constancy condition for the energy of the highest (partially) occupied spinorbital, we can examine the local counterparts of these conditions related to the ground-state density and the Fukui function (FF) [28, 36]. Third indicator tests the fulfillment of the local linearity condition (LLC), Eq. (3). It is the integrated deviation from the LLC, defined as

$$\Delta_{ \pm \delta }^{\text{LLC}} = \left( {J \pm \delta } \right)^{ - 1} \int {\left| {\rho_{ \pm \delta }^{\text{DFA}} \left( {\mathbf{r}} \right) - \left( {\delta \rho_{ \pm 1}^{\text{DFA}} \left( {\mathbf{r}} \right) + \left( {1 - \delta } \right)\rho_{0}^{\text{DFA}} \left( {\mathbf{r}} \right)} \right)} \right|d{\mathbf{r}}} ,\,\,\,\,\rho_{ \pm \delta }^{\text{DFA}} \left( {\mathbf{r}} \right) = \rho \left( {{\mathbf{r}},J \pm \delta ,v} \right) .$$
(13)

(The factor (J ± δ)−1 is omitted in our illustrative example tests as it is irrelevant at comparisons).

The FF (the derivative of the ground-state density, Eq. (3), with respect to the electron number) is independent of the fraction δ (similarly as the chemical potential)

$$f\left( {{\mathbf{r}};J \pm \delta ,v} \right) \equiv \left. {\left( {\frac{{\partial \rho \left( {{\mathbf{r}};N,v} \right)}}{\partial N}} \right)_{v} } \right|_{N = J \pm \delta } = \pm \left( {\rho \left( {{\mathbf{r}};J \pm 1,v} \right) - \rho \left( {{\mathbf{r}};J,v} \right)} \right)$$
(14)

(this form of FF was noted first in Ref. [37]), and the integrated deviation from this local constancy condition (LCC) is defined as

$$\Delta_{ \pm \delta }^{\text{LCC}} = C^{\text{LCC}} \int {\left| {f_{ \pm \delta }^{\text{DFA}} \left( {\mathbf{r}} \right) \mp \left( {\rho_{ \pm 1}^{\text{DFA}} \left( {\mathbf{r}} \right) - \rho_{0}^{\text{DFA}} \left( {\mathbf{r}} \right)} \right)} \right|d{\mathbf{r}}} ,\;C^{\text{LCC}} = \left( {\int {\left| {\left( {\rho_{ \pm 1}^{\text{DFA}} \left( {\mathbf{r}} \right) - \rho_{0}^{\text{DFA}} \left( {\mathbf{r}} \right)} \right)} \right|d{\mathbf{r}}} } \right)^{ - 1}$$
(15)

where \(f_{ \pm \delta }^{\text{DFA}} \left( {\mathbf{r}} \right) \equiv \left. {\left( { \partial \rho_{ \pm \delta }^{\text{DFA}} \left( {\mathbf{r}} \right)/ \partial \rho_{ \pm \delta }^{\text{DFA}} \left( {\mathbf{r}} \right) \partial N. \partial N} \right)} \right|_{N = J \pm \delta }\). (The factor C LCC is omitted in our illustrative example tests as it is irrelevant at comparisons). Note that (C LCC)−1 equals to one when the FF is not negative, but it is larger than one (even two) when FF is negative in some regions as in cases of redox-induced electron transfer [38, 39]. This ΔLCC is the last indicator for testing the quality of functionals, proposed in this work. It should be noted that if the integrand in Eqs. (13) and (15) is taken without absolute value operation, the result of integration is exactly zero. According to the proposed definitions of indicators, high quality of a particular DFA is demonstrated when all inequalities

$$\left| {\Delta_{ \pm \delta }^{GLC} } \right| \ll 1, \, \left| {\Delta_{ \pm \delta }^{GCC} } \right| \ll 1,\Delta_{ \pm \delta }^{LGC} \ll 1,\Delta_{ \pm \delta }^{LCC} \ll 1,$$
(16)

are satisfied in the whole range 0 < δ < 1.

3 Illustrative examples and discussion

In this article, the performance of SVWN5 [40] (a LDA-type functional), PBE [41] (a GGA-type functional), B3LYP [42, 43] (a hybrid functional), CAM-B3LYP [44] (a long-range corrected functional), rCAM-B3LYP [11] (a long-range corrected functional with improved description of systems with fractional numbers of electrons), and HF method was tested. To eliminate complication due to the spatial degeneracy, our test set is restricted only to atomic systems with the total angular momentum equal zero, namely to He, Li, Be, Na, Mg atoms. All calculations were performed with cc-pVTZ basis set [45–47] in spin-unrestricted formalism. The FF for the fractional electron system [48], f DFA±δ (r), was calculated using QM4D program [49] according to methods described in Ref. [48]. Only decreasing electron number process was considered for all systems (lower sign chosen in Eqs. (2), (3), (9)–(15), δ = 1 means a cationic form of an atom). All results are presented per thousand. The results for atoms are accessible from Supporting Information.

To present our illustrative example results in a condensed form, the average value of an indicator over a set of four atoms, Li, Na, Be, Mg is displayed for each deviation. The results for He are not included, as significantly differing from those of other atoms in many cases: for ΔGLC at all methods (see Table 1), for remaining indicators at rCAM-B3LYP. This peculiarity stems from the fact that the two-particle aspects of electron structure in comparison with the single-particle ones are more pronounced for the He atom than for remaining atoms. In all Figures, lines joining calculated points are drawn for eye guidance only.

Table 1 The average (over the set: Li, Na, Be, Mg) of the test indicator for the global linear condition Δ GLC−δ
Full size table

The results of the global linearity condition testing, Eq. (11), are collected in Table 1 (zero values at δ = 0 and δ = 1 are omitted) and in Fig. 1. As expected, the standard local-potential functionals (SVWN5 and PBE) yield the largest, positive \(\Delta_{{_{ - \delta } }}^{{^{\text{GLC}} }}\). The HF method yields the weakest, negative \(\Delta_{{_{ - \delta } }}^{{^{\text{GLC}} }}\). The B3LYP functional combines these two effects, yielding \(\Delta_{{_{ - \delta } }}^{{^{\text{GLC}} }}\) slightly smaller than those generated by SVWN and PBE functionals, but still presents a concave curve. A significant improvement is observed for functionals which contain some Coulomb-attenuated exchange (CAM). The rCAM-B3LYP functional yields the smallest GLC deviation. Reasonable explanation of observed tendencies is available, especially of the fact that extreme deviations are observed at δ = 0.5 for all methods. Very weak GLC deviation for the HF method is the effect of the testing set (only atomic systems without spatial degeneracy) used in this article. All the occurring open shell systems at integral N are not spatially degenerate and the HF method for such systems is self-interaction free [31].

Fig. 1
figure 1

The average of \(\Delta_{ - \delta }^{\text{GLC}}\), Eq. (11), vs. fraction δ

Full size image

At first, it should be noted that the energy expression, Eq. (7), has no explicit N dependence. This dependence is implicit trough the dependence on {n p } and the condition ∑ p n p  = N. Namely, the trace of the first-order density matrix (DM), d 1, is equal to electron number N

$$d_{1} \left( {{\mathbf{x}};{\mathbf{x}}{\prime };\left\{ {n_{p} } \right\},\left\{ {\chi_{p} } \right\}} \right) \equiv \sum\limits_{p} {n_{p} \chi_{p}^{*} \left( {\mathbf{x}} \right)\chi_{p} \left( {{\mathbf{x}}{\prime }} \right)} ,\quad {\text{Tr}}d_{1} = N .$$
(17)

In the direct extension of the HF method, the hyper-HF one, the approximate second-order DM, d 2, is constructed from the first-order DM as the following determinant

$$d_{2} \left( {{\mathbf{x}}_{1} ,{\mathbf{x}}_{2} ;{\mathbf{x}}{\prime }_{1} ,{\mathbf{x}}{\prime }_{2} ;d_{1} } \right) \equiv \frac{1}{2}\left( {d_{1} \left( {{\mathbf{x}}_{1} ,{\mathbf{x}}{\prime }_{1} } \right)d_{1} \left( {{\mathbf{x}}_{2} ,{\mathbf{x}}{\prime }_{2} } \right) - d_{1} \left( {{\mathbf{x}}_{1} ,{\mathbf{x}}{\prime }_{2} } \right)d_{1} \left( {{\mathbf{x}}_{2} ,{\mathbf{x}}{\prime }_{1} } \right)} \right) = d_{2} \left[ {N,v} \right] ,$$
(18)

(it is exact for the single Slater determinant). This d 2 leads to relations [50] \(\text{Tr} d_{2} \left[ {N,v} \right] \ge \frac{1}{2}\left( {N - 1} \right)N\), and [31]

$$\Delta E_{{{\text{ee}}, \pm \delta }}^{\text{HF}} \equiv E_{{{\text{ee}}, \pm \delta }}^{\text{HF}} - \left( {\delta E_{{{\text{ee}}, \pm 1}}^{\text{HF}} + \left( {1 - \delta } \right)E_{{{\text{ee}},0}}^{\text{HF}} } \right) \ge 0 .$$
(19)

where \(E_{{{\text{ee}}, \pm \delta }}^{\text{HF}} = \text{Tr} d_{2} \left[ {J \pm \delta ,v} \right]\hat{V}_{ee}\)(the electron–electron energy due to ‘generalized Fock’ operator).

In the case of atoms considered in this work, the above deviation from the linearity for the electron–electron interaction energy is [31]

$$\Delta E_{{{\text{ee}}, \pm \delta }}^{\text{HF}} = \delta \left( {1 - \delta } \right)J_{ff} \ge 0 ,$$
(20)

(here J ff is the Coulomb integral for the frontier orbital) with the maximum at δ = 0.5. The shape of this derivation, Eq. (21), was already observed in Fig. 1 of Ref. [51]. When combined with the negative denominator in Eq. (11), it leads to negative \(\Delta_{{_{ - \delta } }}^{{^{\text{GLC}} }}\) for the HF method. Contrary to the hyper-HF method, compensation of the self-interaction errors of E es and E x for DFA methods is not complete or absent. Therefore, the deviation from linearity for E es, Eq. (5), is the dominating deviation. It can be easily evaluated as

$$\begin{gathered} \Delta E_{{{\text{es}}, \pm \delta }}^{\text{DFA}} \equiv E_{\text{es}} \left[ {\rho_{ \pm \delta }^{\text{DFA}} } \right] - \left( {\delta E_{\text{es}} \left[ {\rho_{ \pm 1}^{\text{DFA}} } \right] + \left( {1 - \delta } \right)E_{\text{es}} \left[ {\rho_{0}^{\text{DFA}} } \right]} \right) \\ = - \delta \left( {1 - \delta } \right)E_{\text{es}} \left[ {\rho_{ \pm 1}^{\text{DFA}} - \rho_{0}^{\text{DFA}} } \right] \le 0. \\ \end{gathered}$$
(22)

When combined with the negative denominator in Eq. (11), it leads to positive \(\Delta_{{_{ - \delta } }}^{{^{\text{GLC}} }}\) for DFA methods, with the maximum at δ = 0.5.

It should be pointed out that Eqs. (20) and (22) describe the GLC discrepancy of the \(E_{\text{ee}}^{\text{HF}}\) and \(E_{\text{es}}^{\text{DFA}}\) only approximately due to the following inconsistencies. All three terms \(E_{{{\text{ee}}, \pm \delta }}^{\text{HF}} ,\;\;E_{{{\text{ee}}, \pm 1}}^{\text{HF}} \;{\text{and}}\;E_{{{\text{ee}},0}}^{\text{HF}}\) are constructed from a common set of orbitals, determined at N = J ± δ, appropriate for the first term only, while to be consistent, the orbitals determined at (J ± 1) and J should be applied for remaining two terms. In Eq. (22), a combination of \(\rho_{ \pm 1}^{\text{DFA}}\) and \(\rho_{0}^{\text{DFA}}\) is taken for \(\rho_{ \pm \delta }^{\text{DFA}}\) [like in Eq. (3)], while this density should be constructed from orbitals determined at N = J ± δ.

The results of the global constancy condition testing, Eq. (12), are presented in Fig. 2. Since the constancy of \(\mu \left[ N \right] = \left( { \partial E\left[ N \right]/ \partial N} \right) = - \partial E_{ - \delta } /\partial \delta\) is tested, the discrepancy for the HF method is accounted mainly by \(\partial \Delta E_{{{\text{ee,}} - \delta }}^{\text{HF}} / \partial \delta = - \left( {1 - 2\delta } \right)J_{\text{ff}}\), see Eq. (20). When combined with the negative denominator in Eq. (12), this produces \(\Delta_{{_{ - \delta } }}^{{^{\text{GCC}} }}\) linearly decreasing, zeroing at δ = 0.5. Similarly, for DFA methods, Eq. (20) should be replaced by Eq. (22), leading to \(\Delta_{{_{ - \delta } }}^{{^{\text{GCC}} }}\) linearly increasing, zeroing at δ = 0.5, and quite large at end points of δ range, with the smallest absolute values for the rCAM-B3LYP and HF methods. The zeroing at δ = 0.5 is consistent with the transition state approach in which the ionization energy is determined by the highest occupied orbital energy of the intermediate state, \(E\left[ J \right] - E\left[ {J - 1} \right] \approx \varepsilon_{\text{HOMO}} \left[ {J - 0.5} \right]\) [52], see also Ref. [53].

Fig. 2
figure 2

The average of \(\Delta_{ - \delta }^{\text{GCC}}\), Eq. (12), vs. fraction δ

Full size image

As indicated in Eq. (9), the region of μ constancy is limited to 0 < δ < 1. Nevertheless, Fig. 2 includes results at δ = 0. With \(\varepsilon_{f}^{\text{DFA}} \left[ {J - \delta } \right]\) defined to be the frontier (highest occupied) orbital energy, the value of \(\Delta_{0}^{{^{\text{GCC}} }}\) equals \({ \lim } _{\delta \to 0} \Delta_{ - \delta }^{\text{GCC}}\). It is consistent with \(\Delta_{ - 0.01}^{{^{\text{GCC}} }}\) for all atoms. However, the GCC indicator calculated formally at δ = 1, \(\Delta_{ - 1}^{{^{\text{GCC}} }}\), does not equal \({ \lim } _{\delta \to 1} \Delta_{ - \delta }^{\text{GCC}}\) and it differs significantly from \(\Delta_{ - 0. 9 9}^{\text{GCC}}\) (see Supporting Information). This difference is due to different calculational procedures for determination of \(\varepsilon_{f}^{\text{DFA}} \left[ {J - \delta } \right]\): at δ < 1, the KS equations are solved self-consistently in the space of J occupied orbitals, while at δ = 1—in the space of (J − 1) occupied orbitals. The lowest virtual orbital energy is to be taken as \(\varepsilon_{f}^{\text{DFA}} \left[ {J - 1} \right]\).

The results of the local linearity condition testing, Eq. (13), are shown in Fig. 3. The shapes of plots are similar to those in Fig. 1, but the maximum discrepancies are about 50 times larger. This similarity can be understood with the help of the DFT mapping ρ → E[ρ] connecting the ground-state density and energy. The maximal discrepancy concerning the density in Fig. 3 correlates with the maximal discrepancy concerning the energy in Fig. 1. Of all DFA methods, the smallest discrepancy is observed for rCAM-B3LYP \(( {\Delta_{{_{{{ - 0} . 5}} }}^{\text{LLC}} = 0.0 1 5})\). The HF method is even better \(( {\Delta_{{_{{{ - 0} . 5}} }}^{\text{LLC}} = 0.0 0 6})\).

Fig. 3
figure 3

The average of \(\Delta_{ - \delta }^{\text{GCC}}\), Eq. (13), vs. fraction δ

Full size image

The results of the local constancy condition testing are shown in Fig. 4. The shapes of plots are similar to those in Fig. 2, provided \(\left| {\Delta_{ - \delta }^{\text{GCC}} } \right|\) would be plotted there instead of \(\Delta_{ - \delta }^{\text{GCC}}\). The magnitude of discrepancies is practically the same in both figures. This similarity can be viewed as an analog of similarity between the \(\Delta_{ - \delta }^{\text{GLC}}\) and \(\Delta_{ - \delta }^{\text{LLC}}\) plots, but now applied to derivatives with respect to N. Again, the performance of the rCAM-B3LYP functional is the best among tested approximate functionals.

Fig. 4
figure 4

The average of \(\Delta_{ - \delta }^{\text{GCC}}\), Eq. (15), vs. fraction δ

Full size image

4 Conclusions

In concluding, our illustrative examples’ calculations demonstrate that the performance of the rCAM-B3LYP functional is the best among the tested functionals. Moreover, satisfaction of inequalities (16) is quite reasonable \(\left( {\Delta_{ - 0.5}^{\text{GLC}} = 5 \times 10^{ - 5} ,\Delta_{ - 0.3}^{\text{GCC}} = - 4. 5 \times 10^{ - 4} ,\Delta_{ - 0.5}^{\text{LLC}} = 1. 5 \times 10^{ - 2} ,\Delta_{ - 0.5}^{\text{LCC}} = 3. 8 \times 10^{ - 3} } \right)\). The especially small \(\Delta^{\text{GLC}}\) is due to the fact that the rCAM-B3LYP functional was constructed to minimize just the discrepancy from the GLC. Another general conclusion can be reached: when the chemical potentials and the Fukui function are to be determined, the best results are found at δ = 0.5 for all tested DFAs. It is surprising that the performance of the extended HF method is as good (or even better) as that of DFT with the rCAM-B3lYP functional. Probably, this is due to the particular set of systems (atoms) chosen for tests. In any case, the hyper-HF method should not be used for large systems because important correlation effects are missing.