GW approximation of Hedin's equations

Green's functions

The GW method can be understood in terms of the following eigenvalue equation

$(T+V_{ext}+V_h){\varphi }_{n\mathbf{𝐤}}(\mathbf{𝐫})+\int d\mathbf{𝐫}\mathrm{\Sigma }(\mathbf{𝐫},{\mathbf{𝐫}}^{\prime },\omega =E_{n\mathbf{𝐤}}){\varphi }_{n\mathbf{𝐤}}({\mathbf{𝐫}}^{\prime })=E_{n\mathbf{𝐤}}{\varphi }_{n\mathbf{𝐤}}(\mathbf{𝐫})$

Here $T$ is the kinetic energy, $V_{ext}$ the external potential of the nuclei, $V_h$ the Hartree potential and $E_{n\mathbf{𝐤}}$ the quasiparticle energies with orbitals ${\varphi }_{n\mathbf{𝐤}}$. In contrast to DFT, the exchange-correlation potential is replaced by the many-body self-energy $\mathrm{\Sigma }$ and should be obtained together with the Green's function $G$, the irreducible polarizability $\chi$, the screened Coulomb interaction $W$ and the irreducible vertex function $\mathrm{\Gamma }$ in a self-consistent procedure. For completeness, these equations are

$G(1,2)=G_0(1,2)+\int d(3,4)G_0(1,3)\mathrm{\Sigma }(3,4)G(4,2)$

$\chi (1,2)=\int d(3,4)G(1,3)G(4,1)\mathrm{\Gamma }(3,4;2)$

$W(1,2)=V(1,2)+\int d(3,4)V(1,3)\chi (3,4)W(4,2)$

$\mathrm{\Sigma }(1,2)=\int d(3,4)G(1,3)\mathrm{\Gamma }(3,2;4)W(4,1)$

$\mathrm{\Gamma }(1,2;3)=\delta (1,2)\delta (1,3)+\int d(4,5,6,7)\frac{\delta \mathrm{\Sigma }(1,2)}{\delta G(4,5)}G(4,6)G(7,5)\mathrm{\Gamma }(6,7;3)$

Here the common notation $1=({\mathbf{𝐫}}_1,t_1)$ was adopted and $V$ denotes the bare Coulomb interaction. Note, that these equations are exact and provide an alternative to the Schrödinger equation for the many-body problem. Nevertheless, approximations are necessary for realistic systems. The most popular one is the GW approximation and is obtained by neglecting the equation for the vertex function and using the bare vertex instead:

$\mathrm{\Gamma }(1,2;3)=\delta (1,2)\delta (1,3)$

This means that the equations for the polarizability and self-energy reduce to

$\chi (1,2)=G(1,2)G(2,1)$

$\mathrm{\Sigma }(1,2)=G(1,2)W(2,1)$

while the equations for the Green's function and the screened potential remain the same.

However, in practice, these equations are usually solved in reciprocal space in the frequency domain

$W_{\mathbf{𝐆}{\mathbf{𝐆}}^{\prime }}(\mathbf{𝐪},\omega )={[{\delta }_{\mathbf{𝐆}{\mathbf{𝐆}}^{\prime }}-{\chi }_{\mathbf{𝐆}{\mathbf{𝐆}}^{\prime }}(\mathbf{𝐪},\omega )V_{\mathbf{𝐆}{\mathbf{𝐆}}^{\prime }}(\mathbf{𝐪})]}^{-1}{V}_{\mathbf{𝐆}{\mathbf{𝐆}}^{\prime }}(\mathbf{𝐪})$

$G_{\mathbf{𝐆}{\mathbf{𝐆}}^{\prime }}(\mathbf{𝐪},\omega )={[{\delta }_{\mathbf{𝐆}{\mathbf{𝐆}}^{\prime }}-{\mathrm{\Sigma }}_{\mathbf{𝐆}{\mathbf{𝐆}}^{\prime }}(\mathbf{𝐪},\omega )G_{\mathbf{𝐆}{\mathbf{𝐆}}^{\prime }}^{(0)}(\mathbf{𝐪})]}^{-1}{G}_{\mathbf{𝐆}{\mathbf{𝐆}}^{\prime }}^{(0)}(\mathbf{𝐪})$

In principle Hedin's equations have to be solved self-consistently, where in the first iteration $G^{(0)}$ is the non-interacting Green's function

$G^{(0)}(\mathbf{𝐫},{\mathbf{𝐫}}^{\prime },\omega )=\sum _{n\mathbf{𝐤}}\frac{{\varphi }_{n\mathbf{𝐤}}^{*(0)}(\mathbf{𝐫}){\varphi }_{n\mathbf{𝐤}}^{(0)}({\mathbf{𝐫}}^{\prime })}{\omega -E_{n\mathbf{𝐤}}^{(0)}}$

with ${\varphi }_{n\mathbf{𝐤}}^{(0)}$ being a set of one-electron orbitals and $E_{n\mathbf{𝐤}}^{(0)}$ the corresponding energies. Afterwards the polarizability ${\chi }^{(0)}$ is determined, followed by the screened potential $W^{(0)}$ and the self-energy ${\mathrm{\Sigma }}^{(0)}$. This means that GW calculations require a first guess for the one-electron eigensystem, which is usually taken from a preceding DFT step.

In principle, one has to repeat all steps by the updating the Green's function with the Dyson equation given above in each iteration cycle until self-consistency is reached. In practice, this is hardly ever done due to computational complexity on the one hand (in fact fully self-consistent GW calculations are available as of VASP 6 only).

On the other hand, one observes that by keeping the screened potential $W$ in the first iteration to the DFT level one benefits from error cancelling, which is the reason why often the screening is kept on the DFT level and one aims at self-consistency in Green's function only.

Following possible approaches are applied in practice and selectable within VASP with the ALGO tag.

Single Shot: G₀W₀

Performing only one GW iteration step is commonly referred to the G₀W₀ method. Here the self-energy ${\mathrm{\Sigma }}^{(0)}$ is determined and the corresponding eigenvalue equation is solved. Formally, this is a five step precedure

• Determine the independent particle polarizability ${\chi }_{\mathbf{𝐆}{\mathbf{𝐆}}^{\prime }}^{(0)}(\mathbf{𝐪},\omega )$
• Determine the screened Coulomb potential $W_{\mathbf{𝐆}{\mathbf{𝐆}}^{\prime }}^{(0)}(\mathbf{𝐪},\omega )$
• Determine the self-energy ${\mathrm{\Sigma }}^{(0)}(\mathbf{𝐫},{\mathbf{𝐫}}^{\prime },\omega )$
• Solve the eigenvalue equation $(T+V_{ext}+V_h){\varphi }_{n\mathbf{𝐤}}(\mathbf{𝐫})+\int d\mathbf{𝐫}{\mathrm{\Sigma }}^{(0)}(\mathbf{𝐫},{\mathbf{𝐫}}^{\prime },\omega =E_{n\mathbf{𝐤}}^{(1)}){\varphi }_{n\mathbf{𝐤}}({\mathbf{𝐫}}^{\prime })=E_{n\mathbf{𝐤}}^{(1)}{\varphi }_{n\mathbf{𝐤}}(\mathbf{𝐫})$ for the quasi-particle energies $E_{n\mathbf{𝐤}}^{(1)}$.

To save further computation time, the self-energy is linearized with a series expansion around the Kohn-Sham eigenvalues ${ϵ}_{n\mathbf{𝐤}}$

${\mathrm{\Sigma }}^{(0)}(\mathbf{𝐫},{\mathbf{𝐫}}^{\prime },\omega )\approx {\mathrm{\Sigma }}^{(0)}(\mathbf{𝐫},{\mathbf{𝐫}}^{\prime },{ϵ}_{n\mathbf{𝐤}})+{\frac{\partial {\mathrm{\Sigma }}^{(0)}}{\partial \omega }(\mathbf{𝐫},{\mathbf{𝐫}}^{\prime },\omega )|}_{\omega ={ϵ}_{n\mathbf{𝐤}}}(\omega -{ϵ}_{n\mathbf{𝐤}})$

and the renormalization factor $Z_{n\mathbf{𝐤}}^{(0)}={[1-\mathrm{R}\mathrm{e}\left({\frac{\partial {\mathrm{\Sigma }}^{(0)}}{\partial \omega }(\mathbf{𝐫},{\mathbf{𝐫}}^{\prime },\omega )|}_{\omega ={ϵ}_{n\mathbf{𝐤}}}\right)]}^{-1}$ is introduced. This allows to obtain the G₀W₀ quasi-particle energies from following equation

$E_{n\mathbf{𝐤}}^{(1)}={ϵ}_{n\mathbf{𝐤}}+Z_{n\mathbf{𝐤}}^{(0)}\mathrm{R}\mathrm{e}[⟨{\varphi }_{n\mathbf{𝐤}}|-\frac{\mathrm{\Delta }}{2}+V_{ext}+V_h+{\mathrm{\Sigma }}^{(0)}(\omega ={ϵ}_{n\mathbf{𝐤}})-{ϵ}_{n\mathbf{𝐤}}|{\varphi }_{n\mathbf{𝐤}}⟩]$

The G₀W₀ method avoids the direct computation of the Green's function and neglects self-consistency in $G$ completely. In fact, only the Kohn-Sham energies are updated from ${ϵ}_{n\mathbf{𝐤}}\to E_{n\mathbf{𝐤}}^{(1)}$, while the orbitals remain unchanged. This is the reason why the G₀W₀ method is internally selected as of VASP6 with ALGO =EVGW0 ("eigenvalue GW") in combination with NELM=1 to indicate one single iteration, even though the method is commonly known as the G₀W₀ approach. To keep backwards-compatibility, however, ALGO=G0W0 is still supported in VAS6.

Note that avoiding self-consistency might seem a drastic step at first sight. However, the G₀W₀ method often yields satisfactory results with band-gaps close to experimental measurements and is often employed for realistic band gap calculations.

Partially self-consistent: GW₀ or EVGW₀

The G₀W₀ quasi-particle energies can be used to update the poles of the Green's function in the spectral representation $G^{(i)}(\mathbf{𝐫},{\mathbf{𝐫}}^{\prime },\omega )=\sum _{n\mathbf{𝐤}}\frac{{\varphi }_{n\mathbf{𝐤}}^{*(0)}(\mathbf{𝐫}){\varphi }_{n\mathbf{𝐤}}^{(0)}({\mathbf{𝐫}}^{\prime })}{\omega -E_{n\mathbf{𝐤}}^{(i)}}$ which in turn can be used to update the self-energy via ${\mathrm{\Sigma }}^{(0)}=G^{(i)}{W}^{(0)}$. This allows to form a partial self-consistency loop, where the screening is kept on the DFT level. The method is commonly known as GW₀, even though only eigenvalues are updated:

• Determine the independent particle polarizability ${\chi }_{\mathbf{𝐆}{\mathbf{𝐆}}^{\prime }}^{(0)}(\mathbf{𝐪},\omega )$
• Determine the screened Coulomb potential $W_{\mathbf{𝐆}{\mathbf{𝐆}}^{\prime }}^{(0)}(\mathbf{𝐪},\omega )$ and keep it fixed in the following
• Determine the self-energy ${\mathrm{\Sigma }}^{(j)}(\mathbf{𝐫},{\mathbf{𝐫}}^{\prime },\omega )=G^{(j)}{W}^{(0)}$.
• Update quasi-particle energies $E_{n\mathbf{𝐤}}^{(j+1)}={ϵ}_{n\mathbf{𝐤}}+Z_{n\mathbf{𝐤}}^{(j)}\mathrm{R}\mathrm{e}[⟨{\varphi }_{n\mathbf{𝐤}}|-\frac{\mathrm{\Delta }}{2}+V_{ext}+V_h+{\mathrm{\Sigma }}^{(j)}(\omega =E_{n\mathbf{𝐤}}^{(j)})-{ϵ}_{n\mathbf{𝐤}}|{\varphi }_{n\mathbf{𝐤}}⟩]$. In the first iteration use $E_{n\mathbf{𝐤}}^{(0)}={ϵ}_{n\mathbf{𝐤}}$

The last two steps are repeated until self-consistency is reached. The GW₀ method is computationally slightly more expensive than the single-shot approach, but yields often excellent agreement with experimentally measured band gaps while being computationally affordable at the same time.

Note that the GW₀ and its single-shot approach do not allow for updates in the Kohn-Sham orbitals ${\varphi }_{n\mathbf{𝐤}}$, merely the eigenvalues are updated. Furthermore, the name GW₀ indicates an update in the Green's function as a solution of the Dyson equation, while the used spectral representation of the Green's function above is strictly speaking correct only in the single-shot approach. Since VASP6 allows to update the Green's function from the solution of the corresponding Dyson equation, the commonly known GW₀ method is also selectable with ALGO=EVGW0 ("eigenvalue GW") and the number of iteration is set with NELM.

Self-consistent Quasi-particle approximation: scQPGW₀

In addition to eigenvalues one can use the self-consistent Quasi-particle GW0 approach (scQPGW0) to update the orbitals ${\varphi }_{n\mathbf{𝐤}}\to {\psi }_{n\mathbf{𝐤}}^{(j)}$ as well. This approach was presented first by Faleev et. al, and used a hermitized self-energy ${\mathrm{\Sigma }}^{\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{m}}=\frac{\mathrm{\Sigma }+{\mathrm{\Sigma }}^{†}}{2}$ in the eigenvalue equation to determine both, quasi-particle energies $E_{n\mathbf{𝐤}}$ and corresponding orbitals ${\psi }_{n\mathbf{𝐤}}$.

In contrast to the Faleev approach one may consider a generalized eigenvalue problem instead that is obtained consistently from the linearization of the self-energy $\mathrm{\Sigma }(E^{(j+1)})\approx \mathrm{\Sigma }(E^{(j)})+\xi (E^{(j)})(E^{(j+1)}-E^{(j)})$ (where $\xi (E^{(j)})=\partial \mathrm{\Sigma }(E^{(j)})/\partial E^{(j)}$) and reads

$\underset{\mathbf{𝐇}(E_{n\mathbf{𝐤}}^{(j)})}{\underbrace{[T+V_{ext}+V_h+\mathrm{\Sigma }\left(E_{n\mathbf{𝐤}}^{(j)}\right)-\xi \left(E_{n\mathbf{𝐤}}^{(j)}\right)E_{n\mathbf{𝐤}}^{(j)}]}}|{\psi }_{n\mathbf{𝐤}}^{(j+1)}⟩=E_{n\mathbf{𝐤}}^{(j+1)}\underset{\mathbf{𝐒}(E_{n\mathbf{𝐤}}^{(j)})}{\underbrace{[1-\xi (E_{n\mathbf{𝐤}}^{(j)})]}}|{\psi }_{n\mathbf{𝐤}}^{(j+1)}⟩$

The resulting Hamiltonian $\mathbf{𝐇}$ and overlap matrix $\mathbf{𝐇}$ are non-hermitian in general, implying that the resulting orbitals $|{\psi }_{n\mathbf{𝐤}}^{(j+1)}⟩$ are not normalized to 1. Therefore, VASP determines the hermitian parts $H=\frac{\mathbf{𝐇}+{\mathbf{𝐇}}^{†}}{2},S=\frac{\mathbf{𝐒}+{\mathbf{𝐒}}^{†}}{2}$ and diagonalizes following matrix instead

$S^{-1/2}H{S}^{-1/2}=U\mathrm{\Lambda }{U}^{†}$

The resulting diagonal matrix $\mathrm{\Lambda }$ contains the new quasi-particles, while the unitary matrix $U$ determine the new orbitals ${\psi }_{n\mathbf{𝐤}}^{(j+1)}=\sum _m{U}_{nm}{\psi }_{m\mathbf{𝐤}}^{(j+1)}$. The method can be selected in VASP with ALGO=QPGW0. See here for more information.

Low-scaling GW: The Space-time Formalism

The GW implementations in VASP described in the papers of Shishkin et al. avoid storage of the Green's function $G$ as well as Fourier transformations between time and frequency domain entirely. That is, all calculations are performed solely on the real frequency axis using Kramers-Kronig transformations for convolutions in the equation of $\chi$ and $\mathrm{\Sigma }$ in reciprocal space and results in a relatively high computational cost that scales with $N^4$ (number of electrons).

The scaling with system size can, however, be reduced to $N^3$ by performing a so-called Wick-rotation to imaginary time $t\to i\tau$.

Following the low scaling ACFDT/RPA algorithms the space-time implementation determines first, the non-interacting Green's function on the imaginary time axis in real space

$G(\mathbf{𝐫},{\mathbf{𝐫}}^{\prime },i\tau )=-\sum _{n\mathbf{𝐤}}{\varphi }_{n\mathbf{𝐤}}^{*(0)}(\mathbf{𝐫}){\varphi }_{n\mathbf{𝐤}}^{(0)}({\mathbf{𝐫}}^{\prime })e^{-({ϵ}_{n\mathbf{𝐤}}-\mu )\tau }[\mathrm{\Theta }(\tau )(1-f_{n\mathbf{𝐤}})-\mathrm{\Theta }(-\tau )f_{n\mathbf{𝐤}}]$

Here $\mathrm{\Theta }$ is the step function and $f_{n\mathbf{𝐤}}$ the occupation number of the state ${\varphi }_{n\mathbf{𝐤}}^{(0)}$. Because the Green's function is non-oscillatory on the imaginary time axis it can be represented on a coarse grid ${\tau }_m$, where the number of time points can be selected in VASP via the NOMEGA tag. Usually 12 to 16 points are sufficient for insulators and small band gap systems.

Subsequently, the irreducible polarizability is calculated from a contraction of two imaginary time Green's functions

$\chi (\mathbf{𝐫},{\mathbf{𝐫}}^{\prime },i{\tau }_m)=-G(\mathbf{𝐫},{\mathbf{𝐫}}^{\prime },i{\tau }_m)G({\mathbf{𝐫}}^{\prime },\mathbf{𝐫},-i{\tau }_m)$

Afterwards, the same compressed Fourier transformation as for the low scaling ACFDT/RPA algorithm is employed to obtain the irreducible polarizability in reciprocal space on the imaginary frequency axis $\chi (\mathbf{𝐫},{\mathbf{𝐫}}^{\prime },i{\tau }_m)\to {\chi }_{\mathbf{𝐆}{\mathbf{𝐆}}^{\prime }}(\mathbf{𝐪},i{\omega }_n)$.

The next step is the computation of the screened potential

$W_{\mathbf{𝐆}{\mathbf{𝐆}}^{\prime }}(\mathbf{𝐪},i{\omega }_m)={[{\delta }_{\mathbf{𝐆}{\mathbf{𝐆}}^{\prime }}-{\chi }_{\mathbf{𝐆}{\mathbf{𝐆}}^{\prime }}(\mathbf{𝐪},i{\omega }_m)V_{\mathbf{𝐆}{\mathbf{𝐆}}^{\prime }}(\mathbf{𝐪})]}^{-1}{V}_{\mathbf{𝐆}{\mathbf{𝐆}}^{\prime }}(\mathbf{𝐪})$

followed by the inverse Fourier transform $W_{\mathbf{𝐆}{\mathbf{𝐆}}^{\prime }}(\mathbf{𝐪},i{\omega }_n)\to \chi (\mathbf{𝐫},{\mathbf{𝐫}}^{\prime },i{\tau }_m)$ and the calculation of the self-energy

$\mathrm{\Sigma }(\mathbf{𝐫},{\mathbf{𝐫}}^{\prime },i{\tau }_m)=-G(\mathbf{𝐫},{\mathbf{𝐫}}^{\prime },i{\tau }_m)W({\mathbf{𝐫}}^{\prime },\mathbf{𝐫},i{\tau }_m)$

From here, several routes are possible. Apart from the single-shot, EVG₀ and QPEVG₀ approximations discussed above, which can also be done on the imaginary axis, the space-time formulation allows to solve the full Dyson equation for $G(\mathbf{𝐫},{\mathbf{𝐫}}^{\prime },i\tau )$ with decent computational cost.

---