Accurate Ionization Potentials and Electron Affinities of Acceptor

Jan 5, 2016 - Journal of Chemical Theory and Computation 2016 12 (2), 595-604. Abstract | Full ... The LDA-1/2 method implemented in the exciting code...
1 downloads 0 Views 2MB Size
Article pubs.acs.org/JCTC

Accurate Ionization Potentials and Electron Affinities of Acceptor Molecules II: Non-Empirically Tuned Long-Range Corrected Hybrid Functionals Lukas Gallandi,† Noa Marom,‡ Patrick Rinke,§ and Thomas Körzdörfer*,† †

Institut für Chemie, Universität Potsdam, Karl-Liebknecht-Straße 24-25, 14476 Potsdam, Germany Physics and Engineering Physics, Tulane University, New Orleans, Louisiana 70118, United States § COMP/Department of Applied Physics, Aalto University School of Science, P.O. Box 11100, FI-00076 Aalto, Finland ‡

S Supporting Information *

ABSTRACT: The performance of non-empirically tuned long-range corrected hybrid functionals for the prediction of vertical ionization potentials (IPs) and electron affinities (EAs) is assessed for a set of 24 organic acceptor molecules. Basis setextrapolated coupled cluster singles, doubles, and perturbative triples [CCSD(T)] calculations serve as a reference for this study. Compared to standard exchange-correlation functionals, tuned long-range corrected hybrid functionals produce highly reliable results for vertical IPs and EAs, yielding mean absolute errors on par with computationally more demanding GW calculations. In particular, it is demonstrated that long-range corrected hybrid functionals serve as ideal starting points for non-self-consistent GW calculations.

1. INTRODUCTION The accurate prediction of ionization potentials (IPs) and electron affinities (EAs) from first principles has been a longstanding challenge for computational materials science.1−7 In particular in the fields of molecular electronics and organic photovoltaics, the energetics of charged excitations are key factors determining the functionality and efficiency of devices. Given the plethora of compounds that could possibly be of interest for such applications, employing computational data mining approaches to identify the most promising candidates for further experimental studies is highly desirable. However, such an approach requires an accurate, yet computationally efficient, electronic structure method. In terms of its numerical efficiency, density functional theory (DFT) stands out from all commonly available electronic structure methods.8,9 Its superior accuracy-to-numerical-costs ratio is one of the reasons why today it has become the most commonly used approach for large-scale data mining projects based on electronic structure theory.10−12 Another reason is that, at least in principle, DFT provides easy access to charged excitation energies. According to the IP-theorem,2,3 the eigenvalue corresponding to the highest occupied molecular orbital (HOMO) of exact Kohn−Sham9 (KS) or generalized Kohn−Sham13 (GKS) theory equals the first vertical IP, i.e., the total energy difference between the neutral and cationic molecule for a fixed geometry. There exists, however, no corresponding theorem for all other eigenvalues, including the lowest unoccupied molecular orbital (LUMO). Yet, the IPtheorem can also be used to determine the vertical EA, at least for finite systems such as molecules and clusters, simply by © 2016 American Chemical Society

calculating the HOMO eigenvalue of the anion. It is a wellknown fact, however, that commonly used semilocal and global hybrid functionals do not obey the IP-theorem, a problem that has been attributed to the self-interaction error14 (SIE) or localization/delocalization error15,16 of DFT.14,17−23 As a consequence, the HOMO eigenvalue obtained from semilocal or global hybrid functionals is often only a very poor approximation to the first ionization potential. The most commonly used approach to avoid this problem is the Δ-selfconsistent-field (ΔSCF) method, in which the vertical IP (EA) is determined from the total energy difference between the neutral and cationic (anionic) molecules. In fact, ΔSCF values for IPs and EAs are widely considered to be much more reliable than interpreting the KS eigenvalues directly.24−26 In a series of recent publications by several groups, it was demonstrated that the IP-theorem can be restored by using long-range corrected (LC) hybrid functionals, in which the range-separation parameter is non-empirically tuned for each system individually.27−52 By comparison to experiments and GW calculations, the HOMO and LUMO eigenvalues obtained from this new class of functionals have been demonstrated to be excellent approximations to IPs and EAs.29−31,33,52 However, a rigorous assessment of the quality of the approach based on a direct comparison to experimental data suffers from several drawbacks. First, well-resolved gas-phase photoelectron measurements of molecules relevant to organic electronics are rare. This is particularly true for the photoelectron spectra of Received: September 11, 2015 Published: January 5, 2016 605

DOI: 10.1021/acs.jctc.5b00873 J. Chem. Theory Comput. 2016, 12, 605−614

Article

Journal of Chemical Theory and Computation molecular anions, which could be used directly to benchmark electron affinities in the gas phase.53 Second, experimental IPs are typically obtained from the peak maxima of the (0−0) peaks of vibrationally broadened photoelectron spectra at some finite temperature, which are clearly different from the vertical ionization energies calculated from electronic structure theory. While the vibrational structure of photoelectron peaks can be considered by virtue of a Franck−Condon analysis,54−56 this approach introduces a number of approximations that might bias the comparison of IPs and EAs such as the determination of the relaxation energy and vibrational frequencies from approximate DFT functionals. To avoid these problems all together, an alternative approach is to compare the predicted vertical IPs and EAs to those obtained from more accurate (and typically also computationally more demanding) electronic structure methods. This approach has been followed in a number of publications,31,33,34,37,47,57 in which the results obtained from IPtuned LC-hybrid functionals were compared to those obtained from many-body perturbation theory in the GW approximation,5,6 where G stands for the Green’s function and W is the screened Coulomb interaction between quasiparticles, i.e., electrons or holes with their polarization cloud. However, GW calculations for molecules suffer from a number of problems on their own, which are discussed in detail in part III of this study.58 One of these problems is that, due to the high computational costs and limited accuracy of fully self-consistent GW calculations, one often only calculates the quasiparticle energies as one-shot corrections to KS or GKS eigenvalues. The quality of the results obtained from such G0W0 calculations is mainly determined by the quality of the DFT starting point.59−66 For the case of molecules, using different starting points can lead to differences in the predicted vertical IPs and EAs as large as 1 eV. Hybrid functional starting points typically improve the accuracy of G0W0 quasiparticle energies,63,67 and non-empirical tuning procedures for determining the optimal fraction of HF-exchange in global hybrid functionals have been suggested.65,66 In a recent benchmark for G0W0 methods, using the range-separated hybrid functional CAM-B3LYP as a starting point was demonstrated to yield very accurate results.67 Given their accuracy for the prediction of quasiparticle energies, tuned long-range corrected hybrid functionals evidently lend themselves as a potentially improved starting point for G0W0.57 This work is part of a collaborative effort of several groups, whose primary goal is to compare the accuracy of different families of electronic structure methods for the prediction of vertical IPs and EAs.58,68,69 To this end, the vertical IPs and EAs for a test set consisting of 24 organic acceptor molecules (Figure 1) were determined using basis set-extrapolated coupled cluster singles doubles with perturbative triples (CCSD(T)) calculations.68 Here, we will use these results to benchmark the quality of the vertical IPs and EAs predicted from the non-empirically tuned LC-hybrid functional LCωPBE.70,71 In a number of previous studies, we have found that, after the tuning of the range-separation parameter, the IPs and EAs obtained from different LC-hybrid functionals are very similar. We thus anticipate the LC-ωPBE results to be representative also for other LC-hybrid functionals. To compare with the performance of standard exchangecorrelation functionals, we also assess a prototypical semilocal functional, PBE,72,73 and a global hybrid functional with 25% HF-exchange, PBE0.74 As it is already well established in the literature that HOMO and LUMO eigenvalues obtained from

Figure 1. Benchmark set of molecules.

semilocal and standard global hybrid functionals are poor approximations to vertical IPs and EAs, respectively, we here focus on IPs and EAs obtained from a ΔSCF approach. It is found that tuned LC-ωPBE provides mean absolute errors (MAEs) that clearly improve upon standard exchangecorrelation functionals and are on par with commonly used GW methods, which are studied and compared in detail in ref 58. Further, we analyze the performance of standard and tuned LC-ωPBE as an ideal starting point for G0W0.

2. METHODOLOGY The main idea underlying LC-hybrid functionals is the separation of the Coulomb operator into short-range and long-range components, which is typically achieved by virtue of the standard error function 1 − erf(ωr ) erf(ωr ) 1 = + r r  r    SR

(1)

LR

This range-separation ansatz allows the combination of the semilocal or global hybrid description of short-range (SR) interactions with the correct asymptotics of full Hartree−Fock exchange for long-range (LR) interactions. Specifically, the LCωPBE functional combines the SR ωPBE functional75,76 with LR HF-exchange and semilocal PBE correlation,72 i.e., ExcLC − ωPBE = ExωPBE,SR + ExHF,LR + EcPBE

(2)

Here, we follow an approach in which the range-separation parameter ω is found for each system independently by virtue 606

DOI: 10.1021/acs.jctc.5b00873 J. Chem. Theory Comput. 2016, 12, 605−614

Article

Journal of Chemical Theory and Computation of a non-empirical tuning procedure.27,33 Different tuning procedures have been proposed for the accurate prediction of IPs and fundamental gaps.27,29,31,33 The basic idea of IP-tuning is to minimize the difference between the HOMO eigenvalue and the total energy difference between the electronic groundstates of the neutral (N electrons) and cationic (N−1 electrons) molecules at a fixed geometry, i.e., ΔIP = | − ϵωHOMO − [Egs(ω , N ) − Egs(ω , N − 1)]|

3. NUMERICAL DETAILS We implemented the LC-ωPBE functional in the all-electron numerical atom-centered orbital (NAO) code FHI-aims,77−79 with which all DFT and GW calculations that entered the benchmark were performed. The NAO basis in FHI-aims is constructed from a set of atomic-orbital basis functions for the core and valence electrons, followed by a hierarchy of additional basis functions that are categorized as tier 1−4. Although DFT calculations require much smaller basis sets than GW, we have used the same basis set (tier 4) and grid settings (tight) for all calculations to allow for an unbiased comparison to the GW calculations performed in this work and in part III of this study.58 For GW calculations, this basis set gives energies converged to within 0.1 eV (see part III of this study58 for an evaluation of basis set convergence of GW calculations in FHIaims). We used the molecular geometries that were employed in part I of this study.68 The IP- and gap-tuning procedures to determine the optimal value for the range-separation parameter ω were carried out using linear interpolation and a simple inverse parabolic interpolation, respectively. In our implementation, the tuning procedure is automated using a script that acts as a wrapper around FHI-aims. Typically, this procedure converges quickly after only a few self-consistency cycles. For all tuning runs, we used a tier 3 basis set, which we have verified to be sufficiently large to converge the ω-value up to 10−3 a0−1. In terms of the bare computation time, the tuning procedure is, hence, less costly than a G0W0 calculation, although it requires several self-consistency steps to converge. In practice, the tuning procedure can be applied to much larger systems than GW since it requires significantly less memory and since its results are much less sensitive to the size of the basis set. A detailed description of the hybrid-DFT and GW implementation in FHI-aims is provided elsewhere.79 The fractional particle curves for the evaluation of the delocalization error were calculated with a local development version of the PSI4 code.80

(3)

In this way, the IP-tuning procedure enforces the IP-theorem, i.e., the HOMO eigenvalue predicted by IP-tuned functionals is as close as possible to the ΔSCF value for the IP. Similarly, one can tune ω such that the LUMO eigenvalue of the neutral molecule equals the ΔSCF-predicted EA, i.e., ΔEA = | − ϵωLUMO − [Egs(ω , N + 1) − Egs(ω , N )]|

(4)

In an attempt to improve the performance of the functional for the prediction of the fundamental gap (IP-EA), it has been suggested to tune ω such that ΔGAP =

Δ2IP + Δ2EA

(5)

27,34

is minimized. This procedure will be referred to as gaptuning in the following. Specifically, we will apply the IP- and gap-tuning procedures to the LC-ωPBE functional; the resulting methods will be denoted as LC-ωPBE(IP) and LCωPBE(gap), respectively. As the gap-tuning procedure involves SCF calculations for the neutral, cationic, and anionic molecules, it is computationally more expensive than the simpler IP-tuning, for which only the neutral and cationic states need to be considered. In this work, the IPs and EAs for LC-ωPBE(IP) and LCωPBE(gap) as well as for the standard LC-ωPBE (with a fixed range-separation parameter of ω = 0.4 a−1 0 ) are predicted from the HOMO and LUMO eigenvalues, respectively. This is different from PBE and PBE0, for which meaningful IP and EA results can only be obtained from ΔSCF. Of course, for LCωPBE(IP) the ΔSCF value for the IP equals the HOMO eigenvalue as a consequence of the IP-tuning. Strictly speaking, this does not hold for the LUMO/EA and also not for the IPs and EAs predicted from LC-ωPBE(gap). However, as will be discussed in more detail below, also in these cases the difference between HOMO/LUMO eigenvalues and ΔSCF IPs and EAs are small and, in most cases, negligible. In the G0W0 approach,5,6 the quasiparticle energies, ϵQP i , are obtained from a one-shot perturbative correction to the GKS eigenvalues, ϵGKS , as i GW ϵQP = ϵGKS + ⟨ϕiGKS|Σ̂ xc (ϵQP ̂ |ϕiGKS⟩ i i i ) − vxc

4. RESULTS The errors for the IPs, EAs, and fundamental gaps with respect to CCSD(T) of PBE(ΔSCF), PBE0(ΔSCF), LC-ωPBE, LCωPBE(IP), LC-ωPBE(gap), G0W0@LC-ωPBE, G0W0@LCωPBE(IP), and G0W0@LC-ωPBE(gap) are presented as heat maps in Figure 2 and as histograms in Figures 3−5. All numbers are tabulated in the Supporting Information,81 where we also provide the IP- and gap-tuned ω-values for all molecules in the test set. The performance of the individual methods for IPs, EAs, and fundamental gaps is evaluated in terms of their mean absolute errors (MAEs) as well as the standard deviation (STD) of the errors with respect to the CCSD(T) benchmark. These numbers are provided in Figures 3−5. While the MAEs give an indication of the overall performance of each method, the STDs provide additional insight such as systematic deviations from the reference data. 4.1. Ionization Potentials. The semilocal functional PBE yields the largest MAE (0.67 eV) and STD (0.29 eV) among all studied methods. Generally speaking, IPs are significantly underestimated by as much as 1 eV. This might come as a surprise since ΔSCF approaches are often considered to provide quite reliable results, even when using semilocal functionals. The results can, however, be significantly improved by incorporating a fraction of HF-exchange as in the PBE0

(6)

where ϕGKS are the GKS orbitals, Σ̂GW is the exchangei xc correlation self-energy in the GW approximation calculated from the GKS orbitals and eigenvalues, and v̂xc is the exchangecorrelation operator, which for the case of hybrid functionals combines a local exchange-correlation potential with the nonlocal HF operator v̂HF x . Specifically, we use vxĉ = vxωPBE,SR + vxHF,LR + vcPBE ̂

(7)

The GKS orbitals and eigenvalues are determined from a standard, IP-tuned, or gap-tuned LC-ωPBE calculation. We will refer to these methods as G0W0@LC-ωPBE, G0W0@LCωPBE(IP), and G0W0@LC-ωPBE(gap), respectively. 607

DOI: 10.1021/acs.jctc.5b00873 J. Chem. Theory Comput. 2016, 12, 605−614

Article

Journal of Chemical Theory and Computation

Figure 2. Heatmap representation of the errors in the vertical IPs (top), EAs (middle), and fundamental gaps (bottom) predicted by the different methods as compared to CCSD(T).

Figure 3. Histograms of the errors in the vertical IPs predicted by the different methods as compared to the CCSD(T) benchmark. The mean absolute error (MAE) and the standard deviation (STD) of the error for each method are also shown.

and 0.24 eV) and STDs (0.23 and 0.21 eV), both clearly improving upon PBE, PBE0, and standard LC-ωPBE in terms of the MAE. Similarly to PBE and PBE0, the tuned LC-ωPBE approaches tend to underestimate IPs, although the underestimation is less pronounced. The STD, however, is somewhat larger than in PBE0 and also in the standard LC-ωPBE, which is certainly surprising given the fact that these functionals have been specifically tuned for each system independently. Comparison to part III58 demonstrates that the accuracy of IP-tuned LC-hybrid functionals can keep up with and even surpass most of the commonly used GW methods, including

functional, reducing the MAE and STD to 0.35 and 0.13 eV, respectively. The improvement of PBE0 over PBE has been attibuted to the reduction of self-interaction errors by the inclusion of HF-exchange.74,82 The histograms shown in Figure 3 illustrate that PBE0(ΔSCF) also systematically underestimates IPs, but not nearly as much as PBE(ΔSCF). Examining the tuned LC-ωPBE results, it is apparent that the tuning procedure significantly improves the accuracy of the GKS HOMO eigenvalues. Also, we find that the difference in the tuning procedures does not significantly affect the quality of the results. Both tuning procedures yield similar MAEs (0.23 608

DOI: 10.1021/acs.jctc.5b00873 J. Chem. Theory Comput. 2016, 12, 605−614

Article

Journal of Chemical Theory and Computation

Figure 4. Histograms of the errors in the vertical EAs predicted by the different methods as compared to the CCSD(T) benchmark. The mean absolute error (MAE) and the standard deviation (STD) of the error for each method are also shown.

non-tuned variants of G0W0@LC-ωPBE are small and, hence, should not be overinterpreted. Also, it important to note that due to the fact that we compare to a CCSD(T) benchmark all molecules in our test set are small. For small and/or saturated molecules, the standard range-separation parameter is often close to the IP-tuned one. This is, however, not the case for large, conjugated systems.30 Hence, the favorable cancelation of errors seen for this test set of molecules might not be effective for larger π-conjugated sytems. Finally, we note that the LC-ωPBE(IP) and G0W0@LCωPBE(IP) results in Figure 3 indicate the existence of four outliers for which the method predicts too large IPs. An explanation for these outliers is provided in the Further Discussion section below. 4.2. Electron Affinities. Quite unexpectedly, the MAEs found for the EAs are significantly smaller than for the IPs for all methods. In fact, all methods result in very similar MAEs of roughly 0.2 eV (Figure 4). While the smallest MAE (0.16 eV) is found for the two tuned LC-ωPBE approaches, the largest MAE (0.22 eV) is found for PBE0(ΔSCF). Also the STDs for all methods are significantly reduced with respect to what was found for the IPs. This is particularly apparent for the G0W0 approaches, which show an impressively small STD of 0.05 and 0.06 eV, respectively. This means that there is barely any variation in the errors with respect to CCSD(T). As a consequence, the MAE of 0.2 eV found for G0W0@LCωPBE(IP) and G0W0@LC-ωPBE(gap) translates to a systematic overestimation of EAs by exactly this error. Generally speaking, all methods tend to overestimate the EA compared to CCSD(T), which is the opposite of what was found for the IPs. Again, the details of the tuning procedure do not have a significant influence on the results. In particular, LC-ωPBE(gap) and G0W0@LC-ωPBE(gap) do not clearly improve upon their IP-tuned counterparts, despite the fact that the EA was considered explicitely in the tuning procedure. In fact, the MAE for the non-tuned LC-ωPBE calculations is equal to those obtained from the two different tuning procedures. Clearly, this is again related to the fact that we here study a test set of small, π-conjugated organic molecules for which the tuning procedure leads to range-separation parameters that do not differ much from the standard range-separation parameter. One molecule that is worth exploring in more detail is TCNQ. As shown in Figure 1, in particular the tuned LCωPBE and G0W0@LC-ωPBE methods overestimate the EA of

partially and fully self-consistent GW methods. This is an impressive result that clearly underlines the significance of the tuning concept. The low MAEs obtained for both tuned LC-ωPBE variants can be even further reduced by applying a G0W0 correction to the GKS eigenvalue spectrum. This way, one reaches MAEs as low as 0.13 eV for both G0W0@LC-ωPBE(IP) and G0W0@LCωPBE(gap), respectively. The excellent performance of the tuned G0W0@LC-ωPBE methods is further supported by the small STDs of 0.18 eV. Generally speaking, IPs obtained from G0W0@LC-ωPBE(IP) and G0W0@LC-ωPBE(gap) are still somewhat underestimated as compared to the CCSD(T) benchmark. Nevertheless, the accuracy of both methods is unrivaled among all non-self-consistent, partially self-consistent, and fully self-consistent GW schemes studied in part III.58 While G0W0@PBE, G0W0@PBE0, and G0W0@HF yield MAEs of 0.65, 0.29, and 0.43 eV, respectively, an MAE of 0.60 eV is found for fully self-consistent GW. The performance of the LCωPBE starting point for IPs cannot even be challenged when going beyond the GW approximation by including higher order correction terms in many-body perturbation theory. We interpret this finding as the consequence of a favorable cancelation of errors between the mean-field starting point and the neglect of vertex corrections in the GW approximation.63,83 This is further supported by the finding that for the test set of small molecules considered here the non-tuned variant of G0W0@LC-ωPBE with a fixed range-separation parameter of 0.4 a−1 0 results in an even lower MAE (0.09 eV) and a smaller STD (0.06 eV) than what is obtained when using tuned rangeseparation parameters. Given the fact that the tuned LC-ωPBE functionals yield much more accurate IPs (and also EAs, see below) than the non-tuned variant, this is certainly surprising. Again, these results underline that the search for an optimal starting point for G0W0 has to rely on a favorable cancelation between the errors in the mean-field starting point and the errors introduced by using the GW approximation. As a consequence, a range-separation parameter that yields optimal IPs and EAs from LC-ωPBE may not necessarily yield optimal results when used as a starting point for G0W0@LC-ωPBE. We have observed a similar behavior in the past for global hybrids, where the optimal admixture of exact exchange in PBE0 is different than it is in [email protected],66 In fact, however, the differences in the MAEs found between the different tuned and 609

DOI: 10.1021/acs.jctc.5b00873 J. Chem. Theory Comput. 2016, 12, 605−614

Article

Journal of Chemical Theory and Computation

Figure 5. Histograms of the errors in the vertical gaps predicted by the different methods as compared to the CCSD(T) benchmark. The mean absolute error (MAE) and the standard deviation (STD) of the error for each method are also shown.

looking at the gap can, indeed, provide for a more stringent test case than benchmarking the individual IP and EA values. Although the large MAE of LC-ωPBE (0.54 eV) is reduced by the tuning procedures, the MAEs are still as large as 0.38 eV. The G0W0@LC-ωPBE method yields the smallest errors of all methods with an MAE 0.09 eV and an STD of 0.11 eV, respectively. Comparison of the tuned LC-ωPBE approaches demonstrates that the details of the tuning procedure are not important, which means that a standard IP-tuning is typically sufficient and the computationally more expensive gap-tuning is not necessary (note that possible exceptions from this finding will be discussed in the next section). Both tuned G0W0@LCωPBE methods systematically underestimate the gap by roughly 0.25 eV, which is a direct consequence of the overestimation of the EA shown in Figure 4. Especially for the case of G0W0@LC-ωPBE(IP), one can clearly identify four outliers for which the gap is significantly overestimated. The reason for the occurrence of these outliers will be discussed in the next section.

TCNQ significantly. This indicates that a correct description of the EA of TCNQ is very challenging and probably requires the inclusion of correlation effects that go beyond what is included in the DFT and G0W0 methods studied here. In fact, the challenge of accurately predicting the electron affinity of TCNQ has been stressed earlier.84 In ref 58, it is demonstrated that the inclusion of vertex corrections to the GW approximation in the form of second-order screened exchange (SOSEX)85 leads to a satisfactory agreement with the CCSD(T) benchmark. This identifies TCNQ as a valuable test system for beyond GW approaches to the calculation of EAs. Yet, three of the methods studied here predict reasonable EAs for TCNQ, that is, PBE(ΔSCF), LC-ωPBE, and G0W0@ LC-ωPBE. While the good result obtained from LC-ωPBE can be regarded as an outlier from the typically rather poor performance of this method to predict EAs, the excellent performance of G0W0@LC-ωPBE with an error of only 0.06 eV with respect to the CCSD(T) benchmark is certainly remarkable. Another interesting aspect of Figure 4 is that, generally speaking, the G0W0 correction to the tuned LC-ωPBE eigenvalues slightly worsens the agreement with the CCSD(T) results. Both this finding and the excellent performance of G0W0@LC-ωPBE for TCNQ can only be rationalized by a systematic cancelation of errors between the non-selfconsistency and the neglect of vertex corrections in the GW approximation.63,83 This is further supported by the finding that fully self-consistent GW overestimates EAs significantly with an MAE of 0.61 eV.58 As the tuned LC-ωPBE functionals already provide a very good estimate for IPs and EAs, the desired partial cancelation of errors deteriorates. 4.3. Fundamental Gaps. For applications in molecular electronics and photovolatics, one is sometimes more interested in the fundamental gap (IP-EA) than the individual IP and EA values. Figure 5 shows the histogram representation of the errors in the vertical fundamental gaps when compared to the CCSD(T) benchmark. All methods studied in this work systematically underestimate the gap. The underestimation is particularly severe for the PBE(ΔSCF) approach, which shows an MAE of 0.89 eV and an STD of 0.37 eV, where the largest error is found for dinitro-benzonitrile (1.51 eV). The underestimation is somewhat reduced by using PBE0(ΔSCF), but the MAE of 0.56 eV is still large. This underlines that

5. FURTHER DISCUSSION The results presented above demonstrate that the nonempirically tuned LC-ωPBE and G0W0@LC-ωPBE methods predict highly accurate IPs and EAs. In particular after the G0W0 correction is applied, the variation in the errors becomes very small, thus underlining the generality of the approach. Generally speaking, all methods underestimate IPs and gaps and at the same time overestimate EAs. However, in the test set consisting of 24 molecules, we could also identify four outliers, for which in particular the LC-ωPBE(IP) and G0W0@LCωPBE(IP) methods overestimated the IPs and fundamental gaps: benzoquinone, naphtalenedione, maleic anhydride, and phthalimide (see heatmaps in Figure 2). In order to explain these results, it is helpful to recapitulate the motivation for using the IP-tuning procedure from a slightly different perspective. As was pointed out in several previous publications on this topic, the IP-tuning approach can be interpreted as a minimization of the delocalization error at the HOMO level.29,34,39,44,50,86 The delocalization error is inherent to standard semilocal and global hybrid functionals and can be quantified as the deviation of the total energy curve for fractional particle numbers from piecewise linearity, that is, 610

DOI: 10.1021/acs.jctc.5b00873 J. Chem. Theory Comput. 2016, 12, 605−614

Article

Journal of Chemical Theory and Computation

accordingly reduced when inserting an electron into the LUMO and increased when taking an electron out of the HOMO. The reason why benzonitrile and mDCNB behave so differently from benzoquinone and maleic anhydride can be understood by inspection of the isosurface plots of the HOMO and LUMO depicted in Figure 7. While the HOMO and

ΔE(N + δ) = E(N + δ) − E(N ) − [E(N + 1) − E(N )]δ (8)

As the HOMO eigenvalue equals the slope of the total energy evolution with fractional particle number δ, IP-tuning guarantees that the initial slope of the total energy curve at δ = 0− (that is, when going from the neutral to the cation state) equals the vertical IP (that is, the total energy difference between neutral and cation). Hence, the tuned functional obeys the piecewise linearity condition when taking out an electron from the HOMO. This is illustrated in Figure 6 (top), which

Figure 7. HOMO and LUMO of benzonitrile, mDCNB, benzoquinone, and maleic anhydride obtained from IP-tuned LC-ωPBE.

LUMO of benzonitrile and mDCNB are π-orbitals, the HOMO orbitals of benzoquinone and maleic anhydride are n-orbitals that include the lone electron pairs of the oxygen atoms. σ- and n-orbitals typically suffer from larger self-interaction errors than π-orbitals.20,87 As a consequence, these orbitals require larger amounts of HF-exchange or, in the case of LC-hybrids, larger ω-values. Hence, a balanced description of both π-orbitals and σ- or n-orbitals has been found to be problematic when using LC-hybrid functionals with no additional global amount of HFexchange such as the LC-ωPBE functional used here.37,57 By studying outer valence photoelectron spectra, it was demonstrated41,57 that the imbalanced description of occupied σ and π-orbitals when using tuned LC-hybrid functionals can be reduced significantly by including a fraction of HF-exchange in the SR part of the functional via

Figure 6. Deviation of the total energy from linearity as a function of fractional particle number δ for benzonitrile, mDCNB, benzoquinone, and maleic anhydride when inserting one electron into the HOMO (left) and LUMO (right) using IP- (upper) and gap-tuned (lower) LC-ωPBE.

shows the deviation from linearity of the total energy obtained from LC-ωPBE(IP), ΔE(N + δ), as a function of the fractional particle number δ for four of the molecules in the test set, i.e., benzonitrile, benzoquinone, mDCNB, and maleic anhydride. While benzoquinone and maleic anhydride are two of the above-mentioned outliers, benzonitrile and mDCNB are representative for the other 20 molecules for which LCωPBE(IP) and G0W0@LC-ωPBE(IP) produce consistent results. We have chosen those four molecules to point out an important difference between the outliers and the other molecules in the test set. IP-tuning guarantees that the piecewise linearity condition is obeyed when taking out an electron from the HOMO, which is why all fractional-particle energy curves in Figure 6 are perfectly flat for this case (top left). For the case of benzonitrile and mDCNB, IP-tuning also leads to an almost perfect linearity of the total energy when inserting an electron into the LUMO (top right), indicating that the IP-tuned ω-value not only minimizes ΔIP (eq 3) but, at least approximately, also ΔEA (eq 4). This is, however, not the case for benzoquinone and maleic anhydride. Here, the IP-tuned ω-values (0.391 a−1 0 and 0.423 a−1 0 , respectively) do not minimize ΔEA. Also, the ω-values obtained from the gap-tuning (0.346 a−1 and 0.372 a−1 0 0 , respectively) are substantially different from the IP-tuned ones. As a consequence, LC-ωPBE(gap) does not obey the piecewise linearity condition, as shown in Figure 6 (bottom). Instead, by minimizing Δgap (eq 5), gap-tuning provides for a trade-off between a minimal value for ΔIP and a minimum value for ΔEA. Compared to LC-ωPBE(IP), the deviation from linearity is

ExcLC − ωPBEh = αExHF,SR + (1 − α)ExωPBE,SR + ExHF,LR + EcPBE

(9)

Specifically, the authors suggested to tune the fraction of HFexchange (α) based on the piecewise linearity condition.57 For systems for which the tuning procedure is not applicable, it was suggested to empirically set α to 0.2.41 To test how such an approach would affect our benchmark, we included 20% HFexchange in our LC-ωPBE functional and studied how it affects the errors for the IPs, EAs, and gaps of benzonitrile, mDCNB, benzoquinone, and maleic anhydride. The respective numbers are provided in Table 1. It is found that the inclusion of SR HFexchange does not significantly affect the overall performance of the tuned LC-ωPBE and G0W0@LC-ωPBE methods. While the error in gaps of benzoquinone and maleic anhydride are slightly reduced for LC-ωPBE, the gap-errors for benzonitrile and mDCNB are slightly increased. Hence, we anticipate that the inclusion of SR HF-exchange in the LC-ωPBE functional, which was demonstrated to be essential to obtain trustworthy outer valence photoelectron spectra below the IP,41,57 does not have a significant effect on the results for the vertical IP and EA presented in this work. This also means that the description of the outer valence photoelectron spectra can be improved without sacrificing accuracy at the level of the HOMO and LUMO. 611

DOI: 10.1021/acs.jctc.5b00873 J. Chem. Theory Comput. 2016, 12, 605−614

Article

Journal of Chemical Theory and Computation

Table 1. IP-Tuned Range-Separation Parameters and Errors in the Vertical IP, EA, and Fundamental Gap Predicted by LCωPBEh(IP) and G0W0@LC-ωPBEh(IP) for Four Selected Molecules Compared to CCSD(T) When No HF-Exchange Is Employed in the SR-Part of the Functional (α = 0, left) and When 20% HF-Exchange Is Used in the SR (α = 0.2, right) α=0 ω LC-ωPBEh

G0W0@LC-ωPBEh

benzonitrile mDCNB benzoquinone maleic anhydride benzonitrile mDCNB benzoquinone maleic anhydride

[a−1 0 ]

0.278 0.270 0.391 0.423 0.278 0.270 0.391 0.423

α = 0.2

IPerr

EAerr

gaperr

−0.14 −0.18 0.42 0.39 0.00 0.01 0.40 0.37

0.13 0.15 0.01 −0.02 0.17 0.22 0.13 0.09

−0.28 −0.33 0.41 0.42 −0.17 −0.21 0.26 0.27

6. CONCLUSIONS In summary, we have assessed the performance of the standard and non-empirically tuned LC-hybrid functional LC-ωPBE for the prediction of vertical IPs and EAs by comparison to a recently derived benchmark of CCSD(T) calculations on 24 organic acceptor molecules. We found that tuned LC-ωPBE provides mean absolute errors that are on par with computationally more demanding GW methods. The predicted IPs can be further improved by applying a one-shot G0W0 correction to the LC-ωPBE eigenvalues, leading to mean absolute errors well below 0.2 eV, thus outperforming commonly used GW approaches. It was further demonstrated that the details of the tuning procedure (IP or gap-tuning) are typically not decisive for the prediction of accurate IPs, EAs, and gaps. The only exception are systems in which the HOMO and LUMO orbitals are of different characters (π vs σ/n). However, even for these challenging systems, tuned LC-ωPBE still yields errors that are on par with commonly used G0W0 calculations. With the electron affinity of the TCNQ molecule, we have also identified a problem that turns out to be very challenging to DFT and standard GW approaches and requires beyond-GW approaches or correlated wave function methods. Overall, our study demonstrates that non-empirically tuned range-separated hybrid functionals provide for a low-cost alternative to GW calculations for the prediction of charged excitation energies in molecules with no loss in accuracy. It was further demonstrated that the accuracy for the prediction of IPs and EAs can still be improved by applying a one-shot GW to the tuned LC-ωPBE eigenvalues, leading to even smaller MAEs and STDs. The most accurate (smallest MAE) and most consistent (smallest STD) results, however, are found in G0W0 calculations that are based on a LC-hybrid starting point that uses a fixed range-separation parameter of 0.4 a−1 0 . We attribute this finding to a favorable cancelation of errors between errors in the mean-field starting point and the missing vertex corrections in the GW approximation, as well as to the fact that our test set of molecules consists of small, organic molecules, for which the fixed range-separation parameter of 0.4 a−1 yields reasonable IPs and HOMO/ 0 LUMO gaps even without non-empirical tuning. Hence, it can be expected that the favorable cancelation of errors might not be effective for larger π-conjugated systems. For such systems, non-empirical tuning of the range-separation parameter might still be required to obtain reliable results from G0W0@LCωPBE. It should, however, be noted that the tuning procedure described here cannot straightforwardly be applied to periodic

ω

[a−1 0 ]

0.228 0.220 0.314 0.331 0.228 0.220 0.314 0.331

IPerr

EAerr

gaperr

−0.18 −0.20 0.42 0.31 −0.03 0.03 0.37 0.35

0.14 0.17 0.06 0.05 0.15 0.24 0.06 0.12

−0.32 −0.37 0.36 0.26 −0.18 −0.21 0.30 0.23

systems such as solids and interfaces for two reasons. First, the required calculations for cations are not trivial with periodic boundary conditions. Second, the piecewise linearity condition is obeyed even for local and semilocal functionals in the solid state.88,89 In other words, when periodic boundary conditions are applied, all exchange-correlation functionals obey the IPtheorem, but their HOMO eigenvalues may differ. In fact, it was found that, even in large but finite systems, the KS gap approaches the fundamental gap calculated from ΔSCF,90,91 a result that questions the validity of the tuning procedure for systems that are (much) larger than the ones studied here. In light of these unresolved questions, the exploration of alternative approaches to determine the optimal rangeseparation parameter for periodic as well as for large finite systems remains an interesting topic for future research.



ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.5b00873. All relevant numbers from the CCSD(T) benchmark; all calculated IPs, EAs, and band-gaps, as well as the corresponding errors, mean errors, and mean absolute errors for all molecules and all methods studied; all IPand gap-tuned range-separation parameters for the nonempirically tuned LC-ωPBE calculations. (PDF)



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS The authors thank F. Caruso, X. Ren, and V. Blum for helpful discussions. Work at the University of Potsdam was supported by the Deutsche Forschungsgemeinschaft (#KO-4876/1-1). Work at Tulane University was supported by the Louisiana Alliance for Simulation-Guided Materials Applications (LASiGMA), funded by the National Science Foundation (NSF) award number #EPS-1003897. Work at Aalto University was supported by the Academy of Finland through its Centres of Excellence Program (No. 251748 and 284621). Computer time was provided by the National Energy Research Scientific Computing Center (NERSC), which is supported by the Office of Science of the U.S. Department of Energy under contract DE-AC02-05CH11231. 612

DOI: 10.1021/acs.jctc.5b00873 J. Chem. Theory Comput. 2016, 12, 605−614

Article

Journal of Chemical Theory and Computation



(37) Körzdörfer, T.; Parrish, R. M.; Marom, N.; Sears, J. S.; Sherrill, C. D.; Brédas, J.-L. Phys. Rev. B: Condens. Matter Mater. Phys. 2012, 86, 205110. (38) Refaely-Abramson, S.; Sharifzadeh, S.; Jain, M.; Baer, R.; Neaton, J. B.; Kronik, L. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 88, 081204. (39) Srebro, M.; Autschbach, J. J. Chem. Theory Comput. 2012, 8, 245−256. (40) Sutton, C.; Körzdörfer, T.; Coropceanu, V.; Brédas, J.-L. J. Phys. Chem. C 2014, 118, 3925−3934. (41) Egger, D. A.; Weissman, S.; Refaely-Abramson, S.; Sharifzadeh, S.; Dauth, M.; Baer, R.; Kümmel, S.; Neaton, J. B.; Zojer, E.; Kronik, L. J. Chem. Theory Comput. 2014, 10, 1934−1952. (42) Sutton, C.; Körzdörfer, T.; Gray, M. T.; Brunsfeld, M.; Parrish, R. M.; Sherrill, C. D.; Sears, J. S.; Brédas, J.-L. J. Chem. Phys. 2014, 140, 054310. (43) Zhang, C.-R.; Coropceanu, V.; Sears, J. S.; Brédas, J.-L. J. Phys. Chem. C 2014, 118, 154−158. (44) Körzdörfer, T.; Brédas, J.-L. Acc. Chem. Res. 2014, 47, 3284− 3291. (45) Körzdörfer, T.; Parrish, R. M.; Sears, J. S.; Sherrill, C. D.; Brédas, J.-L. J. Chem. Phys. 2012, 137, 124305. (46) Sun, H.; Zhang, S.; Sun, Z. Phys. Chem. Chem. Phys. 2015, 17, 4337−4345. (47) Foster, M. E.; Wong, B. M. J. Chem. Theory Comput. 2012, 8, 2682−2687. (48) Agrawal, P.; Tkatchenko, A.; Kronik, L. J. Chem. Theory Comput. 2013, 9, 3473−3478. (49) Karolewski, A.; Kronik, L.; Kümmel, S. J. Chem. Phys. 2013, 138, 204115. (50) Autschbach, J.; Srebro, M. Acc. Chem. Res. 2014, 47, 2592−2602. (51) de Queiroz, T. B.; Kümmel, S. J. Chem. Phys. 2014, 141, 084303. (52) Sun, H.; Autschbach, J. J. Chem. Theory Comput. 2014, 10, 1035−1047. (53) Wang, X.-B.; Wang, L.-S. Annu. Rev. Phys. Chem. 2009, 60, 105− 126. (54) Thomas, T. D.; Saethre, L. J.; Sorensen, S. L.; Svensson, S. J. Chem. Phys. 1998, 109, 1041−1051. (55) Malagoli, M.; Coropceanu, V.; da Silva Filho, D. A.; Brédas, J. L. J. Chem. Phys. 2004, 120, 7490−7496. (56) Gallandi, L.; Körzdörfer, T. J. Chem. Theory Comput. 2015, 11, 5391−5400. (57) Refaely-Abramson, S.; Sharifzadeh, S.; Govind, N.; Autschbach, J.; Neaton, J. B.; Baer, R.; Kronik, L. Phys. Rev. Lett. 2012, 109, 226405. (58) Knight, J.; Wang, X.; Gallandi, L.; Ren, X.; Rinke, P.; Körzdörfer, T.; Marom, N.; Dolgounitcheva, O.; Ortiz, J. V. Accurate Ionization Potentials and Electron Affinities of Acceptor Molecules for Organic Photovoltaics III: Benchmarking the Accuracy of GW Methods. J. Chem. Theory Comput. 2016, DOI: 10.1021/acs.jctc.5b00871. (59) Rinke, P.; Qteish, A.; Neugebauer, J.; Freysoldt, C.; Scheffler, M. New J. Phys. 2005, 7, 126. (60) Fuchs, F.; Furthmüller, J.; Bechstedt, F.; Shishkin, M.; Kresse, G. Phys. Rev. B: Condens. Matter Mater. Phys. 2007, 76, 115109. (61) Marom, N.; Moussa, J. E.; Ren, X.; Tkatchenko, A.; Chelikowsky, J. R. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 84, 245115. (62) Faber, C.; Attaccalite, C.; Olevano, V.; Runge, E.; Blase, X. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 83, 115123. (63) Marom, N.; Caruso, F.; Ren, X.; Hofmann, O. T.; Körzdörfer, T.; Chelikowsky, J. R.; Rubio, A.; Scheffler, M.; Rinke, P. Phys. Rev. B: Condens. Matter Mater. Phys. 2012, 86, 245127. (64) Caruso, F.; Rinke, P.; Ren, X.; Scheffler, M.; Rubio, A. Phys. Rev. B: Condens. Matter Mater. Phys. 2012, 86, 081102. (65) Körzdörfer, T.; Marom, N. Phys. Rev. B: Condens. Matter Mater. Phys. 2012, 86, 041110. (66) Atalla, V.; Yoon, M.; Caruso, F.; Rinke, P.; Scheffler, M. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 88, 165122. (67) Bruneval, F.; Marques, M. A. L. J. Chem. Theory Comput. 2013, 9, 324−329.

REFERENCES

(1) Koopmans, T. Physica 1934, 1, 104. (2) Janak, J. F. Phys. Rev. B: Condens. Matter Mater. Phys. 1978, 18, 7165. (3) Perdew, J. P.; Parr, R. G.; Levy, M.; Balduz, J. L., Jr. Phys. Rev. Lett. 1982, 49, 1691. (4) Chong, D. P.; Gritsenko, O. V.; Baerends, E. J. J. Chem. Phys. 2002, 116, 1760−1772. (5) Hedin, L. Phys. Rev. 1965, 139, A796. (6) Hybertsen, M. S.; Louie, S. G. Phys. Rev. B: Condens. Matter Mater. Phys. 1986, 34, 5390. (7) Nooijen, M.; Bartlett, R. J. J. Chem. Phys. 1995, 102, 3629−3647. (8) Hohenberg, P.; Kohn, W. Phys. Rev. 1964, 136, B864. (9) Kohn, W.; Sham, L. Phys. Rev. 1965, 140, A1133. (10) Jain, A.; Ong, S. P.; Hautier, G.; Chen, W.; Richards, W. D.; Dacek, S.; Cholia, S.; Gunter, D.; Skinner, D.; Ceder, G.; Persson, K. A. APL Mater. 2013, 1, 011002. (11) Fischer, C. C.; Tibbetts, K. J.; Morgan, D.; Ceder, G. Nat. Mater. 2006, 5, 641−646. (12) Curtarolo, S.; Hart, G. L. W.; Nardelli, M. B.; Mingo, N.; Sanvito, S.; Levy, O. Nat. Mater. 2013, 12, 191−201. (13) Seidl, A.; Görling, A.; Vogl, P.; Majewski, J. A.; Levy, M. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 53, 3764−3774. (14) Perdew, J. P.; Zunger, A. Phys. Rev. B: Condens. Matter Mater. Phys. 1981, 23, 5048. (15) Ruzsinszky, A.; Perdew, J. P.; Csonka, G. I.; Vydrov, O. A.; Scuseria, G. E. J. Chem. Phys. 2006, 125, 194112. (16) Mori-Sánchez, P.; Cohen, A. J.; Yang, W. J. Chem. Phys. 2006, 125, 201102. (17) Patchkovskii, S.; Autschbach, J.; Ziegler, T. J. Chem. Phys. 2001, 115, 26−42. (18) Körzdörfer, T.; Kümmel, S.; Mundt, M. J. Chem. Phys. 2008, 129, 014110. (19) Kümmel, S.; Kronik, L. Rev. Mod. Phys. 2008, 80, 3−60. (20) Körzdörfer, T.; Kümmel, S.; Marom, N.; Kronik, L. Phys. Rev. B: Condens. Matter Mater. Phys. 2009, 79, 201205. (21) Körzdörfer, T.; Mundt, M.; Kümmel, S. Phys. Rev. Lett. 2008, 100, 133004. (22) Messud, J.; Dinh, P. M.; Reinhard, P.-G.; Suraud, E. Phys. Rev. Lett. 2008, 101, 096404. (23) Pederson, M. R.; Ruzsinszky, A.; Perdew, J. P. J. Chem. Phys. 2014, 140, 121103. (24) Curtiss, L. A.; Redfern, P. C.; Raghavachari, K.; Pople, J. A. J. Chem. Phys. 1998, 109, 42−55. (25) de Oliveira, G.; Martin, J. M. L.; de Proft, F.; Geerlings, P. Phys. Rev. A: At., Mol., Opt. Phys. 1999, 60, 1034−1045. (26) Rienstra-Kiracofe, J. C.; Tschumper, G. S.; Schaefer, H. F.; Nandi, S.; Ellison, G. B. Chem. Rev. 2002, 102, 231−282. (27) Stein, T.; Kronik, L.; Baer, R. J. Am. Chem. Soc. 2009, 131, 2818−2820. (28) Stein, T.; Kronik, L.; Baer, R. J. Chem. Phys. 2009, 131, 244119. (29) Baer, R.; Livshits, E.; Salzner, U. Annu. Rev. Phys. Chem. 2010, 61, 85−109. (30) Körzdörfer, T.; Sears, J. S.; Sutton, C.; Brédas, J.-L. J. Chem. Phys. 2011, 135, 204107. (31) Refaely-Abramson, S.; Baer, R.; Kronik, L. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 84, 075144. (32) Karolewski, A.; Stein, T.; Baer, R.; Kümmel, S. J. Chem. Phys. 2011, 134, 151101. (33) Stein, T.; Eisenberg, H.; Kronik, L.; Baer, R. Phys. Rev. Lett. 2010, 105, 266802. (34) Kronik, L.; Stein, T.; Refaely-Abramson, S.; Baer, R. J. Chem. Theory Comput. 2012, 8, 1515−1531. (35) Moore, B., II; Autschbach, J. J. Chem. Theory Comput. 2013, 9, 4991−5003. (36) Sears, J. S.; Körzdörfer, T.; Zhang, J.-L.; Brédas, C.-R. J. Chem. Phys. 2011, 135, 151103. 613

DOI: 10.1021/acs.jctc.5b00873 J. Chem. Theory Comput. 2016, 12, 605−614

Article

Journal of Chemical Theory and Computation (68) Richard, R. M.; Marshall, M. S.; Dolgounitcheva, O.; Ortiz, J. V.; Brédas, J.-L.; Marom, N.; Sherrill, C. D. Accurate Ionization Potentials and Electron Affnities of Acceptor Molecules I. CCSD(T) Reference Data. J. Chem. Theory Comput. 2016, DOI: 10.1021/acs.jctc.5b00875. (69) Dolgounitcheva, O.; Díaz-Tinoco, M.; Zakrzewski, V. G.; Richard, R. M.; Marom, N.; Sherrill, C. D.; Ortiz, J. V. Accurate Ionization Potentials and Electron Affinities of Acceptor Molecules IV: Electron-Propagator Methods. J. Chem. Theory Comput. 2016, DOI: 10.1021/acs.jctc.5b00872. (70) Vydrov, O. A.; Heyd, J.; Krukau, A. V.; Scuseria, G. E. J. Chem. Phys. 2006, 125, 074106. (71) Vydrov, O. A.; Scuseria, G. E. J. Chem. Phys. 2006, 125, 234109. (72) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865. (73) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1997, 78, 1396. (74) Adamo, C.; Barone, V. J. Chem. Phys. 1999, 110, 6158. (75) Heyd, J.; Scuseria, G. E.; Ernzerhof, M. J. Chem. Phys. 2003, 118, 8207−8215. (76) Heyd, J.; Scuseria, G. E. J. Chem. Phys. 2004, 120, 7274−7280. (77) Blum, V.; Gehrke, R.; Hanke, F.; Havu, P.; Havu, V.; Ren, X.; Reuter, K.; Scheffler, M. Comput. Phys. Commun. 2009, 180, 2175− 2196. (78) Havu, V.; Blum, V.; Havu, P.; Scheffler, M. J. Comput. Phys. 2009, 228, 8367−8379. (79) Ren, X.; Rinke, P.; Blum, V.; Wieferink, J.; Tkatchenko, A.; Sanfilippo, A.; Reuter, K.; Scheffler, M. New J. Phys. 2012, 14, 053020. (80) Turney, J. M.; Simmonett, A. C.; Parrish, R. M.; Hohenstein, E. G.; Evangelista, F. A.; Fermann, J. T.; Mintz, B. J.; Burns, L. A.; Wilke, J. J.; Abrams, M. L.; Russ, N. J.; Leininger, M. L.; Janssen, C. L.; Seidl, E. T.; Allen, W. D.; Schaefer, H. F.; King, R. A.; Valeev, E. F.; Sherrill, C. D.; Crawford, T. D. Wiley Interdisciplinary Reviews: Computational Molecular Science 2012, 2, 556−565. (81) Supporting Information. (82) Henderson, T. M.; Izmaylov, A. F.; Scalmani, G.; Scuseria, G. E. J. Chem. Phys. 2009, 131, 044108. (83) Delaney, K.; García-González, P.; Rubio, A.; Rinke, P.; Godby, R. W. Phys. Rev. Lett. 2004, 93, 249701. (84) Milian, B.; Pou-AméRigo, R.; Viruela, R.; Orti, E. Chem. Phys. Lett. 2004, 391, 148−151. (85) Ren, X.; Marom, N.; Caruso, F.; Scheffler, M.; Rinke, P. Phys. Rev. B: Condens. Matter Mater. Phys. 2015, 92, 081104. (86) Stein, T.; Autschbach, J.; Govind, N.; Kronik, L.; Baer, R. J. Phys. Chem. Lett. 2012, 3, 3740−3744. (87) Körzdörfer, T. J. Chem. Phys. 2011, 134, 094111. (88) Mori-Sánchez, P.; Cohen, A. J.; Yang, W. Phys. Rev. Lett. 2008, 100, 146401. (89) Vlček, V.; Eisenberg, H. R.; Steinle-Neumann, G.; Kronik, L.; Baer, R. J. Chem. Phys. 2015, 142, 034107. (90) Lany, S.; Zunger, A. Phys. Rev. B: Condens. Matter Mater. Phys. 2008, 78, 235104. (91) Ö ğüt, S.; Chelikowsky, J. R.; Louie, S. G. Phys. Rev. Lett. 1997, 79, 1770−1773.

614

DOI: 10.1021/acs.jctc.5b00873 J. Chem. Theory Comput. 2016, 12, 605−614