Quantifying the “Subtle Interplay” between Intermolecular and Molecule–Substrate Interactions in Molecular Assembly on Surfaces

The effect of the concurrent action of intermolecular and molecule–substrate interactions on the two-dimensional (2D) self-assembly of organic molecules on solid surfaces is investigated in a combined experimental and theoretical effort. Scanning tunneling microscopy measurements of terephthalic acid on the Cu(111) surface, a model system where the interplay between the two interactions is particularly evident, are used to develop a general, simple, and computationally inexpensive model that quantitatively accounts for the experimental observations. The model, related to the well-known Frenkel–Kontorova model, offers a comprehensive description of the “subtle interplay” between intermolecular and molecule–substrate interactions and provides a qualitative and quantitative predictive capability in the design and fabrication of 2D molecular nanostructures at surfaces.

Molecules were modelled as particles in a 1-D periodic potential (representing the substrate), with periodic boundary conditions. The interactions present in this system are: (i) intermolecular interactions, described with a pairwise Morse-like potential (representing the intermolecular dimeric hydrogen bond) and (ii) the molecules' interactions with the periodic substrate adsorption potential. The model is one-dimensional, i.e. we considered only a 1-D chain of "molecules" and ignored any side interactions with molecules in neighboring chains.
The interaction of molecule i (located at coordinate x i ) with the substrate can be described by the cosine potential (x i ) with periodicity B and amplitude : The intermolecular interaction between molecule i and a neighbor molecule j can be described by the Morse-like potential: where x = |x ix j | is the distance between molecules i and j, A is the ideal gas-phase (equilibrium) intermolecular distance,  is the intermolecular interaction energy, and a describes the width of the potential energy well.
The amplitude of the cosine potential  was set as the unit of energy and the periodicity of the surface adsorption potential B was set as the unit of length. Therefore, the depth of the Morse potential is described by the relative energy  , and the intermolecular distances are expressed in terms of B  in particular, the ideal gas-phase intermolecular distance is given by ; molecular coordinates are also rescaled in units of B.
Moreover, the molecule-substrate cosine potential was shifted downwards by , thus eliminating the second term in equation (1). As a result, the intermolecular and moleculesubstrate interactions become: g(X i ) = cos(2X i ) Thus, for a system of N molecules with coordinates X i , the overall potential energy U is: S3 cos (5) and the force F i acting on molecule i is: The system was set to be periodic with periodicity P, corresponding to the number of surface adsorption sites considered.
The systematic procedure for finding the absolute energy minimum arrangement can be described as follows.
1. The number of molecules N is chosen (starting with N = 1).
2. The number of surface adsorption sites P is chosen. P is varied between some minimum and maximum values, P min and P max , which are chosen so as to accommodate the number of molecules N and therefore depend on L.
3. For each P between P min and P max , two starting arrangements of molecules are constructed: i. N molecules initially equally-spaced at distances equal to P/N (i.e. molecular positions controlled by the substrate periodicity); ii. N-1 molecules with initial intermolecular separations equal to the ideal gas phase distance L (i.e. molecular positions controlled by the intermolecular interaction) and the remaining molecule placed to accommodate the mismatch between (N-1)L and P.
4. The equilibrium structure according to the potential energy (5) for N molecules in P lattice sites is found by steepest descent minimization starting from both arrangements (i) and (ii). Typically, N was varied from 1 to 10, and by N = 10 the minimum energy per molecule was found to have converged to an approximately constant value; only in few cases was there the need to increase N beyond 10 to achieve convergence.
This procedure gives the lowest-energy commensurate superstructure (N molecules in a periodic landscape of P surface adsorption sites).
The fact that the energy per molecule became approximately constant for large values of N suggests that, after a certain threshold, there is no one preferred value of N; instead, several superstructures with different N (and correspondingly different P) are likely to coexist.
The described procedure was performed for a wide range of L and D values, to describe an ample set of structural and energetic mismatches between the intermolecular and molecule-substrate interactions. The value of a (the curvature of the Morse potential) was set to be 2.0 so as to reproduce the curvature of the intermolecular energy vs distance curve (the function plotted in Figure 3 in the main text) and was not further varied.

Molecular mechanics calculations
To obtain the ideal (gas-phase) TPA-TPA distance, the lattice parameter a along the 1- were varied (specifically, the components b x and b y of the vector b were varied, also in steps of 0.05 Å), while the parameter a was kept fixed at its 1-D value. Finally, the values of a, b x and b y were varied simultaneously around their minimum-energy values, to check that the 1-D value of the parameter a was unchanged in the 2-D system.

Density-functional theory calculations
Density-functional theory calculations of bulk Cu and Au metals, metal surfaces and surface-adsorbate systems were done using the GGA Perdew-Burke-Ernzerhof (PBE) functional, [S1] as outlined in the Computational Methods section of the main text. Test calculations of Cu and Au bulk showed that the k-grid cutoff of 15 Å (9×9×9 k-point grid) and mesh cutoff of 150 Ry gave converged values of Cu and Au lattice parameters and bulk cohesive energies. The optimized lattice constant for bulk Cu was 3.68 Å (cf. experimental value of 3.615 Å [S2] ), and for bulk Au 4.17 Å (cf. 4.079 Å in experiment [S3] ). (111)-oriented slabs were converged already at the two-layer thickness, with the surface energy of 1.62 J m -2 (Cu) and 0.80 J m -2 (Au), in very good agreement with earlier calculations. [S4, S5] Nevertheless, a thicker 3-layer metal slabs with a fixed bottom layer was used for the presented calculations.
To compare the energies of adsorption of TPA, benzoic acid and benzene on Cu (111) and Au (111), high-symmetry adsorption sites were considered first: the molecules were placed with the aromatic ring above the hexagonal close-packed (hcp) hollow site on Cu (111) and above the face-centered cubic (fcc) hollow site on Au(111). Then, to model the potential energy surfaces (PES) of TPA adsorption, the molecule was moved laterally in steps of 0.25 Å (or 0.5 Å on Au(111)) to cover the whole surface unit cell of Cu(111) or Au(111). The metal atoms in the lowest layer of the slab were fixed, while the molecules were allowed to S6 optimise both their geometry and position; however, there was almost no change in the molecules' lateral positions during optimisation because of the very flat PES.

Additional theoretical results
3.1.

2-D structure of TPA
By scanning through a large range of b x and b y components of the b lattice vector (which is equivalent to simultaneously varying the length of b and the angle, see Figure S1), only one energy minimum was found. The minimum corresponds to a = 9.64 Å, b = 7.28 Å and  = 56.1 (b x =4.06 Å and b y = 6.04 Å). This geometry is in agreement with our previous studies [S6] and with experimental studies (cited in Ref. [S6] ). The minimum-energy 2-D structure is only 0.10 eV/molecule lower in energy than the 1-D chain of TPA. Thus, interchain interactions are very weak, much weaker than the energy cost of contracting the TPA-TPA distance along the chain to 9.3 Å or expanding it to 10.4 Å (see Figure 3 in the main text).

Potential energy surfaces of TPA on Cu(111)
PES were calculated for TPA adsorbed on Cu (111) with the molecule's long axis oriented along the direction, along the direction, and at a 15 angle to the direction.
The molecule was placed on the metal surface and moved laterally in steps of 0.25 Å to cover the whole unit cell of the Cu (111) surface. Only the metal atoms in the lowest layer of the slab were fixed; the molecule was allowed to optimize both its geometry and position, but its lateral position did not change during optimization because of the very flat PES. In this way, a 2.5 × 4.5 Å rectangular 2-D grid of TPA positions on Cu(111) was explored. Adsorption energies were calculated as described in the Methods section, and the BSSE correction obtained for adsorption on the hcp site was applied throughout.
We are not aware of other computational studies systematically exploring potential energy surfaces of similar molecules (benzene or its derivatives) on surfaces of noble metals. Several computational studies explored benzene adsorption on (111) surfaces of noble metals using DFT-D or van der Waals-corrected functionals, [S7-S9] but these were typically limited to considering a number of high-symmetry sites such as the hcp and fcc hollow, bridge and top sites. In these studies, hcp is typically found to be the lowest-energy site, although there is some variation depending on the computational method.
The PES for the orientation is presented in Figure 4 of the main text of the article.
One energy minimum can be seen in the upper right corner of the PES plot (the same minimum in all four corners, connected by the surface lattice vectors); another minimum is located in the central region of the PES. These minima correspond to the molecule adsorbed with the center of its aromatic ring above the hexagonal close-packed (hcp) hollow site, as shown in Figure 4a. The second minimum is not perfectly centeredthis is most likely caused by the finite step used for calculating the PES; it points to the limit of accuracy of our calculations (few tenths of Ångstrom and few hundredths of eV This could be a further cause of the slight asymmetry observed in the PES. The highest maxima (in the lower central region of the plot) correspond to the molecule adsorbed above a Cu atom. The face-centred cubic (fcc) hollow site is neither a minimum nor a maximum, according to our results. Figure S2 shows line scans of this PES along the y axis (i.e. along the direction).
The results for the orientation and for the orientation at a 15 angle to the direction are presented below in Figs. S3 and S4.   According to these calculations, all orientations of TPA on Cu(111) are very similar in energy, and it is not possible to identify one preferred orientation. This is in agreement with earlier computational studies (both DFT-D and van der Waals density functionals, vdW-DF) of benzene adsorption on noble metals, where the differences between the and orientations are of the order of tens of meV. [S7-S9] However, scanning tunneling microscopy (STM) images (see main text) point to the orientation as the predominant orientation; therefore, this orientation was used in our further analysis.

Potential energy surface of TPA on Au(111)
The PES for adsorption of TPA on Au (111)  S5. Figure S6 shows the corresponding PES. Adsorption energies vary between 1.19 (for the lowest-energy adsorption sites) and 1.14 eV (for the highest maximum of the PES). The lowest-energy configurations (purple regions in the PES) are not found in correspondence of high-symmetry sites but close to bridge sites, see Figure S5. The variation in the adsorption energies is however very small, leading to a corrugation of the PES of only up to 0.05 eV, in agreement with the literature (DFT-D and vdW-DF) results for benzene on Au(111). [S8, S9] Notably, this PES corrugation is much smaller than for TPA adsorption on Cu(111) (which is ~0.20 eV). Figure S7 shows that, similar to TPA on Cu(111), a cut through the PES of TPA/Au(111) along the direction of the TPA chains has a roughly sinusoidal character. The curves look less regular than those for TPA on Cu(111) in Figure S2, because of the very flat PES of TPA/Au(111) with a corrugation that is close to the accuracy of DFT. However, each of the lines in Figure S7 has a discernible wave-like shape.   Figure S6). The notation is the same as in Figure S2 for TPA on Cu(111): the y direction corresponds to the direction, i.e. the direction of the TPA chains; each curve corresponds to a different x-coordinate value. The curves have been offset by 0.05 eV for clarity.

S13
These results for the two surfaces show the generality of the periodic function (x i ) used to model the molecule-substrate interaction. The sinusoidal character is essentially given by the position of the PES minima of TPA on the surface and the situation is expected to be qualitatively similar for any (111)-terminated metal surface, as it is essentially a reflection of the geometric structure of the metal and of the molecule. We expect the same behavior also for the full 2-D PES which should be qualitatively well-described by a 2-D sinusoidal function. As a consequence, a 2-D generalization of our self-assembly model can be devised along very similar lines to the 1-D case.

Benzene and benzoic acid adsorption on Cu(111) and Au(111)
To compare the contributions of carbon (aromatic ring) and oxygen atoms to the TPA binding energy on Cu(111) and Au(111), the adsorption of benzene and benzoic acid was calculated on these two metal surfaces. The structures are shown in Figure S8, and the energies are compared in Table ST1. All the molecules were adsorbed with the aromatic ring above the hcp site on Cu(111), and above the fcc site on Au(111). Note that two adsorption configurations of benzoic acid were considered, corresponding to two different positions of the carboxylic group.  Figure S3 for Cu and in Figure S5 for Au. S15 Table ST1 Adsorption energies of TPA, benzoic acid and benzene on Cu(111) and Au(111).

Molecule
Adsorption energies on Cu (111)  It is difficult to compare these energies with published theoretical values for benzene on metals, because these latter vary rather strongly depending on the methodology chosen for treating dispersion interactions (see Ref. [S9] for an overview of published results). However, our adsorption energy for benzene on Cu(111), 0.56 eV, is close toeven if slightly lower thanthe experimental value of 0.71 eV. [S9] Our value for benzene on Au(111), 0.83 eV, is also close but slightly larger than the experimental value of 0.76 eV [S4]. [S9] These differences are determined by the C 6 parameter values taken from Ref. [S10] and Ref. [S11] for interactions with Cu and Au, respectively.
Adsorption energies have been broken down into contributions from carbon and oxygen atoms (assuming that hydrogen contributions are very small and can be grouped with their carbon atoms). The results shown in table ST1 show that carbon atoms interact with both Cu(111) and Au (111) Figure S9 shows a histogram obtained from a detailed analysis of hundreds of STM images, where the number of molecules in between the dark spots (corresponding to the long intermolecular spacings) is counted. Only images with sufficiently high resolution were used.
The histogram presents a wide distribution with an extended tail towards longer chain lengths. It should be noted that the measurement process favors counting of short chains as longer chains tend to extend beyond the edges of an STM image. Irrespectively of this, Figure S9 clearly does not show any evidence of a strongly preferred repetition length of the larger molecular separations along the TPA chains. Figure S9. Histogram counting the number of TPA molecules between the dark spots in the direction of the dimeric hydrogen bonds.