We show that a rigid scissors-like GW self-consistency approach, labeled here Δ¯GW0, can be trivially implemented at zero additional cost for large scale one-shot G0W0 calculations. The method significantly improves one-shot G0W0 and for large systems is very accurate. Δ¯GW0 is similar in spirit to evGW0 where the self-consistency is only applied on the eigenvalues entering Green’s function, while both W and the eigenvectors of Green’s function are held fixed. Δ¯GW0 further assumes that the shift of the eigenvalues is rigid scissors-like so that all occupied states are shifted by the same amount and analogously for all the unoccupied states. We show that this results in a trivial modification of the time-dependent G0W0 self-energy, enabling an a posteriori self-consistency cycle. The method is applicable for our recent stochastic-GW approach, thereby enabling self-consistent calculations for giant systems with thousands of electrons. The accuracy of Δ¯GW0 increases with the system size. For molecules, it is up to 0.4-0.5 eV away from coupled-cluster single double triple (CCSD(T)), but for tetracene and hexacene, it matches the ionization energies from both CCSD(T) and evGW0 to better than 0.05 eV. For solids, as exemplified here by periodic supercells of semiconductors and insulators with 6192 valence electrons, the method matches evGW0 quite well and both methods are in good agreement with the experiment.

The GW approximation1 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–15GW 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),

Gr,r,t,t=G0r,r,t,t+G0r,r1,t,t1×Σr1,r2,t1t2Gr2,r,t2,tdr1dt1dr2dt2,
(1)

where Σr1,r2,t is the time-dependent self-energy. Reference Green’s function is typically16,17 given by the Kohn-Sham18 (KS) density function theory (DFT).19 In GW, the self-energy Σ is approximated as

Σr,r,t=iGr,r,tWr,r,t+,
(2)

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,17G0W0, where the right hand side of Eq. (2) becomes G0r,r,tW0r,r,t and W0 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 (Eg) 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 G0W0;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-energy20,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 Eg35–37 and the ionization potentials15,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 (evGW0), and it was applied successfully to bulk systems8,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)35 and non-negligible starting point dependence,41,42 evGW0 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 G0W0, there is an issue of scaling and overall cost. Even one-shot G0W0 has usually been quite expensive, scaling as O(Ne3)O(Ne4) with the number of electrons Ne, and self-consistency methods at least multiply the cost by the number of iterations. Thus, in spite of their potential, G0W0 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 G0W0 with moderate computational cost.43–45 The approach has been verified to give the same results as deterministic converged-basis (or converged plane-wave) G0W0 simulations and can handle systems with more than 10 000 electrons.45 But as stochastic G0W0 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 evGW0 whereby the energies of the HOMO and LUMO are shifted rigidly (a scissors approach, labeled below Δ¯GW0) can be easily applied on top of G0W0, without any additional cost. Thus stochastic GW can be enhanced beyond G0W0; 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 G0W0.

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 evGW0 predictions and both methods approach for large systems the more exact coupled-cluster single double triple (CCSD(T)) results.

Next, we perform stochastic G0W0 calculations for periodic semiconductors and insulators using large supercells with 6192 valence electrons. For solids, Δ¯GW0 consistently increases Eg toward better agreement with evGW0 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 G0W0 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 (ε0) are corrected by the QP shift (Δ) using a fixed point equation,

εQPi=ε0i+Δi,
(3)

where

Δi=Σ̃iω=εQPiϕivxcϕi.
(4)

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

Σ̃iω=ϕi|Σ^t|ϕieiωtdt,
(5)

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

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

iG0r,r,t=Trrreiĥ0tθtθβĥ0μθtθβμĥ0,
(6)

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 ĥ0 is the KS Hamiltonian,

ĥ0=122+vext+vH+vxc,
(7)

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 W0.

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

ϕiΣ^0tϕi=iϕirG0r,r,tW0r,r,t+ϕirdrdr.
(8)

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

In the evGW0 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

iGr,r,t=Trrreiĥ0+Δ^t×θtθβĥ0μθtθβμĥ0,
(9)

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

As common in evGW0 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

Δ^iϕiReΔiϕi.
(10)

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 evGW0.

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

iGr,r,t=Trrreiĥ0+Δ¯ĥ0t×θtθβĥ0μθtθβμĥ0.
(11)

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 operator46,47 which shifts the occupied and unoccupied states down and up in energy, respectively,

ReΔ¯ε0iΔHε0iε0HΔLε0iε0L.
(12)

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

Δ¯tΔHt<0ΔLt>0.
(13)

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

ϕi|Σ^t|ϕi=eiΔ¯ttϕi|Σ^0t|ϕi.
(14)

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.

FIG. 1.

The self-consistency cycle (full arrows). In the first step, G0W0 is used to calculate the self-energy Σ0t and the corresponding QP energies εQP for the HOMO and LUMO states (dashed arrow). The shift of the occupied and unoccupied states (Δ¯) is calculated from the QP HOMO and LUMO energies through Eq. (12). The updated time-dependent self-energy Σt=eiΔ¯tΣ0t is obtained via Eq. (13), which leads to new QP HOMO/LUMO energies. The cycle is repeated several times until reaching self-consistency.

FIG. 1.

The self-consistency cycle (full arrows). In the first step, G0W0 is used to calculate the self-energy Σ0t and the corresponding QP energies εQP for the HOMO and LUMO states (dashed arrow). The shift of the occupied and unoccupied states (Δ¯) is calculated from the QP HOMO and LUMO energies through Eq. (12). The updated time-dependent self-energy Σt=eiΔ¯tΣ0t is obtained via Eq. (13), which leads to new QP HOMO/LUMO energies. The cycle is repeated several times until reaching self-consistency.

Close modal

Note that this form of self-consistency is trivial and is a postprocessing step with no additional cost, unlike previous uses46,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 G0W0 method43–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.48 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)49 that the potentials are not periodic. The exchange-correlation interaction is described by the local density approximation (LDA)50 with Troullier-Martins pseudopotentials;51 the DFT eigenvalues are converged up to <10 meV with respect to the spacings of the real space grids (given in Table I).

TABLE I.

Ionization potentials I for small molecules by different methods. The G0W0 and Δ¯GW0 estimates were obtained using an LDA starting point and calculated through a stochastic implementation,43,44 and the statistical errors are reported in parentheses; evGW0 and qpGW refer to eigenvalue self-consistency and full self consistency, respectively, based on a PBE starting point and for the first three molecules (the non-acenes) are not extrapolated to the infinite basis-set limit. Experimental values are from Ref. 56. For the acenes, the CCSD(T) results are estimates extrapolated to the infinite basis-set limit.42 The last two columns show the difference between the partially self consistent values (for the present Δ¯GW0 and for evGW0) and the starting G0W0 energies. Molecular geometries were taken from Refs. 56 and 42.

I (eV)
Systemh (a0)LDAG0W0Δ¯GW0evGW0qpGWCCSD(T)Expt.δ(Δ¯GW0)δ(evGW0)
Nitrogen 0.35 10.44 15.08 (0.05) 15.93 (0.05) 15.32a 16.01a 15.57b 15.58 0.85 0.43 
Ethylene 0.35 6.92 10.50 (0.04) 10.87 (0.04) 10.24a 10.63a 10.67b 10.68 0.37 −0.09 
Urea 0.30 6.10 9.53 (0.08) 10.48 (0.08) 9.81a 10.45a 10.05b 10.28 0.95 0.49 
Naphthalene 0.35 5.71 8.10 (0.09) 8.39 (0.09) 8.15c … 8.25c 8.14 0.29 0.19 
Tetracene 0.35 4.89 6.79 (0.08) 6.94 (0.08) 6.84c … 7.02c 6.97 0.15 0.16 
Hexacene 0.35 4.52 6.15 (0.06) 6.33 (0.06) 6.19c … 6.32c 6.33 0.18 0.15 
I (eV)
Systemh (a0)LDAG0W0Δ¯GW0evGW0qpGWCCSD(T)Expt.δ(Δ¯GW0)δ(evGW0)
Nitrogen 0.35 10.44 15.08 (0.05) 15.93 (0.05) 15.32a 16.01a 15.57b 15.58 0.85 0.43 
Ethylene 0.35 6.92 10.50 (0.04) 10.87 (0.04) 10.24a 10.63a 10.67b 10.68 0.37 −0.09 
Urea 0.30 6.10 9.53 (0.08) 10.48 (0.08) 9.81a 10.45a 10.05b 10.28 0.95 0.49 
Naphthalene 0.35 5.71 8.10 (0.09) 8.39 (0.09) 8.15c … 8.25c 8.14 0.29 0.19 
Tetracene 0.35 4.89 6.79 (0.08) 6.94 (0.08) 6.84c … 7.02c 6.97 0.15 0.16 
Hexacene 0.35 4.52 6.15 (0.06) 6.33 (0.06) 6.19c … 6.32c 6.33 0.18 0.15 
a

Reference 39.

b

Reference 55.

c

Reference 42.

Six molecules were tested, ranging in size from N2 to hexacene. In all cases, stochastic G0W043–45 was used to calculate the self-energy. Δ¯GW0 converges in 3-4 iterations after the initial G0W0 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 G0W0 underestimates I.24,42,44

FIG. 2.

Graphical solution to the QP equation for the IP of hexacene using the self-energies from G0W0 and Δ¯GW0. The gray line is the diagonal frequency, and the arrows indicate solutions to the fixed point equation [Eq. (3)] that are the QP energies. The Δ¯GW0 result is visibly shifted toward a lower (more negative) energy than the one-shot solution.

FIG. 2.

Graphical solution to the QP equation for the IP of hexacene using the self-energies from G0W0 and Δ¯GW0. The gray line is the diagonal frequency, and the arrows indicate solutions to the fixed point equation [Eq. (3)] that are the QP energies. The Δ¯GW0 result is visibly shifted toward a lower (more negative) energy than the one-shot solution.

Close modal

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

The G0W0 column shows our stochastic G0W0 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.44 

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 evGW0, where, as mentioned, Green’s function is also converged (leaving W0 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 evGW0 results.44 Taken together, this still leaves about an error of ≈0.4 eV between our results and evGW0 for the small molecules.

The difference between the present Δ¯GW0 and evGW0 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 G0W0) for our calculation and evGW0. For the small molecules, Δ¯GW0 gives a correction which is higher than evGW0 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 evGW0 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 evGW0 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 evGW0. This is probably due to error cancellation; i.e., evGW0 based on a local potential somewhat underestimates the CCSD(T) results, and Δ¯GW0 generally somewhat overestimates evGW0 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 evGW0 approximation and our Δ¯GW0 further approximation to evGW0. 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 ε0(i). The large errors for small molecules stem from the strong dependence of Δ(i) on ε0. 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, N2, G0W0 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 N2.53 

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 evGW0 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 evGW0 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.45 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, Eg = ε(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.57 

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

TABLE II.

Fundamental bandgaps (Eg) for a sample of solids. For each system, a 6 × 6 × 6 conventional super-cell (with 6192 valence electrons) was used.

Eg (eV)
Systemh (a0)LDAG0W0 (present)G0W0aΔ¯GW0evGWaqpGWaExpt.δ(Δ¯GW0)δ(evGW0)
Si 0.446 0.50 1.24 (0.04) 1.12 1.32 (0.04) 1.20 1.28 1.17b 0.08 0.08 
SiC 0.293 1.36 2.40 (0.04) 2.27 2.50 (0.04) 2.43 2.64 2.42c 0.10 0.16 
AlP 0.368 1.46 2.45 (0.04) 2.44 2.58 (0.04) 2.59 2.77 2.51d 0.13 0.15 
0.336 4.15 5.60 (0.04) 5.50 5.71 (0.04) 5.68 5.99 5.48e 0.11 0.18 
BN 0.380 4.35 6.11 (0.04) 6.10 6.28 (0.04) 6.35 6.73 6.1–6.4f 0.17 0.25 
Eg (eV)
Systemh (a0)LDAG0W0 (present)G0W0aΔ¯GW0evGWaqpGWaExpt.δ(Δ¯GW0)δ(evGW0)
Si 0.446 0.50 1.24 (0.04) 1.12 1.32 (0.04) 1.20 1.28 1.17b 0.08 0.08 
SiC 0.293 1.36 2.40 (0.04) 2.27 2.50 (0.04) 2.43 2.64 2.42c 0.10 0.16 
AlP 0.368 1.46 2.45 (0.04) 2.44 2.58 (0.04) 2.59 2.77 2.51d 0.13 0.15 
0.336 4.15 5.60 (0.04) 5.50 5.71 (0.04) 5.68 5.99 5.48e 0.11 0.18 
BN 0.380 4.35 6.11 (0.04) 6.10 6.28 (0.04) 6.35 6.73 6.1–6.4f 0.17 0.25 
a

Reference 8.

b

Reference 63.

c

Reference 64.

d

Reference 65.

e

Reference 66.

f

Reference 67.

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 G0W0. 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 Eg) and similar in magnitude for all states.58 The underlying assumption that all states in either the occupied or unoccupied subspace are shifted (approximately) by a constant is fulfilled well.

FIG. 3.

Graphical solution to the QP equation for the G0W0 and Δ¯GW0 self-energies, for a diamond solid simulated by a 6 × 6 × 6 supercell. The gray line denotes the frequency, and the intersections with the black and red lines are the solution to the fixed point equation [Eq. (3)]. The full and dashed lines are for the top valence and bottom conduction states, respectively. Both axes are shifted so that zero is associated with the original LDA top of the valence band.

FIG. 3.

Graphical solution to the QP equation for the G0W0 and Δ¯GW0 self-energies, for a diamond solid simulated by a 6 × 6 × 6 supercell. The gray line denotes the frequency, and the intersections with the black and red lines are the solution to the fixed point equation [Eq. (3)]. The full and dashed lines are for the top valence and bottom conduction states, respectively. Both axes are shifted so that zero is associated with the original LDA top of the valence band.

Close modal

Gaps predicted by Δ¯GW0 are similar to previously published evGW results.59 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.60 More importantly for our purpose, here are the energy differences due to the self-consistency (i.e., the differences of either evGW0 or Δ¯GW0 energies from G0W0). 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 Eg due to zero-point motion,61 which decreases the gap by as anywhere from 0.05 eV (silicon) to 0.37 eV (diamond).62 Yet it is encouraging to see that Δ¯GW0 enlarges Eg 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 G0W0 calculation. In conjunction with the nearly linear scaling stochastic G0W0, 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 G0W0 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, evGW0. For small systems, both approximations have faults (dependence on initial DFT Hamiltonian in the final results of evGW0, 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 evGW0, 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.68 

2.
F.
Aryasetiawan
and
O.
Gunnarsson
,
Rep. Prog. Phys.
61
,
237
(
1998
).
3.
M. M.
Rieger
,
L.
Steinbeck
,
I.
White
,
H.
Rojas
, and
R.
Godby
,
Comput. Phys. Commun.
117
,
211
(
1999
).
4.
L.
Steinbeck
,
A.
Rubio
,
L.
Reining
,
M.
Torrent
,
I. D.
White
, and
R. W.
Godby
, “
Enhancements to the GW space-time method
,”
Comput. Phys. Commun.
125
,
105
(
2000
).
5.
G.
Onida
,
L.
Reining
, and
A.
Rubio
,
Rev. Mod. Phys.
74
,
601
(
2002
).
6.
A.
Rubio
and
S. G.
Louie
, “
Quasiparticle and optical properties of solids and nanostructures: The GW-BSE approach
,” in
Handbook of Materials Modeling
, edited by
S.
Yip
(
Springer
,
Dordrecht, New York
,
2005
), p.
215
.
7.
C.
Friedrich
and
A.
Schindlmayr
,
Computational Nanoscience: Do It Yourself!
in
NIC Series
, edited by
J.
Grotendorst
,
S.
Blügel
, and
D.
Marx
(
John von Neumann Institute for Computing
,
2006
), Vol. 31, pp.
335–355
.
8.
M.
Shishkin
and
G.
Kresse
,
Phys. Rev. B
75
,
235102
(
2007
).
9.
P. E.
Trevisanutto
,
C.
Giorgetti
,
L.
Reining
,
M.
Ladisa
, and
V.
Olevano
,
Phys. Rev. Lett.
101
,
226405
(
2008
).
10.
C.
Rostgaard
,
K. W.
Jacobsen
, and
K. S.
Thygesen
,
Phys. Rev. B
81
,
085103
(
2010
).
11.
I.
Tamblyn
,
P.
Darancet
,
S. Y.
Quek
,
S. A.
Bonev
, and
J. B.
Neaton
,
Phys. Rev. B
84
,
201402
(
2011
).
12.
M.
van Setten
,
F.
Weigend
, and
F.
Evers
,
J. Chem. Theory Comput.
9
,
232
(
2012
).
13.
G.
Stefanucci
and
R.
van Leeuwen
,
Nonequilibrium Many-Body Theory of Quantum Systems: A Modern Introduction
(
Cambridge University Press
,
2013
).
14.
M.
Govoni
and
G.
Galli
,
J. Chem. Theory Comput.
11
,
2680
(
2015
).
15.
F.
Kaplan
,
M. E.
Harding
,
C.
Seiler
,
F.
Weigend
,
F.
Evers
, and
M. J.
van Setten
,
J. Chem. Theory Comput.
12
,
2528
(
2016
).
16.
M. S.
Hybertsen
and
S. G.
Louie
,
Phys. Rev. Lett.
55
,
1418
(
1985
).
17.
M. S.
Hybertsen
and
S. G.
Louie
,
Phys. Rev. B
34
,
5390
(
1986
).
18.
W.
Kohn
and
L. J.
Sham
,
Phys. Rev.
140
,
A1133
(
1965
).
19.
P.
Hohenberg
and
W.
Kohn
,
Phys. Rev.
136
,
B864
(
1964
).
20.
S. V.
Faleev
,
M.
Van Schilfgaarde
, and
T.
Kotani
,
Phys. Rev. Lett.
93
,
126406
(
2004
).
21.
F.
Caruso
,
P.
Rinke
,
X.
Ren
,
M.
Scheffler
, and
A.
Rubio
,
Phys. Rev. B
86
,
081102
(
2012
).
22.
F.
Bruneval
and
M. A.
Marques
,
J. Chem. Theory Comput.
9
,
324
(
2012
).
23.
F.
Bruneval
,
J. Chem. Phys.
136
,
194107
(
2012
).
24.
M. J.
van Setten
,
F.
Caruso
,
S.
Sharifzadeh
,
X.
Ren
,
M.
Scheffler
,
F.
Liu
,
J.
Lischner
,
L.
Lin
,
J. R.
Deslippe
,
S. G.
Louie
 et al,
J. Chem. Theory Comput.
11
,
5665
(
2015
).
25.
J. W.
Knight
,
X.
Wang
,
L.
Gallandi
,
O.
Dolgounitcheva
,
X.
Ren
,
J. V.
Ortiz
,
P.
Rinke
,
T.
Körzdörfer
, and
N.
Marom
,
J. Chem. Theory Comput.
12
,
615
(
2016
).
26.
A.
Stan
,
N. E.
Dahlen
, and
R.
Van Leeuwen
,
Europhys. Lett.
76
,
298
(
2006
).
27.
A.
Stan
,
N. E.
Dahlen
, and
R.
Van Leeuwen
,
J. Chem. Phys.
130
,
114105
(
2009
).
28.
X.
Blase
and
C.
Attaccalite
,
Appl. Phys. Lett.
99
,
171909
(
2011
).
29.
J.
Deslippe
,
G.
Samsonidze
,
D. A.
Strubbe
,
M.
Jain
,
M. L.
Cohen
, and
S. G.
Louie
,
Comput. Phys. Commun.
183
,
1269
(
2012
).
30.
H.-V.
Nguyen
,
T. A.
Pham
,
D.
Rocca
, and
G.
Galli
,
Phys. Rev. B
85
,
081101
(
2012
).
31.
F.
Caruso
,
P.
Rinke
,
X.
Ren
,
A.
Rubio
, and
M.
Scheffler
,
Phys. Rev. B
88
,
075105
(
2013
).
32.
P.
Koval
,
D.
Foerster
, and
D.
Sánchez-Portal
,
Phys. Rev. B
89
,
155417
(
2014
).
33.
L.-W.
Wang
,
Phys. Rev. B
91
,
125135
(
2015
).
34.
B.
Holm
and
U.
von Barth
,
Phys. Rev. B
57
,
2108
(
1998
).
35.
F.
Bruneval
and
M.
Gatti
,
First Principles Approaches to Spectroscopic Properties of Complex Materials
(
Springer
,
2014
), pp.
99
135
.
36.
M.
van Schilfgaarde
,
T.
Kotani
, and
S.
Faleev
,
Phys. Rev. Lett.
96
,
226402
(
2006
).
37.
T.
Kotani
,
M.
van Schilfgaarde
, and
S. V.
Faleev
,
Phys. Rev. B
76
,
165106
(
2007
).
38.
F.
Kaplan
,
F.
Weigend
,
F.
Evers
, and
M.
van Setten
,
J. Chem. Theory Comput.
11
,
5152
(
2015
).
39.
F.
Caruso
,
M.
Dauth
,
M. J.
van Setten
, and
P.
Rinke
,
J. Chem. Theory Comput.
12
,
5076
(
2016
).
40.
J. E.
Northrup
,
M. S.
Hybertsen
, and
S. G.
Louie
,
Phys. Rev. Lett.
59
,
819
(
1987
).
41.
N.
Marom
,
F.
Caruso
,
X.
Ren
,
O. T.
Hofmann
,
T.
Körzdörfer
,
J. R.
Chelikowsky
,
A.
Rubio
,
M.
Scheffler
, and
P.
Rinke
,
Phys. Rev. B
86
,
245127
(
2012
).
42.
T.
Rangel
,
S. M.
Hamed
,
F.
Bruneval
, and
J. B.
Neaton
,
J. Chem. Theory Comput.
12
,
2834
(
2016
).
43.
D.
Neuhauser
,
Y.
Gao
,
C.
Arntsen
,
C.
Karshenas
,
E.
Rabani
, and
R.
Baer
,
Phys. Rev. Lett.
113
,
076402
(
2014
).
44.
V.
Vlček
,
E.
Rabani
,
D.
Neuhauser
, and
R.
Baer
,
J. Chem. Theory Comput.
13
,
4997
(
2017
).
45.
V.
Vlček
,
W.
Li
,
R.
Baer
,
E.
Rabani
, and
D.
Neuhauser
,
Phys. Rev. B
98
,
075107
(
2018
).
46.
M. R.
Filip
and
F.
Giustino
,
Phys. Rev. B
90
,
245145
(
2014
).
47.
X.
Qian
,
P.
Umari
, and
N.
Marzari
,
Phys. Rev. B
91
,
245105
(
2015
).
48.

Actually, a first order polynomial could also be incorporated, with some additional storage steps, as will be detailed in a future publication.

49.
G. J.
Martyna
and
M. E.
Tuckerman
,
J. Chem. Phys.
110
,
2810
(
1999
).
50.
J. P.
Perdew
and
Y.
Wang
,
Phys. Rev. B
45
,
13244
(
1992
).
51.
N.
Troullier
and
J. L.
Martins
,
Phys. Rev. B
43
,
1993
(
1991
).
52.
E.
Maggio
,
P.
Liu
,
M. J.
van Setten
, and
G.
Kresse
,
J. Chem. Theory Comput.
13
,
635
(
2017
).
53.

The G0W0 shifts for hexacene are −2.73 ± 0.08 and −1.63 ± 0.06 eV for the bottom valence and HOMO states.

54.
L.
Hung
,
H.
Felipe
,
J.
Souto-Casares
,
J. R.
Chelikowsky
,
S. G.
Louie
, and
S.
Öğüt
,
Phys. Rev. B
94
,
085125
(
2016
).
55.
K.
Krause
,
M. E.
Harding
, and
W.
Klopper
,
Mol. Phys.
113
,
1952
(
2015
).
56.
R. D.
Johnson
,
NIST Computational Chemistry Comparison and Benchmark Database NIST Standard Reference Database, Number 101
, 3rd ed. (
NIST
,
2016
), available at http://cccbdb.nist.gov/.
57.

Of course, just like for molecules, our self-consistency scheme is applicable to any G0W0 implementation that yields the full-frequency or full-time matrix element of the self-energy.

58.

For BN, we observe that in the zeroth iteration (G0W0) the bottom valence state and εH are shifted by −0.73 ± 0.04 and −0.21 ± 0.03 eV, respectively. However, for Si, the shifts are only 0.13 ± 0.03 and −0.16 ± 0.02 eV.

59.
Note that referenced values are obtained with PBE starting point, but
Klimeš
 et al [
Phys. Rev. B
90
,
075125
(
2014
)] and
Grüneis
 et al [
Phys. Rev. Lett.
112
,
096401
(
2014
)] show that difference between LDA and PBE starting points is negligible for the solids considered here.
[PubMed]
60.
J.
Klimeš
,
M.
Kaltak
, and
G.
Kresse
,
Phys. Rev. B
90
,
075125
(
2014
).
61.
F.
Giustino
,
Rev. Mod. Phys.
89
,
015003
(
2017
).
62.
M.
Cardona
,
Solid State Commun.
133
,
3
(
2005
).
63.
W.
Bludau
,
A.
Onton
, and
W.
Heinke
,
J. Appl. Phys.
45
,
1846
(
1974
).
64.
D.
Bimberg
,
M.
Altarelli
, and
N.
Lipari
,
Solid State Commun.
40
,
437
(
1981
).
65.
M.
Lorenz
,
R.
Chicotka
,
G.
Pettit
, and
P.
Dean
,
Solid State Commun.
8
,
693
(
1970
).
66.
C. J.
Wort
and
R. S.
Balmer
,
Mater. Today
11
,
22
(
2008
).
67.
M. E.
Levinshtein
,
S. L.
Rumyantsev
, and
M. S.
Shur
,
Properties of Advanced Semiconductor Materials: GaN, AIN, InN, BN, SiC, SiGe
(
John Wiley & Sons
,
2001
).
68.
J.
Towns
,
T.
Cockerill
,
M.
Dahan
,
I.
Foster
,
K.
Gaither
,
A.
Grimshaw
,
V.
Hazlewood
,
S.
Lathrop
,
D.
Lifka
,
G. D.
Peterson
 et al,
Comput. Sci. Eng.
16
,
62
(
2014
).