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 in using natural sun light for an energy transition toward 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−7] or to design new emitters[8,9] for light and flexible organic light emitting diodes (OLEDs)[10−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 addressing the eigenvalue problem either by solving the Schrödinger equation using a multiconfigurational ansatz[13] or the single-particle Kohn–Sham equations in the time-dependent formalism of the 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 us 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 electron–hole pair 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., subsystems TDDFT,[19] 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. On the other hand, Real-Time Propagation TDDFT (P-TDDFT) provides an efficient alternative to these methods. It has been demonstrated that P-TDDFT is an excellent platform for studies of such properties and processes such as molecular[24] and electron[25] dynamics, linear and nonlinear optics,[26] transport properties,[27] single and triplet excitations,[28,29] dynamical hyperpolarizabilites,[30] and exciton decay dynamics.[31] But probably, the most useful advantage that P-TDDFT offers is the possibility to obtain all frequency excitations at the same cost having converged just the occupied KS states for the ground state.[32,33] Also, propagation of the KS states can be highly parallelized enabling us to compute optical properties for up to several thousands of atoms.[34,35]

These features make P-TDDFT the suitable theoretical framework to study photoactive complex systems such as natural LHC and OLEDS. In such cases, it is important to take into account the environment effects, which is commonly done by adding polarizable force fields[36] or by other methods.[16]

It is well known that the absorption cross section can be computed from the propagation of the KS states in the linear-response regime. Also, the transition densities and plasmons can be qualitatively evaluated from the Fourier transform of the time-dependent induced density.[37−40] Recently Schelter et at.[41] 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 a quantum mechanical solution for the response functions after a boost excitation have the cardinal sine form, and they exemplified their method by performing real-time propagation TDDFT calculations. In that paper, however, they did not address 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 it 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[34] 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 us 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[42] and light-emitting layers of OLEDs.[43]

As a proof of concept of the method here presented, we illustrate its validity 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 us to decompose the absorption spectra into individual contributions.

In the next 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 momentwhere t0 is the initial time just before the application of the time-dependent external field, (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 polarizability 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λ(t) = r̂λ(t), respectively (where e refers to the electron charge), Θ(t – t′) is the Heaviside function, and |Ψ0⟩ is the wave function at the initial moment of time (t0).

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]where k corresponds to the magnitude vector of the electric field, which should be sufficiently small to ensure the linear-response regime.

Then, the polarizability tensor can be written in frequency space (ω) in terms of the Fourier transform of the time-dependent dipole moment as

Note that if k is sufficiently small so that the system is in the linear-response regime, the polarizability tensor does not depend on it.

In the time-dependent Kohn–Sham scheme, the time-dependent dipole moments can be easily obtained from the time-dependent electronic density (ρ(,t)))wherewhere Φ(,t) are the KS orbitals and the summation runs over the Nocc 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, aswhere the |Ψ⟩ is the wave function of the nth excited state, Ω is the corresponding excitation energy, and f 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 ().

In addition, from the imaginary part of the dynamic polarizability tensor, we can obtain the photoabsorption cross section, which is gives the absorption spectrum of the systemwhere c is the speed of light and 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).[34] This procedure is 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.[44] 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 Bader.[45] Among other topological properties, QTAIM enables us 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 photoabsorption spectrum by applying eqs –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 worth obtaining transition dipole moments, and transition densities, ρ() = ⟨Ψ|n̂()|Ψ0⟩, where is the position operator and n̂() the density operator. From eqs –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.[37−39,46]

Recently, Schelter et at.[41] 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 aswhere χ is the density–density response function[14]

The density operator here is defined in second quantization as n̂ = â†â, and the subindex H0 stands for the zero-order Hamiltonian, i.e., with no external perturbation.

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 (Ψ) is kept in the density–density response function, since the transition density is defined by ρ() = ⟨Ψ|n̂()|Ψ0⟩.

We can now make use of the “fluctuation–dissipation theorem”[15] that relates the density–density response function (χ) with the corresponding dynamical structure factor (S) at zero temperature

Since we assume T = 0 K, the dynamical structure factor fulfills S(ω) = 0 for ω < 0 and S(−ω) = 0 for ω > 0. Then, we do not consider stimulated emission and focus on the absorption spectrum (i.e., consider just ω > 0). By insterting eq into eq , where the system is perturbed by an instantaneous electric field, and using eq , we obtain the following equation

Taking into account that each component of the transition moment of a given excitation can be computed from the corresponding transition density (i.e., μ = ∫ rνρ()d), one getswhere k is the magnitude of the external electric field (eq ), and μ 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 μ ≠ 0

Finite Propagation Times in P-TDDFT

As stated in eq , we first need to know the transition moment for a given excitation in order to recover the transition density.

To obtain the photoabsorption spectrum from finite propagation times in P-TDDFT, a damping function is commonly added in the Fourier transform of the polarizability tensor determined by eq (Supporting Information). It induces an artificial decay of the excited population and allows us to get a smooth spectrum removing artificial peaks. In particular, in our calculations, we introduce a Gaussian damping function (e–η). 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 (Supporting 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 eq . The frequency-dependent diagonal elements of the diagonalized polarizability tensor give the information about the probability of orthogonal excitations (i.e., transitions characterized by orthogonal transition dipoles) and enables us to distinguish between quasidegenerate states.

The area under the trace of the diagonalized dynamic polarizability tensor is proportional to the square of the transition dipole moment (eq ), 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.[47] 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 (σ) and the Gaussian center, respectively (ϵ).

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 eq around the peak at Ω, we obtain the analytical expression for the transition densitywhere σ is precisely the Gaussian fitted half-width at half-maxima around Ω.

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 photoactive molecule.

The exciton dynamics is governed by the exciton coupling determined by the Frenkel Exciton Hamiltonianwhere |n⟩ is the excited state located on the site n, ϵ are the site energies (or vertical excitation energies), and V 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[48,49]where g is the exchange-correlation kernel, and ω is the resonant excitation energy for states n and m for symmetric systems; otherwise, the average of the acceptor (ω) and donor (ω) excitation frequencies is used. Note that this expression was derived considering the interaction between chromophores as a perturbation. It is not accurate when the chromophores are so close that their electronic wave functions strongly interact, mostly due to neglect of the charge-transfer effects.

The Coulombic term inside the kernel function is usually attributed to the Förster resonant energy transfer (FRET).[50] 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.[51] 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 ).

Figure 1

Schematic representation of the excitation energy transfer for (a) Förster and (b) Dexter mechanisms.

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)[52] and the transition density cube (TDC) method.[53] 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 aswith the orientation factor κ = · – 3(·)(·), where are the unit vectors describing directions of transition dipole moments for the donor (D) and acceptor (A), is the distance vector connecting both centers of mass, and is the unit vector along the same direction.

This approximation does not take into account the delocalization of the transition densities over the molecules. In other words, different regions of the electronic density of the acceptor molecule feel different potentials generated by the excited density of the donor molecule due to the spacial arrangement.

The use of the IDA approximation requires the knowledge of the transition dipole moments both for donor and acceptor molecules. In Section , we compare exciton couplings computed using different transition dipole moments obtained (i) from consideration of isolated molecules, named here as V(IM) (where IM stands for isolated molecule) and (ii) from the proposed local density analysis, V(LDT) (where LDT stands for local density treatment, not to be confused with the LDA exchange-correlation functional).

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 (eq ) by evaluating the Coulomb and exchange-correlation 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(TDCM).

Procedure and Computational Details

Summarizing, local transition densities for a complex multichromophore system can be obtained following the computational steps shown in Figure . After preforming a ground state (GS) calculation, the local domains are determined to get the densities belonging to each molecule. We have checked that since we stay in a linear-response regime during the time-dependent (td) propagation the electronic charge inside each domain does not vary significantly (less than 0.005%, Figure S1); hence, there is no need to update the grid point assignment to different spatial domains during the time-propagation calculations. Three P-TDDFT calculations are executed by perturbing the initial state with an electric field polarized in each of the Cartesian directions, recording the dipole moments and saving the time-dependent density every 20 steps of the propagation. Then, the dynamic polarizability tensor for each domain is obtained and diagonalized to get the local contribution to the absorption spectrum as well as the transition dipole moments. The frequency-dependent response densities are obtained by the Fourier transform of the stored time-dependent electron density. Finally the transition densities for each molecule and excitation can be reconstructed, and the exciton coupling can be computed using the TDC method.

Figure 2

Scheme of the procedure used to get transition dipole moments and transition densities from P-TDDFT calculations for multichromophore systems.

OCTOPUS is a quantum chemistry/physics code specially designed to solve ground-state and time-dependent (TD) DFT problems.[54−56] In OCTOPUS, the basic quantities (Hamiltonian, potentials and single-particle wave functions) are considered on a discretized real-space grid. Also, OCTOPUS is highly parallelized in both grid points and Kohn–Sham states, which enables fast evaluation of nonlinear 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 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.[57,58] The 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.[59]

In order to be sure that both types of calculations (using ORCA and OCTOPUS) are consistent, we compared the ground-state eigenvalues (Table S1) and excited-state energies (Table ). It can be seen that in all cases the energy difference between the methods is smaller than 0.05 eV.

Table 1

Excitation Energies (ϵ, in eV), Components of the Transition Dipole Moment (μ, μ, μ, in Å), and Oscillator Strength (f) for the Lowest Bright Excited States of the Isolated Benzaldehyde Molecule Determined by Gaussian fitting of the First Eigenvalue of the Diagonalized Polarizability Tensor

Method	ϵn	μx	μy	μz	fn
P-TDDFT	4.21	0.1972	0.0055	0.0439	0.0150
LR-TDDFT	4.23	0.1982	0.0042	0.0440	0.0152
P-TDDFT	4.78	–0.0091	–0.7153	0.0914	0.2175
LR-TDDFT	4.80	–0.0085	–0.7178	0.0915	0.2200

Further computational details of the geometry optimization, ground-state calculations, and TDDFT calculations using linear-response and time-propagation schemes are described in the Supporting Information. The calculated numerical data that support our study are available in the “NOMAD repository” (http://dx.doi.org/10.17172/NOMAD/2019.02.27-2).

Applications and Discussion

Isolated Molecule

Let us exemplify the applicability of the formalism described above by computing the excitation properties for a single benzaldehyde molecule (Figure a). Figure b shows the benzaldehyde UV–visible absorption spectrum computed solving the Casida’s equation within linear-response TDDFT.[18] It presents two major optically allowed excitations at 4.2 and 4.8 eV.

Figure 3

(a) Benzaldehyde molecule orientation. Brown, white, and red spheres represent carbon, hydrogen, and oxigen atoms, respectively. (b) Absorption spectrum of the isolated benzaldehyde molecule computed using P-TDDFT (purple line) and LR-TDDFT (yellow impulses). The decomposition of P-TDDFT results in the contributions from different eigenvalues of the diagonalized polarizability tensor indicated using green, blue, and orange lines.

Figure b also shows how the absorption spectra obtained from the trace of the dynamic polarizability tensor can be decomposed by performing a respective diagonalization. In the frequency regions where just single-degenerate bright states are present [4–5.5 eV], we can see that the first eigenvalue of this diagonalization matches the total contribution of the spectrum. However, in the case when almost degenerate states are present [5.5–6 eV], we observe that although the trace of the dynamic polarizability shows just one broad peak the first eigenvalue splits this band, in good agreement with the LR-TDDFT results.

In Table , we give the values for the excitation energy and transition dipole moments obtained by the Gaussian fitting procedure of the dynamical polarizability tensor in the framework of P-TDDFT as implemented in the OCTOPUS package.[35,55] The values obtained from the solution of the Casida’s equation as implemented in the ORCA package[57] are 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 shows the isocontours for the two bright states in Table . The first transition density (Figure a) is obtained using the eigenvectors of Casida’s equation. The transition densities in Figure b and c are obtained from the time-propagation of the Kohn–Sham orbitals after applying an initial perturbation keeping the system in the linear regime (P-TDDFT). From eq , it can be seen that we can obtain the transition density from any propagation where the projection of the transition dipole to the polarization direction of the electric field does not vanish. For the first excitation of the benzaldehyde molecule, the transition dipole moment is almost aligned along the x axis (Table and Figure a). Although, in principle, the components of the transition dipole moment for the other directions are not exactly zero, they take small values, and their use for the calculations of the transition density can lead to errors. To demonstrate this, we compared transition densities for the propagation run with the electric field along the x axis ( = (1,0,0)) (Figure b) and the average over the propagation runs for all three axes x, y, and z (Figure c). We can see that the propagation with the electric field along the y axis ( = (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 d and e 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 ).

Figure 4

(a–c)Transition densities (isovalue of 0.005 Å–3) for the first bright state of the isolated benzaldehyde molecule located at 4.2 eV obtained from (a) LR-TDDFT and (b) P-TDDFT with the electric field along the x axis, which corresponds to the maximum component of the transition dipole moment and (c) average over P-TDDFT runs with the electric field along x, y, and z axes. (d, e) Transition densities (isovalue of 0.0075 Å–3) for the second bright state located at 4.8 eV obtained from (d) LR-TDDFT and (e) P-TDDFT with the electric field along the y axis, which corresponds to the maximum component of the transition dipole moment.

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 eq 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 a). 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 photoabsorption cross sections (Figure b). The effect of the intermolecular distance is immediately observed in the photoabsorption 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.

Figure 5

(a) Geometry of benzaldehyde dimers. Red, yellow, and blue vectors represent the x, y, z Cartesian axes, respectively. We call “Bnzld Left” the benzaldehyde molecule in gold with a smaller z coordinate (left) and “Bnzld Right” the purple molecule (right). (b) Comparison of the simulated absorption spectra for benzaldehyde dimers for different separations of the molecules along the z axis. (c) Spectrum decomposition using the local density analysis for the dimer with the separation along the z axis of 4 Å. (d) Exciton coupling evaluated for the strongest excitation around 4.8 eV for different benozaldehyde dimers.

The contributions of each individual molecule to the total spectrum can be obtained performing the local density analysis. Figure c shows the contributions of each molecule for the dimer with the separation along the z axis of 4 Å. Due to the small torsion with respect to the xy plane and z axis, similar atoms of the molecules are not perfectly faced in and have a different environment. This results in a slightly different excitation energy for each molecule, shown in the figure by dashed lines. After the reconstruction of the individual transition densities with account of the effect of the other molecule, we compute the excitation coupling in the IDA approximation (eq ) and using the TDC method (eq 24). Figure d 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.[60] This demonstrates the relevance of counting the 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′A1) = (−0.0052 Å, −0.7293 Å, 0.0940 Å)] and transition dipole moments computed from P-TDDFT (V(LDT)) (Table ) 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 2

Comparison of Excitation Energies (ϵ, in eV), Components of the Transition Dipole Moment (μ, μ, μ, in Å), and Oscillator Strength (f) for Molecules in Benzaldehyde Dimers with Different Separations R (in Å) along the z Axis

R	ϵ	μx	μy	μz	fosc
4	4.85	–0.0128	–0.5974	0.0722	0.1537
	4.82	–0.0102	–0.5994	0.0963	0.1556
6	4.81	0.0088	0.6651	–0.0897	0.1896
	4.81	0.0071	0.6678	–0.0902	0.1911
8	4.80	–0.0063	–0.6930	0.0924	0.2053
	4.80	–0.0059	–0.6934	0.0919	0.2055
10	4.79	0.0051	0.7048	–0.0928	0.2118
	4.79	–0.0052	–0.7050	0.0924	0.2119
12	4.79	–0.0049	–0.7105	0.0929	0.2152
	4.79	–0.0049	–0.7103	0.0927	0.2151
16	4.79	–0.0047	–0.7153	0.0930	0.2181
	4.79	–0.0047	–0.7152	0.0929	0.2180
20	4.79	–0.0047	–0.7276	0.0944	0.2256
	4.79	–0.0047	–0.7277	0.0943	0.2256
24	4.78	–0.0048	–0.7284	0.0943	0.2256
	4.78	–0.0048	–0.7286	0.0944	0.2257

Table 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) (eq 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 semilocal GGA exchange-correlation functional such as PBE.[61,62] We expect that the use of long-range XC functionals[63−66] 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 Coulomb (Vh) and Exchange-Correlation (VXC) Terms to V(TDCM) for Dimers of Benzaldehyde with Different Separations R (in Å) along the z Axis

R	V(IM)	V(LDT)	V(TDCM)	Vh	VXC
4.0	933.313	618.470	329.559	328.114	1.445
6.0	276.537	231.079	168.345	168.3662	0.0210
8.0	116.664	105.256	87.874	87.874	0.000
10.0	59.732	55.733	49.392	–49.392	0.000
12.0	34.567	32.746	30.170	30.170	0.000
16.0	14.583	14.034	13.353	13.353	0.000
20.0	7.467	7.420	7.640	7.640	0.000
24.0	4.321	4.2745	4.470	4.470	0.000

Complex Multichromophore System

Let us now consider a complex system with more than two photoactive molecules. We have created a cluster consisting of 14 randomly oriented benzaldehyde molecules. To ensure a plausible disposition, during the cluster generation, we avoid interatomic distances between atoms of different molecules smaller than the sum of the covalent radii plus 30% of the corresponding radii. Figure a 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 Å.

Figure 6

(a) Geometrical distribution of 14 benzaldehyde molecules in the cluster. (b) Spectrum decomposition for the 14 benzaldehyde cluster.

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 (Figure b). Also, as described above, the transition dipole moments for a each peak are obtained by fitting the diagonalized dynamic polarizability tensor (Table S3). Finally, we compute the transition density for a given excitation for each molecule of the cluster defined for the corresponding charge density domain using eq . The computed transition densities at the most intense peak below 5 eV are influenced by the presence of the embedding environment (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 eq . Tables S2–S4 show this comparison, which corresponds to a mean deviation of 0.02 Å in the transition dipole moments and a 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 of exciton coupling are specified in Tables S5 and S6. Figure shows the differences in the exciton couplings computed using TDCM (eq 24) and LDT (eq ), i.e., the transition densities from eq and the transition dipole moments from the fit and diagonalization of the dynamic polarizability tensor (eq ). The difference between both methods is not homogeneous and strongly depends on the specific pair of molecules.

Figure 7

(a) Correlation between the values computed using the transition density cube method (V(TDCM)) versus the values computed from the ideal dipole approximation using the transition dipole moments found for local domains from the fit of the diagonalized dynamic polarizability tensor (V(LDT)). (b) Difference between V(LDT) and V(TDCM) as a function of the center-of-mass distance.

We can see from Figure b 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 systems at the ab initio level of theory allows for accurate calculations of the exciton dynamic properties. Although our approach suffers from intrinsic limitations of some exchange-correlation functionals in the TDDFT method such as the proper description of charge-transfer states, it can provide a more accurate alternative to other approaches for weakly interacting large systems by including explicitly the quantum mechanical environment effects.

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 of the corresponding transition moment.

More interesting, we introduced the local density analysis to determine the transition properties for multichromophore systems. This procedure allows us 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 for disentangling 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 explicit terms due to the polarizable environment[48,67] 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.