Accurate optical spectra through time-dependent density functional theory based on screening-dependent hybrid functionals

We investigate optical absorption spectra obtained through time-dependent density functional theory (TD-DFT) based on nonempirical hybrid functionals that are designed to correctly reproduce the dielectric function. The comparison with state-of-the-art $GW$ calculations followed by the solution of the Bethe-Sapeter equation (BSE-$GW$) shows close agreement for both the transition energies and the main features of the spectra. We confront TD-DFT with BSE-$GW$ by focusing on the model dielectric function and the local exchange-correlation kernel. The present TD-DFT approach achieves the accuracy of BSE-$GW$ at a fraction of the computational cost.

Kohn-Sham density functional theory (DFT) is notoriously bad for describing band gaps in semiconductors due to the lack of the derivative discontinuity in semi-local functionals [1,2].Thus, most of the ab inito methods for optical absorption calculations based on DFT electronic structures have to address two important problems.First, it is necessary to correct the band gap.Second, the interaction between electrons and holes has to be taken into account.The state-of-the-art approach to improve the band gap is the GW approximation [3][4][5][6], whereas the electron-hole interaction can be included by solving the Bethe-Salpeter equation (BSE) [7][8][9][10][11].The combined GW +BSE approach has been shown to give very accurate results compared to experiment, but the main drawback is the scaling, which makes it computationally very challenging for large size systems.Gross developed an alternative approach based on the time-dependent electron density, typically referred to as time-dependent density functional theory (TD-DFT) [7,12,13].In this approach, the time-dependent Kohn-Sham equations include a time-dependent exchangecorrelation (xc) potential v xc and its variation with the time-dependent density, also known as the exchangecorrelation kernel f xc .The exact v xc and f xc are unknown, but several approximations have been introduced.Local approximations to f xc lack the correct long wavelength limit, f xc (q → 0) ∝ 1/q 2 , responsible for the correct description of the electron-hole interaction.Therefore, local approximations are unable to capture excitonic effects, but can perform well for metallic systems [14][15][16][17].The correct asymptotic behavior is recovered in the so-called "nanoquanta" kernel [18][19][20] derived from the BSE to capture excitonic effects.Hence this approach produces accurate spectra for solids, but remains computationally as expensive as solving the BSE.Hybrid functional calculations with non-local Fock exchange can be used to improve band gaps.Moreover, since the long wavelength limit is accounted for, it can be expected that these functionals could be used for calculating optical spectra.Various hybrid functionals have been tested in TD-DFT and it has been shown that a good performance can be achieved in molecules [21,22].However, a good description in solids requires the consideration of the screening in the exchange interaction [23,24].Such a screened interaction was found to be crucial for the correct description of optical spectra [23,25].Several different approximations for the screening of the non-local exchange interaction have been investigated [26][27][28][29].The results suggest that hybrid functionals yield spectra comparable to BSE-GW provided the adopted fraction of Fock exchange accounts for the screening in the long range.In particular, Wing et al. obtained good results with screened range-separated hybrid functionals [29].However, the correct screening in the short and medium range was not imposed but rather followed from the empirical setting of the hybrid functional parameters in their TD-DFT approach.For instance, their choice of taking 25% of Fock exchange in the short range does not describe the physically correct behavior of the screening.More importantly, the range separation parameter was empirically tuned so that the calculated band gaps matched the GW ones.Recently, Chen et al. [30] developed a nonempirical hybrid functional scheme, in which all the parameters are taken from the static screening without tuning.The method showed very accurate electronic structures and band gaps for a large number of semiconductors and insulators.The advantage of this approach is that it accurately accounts for the wave-vector dependent screening: at short-range the exchange interaction is only weakly screened, whereas in the long-range it is reduced by the static dielectric constant.In this work, we investigate the performance of hybridfunctional TD-DFT for optical absorption calculations through the comparison with state-of-the-art BSE-GW .In the TD-DFT scheme, we employ hybrid functionals that have been designed to reproduce the correct screen-ing properties through a self-consistent procedure [30].We show that this scheme provides absorption spectra with an accuracy comparable to that of BSE-GW without tuning parameters.In particular, we show that equivalent descriptions of the screening in TD-DFT and BSE-GW result in very similar absorption spectra.Following recent work on dielectric-dependent hybrid functionals (DDH) [30][31][32], we use in this work the explicit form of the exchange-correlation potential given by where µ is the range-separation parameter, V PBE x and V PBE c are the PBE exchange and correlation potentials [33].Here, V Fock x is the Fock exchange operator given by where ψ nk are one-electron Bloch states, the sum over k is over N k k points of the Brillouin zone, and the sum over n is over all the occupied bands.In Eq. ( 1), the Fock exchange interaction is thus modified by the function In the GW approximation, the Coulomb interaction is screened by the dielectric function which is a frequency dependent tensor −1 (r, r , ω) [6].Thus, c DDH x (|r − r |) in Eq. ( 3) corresponds to a model inverse dielectric function −1 model (|r − r |) that neglects the dynamical screening (ω = 0) and the off-diagonal elements.In the approach of Ref. [30], the parameters in Eq. ( 3) are determined self-consistently.In the long-wave limit, the interaction is set to 1/( ∞ r), where the dielectric constant ∞ is calculated using the random-phase approximation with vertex corrections.The parameter µ is obtained by fitting the model to the calculated dielectric function.In Fig. 1, the model dielectric functions associated with the DDH are compared to the diagonal elements of the dielectric matrix at the Γ point at zero frequency.The dielectric matrix is obtained within partially self-consistent GW with vertex corrections [cf.Supplemental Material (SM) [34]].In all the cases, we find the model dielectric function to be in good agreement with the calculated dielectric function.The excitation spectra in both BSE and TD-DFT are obtained by solving an eigenvalue problem, referred to as the Bethe-Salpeter and Casida equation, respectively [7,35]: where submatrices A and B read  3) with parameters taken from Ref. [30].
with the indices i, j and a, b referring to occupied and unoccupied states, respectively.The excitation frequencies of the system are given by Ω. X and Y are the two-body electron-hole eigenstates in the transition basis ψ a (r)ψ * i (r ) and ψ i (r)ψ * a (r ).Matrix A includes two terms, the energy of the direct transition from occupied to unoccupied states and the electron-hole interaction described by the kernel K [36,37].Equation ( 4) is non-Hermitian, which makes it difficult to solve using standard eigenvalue solvers [38,39].A common practice to avoid this difficulty is to neglect the coupling between excitations and de-excitations by setting B to zero.This approximation is known as the Tamm-Dancoff approximation.The distinction between BSE-GW and TD-DFT approaches results, on the one hand, from the origin of the one-particle eigenfunctions and energies and, on the other hand, from the type of the interaction kernel K.
To make the comparison between different methods more transparent, we provide in Fig. 2 the Feynman diagrams corresponding to the various irreducible polarizabilities χ discussed in this work.In BSE-GW , the orbitals and energies are derived from a preceding GW calculation and the kernel consists of a Hartree term V and a screened exchange term W [37]: The Hartree term describes the bare Coulomb interac- tion and is the same in all the approximations considered here.It can be included straightforwardly in a two-point formulation involving the polarizability χ.The exchange term, however, requires calculating four-point integrals, which drastically increases the complexity of the problem.The screening of the exchange interaction is determined by the frequency-dependent dielectric function obtained from GW and is represented by a vertical wiggly line in the diagrams.However, as it is shown in Refs.[40][41][42], the dynamical effects can often be neglected in BSE calculations.
In the TD-DFT approach, the electron energies and wave functions are obtained from a semilocal or hybridfunctional calculation.The interaction kernel consists of three terms, a Hartree and an exchange term like in the BSE, and an additional local exchange-correlation interaction f loc xc [43]: where The screening of the exchange interaction in TD-DFT is provided by a constant c  [44] for C, from [45] for Si, from [46] for SiC, from [47] for Ar, from [48] for NaCl, and from [49] for MgO.
obtained from the eigenvalue problem in Eq. ( 4).In particular, the BSE-GW calculations are based on partially self-consistent GW using the "nanoquanta" vertex corrections f xc in the polarizability χ [50].These two approaches are tested on a set of materials possessing a wide range of band gaps.The corresponding spectra are given in Fig. 3. Our calculations show that both approaches agree well with experiment and that TD-DDH reproduces all the spectral features with the correct oscillator strengths.In the case of C, Si and SiC, the spectra are nearly on top of each other.For Ar, NaCl, and MgO, the relative positions of the main features in the spectra are found to be shifted slightly.In principle, this shift can result from differences in the band structure and in the screening.Our analysis indicates that the dominant effect is due to the energy transition terms in Eq. ( 5).From Table I, we notice that the calculated band gaps differ by less than 0.15 eV for Si, SiC, and C, but that the disagreement is more substantial for NaCl, Ar, and MgO.Overall, when compared to experimental values corrected for the coupling to phonons, DDH and GW band gaps show mean average errors of 0.11 and 0.22 eV, respectively.These errors are consistent with the current accuracy of ab initio methods [51,52], indicating that the agreement with experiment should be considered excellent for both schemes.As far as the screening is concerned, we show below that the small discrepancies observed in Fig. 1 hardly change the spectra.[53], with a correction of 0.06 eV from Ref. [54].b Ref. [55], with a correction of 0.11 eV from Ref. [56].c Ref. [57], with a correction of 0.37 eV from Ref. [54].d Ref. [48], with a correction of 0.17 eV from Ref. [58].e Ref. [59], with a correction of 0.03 eV calculated using the method described in Refs.[60,61].f Ref. [62], with a correction of 0.53 eV from Ref. [63].
To compare the screening in BSE-GW and TD-DDH, we show in Fig. 4 the absorption spectra of diamond calculated in various approximations using the same energies and wave functions, which are taken from a G 0 W 0 calculation.We start our analysis from a BSE calculation in which the full static inverse dielectric matrix is used (W full ).In particular, we show that the off-diagonal elements of this matrix barely have any effect on the calculated spectrum (W diag ).Next, we replace the inverse dielectric function with the model c DDH x (|r − r |) and find no discernible difference in the spectrum (W m ).This treatment of the screening is equivalent to that of TD-DDH, where the local exchange-correlation kernel is neglected, i.e. f loc xc = 0, also referred to as model BSE (mBSE) [64].To restore the full TD-DDH screening, we include the f loc xc and still obtain essentially the same spectrum (W m + f loc xc ).Hence, these results indicate that the model screening in TD-DDH gives an accurate description of the screening in W = −1 V and that the effect of f loc xc is negligible for extended systems.The numerical complexity of Eq. ( 4) is the same in BSE and TD-DDH.However, the preceding GW calculations required in BSE-GW involve a high computational cost, which scales like N 4 in the number of electrons N in most GW implementations instead of like N 3 in TD-DFT.Additionally, in the calculation of the Green's function in GW , the convergence with respect to the number of unoccupied states and the number of frequency points has to be controlled carefully, which significantly increases the complexity of the calculations.Note that the static dielectric constant only needs to be determined at the Γ point of the Brillouin zone and that it converges quickly with respect to the number of included orbitals.Furthermore, the hybrid-functional approach only requires a model static dielectric function, for which the static limit can be obtained rather efficiently [31].Thus, the hybrid functional approach opens the way to more efficient numerical schemes that can circumvent the calculation of the full dielectric matrix.
In conclusion, we have shown that time-dependent calculations using the parameter-free DDH functional yield optical absorption spectra with an accuracy comparable to BSE-GW .The success of this approach originates from the use of a model dielectric function that gives a physically motivated description of the screened exchange interaction over the full spatial range.Notably, the computational complexity of the method is drastically reduced compared to BSE-GW , as it eliminates the need for preceding GW calculations.This will allow one to consider larger and more complex systems than hitherto possible.Support from the Swiss National Science foundation is acknowledged under Grant No. 200020-172524.The calculations have been performed at the Swiss National Supercomputing Centre (CSCS) (grant under project ID s879) and at SCITAS-EPFL.
FIG.1.Inverse dielectric functions vs. wave vector G for Si, C, SiC, Ar, NaCl and MgO, as calculated in GW and given by the model in Eq. (3) with parameters taken from Ref.[30].

FIG. 4 .
FIG. 4. Optical absorption spectrum of diamond calculated with BSE and full −1 from GW (W full ), BSE and diagonal −1 (W diag ), BSE and model −1 model (W m ), TD-DDH and model −1 model with f loc xc (W m + f loc xc ).All spectra are based on energies and wave functions from a G0W0 calculation.A 6 × 6 × 6 k-point grid is used.

TABLE I .
[33] gaps (in eV) obtained with the semilocal PBE functional[33], DDH, and GW .The experimental values are augmented by theoretical corrections resulting from the coupling to phonons..23 a 2.53 b 5.85 c 9.14 d 14.33 e 8.36 f a Ref.