A Nonorthogonal Configuration Interaction Approach to Singlet Fission in Perylenediimide Compounds

Perylenediimide molecules constitute a family of chromophores that undergo singlet fission, a process in which an excited singlet state converts into lower energy triplets on two neighboring molecules, potentially increasing the efficiency of organic solar cells. Here, the nonorthogonal configuration interaction method is applied to study the effect of the different crystal packing of various perylenediimide derivatives on the relative energies of the singlet and triplet states, the intermolecular electronic couplings, and the relative rates for singlet fission. The analysis of the wave functions and electronic couplings reveals that charge transfer states play an important role in the singlet fission mechanism. Dimer conformations where the PDI molecules are at large displacements along the long axis and short on the short axis are posed as the most favorable for singlet fission. The role of the substituent at the imide group has been inspected concluding that, although it has no effect in the energies, for some conformations it significantly influences the electronic couplings, and therefore, replacing this substituent with hydrogen may introduce artifacts in the computational modeling of the PDI molecules.


INTRODUCTION
Singlet fission (SF), a multiple exciton generation process occurring in organic molecules, turns out to be a promising process to increase the conversion efficiency of photovoltaic devices. 1,2SF takes place in organic chromophores in which the excitation energy of the lowest singlet excited state is approximately twice the excitation energy of the lowest triplet, E(S 1 ) ≳ 2E(T 1 ) . 3,4Hence, SF is thermodynamically favorable when the energy difference E(S 1 ) − 2E(T 1 ) is positive and small or slightly negative (by less than ∼0.2 eV), i.e., for slightly exothermic or isoenergetic processes.Large exothermicity induces heat energy losses, which results in a decrease of SF efficiency.Moreover, to avoid possible triplet−triplet deactivation channels, the excited triplet and quintet states should lie at least two times higher in energy than the low-lying triplet state, E(T 2 /Q 1 ) > 2E(T 1 ).Optimal chromophore molecules require good light-harvesting properties, with intense absorption bands in the visible spectra, and be chemical stable. 5For these reasons, π-conjugated organic molecules are potentially suitable SF candidates.
The singlet fission process starts by absorption of light generating a singlet excited state in a molecule, S 1 , that combines with a neighboring molecule in the ground state, S 0 , producing two triplet states T 1 , one on each of the two molecules, overall coupled to singlet spin, 1 TT.This process is radiationless and spin-allowed and hence, is expected to be fast, usually in the picosecond or subpicosecond time scale, thus competing with vibrational relaxations. 6,7In the next step, the resulting singlet state 1 TT can be separated into two decoupled triplet states, thus doubling the number of electron−hole pairs in the system.Despite the fact that SF usually takes place between neighboring molecules in a molecular crystal, it can also occur on a single molecule, so-called intramolecular singlet fission. 8,9lthough SF is a complex, multistep process, and is not yet fully understood, the basic mechanism has been rationalized in terms of a direct or indirect mechanism. 3The direct mechanism involves a straight two-electron process from a S 0 S 1 state, accessed by light absorption, to a 1 TT state, while in the indirect mechanism, the process is mediated by charge transfer states between neighboring chromophores.The involvement of charge transfer states can be envisaged either as virtual states facilitating the 1 TT formation via a superexchange mechanism or as real intermediate states in a sequential two one-electron transfer step process. 7ince SF implies energy transfer between neighboring chromophores, the process depends not only on the nature of the particular chromophore but also on the relative orientation of the monomers in the molecular crystal. 10,11he energetic condition is that the SF process should be slightly exoergic or isoergic.However, this is not the only parameter that controls the SF rate and yield; intermolecular electronic couplings between neighboring molecules also play a key role in the process.In fact, electronic couplings, which allow for nonadiabatic transitions between different states, are markedly dependent on the particular arrangement of the chromophores. 12inglet fission has been observed in different organic materials showing a slightly exothermic or isoenergetic process, such as pentacene, 13−16 1,3-diphenylisobenzofuran (DPIBF), 17,18 carotenoids, 19 and their derivatives, but also in materials as tetracene, 20,21 rubrene, 22,23 perylene, 24,25 and related derivatives, where the process is slightly endothermic.
Perylene-3,4:9,10-bis(dicarboximide) (PDI) molecule (Figure 1) and its derivatives are chromophores that form ordered π-stacked structures where the crystal packing depends on the different functional groups at the N atom of the imide groups (R 1 ).They constitute a vast family of versatile dyes for photophysical applications due to their good light-harvesting properties, showing intense absorptions in the visible spectra and thermal and photochemical stability.Moreover, PDI derivatives are suitable chromophores for singlet fission since the energy of the singlet excited state is approximately the sum of the energy of two triplets, fulfilling the thermodynamic condition. 25,26These characteristics make PDI derivatives ideal systems to explore the optical properties by changing the substituents at the imide groups and, consequently, the crystal packing.This opens the possibility of modifying the structure to tune their photophysical properties at will, eventually targeting the control and increase of the SF efficiency.In fact, the effect of crystal packing in the SF rates and yields has been experimentally studied by Le et al. 25,27 by transient absorption and time-resolved photoemission spectroscopy in six different PDI derivatives, showing that triplet excitons can be produced at the picosecond time scale with a maximum yield of 178% when the PDI molecules are arranged cofacially with a slip displacement along the long x-axis of ∼3.0 Å and minimal short y-axis displacement (see Figure 2 for the definition of the axes).However, the singlet fission created triplets are found to decay to the ground state over tens of nanoseconds.In the theoretical studies of the electronic couplings and SF rates 28 based on the restricted active space configuration interaction with double spin-flip method, 29 it has been postulated that displacements ≥2.5 Å along both axis are more favorable to prompt singlet fission.On the other hand, Mirjani et al. 30 suggested a zero displacement in the short axis and a ∼2.8 Å slip in the long direction as the optimal region for SF.Moreover, the process takes place through a charge transfer enhanced pathway where the charge transfer states act as virtual intermediate states.
In the present work, we study a series of PDI aggregates to analyze the dependence of the two key parameters for singlet fission, SF energy, and intermolecular electronic coupling, on the stacking structure of these molecules in the crystal.From their magnitudes, relative rates are estimated by a simple kinetic model. 28This allows one to find dimer arrangements that are optimal for SF, ultimately helping the design of new advantageous materials.To do so, we employ a nonorthogonal configuration interaction approach, 31,32 based on monomer multiconfigurational wave functions, that fully accounts for orbital relaxation and nondynamical electron correlation.The remaining, mostly dynamical, electron correlation, can be included by a second-order perturbation theory energy correction. 33This methodology has the advantage of yielding very short wave functions expressed in terms of many-electron basis functions with a clear diabatic character that can be directly assigned to specific electronic states such as local excitations or charge transfer states.This vastly facilitates the physical interpretation and analysis of the process under study.Moreover, it allows quantification of the effect of the dynamical electron correlation in the properties of interest.Compared to other theoretical approaches previously applied in the computation of electronic couplings, such as simple models considering only the highest occupied and lowest unoccupied orbitals of the interacting molecules 3,34,35 and density functional theory (DFT) based methods, 36−38 the present approach is based on ab initio multiconfigurational wave functions, which allows researchers to treat all excited states on the same footing and is free of any dependence on the particular choice of the functional.
The paper is organized as follows: the computational approach employed to evaluate the energies and electronic couplings of the PDI derivatives studied is presented in Section 2. Results of the relative energies and character of the relevant electronic states for the PDI monomer and the studied dimers are presented in Section 3, followed by the variation of the energies, electronic couplings, and relative rates with the slip displacement between PDI monomers, from which the optimal geometries for SF can be foreseen.Finally, the effect of the

The Journal of Physical Chemistry A
substituent of the imide group on the energies and electronic couplings has been analyzed.

COMPUTATIONAL APPROACH
The geometry of the ground state, S 0 , of the PDI monomer with the R 1 and R 2 substituents modeled by H atoms (see Figure 1), has been optimized by Density Functional Theory (DFT) using different exchange-correlation functionals and the D3 van der Waals dispersion energy correction by Grimme et al. 39 The functionals applied include two GGA functionals, PBE-D3 40 and BP86-D3, 41,42 two hybrid functionals, B3LYP-D3 43,44 and PBE0-D3, 45 and two range-separated hybrid functionals, ωB97X-D3 46 and CAM-B3LYP. 47Gaussian type basis sets of triple-ζ plus polarization quality 48 (def2-TZVP) were applied for all atoms.All DFT calculations were performed with the Orca 4.2.0 code. 49In all cases, the optimized geometry of the PDI molecule for the S 0 ground state is planar and has a D 2h symmetry.Differences in interatomic distances and angles among the different functionals are small and do not significantly affect the molecular geometry.The largest internuclear distances are found for the GGA functionals, which differ at most by 0.02 Å compared to the ones obtained applying hybrid functionals, while the angles vary by less than 0.5 degrees.The vertical energies of the lowest singlet and triplet states, S 1 and T 1 , have been computed by Time Dependent DFT (TD-DFT) for each optimized structure (Table S1 of the Supporting Information).Results show dispersion of the excitation energies varying the functional, with all values being blue-shifted compared to the available experimental data in solution.The E(S 1 ) − 2E(T 1 ) energy difference is negative in all cases, meaning that the thermodynamic criterion for SF is not satisfied since the process is endoergic by more than 0.2 eV.The structures of the PDI monomer with three different substituents at the R 1 position, a methyl (C1), an ethyl (C2) and an ethylphenyl (EP) group, have been optimized applying the BP86-D3 functional and the same basis set as for the model PDI with R 1 = H.
At the BP86-D3 optimized geometry, state specific Complete Active Space Self Consistent Field (CASSCF) and CAS second-order perturbation theory (CASPT2) calculations on the monomer in vacuum were performed with Open-Molcas. 50,51An active space containing 8 π orbitals and 8 electrons has been considered for the PDI monomer.ANO-RCC basis sets 52 with three different contraction schemes were applied.First, Basis 1 consists of a (5s,4p,2d,1f) contraction for C, N and O atoms and (3s,2p,1d) for H. Basis 2 is composed of a (4s,3p,1d) contraction for C, N and O atoms and (3s,1p) for H. Finally, in Basis 3 a smaller (2s) contraction for the H atoms is applied, while for the rest of the atoms a (4s,3p,1d) contraction is maintained.All electrons were included in the perturbational treatment of the dynamic correlation except the deep-core electrons (1s 2 for C, N and O).The standard CASPT2 zeroth-order Hamiltonian with ionization potentialelectron affinity (IPEA) shift set to 0.25 au, CASPT2-0.25,was used as well as the CASPT2-0 without IPEA shift.As can be seen in Table S2 of the Supporting Information, the size of the basis set does not markedly affect the relative energies of the lowest singlets and triplet states.ANO Basis sets 2 and 3 lead virtually to the same values, which differ by less than 0.1 eV when applying a def2-TZVP basis set.For this reason, in this article, only the results obtained with Basis 2 are reported.Graphical representations of the active orbitals of the monomer are included in Figure S1 of the Supporting Information together with a short description of the computation of the monomer wave functions.
To avoid the so-called IPEA dilemma, 53 N-Electron Valence State Second-Order Perturbation Theory (NEVPT2) 54 calculations, as implemented in OpenMolcas, have been performed.In contrast to CASPT2, NEVPT2 applies a Dyall Hamiltonian as a zeroth-order Hamiltonian, which is intruder state free and does not need parameters like the IPEA shift.The two formulations of NEVPT2, strongly contracted (SC-NEVPT2) and partially contracted (PC-NEVPT2), have been applied.However, in the present implementation in OpenMolcas, NEVPT2 calculations can be performed only with a Density Matrix Renormalization Group (DMRG) wave function as a reference, instead of a standard CASSCF wave function.Hence, DMRG calculations were carried out for all states starting from the converged CASSCF(8,8) monomer wave functions.A value of 2000 for the number of renormalized states, m, and a maximum number of 5 sweeps were sufficient to reproduce the CASSCF (8,8) energies.
The nonorthogonal configuration interaction for fragments (NOCI-F) method as implemented in the GronOR code 32,55 has been used to study the dimers.PDI monomer CASSCF- (8,8) wave functions (computed at the BP86-D3 optimized geometry with Basis 2) for the two lowest singlets, S 0 and S 1 , the lowest triplet, T 1 , and the doublet cationic, D + , and anionic, D − , states have been employed to construct the dimer manyelectron basis functions (MEBFs) as antisymmetrized spinadapted products of the monomer wave functions, resulting in six singlet dimer MEBFs (S 0 S 0 , S 0 S 1 , S 1 S 0 , 1 TT, D + D − , and D − D + ).The NOCI energies and wave functions are obtained by solving the nonorthogonal 6 × 6 Hamiltonian matrix.The electronic coupling between the initial (i) and final (f) states involved in the singlet fission process is defined by the following the expression: 34 where the initial state is a linear combination of MEBFs representing the excited singlet state and the final state is a linear combination of MEBFs corresponding to the singletcoupled double triplet (see section 3.2.1).This method has the advantage that each electronic state of the monomer is expressed in its own optimal molecular orbitals, hence including full orbital relaxation from the start, and leads to compact wave functions, facilitating the analysis and interpretation of the results.However, this implies that the molecular orbitals of the Slater determinants that contribute to the wave function expansion are nonorthogonal to each other, making the calculations computationally much more expensive than the standard orthogonal approaches.Fortunately, the program was optimized to operate in massively parallel supercomputers, allowing NOCI calculations to be efficient and feasible for relatively large molecular systems.To reduce the dimension of the two-electron integrals, the MEBFs are expressed in a common molecular orbital basis 33 instead of the atomic orbital basis used in the optimization of the fragment wave functions.The default threshold of 10 −5 has been used for the removal of linearly dependent basis functions in the common molecular orbital basis of the dimer.The same threshold has been used to select the determinant pairs that are evaluated in the calculation of the nonorthogonal matrix The Journal of Physical Chemistry A elements.These values have been shown to be sufficient to get reliable energies and electronic couplings. 33,56he multiconfigurational CASSCF wave functions used to construct the MEBFs account mainly for nondynamic molecular electron correlation.The energy effects of the remaining, mainly dynamic, electron correlation are considered by shifting the diagonal elements of the NOCI matrix by a second-order perturbation theory correlation energy correction on the MEBFs. 56This correction is calculated from the second-order perturbation electron correlation energy for each electronic state of the monomer.Hence, from now on, we will refer to those NOCI calculations obtained by diagonalizing the Hamiltonian matrix where the MEBFs are constructed from CASSCF wave functions by NOCI(CASSCF).When the diagonal Hamiltonian matrix elements are shifted to include dynamic electron correlation, the results will be labeled with NOCI(CASPT2-0), NOCI(CASPT2-0.25),NOCI(SC-NEVPT2), or NOCI(PC-NEVPT2) depending on the specifics of the second order perturbation approach used to estimate the correlation energy correction.

PDI Monomer.
The crystal packing of the different PDI derivatives (Figure 1) depends on the particular functional group at the N atom of the imide groups (R 1 ).However, the absorption spectra of different PDI derivatives in solution have been shown to be almost indistinguishable from each other, 25 indicating that the substituents do not affect the electronic properties of the PDI molecule.Hence, for simplification, the PDI molecules are usually modeled by substituting both R 1 and R 2 by H atoms. 28 Here, the relative energies of the different states that are potentially involved in the SF process (S 1 , T 1 , D + and D − ) have been calculated for four different PDI derivatives.First, the simplified PDI monomer with H atoms bonded to the imide groups and, in the second place, three different substituents in the R 1 position have been considered: a methyl (C1), an ethyl (C2), and an ethylphenyl group (EP).
The relative energies of the different states with respect to the S 0 ground state are listed in Table 1.Results show that inclusion of the electron correlation by second-order perturbation theory has an important effect on the excitation energies, which largely stabilizes the S 1 and T 1 states.Compared to the experimental estimation of the excitation energies, obtained from the solution-phase absorption spectra of six PDI derivatives, 25,26 2.34 eV for S 1 and 1.19 eV for T 1 , CASPT2-0.25 tends to overestimate the values while CASPT2-0 underestimates them.SC-NEVPT2 and PC-NEVPT2 give similar values; in both cases, excitation energies lie between the CASPT2-0.25 and CASPT2-0 results.The vertical excitation energies reported in Table 1 for the S 1 and T 1 states turn out to be smaller than those computed by DFT calculations applying the ωB97X-D functional, with values for the vertical S 0 to S 1 excitation of 2.84 57 and 2.90 eV, 28 depending on the basis sets, and 1.55 eV for the T 1 excitation. 28Restricted active space configuration interaction with double spin-flip calculations lead to even higher excitations, with S 1 ∼ 0.5 eV and T 1 ∼ 0.3 eV blue-shifted compared to the ωB97X-D results. 28Here, except for the CASPT2-0.25 relative energies, the other second-order electron correlation corrections lead to lower excitation energies than experiment, with the PC-NEVPT2 showing better agreement with spectroscopic data, 25,26 with values underestimated by about 0.4 eV for the S 1 excitation and less than 0.2 eV for the T 1 state.The E(S 1 ) -2 E(T 1 ) energy difference is negative, i.e., singlet fission is not thermodynamically favorable, by −0.6 eV at the CASPT2-0.25 level of calculation and by −0.2 eV for SC-NEVPT2 and PC-NEVPT2, while the process is exoergic by 0.2 eV at the CASPT2-0 level.As can be seen in Table 1, the relative energies obtained with the three substituents, methyl, ethyl, and ethylphenyl, do not significantly differ from those obtained with the model monomer with H, meaning that the states involved in the SF process are not affected by these substituents at the R 1 position.The most remarkable difference is found for the ionized doublet D + state, which is stabilized by ∼0.2 eV when The Journal of Physical Chemistry A the H at the R 1 position is replaced by substituents with larger electron donating character, such as the alkyl groups.

PDI Dimer.
To study the singlet fission process in different PDI derivatives, a minimal model composed of two neighboring monomers has been applied, and various dimer conformations have been considered.−60 These include C1 (R 1 = methyl), C3−I and C3−II (R 1 = propyl) showing two different crystal structures, C7 (R 1 = heptyl), C8 (R 1 = octyl), EP (R 1 = ethylphenyl), and MO (R 1 = methoxyphenyl), while R 2 = H in all cases (see Figure 1).For these materials, the model dimers were first constructed by placing the two monomers with R 1 = H arranged as in the respective crystal structures.The interplane distances (dz) between monomers are 3.40 Å for C1 and C3−I, 3.41 Å for C3−II, 3.50 Å for C7, 3.28 Å for C8, 3.48 Å for EP, and 3.46 Å for MO.Furthermore, to analyze the dependence of the electronic properties of the PDI dimers with the dx and dy slip-stack displacements (Figure 2), 23 configurations have been generated in the region where most of the experimental PDI derivatives are found.All those configurations are placed at dz = 3.48 Å, the one corresponding to the EP derivative, which is one of the candidates that has been suggested to exhibit SF.Notice that for most of the PDI derivatives, the separation in the perpendicular z-direction is very similar, around 3.40−3.50Å. 58−60 Figure 3 shows the whole set of conformations scanned, disclosing the geometries extracted from the crystal structures.
For the entire set of dimer geometrical conformations, NOCI calculations have been performed in order to construct the NOCI wave functions of the six lowest electronic states, analyze their character, and determine both the relative energies between the different dimer states and the electronic couplings.As commented above, one of the main advantages of the NOCI approach is the simple interpretation of the wave functions, since they are written in a short expansion of MEBFs with recognizable character.For all of the dimers, the NOCI wave functions have been analyzed.As an example, Table 2 shows the expansion of the NOCI wave functions in terms of the MEBFs for the PDI dimer in the EP orientation, computed at the CASSCF level and by diagonalizing the NOCI matrix including the electron correlation correction by CASPT2-0.25,CASPT2-0, and PC-NEVPT2 (SC-NEVPT2 results are not reported since they are virtually identical with those obtained by PC-NEVPT2).The relative energies of the NOCI states are also collected in Table 2. Similar NOCI wave functions are obtained for the different calculations performed (with and without dynamic electron correlation correction), although the energetic ordering of the states can be slightly different.As shown in Table 2, the NOCI wave functions show important mixing between the different MEBFs, except for the lowest electronic state, Ψ 1 , which corresponds to the ground state S 0 S 0 configuration.For instance, for the NOCI(CASSCF), Ψ 3 is a combination of S 0 S 1 and S 1 S 0 MEBFs with similar weights, so-called excitonic resonance (ER) states, while Ψ 2 is a combination of ER states, with charge transfer resonance (CR) states, mixtures of D + D − and D − D + charge transfer states between the two monomers.Ψ 4 has a large contribution from the 1 TT state, with important weights of CR states.Ψ 5 and Ψ 6 have large CR contributions, with sizable mixing with 1 TT or ER states, respectively.These large mixtures between different states were already pointed out by Farag and Krylov 28 in their study of SF in PDI dimers, based on restricted active space configuration interaction method with double spin-flip (RAS-2SF).This is in contrast to that found in other SF systems like acenes, diphenylisobenzofuran derivatives or diazadiborinine, 10,11,56,61 where the wave functions appear to be more localized in few dimer states.On the other hand, the CASPT2-0 corrected NOCI wave functions slightly differ from the NOCI(CASSCF), NOCI(CASPT2-0.25) and NOCI-(NEVPT2) wave functions.In particular, the lowest excited state, Ψ 2 , is mainly of 1 TT character, with small contributions from CR states.The lack of IPEA correction in NOCI-(CASPT2-0) overstabilizes the 1 TT state, which now lies below the S 0 S 1 ± S 1 S 0 states.The energy gap of the 1 TT and D 0 D 1 ± D 1 D 0 states also increases, which significantly reduces the mixing among them in the final NOCI wave functions.
In some cases, two states with sizable and almost equal contributions from the 1 TT MEBF are found.For instance, in the NOCI(CASSCF) wave functions for the EP dimer (see Table 2), Ψ 4 and Ψ 5 show similar coefficients in the 1 TT MEBF.Hence, in this case, the state that is lowest in energy has been considered to represent the situation where the two monomer triplet states are coupled to singlet.However, sometimes a higher electronic state shows a larger coefficient in the 1 TT MEBF than a lower-lying state, as found in the NOCI(CASPT2-0.25)wave functions, where Ψ 6 shows a somewhat larger coefficient in the 1 TT MEBF than Ψ 4 .In these cases, as the difference in weights is not substantial, the lowest Ψ 4 state has been chosen as the reference for the 1 TT state.Farag and Krylov 28 found a similar behavior for some PDI dimers (EP, MO, C3−I, C7 and C8).However, in their calculations, the higher state had a much larger 1 TT character that the lower one and the two possible 1 TT states were explored in their work.
Among the different MEBFs considered in the NOCI, the CT configurations imply the movement of electrons among the monomers.Therefore, these configurations are, in principle, most susceptible to polarization effects induced by the environment.Bellinger et al. have reported a detailed study of the relative energies of the Frenkel excitons (S 0 S 1 ± S 1 S 0 ) and the CT states as a function of different representations of the environment of a PDI dimer. 62As the most relevant The Journal of Physical Chemistry A conclusion in the scope of the present study, they revealed that the relative energy of the symmetric CT states, i.e., those that are mainly described by the D + D − ± D + D − configurations, is hardly affected by inclusion of solvent effects.They found a substantial lowering of the energy only when the CT state distorts to its optimal geometry, associated with a collapse of the wave function in either the D + D − or D − D + configuration.As shown in Table 2, the contribution of the CT configurations to the electronic states is symmetric in all cases, which means that the polarization effects in the environment are expected to be small.Moreover, it should be realized that the CT states are never populated in the SF process in PDI, they only contribute to the multiconfigurational wave function of the states involved in the SF process.Therefore, here, geometry distortions are not at play.
The singlet fission process starts by exciting one PDI molecule from the ground state to the lowest singlet excited state.In a dimer, this situation is described as an exciton resonance state of S 0 S 1 ± S 1 S 0 character.Calculations of the NOCI oscillator strengths from the initial S 0 S 0 ground state for the various PDI dimer conformations, show that the state that carries intensity by optical absorption corresponds to S 0 S 1 + S 1 S 0 with an oscillator strength of ∼1.6 au The S 0 S 1 − S 1 S 0 combination is the dark state, showing an oscillator strength of  The Journal of Physical Chemistry A the order of 1 × 10 −5 au.The bright state turns out to be the lowest electronic state for dimer conformations with small displacements along the short axes, dy ≤0.41 Å, (EP and MO among them), while for the rest of the conformations explored, the S 0 S 1 − S 1 S 0 is lower in energy (see as an example Table S3 of the Supporting Information where the NOCI wave functions for the C7 dimer conformation are detailed).The energetic difference between these two states for all conformations is usually small, on the order of 0.1 eV.Although the S 0 S 1 + S 1 S 0 state is the bright state, the electronic coupling of this state with the state of mainly 1 TT character is small.Instead, the coupling between the S 0 S 1 − S 1 S 0 and 1 TT states is much larger (as will be shown below), thus favoring the SF process.As the two states, S 0 S 1 + S 1 S 0 and S 0 S 1 − S 1 S 0 , are singlets and lie close in energy, it is expected that both are accessible after irradiation.Hence, the S 0 S 1 − S 1 S 0 state, which shows larger coupling with the 1 TT state, is taken as a reference to compute the SF energies and electronic couplings.
In Table 3 the SF energy, defined as E SF = E( 1 TT) − E(S 0 S 1 − S 1 S 0 ) is reported.A positive value stands for an endoergic process, while a negative value means an exoergic process, with the 1 TT state lower in energy than the initial S 0 S 1 − S 1 S 0 state.Results in Table 3 reveal that NOCI(CASSCF) calculations show an endoergic process for all of the systems.Inclusion of correlation energy correction in the NOCI calculations by the standard CASPT2-0.25method increases these values while NOCI(CASPT2-0) leads to exoergic energies of around 0.1 eV, which means that singlet fission would be favorable.SC-NEVPT2 and PC-NEVPT2 energy corrections in the NOCI calculations have also been considered, both resulting in very similar values.Singlet fission NOCI(NEVPT2) energies are in between the values obtained by NOCI(CASPT2) with and without IPEA shift.These results are in line with those obtained for the monomer, showing that NOCI(CASPT2-0) tends to overstabilize the triplet states in the monomers.Comparing the singlet fission energy of the different dimers with the corresponding energy difference on the monomer, 2E(T 1 ) − E(S 1 ), it turns out that the process is energetically more favorable in the EP and MO dimers than in the monomer, while it is less favorable for the rest of conformations listed in Table 3. From these results, the conformations corresponding to the EP and MO systems fulfill the thermodynamic requirement, with E SF negative or slightly positive, depending on the particular electron correlation correction included, and hence, from the energetic point of view, they are the most likely candidates to undergo SF.
To investigate the dependency of the singlet fission energy on the dx and dy slip-stack displacements, apart from the seven PDI experimental structures, as commented on before, 23 additional dimer conformations have been studied.Figure 4 shows the variation of the singlet fission energy with the displacements along the long x and the short y axes as heat maps to facilitate the interpretation.Small positive or negative values favor the singlet fission process.It can be seen that all levels of the NOCI calculations, CASSCF, CASPT2-0.25,CASPT2-0 and PC-NEVPT2, coincide in pointing to large x displacements (dx ≥ 2.6 Å) and short y values (dy ≤ 0.8 Å) as the energetically most favorable region for SF (blue-violet The Journal of Physical Chemistry A region in Figure 4).From the seven PDI derivatives reported experimentally, the EP and MO dimers display slip distances in dx and dy that meet this criterion.
3.2.1.Electronic Couplings.The electronic coupling between the S 0 S 1 − S 1 S 0 state and the state where the two monomer triplets are coupled to a singlet, 1 TT, has been computed by equation 1 in two different ways.First, it is a direct electronic coupling between the S 0 S 1 − S 1 S 0 and the 1 TT MEBFs.In the second place, the coupling is computed allowing for a charge transfer enhanced mechanism, in which the effect of the charge transfer excitations between the two monomers is included.This is done by first diagonalizing the S 0 S 1 , S 1 S 0 , D + D − , and D − D + 4 × 4 Hamiltonian matrix leading to new MEBFs.Two eigenvectors are dominated by the S 0 S 1 ± S 1 S 0 combination of the original MEBFs, now dressed with D + D − ± D − D + contributions.Similarly, a dressed 1 TT state is obtained from a 3 × 3 Hamiltonian diagonalization including the 1 TT, D + D − and D − D + MEBFs.After reexpressing the NOCI Hamiltonian and overlap matrix in this new set of MEBFs, the total electronic couplings between the S 0 S 1 − S 1 S 0 and 1 TT states can be calculated using equation 1.The total coupling accounts for both, the direct coupling and the charge transfer enhanced contribution to the coupling.Charge transfer effects on the coupling are expected to be important since the final NOCI wave functions showed large mixing of states of different nature; particularly, 1 TT and S 0 S 1 ± S 1 S 0 substantially mix with CR states (see Table 2).The importance of the involvement of charge transfer states in the singlet fission process has been previously pointed out in tetracene 63 and PDI derivatives. 30,64n Table 4 the calculated electronic couplings are reported.The total coupling is computed by NOCI(CASSCF) calculations and including the electron correlation energy correction as commented before.In the latter case, the charge transfer dressed MEBFs are constructed by diagonalizing the shifted Hamiltonian.Instead, results of the direct coupling are only shown for the NOCI(CASSCF) calculations since inclusion of the electron correlation energy correction, either by CASPT2-0.25,CASPT2-0 or NEVPT2, does not significantly affect the couplings when charge transfer states are not incorporated.
As can be extracted from Table 4, direct couplings are small in all systems.A somewhat larger value is found for the C8 dimer, for which the interplanar distance is a bit shorter than that in the other dimers.Inclusion of charge transfer configurations to describe the 1 TT and S 0 S 1 − S 1 S 0 states largely increases the magnitude of the coupling, by at least 1 order of magnitude.NOCI(CASSCF) values are similar to those obtained by NOCI(CASPT2-0.25) and NOCI-(NEVPT2), while NOCI(CASPT2-0) total electronic couplings are significantly smaller as a consequence of the lesser mixing of the 1 TT with the charge transfer MEBFs.Note that the systems for which the thermodynamic criterion was more favorable to experience SF (EP and MO) do not show the  The Journal of Physical Chemistry A largest coupling.The largest couplings are found for dimer conformations corresponding to the C7 and C8 PDI compounds, while the smaller is observed for the C3−II PDI conformation.
In Figure 5, the electronic couplings, both direct and total, between the 1 TT and S 0 S 1 − S 1 S 0 states are plotted as a function of the dx and dy displacements.The region with larger electronic couplings differs from that found for the favorable E SF energy.Although the magnitude of the couplings is very different, the heat map for the NOCI(CASSCF) direct coupling closely resembles the one obtained for the total coupling NOCI(CASPT2-0) calculations, with a small region of maximum values at dx ≥ 3.0 Å and dy ∼ 1.0−1.2Å.This can be explained by the earlier discussed tendency of NOCI(CASPT2-0) to overstabilize the 1 TT state, leading to very little interaction with charge transfer states.This region of high couplings is also present when the charge transfer effects are incorporated in the total electronic couplings, computed by NOCI(CASSCF), NOCI(CASPT2-0.25) and NOCI(PC-NEVPT2).However, in these calculations, a new region of large couplings is found at intermediate dx distances and dy ∼1.0−1.2Å. NOCI(CASSCF) and NOCI(PC-NEVPT2) total coupling heat maps turn out to be very similar, with a region of large couplings that expands from dx ∼ 2.2−2.6 Å.The difference between these two heat maps and the one obtained by a CASPT2-0.25 energy correction is due to the fact that different dimer conformations studied at dx = 2.35 Å do present two states with significant 1 TT character.In the CASSCF and NEVPT2 NOCI calculations the lowest state shows a sizable 1 TT contribution and has been taken as reference although a higher state with slightly larger 1 TT coefficient exists.Instead, in the NOCI wave functions obtained including the CASPT2-0.25 correction, the coefficient of the lowest state with 1 TT character is too small (≤0.3) to be considered a reasonable description of the 1 TT state, and therefore, a higher state with larger 1 TT character has been considered to evaluate the electronic couplings.Hence, the higher relative energy of 1 TT leads to small couplings in this region.Overall, the results show that charge transfer states play an important role in the mechanism of SF in PDI dimers since the mixing of the 1 TT with CR states significantly strengthens the electronic coupling.Moreover, heat maps reveal that, regardless of the method, the couplings are more favorable at relatively large dx and intermediate dy displacements.
Up to this point, the outcome of the two main properties to describe the SF process slightly differ.The thermodynamic criterion for singlet fission is accomplished for large displacements in the x direction and short y shifts, while the electronic coupling is bigger for rather large dx and intermediate dy displacements.However, none of these two properties by itself can be taken as unique reliable descriptor for optimal SF; 25,28 instead, a combination of both relative energies and electronic coupling would be more suitable to establish the more

The Journal of Physical Chemistry A
favorable conformations for SF to occur.In this sense, kinetic models allow combining both descriptors and become a better tool to assess the feasibility of the process.

Relative Singlet Fission
Rates.Here, we follow the kinetic model proposed by Krylov and co-workers, 10,28,65 according to which singlet fission occurs in two steps: the first step corresponds to the transition from the initial photoexcited state, here S 0 S 1 − S 1 S 0 , to the 1 TT, while the second step is the decoupling and separation of the two triplets.The rate of the first step of the SF process, r 1 , can be estimated by the expression used in reference 28: where NAC reads for nonadiabatic coupling between the S 0 S 1 − S 1 S 0 and 1 TT states.Here, it is computed as the total intermolecular electronic coupling.The parameter α is set to 0.5 as suggested by Kolomeisky et al., 66 assuming that the transition state is placed halfway along the reaction coordinate between reactants and products, β = kT 1 with k the Boltzmann constant, T the temperature, typically T = 300 K, and E SF the SF energy as defined before.Notice that the rate, r 1 , leads to faster processes for large electronic couplings and exothermic or small endothermic SF energies and therefore combines the two main descriptors of the SF process.This approach, however, only allows to compute relative rates in a series of similar compounds, but cannot provide absolute rate values.
Here, we have computed the rates with respect to the lowest value and represent the decimal logarithm of the relative singlet fission rates as a function of the x and y variation; see Figure 6.The relative rates have been calculated from the total electronic couplings and singlet fission energies, computed by NOCI(CASSCF), NOCI(CASPT2-0.25),NOCI(CASPT2-0), and NOCI(PC-NEVPT2) calculations, as previously described.From Figure 6, it can be seen that all of the methods result in similar heat maps.Larger values mean kinetically more favored SF process.These can be found in a region with dx ∼2.6−2.9Å and dy below ∼0.9 Å.It is also remarkable that the rate becomes significant again at dx ∼ 3.2 Å and dy ≤ 0.5 Å, which is compatible with the conformation of the EP dimer.
The region around dx ∼ 3.05 Å and dy ≤ 0.5 Å shows a slow rate, which can be explained by the low electronic couplings in this region, as shown in Figure 5. Regarding the seven derivatives experimentally characterized, the relative rates for the formation of the 1 TT state obtained here posit the EP conformation as the more favorable for singlet fission, followed by the MO conformation, C1, C31, and C8 with similar values, and C32 as the most unfavorable arrangement.These results are in accordance with the relative rates reported by Farag et al., 28 showing the MO derivative as the kinetically most favorable for SF, closely followed by the EP, while the relative rates of the C1, C8 and C31 derivatives remain very similar.These results also concur with the outcomes from timeresolved emission experiments 25 that point to the MO and EP PDI derivatives as the systems undergoing fastest SF, the former within 231 ps and 173% yield, and the second around 260 ps and 178% SF yield, while the rate is 1 order of magnitude slower for the C8 PDI derivative.All in all, the results suggest that small displacements along the short axis y and large displacements along the long axis x are more favorable for the formation of the 1 TT state and, thus, to singlet fission.These findings are in agreement with previous works by Mirjani et al., 30 who suggested a dy = 0 and dx ∼ 2.8 Å as the optimal geometry, and by Le et al., 25 who indicated a slip along the long axis dx ∼ 3.0 Å and minimal short axis displacement.In contrast, Farag et al. 28 postulate displacements ≥2.5 Å for both x and y axes.
The rate of the second step of the SF process, the separation of the dimer 1 TT state into two decoupled triplets, is related to the 1 TT binding energy, E b , that can be estimated by the quintet state approximation as the energy difference between the dimer quintet state, 5 TT, and the 1 TT state, E b = E( 5 TT) − E( 1 TT). 7Notice that the 5 TT state has a pure diabatic TT character and therefore, gives a reasonable estimate of the energy of two isolated triplets, while the dimer 1 TT state mixes with charge resonance states, as discussed before.As shown in Table 5, the separation of the 1 TT state into two triplets is not energetically favored, and the energy needed varies depending on the conformation, from ∼60 meV for the C3−II conformation to ∼320 meV for the EP structure.As previously found for the relative energies, here the binding energies obtained by NOCI(NEVPT2) lie in between the values obtained by NOCI(CASPT2-0.25) and NOCI(CASPT2-0), the later leading to the lowest values.The mixing with charge resonance states stabilizes the 1 TT state, favoring exciton trapping with respect to triplet separation.Inclusion of the CASPT2-0 correction tends to localize the 1 TT state, which shows much less mixing with charge transfer states, and therefore, the energy needed to separate the triplets becomes smaller.
Comparing the various derivatives, it can be observed that conformations with a large rate for the 1 TT formation, such as EP or MO, tend to display energetically unfavorable triplet separation, and consequently, the charge generating process could be hindered by the second step.In contrast, C3−II shows the lowest rate since the interaction of the 1 TT MEBF with charge resonance states is limited, and as a result, triplet separation is energetically more likely.Clearly, to describe the mechanism of the triplet separation process, the estimation of the 1 TT binding energy is not sufficient, and the different competing relaxation pathways have to be considered, which is out of the scope of the present work.Still, the fact that the The Journal of Physical Chemistry A triplet separation is not energetically favored is in accordance with the fast decay of the 1 TT to the ground state observed experimentally. 25.2.3.Effect of the R 1 Substituent in the Dimers.As shown before, the specific substituent at the imide group, R 1 , has no influence on the relative energies of the electronic states of the PDI monomer.Based on this fact, most computational studies model the R 1 substituent by a hydrogen atom using the relative orientation of the monomers as found in the crystal structure to differentiate the derivatives.Here, we quantitatively explore the effects of the R 1 substituent in the properties analyzed, i.e., the singlet fission energy, the electronic couplings and the relative rates for the formation of the 1 TT state.To start, we have studied the PDI dimer in the relative orientation of the EP PDI derivative considering four possible R 1 groups: an H atom, a methyl, an ethyl, and the real ethylphenyl group.The relative energies of the six NOCI states for the four models of the EP dimer are very similar, with only small differences of ∼0.05 eV (see Table S4 of the Supporting Information).Singlet fission energies, E SF , are also similar in all cases, with values of 0.18, 0.43, −0.16, and 0.04 eV for R 1 = H at the NOCI(CASSCF), NOCI(CASPT2-0.25),NOCI(CASPT2-0), and NOCI(PC-NEVPT2) level of calculation and 0.19, 0.46, −0.15, and 0.10 eV, respectively, for R 1 = ethylphenyl.This behavior is general, as shown in the heat maps of the E SF with respect to the x and y displacements for the dimer represented with R 1 = ethyl reported in Figure S2 of the Supporting Information as an example, showing similar trends as those obtained by R 1 = H (see Figure 4), thus indicating that the substituent does not significantly affect the singlet fission energy.
In contrast, the electronic couplings, as reported in Table 6, are more sensitive to the specific substituent, and significant differences appear with values that are around two times larger when R 1 is substituted by the real ethylphenyl group in the EP structure in place of H. From Table 6 it can be observed that the electronic couplings for the EP dimer are very similar for H and methyl at the R 1 position; however, for the ethyl substituent, which is closer to the ethylphenyl group, an important increment of the couplings is observed.The couplings obtained by the ethyl substituent are much closer to those obtained for the real ethylphenyl system, which is an indication that the ethyl group, different from the H or methyl, better represents the EP derivative.This is an interesting result since it shows that, although the substituent does not affect the singlet fission energy, it has an impact on the electronic couplings.
In light of these results, we computed the electronic couplings of all of the PDI derivatives building the dimers based on a monomer with an ethyl group at the R 1 position.The aim is to discern whether the ethyl substituent, which apparently better mimics alkyl groups as propyl (C3−I and C3−II dimers), heptyl (C7) or octyl (C8) than the bare Hatom, has an effect on the electronic couplings in these PDI derivatives.The optimized structure of the monomer with R 1 = ethyl was used to construct the dimers.Comparing the electronic couplings obtained with R 1 = H (Table 4) and R 1 = ethyl (Table 7) one can see that for the C1, C3−I, C3−II and C7 dimer conformations the variation is small; the values do not change by more than 5−7 meV, and for C8 the changes are also moderate, with an 11 meV decrease being the maximum change.These small changes do not alter the relative values of the couplings.However, for the EP conformation, replacing the H-atom at the R 1 position by an ethyl group remarkably increases the magnitude of the couplings, ∼30 meV, practically doubling the values obtained for R 1 = H, while for the MO conformation the coupling decreases by more than ∼15 meV.This reverse trend exchanges the size of the electronic couplings, which are larger for the MO conformation compared to EP with R 1 = H and smaller when R 1 = ethyl.The origin of this effect must be ascribed to crystal packing effects because the same replacement does not affect the electronic structure of the monomer, as illustrated in Table 1.The EP and MO conformations show the lowest values for the dy displacement, 0.41 Å for EP and 0.39 Å for MO, while dx is 3.20 Å for EP and 2.67 Å for MO (Figure 3).A short dy displacement puts the two monomers nearly on top of each other, which combined with a large dx displacement places the ethyl substituent of one monomer close to the backbone of the second monomer (Figure 7 a).Analyzing the structure, one Table 6.Direct and Total Electronic Couplings (in meV) between 1 TT and S 0 S 1 − S 1 S 0 for the PDI Dimer at the EP Conformation with R 1 : H, Methyl (C1), Ethyl (C2), and Ethylphenyl (EP) The Journal of Physical Chemistry A can observe interactions between the H atoms of the ethyl group and the C atoms that are involved in the active orbitals of the S 1 and T 1 states.These contacts are also present in the real EP structure but are missing when either an H or a methyl group is used to model the R 1 substituent (see Figure 7, parts b and c).These interactions are stronger for the EP conformation, which shows the largest dx displacement and, consequently, the largest variation on the electronic coupling; somewhat smaller for the MO conformation, with similar dy but smaller dx, thus diminishing the overlap region of the ethyl of one monomer with the backbone carbon active orbitals of the second monomer.This hydrogen-backbone interaction is only of relative importance for the C8 conformation, with a similar dx displacement as the EP dimer but much larger dy slip-stacking displacement.For the remaining conformations, either with shorter dx displacement (see Figure 7d) or with larger lateral dy displacement (see Figure 7e), no significant intermonomer contacts are found, and the systems can be adequately modeled by an H atom instead of an alkyl group, as corroborated by the electronic couplings reported on Tables 4  and 7. Hence, for conformations with large dx and short dy displacements, due to the interactions of the substituents of one monomer with the carbon backbone of the neighboring monomer, the simplification of the substituent by an H atom is not the most appropriate one to evaluate electronic couplings.Heat maps of the electronic couplings for the scanned x and y slip-stacking displacements of the dimer computed with R 1 = ethyl are reported in Figure S3 of Supporting Information.Since the effect of substituting the R 1 group only influenced some conformations, the heat maps only slightly changed when replacing the H atom with an ethyl.Two regions of high coupling are clearly defined, one at intermediate dx and dy displacements and a small area at large dx and intermediate dy positions.Likely, 1 TT formation rates, which depend on the value of the electronic coupling, are also affected when the H is substituted at R 1 by an ethyl.The corresponding heat maps are presented in Figure S4 of Supporting Information, displaying also only slight changes.Again, the EP conformation emerges as the most favorable arrangement for singlet fission, followed by the MO conformation, and C32 as the less favorable geometrical disposition.
Overall, the R 1 substituent chosen to model the PDI derivatives, does not significantly affect properties like the energies of the states of the monomer, the energies of the NOCI states in the dimers and the E SF energies, however, for some conformations, it has a quantitative influence in the electronic couplings, which is a much smaller magnitude.This is an important finding since it is common practice to simplify the R 1 substituent by an H atom to model chromophores and their derivatives uniquely on the basis of relative energies of the monomers.

CONCLUSIONS
Singlet fission energies, intermolecular electronic couplings, and relative 1 TT formation rates have been computed by nonorthogonal configuration interaction (NOCI) calculations in a series of PDI dimers.The NOCI wave functions are written in terms of many-electron basis functions associated with physically sound electronic distributions in the dimers.This approach accounts for orbital relaxation and nondynamical electron correlation effects, while the remaining dynamical electron correlation has been introduced by a second-order perturbation theory energy correction.Inclusion of dynamical electron correlation effects is essential to properly describe the relative energies of the states involved in the SF process.

The Journal of Physical Chemistry A
Results obtained for various PDI dimer conformations show that the 1 TT state generated in the first step of the singlet fission process strongly mixes with other states, especially with charge transfer states.This mixing favors the electronic coupling between the initial S 0 S 1 − S 1 S 0 state and the 1 TT state.NOCI(CASPT2-0) corrected wave functions show a tendency to overstabilize the 1 TT state leading to a NOCI wave function with small contributions of charge transfer MEBFs, and consequently, the couplings computed including a CASPT2-0 electron correlation correction are systematically smaller than those obtained by NOCI(CASPT2-0.25) and NOCI(NEVPT2), where the 1 TT state strongly mix with charge transfer states.
The geometrical conformations that best fulfill the thermodynamic criterion for singlet fission, i.e., being an exoergic or isoergic process, do not necessarily lead to larger electronic couplings.In fact, none of these criteria by itself can be taken as a unique reliable descriptor for optimal SF, instead, a combination of both relative energy and electronic coupling is important to establish the more favorable conformations for SF to occur.Computed relative rates for the 1 TT formation points to dimer conformations with small displacements along the short axis and large shifts in the long axis as the most favorable region for the generation of the 1 TT state.EP and MO PDI derivatives turn out to be the most suitable PDI compounds for SF.However, estimates of the 1 TT binding energy based on the quintet state approximation reveal that the separation into two triplets is energetically unfavorable, thus inhibiting carrier diffusion.A more complete description of the triplet separation can be obtained by adding an extra monomer to the model and comparing the energy of the T 1 T 1 S 0 state to the T 1 S 0 T 1 state.Moreover, the NOCI calculation on the trimer provides the electronic coupling between the two states.This is subject of ongoing research.
The substituent at the imide functional group does not appreciably influence the relative energies of the different states; however, it has a significant effect on the values of the electronic coupling for some conformations where interactions between the substituent of one monomer and the carbon backbone of the nearest monomer come into play.This is the case for the EP and MO PDI derivatives and, to a lesser extent, the C8 PDI compound.Hence, the approximation of replacing the substituent with a H atom seems to be more appropriate to compute energies than to evaluate the electronic coupling.

Figure 2 .
Figure 2. PDI dimer displaying the displacement along the long (x) and short (y) axes.

Figure 3 .
Figure3.Slip-stacking displacements dx and dy for the PDI dimers considered.Labeled points in red represent geometries that have been taken from reported crystal structures.58−60

Figure 4 .
Figure 4. Singlet fission energy (eV), E SF , computed by NOCI calculations, as a function of the dx and dy displacements for PDI dimers with R 1 = H.

Figure 5 .
Figure 5. Electronic coupling (meV) between the 1 TT and S 0 S 1 − S 1 S 0 states, computed by NOCI calculations, as a function of the dx and dy displacements between two PDI molecules with R 1 = H.

Figure 6 .
Figure 6.log of the relative singlet fission rates, computed by NOCI calculations, as a function of the dx and dy displacements for PDI dimers with R 1 = H.
) containing a representation of the active orbitals of the PDI monomer and the heat maps for the PDI dimer with an ethyl substituent and Tables S1−S4, including relative energies and wave function that support the findings of this study (PDF)■ AUTHOR INFORMATION Corresponding Author C. Sousa − Departament de Cieǹcia de Materials i Química Física and Institut de Química Teorica i Computacional, Universitat de Barcelona, 08028 Barcelona, Spain;

Table 2 .
MEBF Coefficients and Relative Energies (eV) of the Six NOCI Wavefunctions for a PDI Dimer at the EP Conformation with R 1 = H a aIn bold, coefficients larger than 0.15.

Table 3 .
Singlet Fission Energies (in eV), E SF , for the PDI Dimers (with R 1 = H) for Which the Crystal Structure Is Available

Table 4 .
1irect and Total Electronic Couplings (in meV) between1TT and S 0 S 1 − S 1 S 0 for the PDI Dimers (with R 1 = H) for Which the Crystal Structure Is Available

Table 5 .
1TT Binding Energies (in eV), E b , for the PDI Dimers for Which the Crystal Structure Is Available

Table 7 .
1irect and Total Electronic Couplings (in meV) between1TT and S 0 S 1 − S 1 S 0 for the PDI Dimers for Which the Crystal Structure Is Available with R 1 = Ethyl