Development and Implementation of Excited-State Gradients for Local

Subscriber access provided by Columbia University Libraries

Spectroscopy and Excited States

Development and Implementation of ExcitedState Gradients for Local Hybrid Functionals Robin Grotjahn, Filipp Furche, and Martin Kaupp J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.9b00659 • Publication Date (Web): 29 Aug 2019 Downloaded from pubs.acs.org on August 31, 2019

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

Development and Implementation of Excited-State Gradients for Local Hybrid Functionals Robin Grotjahn, Filipp Furche, and Martin Kaupp∗, Technische Universit¨ at Berlin, Institut f¨ ur Chemie, Theoretische Chemie/Quantenchemie, Sekr. C7, Straße des 17.

Juni 135, D-10623, Berlin, Germany University of California, Irvine, Department of Chemistry, 1102 Natural Sciences II, Irvine, California

92697-2025, United States E-mail: [email protected]

Abstract

strongly on the choice of the XC functional. In particular, benchmark studies 19–22 show a dependence of the mean absolute errors (MAEs) of AEEs on the amount of exactexchange (EXX) admixture in hybrid functionals. However, a constant admixture of EXX delivers mixed accuracy for states of different character, such as charge-transfer, Rydberg and valence states. For instance, the popular B3LYP functional gives fairly accurate results for vertical valence excitations with MAEs of 0.22 eV 23 and for AEEs of lowlying states (0.19 eV 24 ) but fails for Rydberg (1.11–eV 23 ) and charge-transfer (1.36 eV 23 ) excitations. While chargetransfer excitation energies are usually improved by increased EXX admixture, 20,25 valence excitation energies are deteriorated. 20,26 This motivates the consideration of new classes of hybrid functionals for the calculation of AEEs. Range-separated hybrid functionals have been successful for excitations with charge-transfer character. 27 Also, recent work by Truhlar and coworkers 28 showed that after proper reoptimization their applicability is extended to include valence and Rydberg excitations. Local hybrid functionals 29,30 (local hybrids, LHs) are another logical extension of the concept of hybrid functionals. LHs use a local mixing function (LMF) governing the real-space position-dependent admixture of EXX. This allows for adaptive compensation of the self-interaction error in different electronic regions of the molecule. In spite of the simplicity of their basic idea and the still relatively low level of empiricism of current LHs, recent self-consistent implementations for ground-state SCF 31 and gradients 32 as well as vertical excitation energies 33 (VEEs) have demonstrated competitive performance. In particular, valence triplet VEEs are remarkably improved with MAEs of 0.16 eV (B3LYP: 0.45 eV) for the Thiel test set, 26,34 as well as for a more demanding set by Michl and coworkers 35 with MAEs of 0.50 eV (B3LYP: 0.74 eV). 36 For the calculation of AEEs it is, of course, possible to optimize the excited state with an already implemented (hybrid) functional and then calculate single-point LH energies using the existing TDDFT implementation. In fact, such a two-pronged approach is generally well-established and errors are usually under control, since AEEs depend quadratically on errors in the underlying structure. 19,37 However, VEEs depend linearly on the structures, and accurate excited-state structures are thus important for the calculation of fluorescence spectra. Moreover, analytical derivatives of excited-state energy surfaces are crucial ingredients for nonadiabatic molecular dynamics simulations. 38

Local hybrid functionals are a relatively recent class of exchange-correlation functionals that use a real-space dependent admixture of exact exchange. Here, we present the first implementation of TDDFT excited-state gradients for these functionals. Based on the ansatz of a fully variational auxiliary Lagrangian of the excitation energy, the working-equations for the case of a local hybrid functional are deduced. For the implementation, we derive the third-order functional derivatives used in the hyper-kernel and kernel-gradients following a semi-numerical integration scheme. The first assessment for a test set of small molecules reveals competitive performance for excited-state structural parameters with typical MAEs of 1.2 pm (PBE0: 1.4 pm) for bond lengths as well as for vibrational frequencies with typical MAEs of 81 cm−1 (PBE0: 76 cm−1 ). Excellent performance was found for adiabatic triplet excitation energies with typical MAEs of 0.08 eV (PBE0: 0.32–eV). In a detailed case analysis of the first singlet and triplet excited states of formaldehyde the conceptional (dis-)advantages of the local hybrid scheme for excited-state gradients are exposed.

1

Introduction

Excited-state potential-energy surfaces and adiabatic excitation energies (AEEs) are critical for the study of a wide range of photochemical and photophysical processes such as phosphorescence, 1 photoinduced reactions 2 , emission spectroscopy, 3 or singlet fission. 4 Time-dependent density functional theory, (TDDFT) within the adiabatic approximation, 5 i.e., the static limit of the frequency-dependent exchange-correlation (XC) kernel, enables quantum chemists to routinely calculate excited-state properties at affordable costs. Excited-state TDDFT gradients, required for efficient excited-state optimizations, were first implemented by Van Caillie et al. in the CADPAC package. 6,7 Here, we follow closely the implementation by Furche et al. 8 in the TURBOMOLE program package. 9 Further approaches have been implemented in various quantum chemistry packages and have been extended, e.g., using the Tamm-Dancoff approximation, 10,11 frequency domain approaches, 12 seminumerical integration techniques, 13 polarizable continuum solvent models, 14 QM/MM techniques, 15,16 long-range corrected functionals, 17 and spin-flip TDDFT. 18 In general, the accuracy of such calculations depends

ACS Paragon Plus Environment 1

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

LHs have been found to be well suited to treat the ground states of mixed-valence (MV) systems, including their VEEs, where a proper balance between reduced self-interaction errors and the simulation of left-right correlation is crucial. 39–43 It is expected that this usefulness will also extend to excited-state electronic structures featuring a similarly sensitive balance as ground-state MV systems. For the calculation of excited-state vibrational frequencies, it is essential to locate stationary points on the excitedstate surface, so that the harmonic approximation can be exploited. Excited-state vibrational frequencies and related quantities like vibronic couplings are vital in many applications 44 and also tend to be influenced significantly by the EXX admixture to the functional. 45,46 Hence, an excitedstate gradient implementation of LHs is needed if we want to make full use of the potential of this class of functionals in photochemistry and -physics.

γσσ (r) is the squared density gradient γσσ0 (r) = ∇ρσ (r) · ∇ρσ0 (r) ,

2.1

2.2

(7)

Excited-State Gradients

We will closely follow the exposition of Ref. 8 on LDA, GGA, TDHF, and GH TDDFT excited-state gradients in deriving the corresponding LH excited-state gradients, and we will focus particularly on the modifications needed to treat LHs. We start from the variational formulation of TDDFT, 58 i.e. excitation energies Ω are stationary points of the functional 1 X [(X + Y )iaσ (A + B)iaσjbσ0 (X + Y )jbσ0 2 0 iaσjbσ

+(X − Y )iaσ (A − B)iaσjbσ0 (X − Y )jbσ0 ] " # X (X + Y )iaσ (X − Y )iaσ , +Ω

Local Hybrid Functionals

In contrast to the constant EXX admixture of global hybrid functionals (global hybrids, GHs), LHs mix Hartree-Focklike exact-exchange energy densities with semi-local ones in a position-dependent way governed by the LMF aσ (r), leading to the expression XZ XZ LH 3 3 Exc = aσ (r) · eex [1 − aσ (r)] · esl x,σ (r)d r + x,σ (r)d r σ

+

1X σ Dpq ∇φpσ (r) · ∇φqσ (r) . 2 pq

τσ (r) =

Theory

3 esl c (r)d r ,

where (X + Y ) is the symmetric and (X − Y ) the antisymmetric transition-density matrix, and (A±B) are the orbital rotation Hessians 0

0

σσ (A − B)iaσjbσ0 = δij δab δσσ0 (σa − σi ) + f−iajb ,

0

(10)

where the kernel integrals are conveniently written as a second derivative w.r.t. density matrices (cf. Sec. S1) (3) 0

σσ fpqrs =

σ Dpq

is a matrix element of the ground-state density matrix in the MO basis. We will use a, b, c, . . . for virtual, i, j, k, . . . for occupied, and p, q, r, . . . for any (arbitrary) orbitals, while σ, σ 0 , σ 00 , . . . are spin labels. Using a constant for the LMF aσ (r) gets us back to the corresponding equations for a GH. While a number of LMFs have been suggested in the literature, 29,30,47–57 and the development of improved LMFs is an active area of research, here we focus exclusively on the simplest type of t-LMF for illustration. 29,48 The t-LMF is defined as the ratio of the von-Weizs¨ acker kinetic energy (r) and the conventional Pauli kineticdensity τσW (r) = γ8σσ ρ(r) energy density τσ (r), scaled by an adjustable prefactor b, γσσ (r) , 8 ρσ (r)τσ (r)

0

0

σσ σσ σσ = fpqrs ± fpqsr , f±pqrs

expressed in terms of Kohn-Sham MOs {φpσ (r)}:

aσ (r) = b ·

(9b)

and σp is the p-th eigenvalue corresponding to MO φpσ (r), RR σσ0 σσ 0 = wiajb (r, r0 )d3 rd3 r0 is the Coulombic part. For Jiajb σσ 0 we adopt the the matrix elements of the XC kernel f±iajb notation introduced in Ref. 59

where esl x/c,σ (r) are (semi-)local exchange/ correlation energy densities and eex x,σ (r) is the exact exchange energy density Z 1X σ σ σσ eex (r) = − D D wpqrs (r, r0 )d3 r0 (2) x,σ ps rq 2 pqrs

φpσ (r)φqσ (r)φrσ0 (r0 )φsσ0 (r0 ) σσ wpqrs . (r, r0 ) = |r − r0 |

0

σσ σσ (A + B)iaσjbσ0 = δij δab δσσ0 (σa − σi ) + 2Jiajb + f+iajb (9a)

(1)

0

(8)

iaσ

σ

Z

(6)

and the kinetic-energy density is defined as

G[X, Y, Ω] =

2

Page 2 of 18

∂ 2 Exc . σ ∂D σ 0 ∂Dpq rs

(11)

While for conventional functionals, the two kernel terms in Eq. 10 are identical and may be summarized yielding σσ 0 σσ 0 σσ 0 f+pqrs = 2 fpqrs and f−pqrs = 0, this does not hold for LHs, σσ 0 σσ 0 σσ 0 i.e. generally f+pqrs 6= 2 fpqrs and f−pqrs 6= 0. The reason is that weighting of the exact-exchange energy density terms by the LMF precludes the required interchange of MO indices. 34 This observation will be revisited further below. The excited-state gradients can be expressed as the derivative of Ω w.r.t. a change ξ in the set of nuclear coordinates Ωξ = Gξ [X, Y, Ω]. Since the involved MOs are usually expanded in an atomic-orbital basis {χµ (r)} with MO coefficients Cµpσ , X φpσ (r) = Cµpσ χµ (r) , (12)

(4)

µ

where ρσ (r) is the ground-state electron density X σ ρσ (r) = Dpq φpσ (r)φqσ (r) ,

ξ

the evaluation of G via the chain-rule of differentiation gives (5) Ωξ =

pq


∂G ∂G ∂Cµpσ dG = + , dξ ∂ξ ∂Cµpσ ∂ξ

(13)

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


where we omitted index summations for clarity. Since the handling of MO coefficient derivatives is rather involved, a fully variational Lagrangian L has been suggested, 8 X L [X, Y, Ω, C, Z, W ] =G[X, Y, Ω] + Ziaσ Fiaσ

equation, X

−

Wpqσ (Spqσ − δpq ) ,

(14)

Riaσ =

X

pqσ p≤q

b

∂L = Fiaσ = 0. ∂Ziaσ

j

+ − [T ] [(X − Y )] + Hiaσ +(X − Y )jaσ Hjiσ X σσ 0 σ 00 + g+iajbkc (X + Y )jbσ0 (X + Y )kcσ00

(15)

jbσ 0 kcσ 00

Spqσ are matrix elements of the overlap matrix, which are fixed by the Lagrange multiplier W such that orthonormality of the MOs is preserved ∂L = Spqσ − δpq = 0 . ∂Wpqσ

where we used for arbitrary matrices M the linear transformations i X h σσ0 + σσ 0 Hpqσ [M ] = 2 Jpqrs + f+pqrs Mrsσ0 , (24a) rsσ 0 − Hpqσ [M ]

(18)

As shown in detail in Ref. 8, the determination of Z and W starting from Eq. 17 involves several derivatives w.r.t. the MO coefficients. While most of them are unaffected by the more general LH scheme and can be adopted for this work, the derivatives of the exchange-correlation part need closer scrutiny (cf. Sec. S2). Similar to the particularities σσ 0 (Eq. 10) for LHs it encountered in the definition of f±pqrs xc turns out, that derivatives of potential terms Vrsσ 0 lead to two XC kernel terms with a flipped index pair,

µνσ

0

00

+

00

1 2

1 + 2

µνσ

ξ Jµνκλ [Pµνσ Dκλσ0 + (X + Y )µνσ (X + Y )κλσ0 ]

X

(ξ)

f+µνσκλσ0 (X + Y )µνσ (X + Y )κλσ0

µνσκλσ 0

X

(ξ)

f−µνσκλσ0 (X − Y )µνσ (X − Y )κλσ0 , (25)

µνσκλσ 0

where all matrices are expressed in their atomic-orbital representation, and hµν is the core Hamilton matrix.

(20) i.e. the third functional

∂ 3 Exc , σ ∂D σ 0 ∂D σ 00 ∂Dpq rs tu

(24b)

µνσκλσ 0

(19)

σ 0 σ 00 X ∂frstu σσ 0 σ 00 σσ 0 σ 00 σσ 0 σ 00 σσ 0 σ 00 Cµqσ = gqprstu + gpqrstu = g+pqrstu 6= 2 gpqrstu , ∂C µpσ µ

σσ σ gpqrstu =

µνσ

X

+

σ σ and derivatives of kernel terms fpqrs lead to two XC hyperkernel terms with a flipped index pair,

σσ 0 σ 00 gpqrstu ,

0

σσ Mrsσ0 , f−pqrs

and where T is the unrelaxed difference density defined as in Ref. 8. Note that the six-index XC hyper-kernel is never calculated as such and is instead directly contracted with the transition density matrices in the atomic-orbital basis as outlined in Sec. 3.3. Having solved for Z allows to use the remaining relations gained from Eq. 17 to obtain the relaxed one-particle difference density matrix P = T + Z and solve for W . Local hybrid specific changes in the solutions for W are similar to those in the definition of R (Eq. 23) and are thus relegated to Appendix A. In the final excited-state gradient equation (derived proceeding from Eq. 18) the XC kernel derivative is affected by the LH scheme following an entirely analogous reasoning, X ξ X ξ X xc (ξ) Ωξ = hµν Pµνσ − Sµνσ Wµνσ + Vµνσ Pµνσ

=0

X ∂V xc 0 σσ 0 σσ 0 σσ 0 σσ 0 rsσ + fpqrs = f+pqrs 6= 2 fpqrs , Cµqσ = fqprs ∂Cµpσ µ

X rsσ 0

the excited-state gradients may be evaluated as an explicit derivative w.r.t. ξ (hereafter indicated by superscript (ξ)),

0

=

(17)

dL ∂L ∂L ∂Cµpσ ∂L = + = . dξ ∂ξ ∂Cµpσ ∂ξ ∂ξ | {z }

00

σσ σ g−iajbkc (X − Y )jbσ0 (X − Y )kcσ00 , (23)

jbσ 0 kcσ 00

(16)

∂L = 0, ∂Cµpσ

0

X

+

If Z and W are chosen such that L becomes stationary w.r.t. the MO coefficients,

where the XC hyper-kernel derivative, is defined as

+ (X + Y )ibσ Habσ [(X + Y )]

− +(X − Y )ibσ Habσ [(X − Y )] X + [(X + Y )] − (X + Y )jaσ Hjiσ

where Fiaσ are matrix elements of the occupied-virtual block of the Fock matrix, which are enforced to be zero by the Lagrange multiplier Z,

Ωξ =

(22)

becomes

iaσ

X

Zjbσ0 (A + B)iaσjbσ0 = −Riaσ ,

jbσ 0

3

Implementation

Our implementation of excited-state gradients for LH functionals into a development branch of the TURBOMOLE 7.3 package is based on the existing EGRAD program, 8,60–65 which so far has allowed use of the local density approximation (LDA), generalized gradient approximation (GGA) and hybrid-GGA functionals. Routines for the initial determination of the excitation energy Ω and excitation vectors (X ± Y ) were readily accessible from the existing ESCF im-

(21)

and where we used a sign convention analogous to Eq. 10. The LH XC hyper-kernel will be derived explicitly in Section 3.3. Then the right-hand side R of the so-called Z vector



plementation for LH vertical excitation energies. 33 For the subsequent solving of the Z-vector equation, the right-hand side R has to be set up and new routines for the calculation of the LH hyper-kernel contracted with two transition density matrices have to be developed (cf. Eq. 23). For the calculation of Coulomb and XC kernel integrals (Eqs. 24a24b) the existing RI-J 62 and LH kernel routines were used. The final evaluation of the excited-state gradients (Eq. 25) requires an implementation of the gradient of the LH potential contracted with the relaxed one-particle difference density matrix and of the gradient of the LH kernel contracted with two transition density matrices, while the remaining contributions are again covered by the already available implementation. For the evaluation of LH derivatives, several non-standard integrals, i.e. LMF-weighted EXX integrals, have to be computed. The adaptation of the semi-numerical COSX scheme, 13,66 is ideally suited for this task and proved to be efficient and reliable in the SCF, 31 the TDDFT, 33 and the ground-state gradient 32 implementations of LHs. We will first introduce the general idea of this scheme and some convenient shorthand notations used in the subsequent equations for the above-mentioned LH derivatives.

3.1

clear documentation, here we introduce some helpful shorthand notations. The basis function vector   χ1 (ri )  χ2 (ri )    (30) X(ri ) =   ..   . χNbasis (ri ) helps to condense orbital summations and to rewrite Eq. 28 conveniently using a matrix-vector notation, X σ σ X Dνκ Dλµ qi · aσ (ri ) · χκ (ri )χλ (ri )Λµν (ri ) =

X

=aσ · X> Dσ ΛDσ X ,

U(1) = (X + Y ) , U

= (X − Y )

(32a) (32b)

That is, its action on an LMF a gives: ∂a ∂a ∂a b ∂µνσ a =Xµ Xν + ∇(Xµ Xν ) · 2∇ρσ + ∇ρσ0 ∂ρσ ∂γσσ ∂γσσ0 1 ∂a + ∇Xµ · ∇Xν . (34) 2 ∂τσ bξ is defined Analogously, the (semi-)local gradient operator ∇ as X (ξ) ∂ bξ = ∇ . (35) Q ∂Q Q∈Q

κλµν

κλµν

(−1)

in the following. For the derivatives of the LMF and (semi-)local energy σ densities w.r.t. some density matrix Dµν we extend the idea b of a semi-local potential operator ∂µνσ . 33 Given a set of (semi-)local quantities Q = {ρσ , γσσ , γσσ0 , τσ }, it is defined as X ∂Q ∂ ∂bµνσ = . (33) σ ∂Q ∂Dµν Q∈Q

It is calculated analytically at each given grid point. Then, an LMF-weighted four-center integral may be approximated as X σ σ ZZ Dνκ Dλµ aσ (r) · wκλµν (r, r0 )d3 rd3 r0 X

(31)

where here and below the dependence on grid points ri and the grid summation is omitted for brevity. To avoid confusion of X with the transition-density matrices (X ± Y ) we will use

where qi is the grid weight of the i-th grid point positioned at ri . We define the Λ matrix (also known as A Matrix in the literature) as Z χµ (r)χν (r) 3 Λµν (ri ) = d r. (27) |ri − r|

σ σ Dνκ Dλµ

σ σ χλ (ri )Dλµ Λµν (ri )Dνκ χκ (ri )

κλµν

Semi-numerical Integration

X

X

qi · aσ (ri ) ·

i

i

qi · aσ (ri ) · χκ (ri )χλ (ri )Λµν (ri ) . (28)

For higher-order derivatives we define a whole range of mixed and partly contracted operators

i

A convenient bonus of this “semi-numerical” scheme is the 4 reduction of the formal Nbasis scaling of the fully analytical 2 four-center integral evaluation to a formal Ngrid · Nbasis scaling. Efficient prescreenings can reduce the overall computational costs and enable an asymptotic linear scaling. 33,66,67 For the calculation of gradients of wκλµν (r, r0 ), we additionally define the gradient of a basis function ∇ξ χµ (r) = −∇χµ (r) and the Λ0 matrix: Z ∇ξ χµ (r0 )χν (r0 ) 3 0 d r . (29) Λ0µν (ri ) = |ri − r0 |

3.2

i

κλµν

In a DFT code the use of numerical quadrature is of course well established for spatial integrations, i.e. some matrix element in the atomic orbital basis Vµν is approximated as Z X Vµν = vµν (r)d3 r ≈ qi · vµν (ri ) , (26)

≈

Page 4 of 18

∂b(1) a =

X

(1) ∂bκλσ0 a · Uκλσ0 ,

κλσ 0

c (1) ∂∂ µνσ a =

X

(1) ∂bµνσ ∂bκλσ0 a · Uκλσ0 ,

κλσ 0

c ∂∂

(1)(1)

a=

X

(1) (1) ∂bησ0 ∂bκλσ00 a · Uησ0 Uκλσ00 ,

ησ 0 κλσ 00

d (1)(1) ∂∂∂ µνσ a

=

X

(1) (1) ∂bµνσ ∂bησ0 ∂bκλσ00 a · Uησ0 Uκλσ00 ,

ησ 0 κλσ 00

d ∇ ξ∂

Short-Hand Notation

(1)

a=

X

bξ ∂bκλσ0 a · U (1) 0 , ∇ κλσ

κλσ 0

Due to the product-rule of differentiation, higher-order LH derivatives are tedious to derive and lengthy to denote. For

[ ∇ ξ ∂∂

(1)(1)

a=

X ησ 0 κλσ 00


bξ ∂bησ0 ∂bκλσ00 a · U (1) 0 U (1) 00 , ∇ ησ κλσ

(36)

Page 5 of 18 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 These follow from exploiting the relations Dα = Dβ for (±) (±) restricted ground states and Uα = Uβ for singlet or

with their explicitly expanded forms given in Supporting Information (Eqs. S11, S17).

(±)

Uα

3.3

Hyper-Kernel and Kernel Gradient of Local Hybrid Functionals

3.4

κλη σ 0 σ 00

h i (1) = − 2 ∂b(1) aσ · Xµ (ΛU(1) σ X)ν + (ΛUσ X)µ Xν c (1)(1) aσ · [Xµ (ΛDσ X)ν + (ΛDσ X)µ Xν ] − ∂∂ i h X (1) (1) c (1)(1) esl −2 ∂bµνσ aσ0 · X> U 0 ΛU 0 X + ∂∂ x,σ 0 σ

σ0

−2

X

−4

X

−4

X

−2

X

+2

X

c (1)(1) aσ0 ∂bµνσ esl x,σ 0 · ∂∂

σ0

i h (1) > σ0 c (1) b(1) esl ∂∂ µνσ aσ 0 · X D ΛUσ 0 X + ∂ x,σ 0

σ0 sl c (1) b(1) aσ0 ∂∂ µνσ ex,σ 0 · ∂

σ0

σ0

0 1 > σ0 d (1)(1) ∂∂∂ X D ΛDσ X + esl µνσ aσ 0 · x,σ 0 2 sl d (1)(1) [1 − aσ0 ] · ∂∂∂ µνσ ex,σ 0

σ0 (1)(1)

d µνσ esl + 2 ∂∂∂ c .

(37)

For the antisymmetric term, most derivative terms cancel each other due to the symmetry of the (semi-)local potential operator on the one hand and due to the antisymmetry of the U(−1) matrix on the other hand: X X σσ0 σ00 (−1) (−1) g−µνκλη Uκλσ0 Uησ00 κλη σ 0 σ 00

=−2

X

h i (−1) (−1) ∂bµνσ aσ0 · X> Uσ0 ΛUσ0 X .

for triplet excitations.

Performance Considerations

Computationally, the calculation of the Λ and Λ0 matrices is the most demanding step. As explained in detail in Ref. 32, this is done using the Rys scheme. 68 Calculation of the Λ0 matrix takes approximately six times the computation time and memory of the Λ matrix calculation, because compared to Λ, Λ0 is not symmetric and each element of Λ0 is in fact a three-dimensional vector. It has already been seen for LH ground-state gradients 32 that this affects somewhat the absolute computational performance compared to the conventional analytical implementation of GH gradients. Still, the semi-numerical implementation showed vastly superior relative scaling with basis-set size. The savings in computational time by the use of S- and P-screenings turned out to be noticeable, reducing the overall costs by around 10% for medium-sized molecules and basis sets. Their impact grows with increasing basis-set and system size. It must be noted, however, that the P-junction screening for negligible density-matrix related terms (e.g. the vector Dσ X) becomes less efficient the more density matrices are considered. As opposed to ground-state gradients, excited-state gradients require the contraction of not only Dσ X with the Λ or Λ0 (1) (−1) matrix but also of Uσ X, Uσ X and Pσ X, making it less likely that no matrix element contributes and thus some Λ matrix calculations could be neglected. The matrix-vector multiplications can become a timecritical step and were thus extensively optimized in the LH TDDFT code. 33 However, compared to the calculation of VEEs that typically involve a considerable number of states and associated transition-density matrices U(±1) , only one state is considered in an excited-state optimization. Hence, the number of executed matrix-vector multiplications is fixed and is typically smaller than in a standard calculation of several VEEs. The efficiency of the matrix-vector multiplications is thus expected to be of secondary importance. Indeed, it was found that approximately 70% of the total computation time of a typical EGRAD run is taken up by the calculation of the Λ and Λ0 matrices. Overall, in its current state the performance of the implementation can be expected to leave room for improvement, and we plan to employ improved prescreenings and efficient parallelization. Hence, timings are not discussed in the main text. To give a rough estimate, an EGRAD run using the conventional analytical implementation of a GH takes about 84 min. for an organic molecule with 49 light atoms and triple-ζ basis set (1380 primitives) to complete on one core of an Intel® Core i7-7700 CPU @ 3.60GHz × 8. The same calculation takes approx. 2.5 times as long (216 min.) for an LH functional. Both calculations take about four to five times longer than their corresponding ground-state calculations. Details can be found in Supporting Information (Tab. S1).

Using the above notation, the hyper-kernel contracted with two symmetric transition-density matrices reads: X X σσ0 σ00 (1) (1) g+µνκλη Uκλσ0 Uησ00

σ

(±)

= −Uβ

(38)

σ0

Note that such simplifications do not hold if the currentdensity response 59 is included for the general case of a τ dependent functional. Results so far have revealed that the additional current-density terms hardly affect VEEs. 34 We thus decided at this point to neglect them for the excitedstate gradients to benefit from the significantly reduced theoretical and computational effort in the calculation of the antisymmetric part. We carefully examine the validity of this approximation in Section 5. With knowledge of the LH hyper-kernel the kernel gradients are derived by formally replacing certain derivative operators. Note that opposed to Λ, the Λ0 matrix is not symmetric which rules out the symmetry-based rearrangements used in the derivation of Eq. 37 and causes a more lengthy form of the non-local EXX-like terms. By the same token, an additional term remains in the antisymmetric kernel gradient. These equations, as well as the equation for the potential gradient contracted with the relaxed one-particle P xc (ξ) difference-density matrix, i.e. µνσ Vµνσ Pµνσ , are given in Appendix B for clarity. We note in passing that the presented equations were derived in the spin-unrestricted formalism but are easily transformed to the corresponding singlet and triplet equations.

4

Computational Details

The accuracy and correctness of the implementation has been validated using the test set composed in Ref. 8. It covers structural parameters, AAEs, and excited-state vibrational frequencies for organic and inorganic molecules with



a wide variety of bond types, excited-state character, and spin multiplicities. The reference values are gas-phase experimental data. We excluded the 1 1 B1 state of CCl2 from our study, as it showed non-real instabilities for all considered LH functionals. If not stated otherwise, all calculations have been performed using def2-QZVPPD basis sets. 69 Compared to the aug-TZVPP basis set used in Ref. 8, this basis set is of quadruple-ζ quality and has been optimized for the calculation of molecular properties. While impractically large for many applications, this basis set is feasible for the given test set, and we thus exploited the benefit of essentially working at the orbital-basis limit. The 2 1 B1 excited state of the H2 O molecule exhibits distinct Rydberg character and thus required additional diffuse f - and g-functions. The corresponding exponents were taken from Dunning’s augcc-pVQZ basis set. 70 The dependence of the semi-numerical EXX integration results on the size of the quadrature grid was carefully checked, and a quantitative evaluation is given in Sec. 5.1. Independent of this evaluation, large and diffuse integration grids (TURBOMOLE grid size 5, diffuse 2) were used in producing the benchmark results, to preclude any errors from the semi-numerical integration. Additionally, neither S- nor P-junctions were used for an integral prescreening in the semi-numerical scheme. For the optimization of the ground- and excited-state structures very tight convergence criteria were used, i.e. ground-state calculations were converged to energy changes 6 10−10 Hartree and density-matrix changes 6 10−9 a.u. Excitation vectors were converged to residual excitationvector norms 6 10−7 a.u. The maximum gradient norms were converged to 6 10−4 a.u. The large basis set, grid size, and strict convergence criteria are by no means recommendations for practical applications but were chosen to avoid error compensation and to hence unveil the true errors of a given functional. For the analysis of the influence of EXX admixtures on the excited states of formaldehyde (cf. Sec. 6.2) slightly looser thresholds (ground state density 6 10−8 a.u., residual excitation-vector norms 6 10−6 a.u.) and grid size 3 (diffuse 1) were used for performance reasons. Harmonic vibrational frequencies of the ground and excited states were calculated using the TURBOMOLE NumForce script, i.e. by means of numerical differentiation (first derivative, central differences) of the gradients. The displacement of the atoms generally breaks the molecular symmetry, and calculations were thus performed in C1 . However, this led to non-real instabilities for the displaced structure with some functionals. Since this was only observed for the 1 2 Σ+ state of NO, the 1 2 Π state of ScO and the 1 4 Π state of VO, i.e. some linear molecules, the calculation of the corresponding harmonic vibrational frequencies was carried out by manually displacing one atom along the molecular axis and exploiting full point group symmetry. We also included derivatives of the quadrature weights to improve the numerical accuracy of gradients in displaced structures during frequency calculations and to improve the convergence of structure optimizations. In general, the resolution of the identity (RI-J) 62 approximation was used for the Coulomb part, with the only exceptions being excited-state gradient calculations of triplet excited states, which were empirically found to be unstable with RI-J. We used standard TURBOMOLE auxiliary basis sets (def2-universal). 71

Page 6 of 18

We studied the performance of the following four LH functionals, which had already been extensively tested for ground-state calculations, 48,72–74 and for vertical excitation energies, 33,34,36 and which are available in TURBOMOLE since release 7.2. Lh07s-SVWN 49,72 is based on a spinchannel s-LMF (aσ (r) = erf (β · sσ (r)), where sσ (r) is the dimensionless density gradient and β = 0.22) and uses Slater exchange 75 and VWN correlation. 76 Lh07t-SVWN 48 uses the same exchange and correlation functionals and a scaled (b = 0.48) spin-channel t-LMF as given in Eq. 4. For Lh12ctSsirPW92 77 and Lh12ct-SsifPW92 77 a common version of the t-LMF (b = 0.646 and b = 0.709, respectively) was used, i.e. the LMF was constructed from the corresponding total spin quantities. Both functionals use Slater exchange and a short-range self-interaction-reduced (sir) or self-interactionfree (sif) version of the PW92 correlation functional. 77,78 Additionally, we investigated the GGA functional PBE 79 and the GH PBE0; 80 the latter both with conventional analytical EXX integration and with the semi-numerical scheme for comparison. The implementation of the τ -dependence as described in Sec. 3 also enabled the testing of the meta-GGA functional TPSS 81,82 and of the meta-GGA GH TPSSh 83 (implemented with semi-numerical EXX integration).

5

Technical Evaluation

As the intricate form of the third-order LH derivatives makes them very prone to errors we carefully checked the correctness of their implementation using finite-difference techniques. We found the numbers produced by our EGRAD implementation to be in excellent agreement with these numerical reference values (typically |Ωξimpl. − Ωξfin.-diff. | < 10−6 a.u.). Some further details can be found in Sec. S5 of Supporting Information. In the following, the impact of the introduced approximations, i.e. the semi-numerical integration and the neglect of the current-density response, is scrutinized.

5.1

Grid Dependence

The results obtained within a semi-numerical scheme naturally depend on the quadrature grid. It was already shown for ground-state calculations 31,32 and VEEs 33 that standard TURBOMOLE DFT grids, i.e. radial Chebyshev and spherical Lebedev grids, 84 are sufficient. We followed this path and were mainly concerned with the adequate grid size, which crucially affects the overall computation time as it determines the number of Λ and Λ0 matrices to be calculated. It is already known that small grids (TURBOMOLE grid size 1) give satisfactory vertical excitation energies with the exception of Rydberg excitations for which slightly larger and more diffuse grids are required. 33 To investigate the grid convergence of the semi-numerical integration for excited-state structural parameters, we performed calculations with a small grid (TURBOMOLE grid size 1) and a medium-sized, slightly more diffuse grid (TURBOMOLE grid size 2, diffuse 1) at the PBE0/ def2-QZVPPD level of theory, both in conventional (analytical) and in seminumerical implementation, compared against analytical results for an extremely large and diffuse grid (TURBOMOLE grid size 5, diffuse 2). For most of the 63 parameters the small grid gives deviations well below 0.1 pm in bond lengths and below 0.1 deg in bond angles, respectively. Such small deviations are negligible compared with typical errors due to functional or basis set and thus can be safely ignored.


Page 7 of 18 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 1: Effects of the current-density response on excitedstate bond lengths (in pm) for diatomic molecules at the Lh07t-SVWN/ def2-QZVPPD level of theory.

For parameters above this threshold, deviations are given in Fig. 1. While the deviation of the data points from the zero line is related to the total grid error, differences between the conventional (squares) and semi-numerical (crosses) results at the same grid size reveal the additional error introduced by the semi-numerical scheme. No systematic over- or underestimation of structural parameters is observed when using smaller grids. For the small grid (red labels), several parameters show notable but small total absolute deviations of around 0.5 pm or 0.5 deg. Somewhat larger effects include a distinct 3.0 deg deviation for the out-of-plane dihedral angle in the 1 3 A00 state of CH2 S and a 1.8 pm deviation for the bond length in the 1 1 Σ+ u state of Mg2 . These deviations have to be put into perspective since the potential energy surfaces of these states are very flat along the respective coordinates. Nonetheless, these cases expose the limitations of the small grid. However, even a moderate increase of the grid size (green labels) eliminates almost all deviations. We note in passing that diffuse grids are of course only required for diffuse states (excited states with notable Rydberg character) and do not necessarily entail an increased number of grid points as a tuning of the

Mol.

State

AsF BeH BeO BF BH CO CO CuH Li2 Mg2 N2 N2 N2 N2 NH NO P2 SiO ScO VO

1 3Π 1 2Π 1 1Π 1 1Π 1 1Π 1 1Π 1 3Π 2 1 Σ+ 1 1 Σ+ u 1 1 Σ+ u 1 1 ∆u 1 1 Πg 1 1 Σ− u 1 3 Πg 3 1 Π 1 2 Σ+ 1 1 Πg 1 1Π 1 2Π 1 4Π

MAE MAE*c AsF 13Π re 1

BeO 1 Π re H2O 11B1 O-H

senu. grid 2 diffuse conv. grid 2 diffuse senu. grid conv. grid

1 1 1 1

kb 149 249 501 673 316 996 1334 173 13 17 1074 1291 1020 1409 561 2904 374 513 520 664

∆re cur re

reno cur

recur

∆re

197.70 133.10 144.35 130.68 121.47 122.91 119.94 158.44 313.20 323.08 127.38 121.16 126.57 120.20 104.42 104.92 198.37 160.17 168.33 163.03

197.08 133.10 143.88 130.33 121.23 122.62 119.80 158.03 307.64 325.82 127.27 121.06 126.47 120.10 104.21 104.67 198.23 159.91 168.19 163.10

0.62 0.00 0.47 0.35 0.24 0.29 0.14 0.41 5.56 −2.74 0.11 0.10 0.10 0.10 0.21 0.25 0.14 0.26 0.14 −0.07

0.31% 0.00% 0.32% 0.27% 0.20% 0.24% 0.12% 0.26% 1.81% −0.84% 0.09% 0.09% 0.08% 0.08% 0.20% 0.24% 0.07% 0.16% 0.08% −0.05%

0.62 0.22

0.28% 0.16%

a

Optimizations without inclusion of the current-density response (reno cur ) were performed using the new gradient implementation, while optimizations with inclusion of the currentdensity response (recur ) were performed using a linear-search algorithm for excited-state single-point energies. b Excited-state force constants (in kg s−2 ) from calculations without current-density response are given as an estimate of the steepness of the potential-energy curves. c Calculated excluding Li2 and Mg2 .

Li2 11Σ+ u re Mg2 11Σ+ u re P2 11Πg re Propynal 21A C-O 2

ScO 1 Π re SiF2 11B1 Si-F 1

SiO 1 Π re

radial mapping parameter is often sufficient. 33 Errors introduced by the semi-numerical scheme are mostly small with the largest deviation being 0.4 deg for the 1 1 B1 state of SiF2 for the small grid (excluding here the two most severe absolute outliers). They are further reduced when switching to the medium-sized grid. We thus conclude, that the errors of the semi-numerical scheme are well controlled and the choice of grid requires no particular additional care for LHs compared to standard DFT functionals.

VO 14Π re 1

CH2O 1 A" Pσ X h i X (P ) − ∂b aσ · X> Dσ ΛDσ ∇ξ X + X> Dσ Λ0 Dσ X

−(X − Y )iaσ (X + Y )ibσ ] X σ + i [(X + Y )iaσ (X + Y )ibσ

σ

i

+(X − Y )iaσ (X − Y )ibσ ] ,

(40)

−

X

−

X

−

X

+

X

i h ∇ξ aσ · X> Dσ ΛPσ X + ∂b(P ) esl x,σ

σ

and

∂b(P ) aσ · ∇ξ esl x,σ

σ

Wiaσ =

X

+ (X + Y )jaσ Hjiσ [(X + Y )]

j

− [(X − Y )] + σi Ziaσ . +(X − Y )iaσ Hjiσ

σ

(41)

(P ) 1 > σ σ sl d ∇ ∂ a · X D ΛD X + e σ ξ x,σ 2 (P ) sl d [1 − aσ ] · ∇ ex,σ ξ∂

σ

B

Gradient of the LH Kernel and Potential

(P ) sl d ec , +∇ ξ∂

where the operators ∂b(P ) , . . . are defined as given in Eq. 36 but using the (symmetric) relaxed one-particle difference density matrix P instead of a transition density matrix.

The (symmetric) kernel gradient part reads: 1 X X (ξ) (1) (1) f+µνσκλσ0 Uµνσ Uκλσ0 2 0 σσ κλµν h i X (1) > (1) 0 (1) =−2 aσ · X> U(1) σ ΛUσ ∇ξ X + X Uσ Λ Uσ X

References (1) Baryshnikov, G.; Minaev, B.; ˚ Agren, H. Theory and Calculation of the Phosphorescence Phenomenon. Chem. Rev. 2017, 117, 6500–6537.

σ

−2

X

h > σ (1) ∂b(1) aσ · X> Dσ ΛU(1) σ ∇ξ X + ∇ξ X D ΛUσ X

σ

i > σ 0> (1) +X> Dσ Λ0 U(1) σ X + X D Λ Uσ X i h X (1)(1) c − ∂∂ aσ · X> Dσ ΛDσ ∇ξ X + X> Dσ Λ0 Dσ X

(2) Ashfold, M. N. R.; Ingle, R. A.; Karsili, T. N. V.; Zhang, J. Photoinduced C–H bond fission in prototypical organic molecules and radicals. Phys. Chem. Chem. Phys. 2019, 21, 13880–13901.

σ

−

X

−

X

h i (1) c (1)(1) esl bξ aσ · X> U(1) ∇ σ ΛUσ X + ∂∂ x,σ

(3) Petsalakis, I. D.; Georgiadou, D. G.; Vasilopoulou, M.; Pistolis, G.; Dimotikali, D.; Argitis, P.; Theodorakopoulos, G. Theoretical Investigation on the Effect of Protonation on the Absorption and Emission Spectra of Two Amine-Group-Bearing, Red “PushPull” Emitters, 4-Dimethylamino-4’-nitrostilbene and 4-(dicyanomethylene)-2-methyl-6-p-(dimethylamino) styryl-4H-pyran, by DFT and TDDFT Calculations. J. Phys. Chem. A 2010, 114, 5580–5587.

σ

c bξ esl ∇ x,σ · ∂∂

(1)(1)

aσ

σ

−2

X

−2

X

h i (1) d b(1) esl ∇ aσ · X> Dσ ΛU(1) ξ∂ σ X+∂ x,σ

σ (1) sl d ∇ ex,σ · ∂b(1) aσ ξ∂

σ

−

X

+

X

[ ∇ ξ ∂∂

(1)(1)

σ

[1 − aσ ] ·

aσ ·

1 > σ X D ΛDσ X + esl x,σ 2

(4) Smith, M. B.; Michl, J. Singlet Fission. Chem. Rev. 2010, 110, 6891–6936.

(1)(1) sl [ ∇ ex,σ ξ ∂∂

σ

[ +∇ ξ ∂∂

(1)(1) sl ec

.

(5) Gross, E. K. U.; Kohn, W. In Density Functional Theory of Many-Fermion Systems; L¨ owdin, P.-O., Ed.; Advances in Quantum Chemistry; Academic Press: Cambridge, Massachusetts, 1990; Vol. 21; pp 255–291.

(42)

The (antisymmetric) kernel gradient part reads: 1 X X (ξ) (−1) (−1) f−µνσ0 κλσ00 Uµνσ0 Uκλσ00 2 00 0 σ σ κλµν X bξ aσ · X> U(−1) =− ∇ ΛU(−1) X σ σ

(6) Caillie, C. V.; Amos, R. D. Geometric derivatives of excitation energies using SCF and DFT. Chem. Phys. Lett. 1999, 308, 249–255.

σ

−2

X

(44)

h i aσ · X> U(−1) ΛU(−1) ∇ξ X + X> U(−1) Λ0 U(−1) X . σ σ σ σ

σ

(43)

(7) Caillie, C. V.; Amos, R. D. Geometric derivatives of density functional theory excitation energies using gradient-corrected functionals. Chem. Phys. Lett. 2000, 317, 159–164.


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


(8) Furche, F.; Ahlrichs, R. Adiabatic time-dependent density functional methods for excited state properties. J. Chem. Phys. 2002, 117, 7433–7447.

(21) Santoro, F.; Jacquemin, D. Going beyond the vertical approximation with time-dependent density functional theory. WIREs Comput. Mol. Sci. 2016, 6, 460–486.

(9) TURBOMOLE V7.3 2018, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989-2007, TURBOMOLE GmbH, since 2007; available from http://www.turbomole.com.

(22) Br´emond, E.; Savarese, M.; Adamo, C.; Jacquemin, D. Accuracy of TD-DFT Geometries: A Fresh Look. J. Chem. Theory Comput. 2018, 14, 3715–3727. (23) Peach, M. J. G.; Benfield, P.; Helgaker, T.; Tozer, D. J. Excitation energies in density functional theory: An evaluation and a diagnostic test. J. Chem. Phys. 2008, 128, 044118.

(10) Liu, F.; Gan, Z.; Shao, Y.; Hsu, C.-P.; Dreuw, A.; Head-Gordon, M.; Miller, B. T.; Brooks, B. R.; Yu, J.G.; Furlani, T. R.; Kong, J. A parallel implementation of the analytic nuclear gradient for time-dependent density functional theory within the Tamm–Dancoff approximation. Mol. Phys. 2010, 108, 2791–2800.

(24) Winter, N. O. C.; Graf, N. K.; Leutwyler, S.; H¨ attig, C. Benchmarks for 0–0 transitions of aromatic organic molecules: DFT/B3LYP, ADC(2), CC2, SOS-CC2 and SCS-CC2 compared to high-resolution gas-phase data. Phys. Chem. Chem. Phys. 2013, 15, 6623–6630.

(11) Hutter, J. Excited state nuclear forces from the Tamm–Dancoff approximation to time-dependent density functional theory within the plane wave basis set framework. J. Chem. Phys. 2003, 118, 3928–3934.

(25) Zhao, Y.; Truhlar, D. G. Density Functional for Spectroscopy: No Long-Range Self-Interaction Error, Good Performance for Rydberg and Charge-Transfer States, and Better Performance on Average than B3LYP for Ground States. J. Phys. Chem. A 2006, 110, 13126– 13130.

(12) Sitt, A.; Kronik, L.; Ismail-Beigi, S.; Chelikowsky, J. R. Excited-state forces within time-dependent densityfunctional theory: A frequency-domain approach. Phys. Rev. A 2007, 76, 054501.

(26) Silva-Junior, M. R.; Schreiber, M.; Sauer, S. P. A.; Thiel, W. Benchmarks for electronically excited states: Time-dependent density functional theory and density functional theory based multireference configuration interaction. J. Chem. Phys. 2008, 129, 104103.

(13) Petrenko, T.; Kossmann, S.; Neese, F. Efficient timedependent density functional theory approximations for hybrid density functionals: Analytical gradients and parallelization. J. Chem. Phys. 2011, 134, 054116.

(27) Goerigk, L.; Grimme, S. Assessment of TD-DFT methods and of various spin scaled CIS(D) and CC2 versions for the treatment of low-lying valence excitations of large organic dyes. J. Chem. Phys. 2010, 132, 184103.

(14) Scalmani, G.; Frisch, M. J.; Mennucci, B.; Tomasi, J.; Cammi, R.; Barone, V. Geometries and properties of excited states in the gas phase and in solution: Theory and application of a time-dependent density functional theory polarizable continuum model. J. Chem. Phys. 2006, 124, 094107.

(28) Verma, P.; Wang, Y.; Ghosh, S.; He, X.; Truhlar, D. G. Revised M11 Exchange-Correlation Functional for Electronic Excitation Energies and GroundState Properties. J. Phys. Chem. A 2019, 123, 2966– 2990.

(15) Zeng, Q.; Liang, W. Analytic energy gradient of excited electronic state within TDDFT/MMpol framework: Benchmark tests and parallel implementation. J. Chem. Phys. 2015, 143, 134104.

(29) Jaramillo, J.; Scuseria, G. E.; Ernzerhof, M. Local hybrid functionals. J. Chem. Phys. 2003, 118, 1068– 1073.

(16) Si, D.; Li, H. Analytic energy gradient in combined time-dependent density functional theory and polarizable force field calculation. J. Chem. Phys. 2010, 133, 144112.

(30) Maier, T. M.; Arbuznikov, A. V.; Kaupp, M. Local hybrid functionals: Theory, implementation, and performance of an emerging new tool in quantum chemistry and beyond. Wiley Interdiscip. Rev.-Comput. Mol. Sci. 2018, 9, e1378.

(17) Chiba, M.; Tsuneda, T.; Hirao, K. Excited state geometry optimizations by analytical energy gradient of long-range corrected time-dependent density functional theory. J. Chem. Phys. 2006, 124, 144106.

(31) Bahmann, H.; Kaupp, M. Efficient Self-Consistent Implementation of Local Hybrid Functionals. J. Chem. Theory Comput. 2015, 11, 1540–1548.

(18) Seth, M.; Mazur, G.; Ziegler, T. Time-dependent density functional theory gradients in the Amsterdam density functional package: geometry optimizations of spin-flip excitations. Theor. Chem. Acc. 2010, 129, 331–342.

(32) Klawohn, S.; Bahmann, H.; Kaupp, M. Implementation of Molecular Gradients for Local Hybrid Density Functionals Using Seminumerical Integration Techniques. J. Chem. Theory Comput. 2016, 12, 4254– 4262.

(19) Send, R.; K¨ uhn, M.; Furche, F. Assessing Excited State Methods by Adiabatic Excitation Energies. J. Chem. Theory Comput. 2011, 7, 2376–2386.

(33) Maier, T. M.; Bahmann, H.; Kaupp, M. Efficient Seminumerical Implementation of Global and Local Hybrid Functionals for Time-Dependent Density Functional Theory. J. Chem. Theory Comput. 2015, 11, 4226– 4237.

(20) Laurent, A. D.; Jacquemin, D. TD-DFT benchmarks: A review. Int. J. Quantum Chem. 2013, 113, 2019– 2039.



Page 14 of 18

(46) Dierksen, M.; Grimme, S. Density functional calculations of the vibronic structure of electronic absorption spectra. J. Chem. Phys. 2004, 120, 3544–3554.

(34) Maier, T. M.; Bahmann, H.; Arbuznikov, A. V.; Kaupp, M. Validation of local hybrid functionals for TDDFT calculations of electronic excitation energies. J. Chem. Phys. 2016, 144, 0741061.

(47) Janesko, B. G.; Scuseria, G. E. Local hybrid functionals based on density matrix products. J. Chem. Phys. 2007, 127, 164117.

(35) Wen, J.; Havlas, Z.; Michl, J. Captodatively Stabilized Biradicaloids as Chromophores for Singlet Fission. J. Am. Chem. Soc. 2015, 137, 165–172.

(48) Bahmann, H.; Rodenberg, A.; Arbuznikov, A. V.; Kaupp, M. A thermochemically competitive local hybrid functional without gradient corrections. J. Chem. Phys. 2007, 126, 011103.

(36) Grotjahn, R.; Maier, T. M.; Michl, J.; Kaupp, M. Development of a TDDFT-Based Protocol with Local Hybrid Functionals for the Screening of Potential Singlet Fission Chromophores. J. Chem. Theory Comput. 2017, 13, 4984–4996.

(49) Arbuznikov, A. V.; Kaupp, M. Local hybrid exchangecorrelation functionals based on the dimensionless density gradient. Chem. Phys. Lett. 2007, 440, 160–168.

(37) Loos, P.-F.; Jacquemin, D. Chemically Accurate 0-0 Energies with Not-so-Accurate Excited State Geometries. J. Chem. Theory. Comput. 2019, 15, 2481–2491.

(50) Arbuznikov, A. V.; Kaupp, M. What can we learn from the adiabatic connection formalism about local hybrid functionals? J. Chem. Phys. 2008, 128, 214107.

(38) Tapavicza, E.; Bellchambers, G. D.; Vincent, J. C.; Furche, F. Ab initio non-adiabatic molecular dynamics. Phys. Chem. Chem. Phys. 2013, 15, 18336–18348.

(51) Janesko, B. G.; Scuseria, G. E. Parameterized local hybrid functionals from density-matrix similarity metrics. J. Chem. Phys. 2008, 128, 084111.

(39) Parthey, M.; Kaupp, M. Quantum-chemical insights into mixed-valence systems: within and beyond the Robin–Day scheme. Chem. Soc. Rev. 2014, 43, 5067– 5088.

(52) Perdew, J. P.; Staroverov, V. N.; Tao, J.; Scuseria, G. E. Density functional with full exact exchange, balanced nonlocality of correlation, and constraint satisfaction. Phys. Rev. A 2008, 78, 052513.

(40) G¨ uckel, S.; Gluyas, J. B. G.; El-Tarhuni, S.; Sobolev, A. N.; Whiteley, M. W.; Halet, J.-F.; Lapinte, C.; Kaupp, M.; Low, P. J. Iron versus Ruthenium: Clarifying the Electronic Differences between Prototypical Mixed-Valence Organometallic Butadiyndiyl Bridged Molecular Wires. Organometallics 2018, 37, 1432–1445.

(53) Haunschild, R.; Janesko, B. G.; Scuseria, G. E. Local hybrids as a perturbation to global hybrid functionals. J. Chem. Phys. 2009, 131, 154112. (54) Arbuznikov, A. V.; Bahmann, H.; Kaupp, M. Local Hybrid Functionals with an Explicit Dependence on Spin Polarization. J. Phys. Chem. A 2009, 113, 11898– 11906.

(41) Kaupp, M.; Karton, A.; Bischoff, F. A. [Al2 O4 ]= , a Benchmark Gas-Phase Class II Mixed-Valence Radical Anion for the Evaluation of Quantum-Chemical Methods. J. Chem. Theory Comput. 2016, 12, 3796–3806.

(55) Johnson, E. R. Local-hybrid functional based on the correlation length. J. Chem. Phys. 2014, 141, 124120.

(42) Klawohn, S.; Kaupp, M.; Karton, A. MVO10: A Gas-Phase Oxide Benchmark for Localization/Delocalization in Mixed-Valence Systems. J. Chem. Theory Comput. 2018, 14, 3512–3523.

(56) Schmidt, T.; Kraisler, E.; Makmal, A.; Kronik, L.; K¨ ummel, S. A self-interaction-free local hybrid functional: Accurate binding energies vis-` a-vis accurate ionization potentials from Kohn-Sham eigenvalues. J. Chem. Phys. 2014, 140, 18A510.

(43) Low, P. J.; G¨ uckel, S.; Gluyas, J. B. G.; Eaves, S. G.; Safari, P.; Yufit, D. S.; Sobolev, A. N.; Kaupp, M. A spectroscopic and computationally minimal approach to the analysis of charge transfer processes in conformationally fluxional mixed-valence and heterobimetallic complexes. Chem. Eur. J. 2019, 25, 8837– 8853.

(57) de Silva, P.; Corminboeuf, C. Local hybrid functionals with orbital-free mixing functions and balanced elimination of self-interaction error. J. Chem. Phys. 2015, 142, 074112. (58) Furche, F. On the density matrix based approach to time-dependent density functional response theory. J. Chem. Phys. 2001, 114, 5982–5992.

(44) Tapavicza, E.; Furche, F.; Sundholm, D. Importance of Vibronic Effects in the UV–Vis Spectrum of the 7,7,8,8-Tetracyanoquinodimethane Anion. J. Chem. Theory Comput. 2016, 12, 5058–5066.

(59) Bates, J. E.; Furche, F. Harnessing the metageneralized gradient approximation for timedependent density functional theory. J. Chem. Phys. 2012, 137, 164105.

(45) Dierksen, M.; Grimme, S. The Vibronic Structure of Electronic Absorption Spectra of Large Molecules: A Time-Dependent Density Functional Study on the Influence of “Exact” Hartree−Fock Exchange. J. Phys. Chem. A 2004, 108, 10225–10237.

(60) Furche, F.; Ahlrichs, R. Erratum: “Time-dependent density functional methods for excited state properties” [J. Chem. Phys. 117, 7433 (2002)]. J. Chem. Phys. 2004, 121, 12772.


Page 15 of 18 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 (75) Slater, J. C. A Simplification of the Hartree-Fock Method. Phys. Rev. 1951, 81, 385–390.

(61) Furche, F.; Rappoport, D. In Density functional methods for excited states: equilibrium structure and electronic spectra; Olivucci, M., Ed.; Theoretical and Computational Chemistry; Elsevier: Amsterdam, Netherlands, 2005; Vol. 16; Chapter III, pp 93–128.

(76) Vosko, S. H.; Wilk, L.; Nusair, M. Accurate spindependent electron liquid correlation energies for local spin density calculations: a critical analysis. Can. J. Phys. 1980, 58, 1200–1211.

(62) Rappoport, D.; Furche, F. Analytical time-dependent density functional derivative methods within the RI-J approximation, an approach to excited states of large molecules. J. Chem. Phys. 2005, 122, 064105.

(77) Arbuznikov, A. V.; Kaupp, M. Importance of the correlation contribution for local hybrid functionals: Range separation and self-interaction corrections. J. Chem. Phys. 2012, 136, 014111.

(63) Rappoport, D.; Furche, F. Lagrangian approach to molecular vibrational Raman intensities using timedependent hybrid density functional theory. J. Chem. Phys. 2007, 126, 201104.

(78) Perdew, J. P. Accurate and simple analytic representation of the electron-gas correlation energy. Phys. Rev. B 1992, 45, 13244–13249.

(64) Send, R.; Furche, F. First-order nonadiabatic couplings from time-dependent hybrid density functional response theory: Consistent formalism, implementation, and performance. J. Chem. Phys. 2010, 132, 044107.

(79) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 1996, 77, 3865–3868. (80) Adamo, C.; Barone, V. Toward reliable density functional methods without adjustable parameters: The PBE0 model. J. Chem. Phys. 1999, 110, 6158–6170.

(65) Furche, F.; Krull, B. T.; Nguyen, B. D.; Kwon, J. Accelerating molecular property calculations with nonorthonormal Krylov space methods. J. Chem. Phys. 2016, 144, 174105.

(81) Tao, J.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E. Climbing the Density Functional Ladder: Nonempirical Meta–Generalized Gradient Approximation Designed for Molecules and Solids. Phys. Rev. Lett. 2003, 91, 146401.

(66) Neese, F.; Wennmohs, F.; Hansen, A.; Becker, U. Efficient, approximate and parallel Hartree–Fock and hybrid DFT calculations. A ‘chain-of-spheres’ algorithm for the Hartree–Fock exchange. Chem. Phys. 2009, 356, 98–109.

(82) Perdew, J. P.; Tao, J.; Staroverov, V. N.; Scuseria, G. E. Meta-generalized gradient approximation: Explanation of a realistic nonempirical density functional. J. Chem. Phys. 2004, 120, 6898–6911.

(67) Laqua, H.; Kussmann, J.; Ochsenfeld, C. Efficient and Linear-Scaling Seminumerical Method for Local Hybrid Density Functionals. J. Chem. Theory Comput. 2018, 14, 3451–3458. (68) Dupuis, M.; Rys, J.; King, H. F. Evaluation of molecular integrals over Gaussian basis functions. J. Chem. Phys. 1976, 65, 111–116.

(83) Staroverov, V. N.; Scuseria, G. E.; Tao, J.; Perdew, J. P. Comparative assessment of a new nonempirical density functional: Molecules and hydrogenbonded complexes. J. Chem. Phys. 2003, 119, 12129– 12137.

(69) Rappoport, D.; Furche, F. Property-optimized Gaussian basis sets for molecular response calculations. J. Chem. Phys. 2010, 133, 134105.

(84) Treutler, O.; Ahlrichs, R. Efficient molecular numerical integration schemes. J. Chem. Phys. 1995, 102, 346– 354.

(70) Kendall, R. A.; Dunning, T. H.; Harrison, R. J. Electron affinities of the first-row atoms revisited. Systematic basis sets and wave functions. J. Chem. Phys. 1992, 96, 6796–6806.

(85) Jensen, P.; Bunker, P. R. The geometry and the inver˜ 1 A2 sion potential function of formaldehyde in the A 3 and a ˜ A2 electronic states. J. Mol. Spectrosc. 1982, 94, 114–125.

(71) Weigend, F. Accurate Coulomb-fitting basis sets for H to Rn. Phys. Chem. Chem. Phys. 2006, 8, 1057–1065.

ˇ Scalmani, G.; Jacquemin, D. Accurate (86) Budz´ ak, S.; Excited-State Geometries: A CASPT2 and CoupledCluster Reference Database for Small Molecules. J. Chem. Theory Comput. 2017, 13, 6237–6252.

(72) Kaupp, M.; Bahmann, H.; Arbuznikov, A. V. Local hybrid functionals: An assessment for thermochemical kinetics. J. Chem. Phys. 2007, 127, 194102.

˜ electronic (87) Huet, T.; Godefroid, M.; Herman, M. The A state of acetylene: Geometry and axis-switching effects. J. Mol. Spectrosc. 1990, 144, 32–44.

(73) Arbuznikov, A. V.; Kaupp, M. Towards improved local hybrid functionals by calibration of exchange-energy densities. J. Chem. Phys. 2014, 141, 204101.

(88) Ventura, E.; Dallos, M.; Lischka, H. The valenceexcited states T1 –T4 and S1 –S2 of acetylene: A highlevel MR-CISD and MR-AQCC investigation of stationary points, potential energy surfaces, and surface crossings. J. Chem. Phys. 2003, 118, 1702–1713.

(74) Theilacker, K.; Arbuznikov, A. V.; Kaupp, M. Gauge effects in local hybrid functionals evaluated for weak interactions and the GMTKN30 test set. Mol. Phys. 2015, 114, 1118–1127.



(89) K´ allay, M.; Gauss, J. Calculation of excited-state properties using general coupled-cluster and configurationinteraction models. J. Chem. Phys. 2004, 121, 9257– 9269. (90) Scott, A. P.; Radom, L. Harmonic Vibrational Frequencies: An Evaluation of Hartree–Fock, Møller– Plesset, Quadratic Configuration Interaction, Density Functional Theory, and Semiempirical Scale Factors. J. Phys. Chem. 1996, 100, 16502–16513. (91) Arbuznikov, A. V.; Kaupp, M. Localized hybrid exchange-correlation potentials for Kohn–Sham DFT calculations of NMR and EPR parameters. Int. J. Quantum Chem. 2005, 104, 261–271. (92) Arbuznikov, A. V.; Kaupp, M.; Bahmann, H. From local hybrid functionals to “localized local hybrid” potentials: Formalism and thermochemical tests. J. Chem. Phys. 2006, 124, 204102. (93) Schattenberg, C. J.; Maier, T. M.; Kaupp, M. Lessons from the Spin-Polarization/Spin-Contamination Dilemma of Transition-Metal Hyperfine Couplings for the Construction of Exchange-Correlation Functionals. J. Chem. Theory. Comput. 2018, 14, 5653–5672. (94) Ipatov, A.; Cordova, F.; Doriol, L. J.; Casida, M. E. Excited-state spin-contamination in time-dependent density-functional theory for molecules with open-shell ground states. J. Mol. Struct.:THEOCHEM 2009, 914, 60–73. (95) Myneni, H.; Casida, M. E. On the calculation of ∆hS 2 i for electronic excitations in time-dependent densityfunctional theory. Comput. Phys. Commun. 2017, 213, 72–91.


Page 16 of 18

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


TOC-Graphic


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

ACS Paragon Plus Environment

Page 18 of 18

Development and Implementation of Excited-State Gradients for Local

Recommend Documents