Toward Efficient GW Calculations Using Numerical Atomic Orbitals

Article pubs.acs.org/JCTC

Cite This: J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Toward Efficient GW Calculations Using Numerical Atomic Orbitals: Benchmarking and Application to Molecular Dynamics Simulations Peter Koval,†,* Mathias Per Ljungberg,† Moritz Müller,† and Daniel Sań chez-Portal†,‡ †

Donostia International Physics Center, Paseo Manuel de Lardizabal 4, 20018 Donostia-San Sebastián, Spain Centro de Física de Materiales, Centro Mixto CSIC-UPV/EHU, Paseo Manuel de Lardizabal 5, 20018 Donostia-San Sebastián, Spain

‡

Downloaded via NOTTINGHAM TRENT UNIV on August 12, 2019 at 18:45:14 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: The use of atomic orbitals in Hedin’s GW approximation provides, in principle, an inexpensive alternative to plane-wave basis sets, especially when modeling large molecules. However, benchmarking of the algorithms and basis sets is essential for a careful balance between cost and accuracy. In this paper, we present an implementation of the GW approximation using numerical atomic orbitals and a pseudopotential treatment of core electrons. The combination of a contour deformation technique with a one-shot extraction of quasiparticle energies provides an efficient scheme for many applications. The performance of the implementation with respect to the basis set convergence and the effect of the use of pseudopotentials has been tested for the 117 closed-shell molecules from the G2/97 test set and 24 larger acceptor molecules from another recently proposed test set. Moreover, to demonstrate the potential of our method, we compute the thermally averaged GW density of states of a large photochromic compound by sampling ab initio molecular dynamics trajectories at different temperatures. The computed thermal line widths indicate approximately twice as large electron−phonon couplings with GW than with standard DFT-GGA calculations. This is further confirmed using frozenphonon calculations.

1. INTRODUCTION Accurate ab initio calculations of electronic excited states are becoming increasingly practicable1−5 due to both the increase of the available computational power and advances in method development. Despite its large computational cost, Hedin’s GW approximation (GWA) has been intensively applied during the last decades. A broad range of applications has been addressed, ranging from a straightforward correction of band gaps6−10 to calculations of electron−phonon couplings11−15 and finite-temperature level broadening.16−20 GWA has been implemented in many electronic structure packages.5,21−37 In the majority of applications, plane-wave (PW) basis sets are used to discretize the spatial degrees of freedom. Among the implementations listed above, there are several ones that use mixed basis sets33,36,38 or localized basis sets.27,30−32,34,35,39−41 The use of local atomic orbitals in conjunction with ab initio pseudopotentials (PP) promises the computational efficiency needed to model the excited states in multiatomic systems41−43 and multiscale simulations.13,41 However, realizations of the GWA comprising both featuresbasis sets of atomic orbitals and PPsare still scarce.31,34,44,45 In this work, we present an implementation of Hedin’s GWA using short-range, numerical atomic orbitals (NAO).46−49 Our implementation is compatible with the popular DFT software suite SIESTA50,51 albeit it can be used in principle with the other similar software such as OPENMX,52 FIREBALL,53 GPAW,54 and PLATO.55 All these packages can use normconserving PPs.56,57 The use of PPs limits the explicit © XXXX American Chemical Society

description of the spatial dependencies to only the valence Kohn−Sham (KS) orbitals. Key features of our GW implementation are the use of atom-centered product basis sets and a contour deformation technique to treat the frequency dependence of the self-energy operator. The quasiparticle energies are obtained in a one-shot manner, i.e., within the so-called G0W0 approximation. The NAOs in SIESTA have a well-defined, tunable spatial cutoff, allowing to exploit locality and achieve simulations of large condensed matter systems43,58,59 with an accuracy comparable to PWs47,60−64 for ground-state DFT calculations. In this work we characterize the accuracy of SIESTA’s NAOs toward the one-shot GWA using two molecular test sets. Another accuracy-limiting factor in PP implementations of the GWA is the usage of PPs constructed for local or semilocal exchange−correlation (xc) functionals even though in the GWA the xc is accounted for by nonlocal, Fock-like operators. This PP inconsistency is typically tolerated in existing realizations of GWA,3,24,31,44,65−67 where the introduced errors are estimated to range from 0.1 to 0.2 eV. We characterize the PP inconsistency by comparing our calculations with an allelectron (AE) reference using similar basis sets, finding discrepancies within the mentioned range. Previously, ab initio molecular dynamics (AIMD) simulations have been used to study the effects of structural fluctuations at finite-temperature on the line shapes of Received: May 6, 2019 Published: June 24, 2019 A

DOI: 10.1021/acs.jctc.9b00436 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

function G0−1(r, r′,ω) is defined. Hedin’s self-energy operator (2) accounts for the nonlocal exchange and correlation operators,3,34,76,78,81 while G0 often originates from an effective one-electron Hamiltonian that only contains a much simpler (frequency-independent) operator Σ0xc to account for exchange and correlation effects. The mean-field propagator G0 can be written as

electronic spectra19,46,68 among other problems.46,51,68,69 The quality of such combined AIMD+excitation simulations depends strongly on the accuracy of the employed excitedstate methods. The computational cost of the excited-state calculations can be substantial and it is often necessary to lower the accuracy requirements to limit the computational cost. Here, the efficiency of our methods allows us to study the effect of temperature on the broadening of quasiparticle levels. By means of constant-temperature AIMD simulations, we found that it is important to use GWA instead of simpler semilocal density functionals to estimate the broadening of quasiparticle levels. In our simulations the broadening of the quasi-particle levels will be connected to the electron−phonon couplings.70−73 In this work we compare the normal-mode coupling strength computed directly, in the frozen-phonon approximation, with that obtained from AIMD simulations. The GW-derived normal-mode coupling strengths are substantially larger than those computed using semilocal DFT, indicating the importance of high-level methods for studying electron− phonon interactions. The rest of the paper is organized as follows. In section 2, we briefly review Hedin’s GWA and a simple theory relating the finite-temperature broadening with the diagonal electron− phonon coupling. The implementation of the one-shot GWA is further detailed in section 3. In section 4, we characterize the accuracy of our implementation by comparing to reference AE GW calculations and experiments. At the end of section 4, we analyze the dependence of the computed electron−phonon strength on the employed electronic structure method, comparing the generalized gradient approximation (GGA) of density functional theory (DFT), the Hartree−Fock approximation (HFA) and the GWA. Finally, we conclude in section 5.

G0−1(ω) = ω − T − Vext − VH − Σ0xc

incorporating the (electron-density-dependent) Hartree potential VH in addition to the kinetic energy operator T and external potential Vext created by the nuclei. The system of eqs 1, 2, and 3 becomes closed and, similarly to other one-particle methods, should be solved to self-consistency.31,39,40,82−84 The self-energy (2) is computationally expensive and one would like to avoid to repeatedly compute it during the self-consistent procedure. Therefore, one tries to improve the initial meanfield propagator (3) by using more accurate approximations to the xc operator Σ0xc. The improvements in the mean-field propagator reduce the size of the corrections described by the self-energy Σ(ω) − Σ0xc, which is now expected to become a small perturbation. The improvements also reduce the importance of self-consistency. In the most favorable case, we compute the self-energy Σ(ω) only once; i.e., we use the noninteracting propagator G0(ω) in place of the interacting G(ω) in the GW self-energy (2). Furthermore, it has been argued in the literature that, in the absence of proper vertex corrections, the results obtained within the pure GWA are not improved and might be even worsen by iteration to selfconsistency, at least in the case of solids.82−84 Thus, although recent studies seem to indicate that the situation might be somewhat different for small molecules,31,39,40,85 it has been customary to stay at the computationally less demanding oneshot GWA level. We will focus on the one-shot GWA in this work. There are indications that the hybrid functionals,9,27,29,67,86,87 or plain Hartree−Fock (HF) Hamiltonian31,40,85 provide a suitable noninteracting Green function G0(ω) in the case of small finite systems. In this work, we focus on the HF starting point. GW calculations for molecules using different levels of selfconsistency9,27,29,31,40,67,85,86 indicate that HF frontier molecular orbitals (MOs) Ψn(r) provide a good approximation to the quasiparticle wave functions, i.e., allow for an accurate representation of the interacting Green function G(ω) in an energy range close to the electron chemical potential

2. THEORY In the GWA, the approximation to the interacting singleparticle Green function G(r,r′,ω) is constructed with a Dyson equation ansatz3,34,74−80 G−1(r,r′,ω) = G0−1(r,r′,ω) + Σ(r,r′,ω) − Σ0xc(r,r′)

(1)

where G0(r,r′,ω) is the electron Green function corresponding to a given mean-field model used as reference. The self-energy Σ(r,r′,ω) is an operator that complements the noninteracting Green function G0(r,r′,ω) and provide a better approximation for the interacting Green function G(r,r′,ω). Σ0xc(r,r′) is an effective static operator approximating the effect of exchange and correlation in the calculation of the mean-field Green function G0(r,r′,ω). In Hedin’s GWA the self-energy is given by Green’s function G and dynamically screened Coulomb interaction W as76,78 Σ(r,r′,ω) =

i 2π

∫ G(r,r′,ω+ω′) W (r,r′,ω′)eiω′η dω′

(3)

G(r,r′,ω) =

∑ n

Ψn(r) Ψ*n (r′) ω − En

(4)

In this diagonal ansatz, only GW poles En vary from the HF eigenvalues E0n, while Ψn(r) still denotes the HF eigenstates. Accepting the diagonal ansatz (4) for the interacting propagator, we can replace the Dyson eq 1 by an eigenvalue equation with fixed eigenvectors Ψn(r)31,85,88

(2)

(HHF + Σc(En))Ψn(r) = En Ψn(r)

where η is a small positive regularization constant. This ansatz gives the name to the approximation. There is a closed form expression for the screened interaction (SI) W(r,r′,ω) in the random-phase approximation (RPA) involving Green’s function G(r,r′,ω) and the bare Coulomb interaction v(r,r′) = |r − r′|−1. Therefore, one can turn eqs 1 and 2 into an explicit computational scheme once the mean-field Green

(5)

and the unknown G0W0 quasiparticle energies En. Multiplying this equation by the eigenvectors Ψn(r) from the left side and integrating the spatial degrees of freedom, we obtain the fixed point equation En0 + ⟨n|Σc(En)|n⟩ = En B

(6) DOI: 10.1021/acs.jctc.9b00436 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article


2.2. GW Corrections of Molecular Levels. Computation of the correlation part of Hedin’s self-energy ⟨n|Σc(ω)|n⟩ in eq 6 turns out to be a bottleneck in GW calculations because of the frequency convolution in the expression for the self-energy (12). Probably the most numerically convenient and mathematically rigorous way of computing the frequency integral in eq 12 is via a contour deformation technique.26,85,88,91 In the complex plane (i.e., allowing ω′ to be complex), the integrand in eq 12 possesses poles due to Green’s function G0(ω+ω′) and due to the SI Wc(ω′), both of which are time-ordered f unctions.74,78,81 Because of the time ordering of Green’s function G0(ω), the poles of the occupied states appear at zn = E0n − ω + iη (i.e., Im zn > 0). The poles due to the unoccupied states appear at zm = E0m − ω − iη (i.e., Im zm < 0) for a given frequency ω. The frequency argument ω at which the self-energy is evaluated changes during the iterative solution of the fixed point eq 6, whereas the poles of the SI Wc(ω′) do not depend on the frequency ω and are located above the real axis in the second quadrant and below the real axis in the fourth quadrant on the complex plane of ω′. Therefore, deforming the integration contour Γ in the complex plane as shown in Figure 1, we can avoid all poles of the SI

This equation can be approximately solved with a Z-factor approach, or, preferably, with an iterative procedure. In this work, we could afford the latter option (see sections 2.2 and 3.2 for details). Note that eq 6 only includes the correlation part of the self-energy Σc(ω) = Σ(ω) − Σx in order to avoid doublecounting of the exchange operator Σx. The exchange operator Σx is already accounted for in the HF eigenvalues E0n. 2.1. Screened Interaction and GW Self-Energy. The screened interaction includes, in addition to the bare Coulomb potential v(r,r′) created at a point r by a charge located at r′, also the induced electric potential vind(r,r′,ω) due to the redistribution of charge in response to the presence of the charge at r′. The SI W (r,r′,ω) = v(r,r′) +

∬ v(r,r″) χ(r″,r‴,ω) v(r‴,r′) dr″ dr‴ (7)

can be constructed via the interacting density response δn(r,ω) function χ (r,r′,ω) ≡ δU (r ′ ,ω) , where δU(r,ω) is an external perturbing potential applied to the system. In RPA, the interacting response function χ (r,r′,ω) = χ0 (r,r′,ω) +

∬ χ0 (r,r″,ω) v(r″,r‴) χ(r‴,r′,ω) dr″ dr‴ (8)

is connected to the noninteracting response function χ0(r,r′,ω), which has an explicit expression in terms of the mean-field energy levels and MOs.89,90 Inserting eq 8 into eq 7 leads to the RPA SI W (r,r′,ω) = v(r,r′) +

∬ v(r, r″) χ0 (r″,r‴,ω) W (r‴,r′,ω) dr″ dr‴ (9)

It is instructive at this point to subtract the bare Coulomb interaction from the SI (9), i.e., to define the induced part of the SI vind(ω). Omitting the arguments from the operator expressions, we define the induced part of the SI vind ≡ W − v and express it via the χ0 operator vind = (δ − vχ0 )−1vχ0 v

(10)

where δ is Dirac δ-function. The splitting of the SI (7) into the bare Coulomb interaction v and the induced part vind(ω) is useful because due to the time-ordering of Green’s function,74,79 the self-energy (2) computed with the bare Coulomb interaction v reduces to a sum over occupied states,78,81 giving rise to the familiar Fock exchange operator Σ x (r,r′) = −

∑ n ∈ occ

Ψn(r) Ψ*n (r′) |r − r′|

Figure 1. Integration contour Γ for computing the integral over the frequency ω′ in eq 12. Depending on the frequency ω for which the self-energy Σc(ω) is computed, several poles of Green’s function G0(ω+ω′) may be located inside of the integration contour Γ (red filled circles).

Wc(ω′) and most of the poles of Green’s function G0(ω+ω′) from appearing inside the integration contour Γ. Integrating along the contour Γ, we split the self-energy Σc(ω) into an integral along the imaginary frequency axis and a sum over residues generated by the poles of Green’s function G0(ω+ω′) inside the contour Γ

(11)

Therefore, the rest of the self-energy (2) is interpreted as the correlation operator Σc(r,r′,ω) =

i 2π

+∞

∫−∞

G0(r,r′,ω+ω′) vind(r,r′,ω′)e iηω′ dω′ (12)

Σc(r,r′,ω) = −

Due to this relation, the induced part of the SI is often called the correlation part of the SI Wc(r,r′,ω) ≡ vind(r,r′,ω). We will also use this nomenclature in the rest of this paper. Having defined the operators of Fock exchange (11) and GW correlation (12), we can set up the quasiparticle eq 5. In the following, we solve the quasiparticle eq 5 within the diagonal approximation, i.e., neglecting the effect of the GW correlation operator (12) on the quasiparticle orbitals Ψn(r).

+

1 2π

+∞

∫−∞

∑ zn(ω) ∈Γ

G0(r,r′,ω + iω′) Wc(r,r′,iω′) dω′

fn Ψn(r) Wc(r,r′,znω) Ψn(r′)

(13)

The set of poles {zωn } that gets into the integration contour Γ depends on the frequency ω at which the correlation energy operator Σc(ω) is computed. Inserting eq 13 into eq 6, we get C


Article

Journal of Chemical Theory and Computation En = En0 −

1 2π

*

+∞

×

where Obs(R,ω) is an instantaneous spectrum for a given geometry snapshot R. Thereby, it is possible to estimate the dependency of the observable Obs(T,ω) on the temperature T of nuclear motion.104,121,122 In eq 17, we assume that the influence of the actual nuclear dynamics, i.e., nonadiabatic effects as well as the effects of quantum motion in the vibrational degrees of freedom, are negligible. Therefore, this treatment only includes the statistical effects due to structural fluctuations and it is only strictly valid at high temperature T. In this work, we will focus on the photoelectron spectra of a molecule in the gas phase. The probability of the photoemission is estimated by the density of states (DoS) of occupied levels, while the probability of the inverse photoemission is estimated by the DoS of unoccupied levels.27,34,39,78,80,81,103 Hedin’s GWA provides a high-quality estimation of the DoS, while GGA and HFA will be contrasted against the GWA reference. Because one-shot, diagonal HF +G0W0 only corrects the single-particle energy levels En[R], the thermal average of the DoS is given by

∬ Ψ*n (r) Ψn(r′)

Ψ (r) Ψ (r′) ∫−∞ ∑ E m+ iω′ m− E 0 Wc(r,r′,iω′) dω′ dr dr′ n

m

+

m

∬ Ψ*n (r) Ψn(r′) ∑

×

fk Ψk(r) Wc(r,r′,zk(En)) Ψ*k (r′) dr dr′

z k(En) ∈Γ

(14)

where we used the Lehmann expression for the HF propagator G0(En+iω′) (cf. eq 4). Treating the last expression, it is beneficial to define the (frequency-dependent) electron−hole matrix elements of the SI I nm(ω) =

∬ Ψ*n (r) Ψm(r) Wc(r,r′,ω) Ψn(r′) Ψ*m(r′) dr dr′ (15) nm

The matrix elements I (ω) facilitate the iterative procedure of solving the fixed point eq 15 because the spatial degrees of freedom get integrated out En =

En0 +

1 − 2π

+∞

∫−∞ ∑ m

S −1

DoS(T ,ω) = (Sπ ) Im ∑ ∑ (ω − En[R s(T )] + iε)−1 s=1

I nm(iω′) dω′ En + iω′ − Em0

(18)

where S is the number of snapshots in the ensemble and ε is the width of quasiparticle levels due to electronic relaxation. The broadening constant ε accounts for an intrinsic level width due to the radiative decay and Auger-like processes.102,104,123,124 We will set the electronic level width ε to a small value in order to reveal the finite-temperature effects associated with structural fluctuations. According to experimental evidence,98,99,101,102,125,126 the thermal agitation motion leads to a Gaussian line shape with a full-width-at-half-maximum (FWHM) proportional to the square root of the finite temperature: FWHM ∼ T . This phenomenology is similar to solvent-induced broadening.127 We rationalize this expectation within the framework of statistical physics, treating the ionic motion classically in analogy with our AIMD simulations. Moreover, we assume a harmonic model of ionic motion in the vicinity of equilibrium and a linear dependence of the quasiparticle levels on the normal-mode displacements Qm. In normal coordinates {Pm, Qm},123,128,129 the Hamiltonian in the harmonic model of the nuclear motion is a sum of independent contributions

nk

∑

fk I (zk(En))

(16)

z k(En) ∈Γ nk

Unfortunately, the calculation of matrix elements I (zk(En)) = ⟨nk|Wc(zk(En))|kn⟩ in the vicinity of the real frequency axis must be repeated during the iterations because the set of poles {zk(En)} depends on the current estimate of the quasiparticle eigenenergy En. We will detail the numerical procedures to realize eqs 10, 15, and 16 in section 3. 2.3. Theory of Finite-Temperature Classical Broadening. The thermal motion of atoms affects the observed properties of organic molecules in a significant, sometimes decisive manner. Chemical reactions,92 charge transport,93,94 isomerization states,95−97 and spectral line shapes98−102 are strongly influenced by the thermal motion. For example, the spectral line shapes in photoabsorption, fluorescence, and photoemission experiments are broadened due to the finitetemperature motion, limiting the resolution of spectroscopic experiments.98,101,103,104 Theoretical modeling of the effects of atomic thermal motion gives rise to combined approaches ranging from simple AIMD68,105−108 to state-of-the-art approaches such as Ehrenfest dynamics,68,107,109,110 nonadiabatic MD,68,107,111−113 and Monte Carlo nuclear ensembles.114−116 In this work we present AIMD simulations within the Born− Oppenheimer approximation (BOA) that can be realized using the distributed version of the SIESTA package.117 In AIMD the nuclei are propagated using Newton’s equations with forces computed ab initio. The temperature is usually controlled by a thermostat,118−120 ensuring a thermal distribution of the molecular configurations in time-averaged samples. To derive the statistical average of a spectroscopic observable, Obs(T,ω), one can take the average over a suitable, temperaturedependent set of S molecular geometry snapshots Rs(T) originating from AIMD Obs(T ,ω) =

1 S

3N

H=

s=1

∑ m=1

Q Pm 2 + ωm 2 m 2 2

2

(19)

where N is the number of atoms and ωm are the frequencies of the normal modes. In the following, we assume the mutual independence of the different quasiparticle levels, which is usually a good approximation if the energy level separation is sufficiently large. Thus, the energy E of each quasiparticle can be written as 3N

E = E(Q) = E0 +

∑ αmQ m m=1

(20)

where E0 is the quasiparticle energy at equilibrium geometry and αm are normal-mode coupling constants. We will also use boldface P and Q to denote sets of 3N normal-mode momenta {Pm} and coordinates {Qm}. The ensemble average of ⟨DoS(ω)⟩ is calculated from the statistical sum123,130

S

∑ Obs(R s(T ),ω)

n

(17) D


Article

Journal of Chemical Theory and Computation ⟨δ(ω−E(Q))⟩ =

In this equation and in the following, we use Einstein’s summation convention over repeated indices unless the indices appear on the left-hand side of an equation. It is customary for atomic orbitals fa(r) to be expressed as products of radial and angular components fa(r) = fa(r) Ylama(r̂), where Ylm(r̂) is a spherical harmonic function131 depending only on the direction of the vector r̂ and fa(r) is a function of the radial distance r. In this work, we employ NAOs as implemented in the SIESTA software suite50,51 to represent the radial functions fa(r). Independently, on the type of radial orbitals fa(r), we have to treat the products of atomic orbitals fa(r) f b(r) in GW calculations. The need for the atomic-orbitals’ products can be seen already in the calculation of the Fock exchange operator (11). Calculating the matrix elements of the Fock exchange operator between NAOs fa(r), we insert the LCAO ansatz (26) into the expression (11) and obtain

∫ dQ dP δ(ω−E(Q )) exp(−H /kT ) ∫ dQ dP exp(−H /kT ) (21)

where k is the Boltzmann constant. Performing the integrations over normal coordinates and momenta P and Q in the last equation (21) with the assumptions (19) and (20), we get the normal Gaussian distribution (see the Supporting Information) 2

⟨δ(ω−E(Q))⟩ = (2πσ 2)−1/2 e−(ω − E0)

/(2σ 2)

(22)

where the square of the standard deviation σ 2 = kTV ep

(23)

is proportional to the temperature T and to the diagonal matrix element of electron−phonon couplings to vibrational modes Vep 13,72,73 3N

V ep =

∑ αm2/ωm2 m=1

Σab x = −

(24)

(27)

f a (r) f b (r) = V μabAμ(r)

(28)

There are several procedures to construct the product basis (PB) functions Aμ(r) and the complementary product vertex 41,43,132,133 coef f icients Vab μ . In this work, we continue to use μ atom-centered product functions A (r) with factorized product vertex coefficients

(25)

ab

V μab = Vν̃ cμν

In this equation, we take into account that the ensembleaveraged quasiparticle energy ⟨E⟩ depends on the finite temperature T.

(29)

ν where Ṽ ab ν is the so-called dominant product vertex and cμ is a matrix of expansion coefficients of dominant products in terms of atom-centered product functions Aμ(r). The utility of the factorization (29) is in a minimal number of nonzero elements in the dominant product vertex Ṽ ab ν as explained in subsections 3.2 and 3.3. In refs 41 and 133 we described carefully our product basis and the advantage of combining an expansion in terms of dominant product functions linked to each atom pair, which are finally expressed in terms of atom-centered product functions Aμ(r). This basis set has been used to compute the optical response at the TDDFT level for systems containing up to several thousands of atoms.43,132,134 Inserting the expansion (28) into the Fock operator (27), we get

3. METHODS The one-shot GW calculation sketched above in section 2 depends on the choice of a basis set, representing the spatial dependency of the involved operators and on the discretization of their frequency dependence. In this work, we rely on a linear combination of atomic orbitals (LCAO) to expand the MOs and their products, and on the contour deformation technique to realize the integration over the frequency variable. The contour deformation technique was covered above in section 2.2. In this section, it remains only to specify the numerical grid used for the quadratures along imaginary frequency axis (16) which is done in subsection 3.4. In subsections 3.1, 3.2, and 3.3, we describe the treatment of real-space variables by means of NAOs and an efficient method of computing the response function. Moreover, we report the methods to realize AIMD simulations and the frozen-phonon approximation for the normal-mode coupling coefficients (20) in the Supporting Information. 3.1. Numerical Atomic Orbitals in GW. The gist of LCAO is to expand the MOs Ψn(r) into a linear combination of localized, atom-centered functions fa(r −Ra) Ψn(r) = Xanf a (r−R a)

∬

f a (r) f a′ (r) f b′ (r′) f b (r′) dr dr′ |r − r′|

In the LCAO formalism, the last equation leads to the notion of four-center integrals. However, there is a well-known way to avoid explicit computation of the four-center integrals, speeding up various quantum chemical methods, by using an auxiliary basis set {Aμ(r)} that is capable of expanding the products of NAOs fa(r) f b(r)

S s=1

Xan′Xbn′

n ∈ occ

For a more detailed derivation we refer the reader to our Supporting Information. The FWHM for the Gaussian distribution (22) is proportional to the standard deviation FWHM = σ 8 ln 2 , which complies with the experimental evidence. In section 4, we will compare the matrix elements of electron−phonon couplings Vep computed using the frozenphonon approximation with those extracted from AIMD. Hence, we equate the second moment of the Gaussian distribution (22) M2(T) = σ2 to the second moment of the AIMD DoS M2[AIMD] = S−1∑Ss=1(Es − ⟨E⟩)2 and obtain V ep = (kTS)−1 ∑ (Es − ⟨E⟩)2

∑

aa′ μν b ′ b Σab x = − Da ′ b′V μ v V ν

(30)

where we used the matrix elements v interaction v μν =

∬

μν

of the bare Coulomb

Aμ(r) Aν (r′) dr dr′ |r − r′|

(31)

∑n∈occXnaXnb.

and the density matrix Dab = 3.2. Calculation of the GW Correction. To calculate the matrix elements of the SI (15), we insert the LCAO (26) and the atomic-orbital product expansion (28) into eq 15

(26) E


Article

Journal of Chemical Theory and Computation μν mn Icnm(ω) = = nm μ Wc (ω)= ν

(32)

where we introduced the product vertex between MOs n ab m = nm μ = Xa V μ Xb

(33)

= nm μ

The MO product vertex is a dense table that would take O(N3) elements of memory. To avoid memory depletion, we store only chunks of this table during computations. Moreover, the factorization (29) is used in order to reduce number of mathematical operations and speedup the calculation. The matrix elements of the SI (10) between the PB functions Aμ(r) entering (32) read40,44 Wcμν(ω) = (δνμ′ − v μμ′χμ0′ ν′ (ω))−1v ν ′ μ ″χμ0″ ν ″ (ω)v ν ″ ν

δμν′

Figure 2. Algorithm to compute the response function for the atomcentered PB using sum-over-states expression (35). The critical part of the algorithm is the computation of the MO product vertex = mn μ . It is computed using the factorization of the atom-centered product vertex according to eq 29.

(34)

μν

where is identity matrix, v are the matrix elements of bare Coulomb interaction (31), and χ0μν(ω) are the expansion coefficients of the noninteracting response function in terms of PB functions Aμ(r). The expansion coefficients χ0μν(ω) can be conveniently calculated using a sum-over-states expression34,41,80,88,133 and the MO product vertex (33) (cf. refs 41 and 133) 0 (ω) = χμν

∑ m,n

most cases. On the contrary, working with each electron−hole pair individually (i.e., using chunks of unit size) would limit the numerical efficiency. For these reasons, the electron−hole pairs are split into chunks of adjustable size. We found that the speed of the matrix operations saturates at the chunk size about 500 while using the MKL 11.2 BLAS library on Intel Xeon processors X5550 with 8 MB of cache memory. The sum-over-states approach is convenient in the context of summation over poles in eq 16 because one can get the SI matrices Wμν({ωi}) for any set of frequencies {ωi} without the need to discretize the whole frequency range. For the computation of matrix elements of the SI (15) along the imaginary axis iω′, the set of frequencies ω′ is known prior to the computation. The quantities along the imaginary axis are rather smooth, and we found that it is accurate and convenient to use a logarithmic grid to organize the corresponding integrations. We discuss our choice of logarithmic-grid parameters in the following. 3.4. Frequency Discretization in Contour Deformation Technique. Since the frequency dependence along both frequency- and time-imaginary axes is featureless and requires long extents and higher density at the origin, we followed the well-known practice of discretizing the frequency dependence with a logarithmic numerical grid. We use Talman’s parametrization of the logarithmic grid138

nm (fn − fm )(= mn μ =ν )

ω − (Em − En)

(35) μ

As a rule of thumb, the size of the product basis {A (r)} is only 4 times larger than the size of the NAO basis {fa(r)}. With this size of the product basis and for the relatively small molecules we aim at in this work, the matrix inversions in eq 34 take a small part of runtime comparing to the computation of the response matrix (35). 3.3. Calculation of Response Function. The contour deformation technique reduces the amount of computation provided we are sufficiently satisfied with the one-shot GW and using only the diagonal part of the self-energy. Sticking to the contour deformation technique and to the atom-centered PB {Aμ(r)}, we have to compute the response function χ0μν(ω) on a fixed grid along the imaginary-frequency axis and in the vicinity of the real-frequency axis at some few frequencies unknown in advance. There are several algorithms to compute the response function along the imaginary-frequency axis26,135 and in the vicinity of the real-frequency axis.44,66,136,137 In this work, we have compared essentially the time-frequency algorithm by Godby et al. for the response function135 and also, following the advice by Blase et al.,85 the sum-over-states expression (35). It turns out that a sophisticated, lowcomplexity algorithm is not justified for any small- and medium-sized molecules for the compact NAO basis sets we use. As a result, we also advocate for the simplest sum-overstates approach in this work. In Figure 2, we sketch an O(N4) algorithm, which is typically 5−10 times faster than the time-frequency domain algorithms for many molecules of practical interest, containing up to 100 atoms at least. As one can see in this scheme, we use a double loop with the inner part going through the frequency variable ω and outer part through the electron−hole pairs {mn}. In the outer loop, we compute chunks of the MO vertex (33) using also the factorization of the atom-centered vertex (29). The factorization (29) is important for computational efficiency as reported in detail in refs 41 and 133. The outer loop must generally operate with chunks of electron−hole space {mn} because the total size of the MO vertex is prohibitively large in

tn = exp(log(tmin) + Δt(n − 1)) ωn = tn(ωmax /tmax ), n ∈ [1 ... Nt ]

(36)

where Δt = (log(tmax) − log(tmin))/(Nt − 1). This parametrization has four independent parameters (tmin, tmax, ωmax, and Nt), and it is suitable to discretize both frequency and time-domain quantities. Such Fourier-transform suitability is convenient for the time-frequency algorithms and it does not interfere with the frequency-domain approach we finally pursue (see subsection 3.3). The choice of the grid parameters is based on the properties of the response function along the imaginary time axis it. Along the imaginary time axis, the response function is given by a sum of exponential functions of negative arguments,135 which makes the choice of the time grid parameters straightforward. The minimal (maximal) time tmin (tmax) is reciprocal to a maximal (minimal) frequency Ωmax (Ωmin) present in response function. We estimate the maximal frequency by the maximal energy difference present in the initial mean-field eigenenergies Ωmax = EN − E1. Analogously, the minimal frequency is F


Article

Journal of Chemical Theory and Computation estimated by the HOMO−LUMO gap Ωmin = ELUMO − EHOMO. For a consistent estimation of tmin and tmax, we introduce a tolerance parameter ϵtol

pV5Z basis omitting higher angular momenta by cc-pV5Zl2, i.e., adding the suffix l2. For both AE basis sets cc-pV5Z and cc-pV5Zl2, the calculations were done with the open-source package MOLGW.143 Calculations with lower cardinal numbers ζ < 5 were also done in our basis set convergence study. The results for the cc-pVQZ are available in our Supporting Information. In order to assess the effect of using GGA PPs, we constructed mixed basis sets consisting of the first-ζ NAOs and other valence and polarization functions acquired from ccpVζZl2 basis sets. The first-ζ NAOs are the solution of the radial atomic KS equation, providing an accurate replacement of the (valence) Gaussian contractions. The mixed, numerical−Gaussian basis sets are denoted as NGζ, matching accurately the reference all-electron basis sets cc-pVζZl2. In the main text of the paper, we compare results obtained with quintuple basis sets, i.e., cc-pV5Zl2 and NG5, while the results for smaller cardinal numbers 2 < ζ < 5 are made available in our Supporting Information. The first-ζ NAO is generated with an energy shift parameter50 of 1 meV. The small energy shift translates into extended NAOs, the radius of which varies from one chemical species to the other, ranging from a minimum of 4.97 to a maximum of 14.91 Bohr for F and Na, correspondingly. The NGζ basis set contains valence-polarization functions that have considerably larger spatial extent than the extent of regular NAOs used in SIESTA ground-state calculations. The short spatial support of regular, purely numerical NAOs allows for efficient calculations, but introduces simultaneously inaccuracies. In order to gauge the accuracy of standard NAO basis for GWA, we present calculations with regular double-ζ-polarized (DZP) basis sets generated with the splitvalence method,50 using the energy shifts 1 and 50 meV. Split basis sets are routinely available in the SIESTA package and provide good results for ground-state calculations. These basis sets consist of basis functions shorter than the NGζ basis set, which positively affects the computing runtime. For example, the split DZP with the energy shift of 50 meV can be 8−15 times faster than the smallest NG3 for GW calculations in the Org24 test set. The HF self-consistent calculations have been performed using our own code,40,41,133 starting from the PBE eigenstates provided by SIESTA. We tested two types of PPs: Troullier− Martins56 and Hamann PPs.57,144 Since the differences (in GW IPs and EAs) between Troullier−Martins and Hamann PPs are minimal, we present only the results for the widely available Troullier-Martins PPs. The initial HF and following GW calculations were done with the atom-centered PB discussed in subsection 3.1. We gave a detailed description of the generation procedure elsewhere.41,133 During the generation procedure, we construct orthogonal linear combinations of the NAOs’ products via a diagonalization-based procedure.145 The importance of a particular linear combination is estimated by its corresponding eigenvalue.145 In the construction of the PB, we discard unimportant linear combinations that have eigenvalues lower than a threshold λPB. In subsections 4.1 and 4.2, the threshold was set to the notably small value λPB = 10−8, providing a highaccuracy PB set. For instance, the PB sets for the NG5 basis set contain a minimum of 164 and a maximum of 256 functions for H/Na and B atomic species, correspondingly. Control calculations with a 10 times larger threshold λPB = 10−7 were done for the G2/97 test set. The maximal absolute deviation of

tmax = −log(ϵtol)/Ωgap tmin = −100log(1 − ϵtol)/Ω max

(37)

The extent of the frequency grid ωmax is set according to asymptotic property of the integrand in eq 13 ωmax = (Ω max 2 /ϵtol)1/4

(38)

This estimation is derived using the frequency dependence of single-pole, model functions for G(iω) and W(iω) along the imaginary axis, assuming the pole is at Ωmax in both cases. For the SI W(iω) and Re G(iω) these model functions will be different by factor 2: W(iω) ∼ 2Ωmax/(ω2 + Ωmax2), and Re G(iω) ∼ Ωmax/(ω2 + Ωmax2). In the limit of large frequencies ω′ ≫ Ωmax and ω′ ≫ ω, the asymptotic expression of the integrand G(ω+iω′) W(iω′) = 2Ωmax2/ω′4, which explains eq 38. The tolerance estimator ϵtol gives approximately the accuracy of the GW corrections computed in the time-frequency approaches provided the number of grid points Nt is sufficiently high. Typically, accurate results are achieved with a moderate number of points Nt = 32−64. A similar number of points was used by other authors.135,139 The number of points Nt only weakly affects the results and does not require any system-size-dependent adjustments. The imaginary offset of the poles’ location zp for computing residues of the SI (see section 2.2) is chosen equal to the minimal frequency along imaginary frequency axis 0zp= ±ωmin.

4. RESULTS Before using NAOs in an AIMD+GW simulation in subsection 4.3, we assess the performance of small basis sets of short-range NAOs with limited angular momentum content, and we gauge the effect of using GGA PPs in subsections 4.1 and 4.2. These assessments are achieved by calculations of ionization potentials (IP) and electron affinities (EA) of two molecular sets. The IP and EA are computed as the negative of HOMO and LUMO of HF+G0W0 energy levels, correspondingly. The first set of molecules comprises all 117 closed-shell molecular species in the G2/97 test set.140 The closed-shell molecules of the G2/97 test set consist of 2−14 atoms. Since many of these small molecules have negative EAs, we analyze only their IPs. We refer to the second set of 24 molecules as the Org24 test set, because of the relevance of these molecules in organic electronics as electron acceptors.4 The molecules of the Org24 test set consist of 8−24 atoms. For the Org24 test set, we evaluate both IPs and EAs. For the G2/97 test set, we fixed the geometry of the molecules to the original geometries of the G2/97 test set141 optimized at MP2 level, whereas for the Org24 test set, we used the geometries specified in the SM of ref 4. We use Dunning’s correlation-consistent, polarized-valence quintuple-ζ basis set142 (cc-pV5Z) as the reference basis for both G2/97 and Org24 test sets. To assess the effect of using GGA PPs, we use the cc-pV5Z basis set without functions of angular momentum l > 2. The omission of basis functions with l > 2 enables a direct comparison with the readily available NAO basis sets in the SIESTA package. We denote the ccG


Article


Figure 3. Discrepancies between ionization potentials (eV) of the closed-shell molecules in the G2/97 test set calculated at the level of HF+G0W0. The mean signed errors are computed with respect to the AE cc-pV5Z reference and summarized on panel a. The mean absolute errors and maximal absolute deviations are shown on panel b in lower left and upper right corners of the heatmap, correspondingly.

Figure 4. Discrepancies between ionization potentials (eV) calculated at the level of HF+G0W0 for the Org24 test set. The mean signed errors are computed with respect to the AE cc-pV5Z reference and shown in panel a. The mean absolute errors and maximal absolute deviations are shown in panel b.

the IPs obtained with λPB = 10−7 is below the accuracy of reported numbers, i.e., ≤10−3 eV. The frequency grid along the imaginary axis is logarithmic with automatically chosen parameters defined by the tolerance estimator ϵtol, as explained in subsection 3.4. The iterative solution of fixed point eq 6 is converging robustly with small values of the tolerance estimator ϵtol = 10−10. The small value of the tolerance estimator ϵtol translates to small minimal frequencies ωmin and large maximal frequencies ωmax. For instance, in the G2/97 test set, the tolerance estimator ϵtol = 10−10 for calculations with NG5 basis translates to ωmin in the range (2.55−5.99) × 10−7 eV for hydrogen fluoride and diatomic hydrogen, correspondingly, while the maximal frequencies ωmax vary in the range 1155−1933 eV for nitrogen and fluorine dimers, correspondingly. It is worth noting that these huge frequency ranges along the imaginary frequency axis can be accurately discretized by a few dozens of frequency grid points owing to the smooth character of the relevant functions. The computed IPs and EAs are available in the Supporting Information. Here we discuss exclusively the statistical deviations within the calculated data sets. For instance, we 1 N examine the mean signed error (MSE = N ∑i ei , where ei = xi − xref i ), the mean absolute error (MAE = 1 N ∑ | e | ), and the maximal absolute deviation (MxAD = i i N max(|ei|)), where N is the size of the test set, xi are values in the examined data set and xref i are reference values.

4.1. Ionization Potentials of the G2/97 Test Set. In Figure 3, we illustrate the MSE, MAE, and MxAD for the G2/ 97 test set. The MSE is shown with a stair chart on panel a. The calculations with AE cc-pV5Z provide reference data. For the sake of the readability of the figures, the basis sets cc-pV5Z, cc-pV5Zl2, DZP[1 meV], and DZP[50 meV] are denoted by cc5 and cc5l2, 1, and 50 meV, respectively. The omission of higher angular momenta in the basis set leads to underestimation of IPs by about 0.31 eV (column cc5l2 on panel a). The cc5l2 and NG5 basis sets are the most similar by constuction and the comparison of the respective calculations allows to gauge the effect of using PPs as compared to AE computations. The MSE between cc5l2 and NG5 is the smallest in this set of discrepancies and amounts to only 0.05 eV. The corresponding MAE (0.08 eV) is larger than MSE, indicating the deviations to both sides of the reference data, but still reflecting a mild effect (in average) of the use of PPs in the IPs of the G2/97 test set. The use of purely numerical NAOs, which are routinely available in the SIESTA package for ground-state calculations leads to simplified models and significantly cheaper simulations at the cost of further deviation from the high-quality AE and PP references. For example, a DZP basis with an energy shift of 1 meV gives rise to MSEs of 0.35 and 0.71 eVs with respect to the NG5 and cc5 reference data. Further limitation of the NAOs’ spatial extent in the DZP[50 meV] basis set causes H


Article


Figure 5. Discrepancies between electron affinities (eV) calculated at the level of HF+G0W0 for the Org24 test set. The mean signed errors are computed with respect to the AE cc-pV5Z reference and shown in panel a. The mean absolute errors and maximal absolute deviations are shown in panel b.

the values reported by other authors.66 This assessment sets the upper bar of the errors due to the incompatibility of the GGA PP and Fock-like exchange−correlation operators of the GWA to be below 0.2 eV. Similarly to the G2/97 test set, the use of short-range, low-ζ basis sets causes underestimation of the IPs and EAs. Panel a in Figure 4 shows that the use of the DZP[1 meV] basis set causes the deviation of IPs with the MSEs of 0.25 and 0.59 eV with respect to the NG5 and cc5 references, correspondingly. In turn, panel a in Figure 5 shows that the EAs deviate by a sligthly larger amount with MSEs of 0.32 and 0.64 eV with respect to the NG5 and cc5 references, correspondingly. These deviations are sligthly smaller in the Org24 test set comparing to the G2/97 test set. The diminished deviations are likely due to the larger atom count in the Org24 test set. Further limitation of the NAOs’ spatial extent in the DZP[50 meV] basis set causes an MSE of 0.22 eV with respect to the DZP[1 meV] basis set for IPs and 0.32 eV for EAs. These deviations became slightly larger for the Org24 test set in comparison to the G2/97 test set results, although they remain smaller than the MSEs between NG5 and DZP[1 meV] basis sets. In panel b of Figures 4 and 5, we show heatmaps of all MAEs and MxADs among cc5, cc5l2, NG5, DZP[1 meV], and DZP[50 meV] basis sets. The MAEs are the same as MSEs except for the pair cc5l2 versus NG5. The same magnitude of the MSEs and MAEs is a sign of the errors’ systematicity. This systematicity is further illustrated with histograms and the calculated error spreads in our Supporting Information. MxADs are inevitably larger than MAEs; however, they became smaller compared to the G2/97 test set. Furthermore, it is worth noting that the deviations for the IPs and EAs occur in the same direction. Namely, worsening of the basis set causes the underestimation of the IPs and EAs by approximately the same amounts in the Org24 test set. Summarizing subsections 4.1 and 4.2, we note that GW calculations with the PP description of core electrons deviate by no more than 0.2 eV from their AE counterparts computed using basis sets of similar quality. Moreover, relatively small (double-ζ polarized) and short-ranged (energy shift 50 meV) basis sets of NAOs provide an acceptable discretization of spatial degrees of freedom for GW calculations. Namely, the computationally affordable DZP[50 meV] basis set causes a rather systematic, rigid increase of quasiparticle energy levels by about 0.9 eV with respect to the computationally expensive

another 0.16 eV deviation with respect to the DZP[1 meV] basis. On panel b of Figure 3, we show a heatmap of all MAEs and MxADs among cc5, cc5l2, NG5, DZP[1 meV], and DZP[50 meV] basis sets. The MAEs shown in the lower left corner of the heatmap are the same as MSEs except for the pair cc5l2 versus NG5 (MAE = 0.08 eV, MSE = 0.05 eV). The same magnitude of the MSEs and MAEs is a sign of the errors’ systematicity. This systematicity is further illustrated with histograms and the calculated error spreads in our Supporting Information as well as on the example of the Org24 test set, in the next subsection 4.2. MxADs are shown in the upper right corner of the heatmap. MxAEs are inevitably larger than MAEs. As a rule of thumb, MxAD are larger than MAE by factor 2, but not exceeding 1.53 eV. 4.2. Larger Acceptor Molecules. In Figures 4 and 5, we show the deviations in calculated IPs and EAs, correspondingly, for the Org24 test set. Figures 4 and 5 are similar to Figure 3 showing the deviations for the G2/97 test set. Similarly to the G2/97 test set, the omission of higher angular momenta leads to underestimation of both IP and EA, as can be read from columns cc5l2. In contrast to the small molecules of the G2/97 test set, the omission of higher angular momenta causes smaller deviations in the Org24 test set. In the case of the G2/97 test set, the MSE is 0.31, while in the case of the Org24 test set, the MSE amounts to 0.28 for IPs and 0.27 eV for EAs. We attribute this reduction in the error to the larger size of the molecules in the Org24 test set. The mean number of atoms in the G2/97 test set is only ∼6 while it is ∼16 in the Org24 test set. The basis functions centered on neighboring atoms overlap significantly and can provide a partial replacement for the missing high-order angular momentum functions. Panels a in Figures 4 and 5 show that, similarly to the G2/97 test set, the MSE between cc5l2 and NG5 are the smallest in the set of discrepancies and amount only 0.05 eV for both IPs and EAs of the Org24 test set. The corresponding MAEs are smaller for the Org24 test set (0.05 and 0.04 eV for IP and EA, correspondingly) than for the G2/97 test set. Similarly, panels b in Figures 4 and 5 corroborate that the MxAD between cc5l2 and NG5 are twice smaller for the Org24 test set comparing with the G2/97 test set and does not exceed 0.14 eV. These small discrepancies confirm the average good performance of norm-conserving GGA PPs for GW calculations, in line with I


Article


Figure 6. Ball and stick model of the fluorinated thiazole diarylethene derivative (11aS,11bS)-4,4,5,5,6,6-hexafluoro-11b-methyl-2-phenyl-11apropyl-5,6,11a,11b-tetrahydro-4H-benzofuro[2′,3′:6,7]indeno[4,5-d]thiazole and its structural formula. We denote the diarylethene derivative as DAE-iii. DAE-iii is photochromic and we model its closed isomer.

Figure 7. Time evolution of the seven highest occupied levels and the five lowest unoccupied levels of the DAE-iii compound for the thermostat temperature T0 = 300 K. The eigenvalues are computed with DFT (PBE functional) and with the HF+G0W0 approximation on panels a and b, correspondingly. Note the different energy scale in the two panels. Light blue horizontal bars indicate the average positions and widths of the respective energy level sequences. The 800 snapshots used in the AIMD average of the DoS are shown. The HOMO−1, HOMO and LUMO, LUMO+1, LUMO+2, and LUMO+3 levels are well separated from the rest of the energy levels for both approximations.

isomers. Born−Oppenheimer dynamics is suitable to generate the relevant trajectories because of the absence of thermochromic effects. We selected a fluorinated thiazole DAE derivative that is shown in Figure 6. This DAE derivative is identical to compound iii in the comprehensive review in ref 146. Therefore, we refer to this compound as DAE-iii in this section. Moreover, we model only the closed (more stable) isomer of DAE-iii in order to reduce the computational demand, as in our previous simulations41 of the optical polarizability at the TDDFT level. We compare two AIMDs characterized by the Nosé thermostat temperatures of T0 = 100 and 300 K. Prior to these AIMDs, we calculated the phonon spectrum of DAE-iii using the VIBRA utility of the SIESTA software suite. The

cc-pV5Z basis set of functions with much larger spatial extent and the angular momentum content. In the following, we will employ the DZP[50 meV] basis in the AIMD+G 0 W 0 simulation to study the thermal agitation effects in the photoelectron spectra of a relatively large diarylethene compound. 4.3. AIMD+GW Simulations of a Diarylethene Derivative. Diarylethenes (DAEs) form a promising family of photochromic compounds that has been studied extensively.95,146 Among a variety of spectroscopical measurements, photoelectron spectroscopy has been used to measure the absolute positions of the electronic levels of the individual DAE isomers.147−151 There is no thermochromic effect in DAEs; i.e., thermal motion causes only adiabatic level broadening,146,151 without switching between open and closed J


Article


Figure 8. AIMD averages of the DoS for the thermostat temperatures T0 = 100 and 300 K computed with DFT (PBE) and HF+G0W0. The red solid lines represent an average over 800 snapshots of the Nosé−Hoover AIMD. The blue dashed lines correspond to the equilibrium configuration with a Gaussian broadening matching the widths of the HOMO thermal broadening of the respective level of theory.

threshold λPB = 10−7. Such PB contains 78 functions for each of the “heavy” atoms (C, N, O, F, S) and 15 functions for the hydrogen atom. Similar to previous calculations, G0 W0 corrections were determined for the 10 highest occupied and the 10 lowest virtual levels. These G0W0-corrected levels were converged within 10−6 Ha in less than 12 iterations for all the 1600 GW calculations. It is instructive to see the HF+G0W0-corrected levels as a function of time. In Figure 7, we show the evolution of the seven highest occupied and five lowest unoccupied levels at the thermostat temperature of T0 = 300 K. The light blue horizontal bars are centered at the mean energy of the respective time sequence ⟨E⟩, and the width of the bars is equal to the respective FWHM. We see that two of the highest occupied and four of the lowest unoccupied levels evolve within nonoverlapping ranges. In panel a we see that the horizontal bars of HOMO−5, HOMO−4, and HOMO−3 overlap, implying the possibility of level crossings (mixing) for these levels in AIMD+GGA simulation. As in AIMD+GGA simulation, the HF+GW energy bands of HOMO−4, HOMO−3, and HOMO−2 do overlap, as can be seen on panel b of Figure 7. The overlapping energy ranges hamper the determination of the levels’ ensemble averages, i.e., their widths and positions, due to eventual loss of the time-continuity of the corresponding mean-field orbitals. In contrast to the individual level widths, the DoS is a cumulative property of all energy levels regardless of the symmetry of the underlying MOs. Therefore, the ensemble average of the DoS is easier to estimate in a wide energy range. In Figure 8, we compare the ensemble-averaged DoSs for PBE (panels a and b) and HF+G0W0 (panels c and d) levels of theory, for both thermostat temperatures T0 = 100 and 300 K. The ensemble-averaged DoS is computed according to eq 18. Because the focus is on the thermal level broadening, we assume a small Lorentzian broadening constant ε = 0.026 eV

minimal and maximal frequencies of normal modes were found to be 44 and 3135 cm−1, which translate to time periods 758 and 10.6 fs, correspondingly. Guided by these time periods, we carried out the AIMD simulations with the time step of 1 fs for a total of 20 000 steps. Each HF+G0W0 calculation is considerably more expensive than a DFT force calculation (10 h versus 3 min). Due to the elevated computational cost of the HF+G0W0, we have to downsample the initial AIMD ensemble. Using again the phonon frequencies as a guide, we decided for a sampling time step of 5 fs and a total of S = 800 snapshots in the DoS average (18), covering 4 ps of the entire AIMD propagation. Because, we accumulated molecular geometries for a total of 20 ps, we have an additional degree of freedom in choosing the timewindow for the AIMD+GW ensemble. Using this freedom, we minimized the difference between MD averages of the instantaneous temperature ⟨T⟩ and temperature fluctuations ⟨ΔT2⟩ and the corresponding ensemble averages for the canonical ensemble.130 For instance, for the thermostat temperature T0 = 100 K, we found the MD averages ⟨T⟩ = 105 K and ⟨ΔT2⟩ = 1.05 T02/CV, where CV = 1.5N (where N = 54 is the number of atoms), computing averages between the 4000th and 8000th AIMD steps. For the thermostat temperature T0 = 300 K, we found ⟨T⟩ = 290 K and ⟨ΔT2⟩ = 1.08 T02/CV, computing averages between the 10000th and 14000th AIMD steps. The parameters of the HF+G0W0 calculation were chosen, lowering the tolerance criteria as compared to the calculations in sections 4.1 and 4.2. In particular, we used only the pure NAO, split DZP[50 meV] basis set, which leads to 10 times faster calculations than the mixed NAO-Gaussian-double-ζ (NG2) basis set. The frequency grid along the imaginary axis contained 48 points, while the frequency tolerance ϵtol = 10−10 is kept equal to the previous choice. Moreover, we used a smaller PB set characterized by a 10 times larger eigenvalue K


Article


Table 1. MD-Averaged Energy Levels (eV) of the Two Highest Occupied and the Four Lowest Unoccupied States for Three Electronic Structure Methods: DFT (PBE Functional), HF, and HF+G0W0a GGA LUMO+3 LUMO+2 LUMO+1 LUMO HOMO HOMO−1

HF+G0W0

HF

0K

100 K

300 K

0K

100 K

300 K

0K

100 K

300 K

−1.5571 −1.8476 −2.5472 −3.7186 −5.1434 −6.1157

−1.6120 −1.8558 −2.6151 −3.7321 −5.1363 −6.1412

−1.6599 −1.9082 −2.6220 −3.7507 −5.1122 −6.1190

3.2104 2.7027 1.7618 0.1911 −7.4654 −9.1348

3.1464 2.7132 1.6671 0.1657 −7.4712 −9.1619

3.0822 2.6627 1.6681 0.1614 −7.4582 −9.1436

1.6868 1.3304 0.3621 −1.1823 −7.3004 −8.6749

1.6213 1.3389 0.2746 −1.2011 −7.2999 −8.7138

1.5630 1.2812 0.2727 −1.2161 −7.2840 −8.7030

a

In the column 0 K, the equilibrium geometry values are given. In the columns labeled 100 and 300 K the AIMD averages for the corresponding thermostat temperature are given.

Table 2. Diagonal Matrix Elements of the Electron−Phonon Coupling Vep (eV), for the Two Highest Occupied and Four Lowest Unoccupied Levelsa GGA LUMO+3 LUMO+2 LUMO+1 LUMO HOMO HOMO−1 deg set

HF+G0W0

HF

0K

100 K

300 K

0K

100 K

300 K

0K

100 K

300 K

0.25 0.35 0.20 0.20 0.28 0.21 0.13

0.28 0.54 0.34 0.20 0.27 0.28 0.12

0.23 0.43 0.26 0.20 0.32 0.21 0.13

0.45 0.59 0.38 0.63 1.21 0.67 0.18

0.49 0.92 0.64 0.56 0.92 0.91 0.18

0.38 0.76 0.48 0.63 1.39 0.60 0.16

0.37 0.55 0.34 0.52 0.83 0.51 0.093

0.44 0.81 0.54 0.49 0.65 0.67 0.11

0.32 0.59 0.43 0.49 0.95 0.53 0.11

The derivatives of the energy levels αm with respect to normal modes Qm are used in the columns 0 K according to eq 24, while the standard deviations of the energy levels extracted from AIMD are used in the colums labeled by 100 and 300 K according to eq 25. In the last row, we state Vep averaged over sets of nearly degenerate states. For the GGA, the set includes HOMO−3, HOMO−4, and HOMO−5, for HF the set includes HOMO−2 and HOMO−3 and for HF+GW the set includes HOMO−2, HOMO−3, and HOMO−4 levels. See the description of the table in the text for the motivation of the sets’ choice. a

(25) show a stronger dependence on the method of computing the electronic excitation energies. In Table 2, we gathered the coupling constants Vep, computed with the frozen-phonon method, according to eq 24, and acquired from the AIMD standard deviations, according to eq 25. In the case of nearly degenerate energy levels, the AIMD and frozen-phonon calculations of the electron−phonon coupling coefficients (24) and (25) can be reconciled by averaging over the set of nearly degenerate energy levels. The averaging should be done for both the AIMD and the frozen-phonon approach by defining an average energy of each snapshot E̅ s = ⟨Es⟩set and average the normal-mode coupling constants α̅ m = ⟨αm⟩set. In both cases, the average is done over the set of nearly degenerate levels. In our Supporting Information, the time evolution of the energy levels is shown in more detail. Inspecting Table 2, we see that the coupling constants estimated by the HF+G0W0 method are a factor of 2−3 larger than for the PBE estimate, while HF always gives the largest coupling. Similar trends were found in the widely studied electron−phonon coupling in C60 fullerene.13,70,72,73 The good agreement, particularly for PBE, between the couplings computed at 0 K within the harmonic approximation and those estimated from the AIMD results justifies a posteriori the use of the harmonic approximation and the linear electron− phonon coupling (20). The small dependence of the mean level position on temperature gives a measure of the deviations from such linear/harmonic behavior. In some cases, the AIMD estimates of the electron−phonon coupling constants Vep do change noticeably, between the thermostats at 100 and 300 K for HF and HF+G0W0 methods and to a lesser extend for PBE. The variation of the AIMD

for the spectra of each individual configuration at both thermostat temperatures. In addition to the ensemble-averaged DoS, we also show the equilibrium DoS plotted using a normal Gaussian distribution (22). The standard deviation σ is chosen to match the FWHM of the temperature broadened HOMO level for each electronic structure method (DFT or HF +G0W0). Inspecting Figure 8, we see that the level width is larger for the higher thermostat temperature (panels b and d) than for the lower thermostat temperature (panels a and c). Within the harmonic approximation the level widths are expected to be proportional to the square root of the thermostat temperature T0 , as seen in eq 23. The accuracy of the square root proportionality is apparent from the estimation of the electron−phonon coupling Vep, as discussed below (discussion of Table 2). In the vicinity of the HOMO and the LUMO, the ensemble-averaged DoS matches the equilibrium DoS, while there are noticeable differences for the lower and higher binding energies. Looking closely at the ensemble-averaged DoS, we see that the mean level positions have a tendency to move to higher energies for the occupied states and lower for the unoccupied states as the temperature is increased. To quantify this faint trend, we gathered the MD averages of the levels’ energies N ⟨En⟩ = Ns−1 ∑s =s 1 Es , n in Table 1. Inspecting Table 1, we see that the HOMO−LUMO gap tends to shrink with increasing temperature for all electronic structure methods, albeit at a rather slow pace (1−2) × 10−4 eV/K. In contrast to the MD averages of the energy levels, the electron−phonon coupling constants Vep defined in eqs 24 and L


Article

Journal of Chemical Theory and Computation estimates of the coupling constants Vep with temperature (less than 50%) are not consistent with the harmonic ansatz (19) and the linearity hypothesis (20). The inconsistency may be due to breaking of the linearity hypothesis (20) because the sampling of trajectories is the same for all the electronic structure methods (DFT, HF, and HF+G0W0), and thus possible anharmonic effects should affect all of them equally. In any case, the studied AIMD+GW simulations point to the importance of going beyond standard semilocal GGA for the accurate computation of electron−phonon couplings and the quantitative thermal averaging of excitation spectra.

we envision that further algorithmic developments of the presented methods will result in a computationally superior approach in the near future. It is very likely that the use of iterative techniques to treat the screened interaction leads to a considerable speed-up of one-shot GW calculations. Establishing a more direct communication in-between subsequent AIMD steps may bring about further improvement.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.9b00436.

5. CONCLUSIONS AND OUTLOOK In this paper, we presented an implementation of Hedin’s GW approximation that is suitable for short-ranged, numerical atomic orbitals (NAO) in electronic structure calculations employing pseudopotentials. The implementation makes use of a contour deformation technique for the calculation of the self-energy and numerical atom-centered basis sets to expand orbital products. The short-ranged, NAOs are designed for fast ground-state calculations, for example, making it possible to perform long large-scale AIMD simulations. In this paper, we explored the suitability of a small double-ζ polarized (DZP) basis set formed by short-ranged NAOs for GW calculations. We did so by comparing with all-electron GW calculations using a rather complete quintuple ζ basis. The short-ranged DZP NAO bases cause a systematic upshift of HOMO and LUMO levels with a magnitude of ∼0.5 eV. These deviations are related to both the low number of orbitals and the short spatial support of these orbitals. However, NAOs are computationally rather efficient, and it is possible to envision that the deviations due to the finite support of the orbitals, when affecting similarly the HOMO and LUMO levels, could be alleviated by a molecularsize-dependent a posteriori correction.152 More accurate results are obtained using combined basis sets composed of NAOs and Gaussian functions. In comparison to the effect of the basis set, the uncertainties due to the use of GGA pseudopotentials are much smaller. We find deviations in the quasiparticle spectra of less than 0.2 eV with respect to the corresponding all-electron references, i.e., comparing against the all-electron basis sets with limited angular momentum content. We demonstrate the efficiency of our approach in combined AIMD+GW simulations. In particular, we study classical finitetemperature effects on the line shapes in photoemission spectra. Increasing the temperature in our AIMD simulations gives rise to a larger broadening of the photoemission lines. A direct analysis of this thermal broadening allows us to investigate the dependence of the diagonal electron−phonon coupling on the employed electronic structure method. We find that the electron−phonon coupling of the frontier orbitals is roughly twice as large in the GWA as compared to semilocal GGA-DFT. Moreover, a comparison of GGA-DFT and GW does not reveal a systematic trend across the quasiparticle levels; i.e., the ratio of the electron phonon couplings Vep GWA/ Vep GGA varies for different quasiparticle states. This emphasizes the importance of going beyond simple DFT functionals, or using correlated methods, for the calculation of electron− phonon couplings and thermal-averaged excitation spectra. The treatment of excited states at the level of GW is computationally still about 2 orders of magnitude more expensive than DFT with pure density functionals. However,

■

Details of the constant temperature molecular dynamics simulations, derivation of the electron−phonon coupling coefficients, details of the frozen-phonon approximation, ionization potentials and electron affinities of the G2/97 and Org24 test sets (including figures of the time evolution of highest and lowest occupied levels, ball and stick model, HOMO and LUMO isosurfaces, and error distribution graphs; tables of frequencies and coupling coefficients, numbers of basis functions per atom, CAS Registry Numbers, IPs, and EAs) (PDF)

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +34 943 01 84 15. Fax: +34 943 01 56 00. ORCID

Peter Koval: 0000-0002-5461-2278 Moritz Müller: 0000-0003-0111-8214 Daniel Sánchez-Portal: 0000-0001-6860-8790 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Prof. Giulia Galli and Dr. T. Anh Pham provided results of their GW calculations in a convenient text form. Prof. Diego Peña Gil discerned the chemical designation of the diarylethene photochromic compound we modeled. P.K. and D.S.P. acknowledge support from Spanish MINECO Grants MAT2016-78293-C6-4-R and RTC-2016-5681-7 (SIESTAPRO). P.K. acknowledges financial support from the Fellows Gipuzkoa program of the Gipuzkoako Foru Aldundia through the FEDER funding scheme of the European Union. M.M. and M.p.L. acknowledge support from the Department of education of the Basque Government, as well as from Euskampus and Donostia International Physics Center.

■

REFERENCES

(1) Ma, Y.; Rohlfing, M.; Molteni, C. Excited states of biological chromophores studied using many-body perturbation theory: Effects of resonant-antiresonant coupling and dynamical screening. Phys. Rev. B: Condens. Matter Mater. Phys. 2009, 80, 241405. (2) Faber, C.; Attaccalite, C.; Olevano, V.; Runge, E.; Blase, X. Firstprinciples GW calculations for DNA and RNA nucleobases. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 83, 115123. (3) van Setten, M. J.; Caruso, F.; Sharifzadeh, S.; Ren, X.; Scheffler, M.; Liu, F.; Lischner, J.; Lin, L.; Deslippe, J. R.; Louie, S. G.; Yang, C.; Weigend, F.; Neaton, J. B.; Evers, F.; Rinke, P. GW 100: Benchmarking G0W0 for Molecular Systems. J. Chem. Theory Comput. 2015, 11, 5665−5687.

M


Article


energy calculations within the GW approximation. J. Comput. Phys. 2015, 286, 1−13. (22) Poncé, S.; Antonius, G.; Boulanger, P.; Cannuccia, E.; Marini, A.; Côté, M.; Gonze, X. Verification of first-principles codes: Comparison of total energies, phonon frequencies, electron−phonon coupling and zero-point motion correction to the gap between ABINIT and QE/Yambo. Comput. Mater. Sci. 2014, 83, 341−348. (23) Wilhelm, J.; Golze, D.; Talirz, L.; Hutter, J.; Pignedoli, C. A. Toward GW Calculations on Thousands of Atoms. J. Phys. Chem. Lett. 2018, 9, 306−312. (24) Giannozzi, P.; Andreussi, O.; Brumme, T.; Bunau, O.; Nardelli, M. B.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Cococcioni, M.; Colonna, N.; Carnimeo, I.; Corso, A. D.; de Gironcoli, S.; Delugas, P.; Jr, R. A. D.; Ferretti, A.; Floris, A.; Fratesi, G.; Fugallo, G.; Gebauer, R.; Gerstmann, U.; Giustino, F.; Gorni, T.; Jia, J.; Kawamura, M.; Ko, H.-Y.; Kokalj, A.; Kücu̧ ̈kbenli, E.; Lazzeri, M.; Marsili, M.; Marzari, N.; Mauri, F.; Nguyen, N. L.; Nguyen, H.-V.; Otero-de-la Roza, A.; Paulatto, L.; Poncé, S.; Rocca, D.; Sabatini, R.; Santra, B.; Schlipf, M.; Seitsonen, A. P.; Smogunov, A.; Timrov, I.; Thonhauser, T.; Umari, P.; Vast, N.; Wu, X.; Baroni, S. Advanced capabilities for materials modelling with Quantum ESPRESSO. J. Phys.: Condens. Matter 2017, 29, 465901. (25) Liu, P.; Kaltak, M.; Klimeš, J.; Kresse, G. Cubic scaling GW: Towards fast quasiparticle calculations. Phys. Rev. B: Condens. Matter Mater. Phys. 2016, 94, 165109. (26) Govoni, M.; Galli, G. Large Scale GW Calculations. J. Chem. Theory Comput. 2015, 11, 2680−2696. (27) Bruneval, F.; Rangel, T.; Hamed, S. M.; Shao, M.; Yang, C.; Neaton, J. B. MOLGW 1: Many-body perturbation theory software for atoms, molecules, and clusters. Comput. Phys. Commun. 2016, 208, 149−161. (28) Ren, X.; Rinke, P.; Blum, V.; Wieferink, J.; Tkatchenko, A.; Sanfilippo, A.; Reuter, K.; Scheffler, M. Resolution-of-identity approach to Hartree−Fock, hybrid density functionals, RPA, MP2 and GW with numeric atom-centered orbital basis functions. New J. Phys. 2012, 14, 053020. (29) van Setten, M. J.; Weigend, F.; Evers, F. The GW-Method for Quantum Chemistry Applications: Theory and Implementation. J. Chem. Theory Comput. 2013, 9, 232−246. (30) Gui, X.; Holzer, C.; Klopper, W. Accuracy Assessment of GW Starting Points for Calculating Molecular Excitation Energies Using the Bethe−Salpeter Formalism. J. Chem. Theory Comput. 2018, 14, 2127−2136. (31) Rostgaard, C.; Jacobsen, K. W.; Thygesen, K. S. Fully selfconsistent GW calculations for molecules. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 81, 085103. (32) Sun, Q.; Berkelbach, T. C.; Blunt, N. S.; Booth, G. H.; Guo, S.; Li, Z.; Liu, J.; McClain, J. D.; Sayfutyarova, E. R.; Sharma, S.; Wouters, S.; Chan, G. K.-L. PySCF: the Python-based simulations of chemistry framework. Wiley Interdisciplinary Reviews: Computational Molecular Science 2018, 8, No. e1340. (33) Nabok, D.; Gulans, A.; Draxl, C. Accurate all-electron G0W0 quasiparticle energies employing the full-potential augmented planewave method. Phys. Rev. B: Condens. Matter Mater. Phys. 2016, 94, 035118. (34) Blase, X.; Duchemin, I.; Jacquemin, D. The Bethe−Salpeter equation in chemistry: relations with TD-DFT, applications and challenges. Chem. Soc. Rev. 2018, 47, 1022−1043. (35) Rohlfing, M.; Louie, S. G. Electron-hole excitations and optical spectra from first principles. Phys. Rev. B: Condens. Matter Mater. Phys. 2000, 62, 4927−4944. (36) Kotani, T.; van Schilfgaarde, M. Fusion of the LAPW and LMTO methods: The augmented plane wave plus muffin-tin orbital method. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 81, 125117. (37) Kutepov, A. L.; Oudovenko, V. S.; Kotliar, G. Linearized selfconsistent quasiparticle GW method: Application to semiconductors and simple metals. Comput. Phys. Commun. 2017, 219, 407−414.

(4) Knight, J. W.; Wang, X.; Gallandi, L.; Dolgounitcheva, O.; Ren, X.; Ortiz, J. V.; Rinke, P.; Körzdörfer, T.; Marom, N. Accurate Ionization Potentials and Electron Affinities of Acceptor Molecules III: A Benchmark of GW Methods. J. Chem. Theory Comput. 2016, 12, 615−626. (5) Vlček, V.; Rabani, E.; Neuhauser, D.; Baer, R. Stochastic GW Calculations for Molecules. J. Chem. Theory Comput. 2017, 13, 4997− 5003. (6) Rohlfing, M.; Krüger, P.; Pollmann, J. Efficient scheme for GW quasiparticle band-structure calculations with applications to bulk Si and to the Si(001)-(2 × 1) surface. Phys. Rev. B: Condens. Matter Mater. Phys. 1995, 52, 1905−1917. (7) Ku, W.; Eguiluz, A. G. Band-Gap Problem in Semiconductors Revisited: Effects of Core States and Many-Body Self-Consistency. Phys. Rev. Lett. 2002, 89, 126401. (8) Liao, P.; Carter, E. A. Testing variations of the GW approximation on strongly correlated transition metal oxides: hematite (α-Fe2O3) as a benchmark. Phys. Chem. Chem. Phys. 2011, 13, 15189−15199. (9) Körzdörfer, T.; Parrish, R. M.; Marom, N.; Sears, J. S.; Sherrill, C. D.; Brédas, J.-L. Assessment of the performance of tuned rangeseparated hybrid density functionals in predicting accurate quasiparticle spectra. Phys. Rev. B: Condens. Matter Mater. Phys. 2012, 86, 205110. (10) Gulans, A.; Kontur, S.; Meisenbichler, C.; Nabok, D.; Pavone, P.; Rigamonti, S.; Sagmeister, S.; Werner, U.; Draxl, C. exciting: a fullpotential all-electron package implementing density-functional theory and many-body perturbation theory. J. Phys.: Condens. Matter 2014, 26, 363202. (11) Lazzeri, M.; Attaccalite, C.; Wirtz, L.; Mauri, F. Impact of the electron-electron correlation on phonon dispersion: Failure of LDA and GGA DFT functionals in graphene and graphite. Phys. Rev. B: Condens. Matter Mater. Phys. 2008, 78, 081406. (12) Grüneis, A.; Attaccalite, C.; Rubio, A.; Vyalikh, D. V.; Molodtsov, S. L.; Fink, J.; Follath, R.; Eberhardt, W.; Büchner, B.; Pichler, T. Electronic structure and electron-phonon coupling of doped graphene layers in KC8. Phys. Rev. B: Condens. Matter Mater. Phys. 2009, 79, 205106. (13) Faber, C.; Janssen, J. L.; Côté, M.; Runge, E.; Blase, X. Electron-phonon coupling in the C60 fullerene within the many-body GW approach. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 84, 155104. (14) Noffsinger, J.; Kioupakis, E.; Van de Walle, C. G.; Louie, S. G.; Cohen, M. L. Phonon-Assisted Optical Absorption in Silicon from First Principles. Phys. Rev. Lett. 2012, 108, 167402. (15) Wright, A. D.; Verdi, C.; Milot, R. L.; Eperon, G. E.; PérezOsorio, M. A.; Snaith, H. J.; Giustino, F.; Johnston, M. B.; Herz, L. M. Electron-phonon coupling in hybrid lead halide perovskites. Nat. Commun. 2016, 7, 11755. (16) Shigeta, Y.; Nagao, H.; Nishikawa, K.; Yamaguchi, K. Density functional theory without the Born−Oppenheimer approximation. II. Green function techniques. Int. J. Quantum Chem. 1999, 75, 875−883. (17) Prendergast, D.; Galli, G. X-Ray Absorption Spectra of Water from First Principles Calculations. Phys. Rev. Lett. 2006, 96, 215502. (18) Swartz, C. W.; Wu, X. Ab initio Studies of Ionization Potentials of Hydrated Hydroxide and Hydronium. Phys. Rev. Lett. 2013, 111, 087801. (19) Drisdell, W. S.; Poloni, R.; McDonald, T. M.; Pascal, T. A.; Wan, L. F.; Pemmaraju, C. D.; Vlaisavljevich, B.; Odoh, S. O.; Neaton, J. B.; Long, J. R.; Prendergast, D.; Kortright, J. B. Probing the mechanism of CO2 capture in diamine-appended metal-organic frameworks using measured and simulated X-Ray spectroscopy. Phys. Chem. Chem. Phys. 2015, 17, 21448−21457. (20) Pham, T. N.; Ono, S.; Ohno, K. Ab initio molecular dynamics simulation study of successive hydrogenation reactions of carbon monoxide producing methanol. J. Chem. Phys. 2016, 144, 144309. (21) Liu, F.; Lin, L.; Vigil-Fowler, D.; Lischner, J.; Kemper, A. F.; Sharifzadeh, S.; da Jornada, F. H.; Deslippe, J.; Yang, C.; Neaton, J. B.; Louie, S. G. Numerical integration for ab initio many-electron self N


Article

Journal of Chemical Theory and Computation (38) Nilsson, F.; Aryasetiawan, F. Recent Progress in First-Principles Methods for Computing the Electronic Structure of Correlated Materials. Computation 2018, 6, 26. (39) Caruso, F.; Rinke, P.; Ren, X.; Scheffler, M.; Rubio, A. Unified description of ground and excited states of finite systems: The selfconsistent GW approach. Phys. Rev. B: Condens. Matter Mater. Phys. 2012, 86, 081102. (40) Koval, P.; Foerster, D.; Sánchez-Portal, D. Fully self-consistent GW and quasiparticle self-consistent GW for molecules. Phys. Rev. B: Condens. Matter Mater. Phys. 2014, 89, 155417. (41) Koval, P.; Barbry, M.; Sánchez-Portal, D. PySCF-NAO: An efficient and flexible implementation of linear response timedependent density functional theory with numerical atomic orbitals. Comput. Phys. Commun. 2019, 236, 188−204. (42) Rossi, T. P.; Zugarramurdi, A.; Puska, M. J.; Nieminen, R. M. Quantized Evolution of the Plasmonic Response in a Stretched Nanorod. Phys. Rev. Lett. 2015, 115, 236804. (43) Marchesin, F.; Koval, P.; Barbry, M.; Aizpurua, J.; SánchezPortal, D. Plasmonic Response of Metallic Nanojunctions Driven by Single Atom Motion: Quantum Transport Revealed in Optics. ACS Photonics 2016, 3, 269−277. (44) Foerster, D.; Koval, P.; Sánchez-Portal, D. An O(N3) implementation of Hedin’s GW approximation for molecules. J. Chem. Phys. 2011, 135, 074105. (45) Hüser, F.; Olsen, T.; Thygesen, K. S. Quasiparticle GW calculations for solids, molecules, and two-dimensional materials. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 87, 235132. (46) Sankey, O. F.; Niklewski, D. J. Ab initio multicenter tightbinding model for molecular-dynamics simulations and other applications in covalent systems. Phys. Rev. B: Condens. Matter Mater. Phys. 1989, 40, 3979−3995. (47) Junquera, J.; Paz, O.; Sánchez-Portal, D.; Artacho, E. Numerical atomic orbitals for linear-scaling calculations. Phys. Rev. B: Condens. Matter Mater. Phys. 2001, 64, 235111. (48) Rossi, T. P.; Lehtola, S.; Sakko, A.; Puska, M. J.; Nieminen, R. M. Nanoplasmonics simulations at the basis set limit through completeness-optimized, local numerical basis sets. J. Chem. Phys. 2015, 142, 094114. (49) Chen, Y.-C.; Chen, J.-Z.; Michaud-Rioux, V.; Shi, Q.; Guo, H. Efficient evaluation of nonlocal operators in density functional theory. Phys. Rev. B: Condens. Matter Mater. Phys. 2018, 97, 075139. (50) Soler, J. M.; Artacho, E.; Gale, J. D.; García, A.; Junquera, J.; Ordejón, P.; Sánchez-Portal, D. The SIESTA method for ab initio order-N materials simulation. J. Phys.: Condens. Matter 2002, 14, 2745−2779. (51) Artacho, E.; Anglada, E.; Diéguez, O.; Gale, J. D.; García, A.; Junquera, J.; Martin, R. M.; Ordejón, P.; Pruneda, J. M.; SánchezPortal, D.; Soler, J. M. The SIESTA method; developments and applicability. J. Phys.: Condens. Matter 2008, 20, 064208. (52) Ozaki, T.; Kino, H. Efficient projector expansion for the ab initio LCAO method. Phys. Rev. B: Condens. Matter Mater. Phys. 2005, 72, 045121. (53) Mendieta-Moreno, J. I.; Walker, R. C.; Lewis, J. P.; GómezPuertas, P.; Mendieta, J.; Ortega, J. FIREBALL/AMBER: An Efficient Local-Orbital DFT QM/MM Method for Biomolecular Systems. J. Chem. Theory Comput. 2014, 10, 2185−2193. (54) Walter, M.; Häkkinen, H.; Lehtovaara, L.; Puska, M.; Enkovaara, J.; Rostgaard, C.; Mortensen, J. J. Time-dependent density-functional theory in the projector augmented-wave method. J. Chem. Phys. 2008, 128, 244101. (55) Kenny, S.; Horsfield, A. Plato: A localised orbital based density functional theory code. Comput. Phys. Commun. 2009, 180, 2616− 2621. (56) Troullier, N.; Martins, J. L. Efficient pseudopotentials for planewave calculations. Phys. Rev. B: Condens. Matter Mater. Phys. 1991, 43, 1993−2006. (57) García, A.; Verstraete, M. J.; Pouillon, Y.; Junquera, J. The psml format and library for norm-conserving pseudopotential data curation and interoperability. Comput. Phys. Commun. 2018, 227, 51−71.

(58) Cole, D. J.; Hine, N. D. M. Applications of large-scale density functional theory in biology. J. Phys.: Condens. Matter 2016, 28, 393001. (59) Ratcliff, L. E.; Mohr, S.; Huhs, G.; Deutsch, T.; Masella, M.; Genovese, L. Challenges in large scale quantum mechanical calculations. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2017, 7, No. e1290. (60) Sánchez-Portal, D.; Artacho, E.; Soler, J. M. Analysis of atomic orbital basis sets from the projection of plane-wave results. J. Phys.: Condens. Matter 1996, 8, 3859. (61) Lee, K.; Yu, J.; Morikawa, Y. Comparison of localized basis and plane-wave basis for density-functional calculations of organic molecules on metals. Phys. Rev. B: Condens. Matter Mater. Phys. 2007, 75, 045402. (62) Louwerse, M. J.; Rothenberg, G. Transferable basis sets of numerical atomic orbitals. Phys. Rev. B: Condens. Matter Mater. Phys. 2012, 85, 035108. (63) Smidstrup, S.; Stradi, D.; Wellendorff, J.; Khomyakov, P. A.; Vej-Hansen, U. G.; Lee, M.-E.; Ghosh, T.; Jónsson, E.; Jónsson, H.; Stokbro, K. First-principles Green’s-function method for surface calculations: A pseudopotential localized basis set approach. Phys. Rev. B: Condens. Matter Mater. Phys. 2017, 96, 195309. (64) Silva, F. W. N.; Cruz-Silva, E.; Terrones, M.; Terrones, H.; Barros, E. B. BNC nanoshells: a novel structure for atomic storage. Nanotechnology 2017, 28, 465201. (65) Trail, J. R.; Needs, R. J. Norm-conserving Hartree−Fock pseudopotentials and their asymptotic behavior. J. Chem. Phys. 2005, 122, 014112. (66) Shishkin, M.; Kresse, G. Implementation and performance of the frequency-dependent GW method within the PAW framework. Phys. Rev. B: Condens. Matter Mater. Phys. 2006, 74, 035101. (67) Körbel, S.; Boulanger, P.; Duchemin, I.; Blase, X.; Marques, M. A. L.; Botti, S. Benchmark Many-Body GW and Bethe−Salpeter Calculations for Small Transition Metal Molecules. J. Chem. Theory Comput. 2014, 10, 3934−3943. (68) Brunk, E.; Rothlisberger, U. Mixed Quantum Mechanical/ Molecular Mechanical Molecular Dynamics Simulations of Biological Systems in Ground and Electronically Excited States. Chem. Rev. 2015, 115, 6217−6263. (69) Kresse, G.; Hafner, J. Ab initio molecular dynamics for openshell transition metals. Phys. Rev. B: Condens. Matter Mater. Phys. 1993, 48, 13115−13118. (70) Antropov, V. P.; Gunnarsson, O.; Liechtenstein, A. I. Phonons, electron−phonon, and electron−plasmon coupling in C60 compounds. Phys. Rev. B: Condens. Matter Mater. Phys. 1993, 48, 7651− 7664. (71) Gunnarsson, O.; Handschuh, H.; Bechthold, P. S.; Kessler, B.; Ganteför, G.; Eberhardt, W. Photoemission Spectra of C60−: ElectronPhonon Coupling, Jahn-Teller Effect, and Superconductivity in the Fullerides. Phys. Rev. Lett. 1995, 74, 1875−1878. (72) Laflamme Janssen, J.; Côté, M.; Louie, S. G.; Cohen, M. L. Electron-phonon coupling in C60 using hybrid functionals. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 81, 073106. (73) Faber, C. Electronic, excitonic and polaronic properties of organic systems within the many-body GW and Bethe-Salpeter formalisms: towards organic photovoltaics. Ph.D. thesis, Université de Grenoble, 2014; https://tel.archives-ouvertes.fr/tel-01369304 (accessed June 13 2019). (74) Fetter, A. L.; Walecka, J. D. Quantum Theory of Many-Particle Systems; Dover: New York, 2003. (75) Mattuck, R. A. Guide to Feynman Diagrams in the Many-body Problem; Dover Books on Physics Series; Dover Publications: New York, 1976. (76) Hedin, L. New Method for Calculating the One-Particle Green’s Function with Application to the Electron-Gas Problem. Phys. Rev. 1965, 139, A796−A823. (77) Gross, E.; Runge, E.; Heinonen, O. Many-Particle Theory; Taylor & Francis: London, 1991. O


Article

Journal of Chemical Theory and Computation (78) Hedin, L. On correlation effects in electron spectroscopies and the GW approximation. J. Phys.: Condens. Matter 1999, 11, R489. (79) Economou, E. Green’s Functions in Quantum Physics; Springer Series in Solid-State Sciences; Springer: Berlin, Heidelberg, 2013. (80) Bechstedt, F. Many-Body Approach to Electronic Excitations: Concepts and Applications; Springer: Berlin, Heidelberg, 2015. (81) Friedrich, C.; Schindlmayr, A. Many-Body Perturbation Theory: The GW Approximation; NIC series; John von Neumann Institute for Computing: Jülich, 2006; Vol. 31, pp 335−355. (82) von Barth, U.; Holm, B. Self-consistent GW0 results for the electron gas: Fixed screened potential W0 within the random-phase approximation. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 8411−8419. (83) Shirley, E. L. Self-consistent GW and higher-order calculations of electron states in metals. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 7758−7764. (84) Holm, B.; von Barth, U. Fully self-consistent GW self-energy of the electron gas. Phys. Rev. B: Condens. Matter Mater. Phys. 1998, 57, 2108−2117. (85) Blase, X.; Attaccalite, C.; Olevano, V. First-principles GW calculations for fullerenes, porphyrins, phtalocyanine, and other molecules of interest for organic photovoltaic applications. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 83, 115103. (86) van Setten, M. J.; Gremaud, R.; Brocks, G.; Dam, B.; Kresse, G.; de Wijs, G. A. Optical response of the sodium Alanate system: GW0BSE calculations and thin film measurements. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 83, 035422. (87) Bruneval, F.; Marques, M. A. L. Benchmarking the Starting Points of the GW Approximation for Molecules. J. Chem. Theory Comput. 2013, 9, 324−329. (88) Pham, T. A.; Nguyen, H.-V.; Rocca, D.; Galli, G. GW calculations using the spectral decomposition of the dielectric matrix: Verification, validation, and comparison of methods. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 87, 155148. (89) Mark, C. In Recent Advances in Density Functional Methods; Delano, C., Ed.; Recent Advances in Computational Chemistry; World Scientific, 1995; Vol. 1; pp 155−192. (90) Marques, M. A. L.; Gross, E. K. U. In A Primer in Density Functional Theory; Fiolhais, C., Nogueira, F., Marques, M. A. L., Eds.; Springer: Berlin, Heidelberg, 2003; pp 144−184. (91) Shao, M.; Lin, L.; Yang, C.; Liu, F.; Da Jornada, F. H.; Deslippe, J.; Louie, S. G. Low rank approximation in G0W0 calculations. Sci. China Math. 2016, 59, 1593−1612. (92) Guo, Q.; Zhou, C.; Ma, Z.; Ren, Z.; Fan, H.; Yang, X. Elementary Chemical Reactions in Surface Photocatalysis. Annu. Rev. Phys. Chem. 2018, 69, 451−472. (93) Bässler, H.; Köhler, A. Unimolecular and Supramolecular Electronics I: Chemistry and Physics Meet at Metal-Molecule Interfaces; Springer: Berlin, Heidelberg, 2012; pp 1−65. (94) Oberhofer, H.; Reuter, K.; Blumberger, J. Charge Transport in Molecular Materials: An Assessment of Computational Methods. Chem. Rev. 2017, 117, 10319−10357. (95) Dürr, H. Photochromism: Molecules and Systems; Elsevier: San Diego, CA, 2003. (96) Kuwahara, S.; Suzuki, Y.; Sugita, N.; Ikeda, M.; Nagatsugi, F.; Harada, N.; Habata, Y. Thermal E/Z Isomerization in First Generation Molecular Motors. J. Org. Chem. 2018, 83, 4800−4804. (97) Simeth, N. A.; Crespi, S.; Fagnoni, M.; König, B. Tuning the Thermal Isomerization of Phenylazoindole Photoswitches from Days to Nanoseconds. J. Am. Chem. Soc. 2018, 140, 2940−2946. (98) Samson, J. A. R. Line Broadening in Photoelectron Spectroscopy. Rev. Sci. Instrum. 1969, 40, 1174−1177. (99) Citrin, P. H.; Eisenberger, P.; Hamann, D. R. Phonon Broadening of X-Ray Photoemission Linewidths. Phys. Rev. Lett. 1974, 33, 965−969. (100) Chen, W.; Wu, X.; Car, R. X-Ray Absorption Signatures of the Molecular Environment in Water and Ice. Phys. Rev. Lett. 2010, 105, 017802.

(101) Patrick, C. E.; Giustino, F. Quantitative Analysis of Valence Photoemission Spectra and Quasiparticle Excitations at Chromophore-Semiconductor Interfaces. Phys. Rev. Lett. 2012, 109, 116801. (102) Vinson, J.; Jach, T.; Müller, M.; Unterumsberger, R.; Beckhoff, B. Quasiparticle lifetime broadening in resonant X-Ray scattering of NH4NO3. Phys. Rev. B: Condens. Matter Mater. Phys. 2016, 94, 035163. (103) Schleife, A.; Varley, J. B.; Fuchs, F.; Rödl, C.; Bechstedt, F.; Rinke, P.; Janotti, A.; Van de Walle, C. G. Tin dioxide from first principles: Quasiparticle electronic states and optical properties. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 83, 035116. (104) Vogel, D. J.; Kilin, D. S. First-Principles Treatment of Photoluminescence in Semiconductors. J. Phys. Chem. C 2015, 119, 27954−27964. (105) Akola, J.; Manninen, M.; Häkkinen, H.; Landman, U.; Li, X.; Wang, L.-S. Photoelectron spectra of aluminum cluster anions: Temperature effects and ab initio simulations. Phys. Rev. B: Condens. Matter Mater. Phys. 1999, 60, R11297−R11300. (106) Müller, M.; Sánchez-Portal, D.; Lin, H.; Brivio, G. P.; Selloni, A.; Fratesi, G. Effect of Structural Fluctuations on Elastic Lifetimes of Adsorbate States: Isonicotinic Acid on Rutile(110). J. Phys. Chem. C 2018, 122, 7575−7585. (107) Curchod, B. F. E.; Martínez, T. J. Ab Initio Nonadiabatic Quantum Molecular Dynamics. Chem. Rev. 2018, 118, 3305−3336. (108) Crespo-Otero, R.; Barbatti, M. Recent Advances and Perspectives on Nonadiabatic Mixed Quantum−Classical Dynamics. Chem. Rev. 2018, 118, 7026−7068. (109) Li, X.; Tully, J. C.; Schlegel, H. B.; Frisch, M. J. Ab initio Ehrenfest dynamics. J. Chem. Phys. 2005, 123, 084106. (110) Alonso, J. L.; Andrade, X.; Echenique, P.; Falceto, F.; PradaGracia, D.; Rubio, A. Efficient Formalism for Large-Scale Ab Initio Molecular Dynamics based on Time-Dependent Density Functional Theory. Phys. Rev. Lett. 2008, 101, 096403. (111) Svoboda, O.; Ončaḱ , M.; Slavíček, P. Simulations of light induced processes in water based on ab initio path integrals molecular dynamics. I. Photoabsorption. J. Chem. Phys. 2011, 135, 154301. (112) Humeniuk, A.; Wohlgemuth, M.; Suzuki, T.; Mitrić, R. Timeresolved photoelectron imaging spectra from non-adiabatic molecular dynamics simulations. J. Chem. Phys. 2013, 139, 134104. (113) Curchod, B. F. E.; Sisto, A.; Martínez, T. J. Ab Initio Multiple Spawning Photochemical Dynamics of DMABN Using GPUs. J. Phys. Chem. A 2017, 121, 265−276. (114) Crespo-Otero, R.; Barbatti, M. Spectrum simulation and decomposition with nuclear ensemble: formal derivation and application to benzene, furan and 2-phenylfuran. Theor. Chem. Acc. 2012, 131, 1237. (115) de Sousa, L. E.; Ribeiro, L. A.; Fonseca, A. L. d. A.; da Silva Filho, D. A. Modeling the Emission Spectra of Organic Molecules: A Competition between Franck−Condon and Nuclear Ensemble Methods. J. Phys. Chem. A 2016, 120, 5380−5388. (116) Kossoski, F.; Barbatti, M. Nuclear Ensemble Approach with Importance Sampling. J. Chem. Theory Comput. 2018, 14, 3173−3183. (117) A first-principles materials simulation code using DFT. https://launchpad.net/siesta (accessed June 14 2019). (118) Nosé, S. A unified formulation of the constant temperature molecular dynamics methods. J. Chem. Phys. 1984, 81, 511−519. (119) Hoover, W. G. Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31, 1695−1697. (120) Schlick, T. Molecular Modeling and Simulation: An Interdisciplinary Guide; Springer: New York, 2010. (121) Bergsma, J. P.; Berens, P. H.; Wilson, K. R.; Fredkin, D. R.; Heller, E. J. Electronic spectra from molecular dynamics: a simple approach. J. Phys. Chem. 1984, 88, 612−619. (122) Garbuio, V.; Cascella, M.; Reining, L.; Sole, R. D.; Pulci, O. Ab Initio Calculation of Optical Spectra of Liquids: Many-Body Effects in the Electronic Excitations of Water. Phys. Rev. Lett. 2006, 97, 137402. (123) Ashcroft, N.; Mermin, N. Solid State Physics; HRW international editions; Holt, Rinehart and Winston: New York, 1976. P


Article

Journal of Chemical Theory and Computation (124) Zhukov, V. P.; Chulkov, E. V.; Echenique, P. M. GW + T theory of excited electron lifetimes in metals. Phys. Rev. B: Condens. Matter Mater. Phys. 2005, 72, 155109. (125) Salaneck, W. R.; Duke, C. B.; Eberhardt, W.; Plummer, E. W.; Freund, H. J. Temperature-Dependent Ultraviolet Photoemission Linewidths of Molecular Solids: Isopropyl Benzene. Phys. Rev. Lett. 1980, 45, 280−283. (126) Weightman, P. X-Ray-excited Auger and photoelectron spectroscopy. Rep. Prog. Phys. 1982, 45, 753. (127) Sevian, H. M.; Skinner, J. L. A molecular theory of inhomogeneous broadening, including the correlation between different transitions, in liquids and glasses. Theor. Chem. Acc. 1992, 82, 29−46. (128) Wilson, E.; Decius, J.; Cross, P. Molecular Vibrations: The Theory of Infrared and Raman Vibrational Spectra; Dover Books on Chemistry Series; Dover Publications: New York, 1980. (129) Giustino, F. Electron-phonon interactions from first principles. Rev. Mod. Phys. 2017, 89, 015003. (130) Landau, L.; Lifshitz, E. Statistical Physics; Pergamon Press Ltd.: Amsterdam, 1980; Vol. 5. (131) Blanco, M. A.; Flórez, M.; Bermejo, M. Evaluation of the rotation matrices in the basis of real spherical harmonics. J. Mol. Struct.: THEOCHEM 1997, 419, 19−27. (132) Barbry, M.; Koval, P.; Marchesin, F.; Esteban, R.; Borisov, A. G.; Aizpurua, J.; Sánchez-Portal, D. Atomistic Near-Field Nanoplasmonics: Reaching Atomic-Scale Resolution in Nanooptics. Nano Lett. 2015, 15, 3410−3419. (133) Koval, P.; Marchesin, F.; Foerster, D.; Sánchez-Portal, D. Optical response of silver clusters and their hollow shells from linearresponse TDDFT. J. Phys.: Condens. Matter 2016, 28, 214001. (134) Barbry, M. Plasmons in Nanoparticles: Atomistic Ab Initio Theory for Large Systems. Ph.D. thesis, University of Basque Country, Donostia-San Sebastián, Spain, 2018; http://cfm.ehu.es/view/files/ MArc_barbry_2-1.pdf (accessed June 14 2019). (135) Rieger, M. M.; Steinbeck, L.; White, I.; Rojas, H.; Godby, R. The GW space-time method for the self-energy of large systems. Comput. Phys. Commun. 1999, 117, 211−228. (136) Blase, X.; Ordejón, P. Dynamical screening and absorption within a strictly localized basis implementation of time-dependent LDA: From small clusters and molecules to aza-fullerenes. Phys. Rev. B: Condens. Matter Mater. Phys. 2004, 69, 085111. (137) Foerster, D.; Koval, P. On the Kohn-Sham density response in a localized basis set. J. Chem. Phys. 2009, 131, 044103. (138) Talman, J. D. NumSBT: A subroutine for calculating spherical Bessel transforms numerically. Comput. Phys. Commun. 2009, 180, 332−338. (139) Caruso, F.; Rinke, P.; Ren, X.; Rubio, A.; Scheffler, M. Selfconsistent GW: All-electron implementation with localized basis functions. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 88, 075105. (140) Curtiss, L. A.; Redfern, P. C.; Raghavachari, K.; Pople, J. A. Assessment of Gaussian-2 and density functional theories for the computation of ionization potentials and electron affinities. J. Chem. Phys. 1998, 109, 42−55. (141) Curtiss, L. A.; Redfern, P. C.; Raghavachari, K.; Pople, J. A. G2/97 Test Set. http://www.cse.anl.gov/OldCHMwebsiteContent/ compmat/g2-97.htm (accessed May 24 2018). (142) Dunning, T. H. Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. J. Chem. Phys. 1989, 90, 1007−1023. (143) Bruneval, F. Accurate many-body perturbation theory calculations of the electronic structure of molecules. http://web. archive.org/web/20181025222541/http://www.cse.anl.gov/ OldCHMwebsiteContent/compmat/g2-97.htm (accessed June 14 2019). (144) Hamann, D. R. Optimized norm-conserving Vanderbilt pseudopotentials. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 88, 085117.

(145) Foerster, D. Elimination, in electronic structure calculations, of redundant orbital products. J. Chem. Phys. 2008, 128, 034108. (146) Irie, M.; Fukaminato, T.; Matsuda, K.; Kobatake, S. Photochromism of Diarylethene Molecules and Crystals: Memories, Switches, and Actuators. Chem. Rev. 2014, 114, 12174−12277. (147) Tanaka, S.; Toba, M.; Nakashima, T.; Kawai, T.; Yoshino, K. Photoelectron Spectroscopy Study of Photoresponsive Change in Electronic Structure of Amorphous Photochromic Diarylethene Film. Jpn. J. Appl. Phys. 2008, 47, 1215. (148) Jakobsson, F. L. E.; Marsal, P.; Braun, S.; Fahlman, M.; Berggren, M.; Cornil, J.; Crispin, X. Tuning the Energy Levels of Photochromic Diarylethene Compounds for Opto-Electronic Switch Devices. J. Phys. Chem. C 2009, 113, 18396−18405. (149) Frisch, J.; Herder, M.; Herrmann, P.; Heimel, G.; Hecht, S.; Koch, N. Photoinduced reversible changes in the electronic structure of photochromic diarylethene films. Appl. Phys. A: Mater. Sci. Process. 2013, 113, 1−4. (150) Nickel, F.; Bernien, M.; Herder, M.; Wrzalek, S.; Chittas, P.; Kraffert, K.; Arruda, L. M.; Kipgen, L.; Krüger, D.; Hecht, S.; Kuch, W. Light-induced photoisomerization of a diarylethene molecular switch on solid surfaces. J. Phys.: Condens. Matter 2017, 29, 374001. (151) Wang, Q.; Frisch, J.; Herder, M.; Hecht, S.; Koch, N. Electronic Properties of Optically Switchable Photochromic Diarylethene Molecules at the Interface with Organic Semiconductors. ChemPhysChem 2017, 18, 722−727. (152) Darghouth, A. A. M. H. M.; Casida, M. E.; Taouali, W.; Alimi, K.; Ljungberg, M. P.; Koval, P.; Sánchez-Portal, D.; Foerster, D. Assessment of Density-Functional Tight-Binding Ionization Potentials and Electron Affinities of Molecules of Interest for Organic Solar Cells Against First-Principles GW Calculations. Computation 2015, 3, 616.

Q


Toward Efficient GW Calculations Using Numerical Atomic Orbitals

Recommend Documents