Environmental Effects via Frozen Density Embedding in Real-Time Time-Dependent Dirac–Kohn–Sham Theory: Solvation of Lead Halides

TL;DR

This paper introduces a new method to study the electronic behavior of heavy-element molecules in complex environments, like solvents, using advanced computational techniques.

Contribution

The novel contribution is extending real-time time-dependent Dirac-Kohn–Sham theory with frozen density embedding to model environmental effects on heavy-element systems.

Findings

The FDE scheme maintains numerical stability in density matrix propagation for active subsystems.

The method was successfully applied to lead halides in γ-butyrolactone solution, showing solvent effects on absorption spectra.

The approach is suitable for studying electron dynamics in realistic systems under various regimes.

Abstract

Accurately describing
the electronic properties of heavy-element
molecular systems in complex environments is essential for advancing
technologies such as optoelectronics and solar cells. However, achieving
accurate predictions remains challenging because both relativistic
and electron correlation effects must be considered equally, along
with interactions involving other species in the complex environment
(e.g., solvent). This paper extends our real-time time-dependent Dirac-Kohn–Sham
(rt-TDDKS) implementation in PyBERTHA-RT to include environmental
effects using the “uncoupled” Frozen-Density-Embedding
(FDE) scheme, where only the active subsystem evolves dynamically
in time. This adaptation utilizes existing FDE functionality within
the PyEmbed module of the PyADF scripting framework. The native Python
APIs of PyBERTHA-RT and PyADF provide an ideal environment for development,
enhancing readability and reusability. We demonstrate that the FDE
potential maintains the numerical stability of the active subsystem’s
density matrix propagation. Illustrative results for lead halides
(PbCl2 and PbI2) in γ-butyrolactone (GBL)
solution show the effects of increasing solvent molecules on absorption
spectra. This case study demonstrates the new implementation’s
applicability to realistic systems, offering a basis for studying
electron dynamics in heavy-element molecules in complex environments
under linear and nonlinear regimes, relevant to perovskite precursor
chemistry.

Figures

Basic real-time DKS-uFDE flow in PyBERTHA-Embed-RT: (a) numerical representation of ρ̃(r, t′ + dt) on the grid; (b) PyBERTHA-Embed classes are used to calculate the embedding potential; (c) time-evolution step.

Geometric structures used in this work, based on a snapshot of molecular dynamics simulations by Kaiser et al.: (a) PbCl2-GBL system; (b) PbCl2-(GBL)2 system; (c) PbCl2-(GBL)3 system; and (d) PbCl2-(GBL)32 system. For both PbI2 and PbCl2 (not shown), the systems (a)–(c) were obtained from the full system (d) by selecting the nearest GBL neighbors to the respective lead halide.

Absorption spectrum (Fourier Transform of the α xx component) for the PbCl2 molecule with/without (upper/lower panel) embedded with one molecule of GBL. Data are reported using different auxiliary basis sets with increasing accuracy (see text for details).

Absorption spectrum (Fourier transform of the α xx component) for PbCl2 embedded with one molecule of GBL. Data are reported using different G-spinor basis sets (namely, dyall.v2z and dyall.v3z; see text for details). The energy range extends up to 6.2 eV (upper panel) and up to 14.2 eV (lower panel), so that the lower-intensity features appearing between 3 and 6.2 eV can be properly visualized.

Effect on the absorption spectrum due to the interaction of PbCl2 with an increasing number of the GBL molecules (1, 2, 3, and 32 molecules in the panel (a)–(d), respectively). The spectrum of the isolated PbCl2 (denoted by “free”) has also been reported for comparison.

Effect on the absorption spectrum due to the interaction of PbI2 with an increasing number of the GBL molecules (1, 2, 3, and 32 molecules in panel (a), (b), (c), and (d), respectively). The spectrum of the isolated PbI2 (denoted by “free”) has also been reported for comparison.

Comparison of the absorption spectrum of PbI2[1GBL] using the uFDE scheme, with PbI2 as the active system and GBL as the environment (solid blue line), and supermolecular calculation (dashed green line). The spectrum of the isolated PbI2 (solid black line) has also been reported for comparison.

Tables

1. Timing (s) for the PbCl2 Molecule with n Molecules of GBL

n	t	t	t
0	5.4
1	7.6	0.02	1.1
1	(10.1)	(0.04)	(1.1)
2	9.0	0.03	1.9
3	9.1	0.05	2.7
32	9.9	0.49	24.3

Full text

Introduction

1

The interaction between matter and radiation underlies a wide range of phenomena, spanning weak-field processes such as photoexcitation, absorption, scattering, light harvesting in dye-sensitized solar cells, ?,? and photoionization, as well as strong-field effects including high harmonic generation, ?,? optical rectification, ?,? multiphoton ionization,? and above-threshold ionization.? The advent of Free Electron Lasers (FELs) and attosecond techniques ?,? has opened new avenues for real-time studies of electron dynamics and chemical reactions, enabling the observation and potential control of electron motion within molecules. These approaches provide direct insights into bond breaking, ?−? ? bond formation,? and ionization processes ?,? on ultrafast time scales.

Real-time time-dependent electronic structure theory, which directly solves the time-dependent equations of electronic motion, is particularly well suited for studying molecular responses and electronic dynamics. This approach is capable of capturing complex nonlinear phenomena and the effects of shaped external fields, making it essential for applications in quantum optimal control theory.? Significant progress in these techniques has been reviewed by Goings et al.,? with real-time time-dependent density functional theory (rt-TDDFT) emerging as a method of choice due to its favorable balance between accuracy and computational efficiency. ?−? ? ? ? ? ? ? ? ? ? ?

Early rt-TDDFT simulations primarily focused on isolated systems, but electronic and optical properties are often strongly influenced by environmental polarization, making isolated-molecule models insufficient. To address this limitation, rt-TDDFT has been combined with QM/MM ?,? and polarizable continuum models (PCM) ?−? ? to include environmental effects at reduced computational cost. These approaches have subsequently been extended to account for nonequilibrium solvent responses ?−? ? ? and nonequilibrium cavity-field polarization in homogeneous dielectrics.? To move beyond classical descriptions of the environment, Koh et al.? integrated rt-TDDFT with block-orthogonalized (BO) Manby-Miller theory,? enabling faster simulations and the explicit inclusion of solvation effects. This framework was later extended by some of the present authors to study environmental influences on X-ray absorption spectra of halogen anions in water.? A quantum mechanical alternative for incorporating environmental effects is provided by the frozen-density embedding (FDE) scheme, a DFT-in-DFT framework that partitions a larger Kohn–Sham system into smaller, coupled subsystems. ?−? ? By treating all subsystems quantum mechanically, FDE not only offers computational advantages but also provides detailed insight into subsystem properties and intersystem couplings.? Neugebauer ?,? introduced coupled FDE, a subsystem TDDFT formulation based on linear response that eliminated some of the approximations present in earlier TDDFT-FDE implementations and allowed great flexibility in controlling which subsystems are allowed to respond to the time-dependent perturbations. More recently, this approach has been further extended ?,? to address charge-transfer excitations, leveraging an exact FDE framework. ?−? ? ? ? Krishtal et al.? developed a real-time DFT subsystem formulation, extending FDE to rt-TDDFT for the first time by updating the embedding potential and evolving Kohn–Sham subsystems simultaneously. This work showed that, at least within an implementation based on plane waves and periodic boundary conditions using pseudopotentials, the method is numerically stable. An “uncoupled” scheme can approximate the full response of the system by restricting the density response to a single active subsystem while freezing all others. Even in this uncoupled formulation, the embedding potential remains time-dependent in real-time TDDFT-FDE and must therefore be recalculated during the propagation. This approach restricts time evolution to the active subsystem while still effectively incorporating environmental effects. Uncoupled FDE has been shown to reproduce supermolecular results in linear-response TDDFT, even in the presence of hydrogen bonding, provided that excitation couplings between subsystems are negligible. ?,? More recently, De Santis et al.? extended the nonrelativistic real-time time-dependent Kohn–Sham formulation based on localized Gaussian basis set functions and implemented it within the Psi4NumPy framework, ?,? including real-time embedded TDDFT schemes that combine the FDE and BOMME approaches.? Additional recent implementations can be found in refs ?,? .

For molecular systems containing heavy elements, the accurate inclusion of relativistic effects is mandatory. The rt-TDDFT scheme has been extended to incorporate high-level relativistic corrections. In particular Repisky et al.? pioneered relativistic real-time TDDFT for atomic and molecular systems using the four-component Dirac-Kohn–Sham (DKS) method (rt-TDDKS), while nearly concurrently Goings et al.? developed X2C Hamiltonian-based electron dynamics for UV/vis spectral calculations. More recently, a rt-TDDKS implementation (PyBERTHA-RT) leveraging modern software engineering practices, including interlanguage communication between Python, C, and Fortran was presented by some of us. ?,? This implementation is based on the recently developed PyBERTHA, ?−? ? ? which serves as the Python API for the relativistic four-component BERTHA code, ?−? ? and has been subsequently been extended to exploit GPU-based acceleration technology.?

Recently, Olejniczak et al.? reviewed strategies to efficiently include environmental effects into accurate relativistic Hamiltonians via embedding theories. While approaches such as the Polarizable Continuum Model (PCM) can account for such effects within real-time Dirac–Kohn–Sham (rt-TDDKS) simulations, the use of quantum embedding approaches such as FDE in combination with rt-TDDKS has not yet been reported. A noteworthy recent development is an implementation in Turbomole, based on a two-component (2c) formalism for the inclusion of spin–orbit coupling, which has also been extended to periodic systems.?

In this work, we extend the rt-TDDKS method formulated with relativistic G-spinor localized basis functions, to the FDE framework. We take advantage of our recent implementation for ground-state DKS calculations (PyBERTHA-Embed),? specifically in its uncoupled variant (uFDE). A unified Python-based framework has been developed to enable efficient interoperability between the rt-TDDKS implementation in PyBERTHA ?,? and the PyADF API? and PyEmbed module, ?,? which are used to define the basic procedures for the embedding scheme. A detailed description of the PyEmbed module has been recently presented by Focke et al.? For clarity, we refer to the original rt-TDDKS implementation as PyBERTHA-RT, and to its extension to the uFDE subsystem framework presented here as PyBERTHA-Embed-RT.

To demonstrate the capability of the new implementation to tackle large-scale systems, we apply it to investigate the influence of a solvent environment composed of γ-butyrolactone (GBL) molecules on the valence spectra of the perovskite precursors PbCl_2_ and PbI_2_. These systems were previously studied by some of us? using frequency-domain TDDFT calculations based on structures obtained from ab initio molecular dynamics (AIMD) simulations. Here, the focus is on assessing the numerical stability of the implementation and its accuracy with respect to both auxiliary fitting and G-spinor basis sets; consequently only structures from a single AIMD snapshot are considered. Although a single structure is insufficient to yield spectra directly comparable to experiment, it allows us to distinguish solvent effects arising from the nearest neighbors to those associated with the second solvation shell and beyond. To this end, we consider models containing 1, 2, and 3 GBL molecules (the latter corresponding to roughly the average GBL coordination number for the first solvation shell? for PbCl_2_ and PbI_2_) as well as the full system containing all 32 GBL molecules in the simulation supercell.?

The manuscript is organized as follows. In Section, we review the FDE formalism, its extension to rt-TDDKS, and the use of density fitting, with emphasis on code interoperability among PyBERTHA,? XCFun,? PyADF,? and the PyEmbed module. ?,?,?

Section presents the results of solvation effects in the absorption spectra of PbCl_2_ and PbI_2_, and Section provides conclusions and discusses perspectives.

Theory and Implementation

2

In this section, we briefly review the theoretical foundations of the FDE scheme, the basic equations of the Dirac-Kohn–Sham (DKS) method, and its extension to the real-time time-dependent DKS (rt-TDDKS) method as implemented in the BERTHA code using its Python API (PyBERTHA). ?,? We also describe the new rt-TDDKS-uFDE method, focusing on the use of the density fitting scheme, previously applied by some of us to the electronic ground state using the DKS-uFDE scheme,? which accelerates the calculation of the FDE potential matrix representation for the active subsystem. Finally, we illustrate the general implementation strategy, which is greatly facilitated by BERTHA’s Python API. This framework provides flexible code development, taking advantage of Python’s reusability and portability, and enables interoperability with other FDE implementations available through the PyADF and PyEmbed modules.

Subsystem DFT and Frozen Density Embedding

Formulation

2.1

The original concept of subsystem DFT can be traced to the work of Senatore and Subbaswamy,? with a formal derivation provided by Cortona.? Within the linear-response framework, Casida and Wesolowski proposed a TDDFT extension? of the FDE scheme. In 2014, Jacob and Neugebauer presented a detailed theoretical analysis,? and very recently published an update, including new areas of application.? We refer interested readers to these comprehensive reviews.

The whole system is divided into N subsystems, and the total density ρ_tot_(** r ) is expressed as the sum of the electron densities of the various subsystems [i.e., ρ_ a _( r **) (a = 1,..., N)]. In the simplest case of a single subsystem of interest (I) interacting with an environment (II), one can consider the total density as partitioned into two contributions as

$$[eqn]$$

The total energy of the system can then be written as

$$[eqn]$$

with the energy of each subsystem (E _ i [ρ i _], with i = I, II) given according to the usual definition in DFT as

$$[eqn]$$

In the above expression, v nuc ^ i ^(** r **) is the nuclear potential due to the set of atoms which defines the i-th subsystem and E nuc ^ i ^ is the related nuclear repulsion energy. T _ s [ρ i _]­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 expression

$$[eqn]$$

with v nuc ^I^ and v nuc ^II^ 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 nuc ^I,II^ term. The nonadditive contributions are defined as

$$[eqn]$$

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) 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 (first introduced by Wesolowski and Warshel in ref ?.) and leads to a set of Kohn–Sham-like equations, one for each subsystem.

$$[eqn]$$

which are coupled by the embedding potential term v emb ^I^(** r ), which carries all dependence on the other fragment’s density. In this equation, v eff ^KS^[ρ_I_]­( r ) 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). Here $[eqn]$ denotes the kinetic energy operator, which in a nonrelativistic framework has the form −∇^2^/2, whereas for a relativistic framework (Dirac-Kohn–Sham theory) is c α·p (see discussion below). We also note that, in the relativistic framework, the FDE expressions above correspond to the case in which an external vector potential is absent. Further details for their generalization can be found in ref ?. In the framework of FDE theory, v emb ^I^( r **) is explicitly given by

$$[eqn]$$

where the nonadditive 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, in practice, prohibited (as orthogonal KS orbitals are only obtained for each subsystems, not for the whole system). 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 ( $[eqn]$ , 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? or the GGA functional PW91k.? These functionals have been shown to be accurate for weakly interacting systems, including hydrogen-bonded systems, whereas their use for subsystems with a larger covalent character is problematic (see ref ? and references therein). The search for more accurate KEDFs is a key aspect for the applicability of the FDE scheme as a general approach, including the partitioning of systems that also involve breaking covalent bonds.? For a comprehensive review of recent efforts in the development of accurate KEDFs, see ref ?.

In general, the set of coupled equations that arise in the FDE scheme for the subsystems must be solved iteratively. Typically, a “freeze-and-thaw” procedure is employed, where the electron density of the active subsystem is determined while keeping the electron densities of the other subsystems frozen. The former is then frozen in turn when the electron densities of the remaining subsystems are relaxed. This process may be repeated multiple times until all subsystem densities have converged. In this case, the FDE scheme can be viewed as an alternative formulation of the conventional KS-DFT approach for large systems, and by construction, it scales linearly with the number of subsystems.

Within the ground-state SCF or real-time propagation frameworks, the implementation of FDE is efficient, as the v _ emb _ ^ I ^(r) contribution can be treated as an effective one-electron potential added to the Kohn–Sham matrix. Unlike linear-response formulations, which require the explicit evaluation of nonadditive embedding kernels, ?,?,? the real-time scheme captures these interactions implicitly through the time-evolution of the potential. When using localized basis functions, the matrix representation of the embedding potential (V ^emb^) can 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. It should be noted that, regardless of whether one or multiple subsystem densities are optimized, the matrix V ^emb^ must be updated during the SCF procedure because it also depends on the density of the active subsystem (see eq).

Thanks to its flexibility, the FDE scheme has been implemented in various forms in many computational packages, ?,?−? ? ? ? ? ? ? ? ? based on plane waves, Slater-type functions, or Gaussian-type functions. FDE has also been implemented to treat subsystems at the fully relativistic four-component level based on the Dirac equation within the DIRAC code,? and can be used with DFT and various wave function methods for both molecular properties and energies involving ground or excited electronic states. ?,?,?,?,?−? ?

Recently, we extended the full four-component relativistic Dirac-Kohn–Sham (DKS) method, as implemented in the BERTHA code, to include environmental and confinement effects via the FDE scheme (DKS-in-DFT uFDE), see ref ?. The implementation takes advantage of the DKS formulation in BERTHA, which uses density fitting algorithms at the core of the computation (i.e., in the evaluation of the embedding potential and its matrix representation in relativistic G-spinor functions); some details will be given below.

As already mentioned, the FDE scheme has been extended to response theory since many years to access various interesting properties, ?,?,?,?,? such as electronic absorption? and NMR shielding, ?,? for which FDE has been shown to perform well, as these properties are often relatively local. In a response formulation, the embedding potential and its derivatives enter the equations. If more than one subsystem is allowed to respond to external perturbations, ?,?,?,? the derivatives of the embedding potential introduce coupling in the subsystems’ response, just as the embedding potential couples the subsystems’ electronic structures in the ground state.

The theoretical background of subsystem TDDFT has also been reviewed in ref ?, which discusses recent advances in both theory and applications. An important development concerns the real-time versions of subsystem TDDFT, introduced by Krishtal et al.? and implemented using plane waves and periodic boundary conditions. An implementation using localized basis functions was presented by De Santis et al.? and more recently in ref ?. In its simplified “uncoupled” form, one considers only the response to the external applied field of a subsystem of interest (and thus the embedding potential and its derivative with respect to this subsystem’s density). While neglecting the environment response may seem a drastic approximation, good performance relative to supermolecular reference data has been obtained for the excitation energies of a chromophore in a solvent or crystal environment, even when only retaining the embedding potential. ?,? We will employ this framework (uFDE) in this contribution.

The Real-Time TDDKS Method and Its Extension

to FDE Based on Density Fitting

2.2

For the theoretical foundation of the Dirac-Kohn–Sham (DKS) methodology, we refer the reader to previous works ?,?−? ? ? ? ? ? ? and the references therein. Here, we summarize only the main aspects of the implementation of an all-electron DKS method based on G-spinor basis sets and density-fitting techniques, as implemented in BERTHA? and its Python API, PyBERTHA, ?,? which are important for understanding the extension of our real-time TDDKS formulation to the uncoupled FDE scheme (rt-TDDKS-uFDE) in this context.

In atomic units, and considering only the longitudinal electrostatic potential, the DKS equation is given by

$$[eqn]$$

where c is the speed of light in vacuum, p is the electron momentum, and

$$[eqn]$$

where σ = (σ_ x , σ y , σ z ), σ q _ is a 2 × 2 Pauli spin matrix and ** I ** is the 2 × 2 identity matrix. The longitudinal interaction term is represented by a diagonal operator borrowed from nonrelativistic theory and made up of: a nuclear potential term v N(** r ), a Coulomb interaction term v H ^(l)^[ρ­( r )], and the exchange-correlation term v XC ^(l)^[ρ­( r **)]. We note that the Breit interaction contributes to the transverse part of the Hartree interaction and is not considered here, as we restrict ourselves to nonhybrid, nonrelativistic functionals of the electron density.

In BERTHA, the spinor solution (Ψ _ i _(r) in eq) is expressed as a linear combination of G-spinor basis functions,? M μ ^ T ^(** r **) (T = L, S with L and S referring to the so-called “large” and “small” components, respectively). The G-spinors do not suffer from the variational problems of kinetic balance (see ref ?. and references therein) and, for the evaluation of multicenter integrals, retain the advantages that have made Gaussian-type functions the most widely used expansion set in nonrelativistic quantum chemistry.

Recently, the DKS method has been extended to real time, solving the real-time time-dependent DKS equation (rt-TDDKS) to investigate both linear and nonlinear properties in molecular systems. The method has been detailed by Repisky and co-workers ?,?,? and also by some of us; ?−? ? here, we summarize only the main points specific to our implementation in the BERTHA code. The time-dependent equation for rt-TDDKS is conveniently expressed in terms of the Liouville–von Neumann (LvN) equation, which in an orthonormal basis reads

$$[eqn]$$

where ** D **(t) and ** H **(t) are the one-electron density matrix and time-dependent DKS matrix, respectively. The time dependence of ** H **(t) arises from the explicitly time-dependent external electric field and, implicitly, from the time dependence of the density matrix ** D **(t) in the Coulomb and exchange-correlation terms. In BERTHA, the DKS matrix is built in atomic (AT) G-spinor basis set? and defined at each time t as

$$[eqn]$$

The ** H ** _ DKS _ ^ AT ^ matrix is defined in terms of ** V ** ^(TT)^, ** S ** ^(TT)^, and Π ^(TT′)^ matrices being, respectively, the G-spinor basis set representation of the local potential, the overlap matrix, and the matrix of the kinetic operator (the label TT = LL, SS and TT′ = LS, SL, with L and S). The local potential ** V ** ^(TT)^ is given by the sum of four terms:

$$[eqn]$$

the nuclear (** v ** ^(TT)^), Coulomb (** J ** ^(TT)^[ρ­(t)]), and exchange-correlation potential (** V ** _ xc _ ^(TT)^[ρ­(t)]) and external time dependent potential (** v ** _ ext _(t)).

The DKS method enables the variational incorporation of both scalar and spin–orbit interactions. However, several approximations are implicitly introduced when the Hamiltonian is used in the form defined by eq. The most significant approximation in the theory arises from the use of exchange-correlation functionals that depend only on the density, rather than on the relativistic four-current; this is reflected in ** V ** _ xc _ ^ TT ^[ρ]. Research on current-dependent functionals in a relativistic context is still at an early stage. In practice, nonrelativistic density functionals are used, even though they were not explicitly designed for relativistic calculations. Specifically, in our DKS implementation, we use the so-called “density-only” approximation, in which the exchange-correlation term depends only on the charge density (and its gradients), and not on other variables, such as spin density or magnetization,? which may also be used to reparametrise the exchange-correlation potential. It is known that such a restriction limits the achievable accuracy in describing electronic transitions involving, for example, spin-flip mechanisms. Another approximation, also implicit in the exchange-correlation potential of the form ** V ** _ xc _ ^ TT ^[ρ­(t,r) ], is the so-called “adiabatic approximation”, which assumes that the exchange-correlation potential depends only on the density at time t and, by definition, cannot include any memory effects.

As already mentioned, the DKS Hamiltonian (** H ** _ DKS _ ^ AT ^) is evaluated in atomic basis set, whereas the LvN eq (eq 10) is represented in the orthonormal ground-state molecular spinors. The two representations are related by the simple transformation ** H **(t)_ DKS _ = ** C ** ^†^ ** H **(t)_ DKS _ ^ AT ^ ** C **, where ** C ** is the matrix of occupied ground-state molecular spinors (see ref ?).

The time evolution of the density matrix is given by

$$[eqn]$$

where ** U **(t, t 0) denotes the time-evolution operator. The use of eq requires selecting a density matrix that represents the initial condition (at t 0) of the molecular system. In most applications, it is customary to set the initial density as that corresponding to the ground state of the system (for details, see refs ?,? ).

Within a finite time interval, solving the Liouville-von Neumann equation requires calculating the DKS matrix at discrete time steps and propagating the density matrix over time. In the current implementation, we use the Magnus propagator, performing matrix exponentiation exactly through matrix diagonalization. Typically, the Magnus expansion is truncated at the first order, and the time integral is evaluated using numerical quadrature, using the midpoint rule. Provided that the time interval Δt is sufficiently short, the time evolution operator can be approximated as

$$[eqn]$$

This approach exhibits an error proportional to (Δt)^3^. The expression in eq coincides with the so-called modified-midpoint unitary transform time-propagation scheme originally introduced by Li et al.? The ** H ** matrix at time t + Δt/2, where no density is available, is obtained using an iterative series of extrapolations and interpolations at each time. If this predictor/corrector procedure is converged in a self-consistent manner, the second-order midpoint Magnus propagator preserves 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 ingredient in preserving the numerical stability of the propagation with a range of algorithms that can be applied in this context.? We use a particularly stable predictor/corrector scheme, originally proposed by Repisky et al.,? and implemented in the PyBERTHA framework. ?,?,?

From the time evolution of the system, one can obtain both linear and nonlinear properties, depending on the specific shape and strength of the external field. As an example, the dipole strength function S­(ω) is related to the imaginary part of the frequency dependent linear polarizability by

$$[eqn]$$

In the linear-response regime, the components of the frequency dependent tensor are related with the Cartesian components of the electric induced dipole moment (μ̅_ p _)­

$$[eqn]$$

In our implementation,? the external field can be chosen as an impulsive perturbation given by E(t) = kδ(t)n, where n is a unit vector representing the orientation of the field and for the δ-function we adopt its analytic representation as proposed in ref ?. Our implementation supports various field envelope shapes (e.g., Gaussian envelope); however, using an impulsive perturbation in time has the advantage of containing all frequencies (E(ω) = k), which are therefore probed simultaneously. In the case where one chooses to perturb the system along a selected Cartesian axis, p, the diagonal components of the polarizability tensor are proportional to the components of the induced dipole moment,

$$[eqn]$$

that can be easily extracted from our real-time simulation via Fourier Transform of the time-dependent electric dipole moment μ⃗(t). Each Cartesian component p (with p = x, y, z) is given by

$$[eqn]$$

To simulate the absorption spectra of a molecule in a sample with randomly distributed orientations, three independent simulations must be carried out, each with the perturbation applied along one of the three Cartesian axes. To extract nonlinear effects, one can apply a strong field with a specific shape.? For evaluating high harmonic generation spectra using PyBERTHA-RT, see, for instance, ref ?.

The rt-TDDKS scheme summarized above is extended here to the subsystem density functional theory framework, specifically to FDE in the uncoupled scheme (rt-TDDKS-uFDE). We therefore consider the active subsystem, which may be subjected to a time-dependent external field, at the DKS level, while the electron density of the environment is kept frozen at its ground state value. Thus, a LvN-type equation can be solved solely for the active system, coupled to the frozen environment considered at the DFT level. Note that, although the electron density of the environment is kept frozen, the embedding potential is nevertheless time-dependent due to the time dependence of the electron density of the active system itself (ρ_I_ in eq is now time-dependent). The only modification to eq is in the definition of the effective Hamiltonian matrix representation, which now refers to the active subsystem (** H ** _ DKS _ ^ AT,I^(t)). The matrix representation of the embedding potential (** V ** ^emb(TT)^(t)), which accounts for the effect of the environment, is added to the local potential ** V ** ^(TT)^ introduced in eq. With this substitution, the time propagation scheme remains unchanged.

As mentioned above, some of us have already presented the integration of the uFDE scheme into the DKS method for the ground state, using an efficient implementation based on density fitting and prototyping techniques, resulting in the PyBERTHA-Embed code.? Here, we use a similar strategy.

The current version of BERTHA takes advantage of both density fitting algorithms ?,?,? and advanced parallelization techniques ?,?,? for the efficient evaluation of Coulomb and exchange correlation contributions to the DKS matrix. In the G-spinor basis set, the “active” system electron density is defined as

$$[eqn]$$

where D _ μν _ ^(TT)^ represents the density in the G-spinor basis and ρ_ μν _ ^(TT)^(** r **) are the G-spinor overlap densities. The label TT = LL, SS indicates that L and S refer to the “Large” and “Small” components of the G-spinor basis, respectively.

This quantity (which is a real scalar function) can be accurately approximated and linearized using a set of N_aux_ auxiliary functions.

$$[eqn]$$

The expansion coefficients d _ t _ are worked out as the solution of the linear system, given by

$$[eqn]$$

where ** A ** is a real and symmetric matrix, representing the Coulomb interaction in the auxiliary basis, A _ st _ = ⟨f _ s _ ∥ f _ t _⟩ while the elements (v _ s ) of the vector ** v ** are the projection of the electrostatic potential on the fitting functions, ⟨f _ s _ ∥ ρ I _⟩.

For the exchange-correlation ** K ** matrix, we adopt a similar strategy by solving for ** z ** in the linear system

$$[eqn]$$

where vector ** w ** is the projection of the exchange-correlation potential (ṽ xc ^(l)^[ρ̃(** r **)]) on the fitting functions

$$[eqn]$$

The elements of vector w, which involve integrals of the exchange-correlation potential, are computed numerically by the integration scheme already implemented in the code.? Once the vectors ** d ** and ** z ** have been worked out, the Coulomb and the exchange-correlation contributions to the DKS matrix can be evaluated in terms of 3-center two-electron integrals I _ s,μν _ ^(TT)^ = ⟨f _ s _ ∥ ρ_ μν _ ^(TT)^⟩

$$[eqn]$$

The extension of this scheme to include the contribution of the FDE potential, v emb ^I^[ρ_I_, ρ_II_]­(** r **), is relatively straightforward. The evaluation of the embedding potential matrix (Ṽ _ μν _ ^ emb(TT)^) can strictly follow a procedure similar to that employed for the exchange-correlation? term, and it is given by

$$[eqn]$$

The matrix elements, Ṽ _ μν _ ^ emb(TT)^, is expressed as a linear combination of the three-index coulomb repulsion integrals as?

$$[eqn]$$

The expansion coefficients (c _ t ) are the elements of vector ** c **, solution of the linear system ** Ac ** = ** g **, whre vector ** g ** is the projection of the embedding potential, v emb ^I^[ρ̃_I, ρ_II_]­(** r **), on the fitting functions.

$$[eqn]$$

The embedding potential contribution is efficiently evaluated in a single step together with the Coulomb and the exchange-correlation ones, and finally added to the DKS matrix, J̃ _ μν _ ^(TT)^ + K̃ _ μν _ ^(TT)^ + Ṽ _ μν _ ^ emb(TT)^ = ∑_ t = 1_ ^ N aux ^ I _ t,μν _ ^(TT)^(d _ t _ + z _ t _ + c _ t _). As in the electronic ground state case, where changes in the active subsystem density require updating ** V ** ^emb^ during the SCF cycles, the time propagation of the electron density introduces a time dependence in ** V ** ^emb^ even though the environment densities remain frozen due to the use of the uncoupled scheme. This means that the embedding potential matrix, ** V ** ^emb^, must be updated at each time step. We will show that the numerical noise associated with constructing the embedding potential in this scheme does not affect the numerical stability of the propagation.

Fast Prototyping and Implementation

2.3

Here, we outline the computational strategy we have adopted to implement the rt-TDDKS-uFDE scheme. The developed Python program PyBERTHA-Embed-RT and the related modules (PyBERTHA, PyBERTHA-Embed and PyBERTHA-RT) are freely available under GPLv3 license at ref ?. A data set collection of computational data including geometries of the systems, numerical data, parameters and job input instructions used to obtain the results of Section, is available and can be freely accessed at the Zenodo repository, see ref ?.

In our workflow, we have combined the real-time TDDKS reference procedure, thus the PyBERTHA-RT module, ?,? already implemented within the BERTHA code thanks to its Python API PyBERTHA, ?,?,? with the FDE computational core that we developed in the PyBERTHA-Embed module. ?,? Both PyBERTHA-RT and PyBERTHA-Embed modules are described in detail in ref ?. and ref ?., respectively.

The PyBERTHA-RT code is based on PyBERTHA, a Python? API that has contributed to improving both the usability and interoperability of the BERTHA code. ?,?,?,? All the basic kernel functions written in FORTRAN are collected in a single Shared Object (SO) capable of performing both serial and parallel calculations using both OpenMP-based operations? or hybrid OpenMP/OpenACC ones (i.e., CPU and GPU).? A FORTRAN module, named bertha_wrapper, contains a class implementing all the methods needed to access the basic quantities, such as energy, density, DKS and overlap matrices, as well as other functions that can be easily implemented. Finally, the main PyBERTHA? module has been developed using the ctypes Python module. This module provides C-compatible data types and allows calling functions collected in shared libraries. Thanks to this Python API, the Magnus propagator scheme in PyBERTHA-RT is written entirely in Python, using the numpy or cupy modules for the required linear algebra operations ?,? (matrix–matrix multiplication and diagonalization).

The PyBERTHA-Embed module is based on PyBERTHA and on general PyEMBMOD module. A specific class inside the PyEMBMOD module allows to well isolate all the FDE data and operations increasing the level of abstraction. The module is used to manage all the required quantities for the generation of the embedding potential, that is v emb[ρ̃_ I , ρ II _(** r **)]. It has been engineered in a such manner that all details of the FDE low-lying implementation are completely transparent from the PyBERTHA side. This has the advantage that all future developments and/or integration of the FDE scheme (g.e. using DKS theory also for the environment DKS-in-DKS FDE) will not affect the PyBERTHA-Embed code, i.e., it will remain completely unchanged. In particular in this first version, the PyEMBMOD module can handle the basic procedures which are based on the use of PyADF, ?,?,? PyEmbed module, ?,?,? and the XCFun library? to evaluate nonadditive exchange-correlation and kinetic energy contributions on user-defined integration grids. This approach gave us both the advantages of the code reusability and clearly made the debugging phase in the development of software straightforward. Finally, we mention that PyADF is not specific to a single program, but works with a number of different quantum chemistry codes, in this context we use it with ADF program.?

An example of the interoperability achieved between different modules is presented in Figure, where a basic workflow of our rt-TDDKS-uFDE implementation in PyBERTHA-Embed-RT is shown. The code accepts several parameters that define the molecular fragments and specify the details of the calculations: geometry, G-spinor basis set, and exchange-correlation functional used in the DKS calculation for the active system. Similarly, it requires the geometry, basis set, exchange-correlation functional, and other details for the environmental system, which, in this context, is handled using the ADF code via PyADF. It also takes the information required for real-time propagation. All of these options are already available in the PyBERTHA-RT code. These include, for example, the definition of the external electric field function, its direction, the total simulation time, and the threshold for the predictor/corrector scheme used in the Magnus propagator. An additional option, which can now be set in a PyBERTHA-Embed-RT run, specifies the frequency at which the embedding potential is updated during propagation. In all numerical examples reported below, we update the embedding potential at each time step of the propagation.

Basic real-time DKS-uFDE flow in PyBERTHA-Embed-RT: (a) numerical representation of ρ̃(r, t′ + dt) on the grid; (b) PyBERTHA-Embed classes are used to calculate the embedding potential; (c) time-evolution step.

The electron density of an active system is obtained via a DKS calculation using PyBERTHA while the electron density, the Coulomb potential and nonadditive terms of the environment are managed using PyADF. These quantities are mapped on a common numerical grid and used in PyEMBMOD to evaluate the embedding potential. We refer the interested reader to refs ?,? , where some of us presented the workflow implemented in PyBERTHA-Embed.

The complete workflow can be summarized as follows: there is an initial phase in which: (i) a single self-consistent DKS calculation is performed on the active system to obtain an initial electron density (ρ̃_ I ); (ii) the PyEMBEDMOD module initializes and manages the calculation of the frozen environment, using PyADF in combination with the ADF code. This module also defines the numerical grid that will be used to represent the key quantities in the procedure. In this phase, the SCF calculation is carried out on the environment systems, and the first two terms of the embedding potential (the nuclear potential and the Coulomb potential of fragment II, see eq), which remain unchanged during the propagation, are mapped onto the numerical grid. Once the electron density of the active system (ρ̃ I (r)) has been calculated at the DKS (or rt-TDDKS) level, then its G-spinor representation is converted to the numerical representation on the grid (i.e., getdensityongrid in Figure). Then we evaluate the corresponding nonadditive embedding potential and return the embedding potential (v _ emb (r _ i )) on the grid as a numpy array (i.e., getpotential in Figure). Afterward, starting from the numerical representation v _ emb (r _ i ) and the numerical grid, we build its G-spinor matrix representation (Ṽ _ μν _ ^ emb(TT)^) by projecting the embedding potential onto the density fitting basis set via numerical integration (see eq and eq) (see setembpotongrid in Figure).

The Ṽ _ μν _ ^ emb(TT)^ matrix is then added to the full ** H ** _ DKS _ ^ AT ^ and used to converge the ground state DKS calculation and subsequently during the time propagation of the rt-TDDKS-uFDE. Note that in this latter case, the procedure described above (see also Figure, panel b) is performed at each time step of the propagation because the electron density of the active fragment (ρ_ I (t)) is time dependent, which is reflected in the time dependence of the embedding potential (v emb ^I^[ρ_I(t), ρ_II_]­(** r **)).

Absorption

Spectra of PbX2-GBL systems (X = Cl, I)

3

Here, we use the newly developed PyBERTHA-Embed-RT code to simulate the absorption spectra of PbCl_2_ and PbI_2_ molecules interacting with an increasing number of solvent molecules, up to 32, specifically γ-butyrolactone (GBL). The solution chemistry of PbCl_2_ and PbI_2_ is relevant to the synthesis of lead-halide perovskites, which are widely investigated for applications in solar cells and other electronic devices.? Solution processing typically involves mixing a lead-halide precursor with a halide salt in a strongly coordinating solvent (e.g., GBL or others), followed by the formation of the perovskite solid phase upon solvent elimination. Chlorine-containing precursors are known to increase perovskite grain size and crystallinity, leading to more homogeneous perovskite thin films and improved charge-carrier transport. The different behaviors of PbCl_2_ and PbI_2_ in solution therefore contribute to controlled perovskite synthesis. A detailed account of these precursors in solution was provided by Kaiser et al.,? combining ab initio molecular dynamics, experimental UV–vis absorption measurements, and frequency-domain TDDFT simulations employing the B3LYP functional. One key finding from this prior study is that PbCl_2_ is generally less solvated than PbI_2_, consistent with a higher lead-halide bonding energy in PbCl_2_: chlorine acts as a stronger ligand and limits coordination by solvent molecules.

Before presenting our results, we emphasize that the primary goal of these simulations is to demonstrate the numerical stability of the new implementation and its potential applicability to complex chemical systems. The present calculations are not intended to enable a strict comparison with experimental UV–vis measurements. Two main factors limit quantitative agreement: (i) the use of GGA exchange-correlation functional (due to current implementation constraints), which is expected to underestimate the transition energies; (ii) the use of a single AIMD snapshots, which neglects nuclear dynamics and thermal averaging. The error associated with approximation (i) may be reduced by using hybrid or range-separated hybrid functionals, which are known to yield more accurate transition energies, while the limitation in approximation (ii) can be alleviated by averaging the absorption spectra over multiple nuclear configurations. Furthermore, as we use FDE in the uncoupled scheme, we neglect the polarization response of the environment upon excitation/perturbation, cannot accurately describe electronically excited states with significant contributions from both the active subsystem and the environment. Within uFDE, the smaller the active subsystem, the less one can expect to reproduce all features of a supermolecular calculation, as only a limited portion of the occupied and virtual spaces is represented.

The geometries used in this study were extracted from a snapshot of the molecular dynamics simulation reported in the literature by some of the authors? and are provided explicitly in the referenced data set.? As an example, Figure shows the structures for PbCl_2_ interacting with 1, 2, 3, and 32 molecules of GBL. For the active systems (PbCl_2_ and PbI_2_) we used DKS level of theory and Dyall’s basis sets, specifically dyall.v2z and dyall.v3z,? which are designed for all-electron relativistic calculations for the whole periodic table. Alongside these orbital basis sets, we used the auxiliary basis set Gen-n2-v1 and evaluated the effect of two more accurate auxiliary fitting basis sets, Gen-n3-v1 and Gen-n3-v2. These auxiliary basis sets were generated through an algorithm we have recently developed in ref ?. For all real-time simulations, we used 20 000 time steps with a time step of 0.15 au for a total time of 3000 au (72 fs). An analytic δ­(t) was used as an external perturbation with a field intensity of 5 × 10^–5^ a.u. In all cases, the PBE exchange-correlation functional was applied within the density-only approximation (see above for more details). The embedding potential, originated by the GBL solvent molecules, has been updated at each time step of the real-time evolution and has been obtained using PyADF (and the PyBerthaEmbed module) in combination with ADF packages? and the TZ2P basis set. The absorption spectra, related to the dipole strength function S(ω) (eq), have been obtained using the Fourier Transform of time-dependent induced dipole moment of the active system molecule using Pad’ approximation,? in combination with a Gaussian damping function with an exponent of 1 × 10^–6^.

Geometric structures used in this work, based on a snapshot of molecular dynamics simulations by Kaiser et al.: (a) PbCl2-GBL system; (b) PbCl2-(GBL)2 system; (c) PbCl2-(GBL)3 system; and (d) PbCl2-(GBL)32 system. For both PbI2 and PbCl2 (not shown), the systems (a)–(c) were obtained from the full system (d) by selecting the nearest GBL neighbors to the respective lead halide.

As previously mentioned, density fitting is central to our implementation. It is used in the evaluation of the Dirac-Kohn–Sham matrix during the initial SCF procedure, throughout the real-time evolution, and in the evaluation of the G-spinor matrix representation of the FDE potential. Therefore, we first assess the impact of using different density fitting auxiliary basis sets of increasing accuracy. Data are reported for PbCl_2_ (active system) interacting with a single molecule of GBL (environment) in Figure (upper panel). For a consistent comparison, we have also included the spectrum obtained for the free PbCl_2_ molecule at the same geometry (see Figure, lower panel). The inclusion of a single molecule of GBL via uFDE scheme reduces the symmetry and increases the number of electronic states, resulting in a much richer spectrum than that of the isolated PbCl_2_. In the low-energy region, the solvent induces a pronounced blue shift of the first transition (approximately 0.5 eV) for the specific configuration considered here. It is noteworthy that the effect of the auxiliary fitting basis sets is almost negligible for the isolated PbCl_2_. For this system, the three auxiliary fitting basis sets produce nearly indistinguishable spectra, and the smallest auxiliary basis set used (Gen-n2-v1) is as accurate as the larger one, Gen-n3-v2. However, for the embedded system, the smallest auxiliary basis set tends to shift the spectrum slightly to higher energy (by about 0.05 eV) compared to the reference data (obtained using Gen-n3-v2). This small discrepancy in the spectrum suggests that, to maintain high accuracy, the inclusion of the embedding potential requires an auxiliary fitting basis with sufficient flexibility to accurately project the embedding potential (see eq).

Absorption spectrum (Fourier Transform of the α xx component) for the PbCl2 molecule with/without (upper/lower panel) embedded with one molecule of GBL. Data are reported using different auxiliary basis sets with increasing accuracy (see text for details).

In Figure, we examine the effect of the G-spinor basis sets on the resulting spectrum for the PbCl_2_ molecule in the presence of one solvent GBL molecule. The effect, passing from dyall.v2z to dyall.v3z, is relatively small at low energies (see upper panel), with differences below 0.05 eV, while at energies above 5 eV (see lower panel), the effect becomes more significant; the spectrum obtained using dyall.v3z basis clearly shows more electronic states due to its increased flexibility. Nevertheless, it is important to note that the numerical stability of our real-time DKS simulation, including the time-dependent embedding potential, is not compromised when using a more flexible and significantly larger principal basis set.

Absorption spectrum (Fourier transform of the α xx component) for PbCl2 embedded with one molecule of GBL. Data are reported using different G-spinor basis sets (namely, dyall.v2z and dyall.v3z; see text for details). The energy range extends up to 6.2 eV (upper panel) and up to 14.2 eV (lower panel), so that the lower-intensity features appearing between 3 and 6.2 eV can be properly visualized.

To demonstrate the applicability of our method to more complex systems, we investigated how the absorption spectra of PbCl_2_ and PbI_2_ evolve as the number of GBL solvent molecules increases, considering n = 1, 2, 3, and 32 (see Figures and ?). All spectra were obtained from the isotropic component of the frequency-dependent dipole polarizability (eq). For each system, three independent simulations were performed, applying perturbations along the three Cartesian directions, and induced dipole components were combined to obtain the isotropic response. The all-electron dyall.v2z basis set, in combination with the Gen-n2-v2 auxiliary fitting set,? was employed. In all cases, the simulations were found numerically stable and well converged.

Effect on the absorption spectrum due to the interaction of PbCl2 with an increasing number of the GBL molecules (1, 2, 3, and 32 molecules in the panel (a)–(d), respectively). The spectrum of the isolated PbCl2 (denoted by “free”) has also been reported for comparison.

Effect on the absorption spectrum due to the interaction of PbI2 with an increasing number of the GBL molecules (1, 2, 3, and 32 molecules in panel (a), (b), (c), and (d), respectively). The spectrum of the isolated PbI2 (denoted by “free”) has also been reported for comparison.

For the PbCl_2_ series, the lowest-energy transition exhibits a blue shift relative to the isolated molecule but does not appear to vary systematically with the number of solvent molecules (0.2, 0.6, 0.15, and 0.35 eV for 1, 2, 3, and 32 GBL molecules respectively). The second feature is separated from the first one by roughly the same energy in the 1GBL and 2GBL systems, but their splitting increases when going from 2 to 3 GBL molecules, and further increases from 3 GBL to 32 GBL molecules. For the third feature, appearing near 3.8 eV in the isolated molecule and which is significantly more intense than the first two, we observe a red shift for the 1GBL and 3GBL and a blue shift for the 2GBL one, roughly in line with how the peak maxima the theoretical results shown in Figure 4 of ref ?. change with the number of GBL molecules. The blue shift of this feature is the most significant for the 32GBL case. Taken together, these changes in spectra highlight the importance of incorporating both first- and outer solvation shells, even if represented through an effective embedding potential, as done here.

In the PbI_2_ series, consistent with experimental observations,? electronic transitions occur at lower energies (Figure) than for PbCl_2_. The first excited state of the isolated PbI_2_ molecule appears around 2.0 eV, compared with approximately 3.0 eV for PbCl_2_. Solvation induces an overall shift of approximately 0.3 eV, with PbI_2_(3GBL) exhibiting a slight red shift. Notably, the lowest-energy transition for PbI_2_(1GBL) nearly coincides with that of PbI_2_(32GBL), but there are significant differences between the isolated and solvated systems for energies above 3.0 eV, not only in terms of peak positions but also in intensities. Given that in all calculations the active subsystem has the same structure, these results again highlight the importance of incorporating molecules in the first and second solvation shells, albeit through an effective potential, as done here.

For completeness, Figure compares the absorption spectrum of PbI_2_[1GBL] obtained using the uFDE (PbI_2_ active; GBL environment) with that from a supermolecular calculation. The agreement is reasonably satisfactory for the lowest-energy states and those above 3.5 eV. In contrast, the spectral region between 2.5 and 3.4 eV shows significant deviations in both peak positions and intensities. These differences are likely due to the approximations inherent in the current uncoupled FDE setup and indicate that the chosen active subsystem is too small to describe electronic states that are delocalized over both the active system and the environmental fragment. This observation suggests that a more realistic description would require the inclusion of at least the first solvent shell, with the remaining environment treated within the uFDE framework. At present, such an approach is computationally demanding and would likely benefit from the GPU-accelerated version of the code presented in ref ?.

Comparison of the absorption spectrum of PbI2[1GBL] using the uFDE scheme, with PbI2 as the active system and GBL as the environment (solid blue line), and supermolecular calculation (dashed green line). The spectrum of the isolated PbI2 (solid black line) has also been reported for comparison.

Before concluding, it is interesting to comment on the computational cost associated with incorporating the FDE scheme in our real-time DKS implementation. The timings reported here were obtained using the dyall.v2z basis in combination with the GEN-n2-v1 auxiliary basis set (with results obtained with the GEN-n2-v3 auxiliary basis reported in parentheses), on an Intel­(R) Xeon­(R) CPU E5–2683 v4­(2.10 GHz) processor using 32 threads. Table summarizes the total wall time per simulation time step. The computational cost associated with the propagation of the active subsystem time propagation (t ^ a ^) is largely independent of the number of environmental molecules, and increases only slightly due to the larger total number of grid points required to evaluate the active system density on the full grid (see t ^ b ^). In our implementation, a global grid is defined for the entire system (active subsystem plus environment), thus, the total number of grid points increases with the number of GBL molecules included in the environment. As previously mentioned, the mapping of the environment’s electron density and Coulomb potential, obtained from ADF, onto the numerical grid is performed once at the beginning of the simulation. The wall time required for this initialization step is also reported in Table (t ^ c ^).

1: Timing (s) for the PbCl2 Molecule with n Molecules of GBL

Conclusions and Perspectives

4

In this work, we presented an extension of the relativistic four-component real-time Dirac-Kohn–Sham (rt-TDDKS) method to include environmental effects via the Frozen Density Embedding (FDE) scheme. This was achieved by developing a unified framework that integrates the Python APIs of PyBERTHA and the PyADF scripting framework, resulting in an efficient and interoperable computational tool, namely PyBERTHA-Embed-RT.

Our implementation uses an uncoupled FDE approach in which only the active subsystem evolves in time, while the environment remains frozen in its ground state. We have shown that incorporating the FDE embedding potential through a density-fitting scheme does not compromise the numerical stability of the time-propagation algorithm. This was validated through systematic tests on the absorption spectrum of PbCl_2_ in the presence of the frozen density embedding potential, which produced robust and well-converged results with respect to both the primary G-spinor and auxiliary fitting basis sets.

As a practical demonstration, we investigated the absorption spectra of lead halide perovskite precursors, PbCl_2_ and PbI_2_, treated as active systems and solvated by increasing numbers of GBL molecules. The simulations successfully captured several solvent-induced effects, including the observed blue shift in the lower-energy bands as the number of solvent molecules increased. These results highlight the capability of the method to provide insights into the solution-phase chemistry of complex systems where both relativistic and environmental effects play an essential role.

By construction, the uncoupled FDE approach implemented here neglects the polarization response of the environment upon perturbation and, therefore, cannot accurately describe electronically excited states with significant contributions arising from changes in the environment. A comparison with a supermolecular calculation for PbCl_2_(1GBL) reveals features in the absorption spectrum that are not reproduced by the uFDE calculations with PbCl_2_ as the active subsystem. This suggests that, for a more accurate description, at least a solvation shell should be included in the active system of the uFDE scheme.

Looking ahead, the present framework provides a solid foundation for several promising extensions. A natural next step is the development of a coupled FDE scheme, in which the environment is allowed to respond dynamically to external fields. Such an approach would enable the study of more complex phenomena, such as intersubsystem charge-transfer excitations. Moreover, the proven stability of the real-time approach makes it particularly suitable for exploring the nonlinear optical properties of molecules in solution, an area of significant scientific interest.

The predictive power of the current implementation is also limited by the approximations used in the exchange-correlation functional, which restrict the quantitative accuracy of the computed absorption spectra relative to experiment. Future developments incorporating exact exchange as well as hybrid or long-range-corrected functionals are expected to improve predictive accuracy and fully exploit the relativistic effects, including spin–orbit coupling, offered by the DKS formalism. Furthermore, extending the present implementation to fully exploit GPU acceleration, which we have already demonstrated in the rt-TDDKS context,? would be a valuable direction for improving performance and enabling larger-scale simulations. In conclusion, PyBERTHA-Embed-RT is a powerful and flexible computational framework and represents an important step toward accurate real-time simulations of molecular systems containing heavy elements embedded in complex environments.