Simulating adsorption of complex molecules using the linearity between interaction energies and tunnelling currents: the case of hexabenzocoronene on a Ag/Pt dislocation network

We present a method for determining the adsorption position of complex molecules on surfaces from first-principles calculations. We use electron transport theory through a vacuum barrier (theory of scanning tunnelling microscopy), and the relationship between the tunnelling current and the interaction energy between surface and tip. This method is especially useful for obtaining reasonable adsorption positions of relatively large molecules (>50 atoms) on reconstructed surfaces, e.g. dislocation networks. The main advantage is that this approach is computationally efficient and does not require mapping of all adsorption sites by separate simulations. The reliability is illustrated by simulating the adsorption of hexabenzocoronene on a model Ag/Pt(111) surface.


Introduction
Molecular self-assembly has attracted considerable interest in recent years [1,2]. The aim is to design and build molecular electronic devices in a controllable and reproducible way. This is essential for mass production of computer elements, thus it has a major impact on the computer industry. A large variety of supramolecular assemblies on solid surfaces has been reported in the literature so far [3]- [5] due to the diversity of molecules and their combinations. Fortunately, the molecule-molecule interaction can be controlled by appropriate functionalisation of the building blocks [6]. As was recently shown, supramolecular structures can be formed driven by molecule-substrate interaction [7]- [9]. In this case, a substrate whose periodicity is comparable to the supramolecular structure is needed. Ideal substrates for this purpose are partial dislocation structures, e.g. Ag/Pt, which were studied extensively in the last decade or so [9]- [12]. In this case, molecule-substrate interactions play the key role which is the focus of this paper. We present a theoretical approach to determine the adsorption position of molecules on noninsulating surfaces. The paper is organised as follows: section 2 describes the applied method and section 3 shows an example of how the approach works applied to the adsorption problem of a hexabenzocoronene (HBC) molecule on a model Ag/Pt(111) surface.

Method
Determining adsorption positions or rather adsorption energies of small molecules on simple metal surfaces from computer simulations is becoming less time consuming due to the improvement of computing power in recent decades. The most advanced methods are parameterfree and are based on density functional theory (DFT) [13,14]. However, calculating big molecules on complicated surfaces has a great computing time requirement. This motivated us to develop a more efficient method based on the theory of electron tunnelling [15]- [17] of which a small summary is found below.
Employing the multiple scattering approach [18] to first-order, the leading term of the tunnelling current between a surface and tip can be written as Here, f denotes the Fermi-function, V is the bias voltage between surface and tip, κ k and κ i are vacuum decays for surface and tip, respectively, and M ik is the usual Bardeen matrix element [19], with real-space wavefunctions ψ and χ for surface and tip, respectively, resulting from DFT calculations, we used the VASP code [20]. Here, the surface integral is taken on the separation surface between surface and tip. Taking the zero bias limit the tunnelling current results at the 3 Bardeen current [18,19], It has to be emphasised that the multiple scattering treatment is necessary for high bias voltages whereas the Bardeen method is a good approximation for the small bias regime. The first-order interaction energy between surface and tip is [18] Interactions between a molecule and the surface can be determined in the following way. First, the molecule has to be attached to the tip which we will refer to as modified tip. We calculate the wavefunctions of the modified tip (χ) and the surface (ψ) separately using DFT, in our case the VASP code [20]. Next, we calculate the interaction energy (E int ) between the molecule and surface at a single (x 0 , y 0 , z 0 ) point without geometrical relaxation, and using the fact that the interaction energy is proportional to the tunnelling current [21], the whole interaction map E int (x, y, z) can be obtained by calculating the current map I (x, y, z) between the surface and modified tip at close to zero bias voltage (V → 0) and employing the equation The maxima of the interaction map reveal the most likely adsorption positions of the molecule on the surface. It has to be noted that this is a very rough approximation since it is assumed that the surface and tip are weakly coupled which is obviously not the case by calculating interaction energy between the molecule and the surface where even chemical bonds can be formed. However, this will change the findings obtained from the weak coupling limit only if the electronic structure of the surface and the adsorbed molecule are substantially changed due to interactions, which is not the case for large organic molecules composed of π-systems. The aim of the method is to obtain a reliable starting point for further analysis, i.e. geometrical relaxation of the molecule. Both from a theoretical point of view, and in experimental work, the relation of the tunnelling current to the interaction energy has been widely discussed, since the original paper deriving a linear relationship was published [21]- [23] (and references therein). At present, experimental evidence suggests that it is close to linear, if only short-range chemical interactions are considered. Since these interactions describe the current variation due to the electronic surface structures, they are also the most important interactions regarding the chemical bonding between a surface and a molecule. We stress here, that the estimate, gained from current simulations, has to be checked for consistency by actual simulations. In this sense, the value of the method is due to the limitation of the adsorption sites, one has to simulate by DFT to a manageable number of sites. While this relation is consistently formulated in terms of Green's functions of the interface and electron scattering, we could not verify the conflicting claim of Julian Chen, that the relation is not linear [22]. Part of the problem, in our view, is the unclear formulation of the proposition from various-and not necessarily consistent-simple models and assumptions based on elementary quantum mechanics. However, even though this is an important issue, it is not the focus of the present paper since it does not influence our presented method here because for this method to succeed it is only necessary that interaction energy and current are monotonously related. The requirement of monotonicity is sufficient and ensures that the maxima of the current map will define the maxima of the interaction map, thus the adsorption positions of molecules can be determined. The above procedure gives us a qualitative picture of the adsorption positions practically on any surface, even on dislocation structures. The simulated adsorption geometries can then be analysed in more detail, i.e. by a full geometrical relaxation of the molecule at the selected positions. The only disadvantage is that the molecule attached to the tip is in a fixed orientation. Changing the orientation induces an additional degree of freedom in the simulation. Therefore, this method is less efficient for functional molecules with high degrees of asymmetry. Assuming reduced degrees of freedom for the molecule, several possible orientations can be identified and modified tips formed accordingly. A search for the maxima of the interaction map can then be performed using modified tips with different molecular orientation as written above then full geometrical relaxation can be carried out. In spite of the above mentioned restrictions on the applicability of the method it is still less time consuming than the straightforward geometrical relaxation of the molecule on the surface, especially if the surface is reconstructed and the molecule is big. We tested the reliability of the method for a flat molecule, see section 3.

HBC molecule on Ag/Pt(111)
Hexa-peri-HBC (C 42 H 18 ) is a well characterised molecule [24] and a good candidate for a building block in self-assembled systems due to its ability to act as molecular anchor on template surfaces with D3h and D6h symmetry [25,26].
The method described in section 2 has been used to determine where the HBC molecule is most likely to adsorb on a model Ag/Pt(111) dislocation structure (1.64 ML Ag coverage, 2D pattern), see figure 1. Here, green circles denote Pt atoms, light blue circles form the full first monolayer of Ag and dark blue circles are Ag atoms in the topmost layer forming bigger and smaller triangles in fcc and hcp sites, respectively. We used different blue colours for visibility purposes only, light and dark blue coloured circles denote Ag atoms. The size difference of fcc-and hcp-stacked triangles is explained by the energy difference between fcc and hcp sites in favour of fcc. This is consistent with the earlier studies [11,12]. Our model set-up was proposed by Brune et al [10], namely that the experimentally observed structure was modelled by fcc and hcp regions in a trigonal network. For computational limitations the full unit cell of 16 × 16 [10] or 24 × 24 [9] could not be calculated. Instead, we simulated a 10 × 10 structure containing 464 atoms. A recent study [9] questions Brune's model [10] for the Ag/Pt(111) dislocation network as careful STM measurements show that the different stacking regions in the topmost layer are almost in-plane and the dark discommensuration lines cannot originate from there, rather from the layer below. Moreover, slight lattice parameter variations were observed in the uppermost Ag layer [9] which could also have a big impact on the structure. This type of dislocation formation could be related to other systems as well [27] but the full understanding is still missing. It has to be emphasised that the aim of our study is to show the applicability of the method in section 2 therefore we used the simplest model as proposed by Brune et al [10].
We considered three different adsorption sites instead of the one described above (x 0 , y 0 , z 0 ) in order to show the reliability of this approach. All chosen sites have three-fold symmetry, hollow: intersection of discommensuration lines, middle of small and big triangles, and they are shown in figure 2. We calculated the interaction energy between HBC and the template for all configurations. These values were obtained as E = −(E Surface + Molecule − E Surface − E Molecule ) using a set of three independent total energy calculations of the adsorbed system (surface + molecule) and the surface and molecule separately. These interaction energy values will be later compared to the tunnelling currents extracted from a three-dimensional current map using the modified tip, under low bias conditions (we used 500 mV, which is enough to obtain electronic transitions from the lowest unoccupied orbitals of the molecule attached to the tip). Let us now calculate the current map. The HBC molecule was attached to a flat tungsten surface in the orientation of figure 2. The scan area of the clean template surface, the topography of it (W tip) and the topography with the modified tip (W tip + HBC) are shown in figure 3. As can be seen, the surface topography significantly differs from the topography measured with the modified tip. We find that the maxima of the latter map represent the most likely adsorption positions. In our case, the hollow position (left part of figure 2) seems to be the best adsorption site. In order to prove this, we compare the calculated interaction energies ( E) with the currents (I ) in the three considered configurations in table 1. As revealed by the interaction energies and the calculated current values, the hollow site is the preferred adsorption position. Comparing the calculated interaction energies with the extracted current values, it can be stated that their fraction is fairly constant. This proves the reliability of the approach.

Summary
We presented a method which allows to identify adsorption positions of molecules on surfaces using the electron tunnelling theory. This method is especially useful for big molecules with reduced degrees of freedom or complicated surfaces or their combinations. The main advantage is time saving of calculations compared to straightforward geometrical relaxation of the molecule. Disadvantage is the variety of molecular orientations which increases with increasing degrees of freedom of the molecule. This method could contribute to speed up theoretical investigations of self-assembled molecular systems on surfaces.