Outer-valence Electron Spectra of Prototypical Aromatic Heterocycles from an Optimally Tuned Range-Separated Hybrid Functional.

Density functional theory with optimally tuned range-separated hybrid (OT-RSH) functionals has been recently suggested [Refaely-Abramson et al. Phys. Rev. Lett.2012, 109, 226405] as a nonempirical approach to predict the outer-valence electronic structure of molecules with the same accuracy as many-body perturbation theory. Here, we provide a quantitative evaluation of the OT-RSH approach by examining its performance in predicting the outer-valence electron spectra of several prototypical gas-phase molecules, from aromatic rings (benzene, pyridine, and pyrimidine) to more complex organic systems (terpyrimidinethiol and copper phthalocyanine). For a range up to several electronvolts away from the frontier orbital energies, we find that the outer-valence electronic structure obtained from the OT-RSH method agrees very well (typically within ∼0.1-0.2 eV) with both experimental photoemission and theoretical many-body perturbation theory data in the GW approximation. In particular, we find that with new strategies for an optimal choice of the short-range fraction of Fock exchange, the OT-RSH approach offers a balanced description of localized and delocalized states. We discuss in detail the sole exception found-a high-symmetry orbital, particular to small aromatic rings, which is relatively deep inside the valence state manifold. Overall, the OT-RSH method is an accurate DFT-based method for outer-valence electronic structure prediction for such systems and is of essentially the same level of accuracy as contemporary GW approaches, at a reduced computational cost.

Introduction

Many electronic properties of molecules and materials are determined by and understood through the energetics of the valence electrons, which are often probed experimentally using photoemission spectroscopy (PES).[1] Via measurement of the kinetic energy of photoemitted electrons, PES provides direct experimental access to the electron ejection energies, the smallest of which is the ionization potential (IP). Calculating PES data from first principles is a long-standing challenge to modern electronic structure methods.[2,3] A state-of-the-art method for obtaining ionization spectra theoretically is many-body perturbation theory (MBPT), which calculates quasi-particle excitation energies via solving the Dyson equation, typically within the GW approximation (where G is the Green function and W is the dynamically screened Coulomb potential).[4−7] However, GW calculations can be computationally demanding. Moreover, in particular for gas-phase computations that are at the center of this work, results can be sensitive to details of the specific GW scheme employed[8−13] and can also be challenging to converge.[11,14]

Density functional theory (DFT),[15,16] in which the ground-state electron density, rather than the many-electron wave function, is the fundamental quantity,[17] is a computationally efficient first principles method for calculating the electronic structure of many-electron systems. DFT is often made practical via the Kohn–Sham (KS) scheme,[18] in which the original many-electron problem is mapped uniquely into a fictitious noninteracting electron system yielding the same electron density. This mapping leads to effective single-particle equations that provide a significant conceptual and computational simplification of the original many-electron problem. However, due to the fictitious nature of the noninteracting electrons, the correspondence of KS eigenvalues with ionization energies measured in an experiment is not at all straightforward.[3,19] It can be shown that for the exact KS potential, the energy of the highest occupied molecular orbital (HOMO) equals the negative of the IP, a result known as the IP theorem.[20−23] Lower-lying eigenvalues do not strictly correspond to electron removal energies.[19] For outer-valence electrons, however, exact DFT eigenvalues may still serve as a useful and even quantitative approximation to electron removal energies.[3,24] We note that in general the simulation of photoemission spectra also requires that the photoionization cross-section is addressed.[1] Its calculation is outside the scope of the present manuscript, which focuses on the correct description of the energetics. Our comparison between the calculated density of states and the photoelectron spectra concentrates, therefore, on peak positions rather than peak intensities. We note, however, that for gas-phase ultraviolet photoemission spectroscopy (UPS)—the standard choice for probing molecular outer-valence states—angle-dependent cross-section effects are irrelevant, owing to orientational averaging, and orbital-dependent cross-section effects are relatively weak.[25] Therefore, their neglect is not expected to have serious consequences.

The KS mapping scheme relies on an exchange-correlation energy-functional of the density, the exact form of which is generally unknown and must be approximated. Common approximate exchange-correlation functionals used within the KS scheme are the local density approximation (LDA) and the generalized gradient approximation (GGA).[15,16] In the former, one assumes that at each point in space the exchange-correlation energy per particle is given by its value for a homogeneous electron gas. In the latter, information on deviations from homogeneity is partly accounted for by considering gradients of the charge density. Unfortunately, for gas-phase molecules, KS eigenvalues obtained through the use of these approximations are not in good agreement with experimental ionization spectra.[3,26] First, the IP theorem is grossly disobeyed and the negative of the HOMO energy usually underestimates the IP severely.[27−30] Even if this difference is accounted for by rigidly shifting the theoretical eigenvalue spectrum,[31,32] the calculated eigenvalue spectrum may still exhibit qualitative failures, notably an erroneous ordering of the electronic levels (see, e.g., refs (33)–[41]).

These two drawbacks can be traced back to two different (yet not unrelated) deficiencies of the above-described approximations: lack of a derivative discontinuity (DD) and orbital self-interaction errors (SIE).[3,26,28,42] The DD is a uniform “jump” in the KS potential, when approaching the integer electron number either from above or from below. This “jump” helps to account for the discontinuity of the chemical potential at integer electron numbers, i.e., for the fact that the electron removal energy is not the same as the electron insertion energy.[16] Part of the discontinuity in the chemical potential is accounted for by the KS kinetic energy term.[20,23] But the kinetic energy contribution to the discontinuity in the chemical potential is generally insufficient, and the remaining discontinuity must be provided by the DD. Because the electron–ion and Hartree energies are continuous in the electron density, the remaining discontinuity can only be attained by a jump in the exchange-correlation potential. However, in standard LDA or GGA functionals, the exchange-correlation potential explicitly depends on the electron density and does not incorporate any orbital dependence (in contrast to, e.g., Fock exchange). As a result, calculations based on these approximations cannot exhibit any DD in the exchange-correlation part of the potential.[43] Instead, they approximately average over it, and as a consequence the KS-HOMO energy is strongly underbinding with respect to the “true” ionization potential.[31,44,45]

The SIE[46] arises from the fact that the classical electron–electron repulsion term (Hartree potential) in the KS equation means that each electron is repelled from the total charge in the system, including a spurious repulsion from itself. Because KS theory is, in principle, exact, whatever error one makes in the Hartree term must then be completely canceled out by the exact exchange-correlation term. Unfortunately, only partial error cancellation is obtained in either LDA or GGA. For strongly localized orbitals, self-interaction may be significant and spuriously destabilize electron energies.[26]

The situation can be improved, if one uses “hybrid” DFT, i.e., functionals that contain exact exchange based expressions employed using a nonlocal Fock operator.[26] We emphasize that while such functionals do not fall within the KS scheme, they are still very much within the DFT framework[3,26,47] through the generalized Kohn–Sham (GKS) scheme.[48] In the KS scheme, many-body effects are incorporated entirely in a multiplicative potential (which is the sum of the Hartree and exchange-correlation potentials). In contrast, in the GKS scheme, many-body effects are incorporated in a combination of a multiplicative potential and a nonlocal operator (nonmultiplicative potential). This is achieved via mapping to a partially interacting electron gas. Generally, in the GKS scheme, the additional nonlocal operator can mitigate the need for a DD in the multiplicative potential.[26] This may reduce errors associated with averaging over the DD in practical approximations. Most hybrid functionals in everyday use are of the global type, i.e., they contain a fixed fraction of Fock exchange.[49−51] In practice, one typically observes that HOMO energies extracted from such hybrid functionals are closer to experimental ionization energies than those obtained from LDA or GGA but still significantly underestimate those observables.[9,32,52] Spectral distortions (including the possibility of an erroneous ordering of the electronic levels) are often mitigated, as the fraction of Fock exchange reduces the SIE.[3,26,33,34] However, the quantitative details of the eigenvalue spectrum typically still depend on the specific choice of the approximate hybrid functional.[53]

One reason for the failure of conventional hybrid functionals to obey the IP theorem is the presence of only a fraction of exact exchange. Because of this, they do not yield the correct ∼1/r asymptotic potential that should be “felt” by an electron at large distances from the molecule, which is especially relevant for describing the ionization process. However, when using the full Fock exchange to correct for that, the delicate balance between exchange and correlation is disrupted, which is highly detrimental—especially for short-range electron–electron interactions that govern chemical bonding.[26] A promising strategy for tackling that problem is offered by the more recent class of range-separated hybrid (RSH) DFT functionals,[54] pioneered by Savin and co-workers.[55] In these functionals, the interelectron Coulomb repulsion term is separated into long-range (LR) and short-range (SR) components via a range-separation parameter γ. The LR term is mapped using full Fock exchange, thereby establishing the correct asymptotic potential. The SR termed is (typically) mapped using a GGA approach, maintaining the compatibility between the exchange and correlation expressions. In this approach, one still needs to determine the range-separation parameter, γ.

Both formal considerations[56] and practical simulations show that aiming at accurate results for a broad range of systems, one typically needs significantly different values of γ. This is taken care of by using optimally tuned RSH (OT-RSH) functionals.[32,47] There, instead of using one and the same range-separation parameter for all systems, γ is tuned for each system such that physically motivated tuning conditions are fulfilled without introducing any empirical parameters. In particular, it has been shown for gas phase molecules that insisting on the HOMO energy being equal to the negative of the IP (i.e., fulfilling the IP theorem) for the neutral and anion species yields highly accurate HOMO energies and HOMO–LUMO gaps when compared to experimentally measured fundamental gaps or to the results of GW calculations.[52,57] More recently, some of us have shown that in addition to the IP and the fundamental gap, the entire higher-lying part of the valence-electron spectrum of gas-phase molecules and molecular solids can be accurately described by the eigenvalues of OT-RSH ground-state DFT calculations.[58,59] In particular, it has been suggested that a more general OT-RSH functional that introduces a fraction of Fock exchange in the SR and simultaneously maintains the full Fock exchange in the LR[60] allows for a more flexible treatment of differently localized molecular orbitals, resulting in an accurate description of more complex organic molecules relevant for applications in organic electronics.[58]

In the present contribution, we further investigate the capabilities of such more general OT-RSH functionals for predicting outer-valence electron spectra of organic molecules in the gas phase. First, we study the prototypical aromatic building blocks—benzene, pyridine, and pyrimidine (see Figure 1a). This choice is motivated by the fact that for such simple systems, existing high-level experimental PES data can serve as useful benchmarks for theory. Moreover, the nitrogen heteroatoms in the azabenzenes can be expected to result in differently localized molecular orbitals in the higher-lying part of the valence electron system and, in particular, close-lying σ and π states.[30] For their accurate description, the OT-RSH approach would need to attain a quantitatively satisfactory balance of self-interaction errors for both.[58] Moreover, these prototypical systems were recently identified as a challenge for both the GW[61] and the OT-RSH[62] methods. Here, we perform both OT-RSH and GW calculations for these systems and find that they overall yield similar (to ∼0.2 eV) eigenvalue energies, both being highly accurate compared to PES experiments. We suggest that a simultaneous reliable prediction of both π and σ orbital energies is indeed within the realm of the OT-RSH functional applied here. Still, we identify one specific molecular orbital that is peculiar to ring-type molecules, which in the OT-RSH calculations displays a significant deviation from experimental results and GW calculations. We analyze the origin of this discrepancy by further computing the spectra of the same systems using conventional hybrid calculations, as well as explicitly self-interaction corrected (SIC) calculations.

Figure 1

Molecules studied in this article. (a) Prototypical aromatic rings—benzene, pyridine, and pyrimidine. (b) More complex representative systems—terpyrimidinethiol and copper phthalocyanine.

With the obtained overall very encouraging results at hand, we proceed toward larger and more complex systems, here chosen to be terpyrimidinethiol and copper phthalocyanine (3N-thiol and CuPc, see Figure 1b). These molecules, which also contain N as a heteroatom, are interesting for novel applications in organic electronics[63−68] but at the same time challenging to assess theoretically due to pronounced differences in SIE among the high-lying orbitals in the valence electron spectrum.[34,68,69] Through a comparison to GW calculations and/or PES experiments, we show that OT-RSH can provide accurate valence-electron spectra also for these more complex organic systems, with an optimal choice for the short-range Fock exchange that is guided by conventional hybrid functional calculations. Our results clearly demonstrate that OT-RSH functionals are a highly promising, state-of-the-art approach for predicting ionization spectra of gas-phase organic molecules.

Theoretical and Methodological Details

Formalism

As mentioned in the Introduction, we examine a generalized RSH form, which allows for different amounts of Fock exchange in the short range and in the long range.[70] We use the range-partitioning expression of Yanai et al., given by[71]Here, r is the interelectron coordinate and α, β, and γ are adjustable parameters. Naturally, this partition is not unique, but the choice of the error function is computationally convenient when using a Gaussian basis for expanding the wave functions of finite systems.[55] Equation 1 defines the range-separation procedure, where the Coulomb operator, 1/r, in the exchange-part of the xc potential is replaced by two complementary terms, which are treated differently. As suggested by the “HF” and “GGA” delimiters in the equation, the first term is treated using Hartree–Fock exchange, the second term is treated using GGA semilocal exchange. Specifically, using the Perdew–Burke–Ernzerhof (PBE)[72] form for the GGA exchange-correlation leads to the following expression for the exchange-correlation energy[73]where the superscripts LR and SR denote that the full Coulomb repulsion, 1/r, has been substituted by the LR and SR Coulomb repulsions, erf(γr)/r and erfc(γr)/r, respectively. PBE correlation is used for the entire range. Yanai et al. viewed α, β, and γ as semiempirical parameters and determined “universal” values for them based on benchmark thermochemistry data. Here, we wish to determine these parameters based on the satisfaction of inherent constraints on the exchange-correlation density functional, without recourse to experimental data.

From eq 2 it is clear that for full long-range Fock exchange, which guarantees the correct asymptotic potential, one condition is that α + β = 1 (note that the fit to thermochemistry data performed by Yanai et al. leads to α + β = 0.65, i.e., to a potential that is not asymptotically correct and will, therefore, run into problems when trying to associate eigenvalues with ionization energies or electron affinities). Equation 2 further shows that α controls the fraction of SR Fock exchange. We restrict the current investigation to α values between 0 and 0.2. This is because in the former limit, the SR behavior is expected to be GGA-like; in the latter limit, the SR behavior is expected to resemble that of a conventional hybrid functional.[58,74] For example, in the PBE0 hybrid functional, which is based on the PBE semilocal functional, a global exact-exchange fraction of 25% is used.[51] In a RSH functional, the SR Fock exchange fraction is indeed expected to be somewhat smaller than in a conventional hybrid functional, because the LR Fock exchange is 100% in the former but only a finite fraction in the latter. The remaining parameter, γ, controls the range-separation. As mentioned in the Introduction, we do not seek a universal value for γ. Instead, we rely on a nonempirical tuning procedure, where γ is adjusted on a per-system basis. One possibility to determine the optimal gamma, γopt, is to choose it such that the IP theorem is satisfied:[32,75]εH(N) is the HOMO energy of the N-electron system and IP(N) is the ionization potential of the N-electron system, determined from total energy difference of the N- and (N – 1)-electron systems. In general, as emphasized by the superscript notation, both εH(N) and IP(N) display a strong γ dependence.

It is often useful to invoke the IP theorem not only for the molecule in its neutral state but also in certain charged states. The condition which then needs to be fulfilled is that the target function J2, given by[32]is minimized. In eq 4, i can in principle adopt any integer number, and one can observe the effect of adding further terms on the residual J value (vide infra). In practice, only values of i that are close to zero are of interest to avoid information from highly charged radical species—a point elaborated below. In particular, including the anion has been found to be highly beneficial for the prediction of fundamental gaps and the energies of charge-transfer excited states.[52,76,77] Note that when restricting i to 0 in eq 4, which we do here if the electron affinity is negative and the molecule does not bind an extra electron, then eq 3 is recovered. In agreement with previous work,[32,52,57−60,78−80] it was repeatedly found that the optimal γ strongly changes from one system to the other, showing that our “per system” tuning approach is indeed necessary.

The above-mentioned optimal-tuning of γ can, in principle, be performed for any choice of the SR Fock exchange parameter, α. In fact, various strategies can be employed to determine α. These are discussed below in detail, along with their pros and cons, in the context of specific computational results.

Computational Details

All PBE calculations (within the KS scheme) and PBE0 and RSH calculations (within the GKS scheme) presented in this article were obtained within the QChem[81] and NWChem[82] codes, using cc-PVTZ[83] basis functions. All geometries were optimized using the PBE functional. As the above-described tuning scheme is based on components taken from well-established density functionals (cf. eq 2) and range-separation is available in many different electronic structure codes, optimal tuning of γ via eq 4 is straightforward to perform. Moreover, the tuning strategy is computationally efficient, as it relies on a series of standard total-energy DFT calculations. We note, however, that when charging the gas-phase molecule as part of the above-explained tuning procedure, the configuration of the cation and/or anion must be identified carefully, as it may affect the results of the calculation.[84]

We illustrate the possible complications associated with different ion configurations by considering pyridine—one of the molecules of Figure 1, which is discussed extensively in the results section. When tuning γ using the above-described approach, it is important to ensure that the character of the HOMO (related to the left-hand term in eq 3) corresponds to the “hole density,” defined as the charge-density difference of the neutral and charged states (and, consequently, related to the right-hand term in eq 3). For pyridine, the self-consistent solution of the GKS equation with the RSH functional was found to lead to two different doublet configurations of the cation, depending on the initial guess used in the procedure. These two configurations correspond to two qualitatively different hole densities (see Figure 2, left part). As shown in the right part of Figure 2, the reason for the different hole densities is that the two cation configurations possess two different LUMO orbitals; i.e., the two cationic ground states represent two different ionization processes. The main difference is that the electron “loss” is, in one case, from a π orbital and, in the other case, from a σ orbital. These two cationic configurations are energetically close, which is consistent with the observation that the HOMO and the HOMO–1 of pyridine are very close in energy (vide infra). Importantly, however, only the “hole density” depicted in Figure 2b—which is associated with the configuration lower in energy, i.e., the true ground state predicted for the cation—corresponds to the HOMO of pyridine (see Figure 2c). Therefore, one has to ensure that the cationic state shown in Figure 2b is indeed the one entering the tuning procedure, in order to retain consistency for the orbital energies and total energies required in eq 3.

Figure 2

Charge-density difference between neutral and cation (left) and LUMO of cation (right), obtained using two different doublet configurations (a and b) of the pyridine cation. The configuration denoted by b is the energetically more stable one. In the charge-density difference plots, red (blue) regions denote areas of electron density depletion (accumulation) as a consequence of the ionization process. (c) HOMO of the neutral pyridine molecule.

For our analysis of the OT-RSH results, we also performed comparative GW calculations, as well as self-interaction-corrected calculations and KS-PBE0 calculations (the latter are defined and explained below). Our GW calculations are based on a standard G0W0 scheme,[14] where quasi-particle energies are computed via a first-order correction to DFT eigenvalues, with no self-consistent update of the starting wave functions. The starting quasi-particle wave function for the G0W0 corrections was obtained from the PBE functional.[72] The static dielectric function is computed within the random-phase approximation and extended to finite frequency via the generalized plasmon-pole (GPP) model of Hybertsen and Louie.[5]

Our G0W0 calculations were performed using the BerkeleyGW package,[85] which employs a plane-wave basis set to compute the dielectric function and self-energy, using a PBE starting point. DFT-PBE calculations were performed within the Quantum Espresso package,[86] which is compatible with BerkeleyGW. The nuclei and core electrons were described by Troullier–Martins relativistic norm-conserving pseudopotentials,[87] which are part of the Quantum Espresso pseudopotential library. Here, one, four, five, and six electrons were explicitly considered as valence electrons for H, C, N, and S, respectively, with cutoff radii (in Bohr) of 1.3, 0.5, 1.0, and 1.7, respectively. We used a plane-wave basis cutoff of 80 Ry for benzene and 120 Ry for pyridine, pyrimidine, and 3N-thiol. These values lead to a total DFT energy converged to <1 meV/atom. To avoid spurious interactions with periodic images, we used a supercell with lattice vectors twice the size necessary to contain 99% of the charge density and, when computing the GW self-energy, the Coulomb interaction was truncated at distances larger than half of the unit cell size. The supercell dimensions, in atomic units, were 35 × 39 × 24; 30 × 20 × 32, 19 × 30 × 30; and 64 × 26 × 15.5 for benzene, pyridine, pyrimidine, and 3N-thiol, respectively.

Our static dielectric function and self-energy were constructed from 4914, 5515, 5071, and 3598 unoccupied states, respectively, for benzene, pyridine, pyrimidine, and 3N-thiol. For the former three prototypical small molecules, this energy range corresponds to 90 eV above the vacuum energy, while for 3N-thiol, it corresponds to 50 eV above the vacuum energy. Fewer states were included for 3N-thiol due to the greater computational expense associated with this rather large system. A static remainder approach was applied to the self-energy to approximately complete the unoccupied subspace.[88] The plane-wave cutoff for the dielectric function was 30 Ry for pyridine and pyrimidine and 24 Ry for 3N-thiol and benzene. We find that these parameters converge the HOMO energies of the prototypical small molecules to less than 0.1 eV. On the basis of the convergence behavior of these molecules and the residual differences that we find for GW and OT-RSH HOMO energies (vide infra), we extrapolate the errors associated with eigenvalues of the 3N-thiol calculation to be less than 0.2 eV.

All SIC calculations were based on the seminal SIC concept of Perdew and Zunger.[46] However, we constructed a spatially local, multiplicative exchange-correlation potential identical for all orbitals in the system, which ensures that the SIC remains within the KS realm.[89,90] This is based on the generalized optimized effective potential (OEP) equation, which extends the original OEP equation to the case of unitarily variant functionals. It is solved using the generalized Krieger–Li–Iafrate (KLI) approximation.[89,91] Unlike the KLI approximation to the standard OEP equation, which can introduce significant deviations for the SIC,[92] the generalized KLI approximation used here has been shown to be an excellent approximation to the generalized OEP. The additional degree of freedom arising from the variance inherent to our procedure can be used to construct a set of orbitals that minimize the total SIC energy of the system, where we applied a complex-valued energy minimizing unitary orbital transformation.[90] For additional insights, we also used the PBE0 functional in conjunction with a local multiplicative potential, constructed—in contrast to traditional GKS schemes—via the KLI[91] approximation for the exact exchange part of the functional. We refer to these calculations as PBE0KS.

All SIC and PBE0KS calculations were performed with the Bayreuth version[93] of the PARSEC real-space code,[94] where we employed a grid-spacing of 0.2 Bohr and Troullier–Martins norm-conserving pseudopotentials.[87]

Finally, for a meaningful comparison between the results of different functionals and/or computational approaches, we tested explicitly that eigenvalues obtained from different codes and basis-set expansions (Gaussian, planewave, real-space) do not differ by more than 0.1 eV for the same underlying functional. Furthermore, we verified by visual inspection[95,96] that eigenvalues calculated from different methods correspond to the same molecular orbitals.

Prototypical Aromatic Rings—Results and Discussion

We start our analysis by considering the prototypical aromatic gas-phase molecules of Figure 1a. The computed peak positions in the eigenvalue spectra (neglecting differences of the photoionization cross sections) should correspond to those in the photoemission spectra (in arbitrary units, denoted as a.u. or arb. units) and are shown in Figure 3. The computed spectra are obtained from OT-RSH eigenvalues with SR Fock exchange fractions of α = 0 and α = 0.2 for (a) benzene, (b) pyridine, and (c) pyrimidine. The results are compared to gas-phase PES spectra from refs (97)–[99] and to the results of our GW calculations. Importantly, the spectra are plotted on an absolute energy scale; i.e., the calculated eigenvalue spectra shown in Figure 3 have not been shifted so as to align the theoretical data with experimental results. To facilitate a visual comparison with the experimental spectra, we broadened the calculated eigenvalue spectra via convolution with a Gaussian (with a standard deviation of 0.4 eV for benzene and 0.2 eV for pyridine and pyrimidine). For all three systems, we find an overall excellent agreement between OT-RSH, our GW calculations, and experimental results. In particular, the HOMO energies calculated from the OT-RSH approach correspond very well to the first peak in the experimental spectrum, indicating the success of the tuning procedure. In addition, the shape of the experimental spectra at higher binding energies is very well reproduced by the OT-RSH spectra. Note that the comparison extends over a relatively large energy range of ∼8 eV below the first IP (much more than the ∼3 eV in previous work on somewhat larger aromatic-based organic molecules[58]). The agreement over such a wide energy range is remarkable, considering that the correspondence between DFT eigenvalues and photoemission energies is expected, on general theoretical grounds, to deteriorate for lower-lying states.[3,24,100] Interestingly, the shape of the OT-RSH spectra for these molecules is only mildly sensitive to the choice of α (which is not always the case—see below).

Figure 3

Valence electron spectra for (a) benzene, (b) pyridine, (c) pyrimidine, as obtained from experiment (benzene and pyridine, ref (97); pyrimidine, refs (99) and (98), reference I and II in the figure), compared with simulated data obtained from GW and OT-RSH calculations (with two different amounts of SR Fock exchange, α), broadened with Gaussians of widths 0.4, 0.2, and 0.2 eV, respectively, to facilitate comparison with experimental results.

Despite the overall agreement, there is one particular feature which does not agree with experimental results for each of the three molecules: For benzene, the pronounced intensity in the OT-RSH spectra around −13 eV is in contrast to the very low intensity seen in the experiment. For pyridine and pyrimidine, mismatches occur around −13 eV and −15 eV, respectively. These discrepancies, discussed in more detail below, are observed for both considered SR choices (α = 0 and α = 0.2). The main deviation between our GW calculations and experiments is found for the benzene molecule, where the feature that in the experiment occurs around −11.5 eV is shifted to higher binding energies. A similar discrepancy (as well as some additional ones for the other rings) has previously been pointed out by Marom et al.[61]

A precise quantitative comparison between theory and experiment beyond the above statements is complicated by several aspects: (i) The peaks in the experiments are not well separated; i.e., to determine peak maxima and, correspondingly, vertical ionization energies, one would need to perform extensive fitting. (ii) While important information on orbital symmetry and localization can be obtained directly from experimental results, it is common practice to provide a detailed assignment of orbital shape and ordering based on comparisons to theory, as was done in, e.g., ref (97). (iii) There are also some deviations between different experimental reports, e.g., for pyrimidine, where the two experimental spectra from refs (98) and (99) (shown in Figure 3c) appear to be shifted with respect to each other by ∼0.2 eV. Therefore, in the following, we choose to compare theoretical data only to the original experimental spectra and avoid the extraction of vertical ionization series from the experiments. Nevertheless, in the interest of putting our OT-RSH results for these prototypical systems into perspective with previous literature findings,[62] we adopt the extracted experimental values for the ionization energies reported by Marom et al.[61] We can then perform a straightforward statistical analysis of mean absolute differences (MADs) between theory and experiment. This yields MADs of ∼0.1 eV for our GW calculations and MADs of ∼0.2 eV for the OT-RSH (α = 0.2) calculations. Ref (62) can be understood to suggest that the accuracy of the OT-RSH method falls short of that of GW. In fact, the MADs of the OT-RSH calculations differ from those of our GW calculations by an extent which is marginal and on par with the experimental resolution. Furthermore, the MAD between OT-RSH and GW results is ∼0.2 eV, i.e., in the same range. Similar conclusions hold for the 3,4,9,10-perylene-tetracarboxylic-dianhydride (PTCDA) and 1,4,5,8-naphtalene-tetracarboxylic-dianhydride (NTCDA) molecules discussed in detail in refs (58) and (62).

For a deeper analysis of the remaining differences between theory and experiment and their origin, we turn to a detailed comparison between the OT-RSH and GW eigenvalues, given in Figure 4. For each eigenvalue, the figure also provides the corresponding single-electron orbital.

Figure 4

Valence eigenvalue spectra of (a) benzene, (b) pyridine, and (c) pyrimidine, obtained from GW, OT-RSH, and the nonoptimally tuned LC- ωPBE and LC-ωPBE0 RSH functionals. The GW spectra shown were obtained from literature data (scanned in with corresponding color coding from ref (61)) with two different starting points, and from our own work. The OT-RSH data were obtained from two different choices for the amount of SR Fock exchange, α.

It has been recently implied[62] that the satisfactory agreement between OT-RSH and GW results, reported in ref (58), worsens if the comparison is made to a GW calculation based on a hybrid-functional DFT starting point, rather than on a PBE starting point. To examine this, the first two columns in Figure 4 show PBE- and PBE0-based G0W0 valence-electron energies computed by Marom et al.[61] for the small aromatic systems. Our PBE-based G0W0 calculations, based on the scheme presented in ref (14) and already shown in Figure 3 to agree well with experiment results, are given in the third column in Figure 4. They are followed by our two (α = 0 and α = 0.2) OT-RSH results. As in Figure 3, all eigenvalues are shown on an absolute energy scale, with no shifting of the computed data. For all three prototypical rings, our PBE-based G0W0 calculations, within a plasmon-pole approximation,[5] are closer to the PBE0-based “full frequency”[101] G0W0 calculations of Marom et al.[61] than to their PBE-based ones. Specifically, the MAD between our PBE-based G0W0 results and the literature ones (averaged over outer valence orbitals of all three molecules) is 0.56 eV. The MAD is reduced to only 0.22 eV upon comparison with the G0W0 results based on PBE0. The difference stems from the fact that these two non-self-consistent, “one-shot” G0W0 results are calculated using different approximations within the GW scheme itself. Here, we use a plasmon-pole approximation, pseudopotentials, and a plane-wave basis, an approach that yields good agreement with measured IPs for small gas-phase molecules;[12,14] ref (61) reports all-electron calculations, with a fully frequency-dependent dielectric function and a finite localized basis. Indeed, the outcome of a GW calculation can depend on more than just the starting point,[9,12,14,102,103] and more specifically, plasmon-pole models have been seen, in other classes of systems such as simple oxides,[104,105] to enhance the magnitude of GW corrections. A detailed discussion of the relative accuracy and precision and the pros and cons of different G0W0 approaches for molecular systems is well outside the scope of this manuscript and will be taken up elsewhere. Here, the only salient point is that the similarity between the second and third columns of Figure 4 validates the comparison made here and in ref (58), to the PBE-based GW methodology of ref (14) (simply referred to as GW hereafter). Generally, differences of 0.1–0.2 eV, owing to numerical precision, are expected in present-day GW approaches. As mentioned above, uncertainties of similar magnitude are expected in many of the experiments relevant in the present context. Beyond numerical precision, owing to the physical approximations inherent in any choice of DFT functional, as well as in a G0W0 or GW approximation to the self energy, energy differences smaller than ∼0.1–0.2 eV are beyond the accuracy of these approaches.

Turning to the detailed comparison between OT-RSH and our GW calculation, given in Figure 4, the following picture emerges. First, the OT-RSH HOMO energies for all three systems are in very good (<0.15 eV) agreement with the respective lowest quasi-hole energy computed from GW. Excellent agreement of OT-RSH HOMO values with reference theoretical or experimental IPs has been previously reported for a variety of systems—see, e.g., refs (52), (57), and (78). Nevertheless, this finding is still not trivial, because the HOMO of benzene consists of degenerate π orbitals, whereas the HOMO of pyridine and pyrimidine is of σ character (n, to be precise). We emphasize the importance of the tuning procedure for obtaining this level of accuracy for the HOMO energy. For comparison, consider the LC-ωPBE[106] or LC-ωPBE0[74] functionals. Both are PBE-based RSH functionals just like the one used here, with (α, γ) pairs of (α = 0, γ = 0.4 Bohr–1) and (α = 0.2, γ = 0.2 Bohr–1), respectively; i.e., they employ a fixed nontuned value of γ. The outer-valence energies obtained from these two functionals are shown in the two right-most columns of Figure 4. For the LC-ωPBE functional, the fixed-γ value is much larger than the optimally tuned one at α = 0 (cf. Figure 3; γopt = 0.287 Bohr–1, γopt = 0.312 Bohr–1, and γopt = 0.353 Bohr–1 for benzene, pyridine, and pyrimidine, respectively), and consequently, the LC-ωPBE spectrum is too low in energy. Complementarily, the LC-ωPBE0 fixed-γ value is smaller than that obtained with optimal γ tuning at α = 0.2 (cf. Figure 3; γopt = 0.238 Bohr–1, γopt = 0.235 Bohr–1, and γopt = 0.278 Bohr–1 for benzene, pyridine, and pyrimidine, respectively), and the spectrum is too high in energy. Even for benzene, where this is a smaller effect, the LC-ωPBE0 are shifted by ∼0.35 eV from our GW data and by ∼0.25 eV from the OT-RSH (α = 0.2) data. For pyrimidine, the effect is more pronounced and the shifts amount to ∼0.55 eV and ∼0.6 eV, respectively. This observation is fully consistent with the discussion of ref (32) (especially Figure 6 therein) and underscores once again the importance of optimal tuning.

With the optimally tuned functional, the predicted deeper-lying parts of the spectra (with one exception, mentioned previously and discussed in more detail below) are in overall excellent agreement with our GW calculations in spite of the large energy range and the fact that for all systems both σ and π orbitals are present. Increasing α from 0 to 0.2 slightly improves the agreement between OT-RSH and GW, but the effect is overall minor for this kind of system; shifts in orbital energy amount to up to ∼0.2 eV, as shown in the Supporting Information (SI)). The reasons for this minor impact of α are discussed below.

One may question the accuracy of the OT-RSH functional for describing states with higher binding-energies because of the absence of piecewise linearity[62] (discussed in more detail below) for highly charged species of the studied system. Such piecewise linearity should indeed be obtained with the exact DFT functional. However, it should be kept in mind that a highly charged radical is much less stable and undergoes considerable charge density relaxation and orbital reordering with respect to the neutral or cationic molecule. As shown above for the case of the pyridine molecule, even for the (only singly charged) cation, it can be difficult to obtain a meaningful ground state, and this issue can be expected to be more severe for more highly charged molecules. Therefore, while indeed piecewise linearity in principle should be observed with the exact functional also for highly charged species, we believe that other issues are more relevant for describing accurately the excited states of the singly charged cation, which are the ones essential for describing the photoemission process.[3,93,107,108]

The one discrepancy between the GW and OT-RSH data (vide supra) involves a particular ring-shaped π orbital (orange-colored in Figure 4), which is specific to cyclic compounds. In the OT-RSH, this orbital appears significantly below its GW position (by ∼0.4 eV for benzene and pyridine and an even worse ∼0.7 eV for pyrimidine, almost independently of α). For all three molecules, this orbital is primarily responsible for the remaining disagreement between the OT-RSH calculations and the experimental spectra in Figure 3.

To understand the OT-RSH results in more detail, and to explore the possible origins of this remaining discrepancy, we performed additional DFT calculations with several different functionals. Of all DFT functionals studied here, the OT-RSH one is the only one capable of obeying the ionization potential theorem of eq 3 and, consequently, providing HOMO energies that are close to the experimental IP.[3,32] Therefore, in the following comparison, all energies are reported relative to the HOMO energy of the respective calculation (which is set to zero). The resulting shifted eigenvalue spectra for the three prototypical aromatic molecules are shown in Figure 5, with the original HOMO energies shown on top.

Figure 5

Shifted eigenvalue spectra of (a) benzene, (b) pyridine, and (c) pyrimidine obtained from different theoretical schemes: a semilocal functional (PBE), a conventional hybrid functional (PBE0) in both the Kohn–Sham (KS) and generalized Kohn–Sham (GKS) scheme, a self-interaction-corrected (SIC) calculation, with and without additional “stretching” of the energy axis, GW calculations based on a PBE starting point, and optimally tuned range-separated hybrid (OT-RSH) calculations with two different short-range exchange parameters, α = 0 and α = 0.2. See text for further details on the computational approaches. The absolute HOMO energy, in electronvolts, is given at the top of each column.

First, we compare our OT-RSH results to those obtained with their “parent” semilocal functional—PBE—and those obtained from PBE0, the most popular “conventional” hybrid functional based on PBE. In the latter functional, semilocal PBE exchange is mixed with 25% nonlocal Fock exchange for the entire interaction range. This leads to an orbital-dependent functional for the exchange-correlation energy,[26] which can then be used within DFT in two different ways. If one wishes to remain within the KS framework, one must take the variational derivative of the exchange-correlation energy with respect to the density, so as to determine the multiplicative KS potential. This can be achieved for an implicit density-functional by solving the (integro-differential) optimized effective potential (OEP) equation.[26,109−111] Here, we solve this equation within the KLI approximation[91] (for details, see the “Computational Details” section) and refer to the result as PBE0KS. Alternatively, one can minimize the total energy with respect to the orbitals. This is the almost universal practice with hybrid functionals, which (as explained in the Introduction) is rigorously justified within the generalized Kohn–Sham scheme. We refer to this result as PBE0GKS.

A first observation is that the shifted PBE0KS eigenvalues are virtually identical to the shifted PBE ones for all systems. This is in agreement with similar observations reported in ref (53). One notable difference between the PBE and PBE0KS results is that the shift needed to align the PBE0KS HOMO energy to the experimental IP energy is substantially smaller than in PBE, albeit still significant. More profound differences occur between the PBE0GKS and the PBE0KS data. While the amount of shift needed for alignment of the HOMO with experimental results is essentially the same for PBE0GKS and PBE0KS, the PBE0GKS eigenvalues are overall “stretched” (= energy rescaled) with respect to the PBE0KS data—as seen most clearly for the case of benzene. For pyridine and pyrimidine, some additional orbital reordering is found, especially in the upper part of the valence band, due to mitigation of self-interaction errors, as discussed below. Such a “stretching” of the energy scale has been observed previously (see, e.g., refs (53), (112), and (113)): Körzdörfer and Kümmel[53] have rationalized it by arguing that, if the differences in the shapes of the KS and GKS orbitals are ignored, then the difference between GKS and corresponding KS eigenvalue can be viewed in terms of first-order perturbation theory involving the difference between Fock and semilocal exchange. They suggested that this first-order correction mimics successfully the first-order correction between KS values and ionization energies, which is known to dominate for outer-valence electrons.[3,24,100] Indeed, the shifted PBE0GKS eigenvalues are in much better agreement with GW than the PBE0KS ones. This underscores the beneficial effect of the nonlocal potential operator that is inherent to the GKS scheme but absent in the KS schemes.

As a next step, the role played by the self-interaction error (SIE) in the description of these systems shall be assessed, also in order to understand whether orbital SIE considerations may explain the discrepancy found above for the ring-shaped π orbital. In order to examine that, we have performed SIC calculations within the KS scheme (for details, see the “Computational Details” section), which results in SIE-free orbitals and eigenvalues. Generally speaking, the shifted SIC spectra are quite similar (but not identical) to the shifted PBE0KS spectra for all systems, although some eigenvalue differences between the two spectra (up to ∼0.3 eV for benzene and ∼0.4 eV for pyridine and pyrimidine) are found. This leads to the above-mentioned orbital reordering found for pyridine and pyrimidine once the SIE is mitigated. We furthermore find that SIC modifies PBE eigenvalues differently for π and σ orbitals, by an average of 0.2 eV, which in particular leads to a reordering of the close-lying HOMO–1 and HOMO–2 of pyrimidine. Notably, these differences in the SIC corrections for π and σ PBE orbital energies are rather similar for all three systems, meaning that the relative differences of SIEs in the σ and π orbitals of benzene are similar to the ones in pyridine and pyrimidine. We further note that, for these particular systems, the magnitude of the SIC modifications is moderate, suggesting a moderate SIE in these systems.

Comparison of the SIC data to our GW calculations reveals immediately that, just like the PBE0KS data, they would benefit from further “stretching.” This is fully consistent with the fact that both SIC and PBE0KS spectra arise from a KS, rather than a GKS, calculation. We can simulate the effect of such “stretching” on the SIC spectrum by multiplying the shifted SIC eigenvalues with a constant (a procedure that is, in fact, often performed heuristically in the comparison of KS data with experimental spectra—see, e.g., refs (114) and (115)). For each SIC spectrum, we chose this stretching factor such that the lowest-energy PBE0GKS and SIC eigenvalues agree. For all three prototypical systems, the resulting stretched SIC eigenvalue spectra are in overall better agreement with the PBE0GKS spectra than the regular (unstretched) SIC spectra and, furthermore, are in very good agreement with the OT-RSH data and with GW. For a few of the pyridine and pyrimidine orbitals, we find some remaining deviations of the stretched SIC from the PBE0GKS and OT-RSH results, which is, however, expected as the three approaches are very different. This agreement shows that a (range-separated or conventional) hybrid can be highly effective in the quantitative mitigation of SIEs, even though it is not rigorously self-interaction free because only a fraction of exact exchange is employed in the short range. Indeed, similar observations have been made for many other molecular systems,[33−35,39,53,58,61,112,116] often by comparison to experimental results, but also by comparison to SIC.[36,68] This rules out that the poor description of the “ring-shaped” (orange-colored) π orbital by OT-RSH is due to a general lack of balance in the description of π and σ orbitals owing to the SIE.

Instead, we attribute the incorrect description of this highly delocalized orbital to the absence of beyond-PBE (and possibly nonlocal) correlation. This is supported by the observation that this π orbital is also overbound in the shifted PBE0GKS data, albeit by a lesser amount of ∼0.2 to ∼0.35 eV. As semilocal exchange is known to mimic static correlation,[117,118] it stands to reason that further removal of some of it, as in the OT-RSH scheme, would worsen the situation. From this point of view, the error in this orbital is then yet one more issue of compatibility between exact exchange and correlation, especially a static one.[26] However, one cannot rule out that even if the correlation issue is overcome, an error would still remain, owing to further dynamic relaxation which is not captured in any eigenvalue-based and thus static description. Indications for such effects have been reported in the literature for small metal clusters,[119] and the response of a very delocalized orbital may be “metallic like.” Overcoming this drawback, e.g., by examining more advanced correlation functionals is, therefore, a challenge for future work.

More Complex Aromatic Heterocycles—Results and Discussion

Encouraged by the overall success of the OT-RSH approach for the prototypical aromatic rings, we now turn to examining the performance of the OT-RSH method for more complex organic molecules: terpyrimidinethiol (3N-thiol) and copper phthalocyanine (CuPc)—see Figure 1b.

Pyrimidinethiols such as 3N-thiol are known to form self-assembled monolayers (SAMs),[120] and they were shown to display a number of interesting phenomena; e.g., they exhibit diode-type current–voltage characteristics in molecular-scale electronic devices.[66,121] In the context of 3N-thiol SAMs, some of us have predicted theoretically that the electronic structure is significantly altered due to collective electrostatic effects in the SAM, leading to a localization of the frontier molecular orbitals and a concomitant pronounced reduction in the HOMO–LUMO gap.[68] 3N-thiols were also suggested as SAMs to strongly reduce or enhance the work function of an underlying metal.[67] These appealing findings render 3N-thiols an interesting candidate for applications in organic and molecular electronics. 3N-thiols are also interesting from a methodological perspective, as they pose a serious challenge to semilocal KS schemes:[68] By means of SIC calculations analogous to the ones performed here, it was shown that close-lying π and σ orbitals in the outer-valence region of 3N-thiol have markedly different SIEs and, as a consequence, are wrongly ordered in (semi)local DFT functionals such as PBE. It was additionally found that the Heyd–Scuseria–Ernezerhof (HSE) SR hybrid functional[122,123] to a large extent reproduces the SIC corrections to the PBE outer-valence eigenvalues. In HSE, the SR components are the same as in PBE0, but in the LR there is no nonlocal exchange. From the photoemission spectroscopy point of view, the behavior of HSE is known to be similar to that of conventional hybrid functionals.[34,124−128]

The above-mentioned findings are supported by the data in Figure 6a, which shows the shifted eigenvalue spectra of 3N-thiol, i.e., aligning all HOMO levels to zero, starting with the literature data[68] for PBE, SIC, and HSE. Here, we focus on the range of ∼3 eV below the HOMO, which we have previously identified for molecules of a similar size as a useful range of accuracy for the RSH-derived eigenvalues.[58] Generally, PBE tends to produce σ orbitals that are underbound owing to significant SIE.[33,34,58,68] Because the HOMO obtained from PBE is of σ nature, shifting according to its energy results in π orbitals that appear to be overbound. The SIC calculations (as above, the spectrum is shown both “as calculated” and “stretched”) correct the SIE and, as discussed above, strongly shift the σ orbitals down in energy with respect to the π ones, leading to extensive reordering of the eigenvalue spectrum. This effect also applies to the HSE spectrum, which needs no stretching as it originates from a GKS calculation. We additionally performed PBE0GKS calculations, also shown in Figure 6a,[129] and found an overall very good agreement with the HSE spectrum and with the “stretched” SIC with an average deviation of the shifted eigenvalues around ∼0.1 eV. This is consistent with the above-mentioned findings, namely that HSE is comparable to conventional hybrid functionals as far as photoemission spectra are concerned.

Figure 6

(a) Shifted eigenvalue spectra of 3N-thiol as obtained from different theoretical schemes (see text for details). PBE, SIC, and HSE values were taken from ref (68). For the OT-RSH calculations, the optimal γ value was 0.217 and 0.187 (in Bohr−1) for α = 0 and α = 0.15, respectively. The absolute HOMO energy, in electronvolts, is given at the top of each column. (b) Unshifted eigenvalue spectra of the same system, as calculated from HSE, PBE0GKS, OT-RSH calculations with the two different α values, and GW.

To the best of our knowledge, an experimental PES spectrum of 3N-thiol is not available in the literature. Therefore, we performed reference GW calculations for which the (again shifted) eigenvalue spectrum is also given in Figure 6a. As seen above for the prototypical aromatic rings, the agreement of the shifted GW spectrum with the shifted hybrid functional (PBE0 or HSE) spectrum is very good. Some discrepancy is observed between ∼−2.0 to ∼−2.5 eV, but in this context it should be kept in mind that here the plotted energy range is much smaller than for the aromatic rings in Figures 4 and 5. In fact, in the present case, the errors do not exceed 0.2 eV, which is within the expected level of agreement between the two methods (vide supra). This confirms the conclusion drawn already from the data of the “prototypical” molecules that the partial correction of SIE and a “stretching” of the spectrum are the essential effects of a hybrid functional.

Despite the excellent agreement of the shifted conventional hybrid (and stretched SIC) data with the GW results, the two differ by a large rigid energy shift. This can be gleaned even from Figure 6a, where the HOMO eigenvalue obtained from each method is given at the top of each column, but is made obvious by Figure 6b, which shows unshifted spectra on an absolute energy scale. This discrepancy is removed by the OT-RSH results, also shown in Figure 6a and b (again with two different choices of SR exact exchange—α = 0 and α = 0.15—the motivation for this particular selection of α values is discussed below). As with the simpler systems, the IPs deduced from our computed OT-RSH and GW eigenvalues are very close in energy (∼0.2 eV) for both choices of α. Consequently, the OT-RSH scheme produces quantitatively meaningful spectra also on an absolute energy scale. This is in sharp contrast to all other DFT methods presented in Figure 6.

A remaining and important issue is the best choice for the short-range Fock-exchange fraction, α. In contrast to the case of the small aromatic rings, where different choices for α have had no qualitative effect and only a small quantitative effect, for the 3N-thiol system the differences between the OT-RSH spectra for the two choices of α are more pronounced. As shown in Figure 6, these differences include qualitative deviations, e.g., a reordering of the HOMO–3 and HOMO–4 orbitals. In particular, we observe that, while the position of the π orbitals is similar in both spectra, the σ orbitals are systematically underbound in the α = 0 calculation, in comparison to the α = 0.15 or the GW calculations. The different dependence of π and σ orbitals on α is underscored in Figure 7, which shows the dependence of orbital eigenvalues on α in the entire range of α = 0 to 0.2, with the range-separation parameter, γ, optimized individually for each choice of α. We relate the underbinding of σ orbitals at small α to the fact that the spatial extent of the σ and π orbitals within the molecular plane is significantly different; i.e., the σ orbitals are clearly more localized than π states and, therefore, should suffer from a larger SIE.[68] The finding that the agreement with our GW calculations is significantly improved for α = 0.15 underscores that, as suggested in ref (58) and discussed in the Introduction, a fraction of Fock exchange in the SR together with full Fock exchange in the LR allows for a mitigation of SIE and therefore a balanced treatment of differently localized molecular orbitals, while retaining the correct absolute position of the energy levels. In principle, a similar statement applies for the small aromatic rings as well. However, there the difference in degree of localization is not the same and, perhaps more importantly, the energy differences between individual eigenvalues are much larger. Therefore, the practical importance of α is smaller in the prototypical small (hetero)cycles, as shown in Figure 4 and in the SI.

Figure 7

Eigenvalues of 3N-thiol, obtained from optimally tuned range-separated hybrid (OT-RSH) calculations, as a function of the short-range exchange fraction, α, with the optimal value of the range-separation parameter, γopt (in Bohr–1) determined for each choice of α.

These considerations raise the question of whether an optimal α can be chosen a priori. Srebro and Autschbach[60] have suggested that it can be achieved by tuning α. Specifically, they sought the value of α that minimizes the curvature of the ideally piecewise-linear total energy versus particle number curve for the addition and removal of one charge from the neutral molecule. This was done while simultaneously determining the optimal γ for each choice of α obtained from eq 4 (with i restricted to 0 and 1). As discussed in ref (58), some of the present authors successfully employed a similar approach to determine an optimal value of α also for the purpose of obtaining an outer-valence eigenvalue spectra, but only if it also involved the removal of a second electron. Specifically, αopt was found to equal 0.2 for PTCDA and NTCDA. This value agrees well with the above considerations. Recently, Stein et al. established a rigorous quantitative equality between deviations from piecewise linearity and deviations from the IP theorem.[130] Therefore, one can equivalently seek α directly by minimizing the target function J from eq 4 without an explicit consideration of fractional densities. Indeed, if this procedure is applied to PTCDA and NTCDA, the same optimal α value is found (see SI). Therefore, we discuss these two tuning procedures (by piecewise linearity and by the IP theorem) together in the following.

Indeed, in cases where an optimal α can be clearly identified, several properties, including the eigenvalue spectrum, have been found to be predicted satisfactorily.[58,60,131] However, in the case of 3N-thiol, we could not employ that strategy: With i restricted to 0 and 1, a similar minimal value of the target function J of eq 4 is obtained across a large range of α values: the residual J only changed by 0.01 eV from α = 0.00 to α = 0.60, and no minimum of J occurred. Removal of the second electron resulted in orbital reordering (see SI for more details).

In the absence of viable alternatives, a possible strategy would be to resort to the generally recommended value[74] of α = 0.2, which indeed has been shown to yield very good results, for the prototypes studied above as well as for other systems.[59] This would also be the case here; i.e., the MAD between the absolute eigenenergies from our GW calculations and OT-RSH with α = 0.2, for the states shown in Figure 7, is ∼0.1 eV. However, as the 3N-thiol results do depend more sensitively on α, it is still interesting, if possible, to determine an optimal α value from additional considerations. In this context, a useful practical observation is that the shifted results of the conventional hybrid functional PBE0 are in good agreement with experimental and/or theoretical reference data for the systems discussed here, as well as for many other cases.[35,58,59,132] Therefore, a pragmatic approach would be to simply tune α so as to obtain agreement between splitting of eigenvalues in the OT-RSH and in the PBE0 spectra. Specifically, one should seek a balance in the relative description of delocalized π and localized σ states. The easiest way to achieve that is to tune α such that the π–σ energy difference between the HOMO and HOMO–1 states is the same in the OT-RSH and the PBE0GKS calculations. This allows us to obtain the same useful level of SIE mitigation as in a conventional hybrid functional, without the need for spectral shifting. It is this approach which has led to the value of α = 0.15 for which the data shown in Figure 6 have been obtained.

Comparing our unshifted GW and OT-RSH (α = 0.15) calculations of Figure 6b, one, indeed, observes a very good agreement between the two methods. The only clear difference is some orbital reordering around ∼−10.8 eV. One should, however, note that the eigenvalues clustered there are very close in energy and that this reordering does not involve deviations greater than ∼0.1 eV between corresponding eigenvalues of the two methods. Overall, the MAD between the GW and OT-RSH (α = 0.15) calculations is a satisfying ∼0.1 eV, with the largest deviation being ∼0.2 eV for the HOMO. For comparison, with OT-RSH (α = 0) the MAD is ∼0.25 eV, with the largest deviation being ∼0.45 eV, where the reduced level of agreement is mostly due to the less accurate description of the σ orbitals.

Encouraged by this further success, we now turn to an even more complex system—copper phthalocyanine (CuPc—see Figure 1b). In molecular-solid form, CuPc is a highly stable organic semiconductor with a broad range of applications, including light-emitting diodes, solar cells, gas sensors, and thin-film transistors.[133] Owing to these applications, there is considerable interest in investigating its electronic structure (see ref (69) and additional references therein). In the present context, CuPc mainly serves as a test case that poses several additional challenges for the OT-RSH method: First, it is an open shell molecule (s = 1/2). Second, the interaction between the d orbitals of the copper atom and the s and p orbitals of the embedding macrocycle result in a highly nontrivial set of localized and delocalized orbitals—an issue elaborated below.

For the CuPc molecule, Marom et al. have previously established that: (1) semilocal functionals, such as PBE, result in eigenvalue spectra that are strongly distorted by severe SIE;[34] (2) these errors are mitigated substantially, though not completely, by the use of hybrid functionals such as PBE0 or HSE;[34,69] (3) severe SIE distortions at the DFT level are partly carried over to GW calculations building on the DFT densities.[69]

To better understand the above claims, consider Figure 8a, which provides a scheme containing the computed frontier eigenvalues and orbitals. We illustrate the severe SIE (claim 1) by considering several selected frontier orbitals. Whereas the a1 orbital and the doubly degenerate e orbitals are highly delocalized on the macrocycle, the spin-split b1 orbitals are quite strongly localized around the copper atom. Configuration 1 (left side of Figure 8a) was obtained from the OT-RSH calculations discussed below but is equivalent to the one obtained from the GW calculations of ref (69). It identifies the delocalized a1 and e orbitals as the HOMO and LUMO, respectively, and places the spin-split b1 and b1 orbitals as HOMO–1 and LUMO+1, respectively. PBE calculations, however, predict configuration 2 (right side of Figure 8a), in which the spin-split b1 and b1 are identified as HOMO and LUMO, respectively. Clearly, the localized b1 orbital is spuriously pushed to higher energies by the SIE, the effect being strong enough to make the b1 orbital become the HOMO. This forces the unoccupied b1 orbital to be spuriously shifted to lower energies in order to maintain the spin-splitting symmetry, to the point of it becoming the LUMO. PBE0 or HSE strongly mitigate this error (claim 2) and yield configuration 1.

Figure 8

(a) Schematic diagram of selected frontier eigenvalues and orbitals for the CuPc molecule, as obtained from an OT-RSH calculation (configuration 1, left) and from a PBE calculation (configuration 2, right). (b) OT-RSH charge-density differences between neutral and cation (left) and LUMO of cation (right), obtained from the “open-shell singlet” configuration of the CuPc cation, corresponding to ionization of netrual configuration 1 (in green), and from the “closed-shell singlet” configuration of the CuPc cation, corresponding to ionization of configuration 2 (in red). In the charge density difference plots, green or red (blue) regions denote areas of electron density depletion (accumulation) as a consequence of the ionization process.

To understand the manifestation of the SIE (and its mitigation) in the simulated photoelectron data, spectra computed with different computational methods are compared to experimental photoemission data in Figure 9. Importantly, in this figure we compare occupied-state eigenvalues to the gas-phase photoelectron spectrum as before, but further compare unoccupied-state eigenvalues to experimental inverse photoemission spectroscopy (IPES). In IPES, photons are emitted from a sample due to its irradiation with fixed-energy electrons, and the energy distribution of the emitted photons is measured, yielding experimental information on virtual states.[1] CuPc is the first system in this article for which such comparison is possible, because for the small aromatic rings the virtual states are unbound and because for 3N-thiol no (regular or inverse) experimental PES data are available. In fact, gas-phase inverse photoemission spectroscopy is nonexistent in general. Therefore, the comparison is to experimental data obtained from thin CuPc films. Due to polarization effects,[59,134,135] the electron affinity of a film is much smaller than that of an isolated molecule, and the computed empty state energies and the measured inverse photoemission spectroscopy data can be compared only up to a rigid shift. Therefore, to allow comparison to experimental results without modifying the computational gas-phase data, the experimental IPES spectra were shifted so as to align the lowest-energy peak with the GW spectrum shown at the topmost computed spectrum in Figure 9.

Figure 9

Simulated DFT and GW spectra (see text for details), obtained from computed energy levels (shown as sticks) by broadening via convolution with a 0.15-eV-wide Gaussian. PBE, PBE0, and GW data were taken from ref (69). Theoretical data are compared to the experimental gas phase photoemission data of Evangelista et al.[141] and to the thin film inverse photoemission data of (a) Murdey et al.[142] (shown with curve fitting results) and (b) Hill et al.[143] The experimental inverse photoemission spectra were shifted so as to align the LUMO peak with the computed GW@PBE0 LUMO peak. On the OT-RSH data, optimal γ values are given in Bohr–1. Eigenvalues corresponding to the e, a1, and b1 orbitals are designated by the same color scheme as in Figure 8.

As discussed above, it is well-known that for either PBE or PBE0 the HOMO and LUMO do not correspond to the ionization potential or the electron affinity, respectively (see, e.g., refs (26) and (32)). This causes an uncontrolled rigid shift of the simulated photoemission curve. To facilitate a meaningful comparison also for these functionals, the PBE and PBE0 spectra of Figure 9 are additionally shown in shifted form. To determine the required rigid shift, the ionization potential is computed as the total energy difference between the neutral species and the cation, and the filled-state eigenvalue spectrum is rigidly shifted such that the HOMO energy coincides with the computed ionization potential. An equivalent procedure is applied to the unoccupied-state eigenvalue spectrum, which is rigidly shifted such that the LUMO energy coincides with the computed electron affinity. Comparison of the PBE- and PBE0-based simulated spectra in Figure 9 (taken from ref (69)) with the experimental data reveals that whereas the shifted PBE spectrum is in poor agreement with experiment for both filled and empty states, the shifted PBE0 spectrum provides for a much improved agreement with experimental results. This highlights both the severe SIE in the PBE calculation and its mitigation by the PBE0 calculation. Better agreement yet, including an accurate placement of the spin-split filled b1 and empty b1 orbitals, is afforded by GW calculations of Marom et al.,[69] based on a PBE0 starting point. It is in significantly better agreement with experimental results than their GW calculations based on a PBE starting point (at least with the particular flavor of GW used in ref (69), which, as discussed earlier, is different from the one we used above for benzene, pyridine, and pyrimidine; we refrained from conducting such GW calculations also for CuPc, because of the computational cost involved). This improved description of the electronic structure of CuPc by the PBE0 based GW calculations substantiates the third above claim, i.e., the possible carry-over of severe SIE to the GW calculation.

It is intriguing to determine how well the OT-RSH approach performs in this more complicated scenario. As for the case of pyridine (see Figure 2), it is important to choose the correct cationic configuration for the tuning procedure. Two cationic configurations can be obtained from the OT-RSH calculation. One is an “open shell singlet,” with a single spin in the b1 orbital and in the a1 orbital, which corresponds to removal of an electron from the a1 HOMO of configuration 1 in Figure 8a. The other is a “closed-shell singlet” cation, which corresponds to removal of an electron from the b1 HOMO of configuration 2 in Figure 8a. Importantly, only for the “open-shell singlet” configuration (which is lower in total energy) is the LUMO consistent with the character of the neutral HOMO, and with the “hole density,” as shown in Figure 8b. Just like in the pyridine example of Figure 2, all tuning must be performed solely with this configuration. We note that multiple stable configurations are well-known to occur in metal-phthalocyanines,[132] which is why we strongly stress the need for correct identification and usage of compatible neutral and ionic configurations in the tuning procedure.[136]

Having ascertained the validity of our tuning procedure, we turn to the spectra simulated using the OT-RSH method, also shown in Figure 9. As for the case of 3N-thiol, the obtained spectra (second and fourth theoretical curves in the figure) are strongly α-dependent. This is shown explicitly in Figure 10, which shows the dependence of orbital eigenvalues on α, with the range-separation parameter, γ, again optimized individually for each choice of α. It is particulary striking that the orbitals most sensitive to α are the highly localized spin-split b1 pair, which is fully consistent with our above discussion on the relation between orbital localization and sensitivity to short-range Fock exchange. Unfortunately, here direct tuning of α based on eq 4, even with i = 0, 1, −1, was not useful because the minimal value of the target function J was very small and very weakly dependent on α within the numerical accuracy of our work (see SI). Thus, while the procedure of ref (58) did not fail, in the sense of providing an incorrect value for α, it did not succeed either. Therefore, we determined α as for the case of 3N-thiol, i.e., based on the success of the shifted PBE0 spectrum. As discussed above, a crucial aspect in the theoretical description of the CuPc spectrum is the energy separation between the delocalized a1 orbital and the filled localized b1 orbital of Figure 8. Therefore, we tune the value of α so as to agree with the PBE0-calculated separation between the eigenvalues for the a1 and b1 orbitals. This yields α = 0.17, which is again close to the “default” value of 0.2, further supporting the latter as a useful choice even without any α-tuning. α = 0.17 was then used to calculate the second theoretical curve in Figure 9. It yields excellent agreement with prior GW@PBE0 calculations (shown in the figure) and the experimental data in both the top valence (i.e., occupied) states and the bottom conduction (i.e., unoccupied) states. In particular, the OT-RSH alignment of the a1, b1, and e states of Figure 8 agrees very well with the prior GW calculations. As orbitals above and below these frontier ones tend to be “clustered” for these larger molecule, we do not discuss them individually. However, it is clear that the agreement between GW and OT-RSH gradually deteriorates as one goes further down or up in the valence or conduction states, respectively. In particular, it appears that agreement of the lower valence states with the GW results could be further assisted by some “stretching” of the energy scale—an observation consistent with the generally expected differences between nonlocal exchange and self-energy operators,[3] as well as with our general guideline of restricting our attention to states within a few electronvolts of the frontier orbitals. Nevertheless, in the absence of detailed orbital information, at the level of the given experimental resolution, the OT-RSH results can be viewed as agreeing with experimental results as well as the GW ones, possibly even relatively deep into the valence band.

Figure 10

Eigenvalues of CuPc, obtained from optimally tuned range-separated hybrid (OT-RSH) calculations, as a function of the short-range exchange fraction, α, with the optimal value of the range-separation parameter, γopt (in Bohr–1), determined for different choices of α.

As an additional observation, we note that the HOMO and LUMO eigenvalues obtained from the OT-RSH (α = 0.17) calculation (as is, without any rigid shifting of the data) are in excellent agreement with the lowest quasi-hole and quasi-electron excitations, respectively, obtained from GW@PBE0. This is fully consistent with the work of Stein et al.,[52] who suggested that OT-RSH gaps can be identified with quasi-particle gaps—an observation that has since been confirmed for a large variety of systems by ourselves, as well as by several other groups.[58,59,78,137−140] A related interesting observation is that the OT-RSH (α = 0) data are rather similar to the GW@PBE calculations from ref (69) (again, for that particular flavor of GW): Both are far better than the PBE data, but in both spectra, the a1–b1 and e–b1 separations are not large enough. This highlights yet again the positive role of SR exchange in mitigating SIEs where they are significant.[58]

Conclusions

In conclusion, we examined the performance of optimally tuned range-separated hybrid functionals for predicting the photoemission spectra of several challenging prototypical (benzene, pyridine, and pyrimidine) and complex (3N-thiol and CuPc) (hetero)cyclic organic molecules. Overall, the resulting spectra agree very well (typically within ∼0.1–0.2 eV) with experimental data and our GW calculations, even for states several electronvolts away from the frontier orbitals. We performed several additional hybrid DFT calculations in both the KS and GKS schemes and found that the inclusion of a nonlocal operator strongly benefits the calculated spectrum. Moreover, our SIC calculations confirmed that self-interaction errors can be efficiently mitigated in OT-RSH functionals, which shows that, with a PBE0-based optimal choice of the short-range fraction of Fock exchange, the OT-RSH method can offer an excellent balance in the description of localized and delocalized states. The sole exception found in the studied systems is a high-symmetry orbital, particular to small aromatic rings. We conclude that the OT-RSH method is a highly accurate DFT method for outer-valence PES prediction for such systems, and its accuracy is comparable to state-of-the-art GW schemes. This success comes at the price of increased computational and conceptual cost that is inherent to the parameter tuning. While this increase in effort is notable compared to standard, nontuned DFT calculations, it is not overwhelming and still lower than the cost of a GW calculation.