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

TL;DR

This paper introduces a real-time propagation TDDFT method combined with a local density decomposition technique to analyze exciton couplings and transition properties in large, multi-chromophore systems, enabling detailed exciton dynamics studies.

Contribution

The authors extend their local density decomposition method to compute exciton dynamic properties using Real-Time Propagation TDDFT, providing a new analytical expression for transition density applicable to complex systems.

Findings

Validated method for transition dipole moments and densities

Analyzed exciton coupling in benzaldehyde dimers and clusters

Mapped exciton interactions in a 14-molecule benzaldehyde system

Abstract

Photo-active systems are characterized by their capacity of absorbing light energy and transforming 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 well-known Casida's equation for single molecules. However, this methodology is not useful in practice when dealing with multichromophore systems. In this work, we extend our local density decomposition method that enables 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 $\pi$-stacked dimers and finally map the exciton coupling for a 14 benzaldehyde cluster.

Tables

Table 1. Table 1 : Excitation energies ( ϵ n subscript italic-ϵ 𝑛 \epsilon_{n} , in eV), components of the transition dipole moment ( μ x , μ y , μ z subscript 𝜇 𝑥 subscript 𝜇 𝑦 subscript 𝜇 𝑧 \mu_{x},\mu_{y},\mu_{z} , in Å) and oscillator strength ( f n subscript 𝑓 𝑛 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.

Method	${ϵ}_n$	${\mu }_x$	${\mu }_y$	${\mu }_z$	$f_n$
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

Table 2. Table 2 : Comparison of excitation energies ( ϵ italic-ϵ \epsilon , in eV), components of the transition dipole moment ( μ x , μ y , μ z subscript 𝜇 𝑥 subscript 𝜇 𝑦 subscript 𝜇 𝑧 \mu_{x},\mu_{y},\mu_{z} , in Å) and oscillator strength ( f o ​ s ​ c subscript 𝑓 𝑜 𝑠 𝑐 f_{osc} ) for the molecules in the benzaldehyde dimers with different separations R 𝑅 R (in Å) along the z axis.

$R$	$ϵ$	${\mu }_x$	${\mu }_y$	${\mu }_z$	$f_{osc}$
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 3. 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 𝑅 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

Equations

$${\mu }_{\nu }(t)={\mu }_{\nu }(t_0)+{\int }_{-\infty }^{+\infty }{\alpha }_{\nu \lambda }(t-t^{\prime }){\mathcal{E}}_{\lambda }(t^{\prime })d{t}^{\prime },$$

$${\alpha }_{\nu \lambda }(t-t^{\prime })=-i\Theta (t-t^{\prime })⟨{\Psi }_0|[{\hat{\mu }}_{\nu }(t-t^{\prime }),v_{\lambda }^{ext}]|{\Psi }_0⟩,$$

$$v^{ext}=-e\text{r}\cdot \text{k}\delta (t_0),$$

$${\alpha }_{\nu \lambda }(\omega )=\frac{1}{k_{\lambda }}{\int }_0^{\infty }dt[{\mu }_{\nu }(t)-{\mu }_{\nu }(t_0)]e^{-i\omega t}.$$

$$\mu (t)=\int \text{r}\cdot \rho (\mathbf{r},t)d\text{r},$$

$$\rho (\mathbf{r},t)=\sum _{i=0}^{N_{occ}}|{\Phi }_i(\mathbf{r},t){|}^2,$$

$$\frac{1}{3}\text{Tr}(\Im m\mathbf{\alpha }(\omega ))=\frac{\pi }{3}\sum _{n=1}^{\infty }\sum _{\nu =1}^3|⟨{\Psi }_0|{\hat{r}}_{\nu }|{\Psi }_n⟩{|}^2\delta (\omega -{\Omega }_n)=\sum _{n=1}^{\infty }\frac{\pi }{2{\Omega }_n}{f}_n\delta (\omega -{\Omega }_n),$$

$$\sigma (\omega )=\frac{4\pi \omega }{3c}Tr[\Im m(\mathbf{\alpha }(\omega )]=\frac{2{\pi }^2}{c}S(\omega )$$

$$\delta \rho (\mathbf{r},t)={\int }_{t_0}^{t}d{t}^{\prime }\int d{\mathbf{r}}^{\prime }{\chi }_{n{n}^{\prime }}(\mathbf{r}t,{\mathbf{r}}^{\prime }{t}^{\prime })v^{ext}({\mathbf{r}}^{\prime },t^{\prime }),$$

$${\chi }_{n{n}^{\prime }}(\mathbf{r}t,{\mathbf{r}}^{\prime }{t}^{\prime })=-i\theta (t-t^{\prime })⟨{\Psi }_0|[{\hat{n}}_{H_0}(\mathbf{r},t),{\hat{n}}_{H_0}({\mathbf{r}}^{\prime },t^{\prime })]|{\Psi }_0⟩.$$

$$\chi_{nn^{\prime}}(\boldsymbol{r},\boldsymbol{r}^{\prime},\omega)=\sum_{I\neq 0}\Bigg{\{}\frac{\langle\Psi_{0}|\hat{n}(\boldsymbol{r})|\Psi_{I}\rangle\langle\Psi_{I}|\hat{n}(\boldsymbol{r}^{\prime})|\Psi_{0}\rangle}{\omega-\Omega_{I}+i\eta}-\frac{\langle\Psi_{0}|\hat{n}(\boldsymbol{r}^{\prime})|\Psi_{I}\rangle\langle\Psi_{I}|\hat{n}(\boldsymbol{r})|\Psi_{0}\rangle}{\omega+\Omega_{I}+i\eta}\Bigg{\}}.$$

$$\Im m\chi_{nn}(\boldsymbol{r},\boldsymbol{r}^{\prime},\omega)=-\pi\big{(}S_{nn}(\boldsymbol{r},\boldsymbol{r}^{\prime},\omega)-S_{nn}(\boldsymbol{r}^{\prime},\boldsymbol{r},-\omega)\big{)}$$

$$S_{nn}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )=n_0(\mathbf{r})n_0({\mathbf{r}}^{\prime })\delta (\omega )+\sum _{n=1}^{\infty }⟨{\Psi }_0|\hat{n}(\mathbf{r})|{\Psi }_n⟩⟨{\Psi }_0|\hat{n}({\mathbf{r}}^{\prime })|{\Psi }_n⟩\delta (\omega -{\Omega }_n).$$

$$\Im m\delta \rho (\mathbf{r},\omega )=\int d{\mathbf{r}}^{\prime }\pi S_{nn}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )e{\mathbf{r}}^{\prime }\cdot \mathbf{k}\delta (\omega -{\Omega }_n).$$

$$\int \mathbf{r}\cdot \mathbf{k}⟨{\Psi }_0|\hat{n}(\mathbf{r})|{\Psi }_n⟩d\mathbf{r}={\mu }_{n,k}\cdot k,$$

$$\Im m\delta \rho (\mathbf{r},\omega )=\pi k\sum _{n=1}^{\infty }{\mu }_{n,k}⟨{\Psi }_0|\hat{n}(\mathbf{r})|{\Psi }_n⟩\delta (\omega -{\Omega }_n).$$

$$⟨{\Psi }_0|\hat{n}(\mathbf{r})|{\Psi }_n⟩\delta (\omega -{\Omega }_n)=\frac{\Im m\delta \rho (\mathbf{r},\omega )}{\pi k{\mu }_{n,k}}.$$

$$D(t)=e^{-{\eta }^2{t}^2},$$

$$\eta =\frac{\sqrt[]{-ln(1e^{-4})}}{\tau }.$$

$${\sigma }_{hwhm}=2\sqrt[]{ln2}\cdot \eta ,$$

$${\rho }_{n0}^T(\mathbf{r})=⟨{\Psi }_0|\hat{n}(\mathbf{r})|{\Psi }_n⟩=\frac{{\sigma }_n}{{\mu }_{n,k}k\sqrt{\pi \ln 2}}\Im m\delta \rho (\mathbf{r},{\Omega }_n),$$

$$\hat{H}=\sum _n{ϵ}_n|n⟩⟨n|+\sum _{n,m}{V}_{nm}|m⟩⟨n|,$$

$$V_{nm}=\int d\text{r}\int d{\text{r}}^{\prime }{\rho }_m^{T\ast }({\text{r}}^{\prime })\left(\frac{1}{|{\text{r}}^{\prime }-\text{r}|}+g_{xc}({\text{r}}^{\prime },\text{r})\right){\rho }_n^T(\text{r})-{\omega }_{n/m}\int d\text{r}{\rho }_m^{T\ast }(\text{r}){\rho }_n^T(\text{r}),$$

$$V_{IDA}^{FRET}={\kappa }^{DA}\frac{|{\mathbf{\mu }}^A||{\mathbf{\mu }}^D|}{|\text{R}{|}^3}$$

$$M_X^{n0}(x,y,z)=V_{\delta }{\int }_z^{z+\delta z}{\int }_y^{y+\delta y}{\int }_x^{x+\delta x}{\Psi }_0^X{\Psi }_n^{X\ast }\text{d}x\text{d}y\text{d}z$$

$$V_s=\int d\text{r}\int d{\text{r}}^{\prime }{M}_A^{n0}({\text{r}}_A^{\prime })\left(\frac{1}{|{\text{r}}^{\prime }-\text{r}|}+g_{xc}({\text{r}}^{\prime },\text{r})\right)M_D^{0n}({\text{r}}_B)$$

Full text

University of the Basque Country] Nano-Bio Spectroscopy Group and ETSF Scientific Development Centre, Department of Materials Physics, University of the Basque Country, CFM CSIC-UPV/EHU-MPC and DIPC, Tolosa Hiribidea 72, E-20018 Donostia-San Sebastián \alsoaffiliation[Max Planck Institute for the Structure and Dynamics of Matter] Theory Department, Max Planck Institute for the Structure and Dynamics of Matter and Center for Free-Electron Laser Science, Luruper Chaussee 149, 22761 Hamburg, Germany

University of the Basque Country] Nano-Bio Spectroscopy Group and ETSF Scientific Development Centre, Department of Materials Physics, University of the Basque Country, CFM CSIC-UPV/EHU-MPC and DIPC, Tolosa Hiribidea 72, E-20018 Donostia-San Sebastián

\abbreviationsIR,NMR,UV

{tocentry}

1 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 addressing the eigenvalue problem either by solving the Schrödinger equation using a multiconfigurational ansatz13 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 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+TB23), 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 as molecular24 and electron25 dynamics, linear and non-linear 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 Besides, propagation of the KS states can be highly parallelized enabling to compute optical properties for up to several thousands of atoms. 34, 35.

These features make P-TDDFT the suitable theoretical framework to study photo-active 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 fields36 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. Besides, the transition densities and plasmons can be qualitatively evaluated from the Fourier transform of the time-dependent induced density.37, 38, 39, 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 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 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 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 complexes42 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.

2 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 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.

2.1 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:

[TABLE]

where $t_{0}$ is the initial time just before the application of the time dependent external field, $\mathcal{E}_{\lambda}(t)$, and the dynamic polarizability element, $\alpha_{\nu\lambda}(t-t^{\prime})$, is the retarded response that relates the $\lambda$ component of the external electric field to the change of the $\nu$ 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:

[TABLE]

where the dipole moment operator and the external potential are described as \hat{\mu}_{\nu}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}=e\hat{r}_{\nu} and \delta v_{\lambda}^{ext}(t)\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}=\hat{r}_{\lambda}\mathcal{E}(t), respectively (where $e$ refers to the electron charge), $\Theta(t-t^{\prime})$ is the Heaviside function and $|\Psi_{0}\rangle$ is the the wave function at the initial moment of time ($t_{0}$).

The polarizability tensor $\boldsymbol{\alpha(\omega)}$ 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:

[TABLE]

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 ($\omega$) in terms of the Fourier Transform of the time-dependent dipole moment as

[TABLE]

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 ($\rho(\boldsymbol{r},t))$):

[TABLE]

where

[TABLE]

where $\Phi_{i}(\boldsymbol{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

[TABLE]

where the $|\Psi_{n}\rangle$ is the wave function of the $n$-th excited state, $\Omega_{n}$ is the corresponding excitation energy and $f_{n}$ is the oscillator strength which gives the transition probability for that particular excitation. This expression relates directly the dynamic polarization tensor ($\boldsymbol{\alpha}(\omega)$) and the transition dipole moment ($\vec{\mu}_{n0}=\langle\Psi_{0}|\hat{\boldsymbol{r}}|\Psi_{n}\rangle$).

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:

[TABLE]

where $c$ is the speed of light and $S(\omega)=\sum_{n=1}^{\infty}f_{n}\delta(\omega-\Omega_{n})$ is the dipole spectral function.

2.2 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 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 Richard Bader.45 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 eqs 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.

2.3 Transition Dipole Moments and Transition Densities from P-TDDFT

In order to properly characterize electronic transitions it is worthy to obtain transition dipole moments, $\vec{\mu}_{n0}=\langle\Psi_{0}|\hat{\boldsymbol{r}}|\Psi_{n}\rangle$, and transition densities, $\rho^{T}_{I0}(\boldsymbol{r})=\langle\Psi_{I}|\hat{n}(\boldsymbol{r})|\Psi_{0}\rangle$, where $\hat{\boldsymbol{r}}$ is the position operator and $\hat{n}(\boldsymbol{r})$ the density operator. From eqs 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. 46, 37, 38, 39.

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 as

[TABLE]

where $\chi_{nn^{\prime}}$ is the density-density response function: 14

[TABLE]

The density operator here is defined in second quantization as $\hat{n}=\hat{a}_{j}^{\dagger}\hat{a}_{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":

[TABLE]

From this expression we can see that the information regarding all transition densities for excitations from the ground state ($\Psi_{0}$) to any excited state ($\Psi_{I}$) is kept in the density-density response function, since the transition density is defined by $\rho^{T}_{I0}(\boldsymbol{r})=\langle\Psi_{I}|\hat{n}(\boldsymbol{r})|\Psi_{0}\rangle$.

We can now make use of the "fluctuation-dissipation theorem" 15, that relates the density-density response function ($\chi_{nn}$) with the corresponding dynamical structure factor ($S_{nn}$) at zero temperature:

[TABLE]

[TABLE]

Since we assume $T=0$ K , the dynamical structure factor fulfills $S_{nn}(\omega)=0$ for $\omega<0$ and $S_{nn}(-\omega)=0$ for $\omega>0$. Then, we do not consider stimulated emission and focus on the absorption spectrum (i.e. consider just $\omega>0$). By insterting eq 12 in eq 9, where the system is perturbed by an instantaneous electric field, and using eq 3, we obtain the following equation:

[TABLE]

Taking into account that the component of the transition moment of a given excitation can be computed from the corresponding transition density (i.e. $\mu_{n,\nu}=\int r_{\nu}\rho_{n0}^{T}(\boldsymbol{r})d\boldsymbol{r}$) , one gets

[TABLE]

where $k$ is the magnitude of the external electric field (eq 3), and $\mu_{n,k}$ corresponds to the projection of the transition dipole moment ($0\rightarrow 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,

[TABLE]

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 $\mu_{n,k}\neq 0$,

[TABLE]

2.3.1 Finite propagation times in P-TDDFT

As stated in eq 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 eq 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^{-\eta^{2}t^{2}}$). We set the parameter $\eta$ such that the damping function reaches the value of $1e^{-4}$ at the end of the propagation time $\tau$:

[TABLE]

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 ($\sigma_{hwhm}$)

[TABLE]

Transition dipole moments can be obtained from a P-TDDFT calculation by diagonalization of the dynamic polarizability tensor given by eq 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 eq 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. 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 ($\sigma_{n}$) and the Gaussian center, respectively ($\epsilon_{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 eq 17 around the peak at $\Omega_{n}$, we obtain the analytical expression for the transition density:

[TABLE]

where $\sigma_{n}$ is precisely the Gaussian fitted half-width at half-maxima around $\Omega_{n}$.

2.4 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 Exciton Hamiltonian:

[TABLE]

where $|n\rangle$ is the excited state located on the site $n$, $\epsilon_{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, 48, 49

[TABLE]

where $g_{xc}$ is the exchange-correlation kernel, and $\omega_{n/m}$ is the resonant excitation energy for states $n$ and $m$ for symmetric systems, otherwise the average of the acceptor ($\omega_{n}$) and donor ($\omega_{m}$) 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 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)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 as

[TABLE]

with the orientation factor $\kappa^{DA}=\boldsymbol{u}^{D}\cdot\boldsymbol{u}^{A}-3(\boldsymbol{u}^{D}\cdot\boldsymbol{u}^{R})(\boldsymbol{u}^{A}\cdot\boldsymbol{u}^{R})$, where $\boldsymbol{u}^{D,A}$ are the unit vectors describing directions of transition dipole moments for the donor (D) and acceptor (A), $\boldsymbol{R}$ is the distance vector connecting both centers of mass, and $\boldsymbol{u}^{R}$ 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 potential 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 the Section 3 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 confuse 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:

[TABLE]

The TDC method can therefore include both Coulomb and exchange-correlation contributions to the exciton coupling (eq 22) 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)$.

2.5 Procedure and Computational Details

Summarizing, local transition densities for a complex multichromophore system can be obtained the following computational steps shown in Figure 2. 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%, see Figure SLABEL:SIfig:ld_gsvstd) ), hence there is no need to update the grid point assignment to different spatial domains during the time-propagation calculations. Three P-TDDFT calculation 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.

OCTOPUS is a quantum chemistry/physics code specially designed to solve ground-state and time-dependent (TD) DFT problems.54, 55, 56 In OCTOPUS the basic quantities (Hamiltonian, potentials and single-particle wave functions) are considered on a discretized real-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 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 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 the calculations (using ORCA and OCTOPUS) are consistent, we compared the ground-state eigenvalues (see Table SLABEL:SItbl:ks_oct_vs_orca) and excited-state energies (see Table 1). It can be seen that in all cases the energy difference between the methods is smaller than 0.05 eV.

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).

3 Applications and Discussion

3.1 Isolated Molecule

Let us exemplify the applicability of the formalism described above by computing the excitation properties for a single benzaldehyde molecule (Figure 3(a)). Figure 3(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 eV and 4.8 eV.

Figure 3(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 to the total contribution of the spectrum. Whereas 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 1 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. 55, 35 The values obtained from the solution of the Casida’s equation as implemented in the ORCA package57 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. The first transition density (Figure 4(a)) is obtained using the eigenvectors of the Casida’s equation. The transition densities in Figures 4(b) and 4(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 20 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 1 and Figure 3(a)). Though 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 ($\boldsymbol{k}=(1,0,0)$) (Figure 4(b)) and the average over the propagation runs for all three axes x,y and z (Figure 4(c)). We can see that the propagation with the electric field along the y axis ($\boldsymbol{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 4(d) and 4(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 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.

3.2 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 5(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 photo-absorption cross-sections (see Figure 5(b)). 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 5(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 23) and using the TDC method (eq 24). Figure 5(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)) [$\epsilon$ = 4.78 eV, $\mu(3^{\prime}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) (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 semi-local GGA exchange-correlation functional such as PBE.61, 62 We expect that the use of long-range XC functionals63, 64, 65, 66 modifies this trend and the XC term should become more significant also at larger distances.

3.3 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 6(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 Å.

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 6(b)). Besides, as described above, the transition dipole moments for a each peak are obtained by fitting the diagonalized dynamic polarizability tensor (see Table SLABEL:SItbl:cluster_mu_fit). 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 20. The computed transition densities at the most intense peak below 5 eV are influenced by the presence of the embedding environment (see Figure SLABEL:SIfig:cluster_tp). 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 eq 5. Tables SLABEL:SItbl:cluster_mu_ida-SLABEL:SItbl:cluster_mu_ftp 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 exciton coupling are specified in Tables SLABEL:SItbl:exc_coupling_lfret-SLABEL:SItbl:exc_coupling_ida. Figure 7 shows the differences in the exciton couplings computed using TDCM (eq 24) and LDT (eq 23), i.e. the transition densities from eq 20 and the transition dipole moments from the fit and diagonalization of the dynamic polarizability tensor (eq 7). The difference between both methods is not homogeneous and strongly depends on the specific pair of molecules.

We can see from Figure 7(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.

4 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. 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 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 environment48, 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.

{acknowledgement}

JJS and IL thank the European Research Council (ERC-2010-AdG-267374), Spanish grant (FIS2016-79464-P), and Grupos Consolidados (IT578-13) and EU-H2020 project “MOSTOPHOS” (n. 646259) for the financial support. JJS gratefully acknowledges the Spanish grant IJCI-2014-22204, and H2020-EINFRA-5-2015 project “NOMAD” (n. 676580) and the funding from the European Union Horizon 2020 research and innovation program under the Marie Sklodowska-Curie Grant Agreement No. 795246-StrongLights. The authors gratefully thank Prof. Angel Rubio for his comments and support.

{suppinfo}

• •

Supplementary Information document includes: Section SLABEL:SIsec:ComputDetails. Computational Details: Extended description the computational details used to obtain the ab initio results. Section SLABEL:SIsec:broadening. 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 SLABEL:SIsec:SupFigs and Section SLABEL:SIsec:SupTabs, supplementary figures and tables that support main text results.