Morphology and Electronic Properties of N,N′-Ditridecylperylene-3,4,9,10-tetracarboxylic Diimide Layered Aggregates: From Structural Predictions to Charge Transport

The morphology of layered aggregates of N,N′-ditridecylperylene-3,4,9,10-tetracarboxylic diimide (PTCDI-C13), a prototypical n-type semiconductor for organic electronic devices, was investigated by molecular dynamics and corroborated by metadynamics simulations. Calculations were targeted to ordered 3D aggregates, differing in the relative orientation of the perylene π-cores and on the degree of interdigitation among contiguous planar layers. Our simulations indicated the noninterdigitated cofacial structure as the thermodynamically most stable form of ordered PTCDI-C13 aggregates, in both bulk crystals and bilayers. Other structures, however, may occur in the growth of PTCDI-C13 under kinetic conditions. Density functional theory calculations were also performed to evaluate the relative total electronic energy of 3D crystals of PTCDI-C13 and related transfer integrals, correlating structure with potential charge-transport properties in devices. The most stable ordered aggregated form of PTCDI-C13 exhibit...


Introduction
Recent progress in organic electronics underlines the need for the development of novel materials with tailored structure at the nanoscale and exceptional charge transport properties. [1][2][3][4] Particularly, organic materials that are able to aggregate in layered stacks of molecules have been spotted as ideal candidates to achieve long-range in-plane crystal ordering, leading to remarkable conductivities in devices. 1,5,6 Among these, small molecules based on functionalized perylene cores, as perylene-tetracarboxylic diimides (PTCDI), display significant electron mobilities, as a result of low LUMO energy, electron delocalization of unoccupied levels and efficient packing in thin films. 6,7 Although the energy of the LUMO level is generally associated to low air stabilities of thin-films based on the PTCDI core, ambient stability of devices can be achieved by co-evaporation of dopants or by polymer interfacial layers. 8,9 The possibility of tuning the functionality of PTCDI side chains has also led to the synthesis of molecular systems where the charge transport properties of the perylene core is coupled to structural and interface properties. 10,11 Namely, the N,N'-ditridecylperylene-3,4,9,10-tetracarboxylic diimide (PTCDI-C13), where the planar core of PTCDI is functionalized with C13 alkyl chains, has been used with success in devices for electronics and bioelectronics applications. 8,[12][13][14][15] The molecular structure of chain-substituted perylene molecules, with an interplay between the rigid π-core and flexible alkyl chains, leads to peculiar morphologies in nanostructures, usually in the form of stable layered aggregates, where packing is driven by π-π stacking interactions. 16 However, the morphology of thin-films formed by molecules based on the PTCDI core can also be affected by the substituents as well as by nature of the substrate, growing conditions and thermal treatments. 17 This variation in morphology may dramatically influence the properties of thin-film in devices, for example transport of charge carriers. 18,19 Structural properties of PTCDI-C13 thin-films have been investigated on a variety of substrates and deposition conditions by scanning probe and X-ray diffraction (XRD) techniques. 20,21 Nevertheless, the intrinsic polycrystalline nature of aggregates and the alterations of morphology with film thickness prevent a precise determination of the packing of PTCDI-C13 in devices at the molecular level. 22,23 In this work, we investigate the structural properties of bulk and bilayer ordered aggregates of PTCDI-C13 by molecular dynamics (MD) simulations. The relationship between morphology and charge transport is evaluated by computing electron transfer integrals at the density functional theory (DFT) level. In particular, we relax putative structures of 3D aggregates of PTCDI-C13 by atomistic MD. Relaxed structures are used to carry out metadynamics 24 calculations in the space of suitable collective variables, identifying global and local minima and computing the related free energy. The structures of the predicted stable and metastable configurations are also used to evaluate total energies and electron transfer integrals in 3D ordered aggregates at the DFT level. The effect on the charge transport properties of the disorder induced by interactions of the organic layer with substrates, as in the case of organic/dielectric interfaces, is further assessed by evaluating transfer integrals of aggregates of PTCDI-C13 in contact with a model surface. It is however worth noting that, although the combination of MD and metadynamics calculations is potentially able to provide information on the free energy surface of molecular aggregates, kinetic aspects of nucleation and growth processes in organic materials and the dynamics of thin-film formation are not explicitly considered in the proposed approach.
As we observed in previous work, 23 the PTCDI-C13 molecule exhibits two low-energy packing configurations in planar two-dimensional (2D) aggregates, with the perylene cores arranged according to a symmetric cofacial (CF) or staggered (ST) configuration, respectively. Moreover, aggregates of molecules with a planar π-core and long tail substituents may exhibit various degrees of interdigitation of alkyl chains, which assists formation of regularly packed structures. 7,16,25 We therefore focused our simulations on bulk and bilayer aggregates of the CF 2D crystal of PTCDI-C13, which is the thermodynamically most stable structure at ambient conditions, considering both interdigitated and non-interdigitated configurations. Furthermore, calculations were also performed on the 3D aggregation of PTCDI-C13 in the ST configuration, which is likely to occur in kinetically-controlled growth conditions. 23 The integration between atomistic MD, metadynamics simulations, DFT total energy calculations, evaluation of electronic transfer integrals and the comparison with available experimental data, guides us in assessing the structure of PTCDI-C13 in nanoaggregates as a function of growing conditions, correlating molecular structure to potential performances in devices.

Theoretical Methods
Molecular mechanics calculations were performed by applying a customized version of the allatom OPLS force-field, 19,26 which was successfully applied in previous work on PTCDI-C13 aggregates. 23 Electrostatic interactions were evaluated by the particle-mesh Ewald (PME) method, using cut-off of 10.0 Å for both Coulomb and van der Waals interactions. The Nose-Hoover thermostat 27,28 was used for simulations in the NVT ensemble, with a time constant of 1.0 ps, and the Parrinello-Rahman barostat 29 added in variable-cell (NPT) simulations, with a time constant of 10.0 ps. Periodic boundary conditions (PBCs) in three dimensions were used for the simulation of all systems. For finite-size systems (e.g. bilayers), a vacuum region of about 50 Å was inserted along the z direction of the box.
All structures were initially relaxed to the nearest energy minimum by steepest descent optimization and subsequently equilibrated at 300 K, in the NPT ensemble, for 4 ns. Stable configurations were cooled to negligible kinetic energies and re-optimized by steepest descent. These minimum-energy configurations were used for further DFT calculations. The effect of thermal treatments on PTCDI-C13 aggregates was simulated by MD annealing from room temperature to 500 K, followed by 2 ns equilibration.
Free energy calculations were performed by metadynamics, 24 using relevant structural parameters as collective variables, as detailed below, adding Gaussian functions to the metadynamics potentials with a frequency of 1000/step, a height of 10 -3 kJ mol -1 and a width sigma of 0.10 kJ mol -1 . MD calculations were performed with the Gromacs (version 5.1) software package, 30 whereas LAMMPS package 31 was used for metadynamics simulations, including the Plumed plugin. 32 The free energy profile evaluated by metadynamics was also compared with the free energy obtained by steered MD simulations, using the same collective variable as pulling coordinate. For either MD and metadynamics simulations the same force-field and simulation parameters were used.
Periodic DFT calculations on bulk model systems were carried out with the SIESTA 33 package (version 4.0), by applying the PBE 34 gradient-corrected approximation to DFT, using a basis set of double-ζ plus polarization quality for orbital expansion and a cut-off of 150 Ry for density expansion. Van der Waals interactions were accounted for by including a Grimme 35 dispersion potential. A reciprocal space mesh of 3x3x1 k-points was used in periodic calculations. Geometries were optimized by relaxing all the nuclear degrees of freedom until convergence to a maximum force of 0.02 eV/Å on atoms and 1 GPa on cell stress.
Transfer integrals were evaluated on dimers extracted from the DFT periodic calculations described above using the charge transfer integrals module of the ADF program package. 36 Electronic couplings for electrons and holes were computed by the PW91 37 exchange-correlation functional and a triple-ζ plus-polarization basis set. This set-up is found to provide results that are in closest agreement with experimental results. 38 The distribution of intermolecular transfer integrals in PTCDI-C13 was also computed for layered three-dimensional crystalline aggregates at the interface with model surfaces. To this end, a periodic model surface was first generated by constructing a 5x5 nm regular square-lattice of dummy atoms with 0.1 nm spacing, and assigning a random value, within selected boundaries, for the z (i.e., orthogonal to the molecular plane) component of the coordinates. The z component was then changed iteratively by averaging the first derivative of neighbouring sites, until convergence to a given total surface root mean square (RMS) roughness. This algorithm allows to obtain a model surface with selected RMS and higher-order moment surface parameters (skewness and kurtosis) depending on the initial random distribution of z coordinates. In agreement with the value commonly observed on the surface of polymer or inorganic dielectric materials for in-plane devices, 39-41 a model surface with a RMS roughness of 0.3 nm was generated, with relatively large skewness and kurtosis values (0.098 and -0.414, respectively) to enhance the effects of local molecule/surface interactions. All atoms of the model surface were assigned non-bonded parameters of graphitic carbon atoms, to mimic the effect of hydrophobic substrate/molecule weak forces, and were kept frozen throughout all simulations. A supercell of the multi-layer crystalline PTCDI-C13 aggregate, with lateral dimension of about 5x5 nm, was first relaxed in vacuum, then put in close contact with the model surface, and subsequently equilibrated at room temperature by MD for 10 ns. This simulation time was found to be sufficient for the formation of a stable interface between the PTCDI-C13 aggregate and the model surface. Transfer integrals were computed on the equilibrated system for all molecule pairs. The effect of annealing on the morphology of PTCDI-C13 aggregates on a model surface was further evaluated by analyzing relevant structural parameters upon annealing at 400 K followed by relaxation for 5 ns.

Results and Discussion
Simulation models were prepared basing on the configurations of the perylene diimide cores taken from equilibrated 2D (monolayer) CF and ST layer phases, as described in ref. 23. Initially, the flexible C13 alkyl chains of PTCDI-C13 were arranged to form, in the direction orthogonal to that of π-stacking, a fully-interdigitated 3D bulk crystal. Model systems were relaxed at room temperature in the NPT ensemble, as described in the method section. The degree of interdigitation was progressively decreased, by changing the cell parameter orthogonal to the π-stacking plane, relaxing again the system by variable-cell MD. Simulations converged to four stable configurations at room temperature, corresponding to a fully-interdigitated and a layon structure for both the CF and the ST crystals, shown in Fig. 1. Similar configurations were also observed by equilibration of ordered bilayers of PTCDI-C13 molecules, as shown in Fig. 2. In bilayers, however, the topmost and bottommost alkyl chains of PTCDI-C13 lean towards the π-cores, in analogy to what is observed in monolayers. 23 The tilted configuration of the C13 chains can thus be considered a structural feature of the topmost layers in ordered PTCDI-C13 aggregates, corresponding to the surfaces exposed to scanning-probe techniques, or at heterointerfaces. Notably, the computed inter-layer distance for the CF lay-on bulk configurations (25.2 Å) agrees remarkably well with the height of the PTCDI-C13 unit cell (25.3 Å) obtained by XRD experiments. 20,21,42 In addition, the slight increase of the inter-layer distance in bilayers (25.6 Å) is compatible with the larger d-spacing observed in thin-films (25.6 -26.1 Å). 20,42 Unlike other substituted perylene derivatives, exhibiting interdigitated stable structures, 7 the relatively short inter-layer distance in PTCDI-C13 thin-films, with respect to the end-to-end length of a relaxed isolated molecule (about 43 Å), can be related to the tilted configuration assumed by molecules on the xy plane. To assess the relative stability of the relaxed structures, and to explore the underlying free energy landscape, further MD simulations were performed. Namely, the equilibrated 3D bulk structures were annealed at high temperatures, as described in the methods section. No relevant structural changes were observed for the CF lay-on structures in the timescale of MD simulations. However, annealing of the relaxed ST lay-on and interdigitated structures at 400 K for 2 ns led to interconversion to the CF lay-on configuration. The interconversion mechanism from ST lay-on and interdigitated structures to the CF lay-on configuration suggests a larger stability for this latter, in complete analogy with the case of PTCDI-C13 monolayers. 23 Further assessments of the equilibrated structures of Fig. 1 and quantitative estimates of the related Gibbs free energy surface were carried out by predictive structural simulations using metadynamics. To this end, the CF interdigitated bulk structure of PTCDI-C13, which, for the above considerations, is expected to constitute a local free energy minimum, was used as the starting configuration for metadynamics. The system was thus propagated with the history-dependent potential (as detailed in the methods section), using the zcomponent of the distance (dz) between the terminal carbon atoms of contiguous alkyl chains as a collective variable (see Fig. 3). The defined collective variable accounts for structural changes from interdigitated to lay-on configurations. The resulting Gibbs free energy surface, as a function of the selected collective variable (see Fig. 3), confirms the occurrence of two main basins of attraction, corresponding to the fully-interdigitated (dz=1.15 nm) and to the lay-on (dz=-0.13 nm) CF configurations, respectively. The CF lay-on structure results more stable than the corresponding fully-interdigitated configuration, with a difference in free energy of about 20 kJ mol -1 at 300 K and a barrier for interconversion that is estimated to be about 45 kJ mol -1 . The main contribution to the observed free energy difference can be ascribed to the coulombic repulsion between alkyl C13 chains, which is significantly larger in the interdigitated configuration. Therefore, from our observations, alkyl chains assist interlayer ordered packing, despite a low propensity to interdigitation. Remarkably, the broader free energy basin corresponding to the interdigitated configuration lies closely to a local minimum, located at dz = 0.83 nm. This minimum corresponds to the partial slipping of interdigitated PTCDI-C13 molecules from each other, sterically controlled by the hindrance of -CH2-and -CH3 groups. The transition state connecting the interdigitated and lay-on structures is located at an intermediate value (dz = 0.37 nm) on the computed free energy surface, with a structure corresponding to an incompletely-interdigitated configuration. A movie representing the transition from the interdigitated to the lay-on structure, extracted from the metadynamics trajectory, is provided in the Supporting Information. Steered-MD simulations of the interconversion between the interdigitated and lay-on forms of PTCDI-C13 confirm essentially the same underlying free energy landscape (see Fig. S1 in the Supporting Information).
More detailed information about the structure and energetics of PTCDI-C13 crystals were obtained by DFT calculations. The structural parameters for the configurations in Fig. 1, optimized at the DFT level and corresponding total DFT electronic energies are reported in Table  1. Calculations indicate the CF lay-on configuration as the most stable, with the ST lay-on configuration lying about 44 kJ mol -1 above in energy. Moreover, interdigitated crystals are generally less stable than the corresponding lay-on structures by more than 65 kJ mol -1 . Table 1. Lattice parameters of PTCDI-C13 crystals (a,b,c: Å; α,β,γ: degrees), computed at the DFT level, and total electronic relative energies (kJ mol -1 ). Therefore, it can be concluded that the CF lay-on configuration corresponds to the global energy minimum configuration for PTCDI-C13 crystals. This picture is confirmed by the excellent agreement between the DFT computed interlayer distance, corresponding to the c value of the cell parameter (25.9 Å) and the value observed by XRD and AFM experiments (between 25.6 and 26.1 Å) in PTCDI-C13 aggregates. 20,42 Other computed CF lay-on lattice parameters are in agreement with experiments. 21 To shed light on the relationship between morphology and charge transport properties of PTCDI-C13, the intermolecular transfer integrals were computed on nearest neighbour dimers extracted from DFT-optimized structures and from equilibrated MD configurations. The resulting transfer integrals are reported in Table 2. Table 2. Intermolecular transfer integrals in meV for electrons (coupling of LUMO orbitals) and holes (coupling of HOMO integrals) for dimers extracted from the DFT optimized PTCDI-C13 crystal structures (a), average values for dimers of the PTCDI-C13 MD crystal structure at 300 K (b), first layer (c) and second layer (d) at the interface with a model surface. The proposed multiscale approach is able to guide us to a better understanding of the structure/properties relationship in molecular aggregates of organic materials. In particular, starting from the molecular modelling of 3D bulk and layered phases, we are able to correlate morphology to charge transport through electronic structure calculations. This, in turn, allows us to speculate on the electronic behaviour of materials in thin-film devices. In the case of crystalline aggregates of PTCDI-C13, the computed values of transfer integrals were found in line with previous calculations on substituted PTCDI systems. 23,43,44 Despite the relatively large transfer integrals, the poor hole conducting properties of PTCDI derivatives are usually ascribed to the deep ionization potential. 19 As expected, computed transfer integrals suggest significant charge transport in the direction of the π-stacking (values shown in Table 2) for crystalline aggregates. The CF structures exhibits comparable transfer integrals for electrons in the lay-on and interdigitated configurations. Although less stable, the interdigitated structures were found to shows higher values of transfer integrals with respect to the lay-on ones. Transfer integrals for the ST structures are generally lower than for the CF configurations, as a result of the tilted πstacking, similarly to the case of monolayers. 23 Beside crystalline ideal structures, the effect of more realistic environments on charge transport properties was assessed by computing the distribution of transfer integrals for a bilayer of the most stable PTCDI-C13 phase (CF lay-on) equilibrated at room temperature and at the interface with a model surface with a morphology resembling that of a typical substrate for the growth of organic thin-films. As expected, thermal disorder broadens the distribution of transfer integrals for both electrons and holes, 45 also slightly lowering the average value by about 22 and 6 meV for electrons and holes, respectively.  (Fig. S2).
The distribution of transfer integrals for holes is further broadened and shifted to lower values for an ordered bilayer of PTCDI-C13 aggregates relaxed on surfaces (See Fig. 4). However, electron transfer integrals are essentially unchanged with respect to the crystal for the first layer of the PTCDI-C13 in contact with the surface, whereas lower values are observed for the second layer. Indeed, in the first PTCDI-C13 layer ordered packing is maintained by the interplay between strong π-stacking interactions and flexibility of alkyl chains, which induce a bulk-like interlayer interaction and adapt to the morphology of the underlying surface. Moreover, in the particular case of interaction with a model surface, the morphology of PTCDI-C13 aggregates is slightly affected by annealing at 400 K, as suggested by computed structural parameters (see Fig.  S3 in Supporting Information). Despite the slight lowering of intermolecular transfer integrals related to thermal and structural disorder, computed values suggest remarkable charge transport properties for multi-layer aggregates of PTCDI-C13 in realistic environments. In particular, transport of electrons exhibits a lower propensity to degradation upon structural modifications with respect to holes.

Conclusions
In summary, an integrated computational approach, based on MD and DFT calculations, allowed to understand the details of the morphology in PTCDI-C13 aggregates, to evaluate the related stability and energetics and to correlate structure to charge transport properties. The most stable ordered bulk crystal structure found corresponds to a CF lay-on configuration. Computed crystal parameters are in excellent agreement with available structural data. The CF lay-on configuration is also prevalent in bilayers of PTCDI-C13, mimicking the situation occurring in thin-films or at interfaces. Our calculations also indicate local energy minima for interdigitated (CF and ST) and for asymmetric ST structures, where the π-cores of perylene diimide arrange in pairs of rotated units. These configurations, however, are more than 40 kJ mol -1 higher than the global energy minimum. Indeed, metadynamics calculations and simulation of heat treatments suggest a likely interconversion path from local minima configurations to the thermodynamically stable CF layon structure. In particular, the formation of interdigitated aggregates of PTCDI-C13 is found to be unfavored with respect to lay-on structures, with C13 alkyl chains assisting intra-layer aggregation. Nevertheless, the possible occurrence of different structures in nanoscale aggregates of PTCDI-C13 confirms experimental evidence, where an interplay between layered and threedimensional growth is observed, as a function of processing conditions. Moreover, the thermodynamically stable CF lay-on configuration exhibits remarkable charge transport properties, as suggested by intermolecular charge transfer integral calculations. A broadening of the distribution of transfer integral values for electrons and holes and a lowering of the average is observed for PTCDI-C13 layers equilibrated at room temperature and at the interface with a model substrate. However, computed values suggest remarkable charge transport properties for the CF lay-on configuration also in partially disordered morphologies, suggesting this configuration as responsible of the observed high conductivities in layered aggregates of PTCDI-C13.

Supporting Information.
Comparison between the free energy computed by metadynamics and steered-MD. Convergence of metadynamics with respect to simulation time. Comparison between structural parameters of unannealed and annealed PTCDI-C13 aggregates on a model surface. Distributions of transfer integrals for electrons and holes of PTCDI-C13 relaxed structures (PDF). A movie extracted from the metadynamics trajectory (MPG).