Electronic band structure from first-principles Green’s function approach: theory and implementations

Abstract

Electronic band structure is one of the most important intrinsic properties of a material, and is in particular crucial in electronic, photo-electronic and photo- catalytic applications. Kohn-Sham Density-functional theory (KS-DFT) within currently available local or semi-local approximations to the exchange-correlation energy functional is problematic for the description of electronic band structure. Many-body perturbation theory based on Green’s function (GF) provides a rigorous framework to describe excited-state properties of materials. The central ingredient of the GF-based many-body perturbation theory is the exchange- correlation self-energy, which accounts for all non-classical electron-electron interaction effects beyond the Hartree theory, and formally can be obtained by solving a set of complicated integro-differential equations, named Hedin’s equations. The GW approximation, in which the self-energy is simply a product of Green’s function and the screened Coulomb interaction (W), is currently the most accurate first-principles approach to describe electronic band structure properties of extended systems. Compared to KS-DFT, the computational efforts required for GW calculations are much larger. Various numerical techniques or approximations have been developed to apply GW for realistic systems. In this paper, we give an overview of the theory of first-principles Green’s function approach in the GW approximation and review the state of the art for the implementation of GW in different representations and with different treatment of the frequency dependence. It is hoped that further methodological developments will be inspired by this work so that the approach can be applied to more complicated and scientifically more interesting systems.

Keywords

electronic band structure / many-body perturbation theory / GW approximation

Introduction

Electronic band structure is one of the most important intrinsic properties of a material that has far-reaching effects on electronic, optical, photo-electronic and photo-catalytic applications. Experimentally electronic properties of a material of interest are often measured in terms of its response to some external perturbation, using either electron or photon as the probe. Among the most widely used spectroscopic techniques are photo-emission, inverse photo-emission and optical absorption spectroscopies [1-3]. In the photo-emission spectroscopy (PES) [1], photons of given energy impinge on the sample and liberate electrons from their occupied levels, either deep core states or valence states; by monitoring the kinetic energy of photoelectrons, the information of electronic occupied states can be probed. The inverse photo-emission spectroscopy (IPS) is essentially a time-reversed process of PES: free electrons of given energy are injected into the sample and fill originally unoccupied states; the redundant energy is released as radiation, whose intensity as a function of energy reveals the information of unoccupied (or conduction band) states. In essence, PES and IPS measure single-particle excitations of the system, with the total electron number N changing to N - 1 or N + 1, which are often termed as quasiparticle (QP) excitations. The optical absorption spectroscopy (OAS), on the other hand, measures neutral excitations, in which electrons are excited from the ground state to excited states, corresponding to creating a hole and an electron simultaneously in a QP picture. If we neglect the interaction between photo-excited electrons and holes, OAS can be regarded as the combination of PES and IPS in a single process. In this paper, we are mainly concerned with the theoretical description of QP excitations from a first-principles perspective.

Kohn-Sham (KS) density functional theory (DFT) [4,5] within the local density or generalized gradient approximation (LDA/GGA) to the exchange-correlation (xc) energy functional has become “the standard model” for first-principles electronic structure calculations of extended systems [6]. By mapping the many-electron interacting system to a fictitious non-interacting (Kohn-Sham) system, which has the same ground state electron density as the interacting one, the original highly complicated many-body problem is transformed to solving a single-particle equation that formally resembles a mean-field theory, but is exact in principles. Within LDA/GGA, KS-DFT can provide accurate descriptions of energetic and structural properties for many materials with feasible computational efforts. As a byproduct, KS-DFT also gives a set of single-particle energies and wave functions, which are, rigorously speaking, auxiliary quantities introduced only to calculate electron density, and therefore do not have any physical meanings. In practice, however, they are often used to interpret electronic quasi-particle band structure of materials as probed by PES/IPS [1]. This practice, however, has to be exercised with caution. Even for many simple sp semiconductors, KS-DFT within LDA/GGA gives band gaps that are systematically underestimated when compared to experiment, and can predict metallic ground state for small-gap semiconductors, e.g., Ge, InN and so on [7].

A completely different theoretical framework for electronic band structure is provided by many-body perturbation theory (MBPT), formulated in terms of one-body Green’s function [8,9], in the GW approximation [10,11]. In this approach, the exchange-correlation self-energy Σ, the central quantity in GF-based MBPT, is obtained as a simple product of single-particle Green’s function (G) and the screened Coulomb interaction (W), hence termed GW. In practice, quasi-particle energies in the GW approximation are often calculated as a first-order correction to eigenenergies of some reference single-particle Hamiltonian H₀, and both G and W are calculated using eigen-energies and eigen-functions of H₀, hence called G₀W₀. Since the middle of 1980s [12,13], the G₀W₀ approach based on the LDA or GGA H₀ has become the method of choice for quasi-particle band structures of various semiconductors [11]. From a practical point of view, the computational demand of GW calculations is much heavier than KS-DFT within LDA/GGA. As a result, various physical or numerical approximations have been developed to facilitate the applications of GW to physically interesting systems. Recent years have seen intensive developments in the methodology of GW with the goal of either going beyond LDA-based G₀W₀ [14-16], or significantly extending the capability of current implementations so that more complex systems can be reached [17-21]. The main goal of this paper is to present an overview of many-body perturbation theory for quasi-particle excitations based on the Green’s function, and review various numerical techniques developed for the implementation of the GW approach. The topic covered in this paper partly overlap with that in Ref. [22], but with different focus. Considering the recent increasing interest to apply the GW for large systems that are currently out of reach of current implementations, it is hoped that this work will inspire further methodological developments so that routine applications of GW and related first-principles methods to more complicated and scientifically more interesting systems become feasible.

The paper is organized as follows. In the next section we present the general theoretical framework in which the GW approximation is derived. In Section 3, main numerical techniques used in the implementation of the GW method are discussed including the basis used to represent the GW equations and the treatment of the frequency dependence by different approaches. Section 4 concludes the paper with a few general remarks.

Many-body perturbation theory and the GW approximation

In this section we present the fundamentals of first principles many-body perturbation theory based on one-body Green’s function. More complete formalisms can be found in Refs. [9,11,23,24].

Green’s function

The Hamiltonian of a many-electron interacting system, which is determined by the external potential V_ext(x) and the number of electrons (N), reads in the second-quantization representation [8]

(1)\widehat{H}={\displaystyle \int d\mathbf{x}{\widehat{\psi }}^{†}(\mathbf{x}t)h_{\mathrm{0}}(\mathbf{x})\widehat{\psi }(\mathbf{x}t)}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \int d\mathbf{x}d{\mathbf{x}}^{\prime }v(\mathbf{x},{\mathbf{x}}^{\prime }){\widehat{\psi }}^{†}(\mathbf{x}t){\widehat{\psi }}^{†}({\mathbf{x}}^{\prime }t)\widehat{\psi }({\mathbf{x}}^{\prime }t)\widehat{\psi }(\mathbf{x}t)}

where we use \mathbf{x}\equiv \{\mathbf{r},\sigma \} to denote both spatial and spin coordinates. \widehat{\psi }(\mathbf{x}t) and {\widehat{\psi }}^{†}(\mathbf{x}t) are the annihilation and creation field operators in the Heisenberg picture, respectively. h_{\mathrm{0}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\nabla }^{\mathrm{2}}+V_{ext}(\mathbf{x}) is the one-body part of the Hamiltonian, and v(\mathbf{x},{\mathbf{x}}^{\mathbf{\text{'}}})\equiv {\displaystyle \frac{\mathrm{1}}{\left|\mathbf{r}-{\mathbf{r}}^{\mathbf{\text{'}}}\right|}}{\delta }_{\sigma ,{\sigma }^{\prime }} is the bare Coulomb interaction, which is spin-independent. Atomic units are used through the paper.

One-body Green’s function G(xt, x′t′), also called the propagator, is then defined as [8]

(2)G(\mathbf{x}t;{\mathbf{x}}^{\prime }{t}^{\prime })=-i\langle N|\mathbf{T}\left[\widehat{\psi }(\mathbf{x}t){\widehat{\psi }}^{†}({\mathbf{x}}^{\prime }{t}^{\prime })\right]|N\rangle

where T is the time-ordering operator, and |N\rangle denotes the ground state of the N-electron interacting system. Based on the Heisenberg equation of motion (EOM) for the field operator \widehat{\psi }(\mathbf{x}t)

(3)i{\displaystyle \frac{\delta \widehat{\psi }(\mathbf{x}t)}{\delta t}}=\left[\widehat{\psi }(\mathbf{x}t),\widehat{H}\right]

we can obtain the equation of motion for the Green’s function

(4)\left[\mathrm{i}{\displaystyle \frac{\delta }{\delta t}}-h_{\mathrm{0}}(\mathbf{x})\right]G(\mathbf{x}t,{\mathbf{x}}^{\prime }{t}^{\prime })-i{\displaystyle \int d{\mathbf{x}}^{″}v(\mathbf{x},{\mathbf{x}}^{″})}\times \langle N|\mathbf{T}\left[{\widehat{\psi }}^{†}({\mathbf{x}}^{″}t)\widehat{\psi }({\mathbf{x}}^{″}t)\widehat{\psi }(\mathbf{x}t){\widehat{\psi }}^{†}({\mathbf{x}}^{\prime }t)\right]|N\rangle =\delta (\mathbf{x}-{\mathbf{x}}^{\prime })\delta (t-t^{\prime })

The last term in the left side of the equation above is a special case of the two-body Green’s function, generally defined as

(5)G_{\mathrm{2}}({\mathbf{x}}_{\mathrm{1}}{t}_{\mathrm{1}},{\mathbf{x}}_{\mathrm{2}}{t}_{\mathrm{2}},{\mathbf{x}}_{\mathrm{3}}{t}_{\mathrm{3}},{\mathbf{x}}_{\mathrm{4}}{t}_{\mathrm{4}})=i^{\mathrm{2}}\langle N|\mathbf{T}\left[\widehat{\psi }({\mathbf{x}}_{\mathrm{1}}{t}_{\mathrm{1}})\widehat{\psi }({\mathbf{x}}_{\mathrm{3}}{t}_{\mathrm{3}}){\widehat{\psi }}^{†}({\mathbf{x}}_{\mathrm{4}}{t}_{\mathrm{4}}){\widehat{\psi }}^{†}({\mathbf{x}}_{\mathrm{2}}{t}_{\mathrm{2}})\right]|N\rangle

Equation (4) can therefore be written as

(6)\left[i{\displaystyle \frac{\delta }{\delta t}}-h_{\mathrm{0}}(\mathbf{x})\right]G(\mathbf{x}t,{\mathbf{x}}^{\prime }{t}^{\prime })+i{\displaystyle \int d{\mathbf{x}}^{″}v(\mathbf{x},{\mathbf{x}}^{″})}\times G_{\mathrm{2}}(\mathbf{x}t,{\mathbf{x}}^{\prime }{t}^{\prime },{\mathbf{x}}^{″}t,{\mathbf{x}}^{″}{t}^{+})=\delta (\mathbf{x}-{\mathbf{x}}^{\prime })\delta (t-t^{\prime })

with t^{+}\equiv t+\eta, where η is an infinitesimal positive number. The equation of motion for G₂, which depends on the three-body Green’s function, can be derived in a similar way. This process can continue, leading to a set of hierarchical equations that are formally exact but useless in practice. To obtain practically feasible approximations, it is necessary to cut off the series or reformulate the equations in a form that is more accessible by approximations.

Equation (6) can be used as the starting point for simple approximations [23] such as the Hartree approximation

(7)G_{\mathrm{2}}(\mathbf{x}t,{\mathbf{x}}^{\prime }{t}^{\prime },{\mathbf{x}}^{″}t,{\mathbf{x}}^{″}{t}^{+})\simeq G(\mathbf{x}t,{\mathbf{x}}^{\prime }{t}^{\prime })G({\mathbf{x}}^{″}t,{\mathbf{x}}^{″}{t}^{+})

or the Hartree-Fock (HF) approximation

(8)G_{\mathrm{2}}(\mathbf{x}t,{\mathbf{x}}^{\prime }{t}^{\prime },{\mathbf{x}}^{″}t,{\mathbf{x}}^{″}{t}^{+})\simeq G(\mathbf{x}t,{\mathbf{x}}^{\prime }{t}^{\prime })G({\mathbf{x}}^{″}t,{\mathbf{x}}^{″}{t}^{+})+G(\mathbf{x}t,{\mathbf{x}}^{″}{t}^{+})G({\mathbf{x}}^{″}t,{\mathbf{x}}^{\prime }{t}^{\prime })

It is, however, difficult to incorporate high order approximations directly based on Eq. (6). Formally it is much more convenient to reformulate Eq. (6) in terms of the exchangecorrelation self-energy.

The self-energy and quasi-particle equation

The exchange-correlation self-energy Σ (or simply called self-energy) can be formally defined as

(9)i{\displaystyle \int d{\mathbf{x}}^{″}v(\mathbf{x},{\mathbf{x}}^{″})}{G}_{\mathrm{2}}(\mathbf{x}t,{\mathbf{x}}^{\prime }{t}^{\prime },{\mathbf{x}}^{″}t,{\mathbf{x}}^{″}{t}^{+})\equiv -V_H(\mathbf{x})G(\mathbf{x}t,{\mathbf{x}}^{\prime }{t}^{\prime })-{\displaystyle \int d{\mathbf{x}}^{″}d{t}^{″}{\displaystyle \sum (\mathbf{x}t,{\mathbf{x}}^{″}{t}^{″})G({\mathbf{x}}^{″}{t}^{″},{\mathbf{x}}^{\prime }{t}^{\prime })}}

where the first term on the right side contains the contribution of classical Coulomb repulsion (the Hartree potential), and all complexities related to the two-body and higher order Green’s functions are now wrapped up in the self-energy operator Σ. The equation of motion for the Green’s function now reads

(10)\left[i{\displaystyle \frac{\delta }{\delta t}}-h_{\mathrm{0}}(\mathbf{x})-V_H(\mathbf{x})\right]G(\mathbf{x}t,{\mathbf{x}}^{\prime }{t}^{\prime })-{\displaystyle \int d{\mathbf{x}}^{″}d{t}^{″}{\displaystyle \sum (\mathbf{x}t,{\mathbf{x}}^{″}{t}^{″})G({\mathbf{x}}^{″}{t}^{″},{\mathbf{x}}^{\prime }{t}^{\prime })=\delta (\mathbf{x}-{\mathbf{x}}^{\prime })\delta (t-t^{\prime })}}

For a time-independent system, both G(xt, x′t′) and Σ(xt, x′t′) depend on only the time difference τ = t-t′. Theoretically it is more convenient to work in the frequency domain, obtained by taking a Fourier transform with respect to τ, which leads to

(11)[\omega -h_{\mathrm{0}}(\mathbf{x})-V_H(\mathbf{x})]G(\mathbf{x},{\mathbf{x}}^{\prime };\omega )-{\displaystyle \int d{\mathbf{x}}^{″}{\displaystyle \sum (\mathbf{x},{\mathbf{x}}^{″};\omega )G({\mathbf{x}}^{″},{\mathbf{x}}^{\prime };\omega )=\delta (\mathbf{x}-{\mathbf{x}}^{\prime })}}

By introducing the non-interacting Green’s function G₀ as the solution of

(12)[\omega -h_{\mathrm{0}}(\mathbf{x})+V_H(\mathbf{x})]G_{\mathrm{0}}(\mathbf{x},{\mathbf{x}}^{\prime };\omega )=\delta (\mathbf{x}-{\mathbf{x}}^{\prime })

Equation (11) can be re-written as

(13)G(\mathbf{x},{\mathbf{x}}^{\prime };\omega )=G_{\mathrm{0}}(\mathbf{x},{\mathbf{x}}^{\prime };\omega )+{\displaystyle \int d{\mathbf{x}}_{\mathrm{1}}d{\mathbf{x}}_{\mathrm{2}}{G}_{\mathrm{0}}(\mathbf{x},{\mathbf{x}}_{\mathrm{1}};\omega )}\times {\displaystyle \sum ({\mathbf{x}}_{\mathrm{1}},{\mathbf{x}}_{\mathrm{2}};\omega )G({\mathbf{x}}_{\mathrm{2}},{\mathbf{x}}^{\prime };\omega )},

which is the Dyson’s equation for the Green’s function G.

The (one-body) Green’s function contains a lot of important information related to the ground and excited states of the system [8]. The poles of G(x, x′; ω) in the complex frequency domain give single-particle excitation properties (N → N±1) as probed by PES/IPS. In addition, the ground state total energy can also be obtained from G(x, x′, ω) [8]. For the description of the single-particle excitation, it is more convenient to work directly with the so-called quasi-particle equation, whose solutions give QP energies and wave functions. Here we present a schematic derivation of the QP equation [11]. According to the Green’s function theory of differential equations [25], G(x, x′, ω) as the solution of Eq. (11) can be formally written as

(14)G(\mathbf{x},{\mathbf{x}}^{\prime };\omega )={\displaystyle \sum _n{\displaystyle \frac{{\mathrm{\Psi }}_n(\mathbf{x};\omega ){\mathrm{\Psi }}_n^{*}({\mathbf{x}}^{\prime };\omega )}{\omega -{ϵ}_n(\omega )}}}

where {\mathrm{\Psi }}_n(\mathbf{x};\omega ) and {ϵ}_n(\omega ) are the eigen-solutions of the following equation

(15)[h_{\mathrm{0}}(\mathbf{x})+V_H(\mathbf{x})]{\mathrm{\Psi }}_n(\mathbf{x};\omega )+{\displaystyle \int d{\mathbf{x}}^{\prime }{\displaystyle \sum (\mathbf{x},{\mathbf{x}}^{\prime };\omega ){\mathrm{\Psi }}_n({\mathbf{x}}^{\prime };\omega )}}=E_n(\omega ){\mathrm{\Psi }}_n(\mathbf{x};\omega )

By definition, ω is real, but formally the domain of ω can be extended to the whole complex plane by analytic continuation, so that the poles of G(x, x′, ω), with \omega ={ϵ}_n(\omega ), correspond to quasi-particle energies. Since QP energies and wave functions are the quantities of main interest, Eq. (15) can be rewritten in the following form, often denoted as the QP equation,

(16)[h_{\mathrm{0}}(\mathbf{x})+V_H(\mathbf{x})]{\mathrm{\Psi }}_n(\mathbf{x})+{\displaystyle \int d{\mathbf{x}}^{\prime }{\displaystyle \sum (\mathbf{x},{\mathbf{x}}^{\prime };{ϵ}_n){\mathrm{\Psi }}_n({\mathbf{x}}^{\prime })}}={ϵ}_n{\mathrm{\Psi }}_n(\mathbf{x})

We note that the term “quasi-particle” here is used in a generalized sense.

Hedin’s equations and the GW approximation

The self-energy operator is itself a highly complicated quantity, and requires further approximations in practice. There are mainly two approaches to approximate Σ in a computationally accessible manner: 1) the Feynman’s diagrammatic approach based on Wick’s theorem [8], and 2) the functional derivative approach [9,10,23]. The former, for example, is used in the so-called electron propagator theory for molecular systems [26], where the screening among electrons is weak, and a many-body perturbation theory with respect to the bare Coulomb interaction up to the second or third order often gives very accurate descriptions as recently demonstrated by Ortiz and coworkers (e.g., [27] and references therein). For extended systems with stronger screening, the second approach, in the form of a set of coupled equations (Hedin’s equations), is more appropriate, and is the subject of this work.

Hedin’s equations, as first formulated in a complete form by Hedin [10] using the functional derivative technique, link the self-energy with the dynamically screened Coulomb interaction W, which, when calculated perturbatively, leads to a many-body perturbation expansion of the self-energy with respect to W instead of v. The derivation of Hedin’s equations can be found in Refs. [9,23,24], and here we just give the final expressions with short explanations on the important quantities that are involved in the equations.

Using the short-hand notation \mathrm{1}\equiv ({\mathbf{x}}_{\mathrm{1}},t_{\mathrm{1}}), Hedin’s equations read

{\displaystyle \sum (\mathrm{1},\mathrm{2})=i{\displaystyle \int d(\mathrm{34})G(\mathrm{1},\mathrm{3})W(\mathrm{4},\mathrm{1})\mathrm{\Gamma }(\mathrm{3},\mathrm{2},\mathrm{4})}}

\mathit{\Gamma }\left(1,2,3\right)=\delta \left(1,2\right)\delta \left(2,3\right)+{\displaystyle \int d\left(4567\right){\displaystyle \frac{\delta {\displaystyle \sum \left(1,2\right)}}{\delta G\left(4,5\right)}}G\left(4,6\right)}\times G\left(7,5\right)\mathrm{\Gamma }\left(6,7,3\right)

W(\mathrm{1},\mathrm{2})=v(\mathrm{1},\mathrm{2})+{\displaystyle \int d(\mathrm{34})v(\mathrm{1},\mathrm{3})P(\mathrm{3},\mathrm{4})W(\mathrm{4},\mathrm{2})}

(17)P\left(1,2\right)=-i{\displaystyle \int d\left(34\right)G\left(1,3\right)\mathit{\Gamma }\left(3,4,2\right)G\left(4,1^{+}\right)}

which, together with the equation of motion for Green’s function (Eq. (10)), form a close set of equations. In Eq. (17), P is the polarization function,

(18)P(\mathrm{1},\mathrm{2})\equiv {\displaystyle \frac{\delta \rho (\mathrm{1})}{\delta V(\mathrm{2})}}

which describes the response of the electron density (ρ) with respect to the total effective potential V = V_H + ø, with ø being the external perturbation. W is the screened Coulomb interaction defined as

(19)W(\mathrm{1},\mathrm{2})={\displaystyle \int d(\mathrm{3}){ϵ}^{-\mathrm{1}}(\mathrm{1},\mathrm{3})v(\mathrm{3},\mathrm{2})}

with the inverse dielectric function ϵ^-1 defined as

(20){ϵ}^{-\mathrm{1}}(\mathrm{1},\mathrm{2})\equiv {\displaystyle \frac{\delta V(\mathrm{1})}{\delta \varphi (\mathrm{2})}}

The dielectric function ϵ (1, 2) is related to the polarization function by

(21)ϵ(\mathrm{1},\mathrm{2})=\delta (\mathrm{1},\mathrm{2})-{\displaystyle \int d(\mathrm{3})v(\mathrm{1},\mathrm{3})P(\mathrm{3},\mathrm{2})}

Г (1, 2, 3) is the vertex function defined as

(22)\mathit{\Gamma }\left(1,2,3\right)\equiv -{\displaystyle \frac{\delta G^{-1}\left(1,2\right)}{\delta V\left(3\right)}}.

Hedin’s equations are a set of complicated integro-differential equations, which are not solvable even for the simplest systems like the homogeneous electron gas [9]. The main use of Hedin’s equations is to serve as the starting point for a many-body perturbation theory in terms of W. Since the screening, the most important electron-electron interaction effect for many extended systems is a built-in feature of W, and even a low-order expansion of the self-energy with respect to W is hopeful to give reasonable description. In practice, the most widely used approach to solve Hedin’s equations is the so-called GW approximation, which can be obtained by taking the zero-order approximation for the vertex function \mathit{\Gamma }\left(1,2,3\right)\simeq \delta \left(1,2\right)\delta \left(2,3\right). In this case, the polarization function is simplified to be the product of two Green’s functions,

(23)P(\mathrm{1},\mathrm{2})=-iG(\mathrm{1},\mathrm{2})G(\mathrm{2},{\mathrm{1}}^{+})

which is often called random-phase approximation (RPA) for the historical reason. The self-energy now becomes a simple product of the Green’s function G and screened Coulomb interaction W

(24){\displaystyle \sum (\mathrm{1},\mathrm{2})}=iG(\mathrm{1},{\mathrm{2}}^{+})W(\mathrm{2},\mathrm{1})

The positive infinitesimal in the equation above is needed for the bare exchange self-energy as a result of the Feynman’s rules [8], G({\mathbf{x}}_{\mathrm{1}}{t}_{\mathrm{1}};{\mathbf{x}}_{\mathrm{2}}{t}_{\mathrm{2}}=t_{\mathrm{1}})\equiv G({\mathbf{x}}_{\mathrm{1}}{t}_{\mathrm{1}};{\mathbf{x}}_{\mathrm{2}}{t}_{\mathrm{1}}^{+}) [11,24]. When represented in the frequency domain, the GW self-energy reads

(25){\displaystyle \sum (\mathbf{x},{\mathbf{x}}^{\prime };\omega )={\displaystyle \frac{i}{\mathrm{2}\pi }}{\displaystyle \int d{\omega }^{\prime }{e}^{i{\omega }^{\prime }\delta }G(\mathbf{x},{\mathbf{x}}^{\prime };\omega +{\omega }^{\prime })}W({\mathbf{x}}^{\prime },\mathbf{x};{\omega }^{\prime })}

GW in practice: Numerical techniques

The G₀W₀ approach

Even within the GW approximation Hedin’s equations are still highly complicated, as illustrated in Fig. 1, requiring a self-consistent calculation of the self-energy from QP energies and wavefuctions as solutions of Dyson’s equation. The latter is mathematically cumbersome to tackle due to the nonhermiticity of the self-energy. In practice, further approximations are often introduced. The most widely used implementation of the GW approximation is the so-called G₀W₀ or one-shot GW approach, in which both G and W are calculated based on eigen-energies є_n and eigen-functions {\psi }_n of some non-interacting reference system,

(26)\left[-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\nabla }^{\mathrm{2}}+V_{ext}(\mathbf{x})+V_H(\mathbf{x})+V_{xc}(\mathbf{x})\right]{\psi }_n(\mathbf{x})={є}_n{\psi }_n(\mathbf{x})

where V_xc(x) is the Kohn-Sham exchange-correlation potential, often in LDA or GGA. G₀ now reads

(27)G_{\mathrm{0}}(\mathbf{x},{\mathbf{x}}^{\prime };\omega )={\displaystyle \sum _n{\displaystyle \frac{{\psi }_n(\mathbf{x}){\psi }_n^{*}({\mathbf{x}}^{\prime })}{\omega -{\widetilde{є}}_n}}}

where {\widetilde{є}}_n\equiv {є}_n+i\eta \mathrm{sgn}({є}_F-{є}_n). The polarization function can then be calculated as

(28)$$\begin{gathered} P_{\mathrm{0}}(\mathbf{x},{\mathbf{x}}^{\prime };\omega )=-{\displaystyle \frac{i}{\mathrm{2}\pi }}{\displaystyle \int G_{\mathrm{0}}(\mathbf{x},{\mathbf{x}}^{\prime };\omega +{\omega }^{\prime })G_{\mathrm{0}}({\mathbf{x}}^{\prime },\mathbf{x};{\omega }^{\prime })d{\omega }^{\prime }} \\ ={\displaystyle \sum _{n,m}{f}_n(\mathrm{1}-f_m){\psi }_n(\mathbf{x}){\psi }_m^{*}(\mathbf{x}){\psi }_n^{*}({\mathbf{x}}^{\prime }){\psi }_m({\mathbf{x}}^{\prime })}\times {\displaystyle \frac{\mathrm{1}}{\omega -{є}_m+{є}_n+i\eta }}-{\displaystyle \frac{\mathrm{1}}{\omega -{є}_m-{є}_n-i\eta }} \\ \equiv {\displaystyle \sum _{n,m}{F}_{nm}(\omega ){\psi }_n(\mathbf{x}){\psi }_m^{*}(\mathbf{x}){\psi }_n^{*}({\mathbf{x}}^{\prime }){\psi }_m({\mathbf{x}}^{\prime })} \end{gathered}$$

where f_n is the occupation number of the n-th state, and F_nm(ω) denotes the occupation and frequency dependent factor. Using P₀, the RPA screened Coulomb interaction can be obtained

(29)\begin{array}{l}{W}_{\mathrm{0}}(\mathbf{x},{\mathbf{x}}^{\prime };\omega )={\displaystyle \int d{\mathbf{x}}^{″}{ϵ}^{-\mathrm{1}}(\mathbf{x},{\mathbf{x}}^{″};\omega )v({\mathbf{x}}^{″},{\mathbf{x}}^{\prime })}\hfill \\ ϵ(\mathbf{x},{\mathbf{x}}^{\prime };\omega )=\delta (\mathbf{x},{\mathbf{x}}^{\prime })-{\displaystyle \int d{\mathbf{x}}^{″}v(\mathbf{x},{\mathbf{x}}^{″})}{P}_{\mathrm{0}}({\mathbf{x}}^{″},{\mathbf{x}}^{\prime };\omega )\hfill \end{array}

It is practically common and convenient to decompose the self-energy into exchange and correlation terms by defining

(30)W_{\mathrm{0}}^c(\mathbf{x},{\mathbf{x}}^{\prime };\omega )=W_{\mathrm{0}}(\mathbf{x},{\mathbf{x}}^{\prime };\omega )-v(\mathbf{x},{\mathbf{x}}^{\prime })

The exchange part, which is just given by the Hartree-Fock exchange potential, reads

(31){\displaystyle {\sum }_{}^{\mathrm{x}}(\mathbf{x},{\mathbf{x}}^{\prime })={\displaystyle \frac{i}{\mathrm{2}\pi }}}{\displaystyle \int G_{\mathrm{0}}(\mathbf{x},{\mathbf{x}}^{\prime };{\omega }^{\prime })v({\mathbf{x}}^{\prime },\mathbf{x})e^{i{\omega }^{\prime }\eta d{\omega }^{\prime }}=-{\displaystyle \sum _n{f}_n{\psi }_n(\mathbf{x})v({\mathbf{x}}^{\prime },\mathbf{x}){\psi }_n^{*}({\mathbf{x}}^{\prime })}}

The GW correlation self-energy is obtained from the frequency integral

(32){\displaystyle {\sum }_{}^c(\mathbf{x},{\mathbf{x}}^{\prime };\omega )={\displaystyle \frac{i}{\mathrm{2}\pi }}{\displaystyle \int G_{\mathrm{0}}(\mathbf{x},{\mathbf{x}}^{\prime };\omega +{\omega }^{\prime })W_{\mathrm{0}}^c({\mathbf{x}}^{\prime },\mathbf{x};{\omega }^{\prime })d{\omega }^{\prime }.}}

The QP energies ϵ_n are calculated by the first order perturbation theory, treating \delta {\displaystyle \sum \equiv {\displaystyle \sum -V_{xc}}} as the perturbation

(33){ϵ}_n={є}_n+Z_n({є}_n)\Re \langle {\psi }_n|\sum ({є}_n)-V_{xc}|{\psi }_n\rangle

where Z_n is the QP renormalization factor,

(34)Z_n(E)={\left[\mathrm{1}-{\left({\displaystyle \frac{\partial }{\partial \omega }}\Re \langle {\psi }_n|\sum (\omega )|{\psi }_n\rangle \right)}_{\omega =E}\right]}^{-\mathrm{1}}

accounting for the frequency dependency of Σ. For sp-semiconductors it has been demonstrated that further improvement can be obtained by introducing partial selfconsistency in the so-called GW₀ approach, in which the energies used for the calculation of the Green’s function are updated by QP energies with fixed W₀ [22,28]. The GW₀ can be implemented with little computational overhead if the intermediate quantities used for Σ^c are stored during the calculation<FootNote>

Jiang, H.; Gomez-Abal, R. I.; Li, X.; Meisenbichler, C.; Ambrosch-Draxl, C.; Scheffler, M. “Fhi-gap: A green-function code based on augmented planewaves”

</FootNote>.

The G₀W₀ equations in the matrix form

To put the GW approach in action, the first step is to convert the equations above into the matrix form. For that purpose we need a set of basis functions that are able to represent the products of any two Kohn-Sham wave functions accurately. Since we are mainly concerned with extended systems, we will include explicitly the dependence on wave-vectors (k) that characterize the translational invariance of the system in the equations from now on unless stated otherwise. We use Ω to denote the volume of the unit cell. The difference between two k vectors is denoted as q. We first start with a general basis set, denoted as \{{\chi }_i^q(\mathbf{x})\}. More technical details on using a particular basis will be discussed in the next subsection. To simplify the notation, we assume \{{\chi }_i^q(\mathbf{x})\} with the same q are orthonormal, although non-orthonormal basis functions can also be used. In addition, we assume that the system under study is spin-unpolarized so that the spin index is dropped.

The most important requirement for the basis functions used in G₀W₀ is to represent the products of two Kohn-Sham wave functions accurately

(35){\psi }_{n\mathbf{k}}(\mathbf{r}){\psi }_{m\mathbf{k}-\mathbf{q}}^{*}(\mathbf{r})={\displaystyle \sum _i{M}_{nm}^i(\mathbf{k},\mathbf{q})}{\chi }_i^{\mathbf{q}}(\mathbf{r})

where M_{nm}^i(\mathbf{k},\mathbf{q}) are expansion coefficients given by

(36)M_{nm}^i(\mathbf{k},\mathbf{q})\equiv {{\int }_V[{\chi }_i^{\mathbf{q}}(\mathbf{r}){\psi }_{n\mathbf{k}-\mathbf{q}}(\mathbf{r})]}^{*}{\psi }_{n\mathbf{k}}(\mathbf{r})d^3\mathbf{r}

V is the total volume of the system, which is related to the volume of the unit cell by V = N_kΩ with N_k being the number of k-points. We denote this basis the product basis to distinguish it from the basis that is used to expand single-particle wave functions. In quantum chemistry, it is also often called the “auxiliary” basis. Any two-point function, g(r, r′), such as the polarizability P₀(r, r′; ω), the bare Coulomb interaction v(r, r′), the dielectric function ϵ (r, r′; ω), and the screened Coulomb interaction W(r, r′; ω), will be expanded in terms of the product basis:

(37)g(\mathbf{r},{\mathbf{r}}^{\mathit{\text{'}}})={\displaystyle \sum _{\mathbf{q}}^{\mathrm{B}\mathrm{Z}}{\sum _{i,j}{g}_{ij}(\mathbf{q}){\chi }_i^{\mathbf{q}}(\mathbf{r})[{\chi }_j^{\mathrm{q}}({\mathbf{r}}^{\mathit{\text{'}}})]}^{*}}

and the expansion coefficients are determined by

(38)g_{ij}(\mathbf{q})\equiv {\displaystyle {\int }_V{\displaystyle {\int }_V{[{\chi }_i^{\mathbf{q}}(\mathbf{r})]}^{*}g(\mathbf{r},{\mathbf{r}}^{\prime }){\chi }_j^{\mathbf{q}}({\mathbf{r}}^{\prime })}}

The diagonal elements of the bare (Fock) exchange self-energy can be written in the product basis as

(39)$$\begin{gathered} {\displaystyle {\sum }_{n\mathbf{k}}^x}\equiv {\displaystyle {\int }_{V}d\mathbf{r}{\displaystyle {\int }_V{\mathbf{r}}^{\prime }{\left[{\psi }_{n\mathbf{k}}\left(\mathbf{r}\right)\right]}^{*}{\displaystyle {\sum }_{}^x\left(\mathbf{r},{\mathbf{r}}^{\prime }\right){\psi }_{n\mathbf{k}}\left({\mathbf{r}}^{\prime }\right)}}} \\ =-{\displaystyle \sum _{\mathbf{q}}^{BZ}{\displaystyle \sum _{i,j}{v}_{ij}\left(\mathbf{q}\right){{\displaystyle \sum _m^{occ}\left[M_{nm}^i\left(\mathbf{k},\mathbf{q}\right)\right]}}^{*}{M}_{nm}^i\left(\mathbf{k},\mathbf{q}\right)}} \end{gathered}$$

The matrix elements of the polarizability (Eq. (28)) can be written as

(40)$$\begin{gathered} P_{ij}\left(\mathbf{q},\omega \right)\equiv {\displaystyle {\int }_V{\displaystyle {\int }_V{\left[{\chi }_i^{\mathbf{q}}\left(\mathbf{r}\right)\right]}^{*}P\left(\mathbf{r},{\mathbf{r}}^{\prime };\omega \right){\chi }_j^{\mathbf{q}}\left({\mathbf{r}}^{\prime }\right)d\mathbf{r}d{\mathbf{r}}^{\prime }}} \\ =2{\displaystyle \sum _k^{BZ}{{\displaystyle \sum _{n,m}{F}_{nm}\left(\mathbf{k},\mathbf{q};\omega \right)M_{nm}^i\left(\mathbf{k},\mathbf{q}\right)\left[M_{nm}^j\left(\mathbf{k},\mathbf{q}\right)\right]}}^{*}} \end{gathered}$$

where the factor of 2 is from the spin degeneracy, and

(41)F_{nm}(\mathbf{k},\mathbf{q};\omega )\equiv f_{n\mathbf{k}}(\mathrm{1}-f_{m\mathbf{k}-q})\left[{\displaystyle \frac{\mathrm{1}}{\omega -{\omega }_{n\mathbf{k},m\mathbf{k}-q}+i\eta }}-{\displaystyle \frac{\mathrm{1}}{\omega +{\omega }_{n\mathbf{k},m\mathbf{k}-q}-i\eta }}\right]

with {\omega }_{n\mathbf{k},m\mathbf{k}-q}\equiv {є}_{m\mathbf{k}-q}-{є}_{nk}. Due to the singularity of the bare Coulomb interaction in reciprocal space as q goes to zero, the dielectric function at q → 0 has to be treated carefully. Mathematically, it is more convenient to use the symmetrized form of the dielectric function [29], which, in the matrix form, can be obtained from P_ij(q, ω) by

(42){ϵ}_{ij}(\mathbf{q},\omega )={\delta }_{ij}-{\displaystyle \sum _{lm}{v}_{il}^{\frac{\mathrm{1}}{\mathrm{2}}}(\mathbf{q})P_{lm}(\mathbf{q},\omega )v{v}_{mj}^{\frac{\mathrm{1}}{\mathrm{2}}}(\mathbf{q})}

The correlation term of the screened Coulomb interaction can then be calculated through

(43)W_{ij}^c(\mathbf{q},\omega )={\displaystyle \sum _{lm}{v}_{il}^{\frac{\mathrm{1}}{\mathrm{2}}}(\mathbf{q})[{ϵ}_{lm}^{-\mathrm{1}}(\mathbf{q},\omega )-{\delta }_{lm}]v_{mj}^{\frac{\mathrm{1}}{\mathrm{2}}}.}

The diagonal matrix elements of the correlation self-energy can be written as

(44)$$\begin{gathered} {\displaystyle {\sum }_{n\mathbf{k}}^{\mathrm{c}}(\omega )}\equiv {{\int }_V\mathrm{d}\mathbf{r}{\displaystyle {\int }_V{\mathbf{r}}^{\mathit{\text{'}}}[{\psi }_{n\mathbf{k}}(\mathbf{r})]}}^{*}{\displaystyle {\sum }_{}^c(\mathbf{r},{\mathbf{r}}^{\mathit{\text{'}}};\omega ){\psi }_{n\mathbf{k}}({\mathbf{r}}^{\mathit{\text{'}}})} \\ ={\displaystyle \sum _{\mathbf{q}}^{BZ}{\displaystyle \sum _m{\displaystyle \frac{\mathrm{i}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle {\int }_{-\infty }^{\infty }\mathrm{d}{\omega }^{\mathit{\text{'}}}{\displaystyle \frac{X_{nm}(\mathbf{k},\mathbf{q};{\omega }^{\mathit{\text{'}}})}{\omega +{\omega }^{\mathit{\text{'}}}-{\widetilde{\mathit{є}}}_{m\mathbf{k}-\mathbf{q}}}}}}} \end{gathered}$$

where the auxiliary quantity X_nm(k, q; ω) is defined as

(45)X_{nm}(\mathbf{k},\mathbf{q};\omega )\equiv {\displaystyle \sum _{ij}{[M_{nm}^i(\mathbf{k},\mathbf{q})]}^{*}{W}_{ij}^c(\mathbf{q},\omega )M_{nm}^i(\mathbf{k},\mathbf{q})}

We can see from the equations above that the main ingredients in the numerical implementation of G₀W₀ are 1) the bare Coulomb matrix, 2) the overlap integrals between the product basis functions and Kohn-Sham wave functions products, and 3) the summation of states in the calculation of P₀ and the self-energy. The treatment of 1) and 2) depends strongly on the nature of the product basis, and 3) is usually the most time-consuming part to calculate.

GW in different basis representations

The first “first-principles” implementation of the GW method by Hybertsen and Louie [12], and Godby et al. [30] used the planewave (PW) representation (also called the reciprocal space representation by many authors) in combination with the pseudopotential (PP) approximation. As we will see soon, the greatest advantage of the planewave representation is the simpleness in terms of implementation. In particular, since the planewave is the most popular basis for Kohn-Sham DFT for extended systems, the extension from DFT to GW is relatively straightforward. In addition, since the quality of the PW basis is controlled by a single parameter, the energy cutoff, the accuracy of the representation can be systematically monitored. Nevertheless, the limitation of the PW representation is also obvious. For systems with d- or f- states not very far way from the Fermi level, it is often necessary to treat these states as valence, but since these states are much more localized than sp states, the number of the planewaves required to represent them is usually quite large, leading to very heavy calculations. For systems with loose structure or a lot of open space such as weakly bonded molecular crystals and surfaces, the PW representation is also not very efficient. There is another more technical drawback of PW: Since PWs are global basis functions, the matrices represented by PWs are usually far from sparse, which makes parallelization of the corresponding GW code much more difficult that those based on local basis functions.

To overcome the limitation of the PW basis functions, a lot of efforts have been invested in the past two decades to develop new techniques that are either more accurate or can treat larger systems than that the PW representation can afford. In the following, we present an overview on the implementations of GW with different basis functions.

The planewaves representation

When using the planewaves

(46){\chi }_i^q(\mathbf{r})\to {\chi }_G^q(r)\equiv {\displaystyle \frac{\mathrm{1}}{\sqrt{V}}}\exp [i(\mathbf{q}+\mathbf{G})·\mathbf{r}]

as the product basis, then the integrals M_{nm}^i(\mathbf{k},\mathbf{q}) become simply the Fourier transform of {\psi }_{nk}(\mathbf{r}){\psi }_{mk-q}^{*}(\mathbf{r}). In addition, since {\psi }_{n\mathbf{k}}(\mathbf{r}) is also represented by plane waves,

(47){\psi }_{n\mathbf{k}}={\displaystyle \sum _{\mathbf{G}}{c}_{n\mathbf{k}};\mathbf{G}{\chi }_G^k(\mathbf{r})}

one obtains

(48)M_{nm}^G(\mathbf{k},\mathbf{q})=V^{-\mathrm{1}/\mathrm{2}}{\displaystyle \sum _{{\mathbf{G}}^{\prime }}{C}_{n\mathbf{k};{\mathbf{G}}^{\prime }}{C}_{m\mathbf{k}-\mathbf{q};G^{\prime }-G}^{*}}

An even more favorable feature of the PW basis is that the bare Coulomb matrix is diagonal in the PW representation

(49)v_{{GG}^{\prime }}(\mathbf{q})={\displaystyle \frac{\mathrm{4}\pi }{{\left|\mathbf{q}+\mathbf{G}\right|}^{\mathrm{2}}}}{\delta }_{{G,G}^{\prime }}

The symmetrized dielectric function now reads

(50){ϵ}_{{GG}^{\prime }}(\mathbf{q},\omega )={\delta }_{{GG}^{\prime }}-{\displaystyle \frac{\mathrm{4}\pi }{\left|\mathbf{q}+\mathbf{G}\right|\left|\mathbf{q}+{\mathbf{G}}^{\prime }\right|}}{P}_{{GG}^{\prime }}(q,\omega ).

When q is at the Г point, i.e., q = 0, the second term in Eq. (50) is divergent when G = 0 and/or G′ = 0. However, this singularity can be easily removed by expanding {\mathrm{P}}_{\mathrm{G}{\mathrm{G}}^{\prime }}\left(\mathbf{q},\omega \right) (q, ω) around q = 0 when G = 0 and/or G′ = 0, which gives the lowest order term \propto {|q|}^{\mathrm{2}} (when both G = 0 and G′ = 0) or \propto |q| (G = 0 or G′ = 0).

The space-time approach

A variant of the PW representation is the space-time approach developed by R. Godby and coworkers [31,32] The basic ideas are as follow.

1. Since both the RPA polarization function and the GW self-energy are simple products in the space and time domain,

(51)\begin{array}{l}{P}_{\mathrm{0}}(\mathbf{r},{\mathbf{r}}^{\prime };\tau )=-i{G}_{\mathrm{0}}(\mathbf{r},{\mathbf{r}}^{\prime };\tau )G_{\mathrm{0}}({\mathbf{r}}^{\prime },\mathbf{r};-\tau )\hfill \\ {\displaystyle \sum (\mathbf{r},{\mathbf{r}}^{\prime };\tau )=iG(\mathbf{r},{\mathbf{r}}^{\prime };\tau )W({\mathbf{r}}^{\prime },\mathbf{r},;\tau )}\hfill \end{array}

it is computationally more efficient to calculate them directly in the space-time representation instead of in the planewave representation.

2. On the other hand, the evaluation of the dielectric function ϵ and the screened Coulomb interaction (W) involves the expensive integration in the real space. It is more convenient to calculate them in the reciprocal space to take advantage of the diagonality of the bare Coulomb interaction. The space-time representation and the planewaves-frequency representation can be converted to each other by the fast-Fourier transform technique.

3. To avoid tackling the complicated structures that Σ, W and G have along the real frequency/time, it is more efficient to construct these functions first on the imaginary time/frequency. Self-energies along the real frequency axis can be obtained by the analytic continuation technique, thanks to the analyticity of these functions in the complex frequency domain.

It is important to bear in mind that since in the space-time approach, quantities like G, P, W and Σ need to be represented in both real-space (usually on a uniform grid) and the reciprocal space, the effective potential of the system under study has to be quite smooth, so that it is necessary to use the pseudo-potential approximation, as in the standard planewaves approach.

The local basis representation

A completely different framework for the implementation of the GW method is to use local atomic-like basis. The latter can be analytic Gaussian-type orbitals (GTO) [33,34] that is the de facto standard in molecular quantum chemistry community [35], or numerical atomic orbitals [36]. The general matrix formulation in the preceding subsection is based on an orthogonal basis representation. When using the local atomic-like (numerical or analytic) basis functions, basis functions are in general not orthogonal, for which particular attention is needed to obtain a robust and efficient implementation.

Local basis functions written in the Bloch form reads

(52){\chi }_{\alpha }^q(\mathbf{r})={\displaystyle \frac{\mathrm{1}}{N_c^{\mathrm{1}/\mathrm{1}}}}{\displaystyle \sum _R{e}^{iq·(R+{\tau }_{\alpha })}{\varphi }_{\alpha }(\mathbf{r}-\mathbf{R}-{\tau }_{\alpha })}

where {\varphi }_{\alpha }(\mathbf{r}) are atomic-like functions centered at τ_a. We use the index α to denote both the positions of the atoms and the atomic quantum numbers that characterize the atomic basis functions. Now we consider a general two-point function g(r, r′), and define the matrix elements

(53){[g]}_{\alpha \beta }(\mathbf{q})\equiv {\displaystyle {\int }_{V}d\mathbf{r}{\displaystyle \underset{V}{\int }d{\mathbf{r}}^{\mathbf{\text{'}}}{\chi }_{\alpha }^{\mathbf{q}*}(\mathbf{r})g(\mathbf{r},{\mathbf{r}}^{\prime }){\chi }_{\beta }^{\mathbf{q}}({\mathbf{r}}^{\prime })}}

Then g(r, r′) can be expanded by the local basis functions as

(54)g(\mathbf{r},{\mathbf{r}}^{\prime })={\displaystyle \sum _q{\displaystyle \sum _{\alpha ,\beta }{\chi }_{\alpha }^q(\mathbf{r}){\langle g\rangle }_{\alpha \beta }(\mathbf{q}){\chi }_{\beta }^{q^{*}}(r)}}

The matrix \langle g\rangle is related to [g] by (in the matrix form)

(55)\langle g\rangle (\mathbf{q})={\mathbf{S}}^{-\mathrm{1}}(\mathbf{q})[g](\mathbf{q}){\mathbf{S}}^{-\mathrm{1}}(\mathbf{q})

where S(q) is the overlap matrix

(56)S_{\alpha \beta }(\mathbf{q})\equiv {\displaystyle \underset{V}{\int }d\mathbf{r}{\chi }_{\alpha }^{q^{*}}(\mathbf{r}){\chi }_{\beta }^q(\mathbf{r})}

To obtain the GW equations in the non-orthogonal basis representation, we can introduce the orthogonalized basis corresponding to {\chi }_{\alpha }^{\mathbf{q}}(\mathbf{r})

(57)|{\tilde{\chi }}_i^q\rangle ={\displaystyle \sum _{\alpha }|{\chi }_{\alpha }^q\rangle S_{\alpha i}^{-\mathrm{1}/\mathrm{2}}(\mathbf{q})}

The orthogonal and non-orthogonal representations are related by

(58)[\tilde{g}](\mathbf{q})={\mathbf{S}}^{-\mathrm{1}/\mathrm{2}}(\mathbf{q})[g](\mathbf{q}){\mathbf{S}}^{-\mathrm{1}/\mathrm{2}}(\mathbf{q})

Using the relation above, one can easily convert the equations in Section 3.2 into their counter-parts in the non-orthogonal basis representation.

A relatively simpler scheme to handle the nonorthogonality of atomic basis functions is to use the orthogonalized basis from the very beginning [34,36]. The basic ingredients such as the bare Coulomb matrix v_ij and the product expansion matrix M_{nm}^i are first calculated in the original non-orthogonal basis, and then are transformed into the orthogonalized form. Then all subsequent treatments are essentially same as that in the orthogonal basis representation. In addition, by using this pre-orthogonalized basis, the linear-dependency of the non-orthogonal basis functions, usually occurring between basis functions centered on different atoms, can be removed from the very beginning by dropping those eigenvectors of the overlap matrix that have vanishing eigen-values [34].

The mixed basis representation

One of the main concerns over using the planewave representation or the space-time approach is that the pseudopotential approximation is necessary. In principles, the error introduced by the use of the PP approximation can be controlled if all relevant atomic states are treated as “valence,” but in practice the situations are more complicated. In KS-DFT, the accuracy of using the PP relies on the validity of two approximations: 1) the frozen-core approximation, i.e., that the states treated as “core” are chemically inert, and 2) the linear approximation for the exchange-correlation potential as the functional of the electron density. In the GW method, there is one more aspect of using the PP approximation: since the GW self-energy depends on wave functions instead of electron density, the use of the pseudo-wave functions can also have significant influences on the accuracy of the GW [N-39]: 1) the use of the PP approximation has much stronger effects on the accuracy of the GW results than it does on KSDFT calculations; and 2) for many elements, it is often necessary to use the pseudopotentials that treat shallow semi-core states as “valence,” and therefore they are much harder than those used in KS-DFT.

To avoid the difficulty of using the PP approach, and in the meanwhile maintain the advantage of the planewave basis that the accuracy can be systematically improved, the full-potential augmented basis approaches have been developed for the implementation of GW [40-43] (see footnote on Page 6). We will use the full-potential linearized-planewaves (FP-LAPW) as an example to describe the main features of this type of approaches (see footnote on Page 6). In the FP-LAPW approach, the space in the unit cell is partitioned into non-overlapping muffin-tin (MT) spheres, centered around each atom (indexed by α, positioned at τ_a in the unit cell), and the interstitial (IS) region. The LAPW basis is given in the IS region as plane-waves

(59){\varphi }_G^k(\mathbf{r})={\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{\Omega }}}}{e}^{i(k+G)·r},\mathbf{r}\in IS

In the MT spheres, it is represented by atomic-like wave functions

(60)

where {\mathbf{r}}^{\alpha }\equiv \mathbf{r}-{\tau }_{\alpha } and R_{MT}^{\alpha } is the radius of α-th MT sphere. ul({\mathbf{r}}^{\alpha }) are the solutions of the radial Schrödinger equation in the spherical potential of the respective MT sphere, taken at a l-dependent reference energy E_l, and {\dot{u}}_l({\mathbf{r}}^{\alpha }) is the first order derivative of ul({\mathbf{r}}^{\alpha }) with respect to the energy at E_l. The augmentation coefficients, A_{\alpha lm}(\mathbf{k}+\mathbf{G}) and B_{\alpha lm}(\mathbf{k}+\mathbf{G}), are determined from the continuity of the basis functions and their first derivatives at the MT sphere boundary.

To implement GW in the FP-LAPW framework, the product basis with similar mixed features should be used to represent the products of KS wave functions accurately, hence termed as the mixed basis (MB). The MB functions can be constructed in the following way: Inside the MT sphere of atom α, the basis functions are atomic like, similar to the product basis originally proposed by Aryasetiawan and Gunnarsson [40].

(61){\chi }_i^{\mathbf{k}}(\mathbf{r})={\displaystyle \frac{\mathrm{1}}{N_c^{\mathrm{1}/\mathrm{2}}}}{\displaystyle \sum _R{e}^{ik\times (R+{\tau }_{\alpha })}}{v}_{\alpha NL}(\left|\mathbf{r}-{\tau }_{\alpha }-\mathbf{R}\right|)Y_{LM}(\widehat{\mathbf{r}})

The radial functions v_{\alpha NL}({\mathbf{r}}^{\alpha }) are often constructed from the products of the radial wave functions used in the LAPW basis (see Eq. (60)). Since the number of the latter is quite big, and in addition, they are not fully linearly independent, the following procedure is often used to construct an optimal set of orthonormal radial functions:

● To reduce the number of product functions, only u_{al}({\mathbf{r}}^{\alpha })’s with l\le l_{\max }^{MB} are considered; in addition, {\dot{u}}_{\alpha l}(\mathbf{r})’s are not taken into account as they are typically one order smaller than u_{\alpha l}(\mathbf{r}) [44,45].

● For each L, all the products of two radial functions u_{\alpha l}({\mathbf{r}}^{\alpha })u_{\alpha l^{\prime }}({\mathbf{r}}^{\alpha }) that fulfill the triangular condition \left|l-l^{\prime }\right|\le L\le l+l^{\prime } are considered.

● The overlap matrix between this set of product radial functions

(62)\mathfrak{O}{u}^{\prime };l_{\mathrm{1}}{l}_{\mathrm{1}}^{\prime }={\displaystyle \underset{\mathrm{0}}{\overset{R_{MT}^{\alpha }}{\int }}{u}_{\alpha l}(\mathbf{r})u_{\alpha l^{\prime }}(\mathbf{r})u_{\alpha l_{\mathrm{1}}}(\mathbf{r})u_{\alpha l_{\mathrm{1}}^{\prime }}(\mathbf{r}){\mathbf{r}}^{\mathrm{2}}d\mathbf{r}}

is diagonalized, yielding the corresponding set of eigenvalues, {\lambda }_N^{MB}, and eigenvectors, \{c_{u^{\prime },N}\}.

● Eigenvectors corresponding to eigenvalues smaller than a certain tolerance, {\lambda }_{\min }^{MB} (typically about 10^-4), are assumed to be linearly dependent and discarded [38]; The remaining eigenvectors, after normalization, form the radial functions

(63)v_{\alpha NL}({\mathbf{r}}^{\alpha })={\displaystyle \sum _{u^{\prime }}c{u}^{\prime },N^{ul}({\mathbf{r}}^{\alpha })u{l}^{\prime }({\mathbf{r}}^{\alpha })}

In the interstitial region, the product basis functions are constructed as ortho-normalized interstitial planewaves (IPW’s),

(64)P_i^{\mathbf{q}}(\mathbf{r})\equiv {\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{\Omega }}}}{\displaystyle \sum _{\mathbf{G}}{\displaystyle \frac{S_{\mathbf{G}i}}{{\lambda }_i^{IS}}}{e}^{i(\mathbf{G}+\mathbf{q})·\mathbf{r}}{\theta }_{IS}(\mathbf{r})}

where

(65){\theta }_{IS}(\mathbf{r})=\{\begin{array}{ll}\mathrm{1}\hfill & \mathbf{r}\in IS,\hfill \\ \mathrm{0}\hfill & otherwise\mathrm{.}\hfill \end{array}

S_{\mathbf{G}i} and {\lambda }_i^{IS} are eigenvectors and eigenvalues of the overlap matrix between IPW’s

(66){\mathfrak{O}}_{\mathbf{G}{\mathbf{G}}^{\mathit{\text{'}}}}={\displaystyle \frac{1}{\mathrm{\Omega }}}{\displaystyle {\int }_{\mathrm{\Omega }}{\theta }_{\mathrm{I}\mathrm{S}}\left(\mathbf{r}\right){\mathrm{e}}^{i\left(\mathbf{G}-{\mathbf{G}}^{\mathit{\text{'}}}\right)\cdot \mathbf{r}}{\mathrm{d}}^3\mathbf{r}\equiv {\displaystyle \frac{1}{\mathrm{\Omega }}}{\mathfrak{I}}_{\mathbf{G}-{\mathbf{G}}^{\mathit{\text{'}}}}}

where {\mathfrak{I}}_{\mathbf{G}} is given by

(67){\mathfrak{I}}_{\mathbf{G}}=\mathrm{\Omega }{\delta }_{\mathbf{G},\mathrm{0}}-{\displaystyle \sum _{\alpha }{V}_{MT}^{\alpha }{e}^{i\mathbf{G}·{\mathbf{r}}_{\alpha }}ℐ(\mathbf{G}{R}_{MT}^{\alpha })}

with V_{MT}^{\alpha }={\displaystyle \frac{\mathrm{4}\pi }{\mathrm{3}}}{(R_{MT}^{\alpha })}^{\mathrm{3}} being the volume of MT sphere α, and

(68)ℐ(x)\equiv {\displaystyle \frac{\mathrm{3}(\sin x-x\cos x)}{x^{\mathrm{3}}}}

To summarize, the mixed basis set is given by

(69)\{{\chi }_j^{\mathbf{q}}(\mathbf{r})\}\equiv \{{\gamma }_{\alpha NLM}^{\mathbf{q}}(\mathbf{r}),P_i^{\mathbf{q}}(\mathbf{r})\}

The main ingredients of the implementation using the mixed basis is the construction of the bare Coulomb matrix v_ij(q) and the evaluation of the production expansion coefficients M_{nm}^i(\mathbf{k},\mathbf{q}). More details can be found in the footnoted reference (see the footnote on Page 6).

Frequency dependence

How to treat the frequency dependence of the dielectric function has significant influences on the accuracy of the results as well as the efficiency of the implementation. Various techniques or approximations have been developed and roughly they fall into the following categories:

● static approximations

● generalized plasmon models

● imaginary frequency plus analytic continuation approach

● the Hilbert transform approach

● the contour deformation approach

In the following we will discuss the essence of each scheme in turn.

Static approximations

The simplest approximation concerning the treatment of frequency is to neglect the frequency dependence at all. The first static approximation to GW is the so-called static Coulomb-hole and screened-exchange (COHSEX) approximation proposed by Hedin [10]. The real part of the self-energy can be written as

(70)Re{\displaystyle \sum \mathbf{r},{\mathbf{r}}^{\prime };\omega =-{\displaystyle \sum _{n\mathbf{k}}^{occ}{\psi }_{n\mathbf{k}}(\mathbf{r}){\psi }_{n{\mathbf{k}}^{\prime }}^{*}({\mathbf{r}}^{\prime })\Re W({\mathbf{r}}^{\prime },\mathbf{r};\omega -{\mathit{є}}_{\mathit{n}\mathbf{k}})-{\displaystyle \sum _{n\mathbf{k}}{\psi }_{n\mathbf{k}}(\mathbf{r}){\psi }_{n\mathbf{k}}^{*}({\mathbf{r}}^{\prime }){\displaystyle \frac{\mathrm{1}}{\pi }}{\displaystyle {\int }_{\mathrm{0}}^{\infty }d{\omega }^{\prime }{\displaystyle \frac{\Im W_c({\mathbf{r}}^{\prime },\mathbf{r};\omega )}{\omega -{\mathit{є}}_{\mathit{n}\mathbf{k}}-{\omega }^{\prime }}}}}}}

where  represents taking the Cauchy principle value. The static COHSEX approximation can be obtained by setting \omega -{є}_{n\mathbf{k}}=\mathrm{0} in Eq. (70), which gives

(71)$$\begin{gathered} {\displaystyle {\sum }_{}^{\mathrm{C}\mathrm{O}\mathrm{H}\mathrm{S}\mathrm{E}\mathrm{X}}\left(\mathbf{r},{\mathbf{r}}^{\mathit{\text{'}}}\right)}={\displaystyle {\sum }_{}^{\mathrm{C}\mathrm{O}\mathrm{H}}\left(\mathbf{r},{\mathbf{r}}^{\mathit{\text{'}}}\right)+{\displaystyle {\sum }_{}^{\mathrm{S}\mathrm{E}\mathrm{X}}\left(\mathbf{r},{\mathbf{r}}^{\mathit{\text{'}}}\right)}} \\ {\displaystyle {\sum }_{}^{\mathrm{C}\mathrm{O}\mathrm{H}}\left(\mathbf{r},{\mathbf{r}}^{\mathit{\text{'}}}\right)}={\displaystyle \frac{1}{2}}\delta \left(\mathbf{r}-{\mathbf{r}}^{\mathit{\text{'}}}\right)\left[W\left(\mathbf{r},{\mathbf{r}}^{\mathit{\text{'}}};0\right)-v\left(\mathbf{r}-{\mathbf{r}}^{\mathit{\text{'}}}\right)\right] \\ {\displaystyle {\sum }_{}^{\mathrm{S}\mathrm{E}\mathrm{X}}\left(\mathrm{r},{\mathbf{r}}^{\mathit{\text{'}}}\right)=-{\displaystyle \sum _{n\mathbf{k}}^{\mathrm{o}\mathrm{c}\mathrm{c}}{\psi }_{n\mathbf{k}}\left(\mathbf{r}\right){\psi }_{n\mathbf{k}}^{*}\left({\mathbf{r}}^{\mathit{\text{'}}}\right)W\left(\mathbf{r},{\mathbf{r}}^{\mathit{\text{'}}};0\right)}} \end{gathered}$$

The static COHSEX tends to overestimate the band gaps of many semiconductors [46]. On the other hand, since the COHSEX self-energy is hermitian, it can be calculated self-consistently, and the resultant single-particle energies and wave functions can be used as the starting point for a full G₀W₀ calculation. This self-consistent COHSEX based G₀W₀ approach has been used recently for a series of systems with remarkable success [15,47-50].

An approach closely related to the static COHSEX approximation is the model GW approach first proposed by Gygi and Baldereschi [51], in which the self-energy correction with respect to the LDA exchange-correlation potential is approximated by

(72)\delta {\displaystyle \sum (\mathbf{r},{\mathbf{r}}^{\prime })}=-\rho (\mathbf{r},{\mathbf{r}}^{\prime })\delta W(\mathbf{r}-{\mathbf{r}}^{\prime })

where ρ(r, r′) is the density matrix function, and δW(r-r′) is a model function for the screened Coulomb interaction correction, whose Fourier transform takes the form

(73)\delta W(\mathbf{q})={\displaystyle \frac{\mathrm{4}\pi }{\mathrm{\Omega }{\mathbf{q}}^{\mathrm{2}}}}[{ϵ}_{SC}^{-\mathrm{1}}(\mathbf{q},\omega =\mathrm{0})-{ϵ}_M^{-\mathrm{1}}(\mathbf{q},\omega =\mathrm{0})]

{ϵ}_{SC}^{-\mathrm{1}} and {ϵ}_M^{-\mathrm{1}} are the inverse dielectric functions for a semiconductor and a metal, respectively, both approximated by some model functions. This model GW approach was used by Massida and coworkers to investigate electronic band structure properties of transition metal oxides [52-54].

Generalized plasmon models

In the random-phase approximation to the polarization, the screening mainly comes from two types of excitations [11], the plasmon excitation and electron-hole excitations. The former corresponds to the collective excitation of all electrons with respect to fixed positive ions. It turns out that for self-energy, the plamson excitation is dominant, and it is therefore physically plausible to approximate the polarization by retaining only the plasmon excitation, leading to the plasmon pole model. In the following, we will mainly use the planewaves representation for the discussion since various plasmon models originate from the homogeneous electron gas (HEG), for which the reciprocal space is the representation of choice.

For the homogeneous electron system, the plasmon pole model gives the polarization with a single pole [11,46]

(74)\Im {ϵ}^{-\mathrm{1}}(\mathbf{q}\mathbf{,}\omega )=A(a)\delta (\omega -{\omega }_p(\mathbf{q}))

where ω(q) is q-dependent plasmon frequency. For inhomogeneous systems, this plasmon pole model is usually generalized into the following form

(75)\Im {ϵ}_{\mathbf{G}{\mathbf{G}}^{\mathbf{\text{'}}}}^{-\mathrm{1}}(\mathbf{q}\mathbf{,}\omega )=A_{\mathbf{G}{\mathbf{G}}^{\mathbf{\text{'}}}}(\mathbf{q})\delta \left(\omega -{\tilde{\omega }}_{\mathbf{G}{\mathbf{G}}^{\prime }}(\mathbf{q})\right)

hence called generalized plasmon pole (GPP) model. Using the Kramers-Kronig’s relation between the real and imaginary part of the inverse dielectric function [3]

(76)\Re {ϵ}^{-\mathrm{1}}(\mathbf{q},\omega )=\mathrm{1}+{\displaystyle \frac{\mathrm{2}}{\pi }}{\displaystyle \underset{\mathrm{0}}{\overset{\infty }{\int }}d{\omega }^{\prime }{\displaystyle \frac{{\omega }^{\prime }{\Im }^{-\mathrm{1}}(\mathbf{q},\omega )}{{\omega }^{\prime \mathrm{2}}-{\omega }^{\mathrm{2}}}}}

one has

(77)\Re {ϵ}_{\mathbf{G}{\mathbf{G}}^{\prime }}^{-\mathrm{1}}(\mathbf{q},\omega )={\delta }_{\mathbf{G}{\mathbf{G}}^{\prime }}+{\displaystyle \frac{\mathrm{2}}{\pi }}{\displaystyle \frac{{\tilde{\omega }}_{\mathbf{G}{\mathbf{G}}^{\prime }}(\mathbf{q})A_{\mathbf{G}{\mathbf{G}}^{\prime }}(\mathbf{q})}{{\tilde{\omega }}_{\mathbf{G}{\mathbf{G}}^{\prime }}^{\mathrm{2}}(\mathbf{q})-{\omega }^{\mathrm{2}}}}

For given G, G′ and q, there are two unknown parameters, A_{G,G^{\prime }}(q) and {\tilde{\omega }}_{{GG}^{\prime }}(\mathbf{q}) that need to be fixed.

Hybertsen-Louie (HL) model: In the Hybertsen-Louie’s GPP model [12], the two parameters are determined by 1) the static dielectric function (ω = 0); and 2) the f-sum rules [55]

(78){\displaystyle {\int }_0^{\infty }\omega S{ϵ}_{{GG}^{\prime }}^{-1}\left(\mathbf{q},\omega \right)d\omega =-{\displaystyle \frac{\pi }{2}}{\mathrm{\Omega }}_{\mathbf{G}{\mathbf{G}}^{\prime }}^2}\left(\mathbf{q}\right)

Where

(79){\mathrm{\Omega }}_{\mathbf{G}{\mathbf{G}}^{\prime }}^{\mathrm{2}}(\mathbf{q})\equiv {\displaystyle \frac{\mathrm{4}\pi (\mathbf{q}+\mathbf{G})\cdot (\mathbf{q}+{\mathbf{G}}^{\prime })}{\left|\mathbf{q}+\mathbf{G}\right|\left|\mathbf{q}+{\mathbf{G}}^{\prime }\right|}}\rho (\mathbf{G}-{\mathbf{G}}^{\prime }).

The final expressions are

(80)\begin{array}{l}{A}_{\mathbf{G}{\mathbf{G}}^{\prime }}(\mathbf{q})={\displaystyle \frac{\pi }{\mathrm{2}}}{[{\delta }_{\mathbf{G}{\mathbf{G}}^{\prime }}-{ϵ}_{\mathbf{G}{\mathbf{G}}^{\prime }}^{-\mathrm{1}}(\mathbf{q},\mathrm{0})]}^{\mathrm{1}/\mathrm{2}}\left|{\mathrm{\Omega }}_{\mathbf{G}{\mathbf{G}}^{\prime }}\right|\hfill \\ {\tilde{\omega }}_{\mathbf{G}{\mathbf{G}}^{\prime }}(\mathbf{q})={\displaystyle \frac{\left|{\mathrm{\Omega }}_{\mathbf{G}{\mathbf{G}}^{\prime }}\right|}{{[{\delta }_{\mathbf{G}{\mathbf{G}}^{\prime }}-{ϵ}_{\mathbf{G}{\mathbf{G}}^{\prime }}^{-\mathrm{1}}(\mathbf{q},\mathrm{0})]}^{\mathrm{1}/\mathrm{2}}}}\hfill \end{array}

Godby-Needs (GN) model: Godby and Needs [56] proposed a different way to determine the parameters in the GPP model: besides the requirement to give the static inverse dielectric function, the second condition they introduced is that the GPP model reproduce the inverse dielectric function at a chosen imaginary frequency, usually ω= iω_p with {\omega }_p\equiv \sqrt{\mathrm{4}\pi \rho (\mathrm{0})} being the classical Drude plasmon frequency. In this case, one obtains

(81)A_{{GG}^{\prime }}(q)={\displaystyle \frac{\pi }{\mathrm{2}}}{\omega }_p^{\mathrm{1}/\mathrm{2}}{[({\delta }_{\mathbf{G}{\mathbf{G}}^{\prime }}-{ϵ}_{\mathbf{G}{\mathbf{G}}^{\prime }}^{-\mathrm{1}}(\mathbf{q},\mathrm{0}))\times ({ϵ}_{\mathbf{G}{\mathbf{G}}^{\prime }}^{-\mathrm{1}}(\mathbf{q},\mathrm{0})-{ϵ}_{\mathbf{G}{\mathbf{G}}^{\prime }}^{-\mathrm{1}}(\mathbf{q},i{\omega }_p))]}^{\mathrm{1}/\mathrm{2}}

{\tilde{\omega }}_{\mathbf{G}{\mathbf{G}}^{\prime }}(\mathbf{q})={\omega }_p^{\mathrm{1}/\mathrm{2}}{\left[{\displaystyle \frac{{ϵ}_{\mathbf{G}{\mathbf{G}}^{\prime }}^{-\mathrm{1}}(\mathbf{q},\mathrm{0})-{ϵ}_{\mathbf{G}{\mathbf{G}}^{\prime }}^{-\mathrm{1}}(\mathbf{q},i{\omega }_p)}{{\delta }_{\mathbf{G}{\mathbf{G}}^{\prime }}-{ϵ}_{\mathbf{G}{\mathbf{G}}^{\prime }}^{-\mathrm{1}}(\mathbf{q},\mathrm{0})}}\right]}^{\mathrm{1}/\mathrm{2}}

von der Linden-Horsch (vdLH) model: In the HL or GN GPP models, for each q there are totally \mathrm{2}{N}_{\mathrm{P}\mathrm{W}}^{\mathrm{2}} (N_PW: the number of planewaves) parameters to determine, which is not always numerically stable. In some cases, unphysical complex plasmon frequencies can arise for very large G. To overcome this difficulty, vor der Linden and Horsch (vdLH) [57] developed a different GPP model that depends on only 2 N_PW parameters, and is numerically robust. The first step in the construction of the vdLH model is to diagonalize the the symmetrized static dielectric function

(82){ϵ}_{\mathbf{G}{\mathbf{G}}^{\prime }}(\mathbf{q},\mathrm{0})={\displaystyle \sum _i{U}_{\mathbf{G}i}(\mathbf{q})d_i(\mathbf{q})U_{{\mathbf{G}}^{\prime }i}^{*}}(\mathbf{q})

where {U_Gi} and {d_i} are eigenvectors and eigenvalues of ϵ(\mathbf{q},\mathrm{0}), respectively. The inverse dielectric function at arbitrary frequency is obtained by assuming that eigenvectors {U_Gi} do not depend on the frequency, and the frequency dependence occurs only for eigenvalues.

(83){ϵ}_{\mathbf{G}{\mathbf{G}}^{\prime }}^{-\mathrm{1}}(\mathbf{q},\omega )={\displaystyle \sum _i{U}_{\mathbf{G}i}(\mathbf{q})}{d}_i^{-\mathrm{1}}(\mathbf{q},\omega )U_{{\mathbf{G}}^{\prime }i}^{*}(\mathbf{q})

and it is assumed that d_i^{-\mathrm{1}}(\mathbf{q},\omega ) takes the plasmon-pole model like the expression

(84)d_i^{-\mathrm{1}}(\mathbf{q},\omega )=\mathrm{1}+{\displaystyle \frac{z_i(\mathbf{q})}{{\omega }^{\mathrm{2}}-{[{\omega }_i(\mathbf{q})-i\eta ]}^{\mathrm{2}}}}

As in the HL GPP model, the parameters z_i(q) and ω_i(q) are determined again from the static limit

(85)d_i^{-\mathrm{1}}(\mathbf{q})=\mathrm{1}-{\displaystyle \frac{z_i(\mathbf{q})}{{\omega }_i^{\mathrm{2}}(\mathbf{q})}}

and the f-sum rules. After some algebraic formulation, one obtains

(86)\begin{array}{l}{z}_i(\mathbf{q})={\displaystyle \sum _{\mathbf{G}{\mathbf{G}}^{\prime }}{U}_{{\mathbf{G}}^{\prime }i}^{*}(\mathbf{q}){\mathrm{\Omega }}_{\mathbf{G}{\mathbf{G}}^{\prime }}^{\mathrm{2}}{U}_{{\mathbf{G}}^{\prime }i}}(\mathbf{q})\hfill \\ {\omega }_i(\mathbf{q})={\left[{\displaystyle \frac{z_i(\mathbf{q})}{\mathrm{1}-d_i^{-\mathrm{1}}(\mathbf{q})}}\right]}^{\mathrm{1}/\mathrm{2}}\hfill \end{array}

Engel-Farid model: A further extension was made by G. E. Engel and B. Farid (EF) [58]: besides taking the f-sum rule into account, the exact behavior at large frequency limit is also taken into account. Based on the formal analysis, the inverse dielectric function at large frequency limit satisfies the following relation

(87)\underset{\omega \to \infty }{\lim }{\omega }^2{ϵ}^{-1}\left(\mathbf{r},{\mathbf{r}}^{\prime };\omega \right)=-{\displaystyle \frac{2}{\pi }}{\displaystyle {\int }_0^{\infty }\omega S{ϵ}^{-1}\left(\mathbf{r},{\mathbf{r}}^{\prime };\omega \right)d\omega }

In the reciprocal space representation and using the f-sum rules (Eq. (78)), one obtains

(88)\underset{\omega \to \infty }{\lim }{\omega }^{\mathrm{2}}{ϵ}_{\mathbf{G}{\mathbf{G}}^{\prime }}^{-\mathrm{1}}(\mathbf{q},\omega )={\mathrm{\Omega }}_{\mathbf{G},{\mathbf{G}}^{\prime }}^{\mathrm{2}}(\mathbf{q})\equiv L_{\mathbf{G},{\mathbf{G}}^{\prime }}(\mathbf{q})

The EF GPP model for the inverse dielectric function is then defined in the planewave representation as (in the matrix form)

(89){ϵ}^{-\mathrm{1}}(\mathbf{q},\omega )={[{\omega }^{\mathrm{2}}{L}^{-\mathrm{1}}+ϵ(\mathbf{q},\mathrm{0})]}^{-\mathrm{1}}

The imaginary frequency plus analytic continuation approach

A relatively simple approach to consider the full frequency dependence of the dielectric function and the correlation selfenergy is to compute them first on the imaginary frequency axis, ω = iu, and then obtain the self-energy on the real axis by analytic continuation, hence termed as the IF+ AC scheme. The polarization function is a smooth function of the imaginary frequency, so that only a small number of discrete frequency points are needed

(90)P_{ij}(\mathbf{q},iu)=\mathrm{2}{\displaystyle \sum _{\mathbf{k}}^{BZ}{\displaystyle \sum _n^{occ}{\displaystyle \sum _m^{unocc}{\displaystyle \frac{-\mathrm{2}({є}_{n\mathbf{k}}-{є}_{n\mathbf{k}-\mathbf{q}})}{u^{\mathrm{2}}+{({є}_{n\mathbf{k}}-{є}_{n\mathbf{k}-\mathbf{q}})}^{\mathrm{2}}}}}}}\times M_{nm}^i(\mathbf{k},\mathbf{q}){[M_{nm}^j(\mathbf{k},\mathbf{q})]}^{*}

(90)Making use of the inversion symmetry of W^c on the imaginary frequency axis, W_{ij}^{\mathrm{c}}(\mathbf{q},\mathrm{i}u)=W_{ij}^{\mathrm{c}}(\mathbf{q},-\mathrm{i}u), the correlation self-energy for the imaginary frequency reads

(91){\displaystyle {\sum }_{n\mathbf{k}}^{\mathrm{c}}(\mathrm{i}u)}={\displaystyle \sum _{\mathbf{q}}^{BZ}{\displaystyle \sum _m{\displaystyle {\int }_{\mathrm{0}}^{\infty }{\displaystyle \frac{\mathrm{d}{u}^{\mathit{\text{'}}}}{\mathrm{2}\mathrm{\pi }}}}{\displaystyle \frac{\mathrm{2}({\mathit{є}}_{\mathit{m}\mathbf{k}-\mathbf{q}}-\mathrm{i}u)X_{nm}(\mathbf{k},\mathbf{q};\mathrm{i}{u}^{\mathit{\text{'}}})}{u^{\mathit{\text{'}}\mathrm{2}}+{({є}_{m\mathbf{k}-\mathbf{q}}-\mathrm{i}u)}^{\mathrm{2}}}}}}.

The integrand in Eq. (91) is peaked around u′ = u for small {є}_{m\mathbf{k}-\mathbf{q}}\mathit{} such that a direct numerical integration is unstable. This can be avoided by adding and subtracting the following analytically integrable term [40]

(92){\displaystyle \underset{\mathrm{0}}{\overset{\infty }{\int }}{\displaystyle \frac{d{u}^{\prime }}{\mathrm{2}\pi }}}{X}_{nm}(\mathbf{k},\mathbf{q};iu){\displaystyle \frac{({є}_{m\mathbf{k}-\mathbf{q}}-iu)}{u^{\prime \mathrm{2}}+{({є}_{m\mathbf{k}-\mathbf{q}}-iu)}^{\mathrm{2}}}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathrm{sgn}({є}_{m\mathbf{k}-\mathbf{q}})X_{nm}(\mathbf{k},\mathbf{q};iu)

which gives

(93){\displaystyle {\sum }_{n\mathbf{k}}^c(iu)={\displaystyle \sum _{\mathbf{q}}^{BZ}{\displaystyle \sum _m\left\{{\displaystyle {\int }_{\mathrm{0}}^{\infty }{\displaystyle \frac{d{u}^{\prime }}{\mathrm{2}\pi }}[X_{nm}(\mathbf{k},\mathbf{q};i{u}^{\prime })-X_{nm}(\mathbf{k},\mathbf{k};iu)]}\times {\displaystyle \frac{\mathrm{2}({є}_{m\mathbf{k}-\mathbf{q}}-iu)}{u^{\prime \mathrm{2}}+{({є}_{m\mathbf{k}-\mathbf{q}}-iu)}^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathrm{sgn}({є}_{m\mathbf{k}-\mathbf{q}})X_{nm}(\mathbf{k},\mathbf{q},iu)\right\}}}}

The integrand is now a smooth function of u for any {є}_{m\mathbf{k}-\mathbf{q}}, and therefore a standard Gaussian quadrature can be used. In practice, a double Gauss-Legendre quadrature [13,59] is often used, in which the semi-infinite integral is divided into two intervals, [\mathrm{0},{\omega }_{\mathrm{0}}] and [{\omega }_{\mathrm{0}},\infty ), and the integration in each interval is carried out by standard Gauss-Legendre quadrature. In most cases, the integration converges quite quickly with respect to the number of discrete frequency points when an appropriate ω₀ is chosen [59] (see footnote on Page 6).

Once the correlation self-energy matrix elements along the imaginary axis are calculated according to Eq. (93), they can be fitted by a function with N_p poles [31]

(94){\displaystyle {\sum }_{n\mathbf{k}}^c(iu)}={\displaystyle \sum _p^{N_p}{\displaystyle \frac{a_{p;n\mathbf{k}}}{iu-b_{p;n\mathbf{k}}}}}

where a_p;nk and b_p;nk are fitting parameters. Eq. (94) is then analytically continued onto the real frequency axis.

The IF+ AC approach to treat the frequency dependence is relatively simple to implement and quite efficient in terms of the computational demand, but it has also such a disadvantage that the error of the analytic continuation cannot be controlled: Increasing the number of imaginary frequencies does not guarantee a convergence to more accurate real-frequency self-energies.

The Hilbert transform approach

To avoid the uncertainty related to the use of the IF+ AC approach, one can treat everything along the real frequency directly. In this case, it is often advantageous to use the spectral representation of G, P, W and Σ^c [11,60].

Polarization matrix: In the general matrix form, the spectral representation of the polarization function reads

(95)$$\begin{gathered} S_{ij}\left(\mathbf{q},\omega \right)\equiv -{\displaystyle \frac{1}{\pi }}S{P}_{ij}\left(\mathbf{q},\omega \right)\mathrm{sgn}\left(\omega \right) \\ =2{\displaystyle \sum _{\mathbf{k}}^{BZ}{{\displaystyle \sum _n^{occ}{\displaystyle \sum _m^{unocc}{M}_{nm}^i\left(\mathbf{k},\mathbf{q}\right)}\left[M_{nm}^j\left(\mathbf{k},\mathbf{q}\right)\right]}}^{*}}\times \delta \left(\left|\omega \right|-{\omega }_{n\mathbf{k},m\mathbf{k}-\mathbf{q}}\right)\mathrm{sgn}\left(\omega \right) \end{gathered}$$

Since only a small number of occupied-unoccupied pairs contribute for each frequency ω, S_ij(q, ω) can be calculated quite efficiently. The full polarization function can then be calculated by the following Hilbert transform

(96)P_{ij}(\mathbf{q},\omega )={\displaystyle {\int }_{-\infty }^{\infty }d{\omega }^{\prime }{\displaystyle \frac{S_{ij}(\mathbf{q},{\omega }^{\prime })}{\omega -{\omega }^{\prime }-i\eta \mathrm{sgn}({\omega }^{\prime })}}}

Correlation self-energies　Correlation self-energies on real-frequency can be calculated directly by using Eq. (44) [60]: one first calculates the auxiliary quantities X_nm(k, ω) on a series of real frequencies, and the correlation selfenergies are then calculated by a numerical integration over real frequency in terms of Eq. (44). This scheme is useful when the full frequency dependent self-energy is required. On the other hand, in G₀W₀ calculation, only the self-energy and its first-order derivative at the corresponding Kohn-Sham orbital energy are needed (See Eq. (33)). Following Ref. [60, first define the Hilbert-transformed screened Coulomb interaction

(97){\tilde{W}}_{ij}^{\pm }(\mathbf{q},\omega )\equiv {\displaystyle \frac{i}{\mathrm{2}\pi }}{\displaystyle {\int }_{-\infty }^{\infty }d{\omega }^{\prime }{\displaystyle \frac{W_{ij}^c(\mathbf{q},{\omega }^{\prime })}{\omega -{\omega }^{\prime }\pm i\eta }}}

from which one obtains

(98){\tilde{X}}_{nm}^{\pm }(\mathbf{k},\mathbf{q};\omega )\equiv {\displaystyle \sum _{ij}[M_{nm}^i{(\mathbf{k},\mathbf{q})}^{*}]{\tilde{W}}_{ij}^{\pm }(\mathbf{q},\omega )M_{nm}^j(\mathbf{k},\mathbf{q})}

Using the fact that W^c is an even function of ω, it is easy to prove that

(99){\tilde{X}}_{nm}^{\pm }(\mathbf{k},\mathbf{q};-\omega )=-{\tilde{X}}_{nm}^{\mp }(\mathbf{k},\mathbf{q};\omega )

We note that it is actually possible to calculate {\tilde{X}}_{nm}^{\pm }(\mathbf{k},\mathbf{q};\omega ) directly

(100){\tilde{X}}_{nm}^{\pm }(\mathbf{k},\mathbf{q};\omega )={\displaystyle \frac{i}{\mathrm{2}\pi }}{\displaystyle {\int }_{-\infty }^{\infty }d{\omega }^{\prime }{\displaystyle \frac{X_{nm}(\mathbf{k},\mathbf{q};{\omega }^{\prime })}{\omega -{\omega }^{\prime }\pm i\eta }}}

which is likely more efficient since the size of X is usually smaller than W_c so that the number of numerical Hilbert transforms can be reduced. Using {\tilde{X}}_{nm}^{\pm }, the diagonal correlation self-energy elements at an arbitrary frequency can be calculated as

(101){\displaystyle {\sum }_{n\mathbf{k}}^c(\omega )={\displaystyle \sum _{\mathbf{q}}{\displaystyle \sum _m{\tilde{X}}_{nm}^{s_m}(\mathbf{k},\mathbf{q};\omega -{є}_{m\mathbf{k}-\mathbf{q}})}}}

where we use the abbreviation s_m\equiv \mathrm{sgn}({є}_{m\mathbf{k}-\mathbf{q}}-{є}_F). In particular, at \omega ={є}_{n\mathbf{k}} (the abbreviation {\mathrm{\Delta }}_{nm}\equiv {є}_{n\mathbf{k}}-{є}_{m\mathbf{k}-\mathbf{q}} is used to simplify the notation)

(102){\displaystyle {\sum }_n^c({є}_{n\mathbf{k}})={\displaystyle \sum _{\mathbf{q}}{\displaystyle \sum _m\mathrm{sgn}({\mathrm{\Delta }}_{nm}){\tilde{X}}_{nm}^{s_m\mathrm{sgn}({\mathrm{\Delta }}_{nm})}(\mathbf{k},\mathbf{q};\left|{\mathrm{\Delta }}_{nm}\right|)}}}

for which Eq. (99) has been used. Since {\tilde{X}}_{nm}^{\pm } are usually calculated only on a set of discrete points {ω_i}, (i = 1, 2, · · ·, N_w), to obtain {\tilde{X}}_{nm}^{\pm } at \omega =\left|{\mathrm{\Delta }}_{nm}\right|, a linear interpolation can be used; assuming , one has

(103){\tilde{X}}^{\pm }({\mathrm{\Delta }}_{nm})={\displaystyle \frac{{\omega }_{i+\mathrm{1}}-\left|{\mathrm{\Delta }}_{nm}\right|}{{\omega }_{i+\mathrm{1}}-{\omega }_i}}{\tilde{X}}^{\pm }({\omega }_i)+{\displaystyle \frac{{\omega }_i-\left|{\mathrm{\Delta }}_{nm}\right|}{{\omega }_i-{\omega }_{i+\mathrm{1}}}}{\tilde{X}}^{\pm }({\omega }_{i+\mathrm{1}})

The derivative of the correlation self-energy at \omega ={є}_{n\mathbf{k}} can be calculated as

(104){\displaystyle \frac{\partial \Re {\displaystyle {\sum }_{n\mathbf{k}}^c({є}_{n\mathbf{k}})}}{\partial {є}_{n\mathbf{k}}}}={\displaystyle \sum _{\mathbf{q}}{\displaystyle \sum _m{\displaystyle \frac{{\tilde{X}}_{nm}^{sm}(\mathbf{k},\mathbf{q};{\omega }_{i+\mathrm{1}})-{\tilde{X}}_{nm}^{sm}(\mathbf{k},\mathbf{q};{\omega }_i)}{{\omega }_{i+\mathrm{1}}-{\omega }_i}}}}

The contour deformation approach

Finally we discuss the contour deformation (CD) approach [11,13] to treat the frequency integration for the calculation of the correlation self-energy. The greatest advantage of the CD approach is that the integration is performed along the imaginary frequency, but the self-energy is calculated directly for real frequencies so that the uncertainty related to the analytic continuation can be avoided.

We start from the matrix form of the correlation self-energy, and to simplify the notation, we drop the dependence on wavevectors (k and q),

(105){\displaystyle {\sum }_n^c(\omega )={\displaystyle \sum _m{\displaystyle \frac{i}{\mathrm{2}\pi }}}{\displaystyle {\int }_{-\infty }^{\infty }d{\omega }^{\prime }{\displaystyle \frac{X_{nm}({\omega }^{\prime })}{\omega +{\omega }^{\prime }-{є}_m-i\eta \mathrm{sgn}({є}_F-{є}_m)}}}}.

As illustrated in Fig. 2, the integrand in the equation above has two types of poles: 1) the pole from the Green’s function, {\omega }^{\prime }={є}_m-\omega +i\eta \mathrm{sgn}({є}_F-{є}_m), which is above (below) the real frequency for occupied (unoccupied) states; 2) the poles from the screened Coulomb interaction, which are above (below) the real axis for ω′<0 (ω′>0). Using the fact W^c(\omega )\propto {\displaystyle \frac{\mathrm{1}}{{\left|\omega \right|}^{\mathrm{2}}}} as \left|\omega \right|\to \infty, it is possible to transform the integration over the real frequency to that over the imaginary frequency (ω = iu) plus the contributions from the poles of the Green’s function by using the contour deformation technique and the residue theorem [25],

{\displaystyle {\sum }_n^c(\omega )={\displaystyle \sum _m\left\{{\displaystyle {\int }_{\mathrm{0}}^{\infty }{\displaystyle \frac{du}{\mathrm{2}\pi }}{X}_{nm}(iu){\displaystyle \frac{\mathrm{2}({є}_m-\omega )}{{({є}_m-\omega )}^{\mathrm{2}}+u^{\mathrm{2}}}}+X_{nm}({є}_m-\omega )[\theta ({є}_m-{є}_F)\theta (\omega -{є}_m)-\theta ({є}_F-{є}_m)\theta ({є}_m-\omega )]}\right\}}}

(106)

In the CD approach, the polarization matrix for both real and imaginary frequencies is needed for the evaluation of the correlation self-energies, so apparently it is more expensive than the IF+ AC approach and the Hilbert transform approach. But it is likely the most accurate approach due to the fact that 1) the numerical integration along the imaginary frequency converges very quickly with respect to the number of mesh points (N_w ~ 10 - 20), in contrast to that along the real frequency for which hundreds of mesh points may be necessary [60]; 2) the correlation self-energy on the real frequency are directly calculated so that it does not rely on the accuracy of the analytic continuation whose error is often difficult to control.

The choice of the reference H₀

In the spirit of the “best G best W” strategy, the reference Hamiltonian H₀ should be chosen to deliver the best possible approximation for G and W in the G₀W₀ framework. A certain dependence of the G₀W₀ results on the reference H₀, although theoretically unsatisfactory, is there unavoidable. When one talks about G₀W₀, it is important to indicate explicitly the starting point used. For extended sp-electron systems, the LDA/GGA KS single particle Hamiltonian is the most popular choice in practice. We note, however, that using the KS H₀ as the reference for GW calculations is mainly based on pragmatic considerations without much formal justifications since the KS single-particle energies in general cannot been identified as the quasi-particle energies, even at an approximate level. Besides the easy availability of LDA/GGA, there are at least two factors that contribute to the popularity of using the LDA/GGA H₀: 1) LDA/GGA can often describe the dispersion (with respect to k-vector) of quasi-particle band structures very well, although the band gaps are significantly underestimated; and 2) The macroscopic dielectric constants of many semiconductors calculated from LDA/GGA eigenenergies and wave functions, which characterizes the strength of screening, are often in good agreement with experiment.

Formally the HF single-particle Hamiltonian is more appropriate as the GW starting point since the HF energies can be related to quasi-particle excitations in terms of Koopmans’s theorem [61], although in a quite crude manner. However, the HF descriptions of metals and semiconductors are usually much less satisfactory than the LDA/GGA ones due to its failure to treat the screening (dynamic correlation) effect properly. For wide-gap insulators or finite molecular systems, in which the screening is usually very weak, the HF H₀ might become more appropriate. However, since the implementation of the HF method for extended systems is much more demanding, most G₀W₀ calculations use LDA/GGA H₀ as the starting point.

Besides LDA/GGA Kohn-Sham single-particle Hamiltonian [12,13], other alternative references are also used to improve the starting point for G₀W₀, including H₀ from the optimized effective potential for exact exchange (OEPx) [62,63], the LDA/GGA plus the Hubbard U correction (LDA/GGA+ U) [64-66], and the hybrid functionals approach [67]. The latter two are especially useful for systems with partially occupied d- or f-shells, often termed as strongly correlated systems, for which the LDA/GGA descriptions are often qualitatively wrong so that they are problematic as the reference for G₀W₀. We refer to Refs. [24,66] for more detailed discussions on this issue.

Concluding remarks

In this paper we present an overview of the theoretical foundation of the GW method and main numerical techniques developed for its implementation. The basic theoretical framework for Green’s function based approaches has been set up for several decades, but they remain highly specialized and limited to relatively simple systems. The situations are changing quickly in recent years thanks to latest methodological developments as well as the rapid increase of computer power. Even more importantly, there is increasingly stronger demand in fundamental and applied research for accurate theoretical descriptions of electronic band structure. For example, efficient solar energy conversion via photo-catalytic or photo-voltaic processes is currently one of the most intensively pursued scientific frontiers, and electronic band structure is one of the most important parameters of working materials used in solar energy conversion. The Green’s function based first-principles methods are therefore highly promising to provide valuable information to interpret existing experimental findings and to achieve rational design of new materials.

There are still several great challenges to meet before it becomes feasible to employ the GF methods routinely for more complicated systems. The first challenge is related to the efficiency of the current GF methods. In spite of the great efforts that have been invested for the implementation of the GW method, including various numerical techniques reviewed in this paper, it is still computationally formidable to routinely treat systems with hundreds of atoms. We expect that some paradigm shift is needed for the implementation of the GF methods that exploits both latest algorithmic advances and new computer architectures. Several interesting developments along this direction have appeared recently [17,19,20,68], but further efforts are obviously needed. We note that in this aspect the experiences from recent developments in linear scaling post-HF quantum chemistry methods [69,70] may give useful inspiration. The second challenge is related to the accuracy of the current GF methods. Even at the GW level, additional approximations are involved in practice. With recent proposal of several approximate self-consistent GW schemes [14-16], it is still under debate what is the optimal way to perform GW calculations. The GW is just the lowest order term in the many-body perturbation expansion of the self-energy with respect to W. For many systems, in particular, open-shell d/f-electron systems, higher order contributions beyond GW are probably necessary. In addition, the current GF methods consider only electronic degrees of freedom, and the screening included in W has only the contribution from electronelectron interaction, but for some materials, especially those used in solar energy conversion, the electron-phonon coupling can play an important role, and the polaronic contribution to screening may become significant [71-73]. How to incorporate such new important physics in the GF methods is still a frontier not well explored.

References

[1] Hüfner, S., Photoelectron Spectroscopy: Prinples and Applications 3rd ed. Berlin: Springer, 2003

[2] Onida, G.; Rubio, A., Rev. Mod. Phys. 2002, 74, 601-659

[3] Yu, P. Y.; Cardona, M., Fundamentals of semiconductors: physics and materials properties 3rd ed. Berlin: Springer, 2001

[4] Parr, R. G.; Yang, W., Density-Functional Theory of Atoms and MoleculesNew York: Oxford University Press, 1989

[5] Dreizler, R. M.; Gross, E. K. U., Density Functional Theory: An Approach to the Quantum Many-Body ProblemBerlin: Springer-Verlag, 1990

[6] Martin, R. M., Electronic Structure: Basic Theory and Practical Methods, Cambridge UK: Cambridge University Press, 2004

[7] Aryasetiawan, F., in Anisimov, V. I., ed., Strong Coulomb Correlations in Electronic Structure Calculations: Beyond the Local Density Approximation Gordon and Breach Science Publishers, 2000 (1)

[8] Fetter, A. L.; Walecka, J. D., Quantum theory of many-particle systems McGraw-Hill, New York, 1971

[9] Hedin, L.; Lundqvist, B. I., Solid State Phys. 1969, 23, 1-181

[10] Hedin, L., Phys. Rev.1965, 139, A796-A823

[11] Aryasetiawan, F.; Gunnarsson, O., Rep. Prog. Phys. 1998, 61, 237-312

[12] Hybertsen, M. S.; Louie, S. G., Phys. Rev. B1986, 34, 5390-5413

[13] Godby, R. W.; Schlüter, M.; Sham, L. J., Phys. Rev. B1988, 37, 10159-10175

[14] Faleev, S. V.; van Schilfgaarde, M.; Kotani, T., Phys. Rev. Lett.2004, 93, 126406

[15] Bruneval, F.; Vast, N.; Reining, L., Phys. Rev. B2006, 74, 045102

[16] Shishkin, M.; Marsman, M.; Kresse, G., Phys. Rev. Lett.2007, 99, 246403

[17] Bruneval, F.; Gonze, X., Phys. Rev. B2008, 78, 085125

[18] Hamann, D. R.; Vanderbilt, D., Phys. Rev. B2009, 79, 045109

[19] Berger, J. A.; Reining, L.; Sottile, F., Phys. Rev. B2010, 82, 2010

[20] Umari, P.; Stenuit, G.; Baroni, S., Phys. Rev. B2010, 81, 115104

[21] Samsonidze, G.; Jain, M.; Deslippe, J.; Cohen, M. L.; Louie, S. G., Phys. Rev. Lett.2011, 107, 186404

[22] Jiang, H.; Gomez-Abal, R.; Rinke, P.; Scheffler, M., Phys. Rev. B2010, 81, 085119

[23] Inkson, J. C., Many-body theory of solids: An IntroductionNew York: Plenum, 1983

[24] Jiang, H.Acta., Acta.Phys. Chim. Sin2010, 26, 1017

[25] Arfken, G. B.; Weber, H. J., Mathematical Methods for Physicists ed. 5th ed. Academic Press, 2001

[26] Linderberg, J., Öhrn Propagators in Quantum Chemistry 2nd ed. John Wiley & Sons, 2004

[27] Zakrzewski, V. G.; Dolgounitcheva, O.; Zakjevskii, A. V.; Ortiz, J. V., Ann. Rep. Comput. Chem2010, 6, 79-94

[28] Shishkin, M.; Kresse, G., Phys. Rev. B2007, 75, 235102

[29] Baldereschi, A.; Tosatti, E., Solid State Commun. 1979, 29, 131-135

[30] Godby, R. W.; Schlüter, M.; Sham, L. J., Phys. Rev. B1987, 36, 6497-6500

[31] Rojas, H. N.; Godby, R. W.; Needs, R. J., Phys. Rev. Lett.1995, 74, 1827-1830

[32] Rieger, M. M.; Steinbeck, L.; White, I. D.; Rojas, H. N.; Godby, R. W., Comput. Phys. Commun.1999, 117, 211-228

[33] Rohlfing, M.; Krüger, P.; Pollmann, J., Phys. Rev. Lett.1995, 75, 3489-3492

[34] Blase, X.; Attaccalite, C.; Olevano, V.prb, 2011, 83: 115103

[35] Helgaker, T.; Jorgensen, P.; Olsen, J., Molecular Electronic-Structure Theory John Wiley & Sons, 2000

[36] Foerster, D.; Koval, P.; Sanchez-Portal, D. J., Chem. Phys.2011, 135, 074105

[37] Gómez-Abal, R.; Li, X.; Scheffler, M.; Ambrosch-Draxl, C., Phys. Rev. Lett.2008, 101, 106404

[38] Li, X., All-Electron G0W0 code based on FP-(L)APW+lo and applications Ph.D. thesis Free University of Berlin, 2008

[39] Li, G. L.; Yin, Z., Phys. Chem. Phys. Chem2011, 13, 2824

[40] Aryasetiawan, F., Phys. Rev. B1992, 46, 13051-13064

[41] Kotani, T.; van Schilfgaarde, M., Solid State Commun.2002, 121, 461-465

[42] Friedrich, C.; Schindlmayr, A.; Blügel, S.; Kotani, T., Phys. Rev. B2006, 74, 045104

[43] Friedrich, C.; Blügel, S.; Schindlmayr, A.prb, 2010, 81: 125102

[44] Aryasetiawan, F.; Gunnarsson, O., Phys. Rev. B1994, 49, 16214-16222

[45] Andersen, O. K., Phys. Rev. B1975, 12, 3060-3083

[46] Aulbur, W. G.; Jönsson, L.; Wilkins, J. W., Solid State Phys.2000, 54, 1-218

[47] Gatti, M.; Bruneval, F.; Olevano, V.; Reining, L., Phys. Rev. Lett.2007, 99, 266402

[48] Vidal, J.; Botti, S.; Olsson, P.; Guillemoles, J.-F.; Reining, L. prl, 2010, 104: 056401

[49] Vidal, J.; Trani, F.; Bruneval, F.; Marques, M. A. L.; Botti, S.prl, 2010, 104: 136401

[50] Botti, S.; Kammerlander, D.; Marques, M. A. L.apl, 2011, 98: 241915

[51] Gygi, F.; Baldereschi, A., Phys. Rev. Lett.1989, 62, 2160-2163

[52] Massidda, S.; Continenza, A.; Posternak, M.; Baldereschi, A., Phys. Rev. Lett.1995, 74, 2323-2326

[53] Massidda, S.; Continenza, A.; Posternak, M.; Baldereschi, A., Phys. Rev. B1997, 55, 13494-13502

[54] Continenza, A.; Massidda, S.; Posternak, M., Phys. Rev. B1999, 60, 15699-15704

[55] Johnson, D. J., Phys. Rev. B1974, 9, 4475-4484

[56] Godby, R. W.; Needs, R. J.prl, 1989, 62: 1169

[57] von der Linden, W.; Horsch, P., Phys. Rev. B1988, 37, 8351-8362

[58] Engel, G. E.; Farid, B., Phys. Rev. B1993, 47, 15931-15934

[59] Jiang, H.; Engel, E. J., Chem. Phys.2007, 127, 184108

[60] Shishkin, M.; Kresse, G., Phys. Rev. B2006, 74, 035101

[61] Szabo, A.; Ostlund, N. S., Modern Quantum ChemistryNew York: McGraw-Hill, 1989

[62] Rinke, P.; Qteish, A.; Neugebauer, J.; Freysoldt, C.; Scheffler, M., N. J. Phys.2005, 7, 126

[63] Rinke, P.; Qteish, A.; Neugebauer, J.; Scheffler, M.phys. stat. sol. (b), 2008, 245: 929

[64] Miyake, T.; Zhang, P.; Cohen, M. L.; Louie, S. G., Phys. Rev. B2006, 74, 245213

[65] Jiang, H.; Gomez-Abal, R. I.; Rinke, P.; Scheffler, M., Phys. Rev. Lett.2009, 102, 126403

[66] Jiang, H.; Gomez-Abal, R. I.; Rinke, P.; Scheffler, M., Phys. Rev. B2010, 82, 045108

[67] Rödl, C.; Fuchs, F.; Furthmüller, J.; Bechstedt, F., Phys. Rev. B2008, 77, 184408

[68] Caramella, L.; Onida, G.; Finocchi, F.; Reining, L.; Sottile, F., Phys. Rev. B2007, 75, 205405

[69] Schütz, M.; Hetzer, G.; Werner, H. J., J. Chem. Phys.1999, 111, 5691

[70] Ayala, P. Y.; Scuseria, G. E. J., Chem. Phys.1999, 110, 3660

[71] Chiodo, L.; Garcia-Lastra, J. M.; Iacomino, A.; Ossicini, S.; Zhao, J.; Petek, H.; Rubio, A.prb, 2010, 82: 045207

[72] Kang, W.; Hybertsen, M. S., Phys. Rev. B2010, 82, 085203

[73] Wang, H.; Wu, F.; Jiang, H. J., PhysChemComm2011, 115, 16180

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg