# Large scale GW calculations

###### Abstract

We present GW calculations of molecules, ordered and disordered solids and interfaces, which employ an efficient contour deformation technique for frequency integration, and do not require the explicit evaluation of virtual electronic states, nor the inversion of dielectric matrices. We also present a parallel implementation of the algorithm which takes advantage of separable expressions of both the single particle Green’s function and the screened Coulomb interaction. The method can be used starting from density functional theory calculations performed with semi-local or hybrid functionals. We applied the newly developed technique to GW calculations of systems of unprecedented size, including water/semiconductor interfaces with thousands of electrons.

Materials Science Division, Argonne National Laboratory, Argonne IL, USA \alsoaffiliationMaterials Science Division, Argonne National Laboratory, Argonne IL, USA \abbreviationsMG,GG

## 1 Introduction

The accurate description of the excited state properties of electrons plays an important role in many fields of chemistry, physics, and materials science^{1}. For example, the interpretation and prediction of photoemission and opto-electronic spectra of molecules and solids rely on the ability to compute transitions between occupied and virtual electronic states from first principles, as well as their lifetimes^{2}.

In particular, in the growing field of materials for energy conversion processes – including solar energy conversion by the photovoltaic effect and solar to fuel generation by water photocatalysis – it has become key to develop predictive tools to investigate the excited state properties of nanostructures and nanocomposites and of complex interfaces^{3, 4, 5}. The latter include solid/solid and solid/liquid interfaces, e.g. between a semiconductor or insulator and water or an electrolyte^{6, 7, 8, 9, 10}.

In the last three decades, Density Functional Theory (DFT) has been widely used to compute structural and electronic properties of solids and molecules^{11, 12, 13, 14, 15}. Although successful in describing ground state and thermodynamic properties, and in many ab initio molecular dynamics studies^{16, 17}, DFT with both semi-local and hybrid functionals has failed to accurately describe excited state properties of several materials^{18}. However, hybrid functionals have brought great improvement for properties computed with semi-local ones, e.g. for defects in semiconductor and oxides^{19, 20, 21, 22}. In particular hybrid functionals with admixing parameters computed self-consistently have shown good performance in reproducing experimental band gaps and dielectric constants of broad classes of systems^{23}. In the case of the electronic properties of surfaces, interfaces (and hence nanostructures), the use of hybrid functionals has in many instances not been satisfactory. Indeed calculations with hybrid functionals yield results for electronic levels that often depend on the mixing parameter chosen for the Hartree-Fock exchange; such parameter is system dependent and there is no known functional yielding satisfactory results for the electronic properties of interfaces composed of materials with substantially different dielectric properties, as different as those of, e.g. water ( = 1.78)^{24} and Si ( = 11.9)^{25} or water and transition metal oxides of interest for photoelectrodes ( = 5-7)^{26}.

The use of many body perturbation (MBPT) starting from DFT single particle states has proven accurate for several classes of systems^{27, 28, 29, 30, 31, 32, 33, 34, 35, 36} and it appears to be a promising framework to describe complex nanocomposites and heterogeneous interfaces. MBPT for the calculations of photoemission spectra in the GW approximation^{37}, and of optical spectra by solving approximate forms of the Bethe Salpeter Equation (BSE)^{38} is in principle of general applicability; however its use has been greatly limited by computational difficulties in solving the Dyson’s equation and the BSE for realistic systems.

Recently we proposed a method to compute quasi particle energies within the approximation (i.e. the non-selfconsistent approximation) that does not require the explicit calculation of virtual electronic states, nor the inversion of large dielectric matrices^{39, 40}. In addition the method does not use plasmon pole models but instead frequency integrations are explicitly performed and there is one single parameter that controls the accuracy of the computed energies, i.e. the number of eigenvectors and eigenvalues used in the spectral decomposition of the dielectric matrix at zero frequency. The method was successfully used to compute the electronic properties of water^{41} and aqueous solutions^{42} and of heterogeneous solids^{5}, including crystalline and amorphous samples^{40}.

However the original method contained some numerical approximations in the calculations of the head and wings of the polarizability matrix; most importantly the correlation self-energy was computed on the imaginary axis and obtained in real space by analytic continuation. Finally, although exhibiting excellent scalability, the method was not yet applied to systems with thousands of electrons, e.g. to realistic interfaces, due to the lack of parallelization in its original implementation.

In this paper we solved all of the problems listed above, by (i) eliminating numerical approximations in the calculation of the polarizability; (ii) avoiding the use of an analytic continuation and using efficient contour deformation techniques; (iii) providing a parallel implementation of the algorithm based on separable forms of both the single particle Green function and the screened Coulomb interaction. The method presented here can be used starting from DFT orbitals and energies obtained both with semi-local and hybrid functionals. We applied our technique to the calculation of the electronic properties of systems of unprecedented size, including water/semiconductor interfaces with more than one thousand electrons. These calculations allow one to accurately study states at heterogeneous interfaces and to define an electronic thickness of solid/liquid interfaces using MBPT.

The rest of the paper is organized as follows. Sec. 2 describes the methodology that we implemented in a computational package called West. Sec. 3 presents results for the ionization potentials of closed and open shell molecules and for the electronic structure of crystalline, amorphous and liquid systems, aimed at verifying and validating the algorithm and code West against previous calculations and measurements. Sec. 4 presents the study of the electronic properties of finite and extended large systems, i.e. nanocrystals and solid/liquid interfaces, of interest to photovoltaic and photocatalysis applications, respectively. Our conclusions are reported in Sec. 5.

## 2 Method

Within DFT, the -th single particle state and energy of a system of interacting electrons is obtained by solving the Kohn-Sham (KS) equation^{11, 12, 13, 14, 15}:

(1) |

where is the KS Hamiltonian, is the kinetic energy operator, and , and are the ionic, Hartree, and exchange-correlation potential operators, respectively. The indexes and label a wavevector within the first Brillouin zone (BZ) and spin polarization, respectively. Here we consider collinear spin, i.e. decoupled up and down spins.

In a fashion similar to Eq. (1) one may obtain quasiparticle (QP) states and QP energies by solving the equation:

(2) |

where the QP Hamiltonian is formally obtained by replacing, in Eq. (1), the exchange-correlation potential operator with the electron self-energy operator ; is the interacting one-particle Green’s function, is the screened Coulomb interaction and is the vertex operator^{28, 43}. All quantities entering the definition of the self-energy are interdependent and can be obtained self-consistently adopting the scheme suggested by L. Hedin^{44, 45, 46}. In the approximation, is set equal to the identity, which yields the following expression for the electron self-energy^{47}:

(3) |

where is the screened Coulomb interaction obtained in the random phase approximation (RPA). Due to the non-locality and the frequency dependence of the electron self-energy, a self-consistent solution of Eq. (2) is computationally very demanding also for relatively small systems, containing tens of electrons, and usually one evaluates QP energies perturbatively:

(4) | |||||

(5) |

We note that enters both the left and right hand side of Eq. (5), whose solution is usually obtained recursively, e.g. with root-finding algorithms such as the secant method. The use of Eq. (5) to evaluate QP energies from KS states and of the corresponding KS wavefunctions is known as
the approximation.

Within , using the Lehmann’s representation, the Green’s function is:

(6) |

where is a small positive quantity and is the Fermi energy. In Eq. (6) we used the subscript KS to indicate that the Green’s function is evaluated using the KS orbitals obtained by solving Eq. (1).

In Ref. [ 39, 40] an algorithm was introduced to compute the self-energy matrix elements of Eq. (5) without the need to evaluate explicitly empty (virtual) electronic states, by using a technique called projective eigendecomposition of the dielectric screening (PDEP). A diagram of the method is reported in Fig. 1. After KS single particle orbitals and energies are obtained using semilocal or hybrid functionals, the screened Coulomb interaction is computed using a basis set built from the eigenpotentials of the static dielectric matrix at zero frequency. In this way entering Eq. (3) is expressed in a separable form, similar to that of in Eq. (6). In the following sections we describe in detail all the steps outlined in Fig. 1. The separable form of is given in Sec. 2.1. Calculation of the polarizability and the spectral decomposition of the dielectric matrix are described in Sec. 2.2 and 2.3, respectively. Matrix elements of and are then obtained without the explicit use of empty electronic states and simultaneously at several frequencies by using a deflated Lanczos technique, described in Sec. 2.4. Finally the frequency integration is carried out by introducing a contour deformation method, as described in Sec. 2.5. The use of the analytic continuation used in the original method of Ref. [ 39, 40] is thus avoided.

### 2.1 Separable form of the screened Coulomb interaction

In order to solve equation Eq. (5) and obtain QP energies, one needs to compute the matrix elements of the electron self-energy between KS states, which in the approximation is given by Eq. (3). The Green’s function may be expressed in a fully separable form using its Lehmann’s representation, Eq. (6). However is not trivially separable; it is given as the sum of two terms:

(7) |

where is the bare Coulomb interaction\bibnoteWe use everywhere atomic Rydberg units. and is a nonlocal and frequency dependent function. Using Eq. (7), we write the self-energy as the sum of two contributions , where only the latter depends on the frequency:

(8) | |||||

(9) |

is the number of occupied states with spin ; is usually called exchange self-energy because it is formally equivalent to the Fock exact exchange operator\bibnoteThe calculation of the Fock exact exchange matrix elements in reciprocal space is done employing the Gygi-Baldereschi method.^{125};

(10) |

is referred to as correlation self-energy. Using Eq.s (8)-(10) the QP Hamiltonian of Eq. (2) may be expressed as:

(11) |

where is the Hartree-Fock potential operator. The ionic potential is treated within the pseudopotential approach\bibnoteAlthough the presented methodology is general, for the purpose of this work we have considered only normconserving pseudopotentials^{126}..

In this work we express in a separable form by adopting the projective dielectric eigendecomposition (PDEP) technique, proposed in Ref. [ 51- 52], and we use a plane wave basis set to express the single particle wave functions and charge density, within periodic boundary conditions:

(12) |

where is a reciprocal lattice vector, and is the unit cell volume. In Eq. (12) all reciprocal lattice vectors such that are included in the summation. Using a plane wave description also for we have

(13) |

where ( is the Kronecker delta) and

(14) |

In Eq. (14) we have introduced the symmetrized reducible polarizability , related to the symmetrized inverse dielectric matrix by the relation:

(15) |

The symmetrized form of the polarizability is

(16) |

The reducible polarizability is related to the irreducible polarizability by the Dyson’s equation, which within the RPA reads:

(17) |

or in terms of symmetrized polarizabilites:

(18) |

Within a plane wave representation each quantity in Eq. (18) is a matrix of dimension , and in principle can be obtained from via simple linear algebra operations. In practice, the matrices of Eq. (18) may contain millions of rows and columns for realistic systems; for example for a silicon nanocrystal with 465 atoms, placed in a cubic box of edge , 1.5 million plane waves are needed in the expansion of Eq. (12) with . It is thus important to find alternative representations of and reduce the number of elements to compute. One could think of a straightforward spectral decomposition:

(19) |

where and are the eigenvectors and eignvalues of , respectively. Unfortunately this strategy is still too demanding from a computational standpoint, as it implies finding eigenvectors and eigenvalues at multiple frequencies.

A computationally more tractable representation may be obtained using the spectral decomposition of at . As apparent from Eq. (18), eigenvectors of are also eigenvectors of ; the latter is easier to iteratively diagonalize than , and the frequency dependency may be dealt with iterative techniques, starting from the solution at , as discussed in Sec. 2.4. Therefore we proceed by solving the secular equation for only at , generating what we call the PDEP basis set , which is used throughout this work to represent the polarizability :

(20) |

here is a matrix of dimension . In general ,^{51, 52} leading to substantial computational savings.\bibnoteThe approach presented here scales as the fourth power of the system size^{40}, as conventional approaches do. However the computational workload of our method represents a substantial improvement over that of conventional approaches, because of a much more favorable pre-factor: the scaling is instead of , where () is the number of occupied (empty) states, in the number of plane waves and is the number of eigenpotentials used in the PDEP expansion of the static dielectric screening. The functions may be computed by solving the Sternheimer equation^{54}, without explicitly evaluating empty (virtual) electronic states. In addition, turns out to be the only parameter that controls the accuracy of the expansion in Eq. (20). The details of the derivation of the PDEP basis set are given in Sec. 2.3. We note that alternative basis sets, based on the concepts related to maximally localized Wannier functions, have been proposed in the literature to reduce the dimensionality of the polarizability matrices^{55}.

By defining , we formally obtain the desired separable form for :

(21) |

The scaling operation used to define is divergent in the long wavelength limit () and for . However such divergence can be integrated yielding:

(22) |

where

(23) |

In Eq. (23) the integration is evaluated on the region of the enclosing the -point (i.e. ).\bibnoteThe divergence of the Coulomb potential present in Eq. (23) can be numerically removed using spherical coordinates. The specific shape of the BZ is taken into account by using a Monte Carlo integration method.

In the limit, we can now write the matrix elements of using: i) the separable form of of Eq. (22) and ii) the expression of , given in Eq. (6), in terms of projector operators:

(24) |

where

(25) |

and are the projector operator over the occupied and unoccupied manyfold of states belonging to k-point and spin , respectively\bibnoteIn the present formulation we considered integer occupation numbers and only the dielectric response given by interband transitions.. Hence we have:

(26) |

where and ( and ) are contributions to the correlation self-energy originating from occupied (empty) states:

(27) |

(28) |

(29) |

(30) |

We have defined . The quantities , , and entering Eq. (26) are now in a form where iterative techniques (see Sec. 2.4) can be applied to obtain the matrix elements of the correlation self-energy. Moreover, because of the completeness of energy eigenstates (), we may compute all quantities in Eq.s (27)-(30) considering only occupied states. The integration over the frequency domain will be discussed in Sec. (2.5).

### 2.2 Polarizability within the random phase approximation

Here we discuss how to compute the polarizability from within the RPA, in the long wavelength limit (), without explicitly evaluating electronic empty states. The Fourier components of the symmetrized irreducible polarizability are given by the Adler-Wiser expression^{58, 59}, which contains an explicit summation over unoccupied states:

(31) | |||||

where the matrix element

(32) |

is often referred to as oscillator strength; it has the following properties:

(33) | |||||

(34) |

Following Ref. [ 60, 61], we partition the polarizability of Eq. (31) into head (), wings (, or ,) and body ( and ) elements. The limit of the body, which we call , is analytic, while the limits of the head and wings are non-analytic, i.e. they depend on the Cartesian direction along which the limit is performed. The long wavelength limits of the head, body and wings of the polarizability matrix are summarized in Table 1.

Using the PDEP basis set we obtain:

(35) |

(36) |

We can now express all the quantities in Table 1 without any explicit summation over empty (virtual) states:

(37) |

(38) |

(39) |

(40) |

Note that the greek letters and identify Cartesian directions, while the roman letters and label the eigevectors of at , i.e. the elements of the PDEP basis set. We have also defined the auxiliary functions and . Within periodic boundary conditions the position operator is ill-defined and is obtained by solving the linear system

(41) |

where the commutator of the KS Hamiltonian with the position operator includes the contribution of the non-local part of the pseudopotential^{60, 61}. Once is obtained, is computed using Eq. (18).

The quantities required to evaluate Eq. (26) are the following^{62}:

(42) |

(43) |

with ; and the matrix elements of . In Eq. (43) the bold symbols denote matrices of dimension .
In order to compute the matrix elements of the correlation self-energy, Eq. (26), we need to evaluate and , namely the head and body of the operator. These are easily obtained via linear algebra operations from , , and .

By replacing explicit summations over unoccupied states with projection operations, Eq.s (37)-(40) may be evaluated using linear equation solvers and (owing to the completeness of the energy eigenstates) the calculation of polarizabilities is carried out without the explicit evaluation of the virtual states. In a similar fashion one obtains the auxiliary functions in Eq. (41) and the PDEP basis set as described in Sec. 2.3.
We note that other approaches were developed in the literature^{63, 64, 65, 66} to improve the efficiency of calculations
by avoiding the calculation of virtual states, or by limiting the number of virtual states to be computed. However such approaches did not make use of the spectral decomposition of the irreducible polarizability to obtain the reducible polarizability, but instead inverted explicity large matrices. Specifically in Reining et al.^{67} the Sternheimer equation was used to obtain the irreducible polarizability without virtual states and then a plasmon pole model was adopted to compute the dielectric response as a function of frequency. In Giustino et al.^{68} the Sternheimer equation was used as well to obtain the irreducible polarizability without computing virtual states; the polarizability matrix was then inverted numerically and either a plasmon pole model or a a Padé expansion were used to treat the frequency dependence.
In our approach we avoided large matrix inversion by using the PDEP basis set to express all polarizability matrices. Finally, we note that an additional advantage of our approach is that Eq.s (37)-(40) may be computed using a deflated Lanczos algorithm for multiple values of the frequency, as discussed in Sec. 2.4. A Lanczos algorithm was also used by Soininen et al.^{69} to iteratively include local field effects in RPA Hamiltonians and avoid explicit inversion of large matrices. However the authors of Ref. [ 69] computed explicitly virtual states.

### 2.3 Projective dielectric eigenpotential (PDEP) basis set

We now describe in detail how to obtain the PDEP basis set introduced in Eq. (20); each function is computed by the iterative diagonalization procedure, summarized in Fig. 2, the procedure is initiated by building an orthonormal set of basis vectors, e.g. with random components. Then Sternheimer equations are solved in parallel, where the perturbation is given by the -th basis set vector . In particular, given a perturbation , the linear variation of the occupied eigenstates of the unperturbed system may be evaluated using the Sternheimer equation^{54}:

(44) |

Eq. (44) may be iteratively solved using e.g. preconditioned conjugate-gradient methods. The linear variation of the density due to the -th perturbation is obtained within density functional perturbation theory^{70, 71} (DFPT) as

(45) |

The matrix elements of the irreducible polarizability in the space spanned by are given by:

(46) |

where is computed using Eq.s (44)-(45) and assuming that . The matrix is then diagonalized to obtain new basis vectors , and the procedure is iterated using e.g. a Davidson algorithm^{72} (See Fig 3). We note that at each iteration, all Sternheimer problems are independent from each other, thus offering the opportunity to carry out embarrassingly parallel calculations. A description of the parallel operations and data layout will be given elsewhere ^{73}. As a result, the algorithm shows a good scalability up to 524288 cores (see Fig. 4).

As an example we show in Fig. 5 the eigenvalues of the matrix obtained with the PDEP algorithm for the water, silane, benzene and sodium chloride molecules, using KS Hamiltonians with different exchange-correlation functionals. The choice of the functional only affects the most screened eigenpotentials, whereas the eigenvalues corresponding to the least screened ones rapidly approach^{51} zero with a decay similar to that predicted by the Lindhard model^{52}. This indicates that the computation of the least screened eigenpotentials might be avoided and carried out using model functions.

If instead of , one wishes to diagonalize , the potential arising from the rearrangements of the charge density in response to the applied perturbation needs to be included in the definition of the perturbation of Eq. (44)^{74, 75, 61}. In a generalized KS scheme the is given by:

(47) |

where is the fraction of exact exchange (EXX) that is admixed to the semilocal exchange potential. The linear variation of the Hartree potential is

(48) |

and those of the exchange and correlation potentials are given by the functional derivatives:

(49) |

The linear variation of the exact exchange potential (EXX) is expressed in terms of variations of the single particle orbitals

(50) |

We note that calculations including require a double self-consistent procedure (see Fig. 2); hence it is computationally more efficient to iteratively diagonalize first and then obtain the reducible polarizabilty by linear algebra operations. \bibnoteBecause of the RPA, Eq. (18) formally corresponds to the solution of when only the Hartree contribution of Eq. (48) is considered in Eq. (47). We recall that both and are Hermitian operators^{77} and because of Eq. (18) they have the same eigenvectors.

### 2.4 The evaluation of and without empty electronic states using a Lanczos algorithm

The calculation of the correlation self-energy, Eq.s (27)-(30), and of the screening, Eq.s (37)-(40), requires the computation of the matrix elements of , defined in Eq. (25), for multiple values of . For each frequency , given two generic vectors and , we define

(51) |

and

(52) |

Eq. (51) can be easily written in terms of the eigenstates and eigenvalues of the KS Hamiltonian:

(53) |

Eq. (52) may be cast as well in terms of the occupied states and energies, by using the relation and writing (called deflated Hamiltonian)

(54) |

The Lanczos alorithm^{78} is used to obtain a set of orthonormal vectors that are used to recast the deflated Hamiltonian in Eq. (54) into a tri-diagonal form:

(55) |

where

(56) |