We present a relativistic correction scheme to improve the accuracy of 1s core-level binding energies calculated from Green’s function theory in the GW approximation, which does not add computational overhead. An element-specific corrective term is derived as the difference between the 1s eigenvalues obtained from the self-consistent solutions to the non- or scalar-relativistic Kohn–Sham equations and the four-component Dirac–Kohn–Sham equations for a free neutral atom. We examine the dependence of this corrective term on the molecular environment and the amount of exact exchange in hybrid exchange–correlation functionals. This corrective term is then added as a perturbation to the quasiparticle energies from partially self-consistent and single-shot GW calculations. We show that this element-specific relativistic correction, when applied to a previously reported benchmark set of 65 core-state excitations [D. Golze et al., J. Phys. Chem. Lett. 11, 1840–1847 (2020)], reduces the mean absolute error (MAE) with respect to the experiment from 0.55 eV to 0.30 eV and eliminates the species dependence of the MAE, which otherwise increases with the atomic number. The relativistic corrections also reduce the species dependence for the optimal amount of exact exchange in the hybrid functional used as a starting point for the single-shot G0W0 calculations. Our correction scheme can be transferred to other methods, which we demonstrate for the delta self-consistent field (ΔSCF) approach based on density functional theory.

Core-level binding energies (BEs), measured by x-ray photoemission spectroscopy (XPS), are element-specific but also depend on the local chemical environment and thus afford access to information about the chemical bonding, oxidation state, and coordination of a given element in a sample.1–3 The energetic differences (chemical shifts) between atomic species of the same type can be smaller than 0.5 eV for second row elements and can be as low as 0.1 eV for carbon 1s excitations.4 The interpretation of an XPS spectrum can be very difficult due to overlapping features or the lack of well-defined reference data.5,6 Highly accurate theoretical tools for the prediction of relative and absolute BEs are, therefore, necessary to guide the experiment and its interpretation. The reliable computation of absolute core-level energies is generally more challenging than the calculation of energy shifts7 and is the focus of this work.

The most common approach to calculating core-level BEs is the delta self-consistent field (ΔSCF) method,8 wherein one computes the total energy difference between the ground and core-ionized states using Kohn–Sham density functional theory (KS-DFT). The best absolute core-level BEs have been obtained with meta-generalized gradient approximation (meta-GGA) functionals, yielding mean deviations of ≈0.2 eV with respect to the experiment for small molecules.9,10 A similar accuracy has been obtained with high-level wavefunction-based delta coupled cluster methods,11–13 albeit at much higher computational cost. The introduction of occupation constraints and the explicit generation of a charged system in Δ-based approaches lead to a plethora of problems.14,15 Most importantly, the application to periodic systems requires further approximations, e.g., neutralizing the unit or supercell.16,17

These problems do not occur in response theories, which avoid the explicit introduction of a core-ionized system and recently emerged as viable alternatives to Δ-based schemes for core-level calculations.18–20 One particularly promising response approach is the GW approximation21 to many-body perturbation theory. GW offers access to the core-state quasiparticle energies, which can be directly related to the core-level BEs. The GW approach is the standard approach to compute valence excitations (band structures) in solid-state physics.22 With the advent of efficient implementations with localized basis sets,23–28GW has also become the method of choice for calculating the BEs of frontier orbitals in molecules and is nowadays routinely applied to large molecular structures with several hundred atoms.27,29–31GW has also been successfully applied to deep valence and semi-core states with excitation energies up to 80 eV.32–34 However, the application of GW for deep core states with BEs larger than 100 eV, which are the ones relevant for the chemical analysis, has rarely been attempted. The first GW studies19,35–37 for deep core states reported partly large deviations of several electronvolts from the experiment. Recently, we have shown that GW can also be successfully applied for relative and absolute 1s molecular core-level BEs when utilizing highly accurate techniques for the integration of the self-energy38 and eigenvalue self-consistency in the Green’s function.20 

The application of GW (or any other method) to absolute core-level binding energies requires an accurate treatment of relativistic effects. Already for the 1s core levels of the p-block second period elements, the magnitude of relativistic effects begins to exceed the size of the chemical shifts. Recently, we have applied the GW method to a benchmark set of 65 core-level BEs of the second period elements C, N, O, and F in small and medium-sized molecules (henceforth, the CORE65 benchmark set), obtaining a mean deviation from the experiment of 0.3 eV.20 Therein, we applied a simple relativistic correction to the GW computed core-level energies. The purpose of the present article is to describe this correction.

Relativistic effects enter the GW formalism via the underlying reference calculation. A fully relativistic one-particle reference is described by a four-component Dirac spinor, or approximately by two-component spinors. Both types of spinors enable the description of noncollinear electronic states. GW equations with an explicit spin-dependent electron–electron interaction were only developed 12 years ago39,40 and provide a framework to properly describe spin–orbit coupling (SOC) effects in GW. Several GW codes emerged over the past years that implement these equations, employing spinors from all-electron two-component Dirac–KS-DFT calculations,41,42 KS-DFT with second-variation SOC43 or KS-DFT with fully relativistic pseudopotentials.44,45 These implementations are in the following referred to as fully relativistic GW approaches.

Fully relativistic GW calculations have been primarily used to compute valence and conduction bands of solids with strong SOC effects. Fully relativistic results have been reported for actinide metals,46,47 transition metal chalcogenides43,44 and dichalcogenides,45,48,49 perovskites,44,50 and bismuth-based topological insulators.51–56 Compared to non-relativistic GW, the fully relativistic approach is computationally at least four times more expensive.42 Alternatively, the SOC can be added as a perturbative correction to the quasiparticle band structure, which was common in earlier GW studies.57 Since this approach is computationally less expensive, it continues to be employed.58–60 However, it has been shown that the SOC post-correction scheme can fail, for example, to describe the band inversion in topological insulators.53 

An explicit treatment of the SOC is typically not necessary for valence excitations of molecules, unless very heavy atoms are involved. Since molecular GW studies have mainly focused on organic semiconductors,22 fully relativistic GW calculations are rare and have been mostly conducted for diatomic molecules,41,44 but recently also for transition metal halide complexes.42 If relativistic effects are considered in molecular GW calculations at all, they are more commonly included by employing one-component reference states that capture only scalar-relativistic effects. This approach has been employed in GW calculations with scalar-relativistic pseudopotentials23,27 or the zeroth-order regular approximation (ZORA).31 

Relativistic considerations for the innermost core-levels differ from those for valence excitations. SOC leads, for core states with the angular quantum number l > 0, to a splitting. Already for 2p states of third-row elements, e.g., sulfur, this splitting lies in the range of 1 eV and increases for fourth-period elements to several tens of eV.61 SOC affected core states require, thus, a noncollinear treatment, and we are only aware of one GW study41 that reports results for deep, spin–orbit-coupled p states. In this study, we focus on organic molecules and small inorganic molecules with elements C, N, O, and F, for which, typically, the 1s excitations in the energy range from 250 eV to 700 eV are measured in XPS. While scalar relativistic effects are heightened in the proximity of the poorly screened nuclear charge, SOC is suppressed for 1s core states. In contrast to valence excitations, the most common scalar-relativistic approximation, ZORA (or, in the variant used in this work and explained in Sec. II D 2, atomic ZORA62,63), performs poorly for the absolute eigenvalues associated with the innermost core levels.64,65 In order to avoid the computational expense of a fully Dirac-like, two-component GW calculation, which is not essential for the physics of 1s core excitations, we derive here an element-specific relativistic corrective term for non-relativistic and scalar-relativistic reference states, which we add in a post-GW perturbative step.

The remainder of this article describes this correction. In Sec. II, we describe the GW formalism, highlighting aspects that are particularly relevant for core-level calculations, and follow this up with an overview of the aspects of relativistic theory relevant to this work. We then describe the methods employed in Sec. III. We present and discuss the results of our correction schemes, which we apply to GW and ΔSCF computed core-level BEs of the CORE65 benchmark set in Sec. IV and finally draw conclusion in Sec. V.

In practice, GW is often performed as a one-shot perturbative approach (G0W0) on top of an underlying mean field theory calculation. Possible mean field theories are Hartree–Fock (HF), KS-DFT, or hybrid DFT, which yield the molecular orbitals (MOs) {ϕn} and eigenvalues {ϵn} used as an input (starting point) for the G0W0 calculation. The G0W0 quasiparticle (QP) energies ϵ n G 0 W 0 are obtained by solving the QP equation,

ϵ n σ G 0 W 0 = ϵ n σ + R e Σ n σ ( ϵ n σ G 0 W 0 ) v n σ x c ,
(1)

and can be related to the BE of state n by B E n σ = ϵ n σ G 0 W 0 , where σ denotes the spin index. The nth diagonal elements of the KS exchange–correlation (XC) potential and self-energy operator Σ are denoted by v n σ x c = ϕ n σ | v x c | ϕ n σ and Σ n σ = ϕ n σ | Σ σ | ϕ n σ , respectively. The self-energy operator is given by

Σ σ ( r , r , ω ) = i 2 π d ω G 0 σ ( r , r , ω + ω ) W 0 ( r , r , ω ) e i ω η ,
(2)

where G 0 σ is the non-interacting KS Green’s function, W0 is the screened Coulomb interaction, and η is a positive infinitesimal. The KS Green’s function is obtained from the KS orbitals and eigenvalues by

G 0 σ ( r , r , ω ) = m ϕ m σ ( r ) ϕ m σ ( r ) ω ϵ m σ i η sgn ( ϵ F ϵ m σ ) ,
(3)

where ϵF is the Fermi energy and the sum runs over both the occupied and virtual KS orbitals. Within the random phase approximation, the screened Coulomb interaction is given by

W 0 ( r , r , ω ) = d r ϵ 1 ( r , r , ω ) v ( r , r ) ,
(4)

with the bare Coulomb interaction v(r, r′) = 1/|rr′| and the dielectric function, ϵ,

ϵ ( r , r , ω ) = δ ( r , r ) d r v ( r , r ) χ 0 ( r , r , ω ) ,
(5)

where the irreducible polarizability χ0 is given in the real-space Adler–Wiser representation66,67 as

χ 0 ( r , r , ω ) = σ i o c c a v i r t ϕ a σ ( r ) ϕ i σ ( r ) ϕ a σ ( r ) ϕ i σ ( r ) × 1 ω ( ϵ a σ ϵ i σ ) + i η + 1 ω ( ϵ a σ ϵ i σ ) + i η .
(6)

The index i runs over occupied orbitals, and a runs over virtual orbitals. Note that we drop the spin index σ in the following for simplicity.

We split the self-energy in a correlation part Σc and an exchange part Σx, where Σ = Σc + Σx. The exchange term is frequency-independent and calculated as the HF-like exchange with the KS orbitals, as described in Ref. 38. Details on the calculation of the correlation part are given in Sec. II B.

Equation (1) is non-linear and can be solved iteratively or approximately by linearization using a Taylor expansion to the first order around ϵn.22 As we pointed out in our previous work,20,38 the linearization error increases rapidly with an increase in BE and can already amount to 0.5 eV for deeper valence states.22 The magnitude of this error is in the range of the chemical shifts expected for 1s excitations. In addition, core-level BEs are an order of magnitude larger than deep valence BEs, potentially leading to even larger linearization errors. We, therefore, always solve the QP equation iteratively.

The accurate frequency integration of the self-energy [Eq. (2)] is one of the major challenges for the calculation of deep core states. A common approach for valence states is to evaluate the self-energy for imaginary frequencies and analytically continue it to the real frequency axis by fitting the self-energy matrix elements to a multipole model. Analytic continuation is employed in many state-of-the-art GW implementations24,27,29,68 and yields accurate results for valence states.69 The structure of Σn for a valence state is typically smooth in the frequency region where the QP solution is expected and is well reproduced by analytic continuation. For core states, the self-energy has a complicated structure with many poles. We showed that analytic continuation becomes numerically unstable in the core region and completely fails to reproduce the self-energy structure.38 

For core states, a more accurate evaluation of the self-energy on the real frequency axis is required. We employ the contour deformation (CD) technique to evaluate the correlation self-energy Σc,22,23,38,70–72 where the numerically unstable integration along the real frequency axis is avoided by extending the integrand to the complex plane. The contours are chosen such that only the poles of G0 are enclosed in the contours, and the contour integral is evaluated using the residue theorem; see Refs. 22 and 38 for details. The integral along the real frequency axis is then evaluated as

Σ c ( r , r , ω ) = 1 2 π d ω G 0 ( r , r , ω + i ω ) W 0 c ( r , r , i ω ) i ϕ i ( r ) ϕ i ( r ) W 0 c ( r , r , | ϵ i ω | + i η ) θ ( ϵ i ω ) + a ϕ a ( r ) ϕ a ( r ) W 0 c ( r , r , | ϵ a ω | + i η ) θ ( ω ϵ a ) ,
(7)

where θ is the Heaviside step function, W 0 c ( r , r , ω ) = W 0 ( r , r , ω ) v ( r , r ) , and i refers to occupied and a to unoccupied orbitals.

The CD technique reproduces the self-energy structure for core excitations exactly and matches the results from the computationally more expensive fully analytic solution of Eq. (2).38 Recently, combining the CD approach with analytic continuation of W0 has been proposed73 and might be an alternative approach for the accurate calculation of the self-energy for core states.74 

In this section, we briefly present the GW variants that we used in this work. By now, many different GW flavors have emerged in practical calculations.22 The most common flavor, G0W0 based on a semi-local DFT starting point, breaks down for core states, as we will detail in the following. We, therefore, need to go beyond the most common approach.

The problem with the conventional GW approach is related to a loss of spectral weight in the QP peak in the spectral function A ( ω ) = 1 / π m I m G m , where m runs over all occupied and virtual states and G = G0 + G0ΣG. For molecular valence states, G0W0 performed on top of a DFT calculation with the Perdew–Burke–Ernzerhof (PBE)75 functional (G0W0@PBE) yields a clearly identifiable QP peak. This peak corresponds to a distinct solution of Eq. (1). Multiple solutions, which would indicate spectral weight transfer to other peaks in the spectral function, have only been observed for a few systems for frontier orbitals.69,76–78 These are rare cases and usually not more than two possible solutions are observed.

The situation is dramatically different for deep core states. The analysis of the spectral functions in our recent work showed that a unique QP solution is not obtained with G0W0@PBE for 1s states.20 The spectral function shows a multitude of peaks with similar spectral weight, but no distinct QP excitation. A spurious transfer of spectral weight from the QP peak to plasmon satellites has been previously observed for deep valence states of transition metal oxides79,80 and semi-core excitations of sodium.33 However, for deep core states, the transfer of spectral weight is far more extreme resulting in a spectral function, where the satellite spectrum and quasiparticle have completely merged. Such a spectral function contradicts the expected physics. Photoemission spectra of molecular 1s excitations show strong QP peaks accompanied by satellites features due to multi-electron excitations such as shake-up processes, which are orders of magnitude smaller than the main excitation.81,82

We showed that the 1s QP peak can be correctly restored by including eigenvalue self-consistency in G, while keeping PBE as the starting point.20 This scheme is referred to as evGW0@PBE. In evGW0, the screened Coulomb interaction is kept fixed at the W0 level, and the Green’s function is recomputed replacing the mean field eigenvalues with the QP energies from Eq. (1). We enforce eigenvalue self-consistency only in G. Inserting the QP energies also in W0 (evGW) reduces the screening, which is not advantageous because the overscreening in W at the PBE level compensates the underscreening due to missing vertex corrections. It has been shown that evGW0 yields bandgaps in good agreement with the experiment,83 while underscreening errors in the evGW scheme lead to too large bandgaps83 and overly stretched spectra.84 

Higher-level self-consistency schemes, such as fully self-consistent GW (scGW)85–87 or quasiparticle self-consistent GW (QSGW),88 are expected to restore the 1s QP peak as well but might not yield better agreement with the experiment than evGW0. It has been shown that scGW overestimates molecular HOMO excitations89 and bandgaps in solids.90 Similar underscreening effects are also expected for core states. A first exploratory study seems to confirm this assumption for QSGW, reporting an overestimation of 2 eV for 1s core states of small molecules.37 

evGW0 is computationally more demanding than a G0W0 calculation because the QP equation is not only solved for the 1s core states of interest but also repeatedly for all occupied and virtual states until convergence in G is reached. We showed that the core-level QP peak can also be restored in a G0W0 calculation by using a XC functional with a high fraction of exact exchange as the starting point.20 We employ the PBEh(α) functional family with an adjustable amount α of HF exact exchange.91 The XC energy Exc is given by

E x c = α E x E X + ( 1 α ) E x P B E + E c P B E , α [ 0,1 ] ,
(8)

where E x E X denotes the HF exchange energy. E x P B E and E c P B E are the PBE exchange and correlation energy, respectively.

In this work, we followed Ref. 20 and used both evGW0 and G0W0@PBEh(α). We analyze how the relativistic corrections we devised affect the two schemes.

It has long been recognized that relativistic effects play a large role in the chemistry of heavy elements.92,93 In this work, we treat the core states of light elements carbon through fluorine, whose BEs are in the range of 250 eV–700 eV, and for which relativistic effects are usually smaller than 1 eV. However, the accuracy required to resolve XPS spectra of 1s excitations of second row elements is in the range of some tenths of an electronvolt and, therefore, on the same order of magnitude as the relativistic effects.

Our relativistic correction scheme for GW is based on two different relativistic KS-DFT methods. The first is the four-component Dirac Kohn–Sham (4c-DKS) approach, further also referred to as the fully relativistic scheme. The second uses the scalar-relativistic atomic ZORA.

1. Fully relativistic Dirac approach

The relativistic description of a non-interacting electron in an external potential V is given by the Dirac equation,94 

h D Ψ = ϵ Ψ ,
(9)

where Ψ = ϕ χ is a four-component Dirac spinor, which is comprised of a large component ϕ and a small component χ, each of which has two components for the spin functions. The Dirac Hamiltonian shifted by the rest energy c2 is given by

h D = V c σ p c σ p 2 c 2 + V ,
(10)

where c is the speed of light, σ is a vector of Pauli spin matrices, and p is the momentum operator. For electron-like states in the non-relativistic limit (c → ∞), the large component ϕ reduces to the wave function of the Schrödinger equation, while the small component vanishes.

For the case with N interacting electrons, the electronic relativistic Hamiltonian is

H D = i N h D ( i ) + i < j N g ( i , j ) ,
(11)

where g(i, j) is the electron–electron interaction. In the non-relativistic case, g(i, j) corresponds to the Coulomb operator, where the interaction between two electrons is instantaneous. When including relativity, this cannot be correct because the Coulomb interaction between electrons involves the exchange of photons traveling at the speed of light. The relativistic electron–electron operator is much more complicated than the non-relativistic one and cannot be written in the closed form. Its perturbation expansion in terms of 1/c2 yields

g ( i , j ) = n = 0 1 c 2 n g n ( i , j ) = 1 | r i r j | + O ( c 2 ) ,
(12)

where g0 is the instantaneous Coulomb interaction and the first order correction g1 is the Breit term,95,96 which introduces magnetic and retardation effects. In our 4c-DKS calculations, we include only the g0 term. The contributions from the Breit term are believed to be small,97–100 and they are, therefore, neglected in most relativistic calculations.

In relativistic KS-DFT,101–105 the XC functional should, in principle, also include relativistic effects and should be formulated in terms of the four-current density.106 The latter is the basic density variable in relativistic KS-DFT. A relativistic generalization of the local density approximation (RLDA)103 has been proposed, as well as a semi-empirical gradient corrected variant (RGGA).102,107 Common practice, however, is to use a non-relativistic XC functional in conjunction with the Dirac kinetic energy,100,106 which is the procedure we follow in this work. Here, we use a 4c-DKS approach with non-relativistic GGA and hybrid GGA functionals.

2. Scalar-relativistic atomic ZORA approach

The computational cost for a fully relativistic 4c-DKS approach is significantly higher than for non-relativistic Schrödinger Kohn–Sham (SKS). The scalar-relativistic atomic ZORA [see Eq. (18)] retains the computational effort of an SKS calculation and has been shown to capture relativistic effects in good agreement with other scalar-relativistic all-electron schemes.62,63,108

The original ZORA scheme is derived by solving one of the two coupled equations in Eq. (9) for the small component χ and inserting it into the other equation, which yields the following (still exact) expression for the un-normalized large component:

σ p c 2 2 c 2 V 1 + ϵ n 2 c 2 V 1 σ p ϕ n + V ϕ n = ϵ n ϕ n ,
(13)

and can be rewritten as

σ p c 2 2 c 2 + ϵ n V σ p ϕ n + V ϕ n = ϵ n ϕ n .
(14)

Expanding the parenthetical term in Eq. (13) as a geometric series yields the regular approximation,105 

σ p c 2 2 c 2 V k = 0 ϵ n V 2 c 2 k σ p ϕ n + V ϕ n = ϵ n ϕ n .
(15)

Retaining only the scalar part of the zeroth order term (k = 0) transforms Eq. (15) to

( T Z O R A + V ) ϕ n = ϵ n ϕ n ,
(16)

where the ZORA kinetic energy TZORA is defined as109 

T Z O R A = p c 2 2 c 2 V p .
(17)

The zeroth order is the only order in this expansion that leaves the Hamiltonian independent of ϵn, enabling a standard numerical treatment of the ensuing KS equations. From Eq. (14), it is apparent that the ZORA is equivalent to the supposition that ϵn ≪ 2c2. This relation is correct for valence states, i.e., ZORA correctly reflects their scalar-relativistic behavior. In contrast, ϵn ≪ 2c2 is a potentially severe approximation for the deep core states. However, the core states are relatively inert, i.e., the chemical shifts are small compared to the absolute value of ϵn. The error introduced for core states by neglecting the higher order terms ( ϵ n / ( V 2 c 2 ) ) k in Eq. (15) is, therefore, largely an atomic property. For a given atomic species, orbital, and orbital occupation, the dependence of this error on different chemical environments is minor. Physical quantities that depend on a difference in energies are, therefore, not strongly affected.

Since the potential V enters in the denominator of Eq. (17), rather than as ϵnV, it is clear that ZORA is gauge dependent, i.e., a constant shift in the electrostatic potential does not lead to a constant shift in the energy. In practice, this dependence on the choice of the zero-energy level of the potential in different chemical environments would be a severe problem. Different methods have been proposed to restore gauge-invariance. One of them is the scaled ZORA,110 where the eigenvalues are rescaled after self-consistency is reached. Scaled ZORA restores almost, but not completely gauge-invariance.

Full gauge-invariance is achieved in the atomic ZORA scheme (aZORA), which we use in this work. We employ the aZORA approach as defined in Eqs. (55) and (56) of Ref. 62 and benchmarked in Refs. 63 and 108. In aZORA, the potential in the denominator of TZORA is replaced with the potential vat(j) of a radially symmetric free atom in empty space, which is calculated once at the outset of the calculation and, therefore, does not depend on the chemical environment in any way. The potential vat(j) is evaluated using boundary conditions with zero potential at infinity (i.e., non-periodic boundary conditions). The dependence on index j indicates that vat(j) corresponds to the potential of the atomic species at the center where the basis function φj is localized. For semilocal DFT calculations with, e.g., the PBE functional, we construct vat(j) using the same functional. For our hybrid DFT calculations with the PBEh family, vat(j) is, instead, constructed using PBE exchange–correlation for simplicity. This approximation is entirely minor for the same reasons as using vat(j) instead of V in the first place. The dominant effect of either V or vat(j) in Eq. (17) with respect to the 2c2 term is the deep electron–nucleus potential, modified by the overall electrostatic screening of the surrounding core electrons. In contrast, differences associated with vxc are negligible compared to the magnitude of 2c2.

The atomic potential vat(j) introduces a dependence on the localized basis functions φj in the kinetic term TaZORA. Therefore, the matrix elements need to be symmetrized to restore Hermiticity, which finally gives

φ i | T a Z O R A | φ j = 1 2 φ i | p c 2 2 c 2 v a t ( j ) p | φ j + 1 2 φ j | p c 2 2 c 2 v a t ( i ) p | φ i .
(18)

While the absolute values of the scaled ZORA eigenvalues are closer to the 4c-DKS reference,110 we expect the relative shifts with respect to 4c-DKS, which are relevant for the proposed correction scheme, to be more consistent with aZORA. The reason is that the latter, unlike scaled ZORA, restores the gauge-invariance completely.

For GW, we have developed three simple correction schemes to account for relativistic effects: (I) Atomic relativistic corrections are added to the QP energies. (II) The aZORA Hamiltonian is used for the underlying DFT calculation, and the obtained KS eigenvalues and MOs are used as a starting point for GW. (III) aZORA is used as in II, and atomic relativistic corrections are added to the QP energies. The atomic corrections are always added as a post-processing step to the converged QP energies and have been obtained as follows.

For scheme I, the atomic relativistic corrections Δ ϵ 1 s , a t S K S are computed as the difference between the non-relativistic SKS 1s eigenvalues ( ϵ 1 s , a t S K S ) and the fully relativistic 4c-DKS 1s eigenvalues ( ϵ 1 s , a t 4 c - D K S ) ,

Δ ϵ 1 s , a t S K S = ϵ 1 s , a t 4 c - D K S ϵ 1 s , a t S K S .
(19)

The label “at” indicates that the calculations are performed for a free neutral atom. For scheme III, we use the atomic corrections Δ ϵ 1 s , a t a Z O R A ,

Δ ϵ 1 s , a t a Z O R A = ϵ 1 s , a t 4 c - D K S ϵ 1 s , a t a Z O R A ,
(20)

evaluating the difference to the aZORA 1s eigenvalues ( ϵ 1 s , a t a Z O R A ) instead to the SKS eigenvalues. The atomic 1s eigenvalues ϵ 1 s , a t S K S , ϵ 1 s , a t a Z O R A , and ϵ 1 s , a t 4 c - D K S are computed self-consistently at the PBE level by solving the radial SKS, aZORA, and 4c-DKS equations, respectively.

All GW and ΔSCF calculations are performed with the all-electron FHI-aims program package,24,62,111 which is based on numerically tabulated atom-centered orbitals (NAOs). Core-level BEs from G0W0, evGW0, and ΔSCF calculations are calculated for the CORE65 benchmark set introduced in Ref. 20, which contains 65 1s binding energies of second row elements (C, N, O, and F) for 31 closed- and one open-shell molecule. The settings for G0W0, evGW0, and ΔSCF are the same as in our previous work20 and are summarized in the following.

The ΔSCF calculations are performed with the PBE0112,113 hybrid functional employing def2 quadruple-ζ valence plus polarization (def2-QZVP)114 basis sets. The all-electron def2-QZVP Gaussian basis sets are treated numerically in FHI-aims for compliance with the NAO scheme. We decontract the def2-QZVP basis sets to enable a full relaxation of the other electrons in the presence of a core-hole; see Ref. 20 for further details and an explanation of the basis set choice.

For the GW calculations, the QP equation [Eq. (1)] is always solved iteratively. For the partially self-consistent evGW0 scheme, we iterate the eigenvalues additionally in G. We use the PBE functional75 as a starting point for the evGW0 calculations. As discussed in Sec. II C, a distinct 1s QP solution does not exist at the G0W0@PBE level, i.e., for the initial evGW0@PBE step. We, therefore, initialize the evGW0 loop with approximate 1s QP solutions obtained, e.g., by linearizing the 1s QP equation or employing a larger value for the broadening parameter η. Note that this is only done for the first evGW0@PBE steps. For G0W0, we employ the PBEh(α) hybrid functionals91 for the underlying DFT calculation, where α indicates the fraction of HF exchange in the functional. The core-level BEs are extrapolated to the complete basis set limit to account for the slow convergence of the GW QP energies with respect to the basis set size.22,27,69,115,116 All empty states that can be resolved by the given basis set are included in the calculation of the self-energy. The extrapolation is performed by a linear regression with respect to the inverse of the total number of basis functions using the Dunning basis set family cc-pVnZ (n = 3–6).117,118 Details are given in Table S3 of the supplementary material, and comprehensive convergence studies are presented in Fig. S1. Furthermore, we use the CD technique38 to compute the GW self-energy. The integral over the imaginary frequency axis in Eq. (7) is computed using modified Gauss–Legendre grids24 with 200 grid points.

Relativistic effects for GW are included in three different ways, as described in Section IIE. For ΔSCF, we account for relativistic effects self-consistently using the aZORA.62 We also apply the atomic relativistic schemes introduced in Sec. II E to ΔSCF for comparison. To obtain the atomic relativistic corrections for Eqs. (19) and (20), the radial DKS, SKS, and aZORA-KS equations are solved self-consistently on numerical real-space grids with the DFTATOM code119 incorporated into FHI-aims.

We investigate the dependence of the relativistic eigenvalue corrections on the molecular environment, XC functional, and basis set using the DIRAC program,120,121 which features a 4c-DKS DFT implementation for the 3D electronic wave function, enabling also molecular calculations. Similar to Eq. (19), we define the molecular corrections as

Δ ϵ 1 s , m o l S K S = ϵ 1 s , m o l 4 c - D K S ϵ 1 s , m o l S K S ,
(21)

where ϵ 1 s , m o l 4 c - D K S are molecular 1s eigenvalues of the 4c-DKS Hamiltonian. The corresponding non-relativistic eigenvalues ϵ 1 s , m o l S K S are here obtained from a 4c-DKS calculation, resetting the speed of light to the non-relativistic limit (c → ∞).

The DIRAC calculations are performed for the molecular structures of the CORE65 benchmark set, excluding the spin-polarized O2 case, using all-electron Dyall basis sets122 of triple-zeta quality and the PBE functional. We define the difference ΔMOL between the molecular and atomic eigenvalue corrections as

Δ M O L = Δ ϵ 1 s , a t S K S Δ ϵ 1 s , m o l S K S .
(22)

The functional dependence of the atomic corrections is assessed for the PBEh(α) hybrid family. We also study the basis set dependence for the Dyall series122 with reference to the fully converged radial solution from DFTATOM,

Δ B A S = Δ ϵ 1 s , a t S K S ( d y a l l ) Δ ϵ 1 s , a t S K S ( r a d i a l ) .
(23)

We first present the atomic relativistic corrections and discuss their dependence on technical and convergence parameter, the XC functional, and the molecular environment. We proceed with a discussion of non-relativistic results for the CORE65 benchmark set and demonstrate how our simple correction schemes, based on these atomic corrections, improve the agreement of the computed absolute 1s BEs to the experiment.

Figure 1(a) shows a sketch of the 1s eigenvalues from the non-relativistic SKS, 4c-DKS, and scalar-relativistic aZORA calculations. The SKS eigenvalues are generally overestimated with respect to the 4c-DKS reference, while aZORA underestimates the 1s eigenvalues by nearly as much. The atomic eigenvalue corrections Δϵ1s,at for SKS [Eq. (19)] and aZORA [Eq. (20)] are given in Table I. The SKS corrections are negative and increase in magnitude with the atomic number, ranging from −4 meV for Li to −13.4 eV for Ar. The aZORA corrections are positive and increase from 4 meV (Li) to 12.4 eV (Ar).

FIG. 1.

(a) Schematic of an energy level diagram of 1s eigenvalues comparing non-relativistic Schrödinger Kohn–Sham (SKS), fully relativistic four-component Dirac–Kohn–Sham (4c-DKS), and scalar-relativistic aZORA. (b) Atomic relativistic Δϵ1s,at for the spherical SKS [Eq. (19)] and aZORA Hamiltonian [Eq. (20)] with respect to the atomic number.

FIG. 1.

(a) Schematic of an energy level diagram of 1s eigenvalues comparing non-relativistic Schrödinger Kohn–Sham (SKS), fully relativistic four-component Dirac–Kohn–Sham (4c-DKS), and scalar-relativistic aZORA. (b) Atomic relativistic Δϵ1s,at for the spherical SKS [Eq. (19)] and aZORA Hamiltonian [Eq. (20)] with respect to the atomic number.

Close modal
TABLE I.

Atomic eigenvalue corrections Δϵ1s,at for the SKS [Eq. (19)] and aZORA Hamiltonian [Eq. (20)] for the second and third period elements. Z indicates the atomic number.

Z Excitation SKS (eV) aZORA (eV)
Li1s  −0.004 245  0.003 683 
Be1s  −0.016 54  0.015 16 
B1s  −0.050 18  0.043 54 
C1s  −0.117 6  0.100 1 
N1s  −0.235 5  0.199 6 
O1s  −0.424 4  0.359 3 
F1s  −0.708 0  0.600 0 
10  Ne1s  −1.113  0.945 6 
11  Na1s  −1.658  1.442 
12  Mg1s  −2.387  2.116 
13  Al1s  −3.360  3.009 
14  Si1s  −4.607  4.161 
15  P1s  −6.180  5.620 
16  S1s  −8.129  7.435 
17  Cl1s  −10.51  9.664 
18  Ar1s  −13.39  12.37 
Z Excitation SKS (eV) aZORA (eV)
Li1s  −0.004 245  0.003 683 
Be1s  −0.016 54  0.015 16 
B1s  −0.050 18  0.043 54 
C1s  −0.117 6  0.100 1 
N1s  −0.235 5  0.199 6 
O1s  −0.424 4  0.359 3 
F1s  −0.708 0  0.600 0 
10  Ne1s  −1.113  0.945 6 
11  Na1s  −1.658  1.442 
12  Mg1s  −2.387  2.116 
13  Al1s  −3.360  3.009 
14  Si1s  −4.607  4.161 
15  P1s  −6.180  5.620 
16  S1s  −8.129  7.435 
17  Cl1s  −10.51  9.664 
18  Ar1s  −13.39  12.37 

The atomic corrections given in Table I are visualized in Fig. 1(b). We observe that the magnitude of the atomic corrections for both SKS and aZORA depends on the fourth power of atomic number Z. This dependence for the Schrödinger equation is known from the relativistic correction to the exact energy of the hydrogenic orbital, whose leading order term in a perturbative expansion scales as the fourth power of Z,94,123–125 and a similar result has been discussed for the ZORA Hamiltonian as well.64 As the 1s orbitals are poorly screened by the outer orbitals, the magnitude of the relativistic correction trends similarly to that of the unscreened hydrogenic orbitals.

The results reported in Table I are obtained from the self-consistent solutions of the radial SKS and 4c-DKS equations. The radial equations enforce a spherical symmetry of the solution. However, most atoms have ground states with non-spherical symmetry. For the second row, this applies to B, C, O, and F. For these elements, the spherical solutions are too high in total energy by several tenths of eV and assume fractional occupation numbers. The 1s eigenvalue corrections Δ ϵ 1 s , a t S K S obtained from the radial SKS and 4c-DKS equations are, therefore, an approximation. To estimate the error introduced by this approximation, we solved the 3D SKS equations for the free neutral atom and compared the 1s eigenvalues of the spherical and non-spherical solution. The spherical solution is also obtained with the 3D equations and is identical to the radial one if we do not break the symmetry and enforce integer occupations. The non-spherical 3D solution is obtained by employing occupation constraints. We find that the difference in the absolute 1s eigenvalues between the spherical and non-spherical solutions is less than 50 meV, see Table S1 of the supplementary material, which is an order of magnitude smaller than the relativistic corrections themselves. The error in the relative values, Δ ϵ 1 s , a t S K S , is expected to be even smaller, and thus, we conclude that the radial approximation is sufficient.

The radial calculations are performed on a numeric real-space grid, which can be easily converged, whereas the 3D calculations rely on relativistic all-electron Gaussian basis sets, potentially introducing a basis set incompleteness error. Figure 2(a) shows the basis set convergence of the Dyall series with respect to the radial solution. At the double-ζ level, the error is within a few meV, and for the relevant 1s states, we reach convergence already at the triple-ζ level; see Fig. 2(a). We use the quadruple-ζ basis set for the calculations shown in Fig. 2(b) and the triple-ζ basis set for the calculations in Fig. 2(c).

FIG. 2.

Basis set convergence and dependence of the atomic corrections on the XC functional and molecular environment. (a) Difference ΔBAS as defined in Eq. (23) between the Dyall all-electron basis set at the double, triple, and quadruple-ζ levels122 and the radial solution on numeric real-space grids. (b) Atomic eigenvalue correction Δ ϵ 1 s , a t S K S [Eq. (19)] computed with the PBEh(α) functional with three different values of α. (c) Difference ΔMOL as defined in Eq. (22) between the molecular eigenvalue correction and the atomic value for C1s, N1s, O1s, and F1s excitations of the CORE65 benchmark set. Bars contain the second and third quartiles, whiskers extend to encompass 95% of the results, and outliers are shown as dots.

FIG. 2.

Basis set convergence and dependence of the atomic corrections on the XC functional and molecular environment. (a) Difference ΔBAS as defined in Eq. (23) between the Dyall all-electron basis set at the double, triple, and quadruple-ζ levels122 and the radial solution on numeric real-space grids. (b) Atomic eigenvalue correction Δ ϵ 1 s , a t S K S [Eq. (19)] computed with the PBEh(α) functional with three different values of α. (c) Difference ΔMOL as defined in Eq. (22) between the molecular eigenvalue correction and the atomic value for C1s, N1s, O1s, and F1s excitations of the CORE65 benchmark set. Bars contain the second and third quartiles, whiskers extend to encompass 95% of the results, and outliers are shown as dots.

Close modal

In Fig. 2(b), we examine the dependence of the atomic eigenvalue correction Δ ϵ 1 s , a t S K S on the fraction α of exact HF exchange in the PBEh(α) functional for α = 0 (PBE), α = 0.25 (PBE0), and α = 0.5. The magnitude of the eigenvalue correction shows a slight dependence on α, increasing by an amount that is proportional to the fraction of exact exchange. At first glance, the α dependence seems more pronounced for heavier elements. However, this is only true for the absolute values: setting the PBE functional as reference and comparing to PBEh(α = 0.5), the magnitude of the atomic correction increases by 12 meV for carbon 1s, which corresponds to 10.2%, and 40 meV for fluorine 1s, which, however, corresponds only to 5.7%. In fact, the α dependence seems to decrease with the atomic number when comparing relative deviations; see Table S6 of the supplementary material for the tabulated values. For all elements listed in Table I, we find that the α dependence is an order of magnitude smaller than the relativistic correction itself. We thus neglect it when applying our relativistic correction schemes to the 1s QP energies from G0W0@PBEh.

In Fig. 2(c), we compare the atomic eigenvalue correction [Eq. (19)] to the molecular eigenvalue correction [Eq. (21)]. For most of the excitations considered, the atomic eigenvalue correction slightly underestimates in magnitude the molecular correction, but the difference between the two is below 5 meV for 49 of the 63 excitations considered, with a maximum deviation of 12.6 meV. The distribution of these differences is similar for the core excitations of different elements and small enough in comparison with the magnitude of the eigenvalue correction to justify the use of the atomic values irrespective of the chemical environment.

The atomic SKS corrections reported in Table I are very similar to the atomic corrections published in Ref. 9 for second row elements B → F. The atomic corrections in Ref. 9 are 10 meV–40 meV larger than ours and were computed by comparing the four-component Dirac–Hartree–Fock (DHF) energies with the non-relativistic HF energies. Our analysis of different PBEh functionals in Fig. 2(b) suggests that these differences must be partly attributed to the exchange treatment. The remaining differences might be due to the usage of non-relativistic basis sets in combination with the 4c-DHF Hamiltonian in Ref. 9. The atomic SKS corrections are also surprisingly similar to the corrections derived for second period elements in an early work from the 1960s based on Pauli perturbation theory of charged two-electron atoms.126,127 Pauli perturbation theory is based on the first order in the expansion of Eq. (13) in terms of 1/c2. It is highly singular in the deep core region128 and has been largely replaced by the regular approximation, which expands Eq. (13) in terms of ϵ n V 2 c 2 . The correspondence worsens when valence electrons are included.129 

In our previous work,20 we briefly discussed the effect of relativistic corrections, comparing non-relativistic and relativistic 1s BEs from evGW0 to experiment. We will now analyze the non-relativistic evGW0 calculations in more detail and additionally include the non-relativistic G0W0@PBEh and ΔSCF results in the discussion.

Figure 3 displays the mean absolute error (MAE) of the absolute 1s BEs with respect to the experiment for the CORE65 benchmark set. The MAEs obtained from non-relativistic evGW0 calculations increase with the atomic number [Fig. 3(c)], and the magnitude of this increase is within the range of the atomic relativistic corrections given in Table I. The distribution of these errors is shown in Fig. 4(a), where the grouping in species is evident. evGW0 systematically underestimates the 1s BEs for all 65 excitations. The non-relativistic ΔSCF calculations underestimate the 1s BEs as well, and the MAEs show a very similar trend with respect to the atomic number; see Fig. 3(c). Comparing the overall MAE for the non-relativistic calculations, we find that the MAE for ΔSCF is 0.83 eV, which is slightly larger than the 0.55 eV MAE for evGW0.

FIG. 3.

Mean absolute error (MAE) of absolute 1s BEs with respect to the experiment for the CORE65 benchmark set. MAE for G0W0@PBEh dependent on the fraction of exact exchange α in the PBEh functional (a) without and (b) with the atomic relativistic correction (RC). (c) MAE for ΔSCF and evGW0@PBE, without and with RC.

FIG. 3.

Mean absolute error (MAE) of absolute 1s BEs with respect to the experiment for the CORE65 benchmark set. MAE for G0W0@PBEh dependent on the fraction of exact exchange α in the PBEh functional (a) without and (b) with the atomic relativistic correction (RC). (c) MAE for ΔSCF and evGW0@PBE, without and with RC.

Close modal
FIG. 4.

Distribution of errors with respect to the experiment for absolute 1s BEs of the CORE65 benchmark set, where E r r o r i = BE i theory BE i exp . Note that the histogram is stacked. Three different relativistic correction (RC) schemes are compared. (a) Non-relativistic evGW0, (b) evGW0 adding the atomic corrections ϵ 1 s , a t S K S to the QP energies, (c) evGW0 using aZORA in the underlying DFT calculation, (d) evGW0 using aZORA as in (c) and adding atomic corrections ϵ 1 s , a t a Z O R A to the QP energies.

FIG. 4.

Distribution of errors with respect to the experiment for absolute 1s BEs of the CORE65 benchmark set, where E r r o r i = BE i theory BE i exp . Note that the histogram is stacked. Three different relativistic correction (RC) schemes are compared. (a) Non-relativistic evGW0, (b) evGW0 adding the atomic corrections ϵ 1 s , a t S K S to the QP energies, (c) evGW0 using aZORA in the underlying DFT calculation, (d) evGW0 using aZORA as in (c) and adding atomic corrections ϵ 1 s , a t a Z O R A to the QP energies.

Close modal

Relativistic effects are also apparent when considering the optimal α for use in a G0W0@PBEh(α) scheme. Figure 3(a) shows the MAE for non-relativistic G0W0@PBEh(α) calculations with respect to the fraction of exact exchange α. These calculations have been carried out for a subset of 43 excitations of the CORE65 benchmark set, for which the mapping between the core state and the atom does not require the analysis of, e.g., MO coefficients. In our previous work,20 we reported the α dependence of the MAE including relativistic effects [Fig. 3(b)] and found that the smallest MAE is obtained for α values around 0.45. For MAEs smaller than 0.45, the BEs are underestimated and for larger α values increasingly overestimated. For the non-relativistic results, we observe a stronger species dependence of the optimal α value. As the non-relativistic Hamiltonian underestimates the core-level BE, increasing the exact exchange reduces the screening, resulting in a larger BE. An increase in α can thus offset the relativistic error. Comparing the MAE from non-relativistic G0W0@PBEh(α) calculations [Fig. 3(a)], we find that the optimal α indeed increases with the atomic number, from 0.44 for C1s excitations to about 0.53 for F1s.

We investigate three simple schemes to account for relativistic corrections (RC) in 1s core-level BEs from GW. Additionally, we discuss the application of these three schemes to ΔSCF. The first approach is to add the atomic corrections Δ ϵ 1 s , a t S K S defined in Eq. (19) to the QP energies and corresponds to the scheme we employed in our latest work.20 We label scheme I “method + RC.” The second is to use aZORA for the underlying DFT calculations and use the aZORA eigenvalues and MOs as the starting point for the GW calculation. We refer to scheme II as “method + aZORA.” In the third scheme, we use scheme II to obtain the QP energies and add the atomic corrections Δ ϵ 1 s , a t a Z O R A defined in Eq. (20) afterward. We label scheme III “method + aZORA + RC.”

For evGW0, we also explored a variant of schemes I and III, where we added the atomic corrections to the DFT eigenvalues instead to the QP energies. These corrected eigenvalues were then used as the starting point for the evGW0 calculation. These pre-correction variants yield with a mean absolute difference of 16 meV (pre-correction variant of scheme I) and 2.2 meV (pre-correction variant of scheme III) BEs that are extremely similar to the ones from evGW0 + RC and evGW0 + aZORA + RC, respectively; see Table S5 and Fig. S2 of the supplementary material for details. Adding the atomic correction as a post-processing step is transferable to non-relativistic GW results obtained from any code, and we thus disregard the pre-correction variant in the following.

Compared with the non-relativistic energies, the evGW0 + RC scheme reduces the error with respect to the experiment, as shown in Fig. 4(b). The errors are more tightly distributed, and the clustering by species is no longer evident. Generally, the BEs are still underestimated. However, the overall MAE is reduced from 0.55 eV to 0.30 eV and is now well within the accuracy required for the chemical analysis. Furthermore, the species dependence in the MAE is largely eliminated; see Fig. 3(c) and Table II. Solely the MAE for the F1s excitations is with 0.44 eV slightly larger than for the other elements. This might be attributable to poor statistics since our benchmark set contains only three F1s excitation.

TABLE II.

Mean absolute error (MAE) and mean error (ME) in eV with respect to the experiment, by species and in aggregate for absolute BEs of the CORE65 benchmark set. The error for excitation i is defined as E r r o r i = BE i theory BE i exp .

ΔSCF ΔSCF evGW0 evGW0
ΔSCF ΔSCF + RC + aZORA + aZORA + RC evGW0 evGW0 + RC + aZORA + aZORA + RC
Core-level MAE ME MAE ME MAE ME MAE ME MAE ME MAE ME MAE ME MAE ME
All  0.83  −0.83  0.57  −0.57  0.33  −0.31  0.46  −0.46  0.55  −0.55  0.30  −0.29  0.18  −0.07  0.32  −0.29 
C1s  0.60  −0.60  0.48  −0.48  0.36  −0.36  0.52  −0.52  0.39  −0.39  0.27  −0.27  0.21  −0.15  0.27  −0.25 
N1s  0.79  −0.79  0.56  −0.56  0.32  −0.32  0.52  −0.52  0.54  −0.54  0.30  −0.30  0.10  −0.07  0.27  −0.27 
O1s  1.11  −1.11  0.69  −0.68  0.32  −0.27  0.64  −0.63  0.71  −0.71  0.32  −0.28  0.19  0.02  0.38  −0.34 
F1s  1.36  −1.36  0.65  −0.65  0.12  0.03  0.57  −0.57  1.15  −1.15  0.44  −0.44  0.10  0.08  0.52  −0.52 
ΔSCF ΔSCF evGW0 evGW0
ΔSCF ΔSCF + RC + aZORA + aZORA + RC evGW0 evGW0 + RC + aZORA + aZORA + RC
Core-level MAE ME MAE ME MAE ME MAE ME MAE ME MAE ME MAE ME MAE ME
All  0.83  −0.83  0.57  −0.57  0.33  −0.31  0.46  −0.46  0.55  −0.55  0.30  −0.29  0.18  −0.07  0.32  −0.29 
C1s  0.60  −0.60  0.48  −0.48  0.36  −0.36  0.52  −0.52  0.39  −0.39  0.27  −0.27  0.21  −0.15  0.27  −0.25 
N1s  0.79  −0.79  0.56  −0.56  0.32  −0.32  0.52  −0.52  0.54  −0.54  0.30  −0.30  0.10  −0.07  0.27  −0.27 
O1s  1.11  −1.11  0.69  −0.68  0.32  −0.27  0.64  −0.63  0.71  −0.71  0.32  −0.28  0.19  0.02  0.38  −0.34 
F1s  1.36  −1.36  0.65  −0.65  0.12  0.03  0.57  −0.57  1.15  −1.15  0.44  −0.44  0.10  0.08  0.52  −0.52 

Scheme I has also been successfully employed for G0W0@PBEh. The range of optimal α is reduced by a factor of two for the G0W0@PBEh(α) + RC scheme vis-á-vis the non-relativistic one; see Figs. 3(a) and 3(b). With the relativistic correction, the value of α that minimizes the MAE ranges from 0.44 for C1s excitations to 0.49 for F1s excitations. This also shows in a slight species dependence of the MAE values we reported for the G0W0@PBEh(α = 0.45) results with RC earlier.20 

Judging by the MAE alone (Table II), the evGW0 + aZORA results are an improvement over the evGW0 + RC scheme. The overall MAE is 0.18 eV. Their distribution [Fig. 4(c)] is more centered. A slight clustering by species is observed, although this is not as obvious as for the non-relativistic results shown in Fig. 4(a). In contrast to the non-relativistic values, the 1s excitations of the lighter elements, such as carbon, tend to be more underestimated than the 1s BEs of oxygen.

The evGW0 + aZORA + RC scheme performs worse than the evGW0 + aZORA approach, with an error distribution similar to evGW0 + RC; see Fig. 4(d). The MAE for the individual species is in the same range as for evGW0 + RC, as is the overall MAE with 0.32 eV; see Table II. The evGW0 + aZORA + RC scheme is the most sophisticated among the three relativistic corrections discussed here: scalar-relativistic effects are included in the MOs used as the starting point for the evGW0 calculation, and the QP energies are corrected with respect to the fully relativistic atomic reference. It is thus surprising that it performs worse than evGW0 + aZORA. The similar performance of evGW0 + RC and evGW0 + aZORA + RC rather implies that the effect of including relativistic effects in the MOs is minimal.

To further investigate this surprising behavior, we visualized the mean errors (MEs) for C1s, N1s, O1s, and F1s in Fig. 5(b). The error is defined as BEtheory − BEexperiment. Thus, negative MEs indicate a systematic underestimation of the BEs with respect to the experiment. For the non-relativistic results, the absolute value of the ME corresponds directly to the MAE, and the ME is increasingly negative with the atomic number. The MEs for evGW0 + RC and evGW0 + aZORA + RC are negative and show almost no species dependence, which is in agreement with our previous analysis of the MAE and the error distribution in Fig. 4. For evGW0 + aZORA, however, we observe a trend that is reverse to the non-relativistic results. The ME increases with the atomic number and becomes even positive for F1s.

FIG. 5.

[(a) and (b)] Mean error with respect to the experiment for 1s BEs from (a) ΔSCF and (b) evGW0 for the CORE65 benchmark set, where the error for excitation i is defined as E r r o r i = BE i theory BE i exp . Three relativistic schemes are compared to the non-relativistic results. (c) Average size of the relativistic correction for 1s BEs from evGW0 for the CORE65 benchmark set with five additional molecules containing third period elements, where ΔBE = BErelativistic − BEnon-relativistic. For evGW0 + RC, ΔBE corresponds to the negative of the atomic corrections Δ ϵ 1 s , a t S K S , given in Table I.

FIG. 5.

[(a) and (b)] Mean error with respect to the experiment for 1s BEs from (a) ΔSCF and (b) evGW0 for the CORE65 benchmark set, where the error for excitation i is defined as E r r o r i = BE i theory BE i exp . Three relativistic schemes are compared to the non-relativistic results. (c) Average size of the relativistic correction for 1s BEs from evGW0 for the CORE65 benchmark set with five additional molecules containing third period elements, where ΔBE = BErelativistic − BEnon-relativistic. For evGW0 + RC, ΔBE corresponds to the negative of the atomic corrections Δ ϵ 1 s , a t S K S , given in Table I.

Close modal

Comparing non-relativistic with the relativistic BEs, we find that the size of the relativistic correction is 2–3 times larger with evGW0 + aZORA than with the other two schemes; see Fig. 5(c). In combination with the upward trend observed for the ME, this implies that aZORA is overestimating the relativistic correction for 1s states. This reflects the well-known tendency of aZORA to overestimate the relativistic correction to core-state eigenvalues.64,65 For the second period elements, where the relativistic error ranges from 0.2 eV to 0.7 eV, the aZORA overcorrection compensates, in part, a chronic underestimation of the BEs. While this error cancellation may seem fortuitous for second row elements, the rapid growth of the relativistic correction with the atomic number implies that evGW0 + aZORA might lead to large errors for 1s BEs of heavier elements. To illustrate this, we analyze the relativistic corrections for a small number of phosphorus and sulfur containing small molecules alongside those for the CORE65 benchmark set: H2S, SO2, PH3, PF3, and PF5 [see Fig. 5(c)]. The relativistic corrections obtained with the evGW0 + aZORA scheme are more than 10 eV larger than with evGW0 + RC and evGW0 + aZORA + RC. Also the difference between evGW0 + aZORA + RC and evGW0 + RC, which is negligible for second period elements, becomes more significant: 2.1 eV for the P1s excitations and 2.6 eV for the S1s. This suggests that the use of a scalar-relativistic reference for the underlying DFT calculation becomes more relevant as the magnitude of relativistic effects increases. However, the effect of the relativistic reference is only about one third of the magnitude of the relativistic correction for these states.

We also applied the three correction schemes to 1s BEs obtained from ΔSCF and plotted the MEs by species in Fig. 5(a). We observe the same trends as for evGW0. Note the similarity between Figs. 5(a) and 5(b). ΔSCF + RC or ΔSCF + aZORA + RC largely eliminate the species dependence of the ME. Although it is not as marked as for evGW0 + aZORA, the MEs also increase slightly with the atomic number for ΔSCF + aZORA. The size of the relativistic correction is two times larger with ΔSCF + aZORA than with ΔSCF + RC or ΔSCF + aZORA + RC, which suggests that the relativistic 1s corrections are also overestimated at the ΔSCF + aZORA level. ΔSCF entails the difference in total energy between the ionized and neutral systems. However, the ZORA overcorrection of core states propagates to the final calculation of the BE in approximately the same manner as in the 1s levels. This can be seen from Slater transition state theory:130,131 In ΔSCF, we calculate the energy difference ℏω = E(N − 1) − E(N), where N is the number of electrons and E(N − 1) and E(N) are the total energies of the core-ionized and neutral systems, respectively. This energy difference is approximately related to the 1s eigenvalue by ω = E / n 1 s + O ( ( δ n 1 s ) 3 ) and ∂E/∂n1s = ϵ1s, where n1s is the occupation number of the 1s state.132 This is an important detail to consider when comparing the performance of XC functionals for ΔSCF calculations of 1s excitations since the relativistic treatment makes a difference. For example, a recent study10 with the SCAN functional uses the ΔSCF approach in combination with scaled ZORA, while an atomic correction scheme was used for a similar benchmark study9 with the Tao-Perdew-Staroverov-Scuseria (TPSS) functional. Both studies report very good agreement with the experiment. However, it is difficult to judge, which functional performs better, since the relativistic effects are not treated on equal footing.

Both schemes that consistently improve the agreement with the experiment, evGW0 + RC and evGW0 + aZORA + RC, chronically underestimate the 1s BEs. This might be attributed to the broadening of the experimental spectra due to vibration effects, while GW yields vertical excitation energies. It has been demonstrated for GW-computed excitations of frontier orbitals that the deviation to the experiment can be reduced to 0.1 eV when fully resolving the vibrational structure based on the Franck–Condon multimode analysis, performed as a post-processing step to the GW calculation.133 This level of accuracy is often not required to resolve most XPS spectra but could be, in principle, reached by applying the same approach.

Relativistic corrections for 1s core-level energies from GW have been derived for non-relativistic and scalar-relativistic starting points. We have investigated three schemes for 1s QP energies from evGW0: A post-GW atomic correction (evGW0 + RC) using a non-relativistic reference, employing an aZORA reference (evGW0 + aZORA) and an aZORA reference along with a post-GW atomic correction (evGW0 + aZORA + RC). All three schemes improve agreement with the experiment. The evGW0 + RC and evGW0 + RC + aZORA schemes reduce the mean absolute error to about 0.3 eV and eliminate the species dependence. The evGW0 + aZORA scheme further reduces the overall MAE to 0.2 eV, but does so inconsistently, due to a species-dependent overcorrection.

The similarity of the results for the evGW0 + RC and evGW0 + aZORA + RC schemes indicates that the use of the scalar-relativistic reference has no significant effect on the result. Of the two, the evGW0 + RC scheme offers the further advantage that it is readily applicable in codes that have not implemented the aZORA Hamiltonian, and this is our recommendation to the practitioner. The application of the scheme entails an a posteriori shift of the non-relativistic evGW0 energies by the Δ ϵ 1 s , a t S K S values in Table I. We have shown that the derived corrections for the non-relativistic reference also consistently improve core-level energies from G0W0 and ΔSCF, suggesting that they are generally applicable to core-level BEs from different theoretical methods.

The evGW0 + RC and evGW0 + aZORA + RC schemes correct the non-relativistic values in a consistent manner and, therefore, form a solid foundation for further development and refinement. Further refinements may include the inclusion of vibrational effects to further improve the accuracy, in particular for C1s excitations, which are generally subject to very small chemical shifts. The development of relativistically corrected GW also paves the way for the accurate calculation of XPS in condensed phase systems, where Δ-approaches are problematic. Relativistic corrections will also improve the accuracy of x-ray absorption spectra from the Bethe–Salpeter equation, which employs the quantities computed in the GW calculations.134,135

See the supplementary material for the comparison of 1s eigenvalues for spherical and non-spherical solutions with the KS equations for neutral atoms (Table S1); plot of basis set dependence of the extrapolation scheme (Fig. S1); convergence of KS 1s eigenvalues and total energies for cc-pVnZ and NAO-VCC-nZ basis set series (Table S2); scalar-relativistic evGW0@PBE + aZORA results for basis set series cc-pVnZ (n = 3–6), extrapolated values, standard errors, and correlation coefficients (Table S3); results and experimental data for the CORE65 benchmark set for evGW0 + aZORA and evGW0 + aZORA + RC (Table S4); non-relativistic ΔSCF results and pre-corrected evGW0 energies (Table S5); difference between the pre- and post-corrected schemes for the CORE65 benchmark set (Fig. S2); and tabulated values of Fig. 2(b) (Table S6).

In pursuit of open materials science,136 we made the results of all relevant calculations available on the Novel Materials Discovery (NOMAD) repository.137,138

The authors would like to thank the CSC–IT Center for Science for providing computational resources. D.G. acknowledges financial support from the Academy of Finland (Grant No. 316168).

1.
P. S.
Bagus
,
E. S.
Ilton
, and
C. J.
Nelin
, “
The interpretation of XPS spectra: Insights into materials properties
,”
Surf. Sci. Rep.
68
,
273
304
(
2013
).
2.
P. S.
Bagus
,
F.
Illas
,
G.
Pacchioni
, and
F.
Parmigiani
, “
Mechanisms responsible for chemical shifts of core-level binding energies and their relationship to chemical bonding
,”
J. Electron Spectrosc. Relat. Phenom.
100
,
215
236
(
1999
).
3.
K.
Siegbahn
,
C.
Nordling
,
G.
Johansson
,
J.
Hedman
,
P. F.
Hedén
,
K.
Hamrin
,
U.
Gelius
,
T.
Bergmark
,
L. O.
Werme
,
R.
Manne
, and
Y.
Baer
,
ESCA Applied to Free Molecules
(
North-Holland Publishing Company
,
Amsterdam; London
,
1969
), pp.
51
136
.
4.
J. J.
Pireaux
,
S.
Svensson
,
E.
Basilier
,
P.-Å.
Malmqvist
,
U.
Gelius
,
R.
Caudano
, and
K.
Siegbahn
, “
Core-electron relaxation energies and valence-band formation of linear alkanes studied in the gas phase by means of electron spectroscopy
,”
Phys. Rev. A
14
,
2133
2145
(
1976
).
5.
A.
Aarva
,
V. L.
Deringer
,
S.
Sainio
,
T.
Laurila
, and
M. A.
Caro
, “
Understanding x-ray spectroscopy of carbonaceous materials by combining experiments, density functional theory and machine learning. Part I: Fingerprint spectra
,”
Chem. Mater.
31
,
9243
9255
(
2019
).
6.
A.
Aarva
,
V. L.
Deringer
,
S.
Sainio
,
T.
Laurila
, and
M. A.
Caro
, “
Understanding x-ray spectroscopy of carbonaceous materials by combining experiments, density functional theory and machine learning. Part II: Quantitative fitting of spectra
,”
Chem. Mater.
31
,
9256
9267
(
2019
).
7.
F.
Viñes
,
C.
Sousa
, and
F.
Illas
, “
On the prediction of core level binding energies in molecules, surfaces and solids
,”
Phys. Chem. Chem. Phys.
20
,
8403
8410
(
2018
).
8.
P. S.
Bagus
, “
Self-consistent-field wave functions for hole states of some Ne-like and Ar-like ions
,”
Phys. Rev.
139
,
A619
A634
(
1965
).
9.
N.
Pueyo Bellafont
,
F.
Viñes
, and
F.
Illas
, “
Performance of the TPSS functional on predicting core level binding energies of main group elements containing molecules: A good choice for molecules adsorbed on metal surfaces
,”
J. Chem. Theory Comput.
12
,
324
331
(
2016
).
10.
J. M.
Kahk
and
J.
Lischner
, “
Accurate absolute core-electron binding energies of molecules, solids, and surfaces from first-principles calculations
,”
Phys. Rev. Mater.
3
,
100801
(
2019
).
11.
X.
Zheng
and
L.
Cheng
, “
Performance of delta-coupled-cluster methods for calculations of core-ionization energies of first-row elements
,”
J. Chem. Theory Comput.
15
,
4945
4955
(
2019
).
12.
A.
Holme
,
K. J.
Børve
,
L. J.
Sæthre
, and
T. D.
Thomas
, “
Accuracy of calculated chemical shifts in carbon 1s ionization energies from single-reference ab initio methods and density functional theory
,”
J. Chem. Theory Comput.
7
,
4104
4114
(
2011
).
13.
S.
Sen
,
A.
Shee
, and
D.
Mukherjee
, “
Inclusion of orbital relaxation and correlation through the unitary group adapted open shell coupled cluster theory using non-relativistic and scalar relativistic Hamiltonians to study the core ionization potential of molecules containing light to medium-heavy elements
,”
J. Chem. Phys.
148
,
054107
(
2018
).
14.
G. S.
Michelitsch
and
K.
Reuter
, “
Efficient simulation of near-edge x-ray absorption fine structure (NEXAFS) in density-functional theory: Comparison of core-level constraining approaches
,”
J. Chem. Phys.
150
,
074104
(
2019
).
15.
P. S.
Bagus
,
C. J.
Nelin
,
X.
Zhao
,
S. V.
Levchenko
,
E.
Davis
,
X.
Weng
,
F.
Späth
,
C.
Papp
,
H.
Kuhlenbeck
, and
H.-J.
Freund
, “
Revisiting surface core-level shifts for ionic compounds
,”
Phys. Rev. B
100
,
115419
(
2019
).
16.
N.
Pueyo Bellafont
,
F.
Viñes
,
W.
Hieringer
, and
F.
Illas
, “
Predicting core level binding energies shifts: Suitability of the projector augmented wave approach as implemented in VASP
,”
J. Comput. Chem.
38
,
518
522
(
2017
).
17.
L.
Köhler
and
G.
Kresse
, “
Density functional study of CO on Rh(111)
,”
Phys. Rev. B
70
,
165405
(
2004
).
18.
J.
Liu
,
D.
Matthews
,
S.
Coriani
, and
L.
Cheng
, “
Benchmark calculations of K-edge ionization energies for first-row elements using scalar-relativistic core-valence-separated equation-of-motion coupled-cluster methods
,”
J. Chem. Theory Comput.
15
,
1642
1651
(
2019
).
19.
V. K.
Voora
,
R.
Galhenage
,
J. C.
Hemminger
, and
F.
Furche
, “
Effective one-particle energies from generalized Kohn–Sham random phase approximation: A direct approach for computing and analyzing core ionization energies
,”
J. Chem. Phys.
151
,
134106
(
2019
).
20.
D.
Golze
,
L.
Keller
, and
P.
Rinke
, “
Accurate absolute and relative core-level binding energies from GW
,”
J. Phys. Chem. Lett.
11
,
1840
1847
(
2020
).
21.
L.
Hedin
, “
New method for calculating the one-particle Green’s function with application to the electron-gas problem
,”
Phys. Rev.
139
,
A796
A823
(
1965
).
22.
D.
Golze
,
M.
Dvorak
, and
P.
Rinke
, “
The GW compendium: A practical guide to theoretical photoemission spectroscopy
,”
Front. Chem.
7
,
377
(
2019
).
23.
X.
Blase
,
C.
Attaccalite
, and
V.
Olevano
, “
First-principles GW calculations for fullerenes, porphyrins, phtalocyanine, and other molecules of interest for organic photovoltaic applications
,”
Phys. Rev. B
83
,
115103
(
2011
).
24.
X.
Ren
,
P.
Rinke
,
V.
Blum
,
J.
Wieferink
,
A.
Tkatchenko
,
A.
Sanfilippo
,
K.
Reuter
, and
M.
Scheffler
, “
Resolution-of-identity approach to Hartree-Fock, hybrid density functionals, RPA, MP2 and GW with numeric atom-centered orbital basis functions
,”
New J. Phys.
14
,
053020
(
2012
).
25.
M. J.
van Setten
,
F.
Weigend
, and
F.
Evers
, “
The GW-method for quantum chemistry applications: Theory and implementation
,”
J. Chem. Theory Comput.
9
,
232
246
(
2013
).
26.
F.
Bruneval
,
T.
Rangel
,
S. M.
Hamed
,
M.
Shao
,
C.
Yang
, and
J. B.
Neaton
, “
MOLGW1: Many-body perturbation theory software for atoms, molecules, and clusters
,”
Comput. Phys. Commun.
208
,
149
161
(
2016
).
27.
J.
Wilhelm
,
M.
Del Ben
, and
J.
Hutter
, “
GW in the Gaussian and plane waves scheme with application to linear acenes
,”
J. Chem. Theory Comput.
12
,
3623
3635
(
2016
).
28.
J.
Wilhelm
and
J.
Hutter
, “
Periodic GW calculations in the Gaussian and plane-waves scheme
,”
Phys. Rev. B
95
,
235123
(
2017
).
29.
J.
Wilhelm
,
D.
Golze
,
L.
Talirz
,
J.
Hutter
, and
C. A.
Pignedoli
, “
Toward GW calculations on thousands of atoms
,”
J. Phys. Chem. Lett.
9
,
306
312
(
2018
).
30.
J.
Wilhelm
,
J.
VandeVondele
, and
V. V.
Rybkin
, “
Dynamics of the bulk hydrated electron from many-body wave-function theory
,”
Angew. Chem., Int. Ed.
58
,
3890
3893
(
2019
).
31.
A.
Stuke
,
C.
Kunkel
,
D.
Golze
,
M.
Todorović
,
J. T.
Margraf
,
K.
Reuter
,
P.
Rinke
, and
H.
Oberhofer
, “
Atomic structures and orbital energies of 61 489 crystal-forming organic molecules
,”
Sci. Data
7
,
58
(
2020
).
32.
M.
Guzzo
,
G.
Lani
,
F.
Sottile
,
P.
Romaniello
,
M.
Gatti
,
J. J.
Kas
,
J. J.
Rehr
,
M. G.
Silly
,
F.
Sirotti
, and
L.
Reining
, “
Valence electron photoemission spectrum of semiconductors: Ab initio description of multiple satellites
,”
Phys. Rev. Lett.
107
,
166401
(
2011
).
33.
J. S.
Zhou
,
J. J.
Kas
,
L.
Sponza
,
I.
Reshetnyak
,
M.
Guzzo
,
C.
Giorgetti
,
M.
Gatti
,
F.
Sottile
,
J. J.
Rehr
, and
L.
Reining
, “
Dynamical effects in electron spectroscopy
,”
J. Chem. Phys.
143
,
184109
(
2015
).
34.
J. J.
Kas
,
J. J.
Rehr
, and
J. B.
Curtis
, “
Particle-hole cumulant approach for inelastic losses in x-ray spectra
,”
Phys. Rev. B
94
,
035156
(
2016
).
35.
S.
Ishii
,
S.
Iwata
, and
K.
Ohno
, “
All-electron GW calculations of silicon, diamond, and silicon carbide
,”
Mater. Trans.
51
,
2150
2156
(
2010
).
36.
T.
Aoki
and
K.
Ohno
, “
Accurate quasiparticle calculation of x-ray photoelectron spectra of solids
,”
J. Phys.: Condens. Matter
30
,
21LT01
(
2018
).
37.
M. J.
van Setten
,
R.
Costa
,
F.
Viñes
, and
F.
Illas
, “
Assessing GW approaches for predicting core level binding energies
,”
J. Chem. Theory Comput.
14
,
877
883
(
2018
).
38.
D.
Golze
,
J.
Wilhelm
,
M. J.
van Setten
, and
P.
Rinke
, “
Core-level binding energies from GW: An efficient full-frequency approach within a localized basis
,”
J. Chem. Theory Comput.
14
,
4856
4869
(
2018
).
39.
F.
Aryasetiawan
and
S.
Biermann
, “
Generalized Hedin’s equations for quantum many-body systems with spin-dependent interactions
,”
Phys. Rev. Lett.
100
,
116402
(
2008
).
40.
F.
Aryasetiawan
and
S.
Biermann
, “
Generalized Hedin equations and σGσW approximation for quantum many-body systems with spin-dependent interactions
,”
J. Phys.: Condens. Matter
21
,
064232
(
2009
).
41.
M.
Kühn
and
F.
Weigend
, “
One-electron energies from the two-component GW method
,”
J. Chem. Theory Comput.
11
,
969
979
(
2015
).
42.
C.
Holzer
and
W.
Klopper
, “
Ionized, electron-attached, and excited states of molecular systems with spin-orbit coupling: Two-component GW and Bethe-Salpeter implementations
,”
J. Chem. Phys.
150
,
204116
(
2019
).
43.
R.
Sakuma
,
C.
Friedrich
,
T.
Miyake
,
S.
Blügel
, and
F.
Aryasetiawan
, “
GW calculations including spin-orbit coupling: Application to Hg chalcogenides
,”
Phys. Rev. B
84
,
085144
(
2011
).
44.
P.
Scherpelz
,
M.
Govoni
,
I.
Hamada
, and
G.
Galli
, “
Implementation and validation of fully relativistic GW calculations: Spin-orbit coupling in molecules, nanocrystals, and solids
,”
J. Chem. Theory Comput.
12
,
3523
3544
(
2016
).
45.
B.
Schuler
,
D. Y.
Qiu
,
S.
Refaely-Abramson
,
C.
Kastl
,
C. T.
Chen
,
S.
Barja
,
R. J.
Koch
,
D. F.
Ogletree
,
S.
Aloni
,
A. M.
Schwartzberg
,
J. B.
Neaton
,
S. G.
Louie
, and
A.
Weber-Bargioni
, “
Large spin-orbit splitting of deep in-gap defect states of engineered sulfur vacancies in monolayer WS2
,”
Phys. Rev. Lett.
123
,
076801
(
2019
).
46.
A.
Kutepov
,
K.
Haule
,
S. Y.
Savrasov
, and
G.
Kotliar
, “
Electronic structure of Pu and Am metals by self-consistent relativistic GW method
,”
Phys. Rev. B
85
,
155129
(
2012
).
47.
T.
Ahmed
,
R. C.
Albers
,
A. V.
Balatsky
,
C.
Friedrich
, and
J.-X.
Zhu
, “
GW quasiparticle calculations with spin-orbit coupling for the light actinides
,”
Phys. Rev. B
89
,
035104
(
2014
).
48.
A.
Molina-Sánchez
,
D.
Sangalli
,
K.
Hummer
,
A.
Marini
, and
L.
Wirtz
, “
Effect of spin-orbit interaction on the optical spectra of single-layer, double-layer, and bulk MoS2
,”
Phys. Rev. B
88
,
045412
(
2013
).
49.
L.
Guo
,
M.
Wu
,
T.
Cao
,
D. M.
Monahan
,
Y.-H.
Lee
,
S. G.
Louie
, and
G. R.
Fleming
, “
Exchange-driven intravalley mixing of excitons in monolayer transition metal dichalcogenides
,”
Nat. Phys.
15
,
228
232
(
2019
).
50.
P.
Umari
,
E.
Mosconi
, and
F.
De Angelis
, “
Relativistic GW calculations on CH3 NH3 PbI3 and CH3 NH3 SnI3 perovskites for solar cell applications
,”
Sci. Rep.
4
,
4467
(
2014
).
51.
I. A.
Nechaev
and
E. V.
Chulkov
, “
Quasiparticle band gap in the topological insulator Bi2Te3
,”
Phys. Rev. B
88
,
165135
(
2013
).
52.
I.
Aguilera
,
C.
Friedrich
,
G.
Bihlmayer
, and
S.
Blügel
, “
GW study of topological insulators Bi2Se3, Bi2Te3, and Sb2Te3: Beyond the perturbative one-shot approach
,”
Phys. Rev. B
88
,
045206
(
2013
).
53.
I.
Aguilera
,
C.
Friedrich
, and
S.
Blügel
, “
Spin-orbit coupling in quasiparticle studies of topological insulators
,”
Phys. Rev. B
88
,
165136
(
2013
).
54.
I.
Aguilera
,
C.
Friedrich
, and
S.
Blügel
, “
Electronic phase transitions of bismuth under strain from relativistic self-consistent GW calculations
,”
Phys. Rev. B
91
,
125129
(
2015
).
55.
I. A.
Nechaev
,
I.
Aguilera
,
V.
De Renzi
,
A.
Di Bona
,
A.
Lodi Rizzini
,
A. M.
Mio
,
G.
Nicotra
,
A.
Politano
,
S.
Scalese
,
Z. S.
Aliev
,
M. B.
Babanly
,
C.
Friedrich
,
S.
Blügel
, and
E. V.
Chulkov
, “
Quasiparticle spectrum and plasmonic excitations in the topological insulator Sb2Te3
,”
Phys. Rev. B
91
,
245123
(
2015
).
56.
M.
Battiato
,
I.
Aguilera
, and
J.
Sánchez-Barriga
, “
Generalized GW + Boltzmann approach for the description of ultrafast electron dynamics in topological insulators
,”
Materials
10
,
810
(
2017
).
57.
M.
Rohlfing
and
S. G.
Louie
, “
Quasiparticle band structure of HgSe
,”
Phys. Rev. B
57
,
R9392
R9395
(
1998
).
58.
W.
Gao
,
W.
Xia
,
Y.
Wu
,
W.
Ren
,
X.
Gao
, and
P.
Zhang
, “
Quasiparticle band structures of CuCl, CuBr, AgCl, and AgBr: The extreme case
,”
Phys. Rev. B
98
,
045108
(
2018
).
59.
K.
Bushick
,
K.
Mengle
,
N.
Sanders
, and
E.
Kioupakis
, “
Band structure and carrier effective masses of boron arsenide: Effects of quasiparticle and spin-orbit coupling corrections
,”
Appl. Phys. Lett.
114
,
022101
(
2019
).
60.
M.-Y.
Zhang
and
H.
Jiang
, “
Electronic band structure of cuprous and silver halides: An all-electron GW study
,”
Phys. Rev. B
100
,
205123
(
2019
).
61.
See https://xdb.lbl.gov/Section1 for x-ray data booklet (accessed on 12 May 2020).
62.
V.
Blum
,
R.
Gehrke
,
F.
Hanke
,
P.
Havu
,
V.
Havu
,
X.
Ren
,
K.
Reuter
, and
M.
Scheffler
, “
Ab initio molecular simulations with numeric atom-centered orbitals
,”
Comput. Phys. Commun.
180
,
2175
2196
(
2009
).
63.
W. P.
Huhn
and
V.
Blum
, “
One-hundred-three compound band-structure benchmark of post-self-consistent spin-orbit coupling treatments in density functional theory
,”
Phys. Rev. Mater.
1
,
033803
(
2017
).
64.
R.
van Leeuwen
,
E.
van Lenthe
,
E. J.
Baerends
, and
J. G.
Snijders
, “
Exact solutions of regular approximate relativistic wave equations for hydrogen-like atoms
,”
J. Chem. Phys.
101
,
1272
1281
(
1994
).
65.
D.
Sundholm
, “
Perturbation theory based on quasi-relativistic Hamiltonians
,” in
Relativistic Electronic Structure Theory
, edited by
P.
Schwerdtfeger
(
Elsevier
,
2002
), Chap. 13, pp.
758
792
.
66.
S. L.
Adler
, “
Quantum theory of the dielectric constant in real solids
,”
Phys. Rev.
126
,
413
420
(
1962
).
67.
N.
Wiser
, “
Dielectric constant with local field effects included
,”
Phys. Rev.
129
,
62
69
(
1963
).
68.
P.
Liu
,
M.
Kaltak
,
J.
Klimeš
, and
G.
Kresse
, “
Cubic scaling GW: Towards fast quasiparticle calculations
,”
Phys. Rev. B
94
,
165109
(
2016
).
69.
M. J.
van Setten
,
F.
Caruso
,
S.
Sharifzadeh
,
X.
Ren
,
M.
Scheffler
,
F.
Liu
,
J.
Lischner
,
L.
Lin
,
J. R.
Deslippe
,
S. G.
Louie
,
C.
Yang
,
F.
Weigend
,
J. B.
Neaton
,
F.
Evers
, and
P.
Rinke
, “
GW100: Benchmarking G0W0 for molecular systems
,”
J. Chem. Theory Comput.
11
,
5665
5687
(
2015
).
70.
R. W.
Godby
,
M.
Schlüter
, and
L. J.
Sham
, “
Self-energy operators and exchange-correlation potentials in semiconductors
,”
Phys. Rev. B
37
,
10159
10175
(
1988
).
71.
X.
Gonze
,
B.
Amadon
,
P.-M.
Anglade
,
J.-M.
Beuken
,
F.
Bottin
,
P.
Boulanger
,
F.
Bruneval
,
D.
Caliste
,
R.
Caracas
,
M.
Côté
,
T.
Deutsch
,
L.
Genovese
,
P.
Ghosez
,
M.
Giantomassi
,
S.
Goedecker
,
D. R.
Hamann
,
P.
Hermet
,
F.
Jollet
,
G.
Jomard
,
S.
Leroux
,
M.
Mancini
,
S.
Mazevet
,
M. J. T.
Oliveira
,
G.
Onida
,
Y.
Pouillon
,
T.
Rangel
,
G.-M.
Rignanese
,
D.
Sangalli
,
R.
Shaltaf
,
M.
Torrent
,
M. J.
Verstraete
,
G.
Zerah
, and
J. W.
Zwanziger
, “
ABINIT: First-principles approach to material and nanosystem properties
,”
Comput. Phys. Commun.
180
,
2582
2615
(
2009
).
72.
M.
Govoni
and
G.
Galli
, “
Large scale GW calculations
,”
J. Chem. Theory Comput.
11
,
2680
2696
(
2015
).
73.
C.
Friedrich
, “
Tetrahedron integration method for strongly varying functions: Application to the gt self-energy
,”
Phys. Rev. B
100
,
075142
(
2019
).
74.
I.
Duchemin
and
X.
Blase
, “
Robust analytic-continuation approach to many-body GW calculations
,”
J. Chem. Theory Comput.
16
,
1742
1756
(
2020
).
75.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
, “
Generalized gradient approximation made simple
,”
Phys. Rev. Lett.
77
,
3865
3868
(
1996
).
76.
J.
Lischner
,
J.
Deslippe
,
M.
Jain
, and
S. G.
Louie
, “
First-principles calculations of quasiparticle excitations of open-shell condensed matter systems
,”
Phys. Rev. Lett.
109
,
036406
(
2012
).
77.
P.-F.
Loos
,
P.
Romaniello
, and
J. A.
Berger
, “
Green functions and self-consistency: Insights from the spherium model
,”
J. Chem. Theory Comput.
14
,
3071
3082
(
2018
).
78.
M.
Véril
,
P.
Romaniello
,
J. A.
Berger
, and
P.-F.
Loos
, “
Unphysical discontinuities in GW methods
,”
J. Chem. Theory Comput.
14
,
5220
5228
(
2018
).
79.
M.
Gatti
,
G.
Panaccione
, and
L.
Reining
, “
Effects of low-energy excitations on spectral properties at higher binding energy: The metal-insulator transition of VO2
,”
Phys. Rev. Lett.
114
,
116402
(
2015
).
80.
Y.-M.
Byun
and
S.
Öğüt
, “
Practical GW scheme for electronic structure of 3D-transition-metal monoxide anions: ScO, TiO, CuO, and ZnO
,”
J. Chem. Phys.
151
,
134305
(
2019
).
81.
R.
Sankari
,
M.
Ehara
,
H.
Nakatsuji
,
A.
De Fanis
,
H.
Aksela
,
S. L.
Sorensen
,
M. N.
Piancastelli
,
E.
Kukk
, and
K.
Ueda
, “
High resolution O 1s photoelectron shake-up satellite spectrum of H2O
,”
Chem. Phys. Lett.
422
,
51
57
(
2006
).
82.
J.
Schirmer
,
G.
Angonoa
,
S.
Svensson
,
D.
Nordfors
, and
U.
Gelius
, “
High-energy photoelectron C 1s and O 1s shake-up spectra of CO
,”
J. Phys. B: At., Mol. Opt. Phys.
20
,
6031
6040
(
1987
).
83.
M.
Shishkin
and
G.
Kresse
, “
Self-consistent GW calculations for semiconductors and insulators
,”
Phys. Rev. B
75
,
235102
(
2006
).
84.
N.
Marom
,
F.
Caruso
,
X.
Ren
,
O. T.
Hofmann
,
T.
Körzdörfer
,
J. R.
Chelikowsky
,
A.
Rubio
,
M.
Scheffler
, and
P.
Rinke
, “
Benchmark of GW methods for azabenzenes
,”
Phys. Rev. B
86
,
245127
(
2012
).
85.
W.-D.
Schöne
and
A. G.
Eguiluz
, “
Self-consistent calculations of quasiparticle states in metals and semiconductors
,”
Phys. Rev. Lett.
81
,
1662
1665
(
1998
).
86.
F.
Caruso
,
P.
Rinke
,
X.
Ren
,
M.
Scheffler
, and
A.
Rubio
, “
Unified description of ground and excited states of finite systems: The self-consistent GW approach
,”
Phys. Rev. B
86
,
081102(R)
(
2012
).
87.
F.
Caruso
,
P.
Rinke
,
X.
Ren
,
A.
Rubio
, and
M.
Scheffler
, “
Self-consistent GW: All-electron implementation with localized basis functions
,”
Phys. Rev. B
88
,
075105
(
2013
).
88.
M.
van Schilfgaarde
,
T.
Kotani
, and
S.
Faleev
, “
Quasiparticle self-consistent GW theory
,”
Phys. Rev. Lett.
96
,
226402
(
2006
).
89.
F.
Caruso
,
M.
Dauth
,
M. J.
van Setten
, and
P.
Rinke
, “
Benchmark of GW approaches for the GW100 test set
,”
J. Chem. Theory Comput.
12
,
5076
5087
(
2016
).
90.
M.
Grumet
,
P.
Liu
,
M.
Kaltak
,
J.
Klimeŝ
, and
G.
Kresse
, “
Beyond the quasiparticle approximation: Fully self-consistent GW calculations
,”
Phys. Rev. B
98
,
155143
(
2018
).
91.
V.
Atalla
,
M.
Yoon
,
F.
Caruso
,
P.
Rinke
, and
M.
Scheffler
, “
Hybrid density functional theory meets quasiparticle calculations: A consistent electronic structure approach
,”
Phys. Rev. B
88
,
165122
(
2013
).
92.
P.
Pyykkö
, “
Relativistic effects in structural chemistry
,”
Chem. Rev.
88
,
563
594
(
1988
).
93.
Relativistic Electronic Structure Theory
, edited by
P.
Schwerdtfeger
(
Elsevier
,
2002
).
94.
P. A. M.
Dirac
, “
The quantum theory of the electron
,”
Proc. R. Soc. London, Ser. A
117
,
610
624
(
1928
).
95.
G.
Breit
, “
The effect of retardation on the interaction of two electrons
,”
Phys. Rev.
34
,
553
573
(
1929
).
96.
W.
Kutzelnigg
, “
Perturbation theory of relativistic effects
,” in
Relativistic Electronic Structure Theory
, edited by
P.
Schwerdtfeger
(
Elsevier
,
2002
), Chap. 12, pp.
664
757
.
97.
L.
Visscher
and
K. G.
Dyall
, “
Relativistic and correlation effects on molecular properties. I. The dihalogens F2, Cl2, Br2, I2, and At2
,”
J. Chem. Phys.
104
,
9040
9046
(
1996
).
98.
L.
Visscher
,
J.
Styszyñski
, and
W. C.
Nieuwpoort
, “
Relativistic and correlation effects on molecular properties. II. The hydrogen halides HF, HCl, HBr, HI, and HAt
,”
J. Chem. Phys.
105
,
1987
1994
(
1996
).
99.
H. M.
Quiney
,
H.
Skaane
, and
I. P.
Grant
, “
Relativistic, quantum electrodynamic and many-body effects in the water molecule
,”
Chem. Phys. Lett.
290
,
473
480
(
1998
).
100.
O.
Fossgaard
,
O.
Gropen
,
M. C.
Valero
, and
T.
Saue
, “
On the performance of four-component relativistic density functional theory: Spectroscopic constants and dipole moments of the diatomics HX and XY (X, Y=F, Cl, Br, and I)
,”
J. Chem. Phys.
118
,
10418
10430
(
2003
).
101.
A. K.
Rajagopal
and
J.
Callaway
, “
Inhomogeneous electron gas
,”
Phys. Rev. B
7
,
1912
1919
(
1973
).
102.
A. K.
Rajagopal
, “
Inhomogeneous relativistic electron gas
,”
J. Phys. C
11
,
L943
L948
(
1978
).
103.
A. H.
MacDonald
and
S. H.
Vosko
, “
A relativistic density functional formalism
,”
J. Phys. C
12
,
2977
2990
(
1979
).
104.
T.
Saue
and
T.
Helgaker
, “
Four-component relativistic Kohn–Sham theory
,”
J. Comput. Chem.
23
,
814
823
(
2002
).
105.
M.
Reiher
and
A.
Wolf
,
Relativistic Quantum Chemistry: The Fundamental Theory of Molecular Science
(
Wiley
,
Weinheim
,
2009
).
106.
E.
Engel
, “
Relativistic density functional theory
,” in
Relativistic Electronic Structure Theory
, edited by
P.
Schwerdtfeger
(
Elsevier
,
2002
), Chap. 10, pp.
253
621
.
107.
M. V.
Ramana
and
A. K.
Rajagopal
, “
Inhomogeneous relativistic electron gas: Correlation potential
,”
Phys. Rev. A
24
,
1689
1695
(
1981
).
108.
K.
Lejaeghere
,
G.
Bihlmayer
,
T.
Björkman
,
P.
Blaha
,
S.
Blügel
,
V.
Blum
,
D.
Caliste
,
I. E.
Castelli
,
S. J.
Clark
,
A.
Dal Corso
,
S.
de Gironcoli
,
T.
Deutsch
,
J. K.
Dewhurst
,
I.
Di Marco
,
C.
Draxl
,
M.
Dułak
,
O.
Eriksson
,
J. A.
Flores-Livas
,
K. F.
Garrity
,
L.
Genovese
,
P.
Giannozzi
,
M.
Giantomassi
,
S.
Goedecker
,
X.
Gonze
,
O.
Grånäs
,
E. K. U.
Gross
,
A.
Gulans
,
F.
Gygi
,
D. R.
Hamann
,
P. J.
Hasnip
,
N. A. W.
Holzwarth
,
D.
Iuşan
,
D. B.
Jochym
,
F.
Jollet
,
D.
Jones
,
G.
Kresse
,
K.
Koepernik
,
E.
Küçükbenli
,
Y. O.
Kvashnin
,
I. L. M.
Locht
,
S.
Lubeck
,
M.
Marsman
,
N.
Marzari
,
U.
Nitzsche
,
L.
Nordström
,
T.
Ozaki
,
L.
Paulatto
,
C. J.
Pickard
,
W.
Poelmans
,
M. I. J.
Probert
,
K.
Refson
,
M.
Richter
,
G.-M.
Rignanese
,
S.
Saha
,
M.
Scheffler
,
M.
Schlipf
,
K.
Schwarz
,
S.
Sharma
,
F.
Tavazza
,
P.
Thunström
,
A.
Tkatchenko
,
M.
Torrent
,
D.
Vanderbilt
,
M. J.
van Setten
,
V.
Van Speybroeck
,
J. M.
Wills
,
J. R.
Yates
,
G.-X.
Zhang
, and
S.
Cottenier
, “
Reproducibility in density functional theory calculations of solids
,”
Science
351
,
aad3000
(
2016
).
109.
E. v.
Lenthe
,
E. J.
Baerends
, and
J. G.
Snijders
, “
Relativistic regular two-component Hamiltonians
,”
J. Chem. Phys.
99
,
4597
4610
(
1993
).
110.
E.
van Lenthe
,
E. J.
Baerends
, and
J. G.
Snijders
, “
Relativistic total energy using regular approximations
,”
J. Chem. Phys.
101
,
9783
9792
(
1994
).
111.
V.
Havu
,
V.
Blum
,
P.
Havu
, and
M.
Scheffler
, “
Efficient O(N) integration for all-electron electronic structure calculation using numeric basis functions
,”
J. Comput. Phys.
228
,
8367
8379
(
2009
).
112.
C.
Adamo
and
V.
Barone
, “
Toward reliable density functional methods without adjustable parameters: The PBE0 model
,”
J. Chem. Phys.
110
,
6158
6170
(
1999
).
113.
M.
Ernzerhof
and
G. E.
Scuseria
, “
Assessment of the Perdew-Burke-Ernzerhof exchange-correlation functional
,”
J. Chem. Phys.
110
,
5029
5036
(
1999
).
114.
F.
Weigend
and
R.
Ahlrichs
, “
Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy
,”
Phys. Chem. Chem. Phys.
7
,
3297
3305
(
2005
).
115.
F.
Bruneval
, “
Ionization energy of atoms obtained from GW self-energy or from random phase approximation total energies
,”
J. Chem. Phys.
136
,
194107
(
2012
).
116.
F.
Bruneval
and
M. A. L.
Marques
, “
Benchmarking the starting points of the GW approximation for molecules
,”
J. Chem. Theory Comput.
9
,
324
329
(
2013
).
117.
T. H.
Dunning
, “
Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen
,”
J. Chem. Phys.
90
,
1007
1023
(
1989
).
118.
A. K.
Wilson
,
T.
van Mourik
, and
T. H.
Dunning
, “
Gaussian basis sets for use in correlated molecular calculations. VI. Sextuple zeta correlation consistent basis sets for boron through neon
,”
J. Mol. Struc.-THEOCHEM
388
,
339
349
(
1996
).
119.
O.
Čertík
,
J. E.
Pask
, and
J.
Vackář
, “
DFTATOM: A robust and general Schrödinger and Dirac solver for atomic structure calculations
,”
Comput. Phys. Commun.
184
,
1777
1791
(
2013
).
120.
DIRAC, a relativistic ab initio electronic structure program, release DIRAC18 (2018), written by T. Saue, L. Visscher, H. J. Aa. Jensen, and R. Bast, with contributions from V. Bakken, K. G. Dyall, S. Dubillard, U. Ekström, E. Eliav, T. Enevoldsen, E. Faßhauer, T. Fleig, O. Fossgaard, A. S. P. Gomes, E. D. Hedegård, T. Helgaker, J. Henriksson, M. Iliaš, Ch. R. Jacob, S. Knecht, S. Komorovský, O. Kullie, J. K. Lærdahl, C. V. Larsen, Y. S. Lee, H. S. Nataraj, M. K. Nayak, P. Norman, G. Olejniczak, J. Olsen, J. M. H. Olsen, Y. C. Park, J. K. Pedersen, M. Pernpointner, R. di Remigio, K. Ruud, P. Sałek, B. Schimmelpfennig, A. Shee, J. Sikkema, A. J. Thorvaldsen, J. Thyssen, J. van Stralen, S. Villaume, O. Visser, T. Winther, and S. Yamamoto (available at , see also http://www.diracprogram.org).
121.
T.
Saue
,
R.
Bast
,
A. S. P.
Gomes
,
H. J. A.
Jensen
,
L.
Visscher
,
I. A.
Aucar
,
R.
Di Remigio
,
K. G.
Dyall
,
E.
Eliav
,
E.
Fasshauer
,
T.
Fleig
,
L.
Halbert
,
E. D.
Hedegård
,
B.
Helmich-Paris
,
M.
Iliaš
,
C. R.
Jacob
,
S.
Knecht
,
J. K.
Laerdahl
,
M. L.
Vidal
,
M. K.
Nayak
,
M.
Olejniczak
,
J. M. H.
Olsen
,
M.
Pernpointner
,
B.
Senjean
,
A.
Shee
,
A.
Sunaga
, and
J. N. P.
van Stralen
, “
The Dirac code for relativistic molecular calculations
,”
J. Chem. Phys.
152
,
204104
(
2020
).
122.
K. G.
Dyall
, “
Relativistic double-zeta, triple-zeta, and quadruple-zeta basis sets for the light elements H–Ar
,”
Theo. Chem. Acc.
135
,
128
(
2016
).
123.
W.
Gordon
, “
Über den stoß zweier punktladungen nach der wellenmechanik
,”
Z. Phys.
48
,
180
191
(
1928
).
124.
C. G.
Darwin
, “
The wave equations of the electron
,”
Proc. R. Soc. London, Ser. A
118
,
654
680
(
1928
).
125.
M.
Eides
,
H.
Grotch
, and
V.
Shelyuto
,
Theory of Light Hydrogenic Bound States
(
Springer
,
2007
), Vol. 222.
126.
C. L.
Pekeris
, “
Ground state of two-electron atoms
,”
Phys. Rev.
112
,
1649
1658
(
1958
).
127.
P. K.
Mukherjee
and
D. P.
Chong
, “
Ab initio calculation of core-electron binding energies in small molecules
,”
Chem. Phys. Lett.
120
,
163
166
(
1985
).
128.
C.
Chang
,
M.
Pelissier
, and
P.
Durand
, “
Regular two-component Pauli-like effective Hamiltonians in Dirac theory
,”
Phys. Scr.
34
,
394
404
(
1986
).
129.
L.
Triguero
,
O.
Plashkevych
,
L. G. M.
Pettersson
, and
H.
Ågren
, “
Separate state vs transition state Kohn-Sham calculations of x-ray photoelectron binding energies and chemical shifts
,”
J. Electron Spectrosc. Relat. Phenom.
104
,
195
207
(
1999
).
130.
J. C.
Slater
, “
A simplification of the Hartree-Fock method
,”
Phys. Rev.
81
,
385
390
(
1951
).
131.
J. C.
Slater
, “
The self-consistent field method for crystals
,” in
Computational Methods in Band Theory
, edited by
P. M.
Marcus
,
J. F.
Janak
, and
A. R.
Williams
(
Springer US
,
Boston, MA
,
1971
), pp.
447
457
.
132.
A. R.
Williams
,
R. A.
deGroot
, and
C. B.
Sommers
, “
Generalization of Slater’s transition state concept
,”
J. Chem. Phys.
63
,
628
631
(
1975
).
133.
L.
Gallandi
and
T.
Körzdörfer
, “
Long-range corrected DFT meets GW: Vibrationally resolved photoelectron spectra from first principles
,”
J. Chem. Theory Comput.
11
,
5391
5400
(
2015
).
134.
E. E.
Salpeter
and
H. A.
Bethe
, “
A relativistic equation for bound-state problems
,”
Phys. Rev.
84
,
1232
1242
(
1951
).
135.
C.
Liu
,
J.
Kloppenburg
,
Y.
Yao
,
X.
Ren
,
H.
Appel
,
Y.
Kanai
, and
V.
Blum
, “
All-electron ab initio Bethe-Salpeter equation approach to neutral excitations in molecules with numeric atom-centered orbitals
,”
J. Chem. Phys.
152
,
044105
(
2020
).
136.
L.
Himanen
,
A.
Geurts
,
A. S.
Foster
, and
P.
Rinke
, “
Data-driven materials science: Status, challenges, and perspectives
,”
Adv. Sci.
6
,
1900808
(
2019
).
137.
L.
Keller
and
D.
Golze
(
2020
), “
gw_core65_rc
,” Dataset in NOMAD repository. .
138.
D.
Golze
(
2020
), “
CORE65 benchmarkset
,” Dataset in NOMAD repository. .

Supplementary Material