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 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, M.; J. Chem. Phys. 2015, 142, 154116]. In the current paper, we extend our recent formulation of the real-time time-dependent Kohn-Sham method based on localized basis set functions and developed within the 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 on 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 the 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 a strong external field inducing nonlinear 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 toward 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 than one per time step, without jeopardizing 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., photoexcitation, absorption, scattering, light harvesting in dye-sensitized solar cells,[1,2] and photoionization to strong-field processes encompassing high harmonic generation,[3,4] optical rectification,[5,6] multiphoton ionization,[7] and the above threshold ionization.[8] Furthermore, the emergence of new free electron laser (FEL) and attosecond methodologies[9,10] opened an area of research in which experiments can probe electron dynamics and chemical reactions in real time and the movement of electrons in molecules may be controlled. These experiments can provide direct insights into bond breaking[11−13]/forming[14] and ionization[15,16] by directly probing nuclear and electron dynamics.

Real-time time-dependent electronic-structure theory, in which the equation of motion is directly solved in the time domain, is clearly the most promising for investigating time-dependent molecular response and electronic dynamics. The recent progress in the development of these methodologies is impressive (see, for instance, a recent review by Li et al.[17]). 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 in 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 the response theory, and one has the possibility of obtaining all frequency excitations at the same cost. Furthermore, the method is suitable to treat complex nonlinear phenomena and external fields with an explicit shape, which is a key ingredient in the quantum optimal control theory.[19]

Several implementations have been presented[20−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 being 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 et al. proposed the first application and implementation of relativistic TDDFT to atomic and molecular systems[35] based on the four-component Dirac Hamiltonian, and almost simultaneously Goings et al.[36] published the development of X2C Hamiltonian-based electron dynamics and its application to the evaluation of UV–vis spectra. Very recently, some of us presented an rt-TDDFT implementation[37,38] based on state-of-the-art software engineering approaches (i.e., including interlanguage communication between high-level languages such as Python, C, FORTRAN, and prototyping techniques). The method, based on the design of an efficient propagation scheme within the Psi4NumPy[39] framework, was also extended to the relativistic four-component framework based on the BERTHA code[40−42] (more specifically based on the recently developed PyBERTHA,[37,38,43,44] i.e., the Python API of BERTHA).

The applications of the rt-TDDFT approach encompass studies of linear[45] and nonlinear 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, there have been many studies in the relativistic and quasi-relativistic framework, ranging from X-ray near-edge absorption[53] to nonlinear optical properties[54] to chiroptical spectroscopy.[55]

Mostly, initial applications of the real-time methodology to chemical systems were largely focused on the electron dynamics and optical properties of 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−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 nonequilibrium solvent response.[60−63] A recent extension also considers nonequilibrium cavity field polarization effects for molecules embedded in a homogeneous dielectric.[64]

Going beyond a classical description for the environment, very recently, Koh et al.[65] combined the rt-TDDFT method with the block-orthogonalized Manby–Miller theory[66] to accelerate the rt-TDDFT simulations; the approach is also suitable for cheaply accounting for the solvation effect on the molecular response. Another fully quantum mechanical approach to include environmental 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 one to partition a larger Kohn–Sham system into a set of smaller, coupled Kohn–Sham subsystems. In addition to the computational advantage, FDE provides physical insights 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 that 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−85]

A DFT subsystem formulation of the real-time methodology has been presented in a seminal work by Pavanello and co-workers[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 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 a 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, henceforth in this article, 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 to as Psi4-rt-PyEmbed.

In Section , we review the fundamentals of FDE and its extension to the rt-TDDFT methodology. In Section , 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 , we report and comment on 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 the presence of an intense external field inducing strong nonlinear effects such as high harmonic generation (HHG) in the active system. Concluding remarks and perspectives are finally given in Section .

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() is represented as the sum of electron densities of the various subsystems [i.e., ρ() (a = 1, ..., N)]. Focusing on a single subsystem, we can consider the total density as partitioned in only two contributions asThe total energy of the system can then be written aswith the energy of each subsystem (E[ρ], with i = I, II) given according to the usual definition in DFT asIn the above expression, vnuc() is the nuclear potential due to the set of atoms that defines the subsystem and Enuc is the related nuclear repulsion energy. Ts[ρ] is the kinetic energy of the auxiliary noninteracting system, which is, within the Kohn–Sham (KS) approach, commonly evaluated using the KS orbitals. The interaction energy is given by the expressionwith vnucI and vnucII being the nuclear potentials due to the set of atoms associated with subsystems I and II, respectively. The repulsion energy for nuclei belonging to different subsystems is described by the EnucI,II term. The nonadditive contributions are defined aswith X = Exc, Ts. 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 ) 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)which are coupled by the embedding potential term vembI(), which carries all dependence on the other fragment’s density. In this equation, veffKS[ρI]() is the KS potential calculated on the basis of the density of subsystem I only, whereas the embedding potential takes into account the effect of the other subsystem (which we consider here as the complete environment). In the framework of the FDE theory, vembI() is explicitly given bywhere the nonadditive exchange correlation and kinetic energy contributions are defined as the difference between the associated exchange-correlation and kinetic potentials defined using ρtot() and ρI(). 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 nonadditive kinetic term (, in eq ) 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 been shown to be accurate in the case of weakly interacting systems including hydrogen-bond systems, whereas their use in 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 arises in the FDE scheme for the subsystems has 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 other 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 a very good starting point, as is the case for ionic species.[91−93]

The implementation of FDE is relatively straightforward, in that the vembI() 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 (Vemb) may be evaluated using numerical integration grids similar to those used for the exchange-correlation 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 Vemb needs to be updated during the SCF procedure because it also depends on the density of the active subsystem (see eq ).

Going beyond the ground state is necessary to access many interesting properties, which for DFT are expressed via the response theory,[76−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 enters 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 can still 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 like 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.

Real-Time Time-Dependent Kohn–Sham Method and Its Extension to FDE

The time-dependent equation for the Kohn–Sham method can be conveniently formulated in terms of the Liouville–von Neumann (LvN) equation. In an orthonormal basis set, the LvN equation readswhere i is the imaginary unit and (t) and (t) are the one-electron density matrix and time-dependent Kohn–Sham matrix, respectively. The above equation holds in both the nonrelativistic and relativistic four-component formulations.[35,40]

In a nonrelativistic framework, the Kohn–Sham matrix ((t)) is defined aswhere and nuc are the one-electron nonrelativistic kinetic energy and intramolecular nuclear attraction terms, respectively. The explicit time dependence of (t) is due to the time-dependent external potential 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 (t) in the Coulomb ([ρ(t)]) and exchange-correlation terms (xc[ρ(t)]).

The propagation in time of the density matrix can be expressed aswhere (t,t0) is the matrix representation of the time-evolution operator.

If we start (initial condition, i.e., initial time t0) with the electronic ground-state density matrix and use as orthonormal basis the ground-state molecular orbitals, (t0) assumes the formwhere 1oo is the identity matrix over the occupied orbital space of size nocc (total number of electrons). The matrix has the dimension of ntot (ntot = nocc + nvirt), which is the total number of 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 time t the Fock and density matrices are related to their AO basis representation simply bywhere the matrix contains the reference MO expansion coefficients. The same coefficients satisfy a similar relation for (t)MOIn a finite time interval, the solution of the Liouville–von Neumann equation consists of 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 the absence of external fields, the time-evolution operator can be expressed by means of a Dyson-like serieswhich in compact notation, using the time-ordering operator , reads asThe time ordering is necessary since (t) at different times does not necessarily commute ([(t),(t′)] ≠ 0). Typically, this time-ordering problem is overcome by exploiting the composition property of the time-evolution operator ((t,t0) = (t,t1) (t1,t0)) and discretizing the time using a small time step. It is clear that the exact time ordering can be achieved only in the limit of an infinitesimal time step. Many different propagation schemes have been proposed[98] in the context of rt-TDDFT. Among others, we mention the Crank–Nicholson,[99] Runge–Kutta,[100] or Magnus[101,102] methods.

The Magnus expansion has found the widest application, in particular, in those implementations that employ localized basis set 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 becomesThis approach, also referred to as the second-order midpoint Magnus propagator, is unitary by construction, provided that is Hermitian. This scheme exhibits an error that is proportional to (Δt)3. The expression in eq coincides with the so-called modified-midpoint unitary transform time-propagation scheme originally introduced by Schlegel et al.[21]

The 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 the absence of a magnetic field. The predictor/corrector scheme is a key factor 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 Repisky et al.,[35] in the interactive quantum chemistry programming environment Psi4NumPy.[37−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, which we will refer to as uFDE-rt-TDDFT). Thus, an LvN-type equation is solved in the space of the active subsystem. The only modification to eq is in the definition of the effective Hamiltonian matrix representation, which now refers to the active subsystem (I(t) = I + nucI + xc[ρI(t)] + [ρI(t)] + ext(t)) and to which the matrix representation of the embedding potential (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 emb is updated at each SCF iteration, the time propagation of the electron density will introduce a time dependence in emb even though the environment densities are kept frozen at their ground-state value (due to the use of the uncoupled scheme).

Thus, the emb matrix needs to be updated during the propagation. In the present implementation, we use atomic centered Gaussian function as the 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 nonlinear 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 jeopardizing 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 bywhere is the matrix representation of the p-th component of the electric dipole moment operator (see also eq ). 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 ), 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 nonlinear properties.

In the linear-response regime, each component of the electric dipole moment, μ(ω), with an external field E in the direction q (with q = x, y, z), is given in the frequency space byThe components depend on the polarizability tensor (α) 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 byIn 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 ((t) = kδ(t) , where 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 nonlinear optical phenomena is HHG in atoms and molecules. HHG occurs via photoemission by the molecular system in a strong field and can also be 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 momentOther 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 uFDE-rt-TDDFT scheme. We devised a multiscale approach where we take advantage of the real-time TDDFT reference procedure, recently implemented within the Psi4Numpy framework (i.e., the Psi4-RT program),[37,38] while the FDE computational core relies on PyADF[86,105] and makes use of its PyEmbed module, which some of us have recently developed.[106,107] PyEmbed provides a Python implementation for computing the interaction energy (eq ) and embedding potential (eq ) from FDE on user-defined integration grids while using the XCFun library[18,87] to evaluate nonadditive exchange-correlation and kinetic energy contributions. With PyEmbed, quantum chemistry codes require only minimal changes: functionality to provide electron densities and its derivatives, as well as the electrostatic potential, over the grid, as well as to read in the embedding potential, and add it as a one-electron operator in the Fock matrix.[91] The PyADF scripting framework provides all of the necessary tools to manage various computational tasks and manipulate the relevant quantities for electronic-structure methods. The resulting Python code, referred to as Psi4-rt-PyEmbed, is available under GPLv3 license in ref (108). A data-set collection of computational results, including numerical data, parameters, and job input instructions used to obtain the absorption spectra of Sections , 4.3, and 4.5, is available in the Zenodo repository and can be freely accessed in ref (109).

Rapid Prototyping and Implementation

Psi4Numpy[39,110] and PyADF[86,105] both provide a Python interface, which greatly simplifies the computational workflow from input data to the results. PyADF is a quantum chemistry scripting framework that provides mechanisms for both controlling the execution of different computational tasks and for managing the communication between these tasks using Python object-oriented programming techniques. As we already mentioned, its built-in classes permit one to handle different aspects involved in the workflow as a single unit. All of the advantages coming from object-oriented programming (i.e., extensibility and inheritance) are readily available and allow us to incorporate a third-party scientific code and directly manipulate quantities coming from different codes (Psi4Numpy) in our case.

The Python HLL (high-level language), among others, permits one to formally express complex algorithms in comparatively few lines of codes. This makes it rather straightforward to let PyADF interact with Psi4Numpy native Python API. For the sake of completeness, we want to finally mention that, to accomplish our goal, we first had to port some of the frameworks (specifically XCFun, PyADF, and PyEmbed) to the new Python 3.0 standard (i.e., we forked the original code of the cited packages, and we made them publicly available at refs (111, 112)).

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 nonadditive 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 nonadditive embedding potential. Thus, the geometry and basis set of the active system (in this specific case, a H2O molecule) are parsed at Line 7 and the ground-state wave-function 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 (Lines 14–20). A ground-state calculation of the environment molecule (i.e., an NH3 molecule in this example) is carried out using the 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, which are represented on a common numerical grid, are obtained using get_density() and get_potential() methods (Lines 25, 27), respectively. The PyEmbed module has all of the methods needed to manage the density of both the reference system and the environment to finally compute the nonadditive embedding potential. Indeed, the embed_eval object is instantiated (Line 34) and the nonadditive embedding potential is evaluated on the numerical grid using get_nad_pot (Line 36), once the densities of both the active system and of the environment have been provided.

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 arising 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 , we present its pictorial workflow.

Figure 1

Working flowchart of the uFDE-RT-TDDFT. 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 a grid function object and used to calculate the embedding potential. The active system density is optimized self-consistently according to eq . 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 a grid function. (b) PyEmbed classes are used to calculate the embedding potential. (c) The embedding potential is expressed on the finite basis set representation (GTOs). (d) The active density matrix is evolved according to the real-time-propagation scheme.

We start describing the out-of-loop section. First, the geometry and basis set of the environment are initialized (the orange leftmost block); thus, the ADF package provides, through a standalone single-point calculation, the electrostatic and nuclear potentials of the environment and its density ρII and a suitable integration grid for later use. At this stage, all of the basis sets and exchange-correlation functionals available in the ADF library can be used. In the next step, the green block, the geometry, and the basis set of the active system are parsed from the 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() on the numerical grid is efficiently accomplished using the molecular orbitals (MOs), which requires the valuation of the localized basis functions at the grid points.

Finally, the PyEmbed module comes into play (the 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 nonadditive kinetic and exchange potentials are obtained. The embedding potential is then calculated from its constituents (i.e., the environment electrostatic and nuclear potentials and the nonadditive contribution as detailed in eq ) and evaluated at each grid point, vemb(). The embedding potential matrix representation in the active subsystem basis set, emb, is calculated numerically on the grid aswhere χμ() are the Gaussian-type basis set functions employed in the active systems (used in Psi4Numpy) evaluated at the grid point, . In the above expression, w are specific integration weights.

In the case of an FDE-rt-TDDFT calculation, the electron density of the active system at the beginning of the propagation (t0 = 0, initial condition) is not the ground-state density of the isolated molecule but 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, indicates 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 resulting from the SCF procedure, including the embedding potential, is the starting point for the real-time propagation. Whereupon, at each time step, we determine the embedding potential corresponding to the instantaneous active density (vemb[ρI(t), ρII]). Again, we need its mapping onto the real-space grid as shown in the first green box (box a). Then, we utilize the methods reported in the rectangular orange box (box b) to calculate the nonadditive part of the embedding potential at each grid point. Finally, we add the nonadditive (kinetic and exchange-correlation) potential to the electrostatic potential of the environment calculated again at each grid point. It should be noted that because the density of the environment is frozen, the corresponding electrostatic potential remains constant during the time propagation. In the next phase, its matrix representation in the localized Gaussian basis functions is obtained as in eq (box c, in Figure ). The active system is evolved (box d) using an effective time-dependent Kohn–Sham matrix, which contains the usual implicit and explicit time-dependent terms, respectively, (J[ρI(t)] + Vxc [ρI(t)]) and vext(t), plus the time-dependent embedding potential (emb[ρI(t), ρII]).

For the sake of completeness, the pseudocode needed to evolve the density using the second-order midpoint Magnus propagator is reported in the SI and relies on the methodology illustrated in Section . 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 assessing 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]

Initial Validation and Numerical Stability

Before going into the details of the numerical comparison between our implementation and the FDE-TDDFT scheme based on the linear-response (ADF-LR) formalism, whether in combination with FDE (ADF-LR-FDE) or not, it is important to first assess the basis set dependence of the calculated excitation energies using the two different approaches. This preliminary study is mandatory because Psi4Numpy (Gaussians) and ADF (Slaters) employ different types of atom-centered basis functions. Due to this difference, perfect numerical agreement between the two implementations cannot be expected, but it is important to quantify the variability of our target observables (the excitation energies of a water molecule) with variations in the basis set.

To simulate the linear-response regime within our Psi4-rt, the electronic ground state of a water molecule, calculated in the absence of an external electric field, was perturbed by an analytic δ-function pulse with a strength of κ = 1.0 × 10–5 a.u. along the three directions, x, y, and z. The induced dipole moment was collected for 9000 time steps with a length of 0.1 a.u. per time step, corresponding to 21.7 fs of simulation. This time-dependent dipole moment was then Fourier-transformed to obtain the dipole strength function S(ω), according to eq and the transition energies. The Fourier transform of the induced dipole moment was carried out by means of Padé approximants.[17,117]

As shown in Table , 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, indicating 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 the 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 Moleculea

excitation energy (eV)
	Psi4-rt			ADF-LR
	D	T	Q	D	T	Q
root 1	6.215	6.227	6.224	6.161	6.189	6.287
root 2	7.512	7.466	7.440	7.454	7.465	7.884
root 3	8.363	8.352	8.344	8.309	8.288	8.427
root 4	9.536	8.953	8.651	8.803	8.482	8.628
root 5	9.644	9.572	9.306	8.945	8.845	10.022

Data obtained using TDDFT based on linear-response implemented in ADF (ADF-LR) and the 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 the text for details).

To assess differences in the presence of an environment, we next tested our uFDE-rt-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 nonadditive kinetic and nonadditive exchange-correlation potentials, respectively. The numerical results are reported in Table . Although, as expected, there is no quantitative agreement on the absolute value of the transitions, the shift Δ (Eiso. – Eemb) shows an acceptable agreement for the lowest transitions (for additional excitation energies, see the Supporting Information).

Table 2

Excitation Energies (in eV) Corresponding to the First Five Low-Lying Transitions of Both the Isolated and Embedded Water Moleculesa

excitation energy (eV)
	Psi4-rt-PyEmbed			ADF-LR-FDE
	isolated	emb.	Δ	isolated	emb.	Δ
(a) double-zeta calculations
root 1	6.215	5.817	0.398	6.161	5.687	0.474
root 2	7.512	6.694	0.818	7.454	6.578	0.876
root 3	8.363	7.892	0.470	8.309	7.782	0.527
root 4	9.536	8.768	0.768	8.803	8.336	0.467
root 5	9.644	9.186	0.458	8.945	8.422	0.523
(b) triple-zeta calculations
root 1	6.227	5.796	0.430	6.189	5.689	0.500
root 2	7.466	6.573	0.893	7.465	6.559	0.905
root 3	8.352	7.848	0.503	8.288	7.734	0.554
root 4	8.953	8.560	0.393	8.482	7.969	0.513
root 5	9.572	8.625	0.948	8.845	8.318	0.527

In the embedded water molecule, an ammonia molecule is used as the environment. Data have been obtained using our new Psi4-rt-PyEmbed implementation and the reference ADF-LR-FDE implementation with (a) aug-cc-pVDZ and AUG-DZP basis sets and (b) aug-cc-pVTZ and AUG-TZ2P basis sets (see the text for details). The shift Δ (Eiso. – Eemb) in the transition energies due to the embedding environment is also reported.

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.

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, 125), which correspond to one snapshot taken from an MD simulation. Different cluster models were taken into consideration, by progressive addition of surrounding water molecules (from 1 to 5 molecules) to the single active water molecule. For the active system water molecule propagation in Psi4-rt-PyEmbed, we use the aug-cc-pVDZ basis set, while for the environment, computed using the ADF code, we use the AUG-DZP basis set. In both cases, we use the BLYP[120,121] exchange-correlation functional, while for the nonadditive kinetic and nonadditive exchange-correlation terms in the generation of the embedding potential, the Thomas–Fermi and LDA functionals are used, respectively. In each case, we use 9000 time steps of propagation, which correspond to a simulation of ≈22 fs (time step of 0.1 a.u.). The corresponding dipole strength functions (S(ω) = 2ω/(3π)Im[α(ω)]) along the z-direction are reported in Figure . Upon the increase of the cluster dimension, the lowest-lying transition shifts within a range of about 1 eV, which is consistent with the results for liquid water.[126] It is worth noting that, for this specific transition, many-body excitonic effects (included via the full coupled FDE-rt-TDDFT scheme by Pavanello et al.) are negligible for the energy but are found to be very important to reproduce its spectral intensity.[126] These results give confidence in the numerical stability of the propagation when the number of molecules in the environment is increased.

Figure 3

Dipole strength function S of the water cluster as a function of the number of surrounding molecules (left panel). Right panel: detailed representation of the low-lying transition. The peak corresponding to the isolated molecule is reported for comparison.

The systematic increase of the size of the environment makes it possible to also consider the actual computational scaling of the Psi4-rt-PyEmbed code for this case. To show this scaling, we carried out a single time step of the real-time propagation and broke down the computational cost into those of the different steps in the workflow, as reported in Figure .

In Table and Figure , 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 s 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 Figure (see for instance the tc column of Table ). The time spent in this evaluation depends on the number of numerical integration points used to represent the potential and can be reduced using special grids for embedding purposes once the environment is large enough.

Table 3

Time Usage in Seconds

	ta	tb	tc	td	te
1	0.007	0.29	0.48	0.77	1.75
2	0.01	0.45	0.74	1.20	2.17
3	0.014	0.61	1.0	1.62	2.58
4	0.015	0.74	1.2	1.96	2.89
5	0.02	0.87	1.42	2.32	3.25

Density on the grid (through MOs).

XCFun (nonadditive potential calculation).

Vemb projection onto the basis set.

Total time for Vemb evaluation.

Total time for an rt-iteration.

Figure 2

Time needed for different tasks vs number of surrounding molecules.

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 isolated and using an explicit water cluster to model solvation. To assess the shift due to the embedding potential, we calculated 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 an MD simulation, where the acetone is surrounded by an environment consisting of 56 water molecules. The uFDE-rt-TDDFT calculation was obtained specifying in our Psi4-rt-PyEmbed framework all of 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 the DZP basis set, while for 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, it is located at 3.96 eV. The full absorption spectrum is reported in Figure .

Figure 4

Absorption spectrum of isolated acetone (left panel) and embedded acetone in a water cluster (right panel).

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 20 000 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 acetone and acetone embedded in a water cluster obtained using both of our Psi4-rt-PyEmbed and ADF-LR-FDE methods. The active system response was calculated at the BLYP level of theory, while the Thomas–Fermi and LDA functionals were employed for the nonadditive 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 , we obtain a good agreement in the absolute values, for both isolated and embedded acetone, and the computed shift is likewise in rather good agreement.

Table 4

Isolated and Embedded in a Water Cluster Acetone n → π* Transition, Reported for Both ADF-LR-FDE and Our Psi4-rt-PyEmbed Code

	iso. (eV)	emb. (eV)	ΔE (eV)
Psi4-rt-PyEmbed	3.734	3.958	0.225
ADF-LR-FDE	3.793	3.975	0.182

FDE-rt-TDDFT in the Nonlinear Regime

A specificity in the real-time approach is that the evolution of the electron density can be driven by a 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 the optical control theory; furthermore, it is possible, employing a 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(ω0t), where ω0 is the carrier frequency. In this work, we employ a cos2 shape for the envelope function[127]where τ is the width of the field envelope. We calculated the response of H2O embedded in a water cluster model made of five water molecules (all of the details about the geometry have been reported in a previous section) to a cos2-shaped laser field with carrier frequency ω0 = 1.55 eV (analogous to a Ti:sapphire laser), intensity I = 1.02 × 1014 W cm–2 (which corresponds to a field E = 0.054 a.u.), and a duration of 20 optical cycles. Each cycle lasts 2π/ω0, and the overall pulse spans over 2250.0 a.u. (i.e., 54 fs). The field was chosen along the molecular symmetry axis (z), and the 6-311++G** basis set and the B3LYP functional were used. The propagation was carried out for a total time of 3500 a.u. without any numerical instabilities.

As shown in Figure , the induced dipole does not follow the applied field adiabatically when a strong field is applied; especially in the last few optical cycles, strong diabatic effects are clearly present. These effects lead to the presence of a residual dipole oscillation. Following a previous work on high harmonic generation (HHG) in the H2 molecule,[127] we extracted the high-order harmonic intensities via the Fourier transform of the laser-driven induced dipole moment (neglecting the remaining part, i.e., for t larger than τ, i.e., 2250 a.u. in the present simulation) asIn Figure , 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 were able to observe relatively well-defined peaks up to the 21st harmonics. We mention that this finding qualitatively agrees with data obtained by Sun et al.[20] (see Figure of ref (20)).

Figure 5

Induced dipole moment in the H2O molecule. The representation of the external field is also reported as a green line.

Figure 6

Upper panel: emission spectrum of the isolated water molecule. Lower panel: emission spectrum of the same water molecule embedded in the (H2O)5 cluster.

An important parameter in the analysis of the HHG spectrum is the value of the energy cutoff (Ecutoff), which is related to the maximum number of high harmonics (Nmax ≈ Ecutoff/ω0). In a semiclassical formulation,[128] which, among others, assumes that only a single electron is active for HHG, Ecutoff ≈ I + 3.17U, where I is the ionization potential of the system and U () is the ponderomotive energy in the laser field of strength E and frequency ω0.[128] In the case of molecular systems, the HHG spectra present more complex features and the above formula is not strictly valid. With the laser parameters used here (E = 0.054 a.u., ω0 = 0.05696 a.u.) and the experimental ionization potential of H2O (I = 0.4637 a.u.), the above formula predicts an Ecutoff value of 1.17601 a.u. (Nmax at about the 21st harmonic), which is remarkably consistent with the HHG spectra we observed here.

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 to the introduction of the embedding potential of the surrounding cluster. The latter is consistent with a shift toward lower ionization energy passing from the free water molecule to a small water cluster observed experimentally.[129]

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 used a water–ammonia complex as a general test-case, where the geometry of the adduct was taken from ref (130) and water is the active subsystem.

In the real-time framework, the embedding potential is, evidently, an implicit time-dependent 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 nonadditive terms, which in turn depend 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 to update the embedding potential, which weakly varies in time.

To investigate such a possible speed-up, we carried out different simulations in which the time interval of the embedding potential updating is progressively increased. The results are reported in Table . Of course, as the number of time steps between consecutive 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 Figure , 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 on the order of 0.02 eV.

Table 5

Time in Seconds as a Function of the Number n of Time Steps between Consecutive Updates of the Embedding Potential

ta	tb	tc	n
0.87		94.84	inf (static)
0.87	2.59	97.52	30
0.85	4.32	99.26	20
0.86	8.56	103.31	10
0.86	85.67	180.97	1

Time for Vemb evaluation.

Total time for Vemb evaluation in the propagation.

Total time needed for 100 real-time iterations.

Figure 7

Left: frequency shift in the S function due to the increasing rate of update of the embedding potential. The peaks corresponding to the isolated water molecule are also reported as red trace. Right: expanded view of the homo–lumo transition.

We also reported the partition between different tasks of the time needed for the calculation of the embedding potential in Table . 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 Figure . 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.

Table 6

Time Usage in Seconds

ta	tb	tc	td
0.01	0.33	0.53	0.87

Density on the grid (through MOs).

XCFun (nonadditive potential calculation).

Vemb projection onto the basis set.

Total time for Vemb evaluation.

Figure 8

Comparison of the S dipole strength function obtained with two different integration grids. The violet trace corresponds to the full supramolecular integration grid, while the green trace corresponds to the active system grid. Right: expanded view of the lowest-lying transition.

We furthermore note the possibility of using a 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 ), which makes it possible to neglect points on which these functions have a small value.

Conclusions and Perspectives

In this work, we focused on the implementation of the frozen-density embedding scheme in the real-time TDDFT. We integrated the Psi4Numpy real-time module we recently developed within the PyADF framework. We devised a real-time FDE scheme in which the active density evolves in 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 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 proved to be reliable also in the case of propagation in the nonlinear regime. As a demonstration, we perturbed with a strong electric field a water molecule surrounded by five water molecules acting as a frozen environment. Numerically stable induced dipole moment and the corresponding emission spectrum were obtained.

Finally, we would 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 nonadditive terms. Significant improvement could be achieved by exploiting the sparsity of the matrix corresponding to that difference. Moreover, the use of a smaller integration grid would probably further improve the procedure. Last but not the least, the effect of relaxation of the environment has to be investigated. In our uncoupled FDE-rt-TDDFT scheme, we were 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, to recover supramolecular excitations, and in the nonlinear 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.