Density-functional theory of linear and nonlinear time-dependent molecular properties

Density-functional theory of linear and nonlinear time-dependent molecular properties

Paweł Sałek

Laboratory of Theoretical Chemistry, The Royal Institute of Technology, SE-10044 Stockholm, Sweden and Department of Chemistry, University of Oslo, P.O. Box 1033 Blindern, N-0315 Oslo, Norway

Olav Vahtras

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

Department of Chemistry, University of Oslo, P.O. Box 1033 Blindern, N-0315 Oslo, Norway Hans Ågren

Laboratory of Theoretical Chemistry, The Royal Institute of Technology, SE-10044 Stockholm, Sweden 共Received 10 June 2002; accepted 4 September 2002兲

We present density-functional theory for linear and nonlinear response functions using an explicit exponential parametrization of the density operator. The response functions are derived using two alternative variation principles, namely, the Ehrenfest principle and the quasienergy principle, giving different but numerically equivalent formulas. We present, for the first time, calculations of dynamical hyperpolarizabilities for hybrid functionals including exchange-correlation functionals at the general gradient-approximation level and fractional exact Hartree–Fock exchange. Sample calculations are presented of the first hyperpolarizability of the para-nitroaniline molecule and of a porphyrin derived push–pull molecule, showing good agreement with available experimental data. © 2002 American Institute of Physics. 关DOI: 10.1063/1.1516805兴

I. INTRODUCTION

The interest in density-functional theory 共DFT兲 has in- creased steadily since Hohenberg and Kohn¹ provided its theoretical justification and Kohn and Sham²its practical for- mulation for calculations. The independent-particle Kohn–

Sham model, which describes correlation in terms of a unique, albeit unknown, functional of the electron density, has had an enormous impact on computational chemistry, and the available approximate exchange-correlation func- tionals have provided a level of accuracy that is adequate in many situations.

The generalization of the Kohn–Sham equations by Runge and Gross³ represents an important milestone in the development of time-dependent DFT 共TDDFT兲. Thus, al- though most applications of DFT have focused on structure and static properties, the interest in time-dependent response properties has grown considerably in recent years.^{4 – 43}At the same time one has witnessed a great success of ab initio wave function quantum chemistry in computing molecular properties for nonvariational and variational wave functions using time-independent analytic gradient theory and time- dependent response theory. Recently, time-dependent proper- ties have been obtained also from the quasienergy ansatz, where properties of any type 共linear, nonlinear, static, dy- namic, magnetic, electric, internal, external兲 are obtained from the solution of a common set of equations, without further parametrizations or manipulations than those applied to the reference wave function. Considering the broad range of systems that nowadays can be treated at an acceptable level of accuracy by DFT, the implementation of such a

‘‘toolbox approach’’ for DFT calculations of molecular prop- erties would have broad ramifications. An implementation should then involve all the common forms of the exchange- correlation functional, also those that include some propor- tion of the Hartree–Fock exchange. In this paper, such a program is carried out by invoking two time-dependent variation principles: one based on the Ehrenfest theorem⁴⁴ and the other based on the quasienergy ansatz,^45,46where a parametrized time-evolution operator acts on a reference state, here the Kohn–Sham determinant. When the electron density is calculated as an expectation value of the density operator, all the changes in the density to a given order in the perturbation are obtained as the terms of a Baker–

Campbell–Hausdorff expansion.

To the best of our knowledge, the implementation pre- sented here is the first complete implementation of DFT frequency-dependent nonlinear properties beyond the local- density approximation 共LDA兲. Thus, in the previous imple- mentation of frequency-dependent hyperpolarizabilities by Gisbergen et al.,¹⁸the calculations are carried out at the adia- batic LDA共ALDA兲level, in which the response of the den- sity is always obtained from LDA functional derivatives, even when the density itself has been calculated in the generalized-gradient approximation 共GGA兲. By contrast, in our implementation, the response of the density is always obtained using the stated DFT functional. We do invoke the adiabatic approximation, however, meaning that the time- dependent functional is assumed to depend on the time- dependent density in the same manner that the time- independent functional depends on the time-independent density.

9630

0021-9606/2002/117(21)/9630/16/$19.00 © 2002 American Institute of Physics

This paper is organized as follows: In Sec. II we present DFT linear and quadratic response functions—in the Ehren- fest form in Sec. II E, and in the quasienergy form in Sec.

II F. Our computer implementation is outlined in Sec. III. In Sec. IV, we present two applications to illustrate the useful- ness of the theory. The first application concerns the static and dynamic hyperpolarizabilities of para-nitroaniline with a detailed analysis of the basis-set and functional dependence of the results; the second application addresses the second harmonic generation coefficient for a porphyrin derived push–pull molecule, which has the potential to be an effec- tive material for this particular property. The paper is sum- marized in Sec. V. Finally, in the Appendices, we provide the explicit formulas for the functional derivatives, together with a variable transformation that has been useful in deriving them.

II. THEORY

A. Linear and quadratic response functions

Response theory is perturbation theory with emphasis on properties rather than on states and energies, aiming to de- scribe how the properties respond to internal or external per- turbations. Consider the expectation value of a time- independent operator Aˆ , which for the unperturbed system is given by 具⁰兩Aˆ兩0典. In the presence of a time-dependent per- turbation Vˆ (t), the expectation value of Aˆ becomes time- dependent and may furthermore be expanded in powers of the perturbation,

具^t兩Aˆ兩t典⫽具^t兩Aˆ兩t典⁽⁰⁾⫹具^t兩Aˆ兩t典⁽¹⁾⫹具^t兩Aˆ兩t典⁽²⁾⫹¯. 共1兲 Each term in this expansion has a Fourier representation of the form

具^t兩Aˆ兩t典⁽¹⁾⫽

冕

^{具具Aˆ ;Vˆ␻1典典␻1exp共⫺i␻1t兲d␻1, 共2兲}

具^t兩Aˆ兩t典⁽²⁾⫽1

2

冕 冕

^{具具Aˆ ;Vˆ␻1,Vˆ␻2典典␻1,␻2}

⫻exp关⫺i共␻1⫹␻2兲t兴d␻1d␻2, 共3兲 and so on, where Vˆ^␻ is the Fourier transform of Vˆ (t),

Vˆ共t兲⫽

冕

⫺⬁

⬁

Vˆ^␻exp共⫺i␻^t兲d␻^. 共4兲 In practice, we work with monochromatic perturbations; the integrals in Eqs.共2兲–共4兲are then replaced by discrete sums.

In this paper, we shall derive and implement expressions for the linear and quadratic response functions 具具^{Aˆ ;Vˆ␻1}典典␻₁

and 具具^{Aˆ ;Vˆ␻1,Vˆ␻2}典典␻1,␻2 within the framework of Kohn–

Sham DFT. In this particular formulation of DFT, the elec- tron density ␳(r,t) is parametrized in terms of a time- dependent reference determinant 兩t典. There are thus no problems associated with the calculation of expectation val- ues of one-electron operators共even differential ones兲, as as- sumed in Eq.共1兲. We note that, in this paper, we shall only consider properties of spin-restricted closed-shell states asso-

ciated with real singlet operators; the generalization to imaginary and triplet operators is fairly straightforward and will be treated in a subsequent paper.

B. Kohn–Sham theory

In Kohn–Sham theory, we assume that the time- dependent density ␳(r,t) is represented in terms of a time- dependent reference Slater determinant兩t典. The Kohn–Sham energy is then written as a functional of this density in the following manner:

E关␳^,t兴⫽T_s关␳兴⫹V_ext关␳^,t兴⫹J关␳兴⫹E_xc关␳兴⫹V_NN. 共5兲 The first term is the kinetic energy evaluated as an expecta- tion value

T_s关␳兴⫽⫺1

2

兺

_i ^{具t兩ⵜi2兩t典, 共6兲}

whereas the second and third terms represent, respectively, the classical Coulomb interactions of the electron density with the external potential and with itself,

V_ext关␳^,t兴⫽V_N关␳兴⫹V关␳^,t兴

⫽

冕

^␳共r,t兲

冉 ^兺

^{I 兩r⫺ZRI I兩⫹v共r,t兲}

冊

^d␶,

J关␳兴⫽1

2

冕 冕

^{␳共r1,tr兲}12^{␳共r2,t兲}

d␶1d␶2. 共7兲 In Eq. 共7兲, we have split the external potential into a static nuclear potential and an explicitly time-dependent perturba- tion. The exchange-correlation functional E_xc关␳兴 in Eq.共5兲 contains all two-electron interactions except the Hartree term J关␳兴, that is, it includes the effects of exchange and correla- tion. In addition, it corrects for the error made in the evalu- ation of the kinetic energy according to Eq.共6兲. The last term in Eq. 共5兲represents the classical nuclear–nuclear repulsion energy.

In the widely used adiabatic approximation, which we have adopted here, the time-dependence of the exchange- correlation energy is contained in the density, that is, the exchange-correlation functional is approximated using the same functional form in the time-dependent and time- independent cases. It is not obvious that this approximation holds for other than slowly varying external fields, but it has been verified that the adiabatic approximation is adequate for calculating excitation energies.¹²

In Kohn–Sham theory, the time evolution of the spin orbitals is governed by the differential equation

关f共r₁,t兲⫹v共r₁,t兲兴␾j共r₁,t兲⫽id␾j共r₁,t兲

dt , 共8兲

where we have introduced the Kohn–Sham operator, f共r₁,t兲⫽h共r₁兲⫹j共r₁,t兲⫹vxc共r₁,t兲, 共9兲 which we choose to define without the explicit perturbation term 共it still depends implicitly on the perturbation through the density兲. The first term in Eq. 共9兲 contains the one- electron combined kinetic and nuclear-attraction operators,

h共r₁兲⫽⫺1 2ⵜ1

2⫹

兺

_{I 兩}r₁⫺^ZIR_I兩, 共10兲

whereas the last two terms are the Coulomb and exchange- correlation potentials, respectively,

j共r₁,t兲⫽

冕

兩^␳r^共₁⫺^r2r^,t₂^兲兩d␶2, 共11兲 vxc共r₁,t兲⫽ ␦^Exc

␦␳共r₁兲

冏

_␳(r

1)⫽␳(r₁,t)

. 共12兲

Note that these terms are themselves functionals of the den- sity.

C. Time evolution of the reference determinant and the density

The Kohn–Sham reference determinant satisfies the equation

共H⫹V兲兩t典^⫽i_dt^d兩t典^{; H⫽}

兺

_i ^{f共ri,t兲, 共13兲}

where the total Kohn–Sham Hamiltonian has been parti- tioned into a term H implicitly dependent on the perturbation and the perturbation itself V 关cf. Eq. 共8兲兴. In our second- quantization formulation of time-dependent Kohn–Sham theory, we adopt an exponential parametrization of the time- evolution operator, representing the time-dependent refer- ence Kohn–Sham determinant in the following manner:

兩t典⫽exp关⫺␬^ˆ共t兲兴兩0典^. 共14兲

Here兩0典 is the unperturbed Kohn–Sham determinant,

兩0典^⫽

兿

_k,␴ ^{ak†␴兩vac典 共15兲}

and␬ˆ (t) the anti-Hermitian operator,

␬^ˆ共t兲⫽

兺

_rs ^{␬rs共t兲Ers⫽}

兺

_rs ^␬rs共t兲

兺

_␴ ^{ar†␴as␴;}

␬rs*_共t兲⫽⫺␬sr共t兲, 共16兲 where a_r^†_␴ and a_s␴ are the creation and annihilation opera- tors, respectively, of orbitals r and s of spin␴. This particu- lar representation of the Kohn–Sham determinant has previ- ously been used in connection with relativistic four- component theory.⁴⁷

The fundamental variational parameters of our theory are the elements of the rotation matrix␬rs(t). As in Hartree–

Fock theory, the nonredundant rotations are those between occupied and unoccupied orbitals. Equation共14兲implies that the individual Kohn–Sham spin–orbitals obey the transfor- mation law,

a_k^†_␴共t兲⫽exp关⫺␬^ˆ共t兲兴a_k^†_␴exp关␬^ˆ共t兲兴 共17兲 and that orthonormality is preserved.

To study the time-dependence of the electron density, we introduce the second-quantized density operator,

␳^ˆ共r兲⫽

兺

pq ␾*p_共r兲␾q共r兲E_pq. 共18兲 Recognizing that the density is an expectation value of this operator and given the parametrization in Eq.共14兲, we obtain the following expression for the time-dependent density:

␳共r,t兲⫽具^t兩␳^ˆ共r兲兩t典⫽具⁰兩exp关␬^ˆ共t兲兴␳^ˆ共r兲exp关⫺␬^ˆ共t兲兴兩0典^. 共19兲 A Baker–Campbell–Hausdorff expansion of the exponential time-evolution operator gives for the density 共and similarly for other operators兲,

␳共r,t兲⫽␳共r,0兲⫹具^0兩关␬^ˆ共t兲,␳^ˆ共r兲兴兩0典

⫹¹2具⁰兩关␬^ˆ共t兲,关␬^ˆ共t兲,␳^ˆ共r兴兴兩0典^⫹O共␬³兲, 共20兲 which will be used frequently in the following sections.

D. Perturbation expansion of the density and Kohn–Sham operators

In the absence of the perturbation, the time-evolution operator produces only a dynamical phase factor, which can- cels out for expectation values. In the presence of the pertur- bation, it is assumed that an expansion in powers of the perturbation exists such that

␬^ˆ共t兲⫽␬^ˆ(1)共t兲⫹␬^ˆ(2)共t兲⫹¯ 共21兲 with Fourier representations

␬^ˆ(1)共t兲⫽

冕

^{␬ˆ␻1exp共⫺i␻1t兲d␻1, 共22兲}

␬^ˆ(2)共t兲⫽1

2

冕 冕

^{␬ˆ␻1,␻2exp关⫺i共␻1⫹␻2兲t兴d␻1d␻2.}

共23兲 For monochromatic perturbations, the integrations are re- placed by summations. Using the expansions Eqs. 共20兲 and 共21兲, we introduce the perturbed density matrices up to sec- ond order as

D_pq⁽⁰⁾⫽具^0兩Epq兩0典^{, 共24兲} D_pq⁽¹⁾⫽具^0兩关␬^ˆ(1),Epq兴兩0典^{, 共25兲} D_pq⁽²⁾⫽具^0兩关␬^ˆ(2),Epq兴⫹¹2关␬^ˆ(1),关␬^ˆ(1),Epq兴兴兩0典^{, 共26兲} which allow us to write the nth order correction to the den- sity as

␳⁽ⁿ⁾共r,t兲⫽

兺

_pq ^{␾*p共r兲␾q共r兲D(n)pq共t兲. 共27兲}

We now expand the second-quantized Kohn–Sham Hamiltonian in orders of the perturbation,

Hˆ⫽

兺

n

Hˆ⁽ⁿ⁾, 共28兲

Hˆ⁽ⁿ⁾⫽

兺

_pq ^{f(n)pqEpq, 共29兲}

where

f_pq⁽ⁿ⁾⫽␦n0h_pq⫹j_pq⁽ⁿ⁾⫹vxc, pq

(n) . 共30兲

The first contribution to f⁽⁰⁾_pq is the one-electron integral over the kinetic-energy and nuclear-attraction operators of Eq.

共10兲,

h_pq⫽

冓

^␾p

^冏

^⫺12ⵜ2⫹

^兺

^{I 兩r⫺ZIRI}

^冏

^␾q

冔

^{. 共31兲}

The electron-repulsion nth order Coulomb interaction inte- grals

j_pq⁽ⁿ⁾⫽

兺

rs

g_pqrsD_rs⁽ⁿ⁾ 共32兲

constitute the second contribution to the Kohn–Sham matrix element f_pq⁽ⁿ⁾ and are obtained from the two-electron inte- grals,

g_pqrs⫽

冓

^{␾p共1兲␾r共2兲}

^冏

^r112

^冏

^{␾q共2兲␾s共2兲}

冔

^共33兲

and from the time-dependent density-matrix elements in Eqs.

共24兲–共26兲. The last contribution to Eq. 共30兲 is the integral over the exchange-correlation potentialvxc关␳兴,

v_{xc, pq}⁽ⁿ⁾ 共t兲⫽具␾p兩v_xc⁽ⁿ⁾共r,t兲兩␾q典^{. 共34兲} Since the potential depends on the density ␳, this contribu- tion to the Kohn–Sham matrix element depends on the per- turbation,

v_xc⁽¹⁾共r,t兲⫽

冕

^␦␦␳^vxc共r^共⬘^r兲^兲␳⁽¹⁾共r⬘^,t兲d␶⬘^, 共35兲

v_xc⁽²⁾共r,t兲⫽

冕

^␦␦␳^vxc共r^共⬘^r兲^兲␳⁽²⁾共r⬘^,t兲d␶⬘

⫹1

2

冕 冕

␦␳共^␦r²⬘^v兲^xc␦␳^共r共^兲r⬙兲␳⁽¹⁾共r⬘^,t兲␳⁽¹⁾

⫻共r⬙^,t兲d␶⬘^d␶⬙^. 共36兲

E. The Ehrenfest method

Let us consider the time-development of the expectation value of an operator of the form,

Qˆ共t兲⫽exp关⫺␬^ˆ共t兲兴Qˆ exp关␬^ˆ共t兲兴. 共37兲 Differentiating the expectation value 具^{t兩Qˆ (t)兩t}典^⫽具^0兩Qˆ兩0典 关see Eq. 共14兲兴 with respect to time and invoking the Schro¨- dinger equation共atomic units兲

冉

^{Hˆ共t兲⫹Vˆ共t兲⫺idtd}

冊

^{兩t典⫽0 共38兲}

we obtain the Ehrenfest theorem

具^t兩iQˆ˙共t兲⫹关Qˆ共t兲,Hˆ共t兲⫹Vˆ共t兲兴兩t典⫽0, 共39兲 which forms the basis for the time-dependent variation principle.⁴⁴Collecting the exponential operators in Eqs.共14兲 and共37兲and using the identity

exp关␬^ˆ共t兲兴

冉

^{dtd exp关⫺␬ˆ共t兲兴Qˆ exp关␬ˆ共t兲兴}

冊

^{exp关⫺␬ˆ共t兲兴}

⫽

冋

^{exp关␬ˆ共t兲兴}

冉

^{dtd exp关⫺␬ˆ共t兲兴}

冊

^,Qˆ

册

^{, 共40兲}

we find that the Ehrenfest theorem may be written in the more convenient form

冓

⁰

^{冏 冋}

^{Qˆ ,exp关␬ˆ共t兲兴}

^冉

^{Hˆ共t兲⫹Vˆ共t兲⫺idtd}

^冊

^{exp关⫺␬ˆ共t兲兴}

^{册 冏}

⁰

冔

^⫽0.

共41兲 This equation holds for any time-independent one-electron operator Qˆ . In particular, it holds for the spin-averaged ex- citation operators E_pq in the expansion of ␬ˆ (t) in Eq. 共16兲. Collecting these operators in the column vector qˆ, we arrive at a set of nonlinear equations from which the time- dependence of␬ˆ (t) may be determined. In the following, we shall use these equations to determine the first- and second- order terms in Eq.共21兲and thereby the linear and quadratic response functions.

1. Linear response

The linear response equations are obtained by expanding Eq.共41兲to first order, yielding a differential equation in time for ␬^ˆ(1),

具⁰兩关qˆ,关␬^ˆ(1),Hˆ(0)兴⫹Hˆ⁽¹⁾兴兩0典⫹i具⁰兩关qˆ,␬^ˆ˙(1)兴兩0典

⫽⫺具⁰兩关qˆ,Vˆ共t兲兴兩0典^. 共42兲 In the frequency domain, this equation becomes an algebraic equation for␬^ˆ␻,

具⁰兩关qˆ,关␬^ˆ␻,Hˆ(0)兴⫹Hˆ␻兴兩0典⫹␻具⁰兩关qˆ,␬^ˆ␻兴兩0典

⫽⫺具⁰兩关qˆ,Vˆ␻兴兩0典^{. 共43兲} We have here introduced the Fourier transform of the per- turbed Kohn–Sham matrix,

Hˆ␻⫽

兺

_pq ^{f␻pqEpq, 共44兲}

where

f_pq^␻⫽j_pq^␻⫹v_{xc, pq}^␻

⫽

兺

rs

g_pqrsD_rs^␻

⫹

冓

^␾p

^{冏 冕}

^{␦␳␦v共rxc⬘兲␳␻共r⬘兲d␶⬘}

^冏

^␾q

冔

^{, 共45兲}

D_rs^␻⫽具⁰兩关␬^ˆ␻,Ers兴兩0典^, 共46兲

␳^␻共r兲⫽具⁰兩关␬^ˆ␻,␳^ˆ共r兲兴兩0典^. 共47兲 To bring out the matrix structure of Eq.共43兲, we may express the first-order parameters in matrix form,

␬^ˆ␻⫽¯E_pq¯

冉

^{␬⯗⯗pq␻}

冊

^{⬅qˆ†␬␻ 共48兲}

such that we formally solve the linear equations

共E⫺␻^S兲␬^␻⫽V^␻. 共49兲

In Eq. 共49兲, the matrix E has the structure of an electronic Hessian,

E⫽具^0兩关qˆ,关Hˆ⁽⁰⁾,qˆ^†兴兩0典^⫹_pqrs

兺

^{具0兩关qˆ,Epq兴兩0典}

冋

^兩r⫺1r⬘兩

⫹␦^vxc共r兲

␦␳共r⬘兲

册

pqrs

具⁰兩关E_rs,qˆ^†兴兩0典^{, 共50兲} where, in the last term, we have introduced two-electron in- tegrals defined by analogy with Eq. 共33兲. In Eq. 共49兲, we have introduced the generalized overlap matrix,

S⫽具⁰兩关qˆ,qˆ^†兴兩0典^, 共51兲

and the perturbation vector on the right-hand side of Eq.

共43兲,

V^␻⫽具^{0兩关qˆ,V␻兴兩0}典^{. 共52兲} From the solution ␬^{ˆ␻ to Eq.} 共49兲, the linear response function is easily obtained for any one-electron operator Aˆ ,

具具^{Aˆ ;Vˆ}典典␻⫽具⁰兩关␬^ˆ␻,Aˆ兴兩0典⫽⫺A^†共E⫺␻^S兲^⫺1V^␻, 共53兲 where

A⫽具⁰兩关qˆ,Aˆ兴兩0典^{. 共54兲}

It should be noted that the working equations of linear response theory resemble Eq.共43兲more closely than Eq.共49兲 since E and S are never constructed explicitly and since, in the iterative algorithms we use, ␬ corresponds to trial vec- tors. Furthermore, the exchange-correlation potentials that we consider are local and never give rise to two-electron integrals—even the GGA functionals are local in this sense.

2. Quadratic response

The quadratic response to a perturbation requires the so- lution of ␬^ˆ(2) from the second-order expansion of Eq.共41兲, 具^0兩关qˆ,关␬^ˆ(2),Hˆ(0)兴⫹Hˆ⁽²⁾兴兩0典^⫹i具^0兩关qˆ,␬^ˆ˙(2)

⫹¹2关␬^ˆ(1),␬^ˆ˙(1)兴兴兩0典

⫽⫺具⁰兩关qˆ,关␬^ˆ(1),Vˆ共t兲⫹Hˆ⁽¹⁾兴⫹¹2关␬^ˆ(1),关␬^ˆ(1),Hˆ(0)兴兴兴兩0典^. 共55兲 With the Fourier transform Eq. 共23兲, this equation can be rearranged to give

具⁰兩关qˆ,关␬^{ˆ␻1,␻2,Hˆ(0)}兴⫹Hˆ^␻1,␻2兴兩0典

⫹共␻1⫹␻2兲具⁰兩关qˆ,␬^ˆ␻1,␻2兴兩0典

⫽⫺Pˆ

12具⁰兩关qˆ,2关␬^{ˆ␻1,Vˆ␻2}⫹Hˆ␻2兴⫹关␬^ˆ␻1,关␬^{ˆ␻2,Hˆ(0)}兴兴

⫹␻2关␬^ˆ␻1,␬^ˆ␻2兴兴兩0典^. 共56兲 Here the operator Pˆ

12 symmetrizes with respect to the fre- quencies ␻1 and␻2,

Pˆ₁₂A共␻1,␻2兲⫽¹2关A共␻1,␻2兲⫹A共␻2,␻1兲兴 共57兲

and the elements of the second-order perturbed Kohn–Sham matrix,

Hˆ^␻1,␻2⫽

兺

_pq ^{f␻1,␻2Epq 共58兲}

are given by

f_pq^␻1,␻2⫽j_pq^␻1,␻2⫹v_{xc, pq}^␻1,␻2 共59兲 with Coulomb and exchange-correlation parts,

j^␻_pq^1,␻2⫽

兺

_rs ^{gpqrsDrs␻1,␻2, 共60兲}

v_{xc, pq}^␻1,␻2⫽

冓

^␾p

^{冏 冕}

^{␦␳␦共vrxc⬘兲␳␻1,␻2共r⬘兲d␶⬘}

^冏

^␾q

冔

⫹Pˆ

12

冓

^␾p

^{冏 冕 冕}

^{␦␳共r␦⬘2兲v␦␳xc共r⬙兲␳␻1共r⬘兲␳␻2}

⫻共r⬙兲d␶⬘^d␶⬙

冏

^␾q

冔

^{, 共61兲}

The solution of Eq. 共56兲 requires the separation of the second-order density matrix elements into first- and second- order parts,

␳^␻1,␻2共r兲⫽

兺

_pq ^{␾*p共r兲␾q共r兲D␻pq1,␻2, 共62兲}

D_rs^␻1,␻2⫽D¯

rs

␻1,␻2⫹Dញ

rs

␻1,␻2, 共63兲

D¯

rs

␻1,␻2⫽Pˆ₁₂具⁰兩关␬^ˆ␻1,关␬^ˆ␻2,Ers兴兴兩0典^, 共64兲 Dញ_rs^␻1,␻2_⫽具⁰兩关␬^{ˆ␻1,␻2,E}rs兴兩0典^, 共65兲 such that the second-order Hamiltonian may be written as the sum of two operators that depend on the first- and second- order parameters, respectively,

Hˆ^␻1,␻2⫽HC^␻1,␻2⫹Hញ^ˆ␻1,␻2. 共66兲 As in the linear case, the equation for the second-order pa- rameters Eq.共56兲may then written in matrix form,

关E⫺共␻1⫹␻2兲S兴␬^␻1,␻2⫽V^␻1,␻2, 共67兲 where all first-order parameters have been collected on the right-hand side,

V^␻1,␻2⫽2 Pˆ

12兵具⁰兩关qˆ,关␬^ˆ␻1,关␬^{ˆ␻2,Hˆ(0)}兴兴兴兩0典

⫹具⁰兩关qˆ,关␬^ˆ␻1,␻2␬^ˆ␻2兴兴兩0典⫹2具⁰兩关qˆ,关␬^␻1,Hˆ␻2

⫹V^␻2兴兴兩0典^⫹具⁰兩关qˆ,HC^␻1,␻2兴兩0典其^{. 共68兲} The quadratic response function may now be calculated as

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

⫹Pˆ

12具⁰兩关␬^ˆ␻1,关␬^ˆ␻2,Aˆ兴兴兩0典^, 共69兲 where the first-order responses ␬^{ˆ␻1 and} ␬^ˆ␻2 are obtained from Eq. 共43兲 and the second-order response ␬^{ˆ␻1,␻2 from} Eq.共56兲. We note that the first term in Eq.共69兲resembles the evaluation of the linear response function,

具⁰兩关␬^␻1,␻2Aˆ兴兩0典⫽⫺A^†␬^␻1,␻2

⫽⫺A^†关E⫺共␻1⫹␻2兲S兴^⫺1V^␻1,␻2. 共70兲 An alternative approach is to solve the adjoint linear re- sponse equation for A at frequency␻1⫹␻2, as done in our implementation.

In short, for a given property Aˆ and two periodic pertur- bations Bˆ and Cˆ with associated frequencies ␻b and␻c re- spectively, the evaluation of the quadratic response equation 具具^{Aˆ ;Bˆ;Cˆ}典典␻_b,␻_c is carried out by first solving the three linear response equations,

␬^A†关E⫺共␻b⫹␻c兲S兴⫽A^†, 共71兲

共E⫺␻bS兲␬^␻b⫽B, 共72兲

共E⫺␻cS兲␬^␻c⫽C, 共73兲

and then evaluating

具具^{Aˆ ;Bˆ;Cˆ}典典␻_b,␻_c⫽␬^{A†V␻b,␻c}⫹Pˆ_bc具^0兩关␬^ˆ␻b,关␬^ˆ␻c,Aˆ兴兴兩0典^. 共74兲 The main complication associated with extending a Hartree–

Fock quadratic response code to the Kohn–Sham DFT is the evaluation of the exchange-correlation contribution to V^␻b,␻c. Section III is devoted to this implementation in de- tail.

F. The quasienergy method

The quasienergy formulation provides an alternative time-dependent variation principle for deriving dynamical response functions. It is arguably more attractive than the Ehrenfest method in that it provides a unified framework for treating variational and nonvariational wave functions, by analogy with time-independent theory, to which it naturally reduces in the limit of static perturbations.^45,46As an addi- tional advantage over the Ehrenfest method, the permuta- tional symmetries共with respect to the exchange of operators兲 becomes manifest in the quasienergy method.

The concept of quasienergy arises naturally in TDDFT—

indeed, it turns out to be nothing but the action integral of Runge and Gross,³where the integration limits are chosen to span a period of the perturbation, scaled by the inverse of the period length. In the following, we discuss the application of the quasienergy method to the calculation of response func- tions in Kohn–Sham DFT. Since the Kohn–Sham quasien- ergy method follows closely the corresponding Hartree–

Fock method, we here consider only those aspects of Kohn–

Sham theory that differ from Hartree–Fock theory.

In Hartree–Fock theory, the quasienergy and its time av- erage are defined as follows:

Q共t兲⫽

冓

^t

^冏

^Hˆ⫺idtd

^冏

^t

冔

^{, 共75兲}

兵^Q其T⫽1 T

冕

⫺T/2

T/2

Q共t兲dt, 共76兲

where T constitutes one period of the time-dependent pertur- bation. In Kohn–Sham theory, there is an additional contri- bution to the quasienergy from the exchange-correlation functional,

Q_xc共t兲⫽E_xc关␳共t兲兴. 共77兲 The perturbation is expanded in its 共discrete兲 Fourier com- ponents, each of which is further expanded in field strengths

⑀_x^␻k coupled to a quantum-mechanical operator,

Vˆ共t兲⫽

兺

_kx ^{exp共⫺i␻kt兲⑀x␻kVˆx␻k. 共78兲}

Included in the summation is the operator Aˆ , the response functions of which we are calculating. We associate a non- zero frequency ␻ but zero perturbation strength with this operator.

The periodic perturbation Eq. 共78兲 induces a change in the quasienergy, which is expanded in orders of the pertur- bation,

Q共t兲⫽Q⁽⁰⁾共t兲⫹Q⁽¹⁾共t兲⫹Q⁽²⁾共t兲⫹¯. 共79兲 The response function of order n is then recovered by differ- entiating the time-averaged perturbed quasienergy Q^(n⫹1)(t) with respect to the frequency-dependent field strengths,

具具^{Aˆ ;Bˆ}典典␻⫽d²兵^Q(2)其^T

d⑀A⫺␻d⑀B␻, 共80兲

具具^{Aˆ ;Bˆ,Cˆ}典典␻₁,␻₂⫽ d³兵^Q(3)其^T

d⑀_A^{⫺␻1⫺␻2d}⑀_B^␻1d⑀_C^{␻2. 共81兲} The frequency associated with Aˆ is set equal to minus the sum of the perturbing frequencies so that the time-averaged quasienergy does not vanish.

To calculate the perturbed quasienergy Eq.共79兲, we must first determine the perturbed orbital-rotation parameters

␬ˆ (t), which are likewise expanded in orders of the perturba- tion,

␬^ˆ(1)共t兲⫽

兺

_kx ^{exp共⫺i␻kt兲⑀x␻k␬ˆx␻k, 共82兲}

␬^ˆ(2)共t兲⫽1

2_klxy

兺

^{exp关⫺i共␻k⫹␻l兲t兴⑀x␻k⑀y␻l␬ˆx␻k␬ˆy␻l. 共83兲}

From the 2n⫹1 rule, it follows that, to calculate the quasien- ergy to order 2n⫹1, we need only determine␬ˆ (t) to order n, the contributions from higher orders being zero. Thus, to calculate the linear and quadratic response functions, we need only determine␬^ˆ(1). Note that, even though the contri- bution of for example␬^{ˆ(2) to Q(2)}vanishes, its contribution to the exchange-correlation part of the quasienergy Q_xc⁽²⁾ is nonzero but cancelled by a similar contribution to Q⁽²⁾

⫺Q_xc⁽²⁾. In our discussion of Q_xc⁽²⁾and Q_xc⁽³⁾we shall therefore ignore all contributions that do not depend on␬^ˆ(1), in accor- dance with the 2n⫹1 rule.

1. Linear response

In linear response theory, the exchange-correlation en- ergy is expanded to second order in␬^ˆ(1),

E¯

xc

(2)⫽

冕

␦␳^␦E共r^xc兲^¯␳⁽²⁾共r,t兲 d␶⫹1

2

冕 冕

␦␳共^␦r兲²␦␳^Exc共r⬘兲

⫻␳⁽¹⁾共r,t兲␳⁽¹⁾共r⬘^,t兲 d␶ ^d␶⬘^. 共84兲 The perturbed densities are given by Eq.共27兲except that the contribution from ␬^{ˆ(2) to}␳⁽²⁾ is ignored in accordance with the 2n⫹1 rule, as indicated by the overbar in E¯

xc

(2)and^¯␳^(2). Next, we expand ␬^ˆ(1)(t) according to Eq.共82兲and obtain

E¯

xc (2)⫽1

2

兺

_kl ^{exp关⫺i共␻k⫹␻l兲t兴}

兺

_xy ^{⑀x␻k⑀␻yl}

⫻

冋 ^冕

^{␦␳␦E(r)xc¯␳␻1,␻2(r) d␶⫹}

^{冕 冕}

^{␦␳(r)␦2␦␳Exc共r⬘兲}

⫻␳^␻1(r)␳^␻2(r⬘^{) d}␶ ^d␶

册

^{, 共85兲}

where we have introduced

␳x␻⫽具^0兩关␬^ˆx␻,␳^ˆ兴兩0典^{, 共86兲}

␳

¯_xy^␻1,␻2⫽Pˆ₁₂具⁰兩关␬^ˆ_x^␻1,关␬^ˆ_y^␻2,␳^ˆ兴兴兩0典^. 共87兲 The first-order parameters␬_x^{␻1 and}␬_y^␻2 are determined by a variational condition on the second-order quasienergy, the result of which is an equation identical to Eq. 共43兲. Differ- entiating the time-averaged exchange-correlation quasien- ergy with respect to the field strengths⑀_a^␻a,⑀_b^␻b, according to Eq. 共80兲, we obtain

具具^{Aˆ ;Bˆ}典典␻xc_b⫽

冕

␦␳^␦E共r^xc兲␳^␻a,␻b共r兲 d␶⫹

冕 冕

␦␳共^␦r兲²␦␳^Exc共r⬘兲

⫻␳^␻a共r兲␳^␻b共r⬘兲 d␶ ^d␶⬘ 共88兲 for the exchange-correlation energy contribution to the linear response function.

2. Quadratic response

The quadratic response function is obtained as the third derivative of the time-averaged quasienergy. The program is then to expand the energy to third order in the first-order parameters,

E¯

xc

(3)⫽

冕

␦␳^␦E共r^xc兲^¯␳⁽³⁾共r,t兲 d␶

⫹

冕 冕

␦␳共^␦r兲²␦␳^Exc共r⬘兲␳⁽¹⁾共r,t兲␳^¯(2)共r⬘^,t兲d␶ ^d␶⬘

⫹1

6

冕 冕 冕

␦␳共r兲␦␳^␦3共^Er⬘^xc兲␦␳共r⬙兲

⫻␳⁽¹⁾共r,t兲␳⁽¹⁾共r⬘^,t兲␳⁽¹⁾共r⬙^,t兲 d␶ ^d␶⬘ ^d␶⬙^{, 共89兲} where

␳

¯⁽³⁾共r,t兲⫽¹6具⁰兩关␬^ˆ(1)共t兲,关␬^ˆ(1)共t兲,关␬^ˆ(1)共t兲,␳^ˆ共r兲兴兴兴兩0典^. 共90兲 Expanding the parameters as in the linear case Eq.共85兲, we obtain

E¯

xc (3)⫽1

6

兺

_klm ^{exp关⫺i共␻k⫹␻l⫹␻m兲t兴}

兺

_{x y z} ^{⑀x␻k⑀y␻l⑀z␻m}

⫻

冋 ^冕

^{␦␳␦E共xcr兲¯␳␻k,␻l,␻m共r兲 d␶}

⫹2 Pˆ

x y z

冕 冕

␦␳共^␦r兲²␦␳^Exc共r⬘兲␳^␻k共r兲^¯␳^␻l,␻m

⫻共r⬘兲d␶ ^d␶⬘⫹

冕 冕 冕

␦␳共r兲␦␳^␦3共^Er⬘^xc兲␦␳共r⬙兲␳^␻k

⫻共r兲␳^␻l共r⬘兲␳^␻l共r⬙兲d␶ ^d␶⬘ ^d␶⬙

册

^{, 共91兲}

where the first- and second-order perturbed densities are given by Eqs. 共86兲 and 共87兲 and the third-order perturbed density by

␳

¯_{x y z}^{␻1,␻2,␻3}共r兲⫽Pˆ₁₂₃具⁰兩关␬^ˆ_x^␻1,关␬^ˆ_y^␻2,关␬^ˆ_z^␻3,␳^ˆ共r兲兴兴兴兩0典^. 共92兲 As in the linear case, time averaging cancels the time- dependent phase factor if the sum of the frequencies is zero.

The exchange-correlation energy contribution to the qua- dratic response function then becomes

具具^{Aˆ ;Bˆ,Cˆ}典典_␻xc_b,␻_b⫽1

6

冋 ^冕

^{␦␳␦E共xcr兲␳␻a,␻b,␻c共r兲 d␶}

⫹3 Pˆ_ABC

冕 冕

_␦␳共^␦r兲²␦␳^Exc共r⬘兲

⫻␳^␻a共r兲␳^␻b␻c共r⬘兲 d␶ ^d␶⬘

⫹

冕 冕 冕

_␦␳共r兲␦␳^␦3共^Er⬘^xc兲␦␳共r⬙兲␳^␻a

⫻共r兲␳^␻b共r兲␳^␻c共r兲d␶ ^d␶⬘ ^d␶⬙

册

^{. 共93兲}

III. IMPLEMENTATION

The implementation of the response module consists of two major parts: The computation of the DFT response con- tribution expressed in terms of commutators of the density matrix and response vectors ␬ⁱ, and the computation of all 19 functional derivatives appearing in the formalism. Tech- nically, these two parts have very different computational characters. The functional-derivative evaluation involves a limited amount of data and a fixed number of floating point operations, including evaluation of power operators and tran- scendental functions共arcsin, etc.兲. The evaluation time for a single grid point may be substantial but does not depend on the system size. The DFT contribution is expressed as a sum of matrices. The evaluation of the prefactors and the sum of the matrices is a memory-intensive operation. The evaluation time of this part of the code per single numerical integration grid point scales formally as the squared number of basis functions but can be reduced to linear scaling by taking ad- vantage of the matrix sparsity.