Benchmark of Excitation Energy Shifts from Frozen-Density

Jun 15, 2018 - Benchmark of Excitation Energy Shifts from Frozen-Density Embedding ... [email protected]., *E-mail: [email protected]...
0 downloads 0 Views 1MB Size
Subscriber access provided by University of South Dakota

Quantum Electronic Structure

Benchmark of excitation energy shifts from Frozen-Density Embedding Theory: introduction of a density-overlap based applicability threshold Alexander Zech, Niccolo Ricardi, Stefan Prager, Andreas Dreuw, and Tomasz A. Wesolowski J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00201 • Publication Date (Web): 15 Jun 2018 Downloaded from http://pubs.acs.org on June 17, 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 43 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

Benchmark of excitation energy shifts from Frozen-Density Embedding Theory: introduction of a density-overlap based applicability threshold Alexander Zech,∗,† Niccol`o Ricardi,∗,† Stefan Prager,‡ Andreas Dreuw,‡ and Tomasz A. Wesolowski∗,† Department of Physical Chemistry, University of Geneva, Geneva, Switzerland, and Interdisciplinary Center for Scientific Computing, University of Heidelberg, Heidelberg, Germany E-mail: [email protected]; [email protected]; [email protected]

Abstract We present a thorough investigation of the errors in results obtained with the combination of Frozen-Density Embedding Theory and the Algebraic Diagrammatic Construction Scheme for the polarization propagator of second order (FDE-ADC(2)). The study was carried out on a set of 52 intermolecular complexes with varying interaction strength, each consisting of a chromophore of fundamental interest and a few small molecules in its environment. The errors emerging in Frozen-Density Embedding Theory based methods originate from: (a) the solver of the quantum many body problem used to obtain the embedded wavefunction (Ψemb A ), (b) the approximation for the explicit density functional for the embedding potential, and (c) the choice of the density ∗

To whom correspondence should be addressed Department of Physical Chemistry, University of Geneva, Geneva, Switzerland ‡ Interdisciplinary Center for Scientific Computing, University of Heidelberg, Heidelberg, Germany †

1

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

Page 2 of 43

representing the environment (ρB (~r)). The present work provides a comprehensive analysis of the errors in the excitation energies based on the last two factors. Furthermore a density-overlap based parameter is proposed to be used as an a priori criterion of applicability.

1

Introduction

Recently, multi-level simulations using the idea of density embedding have become increasingly popular. A common feature of such methods is the charge density dependent operator ˆ , to its environment B in the vˆemb coupling subsystem A, defined by its Hamiltonian H A following Schr¨odinger-like equation: 

 ˆ HA + υˆemb Ψemb = Ψemb A A

(1)

Such methods are based on either Frozen-Density Embedding Theory (FDET) or on the energy-error compensation ansatz. FDET was originally formulated by Wesolowski and Warshel for embedding a non-interacting reference system 1 and subsequently extended to other quantum-mechanical descriptors of the embedded subsystem: a) the wavefunction of the multi-determinantal form 2 and b) the one-particle spinless density matrix and the corresponding energy functional. 3 Error compensation Ansatz based methods were pioneered by Carter and collaborators who introduced it for embedding a system described by an explicitly interacting Hamiltonian. 4 Both approaches (FDET and the energy-error compensation Ansatz) feature a multiplicative embedding operator, which has the form of a local one-particle potential:

υˆemb =

N X i

υemb (~r)δ(~r − ~ri ).

2

ACS Paragon Plus Environment

(2)

Page 3 of 43 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

For FDET the embedding potential reads: Z vemb [ρA , ρB , vB ](~r) = vB (~r) +

nad δExc,T [ρA , ρB ] ρ (~r 0 ) B 0 d~r 0 + ~r − ~r δρA (~r)

(3)

The first two terms in the potential (Eq. 3) represent the classical electrostatic potential generated by the environment. The last term,

nad vxc,T [ρA , ρB ](~r) =

nad [ρA , ρB ] δTsnad [ρA , ρB ] δExc + δρA (~r) δρA (~r)

(4)

has a quantum mechanical origin and is a bi-functional uniquely determined by the pair ρA (~r) and ρB (~r). This implicit functional is defined by virtue of the constrained search formalism. 2 Eq. 3 is exact without additional non-local terms (a suggestion which is made sometimes in the literature 5–7 ) within the formal framework of FDET. Both the exact solutions of FDET equations available for analytically solvable systems 8 as well as practical numerical examples discussed in detail by Unsleber et al. 9 show that additional non-local components of the embedding operator are not needed. nad [ρA , ρB ](~r) is frequently approximated by means of an exIn practical simulations, vxc,T

plicit bi-functional. The domain of currently used inexpensive semi-local approximations for the non-electrostatic component of the FDET embedding potential is limited to cases where the overlap between ρA (~r) and ρB (~r) is small. In practice these approximations are applicable only in the absence of covalent bonds between subsystems, 10–12 although there exist approaches which address this issue by means of going beyond the formal framework of FDET 13 or constructing approximations for the embedding potential without the use of an nad explicit density functional to approximate the potential vxc,T [ρA , ρB ] (~r). In the latter group nad of approaches, an approximation to the potential vxc,T [ρA , ρB ] (~r) is constructed using nu-

merical inversion of the Kohn-Sham equation (see the original works 12,14–16 from the groups of Carter, Roncero, Miller, or Jacob or our recent review 17 ).

3

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

Page 4 of 43

Due to the inhomogeneity of non-additive functionals,

nad Exc,T [ρA , ρB ]

6=

Z

nad ρA (~r) vxc,T [ρA , ρB ] (~r)d~r,

(5)

the Lagrange multiplier () obtained from Eq. 1 is not simply related to the electronic energy of the whole system. The latter reads: D E ˆ nuc EW F [ΨA , ρB ] = ΨA H EAB Ψ [ρA ] + Jint [ρA , ρB ] A A + VB nad +Exc,T [ρA , ρB ] + EvHK [ρB ] + VAnuc [ρB ] B

(6)

E ˆ [ρB ] + VAnuc [ρB ] 6 = ΨA HA + υˆemb ΨA + EvHK B D

[ρB ] denotes the energy of the isolated subsystem B as obtained from the The term EvHK B Hohenberg-Kohn formulation of DFT. For the complete description of the formalism see Ref. 2. Here, we underline that FDET features the embedded wavefunction (Ψemb A ) and the electron density of the environment (ρB (~r)) as independent variables, and that it leads to self-consistent expressions for the total energy, embedding potential and optimal embedded wavefunction. Both the environment density (used as an independent variable in FDET) and the density obtained from Eq. 1 are non-negative functions (ρB (~r) ≥ 0 and ρA (~r) ≥ 0). FDET yields the exact total density (ρItot (~r)) for a given electronic state (I) only for ρB (~r) such that  ∀~r ρItot (~r) > ρB (~r) .

(7)

By virtue of the Hohenberg-Kohn theorem, FDET can yield only the upper limit of the ground state energy of the whole system if ρB (~r) violates the above inequality. Several reviews are available in the literature which deal with various flavors of the density based embedding idea. 18–22 The present work concerns methods based on the formal

4

ACS Paragon Plus Environment

Page 5 of 43 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

nad nad framework of FDET, in which the exact functionals vxc,T [ρA , ρB ] (~r) and Exc,T [ρA , ρB ] are

approximated by means of explicitly density dependent analytical expressions:

nad nad vxc,T [ρA , ρB ] (~r) ≈ v˜xc,T [ρA , ρB ] (~r)

(8)

nad nad Exc,T [ρA , ρB ] ≈ E˜xc,T [ρA , ρB ] .

(9)

and

Careful analysis of the results of FDET-based simulations published in the literature nad indicates that the accuracy of the approximation for vxc,T [ρA , ρB ](~r) deteriorates with in-

creasing overlap between ρA (~r) and ρB (~r). So far, no systematic analysis has been made concerning the relation between the accuracy of the observables derived from FDET and the magnitude of the overlap. The errors in the total energy obtained for real cheminad cal systems using a given approximation for E˜xc,T [ρA , ρB ] (and the corresponding potential nad v˜xc,T [ρA , ρB ](~r) =

˜ nad [ρ ,ρ ] E xc,T A B ) δρA (~ r)

do not provide a direct insight into the accuracy of the po-

tential. The error in the total energy results from accumulated errors in the functional and its functional derivative. The errors in these two quantities, however, were shown to not be correlated. 8 A similar issue was discussed recently for the functionals used in the DFT formulation of the quantum many body problem. 23 The analysis of errors in the density is more useful. Results reported in Refs. 24,25 clearly show that errors in the density increase with decreasing intermolecular distance. Another body of evidence for the correlation between these errors comes from studies on model systems for which FDET equations can be solved analytically. 26 Inspired by these observations, we chose to investigate such correlations in more detail. The magnitude of the overlap relates to well-defined quantities such as the intermolecular distance or to qualitative descriptors characterizing the type of interaction between the chromophore and its environment in the ground and excited states (hydrogen bonding character, type of excitation, etc.). To this end, we investigate environment-induced shifts in

5

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

the excitation energies for various types and strengths of chromophore-environment interactions in several intermolecular complexes. In most of them, the chromophore is hydrogen bonded to one or more molecules represented by means of vB (~r) and ρB (~r) in FDET. The chosen complexes differ in magnitude and topology of the overlap of ρA (~r) and ρB (~r). Such systems represent a typical domain of applicability of approximate FDET-based methods. nad Attempts to use explicit semi-local approximations for v˜xc,T [ρA , ρB ](~r) for embedded species

covalently bound to the environment were shown to fail. 11 The present work does not cover these applications which we believe lie outside of the domain of practical applications of such approximations. In view of these considerations, we developed and tested an overlap-based diagnostic parameter for FDET calculations. We should also point out, that the magnitude of the overlap of ρA (~r) and ρB (~r) – regardless of how it is measured – cannot be expected to be the sole factor responsible for the error in the embedding calculations. The effect of the environment on electronic excitations is frequently interpreted by means of quantities featured in the theory of intermolecular interactions (see the seminal work by McRae 27 ). In methods in which the polarization, mutual induction, overlap, and dispersion contributions are explicitly included, it is possible to determine the significance of such contributions. These contributions are known to be strongly system-dependent. Unfortunately, the one-toone correspondence between the aforementioned contributions to the interaction energy and the components of the FDET embedding potential does not exist. Nevertheless, a strong system-dependency of the FDET errors can be expected as well. It originates, however, from other sources. For instance, the quality of the approximation for the non-additive kinetic potential is irrelevant at small overlaps between ρA (~r) and ρB (~r) because the functional vTnad [ρA , ρB ](~r) (both exact as well as approximate as in the present work) disappears at zero overlap. The quality of the approximation for the non-additive kinetic potential is, on the other hand, crucial for overlaps occuring at equilibrium geometries of the interacting molecules.

6

ACS Paragon Plus Environment

Page 6 of 43

Page 7 of 43 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

For the numerical variant of FDET applied in this work, we have chosen the recent implementation of Linearized FDET 28,29 combined with the Algebraic Diagrammatic Construction (ADC(2)) scheme of the polarization propagator of second order for the calculation of excitation energies. We have used this implementation in view of the numerical soundness of ADC(2) (see Ref. 30) and the recently shown excellent performance of this combination 31,32 in describing environment induced shifts of vertical excitation energies in small intermolecular complexes. Thereby we also aim to extend the database of systems for which the performance of such a combination is tested. We point out that possible conclusions connad cerning applicability of a given approximation for the functional vxc,T [ρA , ρB ](~r) emerging

from such studies as the present one concern also any other method applying approximations nad [ρA , ρB ](~r) including all possible combinations of the FDET embedding potential with for vxc,T

correlated ab intio quantum mechanical methods (see for example Refs. 32–38) as well as methods following the LR-TDDFT approach to excitation energies 39 and its extensions. 40,41 Due to the use of a large dataset, the analysis of ρB -dependence of the excitation energies in this work extends and complements the previous works by ourselves 42 and by Daday et al. 37

2

Methodology

2.1

Molecular test set

A representative set of optimized supermolecular structures was devised to both thoroughly test the FDE-ADC protocol and to investigate the correlation between its error and the newly designed overlap parameter Ω (see Sec. 2.3). All of the intermolecular clusters considered, and their main characteristics, are summarized in Table 1. The test set focuses on hydrogen bonded clusters due to their importance in supramolecular chemistry (particularly in biological systems), their considerable solvatochromic shifts, and the wide range of interaction strengths. The selected set includes structures with: 7

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

ˆ The embedded system acting as hydrogen bond donor, acceptor, or both ˆ Environments differing in size, complexity, and amount of intra-environment interac-

tions ˆ System-environment interactions involving different atoms and functional groups

The choice of systems and their environment is a tradeoff between maintaining a feasible size regarding the computational cost and selecting chromophores of chemical interest. The systems were chosen in order to posses mainly localized transitions which can be treated with FDET, i.e. where neither charge transfer between system A and environment B nor low lying excitations on the environment are relevant. Such excitations may interact with those of the system of interest via excitonic coupling, and this would not be described due to the neglect of the dynamic response of the environment. In the presence of such couplings, they can be explicitly simulated using the extension of the subsystem DFT, 43,44 which is the ground-state formalism, to excited states following the LR-TDDFT approach. 40,41 Systems 1a, 1c-e, 1g-i, 1k-o, and 1t were chosen because they were subject of previous studies 42,45–47 where they were assigned to correspond to experimental measurements. In this subset mono-solvated clusters always have the hydroxy group acting as hydrogen bond donor, doubly coordinated clusters have a nitrogen acting as acceptor and a hydroxy group acting as donor, and triply coordinated clusters have a third molecule indirectly bound through the two other ones. Another subset of geometries was developed in order to test different coordination patterns on 1. This subset includes mono-solvated structures with the nitrogen acting as acceptor (1b, 1f, and 1j) and π systems in the environment (1p-z) with a different type of interaction with the chromophore. A third smaller subset is based on the test systems (8 and 9) used in previous work on FDE-ADC. 31 Two structures from the S66 benchmark set 48 were selected as examples for π π (8b) and NH π (12a) interactions. The remaining structures were developed starting from chromophores of interest such as xanthine (2) and its deprotonated anion (3), diketopyrrolopyrrole (7), pyridiniumyl ben8

ACS Paragon Plus Environment

Page 8 of 43

Page 9 of 43 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

O

5

4 3

6

6

HN 1

5

6

1

H N

N

8 8

1

2

OH

8

O

1

N

7

7 2

3

4

N H

H2N

N

2

N

11

N

3

4

N

R

R=H R=Me

3

4 5

1N

H N 1

5

4

2 3

1

3

H

5

3

O

6

9

7

O

8

1

N

5

6

2

4 4

8

7

2

10

O

H

6

6

2

6

6

O

O

NH

5

1

8

3

4

NH

HN

9

4 5

O

2

7

6

5

10

N+

N9

9

O 1

4

3

2 3=(2)-

2

7

N

5

11

NH2

12

Figure 1: Chromophores considered in this work are (1) 7-hydroxyquinoline, (2) xanthine, (3) xanthinide, (4) 2-aminopurine, (5) 7-methyl-2-aminopurine, (6) pyridiniumyl benzimidazolide, (7) diketopyrrolopyrrole, (8) uracil, (9) benzaldehyde, (10) 4-dimethylaminopyridine, (11) 7-Amino-4-methylcoumarin and (12) benzene. zimidazolide (6), aminopurine (4) and methylaminopurine (5), coumarin 120 (11), and dimethylaminopyridinium cation (10). These chromophores were coordinated with dipolar or H-bonding small molecules in different numbers and positions. It should be noted that the resulting set, albeit focused on hydrogen bonds, spans a very large range of the most important properties: bond length, and consequently interaction strength, dipole moment and polarizability of the chromophore and of the environment. Furthermore, neutral, positively, and negatively charged chromophores and environments were tested.

9

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

Page 10 of 43

Table 1: Chromophore environment clusters considered in this work. D/A corresponds to the donor or acceptor atom(s) of the chromophore involved in the interaction with the environment. Label

Environment

1 7-hydroxyquinoline 1a H2O 1b H2O 1c 2 H2O 1d 3 H2O 1e NH3 1f NH3 1g 2 NH3 1h 3 NH3 1i MeOH 1j MeOH 1k 2 MeOH 1l 3 MeOH 1m NH3 NH3 H2O 1n NH3 H2O NH3 1o NH3 H2O H2O 1p formamide 1q HCOO– 1r formimidamide 1s guanidine 1t HCOOH 1u HCOOH 1v HCOOH 1w HCOOH 1x HCOOH 1y HCOOH 1z HCOOH 2 Xanthine 2a MeCN 2b MeSH 2c MeOH 2d EtOH 2e F3C COO– 2f 3 CH3F 2g 2h 1 2

2.2

3 BrF 3 NH4+

D/A

Label

H-7 N-1 N-1, H-7 N-1, H-7 H-7 N-1 N-1, H-7 N-1, H-7 H-7 N-1 N-1, H-7 N-1, H-7 N-1, H-7 N-1, H-7 N-1, H-7 N-1, H-8 H-7, H-8 H-7 H-7 N-1, H-8 H-7, H-8 N-1 (⊥)2 H-7 (⊥)2 N-1, H-2 H-6, O-7 H-7, C-8

3 Xanthinyl anion 3a HCOOH N-7, H-8 4 Aminopurine 4a H2O N-3, H-9 4b H2O N-1, H-2 4c 2 H2O N-1, H-2, N-3, H-9 4d 2 H2O N-1, H-2a, H-2b, N-3 5 Methyl aminopurine 5a H2O H-2, N-3 6 Pyridiniumyl benzimidazolide 6a 4 NH3 N-1, N-6 6b 5 H2O N-1, N-6, H-7, H-10 6c 2 HCOOH N-1, N-6, H-7, H-11 7 Diketopyrrolopyrrole 7a 2 HCOOH O-1, H-2, O-4, H-5 7b 6 H2O O-1, H-2, O-4, H-5 7c 10 H2O O-1, H-2, O-4, H-5 8 Uracil 8a 5 H2O H-1, O-2, H-3, O-4, H-5 8b1 pyridine π-system 9 Benzaldehyde 9a 2 H2O O-1 10 Dimethylaminopyridinium cation 10a 4 H2O H-1, H-3, H-4a, H-4b, H-5 11 Coumarin 120 11a 6 H2O O-1, O-2, H-7a, H-7b 12 Benzene 12a1 acetamide π-system

O-6, O-6, O-6, O-6, H-7, H-1, O-6, O-2, O-2,

H-7 H-7 H-7 H-7 H-8 O-2, H-7, O-6, O-6,

Environment

D/A

H-3, N-9 N-9 N-9

Taken from Ref. 48 The symbol ⊥ indicates the perpendicular alignment of the molecular planes of chromophore and environment.

Assignment of the excitations

In order to analyze the errors in the excitation energy shifts obtained with a certain embedding protocol the same excited state was calculated in three different manners: ADC(2) for 10

ACS Paragon Plus Environment

Page 11 of 43 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

the isolated chromophore, ADC(2) for the complex, and lastly FDE-ADC(2). Since in these calculations the order of states can be different, the states need to be identified and assigned. A simple assignment can be obtained by comparing oscillator strengths and considering the standard range of solvatochromic shifts. In the majority of the cases a more detailed analysis based on MOs and transition amplitudes was sufficient, while in few difficult instances the correspondence was established by means of tools for wavefunction analysis (libwfa) 49–51 for difference densities, natural transition orbitals, 52,53 detachment, and attachment densities. From these calculations singlet excited states with predominant nπ ∗ or ππ ∗ character could be clearly distinguished and will be henceforth referred to as nπ ∗ and ππ ∗ states, respectively. Quasi-degeneracy may lead to state-mixing and loss of the correspondence between excited states in the three types of calculations. This affects oscillator strengths significantly more than excitation energies because the mixing of dark and bright states can result in socalled ”electronic intensity borrowing“. Concerning Rydberg states, they were not included in the present analysis (see also the discussion of delocalization below). For excitations retained in the set of assigned excitations, one can expect that the errors in the excitation energies are related not only to the magnitude of the overlap but also to the degree of localization of the excited state. Two extreme cases illustrate the relevance of delocalization for expected errors in FDET based calculations. The first case is an excitation (charge transfer from the environment to the chromophore for instance), which leads to a  decrease of electron density somewhere in the environment: ∃~r ∈ B ρItot (~r) − ρ0tot (~r) < 0 . If the density of the isolated environment is used as ρB (~r), it is very likely that the inequality given in Eq. 7 is violated locally somewhere in the environment. As a consequence, FDET cannot yield the correct total density and hence the energy. FDET using ground-state density of the isolated environment as ρB (~r) is, therefore, doomed to fail for such excitanad tions regardless of the quality of the approximation used for vxc,T [ρA , ρB ] (~r). The other

extreme case is an excitation (charge transfer from the chromophore to the environment

11

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

Page 12 of 43

for instance), which leads to an increase of electron density somewhere in the environment:  ∃~r ∈ B ρItot (~r) − ρ0tot (~r) > 0 . FDET does not have any fundamental problem in describing such an excitation, but FDET-based calculations using certain classes of approximations for nad vxc,T [ρA , ρB ] (~r) are expected to be less reliable. Most of the approximations used in practice

(and in the present work) comprise the LDA component (derived from the Thomas-Fermi functional) which is repulsive. Studies on model systems show that this repulsive character is frequently overestimated. 8,54 Using such an approximation in Eq. 1 results in an excessive localization of ρA (~r). Although the excitations in our data set do not include cases of clear charge-transfer between the chromophore and the environment, a partial delocalization of the states always takes place. For this reason, we also investigated the relation between the errors of the excitation energies and the magnitude of the delocalization of the excited state. Each excited state obtained from a supermolecular ADC(2) calculation with a number of amplitudes n is characterized by means of a parameter rw (reduced weight ratio). n X (via )j rw = , a (w ) j i j

(10)

Here wia represents the weight for an MO transition i → a, whereas the reduced weight via is defined as: via =

 1X 2 cµ,i + c2µ,a · wia 2 µ∈A

(11)

in which cµ,i and cµ,a are the MO coefficients in the orthonormal basis of occupied and virtual orbitals respectively, spanning only atomic basis functions of subsystem A. Thus, an excitation fully localized on subsystem A would yield rw = 1, while larger MO coefficients for basis functions on the environment would decrease the reduced weight ratio to the point of rw = 0, which represents a state fully localized on the environment. We would like to stress that this parameter serves as a diagnostic tool and is not useful for practical calculations as it requires a computationally expensive ADC calculation on the full complex. Furthermore, the parameter rw is basis set dependent and will gradually lose 12

ACS Paragon Plus Environment

Page 13 of 43 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

interpretability if the number of diffuse basis functions in the basis set increases. A somewhat similar measure for the extent of delocalization already exists and was developed on the basis of the transition density matrix. 55–58

2.3

Measures of the overlap

In principle, various quantities can be used as a measure of the overlap of ρA (~r) and ρB (~r) . The simplest would be just an integral of the form: Z Ωk =

with k ≥ 1 or

Z Ωp =

(ρA (~r) · ρB (~r))k d~r

p

(ρA (~r) + ρB (~r)) d~r −

Z

ρpA (~r)d~r

(12)



Z

ρpB (~r)d~r

(13)

with p ≥ 1. The above measures are extensive, which is an undesired property for an universal applicability criterion of FDET. The following example demonstrates why the applicability criterion should rather not be based on an extensive measure of the overlap: let us consider a dimer consisting of two molecules (A, B) and some partitioning of the total electron density such that ρA (~r) is localized mainly near molecule A and ρB (~r) near molecule B which leads to some value of Ωp . Now, consider a larger system consisting of two such dimers A-B and A’-B’, but infinitely separated. Obviously, the value of Ωp calculated for ρdimer (~r) = ρA (~r) + ρ0A (~r) A and ρdimer (~r) = ρB (~r) + ρ0B (~r) will be larger by a factor of two despite the fact that the errors B in the embedding potential for the two identical dimers are independent of each other. The same observation concerns Ωk . However, a useful applicability criterion should be based on a system independent threshold. This can be easily achieved by normalization of Ωp . Additionally chosing p = 5/3 in Ωp leads to, up to a multiplicative constant, the Thomas-Fermi density functional for the kinetic energy. The final expression for the overlap measure used

13

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

Page 14 of 43

in the present analysis is thus: R p R ρA (~r)d~r + ρpB (~r)d~r Ω=1− R (ρA (~r) + ρB (~r))p d~r

(14)

with p = 5/3. For any pair ρA (~r) and ρB (~r), Ω is positive and disappears for non-overlapping ρA (~r) and ρB (~r). We also point out that this measure is dimensionless. As a first test Ω has been calculated for five selected systems at nine different intersystem distance values ranging from the equilibrium distance to twice as much (see supporting information). The five values of the overlap parameter and their spread approach zero at about twice the equilibrium distance. One can envision several choices of ρA (~r) and ρB (~r) for the evaluation of Ω depending on different variants and techniques of FDET. In the context of linearized FDET 29 (or conventional FDET) one has access to the initial density ρref r) and a stateA (~ specific (optimized) embedded density ρIA (~r). Beside the initial ρB (~r) also an optimized environment density is available through freeze-and-thaw cycles 44 (alternating the role of the embedded species). In our calculations we chose to compute Ω = Ωref with the reference density ρref r) at MP2/cc-pVDZ level of theory and the corresponding environment density, A (~ both of which are used to construct the embedding potential and are therefore available before the most computationally demanding portion of the calculation. It is worth noticing that Ωref is the same for each state for a given system.

2.4

Computational Details

Concerning the structures of the intermolecular complexes described in Sec. 2.1, the majority of structures were optimized with the TURBOMOLE 6.6 program 59 either at RI-CC2/ccpVTZ 60 (all complexes of 1, 4, and 5) or B3LYP/cc-pVTZ (all complexes of 2, 3, 7, 6, and 10) level of theory. The structures of 8a and 9a were optimized at MP2/cc-pVDZ level of theory using Q-Chem 4.3. 61 All molecular geometries used in this work (except for 8b and 12a) are listed in the supporting information.

14

ACS Paragon Plus Environment

Page 15 of 43 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

In each single-point calculation reported in the results section, ADC(2) for the complex, ADC(2) for the isolated chromophore, FDE-ADC(2) for the embedded chromophore, and Hartree-Fock or Kohn-Sham/PBE for the generation of ρB (~r), the cc-pVDZ basis set 60 was used. In order to isolate the electronic component of the solvatochromic shift, the geometry in all calculations was based on that of the intermolecular complex. For details regarding the embedding calculations combining Linearized FDET with ADC(2), we refer the reader to Refs. 31,32. Here, we specify only the most relevant features for the present analyses. All FDE-ADC(2) calculations in this work were carried out with a developer version of Q-Chem 5.0 employing the so-called supermolecular expansion (SE), i.e. all embedding steps utilize all available basis functions of the supersystem. Such SE calculations are computationally inefficient because the number of basis functions is the same as in conventional ADC(2) calculations for the whole system. However, the use of SE is motivated in the present analysis by the necessity to reduce the number of approximations. Note that if the same set of basis functions is used in both ADC(2) and FDE-ADC(2) calculations, the only sources of error are the violation of Eq. 7 locally for the chosen ρB (~r) and the approximation nad nad nad denotes an approximate explicit bi-functional. [ρA , ρB ](~r), where v˜xc,T [ρA , ρB ](~r) ≈ v˜xc,T vxc,T

Furthermore, we employed two different approximations for the functional vTnad [ρA , ρB ](~r): LDA type (Thomas-Fermi functional 62,63 ) and GGA type (GGA97, the bi-functional introduced in Refs. 24,25 which uses the enhancement factor of the Perdew-Wang91 (PW91) exchange functional 64,65 re-parametrized for the kinetic energy by Lembarki and Chermette 66 ). Throughout this work the PBE expression for the exchange-correlation functional 67 was used nad to approximate the non-additive exchange-correlation potential vxc [ρA , ρB ](~r). Unless ex-

plicitly specified, ρB (~r) used in FDET calculations was obtained from Hartree-Fock calculations for the isolated environment and the Thomas-Fermi functional was applied to evaluate vTnad [ρA , ρB ](~r). All plots were created using the matplotlib (v2.0.0) Python module. 68 The linear regression was carried out with the statsmodels (v0.8.0) Python module 69 and data processing

15

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

was facilitated using the pandas (v0.19.2) Python module. 70

3

Results and Discussion

The results section is organized as follows. Section 3.1 reports the value-by-value comparison between the corresponding excitation energies obtained from the supermolecular (ADC(2)) and embedded (FDET-ADC(2)) calculations. The dataset includes up to nine excited states for each considered complex. The analysis of the errors is made for all excitations in the dataset as well as for two separate subsets including excitations of different nature (nπ ∗ and ππ ∗ ). Moreover, correlations between the errors of the shifts and the magnitude of the shifts are investigated in detail. In the subsequent sections we investigate connections of the data obtained in Section 3.1 with other selected parameters. Section 3.2 concerns the relation between the error of the shift and the a priori available parameter Ωref (see Section 2.3 for details). Such analysis aims at determining the applicability domain of the used approximations defined by means of a well-defined system-independent parameter. Section 3.3 concerns the expected worsening of the performance of FDET for excitations nad [ρA , ρB ](~r)). Here, we take advantage extending to the environment (with approximated vxc,T

of the fact that in our dataset the information about the excited state wavefunction for each considered complex is available at the supermolecular ADC(2) level of theory. Of course, the information about the delocalization is not available in practical embedding calculations. The relation between the errors in the calculated shifts and the localization of the excited state is investigated in the present work not to introduce another predictive reliability threshold but rather to understand the origin of the errors or exceptions from the trends noted for other descriptors. Finally, Section 3.4 concerns another factor which can be expected to lead to errors nad in embedding calculations, i.e. the approximations used for the functional vxc,T [ρA , ρB ](~r).

16

ACS Paragon Plus Environment

Page 16 of 43

Page 17 of 43 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

nad Since the exact density functional vxc,T [ρA , ρB ](~r) is not available, the analysis is made for

pragmatic purposes, namely to determine the sensitivity of the shifts on the choice of two conventional approximations used in the literature for this functional.

3.1

Errors in Excitation Energies

All numerical values discussed and presented graphically in the Results section are available in the supporting information. Only in a few instances the numerical values of individual excitation energies are specified in the main text. The error in FDE-ADC(2) excitation energy is defined as the difference between the excitation energy obtained in FDE-ADC(2) and ADC(2) calculations for the complex: δ = FDE-ADC(2) − ADC(2) .

(15)

In most analyses, we discuss not the excitation energies but their shifts compared to the isolated chromophore. The ultimate goal of any embedding method is, after all, the ability to predict reliably the effect of the environment on the properties of the embedded species. Note that the errors in the excitation energies defined above and the errors in the shifts are the same. In Fig. 2, the errors (δ) are plotted against the shifts (∆), and same scale is used for both axes. The ideal case of a perfect embedding method would therefore correspond to the errors forming a horizontal line at δ = 0. The other extreme, in which the errors are scattered randomly, would show the embedding method to be unusable. If the calculated errors form a straight line, more information can be extracted from such plots and indicate systematic sources of errors. Note that, the slope of that line represents the fitted relative error. The mean absolute error (MAE) in excitation energy amounts to 39.2 meV if evaluated for all 351 computed excitations. The individual MAE for nπ ∗ and ππ ∗ states are 26.6 meV and 45.7 meV, respectively. It is immediately recognizable in Fig. 2, that the errors trend

17

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1.50

nπ ∗ state ππ ∗ state

y1 = -0.04*x + 0.02; R2 = 0.087 y2 = -0.18*x + 0.02; R2 = 0.856

1.25 1.00 0.75

δ [eV]

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 18 of 43

0.50 0.25 0.00 −0.25 −0.50 −0.75 −0.75 −0.50 −0.25

0.00

0.25

0.50

0.75

1.00

1.25

1.50

∆ [eV]

Figure 2: Correlation of the solvent shift ∆ at ADC(2)/cc-pVDZ level of theory and the error of the embedding method δ. The regression was performed separately on the set of nπ ∗ and ππ ∗ states, respectively. A red shift corresponds to ∆ < 0. differently for nπ ∗ and ππ ∗ states. The large majority of nπ ∗ solvent shifts are overestimated independent of the magnitude of the shift, i.e. the errors show a fluctuation even for smaller shifts. The magnitude of the solvent shift for ππ ∗ states, however, is consistently underestimated and follows a linear trend as indicated by the regression line. From the least squares fit, a fitted relative error of 18 % for ππ ∗ solvent shifts is obtained. The resulting coefficients 2 of determination (R2 ) of both types of excitations are very different. While Rππ ∗ is 86 % 2 explaining most of the variability around the mean, Rnπ ∗ is almost zero indicating that the

18

ACS Paragon Plus Environment

Page 19 of 43 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

fit function does not perform significantly better than the simple mean. The errors in nπ ∗ excitations effectively follow a Gaussian-like distribution (see supporting information). The largest red and blue shifts of nπ ∗ states for micro-solvated clusters are in systems 10a (-0.163 eV, S7 ) and 6c (1.484 eV, S3 ). For ππ ∗ excitations the systems exhibiting the largest shifts correspond to 1q (-0.710 eV, S6 ) and 6b (1.204 eV, S4 ). For nπ ∗ states the two largest errors in excitation energies stem from 1v (0.122 eV, S2 ) and 8b (0.105 eV, S5 ). For excited states of ππ ∗ character the two largest errors δ belong to 6c (-0.188 eV, S1 ) and 1q (0.172 eV, S6 ). An important characteristic of 1v is the almost perpendicular arrangement of the molecular planes of the chromophore and environment (](CB -OB -NA -CA ) = −79◦ ). Also, the carboxyl oxygen atom of formic acid is pointing in the opposite direction of 1 in this intermolecular complex. In comparison to the structurally similar systems 1t and 1x, in which both molecular planes are aligned, the error in the S2 excitation energy is notably increased (0.079 and 0.055 eV, respectively). Chromophore and environment in 8b exhibit a π-interaction in the ground state (6.8 kcal/mol 48 ). Upon excitation, such systems often entail the formation of very delocalized exciplex-like states, implying particle exchange between the two subsystems and therefore violation of Eq. 7. The S1 state in 6c shows a strong hypsochromic shift (∆ = 0.91 eV). Also, a substantial charge transfer character is to be expected in the excited state due to the betainic nature of 6. For the first four ππ ∗ excitations in 6c, the electron-hole distance |~re − ~rh | based on the analysis of the transition density matrix is in the range of 3.8 ˚ A to 4.7 ˚ A. System 1q features a negatively charged environment and the shortest hydrogen bond (1.358˚ A) in our benchmark set. The fact that nπ ∗ and ππ ∗ excitations display different trends in errors is intriguing and calls for an explanation. To this end we address this issue in a preliminary but intuitive way by means of difference density plots for a selection of excitations that we deem representative (see supporting information). For each system (1k, 1q, 8b, and 6c) one nπ ∗

19

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

and one ππ ∗ excitation were considered. The difference density plots were generated from FDE-ADC(2) calculations for the selected complexes and compared to the corresponding supermolecular reference. In three out of the four cases the difference densities for nπ ∗ states of the supermolecular reference exhibit considerable contributions on the environment molecules, whereas the corresponding difference densities for ππ ∗ states show small or no density changes at all in the same region. Embedding calculations fail to reproduce this feature of difference density for nπ ∗ states (see SI-Figures and discussion in Sec. 2.2). Moreover, the error in excitation energy for an nπ ∗ state is smaller than that of the corresponding ππ ∗ state in three out of the four complexes, which is in line with the mean absolute errors reported above. Clearly - as the S4 (nπ ∗ ) in system 1k (δ = 0.003 eV) illustrates - neglect of a substantial density change in the environment can only lead to small errors in nπ ∗ excitation energies if other factors compensate for the failure to properly describe density changes in the environment. This error cancellation mechanism seems to be absent or significantly less pronounced in the case of the selected ππ ∗ states giving rise to a more consistent (relative) error. Although it is tempting to accept this as an explanation for the different behaviour of nπ ∗ and ππ ∗ excitations in Fig. 2, one has to be cautious extrapolating from eight cases to the full dataset. So far, the discussion focused on extreme cases of FDE-ADC(2) errors and salient properties of these complexes. In the next sections, it is our intention to find possible correlations between the error and other descriptors. As elucidated in Section 1 and 2.2, we expect the accuracy of FDET-based methods to be affected by the density overlap of the two subsystems (measured by Ωref ) and the extent of delocalization of the (reference) excited state (measured by rw ).

3.2

Overlap parameter Ω

Fig. 3 shows the relation between Ωref and the error in FDE-ADC(2) excitation energies for nπ ∗ states and ππ ∗ states, respectively. The color coding allows to consider both the 20

ACS Paragon Plus Environment

Page 20 of 43

Page 21 of 43

absolute and (qualitatively) the relative error. Concerning the nπ ∗ excitations, Figures 2 nπ ∗ states

ππ ∗ states 1.4

0.175 1.2 0.150 1.0

0.125

|δ| [eV]

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

0.100

0.8

0.075

0.6

0.050 0.4 0.025 0.2 0.000 50

100

150

200

250

50

Ωref [ppm]

100

150

Ωref [ppm]

200

250

|∆| [eV]

Figure 3: Relation of the overlap parameter Ωref and the unsigned error |δ| of the embedding method for nπ ∗ states (left) and ππ ∗ states (right). Color coding is used additionally to indicate the magnitude of the shift |∆| for each excitation. and 3 indicate that the errors (which are nevertheless small) correlate with neither the interaction strength between the chromophore and the environment (measured by ∆) nor with the magnitude of the overlap Ωref . Even for systems with a similar magnitude of the overlap parameter, for instance Ωref ≈ 100 ppm, the errors vary over the whole range of |δ|. One surprising case is system 8b, for which the error is comparatively large in view of the small overlap (Ωref = 49 ppm). Consequently, this case shows the limits of using a parameter solely based on overlap. We attribute this error mainly to environment contributions to the excited state which have been quantified using the rw parameter (see Eq. 10). The three excitations in question here yield rw values between 0.81 and 0.86. In fact, all of the first eight excited states in the complex 8b show an involvement of the environment in the excitation, i.e. all reduced weight ratios rw are smaller than 0.88. For ππ ∗ states (Fig. 3 right), the errors in excitation energy show some correlation with Ωref , even for large overlap values (e.g. S1 and S6 of system 1q with Ωref = 245.9 ppm). One 21

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

also observes an increasing variance of errors with larger values for Ωref . Several excitations for system 6b (Ωref = 131.8 ppm, S1 , S2 , S4 , S5 , S7 ) and 6c (Ωref = 155.3 ppm, S1 , S2 , S8 ) clearly stand out of the general pattern. Nonetheless, their shifts are extremely large, as indicated by the color coding, so that also their corresponding relative errors are consistent with the fitted relative error (cf. Fig. 2). Notably, all systems displayed values of Ωref lower than 245.9 ppm. Larger values ought to be considered to fall within untested domain. The relation between overlap and error for ππ ∗ states suggests an increase of the likelihood of large errors for such excitations.

3.3

Reduced weight ratio

As discussed in Sec. 2.2, FDET does not account properly for charge transfer between the chromophore and its environment. Measuring the locality of a state is, therefore, crucial to investigate the error of this approximation. Fig. 4 shows the magnitude of the excitation energy error relative to the localization parameter rw . Surprisingly, the error does not seem to correlate directly with the reduced weight ratio. However, it is possible that the error due to excitations from and to the environment accounts only for a small part in the total error. Hence, a potential correlation between rw and |δ| could be hiding beneath a larger effect, e.g. increased overlap. Fig. 4 also reveals two different domains of localization. While the majority of ππ ∗ states are well localized, most nπ ∗ states are not. The two lowest occurences of rw for excited states with nπ ∗ character correspond to 1u (0.55, S4 ) and 1x (0.72, S4 ). For ππ ∗ states the two lowest rw values are observed for 8b (0.55, S2 ) and 1u (0.74, S7 ). Introducing the reduced weight ratio as another dimension in the comparison of |δ| and Ωref gives a better indication of its role in the total FDE-ADC error. For nπ ∗ states (see Fig. 5 left) there are now two trends visible: the excited states with the largest error show a high extent of delocalization while those well localized on the chromophore (rw ≈ 1) display the lowest errors. There are though some delocalized states with considerably low errors. Therefore a high value for rw and a low value for Ωref are a sufficient but not a necessary 22

ACS Paragon Plus Environment

Page 22 of 43

Page 23 of 43

0.175

nπ ∗ state ππ ∗ state

0.150

|δ| [eV]

0.125

0.100

0.075

0.050

0.025

0.000 0.6

0.7

0.8

0.9

1.0

rw

Figure 4: Relation of the reduced weight ratio and the error of the embedding method for all states. nπ ∗ states

ππ ∗ states

0.175

0.9

0.150

0.125 0.8

|δ| [eV]

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

0.100

0.075 0.7 0.050

0.025 0.6 0.000 50

100

150

200

250

50

Ωref [ppm]

100

150

Ωref [ppm]

200

250

rw

Figure 5: Absolute error |δ| versus Ω for nπ ∗ states (left) and ππ ∗ states (right). The individual excitations are color-coded by the reduced weight ratio rw . condition to obtain small errors in excitation energy shifts, which reflects the multivariant nature of the error. Interestingly, for system 1u the errors for the S2 and S4 states (Ωref = 47 ppm) are very similar although the respective values for rw differ by more than 0.4. However, the solvent shifts for these two states (0.04 eV and 0.07 eV, respectively) are in the same order of 23

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

Page 24 of 43

magnitude as the mean absolute error for this benchmark set. For the other few delocalized cases below |δ| = 0.02 eV, in particular excited states of 6a, 1h and 8a, the solvent shifts are one order of magnitude larger ranging from 0.18 eV to 0.70 eV. Fig. 5 also illustrates the relation of δ, Ωref and rw for excited states of ππ ∗ character. The lack of correlation between δ and rw (cf. Fig. 4) (at first unexpected) becomes apparent. Since the majority of ππ ∗ states are well localized on the chromophore, it is the value of the overlap parameter Ωref which relates to the error. A few data points cannot be correlated with Ωref and rw , for instance errors of 2g (S2 , Ωref = 28 ppm) and 11a (S1 , S7 , Ωref = 75 ppm) exhibit low values for Ωref and a high degree of localization rw > 0.98. This shows that some residual underlying factors to the errors have not been identified yet.

3.4

Other approximations for vTnad [ρA , ρB ](~r) and ρB (~r)

So far, our discussion has concerned only FDE-ADC(2) calculations employing the PBE nad exchange-correlation functional to approximate vxc [ρA , ρB ](~r), the Thomas-Fermi kinetic

energy functional to approximate vTnad [ρA , ρB ](~r), and ρB (~r) obtained from Hartree-Fock calculations. Each of these factors contributes to the errors and this section collects additional results concerning the contribution of these factors to the overall errors and their trends. To this end, the statistical analysis was repeated using other choices for v˜Tnad [ρA , ρB ](~r) and ρB , that is the GGA97 approximation for v˜Tnad [ρA , ρB ](~r) and Kohn-Sham/PBE calculations to generate ρB (~r). The mean absolute errors for each of the resulting four types of FDE-ADC(2) calculations are collected in Tab. 2 and respective plots in the same fashion as Fig. 2 are available in the supporting information. Neither the choice of method used to generate ρB (~r), nor the approximations for vTnad [ρA , ρB ](~r), have a noticable influence onthe overall quality of the FDE-ADC(2) results for ππ ∗ states. Even the fitted relative errors (see Tab. 3) remain almost unaffected by the applied ap24

ACS Paragon Plus Environment

Page 25 of 43 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

Table 2: Mean absolute errors in excitation energies for ππ ∗ , nπ ∗ states and both combined with respect to four different types of approximations in FDE-ADC(2). v˜Tnad [ρA , ρB ](~r) ρB (~r) LDA HF LDA PBE GGA97 HF GGA97 PBE

MAE() [meV] MAE(ππ∗ ) [meV] MAE(nπ∗ ) [meV] 39.2 45.7 26.6 41.9 48.6 28.7 43.8 44.9 41.7 41.8 48.1 29.7

proximations. A weak dependence on the choice of ρB (~r) is in line with general variational character of FDET and the trends reported in the literature. 32,42 Interestingly, for the nπ ∗ states, switching between the approximations listed in Tab. 2 has a more pronounced effect. Firstly, determining ρB (~r) at Hartree-Fock level of theory leads to a noticeable overall increase of the error if the GGA97 bi-functional is used to approximate vTnad [ρA , ρB ](~r). Secondly, utilisation of a PBE density instead of a Hartree-Fock density for ρB (~r) leads to a shift in the error δ towards negative values, while a change in the approximation to vTnad [ρA , ρB ](~r) from LDA to GGA97 has the opposite effect on the same order of magnitude. The errors are still small and acceptable for practical applications, but the origin of this increase and change of errors remains to be investigated. Table 3: Linear regression parameters δ = a·∆+b for four different types of approximations in FDE-ADC(2) calculations. v˜Tnad [ρA , ρB ](~r) ππ ∗ states LDA LDA GGA97 GGA97 nπ ∗ states LDA LDA GGA97 GGA97

R2

ρB (~r)

a

b [eV]

HF PBE HF PBE

−0.18 −0.20 −0.18 −0.20

0.02 0.02 0.02 0.02

0.856 0.843 0.856 0.850

HF PBE HF PBE

−0.04 −0.11 0.03 −0.05

0.02 0.02 0.03 0.02

0.084 0.537 0.029 0.134

25

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

4

Suggested benchmark set

The set of intermolecular complexes introduced in Sec. 2.1 is large and can be reduced in order to suggest a benchmark set for future statistical analyses not only for embedding methods, but also for solvatochromic shifts of such complexes. Among these studied complexes there are some redundant cases of chemically similar interactions between chromophore and environment leading to almost identical solvatochromic shifts at ADC(2) level of theory. Such redundancies, as well as the two examples exhibiting π-interactions (8b and 12a) were eliminated from the proposed benchmark set. The resulting set, labeled XH-27 (see supporting information), spans 27 individual systems based on 9 chromophores and a total of 121 excitations (if no more than the first five excited states per complex are considered). The mean absolute error of FDE-ADC(2) excitation energies of this subset amounts to 45.8 meV.

5

Conclusions

A comprehensive analysis of the errors of 351 electronic excitation energies in a representative set of 52 embedded clusters shows good performance of the used FDET-based method. It can thus be considered as an attractive alternative to supermolecular calculations. The mean absolute error of the excitation energy (and its environment-induced shift) amounts to only 0.039 eV, whereas the maximum absolute error is 0.188 eV. It occurs for pyridiniumyl benzimidazolide · 2 HCOOH (6c), for which the shift is equal to 0.907 eV (S1 ). The largest absolute errors are observed in systems where the magnitude of the excitation energy shift is significantly larger (0.6 to 1.5 eV) than the error itself. The errors of the excitation energies were analyzed further by investigating their correlation with the magnitude of the shift, the magnitude of the density overlap, and the quantity measuring the delocalization of the excited state. For nπ ∗ excitations, the errors and the magnitudes of the shifts are not correlated in the whole range of variation of the shift which stretches from -0.16 to 1.48 eV. On the other hand, for ππ ∗ excitations (231 out of the 351 26

ACS Paragon Plus Environment

Page 26 of 43

Page 27 of 43 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

total excitations), the calculated shifts are systematically underestimated by about 18 %. No correlation with the errors in the excitation energy for nπ ∗ states with the magnitude of the overlap (Ωref ) is observed. However, the errors weakly correlate, with Ωref for ππ ∗ states, i.e. the larger the overlap the larger the error in the excitation energy. This suggests the usage of Ωref as a reliability/confidence parameter which can be applied a priori in FDET-based calculations. Values of Ωref < 250 ppm characterize all excitations considered. The correlation of Ωref with the error of the excitation energy for ππ ∗ excitations suggests that intermolecular complexes for which Ωref > 250 ppm lie outside of the domain of applicability of the approximations used for the FDET embedding potential. We also emphasize that other choices for the measure of the overlap magnitude can be envisaged. The one used in this work was chosen for the sake of convenience and expected transferability. For charge transfer excitations (from/to the chromophore to/from the environment) as well as for the excitations leading to the change of electron densities delocalized over the whole complex, FDET-based methods are expected to fail. Such cases should not be considered targets of FDET-based calculations at the current stage of development of approximations for non-additive functionals. Indeed, a high degree of localization and low overlap, measured by rw and Ωref , are observed to be a sufficient, although not necessary, condition for small errors in excitation energy. In contrast to Ωref , which can be used a priori, rw is available neither a priori nor a posteriori in embedding methods. The embedding potential of the form given in Eq. 3 is used in several variants of FDET differing in the treatment of the quantum many-body problem. However, all of these variants feature the same form of the embedding potential as a functional of charge densities. The accuracy of the approximations used for the embedding potential determines the shift of any observable due to the interactions with the environment. For this reason, we believe that the conclusions concerning the factors determining the errors of the excitation energy shifts based on comparisons between embedding (FDE-ADC(2)) and supermolecular (ADC(2)) calculations are transferable for other possible methods based on FDET.

27

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

The set of complexes and the localized excitation energies is rather large and includes several systems of similar character. A smaller set proposed in the present work labeled as XH-27 can be used for benchmarking purposes for both supermolecular and embedding methods. It includes a representative set of excitations as well as most common arrangement of interacting molecules.

Supporting Information Available A zip archive is available for download which comprises the following: ˆ all structures (folder): contains all structures used in this work in .xyz format

except for molecular geometries taken from Ref. 48, ˆ XH-27 (folder): contains all structures from the chosen subset denoted as XH-27 in

.xyz format and ˆ supporting info.pdf (file): Supporting information document which includes addi-

tional plots and the benchmarking results listed in tabular form. This material is available free of charge via the Internet at http://pubs.acs.org/.

Acknowledgement This research was supported by grant from Swiss National Science Foundation (Grant No. 200021-152779). AD and SP acknowledge support from the Heidelberg Graduate School ”Mathematical and Computational Methods for the Sciences“ (GSC 220).

References (1) Wesolowski, T. A.; Warshel, A. Frozen density functional approach for ab initio calculations of solvated molecules. J. Phys. Chem. 1993, 97, 8050–8053.

28

ACS Paragon Plus Environment

Page 28 of 43

Page 29 of 43 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

(2) Wesolowski, T. A. Embedding a multideterminantal wave function in an orbital-free environment. Phys. Rev. A 2008, 77, 012504. (3) Pernal, K.; Wesolowski, T. A. Orbital-free effective embedding potential: Densitymatrix functional theory case. Int. J. Quantum Chem. 2009, 109, 2520–2525. (4) Govind, N.; Wang, Y.; da Silva, A.; Carter, E. Accurate ab initio energetics of extended systems via explicit correlation embedded in a density functional environment. Chem. Phys. Lett. 1998, 295, 129–134. (5) Stefanovich, E. V.; Truong, T. N. Embedded density functional approach for calculations of adsorption on ionic crystals. J. Chem. Phys. 1996, 104, 2946–2955. (6) Tamukong, P. K.; Khait, Y. G.; Hoffmann, M. R. Density Differences in Embedding Theory with External Orbital Orthogonality. J. Phys. Chem. A 2014, 118, 9182–9200. (7) Chulhai, D. V.; Jensen, L. Frozen Density Embedding with External Orthogonality in Delocalized Covalent Systems. J. Chem. Theory Comput. 2015, 11, 3080–3088. (8) Wesolowski, T. A.; Savin, A. In Recent Prog. Orbital-Free Density Funct. Theory; Wesolowski, T. A., Wang, Y. A., Eds.; Recent Progress in Computational Chemistry; World Scientific: Singapore, 2013; Vol. 6; pp 277–295. (9) Unsleber, J. P.; Neugebauer, J.; Jacob, C. R. No need for external orthogonality in subsystem density-functional theory. Phys. Chem. Chem. Phys. 2016, 18, 21001–21009. (10) Fux, S.; Kiewisch, K.; Jacob, C. R.; Neugebauer, J.; Reiher, M. Analysis of electron density distributions from subsystem density functional theory applied to coordination bonds. Chem. Phys. Lett. 2008, 461, 353–359. (11) G¨otz, A. W.; Beyhan, S. M.; Visscher, L. Performance of Kinetic Energy Functionals for Interaction Energies in a Subsystem Formulation of Density Functional Theory. J. Chem. Theory Comput. 2009, 5, 3161–3174. 29

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

(12) Fux, S.; Jacob, C. R.; Neugebauer, J.; Visscher, L.; Reiher, M. Accurate frozen-density embedding potentials as a first step towards a subsystem description of covalent bonds. J. Chem. Phys. 2010, 132, 164101. (13) Jacob, C. R.; Visscher, L. A subsystem density-functional theory approach for the quantum chemical treatment of proteins. J. Chem. Phys. 2008, 128, 155102. (14) Zhou, B.; Alexander Wang, Y.; Carter, E. A. Transferable local pseudopotentials derived via inversion of the Kohn-Sham equations in a bulk environment. Phys. Rev. B 2004, 69, 125109. (15) Roncero, O.; de Lara-Castells, M. P.; Villarreal, P.; Flores, F.; Ortega, J.; Paniagua, M.; Aguado, A. An inversion technique for the calculation of embedding potentials. J. Chem. Phys. 2008, 129, 184104. (16) Goodpaster, J. D.; Ananth, N.; Manby, F. R.; Miller, T. F. Exact nonadditive kinetic potentials for embedded density functional theory. J. Chem. Phys. 2010, 133, 084103. (17) Banafsheh, M.; Adam Wesolowski, T. Nonadditive kinetic potentials from inverted Kohn-Sham problem. Int. J. Quantum Chem. 2018, 118, e25410. (18) Wesolowski, T. A. In Comput. Chem. Rev. Curr. Trends; Leszczynski, J., Ed.; Computational Chemistry: Reviews of Current Trends; World Scientific, 2006; Vol. 10; pp 1–82. (19) Libisch, F.; Huang, C.; Carter, E. A. Embedded Correlated Wavefunction Schemes: Theory and Applications. Acc. Chem. Res. 2014, 47, 2768–2775. (20) Jacob, C. R.; Neugebauer, J. Subsystem density-functional theory. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2014, 4, 325–362. (21) Wesolowski, T. A.; Shedge, S.; Zhou, X. Frozen-Density Embedding Strategy for Multilevel Simulations of Electronic Structure. Chem. Rev. 2015, 115, 5891–5928. 30

ACS Paragon Plus Environment

Page 30 of 43

Page 31 of 43 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

(22) Krishtal, A.; Sinha, D.; Genova, A.; Pavanello, M. Subsystem density-functional theory as an effective tool for modeling ground and excited states, their dynamics and manybody interactions. J. Phys. Condens. Matter 2015, 27 . (23) Wasserman, A.; Nafziger, J.; Jiang, K.; Kim, M.-C.; Sim, E.; Burke, K. The Importance of Being Inconsistent. Annu. Rev. Phys. Chem. 2017, 68, 555–581. (24) Wesolowski, T. A.; Chermette, H.; Weber, J. Accuracy of approximate kinetic energy functionals in the model of KohnSham equations with constrained electron density: The FHNCH complex as a test case. J. Chem. Phys. 1996, 105, 9182–9190. (25) Wesolowski, T. A. Density functional theory with approximate kinetic energy functionals applied to hydrogen bonds. J. Chem. Phys. 1997, 106, 8516–8526. (26) Bernard, Y. A.; Dulak, M.; Kami´ nski, J. W.; Wesolowski, T. A. The energy-differences based exact criterion for testing approximations to the functional for the kinetic energy of non-interacting electrons. J. Phys. A Math. Theor. 2008, 41, 055302. (27) McRae, E. G. Theory of Solvent Effects on Molecular Electronic Spectra. Frequency Shifts. J. Phys. Chem. 1957, 61, 562–572. (28) Wesolowski, T. A. Embedding potentials for excited states of embedded species. J. Chem. Phys. 2014, 140, 18A530. (29) Zech, A.; Aquilante, F.; Wesolowski, T. A. Orthogonality of embedded wave functions for different states in frozen-density embedding theory. J. Chem. Phys. 2015, 143, 164106. (30) Dreuw, A.; Wormit, M. The algebraic diagrammatic construction scheme for the polarization propagator for the calculation of excited states. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2015, 5, 82–95.

31

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

(31) Prager, S.; Zech, A.; Aquilante, F.; Dreuw, A.; Wesolowski, T. A. First time combination of frozen density embedding theory with the algebraic diagrammatic construction scheme for the polarization propagator of second order. J. Chem. Phys. 2016, 144, 204103. (32) Prager, S.; Zech, A.; Wesolowski, T. A.; Dreuw, A. Implementation and Application of the Frozen Density Embedding Theory with the Algebraic Diagrammatic Construction Scheme for the Polarization Propagator up to Third Order. J. Chem. Theory Comput. 2017, 13, 4711–4725. (33) Kl¨ uner, T.; Govind, N.; Wang, Y.; Carter, E. Prediction of Electronic Excited States of Adsorbates on Metal Surfaces from First Principles. Phys. Rev. Lett. 2001, 86, 5954– 5957. (34) Aquilante, F.; Wesolowski, T. A. Self-consistency in frozen-density embedding theory based calculations. J. Chem. Phys. 2011, 135, 084120. (35) Libisch, F.; Cheng, J.; Carter, E. A. Electron-Transfer-Induced Dissociation of H 2 on Gold Nanoparticles: Excited-State Potential Energy Surfaces via Embedded Correlated Wavefunction Theory. Zeitschrift f¨ ur Phys. Chemie 2013, 227, 1455–1466. (36) H¨ofener, S.; Gomes, A. S. P.; Visscher, L. Solvatochromic shifts from coupled-cluster theory embedded in density functional theory. J. Chem. Phys. 2013, 139, 104106. (37) Daday, C.; K¨onig, C.; Neugebauer, J.; Filippi, C. Wavefunction in Density Functional Theory Embedding for Excited States: Which Wavefunctions, which Densities? ChemPhysChem 2014, 15, 3205–3217. (38) Dresselhaus, T.; Neugebauer, J.; Knecht, S.; Keller, S.; Ma, Y.; Reiher, M. Selfconsistent embedding of density-matrix renormalization group wavefunctions in a density functional environment. J. Chem. Phys. 2015, 142, 044111.

32

ACS Paragon Plus Environment

Page 32 of 43

Page 33 of 43 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

(39) Wesolowski, T. A. Hydrogen-Bonding-Induced Shifts of the Excitation Energies in Nucleic Acid Bases: An Interplay between Electrostatic and Electron Density Overlap Effects. J. Am. Chem. Soc. 2004, 126, 11444–11445. (40) Casida, M. E.; Wesolowski, T. A. Generalization of the Kohn-Sham equations with constrained electron density formalism and its time-dependent response theory formulation. Int. J. Quantum Chem. 2004, 96, 577–588. (41) Neugebauer, J. Couplings between electronic transitions in a subsystem formulation of time-dependent density functional theory. J. Chem. Phys. 2007, 126, 134116. (42) Humbert-Droz, M.; Zhou, X.; Shedge, S. V.; Wesolowski, T. A. How to choose the frozen density in Frozen-Density Embedding Theory-based numerical simulations of local excitations? Theor. Chem. Acc. 2014, 133, 1405. (43) Cortona, P. Self-consistently determined properties of solids without band-structure calculations. Phys. Rev. B 1991, 44, 8454–8458. (44) Wesolowski, T. A.; Weber, J. Kohn-Sham equations with constrained electron density: an iterative evaluation of the ground-state electron density of interacting molecules. Chem. Phys. Lett. 1996, 248, 71–76. (45) Bach,

A.;

Leutwyler,

S.

Water-chain

clusters:

vibronic

spectra

of

7-

hydroxyquinoline·(H2O)n, n=14. Chem. Phys. Lett. 1999, 299, 381–388. (46) Thut, M.; Tanner, C.; Steinlin, A.; Leutwyler, S. Time-Dependent Density Functional Theory As a Tool for Isomer Assignments of Hydrogen-Bonded Solute·Solvent Clusters. J. Phys. Chem. A 2008, 112, 5566–5572. (47) Thut, M.; Manca, C.; Tanner, C.; Leutwyler, S. Spectral tuning by switching CHO hydrogen bonds: Rotation-induced spectral shifts of 7-hydroxyquinolineHCOOH isomers. J. Chem. Phys. 2008, 128, 024304. 33

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

(48) ez´aˇc, J.; Riley, K. E.; Hobza, P. S66: A Well-balanced Database of Benchmark Interaction Energies Relevant to Biomolecular Structures. J. Chem. Theory Comput. 2011, 7, 2427–2438. (49) Plasser, F.; Wormit, M.; Dreuw, A. New tools for the systematic analysis and visualization of electronic excitations. I. Formalism. J. Chem. Phys. 2014, 141, 024106. (50) B¨appler, S. A.; Plasser, F.; Wormit, M.; Dreuw, A. Exciton analysis of many-body wave functions: Bridging the gap between the quasiparticle and molecular orbital pictures. Phys. Rev. A 2014, 90, 052521. (51) Plasser, F.; Thomitzni, B.; B¨appler, S. A.; Wenzel, J.; Rehn, D. R.; Wormit, M.; Dreuw, A. Statistical analysis of electronic excitation processes: Spatial location, compactness, charge transfer, and electron-hole correlation. J. Comput. Chem. 2015, 36, 1609–1620. (52) Luzanov, A. V.; Sukhorukov, A. A.; Umanskii, V. E. Application of transition density matrix for analysis of excited states. Theor. Exp. Chem. 1976, 10, 354–361. (53) Martin, R. L. Natural transition orbitals. J. Chem. Phys. 2003, 118, 4775–4777. (54) Lastra, J. M. G.; Kaminski, J. W.; Wesolowski, T. A. Orbital-free effective embedding potential at nuclear cusps. J. Chem. Phys. 2008, 129, 074107. (55) Tretiak, S.; Mukamel, S. Density Matrix Analysis and Simulation of Electronic Excitations in Conjugated and Aggregated Molecules. Chem. Rev. 2002, 102, 3171–3212. (56) Luzanov, A. V.; Prezhdo, O. V. Irreducible charge density matrices for analysis of many-electron wave functions. Int. J. Quantum Chem. 2005, 102, 582–601. (57) Luzanov, A. V.; Zhikol, O. A. Electron invariants and excited state structural analysis for electronic transitions within CIS, RPA, and TDDFT models. Int. J. Quantum Chem. 2010, 110, 902–924. 34

ACS Paragon Plus Environment

Page 34 of 43

Page 35 of 43 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

(58) Plasser, F.; Lischka, H. Analysis of Excitonic and Charge Transfer Interactions from Quantum Chemical Calculations. J. Chem. Theory Comput. 2012, 8, 2777–2789. (59) TURBOMOLE V6.6 2014,

a development of University of Karlsruhe and

Forschungszentrum Karlsruhe GmbH, 1989-2007, TURBOMOLE GmbH, since 2007; available from http://www.turbomole.com. (60) Dunning, T. H. Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. J. Chem. Phys. 1989, 90, 1007–1023. (61) Shao, Y.; Gan, Z.; Epifanovsky, E.; Gilbert, A. T.; Wormit, M.; Kussmann, J.; Lange, A. W.; Behn, A.; Deng, J.; Feng, X.; Ghosh, D.; Goldey, M.; Horn, P. R.; Jacobson, L. D.; Kaliman, I.; Khaliullin, R. Z.; Ku´s, T.; Landau, A.; Liu, J.; Proynov, E. I.; Rhee, Y. M.; Richard, R. M.; Rohrdanz, M. A.; Steele, R. P.; Sundstrom, E. J.; Woodcock, H. L.; Zimmerman, P. M.; Zuev, D.; Albrecht, B.; Alguire, E.; Austin, B.; Beran, G. J. O.; Bernard, Y. A.; Berquist, E.; Brandhorst, K.; Bravaya, K. B.; Brown, S. T.; Casanova, D.; Chang, C.-M.; Chen, Y.; Chien, S. H.; Closser, K. D.; Crittenden, D. L.; Diedenhofen, M.; DiStasio, R. A.; Do, H.; Dutoi, A. D.; Edgar, R. G.; Fatehi, S.; Fusti-Molnar, L.; Ghysels, A.; Golubeva-Zadorozhnaya, A.; Gomes, J.; Hanson-Heine, M. W.; Harbach, P. H.; Hauser, A. W.; Hohenstein, E. G.; Holden, Z. C.; Jagau, T.-C.; Ji, H.; Kaduk, B.; Khistyaev, K.; Kim, J.; Kim, J.; King, R. A.; Klunzinger, P.; Kosenkov, D.; Kowalczyk, T.; Krauter, C. M.; Lao, K. U.; Laurent, A. D.; Lawler, K. V.; Levchenko, S. V.; Lin, C. Y.; Liu, F.; Livshits, E.; Lochan, R. C.; Luenser, A.; Manohar, P.; Manzer, S. F.; Mao, S.-P.; Mardirossian, N.; Marenich, A. V.; Maurer, S. A.; Mayhall, N. J.; Neuscamman, E.; Oana, C. M.; Olivares-Amaya, R.; O’Neill, D. P.; Parkhill, J. A.; Perrine, T. M.; Peverati, R.; Prociuk, A.; Rehn, D. R.; Rosta, E.; Russ, N. J.; Sharada, S. M.; Sharma, S.; Small, D. W.; Sodt, A.; Stein, T.; St¨ uck, D.; Su, Y.-C.; Thom, A. J.; Tsuchimochi, T.; Vanovschi, V.; Vogt, L.; Vydrov, O.; Wang, T.; Watson, M. A.; Wenzel, J.; White, A.; Williams, C. F.; Yang, J.; 35

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

Yeganeh, S.; Yost, S. R.; You, Z.-Q.; Zhang, I. Y.; Zhang, X.; Zhao, Y.; Brooks, B. R.; Chan, G. K.; Chipman, D. M.; Cramer, C. J.; Goddard, W. A.; Gordon, M. S.; Hehre, W. J.; Klamt, A.; Schaefer, H. F.; Schmidt, M. W.; Sherrill, C. D.; Truhlar, D. G.; Warshel, A.; Xu, X.; Aspuru-Guzik, A.; Baer, R.; Bell, A. T.; Besley, N. A.; Chai, J.-D.; Dreuw, A.; Dunietz, B. D.; Furlani, T. R.; Gwaltney, S. R.; Hsu, C.-P.; Jung, Y.; Kong, J.; Lambrecht, D. S.; Liang, W.; Ochsenfeld, C.; Rassolov, V. A.; Slipchenko, L. V.; Subotnik, J. E.; Van Voorhis, T.; Herbert, J. M.; Krylov, A. I.; Gill, P. M.; Head-Gordon, M. Advances in molecular quantum chemistry contained in the Q-Chem 4 program package. Mol. Phys. 2015, 113, 184–215. (62) Thomas, L. H. The calculation of atomic fields. Math. Proc. Cambridge Philos. Soc. 1927, 23, 542. (63) Fermi, E. Eine statistische Methode zur Bestimmung einiger Eigenschaften des Atoms und ihre Anwendung auf die Theorie des periodischen Systems der Elemente. Zeitschrift f¨ ur Phys. 1928, 48, 73–79. (64) Perdew, J. P. In Electron. Struct. Solids; Ziesche, P., Eschrig, H., Eds.; Akademie Verlag: Berlin, 1991; pp 11–20. (65) Perdew, J. P.; Wang, Y. Accurate and simple analytic representation of the electron-gas correlation energy. Phys. Rev. B 1992, 45, 13244–13249. (66) Lembarki, A.; Chermette, H. Obtaining a gradient-corrected kinetic-energy functional from the Perdew-Wang exchange functional. Phys. Rev. A 1994, 50, 5328–5331. (67) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865–3868. (68) Hunter, J. D. Matplotlib: A 2D Graphics Environment. Comput. Sci. Eng. 2007, 9, 90–95.

36

ACS Paragon Plus Environment

Page 36 of 43

Page 37 of 43 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

(69) Seabold, S.; Perktold, J. 9th Python Sci. Conf.; 2010. (70) McKinney, W. Data Structures for Statistical Computing in Python. Proc. 9th Python Sci. Conf. 2010; pp 51–56.

37

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

Graphical TOC Entry

Schematic representation of the overlap region of ρA (~ r) and ρB (~ r).

38

ACS Paragon Plus Environment

Page 38 of 43

Page439 of 43 3

1 2 2 N 3 1 4 5 6 7 8 9 10 11 2 12 3 13 14 15 4 16 17 5 18 19 20 21 22 23 24 25 2 3 26 27 28 4 29 30 5 31 32 33

O Journal of Chemical Theory and Computation H 5 N 6 7 HN 1

5 6 7

8

O

1

3

4

N H

H2N

N

N 2

N

N

-

10

R

R=H R=Me

N

4 5

H

8

7

H N 1

5

4

2 3

1

H

5

3 4

6

9

6 7

2

O O 1 N ACS Paragon Plus Environment

10

O

1N

6

6

2

6

6

O

O

NH

5

1

8

3

4

NH

HN

9

7

6

3

4 5

O

2

N+

N9

9

O 11

4

3

2 3=(2)-

1

7

N

5

8

2

OH

8

6 1

8

11

NH2

12

1.50

δ [eV]

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

Journal of Chemical Theory and Computation

Page 40 of 43

nπ ∗ state ππ ∗ state

y1 = -0.04*x + 0.02; R2 = 0.087 y2 = -0.18*x + 0.02; R2 = 0.856

1.25 1.00 0.75 0.50 0.25 0.00 −0.25 −0.50 −0.75 −0.75 −0.50 −0.25

0.00

0.25

0.50

ACS Paragon Plus Environment

∆ [eV]

0.75

1.00

1.25

1.50

nπ ∗ states

Page 41 of 43

1.4 0.175 1.2 0.150 1.0

0.125

|δ| [eV]

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

ππ ∗ states

Journal of Chemical Theory and Computation

0.100

0.8

0.075

0.6

0.050 0.4 0.025 0.2 0.000 50

100

150

Ωref [ppm]

200

250

ACS Paragon Plus Environment

50

100

150

Ωref [ppm]

200

250

|∆| [eV]

0.175

Journal of Chemical Theory and Computation

Page 42 of 43

nπ state ππ ∗ state

0.150

0.125

|δ| [eV]

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



0.100

0.075

0.050

0.025

0.000 0.6

0.7

0.8

ACS Paragon Plus Environment

rw

0.9

1.0

nπ ∗ states

Page 43 of 43

0.175

0.9

0.150

0.125 0.8

|δ| [eV]

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

ππ ∗ states

Journal of Chemical Theory and Computation

0.100

0.075 0.7 0.050

0.025 0.6 0.000 50

100

150

Ωref [ppm]

200

250

ACS Paragon Plus Environment

50

100

150

Ωref [ppm]

200

250

rw