Improving the exchange and correlation potential in density functional approximations through constraints

We review and expand on our work to impose constraints on the effective Kohn Sham (KS) potential of local and semi-local density functional approximations. In this work, we relax a previously imposed positivity constraint, which increased the computational cost and we find that it is safe to do so, except in systems with very few electrons. The constrained minimisation leads invariably to the solution of an optimised effective potential (OEP) equation in order to determine the KS potential. We review briefly our previous work on this problem and demonstrate with numerous examples that despite well-known mathematical issues of the OEP with finite basis sets, our OEP equations are well behaved. We demonstrate that constraining the screening charge of the Hartree, exchange and correlation potential not only corrects its asymptotic behaviour but also allows the exchange and correlation potential to exhibit nonzero derivative discontinuity, a feature of the exact KS potential that is necessary for the accurate prediction of band-gaps in solids but very hard to capture with semi-local approximations.


I. INTRODUCTION
A challenge with common density functional approximations is the imbalance of accuracy between the energy functionals and the corresponding Kohn-Sham (KS) potentials, i.e. the functional derivatives of the energy density-functionals. Although the accuracy and quality of an energy density-functional is often quite good, the resulting KS potential is inferior [1][2][3]. The quest to derive ever more accurate energy density-functionals to obtain moderate improvements on the KS potential may not be the best strategy (it is vulnerable to diminishing returns in the accuracy of the KS potential). We explore different routes to improved accuracy for these calculations.
Previously, we explored the minimisation of potential functionals defined by an energy difference, instead of density fucntionals of the total energy, as a means of improving the quality of the KS potential [4][5][6][7]. The advantage of this approach is that the energy difference is bound from below, even in approximations from finite-order (second) perturbation theory; the latter can then be employed directly to derive accurate exchange-correlation (xc) potentials without the risk of variational collapse [4,6].
In this paper, we review briefly and expand on our work [8][9][10] to improve the performance of local and semi-local density-functional approximations (DFAs), by imposing physical constraints on the single-particle, local, effective (KS) potential, whose orbitals minimise the total energy functional. In Refs. 8-10 we argued that these constraints improve the asymptotic behaviour and overall quality of the KS potential by removing the erroneous effects of self-interactions (SIs). As evidence, we demonstrated that, compared with the results

II. CONSTRAINED MINIMISATION OF DENSITY FUNCTIONAL APPROXIMATIONS
In the constrained minimisation method [8][9][10] we employ the standard total energy expression in DFT, using a density functional approximation (DFA) for the xc energy density- The various quantities have their usual definitions, v en is the external potential, are the noninteracting kinetic energy and Hartree energy density functionals. Following the optimised effective potential method (OEP) [12,13], we set that the KS orbitals satisfy single-particle KS equations employing an effective potential v(r), The total energy is then minimised by imposing constraints, detailed below, on the effective potential v(r). The effective potential v(r) is akin to the Hartree-exchange and correlation (Hxc) potential of KS theory v DFA Hxc (r). However, the constraints we impose correct the asymptotic form of v(r) and alleviate other effects of SIs from it, so in general, v(r) = v DFA Hxc (r).
In the constrained method, we treat the Hxc screening density, or electron repulsion density, ρ scr (r), [52], as the fundamental quantity. It is defined via Poisson's equation from the Laplacian of the difference of the (exact or approximate) KS potential minus the external potential, ∇ 2 v s (r) − v en (r) ; for example, the Hxc screening density of the exact KS potential is given by, Together with the integrated Hxc screening charge Q scr , Q scr = dr ρ scr (r), the Hxc screening density plays a central role in our constrained method to mitigate against the effects of self-interactions. The concept of an effective screening density was first explored in Refs. 14-16 in terms of a screening density for the xc (or exchange only) potential; for the exact xc-potential, the screening charge is Q xc scr = −1 [14][15][16]. It has been used in various applications of the OEP method to fix the freedom of a constant in the OEP solution [14,[17][18][19][20].
The Hxc-screening charge of the exact KS potential satisfies the intuitive sum rule [14][15][16], However, in common DFAs (such as the local density approximation, L(S)DA, and most generalized-gradient approximations (GGAs)) this sum rule is violated and the screening charge is in fact given by Q scr = N . We argue [8] that this violation of the sum rule can be attributed to the presence of SIs, since it implies that any of the electrons of an N -electron system are effectively repelled, via the Hxc potential, by a net charge of N electrons. We note that the sum rule (5), which depends on the screening density and is violated for LDA and common GGAs, is different from the well-known sum rule [21] for the xc hole, dr ρ xc (r, r ) = −1, which is satisfied by LDA and common GGAs. The quantities ρ xc scr (r) and ρ xc (r, r ) are not directly related.
Accordingly, in the constrained minimisation of DFAs [8][9][10] (which we henceforth refer to as the CDFA method), our strategy to mitigate the effects of SIs from the effective potential is to enforce that the KS orbitals satisfy Eq. (2) with the effective potential v(r) represented by the effective screening density v(r) = dr ρ scr (r ) |r − r | , where ρ scr (r) satisfies two constraints: and ρ scr (r) ≥ 0.
The second constraint (8) is physically intuitive, hinting at interpreting ρ scr (r) as the charge density of N − 1 electrons. However, this condition is too restrictive and not satisfied by the exact KS potential.
Nonetheless, the positivity constraint (8) has a double role in the constrained minimisation method. As explained in Refs. 8-10, the CDFA minimization procedure must be solved within the optimized effective potential (OEP) framework [12,13]. Primarily, the positivity constraint allows the mathematical problem of constrained minimisation to remain well posed in the limit of complete orbital and auxiliary basis sets [8,9]; without the positivity constraint, there is nothing to prevent the screening density from separating into a component in the energetically important spatial region near the molecule, with charge Q a scr = N , and a separate component with charge Q b scr = −1 pushed out to infinity (within the basis set limits). Secondly, the solution of the OEP equation in Gaussian basis set codes is a longstanding problem in DFT, which has hindered the widespread adoption of OEPbased methods in practical calculations. Various methods have been developed to overcome these numerical difficulties which typically manifest themselves as spurious oscillations in the effective potential [17,18,22]. With finite orbital and auxiliary basis sets, the positivity constraint (8) offers a simple way to reduce drastically the variational flexibility of ρ scr (r) and of v(r) and thereby helps to overcome mathematical pathologies in the solution of the OEP equation.
In the previous implementation of the CDFA method, the positivity constraint was used in combination with a singular value decomposition (SVD) of the density-density response when this energy term is written as an implicit functional of the density. Alternatively, it can be obtained by minimising the DFT total energy expression, (1), indirectly by searching for the effective potential v(r) in (2) whose KS orbitals minimise the total energy [26,27].
Either way, we obtain the integral OEP equation, where χ v (r, r ) is the density-density response function given by (in a complete orbital basis set) The KS orbitals from (2) are assumed to be real-valued. The right-hand side (RHS) b v (r) depends on the energy term whose functional derivative we take, in our case the Hxc energy If no constraints are imposed, the solution of (9) is trivially v(r) = v DF A Hxc (r) within a constant, since χ v (r, r ) has no null eigenfunctions except the constant function. In Ref. 10 we explain how we impose the normalisation constraint (7) on the effective screening density and demonstrate that the the scheme can be applied for any given DFA, including LDA, GGAs and hybrid functionals.
To understand the effect of finite orbital basis sets on the solution of the OEP equation, we focus on the density-density response function; the analysis below also applies to the RHS b v (r). We split χ v into two terms, the first of which can be represented exactly in the orbital basis, and the second which must be approximated. χ v is given, for λ = 1, by with The sum is over occupied {φ i } and unoccupied {φ a , φ b } KS orbitals (2)  By definition, the complementχ v cannot be represented exactly so we must approximate it. We use theÜnsold approximation [28] together with the completeness relation for the KS orbitals (in much the same manner as the well-known Krieger-Li-Iafrate (KLI) approximation [29,30] and common energy demoninator approximation (CEDA) [31,32] methods), in which caseχ v reduces tō where −∆ is the common energy denominator that replaces i − b in (14), ∆ > 0. In Eq. (15), we omit the final term with the same domain as χ 0 v , because its contribution to χ λ v vanishes for small λ, which is ultimately the limit we seek.
We observe that, as long as ∆ > 0, the value of ∆ does not play a role in the results, since ∆ always appears together with λ, in the ratio λ/∆, and we investigate the limit λ → 0.
We shall also consider the limit λ → ∞, for which the value of positive ∆ does not matter either. It is straightforward to confirm thatχ v is negative semi-definite, like χ 0 v , and that the only null eigenfunction ofχ v is the constant function.
The same procedure is applied for the RHS b v (r) of the OEP equation (10), which yields the following expressions for the terms b 0 v (r) and its complementb v (r), The OEP equation thus takes the following form, To solve this equation in a Gaussian basis set code, the screening density is expanded in an auxiliary basis set and its coefficients can be found by a straightforward matrix inversion.
The screening charge constraint (7) besides mitigating against SI errors is also necessary to fix the freedom of a constant in the effective potential [17] and is enforced using a Lagrange multiplier. The optimization procedure is explained in detail in Ref. 10; the only difference here is that the matrices for the LHS and RHS of the OEP equation now contain the additional complement terms.
Prior to Ref. 23, the finite orbital basis OEP was given by the solution of (18) at λ = 0.
However, this solution leaves the effective potential v 0 (r) indeterminate in the null space of In order to obtain a smooth potential, v 0 (r), one must restrict the freedom of v 0 (r), which has spawned a variety of approaches in the literature. These include, for example, schemes to balance the relative sizes of the orbital and auxiliary basis set [18,19]; regularization techniques to smooth the effective potential [33,34]; and removing the additional freedom of an auxiliary basis set [35]. In our method, rather than restricting the freedom of v 0 (r), we instead solve the OEP equation to find the potential v λ (r) which is defined mathematically to be unique for finite λ.
The main point of Ref. 23 is the observation that the solution of the same equation (18) for any finite λ > 0, even λ tending to zero, determines the effective potential fully, up to a constant. The extension of the response function withχ v amounts to using an effectively complete orbital basis. Numerically, we find that the solution of (18) is smooth for almost any λ > 0 [53], including the limits for small and for large λ, which correspond respectively to the OEP potential in a finite orbital basis, v λ→0 (r), and to its (Unsöld) approximation with a common energy denominator, v ∞ (r). It turns out that for the effective xc potentials in the constrained minimisation method, the two solutions are close to each other.

A. Relaxing the positivity constraint
In Refs. 8-10 we solved the OEP equation for CDFA method, using finite orbital and auxiliary basis sets, with λ = 0. The indeterminacy of the effective potential was restricted by expressing v(r) in terms of the screening density ρ scr (r) in (6) and then constraining the screening charge Q scr (7) as well as the sign of ρ scr (r) (8).
However, the positivity constraint, implemented with a penalty function [10] is a computational bottleneck for the calculation. In a forthcoming paper, we implement the positivity constraint more efficiently, by writing ρ scr (r) = |f scr (r)| 2 , and solving for the screening amplitude f scr (r) [36], which ensures the constrained minimization is mathematically well-posed regardless of basis set size.
In the next part, we investigate the effects of relaxing the positivity constraint on the convergence of the screening potential and screening density. A weak effect, for sufficiently flexible auxiliary basis sets, will justify the relaxation of the positivity constraint and reduce the computational effort. The auxiliary basis sets we use are un-contracted cc-pVXZ [37,38], with X=D,T,Q.
In the rest of the section, we show indicative results for the CDFA method applied to the LDA functional, henceforth denoted by CLDA, where the minimisation was performed under just the constraint for the screening charge, Q scr = N − 1 (7). In order to determine v(r) and ρ scr (r), we employ the extended response function χ λ v (r, r ), in the limit of small λ. We use λ/∆ = 0.01, but the results seem converged and do not change if we reduce λ/∆ by an order of magnitude. The positivity constraint enabled the constrained minimisation problem to remain well posed in the limit of large (complete) orbital and auxiliary basis sets. Consequently, we expect the screening charge to change gradually, as we increase the size of the auxiliary basis. This effect will be stronger for systems with few electrons, since then, the difference between N − 1 and N is largest.
Calculations were performed in the Gaussian basis set code HIPPO [54], with oneand two-electron integrals for the Cartesian Gaussian basis elements calculated using the GAMESS code [39,40]. Basis set data was obtained from the Basis Set Exchange database [41]. the LDA potential v LDA xc (r) is also shown with a blue dashed line. In Figs. 2a-2c, r 2 ρ scr (r) (CLDA screening density multiplied by r 2 ), is shown for the Ne atom and three auxiliary basis sets, un-contracted cc-pVXZ, with X=D,T,Q. In each sub-figure r 2 ρ scr (r) is shown for fixed auxiliary basis and various orbital basis sets: cc-pVXZ, with X=D,T,Q,5. The overall convergence of the xc potential is excellent. The convergence of ρ scr (r) for fixed aux basis and increasing size of orbital basis is also very good.
Before proceeding, it is worth noting that despite not deploying the positivity constraint   basis and increasing orbital basis size is good. Note that for any combination of orbital and auxiliary basis, the xc potential is smooth. The convergence of the screening density for fixed auxiliary basis and increasing size of orbital basis is slower than the other systems. In addition, as the size of the auxiliary basis increases, from 6a to 6b to 6c, the screening density keeps changing considerably. Note specifically the negative part of the screening density in Figs. 6a-6c. In Fig. 6a the negative lump is centred around 2.5 a 0 away from the origin, in Fig. 6b it is centred around 3.0 a 0 away from the origin and in Fig. 6c it has moved to 3.5 a 0 .
This is the effect we discussed in section II. The positivity constraint enables the constrained minimisation problem to remain well posed for large basis sets (here large auxiliary bases). With only the constraint on Q scr enabled and without positivity, it becomes energetically preferable, during the total energy minimisation, to converge to a screening density with the screening charge locally equal to N (=Q LDA scr ), and to shift negative charge density away from the system. The effect is already evident for the moderately large auxiliary bases used in our study, because the difference between N − 1 and N is relatively large for N = 2.  It is evident that both the CLDA xc potential and the CLDA screening density are well converged and these systems do not present a challenge to the constrained minimisation, at least regarding convergence.
In Table I  In concluding this section, we first recall the reasons why our CDFA method was implemented with the positivity constraint (8). This constraint is intuitive if one considers each electron to experience a repulsive electronic density from the other N − 1 electrons, but it also serves two computational purposes: (i) to avoid shifting negative screening density to infinity as the size of the orbital and auxiliary basis sets increase and (ii) as a regularization technique to avoid pathological behaviour of the OEP solution. As we have seen from the good convergence of the screening densities and potentials, the latter reason is no longer necessary with the introduction of the complement terms in the OEP equation (18).
Regarding the first reason (i), we note that for the moderately large auxiliary basis sets we tested, it is safe to carry out constrained minimisations of the DFA total energy under the constraint of the screening charge only, Q scr = N − 1, except for systems with few electrons;   In the next section, we shall argue that the screening charge constraint endows the xc potential of local and semi-local DFAs with exotic qualities, such a finite derivative discontinuity ∆ xc . Although crucial for the accurate prediction of band gaps, ∆ xc is notoriously hard to capture in approximations. Advanced approximations have been proposed which capture this discontinuous behaviour, e.g., Refs. 42-46, however, further development is required for these methods to yield reliable band gaps for all materials.

III. DERIVATIVE DISCONTINUITY OF THE CDFA XC POTENTIAL
The discontinuity of the xc potential is defined by where v N ±ω xc (r) is the xc potential of an ensemble with N ± ω electrons.
The ensemble KS densities with N ± ω electrons are given by, where ρ M ven (r), M = N −1, N, N +1, is the ground state density of the M -electron KS system in the external potential v en (r). We shall use the CLDA KS equation (2), with constraints (7,8).
We seek the derivative discontinuity ∆ xc of the CLDA xc potential from (19) for reference.
In order to obtain ∆ ω xc (r) and then ∆ xc , one must first find the ensemble KS xc potentials with densities ρ N ±ω ven (r) and subtract them. Work is in progress in our group to obtain directly these ensemble KS xc potentials. Here, we use the method of Refs. 47 and 48 to obtain the ensemble KS xc potential by constructing the ensemble density ρ N ±ω ven from separate KS calculations for N , and N ± 1 particles and then inverting ρ N ±ω ven (r) to obtain v N ±ω xc (r).
Let us follow this construction in detail. The two KS ground state densities that build the ensemble density ρ N +ω ven (r) can be written: The notation makes explicit that {φ i [ρ M ](r)} are the KS orbitals of the M -electron system with density ρ M .
In terms of the ensemble KS orbitals {φ i [ρ N +ω ](r)} the ensemble density is given by In addition, from Eqs. 20-23, it is also equal to In general, the ensemble KS orbitals, {φ i [ρ N +ω ](r)} in (24) In the limit ω → 0 + , we have: where Q N + scr = lim ω→0 + Q N +ω scr . This stepwise increase of screening charge obviously causes a discontinuous jump in the constrained xc potential v N +ω xc (r). In the limit ω → 0 + , the jump of the xc potential is where v N + xc (r) = lim ω→0 + v N +ω scr (r). From (19) the jump of the xc potential due to the stepwise increase in the screening charge gives the derivative discontinuity in the CDFA method, We note that Eq. 28 does not require an ensemble calculation, but only the evaluation of the N -electron CDFA xc potential for two values of the screening charge and hence could be employed in practical calculations at a moderate computational cost.
In the last part of the paper, we shall compare ∆ xc from the constrained minimisation method (28) with the result for ∆ xc from (19). We shall calculate the differences in CLDA for various values of ω and investigate the limit of small ω.
Before we continue, we note that in the simple model we have constructed to predict the derivative discontinuity, using the inversion of the ensemble density (19) and with the CDFA method (28), we have restricted the freedom of the Hxc potentials, by the ansatz in (6); the restriction is that v N +ω xc (∞) = 0. Consequently, the derivative discontinuities we obtain with Eq. 19 and Eq. 28 cannot be perfect constant functions but have to vanish at r → ∞. We aim to investigate whether the resulting approximate derivative discontinuity, ∆ xc (r), as a function of r remains flat and almost equal to a constant over the region of the atom or the molecule. Finally, we want to obtain the converged value of the constant in the limit of an infinite basis set.
In order to proceed and construct the ensemble density ρ N +ω ven (r), we need the densities from two KS calculations for N and N +1 particles allowing us to then find the corresponding ensemble xc potential against which ∆ CDFA xc can be compared. We use our CLDA method to obtain the densities ρ N ven (r) and ρ N +1 ven (r), in order to control the screening densities of the constituent xc potentials. One of the integers N , N + 1 is an odd number, corresponding to an open shell system. The LDA exchange energy for open shells contains an error ("ghostexchange error" [49]) in modelling exchange with half the electrons spin-up and half spindown. In a forthcoming publication [49], we propose how to correct this error, still within LDA (not local spin density approximation). Hence, in the KS calculation for an odd number of electrons (either for N or for N + 1), we employ our method to correct for the ghostexchange error, in order to improve the accuracy of the resulting CLDA xc potential and density. Details will be published in Ref. 49.
Once we construct the ensemble density, we invert it to obtain the ensemble KS potential, v N +ω xc (r). For the inversion, we apply the method in Refs. 5, 50. The inversion method [50] requires the a priori selection of a value for the screening charge of the xc potential.
According to (26), for v N +ω xc (r) we set Q N +ω scr = N . The orbital and auxiliary basis sets are un-contracted cc-pVTZ.  In Fig. 10a, the ensemble xc potential, v N +ω xc (r), for ω = 0.1 (with Q N +ω scr = N ) is shown together with v N xc (r) and v N + xc (r), which have screening charges Q N scr = N −1 and Q N + scr = N . In Fig. 10b, the ω-dependent (29) derivative discontinuity, ∆ ω xc (r) = v N +ω xc (r) − v N xc (r), is shown for various values of ω. In the limit of small ω, ∆ ω xc (r) yields the derivative discontinuity using ensembles, ∆ ω→0 xc (r) = ∆ xc (r). The inversion method has some numerical instabilities which are exaggerated when the difference of two potentials is taken. This explains why ∆ xc (r) is not flat for small r. The distance r after which ∆ xc (r) tends to zero depends on the basis set. However, we do not propose this method as a means of computing the derivative discontinuity in practice, but rather to compare with the results of the CDFA method.
The blue line in Fig. 10b shows the CLDA prediction for the derivative discontinuity, ∆ CLDA xc (r), without an ensemble calculation. ∆ CLDA xc (r) remains almost a constant up to a distance of about 2.5 a 0 , beyond which it tends to zero. The differences ∆ ω xc (r) for decreasing ω approach the line of ∆ CLDA xc both in height and in the spatial extent over which ∆ CLDA xc and ∆ ω xc stay almost constant. The value of the constant can be obtained by inspection from Fig. 10b to be approximately 0.35 hartrees, or about 9.5eV. We can obtain the constant more accurately from the shift of the occupied single-particle energy levels between the two xc potentials v N xc (r) and v N + xc (r). For the un-contracted cc-pVTz basis used for the results in Fig. 10b we find ∆ CLDA xc = 9.48eV. See Table II. We calculated the xc derivative discontinuity with our model ∆ CDFA xc (r) (28) using as orbital and auxiliary basis sets the un-contracted cc-pVXz sets, with X=D,T,Q,5. The last row of Table II shows the value of the derivative discontinuity, ∆ CDFA xc , for each basis set.
Each column in Table II shows the eigenvalues of the occupied orbitals in the Ne atom, with the two constrained xc potentials v N xc (r) and v N + xc (r), for a specific choice of orbital and auxiliary basis sets. Each column also shows the shift of each eigenvalue ∆ i = N + i − N i . The average value of these shifts, gives ∆ CLDA xc in the specific basis. Using the un-contracted cc-pVDz orbital and auxiliary basis, the shifts of the orbital eigenvalues are almost the same within 0.03eV. In the un-contracted cc-pVTz orbital and auxiliary basis, the differences in cc-pVDz cc-pVTz cc-pVQz cc-pV5z for the orbitals i = 1s, 2s, 2p. The orbital and auxiliary basis set is un-contracted cc-pVXz, X=D,T,Q,5. The average difference ∆ i per basis set gives ∆ CLDA xc . All energies are in eV. the shifts in each energy level are smaller than 0.01eV. In the two larger basis sets, the differences between the almost constant shifts for each energy level are within 0.01eV. These results are consistent with a near perfectly constant ∆ CLDA xc over the whole spatial region where the electronic density of the Ne atom is appreciable.
For the un-contracted cc-pVTz basis, we performed another check to confirm that the shift between the two xc potentials v N + xc (r), v N xc (r) is almost a constant over a large region of space. We evaluated the overlaps of the occupied orbitals in the two potentials, φ N + i |φ N i , i =1s, 2s, 2p (triply degenerate). We found that the numerical values of all overlaps were indeed very close to one, with the overlap in the worst case differing from one by about ∼ 10 −7 .
In Fig. 11a we show the derivative discontinuity ∆ CLDA xc (r) as a function of r (Eq. 28) for orbital and auxiliary basis sets un-contracted cc-pVXz, X=D,T,Q,5. These functions have a plateau at the origin where the atom lies. The extent of the plateau increases with basis set size and the height decreases and seems to converge. To establish that the discontinuity ∆ CLDA xc (height of the plateau) indeed converges and does not vanish in the limit of infinite basis set, we plot ∆ CLDA xc against the inverse of the number of basis set elements, n bas . The behaviour is fitted well by a straight line with equation ∆ CLDA xc (n bas ) = 5.6 + 160 × (n bas ) −1 .
The extrapolation gives a nonzero derivative discontinuity of 5.6eV for the infinite basis limit.

IV. CONCLUSIONS
A common theme of popular local and semi-local density functional approximations is the imbalance of accuracy between energy density-functionals, which can be quite accurate, and the corresponding effective KS potentials, with inferior accuracy [1][2][3]. We have approached this problem from several directions [4][5][6]. In this paper, we review and expand our work on imposing physical constraints during the energy minimisation in order to yield more a accurate corresponding xc potential [8][9][10]. Specifically, we investigate the relaxation of a constraint that is computationally expensive and find that its omission leads to well behaved results, except for very small systems with only a few electrons. The constraints we impose raise the total energy minimally [8,10] but have a dramatic impact on the quality of the effective KS potential, gifting it with the correct asymptotic behaviour and enabling it to exhibit important non-analytic behaviour (derivative discontinuity) shared by the exact KS potential but elusive from the lower rungs of Jacob's ladder of DFAs where semi-local DFAs reside.