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

We propose a new generalised Kohn-Sham or constrained hybrid method, where the exchange potential is the (equally weighted) average of the nonlocal Fock exchange term and the self-interaction-corrected exchange potential, as obtained from our constrained minimisation method of semi-local approximations. The new method gives an accurate single-particle eigenvalue spectrum with an average deviation between (the negative of) the valence orbital eigenvalues and the experimental ionisation potentials of about 0.5eV, while the deviation of core orbitals is within 2eV. The improvement in the eigenvalue spectrum is achieved with a minimal increase in the total energy.


Introduction
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 ionisation potential (IP) of the system [11]. 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.1eV) to the corresponding vertical IPs of the system. Bartlett and coworkers [19,4] and Baerends and coworkers [18] have independently proven theorems predicting the observed closeness of the negative of the eigenvalues of 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-interactionfree model. Consequently, the occupied orbital eigenvalues are erroneously shifted to higher energies, resulting in an underestimation of the IP by approximately ∼4 eVs.
In our group, we have addressed this error with a constrained minimisation of the total energy in local or semilocal DFAs. The DFA total energy remains unchanged but the constrained minimisation forces the effective "screening" (or "electron repulsion") charge of the HXC potential of the DFA to equal N − 1 [8]. As a result, the negative of the HOMO eigenvalues from our constrained minimisation 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 (LDA, PBE, B3LYP). 1 In Hartree Fock (HF) theory, Koopmans theorem predicts that the negative of the eigenvalues of all the occupied orbitals are equal to minus the corresponding vertical IPs. The theorem neglects orbital relaxation when an electron is ionised and HF theory omits electron correlation altogether. These two errors have opposite signs for occupied orbitals [3] 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 minimisation 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 minimisations are virtually identical to the unconstrained DFA energies (within ∼0.1 meV) [13], so the good quality of the DFA total energies remains intact in the constrained minimisation. The aim of this paper is to consider a hybrid of the constrained minimisation of local/semi-local DFAs and of HF theory, aiming to formulate a generalised KS (GKS) scheme 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 organised as follows: in section II we review briefly 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 section 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 section 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.

Constrained minimisation of the total energy
In our previous work [8,13] the approximate KS potential is constrained in order to remove the effect of SIs from its longrange 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 : 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. Additionally, 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. The first constraint in (3) is incorporated via a Lagrange multiplier in the objective functional In Refs. [8,13] the second constraint in (4) is enforced through the use of a penalty function. However, our recent work [12] 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 minimised with respect to ρ rep . The minimisation can be performed as an extension to the optimised effective potential (OEP) method [15,16]. (When ρ rep is given by Eq. 6, the minimisation is with respect to f . Details of the minimisation of the objective function with respect to the amplitude f in (5) are presented in a separate publication.) At the minimum of G, where the potential v DFA HXC (r) = δ E DFA HXC [ρ]/δ ρ(r) is the DFA HXC potential and is the KS density-density response function. This method can be applied to a generic energy density functional E[ρ]. In Ref. [13] we applied this method to a range of functionals, including a hybrid functional for which the HXC potential is non-local. The constrained minimisation 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 (34). 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 4 [17,2]. Therefore the asymptotic behaviour of the exact v HXC will be dominated by that of the HX potential v HX . The exact, leading-order asymptotic behaviour 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 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 the potential v DFA HXC is now v DFA HX . Table 2 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, similarly 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 minimisation results in a    = 0), for various values of α, for the LDA and PBE functionals and a range of molecules. These results are the average of the differences in total energies for the molecules in Table 2. Energies in eV.
KS potential that is no longer exact for a system with a single electron. The exact behaviour 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.

The Generalised KS Method
Separating the correlation potential from the constrained minimisation allows 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 non-local HF potential in a GKS scheme. A simple generalised KS equation can be constructed using a single parameter α to control the contribution of local and non-local exchange where v F x is the nonlocal Fock exchange potential, v DFA c is the DFA correlation potential and v eff is the self-interaction corrected exchange potential, (difference of v rep and Hartree potentials, v eff = v rep − v H ), with v rep given by Eq. 10. Correspondingly the total energy is given as where T s is the non-interacting kinetic energy, E en is the energy from the nuclear potential, E H is the Hartree energy, E DFA x and E DFA c are the exchange and correlation energy functionals, and E F x 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 LDA 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 optimise the prediction of ionisation energies for this hybrid. For LDA, Table 2 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 lowers further 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 2) 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 [13] the ionisation energy calculated from the constrained method is insensitive to the choice of DFA and therefore an optimisation of the ionisation energies is likely to be valid for a range of exchange functionals.

Results
With our choice of CPBE, the total energies of our hybrid scheme are relatively unchanged for a range of values of α. This allows the optimisation of the hybrid parameter α for the calculation of ionisation energies. However, instead of focusing just on the HOMO eigenvalue [18,20], we optimised the present hybrid method to best approximate the ionisation energies of all occupied orbitals. The vaue of the hybrid parameter was optimised for the simple molecules in Table 2 with molecular geometries taken from the supplementary material of Ref. [14] along with experimental and coupled cluster results. The minimising α can be determined from the results in Fig. 4. These results show that the minimising α varies slightly when considering the HOMO, core, and valence orbitals. The energy eigenvalues of the valence orbitals are optimised 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-andhalf mixing of the original hybrid formulation of Becke [5], 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 selfinteraction corrected analogue of the PBE50 functional [6] that was employed in time-dependent DFT and many-body perturbation theory methods. In DFT applications, the PBE0 hybrid functional [10,1] with 25% HF exchange is a more widely used parameterisation.
The eigenvalue accuracy of the constrained hybrid method with α = 0.5 will now be investigated. Figure 4 demonstrates the agreement of the CHYB method with experimental results for the HOMO, valence and core orbitals. We see Figure 2: The orbital energy "error" (difference of orbital eigenvalue magnitude from experimental IP) as a function of α. We show average error (AE) and absolute average error (AAE) for all the valence orbitals and separately for the HOMO. In the inset, we show the average error and absolute average error for the core orbitals. Results are calculated for the molecules in Table 2 using a cc-pVDZ basis set. a clear improvement of the CHYB method over the CPBE and HF methods. Table 4 compares our new CHYB method with results from unconstrained and constrained PBE, HF and coupled-cluster singles and doubles(CCSD) results. These results demonstrate that the constrained hybrid method has a significantly reduced error over the other DFT methods for calculation of all IPs, and has accuracy approaching CCSD results. Ionisation energies of core orbitals show a significant improvement over PBE and HF results. These ionisation energies are largely independent of the surrounding chemical environment and as such the four groupings of points in Figure 4(inset) correspond (from left to right) to C,N,O,F. As can be seen from table 4, 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 3 with one valence orbital underestimated by 3eV and SiF 4 where every orbital has errors over 1eV. The large error on these molecules requires further investigation, although similarly poor results were also found for the CCSD calculations for O 3 [14]. 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 [13]. An increase in the total energy of ∼ 0.1meV 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 behaviour compared to the unconstrained method. This hybrid will therefore have the same energetic behaviour as the 50% PBE hybrid method.  Table 4 using the cc-pVDZ basis set. Coloured points correspond to the hybrid method, grey points show results for CPBE (crosses) and HF(diamonds), the black line shows the ideal exact agreement between experimental and calculated results. Inset, shows the same results including high energy core orbitals. Table 4 shows the full orbital energies for a selection of molecules, compared with experimental results [7]. 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 1eV. The present method demonstrates a significant improvement in calculating the spectra of these molecules when compared to HF and the CPBE methods.

Conclusions
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 1eV 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 (non-local 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-PBEcorrelation methods, and is increased over an unconstrained hybrid method by ∼ 0.1meV. The accuracy of these results, especially for core orbitals suggest that there is a cancellation of errors between the exchange potential in the constrained minimisation method and HF. For orbitals with large ionisation 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. [13], the IPs calculated from the constrained method is 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 optimised 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 [7,9,14,4] 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.