Towards fully automatized GW band structure calculations: What we can learn from 60.000 self-energy evaluations

We analyze a data set comprising 370 GW band structures composed of 61716 quasiparticle (QP) energies of two-dimensional (2D) materials spanning 14 crystal structures and 52 elements. The data results from PAW plane wave based one-shot G$_0$W$_0$@PBE calculations with full frequency integration. We investigate the distribution of key quantities like the QP self-energy corrections and renormalization factor $Z$ and explore their dependence on chemical composition and magnetic state. The linear QP approximation is identified as a significant error source and propose schemes for controlling and drastically reducing this error at low computational cost. We analyze the reliability of the $1/N_\text{PW}$ basis set extrapolation and find that is well-founded with narrow distributions of $r^2$ peaked very close to 1. Finally, we explore the validity of the scissors operator approximation concluding that it is generally not valid for reasonable error tolerances. Our work represents a step towards the development of automatized workflows for high-throughput G$_0$W$_0$ band structure calculations for solids.


I. INTRODUCTION
In computational materials science, the highthroughput (HT) mode of operation is becoming increasingly popular 1 . The development of automatized workflow engines capable of submitting, controlling and receiving thousands of interlinked calculations [2][3][4] with minimal human intervention has greatly expanded the range of materials, and properties, that can be investigated by a single researcher. Several HT studies have been conducted over the past decade mostly with the aim of identifying new prospect materials for various applications including catalysis 5 , batteries 6,7 , thermoelectrics 8,9 , photocatalysts 10 , transparent conductors 11 , and photovoltaics 12,13 , just to mention some. The vast amounts of data generated by such screening studies have been stored in open databases [14][15][16][17] making them available for further processing, testing and comparison of methods and codes, training of machine learning algorithms etc. With very few exceptions, the HT screening studies and the generation of the materials databases, have been based on density functional theory (DFT) at the level of the generalized gradient approximation (GGA).
While DFT is fairly accurate for structural parameters and other properties related to the electronic ground state, it is well known that electronic band structures, in particular the size of band gaps, are not well reproduced by most xc-functionals 18 . In fact, as of today quantitatively accurate band structures of solids can only be obtained from many-body methods like the GW approximation [19][20][21] , which explicitly accounts for exchange and dynamical screening effects. In fact, while (standard) DFT easily underestimates the experimental band gap by 1 eV or more, the GW method even in its simplest non-selfconsistent G 0 W 0 flavour, is accurate to within 0.4 eV for typical solids [22][23][24] . We note in passing that for partially self-consistent GW 0 22 or when vertex corrections are included 25,26 , the deviation from experiments falls below 0.2 eV, which is comparable to the accuracy of the experimental band gaps. The improved accuracy of the GW method(s) comes at the price of a significantly more involved methodology both conceptually and numerically as compared to DFT. While DFT calculations can be routinely performed by non-experts using codes that despite very different numerical implementations produce identical results 27 , GW calculations remain an art for the expert.
The high complexity of GW calculations is due to several factors including: (i) The basic quantities of the theory, i.e. the Greens function (G) and screened Coulomb interaction (W ) are dynamical quantities which depend on time/frequency. Several possibilities for handling the frequency dependence exists including the formally exact direct integration 23 and contour deformation techniques 28 as well as the controlled approximate analytic continuation methods 29 and the rather uncontrolled but inexpensive plasmon-pole approximations 19 . (ii) The formalism involves infinite sums over the unoccupied bands. While most implementations perform the sum explicitly up to a certain cut-off, schemes to avoid the sum over empty states have been developed 30,31 . (iii) The basic quantities are two-point functions in real space (or reciprocal space) that couple states at different k-points. This leads to large memory requirements and makes it unfeasible to fully converge GW calculations with respect to basis set. Consequently, strategies for extrapolation to the infinite basis set limit must be employed 32,33 . (vi) Unless the GW equations are solved fully self-consistently, which is rarely done and does not improve accuracy 22,26 , there is always a starting point dependence. This has been systematically explored for molecules where it was found that LDA/GGA often comprise a poor starting point whereas hybrids per-form better in the sense that they lead to better agreement with experimental ionization potentials and produce more well defined spectral peaks with higher quasiparticle weights 34,35 . These and other factors imply that GW calculations not only become significantly more demanding than DFT in terms of computer resources, but they also involve more parameters making it difficult to assess whether the obtained results are properly converged or perhaps even erroneous.
Successful application of the HT approach to problems involving excited electronic states, e.g. light absorption/emission, calls for development of automatized and robust algorithms for setting the parameters of manybody calculations such as GW (according to available computational resources and required accuracy level), extrapolating the basis set, and assessing the reliability of the obtained results. The first step towards this goal is to analyse and systematize the data from large-scale GW studies. With a similar goal in mind van Setten et al. compared G 0 W 0 @PBE band gaps, obtained with the plasmon-pole approximation, to the experimental band gaps. They analyzed the correlations between different quantities and concluded that that G 0 W 0 (with plasmonpole approximation) is more accurate than using an empirical correction of the PBE gap, but that, for accurate predictive results for a broad class of materials, an improved starting point or some type of self-consistency is necessary.
In this work we perform a detailed analysis of an extensive GW data set consisting of G 0 W 0 @PBE band structures of 370 two-dimensional semiconductors comprising a total of 61716 QP energies. Our focus is not on the ability of the G 0 W 0 to reproduce experiments, i.e. its accuracy, which is well established by numerous previous studies, but rather on the numerical robustness and reliability of the method and the basis set extrapolation procedure. The calculations employ a plane wave basis set and direct frequency integration; thus the use of projector augmented wave (PAW) potentials represents the only significant numerical approximation. We investigate the distribution of self-energy corrections and renormalization factors, Z, and explore their dependence on the materials composition and magnetic state. By investigating the full frequency dependent self-energy for selected materials we analyse the error caused by the linear approximation to the QP equation and propose methods to estimate and correct this error. We assess the reliability of a plane wave basis set extrapolation scheme finding it to be very accurate with r 2 values above 0.95 in more than 90% of the cases when extrapolation is performed from 200 eV. Finally, we assess the accuracy of the scissors operator approach, and conclude that it should only be used when average (maximal) band energy errors of 0.2 eV (2 eV) are acceptable.

A. The G0W0 data set
The 370 G 0 W 0 calculations were performed as part of the Computational 2D Materials Database (C2DB) project 36 . Below we briefly recapitulate the computational details behind the G 0 W 0 calculations and refer to Ref. 36 for more details. All calculations were performed with the projector augmented wave function code GPAW 37 .
The C2DB database contains around 4000 monolayers comprising both known and hypothetical 2D materials constructed by decorating experimentally known crystal prototypes with a subset of elements from the periodic table 36 . Currently, G 0 W 0 calculations have been performed for 370 materials spanning 14 different crystal structures and 52 different chemical elements. Fig. 1a illustrates the distribution of elements. The number of materials containing a given element is shown below the element symbol. The number of magnetic materials containing the elements is shown in a parenthesis next to the total number.
To give an overview of some of the data analysed in this work, the distribution of the 61716 G 0 W 0 corrections for the six bands around the band gap is shown in figure  1 (b). The distribution for the valence bands is shown in blue and for the conduction in orange. It is usually the case in GW studies that the DFT valence bands are shifted down and the conduction bands are shifted up. A similar behavior is found for the main part of our data, but we also observe a small subset of states for which the correction has the opposite sign. Figure 1 (c) shows a scatter plot of the PBE energies versus the G 0 W 0 energies. We only show energies from -10 eV to 10 eV for clarity. The color of a point shows the Z value. The latter has been truncated to the region [0.5, 1.0] to show the variation of the main part of the distribution. The main observation we can make from this figure is that there is no obvious correlation between the energies and the Z values. This is also verified by the calculated correlation coefficient between E PBE and Z (= 0.27), E G0W0 and Z (= 0.23) and between the G 0 W 0 correction, E G0W0 − E PBE , and Z (= 0.10). We conclude that there is no significant correlation between the energies and Z, meaning that low Z values (which signals a break down of the QP approximation) may occur in any energy range.

B. Quasiparticle weight Z
The quasiparticle weight, Z, gives a rough measure of the validity of the quasiparticle picture, i.e. how well the charged excitations of the interacting electron system can be described by single-particle excitations from the ground state. In the Methods section we prove a physical interpretation of the quasiparticle weight. In the following we analyse the 61716 calculated quasiparticle weights, Z, contained in the C2DB database. As discussed in the Methods section, for the QP approximation to be well-founded Z should be close to 1. We split the Z values into two groups quasiparticle-consistent (QP-c) for Z ∈ [0.5, 1.0] and quasiparticle-inconsistent (QP-ic) for Z ∈ [0.5, 1.0]. We can expect that the QP approximation works well for QP-c states and worse for QP-ic states. Figure 2 shows a histogram of the Z-values (all extrapolated to the infinite plane wave limit) corresponding to the 3 highest valence bands and 3 lowest conduction band of 370 semiconductors. The vast majority of the values are distributed around ≈ 0.75 with only 0.28% lying outside the physical range from 0 to 1 (0.16% are larger than one and 0.12% are negative). We find that 97.5% of states are in QP-c.
It is of interest to investigate if there are specific types of materials/elements that are particularly challenging to describe by G 0 W 0 . Figure 3 shows a barplot of the percentage of QP-ic states in materials containing a given element (note the logarithmic scale). The result of this analysis performed on the non-magnetic (ferromagnetic) materials is shown in blue (orange). For example, a large percentage (about 65%) of the states in Co-containing materials are QP-ic. It is clear that magnetic materials contribute a large fraction of the QP-ic eigenstates. In fact, 0.36% of the non-magnetic states are QP-ic while 22% of the magnetic states are QP-ic. In general it thus seems that the QP approximation with KS eigenstates as the reference state is generally worse for magnetic materials.
Based on the distribution of QP weights in figure 2, it appears that the QP approximation is valid for essen-   tially all the states in the non-magnetic materials and most of the states in the magnetic materials. However, while a QP-c Z value is likely a necessary condition for predicting an accurate QP energy from the linearized QP equation [Eq. (6) in the Method section], it is not sufficient. This is because the assumption behind Eq. (6), i.e. that Σ(ε) varies linearly with ε in the range between the KS energy and the QP energy, is not guaranteed for QP-c states. This is illustrated in figure 4 which shows the full frequency dependent self-energy for three states in the ferromagnetic FeCl 2 . Case (a) is a typical example where the self-energy of a QP-c state (Z = 0.61) varies linearly around ε KS and the 1st order approximation works well. The second case (b) shows an example where the 1st order approximation breaks down for a QPic state (Z = 1.19). The final case (c) illustrates that the 1st order approximation can break down even in cases where Z is very close to 1. Unfortunately, there is no simple way to diagnose such cases from the information available in a standard G 0 W 0 calculation (Σ(ε KS ) and Z). We stress that the example in figure 4(c) is a special case and that in general, the linear approximation is significantly more likely to hold for QP-c states than for QP-ic states (see discussion below).

C. Beyond the linear QP approximation
Under the assumption that the KS wave functions constitute a good approximation to the QP wave functions, so that off-diagonal elements can be neglected, the solution to the QP equation reduces to solving an equation of the form  In this section we investigate different root-finding schemes to estimate the size of the error introduced by the linear approximation and obtain an improved QP energy. With high-throughput computations in mind, a good algorithm provides a reasonable balance between computation time (number of Σ/Z evaluations) and accuracy. To benchmark the different schemes we computed the full frequency dependent self-energy for 3192 states, corresponding to the 3 highest valence bands and 3 lowest conduction bands, for 12 of the 370 2D materials (including two ferromagnetic materials). The self-energy is evaluated on a uniform frequency grid and interpolated using cubic splines. The "true" solution of the QP equation is then determined and used to evaluate the errors of the approximate schemes. We first consider the iterative Newton-Raphson (NR) method where we limit ourselves to 1 and 2 iterations to keep the number of self-energy evaluations and thus the computational cost low. We note that 1 iteration (NR1) is equivalent to the linear approximation. The distribution of the errors is shown in figure 6 (a). Although 87% of the errors from NR1 are below 0.1 eV, the mean absolute error (MAE) is 0.11 eV due to outliers. Most of these errors are significantly reduced by performing one more iteration of Newton-Raphson (NR2), but again large outliers pull the MAE up. If we evaluate the MAE without the outliers (those lying outside the displayed error range), the MAE reduces to only 0.006 eV.
Motivated by the relatively narrow distribution of Z values in figure 2 we consider an empirical solution estimate consisting of replacing the actual Z value with the mean value of the distribution, i.e. we simply set Z = 0.75. This has the advantage of being simple, computationally cheap, and robust in the sense of avoiding outlier Z-values arising from local irregularities in Σ at the KS energy, see figure 4(b). The resulting error distribution is shown in 6 (b). While the central part of   the distribution is slightly broadened compared to the 1st order approximation, the MAE is reduced due to a reduction of outliers (enhanced robustness). As shown in panel (c), the central part of the distribution can be narrowed by applying the empirical approach only for QP-ic states, i.e. when Z ∈ [0.5, 1]. In fact, this approach (empZ@QP-ic) has a MAE equal to that of NR2 but with half the computational cost (two Σ/Z evaluations compared to four). Next, we examine polynomial fitting of the self-energy. We construct second and fourth order polynomials, P n (ω), from the self-energy at energies in a range of ±1 eV around the KS energy. The cost of the second and fourth order fits are equivalent to three and five self-energy evaluations, respectively. In general, the polynomial fits have rather low correlation coefficients of C < 0.9 and are sensitive to the choice of frequency points and self-energy data used for the fit. As a consequence the resulting errors are large (not shown) and the approach is not suitable. We attribute this to our observation that self-energies are often irregular (on the relevant scale of 1 eV) and not well-described by loworder polynomials.
Finally, we consider a scheme that we refer to as ΣdE, which estimates the error as The motivation for this expression is the following. If the linear approximation is exact, then δ vanishes as it should. Moreover, if the self-energy has a non-zero curvature it can show that δ equals the true error to leading order in the curvature. In that sense it is similar to the second order polynomial fit, but with the important difference that whereas the polynomial fit was based on uniformly distributed points, ΣdE uses the value and slope at E KS and the value at E QP, lin .
The errors resulting from Eq. (2) are shown in 7(b). Compared to the linear approximation, the ΣdE reduces the MAE from 0.11 eV to 0.05 eV, at the cost of one additional self-energy evaluation. Interestingly, Eq. (2) systematically overestimates the error. This can be seen from figure 7(a), which shows the ratio of the estimated and true error. A Gaussian fit to the distribution (red curve) has a mean value of α 0 = 1.5 and standard deviation of 0.2. Based on this observation it is tempting to correct for the systematic error using α = α 0 . As can be seen from 7(c), this corrected ΣdE scheme shows excellent performance with an almost four-fold reduction of the MAE from 0.11 eV for the linear approximation to only 0.03 eV at a computational overhead of just one additional self-energy evaluation.

D. Plane wave extrapolation
The self-energy and the derivative of the self-energy (both evaluated at the KS energy) are calculated at three cutoff energies: 170 eV, 185 eV, and 200 eV. These values are then extrapolated to infinite cutoff, or infinite number of plane waves, by assuming a linear dependence on the inverse number of plane waves 38 . An example of this fitting procedure is shown in figure 8 (a). The extrapolation procedure saves computational time while improving the accuracy of the results -provided the extrapolation is sufficiently accurate. Extrapolation can fail if convergence as a function of plane wave cutoff for the given quantity does not follow the expected 1/N PW behaviour in the considered cutoff range.
To validate this approach we investigate the distribution of the r 2 values (coefficient of determination) for all 61716 extrapolations in C2DB. We split them into two cases: extrapolation of the self-energy and extrapolation of the derivative of the self-energy. The distributions are shown as histograms in figure 8 (b). The distributions are clearly peaked very close to 1, and in general it seems that the extrapolation is very good. The distribution for the derivatives is somewhat broader, and the extrapolation is generally less accurate than for the self-energies, which indicates a slower convergence with plane waves than for the self-energies. If we choose r 2 = 0.8 as an acceptable threshold, we find that 1.7% of the r 2 values of the selfenergy extrapolation fall below this criterion while 5.0% are below for the derivative extrapolation. While these numbers might seem large, the problem is readily diagnosed (by the r 2 value) and can be alleviated by using higher plane wave cutoffs.

E. Scissors operator approximation
Within the so-called scissors operator approximation (SOA) it is assumed that the G 0 W 0 correction is independent of band-and k-index. Consequently, the G 0 W 0 correction calculated at e.g. the Γ point is applied to all the eigenvalues thus saving computational time as only one G 0 W 0 correction is required. In figure 9 (a) the idea is illustrated for a generic band. With the notation from the figure, the SOA consists of setting ∆(k) = ∆ (or ∆ nσ (k) = ∆ nσ when more than one band and spin is involved).
To test the accuracy of the SOA, we evaluate the mean absolute error ( | | ) and maximum absolute error (max(| |)) of the band energies obtained with the SOA for each of the 370 materials: and The distribution of these errors are shown in figure 9 (b) and (c). From panel (b) we see that the mean error exceeds 100 meV for about half of all materials -a rather large error, comparable to the target accuracy of the G 0 W 0 method itself. Furthermore, it follows from (c) that the maximum absolute error is often 0.5−1.0 eV. We conclude that while the average error of the SOA might be acceptable, it can produce significant errors for specific bands and should be used with care. A compromise between the SOA and calculating the G 0 W 0 energies at every k-point could be to calculate the G 0 W 0 correction at several k-points and fit ∆ nσ (k) to the resulting values.

III. DISCUSSION
As high-throughput computations are gaining popularity in the electronic structure community it becomes important to establish protocols for performing various types of calculations in an automated, robust and errorcontrolled manner. In this work we have taken the first steps towards the development of automated workflows for G 0 W 0 band structure calculations of solids. With G 0 W 0 representing the state-of-the-art for predicting QP energies in condensed matter systems, such workflows are essential for continued progress in the field of computational materials design.
Based on our detailed analysis of 61716 G 0 W 0 selfenergy evaluations for the eigenstates of 370 twodimensional semiconductors we were able to draw several conclusions relevant to large-scale GW studies. First of all, we found it useful to divide the states into two categories, namely quasiparticle consistent (QP-c) and quasiparticle inconsistent (QP-ic) states defined by Z ∈ [0.5, 1.0] and Z ∈ [0.5, 1.0], respectively. Importantly, we found that the QP energies obtained from the standard linearized QP equation are significantly more accurate for QP-c states than for QP-ic state. Moreover, we found the fraction of QP-ic states to be much larger in magnetic materials (22%) than in non-magnetic materials (0.36%). Thus extra care must be taken when performing G 0 W 0 calculations for magnetic materials; in particular, such materials might require a special treatment in high-throughput workflows.
The mean absolute error (MAE) on the QP energies resulting from the linearized QP equation was found to be 0.11 eV on average. The MAE evaluated separately for QP-c and QP-ic states was 0.04 eV and 0.27 eV, respectively. In comparison, the accuracy of the GW approximation itself (compared to experiments) is on the order of 0.2 eV. It is therefore of interest to reduce or at least estimate the numerical error bar on the QP energies obtained from G 0 W 0 calculations. We found that an empirical scheme, where we set Z = 0.75 (corresponding to the mean of the Z-distribution) for QP-ic states, reduces the MAE from 0.11 eV to 0.06 eV with no computational overhead. Similarly, the method dubbed the corrected ΣdE scheme reduces the MAE to 0.03 eV, at the cost of one additional self-energy evaluation. From these studies it seems natural to accompany the QP energies obtained from G 0 W 0 with estimated error bars derived from one of the these correction schemes. In fact, we have used the empZ@QP-ic method to correct all the GW the band structures in the C2DB database.
Our analysis of the well known and widely used scissors operator approximation shows that the errors introduced on the individual QP energies when averaged over all bands (specifically the 3 highest valence and 3 lowest conduction bands) typically is on the order of 0.1 eV while the maximum error typically exceeds 1 eV. We stress that our scissors operator fits each of the six bands separately using the G 0 W 0 corrections at the Γ-point. Thus the errors introduced by the more standard scissors approximation that fits only the band gap, are expected to be even larger. We conclude that the scissors operator should be used with care and only in cases where errors on specific band energies of 1-3 eV are acceptable.
Finally, the 1/N PW extrapolation scheme was found to be highly reliable for our PAW calculations when applied to cutoff energies in the range 180-200 eV. In fact only 1.7% (5.0%) of the self-energy (derivative of self-energy) extrapolations had an r 2 below 0.8. However, for the purpose of high-throughput studies it may be prudent to store and make available information on the r 2 for the extrapolation so that the quality of the extrapolation can always be examined and improved calculations with higher cutoff can be performed if deemed necessary.

A. G0W0 Calculations
For the materials considered here, DFT calculation using PBE 39 were performed using an 800 eV plane-wave cutoff. Spin-orbit coupling is included by diagonalizing the spin-orbit Hamiltonian in the k-subspace of the Bloch states found from PBE.
Those materials that have a finite gap and up to 5 atoms in the unit cell are selected for G 0 W 0 calculations. The QP energies are calculated for the 8 highest occupied and the 4 lowest unoccupied bands. Three energy cutoffs are used: 170 eV, 185 eV, and 200 eV. The results are then extrapolated to infinite energy, i.e. to an infinite number of plane-waves. This extrapolation is done by expressing the self-energies in terms of the inverse number of plane-waves, 1/N PW , performing a linear fit, and determining the value of the fit at 1/N PW = 0, see Refs. 40 and 41.
The screened Coulomb interaction entering in the selfenergy is calculated using full frequency integration in real frequency space. To avoid effects from the (artificially) repeated layers. A Wigner-Seitz truncation scheme is used for the exchange part of the self-energy 42 and a 2D truncation of the Coulomb interaction is used for the correlation part 38,43 . A truncated Coulomb interaction leads to significantly slower k-point convergence because the dielectric function strongly depends on q around q = 0; this is remedied by handling the integral around q = 0 analytically, see 44,45 . A k-point density of 5.0/Å −1 was used.
The statistical analyses performed here use the data from all spins, k-points, and the three highest occupied bands and the three lowest unoccupied bands. In section II B we consider several examples of the full frequencydependent self-energies for a randomly selected spin, kpoint, and band combination, subject to some requirements on the quasi-particle weight, Z, which are described below.

B. Quasiparticle theory
The G 0 W 0 quasi-particle energies are found by solving the quasi-particle equation (QPE) 22 : Here ψ nkσ is the Kohn-Sham wavefunction for band n, crystal momentum k, and spin σ, H KS is the singleparticle Kohn-Sham Hamiltonian, v xc is the exchangecorrelation potential and Σ is the self-energy. Typically, and in C2DB, the QPE is solved via one iteration of the Newton-Raphson method starting from the KS energy, nkσ , which is equivalent to making a linear approximation of the self-energy. This yields the solution E QP nkσ ≈ nkσ + ZRe [ ψ nkσ |Σ( nkσ ) − v xc |ψ nkσ ] , (6) Z is known as the quasi-particle weight. We define the G 0 W 0 correction, ∆E nkσ , as the difference between G 0 W 0 energies and KS energies E QP nkσ − nkσ at spin σ, crystal momentum k, and band n.

Physical Interpretation of Z
Following Ref. 44 we provide here a physical interpretation of Z. We denote the many-body eigenstates for the N particle system by |Ψ N i , where i is the excitation index. An interesting question is how well the state |Ψ N +1 i can be described as the addition of a single electron to the ground state |Ψ N 0 . In other words, can we find an state φ such that |Ψ N +1 i ≈ c † φ |Ψ N 0 ? The optimal φ is determined from maximizing the overlap, i.e.
If the maximal overlap is close to 1 the excited manybody state is well approximated by a single-particle excitation. It turns out that the square of this maximal overlap is exactly equal to the QP weight Z defined by Eq. (6) if it is evaluated at the true QP energy and with the true QP wave function rather than at the KS energy and with the KS wave function. Furthermore Z can be shown to be equal to the squared norm of the QP wave function, which is defined as For a proof of these results we refer to Ref. 44. In standard G 0 W 0 calculations, the self-energy is evaluated at the KS energy using KS eigenstates. In this case, Z is no longer equal to the exact QP weight but only approximates it. If Z deviates significantly from 1, we can only conclude that either 1) the system is strongly correlated so that the QP approximation fails, or 2) the Kohn-Sham energy and/or wave function are a bad approximation to the true QP energy and/or wave function. In either case we would expect that the G 0 W 0 calculation is problematic and requires special attention.

VII. AUTHOR CONTRIBUTIONS
AR performed the statistical analyses and full, frequency-dependent self-energy calculations. TD performed the G 0 W 0 calculations. All authors interpreted the analyses and wrote the article.

VIII. COMPETING INTERESTS
The authors declare no competing interests.