Generalized Kohn–Sham equations with accurate total energy and single-particle eigenvalue spectrum

Historically, the Kohn–Sham (KS) scheme was proposed to provide an improved approximation for the kinetic energy of an electronic system in its ground state as a functional of the density. As such, the (exact) KS equations were not expected to yield an accurate single-particle eigenvalue spectrum, except for the highest occupied (HOMO) eigenvalue, which is equal to the ionization potential (IP) of the system.¹ 

This view has changed over the years. Numerical inversion of accurate electronic densities to obtain a nearly exact KS potential has revealed that the eigenvalues of the valence KS orbitals are very close in magnitude (within ∼0.1 eV) to the corresponding vertical IPs of the system. Verma and co-workers^2,3 and Van Meer and co-workers⁴ have independently proven theorems predicting the observed closeness of the negative of the eigenvalues of the valence KS orbitals to the IPs of the electronic system. For deeper valence electrons, the coincidence between the negative of KS orbital eigenvalues and IPs deteriorates, and for core electrons, the difference reaches a few tens of electron volts.

Approximate local and semi-local density functionals yield rather accurate total energy predictions, but their corresponding KS potentials are not equally good. The typical KS HOMO eigenvalue error of local and semi-local density functional approximations (DFAs) is about ∼4 eV, relative to the exact KS result. The reason is simple: for an N-electron system, the electrons in the occupied KS orbitals are repelled via the Hartree–exchange–correlation (HXC) part of the approximate KS potential by an effective charge of N electrons rather than the correct N − 1 electrons of a self-interaction-free model. Consequently, the occupied orbital eigenvalues are erroneously shifted to higher energies, resulting in an underestimation of the IP by ∼4 eVs.

In our group, we have addressed this error with a constrained minimization of the total energy in local or semi-local DFAs. The DFA total energy remains almost unchanged, but the constrained minimization forces the effective “screening” (or “electron repulsion”) charge of the HXC potential of the DFA to equal N − 1.⁵ As a result, the negative of the HOMO eigenvalues from our constrained minimization method come closer to the true IPs (still underestimating them on average) with a typical error of about ∼1 eV. It is noteworthy that the resulting IP predictions are similar for the three different DFAs we tried [local-density approximation (LDA), Perdew–Burke–Ernzerhof (PBE), and B3LYP].⁶ 

In Hartree–Fock (HF) theory, Koopmans theorem predicts that the negative of the eigenvalues of all the occupied orbitals is equal to minus the corresponding vertical IPs. The theorem neglects orbital relaxation when an electron is ionized, and HF theory omits electron correlation altogether. These two errors have opposite signs for occupied orbitals⁷ and almost cancel each other. Still, the orbital relaxation error is stronger than the correlation error, and the dominance of relaxation increases for core electrons. The HF eigenvalues typically overestimate the vertical IPs from about ∼1 eV for the least bound electron (HOMO) to a few tens of eVs for core electrons.

The improvement of the HOMO eigenvalue with the constrained minimization of local/semi-local DFAs brings the magnitude of the IP error to the same ballpark as HF theory but with the opposite sign (Fig. 1). At the same time, the total energies from DFA constrained minimizations are virtually identical to the unconstrained DFA energies (within ∼0.1 meV),⁸ so the good quality of the DFA total energies remains intact in the constrained minimization. The aim of this paper is to consider a hybrid of the constrained minimization of local/semi-local DFAs and HF theory, aiming to formulate a generalized KS (GKS) scheme^9,10 with further improved single-particle spectrum over constrained DFA (CDFA) and HF theory for all occupied electrons, preserving at the same time the quality of the DFA total energy.

The rest of the paper is organized as follows: In Sec. II, we briefly review the CDFA method and modify it to aid its integration in a hybrid scheme and to make the implementation of our constraints more efficient. In Sec. III, we present the GKS formalism valid for any mixing between CDFA and HF. We seek the optimal mixing and find it is close to 50%, as expected from Fig. 1. In Sec. IV, we confirm that the predictions for the first IPs are better than those of either HF or CDFA and that the quality of the resulting total energies is of similar quality compared with those of the DFA (and CDFA). Finally, we investigate whether the flexibility of our hybrid potential corrected for self-interactions (SIs) allows for a better prediction for the IPs of all the bound electrons, including those in the core where relaxation effects are expected to become dominant.

In our previous work,^5,8 the approximate KS potential is constrained in order to remove the effect of SIs from its long-range asymptotic decay (without correcting for SIs the total energy). For a finite system of N electrons, the HXC potential should decay as (N − 1)/r, but in many common approximations, it incorrectly decays as N/r, a failure that we attribute to SIs.

To correct the potential for SIs, the KS equations are written as

where v_rep takes the role of the HXC potential, v_HXC. The potential v_rep is written as the Coulomb–Hartree potential of an auxiliary density ρ_rep as follows:

We have called this density the “effective repulsion” or “screening” density in the KS system. Mimicking in a mean field way the absent repulsion of the real electrons, the KS screening density is distinct from the physical density of the interacting electron system and is subject to the following constraints:

The application of these constraints ensures that the potential v_rep has the correct asymptotic decay of the exact HXC potential. In addition, for a single electron, the effective repulsion potential becomes zero everywhere, v_rep(r) = 0, which is the exact one electron result for the HXC potential. Even though we employ v_rep to approximate v_HXC, we avoid writing v_rep = v_HXC because so far, we have not been able to write v_rep satisfying (3, 4) as the functional derivative of an approximate HXC energy. Görling¹¹ was the first to employ the ansatz (2, 3) for the exact-exchange potential, constraining the exchange density to integrate to −1.

The first constraint in (3) is incorporated via a Lagrange multiplier in the objective functional

In Refs. 5 and 8, the second constraint in (4) is enforced through the use of a penalty function. However, our recent work¹² provides an alternative by constructing the density ρ_rep(r) as the modulus square of an effective repulsion, or screening, amplitude,

To obtain the ground state, the objective functional in Eq. (5) is minimized with respect to ρ_rep. The minimization can be performed as an extension to the optimized effective potential (OEP) method.^13,14 [When ρ_rep is given by Eq. (6), the minimization is with respect to f. Details of the minimization of the objective function with respect to the amplitude f in (5) are presented in a separate publication.] At the minimum of G,

where $vHXCDFA(r)=δEHXCDFA[ρ]/δρ(r)$ is the DFA HXC potential and

is the KS density–density response function. The computational cost to solve the OEP in Eq. (7) enforcing (3, 4) is low, scaling similar to a few (∼10) KS-DFA calculations, since the potential $vHXCDFA(r)$ in (7) is local.

This method can be applied to a generic energy density functional E[ρ]. In Ref. 8, we applied this method to a range of functionals, including a hybrid functional, for which the HXC potential is nonlocal. The constrained minimization increased total energies minimally, while improving calculations of the IP from the negative of the HOMO eigenvalue. It was also found that the HOMO eigenvalues were largely independent from the functional used to calculate them.

In the original formulation of the constrained method, the sum of the Hartree, exchange, and correlation potentials is subject to the constraints (3) and (4). Alternatively, these constraints can instead be applied only to the Hartree and exchange (HX) potential. To leading order, the long-range decay of the exact HX potential is (N − 1)/r, while the exact correlation potential decays as ∝−1/r⁴.^15,16 Therefore, the asymptotic behavior of the exact v_HXC will be dominated by that of the HX potential v_HX. The exact, leading-order asymptotic behavior of the potential is, therefore, maintained if the unconstrained correlation potential is included separately to the constrained HX potential. Correlation can be treated separately to the constrained method by constructing the KS equations as

where v_c(r) is the density functional theory (DFT) correlation potential δE_c[ρ]/δρ(r). Treating correlation separately to Hartree–exchange, represented again by v_rep, the equation to determine ρ_rep becomes

identical to Eq. (7), except that the potential $vHXCDFA$ is now $vHXDFA$⁠.

Table I shows that excluding correlation from the constraint results in an upshift of the calculated IP. For LDA, this results in an average upshift of 0.90 eV, and for PBE, this average is 0.42 eV. Typically, calculations of the HOMO energy, similar to unconstrained LDA calculations, underestimate the IP. This upshift of the calculated IP, therefore, after leaving correlation unconstrained, generally improves slightly the quality of the approximation for the IP of the system. However, the removal of correlation from the constrained minimization results in a KS potential that is no longer appropriate for a system with a single electron. The exact behavior for a one electron system is only recovered if the approximate correlation potential is zero for a single electron, which holds for the exact functional.

Separating the correlation potential from the constrained minimization allows for the calculation of a constrained HX component of the KS potential, where the main effects of SIs have been removed. To construct our hybrid method, this potential is combined with the SI-free nonlocal HF potential in a GKS scheme.^9,10 A simple generalized KS equation can be constructed using a single parameter α to control the contribution of local and nonlocal exchanges. The hybrid, single-particle potential corresponding to HXC is written as follows:

where v_rep is the effective repulsion, or screening, potential given by Eq. (10) (with χ_v constructed using the GKS orbitals), $vHXHF$ is the sum of the Hartree plus nonlocal Fock–exchange potentials, and $vcDFA$ is the DFA correlation potential. Correspondingly, the total energy is given as

where T_s is the noninteracting kinetic energy, E_en is the energy from the nuclear potential, E_H is the Hartree energy, $ExDFA$ and $EcDFA$ are the exchange and correlation energy functionals, and $ExF$ is the Fock–exchange energy.

The treatment of correlation separately to exchange allows for the new hybrid method to have a consistent correlation energy for all values of α. The total energy varies almost linearly from the HF plus DFA correlation limit (α = 0) to the constrained DFA energy limit (α = 1) as the hybrid parameter α varies from 0 to 1. We look to choose a DFA for which this change in total energy is small such that total energies are insensitive to the choice of α. This allows the parameter α to be varied to optimize the prediction of ionization energies for this hybrid. For LDA, Table II demonstrates a significant change of 1.11 eV in the total energy as α varies over 0 ≤ α ≤ 1. This large energy change can be explained because the HF energy is lower than the LDA total energy, and adding correlation energy, a negative quantity, on top of HF at α = 0 further lowers the total energy and increases the energy difference between the end points, α = 1 (LDA) and α = 0 (HF + LDA correlation).

The constrained PBE (CPBE) approximation shows a smaller energy shift of about 0.01 eV (Table II) between the constrained PBE total energy (α = 1) and the HF plus PBE correlation energy (α = 0). To exploit the insensitivity of the total energy on the hybrid parameter α, we restricted our investigation to the hybrid method with constrained PBE exchange. Because the energetics of the CPBE hybrid remain consistent throughout the range 0 ≤ α ≤ 1, the effect of varying α on the orbital energies can be investigated while the quality of the energetics is not affected. As we found in our previous work,⁸ the ionization energy calculated from the constrained method is insensitive to the choice of DFA, and therefore, an optimization of the ionization energies is likely to be valid for a range of exchange functionals.

With our choice of CPBE, the total energies of our hybrid scheme are relatively unchanged for a range of values of α. This allows the optimization of the hybrid parameter α for the calculation of ionization energies. However, instead of focusing just on the HOMO eigenvalue,^4,17 we optimized the present hybrid method to best approximate the ionization energies of all occupied orbitals. The value of the hybrid parameter was optimized for the simple molecules in Table I with molecular geometries taken from the supplementary material of Ref. 18 along with experimental and coupled cluster results. The minimizing α can be determined from the results in Fig. 2. These results show that the minimizing α varies slightly when considering the HOMO, core, and valence orbitals. The energy eigenvalues of the valence orbitals are optimized at a slightly larger value of α, i.e., a smaller percentage of HF exchange, while IPs from core orbital eigenvalues are better approximated with a slightly smaller percentage of constrained DFT exchange. A choice that achieves good accuracy for the full range of orbitals without overfitting the results is the value α = 0.5. This value coincides with the half-and-half mixing of the original hybrid formulation of Becke¹⁹ based on a linear interpolation of the adiabatic connection formula. The present 50% hybrid of CPBE and HF exchange (which we denote “constrained hybrid”—CHYB) is the self-interaction-corrected analog of the PBE50 functional²⁰ that was employed in time-dependent DFT and many-body perturbation theory methods. In DFT applications, the PBE0 hybrid functional^21,22 with 25% HF exchange is a more widely used parameterization.

The eigenvalue accuracy of the constrained hybrid method with α = 0.5 will now be investigated. Figure 3 demonstrates the agreement of the CHYB method with experimental results for the HOMO (blue), valence (red), and core (inset) orbitals. We see a clear improvement of the CHYB method over the CPBE and HF methods. Table III shows a comparison of our new CHYB method with results from unconstrained and constrained PBE, PBE0, HF, and coupled-cluster singles and doubles (CCSD) methods. These results demonstrate that the constrained hybrid method has a significantly reduced error over the other DFT methods for the calculation of all IPs and has accuracy approaching CCSD results. Ionization energies of core orbitals show a significant improvement over PBE and HF results. These ionization energies are largely independent of the surrounding chemical environment, and as such, the four groupings of points in Fig. 3 (inset) correspond (from left to right) to C, N, O, and F. As shown in Table III, the error for valence orbitals is, in general, positive; for core orbitals, this error becomes negative. This is likely caused by the HF error increasing more relative to CPBE due to the neglect of relaxation effects. The core orbitals especially are expected to have a large relaxation error and the HF error to increase. Of the molecules investigated with this hybrid, outlying results are for the molecules O₃ with one valence orbital underestimated by 3 eV and SiF₄ where every orbital has errors over 1 eV. The large error on these molecules requires further investigation, although similarly poor results were also found for the CCSD calculations for O₃.¹⁸ As would be expected from the imposition of additional constraints, the calculated total energy is higher than the unconstrained 50% PBE hybrid and in line with the increase in energy when applying the constrained method.⁸ An increase in the total energy of ∼0.1 meV was found for the hydrogen fluoride molecule. This small increase was maintained for a range of interatomic separations, showing that the addition of the constrained method will have minimal impact on the energetic behavior compared to the unconstrained method. This hybrid will, therefore, have the same energetic behavior as the 50% PBE hybrid method.

Table IV shows the full orbital energies for a selection of molecules compared with experimental results. The agreement with experiment of the self-interaction-corrected hybrid functional can clearly be seen throughout these molecules. Core orbitals have the largest absolute error due to their large eigenvalues; however, even for these orbitals, the constrained hybrid method has average errors smaller than 1 eV. The present method demonstrates a significant improvement in calculating the spectra of these molecules when compared to HF, PBE0, and CPBE methods. The computational cost of a constrained hybrid calculation scales similarly to a small number of conventional hybrid calculations.

From the results presented in this paper, it is clear that the constrained hybrid method offers significant improvements in the approximation of orbital energies not just for the HOMO but for the full range of occupied orbitals. The errors in orbital energies, below 1 eV for valence orbitals and within 1–2 eV for core orbitals, are surprisingly small considering the semi-local DFA we employ. We attribute the improvement to the hybrid character of the exchange potential whose equally weighted components are either SI-free (nonlocal Fock term) or corrected from the main errors of SIs (constrained DFA exchange potential). The total energy of this constrained hybrid remains relatively consistent with the PBE and HF plus PBE correlation methods and is increased over an unconstrained hybrid method by ∼0.1 meV. The accuracy of these results, especially for the core orbitals, suggests that there is a cancellation of errors between the exchange potential in the constrained minimization method and HF. For orbitals with large ionization energies, there should be significant orbital relaxation effects, which, nevertheless, are accounted for in the orbital energy spectrum of the constrained hybrid method. As we found in Ref. 8, the IPs calculated from the constrained method are largely independent of the DFA. Therefore, the constrained hybrid method will show similar improvements in the calculation of IPs for a range of DFAs using a mixing parameter of 50%. This mixing is in line with a simple linear interpolation of the adiabatic connection, and while further improvements for specific orbitals may be obtained for a highly optimized value of α, the choice of α = 0.5 is expected to work for a variety of systems and orbitals. Finally, the accuracy of these results is consistent with recent arguments^3,18,23,24 that the DFT Koopmans’ theorem holds approximately for all orbitals, although to achieve the systematic high accuracy of the present results, a hybrid exchange potential is necessary.

T.C.P. and N.I.G. acknowledge financial support from The Leverhulme Trust through a Research Project Grant under Grant No. RPG-2016-005. N.I.G. thanks Professor Rod Bartlett for helpful discussions during his visit to Durham University in early 2019 and acknowledges the Institute of Advanced Study at Durham University for hosting this visit.

The authors state that they have no conflict of interest to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.