A Finite-field Approach for GW Calculations Beyond the Random Phase Approximation

A Finite-field Approach for Calculations Beyond the Random Phase Approximation

He Ma Institute for Molecular Engineering, University of Chicago, Chicago, USA.    Marco Govoni Institute for Molecular Engineering, University of Chicago, Chicago, USA.    Francois Gygi Department of Computer Science, University of California Davis, Davis, USA.    Giulia Galli Institute for Molecular Engineering, University of Chicago, Chicago, USA. gagalli@uchicago.edu

We describe a finite-field approach to compute density response functions, which allows for efficient and calculations beyond the random phase approximation. The method is easily applicable to density functional calculations performed with hybrid functionals. We present results for the electronic properties of molecules and solids and we discuss a general scheme to overcome slow convergence of quasiparticle energies obtained from calculations, as a function of the basis set used to represent the dielectric matrix.


Department of Chemistry, University of Chicago, Chicago, USA. \alsoaffiliationMaterials Science Division, Argonne National Laboratory, Chicago, USA. \alsoaffiliationDepartment of Chemistry, University of Chicago, Chicago, USA. \alsoaffiliationMaterials Science Division, Argonne National Laboratory, Chicago, USA.

1 Introduction

Accurate, first principles predictions of the electronic structure of molecules and materials are important goals in chemistry, condensed matter physics and materials science 1. In the past three decades, density functional theory (DFT) 2, 3 has been successfully adopted to predict numerous properties of molecules and materials 4. In principle, any ground or excited state properties can be formulated as functionals of the ground state charge density. In practical calculations, the ground state charge density is determined by solving the Kohn-Sham (KS) equations with approximate exchange-correlation functionals, and many important excited state properties are not directly accessible from the solution of the KS equations. The time-dependent formulation of DFT (TDDFT) 5 in the frequency domain 6 provides a computationally tractable method to compute excitation energies and absorption spectra. However, using the common adiabatic approximation to the exchange-correlation functional, TDDFT is often not sufficiently accurate to describe certain types of excited states such as Rydberg and charge transfer states 7, especially when semi-local functionals are used.

A promising approach to predict excited state properties of molecules and materials is the many-body perturbation theory (MBPT) 8, 9, 10. Within MBPT, the approximation can be used to compute quasiparticle energies that correspond to photoemission and inverse photoemission measurements; furthermore, by solving the Bethe-Salpeter equation (BSE), one can obtain neutral excitation energies corresponding to optical spectra. For many years since the first applications of MBPT 9, its use has been hindered by its high computational cost. In the last decade, several advances have been proposed to improve the efficiency of MBPT calculations 11, 12, 13, which are now applicable to simulations of relatively large and complex systems, including nanostructures and heterogeneous interfaces 14, 15, 16. In particular, and BSE calculations can be performed using a low rank representation of density response functions 17, 18, 19, 20, whose spectral decomposition is obtained through iterative diagonalization using density functional perturbation theory (DFPT) 21, 22. This method does not require the explicit calculation of empty electronic states and avoids the inversion or storage of large dielectric matrices. The resulting implementation in the WEST code \bibnotewww.west-code.org has been successfully applied to investigate various systems including defects in semiconductors 24, 25, nanoparticles26, aqueous solutions27, 15, 28, and solid/liquid interfaces19 .

In this work, we developed a finite-field (FF) approach to evaluate density response functions, which enters the definition of the screened Coulomb interaction . The FF approach can be used as an alternative to DFPT, and presents the additional advantage of being applicable, in a straightforward manner, to both semilocal and hybrid functionals. In addition, FF calculations allow for the direct evaluation of density response functions beyond the random phase approximation (RPA).

Here we first benchmark the accuracy of the FF approach for the calculation of various density response functions, from which one can obtain the exchange correlation kernel (), defined as the functional derivative of the exchange-correlation potential with respect to the charge density. Then we discuss calculations for various molecules and solids, carried out with either semi-local or hybrid functionals, and by adopting different approximations to include vertex corrections in the self-energy. In the last two decades a variety of methods 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43 have been proposed to carry out vertex-corrected calculations, with different approximations to the vertex function and including various levels of self-consistency between , and . Here we focus on two formulations that are computationally tractable also for relatively large systems, denoted as and . In , is included in the evaluation of the screened Coulomb interaction ; in , is included in the calculation of both and the self-energy through the definition of a local vertex function. Most previous and calculations were restricted to the use of the LDA functional 29, 30, 34, 35, for which an analytical expression of is available. A work from Paier et al. 44 reported results for solids with HSE03 range-separated hybrid functional 45, with the exact exchange part of defined through the nanoquanta kernel 46, 47, 48, 33. In this work semilocal and hybrid functionals are treated in equal footing, and we present calculations using LDA 49, PBE 50 and PBE0 51 functional, as well as a dielectric-dependent hybrid (DDH) functional for solids 52.

A recent study of Thygesen and co-workers 53 reported basis set convergence issues when performing calculations, which could be overcome by applying a proper renormalization to the short-range component of 54, 55, 56. In our work we generalized the renormalization scheme of Thygesen et al. to functionals other than LDA, and we show that the convergence of quasiparticle energies is significantly improved using the renormalized .

The rest of the paper is organized as follows. In Sec. 2 we describe the finite-field approach and benchmark its accuracy. In Sec. 3 we describe the formalism used to perform calculations beyond the RPA, including a renormalization scheme for , and we compare the quasiparticle energies obtained from different approximations (RPA or vertex-corrected) for molecules in the GW100 test set 57 and for several solids. Finally, we summarize our results in Sec. 4.

2 The finite field approach

We first describe the FF approach for iterative diagonalization of density response functions and we then benchmark its robustness and accuracy.

2.1 Formalism

Our calculations are based on DFT single-particle energies and wavefunctions, obtained by solving the Kohn-Sham (KS) equations:


where the KS Hamiltonian . is the kinetic energy operator; is the KS potential that includes the ionic , the Hartree and the exchange-correlation potential . The charge density is given by . For simplicity we suppressed the spin index.

We consider the density response function (polarizability) of the KS system and that of the physical system ; the latter is denoted as when the random phase approximation (RPA) is used. The variation of the charge density due to either a variation of the KS potential or the external potential is given by:


where if and if . The density response functions of the KS and physical system are related by a Dyson-like equation:


where is the Coulomb kernel and is the exchange-correlation kernel.

Within the RPA, is neglected and is approximated by:


In the plane-wave representation (for simplicity we only focus on the point of the Brillouin zone), (abbreviated as ), and the dimensionless response function , obtained by symmetrizing with respect to , is defined as:


In the formalism which we adopted to perform calculations without empty states, one needs to compute a low rank decomposition of :


where and denote eigenvalue and eigenvectors of , respectively. The set of constitute the projective dielectric eigenpotential (PDEP) basis 17, 18, 19, and the accuracy of the low rank decomposition is controlled by , the size of the PDEP basis. and are obtained through iterative diagonalization, e.g. with the Davidson algorithm 58, which requires to evaluate the action of on an arbitrary function :


where and denote forward and inverse Fourier transforms, respectively.

Defining = , the calculations of the real space integral in Eq. 7 is equivalent to solving for the variation of the charge density due to :


In Ref. 19, is solved using DFPT for the case of . In this work we solved Eq. 8 by a finite-field approach. In particular, two SCF calculations are performed under the action of the potentials :


and is computed as:


where a central difference instead of forward/backward difference is performed in Eq. 10 to increase the numerical accuracy of the computed .

Self-consistent solutions of Eq. 8-10 correspond to considering in Eq. 12. If is kept fixed during the SCF iterations, the solution of Eq. 10 corresponds to . If both and are kept fixed during the SCF iterations, the solution of Eq. 10 corresponds to .

The response functions and (see Eq. 4) have the same eigenvectors, and their eigenvalues are related by:


where ’s and ’s are eigenvalues of and , respectively. In general the eiegenvalues and eigenvectors of are different from those of due to the presence of in Eq. 3.

In comparison to DFPT, the finite-field approach adopted here allows for the straightforward calculation of response functions beyond the RPA (i.e. for the calculation of instead of or ), and it can be readily applied to hybrid functionals for which analytical expressions of are not available. We note that finite-field calculations with hybrid functionals can easily benefit from any methodological development that reduces the computational complexity of evaluating exact exchange potentials 59, 60, 61.

Once the PDEP basis is obtained by iterative diagonalization of \bibnoteHere, we defined the PDEP basis to be the eigenvectors of . Alternatively, one may first iteratively diagonalize and define its eigenvectors as the PDEP basis. Then and can be evaluated in the space of the eigenvectors of . This choice is not further discussed in the paper; we only mention that some comparisons for the quasiparticle energies (at the level, see Sec. 3) of selected molecules obtained using either or eigenvectors as the PDEP basis are identical within 0.01 (0.005) eV for the HOMO (LUMO) state., the projection of on the PDEP basis can be performed using the finite field approach as well. Then the symmetrized exchange-correlation kernel can be computed by inverting the Dyson-like equation (Eq. 3):


On the right hand side of Eq. 12 all matrices are and therefore the resulting is also defined on the PDEP basis.

When using orbital-dependent functionals such as meta-GGA and hybrid functionals, the computed from Eq. 12 needs to be interpreted with caution. In this case, DFT calculations for can be performed using either the optimized effective potential (OEP) or the generalized Kohn-Sham (GKS) scheme. In the OEP scheme, is local in space and depends on and , as in the case of semi-local functionals. In the GKS scheme, is non-local and depends on three position vectors. We expect to be almost independent of the chosen scheme, whether GKS or OEP, since both methods yield the same result within first order in the charge density 63. We conducted hybrid functional calculations within the GKS scheme, assuming that for every GKS calculation an OEP can be defined yielding the same charge density; with this assumption the from Eq. 12 is well defined within the OEP formalism.

2.2 Implementation and Verification

We implemented the finite-field algorithm described above by coupling the WEST 19 and Qbox 64 codes in client-server mode, using the workflow summarized in Fig 1. In particular, in our implementation the WEST code performs an iterative diagonalization of by outsourcing the evaluation of the action of on an arbitrary function to Qbox, which performs DFT calculations in finite field. The two codes communicate through the filesystem.

Figure 1: Workflow of finite-field calculations. The WEST code performs an iterative diagonalization of (, , ). In calculations beyond the RPA, is computed from Eq. 12, which requires computing the spectral decomposition of and evaluating in the space of eigenvectors. Finite-field calculations are carried out by the Qbox code, and the communications of and between WEST and Qbox is carried through the filesystem.

To verify the correctness of our implementation, we computed , , for selected molecules in the GW100 set and we compared the results to those obtained with DFPT. Sec. 1 of the SI summarizes the parameters used (, , etc.). In finite-field calculations we optimized the ground state wavefunction using a preconditioned steepest descent algorithm with Anderson acceleration65. The magnitude of was chosen to insure that calculations were performed within the linear response regime (see Sec. 2 of the SI). All calculations presented in this section were performed with the PBE functional unless otherwise specified.

Fig 2a shows the eigenvalues of for a few molecules obtained with three approaches: iterative diagonalization of with the finite-field approach; iterative diagonalization of with either the finite-field approach or with DFPT, followed by a transformation of eigenvalues as in Eq. 11. The three approaches yield almost identical eigenvalues.

Figure 2: Comparison of the eigenvalues(a) and eigenfunctions(b) of obtained from DFPT and finite-field (FF) calculations. Three approaches are used: diagonalization of by DFPT, diagonalization of by FF (denoted by FF(0)) and diagonalization of by FF (denoted by FF(RPA)). In the case of DFPT and FF(0), Eq. 11 was used to obtain the eigenvalues of from those of . In (b) we show the first elements of the and matrices (see Eq. 6).

The eigenvectors of the response functions are shown in Fig 2b, where we report elements of the matrices defined by the overlap between finite-field and DFPT eigenvectors. The inner product matrices are block-diagonal, with blocks corresponding to the presence of degenerate eigenvalues. The agreement between eigenvalues and eigenvectors shown in Fig 2 verifies the accuracy and robustness of finite-field calculations.

Fig 3 shows the eigendecomposition of compared to that of .

Figure 3: Comparison of eigenvalues(a) and eigenfunctions(b) of and obtained from finite-field calculations. In (b), the first elements of the matrices are presented.

As indicated by Fig 3a, including in results in a stronger screening. The eigenvalues of are systematically more negative than those of , though they asymptotically converge to zero in the same manner. While the eigenvalues are different, the eigenvectors (eigenspaces in the case of degenerate eigenvalues) are almost identical, as indicated by the block-diagonal form of the eigenvector overlap matrices (see Fig. 3b).

Finally, can be computed from and according to Eq. 12. Due to the similarity of the eigenvectors of and (identical to that of ), the matrix is almost diagonal. In Sec. 3 of the SI we show the matrix in the PDEP basis for a few systems. To verify the accuracy of obtained by the finite-field approach, we performed calculations with the LDA functional, for which can be computed analytically. In Fig 4 we present for a number of systems the average relative difference of the diagonal terms of the matrices obtained analytically and through finite-field (FF) calculations. We define as


As shown in Fig 4, is smaller than a few percent for all systems studied here. To further quantify the effect of the small difference found for the matrices on quasiparticle energies, we performed calculations ( is included in , see Sec. 3) for all the systems shown in Fig 4, using the analytical and computed from finite-field calculations. The two approaches yielded almost identical quasiparticle energies, with mean absolute deviations of 0.04 and 0.004 eV for HOMO and LUMO levels, respectively.

Figure 4: Average relative differences (see Eq. 13) between diagonal elements of the matrices computed analytically and numerically with the finite-field approach. Calculations were performed with the LDA functional.

3 calculations

3.1 Formalism

In this section we discuss calculations within and beyond the RPA, utilizing computed with the finite-field approach. In the following equations we use 1, 2, … as shorthand notations for , , … Indices with bars are integrated over. When no indices are shown, the equation is a matrix equation in reciprocal space or in the PDEP basis. The following discussion focuses on finite systems; for periodic systems a special treatment of the long-range limit of is required and relevant formulae are presented in Sec. 4 of the SI.

Based on a KS reference system, the Hedin equations 8 relate the exchange-correlation self-energy (abbreviated as ), Green’s function , the screened Coulomb interaction , the vertex and the irreducible polarizability :


We consider three different approximations: the first is the common formulation within the RPA, here denoted as , where and is given by:






The second approximation, denoted as , includes in the definition of . Specifically, is computed from and with Eq. 3:


and is used to construct screened the Coulomb interaction beyond the RPA:


The third approximation, denoted , includes in both and . In particular, an initial guess for is constructed from :


from which one can obtain a zeroth order vertex function by iterating Hedin’s equations once 29:


Then the self-energy is constructed using , and :


where we defined an effective screened Coulomb interaction\bibnoteOne may note that is not symmetric with respect to its two indices, and it can be symmetrized by using in Eq. 31. We found that the symmetrization has negligible effects on quasiparticle energies. We performed calculations for systems as shown in Fig 4 with either symmetrized or unsymmetrized , the mean absolute deviations for HOMO and LUMO quasiparticle energies are 0.006 eV and 0.001 eV respectively.


The symmetrized forms of the three different density response functions (reducible polarizabilities) defined in Eq. 21, 22, 28 are:


Eqs. 29-31 have been implemented in the WEST code 19.

We note that finite-field calculations yield matrices at zero frequency. Hence the results presented here correspond to calculations performed within the adiabatic approximation, as they neglect the frequency dependence of . An interesting future direction would be to compute frequency-dependent by performing finite-field calculations using real-time time-dependent DFT (RT-TDDFT).

When using the formalism, the convergence of quasiparticle energies with respect to turned out to be extremely challenging. As discussed in Ref. 53 the convergence problem originates from the incorrect short-range behavior of . In Sec. 3.2 below we describe a renormalization scheme of that improves the convergence of results.

3.2 Renormalization of

Thygesen and co-workers 53 showed that calculations with computed at the LDA level exhibit poor convergence with respect to the number of unoccupied states and plane wave cutoff. We observed related convergence problems of quasiparticle energies as a function of , the size of the basis set used here to represent response functions, see Sec. 5 of the SI. In this section we describe a generalization of the renormalization scheme proposed by Thygesen and co-workers 54, 55, 56 to overcome the convergence issues.

The approach of Ref. 53 is based on the properties of the homogeneous electron gas (HEG). For an HEG with density , depends only on due to translational invariance, and therefore is diagonal in reciprocal space. We denote the diagonal elements of as where . When using the LDA functional, the exchange kernel exactly cancels the Coulomb interaction at wavevector (the correlation kernel is small compared to for ), where is the Fermi wavevector. For , shows an incorrect asymptotic behavior, leading to an unphysical correlation hole 54, 55. Hence Thygesen and co-workers introduced a renormalized LDA kernel by setting for and for . They demonstrated that the renormalized improves the description of the short-range correlation hole as well as the correlation energy, and when applied to calculations substantially accelerates the basis set convergence of quasiparticle energies.

While within LDA can be computed analytically and at exactly , for a general functional it is not known a priori at which this condition is satisfied. In addition, for inhomogenous systems such as molecules and solids the matrix is not diagonal in reciprocal space. Ref. 53 used a wavevector symmetrization approach to evaluate for inhomogenous systems, which is not easily generalizable to the formalism in this work where is represented in the PDEP basis.

To overcome these difficulties, here we first diagonalize the matrix in the PDEP basis:


where and are eigenvalues and eigenvectors of . Then we define a renormalized as:


Note that for , , therefore is strictly greater or equal to . When applied to the HEG, the is equivalent to in the limit , where the PDEP and plane-wave basis are related by a unitary transformation. Thus, Eq. 33 represents a generalization of the scheme of Thygesen et al. to any functional and to inhomogeneous electron gases. When using , we observed a faster basis set convergence of results than results, consistent with Ref. 53. In Sec. 5 of the SI we discuss in detail the effect of the renormalization on the description of the density response functions and , and we rationalize why the renormalization improves the convergence of results. Here we only mention that the response function may possess positive eigenvalues for large PDEP indices. When the renormalized is used, the eigenvalues of are guaranteed to be nonpositive and decay rapidly toward zero as the PDEP index increase, which explains the improved convergence of quasiparticle energies.

All results shown in Sec. 3.3 were obtained with renormalized matrices, while calculations were conducted without renormalizing , since we found that the renormalization had a negligible effect on quasiparticle energies (see SI Sec. 5).

3.3 Results

In this section we report quasiparticle energies for molecules in the GW100 set 57 and for several solids. Calculations are performed at , and levels of theory and with semi-local and hybrid functionals. Computational parameters (, , etc.) for all calculations are summarized in Sec. 1 of the SI. A discussion of the convergence of quasiparticle energies with respect to these parameters can be found in Ref. 20.

We computed the vertical ionization energy (VIP), vertical electron affinity (VEA) and fundamental gaps for molecules with LDA, PBE and PBE0 functionals. VIP and VEA are defined as and respectively, where is the vacuum level estimated with the Makov-Payne method 67; and are HOMO and LUMO quasiparticle energies, respectively. The results are summarized in Fig 5 and compared with experimental values \bibnoteExperimental data are taken from the WEST GW100 data collection website (www.west-code.org/database). ; VIP values are also compared with existing results from quantum chemistry CCSD(T) calculations 69.

Figure 5: Vertical ionization energy (VIP), vertical electron affinity (VEA) and electronic gap of molecules in the GW100 set computed at , and levels of theory, compared to experimental and CCSD(T) results (black dashed lines) \bibnoteFor \ceKH molecule, calculation for the HOMO converged to a satellite instead of the quasiparticle peak. The spectral function of \ceKH is plotted and discussed in SI Sec. 6 and the correct quasiparticle energy is used here..

The VIP values computed at () level are systematically higher (lower) than the corresponding results. Compared to experiments and CCSD(T) results, results show an improvement over ones, when semilocal functionals (LDA, PBE) are used as starting points; instead, when using the PBE0 functional leads to a slight overestimation of VIP. results underestimate VIP with all functionals tested here, compared to experiments. At the LDA level, the comparison between different approximations is consistent with the observation by Morris et al for He, Be and Ne atoms 35. In general, vertex corrections have larger effects on quasiparticle energies computed with hybrid functionals, as the difference between the results obtained with different approximations are more prominent when starting from PBE0 calculations than starting from semi-local ones.

Different approximations yielded very similar results for the VEA, with () results marginally higher (lower) than ones.

Finally we report , and results for several solids: \ceSi, \ceSiC (4H), \ceC (diamond), \ceAlN, \ceWO3 (monoclinic), \ceSi3N4 (amorphous). We performed calculations starting with LDA and PBE functionals for all solids, and for Si we also performed calculations with the dielectric-dependent hybrid (DDH) functional 52. All solids are represented by supercells with 64-96 atoms (see Sec. 1 of the SI) and only the -point is used to sample the Brillioun zone. In Table 1 we present the band gaps computed with different approximations and functionals. Note that the supercells used here do not yield fully converged results as a function of supercell size (or k-point sampling); however the comparisons between different calculations are sound and represent the main result we are discussing in this section.

System XC
\ceSi LDA 0.55 1.35 1.33 1.24
PBE 0.73 1.39 1.37 1.28
DDH 1.19 1.57 1.50 1.48
\ceC (diamond) LDA 4.28 5.99 6.00 5.89
PBE 4.46 6.05 6.06 5.95
\ceSiC (4H) LDA 2.03 3.27 3.23 3.26
PBE 2.21 3.28 3.23 3.28
\ceAlN LDA 3.85 5.67 5.72 5.66
PBE 4.04 5.67 5.74 5.68
\ceWO3 (monoclinic) LDA 1.68 3.10 3.07 3.15
PBE 1.78 2.97 2.87 3.03
\ceSi3N4 (amorphous) LDA 3.04 4.84 4.92 4.81
PBE 3.19 4.86 4.96 4.83
Table 1: Band gaps (eV) for solids computed by different approximations and exchange-correlation (XC) functionals (see text). All calculations are performed at the -point of supercells with 64-96 atoms (see Sec. 1 of the SI for details).

Overall, band gaps obtained with different approximations are very similar, with differences much smaller than those observed for molecules. To further investigate the positions of the band edges obtained from different approximations, we plotted in Fig 6 the quasiparticle corrections to VBM and CBM, defined as where and are the quasiparticle energy and the Kohn-Sham eigenvalue of VBM/CBM, respectively.

Figure 6: quasiparticle corrections to the valance band maximum (VBM) and the conduction band minimum (CBM). Circles, squares and triangles are , and results respectively; red, blue, green markers correspond to calculations with LDA, PBE and DDH functionals.

Compared to , VBM and CBM computed at the level are slightly lower, while VBM and CBM computed at the level are significantly higher. For Si, obtained with LDA starting point are -0.75/0.06 (), -0.86/-0.08 (), -0.21/0.49 eV () respectively, showing a trend in agreement with the results reported by Del Sole et al (-0.36/0.27, -0.44/0.14, 0.01/0.67 eV) 29, but with an overall overestimate of the band gap due a lack of convergence in our Brillouin zone sampling. The deviation of band edges computed by different approximations is larger with the DDH functional compared to that of semi-local functionals. Overall the trends observed for solids are consistent with those found for molecules, except that for solids the shift of the CBM resembles those of the VBM when vertex corrections are included, while for molecules VEA is insensitive to vertex corrections.

4 Conclusions

In summary, we developed a finite-field approach to compute density response functions (, and ) for molecules and materials. The approach is non-perturbative and can be used in a straightforward manner with both semilocal and orbital-dependent functionals. Using this approach, we computed the exchange-correlation kernel and performed calculations using dielectric responses evaluated beyond the RPA.

We evaluated quasiparticle energies for molecules and solids and compared results obtained within and beyond the RPA, and using DFT calculations with semi-local and hybrid functionals as input. We found that the effect of vertex corrections on quasiparticle energies is more notable when using input wavefunctions and single-particle energies from hybrid functionals calculations. For the small molecules in the GW100 set, calculations yielded higher VIP compared to results, leading to a better agreement with experimental and high-level quantum chemistry results when using LDA and PBE starting points, and to a slight overestimate of VIP when using PBE0 as the starting point. calculations instead yielded a systematic underestimate of VIP of molecules. VEA of molecules were found to be insensitive to vertex corrections. In the case of solids, the energy of the VBM and CBM shifts in the same direction, relative to RPA results, when vertex corrections are included, and overall the band gaps were found to be rather insensitive to the choice of the approximation.

In addition, we reported a scheme to renormalize , which is built on previous work 53 using the LDA functional. The scheme is general and applicable to any exchange-correlation functional and to inhomogeneous systems including molecules and solids. Using the renormalized , the basis set convergence of results was significantly improved.

Overall, the method introduced in our work represents a substantial progress towards efficient computations of dielectric screening and large-scale calculations for molecules and materials beyond the random phase approximation.


We thank Ngoc Linh Nguyen for helpful 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 number DE-AC02-06CH11357. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a DOE Office of Science User Facility supported by the Office of Science of the US Department of Energy under Contract No. DE-AC02-05CH11231, 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.


The Supporting Information contains parameters used for calculations, convergence tests, detailed discussion of matrix and its renormalization, extension of beyond-RPA formalism to solids, and an analysis of the spectral function of \ceKH molecule.

Table of Contents:


  • Onida et al. 2002 Onida, G.; Reining, L.; Rubio, A. Electronic excitations: density-functional versus many-body Green’s-function approaches. Reviews of Modern Physics 2002, 74, 601–659.
  • Hohenberg and Kohn 1964 Hohenberg, P.; Kohn, W. Inhomogeneous Electron Gas. Physical Review 1964, 136, B864–B871.
  • Kohn and Sham 1965 Kohn, W.; Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. Physical Review 1965, 140, A1133–A1138.
  • Becke 2014 Becke, A. D. Perspective: Fifty years of density-functional theory in chemical physics. The Journal of Chemical Physics 2014, 140, 18A301.
  • Runge and Gross 1984 Runge, E.; Gross, E. K. U. Density-Functional Theory for Time-Dependent Systems. Physical Review Letters 1984, 52, 997–1000.
  • Casida 1995 Casida, M. Recent Advances in Density Functional Methods; WORLD SCIENTIFIC, 1995; pp 155–192.
  • Casida and Huix-Rotllant 2012 Casida, M.; Huix-Rotllant, M. Progress in Time-Dependent Density-Functional Theory. Annual Review of Physical Chemistry 2012, 63, 287–323.
  • Hedin 1965 Hedin, L. New Method for Calculating the One-Particle Green’s Function with Application to the Electron-Gas Problem. Physical Review 1965, 139, A796–A823.
  • Hybertsen and Louie 1986 Hybertsen, M. S.; Louie, S. G. Electron correlation in semiconductors and insulators: Band gaps and quasiparticle energies. Physical Review B 1986, 34, 5390–5413.
  • Martin et al. 2016 Martin, R. M.; Reining, L.; Ceperley, D. M. Interacting Electrons; Cambridge University Press, 2016.
  • Umari et al. 2009 Umari, P.; Stenuit, G.; Baroni, S. Optimal representation of the polarization propagator for large-scale GW calculations. Physical Review B 2009, 79.
  • Neuhauser et al. 2014 Neuhauser, D.; Gao, Y.; Arntsen, C.; Karshenas, C.; Rabani, E.; Baer, R. Breaking the Theoretical Scaling Limit for Predicting Quasiparticle Energies: The Stochastic GW Approach. Physical Review Letters 2014, 113.
  • Liu et al. 2016 Liu, P.; Kaltak, M.; Klimeš, J.; Kresse, G. Cubic scaling GW: Towards fast quasiparticle calculations. Physical Review B 2016, 94.
  • Ping et al. 2013 Ping, Y.; Rocca, D.; Galli, G. Electronic excitations in light absorbers for photoelectrochemical energy conversion: first principles calculations based on many body perturbation theory. Chemical Society Reviews 2013, 42, 2437.
  • Pham et al. 2017 Pham, T. A.; Govoni, M.; Seidel, R.; Bradforth, S. E.; Schwegler, E.; Galli, G. Electronic structure of aqueous solutions: Bridging the gap between theory and experiments. Science Advances 2017, 3, e1603210.
  • Leng et al. 2016 Leng, X.; Jin, F.; Wei, M.; Ma, Y. GW method and Bethe-Salpeter equation for calculating electronic excitations. Wiley Interdisciplinary Reviews: Computational Molecular Science 2016, 6, 532–550.
  • Nguyen et al. 2012 Nguyen, H.-V.; Pham, T. A.; Rocca, D.; Galli, G. Improving accuracy and efficiency of calculations of photoemission spectra within the many-body perturbation theory. Physical Review B 2012, 85.
  • Pham et al. 2013 Pham, T. A.; Nguyen, H.-V.; Rocca, D.; Galli, G. GW calculations using the spectral decomposition of the dielectric matrix: Verification, validation, and comparison of methods. Physical Review B 2013, 87.
  • Govoni and Galli 2015 Govoni, M.; Galli, G. Large Scale GW Calculations. Journal of Chemical Theory and Computation 2015, 11, 2680–2696.
  • Govoni and Galli 2018 Govoni, M.; Galli, G. GW100: Comparison of Methods and Accuracy of Results Obtained with the WEST Code. Journal of Chemical Theory and Computation 2018, 14, 1895–1909.
  • Baroni et al. 1987 Baroni, S.; Giannozzi, P.; Testa, A. Green’s-function approach to linear response in solids. Physical Review Letters 1987, 58, 1861–1864.
  • Baroni et al. 2001 Baroni, S.; de Gironcoli, S.; Corso, A. D.; Giannozzi, P. Phonons and related crystal properties from density-functional perturbation theory. Reviews of Modern Physics 2001, 73, 515–562.
  • 23 www.west-code.org.
  • Seo et al. 2016 Seo, H.; Govoni, M.; Galli, G. Design of defect spins in piezoelectric aluminum nitride for solid-state hybrid quantum technologies. Scientific Reports 2016, 6.
  • Seo et al. 2017 Seo, H.; Ma, H.; Govoni, M.; Galli, G. Designing defect-based qubit candidates in wide-gap binary semiconductors for solid-state quantum technologies. Physical Review Materials 2017, 1.
  • Scherpelz et al. 2016 Scherpelz, P.; Govoni, M.; Hamada, I.; Galli, G. Implementation and Validation of Fully Relativistic GW Calculations: Spin–Orbit Coupling in Molecules, Nanocrystals, and Solids. Journal of Chemical Theory and Computation 2016, 12, 3523–3544.
  • Gaiduk et al. 2016 Gaiduk, A. P.; Govoni, M.; Seidel, R.; Skone, J. H.; Winter, B.; Galli, G. Photoelectron Spectra of Aqueous Solutions from First Principles. Journal of the American Chemical Society 2016, 138, 6912–6915.
  • Gaiduk et al. 2018 Gaiduk, A. P.; Pham, T. A.; Govoni, M.; Paesani, F.; Galli, G. Electron affinity of liquid water. Nature Communications 2018, 9.
  • Sole et al. 1994 Sole, R. D.; Reining, L.; Godby, R. W. GW approximation for electron self-energies in semiconductors and insulators. Physical Review B 1994, 49, 8024–8028.
  • Fleszar and Hanke 1997 Fleszar, A.; Hanke, W. Spectral properties of quasiparticles in a semiconductor. Physical Review B 1997, 56, 10228–10232.
  • Schindlmayr and Godby 1998 Schindlmayr, A.; Godby, R. W. Systematic Vertex Corrections through Iterative Solution of Hedin’s Equations Beyond the GW Approximation. Physical Review Letters 1998, 80, 1702–1705.
  • Marini and Rubio 2004 Marini, A.; Rubio, A. Electron linewidths of wide-gap insulators: Excitonic effects in LiF. Physical Review B 2004, 70.
  • Bruneval et al. 2005 Bruneval, F.; Sottile, F.; Olevano, V.; Sole, R. D.; Reining, L. Many-Body Perturbation Theory Using the Density-Functional Concept: Beyond the GW Approximation. Physical Review Letters 2005, 94.
  • Tiago and Chelikowsky 2006 Tiago, M. L.; Chelikowsky, J. R. Optical excitations in organic molecules, clusters, and defects studied by first-principles Green’s function methods. Physical Review B 2006, 73.
  • Morris et al. 2007 Morris, A. J.; Stankovski, M.; Delaney, K. T.; Rinke, P.; García-González, P.; Godby, R. W. Vertex corrections in localized and extended systems. Physical Review B 2007, 76.
  • Shishkin et al. 2007 Shishkin, M.; Marsman, M.; Kresse, G. Accurate Quasiparticle Spectra from Self-Consistent GW Calculations with Vertex Corrections. Physical Review Letters 2007, 99.
  • Shaltaf et al. 2008 Shaltaf, R.; Rignanese, G.-M.; Gonze, X.; Giustino, F.; Pasquarello, A. Band Offsets at the Si/SiO2 Interface from Many-Body Perturbation Theory. Physical Review Letters 2008, 100.
  • Romaniello et al. 2009 Romaniello, P.; Guyot, S.; Reining, L. The self-energy beyond GW: Local and nonlocal vertex corrections. The Journal of Chemical Physics 2009, 131, 154111.
  • Gruneis et al. 2014 Gruneis, A.; Kresse, G.; Hinuma, Y.; Oba, F. Ionization Potentials of Solids: The Importance of Vertex Corrections. Physical Review Letters 2014, 112.
  • Chen and Pasquarello 2015 Chen, W.; Pasquarello, A. Accurate band gaps of extended systems via efficient vertex corrections in GW. Physical Review B 2015, 92.
  • Kutepov 2016 Kutepov, A. L. Electronic structure of Na, K, Si, and LiF from self-consistent solution of Hedin’s equations including vertex corrections. Physical Review B 2016, 94.
  • Kutepov 2017 Kutepov, A. L. Self-consistent solution of Hedin’s equations: Semiconductors and insulators. Physical Review B 2017, 95.
  • Maggio and Kresse 2017 Maggio, E.; Kresse, G. GW Vertex Corrected Calculations for Molecular Systems. Journal of Chemical Theory and Computation 2017, 13, 4765–4778.
  • Paier et al. 2008 Paier, J.; Marsman, M.; Kresse, G. Dielectric properties and excitons for extended systems from hybrid functionals. Physical Review B 2008, 78.
  • Heyd et al. 2003 Heyd, J.; Scuseria, G. E.; Ernzerhof, M. Hybrid functionals based on a screened Coulomb potential. The Journal of Chemical Physics 2003, 118, 8207–8215.
  • Reining et al. 2002 Reining, L.; Olevano, V.; Rubio, A.; Onida, G. Excitonic Effects in Solids Described by Time-Dependent Density-Functional Theory. Physical Review Letters 2002, 88.
  • Marini et al. 2003 Marini, A.; Sole, R. D.; Rubio, A. Bound Excitons in Time-Dependent Density-Functional Theory: Optical and Energy-Loss Spectra. Physical Review Letters 2003, 91.
  • Sottile et al. 2003 Sottile, F.; Olevano, V.; Reining, L. Parameter-Free Calculation of Response Functions in Time-Dependent Density-Functional Theory. Physical Review Letters 2003, 91.
  • Perdew and Zunger 1981 Perdew, J. P.; Zunger, A. Self-interaction correction to density-functional approximations for many-electron systems. Physical Review B 1981, 23, 5048–5079.
  • Perdew et al. 1996 Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Physical Review Letters 1996, 77, 3865–3868.
  • Perdew et al. 1996 Perdew, J. P.; Ernzerhof, M.; Burke, K. Rationale for mixing exact exchange with density functional approximations. The Journal of Chemical Physics 1996, 105, 9982–9985.
  • Skone et al. 2014 Skone, J. H.; Govoni, M.; Galli, G. Self-consistent hybrid functional for condensed systems. Physical Review B 2014, 89.
  • Schmidt et al. 2017 Schmidt, P. S.; Patrick, C. E.; Thygesen, K. S. Simple vertex correction improves GW band energies of bulk and two-dimensional crystals. Physical Review B 2017, 96.
  • Olsen and Thygesen 2012 Olsen, T.; Thygesen, K. S. Extending the random-phase approximation for electronic correlation energies: The renormalized adiabatic local density approximation. Physical Review B 2012, 86.
  • Olsen and Thygesen 2013 Olsen, T.; Thygesen, K. S. Beyond the random phase approximation: Improved description of short-range correlation by a renormalized adiabatic local density approximation. Physical Review B 2013, 88.
  • Patrick and Thygesen 2015 Patrick, C. E.; Thygesen, K. S. Adiabatic-connection fluctuation-dissipation DFT for the structural properties of solids—The renormalized ALDA and electron gas kernels. The Journal of Chemical Physics 2015, 143, 102802.
  • van Setten et al. 2015 van Setten, M. J.; Caruso, F.; Sharifzadeh, S.; Ren, X.; Scheffler, M.; Liu, F.; Lischner, J.; Lin, L.; Deslippe, J. R.; Louie, S. G.; Yang, C.; Weigend, F.; Neaton, J. B.; Evers, F.; Rinke, P. GW100: Benchmarking G0W0 for Molecular Systems. Journal of Chemical Theory and Computation 2015, 11, 5665–5687.
  • Davidson 1975 Davidson, E. R. The iterative calculation of a few of the lowest eigenvalues and corresponding eigenvectors of large real-symmetric matrices. Journal of Computational Physics 1975, 17, 87–94.
  • Gygi 2009 Gygi, F. Compact Representations of Kohn-Sham Invariant Subspaces. Physical Review Letters 2009, 102.
  • Gygi and Duchemin 2012 Gygi, F.; Duchemin, I. Efficient Computation of Hartree-Fock Exchange Using Recursive Subspace Bisection. Journal of Chemical Theory and Computation 2012, 9, 582–587.
  • Dawson and Gygi 2015 Dawson, W.; Gygi, F. Performance and Accuracy of Recursive Subspace Bisection for Hybrid DFT Calculations in Inhomogeneous Systems. Journal of Chemical Theory and Computation 2015, 11, 4655–4663.
  • 62 Here, we defined the PDEP basis to be the eigenvectors of . Alternatively, one may first iteratively diagonalize and define its eigenvectors as the PDEP basis. Then and can be evaluated in the space of the eigenvectors of . This choice is not further discussed in the paper; we only mention that some comparisons for the quasiparticle energies (at the level, see Sec. 3) of selected molecules obtained using either or eigenvectors as the PDEP basis are identical within 0.01 (0.005) eV for the HOMO (LUMO) state.
  • Kümmel and Kronik 2008 Kümmel, S.; Kronik, L. Orbital-dependent density functionals: Theory and applications. Reviews of Modern Physics 2008, 80, 3–60.
  • Gygi 2008 Gygi, F. Architecture of Qbox: A scalable first-principles molecular dynamics code. IBM Journal of Research and Development 2008, 52, 137–144.
  • Anderson 1965 Anderson, D. G. Iterative Procedures for Nonlinear Integral Equations. Journal of the ACM 1965, 12, 547–560.
  • 66 One may note that is not symmetric with respect to its two indices, and it can be symmetrized by using in Eq. 31. We found that the symmetrization has negligible effects on quasiparticle energies. We performed calculations for systems as shown in Fig 4 with either symmetrized or unsymmetrized , the mean absolute deviations for HOMO and LUMO quasiparticle energies are 0.006 eV and 0.001 eV respectively.
  • Makov and Payne 1995 Makov, G.; Payne, M. C. Periodic boundary conditions in ab initio calculations. Physical Review B 1995, 51, 4014–4022.
  • 68 Experimental data are taken from the WEST GW100 data collection website (www.west-code.org/database).
  • Krause et al. 2015 Krause, K.; Harding, M. E.; Klopper, W. Coupled-cluster reference values for the GW27 and GW100 test sets for the assessment of GW methods. Molecular Physics 2015, 113, 1952–1960.
  • 70 For \ceKH molecule, calculation for the HOMO converged to a satellite instead of the quasiparticle peak. The spectral function of \ceKH is plotted and discussed in SI Sec. 6 and the correct quasiparticle energy is used here.
  • 71 For \ceKH molecule, calculation for the HOMO converged to a satellite instead of the quasiparticle peak. The spectral function of \ceKH is plotted and discussed in SI Sec. 6 and the correct quasiparticle energy is used here.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description