Analytical gradients for excitation energies from frozen-density embedding

The formulation of analytical excitation-energy gradients from time-dependent density functional theory within the frozen-density embedding framework is presented. In addition to a comprehensive mathematical derivation, we discuss details of the numerical implementation in the Slater-function based Amsterdam Density Functional (ADF) program. Particular emphasis is put on the consistency in the use of approximations for the evaluation of secondand third-order non-additive kinetic-energy and exchange–correlation functional derivatives appearing in the final expression for the excitationenergy gradient. We test the implementation for different chemical systems in which molecular excitedstate potential-energy curves are affected by another subsystem. It is demonstrated that the analytical implementation for the evaluation of excitation-energy gradients yields results in close agreement with data from numerical differentiation. In addition, we show that our analytical results are numerically more stable and thus preferable over the numerical ones.


Introduction
A full interpretation of excited-state processes requires to have reliable theoretical models for the explanation of experimental data in terms of excited-state geometric and electronic structure at hand. 1 The characterization of excited-state potential energy surfaces is important for understanding a variety of phenomena, such as photochemical reactions, vibrationally resolved UV-Vis spectra, resonance Raman spectra, excited-state dynamics etc. Efficient calculations of potential energy surfaces require an analytical approach for the evaluation of excited-state gradients. The main advantage of such an analytical calculation of gradients in contrast to a numerical approach is that the computational cost does not depend on the number of nuclear degrees of freedom. For a long time, the standard method to calculate nuclear derivatives of excited-state energies for large molecular systems has been the configuration interaction singles method. 2 In 1994, Ortiz derived and implemented gradients of excited states within the random phase approximation. 3 The first implementation of analytical excited-state gradients within Time-Dependent Density Functional Theory (TDDFT) was presented by Van Caillie and Amos. 4 While in that work the expression for the gradients was obtained by straightforward differentiation, Furche and Ahlrichs employed 5 a more convenient Lagrangianbased approach, similar to the one used in ref. 6 for configuration interaction methods. The same technique was used by Hutter in his work, 7 where the calculation of analytical excited-state gradients from the Tamm-Dancoff approximation to TDDFT was implemented within a plane-wave-basis set framework. In 2005, Doltsinis extended 8 plane-wave TDDFT excited-state gradient calculations beyond the Tamm-Dancoff approximation, where in addition, an active space approach 9 and an implicit differentiation scheme were employed. Seth et al. 10 implemented analytical gradients for spin-flip TDDFT excitation in the Amsterdam Density Functional program package (ADF), 11,12 which uses Slatertype functions in contrast to most other quantum chemistry packages. These implementations provided a powerful toolbox for the theoretical characterization of excited-state potentialenergy surfaces of large molecules. Embedding techniques like the Polarizable Continuum Model 13 (PCM), Quantum Mechanics/Molecular Mechanics 14 (QM/MM) or Frozen-Density Embedding (FDE) 15 can increase the range of application for TDDFT gradients. Analytical gradients for the evaluation of excitation-energy gradients have been implemented, e.g. for PCM-TDDFT 16 and QM/MM-TDDFT 17 and other fragmentation or embedding strategies. 18,19 However, no such analytical gradients are available so far for FDE, which partitions the total electron density of a large system into that of an active region and the frozen density of an environment. The TDDFT extension of FDE has been presented in ref. 20 and 21, and several challenging applications have been presented since then (see ref. [22][23][24][25] or the review in ref. 26 for some examples). Following the strict version of ''frozen-density embedding'', we will assume here and in the following that the response of the system to an external perturbation is entirely localized in the active system. This corresponds to the case of ''uncoupled FDE'' (FDEu) for excited states in the nomenclature of ref. 27, which is also known as the ''neglect of differential response of the environment'' (NDRE). 28 This can be relaxed in a full subsystem TDDFT treatment (see ref. 27 and the work by Casida and Wesolowski 20 ). For recent reviews, see ref. 29 and 30. Note that analytical ground-state energy gradients for FDE have been reported before, 31 and will thus not be considered here. The development of analytical excited-state gradient theory in TDDFT within FDE will give fast access to potential energy surfaces of the excited states of complex systems.

Theory
The excited-state gradient in TDDFT is calculated as the sum of the ground-state gradient and the gradient of the excitation energy, but since the former is established for FDE 31 we focus here only on the latter part. In order to introduce the notation, the underlying theory of FDE and its TDDFT extension are briefly presented in Sections 2.1 and 2.2, respectively. The formulation of FDEu-TDDFT excitation-energy gradients is concisely established in Section 2.3, which is followed by a derivation of the Z-vector equation for FDEu-TDDFT gradients in Section 2.4. An in-depth derivation including all intermediate steps is presented in Appendix A. Note that our theoretical development is based on the formalism introduced in the work by Seth et al. 10 as well as in that by Furche and Ahlrichs. 5

Frozen-density embedding
The FDE theory 15 is a fundamental framework for embedding schemes within DFT. The total density of the system in the FDE approach is partitioned into an active part and a frozen part r tot (r) = r(r) + r f (r). (1) The active electron density, r(r), can be expressed through the orbitals located in the active fragment of the molecular system, where N is the number of electrons in the active part. The molecular orbitals of this fragment are usually varied during the calculation, and hence its density as well. Contrarily, the density of the frozen part, r f (r), is chosen to be constraint during the calculation. In the context of the related subsystem DFT method, [32][33][34] it consists of orbitals located in the environment, f i f , and both subsystem densities can be optimized. Please note that we do not use an extra label for indices of orbitals, densities, and potentials of the active fragment, while the terms for the frozen fragment are marked with ''f'' as a subscript (or, in the case of potentials, as a superscript).
In order to obtain the ground-state density of the active system, r(r), in the presence of the frozen density, the Kohn-Sham (KS) equations with Constrained Electron Density (KSCED) 15 must be solved, where v eff [r](r) has the same form as the effective potential for the standard Kohn-Sham equation of the isolated active system, and v emb [ r,r f ](r) is responsible for the interaction of the active part with the frozen environment density, þ v nadd kin r; r f ½ ðrÞ: The last two terms represent the so-called non-additive XC potential and non-additive kinetic potential, which read The term in parenthesis on the left hand-side of eqn (4) is denoted as F and can be identified as the time-independent FDE-KS Hamiltonian.

Frozen-density embedding within TDDFT
The FDE formalism has been extended to excited states by Casida and Wesolowski, 20 and a general subsystem TDDFT framework was established in ref. 27. The FDE approach can be considered as a special case of the subsystem-based TDDFT scheme in which local excitations are obtained under the influence of a given environment, e.g. a solvent. The matrix elements of the timeindependent FDE-KS Hamiltonian are given as where f p are orbitals of the active part of the system. The convention that indices i,j, etc. correspond to occupied orbitals, a, b, etc. to virtual orbitals and p,q, etc. refer to general orbitals is adopted for the following derivations. Expanding the timedependent one-electron functions, f p (r,t), in a linear combination of time-independent one-electron functions, f j (r) with time-dependent orbital coefficients, K jp (t), the Time-Dependent KS (TDKS) equations will take the following form 35,36 i @ @t K ¼ FK; (11) where the i-th column of the matrix K contains the expansion coefficients of f i (r,t). F is a matrix representation of the TDKS Hamiltonian operator in the basis of {f i (r)}. Assuming that a time-dependent external field is applied 35 one can derive the main equation of linear-response TDDFT within the FDE approach from eqn (11), where the elements of A and B read (real orbitals are assumed), The response of the XC potential corresponds to the second functional derivative of the exchange-correlation energy, which is called XC kernel and given as Only non-additive terms from the embedding potential, v emb [ r,r f ](r), remain and they are represented as a second functional derivative of the corresponding non-additive energies, which are also called non-additive kernels, Using the fact that in the basis of canonical orbitals elements of matrices A and B change to B ia,jb = (ia|bj) + (ia| f XC |bj) + (ia| f nadd kin |bj) + (ia| f nadd XC |bj).
In the ADF package XC, non-additive XC, and non-additive kinetic kernels in the equations above are by default evaluated within the adiabatic local density approximation (ALDA) for non-hybrid XC potentials. This gives a good compromise between accuracy and efficiency. Moreover, ALDA kernels show smoother behavior for the low density regions than adiabatic generalized gradient approximation (AGGA) kernels. The last statement can be explained by the absence of 1/|rr(r)| b dependencies in kernels, where b is a positive exponent. We note that it is sometimes useful to rearrange eqn (12) by means of the unitary transformation U, into the two coupled equations

TDDFT excitation-energy gradients within FDE
Having the concepts of excitation-energy gradients for regular TDDFT at hand, 5 The quantity G[X,Y,o] has the following form, where A and B are the matrices specified in Section 2.2, X and Y are the solution vectors from eqn (12), C is the matrix of orbital coefficients. Z and W in eqn (25) contain two additional sets of Lagrangian multipliers. The sum P ia Z ia F ia takes care of the conditions F ia = 0 (which implies that KSCED equations are solved, but not necessarily in canonical orbitals) and P pq;p q W pq S pq À d pq À Á ensures orthonormality of the orbitals (S pq is the overlap integral for orbitals p and q). Note that our definition of the vector Z implies that only occupied-virtual components among all Z rs are non-zero. The FDE-TDDFT eigenvalue problem can be expressed as and the orthonormality constraint is given as It is possible to rewrite the Lagrangian as (see Appendix A.3; for simplicity of notation the Lagrangian, L[X,Y,o,C,Z,W], will be written from now on as L) In the equation above, a relaxed difference density matrix P was introduced (see Appendix A.6 for more details), where T is an unrelaxed difference density matrix with elements equal to T ia = T ai = 0.
We now introduce some external perturbation x. In our particular case, x represents some nuclear coordinate R A k (i.e., component k of the position vector of atom A). The derivatives of the Lagrangian with respect to all possible nuclear coordinates x in the system constitute the gradient of the excitation energy with respect to the position of the nuclei. As stressed by Furche and Ahlrichs, the formalism is also valid if x represents some general perturbation, e.g. a component of a static electric field, which gives access to the excited-state dipole moment. 5 In order to differentiate the Lagrangian with respect to a perturbation x, it is useful first to convert it to atomic orbital representation (see also Appendices A. 4  For convenience, we will sometimes drop the notation v emb [ r,r f ], v eff [ r], etc. and write v emb , v eff instead. Converting P pq F pq P pq to atomic orbital basis, we get For convenience, we define the kinetic energy operatorT s ¼ À 1 2 r 2 . Keeping the orbital coefficients constant and taking the derivative with respect to the perturbation x [x as a superscript indicates differentiation over x, while (x) means that MO coefficients are kept constant during differentiation], we get from the equation above, The first term reads 10 ð X mn w mT s w n À Á x P mn ¼ T s r PðxÞ h i : Here, we introduced the density The second term P mn v eff mn ðxÞ reads (here, we introduce a simplified notation in which the integration of a function Ð f ðrÞdr or Ð f ðr; r 0 Þdr dr 0 will be expressed as where v C and v P C are the Coulomb potentials due to r and r P , respectively, and is an element of the ground-state density matrix. Similarly, we get for the third term P mn v emb mn ðxÞ P mn , where v f C is the Coulomb potential due to r f . Converting the third term in eqn (29) to the atomic basis representation and taking its derivative with respect to x leads to X or, in terms of densities, where g XC , g nadd XC , and g nadd kin are given as functional derivatives of the kernels (third functional derivatives of the corresponding energy functionals with respect to the density).
The derivative of the last term in eqn (29) with respect to x gives zero. Summing up all the terms mentioned above together with À P mn S x mn W mn ! the expression for the Lagrangian deriva-

Z-vector equation for FDE-TDDFT
For the evaluation of L x , the corresponding vectors W and Z have to be known (P = Z + T, and T is already determined from the solution of FDE-TDDFT problem). These vectors can be obtained from the so-called Z-vector equation, 5,37 which is derived using the condition, If eqn (44) is multiplied by C mq and summed over index m, we get the following relation, where Q stands for Introducing the abbreviation for some arbitrary M, it is possible to derive the following equations (see ref. 5 and Appendix B), where the explicit expressions for Q pq look like (four cases are distinguished), Subtracting eqn (50) from eqn (49) we get the Z-vector equation This equation can be rewritten as X where R ia = À(Q ia À Q ai ) reads The matrix W is determined from the following equations The Z-vector equation and all relevant variables can be rewritten in terms of densities. 10 The corresponding quantities H ij + [M], Q ij , and Q ia then read, while Q ai and Q ab will remain the same. Note that r M and r (X+Y) are defined in analogy to eqn (37). Finally, one can rewrite R ia in terms of densities

Details of the implementation
An algorithm for the evaluation of excitation-energy gradients within FDEu-TDDFT was implemented in ADF. 11 Our implementation is based on the modules for the calculation of excitation-energy gradients in TDDFT, which were developed by Seth et al. 10 In addition, we created modules for the evaluation of the second functional derivative of the non-additive kinetic energy using the AGGA (denoted as FULL in the formulas below to distinguish from the commonly used ALDA approximation), and the second and third functional derivatives of the non-additive kinetic energy using the Thomas Fermi approximation (denoted as TF in the formulas below). These modules make use of the XCFUN library. 38 As was already mentioned in Section 2.2, the evaluation of the second functional derivative in response calculations is usually done within the ALDA in ADF, independently of the approximation used for the XC potential. This introduces some consistency issues in the evaluation of TDDFT gradients. Every kernel and potential term in the expression for the gradient which originates from the Fock matrix (XC potential and second XC functional derivative) have to be evaluated in the original approximation (LDA, GGA,. . .), while terms which result from the response part are always evaluated in the ALDA (which in the case of the non-additive kinetic energy is the TF approximation). Within the ADF setup, the ALDA-consistent expression for the excited-state gradients thus reads, The difference in the approximations for the XC part has to be taken into account in the Z-vector equation, eqn (57), as well. Special attention should be paid to the elements H ij when M is (X + Y), and when M is T or Z. An additional comment should be made on the evaluation of the non-additive XC kernels. One can significantly simplify the evaluation of the XC kernels and their functional derivatives if the expression for the non-additive XC potential, eqn (7) is recalled, This basically means that instead of the evaluation and summation of XC kernels for active and total densities, one needs to evaluate the XC kernels only for the total density r tot (r).
Concerning the computational scaling of the FDEu-TDDFT gradient calculations, we note that only a small additional step in the calculation of embedding potentials and kernels is needed, compared to the effort in the standard TDDFT gradients. In particular, the orbital space is not increased in FDEu-TDDFT calculations when using monomer basis sets for the active system (which is the default in these types of calculations).
The differences in timing and computational effort of the ground-state DFT step and in TDDFT excitation-energy calculations have been addressed in ref. 23, and we refer the interested reader to this reference for further details.

Computational details
A locally modified version of the ADF package 11,12 was used for all calculations. For ground-state structure optimizations, the Becke-Perdew (BP86) XC functional 39,40 in combination with the triple-z polarized Slater-type (TZP) basis set from the ADF basis set library was employed. The same basis set and the XC functional were used for the excitation energy calculations. For the FDE calculations, we employed the corresponding ADF implementation, 41 making use of the so-called GGA97 generalized-gradient approximation (GGA) for the non-additive kinetic contribution to the embedding potential. 42 This approximation is also known as PW91k, because it has the same functional form for the enhancement factor as used in the exchange functional of Perdew and Wang (PW91). 43

Results and discussions
To test the analytical gradients for excitation energies within the FDEu approach, three different test systems were studied. A first test was performed on the complex of formaldehyde with a lithium ion (charge +1), where the density of the lithium ion was frozen. The structure of the system was obtained in the following way: at first, the formaldehyde structure was optimized with BP86/TZP; then the lithium ion was put on the axis along the CO bond (oriented along the z-axis). The structure is shown in Fig. 1. Such an orientation of the lithium ion was chosen to enforce a rather large effect of the embedding potential in this initial test. This system is thus highly polarized, and the effect of the classical electrostatic part of the embedding potential will be quite pronounced. But in addition, the nonadditive kinetic-energy contributions can be important in such systems, because they are necessary to prevent over-polarization effects.
The excitation energy dominated by the HOMO -LUMO transition was chosen for the calculation of the gradient at different CQO bond distances. Note that the CÁ Á ÁLi + distance is kept constant in our calculations, so that increasing CQO distances correspond to decreasing OÁ Á ÁLi + distances. The HOMO of this Fig. 1 Structure of the CH 2 OÁ Á ÁLi + system and its orientation in coordinate space. system corresponds to a nonbonding (n) orbital, while the LUMO corresponds to a p* orbital. Isosurface plots of the orbitals involved in the transition are shown in Fig. 2.
To check whether the embedding potential and all of its components affect the electronic transitions, the excitation energies for the CH 2 OÁ Á ÁLi + complex and for isolated formaldehyde were calculated for different CQO bond distances. These results are shown in Fig. 3. One can see from this figure that the closer the oxygen atom is to the lithium ion, the more the excitation energy is affected by the embedding potential. Thus, the complex represents a good test system to check the implemented procedure.
The gradients of the excitation energies 10 for isolated formaldehyde for different CQO bond distances were calculated analytically with BP86/TZP. Then for the same distances gradients were calculated numerically for the complex, where the lithium ion was described with FDEu. In addition, the analytical gradients of the excitation energy were calculated for the same system with FDEu. To compare the obtained results, the dependence of the z-component of the gradient for the oxygen atom on the CQO bond distance is plotted for each calculation in Fig. 4.
One can see from the plot that the presence of the lithium ion is simulated very well by the embedding potential, and that it contributes significantly to the value of the excitation-energy gradient. For all displacements of the oxygen atom towards lithium, the values of the z-component of the gradient for oxygen are in perfect agreement.
In the second test system, the lithium ion was substituted by a helium atom, and the OÁ Á ÁHe distance was changed while keeping the CQO distance fixed at 1.7 Å. This is significantly larger than the ground-state equilibrium CQO distance of about 1.2 Å, and was selected in order to observe a significant effect of the embedding contribution (at the equilibrium distance, this effect is rather small). Contrary to the first, the second test system has a small contribution of the classical electrostatic part to the embedding potential. The same excitation as in the previous case was considered, see Fig. 2. The effect of the FDE potential on the excitation energy under study can clearly be seen in Fig. 5.
Again, as in the previous case, the dependence of the z-component of the gradient for the oxygen atom on the OÁ Á ÁHe distance is plotted in Fig. 6 for the CH 2 O molecule alone as well as for the CH 2 OÁ Á ÁHe system (analytical/numerical). We observe a very good agreement between numerical and analytical results. The small remaining discrepancies can be attributed to inaccuracies of the numerical gradient. This can be seen by considering the region from 3.1 Å to 3.3 Å (see inset in Fig. 6). In this region the influence of the He atom on the gradient

View Article Online
should vanish and the result should be very close to the one for the isolated molecule. While the analytical FDE gradient approaches the result for the isolated molecule with good accuracy, the numerical one underestimates the z-component of the gradient by roughly 0.0005 a.u. Å À1 . This shows the advantages of the analytical approach over the numerical one. Note that the error for the numerical approach might be larger in regions where the gradient changes faster, e.g. near 0.4 Å, see Fig. 6. The third test system consists of an acetone molecule and water treated with FDE. We take this system from our previous work on frozen density embedding. 44 The system was optimized with BP86/ TZP imposing C s symmetry, see Fig. 7. Here, we studied the excitation dominated by the HOMO -LUMO (np*) orbital transition. The orbital composition of the excitation considered for the gradient calculation is shown in Fig. 8.
In comparison to the previous tests where only single atoms were considered as frozen systems within FDE, here we deal with a frozen molecule. In this case, both electrostatic and nonadditive kinetic parts of the FDE potential are significant. Structures in which the CQO bond distance in acetone was varied from 1.2 Å to 2 Å were generated for the comparison of numerical and analytical gradients. Similar to our first test system, the carbon-water distance was kept fixed, i.e., longer CQO distances correspond to shorter OÁ Á ÁH 2 O contacts. In Fig. 9 the dependence of the excitation energy on the CQO     bond distance for isolated acetone and for the acetoneÁ Á Áwater system are plotted. One can see that the FDE potential has a slight effect on the excitation energy.
The analytical and numerical excitation energy gradients for C 3 H 6 OÁ Á ÁH 2 O as well as for isolated acetone were calculated for the above-mentioned structures. The z-component of the excitation energy gradient for the oxygen atom in acetone (analytical) and in acetoneÁ Á Áwater (analytical/numerical) are shown in Fig. 10. Again, the results from the analytical approach are in good agreement with the numerical ones.

Summary and conclusions
In this paper we have, for the first time, presented the formulation of analytical excitation-energy gradients in FDE-TDDFT. Furthermore, an algorithm and an implementation for the analytical evaluation of excitation-energy gradients within the FDEu approximation were developed and discussed in detail. We have shown how to incorporate the contributions from the embedding potential into the final expression for excitationenergy gradients consistently with the non-FDE implementation in ADF by Seth et al. 10 Special care has to be taken concerning approximations for non-additive kinetic and XC kernel and their functional derivatives in order to be consistent with the previous ALDA implementation. The kernels arising from the Fock operator should be evaluated in the original approximation, while terms originating from the response part are evaluated in the ALDA. Of course, the implementation is flexible and allows for an extension beyond the ALDA in future work.
The implemented procedure was successfully tested and yields very good agreement with numerical reference results. Three different test systems were considered, which were chosen according to the expected contribution of different parts of the FDE embedding potential. In particular, the correct behavior of the gradient was tested on formaldehyde with a lithium ion as the frozen part and on formaldehyde with a frozen helium atom. The former has a large contribution from classical electrostatics in the FDE potential, while in the latter system the contribution from the non-additive kinetic potential is more important. As another test with a frozen molecular environment system, an acetoneÁ Á Áwater complex was chosen. All calculations on the above-mentioned systems showed good agreement between analytical and numerical results. Moreover, it was demonstrated that results of analytical FDEu-TDDFT gradients are more accurate and thus preferable to the ones obtained by numerical differentiation.
The whole formalism was successfully implemented into the ADF package. In future work, we plan to address more complicated test systems to optimize the performance of the implementation and make FDE calculations beyond the ALDA possible. Also extensions beyond the approximation of a response localized on the active system are desirable if subsystem-based descriptions of excitonically coupled chromophores 45 or explicit solvent response effects 46 are of interest.
A Detailed derivation of the excitationenergy gradients

A.1 Introduction
In this appendix, we provide an in-depth derivation of the expression for excitation-energy gradients in FDE-TDDFT, based on the studies of Furche and Ahlrichs 5 as well as Seth et al. 10 The latter directly focuses on aspects relevant for the Slater-based ADF framework and thus provides the basis for the derivation of excited-state gradients in FDE-TDDFT, which is established here in Sections 2.3 and 2.4.

A.2 Excited state gradients for TDDFT
In FDE-TDDFT, the excited-state energy is represented as a sum of the ground-state energy and the excitation energy. If the gradient of the excited-state energy has to be evaluated, one has to consider the derivatives of these two terms. The derivatives of the former are formulated in ref. 31, and the derivative of the latter can be established in a similar way using Lagrangian techniques. We will only be concerned with the latter here. The variational formulation of FDE-TDDFT (by analogy to TDDFT case 5,47,48  The quantities L and D are super-matrices, The |X,Yi vector in Dirac notation is After differentiation and elementary manipulations, eqn (71) leads to a set of two matrix-vector equations, eqn (23) and (24) given in Section 2.2. The first term of the Lagrangian in eqn (73), hX,Y|L|X,Yi, can be rewritten in the following way, After substitution of this result into eqn (73) and writing the second term in an explicit way, the Lagrangian will take the following form for every combination ia, since one can use (X + Y) ia and (X À Y) ia instead of X ia and Y ia as free variables. These two equations lead to eqn (23) and (24). If the TDDFT eigenvalue problem is solved, then the Lagrangian is stationary with respect to (X + Y) ia and (X À Y) ia but not with respect to the molecular orbital coefficients C mp . This introduces difficulties if one studies the derivative of the Lagrangian with respect to some perturbation x, substituting x with each coordinate of every atom in the system. This is a rather inefficient way to calculate the gradient of the excitation energy. But in 1984, Handy and Schaefer 37 suggested an efficient method to avoid solving these 3M equations, which instead requires to solve only one equation with the same dimension as the CPKS equation. This approach is called the Z-vector method 49 and it is set up in such a way that it does not require derivatives of molecular orbital coefficients. 5,37 It introduces the so-called relaxed difference density matrices 5 P (discussed in detail in Section A.6) and removes any implicit dependence of the Lagrangian on the molecular orbital coefficients. This can be done by introducing a new Lagrangian for the calculation of excitation energies in FDE-TDDFT (in analogy to TDDFT 5 ), and requiring it to be stationary with respect to the all molecular orbital coefficients C mp (for simplicity of notation the Lagrangian, L[X,Y,o,C,Z,W], will be written from now on as L), The Lagrangian multipliers Z ia and W pq ensure implicit relaxation of the ground-state density and also enforce the relaxed density to be a solution of the KS equation. In particular, P ia Z ia F ia takes care of the conditions F ia = 0 (which implies that the KS equations are solved), and P pq;p q W pq S pq À d pq À Á ensures orthonormality of the orbitals, The relations below require that the Lagrangian satisfies the TDDFT eigenvalue problem @L @hX; Yj where eqn (90) is the usual TDDFT eigenvalue problem and eqn (91) is the orthonormality constraint for the eigenvectors |X,Yi.

A.3 Expression for the Lagrangian in terms of MO integrals
Substituting the results of eqn (13) and (14) For ease of notation the unrelaxed difference density matrix T is introduced, (94) where the last equality holds due to the properties of (X + Y) and (X À Y), In eqn (92) the sum of terms containing F ab over two indexes i, j collapses to a sum over one index i because of d ij , which is then absorbed into T ab , see eqn (93). In the same way, the sum of terms containing F ij over two indexes a, b collapses to a sum over one index a (because of d ab ), which is included in the definition of T ij , see eqn (94). Considering eqn (95), the terms with F ij and F ab can be unified as terms with F pq (since F ia = F ai = 0) and eqn (92) can be simplified to Using this form of G[X,Y,o], eqn (86) takes the form, After the introduction of the relaxed difference density matrix P = T + Z, the Lagrangian can be rewritten as Recalling the Fock matrix, eqn (9), and using the definition of the core hamiltonianĥ ¼ À 1 2 r 2 þ v ext ðrÞ, we get and we can finally transform the Lagrangian to the following form,

A.4 Expression for the Lagrangian in terms of AO integrals
The imposed requirement for the Lagrangian to be stationary with respect to the elements of C, cancels out the implicit dependence of L on a perturbation x through the MO coefficients C mp . Hence it is more convenient to convert all terms to atomic orbitals by means of the relation, Noting that and defining, the first term in eqn (102) results in, From the second term, we obtain, where is an element of the ground-state density matrix. The third term reads X Similarly, the first part of the fifth term can be re-written as X pq;p q Grouping everything together, the Lagrangian in terms of atomic integrals is established (except for the last two terms, which will not survive the differentiation as shown in the following section),

A.5 Derivative of the Lagrangian in terms of AO integrals and densities
At this point the dependency of the Lagrangian, L, on x has been clearly established, and it is apparent that only atomicbasis-function integrals are involved in the differentiation with respect to x, i.e., L x = L (x) by means of the conditions in eqn (103). Then, after taking the derivative with respect to x, the following result occurs [x as a superscript indicates differentiation over x, while (x) means that MO coefficients are kept constant during differentiation] In the context of a DFT framework, it is more convenient to write this expression in terms of densities, in particular, for the purpose of implementation in the ADF code. Using the defini- (114) can be rewritten in terms of densities. 10 The Coulomb potential based on this density takes the following form (for ease of notation the integration of a function Ð f ðrÞdr will be defined as Eqn (114) will be analyzed term by term. Defining the kinetic energy operator asT s ¼ À 1 2 r 2 the first term from eqn (114) will take the following form: The sum containing derivatives of the matrix elements of the XC potential read X In the same way one can transform the sum with the derivatives of the embedding potential matrix elements X mn v emb Here, we have explicitly assumed that the external potential from the frozen region does not change, i.e., no perturbations in the positions of frozen atoms are included (although this can easily be generalized). The derivative of the term containing the kernels in eqn (114) The term with the Coulomb integrals which originates from the Fock matrix and from the TDDFT part can be transformed as follows: Grouping together all the terms mentioned above, it is finally possible to write the derivative of the Lagrangian in terms of densities A.6 Derivation of the Z-vector equation For the calculation of the Lagrangian derivative, eqn (114) and (121), the elements of W and Z have to be known. The matrix P contains the difference between the excited-and ground-state density matrices where T is the unrelaxed difference density matrix (''unrelaxed'' part) and is already determined once the TDDFT eigenvalue problem has been solved, since it contains products of transition density matrices (excitation vectors, TDDFT solution factors). The matrix Z takes care of the relaxation of the ground-state orbitals and can be of the same order of magnitude as T. 47 The information in P is complementary to the information contained in the transition vector |X,Yi (the transition density matrix). 47 P is related to the difference of expectation values for excited states and the ground state, while the transition vector is related to mixed matrix elements between the ground and excited states. P can give insight into the re-distribution of the electronic charge due to the excitation process. 47 These Lagrangian multipliers can be obtained from the so-called Z-vector equation, 5,37 which can be derived from the following condition: If eqn (123) is multiplied by C mq and summed over index m, it leads to the following equation: where Q is a matrix with elements, Four possible cases for the orbitals p and q in eqn (124) have to be considered at this point. The first case corresponds to occupied orbitals for p and q, so that we can set p = k r q = l. The right hand-side of eqn (124) can be expanded as follows: Taking into account that W rs = W sr , the following equation holds for any k, X The last equality holds because of S tk = d tk . The second term of the left hand-side of eqn (124) can be transformed by analogy with the derivation mentioned above The first sum on the right hand-side of eqn (128) is zero due to the condition F la = 0. Since d ka = 0 for every a, the second term of eqn (128) also vanishes, and finally eqn (124) changes to For the second case, where only virtual orbitals are considered for p and q, p = c r q = d in eqn (124), the result for the terms containing W rs will be similar to the previous case, namely W cd (1 + d cd ). The second term of the left-hand side in eqn (124) transforms as follows in this case: Since d ci = 0 for every i, the first term on the right-hand side of eqn (130) vanishes and the second term is zero due to the condition F id = 0. Then eqn (124) can be simplified for the case p = c r q = d to The third case concerns the combination of virtual (c) and occupied (k) orbitals, so that p = k, q = c in eqn (124), Taking into account that W rs = W sr , the following equation holds for a given k X s;s!k W ks S cs þ X r;r k W rk S rc ¼ X t W tk S tc þ W kk S kc ¼ W kc : The last equality holds true since S tc = d tc and k a c. The second term on the left-hand side of eqn (124) can then be rewritten in the following way: Since d ka = 0 for every a, the second term vanishes. In terms of canonical orbitals, the equation above changes to, Eqn (124) will change to Q kc + e c Z kc + H kc + [Z] = W kc .
Finally, the fourth case concerns the combination of occupied and virtual orbitals by making the choice p = c, q = k. The righthand side of eqn (124) will give the same result as in the case above (W kc ), while the second term on the left hand side will take the form The last step holds after making use of canonical orbitals and realizing that d ic = 0 for every i (the first term vanishes). Finally eqn (124) takes the form Q ck + e k Z kc = W kc .
Thus, the following set of equations has to be solved in order to get the vectors Z and W, where R ia = À(Q ia À Q ai ). It is common to call eqn (143) the Z-vector equation in the context of TDDFT, and it is usually written in the following form Once this equation is solved with respect to Z, it is possible to determine W from the following equations: To solve the Z-vector equation, eqn (144), explicit expressions for the elements of the vector Q have to be established: These expressions are derived in Section B.

B Elements of the vector Q
Here we present a derivation of the expressions for elements of the vector Q. First, two auxiliary relations will be derived in Section B.1 and then the necessary expressions will be obtained in Section B.2.

B.1 Auxiliary relation for the derivation of elements Q pq
If differentiation of eqn (78) with respect to (X + Y) la [note that (X + Y) ia and (X À Y) ia are considered as independent variables] is followed by multiplication with (X + Y) ka and summation over a, then due to the condition in eqn (79), the following equation is valid: