Accurate Prediction of Band Structure of FeS2: A Hard Quest of Advanced First-Principles Approaches

The pyrite and marcasite polymorphs of FeS2 have attracted considerable interests for their potential applications in optoelectronic devices because of their appropriate electronic and optical properties. Controversies regarding their fundamental band gaps remain in both experimental and theoretical materials research of FeS2. In this work, we present a systematic theoretical investigation into the electronic band structures of the two polymorphs by using many-body perturbation theory with the GW approximation implemented in the full-potential linearized augmented plane waves (FP-LAPW) framework. By comparing the quasi-particle (QP) band structures computed with the conventional LAPW basis and the one extended by high-energy local orbitals (HLOs), denoted as LAPW + HLOs, we find that one-shot or partially self-consistent GW (G 0 W 0 and GW 0, respectively) on top of the Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation with a converged LAPW + HLOs basis is able to remedy the artifact reported in the previous GW calculations, and leads to overall good agreement with experiment for the fundamental band gaps of the two polymorphs. Density of states calculated from G 0 W 0@PBE with the converged LAPW + HLOs basis agrees well with the energy distribution curves from photo-electron spectroscopy for pyrite. We have also investigated the performances of several hybrid functionals, which were previously shown to be able to predict band gaps of many insulating systems with accuracy close or comparable to GW. It is shown that the hybrid functionals considered in general fail badly to describe the band structures of FeS2 polymorphs. This work indicates that accurate prediction of electronic band structure of FeS2 poses a stringent test on state-of-the-art first-principles approaches, and the G 0 W 0 method based on semi-local approximation performs well for this difficult system if it is practiced with well-converged numerical accuracy.

Despite progress towards understanding the origin of the low V OC in pyrite FeS 2 (Rahman et al., 2020), consensus is still not reached on the fundamental band gaps of the two FeS 2 phases. Experimentally, values varying from 0.6 to 2.6 eV have been reported for pyrite, primarily due to differences in sample preparation, measuring technique, and analytical model of spectra used in experimental studies (Ferrer et al., 1990;Ennaoui et al., 1993). Measurements of the pyrite band gap are generally carried out through optical absorption spectroscopy (Schlegel and Wachter, 1976;Kou and Seehra, 1978;Ennaoui et al., 1993), which features the neutral excitation (exciton) instead of the charged one as in the photo-electron spectroscopy (PES). Therefore, the measured excitation energies are in fact coupled to the electron-hole binding. Careful investigation by absorption spectroscopy for the marcasite phase is done only recently and gives an optical gap similar to pyrite, which essentially precludes the possibility of marcasite being the culprit for the low V OC of FeS 2 photovoltaics Wu et al., 2016). Furthermore, even though PES measurements of pyrite FeS 2 have been conducted (Ohsawa et al., 1974;van der Heide et al., 1980;Folkerts et al., 1987;Mamiya et al., 1997;Ollonqvist et al., 1997;Nesbitt et al., 2003), combined studies of direct and inverse PES (IPS) for regions near the Fermi level are rare. Reported relevant works (Folkerts et al., 1987;Mamiya et al., 1997) were done more than 20 years ago and the spectra were not resolved enough to identify a well-defined fundamental band gap.
Difficulties in characterizing band structures of FeS 2 polymorphs are also encountered from the perspective of firstprinciples calculations. Within the framework of density functional theory (DFT) (Hohenberg and Kohn, 1964), calculations with Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation (GGA) (Perdew et al., 1996a) predict pyrite to have a band gap of about 0.3 eV smaller than the experimental value of 0.95 eV as generally accepted (Ennaoui et al., 1993;Schena et al., 2013;Kolb and Kolpak, 2013;Li et al., 2015). Considering the well-known band gap problem of local density approximation (LDA) or GGA (Perdew et al., 1982), orbitaldependent functionals in spirit of generalized Kohn-Sham (GKS) DFT (Seidl et al., 1996;Perdew et al., 2017;Zhang et al., 2020) are also employed to tackle the problem, e.g. PBE plus the Hubbard-U correction (DFT + U) and hybrid functionals (Becke, 1993a;Becke, 1993b;Perdew et al., 1996b). Using an ad hoc U of 2 eV, the PBE + U method is able to reproduce the experimental band gap (Sun et al., 2011;Hu et al., 2012;Li et al., 2018) but meanwhile deteriorates the simulated optical spectra compared to PBE (Choi et al., 2012;Schena et al., 2013). Furthermore, despite the good performance in predicting band gaps for typical semiconducting materials (Heyd et al., 2005;Paier et al., 2006b,a;Marsman et al., 2008), hybrid functionals such as Heyd-Scuseria-Ernzerhof (HSE) method (Heyd et al., , 2006 have been shown to give large band gaps for pyrite of over 2 eV (Muscat et al., 2002;Sun et al., 2011;Choi et al., 2012;Hu et al., 2012;Schena et al., 2013;Liu et al., 2019). There are also works using beyond-DFT methods, particularly, the GW method based on many-body perturbation theory (MBPT) (Hedin, 1965). However, the GW results for the pyrite phase are rather scattered, ranging from 0.3 to 1.1 eV (Choi et al., 2012;Lehner et al., 2012;Kolb and Kolpak, 2013;Schena et al., 2013). It is worth noting that Schena and coworkers conducted the state-of-the-art all-electron G 0 W 0 calculations with the linearized augmented plane-wave (LAPW) basis for both pyrite and marcasite, and report a pyrite band gap only about 0.3 eV (Schena et al., 2013). The GW gap value is smaller than that from PBE, which is rarely observed in GW practices and hence deserves closer investigation.
For GW implementations involving explicit summation of states, it is established recently by a number of works (Friedrich et al., 2006;Friedrich et al., 2011a;Friedrich et al., 2011b;Klimes et al., 2014;Jiang and Blaha, 2016;Nabok et al., 2016;Jiang, 2018;Zhang and Jiang, 2019;Ren et al., 2021) that an accurate description of high-lying empty states is essential to give accurate correlation self-energy operator and consequent QP band structure. In the pseudo-potential framework, one can improve the accuracy by using a norm-conserving potential with specifically tailored projectors at high energies (Klimes et al., 2014;van Setten et al., 2018). In all-electron calculations with the LAPW basis set, local orbitals with large energy parameters (usually 10 1∼2 Ry higher than the Fermi level) are introduced as additional basis functions to remove the linearization error in unoccupied states up high in the conduction band regime (Friedrich et al., 2006(Friedrich et al., , 2011aJiang and Blaha, 2016;Nabok et al., 2016). The LAPW basis extended by these high-energy local orbitals (HLOs), termed as LAPW + HLOs, has succeeded in helping produce accurate QP band structures in good agreement with experiment for a variety of semiconductors (Jiang and Blaha, 2016) including the conventionally challenging systems such as ZnO (Friedrich et al., 2011a;Friedrich et al., 2011b;Stankovski et al., 2011;Jiang and Blaha, 2016;Nabok et al., 2016), d/f-electron oxides (Jiang, 2018) and cuprous and silver halides (Zhang and Jiang, 2019). Particularly, the effects of including HLOs on the QP correction have been demonstrated quantitatively to be larger for states with stronger metal-d characters (Zhang and Jiang, 2019). For the FeS 2 polymorphs with states of significant Fe-3d characters in both valence and low-energy conduction band regimes, GW with LAPW + HLOs is likely to give better description of the QP energies and dispersion relation than that with the standard LAPW basis.
A competitive alternative in the DFT framework to GW for band structure prediction is the doubly screened hybrid (DSH) functional method  in the category of hybrid functionals with system-dependent parameters . Derived from a model dielectric function (Cappellini et al., 1993;Shimazaki and Yoshihiro, 2008), the exchangecorrelation potential in DSH can be regarded as a further approximation to the Coulomb hole and screened exchange (COHSEX) approximation to the GW self-energy, and is able to capture both dielectric and metallic screening in the exchange interaction . It is shown that the DSH can evaluate band gaps of typical sp semiconductors with accuracy comparable to GW with the LAPW + HLOs basis while only at modest computational cost . Furthermore, the one-shot variant DSH0 can outperform fixed-parameter hybrid functionals for band gap predictions in a wide range of materials including narrow-gap semiconductors and transition metal mono-oxides Liu et al., 2020). Hence we consider DSH as a hopeful approach to solve the FeS 2 band gap puzzle within the GKS framework of DFT.
In the present work, we investigate the electronic band structures of the pyrite and marcasite polymorphs of FeS 2 by applying the state-of-the-art all-electron GW method with the LAPW + HLOs basis. For comparison, we examine the results from GW with the standard LAPW basis as well. We also investigate the performances of several hybrid functionals, including PBE0 (Perdew et al., 1996b), HSE06 , Heyd et al., 2006, screened-exchange-PBE hybrid functional (SX-PBE) (Bylander and Kleinman, 1990;Seidl et al., 1996) and DSH , in attempt to obtain insights into the failure of the conventional fixedparameter functionals in predicting the band gap of FeS 2 .

The GW Method
The central task of the GW method is to solve the quasi-particle (QP) equation with the self-energy operator Σ in the frequency domain expressed as (Hedin, 1965) where G is the time-ordered Green's function With ψ nk and ε nk being the wave function and energy of the single-particle state |nk〉 respectively, μ the chemical potential, and δ and η positive infinitesimals. Atomic units are used throughout the paper. The screened Coulomb interaction W writes where v (r, r′) 1/|r − r′| is the bare Coulomb interaction and ε(r, r′; ω) is the microscopic dielectric function calculated at the level of random phase approximation (RPA). In principle, Eqs 1-3 have to be solved self-consistently along with the Dyson equation for the Green's function (Hedin, 1965). However, due to the computational cost and generally unsatisfactory results of the fully self-consistent GW for solids [e.g. Grumet et al. (2018)], one usually turns to the non-self-consistent variant G 0 W 0 . Considering the resemblance of KS and QP wave functions in weakly correlated systems (Hybertsen and Louie, 1986), the selfenergy or QP energy ε QP nk can be computed perturbatively upon the acquisition of Σ from the KS states as where V xc is the KS exchange-correlation potential and Z nk a renormalization factor. One can further perform the socalled energy-only self-consistent GW 0 calculations, where QP energies ε QP nk in place of ε KS nk in Eq. 2 are updated iteratively while W is kept the same as in G 0 W 0 (Shishkin and Kresse, 2007). The GW method has been implemented in various numerical frameworks (Jiang, 2011;Golze et al., 2019). For a detailed explanation of the basic theory and computational techniques used in the present GW implementation, the readers can refer to Jiang et al. (2013).

All-Electron Calculations With HLOs-Extended LAPW Basis
In the all-electron framework with LAPW, KS wave functions are expanded by the LAPW basis (Andersen, 1975;Singh and Nordström, 2006;Blaha et al., 2020) where V α is the region enclosed by the muffin-tin (MT) sphere of atom α centered at r α with radius R α MT , r α r − r α , u αl (E αl ) is the solution of radial KS equation inside V α at chosen energy E αl , _ u αl (E αl ) ≡ zu αl (E)/zE| E E αl , and Y m l is the spherical harmonic function. The coefficients A k+G αlm and B k+G αlm are determined by enforcing that φ LAPW k+G (r) be smooth at the boundary of V α . Local orbitals (LOs) which vanish outside the atomic spheres are proposed to supplement the LAPW basis to better describe the semi-core states (Singh, 1991). Inside the atomic sphere V α , LOs take the following form where E LO,i αl is the energy parameter for the ith LO centered on atom α with angular and azimuthal quantum numbers l and m, respectively.
HLOs fall into the category of LOs with E LO αl typically 10 ∼ 100 Ry above the Fermi level. Such extra LOs have been found to facilitate accurate description of unoccupied states by remedying the linearization error therein when using the LAPW basis (Krasovskii et al., 1994;Krasovskii, 1997;Friedrich et al., 2006;Michalicek et al., 2013). In ground state calculations with LDA/ GGA or hybrid functionals, the error causes no essential difficulties, since only occupied and low-lying unoccupied states are involved which are usually handled in sufficient accuracy with the usual or standard LAPW basis generated as default in popular DFT implementations with LAPW basis (Blaha et al., 2020). However, the error can be detrimental to the numerical accuracy of methods where the summation over unoccupied states is required, e.g. GW and DFT methods with density approximations belonging to the fifth rung of Jacobi ladder (Perdew and Schmidt, 2001) such as the adiabaticconnection dissipation-fluctuation (ACFD) calculation under RPA for ground-state energy (Ren et al., 2012;Cui et al., 2016;Zhang et al., 2018). In these methods, the completeness of summation and quality of unoccupied states play a crucial role. Previous GW studies (Jiang and Blaha, 2016;Jiang, 2018;Zhang and Jiang, 2019;Shen et al., 2020) have suggested that both can be taken into account by including localized orbitals energetically higher than the Fermi level in addition to the standard LAPW basis. HLOs have been shown to effectively improve the optical properties (Krasovskii et al., 1994;Krasovskii, 1997), NMR chemical shifts (Laskowski and Blaha, 2012;Laskowski and Blaha, 2014), GW QP energies (Friedrich et al., 2006;Friedrich et al., 2011a;Jiang and Blaha, 2016;Nabok et al., 2016;Jiang, 2018;Zhang and Jiang, 2019), optimized effective potential (Betzinger et al., 2011(Betzinger et al., , 2012 and RPA correlation energy (Betzinger et al., 2015).
In the current implementation, HLOs are generated systematically by following the way described by Laskowski and Blaha (2012). The quality of LAPW + HLOs is controlled by two parameters besides those for the LAPW basis, namely, the additional number of nodes in the radial function of highest energy local orbital with respect to that of the LAPW function with the same angular quantum number and the maximum angular quantum number of used HLOs, denoted as n LO and l (LO) max , respectively. Generally speaking, the larger n LO and l (LO) max are, the higher the HLOs can reach in the energy space. We use n LO 0 to denote the usual or standard LAPW basis. Since the convergence rate of the QP energy with respect to the two parameters can be different for states featuring distinct atomic characters, careful convergence check is required to obtain numerically accurate GW results.

Hybrid Functionals
Hybrid functionals have been widely used in first-principles simulations of condensed matter for their good balance between performance and computational cost, and have been actively developed to further exploit the potential of its particular functional form. Readers interested in detailed description on the current status of hybrid functional development are directed to several recent reviews (Kümmel and Kronik, 2008;Baer et al., 2010;Maier et al., 2019;Zhang et al., 2020). Here we briefly introduce the general formalism of the range-separated hybrid functionals and the variants relevant to the current study.
The essential ingredient in hybrid functional methods is the exchange-correlation energy E xc or potential V xc composed of non-local orbital-dependent (screened) Hartree-Fock (HF) exchange terms. In the present work, we focus on hybrid functionals with V xc in the range-separated form as where V HF,sr x and V HF,lr x are the short-and long-ranged Fock exchange potentials, respectively, which, using the reduced density-matrix defined as ρ( In Eq. 8 v sr (r, r′; μ) denotes the short-ranged Coulomb interaction of a certain form characterized by screening parameter μ . x denotes collectively the spatial and spin coordinates of an electron, x ≡ (r, σ). V SL,sr x and V SL,lr x are the semi-local (SL) counterparts of the exchange potentials in LDA, GGA or meta-GGA. μ and the mixing ratios α sr and α lr are the adjustable parameters of the hybrid functional form.
Conventionally, the parameters are determined by either theoretical analysis or fitting against some dataset of particular properties, and then applied to other systems as fixed. Famous examples of the fixed-parameter hybrid functionals include PBE0 α sr α lr 1/4 (Perdew et al., 1996b) and the HSE series α sr 1/4, α lr 0, μ 0.2-0.3 Å -1 Heyd et al., 2006). Recently, hybrid functionals with system-dependent parameters are developed by several groups (Shimazaki and Yoshihiro, 2008;Marques et al., 2011;Kronik et al., 2012;Koller et al., 2013;Skone et al., 2014;Chen et al., 2018;Cui et al., 2018). Among different methods, the doubly screened hybrid (DSH) functional has been demonstrated as a competitive candidate for accurate description of band structures of both wide-and narrow-gap semiconductors . The underlying idea of DSH is to approximate the screening effect in solids by employing the Bechstedt model dielectric function (Bechstedt et al., 1992) where ε M is the macroscopic dielectric constant, q TF the Thomas-Fermi wave vector and α an empirical parameter chosen for semiconductors (Cappellini et al., 1993). A screened Coulomb interaction can be derived from this model to take both dielectric and metallic screening into account, leading to parameters in Eq. 7 as Frontiers in Chemistry | www.frontiersin.org September 2021 | Volume 9 | Article 747972 α sr 1, The corresponding short-ranged Coulomb interaction in Eq.
where erfc is the complementary error function. In practice, an initial ε M is required, which can be obtained from the PBE calculation or experimental measurements, to construct the DSH potential and solve the GKS equation. The resulting single-particle states act as the inputs to compute a new ε M , which is in turn used to update the DSH potential. The self-consistent loop stops when ε M is converged. Alternatively, one can break after solving the GKS equation with the initial ε M , leading to the one-shot scheme denoted as DSH0.

Computational Details
The unit cells of pyrite and marcasite FeS 2 used in our calculations are shown in the left panel of Figure 1. The crystal structure of pyrite FeS 2 ( Figure 1A) can be viewed as a faced-centered cubic cell of Fe atoms with S 2 dumbbells occupying the octahedral interstitials and pointing to different <111> crystallographic axes. The anion coordination octahedra (FeS6) are connected only through sharing vertices. In the orthorhombic marcasite phase ( Figure 1C), (FeS6) are connected by sharing edges with the two neighbors along c-axis and linked together through sharing vertices on the aOb plane. In terms of lattice parameters, we use a 5.418 Å, u 0.3850 for pyrite (space group Pa3) and a 4.443 Å, b 5.425 Å, c 3.387 Å, u 0.2005, v 0.3783 for marcasite (space group Pnnm). These values follow the results from X-ray diffraction experiments at ambient conditions (Brostigen and Kjekshus, 1969;Brostigen et al., 1973;Chattopadhyay and Von Schnering, 1985;Zuñiga-Puelles et al., 2019). The corresponding S-S bond lengths in the two polymorphs are 2.16 and 2.21 Å, respectively.
The present all-electron GW calculations are performed by the GW facilities in the GAP2 program (Jiang et al., 2013;Jiang and Blaha, 2016) interfaced to WIEN2K (Blaha et al., , 2020. Results in both G 0 W 0 and GW 0 schemes are presented, where KS orbital energies and wave functions calculated with the PBE (Perdew et al., 1996a) GGA are used as the input to construct one-body Green's function and screened Coulomb interaction. The KS states are obtained by using charge density preconverged under self-consistent field (SCF) calculation with PBE and the standard LAPW basis. The energy criterion for convergence of SCF iterations is set to 10 -8 Rydberg (Ry). 64 (4 × 4 × 4) and 120 (5 × 4 × 6) k points are sampled in the first Brillouin zones of pyrite and marcasite FeS 2 , respectively. All available unoccupied states are considered in the summation of states for screened Coulomb interaction W and self-energy Σ. Mixed product basis is used to describe the wave function products in the two-point functions, e.g. W and Σ (Aryasetiawan and Gunnarsson, 1994;Kotani and van Schilfgaarde, 2002). We choose Q 0.75 and l MB max 3for the interstitial plane wave and MT product basis, respectively [Jiang et al. (2013) for the meanings of these parameters]. LAPW and LOs with E αl < 20 Ry are used to build the MT product basis. Frequency dependence of W is treated explicitly on a 16-point double Gauss-Legendre grids along the positive imaginary axis. Σ on the same grid is calculated and analytically continued to the real axis (Rojas et al., 1995). A rather coarse k/q-point mesh, 2 × 2 × 2 for pyrite and 4 × 2 × 4 for marcasite, is sufficient to converge the direct band gap at the Γ point E g Γ within 0.01 eV. The QP band structure diagrams along particular k-point paths (see the right panel of Figure 1) are calculated by interpolating the QP energies obtained with the above mesh using the Fourier interpolation technique (Pickett et al., 1988).
In terms of the LAPW basis, the usual or standard LAPW basis set is created automatically in the recent version of WIEN2K , which is actually a mixture of the APW + lo basis for the valence states , the ordinary LAPW basis for higher l channels up to l max 10 and additional local orbitals (LOs) for semi-core Fe-3s and Fe-3p states . The convergence with respect to the two HLOs parameters n LO and l (LO) max is investigated, the latter being represented by max is the largest l of valence orbitals for each element. In the present study, l (v) max 2 and l (v) max 1 for Fe and S, respectively. Since the convergence with respect to HLOs parameters are decoupled from the choice of k-point mesh, we choose marcasite with a coarse 2 × 1 × 2 mesh for HLOs convergence test. RKmax ≡ R MT,min K max 7.0 is chosen for the plane-wave cut-off in the interstitial region, where R MT,min is the minimal muffin-tin radius R MT used in the lattice. In the present FeS 2 case, R MT is set to 2.1 Bohr for Fe and 1.9 Bohr for S. Using RKmax 9.0 will reduce the band gap from GW (LAPW + HLOs) by less than 0.03 eV, indicating that adequate accuracy can be delivered with the current RKmax 7.0 setup. Due to limited computational resources, RKmax 6.0 is used for HLOs convergence test. Following Laskowski and Blaha (2012), the linear independence of HLO basis functions is assured by choosing the energy parameters such that the overlap between the HLO radial functions is smaller than a threshold, which is 0.6 in the present work.
For hybrid functional calculations, we consider PBE0 (Perdew et al., 1996b), HSE06 (Heyd et al., , 2006 and screened exchange SX-PBE (Bylander and Kleinman, 1990) methods as well as DSH. All hybrid functional calculations are performed with the projector augmented waves (PAW) method (Blöchl, 1994) implemented in the Vienna ab-initio Simulation Package (VASP) (Kresse and Furthmüller, 1996). The static dielectric function is calculated from the average of diagonal elements of macroscopic dielectric tensor computed by using density functional perturbation theory (DFPT) with local field effect included (Baroni et al., 2001). Apart from 3d and 4s, the 3s and 3p electrons of Fe are also treated explicitly in the valence region. The Thomas-Fermi wave vectors are 2.57 and 2.56 Å -1 for pyrite and marcasite FeS 2 , respectively. The cut-off energy of plane-wave basis for wave function expansion is chosen as 400 eV, which is sufficient to converge E g Γ of both FeS 2 polymorphs within 2 meV. In terms of k-point mesh, 64 (4 × 4 × 4) and 120 (5 × 4 × 6) k points are sampled in the first Brillouin zones of pyrite and marcasite for the self-consistent calculations, respectively. Using a finer 6 × 6 × 6 sampling for pyrite will change the band gap by less than 0.01 eV, and hence we consider the results well converged with respect to the k-point mesh. The energy convergence criterion is chosen to be 10 -6 eV for the SCF iterations.

3.1The GW Results
In this part, we present the electronic band structures of pyrite and marcasite FeS 2 computed by the all-electron GW method. In particular, we analyse the effect of high-energy local-orbitals (HLOs) by comparing the results from GW with the standard LAPW and LAPW + HLOs basis.

Convergence of QP Energies with Respect to HLOs Parameters
To achieve a balance between the computational cost and numerical accuracy of the LAPW + HLOs based GW method, we have to decide an optimized HLOs setup for the FeS 2 polymorphs of interest. That is to say, certain convergence with respect to the two HLOs parameters, namely n LO and Δl LO , must be achieved for the QP band structures of both polymorphs, while the number of basis functions should be kept as few as possible. To simplify the notation, we denote the setup of HLOs by (n LO , Δl LO ) so that (1, 1) indicates a set of HLOs with n LO 1 and Δl LO 1, for example. Since we are most interested in the band gaps (direct and indirect) of the systems, we choose the indirect band gap from the X point to the Γ point, E X-Γ g , as the descriptor for the band structure, and investigate its dependence on the two HLOs parameters for marcasite.
Before discussing the results, we briefly illustrate the appropriateness of this choice. First of all, E X-Γ g is a representative band gap energy for pyrite and marcasite FeS 2 . This is because in both phases, the topmost valence state at the X point, X v , is close to the valence band maximum (VBM) and the bottommost conduction state at the Γ point, Γ c , is the conduction band minimum (CBM) (that is the case for marcasite given the coarse 2 × 1 × 2 k mesh in the convergence study). Second, either X v or Γ c has similar atomic contributions in the two polymorphs, and the effects of HLOs on such states are also similar, as shown in the results for other polymorphs like zinc-blende and wurtzite ZnO (Jiang and Blaha, 2016). Therefore the parameters optimized for marcasite are considered transferable and can be applied to the pyrite polymorph. Last but not least, as we will discuss later, the effects of HLOs on X v and Γ c differ significantly, avoiding considerable error cancellation in change of the QP correction to the band gap upon including HLOs.   Figure 2 summarizes the results of convergence study for the G 0 W 0 @PBE method. E X-Γ g ( Figure 2A) is about 1.0 eV with the standard LAPW basis (n LO 0) and is significantly increased by extending LAPW with HLOs. One can see that the convergence rate of E X-Γ g with respect to n LO differs with different Δl LO , and is faster for lower Δl LO . The reverse is also true, i.e. the convergence with respect to Δl LO is faster when n LO is smaller. It clearly indicates that the convergence with respect to n LO and Δl LO is coupled. Increasing HLOs parameters from (4, 4) to (5, 5) changes E X-Γ g by less than 0.05 eV, indicating that HLOs (4,4) is able to deliver an adequate accuracy. Therefore, unless stated otherwise, HLOs (4, 4), amounting to 196 and 145 HLOs for Fe and S atoms, respectively, is considered optimized and will be used in the subsequent GW calculations denoted by LAPW + HLOs. The energy parameters for HLOs (4, 4) can be found in Table 1.
It is worth noting that the effect of including HLOs on the QP correction to E X-Γ g is different from those on the valence and conduction states. To illustrate this, we show the dependence on n LO and Δl LO of the self-energy corrections to Γ c (Δε Γ c ) and X v (Δε X v ) states in Figures 2B,C, respectively. Both Δε Γ c and Δε X v decrease with increasing n LO or Δl LO , but the former converges much faster than the latter, which agrees with the general trend observed previously (Jiang and Blaha, 2016;Zhang and Jiang, 2019). With the standard LAPW basis, G 0 W 0 @PBE gives Δε Γ c 0.48 eV and Δε X v 0.76 eV, indicating a negative QP correction to the band gap, which is rarely observed in LDA/GGA-based GW calculations of semiconductors (Jiang and Blaha, 2016). When HLOs (5, 5) are included, Δε Γ c decreases by 0.3 eV, much smaller compared to the decreasing of 1.2 eV in Δε X v . Such biased effects of including HLOs on valence and conduction band states can be attributed to the difference in atomic characteristics between the states, and will be further discussed in the following sections.

Quasi-Particle Band Gaps
After having obtained the optimized HLOs, we perform the PBEbased GW calculations for pyrite and marcasite FeS 2 with the LAPW + HLOs basis set, and compare with the PBE method and GW with the standard LAPW basis.
The band gaps of pyrite and marcasite FeS 2 calculated by PBE and GW methods are presented in Table 2. The fundamental band gaps are obtained by computing the band energies along the k-point paths indicated in the right panel of Figure 1. In the PBE reference, pyrite and marcasite are predicted to have indirect fundamental band gaps of 0.70 and 0.83 eV, respectively. Our PBE results are consistent with those from previous all-electron LAPW study (Schena et al., 2013) and close to the recently reported optical band gaps obtained from diffuse reflectance spectroscopy . However, our PBE band gap for pyrite is slightly larger than several reported PBE results (Sun et al., 2011;Kolb and Kolpak, 2013;Lazić et al., 2013;Li et al., 2015;Zhang et al., 2018). This can be attributed to the use of different lattice structures in those studies (Eyert et al., 1998;Lazić et al., 2013;Schena et al., 2013) from the current work. Particularly, geometry optimization by PBE (Eyert et al., 1998;Schena et al., 2013) generally gives a longer S-S dimer, which leads to smaller splitting between bonding and anti-bonding S-3pσ orbitals and a consequent shrink in the band gap.
For GW calculations with the standard LAPW basis, the QP fundamental band gaps by G 0 W 0 @PBE are smaller than the PBE counterparts in both FeS 2 polymorphs. Pyrite FeS 2 is predicted to have a band gap of only 0.06 eV, which is 0.64 eV smaller than that by PBE. The negative QP correction for pyrite band gap has been reported by Schena et al. (2013). The QP fundamental band gap for marcasite predicted by G 0 W 0 @PBE (LAPW) is also smaller than PBE, while the change (0.26 eV) is less dramatic 2 | Fundamental band gap (indicated by "fund.") and other direct and indirect band gaps (unit: eV) for pyrite and marcasite FeS 2 calculated by PBE and GW methods. Results from previous GW studies and experimental measurements are presented for comparison. To simplify the notation, we use "L" and "L + H" to denote the standard LAPW and LAPW + HLOs basis sets, respectively. PBE is used as the starting point for G 0 W 0 and GW 0 calculations unless stated otherwise.

Pyrite Marcasite
Fund.  than that for pyrite. Such negative QP corrections to LDA/GGA band gaps are uncommon in GW studies for closed-shell systems (Klimes et al., 2014;Jiang and Blaha, 2016;van Setten et al., 2017;Zhang and Jiang, 2019) as well as open-shell d/f-electron semiconductors (Jiang, 2018). Switching on self-consistency of the Green's function by GW 0 @PBE further reduces the fundamental band gaps of FeS 2 . In particular, pyrite is predicted to be metallic by GW 0 @PBE, which disagrees qualitatively with its semiconducting nature in experiment (Ennaoui et al., 1993). For other direct and indirect band gaps, those for Γ → Γ and X → Γ in pyrite and marcasite and M → Γ in marcasite are decreased from PBE to G 0 W 0 @PBE (LAPW). The decrease is largest for the Γ → Γ gap in the two phases, 0.71 and 0.86 eV for pyrite and marcasite, respectively. On the other hand, the gaps for X → X in pyrite, Γ → T and X → T in marcasite are increased by 0.28, 0.25, and 0.58 eV, respectively. However, it should be noted that the distinction in signs of corrections to the QP gaps in different channels should not be considered as intrinsic for FeS 2 . Instead, it is an artifact as a result of the incomplete basis, which we will discuss in details below. Now we turn to the LAPW + HLOs-based GW calculations. With the G 0 W 0 @PBE method, including HLOs increases the QP fundamental gap by 0.98 eV for pyrite and 0.58 eV for marcasite. The resulting G 0 W 0 @PBE band gaps are 1.04 and 1.15 eV for pyrite and marcasite, respectively. In contrast, all band gaps investigated are increased by G 0 W 0 with LAPW + HLOs compared to their PBE counterparts. We note that HLOs have distinct effects among band gaps for different channels. Once the HLOs are included, band gaps for channels with the conduction state at the Γ point are increased by about 1 eV. On the other hand, the QP correction to the X → X band gap in pyrite increases by only 0.18 eV. Moreover, the gaps for Γ → T and X → T in marcasite even decrease. With the LAPW + HLOs basis, using GW 0 to switch on partial self-consistency further increases the band gaps, but the change is moderate and no more than 0.1 eV.
As explained at the beginning, the fundamental band gap of FeS 2 has been controversial in the recent decades, partly due to the widely varying experimental values (Ennaoui et al., 1993). In the present study, the GW 0 @PBE method with the LAPW + HLOs predicts that pyrite and marcasite have indirect fundamental band gaps of 1.14 and 1.16 eV, respectively. The GW 0 gap of pyrite is slightly larger than the generally accepted experimental value of 0.95 eV (Ennaoui et al., 1993). Furthermore, the fact that the two polymorphs have almost identical band gaps is consistent with the optical measurements by Sánchez et al. (2016), although our predicted band gaps are about 0.3 eV larger. However, it should be noted that one must take exciton binding energy E B into account for a meaningful comparison between the QP fundamental band gap and experimentally measured optical gap. The difference between the fundamental and optical gaps can be significant when the exciton is localized, i.e. of Frenkel type (Fox, 2010). On the other hand, while it is more straightforward to compare the QP gap with spectral data from direct and inverse PES (Folkerts et al., 1987;Mamiya et al., 1997), the resolutions of available measurements for pyrite FeS 2 are too low to extract a meaningful gap value for comparison. Moreover, to the best knowledge of the authors, no data of combined PES/IPS measurements are available for marcasite. Therefore, further experimental studies are required to determine and verify the band gaps of the FeS 2 polymorphs.
To close this part, we highlight that the present work resolve two issues reported in previous GW studies in terms of QP band structures of FeS 2 . First, Schena et al. (2013) performed a G 0 W 0 @ PBE study on pyrite and marcasite FeS 2 with similar HLOsextended LAPW basis. The fundamental band gap of pyrite was estimated as about 0.3 eV, by which the authors claimed to explain the low V OC encountered in the pyrite solar cell. However, according to our convergence study, such a small band gap is likely to result from inadequate convergence with respect to HLOs. More specifically, the largest angular momentum of HLOs l (LO) max used in Schena et al. (2013) is 3, i.e. f orbital, while l (LO) max 2 + 4 6 (i orbital) is used in the optimized HLOs of the present work. As a result, the highest energy covered by HLOs in Schena et al. (2013) (800 eV) is much smaller than that used in the present work (about 1800 eV). Second, fully self-consistent GW (scGW) and quasi-particle selfconsistent GW (QPscGW or QSGW) calculations have also been carried out to study the band structure of pyrite, and give apparently satisfactory results (Lehner et al., 2012;Kolb and Kolpak, 2013). However, variants of self-consistent GW without taking the vertex function into account tend to overestimate the band gaps of typical semiconductors, as indicated by several works (Shishkin and Kresse, 2007;Deguchi et al., 2016;Cao et al., 2017;Grumet et al., 2018). Thus the error cancellation between the general tendency of overestimating band gaps of semiconductors and the numerical inaccuracy in the LAPW basis or the use of conventional pseudo-potentials could contribute to the apparent agreement between the generally accepted band gap and the self-consistent GW results. Of course, without looking into computational details of previous self-consistent GW calculations, this is just our speculation. Further investigations are needed to fully clarify this issue. We also note that similar LAPW + HLOs calculation has been conducted for the pyrite phase by Ouarab and Boumaour (2017) and gives a band gap (0.97 eV) close to ours.

Quasi-Particle Band Structure
To further illustrate the significance of HLOs in applying the GW methods to FeS 2 , we present the QP band structures of pyrite and marcasite FeS 2 calculated from the G 0 W 0 @PBE method with the LAPW + HLOs basis set, and compare the results to those with the standard LAPW basis. Figure 3 shows the electronic band structures of the two FeS 2 phases from different methods. Note that the bands are aligned to the CBM at the Γ point for a better view of QP correction to the valence states. With PBE, pyrite ( Figure 3A) is found to be an indirect band gap material with the CBM located at the Γ point and the VBM near the X point along Γ-X. The top valence bands within 1 eV below the VBM are dominated by the localized Fe-3d states, also manifested by their flat dispersion. The dispersive bands about 2 eV below the VBM are mainly composed of S-3p states and well separated from the Fe-3d (t 2g ) valence bands. In the conduction band region, the lowermost conduction bands are also largely composed of Fe-3d (e g ), except for the states close to the Γ point with predominant S-3p characters. Particularly, the CBM Γ c state is exclusively formed by the σ anti-bonding overlapping of S-3p orbitals in the S-S dimer (see projected bands in Figure 4A). Valence and conduction bands with strong Fe-3d characters are separated by about 2 eV. For marcasite, an indirect band gap is also observed, with the CBM located at the T point (T c ) and the VBM along Γ-X (Δ v ). Both states at the VBM and CBM of marcasite are of dominant Fe-3d characters ( Figure 4D), in contrast to pyrite where CBM is of pure S-3p characters. The wider Fe-3d valence bands overlap with the S-3p bands near about 1.5 eV below the VBM, which indicates stronger covalent bonding between Fe and S in marcasite than in pyrite.
Then we compare the QP band structures obtained from G 0 W 0 @PBE with the LAPW and LAPW + HLOs basis (Figure 3). With the standard LAPW basis, G 0 W 0 @PBE predicts pyrite almost as a semimetal with a nearly vanishing band gap ( Figure 3A). Dispersion of the conduction band around the Γ point and the separation between the Fe-3d and S-3p valence bands are enhanced compared to the PBE reference. For marcasite ( Figure 3B), although a noticeable gap (0.57 eV) is predicted by G 0 W 0 @PBE, the band edges are different from those in PBE: the CBM is located at the Γ point (Γ c ) and the VBM in the middle of the Z-Γ path (Λ v ). The change in the nature of band edges from semilocal functional to GW method is also observed by Schena et al. (2013). Once HLOs are included in the basis set, QP band gaps of both phases are dramatically enlarged. The fundamental gaps of pyrite and marcasite are 1.04 and 1.15 eV, respectively, which are 0.2 ∼ 0.3 eV larger than the optical gaps from absorption spectra . Band edges of marcasite by GW are also recovered to those by PBE. The comparison indicates that both negative QP corrections to band gaps and change of band edges in GW (LAPW) are indeed artifacts due to the inadequate numerical accuracy of the basis set.
To better understand how the HLOs basis functions influence the QP band structures of FeS 2 , we scrutinize the QP correction to Kohn-Sham state Δε, defined by the difference between the QP energy ε QP and the KS energy ε KS , i.e. Δε ≡ε QP − ε KS . For pyrite, with the standard LAPW basis, Δε to the CBM is smaller than those to the valence Fe-3d t 2g and conduction Fe-3d e g states as shown in Figure 4B. Particularly, Δε for the VBM is about 0.7 eV greater than that for the CBM. This leads to a up-shift of Fe-3d states with respect to the CBM on a whole. Extending LAPW with HLOs reduces Δε for all states, but the reduction in Δε to the VBM is more than that to the CBM by about 1.0 eV, resulting in the sign change of the QP correction to the band gap. Similar conclusion can be drawn from Δε in the marcasite phase ( Figure 4E). With the standard LAPW method, Δε to T c exceeds that to Γ c by more than 1.1 eV. Consequently, Γ c drops down below T c and becomes the CBM, as we have seen in Figure 3B. Upon including HLOs, Δε to T c is reduced more significantly than Δε to Γ c such that T c recovers the conduction band edge as in PBE.
Such biased effects of HLOs are clearly associated with the atomic characteristics of Kohn-Sham states, as we have demonstrated in the GW calculations of cuprous and silver halides (Zhang and Jiang, 2019). In Figures 4C,F, we plot the difference between Δε computed by G 0 W 0 with LAPW + HLOs and LAPW against the weight of Fe-d characters of the Kohn-Sham orbitals |ψ nk 〉, w Fe−d nk , defined by where ϕ Fe i l 2,m represents the pre-defined atomic function centered on the ith Fe atom featuring spherical harmonic function Y m 2 . The negative difference implies that including HLOs generally brings down Δε. Moreover, the difference is more dramatic for states with larger w Fe−d nk , indicating that numerical error is more significant for states with stronger Fe-d characters in GW calculations with the incomplete LAPW basis.

GW Density of States
To end this section, we present the GW calculated density of states (DOS) of FeS 2 polymorphs in Figure 5. The results for pyrite FeS 2 are shown in Figure 5A. Due to different definitions of the Fermi level in theoretical results and experimental spectral data, we have shifted the experimental data to match up the highest valence peak near the Fermi level. With this alignment, the overall DOS from G 0 W 0 (LAPW + HLOs) agrees well with the energy distribution curves (EDCs) from the PES experiments. The width of the valence Fe-3d band and separation between the Fe-3d and S-3p valence bands are consistent with the UPS experiment by Mamiya et al. (1997) and the XPS experiment by Folkerts et al. (1987). The location of the first peak in the conduction band region is also in good agreement with the BIS data (Folkerts et al., 1987). Interestingly, although G 0 W 0 (LAPW) underestimates the band gap severely, the location of the first peak in the conduction region is almost identical to that by G 0 W 0 (LAPW + HLOs), probably due to the error cancellation between QP corrections to the valence and conduction Fe-3d bands. However, such fortuitous cancellation does not hold in the valence region as inferred by the too deep S-3p band in the G 0 W 0 (LAPW) results. Figure 5B shows the calculated DOS for marcasite. Regardless of the theoretical method used, the valence Fe-3d band of marcasite has larger width than that of pyrite, indicating a stronger Fe-S interaction in the marcasite phase. In the conduction region, a sharp peak is observed with the G 0 W 0 (LAPW) method, while only a plateau is found with G 0 W 0 (LAPW + HLOs). However, the sharp peak is actually an artifact of wrongly pushed up Fe-3d conduction bands due to the inaccuracy of the standard LAPW basis as explained above. Energy distribution curves (EDCs) of pyrite extracted from photoelectron spectroscopy (PES) are presented for comparison. Each dataset is normalized with respect to its highest peak. Theoretical data are aligned to its valence band maximum as energy zero. The XPS + BIS and UPS + BIS data for pyrite are obtained from Folkerts et al. (1987) and Mamiya et al. (1997), respectively. To take into account the different definitions of the Fermi level in theory and experiment, a rigid shift of 0.60 and 0.40 eV are employed for the EDCs from XPS + BIS and UPS + BIS, respectively, to match the highest valence peaks below the Fermi level.

Results From Hybrid Functionals
As mentioned in the introduction, previous theoretical studies found that various hybrid functionals, which are typically able to describe the band gaps of semiconductors quite accurately, performed badly for FeS 2 . In this section, we look into this issue and present results by several hybrid schemes including the DSH functional with systemdependent parameters.

Band Gaps by Hybrid Functionals
Band gaps computed by different hybrid functionals are collected in Table 3. The widely used PBE0 and HSE06 functionals have been reported to predict fundamental gaps of pyrite and marcasite FeS 2 larger than 2 eV in the literature (Sun et al., 2011;Choi et al., 2012;Hu et al., 2012;Schena et al., 2013;Liu et al., 2019), which is confirmed by our results. DSH, the hybrid functional with system-tuned parameters, does not improve the prediction over PBE0 and HSE06. This is surprising, given that DSH has been previously shown to outperform several other hybrids in evaluating band structures for wide-and narrow-gap systems Liu et al., 2020), including PBE0, HSE06 and the dielectric-dependent hybrid (DDH) functionals (Marques et al., 2011;Skone et al., 2014). SX-PBE screened exchange functional gives band gaps of FeS 2 significantly smaller than the hybrids mentioned above, but the gaps are still larger than those from GW with the LAPW + HLOs basis ( Table 2) by about 0.5 eV.
Considering that the one-shot DSH, i.e. DSH0, may outperform the self-consistent scheme in some transition metal compounds Liu et al., 2020), we also employ DSH0 to calculate the two FeS 2 polymorphs. The macroscopic dielectric constant calculated with PBE ε M PBE is 20.6 for pyrite, which agrees well with ε M PBE 21 from a previous study (Choi et al., 2012). DSH0 with ε M PBE predicts smaller band gaps than DSH, but the values are still above 2 eV. Meanwhile, DSH0 with experimentally obtained ε M 10.9 (Husk and Seehra, 1978) gives the pyrite band gap of 2.72 eV. In contrast, a modified HSE functional (MHSE) with HSE06 screening parameter and 10% hybrid ratio, which is roughly equal to the inverse of the experimental dielectric constant, as suggested by Liu et al. (2019), gives band gaps close to the GW 0 (LAPW + HLOs) result. The MHSE results agree with those by Liu et al. (2019) and seem to verify the suggestion by Schena et al. (2013) of using 1/ε M as the hybrid ratio in the HSE-type screened hybrid functional.

Band Structures by Hybrid Functionals
As summarized above, the investigated hybrid functionals except for MHSE fail to give reasonable predictions for the band gaps of pyrite and marcasite FeS 2 . In this section, we take a close look at the band structures computed from these methods to understand the failure.
The band structures for pyrite calculated from selected hybrid functionals are shown in the upper panel of Figure 6. With PBE0 and HSE06 ( Figures 6A,B), the fundamental band gap is a direct one with both VBM and CBM located at the Γ point. An indirect TABLE 3 | Fundamental band gap (indicated by "fund.") and other direct and indirect band gaps (unit: eV) for pyrite and marcasite FeS 2 calculated by different hybrid functionals. Results from other theoretical studies and experimental measurements are presented as comparison.

Methods
Pyrite Marcasite   Sánchez et al. (2016), optical gap at room temperature using diffuse reflectance spectroscopy.
Frontiers in Chemistry | www.frontiersin.org September 2021 | Volume 9 | Article 747972 fundamental gap is obtained by SX-PBE and DSH ( Figures 6C,D), but the VBM is different from that in PBE or the GW method ( Figure 3A). In addition, compared to the GW (LAPW + HLOs) results, the separation between valence Fe-3d and S-3p bands is reduced and the splitting between the valence and conduction Fe-3d bands is significantly increased by the hybrid functionals. We note that both features can be understood tentatively as a result of increased ligand field strength from the perspective of ligand field theory. This indicates an overestimated interaction between the ligand S-3pσ and Fe-3d orbitals in the selected hybrid functionals than that in PBE. The overestimation is most significant in the DSH method ( Figure 6D), where the state of predominant S-3pπ characters along the M-Γ path becomes the VBM and conduction Fe-3d bands are raised beyond 6 eV above the Fermi level. We can observe similar features in marcasite band structures from hybrid functionals, as shown in the lower panel of Figure 6. In the valence band region, the S-3p bands are pushed up relatively to Fe-3d bands compared with PBE and GW. The increase is so significant that the VBM along Γ-X, which is mainly of Fe-3d in PBE and GW, is now of predominant S-3p characters. This also leads to a considerable overlap between the two sets of bands in the energy window 1 ∼ 3 eV below the Fermi level. The conduction bands are also shifted to higher energies. However, the shifts are larger for the conduction Fe-3d bands than for S-3p. For the DSH method as an extreme case, the Fe-3d bands are raised up too high and even separated from the S-3p bands.
The radical failure of DSH invites a close inspection of feasibility of DSH for FeS 2 . As a preliminary exploration to the possible cause, we make a direct comparison between the  inverse static dielectric function used in the DSH with ε M PBE and that from the RPA calculation with the LAPW + HLOs basis in pyrite FeS 2 as a function of the length of wave vector in the longrange limit, i.e. q → 0. The inverse dielectric function corresponding to DSH reads (Liu et al., 2020) We note that Eq. 13 differs from the inverse of Eq. 9 because in the derivation of DSH, the exponential function is replaced by erfc  for more details]. As shown in Figure 7, while DSH overestimates ε −1 and underestimates the screening in the short-wavelength region, i.e. near |G| 0, DSH0 model dielectric function with ε M PBE closely resembles that from RPA calculation, which is similar to the observation by Liu et al. (2020) in transition metal oxides. Hence we consider that the screening effect is reasonably captured in DSH0. Further investigation is needed to understand the cause for the failure of DSH for FeS 2 .

CONCLUSION
In the present study, we have investigated the electronic band structures of two FeS 2 polymorphs, namely pyrite and marcasite, by using methods in different frameworks. With the all-electron many-body GW method implemented in the LAPW framework, we find that by using GW 0 @PBE with the LAPW + HLOs basis, pyrite and marcasite are predicted to have indirect fundamental band gaps of 1.14 and 1.16 eV, respectively. The closeness of band gaps for the two polymorphs agrees with the experimental observation . The pyrite band gap from GW 0 @PBE with LAPW + HLOs is very close to the generally accepted experimental value (Ennaoui et al., 1993) and the corresponding density of states also agrees well with energy distribution curves obtained from the photoelectron spectroscopy measurements (Folkerts et al., 1987;Mamiya et al., 1997). In contrast, with the standard LAPW basis, PBEbased G 0 W 0 and GW 0 both lead to negative QP correction to the PBE fundamental gap, which is rarely observed in LDA/GGAbased G 0 W 0 and GW 0 treatments of semiconductors. The splitting between Fe-3d and S-3p valence bands of pyrite is also significantly overestimated compared to experiment. These artifacts exist not only in calculations with the standard LAPW basis, but also in those with LAPW basis extended by an inadequately converged HLOs (Schena et al., 2013). Therefore in order to eliminate such artifacts, it is instrumental to carefully converge the fundamental band gap with respect to the two controlling parameters, namely n LO and Δl LO . We have further studied electronic band structures of FeS 2 polymorphs with different hybrid functionals, including PBE0, HSE06, the screened exchange SX-PBE and the recently developed DSH functional with system-tuned hybridization parameters. We find that all those methods overestimate the band gaps of the two polymorphs by 0.5 ∼ 1.9 eV compared to the results obtained from G 0 W 0 (LAPW + HLOs). The overestimation by PBE0 and HSE06 as reported in the literature is reproduced in this work. Furthermore, either self-consistent or one-shot DSH method fails to improve over the conventional fixed-parameter hybrid functionals. By comparing the model dielectric function used in DSH with that from RPA calculation with LAPW + HLOs in pyrite, we point out that the failure of DSH may not be caused by the insufficiency of the dielectric model used and therefore requires further investigation. Our investigations clearly show that accurate prediction of electronic band structures of FeS 2 polymorphs poses a stringent test on the state-of-the-art firstprinciples approaches, and the GW method based on semilocal density approximation performs well for this difficult system if it is practiced with well-converged numerical accuracy.
Finally, we note that further work in the following aspects can be done to shed more light onto the band gap problem of FeS 2 in terms of GW and hybrid functional calculations. For one thing, it is possible to build the screened Coulomb interaction W using the KS states from the LAPW calculations and calculate the selfenergy Σ with G from LAPW + HLOs. One can compare it with GW using LAPW to see whether it is the inaccurate band summation in W or G to blame. For another, replacing the PBE with the hybrid functional as starting point will be worthwhile to evaluate the dependence of G 0 W 0 /GW 0 results on initial input for FeS 2 . Particularly, considering the severe overestimation of FeS 2 band gaps by the hybrid functionals, it is of great interest to see whether G 0 W 0 /GW 0 can produce a negative QP correction to the gap from hybrid functional calculations such that the experimental gap is approached from above.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
M-YZ performed the calculations, conducted the analysis and wrote the manuscript. HJ revised the manuscript and supervised all the work. All authors listed contributed to the article and approved it for publication.