Morphology, dynamic disorder, and charge transport in an indoloindole-based hole-transporting material from a multi-level theoretical approach

An exponential effort in the design of the hole-transporting materials (HTMs) during the last decade is motivated by their key role as p -type semiconductors for (opto)electronics. Although structure-property relationships have been successfully rationalized to decipher optimal site substitutions, aliphatic chains length or efficient aromatic cores for enhanced charge conduction, the impact of molecular shape, material morphology and dynamic disorder is generally overlooked. In this work, we characterize by means of a multi-level theoretical approach the charge transport properties of a novel planar small-molecule HTM based on the indoloindole aromatic core (IDIDF), and compare it with spherical-like spiro-OMeTAD. Hybrid DFT calculations predict moderate band dispersions in IDIDF associated to the main transport direction characterized by  -  stacked molecules, both between the indoloindole cores and the thiophene groups. Strongly coupled dimers show relevant non-covalent interactions, indicating that NCI surfaces are a necessary but not exclusive requirement for large electronic couplings. We evidence remarkable differences in the site energy standard deviation and electronic coupling distributions between the conduction paths of IDIDF and spiro-OMeTAD. Despite the spherical vs planar shape, however, theoretical calculations predict in the static crystal strong direction-dependent charge transport in the two HTMs, with ca. one-order-of-magnitude higher mobility (  ) for IDIDF. The dynamical disorder promoted by finite temperature effects in the crystal leads to a reduction of the hole transport properties in both HTMs, with maximum  values of 2.42 and 4.2 × 10 −2 cm 2 V −1 s −1 for IDIDF and spiro-OMeTAD, respectively, as well as a significant increase in the transport anisotropy in the latter. Finally, the impact of the material amorphousness in the hole mobility is analysed by modelling a fully random distribution of HTM molecules. An average (lower-bound) mobility of 1.1 × 10 −3 and 4.9 × 10 −5 cm 2 V −1 s −1 is predicted for planar IDIDF and spherical spiro-OMeTAD, respectively, in good accord with the experimental data registered in thin-film devices. Our results demonstrate the strong influence of molecular shape, dynamic structural fluctuations, and crystal morphology on the charge transport, and pose indoloindole-based HTMs as promising materials for organic electronics and photovoltaics.


Introduction
Hole transport materials (HTMs) are one of the crucial components of photovoltaic devices with dye-sensitized (DSC) and perovskite-sensitized (PSC) active layers.HTMs are used to transport photogenerated holes in the active layer to the electrodes, and to prevent mechanical activelayer/electrode contacts.Compared to the inorganic analogues, organic HTMs as a p-type semiconductor have the advantages of being biodegradable, low cost, and easily processed. 1To achieve high power conversion efficiencies (PCEs) and long-term device stabilities, however, novel HTM with a proper structural and chemical design, appropriate energy level offsets, high crystallinity and high hole mobility (> 10 −3 cm 2 V −1 s −1 ) are required. 24][5][6] Among them, spiro-OMeTAD [2,2′,7,7′-tetrakis(N,N-di-p-methoxyphenylamine)9,9′-spirobifluorene], consisting of two perpendicular fluorene moieties with methoxyphenyl-amine substituents (see Figure 1), is arguably one of the most recognized organic HTMs, boosting the rapid development of PSCs exceeding PCEs of 21%. 7The spiro-type central structure confers an orthogonal configuration, which acts in favour of morphological stability.Meanwhile, the peripheral arylamine units in spiro-OMeTAD allow suitable energy level alignment, and the pendant methoxy groups (-OMe) facilitate solution processability.High mobilities for single-crystal field-effect transistors have been reported around 10 −3 cm 2 V −1 s −1 , three orders of magnitude higher than that found in thin-film devices; 8 thus, crystallinity plays a key role in device engineering for optimal performance.The substituent position of −OMe group in spiro-OMeTAD is found to have strong impacts on material properties, such as the electron-donor/withdrawing character, allowing tunability of the optical bandgap. 9On the other hand, the adjustment of the electron donating ability through the pendant groups is proved to be another efficient way to fine-tune the energy levels, e.g.substituting the −OMe groups by ethyl, N,N-dimethylamino or methylsulfanyl moieties. 10Still, spiro-OMeTAD HTMs present several disadvantages such as the requirement of dopants and additives to increase charge transport at the expense of enhanced degradation.Moreover, the complex multistep procedure for its synthesis and purification hinders commercialization due to exceptionally high costs.
HTMs based on other spiro-type structures have been developed, including heteroyclic moieties based on nitrogen, 11-13 oxygen [14][15][16] and sulphur. 17On the other hand, non-spiro core-based HTMs containing thiophenes, 18-20 triazatruxenes, 21 azulenes, 22 tetrathiafulvalene, 23 carbazole 24,25 or phenothiazine, 26 among others, have shown promising PCEs ranging 10−23%, including the poly[bis(4-phenyl)(2,4,6-trimethylphenyl)amine (PTAA) and poly(3hexylthiophene-2,5-diyl) (P3HT) systems, which hold record PCEs of ca.22% among polymeric HTMs. 27,28 linear fluorinated indolo [3,2-b]indole core unit-based HTM (IDIDF) was reported by Park et al. (Figure 1), 29 showing one order of magnitude higher hole mobility (1.7 × 10 −3 cm 2 V −1 s −1 ) compared to referable spiro-OMeTAD in space-charge-limited-current measurements, and superior photovoltaic performance for perovskite solar cells.The high charge-carrier transport was explained by the tight and parallel molecular arrangement stacked by strong π-π interactions, which contrasts to the spherical shape of spiro-OMeTAD.The indoloindole core was previously used to characterize the nonfluorinated 4H4TIDID analogue, which presented an extraordinary field-effect mobility of 0.97 cm 2 V −1 s −1 for a vacuum-deposited polycrystalline device.30 More recently, a family of IDIDF HTMs has been synthesized with varying alkyl chain position and length, leading to PCEs as high as 23% in non-doped and 24% in doped conditions, as well as outstanding thermal and moisture stabilities.31 Despite the promising characteristics of this new class of indoloindole HTMs, a detailed analysis on the importance of a well-packed - stacking and the dynamic disorder, as well as the impact of mesoscale ordering and anisotropy, in their charge-carrier transport properties is missing so far.
From the experimental point of view, manufacturing variability, layer homogeneity and/or crystallinity, interface morphology, structural orientation, presence of trap defects or stability issues, among others, inhibit precise comparison of the hole mobility in a series of HTMs.In this regard, theoretical chemistry and computational modelling has proved crucial to deepen into the understanding of charge transport and the key factors affecting it, boosting the development of next-generation hole-transport materials.For example, studies employing state-of-the-art ab-initio calculations, multiscale modelling and/or kinetic Monte Carlo simulations shed light on the impact of anisotropy, 32 mesoscale ordering, 33 amorphousness, 34 substituent connectivity 35,36 or aliphatic side chain engineering in the charge transport of HTMs, allowing guidance prior experimentation.The role of disorder, nature of the semiconductor-dielectric interface, and presence of non-Ohmic contacts can also be assessed theoretically, allowing for the first time a robust molecular-scale description of a field-effect transistor operation. 37Finally, big-data science and machine-learning algorithms are being applied to high-throughput screening of HTMs for the discovery of new candidates with unprecedented performances or develop costeffective theoretical approaches for large-scale modelling. 38,39 this work, we perform a detailed analysis on the charge transport properties of a recently reported indoloindole-based hole-transporting material (IDIDF), and compare it throughout with referable spiro-OMeTAD.First, the crystal and electronic structure is calculated, and the different paths for hole transport are characterized.A multi-level theoretical protocol combining ab-initio calculations, molecular dynamics and kinetic Monte Carlo simulations was applied to predict hole mobilities in the static crystal, and to determine the effect of dynamic disorder by means of finite temperature structural fluctuations.Non-covalent interactions were analysed for most representative conduction paths.The impact of the molecular shape and the selfassembly crystal packing along the efficient charge-transport paths was disentangled by comparing the structural isotropy, site energy and electronic coupling distributions on the planar-like IDIDF vs the spherical-like spiro-OMeTAD.Finally, random amorphous phases of the two HTMs were modelled to unravel the impact of crystallinity in their charge transport properties.

Computational details
Minimum-energy crystal structures were obtained under the density functional theory (DFT) framework by using the Gradient-Generalized Approximation through the PBEsol functional, 40 a 4×4×4 k-point grid mesh, and the "light" tier-1 numerically atom-centred orbital basis set, as implemented in FHI-aims. 41Single-point electronic structure calculations were performed at the hybrid meta-GGA HSE06 functional using the same k-point grid mesh to explore the Brillouin zone.A k-path following the high-symmetry points of was chosen according to the centrosymmetric triclinic space group of spiro-OMeTAD and IDIDF. 1 P Effectives masses for holes and electrons were obtained according to the finite differences method as implemented in the effmass code. 42Frontier crystal and molecular orbitals were represented by means of the VESTA 43 and Chemcraft software. 44e shape of the HTMs was evaluated by employing an isotropic index (ISO) calculated as the fraction between the smallest and the largest eigenvalue of the inertia tensor.To do so, the centre of mass was first computed, and then the coordinates were translated to that point.Subsequently, the inertia tensor was calculated, and the eigenvalues were obtained.Using this parameter, a value of ISO = 1 indicates a fully isotropic system, such as methane or fullerene, and an ISO approaching 0 is related to a highly anisotropic structure.
The non-covalent interactions involved in relevant dimers of the HTM materials were analysed through the reduced density gradient (RDG), or NCI index, by means of the NCIPLOT package. 45romolecular densities, which are accurate enough for a semi-quantitative analysis of the noncovalent interactions, were used to save computational resources.Graphical representations of the NCI index were performed by means of the VMD software using standard cutoffs of  = 0.5 a.u. and RDG = 1.0 a.u.for the calculation, and of RDG = 0.3 a.u.for the isosurface visualization.
According to the high localization of the polaron generated in a single molecule, we applied the Marcus theory for calculating the rate constants between dimers on each relevant direction through the formula: (2) where is the reorganization energy, is the Boltzmann constant, is the temperature (in  B k T this case, 298.15 K), V is the electronic coupling, and is the energy difference between the E hopping sites.
The reorganization energy was obtained by means of the four-point approximation: (3) where is the energy of the state X (0 = neutral or + = cation) at the minimum-energy X Y E geometry of Y (g0 = neutral or g+ = cation).We assume that the external reorganization energy in a crystal is negligible compared to the internal . 46,47The molecular geometries and the corresponding energies at the different states were obtained using the Boese and Martin's τdependent 2004 hybrid functional BMK 48 and the 6-31G(d,p) basis set in gas phase.Note that J.-L.Brédas and coworkers 8 evidenced a large underestimation of the reorganization energy at the B3LYP level, caused by the well-known delocalization error (overestimation of conjugation) inherent to this functional, whereas a negligible impact of the functional is found on the electronic coupling. 49e electronic coupling V between dimers was estimated according to the dimer Fock operator with frontier orbitals of two neighbouring monomers based on a dimer projection (DIPRO). 50In this approximation, the adiabatic molecular orbitals of a dimer are expressed as a linear combination of the molecular orbitals of the isolated monomers, and the monomer diabatic energies and the electronic coupling are obtained after solving the corresponding secular equation.The wavefunction used to compute the electronic coupling was obtained at the B3LYP/6-31G(d,p) level of theory.
The mobility of charge carriers was assessed by considering the most relevant paths.In the crystalline phases, we considered a 5x5x5 supercell for IDIDF (125 molecules) and 3x3x3 for spiro-OMeTAD (54 molecules), and the crystalline structure was minimized using the experimental crystal cell parameters at the molecular mechanics level under the scalable molecular dynamics NAMD program package. 51The web-based graphical interface CHARMM-GUI was employed for the automatic parameterization of the HTMs using the CHARMM generalized force-field (CGenFF). 52A kinetic Monte Carlo program was coded to compute the hole mobility according to the direct method based on the Gillespie algorithm (see Supplementary Information for further details). 53All possible paths were considered for firstnearest and second-nearest neighbours.Electronic couplings between dimers higher than a threshold of 1 meV were chosen as non-negligible charge transport paths.For kMC diffusion convergence, 5000 simulations of 5 and 1000 ns for crystalline and amorphous phases, respectively, were used, whereas an electric field of 10 4 V cm −1 was applied in all the directions of ab, ac and bc crystallographic planes for static and dynamic crystal cells, or xy, xz and yz coordinates planes in the case of the amorphous materials.
Dynamic disorder due to finite temperature structural fluctuations on the charge transport properties of our HTMs were considered by performing NVT molecular dynamics (MD) simulations using the experimental crystal lattice parameters and the CGenFF force-field.MD simulations were performed along 500 ps with a time step of 2 fs and a temperature of 298 K. Crystal structure snapshots were extracted every 1 ps of the last 200 ps in a total of 200 geometries.The previously selected dimers with relevant electronic couplings for the static crystal structure were analysed during the dynamics to obtain a Gaussian distribution of energy sites ( ± ) and couplings ( ± ).These parameters were set as an input to perform kMC   simulations and predict the dynamic disorder effect of temperature in the mobility (see the Electronic Supplementary Information, ESI).Anisotropy in charge transport was elucidated my computing the mobility along the ab, ac and bc crystallographic planes in steps of 10º.Anisotropy plots of mobility were generated by means of the matplotlib library in python. 54nally, amorphous phases were modelled starting from a random molecular packing including 54 and 100 molecules for spiro-OMeTAD and IDIDF materials, respectively, using the PackMol software. 55,56Lattice parameters were fully optimized in an orthorhombic shape until convergence.In this case, long MD simulations of up to 200 ns using a time step of 2 fs were performed to explore the large conformational space and randomness of molecular positions (snapshot extraction every 1 ns).In contrast to the dynamic disorder analysis, all first-nearest neighbours were considered for dimers, and joined in a unique distribution.For the kMC simulations, a regular grid of equidistant cubic-like disposition of the molecules with a mean distance as extracted from the molecular dynamics of 200 ns were assumed, with the site energies and electronic couplings derived from normal distributions as extracted along the dynamics.

Static crystals
Minimum-energy crystal structures of indoloindole-based IDIDF and referable spiro-OMeTAD HTMs were obtained upon full ionic and unit cell relaxation at the PBEsol level (Figure 2).Theoretical calculations predict lattice parameters very close to the experimental ones, within a maximum deviation of 1.7 and 1.6% for IDIDF and spiro-OMeTAD, respectively (see Table S1  and S2).Inclusion of the van-der-Waals (vdW) dispersion correction to account for the noncovalent interactions leads to an over-compression of the unit cell, resulting in larger discrepancies with respect to the experimental data.Although its better performance is probably due to error cancelation (neglect of thermal effects and dispersion interactions), the dispersion non-corrected methodology was selected for subsequent electronic structure analysis.As shown by the highest-occupied molecular orbital of the molecular systems (Figure S2), the hole in IDIDF is generated over the full -conjugated backbone constituted by the indoloindole core and the peripheral bithiophene moieties, with negligible participation of the aliphatic chains.On the other hand, the hole in spiro-OMeTAD is delocalized over the fluorene core and extends through the methoxyphenyl-amine pendant groups.Crystal orbital representation of the VBM and CBM of IDIDF indicates a certain delocalization of the electronic structure (Figure S3), suggesting a limit case to consider a hopping mechanism.Similarly, although presenting a flat band structure, the VBM in spiro-OMeTAD is partially shared between proximal fluorenes of vicinal moieties (Figure S3).The localized nature of the hole charge carrier in both HTMs is however supported by the single-molecule localization of the spin density upon considering hole formation (IDIDF + and spiro-OMeTAD + ) in the most interacting dimers (Figure S4).

Faraday Discussions Accepted Manuscript
To further confirm the viability of applying a hopping transport regime in our materials, the reorganization energy () of spiro-OMeTAD and IDIDF was computed at the BMK/6-31(d,p) level, and compared it with a rough estimation of the electronic coupling extracted from the VBM depth.Theoretical calculations indicate that spiro-OMeTAD presents a reorganization energy of 302 meV, in good accord with previous reports. 49On the other hand, the highly conjugated structure of IDIDF is predicted with a similar  (313 meV) compared to referable spiro-OMeTAD, supporting its suitability for hole conduction.In both cases,  is notably larger than the electronic coupling estimated preliminarily at 0.13 and 67 meV in spiro-OMeTAD and IDIDF, respectively, from the VBM depth.These results therefore confirm the valid assumption of a hopping regime and the application of Marcus theory to study the charge transport properties in these HTMs. 57e different shape of the spherical-like spiro-OMeTAD and the more planar and longitudinal IDIDF molecules may be translated into a different crystal morphology and hole transport directionality.To measure the level of geometrical isotropy in these two HTMs, we calculated the ratio between the smallest and largest eigenvalue of the inertia tensor, hereafter ISO index (values close to 1 indicate spherical structures and ISO close to 0 imply strong anisotropy; see Computational Details).In particular, we predict an ISO index of 0.688 and 0.091 for spiro-OMeTAD and IDIDF, respectively, which emphasizes the high level of structural anisotropy of the latter compared to the spiro HTM.This level of anisotropy, accompanied by the consequent self-assembling morphology in the crystal, is translated into preferential directions for charge transport.For example, in the crystal structure of spiro-OMeTAD, we predict up to 4 directions where relevant electronic couplings between neighbours must be considered for subsequent mobility calculations (Figure 4).In contrast, IDIDF presents only two directions with important electronic coupling (direction with V = 63 meV, and direction with V = 17 meV).Note 100 211 that one of the main transport directions in IDIDF is predicted between second neighbours, suggesting the necessity of exploring further than first neighbours in highly anisotropic HTMs.Non-covalent interactions (NCI) index was calculated to unveil the relation between the weak interacting forces governing the crystal assembly and the coupling between vicinal molecules.In IDIDF, the main conduction path is characterized by dimers that interact with a large 100 surface of NCIs, as evidenced by the graphical representation of the reduced density gradient in Figure 5 (see Computational Details).The NCI surface confirms the strong interaction between the aromatic cores of vicinal IDIDF entities, with - distances of 3.3−3.6Å between indoloindole cores and of 3.8−4.0Å between thiophene moieties, as well as the proximity between the aliphatic chains (ca.2.8 Å) that further stabilize the packing through CH•••HC aliphatic interactions (Figure 5a).On the other hand, the NCI analysis for the dimer predicts short 211 contacts of 3.8 Å between the thiophene moieties, and the interaction extends through the aliphatic chains with the peripheral benzene rings of the fluoro-indole groups through CH••• contacts at 2.8 Å (Figure 5b).Otherwise, the main interacting dimers in the spiro-OMeTAD crystal are calculated with a large NCI surface (Figure 5c and 5d).Note that the unit cell of spiro-OMeTAD contains two inequivalent molecules so we compute two different types of dimers along a direction ( with next cell and within the same cell; Figure 4).Despite the large 100 000 differences in terms of V, these dimers present similar NCI surfaces located mainly between the aromatic rings of the interacting monomers, with - distances as short as 3.2−4.0Å.The other directions with weak electronic couplings are characterized by relatively small NCI surfaces involving the peripheral methoxyphenyl groups (Figure 5e and S5).

Dynamic disorder
To shed light on the dynamic disorder promoted by structural fluctuations at room temperature on the charge transport, we performed NVT molecular dynamics on the minimum-energy crystal structure of spherical-like spiro-OMeTAD and planar IDIDF up to 500 ps.The timestep was chose to be 2 fs and snapshot extraction was made every 1 ps for the last 200 ps to have vibrational motion resolution.The electronic coupling was sampled for the most relevant dimers obtained in the previous section for the static crystals.First, we recomputed the isotropic ISO parameter along the 200 snapshots of the dynamics simulation.Theoretical calculations predict normal distributions with a mean ISO value of 0.657 ( = 0.012) and 0.089 ( = 0.003) for spiro-OMeTAD and IDIDF, respectively (Figure S6), to be compared with 0.688 and 0.091 obtained in the static crystal.These results indicate that, whereas the IDIDF molecule maintains its large anisotropy, a significant increase in anisotropy is predicted in spherical-like spiro-OMeTAD due to the roomtemperature structural deformations.
The site energy and electronic coupling distribution along the dynamics between the relevant monomer pairs was calculated and fitted to a Gaussian function.In IDIDF, the site energies nicely correlate to a normal distribution, with a mean value centred in zero (Figure 6a in range of 399−481 meV), which are notoriously larger than for IDIDF (Figure 6b).The higher  values obtained in spiro-OMeTAD indicates stronger structural fluctuations along the dynamics, which is in good accord with the dispersion trends of the isotropic index.In terms of electronic couplings, distributions with [ = 16,  = 16 meV] for dimer and [ = 4,  = 6 meV] for V 100 V dimer were obtained, whereas the electronic couplings in the and directions are 000 010 111 practically vanished upon dynamic disorder (Figure 6b).These results not only indicate a lower charge transport but also an increased mobility anisotropy upon considering thermal fluctuations.Hole mobilities were estimated by means of a homemade kinetic Monte Carlo code, which applies a Marcus-type rate for a list of preferential directions using a stochastic method (see Computational Details and ESI).First, the mobilities in the static crystals (0 K) of IDIDF and spiro-OMeTAD were calculated along the main crystallographic planes, according to the preferred conduction paths and the electronic couplings obtained for the most relevant dimers (V larger than 1 meV).Theoretical calculations predict a high level of anisotropy in the case of IDIDF, with a maximum  value of 6.81, 0.51 and 2.94 cm 2 V −1 s −1 in the a, b and c axes, respectively (Figure 7a and Table S3).Note that the maximum mobility does not fall exactly at the a axis due to the non-negligible contribution of the path.Similarly, the spherical spiro-OMeTAD also shows 211 a high level of anisotropic hole conduction (Figure 7b), as evidenced in previous reports. 8,58In this case, maximum  values of 0.30, 0.12 and 0.09 cm 2 V −1 s −1 , are predicted in a, b and c directions, respectively.These results for the static crystals indicate a significantly more efficient charge transport (ca.one order of magnitude) in the planar IDIDF compared to spiro-OMeTAD.
Dynamic effects were considered in the calculation of hole mobilities through the normal distributions of site energies and electronic couplings predicted along the molecular dynamics at room temperature (298 K). kMC simulations indicate a notable reduction of the charge transport by a factor of 3 in IDIDF and of 8 in the case of spiro-OMeTAD due to the finite temperature structural fluctuations.This leads to maximum  values of 2.42 cm 2 V −1 s −1 for IDIDF and of 4.2 × 10 −2 cm 2 V −1 s −1 for spiro-OMeTAD (Table 1).These values are in good correlation with the experimentally reported mobilities of 0.97 cm 2 V −1 s −1 for the non-fluorinated analogue 4H4TIDID and 1.3 × 10 −3 cm 2 V −1 s −1 for spiro-OMeTAD in single-crystal devices. 8,30Whereas the mobility profiles along the three crystallographic planes are maintained in IDIDF upon considering dynamic disorder (Figure 7a), the anisotropy in mobility is increased for spiro-OMeTAD, especially in the ab plane (Figure 7b), as anticipated by the electronic coupling analysis.

Amorphous phase
Perfect crystal morphology of HTMs leads to a theoretical upper bound mobility compared to the experimental data, allowing a fair comparison between theory and experiment only when single-crystal devices can be fabricated.In all other cases, a certain amount of paracrystallinity, semicrystallinity or even amorphousness are present in the HTM materials constituting the (opto)electronic device.In order to compare planar IDIDF vs spherical spiro-OMeTAD in the lower-bound case scenario of having an amorphous material, we modelled large unit cells of up to 100 and 54 molecular structures for IDIDF and spiro-OMeTAD, respectively, where the molecules were placed randomly, and the size of the orthorhombic cell was varied until convergence.Then, a long molecular dynamics simulation of 200 ns in an NPT ensemble was performed.
Representative snapshots of the amorphous HTMs are displayed in Figure 8a for IDIDF and Figure S7 for spiro-OMeTAD.The amorphous phase of the IDIDF HTM is characterized by a density of 1.11 g/cm 3 , slightly smaller than that of the perfect crystal (1.22 g/cm 3 ).Unexpectedly, and despite its spherical-like shape, the amorphous phase of the highly isotropic spiro-OMeTAD is predicted with a converged density of 1.02 g/cm 3 , significantly smaller than that of the perfect crystal (1.25 g/cm 3 ).This result may be explained by the unachieved intercalation between fluorene moieties of vicinal spiro-OMeTAD molecules in the amorphous phase, hindered by the bulky methoxyphenylamine groups.Radial distribution functions (g(r)) of the central molecule centroid with respect to the rest of the molecules evidence for IDIDF a shift to smaller distances (ca. 3 Å) of the shortest contacts compared to the crystal (3.5 Å) due to CH•••(indoloindole) interactions.Whereas in the perfect crystal a g(r) profile with three maxima is predicted at 4, 8 and 12 Å, the amorphous phase is predicted with a smoother g(r) practically stabilized at 5 Å (Figure 8b).At this distance, 0.20 and 0.13 molecules are integrated for IDIDF crystal and amorphous phases, respectively, whereas 2.71 and 1.66 molecules are contained in a sphere of 10 Å.On the other hand, the crystal and amorphous phases of spiro-OMeTAD present a g(r) profile with a steep slope from 3.5 to 5 Å, and then a progressive increase up to ca. 11 Å (Figure S7b).These distances are correlated with the internal (the spirobifluorene core) and external (the whole molecule) radius calculated at 4.3 and 11.2 Å, respectively.In spiro-OMeTAD, 1.76 and 0.78 molecules are integrated at 10 Å for the crystal and the amorphous phase, respectively, which is in line with the lower density predicted for the latter.Site energy and electronic coupling distributions were obtained for the nearest neighbours taking the central molecule as a reference in the amorphous phase of spiro-OMeTAD and IDIDF along the molecular dynamics simulations.In contrast to the dynamic disorder calculations, the long dynamics performed along 200 ns with geometry extractions each 1 ns inhibits vibrational resolution but instead we reach the time scale of conformational changes.Site energies nicely fit to a normal distribution with parameters = −87;  = 1881 meV] and [ = −41,  = 3232 [  meV] for spiro-OMeTAD and IDIDF, respectively.However, electronic couplings are highly centred at zero with a normal distribution that slightly overestimates the probability of having moderate V values (Figure S8): [ = 0.03,  = 1.62 meV] and [ = 0.14,  = 4.63 meV] for spiro-V V OMeTAD and IDIDF, respectively.
The stabilizing intermolecular NCI surfaces were integrated along the nearest-neighbour dimers extracted for the central molecule from the molecular dynamics of the amorphous phase, and compared with the value of the electronic coupling.Theoretical calculations indicate that there is not a clear correlation between the integrated volume of  associated to the attractive NCIs and V (Figure 9).However, we predict a triangular profile in IDIDF (Figure 9a), suggesting that NCI surfaces are important but not enough for predicting the electronic coupling trends (see e.g. the NCI plots for characteristic dimers A, B and C in Figure S9).In contrast, spiro-OMeTAD does not present strong electronic couplings although large NCI surfaces are integrated in several dimers (Figure 9b).Not only the proximity of the molecules but also the nature of this noncovalent interaction and the wavefunction symmetry must therefore be considered for predicting trends in V. Finally, hole mobilities were calculated from the kMC simulations and compared them with the previous data obtained for the static (0 K) and dynamic (298 K) crystals.Note that in the amorphous phase, all the first nearest neighbours were considered in the normal distribution for energy sites and electronic couplings (Figure S8).Theoretical calculations indicate an isotropic hole mobility profile in the amorphous phase of both spiro-OMeTAD and IDIDF, as expected; the average values of  are summarized in Table 1.Compared to the crystal at 298 K, amorphousness leads to a decrease in hole transport of ca. 3 orders of magnitude for the two HTMs, with mean  values of 1.1 × 10 −3 and 4.9 × 10 −5 cm 2 V −1 s −1 for IDIDF and spiro-OMeTAD, respectively (Table 1).These data, which are in good accord with the experimentally reported mobilities of 1.7 × 10 −3 and 4 × 10 −5 cm 2 V −1 s −1 for IDIDF and spiro-OMeTAD in thin-film devices, 29,59 confirm the efficient charge transport of the indoloindole-based HTM for future application in organic electronics and photovoltaics.Table 1.Maximum hole mobilities (cm 2 V −1 s −1 ) calculated for the static crystal (0 K), including dynamic disorder (298 K), and in the amorphous phase for IDIDF and spiro-OMeTAD HTMs.Experimental data available is included for comparison (SC and TF refer to values reported from single-crystal and thin-film devices, respectively).

Conclusions
Herein, we report on the charge transport properties of a novel hole-transporting material based on the indoloindole core with peripheral thiophene groups by using a multi-level theoretical approach.The impact of its anisotropic shape, the dynamic disorder and the morphology of the material are analysed in detail in comparison with referable spiro-OMeTAD.First-principles theoretical calculations evidence in IDIDF a dispersive valence band maximum with small effective masses of the order of 1 m 0 , associated to the main transport direction characterized by the - stacking of the aromatic backbone.Extensive non-covalent interactions are calculated for the most-interacting dimers, but their nature and wavefunction symmetry considerations must be examined for electronic coupling predictive power.Molecular dynamics and kinetic Monte Carlo simulations were performed to unveil the impact of the dynamic disorder in the charge transport of our HTMs, leading to  values of 2.42 cm 2 V −1 s −1 for IDIDF and 4 × 10 −2 cm 2 V −1 s −1 for spiro-OMeTAD (3 and 10 times lower mobilities, respectively, than in the static crystal), in good accord with reported experimental data.Despite the different molecular shape, we predict strong mobility anisotropy in both the spherical-like spiro-OMeTAD and the planar IDIDF, which is further enhanced upon dynamic structural fluctuations.Finally, amorphous phases were modelled with random distribution of the HTM molecules, leading to lower-bound mobilities of 1.1 × 10 −3 and 4.9 × 10 −5 cm 2 V −1 s −1 for IDIDF and spiro-OMeTAD, respectively.These results confirm the excellent charge transport properties of the indoloindole-based HTMs for future applications in organic electronics and photovoltaics.

Figure 2 .
Figure 2. Minimum-energy crystal structure calculated at the PBEsol level of theory for IDIDF (left) and spiro-OMeTAD (right).Color coding: C in black, S in yellow, N in blue, O in red and F in cyan.H atoms are omitted for clarity.

Figure 3 .
Figure 3. Band structure calculated at the HSE06 level of IDIDF (left) and spiro-OMeTAD (right) on the PBEsol-optimized crystal geometries.Only two k-segments are displayed for simplicity (see Figure S1 for full band structure diagrams).Relevant effective masses for holes (in the VBM) and for electrons (in the CBM) are indicated.

Figure 4 .
Figure 4. Schematic representation of the main directions of hole transport in a) spiro-OMeTAD and b) IDIDF.The location of each molecule is indicated by the position of its centroid.Two colours for the centroids in spiro-OMeTAD manifest the two inequivalent molecules within the unit cell.The electronic couplings between neighbours connected by arrows are summarized.Note that the IDIDF representation requires a larger supercell due to the relevant secondneighbour contacts.Notation of ABC direction indicates the displacement along the crystallographic axes a, b and c in unit cells.
) and a  deviation slightly larger for the dimer ( = 179 meV) compared to dimer ( = 131 meV).211 100 The mean electronic coupling, however, significantly differs between the two dimers: [ = 36, V  = 31 meV] for dimer , and [ = 9,  = 8 meV] for dimer (Figure 6a).The obtained 100 V 211 V for dimers and are therefore smaller than those calculated in the static crystal (63 and 100 211 17 meV, respectively), suggesting a detrimental character of the dynamic disorder in the charge transport.Moving to spherical-like spiro-OMeTAD, we predict again a set of well-fitted Gaussian distributions for site energies in the four relevant types of dimers.In this case, despite not all distributions are centred at zero ( from −204 to 425 meV), they present similar dispersions (  Faraday Discussions Accepted Manuscript Open Access Article.Published on 21 August 2023.Downloaded on 9/23/2023 1:35:43 PM.This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.View Article Online DOI: 10.1039/D3FD00144J

Figure 6 .
Figure 6.Normal distributions for the site energy (top) and dimer electronic coupling (bottom) of IDIDF (a) and spiro-OMeTAD (b) obtained along the molecular dynamics simulations for the main charge transport directions.

Figure 7 .
Figure 7. Anisotropic plots of the hole mobility (cm 2 V −1 s −1 ) calculated for the static crystals (0 K) and considering dynamic disorder (298 K) in IDIDF (a) and spiro-OMeTAD (b) along the main crystallographic planes generated by axes a (black), b (red) and c (green).

Figure 8 .
Figure 8. a) Representative snapshot of the amorphous IDIDF phase.b) Radial distribution function (g(r)) of the centroid of a central molecule with respect to the other molecules in crystalline and amorphous IDIDF.

Figure 9 .
Figure 9. Absolute value of the electronic coupling vs the integrated volume of  associated to the weakly stabilizing non-covalent interactions for IDIDF (a) and spiro-OMeTAD (b) first nearest neighbouring dimers as obtained during the NPT molecular dynamics in the amorphous phase.The representation of the RDG index for representative dimers A, B and C of IDIDF are displayed in Figure S9.