Generalization of Dielectric-Dependent Hybrid Functionals to Finite Systems

The accurate prediction of electronic and optical properties of molecules and solids is a persistent challenge for methods based on density functional theory. We propose a generalization of dielectricdependent hybrid functionals to finite systems where the definition of the mixing fraction of exact and semilocal exchange is physically motivated, nonempirical, and system dependent. The proposed functional yields ionization potentials, and fundamental and optical gaps of many, diverse molecular systems in excellent agreement with experiments, including organic and inorganic molecules and semiconducting nanocrystals. We further demonstrate that this hybrid functional gives the correct alignment between energy levels of the exemplary TTF-TCNQ donor-acceptor system.


I. INTRODUCTION
Current first-principles methods to predict the electronic and optical properties of solids and nanostructures are mainly based on density functional theory (DFT) [1,2] and many-body perturbation theory (MBPT) [3][4][5][6].In the case of molecules, methods adopted for their optoelectronic properties include ab initio quantum chemistry (QC) techniques [7][8][9][10] and hybrid density functionals [11][12][13][14].While QC and MBPT methods have proven accurate for several classes of molecules and solids [15,16], they remain computationally more demanding than DFT-based techniques; they are limited to smaller-sized systems, and they have rarely been applied in conjunction with ab initio molecular dynamics or Monte Carlo simulations.
While calculations of electronic excitation energies of solids rely mostly on single-particle energies, obtained within either DFT or MBPT, those for molecules encompass single-particle energies and total energy differences.A popular and rather accurate method, especially for small molecules, is the so-called delta self-consistent field method (ΔSCF) [17]; ionization potentials and electron affinities are obtained as total energy differences of the neutral and charged molecule, where an electron is either removed or added.Total energies may be obtained at different levels of theory, either with DFT or QC techniques.However, the ΔSCF method is limited to obtaining bound frontier orbital energies, and it cannot be easily generalized to solids [18].In addition, it is not suitable for the calculation of energy-level alignment between different molecules or nanoparticles in close proximity with each other, and hence to study charge transfer processes.
Therefore, the search for density functionals yielding accurate electronic properties and, in particular, singleparticle energies of broad classes of condensed and finite systems is an active field of research, impacting applications ranging from materials discovery to the interpretation of spectroscopic and transport measurements.Functionals adopted by the condensed matter physics community in the early days of DFT applications to solids were strictly dependent on the electronic density [19,20]; those employed by the chemistry community were instead often defined using linear combinations of Hartree-Fock exchange and density-dependent exchange functionals, which are referred to as hybrid functionals.The latter have also seen many implementations and uses [11,[21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38] in condensed matter physics.One successful example of hybrid functionals is that of the optimally tuned rangeseparated hybrid [39], defined using parameters selfconsistently tuned to enforce the validity of Koopman's theorem [40] for molecules.However, the generalization of these functionals to solids is not a simple task [41].Another interesting example is that of Koopman's-complaint functionals [42].Motivated by the COHSEX approximation used within MBPT [43], an alternative class of hybrid functionals for solids, referred to as dielectric-dependent hybrids, were introduced by equating the optimal mixing fraction of exact exchange with the inverse of the bulk dielectric constant [44][45][46][47][48][49][50][51].The generalization of dielectric hybrid functionals to finite systems is not straightforward, as there is no unique definition of an average dielectric constant for molecules or nanoparticles.In fact, the concept of a dielectric constant of finite systems is ill defined since the ratios of applied and induced fields differ, depending on the location and form of the probe.We note that several definitions of dielectric constants of nanocrystals have been proposed in the literature, some dependent on volume, yielding different results for the same system [52,53], depending on the chosen definition of the volume.
Here, we present a generalization of dielectric-dependent hybrid functionals to molecules and nanostructures, where we overcome the problem of introducing a dielectric constant of finite systems, by using a mixing fraction of exact and semilocal exchange defined as the ratio of the screened and unscreened exchange energies.Such a ratio reduces to the inverse dielectric constant in the bulk limit as shown below.We show that the proposed functional yields ionization potentials, and fundamental as well as optical gaps of molecular systems, in excellent agreement with experiments, including organic and inorganic molecules and semiconducting nanocrystals.We further demonstrate that this hybrid functional yields the correct alignment between energy levels of the exemplary Tetrathiafulvalene-Tetracyanoquinodimethane (TTF-TCNQ) donor-acceptor system.

II. THEORY
As mentioned above, the concept of a dielectric constant is ill defined for finite systems and so is the volume.Hence, in order to generalize dielectric-dependent hybrid functionals to molecules, and, in general, to finite systems, we do not attempt to generalize the concept of a dielectric constant used for condensed systems.Rather, we start from approximate forms of the self-energy obtained within MBPT, at the GW level [43], and we define a nonempirical system-dependent fraction of exact and semilocal exchange.
To obtain the mixing fraction of exact exchange, we approximate the quasiparticle GW self-energy assuming that the screened Coulomb interaction W is static (i.e., frequency independent): Here, γ is a small positive parameter, Wðr 0 ; rÞ ¼ R dr 00 ϵ −1 ðr 0 ; r 00 Þvðr 00 ; rÞ is the statically screened Coulomb interaction, vðr 00 ; rÞ is the Coulomb potential, ϵ −1 ðr 0 ; r 00 Þ is the static inverse dielectric matrix, and Gðr; r 0 ; τÞ is the Green's function.Using the Kohn-Sham wave functions and eigenvalues, the Green's function can be written as where η is a small positive quantity, ε F is the Fermi energy, ϕ j is the jth Kohn-Sham single-particle wave function, and ε j is the jth Kohn-Sham eigenvalue.Here, ϵ −1 may be evaluated at different levels of theory; in this paper, we used the random phase approximation [4], with where χ is the full polarizability.Within RPA, χ is given by χ ¼ ð1 − vχ 0 Þ −1 χ 0 , where χ 0 is the independent-particle polarizability given by with f i the occupation factor.We adopted a recently developed method to compute ϵ −1 , which does not require the explicit calculation of unoccupied single-particle electronic states nor any explicit matrix inversion [54][55][56].
Under the static approximation, Eq. ( 1) becomes [57] Σ GW ðr; where is the screened exchange self-energy and N is the number of occupied electronic states.We then assume that the nonlocal screened exchange term is proportional to the Hartree-Fock exchange, as obtained using KS singleparticles: where Σ X ðr; r 0 Þ is given by From Eq. ( 6), we define the state-dependent screened exchange constant α SX i : Finally, we define the global screened exchange constant (SXC) that will be used as the mixing fraction of exact exchange: Note that α SX [Eq.( 9)] does not depend on any arbitrary volume definitions.The exchange correlation potential entering the generalized Kohn-Sham formalism [58] is obtained as where v x and v c are semilocal exchange correlation potentials, here chosen to be described within the PBE [20] approximation.Hereafter, we refer to this functional as the SX hybrid functional.We note that Eq. ( 10) could be evaluated self-consistently, along the lines proposed in Ref. [48].We also note that the use of Eq. ( 10) with the definition given in Eq. ( 9) is not size consistent.Size consistency will be addressed in future investigations.
In the case of condensed systems, if every element of the inverse dielectric matrix ϵ −1 G;G 0 is replaced by the inverse of the macroscopic dielectric constant ϵ −1 ∞ , the screened exchange constant reduces to the inverse of the macroscopic dielectric constant, and the SX functional reduces to the dielectric-dependent hybrid functional.

III. COMPUTATIONAL METHODS
In this section, we summarize the properties computed in this work using the SX functional along with the computational details of our calculations.We report electronic properties corresponding to charged and neutral excitations, in particular, vertical ionization potentials and fundamental gaps, typically measured in photoemission experiments, and absorption gaps and spectra.For both types of excitations, we first determined the ground-state total energy and single-particle energies within DFT, at the PBE level of theory.The screened exchange constant [Eq.( 8)] was then evaluated for the highest occupied molecular orbital and used to define the mixing fraction of exact exchange entering the global hybrid definition [Eq.(10)].
Vertical ionization potentials and electron affinities of the G2 set and of TTF and TCNQ were computed using the single-particle energies obtained from the self-consistent solution of the Kohn-Sham equations with the global hybrid defined by Eq. (10).The fundamental gaps of the Si nanoparticles and the vertical ionization potentials of Thiel's test set were obtained in the same way.Absorption properties (i.e., optical gaps) of the Thiel's test set were instead computed using time-dependent density functional theory, as implemented in Quantum ESPRESSO in the turboTDDFT code [59].We used the latter code to solve the TDDFT equations in the adiabatic approximation, where the exchange-correlation potential was described by the functional of Eq. (10).We note that the evaluation of empty electronic states is not necessary in TDDFT calculations, when following the procedure of Ref. [59].
As mentioned previously, all results reported in this work were obtained using the mixing parameter from Eq. ( 8), where i is the index of the highest occupied molecular orbital (HOMO).The values of α SX i¼HOMO and α SX from Eq. ( 9) were found to be very similar for all the systems studied here; for example, we found a maximum difference (α SX i¼HOMO − α SX of 0.06 for the G2=97 test set [60] described below, with an average difference of 0.02.In the case of silicon nanoparticles, the difference between α SX i¼HOMO and α SX decreases as the size of the system increases, leading to α SX i¼HOMO ≈ α SX for Si 86 H 76 .We note that the application of the SX functional consists of a single hybrid calculation for neutral systems, unlike the functionals of Refs.[39,40], which require multiple calculations of neutral and ionized systems.In addition, the use of SX is computationally less demanding than GW calculations, as frequency-dependent dielectric matrices are not required.Once the screened exchange constant is obtained (using, e.g., the WEST code [56]), SX calculations of absorption and photoemission properties may be performed by any DFT code, where global hybrid functionals and TDDFT methods are available.

IV. RESULTS
Below, we discuss the optoelectronic properties of inorganic and organic molecules as obtained using the SX functional; the latter is aimed at the accurate determination of response properties of finite systems rather than at ground-state properties such as total energies.First, we study the ionization potentials and fundamental gaps of several molecules and nanoparticles and compare our results with those obtained with MBPT techniques and several other global hybrid methods.Several photoemission spectra of selected systems are also included in Ref. [61].With the exception of the deepest excitations, we found good agreement with both G 0 W 0 results and experiment over a wide range of energies (see Ref. [61]).We then used Eq. ( 10) to compute optical properties within the framework of time-dependent density functional theory, and we show below that, remarkably, the SX functional yields excellent results for both photoemission and optical gaps of a variety of molecules.
A common benchmark for GW calculations [6] is the evaluation of vertical ionization potentials of different subsets of the G2=97 data set [62][63][64].Here, we chose a subset consisting of 36 closed-shell [65] molecules [56].The calculated vertical ionization potentials, i.e., the negative of the highest occupied molecular orbital energies obtained using Eq. ( 10) and the mean absolute error (MAE) with respect to the experimental values, are presented in Fig. 1 and Table I.All GW calculations were performed using the WEST code and the hybrid functional calculations using the Quantum ESPRESSO software package [56,66]; we employed norm-conserving pseudopotentials.A wave function energy cutoff of 85 Ry was used for all molecular systems, except Si 35 H 36 , Si 66 H 64 , and Si 87 H 64 , for which we used 25 Ry.Molecular geometries for the G2 test set and experimental vertical ionization potentials were taken from the NIST Computational Chemistry Comparison and Benchmark DataBase [67].Molecular geometries for Thiel's test set were taken from Ref. [68].Single-particle energies were corrected using the Makov-Payne correction scheme [69].Nanoparticle geometries, from Refs.[70,71], are provided in Ref. [61].
In most cases, results obtained with the SX functional lie between those obtained at the G 0 W 0 level of theory, with either PBE or SX wave functions, and they are in very good agreement with the experimental values and the G 0 W 0 @PBE0 results.In order to compare the results of the different levels of theory in a statistically meaningful way, we used a Wilcoxon signed-rank test [73] as implemented in the R package [74], and we compared the absolute errors of the SX functional and the popular G 0 W 0 @PBE method; we found a p-value [61] less than 0.043 and a 95% confidence interval of [0.01, 0.28] eV, indicating that the SX functional is more accurate than G 0 W 0 @PBE, on average, for molecules of which this subset is representative.We note that popular hybrid functionals such as PBE0 and B3LYP have fixed mixing fractions (α) below 0.50, rather different from the average value provided by Eq. ( 8) for the molecules studied here, i.e., α ¼ 0.72 with a range of Δα ¼ 0.19.This average mixing fraction is similar to the optimal mixing fraction reported in a previous study for the G3=99 test set [75].In addition, Eq. ( 9) also yields similar mixing fractions as those of Ref. [76] and provides a physical justification to the results reported by Atalla et al. [76], which were obtained by self-consistently fitting hybrid functional results to those of GW calculations.
A smaller fraction of exact exchange than given by Eq. ( 8) severely overestimates the amount of screening present in molecular systems, which in turn results in the poor performance, on this benchmark, of global hybrids such as B3LYP and PBE0, as well as the range-separated hybrid HSE06.However, we note that these and some semilocal functionals may yield more accurate predictions when total energy differences of neutral and ionized systems are considered [60,62], instead of the eigenvalues of the neutral molecule; this is especially so for B3LYP, originally designed to yield total energies not accurate response properties from orbital energies.For example, ΔSCF calculations at the PBE level of theory yield a MAE of 0.24 eV for the G2 test set [62].
Small mixing fractions have also been shown to yield incorrect energy-level alignments for several molecular systems.For example, it was found that functionals with α < 0.30 predict spurious amounts of charge transfer for the TTF (C 6 H 4 S 4 ) and TCNQ (C 12 H 4 N 4 ) donor-acceptor system [76].A correct charge transfer is obtained if the difference (Δ) between the lowest unoccupied molecular orbital (LUMO) of TCNQ and the HOMO of TTF is positive, in the limit of infinite separation between the two molecules.The mixing fractions predicted by Eq. ( 8) for TTF and TCNQ are 0.63 and 0.61, respectively.Evaluated at the average mixing fraction (0.62), the SX functional ΔSCF(PBE) G 0 W 0 @PBE G 0 W 0 @PBE0 SX G 0 W 0 @SX FIG. 1. Vertical ionization potentials (VIP) of the G2=97 subset consisting of 36 closed-shell molecules computed as the highest occupied molecular orbital energy using the SX hybrid functional (SX), G 0 W 0 starting from PBE wave functions (G 0 W 0 @PBE; data from Ref. [56]) and SX wave functions (G 0 W 0 @SX), and total energy differences [ΔSCFðPBEÞ; data from Ref. [62]].TABLE I. MAE in the vertical ionization potential of 36 closedshell molecules predicted by the highest occupied molecular orbital energy computed using the SX functional (SX) (see text), G 0 W 0 starting from PBE0 wave functions (G 0 W 0 @PBE0), the SX hybrid functional (G 0 W 0 @SX), and PBE (G 0 W 0 @PBE).We also show the results for Hartree-Fock exchange with PBE correlation (EXXc), B3LYP [72], PBE0 [21], and HSE06 [22] for comparison.predicted Δ > 0, with ε TCNQ LUMO ¼ −4.32 and ε TTF HOMO ¼ −6.15 eV.The mixing fraction that minimizes the deviation from the straight-line error was shown to be 0.70 for TCNQ and 0.8 for TTF [77,78], not dissimilar from the values predicted by Eq. ( 8).Benzene is an additional example where a larger mixing fraction (α > 0.25) was reported to be advantageous.In Ref. [79], a mixing fraction of ¼ 0.70 for benzene was found to minimize the many-electron self-interaction error.Using Eq. ( 8), we found α SX i¼HOMO ¼ 0.67.We now investigate how different levels of theory, including Eq. ( 10), perform as the size of the system is varied.In Fig. 2, we present the fundamental gaps of several hydrogen passivated silicon nanoparticles predicted using the SX functional.For most of the molecules in Fig. 1, the results obtained using Eq. ( 10) lie between those of G 0 W 0 @PBE and G 0 W 0 @SX; these findings are similar to those of Fig. 2, with the addition that the quasiparticle corrections to the SX eigenvalues decrease as the size of the system increases.In contrast, we found the fundamental gaps predicted by ΔSCFðPBEÞ to fall consistently below those of G 0 W 0 @PBE by an average value of 0.67 eV.Also provided in Fig. 2 are the mixing values predicted by Eq. ( 8), which increase as the size of the system decreases.This variation is related to the magnitude of electronic screening [70,71] and of the quasiparticle corrections, consistent with predictions of previous works [56,80,81]; we note that the SX hybrid proposed here is the first global hybrid functional to yield such an a priori prediction of the magnitude of quasiparticle energy corrections as a function of system size.

Method MAE (eV)
The improvement of the accuracy of adiabatic timedependent density functional theory is an area of active research in condensed matter physics and physical chemistry [59,[82][83][84][85][86][87]; hence, we explored the application of Eq. ( 10) to time-dependent density functional theory (TDSX), and we show below that the SX functional can correctly reproduce not only ionization energies and fundamental gaps but also optical gaps.In Fig. 3, we compare the first lowest optical singlet excitations predicted by TDSX with the revised best theoretical estimates for Thiel's test set, as obtained from reference literature values and quantum chemistry calculations, e.g., coupled cluster [88].Thiel's test set consists of 28 molecules and was recently used for benchmarking state-of-the-art many-body perturbation theory techniques based on the Bethe-Salpeter-equation (BSE) and the GW approximation [89,90].Jacquemin et al. [91] found TDPBE0 (timedependent density functional calculations with the PBE0 functional) to be exceptionally accurate for this set, with the lowest mean absolute error compared to 28 exchange correlation functionals.The mean absolute errors of the lowest optical singlet transitions predicted by TDSX are presented in Fig. 4, which shows that the SX functional performs almost as well as TDPBE0 for optical transitions.Also presented in Fig. 4 are the vertical ionization potentials (i.e., the highest occupied molecular orbital energies) for the same molecules.Importantly, the SX functional captures the physics of both quasiparticle and excitonic effects, yielding accurate photoemission and optical properties.
A statistical analysis of the results for both the photoemission and absorption properties presented in Fig. 4 is provided in Ref. [61].Although TDPBE0 appears to outperform almost every other technique to describe optical excitations, the single-particle energies of the same functional (PBE0) fail to describe vertical ionization potentials, predicted using the SX hybrid functional (see text), G 0 W 0 starting from PBE wave functions (G 0 W 0 @PBE), and SX wave functions (G 0 W 0 @SX).Fundamental gaps obtained from total energy differences [ΔSCFðPBEÞ] are included for comparison.TDPBE0 BSE@G 0 W 0 @PBE0 TDSX FIG. 3. The first lowest optical singlet excitations predicted by the TDSX compared to the revised best theoretical estimates (see text) of Thiel's test set for 28 molecules.Values predicted for the same transitions by TDPBE0 (time-dependent density functional calculations with the PBE0 functional) and BSE@G 0 W 0 @PBE0 (results from the solution of the BSE, starting from single-particle energies computed at the G 0 W 0 level of theory, with DFT results obtained with the PBE0 functional) are taken from Ref. [89].
as seen in the G2=97 and the Thiel's test sets.This result suggests that TDPBE0 underestimates both quasiparticle and excitonic effects, resulting in a fortuitous cancellation of errors for these systems.We find that G 0 W 0 @PBE0 is accurate for vertical ionization potentials, while the mean absolute error of optical excitations obtained from the solution of the BSE is larger than that of TDSX.Remarkably, the SX functional captures the physics of both quasiparticle and excitonic effects, yielding accurate photoemission and optical properties.

V. CONCLUSION
In summary, we generalized the concept of dielectricdependent functionals to finite systems by defining a global hybrid functional (which we called the SX functional) with a mixing parameter between exact and local exchange that is nonempirical and system dependent.We showed that the SX functional yields accurate results for both photoemission and optical properties of a variety of organic and inorganic molecules; in addition, it yields the most accurate combined results (charged and neutral excitations) compared to other global hybrids and methods based on MBPT, while being computationally cost-effective as compared to G 0 W 0 calculations.We also demonstrated that the SX functional gives the correct alignment between energy levels of the exemplary TTF-TCNQ donor-acceptor system, which many first-principles or semiempirical DFT-based schemes proposed in the literature failed to reproduce.Finally, we presented a detailed statistical analysis of the data discussed in our paper, in order to assess the accuracy of the proposed functional for molecules beyond the test sets investigated in this work.TDSX BSE@G 0 W 0 @PBE0 TDPBE0 SX G 0 W 0 @PBE0 PBE0 FIG. 4. MAE for (from left to right) the ionization potentials (photoemission) of Thiel's test set (with the exception of propanamide because of the lack of experimental data) predicted by the highest occupied molecular orbital energy computed using PBE0, G 0 W 0 @PBE0, and the screened exchange constant functional (SX); the lowest optical singlet transitions predicted by the SX functional [TDSX, i.e., time-dependent density functional calculations using the SX functional, Eq. ( 10)], TDPBE0 and PBE0 with many-body perturbation theory corrections (BSE@G 0 W 0 @PBE0) TDPBE0 and BSE@G 0 W 0 @PBE0 are taken from Ref. [89].