Combining the Δ-Self-Consistent-Field and GW Methods for Predicting Core Electron Binding Energies in Periodic Solids

For the computational prediction of core electron binding energies in solids, two distinct kinds of modeling strategies have been pursued: the Δ-Self-Consistent-Field method based on density functional theory (DFT), and the GW method. In this study, we examine the formal relationship between these two approaches and establish a link between them. The link arises from the equivalence, in DFT, between the total energy difference result for the first ionization energy, and the eigenvalue of the highest occupied state, in the limit of infinite supercell size. This link allows us to introduce a new formalism, which highlights how in DFT—even if the total energy difference method is used to calculate core electron binding energies—the accuracy of the results still implicitly depends on the accuracy of the eigenvalue at the valence band maximum in insulators, or at the Fermi level in metals. We examine whether incorporating a quasiparticle correction for this eigenvalue from GW theory improves the accuracy of the calculated core electron binding energies, and find that the inclusion of vertex corrections is required for achieving quantitative agreement with experiment.


INTRODUCTION
The energy required to remove a core electron from a particular atom depends on the atom's chemical environment. In core level X-ray Photoelectron Spectroscopy (XPS), this dependence can be exploited to identify the chemical environments that are present in the sample. XPS is particularly well suited for the analysis of complex surfaces, and it plays an important role in the study of heterogeneous catalysis, 1−4 corrosion, 5−7 environmental degradation, 8−10 or the manufacture of surface coatings. 11−14 However, the interpretation of XPS spectra is challenging, which has motivated the development of computational techniques for calculating core electron binding energies from first principles. 15−30 For the prediction of absolute core electron binding energies in periodic solids, two kinds of methods have emerged. In the total energy difference method based on density functional theory (DFT), also known as the Δ-Self-Consistent-Field (ΔSCF) method, the core electron binding energy is calculated as the difference between total energies from two separate calculations: one for the system with a core hole, and one for the system without it. 31 −33 In contrast, in the GW method, the core electron binding energy is calculated as the GW eigenvalue of the relevant core eigenstate. 34,35 In a typical GW calculation, ground state orbitals and orbital eigenvalues are first obtained using DFT, and next, GW corrections to the eigenvalues are obtained by applying the GW method in a "one-shot" (G 0 W 0 ), or partly self-consistent manner. Direct GW calculations of core electron binding energies involve some additional complications, when compared to GW calculations of valence states. Issues such as the treatment of the frequency-dependent self-energy, basis set convergence and extrapolation, starting point dependence, and the role of (partial) self-consistency have been discussed extensively in recent works. [25][26][27]36,37 In brief, very promising results have recently been obtained for molecular systems (mean absolute error <0.3 eV) [ref 36], whereas somewhat larger mean absolute errors (0.53 and 0.57 eV in references 34 and 35, respectively) have been observed in the few preliminary studies of periodic solids published thus far.
In this work, we examine the formal relationship between the ΔSCF and GW methods and combine the two approaches by establishing the link between total energy differences and energy eigenvalues. In addition, we examine how this insight can be exploited to improve the accuracy of calculated binding energies.

ΔSCF METHOD FOR PERIODIC SOLIDS
When calculating or measuring core electron binding energies in solids, a well-defined point of reference must be used. In experimental XPS, the sample Fermi level is typically used as the zero of the energy scale. However, as discussed in ref 31, this choice is not well suited for theoretical calculations of core electron binding energies in insulators, as the position of the Fermi level within the band gap is not in general known a priori, and it depends strongly on extrinsic factors, such as the concentration of defects or impurities in the sample. Therefore, in recent computational studies, the energy of the highest occupied state, i.e., the Fermi level in metals and the valence band maximum (VBM) in insulators, has been used as the point of reference instead. 31,34,35 For total energy difference methods, this means that the core electron binding energy is defined as the difference between two total energy differences: the ΔSCF result for the core where E B is the calculated core electron binding energy relative to the VBM in insulators or the Fermi level in metals, E N,ground is the ground state total energy, E N−1,ground is the total energy of the system with one electron removed from the highest occupied state, and E N−1,ch is the total energy of the system with a core hole. This formalism was used to calculate absolute core electron binding energies in solids in ref 31. It was shown that core electron binding energies from periodic ΔSCF calculations based on DFT with the SCAN functional 38 were in good agreement with experimental values. In particular, the mean absolute error was just 0.24 eV for a small test set of 15 core electron binding energies. However, in some cases, significantly larger errors were observed, e.g., the C 1s binding energy in diamond was overestimated by 0.39 eV, and the Be 1s and O 1s binding energies in BeO were in error by 0.79 and 1.16 eV, respectively. In ref 31, it was speculated that these errors arise from the inability of DFT to accurately predict the position of the VBM in wide band gap insulators.
2.1. The VBM Energy in Density Functional Theory. In this study, we investigate this matter further. At first, from eq 1, it would seem that the VBM energy in fact never needs to be explicitly calculated for obtaining the core electron binding energy. However, the second term in the brackets before simplification does correspond to a total energy difference calculation of the VBM energy. The relationship between the term (E N−1,ground − E N,ground ) and the VBM Kohn−Sham eigenvalue in DFT, ϵ max , has been previously discussed, e.g., in refs 39 and 40. In particular, as explained in ref 39, the energy difference between a pure material and a material with a single hole becomes equal to the VBM Kohn−Sham eigenvalue in the limit of a dilute hole gas. In real calculations using finite supercells, however, this energy difference only slowly converges to the infinite limit as the system size is increased. Formally, where n is the number of atoms per supercell, and ϵ max is the energy of the highest occupied state, i.e., VBM eigenvalue in insulators, or the eigenvalue at the Fermi level in metals. The preceding discussion pertains to DFT with real (approximate) exchange-correlation functionals. In exact DFT, the equality in eq 2 holds at any supercell size. Equation 2 shows that in solids, at the limit of infinite supercell size, VBM energies calculated as total energy differences must have exactly the same shortcomings as Kohn−Sham eigenvalues.

Alternative Formalism for Periodic ΔSCF Calculations of Core Electron Binding Energies. Equation 2
allows us to write an alternative expression for the core electron binding energy, by replacing the term ( In the limit of infinite supercell size, eq 1 and eq 3 should give the same result, but for finite supercells the calculated core electron binding energies differ. A numerical verification of eqs 2 and 3 is presented next.

Numerical Verification of Equation 2.
We have calculated IE ΔSCF , defined as E N−1,ground (n) − E N,ground (n), and IE ϵ , defined as −ϵ max for all of the 10 solids�Li, Be, Na, Mg, graphite, BeO, hex-BN, diamond, β-SiC, and Si�and all of the supercells considered in ref 31, using DFT with both the SCAN and the PBE functionals. 38,41 As an example, the results for diamond obtained using the SCAN functional are shown in Figure 1. In Figure 1a, IE ΔSCF − IE ϵ is plotted against the number of atoms per supercell (n), and in Figure 1b, the same quantity is plotted against the inverse cube root of n, as is done when extrapolating core electron binding energies to the infinite supercell limit. Figure 1a,b shows that IE ΔSCF − IE ϵ indeed slowly approaches zero as the size of the supercell increases. Similar behavior is also observed for the other materials, using both PBE and SCAN�the detailed results are provided in the SI.

Numerical Verification of Equation 3.
Next, the core electron binding energies calculated using eq 1 and eq 3 are compared in Figure 2. In Figure 2, calculated core electron binding energies in one insulator, diamond, and one metal, Na, are shown as a function of supercell size. In each plot, the infinite supercell limit lies at the y-axis intercept. Figure 2a shows calculated C 1s binding energies in diamond from eq 1 and eq 3. The extrapolated values, 284.43 and 284.36 eV, respectively, differ by 0.07 eV�this is attributed to uncertainties in extrapolation and errors caused by finite kpoint sampling. While not negligible, this difference is less than half of the average error in the calculated binding energies and of the same magnitude as the precision with which experimental binding energies are typically reported. In Figure   Figure 1. Numerical validation of eq 2 for diamond. In panel (a), the difference between the first ionization energy calculated using the total energy difference method and the negative eigenvalue of the highest occupied state is plotted against the number of atoms in the supercell. As the size of the supercell increases, the difference slowly tends toward zero. In panel (b), the same quantity is plotted against the inverse cube root of the number of atoms per supercell. 2b, the calculated Na 1s binding energies in Na metal from the two equations are compared. In this case, and similarly for other metals, for sufficiently large supercells, both equations yield core electron binding energies that are converged to the limiting value. For the Na 1s binding energy in Na metal, the limiting values from eq 1 and eq 3 differ by less than 0.01 eV.
Further numerical verification of eq 3 is provided in Table 1, where a comparison of the extrapolated results from eq 1 and eq 3 is provided for all of the 15 core levels considered in ref 31. The same calculations have been performed using both the PBE and SCAN exchange-correlation functionals. In summary, the two equations yield very similar results. For both SCAN and PBE, the root mean squared deviation between the calculated binding energies from the two equations is just 0.07 eV.

Localized vs Delocalized Hole States.
It is important to emphasize that an identity similar to eq 2 does not hold for the core electrons, i.e., ,ground core (4) provided that the core hole is properly localized in the calculation of E N−1,ch . This is numerically illustrated in Figure  3. This fundamental difference arises due to the fact that in valence ionization an electron is removed from a delocalized state, and as the size of the simulation cell increases, the change in the local potential experienced by all the remaining electrons slowly tends toward zero. In contrast, in core ionization, an electron is removed from a localized state, and in the vicinity of the atom with a core hole, the remaining electrons experience a large change in local potential regardless of the size of the supercell. Here, the terms "localized" and "delocalized" refer to the spatial distribution of a Kohn−Sham state relative to the simulation cell (that in general contains many unit cells of the solid). A localized core hole is centered around exactly one atom, regardless the size of the simulation cell, thus breaking Figure 2. A comparison of calculated core electron binding energies from eq 1 and eq 3 for the C 1s level in diamond (panel a) and the Na 1s core level in sodium metal (panel b). For finite supercells, eqs 1 and 3 can give different results. However, at the limit of infinite supercell size, the calculated binding energies from eq 1 and eq 3 converge to the same limiting value. The results are shown for two sets of calculations, one using the exchange-correlation functional PBE, and the other using the exchangecorrelation functional SCAN. All energies are given in eV. Figure 3. Difference between the calculated core electron binding energy from a total energy difference calculation and the negative eigenvalue of the core orbital, as a function of supercell size. Results for the C 1s core level in diamond are shown in panel (a), and results for the Na 1s core level in sodium metal are shown in panel (b). In contrast to the behavior observed for the first ionization energy (Figure 1), for core electron binding energies, the difference does not approach zero with increasing supercell size. This is due to the localized nature of the core hole, as opposed to the delocalized nature of the hole in the valence band. the translational symmetry in the system. In contrast, a hole in a delocalized state is evenly distributed over all symmetryequivalent atoms in the simulation cell.

SIGNIFICANCE OF EQUATION 3
Conceptually, eq 3 highlights that the accuracy of core electron binding energies from periodic ΔSCF calculations depends on the accuracy of ϵ max , i.e., the DFT eigenvalue of the highest occupied state. However, DFT is widely known to underestimate band gaps in solids, and more advanced theories such as the GW approximation yield significant corrections to both the VBM and CBM (conduction band minimum) energies predicted by DFT. It is therefore reasonable to consider whether it is possible to improve the accuracy of calculated core electron binding energies in insulating solids by adding a quasiparticle correction to ϵ max in eq 3.
In other words, provided that a consistent point of reference can be established, one could try to calculate (E N−1,ch − E N,ground ) using a method that is optimal for predicting core electron binding energies, and ϵ max using a method that is optimal for modeling the removal of valence electrons, and combine the two values to obtain a "theoretical best estimate" core electron binding energy referenced to the VBM (or E F in metals).

COMBINING THE ΔSCF AND GW APPROACHES
In this work, we attempt to combine core electron binding energies calculated using the ΔSCF method with VBM energies calculated using the G 0 W 0 approach. In particular, we have performed the following calculations.
(i) We have calculated (E N−1,ch − E N,ground ), as well as ε max DFT , for all of the materials, core levels, and supercells considered in ref 31, using DFT with two different functionals: PBE and SCAN. These calculations have been performed in the allelectron electronic structure code FHI-aims. 42  = + + . The total energies E N−1,ch and E N,ground , as well as ϵ max , have all been calculated in FHI-aims, using the same structures, physical settings (functional and treatment of relativistic effects), and numerical settings (basis sets, integration grids, etc.).
(iv) We have also attempted to combine a G 0 W 0 correction with the core electron binding energies calculated using the SCAN functional. For technical reasons, and due to the limited current knowledge about the performance of DFT with the SCAN functional as a starting point for perturbative GW calculations, we have not at present calculated G 0 W 0 corrections to the VBM (or Fermi level) eigenvalues from SCAN. Instead, we have chosen to test a strategy where the G 0 W 0 @PBE correction is combined with core electron binding energies from ΔSCF calculated using the SCAN functional. This requires an additional step, because the correction is defined relative to the PBE eigenvalue of the highest occupied state, not the SCAN eigenvalue. Therefore, we also have to correct for the difference between ϵ max PBE and ϵ max SCAN , and the c o r r e c t e d b i n d i n g e n e r g i e s a r e o b t a i n e d a s is evaluated at the optimized density from PBE. In order to assess the severity of this approximation, we have compared ΔE PBE@SCAN with ΔE SCAN@PBE , for each of the materials considered, i.e., the

Journal of Chemical Theory and Computation
pubs.acs.org/JCTC Article differences between the SCAN and PBE eigenvalues at the relaxed density from either functional. We have found that ΔE PBE@SCAN ≈ −ΔE SCAN@PBE in all cases, with all differences in the absolute values being less than 0.02 eV. Thus, we obtain (i) uncorrected core electron binding energies from eq 3 using PBE and SCAN: E B PBE and E B SCAN , and (iii, iv) core electron binding energies that have been recalibrated to the position of the highest occupied state predicted by the G 0 W 0 @PBE method: E E B PBE, G 0 W 0 @PBE , and E E B SCAN, G 0 W 0 @PBE . The initial results obtained using this approach are disappointing. In fact, as shown in Tables 2  and 3, including the correction for ϵ max from G 0 W 0 theory worsens the agreement with experiment considerably.
For PBE, the mean absolute error (MAE) increases from 0.66 to 1.45 eV, and for SCAN, the MAE increases from 0.24 to 0.56 eV. In particular, we find that if the G 0 W 0 correction for the highest occupied state is included, the calculated binding energies are too low, as compared to experiment, in all cases. This means that the mean signed errors (MSE) are equal in magnitude to the mean abolute errors:

Effect of Vertex Corrections in GW.
In ref 43, it was argued that while the G 0 W 0 method is highly accurate for band gaps in periodic solids, it relies partly on error cancellation, and that the absolute band energies predicted by G 0 W 0 are considerably less accurate. An improved methodology, termed G 0 W 0 Γ, was proposed, in which so-called vertex corrections derived from the renormalized adiabatic local density approximation (rALDA) kernel are included. It was shown that, as compared to G 0 W 0 , the band gaps predicted by G 0 W 0 Γ are largely unchanged, whereas the absolute positions of the band edges are shifted upward by approximately 0.6 eV in the examples considered.
We have examined whether using the G 0 W 0 Γ@PBE correction to the energy of the highest occupied state, instead of the G 0 W 0 @PBE correction, improves the results. The 0 @PBE results is just 0.19 eV, which is smaller than the MAE of the uncorrected binding energies. Overall, the G 0 W 0 Γ correction improves the accuracy of the calculated binding energies in nonmetals: the MAE of the corrected binding energies is 0.19 eV, as compared to 0.35 eV for the pure ΔSCF results with SCAN. In particular, the G 0 W 0 Γ correction significantly improves the results for the difficult cases of diamond and BeO�the errors in the C 1s, Be 1s, and O 1s binding energies are reduced to 0.04 eV, −0.46 eV, and −0.11 eV, respectively, compared to 0.32, 0.78, and 1.13 eV for the E B SCAN values. In the metallic systems considered in this work, the accuracy of the original ΔSCF results with SCAN is already very high: MAE = 0.08 eV; the G 0 W 0 Γ correction makes the agreement somewhat worse, although the MAE remains relatively small at 0.20 eV.

CONCLUSIONS
In summary, this study establishes a direct link between the two fundamentally different strategies that can be employed for calculating core electron binding energies: total energy difference methods, and eigenvalue methods. Formally, this is expressed as the equivalence of eqs 1 and 3 in the limit of infinite supercell size. The results indicate that combining a technique that is known to yield accurate absolute core electron binding energies in free molecules (ΔSCF with SCAN) with an approach that yields accurate band energies of valence states (G 0 W 0 Γ) is a viable strategy for calculating core ) shifts it further onto a scale where the zero is defined by the position of the VBM predicted by G 0 W 0 @PBE (G 0 W 0 Γ@PBE). All energies are given in eV. electron binding energies in solids, referenced to the energy of the highest occupied state. Nevertheless, the smallness of the data set (only 15 binding energies) means that additional and more extensive tests are required to properly evaluate the accuracy of the SCAN + E G W @PBE 0 0 approach. In more general terms, we have demonstrated the importance of accurately predicting the position of the VBM in calculations of core electron binding energies, whenever the VBM is used as a point of reference. This includes not only calculations of periodic solids, but also calculations of surface species adsorbed onto a substrate with a band gap. We have found that using the conventional G 0 W 0 approach to predict the VBM energy gives unsatisfactory results. In contrast, using VBM energies predicted by the G 0 W 0 Γ approach, in which vertex corrections are included, yields excellent agreement between the calculated and experimental core electron binding energies. Other strategies for going beyond the G 0 W 0 @PBE level of theory, such as using a different mean-field starting point, or including partial self-consistency in GW, may give similar improvements, 36,44 and will be investigated in future studies. As an alternative with lower computational cost, hybrid functionals could be used to predict the VBM energy. This could be useful in cases where performing a GW calculation of the full unit cell of the material is prohibitively expensive.

COMPUTATIONAL METHODS
All of the ΔSCF calculations were performed using the allelectron electronic structure code FHI-aims. 42,45,46 The results of the calculations reported in ref 31, based on the SCAN functional, have been reused in this work to calculate core electron binding energies based on eq 1 and eq 3. In addition, similar calculations, using the same structures, settings, and numerical parameters, have been run using the PBE functional. Full details are provided in the Supporting Information of ref 31. GW and GWΓ calculations were run using GPAW. 47−49 In these calculations, the valence electrons are modeled using a plane wave basis set, and the effect of core electrons is treated using the projector-augmented wave formalism, as described in refs 47 and 48. In the ground state DFT calculations in GPAW, a plane wave cutoff of 800 eV was employed. The structures and the k-point grids used are given in the Supporting Information. Occupation smearing based on the Fermi−Dirac distribution with a width of 0.001 eV was applied in all cases. In the GW and GWΓ calculations, a nonlinear frequency grid defined by the values ω 2 = 20 eV and Δω 0 = 0.02 eV was used, where Δω 0 is the frequency spacing at ω = 0 and ω 2 is the frequency at which the spacing has increased to 2Δω 0 . For GWΓ, vertex corrections were calculated using the rAPBE kernel. GW and GWΓ calculations were performed at three values of E cut : 300, 350, and 400 eV, where E cut is the plane wave cutoff, and converged values were obtained by using a 1/ E cut 3/2 extrapolation. ■ ASSOCIATED CONTENT
Binding energies from eq 1 and eq 3 and extrapolation plots, numerical verification of eq 2 for all materials considered in this work, structures and k-point grids used in the GW and GWΓ calculations, and extrap-