Environmental effects with Frozen Density Embedding in Real-Time Time-Dependent Density Functional Theory using localized basis functions

Frozen Density Embedding (FDE) represents a versatile embedding scheme to describe the environmental effect on the electron dynamics in molecular systems. The extension of the general theory of FDE to the real-time time-dependent Kohn-Sham method has previously been presented and implemented in plane-waves and periodic boundary conditions (Pavanello et al. J. Chem. Phys. 142, 154116, 2015). In the current paper, we extend our recent formulation of real-time time-dependent Kohn-Sham method based on localized basis set functions and developed within the Psi4NumPy framework (De Santis et al. J. Chem. Theory Comput. 2020, 16, 2410) to the FDE scheme. The latter has been implemented in its"uncoupled"flavor (in which the time evolution is only carried out for the active subsystem, while the environment subsystems remain at their ground state), using and adapting the FDE implementation already available in the PyEmbed module of the scripting framework PyADF. The implementation was facilitated by the fact that both Psi4NumPy and PyADF, being native Python API, provided an ideal framework of development using the Python advantages in terms of code readability and reusability. We demonstrate that the inclusion of the FDE potential does not introduce any numerical instability in time propagation of the density matrix of the active subsystem and in the limit of weak external field, the numerical results for low-lying transition energies are consistent with those obtained using the reference FDE calculations based on the linear response TDDFT. The method is found to give stable numerical results also in the presence of strong external field inducing non-linear effects.

Psi4NumPy framework to the FDE scheme. The latter has been implemented in its "uncoupled" flavor (in which the time evolution is only carried out for the active subsystem, while the environment subsystems remain at their ground state), using and adapting the FDE implementation already available in the PyEmbed module of the scripting framework PyADF. The implementation was facilitated by the fact that both Psi4NumPy and PyADF, being native Python API, provided an ideal framework of development using the Python advantages in terms of code readability and reusability.
We employed this new implementation to investigate the stability of the time propagation procedure, which is based an efficient predictor/corrector second-order midpoint Magnus propagator employing an exact diagonalization, in combination with the FDE scheme. We demonstrate that the inclusion of the FDE potential does not introduce any numerical instability in time propagation of the density matrix of the active subsystem and in the limit of weak external field, the numerical results for low-lying transition energies are consistent with those obtained using the reference FDE calculations based on the linear response TDDFT. The method is found to give stable numerical results also in the presence of strong external field inducing non-linear effects. Preliminary results are reported for high harmonic generation (HHG) of a water molecule embedded in a small water cluster. The effect of the embedding potential is evident in the HHG spectrum reducing the number of the well resolved high harmonics at high energy with respect to the free water. This is consistent with a shift towards lower ionization energy passing from an isolated water molecule to a small water cluster. The computational burden for the propagation step increases approximately linearly with the size of the surrounding frozen environment. Furthermore, we have also shown that the updating frequency of the embedding potential may be significantly reduced, much less that one per time step, without jeopardising the accuracy of the transition energies.

Introduction
The last decade has seen a growing interest in the electron dynamics taking place in molecules subjected to an external electromagnetic field. Matter-radiation interaction is involved in many different phenomena ranging from weak-field processes, i.e., photo-excitation, absorption and scattering, light harvesting in dye sensitized solar cells 1,2 and photo-ionization, to strong-field processes encompassing high harmonic generation, 3,4 optical rectification, 5,6 multiphoton ionization 7 and above threshold ionization. 8  Among different approaches, because of its compromise between accuracy and efficiency, the real-time time-dependent density functional theory (rt-TDDFT) is becoming very popular.
The main obstacle to implementing the rt-TDDFT method involves the algorithmic design of a numerically stable and computationally efficient time evolution propagator. This typically requires the repeated evaluation of the effective Hamiltonian matrix representation (Kohn-Sham matrix), at each time step. Despite the difficulties to realizing a stable time propagator scheme, there are very appealing features in a real-time approach to TDDFT, such as the absence of explicit exchange-correlation kernel derivatives 18 or divergence problems appearing in response theory and one has the possibility to obtain all frequency excitations at the same cost. Furthermore, the method is suitable to treat complex non-linear phenomena and external fields with an explicit shape, which is a key ingredient for the quantum optimal control theory. 19 Several implementations have been presented, [20][21][22] after the pioneering work of Theilhaber 23 and Yabana and Bertsch. 24 Many of them rely on the real space grid methodology 24 with Siesta and Octopus as the most recent ones. 25,26 Alternative approaches employ plane waves such as in Qbox 27 or QUANTUM ESPRESSO 28,29 and analytic atom centered Gaussian basis implementations (i.e., Gaussian, 30,31 NWChem, 32 Q-Chem 33,34 ) also have gained popularity. The scheme has been also extended to include relativistic effects at the highest level. Repisky  was also extended to the relativistic four-component framework based on the BERTHA code, [40][41][42] (more specifically based on the recently developed PyBERTHA, 37,38,43,44 that is the Python API of BERTHA).
The applications of the rt-TDDFT approach encompass studies of linear 45 and non-linear optical response properties, 25,46 molecular conductance, 47 singlet-triplet transitions, 48 plasmonic resonances magnetic circular dichroism, 49 core excitation, photoinduced electric current, spin-magnetization dynamics 50 and Ehrenfest dynamics. 51,52 Moreover, many studies in the relativistic and quasi relativistic framework appeared, ranging from X-ray near-edge absorption, 53 to nonlinear optical properties , 54 to chiroptical spectroscopy. 55 Most part of initial applications of real-time methodology to chemical systems were largely focused on the electron dynamics and optical properties of the isolated target systems. However, it is widely recognized that these phenomena are extremely sensitive to the polarization induced by the environment, such that the simulation on an isolated molecule is usually not sufficient even for a qualitative description. A number of studies aiming at including the effect of a chemical environment within rt-TDDFT have appeared in the literature.
They are based on the coupling of rt-TDDFT with the QM/MM approach which includes the molecular environment explicitly and at a reduced cost using classical mechanical description 31,56 or in a polarizable continuous medium (PCM), where the solvent degrees of freedom are replaced by an effective classical dielectric. [57][58][59] One of the challenges, however, in the dynamical description of the environment is that the response of the solvent is not instantaneous, thus these approaches have been extended to include the non-equilibrium solvent response. 60-63 A recent extension considers also non-equilibrium cavity field polarization effects for molecules embedded in an homogeneous dielectric. 64 Going beyond a classical description for the environment, very recently, Koh et al. 65 have combined the rt-TDDFT method with block-orthogonalized Manby-Miller theory 66 to accelerate the rt-TDDFT simulations, the approach is also suitable for cheaply accounting the solvation effect on the molecular response. Another fully quantum mechanical approach to include environment effects in the molecular response property is based on the frozen-density embedding (FDE) scheme. 67-69 FDE is a DFT-in-DFT embedding method that allows to partition a larger Kohn-Sham system into a set of smaller, coupled Kohn-Sham subsystems.
Additional to the computational advantage, FDE provides physical insight into the properties of embedded systems and the coupling interactions between them. 70 For electronic ground states, the theory and methodology were introduced by Wesołowski and Warshel, 71 based on the approach originally proposed by Senatore and Subbaswamy,72 and later Cortona, 73 for solid-state calculations. It has been further generalized 74,75 and directed to the simultaneous optimization of the subsystem electronic densities. Within the linear-response formalism Casida and Wesołowski put forward a formal TDDFT generalization 76 of the FDE scheme. Neugebauer 77,78 then introduced coupled FDE, a subsystem TDDFT formulation which removed some of the approximations made in the initial TDDFT-FDE implementations. Recently, the approach has been further extended 79,80 to account for charge-transfer excitations, taking advantage of an exact FDE scheme. [81][82][83][84][85] A DFT subsystem formulation of the real-time methodology has been presented in a seminal work by Pavanello and coworkers 70 together with its formulation within the FDE framework. They showed that the extension of FDE to rt-TDDFT can be done straightforwardly by updating the embedding potential between the systems at every time step and evolving the Kohn-Sham subsystems in time simultaneously. Its actual implementation, based on the use of plane-waves and ultrasoft pseudopotentials, 29,70 showed that the updating of the embedding potentials during the time evolution of the electron density does not affect the numerical stability of the propagator. The approach may be approximated and devised in the so called "uncoupled" scheme where the density response to the external field is limited to one active subsystem while keeping the densities of the other subsystems frozen in time. Note that also in this uncoupled version the embedding potential is time-dependent and needs to be recomputed and updated during the time propagation. However, the propagation scheme is restricted to the active subsystem and the approach is promising to include environmental effects in the real-time simulation. Numerous applications within the context of the linear response TDDFT showed that an uncoupled FDE is sufficient for reproducing supermolecular results with good accuracy even in the presence of hydrogen bonds as long as there are no couplings in the excitations between the systems.
In this work we extend rt-TDDFT based on localized basis functions to the FDE scheme in its uncoupled version (uFDE-rt-TDDFT), taking advantage of modern software engineering and code reusability offered by the Python programming language. We devised an unified framework based on Python in which the high interoperability allowed the concerted and efficient use of the recent rt-TDDFT procedure, which some of us have implemented in the framework of the Psi4Numpy API, 37,38 and the PyADF API. 86 The rt-TDDFT procedure has served as the main interface where the PyADF methods, which gave direct access to the key quantities necessary to devise the FDE scheme, can be accessed within a unified framework.
Since in this work we introduced a new flavor of the rt-TDDFT Psi4Numpy-based program, to avoid confusions, from now on, we will refer to the aforementioned rt-TDDFT based on Psi4Numpy as Psi4-rt, while its extension to the FDE subsystem framework will be referred as Psi4-rt-PyEmbed.
In Section 2 we review the fundamentals of FDE and its extension to rt-TDDFT methodology. In Section 3 computational details are given with a specific focus on the interoperability of the various codes we merged and used: Psi4Numpy, 39 XCFun 87 and PyADF, 86 including the PyEmbed module recently developed by some of us. In Section 4 we report and comment the results of the calculations we performed on excitation transitions for different molecular systems, including: a water-ammonia complex, a water cluster and a more extend acetone-in-water cluster case. Finally, we give some preliminary results about the applicability and numerical stability of the method in presence of intense external field inducing strong non-linear effects as High Harmonic Generation (HHG) in the active system. Concluding remarks and perspectives are finally given in Section 5.

Theory
In this section we briefly review the theoretical foundations of the FDE scheme and its extension to the rt-TDDFT methodology. As mentioned above, a previous implementation was presented by Pavanello et al. 70 using plane waves and ultrasoft pseudopotentials. We refer the interested reader to this seminal work for a general theoretical background, and for additional details of the FDE-rt-TDDFT formal derivation.

Subsystem DFT and Frozen Density Embedding formulation
In the subsystem formulation of DFT the entire system is partitioned into N subsystems, and the total density ρ tot (r) is represented as the sum of electron densities of the various subsystems [i.e., ρ a (r) (a = 1, .., N )]. Focusing on a single subsystem, we can consider the total density as partitioned in only two contributions as ρ tot (r) = ρ I (r) + ρ II (r). (1) The total energy of the system can then be written as with the energy of each subsystem (E i [ρ i ], with i = I, II) given according to the usual definition in DFT as In the above expression, v i nuc (r) is the nuclear potential due to the set of atoms which defines the subsystem and E i nuc is the related nuclear repulsion energy. T s [ρ i ] is the kinetic energy of the auxiliary non-interacting system, which is, within the Kohn-Sham (KS) approach, commonly evaluated using the KS orbitals. The interaction energy is given by the expression: with v I nuc and v II nuc the nuclear potentials due to the set of atoms associated with the subsystem I and II, respectively. The repulsion energy for nuclei belonging to different subsystems is described by the E I,II nuc term. The non-additive contributions are defined as: with X = E xc , T s . These terms arise because both exchange-correlation and kinetic energy, in contrast to the Coulomb interaction, are not linear functionals of the density.
The electron density of a given fragment (ρ I or ρ II in this case) can be determined by minimizing the total energy functional (Eq.2) with respect to the density of the fragment while keeping the density of the other subsystem frozen. This procedure is the essence of the FDE scheme and leads to a set of Kohn-Sham-like equations (one for each subsystem) where the non-additive exchange-correlation and kinetic energy contributions are defined as the difference between the associated exchange-correlation and kinetic potentials defined using ρ tot (r) and ρ I (r). For both potentials, one needs to account for the fact that only the density is known for the total system so that potentials that require input in the form of KS orbitals are prohibited. For the exchange-correlation potential, one may make use of accurate density functional approximations and its quality is therefore similar to that of ordinary KS.
The potential for the non-additive kinetic term ( δT nadd δρ I (r) , in Eq.7) is more problematic as less accurate orbital-free kinetic energy density functionals (KEDFs) are available for this purpose. Examples of popular functional approximations applied in this context are the Thomas-Fermi (TF) kinetic energy functional 88 or the GGA functional PW91k. 89 These functionals have shown to be accurate in the case of weakly interacting systems including hydrogen bond systems, whereas their use for subsystems interacting with a larger covalent character is problematic (see Ref. 81 and references therein). The research for more accurate KEDFs is a key aspect for the applicability of the FDE scheme as a general scheme, including the partitioning of the system also breaking covalent bonds. 90 In general, the set of coupled equations that arise in the FDE scheme for the subsystems have to be solved iteratively. Typically, one may employ a procedure of "freeze-and-thaw" where the electron density of the active subsystem is determined keeping frozen the electron density of the others subsystems, which is then frozen when the electron density of the other subsystems is worked out. This procedure may be repeated many times until all subsystems' densities are converged. In this case the FDE scheme can be seen as an alternative formulation of the conventional KS-DFT approach for large systems (by construction it scales linearly with the number of subsystems). The update of the density for (part of) the environment can be important when trial densities obtained from isolated subsystems are not are not a very good starting point, as is the case for ionic species. [91][92][93] The implementation of FDE is relatively straightforward, in that the v I emb (r) potential is a one-electron operator that needs to be added to the usual KS hamiltonian. When using localized basis functions, the matrix representation of the embedding potential (V emb ) may be evaluated using numerical integration grids similar to those used for the exchangecorrelation term in the KS method. This contribution is then added to the KS matrix and the eigenvalue problem is solved in the usual self-consistent field manner.
We note that, irrespective of whether one or many subsystem densities are optimized, the matrix V emb needs to be updated during SCF procedure because it also depends on the density of the active subsystem (see Eq.7).
Going beyond the ground state is necessary to access many interesting properties, which for DFT are expressed via response theory, [76][77][78]94,95 such as electronic absorption 96 or NMR shielding, 93,97 and for which FDE has been shown to work properly since these are quite often relatively local. In a response formulation, the embedding potential as well as its derivatives enter the equations and, if more than one subsystem is allowed to react to the external perturbations, 77,78,94,95 the derivatives of the embedding potential introduce the coupling in the subsystems' response (as the embedding potential introduces the coupling of the subsystems' electronic structure in the ground state).
While such couplings in response may be very important in certain situations, such as for strongly interacting systems 79,80 or for extensive properties, 78 disregarding them still can provide a very accurate picture, notably for localized excited states. 91,96 In this simplified "uncoupled" framework, one considers only the response of the subsystem of interest (and thus the embedding potential and its derivative with respect to this subsystem's density).
While neglecting environment response may seem a drastic approximation, good performance relative to supermolecular reference data has been obtained for excitation energies of a chromophore in a solvent or a crystal environment, even when only retaining the embedding potential. 91 We will therefore employ this framework in the following.

The Real-Time Time-Dependent Kohn-Sham method and its extension to FDE
where i is the imaginary unit and D(t) and F (t) are the one-electron density matrix and time-dependent Kohn-Sham matrix, respectively. The above equation holds in both the non-relativistic and relativistic four-component formulations. 35,40 In a non-relativistic framework, the Kohn-Sham matrix (F (t)) is defined as where T and v nuc are the one-electron non-relativistic kinetic energy and intramolecular nuclear attraction terms, respectively. The explicit time dependence of F (t) is due to the time-dependent external potential v ext (t), which accounts for the interaction of the molecular system with an applied external electric field. Even in the absence of an external field the Fock operator is implicitly dependent on time through the density matrix D(t) in the The propagation in time of the density matrix can be expressed as where U (t, t 0 ) is the matrix representation of the time-evolution operator.
If we start (initial condition, i.e. initial time t 0 ) with the electronic ground state density matrix and use as orthonormal basis the ground state molecular orbitals, D(t 0 ) assumes the form where 1 oo is the identity matrix over the occupied orbital space of size n occ (total number of electrons). The D matrix has the dimension of n tot (n tot = n occ + n virt ) that is total number of the basis functions.
In our implementation, which uses a basis set of atomic centered (AO) Gaussian-type functions, the ground-state molecular orbitals are conveniently used as the reference orthonormal basis and at the time t the Fock and density matrices are related to their AO basis representation simply by: where the C matrix contains the reference MO expansion coefficients. The same coefficients satisfy a similar relation for D(t) M O : In a finite time interval, the solution of the Liouville-von Neumann equation consists in the calculation of the Fock matrix at discrete time steps, and in propagating the density matrix in time.
In the most general case, where the Fock operator depends on time even in absence of external fields, the time-evolution operator can be expressed by means of a Dyson-like series: which in compact notation, using the the time ordering operatorT , reads as: The time ordering is necessary since F (t) at different times do not necessarily commute . Typically, this time-ordering problem is overcome by exploiting the composition property of time-evolution operator (U (t, t 0 ) = U (t, t 1 )U (t 1 , t 0 )) and discretiz- The Magnus expansion has found the widest application, in particular, in those implementations that employ localized basis sets functions, for which matrix exponentiation can be performed exactly via matrix diagonalization. Typically, the Magnus expansion is truncated to the first order evaluating the integral over time using numerical quadrature, provided that the time interval ∆t is sufficiently short. Using the midpoint rule the propagator becomes This approach, also referred as second-order midpoint Magnus propagator, is unitary by construction, provided that F is hermitian. This scheme exhibits an error which is proportional to (∆t) 3 . The expression in Eq.15 coincides with the so-called modified-midpoint unitary transform time-propagation scheme originally introduced by Schlegel et al. 21 The F matrix at time t + ∆t/2, where no density is available, can be obtained using an iterative series of extrapolations and interpolations at each time. Note that, if this predictor/corrector procedure is converged in a self-consistent manner the second-order midpoint Magnus propagator preserves the time reversal symmetry, which is an exact property of the equation of motion in absence of magnetic field. The predictor/corrector scheme is a key ingredient in preserving the numerical stability of the propagation with a range of algorithms that can be applied in this context. 103 We have recently implemented a particularly stable predictor/corrector scheme, originally proposed by Repiskyet al., 35 in the interactive quantum chemistry programming environment Psi4NumPy. [37][38][39] The methodology that we have described above can be straightforwardly extended to the subsystem density functional theory framework and in particular to FDE (FDE-rt-TDDFT). 70 In the present work we consider one active subsystem and keep frozen the density of the environment along the time propagation (uncoupled scheme, to which we will refer as uFDE-rt-TDDFT). Thus, a LvN type equation is solved in the space of the active subsystem. The only modification to Eq.8 is in the definition of the effective hamil-tonian matrix representation which now refers to the active subsystem (F +v ext (t)) and to which the matrix representation of the embedding potential (V emb (t)) is added to take into account the effect of the environment. The propagation scheme itself remains unaltered.
As in the case for the ground state, in which the change of the active subsystem density requires that V emb is updated at each SCF iteration, the time propagation of the electron density will introduce a time dependence in V emb even though the environment densities are kept frozen at their ground state value (due to the use of the uncoupled scheme).
Thus, the V emb matrix needs to be updated during the propagation. In the present implementation we use atomic centered Gaussian function as basis set for the active subsystem and evaluate the V emb µν matrix elements numerically. 91 We will show that the numerical noise associated with the construction of the embedding potential introduced by this scheme does not affect the numerical stability of the density matrix propagation in the linear and non-linear regimes. In the following sections we will also demonstrate, for a specific application, that the updating frequency of the embedding potential may be significantly reduced (much less than one per time step used to solve the LvN equation) without jeopardising the accuracy.
As usual, the key quantity in a real time simulation is the time-dependent electric dipole moment µ(t). Each Cartesian component p (with p = x, y, z) is given by where P p is the matrix representation of the p-th component of the electric dipole moment operator (see also Eq. 16). Since, in our uFDE-rt-TDDFT implementation, the time dependency response of the external field is due only from the active system, in the above expression (Eq.16) all quantities refer to the active subsystem. The vector µ(t) defines the polarization response to all orders and is easily computed by the electronic density at any time, t. From this quantity one can then compute both linear and non-linear properties.
In the linear response regime, each component of electric dipole moment, µ p (ω), with an external field E q in the direction q (with q = x, y, z), is given in frequency space by The components depend on the polarizability tensor (α pq ) through the Fourier-transformation of the q-component of the applied field. The dipole strength function S(ω) is related to the imaginary part of the frequency dependent linear polarizability by In our implementation 37 the perturbation can be chosen to be either an impulsive kick or a continuous wave whose amplitude is modulated by an analytic envelope function. Different explicit functional forms are available. 37,38 In the case of an impulsive perturbation (E(t) = kδ(t)n, where n is a unit vector representing the orientation of the field) we adopt the δanalytic representation as proposed in Ref. 35 One of the best-known examples of non-linear optical phenomena is HHG in atoms and molecules. HHG occurs via photo-emission by the molecular system in a strong field and can be also computed from µ(t). 104 In this work we calculate the HHG power spectrum for a particular polarization direction as the Fourier transform of the laser-driven induced dipole moment, Other suitable approaches have been investigated in the literature, 104 but in all cases the key quantity is µ(t).

Computational Details and Implementation
In this section we outline the computational strategy we adopted to implement the uFDErt-TDDFT scheme. We devised a multi-scale approach where we take advantage of the  As an explicit example of the interoperativity achieved between different codes we report in Algorithm 1 some basic directives used to compute those key quantities necessary for our uFDE-rt-TDDFT. The electron density of an active system is obtained via Psi4Numpy while the electron density, the Coulomb potential and non-additive terms of the environment are managed using PyADF. These quantities can be easily mapped on a common numerical grid and used in PyEmbed to evaluate the relative non additive embedding potential. Thus, the geometry and basis set of the active system (in this specific case a H 2 O molecule) are parsed at Line 7 and the ground state wavefunction object is returned by the psi4.energy() method. The corresponding electron density matrix is then obtained as a NumPy array by the h2o_wfn object. The electron density is mapped into a real-space grid representation using a preset numerical grid, and used to populate a suitable object container (Line [14][15][16][17][18][19][20]. A ground state calculation of the environment molecule (that is a NH 3 molecule in this example) is carried out using PyADF run() method (Line 23). In this case we use the adfsinglepointjob method to execute the corresponding ADF calculation. 113 We mention here that PyADF, despite its name, is not specific to this program, but works with a number of different quantum chemistry codes. The density and Coulomb potential resulting from this calculation, that are represented on a common numerical grid, are obtained using get_density() and get_potential() methods (Line 25,27) respectively. The PyEmbed module has all the methods needed to manage the density of both the reference system and environment to finally compute the non-additive embedding potential. Indeed, the embed_eval object is instantiated (Line 34) and the non-additive embedding potential is evaluated on the numerical grid using get_nad_pot (Line 36), once the density of both the active system and of the environment has been provided.

Rapid prototyping and implementation
Algorithm 1 has well illustrated how we can utilize the classes provided by PyADF to obtain a very simple workflow in which we are able to manipulate quantities coming from Psi4Numpy. Thus, we are now in a position to draw the main lines of our uFDE-rt-TDDFT implementation, the Psi4-RT-PyEmbed code. 108 In Figure 1 we present its pictorial workflow.  In the out-of-loop section the density and electrostatic potential of the environment are obtained as grid functions. The active system density matrix is expressed as grid function object and used to calculate the embedding potential. The active system density is optimized self-consistently according to Eq 7. The red star and the arrow pointing at it, symbolize that the out-of-loop blocks of tasks are involved only in the initial stage of the procedure. a) The relaxed active density matrix is exported as grid function. b) PyEmbed classes are used to calculate the embedding potential. c) The embedding potential is expressed on the finite basis set representation (GTO's). d) The active density matrix is evolved according to the real-time propagation scheme.
Algorithm 1 Illustrative Python code to compute active system density (using the Psi4Numpy code), environment density and Coulomb potential (using the ADF code) and non-additive embedding potential via the PyEmbed module. 1: import psi4 2: import pyadf 3: import pyadf.PyEmbed 4: from pyadf.Plot.GridFunctions import GridFunctionFactory 5: from pyadf.Plot.GridFunctions import GridFunctionContainer 6: ... 7: geom,mol = fde_util.set_input('h2o.xyz',basis_set) 8: # psi4 run 9: ene, h2o_wfn = psi4.energy(func,return_wfn=True) 10: # get psi4 h2o density 11: D = np.array(h2o_wfn.Da()) 12: ... 13  We start describing the out-of-loop section. Firstly the geometry and basis set of the environment are initialized (the orange left most block), thus the ADF package provides, through a standalone single point calculation, the electrostatic and nuclear potential of the environment and its density ρ II and a suitable integration grid for later use. At this stage all the basis sets and exchange-correlation functionals available in the ADF library can be used. In the next step, green block, the geometry and basis set of the active system are parsed from input and the ground state density ρ I is calculated using the Psi4Numpy related methods. The right pointing arrow, connecting the last block, sketches the mapping of the density matrix onto the real-space grid representation. The evaluation of ρ I (r) on the numerical grid is efficiently accomplished using the molecular orbitals (MO), which requires the valuation the localized basis functions at the grid points.
Finally, the PyEmbed module comes into play (last block of the out-of-loop section), the real-space electron densities ρ I and ρ II serve as input for the get_nad_pot() method. Thus the non-additive kinetic and exchange potential are obtained. The embedding potential is then calculated from its constituents (i.e. the environment electrostatic and nuclear potential and the non-additive contribution as detailed in Eq. 7) and evaluated at each grid point, v emb (r k ). The embedding potential matrix representation in the active subsystem basis set, V emb , is calculated numerically on the grid as where χ µ (r k ) are the Gaussian-type basis set functions employed in the active systems (used in Psi4Numpy) evaluated at the grid point, r k . In the above expression, w k are specific integration weights.
In the case of a FDE-rt-TDDFT calculation, the electron density of the active system at the beginning of the propagation (t 0 = 0, initial condition) is not the ground state density of the isolated molecule, rather a polarized ground state density. The latter is obtained through a self-consistent-field calculation in the presence of the embedding potential. We adopt the so-called split-scf scheme as described in Ref. 114. It should be noted that the density matrix, corresponding to the optimized ρ I electron density, is the input data for the green block (block a) of the in-loop section. The outgoing red arrow, connecting the out-and the in-loop branch of the diagram, it means that the former is only involved in the early step of the procedure and it will no longer come into play during the time propagation. As mentioned, the optimized density matrix of the active system as resulting from the SCF procedure including the For the sake of completeness, the pseudo code needed to evolve the density using the second-order midpoint Magnus propagator is reported in SI and relies on the methodology illustrated in Section 2.2. We refer the interested readers to our recent work on real-time propagation for further details. 37

Results and Discussion
In the present section we report a series of results mainly devoted to assess the correctness of the uFDE-rt-TDDFT scheme. To the best of our knowledge, this implementation is the first available for localized basis sets. Since our implementation relies on the embedding strategies implemented in PyADF, it appears natural and appropriate to choose as a useful reference the uncoupled FDE-TDDFT scheme, based on the linear response 75,115 and implemented in the ADF program package. 116  As shown in Table 1, convergence can be observed with both Psi4-rt and ADF-LR, in particular for the first low-lying transitions (additional excitation energies are reported in the Supporting Information). For some of the higher energy transitions the convergence is less prominent, pointing to deficiencies in the smaller basis sets. We mention that the results obtained using our Psi4-rt implementation perfectly agree with those obtained using the TDDFT implementation based on linear response implemented in the NWChem code, which uses the same Gaussian type basis set (see Table S1 in SI). Thus, we conclude that most of the deviations from the ADF-LR values can be ascribed to unavoidable basis set differences. A qualitatively similar pattern of differences is to be expected when including the environment effect within the FDE framework. Table 1: Excitation energies (in eV) corresponding to the first five low-lying transitions of the isolated water molecule. Data obtained using TDDFT based on linear response implemented in ADF (ADF-LR) and the our real-time TDDFT implemented (Psi4-rt). The labels (D, T, Q) correspond to data obtained using the Gaussian-type basis sets aug-cc-pVXZ (X = D,T,Q) and Slater-type basis sets AUG-X (X = DZP,TZ2P,QZ4P) which are used in the Psi4-rt and ADF-LR codes, respectively (see text for details). To assess differences in the presence of an environment, we next tested our uFDErt-TDDFT results against ADF-LR-FDE ones. The target system is the water-ammonia adduct, in which the water molecule is the active system that is bound to an ammonia molecule, which plays the role of the embedding environment. In the Psi4-rt-PyEmbed case we employed a contracted Gaussian aug-cc-pVXZ (X=D,T) basis set 118,119 for the active system whereas the basis set used in PyADF for the calculation of the environment frozen density (ammonia) and the embedding potential is the AUG-X (X =DZP,TZ2P) Slater-type set from the ADF library. 116 The ADF-LR-FDE employs the AUG-X (X =DZP,TZ2P) basis sets from the same library. For the real-time propagation of the active system (water), in both the isolated and the embedded case the BLYP 120,121 exchange-correlation functional is used, while the Thomas-Fermi and LDA functionals 122,123 have been employed for the non-additive kinetic and non-additive exchange-correlation potentential, respectively. The numerical results are reported in Table 2. Although, as expected, there is no quantitative agreement on the absolute value of the transitions, the shift ∆ (E iso. − E emb ) shows an acceptable agreement for the lowest transitions (see for additional excitation energies the Supporting Information). From these results, we conclude that our implementation is both stable and numerically correct, with differences between the methods explainable by the intrinsic basis set differences.

The water in water test case
To provide a further test of our implementation, we also computed the absorption spectra of a water molecule embedded in a water cluster of increasing size. The geometries of the different water clusters are taken from Refs. 124    In Table 3, and in Fig.2, we report how the time for the embedding potential calculation is distributed over the different tasks, when the number of surrounding water molecules increases from one to five. It is interesting to note that the time needed to evaluate the embedding potential increases almost linearly, for the limited number of water molecules considered here. The standard real-time iteration time (corresponding to the isolated water molecule) takes less than 1 sec and shows up as a fixed cost in the increasing computation time, while the time spent in the embedding part is dominated by the evaluation of the matrix representation for the active subsystem, e.g step c) of Fig.1 (see for instance t c column of Table 3). The time spent in this evaluation depends on the number of numerical integration points used to represent the potential, and can be reduced by using special grids for embedding purposes once the environment is large enough.

The acetone in water test case
As a further test of the numerical stability of accuracy of the method, we investigated the n → π * transition in the acetone molecule, both isolation and using an explicit water cluster to model solvation. In order to assess the shift due to the embedding potential, we calculate the absorption spectrum of the isolated molecule at the same geometry it has in the cluster model. The geometry for the solvated acetone system was taken from Ref. 91, corresponding to one snapshot from a MD simulation, where the acetone is surrounded by an environment consisting of 56 water molecules. The uFDE-rt-TDDFT calculation has been obtained specifying in our Psi4-rt-PyEmbed framework all the computational details. In particular, the frozen density of the environment is obtained from a ground state calculation using ADF in combination with the PBE functional and DZP basis set, while for the acetone we employ the BLYP functional and the Gaussian def2-svp basis set using the Psi4-rt code. The nonadditive kinetic and exchange-correlation terms of the embedding potential are calculated using the Thomas-Fermi and LDA functionals respectively. For the isolated acetone the n → π * transition is found at 3.73 eV whereas for the embedded molecule is located at 3.96 eV. The full absorption spectrum is reported in Fig. 4.
It is worth noting that, due to its low intensity, this transition is particularly challenging for a real-time propagation framework. To obtain a spectrum up to 11 eV, we carried out a simulation consisting of 20000 time steps and lasting 2000 a.u (48 fs). This relatively long simulation time demonstrates the numerical stability of the approach and its implementation.
As an overall check of our implementation we compare the shift of the n → π * transition  observed between isolated and embedded in a water cluster acetone obtained using both our Psi4-rt-PyEmbed and the ADF-LR-FDE methods. The active system response was calculated at BLYP level of theory, while Thomas-Fermi and LDA functionals were employed for the non-additive kinetic and exchange-correlation terms respectively of the embedding potential in the ADF-LR-FDE calculation. As one can observe by looking at the values reported in Table 4, we obtain a good agreement in the absolute values, both isolated and embedded acetone, and the computed shift is likewise in rather good agreement.

FDE-rt-TDDFT in the non-linear regime
A specificity in the the real-time approach is that the evolution of the electron density can be driven by an real-valued electric field whose shape can be explicitly modulated. Realistic laser fields can be modeled by a sine function of ω 0 frequency using any physically meaningful enveloping function. Using an explicit external field is a key tool in optical control theory, furthermore it is possible, employing high intensity field, to study phenomena beyond linear-response, i.e hyperpolarizability coefficients and high harmonic generation in molecules. The latter point will be detailed in the following section.
In this section we demonstrate that the uFDE-rt-TDDFT scheme gives stable numerical results not only in the perturbative regime, as shown above, but also in the presence of intense fields. Physically meaningful laser fields are adequately represented by sinusoidal pulse of the form E(t) = f (t) sin(ω 0 t) where ω 0 is the carrier frequency. In this work we employ, a cos 2 shape for the envelope function: 126 simulation) as In Fig. 6 we report the base-10 logarithm of the spectral intensity for the embedded water molecule and we compare it to the HHG of the isolated water calculated that has the same geometry it has in the cluster model. In the case of the isolated water molecule we are able to observe relatively well defined peaks up to the 21th harmonics. We mention that this finding qualitatively agrees with data obtained by Sun et al. 20 (see Figure 3 of Ref. 20 ).
An important parameter in the analysis of the HHG spectrum is the value of the energy cutoff (E cutof f ), which is related with the maximum number of high harmonics (N max ≈ E cutof f /ω 0 ). In a semiclassical formulation, 127  For the water molecule embedded in the cluster the same boundary can be approximately found corresponding to the 16th harmonic. The peaks at higher energies have a very small intensity and are much less resolved above the 16th harmonic. The flattening of the HHG intensity pattern is therefore solely due the introduction of the embedding potential of the surrounding cluster. The latter is consistent with a shift towards lower ionization energy passing from free water molecule to a small water cluster observed experimentally. 128

Computational constraints
Before concluding this work it may be interesting to put forward some assessments in terms of time statistics, to be used as a basis for optimizing the computation time and speed-up any uFDE-rt-TDDFT calculations. We are using a water-ammonia complex as a general test-case, where the geometry of the adduct has been taken from Ref. 129 and the water is the active subsystem.
In the real-time framework the embedding potential is, evidently, an implicit timedependent quantity. Since in the uncoupled FDE framework the density of the environment is kept frozen, the embedding potential depends on time only through the relatively small contributions given by the exchange-correlation and kinetic non-additive terms, which in turn it depends on time only through the density of the active subsystem. The electrostatic potential, due to the frozen electron density and nuclear charges of the environment, is the leading term in the overall potential. Thus, it may be reasonable to choose a longer time step for the update of the embedding potential, which is weakly varying in time.  4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 Table 5. Of course, as the number of time steps between consec-utive updates is increased (i.e. the embedding potential is updated less often), the total time needed to perform the full simulation goes down, as the time spent in computing the embedding potential decreases. The update rate of the embedding potential during the propagation affects to some extent the position of the peaks in the absorption spectrum. As can be seen in Fig. 7 the different traces corresponding to dipole strength functions calculated with different update rates, do not differ significantly and tend to coalesce as the number of time steps between consecutive updates decreases below 30 time steps. In particular, in the case of the lowest-energy transition, the energy shift corresponding to a quite long update period (roughly 300 time steps) is of the order of 0.02 eV. We also reported the partition between different tasks of the time needed for the calculation of the embedding potential in Table 6. As seen before, the calculation of the embedding potential is largely dominated by the projection to the basis set of the embedding potential from the numerical-grid representation. Therefore, some preliminary tests in reducing the number of grid points were carried out, and the results are presented in Fig. 8. It can be seen that there is no significant modification in the peak positions due to the use of a coarser integration grid: the overall spectrum is essentially stable and no artifacts are introduced. We furthermore note the possibility to use small grid localized solely on the active system by utilizing the fact that the embedding potential is projected on the localized basis set functions of the active system (see Eq.20), which makes it possible to neglect points on which these functions have a small value.

Conclusions and perspectives
In this work we have focused on the implementation of the Frozen Density Embedding scheme in the real-time TDDFT. We have integrated the Psi4Numpy real-time module we recently developed within the PyADF framework. We have devised a real-time FDE scheme in which the active density is evolved under the presence of the embedding potential. This implementation relies on a multiscale approach, since the embedding potential is calculated by means of PyADF, while the propagation is carried out by Psi4Numpy. We tested the implementation on a simple water cluster showing that the time needed for the propagation scales linearly with the cluster size. We studied many low-lying transitions in the case of a water molecule embedded in ammonia, and we showed that the shift of excitation energies with respect to the isolated water molecule is in good agreement with the results obtained using linear response FDE TDDFT implemented in ADF. Finally, we tackled a challenging case for rt-TDDFT, by computing the lowest-energy transition of acetone, which features an extremely low intensity. The corresponding signal can be identified in the computed spectrum, and we evaluated the solvatochromic shift due to the presence of a surrounding water cluster. We obtained a frequency shift of 0.225 eV, close to the reference value, 0.182 eV, from LR-FDE TDDFT as implemented in ADF. The scheme we developed has proven to be reliable also in the case of propagation in the non-linear regime. As a demonstration, we perturbed with a strong electric field a water molecule surrounded by five water molecules acting as frozen environment. Numerically stable induced dipole moment and corresponding emission spectrum were obtained.
Finally, we like to state that the present work provides an excellent framework for future developments. It is for instance possible and desirable to optimize the embedding potential construction. In our implementation (i.e. Psi4-rt-PyEmbed) the projection onto the basis set of the embedding potential from the numerical grid representation dominates the computational burden. The change of the embedding potential matrix in time, (i.e the difference at two consecutive time steps) depends on the relatively small contributions given by the exchange-correlation and kinetic non-additive terms. Significant improvement could be achieved by exploiting the sparsity of the matrix corresponding to that difference. Moreover, the use of smaller integration grid would probably further improve the procedure. Last but not least, the effect of relaxation of the environment has to be investigated. In our uncoupled FDE-rt-TDDFT scheme we are able to study local transitions within a given subsystem, and particularly those of the active system under the influence of the embedding potential due to the frozen environment. Thus we neglect transitions involving the environment and those due to the couplings of the subsystems. Relaxing the environment can be crucial both in the linear-response framework, in order to recover supramolecular excitations, and in the non-linear regime where a polarizable environment could heavily affect the hyperpolarizabilities of the target system. The limit of the uncoupled FDE scheme can be overcome by carrying out a simultaneous propagation of subsystems 70 and the computational framework developed in the present work represents an important step in that direction.