On the Subsystem Formulation of Linear-Response Time-Dependent DFT
aa r X i v : . [ phy s i c s . c h e m - ph ] J un On the Subsystem Formulation ofLinear Response Time-DependentDFT
Michele Pavanello Department of Chemistry, Rutgers University,Newark, NJ 07102, USADate: October 16, 2018Status: Accepted by J. Chem. Phys. on May 3, 2013 E-mail: [email protected] bstract A new and thorough derivation of linear-response subsystem TD-DFT is presented and analyzed in detail. Two equivalent derivationsare presented and naturally yield self consistent subsystem TD-DFTequations. One derivation reduces to the subsystem TD-DFT for-malism of Neugebauer [J. Neugebauer,
J. Chem. Phys. , 134116(2007)]. The other yields Dyson type equations involving three typesof subsystem response functions: coupled, uncoupled and Kohn–Sham.The Dyson type equations for subsystem TD-DFT are derived herefor the first time.The response function formalism reveals previously hidden quali-ties and complications of the subsystem formulation of TD-DFT com-pared with the regular TD-DFT of the supersystem. For example,analysis of the pole structure of the subsystem response functionsshows that each function contains information about the electronicspectrum of the entire supersystem. In addition, comparison of thesubsystem and supersystem response functions shows that, while thecorrelated response is subsystem additive, the Kohn–Sham responseis not. Comparison with the non-subjective Partition DFT theoryshows that this non-additivity is largely an artifact introduced by thesubjective nature of the density partitioning in subsystem DFT. Introduction
When modeling systems that contain a large number of electrons, even theKohn–Sham Density Functional Theory (KS-DFT) approach [1] has its lim-its. In the past decade many approximations [2–4] made it possible to mas-sively reduce KS-DFT complexity for spatially extended molecules. However,the large pre-factor of such scaling laws left the calculation of most realistic,fully-solvated systems still prohibitive [5, 6].Reducing the computational complexity of KS-DFT by partitioning thetotal electron density of a system into subsystem contributions has been anappealing idea since the early works of Gordon and Kim [7, 8]. However, thesuccess of KS-DFT seemed to have rendered partitioning methods unneces-sary. This is evident from the Quantum Chemistry literature of the 70s and80s, where partitioning methods were frequent only to high-end wave func-tion methods, and interactions between subsystems were treated with varioustypes of perturbation theory [9]. Despite two successful applications of den-sity partitioning techniques, first by Senatore and Subbaswamy [10], and thenby Cortona [11], revival of these methods is due to a paper by Wesolowskiand Warshel published in 1993 [12]. Presently, subsystem DFT is being de-veloped by many research groups worldwide [13–21]. Successful applicationsof subsystem DFT are reported for applications related to the ground state,such as analysis of electron densities [22], and spin densities [23]; and forcalculations of charge and excitation energy transfer parameters [15, 24, 25];as well as for electronic spectra and molecular properties [20, 26–30].The time-dependent extension of susbsystem DFT has been pioneeredby Casida and Wesolowski [31]. However, Neugebauer [32] is credited forderiving working equations for the solution of the subsystem time-dependent3FT (subsystem TD-DFT, hereafter) and for applying subsystem TD-DFTto determine excitation energies [20, 28, 32], charge/exciton couplings [15, 24,25], and molecular properties [33]. Similarly to subsystem DFT, subsystemTD-DFT is developed to take full advantage of the subsystem nature ofthe majority of real life systems. Solvated systems are a typical exampleof this. An early success story of subsystem TD-DFT is the calculationof solvatochromic shifts [34]. More recently, the electronic spectra of lightharvesting complexes model systems containing more that 1000 atoms hasbeen calculated with this method [28].Linear-response TD-DFT has been formulated by many authors in manypublications. This has frequently offered a chance to discuss its limitationsand to offer possible solutions. This has not been the case for subsystemTD-DFT. Partly because this field is relatively new.The large body of work on subsystem DFT and TD-DFT shows theirusefulness and importance. However, a work that aims at clarifying therelationship between subsystem and supersystem TD-DFT, and analyzingwhat the density partitioning does to the collective time-dependent responseof the system is long overdue. This work, aims at filling this gap providingtwo derivations of subsystem TD-DFT. The new derivations are amenable toa deeper analysis and understanding of the theory of subsystem TD-DFT.For example, Dyson type equations for susbsystem TD-DFT are derived herefor the first time.This work is organized as follows. In the next section, a theoretical back-ground is given on KS-DFT, subsystem DFT, and linear-response TD-DFT.In Section 3, a rigorous derivation of linear-response subsystem TD-DFT iscarried out. In Section 4 an alternative derivation of subsystem TD-DFT ispresented in terms of subsystem response functions. Section 5 is devoted to4he comparison of subsystem TD-DFT with TD-DFT of the supersystem. InSection 6, conclusions are drawn.
KS-DFT can be summarized by the following equation, known as the KSequation, in canonical form, (cid:20) − ∇ + v eff ( r ) (cid:21) φ k ( r ) = ε k φ k ( r ) , (1)where v eff is the effective potential that the one-particle KS orbitals, φ k ,experience, and ε k are the KS orbital energies. The spin labels have beenomitted for sake of clarity, as throughout this work only the spin restrictedcase is considered without loss of generality of the derivations. The electrondensity is simply ρ ( r ) = 2 P occ i | φ i ( r ) | .The effective potential, v eff , is given by v eff ( r ) = v appl ( r ) + v eN ( r ) + v Coul ( r ) + v xc ( r ) , (2)with v appl being an externally applied potential, v eN the electron–nucleus at-traction potential, v Coul the Hartree potential, and v xc the exchange–correlation(XC) potential [1].Subsystem DFT is based on the idea that an electronic (molecular) sys-tem can be more easily approached if it is partitioned into many smallersubsystems. In mathematical terms, this is done by partitioning the electron5ensity as follows [10, 11] ρ ( r ) = N S X I ρ I ( r ) , (3)with N S being the total number of subsystems.The ultimate goal is to represent the subsystems as a set of N coupledKohn–Sham systems. Hence, the subsystem densities must be non-negative,must integrate to a preset number of electrons, i.e. R ρ I ( r )d r = N I , andmust be v -representable. In this context, it is perfectly legitimate to wonderwhat then constitutes a subsystem. The three requirements (constraints)mentioned above constitute the only theoretical prescription. It is remarkablethat this prescription does not invoke any real space partitioning.Therefore, subsystem densities can, in principle, strongly overlap and canbe highly delocalized. In practical calculations, however, the subsystem den-sities are constructed from subsystem molecular orbitals which are expandedin terms of localized atomic orbitals, often centered to atoms belonging toonly one molecular fragment in a system (monomer basis set). In practicalcalculations, the latter approximation and the use of local and semilocal non-additive kinetic energy functionals define the subsystems as non-covalentlybound molecules.Self consistent solution of the following coupled KS-like equations (alsocalled KS equations with constrained electron density [35]) yield the set ofsubsystem KS orbitals, i.e. (cid:20) − ∇ + v I eff ( r ) (cid:21) φ Ik ( r ) = ε Ik φ Ik ( r ) , with I = 1 , . . . , N S (4)6ith the effective subsystem potential given by v I eff ( r ) = v appl ( r ) + v I eN ( r ) + v I Coul ( r ) + v I xc ( r ) | {z } same as regular KS − DFT + v I emb ( r ) . (5)In the above it is clear that if an applied potential, v appl , acts on the totalsystem, every subsystem will experience that same potential. In the so-calledFrozen Density Embedding (FDE) formulation of subsystem DFT [12,35], theunknown potential above, v emb , is called embedding potential and is givenby v I emb ( r ) = N S X J = I "Z ρ J ( r ′ ) | r − r ′ | d r ′ − X α ∈ J Z α | r − R α | ++ δT s [ ρ ] δρ ( r ) − δT s [ ρ I ] δρ I ( r ) + δE xc [ ρ ] δρ ( r ) − δE xc [ ρ I ] δρ I ( r ) . (6)Throughout this work, “subsystem DFT” is used as a synonym of FDE.The density of the supersystem is thus found using Eq.(3) and Eq.(4) as ρ ( r ) = 2 P N S I P occ I i (cid:12)(cid:12) φ Ii ( r ) (cid:12)(cid:12) . The time-dependent KS equation, (cid:20) − ∇ + v eff ( r , t ) (cid:21) φ k ( r , t ) = i ∂φ k ( r , t ) ∂t , (7)relates the time dependent KS orbitals, φ k ( r , t ), and the correlated den-sity, ρ ( r , t ) = 2 P occ k | φ k ( r , t ) | , with the externally applied, time-dependentperturbation, v appl ( r , t ). When starting from the ground state density, the7ime-dependent KS potential is defined according to Eq.(2) as v eff ( r , t ) = v appl ( r , t ) + v eN ( r ) + v Coul ( r , t ) + v xc ( r , t ) , (8)The applied potential constitutes the only time-dependent perturbation caus-ing the density ρ to become a time-dependent function [31, 36]. With theexception of v eN , which is considered static (the nuclei are assumed to bestill in the time the perturbation is applied), the other potential terms partof the effective KS potential are dependent on time, but only as a result ofthe perturbation v appl ( r , t ) through their density dependence.Linear-response TD-DFT is based on the assumption that the densityresponse to the external weak perturbation is given by the following linear-response integral equations [37] δρ ( r , t ) = Z χ ( r , r ′ , t − t ′ ) δv appl ( r ′ , t ′ )d r ′ d t ′ (9)= Z χ ( r , r ′ , t − t ′ ) δv eff ( r ′ , t ′ )d r ′ d t ′ , (10)where δv eff ( r ′ , t ) = δv appl ( r ′ , t ) + δv ind ( r ′ , t ) , (11)The induced potential, δv ind , is expressed in linear-response as well, namely δv ind ( r ′ , t ) = Z t ′ = tt ′ = t Z (cid:20) δ ( t − t ′ ) | r ′ − r ′′ | + δv xc ( r ′ , t ) δρ ( r ′′ , t ′ ) (cid:21) δρ ( r ′′ , t ′ )d r ′′ d t ′ . (12)The quantity δv xc ( r ′ ,t ) δρ ( r ′′ ,t ′ ) is called the XC kernel, f xc ( r ′ , r ′′ , t − t ′ ). The func-tions χ ( r , r ′ , t − t ′ ) and χ ( r , r ′ , t − t ′ ) are the correlated and the simplifiedKS response functions (or simply “correlated response” and “KS response”,respectively). Eq.(9) constitutes the definition of linear-response TD-DFT,8nd Eq.(10) derives from Eq.(9) from the Runge–Gross theorem (Theorem4 of Ref. [38]).As it is more convenient to write the working equations in the frequencydomain, by virtue of the convolution theorem, the above equation can berewritten as δv ind ( r ′ , ω ) = Z (cid:20) | r ′ − r ′′ | + f xc ( r , r ′ , ω ) (cid:21) δρ ( r ′′ , ω )d r ′′ (13)For practical calculations, Eq.(10) is the most important. This is because,in the adiabatic approximation, it involves quantities that can be extractedfrom the ground state KS system, such as the KS response function (givenhere in Fourier transform) χ ( r , r ′ , ω ) = occ X i virt X a ω ia ω ia − ω φ i ( r ) φ a ( r ) φ i ( r ′ ) φ a ( r ′ ) , (14)where φ i and φ a are occupied and virtual KS orbitals. In a simpler notation,omitting the integral signs and the variable dependence, practical calcula-tions of the density response are carried out by self consistently solving thefollowing [37] h(cid:0) χ (cid:1) − − f i δρ = δv appl , (15)with f ( r , r ′ , ω ) = 1 | r − r ′ | + f xc ( r , r ′ , ω ) . (16)As Eq.(15) must hold for any δv appl , comparison with Eq.(9) yields( χ ) − = (cid:0) χ (cid:1) − − f, (17)also known as the Dyson equation for the response function [37, 39, 40].9 Subsystem TD-DFT
This section is devoted to the derivation of linear-response subsystem TD-DFT. Even though this theory has been first derived by Neugebauer [32],here it is presented in a different mathematical formalism which makes useof subsystem response functions. The derivations and analyses presented inthis section are important as they pave the road to the formalism presentedthe subsequent sections.
Following the usual decomposition of the density change in subsystem TD-DFT [20, 25, 32], the total electron density change of the system, δρ , due toan external perturbation is given exactly by δρ ( r , t ) = N S X I δρ I ( r , t ) , (18)where δρ I is the density change of the single subsystem I . The subsystemdensity change is, in all respects, equivalent to a regular TD-DFT densitychange. For example, one could think that the single subsystem densitychanges, δρ I ( r , t ), may involve inter-subsystem charge transfer type changesso that they may not integrate to zero. In fact, R δρ I ( r , t )d r = 0 alwaysbecause one of the defining constraints of the subsystems is that they mustbe made of a fixed number of electrons. Charge transfer type excitationsare naturally accounted for in this theory, as no real-space constraints areimposed on the subsystem densities.Similarly to Eq.(11), let us consider the effective time-dependent pertur-bation on subsystem I , δv I eff ( r , t ), as being a functional of all the subsystem10ensities, and defined as follows δv I eff ( r , t ) = δv appl ( r , t ) | {z } perturbation + δv I ind ( r , t ) | {z } induced potentialon subsystem I . (19)Similarly to Eq.(5), in the above it is assumed that the applied potential actson the entire system and therefore it is the same applied potential, δv appl ,that interacts with all the subsystems.The induced potential, δv I ind , can be defined in terms of functional deriva-tives of the subsystem KS potential given in Eq.(5) [20, 31, 32]. Defining K IJ ( r , r ′ , t − t ′ ) = δv I ind ( r , t ) δρ J ( r ′ , t ′ ) , (20)expressing all quantities in Fourier transform, and applying the convolutiontheorem, we get δv I ind ( r , ω ) = N S X J Z K IJ ( r , r ′ , ω ) δρ J ( r ′ , ω )d r ′ (21) K IJ ( r , r ′ , ω ) = 1 | r − r ′ | + f xc ( r , r ′ , ω ) + f T ( r , r ′ , ω ) − f I T ( r , r ′ , ω ) δ IJ , (22)where the kinetic kernels, expressed in the time domain, are defined as f T ( r , r ′ , t − t ′ ) = δ T s [ ρ ] δρ ( r , t ) δρ ( r ′ , t ′ ) , (23) f I T ( r , r ′ , t − t ′ ) = δ T s [ ρ I ] δρ I ( r , t ) δρ I ( r ′ , t ′ ) . (24)Eq.(22) was derived in Ref. [32] and is found by noticing that δ T s [ ρ ] δρ ( r ) δρ J ( r ′ ) = δ T s [ ρ ] δρ ( r ) δρ ( r ′ ) after applying the chain rule. Taking the partial functional deriva-11ive with respect to a single subsystem density of functionals of the totalsupersystem density is equivalent to taking the derivative with respect tothe total density [31, 32].Similarly to Eq.(9) and Eq.(10), with the aid of the Runge-Gross theorem,the time-dependent subsystem density can be obtained self consistently as δρ I ( r , ω ) = Z χ cI ( r , r ′ , ω ) δv appl ( r ′ , ω )d r ′ (25)= Z χ I ( r , r ′ , ω ) δv I eff ( r ′ , ω )d r ′ , (26)where χ I is the KS response of the subsystem to the external perturbation, χ cI is the correlated “coupled” subsystem response function, and δv I eff is givenby Eqs.(19–21).The above equations hold a great deal of information, e.g. the subsystemtime-dependent density can be obtained from the simplified subsystem KSresponse function and the effective time-dependent potential. Eqs.(19–21)can be used in Eq.(26), yielding δρ I ( r , ω ) = Z χ I ( r , r ′ , ω ) δv appl ( r ′ , ω )d r ′ + Z χ I ( r , r ′ , ω ) X J K IJ ( r , r ′ , ω ) δρ J ( r ′′ , ω )d r ′ d r ′′ . (27)From the above equation it becomes clear that a subsystem density responseis coupled to the density responses of other subsystems through the inducedpotential [32]. This is a key piece of information, as the subsystem densityresponse will appear in the expressions of all the other subsystem densityresponses. This is a picture of dynamic coupling between subsystems thatreveals how the labeling of the subsystem time-dependent quantities is just a12ormality. This analysis uncovers the fact that the dynamic response of thesupersystem is collective and generally not subsystem additive.Grouping the terms in Eq.(27) involving δρ I on the lhs and expressing δρ I ( r , ω ) in terms of integrals of suitable Dirac deltas yields Z (cid:2) δ ( r − r ′ ) δ ( r ′ − r ′′ ) − χ I ( r , r ′ , ω ) K II ( r ′ , r ′′ , ω ) (cid:3) δρ I ( r ′′ , ω )d r ′ d r ′′ == Z χ I ( r , r ′ , ω ) δv appl ( r ′ , ω )d r ′ + Z χ I ( r , r ′ , ω ) X J = I K IJ ( r ′ , r ′′ , ω ) δρ J ( r ′′ , ω )d r ′ d r ′′ . (28)Eq.(28) can be rearranged by acting on the left by ( χ I ) − ( r ′′′ , r , ω ) and inte-grating over d r , using the relation R ( χ I ) − ( r ′′′ , r , ω ) χ I ( r , r ′ , ω )d r = δ ( r ′′′ − r ′ ), Z h(cid:0) χ I (cid:1) − ( r ′′′ , r ′ , ω ) δ ( r ′ − r ′′ ) − δ ( r ′′′ − r ′ ) K II ( r ′ , r ′′ , ω ) i δρ I ( r ′′ , ω )d r ′ d r ′′ == δv appl ( r ′′′ , ω ) + Z δ ( r ′′′ − r ′ ) X J = I K IJ ( r ′ , r ′′ , ω ) δρ J ( r ′′ , ω )d r ′ d r ′′ . (29)After integration over r ′ and substitution of r ′′′ → r and r ′′ → r ′ the followingis obtained Z h(cid:0) χ I (cid:1) − ( r , r ′ , ω ) − K II ( r , r ′ , ω ) i δρ I ( r ′ , ω )d r ′ = δv appl ( r , ω )+ (30)+ Z X J = I K IJ ( r , r ′ , ω ) δρ J ( r ′ , ω )d r ′ . We now define the inverse of the “uncoupled” subsystem response functionas ( χ uI ) − ( r , r ′ , ω ) = (cid:0) χ I (cid:1) − ( r , r ′ , ω ) − K II ( r , r ′ , ω ) . (31)13ealizing that Eq.(30) holds for every subsystem, the following N S × N S matrix vector equation can be formally constructed M δρ = δv appl , (32)with M = ( χ uI ) − − K . . . − K (cid:0) χ uN S (cid:1) − , (33)and δρ = δρ I ... δρ N S , (34)where K is the matrix composed of the K IJ kernels. If the matrix in Eq.(33)is invertible, then the poles of M − occur at the true excitation energies ofeach subsystem, and hence of the total supersystem. Eq.(32) can be con-sidered the subsystem DFT equivalent of Eq.(15). The matrix formulationabove yields the coupled subsystem response function defined in Eq.(26) as δρ I = (cid:0) M − (cid:1) II δv appl , (35)thus a formal relationship is χ cI = (cid:0) M − (cid:1) II , (36)and δρ = Tr (cid:2) M − (cid:3) δv appl . (37)14quations (35–37) are well suited to be used in practical calculations. Thisis because, in practice, all the operators (response functions and kernels) areexpressed in a matrix form. However, in practical calculations, frequencyindependent kernels (adiabatic approximation) are usually adopted.It would be very useful to express the above equations completely in termsof subsystem response functions eliminating the δρ and v appl dependence,as that would lead to a Dyson-type equations formalism relating the time-dependent correlated and KS response of the total system with the ones ofthe subsystems. The following section provides precisely such a derivation. Using the definition in Eq.(31), Eq.(30) can be rearranged as follows Z ( χ uI ) − ( r , r ′ , ω ) δρ I ( r ′ , ω )d r ′ = δv appl ( r , ω )+ Z X J = I K IJ ( r , r ′ , ω ) δρ J ( r ′ , ω )d r ′ , (38)and expressing the subsystem density changes in terms of the subsystemresponse functions and the applied potential [using Eq.(25)], the above canbe simplified to Z ( χ uI ) − ( r , r ′ , ω ) χ cI ( r ′ , r ′′ , ω ) δv appl ( r ′′ , ω )d r ′′ d r ′ == Z " δ ( r − r ′ ) δ ( r ′ − r ′′ ) + X J = I K IJ ( r , r ′ , ω ) χ cJ ( r ′ , r ′′ , ω ) δv appl ( r ′′ , ω )d r ′′ d r ′ . (39)15he above equation must hold for any δv appl ( r ′′ , ω ), and specifically for δv appl ( r ′′ , ω ) = δ ( r ′′ − ˜ r ) f ( ω ), where f ( ω ) is any non-zero function of thefrequency. Integration over d r ′′ and simplification of the f ( ω ) term yields Z ( χ uI ) − ( r , r ′ , ω ) χ cI ( r ′ , ˜ r , ω )d r ′ = δ ( r − ˜ r )+ Z X J = I K IJ ( r , r ′ , ω ) χ cJ ( r ′ , ˜ r , ω )d r ′ . (40)Thus, after applying χ uI ( r ′′ , r , ω ) and integration over d r , the following Dyson-type equation is obtained, in simplified notation, χ cI = χ uI + N S X J = I χ uI K IJ χ cJ . (41)The above equation provides a general Dyson equation relating the uncoupledand the coupled subsystem response functions, and where is is clear thatthe coupling between subsystem responses is mediated by the off-diagonalelements of the kernel matrix K which contains exchange-correlation termsas well as kinetic energy terms.Dyson equations for the response functions involving only the kernelsand the KS response functions are derived starting from Eq.(30) and read asfollows χ uI = χ I + χ I K II χ uI , (42) χ cI = χ I + χ I N S X J K IJ χ cJ . (43)Similarly to regular TD-DFT, this formulation shows that the uncoupledresponse in Eq.(42) is similar to the one of the isolated subsystem, albeit asmall correction in the kernel due to the second functional derivative of the16on-additive kinetic energy functional.From the above equations, it is evident that if the poles of the responsefunction of subsystem I are well separated from the ones of the other subsys-tem response functions, then the poles of each subsystem response containthe ones of all other subsystems. This is a particularly interesting result, asit shows that formally the correlated response function of a single subsystemcontains information about the electronic spectrum of the entire supersys-tem. Obviously, in the limit of infinitely separated subsystems, K IJ ( r , r ′ , ω )will be identically zero when r and r ′ span regions of space occupied by dif-ferent subsystems. Thus, the above observation needs to be taken with agrain of salt as it is valid only if the subsystems are spatially close to eachother.The limiting case of infinite subsystem separation seems to simplify theformalism introducing some degree of subsystem additivity. However, theapproximations (such as the adiabatic approximation) usually employed inpractical implementations of this theory will likely break down in this limitingcase. Retardation effects (finite speed of interactions between subsystems)are completely neglected in practice and it is expected that they will stronglyinfluence the subsystem dynamical coupling when the subsystems are sepa-rated by large distances.Another interesting outcome of this formalism is that when two sub-systems have poles at the same frequencies in the isolated case (or in theuncoupled case), then this degeneracy must disappear in the coupled caseotherwise the response function would feature an unphysical “double pole”.This implies that the above formalism is coherent with the existence of Davy-dov splittings in dimeric systems [32, 41].Even though Eq.(41) is aesthetically pleasing, it is not suitable for prac-17ical calculations. The scheme developed in Eqs.(33–37) is recovered byrewriting Eq.(41) as χ uI = χ uI " ( χ uI ) − χ cI − N S X J = I K IJ χ cJ , (44)which leads to 1 = M χ c , (45)where M is the same matrix defined in Eq.(33). Similarly as before, fromEq.(45), the poles of M − are also the poles of the coupled response func-tion. The subsystem TD-DFT equations derived by Neugebauer [32] arereadily recovered in this formalism by rewriting the above equation in termsof occupied-virtual KS orbital products and applying the adiabatic approx-imation (i.e. f xc ( r , r ′ , t − t ′ ) = f xc ( r , r ′ , t − t ′ ) δ ( t − t ′ ), and similarly for thekinetic energy kernels). Comparison of subsystem DFT with regular KS-DFT is straightforward. Insubsystem DFT one has to solve coupled KS-like equations, where the cou-pling term is conveniently expressed as a potential term, v emb , added to theKS effective potential of the isolated subsystem. This means that the twoformalisms involve similar algorithms for practical calculations.As it will be clear from the following derivations, this is not the casefor the time-dependent extensions. In the following, Dyson-type equationswill make it possible to directly compare subsystem DFT with the TD-DFT18f the supersystem. The correlated response of the supersystem has a sim-ple relationship to the subsystem correlated responses, taking the functionalderivative with respect to δv appl of both sides of Eq.(18) one obtains χ = N S X I χ I . (46)Conversely, due to the non-uniqueness of the density partitioning in Eq.(3),the “simplified” KS response of the supersystem has no “simple” relationshipwith the subsystem KS responses. In order to find a relationship between the subsystem and the supersystemKS response functions, let us manipulate Eq.(25) by inverting the subsystemresponse functions, one by one, in the following manner (cid:0) χ I (cid:1) − δρ I = δv I eff , (47)where we have omitted the integration symbols for sake of a lighter notation.An important difference between the induced potential used in subsystemTD-DFT defined in Eq.(22) and the one used in the TD-DFT of the super-system defined in Eq.(12) resides in the kinetic energy kernels. The two canbe related, δv I eff ( r , ω ) = δv eff ( r , ω ) + δv I T ( r , ω ) , (48)19here, using Eq.(21) and Eq.(22), we define the kinetic energy part of thesubsystem kernel as δv I T ( r , ω ) = N S X J Z (cid:0) f T ( r , r ′ ω ) − f I T ( r , r ′ ω ) (cid:1) δρ J ( r ′ , ω )d r ′ = Z f T ( r , r ′ , ω ) δρ ( r ′ , ω ) − Z f I T ( r , r ′ , ω ) δρ I ( r ′ , ω ) , (49)where Eq.(18) has been used for the first term of the rhs. Using Eq.(49) inEq.(47), we obtain h(cid:0) χ I (cid:1) − + f I T i δρ I = (cid:0) f T χ (cid:1) δv eff , (50)where the number 1 above is intended to be the identity in functional space,i.e. a Dirac delta in the position representation. Inverting the operator on thelhs of the above equation and summing over all the subsystems, we obtainthe following N S X I δρ I = δρ = N S X I (cid:26)h(cid:0) χ I (cid:1) − + f I T i − (cid:0) f T χ (cid:1)(cid:27) δv eff . (51)At this point there are several algebraically non-equivalent ways to proceed.Two routes are considered here: the first one leading to an exact expression,and the second one leading to expressions suited for approximations.20 .1.1 Exact expression Extracting χ from the braces of Eq.(51), and realizing that χ δv eff = δρδρ = N S X I (cid:26)h(cid:0) χ I (cid:1) − + f I T i − h(cid:0) χ (cid:1) − + f T i(cid:27) χ | {z } Identity operator δρ. (52)By defining (cid:0) χ TI (cid:1) − = h ( χ I ) − + f I T i , Eq.(52) leads to χ = " N S X I (cid:0) χ TI (cid:1) − − − f T . (53)It should be noted that the KS supersystem considered here is the trueKS supersystem. Other formulations of subsystem TD-DFT [31], instead,considered a supersystem treated with Thomas–Fermi theory. A first approximation can be reached directly from Eq.(53) in the limit ofvanishing kernels, namely (cid:0) χ (cid:1) − = N S X I (cid:0) χ I (cid:1) − . (54)However, a different approximation can be make by first taking the func-tional derivative with respect to δv eff on both sides of Eq.(51), namely χ = N S X I χ I (cid:2) f I T χ I (cid:3) − (cid:0) f T χ (cid:1) , (55)21hich can be arranged to χ = [1 − f T S ] − S, (56)with S = P N S I χ I (cid:2) f I T χ I (cid:3) − . The above inverse operations expression canbe approximated with linear expansions, in the limit of small f I T χ I and small f T χ I , to χ ≃ N S X I χ I − N S X I χ I f I T χ I + N S X IJ χ I f T χ J , (57)featuring an interesting resemblance to the Dyson equation for the responsefunction. The derivations in the preceding section stand out as being too complicatedfor just the KS response, paradoxically in this context known as the “simpli-fied” response. What is the significance of such a complicated relationshipbetween the supersystem and the subsystem KS response functions? Whatis puzzeling in Eqs.(53–57) is that the KS response of the supersystem con-tains terms coupling KS responses of different subsystems. This is not avery good property of this theory, as subsystem additivity is sought in thedensity, in the correlated response in Eq.(46) and it is expected to appear inthe KS density response as well.This apparent artifact is due to the non-uniqueness and subjectivity ofthe density partitioning employed in Eq.(3). An indication of this artificialbehavior of the subsystem KS responses can be easily shown by consideringa more refined version of subsystem DFT known as partition DFT (PDFT).22n PDFT theory [42–45] the effective subsystem time-dependent potential is δv I eff = δv appl + δv ind + δv p , (58)where δv p is the change in the partition potential (a quantity shared byall subsystems and thus unique ). The above equation can be rearrangedsimilarly to the step carried out between Eq.(50) and Eq.(51), to yield h(cid:0) χ I (cid:1) − − f p i δρ I = (cid:0) χ (cid:1) − δρ, (59)with f p = δv p δρ . The above equation is rearranged to χ = N S X I h(cid:0) χ I (cid:1) − − f p i − (60)which can be approximated assuming small f p χ I by χ ≃ N S X I χ I − N S X I χ I f I p χ I . (61)The last two equations feature no cross terms coupling the subsystem KSresponses. Thus, PDFT provides a more intuitive time-dependent behaviorof the subsystems and is completely free of artifacts due to the non-uniquepartitioning appearing in regular subsystem DFT.Non-orthogonality also plays a role. For example, a subsystem additiveKS response function is expected to be a good approximation to the su-persystem KS response function in two limit cases: small electron densityoverlap between subsystems, and orthogonality between orbitals belongingto different subsystems. This is because in these cases, the non-additive ki-23etic energy functional is close to being identically zero and the treatmentbecomes similar to the PDFT case. In this work, the theory of linear-response subsystem TD-DFT is derivedin a complete way and analyzed in detail. For the first time, Dyson equa-tions involving subsystem response functions are derived for linear-responsesubsystem TD-DFT. Three types of subsystem response functions are con-sidered: coupled, uncoupled and KS. The coupled and uncoupled are exactand approximated correlated subsystem responses, respectively.It is found that, for non-infinitely separated subsystems, the pole struc-ture of a correlated (coupled) subsystem response function contains the ex-citations of the entire supersystem. This shows that if an applied potentialis in resonance with an electronic transition of one subsystem, the electronicresponse of another subsystem will also be strongly affected. This behav-ior generally does not fit a picture of “localized excitations” but instead isconsistent with the idea that the response of a collection of subsystems iscollective, and generally delocalized.Localization of the excitations may take place whenever the kernel cou-pling the subsystem’s excitations K IJ is small, which is often the case be-cause in practical calculations the subsystems are chosen to be non-bondedmolecules. However, a local picture of the time-dependent response of asystem is not generally accurate.The formalism presented here, specifically the Dyson equations, shows aremarkable similarity with the set of coupled equations one needs to solve forwhen considering a molecule interacting with a polarizable force field [46], or24he recent model for including non-local correlation in DFT by Tkatchenkoand coworkers [47]. In the two cases mentioned, the local polarizabilitiesare affected by the presence of other polarizabilities centered on differentatoms through Dyson equations, in all respects, similar to Eq.(41). Interest-ingly, the Hamiltonian that couples polarizabilities is the dipole Hamiltonian,whereas here the full Coulomb kernel is considered.Another interesting aspect that was uncovered in this work is that, whilethe correlated response of the supersystem is given by a simple sum of subsys-tem response functions (subsystem additive), the KS response is not. The ki-netic energy kernels are responsible for this non-additivity. Series expansionsreveal that the non-additive portions include terms coupling KS responsesof different subsystems. This non-additivity is also a feature of another par-titioning technique, PDFT. However, in PDFT the non-additive terms donot couple KS responses of different subsystems. This indicates that theunwanted cross-subsystem non-additivity occurring in subsystem TD-DFTis entirely an artifact stemming from the subjective nature of the densitypartitioning. Acknowledgements
I thank Johannes Neugebauer, Neepa Maitra, Ruslan Kevorkyants, and HenkEshuis for illuminating discussions. This work is supported by startup fundsof the Department of Chemistry and the office of the Dean of FASN of RutgersUniversity-Newark. 25 eferences [1] W. Kohn and L. J. Sham, Phys. Rev. , 1133 (1965).[2] G. E. Scuseria, J. Phys. Chem. A , 4782 (1999).[3] S. Goedecker, Rev. Mod. Phys. , 1085 (1999).[4] D. R. Bowler and T. Miyazaki, Rep. Prog. Phys. , 036503 (2012).[5] P. Carloni, U. Rothlisberger, and M. Parrinello, Acc. Chem. Res. ,455 (2002).[6] K. Burke, J. Chem. Phys. , 150901 (2012).[7] R. G. Gordon and Y. S. Kim, J. Chem. Phys. , 3122 (1972).[8] Y. S. Kim and R. G. Gordon, J. Chem. Phys. , 1842 (1974).[9] B. Jeziorski and W. Ko los, Perturbation approach to the study of weakintermolecular interactions, in Molecular Interactions, Vol. 3 , edited byH. Ratajczak and W. J. Orville-Thomas, pages 1–46, Wiley, Chichester,1982.[10] G. Senatore and K. R. Subbaswamy, Phys. Rev. B , 5754 (1986).[11] P. Cortona, Phys. Rev. B , 8454 (1991).[12] T. A. Wesolowski and A. Warshel, J. Phys. Chem. , 8050 (1993).[13] P. de Silva and T. A. Wesolowski, J. Chem. Phys. , 094110 (2012).[14] X. Hu, Y. Jin, X. Zeng, H. Hu, and W. Yang, Phys. Chem. Chem. Phys. , 7700 (2012). 2615] M. Pavanello, T. Van Voorhis, L. Visscher, and J. Neugebauer, J. Chem.Phys. (2013), Submitted, available at http://arxiv.org/abs/1211.4880.[16] A. S. P. Gomes and C. R. Jacob, Annu. Rep. Prog. Chem., Sect. C:Phys. Chem. , 222 (2012).[17] S. H¨ofener, A. S. P. Gomes, and L. Visscher, J. Chem. Phys. , 044104(2012).[18] S. Laricchia, E. Fabiano, and F. D. Sala, J. Chem. Phys. , 014102(2012).[19] J. D. Goodpaster, T. A. Barnes, and T. F. Miller, III, J. Chem. Phys. , 164108 (2011).[20] J. Neugebauer, Phys. Rep. , 1 (2010).[21] M. Iannuzzi, B. Kirchner, and J. Hutter, Chem. Phys. Lett. , 16(2006).[22] S. Fux, K. Kiewisch, C. R. Jacob, J. Neugebauer, and M. Reiher, Chem.Phys. Lett. , 353 (2008).[23] A. Solovyeva, M. Pavanello, and J. Neugebauer, J. Chem. Phys. ,194104 (2012).[24] M. Pavanello and J. Neugebauer, J. Chem. Phys. , 234103 (2011).[25] J. Neugebauer, C. Curutchet, A. Munoz-Losa, and B. Mennucci, J.Chem. Theory Comput. , 1843 (2010).[26] A. S. P. Pawel Tecmer, Henk van Lingen and L. Visscher, J. Chem.Phys. , 084308 (2012). 2727] C. K¨onig, N. Schl¨uter, and J. Neugebauer, J. Chem. Phys. (2012),Accepted.[28] C. K¨onig and J. Neugebauer, Phys. Chem. Chem. Phys. , 10475(2011).[29] R. E. Bulo, C. R. Jacob, and L. Visscher, J. Phys. Chem. A , 2640(2008).[30] C. R. Jacob and L. Visscher, J. Chem. Phys. , 194104 (2006).[31] M. E. Casida and T. A. Wesolowski, Int. J. Quantum Chem. , 577(2004).[32] J. Neugebauer, J. Chem. Phys. , 134116 (2007).[33] J. Neugebauer, J. Chem. Phys. , 084104 (2009).[34] J. Neugebauer, M. J. Louwerse, E. J. Baerends, and T. A. Wesolowski,J. Chem. Phys. , 094115 (2005).[35] T. A. Wesolowski, One-electron Equations for Embedded Electron Den-sity: Challenge for Theory and Practical Payoffs in Multi-Level Mod-eling of Complex Polyatomic Systems, in Computational Chemistry:Reviews of Current Trends , edited by J. Leszczynski, volume 10, pages1–82, World Scientific, Singapore, 2006.[36] M. E. Casida and M. Huix-Rotllant, Annu. Rev. Phys. Chem. , 287(2012).[37] M. E. Casida, Time-Dependent Density Functional Response Theoryfor Molecules, in Recent Advances in Density Functional Methods Part , edited by D. P. Chong, pages 155–192, World Scientific, Singapore,1995.[38] E. Runge and E. K. U. Gross, Phys. Rev. Lett. , 997 (1984).[39] E. K. U. Gross, J. F. Dobson, and M. Petersilka, Density FunctionalTheory of Time-Dependent Phenomena, in Density Functional TheoryII , edited by J. D. Dunitz, K. Hafner, K. N. Houk, S. Ito, J.-M. Lehn,K. N. Raymond, C. W. Rees, J. Thiem, and F. V¨ogtle, Topics in CurrentChemistry, pages 81–172, Springer, 1996.[40] M. Petersilka, U. J. Gossmann, and E. K. U. Gross, Phys. Rev. Lett. , 1212 (1996).[41] A. S. Davydov, Theory of Molecular Excitons , McGraw-Hill, New York,1962.[42] R. Tang, J. Nafziger, and A. Wasserman, Phys. Chem. Chem. Phys. ,7780 (2012).[43] P. Elliott, K. Burke, M. H. Cohen, and A. Wasserman, Phys. Rev. A , 024501 (2010).[44] M. H. Cohen, A. Wasserman, R. Car, and K. Burke, J. Phys. Chem. A , 2183 (2009).[45] P. Elliott, M. H. Cohen, A. Wasserman, and K. Burke, J. Chem. TheoryComput. , 827 (2009).[46] B. Mennucci, Phys. Chem. Chem. Phys. , 6583-6594 (2013).[47] A. Tkatchenko, R. A. DiStasio, R. Car and M. Scheffler, Phys. Rev.Lett.108