Skip to main content

REVIEW article

Front. Chem., 09 July 2019
Sec. Theoretical and Computational Chemistry
Volume 7 - 2019 |

The GW Compendium: A Practical Guide to Theoretical Photoemission Spectroscopy

  • Department of Applied Physics, Aalto University, School of Science, Espoo, Finland

The GW approximation in electronic structure theory has become a widespread tool for predicting electronic excitations in chemical compounds and materials. In the realm of theoretical spectroscopy, the GW method provides access to charged excitations as measured in direct or inverse photoemission spectroscopy. The number of GW calculations in the past two decades has exploded with increased computing power and modern codes. The success of GW can be attributed to many factors: favorable scaling with respect to system size, a formal interpretation for charged excitation energies, the importance of dynamical screening in real systems, and its practical combination with other theories. In this review, we provide an overview of these formal and practical considerations. We expand, in detail, on the choices presented to the scientist performing GW calculations for the first time. We also give an introduction to the many-body theory behind GW, a review of modern applications like molecules and surfaces, and a perspective on methods which go beyond conventional GW calculations. This review addresses chemists, physicists and material scientists with an interest in theoretical spectroscopy. It is intended for newcomers to GW calculations but can also serve as an alternative perspective for experts and an up-to-date source of computational techniques.

1. Introduction

Electronic structure theory derives from the fundamental laws of quantum mechanics and describes the behavior of electrons—the glue that shapes all matter. To understand the properties of matter and the behavior of molecules, the quantum mechanical laws must be solved numerically because a pen and paper solution is not possible. In this context, Hedin's GW method (Hedin, 1965) has become the de facto standard for electronic structure properties as measured by direct and inverse photoemission experiments, such as quasiparticle band structures and molecular excitations.

Electronic structure theory covers the quantum mechanical spectrum of computational materials science and quantum chemistry. The fundamental aim of computational science is to derive understanding entirely from the basic laws of physics, i.e., quantum mechanical first principles, and increasingly also to make predictions of new properties or new materials and new molecules for specific tasks. The exponential increase in available computer power and new methodological developments are two major factors in the growing impact of this field for practical applications to real systems. As a result of these advances, computational science is establishing itself as a viable complement to the purely experimental and theoretical sciences.

Hedin published the GW method in 1965, in the same time period as the foundational density-functional theory (DFT) papers (Hohenberg and Kohn, 1964; Kohn and Sham, 1965). While DFT has shaped the realm of first principles computational science like no other method today, GW's fame took a little longer to develop1. DFT's success has been facilitated by the computational efficiency of the local-density (Kohn and Sham, 1965) or generalized gradient approximations (Becke, 1988; Lee et al., 1988; Perdew et al., 1996a) (LDA and GGA) of the exchange-correlation functional that make DFT applicable to polyatomic systems containing up to several thousand atoms. GW, however, only saw its first applications to realistic materials 20 years after its inception (Hybertsen and Louie, 1985; Godby et al., 1986), due to its much higher computational expense. Soon after it was realized that GW can overcome some of the most notorious deficiencies of common density functionals such as the self-interaction error, the absence of long-range polarization effects and the Kohn-Sham band-gap problem.

The GW approach is now an integral part of electronic structure theory and readily available in major electronic structure codes. It is taught at summer schools along side DFT and other electronic structure methods. This review is intended as a tutorial that complements showcases of GW's achievements with a practical guide through the theory, its implementation and actual use. For GW novices, the review offers a gentle introduction to the GW concept and its application areas. Regular GW users can consult this review as a handbook in their day-to-day use of the GW method. For seasoned GW users and GW experts it might serve as a reference for key applications and some of the subtler points of the GW framework.

In this review we take a more practical approach toward the GW method. We will recap the basic theory starting from theoretical spectroscopic view point as a probe of the electronic structure. Aiming at GW practitioners, we will illustrate how the GW approach emerges from the theoretical spectroscopy framework as a practical scheme for electronic structure calculations. A more in-depth discussion of the theoretical foundations of many-body theory can be found in textbooks (e.g., Fetter and Walecka, 1971; Szabo and Ostlund, 1989; Gross et al., 1991; Bechstedt, 2014; Martin et al., 2016), while the GW theory itself is covered in excellent early reviews (Hedin and Lundqvist, 1970; Aryasetiawan and Gunnarsson, 1998; Hedin, 1999). GW began to flourish at the beginning of the 21st century and two reviews succinctly summarized the state of the field until then (Aulbur et al., 2000; Onida et al., 2002). Our review bridges the ensuing gap of almost 20 years, after which only several more specialized reviews addressed different aspects and applications of GW calculations (Rinke et al., 2008a; Giantomassi et al., 2011; Ping et al., 2013; Bruneval and Gatti, 2014; Faber et al., 2014; Marom, 2017; Gerosa et al., 2018a; Kang et al., 2019), and complements a recent review (Reining, 2017).

A considerable part of our review is devoted to the presentation of different GW implementations. We will discuss the practical considerations that GW users have to make when they decide on a particular GW implementation or code for their work. Moreover, we will guide the reader through computational decisions that might affect the accuracy of their GW calculations and illustrate them with concrete examples from the GW literature. An important aspect in this regard is the issue of self-consistency in GW, which we cover in detail.

Although the GW method might be best known for its success in predicting the band gaps of solids, we will present its diversity and discuss a range of different quantities that can be computed with the GW method. Since no method is perfect, we will conclude with a critical outlook on the challenges faced by the GW method and discuss ways to go beyond GW.

We conclude this introduction with a quote from H. J. Monkhorst, who wrote in a laudation in 2005: It is therefore my conviction that, rather sooner then later, we will see a resurgence of the precise many-body approach to solid-state theory as we envisioned. Almost assuredly the GW method will be the tool of choice (Monkhorst, 2005). In 2019, we can say that Monkhorst was right.

2. Theoretical Spectroscopy

2.1. Direct and Inverse Photoemission Spectroscopy

Spectroscopic measurements are an important component in the characterization of materials. Any spectroscopic technique perturbs the system under investigation and promotes it into an excited state. Experimentally, the challenge then lies in the correct interpretation of the system's response. From a theoretical point of view, however, the challenge is to find (or develop) a suitable, accurate and, most of all, computationally tractable approach to describe the response of the system. Experimental and theoretical spectroscopy are often complementary and, when combined, they are a powerful approach.

Many-body perturbation theory (MBPT) provides a rigorous and systematic quantum mechanical framework to describe the spectral properties of a system that connects central quantities like the Green's function, the self-energy, and the dielectric function with each other. The poles of the single-particle Green's function, the central object in MBPT, correspond to the electron addition and removal energies probed in direct and inverse photoemission, which is explained in detail at the end of section 2.1. In contrast, information about neutral excitations probed in optical or energy loss spectra can be extracted from the dielectric function. In this review, we will not address optical properties or other neutral excitations and instead focus on the single-particle Green's function and its connection to direct and inverse photoemission spectroscopy.

In photoelectron spectroscopy (PES) (Plummer and Eberhardt, 1982; Himpsel, 1983; Kevan, 1992), electrons are ejected from a sample upon irradiation with visible or ultraviolet light (UPS) or with X-rays (XPS), as sketched in Figure 1A. The energy of the bound state ϵs can be reconstructed from the photon energy , the work function Φ and the kinetic energy Ekin of the photoelectrons that reach the detector2

IPs=ϵs=hνEkinΦ, for ϵs<EF.    (1)

The ionization potential IPs is defined as the energy that is required to remove an electron from the bound initial state s of the neutral sample, where the energy of the state is below the Fermi level (EF). It is always a positive number and related to ϵs as shown in Equation (1).


Figure 1. Schematic of the photoemission (PES) and inverse photoemission (IPES) process. In PES (A) an electron is excited by an incoming photon from a previously occupied valence state (lower shaded region) into the continuum (gray shaded region, starting above the vacuum level Evac). In IPES (B) an injected electron with kinetic energy Ekin undergoes a radiative transition into an unoccupied state (upper shaded region) thus emitting a photon in the process.

By inverting the photoemission process, as schematically shown in Figure 1B, the unoccupied states can be probed. This technique is commonly referred to as inverse photoemission spectroscopy (IPES) or bremsstrahlung isochromat spectroscopy (BIS) (Dose, 1985; Smith, 1988; Fuggle and Inglesfield, 1992). In IPES, an incident electron with energy Ekin is scattered in the sample emitting bremsstrahlung. Eventually it will undergo a radiative transition into a lower-lying unoccupied state, emitting a photon that carries the transition energy . The energy of the final, unoccupied state ϵs can be deduced from the measured photon energy according to

EAs=ϵs=Ekinhν+Φ, for ϵsEF.    (2)

EA denotes the electron affinity, which we define as the energy needed to detach an electron from a negatively charged species and which is thus the negative of ϵs. EA can be a positive or negative number. It is negative when the additional electron is in an unbound state, and positive when the electron is bound.

The experimental observable in photoemission spectroscopy is the photocurrent, which is the probability of emitting an electron with the kinetic energy Ekin within a certain time interval. It is related to the intrinsic spectral function A(r, r′, ω) of the electronic system, given by the imaginary part of the single-particle Green's function3: (Almbladh and Hedin, 1983; Onida et al., 2002)

A(r, r,ω)=1πIm G(r, r,ω)sgn(EFω),    (3)

where ω denotes an energy (frequency). The single-particle Green's function, G(r, r′, ω), is the probability amplitude that a particle created or destroyed at r is correlated with the adjoint process at r′ − it will be discussed in detail later. The actual dependence of the photocurrent on the spectral function is quite complicated because the coupling to the exciting light and electron loss processes in the sample, as well as surface effects, have to be taken into account. To our knowledge, no comprehensive theory yet exists for this relation and we therefore proceed with the discussion of the spectral function and will return to the photocurrent later.

The energies ϵs in Equations (1) and (2) are the removal and addition energies of the photoelectron, respectively, and we refer to the transition amplitudes from the N to the N ± 1-body states as ψs(r) (see also section B.1 in Appendix B):

ϵs=E(N)E(N1,s)ψs(r)=N1,s|ψ^(r)|N} forϵs<EF    (4)
ϵs=E(N+1,s)E(N)ψs(r)=N|ψ^(r)|N+1,s} forϵsEF    (5)

The states |N, s〉 are many-body eigenstates (wave functions in real space) of the N-electron Schrödinger equation Ĥ|N, s〉 = E(N, s)|N, s〉, Ĥ is the many-body Hamiltonian and E(N, s) = 〈N, s|Ĥ|N, s〉 is the corresponding total energy. The field operator ψ^(r) annihilates an electron at point r from the many-body states |N〉 or |N + 1〉. The representation given in Equations (4) and (5) is particularly insightful because it allows a direct interpretation of ϵs as the photoexcitation energy from the N-particle ground state with total energy E(N) into an excited state s of the (N − 1)-particle system with total energy E(N − 1, s) upon removal of an electron in the photoemission process. Similarly, the addition energy that is released in the radiative transition in inverse photoemission is given by the total energy difference of the excited (N+1)-particle system and the ground state.

To build a practical scheme for calculating the energies in Equations (4) and (5) we introduce the definition of the single-particle Green's function4

G(r,σ,t,r,σt)=iN|T^{ψ^(r,σ,t)ψ^(r,σ,t)}|N    (6)

where T^ is the time ordering operator for the times t and t′ and σ the spin. T^ arranges the field operators so that the earlier time is to the right and acts on the ground state |N〉 first. G allows for both time orderings: t > t′ or t′ > t. This definition of the Green's function is particularly insightful because it illustrates the process of adding and removing electrons from the system, as done in photoemession spectroscopy. Assuming the time-ordering is as shown in Equation (6), ψ^(r,σ,t) will create an electron with spin σ′ at time t′ in point r′. This electron will then propagate through the system, until it is annihilated by ψ^(r,σ,t) at a later time t in position r. The Green's function is therefore also often called a propagator. We will return to this propagator picture of G in later sections of this review.

To make contact with Equations (4) and (5), we need to Fourier transform the Green's function from the time to the energy axis. For a time-independent Hamiltonian this then produces the spectral or Lehman representation of G (Fetter and Walecka, 1971; Gross et al., 1991)

G(r,r,ω)=limη0+sψs(r)ψs*(r)×  ×[Θ(ϵsEF)ω(ϵsiη)+Θ(EFϵs)ω(ϵs+iη)]    (7)

where we have assumed a spin paired system and summed over the spin quantum number shown in Equation (6). The two terms in brackets are for the two time orderings in G. Θ is the Heaviside step function, which is zero for negative arguments and one for positive arguments5. It kills any processes that do not obey the correct time ordering, as determined by the created/annihilated particle's energy relative to EF. This representation illustrates that the many-body excitations of the system that are associated with the removal or addition of an electron are given by the poles of the single-particle Green's function. The diagonal spectral function

A(r,r,ω)=1πIm G(r,r,ω)sgn(EFω)    (8)
=sψs(r)ψs*(r)δ(ωϵs)    (9)

then assumes the intuitive form of a (many-body) density of states.

2.2. The Quasiparticle Concept

In periodic solids, the crystal has special crystallographic directions so that spectra are direction dependent, with the direction indexed by a wave vector k. By varying the direction of the incident beam relative to the crystallographic axes (ai), one can map the k dependent PES spectra, as shown in Figure 2A. This technique is called angle resolved photoemission spectroscopy (ARPES). Figure 2C shows data from a typical ARPES measurement and Figure 2B a schematic of the spectral function at a single k-point. The spectra usually exhibit distinct peaks that are attributed to particle-like states but have a finite width. There can also be additional, broader peaks away from the main peak called satellites. However, the spectral function in Equation (9) contains only Dirac-delta functions which appear as infinitely sharp peaks. The broadening of the spectral function comes from the sum of many delta functions close in energy, which merge to form a peak of finite width. If the contributing delta functions are closely packed around one energy, the peak is attributed to a quasiparticle (Landau et al., 1980).


Figure 2. (A) Schematic representation of an ARPES experiment. By varying the angles θ and ϕ with respect to the crystallographic axes (ai), the measured spectrum is direction, or k, dependent. In practice, the detector angle is usually varied with respect to a fixed beam. (B) A typical spectral function features a sharp peak attributed to the quasiparticle, an incoherent background, and satellites away from the single particle peak. (C) ARPES data of the upper valance bands of ZnO (Kobayashi et al., 2009). The corresponding G0W0 band structure of ZnO is shown in Figure 27.

To further motivate the association of quasiparticles with particle-like excitations it is insightful to consider non-interacting electrons. In that case, the spectral function consists of a series of delta peaks

Ass(ω)=ψs|A(ω)|ψs=δssδ(ωϵs),    (10)

each of which corresponds to the excitation of a non-interacting particle, see Appendix A for the integral notation used in Equation (10). The many-body states |N〉 and |N ± 1〉 become single Slater determinants so that the exact excited states are characterized by a single creation or annihilation operator acting on the ground state. The excitation energies ϵs and the wave functions ψs(r) are thus the eigenvalues and eigenfunctions of the single-particle Hamiltonian.

When the electron-electron (or electron-ion) interaction is turned on, the exact eigenstates |N, s〉 are no longer single Slater determinants. As a consequence, the matrix elements of the spectral function Ass(ω) will contain contributions from many non-vanishing transition amplitudes. If these contributions merge into a clearly identifiable peak that appears to be derived from a single delta-peak broadened by the electron-electron interaction, this structure can be interpreted as a single-particle like excitation—the quasiparticle. The broadening of the quasiparticle peak in the spectral function is associated with the lifetime τ of the excitation due to electron-electron scattering, whereas the area underneath the peak is interpreted as the renormalization Z of the quasiparticle. This renormalization factor quantifies the reduction in spectral weight due to the electron-electron interaction compared to an independent electron, though the total spectral weight is conserved. We can combine these various arguments and say that the quasiparticle peak for state s will exhibit the following shape:

Ass(ω)1π|Zsω(ϵs+iΓ)|.    (11)

In contrast to the exact energies of the many-body states, which are poles of the Green's function on the real axis, the quasiparticle poles reside in the complex plane and are no longer eigenvalues of the N-body Hamiltonian. The real part of this complex energy is associated with the energy of the quasiparticle excitation and the imaginary part with its inverse lifetime Γ = 2/τ.

To develop a more intuitive understanding of quasiparticles, it is insightful to adopt a real-space picture. The quasiparticle concept can be explained by analogy with a crowd of people, as shown in Figure 3. Picture a group of people, such as at a concert or festival, all crowded into the same area. Not wanting to get too close to each other to preserve their own space, people in the crowd interact with each other. If one person gets too close to another, their mutual repulsion eventually takes over and separates them again. The exact description of the crowd requires the location of each individual person at all times. This is a very difficult task because of the constant interactions, or repulsions, between individual people. This collection of people and their occasional fluctuations are grouped together and labeled the ground state.


Figure 3. Top: Depiction of the quasiparticle concept. (A) A crowd of people is analogous to the electronic ground state. A new person (that represents an additional electron) enters the crowd in (B). The new person begins to interact with other people who, in turn, interact back with the new person in (C) and form a polarization cloud. An effective, or renormalized, object, the quasi-person, moves through the crowd in (D). Even though it is an interacting system, the many-person state in (D) can still be connected to, or identified by, the single person added to the crowd. This connection allows us to identify the quasi-person. Bottom: Schematic representation of photoemission spectroscopy.

A new person arrives and pushes their way into the crowd. We can think of this new person as the electron in inverse photoemission that is injected into the system. The new person enters in a specific direction with a certain energy. As they enter the group, they repeatedly interact with other people as they continue their trajectory, as shown in Figure 3C. These repeated interactions repel people in their immediate vicinity and form a small halo of free space around the incoming person. People seem to move out of their way on their journey, forming a polarization cloud created by the absence of other people around them. The intruder's motion and their polarization cloud can be taken together to form a new composite object, a quasi-person, which appears as a slowed-down version of the newcomer. From far away, one does not need to describe the precise motion of all N + 1 people in the group, but only the motion of this quasi-person propagating through the crowd.

By analogy, a quasiparticle can then be considered a combination of an additional electron or hole in the system that interacts with its surrounding polarization cloud. The situation corresponding to photoemission spectroscopy is depicted in the bottom panel of Figure 3. As time increases, the bare hole left by removal of the interaction is screened. The quasiparticle therefore embodies an electron state with the perturbation of its own surrounding. The feedback via interactions of the particle with surrounding electrons is termed the self-energy. Over time, the propagating quasiparticle can decay into many different elementary excitations, giving it a finite lifetime. Essential quasiparticle properties are dispersion, lifetime, weight, and satellite spectrum. The latter arises from the collective excitations in the medium.

2.3. Comparison to Experimental Spectra

We have now identified quasiparticles as one possible source for peaks in experimental photoemission spectra. Before we introduce the GW approximation as a tractable computational approach for calculating quasiparticle energies, we will first address the photocurrent, which is the quantity measured in direct photoemission experiments. Then we will briefly discuss the reconstruction of the band structure information, as well as other sources of peaks in spectrum.

Establishing rigorous links between the spectral function and the photocurrent is still a challenge for theory (Hedin, 1999; Lee et al., 1999; Minár et al., 2011, 2013). The photocurrent Jk() is the probability per unit time of emitting a photoelectron with momentum k and energy Ekin, k due to an incident photon with the energy . The spectral function defined in Equation (3) describes the removal of an electron from the sample, but does not include intermediate steps on the way to the detector where the electron loses energy. Therefore, it does not correspond to Jk(). However, the spectral function can be related to the photocurrent by using the sudden approximation (Hedin, 1999; Hedin and Lee, 2002) assuming that the ejected photoelectron is immediately decoupled from the sample. Jk() is then given by (Hedin, 1999)

Jk(hν)= s,sΔksAss(Ekin,khν)Δsk,    (12)

where Ass(ω)=ψs|A(ω)|ψs are matrix elements of the spectral function defined in Equation (3). Δk, s are matrix elements of the dipole operator which describe the coupling to photons. The dipole matrix elements capture the promotion of the electron to a highly excited state (often assumed to be a plane wave), i.e., they describe the transition between the initial and the final electron state. In this final state, the electron travels to the detector. On the way, it crosses the surface of the sample, which adds a further perturbation to its path and its energy. The transition matrix elements affect the amplitudes of the spectrum and add selection rules that give rise to the suppression or enhancement of certain peaks. In practice, one compares only matrix elements of the spectral function to the experiment disregarding the effects of the dipole matrix. Furthermore, it is often assumed that only the diagonal elements of the spectral function are dominating.

For the reconstruction of the band structure, e.g., with ARPES, the comparison between theory and experiment is hampered from the experimental side. In ARPES studies of crystalline materials, the emitted photons or electrons inevitably have to pass the surface of the crystal to reach the detector. Therefore, information about their transverse momentum k is lost. This is because the crystal's translational symmetry is broken at the surface and only the in-plane momentum k is conserved. To reconstruct the three-dimensional band structure of the solid from experimental data, assumptions are often made about the dispersion of the final states (Plummer and Eberhardt, 1982; Himpsel, 1983; Hora and Scheffler, 1984; Dose, 1985; Smith, 1988). Ab initio calculations as described in this article can aid in the assignment of the measured peaks. Either way, some layer of interpretation between theoretical and experimental band structures is required.

Apart from quasiparticle excitations, a typical photoemission experiment provides a rich variety of additional information. In core-electron emission for instance, inelastic losses or multi-electron excitations such as shake-ups and shake-offs lead to satellites in the spectrum. Satellites can also appear in the valence region. The outgoing photoelectron or the hole left behind can, for example, excite other quasiparticles like plasmons, phonons or magnons. This gives rise to additional peaks, the so-called plasmon or magnon satellites or phonon side bands, that are typically separated from the quasiparticle peak by multiples of the plasmon, magnon or phonon energy. The broad peak in Figure 4, which shows integrated spectra and therefore has no k dependence, near −40 eV is an example of a satellite. Satellites are collective effects that are not described within the quasiparticle picture.


Figure 4. X-ray photoemission spectrum with 800 eV incident energy compared to two calculated spectra. The red line shows the evGW0@LDA spectrum (see section 5), whereas the blue spectrum contains additional vertex corrections in form of a cumulant expansion (see section 11). The evGW0@LDA+C* spectrum contains the addition of the Shirley background (shown by the black dashed line) and loss effects of the outgoing photoelectron. Data retrieved from Guzzo et al. (2011), where the GW results are labeled as G0W0. However, self-consistency in the eigenvalues was in fact applied, which is in our notation evGW0 (Private Communication).

3. Hedin's GW Equations

Having introduced the general Green's function framework and quasiparticle concept, we are prepared to consider the concrete formalism for GW. GW is an approximation to an exact set of coupled integro-differential equations called Hedin's equations (Hedin, 1965), the full derivation of which can be found in Appendix B.

We build up Hedin's equations from a perturbation theory perspective. We can conveniently represent the perturbation expansion for G that we introduced in the time domain in Equation (6) and in the energy domain in Equation (7) with the Feynman diagram technique. Feynman diagrams are a pictorial way of representing many-body and Green's function theory. We cover the necessary basics in this section and refer the interested reader to an excellent book on Feynman diagrams (Mattuck, 1992).

The perturbation expansion begins with the noninteracting Green's function, denoted G0. G0 is the probability amplitude for a noninteracting particle to propagate from one spacetime point to another. In the diagrammatic technique, G0 is represented by a solid line with an arrow. The ends of G0 indicate spacetime points. The generic notation 1 = (r1, t1, σ1) refers to the spatial coordinate r1, time t1, and spin variable σ1. G0, shown in Figure 5, is one of the basic building blocks for the perturbation expansion.


Figure 5. The most basic pieces of diagrammatic perturbation theory are G0 and v. From these, all other quantities can be built. The interaction v(1, 2) is instantaneous. Therefore, the dashed line is perpendicular to the time axis. The arrows in G0 and G point in only one direction, but both time orderings are included.

G is the probability amplitude for the interacting system that a particle creation at 2 is correlated with a particle annihilation at 1. The rules of quantum mechanics dictate that we must sum over all possible paths for the particle to move from 2 to 1 − this generates the exact G, which is represented by a bold, or double, line with an arrow in the diagram language. Every possible process between 1 and 2 contributes a different amplitude to the total G. The different processes which connect 1 to 2 depend on interactions with other particles in the system at the times between t1 and t2. Without these interactions, the problem is already solved with G0.

Recall from the definition of G in Equation (6) that G contains two time orderings. The second time ordering implies that the annihilation process may come before the creation. Remember that the field operators in Equation (6) act on the interacting ground state. In the ground state, there is some charge for the annihilation operator ψ^(r) to “act” on, even without any preceding creation process, so that the reverse time ordering in G makes sense. Feynman diagrams do not explicitly show both time orderings in G, but it is important to remember that G and G0 lines implicitly contain both time orderings.

The times between t1 and t2 are called internal times. We can add up all the processes contributing to G in a certain order depending on the number of times the particle interacts with other particles in the system. These interactions occur only at internal times, and the number of internal interactions is the order of the diagram. To efficiently represent all these internal interactions, we use a dashed line to represent the interaction in the diagram. At a given order n, we construct all possible processes, or diagrams, which connect n interaction lines with G0 lines at 1 and 2. We connect all of the dashed lines appearing at internal times with additional G0 lines. There is a very specific set of rules for how these arrangements can be done. Wick's theorem defines how to contract these pieces (Fetter and Walecka, 1971). A simple principle is enough to demonstrate the idea, however. Because the Coulomb interaction is a two-body operator, each dashed line must have two G0 lines at each end. To compute the exact G, the expansion must be taken to infinite order, n → ∞, adding up all possible processes along the way. The process of building up all diagrams in the perturbation expansion is shown in Figure 6. A few example diagrams, as well as a couple of forbidden diagrams, are shown in Figure 7.


Figure 6. The exact G contains amplitudes from all possible paths between 1 and 2. Amplitudes from all of these paths are represented by the rectangle placed between the field operators, the action of which is represented by * symbols. These terms can be calculated order-by-order with perturbation theory. At a given order n, we must connect n interaction lines at internal times in all possible − and allowed − ways. Concrete examples of diagrams are in Figure 7. All terms of the topology which can be inserted between two G0 lines form the reducible self-energy.


Figure 7. At first order, n = 1, there are only two possible self-energy diagrams. These are the diagrams of the Hartree-Fock approximation, the direct electrostatic interaction (left) and exchange (right). Two possible n = 2 diagrams are also shown (there are others). The bottom diagrams are forbidden because they do not have two G0 lines at each end of the interaction lines. When drawing the diagrams, a certain degree of flexibility is allowed and they must be interpreted carefully. For example, the curved interaction lines above must still be treated as instantaneous in a calculation.

To go further with our analysis, we must dissect the internal structure of the diagrams and separate it into pieces. By considering the possible topologies of internal parts allowed by the contraction rules, we can group the parts into different categories. Here, the “topology” of the piece is determined by the number of G0 lines and interactions at two different times, without considering the internal structure between the two chosen times. Those parts which have two G0 lines sticking out are called a self-energy diagram. The full self-energy (Σ) can be inserted between two G0 lines to form G (one must also include the separate G0 term). Topologies which connect two G0 legs to an interaction line are labeled a vertex. Summing over all pieces with this topology creates the full vertex (Γ), which depends on three spacetime points. Finally, the diagram parts which end in two dashed lines sum up to the effective, or screened, interaction (W).

Conceptually, the vertex is the most difficult to understand. Figure 8 demonstrates the effect of the vertex in a specific example. The diagram shown in Figure 8A is meant to contain the exact vertex, Γ. Γ has three corners and can be inserted where two G0 lines meet an interaction line. By simply letting these three pieces meet without any internal structure, we replace Γ with a single spacetime point, as shown in Figure 8B. Alternatively, we could allow the vertex to include the curved interaction line shown in Figure 8C. In that case, the vertex has internal structure.


Figure 8. The exact vertex Γ, shown in (A), can be replaced with approximations to simplify the calculation. The approximation in (B) is referred to as a “single spacetime point” because the vertex has no internal structure. In contrast, the vertex in (C) has internal structure. The diagram shown here is only an example to demonstrate the role of Γ and does not correspond to the exact self-energy or the GW self-energy.

Hedin's equations can be interpreted as the self-consistent formulation of these topologically distinct building blocks. While Hedin followed a formal and systematic derivation, a heuristic motivation is to group all diagrams of a certain topology together and replace them with a single dressed, or renormalized, object with the same topology. A critical aspect of this replacement is their energy dependence. By replacing many diagrams of perturbation theory with a single object of the same shape, we reduce the number of objects to be computed. However, the information apparently missing due to the reduction in objects is encoded in the energy dependence of the dressed quantity. The final result, Hedin's equations (Appendix B), are a compact and self-contained set of five integro-differential equations. Despite the reduction in the number of objects to be treated compared to the perturbation expansion, the functional differential equations coupling these pieces are extremely difficult to solve exactly.

Hedin recognized this difficulty and suggested the GW approximation. As mentioned above, the vertex is the building block which has a single interaction connected to two G0 lines. Unlike the other building blocks, the vertex depends on three spacetime points instead of two, making it the most difficult to compute. To simplify the theory, Hedin suggested replacing Γ with a single spacetime point. In Hedin's equations, the exact self-energy is Σ = iGWΓ. With the replacement Γ(1, 2, 3) = δ(1, 2)δ(1, 3), Hedin's approximation gives Σ = iGW, hence the name of the GW approximation. In this approximation, Hedin's equations are

G(1,2)=G0(1,2)  + G0(1,3)Σ(3,4)G(4,2)d(3,4)    (13)
Γ(1,2,3)=δ(1,2)δ(1,3)    (14)
χ0(1,2)=iG(1,2)G(2,1)    (15)
W(1,2)=v(1,2)+ v(1,3)χ0(3,4)W(4,2)d(3,4)    (16)
Σ(1,2)=iG(1,2)W(1+,2)    (17)

where the Hartree potential is included in the solution for G06.

These are the GW equations, which are translated into the diagram language in Figure 9. There is one final important point regarding the reducibility of the quantities in Hedin's equations. Σ and χ0 in Hedin's equations are irreducible, which means that they cannot be broken into smaller pieces with the same topology. To generate the full, or reducible, quantity from its irreducible part, the irreducible component is iterated in a series similar to the perturbation expansion for G. Series of this type are commonly called Dyson series, and Dyson's equation refers to the equation for G shown in Equation (13) (W in Equation (16) also obeys a Dyson series). Dyson's equation is of great importance in many-body physics, and we return to it in later sections in the context of self-consistent GW. It is common in the literature to use the same symbol for both reducible and irreducible components with the same topology, especially when discussing the self-energy. Almost always, the symbol refers to the quantity as it appears in Hedin's equations. When discussing or calculating the self-energy, this implies that one is interested in the irreducible self-energy.


Figure 9. Diagrammatic representation of Equations (13–17). The GW approximation reduces the self-energy to a product of G with W. The first equation (Dyson's equation) has a G line on the left- and right-hand sides. This equation can be iterated, inserting G0 + G0ΣG in place of each G on the RHS, forming the Dyson series. The same iterative procedure for W forms its own Dyson series.

In Hedin's equations, the screened Coulomb interaction W plays a central role. Screening is based on the simple idea that charges in the system rearrange themselves to minimize their interaction. In polarizable materials, screening is significant, and the effective interaction is noticeably weaker than the bare one. W is also dependent on frequency or a time difference. The frequency dependence of W is critical to both the physics and the numerical implementation of GW. Even though the bare interaction is instantaneous, W is time difference dependent because it is built from repeated bare interactions at different times. This series of bare interactions to form W can be built by iterating the fourth line in Figure 9. The underlying G0 lines which connect the bare interactions in the W expansion are themselves dependent on a time difference, so that if we vary the initial or final times the entire expansion changes magnitude. This series of repeated bare interactions connected by G0 is the microscopic mechanism for the quasiparticle screening concept developed in section 2.2. The frequency dependence of W is what allows the system to relax and screen the quasiparticle. The GW self-energy is similar to the bare exchange in Hartree-Fock theory, which can be written as the product of G with v. Given the similarity between the GW self-energy and bare exchange, GW can be thought of as a dynamically screened version of Hartree-Fock.

The GW equations should still be solved self-consistently since all four quantities are coupled to each other. As with other nonlinear equations, including the equations of mean-field theories like Kohn-Sham DFT or Hartree-Fock theory, the GW equations can be solved by iteration. In principle, the prescription is clear. Start from a given G0 and iterate Equations (13–17) to self-consistency (scGW). However, remarkably few fully self-consistent solutions of the GW equations have been performed in the last 50 years. The first calculations for the homogeneous electron gas (HEG) were reported at the turn of the previous century (Holm and von Barth, 1998; Holm, 1999; García-González and Godby, 2001) and reported worse agreement with experiment on quasiparticle band widths and satellite structure compared to non-self-consistent calculations. They were quickly followed by calculations for real solids, like silicon and sodium (Schöne and Eguiluz, 1998; Ku and Eguiluz, 2002). Self-consistency was then dropped for several years because of its high computational expense and the success of non-self-consistent approximations. More recent scGW studies for atoms (Delaney et al., 2004; Stan et al., 2006, 2009), molecules (Rostgaard et al., 2010; Caruso et al., 2012a, 2013a,b; Marom et al., 2012), conventional solids (Kutepov et al., 2009; Grumet et al., 2018) and actinides (Kutepov et al., 2012) have been reported. In practice, non-self-consistent calculations are much more common, and even self-consistent GW calculations come in different types. scGW is discussed in more detail in section 5.

4. The G0W0 Approach: Concept and Implementation

4.1. The G0W0 Equations

The lowest rung in the hierarchy of GW approximations is the widely used G0W0 approach. Starting from a mean-field Green's function, G0W0 calculations correspond to the first iteration of Hedin's equations. We denote the self-energy of such single-shot perturbation calculations Σ0. Since we always refer to the single-shot self-energy in section 4, we drop the label. Furthermore, we define the single-particle Hamiltonians

ĥ0=122+vext    (18)
ĥ=122+vext+vH=ĥ0+vH    (19)
ĥMF=122+vext+vH+vσMF=ĥ+vσMF    (20)

where vext is the external potential, vH is the Hartree potential, and vσMF is the mean-field (MF) exchange-correlation potential. The spin channel is denoted by σ. Possible mean-field Hamiltonians are the Kohn-Sham (KS) or Hartree-Fock (HF) Hamiltonians.

From Dyson's equation for G, one can derive an effective single-particle eigenvalue problem referred to as the quasiparticle (QP) equation. The solutions of the QP equation are then given by

h^MF(r)ψsσ(r)drvσMF(r,r)ψsσ(r)+  drΣσ(r,r,ϵsσ)ψsσ(r)=ϵsσψsσ(r) .    (21)

The self-energy is calculated with a G0 chosen to match the initial mean-field calculation based on ĥMF. The solution of Equation (21) provides the QP energies {ϵ} and wave functions {ψ}.

Most commonly, the QP wave functions are approximated with the eigenfunctions {ϕsσ0} of the mean-field Hamiltonian. Projecting each side of Equation (21) onto ϕsσ0 yields a set of QP equations

ϵsσ=ϵsσ0+ϕsσ0|Σσ(ϵsσ)vσMF|ϕsσ0,    (22)

where {ϵsσ0} are the eigenvalues of ĥMF. Solving Equation (22), the QP energy ϵ is obtained by correcting the mean-field eigenvalue ϵsσ0.

To solve Equation (22), we have to calculate the G0W0 self-energy Σσ,

Σσ(r,r,ω)=i2πdωeiωηG0σ(r,r,ω+ω)W0(r,r,ω)    (23)

where ω is the frequency at which the self-energy is computed. Equation (23) is the frequency space version of Equation (17) for the GW self-energy. The Green's function G0σ stems from the aforementioned mean-field Hamiltonian and is given by

G0σ(r,r,ω)=mϕmσ0(r)ϕmσ0*(r)ωϵmσ0iηsgn(EFϵmσ0).    (24)

W0 in Equation (23) is the screened Coulomb interaction in the random-phase approximation (RPA)

W0(r,r,ω)= drε1(r,r,ω)v(r,r),    (25)

with the bare Coulomb interaction v(r, r′) = 1/|rr′| and the dynamical dielectric function ε. The latter is given by

ε(r,r,ω)=δ(r,r)drv(r,r)χ0(r,r,ω).    (26)

In G0W0, the irreducible polarizability χ0,

χ0(r,r,ω)=      i2πσdωG0σ(r,r,ω+ω)G0σ(r,r,ω),    (27)

simplifies to the Adler-Wiser expression (Adler, 1962; Wiser, 1963)

χ0(r,r,ω)=           σioccavirt{ϕiσ0*(r)ϕaσ0(r)ϕaσ0*(r)ϕiσ0(r)ω(ϵaσ0ϵiσ0)+iη                                ϕiσ0(r)ϕaσ0*(r)ϕaσ0(r)ϕiσ0*(r)ω+(ϵaσ0ϵiσ0)iη},    (28)

where the index i denotes an occupied and a an unoccupied (also called virtual) single-particle orbital.

For numerical convenience as well as insight into the underlying physics, the G0W0 self-energy is often split into a correlation part Σσc,

Σσc(r,r,ω)=i2πdωG0σ(r,r,ω+ω)W0c(r,r,ω),    (29)

where W0c is defined as

W0c(r,r,ω)=W0(r,r,ω)v(r,r),    (30)

and an exchange part

Σσx(r,r)=i2πdωeiωηG0σ(r,r,ω+ω)v(r,r)    (31)
=ioccϕiσ0(r)ϕiσ0*(r)v(r,r).    (32)

Note that the exponential factor in Equation (31) is necessary to close the integration contour, whereas W0c(r,r,ω) falls of quickly with increasing frequency and we can take the zero limit of η before integrating. For a derivation of Equation (32) see van Setten et al. (2013). We introduce the following notation for the (s, s)-diagonal matrix elements of the self-energy,

Σsσ(ω)=ϕsσ0|Σσ(ω)|ϕsσ0.    (33)

The same notation is also used for matrix elements of the mean-field potential vsσMF=ϕsσ0|vσMF|ϕsσ07.

In the literature, G0 is often referred to as the “non-interacting” Green's function. However, this is technically only correct if G0 is constructed from an initial calculation based on ĥ0. This is the definition of G0 in formal many-body theory. However, often times in the theoretical literature, the Hartree potential is included in the G0 solution and excluded from the self-energy. This is the case of starting the calculation from ĥ in Equation (19). For G0W0 in practice, we usually start from ĥMF, which implies that we start from a mean-field Green's function rather than a non-interacting one. Conceptually, such a mean-field G0 is closer to the interacting G than the true G0. This is precisely why the mean-field G0 serves as such a useful starting point for GW calculations—it is closer to a self-consistent solution for G than a true non-interacting G0 is. When consulting literature references, keep in mind that GW calculations most likely refer to a mean-field G0.

4.2. Procedure

G0W0 calculations are usually performed on top of KS-DFT or HF calculations. A flowchart for a typical G0W0 calculation starting from a KS-DFT Hamiltonian is shown in Figure 10. Note that details of the flowchart depend on the treatment of the frequency dependence discussed in section 4.3. Figure 10 starts with the KS energies {ϵsKS}, KS orbitals {ϕsKS} and the exchange-correlation potential vxc from a DFT calculation. The exchange part of the self-energy Σsx is directly computed from the DFT orbitals. For the correlation term Σsc, the frequency integral over G0 and W0 must be computed, see Equation (29). If the integral is evaluated numerically, W0 is computed for a set of frequencies {ω}. The procedure to obtain W0 is as follows: First, the irreducible polarizability χ0 (Equation (28)) is computed with the KS energies and orbitals. Second, χ0 is used to calculate the dielectric function ε (Equation (26)). From the inverse of ε and the bare Coulomb interaction v, we finally obtain the correlation part of the screened Coulomb interaction, see Equations (25) and (30).


Figure 10. Flowchart for a G0W0 calculation starting from a KS-DFT calculation. The KS energies {ϵsKS} and orbitals {ϕsKS} are used as input for the G0W0 calculation. For the full expressions of χ0, ε and W0c see Equations (25–28) and (30). The spin has been omitted for simplicity.

Since the QP energies appear on both sides of Equation (22), an iterative procedure is required. More precisely, the correlation term of the self-energy depends on ϵs and must be updated at each step. Note that only G0 is a function of the QP energy, while W0c depends solely on the frequencies of the integration grid. Therefore, W0c can be pre-computed before entering the QP cycle, as displayed in Figure 10.

The correlation self-energy Σsc is a complex quantity. However, the imaginary part of Σsc is generally small for frequencies around the QP energy8, see Figures 11A,B, where Σsc(ω) is plotted for the highest occupied molecular orbital (HOMO) of the water molecule. To solve Equation (22), often only the real part of Σsc is used, which simplifies the matrix algebra to real operations and reduces the computational cost.


Figure 11. (A) Real and (B) imaginary part of the self-energy Σc(ω). Displayed is the diagonal matrix element Σsc=s|Σc(ω)|s for the HOMO of the water molecule. The gray-dashed line at ≈ − 12.0 eV indicates the QP solution ϵs. (C) Spectral function Ass(ω) computed from Equation (37). The PBE functional is used as starting point in combination with the cc-pV4Z basis set. Further computational details are given in Appendix C.

A common technique to avoid the re-calculation of the self-energy at each iteration step of the QP cycle is the linearization of Equation (22) (Giantomassi et al., 2011; Liu et al., 2016; Wilhelm et al., 2016). Assuming that the difference between QP and mean-field energies is relatively small, the matrix elements Σsc can be Taylor expanded to first-order around ϵs0:

ϵs=ϵs0+Zsϕs0|Σ(ϵs0)vMF|ϕs0    (34)
Zs=[1ddωϕs0|Σ(ω)|ϕs0ω=ϵs0]1.    (35)

The self-energy matrix elements are now only required at the mean-field energies ϵs0. Zs is known as the QP renormalization factor, because it measures how much spectral weight the QP peak carries (see also Equation (11) in section 2.2). The QP solution (main peak) is characterized by large Zs values, which lie around 0.7 to 0.8 for simple insulators, semiconductors and metals (Aulbur et al., 2000; Laasner, 2014) and around 0.9 for the molecules in Figure 12. Small Zs values indicate satellite features.


Figure 12. Error introduced by linearizing the QP equation, ΔZshot=|ϵsiter-ϵsZshot|, where ϵsiter has been obtained from the iterative procedure and ϵsZshot from Equation (34). “HOMO-x” indicates deeper valence states. The PBE functional is used as starting point in combination with the cc-pV4Z basis set. Further computational details are given in Appendix C.

The linearization error depends on the state s. The deviation from the full iterative solution usually is in the range of 0.1 eV for the HOMO, as shown in Figure 12 for a set of small molecules. The Taylor expansion of the QP equation becomes less and less accurate for larger binding energies because the absolute distance between DFT eigenvalues and QP energies increases (i.e., the G0W0 correction increases). For the deeper valence states, the linearization error is already as large as 0.5 eV (see Figure 12).

Another alternative to iterating the QP equation is to find a graphical solution. As shown in Figure 11A, the real part of the self-energy matrix elements is computed and plotted on a fine grid of real frequencies {ω} around the expected solution. All intersections of the straight line ω+vsxc-Σsx-ϵs0 with ReΣsc(ω) are then solutions of Equation (22). The intersection with the largest spectral weight Zs is the QP solution and is characterized by a small slope of ReΣsc.

Another way to calculate the QP excitations is to compute the diagonal elements of the spectral function, Ass(ω), for a set of frequencies as shown in Figure 11C. This is the most accurate procedure to obtain QP energies among the methods discussed here. Ass(ω) is computed from the complex self-energy (Σsc=ReΣsc+iImΣsc),

Ass(ω)=1πImGss(ω)sgn(EFω)    (36)
=1πIm[(ωϵs0(Σsc(ω)+Σsxvsxc))1] ×sgn(EFω),    (37)

where we employed Equation (3), the Dyson equation, G = G0 + G0ΣG, and used only the diagonal matrix elements of Σ. Figure 11C confirms that the solution at around ≈ −12.0 eV is the main solution. The spectral weight of the other solutions, e.g., the satellite peaks in the frequency range −30 to −25 eV, is indeed very small.

The aforementioned iterative procedure is computationally far more efficient than the graphical solution or the calculation of Ass(ω). The number of required QP cycles NQP typically ranges between 5 and 15 and the self-energy only has to be computed for NQP many frequencies. However, the spectral function takes also the imaginary part of the self-energy into account. This is essential for the accurate computation of satellite features in the GW spectrum (Zhou et al., 2015; Reining, 2017). Satellites fall usually in a region where ReΣsc has poles, as demonstrated in Figure 11A. In these regions, the imaginary part ImΣsc exhibits complementary peaks and is non-zero (Kramers-Kronig relation), see Figure 11B. Note that the graphical solution indicates the expected range of the satellite peaks, but does not predict their positions accurately because the imaginary part is omitted.

4.3. Frequency Treatment

The frequency integration in Equation (23) is one of the major difficulties in a G0W0 calculation since both functions that are integrated, G0 and W0, have poles infinitesimally above and below the real frequency axis. In principle, a numerical integration of Equation (23) is possible, but potentially unstable since the integrand needs to be evaluated in regions in which it is ill-behaved. However, a toolbox of approximate and exact alternatives is available. The most frequently used methods are summarized in the following.

4.3.1. Plasmon-Pole Models

The simplest way to calculate the frequency integral is to approximate the frequency dependence of the dielectric function ε and thus the screened Coulomb interaction W0 by a plasmon pole model (PPM) (Hybertsen and Louie, 1986). The PPM approximation takes advantage of the fact that ε−1 is usually dominated by a pole at the plasma frequency ωp (Hybertsen and Louie, 1986). This pole corresponds to a collective charge-neutral excitation (a plasmon) in the material. Assuming that only one plasmon branch is excited, the shape of ε can be modeled by a single-pole function

Reε1(ω)1+Ω2ω2ω~2,    (38)

where Ω and ω~ are two parameters in the model, whose squares are proportional to ωp2, see Giantomassi et al. (2011). ε, Ω and ω~ are matrices typically expressed in a plane wave basis because PPMs are mostly used for periodic systems. Note that Equation (38) holds for each matrix element and that we take the square of the matrix elements in Equation (38) and not the square of the matrix itself. Using a model function for ε−1, the expression for W0 is greatly simplified resulting in an analytic expression for the self-energy, see Deslippe et al. (2012).

The two parameters, Ω and ω~, can be determined in several ways leading to different flavors of the PPM approximation (Giantomassi et al., 2011). The most common PPMs are the Hybertsen-Louie (HL) (Hybertsen and Louie, 1986) and the Godby-Needs (GN) (Godby and Needs, 1989) model. The parameters in the HL model are obtained by requiring that the PPM reproduces the value of ε−1 in the static limit (ω = 0) and that the so-called f-sum rule is fulfilled. The f-sum rule is a generalized frequency sum rule relating the imaginary part of ε−1 to ωp and the electron density in reciprocal space (Johnson, 1974). In the HL model, the low and high real frequency limits are exact and ε has to be calculated explicitly only at ω = 0. The parameters of the GN PPM are determined by calculating ε at ω = 0 and an imaginary frequency point iωp, where ωp is typically chosen to be close to the plasma frequency ωp. The latter corresponds to the energy of the plasmon peak in the electron energy loss spectra (EELS) and can be obtained from experiment. Alternatively, ωp follows from the average electronic density ρ0 per volume, wp=4πρ0 (Giantomassi et al., 2011).

A comparison between different PPMs, namely the HL, GN, Linden-Horsch (von der Linden and Horsch, 1988), and Engel-Farid (Engel and Farid, 1993) model, can be found in Larson et al. (2013) and Stankovski et al. (2011). There it was shown that the GN model best reproduces the inverse dielectric function and the corresponding QP energies of reference calculations with an exact full-frequency treatment, such as the contour deformation discussed in section 4.3.2. However, even the accuracy of the GN model decreases further away from the Fermi energy, i.e., for low-lying occupied and high-lying unoccupied states (Cazzaniga, 2012; Laasner, 2014).

While PPMs made the first G0W0 calculations tractable (Hybertsen and Louie, 1985, 1986), full frequency calculations are now the norm, because the effects of the plasmon-pole approximation on the overall accuracy of the calculation are often hard to judge (Stankovski et al., 2011; Miglio et al., 2012). Moreover, the imaginary part of the self-energy becomes non-zero only at the plasmon poles, which implies that QP lifetimes cannot properly be calculated with PPMs, see Equation (98) and section 6.3. However, PPMs are still used in large scale G0W0 calculations (Deslippe et al., 2012), for example for solids (Jain et al., 2011; Reyes-Lillo et al., 2016), surfaces (Löser et al., 2012), 2D materials (Dvorak and Wu, 2015; Qiu et al., 2016; Drüppel et al., 2018), graphene nanoribbons (Wang et al., 2016; Talirz et al., 2017) or polymers (Hogan et al., 2013; Lüder et al., 2016).

The application of PPMs to molecules is conceptually less straightforward because the dominant charge neutral excitations in molecules are not necessarily collective. This raises the question of how to define the plasma frequency ωp of a molecule. Nevertheless, PPMs have also been used in benchmark studies for molecules, where mean absolute deviations of 0.5 eV from accurate frequency integration methods were reported (van Setten et al., 2015).

The plasmon-pole model can be extended to an arbitrary number of poles, as proposed by Rehr and coworkers (Soininen et al., 2003, 2005; Kas et al., 2007). If many frequencies are required to determine the parameters in the model, the computational cost for the evaluation of Σ is not necessarily reduced compared to full-frequency methods. However, multi-pole models are also well-defined for finite systems since the existence of a distinct plasmon peak is no longer an inherent assumption of the model (Kas et al., 2007).

4.3.2. Contour Deformation

The contour deformation (CD) approach is a widely used, full-frequency integration technique for the calculation of Σc(ω) (Godby et al., 1988; Lebègue et al., 2003; Kotani et al., 2007a; Gonze et al., 2009; Blase et al., 2011; Govoni and Galli, 2015; Golze et al., 2018). In the CD approach, the real-frequency integration is carried out by using the contour integral, see Figure 13. By extending the integrand to the complex plane, the numerically unstable integration along the real-frequency axis, where the poles of G0 and W0 are located, is avoided.


Figure 13. Contour deformation technique: Integration paths in the complex plane to evaluate Σc(ω). Γ+ and Γ are the integration contours, which are chosen such that the poles of G0, but not the poles of W0 are enclosed. Γ+ encircles the upper right and Γ the lower left part of the complex plane. ω′ denotes frequencies of the integration grid and ω the frequency at which Σc is calculated.

The integral along the contours shown in Figure 13 has four terms: an integral along the real (Re) and the imaginary axis (Im) and along the arcs.

i2π dωG0(ω+ω)W0(ω) =Re +Im +arcΓ+ +arcΓ    (39)

The contour integral is evaluated by taking the contours to infinity, which implies that the radius of the arcs is infinite. For infinitely large ω′, G0(ω+ω)W0(ω) vanishes and the integral along the arcs of contours Γ+ and Γ is zero. Therefore, we can compute the real-frequency integral by subtracting the imaginary-frequency integral from the contour integral. After rearranging Equation (39) and using Equation (23), we obtain

Σ(r,r,ω)=i2π dωG0(r,r,ω+ω)W0(r,r,ω)   12πdωG0(r,r,ω+iω)W0(r,r,iω),    (40)

where the second term is the integral along the imaginary axis.

The contours Γ+ and Γ are chosen such that only the poles of G0 fall into DΓ+ and DΓ-, which denote the subsets of the complex plane encircled by Γ+ and Γ, respectively. The location of the poles of G0(ω+ω) depends on the frequency ω at which the self-energy is computed. Recalling Equation (24), the poles of G0 lie at the complex frequencies

ωm=ϵm0ω+iη sgn(EFϵm0).    (41)

For ω < EF, these poles can enter only DΓ+ and must arise from occupied states. Our example in Figure 13 displays a case were ω < ϵ(HOMO − 1). Two poles, namely [ϵ(HOMO) − ω] and [ϵ(HOMO − 1)−ω], fall into DΓ+. For an even smaller ω, more poles from deeper occupied states will shift into DΓ+. Conversely, for ω > EF, the poles from the unoccupied states will enter DΓ-.

We can now calculate the residues of the poles that are in DΓ+ or DΓ-. Employing the residue theorem, the contour integral is then replaced by a sum over these residues:

i2π dωG0(ω+ω)W0(ω) =ωmDΓ+Res{G0(ω+ω)W0(ω),ωm} +ωmDΓRes{G0(ω+ω)W0(ω),ωm}.    (42)

The integral along the imaginary frequency axis is smooth (Rieger et al., 1999; Giantomassi et al., 2011) and the integration is performed numerically. The size of the frequency grid for the numerical integration needs to be carefully converged. For more details and a derivation of the final CD equations see e.g., Golze et al. (2018) or Govoni and Galli (2015).

4.3.3. Analytic Continuation

Analytic continuation (AC) from the imaginary to the real frequency axis is another method in our toolbox that enables an integration over the full-frequency range. The AC technique exploits the fact that the integral of the self-energy along the imaginary frequency axis,

Σc(r,r,iω)=12πdωG0(r,r,iω+iω) ×W0(r,r,iω),    (43)

is smooth and easy to evaluate, unlike the integral along the real-frequency axis. However, the QP energies and spectral functions are measured for real frequencies. To return from the imaginary to the real frequency axis, the procedure is as follows: The self-energy is first calculated for a set of imaginary frequencies {} and then continued to the real-frequency axis by fitting the matrix elements Σsc(iω) to a multipole model. A common analytic form is, e.g., the so-called 2-pole-model (Rojas et al., 1995; Rieger et al., 1999)

Σsc(iω)j=12as,jiω+bs,j+as,0,    (44)

which has been widely used for G0W0 calculations of materials (Rieger et al., 1999; Friedrich et al., 2010; Pham et al., 2013) and molecules (Ke, 2011; van Setten et al., 2015; Wilhelm et al., 2016). The unknown complex coefficients as, j and bs, j are determined by a nonlinear least-squares fit. From the identity theorem of complex analysis we know that the analytic form of a complex differentiable function on the real and imaginary axis are identical. Therefore, we can finally calculate the self-energy in the real-frequency domain by replacing with ω in Equation (44).

An alternative multipole model function is the popular Padé approximant (van Setten et al., 2015; Liu et al., 2016; Wilhelm et al., 2018), which is more flexible but contains more parameters. In the Padé approximation, the complex fitting coefficients are not obtained by a nonlinear least-squares fit, but recursively from the matrix elements Σsc(iω) and the imaginary frequencies {} (Vidberg and Serene, 1977).

In Figure 14 we compare the real self-energy matrix elements ReΣsc obtained from the AC approach to an implementation on the real-frequency axis such as the CD method. The Padé approximation reproduces the self-energy exactly in the frequency range around the QP energy of the HOMO. The deviation in the HOMO-QP energy is smaller than 10−4 eV with respect to the CD results. By using a Padé approximant with a large number of parameters, even some features of the pole structure at higher and lower frequencies are reproduced, as shown in Figure 14. The 2-pole model, on the contrary, is significantly less accurate and yields an error of around 0.1 eV in the first ionization potential. For a comprehensive comparison between the Padé and 2-pole model see van Setten et al. (2015).


Figure 14. G0W0@PBE self-energy matrix elements for the HOMO of benzene obtained with different frequency integration techniques: contour deformation (CD) and analytic continuation (AC) using the Padé model with 128 parameters and the 2-pole model. See Appendix C for further computational details.

The reliability of the AC approach is limited to valence excitations because the self-energy structure of deeper states shows poles closer to the QP solution. Our recent work (Golze et al., 2018) showed that the AC technique fails drastically to describe the complicated features of the self-energy for core states resulting in errors of 10-20 eV for the core-level binding energies. Furthermore, satellite features are difficult to obtain. As discussed in section 4.2, satellites lie in regions in which ReΣsc has poles. As evident from Figure 14, these poles can only partly be reproduced by the AC.

The convergence parameters for the AC approach are the number of frequency points {}, for which the self-energy is computed, and the size of the frequency grid for the numerical integration over ω′. In practice, the same grid is often employed for {} and {′} (Ke, 2011; Wilhelm et al., 2016).

4.3.4. Fully Analytic Approach

The integral in Equation (23) can be carried out fully analytically. In this case, the Adler-Wiser sum-over-states representation of the polarizability introduced in section 4.1 is not used. Instead we start from the reducible polarizability χ(ω). In the spectral representation, χ(ω) is given as sum of its poles n in the complex plane

χ(r,r,ω)=         nρn(r)ρn*(r)(1ω+iηΩn1ωiη+Ωn).    (45)

The pole positions Ωn correspond to charge neutral excitation energies and ρn(r) denotes transition densities. Equation (45) would be exact for the exact Ωn and ρn(r). Both quantities are obtained by solving a conventional eigenvalue problem. The equations that are solved are identical to the Casida equations (Casida, 1995b), except that for G0W0 the exchange-correlation kernel is omitted (otherwise it would be time-dependent density-functional theory).

The reducible polarizability χ(ω) can then be expanded in terms of χ0 in a Dyson series

χ(ω)=χ0(ω)+χ0(ω)vχ(ω)    (46)

and we can thus rewrite W0 given in Equation (25) in terms of χ(ω)

W0(r,r,ω)=v(r,r)+ drdrv(r,r)χ(r,r,ω)v(r,r).    (47)

Inserting Equation (45) yields a pole expansion for W0. The self-energy integral can then be solved analytically and we obtain a closed expression for Σsc(ω):

Σsc(ω)=mnϕs0ϕm0|Pn|ϕm0ϕs0ωϵm0+(Ωniη)sgn(EFϵm0),    (48)

where Pn(r,r)=ρn(r)ρn(r). More precisely, Σsc(ω) also becomes a sum over the poles Ωn. Equation (48) is therefore similar to the PPM approximation, except that we sum over the exact poles of W0 and not over the poles of a W0-model function. A detailed description of the fully-analytic frequency treatment can be found in van Setten et al. (2013). Equivalent expressions are also given in Hedin's review article from 1999 (Hedin, 1999) and were applied in Tiago and Chelikowsky (2006), Bruneval (2012), and Bruneval et al. (2016).

4.3.5. Comparison of Accuracy and Computational Cost

In the previous sections three full-frequency integration techniques have been introduced: the CD, the AC and the fully-analytic approach. The CD and fully analytic method compute the self-energy directly for real frequencies. By design, the fully analytic approach is in principle the most exact one since it is parameter-free, except for the dependence on the basis set and the broadening parameter η. However, the same accuracy can already be achieved with the CD using a moderately-sized numerical integration grid for the imaginary frequency term (Golze et al., 2018). In the AC approach the self-energy is calculated on the imaginary frequency axis, which is fairly featureless. The accuracy of the AC approach depends on the features of the self-energy on the real axis and on the flexibility of the model function, which continues the self-energy to real frequencies.

Generally, QP energies of valence states are well reproduced (van Setten et al., 2015), while the AC is likely to fail for deeper states as discussed in section 4.3.3. In the PPM approximation, the self-energy integral is simplified by introducing a model function for W0. The accuracy is therefore determined by the chosen model function and generally difficult to estimate.

The fully-analytic approach is the computationally most expensive method in our toolbox since solving the eigenvalue problem to obtain the poles of χ is an O(N6) step, where N defines the size of the system. The scaling of the CD and AC approach is generally lower, but depends on the details of the implementation. The CD method requires more computational resources than the AC methods due to the additional sum over the residues of the poles of G0. The overhead is relatively small for QP energies of valence states, but increases for deeper states due to the steady increase of the number of residues. The PPM is computationally the most efficient method, because the dielectric function ε used to compute W0 has to be calculated only at a few frequency points to determine the parameters of the PPM.

4.4. Basis Sets

In any GW implementation, the QP wave functions {ψs} and also the mean-field orbitals {ϕs0} are expanded in a set of normalized basis functions {φj}. Since in G0W0 the QP wave functions are approximated by the KS-DFT or HF ones, we expand in practice only ϕs0,

ϕs0(r)=jcsjφj(r),    (49)

where csj are the expansion coefficients that have to be determined. Performing the G0W0 calculation in a basis transforms the expression for W0, ε, and χ0 into matrix equations suitable for implementation in computer codes. The basis set choice is often guided by the type of system under investigation. In the following we will introduce the most common basis sets with brief comments on their suitability.

4.4.1. Plane Waves

For periodic systems, the energy spacing between discrete energy levels can vanish, in which case the single-particle eigenvalues form bands. According to Bloch's theorem (Bloch, 1929), the single-particle states can be written as Bloch waves

ϕnk0(r)=eik·runk(r),    (50)

where k is a wave vector in the first Brillouin zone and n is the band index. The index s that we had used in Equation (49) and throughout to label states now becomes the compound index nk. The functions unk(r) have the periodicity of the lattice and can be expanded in plane waves {φG},

unk(r)=Gcnk(G)φG(r)    (51)
φG(r)=1ΩeiG·r    (52)

where Ω is the volume of the periodic cell and G is a reciprocal lattice vector. Reciprocal lattice vectors G are given by G · t = 2πn, where n is a positive integer and t is a translation vector of the unit cell. G2 is directly proportional to the kinetic energy E of a free electron. The size of the basis set is characterized by the largest G vector and is usually given in terms of the energy E that corresponds to the largest reciprocal lattice vector, E=Gmax2/2. All G vectors with equal or smaller energies are included in the basis set.

The first GW calculations (Hybertsen and Louie, 1985, 1986; Godby et al., 1986) were performed for solids with plane wave basis sets. Also today plane waves are common in state-of-the-art GW implementations, see Table 1 for a list of GW codes. The real-space representation of G0, W0, ε, and χ0 given in Equations (24–28) can be easily transformed into a basis of plane waves by Fourier transforms. For expressions of these quantities in plane waves see, e.g., Hüser et al. (2013b).


Table 1. Selection of G0W0 codes and large program packages with G0W0 implementations and corresponding basis sets.

Plane wave basis sets are suitable for describing the slowly varying electron density in the valence region, where only the valence orbitals are non-zero. However, the valence wave functions tend to oscillate rapidly close to the nuclei due to orthogonality constraint with respect to the core orbitals. Representing these oscillations requires a large number of plane waves. Plane waves are therefore used in combination with pseudopotentials or the projector-augmented-wave methods (Blöchl, 1994) to approximate the effect of the core electrons. We will introduce the pseudopotential concept in section 4.4.4 and return to plane waves in the context of the projector augmented wave scheme in section 4.4.5.

4.4.2. Localized Basis Sets

While plane waves are mostly used for periodic systems, they can in principle also be used for finite systems by placing, e.g., the molecule in a sufficiently large unit cell to avoid spurious interactions with the neighboring cells. However, large unit cells require a very large plane wave basis set and are therefore computationally expensive. Molecular systems can be more efficiently described by atom-centered localized basis sets.

The most common basis functions of this type are Gaussian basis sets

φα,l,m(r)=NlrlYl,m(θ,ϕ)eαr2,    (53)

where Nl is a normalization constant and Yl, m(θ, ϕ) are spherical harmonic functions given in spherical coordinates (r, θ, ϕ). A Gaussian type orbital is characterized by the exponent α and the angular and magnetic quantum numbers l and m, which are dictated by the basis set selection. The design of Gaussian basis sets requires careful optimization regarding the number of functions, their respective angular momentum and exponents α. In quantum chemistry, Gaussian basis sets are widely used and ample experience exists in designing suitable basis sets for correlated methods such as coupled cluster theory. These Gaussian basis sets can then also be used in GW calculations.

Another type of localized basis functions used in GW calculations are numeric atom-centered orbitals (NAOs),

φl,m,μ(r)=Nluμ(r)rYl,m(θ,ϕ)    (54)

where uμ(r) are radial functions that are not restricted to any particular shape. The radial part of NAOs is tabulated on dense grids and is fully flexible. Gaussian radial functions can be considered as special types of this general NAO form.

Slater type functions, which posses an exponential decay at long range and a cusp at the position of the nuclei, have been also used in GW calculations (Stan et al., 2009). However, this basis set type is less common.

Local basis functions, in particular NAOs that derive from atomic orbitals, are well suited to describe rapid oscillations of wave functions near the nucleus. They are therefore the obvious choice for QP calculations of core and semi-core states.

4.4.3. Augmented Basis Sets

Augmented plane waves (APW) are another basis set type that includes the rapidly varying oscillations near the nuclei. APW methods use the so-called muffin tin approximation, which is a physically motivated approximation to the shape of the potential in solid state systems (Slater, 1937; Martin, 2004). The shape of the potential resembles a muffin tin: it is peaked at the nuclei and predominantly spherical close to it, while it is flat in between. Therefore, real space is partitioned into non-overlapping (muffin-tin) spheres ΩMT, a centered around each nuclei a and interstitial regions ΩI between these spheres. The valence wave functions are then expanded in localized NAO-like functions (Equation (54)) inside the spheres and plane waves in the interstitial regions.

By construction, the APW basis sets produce wave functions with a discontinuity in the first derivative at the muffin-tin boundaries. The linear APW (LAPW) was proposed to guarantee that the solution in the muffin-tin matches continuously and differentiably onto the plane wave part in the interstitial region (Andersen, 1975). With this extension, the explicit form of the LAPW basis functions is


where ula(r) and its derivative u˙la(r) are radial functions centered at the atom a. The coefficients Alm and Blm are determined such that continuity in value and derivative of the basis functions at the muffin-tin boundaries is ensured.

LAPW basis sets can be extended by additional local orbitals, LAPW+lo, that are completely localized in the muffin-tin spheres and go to zero at the boundaries. Inclusion of such local orbitals significantly improves the variational freedom, e.g., the description of d and f electrons (Singh, 1991). It has furthermore been shown that these local orbitals are particularly important for the unoccupied state convergence in GW calculations (Friedrich et al., 2011; Jiang and Blaha, 2016; Jiang, 2018).

A general form of the LAPW method are full-potential LAPW (FLAPW) methods that make no approximations on the shape of the potential (Wimmer et al., 1981) and which are nowadays standard in LAPW codes. Recently a number of FLAPW GW codes have emerged (Friedrich et al., 2006, 2010, 2012; Jiang et al., 2013; Gulans et al., 2014).

Linear muffin-tin orbital (LMTO) schemes are very similar to LAPW basis sets, except that the basis functions in the interstitial region are not plane waves (Andersen, 1975), but for example smooth Hankel functions (Methfessel et al., 2000).

In these augmented basis sets it is straightforward to include core and semicore states in the Green's function G0 (Equation (24)) and the polarizability χ0 (Equation (28)) and therefore in the self-energy. This, in principle, improves the description of QP excitations of valence states and band gaps, even though it has been found that the difference to carefully adjusted plane wave-based projector augmented-wave (PAW) calculations (see section 4.4.5) is typically less than 100 meV (Nabok et al., 2016). However, the same study reported larger differences for deep-lying and very localized d and f states (Nabok et al., 2016). Core excitations are in principle also accessible with FLAPW basis sets. However, these have not been thoroughly investigated yet.

For local and semi-local DFT functionals, the (F)LAPW basis sets have become the ultimate accuracy reference, closely followed by NAOs (Lejaeghere et al., 2014, 2016). For G0W0, first steps in systematically benchmarking solids were made only recently (van Setten et al., 2017). For molecules, G0W0 benchmark calculations emerged during the last years and we will discuss them in section 9.3. The jury is therefore still out on which basis set is most accurate for solids.

4.4.4. Pseudopotentials

GW calculations can be grouped in two categories: those that take all electrons of the system into consideration and those that partition into valence and core electrons. In this latter case, only the valence electrons enter the GW (and the preceding DFT) calculation explicitly, whereas the effect of the core electrons is taken into account only indirectly, for example through a pseudopotential. Such core-valence partitioning is motivated by the observation that deep core states are relatively inert and do not contribute to chemical bonding. The advantage of using a partitioning scheme is that the electron number in the GW calculation is reduced, which decreases the computational cost. An obvious disadvantage is that the core electrons may have an effect on the valence electrons, which will be difficult to include appropriately in the GW calculation and then may lead to incorrect results.

Pseudopotentials have been the default way to partition electrons (Martin, 2004; Marx and Hutter, 2009). In a pseudopotential, the core electrons are removed and the Coulomb potential of the nucleus and the inner-shell electrons is replaced by a smooth effective potential that acts on the valence electrons. The potentials are generated from calculations of isolated atoms. They are constructed such that the wave function of the valence electrons match those of an all-electron calculation outside the core region or outside a chosen radius around the nuclei. Inside the core region, the functions are smooth and nodeless. Additional norm-conservation criteria (Hamann et al., 1979; Bachelet et al., 1982), which preserve the orthonormality condition for the pseudo wave function, are usually applied (Troullier and Martins, 1991; Goedecker et al., 1996; Fuchs and Scheffler, 1999). The resulting potential is finite at the origin of the atom and shallow. Pseudopotentials are mostly used in combination with plane waves since the smooth and shallow potentials greatly reduce the required plane wave cutoff and make plane wave GW calculations with these basis sets feasible. In addition, pseudopotentials have been used for GW calculations with localized functions to reduce the basis set size (Blase et al., 2011; Wilhelm et al., 2016).

The majority of pseudopotential development took place in DFT (Marx and Hutter, 2009). Optimizing the parameters in the pseudopotential to ensure transferability is a complex task and requires thorough testing (Shirley et al., 1990; Goedecker et al., 1996). Transferability means that one and the same pseudopotential should be adequate for an atom in different chemical environments. The parameters of pseudopotentials are precomputed, similar to localized basis sets, and then tabulated for download in libraries like the Pseudo-Dojo (García et al., 2018; van Setten et al., 2018).

In GW, the consistency between pseudopotential and all-electron calculations will almost inevitably be violated. To be fully consistent, the DFT pseudopotentials would have to be cast aside and GW pseudopotentials be used. However, no such GW pseudopotentials have been developed until now, due to the complexity of the GW self-energy, which does not lend itself easily to pseudoization. Even if we had GW pseudopotentials, they would then have to first be used in the preceding DFT calculation, in which they would break the DFT core-valence consistency. Unless we perform fully self-consistent GW calculations, we are stuck with an inconsistency dilemma.

Early efforts toward GW pseudopotentials introduced core polarization effects into DFT pseudopotentials (Shirley and Martin, 1993a; Lee and Needs, 2003). By extending the GW formalism to include core contributions in the dielectric screening and the self-energy, such core-polarization potentials have also been tried successfully in the GW method (Shirley et al., 1997). However, developments in this direction did not continue. The default procedure today for plane wave GW codes is to use only well tested DFT pseudopotentials for the required elements (Govoni and Galli, 2018). Care has to be taken that the scattering states (i.e., the unoccupied states) of the pseudopotential are described well and do not introduce ghost states (Gonze et al., 1990). If no good pseudopotentials are available, it is recommended to either generate customized pseudopotentials, use the PAW method or employ genuine all-electron basis sets. Pseudopotential approaches have to be employed with care in particular for materials with localized d and f electrons. Specific issues in GW calculation of these materials are discussed in section 6.

4.4.5. Projector Augmented-Wave Method (PAW)

The PAW method is commonly used in plane wave G0W0 implementations, see Table 1. It enables computational feasibility and ensures transferability between different chemical environments. The PAW method has been derived by Blöchl combining ideas from the pseudopotential method and LAPW basis sets (Blöchl, 1994). The idea is to express the KS all-electron wave function ϕs0 for state s in terms of a smooth auxiliary function ϕ~s0 and correction terms, which restore the oscillating behavior in the core region. Note that for Bloch waves, the label s contains the k and band index n.

To construct ϕ~s0, we define a linear transformation T^ which establishes a connection between ϕs0 and ϕ~s0,

|ϕs0=T^|ϕ~s0.    (55)

Since the all-electron wave function is already smooth at a certain distance from the nuclei, we partition the space similarly to LAPW schemes: in atom-specific augmentation regions Ωa around the nuclei, where a is the atom index, and an interstitial region ΩI. The augmentation regions are characterized by the cutoff radii rca, which should be chosen such that the augmentation spheres do not overlap. Outside the augmentation regions, ϕ~s0 should be identical to the all-electron wave function. T^ should thus modify ϕs0 only in Ωa and we define

T^=1+aT^a,    (56)

where the atom-centered transformation, T^a, acts only within Ωa.

The transformation operator is derived by introducing atom-centered functions as in LAPW, which is described in detail in Martin (2004) and Rostgaard (2009). The all-electron wave function can then be rewritten as

ϕs0(r)=ϕ˜s0(r)+a(ϕsa(r)ϕ˜sa(r)),    (57)

where the atom-centered hard and smooth auxiliary wave functions are denoted by ϕsa and ϕ~sa, respectively. “Hard” refers to rapidly varying functions in the core region. The concept of the PAW scheme is visualized in a simplified way in Figure 15. By adding ϕa to ϕ~s0 we obtain the oscillating behavior in the core region, but we have to subtract the smooth function ϕ~a to cancel the contribution of ϕ~s0 in Ωa. That implies that the following conditions must hold

ϕs0(r)=ϕ˜s0(r)ϕsa(r)=ϕ˜sa(r)}rΩI ϕs0(r)=ϕsa(r)ϕ˜s0(r)=ϕ˜sa(r)}rΩa.

The atom-centered auxiliary wave functions can be expanded in a finite set of local basis functions {φja} and {φ~ja} and a set of projector functions p~ja, where ‘~’ indicates again smooth functions. These expansions are given by

ϕsa(r)=jφja(r)p˜ja|ϕ˜s0    (58)
ϕ~sa(r)=jφ~ja(r)p~ja|ϕ~s0.    (59)

The variational object in a PAW calculation is ϕ~s0. The latter is expanded using, e.g., a plane wave basis set, for which a low energy cutoff can be used due to its smoothness. The local basis sets and projector functions needed to compute the second and third terms in Equation (57) are tabulated for each element of the periodic table. For specific choices of these basis sets see, e.g., Kresse and Joubert (1999) and Rostgaard (2009). Details regarding the practical implementation within a plane wave code are given in Kresse and Joubert (1999) and for real-space grid codes in Mortensen et al. (2005) and Enkovaara et al. (2010).


Figure 15. Schematic representation of the projector augmented wave (PAW) scheme. The all-electron wave function ϕ is constructed from the smooth auxiliary function ϕ~ and corrections from the hard and smooth atom-centered auxiliary wave functions ϕa and ϕa~, respectively.

GW calculations within the PAW schemes employ usually the frozen core approximation (Shishkin and Kresse, 2006a; Liu et al., 2016). The core states are localized at the atoms and confined in Ωa. In the frozen core approximation we assume that the KS core states are identical to the atomic core states α, i.e., ϕs0=φαa. In this approximation, the decomposition given in Equation (57) is not used for the core states. However, the effect of the core on the valence states is correctly described.

The accuracy of the expansion in Equation (57) depends on the completeness of the set of localized basis and projector functions ({φja}, {φ~ja} and {p~ja}). Achieving completeness is easy for occupied states, but practically impossible if s corresponds to a high-energy empty state, which has been discussed by Klimeŝ et al. (2014). However, for a GW calculation, many of these high-energy empty states need to be included, which is explained in more detail in section 4.5. These incompleteness issues lead to a violation of the norm-conservation for the unoccupied states, which can be the source of substantial errors, in particular for elements with d and f electrons. This error can be avoided using norm-conserving instead of the standard PAW potentials for GW calculations (Klimeŝ et al., 2014).

4.5. Basis Set Convergence

The first criteria for a reliable G0W0 calculation is that the underlying DFT or HF calculation is converged. This convergence has to be checked for all basis set types. The second convergence criteria is the size of the basis set in the G0W0 calculation itself. In quantum chemistry it is well established that correlated electronic structure methods converge slowly with respect to the number of basis functions (Kendall et al., 1992; Kutzelnigg and Morgan, 1992; Klopper et al., 1999). The same has been also observed for G0W0 (Shih et al., 2010; Friedrich et al., 2011; Bruneval, 2012; Yan et al., 2012; Bruneval and Marques, 2013; van Setten et al., 2013, 2015; Jacquemin et al., 2015; Bruneval et al., 2016; Wilhelm et al., 2016). It has been demonstrated that the convergence rate of G0W0 is comparable to other correlated methods such as second-order Møller-Plesset perturbation theory (MP2) and the coupled cluster singles, doubles and perturbative triples [CCSD(T)] method (Bruneval and Marques, 2013).

Converging G0W0 excitations within a plane wave basis set is straightforward since the cutoff energy can be continuously increased, see Figure 16A. In addition, extrapolation schemes to the complete basis set limit have been reported to reduce the computational cost (Klimeŝ et al., 2014; Maggio et al., 2017; Govoni and Galli, 2018). Conversely, for localized basis sets only a limited number of basis set sizes is available and a steady increase in size as for plane waves is not possible. Therefore, extrapolation techniques must always be used to obtain converged G0W0 energies. This is displayed in Figure 16B, where the first IP of benzene is plotted with respect to the inverse of the basis set size. Shown are results for the Dunning basis set family cc-pVnZ (n=3-6) (Dunning, 1989; Wilson et al., 1996), which was designed to smoothly reach the complete basis set limit. Increasing values of n indicate increasingly large basis sets. The convergence is smooth, but very slow, as shown in Figure 16B.


Figure 16. Basis set convergence for G0W0 calculations. (A) Convergence for a plane wave basis set. Bandgap of wurtzite ZnO dependent on the number of bands and on the corresponding cutoff energy (data from SI of Yan et al., 2012). (B) Convergence and extrapolation procedure for a localized basis set. Ionization potential (IP) for the HOMO of benzene plotted with respect to the inverse of the number of basis functions Nfunc using the cc-pVnZ basis set series. Further computational details are given in Appendix C.

The extrapolation to an infinite number of basis functions is performed by a linear regression with respect to the inverse of the total number of basis functions. The extrapolated value of 9.17 eV in Figure 16 is 0.07 eV larger than the IP obtained at the cc-pV6Z level showing that even the largest basis set cannot converge the G0W0 energies completely. This linear fitting procedure is a well-established scheme to extrapolate G0W0 energies and has been tested in large benchmark studies (van Setten et al., 2015). Alternatively, linear regression has also been performed with respect to Cn-3, where Cn is the cardinal number of the basis set, i.e., 3 for cc-pVTZ, 4 for cc-pVQZ, 5 for cc-pV5Z and 6 for cc-pV6Z. Extrapolation with respect to Cn-3 is well-established for correlated methods (Helgaker et al., 1997). The inverse of the basis set number corresponds roughly to Cn-3. The average difference between the two extrapolation schemes for G0W0 energies is indeed very small with 0.04 eV (van Setten et al., 2015).

A common misconception in the plane wave community is that the number of unoccupied states that enter a G0W0 calculation is a separate convergence parameter. The number of unoccupied states that can be resolved with a given basis set typically grows with the size of that basis set, i.e., the Hilbert space of that basis grows. This implies that more empty states enter the sums in the Green's function (Equation (24)) and the polarizability (Equation (28)). Since it is computationally expensive to generate a large number of unoccupied states in the preceding plane wave DFT or HF calculation, plane wave G0W0 practitioners reduced the number of unoccupied states that enter the GW calculation in order to save computational time (Stankovski et al., 2011; van Setten et al., 2017). Localized basis sets on the other hand are significantly smaller and typically all virtual states are computed even in DFT-only calculations. Figure 16A gives an impression of the scale for the plane wave case. It displays the convergence of the band gap of wurtzite ZnO with respect to the number of bands (states) (Yan et al., 2012). The convergence rate of ZnO is particularly slow (Friedrich et al., 2011; Stankovski et al., 2011) compared to other semiconductors, e.g., silicon (Friedrich et al., 2006). The band gap finally converges at around 30 Ha. At this point almost 3000 bands have been included in the G0W0 calculation.

While it might seem appealing to reduce the number of required unoccupied states to less than 3,000, Figure 16A illustrates that a reduction is not possible due to the slow convergence. Since the number of resolvable, unoccupied states is coupled to the plane wave cutoff (Stankovski et al., 2011; Gao et al., 2016; van Setten et al., 2017), one should always include the maximum number of bands in the G0W0 calculation for a given plane wave cutoff. Such a procedure also greatly simplifies the convergence study since only one and not two parameters need to be converged.

4.6. Elimination of Unoccupied State Summation

The complications around the virtual state convergence raised in the previous section for plane wave basis sets can be bypassed completely by eliminating the explicit summation over unoccupied states in the Green's function G0 (Equation (24)) and the polarizability χ0 (Equation (28)) (Reining et al., 1997; Wilson et al., 2008, 2009; Berger et al., 2010; Giustino et al., 2010a; Umari et al., 2010). A practical method for building a perturbation theory without explicit reliance on unoccupied states was pioneered in the context of density functional perturbation theory (DFPT) (Baroni et al., 1987; Gonze, 1995, 1997). Here we will briefly review how the DFPT concept can be transferred to G0W0. For a general introduction to the DFPT formalism see (Baroni et al., 2001).

The central object in DFPT is the response function that measures the response of a system to a perturbation ΔV. In the G0W0 context, we are interested in the response to the introduction of an additional charge to the system at point r. The additional charge perturbs the charge density of the system. The response function mediates the charge density change and the perturbation. We now wish to calculate the response function without invoking the sum over states expression introduced in Equation (28).

We start with the change in the charge density Δn(r, r′, ω) given for a spin-unpolarized system by Giustino et al. (2010a)

Δn(r,r,ω)=2ioccϕi0*(r)×                                     (Δϕi0(r,r,ω)+Δϕi0(r,r,ω)).    (60)

Here Δϕi0(r,r,±ω) is the frequency-dependent variation of the occupied mean-field state i. Instead of expanding Δϕi0(r,r,±ω) in the basis of unperturbed mean-field states ϕi0(r) it is calculated directly with the Sternheimer equation (Baroni et al., 1987; Giustino et al., 2010a)

(ĥMFϵi0±ω)Δϕi0(r,r,±ω)=(1P^occ)ΔV(r,r)ϕi0(r).    (61)

P^occ is a projection operator on the occupied states, P^occ=iocc|ϕi0ϕi0|.

Sternheimer G0W0 formalisms differ in their choice of ΔV. There are two possible choices for ΔV. The first one is to set the perturbation to the bare Coulomb interaction v(r, r′) (Reining et al., 1997)

ΔV(r,r)=v(r,r).    (62)

This choice is known as non-self-consistent Sternheimer GW. The Sternheimer G0W0 formalism is shown in Figure 17. Quantities that depend on r′ are expanded in a basis {φk(r)}, see Equation (49), and both sides of Equation (61) are projected onto φl(r). This leads to a linear set of equations with a parametric dependence on r and ±ω. Solving the Sternheimer equation for each real-space grid point r and the frequencies ±ω yields Δϕi0(r,r,±ω)) for the occupied state i. From the latter we can evaluate the induced charge density Δn. The dielectric function given in Equation (26) can be rewritten in terms of Δn (Reining et al., 1997; Lambert and Giustino, 2013)

ε(r,r,ω)=δ(r,r)Δn(r,r,ω).    (63)

W0 is then calculated from the inverse of ε according to Equation (25) as usual.


Figure 17. Non-self-consistent Sternheimer approach for obtaining W0 without empty states. Δϕ0 and v have a parametric dependence on the real space point r. Δϕ0 depends additionally on the frequency ω.

The second choice for ΔV is to set it to the screened Coulomb interaction

ΔV(r,r,ω)=W0(r,r,ω),    (64)

leading to the self-consistent Sternheimer GW formalism introduced by Giustino et al. (2010a). In this scheme, the (self-consistent) induced charge density ΔnSC generates a potential ΔVscr, which screens the bare Coulomb potential v due to the perturbative charge in the system. From ΔVscr we can directly calculate W0 as

ΔVscr(r,r,ω)=drΔnSC(r,r,ω)v(rr)    (65)
W0(r,r,ω)=v(r,r)+ΔVscr(r,r,ω).    (66)

It can be shown that Equation (66) is equivalent to Equation (25) (Giustino et al., 2010a). Since W0 appears on the right-hand side of Equation (61), it must be solved self-consistently. In the first step, W0 is initialized with v. From the solutions of the Sternheimer equation, we calculate ΔnSC, ΔVscr and finally W0.

Both schemes yield W0, but the non-self-consistent approach requires fewer steps. However, the dimensions of the dielectric matrix increase with system size and its inversion might become a bottleneck for large systems, in particular when using plane wave basis sets. In this case the self-consistent scheme might be more efficient.

The two schemes discussed so far address the elimination of empty states in W0. Removing the sum over virtual states in G0 is also possible by using a similar strategy as for W0, see Giustino et al. (2010a) for a detailed description. Once G0 and W0 have been obtained, the self-energy is composed as usual and the frequency integration is performed with the methods described in section 4.3. Sternheimer approaches have been implemented for plasmon pole models (Reining et al., 1997), the analytic continuation (Giustino et al., 2010a) and contour deformation (Govoni and Galli, 2015).

The Sternheimer GW approach is primarily used in plane wave implementations (Reining et al., 1997; Nguyen et al., 2012; Lambert and Giustino, 2013; Pham et al., 2013; Govoni and Galli, 2015). We are only aware of one non-plane wave implementation using mixed representations of real space and localized basis sets (Hübener et al., 2012a,b). As discussed in section 4.5, converging a G0W0 calculation with plane waves requires a very large number of empty states. The calculation of all empty states in the preceding DFT or HF calculation is computationally expensive and can easily become a computational and storage bottleneck. In the Sternheimer approach, the preceding DFT step is significantly simplified since only occupied states have to be calculated. For localized basis sets, no such benefit is found in DFPT or Sternheimer since the number of virtual states is typically not that large and rarely a bottleneck (Shang et al., 2018).

Sternheimer G0W0 saves not only computational time in the preceding mean-field calculation, but also by not having to carry out the sums over states in G0 and χ0. However, it concomitantly loses time in the Sternheimer iterations. To our knowledge, a detailed comparison of the computational cost to conventional G0W0 implementations has not been reported yet. To speed up Sternheimer G0W0, projection techniques for representing the dielectric matrix in an optimal, smaller basis (Wilson et al., 2008, 2009; Nguyen et al., 2012; Pham et al., 2013; Govoni and Galli, 2015) and Lanczos-chain algorithms that efficiently obtain the Sternheimer solution over a broad frequency range (Umari et al., 2010) have been developed. Furthermore, all the Sternheimer equations, that need to be solved for each r and ω, are independent from each other facilitating massively parallel implementation (Govoni and Galli, 2015; Schlipf et al., 2019).

The Sternheimer approach does not reduce the basis set size, i.e., the plane wave cutoff or equivalently the size of the real-space grid, nor does it change the formal scaling of G0W0 with respect to system size. However, it is an elegant way to facilitate easier convergence, since the temptation of converging the number of virtual states separately is removed.

A modified Sternheimer ansatz has been developed for FLAPW basis sets which accounts for response contributions outside the Hilbert space spanned by the basis set (Betzinger et al., 2012, 2013, 2015). This modified approach thus allows the basis set size to be decreased, unlike the classical Sternheimer technique. The explicit summation over unoccupied states is not completely removed, but the number of empty states needed for convergence is strongly reduced.

4.7. Starting Point Dependence and Optimization

The results of a G0W0 calculation depend on the wave functions {ϕs0} and the energies ϵs0 that are used as input for the Green's function (G0) and the screened Coulomb interaction (W0). The single-particle wave functions and energies are determined by the choice of the single-particle mean-field Hamiltonian ĥMF, e.g., by the chosen DFT functional. To denote this dependence, we will introduce the notation G0W0@startingpoint.

4.7.1. How Severe Is the Dependence on the Reference State?

Until the early 2000s, the majority of all G0W0 calculations were based on DFT calculations using the local-density approximation (LDA) or a generalized gradient approximation (GGA) (Aryasetiawan and Gunnarsson, 1998; Aulbur et al., 2000; Onida et al., 2002). With the advent of exact-exchange based DFT functionals in the solid state community and the proliferation of G0W0 codes that are based on quantum chemical codes, a more diverse range of starting points became available. It was soon realized that G0W0 can exhibit a pronounced starting-point dependence for semiconductors (Rinke et al., 2005; Fuchs et al., 2007). In the last years, the starting point dependence has also been intensively discussed for molecules (Körzdörfer and Marom, 2012; Marom et al., 2012; Bruneval and Marques, 2013; Gallandi and Körzdörfer, 2015; Caruso et al., 2016; Gallandi et al., 2016).

Figure 18 illustrates the starting-point dependence for the HOMO of the water molecule. For the underlying DFT calculations, the PBE-based hybrid (PBEh) functional family (Perdew et al., 1996b; Adamo and Barone, 1999; Ernzerhof and Scuseria, 1999) was scanned. The PBEh functional family is characterized by an adjustable fraction α of HF exchange. The exchange-correlation energy Exc is therefore α-dependent and given by

Exc=αExEX+(1α)ExPBE+EcPBE, α[0,1],    (67)

where ExEX denotes the HF exchange energy. ExPBE and EcPBE are the PBE exchange and correlation energy, respectively. To illustrate the starting point dependence in G0W0, the mixing parameter α in PBEh was varied from 0 to 1 and then a subsequent G0W0 calculation was performed. The resulting G0W0 HOMO energies shown in Figure 18 span a range of more than 1 eV. Although a 1 eV spread appears large, it is much smaller than the range of the corresponding mean-field energies ϵHOMO0 that decrease from −7 eV to −14 eV with increasing α.


Figure 18. Starting point dependence of G0W0: the left side shows the G0W0 HOMO energy of the water molecule for hybrid functional starting points with different amounts of exact exchange. The HOMO energy in self-consistent GW (scGW) is shown on the right. The dashed line marks the experimental value of 12.62 eV (Page et al., 1988; Lias and Liebman, 2003). All GW values are extrapolated to the exact basis set limit using the cc-pVnZ (n = 3–5) basis sets. Further computational details are given in Appendix C.

The strong dependence of G0W0 on the starting point can be largely attributed to over- and underscreening. From the Adler-Wiser expression for χ0 (Equation (28)) it can be deduced that the screening strength in G0W0 is inversely proportional to the eigenvalue gap of the starting point. Since HF typically overestimates gaps, it will underscreen. PBE, on the other hand, underestimates gaps and therefore overscreens.

Another source of error in the KS orbital energies is the spurious self-interaction term (Perdew and Zunger, 1981). The one-electron self-interaction error (SIE) arises from an incomplete cancellation of the unphysical electrostatic Hartree energy of an electron with itself by the exchange-correlation term. The SIE is more pronounced for localized than delocalized orbitals (Körzdörfer et al., 2009). This can be intuitively understood: an electron in a localized orbital has a stronger self-interaction because its wave function is more confined. As a result, the localized orbitals are destabilized with respect to more delocalized orbitals. This can lead to a wrong ordering of the orbital energies in the underlying DFT calculations, which carries over to the GW spectrum. The SIE can be mitigated by a larger amount of exact exchange, which also restores the correct ordering for the QP energies (Marom et al., 2011, 2012; Körzdörfer and Marom, 2012).

The G0W0 starting point dependence generally lies in the range of 1.0 eV for HOMO energies of molecules (Marom et al., 2012), but increases for deeper states. For solids, the spread can exceed 2 eV for the band gap (Fuchs et al., 2007). This beckons for a judicious choice of the starting point in G0W0 calculations or an elimination of the starting point dependence. The dependence on the preceding mean-field calculation can be eliminated or reduced by employing some form of self-consistency as discussed in section 5 or, as proposed only very recently, by replacing G0 by a renormalized singles Green's function (Jin et al., 2019). In this section we focus on the optimal choice of the starting point. The PBEh family of DFT functionals is convenient for this purpose, since one parameter (the amount of exact exchange α) governs the behavior of the starting point. Several schemes have been developed to find the optimal α value within the PBEh functional family (Körzdörfer and Marom, 2012; Atalla et al., 2013; Pinheiro et al., 2015; Dauth et al., 2016; Bois and Körzdörfer, 2017). We summarize them in the following.

4.7.2. Consistent Starting Point Scheme

Körzdörfer and Marom developed a consistent starting point (CSP) scheme that seeks a PBEh reference state (i.e., starting point) that best resembles the G0W0 spectrum. Splitting both the G0W0 self-energy and the starting hybrid functional into their respective exchange and correlation parts (see Equations (29) and (32)) allows us to rewrite the QP equation (Equation (22)) as follows (Körzdörfer and Marom, 2012; Marom, 2017)

ϵs=ϵs0+(1α)Δvsx+Δvsc    (68)
Δvsx:=ϕs0|ΣxvxPBE|ϕs0    (69)
Δvsc:=ϕs0|ΣcvcPBE|ϕs0.    (70)

vxPBE and vcPBE are the exchange and correlation part of the PBE exchange-correlation potential, respectively. The optimal α is determined so that the shift between G0W0 and PBEh for the occupied states is approximately a constant k

Δvsc+(1α)Δvsxk, socc.    (71)

If Equation (71) is satisfied, the positions of the PBEh orbital energies relative to each other are as close as possible to the G0W0 energies. In this case, the QP correction amounts to a rigid shift of the PBEh spectrum. The value of α for which Equation (71) is satisfied yields the optimal starting point in the CSP scheme. If the PBEh and the G0W0@PBEh spectrum matched perfectly, the constant k would be zero. However, in general it is not possible to find a starting point whose spectrum matches the G0W0 spectrum exactly.

For a given guess of α, Δvsx and Δvsc are calculated according to Equations (69) and (70). Δvsc is plotted as a function of Δvsx for a set of occupied states s as data points, see Figure 19. A straight line fit determines a new α which is used to calculate new DFT eigenvalues and orbitals from PBEh(αnew). From the new eigenvalues and orbitals, a new self-energy is calculated and Equations (69) and (70) are reassessed. This procedure is continued until the α of the linear fit equals the initial α. Then the optimal α has been found.


Figure 19. CSP scheme representative for a small molecule. Δvsc (Equation (70)) is plotted with respect to Δvsx (Equation (69)) for a set of occupied states s. The HOMO and HOMO-1 states are indicated. The new α value is obtained from the slope of the straight line fitted through the red symbols. Data retrieved from Körzdörfer and Marom (2012).

By construction, the PBEh(α) eigenvalues are now, up to the rigid shift k, consistent with the G0W0 spectrum. Typical CSP α values lie in the range of 0.25 − 0.30. The CSP scheme has been tested on several organic molecules that are used in organic electronics and yields good agreement with photoemission spectra in all cases (Körzdörfer and Marom, 2012; Körzdörfer et al., 2012; Marom, 2017).

4.7.3. Deviation From the Straight Line Scheme

A physically more rigorously motivated optimization scheme is based on the deviation from the straight line error (DSLE) (Dauth et al., 2016). In 1982, Perdew and co-workers showed that the total energy E of any many-electron system should change linearly when varying the electron number continuously from N to N − 1 electrons (Perdew et al., 1982),

E(f)=(1f)E(N1)+fE(N)  f[0,1].    (72)

The function E(f) is a piecewise linear function of the occupation number f, with cusps at every integer value of f, see Figure 20. Standard DFT functionals, however, violate this piecewise linearity condition and yield energies that deviate from the straight line at fractional occupation numbers f (Mori-Sánchez et al., 2006; Ruzsinszky et al., 2006; Kraisler and Kronik, 2013). The straight-line condition applies not only to DFT, but to any total energy method (we will address GW total energies in section 10). Therefore, we can use the DSLE to find an optimal starting point for G0W0.


Figure 20. Schematic representation of the straight line condition for total energies E (left) and derivatives ∂E/∂f (right). f is the occupation number. The DSLE is shown for three different cases: convex (blue), concave (red), and mixed curvature (green). Reprinted with permission from Dauth et al. (2016). Copyright (2016) by the American Physical Society.

The slope of the total energy as a function of occupation gives the Kohn-Sham eigenvalue or, in the GW case, the quasiparticle excitation energy for a given electron number,

Ef|f=1=E(N)E(N1)=ϵHOMO,N    (73)
Ef|f=0=E(N)E(N1)=ϵLUMO,N1,    (74)

where ϵHOMO, N is the QP energy of the HOMO for the neutral system (N). ϵLUMO, N−1 denotes the QP energy of the lowest unoccupied orbital (LUMO) for the charged system (N − 1). It is evident from Equations (73) and (74) that the slopes must be identical for f = 0 and f = 1 in an exact theory. In other words, a necessary condition for piecewise linearity is that the energy for removing an electron from the neutral system equals the energy for adding an electron to the positively charge system, i.e., the IP for the neutral systems and the electron affinity (EA) of the charged system should be equal. The difference

=ϵLUMO,N-1+ϵHOMO,N    (76)

should thus be zero and if it is not zero it quantifies the deviation from the straight line error ΔDSLE.

The idea is now to find a PBEh starting point for which G0W0@PBEh minimizes ΔDLSE. The optimal α value for PBEh can be found by the following procedure: We select a set of PBEh functionals with α values between 0 and 1. Then two G0W0 calculations are performed for each starting point. One for the neutral system that yields ϵHOMO, N and a second one for the cation to obtain ϵLUMO, N-1. Following Equation (76), we calculate the difference between these two energies to estimate the deviation from the straight line condition. The PBEh(α) functional that yields the smallest ΔDSLE will be the optimal starting point.

This DSLE scheme has been tested for small molecular systems, where it has been shown that the optimal α values are distributed around 0.35 − 0.40 for the first IP (Caruso et al., 2016; Dauth et al., 2016). The reported deviation from the CCSD(T) reference is smaller than 0.2 eV (Caruso et al., 2016). The drawback of the DSLE scheme is that the optimization is restricted to the HOMO. For other states, the straight line condition could still be formulated, but the corresponding G0W0 calculations could not be performed because the electron occupation function would no longer correspond to an equilibrium distribution (see section 11.3 for GW calculations out of equilibrium). If we removed an electron from a lower lying occupied state, the sums over occupied and virtual orbitals in the polarizability χ0 would no longer be rigorously defined, i.e., the energy differences in the denominator in Equation (28) would exhibit the wrong sign.

4.7.4. IP-Theorem Schemes

Several other schemes were developed that are, in spirit, similar to the CSP and DSLE optimization approaches, but are explicitly based on the ionization potential (IP) theorem. The latter states that in exact DFT the negative of the KS orbital can be strictly assigned to the first ionization potential IPHOMO (Levy et al., 1984; Almbladh and von Barth, 1985)


This statement is not true for any other KS state and not valid for approximate DFT functionals. Atalla et al. proposed a scheme that exploits the IP-theorem and minimizes the G0W0 correction for the HOMO level with respect to the amount of exact exchange α in a PBEh starting point (Atalla et al., 2013),

ΔvHOMOc+(1α)ΔvHOMOx=0.    (78)

This approach is similar to the CSP scheme in Equation (71). The difference is that we find the PBEh functional whose HOMO energy matches that of G0W0@PBEh for the same α, whereas CSP looks for the closest spectral match between PBEh and G0W0@PBEh. HOMO excitations obtained from this IP-theorem-tuned scheme agree reasonably well with CCSD(T) reference data (Bois and Körzdörfer, 2017), but are not expected to reproduce the whole excitation spectrum properly (Atalla et al., 2013). They generally yield large αs (around 0.8) and produce HOMOs and HOMO-LUMO gaps that are too large (i.e., underscreened).

Finding a PBEh(α) functional that fulfills the IP-theorem by enforcing consistency with the G0W0 spectrum is one option. An alternative approach is to IP tune the hybrid functionals themselves (Bois and Körzdörfer, 2017) by minimizing

ΔIP=|ϵHOMOKS(α)(E(N,α)E(N1,α))|    (79)

with respect to α. These IP-tuned hybrids already give accurate KS-HOMO energies. Recent benchmark studies for molecular systems showed that G0W0 corrections on top of ΔIP-tuned functionals provide spectral properties in good agreement with experiment for the whole excitation spectrum (Refaely-Abramson et al., 2012; Egger et al., 2014; Gallandi and Körzdörfer, 2015; Gallandi et al., 2016; Knight et al., 2016; Bois and Körzdörfer, 2017). In particular, EAs are well reproduced with a mean absolute devation (MAD) smaller than 0.2 eV from the CCSD(T) reference, while the MAD reported for IPHOMO is 0.1 eV (Knight et al., 2016).

4.8. Computational Complexity and Cost

Of all the GW schemes described in this review, G0W0 is the computationally most efficient one. Only the diagonal elements of the self-energy are needed and the Green's function that enters is always G0. The fully interacting G, on the contrary, depends on the full self-energy and can only be obtained by iterating Dyson's equation G = G0 + G0ΣG.

The computational complexity of G0W0 depends on the frequency integration method and design of the algorithm. The most accurate integration technique, the fully-analytic approach, requires the solution of the full Casida equations, which is an O(N6) step with respect to the system size N, see section 4.3.5. In the canonical implementation, the computational cost is reduced to O(N4). Different implementations with N4 complexity have been developed employing a variety of numerical techniques specific for the respective basis set (Shishkin and Kresse, 2006a; Blase et al., 2011; Deslippe et al., 2012; Ren et al., 2012a; Govoni and Galli, 2015; Wilhelm et al., 2016). For example, the O(N4) algorithm proposed by Ren et al. employs localized basis functions and the AC method (Ren et al., 2012a). The computational and memory costs for the four-center electron repulsion integrals (4c-ERIs) are reduced by refactoring the latter in two- and three-center Coulomb integrals using a resolution-of-the-identity (RI) approach with a so-called Coulomb metric (Vahtras et al., 1993). The accuracy of this algorithm has been validated for valence excitations and EAs by comparing to the fully analytic approach (van Setten et al., 2015). However, for core states it was recently shown that AC fails and that CD is required (Golze et al., 2018). The computational complexity of CD remains unchanged for valence states, but increases to O(N5) for the deep states.

The O(N4) scaling and overall cost of canonical GW implementations restrict the tractable system size and prohibit the study of many systems that are relevant in the chemistry and physics community, such as solid-liquid interfaces, molecules in solution, complex alloys, nanostructures or hybrid interfaces, that require large simulation cells with hundreds to thousands of atoms. To make G0W0 calculations feasible for larger systems, the scaling and computational complexity have been scrutinized. Developments have proceeded in two directions: (1) reducing the prefactor, i.e., the overall computational cost at O(N4) scaling, or (2) reducing the scaling.

The prefactor has been reduced by low-rank approximations of χ0, which map χ0 onto a smaller basis (Wilson et al., 2008, 2009; Umari et al., 2009; Govoni and Galli, 2015; Del Ben et al., 2019). Another approach is the elimination of the sum over empty orbitals in χ0 and in G0 (Giustino et al., 2010a; Umari et al., 2010; Lambert and Giustino, 2013; Pham et al., 2013; Govoni and Galli, 2015) by solving the Sternheimer equation (Sternheimer, 1954), which we discussed in section 4.6. Others developed techniques to reduce the number of unoccupied states (Bruneval and Gonze, 2008; Bruneval, 2016).

The prefactor can also be controlled by choosing an optimal basis set for the respective system under study. In the last years, several algorithms for localized basis sets have been developed (Blase et al., 2011; Ke, 2011; Ren et al., 2012a; van Setten et al., 2013; Bruneval et al., 2016; Wilhelm et al., 2016). These basis sets are generally smaller than traditional plane wave basis sets and considerably more efficient for molecules. However, the development of reliable G0W0 algorithms for periodic systems based on localized basis sets is still underway (Wilhelm and Hutter, 2017).

The reduction of the exponent to O(N3) complexity has been addressed in different ways. Foerster et al. developed a cubic-scaling G0W0 algorithm using Gaussian basis sets and exploiting locality in the electronic structure, albeit with a high prefactor (Foerster et al., 2011). Recently, two cubic-scaling algorithms have been devised (Liu et al., 2016; Wilhelm et al., 2018), which are both variants of the O(N3) GW space time method (Rojas et al., 1995). The key step of these algorithms is the computation of the irreducible polarizability in imaginary time, χ0(it) = −iG0(it)G0(−it) and the subsequent transformation to imaginary frequencies . The time-ordered non-interacting Green's function in imaginary time is given by

G0(r,r,it)={iioccϕi0(r)ϕi0*(r)exp(ϵi0t),t<0,iavirtϕa0(r)ϕa0*(r)exp(ϵa0t),t>0.    (80)

Inserting G0(it), the summation over occupied and virtual states is now decoupled in χ0(it) and can be performed separately, which is fundamental for the reduction to O(N3) complexity.

Liu et al. based their cubic-scaling algorithm on a plane wave basis set in combination with a PAW scheme and reported a linear-scaling with the number of k points used to sample the Brillouin zone (Liu et al., 2016). In combination, this paves the way for GW calculations of large periodic systems. Wilhelm et al. employed a Gaussian basis set and exploited sparse matrix algebra by using an overlap metric for the RI approximation (RI-SVS) to refactor the 4c-ERIs (Vahtras et al., 1993). The step with the largest prefactor, the computation of χ0, is reduced from O(N4) to O(N2) in this algorithm, while the other operations scale with N3. A comparison between the O(N4) algorithm developed by Ren et al. (2012a) and the low-scaling algorithm is shown in Figure 21 for graphene nanoribbons. The canonical algorithm is restricted to system sizes of less than 500 atoms, while systems with more than 1,600 atoms can be addressed with the low-scaling implementation. These are some of the largest G0W0 calculations with high accuracy and full-frequency integration reported so far. The mean absolute deviation with respect to the canonical reference implementation in FHI-aims (Ren et al., 2012a) is less than 35 meV for the GW100 test set (Wilhelm et al., 2018), which is discussed more in detail in section 9.3.


Figure 21. Scaling of state-of-the-art G0W0 implementations with respect to system size using graphene nanoribbons as a benchmark system. (A) Smallest graphene nanoribbon unit with 114 atoms. (B) Comparison of the scaling of the canonical G0W0 (Wilhelm et al., 2016) and the low-scaling implementation (Wilhelm et al., 2018). The latter requires operations of at most O(N3) complexity (red diamonds). Dashed lines represent least-square fits of exponent and prefactor. Data retrieved from Wilhelm et al. (2018). Both algorithms are implemented in the CP2K program package.

An actual linear scaling algorithm was devised within the framework of stochastic GW (Neuhauser et al., 2014) and applied to silicon clusters with 1,000 atoms. However, the verification of its general reliability is still the subject of ongoing research (Vlček et al., 2017).

4.9. Practical Guidelines

In summary, the following points should be considered when conducting G0W0 calculations:

1. Frequency integration technique

A sufficiently accurate method for the frequency integration of Equation (23) has to be chosen. The required precision depends on the system and in particular on the states of interest, see section 4.3.

2. Basis set choice

The decision should be guided by the system of interest. Localized basis sets are generally more efficient for finite systems, while plane wave G0W0 codes are currently superior for extended systems (see detailed discussion in section 4.4).

3. Basis set convergence

GW calculations have to be carefully converged with respect to basis set size. Extrapolation procedures to the complete basis set limits might be required as demonstrated in section 4.5.

4. Starting point

The QP energies depend strongly on the functional of the preceding DFT calculation as shown in section 4.7. While for solid state systems GGA functionals are often suitable starting points (see also section 6), hybrid functionals perform better for molecules. A judicious choice of the starting point is necessary.

5. Convergence of technical parameters

G0W0 calculations usually require the convergence of a few additional parameters, which are strongly implementation dependent. Such a parameter is, e.g., the size of the integration grid for the imaginary frequency terms in the CD and AC approach. G0W0 practitioners should always carefully check their GW code to ensure the robustness and convergence of all available settings and parameters

The G0W0 approximation provides computationally efficient access to the whole QP spectrum. Despite these appealing features G0W0 has certain drawbacks. The most severe is the dependence on the starting point discussed in section 4.7. Furthermore, the ground state energy and density cannot be computed. In sections 5 and 10, we will show how these drawbacks can be overcome.

5. Beyond G0W0: Self-consistency Schemes

5.1. Fully Self-consistent GW

To go beyond G0W0, one must include some level of self-consistency in Hedin's GW equations. The conceptually purest approach to GW is to perform full self-consistency in the GW equations, denoted as scGW. It is also the most expensive. As introduced in section 3, all four quantities are iterated until self-consistency in the Green's function is achieved. Until now, self-consistent GW is the rarest form of GW for reasons of computational expense and conceptual controversy (see below), although that is slowly changing.

The first scGW calculation was performed for the homogeneous electron gas (HEG) by Holm and von Barth (1998), after the same authors had previously applied partial self-consistency (scGW0) (von Barth and Holm, 1996). Later studies were extended to the 2D HEG (García-González and Godby, 2001). scGW deteriorates the spectral properties of the HEG compared to G0W0. This deterioration manifests itself in a quasiparticle bandwidth that is larger than the free electron one and a broad and featureless satellite spectrum. Both results contradict experimental evidence for alkali metals which are HEG-like (von Barth and Holm, 1996). Also, band gaps of simple semiconductors are greatly overestimated by scGW (Schöne and Eguiluz, 1998; Grumet et al., 2018). scGW calculations for atoms (Stan et al., 2009) and molecules (Rostgaard et al., 2010; Caruso et al., 2012a, 2013a,b; Marom et al., 2012) show improvements over G0W0 for the first ionization energies and transport properties (Strange et al., 2011) of finite systems. With regard to the whole spectrum, however, scGW is usually outperformed by G0W0 with a judicious starting-point choice (Marom et al., 2012; Caruso et al., 2013a; Knight et al., 2016).

scGW is computationally more demanding than G0W0 because the full Green's function must be stored and calculated (Caruso et al., 2013a), increasing memory and computation requirements. In G0W0, the full Green's function is only required in O(N3) schemes (see section 4.8). Other implementations make use of the fact that intermediate quantities can be expressed in terms of the mean-field wave functions and eigenvalues, which reduces the computational complexity (see section 4.1). Furthermore, iterations of the GW equations for scGW are expensive. χ0, W, and Σ must be computed at each iteration, with a potentially high computational time for even a single evaluation of Σ.

Conceptually, the additional self-consistency in the Green's function adds more reducible diagrams compared to G0W0, as Figure 22 illustrates. In G0W0, the bare Coulomb interaction is screened by a series of sequentially interacting electron-hole pairs, or “bubbles.” In scGW, this structure is preserved, but the bubbles are now composed of interacting Green's function lines instead of non-interacting G0 lines. This effect is a general feature of iterating Green's function diagrams. By iterating diagrams for a given approximation, initial G0 lines at internal times become interacting G lines. Already after the first iteration of the cycle (G1 in Figure 22), the Green's function contains sequential self-energy insertions.


Figure 22. G0W0 and scGW in terms of Feynman diagrams. In G0W0, the irreducible self-energy is constructed from G0 and W0. The Green's function updated with the lowest order self-energy, G1 (shown as the bold Green's function line), contains an infinite series of Σ insertions. In fully self-consistent GW, the starting point dependence is removed and all quantities in the diagrammatic expansion are fully dressed. Here, we assume a true G0 starting point instead of a mean-field G0 so that subtraction of vMF is not necessary to include in the diagrams.

Let us look more closely at how this occurs. Recall from Equation (13) that Dyson's equation is

G(1,2)=G0(1,2)  +G0(1,3)Σ(3,4)G(4,2)d(3,4).    (81)

The first guess at the full G, labeled G1, would then be

G1(1,2)=G0(1,2)  +G0(1,3)Σ0(3,4)G1(4,2)d(3,4),    (82)

where we have inserted G1(4, 2) in place of G(4, 2) on the right-hand-side (RHS). Σ0 is the first estimate to the self-energy, evaluated with G0 wherever G lines enter the self-energy diagram. G1 appears on both sides of Equation (82) − by replacing G1(4, 2) on the RHS with the entire RHS, one can generate a reducible diagram for G1 that is O(Σ02). At this point, the series for G1 contains three parts: G0, a term of order O(Σ0), and a term of order O(Σ02). These contributions to G1 are shown in Figure 22. By further iterating G1 on the RHS, one can generate all reducible diagrams which contribute to G1. Despite the infinite number of reducible diagrams generated by this prescription, G1 is still computed only with the G0W0 self-energy because we have not updated Σ0. This example also demonstrates why it is conceptually much simpler to work only with the irreducible self-energy and avoid this infinite, reducible series. Indeed, iterating Equation (82) to find G1 while keeping Σ0 fixed is equivalent to generating the entire reducible series for G1.

In scGW, the G0W0 calculation of Σ0 to build G1 is only the first step. Next, we update Σ0 to Σ1 by inserting G1 into the self-energy diagram. The diagram contains one obvious G line (Σ = iGW), but contains more that are hidden in the polarizability entering W. By updating the polarizability with G1 in place of G0, the diagrams contained in G1 enter the screened Coulomb interaction. Just as before, we can generate all reducible diagrams contributing to the updated Green's function (G2) by iterating the Dyson series

G2(1,2)=G0(1,2)+G0(1,3)Σ1(3,4)G2(4,2)d(3,4)    (83)

for a fixed Σ1. Continuing to update Σ and iterate G introduces more and more reducible diagrams. The solution is self-consistent when G entering Σ is the same as G from iterating Dyson's equation.

In real scGW calculations, the procedure is slightly different. G and Σ are updated together instead of iterating to find Gi+1 at a fixed Σi. After the first iteration of Equation (82), the updated − but not yet self-consistent − G1 is inserted into Σ. This way, Σ is updated at each iteration along with G. The combined iterations are much more efficient because G and Σ converge together. Bear in mind that the efficient method of updating G and Σ at each step does not form the same easy-to-interpret series for G1 as in Figure 22. Note that even after one iteration to find G and Σ, we would already go beyond G0W0.

Based on Figure 22, the fully dressed Green's function and screened Coulomb interaction in scGW can be interpreted as double renormalizations of G0 and W0 through the two Dyson's equations in the GW equations. However, the third Dyson equation that the vertex function (Equation (B26) in Appendix B) would introduce is missing from the GW equations. The absence of the vertex function has important consequences. Figure 23 shows two of the reducible diagrams that enter W in scGW, but that are not present in W0. The diagram on the left shows the polarization bubble with the insertion of one interaction, or scattering event, in each arch. It is part of a series of sequential scattering events and builds additional interactions into the screened Coulomb interaction. The other diagram, however, is problematic. After the creation of the first electron-hole pair, both the electron and the hole interact with a new electron-hole pair. The two new electron-hole pairs are composed of the same Green's function lines as the initial electron-hole pair, even though the initial pair still exists. Therefore, the two later pairs do not account for the fact that the initial electron-hole pair has already been created − they should somehow omit the pair already created. The electron (or hole) thus interacts or correlates with itself. The problems that have been identified9 for scGW can be attributed to diagrams like the one on the RHS of Figure 23 (Romaniello et al., 2009).


Figure 23. Two of the diagrams of the screened Coulomb interaction in scGW that are not present in G0W0.

We would also like to briefly comment on partial self-consistency in the Green's function. The above discussion points to the screened Coulomb interaction as the major source of imbalance between self-consistency and the missing vertex corrections. This imbalance can be partially fixed by keeping W fixed at the W0 level and iterating only the Green's function to self-consistency (scGW0). This is shown schematically in Figure 24. This approximation is partially motivated by a “best G, best W” philosophy that can improve agreement with experiment in certain situations.


Figure 24. Schematic of Hedin's full set of equations (A) and Hedin's GW approximation (B–D). In (A), all five quantities are iterated to self-consistency. In (B), self-consistent GW (scGW), Γ is set to a single spacetime point and the remaining four quantities are determined self-consistently. Eigenvalue self-consistent GW shown in (C), evGW, updates only the quasiparticle energies while leaving the wave functions unchanged. In the scGW0 or evGW0 procedures shown in (D), one iterates G to self-consistency in Dyson's equation but does not update χ0 or W.

Before leaving the discussion of self-consistent GW, we introduce one of the most technical and modern aspects of Green's function theory being researched: the existence of multiple solutions for G from a single Dyson equation (Tandetzky et al., 2015). This issue has been studied in detail for the zero-dimensional one-point model (OPM) (Lani et al., 2012; Berger et al., 2014; Tarantino et al., 2017) and has been produced numerically (Kozik et al., 2015; Gunnarsson et al., 2017; Vučičević et al., 2018). In the analytic OPM, there exist two interacting G which can be mapped from the same G0 (Rossi and Werner, 2015; Stan et al., 2015). One of these solutions is the physical G for all values of interaction strength. Here, the physical solution is characterized by a smooth connection to G0, unlike the unphysical solution for G which diverges at zero interaction strength. The reverse map, from G to G0, has two solutions for G0 which must be disentangled at a certain interaction strength. At this point, the physical G0 switches between the two solutions, so that solving the problem for all interaction strengths requires switching solution methods at this point. Otherwise, one would obtain a physical G0 for some interaction strengths and an unphysical G0 for others. In the OPM, this is now understood. However, it is not well understood if or how this same phenomenon emerges in more realistic systems.

5.2. Eigenvalue Self-consistency and Level Alignment

There are a few strategies to reduce the expense of scGW while still including more physics than G0W0. The simplest form of performing approximate self-consistency in GW is to iterate in the eigenvalues (evGW). After completion of the G0W0 loop, the real parts of the quasiparticle energies obtained from Equations (22) or (34) are reinserted into the non-interacting Green's function G0 (Equation (24)) in place of the starting eigenvalues. Through G0, the change in eigenvalues permeates through W0 to the self-energy and eventually to the quasiparticle energies (evGW). After iterating until the input quantities equal the output, the equations are self-consistent in the eigenvalues.

Eigenvalue self-consistency was already proposed in the first G0W0 calculation for real materials (Hybertsen and Louie, 1986) and has since been applied frequently (see, e.g., Shishkin and Kresse, 2006b for a more in-depth analysis). However, since only the real part of the quasiparticle energies is used and the wave functions are not updated self-consistently, the starting point dependence cannot be eliminated entirely. For example, a study for azabenzenes demonstrates that although the starting point dependence is reduced from 1.4 eV in G0W0, it cannot be lowered beyond 0.4 eV (Marom et al., 2012).

For molecular systems, it has been shown that eigenvalue self-consistency improves the HOMO-LUMO gaps, which are then in good agreement with experiment (Blase et al., 2011; Wilhelm et al., 2016). However, examining the entire spectrum reveals that evGW does not lead to consistent improvements over G0W0. The evGW spectra are overly stretched with respect to the experimental spectra, such that large deviations (on the order of 1 eV) from experiment occur for lower lying states. Moreover, for most systems, the orbital ordering deviates from experimental observations (Marom et al., 2012). This is in line with observations for semiconductors and insulators, that find band gaps to be considerably overestimated in evGW (Shishkin and Kresse, 2006b). The reason for this overestimation in solids lies in the fact that the insertion of the quasiparticle energies into the screened Coulomb interaction leads to an underscreening, which should be compensated by the missing vertex corrections, as discussed previously.

Just as scGW0 ameliorates problems in self-consistent GW, one can perform eigenvalue self-consistency only in G to circumvent underscreening errors. Iterating the eigenvalues in only G produces an evGW0 scheme that gives band gaps in good agreement with experiment (Shishkin and Kresse, 2006b). However, for open shell systems, it was observed that eigenvalue self-consistency in G strongly affects the calculated multiplet splittings (Lischner et al., 2012).

The effect of eigenvalue self-consistency in G on the self-energy is demonstrated in Figure 25 for an evGW0 calculation. Compared to G0W0, the structure of the self-energy is almost identical, but shifted to lower energies. The Green's function in the eigenvalue self-consistent GW0 scheme is given by

Gevσ(r,r,ω)=mϕmσ0(r)ϕmσ0*(r)ωϵmσiηsgn(EFϵmσ)    (84)


ϵmσ=ϵmσ0+Δϵmσ    (85)
Δϵmσ=Σmσ(ϵmσ)vmσMF.    (86)

Inserting the GW corrections Δϵ into the Green's function results in a shift of the poles in the self-energy, see Equation (48). On the real axis, the poles of Σsc are located at ωiσn=ϵiσ0+Δϵiσ-Ωnσ and ωaσn=ϵaσ0+Δϵaσ+Ωnσ, where i indicates again occupied and a virtual states. Starting from a GGA functional, the correction Δϵ is negative for occupied and positive for virtual states. Compared to a G0W0 scheme, the poles ωiσn are now located at lower and the poles ωaσn at higher frequencies.


Figure 25. Self-energy matrix elements for the HOMO of a single water molecule obtained with G0W0, evGW0 and level-aligned G0W0. In all three cases PBE is used as starting point. The inlet shows the graphical solutions of the QP equation. See Appendix C for further computational details.

A simplified version of evGW0 was originally suggested by Hedin (Hedin, 1965) and has been revisited a few times by others (Pollehn et al., 1998; Hedin, 1999; Martin et al., 2016). Instead of using an individual shift Δϵ for each state m, a global shift ΔE is employed:

GΔEσ(r,r,ω)=          mϕmσ0(r)ϕmσ0*(r)ω(ϵmσ0+ΔEσ)iηsgn(EFϵmσ0),    (87)

where G0σ(ω-ΔEσ)=GΔEσ(ω). The QP equation (Equation (22)) then transforms into

ϵsσ=ϵsσ0+Σsσ(ϵsσΔEσ)vsσMF.    (88)

For metals, the shift ΔE is chosen in such a way that the G0W0 Fermi energy aligns with that of the starting point calculation, i.e., with the Fermi level of G0. For systems with a energy gap, the highest occupied state is aligned, i.e., the valence band maximum for solids or the HOMO for finite systems. The latter is motivated by “DFT Koopman's theorem,” which states that only the KS energy of the HOMO can be rigorously assigned to the ionization potential when starting from an exact DFT functional (Levy et al., 1984; Almbladh and von Barth, 1985). In that case ΔE would be zero.

The shift can be determined by demanding self-consistency for the highest occupied state

ϵHOMO,σ=ϵHOMO,σ0+ΔEσ.    (89)

Inserting Equation (89) into Equation (88) yields the explicit expression

ΔEσ=ΣHOMO,σ(ϵHOMO,σ0)vHOMO,σMF.    (90)

Adjusting the energy scale of G0σ by ΔE translates to a rigid shift of the self-energy as shown in Figure 25. The results are very similar to evGW0 in the frequency range where the quasiparticle solution is expected.

The ΔE scheme is less frequently used for the calculation of quasiparticle energies than eigenvalue self-consistent schemes. However, it has been shown that it substantially improves satellite spectra (Pollehn et al., 1998). The same holds for the evGW0 scheme, which has been employed to calculate satellite spectra of VO2 (Gatti et al., 2015) and bulk sodium (Zhou et al., 2015).

5.3. Self-consistency via a New Ground State

Building on the idea of iterating in the quasiparticle energies, one can go one step further and also incorporate wave function changes. An elegant way to achieve this is to find the variationally best mean-field potential to a given self-energy (Godby et al., 1986, 1987b; Casida, 1995a; van Schilfgaarde et al., 2006; Kotani et al., 2007b). This mean-field potential can then be used to generate new eigenvalues and eigenfunctions to construct a new G0 for the next iteration of the GW cycle.

If the new potential is local, this iteration can be formalized exactly in the optimized effective potential (OEP) framework (Casida, 1995a; Kümmel and Kronik, 2008), which is equivalent to the Sham-Schlüter equation (Godby et al., 1986, 1987b). The OEP framework and the Sham-Schlüter equation only guarantee that the density generated by the local potential matches the GW density. The eigenvalue spectrum of the local potential will not be the same as the GW spectrum. For explicitly non-local potentials, no formally exact match between the GW self-energy and the potential has been found because the self-energy is non-local and frequency dependent, while the constructed potential is non-local but static.

An approximate non-local potential can be found by introducing the GW Hamiltonian ĥGW(ω)=ĥ0+vH+ΣGW(ω). The mean-field Hamiltonian ĥMF that best reproduces the effects of ΣGW is defined as ĥMF=ĥ0+vH+vMF, see also Equations (18–20) for the definitions of the Hamiltonians. vMF can then be obtained by minimizing ||ĥGWĥMF|| (van Schilfgaarde et al., 2006; Kotani et al., 2007b). An approximate minimization finally yields an analytic expression for the (static and Hermitian) mean-field potential (Faleev et al., 2004; van Schilfgaarde et al., 2006; Kotani et al., 2007b)

vijMF=12[[ReΣ(ϵi)]ij+[ReΣ(ϵj)]ij],    (91)

where “Re” signifies here the Hermitian part of Σ(ϵk)

[ReΣ(ϵk)]ij=12[Σ(ϵk)+Σ(ϵk)]ij.    (92)

The quasiparticle energies for the Green's function G are then given by the self-consistent G0 that follows from vijMF (van Schilfgaarde et al., 2006). Satellites or the incoherent part of the spectral function are not captured by this approximation. This is why the scheme is commonly referred to as quasiparticle self-consistent GW (QSGW). Reports of a starting point dependence for metal oxides (Liao and Carter, 2011; Isseroff and Carter, 2012) have not yet been reproduced by other groups with a different implementation. In general, the QSGW scheme converges to a unique solution.

An alternative definition for a non-local mean-field potential is given by the static Coulomb hole plus screened exchange (COHSEX) approximation to GW (Hedin, 1965; Hedin and Lundqvist, 1970):

vσMF,COHSEX(r,r)=ΣσCOH(r,r)+ΣσSEX(r,r).    (93)

The screened exchange (SEX) term is defined in analogy to the exact-exchange self-energy in Equation (32) but with the statically screened Coulomb interaction instead of the bare one

ΣσSEX(r,r)=ioccϕiσ(r)ϕiσ*(r)W(r,r,ω=0),    (94)

where ϕ are eigenfunctions of the COHSEX mean-field Hamiltonian. The static Coulomb hole (COH) term, on the other hand, becomes local in space

ΣσCOH(r,r)=δ(rr)[W(r,r,ω=0)v(r,r)].    (95)

The statically screened Coulomb interaction W(r, r′, ω = 0), which enters in Equations (94) and (95), is obtained by inserting the COHSEX eigevalues and eigenfunctions in Equations (25–28) for ω = 0. Like in QSGW, vσMF,COHSEX(r,r) produces new eigenvalues and eigenfunctions, which yield a new self-energy. The COHSEX equation can then be iterated until self-consistency is achieved.

COHSEX can also serve as an improved starting point compared to KS-DFT for a perturbative G0W0 calculation. After completing a COHSEX calculation, one can use the self-consistent COHSEX eigenvalues and wave functions for a perturbative G0W0 calculation with the full, dynamical W. In the case of VO2 (Gatti et al., 2007), G0W0@LDA fails to open the band gap. On the other hand, G0W0@COHSEX opens a band gap, in agreement with experiment. The improved COHSEX starting point is especially important for materials with localized electrons (Aguilera et al., 2011). When comparing to benchmark coupled cluster data on organic molecules, the COHSEX starting point decreases the mean absolute error of G0W0 compared to G0W0@PBE (Knight et al., 2016). In Ge under pressure, G0W0@COHSEX predicts a direct gap at the Γ point while G0W0@LDA predicts band overlap (Jain et al., 2014).

Self-consistency in GW is a topic that is still being researched. While the results from true self-consistent GW are usually in worse agreement with experiment than G0W0 (Schöne and Eguiluz, 1998; Shishkin and Kresse, 2006b; Grumet et al., 2018), studies of self-consistency are still necessary to advance the field. We can learn about shortcomings of the theory or assess challenging materials. Most importantly, self-consistent GW implementations are a necessary foundation to go beyond GW in the future, as discussed in section 11.

6. Solids

Solids were the first testbed of GW theory in real materials. The basic quantity to compute in solids is the band structure. Unlike in molecules with single-particle states given by molecular orbitals, single particles in solids occupy Bloch waves, defined in Equation (50) and indexed by their wave vector k. Diagonalizing the Hamiltonian at each k gives its own set of single-particle eigenvalues. One can conveniently visualize the eigenvalues at different k-points by varying k continuously along some path, placed on the x-axis, and plotting the eigenvalues on the y-axis. Eigenvalues change continuously with k, forming separate bands of states. The collection of single-particle states in bands is similar to the grouping of different combinations of bonding and anti-bonding states in molecules or polymers. Quasiparticle Hamiltonians are also k-dependent, and energies at all k form a band structure which can be compared to angle resolved PES and IPES spectra as the incident momentum is varied (see Figure 2).

Computing quasiparticle band structures with GW gives the quasiparticle band gap, analogous to the HOMO-LUMO gap in molecules. The first G0W0 calculations for real materials (Strinati et al., 1980, 1982; Hybertsen and Louie, 1985, 1986) focused on semiconductors (Si, Ge) and insulators (diamond, LiCl). G0W0 calculations in semiconductors and insulators give a uniform improvement in band gap over estimates with either Hartree-Fock or Kohn-Sham eigenvalues, as shown in Figure 26A. This improvement was the first major success of the GW theory.


Figure 26. (A) Band gaps of semiconductors and insulators computed with PBE, G0W0, and evGW0 in the all-electron, linearized augmented plane wave (LAPW) framework. Data taken from Jiang and Blaha (2016). (B) Band gaps of semiconductors and insulators computed with PBE, G0W0, QSGW, and scGW in the projector-augmented-wave (PAW) framework. Data taken from Grumet et al. (2018).

The success of GW applied to semiconductors continued with other studies (Godby et al., 1986, 1987a,b, 1988; Blase et al., 1995). Many common semiconductors lack semicore states and are well described by pseudopotentials (see section 4.4.4). These factors reduce the computational complexity for GW calculations and made early, realistic GW calculations of semiconductors feasible. Screening in simple semiconductors can also be approximated by model dielectric functions like plasmon-pole models (see section 4.3.1), eliminating the need for a numerical evaluation of the self-energy integral.

It was later realized that semicore d-electrons in semiconductors such as GaN, ZnO, ZnS, ZnSe, or CdS (Rohlfing et al., 1995b, 1997a, 1998) and metals such as Cu (Marini et al., 2002) and Au (Rangel et al., 2012) have a strong influence on GW calculations. Due to the strong overlap of the atomic d functions with the atomic s and p functions in the same shell, the exchange self-energy is very sensitive to the inclusion of semicore states. If only the semicore d states are explicitly included as valence state in the G0W0 calculation, while the s and p states in the same shell are frozen in the core of a pseudopotential, the subsequent G0W0 calculation will produce an incorrect band gap (Rohlfing et al., 1995b). This problem can be solved by explicitly including the entire shell as valence in the G0W0 calculation (Rohlfing et al., 1995b; Luo et al., 2002; Tiago et al., 2003a; Fleszar and Hanke, 2005), by using exact-exchange pseudopotentials and exact-exchange starting points (Rinke et al., 2005, 2008a; Qteish et al., 2006) or by all-electron calculations (Friedrich et al., 2006, 2010; Shishkin and Kresse, 2006a; Gulans et al., 2014, Jiang and Blaha, 2016).

6.1. Band Gaps

Figure 26A shows the quasiparticle band gap computed with G0W0 and evGW0 with a modern all-electron LAPW code (Jiang and Blaha, 2016) for several different semiconductors and insulators. Perfect agreement between theory and experiment would place all data points on the dashed blue line. Generally, Kohn-Sham eigenvalues based on a multiplicative (local or semi-local) exchange-correlation potential (here PBE) underestimate the band gap and Hartree-Fock eigenvalues overestimate (not shown in Figure 26). Eigenvalue self-consistency (evGW0) improves the agreement with experiment even further than G0W0, when starting from a local or semi-local DFT calculation.

Figure 26B compares band gaps computed with different self-consistency schemes (Grumet et al., 2018) for a different set of semiconductors and insulators than in panel (A). G0W0@PBE again provides good agreement with experiment (i.e., the red squares are close to the diagonal). Van Schilfgaarde's QSGW scheme (van Schilfgaarde et al., 2006) and fully self-consistent GW calculations consistently overestimate band gaps.

With the predictive accuracy of G0W0 band gaps validated, we provide a few examples in which GW calculations helped to resolve band gap controversies. One case is InN. In the early 2000s, alloys of GaN and InN were revolutionizing light-emitting diode (LED) technology. However, the band gap of InN was believed to be almost 2 eV (Butcher and Tansley, 2005), which would have severely limited the usefulness of InGaN alloys to tune the emission of LEDs. Through G0W0 calculations and more refined experiments, the real value of the InN band gap was found to be 0.7 eV (Bechstedt and Furthmüller, 2002; Furthmüller et al., 2005; Rinke et al., 2006), paving the way for the LEDs we know today. Another example is hybrid perovskites that have triggered a new boom in the emergent photovoltaic materials field. The prototypical material is methylammonium lead triiodide (CH3NH3PbI3 or in short MAPI). Unusually, local or semi-local DFT calculations already predict a band gap in good agreement with measurements, which had caused initial confusion in the field. However, when spin-orbit effects, which are particularly strong in this materials class, are incorporated in the DFT calculations, the band gap becomes significantly underestimated again. G0W0 and QSGW calculations that include spin-orbit effects then predict the correct band gap (Brivio et al., 2014; Umari et al., 2014).

Finally, we consider high pressure physics. At high pressures (~ 100 GPa), many materials experience band gap closure and transition from an insulator to a metal. There can also be many competing structural phases, each with their own metallization pressure, that are difficult to disentangle in experiments. GW is an excellent tool to theoretically predict the metallization pressure for different structural phases and help interpret experimental results (Khairallah and Militzer, 2008; Tse et al., 2008; Ramzan et al., 2010; Jin et al., 2016; Yang, 2017). Solid hydrogen is a noteworthy example of metallization at high pressure, first predicted in 1935 (Wigner and Huntington, 1935). The metallic hydrogen puzzle is an exceptionally difficult one that is still not fully understood, but GW calculations help corroborate experimental measurements and support the existence of certain structural phases (Lebègue et al., 2012; Dvorak et al., 2014; McMinis et al., 2015).

6.2. Band Structures and Band Parameters

With GW, one can compute much more than only the band gap of a solid. A typical band structure computed with GW is shown in Figure 27 for ZnO. To visualize the results, k is allowed to vary on a linear path through the Brillouin zone. The quasiparticle energy as a function of k is also called the dispersion for the system. The G0W0 band structure for ZnO (red lines in Figure 27) is superimposed on the experimental photoemission results shown already in Figure 2. Experiment and G0W0 agree very well both in terms of band positions as well as band curvatures.


Figure 27. G0W0 band structure of ZnO superimposed on experimental ARPES data (Yan et al., 2012). The experimentally measured lifetimes of the states are indicated by the shading, with white shading indicating long lifetime. The G0W0 calculations are based on the optimized effective potential approach for exact exchange mentioned in section 5.3 that includes LDA correlation [OEPx(cLDA)]. Reprinted with permission from Yan et al. (2011). Copyright (2011) by IOP Publishing Ltd.

Another example of a G0W0 band structure is shown in Figure 28 for K2Sn3O7, a new prospective ion conductor or transparent conductor (McAuliffe et al., 2017). The unoccupied states in the PBE band structure have been shifted up for the purposes of plotting so that the bottom of the conduction bands coincides in PBE and G0W0@PBE. This removes the PBE band gap problem from the comparison and makes it easier to spot differences in band curvatures. For the valence bands the PBE and G0W0@PBE band structures agree remarkably well for this material. Toward lower energies the deviations between the band structures become larger with G0W0@PBE generally giving lower band energies than PBE. This downward shift leads to a band width widening in G0W0 compared to PBE. For the conduction bands, the difference between PBE and G0W0@PBE is more pronounced. The band curvatures in G0W0@PBE are much steeper than in PBE, which subsequently leads to a significant underestimation of the PBE bands around the X, S, U, and R points in the Brillouin zone. K2Sn3O7 is another example of a material whose band gap and band structure were not known. The G0W0@PBE band gap amounts to 3.15 eV (McAuliffe et al., 2017), which now provides a reference value for this new material.


Figure 28. G0W0 band structure of K2Sn3O7 (McAuliffe et al., 2017). The main panel illustrates the difference between the PBE (dark gray lines) and the G0W0@PBE (red lines) band structure. The unoccupied states of the PBE band structure have been shifted up in energy for better visibility so that the bottom of the conduction bands coincide in both band structures. The right panel shows the G0W0@PBE density of states (DOS) resolved into s, p and d angular momentum channels.

From the band structure, one can access the band gap, band widths, and estimate effective masses. If one models the dispersion at the band edges as parabolic, as is the case for a free particle, one can extract an effective mass from the band structure. The effective mass, labeled m*, is

m*=2[d2Edk2]1    (96)

so that the band edge dispersion is

E=2k22m*    (97)

to mimic a free particle. In a real crystal, the effective mass is a tensor, not just a scalar. The effective mass model is closely related to the quasiparticle concept, and the renormalization factor Zs (Equation (35)) is one factor contributing to m*. The quasiparticle effective mass is usually heavier than the free electron mass because of the drag induced by the surrounding electrons. Static mean-field theories like Kohn-Sham DFT also give an estimate of effective mass from their band structure. The GW band structure is typically more “curved” or concave than the Kohn-Sham structure, as shown in Figure 28, which means that GW quasiparticles are “lighter” than the KS particles. In silicon and methylammonium lead iodide, the GW level of theory is necessary to predict effective masses in good agreement with experiment (Filip et al., 2015; Poncé et al., 2018).

One can also compute the single-particle density of states (DOS), which in solids is analogous to taking horizontal slices through the band structure. This gives the total number of available states at the energy of that slice. In Green's function theory, the concept of the single-particle DOS is replaced by the spectral function. The spectral function is k-dependent, so that the total effective DOS is obtained by adding up the spectral functions at all k. However, computing the spectral function requires the solution of Dyson's equation, which is often not practical computationally. Instead, a G0W0 quasiparticle DOS is computed by summing up artifically broadened Gaussian peaks centered around each G0W0 energy. The right panel of Figure 28 shows the G0W0@PBE DOS for K2Sn3O7 (McAuliffe et al., 2017). In addition, this DOS is projected on the atomic angular momentum channels s, p, and d. Such information is usually extracted from DFT calculations and illustrates that both the valence band and the conduction band of K2Sn3O7 is largely made up of p states.

Band parameters like effective masses are important characteristics of semiconductors and are key parameters for the semiconductor industry. G0W0 effective masses are more accurate than those computed with DFT. Effective masses are either extracted directly from the GW band structure by fitting Equation (96) to a fine band structure path (Schleife et al., 2009) or by fitting an effective k · p Hamiltonian to GW quasiparticle energies (Rinke et al., 2008b). In this way, important band parameters have been computed for silicon and silicon under strain (Bouhassoune and Schindlmayr, 2010; Poncé et al., 2018), GaAs (Cheiwchanchamnangij and Lambrecht, 2011), AlN, GaN, and InN (Rinke et al., 2006, 2008b; Svane et al., 2010; Yan et al., 2011), MgO, ZnO, and CdO (Schleife et al., 2009; Yan et al., 2012) and more recently for perovskites and hybrid perovskites (Filip et al., 2015). Such band parameters can then be used directly in device simulations to model, for example, charge carrier flows (Kivisaari et al., 2017). If one is interested in charge carrier mobilities and charge carrier densities, scattering due to phonons and impurities has to be taken into account (Kioupakis et al., 2010; Poncé et al., 2018).

GW can further be used as one of the final steps in high-throughput screening studies for new materials. In a search for transparent p-type conductors, G0W0 calculations provided accurate band gaps and effective masses that screened out the final candidates (Hautier et al., 2013).

6.3. Lifetimes

Unlike mean-field theories, GW also allows one to compute the lifetimes of states from first principles. The lifetime of the quasiparticle is the characteristic time over which the added particle decays into surrounding degrees of freedom. States which are “closer” to exact eigenstates of the system have longer lifetime. The lifetime of a quasiparticle with corresponding energy ϵs is directly related to the non-Hermiticity of the self-energy and the magnitude of its imaginary part,

τs1=2|ImΣ(ϵs)|.    (98)

τ can be inferred from experimental spectra by its relation to the quasiparticle peak width, Γ, as τ−1 = Γ/2 (not to be confused with the vertex function Γ).

From simple arguments in Fermi liquid theory, lifetimes decrease as particle energy increases because the available phase space for scattering at a fixed energy grows with increasing energy. Studies of quasiparticle lifetimes in the GW approximation for metals (Cu, Ag, Au) show good agreement with experiment at high energies (Keyling et al., 2000; Bacelar et al., 2002; Marini et al., 2002; Yi et al., 2010). At low energies, however, the agreement is not perfect. For example, GW calculations cannot account for the sudden increase in experimental lifetimes of electrons in Cu at energies below 2 eV (Keyling et al., 2000; Yi et al., 2010). These failures are attributed to the localized, short-range interactions of d-electrons in metals that are not described well by GW.

6.4. More Challenging Solids

As computational power increased and the success of GW became more widely known, studies were extended to more challenging materials like oxides, or d- and f-electron compounds. These materials are both a theoretical challenge for the GW approximation and are numerically more difficult to compute. Broadly speaking, these materials suffer from a severe mean-field starting point problem and/or contain localized electrons which are not well described by GW. Accordingly, studies of these materials required advances in the treatment of core electrons and the evaluation of the self-energy. The first studies of metals focused on the alkali metals Na and Al (Northrup et al., 1987; Surh et al., 1988). Metals served as a valuable test on the effects of self-consistency and vertex corrections (Mahan and Sernelius, 1989; Shirley, 1996). Eventually, studies moved into oxides and materials with d-electrons (Aryasetiawan, 1992; Aryasetiawan and Gunnarsson, 1995; Massidda et al., 1995, 1997, Aryasetiawan and Karlsson, 1996).

Already in the early nineties of the previous century Aryasetiawan tackled ferromagnetic nickel (Ni) with G0W0 (Aryasetiawan, 1992). He found the quasiparticle band structure and the valence bandwidth to be in good agreement with angle-resolved photoemission data. However, the exchange splittings are not well reproduced by G0W0 and a satellite at 6 eV is missing. Later calculations for gadolinium (Gd) revealed similar observations (Aryasetiawan and Karlsson, 1996; Aryasetiawan, 1997). For Gd, satellites were seen in the G0W0 spectrum, but their spectral weight does not match experiment.

The previous millennium concluded with early explorations into transition metal oxides such as nickel oxide (NiO) and manganese oxide (MnO) (Massidda et al., 1995, 1997; Aryasetiawan and Karlsson, 1996). They, as well as iron and cobalt oxide (FeO and CoO, respectively), were then revisited with GW in the 2000s (Li et al., 2005; Kobayashi et al., 2008; Rödl et al., 2008; Rödl et al., 2009; Jiang et al., 2010b). These oxides present a challenge to G0W0 calculations because local and semi-local DFT starting points produce metallic states that then cannot be corrected into semiconductors by G0W0. Instead, DFT+U and hybrid functionals were explored as alternative starting points (Rödl et al., 2009; Jiang et al., 2010b). The resulting G0W0 DOSs are in good agreement with direct and inverse photoemission measurements for the low temperature magnetically ordered phases. However, the GW method cannot describe the DOS in the paramagnetic phase nor the transition to the paramagnetic phase.

The situation is similar in the less correlated copper oxide (Cu2O) (Bruneval et al., 2006). G0W0@LDA again fails to give a proper account of the band structure, while QSGW provides good agreement with ARPES measurements. CuO poses more of a problem, as no starting point or self-consistency scheme produces a satisfying band gap or density of states (Rödl et al., 2015, 2017).

The early 2000s saw other oxides gain rapid interest, as the semiconductor industry sought a replacement for silicon dioxide (SiO2) in silicon-based microelectronic technology. To prevent gate leakage in ever-shrinking transistors, gate materials with a higher dielectric constant (k) than SiO2 were required. Eventually hafnium dioxide won the race. During the development period, the electronic structure, in particular the band gap and the band offsets of so called high-k materials were of enormous interest (Shaltaf et al., 2008; Grüning et al., 2010; Jiang et al., 2010a; Sklénard et al., 2018). G0W0 calculations of the closely related compounds zirconium oxide (ZrO2) and hafnium oxide (HfO2) were performed (Grüning et al., 2010; Jiang et al., 2010a). Plane-wave and FLAPW G0W0 agree very well with each other for these materials. The all-electron calculations investigated the effect of the Hf f-electrons and found that they do not change the self-energy corrections in these materials (Jiang et al., 2010a). The final band gap of monoclinic HfO2, however, is still under debate. It was initially believed to lie around 5.8 eV and is now thought to be in excess of 6.3 eV (Sklénard et al., 2018). What remains a challenge in strongly polarizable materials such as high-k dielectrics, and could thus potentially explain remaining discrepancies between GW and experiment, is how to include ionic screening (i.e., screening due to nuclear motion) consistently in the dielectric function of a GW calculation.

The list of interesting metal oxides and metallic, semiconducting, or insulating solids is long and the number of GW calculations is steadily growing. Recent flagship applications even include defects, surface effects and solvents in their comparison to experiment (Gerosa et al., 2018b). For a recent review on the performance of different GW variants to metal oxides, we refer to (Bruneval and Gatti, 2014; Gerosa et al., 2018a).

While the f-electrons are relatively inert in HfO2, they assume a much more prominent role in lanthanide and actinide metals and oxides. With the exception of early explorations into Gd, GW calculations for f-electron compounds have only emerged fairly recently (Chantis et al., 2007; Jiang et al., 2009, 2012; Richter et al., 2011; Kutepov et al., 2012; Sakuma et al., 2012; Jiang, 2018). These calculations are almost always performed with DFT+U starting points or some form of self-consistency, as local or semi-local DFT provides a poor description of the electronic structure.

QSGW calculations for the rare-earth metals Gd and Er and the rare-earth monopnictides GdN, EuN, YbN, GdAs, and ErAs place the occupied 4f states in agreement with photoemission measurements, but then overestimate the position of the unoccupied f states (Chantis et al., 2007). Also, upon closer inspection, multiplet splittings are not reproduced with GW and require a beyond GW treatment (Richter et al., 2011). For the lanthanide sesquioxide (Ln2O3) series, G0W0@LDA+U calculations reproduce the relative positions of the occupied and unoccupied lanthanide f states across the series and confirm the experimental conjecture derived from phenomenological arguments (Jiang et al., 2009; Jiang, 2018).

Cerium (Ce) is another paradigmatic material. With only one f electron per Ce atom, it should still be relatively easy to describe, but Ce turns out to be an intricate material full of surprises. The phase diagram exhibits an unusual iso-structural phase transition. Both the α and the γ phase have an fcc crystal structure, but the α phase has a smaller equilibrium volume (Amadon et al., 2006; Bieder and Amadon, 2014; Devaux et al., 2015). The different localization of the f-electrons in the two phases is believed to be the driving force for the phase transition (Casadei et al., 2012, 2016; Devaux et al., 2015). Ce is traditionally thought to be a strongly correlated material that belongs to the realm of dynamical mean-field theory (DMFT), see section 11 for further details. However, the α and γ phases are already captured by the random-phase approximation (see section 10 for details). The G0W0 spectral function of the α phase is in good agreement with photo and inverse photoemission spectra (Sakuma et al., 2012). However, in the more correlated γ phase, G0W0 produces a peak at the Fermi level that is absent in the experimental spectra (Sakuma et al., 2012).

In conclusion of this section, we would like to reiterate that materials with d- or f-electrons remain one of the most challenging applications of GW. As a matter of principle, GW cannot yield an insulator with an odd number of electrons per unit cell. Such Mott insulators (Mott, 1968) are a manifestation of strong electronic correlation. Modern approaches to describing such strongly-correlated, localized states often combine GW with either a phenomenological or first-principles treatment of d- or f-electron correlation, a topic we discuss further in section 11.

6.5. Defects in Solids

So far, we have primarily discussed the performance of GW for computing the band gap in solids, which does not depend on the absolute values of the band edges. However, the locations of the valence band maximum (VBM) and conduction band minimum (CBM) are essential for understanding defect level alignment in solids. The conceptual problem is illustrated in Figure 29. Assume an initial LDA calculation and then a G0W0 calculation of the band structure for a system with a defect level in the gap. We assume that the G0W0 band gap is in good agreement with experiment, but what about the defect level? The position of the defect level relative to the band edges and Fermi energy is critical for determining its occupancy when the system is put in contact with an electron reservoir. The band gap alone is no longer enough to assess the accuracy of the calculation. Similar to the defect problem, band alignment at semiconductor heterojunctions depends on the absolute position of the levels, a problem which is discussed more in section 7.


Figure 29. For systems with mid-gap defect levels, computing the band gap alone is not enough to test the material for potential applications. The position of the defect level can also be computed with GW. In these cases, the absolute position of the VBM, CBM, and defect level are important.

The accuracy of GW for defect levels is still under investigation (Hedström et al., 2002, 2006; Weber et al., 2007; Ma and Rohlfing, 2008; Bruneval, 2009; Rinke et al., 2009; Bockstedte et al., 2010; Ma et al., 2010a) and its comparison to hybrid functionals is summarized in the review of Chen and Pasquarello (2015b). With a suitable choice of reference values to align the calculation with experiment, the accuracy of G0W0 is similar to that of hybrid functionals for predicting defect energy levels (Chen and Pasquarello, 2015b). Their major difference in performance can be attributed to their shift in the VBM, which has a direct effect on the defect level alignment and the calculated ionization potential. Hybrid functionals tend to symmetrically shift the VBM and CBM, while G0W0 mostly shifts the VBM down in energy which can worsen agreement with experiment for ionization potentials (Chen and Pasquarello, 2015b).

6.6. Outlook on Solids

As large-scale GW implementations became more common and parallelism increased, GW calculations became an indispensable tool for ab-initio predictions in solids. Today, there are too many GW calculations for solids to count. Even so, comparing GW calculations to experiment must be done with some care because there are additional effects in the experiment that are not included in ordinary GW. For example, electron-phonon coupling can have a significant effect on the band gap in real materials (Giustino et al., 2010b; Cannuccia and Marini, 2011; Botti and Marques, 2013; Antonius et al., 2014; Kawai et al., 2014). The effect of the electron-phonon interaction can also be described by a self-energy and calculated with perturbation theory (Cederbaum and Domcke, 1974; Smondyrev, 1986). Experimental spectroscopies are also surface sensitive techniques, as mentioned briefly in section 2.3, which means that the measured band structure may not correspond perfectly to the bulk states. These considerations aside, the impressive success of GW in solids encouraged studies of other systems, including surfaces and molecules.

7. Surfaces

The application of GW to surfaces and interfaces is not as common because these systems tend to have large unit cells with a number of atoms beyond the tractability of many GW codes. However, what makes surfaces particularly interesting from the GW perspective is a long-range polarization effect contained in the screened Coulomb interaction that is absent for bulk materials: the image effect. As illustrated schematically in Figure 30, an additional charge (hole created in the photoemission process or added electron in inverse photoemission) outside a surface induces an image charge in the surface (Deisz et al., 1993). This charge gives rise to an additional potential, the image potential, that renormalizes the energy of the electron or hole state. For metallic and dielectric surfaces it is easy to show from simple electrostatic considerations that the image potential should decay with the inverse distance from the surface. For other geometries, e.g., quantum dots or nanostructures, this decay behavior is modified (Rinke et al., 2004).


Figure 30. Illustration of the image effect. (B) Shows the image charge and image potential induced by an additional electron (e.g., anionic charge on a molecule) outside a surface. (A) Provides a graphic illustration how the image potential of a germanium (Ge) surface could be probed with a carbon monoxide (CO) test molecule. By adding thicker and thicker sodium chloride (NaCl) layers between CO and Ge, the CO molecule moves along the Ge image potential. The resulting CO gap will then depend on the NaCl layer thickness, which is indeed the case as (C) illustrates. Subfigure (A) adapted from Freysoldt et al. (2009) under the terms of the Creative Commons Attribution 3.0 License. Data for (C) obtained from Freysoldt et al. (2009).

That the GW self-energy encompasses the image effect was first shown by extracting the image potential from the GW self-energy for the Al(111) surface (White et al., 1997). Later, image resonances (Fratesi et al., 2003) and image states for semiconductors, insulators (Rohlfing et al., 2003; Kutschera et al., 2007) and nanoclusters (Rinke et al., 2004) were calculated with GW. Freysoldt et al. (2009) showed that the image potential can also be probed by monitoring the excitation energies of a test molecule (see Figure 30). The test molecule (CO) can be moved along the image potential by introducing insulating spacer layers between molecule and surface. The energy of the CO states gets renormalized stronger the closer it is to the surface, i.e., the smaller the spacer layer is. Freysoldt et al. (2009) also showed that the energy of semi-core states in different NaCl layers is affected by the image potential in the same way, a result that was later corroborated by Strange and Thygesen (2012) in a model study.

The most significant effect of the image potential is that it renormalizes the energy of adsorbates such as organic molecules (Freysoldt et al., 2009; García-Lastra et al., 2009; Thygesen and Rubio, 2009; Puschnig et al., 2012). The energetic position of molecular states near or on the surface is different from the molecule in the gas phase. During the excitation process, an electron or hole is added at the molecule. The additional correlation energy due to the polarization of the surface further stabilizes the added charge. As result, occupied states move up in energy and unoccupied states down and the HOMO-LUMO gaps reduce consequently in size, see Figure 30C. The renormalization depends on the dielectric constant of the surface. The larger the dielectric constant, the larger the renormalization. Already for surfaces of insulators the HOMO-LUMO gap renormalization is of the order of 1 eV and can reach more than 3 eV for metallic surfaces (García-Lastra et al., 2009; Thygesen and Rubio, 2009).

Apart from the HOMO-LUMO gap, the position of adsorbate states relative to the substrate's Fermi level or relative to the band edges is of significant interest in surface and interface science. This relative positioning of adsorbate to substrate states is commonly referred to as level alignment. GW calculations are currently considered to be the holy grail for an accurate determination of the level alignment. However, due to the aforementioned computational reasons (i.e., very large supercells) most GW level alignment calculations reported in the literature are not converged. Careful GW cluster calculations (Patrick and Giustino, 2012; Wippermann et al., 2014; Govoni and Galli, 2015) and very large scale GW calculations report good agreement with experiment. For physisorbed molecules, whose electronic states do not couple strongly to the substrate, the GW self-energy can be split into a surface and a molecular contribution. Such a simplified GW polarization model has been developed by Neaton et al. (2006) for weakly interacting molecules at metallic surfaces. This model has been used to compute GW estimates for the level alignment of amine-gold junctions and interfaces (Quek et al., 2007; Tamblyn et al., 2011) as well as π-conjugated polymers at Au(111) (Giovannantonio et al., 2018).

GW calculations for surfaces and interfaces are not only challenging because of the large supercells. An additional complication is the vacuum spacing in the common repeated slab model. In GW calculations that apply periodic boundary conditions, the surface is modeled as a slab of material that is periodic in two dimensions and finite in the third. The rest of the supercell in the direction perpendicular to the surface is filled with vacuum. Since the periodic boundary conditions apply also in the dimension perpendicular to the surface, the final system is a heterostructure of repeated blocks of material and vacuum (see Figure 31). GW now couples these repeated slab images because the GW interaction is long-ranged. The image potential decays only with the inverse distance between the slabs (see image effect discussion above) and not exponentially fast, as local or semi-local DFT functionals do. As a result, image potential tails generated by the repeated slab images reach into the surface region we would like to model with the slab model (see Figure 31) and obscure the actual image potential.


Figure 31. The image potential of a repeated slab system (B) differs from that of an isolated surface (A). The dashed lines in (C) mark the difference that can be computed with a suitable correction scheme (Freysoldt et al., 2008). As the charge moves across the interface, the ratio of dielectric constants for the “charged” and “uncharged” regions changes. As a result, the image potential changes sign.

In a GW calculation, the image potential is always present, even if we are not explicitly interested in image states. Due to the long range of the interaction, the vacuum spacing cannot be converged out in any GW implementation that has to place basis functions in the vacuum region (as for example plane waves) (Freysoldt et al., 2008; Hüser et al., 2013a). Two prevalent solutions to this problem have emerged: (1) to cut the interaction range and use an effective short-range interaction or (2) to apply post-processing corrections. The easiest way to limit the range is to impose a spherical cutoff on the Coulomb interaction every time it is used in the GW equations (Onida et al., 1995a; Spataru et al., 2004b; Ismail-Beigi, 2006; Rozzi et al., 2006). The largest disadvantage of this approach is that the spherical cutoff also limits the range of the GW interaction inside the material and in the two directions parallel to the surface. The cutoff radius should therefore at least be as large as the slab is thick. This implies that the vacuum separation should at least be equal to the slab thickness, which increases the computation time again for thicker slabs.

A computationally more efficient way is to apply post-processing corrections to a normal GW calculation that does not modify the range of the Coulomb interaction (Freysoldt et al., 2008). Care has to be taken, however, that the GW implementation correctly includes the dielectric tensor (Freysoldt et al., 2007). Otherwise, the GW calculation will not converge with respect to k-points (Freysoldt et al., 2007; Hüser et al., 2013a). Such a post-processing correction has been derived from an electrostatic model (Freysoldt et al., 2008) and is depicted in Figure 31. The true image potential is shown for two scenarios in Figure 31A: for a charge located outside or inside the slab. As the charge moves from outside the slab to inside, the image potential changes sign, as now the dielectric constant in the region where the charge resides (i.e. in the slab) is larger than where the charge is not (i.e. in the vacuum). Figure 31B shows the image potential for a periodic array of slabs in the repeated slab approach. It is notably different from the image potential of a single slab in Figure 31A. The correction derived by Freysoldt et al. (2008) is shown as black dashed lines in Figure 31C and restores the correct behavior for a single slab. The corrections can be several tenth of eV large and yield converged results already for small vacuum thicknesses (Freysoldt et al., 2008).

At surfaces, the DFT wave functions are sometimes poor approximations of certain surface states and image states. In such cases, it is desirable to calculate quasiparticle wave functions. This can be done by solving the full quasiparticle equation (Equation (21)) in a suitable basis. If this solution is performed iteratively in energy, new quasiparticle states, such as image states, can be found that are absent from the DFT spectrum. Examples where the quasiparticle wave functions differ notably from the LDA or PBE wave functions are GaAs(110) (Pulci et al., 1999) and the C(111) surface (Marsili et al., 2005) as well as image states (White et al., 1997; Rohlfing et al., 2003; Kutschera et al., 2007).

Early GW calculations for surfaces focused on surface states of simple semiconductors such as silicon (Rohlfing et al., 1995a, 1997b; Rohlfing and Louie, 1999; Hahn et al., 2001; Weinelt et al., 2004), germanium (Rohlfing et al., 1996), silicon carbide (Rohlfing, 2001), gallium phosphite (Schmidt et al., 1999), indium phosphite (Schmidt et al., 2000; Hedström et al., 2006) and insulators such as diamond (Marsili et al., 2005), lithium fluoride (Wang et al., 2003) and sodium chloride (Freysoldt et al., 2009). Frequently, the GW quasiparticle energies are taken as input for optical absorption or reflectance anisotropy spectroscopy (RAS) studies (Pulci et al., 1998; Schmidt et al., 2000; Hahn et al., 2001). The surface band structure and dispersion of surface states is in good agreement with available photoemission studies. Also, computed optical and RAS spectra agree well with experimental spectra for these systems. Later calculations for more complex surfaces or surface adsorbates have to be taken with a grain of salt, since they may not be fully converged with respect to all computational parameters, unless plasmon pole models, other model dielectric functions, or cluster models were used (Giorgi et al., 2011; Patrick and Giustino, 2012; Alves-Santos et al., 2014).

8. Two-Dimensional Materials

Research in two-dimensional materials developed rapidly after the isolation of graphene in 2004 (Novoselov et al., 2004). The crystal structure and Brillouin zone of graphene are shown in Figure 32. Two-dimensional materials have gained great fame for their interesting electronic structures, which include phenomena like Dirac fermions and topological insulators (Geim and Novoselov, 2007; Neto et al., 2009; Bhimanapati et al., 2015). Models of these effects are largely in the single particle − or single quasiparticle − picture. GW serves an important purpose to parameterize such models from a fully ab-initio perspective.


Figure 32. (A) Graphene has two hexagonal sublattices (A and B) in its honeycomb structure with translation vectors a1 and a2. (B) The Brillouin zone is hexagonal with two symmetry inequivalent corners labeled K and K′. (C) Near the Dirac points at K and K′, the dispersion is linear. The band structure is computed at the PBE level and taken from the Computational 2D Materials Database (Haastrup et al., 2018) with Fermi energy set to zero.

Two-dimensional materials often exist at the size and interaction strength that is ideally suited for GW. They are too large (the required Brillouin zone sampling is too dense) for more expensive wave function or beyond-GW Green's function methods, but their correlation is usually weak enough that GW gives a good description of their electronic structure. Similar to GW calculations on surfaces (section 7) or molecules (section 9), two-dimensional materials can show enhanced interaction effects from reduced dimensionality and decreased screening compared to bulk solids. Technical aspects of GW calculations of two-dimensional materials include the truncation of the screened Coulomb interaction between layers (similar to surfaces) and slow convergence with respect to k-points (Qiu et al., 2016; Rasmussen et al., 2016; Thygesen, 2017).

The band structure of graphene (and many other two-dimensional materials) is characterized by a zero band gap and linear dispersion near the Fermi energy, E(q)± ≈ ±vF|q| where the + (−) sign refers to electrons (holes) and q is the wave vector relative to the K or K′ points of the Brillouin zone, see Figure 32. vF is called the Fermi velocity and is the slope of the dispersion at the band edges. This linear dispersion is strikingly different than the parabolic dispersion in Equation (97), which is the case for most materials. Not long after its discovery, GW was applied to graphene to calculate the band structure and vF from first principles (Trevisanutto et al., 2008; Park et al., 2009; Siegel et al., 2011). Compared with calculations based on the local density approximation, GW preserves band closure at the Fermi energy and increases the Fermi velocity by ~17% to give a value of 1.1 × 106 m/s which is in good agreement with experiment. These studies also found kinks which appear in the low energy band structure from electron-phonon coupling and doping level dependent kinks of purely electronic origin.

Replacing carbon with a different group IV element creates a family of graphene-like materials. By preserving the honeycomb lattice of graphene, the materials still host Dirac fermions, but their chemistry and Fermi velocities depend on the specific element. For example, GW calculations of the Fermi velocity of planar silicon, called silicene, give a value of ~ 7.7 × 105 m/s (Huang et al., 2013). Because of silicon's tendency for sp3 hybridization, silicene also has a buckled structure which preserves linear dispersion at the GW level of theory (Wei et al., 2013). As with graphene, the electronic structures of silicene and germanene (monolayer Ge) subject to hydrogenation, strain, and hybridization with other materials have been studied with GW (Wei and Jacob, 2013b; Drissi and Ramadan, 2015a,b; Wang et al., 2015; Yan et al., 2015; Wang and Wu, 2017).

As one goes down the group IV elements, they become heavier; this has great significance for spin-orbit coupling (SOC) in two-dimensional materials. A SOC induced band gap in two-dimensional materials is critical to the topological character of their electronic structure. Stanene is a group IV monolayer (Sn) that has a sizable band gap due to SOC (Lü et al., 2012; Lu et al., 2017).

A number of functionalizations or structural modifications to graphene have been proposed for modifying its electronic structure. Much of the research in functionalized graphene is directed toward achieving semiconducting graphene or, more generally, two-dimensional semiconductors. As interesting as Dirac fermions are, semiconducting layers are necessary to build many layered electronic devices like field-effect transistors. For example, passivating graphene with hydrogen transforms it from a sheet of sp2 bonded carbon to sp3. The passivated structure, called graphane, has a GW band gap of ~ 5 eV (Lebègue et al., 2009; Leenaerts et al., 2010; Karlický and Otyepka, 2013; Hadipour and Jafari, 2015) and could be useful as a two-dimensional semiconductor. Other passivated graphenes also open a band gap (Klintenberg et al., 2010; Wei and Jacob, 2013a). One can also apply strain, poke holes, or form other planar carbon allotropes by rearranging carbon bonds, many of which open an appreciable band gap (~1 eV) in graphene at the GW level (Appelhans et al., 2010a,b; Liang et al., 2012; Nisar et al., 2012; Dvorak and Wu, 2015).

To fill the need for two-dimensional semiconductors, one can move away from graphene and consider materials that are intrinsically semiconducting. Elements from the third and fifth groups of the periodic table (III-V compounds) often form a semiconducting monolayer, see Figure 33. This is largely because of the A-B sublattice imbalance in these materials, which opens a band gap at the tight-binding level of theory (Wallace, 1947). A monolayer of hexagonal boron nitride (hBN) is one possibility, with a GW band gap of ~ 7.5 eV (Wirtz et al., 2005; Berseneva et al., 2013). Because graphene and hBN have similar lattice constants, they can be layered or hybridized easily, which gives additional tunability of the electronic properties (Bernardi et al., 2012). GaAs is another example, with a GW band gap of ~ 3 eV (Fakhrabad et al., 2014). Other III-V monolayers are also stable and have been studied with GW (Şahin et al., 2009; Wang and Shi, 2010; Fakhrabad et al., 2015; Prete et al., 2017). Phosphorene is a somewhat unusual case, as a monatomic group V material. This is reflected in its unusual structure, which has armchair-like vertical buckling, shown in Figure 33C. Phosphorene is attractive because it has a smaller band gap than many other two-dimensional semiconductors, computed to be ~ 2 eV with GW (Tran et al., 2014; Rudenko et al., 2015; Li et al., 2016b; Steinkasserer et al., 2016; Ferreira and Ribeiro, 2017), which is well-suited for applications.


Figure 33. Top and side views of (A) hexagonal boron nitride, (B) hydrogenated silicene (silicane), (C) phosphorene, (D) and 2H-MoS2. Structures taken from the Computational 2D Materials Database (Haastrup et al., 2018).

Finally, we get to the transition metal dichalcogenides (TMDs). TMDs have the chemical formula MX2 where M is a transition metal and X is a chalcogen, commonly S, Se, or Te. In their stable two-dimensional phase, TMDs usually form a three-layered structure with the transition metal atoms in a central layer between the chalcogens (called the 2H phase). MoS2, MoSe2, WS2, and WSe2 are all semiconductors with G0W0@LDA band gaps from 2.0−2.5 eV when including SOC (Rasmussen and Thygesen, 2015). The band structure of MoS2 is shown in Figure 34. TMDs feature unusual electronic structures derived from strong SOC and lack of inversion symmetry (see Manzeli et al., 2017). GW calculations at either the perturbative or partially self-consistent levels improve the agreement with experiment for fundamental band gaps (Cheiwchanchamnangij and Lambrecht, 2012; Komsa and Krasheninnikov, 2012; Ramasubramaniam, 2012; Espejo et al., 2013; Hüser et al., 2013a; Molina-Sánchez et al., 2013; Shi et al., 2013; Debbichi et al., 2014; Ugeda et al., 2014; Qiu et al., 2016; Robert et al., 2016; Lee et al., 2017). However, conclusions from different GW studies on the magnitude and character (direct or indirect) of the band gap in MoS2 are not entirely consistent. Depending on the level of self-consistency, truncation of Coulomb interaction, treatment of frequency dependence, and k-point sampling, the GW quasiparticle band gap of MoS2 can vary by ~ 0.44 eV (Qiu et al., 2016).


Figure 34. Band structure of MoS2 in the 2H phase at the PBE (black) and G0W0 (red) levels. The G0W0 bands include spin-orbit coupling but the PBE bands do not. The Fermi energies for each case are indicated by horizontal dotted lines. Data taken from the Computational 2D Materials Database (Haastrup et al., 2018).

The MoS2 case study highlights the importance of carefully converging GW calculations and the difficulties of two-dimensional materials, in particular. In two-dimensional semiconductors, the dielectric function is a linear function of q which results in very slow k-point convergence (Rasmussen et al., 2016). TMDs are also commonly stacked in layered materials called van der Waals heterostructures, which allow one to tune the electronic structure for device applications (Zhang et al., 2016; Arora et al., 2017; Winther and Thygesen, 2017). GW allows one to predict band alignment in these heterostructures from first principles (Ganesan et al., 2016).

9. Molecules

The application of GW to molecules is a relatively new field of research that has developed rapidly over the last decade. The electronic screening is much weaker in molecules than in extended systems. The low charge density in molecules does not naturally fit a screening interpretation of correlation which is intrinsic to GW and replacing the bare Coulomb potential with the dynamically screened Coulomb interaction W might not be the obvious choice. Even so, a rigorous test of GW for the He atom, with only two electrons, found excellent agreement with numerically exact results (Li et al., 2017). In addition, the first exploratory G0W0 studies on molecular systems revealed that the inclusion of screening at the GW level substantially improves electron removal and addition energies (Grossman et al., 2001; Niehaus et al., 2005; Dori et al., 2006; Ma et al., 2009, 2010b; Rostgaard et al., 2010; Blase et al., 2011; Ke, 2011).

9.1. First Ionization Potentials and Electron Affinities

In molecules, the single-particle states {ϕs0} correspond to molecular orbitals (MO) with discrete energies. The energy to remove an electron from an MO is referred to as ionization potential. The negative of the electron affinity (EA) corresponds to the energy needed to add an electron to the LUMO of the neutral system (− EALUMO = ϵLUMO), see also Equations (1) and (2). G0W0 provides access to both quantities. Furthermore, we can calculate the fundamental gap from the first ionization potential, IPHOMO, and the electron affinity

Δfgap=IPHOMOEALUMO.    (99)

The fundamental gap should not be confused with the optical gap Δogap, which is the energy needed for the charge neutral excitation from the HOMO to the LUMO. The optical gap is lower in energy than Δfgap and can not be obtained from GW. It defines the threshold for photons to be absorbed and for the formation of a bound electron-hole pair (exciton). Conversely, the fundamental gap is the energy threshold for the formation of a separate electron-hole pair, which is not bound together. It can be considered as the molecular equivalent to the band gap, see also Bredas (2014) and Baerends et al. (2013).

GW has been mainly applied to compute the IP for the HOMO and the electron affinity for π-conjugated molecules with potential for organic photovoltaic applications (Blase et al., 2011; Faber et al., 2011, 2012; Ke, 2011; Gallandi and Körzdörfer, 2015; Gallandi et al., 2016; Knight et al., 2016; Wilhelm et al., 2016; Marom, 2017). Examples for relevant π-conjugated organic molecules are linear acenes (linearly fused benzene rings), quinones, aromatic nitriles, anhydrides, porphyrins, and thiophene polymers. These classes of molecules are particularly suited as organic semiconductor because their EA is often positive10 (Richard et al., 2016), i.e., they are electron acceptors and their fundamental gap is much smaller than in inorganic molecules. For example, smaller acenes have gaps between 6.0 − 7.0 eV (Richard et al., 2016), whereas the fundamental gap of a small inorganic molecule like water is larger than 14.0 eV (van Setten et al., 2015).

The fundamental gap, IPHOMO and EA are critical parameters for the charge transport in organic semiconductors. Over the last years it has been shown that GW predicts these properties well. Using an appropriate starting point (see section 4.7), the reported mean absolute deviations (MADs) of IPHOMO and EA are less than 0.2 eV from the CCSD(T) reference (Gallandi et al., 2016; Knight et al., 2016). The MAD of IPHOMO with respect to experiment can be even reduced to <0.1 eV when including also vibrational effects in the GW spectra (Gallandi and Körzdörfer, 2015).

The electronic properties of π-conjugated molecular structures can be tuned by, e.g., increasing the chain length. It has been shown that GW correctly predicts the decrease of IPHOMO in trans-polyacetylene with increasing chain length (Pinheiro et al., 2015; Bois and Körzdörfer, 2017). Similar GW studies were conducted for band gaps of linear acenes (Wilhelm et al., 2016). The photovoltaic properties can be further modulated by using two different organic semiconductors in the cell: a molecule with a low IPHOMO (electron donor) and molecule with electron-acceptor character, i.e., with a high EA (Kippelen and Brédas, 2009). The level alignment of such donor-acceptor systems has been studied with GW for tetrathiafulvalene (TTF) and tetracyanoethylene (TCNE) or tetracyanoquinodimethane (TCNQ) dimers, demonstrating the importance of well-chosen starting points or self-consistent schemes (Caruso et al., 2014; Gallandi and Körzdörfer, 2015).

The accurate prediction of charged excitations is not only important for organic semiconductors, but also for DNA and RNA nucleobases in order to study their damage following exposure to ionizing radiation. IPs and EAs for these molecules have been reported at the G0W0 level in good agreement with experiment and quantum chemistry methods (Faber et al., 2011; Qian et al., 2011; Gallandi and Körzdörfer, 2015).

9.2. Ionization Spectra

The GW approximation has also been applied to calculate excitations of deeper valence states for small organic molecules (Körzdörfer and Marom, 2012; Marom et al., 2012; Caruso et al., 2013a; Egger et al., 2014; Ren et al., 2015) and also medium-sized π-conjugated molecules (Dori et al., 2006). An example is shown in Figure 35, where the ionization spectrum of pyridine is displayed for the first 12 valence states. Compared are the G0W0@PBE0 spectrum and the experimental PES. The positions of the peaks are in good agreement, in particular for the first three valence excitations. Benchmark studies for azabenzenes showed that a HF starting point yields distorted spectra, while hybrid DFT functionals and self-consistent schemes yield spectra that agree well with experiment (Marom et al., 2012). However, it has been found that the energy spacings and positions are not always reproduced satisfactorily. For example for benzene, the spacing of the HOMO-1 and HOMO-2 is vanishingly small for all starting points and also scGW (Ren et al., 2015). The exact spacing is larger than 0.5 eV. It has been demonstrated that a beyond GW scheme, so-called “vertex corrections,” are necessary to separate these two peaks (Ren et al., 2015).


Figure 35. Ionization spectrum of pyridine. G0W0@PBE0 QP energies compared to the experimental photoemission spectrum (Liu et al., 2011). The calculated spectrum has been artificially broadened; the position of the QP energies is indicated with vertical bars. All QP energies are extrapolated using the cc-pVnZ (n=3–6) basis sets, see Appendix C for further computational details. The QPs of the first valence states are colored in red, green and blue.

The deeper valence states are generally less valuable for characterization and chemical analysis. Core excitation energies, on the other hand, are a powerful tool to investigate the chemical structure of complex molecules and materials. They are element-specific, but are also sensitive to the atomic environment, such as covalent bonding, hybridization or the oxidation state (Siegbahn et al., 1969; Bagus et al., 1999, 2013). The application of GW to core states is more difficult than to valence states, as we will explain in more detail in section 11. Core excitations in GW are an emerging research field and appropriate numerical algorithms have only been developed recently (Golze et al., 2018).

Lastly, we will briefly address peak broadening in GW spectra. The G0W0 spectrum in Figure 35 has been artificially broadened to facilitate comparison with experiment. This broadening mimics vibrational, experimental resolution and finite lifetime effects. With regard to electronic lifetimes, also quasiparticle (QP) excitations in molecules have finite lifetimes accompanied by a finite broadening. Such a finite broadening would be revealed in the full spectral function A(ω), see Equation (3). The peaks close to the Fermi energy are usually sharp delta-like peaks, while higher energy excitations may decay through the formation of electron-hole pairs or collective excitations resulting in broader peaks, see Caruso et al. (2013a) for a detailed discussion of lifetimes of quasiparticles in molecules.

9.3. The GW100 Benchmark Set

An important aspect in electronic structure theory is benchmarking. Benchmark sets are very common in quantum chemistry, but have not found their way into GW until recently. Molecules offer a distinct advantage compared to solids for benchmarking because accurate reference energies can be computed with high-level quantum chemical methods. For this purpose, sets of small molecules are beneficial since they are computationally tractable. Moreover, they provide diversity in the electronic structure due to different types of covalent bonding.

The first systematic benchmarks were performed using a small set of 34 molecules (Rostgaard et al., 2010; Bruneval and Marques, 2013). Van Setten et al. took this idea further and proposed the GW100 benchmark set (van Setten et al., 2015), which is currently the largest and most popular GW benchmark set. It contains 100 molecules that feature a variety of elements from the periodic table. The original GW100 paper reports HOMO and LUMO quasiparticle energies computed at the G0W0@PBE level and the corresponding experimental references. Van Setten et al. used the test set for a quantitative comparison of the different GW methodologies implemented in the program packages Turbomole, FHI-aims and BerkeleyGW. They compared the performance of different basis sets (plane wave vs. localized), handling of core and valence electrons (all-electron vs. pseudopotentials) and different frequency integration techniques. The codes with localized basis sets (FHI-aims and Turbomole) agree to a precision of 1 meV for most molecules. The deviation of the BerkeleyGW plane wave code to the basis-set-extrapolated FHI-aims and Turbomole results is in the range of 200 meV. These numbers refer to the IPs obtained from full-frequency integration techniques available in all three codes. Based on this, van Setten et al. identified the basis set size as one important aspect for the accuracy of GW calculations.

The test set served later as a benchmark for the PAW G0W0 implementation in VASP (Maggio et al., 2017). This comparison established that the carefully converged PAW plane wave G0W0 calculations agree very well with the extrapolated results from the localized basis set codes. The MAD from the FHI-aims reference values is 60 meV. GW100 investigations with the WEST code gave similar results and highlighted the need for a re-evaluation of the pseudopotentials for some elements (Govoni and Galli, 2018). Moreover, the GW100 test set has been used to validate the accuracy of the low-scaling GW algorithm based in CP2K (Wilhelm et al., 2018). A comparison between the different codes is reported in Figure 36. Extrapolated values are represented in Figure 36A comparing plane wave codes to FHI-aims, whereas the comparison in Figure 36B is restricted to codes with localized functions. A list of all codes that ran the GW100 benchmark can be found in GW100 (2018).


Figure 36. GW100 benchmark comparing IPHOMO energies computed at the G0W0@PBE level. FHI-aims is set as reference: ΔIPHOMO = IPHOMO(FHI-aims)−IPHOMO(X). (A) Comparison of extrapolated/converged results for VASP (Maggio et al., 2017), WEST (Govoni and Galli, 2018), BerkeleyGW (van Setten et al., 2015). Shown are the results from full-frequency treatments and iterative solutions of the QP equation. (B) Comparison of localized basis set codes using the Gaussian basis set def2-QZVP (Weigend and Ahlrichs, 2005) for Turbomole (no-RI) (van Setten et al., 2015) and the N3 implementation in CP2K (Wilhelm et al., 2018). Note that BN, O3, MgO, BeO, and CuCN are excluded for WEST, VASP and CP2K and that the BerkeleyGW and Turbomole data contain only a subset of 19 and 70 molecules, respectively. Box plot: Outliers represented by dots; boxes indicate the “interquartile range” measuring where the bulk of the data are.

The GW100 test set was not only used to validate the reliability of numerical techniques in G0W0 implementations. It has been also used for a comprehensive assessment of different self-consistent GW methodologies: scGW, QSGW and scGW0 (Caruso et al., 2016). The results were compared to CCDS(T) at the polarized triple-zeta level reporting the smallest discrepancies for QSGW. A comparison of basis set extrapolated CCSD(T) and GW schemes was performed shortly afterwards for a smaller, more specialized benchmark set of 24 organic electron-acceptor molecules, where G0W0 based on long-range corrected hybrid functionals emerged as the best GW method (Gallandi et al., 2016; Knight et al., 2016). Since then, also equation of motion (EOM) coupled cluster benchmark sets have been published that provide reference spectra (and not just HOMO or LUMO energies) for molecules (Lange and Berkelbach, 2018; Ranasinghe et al., 2019).

9.4. Molecular Crystals

Modern applications of GW comprise not only isolated molecules, but also molecules in the condensed-phase, such as organic molecular crystals. These materials are composed of weakly bonded molecular units held together by, e.g., van-der-Waals interactions, dipole-dipole interactions or hydrogen bonds. Here, we summarize only some key applications of GW to molecular solids. A more comprehensive discussion can be found in the specialized review by Kronik and Neaton (2016).

Molecular solids exhibit a band gap renormalization similar to molecular adsorbates discussed in section 7. The band gap of molecular solids is significantly smaller than the fundamental gap of the isolated molecules (Sato et al., 1981). As for molecular adsorbates, the gap renormalization is a direct consequence of polarization effects. It is also present when there is no wave-function overlap between neighboring molecular units. If an electron is added to or removed from a certain molecule, the new charge carrier is screened not only by the molecule it was added to, but also by the surrounding molecules. This renormalization effect is shown in Figure 37 for the benzene crystal. The HOMO level moves up in energy with respect to its position in the gas phase molecule, whereas the LUMO moves down resulting in a gap reduction.


Figure 37. Fundamental gaps of gas-phase benzene and band gap of the benzene crystal (space group Pbca). PBE was used as starting point for the G0W0 calculations. Data retrieved from Refaely-Abramson et al. (2013).

The gap renormalization typically lies in the range of 2–6 eV (Kronik and Neaton, 2016) and has been studied with GW for benzene (Refaely-Abramson et al., 2013), corannulene-based materials (Zoppi et al., 2011), C60 (Refaely-Abramson et al., 2013), pentacene (Sharifzadeh et al., 2012; Refaely-Abramson et al., 2013), perylenetetracarboxylic dianhydride (PTCDA) (Sharifzadeh et al., 2012), octaethylporphyrin (H2OEP) (Marsili et al., 2014), 6,13-bis(triisopropylsilylethynyl)-pentacene (TIPS-pentacene) (Sharifzadeh et al., 2015) and oligoacenes (Rangel et al., 2016). The gap reduction is not captured by standard DFT calculations (Refaely-Abramson et al., 2013), see also Figure 37. In fact, the DFT gap remains almost unchanged when transitioning from the gas to the crystalline phase because the long-range polarization effects responsible for the gap renormalization are not included in conventional DFT functionals.

The molecular orbitals of molecular crystals resemble those of an isolated molecule. However, the overlap between neighboring molecules is not zero resulting in a k dependence (dispersion) of the energy levels. Starting with early work on C60 (Shirley and Louie, 1993), GW band structures have been reported for a wide range of organic crystals (Tiago et al., 2003b; Sharifzadeh et al., 2012, 2015; Refaely-Abramson et al., 2013, 2015; Fonari et al., 2014; Yanagisawa and Hamada, 2017; Cocchi et al., 2018; Rangel et al., 2018). As for inorganic semiconductors, GW opens the band gap and increases the band width with respect to DFT. For example, GW band widths reported for pentacene (Tiago et al., 2003b; Sharifzadeh et al., 2012), PTCDA (Sharifzadeh et al., 2012), rubrene (Yanagisawa et al., 2013), or picene (Yanagisawa et al., 2014) are larger by more than 15%. The bands of molecular crystals are relatively flat compared to inorganic semiconductors (see section 6). For example, GW-computed band widths for pentacene are only 0.4 eV for the valence and 0.7 eV for the conduction band (Sharifzadeh et al., 2012).

Molecular crystals are an ideal testbed for GW embedding schemes since the band gap of molecular solids is mainly determined by polarization effects and significantly less by dispersion. In the spirit of quantum mechanics/molecular mechanics (QM/MM) embedding schemes the molecular crystal is partitioned into a small part that is calculated with GW and a much larger MM part. In the embedding scheme proposed by Blase and co-workers, the small part to which GW is applied consists of one or more molecules, while a continuum polarization model is used to include the response of the MM system (Duchemin et al., 2016; Li et al., 2016a). They reported GW/MM gaps for pentacene and perfluoropentacene that are in close agreement with the bulk reference (Li et al., 2018). Such embedding schemes are often computationally more efficient than periodic boundary condition calculations, especially for local orbital basis set codes.

10. Total Energy and the Electronic Ground State

In addition to the quasiparticle spectrum, the Green's function also provides information on the electronic ground state. Both the ground state density and the ground state total energy are accessible. However, very few studies have explored ground state properties with GW. Since this review mainly addresses spectroscopic properties, we will only briefly address GW ground state calculations here.

10.1. Electron Density

The ground-state density n(r) follows directly from the Green's function (Fetter and Walecka, 1971)

n(r)=iσGσ(r,r,t=0).    (100)

The total electron number contained in G can be obtained through integration of the density. For a self-consistent G that has been obtained from a converged solution of Dyson's equation (Schindlmayr, 1997), this number should then equal the total number of electrons N in the system. Also scGW0 satisfies this particle number conservation law, but all other approximate self-consistency schemes as well as G0W0 violate particle number conservation.

Figure 38 shows density differences compared to the Hartree-Fock method for PBE, coupled cluster singles-doubles (CCSD), and self-consistent GW for the CO molecule. Overall the pattern is similar. All three methods remove charge from the bonding region and the top of the oxygen atom and focus it on the carbon atom and a p orbital of the oxygen atom. The charge density difference pattern between CCSD, a high-level quantum chemistry method, and scGW is very similar. This indicates that the GW density is of high quality.


Figure 38. Density difference for the CO molecule between Hartree-Fock (HF) and PBE (left), coupled cluster singles-doubles (CCSD) and self-consistent GW (right). Charge depletion in the three methods is encoded by blue and charge accumulation by red colors. The same computational settings as in Caruso et al. (2013a) have been used.

From the density, the dipole moment of CO can be calculated. In PBE the dipole moment amounts to 0.2 Debye, in HF to −0.13 Debye and from scGW we obtain 0.07 Debye (Caruso et al., 2012a). The CCSD dipole moment is 0.06 Debye (Caruso et al., 2012b). All values were computed at the equilibrium bond-length of the respective method and the experimental dipole moment is 0.11 Debye (NIST, 2019). CCSD and scGW again agree closely and also match experiment reasonably well, whereas PBE overestimates the dipole moment and HF gives the wrong sign. The good agreement between scGW and CCSD and experiment is further testimony for the quality of the GW density.

Since fully self-consistent GW calculations are numerically quite involved and can currently only be performed for small systems, DFT densities are still used in the majority of GW studies. However, in situations in which the underlying DFT Kohn-Sham spectrum has the wrong order of states, erroneous charge transfer can occur in the DFT calculation. This is, for example, frequently the case in molecular complexes, if the HOMO of one molecule erroneously ends up above the LUMO of another. The corresponding G0 will not reflect the true ground state density of the complex and the subsequent G0W0 calculation will be wrong. G0W0 itself cannot rectify this situation because it has no access to the density. Only self-consistent schemes can correct the density and the Green's function. Examples of such molecular complexes are dimers of tetrathiafulvalene (TTF) with tetracyanoethylene (TCNE), tetracyanoquinodimethane (TCNQ) and p-chloranil. In all cases, scGW stops the erroneous charge transfer that occurs in PBE and in hybrid functionals with a low amount of exact exchange (Caruso et al., 2014). The resulting charge density reflects the molecular charge densities that are slightly perturbed where the molecules are closest to each other.

10.2. Total Energy

The total electronic energy can be obtained from the single-particle Green's function G via the Galitskii-Migdal (GM) formula: (Galitskii and Migdal, 1958; Fetter and Walecka, 1971)

EGM=iσ dr dtlimrrtt+[it+h^0]Gσ(rt,rt),    (101)

where ĥ0 contains the kinetic energy operator and the external potential. This equation can be recast into a more familiar looking form (Strinati, 1988; Caruso et al., 2013a)

Etot[G]=T[G]+Eext[G]+EH[G]+Exc[G],    (102)

in which T denotes the kinetic energy, Eext the external potential energy, and EH the Hartree energy. The exchange-correlation (xc) energy

Exc[G]=0dω2πTr{Σ(iω)G(iω)},    (103)

is given by the self-energy, Σ, and the Green's function. Equation (102) is appealing because it contains the same terms as the DFT total energy. Notable differences are that the kinetic energy is the fully interacting kinetic energy and not that of an auxiliary non-interacting system. Correspondingly, the exchange-correlation energy is purely due to electronic exchange and correlation and does not need to also approximate the difference between the interacting and the non-interacting kinetic energy as in Kohn-Sham DFT.

The GW total energy is closely related to the popular random-phase approximation (RPA) in DFT (Langreth and Perdew, 1977; Hesselmann and Görling, 2010; Eshuis et al., 2012; Ren et al., 2012b). The xc energy in GW and RPA can be represented in terms of topologically identical Feynman diagrams (Hellgren and von Barth, 2007; Caruso et al., 2013b) and thus have a total energy expression with the same functional dependence on the Green's function (Klein, 1961; Dahlen et al., 2006b; Hellgren and von Barth, 2007). However, the RPA energy is evaluated with a non-interacting Green's function (originating from a local Kohn-Sham potential) and the GW energy with a fully interacting Green's function. In fact, the Dyson equation results as stationary equation from the optimization of the GW total energy with respect to the Green's function in the Klein or Luttinger-Ward functionals.

Early GW calculations for the homogeneous electron gas found the total energy to be in good agreement with Quantum Monte Carlo calculations (von Barth and Holm, 1996; Holm and von Barth, 1998; Holm, 1999; Holm and Aryasetiawan, 2000; García-González and Godby, 2001). GW also captures van der Waals interactions as exemplified by the total energy curve between two jellium slabs (García-González and Godby, 2002) and by changes in the GW density of the argon dimer (Ferri et al., 2015). More recently it was shown that the lattice constants and bulk moduli of simple solids agree much better between experiment and GW than with LDAs, GGAs or HF (Kutepov et al., 2009). However, GW total energy calculations for atoms (see Figure 39) and small molecules show the opposite (Stan et al., 2009; Caruso et al., 2012a, 2013a). Presumably due to the low amount of screening self-consistent GW calculations are outperformed by high-level quantum chemistry methods and even simple DFT functionals.


Figure 39. Total energy of atoms computed with three different GW variants for atoms and small molecules plotted as a difference to the essentially exact Configuration Interaction (CI) results. Data retrieved from Caruso et al. (2012a).

Further analysis (Hellgren et al., 2015) reveals that the difference between G0W0@HF and scGW can be ascribed to the difference in the kinetic energy (termed here kinetic correlation in analogy with DFT) because their Coulomb correlation energies are almost identical for the small systems shown in Figure 39. Conversely, the difference between G0W0@PBE and scGW is almost entirely due to Coulomb correlation. Both Coulomb and kinetic correlation are large, as illustrated in Figure 39. Once included, the remaining difference between scGW and full configuration interaction (the essentially exact solution) must be due to missing vertex corrections. This contribution is much smaller than the two correlation contributions.

11. Current Challenges and Beyond GW

11.1. Challenges

As successful as the GW approximation is for describing quasiparticle excitations, there are still technical and theoretical challenges to overcome. Core-level spectroscopy is a valuable tool for chemical analysis and characterizing materials. The operating principle is the same as PES and IPES discussed in section 2, though at higher incident energies using X-rays. The technique is then referred to as X-ray photoelectron spectroscopy (XPS). Core levels of the same type, for example, different carbon 1s states, are element-specific, but are also sensitive to the local chemical environment, i.e., bonding, hybridization or the oxidation state (Egelhoff, 1987; Bagus et al., 1999). However, these so-called chemical shifts are for second-row elements often smaller than 1 eV (Siegbahn et al., 1969). The energetic differences are particularly minute for carbon with XPS peaks that are separated by less than 0.5 eV. Such spectra are hard to resolve and interpret. Theoretical spectroscopy can be a valuable tool to aid the interpretation of experimental results.

For core levels, however, the simple, single quasiparticle picture can break down. The incident photon in PES may produce spectral features away from the single-particle peak. If the additional peak is broad, these so-called satellites can be attributed to the collective excitation of the system after the electron is excited. If the additional peak is narrow or, equivalently, has a long lifetime, the electron has spectral weight divided between multiple particle-like eigenstates of the system (Golze et al., 2018). This effect can also appear when probing the multiplet structure of open-shell systems (Lischner et al., 2012). In these cases, the quasiparticle equation can have multiple solutions, making both the GW calculation and interpretation of the result more difficult. The problem also appears for more conventional valence states of small molecules, and recent work has shown that these multiple solutions lead to unphysical discontinuities in quasiparticle energies and that evGW can exacerbate the problem (Loos et al., 2018; Véril et al., 2018). Just as for experimental spectroscopy, the sensitivity of core states to the local environment makes the GW calculation more challenging than for conventional valence states. Due to its value for chemical analysis and dearly needed support for XPS experiments, GW for core levels can yield useful insight and is an ongoing topic of research (Zhou et al., 2015; Golze et al., 2018).

Spin dependence in GW calculations is important for understanding magnetic systems and is critical to the electronic structure of topological insulators. Already in the case of collinear spin, when the spin quantum number is either up or down, spin polarization has an effect on the excitation spectrum of MnO (Rödl et al., 2008). By including spin-orbit coupling (SOC) in the one electron Hamiltonian, single-particle states become noncollinear and can no longer be decomposed into up or down. Noncollinear calculations are important in relativistic systems with strong SOC or when describing magnetic effects (Sakuma et al., 2011; Kutepov et al., 2012; Ahmed et al., 2014; Kühn and Weigend, 2015; Scherpelz et al., 2016). For materials with heavy elements, energy shifts due to spin-orbit coupling must be included for good agreement with experiment on band gaps (Scherpelz et al., 2016). Topological insulators commonly contain heavy elements (Se, Te, Bi, Sb) and depend on spin-orbit coupling for band inversion (Aguilera et al., 2013a,b, 2015a; Nechaev and Chulkov, 2013; Nechaev et al., 2015). For a detailed review of GW+SOC calculations, see Aguilera et al. (2015b). To describe spin-dependent interactions between particles, one must generalize Hedin's equations beyond the Coulomb interaction, which has no spin dependence. This generalization was recently completed (Aryasetiawan and Biermann, 2008, 2009) and allows one to treat magnetic dipole-dipole interactions, for example.

11.2. Quantum Chemistry

Quantum chemistry offers an established, albeit expensive, route to compute particle addition/removal energies in molecules. Ionization energies and electron affinities can be computed as the difference of total energies between the neutral molecule and the ion. In fact, GW calculations on small systems are often compared with coupled cluster results as a benchmark. A direct comparison between correlated wave-function and Green's function methods to determine the level of correlation described by each is somewhat challenging. In certain cases, it is possible: recent work compares diagrams included in GW with those included in equation-of-motion coupled cluster theory (Lange and Berkelbach, 2018).

Generally, nonperturbative wave-function methods are considered beyond GW, even if they rely on an ansatz or other approximation. In Green's function embedding theories, quantum chemistry (either full or truncated configuration interaction) can be used as a high accuracy Green's function solver in a subspace (Zgid et al., 2012). After computing the subspace wave function, one directly computes the amplitudes in Equation (5) for the subspace Green's function. With G and G0 in hand, it is then trivial to compute the self-energy (Pavlyukh and Hübner, 2007). In this subspace, the Green's function is computed from accurate many-body wave functions so that correlation is treated beyond GW. The subspace Green's function can be self-consistently iterated with the remaining degrees of freedom described at the GW level of theory (Martin et al., 2016). Other routes to combine GW with quantum chemistry are an emerging field. A newly developed method combines GW with configuration interaction by embedding a wave function calculation inside of a Green's function calculation (Dvorak et al., 2018; Dvorak and Rinke, 2019). These developments offer valuable insight to merge these disciplines in the future.

Green's functions are also directly studied in quantum chemistry, where they are more commonly called propagators. There certainly is some overlap between the two communities in their treatment of GW or GW-like approximations. Because we primarily focus on GW and Hedin's equations in physics, we refer the interested reader to the work of Cederbaum and Domcke (2007) and Ortiz (2012) for a perspective of propagators in chemistry.

11.3. Non-equilibrium Green's Functions

The GW approach has also been applied to systems in strong external fields. These include quantum transport calculations (Thygesen and Rubio, 2008; Spataru et al., 2009) and semiconductors in strong laser fields (Spataru et al., 2004a). The problem of describing quantum transport is similar to that of level alignment at a molecule/metal interface discussed in section 7. First, the alignment of molecular states in the contact region relative to the Fermi level of the metal leads determines the overall conductance. Second, for applied biases, charge will flow from the lead into the molecule or molecules in the contact region. This charge flow will alter the electron density of the system and therefore the quantum mechanical interactions.

Self-consistent GW calculations (Thygesen and Rubio, 2008) take charge transfer and the associated change in screening (e.g., image effect) and the many-body interactions into account correctly. scGW is an appropriate tool for finite, small bias quantum transport calculations, as benchmarked for instance for thiol- and amine-linked benzene/gold (Strange et al., 2011) and alkane/gold junctions (Strange and Thygesen, 2011). Strong correlation effects in quantum transport, such as the Coulomb blockade or the Kondo effect, can, however, not be captured with the GW approach and require a beyond GW treatment of correlation (Spataru et al., 2009; Thoss and Evers, 2018).

For strong external biases in quantum transport and systems in other strong external fields, the electron distribution is perturbed so strongly that it can no longer be described by an equilibrium Fermi function at a finite temperature. In such non-equilibrium cases, the Green's function theory described in this review article is not applicable anymore. Non-equilibrium scenarios can be incorporated into the Green's function formalism, by switching to non-equilibrium Green's functions defined on the Keldysh contour (Dahlen et al., 2006a). These non-equilibrium Green's functions obey the Kadanoff-Baym and not Hedin's equations. The Keldysh contour formalism goes beyond the scope of this review article and we refer the interested reader to two excellent recent books (Stefanucci and Leeuwen, 2013; Karlsson and Leeuwen, 2018).

In one application of this non-equilibrium formalism highly excited GaAs was investigated. It had been hypothesized that GaAs could become metallic if enough electrons could be promoted from the valence to the conduction band with a strong laser. Non-equilibrium GW calculations showed that the band gap was indeed decreasing with increasing laser power, but would not close completely, falsifying the hypothesis (Spataru et al., 2004a).

11.4. Vertex Corrections

To go beyond the GW approximation, one must include vertex corrections. The full set of Hedin's equations, including the vertex, are shown diagramatically in Figure 40, which can be directly compared to Figure 9. The mathematical equations are in Appendix B.2. By comparison to the GW diagrams, we see that treating the vertex, Γ, beyond a single spacetime point significantly complicates the equations, as demonstrated for a single diagram in Figure 8. The vertex contributions beyond Γ(1, 2, 3) = δ(1, 2)δ(1, 3) are commonly called vertex corrections. The exact vertex requires the functional derivative δΣ/δG. This functional derivative cannot be computed numerically and must be derived analytically. Any resulting vertex is extremely expensive to compute because it now depends on three spatial, spin, and time indices. There are a few reasonable options for reducing the expense: using an approximate Σ from a different theory to simplify the derivative, using a diagrammatic but simplified Γ, or using an exact Γ in only a small subspace. As in GW, Γ can be selectively iterated or computed in a single-shot way to further lower the cost.


Figure 40. Diagrammatic representation of Hedin's equations. All 5 quantities are coupled to all others. Here, we omit the Hartree potential from the G diagram, though it must also be included when translating the diagrams to the equations in Appendix B.

The earliest approaches to vertex corrections used an approximate vertex based on the local density approximation (LDA) to Kohn-Sham DFT (Hybertsen and Louie, 1986; Del Sole et al., 1994). In these approximations, the LDA exchange-correlation functional is used in place of the self-energy to compute the functional derivative. The approximate vertex enters the polarizability and interaction as

fxc(1,2)=δvxc(1)δn(2)W˜=v[1χ0(v+fxc)]1    (104)

where n is the electron density and fxc determines the vertex correction. The advantage of the LDA is that fxc can be calculated analytically.

These approaches are computationally much lighter than the true vertex and, for that reason, are still used (Schmidt et al., 2017). Approximate vertex corrections can also be extended beyond the LDA to recover a more realistic behavior (Chen and Pasquarello, 2015a). The inclusion of fxc has its roots in time-dependent density functional theory (TDDFT) and is somewhat inconsistent with the Green's function formalism. The final results of such calculations can improve band gaps compared to G0W0 (Chen and Pasquarello, 2015a) or band centers compared to G0W0 (Schmidt et al., 2017). Shaltaf et al. found that a GGA-based vertex correction had a negligible effect on band offsets in the Si/SiO2 interface (Shaltaf et al., 2008).

Diagrammatic vertex corrections, instead of those based on a density functional, are a more consistent and formal extension to GW. One can build upon GW in a fully diagrammatic framework by including vertex corrections in the perturbative, single-shot self-energy. These methods are analogous to the G0W0 approach in that the vertex correction is computed only one time, and the quasiparticle energies are usually computed in a diagonal approximation. For example, second order screened exchange (SOSEX) and related approximations include a subset of exchange-type diagrams which are a vertex correction to GW (Shirley and Martin, 1993b; Bobbert and van Haeringen, 1994). In the approximation of Ren et al. (2015), the diagonal matrix elements of GW+SOSEX are

     ϕs0|ΣcGW(iω)+ΣcSOSEX(iω)|ϕs0=12πdω qrp(fqfr)sr||qpqp|W(iω)|sr(iω+iωϵp0)(iω+ϵq0ϵr0),    (105)

where fq and fr are Fermi occupation factors and 〈sr| |qp〉 is defined in Appendix A. Notice extra matrix elements in the numerator and factors in the denominator compared to the equation for the GW self-energy in Figure 10. SOSEX-type approximations generally improve upon GW band gaps in molecules (Ren et al., 2015). In the perturbative approach, these calculations are relatively lightweight but have the same starting point dependence of G0W0.

A systematic bridge between diagrammatic vertex corrections and TDDFT was developed by Del Sole, Reining, and others (Streitenberger, 1984a,b; Reining et al., 2002; Adragna et al., 2003; Del Sole et al., 2003; Sottile et al., 2003; Bruneval et al., 2005). In this approach, the kernel to construct the irreducible polarizability is cast as only a two-point function. This is in contrast to the exact vertex, Γ, which depends on four spacetime coordinates to compute, making it much more expensive (four coordinates for the derivative δΣ(1, 2)/δG(4, 5), see Appendix B.2). This two-point kernel can only be used inside of W. Outside of W, the three-point nature of the vertex is unavoidable. The simplified many-body approach retains the simplicity of a TDDFT kernel but has its foundation in many-body theory. By adopting the GW approximation to Σ, an approximate, analytic fxcQP exists. Calculations of the dielectric function in Si and GaAs show good agreement between the fxcQP approach and a solution for the full vertex (Adragna et al., 2003).

More recent work has included diagrammatic vertex corrections to solve Hedin's equations at some level of self-consistency, though still at an approximate level (Grüneis et al., 2014; Kutepov, 2016, 2017; Maggio and Kresse, 2017). The greatest conceptual and computational difficulty to these calculations is how to update Γ. Because Γ enters in both χ0 and Σ, and because full self-consistency is so expensive, it is advantageous to update Γ in only one portion of the calculation. For example, one could evaluate the interaction W with a diagrammatic Γ once at the beginning of the calculation. Keeping W fixed, only G and Σ are updated through Dyson's equation in the self-consistency cycle. While this procedure is only partially self-consistent, it incorporates a diagrammatic Γ while keeping the GW level of complexity through the self-consistency cycle. When applied to semiconductors and insulators, and with some practical restrictions on Γ, solutions of Hedin's full equations give noticeably better band gaps than GW (Kutepov, 2016, 2017). Full, self-consistent solutions of Hedin's exact equations remain out of reach in real systems, and even partially self-consistent schemes are a technical challenge.

11.5. Optical Properties

Calculations of the many-body vertex have another application beyond particle addition/removal energies. Optical processes, such as photon absorption, can also be modeled in the Green's function formalism. In such a case, the relevant correlation function is the time-ordered two-particle correlation function, L. L obeys a Dyson equation like G, except that the role of the self-energy is instead played by δΣ/δG. The Dyson series for the full vertex to determine L is called the Bethe-Salpeter equation (BSE) in physics (Salpeter and Bethe, 1951; Held et al., 2011). BSE calculations describe a different process than ordinary GW, so they are not beyond GW in the same sense as including a vertex in the self-energy. Even so, they are closely related. The common implementation of the BSE for materials relies on the GW approximation to the self-energy. In these GW/BSE calculations, the excited electron and hole instantaneously interact via W.

The effective two-particle Hamiltonian for correlated optical excitations (Strinati, 1982, 1984; Rohlfing and Louie, 2000), called excitons, is

(ϵaϵi)AiaP+iaia|K|iaAia=ΩPAiaia|K|ia=Wiaia+viiaa    (106)

where i and a denote again occupied and empty states, respectively, ϵi/a is a quasiparticle energy, AiaP is the Pth exciton wave function in the single-particle basis, and ΩP is the Pth excitation energy. Equation (106) makes the common Tamm-Dancoff approximation (TDA), which ignores backward propagating electron-hole pairs that are present in the exact BSE. Matrix elements of K, δΣ/δG, include a direct screened interaction (W) and repulsive exchange (v). Schematic and diagrammatic representations of the optical process are shown in Figure 41. BSE calculations can be considered a first iteration of Γ in Hedin's equations to go beyond GW for particle addition/removal energies, if the resulting polarizability is reinserted into W.


Figure 41. (A) Schematic representation of optical absorption. (B) Diagrammatic representation of GW/BSE. The electron and hole are represented by G lines, computed in the GW approximation, and their direct interaction is through the screened Coulomb interaction.

The first BSE calculations included only the bare electron-hole exchange (Hanke and Sham, 1975) in a semi-empirical basis and were then extended to include the direct, screened interaction (Hanke and Sham, 1979, 1980). Ab-initio GW/BSE calculations focused on semiconductors like Si, GaAs, and Li2O where GW/BSE produces optical absorption spectra and exciton binding energies in close agreement with experiment (Onida et al., 1995b; Albrecht et al., 1997, 1998a,b; Benedict et al., 1998a,b; Rohlfing and Louie, 1998, 2000; Benedict and Shirley, 1999). Similar to the proliferation of GW since its early successes, GW/BSE has been applied extensively to solids (Schleife et al., 2011, 2018; Rinke et al., 2012; Erhart et al., 2014), molecules (Körbel et al., 2014; Bruneval et al., 2015; Jacquemin et al., 2015), surfaces (Palummo et al., 2004), and two-dimensional materials (Komsa and Krasheninnikov, 2012; Ramasubramaniam, 2012; Hüser et al., 2013a; Qiu et al., 2013, 2016; Shi et al., 2013; Ugeda et al., 2014; Dvorak and Wu, 2015). As with Dyson's equation, equations with the Bethe-Salpeter form appear in different contexts in many-body theory. For example, a Bethe-Salpeter equation can also describe spin-flip excitations in magnetic materials (Müller et al., 2016).

11.6. T-Matrix

The framework of Hedin's equations, and GW in particular, places great emphasis on the screened Coulomb interaction. Indeed, many of the approximate schemes presented here frame the exact vertex as a correction (hence the term “vertex correction”) to a GW calculation of the self-energy. In certain systems, it may be necessary to abandon this picture entirely. For example, in systems with low electron density and a similarly low number of electron-hole screening channels, as in very small atoms or molecules, screening of the Coulomb interaction may be insignificant. Roughly speaking, diagrams of the vertex type could be more important than screening diagrams included in GW. For these systems, we should adopt a different formalism which does not rely on screening and directly emphasizes other correlation channels. One such formalism is the T-matrix (Zhukov et al., 2004; Romaniello et al., 2012; Zhang et al., 2017), in which the self-energy is written as a product of G with a four-point kernel, T,

Σ(1,2)=id3 d4 G(4,3)T(1,3,2,4).    (107)

The precise form of T depends on the choice of the particle-particle (pp) or particle-hole (ph) T-matrix, which determines the channels that are correlated alongside a third propagating channel. T obeys its own Dyson series and physically corresponds to repeated interactions, or scattering, between the selected channels (pp or ph). There are advantages of the T-matrix approach: it is exact up to second order in the bare interaction and includes many more exchange diagrams than GW, making it useful for magnetic systems. At first glance, the ph T-matrix may sound like GW. However, the two approximations correlate different particle and hole channels in the self-energy diagram. A schematic representation of the correlated channels in GW and T-matrix is shown in Figure 42. Notice the different topologies of G0 lines correlated in GW and ph-T. The T-matrix approach has been applied to understand the role of spin-flip excitations in metals (Zhukov et al., 2004, 2005, 2006; Müller et al., 2018; Młyńczak et al., 2019), double ionizations and Auger spectroscopy (Noguchi et al., 2007, 2008, 2010), as well as satellites in metals (Springer et al., 1998).


Figure 42. Schematic representation of diagrams included with (A) GW, (B) particle-hole T-matrix, and (C) particle-particle T-matrix. In each case, channels going into the box are correlated further with additional interactions.

11.7. Cumulant Expansion

One long-standing problem with the GW approximation is its description of plasmon satellites in, for example, Si, Na and Al (Aryasetiawan et al., 1996; Guzzo et al., 2011, 2014; Zhou et al., 2015). Plasmon satellites are peaks in the spectral function which are attributed not to a single quasiparticle, but to the coupling between a hole (in the particle removal case) and the collective excitation of the remaining electrons. This coupling leads to a quasiparticle peak in the spectral function and a series of progressively weaker plasmon replicas separated by the plasmon energy.

A proven route to improve the plasmon description compared to GW is the cumulant expansion to the Green's function, which has been tested on Na, Al, graphene, Si, and the electron gas (Aryasetiawan et al., 1996; Guzzo et al., 2011; Gatti and Guzzo, 2013; Lischner et al., 2013, 2014; Caruso and Giustino, 2015, 2016; Caruso et al., 2015; Zhou et al., 2018). Based on an exponential ansatz, somewhat analogous to the coupled cluster expansion for the wave function, the cumulant Green's function for a hole takes the form (Aryasetiawan et al., 1996; Guzzo et al., 2011; Lischner et al., 2013)

Gs(t)=Θ(t)eiϵs0t+Cs(t)    (108)

where ϵs0 is the mean-field energy that enters G0 for state s and Cs(t) is the cumulant. The exact form of the cumulant Cs(t) depends on the chosen approximation to the self-energy. If one Taylor expands Equation (108) in powers of Cs and compares it to Dyson's equation with the GW self-energy, an approximate closed form for the cumulant exists. This is called the GW+C method. The cumulant includes vertex corrections beyond the GW self-energy at the same computational expense as ordinary GW. These vertex corrections generally improve the description of satellites over GW when compared with experiment.

The cumulant appears to be a tremendous success − it miraculously provides vertex corrections for the same cost as GW. However, it does not improve the description of valence quasiparticle energies. The quasiparticle energy is still determined by ordinary GW. Furthermore, the cumulant ansatz in Equation (108) relies on the separation of electron and hole branches of the Green's function and produces satellites only below the Fermi energy. In general, this separation is not correct, and it becomes a worse approximation closer to the Fermi energy (Guzzo et al., 2014; Martin et al., 2016). The formal connection between GW and the cumulant is presented in Gumhalter et al. (2016).

As a case study of the cumulant, we highlight the study of doped graphene by Lischner et al. (2013). The spectral function of doped graphene on a SiC(0001) surface displays a quasiparticle peak and satellite, shown in Figure 43. With ordinary G0W0, the splitting between the quasiparticle and satellite is 0.44 eV, which overestimates the experimental value of 0.3 eV. By including vertex corrections to the hole Green's function with the cumulant, the GW+C calculation reduces the splitting to 0.27 eV. GW+C also redistributes spectral weight away from the quasiparticle and to the satellite. Additionally, GW+C eliminates a spurious plasmaron − coupling between a hole and plasmon − solution that appears in GW.


Figure 43. Spectral function A(ω) of graphene on SiC at the Dirac point for electron doping density of n = −5.9 × 1013 cm−2. Data taken from Lischner et al. (2013).

11.8. Local Vertex

The treatment of localized electrons in physics has become its own subfield (Held et al., 2011; Hirayama et al., 2017; Tomczak et al., 2017). In strongly-correlated physics, localized d- or f-electrons usually indicate a need to go beyond GW (Nohara et al., 2009; Gatti and Guzzo, 2013; Sakuma et al., 2013), with transition metal oxides being typical test cases. Much of the previous discussion applies just as well to localized electrons as any others − Hedin's equations are exact. However, the localized nature of these states lends themselves to model Hamiltonians, particularly the Hubbard model (Hubbard, 1963), which describes localized interactions by a parameter U. U is a measure of the repulsive interaction, or energetic cost, for electrons occupying the same spatial orbital. When combined with the LDA, the LDA+U method can improve upon the poor description of localized states by mean-field and GW theories (Jiang, 2018).

In the Green's function formalism, including diagrams beyond GW is made tractable by approximating the true Coulomb interaction by U. The combination of GW with dynamical mean-field theory (DMFT) (Georges and Kotliar, 1992; Kotliar et al., 2006), the GW+DMFT method (Biermann et al., 2003; Biermann, 2014), is a rigorous way of combining diagrams of higher order with the GW approximation. GW+DMFT describes long-range correlation with GW and local d- or f-electron correlation via the Hubbard U. GW+DMFT correlates a small set of states (the d- or f-electrons) using a non-perturbative vertex, called a “local” vertex since single site DMFT only includes diagrams which are local in orbital space. The first studies of SrVO3 with GW+DMFT demonstrated its potential for predicting spectral properties of strongly-correlated solids (Tomczak et al., 2012, 2014). The GW+DMFT method continues to be developed (Biermann, 2014; Choi et al., 2016).

12. Conclusion

Photoelectron and inverse photoelectron spectroscopy will remain some of the most powerful probes of matter available to scientists. While experimental spectroscopy gives the “answer” in the form of the measured spectrum, it may not give the full understanding of the underlying physics. In this regard, theoretical spectroscopy plays a huge role as a complement to experimental techniques.

We have introduced the Green's function formalism and many-body theory from a perturbation theory perspective. The formalism is exact, in principle, and allows one to calculate both ground and excited state properties. From the Feynman diagram construction, we have given a heuristic motivation for Hedin's equations, which are themselves nonperturbative. Hedin's equations place emphasis on the screened Coulomb interaction. The intuitive nature of screening—the simple idea that charges rearrange themselves, or respond, to an added charge—is the major reason behind the appeal and success of the GW approximation. The frequency dependence of the screened Coulomb interaction is largely what sets GW apart from density functional or Hartree-Fock theories.

Impressive advances in code development and computing resources have pushed GW calculations to a new scale. At the computationally lowest level of theory, G0W0 calculations remain the most widely used and can be routinely applied to systems with hundreds of atoms. Within the G0W0 approach, we have outlined the practical considerations presented to the user before performing any calculation: starting point, basis set, evaluation of the self-energy, and convergence are all for the user to decide. For a broad class of systems, G0W0 already gives electron addition and removal energies in good agreement with experiment. These successes give GW an impressive ranking in computational value, or accuracy for computational cost, on any list of first-principles electronic structure methods. The versatility of GW assures that new applications in physics, chemistry, and materials science will continue to emerge in the future.

Author Contributions

All authors listed have made a substantial contribution to the work. DG focused on numerical aspects, MD on theoretical concepts. PR contributed to all parts.


This work is supported by the Academy of Finland through grant Nos. 284621, 305632, 316347, and 316168.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.


The authors acknowledge the CSC-IT Center for Science, Finland, for generous computational resources and the Aalto University School of Science Science-IT project for computational resources. We would like to thank Fabio Caruso, Jan Wilhelm, Leeor Kronik, and Lucia Reining for fruitful discussions and valuable feedback.

Appendix A: Integral Notation

We adopt the following integral notation

ϕi|Ô|ϕj= drdrϕi*(r)Ô(r,r)ϕj(r)    (109)
ϕiϕj|Ô|ϕkϕl=drdrϕi*(r)ϕj*(r)Ô(r,r)  ×ϕk(r)ϕl(r)    (110)

where Ô(r, r) is an operator that depends on the spatial variables r = (x, y, z) and r = (x′, y′, z′). Furthermore the following notation for Coulomb integrals are used

ϕiϕj|ϕkϕl=ij|kl =drdrϕi*(r)ϕj*(r)v(r,r)ϕk(r)ϕl(r),    (111)

where v(r, r) = 1/|rr| is the Coulomb operator. The antisymmetrized Coulomb integrals are defined as

ϕiϕj||ϕkϕl=ij||kl =ϕiϕj|ϕkϕlϕiϕj|ϕlϕk.    (112)

Appendix B: The Many-Body Problem

In first principles electronic structure theory, the aim is to solve the Schrödinger equation. For simplicity we focus on the non-relativistic time-independent Schrödinger equation. For a system of N electrons and M nuclei, the Schrödinger equation is given by

ĤΨ=EΨ    (B1)

with the many-body Hamiltonian

H^=i=1Ni22+12ijN1|rirj|i=1Na=1MZa|riRa|H^elec  a=1Ma22Ma+a=1Mb>aMZaZb|RaRb|.    (B2)

ri and Ra are the positions of the electrons and the nuclei, respectively, and Za is the charge of the nuclei.

To make this system of coupled electrons and nuclei more tractable, the Born-Oppenheimer (BO) approximation of clamped nuclei is frequently introduced. In this case, we only need to consider the electronic Hamiltonian by itself,

ĤelecΨelec=EelecΨelec,    (B3)

in which the many-electron wave function Ψelec depends parametrically on the position of the nuclei. The electronic Hamiltonian is then given by

Ĥelec=i=1N[i22+vext(ri)]+12ijN1|rirj|    (B4)
=i=1Nĥ0(ri)+12ijNv(ri,rj).    (B5)

We use v to denote the bare Coulomb interaction and the external potential is the same for every electron

vext(r)=a=1MZa|rRa|.    (B6)

The kinetic energy and external potential are grouped together in ĥ0(ri) in Equation (B5). Approaches to solve the electronic Schrödinger equation (B3) can be largely grouped according to the basic variable they work with: the many-body wave function, the density-matrix, the density, or Green's functions. Each choice has its advantages and disadvantages and no consensus has been reached on the optimal choice. Green's functions have a natural connection, however, to the particle addition/removal problem and theoretical spectroscopy.

B.1. Green's Function Formalism

To derive equations for the single-particle Green's function that are more amenable to approximations than definitions in Equations (6) and (7), we start from the equation of motion for the field operators, which relates the time derivative of ψ^ to the commutator of ψ^ and Ĥelec (Fetter and Walecka, 1971; Gross et al., 1991):

itψ^(x,t)=[ψ^(x,t),Ĥelec],    (B7)

where x contains also the spin variable σ, i.e., x = (r, σ). Evaluating the commutator in the Heisenberg picture and applying the anti-commutation relations yields for the equation of motion

itψ^(x,t)=              [ĥ0(r)+ψ^(x,t)v(r,r)ψ^(x,t)dx]ψ^(x,t).    (B8)

The equation of motion for the Green's function (6) then follows from (B8)

[itĥ0(r)]G(xt,xt)=                  δ(tt)δ(xx)idxv(r,r)×                  N|T^[ψ^(x,t)ψ^(x,t)ψ^(x,t)ψ^(x,t)]|N    (B9)

The term under the integral contains the two-particle Green's function, G2(xt,xt,xt+,xt), and includes all two-body correlations in the system. In order to calculate the one-particle Green's function we would therefore require the equation of motion for the two-particle Green's function, which in turn introduces the three-particle Green's function. Applied iteratively, this procedure creates an infinite series of higher order Green's functions and thus describes all the possible many-body interactions in the system. In practice, however, the resulting recurrence relation for the nth order Green's function is impossible to solve for large n. Instead we introduce the non-local, time-dependent self-energy Σ̄(xt,xt)

idxv(r,r)G2(xt,xt,xt+,xt)                  dtdxΣ̄(xt,xt)G(xt,xt).    (B10)

Analogous to other electronic structure methods, we separate out the most dominant term, the Hartree potential

vH(r)=drv(r,r)N|ψ^(r,t)ψ^(r,t)|N    (B11)
=drv(r,r)n(r),    (B12)

where n(r) is the electron density. With this definition, the equation of motion for the Green's function (B9) adopts the much more convenient form of an integral equation involving the self-energy, Σ=Σ̄-vH:

[itĥ0(r)vH(r)]G(xt,xt)=δ(tt)δ(rr)                       +dtdxΣ(xt,xt)G(xt,xt).    (B13)

If G0 now denotes the Green's function that is a solution to ĥ=ĥ0+vH(r) (the kinetic energy, the external potential, and the Hartree potential), then Equation (B13) can be further rewritten as

G(1,2)=G0(1,2) +G0(1,3)Σ(3,4)G(4,2)d(3,4).    (B14)

Here we changed to the abbreviated notation (1, 2, …) for the set of position, time and spin variables (x1t1, x2t2, …). Accordingly ∫ d(1) is a shorthand notation for the integration in all three variables of the corresponding triple(s). In this context 1+ implies the addition of a positive infinitesimal to the time argument 1. Equation (B14) is again Dyson's equation (Dyson, 1949a,b) (see also Equations (13) and (82)) and links the non-interacting Green's function, G0, to the interacting one, G. Dyson's equation gives a physical interpretation to the self-energy instead of simply its definition by Equation (B10). The self-energy quantifies the difference between a bare and a fully interacting electron, or quasielectron. This brings us back to our phenomenological consideration in section 2. An additional electron or hole drags a dynamically adjusting polarization cloud through the system that alters its energy. Hence the name self-energy. If instead G0 is the solution to the mean-field Hamiltonian ĥMF [e.g. of Kohn-Sham density-functional theory (Kohn and Sham, 1965)] then the self-energy, Σ, embodies the difference between a quasielectron and an electron in the static mean-field.

At this stage in the derivation, the self-energy as well as the Green's function are still exact and contain the electron-electron interaction to all orders. A full solution of Equations (B13) and (B14) is not tractable and approximations are required.

B.2 Hedin's Equations

In 1965 Hedin expanded the Green's function and the self-energy in terms of the screened instead of the bare Coulomb interaction (Hedin, 1965). Introducing a small perturbing field φ that will later be set to zero, the operator identity due to Schwinger (1951)

G2(1,3,2,3+)=G(1,2)G(3,3+)δG(1,2)δφ(3)    (B15)

can be used to eliminate the two particle Green's function from Equation (B10) for the self-energy. Multiplying the resulting equation with G−1 leads to the following expression for the self-energy

Σ¯(1,2)=δ(1,2) d(3)v(1,3)G(3,3+)  i d(3,4)v(1,3)G(1,4)δG1(4,2)δφ(3) =δ(1,2)vH(1)  i d(3,4)v(1,3)G(1,4)δG1(4,2)δφ(3),    (B16)

where Σ̄ is defined to include the Hartree potential vH, unlike Σ. Using the following definitions

total potential:

V(1)=φ(1)+vH(1)    (B17)

3-point vertex:

Γ(1,2,3)=δG1(1,2)δV(3)    (B18)

dielectric function:

ε1(1,2)=δV(1)δφ(2)    (B19)

screened Coulomb interaction:

W(1,2)=d(3)ε1(1,3)v(3,2)    (B20)

irreducible polarizability:

P(1,2)=iδG(1,1+)δV(2)=δn(1)δV(2)    (B21)

reducible polarizability:

χ(1,2)=iδG(1,1+)δφ(2)=δn(1)δφ(2)    (B22)

we arrive at Hedin's equations (Hedin, 1965)

P(1,2)=id(3,4)G(4,2)G(2,3)Γ(3,4,1)    (B23)
W(1,2)=v(1,2)+d(3,4)v(1,3)P(3,4)W(4,2)    (B24)
Σ(1,2)=id(3,4)G(1,4)W(1+,3)Γ(4,2,3)    (B25)
Γ(1,2,3)=δ(1,2)δ(1,3)  +d(4,5,6,7)δΣ(1,2)δG(4,5)G(4,6)G(7,5)Γ(6,7,3).    (B26)

Dyson's equation (B14) closes this set of integro-differential equations, which is shown pictorially in Figure 24A.

The benefit of Hedin's equations is that the self-energy is given in terms of an effective, or screened, rather than the bare Coulomb interaction

W(1,2)=d(3)ε1(1,3)v(3,2).    (B27)

The screening

ε1(1,2)=δ(1,2)+d(3) v(1,3)χ(3,2)    (B28)

follows from the reducible polarizability or density-density response function

χ(1,2)=χ0(1,2)+d(3,4)χ0(1,3)v(3,4)χ(4,2)               =|T^[δn^(1)δn^(2)]|N .    (B29)

Here δn^(1) is a density fluctuation δn^(1)=n^(1)-n(1) of the density around its average value, where n^(1)=iNδ(r1-ri). For a polarizable medium, like a solid, screening is large and the screened Coulomb interaction will differ considerably from the bare one. It therefore makes sense to build a perturbation series on W instead of v. Hedin's equations are physically transparent in this sense. Electron-hole pairs are created in the polarizability (one Green's function for the electron, a separate Green's function for the hole). They interact through the vertex function, which is determined by the change in potential upon excitation. The polarizability, in turn, determines the dielectric function, which screens the Coulomb interaction. The self-energy quantifies the energy contribution that the added electron or hole experiences through the interaction with its surrounding.

Appendix C: Computational Details

Figures 11, 12, 14, 16B, 18, 25, and 35 present original content. All calculations are performed with the FHI-aims program package (Blum et al., 2009; Havu et al., 2009), which expands the all-electron KS equations in numeric-atom-centered orbitals (NAOs), see section 4.4. The structures have been optimized at the DFT level using the Perdew-Burke-Ernzerhof (PBE) functional (Perdew et al., 1996a) to model the XC potential and NAOs of tier 2 quality (Blum et al., 2009) to represent core and valence electrons. Dispersion corrections are accounted for by employing the Tkatchenko-Scheffler van der Waals correction (Tkatchenko and Scheffler, 2009).

G0W0 calculations have been performed with the contour-deformation (CD) technique (Golze et al., 2018) if not indicated otherwise. Calculations with the analytic continuation (AC) (Ren et al., 2012a) have been conducted for the G0W0 self-energies in Figure 14 and for the scGW result in Figure 18. For both methods, CD and AC, a modified Gauss-Legendre grid with 200 grid points is used for the numerical integration of the integral over the imaginary frequency axis. In case of the AC approach, the same set of grid points {} is employed to calculate Σsc(iω), which has been fitted to a Padé approximant with at least 16 parameters (Vidberg and Serene, 1977), unless stated otherwise.

Quasiparticle energies have been computed by iteratively solving Equation (22). Furthermore, all quasiparticle energies have been extrapolated to the complete basis set limit using the Dunning basis set family cc-pVnZ (n = 3 − 6) (Dunning, 1989; Wilson et al., 1996). The cc-pVnZ basis sets are all-electron Gaussian basis sets, which can be considered as a special case of an NAO and are treated numerically in FHI-aims. The extrapolation has been performed by a linear regression against the inverse of the total number of basis functions.


1. ^The theory directly comparable to DFT is the Luttinger-Ward formalism, not the GW approximation. Here, we focus on the chronological development of practical electronic structure calculations with Kohn-Sham DFT or GW.

2. ^Throughout this article the energy zero is chosen to be the top of the valence bands for extended systems and the vacuum level for finite systems.

3. ^Atomic units 4πϵ0 = h = e = me = 1, where e and me are the charge and mass of an electron, respectively, will be used in the remainder of this article.

4. ^We consider only the zero temperature G and assume μ = EF.

5. ^For the remainder of the review, η is always assumed to be a positive infinitesimal.

6. ^1+ means the time t1 evaluated at an infinitesimally later time. Such time infinitesimals appear in order to define the time-ordering for quantities meant to be evaluated in the instantaneous limit. The ground state density, for example, is time independent but must be written as n(1) = −iG(1, 1+) so that the time-ordering makes sense.

7. ^For simplicity, we will drop spin variables in the following parts of section 4.

8. ^While the imaginary part might be small, it is nonetheless important as its inverse is proportional to the lifetime of the state.

9. ^It has also been pointed out that the exact W from RPA does not satisfy the f-sum rule.

10. ^Note that different sign conventions are used for EA in literature. We define EA as the energy required to detach an electron from a negatively charged species. If EA is defined as the energy required to add an electron to a neutral atom, the sign swaps. EA refers always to the LUMO. The label for the state is therefore dropped in the following.


Adamo, C., and Barone, V. (1999). Towards reliable density functional methods without adjustable parameters: the PBE0 model. J. Chem. Phys 110:6158. doi: 10.1063/1.478522

CrossRef Full Text | Google Scholar

Adler, S. L. (1962). Quantum theory of the dielectric constant in real solids. Phys. Rev. 126:413. doi: 10.1103/PhysRev.126.413

CrossRef Full Text | Google Scholar

Adragna, G., Del Sole, R., and Marini, A. (2003). Ab initio calculation of the exchange-correlation kernel in extended systems. Phys. Rev. B 68:165108. doi: 10.1103/PhysRevB.68.165108

CrossRef Full Text | Google Scholar

Aguilera, I., Friedrich, C., Bihlmayer, G., and Blügel, S. (2013a). GW study of topological insulators Bi2Se3, Bi2Te3, and Sb2Te3: beyond the perturbative one-shot approach. Phys. Rev. B 88:045206. doi: 10.1103/PhysRevB.88.045206

CrossRef Full Text | Google Scholar

Aguilera, I., Friedrich, C., and Blügel, S. (2013b). Spin-orbit coupling in quasiparticle studies of topological insulators. Phys. Rev. B 88:165136. doi: 10.1103/PhysRevB.88.165136

CrossRef Full Text | Google Scholar

Aguilera, I., Friedrich, C., and Blügel, S. (2015a). Electronic phase transitions of bismuth under strain from relativistic self-consistent GW calculations. Phys. Rev. B 91:125129. doi: 10.1103/PhysRevB.91.125129

CrossRef Full Text | Google Scholar

Aguilera, I., Nechaev, I. A., Friedrich, C., Blügel, S., and Chulkov, E. V. (2015b). “Chapter 7: Many-Body Effects in the Electronic Structure of Topological Insulators,” in Topol. Insulat. eds F. Ortmann, S. Roche, and S. O. Valenzuela (John Wiley & Sons, Ltd.), 161–189.

Google Scholar

Aguilera, I., Vidal, J., Wahnón, P., Reining, L., and Botti, S. (2011). First-principles study of the band structure and optical absorption of CuGaS2. Phys. Rev. B 84:085145. doi: 10.1103/PhysRevB.84.085145

CrossRef Full Text | Google Scholar

Ahmed, T., Albers, R. C., Balatsky, A. V., Friedrich, C., and Zhu, J.-X. (2014). GW quasiparticle calculations with spin-orbit coupling for the light actinides. Phys. Rev. B 89:035104. doi: 10.1103/PhysRevB.89.035104

CrossRef Full Text | Google Scholar

Albrecht, S., Onida, G., and Reining, L. (1997). Ab initio calculation of the quasiparticle spectrum and excitonic effects in Li2O. Phys. Rev. B 55, 10278–10281. doi: 10.1103/PhysRevB.55.10278

CrossRef Full Text | Google Scholar

Albrecht, S., Reining, L., Sole, R. D., and Onida, G. (1998a). Ab Initio calculation of excitonic effects in the optical spectra of semiconductors. Phys. Rev. Lett. 80, 4510–4513. doi: 10.1103/PhysRevLett.80.4510

CrossRef Full Text | Google Scholar

Albrecht, S., Reining, L., Sole, R. D., and Onida, G. (1998b). Excitonic effects in the optical properties. Phys. Status Solidi A 170, 189–197. doi: 10.1002/(SICI)1521-396X(199812)170:2<189::AID-PSSA189>3.0.CO;2-3

CrossRef Full Text | Google Scholar

Almbladh, C. O., and Hedin, L. (1983). “Beyond the one-electron model. Many-body effects in atoms, molecules, and solids,” in Handbook on Synchroton Radiation, Vol. 1, ed E. E. Koch (Amsterdam: North-Holland), 607–904.

Google Scholar

Almbladh, C. O., and von Barth, U. (1985). Exact results for the charge and spin densities, exchange-correlation potentials, and density-functional eigenvalues. Phys. Rev. B 31:3231. doi: 10.1103/PhysRevB.31.3231

PubMed Abstract | CrossRef Full Text | Google Scholar

Alves-Santos, M., Jorge, L. M. M., Caldas, M. J., and Varsano, D. (2014). Electronic structure of interfaces between thiophene and TiO2 nanostructures. J. Phys. Chem. C 118, 13539–13544. doi: 10.1021/jp407275e

CrossRef Full Text | Google Scholar

Amadon, B., Biermann, S., Georges, A., and Aryasetiawan, F. (2006). The α−γ transition of cerium is entropy driven. Phys. Rev. Lett. 96:066402. doi: 10.1103/PhysRevLett.96.066402

PubMed Abstract | CrossRef Full Text | Google Scholar

Andersen, O. K. (1975). Linear methods in band theory. Phys. Rev. B 12, 3060–3083. doi: 10.1103/PhysRevB.12.3060

CrossRef Full Text | Google Scholar

Antonius, G., Poncé, S., Boulanger, P., and Côté, M. X. (2014). Many-body effects on the zero-point renormalization of the band structure. Phys. Rev. Lett. 112:215501. doi: 10.1103/PhysRevLett.112.215501

CrossRef Full Text | Google Scholar

Appelhans, D. J., Carr, L. D., and Lusk, M. T. (2010a). Embedded ribbons of graphene allotropes: an extended defect perspective. New J. Phys. 12:125006. doi: 10.1088/1367-2630/12/12/125006

CrossRef Full Text | Google Scholar

Appelhans, D. J., Lin, Z., and Lusk, M. T. (2010b). Two-dimensional carbon semiconductor: density functional theory calculations. Phys. Rev. B 82:073410. doi: 10.1103/PhysRevB.82.073410

CrossRef Full Text | Google Scholar

Arora, A., Drüppel, M., Schmidt, R., Deilmann, T., Schneider, R., Molas, M. R., et al. (2017). Interlayer excitons in a bulk van der Waals semiconductor. Nat. Commun. 8:639. doi: 10.1038/s41467-017-00691-5

CrossRef Full Text | Google Scholar

Aryasetiawan, F. (1992). Self-energy of ferromagnetic nickel in the GW approximation. Phys. Rev. B 46, 13051–13064. doi: 10.1103/PhysRevB.46.13051

PubMed Abstract | CrossRef Full Text | Google Scholar

Aryasetiawan, F. (1997). GW method for 3d and 4f systems. Physica B 237–238, 321–323. doi: 10.1016/S0921-4526(97)00196-8

CrossRef Full Text | Google Scholar

Aryasetiawan, F., and Biermann, S. (2008). Generalized Hedin's equations for quantum many-body systems with spin-dependent interactions. Phys. Rev. Lett. 100:116402. doi: 10.1103/PhysRevLett.100.116402

PubMed Abstract | CrossRef Full Text | Google Scholar

Aryasetiawan, F., and Biermann, S. (2009). Generalized Hedin equations and GW approximation for quantum many-body systems with spin-dependent interactions. J. Phys. Condens. Matter 21:064232. doi: 10.1088/0953-8984/21/6/064232

CrossRef Full Text | Google Scholar

Aryasetiawan, F., and Gunnarsson, O. (1995). Electronic structure of NiO in the GW approximation. Phys. Rev. Lett. 74:3221. doi: 10.1103/PhysRevLett.74.3221

PubMed Abstract | CrossRef Full Text | Google Scholar

Aryasetiawan, F., and Gunnarsson, O. (1998). The GW method. Rep. Prog. Phys. 61:237. doi: 10.1088/0034-4885/61/3/002

CrossRef Full Text | Google Scholar

Aryasetiawan, F., Hedin, L., and Karlsson, K. (1996). Multiple plasmon satellites in Na and Al spectral functions from ab initio cumulant expansion. Phys. Rev. Lett. 77, 2268–2271. doi: 10.1103/PhysRevLett.77.2268

PubMed Abstract | CrossRef Full Text | Google Scholar

Aryasetiawan, F., and Karlsson, K. (1996). GW spectral functions of Gd and NiO. Phys. Rev. B 54:5353. doi: 10.1103/PhysRevB.54.5353

PubMed Abstract | CrossRef Full Text | Google Scholar

Atalla, V., Yoon, M., Caruso, F., Rinke, P., and Scheffler, M. (2013). Hybrid density functional theory meets quasiparticle calculations: a consistent electronic structure approach. Phys. Rev. B 88:165122. doi: 10.1103/PhysRevB.88.165122

CrossRef Full Text | Google Scholar

Aulbur, W. G., Jönsson, L., and Wilkins, J. W. (2000). “Quasiparticle calculations in solids,” in Solid State Physics, Vol. 54, eds H. Ehrenreich and F. Spaepen (San Diego, CA: Academic Press), 1–218.

Bacelar, M. R., Schöne, W.-D., Keyling, R., and Ekardt, W. (2002). Lifetime of excited electrons in transition metals. Phys. Rev. B 66:153101. doi: 10.1103/PhysRevB.66.153101

CrossRef Full Text | Google Scholar

Bachelet, G. B., Hamann, D. R., and Schlüter, M. (1982). Pseudopotentials that work: from H to Pu. Phys. Rev. B 26, 4199–4228. doi: 10.1103/PhysRevB.26.4199

CrossRef Full Text | Google Scholar

Baerends, E. J., Gritsenko, O. V., and van Meer, R. (2013). The Kohn Sham gap, the fundamental gap and the optical gap: the physical meaning of occupied and virtual Kohn Sham orbital energies. Phys. Chem. Chem. Phys. 15, 16408–16425. doi: 10.1039/c3cp52547c

PubMed Abstract | CrossRef Full Text | Google Scholar

Bagus, P. S., Illas, F., Pacchioni, G., and Parmigiani, F. (1999). Mechanisms responsible for chemical shifts of core-level binding energies and their relationship to chemical bonding. J. Electron Spectrosc. Relat. Phenom. 100, 215–236. doi: 10.1016/S0368-2048(99)00048-1

CrossRef Full Text | Google Scholar

Bagus, P. S., Ilton, E. S., and Nelin, C. J. (2013). The interpretation of XPS spectra: insights into materials properties. Surf. Sci. Rep. 68, 273–304. doi: 10.1016/j.surfrep.2013.03.001

CrossRef Full Text | Google Scholar

Baroni, S., Giannozzi, P., and Testa, A. (1987). Green's-function approach to linear response in solids. Phys. Rev. Lett. 58, 1861–1864. doi: 10.1103/PhysRevLett.58.1861

PubMed Abstract | CrossRef Full Text | Google Scholar

Baroni, S., Gironcoli, S. D., Corso, A. D., and Giannozzi, P. (2001). Phonons and related crystal properties from density-functional perturbation theory. Rev. Mod. Phys. 73:515. doi: 10.1103/RevModPhys.73.515

CrossRef Full Text | Google Scholar

Bechstedt, F. (2014). Many-Body Approach to Electronic Excitations: Concepts and Applications. Springer Series in Solid-State Sciences. Berlin; Heidelberg: Springer.

Google Scholar

Bechstedt, F., and Furthmüller, J. (2002). Do we know the fundamental energy gap of InN? J. Cryst. Growth 246:315. doi: 10.1016/S0022-0248(02)01756-6

CrossRef Full Text | Google Scholar

Becke, A. (1988). Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A 38:3098. doi: 10.1103/PhysRevA.38.3098

PubMed Abstract | CrossRef Full Text | Google Scholar

Benedict, L. X., and Shirley, E. L. (1999). Ab initio calculation of ϵ2(ω) including the electron-hole interaction: application to GaN and CaF2. Phys. Rev. B 59, 5441–5451. doi: 10.1103/PhysRevB.59.5441

CrossRef Full Text | Google Scholar

Benedict, L. X., Shirley, E. L., and Bohn, R. B. (1998a). Optical absorption of insulators and the electron-hole interaction: an ab initio calculation. Phys. Rev. Lett. 80:4514–4517. doi: 10.1103/PhysRevLett.80.4514

CrossRef Full Text | Google Scholar

Benedict, L. X., Shirley, E. L., and Bohn, R. B. (1998b). Theory of optical absorption in diamond, Si, Ge, and GaAs. Phys. Rev. B 57, R9385–R9387. doi: 10.1103/PhysRevB.57.R9385

CrossRef Full Text | Google Scholar

Berger, J. A., Reining, L., and Sottile, F. (2010). Ab initio calculations of electronic excitations: collapsing spectral sums. Phys. Rev. B 82:041103. doi: 10.1103/PhysRevB.82.041103

CrossRef Full Text | Google Scholar

Berger, J. A., Romaniello, P., Tandetzky, F., Mendoza, B. S., Brouder, C., and Reining, L. (2014). Solution to the many-body problem in one point. New J. Phys. 16:113025. doi: 10.1088/1367-2630/16/11/113025

CrossRef Full Text | Google Scholar

Bernardi, M., Palummo, M., and Grossman, J. C. (2012). Optoelectronic properties in monolayers of hybridized graphene and hexagonal boron nitride. Phys. Rev. Lett. 108:226805. doi: 10.1103/PhysRevLett.108.226805

PubMed Abstract | CrossRef Full Text | Google Scholar

Berseneva, N., Gulans, A., Krasheninnikov, A. V., and Nieminen, R. M. (2013). Electronic structure of boron nitride sheets doped with carbon from first-principles calculations. Phys. Rev. B 87:035404. doi: 10.1103/PhysRevB.87.035404

CrossRef Full Text | Google Scholar

Betzinger, M., Friedrich, C., and Blügel, S. (2013). Precise response functions in all-electron methods: Generalization to nonspherical perturbations and application to NiO. Phys. Rev. B 88:075130. doi: 10.1103/PhysRevB.88.075130

CrossRef Full Text | Google Scholar

Betzinger, M., Friedrich, C., Görling, A., and Blügel, S. (2012). Precise response functions in all-electron methods: application to the optimized-effective-potential approach. Phys. Rev. B 85:245124. doi: 10.1103/PhysRevB.85.245124

CrossRef Full Text | Google Scholar

Betzinger, M., Friedrich, C., Görling, A., and Blügel, S. (2015). Precise all-electron dynamical response functions: application to COHSEX and the RPA correlation energy. Phys. Rev. B 92:245101. doi: 10.1103/PhysRevB.92.245101

CrossRef Full Text | Google Scholar

Bhimanapati, G. R., Lin, Z., Meunier, V., Jung, Y., Cha, J., Das, S., et al. (2015). Recent advances in two-dimensional materials beyond graphene. ACS Nano 9, 11509–11539. doi: 10.1021/acsnano.5b05556

PubMed Abstract | CrossRef Full Text | Google Scholar

Bieder, J., and Amadon, B. (2014). Thermodynamics of the α-γ transition in cerium from first principles. Phys. Rev. B 89:195132. doi: 10.1103/PhysRevB.89.195132

CrossRef Full Text | Google Scholar

Biermann, S. (2014). Dynamical screening effects in correlated electron materials a progress report on combined many-body perturbation and dynamical mean field theory: GW + DMFT. J. Phys. Condens. Matter 26:173202. doi: 10.1088/0953-8984/26/17/173202

PubMed Abstract | CrossRef Full Text | Google Scholar

Biermann, S., Aryasetiawan, F., and Georges, A. (2003). First-principles approach to the electronic structure of strongly correlated systems: combining the GW approximation and dynamical mean-field theory. Phys. Rev. Lett. 90:086402. doi: 10.1103/PhysRevLett.90.086402

PubMed Abstract | CrossRef Full Text | Google Scholar

Blase, X., Attaccalite, C., and Olevano, V. (2011). First-principles GW calculations for fullerenes, porphyrins, phtalocyanine, and other molecules of interest for organic photovoltaic applications. Phys. Rev. B 83:115103. doi: 10.1103/PhysRevB.83.115103

CrossRef Full Text | Google Scholar

Blase, X., Rubio, A., Louie, S. G., and Cohen, M. L. (1995). Quasiparticle band structure of bulk hexagonal boron nitride and related systems. Phys. Rev. B 51, 6868–6875. doi: 10.1103/PhysRevB.51.6868

PubMed Abstract | CrossRef Full Text | Google Scholar

Bloch, F. (1929). Über die quantenmechanik der elektronen in kristallgittern. Z. Phys. 52, 555–600. doi: 10.1007/BF01339455

CrossRef Full Text | Google Scholar

Blöchl, P. E. (1994). Projector augmented-wave method. Phys. Rev. B 50, 17953–17979. doi: 10.1103/PhysRevB.50.17953

PubMed Abstract | CrossRef Full Text | Google Scholar

Blum, V., Gehrke, R., Hanke, F., Havu, P., Havu, V., Ren, X., et al. (2009). Ab initio molecular simulations with numeric atom-centered orbitals. Comput. Phys. Commun. 180, 2175–2196. doi: 10.1016/j.cpc.2009.06.022

CrossRef Full Text | Google Scholar

Bobbert, P. A., and van Haeringen, W. (1994). Lowest-order vertex-correction contribution to the direct gap of silicon. Phys. Rev. B 49, 10326–10331. doi: 10.1103/PhysRevB.49.10326

PubMed Abstract | CrossRef Full Text | Google Scholar

Bockstedte, M., Marini, A., Pankratov, O., and Rubio, A. (2010). Many-body effects in the excitation spectrum of a defect in SiC. Phys. Rev. Lett. 105:026401. doi: 10.1103/PhysRevLett.105.026401

PubMed Abstract | CrossRef Full Text | Google Scholar

Bois, J., and Körzdörfer, T. (2017). Size-dependence of nonempirically tuned DFT starting points for G0W0 applied to π-conjugated molecular chains. J. Chem. Theor. Comput. 13, 4962–4971. doi: 10.1021/acs.jctc.7b00557

PubMed Abstract | CrossRef Full Text | Google Scholar

Botti, S., and Marques, M. A. L. (2013). Strong renormalization of the electronic band gap due to lattice polarization in the GW formalism. Phys. Rev. Lett. 110:226404. doi: 10.1103/PhysRevLett.110.226404

PubMed Abstract | CrossRef Full Text | Google Scholar

Bouhassoune, M., and Schindlmayr, A. (2010). Electronic structure and effective masses in strained silicon. Phys. Status Solidi C 7, 460–463. doi: 10.1002/pssc.200982470

CrossRef Full Text | Google Scholar

Bredas, J.-L. (2014). Mind the gap! Mater. Horiz. 1, 17–19. doi: 10.1039/C3MH00098B

CrossRef Full Text | Google Scholar

Brivio, F., Butler, K. T., Walsh, A., and van Schilfgaarde, M. (2014). Relativistic quasiparticle self-consistent electronic structure of hybrid halide perovskite photovoltaic absorbers. Phys. Rev. B 89:155204. doi: 10.1103/PhysRevB.89.155204

CrossRef Full Text | Google Scholar

Bruneval, F. (2009). GW Approximation of the many-body problem and changes in the particle number. Phys. Rev. Lett. 103:176403. doi: 10.1103/PhysRevLett.103.176403

PubMed Abstract | CrossRef Full Text | Google Scholar

Bruneval, F. (2012). Ionization energy of atoms obtained from GW self-energy or from random phase approximation total energies. J. Chem. Phys. 136:194107. doi: 10.1063/1.4718428

PubMed Abstract | CrossRef Full Text | Google Scholar

Bruneval, F. (2016). Optimized virtual orbital subspace for faster GW calculations in localized basis. J. Chem. Phys. 145:234110. doi: 10.1063/1.4972003

PubMed Abstract | CrossRef Full Text | Google Scholar

Bruneval, F., and Gatti, M. (2014). “Quasiparticle self-consistent GW method for the spectral properties of complex materials,” in First Principles Approaches to Spectroscopic Properties of Complex Materials, eds C. D. Valentin, S. Botti, and M. Cococcioni (Berlin; Heidelberg: Springer), 99–135.

PubMed Abstract | Google Scholar

Bruneval, F., and Gonze, X. (2008). Accurate GW self-energies in a plane-wave basis using only a few empty states: towards large systems. Phys. Rev. B 78:085125. doi: 10.1103/PhysRevB.78.085125

CrossRef Full Text | Google Scholar

Bruneval, F., Hamed, S. M., and Neaton, J. B. (2015). A systematic benchmark of the ab initio Bethe-Salpeter equation approach for low-lying optical excitations of small organic molecules. J. Chem. Phys. 142:244101. doi: 10.1063/1.4922489

PubMed Abstract | CrossRef Full Text | Google Scholar

Bruneval, F., and Marques, M. A. L. (2013). Benchmarking the starting points of the GW approximation for molecules. J. Chem. Theor. Comput. 9, 324–329. doi: 10.1021/ct300835h

PubMed Abstract | CrossRef Full Text | Google Scholar

Bruneval, F., Rangel, T., Hamed, S. M., Shao, M., Yang, C., and Neaton, J. B. (2016). MOLGW1: many-body perturbation theory software for atoms, molecules, and clusters. Comput. Phys. Commun. 208, 149–161. doi: 10.1016/j.cpc.2016.06.019

CrossRef Full Text | Google Scholar

Bruneval, F., Sottile, F., Olevano, V., Del Sole, R., and Reining, L. (2005). Many-body perturbation theory using the density-functional concept: beyond the GW approximation. Phys. Rev. Lett. 94:186402. doi: 10.1103/PhysRevLett.94.186402

PubMed Abstract | CrossRef Full Text | Google Scholar

Bruneval, F., Vast, N., Reining, L., Izquierdo, M., Sirotti, F., and Barrett, N. (2006). Exchange and correlation effects in electronic excitations of Cu20. Phys. Rev. Lett. 97:267601. doi: 10.1103/PhysRevLett.97.267601

CrossRef Full Text | Google Scholar

Butcher, K. S. A., and Tansley, T. L. (2005). InN, latest development and a review of the band-gap controversy. Superlatt. Microstruct. 38:1. doi: 10.1016/j.spmi.2005.03.004

CrossRef Full Text | Google Scholar

Cannuccia, E., and Marini, A. (2011). Effect of the quantum zero-point atomic motion on the optical and electronic properties of diamond and trans-polyacetylene. Phys. Rev. Lett. 107:255501. doi: 10.1103/PhysRevLett.107.255501

PubMed Abstract | CrossRef Full Text | Google Scholar

Caruso, F., Atalla, V., Ren, X., Rubio, A., Scheffler, M., and Rinke, P. (2014). First-principles description of charge transfer in donor-acceptor compounds from self-consistent many-body perturbation theory. Phys. Rev. B 90:085141. doi: 10.1103/PhysRevB.90.085141

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Caruso, F., and Giustino, F. (2015). Spectral fingerprints of electron-plasmon coupling. Phys. Rev. B 92:045123. doi: 10.1103/PhysRevB.92.045123

CrossRef Full Text | Google Scholar

Caruso, F., and Giustino, F. (2016). The GW plus cumulant method and plasmonic polarons: application to the homogeneous electron gas. Eur. Phys. J. B 89:238. doi: 10.1140/epjb/e2016-70028-4

CrossRef Full Text | Google Scholar

Caruso, F., Lambert, H., and Giustino, F. (2015). Band structures of plasmonic polarons. Phys. Rev. Lett. 114:146404. doi: 10.1103/PhysRevLett.114.146404

PubMed Abstract | CrossRef Full Text | Google Scholar

Caruso, F., Rinke, P., Ren, X., Rubio, A., and Scheffler, M. (2013a). Self-consistent GW: all-electron implementation with localized basis functions. Phys. Rev. B 88:075105. doi: 10.1103/PhysRevB.88.075105

CrossRef Full Text | Google Scholar

Caruso, F., Rinke, P., Ren, X., Scheffler, M., and Rubio, A. (2012a). Unified description of ground and excited states of finite systems: the self-consistent GW approach. Phys. Rev. B 86:081102(R). doi: 10.1103/PhysRevB.86.081102

CrossRef Full Text | Google Scholar

Caruso, F., Rinke, P., Ren, X., Scheffler, M., and Rubio, A. (2012b). Unpublished result from calculations for Ref. PRB 86, 081102. CCSD dipole moment calculated at the aug-cc-pVTZ level using the settings described in SI of PRB 86, 081102.

Caruso, F., Rohr, D. R., Hellgren, M., Ren, X., Rinke, P., Rubio, A., et al. (2013b). Bond breaking and bond formation: how electron correlation is captured in many-body perturbation theory and density-functional theory. Phys. Rev. Lett. 110:146403. doi: 10.1103/PhysRevLett.110.146403

PubMed Abstract | CrossRef Full Text | Google Scholar

Casadei, M., Ren, X., Rinke, P., Rubio, A., and Scheffler, M. (2012). Density-functional theory for f-electron systems: the α-γ phase transition in cerium. Phys. Rev. Lett. 109:146402. doi: 10.1103/PhysRevLett.109.146402

PubMed Abstract | CrossRef Full Text | Google Scholar

Casadei, M., Ren, X., Rinke, P., Rubio, A., and Scheffler, M. (2016). Density functional theory study of the α-γ phase transition in cerium: Role of electron correlation and f-orbital localization. Phys. Rev. B 93:075153. doi: 10.1103/PhysRevB.93.075153

CrossRef Full Text | Google Scholar

Casida, M. E. (1995a). Generalization of the optimized-effective-potential model to include electron correlation: a variational derivation of the Sham-Schlüter equation for the exact exchange-correlation potential. Phys. Rev. A 51:2005. doi: 10.1103/PhysRevA.51.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Casida, M. E. (1995b). “Time-dependent density functional response theory for molecules,” in Recent Advances in Density Functional Methods, ed D. P. Chong (Singapore: World Scientific), 155–192.

Google Scholar

Cazzaniga, M. (2012). GW and beyond approaches to quasiparticle properties in metals. Phys. Rev. B 86:035120. doi: 10.1103/PhysRevB.86.035120

CrossRef Full Text | Google Scholar

Cederbaum, L. S., and Domcke, W. (1974). On the vibrational structure in photoelectron spectra by the method of Green's functions. J. Chem. Phys (New York, NY), 60, 2878–2889. doi: 10.1063/1.1681457

CrossRef Full Text | Google Scholar

Cederbaum, L. S., and Domcke, W. (2007). “Theoretical aspects of ionization potentials and photoelectron spectroscopy: a green's function approach,” in Advances in Chemical Physics, eds I. Prigogine and S. A. Rice (New York, NY: John Wiley & Sons, Ltd), 205–344.

Chantis, A. N., van Schilfgaarde, M., and Kotani, T. (2007). Quasiparticle self-consistent GW method applied to localized 4f electron systems. Phys. Rev. B 76:165126. doi: 10.1103/PhysRevB.76.165126

CrossRef Full Text | Google Scholar

Cheiwchanchamnangij, T., and Lambrecht, W. R. L. (2011). Band structure parameters of wurtzite and zinc-blende GaAs under strain in the GW approximation. Phys. Rev. B 84:035203. doi: 10.1103/PhysRevB.84.035203

CrossRef Full Text | Google Scholar

Cheiwchanchamnangij, T., and Lambrecht, W. R. L. (2012). Quasiparticle band structure calculation of monolayer, bilayer, and bulk MoS2. Phys. Rev. B 85:205302. doi: 10.1103/PhysRevB.85.205302

CrossRef Full Text | Google Scholar

Chen, W., and Pasquarello, A. (2015a). Accurate band gaps of extended systems via efficient vertex corrections in GW. Phys. Rev. B 92:041115. doi: 10.1103/PhysRevB.92.041115

CrossRef Full Text | Google Scholar

Chen, W., and Pasquarello, A. (2015b). First-principles determination of defect energy levels through hybrid density functionals and GW. J. Phys. Condens. Matter 27:133202. doi: 10.1088/0953-8984/27/13/133202

PubMed Abstract | CrossRef Full Text | Google Scholar

Choi, S., Kutepov, A., Haule, K., van Schilfgaarde, M., and Kotliar, G. (2016). First-principles treatment of Mott insulators: linearized QSGW+DMFT approach. NPJ Quantum Mat. 1:16001. doi: 10.1038/npjquantmats.2016.1

CrossRef Full Text | Google Scholar

Cocchi, C., Breuer, T., Witte, G., and Draxl, C. (2018). Polarized absorbance and Davydov splitting in bulk and thin-film pentacene polymorphs. Phys. Chem. Chem. Phys. 20, 29724–29736. doi: 10.1039/C8CP06384B

PubMed Abstract | CrossRef Full Text | Google Scholar

Dahlen, N. E., Stan, A., and Leeuwen, R. (2006a). Nonequilibrium Green function theory for excitation and transport in atoms and molecules. J. Phys. Conf. Ser. 35, 324–339. doi: 10.1088/1742-6596/35/1/030

CrossRef Full Text | Google Scholar

Dahlen, N. E., van Leeuwen, R., and von Barth, U. (2006b). Variational energy functionals of the Green function and of the density tested on molecules. Phys. Rev. A 73:012511. doi: 10.1103/PhysRevA.73.012511

CrossRef Full Text | Google Scholar

Dauth, M., Caruso, F., Kümmel, S., and Rinke, P. (2016). Piecewise linearity in the GW approximation for accurate quasiparticle energy predictions. Phys. Rev. B 93:121115. doi: 10.1103/PhysRevB.93.121115

CrossRef Full Text | Google Scholar

Debbichi, L., Eriksson, O., and Lebègue, S. (2014). Electronic structure of two-dimensional transition metal dichalcogenide bilayers from ab initio theory. Phys. Rev. B 89:205311. doi: 10.1103/PhysRevB.89.205311

CrossRef Full Text | Google Scholar

Deisz, J. J., Eguiluz, A. G., and Hanke, W. (1993). Quasiparticle theory versus density-functional theory at a metal surface. Phys. Rev. Lett. 71, 2793–2796. doi: 10.1103/PhysRevLett.71.2793

PubMed Abstract | CrossRef Full Text | Google Scholar

Del Ben, M., da Jornada, F. H., Canning, A., Wichmann, N., Raman, K., Sasanka, R., et al. (2019). Large-scale GW calculations on pre-exascale HPC systems. Comput. Phys. Commun. 235, 187–195. doi: 10.1016/j.cpc.2018.09.003

CrossRef Full Text | Google Scholar

Del Sole, R., Adragna, G., Olevano, V., and Reining, L. (2003). Long-range behavior and frequency dependence of exchange-correlation kernels in solids. Phys. Rev. B 67:045207. doi: 10.1103/PhysRevB.67.045207

CrossRef Full Text | Google Scholar

Del Sole, R., Reining, L., and Godby, R. W. (1994). GWΓ approximation for electron self-energies in semiconductors and insulators. Phys. Rev. B 49:8024. doi: 10.1103/PhysRevB.49.8024

PubMed Abstract | CrossRef Full Text | Google Scholar

Delaney, K., García-González, P., Rubio, A., Rinke, P., and Godby, R. W. (2004). Comment on “band-gap problem in semiconductors revisited: effects of core states and many-body self-consistency”. Phys. Rev. Lett. 93:249701. doi: 10.1103/PhysRevLett.93.249701

PubMed Abstract | CrossRef Full Text | Google Scholar

Deslippe, J., Samsonidze, G., Strubbe, D. A., Jain, M., Cohen, M. L., and Louie, S. G. (2012). BerkeleyGW: a massively parallel computer package for the calculation of the quasiparticle and optical properties of materials and nanostructures. Comput. Phys. Commun. 183, 1269–1289. doi: 10.1016/j.cpc.2011.12.006

CrossRef Full Text | Google Scholar

Devaux, N., Casula, M., Decremps, F., and Sorella, S. (2015). Electronic origin of the volume collapse in cerium. Phys. Rev. B 91:081101. doi: 10.1103/PhysRevB.91.081101

CrossRef Full Text | Google Scholar

Dori, N., Menon, M., Kilian, L., Sokolowski, M., Kronik, L., and Umbach, E. (2006). Valence electronic structure of gas-phase 3,4,9,10-perylene tetracarboxylic acid dianhydride: experiment and theory. Phys. Rev. B 73:195208. doi: 10.1103/PhysRevB.73.195208

CrossRef Full Text | Google Scholar

Dose, V. (1985). Momentum-resolved inverse photoemission. Surf. Sci. Rep. 5:337. doi: 10.1016/0167-5729(85)90006-8

CrossRef Full Text | Google Scholar

Drissi, L., and Ramadan, F. (2015a). Excitonic effects in GeC hybrid: many-body Green's function calculations. Physica E 74, 377–381. doi: 10.1016/j.physe.2015.07.030

CrossRef Full Text | Google Scholar

Drissi, L., and Ramadan, F. (2015b). Many body effects study of electronic and optical properties of silicene graphene hybrid. Physica E 68, 38–41. doi: 10.1016/j.physe.2014.12.009

CrossRef Full Text | Google Scholar

Drüppel, M., Deilmann, T., Noky, J., Marauhn, P., Krüger, P., and Rohlfing, M. (2018). Electronic excitations in transition metal dichalcogenide monolayers from an LDA+GdW approach. Phys. Rev. B 98:155433. doi: 10.1103/PhysRevB.98.155433

CrossRef Full Text | Google Scholar

Duchemin, I., Jacquemin, D., and Blase, X. (2016). Combining the GW formalism with the polarizable continuum model: A state-specific non-equilibrium approach. J. Chem. Phys. 144:164106. doi: 10.1063/1.4946778

PubMed Abstract | CrossRef Full Text | Google Scholar

Dunning, T. H. (1989). Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. J. Chem. Phys. 90, 1007–1023. doi: 10.1063/1.456153

CrossRef Full Text | Google Scholar

Dvorak, M., Chen, X.-J., and Wu, Z. (2014). Quasiparticle energies and excitonic effects in dense solid hydrogen near metallization. Phys. Rev. B 90:035103. doi: 10.1103/PhysRevB.90.035103

CrossRef Full Text | Google Scholar

Dvorak, M., Golze, D., and Rinke, P. (2018). A quantum embedding theory in the screened Coulomb interaction: Combining configuration interaction with GW/BSE. arXiv:1810.12005.

Google Scholar

Dvorak, M., and Rinke, P. (2019). Dynamical configuration interaction: quantum embedding that combines wave functions and Green's functions. Phys. Rev. B 99:115134. doi: 10.1103/PhysRevB.99.115134

CrossRef Full Text | Google Scholar

Dvorak, M., and Wu, Z. (2015). Tunable many-body interactions in semiconducting graphene: giant excitonic effect and strong optical absorption. Phys. Rev. B 92:035422. doi: 10.1103/PhysRevB.92.035422

CrossRef Full Text | Google Scholar

Dyson, F. J. (1949a). The radiation theories of tomonaga, tchwinger, and feynman. Phys. Rev. 75:486. doi: 10.1103/PhysRev.75.486

CrossRef Full Text | Google Scholar

Dyson, F. J. (1949b). The S matrix in quantum electrodynamics. Phys. Rev. 75:1736. doi: 10.1103/PhysRev.75.1736

CrossRef Full Text | Google Scholar

Egelhoff, W. F. (1987). Core-level binding-energy shifts at surfaces and in solids. Surf. Sci. Rep. 6, 253–415. doi: 10.1016/0167-5729(87)90007-0

CrossRef Full Text | Google Scholar

Egger, D. A., Weissman, S., Refaely-Abramson, S., Sharifzadeh, S., Dauth, M., Baer, R., et al. (2014). Outer-valence electron spectra of prototypical aromatic heterocycles from an optimally tuned range-separated hybrid functional. J. Chem. Theor. Comput. 10, 1934–1952. doi: 10.1021/ct400956h

PubMed Abstract | CrossRef Full Text | Google Scholar

Engel, G. E., and Farid, B. (1993). Generalized plasmon-pole model and plasmon band structures of crystals. Phys. Rev. B 47, 15931–15934. doi: 10.1103/PhysRevB.47.15931

PubMed Abstract | CrossRef Full Text | Google Scholar

Enkovaara, J., Rostgaard, C., Mortensen, J. J., Chen, J., Dułak, M., Ferrighi, L., et al. (2010). Electronic structure calculations with GPAW: a real-space implementation of the projector augmented-wave method. J. Phys. Condens. Matter 22:253202. doi: 10.1088/0953-8984/22/25/253202

PubMed Abstract | CrossRef Full Text | Google Scholar

Erhart, P., Schleife, A., Sadigh, B., and Åberg, D. (2014). Quasiparticle spectra, absorption spectra, and excitonic properties of NaI and SrI2 from many-body perturbation theory. Phys. Rev. B 89:075132. doi: 10.1103/PhysRevB.89.075132

CrossRef Full Text | Google Scholar

Ernzerhof, M., and Scuseria, G. E. (1999). Assessment of the Perdew-Burke-Ernzerhof exchange-correlation functional. J. Chem. Phys. 110:5029. doi: 10.1063/1.478401

CrossRef Full Text | Google Scholar

Eshuis, H., Bates, J. E., and Furche, F. (2012). Electron correlation methods based on the random phase approximation. Theor. Chem. Acc. 131:1084. doi: 10.1007/s00214-011-1084-8

CrossRef Full Text | Google Scholar

Espejo, C., Rangel, T., Romero, A. H. X., and Rignanese, G.-M. (2013). Band structure tunability in MoS2 under interlayer compression: A DFT and GW study. Phys. Rev. B 87:245114. doi: 10.1103/PhysRevB.87.245114

CrossRef Full Text | Google Scholar

Faber, C., Attaccalite, C., Olevano, V., Runge, E., and Blase, X. (2011). First-principles GW calculations for DNA and RNA nucleobases. Phys. Rev. B 83:115123. doi: 10.1103/PhysRevB.83.115123

CrossRef Full Text | Google Scholar

Faber, C., Boulanger, P., Attaccalite, C., Duchemin, I., and Blase, X. (2014). Excited states properties of organic molecules: from density functional theory to the GW and Bethe-Salpeter Green's function formalisms. Philos. Trans. R. Soc. A 372:20130271. doi: 10.1098/rsta.2013.0271

PubMed Abstract | CrossRef Full Text | Google Scholar

Faber, C., Duchemin, I., Deutsch, T., and Blase, X. (2012). Many-body Green's function study of coumarins for dye-sensitized solar cells. Phys. Rev. B 86:155315. doi: 10.1103/PhysRevB.86.155315

CrossRef Full Text | Google Scholar

Fakhrabad, D. V., Shahtahmasebi, N., and Ashhadi, M. (2015). Optical excitations and quasiparticle energies in the AlN monolayer honeycomb structure. Superlatt. Microstruct. 79, 38–44. doi: 10.1016/j.spmi.2014.12.012

CrossRef Full Text | Google Scholar

Fakhrabad, D. V., Shahtamasebi, N., and Ashhadi, M. (2014). Quasiparticle energies and optical excitations in the GaAs monolayer. Physica E 59, 107–109. doi: 10.1016/j.physe.2013.12.019

CrossRef Full Text | Google Scholar

Faleev, S. V., van Schilfgaarde, M., and Kotani, T. (2004). All-electron self-consistent GW approximation: application to Si, MnO, and NiO. Phys. Rev. Lett. 93:126406. doi: 10.1103/PhysRevLett.93.126406

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferreira, F., and Ribeiro, R. M. (2017). Improvements in the GW and Bethe-Salpeter-equation calculations on phosphorene. Phys. Rev. B 96:115431. doi: 10.1103/PhysRevB.96.115431

CrossRef Full Text | Google Scholar

Ferri, N., DiStasio, R. A., Ambrosetti, A., Car, R., and Tkatchenko, A. (2015). Electronic properties of molecules and surfaces with a self-consistent interatomic van der Waals density functional. Phys. Rev. Lett. 114:176802. doi: 10.1103/PhysRevLett.114.176802

PubMed Abstract | CrossRef Full Text | Google Scholar

Fetter, A. L., and Walecka, J. D. (1971). Quantum Theory of Many-Particle Systems. New York, NY: Courier Dover Publications.

Google Scholar

Filip, M. R., Verdi, C., and Giustino, F. (2015). GW band structures and carrier effective masses of CH3NH3PbI3 and hypothetical perovskites of the type APbI3: A = NH4, PH4, AsH4, and SbH4. J. Phys. Chem. C 119, 25209–25219. doi: 10.1021/acs.jpcc.5b07891

CrossRef Full Text | Google Scholar

Fleszar, A., and Hanke, W. (2005). Electronic structure of IIB-VI semiconductors in the GW approximation. Phys. Rev. B 71:045207. doi: 10.1103/PhysRevB.71.045207

CrossRef Full Text | Google Scholar

Foerster, D., Koval, P., and Sánchez-Portal, D. (2011). An O(N3) implementation of Hedin's GW approximation for molecules. J. Chem. Phys. 135:074105. doi: 10.1063/1.3624731

PubMed Abstract | CrossRef Full Text | Google Scholar

Fonari, A., Sutton, C., Brédas, J.-L., and Coropceanu, V. (2014). Impact of exact exchange in the description of the electronic structure of organic charge-transfer molecular crystals. Phys. Rev. B 90:165205. doi: 10.1103/PhysRevB.90.165205

CrossRef Full Text | Google Scholar

Fratesi, G., Brivio, G. P., Rinke, P., and Godby, R. (2003). Image resonance in the many-body density of states at a metal surface. Phys. Rev. B 68:195404. doi: 10.1103/PhysRevB.68.195404

CrossRef Full Text | Google Scholar

Freysoldt, C., Eggert, P., Rinke, P., Schindlmayr, A., Godby, R. W., and Scheffler, M. (2007). Dielectric anisotropy in the GW space-time method. Comput. Phys. Commun. 176:1. doi: 10.1016/j.cpc.2006.07.018

CrossRef Full Text | Google Scholar

Freysoldt, C., Eggert, P., Rinke, P., Schindlmayr, A., and Scheffler, M. (2008). Screening in 2D: GW calculations for surfaces and thin films using the repeated-slab approach. Phys. Rev. B 77:235428. doi: 10.1103/PhysRevB.77.235428

CrossRef Full Text | Google Scholar

Freysoldt, C., Rinke, P., and Scheffler, M. (2009). Controlling polarization at insulating surfaces: quasiparticle calculations for molecules adsorbed on insulator films. Phys. Rev. Lett. 103:056803. doi: 10.1103/PhysRevLett.103.056803

PubMed Abstract | CrossRef Full Text | Google Scholar

Friedrich, C., Betzinger, M., Schlipf, M., Blügel, S., and Schindlmayr, A. (2012). Hybrid functionals and GW approximation in the FLAPW method. J. Phys. Condens. Matter 24:293201. doi: 10.1088/0953-8984/24/29/293201

PubMed Abstract | CrossRef Full Text | Google Scholar

Friedrich, C., Blügel, S., and Schindlmayr, A. (2010). Efficient implementation of the GW approximation within the all-electron FLAPW method. Phys. Rev. B 81:125102. doi: 10.1103/PhysRevB.81.125102

CrossRef Full Text | Google Scholar

Friedrich, C., Müller, M. C., and Blügel, S. (2011). Band convergence and linearization error correction of all-electron GW calculations: the extreme case of zinc oxide. Phys. Rev. B 83:081101. doi: 10.1103/PhysRevB.83.081101

CrossRef Full Text | Google Scholar

Friedrich, C., Schindlmayr, A., Blügel, S., and Kotani, T. (2006). Elimination of the linearization error in GW calculations based on the linearized augmented-plane-wave method. Phys. Rev. B 74:045104. doi: 10.1103/PhysRevB.74.045104

CrossRef Full Text | Google Scholar

Fuchs, F., Furthmüller, J., Bechstedt, F., Shishkin, M., and Kresse, G. (2007). Quasiparticle band structure based on a generalized Kohn-Sham scheme. Phys. Rev. B 76:115109. doi: 10.1103/PhysRevB.76.115109

CrossRef Full Text | Google Scholar

Fuchs, M., and Scheffler, M. (1999). Ab initio pseudopotentials for electronic structure calculations of poly-atomic systems using density-functional theory. Comput. Phys. Commun. 119:67. doi: 10.1016/S0010-4655(98)00201-X

CrossRef Full Text | Google Scholar

Fuggle, J. C., and Inglesfield, J. E. (eds.). (1992). Unoccupied Electronic States. Berlin: Springer-Verlag.

Furthmüller, J., Hahn, P. H., Fuchs, F., and Bechstedt, F. (2005). Band structures and optical spectra of InN polymorphs: Influence of quasiparticle and excitonic effects. Phys. Rev. B 72:205106. doi: 10.1103/PhysRevB.72.205106

CrossRef Full Text | Google Scholar

Galitskii, V., and Migdal, A. (1958). Application of quantum field theory methods to the many body problem. Sov. Phys. JETP 7:96.

Google Scholar

Gallandi, L., and Körzdörfer, T. (2015). Long-range corrected DFT meets GW: vibrationally resolved photoelectron spectra from first principles. J. Chem. Theor. Comput. 11, 5391–5400. doi: 10.1021/acs.jctc.5b00820

PubMed Abstract | CrossRef Full Text | Google Scholar

Gallandi, L., Marom, N., Rinke, P., and Körzdörfer, T. (2016). Accurate ionization potentials and electron affinities of acceptor molecules II: non-empirically tuned long-range corrected hybrid functionals. J. Chem. Theor. Comput. 12, 605–614. doi: 10.1021/acs.jctc.5b00873

PubMed Abstract | CrossRef Full Text | Google Scholar

Ganesan, V. D. S., Linghu, J., Zhang, C., Feng, Y. P., and Shen, L. (2016). Heterostructures of phosphorene and transition metal dichalcogenides for excitonic solar cells: a first-principles study. Appl. Phys. Lett. 108:122105. doi: 10.1063/1.4944642

CrossRef Full Text | Google Scholar

Gao, W., Xia, W., Gao, X., and Zhang, P. (2016). Speeding up GW calculations to meet the challenge of large scale quasiparticle predictions. Sci. Rep. 6:36849. doi: 10.1038/srep36849

PubMed Abstract | CrossRef Full Text | Google Scholar

García, A., Verstraete, M. J., Pouillon, Y., and Junquera, J. (2018). The PSML format and library for norm-conserving pseudopotential data curation and interoperability. Comput. Phys. Commun. 227, 51–71. doi: 10.1016/j.cpc.2018.02.011

CrossRef Full Text | Google Scholar

García-González, P., and Godby, R. W. (2002). Many-body GW calculations of ground-state properties: quasi-2D electron systems and van der Waals forces. Phys. Rev. Lett. 88:056406. doi: 10.1103/PhysRevLett.88.056406

PubMed Abstract | CrossRef Full Text | Google Scholar

García-González, P. P., and Godby, R. W. (2001). Self-consistent calculation of total energies of the electron gas using many-body perturbation theory. Phys. Rev. B 63:075112. doi: 10.1103/PhysRevB.63.075112

CrossRef Full Text | Google Scholar

García-Lastra, J. M., Rostgaard, C., Rubio, A., and Thygesen, K. S. (2009). Polarization-induced renormalization of molecular levels at metallic and semiconducting surfaces. Phys. Rev. B 80:245427. doi: 10.1103/PhysRevB.80.245427

CrossRef Full Text | Google Scholar

Gatti, M., Bruneval, F., Olevano, V., and Reining, L. (2007). Understanding correlations in vanadium dioxide from first principles. Phys. Rev. Lett. 99:266402. doi: 10.1103/PhysRevLett.99.266402

PubMed Abstract | CrossRef Full Text | Google Scholar

Gatti, M., and Guzzo, M. (2013). Dynamical screening in correlated metals: spectral properties of SrVO3 in the GW approximation and beyond. Phys. Rev. B 87:155147. doi: 10.1103/PhysRevB.87.155147

CrossRef Full Text | Google Scholar

Gatti, M., Panaccione, G., and Reining, L. (2015). Effects of low-energy excitations on spectral properties at higher binding energy: the metal-insulator transition of VO2. Phys. Rev. Lett. 114:116402. doi: 10.1103/PhysRevLett.114.116402

CrossRef Full Text | Google Scholar

Geim, A. K., and Novoselov, K. S. (2007). The rise of graphene. Nat. Mater. 6, 183–191. doi: 10.1038/nmat1849

PubMed Abstract | CrossRef Full Text

Georges, A., and Kotliar, G. (1992). Hubbard model in infinite dimensions. Phys. Rev. B 45, 6479–6483. doi: 10.1103/PhysRevB.45.6479

PubMed Abstract | CrossRef Full Text | Google Scholar

Gerosa, M., Bottani, C. E., Valentin, C. D., Onida, G., and Pacchioni, G. (2018a). Accuracy of dielectric-dependent hybrid functionals in the prediction of optoelectronic properties of metal oxide semiconductors: a comprehensive comparison with many-body GW and experiments. J. Phys. Condens. Matter 30:044003. doi: 10.1088/1361-648X/aa9725

PubMed Abstract | CrossRef Full Text | Google Scholar

Gerosa, M., Gygi, F., Govoni, M., and Galli, G. (2018b). The role of defects and excess surface charges at finite temperature for optimizing oxide photoabsorbers. Nat. Mater. 17, 1122–1127. doi: 10.1038/s41563-018-0192-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Giantomassi, M., Stankovski, M., Shaltaf, R., Grüning, M., Bruneval, F., Rinke, P., et al. (2011). Electronic Properties of Interfaces and Defects from Many-Body Perturbation Theory: Recent Developments and Applications. Weinheim: John Wiley & Sons, Ltd., 33–60.

Google Scholar

Giorgi, G., Palummo, M., Chiodo, L., and Yamashita, K. (2011). Excitons at the (001) surface of anatase: Spatial behavior and optical signatures. Phys. Rev. B 84:073404. doi: 10.1103/PhysRevB.84.073404

CrossRef Full Text | Google Scholar

Giovannantonio, M. D., Urgel, J. I., Beser, U., Yakutovich, A. V., Wilhelm, J., Pignedoli, C. A., et al. (2018). On-surface synthesis of indenofluorene polymers by oxidative five-membered ring formation. J. Am. Chem. Soc. 140, 3532–3536. doi: 10.1021/jacs.8b00587

PubMed Abstract | CrossRef Full Text | Google Scholar

Giustino, F., Cohen, M. L., and Louie, S. G. (2010a). GW method with the self-consistent Sternheimer equation. Phys. Rev. B 81:115105. doi: 10.1103/PhysRevB.81.115105

CrossRef Full Text | Google Scholar

Giustino, F., Louie, S. G., and Cohen, M. L. (2010b). Electron-phonon renormalization of the direct band gap of diamond. Phys. Rev. Lett. 105:265501. doi: 10.1103/PhysRevLett.105.265501

PubMed Abstract | CrossRef Full Text | Google Scholar

Godby, R. W., and Needs, R. J. (1989). Metal-insulator transition in Kohn-Sham theory and quasiparticle theory. Phys. Rev. Lett. 62, 1169–1172. doi: 10.1103/PhysRevLett.62.1169

PubMed Abstract | CrossRef Full Text | Google Scholar

Godby, R. W., Schlüter, M., and Sham, L. J. (1986). Accurate exchange-correlation potential for silicon and its discontinuity on addition of an electron. Phys. Rev. Lett. 56, 2415–2418. doi: 10.1103/PhysRevLett.56.2415

PubMed Abstract | CrossRef Full Text | Google Scholar

Godby, R. W., Schlüter, M., and Sham, L. J. (1987a). Quasiparticle energies in GaAs and AlAs. Phys. Rev. B 35, 4170–4171. doi: 10.1103/PhysRevB.35.4170

PubMed Abstract | CrossRef Full Text | Google Scholar

Godby, R. W., Schlüter, M., and Sham, L. J. (1987b). Trends in self-energy operators and their corresponding exchange-correlation potentials. Phys. Rev. B 36, 6497–6500. doi: 10.1103/PhysRevB.36.6497

PubMed Abstract | CrossRef Full Text | Google Scholar

Godby, R. W., Schlüter, M., and Sham, L. J. (1988). Self-energy operators and exchange-correlation potentials in semiconductors. Phys. Rev. B 37:10159. doi: 10.1103/PhysRevB.37.10159

PubMed Abstract | CrossRef Full Text | Google Scholar

Goedecker, S., Teter, M., and Hutter, J. (1996). Separable dual-space Gaussian pseudopotentials. Phys. Rev. B 54, 1703–1710. doi: 10.1103/PhysRevB.54.1703

PubMed Abstract | CrossRef Full Text | Google Scholar

Golze, D., Wilhelm, J., van Setten, M. J., and Rinke, P. (2018). Core-level binding energies from GW: an efficient full-frequency approach within a localized basis. J. Chem. Theor. Comput. 14, 4856–4869. doi: 10.1021/acs.jctc.8b00458

PubMed Abstract | CrossRef Full Text | Google Scholar

Gonze, X. (1995). Adiabatic density-functional perturbation theory. Phys. Rev. A 52, 1096–1114. doi: 10.1103/PhysRevA.52.1096

PubMed Abstract | CrossRef Full Text | Google Scholar

Gonze, X. (1997). First-principles responses of solids to atomic di