FEAST nonlinear eigenvalue algorithm for $GW$ quasiparticle equations

The $G_{0}W_{0}$ approximation quasiparticle equation in atomic units is given as

$$\begin{split}\hat{h}\psi(r)&+\int^{+\infty}_{-\infty}\Sigma(r,r^{\prime},\varepsilon)\psi(r^{\prime})dr^{\prime}=\varepsilon\psi(r),\\ &\mbox{with}\;\;\;\hat{h}=-\frac{1}{2}\nabla^{2}+v_{ext}(r)+v_{H}(r),\end{split}$$

(1)

where $v_{ext}$, $v_{H}$ correspond to the external and the Hartree potentials, respectively. The term $\Sigma(r,r^{\prime},\varepsilon)$ is the exchange-correlation self-energy which contains the many-body quantum effects information. Eq. (1) is clearly a non-local and non-linear eigenvalue equation since the self-energy term contains all of these two difficult-to-handle features. The self-energy in one-shot $GW$ approximation is given by the convolution integral of the Green’s function ($G_{0}$) and the screened Coulomb interaction ($W_{0}$) in energy domain.

$$\Sigma(r,r^{\prime},\varepsilon)=\frac{\mathbf{i}}{2\pi}\int^{+\infty}_{-\infty}e^{\mathbf{i}\omega\eta}G_{0}(r,r^{\prime},\varepsilon+\omega)W_{0}(r,r^{\prime},\omega)d\omega.$$

(2)

An exponential term $e^{\mathbf{i}\omega\eta}$ with a positive infinitesimal $\eta$ is inserted into the integral to make sure that the integral is convergent in the complex plane and physically meaningful. The $G_{0}$ means mean-field one-body Green’s function that could be obtained by a mean-field Hamiltonian, such as DFT, and $W_{0}$ is the screened Coulomb interaction which is determined by the mean-field one-body Green’s function $G_{0}$ again and given by

$$W(r,r^{\prime},\omega)=\int^{+\infty}_{-\infty}\epsilon^{-1}(r,r^{\prime\prime},\omega)v(r^{\prime\prime},r^{\prime})dr^{\prime\prime},$$

(3)

where $v(r,r^{\prime})=1/|r-r^{\prime}|$ is bare Coulomb interaction, and the dielectric function $\epsilon(r,r^{\prime},\omega)$ is calculated as follow.

$$\epsilon(r,r^{\prime},\omega)=\delta(r,r^{\prime})-\int^{+\infty}_{-\infty}v(r,r^{\prime\prime})\chi_{0}(r^{\prime\prime},r^{\prime},\omega)dr^{\prime\prime}.$$

(4)

Using the Randon Phase Approximation (RPA), the irreducible polarizability, $\chi_{0}$, in Eq. (4), is the convolution integral of two mean-field one-body Green’s functions usually with a factor 2 to represent the spin - $\sigma$ for the closed shell systems:

$$\chi_{0}(r,r^{\prime},\omega)=-\frac{\mathbf{i}}{2\pi}\sum_{\sigma}\int^{\infty}_{-\infty}e^{\mathbf{i}\omega^{\prime}\eta}G_{0}(r,r^{\prime},\omega+\omega^{\prime})G_{0}(r^{\prime},r,\omega^{\prime})d\omega^{\prime}.$$

(5)

It can be seen from Eq. (2) to (5) that the self-energy is a energy domain integral that is solely determined by two physical quantities, the bare Coulomb interaction $v$ and the Green’s function $G_{0}$. Therefore, the construction of the Green’s function and the evaluation of the energy integral are essential for $GW$ calculations. For Green’s function, typically, the most commonly used method is the eigenfunction decomposition, another way is to construct the mean-field Hamiltonian – $\hat{H}_{0}$, compute the resolvent (inverse) directly and add the remainder imaginary term.

$$\begin{split}G_{0}(r,r^{\prime},\varepsilon)&=\sum_{n}\frac{\psi_{n}(r)\psi_{n}^{*}(r^{\prime})}{\varepsilon-E_{n}-\mathbf{i}\eta\;sgn(E_{F}-E_{n})}\\ &=(\varepsilon+\mathbf{i}\eta-\hat{H}_{0})^{-1}+R(\mathbf{i}\eta),\end{split}$$

(6)

where $E_{n}$ and $\psi_{n}$ are the eigenvalues and eigenfunctions of $\hat{H}_{0}$, $E_{F}$ is the Fermi energy level of the system $\hat{H}_{0}$ and $R(\mathbf{i}\eta)$ could be computed using the occupied orbitals only.

Using eigenfunction decomposition, Eq. (5) can also be equivalently rewritten as the Lehmann representation. Alternatively, Eq. (5) can be solved analytically and expressed by the Green’s function and the occupied orbitals.

$$\begin{split}\chi_{0}(r,r^{\prime},\omega)=\sum_{\sigma}\sum_{i}^{occ}&\sum_{a}^{unoc}\Bigl{[}\frac{\psi^{*}_{i}(r)\psi_{a}(r)\psi_{i}(r^{\prime})\psi^{*}_{a}(r^{\prime})}{\omega-(E_{a}-E_{i})+\mathbf{i}\eta)}-\frac{\psi_{i}(r)\psi^{*}_{a}(r)\psi^{*}_{i}(r^{\prime})\psi_{a}(r^{\prime})}{\omega+(E_{a}-E_{i})-\mathbf{i}\eta)}\Bigr{]}\\ =\sum_{\sigma}\sum_{i}^{occ}&\psi_{i}(r)\psi_{i}^{*}(r^{\prime})\Bigl{(}G_{0}(r,r^{\prime},E_{i}+\omega)+G_{0}(r,r^{\prime},E_{i}-\omega)\Bigr{)}.\end{split}$$

(7)

In practice, exchange and correlation parts of the self-energy are treated separately, $\Sigma(\varepsilon)=\Sigma^{X}+\Sigma^{C}(\varepsilon)$ , using the following two expressions for the exchange part $\Sigma^{X}$, and the correlation part $\Sigma^{C}$:

$$\begin{split}\Sigma^{X}(r,r^{\prime})&=\frac{\mathbf{i}}{2\pi}\int G_{0}(r,r^{\prime},\varepsilon+\omega^{\prime})v(r,r^{\prime})d\omega^{\prime}\\ &=-\sum_{j}^{occ}\frac{\psi_{j}(r)\psi^{*}_{j}(r^{\prime})}{|r-r^{\prime}|},\end{split}$$

(8)

$$\begin{split}\Sigma^{C}(r,r^{\prime},\varepsilon)&=\frac{\mathbf{i}}{2\pi}\int G_{0}(r,r^{\prime},\varepsilon+\omega)(W_{0}(r,r^{\prime},\omega)-v(r,r^{\prime}))d\omega\\ &=\frac{\mathbf{i}}{2\pi}\int G_{0}(r,r^{\prime},\varepsilon+\omega)W_{0}^{C}(r,r^{\prime},\omega)d\omega.\end{split}$$

(9)

It could be noted that the exchange part (8) has an analytical solution and becomes energy independent, leaving only the correlation part (9) energy dependent (for simplicity, the exponential term and the infinity notations will be omitted from now on).

Several different approaches are available to carefully evalaute the correlation component, including: Analytical Continuation (AC) Rojas et al. (1995); Ren et al. (2012a, b); Wilhelm et al. (2016), Contour Deformation (CD) Godby et al. (1988); Gonze et al. (2009); Blase et al. (2011); Govoni and Galli (2015); Golze et al. (2018, 2020), Plasmon-Pole Models (PPM) Hybertsen and Louie (1986); Soininen et al. (2003, 2005); Kas et al. (2007), and fully analytic method using the Casida equation Casida (1995); Tiago and Chelikowsky (2006); Bruneval (2012); Gao and Chelikowsky (2019). In this work, the CD method is utilized as the main way to calculate the energy integral. The CD method converts the integral on real axis into a contour integral minus an integral on imaginary axis, thus avoiding instability on the real axis due to the poles of Green’s functions:

$$\begin{split}&\Sigma^{C}(r,r^{\prime},\varepsilon)=\frac{\mathbf{i}}{2\pi}\oint G_{0}(r,r^{\prime},\varepsilon+\omega)W_{0}^{C}(r,r^{\prime},\omega)d\omega-\frac{1}{2\pi}\int G_{0}(r,r^{\prime},\varepsilon+\mathbf{i}\omega)W_{0}^{C}(r,r^{\prime},\mathbf{i}\omega)d\omega\\ &=\pm\sum_{m}^{poles}\psi_{m}(r)\psi_{m}^{*}(r^{\prime})W_{0}^{C}(r,r^{\prime},E_{m}-\varepsilon\pm\mathbf{i}\eta)-\frac{1}{2\pi}\int G_{0}(r,r^{\prime},\varepsilon+\mathbf{i}\omega)W_{0}^{C}(r,r^{\prime},\mathbf{i}\omega)d\omega.\end{split}$$

(10)

The summation notation in the first term of Eq. (10) is sum over all the residues of the poles that lie in the contour, and for the second term, the Gauss quadrature could be used to evaluate the imaginary axis integral.

Now, substitute equation $\Sigma(\varepsilon)=\Sigma^{X}+\Sigma^{C}(\varepsilon)$ into Eq. (1) , we get

$$\hat{h}\psi_{i}(r)+\int\Bigl{(}\Sigma^{X}(r,r^{\prime})+\Sigma^{C}(r,r^{\prime},\varepsilon_{i})\Bigr{)}\psi_{i}(r^{\prime})dr^{\prime}=\varepsilon_{i}\psi_{i}(r),$$

(11)

where $\Sigma^{X}$ and $\Sigma^{C}$ is given by Eq. (8) and Eq. (10).

Then, after discretization, the nonlinear eigenvalue equation (1) can be written in the following matrix form:

$$(\mathbf{h}+\mathbf{\Sigma}^{X}+\mathbf{\Sigma}^{C}(\varepsilon))\mathbf{\Psi}=\varepsilon\mathbf{S}\mathbf{\Psi},$$

(12)

There are several ways to discretize and construct the system of matrices Eq. (12), real space finite element method (FEM) is used in this work NES (2022); Kestyn and Polizzi (2020); Li et al. (2024); Li and Polizzi (2024). $\mathbf{S}$ represents then the FEM symmetric positive definite mass matrix.

One way to approximately solve this nonlinear eigenvalue equation is to use the DFT wavefunctions as the eigenvector solutions, $\mathbf{\Psi}\approx\mathbf{\Psi}^{DFT}$. The eigenvalues can then be found graphically by plotting the two expressions on the left and right hand sides of the equal sign of the following expression:

$$\mathbf{A}=\Bigl{\{}~{}\varepsilon\in R~{}|~{}\bra{\mathbf{\Psi}^{DFT}}\mathbf{\Sigma}^{C}(\varepsilon)\ket{\mathbf{\Psi}^{DFT}}=\varepsilon-\bra{\mathbf{\Psi}^{DFT}}\mathbf{h}+\mathbf{\Sigma}^{X}\ket{\mathbf{\Psi}^{DFT}}~{}\Bigr{\}},$$

(13)

where the intersections are the graphical solutions defined by the set $\mathbf{A}$.

Eigenvalue problems in which the coefficient matrices depend nonlinearly on the eigenvalues arise in a variety of applications Betcke et al. (2013); Güttel and Tisseur (2017). Eigenvalues $\lambda$ associated with eigenvectors $\mathbf{x}$ are then solutions of the following general form:

$$\mathbf{T}(\lambda)\mathbf{x}=0,$$

(14)

which includes the common linear eigenvalue problem as a special case, letting $\mathbf{T}(z)=z\mathbf{B}-\mathbf{A}$, as well as the broader polynomial case letting $\mathbf{T}(z)=\sum_{m=0}^{M}z^{m}\mathbf{A}_{m}$. For the QP equation (12), we note that

$$\mathbf{T}(z)=z\mathbf{S}-(\mathbf{h}+\mathbf{\Sigma}^{X}+\mathbf{\Sigma}^{C}(z)).$$

(15)

The FEAST algorithm is a contour integral method initially developed to address linear eigenvalue problems Polizzi (2009); FEA (2020). It has since been adapted for polynomial eigenvalue problems Gavin et al. (2018) and, more recently, for general nonlinear eigenvalue problems (14) Brenneck and Polizzi (2020). The latter can be considered an advancement over Beyn’s contour integration method Beyn (2012), as it incorporates several enhancements from the FEAST algorithm, such as residual inverse iterations for improved convergence and mixed-precision arithmetic. Essentially, the FEAST algorithm and other similar contour integral methods involve converting an eigenvalue problem into solving a series of independent linear systems $\mathbf{T}(z_{j})^{-1}\mathbf{y}$ (i.e. $\mathbf{T}(z_{j})\mathbf{x}=\mathbf{y}$) at particular points $z_{j}$ in the complex plane. When applied to solving the genereral non-linear system (14), the FEAST algorithm involve performing the following contour integrations along $\Gamma$ for $k=0$ and $k=1$:

$$\mathbf{Q}_{k}=\frac{1}{2\pi\mathbf{j}}\oint_{\Gamma}z^{k}\Big{(}\mathbf{X}-\mathbf{T}^{-1}(z)\mathbf{R}\Big{)}{(z\mathbf{I}-\mathbf{\Lambda})}^{-1}\,dz,$$

(16)

where ${\bf X=\{x_{1},\dots\,x_{m}\}}$ represents the eigenvector subspace, and ${\bf R}=\{\mathbf{T}(\lambda_{1}){\bf x}_{1},\dots,\mathbf{T}(\lambda_{m}){\bf x}_{m}\}$ regroups the residuals of the $m$ eigenpairs located within $\Gamma$. In practice, a numerical quadrature is needed to compute the contour integral (16) and the resulting FEAST algorithm is presented in Figure 1. It is important to note the following: (i) while the quality of numerical integration, the number of quadrature nodes, and the size $m_{0}$ of the search subspace can impact the convergence rate, they do not affect the accuracy of the final solutions; (ii) in practice, using a relatively small number of quadrature nodes (around 10) is generally sufficient to achieve convergence to machine precision with very few iterations; (iii) the linear systems that need to be solved at the complex points $z_{j}$ are independent of one another and can be solved in parallel; (iv) importantly, these linear systems can be solved with lower accuracy (using lower arithmetic precision or modest relative residuals, such as $\sim 10^{-2}$ with iterative methods) without affecting the convergence rate of the eigenpairs to machine precision.

Finally, when applying FEAST for solving the QP Eq. (12), it requires replacing the $\varepsilon$ in Eq. (10) with the complex number $z$ ($=\varepsilon+\mathbf{j}\gamma$) at the quadrature nodes. However, due to the nonlinear properties, the process of incorporating FEAST contour integral projection (15) and (16) to $\mathbf{\Sigma}^{C}$ (also $\mathbf{\Sigma}^{X}$) involves the calculations of residues of the poles with the complex energy $z$ deviated from the real axis, which could lead to ill-defined $\mathbf{\Sigma}^{C}(z)$ of Eq. 10 (also $\mathbf{\Sigma}^{X}$ of Eq. (8)). To avoid the complication caused by the confusion between the imaginary parts of FEAST nodes and the imaginary numbers present in $\mathbf{\Sigma}^{C}$ (namely the terms $z+\mathbf{i}\omega$ and $E_{m}-z\pm\mathbf{i}\eta$), we propose handling these two types of imaginary components separately. As a result, hypercomplex number (more precisely, quaternion) should be introduced in the calculations. Moreover, to make the basis of quaternion complete, another dummy imaginary part ($\mathbf{k}\zeta$, where $\zeta=0$) needs to be added. The energy domain of the self-energy is then a hypercomplex space defined by the quaternion:

$$q=\varepsilon+\mathbf{i}\omega+\mathbf{j}\gamma+\mathbf{k}\zeta,\;\;\;q\in\mathbb{H}.$$

(17)

In practice, the domain could be reduced to complex plane for some cases. In the context of considering only real part of the self-energy $\mathbf{\Sigma}^{C}$ in $GW$ complex plane, by taking advantage of the symmetry, only the Green’s function in the second term of Eq. (10) needs the quaternion operations, the first and third imaginary parts will be cancelled out after the $\mathbf{i}\omega$ integration in the end, only left with the imaginary part of $\mathbf{j}\gamma$. Stated otherwise, in the context of considering real part of the self-energy, for Eq. (10), the first term is the normal complex number calculation in $z$ complex plane smearing out the $\pm\mathbf{i}\eta$ term, and the real part of $z$ is used to determine whether if the poles lie in the contour without the need of considering the imaginary part of $z$. For the second term, after computing the Green’s function $G_{0}$ in the hypercomplex space $\mathbb{H}$, only the real and second imaginary part $\mathbf{j}$ are taken from the results to turn it back to the $z$ complex plane, and multiply it with the real part of $W_{0}^{C}$. Consequently, Eq. (10) can be rewritten in FEAST algorithm as:

$$\begin{split}\Sigma^{C}(r,r^{\prime},z)=&\pm\sum_{m}^{poles}\psi_{m}(r)\psi_{m}^{*}(r^{\prime})W_{0}^{C}(r,r^{\prime},E_{m}-z)\\ &-\frac{1}{2\pi}\int\Bigl{[}\Re\Bigl{(}G_{0}(r,r^{\prime},q)\Bigr{)}+\Im^{\{\mathbf{j}\}}\Bigl{(}G_{0}(r,r^{\prime},q)\Bigr{)}\Bigr{]}\cdot\Re\Bigl{(}W_{0}^{C}(r,r^{\prime},\mathbf{i}\omega)\Bigr{)}d\omega,\end{split}$$

(18)

where $\Re$ and $\Im^{\{\mathbf{j}\}}$ stand for the real and the $\mathbf{j}$ imaginary part.

Combining the methods mentioned above, we performed $G_{0}W_{0}$ calculations on a few molecules. In principle, nonlinear FEAST can solve all the occupied states and the LUMO state at the same time, as long as the contour is chosen appropriately. However, it has to be noted that due to the pronounced multisolution behavior of the core electrons in $GW$ Golze et al. (2018, 2020), contouring energy states that are lower than HOMO can sometimes cause the calculation results to not converge. This could be fixed by using different mean-field Hamiltonian other than DFT, such as Hartree-Fock. In this article, we take the advantage of FEAST to choose the contour appropriately and solve the HOMO and LUMO states using DFT-PBE as the starting point.

The following two expressions are used to evaluate the solutions of the nonlinear eigenvalue equation:

$$\begin{split}E_{gs}&=\bra{\psi^{DFT}}\hat{h}+\hat{\Sigma}^{X}+\hat{\Sigma}^{C}(E_{gs})\ket{\psi^{DFT}},\\ E_{nf}&=\bra{\psi^{GW}}\hat{h}+\hat{\Sigma}^{X}+\hat{\Sigma}^{C}(E_{nf})\ket{\psi^{GW}},\end{split}$$

(19)

where $E_{gs}$ means energy obtained by graphical solution method using DFT wavefunction approximations, which is equavalent to (13), $E_{nf}$ means energy obtained by nonlinear FEAST algorithm. The only difference is the use of different wavefunctions, while the Hamiltonians keep unchanged. As a result, if the energies are different, it would be due to the DFT wavefunction approximation made in the GW calculations.

From the Table 1 and 2, it can be seen that the wavefunctions of DFT and $G_{0}W_{0}$ are almost identical at the HOMO level, but the difference is very pronounced at the LUMO level. The deviations in LUMO level are roughly one order larger than that in HOMO level. It is also found that LUMO energy states obtained by nonlinear eigenvalue algorithm are lower than the graphical solutions using DFT wavefunction approximations, which means DFT wavefunction approximations would slightly overestimate the HOMO-LUMO gap. From these results, it could be carefully extrapolated that all the DFT orbitals of occupied states are very close to $GW$’s, but the unoccupied states are different.

To summarize, a nonlinear eigenvalue algorithm for solving $G_{0}W_{0}$ QP equations is presented. We utilized the FEAST algorithm, which is capable of handling both linear and nonlinear eigenvalue problems, to solve the $GW$ QP equation and compute the $GW$ eigenfunctions. By introducing hypercomplex numbers into the FEAST spectral projection process, we effectively circumvent the complications associated with two imaginary numbers. Additionally, to simplify quaternion calculations, we reduce the eigenvalue domain back to the complex plane, making $GW$ calculations more accessible and affordable. Our method was applied to several molecules, revealing that the $G_{0}W_{0}$ LUMO wavefunctions differ from DFT wavefunctions, which challenges the previously accepted consensus. This discrepancy highlights the importance of considering nonlinear eigenvalue problems in $GW$ approximations, an area that may often be overlooked within the community but could greatly benefit from the recent progress made in numerical nonlinear algorithms, such as FEAST. We hope our findings will stimulate interest and lead to the calculation of nonlinear eigenvalue equation in future quantum many-body studies, paving the way for more accurate and comprehensive condensed matter physics analyses.