Front. Chem. Frontiers in Chemistry Front. Chem. 2296-2646 Frontiers Media S.A. 736591 10.3389/fchem.2021.736591 Chemistry Original Research Low-Order Scaling Quasiparticle Self-Consistent GW for Molecules Förster and Visscher Low-Order Scaling qsGW for Molecules Förster Arno * Visscher Lucas Theoretical Chemistry, Vrije Universiteit, Amsterdam, Netherlands

Edited by: Dorothea Golze, Aalto University, Finland

Reviewed by: Jan Wilhelm, University of Regensburg, Germany

Ivan Duchemin, Commissariat à l'Energie Atomique et aux Energies Alternatives (CEA), France

 *Correspondence: Arno Förster, a.t.l.foerster@vu.nl

This article was submitted to Theoretical and Computational Chemistry, a section of the journal Frontiers in Chemistry

 03 09 2021 2021 9 736591 05 07 2021 13 08 2021 Copyright © 2021 Förster and Visscher. 2021 Förster and Visscher

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Low-order scaling GW implementations for molecules are usually restricted to approximations with diagonal self-energy. Here, we present an all-electron implementation of quasiparticle self-consistent GW for molecular systems. We use an efficient algorithm for the evaluation of the self-energy in imaginary time, from which a static non-local exchange-correlation potential is calculated via analytical continuation. By using a direct inversion of iterative subspace method, fast and stable convergence is achieved for almost all molecules in the GW100 database. Exceptions are systems which are associated with a breakdown of the single quasiparticle picture in the valence region. The implementation is proven to be starting point independent and good agreement of QP energies with other codes is observed. We demonstrate the computational efficiency of the new implementation by calculating the quasiparticle spectrum of a DNA oligomer with 1,220 electrons using a basis of 6,300 atomic orbitals in less than 4 days on a single compute node with 16 cores. We use then our implementation to study the dependence of quasiparticle energies of DNA oligomers consisting of adenine-thymine pairs on the oligomer size. The first ionization potential in vacuum decreases by nearly 1 electron volt and the electron affinity increases by 0.4 eV going from the smallest to the largest considered oligomer. This shows that the DNA environment stabilizes the hole/electron resulting from photoexcitation/photoattachment. Upon inclusion of the aqueous environment via a polarizable continuum model, the differences between the ionization potentials reduce to 130 meV, demonstrating that the solvent effectively compensates for the stabilizing effect of the DNA environment. The electron affinities of the different oligomers are almost identical in the aqueous environment.

 GW approximation convergence acceleration analytical continuation quasiparticle quasiparticle self-consistent GW DNA photodamage theoretical spectroscopy Nederlandse Organisatie voor Wetenschappelijk Onderzoek10.13039/501100003246 1 Introduction

The GW approximation (GWA) to Hedin’s equations (Hedin, 1965) is a popular approach to calculate charged excitations in molecular systems. Recent applications include the calculation of band gaps and elucidation of charge-transfer in organic donor-acceptor compounds (Blase and Attaccalite, 2011; Blase et al., 2011; Caruso et al., 2014), applications to dye-sensitized solar cells (Marom et al., 2011; Faber et al., 2012; Umari et al., 2013; Marom et al., 2014; Mowbray and Migani, 2015), electronic level alignment in photocatalytic interfaces (Migani et al., 2013, 2014), core-ionization spectra of medium sized molecules (Van Setten et al., 2018; Golze et al., 2018, 2020) or photo-electron spectra of transition metal oxides (Berardo et al., 2017; Hung et al., 2017; Shi et al., 2018; Rezaei and Ögüt, 2021). Combined with the Bethe-Salpeter equation (BSE) formalism (Salpeter and Bethe, 1951; Strinati, 1988) the GWA has been used to calculate optical spectra of Cyanins (Boulanger et al., 2014), the Bacteriochlorin molecule (Duchemin et al., 2012) or Bacteriochlorophylls and Chlorophylls (Hashemi and Leppert, 2021). At the same time, the GWA has been implemented into an increasing number of molecular electronic structure codes (Ke, 2011; Caruso et al., 2012; Caruso et al., 2013; Ren et al., 2012; Van Setten et al., 2013; Kaplan et al., 2015, 2016; Bruneval et al., 2016; Wilhelm et al., 2016; Tirimbò et al., 2020b). Traditionally, these implementations use localized basis functions and the resolution-of the identity or density fitting approximation (Baerends et al., 1973; Whitten, 1973; Dunlap et al., 1979) within the global Coulomb metric (RI-V) (Vahtras et al., 1993), leading to a scaling of N 4 with system size. Systems of around 100 atoms are within reach on standard hardware (Knight et al., 2016), while highly parallel implementations enable applications to systems with more than 300 atoms on modern supercomputers (Wilhelm et al., 2016; Wilhelm et al., 2018; Wilhelm et al., 2021).

Over the last years, many algorithms with reduced asymptotic scaling with system size have been proposed. These are usually based on the space-time approach by Godby and coworkers (H. N. Rojas et al., 1995; Rieger et al., 1999). The original space-time method is based on the observation that it is much simpler to solve the Dyson equations in the GWA in reciprocal space and imaginary frequency while the kernels of these Dyson equations are most easily evaluated in real space and imaginary time, reducing the asymptotic scaling of the GWA to N 3. Building on earlier work by Almlöf (Almlöf et al., 1982), Kresse, Kaltak and coworkers could significantly reduce the prefactor of these calculations by using non-uniform spaced grids in imaginary time and imaginary frequency and an efficient way to switch between both domains (Kaltak et al., 2014a; Kaltak et al., 2014b; Kaltak and Kresse, 2020). Over the last years, there has been a surge of new GW implementations based on the space-time method for periodic (Kutepov et al., 2012; Chu et al., 2016; Liu et al., 2016; Kutepov et al., 2017; Grumet et al., 2018; Kutepov, 2020; Singh and Wang, 2020; Foerster and Gueddida, 2021) and finite (Wilhelm et al., 2018; Koval et al., 2019; Förster and Visscher, 2020; Duchemin and Blase, 2021; Wilhelm et al., 2021) systems. Other recent examples of low-order scaling implementations include the spectral function based approach by Foerster et al. (2011), the time-shredded propagator formalism by Ismail-Beigi and coworkers (Kim et al., 2020), stochastic GW developed by Neuhauser et al. (2014), Vlček et al. (2017), Vlček et al. (2018), Weng and Vlcek (2021), and also a fragment molecular orbital based implementation (Fujita et al., 2019).

For molecular systems, diagonal approximations to the self-energy are commonly made. They rely on the assumption that the wave function of generalized Kohn-Sham (KS) density functional theory (DFT) is similar to the GW wave function. One then evaluates corrections to the DFT single orbital energies by calculating the diagonal elements of the self-energy matrix Σ. The most economical way to calculate these corrections is the one-shot G 0 W 0 approach which heavily depends on the mean-field starting point. Extensive benchmarks (Marom et al., 2012; Bruneval and Marques, 2013; Caruso et al., 2016; Knight et al., 2016) have provided substantial evidence that hybrid functionals with a rather large amount of exact exchange or long-range corrected hybrids are usually a suitable starting point. In addition, non-empirical procedures to select an optimal starting point for a given system have been proposed (Gallandi and Körzdörfer, 2015; Dauth et al., 2016; Bois and Körzdörfer, 2017). Finally, in eigenvalue-only self-consistent GW (evGW) the QP energies are updated until they are stationary, removing the starting point dependence to a large extent.

QP energies calculated following these strategies are almost always more accurate than fully self-consistent GW (scGW) calculations for molecules. As discussed by Kotani, van Schilfgaarde and Valeev, QP approximations, i.e. approximations in which satellites are neglected, emphasize the importance of the Ward identity (Ward, 1950) in the long-range and low-frequency limit. The Ward identity demands ’Z-factor cancellation’ (Kotani et al., 2007) between the three-point Vertex and the renormalized electron propagator. Z is the QP renormalization factor. In QP approximations, neither the vertex is included nor is the propagator renormalized, and the effect of both approximations cancel in the above-mentioned limit. This limit can be expected to be of particular importance for weakly correlated molecules to which the GWA is frequently applied.

As opposed to diagonal approximations, scGW is strictly starting point independent and also allows to calculate 1-particle reduced density matrices (1RDM) including electron correlation effects from first principles. Most importantly, it does not contain any adjustable parameters. Another method which also offers these advantages is the QP self-consistent GW (qsGW) method by Kotani, van Schilfgaarde and Faleev. (Van Schilfgaarde et al., 2006; Kotani et al., 2007). qsGW can be seen as a non-empirical procedure to find an optimal starting point for a G 0 W 0 calculation. This is accomplished by mapping the GW self-energy self-consistently to a non-local, Hermitian, and static exchange-correlation potential. This potential has been shown to be optimal in a variational sense (Ismail-Beigi, 2017). Diagonalization of the resulting mean-field Hamiltonian yields eigenvectors and eigenvalues from which a new non-interacting Green’s function is obtained. This self consistent field (SCF) procedure is reminiscent of generalized KS theory, with the notable difference that the exchange-correlation potential is not a functional of the 1RDM but rather of the non-interacting single-particle Green’s function. qsGW is starting point independent and fulfills the Ward identity in the low frequency and long range limit.

In canonical implementations (Ke, 2011; Bruneval, 2012; Koval et al., 2014; Kaplan et al., 2015; Kaplan et al., 2016), the need to calculate the off-diagonal elements of the self-energy matrix and the fact that it is typically more difficult to converge make qsGW typically an order of magnitude more expensive than evGW (Gui et al., 2018) which in turn is typically 5–10 times more expensive than G 0 W 0 due to the requirement of self-consistency. Moreover, low-order scaling implementations for molecules are typically restricted to diagonal approximations only (Wilhelm et al., 2018; Förster and Visscher, 2020; Duchemin and Blase, 2021; Wilhelm et al., 2021) To fill this gap, we extend the recently developed low-order scaling diagonal GW implementation in ADF (Baerends et al., 2020; Förster and Visscher, 2020, Förster and Visscher, 2021) to qsGW. We evaluate the qsGW self-energy as a direct product in imaginary time, in the same way as in the diagonal approximation. Even though the qsGW self-energy is static, for larger systems evaluation of the self-energy at an array of imaginary time points is more efficient than its evaluation at a single real frequency point. The procedure is similar to the linearized qsGW method by Kutepov and coworkers (Kutepov et al., 2017) which is also based on the imaginary time formalism and in which the self-energy is averaged over all frequencies. However, in our implementation, we only average over frequencies for the off-diagonal elements but retain the optimum exchange-correlation potential on the diagonal. We achieve stable and rapid convergence of the SCF procedure by a suitable implementation of the direct inversion in the iterative subspace (DIIS) Pulay (1980) approach. Most importantly, the proposed algorithm is easy to implement and only requires to combine the qsGW approach with the space-time implementation for the self-energy and an efficient method to evaluate the exact exchange-contribution to the Fock matrix.

This work is organized as follows: In section 2 we first recapitulate the qsGW procedure and describe some aspects of our implementation. We focus on the implementation of the DIIS and on the analytical continuation (AC) of the self-energy. In section 3, we confirm the correctness of our implementation by comparison to ionization potentials (IP) (Kaplan et al., 2016) from TURBOMOLE (Balasubramani et al., 2020) and investigate the convergence of the SCF equations. We also illustrate the computational performance of our implementation with a proof-of principle application to large DNA oligomers. In section 4 we summarize and conclude this work.

 2 Methods

In this section, we review the qsGW method and comment on our implementation, focusing on the AC of the self-energy as well as our approach to accelerate convergence of the SCF procedure. Greek lowercase letters μ, ν … label atomic orbitals (AO) and run from 1 to n AO . Latin lowercase letters p, q, r, … label general MOs and run from 1 to n MO . i, j, k (a, b, c) label occupied (virtual) MOs and run from 1 to N occ (N virt ). Latin symbols without labels denote tensors in some basis which will always be clear from the context.

2.1 QP Self-Consistent GW

The GWA is an approximation to the self-energy appearing in Dyson’s equation (Dyson, 1949), ∑ r Σ p r ( ω p ) U r q ( ω p ) = ω p − ϵ p U p q ( ω p ) .

We mostly work in a basis of molecular orbitals (MO), ϕ p ( n ) ( r ) = ∑ μ χ μ ( r ) b μ p ( n ) , where the χ μ are AOs. Dysons’s equation is non-linear and will be solved via a fixed point iteration. The superscript (n) means that we are in the nth iteration of a SCF procedure. The self-energy Σ is non-Hermitian and energy dependent. Thus, U is complex and energy dependent as well. We will neglect spin in the following.

The ε p are obtained from solving the generalized KS problem, ∑ ν H μ ν ( 0 ) b ν p ( 0 ) = ∑ ν S μ ν b ν p ( 0 ) ϵ p ( 0 ) , H K S = T + V e x t + V H x c [ P ] , where V Hxc is the sum of exchange-correlation potential V xc and Hartree potential V H , being functionals of the 1RDM P and the electron density, respectively. T and V ext are kinetic energy and external potential, respectively. S is the overlap matrix of AOs and b defines a transformation from AO to MO basis, M p q = b p μ M μ ν b † ν q . In the AO basis, P is given as P μ ν = 2 ∑ i N o c c b μ i b † i ν . We also define the Hamiltonian of the Hartree approximation, H H = H K S − V x c . The Green’s function G 0 corresponding to the non-interacting Hamiltonian is diagonal in the MO basis with G 0 p p ( ω ) = i ω − ϵ p − 1 . We can then expand Σ in terms of G 0 as follows (Hedin, 1965; Martin et al., 2016), Σ ( ω ) = ( G 0 * W 0 ) ( ω ) + … , and in the GWA the expansion is truncated after first order. W 0 is the screened Coulomb interaction, calculated in the bubble approximation (Onida et al., 2002) from G 0 (Hedin, 1965). Without further approximations to Σ, one typically avoids solving (Eq. 1) but instead calculates the interacting Green’s function G by inversion of G ( ω ) − 1 = G 0 ( ω ) − 1 − Σ ( ω ) . From there one proceeds by building the self-energy (Eq. 8) but replaces G 0 by G, and W 0 by W and repeats this procedure until self-consistency is reached. In more approximate GW schemes, one avoids solving (Eq. 9). In diagonal approximations to Dysons’s equation, one assumes Σ to be diagonal. In that case, U in (Eq. 1) is unity for all ω and (Eq. 1) reduces to a set of independent non-linear equations for ω. In qsGW on the other hand, one does not make the diagonal approximation but Σ is mapped to a Hermitian and frequency-independent exchange-correlation potential V x c q s G W . For this mapping, it is convenient to define W 0 ( ω ) = V c + W ̃ 0 ( ω ) , with V c being the bare Coulomb potential. The self-energy can then be decomposed into a static and dynamic part Σ ( ω ) = Σ x + ( G 0 * W ̃ 0 ) ( ω ) = Σ x + Σ c ( ω ) . Σ x is the Fock exchange potential, V x q s G W = Σ x , and following Kotani et al. (2007), the correlation part of V xc is obtained from Σ by taking one of the real symmetric definitions V c q s G W p q = 1 2 R e Σ c p q ( ϵ p ) + R e Σ c p q ( ϵ q ) , or V c q s G W p q = δ p q R e Σ c p q ( ϵ p ) + ( 1 − δ p q ) R e Σ c p q ( ω = 0 ) . There are formal reasons why (Eq. 12) should be preferred over (Eq. 13). Constructing the qsGW Hamiltonian via (Eq. 12) minimizes the length of the gradient of the Klein functional (Klein, 1961) with respect to G 0 (Ismail-Beigi, 2017) and can be seen as an optimized effective non-local potential. The approach bears strong resemblance to what is usually referred to as the optimized effective potential (OEP) method (Talman and Shadwick, 1976). Another possibility is to linearize the self-energy around the chemical potential. This has been implemented by Kutepov et al. (2017). Physically, it is equivalent to taking the static limit of the self-energy, or averaging over frequencies. We will discuss in more detail below that such an approach has advantages with regards to numerical stability. However, we think that one should use the optimum potential at least for the diagonal elements. (Eq. 13) is a hybrid between (Eq. 12) and Σ (ω = 0) which retains the optimum potential on the diagonal. Employing (Eq. 13) can be justified if one assumes that the effect of using the optimum potential as opposed to Σ (ω = 0) will cancel out to a large extent for the off-diagonal elements. We provide numerical evidence later on that this is indeed true. Also an approach using Löwdin’s orthogonalization has been proposed to construct the QP Hamiltonian (Sakuma et al., 2009) but that construction is not considered here.

With these simplifications, we can now solve (Eq. 1) self-consistently. In each iteration, we solve ∑ r H p r q s G W ( n + 1 ) U r q ( n + 1 ) = ω p ( n + 1 ) U p q ( n + 1 ) , with H q s G W ( n + 1 ) = H H + Δ V H ( n + 1 ) + V x c q s G W ( n + 1 ) and V x c q s G W ( n + 1 ) = V x [ P ( n ) ] + V c q s G W [ G 0 ( n ) ] . In each iteration, H qsGW is expressed in the basis in which G 0 ( n ) is diagonal. That is, at the n + 1st iteration, H qsGW is expressed in terms of the ϕ i ( n ) and unless self-consistency has been reached, U (n) will not be unity and defines a rotation of the molecular orbitals. We now set b μ p ( n + 1 ) = ∑ q b μ q ( n ) U q p ( n + 1 ) ϵ p ( n + 1 ) = ω p ( n + 1 ) ∀ p and evaluate G 0 ( n + 1 ) via (Eq. 7) which in turn is used to evaluate (Eq. 11) and finally (Eq. 12) or (Eq. 13). P (n+1) is then evaluated from (Eq. 5) and the change in the Hartree-potential is calculated as Δ V H ( n + 1 ) = V H [ Δ P ( n + 1 ) ] , with Δ P ( n + 1 ) = P ( n + 1 ) − P ( n ) . The cycle is repeated until self-consistency is reached.

 2.2 Implementation

As already stressed in the introduction, for the qsGW implementation no modifications of the code described in Förster and Visscher (2020) for the calculation of the self-energy are needed. A description of the algorithm can be found in Förster and Visscher (2020) and in Förster and Visscher (2021) we reported important modification of our original implementation, increasing accuracy and robustness. The only points we discuss hered are related to the convergence and stability of the self-consistent field (SCF) procedure.

2.2.1 Analytical Continuation

In space-time implementations of the GWA, the self-energy is evaluated in imaginary time and then Fourier transformed to the imaginary frequency axis. In ADF, the self-energy is calculated in the AO basis on a non-uniform grid of imaginary time points. After transformation to the reference basis [the MO basis from the generalized KS calculation in the first iteration and the basis defined by (Eq. 17) later], the self-energy matrix is Fourier transformed to a non-uniform grid in imaginary frequency space. For the implementation of this transformation, we refer to Kaltak et al. (2014b) and to the appendix of Förster and Visscher (2021). Since the non-uniform grids depend on the QP energies used to build G 0 we also need to recalculate these grids at the beginning of each qsGW iteration to ensure independence of the results from the initial guess.

After this transformation, Σ is known on a discrete set of points W = i ω β β = 1 , N ω on the imaginary frequency axis. However, to evaluate Eq. 13, we need to know the self-energy on the real frequency axis at the positions of the QP energies ϵ p ( n ) . To this end, we seek to find a function f which is analytic in the largest possible domain A ⊂ C and coincides with Σ in W . For a meromorphic function (as the self-energy) which is known on the whole imaginary axis, it is always possible to find such a function so that A = C , but since we only know the self-energy on a small subset of points, only an approximate solution can be found. The problem here is, that the AC is exceptionally ill-conditioned, i.e. numerical noise in the input data might significantly affect the output (Shinaoka et al., 2017).

Among the many developed algorithms [see for instance Levy et al. (2017) for an overview], the construction of a continued fraction (Vidberg and Serene, 1977; Beach et al., 2000) via a Padé approximant is most common in implementations of the GWA. While in many codes Thiele’s reciprocal difference method is implemented, (Liu et al., 2016; Grumet et al., 2018; Foerster and Gueddida, 2021), ADF, implements the variant by Vidberg and Serene (Vidberg and Serene, 1977), which for example has also been implemented by Kutepov (Kutepov, 2020). In the latter variant, the coefficients of the continued fraction are calculated while the former method returns the value of the continued fraction (Beach et al., 2000). While it has been claimed that the former variant is numerically more stable (Liu et al., 2016), we did not experience any numerical issues with our implementation for diagonal self-energies. This procedure typically yields good results for states close to the HOMO-LUMO gap while it becomes unreliable for core states (Golze et al., 2018, 2020). Exceptions are cases for which the self-energy has a pole close to the position of the QP energy (Govoni and Galli, 2018). Partial self-consistency in G pushes the poles away from the QP peak (Golze et al., 2019), and consequently, these issues should not be present in qsGW as well. This is different from situations in which the independent QP picture breaks down and the spectral weight of a single excited electrons is distributed between multiple peaks. The former is a purely numerical issue while the latter is caused by strong correlation and can not be overcome by partial self-consistency. It has also been shown in Wilhelm et al., 2021 that AC yields accurate results for semi-core and inner valence states in case the real part of the self-energy does not have poles in the vicinity of the QP solutions.

If one is only interested in accurate valence states, AC via Padé approximants is not problematic for G 0 W 0 where (Eq. 1) reduces to a set of N independent non-linear equations where N is the number of MOs. In evGW, the situation is only slightly different. The N equations are still independent, but information from all QP energies enters the polarizability so that there is an implicit dependence of the QP energies on each other. In practice, this is also not an issue since the numerical errors are typically orders of magnitude smaller than the absolute values of the QP energies.

The situation is different for qsGW. The Off-diagonal elements of Σ c are often equal to or very close to zero (Kaplan et al., 2015) and generally small compared to the diagonal elements. For these off-diagonal elements, numerical errors from AC can be orders of magnitudes larger than the values of the off-diagonal elements. Since there are many of them, this might significantly alter the solutions of Eq. 14. Due to the non-linear nature of the QP equations, this can complicate convergence of the SCF procedure or even lead to erroneous results. The development of more reliable methods for AC is a very active field of research (Bergeron and Tremblay, 2016; Levy et al., 2017; Otsuki et al., 2017; Gull et al., 2018; Fournier et al., 2020; Fei et al., 2021) and it would certainly be interesting to investigate whether other techniques are more suitable for qsGW. For now, we restrict ourselves to the techniques of Padé-approximants. To ensure numerical stability, two aspects need to be considered:

First, it seems reasonable to assume that AC close to the Fermi energy is also more reliable for the off-diagonal elements of Σ. To this end, using (Eq. 13) to construct the exchange-correlation potential seems to be more suitable for our implementation than (Eq. 12). As we will see later on, both constructions of the exchange-correlation potential lead to similar results, but using (Eq. 13), the SCF procedure is significantly easier to converge. In fact, applying the same reasoning one could justify to use Σ (ω = 0) (Kutepov et al., 2017) instead. However, as we will show below, using (Eq. 13) is sufficiently numerically stable.

Second, after evaluating Eq. 13 or (Eq. 12), numerical noise needs to be removed rigorously from V c q s G W . At self-consistency, the off-diagonal elements of V c q s G W need to be zero: In the n + 1 the iteration, V c q s G W is expressed in the basis which diagonalizes the operator defined in (Eq. 15) in the nth iteration. At self-consistency b (n+1) = b (n), which will not be the case when the off-diagonal elements of V c q s G W will be different from zero. In our present implementation, we set all values with magnitude smaller than 1e −6 to zero. This cut-off is of the order of the numerical noise introduced by the AC. As we will show later on, despite this drastic cut-off the HOMO and LUMO energies can be converged to a degree that the QP energies are converged within a few meV.

 2.2.2 Convergence Acceleration

As outlined so far, in each iteration of the self-consistency cycle the previous qsGW Hamiltonian is replaced by the new one, similar to the Roothaan algorithm for the Hartree-Fock (HF) equations. For Hartree-Fock, it is well known, that such a procedure can be numerically unstable (Cances and Le Bris, 2000) and convergence difficulties are encountered already for the simplest molecules (Koutecký and Bonačić, 1971; Bonačić-Koutecký and Koutecký, 1975). Also in many GW implementations, convergence has been shown to be much slower than with a simple linear mixing scheme (Caruso et al., 2013; Kaplan et al., 2016). While the latter seems to work reasonably well for evGW (Gui et al., 2018), it seems that there is room for improvement for qsGW (Gui et al., 2018). An iterative fixed point procedure of the general form G 0 ( m ) 0 ≤ m ≤ n + 1 → H ̃ q s G W n + 1 → ϵ ( n + 1 ) , b ( n + 1 ) is clearly a better option. A practical way to implement this is to replace (Eq. 14) by ∑ r H ̃ p r q s G W ( n + 1 ) U r q ( n + 1 ) = ω p ( n + 1 ) U p q ( n + 1 ) , with H ̃ q s G W ( n + 1 ) = ∑ m = n − n 0 n + 1 α m H q s G W ( m ) , where ∑ m = n − n 0 n α m = 1 , needs to be fulfilled and n 0 is the maximum number of previous iterations taken into account. We determine the expansion coefficients α m using Pulay’s DIIS method (Pulay, 1980). In the DIIS method, we seek to minimise the residual error r ( n + 1 ) = ∑ m = n − n 0 n α m r ( m ) , subject to the constraint Eq. 23. One might additionally require the α m to be positive (what is usually called EDIIS) (Kudin et al., 2002) but we did not find any improvement over the simple DIIS. Different implementations of DIIS differ in the definition of the residual error. Since G 0 uniquely determines H qsGW , we would ideally define r ( n + 1 ) = G 0 ( n + 1 ) − G 0 ( n ) , however, storage (or recalculation) of this quantity for n 0 iterations is inefficient. Therefore, one can use r ( n + 1 ) = P ( n + 1 ) − P ( n ) , which is related to the time-ordered Green’s function by taking the limit τ → 0− (τ is the difference between both time arguments). In this work, we have used a different definition for the residual which is, however, identical to (Eq. 26). 1

Technically, in the nth iteration we solve (Eq. 14) and evaluate the corresponding b (n) from which we calculate P (n) and Q (n). We check for convergence by evaluating the Frobenius norm of the residual (Eq. 26), N F = 1 N M O 2 ∑ μ ν r μ ν ( n + 1 ) 2 , and terminate the SCF as soon as N F < ε SCF for two subsequent iterations. As we will show later on, ε SCF = 1e −7 leads to QP energies which are converged within a few meV for all systems in the GW100 database (Van Setten et al., 2015). Subsequently, we store r (n+1) and H q s G W ( n + 1 ) and determine the expansion coefficients α m using the DIIS method, setting n 0 = 10. Finally, we solve (Eq. 21) and use the resulting U to evaluate (Eq. 17).

 2.3 Computational Details

All calculations have been performed with a locally modified development version of ADF2020 using the implementation as described Förster and Visscher (2020) and using the updated imaginary frequency grids as described in Förster and Visscher (2021).

2.3.1 GW100

We use the same structures as in for our previous benchmarks (Förster and Visscher, 2020; Förster and Visscher, 2021). We use the non-augmented TZ3P and QZ6P basis sets described in Förster and Visscher (2021). Complete basis set (CBS) limit extrapolated results are obtained as described in Förster and Visscher (2021). In all calculations, we set the numericalQuality key to Good. Exceptions are a few systems for which we observed inconsistencies with the Good fit set: For Pentasilane, Na2, Na4, and Na6, we used the Excellent fit set, and for the nucleobases we used the VeryGood fitset. We used 32 imaginary time and 32 imaginary frequency points each [We refer to the explanations in the appendix of Förster and Visscher (2021)]. For all TZ3P calculations, we set Dependency Bas = 1e−3 and for QZ6P we set Dependency Bas = 5e−3 in the AMS input as described in Förster and Visscher (2020). All calculations using augmented basis sets (aug-TZ3P and aug-QZ6P) have been performed in the same way, but using the Excellent auxiliary fit set and numericalQuality VeryGood. No relativistic effects have been taken into account.

 2.3.2 DNA Fragments

The structures of the DNA fragments have been taken from Doser et al. (2009). We performed qsGW calculations using the TZ2P (Van Lenthe and Baerends, 2003), TZ3P and QZ6P basis sets, starting from a PBE0 (Adamo and Barone, 1999; Ernzerhof and Scuseria, 1999) initial guess. We set the numerical quality to VeryGood, but used the Good fitset, with the exception of the QZ6P calculations were we also used the VeryGood fitset. We also set MBPT. ThresholdQuality = Normal. In Förster and Visscher (2020) we have shown that these thresholds are sufficient to converge quasi-particle energies within a few 10 meV. 16 grid points in imaginary time and imaginary frequency have been used. Solvent effects have been accounted for exclusively on the KS level using the conductor like screening model (COSMO) (Klamt and Schüürmann, 1993; Klamt, 1995; Klamt and Jonas, 1996) as implemented in ADF (Pye and Ziegler, 1999) using the BLYP (Becke, 1988; Lee et al., 1988; Miehlich et al., 1989) functional with D3 dispersion correction (Grimme et al., 2010) with Becke-Johnson damping (Grimme et al., 2011) and the TZ2P basis set. Numericalquality Good has been used. The solvent correction ΔE s is then obtained as Δ E s = E s ( + ) − E s ( 0 ) , i.e. as the difference between the solvent contributions to the bonding energies of the oxidized species and the neutral species both at the equilibrium geometry of the neutral species.

 3 Results 3.1 Benchmarks 3.1.1 Comparison of Exchange-Correlation Potentials in qs<italic>GW</italic>

We already noticed in section 2 that the correlated part of the exchange-correlation potential of qsGW can be defined in different ways. Here we compare the two most common ways to construct this quantity (Kotani et al., 2007; Shishkin et al., 2007; Shishkin and Kresse, 2007; Kaplan et al., 2016) (Eq. 12 and Eq. 13) for a subset of molecules from the GW100 database. The data is shown in the supporting information and shows that the exchange-correlation potential obtained from (Eq. 12) is significantly harder to converge than the one from (Eq. 13). An example of the convergence behaviour of both variants is shown in Figure 1. Figure 1 plots log 10 r with r defined in Eq. 26 against the number of iterations with two different initial guesses for Methane. We see, that using (Eq. 13), the SCF rapidly converges towards a fixed point, while log 10 r always remains much larger than −6 for (Eq. 12). On the other hand, for the 10 converged calculations differences in the final QP energies are small; for both, IPs and EAs, both variants differ by only 20 meV on average, i. e the error introduced by averaging over the off-diagonal elements of the self-energy are small. For this reason, we decided to use the correlation potential as defined in (Eq. 13) in all subsequent calculations.

Convergence of the qsGW SCF for Methane for different initial guesses and constructions of the correlation potential.  log 10 r, r defined in Eq. 26, is plotted against the number of iterations.

 3.1.2 Self Consistent Field Convergence

Next, we comment on the convergence of the qsGW SCF procedure. To this end, we compare IPs and electron affinities (EA) for the molecules in the GW100 database for 3 different starting points, PBE (Perdew et al., 1996a; Perdew et al., 1996b), PBE0, and HF. At self-consistency, the QP energies should be independent from the initial guess and their differences will thus provide information about the obtained convergence of the QP energies for a given ε SCF . In all calculations we set ε SCF = 1e −7 and restrict all calculations to a maximum of 30 iterations.

Independent of the starting point, we could not reach convergence for Mgo, BeO, BN, Cu2, and CuCN with our DIIS implementation. Employing a linear mixing procedure as implemented in Wilhelm et al. (2021) with α = 0.35 we could reach convergence for these systems, albeit with a large number of iterations. These systems are problematic for GW approaches since the single the spectral weight of the single excited electron is distributed between multiple peaks (Govoni and Galli, 2018). qsGW relies on the validity of the single QP picture. In situations, in which the quasi-particle equations might have multiple solutions (Govoni and Galli, 2018; Golze et al., 2019) corresponding to the same non-interacting state, different solutions may be found in different iterations of the qsGW SCF procedure. qsGW should select the solution with largest QP weight (Ismail-Beigi, 2017) but in situations where there are at least two solutions with (almost) equal QP weight, the “physical” solution might change in each iteration. In such cases, the DIIS algorithm tries to minimize the residual SCF error by interpolating between different solutions and no fixed point of the map (Eq. 20) is found. On the other hand, linear mixing results in a smooth but slow convergence pattern, if only α is chosen small enough to make sure that in all iterations the same solution is found. We do not know, how to best solve this issue but we do not consider it to be a major concern as such convergence problems are only encountered for systems in which the single QP picture is not valid. This then merely signals that qsGW is not an appropriate level of theory.

Figure 2 shows mean absolute deviations (MAD) as well as maximum absolute deviations of the IPs and EAa obtained from different starting points. With MAD of 6 and 2 meV, respectively, EAs are better converged than IPs. Also the maximum error is about twice as small for EAs than for IPs. These differences are related to the AC procedure which gives smaller errors for unoccupied states with usually featureless self-energy matrix elements. The maximum error never exceeds 50 meV and is of the same order of magnitude than the experimental resolution of photoionization experiments (Knight et al., 2016) of the typical basis set errors of GW calculations after extrapolation. (Knight et al., 2016; Maggio et al., 2017; Govoni and Galli, 2018; Bruneval et al., 2020; Förster and Visscher, 2021). The distribution of iterations required for convergence is displayed in Figure 3. This includes the 5 problematic cases discussed above. The calculations on average converge in around 10 iteration, with little dependence on the initial guess.

Mean absolute deviations (A) and maximum absolute deviations (B) of qsGW IPs (upper triangle) and EAs (lower triangle) obtained with different initial guesses for the GW100 database. All values are in meV.

Number of iterations needed to attain convergence of the SCF for different initial guesses.

 3.1.3 Comparison of Ionization Potentials for the GW100 Database

We now compare the IPs from our algorithm to the ones obtained with the TURBOMOLE code for GW100. The TURBOMOLE results have been obtained with the GTO-type def2-TZVPP basis sets. For some systems, TURBOMOLE results are not available and we exclude these from our discussion. We use the TZ3P basis sets which we have shown to give comparable results to def2-TZVP for GW100 (Förster and Visscher, 2021). However, quantitative accuracy can not be expected.

The deviations to TURBOMOLE are shown in Figure 4. The average deviation between both codes is close to zero, and with one exception, for all IPs deviations are considerably smaller than 300 meV, with the deviations for the majority of systems being smaller than 100 meV. Thus, our results are qualitatively similar and deviations can be attributed to different basis set errors and different constructions of the qsGW exchange-correlation potential. The IP of Cyclooctatetrane is the only exception. Here, TURBOMOLE gives an IP of 9.30 eV, while the ADF IP is with 8.38 eV nearly 1 eV smaller. For different starting points, we obtained the same result within an accuracy of only a few meV, indicating that our IP is well converged. The TURBOMOLE qsGW IPs on average overestimate the CCSD(T) reference values for GW100 by Klopper and coworkers (Krause and Klopper, 2017) in the same basis set by only a little more than 100 meV, while the deviation for Cycloocatetrane is nearly 1 eV. The CCSD (T) IP for this system, is 8.35 eV, which is in very good agreement with our value. These numbers indicate that our IP is reasonable, despite the large deviation to TURBOMOLE.

Distribution of deviations (in eV) of the IPs from TURBOMOLE and with our implementation.

Ideally, we would also like to compare our EAs against literature data, however, with only one exception (were optimized structures do not seem to be available) (Ke, 2011), we are not aware of any published EAs for molecular systems.

 3.1.4 Basis Set Limit Extrapolated Ionization Potentials and Electron Affinities for the GW100 Database

In the supporting information, we report CBS limit extrapolated EAs and IPs for the GW100 database. The qsGW QP energies seem to converge faster to the CBS limit than their G 0 W 0 counterparts. Going from TZ3P to QZ6P, the basis set incompleteness error reduces by 80 meV on average, while for G 0 W 0 @PBE, we found an average reduction of 130 meV (Förster and Visscher, 2021). Self-consistent approaches might converge faster than G 0 W 0 - Caruso et al. have already observed that scGW converges faster to the CBS limit than G 0 W 0 (Caruso et al., 2013). For the EAs, the average differences are much larger which is also due to the many systems with negative EA in the GW100 database. For these systems CBS limit extrapolation is not reliable without adding diffuse functions. Repeating these calculations with augmented basis sets (Förster and Visscher, 2021) yields smaller differences between the aug-TZ3P and aug-QZ6P basis sets. (Förster and Visscher, 2021). In Table 1, these differences are shown for the series of linear alkanes from Methane to Butane (for more numbers we refer to the supporting information). On both the TZ and QZ level the augmented basis sets give a much higher EA. Also, the differences between aug-TZ3P and aug-QZ6P are with in between 150 and 200 meV modest, while they are huge for the non-augmented basis sets. Also the extrapolated values are much smaller using the augmented basis sets. The effect of augmentation is also profound for other systems. For example, using the non-augmented basis sets, the EA of carbontetrachloride is negative (−0.27 eV). Using the augmented basis sets, it becomes positive (0.17 eV) which is in much better agreement with experiment (0.80 ± 0.34 eV) (Staneke et al., 1995).

Comparison of electron affinities for linear alkanes from Methane to Butane using augmented, and non-augmented basis sets.

Non-augmented Augmented
Name TZ3P QZ6P Extrap Aug-TZ3P Aug-QZ6P Extrap
Methane −2.30 −−1.62 −0.78 −0.79 −0.58 −0.26
Ethane −2.27 −1.56 −0.65 −0.72 −0.57 −0.35
Propane −2.23 −1.51 −0.56 −0.72 −0.55 −0.30
Butane −2.24 −1.50 −0.52 −0.71 −0.55 −0.30
3.2 Application to DNA Fragments

Oxidation of DNA is related to genetic damage and to investigate the mechanisms behind these processes quantum chemically, electron addition and removal energies need to be computed with high accuracy. A necessary first step for such studies is the selection of appropriate model system which should represent DNA under physiological conditions as accurately as possible while still being computationally feasible. As an illustrative example how the new qsGW implementation can be used effectively in practice, we investigate the dependence of IP and EA of oligomers of Adenine-Thymine (AT) base pairs on the oligomer size.

The calculated charged excitations are shown in Table 2 for different basis sets and fragment sizes between 1 and 4 AT pairs (We refer to these systems as ATx, were x denotes the number of AT base pairs). These systems are shown in Figure 5. For all fragments, we calculated the IPs with the TZ2P and TZ3P basis set with 1d1f, and 2d1f shells of polarization functions for second and third row atoms (and analogously for other atoms). We see, that going from TZ2P to TZ3P only has a small effect on the IPs and EAs, reducing the basis set incompleteness error by only a few 10 meV. These calculations with two rather similar basis sets are necessary to rule out the possibility that a result is simply an artefact of a chosen basis set. Going from TZ3P to QZ6P, the IP of the AT1+B increases by modest 60 meV, while the EA reduces by 180 meV. Based on the TZ3P and QZ6P calculations, we can estimate the QP energies at the CBS limit by extrapolation. Comparing the TZ3P results to the extrapolated ones, we find a basis set limit incompleteness error of 140 meV for the IP and of 420 meV for the EA of AT1. For AT1, we find a similar basis set limit incompleteness error of 80 meV for the IP and of 340 meV for the EA.

Ionization potentials (IPs) and electron affinities (EAS) of DNA fragments consisting of different numbers of adenine-thymine base pairs calculated with different basis sets and contributions of solvent from ΔBLYP calculations. Extra denotes extrapolation to the CBS limit based on TZ3P and QZ6P calculations and numbers in parentheses are obtained by adding the difference between ϵ i C B S − ϵ i T Z 3 P to the result obtained at the TZ3P level. Δsol. has been calculated using COSMO. All values are in eV.

IP EA
Calculation AT1 AT1+B AT2 AT4 AT1 AT1+B AT2 AT4
TZ2P — 7.84 7.34 6.94 — −0.84 −0.65 −0.45
TZ3P 8.47 7.90 7.35 6.97 −0.41 −0.80 −0.63 −0.40
QZ6P 8.50 7.96 — — −0.26 −0.62 — —
Extra 8.55 8.04 (7.49) (7.11) 0.07 -0.38 (−0.21) (0.02)
Δsol −1.82 -0.99 −0.52 −0.01 1.55 — 1.87 1.62
ϵ + Δsol 6.73 7.05 6.97 7.10 1.62 — 1.66 1.64

DNA model systems used in this work.

On standard hardware, calculations on the QZ level are not feasible for AT4 and already for AT2, the QZ calculation is cumbersome. This is not only due to the large number of diffuse AOs which make makes it difficult to exploit distance-based cut-offs (Förster and Visscher, 2020) but also due to the large auxiliary basis sets which are required to make the calculations numerically stable. However, we can estimate the CBS limit based on the differences between the QP energies at the CBS limit and the largest affordable basis set for the larger systems for the smaller fragments. This is justified with the observations made in Förster and Visscher (2020) for G 0 W 0 were we found the basis set incompleteness error on average to decrease with increasing system size but only to a certain extent since basis functions are localised. Based on this assumption, we correct the IPs and EAs of AT2 and AT4 on the TZ3P level by the basis set limit incompleteness error found for AT1+B. (140 and 420 meV, respectively). There is of course a small uncertainty due to the different basis set errors for AT1 and AT1+B. For the extrapolation itself, we assume the error to be rather small for the IP, since the difference between TZ3P and QZ6P are rather small. For the EAs, the error might be larger. Still, we can safely assume, that the basis set errors for AT2 and AT4 are below 100 meV.

The energy required to remove or add an electron from a DNA oligomer in vacuum is strongly size dependent: The vertical IP in vacuum decreases rapidly with increasing oligomer size, with a difference of almost 1 eV between AT1 and AT4. For the EA, a difference of 0.4 eV is found. The IPs of the solvated DNA oligomers, on the other hand, are almost independent of the number of base pairs. When an electron is removed from the oligomer, the surrounding cloud of electrons stabilizes the resulting hole. Increasing the oligomer size thus reduces the IP potential since the hole becomes more and more stabilized. In the aqueous environment, the solvent plays the same role and consequently, the inclusion of water via the COSMO effectively compensates for the effect of the DNA environment. Of course, the comparison is slightly skewed since the DNA environment and the solvent are not treated at the same level of theory. However, there is some evidence that COSMO and other polarizable continuum models are fairly accurate in describing the dielectric screening properties of water (Deglmann and Schenk, 2012).

The IP of AT1+B, AT2, and AT4, all agree within 130 meV. In light of possible basis set errors and errors of the qsGW method itself, the difference is well within the error margin of our method. Only for AT1 we obtain a significantly lower IP, which indicates that the DNA backbone apparently plays an important role in stabilizing ionized DNA oligomers. For the EAs, we arrive at the same conclusion. The differences between the considered systems are even smaller, the aqueous EAs of AT1, AT2 and AT4 being with 1.62, 1.66, and 1.64 eV in excellent agreement. Recently, Pluhařová et al. (2011), Pluhařová et al. (2013), Pluhařová et al. (2015) also concluded that the effect of the DNA environment on the IPs of individual aqueous nucleobases seems to be modest. On the BMK (Boese and Martin, 2004)/6–31G* level of theory, they obtained an IP of 7.24 eV for a fragment of 2 solvated AT base pairs including backbone from the Dickerson dodecamer, but for the isoltaed AT base pair, they obtained and IP of 7.58 eV. The first number is in good agreement with ours, while the second one differs from our result for AT1 by almost 1 eV. However, the difference of only 340 meV between both fragments is of the same order as our difference between the IPs of AT1 and AT2 of 260 meV. Thus, our conclusions regarding the role of the explicit inclusion of the DNA environment on the calculated IPs are very similar.

Finally, we shortly discuss the compute times of the qsGW calculations for the DNA fragments. A detailed timing analysis for the evaluation of the self-energy in ADF has already been performed in Förster and Visscher (2020). The asymptotic scaling of qsGW will be the same as for G 0 W 0: The only additional cubic step is the diagonalization of the Hamiltonian in each iteration. The LU factorization of each of the N ω N fit × N fit matrices in each iteration to calculate the screened interaction (Förster and Visscher, 2020) requires roughly 2 3 N f i t 3 FLOPS, while the dominant step in the single diagonalization of the N bas × N bas matrix in each iteration requires 4 3 N b a s 3 FLOPS. Since we have N bas ≈ 5 × N fit in a typical calculation, the compute time for diagonalization is negligible. Of course, a qsGW calculation requires multiple iterations and is consequently slower than a G 0 W 0 calculation. For the DNA fragments, all calculations required between 6 and 8 iterations to converge. This is considerably faster than the average number of iterations found for GW100, where we have already observed that convergence is typically faster for organic systems. We have set the converge threshold for all calculations in this section to log 10 (ε SCF ) = −8, as opposed to −7 for GW100. However, the increasing sparsity of G 0 (τ → 0+) and G 0 (τ → 0−) with increasing system size is also responsible for this fast convergence.

The largest calculation here is the one for AT4 using the TZ3P basis set. The system has 260 atoms and 1,220 electrons. We used 6,374 MOs and 33,678 auxiliary fit functions. The calculation took 6 iterations to converge and has been performed on 16 cores of a single Dual AMD EPYC 7302@3.0GHz, 2x RTX2070 machine with 256 GB of memory. On average, a single iteration took a little more than 15 h, or 243 core hours.

 4 Conclusion

As opposed to GW calculations with diagonal self-energy, qsGW is a general, parameter-free, and starting point independent method for the calculation of QP energies. While qsGW is known to severely overestimate band gaps and IPs in three-dimensional (3D) materials (Shishkin et al., 2007; Tal et al., 2021) there is evidence that qsGW is more accurate for molecules (Caruso et al., 2016; Kaplan et al., 2016). In canonical implementations, qsGW is usually a magnitude slower than evGW (Gui et al., 2018) and so far, low-order scaling implementations for molecular systems have focused on diagonal approximations to GW (Wilhelm et al., 2018, Wilhelm et al., 2021; Förster and Visscher, 2020; Duchemin and Blase, 2021). To fill this gap, we have presented a low-order scaling implementation of qsGW for molecular systems and demonstrated its accuracy and robustness. In a proof-of-principle application to DNA fragments we have showcased the capabilities of the new implementation for systems of practical interest (Pluhařová et al., 2015; Balanikas et al., 2020). We have shown, that IPs and EAs of the considered DNA fragments in vacuum are strongly size-dependent. Upon taking into account the effect of the aqueous environment, the QP energies become almost independent of the system size. This confirms the results of previous DFT studies. (Pluhařová et al., 2015, Pluhařová et al., 2013). For the largest of the considered fragments with 1,220 electrons, the respective qsGW calculation with more than 6,300 spherical AOs converged within 6 iterations in less than 4 days on a single compute node with 16 cores.

All in all, the herein presented implementation is a necessary stepping stone towards accurate ab initio studies of the spectroscopic properties of large molecules in realistic environments, relevant to organic optoelectronics or biochemistry. To be able to also study optical properties of large systems, it needs to be combined with an implementation of the BSE formalism. Our implementation does not allow to take into account solvent effects directly. In the present work, we have done that via a ΔDFT calculation and obtained consistent results. However, it would be desirable to take into account environmental effects more directly by combining qsGW with COSMO (or a PCM) (Duchemin et al., 2016; Li et al., 2018) and/or molecular mechanics calculations (Tirimbò et al., 2020a; Tirimbò et al., 2020b).

Another issue in practice is the slow convergence of the QP energies to the CBS limit. This is especially true for algorithms like the present one which exploit sparsity in the AO basis. It is encouraging that this convergence is seemingly faster than for qsGW than G 0 W 0. This doesn’t eliminate the need for basis set limit extrapolation, but the extrapolation schemes become more reliable with decreasing basis set error. Basis set errors for large systems can also be accurately estimated based on results for smaller, chemically similar systems, as exemplified in this work.

 Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

 Author Contributions

All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

 Funding

This research received funding from Netherlands Organisation for Scientific Research (NWO) in the framework of the Innovation Fund for Chemistry and from the Ministry of Economic Affairs in the framework of the TKI/PPS-Toeslagregeling (award number 731.017.417).

 Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

 Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

 Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchem.2021.736591/full#supplementary-material

We experimented with different residuals: In the implementation used in this work, we have used r ( n + 1 ) = P ( n + 1 ) + Q ( n + 1 ) − P ( n ) + Q ( n ) .

Here, Q is defined like P, but with the summation spanning the virtual orbital space, Q μ ν = ∑ a N v i r t b μ a b † a ν .

This was based on the intuitive assumption that convergence could be improved by including information about the virtual orbitals in the residual. However, note that S − 1 = 1 2 P + Q , so that we obtain P + Q = 1 2 P + S − 1 . Apart from the factor of 1/2 (Eq. 26) is therefore completely equivalent to this expression. We thank one of the reviewers for pointing this out

 References Adamo C. Barone V. (1999). Toward Reliable Density Functional Methods without Adjustable Parameters: The PBE0 Model. J. Chem. Phys. 110, 6158–6170. 10.1063/1.478522 Almlöf J. Faegri K. Korsell K. (1982). Principles for a DirectSCFapproach toLICAO-MOab-Initiocalculations. J. Comput. Chem. 3, 385–399. 10.1002/jcc.540030314 Baerends E. J. Ellis D. E. Ros P. (1973). Self-Consistent Molecular Hartree-Fock-Slater Calculations I. The Computational Procedure. Chem. Phys. 2, 41–51. 10.1016/0301-0104(73)80059-x Baerends E. Ziegler T. Atkins A. Autschbach J. Baseggio O. Bashford D. (2020). [Dataset]. ADF2020, Locally Modified Development Version. Balanikas E. Banyasz A. Douki T. Baldacchino G. Markovitsi D. (2020). Guanine Radicals Induced in DNA by Low-Energy Photoionization. Acc. Chem. Res. 53, 1511–1519. 10.1021/acs.accounts.0c00245 Balasubramani S. G. Chen G. P. Coriani S. Diedenhofen M. Frank M. S. Franzke Y. J. (2020). TURBOMOLE : Modular Program Suite for Ab Initio Quantum-Chemical and Condensed-Matter Simulations. J. Chem. Phys. 152, 184107. 10.1063/5.0004635 Beach K. S. D. Gooding R. J. Marsiglio F. (2000). Reliable Pade Analytical Continuation Method Based on a High-Accuracy Symbolic Computation Algorithm. Phys. Rev. B 61, 5147–5157. 10.1103/physrevb.61.5147 Becke A. D. (1988). Density-fUnctional Exchange-Energy Approximation with Correct Asymptotic Behavior. Phys. Rev. A. 38, 3098–3100. 10.1063/1.174983510.1103/physreva.38.3098 Berardo E. Kaplan F. Bhaskaran-Nair K. Shelton W. A. Van Setten M. J. Kowalski K. (2017). Benchmarking the Fundamental Electronic Properties of Small TiO2 Nanoclusters by GW and Coupled Cluster Theory Calculations. J. Chem. Theor. Comput. 13, 3814–3828. 10.1021/acs.jctc.7b00538 Bergeron D. Tremblay A. M. (2016). Algorithms for Optimized Maximum Entropy and Diagnostic Tools for Analytic Continuation. Phys. Rev. E 94, 1–25. 10.1103/PhysRevE.94.023303 Blase X. Attaccalite C. (2011). Charge-transfer Excitations in Molecular Donor-Acceptor Complexes within the many-body Bethe-Salpeter Approach. Appl. Phys. Lett. 99, 12–14. 10.1063/1.3655352 Blase X. Attaccalite C. Olevano V. (2011). First-principles GW Calculations for Fullerenes, Porphyrins, Phtalocyanine, and Other Molecules of Interest for Organic Photovoltaic Applications. Phys. Rev. B 83, 1–9. 10.1103/PhysRevB.83.115103 Boese A. D. Martin J. M. (2004). Development of Density Functional for Thermochemical Kinetics. J. Chem. Phys. 121, 3405–3416. 10.1063/1.1774975 Bois J. Körzdörfer T. (2017). Size-Dependence of Nonempirically Tuned DFT Starting Points for G0W0 Applied to π-Conjugated Molecular Chains. J. Chem. Theor. Comput. 13, 4962–4971. 10.1021/acs.jctc.7b00557 Bonačić-Koutecký V. Koutecký J. (1975). General Properties of the Hartree-Fock Problem Demonstrated on the Frontier Orbital Model - I. Relation Among Various Types of Extrema. Theor. Chim. Acta 36, 149–161. 10.1007/BF00572556 Boulanger P. Jacquemin D. Duchemin I. Blase X. (2014). Fast and Accurate Electronic Excitations in Cyanines with the many-body Bethe-Salpeter Approach. J. Chem. Theor. Comput. 10, 1212–1218. 10.1021/ct401101u Bruneval F. Marques M. A. (2013). Benchmarking the Starting Points of the GW Approximation for Molecules. J. Chem. Theor. Comput. 9, 324–329. 10.1021/ct300835h Bruneval F. Rangel T. Hamed S. M. Shao M. Yang C. Neaton J. B. (2016). MOLGW 1: Many-body Perturbation Theory Software for Atoms, Molecules, and Clusters. Comput. Phys. Commun. 208, 149–161. 10.1016/j.cpc.2016.06.019 Bruneval F. Maliyov I. Lapointe C. Marinica M.-C. (2020). Extrapolating Unconverged GW Energies up to the Complete Basis Set Limit with Linear Regression. J. Chem. Theor. Comput. 16, 4399–4407. 10.1021/acs.jctc.0c00433 Bruneval F. (2012). Ionization Energy of Atoms Obtained from GW Self-Energy or from Random Phase Approximation Total Energies. J. Chem. Phys. 136, 194107. 10.1063/1.4718428 Cances E. Le Bris C. (2000). On the Convergence of SCF Algorithms for the Hartree-Fock Equations. ESAIM Math. Model. Numer. Anal. 34, 749–774. 10.1051/m2an:2000102 Caruso F. Rinke P. Ren X. Scheffler M. Rubio A. (2012). Unified Description of Ground and Excited States of Finite Systems: The Self-Consistent GW Approach. Phys. Rev. B 86, 1–5. 10.1103/PhysRevB.86.081102 Caruso F. Rinke P. Ren X. Rubio A. Scheffler M. (2013). Self-consistent GW: All-Electron Implementation with Localized Basis Functions. Phys. Rev. B 88, 1–18. 10.1103/PhysRevB.88.075105 Caruso F. Atalla V. Ren X. Rubio A. Scheffler M. Rinke P. (2014). First-Principles Description of Charge Transfer in Donor-Acceptor Compounds from Self-Consistent Many-Body Perturbation Theory. Phys. Rev. B - Condens. Matter Mater. Phys. 90, 085141. 10.1103/PhysRevB.90.085141 Caruso F. Dauth M. Van Setten M. J. Rinke P. (2016). Benchmark of GW Approaches for the GW100 Test Set. J. Chem. Theor. Comput. 12, 5076–5087. 10.1021/acs.jctc.6b00774 Chu I. H. Trinastic J. P. Wang Y. P. Eguiluz A. G. Kozhevnikov A. Schulthess T. C. (2016). All-electron Self-Consistent GW in the Matsubara-Time Domain: Implementation and Benchmarks of Semiconductors and Insulators. Phys. Rev. B 93, 1–12. 10.1103/PhysRevB.93.125210 Dauth M. Caruso F. Kümmel S. Rinke P. (2016). Piecewise Linearity in the GW Approximation for Accurate Quasiparticle Energy Predictions. Phys. Rev. B 93, 121115(R). 10.1103/PhysRevB.93.121115 Deglmann P. Schenk S. (2012). Thermodynamics of Chemical Reactions with COSMO-RS: The Extreme Case of Charge Separation or Recombination. J. Comput. Chem. 33, 1304–1320. 10.1002/jcc.22961 Doser B. Lambrecht D. S. Kussmann J. Ochsenfeld C. (2009). Linear-scaling Atomic Orbital-Based Second-Order Møller-Plesset Perturbation Theory by Rigorous Integral Screening Criteria. J. Chem. Phys. 130, 064107. 10.1063/1.3072903 Duchemin I. Blase X. (2021). Cubic-Scaling All-Electron GW Calculations with a Separable Density-Fitting Space-Time Approach. J. Chem. Theor. Comput. 17, 2383–2393. 10.1021/acs.jctc.1c00101 Duchemin I. Deutsch T. Blase X. (2012). Short-range to Long-Range Charge-Transfer Excitations in the Zincbacteriochlorin-Bacteriochlorin Complex: A Bethe-Salpeter Study. Phys. Rev. Lett. 109, 1–6. 10.1103/PhysRevLett.109.167801 Duchemin I. Jacquemin D. Blase X. (2016). Combining the GW Formalism with the Polarizable Continuum Model: A State-specific Non-equilibrium Approach. J. Chem. Phys. 144, 164106. 10.1063/1.4946778 Dunlap B. I. Connolly J. W. Sabin J. R. (1979). On Some Approximations in Applications of Xα Theory. J. Chem. Phys. 71, 3396–3402. 10.1063/1.438728 Dyson F. J. (1949). The S Matrix in Quantum Electrodynamics. Phys. Rev. 75, 1736–1755. 10.1103/PhysRev.75.1736 Ernzerhof M. Scuseria G. E. (1999). Assessment of the Perdew–Burke–Ernzerhof Exchange-Correlation Functional. J. Chem. Phys. 110, 5029. 10.1063/1.478401 Faber C. Duchemin I. Deutsch T. Blase X. (2012). Many-body Green’s Function Study of Coumarins for Dye-Sensitized Solar Cells. Phys. Rev. B - Condens. Matter Mater. Phys. 86, 1–7. 10.1103/PhysRevB.86.155315 Fei J. Yeh C. N. Gull E. (2021). Nevanlinna Analytical Continuation. Phys. Rev. Lett. 126, 56402. 10.1103/PhysRevLett.126.056402 Förster A. Visscher L. (2020). Low-Order Scaling G0W0 by Pair Atomic Density Fitting. J. Chem. Theor. Comput. 16, 7381–7399. 10.1021/acs.jctc.0c00693 Förster A. Visscher L. (2021). GW100: A Slater-Type Orbital Perspective. J. Chem. Theory Comput. 17 (8), 5080–5097. 10.1021/acs.jctc.1c00308 Foerster D. Gueddida S. (2021). A Low Resources Space Time Approach to the GW Approximation. Comput. Mater. Sci. 187, 110078. 10.1016/j.commatsci.2020.110078 Foerster D. Koval P. Snchez-Portal D. (2011). An O(N3) Implementation of Hedins GW Approximation for Molecules. J. Chem. Phys. 135, 074105. 10.1063/1.3624731 Fournier R. Wang L. Yazyev O. V. Wu Q. S. (2020). Artificial Neural Network Approach to the Analytic Continuation Problem. Phys. Rev. Lett. 124, 56401. 10.1103/PhysRevLett.124.056401 Fujita T. Noguchi Y. Hoshi T. (2019). Charge-transfer Excited States in the Donor/acceptor Interface from Large-Scale GW Calculations. J. Chem. Phys. 151, 114109. 10.1063/1.5113944 Gallandi L. Körzdörfer T. (2015). Long-range Corrected DFT Meets GW: Vibrationally Resolved Photoelectron Spectra from First Principles. J. Chem. Theor. Comput. 11, 5391–5400. 10.1021/acs.jctc.5b00820 Golze D. Wilhelm J. Van Setten M. J. Rinke P. (2018). Core-Level Binding Energies from GW: An Efficient Full-Frequency Approach within a Localized Basis. J. Chem. Theor. Comput. 14, 4856–4869. 10.1021/acs.jctc.8b00458 Golze D. Dvorak M. Rinke P. (2019). The GW Compendium: A Practical Guide to Theoretical Photoemission Spectroscopy. Front. Chem. 7, 377. 10.3389/fchem.2019.00377 Golze D. Keller L. Rinke P. (2020). Accurate Absolute and Relative Core-Level Binding Energies from GW. J. Phys. Chem. Lett. 11, 1840–1847. 10.1021/acs.jpclett.9b03423 Govoni M. Galli G. (2018). GW100: Comparison of Methods and Accuracy of Results Obtained with the WEST Code. J. Chem. Theor. Comput. 14, 1895–1909. 10.1021/acs.jctc.7b00952 Grimme S. Antony J. Ehrlich S. Krieg H. (2010). A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 132, 154104. 10.1063/1.3382344 Grimme S. Ehrlich S. Goerigk L. (2011). Effect of the Damping Function in Dispersion Corrected Density Functional Theory. J. Comput. Chem. 32, 1456–1465. 10.1002/jcc.21759 Grumet M. Liu P. Kaltak M. Klimeš J. Kresse G. (2018). Beyond the Quasiparticle Approximation: Fully Self-Consistent GW Calculations. Phys. Rev. B 98, 1–9. 10.1103/PhysRevB.98.155143 Gui X. Holzer C. Klopper W. (2018). Accuracy Assessment of GW Starting Points for Calculating Molecular Excitation Energies Using the Bethe-Salpeter Formalism. J. Chem. Theor. Comput. 14, 2127–2136. 10.1021/acs.jctc.8b00014 Gull E. Iskakov S. Krivenko I. Rusakov A. A. Zgid D. (2018). Chebyshev Polynomial Representation of Imaginary-Time Response Functions. Phys. Rev. B 98, 1–10. 10.1103/PhysRevB.98.075127 Hashemi Z. Leppert L. (2021). Assessment of the Ab Initio Bethe - Salpeter Equation Approach for the Low-Lying Excitation Energies of Bacteriochlorophylls and Chlorophylls. J. Phys. Chem. A. 125, 2163–2172. 10.1021/acs.jpca.1c01240 Hedin L. (1965). New Method for Calculating the One-Particle Green’s Function with Application to the Electron-Gas Problem. Phys. Rev. 139, A796. 10.1103/physrev.139.a796 Hung L. Bruneval F. Baishya K. Öğüt S. (2017). Benchmarking the GW Approximation and Bethe-Salpeter Equation for Groups IB and IIB Atoms and Monoxides. J. Chem. Theor. Comput. 13, 2135–2146. 10.1021/acs.jctc.7b00123 Ismail-Beigi S. (2017). Justifying Quasiparticle Self-Consistent Schemes via Gradient Optimization in Baym-Kadanoff Theory. J. Phys. Condens. Matter 29, 385501. 10.1088/1361-648x/aa7803 Kaltak M. Kresse G. (2020). Minimax Isometry Method: A Compressive Sensing Approach for Matsubara Summation in many-body Perturbation Theory. Phys. Rev. B 101, 205145. 10.1103/PhysRevB.101.205145 Kaltak M. Klimeš J. Kresse G. (2014a). Cubic Scaling Algorithm for the Random Phase Approximation: Self-Interstitials and Vacancies in Si. Phys. Rev. B 90, 1–11. 10.1103/PhysRevB.90.054115 Kaltak M. Klimeš J. Kresse G. (2014b). Low Scaling Algorithms for the Random Phase Approximation: Imaginary Time and Laplace Transformations. J. Chem. Theor. Comput. 10, 2498–2507. 10.1021/ct5001268 Kaplan F. Weigend F. Evers F. Van Setten M. J. (2015). Off-diagonal Self-Energy Terms and Partially Self-Consistency in GW Calculations for Single Molecules: Efficient Implementation and Quantitative Effects on Ionization Potentials. J. Chem. Theor. Comput. 11, 5152–5160. 10.1021/acs.jctc.5b00394 Kaplan F. Harding M. E. Seiler C. Weigend F. Evers F. Van Setten M. J. (2016). Quasi-Particle Self-Consistent GW for Molecules. J. Chem. Theor. Comput. 12, 2528–2541. 10.1021/acs.jctc.5b01238 Ke S. H. (2011). All-electron GW Methods Implemented in Molecular Orbital Space: Ionization Energy and Electron Affinity of Conjugated Molecules. Phys. Rev. B 84, 4–7. 10.1103/PhysRevB.84.205415 Kim M. Martyna G. J. Ismail-Beigi S. (2020). Complex-time Shredded Propagator Method for Large-Scale GW Calculations. Phys. Rev. B 101, 1–24. 10.1103/physrevb.101.035139 Klamt A. Jonas V. (1996). Treatment of the Outlying Charge in Continuum Solvation Models. J. Chem. Phys. 105, 9972–9981. 10.1063/1.472829 Klamt A. Schüürmann G. (1993). COSMO: A New Approach to Dielectric Screening in Solvents with Explicit Expressions for the Screening Energy and its Gradient. J. Chem. Soc. Perkin Trans. 2, 799–805. 10.1039/P29930000799 Klamt A. (1995). Conductor-like Screening Model for Real Solvents: A New Approach to the Quantitative Calculation of Solvation Phenomena. J. Phys. Chem. 99, 2224–2235. 10.1021/j100007a062 Klein A. (1961). Perturbation Theory for an Infinite Medium of Fermions. Phys. Rev. 121, 950–956. 10.1103/PhysRev.112.99410.1103/physrev.121.950 Knight J. W. Wang X. Gallandi L. Dolgounitcheva O. Ren X. Ortiz J. V. (2016). Accurate Ionization Potentials and Electron Affinities of Acceptor Molecules III: A Benchmark of GW Methods. J. Chem. Theor. Comput. 12, 615–626. 10.1021/acs.jctc.5b00871 Kotani T. Van Schilfgaarde M. Faleev S. V. (2007). Quasiparticle Self-Consistent GW Method: A Basis for the Independent-Particle Approximation. Phys. Rev. B 76, 1–24. 10.1103/PhysRevB.76.165106 Koutecký J. Bonačić V. (1971). On Convergence Difficulties in the Iterative Hartree-Fock Procedure. J. Chem. Phys. 55, 2408–2413. 10.1063/1.1676424 Koval P. Foerster D. Sánchez-Portal D. (2014). Fully Self-Consistent GW and Quasiparticle Self-Consistent GW for Molecules. Phys. Rev. B - Condens. Matter Mater. Phys. 89, 1–19. 10.1103/PhysRevB.89.155417 Koval P. Ljungberg M. P. Müller M. Sánchez-Portal D. (2019). Toward Efficient GW Calculations Using Numerical Atomic Orbitals: Benchmarking and Application to Molecular Dynamics Simulations. J. Chem. Theor. Comput. 15, 4564–4580. 10.1021/acs.jctc.9b00436 Krause K. Klopper W. (2017). Implementation of the Bethe-Salpeter Equation in the TURBOMOLE Program. J. Comput. Chem. 38, 383–388. 10.1002/jcc.24688 Kudin K. N. Scuseria G. E. Cancès E. (2002). A Black-Box Self-Consistent Field Convergence Algorithm: One Step Closer. J. Chem. Phys. 116, 8255–8261. 10.1063/1.1470195 Kutepov A. Haule K. Savrasov S. Y. Kotliar G. (2012). Electronic Structure of Pu and Am Metals by Self-Consistent Relativistic GW Method. Phys. Rev. B - Condens. Matter Mater. Phys. 85. 10.1103/PhysRevB.85.155129 Kutepov A. L. Oudovenko V. S. Kotliar G. (2017). Linearized Self-Consistent Quasiparticle GW Method: Application to Semiconductors and Simple Metals. Comput. Phys. Commun. 219, 407–414. 10.1016/j.cpc.2017.06.012 Kutepov A. L. (2020). Self-consistent GW Method: O(N) Algorithm for Polarizability and Self Energy. Comput. Phys. Commun. 257, 107502. 10.1016/j.cpc.2020.107502 Lee C. Yang W. Parr R. G. (1988). Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B 37, 785–789. 10.1103/PhysRevB.37.785 Levy R. LeBlanc J. P. Gull E. (2017). Implementation of the Maximum Entropy Method for Analytic Continuation. Comput. Phys. Commun. 215, 149–155. 10.1016/j.cpc.2017.01.018 Li J. D’Avino G. Duchemin I. Beljonne D. Blase X. (2018). Accurate Description of Charged Excitations in Molecular Solids from Embedded many-body Perturbation Theory. Phys. Rev. B 97, 1–13. 10.1103/PhysRevB.97.035108 Liu P. Kaltak M. Klimeš J. Kresse G. (2016). Cubic Scaling GW: Towards Fast Quasiparticle Calculations. Phys. Rev. B 94, 1–13. 10.1103/physrevb.94.165109 Maggio E. Liu P. Van Setten M. J. Kresse G. (2017). GW100: A Plane Wave Perspective for Small Molecules. J. Chem. Theor. Comput. 13, 635–648. 10.1021/acs.jctc.6b01150 Marom N. Moussa J. E. Ren X. Tkatchenko A. Chelikowsky J. R. (2011). Electronic Structure of Dye-Sensitized TiO 2 Clusters from many-body Perturbation Theory. Phys. Rev. B - Condens. Matter Mater. Phys. 84, 1–15. 10.1103/physrevb.84.245115 Marom N. Caruso F. Ren X. Hofmann O. T. Körzdörfer T. Chelikowsky J. R. (2012). Benchmark of GW Methods for Azabenzenes. Phys. Rev. B 86, 1–16. 10.1103/PhysRevB.86.245127 Marom N. Körzdörfer T. Ren X. Tkatchenko A. Chelikowsky J. R. (2014). Size Effects in the Interface Level Alignment of Dye-Sensitized TiO2 Clusters. J. Phys. Chem. Lett. 5, 2395–2401. 10.1021/jz5008356 Martin R. Reining L. Ceperley D. (2016). Interacting Electrons: Theory and Computational Approaches. Cambridge: Cambridge University Press. 10.1017/CBO9781139050807 Miehlich B. Savi A. Stoll H. Preuss H. (1989). Results Obtained with the Correlation Energy Density Functionals of Becke and Lee, Yang and Parr. J. Pediatr. Urol. 157, 200–206. 10.1016/0009-2614(89)87234-3 Migani A. Mowbray D. J. Iacomino A. Zhao J. Petek H. Rubio A. (2013). Level Alignment of a Prototypical Photocatalytic System: Methanol on TiO2(110). J. Am. Chem. Soc. 135, 11429–11432. 10.1021/ja4036994 Migani A. Mowbray D. J. Zhao J. Petek H. Rubio A. (2014). Quasiparticle Level Alignment for Photocatalytic Interfaces. J. Chem. Theor. Comput. 10, 2103–2113. 10.1021/ct500087v Mowbray D. J. Migani A. (2015). Using G0W0 Level Alignment to Identify Catechol’s Structure on TiO2(110). J. Phys. Chem. C 119, 19634–19641. 10.1021/acs.jpcc.5b05392 Neuhauser D. Gao Y. Arntsen C. Karshenas C. Rabani E. Baer R. (2014). Breaking the Theoretical Scaling Limit for Predicting Quasiparticle Energies: The Stochastic GW Approach. Phys. Rev. Lett. 113, 1–5. 10.1103/PhysRevLett.113.076402 Onida G. Reining L. Rubio A. (2002). Electronic Excitations: Density-Functional versus many-body Green’s-function Approaches. Rev. Mod. Phys. 74, 601. 10.1103/RevModPhys.74.601 Otsuki J. Ohzeki M. Shinaoka H. Yoshimi K. (2017). Sparse Modeling Approach to Analytical Continuation of Imaginary-Time Quantum Monte Carlo Data. Phys. Rev. E 95, 4–9. 10.1103/PhysRevE.95.061302 Perdew J. P. Burke K. Ernzerhof M. (1996a). Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 77, 3865–3868. 10.1103/physrevlett.77.3865 Perdew J. P. Ernzerhof M. Burke K. (1996b). [ERRATA] Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 77, 3865–3868. 10.1103/physrevlett.77.3865 Pluhařová E. Jungwirth P. Bradforth S. E. Slavíček P. (2011). Ionization of Purine Tautomers in Nucleobases, Nucleosides, and Nucleotides: From the Gas Phase to the Aqueous Environment. J. Phys. Chem. B 115, 1294–1305. 10.1021/jp110388v Pluhařová E. Schroeder C. Seidel R. Bradforth S. E. Winter B. Faubel M. (2013). Unexpectedly Small Effect of the DNA Environment on Vertical Ionization Energies of Aqueous Nucleobases. J. Phys. Chem. Lett. 4, 3766–3769. 10.1021/jz402106h Pluhařová E. Slavíček P. Jungwirth P. (2015). Modeling Photoionization of Aqueous DNA and its Components. Acc. Chem. Res. 48, 1209–1217. 10.1021/ar500366z Pulay P. (1980). Convergence Acceleration of Iterative Sequences. The Case of Scf Iteration. Chem. Phys. Lett. 73, 393–398. 10.1016/0009-2614(80)80396-4 Pye C. C. Ziegler T. (1999). An Implementation of the Conductor-like Screening Model of Solvation within the Amsterdam Density Functional Package. Theor. Chem. Acc. 101, 396–408. 10.1007/s002140050457 Ren X. Rinke P. Blum V. Wieferink J. Tkatchenko A. Sanfilippo A. (2012). Resolution-of-identity Approach to Hartree-Fock, Hybrid Density Functionals, RPA, MP2 and GW with Numeric Atom-Centered Orbital Basis Functions. New J. Phys. 14, 053020. 10.1088/1367-2630/14/5/053020 Rezaei M. Ögüt S. (2021). Photoelectron Spectra of Early 3d-Transition Metal Dioxide Molecular Anions from GW Calculations. J. Chem. Phys. 154, 094307. 10.1063/5.0042106 Rieger M. M. Steinbeck L. White I. D. Rojas H. N. Godby R. W. (1999). GW Space-Time Method for the Self-Energy of Large Systems. Comput. Phys. Commun. 117, 211–228. 10.1016/S0010-4655(98)00174-X Rojas H. N. Godby R. W. Needs R. J. (1995). Space-Time Method for Ab Initio Calculations ofSelf-Energies and Dielectric Response Functions of Solids. Phys. Rev. Lett. 74, 1827–1831. 10.1103/physrevlett.74.1827 Sakuma R. Miyake T. Aryasetiawan F. (2009). Effective Quasiparticle Hamiltonian Based on Löwdin’s Orthogonalization. Phys. Rev. B - Condens. Matter Mater. Phys. 80, 1–8. 10.1103/PhysRevB.80.235128 Salpeter E. E. Bethe H. A. (1951). A Relativistic Equation for Bound-State Problems. Phys. Rev. 84, 1232–1242. 10.1103/PhysRev.84.1232 Shi B. Weissman S. Bruneval F. Kronik L. Öğüt S. (2018). Photoelectron Spectra of Copper Oxide Cluster Anions from First Principles Methods. J. Chem. Phys. 149, 064306. 10.1063/1.5038744 Shinaoka H. Otsuki J. Ohzeki M. Yoshimi K. (2017). Compressing Green’s Function Using Intermediate Representation between Imaginary-Time and Real-Frequency Domains. Phys. Rev. B 96, 1–8. 10.1103/PhysRevB.96.035147 Shishkin M. Kresse G. (2007). Self-consistent GW Calculations for Semiconductors and Insulators. Phys. Rev. B - Condens. Matter Mater. Phys. 75, 1–9. 10.1103/PhysRevB.75.235102 Shishkin M. Marsman M. Kresse G. (2007). Accurate Quasiparticle Spectra from Self-Consistent GW Calculations with Vertex Corrections. Phys. Rev. Lett. 99, 14–17. 10.1103/PhysRevLett.99.246403 Singh Y. Wang L. W. (2020). Analysis of diagonal G and Subspace W Approximations within Fully Self-Consistent GW Calculations for Bulk Semiconducting Systems. Phys. Rev. B 101, 1–14. 10.1103/PhysRevB.101.235157 Staneke P. O. Groothuis G. Ingemann S. Nibbering N. M. (1995). Formation, Stability and Structure of Radical Anions of Chloroform, Tetrachloromethane and Fluorotrichloromethane in the Gas Phase. Int. J. Mass. Spectrom. Ion Process. 142, 83–93. 10.1016/0168-1176(94)04127-S Strinati G. (1988). Application of the Green’s Functions Method to the Study of the Optical Properties of Semiconductors. La. Riv. Del. Nuovo Cim. Ser. 11, 1–86. 10.1007/BF02725962 Tal A. Chen W. Pasquarello A. (2021). Vertex Function Compliant with the Ward Identity for Quasiparticle Self-Consistent Calculations beyond GW. Phys. Rev. B 103, 161104. 10.1103/PhysRevB.103.L161104 Talman J. D. Shadwick W. F. (1976). Optimized Effective Atomic central Potential. Phys. Rev. A. 14, 36–40. 10.1103/physreva.14.36 Tirimbò G. De Vries X. Weijtens C. H. Bobbert P. A. Neumann T. Coehoorn R. (2020a). Quantitative Predictions of Photoelectron Spectra in Amorphous Molecular Solids from Multiscale Quasiparticle Embedding. Phys. Rev. B 101, 035402. 10.1103/PhysRevB.101.035402 Tirimbò G. Sundaram V. Çaylak O. Scharpach W. Sijen J. Junghans C. (2020b). Excited-state Electronic Structure of Molecules Using many-body Green’s Functions: Quasiparticles and Electron-Hole Excitations with VOTCA-XTP. J. Chem. Phys. 152, 114103. 10.1063/1.5144277 Umari P. Giacomazzi L. De Angelis F. Pastore M. Baroni S. (2013). Energy-level Alignment in Organic Dye-Sensitized TiO2 from GW Calculations. J. Chem. Phys. 139, 014709. 10.1063/1.4809994 Vahtras O. Almlöf J. Feyereisen M. W. (1993). Integral Approximations for LCAO-SCF Calculations. Chem. Phys. Lett. 213, 514–518. 10.1016/0009-2614(93)89151-7 Van Lenthe E. Baerends J. E. (2003). Optimized Slater-type Basis Sets for the Elements 1–118. J. Comput. Chem. 24, 1142–1156. 10.1002/jcc.10255 Van Schilfgaarde M. Kotani T. Faleev S. (2006). Quasiparticle Self-Consistent GW Theory. Phys. Rev. Lett. 96, 1–4. 10.1103/PhysRevLett.96.226402 Van Setten M. J. Weigend F. Evers F. (2013). The GW-Method for Quantum Chemistry Applications: Theory and Implementation. J. Chem. Theor. Comput. 9, 232–246. 10.1021/ct300648t Van Setten M. J. Caruso F. Sharifzadeh S. Ren X. Scheffler M. Liu F. (2015). GW100: Benchmarking G0W0 for Molecular Systems. J. Chem. Theor. Comput. 11, 5665–5687. 10.1021/acs.jctc.5b00453 Van Setten M. J. Costa R. Viñes F. Illas F. (2018). Assessing GW Approaches for Predicting Core Level Binding Energies. J. Chem. Theor. Comput. 14, 877–883. 10.1021/acs.jctc.7b01192 Vidberg H. J. Serene J. W. (1977). Solving the Eliashberg Equations by Means of N-point Padé Approximants. J. Low Temp. Phys. 29, 179–192. 10.1007/bf00655090 Vlček V. Rabani E. Neuhauser D. Baer R. (2017). Stochastic GW Calculations for Molecules. J. Chem. Theor. Comput. 13, 4997–5003. 10.1021/acs.jctc.7b00770 Vlček V. Li W. Baer R. Rabani E. Neuhauser D. (2018). Swift GW beyond 10,000 Electrons Using Sparse Stochastic Compression. Phys. Rev. B 98, 0–10. 10.1103/PhysRevB.98.075107 Ward J. C. (1950). An Identity in Quantum Electrodynamics. Phys. Rev. 78, 182. 10.1103/physrev.78.182 Weng G. Vlcek V. (2021). Efficient Treatment of Molecular Excitations in the Liquid Phase Environment via Stochastic many-body Theory, 1–9. arXiv:2105.14374. Whitten J. L. (1973). Coulombic Potential Energy Integrals and Approximations. J. Chem. Phys. 58, 4496–4501. 10.1063/1.1679012 Wilhelm J. Del Ben M. Hutter J. (2016). GW in the Gaussian and Plane Waves Scheme with Application to Linear Acenes. J. Chem. Theor. Comput. 12, 3623–3635. 10.1021/acs.jctc.6b00380 Wilhelm J. Golze D. Talirz L. Hutter J. Pignedoli C. A. (2018). Toward GW Calculations on Thousands of Atoms. J. Phys. Chem. Lett. 9, 306–312. 10.1021/acs.jpclett.7b02740 Wilhelm J. Seewald P. Golze D. (2021). Low-scaling GW with Benchmark Accuracy and Application to Phosphorene Nanosheets. J. Chem. Theor. Comput. 17, 1662–1677. 10.1021/acs.jctc.0c01282