Simple eigenvalue-self-consistent $Δ¯GW$

The GW approximation¹ to many-body perturbation theory is often used to calculate electron removal or addition energies and the related (inverse) photoemission spectra of molecules, nanostructures, and bulk materials.^2–15 GW is part of a family of methods that describe the probability amplitude of a quasiparticle (QP) to propagate between two space-time points $r,t$ and $r′,t′$ with Green’s function $Gr,r′,t,t′$⁠, the poles of which are the QP energies. Green’s function is obtained perturbatively from reference (non-interacting) Green’s function, $G0r,r′,t,t′$⁠, via a Dyson equation (all equations use atomic units),

where $Σr1,r2,t$ is the time-dependent self-energy. Reference Green’s function is typically^16,17 given by the Kohn-Sham¹⁸ (KS) density function theory (DFT).¹⁹ In GW, the self-energy Σ is approximated as

where $Wr,r′,t$ is the screened Coulomb interaction, usually evaluated within the random phase approximation (RPA).

A solution of the equations above requires in principle a self-consistent procedure since both $Σr,r′,t$ and $Wr,r′,t+$ depend on $Gr,r′,t$⁠. In practice, the self-consistency is often abandoned and the most common GW treatment is based on “one-shot” scheme,^16,17 G₀W₀, where the right hand side of Eq. (2) becomes $G0r,r′,tW0r,r′,t$ and W₀ is obtained from a random phase approximation that uses the KS eigenstates. The one-shot approach improves significantly the KS-DFT results, yet it depends on the choice of the reference system and often underestimates the QP gaps (E_g) and the ionization potentials (I).^8,20–25 A fully self-consistent solution is computationally extremely demanding,^{10,21,26–33} and in many situations, it yields results that are worse than G₀W₀;^8,34,35 heuristically, since GW is an approximation, an iteration of that approximation may enhance the deviation from the exact result.

To simplify the problem, a static and Hermitian approximation to the self-energy^20,36 is sometimes used in the so-called QP self-consistent GW (qpGW), which iteratively updates the Σ and QP wave-functions (i.e., the Dyson orbitals). The qpGW approach is still computationally expensive and cannot be applied for large systems; furthermore, it tends to overestimate E_g^35–37 and the ionization potentials^15,38,39 due to an overly strong screened Coulomb interaction.^15,35,36 Alternatively, the W term is kept “frozen” and self-consistency is sought only in Green’s function.^{8,27,28,35,38} This method is termed eigenvalue self-consistent GW (evGW₀), and it was applied successfully to bulk systems^8,35,40 and to organic molecules.^25,27,28 Despite the lack of self-consistency in W, which leads to overscreening (since it is based on an underestimated KS DFT gap)³⁵ and non-negligible starting point dependence,^41,42 evGW₀ provides a reasonable route for improving the “one-shot” predictions or gaps in extended systems, notwithstanding the fact that for fairly small molecular systems the starting-point dependence is very large.

While overall self-consistency methods provide an improvement over one-shot G₀W₀, there is an issue of scaling and overall cost. Even one-shot G₀W₀ has usually been quite expensive, scaling as $O(Ne3)−O(Ne4)$ with the number of electrons N_e, and self-consistency methods at least multiply the cost by the number of iterations. Thus, in spite of their potential, G₀W₀ and its self-consistent improvements were only applicable to fairly small systems unless drastically loose convergence criteria are made and/or huge computational efforts are used, and even then large systems with many thousands of electrons were not feasible.

Recently, however, we developed an effectively linear-scaling stochastic approach to G₀W₀ with moderate computational cost.^43–45 The approach has been verified to give the same results as deterministic converged-basis (or converged plane-wave) G₀W₀ simulations and can handle systems with more than 10 000 electrons.⁴⁵ But as stochastic G₀W₀ obtains directly only a column of the self-energy matrix rather than the full-matrix, most self-consistent improvements cannot be applied to it.

Here, we enhance the power of the stochastic approach, by showing that a variant of evGW₀ whereby the energies of the HOMO and LUMO are shifted rigidly (a scissors approach, labeled below $Δ¯GW0$⁠) can be easily applied on top of G₀W₀, without any additional cost. Thus stochastic GW can be enhanced beyond G₀W₀; the combined approach is labeled stochastic $Δ¯GW0$⁠.

Specifically, we first develop the method showing how the application of a scissors operator for Green’s function is an a posteriori step. The $Δ¯GW0$ formulation requires only the knowledge and use of two quasiparticle states, rather than the whole quasiparticle spectrum. As long as one has access to the relevant two diagonal matrix elements (usually the HOMO and LUMO elements) of the self-energy at all frequencies or all times, then, irrespective of the system size, the computational cost of the self-consistency step is negligible (i.e., seconds on a single-core machine), so $Δ¯GW0$ costs no more than G₀W₀.

Equally important, we show that stochastic $Δ¯GW0$ improves with increasing system size. First, we show for large molecules that the ionization energies from $Δ¯GW0$ converge to the evGW₀ predictions and both methods approach for large systems the more exact coupled-cluster single double triple (CCSD(T)) results.

Next, we perform stochastic G₀W₀ calculations for periodic semiconductors and insulators using large supercells with 6192 valence electrons. For solids, $Δ¯GW0$ consistently increases E_g toward better agreement with evGW₀ and the experiment (with resulting mean absolute error relative to the experiment of better than 0.1 eV). In all cases, the self-consistency is reached in very few iterations without any additional cost beyond the G₀W₀ calculation.

The paper is organized as follows: Sec. II presents the theory, Sec. III presents the results for molecules and solids, and conclusions follow in Sec. IV.

The QP energy of the ith state $εQPi$ is calculated using the usual form of the perturbative GW approximation in which the Kohn-Sham eigenvalues (ε⁰) are corrected by the QP shift (Δ) using a fixed point equation,

where

Here, $vxc$ is the Kohn-Sham exchange-correlation potential for the DFT density, $Σ̃iω$ is the Fourier transform of the matrix-element of $Σ^(t)$⁠,

and $Σ^(t)$ is given by Eq. (2).

Starting from a KS DFT reference point, the initial self-energy is constructed from the KS propagator

where Tr denotes a trace over all KS states, μ is the chemical potential, θ is the Heaviside step function that guarantees forward and backward time propagation for particles and holes, respectively, and ĥ⁰ is the KS Hamiltonian,

where we introduced the kinetic energy and the external and Hartree potentials. In the rest of the paper, we employ real time-dependent Hartree propagation to calculate the screened Coulomb interaction;^43,44 this is equivalent to using the RPA approximation for W₀.

In the time-domain, the self-energy matrix element for the ith state is

Finally, after Fourier transformation combined with time-ordering,^43,44 the “one-shot” QP energy is calculated through Eq. (4).

In the evGW₀ procedure, Green’s function is reconstructed in each iteration, employing the QP energies from the previous iteration. In the time domain, this corresponds to writing the propagator as

where $Δ^$ contains all the many-body contributions to the spectrum.

As common in evGW₀ self-consistency,^8,35,38 for the purpose of Eq. (9), one ignores the non-Hermiticity of the self-energy and treats it as diagonal so that in the KS basis the correction operator is approximated as

Hence, all the KS energies in the exponent in Eq. (9) are shifted to the QP energies obtained from Eq. (3). Equation (10) is the fundamental equation of evGW₀.

Since $Δ^$ in Eq. (9) is by construction diagonal in the KS basis set, we view it as a function of the KS Hamiltonian $Δ¯ĥ0$ that interpolates all QP shifts. Green’s function is therefore

This simple expression allows for a further approximation described below that significantly reduces the computational cost associated with self-consistent treatment.

In many cases, $Δ¯$ is well described by a low degree polynomial with discontinuity only at the bandgap energies, $εH$ and $εL,$ corresponding to the highest occupied (H) and lowest unoccupied (L) states, respectively. The zeroth order term in this polynomial is a scissors operator^46,47 which shifts the occupied and unoccupied states down and up in energy, respectively,

We call this approximation $Δ¯GW0$ and use it in Sec. III for molecules and periodic systems. Furthermore, rather than choosing highest occupied and lowest unoccupied states, $ReΔ¯$ may be taken based on an average shift of all occupied and unoccupied states, respectively.

Combining Eqs. (12) and (11) leads to modified Green’s function which acquires an additional phase shift that is different for positive and negative times. Therefore, in the time domain, we define

In each iteration, the updated self-energy matrix element is then calculated as

Next, the self-energy matrix element is transformed to the frequency domain and used in Eq. (3) to calculate a new estimate of the QP energy. The new QP energy is used iteratively to update Eqs. (12) and (13). The full cycle is illustrated in Fig. 1.

Note that this form of self-consistency is trivial and is a postprocessing step with no additional cost, unlike previous uses^46,47 of the scissors-operator in GW which require repeated evaluations of the self-energy. Furthermore, the $Δ¯GW0$ approach is applicable to any implementation which yields $Σ̃iω$ for the HOMO and LUMO states. Thus, this method is naturally suited for the stochastic G₀W₀ method^43–45 which provides the self-energy over the full-time domain and therefore over a wide range of frequencies (spanning several hundred eV). Note that while in principle one could use higher order polynomials of Δ, and such polynomials would generally require going beyond a post-processing method, i.e., would require additional GW calculations.⁴⁸ Therefore, the simple equation (14), for a scissors-type expression, is our main result.

We first test our approach on ionization potentials I [taken as the negative of the HOMO energy, −ε(H)] for a set of small molecules listed in Table I. A ground state DFT calculation is performed using a Fourier real-space grid, ensuring (through the Martyna-Tuckerman approach)⁴⁹ that the potentials are not periodic. The exchange-correlation interaction is described by the local density approximation (LDA)⁵⁰ with Troullier-Martins pseudopotentials;⁵¹ the DFT eigenvalues are converged up to <10 meV with respect to the spacings of the real space grids (given in Table I).

Six molecules were tested, ranging in size from N₂ to hexacene. In all cases, stochastic G₀W₀^43–45 was used to calculate the self-energy. $Δ¯GW0$ converges in 3-4 iterations after the initial G₀W₀ calculation; the self-energy curves are illustrated for hexacene in Fig. 2. Since our self-consistency procedure is performed a posteriori, its computational cost is negligible (less than a second on a single core machine). In all cases, $Δ¯GW0$ increases the ionization potentials above the one-shot values. This is a desired outcome since G₀W₀ underestimates I.^24,42,44

Table I shows the ionization energies at several levels, including G₀W₀, $Δ¯GW0$⁠, deterministic self-consistent GW on two levels, CCSD(T), and also the experiment.

The G₀W₀ column shows our stochastic G₀W₀ values that are fully converged with grid size (i.e., in a plane-wave language, corresponding to a fully converged plane-wave cutoff) and are therefore close to the converged deterministic plane-wave basis results benchmarked in Ref. 52 as well as basis set extrapolated values obtained with deterministic codes, with differences of typically less than 0.1 eV, which are mostly due to the use of an LDA functional rather than PBE exchange-correlation functional.⁴⁴ 

Our partially self-consistent $Δ¯GW0$ is then compared in the table to other partially self-consistent techniques in the literature, suitable for small molecules. The most direct comparison should be to evGW₀, where, as mentioned, Green’s function is also converged (leaving W₀ frozen at its initial value), for an initial Hamiltonian with a local potential. There are large deviations for the small molecules (0.6–0.7 eV), and there are several contributing factors to these differences. First, for the three small molecules about 0.1-0.2 eV difference is due to the non-completeness of the basis set (see Ref. 24 for the effect of basis-set extrapolation). For nitrogen and urea, another factor of ≈0.1 eV is due to our use of LDA rather the PBE in the quoted evGW₀ results.⁴⁴ Taken together, this still leaves about an error of ≈0.4 eV between our results and evGW₀ for the small molecules.

The difference between the present $Δ¯GW0$ and evGW₀ is made more explicit by the last two columns in Table I, which shows the degree of correction due to self-consistency (i.e., the difference between the final result and the one-shot G₀W₀) for our calculation and evGW₀. For the small molecules, $Δ¯GW0$ gives a correction which is higher than evGW₀ by more than 0.4 eV. But for larger systems, the deviation is getting much smaller. Specifically, the deviation from $Δ¯GW0$ decreases to 0.1–0.25 eV for the three larger molecules (three acenes) where the deterministic basis-set results for the acenes were already extrapolated to the infinite basis-set limit. A small portion of the difference is due to our use of LDA rather than PBE. Overall the much better agreement of the rigid-shift approach with evGW₀ for large systems is due to the improved assumption of a rigid shift.

A separate question from the agreement of the rigid-shift approximation with evGW₀ is how accurate the overall method is, i.e., does the method agree with more exact non-GW approaches. Here, the comparison with the CCSD(T) ionization energy (extrapolated to the limit of infinite basis-set) is quite encouraging. For the small molecules, there is a large deviation from CCSD(T), but it decreases dramatically to 0.15 eV and below once we go to the larger acenes.

Interestingly, $Δ¯GW0$ is generally closer to the CCSD(T) results than evGW₀. This is probably due to error cancellation; i.e., evGW₀ based on a local potential somewhat underestimates the CCSD(T) results, and $Δ¯GW0$ generally somewhat overestimates evGW₀ due to the rigid shift assumption. Nevertheless, Table I does show that for the ionization energy, as the system size grows the essentially exact CCSD(T) is better and better predicted by both the evGW₀ approximation and our $Δ¯GW0$ further approximation to evGW₀. Likewise, $Δ¯GW0$ agrees well with the experiment for the larger molecules (the three acenes), which makes sense since except for the effect of vibrational energy change CCSD(T) should reproduce the experiment.

To understand the improvement of ionization energies from $Δ¯GW0$ with the system size, we examine the dependence of the state dependence shift, Δ(i), on the state energy ε⁰(i). The large errors for small molecules stem from the strong dependence of Δ(i) on ε⁰. Namely, the QP shift is far from being constant for all occupied (or unoccupied) states, so Eq. (12) is not fulfilled. For the smallest molecule we tested, N₂, G₀W₀ shifts the lowest valence state by −7.34 ± 0.06 eV, while the HOMO energy is decreased on this level by only −4.63 ± 0.05 eV; i.e., the difference in QP shifts across the valence band is very high, 2.60 eV. But for larger systems, the scissors operator approximation in Eq. (12) is getting gradually better. Thus, for hexacene, the difference in QP shifts between the bottom valence and the HOMO state is only 0.87 eV, three times smaller than that for N₂.⁵³ 

The improvement of the rigid-shift approximation with the system size is important since it implies that for large system sizes, for which other self-consistent techniques are too expensive, the $Δ¯GW0$ approach would be sufficiently accurate.

We conclude this section by a cautionary note. The success of both evGW₀ based on a local potential and $Δ¯GW0$ is for the ionization potential. For the LUMO, i.e., for the electron affinity, the results of both evGW₀ and $Δ¯GW0$ are not as good (deviations of more than 0.5 eV from the exact results even for the larger systems). For the LUMO, a time dependent DFT (TDDFT)-based kernel may give better one-shot results,^43,54 and the application of rigid-shift self-consistency for the LUMO would be described in a future publication.

Next we study self-consistency for periodic solids, using an extension of the stochastic GW approach to treat periodic boundary conditions.⁴⁵ The self-consistent $Δ¯GW0$ method is demonstrated on large supercells with 1712 atoms (corresponding, for the system we study below, to 6 × 6 × 6 conventional cells with 6912 valence electrons). The stochastic approach handles well such large system sizes,^43–45 yielding again self-energies for selected electronic states. Here we are interested in the fundamental bandgap, E_g = ε(L) − ε(H), so we need both the highest occupied and lowest unoccupied states. For large extended systems, the complications associated with a TDDFT vs. RPA based kernel are tamed in comparison with molecules, so it makes sense to study the gap, based solely on RPA.⁵⁷ 

Table II shows for several solids the one-shot and $Δ¯GW0$ and compares them again to previous self-consistent GW calculations and to experiment. The results demonstrate first that the one-shot correction yields bandgaps that are mostly (though not always) lower than experimental values, in agreement with previous calculations.^8,20,35,36

Moving to $Δ¯GW0$⁠, in all cases studied, self-consistency is quickly achieved within 2 or 3 iterations and the self-energy curves are illustrated in Fig. 3. The fundamental bandgaps are enlarged by as much as 0.17 eV compared to G₀W₀. For all the solids in Table II, the difference between the QP shifts for the bottom and top valence states is small (correlating slightly with E_g) and similar in magnitude for all states.⁵⁸ The underlying assumption that all states in either the occupied or unoccupied subspace are shifted (approximately) by a constant is fulfilled well.

Gaps predicted by $Δ¯GW0$ are similar to previously published evGW results.⁵⁹ Much of the difference is due to the use of distinct atomic potentials in previous calculations, which are known to cause energy differences up to ∼0.2 eV.⁶⁰ More importantly for our purpose, here are the energy differences due to the self-consistency (i.e., the differences of either evGW₀ or $Δ¯GW0$ energies from G₀W₀). These differences (which are less sensitive to the particular type of potentials used) are very close, to within about 0.04 eV in average.

In Table II, we also list experimental bandgaps. Such comparison has to be performed with caution since our calculations do not account for renormalization of E_g due to zero-point motion,⁶¹ which decreases the gap by as anywhere from 0.05 eV (silicon) to 0.37 eV (diamond).⁶² Yet it is encouraging to see that $Δ¯GW0$ enlarges E_g and that our predictions are quite close to experimental values.

We derived here a general form of Green’s function self-consistency in the time domain and introduced its simplified form, which we label $Δ¯GW0$⁠. The starting point is the scissors correction, where the differences between Kohn-Sham eigenvalues and quasiparticle energies are approximately just two constants, one for occupied and the other for unoccupied states. We approximate this scissors-like correction by the corrections to the HOMO and LUMO energies. But unlike previous self-consistency schemes (even those relying on scissors-like operators), our approach does not require repeated many-body calculation and uses only two self-energy matrix elements. Practically, it is merely an a posteriori treatment of the time-dependent self-energy matrix. Yet, it captures the major changes in the structure of Green’s function and hence of the self energy for large systems.

$Δ¯GW0$ has essentially no additional computational cost beyond that of a one-shot G₀W₀ calculation. In conjunction with the nearly linear scaling stochastic G₀W₀, it can treat extremely large systems with thousands of electrons. The combined method is best labeled as stochastic $Δ¯GW0$⁠.

We tested stochastic $Δ¯GW0$ on molecules and on periodic semiconductors and insulators with large periodic supercells with 6192 electrons. The predicted ionization potentials and fundamental bandgaps are overall better than one-shot G₀W₀ values when compared to high-level methods and/or experiments. Our simplified self-consistency treatment is especially appropriate for large molecules and periodic systems, for which no other computationally affordable approach existed up to now.

$Δ¯GW0$ is an approximation on top of another, successful approximation, evGW₀. For small systems, both approximations have faults (dependence on initial DFT Hamiltonian in the final results of evGW₀, compounded by strong state dependence of the self-energy diagonal elements violating the underlying assumption on $Δ¯GW0$⁠). But the purpose of the new approach is to address very large systems, for which these two issues are much less severe, and for which $Δ¯GW0$ is automatically suited, since, unlike even evGW₀, it does not require the computation of new operators when the quasi-particle energies are updated.

The stochastic partially self-consistent $Δ¯GW0$ approach presented here is both accurate for large systems and efficient, opening the door to many future applications in chemistry, physics, and nano- and material sciences.

V.V., E.R., and D.N. were supported by the Center for Computational Study of Excited State Phenomena in Energy Materials (C2SEPEM) at the Lawrence Berkeley National Laboratory, which is funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division, under Contract No. DEAC02-05CH11231 as part of the Computational materials Sciences Program. R.B. is grateful for support from the Binational Science Foundation, under Grant No. 2015687, and for support from the Israel Science Foundation—FIRST Program, under Grant No. 1700/14. The calculations were performed as part of XSEDE computational Project No. TG-CHE170058.⁶⁸