冕 冕

Density functional theory of nonlinear triplet response properties with applications to phosphorescence

Ingvar Tunell, Zilvinas Rinkevicius, Olav Vahtras, and Paweł Sałek

Laboratory of Theoretical Chemistry, The Royal Institute of Technology, SE-10691, Sweden Trygve Helgaker

Department of Chemistry, University of Oslo, P.O. Box 1033 Blindern, N-0315 Oslo, Norway Hans A˚ gren^a)

Laboratory of Theoretical Chemistry, The Royal Institute of Technology, SE-10691, Sweden 共Received 15 July 2003; accepted 10 September 2003兲

We present density functional response theory generalized to triplet excitations. A method based on an exponential parametrization of the spin-dependent density operator is derived for the evaluation of linear and quadratic response functions for spin-dependent perturbations. The developed methodology is applicable to commonly available functionals, also hybrid functionals including exchange–correlation functionals at the general gradient-approximation level and fractional exact Hartree–Fock exchange. Illustrative calculations are presented for singlet–triplet transition moments and phosphorescence lifetimes, providing numerical data on these quantities for the first time using time-dependent density functional theory. © 2003 American Institute of Physics.

关DOI: 10.1063/1.1622926兴

I. INTRODUCTION

Although the overwhelming body of quantum methods still is based on the Born–Oppenheimer, nonrelativistic, Hamiltonian picture, with applications to singlet states, one has increasingly realized the important role of spin–orbit coupling for spectroscopy, reactivity, and catalysis even for light organic molecules. The incitement to transcend the tra- ditional approach has become particularly evident when con- sidering triplet or higher multiplicity states, since many prop- erties of these are intimately connected to spin–orbit coupling. A thesis of this work is that triplet state properties can best be studied theoretically through response theory.

This contention follows from a successful implementation some 10 years ago, which allowed the calculations of triplet state and spin–orbit related properties based on linear and quadratic response functions for singlet and triplet perturba- tions when no permutational symmetry in the two-electron operators is assumed and from which various triplet as well as singlet response properties could be derived. Several for- mal and computational advantages followed from this imple- mentation: The spin–orbit coupling matrix elements between singlet and triplet states could be evaluated as residues of 共multiconfiguration兲linear response functions, and therefore automatically become determined between orthogonal and noninteracting states, and in a way in which the sum-over- state values are obtained implicitly rather than by summa- tions of explicitly determined states. Furthermore, nonlinear triplet properties, like phosphorescence, were effectively ob- tained as spin-forbidden transition intensities determined from the spin–orbit coupling induced transitions between two electronic states of different multiplicity and obtained as

residues of quadratic response functions. The potential of the response theory for triplet state properties and its range of applications has been illustrated by a series of investigations covering different molecular phenomena, such as second- order energy contributions, intensity rearrangement in elec- tron spectra, predissociative lifetimes, triplet bands in ab- sorption spectra, intersystem crossings and reactivity, external heavy atom effects, spin–orbit coupling contribu- tions to NMR chemical shifts, phosphorescence spectra and radiative lifetimes of triplet states, to mention some of them.

A review of these applications has been given in Ref. 1.

The phenomenon of phosphorescence turned out to be a particularly rewarding implementation of response theory² for which the summation of perturbation theory contributions to the sum-over-state expression often is slowly converging.

A general spin-dependent response function implementation, given for self-consistent field and multiconfiguration SCF wave functions, was successfully applied to the phosphores- cence lifetime of small molecules like formaldehyde, mo- lecular nitrogen, oxygen, ozone, benzene, and naphthalene, to various azabenzenes and azanaphthalens and for studies of external heavy atom effects; see Ref. 1. Notwithstanding the success of these early applications of response theory to the phosphorescence phenomenon, there are two main obstacles in comparison with ‘‘ordinary’’ singlet excitation response theory or properties and which have prohibited further appli- cations to a larger class of molecules. The first is the involve- ment of the two-electron 共Breit–Pauli兲 spin–orbit operator, leading to the computation and, in particular, storing of a large number of two-electron integrals. The second limitation derives from the fact that the spin–orbit coupling involved describes the interaction of the electron spin with the elec- tron orbital angular momentum, which means that triplet

a兲Electronic mail: [email protected]

11024

0021-9606/2003/119(21)/11024/11/$20.00 © 2003 American Institute of Physics

perturbations must be considered. This brings about instabili- ties towards these triplet excitations for the commonly ap- plied random phase approximation 共response for Hartree–

Fock reference wave functions兲, which rather unpredictably appear for different systems. It seems though that even quite moderately correlated wave functions, e.g., MCSCF wave functions with small active spaces, overcome the triplet in- stabilities, but this remedy of the problem leads instead to a destruction of the good computational scaling of the HF/RPA approach and therefore to a limitation to small species.

The steep growth of the number of spin–orbit integrals owing to the spin–own–orbit and spin–other–orbit interac- tions is perhaps an even more severe obstacle towards appli- cations to larger systems. In the semiempirical context this problem has commonly been circumvented by simply ignor- ing the two-electron part of the spin–orbit operator in addi- tion to a neglect of the off-center terms.³ More modern ef- fective one-electron approaches include two-electron terms through so-called nuclear screening, through atomic mean- field approximations^4,5 共AMFI兲, or through the Ziegler ap- proximation of two-electron SO operators which includes spin–same-orbit as well as exchange–correlation potential SO operators, but neglects spin–other-orbit SO operators.⁶If such an effective one-electron approximation of the spin–

orbit coupling operator is crossed with an effective one- electron theory like density functional theory for response calculations, one can envisage a methodology that over- comes both the obstacles reviewed above and which can be predicted to be both accurate and applicable to phosphores- cence and other triplet state properties. With the present pa- per we intend to present a theory that accomplishes this goal, and illustrate it with a few initial calculations. The theory constitutes a nontrivial extension including triplet state op- erators of a nonlinear density functional theory for singlet state properties based on an explicit exponential parametri- zation of the density operator, recently advanced by the present authors.⁷

II. THEORY

The present work starts out from a density functional theory approach based on an explicit exponential parametri- zation of the density operator, and where the nonlinear re- sponse functions are derived using the Ehrenfest and the quasienergy variation principles, giving different but numeri- cally equivalent formulas.⁷The theory was implemented for dynamical nonlinear property calculations for all commonly available functionals, also hybrid functionals including exchange–correlation functionals at the general gradient- approximation level and fractional exact Hartree–Fock exchange.⁷The goal of the present work is to generalize this theory to triplet excitations and to triplet nonlinear response properties. We outline the derivations for linear and quadratic response functions and emphasize the generalizations that are required for spin-dependent perturbations.

The linear and quadratic response functions are defined by the expansion of an expectation value to second order

具^A典⫽具^A典0⫹

冕

^{d␻具具A;V典典␻e⫺i␻t}

⫹ ¹2

冕 冕

^{d␻1d␻2具具A;V,V典典␻1␻2e⫺i(␻1⫹␻2)t}

⫹¯, 共1兲

where A is a property and V is a perturbation. The starting point of our formalism is an exponential parametrization of the time evolution operator acting on a Kohn–Sham 共KS兲 determinant 兩0典

兩t典^{⫽e⫺␬ˆ (t)兩0}典^{, 共2兲}

where␬ˆ is an anti-Hermitian operator

␬^ˆ⫽_pq

兺

_␴ ^{␬pq␴a†p␴aq␴, 共3兲}

the matrix elements of which form the variational parameters of the theory and describe nonredundant rotations between occupied and virtual orbitals. The sum runs over orbitals p, q, and spin projections␴^{and a}p†␴

(a_q␴) is the standard cre- ation共annihilation兲operator of an electron in orbital␾p(␾q) with spin␴. The density of electrons with spin␴is evaluated as an expectation value of the␴-density operator

␳␴共r,t兲⫽具⁰兩e^␬ˆ␳^ˆ␴共r兲e^⫺␬ˆ兩0典^, 共4兲 which in turn is expressed in terms of the KS orbitals␾p as

␳^ˆ␴共r兲⫽

兺

_pq ^{␾p共r兲*␾q共r兲ap†␴aq␴. 共5兲}

We are considering a closed-shell spin-restricted theory with the same set of spatial orbitals for ␣^and␤spins. Neverthe- less, we treat ␳␣ and␳␤ as independent variables, because this is required in the treatment of spin-dependent perturba- tions. The equivalence between alpha and beta densities is imposed at the end.

A variation in the spin density energy functional

␦^E关␳␣,␳␤兴⫽

兺

_␴

冕

^d␶␦␳^␦␴^E共r兲␦␳_␴共r兲

⫽

兺

_␴

冕

^d␶␦␳^␦␴^E共r兲具⁰兩关␦␬^{ˆ ,}␳^ˆ_␴共r兲兴兩0典

⬅具⁰兩关␦␬^{ˆ ,Hˆ}兴兩0典^, 共6兲

defines the Hamiltonian operator

Hˆ⫽

兺

_␴

冕

^{d␶␳ˆ␴共r兲}␦␳^␦␴^E共r兲. 共7兲 The matrix elements of Hˆ form the spin-dependent KS ma- trix which determines the KS orbitals as well as the time evolution of the KS orbitals when the time-dependent pertur- bation is included.

In the presence of a time-dependent perturbation Vˆ (t), the Ehrenfest theorem applied to an arbitrary operator qˆ gives⁸

具⁰兩

冋

^qˆ,e␬ˆ

冉

^{Hˆ⫹Vˆ⫺idtd}

冊

^e⫺␬ˆ

册

^{兩0典⫽0, 共8兲}

and qˆ is chosen to be a column vector with nonredundant orbital rotation operators. The Hamiltonian Hˆ 共which de- pends on the density and hence implicitly on the perturbation V) and the parameters of ␬ˆ are expanded in orders of the perturbation strength

␬^ˆ⫽n

兺

⫽1 ␬^{ˆ(n) ; Hˆ}⫽n

兺

⫽0

Hˆ⁽ⁿ⁾. 共9兲

␬⁽ⁿ⁾ are determined by requiring that 共8兲 is to hold for all orders. We note that␣^and␤equivalence implies that

␦^E

␦␳␣⫽ ␦^E

␦␳␤, 共10兲

and that the zero-order Hamiltonian may be written as H⁽⁰⁾⫽

冕

^{d␶共␳ˆ␣共r兲⫹␳ˆ␤共r兲兲}␦␳^␦␣^E共r兲, 共11兲 i.e., as a singlet spin-conserving operator.

A. Linear response

Expanding共8兲to first order, we obtain in the frequency domain

具⁰兩关qˆ,关␬^ˆ␻,H(0)兴⫹H^␻⫹V^␻⫹␻␬^ˆ␻兴兩0典⫽0, 共12兲 where the Fourier transforms are defined by

␬^ˆ(1)⫽

冕

^{d␻␬ˆ␻e⫺i␻t, 共13兲}

Hˆ⁽¹⁾⫽

冕

^{d␻Hˆ␻e⫺i␻t, 共14兲}

Vˆ⫽

冕

^{d␻Vˆ␻e⫺i␻t. 共15兲}

The first-order Fourier transformed Hamiltonian of Eq. 共14兲 is

Hˆ␻⫽_␴␴

兺

⬘

冕 冕

^{d␶d␶⬘␳ˆ␴共r兲}␦␳␴共r^␦兲␦␳^2E␴⬘共r⬘兲␳_␴^␻_⬘共r⬘兲, 共16兲 where␳_␴^␻(r) is the Fourier transform of the first-order den- sity

␳_␴^␻共r兲⫽具⁰兩关␬^ˆ␻,␳^ˆ_␴共r兲兴兩0典^. 共17兲 It is instructive to look in detail at the role of the spin rank of the perturbation V. In practice, we consider a discrete sum of monochromatic frequencies

Vˆ⫽

兺

_k ^{Vˆ␻ke⫺i␻kt, 共18兲}

and each term is assumed to have a specific spin rank. We will repeatedly use the notation⁹

Vˆ␻k⫽

兺

_pq ^{vpq␻k共a†p␣aq␣⫹Skap␤† aq␤兲, 共19兲}

where S_kis a sign factor which implies a singlet operator for S_k⫽⫹1 and a triplet operator for S_k⫽⫺1. It is only neces- sary to consider the zero component of the triplet operator since the other components may be reduced to this form with the Wigner–Eckart theorem. A closer inspection of 共12兲 shows that meaningful solutions for a particular frequency can only be obtained if the excitation operators and ␬^have the same rank as V. This gives

␬^␻_pq^k_␤⫽S_k␬^␻_pq^k_␣^, 共20兲

␳_␤^␻k共r兲⫽S_k␳_␣^␻k共r兲, 共21兲 and assuming equivalence between ␣ ^and␤ derivatives the first-order Hamiltonian may be written as

Hˆ␻k⫽

冕 冕

^{d␶d␶⬘共␳ˆ␣共r兲⫹Sk␳ˆ␤共r兲兲}

冉

^␦␳␣共r^␦兲²␦␳^E␣共r⬘兲

⫹S_k ␦^2E

␦␳␣共r兲␦␳␤共r⬘兲

冊

^{␳␣␻k共r⬘兲. 共22兲}

We observe that the first-order Hamiltonian has the same spin rank as the perturbation. Now, we have all definitions in place for solving the first-order response equation共12兲for a particular frequency and perturbation; the linear response function in Eq.共1兲is evaluated as

具具^A;V典典␻k⫽具^0兩关␬^ˆ␻k,Aˆ兴兩0典^{. 共23兲}

B. Quadratic response

To evaluate the quadratic response function with respect to two perturbations V^␻1and V^␻2, we need in addition to the first-order parameters corresponding to these perturbations/

frequencies the second-order parameters obtained from the expansion of 共8兲to second order

具⁰兩关qˆ,Hˆ␻1,␻2⫹2 P₁₂关␬^{ˆ␻1,Hˆ␻2}⫹Vˆ␻2兴⫹共␻1⫹␻2兲␬^ˆ␻1,␻2

⫹P₁₂关␬^ˆ␻1,关␬^{ˆ␻2,Hˆ(0)}兴兴⫹共␻2⫺␻1兲关␬^ˆ␻1,␬^ˆ␻2兴兴兩0典^⫽0, 共24兲 where P₁₂is the symmetrizing operator

P₁₂f共1兲g共2兲⫽f共1兲g共2兲⫹f共2兲g共1兲

2 , 共25兲

and the second-order terms have the transforms

␬^ˆ(2)⫽1

2

冕 冕

^{d␻1d␻2␬ˆ␻1,␻2e⫺i(␻1⫹␻2)t, 共26兲}

Hˆ⁽²⁾⫽1

2

冕 冕

^{d␻1d␻2Hˆ␻1,␻2e⫺i(␻1⫹␻2)t. 共27兲}

Note that the second-order Hamiltonian has second-order contributions from first- and second-order densities

Hˆ␻1,␻2⫽

兺

␴␴⬘

冕 冕

^{d␶d␶⬘␳ˆ␴共r兲}␦␳␴共r^␦兲␦␳^2E␴⬘共r⬘兲

⫻␳_␴^␻_⬘^1,␻2共r⬘兲

⫹P_12␴␴

兺

_⬘␴⬙

冕 冕 冕

^{d␶d␶⬘d␶⬙␳ˆ␴共r兲}

⫻ ␦^3E

␦␳_␴共r兲␦␳_␴⬘共r⬘兲␦␳_␴⬙共r⬙兲␳_␴^␻1共r⬘兲␳_␴^␻_⬘²共r⬙兲, 共28兲 and the second-order densities have contributions from first- and second-order parameters

␳_␴^␻1,␻2共r兲⫽具⁰兩关␬^ˆ␻1,␻2,␳^ˆ␴共r兲兴

⫹P₁₂关␬^ˆ␻1,关␬^ˆ␻2,␳^ˆ␴共r兲兴兴兩0典^. 共29兲 Just as for linear response, we note that Eq.共24兲implies that the spin rank of the second-order operators must be the same as the excitation operators and the combination of the two first-order operators

␬^␻_pq¹_␤^,␻2⫽S₁S₂␬_pq^␻1_␣^,␻2, 共30兲

␳_␤^␻1,␻2共r兲⫽S₁S₂␳_␣^␻1,␻2共r兲, 共31兲 and the ␣ ^and ␤ equivalence implies that the second-order Hamiltonian may be written as

Hˆ^␻1,␻2⫽

冕 冕

^{d␶d␶⬘共␳ˆ␣共r兲⫹S1S2␳ˆ␤共r兲兲}

冉

^␦␳␣共r^␦兲²␦␳^E␣共r⬘兲

⫹S₁S₂ ␦^2E

␦␳␣共r兲␦␳␤共r⬘兲

冊

^{␳␣␻1␻2共r⬘兲}

⫹

冕 冕 冕

^{d␶d␶⬘d␶⬙共␳ˆ␣共r兲⫹S1S2␳ˆ␤共r兲兲}

⫻

冉

^␦␳␣共r兲␦␳^␦␣³共^Er⬘兲␦␳␣共r⬙兲⫹共S₁⫹S₂⫹S₁S₂兲

⫻ ␦^3E

␦␳␣共r兲␦␳␣共r⬘兲␦␳␤共r⬙兲

冊

^{␳␣␻1共r⬘兲␳␣␻2共r⬙兲.}

The ingredients have now been defined for evaluating the quadratic response functions from the solutions of Eqs. 共12兲 and共24兲, which become

具具^A;V,V典典␻₁,␻₂⫽具⁰兩关␬^␻1,␻2,A兴

⫹P₁₂关␬^␻1,关␬^␻2,Aˆ兴兴兩0典^. 共32兲 The high complexity of the first- and second-order Hamilto- nians of Eqs.共16兲and共28兲is hidden in the general formal- ism of functional calculus. For general GGA-type function- als which are expressed as integrals over energy densities that are functions of densities and density gradients, we de- vote the Appendix to the Hamiltonian operators to second order.

III. APPLICATION

In this section we recapitulate the origin of the phospho- rescence effect. The transition between the two states of dif- ferent multiplicities 共singlet and triplet兲 is realized by the spin–orbit coupling perturbing the two states producing non- pure initial and final spin states. The expansion describing these spin contaminations can then be written as

␮⫽具^S0兩r兩T₁典⫽

兺

_n ^{具S0兩r兩S}Eⁿ_T^{典具Sn兩HˆSO兩T1典}

1⫺E_S

n

⫹

兺

n

具^S0兩Hˆ

SO兩T_n典具^Tn兩r兩T₁典 E_S

0⫺E_T

n

, 共33兲 which can then be inserted in the regular expression for phosphorescence lifetimes

1

␶^{⫽const. 共⌬E兲3}␮^2. 共34兲

The use of共33兲is usually referred to as the ‘‘sum-over- states’’ method. The convergence for such intermediate state summations is often slow already for linear properties, and even more so for nonlinear properties such as phosphores- cence moments, where individual contributions of intermedi- ate states are of arbitrary sign. This has been demonstrated by, e.g., Langhoff and Davidson¹⁰in CI calculations of phos- phorescence moments. Thus, ‘‘few-state-models’’ based on a truncated summation are risky in many cases. On the other hand, when such models are applicable they may provide an interpretation of the property value and eventually a

‘‘structure-property’’ relation. Moreover, summations be- come very cumbersome for well-parametrized methods, since the summation then may simply include too many terms to be feasible. In response theory the sum-over-state expression is implicitly obtained by solving equations sys- tems without explicitly addressing the excited states. In this way the sum-over-state value is readily obtained in cases where the explicit summation is impossible.

From the quadratic response function 具具^{A;V␻1,V␻2}典典␻₁⫹i⑀,␻₂⫹i⑀, we can extract many different molecular properties by evaluating its single and double resi- dues. The single residue can be written as

lim

␻2→␻f⫺i⑀共␻2⫺␻f⫹i⑀兲具具^{A;V␻1,V␻2}典典⫺␻₁⫹i⑀,␻₂⫹i⑀

⫽

兺

_k

冋

^{具0兩A兩k典关具⫺k兩␻V⫺␻1⫹1␻兩ff典⫺⫺␻␦kk f⫹具0i⑀兩V⫺␻1兩0典兴}

⫺具⁰兩V^⫺␻1兩k典关具^k兩A兩f典⫺␦k f具⁰兩A兩0典兴

⫺␻1⫹␻k⫹i⑀ 具^f兩V^␻f兩0典^, 共35兲 where we have used the notation of Ref. 2.

In this case we are interested in the phosphorescence matrix elements ␮. By letting ␻1⫽0, A⫽r, and V^⫺␻1

⫽Hˆ_SO we can identify the right-hand side of Eq. 共35兲with Eq. 共33兲. That is, the residue of the response function 具具^r;HˆSO,V^␻2典典0,␻₂ gives the phosphorescence matrix ele-

ments. Just as we, from a response function, get the phos- phorescence matrix elements we can also get individual components of the moment for each state in Eq. 共33兲 by evaluating the corresponding residues. Naturally, this is com- putationally very demanding, requiring solving linear re- sponse equations for the具^S0兩r兩S_n典 ^and具^S0兩Hˆ

SO兩T₁典 ^compo- nents and quadratic response equations for the具^{Sn兩HˆSO兩T1}典 and 具^Tn兩r兩T1典 components. However, by computing these components we can examine the convergence of the ‘‘sum- over-state’’ expression, Eq.共33兲.

IV. CALCULATIONS AND RESULTS

Excitation energies and phosphorescence lifetimes have been calculated for a selected set of small molecules. The choice of molecules was motivated by the ability to perform high-accuracy ab initio calculations such as full CI共FCI兲or large complete/restricted active space self-consistent field 共CASSCF/RASSCF兲^11,12calculations. While DFT in itself is a method more suitable for larger systems than the men- tioned high-accuracy methods, we have instead of the full Breit–Pauli spin-operator used AMFI, which constitutes an approximation to the two-electron spin–orbit contribution.

This improves the scaling of the TD-DFT method further.

Following the original work by Heß et al.,^4,5 this approach has been proven successful for a number of properties,^4,5,13,14 including phosphorescence.¹⁵ Furthermore, no two-electron operators are well defined within the DFT implementations to our knowledge, although the foundation for such has been derived.¹⁶

The geometry optimizations were performed with the

GAUSSIAN 98¹⁷ software package, and all other properties with theDALTON¹⁸ quantum chemistry program.

For lithium hydride共LiH兲and formaldehyde (H₂CO) we used the experimental geometries from Herzberg¹⁹ and Raynes,²⁰respectively. The basis sets used for LiH excitation energies and transition moments were the correlation consis- tent valence double, triple, and quadruple zeta basis sets^21,22 共cc-pVXZ, X⫽D, T, and Q兲and the Pople split valence^23–27 and split valence triple zeta basis sets^25–28with diffuse func- tions 关6-31⫹⫹G** and 6-311⫹⫹G(3d f ,3pd)] in both the B3LYP and FCI calculations.

For beryllium dihydride (BeH₂) the structure was opti- mized using the unrestricted B3LYP^29,30 共UB3LYP兲 func- tional with the cc-pVQZ basis set. The TD-DFT calculations were then performed with the same basis sets as for lithium hydride. As reference we performed large RAS calculations since FCI would be computationally too demanding for all but the smallest basis sets. The RAS space was set up with an RAS2 space of 3a₁,1b₁,1b₂, and 1a₂ orbitals and an RAS3 space including all remaining virtual orbitals of the basis set. We allowed up to two electrons in this RAS3, that is, single and double excitations are accounted for. This col- lects a major part of the dynamical correlation energy which was verified by an FCI calculation with the smallest basis set used.

For formaldehyde we used the atomic natural orbital ba- sis sets共ANO兲of Widmark et al.³¹ with the following con- tractions: C,O: (14s9 p4d3 f )/关4s3 p2d1 f兴, H: (8s4 p3d)/

关3s2p1d兴 labeled ANO-A and C,O: (14s9 p4d3 f )/

关3s2 p1d兴, H: (8s4 p3d)/关2s1 p兴 labeled ANO-B, aug- mented and nonaugmented cc-pVTZ and cc-pVQZ,^21,22 where the TD-DFT calculations were performed using Dirac–Vosko–Wilk–Nusair³² 共LDA兲, Becke–Lee–Yang–

Parr,^30,33共BLYP兲and Becke3–Lee–Yang–Parr^29,30共B3LYP兲 functionals. Similar to in Ref. 2, we perform calculations of the singlet and triplet intermediate state components in Eq.

共33兲. In Ref. 2 the calculations were carried out at the SCF level on the singlet state geometry. In this work we account for electron correlation, performing both DFT and CASSCF calculations on the triplet state geometry more in accordance with the phosphorescence phenomenon.

A. Lithium hydride

An accurate test of a new phosphorescence method is aggravated by the fact that only small, light, molecules are accessible by the benchmarking FCI method, while these molecules at the same time possess exceedingly small singlet–triplet transition moments, and transition moments corresponding to lifetimes over a second are extremely sen- sitive to the parametrization and difficult to converge. De- spite this unfortunate circumstance we have chosen LiH and BeH₂ as the first two test cases for the new DFT phospho- rescence method.

In LiH we study the two lowest triplet states, ³⌺^⫹ and

3⌸, although the former state is strictly repulsive, thus with- out measurable phosphorescence. However, we can still gather information by considering the transition moments.

All properties have been calculated at a lithium–hydrogen distance of 1.5953 Å.¹⁹The results of the DFT-LR excitation energy calculations are listed in Table I and the transition moments from the DFT-QR calculations are given in Table II and III. Using the FCI calculations as reference, DFT-LR gives excitation energies as well as other methods, such as CASSCF with moderate active spaces but at a fraction of the cost. Comparing with FCI the energies are underestimated by 10%–20%. By examining the transition moments we can, for the ³⌺^⫹ state, see a very good agreement between the TD- DFT and FCI results already for moderately sized basis sets.

As noted above, due to the very small transition moments it is difficult to reach convergence prior to the basis set limit of the FCI method. In view of that fact, the comparison be- tween the B3LYP and FCI values for the ³⌸ state is also reassuring for the DFT method, although rather large basis sets are needed.

TABLE I. Excitation energies in a.u. of LiH³⌺^⫹and³⌸.

Basis set

Excitation energies

3⌺^⫹/B3LYP ³⌺^⫹/FCI ³⌸/B3LYP ³⌸/FCI

cc-pVDZ 0.0981 0.1138 0.1416 0.1517

cc-pVTZ 0.0974 0.1187 0.1417 0.1552

cc-pVQZ 0.0979 0.1195 0.1418 0.1557

6-31⫹⫹G** 0.0995 0.1120 0.1467 0.1533

6-311⫹⫹G(3d f ,3pd) 0.0971 0.1189 0.1418 0.1551

B. Beryllium dihydride

Beryllium dihydride has a C_2v structure in the triplet states. The geometry, a beryllium–hydrogen distance of 1.4318 Å, and an HBeH angle of 42.5°, as well as the ener- gies of the triplet states are similar to values of previous works.³⁴The three lowest triplet states,³B₂, ³A₂, and ³A₁, have been examined. The energies and the phosphorescence lifetimes are presented in Table IV. The DFT-LR excited state energies are about 10%–20% lower than for the RAS calculations, and the phosphorescence lifetimes are consis- tently longer. As for LiH the transition moments are very small, and as the lifetimes are well over a second, the agree- ment between DFT and RASSCF must be considered as very satisfactory.

C. Formaldehyde

In formaldehyde we compare our results to the CASSCF phosphorescence calculations of Vahtras et al.²in 1992. The same experimental geometry is used: Carbon–hydrogen dis- tance 1.10 Å, carbon–oxygen 1.28 Å, HCH angle 116.5°, and hydrogen out-of-plane angles 38.5°. The excitation en- ergies and phosphorescence lifetimes of the ³A⬙ ^{state are} given in Table V. The B3LYP lifetimes are roughly three times as long as the multiconfigurational results. We can ex- pect the real lifetimes to lie in between the DFT and CASSCF results. The excitation energy at the triplet state geometry is unknown, while calculations and experiments on

the singlet state geometry imply that the excitation energies are underestimated by about 10% in the DFT-LR calcula- tions, while CASSCF overestimates it by around 10%. For the lifetimes the basis set effect is minor, while the different functionals give some variations, albeit much smaller than the variation between DFT and CASSCF. As the CASSCF method does not include dynamical correlation effects it can- not serve the same benchmarking purpose as FCI or exten- sive RASSCF did for the LiH and BeH₂ calculations dis- cussed above. We interpret the difference of a factor 3 to 4 between the DFT and CASSCF results for H₂CO as an effect of inclusion of dynamical correlation in the DFT method.

In Tables VI, VII, VIII, and IX we present the conver- gence of the largest component of the phosphorescence ma- trix element (r⫽z and H_SO^z ). We see no change in the con- vergence for correlated wave functions. However, this could be expected since singlet and triplet states contributing to␮ have a single excitation character and therefore relatively small contributions to␮from doubly excited states. For the B3LYP/ANO-B calculation the phosphorescence matrix ele- ment is 1.945 and with 11 excited states we get共by looking at Tables VI and VII兲 ⫺0.8988⫺2.1763⫽⫺3.0751, where the sign can be neglected as only representing an arbitrary phase factor. In the same way for CASSCF/ANO-B we have a moment of 2.149, while the 11 states accounted for sums up to 1.0430⫺(⫺1.4358)⫽2.4788共Tables VIII and IX兲.

TABLE III. FCI transition dipole moment components for LiH.共Units of 10^⫺3a.u.)

Basis set Polarization ³⌺^⫹(x) ³⌺^⫹(y ) ³⌸(x) ³⌸(z)

cc-pVDZ x 0.003 ⫺0.015

y ⫺0.003

z ⫺0.008

cc-pVTZ x 0.003 0.018

y ⫺0.003

z 0.011

cc-pVQZ x ⫺0.003 0.021

y 0.003

z 0.011

6-31⫹⫹G** x 0.002 0.011

y ⫺0.002

z 0.006

6-311⫹⫹G** x ⫺0.003 ⫺0.019

y 0.003

z ⫺0.011

TABLE IV. Energies and lifetimes of³B₂,³A₂, and³A₁states in BeH₂.共Energies in a.u. and lifetimes in s.兲

Basis set

3B₂/B3LYP ³B₂/RAS^{a 3}A₂/B3LYP ³A₂/RAS^{a 3}A₁/B3LYP ³A₁/RAS^a

E ␶ ^E ␶ ^E ␶ ^E ␶ ^E ␶ ^E ␶

cc-pVDZ 0.0983 165.56 0.1190 57.104 0.1579 45.579 0.1741 25.273 0.1732 13.172 0.2002 8.0690

cc-pVTZ 0.1006 86.271 0.1250 35.275 0.1589 28.317 0.1783 20.769 0.1735 8.5271 0.2051 6.1719

cc-pVQZ 0.1005 62.798 ¯ ¯ 0.1586 12.546 ¯ ¯ 0.1733 7.1526 ¯ ¯

6-31⫹⫹G** 0.1006 147.36 0.1256 49.966 0.1609 29.971 0.1866 21.842 0.1741 11.421 0.2060 7.6714 6-311⫹⫹G(3d f ,3pd) 0.1004 115.91 0.1330 37.643 0.1583 25.449 0.1945 21.115 0.1732 8.5299 0.2133 6.0298

aSee the text for a description of the RAS space.

TABLE II. DFT共B3LYP兲transition dipole moment components for LiH.

共Units of 10^⫺3a.u.)

Basis set Polarization ³⌺^⫹(x) ³⌺^⫹(y ) ³⌸(x) ³⌸(z)

cc-pVDZ x ⫺0.002 ⫺0.010

y 0.002

z ⫺0.009

cc-pVTZ x ⫺0.003 0.015

y 0.003

z 0.013

cc-pVQZ x 0.003 ⫺0.019

y ⫺0.003

z ⫺0.013

6-31⫹⫹G** x 0.001 ⫺0.007

y ⫺0.001

z ⫺0.007

6-311⫹⫹G** x ⫺0.003 0.017

y 0.003

z 0.013

TABLE V. Excitation energies and phosphorescence lifetimes of the³A⬙^{state of H}2CO.共Energies in a.u. and lifetimes in s.兲

Basis set

LDA BLYP B3LYP CASSCF

E ␶ ^E ␶ ^E ␶ ^E ␶

ANO-A 0.0757 0.0981 0.0750 0.1142 0.0771 0.0735 0.1065^a 0.0238^a ANO-B 0.0765 0.0966 0.0757 0.1107 0.0777 0.0731 0.1078^a 0.0227^a

cc-pVTZ 0.0760 0.1104 0.0753 0.1316 0.0771 0.0850 0.1072 0.0257

aug-cc-pVTZ 0.0766 0.0918 0.0759 0.1047 0.0780 0.0707 0.1079 0.0222

cc-pVQZ 0.0764 0.0994 0.0757 0.1168 0.0776 0.0765 0.1079 0.0235

aug-cc-pVQZ 0.0765 0.0856 0.0758 0.0973 0.0778 0.0660 0.1080 0.0211

aFrom Ref. 2. Active space: 12 electrons in 8a⬘^4a⬙共full valence space兲.

TABLE VI. B3LYP/ANO-B singlet states contribution关first term of Eq.共33兲兴. State 具^S0兩z兩Sn典 具^Sn兩HˆSO

z 兩T1典^{a ES}n Term 1^a Sum^a

1 0.2167 0.1180 0.2708 ⫺0.1324 ⫺0.1324

2 ⫺0.5222 0.0162 0.2721 0.0435 ⫺0.0889

3 0.7899 0.2235 0.3072 ⫺0.7695 ⫺0.8584

4 ⫺0.1534 ⫺0.0247 0.3484 ⫺0.0140 ⫺0.8724

5 0.2259 0.0340 0.3855 ⫺0.0249 ⫺0.8973

6 0.3517 0.0176 0.3930 ⫺0.0196 ⫺0.9170

7 0.1710 0.0143 0.4105 ⫺0.0074 ⫺0.9243

8 0.1529 ⫺0.0016 0.4301 0.0007 ⫺0.9236

9 0.4781 0.0219 0.4508 ⫺0.0281 ⫺0.9517

10 ⫺0.3227 ⫺0.0044 0.4524 ⫺0.0038 ⫺0.9555

11 0.9122 ⫺0.0240 0.4635 0.0567 ⫺0.8988

aUnits of 10^⫺3a.u.

TABLE VII. B3LYP/ANO-B triplet states contribution关second term of Eq.共33兲兴.

State 具^S0兩Hˆ_SO^z 兩Tn典^a 具^Tn兩z兩T1典 ^ET_n Term 2^a Sum^a

1 0.2772 0.5679 0.0777 2.0251 2.0251

2 0.0206 ⫺0.0105 0.2349 ⫺0.0009 2.0242

3 0.0054 0.0916 0.2747 0.0018 2.0260

4 0.0718 0.4640 0.2980 0.1118 2.1379

5 ⫺0.0025 0.4027 0.3164 ⫺0.0032 2.1347 6 ⫺0.0260 ⫺0.8942 0.3478 0.0669 2.2015 7 ⫺0.0076 0.2012 0.3790 ⫺0.0040 2.1975

8 0.0001 ⫺0.0394 0.4075 ⫺0.0000 2.1975

9 0.0200 0.0561 0.4356 0.0026 2.2001

10 0.0546 ⫺0.1022 0.4428 ⫺0.0126 2.1875

11 ⫺0.0278 0.1796 0.4475 ⫺0.0112 2.1763

aUnits of 10^⫺3a.u.

TABLE VIII. CASSCF/ANO-B singlet states contribution关first term of Eq.共33兲兴. State 具^S0兩z兩Sn典 具^Sn兩HˆSO

z 兩T1典^{a ES}n Term 1^a Sum^a

1 ⫺0.0000 0.0000 0.2665 0.0000 0.0000

2 ⫺0.4777 0.1242 0.2918 0.3224 0.3224

3 0.9306 ⫺0.1528 0.3196 0.6714 0.9938

4 ⫺0.1585 0.0899 0.3247 0.0657 1.0595

5 ⫺0.0000 0.0000 0.3884 0.0000 1.0595

6 ⫺0.0730 0.0172 0.4086 0.0042 1.0637

7 ⫺0.2722 ⫺0.0290 0.4277 ⫺0.0247 1.0390

8 ⫺0.3481 ⫺0.0392 0.4508 ⫺0.0398 0.9992

9 ⫺0.1638 0.0847 0.4659 0.0387 1.0379 10 ⫺0.1390 ⫺0.0175 0.4865 ⫺0.0064 1.0315

11 0.1510 ⫺0.0294 0.4931 0.0115 1.0430

aUnits of 10^⫺3a.u.

D. Discussion

The present work represents the first derivation and implementation of nonlinear triplet excitation properties within the framework of density functional theory. The ap- proach is based on an explicit exponential parametrization of the density operator, with the共nonlinear兲response functions derived using the Ehrenfest and the quasienergy variation principles, giving different but numerically equivalent for- mulas. The theory was implemented for dynamical nonlinear property calculations for all common functionals as well as hybrid functionals including exchange–correlation function- als at the general gradient-approximation level and fractional exact Hartree–Fock exchange. The goal with the present work was to generalize the corresponding theory for singlet excitations⁷ to triplet excitations and to triplet nonlinear re- sponse properties, thereby reaching a set of important mo- lecular properties that are associated with the triplet state by time-dependent density functional theory.

Although the ultimate goal of the theory is to make trip- let state property calculations accessible for large molecules, thereby solving the triplet instability problem for the Hartree–Fock-based random phase approximation as well as improving upon the precision of that commonly used method in general, we confine ourselves for illustration purposes in this first work to a set of small molecules focusing on the phosphorescence effect. DFT phosphorescence transition moments and triplet state radiative lifetimes have thus been evaluated for the LiH, BeH₂, and H₂CO molecules.

The small molecules available for the benchmarking FCI and RAS calculations are unfortunately at the same time too light to produce appreciable triplet transition moments and short phosphorescence lifetimes. The comparison between FCI and the DFT method is therefore aggravated by the very sensitive parametric dependence. The comparison of phos- phorescence values with FCI, and also with RASSCF, is nev- ertheless favorable for DFT; considering the very long life- times satisfactory agreement is reached already for double or triple zeta basis sets for the LiH and BeH₂ molecules.

Evaluating the calculated data in some detail, one finds that there is a general trend of DFT to somewhat overesti- mate the phosphorescence lifetimes. However, considering Eq. 共34兲 the phosphorescence lifetime is a very sensitive quantity; an underestimation in energy of approximately

10% gives an increase in lifetime of about a factor of 1.37.

With decent basis sets the transition moments show good agreement, although at this stage too few systems have been studied to establish this as being the general case. The AMFI approximation has been found to be in good agreement with calculations¹⁵using the full Breit–Pauli spin–orbit operator.

This was also verified for the CASSCF calculations on form- aldehyde presented here.

For formaldehyde we explored the effects of different functionals as we believe that most of the remaining discrep- ancies of DFT can be attributed to these; here, we explored the B3LYP, BLYP, and LDA functionals. It seems desirable to include the effects of asymptotically corrected 共AC兲 po- tentials in future works. In the outer region the electron should only feel the 1/r field of the atom it leaves behind while most employed potentials decay much faster, some- thing that the AC functional should remedy according to dis- cussions by van Leeuwen et al.³⁵ As the response calcula- tions involve numerous higher lying states this correction may give a significant improvement to the excitation ener- gies as have been shown by, e.g., Cai et al.³⁶ for singlet excitations.

For small molecules the inclusion of diffuse basis func- tions seems to be more important than the precise description of the inner regions. In formaldehyde the aug-cc-pVTZ cal- culations give better results than with the cc-pVQZ basis set.

At a first glance the need for a large basis set seems to be more important for the DFT method than for other correlated methods such as, e.g., CASSCF. In particular, the longer components of the lifetimes are much more sensitive to basis set effects than are their slow counterparts. In the millisecond region of formaldehyde we do not notice any major basis set effect, while in the hydrides, where the lifetimes are on the order of 1 – 10² s the basis set effect is substantial. It is there- fore our recommendation to use large basis sets if the life- times can be expected to be larger than seconds. From a few preliminary calculations on larger systems involving 20–25 atoms 共to be reported separately兲, it seems that the conclu- sions also carry over to these systems. This should not prove to be an obstacle as implementations of linear scaling tech- niques in DFT by the present authors now support calcula- tions with several thousand basis functions. As the first ap- plication of phosphorescence using the DFT-QR formalism,

TABLE IX. CASSCF/ANO-B triplet states contribution关second term of Eq.共33兲兴.

State 具^S0兩Hˆ_SO^z 兩Tn典^a 具^Tn兩z兩T1典 ^ETn Term 2^a Sum^a 1 ⫺0.2887 0.5868 0.1078 ⫺1.5718 ⫺1.5718

2 0.0000 ⫺0.0000 0.2534 ⫺0.0000 ⫺1.5718

3 ⫺0.0186 ⫺0.0295 0.2786 0.0020 ⫺1.5698

4 ⫺0.0073 0.0668 0.3166 ⫺0.0015 ⫺1.5714

5 0.0770 0.5572 0.3442 0.1247 ⫺1.4467

6 ⫺0.0256 ⫺0.4314 0.3539 0.0312 ⫺1.4155

7 ⫺0.0175 0.0995 0.3880 ⫺0.0045 ⫺1.4200

8 ⫺0.0151 0.2409 0.4441 ⫺0.0082 ⫺1.4282

9 ⫺0.0041 0.0328 0.4513 ⫺0.0003 ⫺1.4285

10 0.0499 ⫺0.2524 0.4821 ⫺0.0261 ⫺1.4546

11 ⫺0.0446 ⫺0.2088 0.4951 0.0188 ⫺1.4358

aUnits of 10^⫺3a.u.