Tuned Range-separated Density Functional Theory and Dyson Orbital

Sep 21, 2018 - User Resources. About Us · ACS Members · Librarians · ACS Publishing Center · Website Demos · Privacy Policy · Mobile Site ...
0 downloads 0 Views 9MB Size
Subscriber access provided by University of Sunderland

Spectroscopy and Excited States

Tuned Range-separated Density Functional Theory and Dyson Orbital Formalism for Photoelectron Spectra Tobias Möhle, Olga S. Bokareva, Gilbert Grell, Oliver Kühn, and Sergey I. Bokarev J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00707 • Publication Date (Web): 21 Sep 2018 Downloaded from http://pubs.acs.org on September 25, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Tuned Range-separated Density Functional Theory and Dyson Orbital Formalism for Photoelectron Spectra T. Möhle,† O.S. Bokareva,‡ G. Grell,† O. Kühn,† and S.I. Bokarev∗,† Institut für Physik, Universität Rostock, Albert-Einstein-Str. 23-24, 18059 Rostock, Germany, and Department of Physical Chemistry, Kazan Federal University, Kremlevskaya str. 18, 420008, Kazan, Russia E-mail: [email protected]

Abstract



Photoelectron spectroscopy represents a valuable tool to analyze structural and dynamical changes in molecular systems. Comprehensive interpretation of experimental data requires, however, involvement of reliable theoretical modeling. In this communication, we present a protocol based on the combination of well-established linear-response timedependent density functional theory and Dyson orbital formalism for the accurate prediction of both ionization energies and intensities. Essential here is the utilization of the optimallytuned range-separated hybrid density functionals, improving the ionization potentials not only of frontier, but also of the deeper lying orbitals. In general, the protocol provides accurate results as illustrated by comparison to experiments for several gas-phase molecules, belonging to different classes. Further, we analyze possible pitfalls of this approach, and namely, discuss the ambiguities in the choice of optimal range-separation parameters, the influence of the stability of the ground state, and the spin contamination issues as possible sources of inaccuracies. ———————————————————

1

Introduction

Photoelectron spectroscopy is one of the most widely used spectroscopic methods to address properties of matter in different states of aggregation. 1,2 However, experimental photoelectron spectra (PESs) of complex systems often exhibit a quite complicated structure and thus require reliable theoretical modeling for insightful interpretation. The computational approaches suggested for this purpose during the last decades both in frequency (energy) and time domains strongly differ in their complexity. Focusing on the frequency domain, in practice especially for extended systems, PESs are often predicted by calculating the single-particle molecular orbital (MO) energies, assuming that they correspond to the ionization energies in the spirit of Koopmans’ theorem. 3 To obtain these energies, Kohn-Sham DFT is frequently utilized due to its computational efficiency and accuracy. In case of DFT, only the eigenenergy of the highest occupied MO (HOMO) can be directly connected to the lowest ionization potential (IP), 4,5 whereas energies of the deeper lying MOs still retain their physical sense as approximate estimates of higher IPs. 6 In DFT, MOs represent quasiparticle wave functions, due to



To whom correspondence should be addressed University of Rostock ‡ Kazan University †

ACS Paragon Plus Environment

1

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

2

Journal of Chemical Theory and Computation

the implicit account for electron correlation effects via an exchange-correlation functional. However, by construction, this Koopmans’like scheme does not account for truly manyelectron effects such as combination or shake up transitions. Moreover, the PES intensities are chosen to be uniform for all transitions, 7–10 representing a quite rough approximation. Despite of its simplicity, this method provides a reasonable estimation of PESs for large systems, where more advanced methods become computationally unfeasible. Another quite popular frequency-domain scheme is the Green’s function based approach, representing a very powerful and more reliable alternative for the prediction of PES. 11,12 An important advantage over most other quantum chemical methods is that the transition energies are obtained directly as poles of the Green’s function rather than differences of the total energies of the initial and final states. Although formally only the initial non-ionized system is considered, many-electron ionization effects are included in this formalism. The PES intensities are commonly taken as residues of the respective poles and, thus, are independent on the kinetic energy of the photoelectron. This leads to the so-called sudden approximation (SA) and makes a natural connection to the Dyson orbital (DO) theory discussed below. 13 Neglecting this dependence might be a critical issue for low kinetic energies of outgoing electrons, possibly giving inaccurate results. To go beyond the SA-level and obtain PESs that depend on the photon energy, the theoretical approach needs to account for the photoelectron explicitly. Besides a large group of methods that treat the photoelectron and the bound electrons in a uniform way (see Refs. 14– 22 to name but a few), the Dyson orbital (DO) formalism 23 is a well-established scheme for the computation of PES. In this scheme, the wave function of the final state is constructed as a properly antisymmetrized product of the bound N − 1-electron counterpart Ψf + and the free electron function ψ el (Fig. 1). Thus, in this approach correlation effects between the outgoing electron and the ionic remainder are neglected, leading to an efficient protocol for the compu-

i

Page 2 of 20

f f+

Energy

1 INTRODUCTION

S

S- ½ S+½

e-

Figure 1: Schematic representation of the photoionization process: the initial N -electron state i with spin S is ionized by a photon with energy hν. The final state is represented by the photoelectron in a continuum state and a bound N − 1-electron state f + of either spin S + 1/2 or S − 1/2, respectively. tation of PESs. As a side remark, this scheme excludes the effects of interchannel coupling and further it cannot be applied to strong-field phenomena. The DOs themselves as well as the respective binding energies EB = Ef + −Ei are calculated from the bound wave functions and energies of the initial, i, and ionized, f +, systems. This implies that the calculation of the ground state of i and the ground and excited states of f + are required, see Fig. 1. The respective wave functions can in principle be obtained with any quantum-chemistry method. 24–27 The DO provides an intuitive interpretation of the transition between many-body states, being a wave function of the outgoing quasiparticle (electron) before ionization. If supplied with a proper representation of the photoelectron, the DO can be integrated providing transition matrix elements as a function of the electron kinetic energy. On the other hand, the need of computing every final state at the “excited state” level makes this method computationally more demanding compared to the Kohn-Sham orbital energy and even Green’s function approaches. DFT and its linear-response time-dependent (TDDFT) extension are the preferred theoreti-

ACS Paragon Plus Environment

2 COMPUTATIONAL PROTOCOL

Page 3 of 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

cal approaches for studying ground and excited states of medium-sized and large molecules. 28 The computational efficiency and the fact that no pre-knowledge on the system under study is required (such as setting active spaces) are the key benefits of the (TD)DFT method, making it very attractive for numerous applications. However, this does not liberate the user from the need of choosing an appropriate density functional, what could be a challenge on its own. 29 The popular GGA functionals include only an approximate description of the exchange-correlation potential, leading to a self-interaction error and violation of the derivative discontinuity condition. 5 These problems manifest themselves in wrong IPs, electronic polarisabilities, and energies of the chargetransfer electronic states. 30,31 This fact complicates PES calculations regardless whether they are evaluated at the DFT-based Koopmanns’, Green’s function, or TDDFT level. To overcome this drawback, range-separated hybrid functionals have been proposed. 32 They are beneficial for calculating PES since they accurately predict IPs 27,33 what motivates their choice as a main method for the current study. In this work, we address the performance of the PES computational protocol combining the DO-formalism with the TDDFT method employing optimally-tuned range-separated hybrid (OTRSH) density functionals. This procedure combines the accuracy of the IPs, provided via first principles functional tuning, with the explicit calculation of the photoionization dipole matrix elements, aiming at improved spectral intensities. The conclusions are exemplified for small model molecules representing different classes of compounds: water, benzene as a simple aromatic system, as well as S8 and Cu4– as examples of non-metallic and metallic clusters. The paper is organized as follows: First, the theoretical methodology for calculation of PES making use of DOs is briefly described in Section 2.1. The protocol and the respective computational details including a summary of the optimal tuning procedure are presented in Sections 2.2 and 2.3. We proceed with the results and discussion in Section 3 and finally give con-

3

clusions in Section 4.

2 2.1

Computational protocol Dyson orbital formalism

In this publication, PES cross sections are calculated in the framework of the DO formalism 26,34 analogously to previous implementations for linear-response TDDFT. 35,36 In the weak-field limit, the N -electron final state |Ψf i = A{|Ψf + i|ψ el (εk )i} can be represented as the properly antisymmetrized product of the N − 1-electron ionic state |Ψf + i with energy Ef + and the state of photoelectron |ψ el (εk )i having kinetic energy εk . Assuming the orthogonality between the wave function ψ el (εk ) and the molecular orbitals of the initial state, the expression can be rewritten in terms of |ψ el (εk )i and a quasi-one-electron (unnormalized) state |ΦDO if + i which is called Dyson orbital (see Eq. 1). Thus, the cross section of transitions between initial |Ψi i and final |Ψf i states is expressed as X σ(εk ) ∝ |hΨf |ˆ µ|Ψi iN |2 ρ(εk ) f

X 2 hψ el (εk )|ˆ ρ(εk ). ≈ µ|ΦDO i 1 if +

(1)

f

Here, we imply that the signal is integrated over the full solid angle, thus, eliminating the directional dependence. The DO represents an N − 1-electron integral over initial and final √ DO ionic states Φif + = N hΨf + |Ψi iN −1 as denoted by the subscript N − 1. Hence, one comprises the system and electron structure method information in the DO, being a reduced quantity. In other words, it corresponds to the photoelectron wave function before ionization, including orbital relaxation and many-electron effects. The N -electron integration in the first part of Eq. 1 is, thus, reduced to a one electron integral. The density of states ρ(εk ) depends on the binding energy, Ef + − Ei , as well as the lineshape of the PES bands and can be represented in simplest case as a δ-distribution ρ(εk ) = δ(Ef + + εk − Ei − ~ω), where ~ω is the energy of an incoming photon.

ACS Paragon Plus Environment

2 COMPUTATIONAL PROTOCOL

4

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Working in the framework of Kohn-Sham DFT, the initial (ground) state of the nonionized system can be taken as a single KohnSham determinant Θ0 , where the constituent orbitals are relaxed self-consistently. (Note the approximate nature of such a representation. 37 ) One performs another orbital optimization for the ground state of the ionized system, where an electron is removed from the HOMO with index k. Therefore, the DO for the groundground transition reads (omitting indices i and f + for simplicity) √ (2) ΦDO = N hΘk |Θ0 iN −1 . In turn, other transitions in Fig. 1 correspond to the excited electronic states of the ionized system with ground state determinant Θk . Within linear-response TDDFT in Tamm-Dankoff approximation 38 (TDA), these final ionized states f + have the form of a P configuration interactionlike expansion Ψf + = j,b Xjb Θbj . Here, Xjb are the elements of the eigenvector of the respective Casida equation AX = ω X. Further, these excited state determinants can be classified in terms of single hole (1h) Θkl and two holes one particle (2h1p) determinants Θakl . This leads to the expression for the DO: ( X √ ∗ Xlk hΘkl |Θ0 iN −1 ΦDO = N l

) +

X

∗ Xla hΘakl |Θ0 iN −1

,

(3)

l,a

where the index k again corresponds to the HOMO of the initial state and the summation indices l and a denote occupied and virtual orbitals of the ionized system, respectively. Note that the MOs of the initial and ionized states are not necessarily orthogonal and, thus, the overlaps hΘakl |Θ0 iN −1 are in general non-zero. The evaluation of Eq. 3 is done via an expan-

Page 4 of 20

sion of determinants in terms of MOs, ϕi , D E X 1 P(1) p ˜ hΘ|Θ0 i = ··· (−1) ϕ˜ ϕ0 P∈SN

D E P(N −1) P(N ) × ϕ˜N −1 ϕ0 · ϕ0 , (4) where quantities with tilde correspond to the final state and P is a permutation of N orbital indices. Further, MO overlaps can be obtained from MO-coefficient matrices C and overlap matrix of atomic basis functions SAO as ˜ T SAO C)ij . Thus, quantities X, C, hϕ˜i |ϕj i = (C ˜ C,and SAO are taken from the quantum chemical calculations and the DOs are computed with the approach described in Ref. 26. The obtained DO can be used directly to approximately estimate the cross section as being 2 proportional to ||ΦDO if + || , which is the so-called sudden approximation. 39 Its utilization implies the neglect of the kinetic energy dependence of the cross section. It should be noted that the respective intensities need to be weighted according to the multiplicity of the final |Ψf + i state to account for its spin degeneracy. Starting with a singlet determinant Θ0 this does not make any difference to the PES apart from a multiplication by a factor of 2 since all possible final states are doublets. However, it modifies the spectrum if an openshell initial state with total spin S 6= 0 is considered. In this case, transitions to both S +1/2 and S − 1/2 final state manifolds are possible as both a spin-up or spin-down electron can be removed from the system. In case of non-singlet states, unrestricted TDDFT needs to be performed and thus the resulting wave function is in general spincontaminated. 40 In this work, we follow the suggestion of Wittbrodt et. al. 41 and use the spin contaminated wave functions and energies; see also discussion in Section 3.4. The respective transition amplitudes are weighted by the multiplicity as if the wave functions had pure spin.

ACS Paragon Plus Environment

2 COMPUTATIONAL PROTOCOL

Page 5 of 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

2.2

5

Journal of Chemical Theory and Computation

OTRSH

Although simply adding a constant HartreeFock-like orbital-dependent exact exchange mitigates the self-interaction error, often it leads to SCF stability issues. 42–44 A more gentle choice is to introduce exact exchange predominantly at large interelectron distances. 32,45 The generalized form of the short-/long-range partitioning of the Coulomb interaction proposed by Yanai et al. 46 reads as: 1 − [α + β · Γ(ωr12 )] α + β · Γ(ωr12 ) 1 = + , r12 r12 r12 (5) where Γ(ωr12 ) is a smooth switching function, varying from 0 to 1. It damps the exchange contribution from the parent density functional (first summand in Eq. 5) and complements it with exact Hartree-Fock exchange (second summand) which has the correct asymptotic behavior. Here, ω sets the switching rate and its inverse corresponds to a characteristic interelectron distance, at which one enters the asymptotic regime. The dimensionless parameter α scales a global r12 -independent exact-exchange contribution (α ∈ [0, 1]). In this work, we choose β = 1 − α to ensure that exact exchange cancels the self-interaction exactly. Certain standard values for the rangeseparation parameter ω have been established using training sets consisting of atoms or small molecules. 47,48 However, the extension of the electron density sets the characteristic interelectron distance. Thus, when dealing with notably different molecules, a more reasonable approach is to choose a system-dependent ω. 30 Tuned functionals (both for 1D and 2D optimization considering only ω and (α, ω)-pair, respectively) have been successfully applied in a number of works, yielding improved descriptions of various molecular properties related to fundamental and optical gaps such as IPs, 49,50 charge-transfer and Rydberg transition energies, 51–53 optical rotation, hyperfine couplings, etc, see e.g. Refs. 7,33,54–58. The asymptotic decay of the electron density is connected to the IP as limr→∞ ρ(r) ∝ √ exp (−2 2 IPr). That is why it is logical to select the parameters in Eq. 5 to repro-

duce the IP. To perform a tuning of the functional for a particular system, the so-called ∆SCF method 51,59,60 has been applied. In this method, the IP is calculated as the difference between ground state (gs) energies of the systems with N and N − 1 electrons, i.e. α,ω α,ω IPα,ω (N ) = Egs (N − 1) − Egs (N ).

(6)

This yields the tuning condition α,ω min J(α, ω) = min |εα,ω (N )| . HOMO (N ) + IP α,ω

α,ω

(7) This procedure provides a manifold of (α, ω)pairs where J is minimal (Fig. 2a)) and one requires a further criterion to select optimal values. For the exact exchange-correlation density functional, the energy E(N ) must vary linearly for fractional electron numbers between integer N s. 5 This, however, does not hold true for many functionals. Accordingly, segments of E(N ) have a certain curvature (could be both positive and negative) referred to as (de)localization error, 61 for a detailed discussion see Ref. 58. Tuned range-separated hybrid functionals result usually in small E(N ) curvatures, ensuring a small delocalization error. However, to assist an unambiguous choice of optimal parameters, we have used the curvatures of the E(N )-dependence for fractional charges and have chosen, whenever possible, the (α, ω)-pair whose curvature is closest to being piecewise linear. As a particular measure of the curvature, we used Z N Z N ∆= dn∆(n) = dn |E(n) − E(N ) N −1

N −1

− [E(N ) − E(N − 1)](N − n) | m

(8)

with m = 1, 2 to indicate the deviation from the requested linear dependence, see also Fig. 2b). Both values of m yielded very similar measures of linearity.

2.3

Computational Details

In the present paper, the tuning of the rangeseparation parameters α and ω has been done for the generalized form (Eq. 5) of the LC-

ACS Paragon Plus Environment

2 COMPUTATIONAL PROTOCOL

6

Journal of Chemical Theory and Computation

0.8 0.7

2.8

a)

J (α, ω)

2.4

[meV]

0.6

2.0

0.5 α

BLYP functional. 62 An error function kernel Γ(ωr12 ) = erf(ωr12 ) has been applied for the exchange part complementing regular DFT. The (α, ω)-optimization has been performed using the LANL2DZ ECP basis set for Cu and the 6-31G(d) basis set for all other atoms. The step size for varying α has been 0.05 and for ω 0.01 bohr−1 . The geometries are optimized using the LCBLYP functional with standard parameters. Since the constraint to full symmetry of the respective point-groups was found to lead to negative stability matrix eigenvalues 63 for all parameters, the optimization is done without symmetry-constraints. The gas phase ground state DFT and excited state TDA calculations have been conducted employing the optimized α and ω and the def2TZVP 64 basis set. To describe the extended anionic wave function of Cu4– , the def2-QZVPP 64 basis has been employed. To cover the relevant binding energy ranges, 20 excited doublet states were taken into account for the H2 O+ molecule, 150 for C6 H6+ and 280 states for S8+ . For Cu4 , both singlet and triplet states of the neutral system were taken into account, with 70 states in each manifold. The intensities of the particular ionization channels have been weighted with the respective spin-degeneracy factors of the final states. All calculations have been performed with the Gaussian 09 program package; 65 the calculations with partial charges have been done with the NWChem program suite. 66 The atomic overlap matrix, Kohn-Sham MO and CI-coefficients have been used for the calculation of DOs as described in Ref. 26. The photoelectron has been represented by a plane wave expanded in Coulomb waves that describe the unbound electron, assuming a spherical Coulomb-potential. 67 An exclusion is the negatively charged Cu4– cluster, where the final state of the system is uncharged and, thus, spherical waves should better describe the outgoing photoelectron wave function ψ el (εk ). Note that both representations in terms of Coulomb and spherical waves are approximations; this especially holds for low-symmetry molecules. The maximum angular momentum

1.6

0.4

1.2

0.3 0.2

0.8

0.1

0.4

0.0

0.0 0.1

2.0 b)

0.2 0.3 −1 ω[bohr ] (0.10, 0.29) (0.45, 0.19)

0.4

0.5

(0.30, 0.24) (0.25, 0.25)

1.5 ∆(n) [meV]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 20

1.0 0.5 0.0 −0.5 −1.0 −1.5 1.0

0.8

0.6 0.4 N −n

0.2

0.0

Figure 2: a) J(α, ω), Eq. 7, for the example of the S8 molecule. Red points denote minima of J(α, ω) at constant α values. b) The deviation from the linear dependence, ∆(n), see Eq. 8, for the S8 molecule. in the expansion of ψ el (εk ) is chosen to be l = 10. The integration of the DO with the ψ el (εk ) is performed using the ezDyson 4.0 program. 68 It uses a regular (rectangular) grid which has been constructed in each case by adding 4 Å to every side of the molecule with a density of 30 points per Ångström. Due to the diffuse character of the Cu4– DOs, the box size is chosen 9Å larger than the size of the cluster in this case. In all cases, the results have been averaged over all possible polarization orientations of the incoming light, accounting for free tumbling of the molecules in the gas phase and/or the unpolarized character of the ionizing radiation. The photoelectron signal has also been integrated with respect to the angular dependence.

ACS Paragon Plus Environment

3 RESULTS AND DISCUSSION

Page 7 of 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

3

Results and discussion

The accuracy of the theoretical PES obtained with the TDDFT DO method depends on several points addressed in this section. First, we discuss the design of the optimally-tuned functional in Section 3.1 and the closely connected problem of the ground state stability of the ionic remainder in Section 3.2. The actual PESs employing different approximations are presented in Section 3.3. Finally, we comment on the spin purity of the excited states in Section 3.4.

3.1

Optimal tuning

In order to optimize the range-separation parameters for each compound, single-point calculations on the (α, ω) grid have been done. A typical example of a 2D plot of J(α, ω) is shown in Fig. 2a) for the S8 molecule. The values of optimal range-separation parameters, fulfilling Eq. 7 and the lowest curvature condition, are given in Table 1 for all considered systems. Table 1: Optimized range separation parameters, α and ω, for different molecules. Compound H2 O C 6 H6 S8 Cu4– Cu4– a) Cu4– b)

α 0.00 0.25 0.10 0.00 0.20 0.50

ω [bohr−1 ] 0.59 0.26 0.29 0.25 0.22 0.16

a), b)

Alternative sets of parameters corresponding to minimal piecewise curvature. Note that in case of the closed-shell systems, benzene, water, and S8 , there is no uncertainty with respect to the choice of the multiplicity of the N − 1-electron system. However, when N is odd, corresponding to a doublet initial state, singlet and triplet multiplicities of the ionized species are possible. For Cu4– , only the singlet case leads to a distinct minimum valley for J (see Supplement), whereas for the triplet case such a valley is absent. Therefore, the set of parameters for this “singlet” case is considered further.

7

Importantly, the criterion of smallest piecewise curvature does not in general provide a tool to unambiguously select the best (α, ω)-pair as this curvature is very small for OTRSH functionals. An example is given in Fig. 2b) for the S8 system, with curvatures being three orders of magnitude smaller than for standard density functionals. 69 Further, this criterion yields ambiguous values for Cu4– , having three equivalent points for which the curvature is very low (Table 1). Even the additional requirement of the variational stability of the ionic ground state (see Section 3.2) does not lead to a solution. To resolve this ambiguity, we have chosen the pair for which the resulting IP is closest to the experimental value.

3.2

Ionic ground state stability

Along with the elimination of self-interaction and delocalization errors, one needs to pay attention on the stability of the ground-state DFT solution, which is a crucial issue for TDDFT to predict reliable excitation energies. 43 This is especially important for calculation of the excited states of the ion, when the molecule has a high symmetry, resulting in degenerate MOs and a singlet closed-shell unionized ground state. Removing an electron from one of such MOs leads to an ambiguity, in which of the degenerate MOs the hole is residing. Since DFT with approximate functionals is a single-reference method, it cannot tackle such situations which require in general a multi-reference treatment, what is reflected in the slight artificial lowering of the energy of one of the orbitals from a degenerate pair. Thus, performing the SCF procedure for the ionic ground state taking orbital relaxation into account, results in the lowest excited states lying very close to the ground one. This situation can lead to, e.g., triplet instabilities and imaginary excitation energies for the lowest electronic states of the ionized system. The same problem may arise even for low-symmetry molecules if the frontier occupied MOs are very close in energy. Note that usually the ground states of the unionized systems are stable with this respect. Even though TDA calculations are not af-

ACS Paragon Plus Environment

3 RESULTS AND DISCUSSION

8

Journal of Chemical Theory and Computation

0.7

S(α, ω)

b)

S(α, ω)

0.00 -0.16

0.22

0.6

0.5

0.17

0.5

-0.32

0.4

0.13

0.4

-0.48

0.3

0.09

0.3

-0.65

0.2

0.05

0.2

-0.81

0.1

0.01

0.1

-0.97

-0.03

0.0

[eV]

[eV]

α

0.6

0.7

0.26

a)

α

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 20

0.0

0.1

0.2 0.3 ω [bohr−1]

0.4

0.5

0.1

0.2 0.3 ω [bohr−1]

0.4

0.5

-1.14

Figure 3: Lowest eigenvalues S(α, ω) of the stability matrix (orbital Hessian) for the ground states of a) C6 H6+ and b) S8+ . Note that imaginary values are denoted as negative in the plots. The minimum values of J are depicted by red circles. For clarity, the zero value of S is denoted by dashed lines. fected by such instabilities, 43,44 parameters giving negative eigenvalues of the stability matrix still indicate unphysical situations, making the stability analysis an important validity check in the present protocol. In Fig. 3, the lowest eigenvalue of the stability matrix is shown as a function of α and ω for the C6 H6+ and S8+ ions. The zero eigenvalues are denoted by a dashed line. The stable regions do not coincide with the valley of J(α, ω) (depicted with red points). Fortunately, for both cases values of optimal (α, ω)-pairs (Table 1) lie within the positive stability regions, although for S8+ it is quite close to zero (about 8·10−5 eV). For other test systems, H2 O+ and Cu4 , the lowest stability matrix eigenvalues are positive in the full range of α and ω.

3.3

Photoelectron spectra

The accuracy of PES calculations depends on two factors: i) description of the photoionization itself and ii) reliability of the underlying electronic structure method. In this section, we discuss both issues for test systems staying within the TDA DFT technique. The highest level in the hierarchy of photoionization description employed in this work is the explicit integration of the dipole matrix element between the photoelectron, in the form of

Coulomb-waves expansion, and the DO (Eq. 1). This result is denoted as ’Full’ calculation below. Further moving down the ladder, the dependence on the photoelectron wave function ψel (εk ) is neglected resulting in the sudden approximation, denoted as SA. It is justified in the limit of high kinetic energies. 39 Finally, the simplest method is based on the Koopmans’-like approach where the peak positions are set at the binding energies of the MOs of the initial state, assuming uniform intensities. This approach, called MO density of states (MO-DOS), corresponds to the SA in the frozen orbital limit. Thus, no orbital relaxation and combination transitions are taken into account. Nevertheless, when speaking about Kohn-Sham MO energies, one should view them as quasi-particle energies that contain electron correlation in the initial state implicitly. To study the effect of range-separation, we also compare the OTRSH results to those obtained with conventional BLYP and B3LYP functionals. The former has been chosen since the OTRSH scheme used in this work is based on the non-hybrid BLYP and the effect of range-separation can be clearly extracted from the comparison. The B3LYP functional has a constant HF-exchange contribution of 20% and is a popular and robust functional. The stick spectra from all approaches are

ACS Paragon Plus Environment

3 RESULTS AND DISCUSSION

Page 9 of 20

1.2 a) Intensity [arb. u.]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

9

Journal of Chemical Theory and Computation

b)

0.9 0.6 0.3 0.0

1.5

2.0

2.5

3.0

3.5

c)

4.0

1.5

2.0

2.5

3.0

3.5

4.0

d) 1.5

2.0 2.5 3.0 Binding energy [eV]

3.5

a) Exp. Full SA MO-DOS

b) Exp. OTRSH B3LYP BLYP

c)-e) alpha beta

-1.0

0.0

1.0

2.0

3.0

4.0

1.0 2.0 3.0 Binding energy [eV]

4.0

e) -1.0

0.0

Figure 4: PES of Cu4– computed with different schemes as compared with experimental data (black line). 70 Legends for the respective panels are given at the bottom left. The spectra in panel a) are based on OTRSH (TD)DFT-calculations for different levels of theory to compute transition energies and intensities. In panel b), PESs computed with the integration of Dyson orbitals (Full) are shown; they differ in the exchange-correlation functional being used. Panels c), d) and e) show the energetic positions of the initial state MO energies (thus corresponding to the stick-spectra of MO-DOS) for OTRSH, B3LYP, and BLYP, respectively. Note the different energy ranges of the panels. broadened using a pseudo-Voigt profile with a uniform width for all transitions to resemble the respective experimental data (water is the only exception, see Supplement). Moreover, each PES is normalized individually according to the spectral feature with highest intensity. Copper cluster We start the discussion with the comparison of approaches for the Cu4– system as presented in Fig. 4. The choice of this cluster is justified since the observed effects follow the same trends as for other systems and it contains few well-resolved lines, making the comparison easier. On the other hand, Cu4– is peculiar as it has a doublet ground state. The ionization occurs via two spin channels and leads to the appearance of additional bands. Fi-

nally, in contrast to other systems, the final singlet and triplet states are almost free of spin contamination (see Section 3.4 for a detailed discussion), leading to more reliable results. The PES comprises four main one-electron transitions with binding energies below 3.7 eV. The first three DOs represent p-type combinations of atomic orbitals, cf. insets in Fig. 4a). Here, the second and third peaks differ mainly in the spin of the final state, whose relative weights of 3 and 1 explain the intensity ratio of the computed spectra. The fourth peak corresponds to a combination of atomic d-orbitals. There is also an intense feature at around 4.0 eV which stems from a number of combination transitions. However, the low photon energy of 4.66 eV employed in the experiment 70 leads to a

ACS Paragon Plus Environment

3 RESULTS AND DISCUSSION

10

Journal of Chemical Theory and Computation

quite low kinetic energy of the photoelectrons; thus, presumably this region can not be treated reliably neither experimentally nor within our theoretical model. This may also explain the poor agreement of the intensities of the peak at 3.5 eV. In general, the shape of the spectrum for lower binding energies is quite well reproduced with DO-based theoretical methods, see Fig. 4a). To note is the excellent agreement of the positions of the first two peaks with experiments. 70–73 Remarkably, the peak at 2.7 eV corresponding to the doublet→singlet transition is absent in the MO-DOS spectrum. Instead, it shows up as a shoulder of the peak at 2.3 eV. Note that here we have summed up α and β MO-DOS with respective weights of 1 and 3 but the energy distance between respective spin-orbitals is notably smaller than that observed in the experiment. Thus, for non-singlet initial states an explicit excited state calculation could be necessary for qualitative reproduction of the PES since MO-DOS approach could be not accurate enough. However, the ratio of intensities for the first two peaks at 1.4 and 2.3 eV is well-reproduced for all approaches including MO-DOS. Both SA and Full calculations underestimate the height of the third peak as the experimental intensity ratio for the second and third peaks is deviating from 3:1 which almost strictly holds for the theoretically predicted spectrum. The SA leads to results which are fairly close to the Full calculation for the bands below 3.0 eV. Even though the SA is expected to break down due to the low kinetic energy, the Full calculation should still be valid in this region. It is thus surprising that the peak at 3.5 eV is fairly well reproduced by the SA but is almost absent in the full spectrum shown in panel a). The reason for this underestimation most probably lies in the representation of the photoelectron as a spherical wave, emitted from the center of the respective non-spherically symmetric DO. The PESs on the MO-DOS level are shown only as stick spectra in the panels c) OTRSH, d) B3LYP and e) BLYP of Fig. 4. It is apparent from these spectra that the choice of the functional is very critical here. While the stan-

Experiment (0.00,0.25) (0.10,0.23)

(0.20,0.22) (0.35,0.19) (0.50,0.16)

Intensity [arb. u.]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 20

1.0

1.5

2.0 2.5 3.0 3.5 Binding energy [eV]

4.0

Figure 5: SA PES of Cu4– computed with different optimal (α, ω)-pairs. dard functionals predict an autoionizing character for the ground state and possess also a qualitatively wrong energetic structure for the lower lying orbitals, the energies predicted by the OTRSH functional fit the most important transitions (see also the light blue (dashed) line in panel a)). Hence, the energetic structure of all electronic states is strongly improved within the OTRSH-formalism as it has also been reported before, see, e.g., Ref. 27. This effect can be observed in most cases and shows that this scheme is not only a fit of the first IP as it might look at the first glance; the quasi-particle energies of deeper lying orbital are also improved in this approach. Despite such dramatic differences between the functionals in the MO-DOS, it is astonishing that the explicit TDA calculation and the DO formalism significantly improve PESs as shown in Fig. 4, panel b). Interestingly, the explicitly calculated IPs are in very good agreement with the experiment for all three functionals, although the intensity of the first peak at 1.4 eV is notably underestimated when using BLYP and B3LYP functionals. The general shape of the PES computed with B3LYP has slightly poorer quality for lower binding energy peaks if compared to that obtained with OTRSH, although

ACS Paragon Plus Environment

3 RESULTS AND DISCUSSION

Water The PESs of gaseous water obtained with the OTRSH functional and shown in panel 6a) are in very good agreement with the experiment 74 for all levels of theory. Note that, in contrast to other systems, the spectra in panel a) are normalized with respect to the area under the first peak since the experimental data show vibrational progressions that are not taken into account in the theoretical procedure. Moreover, since the linewidths of the

1.0

a)

Exp. MO-DOS SA Full

0.9 Intensity [arb. u.]

it outperformes the latter for features above 3 eV. BLYP reproduces the positions of the first two peaks that essentially resemble the singlettriplet energy gap. However, starting from the third peak and especially at higher binding energies, it is in very bad agreement with the experiment. There is no truly unambiguous criterion which (α, ω)-pair to chose in general case and especially for Cu4– as discussed in Section 3.1. Even the additional criterion of stability of the ionic solution does not apply here as it is satisfied for all considered values of range-separated parameters. Therefore, we have calculated PESs for this system using several points in the J(α, ω) ≈ 0 valley (Fig. 5). One can see that different values of parameters lead to an overall shift of the first three features in the PES. The lines at higher binding energies are shifted in the opposite direction and the structure of the respective transitions also substantially changes as is apparent from the corresponding stick spectra. As mentioned in Section 3.1, for our main analysis we have chosen α and ω out of three possibilities (Table 1) such that the computed IP coincides with the experiment. In cases, when such a selection cannot be justified on the basis of experiments, other criteria are required and finding them is the task for future research. The observations made for Cu4– are also supported by the results obtained for H2 O, C6 H6 and S8 which have a singlet initial state. The respective spectra obtained with the OTRSHformalism on different levels of theory are shown in Fig. 6 and compared to experimental data. The most important DOs are shown at the respective peaks.

0.8 0.3 0.2 0.1 0.0

8

1.4 Intensity [arb. u.]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

11

Journal of Chemical Theory and Computation

12

16

20

32

36

b)

1.2 1.0 0.8 0.6 0.4 0.2 0.0

8

10

12

14

16

18

20

1.2 c) Intensity [arb. u.]

Page 11 of 20

1.0 0.8 0.6 0.4 0.2 0.0

8

9

10 11 12 13 Binding energy [eV]

14

15

Figure 6: PESs of a) water with a photon energy of 180 eV 74 (for broadening, see Supplement), b) benzene with a photon energy of 21.22 eV 75 and c) Sulfur vapor at 100◦ C (S8 ) with a photon energy of 21.22 eV 76 as computed with the MO-DOS, SA and Full schemes and compared to experiment.

ACS Paragon Plus Environment

3 RESULTS AND DISCUSSION

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

transitions differ strongly in the PES of water, we used different broadening parameters for different peaks, see Supplement. The IPs of the first three peaks are in excellent agreement with the experiment, albeit MO-DOS gives positions of the first and fourth bands which are slightly deviating from the TDA ionization energies. The position of the feature at 33 eV is overestimated by all approaches by about 1 eV. This is not surprising since highly excited states are predicted less reliably with TDA. It should be noted that several transitions contribute to this region but with very low intensity. The energetic positions of the MO-DOS spectral features for the B3LYP and BLYP functionals are sufficiently different from those of the OTRSH one (see Supplement). Due to the high photon energy of 180 eV, SA could be expected to hold, meaning that the results of SA and Full calculation should almost coincide. This is, however, true only for the third peak at around 19 eV. (The height of the first peak for all theoretical approaches is normalized to one; thus intensities coincide by definition.) For the second and fourth peaks at 15 eV and 33 eV, full integration gives notably lower intensities. This fact might be connected to the approximate representation of the photoelectron. Benzene The PES of benzene, shown in panel 6b), is an example of a system where orbital relaxation upon ionization plays an important role, leading to clear differences between the results obtained with MO-DOS (blue-green, dashed) and SA (wine, dash-dotted). The DObased spectra (SA and Full) also show some clear differences between each other. The Full calculation leads to a notably better agreement with the experiment for the intensity of the first peak. However, it also strongly overestimates the heights of the features at binding energies above 16 eV. These differences between the SA and Full approaches can be explained by the substantially different nature of the DOs, e.g., they have different effective angular momentum as can be seen from the insets in panel 6b). This fact leads to qualitative differences intro-

12

Page 12 of 20

duced by the explicit integration if compared to SA. In general, among all levels of theory, the SA seems to reproduce the shape of the spectrum best. This is especially apparent for the peak pair at 14-16 eV, where the Full calculation yields a reverse order of intensities. Regarding the energetic structure, the first IP is well reproduced, but the double degeneracy of the respective states is lifted due to orbital relaxation in the ionized state, leading to an overall broader lineshape. Apart from the first peak, the computed energies of other features are overestimated by up to 1 eV. In addition, the 14-16 eV peak pair is broader than in the experiment. Further, the MO-DOS of BLYP and B3LYP have a substantial offset of 3.4 and 2.4 eV, respectively, but the overall structure is predicted fairly well, see Supplement. Octasulfur The spectra of S8 demonstrate only minor differences between the methods (Fig. 6c)). The good agreement of the MO-DOS spectrum with the PES computed with the DO formalism shows that for S8 electron relaxation is not important. Moreover, the shape of the DOs is very similar across the spectrum as they all are combinations of sulfur 3p orbitals, explaining the minor differences between the Full integration and SA results. Despite the optimal tuning, the first IP is overestimated if compared to experiment. However, its experimental value is also under debate. 76–81 Similar to the PES of benzene, the computed spectrum is energetically stretched towards higher binding energies. Further, the height of the feature at 11.5 eV (in experiment) is notably underestimated by all computational schemes, whereas the intensities of the groups of bands at 9-10.2 and 12.5 eV are predicted fairly well. While comparing to experiment one should keep in mind that it has been conducted on sulfur vapor at 100 ◦ C. 76 Thus, it is expected that the vapor is in thermal equilibrium and contains only 67.3% of S8 . 82 Some of the differences between experimental and theoretical data thus can be attributed to this fact. Note that another experimental PES taken under slightly different conditions is available 77 whose inten-

ACS Paragon Plus Environment

4 CONCLUSIONS

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Intensity [arb. u.]

Page 13 of 20

Journal of Chemical Theory and Computation

2.4

1.0 Exp. Full 0.8

2.0 1.6

0.6

1.2

0.4

0.8

0.2

0.4

0.0

9

12 15 18 Binding energy [eV]

21

0.0

Figure 7: The PES of benzene computed with the Full scheme using the OTRSH functional. The color of the sticks represents the spin contamination hS 2 i − 0.75 of the respective (doublet) final states. sities are in worse agreement with our calculations. Thus, octasulfur appears to be the most problematic system out of those chosen for the current work, accounting for the possible ambiguity in experiment.

3.4

Spin purity and other issues

Before concluding, we would like to point to some problems that might be encountered within the presented protocol. First of all, linear-response TDDFT or TDA are intrinsically lacking double-electron excitations. This fact can lead to a poor agreement in the region of higher binding energies, where double excitations start to play role. In most cases, this limitation should have only a minor influence, since double-excitations naturally have only very small dipole intensities in photoelectron spectroscopy. However, it might become critical if ionization from an electronically excited state takes place, e.g., as is relevant for optical pump/photoelectron probe transient spectroscopy. 83–85 In this case, 2h1p combination transitions cannot be predicted without taking double excitations into account. This problem could even occur for main 1h transitions when the occupied and virtual orbitals change their order upon ionization, being a case of extremely

13

strong orbital relaxation. Thus, such situations are by construction out of reach for the present DO-based protocol. Another reason for poor agreement in the energetic region of inner-valence ionized states may be spin contamination. In a sense, it is unavoidable, since either non-ionized or ionized or even both (e.g. for doublet→triplet channel) state manifolds need to be treated with unrestricted TDDFT. Apart from errors in energy, spin contamination has a direct influence onto the intensities which are weighted with the spin degeneracies. The quite severe case of the benzene molecule is illustrated in Fig. 7. In this figure, the color of the sticks denotes the deviation of the expectation value hS 2 i for the doublet final manifold from its exact value of 0.75. The region of largest differences from the experiment coincides with the region of largest spin contamination which among other reasons could explain the poor agreement. However, to the best of our knowledge, a consensus on the proper treatment of spin contamination in the TDDFT case has not yet been reached. 40,41,86–88 Thus, in this work we ignored the contamination and treated all states as being spin-pure. However, appropriate treatment of this issue in the context of PES requires further investigation.

4

Conclusions

The developed computational protocol combines the Dyson orbital formalism to calculate PESs with the Tamm-Dankoff approximation employing the optimally-tuned range-separated hybrid density functional. Such a combination should in principle lead to semi-quantitative prediciton of the PESs of medium-sized systems relying on: i) computational efficiency of the DFT-based techniques, ii) accuracy of ionization potentials ensured via a first principles tuning procedure, and iii) explicit calculation of the bound-to-continuum transition matrix elements for intensities. For gas-phase molecules from different classes, it leads to results superior to the quite popular approximation of PES from the density of the molecular orbitals.

ACS Paragon Plus Environment

4 CONCLUSIONS

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

However, it requires larger computational effort, thus, limiting the size of the studied systems. Further, this approach provides a more consistent route to predict spectra of systems with non-singlet ground states, where different ionization spin-channels are possible. In this case, the respective state manifolds with different multiplicities can be explicitly calculated as exemplified for the Cu4– cluster. More specifically, concerning the density functional, the OTRSH-scheme leads to much better energetic positions in the MO-DOS scheme which has been also reported before. For the DO formalism, this scheme also leads to an improvement of the agreement with experiments. However, due to explicit calculation of the electronic states accounting for orbital relaxation and differential correlation upon ionization, the dependence on the functional is sufficiently reduced. Even standard functionals like BLYP and B3LYP might lead to reasonable PESs. An important note is that finding the optimal range-separation parameters cannot be unambiguously done in the general case. Even the additional criteria ensuring the piecewise linearity of the energy vs. fractional orbital occupations and ionic ground state stability cannot resolve this unambiguity. Thus, the formulation of additional criteria is required for the OTRSH functional. To mention is also a strong influence of this ambiguity on the shape of PESs which makes the search for such conditions warrant. Finally, the DO formalism or a similar theoretical method needs to be applied in general to achieve appropriate quality of the spectra and to prevent misinterpretation due to missing peaks and too simplistic models for the intensities. However, our study suggests that the Full DO scheme, including the costly integration to obtain the dipole moment matrix elements, can be omitted employing the more robust sudden approximation, unless a more accurate model for photoelectron is applied. Acknowledgement The authors would like to acknowledge financial support from the following sources: Landesgraduiertenförderung of Mecklenburg-Vorpommern (T.M.), Grant No. 14.Y26.31.0019 from Ministry of Education and

14

Page 14 of 20

Science of Russian Federation (O.S.B.), and Deutsche Forschungsgemeinschaft grant No. BO 4915/1-1 (S.I.B.).

Supporting Information Supporting information file contains description of the broadening parameters of the spectral lines, J(α, ω) for the Cu4– system assuming singlet and triplet ionized states, comparison of MO-DOS and Full spectra for different density functionals for water, benzene, and octasulfur systems, comparison of octasulfur spectrum to different experimental data.

References (1) Stefan Hüfner, Photoelectron troscopy; Springer: Berlin, 1996.

Spec-

(2) Winter, B.; Faubel, M. Photoemission from Liquid Aqueous Solutions. Chem. Rev. 2006, 106, 1176–1211. (3) Koopmans, T. Über die Zuordnung von Wellenfunktionen und Eigenwerten zu den einzelnen Elektronen eines Atoms. Physica 1934, 1, 104–113. (4) Almbladh, C.-O.; von Barth, U. Exact Results for the Charge and Spin Densities, Exchange-Correlation Potentials, and Density-Functional Eigenvalues. Phys. Rev. B 1985, 31, 3231–3244. (5) Perdew, J. P.; Parr, R. G.; Levy, M.; Balduz, J. L. Density-Functional Theory for Fractional Particle Number: Derivative Discontinuities of the Energy. Phys. Rev. Lett. 1982, 49, 1691–1694. (6) 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. (7) Egger, D. A.; Weissman, S.; RefaelyAbramson, S.; Sharifzadeh, S.; Dauth, M.; Baer, R.; Kümmel, S.; Neaton, J. B.; Zojer, E.; Kronik, L. Outer-Valence Electron

ACS Paragon Plus Environment

4 CONCLUSIONS

Page 15 of 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Spectra of Prototypical Aromatic Heterocycles from an Optimally Tuned RangeSeparated Hybrid Functional. J. Chem. Theor. Comp. 2014, 10, 1934–1952. (8) Yepes, D.; Seidel, R.; Winter, B.; Blumberger, J.; Jaque, P. Photoemission Spectra and Density Functional Theory Calculations of 3d Transition Metal–Aqua Complexes (Ti–Cu) in Aqueous Solution. J. Phys. Chem. B 2014, 118, 6850–6863.

15

(16) Colle, R.; Embriaco, D.; Massini, M.; Simonucci, S.; Taioli, S. Ab Initio Calculation of the C1s Photoelectron Spectrum of C2 H2 . Nucl. Instr. Meth. Phys. Res. B 2004, 213, 65–70. (17) Ponzi, A.; Sapunar, M.; Angeli, C.; Cimiraglia, R.; Došlić, N.; Decleva, P. Photoionization of Furan from the Ground and Excited Electronic States. J. Chem. Phys 2016, 144, 084307.

(9) Brumboiu, I. E.; Prokopiou, G.; Kronik, L.; Brena, B. Valence Electronic Structure of Cobalt Phthalocyanine from an Optimally Tuned Range-Separated Hybrid Functional. J. Chem. Phys. 2017, 147, 044301.

(18) Tao, L.; McCurdy, C. W.; Rescigno, T. N. Grid-Based Methods for Diatomic Quantum Scattering Problems: A FiniteElement Discrete-Variable Representation in Prolate Spheroidal Coordinates. Phys. Rev. A 2009, 79, 012719.

(10) Gao, C.-Z.; Wopperer, P.; Dinh, P. M.; Suraud, E.; Reinhard, P.-G. On the Dynamics of Photo-Electrons in C60 . J. Phys. B 2015, 48, 105102.

(19) Larsson, H. R.; Bauch, S.; Sørensen, S.; Bonitz, M. Correlation Effects in StrongField Ionization of Heteronuclear Diatomic Molecules. Phys. Rev. A 2016, 93, 013426.

(11) Von Niessen, W.; Schirmer, J.; Cederbaum, L. S. Computational Methods for the One-Particle Green’s Function. Comp. Phys. Rep. 1984, 1, 57–125. (12) Schirmer, J.; Cederbaum, L. S.; Walter, O. New Approach to the One-Particle Green’s Function for Finite Fermi Systems. Phys. Rev. A 1983, 28, 1237–1259. (13) Kadanoff, L. P.; Baym, G. Quantum Statistical Mechanics: Green’s Function Method in Equilibrium and Nonequilibrium Problems; Addison-Wesley Publishing Co. Inc.: Redwood City, Calif, 1989. (14) Pohl, A.; Reinhard, P.-G.; Suraud, E. Towards Single-Particle Spectroscopy of Small Metal Clusters. Phys. Rev. Lett. 2000, 84, 5090–5093. (15) Bachau, H.; Cormier, E.; Decleva, P.; Hansen, J. E.; Martín, F. Applications of B -Splines in Atomic and Molecular Physics. Rep. Prog. Phys. 2001, 64, 1815– 1943.

(20) Dauth, M.; Graus, M.; Schelter, I.; Wießner, M.; Schöll, A.; Reinert, F.; Kümmel, S. Perpendicular Emission, Dichroism, and Energy Dependence in Angle-Resolved Photoemission: The Importance of The Final State. Phys. Rev. Lett. 2016, 117, 183001. (21) Marante, C.; Klinker, M.; Corral, I.; González-Vázquez, J.; Argenti, L.; Martín, F. Hybrid-Basis Close-Coupling Interface to Quantum Chemistry Packages for the Treatment of Ionization Problems. J. Chem. Theor. Comp. 2017, 13, 499–514. (22) De Giovannini, U.; Hübener, H.; Rubio, A. A First-Principles Time-Dependent Density Functional Theory Framework for Spin and Time-Resolved AngularResolved Photoelectron Spectroscopy in Periodic Systems. J. Chem. Theor. Comp. 2017, 13, 265–273. (23) Pickup, B. T. On the Theory of Fast Photoionization Processes. Chem. Phys. 1977, 19, 193–208.

ACS Paragon Plus Environment

4 CONCLUSIONS

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(24) Oana, C. M.; Krylov, A. I. Cross Sections and Photoelectron Angular Distributions in Photodetachment from Negative Ions Using Equation-of-Motion CoupledCluster Dyson Orbitals. J. Chem. Phys. 2009, 131, 124114. (25) Decleva, P.; Fronzoni, G.; Kivimäki, A.; Álvarez Ruiz, J.; Svensson, S. Shake-up Transitions in S 2p, S 2s and F 1s Photoionization of the SF6 Molecule. J. Phys. B 2009, 42, 055102. (26) Grell, G.; Bokarev, S. I.; Winter, B.; Seidel, R.; Aziz, E. F.; Aziz, S. G.; Kühn, O. Multi-Reference Approach to the Calculation of Photoelectron Spectra Including Spin-Orbit Coupling. J. Chem. Phys 2015, 143, 074104. (27) Kronik, L.; Kümmel, S. In First Principles Approaches to Spectroscopic Properties of Complex Materials; Di Valentin, C., Botti, S., Cococcioni, M., Eds.; Springer Berlin Heidelberg: Berlin, Heidelberg, 2014; Vol. 347; pp 137–191. (28) Jones, R. O. Density Functional Theory: Its Origins, Rise to Prominence, and Future. Rev. Mod. Phys. 2015, 87, 897–923. (29) Leang, S. S.; Zahariev, F.; Gordon, M. S. Benchmarking the Performance of TimeDependent Density Functional Methods. J. Chem. Phys. 2012, 136, 104101. (30) Baer, R.; Livshits, E.; Salzner, U. Tuned Range-Separated Hybrids in Density Functional Theory. Annu. Rev. Phys. Chem. 2010, 61, 85–109. (31) Kümmel, S.; Kronik, L. OrbitalDependent Density Functionals: Theory and Applications. Rev. Mod. Phys. 2008, 80, 3–60. (32) Savin, A. Recent Advances in Computational Chemistry; World Scientific, 1995; Vol. 1; pp 129–153. (33) Refaely-Abramson, S.; Sharifzadeh, S.; Govind, N.; Autschbach, J.; Neaton, J. B.;

16

Page 16 of 20

Baer, R.; Kronik, L. Quasiparticle Spectra from a Nonempirical Optimally Tuned Range-Separated Hybrid Density Functional. Phys. Rev. Lett. 2012, 109, 226405. (34) Gozem, S.; Gunina, A. O.; Ichino, T.; Osborn, D. L.; Stanton, J. F.; Krylov, A. I. Photoelectron Wave Function in Photoionization: Plane Wave or Coulomb Wave? J. Phys. Chem. Lett. 2015, 6, 4532–4540. (35) Walter, M.; Häkkinen, H. Photoelectron Spectra from First Principles: From the Many-Body to the Single-Particle Picture. New J. Phys. 2008, 10, 043018. (36) Humeniuk, A.; Wohlgemuth, M.; Suzuki, T.; Mitrić, R. Time-Resolved Photoelectron Imaging Spectra from Non-Adiabatic Molecular Dynamics Simulations. J. Chem. Phys. 2013, 139, 134104. (37) Onida, G.; Reining, L.; Rubio, A. Electronic Excitations: Density-Functional versus Many-Body Green’s-Function Approaches. Rev. Mod. Phys. 2002, 74, 601– 659. (38) Hirata, S.; Head-Gordon, M. TimeDependent Density Functional Theory within the Tamm–Dancoff Approximation. Chem. Phys. Lett. 1999, 314, 291 – 299. (39) Åberg, T. Theory of X-Ray Satellites. Phys. Rev. 1967, 156, 35–41. (40) Ipatov, A.; Cordova, F.; Doriol, L. J.; Casida, M. E. Excited-State SpinContamination in Time-Dependent Density-Functional Theory for Molecules with Open-Shell Ground States. J. Mol. Struct. (THEOCHEM) 2009, 914, 60–73. (41) Wittbrodt, J. M.; Schlegel, H. B. Some Reasons Not to Use Spin Projected Density Functional Theory. J. Chem. Phys. 1996, 105, 6574–6577.

ACS Paragon Plus Environment

4 CONCLUSIONS

Page 17 of 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(42) Peach, M. J. G.; Tozer, D. J. Overcoming Low Orbital Overlap and Triplet Instability Problems in TDDFT. J. Phys. Chem. A 2012, 116, 9783–9789. (43) Casida, M.; Huix-Rotllant, M. Progress in Time-Dependent Density-Functional Theory. Annu. Rev. Phys. Chem. 2012, 63, 287–323. (44) Peach, M. J. G.; Williamson, M. J.; Tozer, D. J. Influence of Triplet Instabilities in TDDFT. J. Chem. Theor. Comput. 2011, 7, 3578–3585. (45) Leininger, T.; Stoll, H.; Werner, H.-J.; Savin, A. Combining Long-Range Configuration Interaction with Short-Range Density Functionals. Chem. Phys. Lett. 1997, 275, 151–160. (46) Yanai, T.; Tew, D. P.; Handy, N. C. A New Hybrid Exchange-Correlation Functional Using the Coulomb-Attenuating Method (CAM-B3LYP). Chem. Phys. Lett. 2004, 393, 51–57. (47) Tawada, Y.; Tsuneda, T.; Yanagisawa, S.; Yanai, T.; Hirao, K. A Long-RangeCorrected Time-Dependent Density Functional Theory. J. Chem. Phys. 2004, 120, 8425–8433. (48) Song, J.-W.; Hirosawa, T.; Tsuneda, T.; Hirao, K. Long-Range Corrected Density Functional Calculations of Chemical Reactions: Redetermination of Parameter. J. Chem. Phys. 2007, 126, 154105. (49) Refaely-Abramson, S.; Baer, R.; Kronik, L. Fundamental and Excitation Gaps in Molecules of Relevance for Organic Photovoltaics from an Optimally Tuned Range-Separated Hybrid Functional. Phys. Rev. B 2011, 84, 075144. (50) Borpuzari, M. P.; Kar, R. A New Nonempirical Tuning Scheme with Single SelfConsistent Field Calculation: Comparison with Global and IP-Tuned RangeSeparated Functional. J. Comput. Chem. 2017, 38, 2258–2267.

17

(51) Stein, T.; Kronik, L.; Baer, R. Reliable Prediction of Charge Transfer Excitations in Molecular Complexes Using TimeDependent Density Functional Theory. J. Amer. Chem. Soc. 2009, 131, 2818–2820. (52) Bokareva, O. S.; Grell, G.; Bokarev, S. I.; Kühn, O. Tuning Range-Separated Density Functional Theory for Photocatalytic Water Splitting Systems. J. Chem. Theor. Comput. 2015, 11, 1700–1709. (53) Bokarev, S. I.; Bokareva, O. S.; Kühn, O. A Theoretical Perspective on Charge Transfer in Photocatalysis. The Example of Ir-Based Systems. Coord. Chem. Rev. 2015, 304-305, 133–145. (54) Körzdörfer, T.; Sears, J. S.; Sutton, C.; Brédas, J.-L. Long-Range Corrected Hybrid Functionals for π-Conjugated Systems: Dependence of the RangeSeparation Parameter on Conjugation Length. J. Chem. Phys. 2011, 135, 204107. (55) Karolewski, A.; Kronik, L.; Kümmel, S. Using Optimally Tuned Range Separated Hybrid Functionals in Ground-State Calculations: Consequences and Caveats. J. Chem. Phys. 2013, 138, 204115. (56) Srebro, M.; Autschbach, J. Tuned RangeSeparated Time-Dependent Density Functional Theory Applied to Optical Rotation. J. Chem. Theor. Comp. 2012, 8, 245–256. (57) Srebro, M.; Autschbach, J. Does a Molecule-Specific Density Functional Give an Accurate Electron Density? The Challenging Case of the CuCl Electric Field Gradient. J. Phys. Chem. Lett. 2012, 3, 576–581. (58) Autschbach, J.; Srebro, M. Delocalization Error and “Functional Tuning” in Kohn–Sham Calculations of Molecular Properties. Acc. Chem. Res. 2014, 47, 2592–602.

ACS Paragon Plus Environment

4 CONCLUSIONS

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(59) Livshits, E.; Baer, R. A Density Functional Theory for Symmetric Radical Cations from Bonding to Dissociation. Phys. Chem. Chem. Phys. 2007, 9, 2932– 41. (60) Stein, T.; Kronik, L.; Baer, R. Prediction of Charge-Transfer Excitations in Coumarin-Based Dyes Using a RangeSeparated Functional Tuned from First Principles. J. Chem. Phys. 2009, 131, 244119. (61) Mori-Sánchez, P.; Cohen, A. J.; Yang, W. Localization and Delocalization Errors in Density Functional Theory and Implications for Band-Gap Prediction. Phys. Rev. Lett. 2008, 100, 146401. (62) Iikura, H.; Tsuneda, T.; Yanai, T.; Hirao, K. A Long-Range Correction Scheme for Generalized-Gradient-Approximation Exchange Functionals. J. Chem. Phys 2001, 115, 3540–3544. (63) Bauernschmitt, R.; Ahlrichs, R. Stability Analysis for Solutions of the Closed Shell Kohn–Sham Equation. J. Chem. Phys. 1996, 104, 9047–9052. (64) 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. (65) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. J. A.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.;

18

Page 18 of 20

Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, O.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, Revision D.01 ; 2009. (66) Valiev, M.; Bylaska, E. J.; Govind, N.; Kowalski, K.; Straatsma, T. P.; Dam, H. J. J. V.; Wang, D.; Nieplocha, J.; Apra, E.; Windus, T. L.; de Jong, W. A. NWChem: A Comprehensive and Scalable OpenSource Solution for Large Scale Molecular Simulations. Comp. Phys. Comm. 2010, 181, 1477–1489. (67) Landau, L. D.; Lifshitz, E. M. Quantum Mechanics: Non-Relativistic Theory; Pergamon: Oxford, 1977. (68) Gozem, S.; Krylov, A. I. ezDyson. 2017. (69) Dauth, M.; Caruso, F.; Kümmel, S.; Rinke, P. Piecewise Linearity in the GW Approximation for Accurate Quasiparticle Energy Predictions. Phys. Rev. B 2016, 93, 121115. (70) Vetter, K.; Proch, S.; Ganteför, G. F.; Behera, S.; Jena, P. Hydrogen Mimicking the Properties of Coinage Metal Atoms in Cu and Ag Monohydride Clusters. Phys. Chem. Chem. Phys. 2013, 15, 21007. (71) Ho, J.; Ervin, K. M.; Lineberger, W. C. Photoelectron Spectroscopy of Metal Cluster Anions: Cu n , Ag n , and Au n . J. Chem. Phys. 1990, 93, 6987–7002. (72) Taylor, K. J.; Pettiette-Hall, C. L.; Cheshnovsky, O.; Smalley, R. E. Ultraviolet

ACS Paragon Plus Environment

4 CONCLUSIONS

Page 19 of 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Photoelectron Spectra of Coinage Metal Clusters. J. Chem. Phys. 1992, 96, 3319– 3329. (73) Leopold, D. G.; Ho, J.; Lineberger, W. C. Photoelectron Spectroscopy of Massselected Metal Cluster Anions. I. Cu n , n =1–10. J. Chem. Phys. 1987, 86, 1715– 1726. (74) Winter, B.; Weber, R.; Widdra, W.; Dittmar, M.; Faubel, M.; Hertel, I. V. Full Valence Band Photoemission from Liquid Water Using EUV Synchrotron Radiation. J. Phys. Chem. A 2004, 108, 2625–2632. (75) Liu, S.-Y.; Alnama, K.; Matsumoto, J.; Nishizawa, K.; Kohguchi, H.; Lee, Y.-P.; Suzuki, T. He I Ultraviolet Photoelectron Spectroscopy of Benzene and Pyridine in Supersonic Molecular Beams Using Photoelectron Imaging. J. Phys. Chem. A 2011, 115, 2953–2965. (76) Richardson, N. V.; Weinberger, P. The Electronic Structure of the S8 Molecule. J. El. Spec. Rel. Phen. 1975, 6, 109–116. (77) Boschi, R.; Schmidt, W. The Photoelectron Spectrum and Structure of Sulfur in the Gas Phase at 140 ◦ C. Inorg. Nuc. Chem. Lett. 1973, 9, 643–648. (78) Berkowitz, J.; Lifshitz, C. Photoionization of High-Temperature Vapors. II. Sulfur Molecular Species. J. Chem. Phys. 1968, 48, 4346–4350. (79) Rosinger, W.; Grade, M.; Hirschwald, W. Electron Impact Induced Excitation Processes Involving the Sulfur Clusters S2 to S8 . Ber. Bunsenges. Phys. Chem. 1983, 87, 536–542. (80) Berkowitz, J.; Chupka, W. A. Vaporization Processes Involving Sulfur. J. Chem. Phys. 1964, 40, 287–295. (81) Hagemann, R. Determination de La Chaleur de Formation de S2 O Par Spectrometrie de Masse. C. R. Chim. 1962, 255, 1102.

19

(82) Rau, H.; Kutty, T.; Guedes De Carvalho, J. Thermodynamics of Sulphur Vapour. J. Chem. Thermodyn. 1973, 5, 833–844. (83) Engel, N.; Bokarev, S. I.; Moguilevski, A.; Raheem, A. A.; Al-Obaidi, R.; Möhle, T.; Grell, G.; Siefermann, K. R.; Abel, B.; Aziz, S. G.; Kühn, O.; Borgwardt, M.; Kiyan, I. Y.; Aziz, E. F. Light-Induced Relaxation Dynamics of the Ferricyanide Ion Revisited by Ultrafast XUV Photoelectron Spectroscopy. Phys. Chem. Chem. Phys. 2017, 19, 14248–14255. (84) Raheem, A. A.; Wilke, M.; Borgwardt, M.; Engel, N.; Bokarev, S. I.; Grell, G.; Aziz, S. G.; Kühn, O.; Kiyan, I. Y.; Merschjann, C.; Aziz, E. F. Ultrafast Kinetics of Linkage Isomerism in Na2 [Fe(CN)5 NO] Aqueous Solution Revealed by Time-Resolved Photoelectron Spectroscopy. Struct. Dyn. 2017, 4, 044031. (85) Moguilevski, A.; Wilke, M.; Grell, G.; Bokarev, S. I.; Aziz, S. G.; Engel, N.; Raheem, A. A.; Kühn, O.; Kiyan, I. Y.; Aziz, E. F. Ultrafast Spin Crossover in [FeII (bpy)3 ]2+ : Revealing Two Competing Mechanisms by Extreme Ultraviolet Photoemission Spectroscopy. ChemPhysChem 2017, 18, 465–469. (86) Myneni, H.; Casida, M. E. On the Calculation of ∆hS 2 i for Electronic Excitations in Time-Dependent Density-Functional Theory. Comput. Phys. Commun. 2017, 213, 72–91. (87) Yamaguchi, K.; Jensen, F.; Dorigo, A.; Houk, K. A Spin Correction Procedure for Unrestricted Hartree-Fock and MøllerPlesset Wavefunctions for Singlet Diradicals and Polyradicals. Chem. Phys. Lett. 1988, 149, 537–542. (88) Hratchian, H. P. Communication: An Efficient Analytic Gradient Theory for Approximate Spin Projection Methods. J. Chem. Phys. 2013, 138, 101101.

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table of Contents Graph 88x35mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 20 of 20