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 .
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 .
the list of organic films displaying polymorphism 
 and sips has grown 
 considerably in length . 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 , shows thin film phases , some of them 
 with hexaphenyl molecules arranged perpendicularly to the herringone 
 layers .
on the other hand , rubrene ( 5,6,11,12-tetraphenyltetracene ) 
 bulk crystal polymorphs have been obtained by changing growth conditions or ambient pressure , but crystalline 
 thin films show invariably the orthorhombic structure .
pentacene shows a rich bulk polymorphism based on different 
 triclinic 
 crystal structures in space group p1 and characterized 
 by herringbone packings of molecules .
the first triclinic phase ( polymorph i ) has been described in 1961 and appeared again only recently . for many years
thin films also show a variety of polymorphs , some of them displaying 
 small differences in their diffraction pattern , quite often coexisting 
 as interpenetrating microdomains or at 
 different depth in the films .
growth parameters or substrate nature such as sio2 , silicon , h - terminated or 
 passivated silicon , metals , alkali metal halides , and polymer surfaces provide a thin film 
 phase with upright molecules possibly with polymorph ii co - present . 
 moreover , an orthorhombic phase has been observed in low thickness films grown on 
 sio2 or polymers . 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 . 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 .
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 ,
while for inorganic surfaces 
 a rich literature and a surface structure database are available , reports on surface reconstructions of molecular 
 solids are sparse .
we are aware of only a few examples concerning pyrene 
 ( 001 ) , benzyl ( 001 ) , and benzophenone surfaces , 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 . the recent study of morisaki et al .
, 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 ) . 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 .
recent works exploring this 
 approach include studies of growth of thin films and epitaxial structures , bulk crystals , and polymorphism . 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 .
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 can not 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 .
the crystal structure of nnn has been solved 
 recently for a thin film phase grown on thermally oxidized silicon .
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 .
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 1 .
this belongs to 
 the monoclinic p21/a space 
 group , with two molecules per unit cell .
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 .
left : unit cell content of nnn thin film and 
 bulk phase ( room temperature ) . 

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 .
the last three columns collect unit cell volumes , 
 densities , and tilt angles (  ) .
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 

the 
 parameters for bonding and nonbonding interactions were equal to the 
 standard ones for aromatic carbons and hydrogens .
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 . 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 .
for the latter , we performed constrained geometry optimizations 
 using the b3lyp density functional with 
 the 6 - 311 g * * basis set and an empirical dispersion correction ( b3lyp - d3/6 - 311 g * * ) , followed by an analogous 
 single point calculation with a larger basis including diffuse functions 
 ( b3lyp - d3/6 - 311++g * * ) .
it is qualitatively 
 similar to that of biphenyl , but it 
 is less symmetric due to the lower symmetry of the naphthyl units .
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 ) .
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 3 . 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 .
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 .
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 .
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 , 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 ) .
electrostatic interactions were treated via 
 the particle - mesh - ewald method 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 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 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 .
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 .
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
table 1 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)p / v  180  10 k is comparable with that of other organic semiconducting 
 materials , such as naphthalene ( 251(9 ) 
  10 k ) and anthracene ( 181(3 )  10 k ) .
further increasing the temperature 
 to 600 k promotes a structural 
 rearrangement associated with a solid - to - solid phase transition .
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
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 / cm ) is 
 about 10% smaller than the experimental one at 300 k. accordingly , 
 the distance between adjacent ( 001 ) planes increases from 1.939 to 2.077(7 )  .
evolution of the  
 angle of all molecules belonging to the 
 bulk phase at 600 k over the first nanosecond .
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 ) .
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 2 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 , 

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 1 ) .
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 1 ) .
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 4 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 ,
after 1 ns , all molecules have changed 
 their orientation ( see structure b in figure 4 ) .
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 1 . to further confirm 
 the stability of the new structure
, we also performed a long npt run ( 5 ns ) using the parrinello 
in contrast to berendsen s , this barostat correctly samples 
 the canonical ensemble but requires a well equilibrated system to 
 avoid numerical instabilities .
also in 
 this case , the parameters collected in table 1 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 ) .
the presence of sharp peaks at both short and long 
 range at all temperatures confirms the absence of liquid crystalline 
 phases .
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 5 also compares the torsional angle distributions 
 ( 1 and 2 ) at 550 and 600 k , over 
 the first nanosecond of the simulation .
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. 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 ) . 

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
figure 3 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 6 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 , where surface 
 relaxation was observed in the first monolayer of single crystals .
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 6 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 4 ) .
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 7 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 .
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 .
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 8 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 6 .
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 .
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 6 ) 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 9 , 
 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 2 ) . unit cell content 
 for the bulk at 600 k ( a ) and ( 001 ) s0 surface 
 at 550 k ( b ) .
to describe the 
 relative positions of the molecules and the structure 
 of the defects , we performed a voronoi analysis .
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 10 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 10 ) , 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 11 ) seem to have a more dynamic , amoeba - like 
 structure .
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 .
representative 
 voronoi diagrams for the s2-dist surface at 550 
 k. dots represent the molecules centers of mass .
representative voronoi diagrams for the s2-prox surface at 550 
 k. dots represent the molecules centers of mass .
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 , 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 , along 
 with a characterization of the change in electronic properties ( e.g. , 
 charge carrier mobilities ) 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 .
 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
Model and Computational Methods
Results and Discussion
Summary and Conclusions