Accurate Valence Ionization Energies from KohnâSham Eigenvalues

Article pubs.acs.org/JCTC

Accurate Valence Ionization Energies from Kohn−Sham Eigenvalues with the Help of Potential Adjustors Adrian Thierbach,† Christian Neiss,† Lukas Gallandi,‡ Noa Marom,¶ Thomas Körzdörfer,‡ and Andreas Görling*,† †

Lehrstuhl für Theoretische Chemie, Universität Erlangen-Nürnberg, Egerlandstrasse 3, D-91058 Erlangen, Germany Computational Chemistry, University of Potsdam, D-14476 Potsdam, Germany ¶ Materials Science and Engineering, Carnegie Mellon University, 5000 Forbes Avenue, Pittsburgh, Pennsylvania 15213-3890, United States ‡

S Supporting Information *

ABSTRACT: An accurate yet computationally very efficient and formally well justified approach to calculate molecular ionization potentials is presented and tested. The first as well as higher ionization potentials are obtained as the negatives of the Kohn−Sham eigenvalues of the neutral molecule after adjusting the eigenvalues by a recently [Görling Phys. Rev. B 2015, 91, 245120] introduced potential adjustor for exchangecorrelation potentials. Technically the method is very simple. Besides a Kohn−Sham calculation of the neutral molecule, only a second Kohn−Sham calculation of the cation is required. The eigenvalue spectrum of the neutral molecule is shifted such that the negative of the eigenvalue of the highest occupied molecular orbital equals the energy difference of the total electronic energies of the cation minus the neutral molecule. For the first ionization potential this simply amounts to a ΔSCF calculation. Then, the higher ionization potentials are obtained as the negatives of the correspondingly shifted Kohn−Sham eigenvalues. Importantly, this shift of the Kohn−Sham eigenvalue spectrum is not just ad hoc. In fact, it is formally necessary for the physically correct energetic adjustment of the eigenvalue spectrum as it results from ensemble density-functional theory. An analogous approach for electron affinities is equally well obtained and justified. To illustrate the practical benefits of the approach, we calculate the valence ionization energies of test sets of small- and medium-sized molecules and photoelectron spectra of mediumsized electron acceptor molecules using a typical semilocal (PBE) and two typical global hybrid functionals (B3LYP and PBE0). The potential adjusted B3LYP and PBE0 eigenvalues yield valence ionization potentials that are in very good agreement with experimental values, reaching an accuracy that is as good as the best G0W0 methods, however, at much lower computational costs. The potential adjusted PBE eigenvalues result in somewhat less accurate ionization energies, which, however, are almost as accurate as those obtained from the most commonly used G0W0 variants. (ADC) scheme20 are currently considered to be state of the art in predicting accurate quasi-particle energies for molecules with some tens to some hundreds of electrons. In a recent benchmark study for organic molecular acceptors,21−24 basisset extrapolated third-order ADC [ADC(3)] has been demonstrated to predict vertical IPs and EAs with high accuracy, yielding mean absolute errors (MAEs) of 0.12 and 0.16 eV, respectively.24 Similar MAEs were found from the best nonself-consistent GW calculations,23 which hereafter will be referred to as G0W0. However, by involving a one-shot correction to eigenvalues obtained from density functional theory (DFT) or Hartree−Fock (HF) calculations, the accuracy of G0W0 calculations depends heavily on the underlying method.25−31 The starting-point dependence of G0W0 can not only lead to significant errors in absolute IPs and

1. INTRODUCTION The efficient and yet accurate prediction of ionization potentials (IPs) and electron affinities (EAs) of atoms, molecules, clusters, and solids has been a long-standing challenge for electronic structure theory. Theoretical calculations are frequently used in combination with photoelectron spectroscopy (PES) and inverse photoelectron spectroscopy measurements, yielding far-reaching physical insight by connecting PES peaks with molecular orbitals.1−10 Also recent computational screening approaches to find new materials for photovoltaics and molecular electronics have profited significantly from efficient and reliable ways to predict accurate IPs, EAs, and fundamental gaps, since these are key factors determining the energetics of the creation, transport, and annihilation of charge carriers.11−16 Ab initio methods based on many-body perturbation theory such as the GW approximation17,18 or polarization propagator methods19 such as the Algebraic Diagrammatic Construction © 2017 American Chemical Society

Received: May 16, 2017 Published: August 7, 2017 4726

DOI: 10.1021/acs.jctc.7b00490 J. Chem. Theory Comput. 2017, 13, 4726−4740

Article

Journal of Chemical Theory and Computation

energy differences obtained with nonconventional KS methods that treat exchange exactly and correlation via the random phase approximation.53,56 In this work, however, we will focus on finite systems, in particular molecules. For molecules it has a long tradition to identify the negatives of the KS or GKS eigenvalues with IPs or EAs, as exemplified by refs 1−7, 9, 10, and 57−61. However, it is immediately clear that such an identification is fundamentally flawed for formal reasons. For a neutral molecule with a fixed integer electron number N, the KS or GKS eigenvalue spectrum is only defined up to an additive constant, because the exchange-correlation potential is defined only up to an additive constant; see below or ref 52 for details. In an ensemble KS or GKS formalism, two exchange-correlation potentials and, thus, two sets of eigenvalues differing by a constant are associated with the neutral N-electron molecule.62−65 One exchange-correlation potential, vN− xc , is obtained if the limit of N electrons is approached from the electron deficient side and leads to an eigenvalue spectrum whose HOMO eigenvalue equals the negative of the IP in case of the exact formalism. The other exchange-correlation potential, vN+ xc , results if the limit of N electrons is approached from the electron abundant side and leads to an eigenvalue spectrum whose LUMO eigenvalue equals the negative of the EA in case of the exact formalism. N+ N− The difference between vxc and vxc is the derivative discontinuity of the exchange-correlation potential.62−65 Unless the derivative discontinuity was zero, it is impossible that a single KS or GKS eigenvalue spectrum can be identified both with IPs and EAs at the same time. In ref 52 it was suggested to construct an approximate quasiparticle spectrum, i.e., a spectrum containing the negatives of the IPs and EAs, from the two different KS or GKS eigenvalue spectra associated with the neutral molecule. This potential adjustor approach combines the eigenvalues of the occupied orbitals obtained with the exchange-correlation potential vN− xc with the eigenvalues of the unoccupied orbitals obtained with the exchange-correlation potential vN+ xc . In an exact KS or GKS formalism, the resulting spectrum would yield the exact first IP and EA as the negative of the HOMO and LUMO eigenvalues, respectively. Then, the other KS or GKS eigenvalues are well-defined approximations to IPs and EAs.37,66 N+ The question that remains is how to obtain vN− xc and vxc and the associated orbital eigenvalues in case of approximate KS or GKS approaches in practice. Taking straightforwardly the functional derivative of the approximate exchange-correlation energy with respect to the electron density yields only one exchange-correlation potential ṽxc, which neither equals vN− xc nor vN+ xc . In ref 52 it was shown that, with the help of potential N+ adjustors, the potential ṽxc can be shifted to obtain vN− xc or vxc . The potential adjustors have a sound formal justification. In the N+ exact KS formalism, they can be used to obtain vN− xc or vxc from any exchange-correlation potential vx̅ c that differs from vN− xc or vN+ xc by an arbitrary constant. Indeed the potential adjustors N+ enable an exact construction of vN− xc of vxc exclusively from quantities defined in the KS formalism for fixed particle numbers, that is, without the need to resort to quantities accessible only in the ensemble KS formalism. The potential adjustors are straightforwardly obtained from the ΔSCF values of the first IP and the first EA.52 The required N+ KS or GKS eigenvalue spectra associated with vN− xc and vxc are accessible by shifting the KS or GKS eigenvalue spectrum of the neutral molecule such that its HOMO eigenvalue equals the

EAs, but it has also been demonstrated to cause an incorrect ordering of the valence orbitals for some systems.26,27,31−33 This error can be traced back to orbital self-interaction errors34−36 in the underlying DFT calculations and, therefore, can be addressed by using DFT starting points with improved orbital ordering such as global hybrid functionals with a suitable amount of HF-exchange,27,37 tuned long-range corrected hybrid functionals,38−40 or self-interaction-corrections.6,34,35 The unfavorable scaling of GW methods [N4] and ADC(3) [N6] with system size [N] and, in particular, the high memory needs of such calculations, however, limit their applicability to small- and medium-sized systems. An alternative approach that is computationally less expensive than ADC(3) and GW methods is provided by the Kohn−Sham (KS)41 or generalized KS (GKS)42 formalisms of DFT.43−45 Here, the first IP as well as the first EA can be calculated as the energy difference between two electronic ground states energies, i.e., of the cationic and the neutral system or the neutral and the anionic system, respectively. For molecules or clusters, calculations of such energy differences are straightforward with KS and GKS methods, formally well justified, and run under the name ΔSCF (Δ-Self-ConsistentField) calculations. Even though ΔSCF calculations are not size-consistent for conventional KS or GKS methods, such calculations are often more accurate than many GW variants for first IPs,22,23 as will also be demonstrated below. With the term conventional KS and GKS methods we here denote methods based on the local density approximation (LDA) and generalized gradient approximation (GGA) for the exchangecorrelation potential or methods relying on meta GGA or hybrid functionals.43−45 KS methods employing the optimized potential method46−51 for the construction of the exchangecorrelation potential, on the other hand, do not fall under this category. Contrary to a common belief, ΔSCF calculations can also easily be carried out for infinite periodic systems, e.g., crystalline solids, on firm formal ground.52−54 Indeed, in the case of infinite periodic systems, ΔSCF calculations reduce to additional energy evaluations with one electron added or removed using the orbitals of the neutral system. This means no SCF calculations besides the one for the neutral system are necessary; only additional energy evaluations for the electronic system with one electron added or removed are required. The reason for this is that, in an infinite periodic system, the addition or removal of a single electron does not lead to a relaxation of the orbitals. However, for LDA, GGA, and hybrid methods the total energy differences between the neutral and the cationic system and between the neutral and the anionic system equal the negatives of the eigenvalues of the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbitals (LUMO), respectively.52,54,55 Typically, in periodic systems it is the fundamental quasiparticle band gap, and not the individual IP and EA, that is of primary interest. If IP and EA are calculated by energy differences, and the latter equals the negative of the HOMO and LUMO eigenvalues, as in LDA, GGA, or hybrid methods, then the quasiparticle band gap is given by the KS or GKS HOMO− LUMO eigenvalue gap. For LDA or GGA methods, such HOMO−LUMO eigenvalue gaps typically are in poor agreement with experimental values, whereas some hybrid methods yield band gaps of much better accuracy. Furthermore, it was recently shown that fundamental quasiparticle gaps of semiconductors can be calculated with high accuracy from total 4727


Article

Journal of Chemical Theory and Computation negative of the ΔSCF value for the first IP or, respectively, such that the LUMO eigenvalue equals the negative of the ΔSCF value of the first EA. The approximate quasiparticle spectra resulting from the described shifts of the original KS or GKS eigenvalues shall be denoted as potential adjusted KS (paKS) or GKS (paGKS) spectra, respectively. Importantly, the shift of the original KS or GKS eigenvalues corresponds to adding the potential adjustors defined in ref 52 to the usual KS or GKS exchange-correlation potential ṽxc obtained when straightforwardly taking the functional derivative of the exchangecorrelation energy with respect to the electron density, see sections 2 and 3 for details. Effectively, this means that ΔSCF calculations of the first IP and the first EA additionally give access also to approximate higher IPs and EAs via the shifted eigenvalue spectra. The crucial question is how good this disarmingly simple, yet formally well-justified approach, which was propagated in ref 52 on formal grounds but not applied to any electronic system, would perform in practice. To answer this question we applied the approach to various molecules, that is we used the approach to calculate IPs of valence as well as core electrons for a set of molecules and compared the results with G0W0 quasiparticle spectra as well as experimental gas-phase photoelectron spectra. The paper is organized as follows: In the next two sections the formalism behind the paKS and paGKS methods is discussed. Then follows a section which compares the methods to G0W0 calculations and experimental data. Finally conclusions are drawn.

δTs = −

δTs δρ(r)

=−

∫ d rvs̅(r)δρ(r)

vH̅ (r) =

(4)

δU δρ(r)

(6)

with the Coulomb energy U, and the exchange-correlation potential vxc̅ (r) =

δExc δρ(r)

(7)

with the exchange-correlation energy Exc. The Hartree and the exchange-correlation potential are, again, only defined up to an additive constant in the KS formalism for a fixed integer electron number. That is vH̅ and vx̅ c represent sets of potentials obtained by adding arbitrary constants to one member of the set, respectively. In order to determine the absolute energetic positions of the potentials vs̅ , vH̅ , and vx̅ c, we have to turn to the ensemble KS formalism. A straightforward definition for the Coulomb energy in the ensemble KS formalism is U[ρ] = (1/2)∫ drdr′ ρ(r)ρ(r′)/|r − r′|, which can be evaluated for any electron density, including those that integrate to noninteger electron numbers. The Hartree potential is then given by the classical electrostatic potential vH(r) = ∫ dr′ ρ(r′)/|r − r′| of the electron density. In the standard ensemble KS formalism, the Hartree potential vH(r) is defined without the freedom of adding an arbitrary constant, that is, it is constructed such that it falls to zero for large |r|. Determining the exchange-correlation potential in the ensemble KS formalism is less straightforward, since it is known to have a derivative discontinuity at integer electron numbers.62−65 Because typically exchange-correlation potentials are required for an integer number of electrons, two limits have to be considered for the exchange-correlation potential, (i) coming from the electron deficient side corresponding to the limit

(1)

(2)

which is given by the functional derivative of the ’noninteracting’ kinetic energy Ts with respect to the electron density ρ(r). This functional derivative is only defined up to an arbitrary additive constant ν in the KS formalism for a fixed integer electron number. In other words, as indicated by the bar, vs̅ stands for a set of potentials that can be obtained by adding an arbitrary constant to one member of the set. Importantly, the fact that the KS potential is defined only up to an additive constant is a consequence of the fact that any change δρ of the electron density ρ(r) has to integrate to zero, i.e.,

∫ drδρ = 0

∫ d r[vs̃(r) + ν]δρ(r)

vs̅ (r) = vext(r) + vH̅ (r) + vxc̅ (r) (5) of the external potential vext(r), usually the potential of the nuclei, the Hartree potential

with the effective KS potential vs̅ (r) = −

=−

with ν, again, denoting an arbitrary constant. The three righthand sides of eq 4 equal each other due to eq 3, which, of course, can be multiplied by any constant ν. Eq 4 demonstrates that each member of the set of potentials vs̅ equally well acts as a functional derivative of the kinetic energy Ts with respect to the electron density or, vice versa, that the functional derivative of Ts is to be identified with the whole set of potentials vs̅ . It is important to realize that this is a direct consequence of the fact that the electron number is fixed. Hence, in this case, the undefined additive constant in the KS potential has nothing to do with a Lagrange multiplier or a chemical potential. The KS potential is given as the sum

2. ENERGETIC ADJUSTMENT OF POTENTIALS For a start, we consider the Kohn−Sham (KS) formalism for a fixed integer electron number N. For simplicity, we consider neutral nonspin-polarized molecules with real-valued KS orbitals. In this case, the electron number N is even. The KS equation for the KS orbitals φi and their eigenvalues εi[vs̅ ] is given by ⎡ ∇2 ⎤ + vs̅ (r)⎥φi(r) = εi[vs̅ ]φi(r) ⎢− ⎣ 2 ⎦

∫ d rvs̃(r)δρ(r)

(3)

vxcN −(r) = lim

in order to keep the electron number fixed. If we pick one arbitrary member ṽs out of the set vs̅ , then the change δTs of the noninteracting kinetic energy Ts due to a change δρ in the electron density, is given by

q→N−

δExc δρ(r)

(8)

and (ii) coming from the electron abundant side corresponding to the limit 4728


Article

Journal of Chemical Theory and Computation vxcN +(r) = lim

q→N+

ΔNxc+[vxc̅ ] = E0N + 1 − E0N − εLUMO[vext + vH + vxc̅ ]

δExc δρ(r)

(9)

N+1 Here, EN0 , EN−1 are the electronic ground state 0 , and E0 energies of the N-, the (N − 1)-, and the (N + 1)-electron systems, that is, in order to calculate the potential adjustor, selfconsistent KS calculations of the cation and the anion have to be carried out. εHOMO[vext + vH + vx̅ c] and εLUMO[vext + vH + vx̅ c] denote the eigenvalues of the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO) of the neutral molecule, respectively. Any member of the set vx̅ c of exchange-correlation potentials can be chosen and N+ yields the same adjusted potentials vN− xc and vxc . The only requirement is consistency, that is, the eigenvalues and potential adjustors are to be calculated from the same member of the vx̅ c. In contrast to approaches (i) and (ii), the third approach to N+ calculate vN− xc and vxc is well-suited for LDA or GGA functionals for several reasons. First, for any finite system, different N+ potentials vN− xc and vxc are obtained, that is, the LDA and GGA exchange-correlation potentials exhibit a derivative discontinuity as they should. Second, it is known that, for molecules, energy differences EN0 − EN−1 and EN+1 − EN0 , which correspond 0 0 to the negatives of the first ionization potential and the first electron affinity, can be calculated with good accuracy by ΔSCF calculations using LDA or GGA functionals. A third advantage of approach (iii), which is beneficial not only for LDA and N+ GGA functionals but generally, is that vN− xc and vxc are obtained exclusively from quantities that are accessible within the KS formalism for a fixed integer electron number, despite the fact N+ that the potentials vN− xc and vxc are defined in the ensemble KS formalism. Finally, it should be kept in mind that eqs 10−13 are exact equations. Approximations are only introduced by the approximate functionals they are applied to. If the energies EN−1 and EN+1 in eqs 11 and 13 for the 0 0 N− potential adjustors Δxc and ΔN+ xc are not calculated selfconsistently but by evaluating EN−1 and EN+1 with orbitals of the 0 0 N-electron system, then this results in an approximate adjustment of exchange-correlation potentials suggested in refs 69 and 70. The evaluation of all required energies with the orbitals of the N-electron system means that relaxation effects upon ionization or upon addition of an electron are neglected.

where q is the, generally noninteger, number of electrons, and N is an integer. There are at least three different possibilities to N+ determine the potentials vN− xc and vxc . For the exact exchangecorrelation functional, all three approaches are equivalent. For approximate functionals, however, they lead to different sets of eigenvalues. The three possibilities are as follows: N+ (i) The potentials vN− xc and vxc are obtained by taking the functional derivatives of the exchange-correlation energy of the ensemble KS formalism directly according to eqs 8 and 9. This N+ is the obvious route to the potentials vN− xc and vxc , typically followed in standard KS and GKS calculations. For this approach to yield meaningful eigenvalue spectra, however, reliable approximations for the exchange-correlation functional of the ensemble KS formalism are required. Exchangecorrelation functionals within the local density approximation (LDA) or the generalized gradient approximation (GGA), for example, can be evaluated for any electron density including those integrating to noninteger electron numbers. Therefore, in practice, LDA and GGA functionals are often considered as approximate ensemble exchange-correlation functionals. This approximation, however, is an unphysical one, because the functional derivative of LDA and GGA energies with respect to the electron density does not exhibit an integer discontinuity. This means that the common practice to calculate LDA or GGA exchange-correlation potentials is fundamentally flawed and, therefore, questionable. (ii) For the case of finite nonperiodic electronic systems, the potential vN− xc can be obtained from any member of the set of exchange-correlation potentials vx̅ c of the KS formalism with a fixed integer electron number by simply adding a constant such that the potential vanishes at infinity. The justification is that N− the exact potential vxc vanishes far away from a finite nonperiodic electronic system.67,68 For the occupied eigenvalue spectra obtained from LDA or GGA functionals, this convention has the same outcome as approach (i) discussed N+ above. However, apart from the fact that only vN− xc but not vxc is accessible in this way, this approach is, again, questionable. This is due to the incorrect asymptotic behavior of semilocal functionals: While the exact exchange-correlation potential decays as 1/r, with r being the distance from the center of charge, LDA and GGA potentials decay exponentially and, hence, much too fast. In other words, the adjustment of the potential energy relies on a qualitatively wrong asymptotic behavior of the potential itself, which, for obvious reasons, is not advisable. N+ (iii) Following ref 52, the potentials vN− xc and vxc can be determined by well-defined potential adjustors, i.e., the N+ potentials vN− xc and vxc are constructed from any member of the set of exchange-correlation potentials vx̅ c according to vxcN −(r) = vxc̅ (r) + ΔNxc−[vxc̅ ]

3. ENERGETIC ADJUSTMENT OF EIGENVALUE SPECTRA The adjustment of the exchange-correlation potential induces a shift in the corresponding KS eigenvalue spectrum. The KS equations now read ⎡ ∇2 ⎤ + vext(r) + vH(r) + vxc̅ (r) + ΔNxc−[vxc̅ ]⎥φi(r) ⎢− ⎣ 2 ⎦ ⎡ ∇2 ⎤ = ⎢− + vs̅ (r) + ΔNxc−[vxc̅ ]⎥φi(r) ⎣ 2 ⎦

(10)

with the potential adjustor ΔNxc−[vxc̅ ]

=

E0N

−

E0N − 1

= [εi[vs̅ ] + ΔNxc−[vxc̅ ]]φi(r) − εHOMO[vext + vH + vxc̅ ]

(11)

= εi N −φi(r)

and according to vxcN +(r)

= vxc̅ (r) +

(13)

(14)

ΔNxc+[vxc̅ ]

if the exchange-correlation potential is adjusted with the potential adjustor ΔN− xc [vx̅ c] provided in eq 11. Then, the adjusted KS eigenvalues εN− are given by i

(12)

with the potential adjustor 4729


Article

Journal of Chemical Theory and Computation εi N − = εi[vs̅ ] + ΔNxc−[vxc̅ ]

⎡ ∇2 ⎤ + vext(r) + vH(r) + vxc̅ (r) + ΔNxc+[vxc̅ ]⎥φi(r) ⎢− ⎣ 2 ⎦

(15)

Eq 15 yields the same adjusted eigenvalues εN− for all members i are uniquely of the set vx̅ c. The adjusted eigenvalues εN− i defined, because the addition of a constant to a potential of the set vx̅ c leads to a corresponding change in the opposite direction in the potential adjustor ΔN− xc [vx̅ c], see eq 11, and therefore in summary has no effect on eq 15. With the help of eq 11 for ΔN− xc [vx̅ c], we obtain for the HOMO eigenvalue

⎡ ∇2 ⎤ = ⎢− + vs̅ (r) + ΔNxc+[vxc̅ ]⎥φi(r) ⎣ 2 ⎦ = [εi[vs̅ ] + ΔNxc+[vxc̅ ]]φi(r) = εi N +φi(r) (17)

where the adjusted KS eigenvalues

N− εHOMO = εHOMO[vs̅ ] + ΔNxc−[vxc̅ ]

are given by

εi N + = εi[vs̅ ] + ΔNxc+[vxc̅ ]

= εHOMO[vs̅ ] + E0N − E0N − 1 − εHOMO[vs̅ ] = E0N − E0N − 1

εN+ i

(18)

Again, eq 18 yields the same adjusted eigenvalues εN+ for all i members of the set vx̅ c. That means the adjusted eigenvalues N− εN+ i , like the eigenvalues εi , are uniquely defined. For the LUMO eigenvalue we obtain

(16)

Consequently, the HOMO eigenvalue equals the negative of the first IP given by EN0 − EN−1 0 , as it should. In the case of LDA or GGA functionals the negative of the adjusted HOMO eigenvalue equals the first ionization potential calculated by a ΔSCF calculation with the LDA or GGA functional. This means eq 16 holds true not only for the exact functional but also for LDA or GGA functionals, provided the exchangecorrelation potentials and, subsequently, the KS eigenvalues are correctly adjusted energetically. Technically, the adjustment of the KS eigenvalue spectra can be achieved by calculating the first ionization potential via a ΔSCF calculation and by then shifting the original eigenvalue spectrum such that the negative of the HOMO eigenvalue equals the first ionization potential. As already mentioned, the crucial point is that this is not just an ad hoc shift, but it is the formally correct adjustment of the KS eigenvalue spectrum for the case of approaching the integer electron number N from below, which is the case relevant for ionization. For the first ionization potential, the approach to determine ionization energies from energetically adjusted KS eigenvalues is equivalent to the well-known ΔSCF approach. However, having at hand an eigenvalue spectrum that is properly adjusted energetically, we propose to approximate higher ionization potentials by the negatives of the energetically lower occupied KS eigenvalues. At this stage, it is important to emphasize that, as mentioned above, it is well-known in the literature that DFT eigenvalue spectra that have been shifted to align the HOMO eigenvalue with the IP often yield good approximations to experimental valence-electron photoelectron spectra. This is particularly true when using standard global hybrid functionals within the GKS framework, since the inclusion of HF-exchange partially corrects for orbital self-interaction errors and introduces a stretching of the eigenvalue spectrum that mimics screening effects.36,37 A fundamentally new aspect of this work is that it provides a formal and rigorous justification for this shift of the eigenvalue spectrum, which is shown to result directly from the formally correct treatment of semilocal functionals in ensemble DFT. Of course, the potential adjustors approach can also be applied to the spectrum of unoccupied KS and GKS eigenvalues. If the exchange-correlation potential is adjusted with the potential adjustor ΔN+ xc [vx̅ c] from eq 13, then the KS equations read

N+ εLUMO = εLUMO[vs̅ ] + ΔNxc+[vxc̅ ]

= εLUMO[vs̅ ] + E0N + 1 − E0N − εLUMO[vs̅ ] = E0N + 1 − E0N

(19)

with eq 13 for ΔN+ xc [vx̅ c]. This means the LUMO eigenvalue equals the negative of first electron affinity given by EN+1 − EN0 . 0 In the case of LDA or GGA functionals, the negative of the adjusted LUMO eigenvalue equals the first electron affinity calculated by a ΔSCF calculation with the LDA or GGA functional. Analogous to the case of ionization, this means eq 19 holds true not only for the exact functional but also for LDA or GGA functionals provided the exchange-correlation potentials and subsequently the KS eigenvalues are correctly adjusted energetically. It suggests itself to approximate higher electron affinities by the negatives of energetically higher KS eigenvalues εN+ i . In summary, two KS eigenvalues spectra containing the eigenvalues εN− and εN+ i i , respectively, are associated with the charged excitation energies, i.e., the quasiparticle energies, of an electronic system. The occupied spectrum of the eigenvalues εN− is related to ionization, while the unoccupied spectrum of i the eigenvalues εN+ is related to electron attachment. The i approximate full quasiparticle spectrum is then obtained by combining the occupied eigenvalues from the first spectrum with the unoccupied eigenvalues from the second. For the case of the exact functional, the HOMO and LUMO eigenvalues of this combined spectrum are exact quasiparticle energies, that is they yield the exact first ionization potential and electron affinity. All other eigenvalues of this combined spectrum can be interpreted as approximate quasiparticle energies. The fact that two KS eigenvalue spectra are associated with an electronic system that need to be combined for an approximate quasiparticle spectrum can be rationalized based on the ensemble KS formalism, more precisely from the presence of the derivative discontinuity in the exchange-correlation functional. The correct energetic adjustment of LDA and GGA exchange-correlation potentials introduces a derivative discontinuity for finite systems and, therefore, leads to two separate KS eigenvalue spectra for occupied and unoccupied states. As mentioned in the Introduction, such approximate quasiparticle spectra shall be denoted potential adjusted KS (paKS) or GKS (paGKS) spectra. 4730


Article


Figure 1. Deviations of paB3LYP IPs from experimental IPs80,81 for different basis sets.

4. IONIZATION POTENTIALS OF SMALL- AND MEDIUM-SIZED MOLECULES 4.1. Test Sets for Ionization Potentials. We tested the above-described approach to calculate IPs via paKS or paGKS eigenvalues for a test set of 113 IPs of 12 small- and mediumsized molecules, called TS12, in the following. The molecules contained in TS12 are N2(5), F2(6), CO2(8), H2O(4), NH3(4), SiH4(4), methane(4), ethylene(5), benzene(15), naphthalene(17), anthracene(23), and naphthacene(18). (The number of considered IPs is given in parentheses after the name of each molecule.) The test set TS12 is a subset of the GW27 set of van Setten et al.71 A second test set, which is used for comparisons with G0W0 methods and is called TS13, differs from TS12 in that it does not contain naphthacene but additionally considers formaldehyde and butadiene. Moreover, TS13 does not contain the fourth IP of N2 and the first IP of CO2 because for these IPs the G0W0 calculations based on PBE orbitals and eigenvalues did not converge for the convergence factor η = 0.001 hartree which is the default value in the employed G 0 W 0 implementation71 of the Turbomole program package, see below for further details of the G0W0 calculations. That is, TS13 contains 109 IPs of the 13 molecules N2(4), F2(6), CO2(7), H2O(4), NH3(4), SiH4(4), methane(4), ethylene(5), benzene(15), naphthalene(17), anthracene(23), formaldehyde(5), and butadiene(11). The geometries of the molecules of TS12 were adopted from the Supporting Information of van Setten et al.,71 while the geometries for the test TS13 were obtained by geometry optimizations with the B3LYP functional72,73 and def2-TZVP basis sets.74 Further details on the geometries can be found in the SI. The small differences in geometries (bond lengths typically deviate by less than 2%) of the two test sets only lead to negligible changes in the calculated IPs showing that the geometry dependence of the IPs is weak. 4.2. Computational Details. All calculations were performed with the Turbomole package.75 The resolution-ofidentity (RI-J) approximation was applied to speed up the computation of the Coulomb integrals.76 In all cases basis sets of the def2 type contained in the Turbomole basis set library77 and corresponding auxiliary basis sets for the RI74 were chosen. If not stated otherwise we used for paKS and paGKS calculations the def2-TZVP basis set. Our approach was compared to G0W0 methods. To this end, G0W0 calculations were carried out with the recent G0W0 implementation in Turbomole71 with the def2-TZVPP basis set.77 This basis set was found to provide well converged quasiparticle energies.71 The convergence factor η was set to 0.001 hartree, and all singlet TD-DFT excitations have been taken into account to calculate the screened Coulomb

interaction W0. The G0W0 quasiparticle energies zn have been calculated using the equation zn = εn + ⟨φn|Σ[G KS](εn) − vxc|φn⟩

(20)

That is we did not use the linearized quasiparticle equation zn = εn + Zn⟨φn|Σ[G KS](εn) − vxc|φn⟩

(21)

with the normalization factor Zn ⎡ ⎤−1 ∂Σ[G KS](E) Zn = ⎢1 − ⟨φn| |E = εn |φn⟩⎥ ⎣ ⎦ ∂E

(22)

which is obtained by Taylor-expanding Σ[GKS](E) around εn up to first order and which is evaluated by default in Turbomole.71 The linearized formula 21, however, lead to distinctively worse results than the ordinary formula 20. The G0W0 quasiparticle energies of ref 23 calculated with the program FHI-AIMS,78,79 to which we compare in Section 4.6, also were obtained with formula 20. The calculations in FHIAIMS employ a basis set of numeric atom-centered orbitals of rank tier 4. For a better comparison with these results the paKS and paGKS calculations were carried out with the larger def2QZVP basis set.77 4.3. Basis Set Dependence of paKS and paGKS Ionization Potentials. Figure 1 shows the deviation ΔIP ≡ IP(calc) − IP(exp) of IPs calculated with the paGKS approach paB3LYP based on the B3LYP functional72,73 to experimental values80,81 for water and benzene for different basis sets.77 We observe that the basis set def-SVP, a split valence basis set including polarization functions, clearly gives nonconverged results, whereas basis sets of TZVP quality and better yield very similar results. A similar behavior was also found for higher IPs. The same picture emerges if all molecules of the TS12 set are considered, of which the statistical analysis is summarized in Table 1: TZVP, TZVPP, and QZVP basis sets yield very similar results; the SVP results, however, deviate somewhat. The Table 1. Basis-Set Dependence of the Accuracy of paB3LYP IPs for the TS12 Test Set with Respect to Experimental Values, See the SIa TS12

def2SVP

def2TZVP

def2TZVPP

def2QZVP

mean signed error mean absolute error root mean squared error

−0.56 0.64 0.86

−0.49 0.60 0.81

−0.49 0.60 0.81

−0.49 0.60 0.81

a

The mean values are taken over all molecules, and all IPs that are below. All values are in eV.

4731


Article

Journal of Chemical Theory and Computation statistical measures are defined as follows (N is the number of considered ionization potentials): mean signed error =

mean absolute error =

root mean squared error =

N

1 N

− IP(exp) ∑ IP(calc) i i

1 N

− IP(exp) | ∑ |IP(calc) i i

i=1 N i=1

1 N

N 2 − IP(exp)) ∑ (IP(calc) i i i=1

We observed a similar basis set dependence for paPBE and paLDA calculations, i.e., paKS calculations based on PBE or LDA. Because converged results are obtained with the def2TZVP basis set, it was used for the rest of the paper if not stated otherwise. 4.4. Dependence of paKS and paGKS Ionization Potentials on the Exchange-Correlation Functional. Next, we investigate the influence of different exchangecorrelation functionals on paKS or paGKS IPs. We compare the LDA, the GGA functional due to Perdew, Burke, and Ernzerhof (PBE), 82 and the global hybrid functionals B3LYP72,73 and PBE0.83 To begin with, we consider three typical molecules of our test set: water, the nitrogen molecule, and benzene. The individual results are depicted in Figure 2. In case of water the best results both for the lowest three IPs and the energetically much higher fourth IP are obtained by applying the hybrid functional PBE0. The hybrid functional B3LYP and the semilocal functional PBE, however, deliver results of a similar quality except for the fourth IP, which is underestimated slightly by paB3LYP and by more than 1 eV by paPBE compared to the experimental value.80 The LDA exhibits larger errors than PBE. Nevertheless, the errors for the first three IPs are still smaller than 1 eV. In summary, only paPBE0 and paB3LYP are able to describe all four IPs of the water molecule within an error bar of ±0.5 eV. Considering the N2 molecule, the lowest three paPBE and paLDA IPs lie within ±0.3 eV of the experimental data,84 see Figure 2. paPBE0 and paB3LYP yield similarly accurate results except for the third IP, which exhibits a deviation of about 0.7 and 0.5 eV, respectively, compared to the experimental value. It is well-known that Hartree−Fock eigenvalues of N2 are in the wrong order with respect to the first and second IP.84 This problem does not occur for the KS and GKA spectra employed here. In case of the fourth IP, which corresponds to an excitation from the 2a1g orbital, all functionals fail more or less. The smallest error, which, however, amounts to roughly 2.5 eV, is obtained by paPBE0. In benzene (experimental values from ref 81) paLDA overestimates low IPs by approximately 1 eV, while higher IPs are underestimated by roughly 1 eV. The curve of the paPBE results resembles the corresponding paLDA curve but is shifted by nearly 1 eV to lower energies. Thus, the first paPBE IP, which is identical to the ΔSCF value, reproduces very well the experimental value but systematically larger errors are obtained for higher IPs. In contrast, the B3LYP functional and even more so the PBE0 functional appears to be more unbiased: most IPs are underestimated but the errors, even for higher IPs, always remain smaller than 1 eV. For the first five IPs the deviations from the experimental values are smaller than 0.5 eV.

Figure 2. Deviations of paKS and paGKS IPs from experimental IPs80,81,84 for different density functionals.

In Table 2 statistical errors of the paLDA, paPBE, paB3LYP, and paPBE0 IPs of the TS12 test set with respect to experimental values, see SI, are displayed. In agreement with the above findings for N2, H2O, and benzene, paPBE0 yields the best match with experimental data in general, followed by paB3LYP, paPBE, and paLDA. If all IPs up to 40 eV are considered, the mean absolute error of PBE0 is 0.45 eV, while paB3LYP, paLDA, and PBE yield mean absolute errors of 0.60, 0.91, and 1.07 eV, respectively. Interestingly, if one considers the mean signed errors for all three considered ranges of IPs, it shows that the PBE functional yields a significantly larger error than LDA, although the mean absolute errors are comparable. This is a manifestation of the fact that paPBE systematically underestimates IPs, which we already observed above for benzene. paLDA tends to overestimate the first ionization potentials and to underestimate the higher ones. This can be clearly seen from the decline of the mean signed error from the first IP to all IPs up to 40 eV. The PBE0 functional clearly yields the smaller errors; however, it also tends to underestimate IPs, especially higher ones. 4.5. Comparison with G0W0. In this section, paB3LYP, paPBE0, and paPBE IPs are compared to G0W0 quasiparticle 4732


Article

Journal of Chemical Theory and Computation Table 2. Accuracy of the IPs of the TS12 Test Set with Respect to Experimental Values, See the SI, for Different Density Functionalsa TS12

paLDA

paPBE

First IPs 0.12 −0.21 0.30 0.28 0.41 0.34 First Three IPs mean signed error 0.01 −0.31 mean absolute error 0.48 0.50 root mean squared error 0.67 0.63 All IPs up to 40 eV mean signed error −0.71 −1.00 mean absolute error 0.91 1.07 root mean squared error 1.22 1.34 mean signed error mean absolute error root mean squared error

paB3LYP

paPBE0

−0.13 0.21 0.27

−0.10 0.16 0.22

−0.14 0.37 0.47

−0.06 0.30 0.42

−0.49 0.60 0.81

−0.32 0.45 0.67

a The mean values are taken over all molecules and (i) only the first IPs, (ii) the first three different IPs, and (iii) all IPs that are below 40 eV. All values are in eV.

energies. Three G0W0 approaches with different DFT starting points are included in the comparison, that is, G 0 W 0 calculations based upon PBE, B3LYP, and PBE0 denoted as G0W0@PBE, G0W0@B3LYP, and G0W0@PBE0, respectively. Considering again the water molecule first, see Figure 3, we note that for the first three IPs paB3LYP and paPBE0 deliver the best results with respect to experimental data. The IPs of the G0W0 approaches are similarly accurate, with G0W0@PBE being the best for the first two IPs, whereas G0W0@PBE0 turns out to be the best for the third IP. For the fourth IP, the highest valence IP, the methods yield significantly different results. G0W0@PBE overestimates this IP by 2 eV, G0W0@B3LYP by 0.5 eV, and G0W0@PBE0 by 0.2 eV, while paB3LYP and paPBE0 underestimate it by 0.7 and 0.5 eV, respectively. In summary paPBE0 delivers the best result closely followed by paB3LYP, G0W0@PBE0, and G0W0@B3LYP. For the nitrogen molecule findings differ from those of water. The G0W0 methods G0W0@PBE, G0W0@B3LYP, and G0W0@ PBE0 deliver very similar and accurate results for the first three IPs (see Figure 3). The first paPBE0 IP is somewhat too high, whereas the second IP is in perfect agreement with the experimental value. The first two paB3LYP IPs deviate slightly more from the experimental values than the G0W0 IPs. For the third IP paB3LYP and paPBE yield almost identical values being about 0.5 eV too high. The mean absolute deviation of the G0W0 IPs is a bit smaller than the corresponding mean absolute deviation of the paB3LYP IPs and less so of the paPBE0 IPs. All methods are able to describe the first three IPs within a maximum deviation of 0.6 eV. For the benzene molecule, all three G0W0 methods behave almost identical for IPs up to 20 eV. Overall, the three G0W0 methods tend to be slightly more accurate in that region than paB3LYP and paPBE0 with PBE0 being a tiny bit more accurate than paB3LYP. In contrast, the two valence IPs over 20 eV are much better described by the paB3LYP and the paPBE0 approach. While G0W0 based on PBE is the best G0W0 method for the IP at roughly 22.5 eV, G0W0 based on B3LYP is the best G0W0 method for the IP at roughly 25.8 eV. In Table 3 the IPs for the test set TS13 (13 different molecules with a total of 109 IPs) obtained from paPBE, paB3LYP, paPBE0, G0W0@PBE, G0W0@B3LYP, and G0W0@ PBE0 are compared against experimental reference data. (The

Figure 3. Deviations of paB3LYP, paPBE0, G0W0@PBE, G0W0@ B3LYP, and G0W0@PBE0 IPs from experimental IPs.80,81,84

results for the individual molecules can be found in the SI.) For the first IP, the approaches paB3LYP, paPBE, and G0W0@PBE yield results of similar accuracy. paPBE0 exhibits slightly smaller errors, whereas the G0W0@B3LYP and G0W0@PBE0 errors are marginally larger. If the three first IPs of each molecule are considered, then paPBE0 IPs are most accurate, closely followed by paB3LYP, G0W0@PBE, and G0W0@B3LYP IPs, whereas G0W0@PBE0 and paPBE IPs are slightly less accurate. For all six methods the mean absolute error is clearly below 0.5 eV. If, for each molecule, all IPs up to 40 eV are taken into account, then again paPBE0 IPs are most accurate, followed by G0W0@PBE, paB3LYP, and G0W0@B3LYP IPs, while G0W0@PBE0 and paPBE IPs are somewhat less accurate. If all IPs up to 40 eV are considered, both the mean absolute and the root mean squared error are larger than for the first three IPs. For the paPBE0 the mean absolute error increases from 0.25 to 0.34 eV. The paB3LYP, G0W0@PBE, G0W0@ B3LYP, and G0W0@PBE0 mean absolute errors increase from roughly 0.3 to roughly 0.5 eV. For paPBE IPs the mean absolute error increases from about 0.4 to 0.9 eV. In summary, 4733


Article


Table 3. Deviations of paPBE, paPBE0, paB3LYP, G0W0@PBE, G0W0@B3LYP, and G0W0@PBE0 IPs of the TS13 Test Set Compared to Experimental Values, See the SI for Individual IPsa TS13

paPBE

paB3LYP


−0.16 0.21 0.26

−0.09 0.20 0.23


−0.30 0.42 0.56

−0.13 0.29 0.37


−0.82 0.87 1.08

−0.34 0.45 0.58

paPBE0 HOMO IPs −0.07 0.17 0.20 First Three Different IPs −0.06 0.25 0.34 All IPs up to 40 eV −0.19 0.34 0.48

G0W0@PBE

G0W0@B3LYP

G0W0@PBE0

−0.01 0.20 0.25

−0.07 0.25 0.29

−0.11 0.23 0.29

−0.03 0.28 0.41

−0.14 0.30 0.43

−0.19 0.32 0.46

−0.15 0.39 0.55

−0.25 0.47 0.79

−0.35 0.52 1.02

The mean values are taken over all molecules, and (i) only the first IPs, (ii) the first three different IPs, and (iii) all IPs that are below 40 eV. All values are in eV. a

single Gaussian functions obtained from the genetic algorithm is shown at the top of Figures 4−8. We then used the same

the paPBE0 method yields the most accurate results for all considered ranges of IPs. The next accurate methods are G0W0@PBE and paB3LYP, which are just slightly worse, followed by G0W0@B3LYP and G0W0@PBE0. However, the difference in accuracy between the best four approaches is moderate, in particular for the first two considered ranges of IPs. paPBE and, to a lesser extent, G0W0@PBE0 yield slightly less accurate results. The potential adjustors for the molecules of the test test TS13 in the PBE, B3LYP, and PBE0 case can be found in the SI. Furthermore, IPs obtained directly as negatives of the nonadjusted PBE, B3LYP, and PBE0 eigenvalues are listed in the SI. The deviations of these IPs from experimental values are about an order of magnitude larger than that of IPs obtained with the help of potential adjustors. 4.6. Medium-Sized Electron Acceptor Molecules. Finally we applied the paKS and paGKS methods to somewhat larger molecules. To that end, a set of electron acceptor molecules was considered whose IPs were investigated in ref 23, where experimental photoelectron spectra (PES)85 are compared with spectra calculated by various G0W0 variants. Here, we compare the results of paPBE, paPBE0, and paB3LYP with the G0W0 spectra calculated in ref 23 as well as to the spectrum obtained from the eigenvalues of a tuned long-range corrected hybrid functional. As for the latter, we followed the procedure proposed in ref 40 for the accurate prediction of outer-valence photoelectron spectra, that is, we included a fixed fraction of 20% HF-exchange in the short-range part of the long-range corrected hybrid functional LC-ωPBEh86 and then tuned the range-separation parameter to match the IPtheorem.38,40,87 In order to allow for a better visual comparison of the quasiparticle spectra, we developed a procedure in which we aim to reproduce the experimental spectra as well as possible by a representation with Gaussian functions with variable heights and broadenings. This is supposed to mimic the different ionization cross sections, Franck−Condon factors, and vibrational broadenings of the experimental spectrum, which we do not calculate explicitly in this work. To achieve this goal, we first determined the number and position of the quasiparticle states in the experimental spectrum and then used a genetic algorithm to find the optimal peak heights and broadenings for the single Gaussian functions representing the corresponding peaks. The experimental spectrum as well as the best-fit from

Figure 4. Experimental gas-phase photoelectron spectrum of phenazine85 and best fit from single Gaussians compared to calculated spectra. Peak heights and broadenings obtained from the best-fit to experiment. Errors of the different methods (as defined in the text) are provided in brackets. The dashed blue lines represent the position of the first IP from a CCSD(T) calculation.88

Gaussian peak heights and broadenings, strictly following the ordering of the quasiparticle energies as provided by the different theoretical methods, to determine the spectra from the calculated quasiparticle energies obtained from the GW calculations performed in ref 23 as well as the LC-ωPBEh,86 paPBE, paPBE0, and paB3LYP calculations performed in this work. This procedure allows for a clear and direct visual comparison of the calculated spectra with experiment. 4734


Article


estimate the HOMO IP only by 0.04 and 0.15 eV (second energy in brackets), respectively. 4.6.2. Benzoquinone. Benzoquinone is known to be a notorious case for which the prediction of the first IP and the outer-valence photoelectron spectrum has been known to be particularly problematic.22,23 This has to do with the fact that the HOMO of benzoquinone is an n-type orbital that includes the free electron pair of the oxygen atoms. n-Type orbitals typically suffer from larger self-interaction errors than π-type orbitals, and finding a balanced treatment of the errors in both type of orbitals has been found to be problematic for semilocal DFT34,36,39 as well as for G0W023,26,31 approaches based on a semilocal starting point. The experimental ionization spectrum of benzoquinone is shown at the top of Figure 5: two double peaks are followed by

In an attempt to further analyze the deviation from the simulated spectra from the best Gaussian fit, two different error measures are provided in the plots. First, in order to quantify the accuracy of the relative quasiparticle spectrum, we calculate the integral over the magnitude of the difference between the various simulated spectra and the best-fit to the experiment, respectively, with the HOMO peak being set to zero in all spectra. In Figures 4−8, the resulting energy is given by the first number in brackets. The corresponding plots of the difference spectra are provided in the SI. With the second number in brackets we also provide the deviation of the HOMO energy from the experimental value, which measures the accuracy of the absolute ionization energies provided by the method at hand. However, it is important to realize that, by directly comparing calculated vertical IPs with experiment, we neglect the influence of zero-point vibrational corrections, relaxation energies, Franck−Condon integrals, and temperature effects. Hence, we also included theoretical reference values for the vertical first IP, calculated from coupled cluster singles douples perturbative triples (CCSD(T)),88 represented by the dashed blue lines in Figures 4−8. 4.6.1. Phenazine. The experimental PES of phenazine (see Figure 4) can be divided into two regions: the first one starts at the low-energy edge with a well-resolved double-peak followed by two peaks of which the first one is rather broad. Then at approximately 11 eV, after a small gap a second absorption region appears with several broad, poorer resolved peaks. The paPBE data lead to different assignments of the lower IPs than the paB3LYP and paPBE0 results. paPBE suggests that the first absorption peak and its shoulder might be assigned to two different IPs. This, however, is in contrast to all the other approaches shown in Figure 4. It is more probable that the shoulder of the first sharp peak in the PES spectrum is due to a vibrational fine structure. Assuming this, the paB3LYP and paPBE0 spectra reproduce the experimental PES very well, both qualitatively and quantitatively. The clear differences between the paB3LYP/paPBE0 and paPBE spectra and the corresponding assignment of IPs, however, indicate the importance of using a hybrid functional for this kind of molecules. While the paPBE spectrum clearly suffers from an incorrect orbital ordering due to large orbital self-interaction errors,34,35 the paB3LYP and LC-ωPBEh(IP) approaches yield spectra in quite good agreement with experiment. In fact, the good performance of standard global hybrid functionals such as B3LYP for relative outer-valence photoelectron spectra has been noted and discussed earlier in the literature.34−37,89−91 It is therefore not surprising that the combination of such functionals with the rigorous potential adjustors approach leads to spectra that compare very well to experimental outer-valence photoelectron spectra, even for challenging systems such as phenanzine. This is also reflected in the integral over the magnitude of the difference to the experimental spectrum (first energy in brackets) calculated after shifting the spectra such that the HOMO IPs are aligned. Values for this integral of 0.67 and 0.74 eV for paPBE0 and paB3LYP spectra show that the latter methods indeed exhibit an (relative) accuracy similar to that of the most accurate G0W0 methods in representing the best-fit spectrum. The only method that stands out is paPBE, which shows a much larger error (2.15 eV). With respect to the absolute IPs, the most accurate methods in direct comparison to experiment are G0W0@LC-ωPBE(IP) and LC-ωPBEh(IP), which under-

Figure 5. Experimental gas-phase photoelectron spectrum of benzoquinone92 and best fit from single Gaussians compared to calculated spectra. Peak heights and broadenings obtained from the best-fit to experiment. Errors of the different methods (as defined in the text) are provided in brackets. The dashed blue lines represent the position of the first IP from a CCSD(T) calculation.88

a broad band starting approximately at 13 eV. In this case, the paB3LYP/paPBE0 and paPBE approaches yield somewhat more similar spectra than in the case of phenazine; in particular the assignment of the IPs is identical for both methods. However, the paPBE approach clearly underestimates the lowest IP by 0.88 eV, and the energetic distance between the first and the second resolved peaks is too large compared to experiment. paB3LYP and paPBE0 are doing much better again, which is also reflected in smaller difference integrals (0.59 and 0.76 eV) than in all other methods. The LCωPBEh(IP) and G0W0@LC-ωPBE(IP) approaches predict no clear separation of two peaks in the low-energy region, while the remaining G0W0 methods are able to describe this splitting somewhat better. G0W0+SOSEX@PBE shows good perform4735


Article

Journal of Chemical Theory and Computation ance for the higher energy peaks and the best HOMO energy of all methods, but the deviations for the broad peak at lower energies lead to a rather high difference integral of 1.35 eV. The overall shape of the experimental spectrum is best reproduced by paB3LYP and paPBE0. 4.6.3. Cl4-Benzoquinone. Substituting the hydrogen by chlorine atoms in benzoquinone essentially introduces an additional peak in the experimental PES at approximately 12− 13 eV. The two low energy bands are slightly shifted to lower energies. In this case the best agreement between experiment and theory is provided by the G0W0+SOSEX@PBE approach (1.06 eV), followed by G0W0@LC-ωPBE(IP) (1.31 eV), G0W0@PBE0 (1.46 eV), and LC-ωPBEh(IP) (1.47 eV). The latter two also provide the best results for the HOMO energy. The paB3LYP (1.89 eV) and paPBE0 (1.76 eV) spectra are, in this case, less accurate than for the molecules discussed above but still of similar accuracy as G0W0@PBE (1.86 eV) and, again, significantly better than paPBE (2.30 eV). Interestingly the onset of the paB3LYP spectrum, i.e., the position of the first IP, perfectly matches the CCSD(T) value. The paPBE method deviates more from the experiment, making an assignment of peaks difficult. 4.6.4. Maleic Anhydride. For the case of maleic anhydride, the experimental PES can be reproduced qualitatively by all methods shown in Figure 7, except paPBE. The paPBE method swaps the order of the first IPs, which gives rise to a different line shape in the low energy region of the spectrum. Concerning the quantitative agreement of peak positions with experiment the G0W0@PBE (0.88 eV) and G0W0@PBE0 (0.90

Figure 7. Experimental gas-phase photoelectron spectrum of maleic anhydride94 and best fit from single Gaussians compared to calculated spectra. Peak heights and broadenings obtained from the best-fit to experiment. Errors of the different methods (as defined in the text) are provided in brackets. The dashed blue lines represent the position of the first IP from a CCSD(T) calculation.88

eV) spectra are the best, closely followed by the paB3LYP and paPBE0 spectra (1.21 eV each). 4.6.5. Nitrobenzene. For nitrobenzene, see Figure 8, similar conclusions as for maleic anhydride can be drawn: all methods except paPBE reproduce the experimental PES reasonably well. For the paPBE approach again the lowest IPs are especially problematic; they are predicted too low in energy and in incorrect order. The absolute position of the first IP is best described by paB3LYP and paPBE with a deviation from the experimental value of only 0.02 eV. The separation of the first IPs is obtained most accurate with the G0W0@LC-ωPBE(IP), the G0W0+SOSEX@PBE0, and the G0W0@PBE0 approach, while the first two peaks are obtained too narrow to each other by paB3LYP, paPBE0, and G0W0@PBE. The best relative agreement with the experimental PES is found for LCωPBEh(IP) (0.72 eV) and G0W0+SOSEX@PBE (0.73 eV), followed by G0W0@PBE0 (1.14 eV) and paPBE0 (1.25 eV).

5. CONCLUSIONS We showed that the paKS and paGKS methods presented here, in particular the paPBE0 and the paB3LYP approach, enable an accurate yet computationally efficient calculation of the first as well as higher ionization potentials and thus of photoelectron spectra. For the test set TS13 of 109 IPs of 13 small- and medium-sized molecules we compared the accuracy of the paPBE, paB3LYP, and paPBE0 methods with G0W0@PBE, G0W0@B3LYP, and G0W0@PBE0 methods, i.e., G0W0 methods based on PBE, B3LYP, and PBE0 orbitals and eigenvalues, respectively, and found that the accuracy can be ranked as

Figure 6. Experimental gas-phase photoelectron spectrum of Cl4benzoquinone93 and best fit from single Gaussians compared to calculated spectra. Peak heights and broadenings obtained from the best-fit to experiment. Errors of the different methods (as defined in the text) are provided in brackets. The dashed blue lines represent the position of the first IP from a CCSD(T) calculation.88 4736


Article


spectra of comparable, in some cases even superior, accuracy than those obtained from more costly G0W0 approaches. The paPBE approach yielded spectra of a somewhat lower accuracy compared to the other methods. The PBE method and subsequently also the paPBE and G0W0@PBE approaches suffer from large orbital self-interaction errors and, therefore, often predict wrong orbital orderings and qualitatively wrong relative eigenvalue spectra. In accordance with earlier findings, hybrid functionals such as paPBE0 and paB3LYP give much more accurate orbital orderings and, therefore, also a generally more reliable prediction of outer-valence photoelectron spectra. In summary, we showed that the paKS and paGKS methods presented here, in particular the paPBE0 and paB3LYP methods, can be an attractive alternative to common G0W0 methods for the calculation of molecular IPs and photoelectron spectra, especially for large molecules, because they yield similar accuracy with much less computational effort. In addition, they can be carried out with every KS/GKS code without the need to make any changes in the code. The paKS and paGKS methods can be used to calculate electron affinities in a completely analogous way. The accuracy of such electron affinities remains to be investigated.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.7b00490. Geometries of the molecules of the test sets TS12 and TS13, ionization potentials obtained directly from unshifted PBE, B3LYP, and PBE0 eigenvalues for the molecules of test set TS13, potential adjustors for all molecules of the test set TS13 for all considered functionals, figures with the individual paB3LYP, paPBE0, G0W0@PBE, G0W0@B3LYP, and G0W0@ PBE0 ionization energies of all molecules of the test set TS13, and plots of the difference between the best fit to the experimental spectra and the calculated spectra after aligning the ionization potentials of the HOMOs for the considered medium-sized acceptor molecules (PDF)

Figure 8. Experimental gas-phase photoelectron spectrum of nitrobenzene95 and best fit from single Gaussians compared to calculated spectra. Peak heights and broadenings obtained from the best-fit to experiment. Errors of the different methods (as defined in the text) are provided in brackets. The dashed blue lines represent the position of the first IP from a CCSD(T) calculation.88

paPBE0 > paB3LYP > G0W0@PBE ≈ G0W0@B3LYP > G0W0@PBE0 > paPBE

The paPBE0 and to a lesser degree also the paB3LYP method turn out to be more accurate than commonly used G0W0 methods for the test set of IPs considered in this work. paKS or paGKS methods just require two self-consistent KS or GKS calculations and a simple shifting of the spectrum of occupied eigenvalues such that the HOMO eigenvalue equals the negative of the first IP calculated by ΔSCF and, thus, are computationally much more efficient than G0W0 methods. Despite their conceptual simplicity and computational efficiency, paKS and paGKS methods have a solid formal justification in the ensemble KS formalism as discussed in sections 2 and 3 and, in more detail, in ref 52. The good performance of paKS or paGKS approaches can be traced back to the fact that total energy differences from KS or GKS methods are much more accurate than the energetic position of KS or GKS eigenvalues.22 By shifting the eigenvalue spectrum, which corresponds to an energetic adjustment of the exchangecorrelation potential as required by the ensemble KS formalism, the accuracy of KS or GKS total energy differences is transferred to the KS or GKS eigenvalue spectrum. We then calculated photoelectron spectra for five mediumsized electron acceptor molecules with the paPBE, paPBE0, and paB3LYP methods and compared with experimental spectra and with spectra obtained from various G0W0 variants as well as from eigenvalues of a tuned long-range corrected hybrid functional. Again the paB3LYP and paPBE0 spectra yielded

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Noa Marom: 0000-0002-1508-1312 Andreas Görling: 0000-0002-1831-3318 Funding

This work was supported by the Deutsche Forschungsgemeinschaft (DFG) through the Cluster of Excellence “Engineering of Advanced Materials” (www.eam.uni-erlangen.de) at the Friedrich-Alexander-Universität Erlangen-Nürnberg. Work at Carnegie Mellon University was supported by the National Science Foundation (NSF) Division of Materials Research through grant DMR1554428. Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Binggeli, N.; Chelikowsky, J. R. Photoemission Spectra and Structures of Si Clusters at Finite Temperature. Phys. Rev. Lett. 1995, 75, 493−496.

4737


Article

Journal of Chemical Theory and Computation (2) Hill, I.; Kahn, A.; Cornil, J.; dos Santos, D.; Brédas, J. Occupied and unoccupied electronic levels in organic π-conjugated molecules: comparison between experiment and theory. Chem. Phys. Lett. 2000, 317, 444−450. (3) Coropceanu, V.; Malagoli, M.; da Silva Filho, D. A.; Gruhn, N. E.; Bill, T. G.; Brédas, J. L. Hole- and Electron-Vibrational Couplings in Oligoacene Crystals: Intramolecular Contributions. Phys. Rev. Lett. 2002, 89, 275503. (4) Spataru, C. D.; Ismail-Beigi, S.; Benedict, L. X.; Louie, S. G. Excitonic Effects and Optical Spectra of Single-Walled Carbon Nanotubes. Phys. Rev. Lett. 2004, 92, 077402. (5) Stenuit, G.; Castellarin-Cudia, C.; Plekan, O.; Feyer, V.; Prince, K. C.; Goldoni, A.; Umari, P. Valence electronic properties of porphyrin derivatives. Phys. Chem. Chem. Phys. 2010, 12, 10812− 10817. (6) Dauth, M.; Körzdörfer, T.; Kümmel, S.; Ziroff, J.; Wiessner, M.; Schöll, A.; Reinert, F.; Arita, M.; Shimada, K. Orbital Density Reconstruction for Molecules. Phys. Rev. Lett. 2011, 107, 193002. (7) Puschnig, P.; Reinisch, E.-M.; Ules, T.; Koller, G.; Soubatch, S.; Ostler, M.; Romaner, L.; Tautz, F. S.; Ambrosch-Draxl, C.; Ramsey, M. G. Orbital tomography: Deconvoluting photoemission spectra of organic molecules. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 84, 235427. (8) Kronik, L.; Kümmel, S. Gas-Phase Valence-Electron Photoemission Spectroskopy Using Density Functional Theory. Top. Curr. Chem. 2014, 347, 137−191. (9) Gallandi, L.; Körzdörfer, T. Long-Range Corrected DFT Meets GW: Vibrationally Resolved Photoelectron Spectra from First Principles. J. Chem. Theory Comput. 2015, 11, 5391−5400. (10) Puschnig, P.; Boese, A. D.; Willenbockel, M.; Meyer, M.; Lüftner, D.; Reinisch, E. M.; Ules, T.; Koller, G.; Soubatch, S.; Ramsey, M. G.; Tautz, F. S. Energy Ordering of Molecular Orbitals. J. Phys. Chem. Lett. 2017, 8, 208−213. (11) O’Boyle, N. M.; Campbell, C. M.; Hutchison, G. R. Computational Design and Selection of Optimal Organic Photovoltaic Materials. J. Phys. Chem. C 2011, 115, 16200−16210. (12) Bérubé, N.; Gosselin, V.; Gaudreau, J.; Côté, M. Designing Polymers for Photovoltaic Applications Using ab Initio Calculations. J. Phys. Chem. C 2013, 117, 7964−7972. (13) Rutledge, L. R.; McAfee, S. M.; Welch, G. C. Design and Computational Characterization of Non-Fullerene Acceptors for Use in Solution-Processable Solar Cells. J. Phys. Chem. A 2014, 118, 7939− 7951. (14) Cole, J. M.; Low, K. S.; Ozoe, H.; Stathi, P.; Kitamura, C.; Kurata, H.; Rudolf, P.; Kawase, T. Data mining with molecular design rules identifies new class of dyes for dye-sensitised solar cells. Phys. Chem. Chem. Phys. 2014, 16, 26684−26690. (15) Ornso, K. B.; Garcia-Lastra, J. M.; De La Torre, G.; Himpsel, F. J.; Rubio, A.; Thygesen, K. S. Design of two-photon molecular tandem architectures for solar cells by ab initio theory. Chem. Sci. 2015, 6, 3018−3025. (16) Hachmann, J.; Olivares-Amaya, R.; Jinich, A.; Appleton, A. L.; Blood-Forsythe, M. A.; Seress, L. R.; Roman-Salgado, C.; Trepte, K.; Atahan-Evrenk, S.; Er, S.; Shrestha, S.; Mondal, R.; Sokolov, A.; Bao, Z.; Aspuru-Guzik, A. Lead candidates for high-performance organic photovoltaics from high-throughput quantum chemistry - the Harvard Clean Energy Project. Energy Environ. Sci. 2014, 7, 698−704. (17) Hybertsen, M. S.; Louie, S. G. Electron correlation in semiconductors and insulators: Band gaps and quasiparticle energies. Phys. Rev. B: Condens. Matter Mater. Phys. 1986, 34, 5390−5413. (18) Hedin, L. New Method for Calculating the One-Particle Green’s Function with Application to the Electron-Gas Problem. Phys. Rev. 1965, 139, A796−A823. (19) Ortiz, J. V. Electron propagator theory: an approach to prediction and interpretation in quantum chemistry. WIREs Comp. Mol. Sci. 2013, 3, 123−142. (20) von Niessen, W.; Schirmer, J.; Cederbaum, L. Computational methods for the one-particle green’s function. Comput. Phys. Rep. 1984, 1, 57−125.

(21) Richard, R. M.; Marshall, M. S.; Dolgounitcheva, O.; Ortiz, J. V.; Brédas, J.-L.; Marom, N.; Sherrill, C. D. Accurate Ionization Potentials and Electron Affinities of Acceptor Molecules I. Reference Data at the CCSD(T) Complete Basis Set Limit. J. Chem. Theory Comput. 2016, 12, 595−604. (22) Gallandi, L.; Marom, N.; Rinke, P.; Körzdörfer, T. Accurate Ionization Potentials and Electron Affinities of Acceptor Molecules II: Non-Empirically Tuned Long-Range Corrected Hybrid Functionals. J. Chem. Theory Comput. 2016, 12, 605−614. (23) Knight, J. W.; Wang, X.; Gallandi, L.; Dolgounitcheva, O.; Ren, X.; Ortiz, J. V.; Rinke, P.; Körzdörfer, T.; Marom, N. Accurate Ionization Potentials and Electron Affinities of Acceptor Molecules III: A Benchmark of GW Methods. J. Chem. Theory Comput. 2016, 12, 615−626. (24) Dolgounitcheva, O.; Diaz-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, 12, 627−637. (25) Blase, X.; Attaccalite, C.; Olevano, V. First-principles GW calculations for fullerenes, porphyrins, phtalocyanine, and other molecules of interest for organic photovoltaic applications. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 83, 115103. (26) Marom, N.; Caruso, F.; Ren, X.; Hofmann, O. T.; Körzdörfer, T.; Chelikowsky, J. R.; Rubio, A.; Scheffler, M.; Rinke, P. Benchmark of GW methods for azabenzenes. Phys. Rev. B: Condens. Matter Mater. Phys. 2012, 86, 245127. (27) Körzdörfer, T.; Marom, N. Strategy for finding a reliable starting point for G0W0 demonstrated for molecules. Phys. Rev. B: Condens. Matter Mater. Phys. 2012, 86, 041110. (28) Atalla, V.; Yoon, M.; Caruso, F.; Rinke, P.; Scheffler, M. Hybrid density functional theory meets quasiparticle calculations: A consistent electronic structure approach. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 88, 165122. (29) van Setten, M. J.; Caruso, F.; Sharifzadeh, S.; Ren, X.; Scheffler, M.; Liu, F.; Lischner, J.; Lin, L.; Deslippe, J. R.; Louie, S. G.; Yang, C.; Weigend, F.; Neaton, J. B.; Evers, F.; Rinke, P. GW100: Benchmarking G0W0 for Molecular Systems. J. Chem. Theory Comput. 2015, 11, 5665−5687. (30) Caruso, F.; Dauth, M.; van Setten, M. J.; Rinke, P. Benchmark of GW Approaches for the GW100 Test Set. J. Chem. Theory Comput. 2016, 12, 5076−5087. (31) Marom, N. Accurate description of the electronic structure of organic semiconductors by GW methods. J. Phys.: Condens. Matter 2017, 29, 103003. (32) Marom, N. Electronic structure of copper phthalocyanine from G(0)W(0) calculations. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 84, 195143. (33) Salomon, E.; Marom, N. Electronic structure of CoPc adsorbed on Ag(100): Evidence for molecule-substrate interaction mediated by Co 3d orbitals. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 87, 075407. (34) Körzdörfer, T.; Kümmel, S.; Marom, N.; Kronik, L. When to trust photoelectron spectra from Kohn-Sham eigenvalues: The case of organic semiconductors. Phys. Rev. B: Condens. Matter Mater. Phys. 2009, 79, 201205. (35) Körzdörfer, T.; Kümmel, S.; Marom, N.; Kronik, L. Erratum: When to trust photoelectron spectra from Kohn-Sham eigenvalues: The case of organic semiconductors [Phys. Rev. B 79, 201205 (2009)]. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 82, 129903. (36) Körzdörfer, T. On the relation between orbital-localization and self-interaction errors in the density functional theory treatment of organic semiconductors. J. Chem. Phys. 2011, 134, 094111. (37) Körzdörfer, T.; Kümmel, S. Single-particle and quasiparticle interpretation of Kohn-Sham and generalized Kohn-Sham eigenvalues for hybrid functionals. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 82, 155206. (38) Refaely-Abramson, S.; Sharifzadeh, S.; Govind, N.; Autschbach, J.; Neaton, J. B.; Baer, R.; Kronik, L. Quasiparticle Spectra from a 4738


Article


(60) Baerends, E. J. Perspective on ”Self-consistent equations including exchange and correlation effects” - Kohn W, Sham LJ (1965) Phys Rev A 140:133−1138. Theor. Chem. Acc. 2000, 103, 265− 269. (61) Baerends, E. J.; Gritsenko, O. V.; van Meer, R. The Kohn-Sham gap, the fundamental gap and the optical gap: the physical meaning of occupied and virtual Kohn-Sham orbital energies. Phys. Chem. Chem. Phys. 2013, 15, 16408−16425. (62) Perdew, J. P.; Parr, R. G.; Levy, M.; Balduz, J. L. DensityFunctional Theory for Fractional Particle Number: Derivative Discontinuities of the Energy. Phys. Rev. Lett. 1982, 49, 1691−1694. (63) Perdew, J. P.; Levy, M. Physical content of the exact Kohn-Sham orbital energies: band gaps and derivative discontinuities. Phys. Rev. Lett. 1983, 51, 1884. (64) Sham, L. J.; Schlüter, M. Density-functional theory of the energy gap. Phys. Rev. Lett. 1983, 51, 1888. (65) Perdew, J. P. What do the Kohn-Sham orbital energies mean? How do atoms dissociate? Density Functional Methods in Physics; 1985; DOI: 10.1007/978-1-4757-0818-9_10. (66) Chong, D. P.; Gritsenko, O. V.; Baerends, E. J. Interpretation of the Kohn-Sham orbital energies as approximate vertical ionization potentials. J. Chem. Phys. 2002, 116, 1760−1772. (67) Levy, M.; Perdew, J. P.; Sahni, V. Exact differential equation for the density and ionization energy of a many-particle system. Phys. Rev. A: At., Mol., Opt. Phys. 1984, 30, 2745. (68) Almbladh, C.-O.; von Barth, U. Exact results for the charge and spin densities, exchange-correlation potentials, and density-functional eigenvalues. Phys. Rev. B: Condens. Matter Mater. Phys. 1985, 31, 3231−3244. (69) Kraisler, E.; Kronik, L. Piecewise Linearity of Approximate Density Functionals Revisited: Implications for Frontier Orbital Energies. Phys. Rev. Lett. 2013, 110, 126403. (70) Kraisler, E.; Kronik, L. Fundamental gaps with approximate density functionals: The derivative discontinuity revealed from ensemble considerations. J. Chem. Phys. 2014, 140, 18A540. (71) van Setten, M. J.; Weigend, F.; Evers, F. The GW-Method for Quantum Chemistry Applications: Theory and Implementation. J. Chem. Theory Comput. 2013, 9, 232−246. (72) Becke, A. D. Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648−5652. (73) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab Initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields. J. Phys. Chem. 1994, 98, 11623−11627. (74) Weigend, F. Accurate Coulomb-fitting basis sets for H to Rn. Phys. Chem. Chem. Phys. 2006, 8, 1057−1065. (75) TURBOMOLE V6.2 2010, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989−2007, TURBOMOLE GmbH: since 2007. Available from http://www. turbomole.com. (76) Eichkorn, K.; Treutler, O.; Oehm, H.; Haeser, M.; Ahlrichs, R. Auxiliary Basis Sets to approximate Coulomb Potentials. Chem. Phys. Lett. 1995, 240, 283−290. (77) Weigend, F.; Ahlrichs, R. Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy. Phys. Chem. Chem. Phys. 2005, 7, 3297−3305. (78) Blum, V.; Gehrke, R.; Hanke, F.; Havu, P.; Havu, V.; Ren, X.; Reuter, K.; Scheffler, M. Ab initio molecular simulations with numeric atom-centered orbitals. Comput. Phys. Commun. 2009, 180, 2175− 2196. (79) Ren, X.; Rinke, P.; Blum, V.; Wieferink, J.; Tkatchenko, A.; Sanfilippo, A.; Reuter, K.; Scheffler, M. Resolution-of-identity approach to Hartree-Fock, hybrid density functionals, RPA, MP2 and GW with numeric atom-centered orbital basis functions. New J. Phys. 2012, 14, 053020. (80) Siegbahn, K.; Nordling, C.; Johansson, G.; Hedman, J.; Heden, P.; Hamrin, K.; Gelius, U.; Bergmank, T.; Werme, L.; Manne, R.; Baer, Y. ESCA applied to free molecules; 1969.

Nonempirical Optimally Tuned Range-Separated Hybrid Density Functional. Phys. Rev. Lett. 2012, 109, 226405. (39) Körzdörfer, T.; Parrish, R. M.; Marom, N.; Sears, J. S.; Sherrill, C. D.; Brédas, J.-L. Assessment of the performance of tuned rangeseparated hybrid density functionals in predicting accurate quasiparticle spectra. Phys. Rev. B: Condens. Matter Mater. Phys. 2012, 86, 205110. (40) Egger, D. A.; Weissman, S.; Refaely-Abramson, S.; Sharifzadeh, S.; Dauth, M.; Baer, R.; Kümmel, S.; Neaton, J. B.; Zojer, E.; Kronik, L. Outer-valence Electron Spectra of Prototypical Aromatic Heterocycles from an Optimally Tuned Range-Separated Hybrid Functional. J. Chem. Theory Comput. 2014, 10, 1934−1952. (41) Kohn, W.; Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. 1965, 140, A1133− A1138. (42) Seidl, A.; Görling, A.; Vogl, P.; Majewski, J. A.; Levy, M. Generalized Kohn-Sham schemes and the band gap problem. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 53, 3764. (43) Parr, R. G.; Yang, W. Density-Functional Theory of Atoms and Molecules; Oxford University Press: Oxford, 1989. (44) Dreizler, R. M.; Gross, E. K. U. Density Functional Theory; Springer: Heidelberg, 1990. (45) Koch, W.; Holthausen, M. C. A Chemist’s Guide to Density Functional Theory; Wiley-VCH: New York, 2000. (46) Sharp, R. T.; Horton, G. K. A Variational Approach to the Unipotential Many-Electron Problem. Phys. Rev. 1953, 90, 317. (47) Talman, J. D.; Shadwick, W. F. Optimized effective atomic central potential. Phys. Rev. A: At., Mol., Opt. Phys. 1976, 14, 36. (48) Görling, A. A new KS method for molecules based on an exchange charge density generating the exact local KS exchange potential. Phys. Rev. Lett. 1999, 83, 5459. (49) Ivanov, S.; Hirata, S.; Bartlett, R. J. Exact exchange treatment for molecules in finite-basis-set Kohn-Sham theory. Phys. Rev. Lett. 1999, 83, 5455. (50) Görling, A. Orbital-and state-dependent functionals in densityfunctional theory. J. Chem. Phys. 2005, 123, 062203. (51) Bleiziffer, P.; Heßelmann, A.; Görling, A. Efficient self-consistent treatment of electron correlation within the random phase approximation. J. Chem. Phys. 2013, 139, 084113. Note that in eqs 30, 31, and 48−51 of this reference there is a sign error in the denominators of the right-hand side. They should read 1 − σ(iω) instead of 1 + σ(iω). (52) Gö rling, A. Echange-correlation potentials with proper discontinuities for physically meaningful Kohn-Sham eigenvalues and band gaps. Phys. Rev. B: Condens. Matter Mater. Phys. 2015, 91, 245120. (53) Trushin, E.; Betzinger, M.; Blügel, S.; Görling, A. Band gaps, ionization potentials, and electron affinities of periodic electron systems via the adiabatic-connection fluctuation-dissipation theorem. Phys. Rev. B: Condens. Matter Mater. Phys. 2016, 94, 075123. (54) Perdew, J. P.; Yang, W.; Burke, K.; Yang, Z.; Gross, E. K. U.; Scheffler, M.; Scuseria, G. E.; Henderson, T. M.; Zhang, I. Y.; Ruzsinszky, A.; Peng, H.; Sun, J.; Trushin, E.; Gö rling, A. Understanding band gaps of solids in generalized Kohn-Sham theory. Proc. Natl. Acad. Sci. U. S. A. 2017, 114, 2801−2806. (55) Perdew, J. P. Density-functional theory and the band-gap problem. Int. J. Quantum Chem. 1985, 28, 497−523. (56) Niquet, Y. M.; Gonze, X. Band-gap energy in the random-phase approximation to density-functional theory. Phys. Rev. B: Condens. Matter Mater. Phys. 2004, 70, 245115. (57) Görling, A.; Levy, M. Hardness of molecules and the band gap of solids within the Kohn-Sham formalism: A perturbation-scaling approach. Phys. Rev. A: At., Mol., Opt. Phys. 1995, 52, 4493. (58) Baerends, E. J.; Gritsenko, O. V.; van Leeuwen, R. Effective oneelectron potential in the Kohn-Sham molecular orbital theory. Chemical Applications of Density-Functional Theory; 1996; DOI: 10.1021/bk-1996-0629.ch002. (59) Savin, A.; Umrigar, C. J.; Gonze, X. Relationship of Kohn-Sham eigenvalues to excitation energies. Chem. Phys. Lett. 1998, 288, 391− 395. 4739


Article

Journal of Chemical Theory and Computation (81) Gelius, U.; Allan, C.; Johansson, G.; Siegbahn, H.; Allison, D.; Siegbahn, K. The ESCA Spectra of Benzene and the Iso-electronic Series, Thiophene, Pyrrole and Furan. Phys. Scr. 1971, 3, 237−242. (82) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 1996, 77, 3865. (83) Perdew, J. P.; Ernzerhof, M.; Burke, K. Rationale for mixing exact exchange with density functional approximations. J. Chem. Phys. 1996, 105, 9982−9985. (84) Cederbaum, L.; Hohlneicher, G.; Niessen, W. Improved calculations of ionization potentials of closed-shell molecules. Mol. Phys. 1973, 26, 1405−1424. (85) Maier, J.; Muller, J.; Kubota, T.; Yamakawa, M. Ionisation Energies and the Electronic Structures of the N-oxides of Azanaphthalenes and azaanthracenes. Helv. Chim. Acta 1975, 58, 1641−1648. (86) Vydrov, O. A.; Scuseria, G. E. Assessment of a long-range corrected hybrid functional. J. Chem. Phys. 2006, 125, 234109. (87) Stein, T.; Kronik, L.; Baer, R. Prediction of charge-transfer excitations in coumarin-based dyes using a range-separated functional tuned from first principles. J. Chem. Phys. 2009, 131, 244119. (88) Richard, R.; Marshall, M.; Dolgounitcheva, O.; Ortiz, J.; Bredas, J.; Marom, N.; Sherill, C. Accurate Ionization Potentials and Electron Affinities of Acceptor Molecules I. Reference Data at the CCSD(T) Complete Basis Set Limit. J. Chem. Theory Comput. 2016, 12, 595− 604. (89) Dori, N.; Menon, M.; Kilian, L.; Sokolowski, M.; Kronik, L.; Umbach, E. Valence electronic structure of gas-phase 3,4,9,10-perylene tetracarboxylic acid dianhydride: Experiment and theory. Phys. Rev. B: Condens. Matter Mater. Phys. 2006, 73, 195208. (90) Marom, N. Electronic structure of copper phthalocyanine: A comparative density functional theory study. J. Chem. Phys. 2008, 128, 164107. (91) Marom, N. Density functional theory of transition metal phthalocyanines, I: electronic structure of NiPc and CoPc-selfinteraction effects. Appl. Phys. A: Mater. Sci. Process. 2009, 95, 159− 163. (92) Brundle, C.; Robin, M.; Kuebler, N. Perfluoro effect in photoelectron spectroscopy. II. Aromatic molecules. J. Am. Chem. Soc. 1972, 94, 1466−1475. (93) Dougherty, D.; McGlynn, S. Photoelectron spectroscopy of carbonyls. 1,4-Benzoquinones. J. Am. Chem. Soc. 1977, 99, 3234−3239. (94) Kimura, K. Handbook of HeI Photoelectron Spectra of Fundamental Organic Molecules: Ionization Energies, Ab Initio Assignments, and Valence Electronic Structure for 200 Molecules; Japan Scientific Societies Press; Halsted Press: 1981. (95) Rabalais, J. Photoelectron Spectroscopic Investigation of the Electronic Structure of Nitromethane and Nitrobenzene. J. Chem. Phys. 1972, 57, 960−967.

4740


Accurate Valence Ionization Energies from KohnâSham Eigenvalues

Recommend Documents

Accurate Valence Ionization Energies from KohnâSham Eigenvalues