# How the Interplay among Conformational Disorder, Solvation, Local and Charge Transfer Excitations affects the Absorption Spectrum and Photoinduced Dynamics of Perylene-diimide Dimers: a Molecular Dynamics/Quantum Vibronic Approach

Alekos Segalina, Daniel Aranda, James Green, Vito Cristino, Stefano Caramori, Giacomo Prampolini, Mariachiara Pastore, Fabrizio Santoro

## ► To cite this version:

Alekos Segalina, Daniel Aranda, James Green, Vito Cristino, Stefano Caramori, et al.. How the Interplay among Conformational Disorder, Solvation, Local and Charge Transfer Excitations affects the Absorption Spectrum and Photoinduced Dynamics of Perylene-diimide Dimers: a Molecular Dynamics/Quantum Vibronic Approach. Journal of Chemical Theory and Computation, 2022, 18 (6), pp.3718-3736. 10.1021/acs.jctc.2c00063. hal-03871111

HAL Id: hal-03871111

https://hal.science/hal-03871111

Submitted on 25 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
How the Interplay among Conformational Disorder,
Solvation, Local and Charge Transfer Excitations affects
the Absorption Spectrum and Photoinduced Dynamics of
Perylene-diimide Dimers: a Molecular Dynamics/Quantum
Vibronic Approach

Alekos Segalina,†,⊥ Daniel Aranda,‡,⊥ James A. Green,¶ Vito Cristino,§ Stefano Caramori,§ Giacomo Prampolini,*,‖ Mariachiara Pastore,*,† and Fabrizio Santoro*,‖

†Université de Lorraine & CNRS, LPCT, UMR 7019, F-54000 Nancy, France

‡Instituto de Ciencia Molecular (ICMol), Universidad de Valencia, Catedrático J. Beltrán 2, 46980, Paterna, Valencia, Spain

¶Consiglio Nazionale delle Ricerche, Istituto di Biostrutture e Bioimmagini (IBB-CNR), via Mezzocannone
16, I-80136 Napoli, Italy

§Dipartimento di Scienze chimiche, farmaceutiche ed agrarie, Via Fossato di Mortara 17 44121 Ferrara,
Italy

‖Istituto di Chimica dei Composti Organo Metallici, Consiglio Nazionale delle Ricerche, (ICCOM-CNR),
SS di Pisa, Area della Ricerca, via G. Moruzzi 1, I-56124 Pisa, Italy

⊥Contributed equally to this work

E-mail: giacomo.prampolini@pi.iccom.cnr.it; mariachiara.pastore@univ-lorraine.fr;
fabrizio.santoro@pi.iccom.cnr.it
Abstract

In this contribution we present a mixed quantum classical dynamical approach for the computation of vibronic absorption spectra of molecular aggregates and their nonadiabatic dynamics, taking into account the coupling between local excitations (LE) and charge-transfer (CT) states. The approach is based on an adiabatic (Ad) separation between the soft degrees of freedom (DoFs) of the system and the stiff vibrations, which are described by the quantum dynamics (QD) of wavepackets (WPs) moving on the coupled potential energy surfaces (PES) of the LE and CT states. These PES are described with a Linear Vibronic Coupling (LVC) Hamiltonian, parametrized by an overlap-based diabatization on the grounds of time-dependent density functional theory computations. The WPs time-evolution is computed with the multiconfigurational time-dependent Hartree method, using effective modes defined through a hierarchical representation of the LVC Hamiltonian. The soft DoFs are sampled with classical Molecular Dynamics (MD) and the coupling between the slow and fast DoFs is included by recomputing the key parameters of the LVC Hamiltonians, specifically for each MD configuration. This method, named Ad-MD|gLVC, is applied to a perylene diimide (PDI) dimer in acetonitrile and water solutions, and it is shown to accurately reproduce the change in the vibronic features of the absorption spectrum upon aggregation. Moreover, the microscopic insight offered by the MD trajectories allows for a detailed understanding of the role played by the fluctuation of the aggregate structure on the shape of the vibronic spectrum and on the population of LE and CT states. The nonadiabatic QD predicts an extremely fast ( ~ 50 fs) energy transfer between the two LEs. CT states have only a moderate effect on the absorption spectrum, despite the fact that after photoexcitation they are shown to acquire a fast and non negligible population, highlighting their relevance in dictating the charge separation and transport in PDI-based optical devices.
<figure><img src="image_2.png" /></figure>

Ad-MD|gLVC: a Molecular Dynamics/Quantum Vibronic Approach
# 1 Introduction

The peculiar and tunable optoelectronic properties of Perylene-diimide (PDI) based dyes and their relatively low production costs have made them attractive for several technological applications, including electronic displays, solar cells devices, and phototheranostic drugs for cancer therapy.¹⁻⁵ In addition to the possibility of tethering PDI monomers by using various linkers, it has been shown that PDIs may spontaneously self-assemble in various solutions forming large-size aggregates.⁶⁻¹⁰ In the latter case, the self-assembly is driven by a delicate interplay of non-bonded hydrophilic, hydrophobic and π-stacking interactions.⁶⁻¹⁰ Such an equilibrium is easily altered by modifying, for instance, the PDI substitution patterns at the imide position, hence effectively controlling the dimension and the shape of the aggregates.⁶ The shape of the aggregates, in turn, remarkably impacts the materials optoelectronic properties.¹¹⁻¹⁴

From a theoretical point of view, a simple connection between the structure of the supra-molecular aggregates and the optical properties is possible within the qualitative Kasha's model: cofacial arrangement leading to blue-shifted absorption relative to the monomer (H aggregate) and head-to-tail arrangement leading to red-shifted absorption (J aggregate).¹⁵ However, potentially significant contributions are neglected in the Kasha's model: i) the electron-vibrational coupling; ii) the overlap between the frontier orbitals of the chromophores; iii) the effect of interstate couplings among low energy quasi-degenerate states; iv) the non-bonded electrostatic interactions among the chromophores and with the environmental solvent molecules; v) dynamical effects. As a matter of fact, unlike inorganic crystalline or glassy materials, the optical excitations and charge transport in π-conjugated molecular aggregates involve significant rearrangements of the nuclei positions, yielding strong vibronic coupling between electronic excitations and intramolecular vibrations.¹²,¹⁶⁻¹⁸ Since the seminal works of Fulton and Gouterman,¹⁹,²⁰ where the effect of the vibronic coupling in molecular dimers was first discussed, many later works investigated the impact of vibronic coupling on molecular aggregates, showing, in many cases, its fundamental role to obtain a good agreement with experimental spectra.¹⁸,²¹⁻³⁶

In this framework, the case of the PDI is paradigmatic. The optical absorption spectra of both PDI monomer and aggregates actually exhibit an evident vibronic progression which has been connected to the coupling between the first bright electronic transition and the intramolecular totally symmetric vinyl stretching mode which has an energy of ~ 0.17 eV.¹³,³⁷ The net difference between the aggregates' spectrum and the one of the monomer stands not only in the energetic shift, as argued by Kasha, but also in the change of the relative intensity between the first two vibronic peaks ($R = I_{0→0}/I_{0→1}$), that is indeed used as spectral signature to identify the formation of PDI aggregates.¹⁰,¹²,¹³,³⁸,³⁹ As evidenced by Spano and coworkers,¹²,¹³ a large R factor ($R > 1$) is observed in J aggregate spectra, whereas a small value of R ($R < 1$) corresponds to
the case of H aggregates. The magnitude of this difference depends on the extent (number of chromophores involved), on the excitonic bandwidth, and on the nature of the exciton (pure Frenkel, Charge Transfer-like (CT), or mixed Frenkel/CT)

Furthermore, several studies have claimed that a proper description of the absorption and emission properties of PDI aggregates requires taking into account the interstate coupling, e.g. between Frenkel-like and CT-like states that compose the low-energy quasi-degenerate excited states manifold.<sup>12,14,40–43</sup> In PDI aggregates, when the interstate coupling is large enough, in particular when the Frenkel-like and CT-like states are strongly mixed, the excited-state relaxation may lead to the formation of excimers, self-trapped excitons, and to the population of CT states, which have a major role in dictating the emission properties and the exciton diffusion.<sup>2,12,40,44–46</sup> Moreover, as far as aggregates absorption spectra in solution are concerned, these appear broadened by the structural thermal fluctuations of both the self-assembled monomers and surrounding solvent molecules.<sup>47–50</sup> Indeed, recent studies have shown that the contribution of structural disorder and environmental interactions has a remarkable effect on the spectroscopy of self-assembled aggregates in solution and cannot be neglected to get an accurate description of the vibronic lineshape of the absorption spectra.<sup>35,47,48,51,52</sup> Clearly, the spectral line shapes hold a treasure trove of information regarding the nature of the fundamental electronic excitations, and valuable insights into the way in which the chromophores pack together.

Although, as reported above, a qualitative understanding of the factors ruling the change of the ab-
sorption spectrum of PDI upon aggregation has been reached, viable computational approaches to describe
the interplay among all these factors in a coherent and non-phenomenological way have been lacking until
very recently.<sup>35,52</sup> The scope of the present work is to introduce a novel Mixed Quantum Classical (MQC)
computational protocol, hereafter named "Adiabatic Molecular Dynamics with Generalized Linear Vibronic
Coupling" (Ad-MD|gLVC), aimed at simulating the spectroscopy of aggregates in solution, while taking into
account thermal fluctuations of solute-solute and solute-solvent interactions, vibronic effects, and interstate
coupling for excitonic aggregates. This approach is employed to perform a complete study of the dynamics
and the spectroscopy of the N,N-bis (2-(trimethylammonium)-ethylene) perylene-3,4,9,10,1-tetracarboxylic
acid diimide (PDI) dimer, both in ACN and in aqueous solution. Ad-MD|gLVC remarkably extends the
capabilities of the Ad-MD|gVH method, recently proposed by some of us,<sup>50</sup> introducing the possibility to
treat systems with coupled electronic states. Since in these cases no analytical expression is available for
the time-dependent correlation functions needed to compute the vibronic spectra, Ad-MD|gLVC implements
a direct numerical propagation of the vibronic wavepacket on coupled Potential Energy Surfaces (PESs).
To take into account all the effects discussed above, Ad-MD|gLVC combines extensive classical Molecular
Dynamics (MD) simulations with quantum dynamics (QD) wavepacket propagations. For the former, we use a molecule specific quantum mechanically derived force field (QMD-FF), built according to the JOYCE protocol,<sup>53-55</sup> which was recently benchmarked for the PDI monomer.<sup>49</sup> For the latter, we use the Multi-Configurational Time-Dependent Hartree (MCTDH) QD wavepacket propagations,<sup>56,57</sup> and its Multilayer (ML) generalization,<sup>58</sup> in combination with a Linear Vibronic Coupling (LVC) Hamiltonian<sup>59</sup> parametrized with a fragment-diabatization approach recently developed by some of us.<sup>60</sup> The latter technique allows for defining the adiabatic electronic states of the PDI dimer as a combination of diabatic local excitations (LE) of the PDI units and relative CT states. Similarly to the previous Ad-MD|gVH method, to account for the dynamics on the calculation of the spectra, solute+solvent nuclear degrees of freedom are adiabatically separated into stiff and soft modes, including the fast stiff modes at quantum mechanical (QM) level, while the slow soft modes are treated classically. The MD simulations sample a thermodynamically reliable set of large amplitude soft modes, and, at each snapshot extracted from the MD trajectory, the coupling between soft and stiff modes is introduced by parametrizing, in the sub-space of the stiff coordinates, specific vibronic LVC Hamiltonians describing the coupled states. The final spectrum is eventually computed as a conforma- tional average of the nonadiabatic spectra, obtained by Fourier transforming the QD correlation function of the stiff modes in each individual MD configuration.

We will show that the Ad-MD|gLVC approach permits both the reproduction of steady-state spectra, as well as the investigation of the ultrafast nonadiabatic dynamics of the photoexcited aggregate involving local and CT excitations. While here Ad-MD|gLVC is applied to PDI dimers in acetonitrile and in water, the protocol is general and can be applied to a variety of molecular aggregates, even in more complex environments than an homogeneous solvent. It may be worthy to anticipate that Ad-MD|gLVC has a significant, but still affordable, computational cost. Its strength, however, is that it is non-phenomenological so that each parameter is rigorously connected with a microscopic variable and can therefore be computed from first-principles. Due to these characteristics, Ad-MD|gLVC can allow a detailed atomistic understanding of the factors determining the photophysics of the system, permitting to examine or re-examine the robustness of phenomenological reduced-dimensionality models, and to perform the calculations in complex systems when the assumptions on which such models are based may not hold. The paper is organized as follows: Sections 2 and 3 introduce respectively the methods and the computational details. Section 4 is devoted to the presentation of the results and Section 5 reports the discussion and some concluding remarks.
# 2 The Ad-MD|gLVC Method

## 2.1 General workflow

The computational scheme presented in this work is a generalization to nonadiabatic cases of the MQC protocol named Ad-MD|gVH, briefly sketched in the left side of Figure 1, developed by some of us<sup>50</sup> and recently applied to compute the vibronic absorption spectrum of a solvated PDI monomer.<sup>49</sup> As shown in the

<figure><img src="image_3.png" /><figcaption>**Figure 1:** Comparison between the Ad-MD|gVH scheme (left),<sup>50</sup> adopted in Ref. [49] to study the PDI monomer, and the new MQC protocol (right) presented in this work for self-assembled dimers. In the boxes we highlight the main computational steps, while at right of the parenthesis we indicate the main outputs. The key differences are highlighted with a red circle, whereas the QM and classical level of theory of the calculations is indicated in blueish and greenish colours, respectively.</figcaption></figure>

right panel of Figure 1, when dealing with more than one PDI unit, the computational protocol requires a number of generalizations, the most important of which is evidenced in red, and concerns with the calculation of the final absorption spectrum. In brief, we modify the vibronic calculation engine, substituting the Vertical Hessian (VH) model,<sup>61</sup> based on non-interacting harmonic PESs, with a LVC model of coupled harmonic PESs.<sup>62</sup> As a consequence, we change the computational technique, from the evaluation of analytical time correlation functions to numerical quantum dynamical propagations. More in detail, like the Ad-MD|gVH
scheme,⁵⁰ the new protocol is based on a partition of the nuclear degrees of freedom in stiff coordinates **r** (to be treated at QM level) and much slower soft coordinates **R**, comprising the flexible coordinates of the solute and all the solvent coordinates, to be treated at classical level. The final absorption signal can then be retrieved from an average of quantum vibronic spectra computed in the stiff modes subspace, where the average is taken over the conformational space spanned by soft coordinates sampled by a classical MD run, carried out with purposely tailored QMD-FFs.⁶³,⁶⁴ Concretely, the MQC expression for the spectral lineshape is⁵⁰

$$L^{mqc}(\omega) = \frac{1}{N_{con}} \sum_{\alpha} L^{\alpha, q^r}(\omega) \quad (1)$$

where the sum is made on $N_{con}$ representative snapshots $\alpha$ taken from the classical MD trajectory. Each spectrum $L^{\alpha, q^r}(\omega)$ is computed in reduced-dimensionality normal-coordinates $q^r$ obtained projecting out all the flexible coordinates **R**, but it is specific for the snapshot $\alpha$, because the PES for **r** coordinates are re-computed at the specific value of **R** in that snapshot. In Ad-MD|gVH, the initial and final state PES along $q^r$ are assumed to be harmonic and inter-state electronic couplings are neglected. Each PES is constructed from energies, gradients and Hessians at the **R** value, a procedure that we defined as the generalized Vertical Hessian (gVH) model, since Hessians are computed out-of-equilibrium.⁵⁰ In this framework spectra can be computed straightforwardly since the necessary correlation functions are analytical.⁶⁵ Such an approach is perfectly suited for the PDI chromophore, which possesses a single bright (HOMO→LUMO) transition in the low energy range, well separated from other electronic states. However, the assumptions behind Ad-MD|gVH do not hold when addressing the spectroscopic behaviour of a PDI dimer. In fact, in this case, one needs to consider the couplings between pairs of four different diabatic states, displayed in Figure S2: a local excitation on each monomer, $|L_1\rangle$ and $|L_2\rangle$, and two charge-transfer states, arising from transitions between HOMO and LUMO located on different monomers, i.e., $|CT(1 \to 2)\rangle$ and $|CT(2 \to 1)\rangle$. To this end, as evidenced in the bottom right part of Figure 1, in Ad-MD|gLVC, the gVH model is replaced by a more complex generalized LVC Hamiltonian, parametrized at each snapshot with a diabatization procedure. Because of inter-state couplings, time correlation functions between different diabatic states cannot be computed analytically in the LVC model, but must be obtained by numerical propagation of the vibronic wavepackets on the coupled PESs according to quantum dynamics.⁶⁶ Here we adopt the MCTDH approach,⁵⁶,⁵⁷ as detailed in the following.
## 2.2 Computation of the absorption spectrum with the LVC Hamiltonian

We consider the following LVC Hamiltonian for a coupled set of diabatic electronic states $|d\rangle = (|d_1\rangle, |d_2\rangle, \dots, |d_n\rangle)$

$$H = \sum_{i} \left( K(\mathbf{p}) + V_{ii}^{\mathrm{d}}(\mathbf{q}) \right) |d_i\rangle \langle d_i| + \sum_{i,j>i} V_{ij}^{\mathrm{d}}(\mathbf{q}) \left( |d_i\rangle \langle d_j| + |d_j\rangle \langle d_i| \right) \quad (2)$$

Here $\mathbf{q}$ are the dimensionless normal mode coordinates of the ground electronic state $S_0$, and $\mathbf{p}$ are their conjugate momenta and $i$ and $j$ label the electronic states. The kinetic $K$ and potential $V$ terms of the Hamiltonian are defined as

$$K(\mathbf{p}) = \frac{1}{2} \mathbf{p}^T \Omega \mathbf{p} \quad (3)$$

$$V_{ii}^{\mathrm{d}}(\mathbf{q}) = E_{ii}^{\mathrm{d}}(\mathbf{q}_0) + \lambda_{ii}^{\mathrm{T}} \mathbf{q} + \frac{1}{2} \mathbf{q}^{\mathrm{T}} \Omega \mathbf{q} \quad (4)$$

$$V_{ij}^{\mathrm{d}}(\mathbf{q}) = E_{ij}^{\mathrm{d}}(\mathbf{q}_0) + \lambda_{ij}^{\mathrm{T}} \mathbf{q}, \quad (5)$$

where $\hbar = 1$, $\Omega$ is the diagonal matrix of normal mode frequencies $\omega_k$ for mode $k$, $E_{ii}^{\mathrm{d}}(\mathbf{q}_0)$ are the diabatic energies and $E_{ij}^{\mathrm{d}}(\mathbf{q}_0)$ the electronic coupling constants between diabatic states at the reference geometry $\mathbf{q}_0$. Notice that these latter coupling constants do not appear in the standard LVC approach,⁵⁹ but are expected to exist and be important in Exciton/Charge-Transfer diabatic models. The vectors $\lambda_{ii}$ and $\lambda_{ij}$ (with $j \neq i$) are the gradients of the diabatic PESs and of the inter-state couplings, respectively.

In a TD framework, the absorption spectrum $\epsilon(\omega)$ at 0 Kelvin can be expressed as

$$\begin{aligned} \frac{\epsilon(\omega)}{C_{\epsilon}} &= \omega \sum_{ji} \int_{-\infty}^{\infty} dt e^{i\omega t - \Gamma t^2} \langle 0; d_j | \boldsymbol{\mu}_{gj} e^{-i\mathbf{H}t} \boldsymbol{\mu}_{ig} | d_i; 0 \rangle \\ &= \sum_i \epsilon_{ii}(\omega) + \sum_{i,j \neq i} \epsilon_{ij}(\omega) = \epsilon^{\mathrm{auto}}(\omega) + \epsilon^{\mathrm{cross}}(\omega) \end{aligned} \quad (6)$$

where $\boldsymbol{\mu}_{gj} = \langle g | \boldsymbol{\mu} | d_j \rangle$ are the matrix elements of the electric dipole moment, $0$ is the ground-vibrational state of the ground electronic state $g$, whose energy is set to 0, and $C_{\epsilon}$ contains all the physical constants (its expression and its value to give molar absorptivity in units $M^{-1}cm^{-1}$ can be found in refs. [67] and [68]). Moreover, a quadratic damping term ruled by a parameter $\Gamma$ was introduced, corresponding to a Gaussian broadening in the frequency domain. By construction diabatic states are built so to be virtually independent of the nuclear coordinates. This implies that the elements $\boldsymbol{\mu}_{gj}$ can be considered constant (Condon approximation) and the problem can be reduced to the computation of the time correlation functions

$$\phi_{ij}(t) = \langle d_i(0) | d_j(t) \rangle \quad (7)$$
This step must be performed numerically propagating the initial wavepackets $|d_j(0)\rangle$ on the coupled PESs of the LVC Hamiltonian. The cross-correlation functions $\phi_{ij}(t)$ ($j \neq i$) are required to obtain the terms $\epsilon_{ij}(\omega)$, which modulate the spectral shape, although they do not contribute to the total intensity. Once the different $\phi_{ij}(t)$ are computed, the total correlation function, $\phi_{tot}(t)$, is obtained by a weighted sum where the weight on each term is the scalar product of the electric transition dipole moments of the corresponding states:

$$ \phi_{tot}(t) = \sum_{i,j} \mu_{ig} \cdot \mu_{gj} \phi_{ij}(t) \qquad (8) $$

The final absorption spectrum is eventually obtained by Fourier transforming $\phi_{tot}(t)$ (Eq. 6). The derivation of analytical expressions for the total intensity ($M_0$) and the center of gravity (first moment, $M_1$) of the absorption lineshape predicted by our model is straightforward (see Section S1.3 in the Supporting Information). For an ideal PDI dimer with two degenerate local states, $|L_1\rangle$ and $|L_2\rangle$ and two completely dark CT states, $|CT_{12}\rangle$ and $|CT_{21}\rangle$, $M_1$ becomes

$$ M_1 = E_{L_1}^d(0) + E_{L_1,L_2}^d(0) \cos \alpha \qquad (9) $$

where $E_{L_1}^d(0)$ is the first moment of the monomer spectrum and $\alpha$ is the angle between the transition dipoles of the two states. In a stacked dimer fully eclipsed ($\alpha=0$) the so-called “exciton” electronic coupling $E_{L_1,L_2}^d(0)$ is expected to be positive and therefore the spectrum of the dimer blueshifts with respect to the monomer by the extent of $E_{L_1,L_2}^d(0)$.

## 2.3 Implementation

### Ad-MD|gLVC full route

As displayed in Figure 1, the new Ad-MD|gLVC scheme does not significantly differ from the previous Ad-MD|gVH procedure in the production of a set of reliable snapshots through the MD runs. As already mentioned, the main difference stands in the computational protocol applied to each sampled snapshot $\alpha$. All the steps required by the full (i.e. without any further assumption) Ad-MD|gLVC procedure are sketched in the left side of Figure 2. It should be noticed that, while in Ad-MD|gVH we account for differences in the quadratic terms of the final and initial PES, so that they can have different normal modes and frequencies, in Ad-MD|gLVC we assume that all diabatic states have the same Hessian as the ground state (GS), see Eq. 4. An analogous approximation transforms the VH model into the simpler Vertical Gradient (VG) one.⁶¹ In order to account for the effect of these quadratic terms, we would need to replace the LVC model by the so-called Quadratic Vibronic Coupling (QVC) one,⁶² but this is computationally too expensive, and LVC
is enough to capture the main effects we want to study. Conversely, considering that the straightforward
application of the full Ad-MD|gLVC protocol is also rather time consuming, we have also designed a much
faster simplified procedure, displayed in green in the right branch of the same Figure.

In practice, according to the full Ad-MD|gLVC scheme, for each snapshot α the following steps are
necessary:

1. Reduced QM Hessian and ground state equilibrium geometry $q_0^r$

The energy, gradient and Hessian of the GS are computed at a proper QM level. The flexible coordi-
nates are projected out and considered frozen, normal modes **q**^r describing only the stiff coordinates
are computed and a GS equilibrium geometry along these modes is extrapolated assuming harmonic
approximation. Such geometry **r**^α<sub>eq</sub>, as well as the normal modes **q**^r,α and their frequencies are specific
for the snapshot α.

2. FrD-LVC diabatization

The full LVC Hamiltonian is built up by displacing the $\mathbf{r}_{eq}^{\alpha}$ for positive and negative steps $\Delta_k$ along

<figure><img src="image_4.png" /><figcaption>**Figure 2:** Ad-MD|gLVC workflow of both the full route (left side) and the simplified one (right side). The red boxes indicate the computationally more demanding steps, whereas green boxes indicate the steps where the cost in terms of CPU time is significantly reduced. See text for definition of the employed acronyms and symbols.</figcaption></figure>
each mode $q_k^{r,\alpha}$, performing a TD-DFT computation and a fragment-based diabatization (FrD) with a maximum-overlap criterion implemented in our "in house" code *Overdia*. For this reason we indicate this LVC Hamiltonian with the label FrD-LVC. This step (evidenced in red in Figure 2) is the most time consuming because it requires $2 \times N_{vib}^r + 1$ TD-DFT computations for the dimer, (where $N_{vib}^r$ is the number of stiff normal modes), plus the computation of the necessary overlaps between the transition densities necessary to apply the diabatization scheme.

### 3. Hierarchical Effective Mode Selection (HEMS)

The low- medium-resolution spectrum can be computed in a very efficient way by adopting a Hierarchical transformation of the LVC Hamiltonian in blocks of sequentially coupled effective modes.⁶⁹⁻⁷¹ Such a hierarchy is obtained with a generalization of the Lanczos algorithm.⁷² Few blocks are sufficient to converge the spectrum.

### 4. MCTDH

Nonadiabatic QD propagations of the vibronic wave packets are performed with the MCTDH method and time-dependent correlation functions are computed.

### 5. Absorption spectrum

Absorption spectra for each snapshot are obtained by Fourier transform of the total correlation function and the final average spectrum is obtained according to Eq. 1

#### Ad-MD|gLVC simplified route

As commented in the previous list, the generation of the LVC Hamiltonian is the slow part of the protocol and the repetition for a representative number of snapshots (~100) demands a significant computational effort. Yet, two observations, evidenced in green in Figure 2, allow for a remarkable speed up. In the following we will prove in fact that for the PDI dimer is possible to neglect the linear terms in the inter-state couplings in Eq. 5 (i.e. $\lambda_{ij}^{\alpha} = 0$), and to assume that the diagonal gradients $\lambda_{ii}^{\alpha}$ are independent of the specific snapshot $\alpha$. The first approximation implies that we assume that the inter-state couplings depend on the snapshot $\alpha$ since they depend on the fluctuations of the slow coordinates, but they are independent of the small oscillations of the fast coordinates $q^r$. As far as the second approximation is concerned, different strategies can be conceived to estimate the snapshot-independent gradients $\lambda_{ii}$. The protocol that we test and adopt in this work uses only computations for an isolated PDI monomer, optimized in its ground state in the gas phase. More specifically we assume that the relevant normal modes of the PDI dimer can be represented by the sum of the sets of the normal modes of two isolated monomers in gas phase $\{q^{M1}, q^{M2}\}$, hence computing the gradient of the local excitation of the monomer ($g^L$), and of the CT cation ($g^C$) and anion
(g^A). In this framework, the gradients of the two local excitations |L₁⟩ and |L₂⟩ of the dimer and the ones of the two CT states, CT(1→2) and CT(2→1), can be eventually approximated by the vectors {g^L, 0} and {0, g^L}, and by the vectors {g^C, g^A} and {g^A, g^C}, respectively. It is noteworthy that the above gradients can be computed analytically and therefore much more rapidly than what required for the parameterization of the LVC Hamiltonian according to the full-route (see Section 3.4). In practice, in this approximation only the diabatic vertical energies E<sub>ii</sub>(q<sub>0</sub><sup>r,α</sup>) and inter-state couplings E<sub>ij</sub>(q<sub>0</sub><sup>r,α</sup>) (j ≠ i) need to be recomputed at each snapshot α. This operation is much faster than the generation of a full LVC Hamiltonian as it requires a single TD-DFT computation for the dimer and one for each of the two monomers, resulting in a huge speedup with respect to the full protocol. We notice that such a "single-point" diabatization was first presented by some of us with the name "FrD-EX",<sup>60</sup> and it was used to parameterize an purely electronic excitonic model with only the terms E<sub>ii</sub> and E<sub>ij</sub>. Therefore, in this new scheme in Figure 2, step 2 is replaced by the FrD-EX diabatization. Moreover, the generation of the HEMS is identical for all snapshots and can be performed just once, hence replacing step 1, which was carried out for each snapshot α. More technical details on the different steps of the computations are given in the following section.

**3 Experimental and Computational Details**

**3.1 Experimental**

UV-Visible spectra were acquired in transmission mode with an Agilent-Cary 300 UV-Visible Spectropho-
tometer in the 800-200 nm interval at a scan rate of 600 nm/min, with a spectral bandwidth of 2 nm and
employing a quartz cell with an optical path of 2 mm. (PF<sub>6</sub>)<sub>2</sub>-PDI salts are not soluble in water, whereas
they are soluble in acetonitrile (ACN), where however aggregation cannot be observed even for saturated
solutions. Therefore, in order to observe PDI aggregation, 20 µl aliquots of 5 mM PDI in ACN (PDI@ACN)
were progressively added to 2 ml of water up to reaching a total 100 µl ACN in water. After this point a
further final addition of 50 µl 5 mM PDI@ACN was accomplished. Alternatively, difference spectra reflecting
the spectral changes due to PDI aggregation can be obtained by starting from a mother solution of PDI in
ACN and adding increasing amounts of water. All difference spectra, along with the related experimental
details, can be found in Section S2.1 of the Supporting Information.

**3.2 Molecular Dynamics**

The absorption spectrum of a Cl₂-PDI salt solvated in ACN, which can be recorded experimentally,<sup>10</sup>
was previously simulated by us by means of the Ad-MD|gVH method, taking into account a single PDI
monomer.⁴⁹ Here we instead consider two Cl₂-PDI units and compute Ad-MD|gLVC spectra in both ACN and water. Despite the computational burden, this choice was made, on the one hand, for an internal coherence when comparing the vibronic progressions of dimer and monomer (ACN) and, on the other hand, in view of a more realistic comparison with the experimental spectrum of the aggregates, obtained in water. It is worthwhile to stress that Cl₂-PDI solvation and PDI self-aggregation is observed with MD in both environments, because of the resulting high concentration, notwithstanding the significant number (1·3·10³) of solvent molecules considered. The fact that the counterion dissociates in our simulations makes immaterial the use of Cl⁻ or PF₆⁻, thus we used Chlorine for internal coherence with the MD simulations previously performed for the monomer.⁴⁹ As discussed in detail in the following, negligible differences arose in the final computed spectra between the two solvents, hence, for the sake of brevity, most of the results in water will be reported in Supporting Information, whereas the discussion in the main text refers, unless otherwise stated, to simulations in ACN.

An accurate QMD-FF for a PDI monomer was previously parameterized by some of us,⁴⁹ as briefly summarized in the following. All intra-molecular PDI's parameters were obtained with the JOYCE code.⁴⁹,⁵³,⁵⁵ As far as the PDI's intermolecular terms are concerned, the point-charges were obtained through the RESP protocol by using the Antechamber suite, while the Lennard-Jones (LJ) parameters were transferred from the OPLS libraries.⁴⁹,⁷³,⁷⁴ The accuracy of the description of the intermolecular interactions within the dimer, which were not considered in our previous work on PDI monomer,⁴⁹ is here carefully validated through the analysis of the MM Interaction Potential Energy Surface (IPES) along selected dimer arrangements. For solvents and counterions, the TIP3P model was employed for water, whereas the parameters concerning ACN and Cl⁻ were taken from the OPLS FF.⁷³,⁷⁴ Further details on PDI's QMD-FF can be found in the original paper.⁴⁹

As far as MD simulations are concerned, two systems, PDI₂@ACN and PDI₂@H₂O, were prepared by solvating into a cubic box two PDI units with ~ 1000 ACN or ~ 3000 H₂O molecules, respectively. In both systems, four Cl⁻ anions have been also included to ensure the electroneutrality of the cell. The two PDIs were placed randomly in the simulation boxes at a distance between the centers of mass of about 16 Å. Moreover, two sets of MD runs were carried out for PDI₂@ACN, separately for the *syn* and the *anti* conformer, aimed to rationalize the role of the lateral pendants on self-aggregation. All MD simulations (in periodic boundary conditions) have been computed making use of the GROMACS5.1 Engine. The PDI₂@ACN and PDI₂@H₂O systems were initially minimized to avoid bad contacts and, subsequently, thermally equilibrated for ~2 ns at 300 K in the NVT ensemble. Very long production runs were carried out in the NPT ensemble for 1 µs at 1 atm and 300 K and considering a time step of 1 fs using the LINCS algorithm⁷⁵ to fix the
bond distances involving hydrogen atoms. The temperature and the pressure coupling have been described
through Parrinello-Raman and the v-rescale schemes using coupling constants of 0.1 ps and 1 ps, respectively.
Finally, for the short-range Coulombic and LJ terms the cutoff radius was set to 11 Å, while the long-range
electrostatics interactions were treated using the particle mesh Ewald (PME) procedure.

## 3.3 QM calculations

All QM calculations were carried out at the DFT and TD-DFT level, using GAUSSIAN16 package<sup>76</sup> with the CAM-B3LYP<sup>77</sup> functional and the 6-31G(d) basis set, adopting the Grimme's D3 dispersion correction.<sup>78,79</sup>
The QM interaction energy between two monomers, necessary to benchmark the MM IPES, was computed along the coordinates displayed in Figure 3a, including the solvent effect at C-PCM level to screen the electrostatic repulsion due to the +2 charge of the lateral chains. The same solvent description was applied in an unconstrained optimization of the ground-state PDI-dimer geometry, employed as reference to parametrize the full FrD-LVC Hamiltonian, performing TD-DFT calculations on molecular structures displaced along all the reduced-dimensionality normal modes.

The QM calculations required for the MQC spectrum were carried out for each snapshot α extracted
from the MD runs as described in the following. First, the reference geometry, $q_o^{r,α}$, optimized in the reduced
space of the fast coordinates *r* was computed at QM:MM level as described in Section S1.4 in the Supporting
Information. Therein, we also show for a few selected snapshots, alternative computational protocols to run
the optimization lead to only marginal differences on the determination of the $q_o^{r,α}$ structures and therefore of
the diabatic energies *E*<sub>ii</sub>(0) and couplings *E*<sub>ij</sub>(0) for the PDI dimer. Next, to apply the FrD-EX diabatization
scheme, single-point calculations were carried at TD-DFT level on each re-optimized geometry $q_o^{r,α}$ on both
the monomers and the dimer. Calculations on the monomers are used to compute the reference states
to define the diabatic states, whereas calculations on the dimer were performed to obtain the adiabatic
states which are then rotated to overlap as much as possible the reference states. This operation defines the
diabatic-to-adiabatic transformation, which was computed with a freely-distributed code, *Overdia*, developed
by some of us.<sup>80</sup> Further details are given in the following section.

## 3.4 Fragment Diabatization

In all monomer and dimer calculations a two layer approach was adopted for taking the environment into
account, again considering the electrostatic embedding (EE) of all solvent molecules and counterions within a
radius of 20 Å which are included as QMD-FF point charges (pc). Furthermore, in the monomer calculations,
the electrostatic effect of the presence of the second monomer on the reference states was introduced by
including it in the EEpc shell. We considered four diabatic states for the PDI dimer: two local excitations of the individual monomers and their two associated CT states. Reference states for the local excitations were obtained by reading the corresponding TD-DFT response vectors in the monomer calculations. As reference for the CT states, we defined states made up by an orbital transition from the HOMO of the first monomer to the LUMO of the second one, and viceversa. Afterward, according to a maximum-overlap based fragment diabatization, we determined the combination of the adiabatic excited-states of the dimer that resemble as much as possible the reference states (see Supporting Information for a sketch). This procedure^60 defines an adiabatic-to-diabatic transformation matrix (D). Application to the diagonal matrix of the adiabatic (TD-DFT) energies of the dimer, allows for computing the diabatic vertical energies $E_{ii}(q_0^{r,\alpha})$ and couplings $E_{ij}(q_0^{r,\alpha})$. For the simplified route in Figure 2 (right) this operation needs to be done just once for each snapshot $\alpha$. In fact the inter-state couplings are assumed to be independent of the oscillations of the fast coordinates $q^r$, and the gradients of the diabatic diagonal potentials along $q^r$ are taken to be independent of the snapshot $\alpha$. These gradients are computed analytically with DFT for the cation and the anion of the monomer for the CT states, and with TD-DFT for the monomer local excitation. On the contrary, for application of the full-protocol in Figure 2 (left), the gradients of the diagonal and off-diagonal LVC diabatic potentials, ($\lambda_{ii}^\alpha$ and $\lambda_{ij}^\alpha$, respectively) were obtained by numerical differentiation.^81 To do that, for each snapshot we performed $2 \times N_{vib}^r$ additional calculations, each at a structure displaced with respect to the equilibrium position $q_0^{r,\alpha}$ by a small quantity $\Delta_k = \pm 0.02$ along a single normal coordinate $q_k^\alpha$, and then applied the transformation (D$^\alpha$) to the elements of the matrix of the adiabatic potential energies $V^{ad,\alpha}(\Delta_k)$:

$$V^{d,\alpha}(\Delta_k) = (D^\alpha)^T V^{ad,\alpha}(\Delta_k) D^\alpha \qquad (10)$$

$$\lambda_{ij}^{\alpha}(k) = \frac{V^{d,\alpha}(\Delta_k) - V^{d,\alpha}(-\Delta_k)}{2\Delta_k} \qquad (11)$$

## 3.5 Quantum dynamics of the nuclear wavepackets

Nuclear wavepacket propagations were performed using the MCTDH method^56,57 as implemented in the Quantics code, with the specific settings shown in the SI.^82,83 There are 224 fast coordinates of the PDI dimer, making a straightforward application of MCTDH very challenging. In these cases, fully-converged low-resolution spectra, which account for the effect of all nuclear coordinates, can be obtained exploiting a hierarchical representation of the Hamiltonian, in terms of effective collective coordinates. They are divided in blocks, each comprising a number of coordinates, defined in such a way that the short-time dynamics (the only-one relevant for the low-resolution spectrum) is dominated by a few blocks.^69–72,84 In Figure S5 of the Supporting Information we show that 3 blocks (12 coordinates) already provide converged spectra
for a Gaussian broadening with an half width at half maximum (HWHM) of 0.03 eV, and that including a
fourth block results in very minor changes. Therefore, in the following, 3 blocks are adopted to compute the
spectra of all snapshots. The number of coordinates necessary to accurately describe the dynamics increases
with the time range monitored after the photoexcitation. Hence, in order to ensure fast and fully-converged
results, we used the Multi-Layer (ML) extension of MCTDH,⁵⁸ ML-MCTDH. Indeed, this extension allows
us to include even the effect of all the normal coordinates (44) with sizable couplings. QD simulation in 44
dimensions were actually run to compute the time-evolution of the electronic populations up to t = 100 fs,
discussed in Section 4.5.

4 Results

4.1 Validation of the Ad-MD|gLVC simplified route

As discussed in Section 2, the significant cost of the full Ad-MD|gLVC protocol prompted us to devise the
more computationally effective route displayed in Figure 2. For the specific case of the PDI dimer, the
approximations over which such a simplified procedure is rooted were tested with care before systematically
applying them along the MD trajectories. To this end: (i) we considered a single dimer geometry in the
*anti* conformation, namely the optimized structure in ACN; (ii) we applied the full route by projecting
out ⁴⁹,⁵⁰ the coordinates of the whole flexible lateral chains (all bonds, angles and dihedrals) as well as the
inter-monomer coordinates;⁸⁵ iii) on this reduced dimensionality model (from 450 to 224 normal modes) we
parametrized a full LVC Hamiltonian and computed the "static" spectrum (see Section S2.2 in the Supporting
Information for further details). The absorption spectrum was thereafter also computed with the simplified
route described in Section 2.3, obtained by setting $\lambda_{ij} = 0$ and $\lambda_{ii}^{\alpha} = \lambda_{ii}[q_0] \forall \alpha$. While further details
are available in Section S2.3 in the Supporting Information, the results, displayed in Figure S10, show that
imposing $\lambda_{ij} = 0$ has almost no effect on the spectrum, whereas the second approximation only introduces
modest changes. All dimer spectra presented in the following sections are therefore computed adopting both
the approximations of the simplified route.

4.2 PDI dynamics in solution: Structural properties

After validating the reliability of the LVC simplified route, we go beyond the “static” approximation, in-
troducing the effect of the dynamics of both the PDI dimer and the solvent on the spectrum through MD
simulations. PDI shows two distinct isomers (*syn* and *anti*), depending on the relative position of the flexible
lateral pendants with respect to the plane defined by the *π*-core. As recently shown for the monomer case,⁴⁹
although the two isomers have almost the same energy, the inter-conversion from one isomeric form to the other requires overcoming a rather high energy barrier. For this reason, both the syn and anti isomers were separately considered. Despite the fact that the lateral pendants have a negligible effect on the monomer optical properties,⁴⁹ in the case of aggregates, they are expected to play a significant role, by influencing the stability and the organisation of the supramolecular π stacked aggregate, thus indirectly affecting the final spectral shape.

The first point to assess is therefore the accuracy of the intermolecular FF in describing the aromatic π stacking interactions, which significantly depend on the choice of the FF parameters.⁸⁶ To this end, we compare the MM and the IPESs computed at QM level, along the most relevant inter-monomer coordinates, displayed in Figure 3a. As shown in Figure S1 in the Supporting Information, they are defined according to the molecular reference axes of each PDI aromatic core. Thereby, we define the distance between the center of the two PDIs as ρ, which may be decomposed, as shown in Figure 3a, as the distance between the two PDI planes, R<sub>π</sub>, and the displacement along the short and the long axis: R<sub>s</sub> and R<sub>||</sub>, respectively. Finally, we call α the angle between the two long axes, and β (see Figure S1) the angle between the vectors normal to the PDI planes. It is worthy highlighting here that, when both monomers bear the lateral pendants in syn conformation, the aromatic cores are able to get closer, and a larger portion of the IPES can be explored. Therefore, the FF IPES validation was carried out through the analysis of the above mentioned geometrical descriptors and applied to the syn conformers.

In Figure 3b, we compare IPES's profiles along the selected geometrical descriptors, computed at the QM and MM level, respectively. As far as the distance between the centers of the aromatic planes is concerned, Figure 3b's top panel shows that both approaches predict the lower interaction energy for ρ ~ 3.6 Å, in agreement with the value expected for π-stacked systems.³²,⁸⁷ Hence, all the other IPES scans were then obtained keeping ρ fixed at 3.6 Å. In the second panel of Figure 3b, the QM intermolecular energy scan along the twist angle α shows two nearly degenerate minima, at ~160° and ~30°, whose position is well reproduced by the MM scan, although with opposite relative stabilities. Also, the energy barrier separating the two minima is slightly overestimated at the MM level. Finally, for the IPES scans along the R<sub>s</sub> and R<sub>||</sub> directions, displayed in the bottom panels of Figure 3b, we kept the two PDIs eclipsed, hence α = β = 0° and ρ and R<sub>π</sub> coincide (see Section S1.1 in the Supporting Information for further details). DFT and MM scans both indicate that the attractive interactions along both axes are greatest for displacements of around 2 Å. These minima are slightly shallower at MM level. On overall, the MM forces are expected to reliably reproduce the QM reference ones, confirming the accuracy of our FF in describing the interaction between monomers.
<figure><img src="image_5.png" /><figcaption>**Figure 3:** a) Structure of the PDI dimer, where the arrows highlight the angles, the distances, and the displacements analyzed in this work. b) Comparison between QM torsional relaxed energy scans (red lines) and the relative MM relaxed profiles (blue lines) computed for the PDI dimer in ACN solution.</figcaption></figure>

The validated FF was employed in NPT simulations of the solvated PDI dimer in ambient conditions to retrieve the statistical distribution of aggregate configurations required by the Ad-MD|gLVC method. To investigate the effect of side chains on the shape of the sampled geometrical arrangements, two separate starting configurations were considered, where the lateral pendants of both the PDI monomers were placed in either *anti* or *syn* conformation. The time evolution of the geometrical descriptors defined in Figure 3a was monitored along both trajectories, with reference to data in ACN, as this is the only solvent for
which we run MD simulations considering two different starting configurations. It is important to mention that during the long µs dynamics starting with the separated monomers both bearing pendants in *anti* conformation (*anti-anti* dimer), one monomer undergoes a *anti* to *syn* transition, thus settling in a *anti-syn* arrangement. In Figure S11 in the Supporting Information we show that only one transition is observed in ACN during the µs dynamics and it takes place between 700 and 800 ns. Therefore, as evident by the chain conformational distribution shown in the same Figure, when the monomers both start in the *anti* conformation, roughly 80% the dimer configurations are found in the *anti-anti*, while the rest is in an *anti-syn* arrangement. Nonetheless, since in the other MD run starting with a *syn-syn* configuration the opposite *syn* to *anti* transition is never observed, for the sake of simplicity the two MD runs will be still labeled as *anti* or *syn* trajectories. Most importantly, it should be well stressed that the *anti* → *syn* transition takes place well after the PDI aggregate formation, which takes place much more rapidly. In fact, in agreement with the conclusions drawn from the IPES analysis, Figure S14 in the Supporting Information shows the time evolution of ρ, which clearly indicates that, independently from the starting configuration, the two monomers form a stable aggregate in less than 1 ns.

The top panel in Figure 4a shows the distribution along the *anti* MD trajectory of the distances between the $\pi$ planes ($R_{\pi}$) and between the center of the PDI cores, ($\rho$). The former shows a maximum around 2.6 Å and a long tail at larger distances, whereas the latter shows a more narrowed distribution with a maximum around 3.9 Å. The difference between these two distributions suggests that PDIs are almost never found perfectly superimposed. This is confirmed by the distributions of both the displacements along the short ($R_{S}$) and the long ($R_{\|}$) axes (the middle panel of Figure 4a), which show a maximum around 2.2 Å, and the ones of the $\alpha$ and $\beta$ angles (bottom panel in Figure 4a), which indicate that the two PDI are always twisted ($\alpha \approx 20°-50°$), and slightly tilted ($\beta \approx 0°-15°$). The correlation among the different geometrical descriptors is evidenced in Figures S15-S17 in the Supporting Information, where 2D heat maps are displayed for the $\rho$ distance with the considered angles and displacements, computed along the *anti* trajectory. Although the most probable $\rho$ distance is around 3.9 Å, it is evident that, by displacing, twisting or tilting one monomer with respect to the other, the dimer is able to visit regions where the two planar cores are closer, and a not negligible population is registered also in the 3.6-3.9 Å range.

By comparing the distributions of the considered geometrical parameters obtained in the *anti* and *syn* MD runs, one can assess the role of the lateral pendants in modulating the aggregate's structure. The distributions of $\rho$, $R_{\pi}$, $R_{S}$, and $R_{\|}$, displayed in Figure 4b, are essentially unchanged due to the position of the lateral chains. $R_{\pi}$ has a right skewed distribution with a maximum around 2.6 Å, while $\rho$ has an almost symmetric distribution, slightly narrower compared to the one of Figure 4a, whilst still having a
<figure><img src="image_6.png" /><figcaption>**Figure 4:** Distributions of inter-molecular (PDI-PDI) angles, distances, and displacements as defined in Figure 3 computed along the MD trajectories of the PDI₂@ACN starting with the a) *anti* and b) *syn* isomers.</figcaption></figure>

maximum around 3.9 Å. Similarly, R<sub>S</sub> and R<sub>||</sub> have very left skewed distributions and show maxima around 2.2 Å. Again, as evidenced from the ρ-R<sub>π</sub> correlation in Figure S18, the two PDI cores are rarely found superimposed. The main difference between the two isomers is instead evident from the distributions of both α and β (bottom panels of Figure 4b. At variance with the *anti* trajectory, the *syn* isomers form aggregates which populate more eclipsed (α <20°) dimer arrangements, as evidenced from the distribution in the bottom panel of Figure 4b). In turn, this again reflects in the possibility for the two monomers to reach closer inter-core distances, shown by the cross correlation displayed in Figures S15 and S18 in the
Supporting Information. Conversely, as far as the $\beta$ angle is concerned, the two PDI planes are more likely to be parallel when considering the syn isomer, as demonstrated by the narrowed symmetric distribution of $\beta$ centered around 0° in Figure 4b.

Finally, with respect to ACN, the results in water displayed in Figures S13 and S21 in the Supporting Information, show rather similar distributions of the chain conformation and the inter-molecular descriptors, respectively. The only noticeable difference stands in the anti → syn transition of one of the two PDIs, which occurs more rapidly, in less than 100 ns. Consequently, the two PDI planes are in average more parallel, as evidenced by the distribution of $\beta$ angle, and more likely closer and stable, as illustrated by the narrowed distribution of $\rho$. Considering the results achieved in both solvents, MD simulations indicate that the PDI dimer is able to visit a plethora of accessible configurations, whose shape and incidence are driven by the subtle interplay of the conformation of the lateral pendants, the interaction with the solvent, and thermal fluctuations.

## 4.3 Parametric study of stacked dimers

The structural fluctuations of the dimer discussed above are expected to impact the absorption spectrum and hence the excited state dynamics through the electronic coupling between LE and CT states. While a detailed discussion of the effect of the supramolecular dynamics on the final spectroscopic properties will be given in Sections 4.4 and 4.6, here we first report a simple parametric study, aimed at rationalizing the effect of each of the geometrical descriptors introduced in the previous section on diabatic energies ($E_{ii}$) and couplings ($E_{ij}$). Specifically, we will focus our attention on the dependence of the excitonic coupling, $E_{LE-LE}(0)$, the LE-CT coupling, $E_{LE-CT}(0)$, and the CT-LE diabatic energy gap $\Delta E_{CT-LE}(0)$ upon varying the stacking and sliding distances $R_{\pi}$, $R_{\|}$, $R_{S}$, the co-planar angle $\alpha$ and the tumbling angle $\beta$. Notice that $\beta$ actually is a solid angle and therefore for simplicity we only considered the tumbling angles along the sliding directions, defined here as $\beta(R_{\|})$ and $\beta(R_{S})$. In the qualitative spirit of this extended analysis, the dimers in this section are built with two ideal PDI monomers where lateral chains have been substituted by hydrogens. Each monomer is in its equilibrium geometry ($D_{2h}$ symmetry) and the starting reference consists in two coplanar and perfectly superimposed monomers, from which each inter-molecular displacement is separately applied.

Figure 5 shows the dependence of the coupling terms on the geometrical descriptors, which were varied in the range suggested by the previous MD analysis. According to these results – and with the expected exception of the inter-chromophore distance $R_{\pi}$ – the dependence of the $E_{LE-LE}(0)$ coupling is in general small with the structural parameters, and ranges between 0.08-0.12 eV. On the other hand, the $E_{LE-CT}(0)$
<figure><img src="image_7.png" /><figcaption>**Figure 5:** Dependence of the constant interstate coupling $E_{ij}(0)$ (for $i \neq j$) and the LE-CT energy gap ($\Delta E_{CT-LE}(0)$) with different parameters of the PDI dimer: (a) $R_{\pi}$, (b) $R_{S}$, (c) $R_{\|}$, (d) $\alpha$, (e) $\beta$ (component along $R_{S}$) and (f) $\beta$ (component along $R_{\|}$). The scans were performed on the structure with $R_{\pi} = 4.0$ Å, changing one parameter at time and setting all others zero. Since several choices of $E_{LE-CT}(0)$ are possible, we specify that only the $E_{LE-CT}$ coupling between the states $|L_1\rangle,|CT(1 \to 2)\rangle$ is reported because it is larger in magnitude. By symmetry, the $|L_2\rangle,|CT(2 \to 1)\rangle$ coupling is the same.</figcaption></figure>

coupling is more sensitive to the structural changes of the dimer: it shows large fluctuations as the structure changes on $R_S$ and $\beta(R_\|)$, and a very large decrease with increasing $R_\pi$, but it is almost independent of $\beta(R_S)$. Interestingly CT states are less stable than LE ones but the energy gap decreases at the decrease of the three distances $R_S$, $R_\|$ and, especially, $R_\pi$. For $R_\pi < 3.9$ Å, this gap is smaller than the coupling $E_{LE-CT}(0)$, causing a large mixing of LE and CT states in the adiabatic states of the dimer. The absorption spectra corresponding to dimer structures considered in these scans are shown in Figures S23 of the Supporting Information. The effects are complex as all the couplings change at the same time. On the one hand, it is evident that larger couplings result in a red-shift, a decrease of the intensity of the lowest energy band, and an increase and blue-shift of the most intense peak. On the other hand, when the coupling decreases, the spectra progressively become similar to the one of the monomer, characterized by a maximum of the intensity for the lowest energy band.
4.4 Ad-MD|gLVC Absorption spectrum of the dimer in solution

We now move to the discussion of the PDI dimer spectrum in ACN and in water, computed with the Ad-MD|gLVC method, which accounts for the slow dimer and solvent dynamics, by averaging over 100 configurations of the soft coordinates sampled with the MD runs carried out in each solvent. Figure 6 shows

<figure><img src="image_8.png" /></figure>

**Figure 6:** Top: Ad-MD|gLVC for PDI dimer in ACN (from *anti* and *syn* MD runs) and water (*anti* MD trajectory). Bottom: comparison of Ad-MD|gVG (monomer) and Ad-MD|gLVC (dimer) spectra with the experimental spectra in ACN and water, respectively. The calculated spectra were redshifted by 0.25/0.27 eV and all stick transitions were convoluted with a Gaussian of HWHM = 0.005/0.03 eV for monomer/dimer, respectively. In the top panel, the spectra were scaled to have unit area and then rescaled by the same factor. For bottom panel, the calculated spectra were scaled to match the experimental intensities.
the thermally-averaged spectra of the PDI dimer computed along the *anti* and *syn* trajectories in ACN, and compares them with the spectrum obtained for the monomer in the same solvent through the corresponding Ad-MD|*g*VG method (see Section S2.8 in the Supporting Information for further details). As displayed in the top panel, the computed spectra with *anti* and *syn* trajectories are very similar, and same holds for the *anti* spectra in ACN and water. In the bottom panel of Figure 6, our results are also compared with the experimental ones for the monomer in ACN¹⁰ and for the PDI aggregate, obtained in this work by adding few drops of a PDI@ACN solution to water. Given the similarities, in this panel for the dimer we only report the computed spectrum of the *anti* trajectory in water. To ease comparison, computed spectra of the dimer have all been artificially shifted by 0.27 eV in order to match the maxima of the computed spectra of the dimer and the experimental spectrum of the aggregate (analogously the spectrum of the monomer was shifted by 0.25 eV).

The Ad-MD|*g*LVC spectrum of the dimer reproduces very well the existence of two bands of the experimental lineshape of the aggregate, at 2.3 and 2.5 eV, and the shoulder at about 2.7 eV. It is noteworthy that the blue shift of the maximum of the experimental spectrum from the monomer to the aggregate is reproduced with a remarkable accuracy (with an error of 0.02 eV considering the different shifts applied for the monomer and dimer spectra). The fact that the lowest energy band in the dimer is red-shifted with respect to the monomer is also reproduced, although the extent of the shift is overestimated, indicating that the splitting of the two computed bands is too large. The relative intensity of the lowest energy peak with respect to the maximum peak (i.e. the ratio *R* discussed in the Introduction) is on the contrary underestimated, and slightly more for the *syn* configurations (see top panel of Figure 6). It should be however noticed that the relative intensity of the first band is underestimated even for the monomer, and this inaccuracy was traced back⁴⁹ to the use of CAM-B3LYP, by showing that a standard hybrid functional would produce the opposite effect. We remind that the choice of a long-range corrected functional for the dimer is mandatory to get a proper description of the possible role of the inter-monomer CT states. Nonetheless, a further analysis of the possible causes of such underestimation is reported for the dimer in Section 2.9 in the Supporting Information.

## 4.5 An analysis of the factors determining the spectral shape

The separate vibronic spectra computed for each snapshot sampled along the *anti* and *syn* trajectories in ACN are shown as thin gray lines in the left panel of Figure 7, while the same spectra are reported for water in Figure S29 in the Supporting Information. The total spectrum, displayed in Figure 7 and evidenced through thick lines, clearly results in a broadening and smoothing of the different peaks exhibited by the
<figure><img src="image_9.png" /><figcaption>**Figure 7:** Left panels: Ad-MD|gLVC averaged (red and green lines) spectrum and individual (gray lines) spectra for the *anti* (top) and *syn* (bottom) MD trajectories of the PDI dimer. The insets show the averaged spectrum computed with or without accounting for the $E_{LE-CT}(0)$ coupling. All calculated spectra are red-shifted by 0.27 eV. Spectra convoluted with a Gaussian of HWHM = 0.03 eV. Right panels: Distribution of the $E_{LE-LE}(0)$ coupling along the *anti* (red) and *syn* (green) MD trajectories. The sampling interval employed for $E_{LE-LE}(0)$ is the same for both sets and correspond to the one shown by the red bins; *syn* data are reported with narrower bins only for clarity.</figcaption></figure>

many individual signals, which, however, are not wiped out by averaging, in agreement with what is observed in experiment. Actually, the computed spectra are even slightly too wide with respect to the experiment and this may be attributed to the broadening gaussian functions employed. In fact, despite formally speaking the Ad-MD|gLVC method already includes all the main possible sources of broadening and therefore a phenomenological function would not be required, computational feasibility allows for using only a small number of effective coordinates, and a small broadening with an HWHM of 0.03 eV was applied to overcome this issue.

In order to investigate the effect of the coupling between the bright LE and the dark CT states, in the insets of the left panels of Figure 7 we compare the average spectra computed including or neglecting ($E_{LE-CT}(0) = 0$) the LE-CT contribution. Notably, the differences are quite modest, indicating that CT states have a marginal effect on the shape of the absorption spectra. This is not simply an averaging out effect, as the spectra of the individual snapshots also show only a moderate CT contribution. It is interesting that for this solvated dimer the CT states do not play a large role in the spectral shape, whilst for perylene
aggregates they have been shown to have significant effects.¹³,¹⁴ Of further interest is that this does not mean that the CT states are not populated after photoexcitation, as we will show in the following section. The exciton coupling $E_{LE-LE}$, on the contrary, has a large impact on the spectral shape. The right panel of Figure 7 shows that it ranges from 0.06 to 0.16 eV for the different snapshots of the *anti* trajectory, although most of them (about 75%) exhibit a coupling between 0.10-0.14 eV. The effect of such a coupling on the spectral shape was further investigated for the *anti* trajectory by collecting, in the different panels of Figure 8, the individual and averaged spectra corresponding to exciton couplings falling into the same bins considered in Figure 7. It is apparent that as the exciton coupling increases the lowest-energy band becomes weaker with respect to the maximum band, in agreement with the results reported by Spano, showing that, for H-aggregates, the ratio *R* between the first two vibronic peaks decreases by increasing excitonic coupling.¹²,¹⁸ In addition, the splitting between the two bands also increases with the coupling: for snapshots with $E_{LE-LE} \sim 0.08$ eV the gap between both bands is 0.2 eV while for $E_{LE-LE} > 0.14$ eV it is 0.3 eV.

The different ratio *R* found in the simulated spectrum of the dimer with respect to the one experimentally obtained for a PDI aggregate can therefore be due, at least partially, to the value of the exciton coupling. On the one hand, as reported in Table C of the Supporting Information, the exciton-coupling does not change significantly either upon variation of the functional or if it is computed from the Coulomb coupling of the transition densities⁸⁸,⁸⁹ rather than with the present diabatization scheme, thus confirming the robustness of our estimates. On the other hand, the present results strongly point to a connection between the extent of the exciton coupling and the spectral shape. For instance, according to Figure 7, the distribution of such couplings over the *syn* trajectory is narrower than for the *anti* one, and its tallest bin is shifted at larger values, in agreement with the more aligned dimer arrangements visited by the dynamics (see Figure 4). These findings explain why the ratio *R* is lower for *syn* (larger couplings) and why the fluctuation of the spectral shapes of the different snapshots is smaller (narrower distribution).

A more in-depth microscopic understanding can be achieved by correlating the magnitude of $E_{ij}$ terms ($i \neq j$) with the variation, along the *syn* and *anti* MD trajectories, of the intermolecular descriptors defined in Figure 3a and S1. Here we focus on both the LE-LE and LE-CT couplings, to separately analyze the different role respectively played on the spectral shape. Nonetheless, since LE-CT couplings were found to affect such shape only slightly, we display only the coupling between the $|L_1\rangle$ and $|CT(1 \to 2)\rangle$ states, which for the majority of the snapshots is the largest LE-CT coupling. As evidenced in Figure 9 and Figures S30 to S41 in the Supporting Information, both $E_{LE-LE}$ and $E_{LE-CT}$ are larger for the *syn* aggregate, revealing that the diabatic couplings are generally larger for configurations showing the two PDIs in a more cofacial
<figure><img src="image_10.png" /><figcaption>**Figure 8:** Calculated spectra of a set of snapshots with different values of the $E_{LE-CT}(0)$ coupling. PDI dimer with the lateral chains in *anti* conformation. All the calculated spectra were redshifted by 0.27 eV, and all stick transitions were convoluted with a Gaussian of HWHM = 0.03eV.</figcaption></figure>

This analysis also shows that the magnitude of the LE-CT couplings is more sensitive to the geometric variations of the aggregates observed along the MD trajectories compared to the LE-LE ones. This outcome agrees well with the considerations drawn for the parametric study reported in Section 4.3 for the unsubstituted PDI dimer. In terms of absolute value, the $E_{LE-CT}$ terms may reach higher values (up to ~ 0.26 eV for the *syn* aggregate) compared to the $E_{LE-LE}$, although for most of the snapshots, and therefore for most stable PDI dimer configurations, these terms are significantly smaller.

Up to now we have discussed the existing correlation between the extent of the exciton coupling and the
<figure><img src="image_11.png" /><figcaption>**Figure 9:** 2D heat maps, correlating the intensity of the LE-LE Coupling (top panels) and the LE-CT Coupling (bottom panels) with the values of ρ(Å) and α(degrees) assumed from the dimer. For this plot we consider the snapshots extracted from the *anti* (left panels) and *syn* (right panels) MD trajectories. The colour palette refers to the intensity of the interstate couplings (eV).</figcaption></figure>

ratio *R*. However Figure 6 shows the computed *R* is too large already in the monomer, and such a finding can only be due to an overestimation of the energy gradient of the LE state (i.e. the $\lambda_{LE-LE}$ vector). Looking for a qualitative understanding of the interplay between $\lambda_{LE-LE}$ and $E_{LE-LE}$ in determining the spectral shape, in Section 2.9 of the Supporting Information we examined the predictions of a minimal two-states two-coordinates model able to capture the main features of the spectrum. It includes the two local excitations $L_1$ and $L_2$ and two coordinates $q_1$ and $q_2$ (i.e. the first members of the hierarchy), which describe the displacement of the $L_1$ and $L_2$ equilibrium positions. In such a small model $\lambda_{LE-LE}$ becomes a scalar, $\lambda$. Figure S28 in the Supporting Information shows that *R* decreases at the increase of $E_{L_1,L_2}$. On the contrary, *R* increases at the increase of $\lambda$, as an indirect effect due to the transfer of intensity to the third band at 2.7 eV. The splitting of the two lowest bands increases with both $E_{L_1,L_2}$ and $\lambda$. While confirming that the slight disagreement on *R* with respect to the experiment is related to the values of the exciton coupling, this analysis also suggests that the excessive splitting of the two lowest-energy computed bands arises from
the combination of the high values of the exciton coupling and the overestimation of the displacement of
the equilibrium position of the LE state due to the employment of the CAM-B3LYP functional. Indeed, a
similar over-estimation of excited state gradients/reorganisation energies with the CAM-B3LYP functional
was also noted for a PDI dimer in Ref. [52].

## 4.6 Time evolution of the diabatic electronic populations

Although the main focus of the present paper is on the absorption spectrum of a PDI dimer in solution, our approach also gives direct access to truly QD time-dependent observables, like the electronic populations of the different diabatic states. This is a noteworthy result since the QD simulation of nonadiabatic dynamics in condensed phase is still an open challenge for current research.⁹⁰ In Figure 10 we report the time evolution of the diabatic populations after an instantaneous excitation of the first PDI (state $|L_1\rangle$), for each of the snapshots extracted from the *anti* trajectory and their average. The latter, displayed in Figure 10 with thick red lines, represents the final QD prediction for the photoexcited PDI dimer in ACN solution in the limiting situation in which the environment is considered to be frozen during the excited-state dynamics. Such a limiting case, sometimes known as “static-disorder”, has been considered by several authors in the recent past for different systems.⁹¹⁻⁹⁴ Indeed, here we assume that the soft degrees of freedom – which comprise all solvent coordinates as well as those describing the PDI cores relative arrangement and the flexible dihedrals of the alkyl chains – are so slow with respect to the electronic dynamics they can be considered as frozen in the configuration sampled by the MD run when the photoexcitation occurs. Then, as mentioned in the computational details, QD propagations were performed with the ML-MCTDH method including all the relevant fast coordinates (44).

Figure 10 shows that immediately after photoexcitation of one monomer a fast population transfers occur, and in a short time (~60 fs) the averaged populations (red lines) reach limiting values. The final populations are 35-40% for both the local excitations, and 10-15% for the CT states, with a slightly higher population for CT(1→2), i.e. the CT state with a hole on the PDI monomer which was initially excited. Interestingly, the population of the local excitations shows some quantum beats, more evident in the first 40 fs, due to the coupling with the vibrational motion: when the population of $|L_1\rangle$ is maximum, the one of $|L_2\rangle$ state is minimum, and vice versa. On the contrary, the CT population undergoes a nearly monotonic growth with some weak oscillations due to the vibrations. The differences in the dynamics observed in each snapshot, and individually displayed in Figure 10 with thin gray lines, highlight the sensitivity of the population exchange to the initial configuration of the slow degrees of freedom. In this respect, it is noteworthy that quantum beats survive even after the average has been taken, indicating that the coupling of electronic populations
<figure><img src="image_12.png" /><figcaption>**Figure 10:** Population dynamics of the diabatic states $|L_1\rangle$ (top left), $|L_2\rangle$ (bottom left), $|CT(1 \to 2)\rangle$ (top right) and $|CT(2 \to 1)\rangle$ (bottom right) after initial excitation on $|L_1\rangle$. The dynamics of individual snapshots are shown in gray and the average in red. The results for a specific case in which the populations of CT(1→2) and CT(2→1) are, respectively, large and vanishingly small are colored in blue. ML-MCTDH including 44 effective coordinates.</figcaption></figure>

and vibrational motion is quite strong.

Additionally, despite the fact that the average population for CT states is about 10%, for some specific
snapshots much larger populations are observed, up to 60% in a few fs. Indeed, the CT state populations
are much more sensitive to the specific snapshot than the LE populations. This finding can be rationalized
by the fact that the CT diabatic energies and LE-CT couplings exhibit a larger distribution than the other
energies and couplings (shown by the standard deviation of these values from the snapshots in Table D).
Focusing on the results of a specific snapshot, highlighted in blue in Figure 10, it is interesting to notice
that whereas $|CT(1 \to 2)\rangle$ state is steadily populated, $|CT(2 \to 1)\rangle$ is not. We verified that this behaviour
is systematic in other snapshots: when the $|CT(1 \to 2)\rangle$ is strongly favored, $|CT(2 \to 1)\rangle$ is not populated
at all and vice-versa. This finding can be explained by the fact that, as shown in Table E in the Supporting
Information, for the blue highlighted snapshot, $|CT(1 \to 2)\rangle$ is the most stable state (0.35 eV below $|L_1\rangle$ and
$|L_2\rangle$), whereas $|CT(2 \to 1)\rangle$ is remarkably higher in energy (+0.72 eV with respect to $|L_1\rangle$). Comparison to
gas phase results in the same Table indicates that the stabilization/destabilization of the CT states is caused
by the specific solvent configuration. This finding highlights that the solvent electrostatic field can be strong
enough to impart a preferential direction to the charge-transfer, suggesting that the hole/electron separation
can be very efficiently realized in PDI aggregates, in line with the experimental evidences reported at the
interface between a n-type semiconductor and a water oxidation catalyst.⁹

**5 Discussion and Conclusions**

The general MQC computational method presented in this work aims to simulate the nonadiabatic vibronic
spectra of molecular aggregates, coupling classical MD sampling and QD wavepacket propagations. The
protocol, named Ad-MD|gLVC provides a robust framework to account for the effect of the fluctuations of
both the aggregate and its embedding environment. The slow dynamics of the aggregate's supramolecular
structure, the flexibility of the alkyl chains and the fluctuations of explicit solvent molecules are accounted
for with a classical MD sampling, whereas the effect of the motion of the fast vibrations within PDI cores on
the coupled surfaces of local excitations and charge transfer states are included at QD level (thus preserving
vibronic resolution). Ad-MD|gLVC largely extends the capabilities of Ad-MD|gVH, a MQC approach recently
proposed by some of us for cases with negligible inter-state couplings,⁵⁰ and applied to the monomer of PDI
in ACN in ref. 49.

Here we apply this method to the simulation of the absorption spectrum of a PDI dimer in ACN and
in water. The results reproduce all the trends observed experimentally when comparing the spectra of the
PDI in monomeric and aggregate form. In particular the blue-shift of the maximum and the inversion
of the intensities of the first and second vibronic peaks are nicely reproduced by Ad-MD|gLVC, correctly
predicting that vibronic structure remains visible despite the fluctuations of the soft degrees of freedom and
the environment.

The slight underestimation of the ratio "R" of the intensities of first two vibronic peaks and the over-
estimation of their separation were traced back to the combination of two different factors: a too large
displacement after the local excitation (already observed in the case of the monomer, and attributed to the
use of the CAM-B3LYP functional⁴⁹) and an overestimation of the exciton coupling. The latter suggests that
in the dimer arrangements sampled by MD the two monomers are closer to each other with respect to the
experimental aggregate. Considering the agreement found between the adopted inter-molecular FF and the
reference CAM-B3LYP/D3 interaction energy curves, such closer distances can be more probably attributed
to the different number of PDI units in our model with respect to the experimental aggregate. It is likely in
fact that interactions with neighbours on both sides of a PDI monomer will lead to greater inter-monomer
spacing within the aggregate than in the dimer. In summary, our approach allows to accurately describe the
difference of the spectral shape between monomers and dimers even if, in absolute terms, some inaccuracy
on the relative intensities of the vibronic bands exist even for the monomer. DFT/TD-DFT usually offers a
good compromise between accuracy and computational cost. However, in order to adopt our approach as a
predictive tool, i.e. for species not yet synthesized, it might be beneficial to run a preliminary becnhmark
analysis on similar systems to select the best DFT functional. Alternatively, it is in principle possible to
couple our approach with different electronic-structure methods. This potentiality has already been shown
for pyrene, adopting our code *Overdia* to parameterize an LVC Hamiltonian with multiconfigurational meth-
ods with perturbative corrections, like RASPT2.<sup>95</sup> Extensions in combination with methods of the coupled
cluster family, like the promising domain-based local pair natural orbital similarity transformed equation of
motion-coupled cluster singles and doubles (DLPNO-STEOM-CCSD) are also possible.<sup>96,97</sup>

The Ad-MD|gLVC approach also allowed us to investigate the time evolution of the electronic populations
after the photoexcitation. To this end we assumed that the motion of the soft coordinates is much slower
than the nonadiabatic dynamics of interest, a limit sometimes described as "static disorder". Remarkably,
although our results indicate that the existence of CT states alter only slightly the shape of the absorption
spectrum of the dimer, they also point out that such states are partially populated in short time scales. More
importantly, we show that some specific arrangements of the adjacent PDIs, as well as of the solvent can
induce a marked directionality in the charge transfer, preferentially populating CT state 1→2 with respect
to 2→1, or *vice versa*. This result is interesting for the investigation of the opto-electronic properties of PDI
stacked aggregates which are often adopted in photovoltaics devices, where charge separation and transport
within the aggregates initiate the interfacial charge separation and collection.<sup>9</sup> At the same time, the fast
and effective population transfer between the two local excitations indicate that exciton migration/diffusion
might be also very effective in PDI aggregates.

It should be mentioned that most of the calculations reported in this work have been performed with
the simplified route of the Ad-MD|gLVC method described in Section 2.3, since it allows a remarkable
reduction of the computational time. Its reliability is based on the assumption that the fluctuation of the
dimer structure and solvent cavity mostly affects the vertical transition energies and the inter-state couplings,
whose dependence on the small oscillations of the fast coordinates can be, on the contrary, neglected. Whereas
for PDI₂ these approximations were shown to be accurate, it will be interesting in future work to investigate
if they can be challenged in different systems featuring, for instance, stronger inter-monomeric interactions
like hydrogen bonds.

Alternative routes to speed up the Ad-MD|gLVC calculations can also be envisaged. For instance,
instead of running a different QD simulation for the LVC Hamiltonian parameterized at each snapshot α, it is possible to run a single QD simulation driven by a LVC Hamiltonian averaged over all the snapshots. Of course the two procedures are not expected to be equivalent, and the former approach that we adopted in this manuscript is in principle more correct. We compare results reported in this manuscript with those predicted by the averaged Hamiltonian for both the absorption spectrum and the population dynamics in the last section of the Supporting Information. Very interestingly the spectrum compares nicely, even if not perfectly, although there is the need to adopt a larger phenomenological broadening to reintroduce the effect of the fluctuation of the LVC parameters that is lost in the average procedure. The comparison of the time-evolution of the diabatic populations is also good, although, not surprisingly, those predicted by a single averaged Hamiltonian exhibits larger quantum beats, which are smoothed away averaging the predictions of the snapshot-specific Hamiltonians. A significant difference can however be noted in that the CT population predicted by the average Hamiltonian is remarkably smaller. This is a clear demonstration that, as already pointed out by some of us in ref. 98, the outcome of a nonadiabatic dynamics is not linear with the fluctuation of the parameters of the Hamiltonian. In a situation in which CT states are generally less stable than LE ones, their further destabilization in some snapshots leads to a minor effect on the CT populations. On the contrary, specific snapshots for which CT states are stabilized predict a very large population transfer toward these states, making the CT population averaged over the different dynamics larger than the one predicted by the average Hamiltonian.

In summary, combining classical MD sampling and QD simulations, our fully atomistic and non-phenomenological approach has allowed for an in-depth analysis of the correlation between the instantaneous structure of the dimer, the parameters of the diabatic Hamiltonian and the predicted vibronic spectral shape and population dynamics, giving us the possibility to achieve a detailed microscopic understanding of the structural parameters determining the photophysics of the system. Among the other things, we found that (i) the exciton coupling is a key parameter for determining the spectral shape and it has a coupled dependence on both the stacking distance and the rotation α, being larger for cofacial arrangements of the dimer, and that the orientation of flexible pendants play a role in its modulation. On the contrary, the dependence of the exciton coupling on the twisting of the two monomeric planes is less relevant. Additionally, (ii) we discovered that CT states are not important for the spectrum, although they gain non-negligible populations in the 100 fs time scale, and that such population would be underestimated without a careful description of the fluctuation of the dimer structure. It is worthy to highlight that these aspects might have been missed by more reductionist approaches which compute the effect of one factor per time without considering their interplay. Finally, we determined that (iii) although an effective mode for each monomer can qualitatively explain the
shape of the spectra, a minimal, non-phenomenological model to accurately reproduce the vibronic spectral
shape determined by all the fast coordinates is a chain of 6 sequentially coupled effective modes on each
monomer.

The Ad-MD|gLVC method may be compared to similar recent coupled quantum-classical methodologies,
which aim to reproduce absorption spectra and analyse photoexcited dynamics of perylenes in solution.<sup>35,52</sup>
Firstly, the use of QMD-FFs in the MD simulations allows to accurately consider the effect of the flexible side
chains, a factor neglected in Ref. 52, where the authors assumed that the side groups have little influence
on the relative geometries of the two perylene units. On the contrary, the present results indicate that
some modulation of the final spectra can be connected with the conformational dynamics of the side chains,
which does give rise to different dimer arrangements. This should be noted with the caveat that in Ref. 52
the perylenes were linked by a bridging phosphate group, thus limiting the conformational space sampled
by the side chains with respect to two untethered PDI units. Next, our FrD methodology allows us to
straightforwardly include both LE and CT states, whilst in Ref. 52 only LE were considered. However,
both approaches can account for mutual polarisation effects of closely spaced chromophores on the LE
energies and couplings. Conversely, the effect of CT states was included in Ref. 35, but directly from the
adiabatic TD-DFT states of the dimer, rather than from a diabatisation scheme as in the present work. Using
rigorously defined LE and CT states from a diabatization scheme can in principle give a more straightforward
interpretation of the photoexcited dynamics when adiabatic states have strongly mixed LE/CT character.
Finally, both ourselves and the authors of Ref. 52 include nonadiabatic effects on the spectra, whilst in
Ref. 35 the authors do not, as they calculate the spectra via single point energy computations on MD
snapshots, with the nonadiabatic dynamics investigated separately using a surface hopping approach. In Ref.
52 the nonadiabatic effects on the spectra are included semi-classically and with a density matrix formalism
through the partially linearised density matrix (PLDM) method,<sup>99</sup> whereas we utilise a quantum wavepacket
treatment through the LVC model and MCTDH method. The PLDM method offers the advantage of
more straightforward inclusion of temperature and dissipative effects, whilst the MCTDH method offers the
advantage of numerical exactness, as well as wavepacket approaches typically being more numerically efficient
than density matrix ones.<sup>100</sup> Recent ML-MCTDH calculations on PDI aggregates have included temperature
effects through the thermofield dynamics approach however,<sup>101,102</sup> and this represents a promising future
avenue for comparison of temperature effects on vibronic spectra via wavepacket versus density matrix based
methods.

In conclusion, beside the specific application presented here, we believe that the proposed Ad-MD|gLVC
model will be a powerful tool in the close future to study the optoelectronic and photophysical properties of
diverse molecular aggregates. In fact, although here we considered a dimer, application at longer oligomers is possible, thanks to the effectiveness of the hierarchical representation of the LVC Hamiltonians, and the simplifications of the full protocol introduced in the present paper.

## Supporting Information

Theory and further computational details, additional results on model systems, experimental absorption and difference spectra, supplementary 2D heat maps.

## Acknowledgements

The authors thank Javier Cerezo, Universidad Autónoma de Madrid, for the many helpful discussions on the implementation of Ad-MD|gLVC method, and Marcel Nooijen, University of Waterloo for suggesting to compare the obtained results with those predicted by a simple average Hamiltonian. DA acknowledges Fundación Ramón Areces (Spain) and Generalitat Valenciana/European Social Fund (APOSTD/2021/025) for fundings. M.P and A.S. thank the financial support from the COMETE project (COception in silico des Matériaux pour l'Environnement et l'Energie) co-funded by the European Union under the program 'FEDER-FSE Lorraine et Massif des Vosges 2014-2020.

## References

(1) Würthner, F. Perylene bisimide dyes as versatile building blocks for functional supramolecular architectures. *Chem. Commun.* **2004**, 1564–1579.

(2) Würthner, F.; Saha-Möller, C. R.; Fimmel, B.; Ogi, S.; Leowanawat, P.; Schmidt, D. Perylene Bisimide Dye Assemblies as Archetype Functional Supramolecular Materials. *Chemical Reviews* **2016**, *116*, 962–1052.

(3) Zhang, G.; Zhao, J.; Chow, P. C. Y.; Jiang, K.; Zhang, J.; Zhu, Z.; Zhang, J.; Huang, F.; Yan, H. Non-fullerene Acceptor Molecules for Bulk Heterojunction Organic Solar Cells. *Chemical Reviews* **2018**, *118*, 3447–3507.

(4) Yan, C.; Barlow, S.; Wang, Z.; Yan, H.; Jen, A. K.-Y.; Marder, S. R.; Zhan, X. Non-fullerene acceptors for organic solar cells. *Nature Reviews Materials* **2018**, *3*, 18003.
(5) Yang, Z.; Chen, X. Semiconducting Perylene Diimide Nanostructure: Multifunctional Phototheranostic Nanoplatform. *Acc. Chem. Res.* **2019**, *52*, 1245–1254.

(6) Chen, Z.; Fimmel, B.; Würthner, F. Solvent and substituent effects on aggregation constants of perylene bisimide π-stacks – a linear free energy relationship analysis. *Org. Biomol. Chem.* **2012**, *10*, 5845–5855.

(7) Segalina, A.; Assfeld, X.; Monari, A.; Pastore, M. Computational Modeling of Exciton Localization in Self-Assembled Perylene Helices: Effects of Thermal Motion and Aggregate Size. *J. Phys. Chem. C* **2019**, *123*, 6427–6437.

(8) Grande, V.; Soberats, B.; Herbst, S.; Stepanenko, V.; Würthner, F. Hydrogen-bonded perylene bisimide J-aggregate aqua material. *Chem. Sci.* **2018**, *9*, 6904–6911.

(9) Ronconi, F.; Syrgiannis, Z.; Bonasera, A.; Prato, M.; Argazzi, R.; Caramori, S.; Cristino, V.; Bignozzi, C. A. Modification of Nanocrystalline WO3 with a Dicationic Perylene Bisimide: Applications to Molecular Level Solar Water Splitting. *J. Am. Chem. Soc.* **2015**, *137*, 4630–4633.

(10) Berardi, S. et al. Perylene Diimide Aggregates on Sb-Doped SnO₂: Charge Transfer Dynamics Relevant to Solar Fuel Generation. *J. Phys. Chem. C* **2017**, *121*, 17737–17745.

(11) Brixner, T.; Hildner, R.; Köhler, J.; Lambert, C.; Würthner, F. Exciton Transport in Molecular Aggregates – From Natural Antennas to Synthetic Chromophore Systems. *Advanced Energy Materials* **2017**, *7*, 1700236.

(12) Hestand, N. J.; Spano, F. C. Expanded theory of H-and J-molecular aggregates: the effects of vibronic coupling and intermolecular charge transfer. *Chem. Rev.* **2018**, *118*, 7069–7163.

(13) Hestand, N. J.; Spano, F. C. Molecular aggregate photophysics beyond the Kasha model: novel design principles for organic materials. *Acc. Chem. Res.* **2017**, *50*, 341–350.

(14) Hestand, N. J.; Spano, F. C. Interference between Coulombic and CT-mediated couplings in molecular aggregates: H- to J-aggregate transformation in perylene-based π-stacks. *The Journal of Chemical Physics* **2015**, *143*, 244707.

(15) Kasha, M. Energy transfer mechanisms and the molecular exciton model for molecular aggregates. *Rad. Res.* **1963**, *20*, 55–70.
(16) Kilina, S.; Kilin, D.; Tretiak, S. Light-Driven and Phonon-Assisted Dynamics in Organic and Semi-
conductor Nanostructures. *Chem. Rev.* **2015**, *115*, 5929–5978.

(17) De Sio, A.; Lienau, C. Vibronic coupling in organic semiconductors for photovoltaics. *Phys. Chem. Chem. Phys.* **2017**, *19*, 18813–18830.

(18) Spano, F. C. The spectral signatures of Frenkel polarons in H-and J-aggregates. *Acc. Chem. Res.* **2010**, *43*, 429–439.

(19) Fulton, R. L.; Gouterman, M. Vibronic coupling. I. Mathematical treatment for two electronic states.
*J. Chem. Phys.* **1961**, *35*, 1059–1071.

(20) Fulton, R. L.; Gouterman, M. Vibronic coupling. II. Spectra of dimers. *J. Chem. Phys.* **1964**, *41*, 2280–2286.

(21) Hemenger, R. A theory of optical absorption by aggregates of large molecules. *J. Chem. Phys.* **1977**, *66*, 1795–1801.

(22) Knapp, E.; Scherer, P.; Fischer, S. On the lineshapes of vibronically resolved molecular aggregate
spectra. application to pseudoisocyanin (PIC). *Chem. Phys. Lett.* **1984**, *111*, 481–486.

(23) Pochas, C. M.; Kistler, K. A.; Yamagata, H.; Matsika, S.; Spano, F. C. Contrasting photophysical
properties of star-shaped vs linear perylene diimide complexes. *J. Am. Chem. Soc.* **2013**, *135*, 3056–
3066.

(24) Eisfeld, A.; Briggs, J. The J-and H-bands of organic dye aggregates. *Chem. Phys.* **2006**, *324*, 376–384.

(25) Chen, Z.; Stepanenko, V.; Dehm, V.; Prins, P.; Siebbeles, L. D.; Seibt, J.; Marquetand, P.; Engel, V.; Würthner, F. Photoluminescence and conductivity of self-assembled π-π stacks of perylene bisimide dyes. *Chem. Eur. J.* **2007**, *13*, 436–449.

(26) Schröter, M.; Ivanov, S. D.; Schulze, J.; Polyutov, S. P.; Yan, Y.; Pullerits, T.; Kühn, O. Exciton-vibrational coupling in the dynamics and spectroscopy of Frenkel excitons in molecular aggregates. *Phys. Rep.* **2015**, *567*, 1–78.

(27) Seibt, J.; Marquetand, P.; Engel, V.; Chen, Z.; Dehm, V.; Würthner, F. On the geometry dependence of molecular dimer spectra with an application to aggregates of perylene bisimide. *Chem. Phys.* **2006**, *328*, 354–362.
(28) Seibt, J.; Winkler, T.; Renziehausen, K.; Dehm, V.; Würthner, F.; Meyer, H.-D.; Engel, V. Vibronic transitions and quantum dynamics in molecular oligomers: A theoretical analysis with an application to aggregates of perylene bisimides. *J. Phys. Chem. A* **2009**, *113*, 13475–13482.

(29) Spano, F. C.; Yamagata, H. Vibronic coupling in J-aggregates and beyond: a direct means of determining the exciton coherence length from the photoluminescence spectrum. *J. Phys. Chem. B* **2011**, *115*, 5133–5143.

(30) Ambrosek, D.; Kohn, A.; Schulze, J.; Kühn, O. Quantum chemical parametrization and spectroscopic characterization of the frenkel exciton hamiltonian for a J-aggregate forming perylene bisimide dye. *J. Phys. Chem. A* **2012**, *116*, 11451–11458.

(31) Ritschel, G.; Suess, D.; Möbius, S.; Strunz, W. T.; Eisfeld, A. Non-Markovian Quantum State Diffusion for temperature-dependent linear spectra of light harvesting aggregates. *J. Chem. Phys.* **2015**, *142*, 034115.

(32) Clark, A. E.; Qin, C.; Li, A. D. Q. Beyond Exciton Theory: A Time-Dependent DFT and Franck-Condon Study of Perylene Diimide and Its Chromophoric Dimer. *J. Am. Chem. Soc.* **2007**, *129*, 7586–7595.

(33) Oltean, M.; Calborean, A.; Mile, G.; Vidrighin, M.; Iosin, M.; Leopold, L.; Maniu, D.; Leopold, N.; Chiş, V. Absorption spectra of PTCDI: A combined UV–Vis and TD-DFT study. *Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy* **2012**, *97*, 703–710.

(34) Pan, F.; Gao, F.; Liang, W.; Zhao, Y. Nature of Low-Lying Excited States in H-Aggregated Perylene Bisimide Dyes: Results of TD-LRC-DFT and the Mixed Exciton Model. *J. Phys. Chem. B* **2009**, *113*, 14581–14587.

(35) Mukazhanova, A.; Trerayapiwat, K. J.; Mazaheripour, A.; Wardrip, A. G.; Frey, N. C.; Nguyen, H.; Gorodetsky, A. A.; Sharifzadeh, S. Accurate First-Principles Calculation of the Vibronic Spectrum of Stacked Perylene Tetracarboxylic Acid Diimides. *J. Phys. Chem. A* **2020**, *124*, 3055–3063.

(36) Padula, D.; Santoro, F.; Pescitelli, G. A simple dimeric model accounts for the vibronic ECD spectra of chiral polythiophenes in their aggregated states. *RSC Adv.* **2016**, *6*, 37938–37943.

(37) Kistler, K. A.; Pochas, C. M.; Yamagata, H.; Matsika, S.; Spano, F. C. Absorption, Circular Dichroism, and Photoluminescence in Perylene Diimide Bichromophores: Polarization-Dependent H- and J-Aggregate Behavior. *The Journal of Physical Chemistry B* **2012**, *116*, 77–86.
(38) Clark, A. E.; Qin, C.; Li, A. D. Beyond exciton theory: a time-dependent DFT and Franck-Condon study of perylene diimide and its chromophoric dimer. *J. Am. Chem. Soc.* **2007**, *129*, 7586–7595.

(39) Houari, Y.; Laurent, A. D.; Jacquemin, D. Spectral signatures of perylene diimide derivatives: insights from theory. *J. Phys. Chem. C* **2013**, *117*, 21682–21691.

(40) Gao, F.; Zhao, Y.; Liang, W. Vibronic Spectra of Perylene Bisimide Oligomers: Effects of Intermolecular Charge-Transfer Excitation and Conformational Flexibility. *The Journal of Physical Chemistry B* **2011**, *115*, 2699–2708.

(41) Gisslén, L.; Scholz, R. Crystallochromy of perylene pigments: Interference between Frenkel excitons and charge-transfer states. *Phys. Rev. B* **2009**, *80*, 115309.

(42) Gisslén, L.; Scholz, R. Crystallochromy of perylene pigments: Influence of an enlarged polyaromatic core region. *Phys. Rev. B* **2011**, *83*, 155311.

(43) Austin, A.; Hestand, N. J.; McKendry, I. G.; Zhong, C.; Zhu, X.; Zdilla, M. J.; Spano, F. C.; Szarko, J. M. Enhanced Davydov Splitting in Crystals of a Perylene Diimide Derivative. *The Journal of Physical Chemistry Letters* **2017**, *8*, 1118–1123.

(44) Sung, J.; Kim, P.; Fimmel, B.; Würthner, F.; Kim, D. Direct observation of ultrafast coherent exciton dynamics in helical π-stacks of self-assembled perylene bisimides. *Nature Comm.* **2015**, *6*.

(45) Lim, J. M.; Kim, P.; Yoon, M.-C.; Sung, J.; Dehm, V.; Chen, Z.; Würthner, F.; Kim, D. Exciton delocalization and dynamics in helical π-stacks of self-assembled perylene bisimides. *Chem. Sci.* **2013**, *4*, 388–397.

(46) Neuteboom, E. E.; Meskers, S. C. J.; Meijer, E. W.; Janssen, R. A. J. Photoluminescence of Self-organized Perylene Bisimide Polymers. *Macromolecular Chemistry and Physics* **2004**, *205*, 217–222.

(47) Zuehlsdorff, T. J.; Isborn, C. M. Combining the ensemble and Franck-Condon approaches for calculating spectral shapes of molecules in solution. *The Journal of Chemical Physics* **2018**, *148*, 024110.

(48) Zuehlsdorff, T. J.; Napoli, J. A.; Milanese, J. M.; Markland, T. E.; Isborn, C. M. Unraveling electronic absorption spectra using nuclear quantum effects: Photoactive yellow protein and green fluorescent protein chromophores in water. *J. Chem. Phys.* **2018**, *149*, 024107.
(49) Segalina, A.; Cerezo, J.; Prampolini, G.; Santoro, F.; Pastore, M. Accounting for Vibronic Features through a Mixed Quantum-Classical Scheme: Structure, Dynamics, and Absorption Spectra of a Perylene Diimide Dye in Solution. *Journal of Chemical Theory and Computation* **2020**, *16*, 7061–7077.

(50) Cerezo, J.; Aranda, D.; Avila Ferrer, F. J.; Prampolini, G.; Santoro, F. Adiabatic-Molecular Dynamics Generalized Vertical Hessian Approach: A Mixed Quantum Classical Method to Compute Electronic Spectra of Flexible Molecules in the Condensed Phase. *J. Chem. Theory Comput.* **2020**, *16*, 1215–1231.

(51) Bartlett, A.; Markegard, C. B.; Dibble, D. J.; Gorodetsky, A. A.; Sharifzadeh, S.; Nguyen, H. D. Molecular dynamics simulations of DNA-inspired macromolecules from perylenediimide base surrogates. *Synthetic Metals* **2019**, *253*, 146–152.

(52) Kumar, M.; Provazza, J.; Coker, D. F. Influence of solution phase environmental heterogeneity and fluctuations on vibronic spectra: Perylene diimide molecular chromophore complexes in solution. *The Journal of Chemical Physics* **2021**, *154*, 224109.

(53) Cacelli, I.; Prampolini, G. Parametrization and Validation of Intramolecular Force Fields Derived from DFT Calculations. *J. Chem. Theory Comput.* **2007**, *3*, 1803–1817.

(54) Barone, V.; Cacelli, I.; De Mitri, N.; Licari, D.; Monti, S.; Prampolini, G. Joyce and Ulysses: Integrated and User-Friendly Tools for the Parameterization of Intramolecular Force Fields from Quantum Mechanical Data. *Phys. Chem. Chem. Phys.* **2013**, *15*, 3736–51.

(55) Cerezo, J.; Prampolini, G.; Cacelli, I. Developing accurate intramolecular force fields for conjugated systems through explicit coupling terms. *Theor. Chem. Accounts* **2018**, *137*, 80.

(56) Meyer, F., H.-D. and Gatti, Worth, G. A., Eds. *Multidimensional Quantum Dynamics: MCTDH Theory and Applications*; Wiley-VCH: Weinheim, 2009.

(57) Beck, M. H.; Jäckle, A.; Worth, G. A.; Meyer, H.-D. The multiconfiguration time-dependent Hartree method: A highly efficient algorithm for propagating wavepackets. *Phys. Rep.* **2000**, *324*, 1–105.

(58) Wang, H.; Thoss, M. Multilayer formulation of the multiconfiguration time-dependent Hartree theory. *J. Chem. Phys.* **2003**, *119*, 1289–1299.

(59) Köppel, H.; Domcke, W.; Cederbaum, L. S. *Advances in Chemical Physics*; John Wiley & Sons, Ltd, 1984; pp 59–246.
(60) Green, J. A.; Asha, H.; Santoro, F.; Improta, R. Excitonic Model for Strongly Coupled Multichro-mophoric Systems: The Electronic Circular Dichroism Spectra of Guanine Quadruplexes as Test Cases. *Journal of Chemical Theory and Computation* **2021**, 17, 405-415.

(61) Avila Ferrer, F. J.; Santoro, F. Comparison of vertical and adiabatic harmonic approaches for the calculation of the vibrational structure of electronic spectra. *Phys. Chem. Chem. Phys.* **2012**, 14, 13549-13563.

(62) Köppel, H.; Domcke, W.; Cederbaum, L. The Multi-mode vibronic-coupling approach. Conical Intersections, Electronic Structure, Dynamics and Spectroscopy. 2004; pp 323-368.

(63) Cacelli, I.; Cerezo, J.; Soto, J.; Mitri, N.; Prampolini, G. JOYCE2.10, a Fortran 77 code for intra-molecular force field parameterization. , available free of charge at http://www pi icom cnr it/joyce, last consulted July. 2020.

(64) Cerezo, J.; Prampolini, G.; Cacelli, I. Developing accurate intramolecular force fields for conjugated systems through explicit coupling terms. *Theor. Chem. Acc.* **2018**, 137, 80.

(65) Avila Ferrer, F. J.; Cerezo, J.; Soto, J.; Improta, R.; Santoro, F. First-principle computation of absorption and fluorescence spectra in solution accounting for vibronic structure, temperature effects and solvent inhomogenous broadening. *Comput. Theoret. Chem.* **2014**, 1040-1041, 328-337.

(66) Lami, A.; Santoro, F. In *Computational Strategies for Spectroscopy*; Barone, V., Ed.; John Wiley & Sons, Inc., 2011; Chapter 10, pp 475-516.

(67) Aranda, D.; Santoro, F. Vibronic Spectra of π-Conjugated Systems with a Multitude of Coupled States: A Protocol Based on Linear Vibronic Coupling Models and Quantum Dynamics Tested on Hexahelicene. *Journal of Chemical Theory and Computation* **2021**, 17, 1691-1700, PMID: 33606542.

(68) Liu, Y.; Cerezo, J.; Mazzeo, G.; Lin, N.; Zhao, X.; Longhi, G.; Abbate, S.; Santoro, F. Vibronic Coupling Explains the Different Shape of Electronic Circular Dichroism and of Circularly Polarized Luminescence Spectra of Hexahelicenes. *Journal of Chemical Theory and Computation* **2016**, 12, 2799-2819, PMID: 27120334.

(69) Cederbaum, L. S.; Gindensperger, E.; Burghardt, I. Short-Time Dynamics Through Conical Intersections in Macrosystems. *Phys. Rev. Lett.* **2005**, 94, 113003.
(70) Gindensperger, E.; Burghardt, I.; Cederbaum, L. S. Short-time dynamics through conical intersections in macrosystems. I. Theory: Effective-mode formulation. *The Journal of Chemical Physics* **2006**, *124*, 144103.

(71) Gindensperger, E.; Burghardt, I.; Cederbaum, L. S. Short-time dynamics through conical intersections in macrosystems. II. Applications. *The Journal of Chemical Physics* **2006**, *124*, 144104.

(72) Picconi, D.; Lami, A.; Santoro, F. Hierarchical transformation of Hamiltonians with linear and quadratic couplings for nonadiabatic quantum dynamics: Application to the $\pi\pi^*/n\pi^*$ internal conversion in thymine. *The Journal of Chemical Physics* **2012**, *136*, 244104.

(73) Jorgensen, W. L.; Maxwell, D. S.; Tirado-rives, J. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. *J. Am. Chem. Soc.* **1996**, *7863*, 11225-11236.

(74) Jorgensen, W. L.; Tirado-Rives, J. Potential Energy Functions for Atomic-Level Simulations of Water and Organic and Biomolecular Systems. *Proc. Natl. Acad. Sci. USA* **2005**, *102*, 6665-70.

(75) Hess, B.; Bekker, B.; Berendsen, H.; J.G.E.M., F. LINCS: A linear constraint solver for molecular simulations. *J. Comp. Chem.* **1997**, *18*, 1463-1472.

(76) Frisch, M. J. et al. Gaussian 16 Revision B.01. 2016; Gaussian Inc. Wallingford CT.

(77) Yanai, T.; Tew, D. P.; Handy, N. C. A new hybrid exchange-correlation functional using the Coulomb-attenuating method (CAM-B3LYP). *Chemical physics letters* **2004**, *393*, 51-57.

(78) Grimme, S.; Ehrlick, S.; Goerigk, L. Effect of the Damping Function in Dispersion Corrected Density Functional Theory. *J. Comp. Chem.* **2011**, *32*, 1456-65.

(79) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. *The Journal of chemical physics* **2010**, *132*, 154104.

(80) Santoro, F.; Green, A. J. OVERDIA1.0, a Fortran 90 code for parametrization of model Hamiltonians based on a maximum-overlap diabatisation, available free of charge at http://www.iccom.cnr.it/overdia-en, last consulted March. 2022.
(81) Yaghoubi Jouybari, M.; Liu, Y.; Improta, R.; Santoro, F. The Ultrafast Dynamics of the Two Lowest Bright Excited States of Cytosine and 1-Methyl-Cytosine: A Quantum Dynamical Study. *J. Chem. Theory Comput.* **2020**, *16*, 5792–5808.

(82) Worth, G. A.; Giri, K.; Richings, G. W.; Beck, M. H.; Jäckle, A.; Meyer, H.-D. The QUANTICS Package, Version 1.1, (2015), University of Birmingham, Birmingham, U.K.

(83) Worth, G. Quantics: A general purpose package for Quantum molecular dynamics simulations. *Computer Physics Communications* **2020**, *248*, 107040.

(84) Padula, D.; Picconi, D.; Lami, A.; Pescitelli, G.; Santoro, F. Electronic circular dichroism in exciton-coupled dimers: vibronic spectra from a general all-coordinates quantum-dynamical approach. *The Journal of Physical Chemistry A* **2013**, *117*, 3355–3368.

(85) Wang, L.-P.; Song, C. Geometry optimization made simple with translation and rotation coordinates. *J. Chem. Phys.* **2016**, *144*, 214108.

(86) Prampolini, G.; Livotto, P. R.; Cacelli, I. Accuracy of Quantum Mechanically Derived Force-Fields Parameterized from Dispersion-Corrected DFT Data: The Benzene Dimer as a Prototype for Aromatic Interactions. *J. Chem. Theor. and Comput.* **2015**, *11*, 5182–96.

(87) Ghosh, S.; Li, X.-Q.; Stepanenko, V.; Würthner, F. Control of H- and J-Type $\pi$ Stacking by Peripheral Alkyl Chains and Self-Sorting Phenomena in Perylene Bisimide Homo- and Heteroaggregates. *Chemistry - A European Journal* **2008**, *14*, 11343–11357.

(88) Jurinovich, S.; Cupellini, L.; Guido, C. A.; Mennucci, B. EXAT: EXcitonic analysis tool. 2018.

(89) Iozzi, M. F.; Mennucci, B.; Tomasi, J.; Cammi, R. Excitation energy transfer (EET) between molecules in condensed matter: A novel application of the polarizable continuum model (PCM). *The Journal of Chemical Physics* **2004**, *120*, 7029–7040.

(90) Santoro, F.; Green, J. A.; Martinez-Fernandez, L.; Cerezo, J.; Improta, R. Quantum and semiclassical dynamical studies of nonadiabatic processes in solution: achievements and perspectives. *Phys. Chem. Chem. Phys.* **2021**, *23*, 8181–8199.

(91) Cerezo, J.; Aranda, D.; Avila Ferrer, F. J.; Prampolini, G.; Mazzeo, G.; Longhi, G.; Abbate, S.; Santoro, F. Toward a general mixed quantum/classical method for the calculation of the vibronic
ECD of a flexible dye molecule with different stable conformers: Revisiting the case of 2,2,2-trifluoro-
anthrylethanol. *Chirality* **30**, 730–743.

(92) Reiter, S.; Keefer, D.; de Vivie-Riedle, R. RNA Environment Is Responsible for Decreased Photostability of Uracil. *J. Am. Chem. Soc.* **2018**, *140*, 8714–8720, PMID: 29943578.

(93) Zauleck, J. P. P.; Peschel, M. T.; Rott, F.; Thallmair, S.; de Vivie-Riedle, R. Ultrafast Reactive Quantum Dynamics Coupled to Classical Solvent Dynamics Using an Ehrenfest Approach. *J. Phys. Chem. A* **2018**, *122*, 2849–2857, PMID: 29498853.

(94) Pápai, M.; Abedi, M.; Levi, G.; Biasin, E.; Nielsen, M. M.; Møller, K. B. Theoretical Evidence of Solvent-Mediated Excited-State Dynamics in a Functionalized Iron Sensitizer. *J. Phys. Chem. C* **2019**, *123*, 2056–2065.

(95) Aleotti, F.; Aranda, D.; Yaghoubi Jouybari, M.; Garavelli, M.; Nenov, A.; Santoro, F. Parameterization of a linear vibronic coupling model with multiconfigurational electronic structure methods to study the quantum dynamics of photoexcited pyrene. *The Journal of Chemical Physics* **2021**, *154*, 104106.

(96) Dutta, A. K.; Nooijen, M.; Neese, F.; Izsák, R. Exploring the Accuracy of a Low Scaling Similarity Transformed Equation of Motion Method for Vertical Excitation Energies. *Journal of Chemical Theory and Computation* **2018**, *14*, 72–91, PMID: 29206453.

(97) Berraud-Pache, R.; Neese, F.; Bistoni, G.; Izsák, R. Unveiling the Photophysical Properties of Boron-dipyrromethene Dyes Using a New Accurate Excited State Coupled Cluster Method. *Journal of Chemical Theory and Computation* **2020**, *16*, 564–575, PMID: 31765141.

(98) Cerezo, J.; Liu, Y.; Lin, N.; Zhao, X.; Improta, R.; Santoro, F. Mixed Quantum/Classical Method for Nonadiabatic Quantum Dynamics in Explicit Solvent Models: The $\pi\pi^*/n\pi^*$ Decay of Thymine in Water as a Test Case. *J. Chem. Theory Comput.* **2018**, *14*, 820–832, PMID: 29207245.

(99) Huo, P.; Coker, D. F. Communication: Partial linearized density matrix dynamics for dissipative, non-adiabatic quantum evolution. *The Journal of Chemical Physics* **2011**, *135*, 201101.

(100) Huo, P.; Coker, D. F. Consistent schemes for non-adiabatic dynamics derived from partial linearized density matrix propagation. *J. Chem. Phys.* **2012**, *137*, 22A535.
(101) Brey, D.; Popp, W.; Budakoti, P.; D'Avino, G.; Burghardt, I. Quantum Dynamics of Electron-Hole Separation in Stacked Perylene Diimide-Based Self-Assembled Nanostructures. *The Journal of Physical Chemistry C* **2021**, *125*, 25030-25043.

(102) Takahashi, Y.; Umezawa, H. Thermo Field Dynamics. *International Journal of Modern Physics B* **1996**, *10*, 1755-1805.