Recently, it was shown that the calculation of quasiparticle energies using the G0W0 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 G0W0 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.

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 approximation5,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 calculations13–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 = ϵ−1vc, where vc is the Coulomb interaction). The summation usually converges slowly as a function of the number of virtual orbitals (Nc).

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 G0W0 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, G0W0 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 G0W0 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.

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 equation3,29,30

(ĤεvÎ)|Δψv=P^cΔV^|ψv
(1)

to obtain the linear variation of the vth occupied electronic orbital, |Δψv⟩, induced by the external perturbation ΔV^. In Eq. (1), Î is the identity operator, P^c is the projector onto the unoccupied states, εv and ψv are the vth eigenvalue and eigenvector of the unperturbed Kohn-Sham Hamiltonian Ĥ=K^+V^SCF, respectively, where K^=22 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 as31 

Δn=2vψvΔψv+c.c..
(2)

Equations (1) and (2) can be used to iteratively diagonalize the static symmetrized irreducible density-density response, χ̃0,24,25,28

χ̃0=i=1NPDEPξiλiξi,
(3)

where λi and ξi are eigenvalues and eigenvectors of χ̃0 and NPDEP is the number of eigenvectors of χ̃0, respectively. The eigenvectors ξi are called PDEPs: projective dielectric eigenpotentials throughout the manuscript. The symmetrized irreducible density-density response function is defined as χ̃0=vc1/2χ0vc1/2, where vc is the Coulomb potential.28 Within the RPA, the symmetrized reducible density-density response can be expressed as χ̃=1χ̃01χ̃0, and therefore, the ξi are also eigenvectors of χ̃. 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=ÎP^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 G0W0 calculation from Npw2NvNc to NPDEPNpwNv2, where Nv, Nc, NPDEP, and Npw are numbers of occupied orbitals (valence bands in solids), virtual orbitals (conduction bands in solids), PDEPs, and plane waves, respectively, and importantly NPDEPNpw.

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 χ̃0 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 χ̃0 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

(K^εvÎ)|Δψv=P^cΔV^|ψv.
(4)

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(N3), where N is the number of electrons.3,30

FIG. 1.

First 500 eigenvalues λi of the symmetrized irreducible density-density response function χ̃0 (see the text), for three small molecules: CH4 (blue dots), C2H4 (orange up triangles), and C2H2 (green down triangles). Nv is the number of occupied orbitals.

FIG. 1.

First 500 eigenvalues λi of the symmetrized irreducible density-density response function χ̃0 (see the text), for three small molecules: CH4 (blue dots), C2H4 (orange up triangles), and C2H2 (green down triangles). Nv is the number of occupied orbitals.

Close modal

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

χ̃0=i=1NstdPDEPξiλiξi+j=1NkinPDEPηjμjηj,
(5)

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^=ÎNstdPDEPξiξ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 χ̃0.

FIG. 2.

The workflow used in this work to generate eigenvectors of the dielectric matrix using the Kohn-Sham Hamiltonian (stdPDEP) and using the kinetic operator (kinPDEP). See the text.

FIG. 2.

The workflow used in this work to generate eigenvectors of the dielectric matrix using the Kohn-Sham Hamiltonian (stdPDEP) and using the kinetic operator (kinPDEP). See the text.

Close modal

In our G0W0 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 G0W0 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.

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 G0W0 calculations for a set of closed-shell small molecules, a larger molecule (Buckminsterfullerene C60), and an amorphous silicon nitride/silicon interface (Si3N4/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 SG1540 ONCV41 pseudopotentials, and G0W0 calculations were carried out with the WEST code.28 

We first considered a subset of molecules taken from the G2/97 test set42 and calculated their vertical ionization potential (VIP) and electron affinity (EA) using different numbers of stdPDEPs (NstdPDEP) and kinPDEPs (NkinPDEP). 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 G0W0 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.

FIG. 3.

Extrapolation of G0W0 energy of highest occupied orbital of methane with respect to total number of eigenpotential used (NstdPDEP + NkinPDEP). In this plot, NstdPDEP = 20 and NkinPDEP = 100, 200, 300, 400 for the four points, respectively.

FIG. 3.

Extrapolation of G0W0 energy of highest occupied orbital of methane with respect to total number of eigenpotential used (NstdPDEP + NkinPDEP). In this plot, NstdPDEP = 20 and NkinPDEP = 100, 200, 300, 400 for the four points, respectively.

Close modal
TABLE I.

Vertical ionization potential (eV) obtained at the G0W0@PBE level of theory with different numbers of standard and kinetic PDEPs. (A) 20 stdPDEPs + up to 400 kinPDEPs and extrapolated, (B) 100 stdPDEPs + up to 400 kinPDEPs and extrapolated, and (C) pure stdPDEPs and extrapolated. A detailed discussion of extrapolations of quasiparticle energies can be found in Ref. 9. See also Fig. 3 of the supplementary material.

MoleculeABC
C2H2 11.07 11.06 11.06 
C2H4 10.41 10.40 10.40 
C4H48.80 8.77 8.76 
C6H6 9.17 9.14 9.13 
CH3Cl 11.28 11.26 11.25 
CH3OH 10.58 10.56 10.56 
CH3SH 9.39 9.36 9.36 
CH4 14.01 14.01 14.01 
Cl2 11.51 11.51 11.50 
ClF 12.55 12.55 12.54 
CO 13.51 13.50 13.50 
CO2 13.32 13.31 13.31 
CS 11.00 10.98 10.98 
F2 14.99 14.97 14.97 
H2CO 10.43 10.42 10.42 
H211.82 11.82 11.81 
H2O2 10.87 10.87 10.86 
HCl 12.50 12.50 12.50 
HCN 13.20 13.20 13.20 
Na2 4.95 4.95 4.95 
MoleculeABC
C2H2 11.07 11.06 11.06 
C2H4 10.41 10.40 10.40 
C4H48.80 8.77 8.76 
C6H6 9.17 9.14 9.13 
CH3Cl 11.28 11.26 11.25 
CH3OH 10.58 10.56 10.56 
CH3SH 9.39 9.36 9.36 
CH4 14.01 14.01 14.01 
Cl2 11.51 11.51 11.50 
ClF 12.55 12.55 12.54 
CO 13.51 13.50 13.50 
CO2 13.32 13.31 13.31 
CS 11.00 10.98 10.98 
F2 14.99 14.97 14.97 
H2CO 10.43 10.42 10.42 
H211.82 11.82 11.81 
H2O2 10.87 10.87 10.86 
HCl 12.50 12.50 12.50 
HCN 13.20 13.20 13.20 
Na2 4.95 4.95 4.95 
TABLE II.

Vertical electron affinity (eV) obtained at the G0W0@PBE level of theory with different numbers of standard and kinetic PDEPs. (A) 20 stdPDEPs + up to 400 kinPDEPs and extrapolated, (B) 100 stdPDEPs + up to 400 kinPDEPs and extrapolated, and (C) pure stdPDEPs and extrapolated. We defined the electron affinity as the energy of the first unoccupied state, within the given unit cell used, extrapolated as a function of NstdPDEP + NkinPDEP. See also Ref. 9 and Fig. 3 of the supplementary material.

MoleculeABC
C2H2 −2.42 −2.41 −2.41 
C2H4 −1.75 −1.75 −1.75 
C4H4−0.85 −0.81 −0.80 
C6H6 −1.01 −0.96 −0.96 
CH3Cl −1.17 −1.16 −1.16 
CH3OH −0.89 −0.89 −0.89 
CH3SH −0.88 −0.88 −0.88 
CH4 −0.64 −0.64 −0.64 
Cl2 1.65 1.64 1.65 
ClF 1.28 1.28 1.28 
CO −1.56 −1.57 −1.57 
CO2 −0.97 −0.97 −0.97 
CS 0.49 0.51 0.51 
F2 1.16 1.16 1.16 
H2CO −0.69 −0.68 −0.68 
H2−0.90 −0.90 −0.90 
H2O2 −1.80 −1.79 −1.79 
HCl −1.07 −1.07 −1.07 
HCN −2.08 −2.08 −2.08 
Na2 0.64 0.63 0.63 
MoleculeABC
C2H2 −2.42 −2.41 −2.41 
C2H4 −1.75 −1.75 −1.75 
C4H4−0.85 −0.81 −0.80 
C6H6 −1.01 −0.96 −0.96 
CH3Cl −1.17 −1.16 −1.16 
CH3OH −0.89 −0.89 −0.89 
CH3SH −0.88 −0.88 −0.88 
CH4 −0.64 −0.64 −0.64 
Cl2 1.65 1.64 1.65 
ClF 1.28 1.28 1.28 
CO −1.56 −1.57 −1.57 
CO2 −0.97 −0.97 −0.97 
CS 0.49 0.51 0.51 
F2 1.16 1.16 1.16 
H2CO −0.69 −0.68 −0.68 
H2−0.90 −0.90 −0.90 
H2O2 −1.80 −1.79 −1.79 
HCl −1.07 −1.07 −1.07 
HCN −2.08 −2.08 −2.08 
Na2 0.64 0.63 0.63 

Table III shows our results for the C60 molecule. The structure of C60 (point group Ih) 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 NstdPDEP = 2000). The two sets of calculations starting with NstdPDEP = 100 and NkinPDEP = 200 amounted to savings of 32% and 15%, respectively.

TABLE III.

Quasiparticle energies (eV) of C60 calculated at the G0W0@PBE level of theory. Energy levels are labeled by their symmetry in point group Ih. NstdPDEP = 100 and NstdPDEP = 200 are calculations with 100 and 200 stdPDEPs and up to 400 kinPDEPs and extrapolated. NkinPDEP = 0 is the calculation with pure stdPDEPs and extrapolated (see the text).

Energy levelsNstdPDEP = 100NstdPDEP = 200NkinPDEP = 0G0W0Expt.
t1g −1.70 −1.66 −1.70   
t1u −2.70 −2.77 −2.82 −2.74,a −2.62,b −2.82c −2.69d 
hu −7.32 −7.32 −7.38 −7.31,a −7.21,b −7.37c −7.61,d −7.6e 
gg + hg −8.46, −8.52 −8.46, −8.51 −8.51, −8.56 −8.68,c −8.69c −8.59e 
Energy levelsNstdPDEP = 100NstdPDEP = 200NkinPDEP = 0G0W0Expt.
t1g −1.70 −1.66 −1.70   
t1u −2.70 −2.77 −2.82 −2.74,a −2.62,b −2.82c −2.69d 
hu −7.32 −7.32 −7.38 −7.31,a −7.21,b −7.37c −7.61,d −7.6e 
gg + hg −8.46, −8.52 −8.46, −8.51 −8.51, −8.56 −8.68,c −8.69c −8.59e 
a

Reference 27: with 700 standard eigenpotentials.

b

Ref. 20: with 27387 SAPOs, where SAPO stands for simple approximate physical orbitals.

c

Ref. 45: plane wave energy cutoff of 45 Ry and cell edge of 31.7 bohrs.

d

Reference 46.

e

Reference 47.

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

FIG. 4.

Ball and stick representation of the atomistic structure48 of the Si3N4/Si(100) interface used in our study.

FIG. 4.

Ball and stick representation of the atomistic structure48 of the Si3N4/Si(100) interface used in our study.

Close modal

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 G0W0@PBE calculations for each bulk system separately and obtained quasiparticle energies.

The local density of states is given by

D(ε,z)=2idxLxdyLy|ψi(x,y,z)|2δ(εεi),
(6)

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 interface48 

VBMEFD(ε,z)dε=EFCBMD(ε,z)dε=ΔEFD(ε,z)dε,
(7)

where EF 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.

TABLE IV.

Bandgaps of bulk Si, a −Si3N4, and band offsets (VBO and CBO) of the interface (see Fig. 4). All values are in eV.

Method∖EnergyVBOCBOEgSiEgSi3N4
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 
G0W0 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.78b 1.82–2.83c 1.17d 4.5–5.5e 
Method∖EnergyVBOCBOEgSiEgSi3N4
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 
G0W0 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.78b 1.82–2.83c 1.17d 4.5–5.5e 
a

Reference 48.

b

References 52–54.

c

Estimated by the other three experimental values.

d

Reference 55.

e

References 56–58.

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

V¯(r)=VH(r)+Vloc(r)iV¯at(i)(|rri|),
(8)

where V¯at(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 G0W0 corrections on band offsets, we performed G0W0@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 G0W0 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 G0W0 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%.

TABLE V.

Quasiparticle energies of valence band maximum (VBM) and conduction band minimum (CBM) of bulk silicon and amorphous Si3N4 computed with standard eigenpotentials and by combining standard and kinetic eigenpotentials. Columns NstdPDEP = 1000 and NstdPDEP = 2000 are calculations done with 1000 and 2000 stdPDEPs; column fit is extrapolated results; column NkinPDEP = 400 is the calculation with up to 400 kinPDEPs and extrapolated (see the text).

NstdPDEP = 1000NstdPDEP = 2000FitNkinPDEP = 400
Si VBM 5.70 5.55 5.45 5.53 
Si CBM 7.03 6.91 6.79 6.82 
a −Si3N4 VBM 7.14 7.01 7.01 6.99 
a −Si3N4 CBM 11.99 11.87 11.83 11.83 
NstdPDEP = 1000NstdPDEP = 2000FitNkinPDEP = 400
Si VBM 5.70 5.55 5.45 5.53 
Si CBM 7.03 6.91 6.79 6.82 
a −Si3N4 VBM 7.14 7.01 7.01 6.99 
a −Si3N4 CBM 11.99 11.87 11.83 11.83 

The method introduced in Refs. 9 and 26–28 to compute quasiparticle energies using the G0W0 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 G0W0 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 G0W0 calculations by 10%–50%, depending on the system, with the most savings observed for the largest systems studied here.

See the supplementary material for convergence studies of the G0W0 calculations of the vertical ionization potential of the CH4 molecule, which is taken as a representative example of the molecular systems studied in the main text.

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.

1.
P.
Hohenberg
and
W.
Kohn
, “
Inhomogeneous electron gas
,”
Phys. Rev.
136
,
B864
(
1964
).
2.
W.
Kohn
and
L. J.
Sham
, “
Self-consistent equations including exchange and correlation effects
,”
Phys. Rev.
140
,
A1133
(
1965
).
3.
R. M.
Martin
,
Electronic Structure: Basic Theory and Practical Methods
(
Cambridge University Press
,
2004
).
4.
G.
Onida
,
L.
Reining
, and
A.
Rubio
, “
Electronic excitations: Density-functional versus many-body green’s-function approaches
,”
Rev. Mod. Phys.
74
,
601
659
(
2002
).
5.
R. M.
Martin
,
L.
Reining
, and
D. M.
Ceperley
,
Interacting Electrons
(
Cambridge University Press
,
2016
).
6.
L.
Hedin
, “
New method for calculating the one-particle green’s function with application to the electron-gas problem
,”
Phys. Rev.
139
,
A796
(
1965
).
7.
M. J.
van Setten
,
F.
Caruso
,
S.
Sharifzadeh
,
X.
Ren
,
M.
Scheffler
,
F.
Liu
,
J.
Lischner
,
L.
Lin
,
J. R.
Deslippe
,
S. G.
Louie
 et al, “
GW 100: Benchmarking G0W0 for molecular systems
,”
J. Chem. Theory Comput.
11
,
5665
5687
(
2015
).
8.
E.
Maggio
,
P.
Liu
,
M. J.
van Setten
, and
G.
Kresse
, “
GW100: A plane wave perspective for small molecules
,”
J. Chem. Theory Comput.
13
,
635
648
(
2017
).
9.
M.
Govoni
and
G.
Galli
, “
GW100: Comparison of methods and accuracy of results obtained with the west code
,”
J. Chem. Theory Comput.
14
,
1895
1909
(
2018
).
10.
P.
Scherpelz
,
M.
Govoni
,
I.
Hamada
, and
G.
Galli
, “
Implementation and validation of fully relativistic GW calculations: Spin–orbit coupling in molecules, nanocrystals, and solids
,”
J. Chem. Theory Comput.
12
,
3523
3544
(
2016
).
11.
N. P.
Brawand
,
M.
Vörös
,
M.
Govoni
, and
G.
Galli
, “
Generalization of dielectric-dependent hybrid functionals to finite systems
,”
Phys. Rev. X
6
,
041002
(
2016
).
12.
D.
Golze
,
M.
Dvorak
, and
P.
Rinke
, “
The GW compendium: A practical guide to theoretical photoemission spectroscopy
,”
Front. Chem.
7
,
377
(
2019
).
13.
M. S.
Hybertsen
and
S. G.
Louie
, “
First-principles theory of quasiparticles: Calculation of band gaps in semiconductors and insulators
,”
Phys. Rev. Lett.
55
,
1418
(
1985
).
14.
M. S.
Hybertsen
and
S. G.
Louie
, “
Electron correlation and the band gap in ionic crystals
,”
Phys. Rev. B
32
,
7005
7008
(
1985
).
15.
M. S.
Hybertsen
and
S. G.
Louie
, “
Electron correlation in semiconductors and insulators: Band gaps and quasiparticle energies
,”
Phys. Rev. B
34
,
5390
(
1986
).
16.
G.
Strinati
,
H.
Mattausch
, and
W.
Hanke
, “
Dynamical correlation effects on the quasiparticle Bloch states of a covalent crystal
,”
Phys. Rev. Lett.
45
,
290
(
1980
).
17.
G.
Strinati
,
H.
Mattausch
, and
W.
Hanke
, “
Dynamical aspects of correlation corrections in a covalent crystal
,”
Phys. Rev. B
25
,
2867
(
1982
).
18.
S. L.
Adler
, “
Quantum theory of the dielectric constant in real solids
,”
Phys. Rev.
126
,
413
(
1962
).
19.
N.
Wiser
, “
Dielectric constant with local field effects included
,”
Phys. Rev.
129
,
62
(
1963
).
20.
G.
Samsonidze
,
M.
Jain
,
J.
Deslippe
,
M. L.
Cohen
, and
S. G.
Louie
, “
Simple approximate physical orbitals for GW quasiparticle calculations
,”
Phys. Rev. Lett.
107
,
186404
(
2011
).
21.
F.
Bruneval
and
X.
Gonze
, “
Accurate GW self-energies in a plane-wave basis using only a few empty states: Towards large systems
,”
Phys. Rev. B
78
,
085125
(
2008
).
22.
W.
Gao
,
W.
Xia
,
X.
Gao
, and
P.
Zhang
, “
Speeding up GW calculations to meet the challenge of large scale quasiparticle predictions
,”
Sci. Rep.
6
,
36849
(
2016
).
23.
J.
Soininen
,
J.
Rehr
, and
E. L.
Shirley
, “
Electron self-energy calculation using a general multi-pole approximation
,”
J. Phys.: Condens. Matter
15
,
2573
(
2003
).
24.
H. F.
Wilson
,
F.
Gygi
, and
G.
Galli
, “
Efficient iterative method for calculations of dielectric matrices
,”
Phys. Rev. B
78
,
113303
(
2008
).
25.
H. F.
Wilson
,
D.
Lu
,
F.
Gygi
, and
G.
Galli
, “
Iterative calculations of dielectric eigenvalue spectra
,”
Phys. Rev. B
79
,
245106
(
2009
).
26.
H.-V.
Nguyen
,
T. A.
Pham
,
D.
Rocca
, and
G.
Galli
, “
Improving accuracy and efficiency of calculations of photoemission spectra within the many-body perturbation theory
,”
Phys. Rev. B
85
,
081101
(
2012
).
27.
T. A.
Pham
,
H.-V.
Nguyen
,
D.
Rocca
, and
G.
Galli
, “
GW calculations using the spectral decomposition of the dielectric matrix: Verification, validation, and comparison of methods
,”
Phys. Rev. B
87
,
155148
(
2013
).
28.
M.
Govoni
and
G.
Galli
, “
Large scale GW calculations
,”
J. Chem. Theory Comput.
11
,
2680
2696
(
2015
).
29.
R.
Sternheimer
, “
Electronic polarizabilities of ions from the Hartree-Fock wave functions
,”
Phys. Rev.
96
,
951
(
1954
).
30.
G.
Galli
and
A.
Pasquarello
, “
First-principles molecular dynamics
,” in
Computer Simulation in Chemical Physics
, edited by
M. P.
Allen
and
D. J.
Tildesley
(
Springer
,
The Netherlands
,
1993
), pp.
261
313
.
31.
S.
Baroni
,
S.
De Gironcoli
,
A.
Dal Corso
, and
P.
Giannozzi
, “
Phonons and related crystal properties from density-functional perturbation theory
,”
Rev. Mod. Phys.
73
,
515
(
2001
).
32.
H.
Ma
,
M.
Govoni
,
F.
Gygi
, and
G.
Galli
, “
A finite-field approach for GW calculations beyond the random phase approximation
,”
J. Chem. Theory Comput.
15
,
154
164
(
2018
).
33.
D.
Lu
,
F.
Gygi
, and
G.
Galli
, “
Dielectric properties of ice and liquid water from first-principles calculations
,”
Phys. Rev. Lett.
100
,
147601
(
2008
).
34.
J.
Lindhard
, “
On the properties of a gas of charged particles
,”
Dan. Vid. Selsk Mat.-Fys. Medd.
28
,
8
(
1954
).
35.
D.
Rocca
, “
Random-phase approximation correlation energies from lanczos chains and an optimal basis set: Theory and applications to the benzene dimer
,”
J. Chem. Phys.
140
,
18A501
(
2014
).
36.
D.
Lu
,
Y.
Li
,
D.
Rocca
, and
G.
Galli
, “
Ab initio calculation of van der Waals bonded molecular crystals
,”
Phys. Rev. Lett.
102
,
206411
(
2009
).
37.
P.
Giannozzi
,
S.
Baroni
,
N.
Bonini
,
M.
Calandra
,
R.
Car
,
C.
Cavazzoni
,
D.
Ceresoli
,
G. L.
Chiarotti
,
M.
Cococcioni
,
I.
Dabo
 et al, “
QUANTUM ESPRESSO: A modular and open-source software project for quantum simulations of materials
,”
J. Phys.: Condens. Matter
21
,
395502
(
2009
).
38.
P.
Giannozzi
,
O.
Andreussi
,
T.
Brumme
,
O.
Bunau
,
M. B.
Nardelli
,
M.
Calandra
,
R.
Car
,
C.
Cavazzoni
,
D.
Ceresoli
,
M.
Cococcioni
,
N.
Colonna
,
I.
Carnimeo
,
A. D.
Corso
,
S.
de Gironcoli
,
P.
Delugas
,
R. A.
DiStasio
,
A.
Ferretti
,
A.
Floris
,
G.
Fratesi
,
G.
Fugallo
,
R.
Gebauer
,
U.
Gerstmann
,
F.
Giustino
,
T.
Gorni
,
J.
Jia
,
M.
Kawamura
,
H.-Y.
Ko
,
A.
Kokalj
,
E.
Küçükbenli
,
M.
Lazzeri
,
M.
Marsili
,
N.
Marzari
,
F.
Mauri
,
N. L.
Nguyen
,
H.-V.
Nguyen
,
A. O.
de-la Roza
,
L.
Paulatto
,
S.
Poncé
,
D.
Rocca
,
R.
Sabatini
,
B.
Santra
,
M.
Schlipf
,
A. P.
Seitsonen
,
A.
Smogunov
,
I.
Timrov
,
T.
Thonhauser
,
P.
Umari
,
N.
Vast
,
X.
Wu
, and
S.
Baroni
, “
Advanced capabilities for materials modelling with quantum ESPRESSO
,”
J. Phys.: Condens. Matter
29
,
465901
(
2017
).
39.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
, “
Generalized gradient approximation made simple
,”
Phys. Rev. Lett.
77
,
3865
3868
(
1996
).
40.
M.
Schlipf
and
F.
Gygi
, “
Optimization algorithm for the generation of ONCV pseudopotentials
,”
Comput. Phys. Commun.
196
,
36
44
(
2015
).
41.
D.
Hamann
, “
Optimized norm-conserving Vanderbilt pseudopotentials
,”
Phys. Rev. B
88
,
085117
(
2013
).
42.
L. A.
Curtiss
,
P. C.
Redfern
,
K.
Raghavachari
, and
J. A.
Pople
, “
Assessment of Gaussian-2 and density functional theories for the computation of ionization potentials and electron affinities
,”
J. Chem. Phys.
109
,
42
55
(
1998
).
43.

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.

44.
NIST Computational Chemistry Comparison and Benchmark Databasxe, Russell D. Johnson III NIST Standard Reference Database Number 101, Release 18,
2016
.
45.
X.
Qian
,
P.
Umari
, and
N.
Marzari
, “
First-principles investigation of organic photovoltaic materials C60, C70,[C60]PCBM, and bis-[C60]PCBM using a many-body G0W0-Lanczos approach
,”
Phys. Rev. B
91
,
245105
(
2015
).
46.
E. P.
Linstrom
and
W.
Mallard
,
NIST Chemistry WebBook: NIST Standard Reference Database Number 69
(
National Institute of Standards and Technology
,
Gaithersburg, MD
,
2018
).
47.
D. L.
Lichtenberger
,
K. W.
Nebesny
,
C. D.
Ray
,
D. R.
Huffman
, and
L. D.
Lamb
, “
Valence and core photoelectron spectroscopy of C60, buckminsterfullerene
,”
Chem. Phys. Lett.
176
,
203
208
(
1991
).
48.
T.
Anh Pham
,
T.
Li
,
H.-V.
Nguyen
,
S.
Shankar
,
F.
Gygi
, and
G.
Galli
, “
Band offsets and dielectric properties of the amorphous Si3N4/Si(100) interface: A first-principles study
,”
Appl. Phys. Lett.
102
,
241603
(
2013
).
49.
T.
Yamasaki
,
C.
Kaneta
,
T.
Uchiyama
,
T.
Uda
, and
K.
Terakura
, “
Geometric and electronic structures of SiO2/Si(001) interfaces
,”
Phys. Rev. B
63
,
115314
(
2001
).
50.
C. G.
Van de Walle
and
R. M.
Martin
, “
Theoretical study of band offsets at semiconductor interfaces
,”
Phys. Rev. B
35
,
8154
8165
(
1987
).
51.
R.
Sundararaman
and
Y.
Ping
, “
First-principles electrostatic potentials for reliable alignment at interfaces and defects
,”
J. Chem. Phys.
146
,
104109
(
2017
).
52.
J. W.
Keister
,
J. E.
Rowe
,
J. J.
Kolodziej
,
H.
Niimi
,
T. E.
Madey
, and
G.
Lucovsky
, “
Band offsets for ultrathin SiO2 and Si3N4 films on Si(111) and Si(100) from photoemission spectroscopy
,”
J. Vac. Sci. Technol., B: Nanotechnol. Microelectron.: Mater., Process., Meas., Phenom.
17
,
1831
1835
(
1999
).
53.
V. A.
Gritsenko
,
A. V.
Shaposhnikov
,
W.
Kwok
,
H.
Wong
, and
G. M.
Jidomirov
, “
Valence band offset at silicon/silicon nitride and silicon nitride/silicon oxide interfaces
,”
Thin Solid Films
437
,
135
139
(
2003
).
54.
M.
Higuchi
,
S.
Sugawa
,
E.
Ikenaga
,
J.
Ushio
,
H.
Nohira
,
T.
Maruizumi
,
A.
Teramoto
,
T.
Ohmi
, and
T.
Hattori
, “
Subnitride and valence band offset at Si3N4/Si interface formed using nitrogen-hydrogen radicals
,”
Appl. Phys. Lett.
90
,
123114
(
2007
).
55.
C.
Kittel
,
Introduction to Solid State Physics
(
Wiley
,
2005
).
56.
A. M.
Goodman
, “
Photoemission of electrons and holes into silicon nitride
,”
Appl. Phys. Lett.
13
,
275
277
(
1968
).
57.
J.
Bauer
, “
Optical properties, band gap, and surface roughness of Si3N4
,”
Phys. Status Solidi (A)
39
,
411
418
(
1977
).
58.
S. V.
Deshpande
,
E.
Gulari
,
S. W.
Brown
, and
S. C.
Rand
, “
Optical properties of silicon nitride films deposited by hot filament chemical vapor deposition
,”
J. Appl. Phys.
77
,
6534
6541
(
1995
).

Supplementary Material