Surface Reconstructions in Organic Crystals: Simulations of the Effect of Temperature and Defectivity on Bulk and (001) Surfaces of 2,2':6',2″-Ternaphthalene.

2,2':6',2″-Ternaphthalene (NNN) is a novel, blue-emitting material, suitable for preparation of organic light-emitting diodes. Its crystal structure has been solved recently, but its thermal behavior and surface properties have not yet been explored, partly due to the difficulty in obtaining high quality crystals. In the present study we use classical molecular dynamics to investigate the thermal behavior of bulk and (001) surfaces of NNN. Our bulk simulations indicate the occurrence of a phase transition at about 600 K. The transition is facilitated by the presence of a free (001) surface, since a reconstruction leading to a very similar structure occurs around 550 K in our surface models. This holds for both ideal and defective surface models, containing a small number of vacancies (one or two missing molecules in the outermost layer). In all cases, the process is triggered by thermal motion and involves the reorientation of the molecules with respect to the (001) plane. Both the bulk and surface phases share the monoclinic space group P21/a with a herringbone disposition of molecules. These findings and their implications for the use of NNN in organic electronics are discussed.

Introduction

Years of intensive research on organic electronic devices have taught us that supramolecular interactions and organization have a major influence on their final performance. Connected to the preparation of devices, e.g., by vacuum deposition techniques or with solution methods, the sword of Damocles of polymorphism is always present, definitely a widespread phenomenon when dealing with molecular crystals.[1] Polymorphism in bulk phases is a general issue in several fields, including pharmaceuticals and materials science. In thin films—the suitable form for a realistic exploitation of electronic materials—occurrence of polymorphism is further enhanced by the presence of the substrate surface and/or strong kinetic effects during growth. When these polymorphs are distinct from those observed in bulk crystals, one may speak of surface-induced phases (SIPs).

Since the first discovery of polymorphism in pentancene thin films.[2,3] the list of organic films displaying polymorphism and SIPs has grown considerably in length.[4,5] A brief, nonexhaustive, summary of SIPs limited to the most promising compounds containing only C and H atoms already contains several paradigmatic examples. p-Hexaphenyl (4′,4″,4‴,4′′′′-hexaphenyl), although scarcely prone to giving polymorphic bulk structures,[6,7] shows thin film phases,[8] some of them with hexaphenyl molecules arranged perpendicularly to the herringone layers.[9,10] On the other hand, rubrene (5,6,11,12-tetraphenyltetracene) bulk crystal polymorphs have been obtained by changing growth conditions[11−14] or ambient pressure,[15] but crystalline thin films show invariably the orthorhombic structure.[16−27]

Pentacene shows a rich bulk polymorphism based on different triclinic crystal structures in space group P1̅ and characterized by herringbone packings of molecules.[28−33] The first triclinic phase (polymorph I) has been described in 1961[34,35] and appeared again only recently.[28,36−39] For many years only polymorph II was obtained in single crystal form by different growth techniques.[28−31,33] Thin films also show a variety of polymorphs, some of them displaying small differences in their diffraction pattern, quite often coexisting as interpenetrating microdomains[40] or at different depth in the films.[2,3,41] Growth parameters or substrate nature[42] such as SiO2,[43,44] silicon,[45] H-terminated[42] or passivated silicon,[46] metals,[47−51] alkali metal halides,[52−55] and polymer surfaces[56−58] provide a thin film phase with upright molecules possibly with polymorph II co-present. Moreover, an orthorhombic phase[59−62] has been observed in low thickness films grown on SiO2 or polymers.[63]

Up to now, there is no general or clear-cut agreement about the causes triggering the formation of SIPs. They may arise from kinetic phenomena under the influence of the substrate or of the strain induced by lattice mismatch at the film/substrate interface. Sometimes, a thickness dependent behavior has been observed due to the formation of metastable phases.[65] For pentacene, a clear thickness dependence of polymorphs is found with films grown on oxide substrates: the single crystal phase appears beyond a critical film thickness, and the critical thickness decreases on increasing the substrate temperature.[66−68]

Molecular crystals are soft materials with weak molecular interactions. Large-amplitude molecular motions about the average positions may become possible on going from the constrained environment of the bulk to the surface, possibly triggering surface relaxation/reconstruction, even at low temperatures. Despite these characteristics, judging from the scientific literature, molecular crystals usually undergo very small relaxation with respect to the bulk. While for inorganic surfaces a rich literature and a surface structure database are available,[69] reports on surface reconstructions of molecular solids are sparse.[70−72] We are aware of only a few examples concerning pyrene (001),[73] benzyl (001),[74,75] and benzophenone surfaces,[76] while recent molecular mechanics minimizations (i.e., excluding thermal disorder) on organic semiconductors crystals exhibiting the typical herringbone packing showed that surface relaxation of the low energy crystal faces is in general negligible.[77−80] The recent study of Morisaki et al.,[72] conversely, has shown that surface relaxation in tetracene single crystals may have a deep impact on electron transport properties.

In this respect, a deep study of the structure of molecular crystal surfaces is a fundamental step, especially when studying technologically relevant crystalline films. Structural properties of thin films and film/substrate interfaces can be studied experimentally by ex situ and in situ scattering techniques (grazing incidence X-ray diffraction, X-ray reflectivity, and electron diffraction).[43,81−83] In parallel, computer modeling of crystals and crystal surfaces has substantially grown in importance over the past decade, thanks to more reliable force fields and increased computer power. Molecular dynamics (MD) simulations have provided a deeper understanding of the interactions and dynamical processes in organic semiconductor thin films, sometimes with a level of detail hardly accessible by experimental investigations.[83] Recent works exploring this approach include studies of growth of thin films and epitaxial structures,[84−93] bulk crystals,[94] and polymorphism.[95,96] Particularly pertinent to this work, SIPs have also become the subject of MD simulations centered on the role of substrate surface on polymorph selection during thin film deposition.[97−99]

Owing to the fact that polymorphism can be triggered by several mechanisms and parameters, in our opinion intrinsic properties of surfaces such as relaxation/reconstruction cannot be excluded a priori as a factor responsible for the appearance of polymorphs in thin films. The present study is devoted to 2,2′:6′,2″-ternaphthalene (NNN), a molecule suitable for preparation of organic light-emitting diodes.[100−102] The crystal structure of NNN has been solved recently for a thin film phase grown on thermally oxidized silicon.[82] The same crystal structure has been observed on single crystals obtained by vacuum sublimation. Unfortunately, the low quality of the single crystals hampered an accurate comparison of the two data sets, looking for minor differences in their crystal packing.

This work aims to set up a background knowledge about surface features of NNN (001), the low energy face displayed by thin films, by means of classical MD. The role of MD in this work is to investigate the stability of the ideal (001) surface and of surfaces containing point defects (molecule vacancies), looking for possible surface reconstruction phenomena linked to molecular thermal motion. We deem this preliminary step necessary for a full understanding of possible surface induced polymorphism or phase transformations associated with thin film phases of NNN grown on different substrates, by organic molecular beam deposition or hot wall epitaxy. Indeed, our MD simulations of the NNN (001) surface give evidence of a rich structural behavior, with intramolecular effects coupled to intermolecular interactions and temperature as a determining variable for triggering surface reconstruction.

Model and Computational Methods

To derive models of the NNN (001) surface, we started with the unit cell parameters and crystal structure published in ref (82) and illustrated in Figure . This belongs to the monoclinic P21/a space group, with two molecules per unit cell. The experimental lattice parameters are given in Table . The molecules are arranged in layers, parallel to the ab plane. Within a stack, they adopt a classic herringbone arrangement, their long molecular axes being tilted by 23° with respect to the layer normal.

Figure 1

Left: unit cell content of NNN thin film and bulk phase (room temperature). Right: molecular structure and descriptors of NNN molecules during the MD simulations. The modulus of θ is the tilt angle with respect to the surface normal (c*). The sign of θ is the same as the inner product between the molecular axis and the a lattice vector. τ1 and τ2 are dihedral angles defined by the atoms in violet; τ1 is always on the outermost side of the molecule.

Table 1

Temperature Dependence of the Unit Cell Parameters of NNN, As Obtained from the NPT MD Simulations of the Crystalsa

T (K)	a (Å)	b (Å)	c (Å)	α (deg)	β (deg)	γ (deg)	V (Å3)	d (g/cm3)	θ (deg)
300b	8.148	5.978	19.452	90.00	94.58	90.00	944.4	1.314	–23.0
300	8.197(10)	5.990(7)	19.903(19)	90.00(1)	92.58(16)	90.00(1)	976.1	1.294	–18.97(86)
400	8.270(10)	6.010(10)	19.986(24)	90.00(1)	91.98(21)	90.00(1)	993.0	1.274	–18.47(85)
500	8.358(11)	6.028(13)	20.100(32)	90.00(2)	91.12(31)	90.00(1)	1012.0	1.248	–17.72(85)
550	8.412(13)	6.036(16)	20.189(41)	90.00(2)	90.40(42)	90.00(1)	1025.0	1.233	–17.07(86)
600c	9.048(15)	5.658(5)	20.952(5)	90.00(2)	97.43(3)	90.00(2)	1063.0	1.188	8.62(40)
600d	9.084(40)	5.656(20)	20.957(20)	90.00(17)	97.38(15)	90.00(11)	1070.0	1.180	8.77(18)

Standard deviations are also given in parentheses. The last three columns collect unit cell volumes, densities, and tilt angles (θ). Volumes and densities are calculated for the final trajectory. The experimental values[82] at room temperature are also given for comparison.

Experimental data.

Bulk values collected after the transition at 600 K, extending the first simulation by 1 ns.

Bulk values collected after the transition at 600 K, extending the first simulation by 5 ns, and replacing the Berendsen with the Parrinello–Rahman barostat.

A force field for the NNN molecule was developed from OPLS-AA.[103] Except for inter-naphthalene torsions, the parameters for bonding and nonbonding interactions were equal to the standard ones for aromatic carbons and hydrogens. The equilibrium length of the inter-naphthalene C–C bonds was increased to 1.48 Å in order to match the average value from ab initio geometry optimizations (details are given later), and the stretching and bending force constants involving them were reduced by 30% with respect to the standard OPLS-AA ones, following the prescriptions of DuBay et al.[104] for conjugated oligomers and polymers. Finally, the inter-ring intrinsic torsion potential was obtained from the energies of the conformers of 2,2′-dinaphthalene (NN). In practice, we fitted the molecular mechanics energies to those calculated ab initio with the GAMESS-US code.[105] For the latter, we performed constrained geometry optimizations using the B3LYP density functional[106] with the 6-311G** basis set and an empirical dispersion correction[107] (B3LYP-D3/6-311G**), followed by an analogous single point calculation with a larger basis including diffuse functions (B3LYP-D3/6-311++G**). The resulting torsion potential is shown in Figure . It is qualitatively similar to that of biphenyl,[108] but it is less symmetric due to the lower symmetry of the naphthyl units.

Figure 2

Plot of the potential energy curve obtained by fitting B3LYP-D3/6-311++G** data (squares), for the torsion about the central bond of 2,2′-dinaphthalene. The structures superposed on the plot depict the molecule in the transoid (at 38.7°) and the cisoid (at 138.0°) conformation.

A large supercell consisting of 10 × 10 × 8 unit cells along the a, b, and c axes (1600 molecules overall) was constructed and employed to describe the bulk crystal by MD, with three-dimensional periodic boundary conditions. To model the (001) surface, the supercell periodicity was kept within the ab plane, while, along c, it was elongated from ≈15.5 to 30 nm, thus creating a slab with a wide vacuum space between its two sides. In one set of simulations, the molecules occupying the outermost layer on the lower side of the slab were left untouched, leading to an ideal and defect-free surface (excluding possible surface reconstructions and thermal disorder), while one molecule from the layer on the upper side of the slab was deleted to create a point defect (vacancy). Henceforth, we shall refer to them as “S0” and “S1” surfaces (corresponding to 100% and 99.5% of fractional coverage, respectively). In another set of simulations, we created two vacancies on both sides, by deleting two neighboring (top layer) or two far away molecules (bottom layer). Henceforth, we shall refer to them as “S2-prox” and “S2-dist” surfaces (both corresponding to 99% of fractional coverage). An example of such surfaces, obtained from simulations performed at 500 K, is shown in Figure . In all cases, these surfaces are separated by a sufficiently thick layer of bulk-like crystal (on one side) and of vacuum space (on the other side), so that they are effectively noninteracting and independent.

Figure 3

Snapshots of NNN slabs, from the MD simulations at 500 K. Left: full S0/S1 slab. Here the lower surface (S0) is defect-free, while the upper one (S1) has one vacancy. Right: the S2 slab showing two neighboring vacancies on the upper surface (S2-prox) and two distant vacancies on the lower one (S2-dist). To better visualize the defects, the inner layers have been made transparent, and the molecules surrounding the defects have been highlighted.

All MD simulations were performed using the GROMACS 4.6.5 package.[109] The simulations of the bulk crystals were carried out under constant temperature and pressure conditions (NPT), with a standard duration of 1 ns. The leapfrog algorithm with a time step of 1 fs was used to integrate the equations of motion. The pressure was maintained constant (1 atm) and controlled via anisotropic coupling to the Berendsen barostat, with a coupling constant of 1.0 ps. We chose the Berendsen, rather than the Parrinello–Rahman,[110] barostat due to its better numerical stability. The Parrinello–Rahman barostat (with a coupling costant of 10 ps) was adopted in one NPT simulation to check the stability of the new crystal structure (see later). The temperature was controlled by coupling with a velocity rescaling thermostat,[111] with a time constant of 0.1 ps.

After equilibration of the bulk crystal at each temperature, the c axis of the supercell was increased to 30 nm to create the slabs and one or two molecules at the outer surfaces were removed to create the vacancies, as explained previously. The MD simulations of the slabs were then carried out with fixed lattice parameters, under constant temperature and volume conditions (NVT). Each run had a duration of 500 ps. All short-range nonbonded interactions were cutoff at 1.2 nm. Electrostatic interactions were treated via the particle-mesh-Ewald method[112] with a Fourier grid spacing of 0.12 nm. As detailed in later text, for the bulk structure at 600 K, a phase transition was observed after about 500 ps. In order to equilibrate the final structure, the NPT simulation was extended to 2 ns. For this system, we developed two sets of slabs. The former was generated starting from a trajectory taken at the early stage of the simulation (200 ps), whereas, the latter was obtained after equilibration.

Atomic coordinates from the MD simulations were saved once every picosecond, to be visualized with VMD[113] and analyzed with our own programs to extract information about the structural organization at the NNN (001) surfaces. Relevant observables included the distributions of average distances and angles between the chosen molecules and the mean (001) plane, and intramolecular dihedral angles. See the next section for the definition of these geometric parameters. Voronoi diagrams were also calculated from the center-of-mass of the molecules in the top and the bottom layers, so as to highlight local changes in the structure, especially around the defects.

Results and Discussion

In order to investigate NNN bulk dynamics and obtain suitable starting structures for the surface simulations, we performed preliminary NPT simulations of the bulk crystal phase at 300, 400, 500, 550, and 600 K. Each of these had a standard duration of 1 ns, enough to observe the fast and relatively minor changes in the crystal structure due to thermal expansion. For the NPT simulation at 600 K, the overall simulation time was extended to 2 ns. Table compares the experimental (at room temperature) and the calculated average lattice parameters. In order to compare these parameters with those of the crystallographic unit cell, the average axes lengths of the MD supercell were divided by 10, 10 and 8, respectively. At 300 K, the simulated lattice parameters are in good agreement with the experimental ones, demonstrating the reliability of the force field used in the simulations. The calculated cell parameters change gradually and regularly with temperature up to 550 K. The thermal expansion coefficient, αV = (∂V/∂T)/V ≈ 180 × 10–6 K–1 is comparable with that of other organic semiconducting materials, such as naphthalene[114] (251(9) × 10–6 K–1) and anthracene[115] (181(3) × 10–6 K–1). Both the unit cell volume and the crystal density change linearly with temperature.

Further increasing the temperature to 600 K promotes a structural rearrangement associated with a solid-to-solid phase transition. This is illustrated in Figure (see later text for details). The formation of the new crystal phase was observed starting at 0.5 ns and took the remaining simulation time to complete. In order to calculate the structural properties of the new polymorph, we extended this run by 1 ns. The cell parameters extracted from the extended simulation are reported in Table . Compared to that at 550 K, the new cell shows longer a and c axes and a shorter b axis. The value of β increased from 90.40° to 97.43°, whereas no change was observed for α and γ. The cell volume increase across the phase transition is about three times that expected from thermal expansion only. The density of the new structure (1.188 g/cm3) is about 10% smaller than the experimental one at 300 K. Accordingly, the distance between adjacent (001) planes increases from[82] 1.939 to 2.077(7) Å.

Figure 4

Evolution of the θ angle of all molecules belonging to the bulk phase at 600 K over the first nanosecond. Deeper colors indicate higher densities in the angle distribution. The starting (A) and final (B) structures are also shown on the right.

We characterized the temperature dependence of the NNN crystal structure through the values of the angle defining the molecular orientation with respect to the (001) plane (θ) and the intramolecular dihedral angles determining their conformational state (τ1 and τ2). These structural descriptors are defined in Figure . In the crystal, the NNN molecules are always planar on average. This planarity is a result of intermolecular interactions, as the most favorable conformation for an isolated molecule would be distorted from planarity by about ±40° (see again Figure for the conformational energy). It is worth mentioning that, in addition to the transoid conformation, in the gas phase there is also a cisoid conformation (at about ±140°) with nearly identical energy. At ambient temperature and pressure, there are no molecules in such conformation within the bulk crystal. NNN molecules pack more efficiently when they are planar due to the π–π intermolecular interactions, which provide about 18 kJ/mol (=2 × 9 kJ/mol, since in each NNN molecule there are two such torsions such as those in Figure ). Nonetheless, the MD simulations clearly indicate that even at room temperature there can be appreciable instantaneous deviations from planarity. These deviations reach ±30° at 300 K and, as expected, increase with temperature (±50° at 550 K). Also the orientation of the molecular axis with respect to the ab plane, expressed by the angle θ, changes with temperature, from about −19° at 300 K to about −17° at 550 K (see Table ). This is consistent with the increase of the c axis length discussed earlier.

The phase transition observed at 600 K occurs with a significant molecular rearrangement and involves a substantial change in the average value of θ, with only a moderate broadening in the torsional angles. The entire process is induced by thermal motion and occurs gradually across adjacent molecular layers. Figure shows the evolution of θ over time. During the simulation, the value of θ, starting at −23°, gradually increases on the average up to about −15°. At this point, corresponding to 550 ps, the molecules in a single molecular layer change simultaneously their orientation with respect to the ab plane, from negative values to positive values, centered at about 8°. Shortly after, the same process is observed for an adjacent molecular layer. After 1 ns, all molecules have changed their orientation (see structure B in Figure ). It is worth noting that, during the transition, the molecules quickly swap from negative to positive θ angles, thereby suggesting an energy maximum—actually, a transition state—at θ = 0. Torsional angles (not shown for brevity) are only marginally involved in the transition because the molecules in the new structure remain essentially planar. Nonetheless, during the transition, the maximum instantaneous deviations from planarity increase to ±60°. Starting at 0.7 ns, a small fraction of molecules exhibit a change in conformation, from the trans to a cisoid conformation, with τ values around ±160° (in the bulk, there is no distinction between τ1 and τ2). The extension of this simulation by 1 ns revealed no further change in the structural parameters. These are reported in Table . To further confirm the stability of the new structure, we also performed a long NPT run (5 ns) using the Parrinello–Rahman barostat. In contrast to Berendsen’s, this barostat correctly samples the canonical ensemble but requires a well equilibrated system to avoid numerical instabilities.[94] Also in this case, the parameters collected in Table show no evidence of further structural changes.

To further characterize the structural properties of the new phase, we calculated the radial distribution functions (RDFs) between the molecular centers of mass at different temperatures (300, 550, and 600 K). These are shown in Figure . The presence of sharp peaks at both short and long range at all temperatures confirms the absence of liquid crystalline phases. The data at 600 K were obtained from the 5 ns NPT run. Increasing the temperature from 300 to 550 K broadens the peaks without significantly affecting their positions. The RDF profile for the new phase at 600 K clearly evidences a different distribution in peaks’ location and intensity. Figure also compares the torsional angle distributions (τ1 and τ2) at 550 and 600 K, over the first nanosecond of the simulation. The distribution obtained at 600 K with the Parrinello–Rahman barostat is equivalent to that observed right after the phase transition (see above). Also in this case, the torsional angles never exceeded ±60° for the majority of molecules. The number of molecules with significant changes in conformation was small (one to five molecules) compared to the whole system (1600 molecules), and within the range of the statistical variability to be expected at this temperature. We, therefore, exclude the possibility of a plastic crystal phase at 600 K.

Figure 5

Plot of the pairwise radial distribution functions between the molecular centers of mass at 300, 550, and 600 K (left panel) for the bulk phase. On the right: evolution of the dihedral angles (τ1 and τ2) at 550 (top) and 600 K (bottom). The simulation at 600 K was performed with the Parrinello–Rahman barostat for 5 ns after the phase transition. For comparison purposes, only the first nanosecond is shown.

We now turn to the discussion of the NNN slabs exposing the (001) surfaces. We are especially interested in the temperature dependence of their structure, as we are looking for evidence of single molecule rearrangements or collective surface reconstructions triggered either by thermal motion—as observed earlier for the bulk structure—and/or defects. At 300, 400, and 500 K we do not see any major structural changes with respect to the bulk. Figure shows that, even at 500 K, the structure of the outer molecular layers is closely analogous to that of the bulk crystal. On the contrary, Figure demonstrates that at 550 K there is a clear reorientation of all of the molecules belonging to the outer layers, retaining at the same time a high degree of order. A similar outcome has recently been reported for tetracene,[72] where surface relaxation was observed in the first monolayer of single crystals.

Figure 6

Snapshots from the final trajectories of the MD simulations of the NNN slabs at 550 and 600 K. The snapshots at 600 K refer to slabs generated from the corresponding bulk structure after 200 ps (e.g., prior the phase transition).

The bottom panel of Figure shows the snapshots at 600 K for the slabs obtained from the bulk prior the phase transition (i.e., at 200 ps). These trajectories are closely similar to those obtained for the bulk at the same temperature (i.e., structure B in Figure ). A first indication of this analogy is given by comparing the distances between adjacent (001) planes, equal to 2.077(7) Å for the bulk and 2.082(16) Å for the S0/S1 slab. Below, we shall show that the primitive cells extracted from the bulk and the slab share the same space group.

At 600 K, the molecular rearrangement starts from the outer layers and propagates to the whole interior of the slab. These rearrangements occur on all of the surfaces, demonstrating that they are independent of the presence and type of vacancies. The main difference between the bulk and the slabs is the time required for the transition to occur. The process took about 500 ps to begin in the bulk, but only 50 ps in the slabs. Also in this case, to characterize the structural changes at the surfaces, we calculated the inter- and intramolecular angles θ, τ1, and τ2. At relatively low temperatures (see Figure for 500 K), these parameters for the molecules within the S0 and S1 surfaces fluctuate in a relatively narrow and sharply defined range, similar to those for the bulk at the same temperature. On the contrary, the molecules belonging to the two S2 surfaces have more orientational and conformational freedom, as evidenced by the broader distributions of their θ and τ angles. The S2-dist surface appears to be more disordered than the S2-prox one. However, even in these cases the majority of molecules remains relatively unaffected. Thus, one or two vacancies produce only a local perturbation in the film structure, without triggering any collective rearrangement.

Figure 7

Evolution of the orientations (θ, left-hand side) and of the torsion angles (τ1 and τ2, right-hand side) of the molecules belonging to the surface layers, from the MD simulations at T = 500 K.

At 550 K, all of the slab systems undergo a very fast structural change. Figure shows that the time scale for these rearrangements is always of the order of 100 ps, independently of defectiveness of the surface. The structural change is clearly highlighted by the θ angle, which switches from −15° to +10°. The molecules in the outer layers switch from a tilted to a more vertical arrangement in the opposite direction, as already highlighted in Figure . This rearrangement does not involve a major change in conformation, as most molecules remain essentially planar. However, after the transition has occurred, we do see a broadening of the torsional distributions on the S0 surface, as the more vertical orientation of the molecules allows greater librational freedom.

Figure 8

Evolution of the orientations (θ, left-hand side) and of the torsion angles (τ1 and τ2, right-hand side) of the molecules belonging to the surface layers, from the MD simulations at T = 550 K.

On the three defective surfaces, a small fraction of molecules also shows a change in conformation, from transoid to cisoid (τ1, τ2 ≈ ±160°). We have observed that this occurs especially when the molecules’ ends temporarily emerge from the surface and, as a consequence, the outer naphthyl group is freed from the packing restrictions imposed by the neighboring molecules.

At 600 K, the dynamics of surface reconstruction occurs on a shorter time scale (about 50 ps) than at 550 K. The further propagation of the transformation wavefront to the slab interior (Figure ) takes about 170 ps to complete in the S0/S1 case and only 110 ps in the S2 one. In order to test the reversibility of this process, we performed an additional NVT simulation by cooling the S0/S1 slab from 600 to 300 K. Using a linear temperature gradient, and an overall annealing time of 10 ns, the original tilted herringbone structure was not restored, demonstrating a certain irreversibility for this transformation.

After averaging over space and time the MD configurations for the bulk at 600 K and (001) surface S0 at 550 K, it was possible to obtain the unit cell content after the phase transformation. Both new phases, shown in Figure , share the monoclinic space group P21/a with a herringbone disposition of molecules, differing from the starting one for the length of the unit cell axes and, more relevant, the molecule tilt with respect to the (001) plane, as evidenced by the variation of the monoclinic β angle in the three cases (see Table ).

Figure 9

Unit cell content for the bulk at 600 K (A) and (001) S0 surface at 550 K (B).

Table 2

Comparison of the Experimental Cell Parameters of NNN[82] and Those of the New P21/a Cells Obtained from the Symmetry Analysis of the MD Trajectories

T (K)	a (Å)	b (Å)	c (Å)	α (deg)	β (deg)	γ (deg)	V (Å3)	d (g/cm3)
300a	8.148	5.978	19.452	90.00	94.58	90.00	944.4	1.314
Ab	9.056	5.657	20.945	90.00	97.41	90.00	1064.0	1.187
Bb	8.405	6.031	21.240	90.00	90.34	90.00	1076.6	1.174

Experimental data.

See Figure .

To describe the relative positions of the molecules and the structure of the defects, we performed a Voronoi analysis.[116] Given a set of sites on a plane (in this case, the centers of mass of the molecules at a surface), a cell is assigned to each site, such that the points within a cell are closer to that site than to any other. Figures and 11 show some representative Voronoi diagrams for the surfaces containing two defects, at 550 K. Far away from the defects, the molecules always have a well defined hexagonal close packing associated with the herringbone motif. In the S2-dist case (Figure ), the defects have a well defined structure, defined by six neighboring molecules. We do observe some variation in the defect structure shortly after the reconstruction, as exemplified by the diagram at 150 ps, but after a while this disorder disappears and the original defect pattern is restored. Note that, by the end of the simulation, the defects have moved one step closer to each other. This demonstrates that at this temperature there is some defect mobility, which may eventually allow their coalescence. Two neighboring defects (Figure ) seem to have a more dynamic, amoeba-like structure. Also in this case, we observe some limited diffusion. A full clarification of these points and a quantitative description of defect mobilities would require further, much longer simulations and a dedicated analysis, which are currently under way.

Figure 10

Representative Voronoi diagrams for the S2-dist surface at 550 K. Dots represent the molecules’ centers of mass. The defects have been highlighted in yellow.

Figure 11

Representative Voronoi diagrams for the S2-prox surface at 550 K. Dots represent the molecules’ centers of mass. The defects have been highlighted in yellow.

Summary and Conclusions

We have applied classical MD to characterize the structural properties and thermal behavior of bulk and (001) surfaces of NNN. Our simulations on the bulk structure have revealed the occurrence of a crystal-to-crystal transition at high temperatures, which was previously undiscovered, possibly due to the difficulty in growing high quality crystalline samples of this material. The entire process is triggered by molecular thermal motion, and it occurs on a nanosecond time scale at 600 K. It involves the reorientation of the molecular long axes with respect to the (001) plane, by about 30°. Although showing significant instantaneous deviations from planarity, dihedral angles are only marginally involved in phase transition, as the majority of molecules remains essentially planar.

Simulations performed on slab models of the (001) surface with zero, one, or two vacancies have evidenced surface reconstruction phenomena starting at 550 K. At this temperature, the molecular rearrangement is characterized by the fast (100 ps) reorientation of all molecules belonging to the outer layers, again with a change in the molecules’ tilt angle by about 30°. Further increasing the temperature to 600 K, we observe the propagation of this rearrangement to the slab interior, leading to a structure similar to that of the bulk crystals at the same temperature. Increasing surface defectiveness accelerates the process, but it has a negligible effect on the final structures. The fluctuation of the torsional angles with temperature are similar to those observed for the bulk, with the more defective surfaces showing increased orientational and conformational freedom. The analogy between bulk and slab structures after phase transition is further enforced by the fact the unit cells extracted from our simulations share the monoclinic space group P21/a and have similar lattice parameters, with a herringbone disposition of molecules.

Similar to the results of a recent experimental and computational study of tetracene,[71] surface reconstruction appears to be an intrinsic property of NNN crystals, not induced by a specific foreign substrate. These findings prompt us to undertake in the future further computational work devoted to the MD simulation of true thin films on different substrates such as graphite (0001) and amorphous SiO2,[82] along with a characterization of the change in electronic properties (e.g., charge carrier mobilities[117]) associated with these reconstructions. These studies, together with the present one, will allow us to decouple effects due to the substrate and kinetics effects during thin film growth from features which are inherent in the thermodynamic behavior of the crystalline phase of NNN, correctly assigning the role of the underlying substrate.