Real-Time Propagation TDDFT and Density Analysis for Exciton Coupling Calculations in Large Systems

Photoactive systems are characterized by their capacity to absorb the energy of light and transform it. Usually, more than one chromophore is involved in the light absorption and excitation transport processes in complex systems. Linear-Response Time-Dependent Density Functional (LR-TDDFT) is commonly used to identify excitation energies and transition properties by solving the well-known Casida’s equation for single molecules. However, in practice, LR-TDDFT presents some disadvantages when dealing with multichromophore systems due to the increasing size of the electron–hole pairwise basis required for accurate evaluation of the absorption spectrum. In this work, we extend our local density decomposition method that enables us to disentangle individual contributions into the absorption spectrum to computation of exciton dynamic properties, such as exciton coupling parameters. We derive an analytical expression for the transition density from Real-Time Propagation TDDFT (P-TDDFT) based on Linear Response theorems. We demonstrate the validity of our method to determine transition dipole moments, transition densities, and exciton coupling for systems of increasing complexity. We start from the isolated benzaldehyde molecule, perform a distance analysis for π-stacked dimers, and finally map the exciton coupling for a 14 benzaldehyde cluster.


Introduction
In the last decades the interest to use the natural sun light for an energy transition towards green and clean energy sources has increased. Researchers have focused their investigations on the design of new devices to harvest and use this absorbed light. 1,2 The challenge is to create new light harvesting complexes that can transfer the light energy to a desired reaction center, 3 such as natural light harvesting complexes (LHC), [4][5][6][7] or to design new emitters 8,9 for light and flexible organic light emitting diodes (OLEDs) [10][11][12] with the quantum efficiency close to unity. Such a quantum efficiency means that all the light energy is transferred to the reaction center for the former, or that all input energy is emitted in the form of light without dissipation in OLEDs.
These molecular systems usually contain a large number of chromophores, i.e. molecules that can absorb a photon in the UV-Visible region. Their capacity of absorbing and transferring the corresponding energy are the key factors that determine the efficiency of the system.
From the theoretical point of view, characterization of excited states can be done by solving the eigenvalue problem (Schrödinger equation), either by using a multiconfigurational ansatz 13 or the Kohn-Sham scheme, in the time-dependent density functional theory (TDDFT). 14,15 An excited state is characterized by its transition dipole moment and the corresponding transition density. The former gives the information about the probability of exciting an electron from the ground state and the most efficient polarization direction of the light, while the transition density informs about the change of the electronic density from the ground state to the excited state. Both of these properties are extremely sensitive to the polarization induced by the environment, and the spectroscopic fingerprint of the isolated molecule in vacuum is usually not sufficient even for a qualitative description. 16,17 The common way to obtain the transition dipole moment for a given excitation in the TDDFT framework is to solve the linear-response Casida's equation. 18 For that, one needs to compute the ground state Kohn-Sham (KS) wave function. In addition, this procedure requires a large number of well converged unoccupied (or virtual) KS states in order to have a good representation of the pair-hole orbital transition space. Nevertheless, this procedure has two main drawbacks: i) large molecules / complex systems require a large set of KS orbitals which makes the solution of the Casida's equation unfeasible, ii) high energy excitations are usually not properly described due to the lack of well converged high energy virtual KS orbitals.
Several approaches have been developed to reduce the number of electron-hole pair transitions needed to solve the Casida's equation (e.g. subsytems TDDFT, 19 the simplified Tamm-Dancoff density functional approach, 20 simplified TDDFT, 21 tight-binding approaches to TDDFT (TD-DFTB, 22 TD-DFT+TB 23 ), but they make some approximation on the environment interaction and keep failing for high energy excitations. Real-Time Propagation TDDFT (P-TDDFT) [24][25][26][27] scheme offers the possibility to obtain all frequency excitations at the same cost having converged just the occupied KS states for the ground state. Besides, propagation of the KS states can be highly parallelized enabling to compute optical properties for up to several thousands of atoms. 28,29 These features make P-TDDFT the suitable theoretical framework to study photo-active complex systems such as natural LHC and OLEDS. It is well known that the absorption cross-section can be computed from the propagation of the KS states in a linear-response regime, and recently Kummel et at. 30 provided an analytical form to obtain the oscillator strengths, excitation energies and transition densities from the real time propagation of the electron system. They used the many-body ground state Hamiltonian to show that the a quantum mechanical solution for the response functions after a boost excitation have the sine cardinale form and they exemplified their method by performing real-time propagation TDDFT calculations. In that paper, however, they did not address the way how transition densities can be extracted for randomly oriented distorted chromophores in complex systems.
In this work, we provide another derivation using the linear-response formalism, which makes possible to study exciton couplings of arbitrary complex systems. We present a theoretical description of bright excitations by performing P-TDDFT that can be applicable even to very large systems, all treated entirely at the same atomistic level of theory. Combining the local density analysis we recently introduced 28 and our new derivation described below, the individual transition dipole moments and transition densities are obtained from P-TDDFT. The techniques we propose here permit to compute not just the first excitation but a broad energy range of excited states of molecules taking into account all effects induced by the rest of the system. This methodology provides a powerful tool to study exciton dynamics in light-harvesting complexes 31 and light-emitting layers of OLEDs. 32 We illustrate the validity of our method for the cases of increasing complexity including the benzaldehyde molecule, dimer and a cluster of 14 molecules. In the following we first describe the methodology developed and then discuss the results for the systems listed.

THEORETICAL DEVELOPMENT
In this section we first describe the fundamentals of computation of absorption spectra for finite systems from real-time propagation TDDFT (P-TDDFT). Then, we expose our time-dependent local-density analysis that enables to decompose the absorption spectra into individual contributions.
In the following subsection, we derive analytical expressions for the transition moment and transition densities for a given excitation from P-TDDFT in the linear-response regime.
In the third subsection, both techniques are combined in order to study exciton transfer in multichromophore systems taking into account the actual environment, treated entirely at the same level of theory, i.e. TDDFT.

Absorption Spectra Decomposition
In the dipole approximation, the influence of the electric field applied to a quantum system can be supervised by the time evolution of the dipole moment: where t 0 is the initial time just before the application of the time dependent external field, E λ (t), and the dynamic polarizability element, α νλ (t−t ), is the retarded response that relates the λ component of the external electric field to the change of the ν component of the dipole moment of the system.
In the linear-response theory, the dynamic polarizablity tensor is constructed from the general expression for the retarded response function: 15 where the dipole moment operator and the external potential are described asμ ν = er ν and δv ext λ (t) =r λ E(t), respectively (where e refers to the electron charge), Θ(t − t ) is the Heaviside function and |Ψ 0 is the the wave function at the initial moment of time (t 0 ).
The polarizability tensor α(ω) can be obtained by propagation of the KS states after the application of a perturbative boost. Usually, for finite systems, this external pertubative potential is applied in the form of a delta-pulse of the electric field: 14,15 v ext = −er · kδ(t 0 ), where k corresponds to the magnitude vector of the electric field.
Then, the polarizability tensor can be written in frequency space (ω) in terms of the Fourier Transform of the time-dependent dipole moment as Note here, that in the time-dependent Kohn-Sham scheme, the time-dependent dipole moments can be easily obtained from the time-dependent electronic density (ρ(r, t))): where where Φ i (r, t) are the KS orbitals and the summation runs over the N occ occupied orbitals.
At the same time, by describing the dynamic polarizability tensor in the Lehmann's representation and using the fluctuation-dissipation theorem, 15 the imaginary part of the dynamic polarizability tensor can be related with the oscillator strength of each transition, in atomic units, as where the |Ψ n is the wave function of the nth excited state, Omega n is the corresponding excitation energy and f n n is the oscillator strength which gives the transition probability for that particular excitation. This expression relates directly the dynamic polarization tensor (α(ω)) and the transition dipole moment ( µ n0 = Ψ 0 |r|Ψ n ).
In addition, from the imaginary part of the dynamic polarizability tensor we can obtain the photo-absorption cross-section, which is gives the absorption spectrum of the system: where c is the speed of light and S(ω) = ∞ n=1 f n δ(ω − Ω n ) is the dipole spectral function.

Local Density Analysis for Spectrum Decomposition
The first time we introduced this methodology was in the decomposition of theoretical spectra of the major light harvesting complex II (LHCII). 28 This procedure based on the DFT fundamentals, in particular on the DFT theorem exposed by Hohenberg and Kohn, that establish the electron density as a basic variable from which any observable property of a quantum system can be calculated. 33 Then, we can define the local observable property by splitting the total electronic density into the contributions for different subsystems. This partitioning is performed using the quantum theory of atoms in molecules (QTAIM) introduced by Richard Bader. 34 Among other topological properties, QTAIM enables to assign different regions of space to specific atoms by following the gradient of the electron density.
Therefore, we can decompose our total electron density into the sum of electron densities which belong to different molecules forming the complex system of interest. Using this approach, we can determine the local time-dependent density for each chromophore and compute its contribution to the photo-absorption spectrum by applying equations (5)- (8).
It is important to remark that, since we stay in the linear-response regime, the boost on the initial wave function does not produce a significant change in the electron density distribution. Instead just small fluctuations are observed. For this reason, partitioning of the ground state density is enough to define individual molecular space regions, and to study the time evolution of the electron density inside these domains.

Transition Dipole Moments and Transition Densities from P-TDDFT
In order to properly characterize electronic transitions it is worthy to obtain transition dipole moments, µ n0 = Ψ 0 |r|Ψ n , and transition densities, ρ T I0 (r) = Ψ I |n(r)|Ψ 0 , wherer is the position operator andn(r) the density operator. From equations (5)- (7) we can see that transition dipole moments are related with the imaginary parts of Fourier transforms of time-dependent dipole moments. Analogously, in P-TDDFT, transition densities are related with the imaginary parts of Fourier transforms of induced time-dependent densities. [35][36][37][38] Recently, Kummel et at. 30 derived an analytic expression to compute transition dipole moments and transition densities from P-TDDFT calculations, based on the coefficient expansion of the time-dependent wave function in terms of eigenfunctions of the many-body ground-state Hamiltonian.
In this work, we present a complementary derivation of both properties based only in LR-TDDFT arguments. Let us start with the key object in the linear-response TDDFT, which is the definition of the response density as where χ nn is the density-density response function: 14 The density operator here is defined in second quantization asn =â † jâ i , and the subindex H 0 stands for the zero-order Hamiltonian, i.e. with no external pertubation.
In the frequency space this gives the well-known "Lehmann-representation": From this expression we can see that the information regarding all transition densities for excitations from the ground state (Ψ 0 ) to any excited state (Ψ I ) is kept in the density-density response function, since the transition density is defined by ρ T I0 (r) = Ψ I |n(r)|Ψ 0 .
We can now make use of the "fluctuation-dissipation theorem", 15 that relates the density-density response function (χ nn ) with the corresponding dynamical structure factor (S nn ) at zero temperature: Since we assume T = 0 K , the dynamical structure factor fulfills S nn (ω) = 0 for ω < 0 and S nn (−ω) = 0 for ω > 0. Then, we do not consider stimulated emission and focus on the absorption spectrum (i.e. consider just ω > 0). By replacing equation (12) in equation (9), where the system is perturbed by an instantaneous electric field, and using equation (3), we obtain the following equation: Taking into account that the component of the transition moment of a given excitation can be computed from the corresponding transition density (i.e. µ n,ν = r ν ρ T n (r)dr) , one gets r · k Ψ 0 |n(r)|Ψ n dr = µ n,k · k, where k is the magnitude of the external electric field (equation (3)), and µ n,k corresponds to the projection of the transition dipole moment (0 → n, n excited state) over the polarization direction of the electric field.
Then, the imaginary part of the response density can be expressed in terms of all transition densities scaled by the projection of the corresponding transition moment on the direction of the external field, If we are interested in the excited states that are energetically far one from each other, the transition densities can be recovered from the response densities under an active perturbation, i.e. the one for which the transition is dipole allowed and µ n,k = 0,

Finite propagation times in P-TDDFT
As stated in equation (17), we first need to know the transition moment for a given excitation in order to recover the transition density.
To obtain the photo-absorption spectrum from finite propagation times in P-TDDFT, a damping function is commonly added in the Fourier transform of the polarizability tensor determined by equation (4) (see Supplementary Information). It induces an artificial decay of the excited population and allows to get a smooth spectrum removing artificial peaks. In particular, in our calculations we introduce a Gaussian damping function (e −η 2 t 2 ). We set the parameter η such that the damping function reaches the value of 1e −4 at the end of the propagation time τ : The damping parameter is responsible for the spectral broadening (see Supplementary Information) so that each single excitation is represented by a Gaussian function with half-width at half-maxima (σ hwhm ) Transition dipole moments can be obtained from a P-TDDFT calculation by diagonalization of the dynamic polarizability tensor given by equation (4). The frequency-dependent diagonal elements of the diagonalized polarizability tensor give the information about the probability of orthogonal excitations (i.e. transitions charactherized by orthogonal transition dipoles), and enables to distinguish between quasi-degenerate states.
The area under the trace of the diagonalized dynamic polarizability tensor is proportional to the square of the transition dipole moment (see equation (7)) and the direction of the transition dipole moment is determined by the corresponding eigenvectors.
In order to obtain the magnitudes of the transition dipole moments, we fit the spectra generated from the diagonalized frequency-dependent dynamic polarizability tensor by a sum of Gaussian functions. The Gaussian fit is obtained using the methodology described by Goshtasby et al. 39 Then the broadening of the Fourier Transform and the location of the excitation energies are found from the fitted parameters of the Gaussian half-width at half-maxima (σ n ) and the Gaussian center, respectively ( n ).
Knowing the transition dipole moments and broadening introduced due to the damping function, transition densities can be computed from P-TDDFT on finite times by taking the imaginary part of the Fourier transform of the time-dependent response density applying the same Gaussian damping. Then, integrating equation (17) around the peak at Ω n , we obtain the analytical expression for the transition density: where σ n is precisely the Gaussian fitted half-width at half-maxima around Ω n .

Exciton Coupling Evaluation
Both techniques, the local density analysis and calculation of transition densities from P-TDDFT, can be combined in order to study complex systems containing more than one photo-active molecule.
The exciton dynamics is governed by the exciton coupling determined by the Frenkel where |n is the excited state located on the site n, n are the site energies (or vertical excitation energies) and V nm corresponds to the exciton coupling between excited states on the donor site n, and the acceptor site m. The exciton coupling is determined by the interaction between the transition densities, 40 where g xc is the exchange-correlation kernel, and ω n/m is the difference energy between states n and m.
The Coulombic term inside the kernel function is usually attributed to the Förster resonant energy transfer (FRET). 41 This term stands for the interchange of the exciton population between two sites where the fluorescent energy of the donor is absorbed by the acceptor producing a resonant excitation. The second term in the kernel function corresponds to the exchange-correlation effect, and is related the Dexter electron transfer. 42 In this mechanism, the orbital proximity enables the simultaneous interchange of (i) an excited electron from the donor to an unoccupied orbital of the acceptor, and (ii) an electron from highest occupied molecular orbital (HOMO) of the acceptor to a hole on low lying single occupied molecular orbital (SOMO) of the excited donor ( Figure 1).
There are several ways to compute the exciton coupling between two excited states of different sites. The most commonly used ones are the ideal dipole approximation (IDA) 43 and the transition density cube (TDC) method. 44 The former assumes that both, acceptor and donor molecules, are far enough with no overlap of the electronic density and then the Coulomb interaction can be written in terms of the point dipole approximation as with the orientation factor In the TDC method, the transition densities are discretized into small volume units and then the exciton coupling is computed through the following sum over the volume units: The TDC method can therefore include both Coulomb and exchange-correlation contributions to the exciton coupling (equation (22)) by evaluating the Coulomb and exchangecorrelation potentials generated by the donor molecule in the region embracing the acceptor molecule. The exciton coupling obtained using the TDC method with the local transition densities calculated within our P-TDDFT formalism is hereafter denoted as V (T DCM ).

Procedure and Computational Details
Summarizing, local transition densities for a complex multichromophore system can be ob- OCTOPUS is a quantum chemistry/physics code specially designed to solve ground-state Figure 2: Scheme of the procedure used to get transition dipole moments and transition densities from P-TDDFT calculations for multichromophore systems.
and time-dependent (TD) DFT problems. [24][25][26] In OCTOPUS the basic quantities (Hamiltonian, potentials and single-particle wave functions) are considered on a discretized read-space grid.
Besides, OCTOPUS is highly parallelized in both grid points and Kohn-Sham states, which enables fast evaluation of non-linear propagation equations for the initial Slater determinant during real-time electron dynamics. These features make OCTOPUS the perfect code to implement our derivation of transition densities and exciton couplings from P-TDDFT.
In this work we use also use the ORCA package that allows us efficient computation of excited state properties by solving the Casida's equation using atom localized basis sets. 45,46 ORCA package uses several approaches to reduce the computational cost providing a good compromise against the accuracy of the results, e.g. the resolution of the identity (RI) approximation for the Coulomb potential. 47 Further computational details of the geometry optimization, ground-state calculations and TDDFT calculations using linear-response and time-propagation schemes are described in Supplementary Information. The calculated numerical data that support our study are available in "NOMAD repository" ( http://dx.doi.org/10.17172/NOMAD/2019.02.27-2).

Isolated Molecule
Let us exemplify the applicability of the formalism described above by computing the excitation properties for a single benzaldehyde molecule (Figure 3a). Figure 3b shows the ben-    Table 1: Excitation energies ( n , in eV), components of the transition dipole moment (µ x , µ y , µ z , in Å) and oscillator strength (f n ) for the lowest bright excited states of the isolated benzaldehyde molecule determined by Gaussian fitting of the first eigenvalue of the diagonalized polarizability tensor. In Table 1  are also shown in rows 2 and 4. We can see that both methods give the same results even though they rely on very different approaches and implementations, validating in this way the methodology for calculation of transition dipole moments from P-TDDFT.
Analogously, we can compare the respective transition densities obtained using both methods. Figure 4 shows the isocontours for the two bright states in Table 1  and (c) the average over the propagation runs for all three axes x,y and z. We can see that the propagation with the electric field along the y axis (k = (0, 1, 0)) induces a lot of noise in the calculation of the transition density due to the low excitation probability for this direction. Figure 4d and 4e show the transition densities for the second bright state (around 4.8 eV) obtained from LR-TDDFT and P-TDDFT, respectively. For the latter, we did a propagation run with the electric field along the y axis, which in this case corresponds to the maximum projection of the transition dipole moment (Table 1). Again, the perfect match between the transition densities computed for the first two bright excitations of the isolated benzaldehyde molecule using different methods, LR-TDDFT and P-TDDFT, validates our new methodology.

Dimers
Now we can try to extract the information on the exciton dynamics in a more complex system consisting of two molecules applying the local density analysis procedure to get the local transition densities. Once the transition densities for each chromophore are obtained, the exciton coupling is computed using equation (24).
Let us exemplify this method by considering the exciton coupling for the face-to-face benzaldehyde dimer as a function of its separation (Figure 5a). We place two benzalhehyde molecules at distances along the z axis ranging from 4 Å to 24 Å. For each of these dimers, we perform P-TDDFT calculations to get photo-absorption cross-sections (see Figure 5b).
The effect of the inter-molecular distance is immediately observed in the photo-absorption spectrum. We can see that as the molecules approach each other, a small blue shift develops on the higher peak (around 4.8 eV). There is also a transfer of the oscillator strength at low distances for the low-lying states (around 4.2 eV) and a small peak appears below 4 eV for the very close distance.
The contributions of each individual molecule to the total spectrum can be obtained performing the local density analysis. Figure 5c shows  (23)) and using the TDC method (equation (24)). Figure   5d shows the variation of the exciton couplings obtained using different methods as functions of the separation along the z axis. From comparison of the couplings V(LDT) and V(TDCM) we see that for benzaldehyde dimers, the IDA validity is limited to separations above 10 Å, in good agreement with the results shown by Hoffman et al. 49 This demonstrates the relevance of account of spatial distribution of the transition density in calculations of the exciton coupling.
The deviation between the exciton couplings computed using the transition dipole moments of the isolated molecules (V(IM)) [ = 4.78 eV, µ(3 A 1 ) = (-0.0052 Å, -0.7293 Å, 0.0940 Å)] and transition dipole moments computed from P-TDDFT (V(LDT)) ( Table 2) reveals that the electron excitation nature is slightly modified below 8 Å, i.e. the direction and strength of the transition dipole moment are changed due to the molecule-molecule interaction.   Table 3 shows the values of the computed exciton coupling for the different methods.

Columns 6 and 5 show the contributions of the Coulomb and exchange-correlation (XC)
terms, respectively, to the coupling V(TDCM) (equation (24)). We see that the Förster mechanism coming from the Coulomb interaction governs at the considered distances and it is not until very close distances (below 6 Å) that the Dexter mechanism related to the exchange-correlation effects becomes non-negligible. This fact is also attributed to the use of the semi-local GGA exchange-correlation functional such as PBE. 50,51 We expect that the use of long-range XC functionals 52-55 modifies this trend and the XC term should become more significant also at larger distances. Table 3: Exciton couplings V(IM), V(LDT), V(TDCM) (in cm -1 ) and contributions of the Coulomb (V h ) and exchange-correlation (V XC ) terms to V(TDCM) for dimers of benzaldehyde with different separations R (in Å) along the z axis.

Complex multichromophore system
Lets now consider a complex system with more than two photo-active molecules. We have created a cluster consisting of 14 randomly oriented benzaldehyde molecules. To ensure a plausible disposition, during the cluster generation, we avoid inter-atomic distances between atoms of different molecules smaller than the sum of the covalent radii plus 30% of the corresponding radii. Figure 6a shows the spacial distribution of the 14 molecules in the cluster with the distances between centers-of-mass of the molecules ranging between 5.17 Å and 19.86 Å. After the calculation of the ground state for the full system, the local density analysis is carried out in order to determine the charge density belonging to each molecule. The appropriateness of each local density domain is validated by its total charge. It is checked that the valence charge corresponds to 40 electrons for each molecule, with the maximum deviation found of 0.003 electron charge (which corresponds to the error of 0.00075%). Then we perform three P-TDDFT calculations with the electric field perturbation applied along each of the Cartesian axes. We propagate the system for the total propagation time of 40 fs.
For each local domain, the dynamic polarizability is calculated from the local induced dipole moment, and therefore the local contribution of each molecule to the global absorption spectrum is obtained (see Figure 6b). Besides, as described above, the transition dipole moments for a given peak are obtained by fitting the diagonalized dynamic polarizability tensor (see Table S2). Finally, we compute the transition density for a given excitation for each molecule of the cluster defined for the corresponding charge density domain using equation (20). The computed transition densities at the most intense peak below 5 eV are influenced by the presence of the embedding environment (see Figure S2). The effects due to spurious contamination of near located bright states can be evaluated by comparing the transition dipole moments obtained by the Gaussian fit and by applying the dipole operator over the transition density, analogously to the equation (5). Tables S1-S3 show this comparison, which corresponds to the mean deviation of 0.02 Å in the transition dipole moments, and the mean angle twist less than 6 • .
We compute now the exciton coupling for the entire system using the different approaches described above. A comparison of the values obtained for all pairs exciton coupling are specified in Tables S4-S5. Figure 7 shows the differences in the exciton couplings computed using TDCM (equation (24)) and LDT (equation (23)), i.e. the transition densities from equation (20) and the transition dipole moments from the fit and diagonalization of the dynamic polarizability tensor (equation (7)). The difference between both methods is not homogeneous and strongly depends on the specific pair of molecules. We can see from Figure 7b that the ideal dipole approximation is valid for the intermolecular distances (distances between the centers of mass) larger than 12 Å. Hence the excitation energy transfer between those molecules can be explained by the Förster resonant mechanism. However, at lower distances, this approximation fails as the V(TDCM) and V(LDT) values largely differ. This fact demonstrates the need for the method that enables accurate calculations of transition densities for individual subunits by treating all the system at the same level of theory.

CONCLUSIONS
In this work, we demonstrate that treating all the complex photoactive system at the ab initio level of theory allows accurate calculations of the exciton dynamic properties.
In this paper we derived an analytic method to obtain transition dipole moments and transition densities from time-propagation TDDFT calculation in the linear-response regime (P-TDDFT). We validated that the transition dipole moment can be obtained by Gaussian fitting of the absorption spectra and its direction is recovered from the eigenvectors of the dynamic polarizability diagonalization. We proved that the transition density for a given excitation can be accurately recovered from the Fourier transform of the time-dependent response density using the previous knowledge on the corresponding transition moment.
More interesting, we introduced the local density analysis to determine the transition properties for multichromophore systems. This procedure allows to accurately extract transition dipole moments and transition densities for individual molecules taking into account all effects due to their environment, since the entire system is modeled at the same level of theory. The calculation of these properties enables studies of exciton coupling in complex systems, a key parameter to understand the energy transfer processes taking place, for example, in photosynthetic light-harvesting complexes or the light-emitting layer of OLED devices.
This method also makes possible to disentangle the mechanisms of exciton transfer. The contribution of different interactions can easily be evaluated. In this work, we included the Coulomb and Exchange-Correlation effect (for density-dependent XC functionals), but other mechanisms such as an explicit terms due to the polarizable environment 56,57 can be considered as well.
It is important to mention that this method is implemented into the open-source OCTOPUS package, which provides a platform for performing P-TDDFT in the real-space basis. It is highly parallelized and is therefore suitable for efficient calculations of exciton couplings in large and complex systems.

Supporting Information Available
The following files are available free of charge.
The following file is available free of charge.
• Supplementary Information document includes: Section S1. Computational Details: Extended description the computational details used to obtain the ab initio results.
Section S2. Damping factor effect on the peak broadening: analytical derivation on the broadening of the absorption spectrum peaks due to the introduction of a Gaussian damping function on the short-time Fourier transform. Section S3 and Section S4, supplementary figures and tables that support main text results.