Site-selective adsorption of phthalocyanine on h-BN/Rh(111) nanomesh

Experiment and computer simulations were conducted in order to study the adsorption of the phthalocyanine molecules H2Pc and CuPc on the h-BN/Rh(111) nanomesh. We combine STM investigations with the exploration of the potential energy surface as resulting from density functional theory calculations. Both approaches indicate a pronounced adsorption selectivity in the so called pore regions of the h-BN nanomesh, whereas the adsorption energy landscape in the pore turns out to be very shallow. This is seen by the inability to image the molecule stably at 77 K by scanning tunneling microscopy. Understanding the nature of the binding by rationalizing the site-selectivity and the mobility of the molecules is quite a challenge for both experiment and theory. In particular, we observe that the choice of the functional in the DFT description is crucial to be able to discriminate among adsorption sites that are very close in energy and to resolve low energy barriers. Our study reveals how the shape of the corrugated h-BN layer is the dominant factor that determines the subtle features of the potential energy surface for the adsorption of phthalocyanine.


A. Structural parameters
The modulation of the height of N atoms, of the height of the Rh atoms of the topmost layer of the slab, and of the BN bond length have been calculated along the [21] diagonal of the nanomesh (NM) unit cell. The three di↵erent line colors indicate the results for the three tested models, revPBE-D3 (black), PBE-rVV10 (blue), and vdW-DF (red).  Table II) is much larger than for the two other functionals. We can also see that the pore obtained with the vdW-DF functional is rather shallow (smallest values for z max B and z max N ), which can be explained by the distances of the B and N atoms from the Rh surface which are much larger than for the two other functionals. Actually it is well known that the tendency of vdW-DF is to overestimate the equilibrium distances. 40 Considering the number of N (or B) atoms belonging to the pore and wire regions is a more quantitative way of defining the size of these regions. In order to classify a N atom as belonging to the pore or wire region, we consider its coordinate z N with respect to the average coordinate z Rh of the Rh atoms belonging to the last Rh  Table II (see text). The Fermi energy is set at zero. layer. If |z N z Rh | < z lim,1 , then the N atom is considered to belong to the pore, while if |z N z Rh | > z lim,2 , then the N atom is considered to belong to the wire. From Fig. 2(a) we chose (with some degree of arbitrariness) z lim,1 = 2.5Å and z lim,2 = 3.0Å for revPBE-D3, z lim,1 = 2.4Å and z lim,2 = 3.7Å for PBE-rVV10, and z lim,1 = 3.3Å and z lim,2 = 3.8Å for vdW-DF. The number of N atoms belonging to the pore (N pore N ) and wire regions (N wire N ) are shown in Table II, where we can see that the value N pore N = 77 for PBE-rVV10 is much larger

B. Projected density of states
The (normalized) density of state (DOS) of the h-BN/Rh(111) nanomesh is shown in Figs. 2(a) and 2(b) for the pore and wire regions, respectively. It can be seen that the two functionals lead to DOSs which are very similar, in particular for the pore region. It has already been shown in Refs. 1-3 that the -band splitting of about 1 eV measured experimentally 4 could be well reproduced by the present model of nanomesh (i.e., 13 ⇥ 13 h-BN monolayer on top of 12 ⇥ 12 Rh(111) surface). This shift, coming from the shift between the N-p DOSs of the pore and wire regions can also be observed in the present work.

II. CALCULATIONS OF H2PC IN GAS PHASE
The structure of the molecule optimized in gas phase is planar, since all C-C and C-N have sp2 character. Bond lengths and bond angles do not vary significantly by changing among the three models revPBE-D3, PBE-rVV10, and vdW-DF. (revPBE) 30,31 supplemented by the D3 dispersion correction of Grimme. 32 Two nonlocal vdW functionals were considered: vdW-DF and PBE-rVV10. In vdW-DF, 33 the semilocal component consists of the exchange part of revPBE and the local density approximation for correlation. 34 The dispersion term in vdW-DF was proposed in Ref. 33 and consists of a double integral and is therefore nonlocal and leads to calculations which are more expensive than with pure semilocal functionals. In PBE-rVV10, the semilocal component is represented by PBE, 30 while the nonlocal term is of the same type as for vdW-DF, but has a di↵erent analytical form. 35,36 The STM images were generated for the optimized structures at di↵erent adsorption sites based on the Terso↵-Hamann approximation 37,38 at the isocurrent surface, and details can be found in Ref. 15. To evaluate the adsorption strength, the binding energy is calculated for each adsorption site, by using the total energy of the system, minus the sum of the energies of nanomesh and adsorbate optimized in free-standing state, i.e. E b = E total E NM E ad . The basis-set superposition error was corrected using the Boys and Bernardi counterpoise method 39 for all binding energies reported in this work. Periodic boundary conditions were always applied.

A. Isolated H2Pc
In this section, the results obtained for the H 2 Pc molecule with the functionals revPBE-D3, PBE-rVV10, TABLE I. Bond lengths in H2Pc calculated with various functionals. (C-C)1 and (C-C)2 are the averages of the C-C bond lengths belonging to the benzene ring and outside the benzene ring, respectively. (C-N)1 is the average of the C-N bond lengths belonging to the isoindole group, while (C-N)2 is for those connecting two isoindole groups. The values are inÅ. [in meV/(BN pair)] is the total-energy di↵erence between flat and corrugated (as in the nanomesh) 13 13 h-BN monolayer (the values correspond to corrugated minus flat).
ERh-BN (in eV) is the interaction energy between Rh (four layers) and h-BN monolayer (the geometries of isolated Rh and h-BN were kept fixed at as they are in the nanomesh).  and vdW-DF are briefly compared. Since all the C and N atoms in the H 2 Pc molecule ( Fig. 1) are of sp 2 hybridization, the molecule is flat. Geometry optimization of the molecule in the gas phase gives the results shown in Table I, and as we can see the three functionals give similar results. On average, the di↵erences between the functionals is at most 0.01Å. Note that a C-C bond length is shorter (about 1.40Å) when it belongs to the benzene ring than when it is outside the ring (1.46Å).

B. h-BN/Rh(111)
Turning now to the h-BN/Rh(111) nanomesh, the results obtained with the various functionals are shown in Table II as well as in Figs. 2 and 3. Figure 2(a) shows the height (with respect to the Rh surface) of the N atoms along the long diagonal of the nanomesh unit cell. We can see that the size of the pore is the largest with the PBE-rVV10 functional and that FIG. 4. H2Pc structure optimized in gas phase: H cyan, C yellow, N grey.

III. CALCULATIONS OF H2PC ADSORBED ON THE NANOMESH
The calculations have been restricted to the case of H2Pc, because CuPc requires spin polarized electronic structure optimizations, which makes the convergence slower and increases significantly the computational costs. Every structure optimization involves a simulation cell containing almost 1000 atoms. The system is metallic, which makes the wavefunction optimization slowly converging; typically about 50 iterations of the self consistent cycle are required. The shallow potential energy surface for the adsorption of the molecule makes it di cult to locate the minimum. A typical structure optimization procedure requires the generation of more than 100 configurations. The most stable configuration for the revPBE-D3 model is with H 2 Pc centered in the pore. In order to understand why the preferred adsorption site with revPBE-D3 is the center of the pore while with PBE-rVV10 is o↵-center, we compare how the DFT and dispersion contributions change. With revPBE-D3, E DFT is 0.37 eV less repulsive for poreA than for o↵A1, while E dis is only 0.31 eV less attractive. The di↵erent contributions almost compensate eachother, which makes a 0.07 eV lower adsorption energy for poreA. On the other hand, with PBE-rVV10, the change in the repulsive DFT term between poreA and o↵A1 is only 0.046 eV in favor of the former, while the gain in dispersion of o↵A1 is of 0.227 eV. This can be explained with the fact that at the center of the larger PBE-rVV10 NM pore, the molecule minimize the interaction with the atoms of the rim, and all its atoms are equally distant from the closest BN pairs. Moving the molecule towards the rim there is enough space to let the attractve dispersion interaction to increase, between one lobe of the molecule and the rim, still keeping the repulsive term small. With the revPBE-D3, the molecule is almost as large as the NM pore, and already when centered in the pore, some of its atoms are already quite close to the rim. Pushing the molecule further towards the rim brings one lobe almost beyond the it, thus increasing the attraction term between C atoms and BN pairs, but also inducing an equal increase of the repulsive contribution. In order to confirm the implications of the pore shape, H 2 Pc has been optimized with revPBE-D3 on the fixed NM geometry as obtained with PBE-rVV10. Similarly, we have optimized H 2 Pc with PBE-rVV10 on the fixed revPBE-D3 NM. It turns out that with the PBE-rVV10 NM shape the o↵ center adsorption is always favorite, irrespective of the functional used to calculate the interaction.     With the vdW-DF functional the adsorption energy of H 2 Pc in the does not change significantly depending on the specific site. It is always about 3.06 eV. Due to the choice of the functional and the less corrugated h-BN layer, the distance between BN and molecule is larger and does not change significantly by moving the molecule over the pore. Hence, the two terms of the interaction energy change also little.

B. Electronic structure of adsorbed H2Pc
The charge density di↵erence ⇢(r) is calculated as where ⇢ pc (r) and ⇢ NM (r) are the charge densities of the two subsystems taken alone, but in the same geometry of the complex. These volumetric data structures represent the changes in the electronic charge distribution due to the interaction between molecule and substrate. The rather small ⇢(r) computed for all pc/NM configurations confirm that no chemical bonding is formed. However, some polarization e↵ects are present, as it can be observed in the plots of ⇢(r) isosurfaces displayed in figure 6 for three configurations, poreA, o↵A1 and o↵B1. The red isosurface corresponds to accumulation of charge, and it is mainly located in the area below the pyrrole-like subunits, indicating polarization of the molecular charge towards the substrate. The corresponding charge depletion is illustrated by the green isosurface and it is localized close to atomic centers of the molecule (mainly H and N) and at some N and B of the NM pore underneath the molecule. The electron distribution at the metal, instead, is barely a↵ected by the presence of the molecule. The dipole moment of the H 2 Pc molecule in gas phase and adsorbed on the surface has been calculated using the method of the maximally localized Wannier function (MLWF) ? . Once the localization procedure has converged, the MLWF are associated either to the molecule or to the substrate. The center of the MLWF , i.e. the Wannier center, represents the position of the localized electron pairs, charge -2. The atomic positions are insted the location of the positive core charges, i.e. +5 for N, +4 for C, and +1 for H. The dipole moment is then obtained from the sum over all the particles, ionic cores and Wannier centers, of the product of the position vector times the corresponding charge In gas phase, all the three Cartesian components of the dipole are null. The interaction with the substrate, instead, induces a molecular dipole, with the largest component in the directio perpendicular to the surface, µ z . We have computed the projected density of states (PDOS) for H 2 Pc in gas phase and adsorbed in the o↵A1 configuration. In figure 7 (a) and (b), the black curves are the PDOS on C atoms and N atoms, respectively, where the energies are reported with respect to the center of the homo-lumo gap. The isosurfaces orrespondint to the wavefunctions of the homo-1, homo, lumo and lumo+1 states are diplayed in panel (c). The same PDOS computed for the molecule adsorbed in the o↵A1 configuration are reported in red. In this case the energies are given with respect to the Femri energy. It is noticed that the lumo of H 2 Pc is almost aligned with the Fermi energy of the 7 complex. Otherwise, the structure of the molecolar orbitals seems not to be significantly perturbed by the interaction with the substrate.

9
C. Mobility within the pore At higher temperature the molecules cannot be imaged stably anymore, but that they are pushed around in the pore during the STM scanning. The "double cross" configuration seen the figure originates from the bi-stable position of the molecule in the pore and not from two molecules sitting in the same pore. This later situation is realized at higher molecule coverage as can be seen in the figure 8, where the two molecules in the pore are well resolved. Depending on the tip condition and tunneling parameters the images at 77 K can become very blurred due to the mobility of the molecules in the pore. Despite the mobility the molecules are not pushed out of the pores.