Cocrystal Growth in Organic Semiconductor Thin Films: Simulation of Pentacene, Perfluoropentacene, and Their 1:1 Blend Deposited On Graphite

The understanding of crystal formation in thin films and the precise knowledge of the relation between structure and surface diffusion are two important requirements for the efficient (nano)fabrication of organic electronic devices. Here a computational approach for simulating vapor‐phase deposition is employed to obtain and investigate three types of crystalline thin films on graphite. All systems, namely pentacene, perfluoropentacene, and their 1:1 blend, which forms an alternate cocrystal, are constituted by recumbent molecules in accordance with experimental findings. The contributions of intermolecular interactions and of molecular rearrangements occurring during the deposition are analyzed to rationalize the final morphologies. Then, the generated structures are employed to evaluate the energy barriers that prevent molecular diffusion at terraces and step‐edges, and to study the reorganization of the films upon high‐temperature annealing. The broad agreement with experimental observations and the possibility of evaluating the potential energy surface at the molecular detail render the proposed approach a promising tool to make predictions for other systems.


Introduction
The demand for ecological and efficient electronical devices is becoming more pressing over time and organic ones still represent promising substitutes to their inorganic counterparts. Their main advantages include lower costs, mechanical flexibility, and higher margins in the modulation of the molecular functional groups, which impacts the final performance of the device, at least for some applications. Among the most proficient ones, we mention thin film devices such as organic solar cells, [1] transistors, [2] diodes for lightning and displays. [3] Despite many DOI: 10.1002/adts.202300080 advantages, their minor efficiency compels for continuous improvements, both in the search of new powerful materials [4] and in the study of their chemical-physical properties. In this sense, it is known that the magnitude and anisotropy of charge mobility across organic films is related to molecular packing and orientation, and for molecular systems it increases with the crystallinity of the sample. [5] Thus, manufacturing techniques heavily focus on the optimization of this aspect, taking also advantage from thermally activated diffusion (e.g., annealing), with contrasting and temperature-dependent results, [6,7] probably because experiments are typically conducted on devices and the morphology changes affect multiple parameters (e.g., contact resistance). [8] Concerning thin films fabrication methods, vapor-phase deposition is one of the most extensively exploited techniques: several variants [9,10] of this technology have been developed over the years, making it a useful and flexible tool for a wide range of applications, also suitable for depositing amorphous organic blends [11] or crystalline films, [12] thus going beyond its traditional application to inorganic materials. Among the several opportunities that this approach offers, the formation of organic semiconductors mixed crystals is a promising field that still needs to be explored, since their many proven advantages, such as the improved charge transport, and the large amount of experimental research on the realization of such devices. [13] In this context, atomistic simulations already proved to be a valid option for the prediction and the understanding of the underlying mechanisms of the on-surface crystal growth. [14][15][16][17][18][19][20][21][22] Concerning the deposition of mixed systems, although simulations of solid mixtures of organic molecules are becoming increasingly popular, in particular for better understanding the functioning of organic LEDs and solar cells at the nanoscale, [23,24] we are not aware of any theoretical investigation on the formation of organic semiconductor mixed crystals from vapor phase.
Pentacene (PEN) was one of the first organic semiconductors employed as active material for the construction of organic transistors, [25,26] and, despite being largely investigated throughout the years, it is still considered as a worth prototype for fundamental studies of thin films growth and as a model p-type semiconductor. An extensive number of studies also concerns PEN derivatives, where the award for the most prominent one is held by perfluoropentacene (PFP). Due to the inductive effect of fluorine atoms over the aromatic rings, and to the consequent variations of the electrostatic interactions in the crystal-which have a strong impact on transport levels-PFP shows (with respect to PEN) a change from p-type to n-type charge carrier behavior, which makes these molecules employable together in efficient p-n heterojunctions. [27][28][29][30][31][32][33][34][35] Interestingly, PFP seems to present a different crystal structure compared to PEN, as it was found in experimental depositions on graphite: while at least at low temperature PEN molecules organize similarly as in their bulk phase with a herringbone configuration with the long molecular axes parallel to the surface, [36] PFP molecules appear flat-lying on the graphite surface. [37] Among the experimental depositions of 1:1 blends of the two molecules, it was reported the formation of a mixed crystal on both graphite [32] and SiO 2 . [29,38] On the former substrate, thin films of recumbent molecules were reported, like in the case of pure PFP, while on the latter crystallites are found in two different orientations: at low temperature (250 K) the molecules are again flat-lying on the surface, while higher temperature favors their vertical configuration. [29][30][31][32][33][34][35][36][37][38] While the connection between the type of molecule and its crystalline configuration appears to be established for pure systems-although very sensible to the interaction with the substrate and to temperature-here we aim to elucidate the different parameters influencing the final outcome and to demonstrate the possibility of reproducing the formation of a PFP-PEN cocrystal by means of computer simulations. Also, knowing the importance of a well-organized morphology in the fabrication of organic thin films, [39] we aim to prove the soundness of our approach in making predictions for further similar systems and to gain deeper understanding of specific crystal growth mechanisms, in particular for mixtures. [40] Through a tailored workflow of molecular dynamics (MD) calculations, we simulated three types of vapor-phase deposition experiments on graphite: one for PEN molecules, one for PFP molecules and one for a 1:1 blend of the two species (MIX), obtaining crystalline films for all systems. The discussion of the simulation results is divided into two main sections. In the first, we focus on the molecular aggregation, outlining the specific intermolecular interactions between the deposited molecules and evaluating the general stability of the crystal in terms of cohesion energy. In the second, it is studied the presence of diffusive energy barriers across the crystalline film that hinder the molecular diffusion and the film reorganization upon annealing.

Results and Discussion
In crystalline films on smooth and defect-free graphite surfaces, the molecular long axes of PFP and PEN are expected to be parallel to the substrate, while heating or the presence of structural defects [41] that increase surface roughness favor an organization with an upright molecular orientation, typical of acenes on other rough surfaces (e.g., SiO 2 ).  In fact, while the chemical similarity between the deposited molecules and graphite promotes a recumbent alignment, an upright one allows exposing to air a lower energy surface. The unbalance between the two energies and kinetic factors determines the outcome; however, since at low submonolayer coverages the recumbent orientation is always thermodynamically favorite, for subsequently obtaining the upright one the temperature should be high enough to allow the initially recumbent molecules to overcome a collective energy barrier for reorientation. [43] We decided to keep the system at a relatively low temperature (240 K) for this reason and, more importantly, to make the diffusive processes slower and thus less relevant in the growth mechanism of the crystals. In this way, we aim to the formation of stable and long-living aggregates and avoid the complete molecular desorption from the surface (that we rarely observed even at higher temperatures). Furthermore, considering that low temperatures have already been employed for successful experimental depositions in similar systems, [44,45] we retain that this choice somehow increases the realism of the simulation protocol, that, given the necessarily limited size of the graphite surface, underestimates the importance of on-surface diffusion to the growth.
While each of the simulated depositions produced a different crystalline morphology through a peculiar growth mechanism, Figure 1 shows that the simulation attempts were broadly successful in reproducing the experimental results: all deposited molecules are found in a recumbent disposition with the long molecular axis parallel to the graphite plane and are organized in stacked monolayers (MLs), clearly outdistanced and identifiable.

Thin Film Growth
Before comparing the three systems, we start by discussing the peculiarity of the growth mechanism of each sample specifically and its connection with the adsorption energies and interaction energies between molecular pairs. Figure 1a-c shows the MLs formation as a function of the number of deposited molecules for PEN, PFP, and MIX system. The second pentacene monolayer (ML2) starts to form only near the completion of ML1 (see Figure 1a), when it becomes difficult for the incoming molecules to find enough space to descend on the surface. This behavior resembles a Frank-van der Merwe (layer-by-layer) growth type, [46] where the assembly of crystalline layers happens one at a time. With increasing coverage and thickness, the importance of the interaction with graphite fades out and the growth mechanism evolves toward a 3D one.
A different situation is found for the PFP system, where the ML2 formation begins sooner than for PEN. The PFP-PFP couple exhibits a stronger -stacking interaction, allowing the falling molecules to stick permanently above the deposited ones, also owing to the low deposition temperature. Furthermore, from the visual observation of the sample growth, we could mark out an interesting peculiarity: the molecules on graphite do not permanently stick to each other when interacting side-to-side and do not form any aggregate until a -stacked molecule is deposited on top of them, i.e., a molecule virtually belonging to ML2. We can affirm that the presence of these top-stacking molecules is fundamental for the formation of the aggregates, allowing the underneath molecules to organize and favoring also the in-plane crystallinity (Figure 2b).
Since aggregates form only when stacked molecules are present, the growth mechanism attributed to the PFP system is the Stranski-Krastanov one, or island-plus-layer growth. Indeed,  we could initially observe the formation of several small islands on the surface, which subsequently coalesce into one major domain and a second smaller one with a different planar orientation. We identify the second domain as a crystalline defect of the growth, having also a slightly different morphology due to the obstruction of the first domain which occupies most of the surface. Overall, we find a comparable morphology to the one depicted by Salzmann et al. [37] on graphene-coated quartz, where the PFP molecules organize in a flat-lying configuration in contrast with the herringbone structure found for the bulk phase or in the PEN system.
Stacking interaction and adsorption energies can help to rationalize this different behavior. As shown in Table 1, PEN-PEN couples present face-face stacking interactions weaker than those of a PFP-PFP pair, that tend to remain stacked. Actually, it must be noticed that the adsorption energy as well is higher for PFP, thus suggesting a greater affinity for the surface compared to PEN, analogously to what was found experimentally on a similarly flat support, hexagonal boron nitride. [47] Indeed, for a more complete picture of the growth mechanism, energy barriers, that will be discussed later on, must also be taken into account.
Another relevant difference emerging from the simulated samples consists in the molecular orientation. While PEN molecules in the first monolayer always adopt a planar orientation on the surface, once the ML2 begins to form, both the underlying ML1 and all above molecules acquire a tilted configuration and arrange in the classical herringbone packing (Figure 2a). Indeed, this arrangement-as reported by Goetzen et al. [36] and found also in the bulk polymorphs of PEN and PFP [48] maximizes the quadrupolar interactions. [49,50] Considering Table 1. Intermolecular interaction energies (in kcal mol −1 ) for two parallel molecules displaced along their molecular inertia axes, and adsorption energy of PEN and PFP on graphite. Note that for a fair comparison of the two types of energies it should be considered that the former are bimolecular quantities while the latter is unimolecular.  (Table 1), which makes the detachment from the graphite surface easier. As a result, we find that the deposition simulation of pentacene mixture on graphite is in accordance with the experimental observations on graphite [36][37][38][39][40][41][42][43][44][45][46][47][48][49][50][51] and graphene. [52] In the heterogeneous system (MIX), knowing that the deposited species at each step was randomly chosen throughout the simulation, homogeneous, and heterogeneous regions, either amorphous or crystalline, could in principle turn out; however, since the final morphology is almost entirely arranged in a mixed 1:1 crystal (Figure 2c), we surmize that this phase is thermodynamically favored and that the molecular diffusive exchange is fast enough to allow its formation. In completely crystalline regions, each PEN molecule is adjacent to a total of six neighboring PFP molecules (two along each molecular axis) and vice versa. The MIX system has several similarities with the PFP one: all molecules are flat-lying, and the growth mechanism works through the early formation of crystalline islands, as we can infer from Figure 1, where ML2 appears quickly thanks to the strong -stacking component. This Stranski--Krastanov type of growth was reported also by D'Avino et al. for the deposition of the mixed crystal on HOPG. [32]

Crystalline Features
The crystallinity of the final morphologies was characterized through the calculation of the unit cell parameters, reported in Table 2 along with the parameters obtained experimentally. The good agreement between the simulated and experimental parameters proves that the employed methodology qualitatively succeeds in simulating the film growth and structure, and then can be reliably employed to make previsions in similar systems and in the calculation of further physical observables discussed in the next sections.
The new crystalline parameters were used to replicate the unit cell and calculate the cohesion energies, which we defined and calculated as the energy loss upon extraction of one molecule from a bulk crystal. The values of the on-graphite grown crystals are reported in Table 3 along with the ones for the bulk crystals (for calculation details, see the Experimental Section). Significant differences are found between the cohesive energies of the crystals grown on surface, which singularly agree with the pre- Table 2. Crystalline cell sides (Å), angles (°), and volume (Å [3] ) of PFP, PEN, and MIX systems on graphite, obtained experimentally and computationally. Experimental PEN parameters refer to the Siegrist phase, the one obtained on graphite by Goetzen et al. [36] All unit cells contain two molecules.  viously attributed growth mechanisms. The higher cohesion energy of PFP with respect to PEN implies a greater stability of the growing PFP nanocrystals and therefore should favor an island type of growth with respect to a layer-by-layer one, in line with the results for the evolution of the coverage shown in Figure 1. This driving force is even more relevant for the heterogeneous system, that, like PFP, displays a more 3D type of growth with respect to PEN. Remarkably, withdrawing a molecule from the MIX system is more costly than removing it from the homogenous crystals. The difference is conspicuous for PFP (≈7.5 kcal mol −1 ), and lower for PEN (≈2.5 kcal mol −1 ), for a total enthalpic stabilization of the cocrystal, with respect to phase separation in pure PFP and PEN crystals, of about 10 kcal mol −1 . This higher stabilization of the mixed interaction is qualitatively confirmed by the enhanced thermal stability reported experimentally for the mixed system with respect to pure ones in thick [30] and monolayer films [53] grown on other surfaces.

Energy Barriers for Molecular Diffusion
Molecular diffusion is one of the factors influencing the formation of crystalline aggregates, allowing the molecules to spread over the system and interact with a higher number of molecules. [56] Even though the diffusion can be hindered by the energy barriers present in the growing sample, experimental depositions are usually performed at temperatures and deposition rates that allow the molecules to override these barriers. [57] We studied two types of diffusive energy barriers: planar, for the diffusion on the monolayer surface, and orthogonal, responsible for the descent and the climb between the MLs. www.advancedsciencenews.com www.advtheorysimul.com

Energy Barriers for Diffusion on Terraces
It is important to stress that the process of diffusion on formed flat terraces, while considered an important factor for adatoms in inorganic crystal growth, [58] it is often overlooked in kinetic models of organic crystal growth (see, e.g., ref. [59]), probably because in the typical case of films of standing molecules, barriers are low [60] and the process is not the rate determining step ruling the growth mechanism. First, we focus on the homogeneous systems. For PEN, Figure 3 shows that the free energy profile-obtained by displacing an extra molecule over a small region on ML3 and Boltzmann-averaging over all possible rotations-is structured with minima localized on molecular cores and barriers separating them, especially along the direction of the long molecular axes (≈7 kcal mol −1 ), while smaller barriers are present in the direction of the medium molecular axes. The PFP energy map appears flatter but qualitatively very similar, with the minima located at the molecular centers and the barriers at molecular edges. Interestingly, this situation facilitates the diffusion on neighboring molecules in one direction, while the other one shows more significant barriers (≈4 kcal mol −1 ). The easy diffusion of perfluoropentacene molecules along one direction should assist the previously hypothesized "stacking effect" of a top loaded molecule for the generation of the PFP aggregates (see also Figure 2b). For the MIX system we report the energy maps relative to the separate diffusion of both PFP and PEN molecules in the same selected region. As expected, the energy minima match the alternate interactions with like and unlike molecules, and switching the species of the probe molecule determines an inverted energy map, owing to the large and opposite molecular quadrupole. We also find a steeper energy profile compared to the homogeneous surfaces, with barriers of around 8 kcal mol −1 and the presence of low energy channels between the minima. Accordingly, we can assume that the diffusion on the planar crystalline surface is generally arduous, but the presence of these channels permits the molecular exchange that yields the alternate configuration in the mixed crystal. On the other hand, the presence of deep minima in correspondence of molecules of unlike chemical species is the thermodynamic driving force for the mixed crystal formation and once the molecules reach one of those minima, they would favorably linger in it despite the presence of these channels. The strict link between surface structure and molecular diffusion exhibited by all systems hints to a general phenomenon that could occur for organic molecules on organic or inorganic crystalline surfaces, as highlighted by recent studies. [61][62][63]

Step-Edge Energy Barriers
The energy barriers that a molecule has to overcome to cross a step edge between two terraces are important parameters in the growth and nucleation of thin films from vapor, albeit difficult to measure. [64] As in most cases, and in particular for terraces of standing elongated molecules, the downward crossing barrier (the Ehrlich-Schwoebel barrier [65,66] ) is lower than 10 kcal mol −1 while that for the upward process is at least twice as large; [67] the latter is rarely discussed in the literature. However, for films of flat-lying molecules, since the step edge height is smaller, also the Table 4.
Step-edge energy barriers (kcal mol −1 ) for PFP, PEN and MIX crystal terraces, calculated for the steps with the molecular long axes parallel to the step line. For PEN, we distinguish two different barriers depending on the tilt of the molecular plane at the step, while for MIX we discriminate homogeneous (HO) and heterogeneous (HE) barriers. Capital letters correspond to initial and final positions along the diffusion paths, as depicted in reverse process might be possible at least in principle, therefore we evaluated both downward and upward barriers from the free energy maps at the steps. For recumbent PEN and PFP molecules one can further distinguish between two types of step-edges in which either the long or the medium molecular side is parallel to the step. Since we found the same qualitative trends but higher barriers for the latter case, we discuss here only the former, providing the results for the other type of step in the Supporting Information.
In addition, two situations are possible for the PEN system, since, owing to the herringbone packing, the molecular plane has two possible tilts at the step. Indeed, we observe that crossing the step is more likely to happen when molecules are leaning toward the surface ("⇄ ∕" in Table 4), and much easier downward than upward. This finding confirms the facility of the molecular descent and the consequent layer-by-layer growth found in PEN system. Instead, the PFP molecules remain for a longer time sticked on top of the others, due to their higher downward step-edge barrier, and consequently follow an island-plus-layer growth regime.
In the MIX maps the picture is more complicated: we can distinguish two types of possible displacements, through the so called homogeneous and heterogeneous channels, meaning that a molecule can cross a step constituted either by molecules of the same specie or of the other one. A schematic picture of the two channels is showed in Figure 4c,d. We established that this mainly happens through the heterogeneous channels (B→C) due to the lower descent barrier (Table 4). Instead, proceeding along the homogeneous channel implies that the probe molecule moves across molecules of the same species at the edges. In this case the starting and ending points (the minima of the potential energy surface) are not on the two molecules defining the step (of the same species), but on the adjacent molecules of the other species (points A and D). Since these spots are relatively far away from the step, they also are almost isoenergetic, and therefore the downward and upward barriers are accordingly very similar. To conclude the discussion, it is worth noting that despite the downward barriers are generally lower than the upward ones for our systems, in deposition experiments it is most often found that film roughness increases with coverage, owing to a predominance of the flux of incoming molecules on the terraces surface over the downward flux at terrace edges, which in turn is affected not only by the step-edge barrier but also by the barrier for diffusion over the terraces.

Thermal Annealing
The annealing procedure is employed experimentally to possibly modulate the crystallinity of organic semiconductors films at a judiciously chosen temperature, however not known beforehand. [6,68] Here we simulated this process at 500 K for 5 ns and cooled back the samples at 240 K for 3 further ns, to assess if the molecules can overcome the previously discussed diffusive barriers and to establish the magnitude of the consequent molecular displacement, with the objective of providing hints on the morphology alteration in relation to the analyzed growth mechanisms. Thermal diffusion mainly concerned the upper, incomplete MLs of the systems since the top molecules are less obstructed and can diffuse in all directions. Also, the upward diffusion became possible, but was overridden by the more probable molecular descent. The annealing step eventually caused a reduction of the total film thickness and roughness.
The three films varied their MLs population to a different extent (Figure 5). Although the high temperature allowed the molecules of PEN and PFP samples to override the diffusive energy barriers and distribute in the inferior MLs, the MIX system showed an inferior population exchange from the higher MLs to the lower ones. This result seems surprising since here the descent barrier, at least through the heterogeneous channel, is the lowest among the systems. We believe that the origin of this result is that the displacement across monolayers is unlikely because of the difficulty of the molecules in reaching the MLs borders, since it happens only through specific channels (Figure 4). In contrast, we mark out the tendency to eliminate the crystalline vacancies. The high cohesion energies (Table 3) explain this reorganization and the crystallite stability, which is not breached by the annealing process. A direct consequence of ML population shift upon annealing is the variation of the topography of the system, which in turn affects the electrostatic landscape at the film surface, an important quantity in organic electronics, since it affects charge motion in transistors [69] and charge separation in solar cells, [70,71] processes strongly connected to the global efficiency of those electrical devices.
To quantify the variations of both quantities, in Table 5 we report the standard deviations of film height (i.e., roughness) and electrostatic potential before and after the annealing process. While topography roughness clearly decreases upon annealing, the corresponding electrostatic roughness seems to be less affected by this high temperature treatment. Considering that the consequences of annealing are qualitatively similar for all systems, in Figure 6 we highlight the changes for the MIX system only. First, we notice the topography variations and the evidence of downward displacement of some molecules.  Moreover, the electrostatic potential looks correlated not only to the film topography, but also to the chemical nature of the upper molecules in the film. In fact, the opposite charge distribution and quadrupole of PEN and PFP determines an inverted potential for the two molecules, whose magnitude is enhanced when those molecules are located on top of terraces.

Conclusions
We simulated the vapor-phase deposition of pentacene, perfluoropentacene, and a blend of the two molecules on graphite, replicating previous experimental results and deepening the study of the key factors involved in the thin film growth. The depositions were carried out at low temperature to allow the growth mechanism to be mainly driven by intermolecular interactions rather than by kinetic factors. This approach allowed the individuation of peculiar characteristics for the growth of each of the simulated systems, which overall produce crystalline structures formed by recumbent molecules through a layer-by-layer growth for pentacene, and an island plus layer growth for perfluoropentacene and for the 1:1 mixed system. Through the evaluation of the free energy barriers for molecular diffusion, both on terraces and at terrace steps, we found a very strict correlation to the crystal morphology, that could impart anisotropy to the diffusion process. Moreover, we assessed that for these thin films of recumbent aromatic molecules, the diffusion barriers on terraces are similar to those at step edges, and that upward crossing at the edge is also possible.
More in general, the study demonstrates that the proposed methodology for simulating vapor deposition is an efficient tool to investigate at the nanoscale the processes of nucleation and growth of organic crystals, including mixed ones, thanks to the possibility of a visual and mathematical analysis of simulated molecular systems, which is just not possible or more limited with analogue experimental techniques.

Experimental Section
In all simulations, both graphite and pentacene atoms were described with the intermolecular and geometrical parameters defined in the General Amber Force Field (GAFF); [72] while graphite atoms are chargeless, the Electrostatic Potential (ESP) atomic charges of the isolated PEN and PFP molecules were calculated at density functional theory level with the PBE0 hybrid functional and the cc-pVTZ basis set, using the Gaussian16 program. [73] Molecular dynamics simulations were performed with the open-source program NAMD. [74] Periodic boundary conditions (PBCs) were applied. The Particle Mesh Ewald (PME) method [75] was used to compute the electrostatic forces with PBCs, using a cut-off of 12 Å for the calculation in the direct space, as well as for truncating Lennard-Jones interactions.
Deposition on Graphite: All the deposition simulation experiments were performed on the same template graphite surface, composed by four alternated layers. This surface was constructed with the aim of containing a high number of PFP molecules which would completely occupy the surface, in a hypothetical 100% crystalline phase. Therefore, the graphite surface was built approximately matching the surface size (a and b axes) of the ab experimental cell size of PFP on quartz-coated graphene, calculated by Salzmann et al. (a = 15.13, b = 8.94, c = 6.51). [37] The initial cell of graphite (a = 2.456, b = 4.254, c = 6.696), which contains two layers, was replicated respectively 48 and 32 times along the ab surface axes. Moreover, the surface was replicated 2 times along the c axis to obtain four layers. The resulting sample had a surface size of 117.888 × 136.128 Å. [2] Each monolayer could contain ≈118 flat-lying PFP molecules.
The NAMD package in combination with a Python script was used to simulate subsequently the vapor deposition of each molecule of PFP and PEN on the graphite surface. The atomic coordinates of the two underneath graphite layers were kept fixed in all the simulations. Volume was kept constant and atomic velocities were rescaled every 100 fs to maintain the temperature constant (240 K) as well. After each deposition interval (1 ns), a new molecule was introduced in the simulated system, at 40 Å above the highest graphite layer with a random position and orientation. For the mixed PEN/PFP deposition, the new deposited molecule is randomly chosen between the two species to better simulate the experimental conditions. PBCs were applied along x and y axes (a, b graphite axes), while the orthogonal dimension was set at 400 Å to simulate the vacuum. A total number of 450 molecules was deposited for each system. After the deposition, the unit cell parameters of the "on-graphite" polymorphs were extracted from a chosen crystalline region and averaged over 25 configurations corresponding to the last 200 ps of the simulation.
Cohesive Energy Calculation: First, the newly obtained crystalline cell of the systems, obtaining a supercell of 90 molecules was replicated. Bulk polymorphs initial coordinates and cell parameters were downloaded from the Cambridge Structural Database. After a conjugated gradient minimization at fixed volume and cell dimensions, the cohesive energy was calculated as a sum of internal energies as follows Where E TOT is the energy of the supercell, E TOT−MOL is the energy of the supercell without one molecule and E MOL is the energy of the withdrawn molecule inside the otherwise empty supercell. Free Energy Calculations: The adaptive biasing force (ABF) method [76] was applied for the evaluation of the adsorption energy and the potential energy of the crystal surfaces, along suitable collective variables (colvars). These calculations were performed at the same temperature, volume, and pressure of the deposition ones.
1) Absorption Energy Barrier: the ABF method was applied on a simple system containing the graphite template and one molecule of either PFP or PEN. The chosen colvar corresponded to the orthogonal distance from the surface, and the absorber molecule was free to explore the space along that direction. To delimit the colvar space, dummy potentials barriers were set at the borders of the colvar region, extending until 30 Å from the highest layer of graphite. The adsorption free energy was therefore calculated as the difference between the minimum value of the free energy and the value at 30 Å. 2) Mapping of the Potential Energy Surfaces: one molecule was let free to move on a crystalline portion of each sample delimited by two colvars corresponding to the positions of the molecular center of mass along x and y, leaving unconstrained the z coordinate unconstrained. All other molecules were maintained at their original position by fixing one atom per molecule. For on-terrace energy barriers, the colvar space comprehended a 3×3 molecules region in PEN and PFP, and a 4×3 region in MIX, while for step-edge barriers 4×2 regions at the interface between two MLs where chosen. To calculate the step-edge energies reported in Table 4, the maps of Figure 5 were averaged along the y axis (i.e., orthogonally to the climb/descent direction), obtaining the energy profiles from which the barriers were extracted. For MIX, the energy maps in two regions were separated and obtained both the homogeneous and heterogeneous averaged profiles.
Annealing: The final obtained morphology underwent a heating process with a sudden increase of temperature to 500 K. The system was kept at 500 K for 5 ns, before it was cooled down at 240 K again and equilibrated for 3 ns. Properties were calculated for the last configuration of this step.
Electrostatic Potential Maps: For each of the three samples, we calculated the potential energy of the system by inserting a test point charge of +1 e in the system with PBC. In each calculation, the charge was positioned on a different point of a 3d grid (16 166 points) corresponding to the calculated topography of the final configuration, which was set 3 Å above the center of mass of the exposed atoms. The potential energy of the system with no charge was subtracted to the calculated energy in each point, resulting in the electrostatic potential map.

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.