Recently, it was shown that the calculation of quasiparticle energies using the *G*_{0}*W*_{0} approximation can be performed without computing explicitly any virtual electronic states, by expanding the Green function and screened Coulomb interaction in terms of the eigenstates of the static dielectric matrix. Avoiding the evaluation of virtual electronic states leads to improved efficiency and ease of convergence of *G*_{0}*W*_{0} calculations. Here, we propose a further improvement of the efficiency of these calculations, based on an approximation of density-density response functions of molecules and solids. The approximation relies on the calculation of a subset of eigenvectors of the dielectric matrix using the kinetic operator instead of the full Hamiltonian, and it does not lead to any substantial loss of accuracy for the quasiparticle energies. The computational savings introduced by this approximation depend on the system, and they become more substantial as the number of electrons increases.

## I. INTRODUCTION

Devising accurate and efficient methods to predict the electronic properties of molecules and condensed systems is an active field of research. Density functional theory (DFT) has been widely used for electronic structure calculations.^{1–3} However, the exact form of the exchange-correlation functional is unknown, and therefore, DFT results depend on the choice of approximate functionals. Improvement over DFT results may be obtained by using many-body perturbation theory (MBPT).^{4,5} A practical formulation of MBPT for many electron systems was proposed by Hedin,^{6} where the self-energy Σ is written in terms of Green’s function *G* and the screened Coulomb interaction *W*.

The *GW* approximation^{5,6} has been successful in the description of the electronic properties of several classes of materials and molecules;^{7–12} however, the computational cost of *GW* calculations remains rather demanding and many complex systems cannot yet be studied using MBPT. Hence, algorithmic improvements are required to apply MBPT to realistic systems. One of the most demanding steps of the original implementation of *GW* calculations^{13–17} involves an explicit summation over a large number of unoccupied single particle electronic orbitals, which enter the evaluation of the dielectric matrix *ϵ*^{18,19} defining the screened Coulomb interaction *W* (*W* = *ϵ*^{−1}$vc$, where $vc$ is the Coulomb interaction). The summation usually converges slowly as a function of the number of virtual orbitals (*N*_{c}).

In recent years, several approaches have been proposed to improve the efficiency of *GW* calculations. For example, in Ref. 20, it was suggested to replace unoccupied orbitals with approximate physical orbitals (SAPOs); the author of Ref. 21 simply truncated the sum over empty states entering the calculation of the irreducible density-density response function and assigned the same, average energy to all the empty states higher than a preset value; in a similar fashion, in Ref. 22, an integration over the density of empty states higher than a preset value was used. Other approaches adopted sophisticated algorithms to invert the dielectric matrix; e.g., in Ref. 23, they employed a Lanczos algorithm. Recently, an implementation of *G*_{0}*W*_{0} calculations avoiding altogether explicit summations over unoccupied orbitals, as well as the necessity to invert dielectric matrices, has been proposed,^{24–28} based on the spectral decomposition of density-density response functions in terms of eigenvectors (also known as projective dielectric eigenpotentials, PDEPs). In spite of the efficiency improvement introduced by such formulation, *G*_{0}*W*_{0} calculations for large systems remain computationally demanding.

In this paper, we propose an approximation to the projective dielectric technique,^{24} which in many cases leads to computational savings of *G*_{0}*W*_{0} calculations of 10%–50%, without compromising accuracy. The rest of the paper is organized as follows: we describe the proposed methodology in Sec. II, and then, we present results for several systems in Sec. III, followed by our conclusions.

## II. METHODOLOGY

We compute the density-density response function of solids and molecules within the framework of the random phase approximation (RPA), using projective dielectric eigenpotentials (PDEP).^{24,25,28} The accuracy of this approach has been extensively tested for molecules and solids.^{25} The technique relies on the solution of Sternheimer’s equation^{3,29,30}

to obtain the linear variation of the *v*th occupied electronic orbital, |Δ$\psi v$⟩, induced by the external perturbation $\Delta V^$. In Eq. (1), *Î* is the identity operator, $P^c$ is the projector onto the unoccupied states, $\epsilon v$ and $\psi v$ are the *v*th eigenvalue and eigenvector of the unperturbed Kohn-Sham Hamiltonian $\u0124=K^+V^SCF$, respectively, where $K^=\u2212\u220722$ is the kinetic energy, $V^SCF$ is the self-consistent potential. For each perturbation, the first order response of the density Δ*n* can be obtained as^{31}

Equations (1) and (2) can be used to iteratively diagonalize the static symmetrized *irreducible* density-density response, $\chi \u03030$,^{24,25,28}

where *λ*_{i} and *ξ*_{i} are eigenvalues and eigenvectors of $\chi \u03030$ and *N*_{PDEP} is the number of eigenvectors of $\chi \u03030$, respectively. The eigenvectors *ξ*_{i} are called PDEPs: projective dielectric eigenpotentials throughout the manuscript. The symmetrized irreducible density-density response function is defined as $\chi \u03030=vc1/2\chi 0vc1/2$, where $vc$ is the Coulomb potential.^{28} Within the RPA, the symmetrized *reducible* density-density response can be expressed as $\chi \u0303=1\u2212\chi \u03030\u22121\chi \u03030$, and therefore, the *ξ*_{i} are also eigenvectors of $\chi \u0303$. The projective dielectric technique has also been recently applied beyond the RPA using a finite field method.^{32}

When solving Sternheimer’s equation, it is not necessary to compute explicitly the electronic empty states because one can write $P^c=\xce\u2212P^v$, since the eigenvectors of *Ĥ* form a complete basis set ($P^v$ is the projector onto the occupied states). The use of Eq. (3) significantly reduces the cost of *G*_{0}*W*_{0} calculation from $Npw2NvNc$ to $NPDEPNpwNv2$, where $Nv$, *N*_{c}, *N*_{PDEP}, and *N*_{pw} are numbers of occupied orbitals (valence bands in solids), virtual orbitals (conduction bands in solids), PDEPs, and plane waves, respectively, and importantly *N*_{PDEP} ≪ *N*_{pw}.

The application of the algorithm described above to large systems is hindered by the cost of solving Eq. (1). However, we note that the eigenvalues of $\chi \u03030$ rapidly converge to zero,^{9,24,25,33} (an example is shown in Fig. 1). In addition, as shown in Ref. 25, the eigenvalue spectrum of the dielectric function for eigenvectors higher than the first few is similar to that of the Lindhard function.^{34} Hence, we propose to compute the PDEPs of $\chi \u03030$ corresponding to the lowest eigenvalues using Eqs. (1) and (2) and to compute the remaining ones with a less costly approach. Inspired by the work of Rocca,^{35} we approximate the eigenpotentials corresponding to higher eigenvalues with kinetic eigenpotentials, which are obtained approximating the full Hamiltonian entering Eq. (1) with the kinetic operator($K^$)^{35,36}

We note that the application of the kinetic energy ($K^$) amounts simply to computing a sum over the plane-wave expansion coefficients multiplied by the square of the **G**-vectors (we are using a plane-wave basis set). Numerous additional operations involving the solution of a Poisson equation and integrals in real space are necessary to apply the full Hamiltonian in Eq. (1), which includes the Hartree and exchange correlation potentials. Furthermore, the evaluation of the nonlocal pseudopotentials [not needed when solving Eq. (4)] is an expensive operation in principle of *O*(*N*^{3}), where *N* is the number of electrons.^{3,30}

In the following, we refer to the eigenpotentials from Eq. (1) as standard PDEPs (stdPDEP, *ξ*_{i}, *i* = 1, …, *N*_{stdPDEP}) and those from Eq. (4) as kinetic PDEPs (kinPDEP, *η*_{i}, *i* = 1, …, *N*_{kinPDEP}), and we rewrite the irreducible density-density response function as

where *ξ*_{i} and *η*_{j} are standard and kinetic PDEPs, respectively, and *λ*_{i} and *μ*_{j} are their corresponding eigenvalues. The procedure to generate stdPDEPs and kinPDEPs is summarized in Fig. 2. We note that during the construction of kinetic PDEPs, the projection operator $P^=\xce\u2212\u2211NstdPDEP\xi i\xi i$ was applied so as to satisfy the orthonormality constraint, ⟨*ξ*_{i}|*ξ*_{j}⟩ = *δ*_{ij}, ⟨*η*_{i}|*η*_{j}⟩ = *δ*_{ij} and ⟨*ξ*_{i}|*η*_{j}⟩ = 0, *∀*(*i*, *j*); in addition, we applied $vc1/2$ to perturbations to yield a *symmetrized* irreducible response function $\chi \u03030$.

In our *G*_{0}*W*_{0} calculations, both the static Green’s function and the statically screened Coulomb interaction are written in the basis of eigenpotentials of the dielectric matrix. Frequency integration is performed using a contour deformation algorithm. A detailed description of the implementation of *G*_{0}*W*_{0} calculations in the basis of eigenpotentials can be found in Refs. 26–28. One should also note that Eqs. (1) and (2) apply only to semiconductors, but this formalism may be generalized to metallic systems.

## III. VALIDATION AND RESULTS

We now turn to discussing results for molecules and solids obtained by using a combination of standard and kinetic PDEPs. To examine the efficiency and applicability of the approximation proposed in Sec. II, we performed *G*_{0}*W*_{0} calculations for a set of closed-shell small molecules, a larger molecule (Buckminsterfullerene C_{60}), and an amorphous silicon nitride/silicon interface (Si_{3}N_{4}/Si) with a total of 1152 valence electrons. All Kohn-Sham eigenvalues and eigenvectors were obtained with the QuantumEspresso package,^{37,38} using the PBE approximation,^{39} SG15^{40} ONCV^{41} pseudopotentials, and *G*_{0}*W*_{0} calculations were carried out with the WEST code.^{28}

We first considered a subset of molecules taken from the G2/97 test set^{42} and calculated their vertical ionization potential (VIP) and electron affinity (EA) using different numbers of stdPDEPs (*N*_{stdPDEP}) and kinPDEPs (*N*_{kinPDEP}). We chose a plane wave cutoff of 85 Ry and a periodic box of edge 30 Bohr. For all molecules, we included either 20 or 100 stdPDEPs in our calculations, and then added 100, 200, 300, 400 kinPDEPs in subsequent calculations, after which an extrapolation was performed ($a+bNstdPDEP+NkinPDEP$) to find converged values (one example is shown in Fig. 3). These results are given in the second (A) and third columns (B) of Tables I and II.^{43} The reference results reported in the last column (C) of the two tables were obtained with 200, 300, 400, 500 stdPDEPs, and an extrapolation was applied. We found that including only 20 stdPDEPs yields quasiparticle energies accurate within 0.1 eV relative to the reference *G*_{0}*W*_{0} values obtained using only standard eigenpotentials (see also Fig. 2 of the supplementary material). The two data sets starting from 20 or 100 stdPDEPs enabled us to save 40% and 10% of computer time compared to the time usage needed with only standard eigenpotentials.

Molecule . | A . | B . | C . |
---|---|---|---|

C_{2}H_{2} | 11.07 | 11.06 | 11.06 |

C_{2}H_{4} | 10.41 | 10.40 | 10.40 |

C_{4}H_{4}S | 8.80 | 8.77 | 8.76 |

C_{6}H_{6} | 9.17 | 9.14 | 9.13 |

CH_{3}Cl | 11.28 | 11.26 | 11.25 |

CH_{3}OH | 10.58 | 10.56 | 10.56 |

CH_{3}SH | 9.39 | 9.36 | 9.36 |

CH_{4} | 14.01 | 14.01 | 14.01 |

Cl_{2} | 11.51 | 11.51 | 11.50 |

ClF | 12.55 | 12.55 | 12.54 |

CO | 13.51 | 13.50 | 13.50 |

CO_{2} | 13.32 | 13.31 | 13.31 |

CS | 11.00 | 10.98 | 10.98 |

F_{2} | 14.99 | 14.97 | 14.97 |

H_{2}CO | 10.43 | 10.42 | 10.42 |

H_{2}O | 11.82 | 11.82 | 11.81 |

H_{2}O_{2} | 10.87 | 10.87 | 10.86 |

HCl | 12.50 | 12.50 | 12.50 |

HCN | 13.20 | 13.20 | 13.20 |

Na_{2} | 4.95 | 4.95 | 4.95 |

Molecule . | A . | B . | C . |
---|---|---|---|

C_{2}H_{2} | 11.07 | 11.06 | 11.06 |

C_{2}H_{4} | 10.41 | 10.40 | 10.40 |

C_{4}H_{4}S | 8.80 | 8.77 | 8.76 |

C_{6}H_{6} | 9.17 | 9.14 | 9.13 |

CH_{3}Cl | 11.28 | 11.26 | 11.25 |

CH_{3}OH | 10.58 | 10.56 | 10.56 |

CH_{3}SH | 9.39 | 9.36 | 9.36 |

CH_{4} | 14.01 | 14.01 | 14.01 |

Cl_{2} | 11.51 | 11.51 | 11.50 |

ClF | 12.55 | 12.55 | 12.54 |

CO | 13.51 | 13.50 | 13.50 |

CO_{2} | 13.32 | 13.31 | 13.31 |

CS | 11.00 | 10.98 | 10.98 |

F_{2} | 14.99 | 14.97 | 14.97 |

H_{2}CO | 10.43 | 10.42 | 10.42 |

H_{2}O | 11.82 | 11.82 | 11.81 |

H_{2}O_{2} | 10.87 | 10.87 | 10.86 |

HCl | 12.50 | 12.50 | 12.50 |

HCN | 13.20 | 13.20 | 13.20 |

Na_{2} | 4.95 | 4.95 | 4.95 |

Molecule . | A . | B . | C . |
---|---|---|---|

C_{2}H_{2} | −2.42 | −2.41 | −2.41 |

C_{2}H_{4} | −1.75 | −1.75 | −1.75 |

C_{4}H_{4}S | −0.85 | −0.81 | −0.80 |

C_{6}H_{6} | −1.01 | −0.96 | −0.96 |

CH_{3}Cl | −1.17 | −1.16 | −1.16 |

CH_{3}OH | −0.89 | −0.89 | −0.89 |

CH_{3}SH | −0.88 | −0.88 | −0.88 |

CH_{4} | −0.64 | −0.64 | −0.64 |

Cl_{2} | 1.65 | 1.64 | 1.65 |

ClF | 1.28 | 1.28 | 1.28 |

CO | −1.56 | −1.57 | −1.57 |

CO_{2} | −0.97 | −0.97 | −0.97 |

CS | 0.49 | 0.51 | 0.51 |

F_{2} | 1.16 | 1.16 | 1.16 |

H_{2}CO | −0.69 | −0.68 | −0.68 |

H_{2}O | −0.90 | −0.90 | −0.90 |

H_{2}O_{2} | −1.80 | −1.79 | −1.79 |

HCl | −1.07 | −1.07 | −1.07 |

HCN | −2.08 | −2.08 | −2.08 |

Na_{2} | 0.64 | 0.63 | 0.63 |

Molecule . | A . | B . | C . |
---|---|---|---|

C_{2}H_{2} | −2.42 | −2.41 | −2.41 |

C_{2}H_{4} | −1.75 | −1.75 | −1.75 |

C_{4}H_{4}S | −0.85 | −0.81 | −0.80 |

C_{6}H_{6} | −1.01 | −0.96 | −0.96 |

CH_{3}Cl | −1.17 | −1.16 | −1.16 |

CH_{3}OH | −0.89 | −0.89 | −0.89 |

CH_{3}SH | −0.88 | −0.88 | −0.88 |

CH_{4} | −0.64 | −0.64 | −0.64 |

Cl_{2} | 1.65 | 1.64 | 1.65 |

ClF | 1.28 | 1.28 | 1.28 |

CO | −1.56 | −1.57 | −1.57 |

CO_{2} | −0.97 | −0.97 | −0.97 |

CS | 0.49 | 0.51 | 0.51 |

F_{2} | 1.16 | 1.16 | 1.16 |

H_{2}CO | −0.69 | −0.68 | −0.68 |

H_{2}O | −0.90 | −0.90 | −0.90 |

H_{2}O_{2} | −1.80 | −1.79 | −1.79 |

HCl | −1.07 | −1.07 | −1.07 |

HCN | −2.08 | −2.08 | −2.08 |

Na_{2} | 0.64 | 0.63 | 0.63 |

Table III shows our results for the C_{60} molecule. The structure of C_{60} (point group *I*_{h}) was also taken from the NIST computational chemistry database,^{44} (optimized with the *ω*B97X-D functional and cc-pVTZ basis sets) and no further optimization was carried out. We used the PBE exchange-correlation functional, a plane wave energy cutoff of 40 *Ry* and cell size of 40 bohrs, the same as used in Ref. 27. We performed two groups of calculations starting with 100 and 200 standard eigenpotentials, respectively. For both calculations, we computed quasiparticle energies by adding 100, 200, 300, 400 kinetic eigenpotentials, and the extrapolation was done in the same manner. As seen in Table III, the results obtained with 200 standard eigenpotentials and additional kinetic eigenpotentials differ at most by 0.1 eV from those computed with standard eigenpotentials (extrapolated up to *N*_{stdPDEP} = 2000). The two sets of calculations starting with *N*_{stdPDEP} = 100 and *N*_{kinPDEP} = 200 amounted to savings of 32% and 15%, respectively.

Energy levels . | N_{stdPDEP} = 100
. | N_{stdPDEP} = 200
. | N_{kinPDEP} = 0
. | G_{0}W_{0}
. | Expt. . |
---|---|---|---|---|---|

t_{1g} | −1.70 | −1.66 | −1.70 | ||

t_{1u} | −2.70 | −2.77 | −2.82 | −2.74,^{a} −2.62,^{b} −2.82^{c} | −2.69^{d} |

h_{u} | −7.32 | −7.32 | −7.38 | −7.31,^{a} −7.21,^{b} −7.37^{c} | −7.61,^{d} −7.6^{e} |

g_{g} + h_{g} | −8.46, −8.52 | −8.46, −8.51 | −8.51, −8.56 | −8.68,^{c} −8.69^{c} | −8.59^{e} |

Energy levels . | N_{stdPDEP} = 100
. | N_{stdPDEP} = 200
. | N_{kinPDEP} = 0
. | G_{0}W_{0}
. | Expt. . |
---|---|---|---|---|---|

t_{1g} | −1.70 | −1.66 | −1.70 | ||

t_{1u} | −2.70 | −2.77 | −2.82 | −2.74,^{a} −2.62,^{b} −2.82^{c} | −2.69^{d} |

h_{u} | −7.32 | −7.32 | −7.38 | −7.31,^{a} −7.21,^{b} −7.37^{c} | −7.61,^{d} −7.6^{e} |

g_{g} + h_{g} | −8.46, −8.52 | −8.46, −8.51 | −8.51, −8.56 | −8.68,^{c} −8.69^{c} | −8.59^{e} |

We now turn to a more complex system, amorphous silicon nitride interfaced with a silicon surface [Si_{3}N_{4}/Si(100)], whose structure was taken from Ref. 48 (see Fig. 4). This interface is representative of a heterogeneous, low dimensional system.

We computed band offsets (BO) by employing two different methods. The first one is based on the calculation of the local density of electronic states (LDOS);^{48,49} the second one is based on the calculation of the average electrostatic potential which is then used to set a common zero of energy on the two parts of the slab representing the two solids interfaced with each other.^{50} The average electrostatic potential was fitted with the method proposed in Ref. 51. We used a plane wave energy cutoff of 70 Ry. We also performed *G*_{0}*W*_{0}@PBE calculations for each bulk system separately and obtained quasiparticle energies.

The local density of states is given by

where *z* is the direction perpendicular to the interface, *ψ*_{i}(*x*, *y*, *z*) is the wavefunction, and the factor 2 represents spin degeneracy. We computed the variation of the valence band maximum (VBM) and conduction band minimum (CBM) as a function of the direction (*z*) perpendicular to the interface^{48}

where *E*_{F} is the Fermi energy and Δ is an constant that is chosen to be 0.003.^{48} We follow a common procedure adopted to describe the electronic structure of interfaces described in Refs. 48 and 49. The band offsets (see Table IV) at the PBE level of theory were determined to be 0.83 eV and 1.49 eV for the valence band and conduction band, respectively, which are in agreement with the results of 0.8 eV and 1.5 eV reported in Ref. 48.

Method∖Energy . | . | VBO . | CBO . | $EgSi$ . | $EgSi3N4$ . |
---|---|---|---|---|---|

PBE | LDOS | 0.83 | 1.49 | 0.67 | 3.19 |

Potential | 0.89 | 1.63 | 0.76 | 3.19 | |

Ref ^{a} | 0.8 | 1.5 | 0.7 | 3.17 | |

G_{0}W_{0} | LDOS | 1.41 | 1.88 | 1.29 | 4.77 |

Potential | 1.46 | 2.02 | 1.29 | 4.77 | |

Ref ^{a} | 1.5 | 1.9 | 1.3 | 4.87 | |

Expt. | 1.5–1.78^{b} | 1.82–2.83^{c} | 1.17^{d} | 4.5–5.5^{e} |

Method∖Energy . | . | VBO . | CBO . | $EgSi$ . | $EgSi3N4$ . |
---|---|---|---|---|---|

PBE | LDOS | 0.83 | 1.49 | 0.67 | 3.19 |

Potential | 0.89 | 1.63 | 0.76 | 3.19 | |

Ref ^{a} | 0.8 | 1.5 | 0.7 | 3.17 | |

G_{0}W_{0} | LDOS | 1.41 | 1.88 | 1.29 | 4.77 |

Potential | 1.46 | 2.02 | 1.29 | 4.77 | |

Ref ^{a} | 1.5 | 1.9 | 1.3 | 4.87 | |

Expt. | 1.5–1.78^{b} | 1.82–2.83^{c} | 1.17^{d} | 4.5–5.5^{e} |

As mentioned above, another method to obtain the valence band offset (VBO) and conduction band offset (CBO) is to align energy levels with respect to electrostatic potentials. Following Ref. 51, the electrostatic potential was computed as

where $V\xafat(i)$ is the potential near the core region obtained from neutral atom calculations. With this method, VBO and CBO at the PBE level are found to be 0.89 eV and 1.63 eV.

To compute *G*_{0}*W*_{0} corrections on band offsets, we performed *G*_{0}*W*_{0}@PBE calculations of bulk silicon and amorphous silicon nitride. In Table V, quasiparticle corrections to Kohn-Sham energies of bulk silicon and amorphous silicon nitride are shown. The second and third columns are computed with 1000 and 2000 standard eigenpotentials. The fitted *G*_{0}*W*_{0} reference results are extrapolated with 500, 1000, 1500, and 2000 standard eigenpotentials. To test the accuracy of kinetic eigenpotentials, we started with 400 stdPDEPs and added 100, 200, 300, 400 kinPDEPs, after which the same extrapolation was applied. We calculated VBO and CBO at the *G*_{0}*W*_{0} level by applying quasiparticle corrections on PBE results. After applying quasiparticle corrections on LDOS results, VBO and CBO are 1.41 eV and 1.88 eV, while the VBO and CBO are found to be 1.46 eV and 2.02 eV after applying corrections to results based on electrostatic potential alignment. Both of them are close to the range of experimental results of 1.5–1.78 eV and 1.82–2.83 eV. Time saving when using kinetic eigenpotentials to obtain quasiparticle corrections was approximately ∼50%.

. | N_{stdPDEP} = 1000
. | N_{stdPDEP} = 2000
. | Fit . | N_{kinPDEP} = 400
. |
---|---|---|---|---|

Si VBM | 5.70 | 5.55 | 5.45 | 5.53 |

Si CBM | 7.03 | 6.91 | 6.79 | 6.82 |

a −Si_{3}N_{4} VBM | 7.14 | 7.01 | 7.01 | 6.99 |

a −Si_{3}N_{4} CBM | 11.99 | 11.87 | 11.83 | 11.83 |

. | N_{stdPDEP} = 1000
. | N_{stdPDEP} = 2000
. | Fit . | N_{kinPDEP} = 400
. |
---|---|---|---|---|

Si VBM | 5.70 | 5.55 | 5.45 | 5.53 |

Si CBM | 7.03 | 6.91 | 6.79 | 6.82 |

a −Si_{3}N_{4} VBM | 7.14 | 7.01 | 7.01 | 6.99 |

a −Si_{3}N_{4} CBM | 11.99 | 11.87 | 11.83 | 11.83 |

## IV. CONCLUSION

The method introduced in Refs. 9 and 26–28 to compute quasiparticle energies using the *G*_{0}*W*_{0} approximation avoids the calculation of virtual electronic states and the inversion and storage of large dielectric matrices, thus leading to substantial computational savings. Building on the strategy proposed in Refs. 24 and 25 and implemented in the WEST code,^{28} here, we proposed an approximation of the spectral decomposition of dielectric matrices that further improve the efficiency of *G*_{0}*W*_{0} calculations. In particular, we built sets of eigenpotentials used as a basis to expand the Green function and the screened Coulomb interaction by solving two separate Sternheimer equations: one using the Hamiltonian of the system to obtain the eigenvectors corresponding to the lowest eigenvalues of the response function, and the other equation using just the kinetic energy operator to obtain the eigenpotentials corresponding to higher eigenvalues. We showed that without compromising much accuracy, this approximation reduces the cost of *G*_{0}*W*_{0} calculations by 10%–50%, depending on the system, with the most savings observed for the largest systems studied here.

## SUPPLEMENTARY MATERIAL

See the supplementary material for convergence studies of the *G*_{0}*W*_{0} calculations of the vertical ionization potential of the CH_{4} molecule, which is taken as a representative example of the molecular systems studied in the main text.

## ACKNOWLEDGMENTS

We thank He Ma, Ryan L. McAvoy, and Ngoc Linh Nguyen for useful discussions. This work was supported by MICCoM, as part of the Computational Materials Sciences Program funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division, through Argonne National Laboratory, under Contract No. DE-AC02-06CH11357. This research used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract No. DE-AC02-06CH11357, and resources of the University of Chicago Research Computing Center.

## REFERENCES

_{0}W

_{0}for molecular systems

Here, we did not attempt to solve the problem on how to accurately compute resonant molecular energy levels: our goal is to compare results obtained with solutions of the Sternheimer equation using the full Hamiltonian [Eq. (1)] and approximate solutions using only the kinetic operator [Eq. (4)]. As long as the results obtained with the two procedures agree, we consider the results obtained from Eq. (4) as accurate.

_{60}, C

_{70},[C

_{60}]PCBM, and bis-[C

_{60}]PCBM using a many-body G

_{0}W

_{0}-Lanczos approach

_{60}, buckminsterfullerene

_{3}N

_{4}/Si(100) interface: A first-principles study

_{2}/Si(001) interfaces

_{2}and Si

_{3}N

_{4}films on Si(111) and Si(100) from photoemission spectroscopy

_{3}N

_{4}/Si interface formed using nitrogen-hydrogen radicals

_{3}N

_{4}