Exploiting a derivative discontinuity estimate for accurate G0W0 ionization potentials and electron affinities

The GW approximation has become an important tool for predicting charged excitations of isolated molecules and condensed systems. Its popularity can be attributed to many factors, including a favorable scaling and relatively good accuracy. In practical applications, the GW is often performed as a one-shot perturbation known as G0W0 . Unfortunately, G0W0 suffers from a strong starting point dependence and is often not as accurate as one would need. Self-consistent GW methodologies alleviate these problems but come with a marked increase in computational cost. In this manuscript, we propose the use of an estimate of the exchange-correlation derivative discontinuity to provide a remarkably good starting point for G0W0 calculations, yielding ionization potentials and electron affinities with eigenvalue self-consistent GW quality at no additional cost. We assess the quality of the resulting methodology with the GW100 benchmark set and compare its advantages over other similar methods.


I. INTRODUCTION
Density functional theory (DFT) [1,2] is now recognized as one of the more important tools for determining the electronic structure of molecules and extended systems.It is well-known, however, that common density functional approximations (DFAs) possess some shortcomings attributable to the limitations and constraints of such approximations [3].
For example, it has been known for decades that the Kohn-Sham eigenvalue of the highest occupied molecular orbital (HOMO) is the negative of the vertical ionization potential (VIP), −I [4][5][6].Common DFAs, however, yield HOMO eigenvalues that deviate substantially from −I due to the incorrect continuous behavior of the energy as a function of the number of electrons N , instead of the correct piecewise linear relation [4,5,7].The eigenvalue of the lowest unoccupied molecular orbital (LUMO) also shows important deviations from the negative of the vertical electron affinity (VEA), −A, due to the deviation from the linearity condition [8,9].These undesirable behaviors have been identified as fractional charge errors which translate into either delocalization (convex deviations) or localization (concave deviations) errors [9][10][11].It has been found that DFAs based on the local density approximation (LDA), generalized gradient approximation (GGA), or meta-GGA show a rather large convex curve for the energy as a function of N , while Hartree-Fock (HF) yields a concave deviation, a fact that has been exploited in the development of several DFAs including optimally-tuned range-separated hybrids [12][13][14][15][16][17][18][19].Unfortunately, the errors in approximating VIP and VEA through the frontier KS eigenvalues do not cancel out as most of the time VIP is underestimated and VEA is overestimated, usually resulting in the significant underestimation of the fundamental gap E g = I − A in solids (bandgap) and molecules (hardness [20]).This has been known as the bandgap problem (of DFT) and affects the accuracy of several response properties and chemical reactivity descriptors [21] that depend on E g .KS eigenvalues from orbitals below HOMO and above LUMO have also been used as approximate addition and removal energies [15,17,[22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40].There is no rigorous theoretical connection between the KS eigenvalues and the quasiparticle energies, but there are practical benefits of a DFA that yields eigenvalues in reasonably close agreement with quasiparticle energies.On the one hand, these KS eigenvalues could be directly used to try to understand and predict the properties of novel materials in a computationally efficient way.On the other hand, both KS eigenvalues and orbitals can be used as input for other many-body methods known to solve, or at least alleviate, the bandgap problem.Prominent among these methods is the GW approximation (GWA) to the self-energy [41][42][43].The GWA is usually performed perturbatively as a one-shot correction, G 0 W 0 [44,45], which can be significantly affected by the quality of the starting orbitals and eigenvalues [46,47].Fully self-consistent GW (scGW ) [48,49] calculations provide one way to eliminate the starting-point dependence at a steep computational cost.Lower-scaling partial self-consistent GW approaches, like quasiparticle self-consistent GW (qsGW ) [50,51] and eigenvalue self-consistent GW (evGW ) [52,53], reduce the starting-point dependence but still have much higher computational costs than G 0 W 0 .It is therefore highly desirable to find a scheme to achieve accurate G 0 W 0 @DFA with reasonable computational costs and robust accuracy.
This route is similar in spirit to the ∆GW method from Vlček and cols.[64], but here we use a shift before the GW calculation as the NCAP DFAs exploit the fact that the derivative discontinuity (DD) implied in some DFAs as N crosses integer values can be directly obtained from the non-zero constant needed to make the exchange-correlation (XC) potential, v xc , vanish at infinity.We will show that NCAP eigenvalues corrected by the DD shift (NCAP-DD) provide a remarkably good starting point for G 0 W 0 calculations, yielding quasiparticle energies on par with the more expensive partially self-consistent evGW @GGA and evGW @mGGA.We will begin our presentation by introducing the theory behind the development of the NCAP functionals, followed by a basic introduction to the GWA.Then, we will use the GW100 [65] benchmark set to showcase the accuracy that can be achieved by the G 0 W 0 @NCAP-DD method.

II. THEORY A. Derivative discontinuity
The exact behavior of the energy, E, as a function of the number of electrons is given by lines connecting the integer values of N [4,5,7], generally with different slopes on the electron-deficient and electron-abundant sides of N .This piecewise linear behavior implies that for the exact functional, where the superscripts indicate that the derivative is obtained from the electron abundant (+) or deficient side (−).Equation (1), in turn, implies that the functional derivative of E xc with respect to the electron number density, n(r), is also discontinuous as N crosses an integer value [5,6,66,67], i.e.
It is customary to assume that the DD ∆ xc is constant [68] and that with ε + H as the HOMO eigenvalue determined with v + xc (r), and ε − L as the corresponding LUMO eigenvalue determined with v − xc (r).Furthermore, the ionization potential theorem states that [6] ε where ε − H is the HOMO eigenvalue determined with v − xc (r).Several relations can be obtained from the previous set of equations, including the confirmation that the KS gap obtained from XC DFAs that do not incorporate the DD of the XC potential will underestimate the fundamental gap Moreover, the shifts needed to pair the Kohn-Sham eigenvalues to the negative VIP and VEA can be related to the full DD via The takeaway piece of this analysis is that the full DD can be decomposed into separate shifts for the occupied and unoccupied spaces, a property exploited by Carmona-Espíndola and cols.[62,63] B. The NCAP exchange functional The NCAP exchange enhancement factor [61] given in terms of the reduced density gradient s(r) = |∇n(r)|/2(3π 2 ) 1/3 n 4/3 (r), was purposely designed to generate an exchange potential that goes asymptotically to a constant plus a term that decays as −c/r.The values of the constants appearing in Equation (8) were obtained by satisfaction of different constraints or norms and are given in Table I It can be shown that for a density whose asymptotic behavior is given by the NCAP exchange potential will tend to The ionization potential theorem requires that the potential decays to zero, so one needs to add the constant to the NCAP potential and the HOMO eigenvalue obtained directly from the self-consistent field KS calculation, that is Equation ( 11) can be written in terms of ε NCAP H to obtain a quadratic equation in v DD x whose solutions can be related to ∆ − xc and ∆ + xc , and can be used to obtain the approximations We will use v DD− x and v DD+ x to shift the entire occupied and unoccupied eigenvalues to yield the G 0 W 0 @NCAP-DD approach.

C. GW Approximation
The GWA to the self-energy utilizes the one-body Green's function, G(r 1 , r 2 , ω), to describe the spectral properties of materials and molecules as its poles in the complex plane are the true addition and removal energies of the system [41-43, 69, 70].The success of the GWA depends on the ability to compute a good approximation to G without an explicit need for the many-body wave function.Hedin's closed-set of equations (Figure 1) is, in principle, an exact route that achieves precisely this goal.Different levels of approximation can be introduced into Hedin's equations, but the most natural one is neglecting all corrections to the bare vertex Γ 0 (1; 2, 3) = δ(1, 2)δ(1, 3) in the calculation of the polarizability P (1, 2) and the self-energy Σ(1, 2).Doing so yields (in simplified notation) Σ = iGW , which gives the name to this particular approximation.
Formally, Hedin's equations must be solved self-consistently.Approaches that do not do so introduce another level of approximation.In particular, the G 0 W 0 method goes once through the loop shown in Figure 1, entering with some mean-field (here assumed to be the KS-DFT solution) orbitals {ϕ KS } and corresponding eigenvalues {ε KS } to compute the "zeroth-order" Green's function G 0 and screened Coulomb interaction W 0 .The method ends once the desired quasiparticle energy ε G 0 W 0 p is obtained by solving the non-linear quasiparticle equation given by Note that only the diagonal elements of the self-energy are important for the G 0 W 0 method, as the orbitals are not updated.The quality of the G 0 W 0 quasiparticle energies depends significantly on the starting KS solution due to the lack of self-consistency.A partial self-consistent approach that only enforces the self-consistency of the quasiparticle energies, known as evGW , alleviates the starting-point dependence and usually improves the quality of the quasiparticle energies [52,53].In evGW , one loops around Hedin's equations several times without updating the orbitals.As in G 0 W 0 , only the diagonal of the self-energy is needed.The updated addition and removal energies are used to compute an improved approximation for G in each loop and the method ends once the quasiparticle energies stop changing.The evGW method can be several times more expensive than the one- shot G 0 W 0 approximation even when only one additional loop through Hedin's equations is performed.This computational burden follows from the need to solve, in principle, the quasiparticle equations for all orbitals instead of solving them for only a handful of states 1 .It is therefore customary to utilize the so-called "scissor-shift operator" method, wherein a limited number of quasiparticle energies above and below the Fermi level are updated.In contrast, the remaining quasiparticle energies are rigidly shifted using the information of their respective space (occupied/unoccupied).Unfortunately, the validity of this shift has been recently assessed to be very limited [71].
There are several ways to implement the solution of Hedin's equations, each with its own advantages and disadvantages (including numerical stability and computational cost).
We have recently implemented a scalable approach [72], based on the contour deformation (CD) approach [73], in the open-source quantum chemistry package NWChem [74,75].
Our implementation uses the density fitting technique [76,77] (also known as resolution-ofthe-identity technique) to avoid computing and manipulating four-center electron repulsion integrals, and the MINRES solver [78] to avoid expensive, and numerically sensitive, ma-trix inversions.In the molecular CD-GW method, it is customary to use the Adler-Wiser representation of the irreducible polarizability [79, 80] to obtain the dynamical dielectric function and, from its inverse, the screened Coulomb interaction In the preceding equations, σ labels the spin, i labels occupied orbitals, and a labels unoccupied ones.Note that W 0 , Equation ( 23), was conveniently split into an exchange part (first term) and a correlation part (second term).The convolution of G 0 with the exchange part of W 0 is performed analytically to yield exchange-like potential The convolution with the correlation part is performed semi-analytically following the integration paths shown in Figure 2, wherein the contours are taken to infinity and chosen to exclude the poles of W 0 .At infinitely large ω ′ , the product G 0 (r, r ′ , ω + ω ′ )W 0 (r, r ′ , ω ′ ) vanishes, so the real-frequency integral can be computed by subtracting the imaginaryfrequency integral from the contour integrals.The integration along the imaginary axis is performed numerically, whereas the contour integrals are replaced by a sum over the residues enclosed within the contours.More details about the CD-GW approach implemented in NWChem can be found in Reference 72.

III. COMPUTATIONAL METHODS
The G 0 W 0 @NCAP-DD approach was assessed using the GW100 [65] dataset to compare vertical ionization potentials and electron affinities for a series of small molecules.All calculations were performed using the current release of NWChem (version 7.2.2)2 .

FIG. 2. Schematic representation of the contour deformation technique to obtain Σ c . The contours
are chosen to contain the poles of G 0 but not the poles of W 0 .The integral over the real frequencies ω ′ is obtained using a numerical integration over the imaginary axis and the residue theorem.
The polarization-consistent pcSseg-3 basis set [81] was selected due to its balanced description between valence addition/removal energies and core removal energies within the GWA [82].Since there is no optimal fitting basis to pair with pcSseg-3, we automatically generated it following the procedure detailed in Reference 83.We used the def2-QZVP/def2universal-JKFIT basis set [84,85] combination for all those elements for which the pcSseg-3 basis set is not defined.
For CD-GW , the integral over the imaginary frequencies was discretized using a Gauss-Legendre quadrature with 200 points and the imaginary infinitesimal was set to 0.001.The quasiparticle energies were obtained by solving the corresponding non-linear quasiparticle equations using Newton's fixed-point approach.To reduce numerical noise as much as possible (especially for evGW calculations), we used the SIMINT electron repulsion integral package [86] with no prescreening of the primitive products, a Schwarz-screening threshold of 10 −16 , tight SCF convergence of 10 −9 Hartrees, and a xfine numerical integration grid 3 .
Following Section II, the calculations labeled as NCAP-DD used a shift based on the 3 See NWChem documentation for more details about these settings.
estimated NCAP derivative discontinuity given by Equation ( 13), The shift was removed while solving the quasiparticle equations, i.e.
Partially self-consistent evGW calculations were performed by computing the screened Coulomb interaction with the so-called spectral decomposition approach (see Reference 72 for more details).This method scales as O(N 6 ) but, many times, is more stable than the CD-GW method.

A. Vertical Ionization Potentials
We will start by comparing VIPs obtained with G 0 W 0 and evGW using three different starting points: NCAP-DD, NCAPR-DD, and PBE. Figure 3 shows that all but one evGW VIPs differ less than 0.1 eV, most of which lie between ±0.05 eV (shaded area).The largest difference, 0.11 eV, is observed for the VIP of NaCl predicted by evGW @NCAP-DD and evGW @NCAPR-DD (-0.058 eV and +0.053 eV deviations relative to evGW @PBE, respectively).The order of these deviations is in line with those reported in the literature for VIPs when comparing "pure" DFAs only [87].DFAs with some admixture of exact exchange show slightly different behaviors that can push the overall starting-point dependence to values as high as 0.5 eV [88].
Having established the equivalency of the evGW ending points for the three starting points, we now turn to compare the VIP shift going from G 0 W 0 to self-consistent evGW for the same starting points.Figure 4 shows that G 0 W 0 @PBE VIPs underestimate the evGW @PBE values by as much as 2 eV.In contrast, G 0 W 0 @NCAP-DD and G 0 W 0 NCAPR-DD VIPs are much closer to the corresponding evGW values.The modest revisions introduced in NCAPR [63] translate into a non-negligible improvement for the G 0 W 0 VIPs as the mean absolute deviation diminishes from 0.149 eV for G 0 W 0 @NCAP-DD to just 0.061 eV for G 0 W 0 @NCAPR-DD.It is worth mentioning, however, that the average improvement comes with a slight degradation in the dispersion of the errors, with the standard deviation increasing from 0.059 eV (NCAP-DD) to 0.088 eV (NCAPR-DD).
Based on these results, it is possible to conclude that G 0 W 0 @NCAP-DD and, especially, G 0 W 0 @NCAPR-DD VIPs are no different than those obtained with evGW starting from other KS solutions using "pure" DFAs.

B. Vertical Electron Affinities
A similar analysis can be performed for the electron affinities of the molecules contained in the GW100 dataset.Figure 5 shows that the evGW starting-point dependence is stronger for VEAs than for VIPs, as some values deviate as much as 1 eV.These rather large differences might be due to the qualitatively different behavior of the exchange potentials of the NCAP DFAs, which decay as −c/r, from other standard GGA potentials which decay much faster.
The shapes of the orbitals in the unoccupied space, including the LUMO, seem to be more affected by the nearly correct asymptotic decay of the NCAP potential.However, further explorations are needed to conclusively determine the origin of the VEA deviations.We note that a closely related DFA with a correct asymptotic potential exhibits similar effects on the unoccupied part of the spectrum, although the effect was attributed to an odd behavior in its potential at intermediate distances [89].
Unfortunately, there is a lack of high-level reference VEAs values for all the molecules contained in the GW100 dataset.Comparison to the experiment is also not straightforward as, most of the time, the experiments measure adiabatic electron affinities (with rather large error bars).A further complication is that the majority of VEAs in the GW100 dataset are negative, meaning that the corresponding anions are unbound and, consequently, that experimental data is difficult to obtain.Because of this, we will set aside the question of which evGW @DFA limit is better for VEAs.
We can, nevertheless, assess the quality of the G 0 W 0 VEAs relative to the corresponding evGW values.Figure 6 shows that G 0 W 0 @NCAP and G 0 W 0 @NCAP-DD VEAs are very close to the respective evGW results.The mean absolute deviations are just 0.049 eV and 0.021 eV for G 0 W 0 @NCAP and G 0 W 0 @NCAP-DD, respectively.Different from the VIP case, G 0 W 0 @PBE VEAs are relatively close to the corresponding evGW @PBE VEAs but systematically overestimated.It is well-known that updating W in (partial) self-consistent GW calculations often results in underscreening since quasiparticle energy differences appear in the denominator of W [90,91].This underscreening leads to photoelectron spectra that appear elongated, with overestimated ionization potentials and underestimated electron affinities.The fact that G 0 W 0 @NCAP-DD overestimates the VIPs and underestimates the VEAs relative to evGW @NCAP-DD might be problematic because this will translate into overshooting an already overestimated fundamental gap.This effect is almost absent in G 0 W 0 @NCAPR-DD.

C. Comparison to other methods
One of the main advantages of using the NCAP-DD solution as a starting point is that Equation 13 is valid for any molecular system.In contrast, the range separation and exact exchange mixing parameters must be changed from system to system when using optimally tuned range-separated hybrids (OTRSH), often needing several self-consistent calculations to achieve the desired outcome.Furthermore, the NCAP ground-state solution can be obtained very efficiently by using density fitting to approximate the Coulomb potential.
This is very important because the ground-state solution is often the bottleneck of the whole G 0 W 0 @DFA approach when only a handful of frontier quasiparticle energies are of interest.As a consequence, G 0 W 0 @NCAP-DD calculations can be significantly faster than OTRSH and G 0 W 0 @OTRSH ones.We also note the existence of a one-shot tuning, for global and range-separated hybrids [92], which should be more computationally efficient than the optimal tuning in OTRSHs.This novel one-shot approach, however, still needs tuning parameters for each individual system.
The same non-universality is found for the screened Koopman's compliant (KC) DFAs, where the screening parameter should be obtained for each orbital.In the original implementation, each parameter was obtained from separate calculations on the neutral and singly-ionized systems [30].Approximate screening parameters can also be obtained using perturbation theory [93], reducing the number of calculations needed.
Another related method is the ∆GW approach from Vlček and cols.[64].However, ∆GW needs a G 0 W 0 calculation to estimate the size of the occupied and unoccupied shifts.In contrast, NCAP-DD provides the shifts before entering the GW calculation itself.However, we note that both approaches need only one iteration through Hedin's cycle to achieve the desired outcome.Another key point is that the ∆GW method offers ionization potentials Comparison of the performance of several electronic structure methods for the prediction of vertical ionization potentials.The data for OTRSH and G 0 W 0 @OTRSH was obtained from Reference [60], for KI, pKIPZ, and KIPZ Koopman's-compliant functionals from Reference [58], and for CCSD(T) from Reference [94].
that deviate several tenths of eV from the evGW 0 target.We haven't tested our approach in the evGW 0 framework, but we expect to achieve the same fidelity as for evGW , namely less than 0.1 eV difference between the G 0 W 0 and evGW 0 .
For completeness, Figure 7 compares the accuracy of the VIPs obtained with an OTRSH DFA [60], several KC DFAs [58], and several flavors of GW calculations, relative to a CCSD(T) reference [94].The methods are ordered by increasing mean absolute error.Is not surprising that several system-dependent approaches are the best performing of the set, but the evGW approaches, including the G 0 W 0 @NCAP-DD family, perform rather well.
It is unclear what performance would the eVGW @OTRSH, but the VIPs will likely be overestimated due to the underscreening attained while iterating W .
Finally, we will only mention that the renormalized singles approach goes beyond the G 0 W 0 @NCAP-DD family (and the related ∆GW ), as the underlying KS orbitals are also updated before the full GW calculation.

V. CONCLUSIONS AND FUTURE OUTLOOK
We have shown that an estimate of the derivative discontinuity (DD) present in the NCAP family of functionals can be used, at no computational cost, to provide a remarkably good starting point for G 0 W 0 calculations.The DD estimate is split into two shifts, one for the occupied part of the spectrum and one for the unoccupied part of the spectrum, according to the procedure detailed by Carmona-Espíndola and cols. in References [62] and [63].The shifted eigenvalues are then used to compute an estimate of the Green's function G and the screened Coulomb interaction W that are in agreement with those obtained in partially self-consistent evGW methods.This procedure leads to a one-shot G 0 W 0 approach that yields vertical ionization potentials and electron affinities with evGW quality, saving valuable computational time.
The best performance overall was obtained using the G 0 W 0 NCAPR-DD approach, showing mean absolute deviations, relative to evGW @NCAPR-DD, of 0.06 eV and 0.02 eV for ionization potentials and electron affinities, respectively.The G 0 W 0 @NCAP-DD performance is not far behind, although it can aggravate the fundamental gap overestimation seen in evGW approaches.
Vertical electron affinities show a strong starting-point dependence, even with the selfconsistent evGW method.Explorations to obtain the root of these differences are being conducted.
The G 0 W 0 @NCAP-DD approach is naturally well-suited to provide the quasiparticle energies and screened Coulomb interaction needed to obtain neutral excitations via the Bethe-Salpeter equation formalism.

VII. SUPPLEMENTAL MATERIAL
Tables with names of the systems contained in the GW100 benchmark dataset as well as individual ionization potentials and electron affinities.

FIG. 1 .
FIG. 1. Schematic representation of Hedin's equations, showing the omission of vertex effects, Γ in the GW approximation.The entry point is the calculation of G with some mean-field approximation.The quasiparticle energies can be obtained using the exchange-correlation self-energy effective potential.Alternatively, a better approximation to G can be obtained, and the circuit iterated until self-consistency.

FIG. 3 .
FIG.3.Deviation of the self-consistent evGW vertical ionization potentials, in eV, relative to the evGW @PBE value.The x-axis is ordered according to the list of molecules shown in the Supplementary Material.

FIG. 4 .
FIG.4.Vertical ionization difference, in eV, between G 0 W 0 @DFA and the corresponding self-consistent evGW @DFA one.The left plot shows individual differences with the x-axis ordered according to the list of molecules shown in the Supplementary Material.The violin plots on the right show the aggregated statistics with the mean appearing as a horizontal line and the median as a white dot.

FIG. 5 .
FIG.5.Deviation of the self-consistent evGW vertical electron affinities, in eV, relative to the evGW @PBE value.The x-axis is ordered according to the list of molecules shown in the Supplementary Material.

FIG. 6 .
FIG.6.Vertical electron affinity difference, in eV, between G 0 W 0 @DFA and the corresponding self-consistent evGW @DFA one.The left plot shows individual differences with the x-axis ordered according to the list of molecules shown in the Supplementary Material.The violin plots on the right show the aggregated statistics with the mean appearing as a horizontal line and the median as a white dot.
This work was supported by the Center for Scalable Predictive methods for Excitations and Correlated phenomena (SPEC) which is funded as part of the Computational Chemical Sciences (CCS) program, under FWP 70942, by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences and Biosciences at Pacific Northwest National Laboratory (PNNL).PNNL is a multi-program national laboratory operated by Battelle Memorial Institute for the United States Depart-ment of Energy under DOE contract number DE-AC05-76RL01830.

TABLE I .
Constants defining the NCAP and NCAPR exchange enhancement factors.