Koopmans meets Bethe-Salpeter: Excitonic optical spectra from GW

Subscriber access provided by AUBURN UNIV AUBURN

Spectroscopy and Excited States

Koopmans meets Bethe-Salpeter: Excitonic optical spectra from GW-free simulations. Joshua David Elliott, Nicola Colonna, Margherita Marsili, Nicola Marzari, and Paolo Umari J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b01271 • Publication Date (Web): 18 Apr 2019 Downloaded from http://pubs.acs.org on April 25, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Koopmans meets Bethe-Salpeter: Excitonic optical spectra from GW -free simulations. Joshua D. Elliott,∗,†,§ Nicola Colonna,∗,‡ Margherita Marsili,¶ Nicola Marzari,‡ and Paolo Umari

†,§

†Dipartimento di Fisica e Astronomia, Università Degli Studi di Padova, via Marzolo 8, I-35131 Padova, Italy ‡Theory and Simulation of Materials (THEOS) and National Centre for Computational Design and Discovery of Novel Materials (MARVEL), École Polytechnique Fédérale de Lausanne, 1015 Lausanne, Switzerland ¶Laboratoire des Solides Irradiés, École Polytechnique, CNRS F-91128 Palaiseau, France §CNR-IOM DEMOCRITOS, Consiglio Nazionale delle Ricerche--Istituto Officina dei Materiali, c/o SISSA, Via Bonomea 265, 34136, Trieste, Italy E-mail: [email protected]; [email protected]

Abstract The Bethe-Salpeter Equation (BSE) can be applied to compute from first-principles optical spectra that include the effects of screened electron-hole interactions. As input, BSE calculations require single-particle states, quasiparticle energy levels and the screened Coulomb interaction, which are typically obtained with many-body perturbation theory, whose cost limits the scope of possible applications. This work tries to address this practical limitation, instead deriving spectral energies from Koopmanscompliant functionals and introducing a new methodology for handling the screened Coulomb interaction. The explicit calculation of the W matrix is bypassed via a direct

1

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

minimization scheme applied on top of a maximally localised Wannier function basis. We validate and benchmark this approach by computing the low-lying excited states of the molecules in Thiel’s set, and the optical absorption spectrum of a C60 fullerene. The results show the same trends as quantum chemical methods and are in excellent agreement with previous simulations carried out at the TD-DFT or G0 W0 -BSE level. Conveniently, the new framework reduces the parameter space controlling the accuracy of the calculation, thereby simplifying the simulation of charge-neutral excitations, offering the potential to expand the applicability of first-principles spectroscopies to larger systems of applied interest.

1

Introduction

Charge-neutral excitations play a pivotal role in many environmental and technologically relevant applications, including photovoltaics 1,2 and photocatalysis. 3–5 Although in principle quantum-mechanical simulations of excited states offer the potential development of novel technologies, the current trade-off between accuracy and computational viability acts as a barrier to the meaningful investigation of realistic systems of interest. In a simple picture of a charge-neutral excitation, the absorption of a photon promotes an electron to a higher-energy state, creating a hole in the valence manifold. It is well known that describing these excitations with density functional theory (DFT) is incorrect. 6 This is due, both to its inability to describe charged excitations (the domain of many-body perturbation theory) and to capture excitonic effects, which are a consequence of two-body interactions between the excited electron-and-hole pair, that renormalize the energy levels and mix the single-particle transitions. Many-body quantum chemistry methods based on correlated wavefunctions can accurately reproduce experimental data; however, these methods scale poorly with the number of atoms and carry a high computational cost, and viable applications are typically limited to just several tens of atoms. Besides wavefunction-based methods, two additional ab-initio 2


Page 2 of 34

Page 3 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 1: Schematic representation of the usual and the proposed path to the solutions of the Bethe-Salpeter equation. (a) the most common approach, where eigenvalues and charge screening are derived from the GW approximation (b) the method introduced in this work where eigenvalues are computed via Koopmans-compliant DFT and the action of the screened Coulomb potential is calculated directly via iterative minimization. approaches are commonly employed in the computation of charge neutral excitations: the time dependent extension of density functional theory 7,8 (TD-DFT), and, within the Green’s function formalism, the Bethe-Salpeter equation (BSE). 6,9 TD-DFT and BSE simulations both provide better accuracy-viability trade-offs than wavefunction methods. In particular, TD-DFT usually performs well for small molecules, 10 yet it is known to fail in describing Rydberg states and charge-transfer excitations due to an incorrect description of the asymptotic long range exchange. 11,12 Additionally, in solids the long-range exchange should be modulated by the dielectric constant of the system 13,14 making the application of TD-DFT to extended systems more challenging. On the other hand, simulations based on the BSE have been shown to reproduce optical properties of molecules 15–18 and extended systems. 6,19 Whilst attractive in principle, BSE calculations can be more demanding in practice as they typically involve several steps, each of which can require convergence of several parameters to ensure accuracy in the final result. For example, the most common approach is based on a combination of DFT and G0 W0 calculations, 20 as described in Figure 1a. In the 3



first step, a DFT calculation is used to build a set of single particle Kohn-Sham states, ψi . The second step constructs the frequency-dependent screened Coulomb interaction matrix, W (r, r0 ; ω). Next comes the quasiparticle equation, which is similar to the DFT Kohn-Sham equations, but includes the self-energy operator in place of the exchange-correlation potential; this is solved for the set of quasiparticle energy levels, GW . Only in the final step the i BSE is solved for excitonic states, making explicit use of the ψi ’s, GW ’s and the screened i Coulomb interaction matrix at zero frequency, W (r, r0 ; ω = 0). ’s, inThe GW steps (blue in Figure 1a), which provide W (r, r0 ; ω = 0) and the GW i volve extensive and time-consuming convergence testing; in addition they also have an unfavourable O(N 4 ) scaling. These two factors, combined with the intrinsic complexity of the method, limit the size of system to which the GW and GW -BSE methods can be applied. Recently, Koopmans-compliant (KC) functionals have emerged as a reliable approach for calculating quasiparticle energies 21–31 and as an efficient alternative to more complex electronic structure methods. 32 Here, we consider bypassing the GW step entirely by following the approach described in Figure 1b, where the KI functional 27,30–32 is used for the computation of quasiparticle energies, KI i . For the static screened Coulomb interaction, we eliminate the need to compute the W (r, r0 ; ω = 0) matrix explicitly. Instead, we improve on established methods 33 by introducing a new strategy, which applies the action of the screened Coulomb interaction directly onto a basis of Wannier functions, wv (r), via iterative minimization within the BSE step. Not only does this approach remove the need for explicit calculation and storage of the W -matrix elements, but it is also parameter free and therefore drastically reduces the computational time required for convergence testing and could expedite future applications in high-throughput materials modelling. We validate and benchmark this GW -free BSE implementation by computing singlet excitations for the widely studied Thiel’s set 34 and comparing our results with quantum chemical methods, TD-DFT and GW -BSE. Finally, as an example of larger molecular system we compute the absorption spectrum of C60 , which is found to be in good agreement with

4


Page 4 of 34

Page 5 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


experimental measurement and with previous results from GW -BSE.

2

Method

Within many-body perturbation theory, the correlated behavior of electron-hole pairs is described by the Bethe-Salpeter equation: The BSE is a Dyson-like equation for the twoparticle correlation function. 9 When a static approximation for screening is used, the BSE may be recast in the form of an eigenvalue problem involving a two-body Hamiltonian-like ˆ eh . 6,19 Within the Tamm-Dancoff approximation, 6 and expressed in the transition operator H ˆ eh is given by space of the electron-hole (eh) pair, H

eh ˆ vv ˆ ˆx ˆd H 0 cc0 = Dvv 0 cc0 + Kvv 0 cc0 + Kvv 0 cc0 .

(1)

ˆ The indices v & v 0 and c & c0 run over occupied and unoccupied single-particle states, and D, ˆ x and K ˆ d are the diagonal, electron-hole exchange and direct screened Coulomb operators K that we describe in greater detail below. The eigenvalues and eigenstates of the two-body Schrödinger-like equation, X

eh ˆ vv H 0 cc0 Aα,v 0 c0 = Eα Aα,vc ,

(2)

v 0 c0

represent charge-neutral excitation energies Eα and excitonic (excited) states Aα,vc . In the language of second quantization a generic excited state in the Tamm-Dancoff approximation may be expressed as a linear combination of independent particle transitions:

|Θi =

X

Avc a ˆv a ˆ†c |Ψ0 i,

(3)

vc

where |Ψ0 i is the ground-state wavefunction, and a ˆ† and a ˆ are particle creation and annihilation operators. The diagonal operator of (1) is defined as,

Dvv0 ,cc0 = (c − v )δvv0 δcc0 , 5


(4)


Page 6 of 34

ˆ qp . The excitation energy where v and c are eigenvalues of the quasiparticle Hamiltonian H Eα is renormalized according to the strength of the electron-hole interactions through the x d electron-hole exchange (Kvv potentials. 6 These are 0 cc0 ) and direct screened Coulomb Kvv 0 cc0 respectively, x Kvv 0 cc0

ZZ =

dx dx0 ψv∗ (x)ψc (x) vc (r, r0 )ψv0 (x0 )ψc∗0 (x0 )

(5)

dx dx0 ψv∗ (x)ψv0 (x) W (r, r0 )ψc (x0 )ψc∗0 (x0 ),

(6)

and d Kvv 0 cc0

ZZ =−

where the ψi ’s are single particle states, x = (r, σ) is a combined space and spin coordinate, and vc (W ) is the bare (static screened) Coulomb operator. The solution of the BSE as described is highly demanding; in recent years several methodological advancements have been introduced to lower the computational cost of performing such BSE calculations. 15,35,36 We will briefly describe our approach, which is implemented in the Gwl code, 37,38 part of the Quantum ESPRESSO distribution. 39,40 We consider only the case of time-reversal symmetric systems, such that the single-particles states can be taken to be real and a complex-conjugate notation can therefore be neglected. Moreover, we neglect relativistic effects such that ψ(x) factorizes in the product of a spatial function and a spin function and we can retain only the spatial component of the x-coordinate. Adopting techniques based on density functional perturbation theory 41 (DFPT), the explicit calculation of unoccupied single-particle states can be removed, 15,36 and utilizing the so-called batch representation a set of Nv single particle functions 15,36,42

ξv (r) =

X

Avc ψc (r),

(7)

c

are defined for the description of excitonic states such that, in the spatial representation, Z Z |Θi =

† 0 0 ˆ dr dr Ψ(r)Ψ (r )Θ(r, r ) |Ψ0 i, 0ˆ

6


(8)

Page 7 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


with Θ(r, r0 ) =

X

ψv (r)ξv (r0 ).

(9)

v

ˆ † and Ψ ˆ represent the creation and annihilation field operators. In this formulation, the Ψ ˆ |{ξv }i, |{ξ 00 }i = K ˆ x 0 |{ξv }i and |{ξ 000 }i = K ˆ d 0 |{ξv }i actions of the operators |{ξv0 }i = D v v v,v v,v can be written, elementwise, as 36

ξv0 (r) =

ξv00 (r) =

Z

ˆ qp − v I ξv (r), H dr0 Pc (r, r0 )ψv (r0 )

(10) X Z

00 0 00 00 00 dr vc (r , r )ψv0 (r )ξv0 (r ) ,

(11)

v0

ξv000 (r)

Z = −

0

0

dr Pc (r, r )

X

0

ξv0 (r )

Z

dr W (r , r )ψv0 (r )ψv (r ) , 00

0

00

00

00

(12)

v0

where we introduce the operator Pc for projections onto the manifold of empty states, Pc (r, r0 ) =

X

ψc (r)ψc (r0 ) = δ(r − r0 ) −

X

ψv (r)ψv (r0 ).

(13)

v

c

The advantage of this approach is that explicit sums over the empty-state manifold, which enter into the evaluation of W and constitute a significant computational bottleneck, are removed. 15 A further speed-up (and improvement to overall scaling 36 ) can be obtained by noting that, unlike ξv00 that is evaluated in reciprocal space, it is more convenient to compute ξv000 in real space, using a maximally localised Wannier Function (MLWF) representation for the occupied single-particle states. A given valence state is transformed to a MLWF via a unitary rotation, 43 wv (r) =

X

Uvv0 ψv0 (r),

(14)

v0

by the matrix U . In the present implementation we use the algorithm proposed by Gygi et

7



Page 8 of 34

al. 44 The batches, ξv , are also reversibly transformable in the MLWF representation: ξ˜v (r) =

X

Uvv0 ξv (r).

(15)

v0

Thus, we can rewrite Eqn. (12) as

ξ˜v000 (r)

Z =−

0

0

dr Pc (r, r )

X

ξ˜v0 (r ) 0

Z

dr W (r , r )wv0 (r )wv (r ) . 00

0

00

00

00

(16)

v0

Exploiting locality, we define a threshold for which a given pair of MLWFs overlap, Z s

Koopmans meets Bethe-Salpeter: Excitonic optical spectra from GW

Recommend Documents