Quantum-Classical Protocol for Efficient Characterization of Absorption Lineshape and Fluorescence Quenching upon Aggregation: The Case of Zinc Phthalocyanine Dyes

A quantum-classical protocol that incorporates Jahn–Teller vibronic coupling effects and cluster analysis of molecular dynamics simulations is reported, providing a tool for simulations of absorption spectra and ultrafast nonadiabatic dynamics in large molecular photosystems undergoing aggregation in solution. Employing zinc phthalocyanine dyes as target systems, we demonstrated that the proposed protocol provided fundamental information on vibronic, electronic couplings and thermal dynamical effects that mostly contribute to the absorption spectra lineshape and the fluorescence quenching processes upon dye aggregation. Decomposing the various effects arising upon dimer formation, the structure–property relations associated with their optical responses have been deciphered at atomistic resolution.

1 Additional Computational Details

The St-LVC Model and Jahn-Teller Effect
In a tetragonal D 4h system with a pair of degenerate states of the E g or E u symmetry, the E ⊗ b JT Hamiltonian for a single pair of symmetry breaking coordinates is defined as where H 0 is the harmonic oscillator for non JT active coordinates and k 1 and k 2 are gradient and coupling parameters, respectively. Q 1 and Q 2 are modes that carry all gradient or coupling, respectively. Thus, one may find k 1 Q 1 and k 2 Q 2 expressions identical to λ T ii q and λ T ij q in the LVC Hamiltonian for each pair of symmetry-breaking modes. In a D 4h system, these Q 1 and Q 2 modes respectively correspond to the symmetry-breaking b 1g and b 2g coordinates that operate JT distortions for E g electronic states and to b 1u and b 2u coordinates for E u . In the particular case of Zn Pc, the symmetry of the bright pair of electronic states is E u . In the lowest order approach, the JT Hamiltonian is linear on the symmetry breaking coordinates and this suggest that we can safely use our LVC Hamiltonian to properly introduce the JT effect on the spectra, as it shows linear dependency on q as well.
To obtain the linear coupling k 1 and k 2 , identified as λ ij in the notation of the main text, we used the procedure detailed in Ref. S1. First, we displace the geometry of truncated model of Zn PcF12-SH monomer along each dimensionless normal coordinate q α by a small amount ∆ α = ±0.02. Then, for each displaced geometry we compute the adiabatic states |a(∆ α )⟩ and the matrix elements of their overlap with the reference states at equilibrium geometry, i.e. S ij (∆ α )=⟨R i (0) |a j (∆ α )⟩. At each step, this provides us with a new transformation matrix D(±∆ α ), which can be applied to the diagonal matrix of adiabatic energies of the Zn Pc monomer to obtain the diabatic matrix of the potential terms at S-3 displaced geometry Once we get the off-diagonal terms V d ij (q α ), we can apply them to perform a numerical differentiation to compute the linear coupling terms λ ij

Adiabatic-to-Diabatic Transformation for Dimers
In the followed approach, we define diabatic states of truncated Zn PcF12-SH dimer on the basis of reference states |R F rags ⟩ of either the adiabatic states of the isolated fragments/monomers of the dimer (for LEs) or one electron transitions between orbitals on different fragments (for CT states). S1,S2 The diabatic states are then obtained through the following adiabatic-to-diabatic transformation |d⟩ = |a Dimer ⟩ D = |a Dimer ⟩ S T (SS T ) −1/2 (4) where S = ⟨R F rag |a Dimer ⟩ is the overlap of the reference states of the fragments with the adiabatic states of the Zn Pc dimer. Performing the latter operation at the reference geometry yields the transformation matrix D(0). This can be applied to the diagonal matrix containing the energies of the adiabatic states of the Zn Pc dimer to obtain the diabatic matrix, which contains the diabatic energies E d ii (0) on diagonal and inter-state couplings E d ij (0) on off-diagonal.
The correlation functions are obtained from the following expression: where |d i (0)⟩ is the initial wavefunction of the state i and |d j (t)⟩ is the time-evolved wavefunction of state j. Notice that the auto terms are obtained when j = i. Once all the correlation functions are collected, we calculate the "total" correlation function of the systems as the sum of all the correlation functions, each of them weighted by the scalar product of the transition dipoles the corresponding states, i.e.: Finally, the total spectrum is obtained by Fourier transform of the total correlation function ϕ tot . S-5

Diabatic States
The diabatic states for the truncated Zn PcF12-SH dimers considered in this study were defined according to the scheme shown in Figure S1. Figure S1: Sketch of the electronic excitations on the truncated Zn PcF12-SH dimer. The orbital from which the electron is coming from (and the hole is forming) upon excitation is shown in blue, while the orbital to where the electron arrives is highlighted in red. The states where the final orbital is LUMO+1 are labeled with an apostrophe. The diabatic states |L 1 ⟩, |L ′ 1 ⟩, |L 2 ⟩ and |L ′ 2 ⟩ used on the LVC Hamiltonian are linear combinations of the local excitations |LE 1 ⟩, |LE ′ 1 ⟩, |LE 2 ⟩ and |LE ′ 2 ⟩ as defined by the TD-DFT vectors. S-6

ML-MCTDH Settings
For the monomer of truncated Zn PcF12-SH with 213 fast coordinates, a straightforward application of MCTDH method, including only 2 blocks with 6 coordinates (Q, 3 per block) provided us with the converged spectrum for a Gaussian function of the half width at halfmaximum (HWHM) of 0.04 eV (see next section for convergence tests). The first block includes three coordinates, of which two of them respectively carry all the gradients for states |L⟩ and |L ′ ⟩, and the third contains the linear couplings due to JT effects. For dimer instead, with a total number of 426 normal modes and inclusion of the JT-associated λ ij couplings coming from St-LVC diabatization on the monomer, we adopted its ML extension, i.e. ML-MCTDH, which allows for inclusion of more coordinates and ensures getting fullyconverged spectra in a fast and effective way. Accordingly, we adopted 3 blocks of 10Q with a effective collection of 30Q, provided with a hierarchical representation of the LVC Hamiltonian based on HEMS, to compute the spectra with a HWHM of 0.04 eV for the dimer structures taken from MD clustering and quantum mechanical optimization. The first block contains 10Q, each of them containing all gradients the gradient for the 8 diabatic states , CT (2 → 1) and CT ′ (2 → 1), respectively and two additional coordinates for the inter-state JT coupling of states |L 1 ⟩ with |L ′ 1 ⟩ and |L 2 ⟩ with |L ′ 2 ⟩. Mode combination was also utilized for modes with a similar character.
Under this setup, the non-adiabatic QD propagations of the vibronic wave packets on the coupled PESs were performed, adopting the MCTDH and ML-MCTDH for the truncated monomer and dimer models of Zn PcF12-SH, respectively, to compute time-dependent correlation functions in a short dynamics up to 100 fs. This time scale allows to capture and introduce the most important non-adiabatic features on the spectra in a fully-converged low-resolution way. The graphical representation of the MCTDH and ML-MCTDH trees are shown in Figures S2 and S3, respectively.
S-7 Figure S2: Graphical representation of the MCTDH tree adopted in the QD computations of the truncated Zn PcF12-SH monomer (regioisomer C) for a set of 12 normal coordinates. The numbers close to the lines indicate the single particle functions used to combine or the number of primitive functions to describe each normal coordinate. Figure S3: Graphical representation of the ML-MCTDH tree adopted in the QD computations of the truncated Zn PcF12-SH dimers for a set of 30 normal coordinates. The numbers close to the lines indicate the single particle functions used to combine each mode or the number of primitive functions to describe each normal coordinate. Figure S4 shows that the LVC spectrum of the monomer can be considered as converged with 2 blocks (6Q), since increasing to 3 (9Q) or even 4 blocks (12Q) results on very minor changes on the spectral shape.

Convergence Tests
For the dimer, initially we have 10 collective coordinates per block (8 gradients and the linear coupling between the two local excitations of each monomer). In Figure S5 we compared the LVC spectra obtained for OM dimer model with 2 blocks (20Q), 3 blocks (30Q) and 4 blocks (40Q). The spectra are similar, indicating only minor changes in the low-energy band. Therefore, we considered that adopting 3 hierarchy blocks with a effective collection of 30Q is sufficient to reproduce the optical properties of the Zn Pc dimers with this resolution. Figure S6 compares the spectra obtained for the OM dimer model of Zn PcF12-SH where only the four LE states were included, using the gradients from the monomer bright state computed at the TD-DFT level or those obtained from the LVC by displacing along each normal coordinate. We employed 3 blocks of 6Q (4 gradients and the linear coupling between the two local excitations of each monomer) with a effective collection of 18Q. Both provide similar spectral shapes but with some differences on the total width and on the intensity of the maximum. To be consistent with our calculations on the monomer, we used the LVC gradients for the rest of dimer calculations. Notice that the matrix is symmetric and only the upper half is reported. The term on the ij, j ≥ i position represents α λ ij (α), therefore, the terms in the diagonal are the excited state gradients and the off-diagonal term is the excited state coupling.  S-11  2 Synthesis of Zn PcF12-SH

Molecular Dynamics Simulations
We followed the procedure shown in Scheme S1 to synthesis the Zn PcF12-SH compound used in the present study.
Scheme S1: The procedure used in the present study for synthesis of Zn PcF12-SH.
Synthesis of 2(5), 9(12), 16(19)-tri(2,2,3,3-tetrafluoropropoxy)-23-(6-hydroxy hexyl)oxy-phthalocyanine Zn PcF12-OH: 4-[(6-hydroxy hexyl)oxy]phthalonitrile (0.3 g, 1.2 mmol), 4-(2,2,3,3-tetrafluoropropoxy) phthalonitrile (1,50 g, 5,8 mmol), zinc acetate (0,45 g, 2,44 mmol), hexanol (5 mL) and DBU (0.5 mL) were refluxed under argon for 48 h. The green reaction mixture was then poured into hexane (50 mL) and the resulting precipitate filtered. The resulting solid was dissolved in a minimum amount of dichloromethane and precipitated several times from hot n-hexane. The crude product was loaded on a silica gel column, and first eluted by dichloromethane, then a 20/1 dichloromethane / methanol solution to elute the symmetric phthalocyanine, followed by a 10/1 solutionto obtain the The Zn Pc studied in this work (i.e. Zn PcF12-SH), in its truncated model with lateral OCH 3 groups, features 36 regioisomers, resulting from single substitution of one of the nonperipheral or peripheral positions available on each isoindole units (see Figure S7). Therefore, we initially considered all the regioisomers and computed their relative stability and also their excited state properties for the lowest electronic transitions at FC point, as reported in Tables S3-S5. The results show that different regioisomers slightly affect the relative stability of the monomeric forms (see Table S3) and also the energy splitting between the excited states (at most 0.1 eV, see Table S4). For most of the regioisomers, the low-lying excited states S 1 and S 2 feature quasi-degenerate bright states dominated by H→L and H→L+1 transitions, respectively, with more than 90% contributions. The vertical energy gap amounts to 0.04 eV for average values of 1.90 (S 1 ) and 1.94 eV (S 2 ), with standard deviation of ±0.02 eV for both states, as reported in Table S5. The next transition for all regioisomers appears ∼1.5 eV higher than the first two states, and out of Q-band region. Taking into account the full structure, however, the steric restriction imposed by the introduction of the bulky substituents in the non-peripheral positions excludes the regioisomers carrying OCH 3 group at the α positions. In other words, with the precursors we used in the experiment, only four regioisomers A, B, C and D are possible whose vibronic absorption spectra are calculated and discussed in the next section.
S-16 Figure S7: The general scheme used in this study for nomenclature of Zn Pc regioisomers. " * " points to the starting point of nomenclature.

Adiabatic Vibronic Spectra
The experimental absorption spectrum of Zn PcF12-SH in diluted solutions of DMSO, is displayed in Figure S8 in the full frequency range. It exhibits a relatively broad Soret band in 2.76-4.13 eV (300-450 nm) with the maximum at 3.52 eV (352 nm) and another more intense band in 1.65-2.25 eV (550-750 nm) peaking at 1.83 eV (679 nm) which has been assigned to Q-band. The latter intense band is known to be the most important characteristic of the monomeric Zn Pc species. The Q-band is associated with a relatively weaker peak in higher energy side (at ∼2 eV), which has been attributed to vibrational progression of Q-band and also, small contribution from aggregated Zn PcF12-SH species which might be only partially formed in DMSO. Both Soret and Q bonds are indicative of ππ * transitions, as common features of this category of the compounds. Figure S8: The experimental absorption and emission spectra of the diluted solutions of Zn PcF12-SH in DMSO and water in the full frequency range. The intensity of all spectra is normalized to match the maximum intensity of the Q-band in DMSO.
Here we focus on Q-band, since absorption in this region is most affected by aggregation and importantly, fluorescense is registered with excitation falling in this range (1.65-1.80 eV). The spectral shape of the Q band cannot be properly described at pure electronic level, and therefore in the following we run vibronic computations. Figure S9 represents the vibronic absorption spectra computed for the possible regioisomers of truncated Zn PcF12-

S-19
SH monomer in adiabatic approximation with the FC|VG model, and compares them with the experimental absorption spectrum in Q-band region. Although the computed spectra reveal some influences of different regioisomers on the spectral lineshape, the main features of the experimental absorption spectrum are nicely reproduced by all of them, in particular the vibronic structure at ∼2 eV, which is accompanied with a weak progression at ∼2.15 eV. Compared to other regioisomers, the spectrum predicted for regioisomer B exhibits the most appreciable differences with the experiment, indicating the splitting of the lowest peak and an increase of the relative intensity of two peaks. These differences are originated from a notable energy gap of ∼0.1 eV between the two lowest electronic transitions of this regioisomer, as evidenced in Table S4. For the regioisomers A and D, however, the differences are mostly reflected in some slight changes in the relative intensity of the lowest energy peak with respect to the vibronic peak.
Remarkably, the noticeable differences between the regioisomers will wipe out to a large extent by averaging over the four vibronic spectra computed for the regioisomers A-D (see the spectrum shown in orange line in Figure S9), implying various regioisomers can coexist and contribute to the final spectral line shape. Among the possible regioisomers of truncated Zn PcF12-SH, the isomer C delivers the spectrum which resembles more the average and the experiment, concerning the spectral lineshape, spacing and relative intensity of the bands. Notably, substituting the OCH 3 groups with the long lateral chains in regioisomer C does not introduce any significant change on the spectral line shape (except for an extra rigid blue shift of just ∼0.02 eV), as evidenced in Figure S10, confirming the reliability of the truncated model of this regioisomer to simulate the spectroscopy of Zn PcF12-SH.
Furthermore, the regioisomer C features remarkably different optical properties compared to other regioisomers studied here. In fact as evidenced in Table S4,  S-21 Figure S10: Comparison of the adiabatic vibronic FC|VG spectra computed for regioisomer C with that of full structure. The experimental spectrum (in DMSO) in Q-band region is also represented. To facilitate the comparison, all the spectra were normalized to have the same maximum intensity. All stick transitions were convoluted with a Gaussian function with HWHM = 0.04 eV.

S-22
metry, but the π-system still electronically behaves like a D 4h system, exhibiting two nearly degenerate states. The states become less mixed and the virtual orbitals slightly change their character with respect to the symmetric case, but looking at Figure S11, it is evident that like in the symmetric Zn Pc chromophore, both L and L+1 orbitals are still identical with a 90 degree rotation. The system is thus potentially able to undergo JT-like non-adiabatic effects. Accordingly, we further exploit the regioisomer C with the adapted position of the substituents to address this effect on the spectroscopy and photoinduced dynamics of both monomeric and aggregate forms of Zn PcF12-SH, as described in the main text. It is noteworthy that the agreement achieved here for the adiabatic FC|VG spectra with the experiment supports the reliability of harmonic approximation adopted here, implying that the frequency changes and Duschinsky mixing neglected in FC|VG approach play a minor role. This fact is confirmed in Figure S12 by comparing the FC|VG spectrum with the FC|AH one which accounts for these effects. These results also anticipate a negligible role for non-adiabatic interstate couplings in determination of the spectral shape of Zn Pc monomer investigated here, despite of representing two close-lying S 1 and S 2 excited electronic states.
Compared to the vibronic spectrum computed for regioisomer C, the vibronic band in the experimental absorption spectrum features a broad flat region (extended from ∼1.90 to ∼2.02 eV), which can be probably traced back to a partial contribution from self-assembled Figure S12: Comparison of the adiabatic vibronic spectra computed for regioisomer C, adopting different FC|VG and FC|AH approaches. To facilitate the comparison, both the spectra were normalized to match the maximum intensity. All stick transitions were convoluted with a Gaussian function with HWHM = 0.04 eV.
Zn PcF12-SH species. This is clearly evident in Figure S13, if one compares the absorption spectrum with the fluorescence excitation spectra recorded in different wavelengths in DMSO.
In fact, our experiment have revealed that the Zn PcF12-SH in water does not fluorescence, where predominantly exists in self-assembled forms. This suggests formation of aggregated Zn PcF12-SH species in mixture with monomeric forms in DMSO. Since these aggregated species are non-emissive, they do not contribute to the excitation spectra, where the flat region has been apparently removed and the vibronic band appeared sharper than that of the absorption spectrum.

S-24
Figure S13: Comparison of the experimental absorption spectrum in Q-band region with the excitation spectra recorded at different wavelengths for monomer of Zn PcF12-SH in DMSO. All the spectra were normalized to have the same maximum intensity. S-26 Figure S14: The frontier molecular orbitals involved in electronic transitions for the OM dimer, calculated at TD-CAM-B3LYP/6-31G(d,p) level adopting LR-PCM model of solvation in water. The molecular orbitals are represented in both top and side views. Table S7: LVC Hamiltonian in the diabatic representation for the OM dimer. Parametrization at TD-CAM-B3LYP/6-31G(d,p)/LR-PCM(water) level of theory. The matrix is symmetric and only half is reported. Terms in the diagonal represent the vertical energies, E d ii (0), while the off diagonal terms are the constant coupling, E d ij (0). All terms are in eV.
State  Table S8: LVC Hamiltonian in the diabatic representation for the DM dimer. Parametrization at TD-CAM-B3LYP/6-31G(d,p)/LR-PCM(water) level of theory. The matrix is symmetric and only half is reported. Terms in the diagonal represent the vertical energies, E d ii (0), while the off diagonal terms are the constant coupling, E d ij (0). All terms are in eV. State

Species
To investigate the dependence of optical properties of Zn Pc aggregates on the number of molecular units involved in the aggregation, we computed the TD-DFT absorption spectrum of a Zn Pc trimeric species and compared it with that of OM dimer, as shown in Figure S15.
In order to generate a truncated trimer model for the Zn PcF12-SH, we started from OM dimer and added another Zn Pc monomeric unit to reach a H-type trimer structure. This structure was then fully optimized at the CAM-B3LYP-D3/6-31G(d,p) level adopting LR-PCM/water implicit solvation model. To obtain the absorption spectrum, excited states TD-DFT computations were performed employing the same DFT functional and basis set.
The result (while missing vibronic and electronic coupling contributions) shows that trimeric species contributes to the absorption spectrum by altering the relative intensities of the absorption spectrum bands and by adding new peaks with small intensity in the low-energy portion of the Q-band. Therefore, it is expected that accounting for higher order aggregates would increase the relative intensity height of the red tail of the Q-band, while featuring an overall broader and lower band with respect to dimers.
S-30 Figure S15: Comparison of the pure electronic TD-DFT spectra computed for a Zn Pc trimer with that of OM dimer. All the stick transitions were convoluted with a Gaussian function with HWHM = 0.04 eV. No shift has been applied to the computed spectra.
S-31 Figure S16: Comparison of the pure electronic TD-DFT spectra computed for all the Zn Pc dimers considered in this study in Q-band region, with the experimental spectrum recorded for Zn PcF12-SH in water. All the stick transitions were convoluted with a Gaussian function with HWHM = 0.04 eV. For simplicity of comparison, the experimental spectrum was normalized in each panel to match the maximum intensity of theoretical spectrum. No shift has been applied to the spectra. S-32

Definition of Intermolecular Coordinates For Zn Pc Dimers
The intermolecular dimer coordinates employed in this work are defined based on the four important parameters. As indicated in Figure S17a, we defined the d Zn−Zn distance as the distance between the zinc atoms located on the interacting Zn Pc units. The d P −P distance is also defined as an average interatomic distance along a specific axis normal to each Zn Pc ring. In other words, the d P −P is defined as the distance between two imaginary planes, each passing through one of the Zn Pc macrocycle rings with a average distance respect to all the atoms existing in the same ring (the atoms belonging to lateral chains were excluded).
The imaginary planes were defined because the two aromatic macrocycle rings of Zn Pc dimers are not necessary planar and, as indicated in Figure 5c in the main text for the CSs obtained from MD clustering, exhibit significant distortions from planarity. Furthermore, we defined the angles α as the angles between the two zinc atoms with each of the nitrogen atoms located on one of the Zn Pc units, as shown in Figure S17a. The dihedral angle β is also defined according to the Figure S17b Figure S18: Distributions of intermolecular distance between the Zn Pc cores, d Zn−Zn as defined in Figure S17, computed along the MD trajectories in water.

S-34
Figure S19: Distributions of intermolecular angles α for each of the nitrogen atoms located on one of the Zn Pc units, according to the numbering shown in Figure S17, computed along the MD trajectories in water.
S-35 Figure S20: Distributions of the twisting dihedral angle β, as defined in Figure S17, computed along the MD trajectories in water. Figure S21: Distributions of the angle α, averaged over the four angles α 1 , α 2 , α 3 and α 4 (as defined in Figure S17), along the MD trajectories in water. For the sake of clarity, only the values ≤90 • are considered in the distribution.
S-36   S-44 Figure S29: The nonadiabatic vibronic LVC spectra computed for the Zn Pc dimers considered in this study, by including (colored in green) or excluding (colored in red) different contributions from each type of interstate couplings. All vibronic transitions were convoluted with a Gaussian function with HWHM = 0.04 eV. No shift has been applied to the computed spectra.
S-45 Figure S30: Population dynamics of the diabatic states of all the Zn Pc dimers considered in this study, starting photoexcitation from |L ′ 1 ⟩. ML-MCTDH including 30 effective coordinates.

S-46
Table S11: Vertical excitation energies of the minima of the diabatic states according to the LVC models parameterized at different dimer structures with TD-CAM-B3LYP/6-31G(d,p)/LR-PCM(water,nonequilibrium). The zero energy value is different for each dimer and corresponds its ground state. All values are in eV. S-47 Figure S31: Population dynamics of the delocalized states (S 1 -S 8 , top) and diabatic local states (LE + CT, bottom) of CS1 Zn Pc dimer after photoexcitation to the lowest energy dark exciton state S 1 (left) or to the second bright adiabatic state S 4 (right). ML-MCTDH propagations including 30 effective coordinates. The results after a photoexcitation to the bright state S 3 are reported in Figure 10 in the main text.
S-48 S-50 Figure S34: Linear dependence of the total population of the CT states after 100 fs, averaged over initial photoexcitation from the |L 1 ⟩ and |L ′ 1 ⟩ states, on the d Zn−Zn (left panel) and d P −P (right panel) stacking distances for the CS dimers obtained from the MD cluster analysis. S-51