Low-Order Scaling Quasiparticle Self-Consistent GW for Molecules

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.

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⁴ 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³. 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₀W₀ 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₀W₀ 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₀W₀ 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),

$\sum _r{\mathrm{\Sigma }}_{pr}({\omega }_p)U_{rq}({\omega }_p)=\left[{\omega }_p-{ϵ}_p\right]U_{pq}({\omega }_p).(1)$

We mostly work in a basis of molecular orbitals (MO),

${\varphi }_p^{(n)}(\mathit{r})=\sum _{\mu }{\chi }_{\mu }(\mathit{r})b_{\mu p}^{(n)},(2)$

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,

$\sum _{\nu }{H}_{\mu \nu }^{(0)}{b}_{\nu p}^{(0)}=\sum _{\nu }{S}_{\mu \nu }{b}_{\nu p}^{(0)}{ϵ}_p^{(0)},H_{KS}=T+V_{ext}+V_{Hxc}[P],(3)$

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_{pq}=b_{p\mu }{M}_{\mu \nu }{\left[b^{†}\right]}_{\nu q}.(4)$

In the AO basis, P is given as

$P_{\mu \nu }=2\sum _i^{N_{occ}}{b}_{\mu i}{\left[b^{†}\right]}_{i\nu }.(5)$

We also define the Hamiltonian of the Hartree approximation,

$H_H=H_{KS}-V_{xc}.(6)$

The Green’s function G₀ corresponding to the non-interacting Hamiltonian is diagonal in the MO basis with

${\left[G_0\right]}_{pp}(\omega )={\left[i\omega -{ϵ}_p\right]}^{-1}.(7)$

We can then expand Σ in terms of G₀ as follows (Hedin, 1965; Martin et al., 2016),

$\mathrm{\Sigma }(\omega )=(G_0*W_0)(\omega )+\dots ,(8)$

and in the GWA the expansion is truncated after first order. W₀ is the screened Coulomb interaction, calculated in the bubble approximation (Onida et al., 2002) from G₀ (Hedin, 1965). Without further approximations to Σ, one typically avoids solving (Eq. 1) but instead calculates the interacting Green’s function G by inversion of

${\left[G(\omega )\right]}^{-1}={\left[G_0(\omega )\right]}^{-1}-\mathrm{\Sigma }(\omega ).(9)$

From there one proceeds by building the self-energy (Eq. 8) but replaces G₀ by G, and W₀ 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_{xc}^{qsGW}$. For this mapping, it is convenient to define

$W_0(\omega )=V_c+{\tilde{W}}_0(\omega ),(10)$

with V_c being the bare Coulomb potential. The self-energy can then be decomposed into a static and dynamic part

$\mathrm{\Sigma }(\omega )={\mathrm{\Sigma }}_x+(G_0*{\tilde{W}}_0)(\omega )={\mathrm{\Sigma }}_x+{\mathrm{\Sigma }}_c(\omega ).(11)$

Σ_x is the Fock exchange potential, $V_x^{qsGW}={\mathrm{\Sigma }}_x$, and following Kotani et al. (2007), the correlation part of V_xc is obtained from Σ by taking one of the real symmetric definitions

${\left[V_c^{qsGW}\right]}_{pq}=\frac{1}{2}\left[\mathrm{R}\mathrm{e}{\left[{\mathrm{\Sigma }}_c\right]}_{pq}({ϵ}_p)+\mathrm{R}\mathrm{e}{\left[{\mathrm{\Sigma }}_c\right]}_{pq}({ϵ}_q)\right],(12)$

or

${\left[V_c^{qsGW}\right]}_{pq}={\delta }_{pq}\mathrm{R}\mathrm{e}{\left[{\mathrm{\Sigma }}_c\right]}_{pq}({ϵ}_p)+(1-{\delta }_{pq})\mathrm{R}\mathrm{e}{\left[{\mathrm{\Sigma }}_c\right]}_{pq}(\omega =0).(13)$

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₀ (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

$\sum _r{H}_{pr}^{qsG{W}^{(n+1)}}{U}_{rq}^{(n+1)}={\omega }_p^{(n+1)}{U}_{pq}^{(n+1)},(14)$

with

$H^{qsG{W}^{(n+1)}}=H_H+\mathrm{\Delta }{V}_H^{(n+1)}+V_{xc}^{qsG{W}^{(n+1)}}(15)$

and

$V_{xc}^{qsG{W}^{(n+1)}}=V_x[P^{(n)}]+V_c^{qsGW}[G_0^{(n)}].(16)$

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 $\left\{{\varphi }_i^{(n)}\right\}$ and unless self-consistency has been reached, U⁽ⁿ⁾ will not be unity and defines a rotation of the molecular orbitals. We now set

$\begin{array}{cc}\hfill b_{\mu p}^{(n+1)}=& \sum _q{b}_{\mu q}^{(n)}{U}_{qp}^{(n+1)}\hfill \\ \hfill {ϵ}_p^{(n+1)}=& {\omega }_p^{(n+1)}\forall p\hfill \end{array}(17)$

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⁽ⁿ⁺¹⁾ is then evaluated from (Eq. 5) and the change in the Hartree-potential is calculated as

$\mathrm{\Delta }{V}_H^{(n+1)}=V_H[\mathrm{\Delta }{P}^{(n+1)}],(18)$

with

$\mathrm{\Delta }{P}^{(n+1)}=P^{(n+1)}-P^{(n)}.(19)$

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₀ 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 $\mathcal{W}={\left\{i{\omega }_{\beta }\right\}}_{\beta =1,N_{\omega }}$ 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 $\mathcal{A}\subset \mathbb{C}$ and coincides with Σ in $\mathcal{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 $\mathcal{A}=\mathbb{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₀W₀ 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^{qsGW}$. At self-consistency, the off-diagonal elements of $V_c^{qsGW}$ need to be zero: In the n + 1 the iteration, $V_c^{qsGW}$ is expressed in the basis which diagonalizes the operator defined in (Eq. 15) in the nth iteration. At self-consistency b⁽ⁿ⁺¹⁾ = b⁽ⁿ⁾, which will not be the case when the off-diagonal elements of $V_c^{qsGW}$ will be different from zero. In our present implementation, we set all values with magnitude smaller than 1e⁻⁶ 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

${\left\{G_0^{(m)}\right\}}_{0\le m\le n+1}\to {\tilde{H}}^{qsG{W}^{n+1}}\to {ϵ}^{(n+1)},b^{(n+1)}(20)$

is clearly a better option. A practical way to implement this is to replace (Eq. 14) by

$\sum _r{\tilde{H}}_{pr}^{qsG{W}^{(n+1)}}{U}_{rq}^{(n+1)}={\omega }_p^{(n+1)}{U}_{pq}^{(n+1)},(21)$

with

${\tilde{H}}^{qsG{W}^{(n+1)}}=\sum _{m=n-n_0}^{n+1}{\alpha }_m{H}^{qsG{W}^{(m)}},(22)$

where

$\sum _{m=n-n_0}^n{\alpha }_m=1,(23)$

needs to be fulfilled and n₀ 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)}=\sum _{m=n-n_0}^n{\alpha }_m{r}^{(m)},(24)$

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₀ uniquely determines H^qsGW, we would ideally define

$r^{(n+1)}=G_0^{(n+1)}-G_0^{(n)},(25)$

however, storage (or recalculation) of this quantity for n₀ iterations is inefficient. Therefore, one can use

$r^{(n+1)}=P^{(n+1)}-P^{(n)},(26)$

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).¹

Technically, in the nth iteration we solve (Eq. 14) and evaluate the corresponding b⁽ⁿ⁾ from which we calculate P⁽ⁿ⁾ and Q⁽ⁿ⁾. We check for convergence by evaluating the Frobenius norm of the residual (Eq. 26),

$N_F=\frac{1}{N_{MO}^2}\sqrt{\sum _{\mu \nu }{\left[r_{\mu \nu }^{(n+1)}\right]}^2},(27)$

and terminate the SCF as soon as N_F < ε_SCF for two subsequent iterations. As we will show later on, ε_SCF = 1e⁻⁷ 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⁽ⁿ⁺¹⁾ and $H^{qsG{W}^{(n+1)}}$ and determine the expansion coefficients α_m using the DIIS method, setting n₀ = 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, Na₂, Na₄, and Na₆, 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 $\mathrm{\Delta }{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 qsGW

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 ₁₀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 ₁₀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.

FIGURE 1

FIGURE 1. Convergence of the qsGW SCF for Methane for different initial guesses and constructions of the correlation potential.  log ₁₀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⁻⁷ and restrict all calculations to a maximum of 30 iterations.

Independent of the starting point, we could not reach convergence for Mgo, BeO, BN, Cu₂, 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.

FIGURE 2

FIGURE 2. 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.

FIGURE 3

FIGURE 3. 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.

FIGURE 4

FIGURE 4. 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₀W₀ counterparts. Going from TZ3P to QZ6P, the basis set incompleteness error reduces by 80 meV on average, while for G₀W₀@PBE, we found an average reduction of 130 meV (Förster and Visscher, 2021). Self-consistent approaches might converge faster than G₀W₀ - Caruso et al. have already observed that scGW converges faster to the CBS limit than G₀W₀ (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).

TABLE 1

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

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.

TABLE 2

TABLE 2. 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^{CBS}-{ϵ}_i^{TZ3P}$ to the result obtained at the TZ3P level. Δsol. has been calculated using COSMO. All values are in eV.

FIGURE 5

FIGURE 5. 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₀W₀ 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₀W₀: 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 $\frac{2}{3}{N}_{fit}^3$ FLOPS, while the dominant step in the single diagonalization of the N_bas × N_bas matrix in each iteration requires $\frac{4}{3}{N}_{bas}^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₀W₀ 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₁₀ (ε_SCF) = −8, as opposed to −7 for GW100. However, the increasing sparsity of G₀ (τ → 0⁺) and G₀ (τ → 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₀W₀. 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

Footnotes

¹We experimented with different residuals: In the implementation used in this work, we have used

$r^{(n+1)}=\left[P^{(n+1)}+Q^{(n+1)}\right]-\left[P^{(n)}+Q^{(n)}\right].$

Here, Q is defined like P, but with the summation spanning the virtual orbital space,

$Q_{\mu \nu }=\sum _a^{N_{virt}}{b}_{\mu a}{\left[b^{†}\right]}_{a\nu }.$

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}=\frac{1}{2}P+Q$, so that we obtain $P+Q=\frac{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., and Barone, V. (1999). Toward Reliable Density Functional Methods without Adjustable Parameters: The PBE0 Model. J. Chem. Phys. 110, 6158–6170. doi:10.1063/1.478522

CrossRef Full Text | Google Scholar

Almlöf, J., Faegri, K., and Korsell, K. (1982). Principles for a DirectSCFapproach toLICAO-MOab-Initiocalculations. J. Comput. Chem. 3, 385–399. doi:10.1002/jcc.540030314

CrossRef Full Text | Google Scholar

Baerends, E. J., Ellis, D. E., and Ros, P. (1973). Self-Consistent Molecular Hartree-Fock-Slater Calculations I. The Computational Procedure. Chem. Phys. 2, 41–51. doi:10.1016/0301-0104(73)80059-x

CrossRef Full Text | Google Scholar

Baerends, E., Ziegler, T., Atkins, A., Autschbach, J., Baseggio, O., Bashford, D., et al. (2020). [Dataset]. ADF2020, Locally Modified Development Version.

Balanikas, E., Banyasz, A., Douki, T., Baldacchino, G., and Markovitsi, D. (2020). Guanine Radicals Induced in DNA by Low-Energy Photoionization. Acc. Chem. Res. 53, 1511–1519. doi:10.1021/acs.accounts.0c00245

PubMed Abstract | CrossRef Full Text | Google Scholar

Balasubramani, S. G., Chen, G. P., Coriani, S., Diedenhofen, M., Frank, M. S., Franzke, Y. J., et al. (2020). TURBOMOLE : Modular Program Suite for Ab Initio Quantum-Chemical and Condensed-Matter Simulations. J. Chem. Phys. 152, 184107. doi:10.1063/5.0004635

PubMed Abstract | CrossRef Full Text | Google Scholar

Beach, K. S. D., Gooding, R. J., and Marsiglio, F. (2000). Reliable Pade Analytical Continuation Method Based on a High-Accuracy Symbolic Computation Algorithm. Phys. Rev. B 61, 5147–5157. doi:10.1103/physrevb.61.5147

CrossRef Full Text | Google Scholar

Becke, A. D. (1988). Density-fUnctional Exchange-Energy Approximation with Correct Asymptotic Behavior. Phys. Rev. A. 38, 3098–3100. doi:10.1063/1.174983510.1103/physreva.38.3098

PubMed Abstract | CrossRef Full Text | Google Scholar

Berardo, E., Kaplan, F., Bhaskaran-Nair, K., Shelton, W. A., Van Setten, M. J., Kowalski, K., et al. (2017). Benchmarking the Fundamental Electronic Properties of Small TiO2 Nanoclusters by GW and Coupled Cluster Theory Calculations. J. Chem. Theor. Comput. 13, 3814–3828. doi:10.1021/acs.jctc.7b00538

PubMed Abstract | CrossRef Full Text | Google Scholar

Bergeron, D., and Tremblay, A. M. (2016). Algorithms for Optimized Maximum Entropy and Diagnostic Tools for Analytic Continuation. Phys. Rev. E 94, 1–25. doi:10.1103/PhysRevE.94.023303

CrossRef Full Text | Google Scholar

Blase, X., and Attaccalite, C. (2011). Charge-transfer Excitations in Molecular Donor-Acceptor Complexes within the many-body Bethe-Salpeter Approach. Appl. Phys. Lett. 99, 12–14. doi:10.1063/1.3655352

CrossRef Full Text | Google Scholar

Blase, X., Attaccalite, C., and 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. doi:10.1103/PhysRevB.83.115103

CrossRef Full Text | Google Scholar

Boese, A. D., and Martin, J. M. (2004). Development of Density Functional for Thermochemical Kinetics. J. Chem. Phys. 121, 3405–3416. doi:10.1063/1.1774975

PubMed Abstract | CrossRef Full Text | Google Scholar

Bois, J., and 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. doi:10.1021/acs.jctc.7b00557

PubMed Abstract | CrossRef Full Text | Google Scholar

Bonačić-Koutecký, V., and 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. doi:10.1007/BF00572556

CrossRef Full Text | Google Scholar

Boulanger, P., Jacquemin, D., Duchemin, I., and Blase, X. (2014). Fast and Accurate Electronic Excitations in Cyanines with the many-body Bethe-Salpeter Approach. J. Chem. Theor. Comput. 10, 1212–1218. doi:10.1021/ct401101u

PubMed Abstract | CrossRef Full Text | Google Scholar

Bruneval, F., and Marques, M. A. (2013). Benchmarking the Starting Points of the GW Approximation for Molecules. J. Chem. Theor. Comput. 9, 324–329. doi:10.1021/ct300835h

PubMed Abstract | CrossRef Full Text | Google Scholar

Bruneval, F., Rangel, T., Hamed, S. M., Shao, M., Yang, C., and Neaton, J. B. (2016). MOLGW 1: Many-body Perturbation Theory Software for Atoms, Molecules, and Clusters. Comput. Phys. Commun. 208, 149–161. doi:10.1016/j.cpc.2016.06.019

CrossRef Full Text | Google Scholar

Bruneval, F., Maliyov, I., Lapointe, C., and 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. doi:10.1021/acs.jctc.0c00433

PubMed Abstract | CrossRef Full Text | Google Scholar

Bruneval, F. (2012). Ionization Energy of Atoms Obtained from GW Self-Energy or from Random Phase Approximation Total Energies. J. Chem. Phys. 136, 194107. doi:10.1063/1.4718428

PubMed Abstract | CrossRef Full Text | Google Scholar

Cances, E., and Le Bris, C. (2000). On the Convergence of SCF Algorithms for the Hartree-Fock Equations. ESAIM Math. Model. Numer. Anal. 34, 749–774. doi:10.1051/m2an:2000102

CrossRef Full Text | Google Scholar

Caruso, F., Rinke, P., Ren, X., Scheffler, M., and Rubio, A. (2012). Unified Description of Ground and Excited States of Finite Systems: The Self-Consistent GW Approach. Phys. Rev. B 86, 1–5. doi:10.1103/PhysRevB.86.081102

CrossRef Full Text | Google Scholar

Caruso, F., Rinke, P., Ren, X., Rubio, A., and Scheffler, M. (2013). Self-consistent GW: All-Electron Implementation with Localized Basis Functions. Phys. Rev. B 88, 1–18. doi:10.1103/PhysRevB.88.075105

CrossRef Full Text | Google Scholar

Caruso, F., Atalla, V., Ren, X., Rubio, A., Scheffler, M., and 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. doi:10.1103/PhysRevB.90.085141

CrossRef Full Text | Google Scholar

Caruso, F., Dauth, M., Van Setten, M. J., and Rinke, P. (2016). Benchmark of GW Approaches for the GW100 Test Set. J. Chem. Theor. Comput. 12, 5076–5087. doi:10.1021/acs.jctc.6b00774

PubMed Abstract | CrossRef Full Text | Google Scholar

Chu, I. H., Trinastic, J. P., Wang, Y. P., Eguiluz, A. G., Kozhevnikov, A., Schulthess, T. C., et al. (2016). All-electron Self-Consistent GW in the Matsubara-Time Domain: Implementation and Benchmarks of Semiconductors and Insulators. Phys. Rev. B 93, 1–12. doi:10.1103/PhysRevB.93.125210

CrossRef Full Text | Google Scholar

Dauth, M., Caruso, F., Kümmel, S., and Rinke, P. (2016). Piecewise Linearity in the GW Approximation for Accurate Quasiparticle Energy Predictions. Phys. Rev. B 93, 121115(R). doi:10.1103/PhysRevB.93.121115

CrossRef Full Text | Google Scholar

Deglmann, P., and Schenk, S. (2012). Thermodynamics of Chemical Reactions with COSMO-RS: The Extreme Case of Charge Separation or Recombination. J. Comput. Chem. 33, 1304–1320. doi:10.1002/jcc.22961

PubMed Abstract | CrossRef Full Text | Google Scholar

Doser, B., Lambrecht, D. S., Kussmann, J., and 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. doi:10.1063/1.3072903

PubMed Abstract | CrossRef Full Text | Google Scholar

Duchemin, I., and Blase, X. (2021). Cubic-Scaling All-Electron GW Calculations with a Separable Density-Fitting Space-Time Approach. J. Chem. Theor. Comput. 17, 2383–2393. doi:10.1021/acs.jctc.1c00101

CrossRef Full Text | Google Scholar

Duchemin, I., Deutsch, T., and 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. doi:10.1103/PhysRevLett.109.167801

PubMed Abstract | CrossRef Full Text | Google Scholar

Duchemin, I., Jacquemin, D., and Blase, X. (2016). Combining the GW Formalism with the Polarizable Continuum Model: A State-specific Non-equilibrium Approach. J. Chem. Phys. 144, 164106. doi:10.1063/1.4946778

PubMed Abstract | CrossRef Full Text | Google Scholar

Dunlap, B. I., Connolly, J. W., and Sabin, J. R. (1979). On Some Approximations in Applications of Xα Theory. J. Chem. Phys. 71, 3396–3402. doi:10.1063/1.438728

CrossRef Full Text | Google Scholar

Dyson, F. J. (1949). The S Matrix in Quantum Electrodynamics. Phys. Rev. 75, 1736–1755. doi:10.1103/PhysRev.75.1736

CrossRef Full Text | Google Scholar

Ernzerhof, M., and Scuseria, G. E. (1999). Assessment of the Perdew–Burke–Ernzerhof Exchange-Correlation Functional. J. Chem. Phys. 110, 5029. doi:10.1063/1.478401

CrossRef Full Text | Google Scholar

Faber, C., Duchemin, I., Deutsch, T., and 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. doi:10.1103/PhysRevB.86.155315

CrossRef Full Text | Google Scholar

Fei, J., Yeh, C. N., and Gull, E. (2021). Nevanlinna Analytical Continuation. Phys. Rev. Lett. 126, 56402. doi:10.1103/PhysRevLett.126.056402

CrossRef Full Text | Google Scholar

Förster, A., and Visscher, L. (2020). Low-Order Scaling G0W0 by Pair Atomic Density Fitting. J. Chem. Theor. Comput. 16, 7381–7399. doi:10.1021/acs.jctc.0c00693

CrossRef Full Text | Google Scholar

Förster, A., and Visscher, L. (2021). GW100: A Slater-Type Orbital Perspective. J. Chem. Theory Comput. 17 (8), 5080–5097. doi:10.1021/acs.jctc.1c00308

Foerster, D., and Gueddida, S. (2021). A Low Resources Space Time Approach to the GW Approximation. Comput. Mater. Sci. 187, 110078. doi:10.1016/j.commatsci.2020.110078

CrossRef Full Text | Google Scholar

Foerster, D., Koval, P., and Snchez-Portal, D. (2011). An O(N3) Implementation of Hedins GW Approximation for Molecules. J. Chem. Phys. 135, 074105. doi:10.1063/1.3624731

PubMed Abstract | CrossRef Full Text | Google Scholar

Fournier, R., Wang, L., Yazyev, O. V., and Wu, Q. S. (2020). Artificial Neural Network Approach to the Analytic Continuation Problem. Phys. Rev. Lett. 124, 56401. doi:10.1103/PhysRevLett.124.056401

PubMed Abstract | CrossRef Full Text | Google Scholar

Fujita, T., Noguchi, Y., and Hoshi, T. (2019). Charge-transfer Excited States in the Donor/acceptor Interface from Large-Scale GW Calculations. J. Chem. Phys. 151, 114109. doi:10.1063/1.5113944

PubMed Abstract | CrossRef Full Text | Google Scholar

Gallandi, L., and 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. doi:10.1021/acs.jctc.5b00820

PubMed Abstract | CrossRef Full Text | Google Scholar

Golze, D., Wilhelm, J., Van Setten, M. J., and 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. doi:10.1021/acs.jctc.8b00458

PubMed Abstract | CrossRef Full Text | Google Scholar

Golze, D., Dvorak, M., and Rinke, P. (2019). The GW Compendium: A Practical Guide to Theoretical Photoemission Spectroscopy. Front. Chem. 7, 377. doi:10.3389/fchem.2019.00377

PubMed Abstract | CrossRef Full Text | Google Scholar

Golze, D., Keller, L., and Rinke, P. (2020). Accurate Absolute and Relative Core-Level Binding Energies from GW. J. Phys. Chem. Lett. 11, 1840–1847. doi:10.1021/acs.jpclett.9b03423

PubMed Abstract | CrossRef Full Text | Google Scholar

Govoni, M., and Galli, G. (2018). GW100: Comparison of Methods and Accuracy of Results Obtained with the WEST Code. J. Chem. Theor. Comput. 14, 1895–1909. doi:10.1021/acs.jctc.7b00952

PubMed Abstract | CrossRef Full Text | Google Scholar

Grimme, S., Antony, J., Ehrlich, S., and 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. doi:10.1063/1.3382344

PubMed Abstract | CrossRef Full Text | Google Scholar

Grimme, S., Ehrlich, S., and Goerigk, L. (2011). Effect of the Damping Function in Dispersion Corrected Density Functional Theory. J. Comput. Chem. 32, 1456–1465. doi:10.1002/jcc.21759

PubMed Abstract | CrossRef Full Text | Google Scholar

Grumet, M., Liu, P., Kaltak, M., Klimeš, J., and Kresse, G. (2018). Beyond the Quasiparticle Approximation: Fully Self-Consistent GW Calculations. Phys. Rev. B 98, 1–9. doi:10.1103/PhysRevB.98.155143

CrossRef Full Text | Google Scholar

Gui, X., Holzer, C., and 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. doi:10.1021/acs.jctc.8b00014

PubMed Abstract | CrossRef Full Text | Google Scholar

Gull, E., Iskakov, S., Krivenko, I., Rusakov, A. A., and Zgid, D. (2018). Chebyshev Polynomial Representation of Imaginary-Time Response Functions. Phys. Rev. B 98, 1–10. doi:10.1103/PhysRevB.98.075127

CrossRef Full Text | Google Scholar

Hashemi, Z., and 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. doi:10.1021/acs.jpca.1c01240

PubMed Abstract | CrossRef Full Text | Google Scholar

Hedin, L. (1965). New Method for Calculating the One-Particle Green’s Function with Application to the Electron-Gas Problem. Phys. Rev. 139, A796. doi:10.1103/physrev.139.a796

CrossRef Full Text | Google Scholar

Hung, L., Bruneval, F., Baishya, K., and Öğü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. doi:10.1021/acs.jctc.7b00123

PubMed Abstract | CrossRef Full Text | Google Scholar

Ismail-Beigi, S. (2017). Justifying Quasiparticle Self-Consistent Schemes via Gradient Optimization in Baym-Kadanoff Theory. J. Phys. Condens. Matter 29, 385501. doi:10.1088/1361-648x/aa7803

PubMed Abstract | CrossRef Full Text | Google Scholar

Kaltak, M., and Kresse, G. (2020). Minimax Isometry Method: A Compressive Sensing Approach for Matsubara Summation in many-body Perturbation Theory. Phys. Rev. B 101, 205145. doi:10.1103/PhysRevB.101.205145

CrossRef Full Text | Google Scholar

Kaltak, M., Klimeš, J., and Kresse, G. (2014a). Cubic Scaling Algorithm for the Random Phase Approximation: Self-Interstitials and Vacancies in Si. Phys. Rev. B 90, 1–11. doi:10.1103/PhysRevB.90.054115

CrossRef Full Text | Google Scholar

Kaltak, M., Klimeš, J., and Kresse, G. (2014b). Low Scaling Algorithms for the Random Phase Approximation: Imaginary Time and Laplace Transformations. J. Chem. Theor. Comput. 10, 2498–2507. doi:10.1021/ct5001268

PubMed Abstract | CrossRef Full Text | Google Scholar

Kaplan, F., Weigend, F., Evers, F., and 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. doi:10.1021/acs.jctc.5b00394

PubMed Abstract | CrossRef Full Text | Google Scholar

Kaplan, F., Harding, M. E., Seiler, C., Weigend, F., Evers, F., and Van Setten, M. J. (2016). Quasi-Particle Self-Consistent GW for Molecules. J. Chem. Theor. Comput. 12, 2528–2541. doi:10.1021/acs.jctc.5b01238

PubMed Abstract | CrossRef Full Text | Google Scholar

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. doi:10.1103/PhysRevB.84.205415

CrossRef Full Text | Google Scholar

Kim, M., Martyna, G. J., and Ismail-Beigi, S. (2020). Complex-time Shredded Propagator Method for Large-Scale GW Calculations. Phys. Rev. B 101, 1–24. doi:10.1103/physrevb.101.035139

CrossRef Full Text | Google Scholar

Klamt, A., and Jonas, V. (1996). Treatment of the Outlying Charge in Continuum Solvation Models. J. Chem. Phys. 105, 9972–9981. doi:10.1063/1.472829

CrossRef Full Text | Google Scholar

Klamt, A., and 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. doi:10.1039/P29930000799

CrossRef Full Text | Google Scholar

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. doi:10.1021/j100007a062

CrossRef Full Text | Google Scholar

Klein, A. (1961). Perturbation Theory for an Infinite Medium of Fermions. Phys. Rev. 121, 950–956. doi:10.1103/PhysRev.112.99410.1103/physrev.121.950

CrossRef Full Text | Google Scholar

Knight, J. W., Wang, X., Gallandi, L., Dolgounitcheva, O., Ren, X., Ortiz, J. V., et al. (2016). Accurate Ionization Potentials and Electron Affinities of Acceptor Molecules III: A Benchmark of GW Methods. J. Chem. Theor. Comput. 12, 615–626. doi:10.1021/acs.jctc.5b00871

PubMed Abstract | CrossRef Full Text | Google Scholar

Kotani, T., Van Schilfgaarde, M., and Faleev, S. V. (2007). Quasiparticle Self-Consistent GW Method: A Basis for the Independent-Particle Approximation. Phys. Rev. B 76, 1–24. doi:10.1103/PhysRevB.76.165106

CrossRef Full Text | Google Scholar

Koutecký, J., and Bonačić, V. (1971). On Convergence Difficulties in the Iterative Hartree-Fock Procedure. J. Chem. Phys. 55, 2408–2413. doi:10.1063/1.1676424

CrossRef Full Text | Google Scholar

Koval, P., Foerster, D., and 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. doi:10.1103/PhysRevB.89.155417

CrossRef Full Text | Google Scholar

Koval, P., Ljungberg, M. P., Müller, M., and 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. doi:10.1021/acs.jctc.9b00436

PubMed Abstract | CrossRef Full Text | Google Scholar

Krause, K., and Klopper, W. (2017). Implementation of the Bethe-Salpeter Equation in the TURBOMOLE Program. J. Comput. Chem. 38, 383–388. doi:10.1002/jcc.24688

PubMed Abstract | CrossRef Full Text | Google Scholar

Kudin, K. N., Scuseria, G. E., and Cancès, E. (2002). A Black-Box Self-Consistent Field Convergence Algorithm: One Step Closer. J. Chem. Phys. 116, 8255–8261. doi:10.1063/1.1470195

CrossRef Full Text | Google Scholar

Kutepov, A., Haule, K., Savrasov, S. Y., and Kotliar, G. (2012). Electronic Structure of Pu and Am Metals by Self-Consistent Relativistic GW Method. Phys. Rev. B - Condens. Matter Mater. Phys. 85. doi:10.1103/PhysRevB.85.155129

CrossRef Full Text | Google Scholar

Kutepov, A. L., Oudovenko, V. S., and Kotliar, G. (2017). Linearized Self-Consistent Quasiparticle GW Method: Application to Semiconductors and Simple Metals. Comput. Phys. Commun. 219, 407–414. doi:10.1016/j.cpc.2017.06.012

CrossRef Full Text | Google Scholar

Kutepov, A. L. (2020). Self-consistent GW Method: O(N) Algorithm for Polarizability and Self Energy. Comput. Phys. Commun. 257, 107502. doi:10.1016/j.cpc.2020.107502

CrossRef Full Text | Google Scholar

Lee, C., Yang, W., and 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. doi:10.1103/PhysRevB.37.785

PubMed Abstract | CrossRef Full Text | Google Scholar

Levy, R., LeBlanc, J. P., and Gull, E. (2017). Implementation of the Maximum Entropy Method for Analytic Continuation. Comput. Phys. Commun. 215, 149–155. doi:10.1016/j.cpc.2017.01.018

CrossRef Full Text | Google Scholar

Li, J., D’Avino, G., Duchemin, I., Beljonne, D., and Blase, X. (2018). Accurate Description of Charged Excitations in Molecular Solids from Embedded many-body Perturbation Theory. Phys. Rev. B 97, 1–13. doi:10.1103/PhysRevB.97.035108

CrossRef Full Text | Google Scholar

Liu, P., Kaltak, M., Klimeš, J., and Kresse, G. (2016). Cubic Scaling GW: Towards Fast Quasiparticle Calculations. Phys. Rev. B 94, 1–13. doi:10.1103/physrevb.94.165109

CrossRef Full Text | Google Scholar

Maggio, E., Liu, P., Van Setten, M. J., and Kresse, G. (2017). GW100: A Plane Wave Perspective for Small Molecules. J. Chem. Theor. Comput. 13, 635–648. doi:10.1021/acs.jctc.6b01150

PubMed Abstract | CrossRef Full Text | Google Scholar

Marom, N., Moussa, J. E., Ren, X., Tkatchenko, A., and 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. doi:10.1103/physrevb.84.245115

CrossRef Full Text | Google Scholar

Marom, N., Caruso, F., Ren, X., Hofmann, O. T., Körzdörfer, T., Chelikowsky, J. R., et al. (2012). Benchmark of GW Methods for Azabenzenes. Phys. Rev. B 86, 1–16. doi:10.1103/PhysRevB.86.245127

CrossRef Full Text | Google Scholar

Marom, N., Körzdörfer, T., Ren, X., Tkatchenko, A., and Chelikowsky, J. R. (2014). Size Effects in the Interface Level Alignment of Dye-Sensitized TiO2 Clusters. J. Phys. Chem. Lett. 5, 2395–2401. doi:10.1021/jz5008356

PubMed Abstract | CrossRef Full Text | Google Scholar

Martin, R., Reining, L., and Ceperley, D. (2016). Interacting Electrons: Theory and Computational Approaches. Cambridge: Cambridge University Press. doi:10.1017/CBO9781139050807

Miehlich, B., Savi, A., Stoll, H., and Preuss, H. (1989). Results Obtained with the Correlation Energy Density Functionals of Becke and Lee, Yang and Parr. J. Pediatr. Urol. 157, 200–206. doi:10.1016/0009-2614(89)87234-3

CrossRef Full Text | Google Scholar

Migani, A., Mowbray, D. J., Iacomino, A., Zhao, J., Petek, H., and Rubio, A. (2013). Level Alignment of a Prototypical Photocatalytic System: Methanol on TiO2(110). J. Am. Chem. Soc. 135, 11429–11432. doi:10.1021/ja4036994

PubMed Abstract | CrossRef Full Text | Google Scholar

Migani, A., Mowbray, D. J., Zhao, J., Petek, H., and Rubio, A. (2014). Quasiparticle Level Alignment for Photocatalytic Interfaces. J. Chem. Theor. Comput. 10, 2103–2113. doi:10.1021/ct500087v

PubMed Abstract | CrossRef Full Text | Google Scholar

Mowbray, D. J., and Migani, A. (2015). Using G0W0 Level Alignment to Identify Catechol’s Structure on TiO2(110). J. Phys. Chem. C 119, 19634–19641. doi:10.1021/acs.jpcc.5b05392

CrossRef Full Text | Google Scholar

Neuhauser, D., Gao, Y., Arntsen, C., Karshenas, C., Rabani, E., and Baer, R. (2014). Breaking the Theoretical Scaling Limit for Predicting Quasiparticle Energies: The Stochastic GW Approach. Phys. Rev. Lett. 113, 1–5. doi:10.1103/PhysRevLett.113.076402

CrossRef Full Text | Google Scholar

Onida, G., Reining, L., and Rubio, A. (2002). Electronic Excitations: Density-Functional versus many-body Green’s-function Approaches. Rev. Mod. Phys. 74, 601. doi:10.1103/RevModPhys.74.601

CrossRef Full Text | Google Scholar

Otsuki, J., Ohzeki, M., Shinaoka, H., and Yoshimi, K. (2017). Sparse Modeling Approach to Analytical Continuation of Imaginary-Time Quantum Monte Carlo Data. Phys. Rev. E 95, 4–9. doi:10.1103/PhysRevE.95.061302

PubMed Abstract | CrossRef Full Text | Google Scholar

Perdew, J. P., Burke, K., and Ernzerhof, M. (1996a). Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 77, 3865–3868. doi:10.1103/physrevlett.77.3865

PubMed Abstract | CrossRef Full Text | Google Scholar

Perdew, J. P., Ernzerhof, M., and Burke, K. (1996b). [ERRATA] Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 77, 3865–3868. doi:10.1103/physrevlett.77.3865

PubMed Abstract | CrossRef Full Text | Google Scholar

Pluhařová, E., Jungwirth, P., Bradforth, S. E., and 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. doi:10.1021/jp110388v

PubMed Abstract | CrossRef Full Text | Google Scholar

Pluhařová, E., Schroeder, C., Seidel, R., Bradforth, S. E., Winter, B., Faubel, M., et al. (2013). Unexpectedly Small Effect of the DNA Environment on Vertical Ionization Energies of Aqueous Nucleobases. J. Phys. Chem. Lett. 4, 3766–3769. doi:10.1021/jz402106h

CrossRef Full Text | Google Scholar

Pluhařová, E., Slavíček, P., and Jungwirth, P. (2015). Modeling Photoionization of Aqueous DNA and its Components. Acc. Chem. Res. 48, 1209–1217. doi:10.1021/ar500366z

PubMed Abstract | CrossRef Full Text | Google Scholar

Pulay, P. (1980). Convergence Acceleration of Iterative Sequences. The Case of Scf Iteration. Chem. Phys. Lett. 73, 393–398. doi:10.1016/0009-2614(80)80396-4

CrossRef Full Text | Google Scholar

Pye, C. C., and 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. doi:10.1007/s002140050457

CrossRef Full Text | Google Scholar

Ren, X., Rinke, P., Blum, V., Wieferink, J., Tkatchenko, A., Sanfilippo, A., et al. (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. doi:10.1088/1367-2630/14/5/053020

CrossRef Full Text | Google Scholar

Rezaei, M., and Ögüt, S. (2021). Photoelectron Spectra of Early 3d-Transition Metal Dioxide Molecular Anions from GW Calculations. J. Chem. Phys. 154, 094307. doi:10.1063/5.0042106

PubMed Abstract | CrossRef Full Text | Google Scholar

Rieger, M. M., Steinbeck, L., White, I. D., Rojas, H. N., and Godby, R. W. (1999). GW Space-Time Method for the Self-Energy of Large Systems. Comput. Phys. Commun. 117, 211–228. doi:10.1016/S0010-4655(98)00174-X

CrossRef Full Text | Google Scholar

Rojas, H. N., Godby, R. W., and 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. doi:10.1103/physrevlett.74.1827

PubMed Abstract | CrossRef Full Text | Google Scholar

Sakuma, R., Miyake, T., and Aryasetiawan, F. (2009). Effective Quasiparticle Hamiltonian Based on Löwdin’s Orthogonalization. Phys. Rev. B - Condens. Matter Mater. Phys. 80, 1–8. doi:10.1103/PhysRevB.80.235128

CrossRef Full Text | Google Scholar

Salpeter, E. E., and Bethe, H. A. (1951). A Relativistic Equation for Bound-State Problems. Phys. Rev. 84, 1232–1242. doi:10.1103/PhysRev.84.1232

CrossRef Full Text | Google Scholar

Shi, B., Weissman, S., Bruneval, F., Kronik, L., and Öğüt, S. (2018). Photoelectron Spectra of Copper Oxide Cluster Anions from First Principles Methods. J. Chem. Phys. 149, 064306. doi:10.1063/1.5038744

PubMed Abstract | CrossRef Full Text | Google Scholar

Shinaoka, H., Otsuki, J., Ohzeki, M., and Yoshimi, K. (2017). Compressing Green’s Function Using Intermediate Representation between Imaginary-Time and Real-Frequency Domains. Phys. Rev. B 96, 1–8. doi:10.1103/PhysRevB.96.035147

CrossRef Full Text | Google Scholar

Shishkin, M., and Kresse, G. (2007). Self-consistent GW Calculations for Semiconductors and Insulators. Phys. Rev. B - Condens. Matter Mater. Phys. 75, 1–9. doi:10.1103/PhysRevB.75.235102

CrossRef Full Text | Google Scholar

Shishkin, M., Marsman, M., and Kresse, G. (2007). Accurate Quasiparticle Spectra from Self-Consistent GW Calculations with Vertex Corrections. Phys. Rev. Lett. 99, 14–17. doi:10.1103/PhysRevLett.99.246403

PubMed Abstract | CrossRef Full Text | Google Scholar

Singh, Y., and 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. doi:10.1103/PhysRevB.101.235157

CrossRef Full Text | Google Scholar

Staneke, P. O., Groothuis, G., Ingemann, S., and 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. doi:10.1016/0168-1176(94)04127-S

CrossRef Full Text | Google Scholar

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. doi:10.1007/BF02725962

CrossRef Full Text | Google Scholar

Tal, A., Chen, W., and Pasquarello, A. (2021). Vertex Function Compliant with the Ward Identity for Quasiparticle Self-Consistent Calculations beyond GW. Phys. Rev. B 103, 161104. doi:10.1103/PhysRevB.103.L161104

CrossRef Full Text | Google Scholar

Talman, J. D., and Shadwick, W. F. (1976). Optimized Effective Atomic central Potential. Phys. Rev. A. 14, 36–40. doi:10.1103/physreva.14.36

CrossRef Full Text | Google Scholar

Tirimbò, G., De Vries, X., Weijtens, C. H., Bobbert, P. A., Neumann, T., Coehoorn, R., et al. (2020a). Quantitative Predictions of Photoelectron Spectra in Amorphous Molecular Solids from Multiscale Quasiparticle Embedding. Phys. Rev. B 101, 035402. doi:10.1103/PhysRevB.101.035402

CrossRef Full Text | Google Scholar

Tirimbò, G., Sundaram, V., Çaylak, O., Scharpach, W., Sijen, J., Junghans, C., et al. (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. doi:10.1063/1.5144277

PubMed Abstract | CrossRef Full Text | Google Scholar

Umari, P., Giacomazzi, L., De Angelis, F., Pastore, M., and Baroni, S. (2013). Energy-level Alignment in Organic Dye-Sensitized TiO2 from GW Calculations. J. Chem. Phys. 139, 014709. doi:10.1063/1.4809994

PubMed Abstract | CrossRef Full Text | Google Scholar

Vahtras, O., Almlöf, J., and Feyereisen, M. W. (1993). Integral Approximations for LCAO-SCF Calculations. Chem. Phys. Lett. 213, 514–518. doi:10.1016/0009-2614(93)89151-7

CrossRef Full Text | Google Scholar

Van Lenthe, E., and Baerends, J. E. (2003). Optimized Slater-type Basis Sets for the Elements 1–118. J. Comput. Chem. 24, 1142–1156. doi:10.1002/jcc.10255

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Schilfgaarde, M., Kotani, T., and Faleev, S. (2006). Quasiparticle Self-Consistent GW Theory. Phys. Rev. Lett. 96, 1–4. doi:10.1103/PhysRevLett.96.226402

CrossRef Full Text | Google Scholar

Van Setten, M. J., Weigend, F., and Evers, F. (2013). The GW-Method for Quantum Chemistry Applications: Theory and Implementation. J. Chem. Theor. Comput. 9, 232–246. doi:10.1021/ct300648t

CrossRef Full Text | Google Scholar

Van Setten, M. J., Caruso, F., Sharifzadeh, S., Ren, X., Scheffler, M., Liu, F., et al. (2015). GW100: Benchmarking G0W0 for Molecular Systems. J. Chem. Theor. Comput. 11, 5665–5687. doi:10.1021/acs.jctc.5b00453

CrossRef Full Text | Google Scholar

Van Setten, M. J., Costa, R., Viñes, F., and Illas, F. (2018). Assessing GW Approaches for Predicting Core Level Binding Energies. J. Chem. Theor. Comput. 14, 877–883. doi:10.1021/acs.jctc.7b01192

CrossRef Full Text | Google Scholar

Vidberg, H. J., and Serene, J. W. (1977). Solving the Eliashberg Equations by Means of N-point Padé Approximants. J. Low Temp. Phys. 29, 179–192. doi:10.1007/bf00655090

CrossRef Full Text | Google Scholar

Vlček, V., Rabani, E., Neuhauser, D., and Baer, R. (2017). Stochastic GW Calculations for Molecules. J. Chem. Theor. Comput. 13, 4997–5003. doi:10.1021/acs.jctc.7b00770

CrossRef Full Text | Google Scholar

Vlček, V., Li, W., Baer, R., Rabani, E., and Neuhauser, D. (2018). Swift GW beyond 10,000 Electrons Using Sparse Stochastic Compression. Phys. Rev. B 98, 0–10. doi:10.1103/PhysRevB.98.075107

CrossRef Full Text | Google Scholar

Ward, J. C. (1950). An Identity in Quantum Electrodynamics. Phys. Rev. 78, 182. doi:10.1103/physrev.78.182

CrossRef Full Text | Google Scholar

Weng, G., and 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. doi:10.1063/1.1679012

CrossRef Full Text | Google Scholar

Wilhelm, J., Del Ben, M., and Hutter, J. (2016). GW in the Gaussian and Plane Waves Scheme with Application to Linear Acenes. J. Chem. Theor. Comput. 12, 3623–3635. doi:10.1021/acs.jctc.6b00380

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilhelm, J., Golze, D., Talirz, L., Hutter, J., and Pignedoli, C. A. (2018). Toward GW Calculations on Thousands of Atoms. J. Phys. Chem. Lett. 9, 306–312. doi:10.1021/acs.jpclett.7b02740

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilhelm, J., Seewald, P., and Golze, D. (2021). Low-scaling GW with Benchmark Accuracy and Application to Phosphorene Nanosheets. J. Chem. Theor. Comput. 17, 1662–1677. doi:10.1021/acs.jctc.0c01282

CrossRef Full Text | Google Scholar