In MBPT, the central quantity is the one-electron Green’s function. The Green’s function describes the time-propagation of an additional particle in the electronic system: an extra electron for propagation forward in time, or a hole for propagation backward in time. The Green’s function contains a great deal of information. For instance its diagonal is the electronic density, and, most interesting for us, its poles are the ionization energies (Fetter and Walecka, 1971).

Once an approximate Green’s function G ₀ is known, the exact Green’s function G can be obtained thanks to the Dyson equation: G ( ω ) = G 0 ( ω ) + G 0 ( ω ) Δ Σ ( ω ) G ( ω ) , (1) where the spatial indices, later defined as p and q, have been dropped for simplicity.

The operator ΔΣ stands for the self-energy difference. It performs the humongous task of connecting G ₀ to G. If the Hartree-Fock approximation (HF) is used for G ₀, then ΔΣ coincides with the missing correlation part of the self-energy Σ_c.

When a mean-field approximation is selected for G ₀, it can be expressed analytically: G 0 p q ( ω ) = δ p q 2 ω − ϵ p ± i η , (2) where the factor of 2 accounts for spin, p and q are molecular orbital (MO) indices, ϵ _p and ± iη is a vanishing imaginary number that ensures the correct analytic behavior of G ₀ . G ₀ is diagonal in the corresponding MO basis.

In practice, we make the further approximation that the self-energy difference is also diagonal in the MO basis: Δ Σ p q ( ω ) = δ p q Δ Σ p p ( ω ) . (3) This approximation is believed to be very good and is common practice in this field (Golze et al., 2019).

Recasting the Dyson Eq. 1 into G 0 − 1 ( ω ) − Δ Σ ( ω ) G ( ω ) = I , (4) where I is the identity operator, it becomes clear that the diagonal approximation of ΔΣ will induce a diagonal approximation to G, since G and therefore also G 0 − 1 are diagonal in the corresponding molecular basis.

Furthermore, the poles of G correspond to the zeroes of the term in between the brackets in Eq. 4: ω − ϵ p = Δ Σ p p ( ω ) . (5) This equation is named the quasiparticle equation and the highest zero for the p index that corresponds to occupied states is ϵ _HOMO = − IP. The HOMO energies reported in this work are obtained with this procedure, which is often referred to as the “graphical solution” of the quasiparticle equation (Golze et al., 2019).

We can calculate the spectral weight Z associated with a pole of G with Z p ( ω ) = 1 − ∂ Δ Σ p p ∂ ω − 1 . (6) Being a weight, this quantity should range from 0 to 1 and hence ∂ Δ Σ p p ∂ ω should be negative.

Note that the mean-field orbitals indexed by p might not be ordered properly. That is why in practice one needs to consider not only the mean-field HOMO, but also a few states below. This pathological behavior is known to occur for N₂ for instance (von Niessen et al., 1977).

The challenge in MBPT is then to derive approximate expressions for ΔΣ that are both accurate and computationally tractable. Henceforth, we use the Goldstone-Feynman diagram representation to describe the different working approximations. The analytic expressions can be found in the cited references.

In this Section, we follow the traditional approach in quantum chemistry for the so-called post-Hartree-Fock calculations (Szabó and Ostlund, 1996; Helgaker et al., 2000).

Let us start gradually and begin with the formulation of the HF approximation in terms of Goldstone-Feynman diagrams. In the upper panel of Figure 1, we have presented the two Goldstone-Feynman diagrams of HF: the Hartree and the Fock exchange terms. The blue arrows indicate the entry and the exit points. The black arrow is a Green’s function and the red dashed line is the bare Coulomb interaction v. As v is assumed to be instantaneous, we represent it horizontally (so that the vertical axis would be the time axis).

The Hartree diagram (upper left-hand diagram in Figure 1) translates into the following integral: v H ( r ) = ∫ d r ′ ρ ( r ′ ) | r − r ′ | , (7) Where ρ(r) = G(rt, rt⁺) is the electronic density. From the Hartree Goldstone-Feynman diagram, we can immediately see that the Hartree potential is local in space and in time, since the entry and exit points are identical. The exchange diagram (upper right-hand diagram in Figure 1) is non-local in space, but local in time, since its entry and exit points share the same y coordinate.

In regular MBPT, one considers the electron-electron interaction v as the meaningful order parameter that will allow us to derive more and more complex approximations.

The second-order perturbation theory, PT2, considers all the possible Goldstone-Feynman diagrams having two Coulomb interactions. There are only two of those diagrams and they are drawn in the middle panel of Figure 1. The first one accounts for the propagation of an electron (or a hole) interacting with an electron-hole pair. The second one is the so-called second-order exchange (SOX). These two diagrams are still rather simple and can be found in chemistry textbooks (Szabó and Ostlund, 1996).

However, the next level, namely PT3, brings in many new terms. PT3 considers all the possible Goldstone-Feynman diagrams with three Coulomb interactions, which results in the analytic terms reported in the Appendix of Ref. (Cederbaum and Domcke, 1977). The formulas extend over three printed pages and will not be reproduced here. We will instead draw a few instructive Goldstone-Feynman diagrams in the lower panel of Figure 1. PT3 contains some static diagrams (the A-diagrams in Cederbaum’s notation), such as the two first diagrams drawn in the PT3 panel. They can be interpreted as corrections to the Hartree and Fock terms due to a correction to the density and the density-matrix. Besides these, some dynamical diagrams are displayed with two electron-hole pairs, or one interacting electron-hole pair, or a ladder diagram, etc.

The PT3 approximation had been implemented and tested by Cederbaum and coworkers (Cederbaum et al., 1973; Cederbaum and Niessen, 1974; Cederbaum, 1975; Cederbaum and Domcke, 1977; Cederbaum et al., 1978), but never applied to a systematic benchmark, to the best of our knowledge. Those authors noticed that PT3 was not fully satisfactory and proposed the rescaling of some of the terms to form a better estimate of the IP. This empirical rescaling, known as outer valence Green’s function (OVGF) or as electron propagator theory (EPT), is not applied here, as our focus is the MBPT itself.

Considering the huge number of terms in PT3, it is not surprising that PT4 has only rarely been used (Ortiz, 1988).

In condensed-matter physics, it has been realized that the one-ring diagram in PT2 (See Figure 1) was producing an infinite value when evaluated for a gapless system (Mahan, 2000). A renormalized interaction was then introduced then to mitigate this problem (Baym and Kadanoff, 1961; Hedin, 1965).

The upper panel of Figure 2 represents the screened Coulomb interaction W within the random-phase approximation. W is represented with wiggly lines that are not necessarily horizontal in the diagrams, because W is not instantaneous as v is. W is an infinite series of subsequent non-interacting electron-hole pairs.

There exists only one first-order diagram in W: the so-called GW approximation to the self-energy, represented in the middle left-hand panel of Figure 2. As W contains an infinite number of diagrams, the GW approximation cannot be rationalized with the v-based MBPT recapitulated in the previous section. Notice the similarity between the exchange diagram in HF (Figure 1) and the GW diagram: The Coulomb interaction has just been replaced by a non-horizontal W wiggly line.

This single GW diagram has been proven to yield very good results for the homogeneous electron gas (Hedin, 1965; Lundqvist, 1967), and for real periodic solids (Hanke and Sham, 1975; Strinati et al., 1982; Hybertsen and Louie, 1985; Godby et al., 1986). More recently, it has been realized that the same good performance is reached for molecules (Shirley and Martin, 1993; Grossman et al., 2001; Rostgaard et al., 2010; Blase et al., 2011; Bruneval, 2012; Körzdörfer and Marom, 2012; Ren et al., 2012; Sharifzadeh et al., 2012; Bruneval and Marques, 2013; van Setten et al., 2013; Koval et al., 2014; Govoni and Galli, 2015; van Setten et al., 2015; Blase et al., 2016; Knight et al., 2016; Kuwahara et al., 2016; Maggio et al., 2017; Golze et al., 2018; Lange and Berkelbach, 2018; Wilhelm et al., 2018; Golze et al., 2019; Lewis and Berkelbach, 2019; Blase et al., 2020).

Of course, the single GW diagram is just the first of an infinite expansion in W. However, the next diagrams become very complex, very quickly. They are often named “vertex corrections” in the literature. Vertex corrections appear in two different locations in Hedin’s equations (Hedin, 1965) (or equivalently in the diagrams): in W beyond RPA and in the self-energy itself.

Adding more diagrams in W would incorporate the electron-hole interaction that is present in PT3 but not in GW. Lewis and Berkelbach have worked on this point and showed a small effect (Lewis and Berkelbach, 2019). We will test improving GW along that line by using a W interaction calculated within time-dependent Hartree-Fock (TDHF), labeled W _TDHF.

Adding more diagrams in the self-energy would incorporate the SOX diagram and more. For instance, we represent in Figure 2 the so-called second-order screened exchange (SOSEX) of Ren and coworkers (Ren et al., 2015). It is an extension to SOX that considers a screened interaction W together with an unscreened interaction v. The complete second-order diagrams with two W wiggly lines has been very recently considered in Ref. (Wang et al., 2021). The authors conclude that it does not bring large contributions and we will use the simpler SOSEX diagram here.

It should be added that there exist additional low-order diagrams when the Green’s function is not calculated self-consistently. Indeed, the two diagrams in the lower panel of Figure 2 are first-order diagrams in W that give corrections to the Hartree potential and the Fock exchange. Similar diagrams show up in PT3, however with one electron-hole pair only. These two diagrams do not appear in Hedin’s equations, because Hedin’s derivation is obtained considering the self-consistent G. One of us has recently studied these diagrams and highlighted a sizable effect on the IP (Bruneval, 2019a), on the electronic densities (Bruneval, 2019b), and on the total energies (Bruneval et al., 2021). We shall name these diagrams γ ^GW in this work, as they only affect the one-electron reduced-density-matrix.

To summarize the many approximations we have presented above, Figure 3 sketches the different diagram sets used in this study. We see that PT3 contains PT2 and that GW has an overlap with PT2, but misses the SOX diagram. Some diagrams of GW are not present in PT2, nor in PT3: the n-pair diagrams with n > 2. GW + SOSEX entirely contains PT2, but obviously misses many diagrams of PT3. GW _TDDFT captures the 1-interacting-pair diagram of PT3 and adds the further interacting pairs. GW + γ ^GW has the 1-pair inclusion in Hartree and Fock exchange. For instance, the ladder diagram is present in PT3 only.

At this stage, there is no way to judge which approximation is best. Ideally in a perturbation theory, the more diagrams, the better. However, in MBPT, the perturbation is by no means “small” and, in our opinion, only practical calculations on trusted benchmarks are able to draw conclusions. This will be the topic of the next Sections.

In this work we use the set of molecules defined in the GW100 set as our boxing ring. This set came into existence first in a comparison between only three codes (van Setten et al., 2015). In the meantime the developers of many other codes have used the set to test and benchmark their implementations, both for GW and other computational approaches aiming at the calculation of ionization energies and electron affinities (Caruso et al., 2016; Vlček et al., 2017a; Maggio et al., 2017; Wilhelm and Hutter, 2017; Govoni and Galli, 2018; Rodrigues Pela et al., 2018; Colonna et al., 2019; Gao and Chelikowsky, 2019; Brémond et al., 2020; Förster and Visscher, 2020; Gao and Chelikowsky, 2020; Bintrim and Berkelbach, 2021; Duchemin and Blase, 2021; Förster and Visscher, 2021; Wilhelm et al., 2021). At present over a hundred data sets have appeared for the GW100 set.

The GW100 set uses established geometries and keeps them fixed for each set of calculations. In the work on the GW100 set using plane-wave basis sets in the PAW formalism using the VASP code (Maggio et al., 2017), it was noticed that for two molecules, phenol and vinyl bromide, the structure used originally was not correct. From this point the two new structures have been added to the set in order to enable comparison between sets containing only one or both versions. In this work, we use the updated geometries, so that the total number of data points is 100.

For a completely correct comparison of the molecules in the GW100 set, between codes employing different basis sets, an extrapolation to the complete basis set limits is paramount (van Setten et al., 2015; Maggio et al., 2017; Govoni and Galli, 2018). However, the use of more complete basis sets that are necessary for an extrapolation is limited by the numerical scaling of the reference CCSD(T) calculations. Fortunately, for comparisons of different methods “beyond” one-shot GW in codes that are based on Gaussian orbitals, this is not strictly necessary, as long as the same basis set is used consistently. The def2-TZVPP basis set (Weigend and Ahlrichs, 2005) has historically been used for these comparisons (Krause et al., 2015; Caruso et al., 2016). We will hence use this basis set in this work as well. In this work, we refrain from interpreting small differences below 0.1 eV that could be affected by the basis set incompleteness, so that our qualitative conclusions would be equally valid for larger basis sets.

In their work providing CCSD(T) reference values for the GW100 molecules, Krause et al. also used the def2-TZVPP basis set (Krause et al., 2015). Close inspection of these results however shows that in some cases large deviations with the experimental values exist and larger than one would hope for CCSD(T). Moreover, in a number of these cases the discrepancy is larger than the one between GW and experiment. In the present comparison we need especially accurate reference energies and since also three molecular systems of the GW100 set are missing in the data by Krause et al., we start by revisiting the CCSD(T) reference set.

The CCSD(T) benchmark values for the ionization potentials, which have been used by all GW100 studies up to now, were done by the authors of Ref. (Krause et al., 2015). using an unrestricted Hartree-Fock (UHF) reference and no spatial symmetry constraints. In all cases, stability analysis was done at the UHF level to ensure that the UHF solution was indeed the lowest in energy, within their convergence tolerance.

While using the lowest energy UHF solution determined via such a stability analysis, can be a very convenient choice, it may not lead to the most accurate CCSD(T) energy. For example, in the case of S O 2 + , the UHF solution with the lowest energy at Hartree-Fock level, actually can lead to a higher energy at the frozen-core CCSD(T) level than a UHF solution with a higher energy at the Hartree-Fock level (see Table 1). While it is true that in general, a lower CCSD(T) energy does not necessarily mean a better one, the lowest energy in Table 1 is the closest one to our FCI (full configuration interaction, a numerically exact energy within the chosen basis set) estimates (Dattani, 2021), so the lowest CCSD(T) energy is actually the more accurate one in this case. Indeed, S O 2 + was one of the worst cases in the benchmark study of Ref. (Krause et al., 2015), in terms of the disagreement between CCSD(T) and experiment for the ionization energy, and it was a case where the GW calculation matched the experimental ionization energy better than the “benchmark” CCSD(T) calculations did.

In this work we have re-calculated the frozen-core CCSD(T) energies for the entire GW100 set, however we chose to use GAUSSIAN 16 (Frisch et al., 2016) (with default settings) instead of CFOUR, and the default in GAUSSIAN is an RHF (restricted Hartree-Fock) reference for all singlet species (in this paper, all neutral species), and a UHF reference for all species with a higher multiplicity (in this paper, all of the cations). This led to 46 IP values being updated with respect to Ref. (Krause et al., 2015), including the case of OCSe, for which Krause et al. accidentally used sulfur instead of selenium in their calculation. The most noticeable updates are SO₂, MgO, cytosine, and uracil with changes larger than 0.4 eV.

Our revised IPs improve very much the consistency of CCSD(T) with the related method named equation-of-motion coupled-cluster (EOM-IP-CCSD). Indeed, Lange and Berkelbach (Lange and Berkelbach, 2018) have evaluated the IPs for the complete GW100 set within this approximation and found a somewhat good agreement with Krause et al. with an MAE of 0.09 eV. However, this correct MAE is hiding a few terrible outliers, such as SO₂, MgO, cytosine, and uracil.

Now, comparing our updated CCSD(T) to Lange’s EOM-IP-CCSD yields not only an improved MAE of 0.06 eV, but also fixes all the mentioned outliers. The deviations between the updated CCSD(T) and EOM-IP-CCSD never exceed 0.30 eV.

As our updated CCSD(T) set very much improves the consistency across the methods and the comparison to experiment when experimental data are available, we have confidence that our updated values are a genuine improvement. We remind the Reader that all the numerical values are reported in the Supplemental Material.

Figure 4 shows that HF HOMO energies are too deep compared to CCSD(T), with a large spread. PT2 improves very little compared to HF: While the median is closer to zero, the mean-absolute error (MAE) remains almost as large. PT3 is a significant improvement: The median is closer to zero and the spread is reasonable. However there exists a dozen outliers with an error over 1 eV, among which the molecules containing fluorine are over-represented.

Turning to the GW approximation, the situation improves significantly. Not only is the MAE reduced to 0.3 eV, but also the spread is decreased. Furthermore, not a single outlier is identified in the whisker plot! It is striking to see how the computationally simpler GW outperforms PT3, even though PT3 contains many diagrams that GW does not have.

To understand some of the problems with PT3, let us analyze here in greater details the case of beryllium oxide. BeO is one of the worst failures of PT3, with a 2.26 eV deviation from CCSD(T). In Figure 5 we represent the correlation part of the self-energy expectation value (the right-hand side of Eq. 5) and the line ω − ϵ H O M O H F . The intersection between these two curves defines the quasiparticle energy.

In Figure 5, we observe a pathological behavior of PT3: its derivative ∂Σ _c /∂ω is sometimes positive, which is not allowed for the exact self-energy. Remember that Z introduced in Eq. 6 is a spectral weight. A positive slope yields a nonphysical spectral weight that exceeds 1. The PT3 analytic expression contains double poles, such as the C1, D1, C6, D6 terms in the Appendix of Ref. (Cederbaum and Domcke, 1977). These terms can induce this pathological behavior. PT2 and GW only contain single poles as shown for GW in Eq. 47 of Ref. (Bruneval et al., 2016) and have the correct analytic behavior by construction.

As a conclusion, based on the GW100 IP benchmark set, GW is clearly the winner by knock-out on the boxing ring: It shows the best MAE, the narrower distribution of errors, and no outlier. It has, by construction, the correct analytic behavior. Furthermore, the structure of the GW self-energy that contains only electron-hole pairs is perfectly suited for the RI approximation. An N ⁴ scaling is then achieved with the contour deformation integration technique (Mejia-Rodriguez et al., 2021) and numerical methods with better scaling have also been proposed (Foerster et al., 2011; Vlček et al., 2017b; Wilhelm et al., 2018; Duchemin and Blase, 2021). In comparison, PT2 also has N ⁴ scaling due to the infamous “atomic orbital to molecular orbital integral transform” step and PT3 has N ⁵ scaling due to the quintuple MO summations (Cederbaum and Domcke, 1977).

Now a legitimate question would be whether one could improve the GW approximation by adding some of the diagrams shown in Figure 3.

The simplest addition to GW would be to add the SOX diagram of Figure 1. This idea has already been tested by Marom et. al. (Marom et al., 2012) and was not successful according to them. In Figure 6, we confirm their conclusion: the results are better than PT2, but worse than GW alone.

Intuitively, it seems that the effect of the SOX diagram is too strong. That is why the GW + SOSEX proposal is appealing. The SOSEX diagram would temper the bare SOX. And this is precisely what it does: the spread of GW + SOSEX is narrower than that of GW + SOX. However, the results in Figure 6 show that the median and the MAE are still far from zero and that GW alone is still better.

Now let us test the possibility to incorporate the interacting electron-hole pairs, by using the TDHF screened Coulomb interaction W _TDHF. This contribution gives a significant push upwards, so that the median is close to zero. Unfortunately, many outliers appear, mostly the ionic dimers of GW100, such as LiH, LiF, BeO, MgO, FH, KH. Please note that boron nitride, BN, had to be excluded from the benchmark here. Indeed the TDHF calculation failed because of a negative excitation energy. In other words, the HF self-consistent solution reached by MOLGW is not the lowest HF energy. A stability search could solve the problem (Seeger and Pople, 1977), but this implementation is not currently available in MOLGW.

Finally, we evaluate the effect of the first-order correction to the Hartree and Fock exchange terms, as depicted in Figure 2. In agreement with previous work on a smaller benchmark (Bruneval, 2019a), we observe a significant improvement over the GW approximation. The MAE becomes very good and the distribution is well centered around zero. The only worrying point is the existence of two outliers: TiF4 and MgO. While the TiF4 HOMO was already much too negative in GW@HF (−0.62 eV compared to CCSD(T)), MgO is more intriguing. It was very good with GW (−0.08 eV compared to CCSD(T)) and deteriorates very much with GW + γ ^GW . BeO, which is chemically similar to MgO, is quite different in terms of its deviation, with a deviation of only 0.01 eV for GW + γ ^GW with respect to CCSD(T).

Of course, we did not explore all the possible combinations of diagrams beyond GW. However, we can state that with GW being already very good, it is a difficult task to improve over it. Adding diagrams may destroy the subtle balance, which makes GW so successful. Among all the additions we considered, only GW + γ ^GW can be considered as a systematic improvement.