Theoretical insights on morphology and charge transport properties of two-dimensional N,N′-ditridecylperylene-3,4,9,10-tetra carboxylic diimide aggregates

Abstract

The relationship between molecular packing and charge transport properties in bidimensional crystalline aggregates of N,N′-ditridecylperylene-3,4,9,10-tetracarboxylic diimide (PTCDI-C13), a prototypical n-type semiconductor for organic electronics, was investigated by an integrated approach based on atomistic molecular dynamics (MD) and density functional theory simulations. Calculations were performed to assess the morphological and dynamical properties of 2D crystals and nanoaggregates resulting from a set of PTCDI-C13 putative structures, differing in the relative orientation of the π-core between adjacent molecules. MD simulations indicate the cofacial crystal as the thermodynamically most stable phase and suggest the likely occurrence of competing two-dimensional structures at typical conditions for the growth of organic thin-films. In particular, the two most relevant structures found can be related to the kinetic and thermodynamic control of PTCDI-C13 crystal growth, respectively. Moreover, electronic structure calculations indicate a strong dependence of charge transport properties on molecular aggregation in the 2D crystals. Our simulations provide information, at the atomistic level, on the morphology of 2D aggregates of PTCDI-C13 in experimental conditions and shed light on the relationship between material processing and resulting charge transport properties, with a potentially strong impact in the development of devices based on PTCDI-C13.

---

Introduction

The progress of applied research in organic and hybrid electronics relies largely on the development of novel materials, with tailored properties, able to perform targeted, complex functionalities at the molecular level.¹ In this respect, intense research efforts have been carried out, during the past two decades, to develop materials based on small organic molecules, constituting the active component of devices, such as organic light-emitting diodes (OLEDs) and organic field-effect transistors (OFETs).^2,3 Indeed, in these devices the overall performances are generally related to the charge transport properties of the active material.^4,5 Beside the intrinsic properties of individual molecules, including for example electronic energy levels, charge transport in molecular materials also depends on the peculiar morphology and nanoscale aggregation realized in devices.^6,7 These factors, in turn, are related to the device architecture, processing techniques, and physico-chemical interactions with the environment and other materials, occurring for example in heterointerfaces.⁸ In particular, remarkable electronic and optical performances in OLEDs and OFETs can be achieved by active layers constituted by thin-films of crystalline or polycrystalline organic materials,^1,9,10 where the growth of ordered aggregates is a prerequisite for optimal efficiencies. The potential occurrence of polymorphism in the formation of organic crystal aggregates, as reported in several previous studies,^11–15 may impact further the overall device performances. The relationship between the directionality and anisotropy of nanoscale crystalline aggregation and charge transport properties is particularly relevant in planar devices, such as OFETs, where charge accumulation is essentially confined to the organic layers at the semiconductor/dielectric interface.^7,16,17 Accordingly, the properties of OFETs depend critically on aggregation properties and crystallinity in the pseudo-bidimensional layers constituting the active thin-film. Detailed investigations on aggregation phenomena in layered organic semiconductors and the link with charge transport properties constitute thus an essential element for the optimization of devices and the comprehensive understanding of structure/property relationships.

In this work we investigate, by means of atomistic molecular dynamics (MD) simulations, the stability, structure and dynamics of 2-dimensional (2D) crystalline phases of N,N′-ditridecylperylene-3,4,9,10-tetra carboxylic diimide (PTCDI-C13), a n-type organic semiconductor based on the perylene core.^18,19 In addition, density functional theory (DFT) calculations are performed on 2D crystal structures of PTCDI-C13 to assess the relationship between crystal packing and charge transport properties. The outstanding charge mobilities observed in devices¹⁹ are usually ascribed to the peculiar molecular structure of PTCDI-C13, constituted by a perylene diimide π-core functionalized by two extended C13 alkyl chains, which are able to induce a remarkable long-range ordered packing at the solid state.^20–24 Moreover, the interplay between the supramolecular arrangement of alkyl chains and the π–π stacking interactions among perylene diimide cores can be expected to represent the primary phenomenon behind the formation of ordered PTCDI-C13 layered aggregates, observed in several experiments.^22,25,26 Despite the large amount of work performed on the characterization of PTCDI-C13 thin films at heterointerfaces,^{23,25,27–30} the details of crystal formation and packing, at the molecular level, still remain elusive. Indeed, the typical polycrystallinity of organic layered aggregates^22,25,26 hinders the accurate characterization of in-plane packing in PTCDI-C13 thin-films by standard XRD techniques.^19,24,28,31 The most accurate available structural information to date suggests formation of ordered island of PTCDI-C13 with a typical terrace height that is compatible with molecules arranged in a stand-up configuration.^23,28 The resulting morphology of PTCDI-C13 thin-films in best-performing devices consists essentially of compact layers of quasi-2D polycrystalline aggregates,²⁹ with molecules oriented in a direction that is favourable to charge transport. However, details on the peculiar realization of the intermolecular π-stacking, a comprehensive rationalization of structural and dynamical aspects of aggregation of PTCDI-C13 layers and the relationships with growth conditions and processing are still missing. A detailed knowledge of the structure of the involved crystalline aggregates, at the atomistic level, constitutes thus a pivotal field of investigation for the optimization of materials and devices. To shed light on this issue, we investigated three putative crystal structures of PTCDI-C13 layers, differing essentially in the relative orientation of the perylene cores. The arrangement of PTCDI-C13 molecules in each of these tentative structures has been inferred from the known crystal structure of thin-films of similar organic systems. Namely, 2D crystals with perylene cores arranged in a cofacial (CF), staggered (ST) and herringbone (HB) configuration, respectively, commonly observed in perylenes, were considered.^19,32–35

Investigations on 2D structures elucidate the morphology of PTCDI-C13 in monolayers and in layer-by-layer growth, that is, in aggregation conditions that can be considered as the most favourable to charge transport. Indeed, although the nature and morphology of the substrate, as well as growth conditions and material processing, are known to greatly affect the aggregation of organic overlayers,^19,20,23,27 detailed studies on ideal 2D structures can assess the intrinsic propensity of PTCDI-C13 to form ordered layers at interfaces. The analysis of the 2D crystal structures, their morphology and dynamics, and the relationship with the resulting electronic properties allows to understand how these parameters impact on the performances of devices.

Computational methods

Molecular mechanics calculations were carried out by using the all-atom OPLS force field³⁶ augmented by parameters taken from the work of Marcon et al.³⁷ Electrostatic interactions were treated by the particle-mesh Ewald (PME) method and a cutoff of 10.0 and 14.0 Å was used for Coulomb and van der Waals interactions, respectively. The Nosé–Hoover thermostat^38,39 was used for simulations in the canonical (NVT) ensemble, and the Parrinello–Rahman barostat⁴⁰ was added in variable-cell (NPT) simulations, with time constants of 1.0 ps and 10.0 ps, respectively. Periodic boundary conditions (PBCs) in 3 dimensions were used for the simulation of 2D systems, by inserting a vacuum region (about 20 Å) along the z direction. The initial cell parameters of the CF phase (Z = 2) were taken from Tatemichi et al.¹⁹ From these parameters, a 2D supercell was built, by replicating the experimental unit cell along the [100] and [010] directions (see below). Similar model systems were built for the ST and HB phases. Structures were initially relaxed to the local energy minimum by steepest descent optimization and, in analogy with experimental processing,^{19,20,22,23,27,29} annealed by MD at 400 K for 1 ns. Annealed structures were then relaxed to 300 K (1 ns), obtaining configurations on which conformational analysis at room temperature was performed. Equilibration for longer simulation times led to converged structural properties (see ESI†). A further cooling to negligible kinetic energies and final steepest descent relaxation led to optimized structures, used to evaluate structural parameters (shown in Table 1) and as input for DFT calculations. All MD calculations were performed with the Gromacs software package.⁴¹

Table 1 Cell parameters for PTCDI-C13 2D crystals, surface densities, relative energies and transfer integral for holes and electrons computed at the DFT and Force-Field (FF) levels

Structure	a (Å)		b (Å)		γ (°)		Density (mol nm^−2)		ΔE (meV)	τh (meV)	τe (meV)
	DFT	FF	DFT	FF	DFT	FF	DFT	FF	DFT
a The total electronic energy of the CF phase is taken as reference.
CF	8.62	8.78	8.57	8.51	82.74	81.49	2.73	2.71	0^a	92	235
ST	7.36	7.48	8.65	8.62	85.94	84.72	3.15	3.06	244	44	37
HB	7.11	7.81	8.91	8.76	74.50	72.49	3.28	3.07	301	2	46

---

The polymorphs identified by MD calculations were fully optimized in a slab model with 2D periodic boundary conditions. Periodic density functional theory (DFT) calculations were carried out with the CRYSTAL09 package.^42,43 The latter has the noteworthy advantage of the absence of repetition of the slab along the z-direction that avoids artificial electrostatic interactions across the vacuum region affecting other codes. The exchange-correlation contributions to the total energy were treated using the hybrid B3LYP functional with 20% of exact Hartree–Fock exchange.^44,45 The all-electron Gaussian-type basis sets adopted were 6-31(d1) for oxygen,⁴⁶ nitrogen⁴⁷ and carbon⁴⁸ and 31(p1) for hydrogen.⁴⁶ Long-range dispersion interactions were accounted for by a post-DFT dispersive contribution, suggested by Grimme,⁴⁹ to the computed ab initio total energy and gradients. van der Waals radii and dispersion coefficients C₆ were taken from Table 1 of ref. 49. The reciprocal space was sampled by a 2 × 2 × 1 k-point mesh corresponding to 4 k-points of the irreducible Brillouin zone.⁴³ For each optimized structure, the charge transfer integrals in dimers, τ, were computed as:

where J_ij, S_ij and e_ii, e_jj are the transfer integral, the overlap integral and the site energies, respectively:

J_ij = 〈φ_i|h_ks|φ_j〉 S_ij = 〈φ_i|φ_j〉 e_ii = 〈φ_i|h_ks|φ_i〉 e_jj = 〈φ_j|h_ks|φ_j〉

with φ_i the frontier molecular orbital (i.e. HOMO or LUMO) of the i-th molecule and h_ks the Kohn–Sham Hamiltonian of the molecular dimer. In the atomic orbital basis, φ_i and h_ks were provided by the vector of molecular orbital coefficients and the Fock matrix for the converged electron density as provided by the GAMESS code (6-311G(dp) all-electron Gaussian-type basis set).^50,51 The Fock matrix can be easily obtained by:

h_ks = SVΛV⁻¹

where V and Λ are the eigenvectors and eigenvalues matrix of the converged electron density, respectively, and S is the overlap matrix printed by the GAMESS code with the 1-electron integrals (NPRINT = 3 in $CONTRL section). The use of the range-separated CAM-B3LYP functional, which can be crucial for excited-state calculations in charge-transfer complexes,^52,53 led to values of transfer integrals in line with B3LYP calculations (see ESI†). For each crystal structure, the reported τ were computed as the average of values obtained for the two non-equivalent closest dimers in the 2D system.

Results and discussion

The structures of the different 2D crystals of PTCDI-C13 considered were initially investigated by MD simulations and analyzed in terms of structural parameters. To this end, models of the three (CF, ST and HB) structures were initially prepared by inserting dimers in a monoclinic (CF, HB) and orthorhombic (ST) unit cell, as in Fig. 1, and replicated to generate 2D supercells, with a vacuum region along the z direction, containing 24 molecules each. All parameters of the two-dimensional cell (a, b, and γ) were allowed to relax.

Fig. 1 Structure of the PTCDI-C13 molecule (top) and cartoons of the initial two-dimensional crystals of the CF (bottom, left), ST (bottom, center) and HB (bottom, right) crystals.

The supercells of all 2D crystals exhibit remarkable long-range translational crystal order upon MD relaxation at room temperature (see Computational details) and can be reduced to a unit cell with Z = 2. The optimized 2D structures of PTCDI-C13 are shown in Fig. 2. Notably, the computed cell parameters of the CF and ST crystals are compatible with available experimental data.¹⁹

Fig. 2 Optimized structures of the CF (left), ST (center) and HB (right) 2D crystals of PTCDI-C13.

The computed surface densities, which can be related to molecular packing in 2D aggregates, are 2.71 mol nm⁻², 3.06 mol nm⁻² and 3.07 mol nm⁻² for the CF, ST and HB models, respectively (see Table 1). The slightly larger surface density of the ST and HB polymorphs can be ascribed to the peculiar relative orientation of the perylene cores, allowing the minimization of the steric hindrance of the alkyl chains and, in the ST crystal, to strong π-stacking interactions. Inversely, the more symmetric CF structure leads to less effective space filling and, consequently, to a lower surface density. The molecular packing in 2D crystals was further assessed by computing the radial pair distribution function g(r) between centers of mass of PTCDI-C13 molecules. The computed g(r) of the relaxed 2D crystals for models of the three structures considered is shown in Fig. 3a. The most intense peaks of the g(r) (3.75 Å, 4.55 Å and 5.85 Å for the ST, CF and HB structures, respectively) correlate with the crystal packing of PTCDI-C13 molecules along the [100] crystallographic direction, which coincides with the axis of main π–π stacking, perpendicular to the aromatic rings, for the CF and ST structures. Less intense g(r) peaks, between 7.5 and 9.5 Å, correspond to crystal packing along the [010] direction for CF and ST structures, and to a combination of [100] and [010] ordering for the more isotropic HB crystal. The preferential crystal packing of CF and ST structures along the [100] direction agrees with the observed needle-like growth of nanoscale perylene aggregates.^21,25,27 However, charge transport properties of materials based on the perylene diimide unit are also greatly affected by the relative orientation of neighboring π-cores.⁵⁴ The relative orientation of neighboring molecules in 2D PTCDI-C13 aggregates was evaluated in terms of two order parameters, defining the twist (θ) and tilt (φ) angles between the two planes of the perylene cores, respectively, as shown in the inset of Fig. 3b. According to this definition, the pair of angles [θ, φ] results about [0°, 0°], [30°, 0°] and [0°, 30°] for the CF, ST, and HB structures, respectively (see Fig. 3b). The angular structural parameters [θ, φ] can be used to identify the different PTCDI-C13 2D crystals, thus allowing a detailed analysis of morphology and structural modifications in both layers and nanoaggregates.

Fig. 3 Radial distribution function (left) and angular order parameters (right) for the three PTCDI-C13 2D crystals considered (CF, blue; ST, red; HB, green) at 300 K.

To explore the thermal stability of PTCDI-C13 crystals, NPT simulations at higher temperatures were also performed. The structure of the CF crystal is essentially unchanged by relaxation at 500 K for about 2 ns. Also, no relevant structural changes are observed upon annealing of the HB crystal at 500 K. However, variable-cell MD relaxation of the ST structure at 500 K leads to a complete phase transition to the CF crystal, suggesting a viable low-energy dynamical path for the interconversion between these two forms. To further assess the stability and the dynamical behavior of the different 2D structures and in realistic finite-size systems, MD simulations on small planar (about 3 × 3 nm) non-periodic clusters of PTCDI-C13 were also performed in non-periodic conditions. Clusters of the CF crystal were found to be stable under MD at annealing temperatures higher than 500 K on timescales of about 10 ns. However, the analysis of the structure obtained by annealing the CF cluster to 700 K reveals a still relatively ordered arrangement of PTCDI-C13 molecules, with coexistence of aggregates in liquid-crystal CF and ST forms, as shown in Fig. 4a. Likewise, annealing the HB cluster at high temperature (about 700 K) induces molecular aggregation into a similar polycrystalline structure composed of CF and ST phases (see Fig. 4b). Inversely, the ST cluster undergoes a thermally-induced rearrangement to the CF form, as observed for the periodic system, by annealing at 400 K (see Fig. 4c and Movie M1 in the ESI†), thus at temperatures lower than those generally used to anneal PTCDI-C13 thin-films.⁵⁵ Therefore, at high temperatures the free energy basins corresponding to both the CF and ST configurations are sampled, whereas relaxation toward the global minimum configuration is observed at lower temperatures. The remarkable crystal ordering at room temperature in 2D crystals, the coexistence of CF and ST liquid crystal phases at high temperatures and the abrupt interconversion from the ST to the CF structure at relatively low temperatures indicate the CF crystal as the most stable and suggest the occurrence of a relatively low energy barrier for transition from the ST metastable phase. Moreover, the HB crystal constitutes a high-energy metastable form of PTCDI-C13 and, as such, not likely to be observed in 2D aggregates.

Fig. 4 Initial (left) and final (center) structures of non-periodic PTCDI-C13 clusters and angular order parameters for the final structure (right) for (a) relaxation of the CF cluster at 700 K; (b) relaxation of the HB cluster at 700 K; (c) relaxation of the ST cluster at 400 K, inducing a phase transition to the CF phase.

Accordingly, the CF crystal is expected to dominate in the thermodynamically controlled growth of aggregates, whereas both CF and ST forms may coexist in particular environments. These findings agree with recent experiments, where formation of co-facial aggregates of PTCDI-C13 is commonly observed.^19,21,31 It must be pointed out that the energy landscape sampled by MD simulations, referring to 2D layers and clusters in vacuum, may significantly be altered upon interaction with a substrate. However, the propensity to ordered aggregation and the relative stability of the structures considered can be expected to hold also in several practical cases, for example for growth of PTCDI-C13 on smooth surfaces and for relatively weak molecule–substrate interactions.

The factors behind 2D aggregation in PTCDI-C13 crystals were further investigated by modeling nanorod-like⁵⁶ systems with periodicity in one dimension. To this end, MD simulations were performed on quasi-1D supercells of the CF and ST structures, respectively, with periodicity along the [100] direction of the corresponding 2D crystal and a variable number of replica along the orthogonal in-plane non-periodic direction. The mean potential energy at 300 K of PTCDI-C13 CF and ST aggregates composed by 1 to 7 supercells along the non-periodic direction of the system is shown in Fig. 5. For a single stack of PTCDI-C13 molecules, the ST structure is energetically more stable, as a consequence of the reduced steric hindrance of the alkyl chains with respect to the CF phase. However, the stability of the CF structure increases with the size of the system along the [010] direction of the corresponding 2D crystal, as a consequence of the more favourable 2D packing of the perylene cores. The relative stability of CF and ST quasi-1D planar structures, as a function of their width along the [010] direction, can be related to the relative propensity to form aggregates in thermodynamically- and kinetically-controlled growth conditions. Namely, in kinetically-controlled growth of PTCDI-C13 the effect of terms associated to entropy, related mainly to the conformation of flexible alkyl chains, dominates the total free energy of the system. In these conditions, aggregation of needle-like molecular stacks in the ST conformation, characterized by higher flexibility of alkyl chains and consequent lower steric hindrance, can be observed. Inversely, the thermodynamically-controlled growth of aggregates, induced by low initial molecular kinetic energy in the deposition process or by annealing/quenching cycles, is expected to lead to formation of extended 2D layers in the most stable CF form, as also suggested by the phase transition discussed above.

Fig. 5 Difference between the mean potential energies (at 300 K) of the CF and ST nanorods, with periodicity along the main axis only, as a function of the number of rows along the in-plane non-periodic direction.

To assess the energetics and the electronic properties of the PTCDI-C13 structures in more detail, DFT calculations were performed on 2D crystal systems starting from the structures obtained by MD simulations. DFT variable-cell full geometry optimizations indicate the CF crystal as the most stable, with a total electronic energy of 244 and 301 meV per cell lower with respect to the ST and HB crystals, respectively (see Table 1). The DFT optimized lattice parameters for the PTCDI-C13 crystal structures considered, also shown in Table 1, are in good agreement with force-field results, with a slight increase of surface densities. On the basis of the DFT wavefunctions, charge transport properties of the PTCDI-C13 crystals were evaluated by computing transfer integrals τ. Despite the less dense packing, the CF structure exhibits a much larger electron transfer integral (τ_e) with respect to the ST and HB 2D crystals (see Table 1). This difference can be ascribed to the more favourable overlap of LUMO orbitals of PTCDI-C13 molecules stacked along the [100] direction in the CF crystal, as shown by the orbital topologies of Fig. 6. Indeed, even small orbital misalignments are known to impact greatly on the transfer integrals in perylene diimide derivatives.⁴⁹ Inversely, the computed hole transfer integral (τ_h) for the CF structure is about three times lower than the TI for electrons and only two times larger than that of the ST crystal, confirming the low propensity to p-type transport in perylene diimide derivatives,^28,57 as also suggested by the negligible τ_h for the HB structure.

Fig. 6 Topology of the HOMO (bottom) and LUMO (top) orbitals for PTCDI-C13 dimers in the configuration of the CF (left), ST (center) and HB (right) DFT optimized crystal structures.

Conclusions

In summary, an integrated approach based on MD and DFT simulations was applied to investigate the relationship between crystal structure, stability, dynamics and charge transport properties in 2D PTCDI-C13 aggregates. MD simulations indicate a remarkable long-range stability for different putative 2D crystals. The thermodynamically most stable crystal structure corresponds to a face-to-face interaction between the π-cores of the perylene moieties, thus confirming experimental evidences. Moreover, a metastable 2D crystal phase of PTCDI-C13, with perylene cores arranged in a staggered configuration, is also found. The relative stability of the two lowest-energy 2D phases, however, is reversed in quasi-1D nanorods, suggesting a likely role of side-to-side interactions in stabilizing 2D aggregates. The local ordering of aggregates was analyzed in terms of radial and angular structural parameters. The latter, in particular, can be related to the relative orientation of the π-cores. Moreover, transfer integral calculations based on DFT suggested a strong correlation between molecular morphology and charge transport properties. Namely, the most stable 2D structure of PTCDI-C13 is found to be also the most effective for electron transport. Our results suggest a likely dependence of the molecular aggregation in nanostructures with environment, growth conditions and thermal treatments, affecting dramatically charge transport properties in devices. In particular, a controlled aggregation of PTCDI-C13 towards long-range-ordered 2D layers in the thermodynamically stable face-to-face phase is expected to lead to most efficient in-plane charge transport. Work is in progress to extend this study to model 3-dimensional aggregation and layered growth of PTCDI-C13 in hetero-organic interfaces.

Acknowledgements

We acknowledge the CINECA award under the ISCRA initiative, for the availability of high performance computing resources and support. We also acknowledge PRACE for awarding us access to HPC systems based in Poland at Cyfronet. The project has received funding from the NMP-20-2014 – “Widening materials models” program under Grant Agreement No. 646259 (MOSTOPHOS). We thank E.T.C. S.r.l., a company of the SAES Group, for the support to this work.

Notes and references

1. C. Wang, H. Dong, W. Hu, Y. Liu and D. Zhu, Chem. Rev., 2012, 112, 2208–2267.
2. P. M. Beaujuge and J. M. J. Frechet, J. Am. Chem. Soc., 2011, 133, 20009–20029.
3. M. Bendikov, F. Wudl and D. F. Perepichka, Chem. Rev., 2004, 104, 4891–4945.
4. H. Bässler and A. Köhler, Top. Curr. Chem., 2012, 312, 1–65.
5. V. Coropceanu, J. Cornil, D. A. da Silva Filho, Y. Olivier, R. Silbey and J. L. Brédas, Chem. Rev., 2007, 107, 926–952.
6. A. P. Kulkarni, C. J. Tonzola, A. Babel and S. A. Jenekhe, Chem. Mater., 2004, 16, 4556–4573.
7. A. Lorenzoni, M. Muccini and F. Mercuri, RSC Adv., 2015, 5, 11797–11805.
8. J. Wünsche, G. Tarabella, S. Bertolazzi, M. Bocoum, N. Coppede, L. Barba, G. Arrighetti, L. Lutterotti, S. Iannotta, F. Cicoira and C. Santato, J. Mater. Chem. C, 2013, 1, 967.
9. M. Muccini, Nat. Mater., 2006, 5, 605–613.
10. A. Facchetti, Mater. Today, 2007, 10, 28–37.
11. Y. Diao, K. M. Lenn, W. Lee, M. A. Blood-forsythe, J. Xu, Y. Mao, Y. Kim, J. A. Reinspach, S. Park, G. Xue and P. Clancy, J. Am. Chem. Soc., 2014, 136, 17046–17057.
12. A. Troisi and G. Orlandi, J. Phys. Chem. B, 2005, 109, 1849–1856.
13. Y. Diao, B. C.-K. Tee, G. Giri, J. Xu, D. H. Kim, H. a. Becerril, R. M. Stoltenberg, T. H. Lee, G. Xue, S. C. B. Mannsfeld and Z. Bao, Nat. Mater., 2013, 12, 665–671.
14. L. Maini, F. Gallino, M. Zambianchi, M. Durso, M. Gazzano, K. Rubini, D. Gentili, I. Manet, M. Muccini, S. Toffanin, M. Cavallini and M. Melucci, Chem. Commun., 2015, 51, 2033–2035.
15. T. Zykova-Timan, P. Raitei and M. Parrinello, J. Phys. Chem. B, 2008, 112, 13231–13237.
16. W. W. A. Koopman, S. Toffanin, M. Natali, S. Troisi, R. Capelli, V. Biondo, A. Stefani and M. Muccini, Nano Lett., 2014, 14, 1695–1700.
17. M. A. Alam, A. Dodabalapur and M. R. Pinto, IEEE Trans. Electron Devices, 1997, 44, 1332–1337.
18. R. J. Chesterfield, J. C. McKeen, C. R. Newman, P. C. Ewbank, D. a. Da Silva Filho, J. L. Brédas, L. L. Miller, K. R. Mann and C. D. Frisbie, J. Phys. Chem. B, 2004, 108, 19281–19292.
19. S. Tatemichi, M. Ichikawa, T. Koyama and Y. Taniguchi, Appl. Phys. Lett., 2006, 89, 2004–2007.
20. Y. Zhang, H. S. Seo, M. J. An and J. H. Choi, Org. Electron., 2009, 10, 895–900.
21. A. L. Briseno, S. C. B. Mannsfeld, C. Reese, J. M. Hancock, Y. Xiong, S. a. Jenekhe, Z. Bao and Y. Xia, Nano Lett., 2007, 7, 2847–2853.
22. H. G. Jeon, N. Oguma, N. Hirata and M. Ichikawa, Org. Electron., 2013, 14, 19–25.
23. J. Jang, S. Nam, D. S. Chung, S. H. Kim, W. M. Yun and C. E. Park, Adv. Funct. Mater., 2010, 20, 2611–2618.
24. J. Puigdollers, M. Della Pirriera, a. Marsal, a. Orpella, S. Cheylan, C. Voz and R. Alcubilla, Thin Solid Films, 2009, 517, 6271–6274.
25. C. Rolin, K. Vasseur, S. Schols, M. Jouk, G. Duhoux, R. Müller, J. Genoe and P. Heremans, Appl. Phys. Lett., 2008, 93, 033305.
26. a. S. Wan, G. P. Kushto and A. J. Mäkinen, J. Vac. Sci. Technol., A, 2010, 28, 1275.
27. K. Vasseur, C. Rolin, S. Vandezande, K. Temst, L. Froyen and P. Heremans, J. Phys. Chem. C, 2010, 114, 2730–2737.
28. Y. J. Jeong, J. Jang, S. Nam, K. Kim, L. H. Kim, S. Park, T. K. An and C. E. Park, ACS Appl. Mater. Interfaces, 2014, 6, 6816–6824.
29. M. Jang, K. Y. Baek and H. Yang, J. Phys. Chem. C, 2013, 117, 25290–25297.
30. A. S. Wan, J. P. Long, G. Kushto and A. J. Mäkinen, Chem. Phys. Lett., 2008, 463, 72–77.
31. Y. Zhang, H.-S. Seo, M.-J. An and J.-H. Choi, Org. Electron., 2009, 10, 895–900.
32. C. W. Struijk, a. B. Sieval, J. E. J. Dakhorst, M. Van Dijk, P. Kimkes, R. B. M. Koehorst, H. Donker, T. J. Schaafsma, S. J. Picken, a. M. Van de Craats, J. M. Warman, H. Zuilhof and E. J. R. Sudholter, J. Am. Chem. Soc., 2000, 122, 11057–11066.
33. K. A. Kistler, C. M. Pochas, H. Yamagata, S. Matsika and F. C. Spano, J. Phys. Chem. B, 2012, 116, 77–86.
34. R. Schmidt, J. H. Oh, Y. Sen Sun, M. Deppisch, A. M. Krause, K. Radacki, H. Braunschweig, M. Könemann, P. Erk, Z. Bao and F. Würthner, J. Am. Chem. Soc., 2009, 131, 6215–6228.
35. T. N. Krauss, E. Barrena, X. N. Zhang, D. G. de Oteyza, J. Major, V. Dehm, F. Würthner, L. P. Cavalcanti and H. Dosch, Langmuir, 2008, 24, 12742–12744.
36. W. L. Jorgensen, D. S. Maxwell and J. Tirado-Rives, J. Am. Chem. Soc., 1996, 118, 11225–11236.
37. V. Marcon, D. W. Breiby, W. Pisula, J. Dahl, J. Kirkpatrick, S. Patwardhan, F. Grozema and D. Andrienko, J. Am. Chem. Soc., 2009, 131, 11426–11432.
38. S. Nosé, J. Chem. Phys., 1984, 81, 511–519.
39. W. G. Hoover, Phys. Rev. A, 1985, 31, 1695–1697.
40. M. Parrinello, A. Rahman and R. a. Parrinello, Phys. Rev. Lett., 1980, 45, 1196–1199.
41. D. Van Der Spoel, E. Lindahl, B. Hess, G. Groenhof, A. E. Mark and H. J. C. Berendsen, J. Comput. Chem., 2005, 26, 1701–1718.
42. R. Dovesi, R. Orlando, R. Civalleri, C. Roetti, V. R. Saunders and C. M. Zicovich-Wilson, Z. Kristallogr.–Cryst. Mater., 2005, 220, 571.
43. R. Dovesi, V. R. Saunders, C. Roetti, R. Orlando, C. M. Zicovich-Wilson, F. Pascale, B. Civalleri, K. Doll, N. M. Harrison, I. J. Bush, P. D'Arco and M. Llunell, CRYSTAL09 User's Manual, University of Torino, Torino, 2009.
44. A. D. Becke, J. Chem. Phys., 1993, 98, 5648.
45. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785–789.
46. M. Corno, C. Busco, B. Civalleri and P. Ugliengo, Phys. Chem. Chem. Phys., 2006, 8, 2464–2472.
47. C. Gatti, V. R. Saunders and C. Roetti, J. Chem. Phys., 1994, 101, 10686.
48. M. A. Spackman and A. S. Mitchell, Phys. Chem. Chem. Phys., 2001, 3, 1518–1523.
49. S. Grimme, J. Comput. Chem., 2006, 27, 1787–1799.
50. M. W. Schmidt, K. K. Baldridge, J. A. Boatz, S. T. Elbert, M. S. Gordon, J. H. Jensen, S. Koseki, N. Matsunaga, K. A. Nguyen, S. Su, T. L. Windus, M. Dupuis and J. A. Montgomery, J. Comput. Chem., 1993, 14, 1347–1363.
51. H. Tamura, I. Hamada, H. Shang, K. Oniwa, M. Akhtaruzzaman, T. Jin, N. Asao, Y. Yamamoto, T. Kanagasekaran, H. Shimotani, S. Ikeda and K. Tanigaki, J. Phys. Chem. C, 2013, 117, 8072–8078.
52. W. T. Peng and J. D. Chai, Phys. Chem. Chem. Phys., 2014, 1–15.
53. B. M. Wong and T. H. Hsieh, J. Chem. Theory Comput., 2010, 6, 3704–3712.
54. F. May, V. Marcon, M. R. Hansen, F. Grozema and D. Andrienko, J. Mater. Chem., 2011, 21, 9538.
55. S. Tatemichi, M. Ichikawa, S. Kato, T. Koyama and Y. Taniguchi, Phys. Status Solidi RRL, 2008, 2, 47–49.
56. F. Mercuri, D. Costa and P. Marcus, J. Phys. Chem. C, 2009, 113, 5228–5237.
57. J. Roh, J. Lee, C. M. Kang, C. Lee and B. Jun Jung, Phys. Status Solidi RRL, 2013, 7, 469–472.

---

Footnote

† Electronic supplementary information (ESI) available: Convergence of the computed radial distribution function with the MD equilibration time. Comparison between transfer integrals computed with the CAM-B3LYP and B3LYP functionals. Movie of the MD trajectory of the phase transition between two PTCDI-C13 crystal structures. See DOI: 10.1039/c6ra06784k

---