First-order derivative couplings between excited states from adiabatic TDDFT response theory Qi Ou, Gregory D. Bellchambers, Filipp Furche, and Joseph E. Subotnik Citation: The Journal of Chemical Physics 142, 064114 (2015); doi: 10.1063/1.4906941 View online: http://dx.doi.org/10.1063/1.4906941 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/142/6?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Derivative couplings between TDDFT excited states obtained by direct differentiation in the Tamm-Dancoff approximation J. Chem. Phys. 141, 024114 (2014); 10.1063/1.4887256 First-order nonadiabatic coupling matrix elements between excited states: A Lagrangian formulation at the CIS, RPA, TD-HF, and TD-DFT levels J. Chem. Phys. 141, 014110 (2014); 10.1063/1.4885817 Excitation energies with linear response density matrix functional theory along the dissociation coordinate of an electron-pair bond in N-electron systems J. Chem. Phys. 140, 024101 (2014); 10.1063/1.4852195 First-order nonadiabatic couplings from time-dependent hybrid density functional response theory: Consistent formalism, implementation, and performance J. Chem. Phys. 132, 044107 (2010); 10.1063/1.3292571 Transition strengths and first-order properties of excited states from local coupled cluster CC2 response theory with density fitting J. Chem. Phys. 127, 064107 (2007); 10.1063/1.2755778

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Mon, 02 Mar 2015 11:18:02

THE JOURNAL OF CHEMICAL PHYSICS 142, 064114 (2015)

First-order derivative couplings between excited states from adiabatic TDDFT response theory Qi Ou,1 Gregory D. Bellchambers,2 Filipp Furche,2,a) and Joseph E. Subotnik1,b)

1

Department of Chemistry, University of Pennsylvania, Philadelphia, Pennsylvania 19104, USA Department of Chemistry, University of California, Irvine, 1102 Natural Sciences II, Irvine, California 92697-2025, USA 2

(Received 25 November 2014; accepted 14 January 2015; published online 12 February 2015) We present a complete derivation of derivative couplings between excited states in the framework of adiabatic time-dependent density functional response theory. Explicit working equations are given and the resulting derivative couplings are compared with derivative couplings from a pseudo-wavefunction ansatz. For degenerate excited states, i.e., close to a conical intersection (CI), the two approaches are identical apart from an antisymmetric overlap term. However, if the difference between two excitation energies equals another excitation energy, the couplings from response theory exhibit an unphysical divergence. This spurious behavior is a result of the adiabatic or static kernel approximation of time-dependent density functional theory leading to an incorrect analytical structure of the quadratic response function. Numerical examples for couplings close to a CI and for well-separated electronic states are given. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4906941] I. INTRODUCTION

In the past few decades, great efforts in the field of quantum chemistry have been expended to study nonadiabatic processes. Going beyond the Born-Oppenheimer approximation, nonadiabatic processes are ubiquitous and cover many interesting modern topics—including charge transfer, electronic excitation quenching, and spin-forbidden transitions. And, at the bottom, modeling most of these processes requires computing the derivative coupling.1–6 Today, time-dependent density functional theory (TDDFT)7–10 is a mainstay of computational photochemistry. The popularity of TDDFT is owed to a compromise between accuracy and computational efficiency that holds up well in many (although not all) applications. However, the evaluation of TDDFT derivative couplings is complicated by the fact that interacting wavefunctions are inaccessible in TDDFT.11 This led to a plethora of approaches for evaluating both the ground-excited state couplings12–18 and excited-excited state couplings.19–24 In 2010, Send and Furche solved the ground-excited state problem definitively by relating the TDDFT derivative coupling to a residue from linear response theory and calculating the residue in a finite atomic orbital (AO) basis.18 The resulting coupling reduces to the exact expression derived by Chernyak and Mukamel12 in the complete basis set limit, which guarantees convergence to the exact result as better and better approximations to the time-dependent exchange-correlation (XC) potential are used. These developments have enabled efficient TDDFT-based nonadiabatic molecular dynamics simulations for systems in the first excited state.25

a)[email protected] b)[email protected]

0021-9606/2015/142(6)/064114/14/$30.00

In this article, our focus will be exclusively on excited state-excited state couplings. In this case, TDDFT recovers the correct dimensionality of a conical intersection (CI) branching plane26–28 and should thus be even more useful. To our knowledge, there have been three different proposed methods to evaluate TDDFT excited state derivative coupl∂ ings. First, Tavernelli et al. proposed evaluating ⟨ΨI | ∂R ΨJ ⟩ 1 ∂F = Ω J I ⟨ΨI | ∂R |ΨJ ⟩, where F is a generalized Kohn-Sham (KS) Fock operator and ΩJ I is the energy gap between state J and state I. In Refs. 29 and 23, we show that the Tavernelli formalism neither obeys the correct symmetries around a CI nor agrees with the exact Chernyak-Mukamel expression in the limit of infinite basis.30 A second approach is the direct differentiation of pseudowavefunctions that we offered in Ref. 23. This approach is identical to what Li and Liu have recently called the equationof-motion (EOM) TDDFT derivative coupling.22 While Li and Liu hypothesized this approach based on differentiating the RPA particle-hole operator, we began by guessing a TDDFT/time-dependent Hartree-Fock (TDHF) ground state wavefunction of the form  |Ψ0⟩ ≈ *1 + Xˆ I Yˆ I + |ΨDFT⟩. (1) I , Here, Xˆ I and Yˆ I are the excitation operators defined as  Xˆ I ≡ X iI a ai aa† , (2) ia

Yˆ I ≡



YiI a ai aa† .

(3)

ia

The derivative coupling vectors given by our pseudowavefunction ansatz recover the desired behaviors near a CI point and agree with the Chernyak-Mukamel equality in the limit of an infinite basis set near a CI point. 142, 064114-1

© 2015 AIP Publishing LLC

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Mon, 02 Mar 2015 11:18:02

064114-2

Ou et al.

Finally, the third approach is to calculate TDDFT derivative couplings via time-dependent response theory. To date, this is the only known approach that can provide exact couplings from TDDFT. An added advantage is the straightforward treatment of Pulay forces, which is essential when atom-centered basis sets are used. In July 2014, Li and Liu presented an abstract approach from time-dependent response theory for computing the excited states derivative couplings, which is applicable for TDDFT, but they did not present working equations or any numerical investigations of the methodology.22 Our goal in this work is to provide a detailed derivation of TDDFT/RPA derivative couplings from time-dependent response theory and to give working equations that can easily be implemented. (Note that while the present article was under review, Li and Liu have published a similar article exploring numerical examples that are in close agreement with the present manuscript.31) Moreover, we will also compare our response theory derivative couplings with those from a pseudo-wavefunction ansatz. (See also the article by Zhang and Herbert.32) An outline of this paper is as follows. In Sec. II, we present a self-consistent derivation of TDDFT/RPA derivative couplings from time-dependent response theory. In Sec. III, we present a numerical comparison of response theory results with our pseudo-wavefunction results for two cases: (a) two electronic states near a CI point and (b) two well-separated electronic states. In Sec. IV, we conclude. In Appendix A, we provide some necessary definitions, and in Appendix B, we demonstrate the equivalence of response theory and pseudowavefunction derivative couplings near a CI point. Unless otherwise specified, we use lowercase latin letters to denote spin molecular orbitals (MO) (a, b, c, d for virtual orbitals, i, j, k, l, m for occupied orbitals, p, q, r, s, w for arbitrary orbitals), greek letters (α, β, γ, δ, λ, σ, µ, ν) denote AOs. Many-electron excited states are denoted by Ψ (with uppercase latin indices I, J). We use atomic units and set ~ = 1. Q denotes nuclear coordinate and superscript [Q] denotes the differentiation with respect to Q.

J. Chem. Phys. 142, 064114 (2015)

II. THEORY

Fig. 1 gives a summary of how derivative couplings can be calculated with response theory. One calculates the second-order auxiliary coupling from exact many-body wavefunction quantum mechanics with response theory according to a sum-over-states (SOS), and then one calculates the matrix element with TDDFT. By comparing a given residue, one can extract the derivative coupling. One might wonder if the thus obtained couplings are physical, because the time-dependent Kohn-Sham (TDKS) system is fictitious. However, as long as the Chernyak-Mukamel expression for the coupling is recovered in the infinite basis set limit, the couplings converge to the exact result as better and better approximations are used for the exchange-correlation functional. A. Exact many-body wavefunction quantum mechanics (nothing to do with TDDFT) 1. Zeroth-, first-, and second-order response of the exact, many-body wavefunction to a time-dependent field according to perturbation theory

Consider an electronic system with the perturbed Hamiltonian H = H0 + λH1(t),  H1(t) = (V (α)eiω α t +V (α)∗e−iω α t ).

(4) (5)

α

The time-dependent Schrödinger equation is ∂ |Ψ(t)⟩. (6) ∂t Here, |Ψ(t)⟩ is the exact time-dependent wavefunction for the perturbed system (to the second-order) and can be expressed perturbatively as H |Ψ(t)⟩ = i

|Ψ(t)⟩ = |Ψ(0)(t)⟩ + λ|Ψ(1)(t)⟩ + λ2|Ψ(2)(t)⟩ + ···.

(7)

To construct |Ψ(t)⟩ in terms of {|ΨI ⟩}, the eigenstates of H0, we take the following steps. 1. The zeroth-order wavefunction is obtained by turning off the field (V (α) → 0) ∂ (0) |Ψ (t)⟩, ∂t |Ψ(0)(t)⟩ = e−i E0t |Ψ0⟩,

H0|Ψ(0)(t)⟩ = i

(8) (9)

where |Ψ0⟩ is the unperturbed ground state. 2. The first-order wavefunction is obtained by equating all terms linear in λ on both sides of Eq. (6), H1(t)|Ψ(0)(t)⟩ + H0|Ψ(1)(t)⟩ = i

FIG. 1. Summary for evaluating the TDDFT/RPA derivative coupling from time-dependent response theory.

∂ (1) |Ψ (t)⟩. ∂t

(10)

Here, the first-order wavefunction |Ψ(1)(t)⟩ does not contain the contribution from |Ψ0⟩ and can thus be expanded in the basis of eigenstates of H0 as follows:  |Ψ(1)(t)⟩ = e−i E0t b(1) (11) I (t)|ΨI ⟩, I ,0

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Mon, 02 Mar 2015 11:18:02

064114-3

Ou et al.

J. Chem. Phys. 142, 064114 (2015)

where the first-order coefficient b(1) I (t) has the form  b(1) W I(1)(ωα )eiω α t +W I(1)(−ωα )e−iω α t . I (t) =

∂ (0) ∂ (2) Ψ (t)⟩ + ⟨Ψ(0)(t)| Ψ (t)⟩ ∂Q ∂Q ∂ (1) + ⟨Ψ(1)(t)| Ψ (t)⟩ (21) ∂Q

C [Q],(2)(t) ≡ ⟨Ψ(2)(t)| (12)

α

Plugging Eqs. (9) and (11) into Eq. (10), one has  e−i E0t H1(t)|Ψ0⟩ + e−i E0t b(1) I (t)H0 |ΨI ⟩

≡ C1[Q],(2)(t) +C2[Q],(2)(t) +C3[Q],(2)(t) and we define

I ,0

=i

∂ −i E0t e b(1) I (t)|ΨI ⟩. ∂t I ,0 

(13)

C [Q],(2)(t) ≡

∂ (1) b (t) = 0, (14) ∂t I where Ω I = E I − E0 is the excitation energy for state I. Now substituting b(1) I (t) with the expression in Eq. (12) and collecting the coefficient of eiω α t , one obtains the expression for W I(1)(ωα ), ⟨ΨI |V (α)|Ψ0⟩ , Ω I + ωα

(15)

⟨ΨI |V (α)|Ψ0⟩ . Ω I − ωα

(16)

3. Second order: the second-order equation can be written as H1(t)|Ψ(1)(t)⟩ + H0|Ψ(2)(t)⟩ = i

∂ (2) |Ψ (t)⟩, ∂t

(17)

where |Ψ(2)(t)⟩ is the second-order wavefunction with the expansion coefficient b(2) I (t),  |Ψ(2)(t)⟩ = e−i E0t b(2) (18) I (t)|ΨI ⟩,

+ C˜ [Q],(2)(−ωα , −ω β )ei(−ω α −ω β )t .

b(2) I (t) =

αβ

∂ −i E0t  (1) ∂ |Ψ(1)(t)⟩ = e bI (t)|ΨI ⟩ ∂Q ∂Q I ,0 = − iE0[Q]t|Ψ(1)(t)⟩ + e−i E0t



b(1),[Q] (t)|ΨJ ⟩ J

J ,0

+e

−i E 0t



[Q] b(1) J (t)|ΨJ ⟩.

(25)

C3[Q],(2)(t) can be obtained after multiplying by ⟨Ψ(1)(t)|,

+ W I(2)(ωα ,−ω β )ei(ω α −ω β )t

C3[Q],(2)(t) = ⟨Ψ(1)(t)| −iE0[Q]t|Ψ(1)(t)⟩  + e−i E0t ⟨Ψ(1)(t)|b(1),[Q] (t)|ΨJ ⟩ J

+ W I(2)(−ωα , ω β )ei(−ω α +ω β )t

I ,0

(19)

Similar to the first-order case, the following expression for W I(2)(ωα , ω β ) can be obtained by comparing the coefficients of ei(ω α +ω β )t : 1 W I(2)(ωα , ω β ) = Ω I + ωα + ω β   ⟨ΨJ |V (β)|Ψ0⟩ ⟨ΨI |V (β)|ΨJ ⟩ × ΩJ + ω β J ,0  ⟨ΨJ |V (α)|Ψ0⟩ ⟨ΨI |V (α)|ΨJ ⟩ + . Ω J + ωα

(24)

J ,0

W I(2)(ωα , ω β )ei(ω α +ω β )t

 + W I(2)(−ωα ,−ω β )ei(−ω α −ω β )t .

(23)

Let us now show that, in most circumstances, the derivative coupling between state I and state J can be found by evaluating the residue of C˜ [Q],(2)(ωα , ω β ) at poles ωα = Ω I and ω β = −ΩJ . To prove this statement, note that the first two terms in Eq. (21), C1[Q],(2)(t) and C2[Q],(2)(t), have residues at ωα = ±ΩJ and ωα + ω β = ±ΩJ ; as such, these terms are not expected to contribute to the pole of C˜ [Q],(2)(ωα , ω β ). Therefore to isolate the derivative coupling, we will now focus on C3[Q],(2)(t), which involves only the first-order wavefunctions. According to Eq. (11), the nuclear derivative of the first-order wavefunction can be expressed as

I ,0



C˜ [Q],(2)(ωα , ω β )ei(ω α +ω β )t

+ C˜ [Q],(2)(ωα , −ω β )ei(ω α −ω β )t + C˜ [Q],(2)(−ωα , ω β )ei(−ω α +ω β )t

⟨ΨI |H1(t)|Ψ0⟩ + Ω I b(1) I (t) −i

W I(1)(−ωα ) = −

 αβ

Left-multiplying by ⟨ΨI | gives (after relabeling the indices)

W I(1)(ωα ) = −

(22)

+e

−i E 0t

 [Q] ⟨Ψ(1)(t)|b(1) J (t)|ΨJ ⟩.

(26)

I ,0

Inserting the SOS representation for the first-order wavefunction, one has  ( ) (1),[Q] C3[Q],(2)(t) = −iE0[Q]tb(1) (t)e−i E0t J (t) + b I I, J ,0 × b(1)∗ I (t)⟨ΨI |ΨJ ⟩ +



(1) [Q] b(1)∗ I (t)bJ (t)⟨ΨI |ΨJ ⟩.

I, J ,0

(27)

(20)

Similar expression exists for W I(2)(ωα , −ω β ), etc.

According to the orthogonality of {|ΨI ⟩}, C3[Q],(2)(t) =

2. Second-order auxiliary coupling

The second-order auxiliary coupling matrix element is defined as



(1),[Q] [Q] (1) * b(1)∗ (t)e−i E0t I (t) −iE0 tb I (t) + b I I ,0 ,  [Q] + + b(1) (28) J (t)⟨ΨI |ΨJ ⟩ . J ,0 -

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Mon, 02 Mar 2015 11:18:02

064114-4

Ou et al.

J. Chem. Phys. 142, 064114 (2015)

When evaluating the residue of C˜3[Q],(2)(ωα , ω β ) at poles ωα = Ω I and ω β = −ΩJ , one finds that only the last term in C3[Q],(2)(t) (1) ˜ [Q],(2)(ωα , ω β ) can be written as (omitting the contributes. Plugging in the expression of b(1)∗ I (t) and bJ (t), we find that C3 non-contributing part) C˜3[Q],(2)(ωα , ω β ) =

 ⟨Ψ0|V (α)(ωα )|ΨI ⟩ ⟨ΨJ |V (β)(ω β )|Ψ0⟩ ⟨ΨI |ΨJ[Q]⟩ + ···. (Ω − ω )(Ω + ω ) I α J β I, J ,0

Thus, the residue of C˜3[Q],(2)(ωα , ω β ) at poles ωα = Ω I and ω β = −ΩJ is therefore Res[C˜3[Q],(2)(ωα ,ω β );Ω I ,−ΩJ ] = ⟨Ψ0|V (α)(ωα )|ΨI ⟩ ⟨ΨJ |V (β)(ω β )|Ψ0⟩ ⟨ΨI |ΨJ[Q]⟩

(30)

≡ V0I VJ 0⟨ΨI |ΨJ[Q]⟩.

(31)

Note that the factor 2πi is omitted when evaluating all the residues. The final residue of the second-order auxiliary couplings is then just the residue of C˜3[Q],(2)(ωα , ω β ), Res[C˜ [Q],(2)(ωα ,ω β );Ω I ,−ΩJ ] = Res[C˜3[Q],(2)(ωα ,ω β );Ω I ,−ΩJ ] = V0I VJ 0⟨ΨI |ΨJ[Q]⟩.

(32)

Up to this point, all of the above theory is simple timedependent perturbation theory, and we have discussed nothing having to do with TDDFT. B. TDDFT

As shown in Sec. II A 2, a derivative coupling is related to the residue of the second-order auxiliary coupling at certain frequencies according to perturbative time-dependent quantum mechanics. In this section, we will derive an explicit expression for such a residue via TDDFT. 1. First- and second-order TDKS orbitals

According to the usual TDDFT framework, the general eigenvalue equation for a closed electronic system is F|ΨDFT⟩ = EDFT|ΨDFT⟩,

where |ΨDFT⟩ is the non-interacting Kohn-Sham ground state with energy EDFT. Here, |ΨDFT⟩ is |ΨDFT⟩ = |φ1φ2 ...φ n ⟩,

(34)

where |φi ⟩ is the ith non-interacting KS orbital with energy ε i . The Fock operator is 

(35) Fpq = h pq + φ p φ s |φq φr γr(0)s , rs

γr(0)s

where is the time-independent density matrix





γr(0)s = ΨDFT|ar† a s |ΨDFT = φr |φ j φ j |φ s = δrocc s .

(denoted by | φ˜i ⟩) can be expanded as (to second-order) ( ) (2) | φ˜i (t)⟩ = e−iε i t |φi ⟩ + |φ(1) (t)⟩ + |φ (t)⟩ , (37) i i (2) where |φ(1) i (t)⟩ and |φ i (t)⟩ are first- and second-order orbital corrections, respectively. The TDKS density matrix is  γ(t) ˜ = | φ˜i (t)⟩ ⟨φ˜i (t)| (38) i

≡ γ (0) + γ˜ (1)(t) + γ˜ (2)(t),

(39)

where γ˜ (1)(t) =

(

) (1) |φi ⟩ ⟨φ(1) i (t)| + |φ i (t)⟩ ⟨φ i | ,

(40)

i

γ˜ (2)(t) =

(

(2) |φ(2) i (t)⟩ ⟨φ i | + |φ i ⟩ ⟨φ i (t)|

i

) (1) (1) (1) + |φ(1) i (t)⟩ ⟨φ i (t)| + |φ i (t)⟩ ⟨φ i (t)| .

(41)

At this point, we want to express the first- and second-order orbital corrections in terms of first- and second-order density matrices, γ˜ (1)(t) and γ˜ (2)(t). We perform this transformation because the density matrices are the central objects in most TDDFT development.33 Let us define  γ˜ (1)(t) ≡ γ (1)(ωα )eiω α t + γ (1)(−ωα )e−iω α t . (42) α

In the frequency domain, the first-order density matrix response is  γ (1)(ω) = γ (1)(ωα )δ(ω − ωα ), (43) α

γ (ωα ) = (1)

(33)

(29)



X˜ ia (ωα )|φa ⟩ ⟨φi | + Y˜ia (ωα )|φi ⟩ ⟨φa |, (44)

ia

where X˜ ia (ωα ) and Y˜ia (ωα ) are the virt-occ and occ-virt matrix elements of γ (1)(ωα ), respectively.33 If we substitute Eq. (40) for γ˜ (1)(t) and then sandwich everything by ⟨φ j | and |φa ⟩, we find   ⟨φ(1) (t)|φ ⟩ = Y˜ia (ωα )eiω α t + Y˜ia (−ωα )e−iω α t . (45) a i aα

Given that X˜ ia (ωα ) = Y˜ia∗ (−ωα ),34 one arrives at the final expression for the first-order TDKS orbital correction,   (t)⟩ = (46) |φ(1) X˜ ia (ωα )eiω α t + X˜ ia (−ωα )e−iω α t |φa ⟩. i aα

(36)

j

Now, when a time-dependent field is applied (as in Eq. (6)), the system is perturbed and the time-dependent KS orbitals

Our next step is to deal with the second-order density matrix response γ (2) in order to get an expression for the second-order orbital correction. The time-dependence of γ (2) is given by

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Mon, 02 Mar 2015 11:18:02

064114-5

Ou et al.

J. Chem. Phys. 142, 064114 (2015)



γ˜ (2)(t) ≡

γ (2)(ωα , ω β )ei(ω α +ω β )t

αβ

+ γ (2)(ωα ,−ω β )ei(ω α −ω β )t + γ (2)(−ωα , ω β )ei(−ω α +ω β )t + γ (2)(−ωα ,−ω β )ei(−ω α −ω β )t .

(47)

Sandwiching Eq. (41) by occ and virt orbitals, one finds the different components of γ˜ (2)(t): 1. occ-occ (2) γ˜ i(2) j (t) = 2⟨φ i |φ j (t)⟩;

(48)

(2) γ˜ ia (t) = ⟨φ(2) i (t)|φ a ⟩;

(49)

2. occ-virt

Note that the derivative of the exponential part e−iε i t of φ˜i (t) (as shown in Eq. (37)) is not included in the above expression since that term will have no effect on the final result (which [Q],(2) is easy to prove). We will now treat the three terms in CKS separately. [Q]   [Q] 1. T1: j ⟨φ(2) ⟩ and T2: j ⟨φ j |φ(2) ⟩ j (t)|φ j j (t) To evaluate these terms, we plug in the expressions for |φ(2) j (t)⟩ (Eq. (52)) and T1 becomes  1  (2)∗ R[Q] γ˜ i j (t)OiR[Q] + γ˜ a(2)∗ j j (t)Oa j , 2 ij aj

T1 =

where we define the “right” spin-orbital derivative overlap O R[Q] pq the same way as in Ref. 23, [Q] O R[Q] pq ≡ ⟨φ p |φq ⟩

3. virt-occ (2) γ˜ ai (t) = ⟨φa |φ(2) i (t)⟩;

    [Q] = *. Cµ p ⟨µ|+/ * |ν⟩Cνq + ν [Q] Cνq + (58) ν , µ -, ν  

(2) (1) (1) (1) γ˜ ab (t) = ⟨φa |φ(1) i (t)⟩ ⟨φ i (t)|φ b ⟩ + ⟨φ a |φ i (t)⟩ ⟨φ i (t)|φ b ⟩.

From Eq. (51), it is clear that the second-order orbital (2) correction (|φ(2) ˜ ab (t). i (t)⟩) has no explicit dependence on γ According to Eqs. (48) and (49) (or (50)), the expression for |φ(2) i (t)⟩ is  1  (2) (2) |φ(2) (t)⟩ = γ ˜ (t)|φ ⟩ + γ˜ ai (t)|φa ⟩ (52) j i 2 j ij a 1   (2) = γ (ωα , ω β )ei(ω α +ω β )t 2 jα β i j i(ω α −ω β )t + γi(2) j (ωα ,−ω β )e i(−ω α +ω β )t + γi(2) j (−ωα , ω β )e

 i(−ω α −ω β )t + γi(2) |φ j ⟩ j (−ωα ,−ω β )e  (2) i(ω α +ω β )t γai (ωα , ω β )e +

(53)

j



  aγ





(59)

) 1 [Q] Sµν Cνq − Θ[Q] pq 2

(60)

µν

( R[Q] Cµ p Sµν −

µν

A[Q] Cµ p Sµν Cνq − Θ[Q] pq .

(61)

µν

Here, Cµ p denotes the MO coefficients, Θ pq stands for [Q] ≡ ⟨µ|ν⟩[Q] are meaningful an orbital rotation matrix, Sµν R[Q] overlap derivatives contributing to Pulay forces, and Sµν A[Q] R[Q] 1 [Q] 1 [Q] [Q] [Q] ≡ ⟨µ|ν ⟩. Sµν ≡ Sµν − 2 Sµν = 2 (⟨µ|ν ⟩ − ⟨µ |ν⟩) are artificial matrix elements that can be ignored if we introduce electron-translation factors.35 (The detailed derivation of O R[Q] pq can be found in Ref. 35.) The (ωα +ω β ) Fourier coefficient of T1 is  1  (2) T 1= γi j (ωα , ω β )OiR[Q] + γa(2)j (ωα , ω β )OaR[Q] j j , 2 ij aj

) (54) (55)

γ (2) j a (ωα , ω β )

where and correspond to the (2) (ωα + ω β ) Fourier transforms of γ (2) j i (t) and γ j a (t), respectively. Similarly, the (ωα + ω β ) Fourier coefficient for T2 is 1  (2) 1  (2),[Q] γj j (ωα , ω β ) + γ (ωα , ω β )O R[Q] T2= ji 2 j 2 i j ji  R[Q] + γ (2) (63) j a (ωα , ω β )O j a . aj

2. T3:

≡ T1 +T2 +T3.

j

=



R[Q] Cµ p Sµν Cνq

(62)

We are now prepared to evaluate the second-order auxiliary coupling using the perturbed TDKS orbitals ( [Q] [Q],(2) [Q] CKS = ⟨φ(2) ⟩ + ⟨φ j |φ(2) ⟩ j (t)|φ j j (t)

T3 =

µν

γ (2) j i (ωα , ω β )

2. Auxiliary coupling matrix elements for the TDKS determinant

[Q]

[Q] Cµ p Sµν Cνq +

=

(51)

(1) + ⟨φ(1) j (t)|φ j (t)

(57)

(50)

4. virt-virt

aα β (2) + γai (ωα ,−ω β )ei(ω α −ω β )t (2) + γai (−ωα , ω β )ei(−ω α +ω β )t  (2) + γai (−ωα ,−ω β )ei(−ω α −ω β )t |φa ⟩.

(56)



(1) (1) [Q] ⟩ j ⟨φ j (t)|φ j (t)

To evaluated this term, we plug in the expression for the first-order orbital corrections (Eq. (46)) into T3. One finds

∗   ∂  ˜ X˜ j a (ωγ )eiωγ t + X˜ j a (−ωγ )e−iωγ t ⟨φa | X j b (ωδ )eiω δ t + X˜ j b (−ωδ )e−iω δ t φb ⟩ . ∂Q bδ

(64)

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Mon, 02 Mar 2015 11:18:02

064114-6

Ou et al.

J. Chem. Phys. 142, 064114 (2015)

Recall that X˜ ia (ω) = Y˜ia∗ (−ω), and it is easy to write down the (ωα + ω β ) Fourier transform of T3,   ( ) Y˜j a (ωα ) X˜ [Q](ω β ) + Y˜j a (ωα ) X˜ j b (ω β )O R[Q] + Y˜j a (ω β ) X˜ [Q](ωα ) + Y˜j a (ω β ) X˜ j b (ωα )O R[Q] . T3= ja

ab

ja

b

T 1+T 2+T 3   R[Q] 1 γi(2) (ω , ω )O + γ (2),[Q] (ωα , ω β ) = α β j ij jj 2 ij j   R[Q] + Y˜j a (ωα ) X˜ j b (ω β ) + Y˜j a (ω β ) X˜ j b (ωα ) Oab j ab

+ +

(

) R[Q] γa(2)j (ωα , ω β ) − γ (2) j a (ωα , ω β ) Oa j

ja

) ˜ ˜ [Q] ˜ X˜ j[Q] a (ωα )Yj a (ω β ) + X j a (ω β )Yj a (ωα ) .

(66)

ja

As a result of the idempotency of the Kohn-Sham density matrix,33 the virt-virt and occ-occ pieces of the second(2) order density matrix γab (ωα , ω β ) satisfy  (2) γab (ωα , ω β ) = ( X˜ j a (ωα )Y˜j b (ω β ) + X˜ j a (ω β )Y˜j b (ωα )), j

γi(2) j (ωα , ω β )

=−



(67) ˜ ˜ ˜ ˜ ( X j a (ωα )Yia (ω β ) + X j a (ω β )Yia (ωα )).

a

(68) If we apply these expressions to Eq. (66), the total Fourier coefficient becomes  R[Q] T 1+T 2+T 3 = γi(2) j (ωα , ω β )Oi j ij

+ +



(2) R[Q] γab (ωα , ω β )Oab

ab  (

) R[Q] γa(2)j (ωα , ω β ) − γ (2) j a (ωα , ω β ) Oa j

ja

1   ˜ [Q] + X j a (ωα )Y˜j a (ω β ) 2 ja ˜ ˜ ˜ [Q] + X˜ j[Q] a (ω β )Yj a (ωα ) − X j a (ωα )Yj a (ω β ) − X˜ j a (ω β )Y˜j[Q] a (ωα )



(69)

≡ O + XY,

(70)

where we define   (2) R[Q] R[Q] O≡ γi(2) + γab (ωα , ω β )Oab j (ωα , ω β )Oi j ij

+

ab

(

γa(2)j (ωα , ω β ) − γ (2) j a (ωα , ω β )

)

OaR[Q] j ,

(71)

ja

XY ≡

ab

(65)

b

The total (ωα +ω β ) Fourier coefficient of the second-order auxiliary coupling given by TDDFT can thus be expressed as

(

ja

1   ˜ [Q] ˜ X j a (ωα )Y˜j a (ω β ) + X˜ j[Q] a (ω β )Yj a (ωα ) 2 ja  ˜ ˜ [Q] − X˜ j a (ωα )Y˜j[Q] a (ω β ) − X j a (ω β )Yj a (ωα ) .

3. Residues of the coefficient derivative terms

We must now evaluate the residues of the O and XY terms. 1. Residue of XY ˜ α ) and Y˜ (ωα ) We begin with the residue of the X(ω derivatives. Start with the general working equation for RPA,36  A B  ˜ (α) * + + ωα *1 0 + * X(ωα )+ = − *V + . (73) (α) ˜ ,B A,0 −1- ,Y (ωα ) ,V ˜ α ) and Y˜ (ωα ) can be obtained by We note that X(ω   −1 (α) ˜ * X(ωα )+ = − *A B+ + ωα *1 0 + *V + . (74) (α) ˜ ,B A,Y (ωα ) ,0 −1- ,V According to Ref. 33, the inverse of the supermatrix is equivalent to  A B  −1   ( ) 1 * X I + X I YI * + + ωα *1 0 + = Ω I + ωα , YI ,B AI ,0 −1- ( ) 1 * YI + YI X I . + Ω I − ωα , X I (75) Therefore, if V (α) = µ(α) (the dipole operator), X˜ and Y˜ can be expressed as  X jI a µ(α) YjI a µ(α) 0I I0 + *. /, + (76) Ω + ω Ω − ω I α I α I , I a (α)  YjI a µ(α) X µ j I0 + 0I *. /. Y˜j a (ωα ) = + (77) Ω + ω Ω − ω α I α I , I Thus, the residues of X˜ and Y˜ at poles ωα = Ω I and ω β = −ΩJ are X˜ j a (ωα ) =

Res[ X˜ j a (ωα );Ω I ] = YjI a µ(α) I0 ,

(78)

X jI a µ(α) I0 ,

(79)

Res[ X˜ j a (ωα );−ΩJ ] = X jJ a µ(α) 0J ,

(80)

Res[Y˜j a (ωα );−ΩJ ] = YjJ a µ(α) 0J ,

(81)

Res[Y˜j a (ωα );Ω I ] =

and the derivatives of X˜ and Y˜ are I a (α),[Q]   X jI a[Q] µ(α) 0I + X j µ0I  [Q] ˜  X j a (ωα ) = Ω I + ωα I   I a (α),[Q] YjI a[Q] µ(α) I 0 +Yj µ I 0 + Ω I − ωα [Q] X jI a µ(α) 0I Ω I

(72)

[Q]   YjI a µ(α) I 0 ΩI  , − − 2 2 (Ω I + ωα ) (Ω I − ωα )  

(82)

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Mon, 02 Mar 2015 11:18:02

064114-7

Ou et al.

Y˜j[Q] a (ωα ) =

J. Chem. Phys. 142, 064114 (2015) I a (α),[Q]   YjI a[Q] µ(α) 0I +Yj µ0I   Ω I + ωα I   I a[Q] (α) Xj µ I 0 + X jI a µ(α),[Q] I0 + Ω I − ωα



[Q] YjI a µ(α) 0I Ω I (Ω I + ωα )2



=−

(

Res



 (2) R[Q] γab (ωα , ω β )Oab ;Ω I ,−ΩJ

ab

= (83)

  (α)a Ia Ia where µ(α) . We can now write down j a (X j +Yj )Vj 0I = the residue of XY at the poles ωα = Ω I and ω β = −ΩJ , Res[XY;Ω I ,−ΩJ ]  ) 1  ( I a[Q] (α) = Yj µ I 0 +YjI a µ(α),[Q] YjJ a µ(β) I 0 0J 2 ja ( ) J a (β),[Q] + X jJ a[Q] µ(β) X jI a µ(α) 0J + X j µ0J I0 ( ) J a[Q] (β) I a (α) J a (β),[Q] − Yj µ I 0 Yj µ0J +Yj µ0J ( I a[Q] (α) ) I a (α),[Q] J a (β) µI 0 + X j µI 0 − X j µ0J X j  ) 1  ( Ia Ja = X j X j −YjI aYjJ a 2 ja ( ) (β),[Q] × µ(α) − µ(α),[Q] µ(β) I 0 µ0J I0 0J ( + X jI a X jJ a[Q] −YjI aYjJ a[Q] − X jI a[Q] X jJ a ) (α) (β) I a[Q] J a + Yj Yj µ I 0 µ0J . Recalling the orthogonality condition ( ) X jI a X jJ a −YjI aYjJ a = δ I J ,

(84)

(85)

(86)

ja

we find that the first term in Eq. (85) vanishes. Taking the derivative on each side of Eq. (86), one finds (for I , J) ( X jI a[Q] X jJ a −YjI a[Q]YjJ a

(

) R[Q] (α) (β) X jI b X jJ a +YjI aYjJ b Oab µ I 0 µ0J .

= 0,

(87)

(91)

ab j (2) (2) For the virt-occ (γvo (ωα , ω β )) and occ-virt (γov (ωα ,ω β )) pieces of the second-order density matrix, these matrix elements must satisfy33  A B  (2) * + + (ωα + ω β ) * I 0 + *γvo (ωα , ω β ) + (2) ,B A,0 −I- , γov (ωα , ω β ) L˜ vo(ωα , ω β )+ = − * ov . (92) ˜ , L (ωα , ω β )See Appendix B 1 for explicit expressions of L˜ vo(ωα , ω β ) and L˜ ov(ωα , ω β ). Note that L˜ vo(ωα ,ω β ) and L˜ ov(ωα ,ω β ) are bilinear in X˜ j a (ωα ) and Y˜j a (ω β )33 (or some permutations thereof). (2) (2) The residue of γvo (ωα ,ω β ) and γov (ωα , ω β ) at ωα = Ω I and ω β = −ΩJ are γ IvoJ and γ IovJ , with

  −1 vo vo *γ I J + = − *A B+ + ΩJ I * I 0 + * L I J + . (93) ov ov ,B A,γ I J ,0 −I- , L I J ov L vo I J and L I J are the virt-occ and occ-virt Lagrangians, vo ˜ ov L I J = Res[ L˜ vo(ωα , ω β );Ω I ,−ΩJ ] and L ov I J = Res[ L (ωα , ω β );Ω I ,−ΩJ ] (see Appendix B 1 for explicit expressions). Thus, the overall residue of O can be written as Res[O;Ω I ,−ΩJ ]  ( ) X iI a X jJ a +YjI aYiJ a OiR[Q] = − j i ja

+

(

+



) R[Q] X jI b X jJ a +YjI aYjJ b Oab

ab j

ja

+ X jI a X jJ a[Q] −YjI aYj

(90)

i ja

[Q]   X jI a µ(α) I 0 ΩI  , 2 (Ω I − ωα ) 

) J a[Q]

) (β) X iI a X jJ a +YjI aYiJ a OiR[Q] µ(α) j I 0 µ0J ,

  γ IvoJ − γ IovJ a j OaR[Q] j

(β) µ(α) I 0 µ0J .

(94)

ja

(

−X jI a[Q] X jJ a +YjI a[Q]YjJ a

)

ja

=

(

) X jI a X jJ a[Q] −YjI aYjJ a[Q] .

(88)

ja

Plugging Eq. (87) into Eq. (85) and taking out the first term, we get the final result for the residue of XY, Res[XY;Ω I ,−ΩJ ] ( ) (β) = X jI a X jJ a[Q] −YjI aYjJ a[Q] µ(α) I 0 µ0J .

(89)

ja

2. Residue of O The residues of the occ-occ and virt-virt part in O are straightforward. Using Eqs. (68) and (78)–(81), the residues of the occ-occ and virt-virt parts at ωα = Ω I and ω β = −ΩJ are just   R[Q] Res γi(2) (ω , ω )O ;Ω ,−Ω α β I J j ij ij

C. Final expression for TDDFT derivative couplings

Thus, the total residue for the (ωα +ω β ) Fourier transform of the auxiliary couplings obtained from TDDFT is Res[O + XY;Ω I ,−ΩJ ]  ( ) = X jI a X jJ a[Q] −YjI aYjJ a[Q] ja



(

+

(

+



) X iI a X jJ a +YjI aYiJ a OiR[Q] j

i ja

) R[Q] X jI b X jJ a +YjI aYjJ b Oab

ab j

γ IvoJ − γ IovJ

 aj

 (β) OaR[Q] µ(α) j I 0 µ0J .

(95)

ja

Comparing Eqs. (31) and (95), we can write down the Response Theory derivative coupling between state I and

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Mon, 02 Mar 2015 11:18:02

064114-8

Ou et al.

J. Chem. Phys. 142, 064114 (2015)

state J (dRT I J ),



⟨ΨI |ΨJ[Q]⟩

=

(

X iI a X iJ a[Q] −YiI aYiJ a[Q]

)

− +

(

+



X iI a X jJ a +YjI aYiJ a

)

) R[Q] X jI b X jJ a +YjI aYjJ b Oab

ab j

 aj

OaR[Q] j

(96)

ja



dRT IJ.

(97)

Comparing dRT I J with our derivative coupling based on a pseudo-wavefunction ansatz (dPW I J , Eqs. (18) and (24) in Ref. 23), we find that the difference between two approaches is no more than   PW (98) γ IvoJ − γ IovJ a j OaR[Q] dRT I J = dI J + j . ja

This result agrees with the recent work by Li and Liu.22 Plugging in the final expression for dPW I J (Eq. (47) in Ref. 23) and the expression of OaR[Q] given by Eq. (61), one can write j down dRT in a finite AO basis as IJ ⟨ΨI |ΨJ[Q]⟩   XI XJ YI YJ 1  I J ˜ [Q]  * Rµλ Rσν + Rµλ Rσν + ˜ [Q] D µν h µν + = G µνλσ  XI YJ XJ ΩJ I µν +Rµλ Rνσ + RYµλI Rνσ µνλσ ,   ( ) [Q] ˜ IJ IJ IJ  − 1 + D µλ Pσν Π [Q] S P D + D Fα β µα µν βν ν β µνλσ  2 α β µν  R X I R X J + RY I RY J *. νγX I δ βX J νγY I δ βY J +/ 1  [Q] ˜ ..+Rδ β Rνγ + Rδ β Rνγ // Sµν Pµα . X I Y J − G Y I X J / α βγδ 2 µνα βγδ ..+Rνγ Rβδ + Rνγ Rβδ // XI YJ Y I XJ ,+Rβδ Rνγ + Rβδ Rνγ XI XJ YI YJ Rγν Rβδ + Rγν Rβδ

*. X I X J + Y I Y J/ 1  [Q] ˜ ..+Rβδ Rγν + Rβδ Rγν // G − Sµν Pµα . X I Y J Y I X J / α βγδ 2 µνα βγδ ..+Rγν Rδ β + Rγν Rδ β // XI YJ Y I XJ ,+Rδ β Rγν + Rδ β Rγν R X I R X J + RY I RY J *. αγX I δ βX J αγY I δ βY J +/  .+Rδ β Rαγ + Rδ β Rαγ // Ξα βγδλσ [Q] ˜ Pµλ Pνσ .. X I Y J − Sµν Y I X J/ 2 ..+Rαγ Rβδ + Rαγ Rβδ // µνα βγδλσ XI YJ Y I XJ ,+Rβδ Rαγ + Rβδ Rαγ  ) ( 1  [Q] ˜ IJ IJ + Sµν Pµα Pνδ Dγ β + D βγ Gα βγδ 2 µνα βγδ   A[Q] − X iI a X iJ b −YiI aYiJ b CνaCµb Sµν µνiab



 (

) A[Q] X iI a X jJ a −YiI aYjJ a Cνi Cµ j Sµν

µνi j a

+

 ja

γ IvoJ − γ IovJ

  aj

µν

A[Q] CµaCν j Sµν

 aj

+

 1 L a j ΘR[Q] aj . ΩJ I

(99)

After simplifying the orbital response (i.e., those terms involving ΘR[Q] in the last line of Eq. (99); see Appendix aj B 2 for details), and invoking the “z-vector” trick,37 the final TDDFT/RPA derivative coupling between state I and state J given by response theory is then

OiR[Q] j

i ja

γ IvoJ − γ IovJ

γ IvoJ − γ IovJ

ja

ia

(



⟨ΨI |ΨJ[Q]⟩  YI YJ  XI XJ 1  I J ˜ [Q]  * Rµλ Rσν + Rµλ Rσν + ˜ [Q] D µν h µν + G µνλσ =  XI YJ XJ ΩJ I µν +Rµλ Rνσ + RYµλI Rνσ µνλσ ,   ( ) [Q] ˜ IJ IJ − 1 + D µλ Pσν Π [Q] Sµν Pµα D βν + DνI βJ Fα β µνλσ   2 α β µν R X I R X J + RY I RY J *. νγX I δ βX J νγY I δ βY J +/ 1  [Q] ˜ ..+Rδ β Rνγ + Rδ β Rνγ // Sµν Pµα . X I Y J − Gα βγδ Y I X J/ 2 µνα βγδ ..+Rνγ Rβδ + Rνγ Rβδ // XI YJ Y I XJ ,+Rβδ Rνγ + Rβδ Rνγ R X I R X J + RY I RY J *. γνX I βδX J γνY I βδY J +/ 1  [Q] ˜ ..+Rβδ Rγν + Rβδ Rγν // Gα βγδ Sµν Pµα . X I Y J − Y I X J/ 2 µνα βγδ ..+Rγν Rδ β + Rγν Rδ β // XI YJ Y I XJ ,+Rδ β Rγν + Rδ β Rγν R X I R X J + RY I RY J *. αγX I δ βX J αγY I δ βY J +/  .+Rδ β Rαγ + Rδ β Rαγ // Ξα βγδλσ [Q] ˜ − Sµν Pµλ Pνσ .. X I Y J Y I X J/ 2 ..+Rαγ Rβδ + Rαγ Rβδ // µνα βγδλσ XI YJ Y I XJ ,+Rβδ Rαγ + Rβδ Rαγ  ( ) 1  [Q] ˜ IJ Sµν Pµα Pνδ DγI Jβ + D βγ Gα βγδ + 2 µνα βγδ   A[Q] − X iI a X iJ b −YiI aYiJ b CνaCµb Sµν µνiab



 (

) A[Q] X iI a X jJ a −YiI aYjJ a Cνi Cµ j Sµν

µνi j a

+

 ja

γ IvoJ − γ IovJ

  aj

A[Q] CµaCν j Sµν

µν

1  vo [Q] − (γ + γ IovJ )bi Mbi . 2ΩJ I bi I J

(100)

All terms in Eqs. (99) and (100) are defined in Appendix A. This final expression for dRT I J is very similar to the one for dPW except for the orbital response terms and the (rather IJ A[Q] meaningless) Sµν terms. D. Spurious poles of the adiabatic TDDFT derivative couplings

Despite the appeal of response theory, a crucial aspect of Eq. (100) is unphysical. Consider again Eq. (92), where we solve for the relaxed part of the second-order density matrix (2) (2) response, γvo (ωα , ω β ) and γov (ωα , ω β ). Because L˜ vo(ωα , ω β )

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Mon, 02 Mar 2015 11:18:02

064114-9

Ou et al.

FIG. 2. (a) The norm ratio

J. Chem. Phys. 142, 064114 (2015)

∥dRT IJ∥ ∥dPW IJ ∥

and (b) the angle

PW dRT I J ·d I J PW ∥dRT I J ∥ ∥d I J ∥

between CH2NH+2 S1/S2 derivative coupling vectors given by a pseudo-wavefunction ansatz and

response theory at various distances to the CI point.

and L˜ ov(ωα , ω β ) are bilinear in X˜ j a (ωα ) and Y˜j a (ω β ) (or some (2) permutations thereof), the pole structure of γvo (ωα , ω β ) (or (2) γov (ωα ,ω β )) is of the form γa(2)j (ωα , ω β ) =

   X jK a X iK b L˜ vo (ωα , ω β ) +YiK b L˜ ov (ωα , ω β ) bi bi K bi

Ω K + ωα + ω β

+ ··· (101)

=



  X jK a X iK b L vo +YiK b L ov bi bi

I J K bi

(Ω I − ωα )(ΩJ + ω β )(Ω K + ωα + ω β )

+ ···, (102)

ov ˜ vo where, again, L vo I J = Res[ L (ωα , ω β );Ω I ,−Ω J ] and L I J ov ˜ = Res[ L (ωα ,ω β );Ω I ,−ΩJ ] (see Appendix B 1 for explicit (2) expressions). Thus, for adiabatic TDDFT, γvo (ωα , ω β ) (or (2) γov (ωα ,ω β )) contains products of three poles, which is at variance with the pole structure of the exact time-dependent density response—the latter containing at most products of two poles. The existence of such spurious poles appears to have been recognized previously by Dalgaard38 in the case of TDHF response theory. We can expect these terms to be problematic when ±Ω K + Ω I − ΩJ = 0, where Ω K is any excitation energy; at such a geometry, Eq. (92) cannot be inverted and the interstate transition density matrix will diverge. To see this effect, a numerical example will be given in Sec. III B. Note that Li and Liu independently came to the exact same conclusion in their concurrent paper.31

derivative coupling of protonated formaldimine (CH2NH+2 ) near its S1/S2 CI point and we showed that the derivative coupling vectors obtained from the direct differentiating our pseudo-wavefunction ansatz recover all the desired properties along a loop around the CI point: (i) the derivative couplings lie rigorously on the branching plane and are perpendicular to the energy difference gradient direction; (ii) their magnitudes are identical everywhere on the loop in the proper coordinate system; (iii) the path integral gives the Berry’s phase. In this work, we will reexamine this molecule as an example to show that the derivative couplings for TDDFT/RPA response theory are identical to those for a pseudo-wavefunction ansatz near the CI point. We numerically implemented Eq. (100) in Q-Chem39,40 PW and ωB97X/6-31G∗∗ is used to get both dRT I J and d I J between + S1 and S2 of CH2NH2 . Taking the g direction as an example, we gradually increased the distance moved away from the CI point and computed both derivative couplings at various distances. Then, we compared the norm and the direction for

III. NUMERICAL EXAMPLES A. S1/ S2 conical intersection of protonated formaldimine

To investigate the derivative couplings (dRT I J ) derived above, we now study two sample cases: (a) a CI when ΩJ I → 0 and (b) the case of two well separated electronic states (ΩJ I ≫ 0). In a previous study,23 we explored the S1/S2

FIG. 3. Relative energies for S0, S1, and S2 states of CH2NH+2 versus the distance moved from the CI point.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Mon, 02 Mar 2015 11:18:02

064114-10

Ou et al.

J. Chem. Phys. 142, 064114 (2015)

FIG. 4. S1/S4 TDHF derivative coupling elements for LiH obtained from a pseudo-wavefunction ansatz (PWDC) and response theory (RTDC) as a function of bond length. In (a), we use a linear scale, and in (b), we use a logarithm scale. Note that the two methods basically agree except at 1.9 Å, where RTDC diverges. The absolute value of the derivative coupling with respect to Li z coordinate is plotted.

two derivative coupling vectors computed at every distance. As shown in Fig. 2, the agreement between the two approaches is very good up to a large distance from the CI point. The relative error for the norm is less than 5% and the angle between two vectors is less than 1.5◦. We also plot the S1/S2 energy gap with respect to the distance in Fig. 3 for reference. As the distance increased from 0.001 Å to 0.2 Å to the CI point, the energies of S1 and S2 become 4 eV separated. Even with such a large energy gap, the derivative couplings given by these two different approaches are just about indistinguishable. B. S1/ S4 TDHF derivative coupling of LiH

As we have shown in the previous example, the RPA derivative coupling obtained from the original response theory and our direct differentiation of a pseudo-wavefunction ansatz

are practically equivalent near a CI point (and pretty far away too). We now study the opposite case: two well-separated electronic states. LiH is used as the test system and its TDHF S1/S4 derivative coupling at different bond lengths is computed with a 6-31G∗ basis for both approaches. In Fig. 4, we plot the absolute value of the derivative coupling with respect to the z direction of Li atom, using both our direct differentiation of pseudo-wavefunction ansatz and response theory. Relative energies of S0, S1, and S4 are plotted in Fig. 5 as reference. As shown in Fig. 4, when the bond length is increased from 1 Å to 3 Å, both S1/S4 TDHF derivative couplings of LiH usually behave similarly. Near 1.9 Å, however, the two couplings behave completely differently. While dPW 14 has a relatively consistent value in this range, dRT 14 changes dramatically in magnitude. This abnormal behavior of dRT 14

FIG. 5. (a) Relative energies for S0, S1, and S4 states and (b) energy gaps between these three states of LiH as a function of bond length. (Ω10 = Ω1 − Ω0, Ω41 = Ω4 − Ω1.)

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Mon, 02 Mar 2015 11:18:02

064114-11

Ou et al.

J. Chem. Phys. 142, 064114 (2015)

is not hard to explain if we analyze the S0, S1, and S4 relative energies of LiH. As shown in Fig. 5, on the left, the potential energy curves of these three states are smooth and there is no crossing when the bond length changes. However on the right, we plot the two energy gaps, Ω10 and Ω41. One can see that the two energy gaps cross when the bond length ≈1.9 Å, which is exactly the case where the derivative coupling from response theory is expected to be unphysical.

the Department of Energy under Award No. DE-SC0008694; J.E.S. acknowledges a Cottrell Research Scholar Fellowship and a David and Lucile Packard Fellowship.

IV. CONCLUSION

⟨ΨI |ΨJ[Q]⟩   X I R X J + RY I RY J 1  I J ˜ [Q]  * Rµλ σν [Q] µλ σν D µν h µν + =  +R X I RY J + RY I R X J + G˜ µνλσ ΩJ I µν νσ νσ µλ µλ µνλσ ,   ( ) [Q] ˜ IJ IJ − 1 Sµν Pµα D βν + DνI βJ Fα β + D µλ Pσν Π [Q] µνλσ   2 α β µν R X I R X J + RY I RY J *. νγX I δ βX J νγY I δ βY J +/ 1  [Q] ˜ ..+Rδ β Rνγ + Rδ β Rνγ // Sµν Pµα . X I Y J − Gα βγδ Y I X J/ 2 µνα βγδ ..+Rνγ Rβδ + Rνγ Rβδ // XI YJ Y I XJ ,+Rβδ Rνγ + Rβδ Rνγ -

In this paper, we have presented a detailed derivation for the RPA derivative couplings based on response theory. The final expression has been implemented in Q-Chem, and we have compared the resulting derivative couplings with those based on a pseudo-wavefunction ansatz (as obtained in our previous studies41,23). For small energy gaps between the two excited states, i.e., close to a CI, both sets of derivative couplings are identical up to an antisymmetric overlap term, see Appendix B 2. However, when the energy difference between two excited states equals the excitation energy of a third excited state, the derivative couplings from adiabatic TDDFT response theory exhibit a spurious pole. This unphysical behavior appears to be a consequence of an incorrect pole structure of quadratic and higher-order frequency dependent response properties within the adiabatic approximation,8 which neglects frequency dependence in the XC kernel and its derivatives. This incorrect analytical behavior also calls into question the validity of the adiabatic approximation for other non-linear response properties such as state-to-state transition moments and even non-linear polarizabilities. For TDHF response theory, which has the same analytical structure as adiabatic TDDFT, this problem has been discussed before,38 but its consequences for excited-state properties appear to have been largely ignored. Li and Liu have independently come to the same realization.31 An important conclusion of this work is thus that the adiabatic approximation of TDDFT is less robust for higherorder non-linear response properties than it has been generally assumed. The development of practical frequency-dependent XC kernels beyond the adiabatic approximation is also essential for states with double excitation character.42 A remedy of the spurious divergence is to evaluate the relaxed part of the transition density matrix at frequency zero instead of ΩJ I in Eq. (B10). This is exact at a CI, where the couplings are largest, and recovers the previously proposed pseudo-wavefunction approximation up to an antisymmetric overlap term. Thus, in the absence of better XC kernels, the pseudo-wavefunction approximation to the couplings may be the best available option for computing excited state to excited state couplings in a TDDFT framework.

ACKNOWLEDGMENTS

This work was supported (Q.O. and J.E.S.) by NSF CAREER Grant No. CHE-1150851 and (G.D.B. and F.F.)

APPENDIX A: DEFINITIONS FOR TERMS IN EQ. (100)

The final expression for response theory derivative couplings that we implemented and used throughout this paper is Eq. (100),

R X I R X J + RY I RY J *. γνX I βδX J γνY I βδY J +/ 1  [Q] ˜ ..+Rβδ Rγν + Rβδ Rγν // Sµν Pµα . X I Y J − Gα βγδ Y I X J/ 2 µνα βγδ ..+Rγν Rδ β + Rγν Rδ β // XI YJ Y I XJ ,+Rδ β Rγν + Rδ β Rγν R X I R X J + RY I RY J *. αγX I δ βX J αγY I δ βY J +/  .+Rδ β Rαγ + Rδ β Rαγ // Ξα βγδλσ [Q] ˜ − Sµν Pµλ Pνσ .. X I Y J Y I X J/ 2 ..+Rαγ Rβδ + Rαγ Rβδ // µνα βγδλσ XI YJ Y I XJ ,+Rβδ Rαγ + Rβδ Rαγ  ( ) 1 [Q] ˜ IJ + Sµν Pµα Pνδ DγI Jβ + D βγ 2 µνα βγδ  1  vo [Q] ov (γ + γ I J )bi Mbi × Gα βγδ − 2 bi I J   A[Q] − X iI a X iJ b −YiI aYiJ b CνaCµb Sµν µνiab



 (

) A[Q] X iI a X jJ a −YiI aYjJ a Cνi Cµ j Sµν

µνi j a

+

 ja

γ IvoJ − γ IovJ

  aj

A[Q] CµaCν j Sµν .

(A1)

µν

Terms in the above equation are defined as follows. 1. One- and two-electron integrals (a) Total one-electron integral h pq = h0pq + g pq ,

(A2)

where h0pq is the matrix element of the kinetic energy plus the external potential (Eq. (A3)) and g pq is the first derivative of the XC energy functional E xc = dr f xc (r)29 (Eq. (A4)), h0pq ≡ ⟨φ p |h0φq ⟩,

(A3)

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Mon, 02 Mar 2015 11:18:02

064114-12

Ou et al.

J. Chem. Phys. 142, 064114 (2015)



g pq ≡

drφ p (r)

pq

(b)

∂ f xc φq (r). ∂ ρ(r)

The total one-electron-integral derivative for TDDFT can be written as

(A4)

0[Q] [Q] h[Q] µν = h µν + g µν

Total two-electron integral G pq sr = Π pq sr + f pq sr ,

where Π pqsr is the Coulomb term plus whatever fraction of Hartree-Fock exchange is included in the DFT functional (cHF in Eq. (A6)), and f pq sr is the second derivative of the XC functional (Eq. (A6)) Π pqsr ≡ ⟨φ p φq |φ s φr ⟩ − cHF⟨φ p φq |φr φ s ⟩,

[Q] Y[Q] = h0[Q] ˜ µν + g µν µν + g

(A5)

Y[Q] ≡ h˜ [Q] µν + g µν .

(b)

Y[Q] [Q] [Q] , + f µνλσ = f˜µνλσ f µνλσ



Cµ pCν p = Pµν +



∂ 2 f xc φλ(r)φσ (r) ∂ ρ(r)2  [Q] ∂ 2 f xc + dr φ µ (r)φν (r) φλ(r)φσ (r) ∂ ρ(r)2  ∂ 2 f xc (φλ(r)φσ (r))[Q] + drφ µ (r)φν (r) ∂ ρ(r)2  ∂ 3 f xc + drφ µ (r)φν (r) φλ(r) ∂ ρ(r)3 γδ [Q] × φσ (r)Pγδ φγ (r)φδ (r) , (A18)  Y[Q] [Q] f µνλσ ≡ Pγδ Ξ µνλσγδ , (A19)

p

CµaCνa .

(A9)

a

Note that we may express the real-space density as ρ(r) = Pµν φ µ (r)φν (r). (b) The RPA excitation-amplitude matrices  XI Rµν = Cµa X iI aCνi , (A10) 

CµaYiI aCνi .

and Ξ µνλσγδ is the XC functional third derivative,  ∂ 3 f xc Ξ µνλσγδ ≡ drφ µ (r)φν (r) ∂ ρ(r)3 × φλ(r)φσ (r)φγ (r)φδ (r). (A20)

(A11)

ia

(c)

drφ µ (r)φν (r)

γδ

ia

RYµνI =

[Q]

[Q] f˜µνλσ ≡

(A8) 

(A17)

where

m

P˜µν =

Two-electron-integral derivatives The second derivative of the XC functional f [Q] is defined as

(A6)

′′ f pqsr ≡ ⟨φ p φq | f xc |φ s φr ⟩  ∂ 2 f xc φr (r)φ s (r). (A7) = drφ p (r)φq (r) ∂ ρ(r)2 pqsr

2. Density matrices (a) The general density matrices  Pµν = CµmCν m ,

(A16)

The generalized difference-density matrix  IJ D µν = Cµa (X iI a X iJ b +YiJ aYiI b )Cνb

The total two-electron-integral derivative for TDDFT can be written as

iab





Cµi (X iJ a X jI a +YiI aYjJ a )Cν j .

[Q] G[Q] = Π [Q] + f µνλσ µνλσ µνλσ

(A12)

[Q] Y[Q] = Π [Q] + f˜µνλσ + f µνλσ µνλσ  [Q] ≡ G˜ [Q] + Pγδ Ξ µνλσγδ . µνλσ

i ja

3. Derivative terms (a) One-electron-integral derivatives The first derivative of the XC functional g[Q] is defined as [Q] [Q] Y[Q] g µν ≡ g˜ µν + g µν ,

where



[Q] g˜ µν

∂ f xc φν (r) ≡ drφ µ (r) ∂ ρ(r)  [Q] ∂ f xc + dr φ µ (r)φν (r) ∂ ρ(r)  ∂ 2 f xc + drφ µ (r)φν (r) ∂ ρ(r)2 λσ [Q]

× Pλσ (φλ(r)φσ (r))[Q],  Y[Q] [Q] g µν ≡ Pλσ f µνλσ . Q

(A13)

(A21)

γδ

(c)

[Q] The mixed derivative Mbi [Q] Mbi



( µν

∂ 2E ∂ 2E [Q] h[Q] Sµν µν + ∂Θbi ∂h µν ∂Θbi ∂Sµν

∂ 2E G[Q] µνλσ ∂Θ ∂Π bi µνλσ µνλσ   =− CµbCνi +Cµi Cνb h[Q] µν +



)

(A22)

µν

(A14) (A15)

λσ

Here, denotes differentiation with respect to Becke weights.

 − Cµb PνσCλi +Cµi PνσCλb G[Q] µνλσ  + ε bCµbCνi + ε i Cµi Cνb µν

1  + Πα βγδ ( P˜β µ Pδν + P˜βν Pδ µ ) 2 α βγδ  [Q] × (CαbCγi +Cαi Cγb ) Sµν .

(A23)

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Mon, 02 Mar 2015 11:18:02

064114-13

Ou et al.

J. Chem. Phys. 142, 064114 (2015)

APPENDIX B: ORBITAL RESPONSE 1. Definition of Lagrangians

According to Ref. 33, L˜ vo(ωα , ω β ) and L˜ ov(ωα , ω β ) can be expressed as     X˜ id (ωα )Y˜j a (ω β ) + Y˜j a (ωα ) X˜ id (ω β ) G ab j d + X˜ id (ωα ) X˜ j a (ω β ) + X˜ j a (ωα ) X˜ id (ω β ) G adj b L˜ vo bi (ωα , ω β ) = adj





+



+



   X˜ l b (ωα )Y˜j a (ω β ) + Y˜j a (ωα ) X˜ l b (ω β ) G al j i + X˜ l b (ωα ) X˜ j a (ω β ) + X˜ j a (ωα )Xl b (ω β ) G ai jl

a jl

   X˜ j a (ωα )Y˜l d (ω β ) + Y˜j a (ωα ) X˜ l d (ω β ) Ξal j dbi + X˜ j a (ωα ) X˜ l d (ω β ) + Y˜j a (ωα )Y˜l d (ω β ) Ξadjl bi

adjl

   X˜ j a (ωα )Y˜j d (ω β ) + Y˜j d (ωα ) X˜ j a (ω β ) G abdi − X˜ j a (ωα )Y˜l a (ω β ) + Y˜j a (ωα ) X˜ l a (ω β ) G j bl i ,

adj

L˜ ov bi (ωα , ω β ) =



(B1)

a jl

   Y˜id (ωα ) X˜ j a (ω β ) + X˜ j a (ωα )Y˜id (ω β ) G ab j d + Y˜id (ωα )Y˜j a (ω β ) + Y˜j a (ωα )Y˜id (ω β ) G adj b

adj





+



+



   Y˜l b (ωα ) X˜ j a (ω β ) + X˜ j a (ωα )Y˜l b (ω β ) G al j i + Y˜l b (ωα )Y˜j a (ω β ) + Y˜j a (ωα )Yl b (ω β ) G ai jl

a jl

   Y˜j a (ωα ) X˜ l d (ω β ) + X˜ j a (ωα )Y˜l d (ω β ) Ξal j dbi + Y˜j a (ωα )Y˜l d (ω β ) + X˜ j a (ωα ) X˜ l d (ω β ) Ξadjl bi

adjl

   X˜ j a (ωα )Y˜j d (ω β ) + Y˜j d (ωα ) X˜ j a (ω β ) G dbai − X˜ j a (ωα )Y˜l a (ω β ) + Y˜j a (ωα ) X˜ l a (ω β ) Gl b j i .

adj

vo L ov I J and L I J are obtained by evaluating the residues of L (ωα , ω β ) and L˜ ov(ωα , ω β ) at ωα = Ω I and ω β = −ΩJ . In the vo AO basis, L ov I J and L I J can be expressed as

˜ vo

L vo bi =

(B2)

a jl

XI Jd YJ Id Rσν X i + Rσν Yi + G µνλσ CµbCλd * J d Y I Id XJ ,+X i Rνσ +Yi Rνσ µνλσd YJ Ib X I Jb  Yℓ Rσν X ℓ + Rσν + G µνλσ − CµℓCλi * J b Y I Ib XJ +X R +Y R νσ νσ ℓ ℓ , µνλσℓ XI XJ YI YJ  R Rσν + R Rσν + CγbCδi * µλX I Y J µλY I X J + Ξ µνλσγδ +R , νσ Rµλ + Rνσ Rµλ µνλσγδ  IJ + CνbCσi D µλ G µνλσ , (B3)



2. Behavior near a CI

In Ref. 23, we showed that derivative couplings constructed from our pseudo-wavefunction ansatz are consistent with the Chernyak-Mukamel equality and time-dependent response theory near an excited crossing in the limit of an infinite atomic orbital basis.41,23,12 In this subsection, we show that this consistency survives in a finite basis. The key observation is that the only meaningful difference between RT dPW I J and d I J comes from orbital response terms. a. The orbital response in the pseudo-wavefunction derivative coupling

According to Eqs. (47) and (48) in Ref. 23, the orbital response in dPW I J is

µνλσ

L ov bi

=



XJ Id Y I Jd Rσν X i + Rσν Yi

+ G µνλσ CµbCλd * I d Y J Jd XI ,+X i Rνσ +Yi Rνσ µνλσd XJ Ib Y I Jb  Yℓ Rσν X ℓ + Rσν + G µνλσ − CµℓCλi * I b Y J Jb X I ,+X ℓ Rνσ +Yℓ Rνσ µνλσℓ XI XJ YI YJ  Rσν R + Rσν Rµλ + + CγbCδi * X I µλY J Y I X J Ξ µνλσγδ +R R + R R νσ νσ µλ , µλ µνλσγδ  IJ + CνbCσi Dλµ Ω µνλσ . (B4) µνλσ

Finally, we can define a total Lagrangian L bi ≡

L ov + L vo . bi bi



1  1  vo =− . L bi Θ[Q] (L + L obiv )Θ[Q] bi bi ΩJ I bi ΩJ I bi bi

(B5)

This term can be simplified via the coupled-perturbed Hartree Fock (CPHF) equation Θ[Q] aj = −

( aj

=−

∂ 2E ∂Θbi ∂Θa j

) −1

[Q] Mbi

1 [Q] (A+ B)−1 j aib Mbi , 2 aj

(B6) (B7)

where Ma[Q] j refers to the mixed derivative terms (Eq. (A22)). The orbital response in dPW I J finally becomes

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Mon, 02 Mar 2015 11:18:02

064114-14



Ou et al.

J. Chem. Phys. 142, 064114 (2015)

1  1  [Q] (A+ B)−1 L bi Θ[Q] = j aib L bi Mbi . bi ΩJ I bi 2ΩJ I j aib

(B8)

According to Eq. (99), the orbital response in dRT I J is 

γ IvoJ − γ IovJ

aj



+ aj

 1 L a j ΘR[Q] aj . ΩJ I

(B9)

vo vo ov Here, L a j = L ov a j + L a j . γ I J and γ I J are given by

  −1 vo vo *γ I J + = − *A B+ + ΩJ I * I 0 + * L I J + , (B10) ov ov ,B A, γI J ,0 −I- , L I J and one has  1  (A + B)(γ IvoJ + γ IovJ ) + (L vo + L ov) . (B11) γ IvoJ − γ IovJ = − ΩJ I If we substitute Eq. (B11) into Eq. (B9), we find that all L terms are canceled. The orbital response in dRT I J can thus be rewritten as  1 (A+ B) j aib (γ IvoJ + γ IovJ )bi Θ[Q] (B12) aj . Ω J I j aib After applying the CPHF equation (Eq. (B7)), the final expression of the orbital response in dRT I J becomes  1 (A+ B) j aib (γ IvoJ + γ IovJ )bi Θ[Q] aj Ω J I j aib =−

1  vo [Q] (γ + γ IovJ )bi Mbi . 2ΩJ I bi I J

(B13)

c. Equivalence near a CI point

According to Eq. (B10), in the limit of ΩJ I → 0 (i.e., at a CI point), one has   −1 ov γ IvoJ + γ IovJ = −(A + B)−1 L vo (B14) I J + L I J = −(A + B) L. Referring to Eqs. (B13) and (B8), one finds that under such a condition, the orbital responses in dRT I J become exactly the same as the one in dPW . Therefore at a CI point, dRT I J is IJ A[Q] equivalent to dPW up to a factor of S (which can be ignored IJ when applying electron-translation factors35). This proves the equivalence we hypothesized. 1B.

Phys. 120, 7322 (2004). Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997 (1984). 8E. K. U. Gross and W. Kohn, Adv. Quantum Chem. 21, 255 (1990). 9M. Casida, in Recent Advances in Density Funtional Methods, edited by D. P. Chong (World Scientific, 1995), Vol. 1, pp. 155–192. 10M. Marques, N. Maitra, F. Nogueira, E. K. U. Gross, and A. Rubio, Fundamentals of Time-Dependent Density Functional Theory, Lecture Notes in Physics Vol. 837 (Springer, Heidelberg, 2012). 11N. T. Maitra, J. Chem. Phys. 125, 014110 (2006). 12V. Chernyak and S. Mukamel, J. Chem. Phys. 112, 3572 (2000). 13R. Baer, Chem. Phys. Lett. 364, 75 (2002). 14C. Hu, H. Hirai, and O. Sugino, J. Chem. Phys. 127, 064103 (2007). 15U. Werner, R. Mitri´ c, T. Suzuki, and V. Bonaˇci´c-Koutecký, Chem. Phys. 349, 319 (2008). 16C. Hu, H. Hirai, and O. Sugino, J. Chem. Phys. 128, 154111 (2008). 17C. Hu, O. Sugino, H. Hirai, and Y. Tateyama, Phys. Rev. A 82, 062508 (2010). 18R. Send and F. Furche, J. Chem. Phys. 132, 044107 (2010). 19E. Tapavicza, I. Tavernelli, and U. Rothlisberger, Phys. Rev. Lett. 98, 023001 (2007). 20I. Tavernelli, E. Tapavicza, and U. Rothlisberger, J. Mol. Struct.: THEOCHEM 914, 22 (2009). 21I. Tavernelli, B. F. E. Curchod, A. Laktionov, and U. Rothlisberger, J. Chem. Phys. 133, 194104 (2010). 22Z. Li and W. Liu, J. Chem. Phys. 141, 014110 (2014). 23Q. Ou, E. Alguire, and J. E. Subotnik, “Derivative couplings between timedependent density functional theory excited states in the random-phase approximation based on pseudo-wavefunctions: Behavior around conical intersections,” J. Phys. Chem. B (published online). 24X. Zhang and J. M. Herbert, J. Chem. Phys. 141, 064104 (2014). 25E. Tapavicza, G. D. Bellchambers, J. C. Vincent, and F. Furche, Phys. Chem. Chem. Phys. 15, 18336 (2013). 26M. E. Casida, K. C. Casida, and D. R. Salahub, Int. J. Quantum Chem. 70, 933 (1998). 27B. G. Levine, C. Ko, J. Quenneville, and T. J. Martínez, Mol. Phys. 104, 1039 (2006). 28B. G. Levine, J. D. Coe, and T. J. Martínez, J. Phys. Chem. B 112, 405 (2008). 29Q. Ou, S. Fatehi, E. Alguire, Y. Shao, and J. E. Subotnik, J. Chem. Phys. 141, 024114 (2014). 30It is worth noting that although Refs. 19–21 do not produce the correct final derivative coupling, Tavernelli et al. did invoke the correct formalism for deriving such a derivative coupling (i.e., by comparing many-body response theory with TD-DFT). Along the way, however, the authors incorrectly reduced the derivative coupling to a simple one-electron operator between singly excited auxiliary wavefunctions (via Casida’s assignment). 31Z. Li, B. Suo, and W. Liu, J. Chem. Phys. 141, 244105 (2014). 32X. Zhang and J. M. Herbert, J. Chem. Phys. 142, 064109 (2015). 33F. Furche, J. Chem. Phys. 114, 5982 (2001). 34This fact follows because γ ˜ (1)(t) is Hermitian; see Eqs. (43) and (42). 35S. Fatehi, E. Alguire, Y. Shao, and J. E. Subotnik, J. Chem. Phys. 135, 234105 (2011). 36In Ref. 33, the sign of ω α is not consistent with the definitions of X˜ (ω α ) and Y˜ (ω α ). In this work, we adopt the exact same formalism and nomenclature as in Ref. 33 but change the sign in front of ω α in Eq. (73) to be a plus sign to correct this earlier typo. 37N. C. Handy and H. F. Schaefer III, J. Chem. Phys. 81, 5031 (1984). 38E. Dalgaard, Phys. Rev. A 26, 42 (1982). 39J. Kong, M. S. Lee, A. M. Lee, S. R. Gwaltney, T. R. Adams, C. Ochsenfeld, A. T. B. Gilbert, G. S. Kedziora, V. A. Rassolov, D. R. Maurice et al., J. Comput. Chem. 21, 1532 (2000). 40Y. Shao, L. Fusti-Molnar, Y. Jung, J. Kussmann, C. Ochsenfeld, S. T. Brown, A. T. B. Gilbert, L. V. Slipchenko, S. V. Levchenko, D. P. O’Neill et al., Phys. Chem. Chem. Phys. 8, 3172 (2006). 41E. Alguire, Q. Ou, and J. E. Subotnik, “Calculating derivative couplings between time-dependent Hartree–Fock excited states with pseudowavefunctions,” J. Phys. Chem. B (published online). 42N. T. Maitra, F. Zhang, R. J. Cave, and K. Burke, J. Chem. Phys. 120, 5932 (2004). 7E.

b. The orbital response from derivative couplings according to response theory



6H. Lischka, M. Dallos, P. G. Szalay, D. R. Yarkony, and R. Shepard, J. Chem.

H. Lengsfield III, P. Saxe, and D. R. Yarkony, J. Chem. Phys. 81, 4549 (1984). 2P. Saxe, B. H. Lengsfield III, and D. R. Yarkony, Chem. Phys. Lett. 113, 159 (1985). 3B. H. Lengsfield III and D. R. Yarkony, J. Chem. Phys. 84, 348 (1986). 4B. H. Lengsfield III and D. R. Yarkony, Adv. Chem. Phys. 82(2), 1 (1992). 5Conical Intersections: Electron Structure Dynamics and Spectroscopy, edited by W. Domcke, D. R. Yarkony, and H. Köppel (World Scientific Publishing Co. Pte. Ltd., Singapore, 2004).

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Mon, 02 Mar 2015 11:18:02

First-order derivative couplings between excited states from adiabatic TDDFT response theory.

We present a complete derivation of derivative couplings between excited states in the framework of adiabatic time-dependent density functional respon...
1MB Sizes 0 Downloads 5 Views