Structure Prediction for Surface-Induced Phases of Organic Monolayers: Overcoming the Combinatorial Bottleneck.

Structure determination and prediction pose a major challenge to computational material science, demanding efficient global structure search techniques tailored to identify promising and relevant candidates. A major bottleneck is the fact that due to the many combinatorial possibilities, there are too many possible geometries to be sampled exhaustively. Here, an innovative computational approach to overcome this problem is presented that explores the potential energy landscape of commensurate organic/inorganic interfaces where the orientation and conformation of the molecules in the tightly packed layer is close to a favorable geometry adopted by isolated molecules on the surface. It is specifically designed to sample the energetically lowest lying structures, including the thermodynamic minimum, in order to survey the particularly rich and intricate polymorphism in such systems. The approach combines a systematic discretization of the configuration space, which leads to a huge reduction of the combinatorial possibilities with an efficient exploration of the potential energy surface inspired by the Basin-Hopping method. Interfacing the algorithm with first-principles calculations, the power and efficiency of this approach is demonstrated for the example of the organic molecule TCNE (tetracyanoethylene) on Au(111). For the pristine metal surface, the global minimum structure is found to be at variance with the geometry found by scanning tunneling microscopy. Rather, our results suggest the presence of surface adatoms or vacancies that are not imaged in the experiment.

The specific structure of any material constitutes the key to its functionality. Besides the chemical composition of a material, the way its individual constituents arrange, that is, the polymorph it assumes, strongly influences the material’s thermal,[1] mechanical, optical,[2] and electronic[3] properties. Structure determination and prediction is, therefore, the very fundament of material science. Determining polymorphs is particularly challenging for systems with more than one component, such as organic/inorganic interfaces that are prevalent in many applications, ranging from catalysis to organic electronics. Because of the interplay between intermolecular and molecule–substrate interactions, organic molecules on inorganic substrates are prone to form surface-induced phases, giving rise to a particularly rich and intricate polymorphism.[4,5] Such surface-induced phases often contain multiple molecules per unit cell and can display properties that are vastly superior to those of the bulk phase.[6] However, it is commonly a priori not assessable which polymorph the material will assume. Hence, predicting the structure of a material, ideally even before it is synthesized, is of crucial importance to pave the way toward computational materials design.

From a computational point of view, structure prediction can be considered as a global optimization problem, that is, the problem of finding the global minimum of the energy landscape. The principle difficulty in treating nontrivial global optimization problems arises from the huge number of possible minima on the multidimensional potential energy surface (PES), which increases exponentially with the size of the system.[7] In practice, it is furthermore complicated by the fact that polymorphs can be kinetically trapped, that is, also other minima besides the global minimum may play a decisive role for the structure formation at interfaces. Exhaustively sampling the corresponding vast configuration space demands an unfeasible amount of computational resources. Still a variety of different techniques, ranging from Monte Carlo or molecular dynamics-based techniques such as simulated annealing,[8,9] Basin-Hopping,[10−12] or minima hopping[13,14] to evolutionary approaches such as genetic algorithms[15−19] have been successfully developed. Although most of these efforts were dedicated to gas-phase structures or bulk crystals, recently, structure search was extended to organic/inorganic interfaces.[12,20−23] However, these methods often rely on elaborate data fitting or force fields to describe intermolecular interactions.

In this article, we present a powerful and efficient computational structure search algorithm that allows to employ fully converged first-principles calculations throughout. It is specifically designed to explore the potential energy landscape of commensurate organic/inorganic interfaces with large unit cells containing multiple molecules and to efficiently and accurately locate low-energy polymorphs. This allows us to predict global minima (in cases where the experimental growth process is thermodynamically controlled), and furthermore also provides a set of relevant, low-energy structures that are useful to verify the interpretation of experiments. In short, the approach, which we call SAMPLE (surface adsorbate polymorph prediction with little effort) combines a systematic discretization of the PES with an efficient exploration inspired by a Basin-Hopping algorithm. Interfacing SAMPLE with dispersion-corrected density functional theory (DFT), we demonstrate the power and efficiency of our algorithm by application to the strong electron-accepting molecule tetracyanoethylene (TCNE) adsorbed on the Au(111) surface. TCNE is particularly interesting for computational structure search studies, as it forms very different surface-induced phases on various metallic substrates[24−26] with structures that are markedly different to those observed in the bulk.[27,28] Employing the SAMPLE approach to TCNE/Au(111) we find that a “naïve” evaluation of the structure on the basis of the scanning tunneling microscopy (STM) image does not yield the global minimum geometry, and we provide an alternative interpretation that cannot be inferred from experimental data alone.

The SAMPLE Approach

To present our method, we will first provide a general overview of the SAMPLE procedure and introduce the core concepts. Then, the individual steps of the procedure are discussed in detail on the specific example of TCNE adsorbed on Au(111). Because there is “no free lunch”, specialization is key for highly efficient structure search methods. In particular, SAMPLE is designed for commensurate structures of mostly rigid molecules where the adsorption energy dominates over the intermolecular interactions. This has several implications, which we exploit to simplify the structure search problem. We discuss these implications throughout the main text. Furthermore, a discussion on when the related assumption might become invalid is provided in the Supporting Information.

As illustrated in Figure , the SAMPLE approach divides the structure search problem at organic/inorganic interfaces into two main parts: first, a systematic discretization of the configuration space; and second, an efficient exploration of the PES. The outcome is a set of energetically lowest lying polymorphs.

Figure 1

SAMPLE approach: A systematic discretization of the configuration space consisting of an extensive set of supramolecular configurations (left) is followed by an exploration of the PES, where efficiently chosen configurations are relaxed to the corresponding polymorphs, that is, the nearest local minima on the PES (right). The outcome is a set of energetically lowest lying polymorphs. Specifically, Step 1 shows three out of the nine stable local adsorption structures for TCNE/Au(111), that is, a flat-lying TCNE molecule, an upright-standing TCNE molecule with the central C=C bond parallel to the surface, and equivalently perpendicular to the surface. Step 2 depicts two possible configurations generated by the assembly process. In Step 3, exemplary trial moves are illustrated, that is, translation by one primitive surface unit cell (which is indicated by the dashed box), rotation, and exchange by another local adsorption geometry. Step 4 shows a schematic PES (black line) and the corresponding transformation into a set of interpenetrating staircases (red dashed line) by performing local geometry optimizations. Because of the unique labeling of all configurations, a history list containing all visited polymorphs can be provided. For details see main text.

As visualized in Figure , the discretization is performed in two consecutive steps. First, the adsorption geometries that isolated molecules on the surface would assume have to be determined. This allows us to establish which adsorption sites are adopted and in a subsequent assembly process use the substrate as discrete registry to combine these local adsorption geometries into supramolecular configurations. Here, the term “configuration” denotes a specific arrangement of several molecules in the unit cell that serve as a starting point from which later the corresponding polymorph, that is, the closest local minimum on the PES, will be determined via a local geometry optimization. These steps are fully deterministic and performed only once for a given interface. In notable contrast to commonly applied global structure search techniques, thereby an enclosed configuration space with a well-defined and reproducible number of possible configurations is generated a priori.

In principle, one could perform a local geometry optimization for each of these starting configurations to their nearest local minimum and thereby obtain a complete energy ranking. However, the number of polymorphs increases exponentially with the number of molecules in the unit cell. Hence, in many cases, even after the aforementioned discretization procedure, the configuration space is too large to perform an exhaustive search within reasonable time. This demands, as the second part of the SAMPLE procedure, an efficient exploration of the PES.

To sample the PES we explore the configuration space generated in the first part in consecutive Monte Carlo steps by iteratively suggesting new configurations followed by a local geometry optimization. The discretization of the configuration space allows us to define connections between the configurations with the main advantage of performing a selective search along the PES. As illustrated as Step 3 in Figure , we define specific trial moves for the molecules on the surface (that is, translation, rotation, and exchange by a different local adsorption geometry) and thereby set up a neighbor list for each configuration from which the next one is chosen randomly. Each suggested configuration is followed by a local structure relaxation into the closest local minimum. These minima (which can be uniquely identified by the molecular configurations and their positions) and their energies are stored for later reference and analysis. The outcome of the SAMPLE procedure is a set of energetically lowest lying polymorphs.

We note that the SAMPLE approach can be linked to any suitable electronic structure method that provides an accurate energy ranking. For all calculations in this article, we optimize the geometry using PBE+vdWsurf[29] and report adsorption energies using the (supposedly) more accurate many-body dispersions correction scheme.[30] For full details on the methodology, see Computational Methods and the Supporting Information.

Application to TCNE/Au(111)

To convey a more detailed explanation of the SAMPLE approach and to benchmark its efficiency, in the following we will apply it to the specific example of TCNE (see Figure a) adsorbed on Au(111). As mentioned in the introduction, TCNE forms very different surface-induced phases on various metallic substrates.[24−26] In the specific case of Au(111), STM experiments performed at low temperature (T = 7 K, see Supporting Information, Methods) reveal a triangular pattern in a nonorthogonal unit cell containing three TCNE molecules, as shown in Figure b. (Note that these structures have already previously been reported in ref (24).) Interestingly, in contrast to many other conjugated molecules as well as to the adsorption of TCNE on, for example, low-index Ag surfaces[24−26] or Cu(100),[24,31] the molecules do not appear in the expected flat-lying fashion. Rather, from the STM experiment they appear to be tilted onto their sides. This peculiarity makes TCNE/Au(111) an exciting candidate for employing our SAMPLE approach.

Figure 2

(a) Chemical structure of tetracyanoethylene (TCNE, C6N4). (b) Experimental STM topography of TCNE adsorbed on Au(111). The image was taken in constant-current mode (Vs = 0.1 V, I = 5 pA, T = 7 K) and the film was grown at room temperature. TCNE arranges in ordered triangular structures with a unit cell containing three molecules that are tilted onto their sides.

Step 1: Evaluating the Local Adsorption Geometries

All band-structure calculations require a unique set of lattice vectors as input. In other words, any computational structure search for interfaces is limited to structures that are commensurable. As such, each molecule can be assigned a specific adsorption site on the substrate. For SAMPLE, we assume that the adsorption site is mostly independent of the molecular coverage, that is, that the geometry of a molecule with respect to its position on the substrate is only slightly perturbed by the presence of other molecules on the surface. This means, for example, that an individual TCNE molecule might adopt a bridge or an on-top position (besides others). These are then likely also local minima (for each molecule) at high coverage (although the exact position may change, but this will be captured during geometry optimization) and thus suitable starting points for setting up the configuration space. Hence, we systematically discretize the configuration space in two consecutive steps. First, we only consider the metal–molecule interactions of single molecules, while intermolecular interactions are accounted for in a consecutive assembly process including several molecules per unit cell. This procedure significantly reduces the conformational complexity. We note that such a “divide and conquer” approach has also very recently been attempted for interfaces.[32] However, in contrast to ref (23) in our work, the intermolecular interactions are not fitted from gas-phase calculations. Rather, they are directly obtained from first-principles calculations of the molecules on the surface. Indeed, we emphasize that for the present system, which undergoes metal-to-molecule charge transfer, the intermolecular interaction even qualitatively changes between the gas phase and the adsorbed molecules.

Hence, in Step 1 we neglect the impact of intermolecular interactions and determine the stable geometries that an isolated, single molecule on the surface can adopt. In the present case, this is done by performing multiple local geometry optimizations from systematically chosen initial guesses with different molecular orientations and positions relative to the substrate (see Supporting Information, Methods for details). Altogether, we find nine possible local adsorption geometries for TCNE/Au(111), three of which are exemplarily shown in Step 1 of Figure . They differ in their particular adsorption site with respect to the substrate as well as in their orientation, that is, they are flat-lying or upright-standing. Because of the hexagonal lattice of the gold surface and the mirror symmetries to TCNE, each of the nine stable local adsorption geometries has three symmetry-equivalent isomers, which need to be considered in the second assembly step. A comprehensive list of all structures together with their adsorption energies can be found in Figures S1 and S2 in the Supporting Information.

Step 2: Assembling the Configurations

One of the major challenges in any structure search algorithm is the exponential explosion of the configuration space, that is, the exponentially increasing number of configurations with system size. Exemplarily, for a unit cell containing three molecules, considering five degrees of freedom for each molecule on the surface and four different values per degree of freedom, an unfeasible number of 4(3·5) > 1 billion configurations would be generated. Applying our discretization procedure, this problem is greatly mitigated. After evaluating the local adsorption geometries in Step 1, we assemble them into supramolecular configurations, as exemplarily illustrated in Step 2 of Figure . In other words, the systematic preoptimization allows us to “freeze out” the internal coordinates during the assembly step, where then only intermolecular interactions have to be considered.

For this assembly procedure, we use the primitive surface cell of the substrate, that is, the p(1 × 1) (equivalent to √3 x √3 R30) unit cell of Au(111), as registry. In principle, we try all combinations of all possible local adsorption geometries in all rotations on each surface unit cell and discard all configurations that do not have the desired coverage or that are symmetry-equivalent to another configuration. Additionally, we remove configurations that are unphysical, because they interpenetrate or come too close to each other. Specifically, we exclude all configurations where the minimum distance between two atoms of adjacent molecules is smaller than a predefined threshold. Here we use 2.4 Å, but we note that for the present system, the number of configurations is not overly sensitive to the exact choice of this parameter, as shown in Figure S3 in the Supporting Information. For details on our technical implementation, which avoids the configurational explosion associated with this assembly procedure, see Figure S4.

The extent of the generated configuration space, that is, the number of possible configurations, depends on the number of local adsorption geometries that are determined in Step 1 as well as the molecular packing density, that is, the size and shape of the unit cell together with the number of molecules it contains. To obtain the global minimum, structures with different unit cell sizes, shapes, and number of molecules can and should be systematically generated. At this point, it is also possible to include experimental information, if available. For example, the experimental size of the unit cell might be known from various techniques, such as diffraction experiments. For TCNE on Au(111), the structure contains three molecules in a ∼233 Å2 cell (31 surface Au atoms). Because this structure must be commensurate with the gold lattice, there are only five distinct unit cell possibilities, see Figure S5 in the Supporting Information. From these options, the experimental STM shown in Figure b is only consistent with the unit cell (prompting us to neglect the others, because their calculations would only cost computational time but not provide more insight into the physics or the performance of our method). With this input, the assembly procedure generates approximately 200 000 different configurations.

The discretization of the configuration space allows us to introduce a unique labeling for each of these configurations by storing the local adsorption site and position for each molecule in the unit cell using the substrate as registry. This leads to a major advantage in the exploration step, as discussed below.

Step 3: Setting up the Connections between Configurations

To explore the PES in a stepwise Monte Carlo procedure we sample the configuration space generated in the first part in consecutive steps by iteratively suggesting new configurations followed by a local geometry optimization. This is in principle similar to traditional Basin-Hopping, where the complex PES is transformed into a set of staircases by consecutive hops followed by geometry relaxations into the closest local minimum. Hopping between the configurations is performed by randomly choosing trial moves along certain trajectories, typically Cartesian or internal coordinates.[11,12] In contrast to this, in SAMPLE a more efficient approach is pursued by exploiting the fact that each configuration in the discretized configuration space can be systematically connected to a certain number of neighboring configurations by performing selective trial moves. Specifically, the neighborhood for a certain configuration is generated according to the selection rules illustrated in Step 3 of Figure : we allow each molecule in the cell to (i) move to an adjacent space of the substrate, that is, translation by one primitive p(1 × 1) surface cell (that is, shift to the next possible equivalent adsorption site on the substrate), (ii) rotate on the spot to a symmetry-equivalent structure, or (iii) adopt a different adsorption site, that is, exchange to a different local adsorption geometry. Thereby, each configuration obtains an individual neighbor list defining all connections within the configuration space.

Step 4: Exploring the Configuration Space

The iterative procedure to explore the PES is illustrated in the flowchart depicted in Figure . Starting from a certain configuration we suggest a random neighbor according to the neighbor list defined in Step 3. Depending on whether this neighboring configuration has already been visited or not, we would either just look up its energy in a history list or perform a local geometry optimization to the nearest local minimum on the PES. Note that this relaxation is comparably inexpensive as we start from a combination of already optimized local adsorption geometries (convergence is reached after 3–4 geometry steps on average). The suggested polymorph is accepted or rejected on the basis of the Metropolis-Hastings scheme[32]If the new energy is lower than the one of the last accepted polymorphs, it is accepted; in the case of a higher energy, it will be accepted depending on a probability that considers the energy difference to the last accepted polymorph, ΔE, as well as an effective temperature via a Boltzmann factor, β. In the case of rejection, a different neighbor is chosen. The initial temperature was set to 300 K and was varied upon iterations, that is, the temperature was decreased (increased) by 100 K in the case of acceptation (rejection). More details are given in the Supporting Information, Methods. Because of the unique labeling of each polymorph, all visited structures (accepted or rejected) can be added to a history list and need not to be recalculated in case of revisiting. As in traditional Basin-Hopping, this exploration run (as depicted in Figure ) is repeated several times for different starting configurations in order to ensure an unbiased structure search on the PES. However, in SAMPLE the history list of visited structures can be easily transferred from one run to the other, and as we will show in the following this enormously increases the overall performance.

Figure 3

Flowchart representing the iterative Monte Carlo procedure to explore the PES. The iterative process is stopped after a predetermined number of geometry optimizations with DFT. For details, see the main text.

Benchmark, Validation, and Performance

In order to validate the performance of the exploration aspect of the SAMPLE approach, it is useful to reduce the conformation space to a size where a comprehensive search is possible, that is, to a system where we can calculate from first-principles, all configurations that are created in the assembly procedure of the SAMPLE approach and see how efficiently the global minimum is found. Although SAMPLE does not necessarily rely on experimental input at all, it straightforwardly allows one to incorporate information if available. We note that in principle such structural information could also be obtained from any experimental method, such as low energy electron diffraction techniques. For the case of TCNE/Au(111), it suggests itself to incorporate the information that TCNE forms triangular structures (cf. Figure b).

Discarding all conformations where this is not the case (that is, automatically rejecting all trial moves that lead to nontriangular structures in Step 3 of Figure , as described in detail in Figure S6 in the Supporting Information) reduces the complexity from about 200 000 to 144 configurations. The corresponding Sub-PES, that is, the total energy of each of these polymorphs after optimization to their local minimum is shown in Figure a. Each box of this PES represents one of the 144 triangular polymorphs. The global minimum is indicated by the red box.

Figure 4

(a) Sub-PES for TCNE/Au(111) comprising 144 triangular polymorphs. Each of the 144 boxes refers to one of the possible polymorphs generated by the assembly procedure. The color represents the total energy obtained after geometry optimization. Energies are given relative to the global minimum. The global minimum is indicated by a red box, while the polymorph with rank 30 (that is, 0.35 eV higher in energy than the global minimum) is highlighted by the white box. The arrangement of the boxes is described in detail in Figures S8 and S9. (b) Probability to find the global minimum on the PES shown in (a) for an increasing number of DFT evaluations. For details see main text.

On basis of this Sub-PES we evaluate the efficiency of our SAMPLE approach by estimating the probability to find the global minimum within a certain number of first-principles evaluations, that is, the number of geometry optimizations, independent of how many SAMPLE cycles where the energy is just looked up (which is essentially free) have been performed. To ensure an unbiased structure search we start the exploration (as described in Figure ) from a completely random starting configuration on the PES and stop it after a finite number (Nmax) of geometry optimizations. As described above, the peculiarity of our approach is the fact that we explore the PES within a finite configuration space consisting of unique configurations that can be labeled and stored for revisiting. Hence, we provide the energies of already calculated polymorphs in the form of a history list, thereby avoiding expensive recalculation of already visited parts of the PES. Because of the stochastic nature of our approach the whole procedure is repeated 10 000 times to obtain accurate statistics, each time choosing a new random seed for the Metropolis-Hastings evaluation as well as for the starting configuration. Figure b shows the probability to find the global minimum on the Sub-PES as a function of the performed geometry optimizations. When applying the SAMPLE procedure (blue line), the global minimum is found with a probability of approximately 85% already after 40 DFT evaluations (which corresponds to visiting 28% of the PES). In the present case, after 70 DFT evaluations (about 50% of the PES) the global minimum is practically guaranteed to be found. If we do not have a history list (which would be the case in traditional Basin-Hopping that cannot label configurations), the performance is considerably decreased. As shown in the red line in Figure b, the probability to find the global minimum after 40 DFT evaluations only amounts to 50%. Furthermore, for a completely random search, that is, if we did not employ system-specific selection rules and did not have a history list either, the probability would decrease to 25% after 40 DFT evaluations (see green line in Figure b).

Comparison to Experiment

While the evaluation clearly shows that our approach is capable of efficiently determining the global minimum structure on the PES, it is also substantial to validate the results by comparison to the experimentally obtained structure. Applying SAMPLE to TCNE/Au(111), we find a global minimum that is indicated by the red box in the PES of Figure a with the geometric structure shown in Figure a. Although as isolated species, the flat-lying TCNE molecule is found to be the most stable structure indeed among the triangular structures, a configuration consisting of upright-standing molecules is obtained. The agreement between the global minimum and the experimental minimum looks in principle very reasonable.

Figure 5

Geometric structure and overlay to the experimental STM of (a) the predicted global minimum on the in Figure a (red box) shown sub-PES consisting of triangular polymorphs of TCNE/Au(111), (b) the polymorph obtained from a “naïve” starting point where the molecules where placed in the middle of the spots in the STM image, corresponding to rank 30 in our PES evaluation, (c) the predicted global minimum including a gold adatom (pink) in the center of the triangular structures, and (d) the predicted global minimum with a vacancy in the center of the triangular structures.

Yet, interestingly, compared to the experimental picture the TCNE molecules of the global minimum appear to be too close to each other. Indeed, if we start a geometry optimization the same way it would be done without a structure search algorithm, that is, based on the experimental picture we place all TCNE molecules in the middle of the spots in the STM image (Figure b), we arrive at a polymorph that is not the global minimum (shown in Figure b). Rather, this structure is 0.35 eV higher in energy than the global minimum. Importantly, the same polymorph was also found with our SAMPLE approach. It is located at rank 30 in the overall hierarchy of all calculated polymorphs highlighted in Figure a by the white border. Although the fact that SAMPLE finds the same structure as a (potentially) better start from the experimental STM corroborates our approach, the question remains what underlies the apparent discrepancy between our structure search efforts and the experimentally observed result.

It is, in principle, conceivable that the reason for this deviation is rooted in the inability of the employed electronic structure theory method (here PBE+MBD) to correctly reproduce the real potential energy surface. Indeed, it is well-known that the PBE functional generally misjudges the amount of charge transfer,[33−35] and tends to overdelocalize charge at interfaces.[36,37] Furthermore, we neglected the energy contribution from vibrations, which may change the ordering of the free energy relative to the DFT energies. On the other hand, the same method that we employ here has been recently used to gauge the hierarchy of organic bulk crystals in a blind test[38] with outstanding success and also reproduces desorption energies of organic molecules on coinage metals very well.[39] Furthermore, we re-evaluated the adsorption energies of our local adsorption geometries using the meta-GGA functional SCAN[40] and found only minor changes (see Supporting Information, Figure S7). We therefore expect that it is unlikely that the polymorph at rank 30 would actually be the true global minimum if we had used an even more accurate method.

Another possible interpretation of our results is that the experimentally observed structure is not in thermodynamic equilibrium. According to Ostwald’s rule of stages,[41] many crystals go through higher-energy polymorphs before assuming the thermodynamic equilibrium arrangement. Indeed, the preparation of the sample was performed at room temperature, while the measurements were done at 7 K, which could potentially “freeze” the prevalent conformation at that stage. On the other hand, most polymorphs of organic molecules that are experimentally observed exhibit energy differences of less than 0.075 eV,[42] much less than the difference between our global minimum and the rank 30 polymorph. Also, the cool-down procedure prior to insertion of the sample into the liquid helium cooled STM involved a relatively slow (about 40 min) precooling from 300 to 80 K using liquid nitrogen, which makes kinetic trapping unlikely.

A possibly more plausible scenario is that an implicit assumption in the interpretation of the STM data, namely that the system solely consists of TCNE molecules on a pristine Au(111) surface, is misleading. Especially the fact that the molecules seem to appear too close with respect to each other in the calculations (compared to experiment) suggests that the system may contain aspects that are not easily imaged with STM. These could be vacancies[43,44] or surface adatoms[26,45−47] both of which have been observed in several experimental studies. Therefore, we chose also to explore these possibilities and apply SAMPLE to predict the global minimum for these two scenarios. To this aim, we repeated the discretization and exploration procedure for triangular structures including a central adatom or a central vacancy. The corresponding Sub-PESs are shown in Figures S9–S11. The polymorphs of the most stable structures are depicted in Figure c and d. Indeed, the presence of either an adatom or a vacancy increases the relative distances between the molecules to roughly the experimentally observed positions. We note that because the three scenarios (pristine surface, adatom, and vacancy) contain a different number of atoms, their energies are not directly comparable in our calculations (unless an assumption about the chemical potential of Au is made, which is not viable here). However, in particular for the polymorph including the adatom an excellent agreement to the experimental STM is observed, which makes us confident that this adatom structure reflects the actual situation. This situation nicely illustrates that computational structure search is a powerful tool to verify and augment the interpretation of experiments and to pinpoint “hidden” aspects that might not be covered by a particular experimental technique. The concluding interpretation, of course, still lies in the hand of experimentalists, which can confirm–or refute–the options presented by theory.

In conclusion, we have introduced SAMPLE, an innovative computational structure search algorithm to explore the potential energy landscape of organic/inorganic interfaces upon first-principles. Our approach combines a systematic discretization of the PES with an efficient exploration realized as a Monte Carlo method inspired by the Basin-Hopping algorithm. In stark contrast to commonly applied global structure search techniques, in a first assembly process we a priori generate a complete and enclosed configuration space with a well-defined and reproducible number of possible configurations. In a subsequent exploration step, an efficiently chosen set of energetically low-lying polymorphs including the thermodynamic minimum is predicted. Our approach is inspired by traditional Basin-Hopping but is generally more efficient. The discretization is designed such that the resulting configuration is already very close to a local minimum, ensuring that the consecutive local geometry optimization is highly efficient. Furthermore, our method allows to efficiently revisit and cross already known parts of the configuration space without any additional computational effort, as the first assembly process enables the unique labeling of each polymorph. The efficiency is particularly enhanced by applying systematic trial moves according to carefully chosen transition rules, instead of the commonly applied random trial moves. We benchmarked our approach on the example of TCNE/Au(111) by interfacing SAMPLE with dispersion-corrected DFT and validated the predicted global minimum by comparison to the experimental STM data. Besides demonstrating the efficiency and power of the SAMPLE approach, we could moreover provide an alternative interpretation that cannot be inferred from pure experimental data. These aspects endow the SAMPLE approach with significant potential for refined structure determination and prediction of other interface materials far beyond the system discussed here.

Computational Methods

All electronic structure calculations were performed within the framework of DFT using the Fritz-Haber Institute ab initio molecular simulations (FHI-aims) package.[48] The Perdew–Burke–Ernzerhof (PBE)[49] exchange-correlation functional was applied and dispersion correction was included using vdWsurf[29] for geometry optimizations and many body dispersion (MBD)[30] for subsequent single point calculations. When changing from vdWsurf to MBD, we observe qualitative differences in the energetic ordering of the local adsorption geometries as well as the polymorphs spanning the PES (see Figures S2, S8, and S10 for more details). To obtain the local adsorption geometries, we used a (6 × 6) unit cell to avoid spurious interaction between neighboring cells. For all supramolecular polymorph calculations the experimental unit cell was measured from the representative STM image shown in Figure b with an epitaxy matrix of . We employed a tier 1 basis for Au (excluding g and h basis functions) and tier 2 for N and C basis functions. A converged 4 × 4 × 1 k-point grid was applied for all calculations. Geometry optimizations were performed using the trust radius method enhanced version of the Broyden-Fletcher-Goldfarb-Shanno optimization algorithm until the remaining forces were less than 0.1 eV/Å. For the local adsorption structures as well as for the polymorphs we used four layers of gold and relaxed the uppermost two layers together with the monolayer. VMD,[50] VESTA,[51] and Matlab were used for graphical visualization, and python was used to implement SAMPLE. For full details on the applied computational methodology and numerical parameters used in our calculations, see Supporting Information, Methods.

Experimental Section

The experiments were performed in ultrahigh vacuum (UHV) using a home-built STM operated at T = 7 K. TCNE crystals (99% purity) were kept in a small vacuum container that was cleaned by repeated cycles of flushing with Ar gas and pumping. Its high vapor pressure of ∼2 × 10–3 mbar at room temperature[52] permits controlled dosing of TCNE gas from the container into the UHV system through a leak valve. The purity of the TCNE gas was checked with quadrupole mass spectrometry. Prior to TCNE adsorption, the Au(111) single-crystal substrate was cleaned by repeated cycles of Ar sputtering and annealing. Afterward, TCNE was deposited through the leak valve onto the Au(111) substrate that was held at room temperature, leading to submonolayer amounts in fcc regions of the Au(111) herringbone reconstruction.[24] Then the sample was slowly (within about 40 min) precooled to about 80 K before the final transfer into the cryogenic STM. The STM tip was made of PtIr, and topography images were taken in constant-current mode.