$GW$ vertex corrected calculations for molecular systems

TL;DR

This paper implements vertex-corrected $GW$ calculations for small molecules, demonstrating improved accuracy and a new simplified approximation that avoids double counting, aligning well with more expensive wavefunction methods.

Contribution

It introduces a $GW ext{}\Gamma$ scheme with vertex corrections at multiple levels and proposes a simplified $GW^{ m{tc-tc}}$ approximation that maintains accuracy.

Findings

Good agreement with wavefunction methods.

Vertex corrections improve $GW$ accuracy.

Simplified $GW^{ m{tc-tc}}$ avoids double counting.

Abstract

Hedin's scheme is solved with the inclusion of the vertex function ($GW\Gamma$) for a set of small molecules. The computational scheme allows for the consistent inclusion of the vertex both at the polarizability level and in the self-energy. A diagrammatic analysis shows that the self-energy formed with this four-point vertex does not lead to double counting of diagrams, that can be classified as direct "bubbles" and exchange diagrams. By removing the exchange diagrams from the self-energy, a simpler approximation is obtained, called $GW^{\rm{tc-tc}}$. Very good agreement with expensive wavefunction-based methods is obtained for both approximations.

Figures

Tables

Table 1. Table 1: G 0 ​ W 0 subscript 𝐺 0 subscript 𝑊 0 G_{0}W_{0} quasi-particle energies for the highest occupied orbital (HOMO) evaluated on a Hartree-Fock reference state. In the second column the smallest plane-wave cutoff E c subscript 𝐸 𝑐 E_{c} (and the corresponding number of plane waves N pw subscript 𝑁 pw N_{\rm{pw}} ) used in the extrapolation is given. The last column reports the negative of the experimental ionisation potential (IP), vertical values are in italics. Plane wave extrapolated results include also finite box corrections as specified in section 3 . All values in eV.

molecule	$E_c$	($N_{\mathrm{pw}}$)	$G_0{W}_0$	$G_0{W}_0$ [9]	$\Delta$	CCSD(T) [9]	Expt.
			PW	GTO	PW $-$ GTO	GTO
				cc-pVQZ		cc-pVQZ
H${}_{\text{2}}$	300.0	(8620)	-16.72	—	—	-16.39 [74]	-15.43 [75]
Li${}_{\text{2}}$	112.1	(1940)	-5.29	-5.36	0.07	-5.17	-4.73 [76]
N${}_{\text{2}}$	420.9	(10060)	-16.56	-16.48	-0.08	-15.49	-15.58 [77]
P${}_{\text{2}}$	255.0	(9180)	-11.31	-10.57	-0.74	-10.76	-10.62 [78]
Cl${}_{\text{2}}$	262.5	(9780)	-12.08	-12.01	-0.07	-11.62	-11.49 [79]
HF	487.7	(12540)	-16.29	-16.39	0.10	-16.09	-16.12 [80]
LiH	300.0	(8620)	-8.26	-8.2	-0.06	-7.94	-7.90 [81]
NaCl	262.5	(7080)	-9.51	-9.36	-0.15	-9.13	-9.80 [82]
ClF	487.7	(12540)	-13.49	-13.32	-0.17	-12.82	-12.77 [79]
CH${}_{\text{4}}$	414.0	(9800)	-14.95	-14.86	-0.09	-14.4	-13.6 [83]
C${}_{\text{2}}$H${}_{\text{4}}$	414.0	(9800)	-10.91	-10.85	-0.06	-10.69	-10.68 [83]
C${}_{\text{2}}$H${}_{\text{2}}$	414.0	(9800)	-11.73	-11.65	-0.08	-11.42	-11.49 [83]
CH${}_{\text{3}}$Cl	414.0	(9800)	-11.90	-11.74	-0.16	-11.41	-11.29 [84]
CH${}_{\text{3}}$OH	434.4	(10460)	-11.71	-11.69	-0.02	-11.08	-10.96 [85]
CH${}_{\text{3}}$SH	414.0	(9800)	-9.93	-9.81	-0.12	-9.49	-9.44 [86]
SiH${}_{\text{4}}$	300.0	(8620)	-13.40	-13.31	-0.09	-12.82	-12.3 [87]
CO	434.4	(10460)	-15.03	-14.97	-0.06	-14.05	-14.01 [88]
CO${}_{\text{2}}$	434.4	(10460)	-14.35	-14.38	0.03	-13.78	-13.77 [89]
SiO	434.4	(10460)	-12.05	-11.98	-0.07	-11.55	-11.3 [90]
SO${}_{\text{2}}$	434.4	(20380)	-13.20	-13.12	-0.08	-12.41	-12.50 [84]
CS	414.0	(9800)	-12.63	-13.08	0.45	-11.45	-11.33 [91]
CS (*)			-13.19	-13.08	-0.11	-11.45	-11.33 [91]
H${}_{\text{2}}$O	434.4	(10460)	-13.10	-13.04	-0.06	-12.64	-12.62 [84]
H${}_{\text{2}}$O${}_{\text{2}}$	434.4	(10460)	-12.12	-12.13	0.01	-11.49	-11.70 [92]
H${}_{\text{2}}$S	300.0	(6040)	-10.79	-10.67	-0.12	-10.43	-10.50 [93]
HClO	434.4	(10460)	-11.96	-11.83	-0.13	-11.3	-11.12 [94]
HCN	420.9	(10060)	-13.92	-13.86	-0.06	-13.64	-13.61 [95]
NH${}_{\text{3}}$	420.9	(10060)	-11.45	-11.38	-0.07	-10.92	-10.82 [96]
N${}_{\text{2}}$H${}_{\text{4}}$	420.9	(10060)	-10.86	-10.78	-0.08	-10.24	-8.98 [97]
PH${}_{\text{3}}$	300.0	(6040)	-10.89	-10.79	-0.10	-10.49	-10.59 [98]

Table 2. Table 2: HOMO quasi-particle energies on Hartree-Fock reference for all levels of theory considered. Finite box corrections are included for all calculations as specified in section 3 . All values in eV.

molecule	finite basis set				extrapolated results
	HF	$G_0{W}_0$	$G_0{W}_0^{\mathrm{tc}-\mathrm{tc}}$	$G_0{W}_0\Gamma$	$G_0{W}_0$	$G_0{W}_0^{\mathrm{tc}}$	$G_0{W}_0\Gamma$
H${}_{\text{2}}$	-16.17	-16.34	-16.00	-16.25	-16.72	-16.40	-16.52
Li${}_{\text{2}}$	-4.90	-5.26	-5.00	-5.16	-5.29	-5.05	—
N${}_{\text{2}}$	-17.16	-16.12	-15.45	-16.06	-16.56	-15.92	-16.39
P${}_{\text{2}}$	-10.61	-11.10	-10.93	-10.92	-11.31	—	—
Cl${}_{\text{2}}$	-12.06	-11.65	-11.25	-11.52	-12.08	-11.71	-11.80
HF	-17.63	-15.83	-14.98	-15.72	-16.29	-15.43	-16.18
LiH	-8.14	-8.12	-7.39	-7.87	-8.26	-7.50	-7.94
NaCl	-9.57	-9.14	-8.65	-9.06	-9.51	-9.06	-9.32
ClF	-13.50	-13.21	-12.78	-13.14	-13.49	-13.06	-13.33
CH${}_{\text{4}}$	-14.80	-14.65	-14.24	-14.37	-14.95	-14.55	-14.57
C${}_{\text{2}}$H${}_{\text{4}}$	-10.23	-10.69	-10.48	-10.49	-10.91	-10.71	-10.66
C${}_{\text{2}}$H${}_{\text{2}}$	-11.09	-11.50	-11.25	-11.26	-11.73	-11.49	-11.43
CH${}_{\text{3}}$Cl	-11.82	-11.60	-11.18	-11.37	-11.90	-11.50	-11.56
CH${}_{\text{3}}$OH	-12.27	-11.30	-10.68	-11.06	-11.71	-11.10	-11.39
CH${}_{\text{3}}$SH	-9.68	-9.72	-9.37	-9.62	-9.93	-9.60	-9.76
SiH${}_{\text{4}}$	-13.17	-13.12	-12.70	-12.94	-13.40	-12.98	-12.88
CO	-15.15	-14.77	-14.28	-14.71	-15.03	-14.55	-14.89
CO${}_{\text{2}}$	-14.77	-13.94	-13.33	-13.83	-14.35	-13.75	-14.16
SiO	-12.65	-11.72	-10.99	-11.63	-12.05	-11.30	-11.78
SO${}_{\text{2}}$	-13.56	-12.83	-12.31	-12.64	-13.20	—	—
CS	-12.78	-12.37	-11.88	-12.29	-12.63	-12.14	-12.47
H${}_{\text{2}}$O	-13.85	-12.61	-11.87	-12.55	-13.10	-12.30	-12.92
H${}_{\text{2}}$O${}_{\text{2}}$	-13.16	-11.63	-10.81	-11.39	-12.12	-11.31	-11.81
H${}_{\text{2}}$S	-10.45	-10.51	-10.18	-10.39	-10.79	-10.47	-10.60
HClO	-12.14	-11.63	-11.12	-11.41	-11.96	-11.48	-11.66
HCN	-13.31	-13.65	-13.34	-13.40	-13.92	-13.64	-13.61
NH${}_{\text{3}}$	-11.70	-11.10	-10.51	-11.04	-11.45	-10.86	-11.28
N${}_{\text{2}}$H${}_{\text{4}}$	-11.16	-10.47	-9.90	-10.28	-10.86	-10.31	-10.58
PH${}_{\text{3}}$	-10.45	-10.64	-10.37	-10.47	-10.89	-10.62	-10.66
MAE vs CCSD(T)					0.46 $\pm$ 0.16	0.15$\pm$ 0.15	0.24$\pm$ 0.19
MSD vs CCSD(T)					-0.46	0.06	-0.24

Equations

$$G(1,2)$$

$${\Sigma }_{\mathrm{xc}}(1,2)$$

$${\Pi }_0(1,2,3,4)$$

$$W(1,2,3,4)$$

$$\Gamma (1,2,3,4)$$

$$A(1,2,3,4)=\int {\Pi }_0(1,2,5,6)W(5,6,3,4)d5d6.$$

$$I(1,2,3,4)=\frac{\delta {\Sigma }_{\mathrm{xc}}(1,2)}{\delta G(3,4)}.$$

$${\Pi }_0(1,2,3,4)$$

$$\Gamma (1,2,3,4)$$

$${\Pi }_0=L_0+L_0(\mathrm{i}I){\Pi }_0\leftrightarrow {\Pi }_0^{-1}=L_0^{-1}-\mathrm{i}I.$$

$$\Pi ={\Pi }_0+{\Pi }_{0}V\Pi \leftrightarrow {\Pi }^{-1}={\Pi }_0^{-1}-V.$$

$$\Pi =L_0+L_0(V+\mathrm{i}I)\Pi \leftrightarrow {\Pi }^{-1}=L_0^{-1}-\mathrm{i}I-V.$$

$$W=V+V\Pi V.$$

$$\begin{array}{rl}\Gamma W& =(1-(\mathrm{i}I)L_0-(\mathrm{i}I)L_0(\mathrm{i}I)L_0+...)(V+V\Pi V)\\ & =(1+(\mathrm{i}I){\Pi }_0)(V+V\Pi V)\\ & =V+V\Pi V+(\mathrm{i}I)({\Pi }_0+{\Pi }_{0}V\Pi )V\\ & =V+V\Pi V+(\mathrm{i}I)\Pi V=V+(V+\mathrm{i}I)\Pi V.\end{array}$$

$$\Pi$$

$${\Sigma }_{\mathrm{xc}}$$

$$G$$

$${\Sigma }_{\mathrm{x}}(1,3)=\mathrm{i}(GV)(1,3)=\int G(2,4)V(1,2,3,4)d2d4.$$

$$\mathrm{i}I(1,2,3,4)=-\delta (2,4)\delta (1,3)\delta (t_1,t_2)v({\mathbf{r}}_1,{\mathbf{r}}_2),$$

$$\mathrm{i}I(1,2,3,4)=-\delta (2,4)\delta (1,3)\delta (t_1,t_2)w({\mathbf{r}}_1,{\mathbf{r}}_2),$$

$$A_{iajb}^{\prime }=$$

$$B_{iajb}=$$

$$\left(\begin{array}{cc}\mathbf{A}& \mathbf{B}\\ {\mathbf{B}}^{\ast }& {\mathbf{A}}^{\ast }\end{array}\right)|\begin{array}{c}{\mathbf{X}}_{\lambda }\\ {\mathbf{Y}}_{\lambda }\end{array}⟩={\Omega }_{\lambda }|\begin{array}{c}{\mathbf{X}}_{\lambda }\\ {\mathbf{Y}}_{\lambda }\end{array}⟩.$$

$$\mathbf{\Pi }(\omega )=\mathbf{Z}{\left(\omega -\mathbf{\Omega }\right)}^{-1}\mathbf{\Delta }{\mathbf{Z}}^{†},$$

$$\mathbf{Z}=\left(\begin{array}{cc}\mathbf{X}& {\mathbf{Y}}^{\ast }\\ \mathbf{Y}& {\mathbf{X}}^{\ast }\end{array}\right),\mathbf{\Delta }=\left(\begin{array}{cc}1& 0\\ 0& -1\end{array}\right),\mathbf{\Omega }=\left(\begin{array}{cc}\mathrm{diag}\{{\Omega }_{\lambda }\}& 0\\ 0& \mathrm{diag}\{-{\Omega }_{\lambda }\}\end{array}\right).$$

$${\mathbf{Z\Delta Z}}^{†}=\left(\begin{array}{cc}{\mathbf{XX}}^{\ast }& {\mathbf{XY}}^{\ast }\\ {\mathbf{YX}}^{\ast }& {\mathbf{YY}}^{\ast }\end{array}\right)-\left(\begin{array}{cc}{\mathbf{Y}}^{\ast }\mathbf{Y}& {\mathbf{Y}}^{\ast }\mathbf{X}\\ {\mathbf{X}}^{\ast }\mathbf{Y}& {\mathbf{X}}^{\ast }\mathbf{X}\end{array}\right).$$

$$\Pi (\omega )=\sum _{\lambda }\frac{|\begin{array}{c}{\mathbf{X}}_{\lambda }\\ {\mathbf{Y}}_{\lambda }\end{array}⟩⟨{\mathbf{X}}_{\lambda }^{\ast }{\mathbf{Y}}_{\lambda }^{\ast }|}{\omega -{\Omega }_{\lambda }}-\frac{|\begin{array}{c}{\mathbf{Y}}_{\lambda }^{\ast }\\ {\mathbf{X}}_{\lambda }^{\ast }\end{array}⟩⟨{\mathbf{Y}}_{\lambda }{\mathbf{X}}_{\lambda }|}{\omega +{\Omega }_{\lambda }}.$$

$$⟨n|{\Sigma }_p(\omega )|n⟩$$

$$⟨n|{\Sigma }_h(\omega )|n⟩$$

$$⟨n|{\Sigma }_c(\omega )|n⟩$$

$${{\epsilon }_n^{\mathrm{QP}}|}_{\infty }={{\epsilon }_n^{\mathrm{HF}}|}_N+\delta ({\epsilon }_n^{\mathrm{QP}})+C_N({\epsilon }_n^{\mathrm{QP}}).$$

Full text

1 Introduction

In the past decades many theoretical methods have been developed in the attempt to predict and rationalise molecular electronic structures. Coupled cluster approaches are amongst the most widespread reference methods and are based on the exponential Ansatz for the molecule’s wavefunction [1, 2]. The unfavourable scaling of these methods, however, makes them unsuitable for large scale calculations or predictive assessments for a large number of molecular systems.

In contrast to these well-established methods, many-body perturbation theory (MBPT) has emerged as a computational alternative for molecules [3, 4, 5, 6, 7, 8, 9, 10, 11, 12]. The computational scheme first proposed by Hedin [13], and summarised in Figure 1 (a-b), involves several quantities: starting with an initial guess for the Green’s function ($G$) one determines the system’s irreducible polarizability ($\chi_{0}$) and the screened interaction $W$. These quantities alone fully specify the self-energy operator $\Sigma$ in the first iteration of the scheme, which then produces an updated Green’s function, as shown in panel (b). For later iterations the inclusion of vertex effects is in principle required through the vertex function $\Lambda$. The main difficulty in including the vertex is that it involves additional diagrams at each step, since the vertex is formally generated by the functional derivative of the self-energy with respect to the Green’s function. Exploiting exact functional relations [14] can simplify the calculations only to a certain extent, and no purely numerical scheme can be implemented for the original method of Hedin. It is understandable then that the $GW$ approximation (GWA) completely neglects the vertex, thus iteratively applying the scheme in Figure 1-b. This approximation is justified when exchange and correlation effects are completely absent in the reference state and only for the first iteration of the scheme, i.e. when $\Sigma_{\mathrm{xc}}=0$, as for the Hartree case [15]. However, it has been widely used also on Hartree-Fock (HF) or density functional theory (DFT) starting points [16, 9, 17, 18, 19, 20].

Actual GWA calculations vary significantly in details. The most common approximation is to create the starting orbitals using standard density functional theory and perform a single shot $G_{0}W_{0}$ calculation. This works well for extended systems, but for small molecules this approximation tends to yield underbound highest molecular orbitals (HOMO). One way to cure this problem is to perform the calculations self-consistently. An update of the quasi-particle (QP) energies alone produces what is now often called the ev-$G_{n}W_{0}$ or the ev-$G_{n}W_{n}$ schemes, if the updated eigenvalues are included only in the Green’s function or also in the screened interaction, respectively. Fully self-consistent $GW$ calculations [8, 21, 22] require the calculation of the interacting Green’s function $G(\omega)=(\omega-T-V_{\rm{ext}}-\Sigma)^{-1}$, with $T$ the kinetic energy operator and $V_{\rm{ext}}$ the external potential generated by the atomic nuclei and $\Sigma$ being the self-energy. The resulting quasi-particle equations involve an energy dependent non-Hermitian self-energy. To make this problem more manageable, quasi-particle self-consistent $GW$ approaches [23, 24] have been devised that disregard the non-Hermitian component of the self-energy and approximate its energy dependence [25]. Recent studies on selected organic molecules have shown that self-consistency significantly improves the quasi-particle energies[21, 26, 27]. However, the HOMO is now on average too negative; in other words self-consistent approaches tend to over-bind the electrons and yield a too large ionization potential (IP). As we will discuss below the inclusion of the vertex cures this shortcoming already in the first iteration of the scheme.

Apart from the overestimation of the binding energies, self-consistency is also computationally fairly expensive, and other simple means to improve the QP energies have been considered as well. One is the admixture of exact exchange in the reference mean-field ground state calculation. From a fundamental point of view, this should certainly improve the occupied orbitals, since for the standard local and semi-local density functionals that are commonly used as starting point the exchange-correction potential does not properly decay like one over the distance from the nuclei [28]. Also the HF one electron energies are closer to the final $GW$ QP energies, so that an iterative solution of the $GW$ equations should not be necessary. An elegant mean to achieve the correct $1/r$ decay of the nuclear-electron potential is using long-range hybrid functionals [9, 26], but standard hybrid functionals[27] as well as Hartree-Fock starting points have been used as well.

In the present work, we have also decided to use the HF reference point to initiate the $GW$ calculations. The reasons for this choice are summarized below. (i) From a practical point of view, this choice obviously allows for a rather straightforward comparison with wavefunction based methods that usually rely on the HF reference point. Accordingly, we will use CCSD(T) reference data taken from the literature [9]. (ii) There is evidence that HF theory is a reasonably accurate starting point for small molecules: it is well known that, because of the Brillouin condition in HF theory, single excited configurations only contribute in second order in the perturbative expansion of the exact wavefunction, whereas for Brückner orbitals their contribution is exactly zero to all orders. From a diagrammatic standpoint this translates into an identical cancellation of single excited determinants in lowest order for the ground state energy evaluated using either the HF or the Brückner orbitals [29]. Since in localised systems the contribution of correlated higher order excitations is small [30], one can reasonably expect that the HF determinant has a large overlap with the exact many-body wavefunction, similarly to the Brückner orbitals. (iii) The main point of the present work is to investigate how important vertex corrections are. For the HF starting point, the vertex corrections can be evaluated exactly, as they are give by the functional derivative of the exchange potential with respect to the Green’s function. Hence we can unambiguously study the effect of vertex corrections.

This brings us back to a point already alluded to above. $G_{0}W_{0}$@HF, as well as self-consistent $GW$ place the HOMO at too negative binding energies, or in other words, the IP is too large [9]. The reason for this is fairly simple. In the GWA— which implies that the random phase approximation is used to calculate the polarizability —particle-hole ladder diagrams are not included. The resulting $W$ underscreens the interaction between electrons because the energy to create an electron-hole pair is too large. This largely explains why the IP is overestimated. The second issue affecting the GWA stems from the so-called self-screening error [31, 32, 33]. In the $GW$ self-energy, one state can be occupied by two particles, since direct and exchange interactions are not treated on the same footing [34]. This situation is reminiscent of the self-interaction problem that appears at the mean-field level in the Hartree theory. The inclusion of the exact exchange in the Hartree-Fock theory completely cures this shortcoming on the mean-field level but neglects correlation. Including the functional derivative of the exchange with respect to the Green’s function in Hedin’s equations also eliminates the problem on the level of correlations.

To set the stage, we now summarize our computational approach. The initial step is a HF calculation for a set of molecules. We then calculate the standard $G_{0}W_{0}$@HF QP energies for this set, and find as expected that this approximation overestimates the IP. We then present post-$GW$ calculations for the same set of molecules. In these calculations, the vertex corrections are included and given by the functional derivative of the exchange potential $\Sigma_{\mathrm{x}}$ with respect to the HF Green’s function $G(\omega)=(\omega-T-V_{\rm{ext}}-V_{H}-\Sigma_{\mathrm{x}})^{-1}$ (here $V_{H}$ is the Hartree potential). The vertex is included either only in the polarizability ($G_{0}W_{0}^{\rm{tc-tc}}$) or in the polarizability and in the self-energy ($G_{0}W_{0}\Gamma$). In present work, only the diagonal components of the self-energy are calculated, and changes of the orbitals are concomitantly not considered, nor do we perform the calculations self-consistently, since this would in principle require a consistent update of the vertex as already emphasized above.

The inclusion of the vertex in the polarizability is equivalent to calculating the polarizabilty using time-dependent Hartree-Fock [35], and the computational cost is also similar with both methods scaling like $O(N^{6})$. There are also formal similarities to the Bethe-Salpeter equation (BSE) performed on top of a $GW$ reference state. The latter method is commonly adopted for solid state systems [36, 37], but has also been used to study the optical properties of molecular systems [38, 39]. To the best of our knowledge, vertex corrections have been routinely considered only for extended systems [40, 41, 42, 15, 43, 44, 22]. For finite systems, calculations have been performed for atoms [45], for very simple molecules (within the Tamm-Dancoff approximation) [46], or with the SOSEX [47, 26], which approximates the vertex only in second order. A local vertex $\Lambda_{\rm{LDA}}$ [48] has also been used for a set of aromatic molecules, however, this simple two-point vertex behaves quite differently than the four-point many-body vertex used here. Among other things, the local approximation breaks fundamental self-energy symmetries and, albeit these effects are small for extended solids [41], they might be more relevant for very inhomogeneous systems such as molecules.

The work is structured in the following manner. In section 2 we derive simple equations for the vertex corrected $GW$, reducing the set of equations to three equations that are reminiscent of the equation of motion. Contrary to Hedin’s equation we use a four-point notation. In subsection 2.3 we report a diagrammatic analysis of the self-energy thus generated. The computational details are given in section 3, and our results are finally presented and discussed in section 4, where we apply the $GW\Gamma$ method perturbatively to a set of molecules.

2 Theory

2.1 Hedin’s equations

In principle, there are at least two complementary routes to improve on the MBPT approaches: either a self-consistent evaluation of the standard GWA or the inclusion of sub-leading diagrammatic contributions. Practical implementations have been directed primarely towards self-consistency (see section 1 for a discussion), whereas vertex corrections remain largely unaddressed, because of the high implementational as well as computational complexity.

To present our computational approach we will rely on a four point notation, for instance used by Starke and Kresse [49], although in the present work we have decided to rearrange the indices to more easily connect with the available standard literature. Also the specific succession of indices adopted here is easier to memorise. Using a four index notation Hedin’s equations can be written as

[TABLE]

In these equations, integrals over repeated indices are assumed. We note that if two four-point quantities are multiplied, for instance, $\Pi_{0}W$ in Eq. (2.1), then the integral is always over the two intermediate indices, in the equation below 5 and 6, and the order is easy to memorise:

[TABLE]

The advantage of using a four point notation is the following. In Hedin’s equations a numeral, for instance $1$, corresponds necessarily to a space-time and spin point $1=({\bf r}_{1},t_{1},\sigma_{1})$; Hedin’s original equations do not apply in orbital space or momentum space. In the four point notation above, however, one can perform a unitary transformation of the position coordinate to reciprocal space or to any set of orthogonal orbitals without changing the equations. The slight disadvantage of the present ordering of the indices is that the first index and the fourth index in each four point object transform like covariant coordinates (transformed by say $U$), whereas the second and third indices transform like contra-variant coordinates (transformed by $U^{\dagger}$).

Before continuing we briefly reiterate the meaning of the individual objects. The four-point Coulomb interaction in space-time coordinates is $V(1,2,3,4)=\delta(1,2)\delta(3,4)\delta(t_{1},t_{4})v(\mathbf{r}_{1},\mathbf{r}_{4})$ generalising the usual two-point counterpart $v(\mathbf{r}_{1},\mathbf{r}_{4})=\frac{1}{|\mathbf{r}_{1}-\mathbf{r}_{4}|}$. $G_{0}(1,2)$ is the Hartree Green’s function, generated by the corresponding Hartree self-energy: $\Sigma_{0}(1,2)=$ $V_{H}(1,2)=$ $\mathrm{i}\delta(1,2)\delta(t_{1},t_{3})[v(\mathbf{r}_{1},\mathbf{r}_{3})G(3,3^{+})]$ in space-time coordinates, $\Sigma_{\mathrm{xc}}(1,2)$ is the self-energy including exchange as well as correlation terms. $\Pi_{0}$ is the irreducible polarization (propagator); it can not be divided into two individual polarization propagators by cutting a single Coulomb line $V$. $W$ is the screened interaction summing the Coulomb interactions to infinite order. Finally $\Gamma$ is the four-point vertex, which is completely specified by the kernel $I$:

[TABLE]

There is an alternative way to write these equations that avoids the vertex altogether. We first introduce an auxiliary independent two-particle propagator $L_{0}(1,2,3,4)=-G(1,3)G(4,2)$. This allows to rewrite Eqs. (3) and (2.1) as

[TABLE]

We can now combine both equations to a single equation

[TABLE]

Here we have suppressed the indices and the matrix multiplications need to be done according to the rule stated above (Eq. 6). For the inversion, it is understood that the first two and last two indices of each four point quantity are combined to a super index and the matrices are inverted using the two super indices as row and column indices. Likewise, we can introduce the full polarization propagator

[TABLE]

It is a simple matter to see that these two equations can be further combined into a single equation, reading

[TABLE]

The proof is most easily done inspecting the right hand inverted Dyson like equations in the previous three equations. Eq. (10) is precisely the polarization propagator as calculated by means of the “Bethe-Salpeter equation” (BSE) in most solid state codes. In fact, in the preceding lines, we have backtracked the calculations performed by Starke and Kresse, where Hedin’s equations were derived from the equation of motion and the BSE.

We can now proceed with $W$ and the self-energy. By expanding the right hand side of Eq. (2.1) to infinite order in $V$ one can identify $\Pi$ and rewrite Eq. (2.1) as

[TABLE]

The final issue is to rewrite Eq. (2), specifically the four point term $\Gamma W=\Gamma(1,2,3,4)W(3,4,5,6)$ such that it involves only the already calculated polarization propagators. Using Eq. (11) for $W$ and inserting the vertex from Eq. (2.1), we obtain:

[TABLE]

From the first to the second line we have used Eq. (8), and from the third to the fourth line we have used Eq. (9).

In summary, we can rewrite Hedin’s equation in the much more compact form

[TABLE]

where the contraction of the two-point Green’s function and a four-point quantity is here defined as

[TABLE]

In these equations the somewhat arbitrary distinction between the Coulomb kernel $V$ and the remaining interaction kernel $\mathrm{i}I$ has been dropped. These equations are essentially equivalent to the compact equations given by Starke and Kresse but have been derived here from Hedin’s equations instead of the more fundamental equation of motion. We note that the same equation for the self-energy can be also found in Ref. [50] (albeit without derivation) and is seemingly common knowledge in quantum field theoretical publications [51, 52].

A few comments are in place. As already emphasized, $\Pi$ in Eq. (12) is the polarization propagator (or polarizability) that many solid state BSE codes calculate. It describes how test charges are screened by the electronic system, hence it is often referred to as test-charge test-charge (tc-tc) polarizability. The self-energy Eq. (13) describes the effects of the many-body system on an added or removed electron, and Eq. (14) is the related Green’s function. The first term on the right hand side of Eq. (13) is just the exchange potential ($GV$), whereas the second term describes, in principle, all correlation effects. If $iI$ is neglected in the self-energy (13) but included in (12) the approximation is commonly referred to as $GW^{\rm{tc-tc}}$.

Although this set of equations is in principle exact, approximations need to be made for the interaction kernel $I$ specified in Equation 7. In our case, we perform calculations perturbatively on the Hartree-Fock reference state, i.e. we start the Hedin scheme with the sum of the Hartree ($\Sigma_{0}$) and exchange self-energy ($\Sigma_{\mathrm{x}}$): $\Sigma=\Sigma_{0}+\Sigma_{\mathrm{x}}$. Then the interaction kernel is simply the functional derivative of the the exchange self-energy $\Sigma_{\mathrm{x}}$ and thus explicitly given by

[TABLE]

where $v(\mathbf{r}_{1},\mathbf{r}_{2})$ is the bare Coulomb interaction between particles. The polarization propagator is then just the polarizability of time-dependent Hartree-Fock [35]. The approximation employed in this study is also related to the random phase approximation with exchange (RPAx) first proposed by Szabo and Ostlund [53]. There has been a recent revival in interest for the RPAx [54, 55, 56, 57, 58, 59], however, previous studies have focused on the evaluation of correlation energies, rather than the self-energy considered here.

Finally, we comment on the use of the $GW$ reference state instead of the Hartree-Fock starting point employed here. For the $GW$ reference state, the interaction kernel is given by the functional derivative of the $GW$ self-energy with respect to the Green’s function. This makes the interaction kernel more complicated. Specifically, (i) $W$ is frequency dependent, and (ii) since $W$ itself depends on the Green’s function, the derivatives of $W$ with respect to $G$ should be accounted for. The resulting diagrams are sometimes refereed to as “butterfly” diagrams. A common approximation is to neglect the dependence of $W$ on the Green’s functions, and to approximate the frequency dependence of $W$ by an instantaneous interactions. The most common approximation for the interaction kernel is

[TABLE]

with $w(\mathbf{r}_{1},\mathbf{r}_{2})$ being the random phase approximation (RPA) screened Coulomb interaction at $\omega=0$. These approximations are also commonly employed in solid state BSE codes [36, 37, 60, 61, 62], and our present $GW\Gamma$ code can be adopted to this case. We will report on such calculations in forthcoming publications.

2.2 Bethe-Salpeter equation

As already emphasized above, we treat all the interactions as instantaneous in the present case. For the Hartree-Fock reference point this is exact, but it constitutes an approximation, if one were to start from the self-energy of the $GW$ approximation. In full generality, the BSE can be Fourier transformed into the frequency domain and its kernel becomes a function of three independent frequency variables $I=I(\omega,\omega^{\prime},\tilde{\omega})$, with the screened interaction carrying a frequency dependence of $\omega^{\prime}-\tilde{\omega}$ [63]. Using an instantaneous interaction obviously yields a kernel $I$ that does not depend on the frequency, since the Fourier transformation of a $\delta$-function in time is a constant in frequency space. More detailed discussions can be found in literature [64, 65, 66]. Then, in an orbital representation the full kernel $V+\mathrm{i}I$ is given by the matrix elements of the matrices $\mathbf{A^{\prime}}$ and $\mathbf{B}$:

[TABLE]

The BSE reduces to a generalised eigenvalue problem [60], where the eigenvalues $\Omega_{\lambda}$ correspond to the optical transition energies, and the matrix $\mathbf{A}$ is formed from $\mathbf{A^{\prime}}$ adding the independent particle energy differences $\Delta E_{ia}=\epsilon_{a}-\epsilon_{i}$ to the diagonal elements:

[TABLE]

At the $\Gamma-$point the matrix elements of $\mathbf{A}$ and $\mathbf{B}$ can be chosen to be real-valued and since the matrix $\mathbf{A}$ is Hermitian and $\mathbf{B}$ is symmetric [67], the generalised eigenvalue problem (EVP) above can be recast as an Hermitian EVP, for which standard solvers are available. A concise expression for the resulting polarisation propagator’s spectral representation reads [68, 69]:

[TABLE]

with the following definitions for the matrices:

[TABLE]

The matrix $\mathbf{Z}$ contains all the individual eigenvectors $\mathbf{X}_{\lambda}$ and $\mathbf{Y}_{\lambda}$, similarly the matrix $\mathbf{\Omega}$ includes resonant and antiresonant transition energies; the transition energies for the resonant and antiresonant branch are identical (apart from a sign) as a consequence of the downfolding of the original BSE into an Hermitian EVP. For definitiveness we choose $\Omega_{\lambda}>0$. Since the matrix $\omega-\mathbf{\Omega}$ is diagonal in the orbital representation, it can be commuted across; the resulting product $\mathbf{Z}\mathbf{\Delta}\mathbf{Z}^{\dagger}$ can be decomposed into the difference of two matrices:

[TABLE]

The resonant and antiresonant contributions can be more easily singled out by re-writing the matrices above as an external product of the eigenstates of Equation 19. The polarisation propagator then reads:

[TABLE]

This expression can now be inserted into Equation 13 to give the greater and lesser components of the correlation self-energy, corresponding to the propagation of a particle ($p$) and a hole ($h$), respectively:

[TABLE]

$f(\epsilon_{n})$ in the previous equations is the occupancy for the energy level $\epsilon_{n}$ and the summation goes over the set of particle-hole excitations $\lambda$. The composite index $\lambda$ is constructed by considering all possible combinations of the single particle indices $(i,a)$, with $i$ and $a$ belonging to the occupied and unoccupied orbital manifold; in the following $n$ represents a generic orbital. The total correlation part of the self-energy is then obtained by taking the Hilbert transform of the two components:

[TABLE]

Expressions in Eqs. (23) and (24) reduce to the usual RPA limit [70], if the exchange term in $\mathbf{A^{\prime}}$ and $\mathbf{B}$ is omitted.

2.3 Diagrammatic analysis

We start our diagrammatic analysis with the lowest order expressions: these can be obtained by replacing the fully interacting polarisation propagator $\Pi$ in Equation 12 with the non-interacting counterpart $L_{0}$. The corresponding diagrams are shown in Figure 2. In the top row of Figure 2 we are reporting the correlation contributions in Equation 13 that stem from the term labelled with $V$, in the bottom row the terms generated by the kernel $\mathrm{i}I$ are shown. The two columns portray the different orientations of the intermediate state $n$ with respect to the incoming state $n^{\prime}$ (which in the diagram is represented by the external legs which have not been amputated for clarity): in the left column $n$ and $n^{\prime}$ belong to the same manifold (they are both either occupied or unoccupied), in the right column $n$ and $n^{\prime}$ belong to different manifolds. It is obvious, by inspecting the diagrams, that the two contributions (in the top and bottom row, called respectively bubble and exchange in the following) are non equivalent in lowest order and that by cutting a particle-hole pair the diagram does not break up into disjoint pieces. This lower order contribution results in the MP2 self-energy, evaluated for instance in Ref. [44] for extended systems.

If we now switch on the interaction between particles and holes in the polarisation propagator we obtain additional classes of diagrams, as shown in Figure 3 in third order in the interaction, i.e. setting $\Pi=L_{0}(V+\mathrm{i}I)L_{0}$. Two of the third order diagrams belong to the same class as their second order counterparts, these are shown in the left column of Figure 3. Additionally, the ”mixed contributions” shown in the right column represent the insertion of an exchange interaction in the bubble term (top) and the insertion of a bubble in the exchange diagram (bottom). The polarisation propagator, however, is obtained by a resummation to infinite order of the corresponding diagrams via the solution of the BSE. These are represented in Figure 4, where we are showing both orientations for the intermediate particle for completeness. Starting from the top left diagram in Figure 3 it is easy to identify an infinite resummation of this class with the usual RPA approximation, where the so-called bubble diagrams are summed up; this translates into a screened interaction (denoted by a double wavy line) showing up in the corresponding diagram in Figure 4. Likewise, an infinite number of bubble insertions in the exchange diagram of Figure 3 (bottom right) returns an additional screened line. The infinite resummation of the exchange contributions (top right and bottom left in Figure 3), on the other hand, takes into account the interaction between particles and holes in the polarisation propagator and this is represented in Figure 4 by a shaded triangle. We should now point out that in the final self-energy diagrams, the resummation to infinite order of the polarisation propagator generates a frequency dependent, dynamical screened interaction vertex (pictorially defined at the bottom of Figure 4). Considering the resulting diagrams it is easy to see that they result in proper, irreducible self-energy insertions: the presence of bubble and exchange contributions in the kernel assures that two topologically distinct classes are obtained (shown in the top and middle row of Figure 4) and that for either class there is no external polarisation insertion.

As it was mentioned in the introduction, the iteration to self-consistency modifies the diagrammatic structure presented up to this point (see Figure 1-a). This is one reason why we do not consider self-consistency in the present work. As explained at the end of Sec. 2.1, by adopting some further approximations, it is however possible to maintain the simple algebraic structure and extend the present approach to $GW\Gamma$ calculations on top of a $GW$ reference state.

Finally, we mention that the exclusion of the exchange diagrams in the self-energy can be obtained by formally replacing the electron in the intermediate state $n$ with a classical test-charge. In what follows we will refer to this level of theory as to the $GW^{\rm{tc-tc}}$ approximation. This approximation neglects the $\mathrm{i}I$ term in Equation 13, but includes it in Equation 12.

3 Computational procedure

Calculations were performed on a subset of molecules studied in Ref. [9]. For consistency we adopted the same molecular geometries as obtained in Ref.[71] at the MP2/6-31G(d) level of theory. The initial step was always a self-consistent Hartree-Fock calculation. Since the solution of the Bethe-Salpeter equation is very demanding if many virtual orbitals are included, we used natural orbitals to represent the unoccupied states. In our procedure we roughly followed the work of Grüneis et al. [44], but instead of using the MP2 density matrix, we calculated the RPA density matrix using the cubically scaling $GW$ code [72]. We then diagonalized the virtual-virtual sub-block and sorted the natural orbitals by their occupancy. In this step, the number of orbitals is reduced by a factor four compared to the plane wave basis set size specified by the cutoff $E_{c}$ and reported in Table 1. We note that the occupied ($N_{\rm{occ}}$) as well as few unoccupied states are not modified but kept fixed at the level of the HF canonical orbitals. The number of unoccupied orbitals that are kept fixed at their Hartree-Fock level is set to 0.8$\times N_{\rm{occ}}$. By comparing calculations using all canonical orbitals and natural orbitals for selected systems, we found that the reduction of the unoccupied states introduces errors that are smaller than about 10 meV. After determining the natural orbitals, the Hartree-Fock Hamiltonian is again diagonalized using the natural orbitals as basis 111in VASP this is performed by selecting ALGO=SUBROT. Using this basis, the Bethe-Salpeter equation is solved for the polarization propagator, and the diagonal elements of the self-energy are calculated in the orbital basis. The self-energy is evaluated on the real frequency axis sampled with a uniform grid of 0.1 eV spacing. The complex broadening is set to twice the frequency spacing and the resulting self-energy $\langle n|\Sigma_{c}(\omega)|n\rangle$ (shifted by the kinetic, ionic, Hartree and Fock exchange energy contributions) is piecewise linearly interpolated in each interval to determine the quasi-particle energies as the intersection with the bisector line as shown in Figure 5. The quasi-particle shifts ($\delta(\varepsilon^{\rm{QP}})$ in the following) are then obtained by subtracting the corresponding mean-field single particle energies. To minimise the impact of image charges due to the finite simulation box, the Hartree-Fock calculations were repeated increasing the simulation cell up to 25 Å in each linear dimension and then correcting for the residual local potential present at the edge of the cell. The quasi-particle calculations were performed in a smaller simulation box (8 Å in size). The volume dependence of the quasi-particle energies has been carefully analysed in our previous publication [73] and it was found to be neglible for the systems considered here. The quality of this approximation can be assessed for the case of $G_{0}W_{0}$ for the HF reference state by comparison with the localised basis set calculations in Ref. [9] and reported in Table 1.

The final quasi-particle energies were obtained by adding to the cell size converged (HF) single particle energies the quasi-particle shifts and the basis set corrections, according to the expression:

[TABLE]

The basis set correction $C_{N}$ is defined as: $C_{N}(\varepsilon)=\left.\varepsilon\right|_{\infty}-\left.\varepsilon\right|_{N}$, where the subscripts indicate the number of natural orbitals ($N$) or the extrapolated value ($\infty$). In the next section it will be shown how the interplay of quasi-particle energy shift and the basis set correction affects the final quasi-particle levels.

One of the advantages of plane wave basis sets is the full control of the basis set completeness by specifying their kinetic energy cutoff $E_{c}$. 222in VASP this specified by the flag ENCUT To estimate the basis set convergence, we increased $E_{c}$ by a factor 1.3 and 1.5 beyond the value specified in Table 1 and extrapolated assuming that the QP shifts converge like one over the basis set size [16].

4 Results and discussion

In this section, we apply the computational scheme presented above to a subset of the molecules considered in Refs. [9, 71]. As a first check on the accuracy of the computational procedure we compare in Table 1 the HOMO quasi-particle energies at the $G_{0}W_{0}$ level of theory with those obtained by Bruneval and Marques with Gaussian type orbitals (GTOs) using cc-pVQZ basis sets [9].

We immediately point out that we expect differences of the order of 100 meV between PW and GTO $GW$ results. First, pseudopotentials can introduce errors and second the GTO results were not extrapolated to the complete basis set limit. Concerning the pseudpotential error, we know from previous work that flourine containing compounds are affected most strongly with errors of about $\approx$100 meV, owing to the strong degree of localisation for 2$p$ electrons. With the non-normconserving PAW potential used here, we somewhat underestimate the degree of localization for flourine, as explained elsewhere [73]. Concerning the basis set extrapolation of the GTO results, its convergence was studied in Ref. [9] up to the ”correlation-consistent” pV5Z basis set for the specific case of CO; the HOMO energy difference between this basis set and the quadruple zeta employed for the calculations of all the remaining molecules is about $\sim$70 meV. Since the error decays like $1/C_{n}^{3}$, where $C_{n}$ is the cardinal number (i.e. four for quadruple basis sets), the quadruple zeta basis set results reported in Ref. [9] are expected to have errors of the order of $2\times 70~{}$meV$\approx 140~{}$meV. The difference between our plane-wave (PW) results and the GTO ones is typically comparable with this order of magnitude, with the exception of two outliers: CS and P${}_{\text{2}}$. For CS we observe that in going from the HF to the $G_{0}W_{0}$ level there is a rearrangement of the orbital energies, with the HOMO-1 and the HOMO exchanging their positions. We can therefore speculate that the value reported in the literature refers to a lower lying quasi-particle level that corresponds to the HOMO orbital at the mean-field level. For this reason we have included also this quasi-particle energy in Table 1 for comparison. On the other hand, for P${}_{\text{2}}$, we notice how the finite box effects are particularly strong in this case, with a downward shift exceeding 1 eV for the HF HOMO level as the cell edge is increased from 8 to 25 Å. It is, however difficult to draw a solid conclusion on the origin of the observed mismatch, also considering that in Ref. [9] P${}_{\text{2}}$ fails to follow any of the trends observed for the remaining molecules.

We then consider the impact of the vertex on the self-energy and on the quasi-particle energy for the HOMO level. This level is defined as the highest occupied orbital in the (post-)$GW$ approximation. This coincides with the HF HOMO level for all molecular systems except CS and N${}_{\text{2}}$, for which it corresponds to the HOMO-1 level in the HF approximation. Furthermore, the inclusion of the ladder diagrams only in the polarisation propagator produces the so-called $G_{0}W_{0}$ with test charge-test charge interactions (labelled $G_{0}W_{0}^{\rm{tc-tc}}$ in the following). In terms of the diagrams involved, the latter approximation can be obtained from the more general $G_{0}W_{0}\Gamma$ by omitting the exchange diagrams in the self-energy (e.g. those in the bottom row in Figs. 2 - 4 above). For this level of theory we observe the same level alignment between the quasi-particle HOMO level and HF reference as in the $G_{0}W_{0}\Gamma$ calculations, with the sole exception of SiO. For this molecule there is a level crossing as in the case of CS and N${}_{\text{2}}$.

The inclusion of the ladder diagrams in the polarizability is reasonable from a physical standpoint: it introduces an interaction between the ”virtual” particle and hole that are generated as a response of the many-body system to the introduction of a test charge (excitonic effects). In turn, ladder diagrams increase the screening in the many-body system, since an interacting electron-hole pair requires less energy to be generated than a non-interacting pair. This effect will be present in both the $G_{0}W_{0}\Gamma$ and the $G_{0}W_{0}^{\rm{tc-tc}}$, and it shows on the resulting self-energy as a shift of its resonances to higher (lower) energies below (above) the Fermi level (conventionally located mid-gap between the HOMO and LUMO levels). Figure 5 contrasts the HOMO self-energy for all approximations considered in two representative systems (PH${}_{\text{3}}$ and CO): in both cases, there is a reduction of the frequency interval devoid of resonances when the polarizability includes ladder diagrams. Consistently the intersection with the bisector line $y=\omega$ occurs in the $G_{0}W_{0}^{\rm{tc-tc}}$ approximation at higher (less negative) energies, thus reducing the underscreening error of $G_{0}W_{0}$, as shown in Table 2. Not surprisingly, the most sizable changes occur for molecules whose highest occupied orbital (HOMO) has a non-bonding or anti-bonding character, i.e. that has a larger overlap with the LUMO and higher unoccupied states. For instance: H${}_{\text{2}}$O, H${}_{\text{2}}$O${}_{\text{2}}$, HF, CO${}_{\text{2}}$, NH${}_{\text{3}}$, have quasi-particle $G_{0}W_{0}^{\rm{tc-tc}}$ shifts more positive than the $G_{0}W_{0}$ values by more than 600 meV as reported in Figure 6.

On the other hand, in the $G_{0}W_{0}\Gamma$ approximation (shown in Figure 5 by dashed-dotted lines) there is another effect at play: the inclusion of exchange diagrams in the self-energy redistributes the spectral weight of its resonances across the frequency range. This effect, for the vast majority of systems considered here (the only exceptions being the unsaturated hydrocarbons C${}_{\text{2}}$H${}_{\text{4}}$ and C${}_{\text{2}}$H${}_{\text{2}}$), counteracts the upshift of $G_{0}W_{0}^{\rm{tc-tc}}$ in comparison with $G_{0}W_{0}$ quasi-particle energies. In particular the effect is more pronounced for the set of molecules that showed larger quasi-particle shifts from $G_{0}W_{0}$ to $G_{0}W_{0}^{\rm{tc-tc}}$ (see Figure 6). This comparison can be systematically carried out as shown in Figure 6 for the whole range of molecules considered: the impact of vertex corrections in the self-energy, i.e. in going from the $G_{0}W_{0}^{\rm{tc-tc}}$ to the $G_{0}W_{0}\Gamma$ approximation, is system dependent. In particular, systems with a more ionic character (for instance fluorides and chlorides) or with unsaturated chemical bonds (e.g. N${}_{\text{2}}$, CO, CO${}_{\text{2}}$) will have the electron density rather localised around the more electronegative element or in the bonding region respectively. As electrons become more localised correlations among them become pivotal and necessitate the inclusion of the vertex at the self-energy level for an accurate description.

Comparison with previous literature is not straightforward. Our observations are in contrast with the well-established behaviour of the local vertex in the $GW\Lambda_{\rm{LDA}}$ approximation, which results in a rigid shift of the quasi-particle energies [15, 48, 41], regardless of the system’s details. $GW\Gamma$ calculations have also been reported by Shirley and Martin [45] for neutral atoms and ions. Specifically, it has been reported that in ionized atoms the importance of exchange effects in the polarisability exceeds that of correlation effects in the self-energy, whereas for (neutral) atoms correlation effects in the self-energy are more relevant. With the very diverse systems considered here, we can not substantiate this trend.

Now we evaluate the basis set corrections for the various levels of theory. We notice these values are consistenly smaller for $G_{0}W_{0}\Gamma$ than for the $G_{0}W_{0}^{\rm{tc-tc}}$ level of theory (see the finite basis quasi-particle levels and the extrapolated counterparts in Table 2). This observation can be rationalised with a similar diagrammatic analysis as carried out for the paramagnetic electron gas [69]. The argument can be concisely restated as follows: in lowest order there is a cancellation between the bubble and the exchange diagrams in the self-energy (see top and bottom row in Figure 2) because of the sign rule for fermionic loops that assigns a different sign to the two classes of diagrams. This cancellation is only partial in second order because there are twice more ways to assign a spin variable to each propagator line in the bubble diagram than there are for the exchange diagram (this is a consequence of the spin conservation at each interaction line). From these considerations the linear extrapolation of the quasi-particle energy will have a smaller slope by a factor $\frac{1}{2}$ if the exchange diagrams in the self-energy are included. The inclusion of higher order diagrams clearly partially invalidates the previous arguments, since already in third order there are exchange type diagrams that feature the same number of fermionic loops as the bubble diagram (see Figure 3 bottom right). In general terms, we can however still expect a faster convergence over the number of orbitals for the $G_{0}W_{0}\Gamma$ level of theory in comparison with the $G_{0}W_{0}^{\rm{tc-tc}}$ case. This turns out to be the case for the vast majority of the systems considered here. The faster convergence for the $G_{0}W_{0}\Gamma$ results can be seen in Figure 7 for CO and PH${}_{\textrm{3}}$.

Finally, we can compare the resulting basis set extrapolated quasi-particle energies with the CCSD(T) results provided in [9]. We reiterate that this comparison is somewhat problematic, since basis set incompleteness errors were not included in the coupled cluster results, and they are estimated to be over the order of 100 to 150 meV, as already elaborated above. Since CCSD(T) calculations also include exchange diagrams (as all interactions are antisymmetrised) we expect CCSD(T) to have a similar convergence rate as $G_{0}W_{0}\Gamma$ and not as $G_{0}W_{0}$. Neglecting these corrections, it is clear that the standard $G_{0}W_{0}$@HF yields much too large IP’s. The $G_{0}W_{0}^{\rm{tc-tc}}$ approximation significantly improves on these results, with only few outliers with an error exceeding 300 meV. These outliers are two polar molecules, with much too small IP, namely HF and LiH, as well as CS, CO and N${}_{\text{2}}$ where the IP is still way too large. However, the $G_{0}W_{0}^{\rm{tc-tc}}$ approximation is not a systematic approximation, since it includes vertex correction in the polarizability only but neglects them in the self-energy. If we keep in mind, that basis set converged CCSD(T) calculations will result in more negative HOMO energies (about $-100$ meV to $-150$ meV), this could possibly improve agreement between $G_{0}W_{0}\Gamma$ and CCSD(T) and worsen the $G_{0}W_{0}^{\rm{tc-tc}}$ results.

The effect of the vertex in the self-energy is to systematically bring the extrapolated IP again above the CCSD(T) reference values Figure 8. Generally, the $G_{0}W_{0}\Gamma$ IP’s lie between the $G_{0}W_{0}$ and $G_{0}W_{0}^{\rm{tc-tc}}$ values, which is in agreement with the general assumptions in literature. Exceptions to this trend are the unsaturated hydrocarbons, SiH${}_{\text{4}}$ and HCN, whose IP’s, accidentally, reproduce the CCSD(T) values quite well. It is interesting to comment on the biggest outliers in Figure 8. These are N${}_{\text{2}}$, CO and CS. A previous study [74] has established the strong dependence on the used orbitals for carbon monoxide and the nitrogen dimer. Indeed, if full self-consistency is achieved (with the update of the initial Hartree-Fock orbitals), the resulting IP’s are reduced by an excess of 0.8 eV. Unfortunately, similar studies have not been carried out for CS, however we expect a similar behaviour for this system, given the similarities in electronic structures. Therefore, for these molecules we expect that an update of the orbitals, or at least of the orbital’s energies, more than the inclusion of the vertex will be relevant for an estimate of their IP’s.

The mean absolute deviation from the coupled cluster values for the remaining molecules in the set is 0.21 eV for $G_{0}W_{0}\Gamma$ and 0.15 eV for $G_{0}W_{0}^{\rm{tc-tc}}$, scoring a marked improvement in comparison with the $G_{0}W_{0}$ deviation (0.46 eV). We believe that inclusion of basis set corrections in the CCSD(T) results would improve the agreement between $G_{0}W_{0}\Gamma$ and CCSD(T) even further, since this would increase the CCSD(T) IP values by the already quoted 100-150 meV moving them closer to the $G_{0}W_{0}\Gamma$. Obviously this correction would worsen the agreement between CCSD(T) and $G_{0}W_{0}^{\rm{tc-tc}}$, which is to be expected, since $G_{0}W_{0}^{\rm{tc-tc}}$ is not a well balanced diagrammatic theory. In Ref. [26] range-separated density functionals have been used as a starting point of the $G_{0}W_{0}$ calculation, and these scored a similar improvement in quasi-particle energies as the one reported here by applying the vertex. In some sense, this starting point is constructed to reproduce the system’s IP, whereas the method we employ is completely ab initio and, as such, does not require any previous information on the system’s electronic structure.

Another well established class of ab initio approaches is the algebraic diagrammatic construction (ADC). It offers a general framework to generate the self-energy diagrams at any order in perturbation theory. The two-particle one-hole Tamm-Dancoff approximation (2ph-TDA) is the simplest ADC method and it includes bubble and particle-hole ladder diagrams in its resummation. It is quite striking that, albeit this approximation is diagrammatically similar to the $GW\Gamma$ method, its performance for small molecules is unsatisfactory with a mean absolute error (vs CCSD(T)) that exceeds 0.6 eV [99]. The inclusion of all self-energy diagrams up to third order, as it is performed in the ADC(3) approximation, largely improves the IP’s values. It is important to point out that in the 2ph-TDA approach the 2p-2h virtual excitations are not included, whereas these effects are implicitly accounted for by the coupling amplitudes in the full ADC(3) [100]. These excitations are responsible for introducing correlation effects in the HF ground state [68] (which is obviously uncorrelated) and are explicitly built in the $GW$ and $GW\Gamma$ approximation thanks to the inclusion of the so-called resonant-antiresonant coupling. The comparison between our results and those obtained with the electron propagator methods cited above allows us to rationalise the poor performance of the 2ph-TDA. Specifically, the coupling between resonant and antiresonant virtual transistions is a key feature of $GW$ (and $GW\Gamma$) that is implicitly accounted for in ADC(3). This coupling seems to be responsible for a dramatic improvement over the quasi-particle energies obtained at the 2ph-TDA level. On the other hand, the inclusion of other higher order particle-particle ladders in the ADC(3) method appear to have an overall modest impact (with the exceptions discussed above that are poorly described by the HF reference), since its performance on small molecules is comparable to that of the $GW\Gamma$ method.

The comparison with experimental data (Figure 8, bottom panel) proceeds along similar lines as above, with the mean absolute error increasing slightly owing to the impact of adiabatic effects on the measurements for a number of systems considered (these are reported in Table 1 as non-italicised entries).

5 Conclusions

The main topic of this work is the inclusion of exchange diagrams in both the polarization propagator as well as the self-energy using Hedin’s equations. Commonly this is referred to as the $GW\Gamma$ approximation. In Hedin’s equations the interaction kernel is strictly given by the functional derivative of the self-energy employed in the preceding step with respect to the Green’s function. In the present work, we adopt the Hartree-Fock starting point. Then from the derivative of the Hartree-term $V_{\mathrm{H}}$ with respect to the Green’s function, one obtains the standard RPA bubble diagrams, and from the derivative of the exchange potential $\Sigma_{\mathrm{x}}=V_{\mathrm{x}}$ one obtains the ladder diagrams. Hence, applying the common $GW$ approximation— that includes bubbles only but neglects ladder diagrams —is from the outset a flawed prescription for the Hartree-Fock starting point. If it were successful, this could only be related to some sort of error cancellation.

An important aspect of the present work is whether such a cancellation exists. To finally and concisely answer this question, we have performed quasi-particle calculations on top of the Hartree-Fock approximation for a set of 29 small molecules. As already alluded to above, we have included the ladder diagrams related to the derivative of the exchange potential, in both the polarization propagator and the self-energy. The applied equations have been derived from Hedin’s equations and should be fairly easy to implement in any time-dependent Hartree-Fock, Casida or BSE code. We have also argued that, although the final equations (23) essentially only rephrase the equation of motion for the Green’s function, they have not yet been used in quantum chemistry or computational solid state physics. This is somewhat astounding, since the implementation is, in principle, straightforward. All one needs to know are the two-particle eigenstates of the TD-HF equation.

Now let us turn to our results. The first, already often published, observation is that the $G_{0}W_{0}$ approximation is not satisfactory when performed on top of HF. It systematically overestimates the IP’s. This should come as no surprise. In the $G_{0}W_{0}$ approximation, a non-physical polarizability is used, that is calculated without particle-hole ladder diagrams related to the variation of the exchange. Inclusion of the ladder diagrams in the polarization propagator (which is equivalent to using the polarizability from TD-HF) largely rectifies this issue and massively improves agreement with the CCSD(T) reference values. The mean absolute error decreases from about 400 meV in $G_{0}W_{0}$ to 150 meV in this $G_{0}W_{0}^{\mathrm{tc-tc}}$ approximation.

Including the exchange diagrams in the polarizability only is, however, not an entirely well balanced approximation instead, one should include the exchange diagrams also in the self-energy ($G_{0}W_{0}\Gamma$ approximation). The $G_{0}W_{0}\Gamma$ values are always in between the standard $G_{0}W_{0}$ and the $G_{0}W_{0}^{\mathrm{tc-tc}}$ results, as discussed in more detail in the next paragraph. With a mean absolute deviation of 210 meV compared to the CCSD(T) reference results, $G_{0}W_{0}\Gamma$ results are on average slightly worse than the $G_{0}W_{0}^{\mathrm{tc-tc}}$ (MAE 150 meV). This is certainly somewhat unsatisfactory and requires further studies and, in particular, better reference values. With an error of roughly 100 meV for cc-pVQZ basis sets, CCSD(T) values without basis set extrapolation are not sufficiently accurate to allow for an unambigious comparison. In passing, we also emphasize a side result of our study, namely that $G_{0}W_{0}$ converges differently with respect to basis set size than $G_{0}W_{0}\Gamma$. This is a result of the inclusion of the exchange diagrams in the self-energy as shown in Fig. 7. This also means that comparisons at finite basis sets (without basis set corrections) need to be done with caution, since one could observe fortuitous agreement at small basis sets that might not hold up at improved basis sets.

A key result of this study is that changes from $G_{0}W_{0}^{\mathrm{tc-tc}}$ to $G_{0}W_{0}\Gamma$ are very system-dependent: in some cases, $G_{0}W_{0}\Gamma$ recovers the $G_{0}W_{0}$ values, in other cases the final values are closer to $G_{0}W_{0}^{\mathrm{tc-tc}}$. As an example, one can consider the case of ammonia and phosphine: they have a comparable geometry and valence electron configuration, with the obvious difference that phosphorus is one period below nitrogen in the periodic table. In one case (PH${}_{\text{3}}$) the introduction of the vertex in the self-energy has a marginal impact on the IP, whereas for NH${}_{\text{3}}$ the higher electronegativity of the central atom, as discussed in the previous section, leads to increased electron correlations and to a more pronounced difference between $G_{0}W_{0}^{\mathrm{tc-tc}}$ and $G_{0}W_{0}\Gamma$. Albeit these observations are in line with ”chemical intuition”, they are difficult to generalise, let alone if quantitative estimates are required. We conclude that there is simply no shortcut to $GW\Gamma$ or general error cancellation; one needs to include the vertex consistently. All these results are in agreement with previous work on solids, where on average the vertex in the self-energy also compensated the effect of the vertex in the polarization propagator, and as for molecules, the results varied significantly between different solids and even different orbitals [44].

Finally, let us comment on the future perspectives: to us it is clear that the effect of self-consistency on the Green’s function needs to be studied. Unfortunately, this is not a simple matter. If one updates the self-energy and the Green’s function, additional terms come into play in the interaction kernel $I$, and ultimately one is then forced to make compromises between computational efficiency and accuracy [22]. We plan to report a combination of self-consistency and vertex corrections in our forthcoming study.

In summary, we have shown that $G_{0}W_{0}\Gamma$@HF improves upon $G_{0}W_{0}$@HF. This implies that the vertex in Hedin’s equations is important and can not be neglected. Given the simplicity of the present approach, we are confident that it will draw interest in the computational quantum chemistry community, as it constitutes an important fundamental step towards accurate quasi-particle calculations.

Acknowledgements

Funding by the Austrian Science Fund (FWF) within the Spezialforschungs- bereich F41 (SFB ViCoM) is gratefully acknowledged. Computations were predominantly performed on the Vienna Scientific Cluster VSC3.