Full-frequency GW without frequency

Green’s function approaches based on time-dependent many-body perturbation theory provide an economical description of excitation energies and spectral intensities. For the one-particle Green’s function, which describes electron addition and removal processes, the GW approximation to the self-energy¹ performs well for weakly correlated insulators and metals,² which has partially motivated its application to molecules.^3–6 Here and throughout, we are considering the common non-self-consistent G₀W₀ approximation, which we call the GW approximation for simplicity. The size of systems that can be studied with the GW approximation is determined by the implementation, which can be characterized by its asymptotic scaling with the system size N, ranging from O(N³) to O(N⁶) with widely varying prefactors.^7,8

GW implementations can be distinguished based on their handling of a numerically challenging frequency integral, which is relatively uncommon in quantum chemical methods. The earliest works used a generalized plasmon pole model to approximate the dielectric function and, thus, integrate analytically.^2,9,10 More sophisticated approaches treat the full-frequency dependence using numerical integration techniques such as analytic continuation^11–16 and contour deformation.^15–19 These latter methods introduce numerical errors, but ones that, in principle, can be eliminated with an increase in the cost (e.g., the frequency integration grid or the fitting of the self-energy on the imaginary frequency axis). The final class of methods is numerically exact within a given single-particle basis set and requires the explicit enumeration of all neutral excitation energies,^3–5,20,21 typically calculated within the random-phase approximation (RPA). This explicit enumeration, i.e., a sum over states, dominates the cost of such a GW calculation due to its O(N⁶) scaling. This exact handling of the full-frequency dependence is the type that we address in the present work. We note that this class of methods still constructs a frequency-dependent self-energy, which is used to solve the quasiparticle equation for each excitation. In this work, we present a new formulation of the GW approximation by recasting it as an eigenvalue problem in an expanded space, and a frequency variable never appears.

Within Green’s function theories, charged excitation energies E_n, i.e., ionization potentials (IPs) and electron affinities (EAs), are found as the poles of the one-particle Green’s function matrix G(ω) via the eigenvalue problem,

where f = h + J + K is the Fock matrix, h is the kinetic and external potential energy matrix, J is the Hartree matrix, K is the exchange matrix, and Σ_c(ω) is the correlation part of the self-energy matrix. In practice, we typically work in a basis of orbitals ϕ_p(x) that diagonalize a mean-field Green’s function, which serves as the reference and defines the orbital energies ε_p. As usual, the O occupied orbitals will be indexed by i, j, k, l; the V unoccupied orbitals will be indexed by a, b, c, d; and generic orbitals will be indexed by p, q, r, s. For simplicity, we will assume real orbitals. In this basis, we have f = ε + K − V_xc, where V_xc is the exchange–correlation potential matrix. Note that for a HF reference, K − V_xc = 0.

In the GW approximation, the self-energy is calculated to lowest-order in the screened Coulomb interaction W, which gives rise to the aforementioned frequency integral, Σ_c(ω) = (i/2π)∫dω′e^iηω′G(ω + ω′)W_p(ω′), where W_p = W − v is the polarized part of the screened Coulomb interaction. When the polarizability that enters W_p is expressible by a spectral representation, the frequency integration can be performed analytically to yield^3–5,20,21

where $Mpqν=∫dx1dx2ϕp(x1)ϕq(x1)r12−1ρν(x2)$⁠, Ω_ν are neutral excitation energies, and ρ_ν(x) are transition densities.

Although any theory of neutral excitations can be used to calculate the polarizability,²² here, we consider the Tamm–Dancoff approximation (TDA) to the (direct) RPA. Within the TDA, the neutral excitation energies and transition density moments are defined by AX^ν = Ω_νX^ν, where

and $Mpqν=∑iaXiaνpi|qa$⁠. Diagrammatically, such a self-energy has screening due to infinite-order, forward time-ordered ring (or bubble) diagrams. The algebraic form [Eq. (2)] assumes that all eigenvalues and eigenvectors of the A matrix have been calculated, which implies a canonical O(N⁶) scaling, as discussed in the Introduction.

In the GW community, RPA screening is much more commonly implemented without the TDA. Although the frequency-free implementation of the GW approximation that we present here is far simpler to formulate within the TDA, many of the same ideas can be applied for the case of RPA screening, which we discuss in the supplementary material. In particular, we show that a frequency-free formulation with RPA screening exists; however, it is less conducive to cost reductions. In Fig. 1, we show that when a density functional theory (DFT) reference is used, the IPs and EAs obtained with TDA screening exhibit errors that are roughly twice as large as those obtained with RPA screening. However, for molecules, a HF reference significantly outperforms DFT references, as shown previously.²³ Remarkably, we find that the most accurate results of all methods considered are those obtained when a HF reference is combined with TDA screening, achieving an accuracy of about 0.1 eV and empirically justifying our focus on TDA screening. These conclusions are likely to change when the GW approximation is applied to solids, especially metals and small gap semiconductors that are highly polarizable.

In order to make progress on a frequency-free implementation that avoids the explicit sum over states in Eq. (2), we define a vector space of excitations corresponding to one hole (1h), one particle (1p), two holes and one particle (2h1p), and two particles and one hole (2p1h). A vector in this space has elements R = (r_i, r_a, r_i[jb], r_[jb]a). The notation of the 2h1p and 2p1h amplitudes indicates that the j → b excitation is independent of the other particle or hole index, i.e., the amplitudes do not obey any antisymmetry as they do in determinantal approaches. We define a frequency-independent “super-matrix” H,

where C^2h1p = ε^1h ⊕ (−A) and C^2p1h = ε^1p ⊕ A, with matrix elements

This super-matrix can be downfolded into the 1h + 1p space, leading to a frequency-dependent eigenvalue problem of the form Eq. (1), with

It is straightforward to check that this frequency-dependent matrix, arising from the downfolding of the 2h1p and 2p1h spaces, is precisely the correlation part of the GW self-energy; a short proof of the downfolding equivalence is given in the supplementary material. This presentation closely follows the algebraic diagrammatic construction (ADC) literature.^29,30 In particular, the above-mentioned theory, i.e., the GW approximation with TDA screening, is a strict but severe approximation to the ADC(3) and 2p1h-TDA methods.³¹ Diagrammatically, the latter two theories include many vertex corrections beyond the GW approximation, including ladder and exchange diagrams that appear as additional terms in Eq. (6). An analogous approach was also used recently to formulate an efficient renormalized second-order Green’s function theory,³² and similar conceptual ideas were discussed in the context of double excitations in time-dependent density functional theory.³³ 

Importantly, the frequency-independent super-matrix form of the GW approximation enables the use of iterative eigensolvers that lowers the computational scaling. Matrix-vector multiplication is given by HR = σ, with

where all indices correspond to spin-orbitals. For a restricted, closed-shell reference, spin-free equations are straightforward to derive and are given in the supplementary material. Clearly, the these equations exhibit no worse than O(N⁵) scaling (specifically, O²V³ for moderately sized basis sets), which is a significant improvement over the O(N⁶) scaling exhibited by the sum-over-states implementation. Furthermore, because only Coulomb-type electron repulsion integrals are used in the direct TDA (or RPA), the scaling of the most expensive contractions can be easily reduced by density fitting. For example, if the ERIs are approximated as $(pq|rs)≈∑QBpqQBrsQ$⁠, then the worst-scaling O(N⁵) term can be calculated by

which has two steps that scale as O(N_auxOV²) or O(N⁴).

Although we will not show results here, we briefly describe how spectral quantities can also be obtained iteratively with identical scalings. Using a spectral resolution of H, the full Green’s function is given by $Gpq(ω)=∑n(rpnrqn)/(ω−En)$⁠, i.e., the quasiparticle weight is given simply in terms of the 1p + 1h elements of the solution vector R; this formulation naturally precludes the common diagonal approximation Σ_pq(ω) ≈ δ_pqΣ_pp(ω). The matrix H can also be used iteratively (without diagonalization) to calculate the frequency-dependent self-energy Σ(ω) = V^2h1pZ^2h1p(ω) + V^2p1hZ^2p1h(ω), where Z(ω) is the matrix that solves the linear systems of equations, e.g., $ω12h1p−C2h1pZ2p1h(ω)=[V2h1p]†$⁠, which can be solved with iterative methods such as a conjugate gradient or the generalized minimum residual method. Similarly, the Green’s function can be calculated as G(ω) = PZ(ω), where Z(ω) solves $ω1−HZ(ω)=P†$ and P is the matrix that projects onto the 1p + 1h space.

We have implemented the GW techniques described above in the PySCF software package.^24,25 To compare their costs and verify their asymptotic scaling, we have calculated the first IP of a series of linear alkanes in the def2-SVP basis²⁶ up to C₃₇H₇₆, which has 898 basis functions The execution timings of the O(N⁶) sum-over-states, O(N⁵) frequency-free, and O(N⁴) density-fitted frequency-free implementations are shown in Fig. 2; all calculations were performed on a single core of an Intel Xeon Gold 6126 2.6 GHz (Skylake) central processing unit (CPU), and density-fitted calculations used the def2-SVP-JKFIT auxiliary basis set.³⁴ As can be seen, all methods exhibit the expected asymptotic scaling. Comparing the absolute execution times of the sum-over-states and density-fitted frequency-free implementations, we obtained a speed-up of four orders of magnitude for the C₁₀H₂₂ calculation. For our largest system with almost one thousand basis functions, the density-fitted implementation required only two hours on a single core, demonstrating the immense savings available with the advances described here. The use of density fitting was found to introduce a negligible error of around 0.01 eV.

For a rough comparison with the conventional GW approaches that achieve O(N⁴) scaling, we have calculated the first IP of C₂₅H₅₂ (610 basis functions), using analytic continuation and contour deformation, as implemented in PySCF.¹⁶ These two methods required 31 min and 64 min, respectively (with 100 grid points and 50 Padé approximants); the total execution time of our density-fitted matrix-vector implementation was 13 min. A precise comparison among the three methods is complicated because execution times depend on the number of grid points, the number of Padé approximants, and the number of poles that are calculated. In any case, our new implementation that exactly treats the frequency integration is clearly competitive with existing approaches that approximate the frequency integration.

One challenge using iterative eigensolvers on H is that the eigenvalues of typical interest (valence ionization potentials and electron affinities) are interior eigenvalues. Therefore, they must be found using energy-targeting methods such as shift-and-invert or ones that maximize the eigenvector overlap with a given guess vector. We have found the latter to work well, in conjunction with Davidson diagonalization,^35,36 for valence IPs and EAs. However, two simple alternatives exist by introducing additional approximations.

In a first approach, one can make the diagonal approximation to the Green’s function and the self-energy, seeking the self-consistent solution E_n of the algebraic equation f_pp + Σ_pp(ω = E_n) = E_n. This can be solved by iterative diagonalization of a modified matrix H^(p), which has deleted all 1p + 1h rows and columns except that of orbital p. Although the principal eigenvalue of interest is still an interior eigenvalue, this approach eliminates all other quasiparticle energies, which can facilitate energy- or overlap-targeting procedures. In a second approach, one can perturbatively decouple the IP and EA parts of H, which will make the valence IPs and EAs into extremal eigenvalues. For example, for the calculation of IPs, we perturbatively eliminate the 1p and 2p1h subspaces based on their lowest-order influence on the 1h subspace and likewise for EAs. This leads to modified Fock matrix elements,

which are used in place of f_ii and f_aa during the matrix-vector product. This approach has the added benefit of reducing the size of the vector spaces for the IP and EA problems, R^IP = (r_i, r_i[jb]) and R^EA = (r_a, r_[jb]a), and reducing the scaling of the IP matrix-vector product to be O(O³V²). Finally, we point out that the valence IPs can be made into the lowest eigenvalues by negating the matrix, H^IP → −H^IP. With all these changes, the calculation of IPs and EAs within the GW approximation looks quite similar to that within the IP/EA-EOM-CCSD approximation.^37,38 The decoupling of the IP and EA spaces is also common in the ADC literature and referred to as a non-Dyson approach.³⁹ We will use the same language and evaluate the performance of the Dyson (coupled IP and EA) and non-Dyson (perturbatively decoupled IP and EA) GW approximation.

To assess the effect of the diagonal approximation and perturbative decoupling, we used our frequency-free O(N⁵) GW implementation to calculate the first IP of all molecules in the GW100 test set.⁶ In Fig. 3, we compare the non-diagonal, diagonal, and non-diagonal perturbatively decoupled (“non-Dyson”) GW results among themselves and to ΔCCSD(T) results,²⁷ using a HF reference. As shown in Fig. 1, the GW results have a mean absolute error of 0.24 eV, with respect to CCSD(T). On average, the diagonal approximation has negligible effect (less than 0.1 eV), although a maximum deviation of 0.65 eV is observed, indicating the potential importance of off-diagonal elements of the self-energy for some molecules; these results are in agreement with previous studies of off-diagonal contributions to the GW self-energy.⁴⁰ In contrast, perturbative decoupling (the non-Dyson GW approximation) changes the results by 0.41 eV on average (and by as much as 2.44 eV) and increases the mean absolute error from 0.24 eV to 0.51 eV, suggesting that it is a relatively severe approximation. We have performed the same analysis (not shown) for the EA, as well as for a PBE reference; all results are qualitatively similar.

To summarize, we have shown that the typical Dyson equation formulation of the GW approximation can be exactly reformulated as a frequency-independent eigenvalue problem in an expanded space. In addition to providing a new conceptual framework for the GW approximation and related Green’s function theories, the new formulation was used to reduce the computational scaling from O(N⁶) to O(N⁴). Based on our preliminary results, we expect that this frequency-free formulation of the GW approximation will be readily applicable to systems with hundreds or thousands of atoms, likely limited by the memory needed to store three-index quantities.

We anticipate that this eigenvalue formulation of the GW approximation will lead to new methodological developments inspired by using quantum chemical methods with the similar structure. For example, we are exploring the use of partitioning schemes to mitigate the cost of large basis sets,^41–43 the introduction of vertex corrections through the algebraic diagrammatic construction,^29,30 the use of renormalization and compression to perform self-consistent GW,³² and the extension toward strongly correlated systems with multi-reference techniques.^44,45 The above-mentioned directions can be combined with the Galitski–Migdal formula to obtain ground-state energies, which has overlap with previous works that eliminate a frequency integration in the RPA correlation energy.⁴⁶ Finally, the ideas presented in this work and these latter extensions can be applied to treat the frequency dependence of the Bethe–Salpeter equation for neutral excitations.^47,48