Hartree–Fock and Kohn–Sham atomic-orbital based time-dependent response theory

Hartree–Fock and Kohn–Sham atomic-orbital based time-dependent response theory

Helena Larsen,^a)Poul Jo”rgensen, and Jeppe Olsen

Department of Chemistry, University of Aarhus, Langelandsgade 140, 8000 A˚ rhus C, Denmark Trygve Helgaker^b)

Department of Chemistry, University of Cambridge, Lensfield Road, Cambridge, CB2 1EW, United Kingdom

共Received 31 July 2000; accepted 24 August 2000兲

A reformulation of general time-dependent Hartree–Fock and Kohn–Sham response theories that refers strictly to the atomic-orbital basis is presented. It is based on a recently proposed exponential parametrization of the one-electron atomic-orbital density matrix. In the presented formulation, only matrix multiplications and additions of sparse matrices are needed to compute the response functions and linear scaling with system size may, therefore, be obtained. Thus, this formalism is well suited to the computation of dynamic and static properties for large molecules at the Hartree–

Fock and Kohn–Sham density-functional levels of theory. © 2000 American Institute of Physics.

关S0021-9606共00兲30643-2兴

I. INTRODUCTION

Currently, there is an intense and multidisciplinary effort aimed at developing photonic materials and technologies.

Whereas electronic techniques are based on the transport of electrons, photonic techniques use photons for the transport and storage of information. The photonic techniques exploit nonlinear interactions between molecules and electromag- netic fields and computational chemistry may, therefore, con- tribute to the development of photonic materials by supply- ing accurate nonlinear susceptibilities. As the photonic materials typically are large organic molecules or polymers, it is important to develop computational methods that are able to compute nonlinear properties for such molecules. For reviews of nonlinear properties of large molecules, we refer to Ref. 1.

In response function theory we determine the time- development of an observable when the molecular system is subjected to, for example, an external electric or magnetic field. This field may oscillate with a given frequency that causes the wave function and the observed properties to be- come frequency dependent. It provides an efficient method for the calculation of response properties like nonlinear sus- ceptibilities for small and medium sized molecules and much effort has, therefore, been devoted to the development of this technique. In the present paper we present a novel method for the computation of response functions within the self- consistent field 共SCF兲 theories Hartree–Fock and Kohn–

Sham density-functional theory that may scale linearly with the size of the system and thus is suitable for the computa- tion of nonlinear properties of large molecules.

Much work has been done to develop density-functional and ab initio methods that are able to handle very large mo-

lecular systems. Ideally, these methods should scale linearly with the number of basis functions, N. Especially, the SCF theories have been considered, since these methods provide a good compromise between relatively low computational cost and reasonable accuracy—not only for the energy but also for molecular properties. The computational cost of the two- electron integrals is effectively reduced to scale linearly with the system size using the fast multipole method共FMM兲and integral prescreening techniques.^2–7 The N³ scaling of the diagonalization of the density matrix can be reduced using either localized orbitals or by optimizing the density matrix directly without constructing orbitals.^8–14In the methods de- veloped so far, the main concern has been how to optimize the total energy, although a few methods concerning re- sponse properties have already been suggested.¹⁵

The response functions are usually obtained in the molecular-orbital 共MO兲 basis.^16,17 However, since it is not the wave function correction but rather the expectation val- ues that represent the time development of the observables of interest, response functions may be determined directly in the atomic-orbital共AO兲basis. We describe how this may be done expressing response functions in terms of the Fock, overlap and density matrices in the AO basis. Matrix multi- plications and additions involving these sparse matrices al- low linear scaling to be achieved.

The density-based time-dependent response functions are derived exploiting the general exponential parametriza- tion of the AO density matrix introduced by Helgaker et al.

in Refs. 18 and 19. After a brief introduction to time- dependent response theory, we derive a density-based formu- lation of the time-dependent variation principle in Sec. III.

Next, in Sec. IV, we consider the exponential parametriza- tion of the density matrix in the time-dependent case, and in Secs. V–VIII we derive the response functions. Although the derivation allows the construction of response functions to arbitrary order, we restrict our discussion to the linear and

a兲Electronic mail: [email protected]

b兲Permanent address: Department of Chemistry, University of Oslo, P.O.B.

1033 Blindern, N-0315 Oslo, Norway.

8908

0021-9606/2000/113(20)/8908/10/$17.00 © 2000 American Institute of Physics

quadratic response functions. In Sec. IX, we discuss how to solve sets of linear equations needed to obtain the response functions. After a brief discussion of Kohn–Sham theory in Sec. X, we give some concluding remarks in Sec. XI.

II. INTRODUCTION TO TIME-DEPENDENT HARTREE–FOCK THEORY

Consider a molecular system described by a time- independent Hamiltonian, Hˆ

0. When a general field W(t) is applied to the system, the system interacts with the field. The interaction operator may be denoted Vˆ^t. We assume that W(t) vanishes at t⫽⫺⬁. The interaction operator then also vanishes at t⫽⫺⬁ and can be expressed as

Vˆ^t⫽

冕

⫺⬁

⬁

Vˆ␻exp关共⫺i␻⫹␧兲t兴d␻^, 共1兲 where␧ is a positive infinitesimal that ensures V^⫺⬁ is zero.

From the Hermiticity of V^t, it follows that:

共V^␻兲^†⫽V^⫺␻. 共2兲

The Hamiltonian, Hˆ , for the total system becomes

Hˆ⫽Hˆ₀⫹Vˆ^t. 共3兲

In Hartree–Fock theory, ⌿ is represented by a single- determinant wave function

⌿共x₁,x₂,...,x_M,t兲⫽共M !兲^⫺1/2det兩␾1共x₁,t兲␾2共x₂,t兲¯

⫻␾I共x_I,t兲¯␾M共x_M,t兲兩, 共4兲 where ␾I is the Ith occupied molecular spin orbital and where x_I denotes the space and spin coordinates of the Ith electron. For ease of notation, we suppress the explicit coor- dinate and time dependence of␾^and⌿. We assume that the molecular spin orbitals are orthonormal

具␾I兩␾J典^⫽␦IJ. 共5兲

The zero-order wave function is obtained at t⫽⫺⬁ where no perturbation is applied

⌿0⫽共M !兲^⫺1/2det兩␾1 0␾2

0¯␾M

0兩 共6兲

and is assumed to be optimized. This is ensured when the occupied molecular spin orbitals are eigenfunctions of the zero-order Fock operator²⁰

Fˆ

0␾I 0⫽␧I␾I

0. 共7兲

Here␾I

0 is a time-independent molecular spin orbital and␧I

its energy. Fˆ₀ is given by Fˆ

0⫽hˆ₀⫹

兺

_I ^{共JˆI0⫺KˆI0兲, 共8兲}

where hˆ₀ is the one-electron operator hˆ₀⫽⫺1

2ⵜ²⫺

兺

_{A 兩}R_A^Z⫺^Ar兩, 共9兲 and Jˆ_I⁰ and Kˆ

I

0 are the Coulomb and exchange operators, respectively,

Jˆ_I⁰␾J

0共x₁兲⫽

冕

^{共␾I0兲*共x2兲r1}12␾I

0共x₂兲dx₂␾J

0共x₁兲, 共10兲

Kˆ

I 0␾J

0共x₁兲⫽

冕

^{共␾I0兲*共x2兲r1}12␾J

0共x₂兲dx₂␾I

0共x₁兲. 共11兲 The integrations are over the space- and spin coordinates of electron 2.

The time development of the Slater determinant, Eq.共4兲, may be determined by one of the time-dependent variation principles as discussed in the next section. From the time development of the wave function, we may determine the time development of the expectation value of a given opera- tor, Aˆ . This time development is conveniently expressed in terms of response functions

具^Aˆ典⫽具⌿0兩Aˆ兩⌿0典⫹

冕

⫺⬁

⬁

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

⫹1 2

冕

⫺⬁

⬁

冕

⫺⬁

⬁

exp关⫺i共␻1⫹␻2兲t兴

⫻具具^{Aˆ ;Vˆ␻1,Vˆ␻2}典典␻₁,␻₂d␻2d␻1⫹¯, 共12兲 where 具具^{Aˆ ;Vˆ␻1}典典␻₁ and 具具^{Aˆ ;Vˆ␻1,Vˆ␻2}典典␻₁,␻₂ are the linear and quadratic response functions, respectively. The response functions are defined to be symmetric with respect to inter- change of the frequencies—for example,

具具^{Aˆ ;Vˆ␻1,Vˆ␻2}典典␻1,␻2⫽具具^{Aˆ ;Vˆ␻2,Vˆ␻1}典典␻2,␻1. 共13兲 In the following, we discuss how the response functions may be determined using an exponential parametrization of the density matrix in the AO basis.

III. DENSITY-BASED TIME-DEPENDENT HARTREE–FOCK THEORY

In this section, we derive a general equation for the time development of the density matrix. We begin by deriving the variation principle in the MO basis. Next, we consider a density-based formulation that refers strictly to the AO basis.

A. Time-dependent equations in the MO basis

The solutions to the time-dependent Schro¨dinger equa- tion depend upon an overall phase factor, which is redundant in the limit where Vˆ^t becomes static共time-independent兲. To obtain an approximate formulation that comprises both the static limit and the dynamic case, it is pertinent to eliminate the overall phase factor. This may be accomplished by using the Langhoff–Epstein–Karplus time-dependent variation principle in phase-isolated form²¹

Re具␦⌿兩

冉

^{Hˆ⫺i⳵⳵t}

冊

^{兩⌿典⫽0, 共14兲}

where ␦⌿ is an allowed variation in ⌿. Furthermore, the time-dependent wave function should fulfill the orthonormal- ization constraint

具␦⌿兩⌿典⫹具⌿兩␦⌿典⫽0. 共15兲

The phase of⌿ is not defined by Eq. 共15兲. However, since we express our results in terms of response functions 共or equivalently densities兲, which do not depend on the phase factor, it is not necessary to determine the phase explicitly.

Let us consider the terms in Eq.共14兲that contain Hˆ 具␦⌿兩Hˆ兩⌿典^⫹具^⌿兩Hˆ兩␦⌿典^⫽␦具^{⌿兩Hˆ兩⌿}典^{. 共16兲} Expanding in terms of orbitals, we obtain

␦具⌿兩Hˆ兩⌿典⫽

兺

_I ^{兵具␦␾I兩Fˆ⫹Vˆt兩␾I典⫹具␾I兩Fˆ⫹Vˆt兩␦␾I典其,}

共17兲 where Fˆ is the time-dependent Fock operator. It is obtained as a straightforward generalization of Fˆ₀ where the orbitals depend explicitly on time.

Next, we consider the terms in Eq.共14兲that contain the partial derivatives with respect to time. Since ⳵^/⳵^{t behaves} like a one-electron operator, we obtain

Re i具␦⌿兩⌿˙典⫽Re i

冉 ^兺

^{I 具␦␾I兩␾˙I典}

⫹

兺

_I⫽J^{具␦␾I兩␾I典具␾J兩␾˙J典}

冊

^{. 共18兲}

This expression may be simplified using the orthonormaliza- tion constraints on the MOs, Eq. 共5兲

具␦␾I兩␾I典⫹具␾I兩␦␾I典⫽2 Re具␦␾I兩␾I典⫽0, 共19兲 具␾^˙J兩␾J典⫹具␾J兩␾^˙J典⫽2 Re具␾^˙J兩␾J典⫽0. 共20兲 The terms 具␦␾I兩␾I典 ^and 具␾^˙J兩␾J典 are thus pure imaginary and the last sum in Eq. 共18兲therefore vanishes. Combining Eqs. 共17兲and共18兲, we obtain the time-dependent equations for the MOs

Re

兺

I 具␦␾I兩Fˆ⫹Vˆ^t⫺i ⳵

⳵^t兩␾I典⫽0. 共21兲 B. The time-dependent MO equations in matrix

representation

The molecular spin orbitals, ␾I, can be expanded in AOs, ␹␮

␾I⫽

兺

_␮ ^{C␮I␹␮. 共22兲}

Substitution in the first term of Eq.共21兲gives

Re

兺

_I ^{具␦␾I兩Fˆ兩␾I典⫽Re}

兺

_I ^{␦CI†fCI, 共23兲}

where C_I is a column vector containing the elements C_␮I from Eq. 共22兲, and f is the Fock matrix in the AO basis

f_␮␯⫽h_␮␯⫹

兺

_␳␴ ^{D␳␴共g␮␯␳␴⫺g␮␴␳␯兲, 共24兲}

where h_␮␯and g_␮␯␳␴ are given by

h_␮␯⫽

冕

^{␾␮*共x兲}

冉

^⫺12ⵜ2⫺

^兺

^{A 兩RAZ⫺Ar兩}

冊

^{␾␯共x兲dx, 共25兲}

g_␮␯␳␴⫽

冕 冕

^{␾␮*共x1兲␾␳*共x2r兲}12^{␾␯共x1兲␾␴共x2兲}

dx₁dx₂. 共26兲 The second term in Eq.共21兲is treated in the same way; for the last term, we obtain

Re i

兺

_I ^{具␦␾I兩}_⳵^⳵t兩␾I典^{⫽Re i}

兺

_I ^{␦CI†SC˙I, 共27兲}

where

S_␮␯⫽具␹␮兩␹␯典^{. 共28兲}

Thus, Eq.共21兲may be written as Re

兺

I ␦^CI

†

冉

^{f⫹Vt⫺iS⳵⳵t}

冊

^{CI⫽0. 共29兲}

The variations ␦^CI are not independent as the MOs must satisfy the orthonormality condition

C^†SC⫽1, 共30兲

which leads to the constraints

␦^CI

†SC_J⫹C_I^†S␦^CJ⫽0 共31兲

for the first-order variations. Introducing the Hermitian La- grange multipliers,␭IJ, Eqs.共29兲and共31兲may be combined to give the unconstrained Lagrangian equations

Re

兺

_I ^␦CI†

再 ^冉

^{f⫹Vt⫺iS⳵⳵t}

^冊

^CI⫺

^兺

^{J SCJ␭JI}

冎

^{⫽0. 共32兲}

As the variations␦^CIand␦^CI

†may be considered as linearly independent, Eq. 共32兲is satisfied when

冉

^{f⫹Vt⫺iS⳵⳵t}

冊

^CI⫺

^兺

^{J SCJ␭JI⫽0. 共33兲}

A similar equation is obtained for the complex conjugate.

Using standard matrix notation, Eq.共33兲becomes

冉

^{f⫹Vt⫺iS⳵⳵t}

冊

^{C⫽SC␭, 共34兲}

which is the form given in Ref. 20. Note carefully that C is a rectangular matrix containing in column I the expansion coefficients C_␮I of the Ith occupied orbital. In the time- independent limit, Eq.共34兲reduces to the well-known equa- tion fC⫽SC⑀^{, with}␭IJ⫽⑀IJ. Eq.共34兲determines the time- development of the transformation matrix C. Next, we will rewrite this equation to obtain an equation for the time- development of the density matrix.

C. Time-dependent equations for the AO density matrix

Multiplying Eq.共34兲with C^†S from the right, we obtain

再

^{共f⫹Vt兲C⫺iS⳵⳵tC}

冎

^{C†S⫽SC␭C†S. 共35兲}

Likewise, the complex conjugate equations can be rewritten to obtain the conjugate transpose of Eq.共35兲. Subtracting the two sets of equations, we obtain

共f⫹V^t兲CC^†S⫺SCC^†共f⫹V^t兲⫺i共SC˙ C^†S⫹SCC˙^†S兲⫽0.

共36兲 Introducing the one-electron matrix in the AO basis

D⫽CC^†, 共37兲

these equations can be written as

共f⫹V^t兲DS⫺SD共f⫹V^t兲⫽iSD˙ S. 共38兲 Note that the time-dependent variation principle in Eq. 共38兲 refers exclusively to the AO basis and that, in the time- independent limit, Eq.共38兲reduces to the standard SCF con- ditions

f₀DS⫺SDf₀⫽0, 共39兲 where f₀ is the time-independent Fock matrix in the AO basis. Furthermore, as Eq.共38兲defines the time development of the SCF density matrix, it constitutes the SCF Liouville equation in an orthonormal basis.

For future convenience, it is advantageous to rewrite Eq.

共38兲slightly. First note that, for all matrices M

TrE_␮␯M⫽M_␯␮, 共40兲

where E_␮␯ is a unit matrix with elements given by

关E_␮␯兴␳␴⫽␦␮␳␦␯␴. 共41兲

The (␮␯)th element of Eq.共38兲can, therefore, be written as TrE_␮␯兵^{共f⫹Vt兲DS⫺SD共f⫹Vt兲}其^⫽iTrE␮␯共SD˙ S兲. 共42兲 Introducing the S commutator

关A,B兴S⫽ASB⫺BSA, 共43兲

we may write Eq. 共42兲in a short-hand notation as

Tr共f⫹V^t兲关D,E_␮␯兴S⫽iTrE_␮␯共SD˙ S兲. 共44兲 This is the fundamental equation that we shall use to derive the response functions.

IV. PARAMETRIZATION OF THE DENSITY MATRIX In Refs. 18 and 19 we presented a general exponential parametrization of the density matrix in the AO basis. In this section, we briefly review this parametrization and consider how the time development may be introduced. The density matrices in the AO basis fulfill the symmetry, trace, and idempotency conditions, which in the spin–orbital basis are given by

D^†⫽D, 共45兲

TrDS⫽N_e, 共46兲

DSD⫽D. 共47兲

The density matrices in the AO basis are related by the transformation^18,19

D共X兲⫽exp共⫺XS兲D exp共SX兲, 共48兲 where X is an anti-Hermitian matrix. Note that D(X) fulfills Eqs.共45兲–共47兲. Finally, X should comply with the projection relation

X⫽P共X兲⫽PXQ^†⫹QXP^†, 共49兲 where P⫽DS and Q⫽1⫺DS, in order to eliminate redun- dancies, which are nonvanishing choices of X for which D(X)⫽D(0)⫽D. Note that the projectors are constructed from the matrix D(0). The above parametrization of the AO density matrix is equivalent to the one in the MO basis where similar symmetry, trace and idempotency conditions are ful- filled. Furthermore, Eq.共48兲reduces to a standard orthogonal transformation in the MO basis when S⫽1.¹⁸

The transformed density matrix, D(X), may be evalu- ated using the asymmetric Baker–Campbell–Hausdorff 共BCH兲expansion according to Ref. 18

D共X兲⫽D⫺关X,D兴S⫹¹2关X,关X,D兴S兴S⫺¯. 共50兲 Introducing the operator Xˆ by

Xˆ M⫽关X,M兴S, 共51兲

D(X) can be written in the more compact form as

D共X兲⫽_n

兺

_⫽^⬁₀ ^共⫺n!^1兲nXˆⁿD. 共52兲 The time development of the AO density matrix is intro- duced by X(t) which may be written as

X共t兲⫽_␮⬎␯

兺

^{共X␮␯共t兲E␮␯⫺X␮␯*共t兲E␯␮兲⫽}

兺

m

X_m共t兲O_m, 共53兲 where the diagonal elements, X_␮␮, vanish since we have used a phase-isolated variation principle. In Eq. 共53兲, we have introduced the short-hand notation

X_m共t兲⫽

再

^{⫺XX␮␯␮␯*共t共兲t兲 mm⬎⬍00 共54兲}

and similarly

O_m⫽

再

^{QQm†m⫽⫽EE␮␯␯␮ mm⬍⬎0.0 共55兲}

The time evolution of the density matrix is determined by Eq. 共44兲. Using Eqs.共53兲and共55兲, this equation may be expressed as

Tr共f共X兲⫹V^t兲关D共X兲,O_m^†兴S⫽iTrO_m^†SD˙共X兲S, 共56兲 which for t⫽⫺⬁ 共and X⫽0) becomes the standard SCF equation, Eq.共39兲. Note that f(X) depends on X through the density matrix

f_␮␯共X兲⫽h_␮␯⫹G_␮␯共D共X兲兲, 共57兲 where, for convenience, we have introduced

G_␮␯共D共X兲兲⫽

兺

_␳␴ ^{D␳␴共X兲共g␮␯␳␴⫺g␮␴␳␯兲. 共58兲}

V. GENERAL EQUATIONS FOR THE TIME DEVELOPMENT

The time development of the AO density matrix can be determined from Eq. 共56兲. Inserting the expression for the Fock matrix Eq. 共57兲into Eq.共56兲, we obtain

Trh关D共X兲,O_m^†兴S⫹TrG共D共X兲兲关D共X兲,O_m^†兴S

⫹TrV^t关D共X兲,O_m^†兴S⫽i TrO_m^†SD˙共X兲S. 共59兲 Let us first consider the evaluation of the left-hand side of Eq.共59兲. Inserting the expressions for D(X) and X from Eqs.

共52兲and共53兲, the first term can be written as

Trh关D共X兲,O_m^†兴S⫽_n

兺

_⫽^⬁₀ ^共⫺n!^1兲nTrh

⫻

冋冉

^␮⫽

^兿

^{n1 Oˆl␮}

冊

^D,Om†

册

S _␮⫽

兿

ⁿ₁ ^Xl␮共t兲, 共60兲 where summation over the repeated index l_␮ is understood 共the Einstein summation convention兲. Likewise, the second term becomes

TrG共D共X兲兲关D共X兲,O_m^†兴S

⫽_n

兺

_⫽^⬁_{0 k}

兺

_⫽ⁿ_{0 共}n^共⫺^⫺k¹兲^兲!k!ⁿ TrG

冉冉

^␮⫽n

^兿

^{⫺1k Oˆl␮}

冊

^D

冊

⫻

冋冉

^␯⫽

^兿

^{k1 Oˆj␯}

冊

^D,Om†

册

S

冉

^␮⫽n

^兿

^{⫺k1 Xl␮共t兲}

冊冉

^␯⫽

^兿

^{k1 Xj␯共t兲}

冊

^,

共61兲 where terms of identical order in X have been collected.

Defining

E_ml

1l₂¯l_n

关n⫹1兴 ⫽共⫺1兲ⁿ

n! Trh

冋冉

^␮⫽

^兿

^{n1 Oˆl␮}

冊

^D,Om†

册

S

⫹_k

兺

_⫽ⁿ_{0 共}n^共⫺⫺k¹兲^兲!k!ⁿ TrG

冉冉

^␮⫽n

^兿

^{⫺k1 Oˆl␮}

冊

^D

冊

⫻

冋冉

^␮⫽

^兿

^{n⫺nk⫹1 Oˆl␮}

冊

^D,Om†

册

S

, 共62兲

the two first terms of Eq.共59兲can be written as Trh关D共X兲,O_m^†兴S⫹TrG共D共X兲兲关D共X兲,O_m^†兴S

⫽_n

兺

_⫽^⬁₀ ^{Eml关n⫹1l12¯兴 ln}_␮⫽

兿

ⁿ₁ ^Xl␮共t兲. 共63兲

Writing out the expression for E_m^关1兴, we obtain

E_m^关1兴⫽Trf₀关D,O_m^†兴S⫽TrO_m^†共f₀DS⫺SDf₀兲. 共64兲 Thus, E_m^关1兴 is equal to a element of the gradient and vanishes for an optimized state. The elements of the second and third E matrices are given as 共assuming symmetrization of the indices兲

E_mn^关2兴⫽⫺Trf₀关关O_n,D兴S,O_m^†兴S

⫺TrG共关O_n,D兴S兲关D,O_m^†兴S, 共65兲

E_mnk^关3兴 ⫽¹2Trf₀关关O_k,关O_n,D兴S兴S,O_m^†兴S⫹TrG共关O_n,D兴S兲

⫻关关O_k,D兴S,O_m^†兴S⫹¹2TrG共关O_k,关O_n,D兴S兴S兲

⫻关D,O_m^†兴S. 共66兲 Similarly, the last term on the left-hand side of Eq.共59兲may be written as

TrV^t关D共X兲,O_m^†兴S

⫽n

兺

⫽0

⬁ 共⫺1兲ⁿ

n! TrV^t

冋冉

^␮⫽

^兿

^{n1 Oˆl␮}

冊

^D,Om†

册

S _␮⫽

兿

ⁿ₁ ^Xl␮共t兲

⫽_n

兺

_⫽^⬁₀ ^{Vmlt关n1⫹l21¯兴ln}_␮⫽

兿

ⁿ₁ ^Xl␮共t兲, 共67兲 where

V_ml

1l₂¯l_n

t关n⫹1兴 ⫽共⫺1兲ⁿ

n! TrV^t

冋冉

^␮⫽

^兿

^{n1 Oˆl␮}

冊

^D,Om†

册

S

. 共68兲

Writing Eq. 共68兲 out for the matrices V^t关1兴 and V^t关2兴, we obtain

V_m^t关1兴⫽TrV^t关D,O_m^†兴S⫽TrO_m^†共V^tDS⫺SDV^t兲, 共69兲 V_mn^t关2兴⫽⫺TrV^t关关O_n,D兴S,O_m^†兴S. 共70兲 Note that V_m^t关1兴 has the same structure as the gradient, Eq.

共64兲.

Finally, let us consider the term on the right-hand side of Eq. 共59兲

iTrO_m^†SD˙共X兲S⫽iTrO_m^†S

再

^{⳵⳵texp共⫺XS兲D exp共SX兲}

冎

^S.

共71兲 To simplify this equation, we use the relation

⳵

⳵^t兵^exp共⫺MS兲N exp共SM兲其

⫽_n

兺

_⫽^⬁₀ ^共⫺n!^1兲n共MˆⁿN˙兲

⫹n

兺

⫽0

⬁ k

兺

⫽0

n 共⫺1兲^k

k!共n⫺k⫹1兲!兵^Mˆk关N,共Mˆ^n⫺kM˙兲兴S其^{, 共72兲} derived in the Appendix. Noting that D˙⫽0, Eq.共71兲can be written as

iTrO_m^†SD˙共X兲S⫽i_n

兺

_⫽^⬁_{0 k}

兺

_⫽ⁿ₀ k!共n^共⫺^⫺k¹⫹^兲k1兲!

⫻TrO_m^†S兵^Xˆk关D,共Xˆ^n⫺kX˙兲兴S其^{S. 共73兲} Inserting the expression for X from Eq. 共53兲, we obtain

iTrO_m^†SD˙共X兲S

⫽in

兺

⫽0

⬁ k

兺

⫽0

n 共⫺1兲^k

k!共n⫺k⫹1兲!TrO_m^†S

再 ^冉

^␯⫽

^兿

^{k1 Oˆj␯}

^冊

⫻

冋

^D,

冉

^n␮⫽⫺

^兿

^{k⫹21 Oˆl␮}

冊

^Ol1

册

S

冎

^S

^冉

^␯⫽

^兿

^{k1 Xj␯共t兲}

^冊

⫻

冉

^n⫺␮⫽

^兿

^{k⫹21 Xl␮共t兲}

冊

^X˙l1共t兲

⫽i_n

兺

_⫽^⬁₁ ^{Sml关n⫹1l1兴2¯lnX˙l1共t兲}_␮⫽

兿

ⁿ₂ ^Xl␮共t兲, 共74兲 where

S_ml

1l₂¯l_n

关n⫹1兴 ⫽_kⁿ

兺

_⫽^⫺₀¹ k!^共⫺共n⫺^1兲k^k兲!TrO_m^†S

再 ^冉

^␮⫽n

^兿

^{⫺nk⫹1 Oˆl␮}

^冊

⫻

冋

^D,

冉

^␮⫽n

^兿

^{⫺k2 Oˆl␮}

冊

^Ol1

册

S

冎

^{S. 共75兲}

The explicit forms for the first two terms in Eq.共75兲become S_mn^关2兴⫽TrO_m^†S关D,O_n兴SS, 共76兲 S_mnk^关3兴 ⫽⫺TrO_m^†S关O_k,关D,O_n兴S兴SS

⫹¹2TrO_m^†S关D,关O_k,O_n兴S兴SS. 共77兲 Thus, Eq.共59兲may be written as

n

兺

⫽0

⬁

共E_ml

1l₂¯l_n 关n⫹1兴 ⫹V_ml

1l₂¯l_n

t关n⫹1兴 兲_␮⫽

兿

ⁿ₁ ^Xl␮共t兲

⫽i_n

兺

_⫽^⬁₁ ^{Sml关n⫹1l12¯兴 lnX˙l1共t兲}_␮⫽

兿

ⁿ₂ ^Xl␮共t兲, 共78兲 which is formally equivalent to the time-dependent equations previously derived for SCF and multiconfiguration SCF 共MCSCF兲 wave functions in the MO basis. The order-by- order solution of Eq. 共78兲 is, therefore, similar to the one described in Ref. 16. In the following, we summarize how the response equations and response functions are obtained.

VI. FIRST- AND SECOND-ORDER EQUATIONS

The set of parameters X(t) may be expanded in powers of the perturbation

X_n共t兲⫽X_n^共1兲共t兲⫹X_n^共2兲共t兲⫹¯, 共79兲 where the zero-order coefficient vanishes since the reference state is optimized. The parameters, X_n⁽ⁱ⁾(t), may be deter- mined by requiring Eq.共78兲to be valid to each order of the perturbation. Let us consider the expressions needed to de- termine X_n⁽¹⁾(t) and X_n⁽²⁾(t). The terms in Eq.共78兲that may contribute to the evaluation of these coefficients are E_mn^关2兴X_n共t兲⫹E_mnk^关3兴 X_n共t兲X_k共t兲⫺iS_mn^关2兴X˙

n共t兲

⫺iS_mnk^关3兴 X˙_n共t兲X_k共t兲⫹V_m^t关1兴⫹V_mn^t关2兴X_n共t兲⫹¯⫽0. 共80兲

As before summation is carried out over all repeated indices.

Inserting the expansion of X_n(t) from Eq.共79兲and collecting terms of first order, we get

E_mn^关2兴X_n^共1兲共t兲⫺iS_mn^关2兴X˙

n

共1兲共t兲⫽⫺V_m^t关1兴. 共81兲 The second-order equation is given by

E_mn^关2兴X_n^共2兲共t兲⫺iS_mn^关2兴X˙

n 共2兲共t兲

⫽⫺E_mnk^关3兴 X_n^共1兲共t兲X_k^共1兲共t兲⫹iS_mnk^关3兴 X˙

n

共1兲共t兲X_k^共1兲共t兲

⫺V_mn^t关2兴X_n^共1兲共t兲. 共82兲 To solve the first-order equation, Eq.共81兲, we use the Fourier expansion of X_n⁽¹⁾(t)

X_n^共1兲共t兲⫽

冕

^{exp共⫺i␻t兲Xn共1兲共␻兲d␻. 共83兲}

Inserting Eqs. 共1兲and共83兲into Eq.共81兲, we obtain

冕

^{exp共⫺i␻t兲兵共Emn关2兴⫺␻Smn关2兴兲Xn共1兲共␻兲⫹Vm␻关1兴其d␻⫽0,}

共84兲 which implies

X^共1兲共␻兲⫽⫺共E^关2兴⫺␻^S关2兴兲^⫺1V^␻关1兴. 共85兲 Likewise to solve the second-order equation, the Fourier transform

X_n^共2兲共t兲⫽

冕 冕

^{exp关i共␻1⫹␻2兲t兴Xn共2兲共␻1,␻2兲d␻1d␻2}

共86兲 is introduced. Insertion in Eq.共82兲gives

关E_mn^关2兴⫺共␻1⫹␻2兲S_mn^关2兴兴X_n^共2兲共␻1,␻2兲

⫽¹2共⫺E_mnk^关3兴 ⫺E_mkn^关3兴 ⫹␻1S_mnk^关3兴

⫹␻2S_mkn^关3兴 兲X_n^共1兲共␻1兲X_k^共1兲共␻2兲⫺V_mn^␻1关2兴X_n^共1兲共␻2兲

⫺V_mn^␻2关2兴X_n^共1兲共␻1兲, 共87兲 where X_n⁽²⁾(␻1,␻2) is defined to be symmetric in ␻1 and

␻2.

VII. THE STRUCTURE OF E^†2‡ANDS^†2‡

Consider Eq.共85兲in more detail. Ordering the matrices O_m in Eq. 共53兲 in the order 1,2,...,m,⫺1,⫺2,...,⫺m (m

⬎0) and assuming that integrals and density matrices are real, we can write E^关2兴 and S^关2兴 in the following forms:

E^关2兴⫽

冉

^{AB BA}

冊

^{, 共88兲}

S^关2兴⫽

冉

^{⫺⌺⌬ ⫺⌬⌺}

冊

^{, 共89兲}

where

A_mn⫽⫺Tr兵^f0关关Q_n^T,D兴S,Q_m兴S⫺G共关Q_n^T,D兴S兲关D,Q_m兴S其^, 共90兲