Characterizing self-assembled molecular layers on weakly interacting substrates: the role of van der Waals and the chemical interactions

For this study, a single molecule of triazine is adsorbed on a G(3 × 3) cell (see figure 2(b)). This cell corresponds to the smallest and most frequent moiré pattern observed on G/Pt(111). It allows an intermolecular distance larger than the equilibrium distance in isolated layers of triazine molecules, but the size is not enough to completely eliminate the intermolecular interaction. In order to remove this contribution and retain just the molecule–substrate interaction, we will focus on the binding energy, E_bind, defined as

$$\begin{eqnarray}&&{E}_{{bind}}={E}_{{mol}+{sub}}-({E}_{{mol}}+{E}_{{sub}}),\end{eqnarray} \tag{ 1 }$$

where ${E}_{{mol}+{sub}}$ is the energy of the full system in the G(3 × 3) cell and E_mol and E_sub are the energies of the molecular layer and the substrate calculated separately on the G(3 × 3) cell using the geometry obtained during the relaxation of the whole system. Subtracting E_mol, that includes the residual intermolecular interaction, E_bind just provides information about the molecule–substrate interaction.

Zoom In Zoom Out Reset image size

We have also calculated the standard adsorption energy, E_ads, that includes both the binding energy and the intermolecular interaction. The latter obviously depends on the relative orientation between molecules. It is given by

$$\begin{eqnarray}&&{E}_{{ads}}={E}_{{mol}+{sub}}-({E}_{{mol}}^{0}+{E}_{{sub}}^{0}),\end{eqnarray} \tag{ 2 }$$

where the new energy references, E⁰_mol and E⁰_sub, are the energies for a molecule in the gas phase and the isolated substrate respectively, both relaxed to their equilibrium configuration. The single molecule energy, E⁰_mol, is computed in a much larger cell (30 × 30 × 30 ų) to minimize any intermolecular interaction.

If a localized atomic orbital scheme is used, in particular when dealing with calculations performed with the CRYSTAL code, the adsorption energy is affected by the basis set superposition error (BSSE). In order to correct for the BSSE, the counter-poise (CP) method can be used [32]. The CP-corrected energy E^CP_bind, has then the form [33]

$$\begin{eqnarray}&&{E}_{{bind}}^{{CP}}={E}_{{mol}+{sub}}-({E}_{{mol}}^{{Gh}({mol}+{sub})}+{E}_{{sub}}^{{Gh}({mol}+{sub})}),\end{eqnarray} \tag{ 3 }$$

where ${E}_{{mol}+{sub}}$ is the computed total energy, and ${E}_{{mol}}^{{Gh}({mol}+{sub})}$ and ${E}_{{sub}}^{{Gh}({mol}+{sub})}$ are the energies of the molecule and the substrate calculated respectively in the presence of the ghosted atoms of the substrate and of the molecule. In the same way, the corrected adsorption energy E^CP_ads can be written as follows

$$\begin{eqnarray}&&{E}_{{ads}}^{{CP}}={E}_{{mol}+{sub}}-({E}_{{mol}}^{0}+{E}_{{sub}}^{0})+{\rm{\Delta }}{E}_{{mol}}^{{Gh}({mol}+{sub})}+{\rm{\Delta }}{E}_{{sub}}^{{Gh}({mol}+{sub})}\end{eqnarray} \tag{ 4 }$$

being ${\rm{\Delta }}{E}_{{mol}}^{{Gh}({mol}+{sub})}$ and ${\rm{\Delta }}{E}_{{sub}}^{{Gh}({mol}+{sub})}$ the corrected energies of the substrate and the molecule at the geometry adopted in the adsorbate-substrate system that can be expressed

$$\begin{eqnarray}&&{\rm{\Delta }}{E}_{{mol}}^{{Gh}({mol}+{sub})}={E}_{{mol}}-{E}_{{mol}}^{{Gh}({mol}+{sub})},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\rm{\Delta }}{E}_{{sub}}^{{Gh}({mol}+{sub})}={E}_{{sub}}-{E}_{{sub}}^{{Gh}({mol}+{sub})}.\end{eqnarray} \tag{ 6 }$$

We have characterized the potential energy surface (PES) profile by calculating these energies and adsorption distances for several high symmetry adsorption sites of the molecule on the G layer. Apart from the analysis on the plane parallel to the substrate (the xy-plane), we have also studied the energy variation upon changes in the adsorption distance. Starting from the geometry of the ground state of the system we have modified homogeneously the adsorption distance of the molecule letting the atoms relax but freezing the z coordinate of all atoms in the molecule. This procedure allows us to unveil the PES as a function of the molecule–substrate separation as well as to disclose the different contributions (short range, vdW or intramolecular energy changes) to both the binding energy and the energy barriers.

The study of the intermolecular interaction in an isolated layer of triazine molecules can be easily performed by a simulation with a small unit cell including just one molecule. However, the large sizes of the moiré patterns found in the experiments are out of the possibilities of the current implementations of the DFT codes. Our characterization of the intermolecular interaction among triazine molecules adsorbed on different G-terminated substrates is based on simulations of islands formed by three triazine molecules, the smallest combination needed to preserve the triangular symmetry of the SAM lattice. Our calculations have been carried out with a G(6 × 6) cell (see figure 2(b)). With this cell size, there is still a residual interaction between neighboring islands, but it is very small compared to the interaction among molecules inside the island. A larger simulation cell able to include more molecules would be superior to reproduce the intermolecular interactions of the realistic moirés. However, first, the cell size, once the substrate is included, would be too large for DFT calculations; and second, we have tested that the effect of both, the position of the island respect to the substrate and the inclusion of more neighbors (an increase of the molecule coordination) are small and do not change the main conclusions. The validity of this approach is confirmed by the comparison of the energetics and structure of a triazine monolayer with our three-molecule island: the H-bond distances are identical (2.35 Å) for both cases, and the interaction energy per H-bond are very similar (127 meV and 131 meV respectively), with just an expected small change due to the different coordination number (see inset of figure 5(b)).

Our goal is to determine the intermolecular energy as a function of the H-bond distance between the molecules in the island. Our starting point is a full relaxation of the island on top of each of the different substrates to determine the equilibrium structure. Then, we use this optimized geometry to prepare a set of configurations where we displace the molecules from their equilibrium positions in the xy-plane while preserving the symmetry of the island and its orientation with respect to the substrate. Each of these configurations is subsequently relaxed, keeping fixed the xy position of a C and a N atom from each molecule (to enforce the constraint on the H-bond distance) and allowing the rest of the atoms in the molecule and substrate to move freely.

For each configuration, we determine the intermolecular interaction among the molecules on the island, E_intermol, using the equation:

$$\begin{eqnarray}&&{E}_{{intermol}}={E}_{{interac}}-{E}_{{bind}},\end{eqnarray} \tag{ 7 }$$

where we subtract the molecule–substrate interaction, E_bind, from the total interaction energy, E_interac. The total interaction energy is calculated by subtracting from the total energy of the system, ${E}_{{island}+{sub}}$, the energy of the isolated substrate, E_sub, and the energy of the isolated molecules, ${E}_{{{mol}}_{i}}$:

$$\begin{eqnarray}&&{E}_{{interac}}={E}_{{island}+{sub}}-\left(\displaystyle \sum _{i=1}^{3}{E}_{{{mol}}_{i}}+{E}_{{sub}}\right),\end{eqnarray} \tag{ 8 }$$

where the reference energies E_sub and ${E}_{{{mol}}_{i}}$ are calculated using the geometries corresponding to the equilibrium configuration of the whole system.

As the chemical environment of each molecule in the island is different, we have to perform three independent simulations to determine the binding energy of each molecule. Each of these simulations includes only one of the molecules of the island interacting with the substrate. Its contribution to the binding energy is obtained from the total energy ${E}_{{{mol}}_{i}+{sub}}$ by subtracting the total energies of the substrate, E_sub, and the isolated molecule ${E}_{{{mol}}_{i}}$

$$\begin{eqnarray}&&{E}_{{bind}}=\displaystyle \sum _{i=1}^{3}({E}_{{{mol}}_{i}+{sub}}-{E}_{{sub}}-{E}_{{{mol}}_{i}}).\end{eqnarray} \tag{ 9 }$$

For most of our calculations, we have used the PBE functional [38] supplemented by different vdW approaches. Some calculations have been done with XC functionals like optB86b [39] (see below) that include a different exchange contribution and a kernel that accounts for vdW interactions. We have also explored the possible role of an improved description of the short range contributions using hybrid functionals, in particular HSE06 [40], PBE0 [41] and B3LYP [42, 43]. These last calculations were performed using CRYSTAL14 as the convergence with hybrid functionals is more optimized in this code. For these simulations, we carried out only static calculations using the geometry from the PBE+vdW result but optimizing the adsorption distance.

There are two ways to include this interaction in the calculations: either as a correction to the final energy or through a kernel in the electronic XC functional that includes these effects in the self-consistency process.

From the first group of vdW implementations we have used the simplest Grimme approach (PBE-D2) [27] and also an advanced implementation by the same authors (PBE-D3) [28]. For the PBE-D2 calculations, we have used the default values given in [27] for all the chemical species except for the Pt, which is not tabulated. For this metal, we have used the parameters C₆(Pt) = 20 J nm⁶ mol⁻¹ and R₀(Pt) = 1.9 Å which successfully reproduce the G–Pt distance [44]. In the case of the PBE-D3 implementation (for which we used the Becke–Johnson damping [45–47]), the vdW C₆ parameter takes into account the local environment of the atom through its coordination number—which means that it may change during the simulation—and is determined by the program. We have also performed calculations with the TS+SCS [48]—similar to D2, apart from the fact that in this approach the parameters are charge-density dependent, accounting in this way for screening effects—and the MBD approach [25, 29]—which contains both the many-body energy and the screening which are missing in simple pairwise approaches—which are the more sophisticated methods to date. From the second group—the so-called DFT-DF functionals—we have used the Klimes optB86b functional [39] which is able to partially account for screening effects and has proven to work very well for G systems [49].

Table 1 shows the variation on the characteristic distances of our systems depending on the vdW approach used. The mismatch between the C and Pt lattices in the G(3 × 3) moiré of G/Pt(111) is small (0.6%) (see table 1) and we decided to fix the size of the supercell to match the relaxed G lattice calculated with each functional-vdW scheme.

Table 1.  Values of characteristic parameters to describe the substrates that we have simulated for different functionals and vdW flavors. The first two columns correspond to the values for Pt and G lattice parameters and the third an fourth are the G-metal distance and graphite interlayer distance respectively

	A0 (Å)	dnn (Å)	$\langle {d}_{{GPt}}\rangle$ (Å)	$\langle {d}_{{GG}}\rangle$ (Å)
	Bulk Pt	Graphene	G/Pt(111)	Graphite
Exp	3.913 [50]	2.46 [51]	3.30 [52]	3.35 [53]
PBE-D2	3.953	2.458	3.345	3.229
PBE-D3	3.927	2.470	3.317	3.384
PBE-TS+SCS	3.951	2.466	3.325	3.382
PBE-MBD	3.968	2.465	3.427	3.410
optB86b	3.958	2.468	3.361	3.324

5.1.1. Triazine adsorption sites

We have characterized the preferential adsorption sites on a G(3 × 3) cell, in order to make contact with the literature. We have checked that the results for PBE-D2 in the G(3 × 3) cell are consistent to what is found in a bigger G(6 × 6) cell and, in order to minimize the computational cost, we stuck to the G(3 × 3) cell to perform the molecule adsorption characterization. Figure 3 displays the high symmetry adsorption sites considered in our study and the final orientation of the molecule with respect to the G honeycomb lattice. The molecular symmetry directions are aligned with those of the underlying G substrate for the three top and one of the bridge configurations, while they are rotated by 30° in the Cross(R) and Bridge(R) geometries in order to maximize the intermolecular interaction with molecules in neighboring cells. Adsorption energies and adsorption distances for the PBE-D2 approach are shown in table 2. The Bridge(R), followed by the Cross(R), are the most stable configurations due precisely to the intermolecular interaction, in good agreement with previous calculations [16].

Zoom In Zoom Out Reset image size

Table 2.  Characteristic binding energies and mean adsorption distances calculated for different adsorption sites of triazine on G depicted in figure 3 using the PBE [38] functional with the Grimme-D2 [27] vdW implementation. The values of the adsorption energies agree with previous works [16]. The energy differences with respect to C top (minimum energy site) for each site are shown in brackets next to the absolute energies.

Site	Ebinding		Eads		dads
	(meV)		(meV)		(Å)
C top	−404	(0)	−388	(0)	3.15
N top	−377	(27)	−361	(27)	3.20
All top	−344	(60)	−327	(61)	3.28
Bridge(R)	−367	(37)	−494	(−106)	3.21
Cross(R)	−339	(65)	−464	(−76)	3.29

In order to obtain the relevant molecule–substrate interaction for the energetic analysis of the SAMs, we need to remove the intermolecular interaction, which is still important in the G(3 × 3) cell. We follow the procedure outlined in section 3 to obtain for each adsorption site the corresponding binding energy. Table 2 shows that, in terms of the binding energy, the most favorable adsorption configuration is the C top, followed by the N top. In the case of a PAH molecule, benzene would be the comparable molecule in this case, these two most favorable positions would be equivalent—there would be no N atoms but more C atoms instead. In the case of triazine, the charge density is larger in the N atoms, on the contrary of what happens in π–π interacting molecules, which breaks the equivalence between these two sites. The C top binding energy, ∼400 meV/molecule, is similar to the cohesive energy measured for a benzene molecule on HOPG, ∼500 meV/molecule [22], which points out that adsorption in both cases is ruled by the same mechanism: non covalent interactions.

We have studied the adsorption and binding energies of the C top and N top configurations using different XC functionals and vdW implementations. Results are shown in table 3. They all predict the C top to have the largest binding energy, but significant differences in absolute binding energies and adsorption distances are found among the different approaches.

Table 3.  Characteristic binding energies and mean adsorption distances calculated for the C top and N top adsorption sites (see figure 3) of triazine on G using different vdW flavors (DFT-D2 [27], DFT-D3 [28], DFT-DF(optB86b) [39], TS+SCS [48], and MBD [25, 29]) (see section 4). The energy differences with respect to C top (diffusion energy barriers) are shown in brackets.

Functional	Site	Ebinding		dads
		(meV)		(Å)
PBE-D2	C top	−404	(0)	3.15
	N top	−377	(27)	3.20
PBE-D3	C top	−400	(0)	3.30
	N top	−381	(18)	3.34
PBE-TS+SCS	C top	−574	(0)	3.43
	N top	−551	(22)	3.49
PBE-MBD	C top	−360	(0)	3.30
	N top	−342	(17)	3.35
optB86b	C top	−562	(0)	3.19
	N top	−537	(26)	3.25

The PBE-D2 is the one with lower adsorption distances although it yields similar binding energies to PBE-D3 (∼400 meV). The rest of the vdW implementations place the molecule much farther from the substrate. This second approach—PBE-D3—along with the PBE-MBD method result in very similar adsorption distances and, although the binding energies are smaller for the MBD approach, the energy differences are very similar for these two methods. The PBE-TS+SCS and optB86b schemes produce bigger binding energies than the other methods (∼550 meV). However, while the optB86b adsorption distances are somewhere in between the PBE-D2 and PBE-D3 cases, the PBE-TS+SCS is the method that results in the farthest distances between molecule and substrate.

5.1.2. Energy barrier calculation

The calculation of the diffusion energy barrier for triazine on G requires to characterize the PES and to determine the minimum energy path. The results of this study are shown in figure 4, where a triangular path that goes through the high symmetry adsorption positions framed in black in figure 3 is explored. These results suggest that the C top-N top-C top path yields the lower energy barrier. We have confirmed this point by performing calculations with the PBE-D2 approach using the nudged elastic band method [54, 55], and, thus, we restrict the study to determine the energy barrier with the other XC+vdW approaches to the path between these two configurations.

Zoom In Zoom Out Reset image size

According to the results shown in table 3, the biggest barriers are obtained for the cases where the D2 vdW implementation is used. This is due to the known overestimation of the dispersion forces associated with this vdW approach. This behavior pushes the molecule closer to the surface towards a more repulsive chemical interaction range which results into a bigger barrier. The rest of the methods used yield to smaller but similar barriers on the order of 15–25 meV.

In summary, diffusion energy barriers are always underestimated by the calculations regardless of the method used, similarly to what was found in benzene on G. For that case, some DFT approaches predict a barrier of less than 10 meV [24], which is small compared to the experimental measurements (17 ± 12 meV) [56], although compatible with them given the large experimental uncertainty.

In order to understand the interactions that are controlling the molecule–substrate interaction, we have determined the adsorption energy and the short range and vdW contributions as a function of the adsorption distance for the triazine on the C top and N top sites using the PBE-D2 functional (figure 4(b)). The molecule adsorbs at a distance of 3.15 Å (3.20 Å) with respect to the G on the C top (N top) configuration. We can recognize a very similar behavior to what is found in π–π interacting systems: the driving attractive force is the vdW—which is very similar in the two cases—while the short range contribution is responsible for the PES corrugation and determines the difference in equilibrium adsorption distance (∼0.05 Å) between the two sites. Notice that the minimum of both curves is located in the repulsive region of the short range interaction due to the effect of the vdW contribution that pushes the molecule closer to the surface. We will come back to these energy curves in the discussion presented in section 7.

In order to analyze the role of the intermolecular interaction in the SAM formation, we first optimize the structure of the system including both the 3-molecule island and the G substrate on a G(6 × 6) unit cell. The resultant equilibrium geometry is depicted in figure 2(b). The molecules are not exactly on the most stable adsorption configuration for a single adsorbed triazine, the C top site, as the interaction with other molecules induces a small displacement on each molecule from the C top configuration of around 0.3 Å.

We have characterized the strength of the intermolecular interaction following the procedure described in the section 3.2. Figure 5 shows the different energy contributions as a function of the separation (defined as the length of the N–H hydrogen bond) between the three molecules of the island. In figure 5(a), we plot the total interaction energy (E_interac), the binding energy (E_bind)—the sum of the interaction energies of each molecule with the substrate, as defined by equation (9)—and the intermolecular energy calculated as ${E}_{{intermol}}={E}_{{interac}}-{E}_{{bind}}$. Figure 5(b) compares the intermolecular energy per H-bond of the isolated island with the energy of the island adsorbed on G. Although both interactions are very similar, the substrate has the effect of slightly weakening the intermolecular interaction, reducing the H-bond energy from ∼−131 to ∼−119 meV, and increasing the H-bond distance from 2.35 to 2.36 Å.

Zoom In Zoom Out Reset image size

6.1.1. Triazine adsorption sites

6.1.2. Energy barrier calculation

Finally, we analyze the effect of the substrate on the intermolecular interaction. We have followed the same procedure applied for the adsorption on G. Figure 8 presents the energy per H-bond as a function of the bond length for G, graphite, G/Pt(111) and the isolated three-molecule island. Experiments show a change on the intermolecular distance of ∼10 pm between the triazine molecules adsorbed on graphite and on G/Pt(111). Calculations yield distance variations between substrates one order of magnitude smaller (<1 pm) and almost identical bond energies (differences less than 1 meV). Therefore, although the presence of G does induce some small changes in the intermolecular interaction, the substrate underneath seems not to play any role.

Zoom In Zoom Out Reset image size

For all the systems considered, G, graphite and G/Pt(111), the clearly preferred configuration for the adsorption of a single triazine molecule is the C top site. The orientation of the triazine molecule with respect to G (the relative angle between the directions that are equivalent under the three-fold symmetry of the molecule and the honeycomb lattice) defined by the C top configuration is consistent with the molecular orientation observed experimentally on the triazine SAMs formed both on HOPG and G/Pt(111). This suggests that the molecule–substrate interaction, although weak, plays a key role in fixing the SAM orientation with respect to the G.

One would be tempted to conclude that the molecule–substrate interaction controls the SAM formation. If this was the case, we should observe supercells like the $G(\sqrt{7}\times \sqrt{7})$ shown in figure 9(a). This configuration accommodates all the triazine molecules close to C top sites in order to maximize the substrate–molecule interaction, although they are slightly rotated from the perfect C top configuration to favor the correct alignment of the H-bonds. Although retaining the proper orientation, this favorable substrate binding configuration imposes intermolecular N–H distances of ∼2.7 Å that are large compared to the intermolecular distance in the isolated triazine monolayer (2.35 Å). This results in a high cost in intermolecular binding: using the results of figure 8, the intermolecular energy in the $G(\sqrt{7}\times \sqrt{7})$ moiré is reduced by 3 H-bonds/molecule × ∼23 meV/H-bond = 69 meV/molecule.

Zoom In Zoom Out Reset image size

The analysis above points out to a dominant role for the intermolecular interaction. In fact, our calculations show that the strength of molecule–substrate and intermolecular interactions are similar. For G, for example, the binding energy in the C top position (−404 meV with PBE-D2) is of the order of the H-bond contribution per molecule (6/2 H-bonds/molecule × −119 meV/H-bond = −357 meV/molecule). The SAM formation is ruled by the subtle balance between changes in both the binding energy and the intermolecular interaction that are controlled by (i) the corrugation of the binding PES, and (ii) the distance dependence of the intermolecular interaction, and result in large moirés. The latter is the key factor according to our analysis of the equilibrium configuration of the three-molecule island on all of the substrates: we found lateral displacements with respect to the C top site in order to keep the intermolecular distances very close to those of the isolated triazine monolayer. This is also consistent with the experimental results. In the large moiré patterns observed in the experiments, the intermolecular distances are close to its optimum value—with an intermolecular energy loss of less than 3 meV per molecule (see figure 8)—while the molecules are not placed on optimal C top adsorption configurations but distributed along the binding PES that has a corrugation around 60 meV (figure 4). Assuming a uniform distribution, this would result on an average binding energy loss of ∼30 meV/molecule, that is compensated by the more favorable intermolecular interaction. This subtle balance explains the preferential stability of large moiré patterns. Notice that for this system we predict a variation of the adsorption distances of ∼15 pm inside these large moirés (figure 4) and that molecules on non-highly-symmetric sites should be slightly tilted. These features could be confirmed in future STM experiments.

Having understood the role of the molecule–substrate and intermolecular interactions, we address the possible reasons behind the discrepancy found between our calculations and the experimental evidence about the influence of the substrate on the SAM formation. The different moiré patterns found in the experiments for HOPG and G/Pt(111) can be explained phenomenologically in terms of the mismatch between the different lattice periodicities involved [59]. We have analyzed the experimental STM images in order to determine the moiré periodicity (see figures 9(b) (c)), the commensurability and the relative orientation of the molecular and substrate lattices. For graphite, the images can be reproduced with a $2\sqrt{67}\times 2\sqrt{67}-R12.2^\circ$ periodicity (lattice parameter 40.3 Å), where the molecular layer is rotated 19.8° with respect to the graphite lattice (see figures S1 and S2 on the supporting information). The same analysis for the G/Pt(111) case, gives a periodicity of $2\sqrt{79}\times 2\sqrt{79}-R43^\circ$ (lattice parameter 43.7 Å), and a relative angle of 21.2° (see figures S3 and S4 on the supporting information). In graphite, the lattice mismatch between the molecular layer and the G lattice is very small (−0.02%). When Pt is introduced, a new lattice periodicity has to be taken into account. If the moiré pattern for G-triazine were preserved, there would be a mismatch of more than 2% between the metal and the G-triazine system. Thus, it is favorable for the system to choose a different moiré pattern with a unit cell that slightly increases the mismatch of the triazine monolayer with the G (0.05%), and consequently the intermolecular distance, but considerably reduces the mismatch with the Pt underneath (0.2%).

The analysis presented above clearly shows the effect of the Pt support through the interaction of the molecules with the Pt atoms either directly or through the modulation created in the G by the Pt surface. Our calculations should capture this effect: the SAM should be trading the loss of some intermolecular bonding energy to accommodate the adsorption positions to the Pt lattice or to the G(3 × 3)/Pt(111) moiré. Remember that the only trace of the Pt in our calculations on this substrate is a difference of ∼2 meV/molecule in the binding energy among equivalent adsorption sites located on different areas of the G(3 × 3)/Pt(111) moiré.

A simple estimate of the energy difference between the moiré patterns observed on graphite and G/Pt(111) does not support the preferential stability of the latter for the G/Pt(111) case: the larger intermolecular distances would result on a reduction of the intermolecular energy of ∼2 meV/H-bond×3 H-bond/molecule ∼6 meV/molecule (figure 8), that cannot be compensated by the ∼2 meV/molecule adsorption energy corrugation of the G/Pt(111) moiré. This result is consistent with the theoretical underestimation of the diffusion barriers and seems not to be essentially modified by the use of other XC functionals and vdW implementations.

The evolution of the diffusion barriers with the adsorption distance (figure 10) suggests that we can increase significantly those barriers if we push the triazine molecule towards the G layer: a rigid reduction of the adsorption distance for all the sites by ∼0.2 (∼0.3) Å for graphite (G/Pt(111)) would bring those barriers to agreement with the experimental values. Thus, our calculations point out to an underestimation of the molecule–substrate interaction, larger in the presence of Pt, as a possible explanation for the discrepancies with experiments in the value of the diffusion barriers and the influence of the G support.

Zoom In Zoom Out Reset image size

We can analyze the implications of this idea by exploring the evolution of the N top-C top energy difference versus the adsorption distance as is shown in figure 10 for the three different substrates using PBE-D2. According to this test, we can predict the adsorption distance ranges for which the experimental barriers would be reproduced. At the new predicted equilibrium adsorption distances (different for each site), the binding PES for the G/Pt(111) case has a larger corrugation (difference in binding energy among equivalent sites on different areas of the G(3 × 3) moiré) needed (∼10 meV/molecule, see the inset on figure 10) to compensate the intermolecular energy loss (∼6 meV, see above), required to accommodate the triazine SAM to the Pt lattice and stabilize the $2\sqrt{79}\times 2\sqrt{79}-R43^\circ$ moiré observed on the G/Pt(111) case.

Our results clearly show that the effect of the short range interaction should be larger in order to get a good comparison with the experimental evidences. The first candidate to explain this would be that the attraction, provided by the vdW dispersion forces, yielded by our calculations is being underestimated: an increase of the vdW contribution would decrease the adsorption distances towards the values that are required to reproduce the experiments.

Notice that what it is needed is an increase of the vdW force, as it is the one that should compensate the repulsive short range force to fix the adsorption position of the molecules.

Nevertheless, the most accurate vdW implementations that we have tried, D3 and MBD approaches, show similar forces (see figure 11) [13]. Even the D2 scheme, which is well-known to overestimate the vdW contributions [28, 29], is unable to reduce the adsorption distances enough to reproduce the experimental data. Indeed, it would be needed to increase the vdW-D2 by a factor of 2 (3) for the graphite (G/Pt(111)) substrate to get to the required adsorption distances. The good agreement between the most reliable vdW implementations, that have been well tested in other systems [13, 25, 28, 29], lead us to think that assuming errors on the calculation of this contribution of 200%–300% is not reasonable.

Zoom In Zoom Out Reset image size

After discarding the vdW as the prime suspect, we need to consider other sources of disagreement. The next candidate left is the short range interaction. In order to characterize the system with more sophisticated XC functionals we have used the CRYSTAL code which is more optimized than VASP for simulations with hybrids.

To check the comparison between the two codes we first perform a benchmark calculation with PBE-D2. In the top panel of figure 12 the results of the calculations with VASP and CRYSTAL are compared. It can be observed that, as expected, the vdW contribution matched exactly between the two codes (the D2 implementation is purely geometrical) while there is an evident difference in the short range interaction. This discrepancy can be explained in terms of the basis set used: the decay of the Gaussian wave functions of CRYSTAL is more steep than that of the plane waves of VASP. As a consequence, the equilibrium adsorption distances are pushed farther from the substrate in the case of CRYSTAL than in the case of VASP (see figure 12). Once we have clarified this difference between codes let us analyze the results obtained with CRYSTAL for hybrid XC functionals.

Zoom In Zoom Out Reset image size

Three different well-tested functionals have been used in these calculations (HSE06, PBE0 and B3LYP). In the lower panel of figure 12 we compare the N top-C top energy differences between the different functionals including the results for PBE obtained with both VASP and CRYSTAL—the difference between codes previously discussed is evident too. The obvious conclusion from this figure is that there is no significant change in the energy differences originated for the use of different XC functionals: the energy differences are larger than the one calculated with VASP but all the CRYSTAL results are very similar among them.

We are, then, unable to address the source of our lack of interaction because the two main suspects—the vdW and the chemical interactions—separately seem to be not enough to explain the discrepancies with the experiment.