Theoretical study of PTCDA adsorbed on the coinage metal surfaces, Ag(111), Au(111) and Cu(111)

A thorough understanding of the adsorption of molecules on metallic surfaces is a crucial prerequisite for the development and improvement of functionalized materials. A prominent representative within the class of π-conjugated molecules is 3,4,9,10-perylene-tetracarboxylic acid dianhydride (PTCDA) which, adsorbed on the Ag(111), Au(111) or Cu(111) surfaces, shows characteristic trends for work-function modification, alignment of molecular levels with the substrate Fermi energy and binding distances. We carried out density functional theory (DFT) calculations to investigate to what extent these trends can be rationalized on a theoretical basis. We used different density functionals (DF) including a fully non-local van der Waals (vdW) DF capable of describing dispersion interactions. We show that, rather independent of the DF, the calculations yield level alignments and work-function modifications consistent with ultra-violet photoelectron spectroscopy when the monolayer is placed onto the surfaces at the experimental distances (as determined from x-ray standing wave experiments). The lowest unoccupied molecular orbital is occupied on the Ag and Cu surfaces, whereas it remains unoccupied on the Au surface. Simultaneously, the work function increases for Ag but decreases for Cu and Au. Adsorption distances and energies, on the other hand, depend very sensitively on the choice of the DF. While calculations in the local density approximation bind the monolayer consistently with the experimental trends, the generalized gradient approximation in several flavors fails to reproduce realistic distances and energies. Calculations employing the vdW-DF reveal that substantial bonding contributions arise from dispersive interactions. They yield reasonable binding energies but larger binding distances than the experiments.


Introduction
The adsorption of organic molecules on surfaces is a matter of intensive research as many material properties can be decisively modified. For example, in the context of molecular electronics, applying molecular monolayers right at the interface between the electrode and the organic material can be used to tune the electrode work function [1,2] thereby enhancing charge injection [3]. Furthermore, the possibility to change the chemical functionality of the surface enables chemical sensing [4,5], a change in wetting properties or corrosion resistance [6,7].
Much effort has been put into understanding the bond formed between organic molecules and metallic substrates. One of the best characterized molecules of large size is 3,4,9,10perylene-tetracarboxylic acid dianhydride (PTCDA), which has gained fundamental academic importance [8]. In particular, recent photoemission-based experiments carried out for PTCDA adsorbed onto Ag(111), Au(111) and Cu(111) have shown characteristic trends for workfunction change, alignment of molecular levels and vertical binding distances. More specifically, upon monolayer deposition (i) the work function of the metallic surfaces is modified in such a way that initial differences are greatly reduced [9], (ii) the lowest unoccupied molecular orbital (LUMO) of PTCDA becomes occupied for the Ag and Cu surfaces, while it remains empty in the Au case [9] and (iii) the adsorption distance is smallest for Cu, intermediate for Ag and largest for Au [10]- [12]. A coherent explanation for these phenomena is challenging, and support from the theoretical side is required. We performed density functional theory (DFT) calculations for PTCDA adsorbed on the three coinage metals, Ag, Cu and Au and evaluate to what extent the observed trends can be explained theoretically.
Previous investigations based on DFT for PTCDA on Ag(111) have revealed a number of different results showing that a proper description of the present systems is demanding. Studies employing the generalized gradient approximation (GGA) [10], [13]- [15] reported adsorption distances much larger than the ones observed in experiment together with unrealistically small adsorption energies. In [16], the authors found better agreement with experiment, however,  in connection with a strong arching of the molecule. Studies based on the local density approximation (LDA) [17,18] show very high binding energies of about 2.5 eV per molecule and adsorption distances in quite good agreement with experiment. Furthermore, an approach employing exact exchange and the random phase approximation (RPA) yielded reasonable binding energies together with somewhat overestimated binding distances for PTCDA on Ag(111) [19].
Due to these controversial results the present study is divided into several steps. First, the work-function modification and level alignment are investigated by means of the Perdew-Burke-Ernzerhof (PBE) exchange-correlation (xc) functional when placing the monolayer at the experimental distance. A comparison with the experimental results allows us to benchmark these calculations of the electronic structure independent of the evaluation of the bonding distance. Then, the effect of varying the vertical adsorption distance is examined. Here, the focus is laid on the formation of the bond, which gradually builds up as the monolayer approaches the surface. Finally, we explore how the results depend on the choice of xc functionals. In particular, we focus on the LDA and the Perdew-Wang 91 (PW91) functional, and a recently developed van der Waals density functional (vdW-DF) [20,21], which has been designed to include dispersion forces and has been successfully applied to a variety of systems [20]- [24]. The latter allows us to assess how strong the non-local dispersion forces are for the present system and how they impact adsorption energies and distances.

System specifications
In the upper part of figure 1, a top view of a PTCDA monolayer adsorbed on the Ag (left) and Cu (right) surfaces is shown. On the Ag(111) surface, it has been found that PTCDA forms monolayers commensurate with the substrate [25,26]. The blue arrows indicate the 6 1 −3 5 surface unit cell containing two inequivalent PTCDA molecules and 33 atoms per metallic layer. The adsorption site is chosen in agreement with previous investigations [17].
For the Au(111) surface, it has been found that the molecule does not form monolayers commensurate with the substrate but rather exhibits a point-on-line growth on the (22 × √ 3) reconstructed surface [27]- [30]. Inclusion of these effects in the present calculations is hampered by the excessive computational effort arising from the required huge supercells. Hence, we adopt the same unit cell as used in the Ag case as the differences in the lattice constant between Au and Ag are quite small (2.95 Å versus 2.94 Å). This structure can be regarded as a reasonable approximation to the one encountered in experiments. Also on the Cu(111) surface the monolayer has been found to form a herringbone structure [31]. As the lattice parameter, is however, significantly smaller than in the Ag case (2.57 Å versus 2.94 Å) a larger 6 3 −4 6 surface unit cell (figure 1, right) has to be used to accommodate the monolayer.

Methodology
DFT calculations were carried out employing the VASP code [32]. The projector augmented waves (PAW) [33] approach was used allowing for a relatively low kinetic energy cutoff for the plane-wave basis set (about 20 Ry). We use a Monkhorst-Pack [34] 2 × 2 × 1 grid of kpoints in the reciprocal space and a first-order Methfessel-Paxton smearing [35] of 0.2 eV. To avoid spurious electrical fields, a dipole layer is inserted in the vacuum region [36]. As shown in the main part of this work, a proper DFT-based description of adsorption geometries is very challenging for the present system. Hence, several xc functionals were adopted, the LDA as well as the GGA in the PBE [37], PW91 [38] and revised PBE (revPBE) [39] flavor. Calculations with the vdW-DF were carried out following [20,21,40] according to which a non-local contribution to the binding energy is calculated on the basis of the PBE self-consistent charge density. More details about the construction of the functional will be given below.
Partitioning of the density of states (DOS) into contributions from individual molecular orbitals (MODOS) was performed in analogy to previous studies [41,42] using the atomic orbital-based SIESTA code [43]. For these calculations, the PBE xc functional and norm-conserving pseudopotential in the Troullier Martins scheme were adopted including linear core corrections. A double zeta polarization basis and a 200 Ry cutoff value for the mesh, and a 34 meV cutoff value for the localization of the atomic orbitals were employed. Utilizing a different code to perform the partition analysis is possible as the DOS obtained from SIESTA and VASP are very similar.
We use the repeated-slab method to model the metal-organic interface. The twodimensional periodic system is calculated with three-dimensional periodic boundary conditions required by the band-structure program. Hence, the metal is approximated with a finite number of metallic layers and a vacuum gap is inserted between the periodic images to decouple them. We used three metallic layers to perform the structural optimization of the systems in-line with previous investigations [10], [13]- [15], [17]- [19]. In all cases, the size of the unit cell along the surface normal was chosen to be 26 Å, which guarantees a sufficiently large vacuum gap (between 16 and 25 Å) for all studied systems. Note that in section 4.1, we included a fourth layer at the bottom of the metallic slab (subsequent to structural optimization) for a more precise evaluation of work functions and the quantities related to rearrangements in electron density, ρ(z), Q(z) and E vac (z). In all other cases, we report results obtained with the three-layer slab.

Electronic structure at the experimental adsorption distances
When a molecular monolayer adsorbs on a metallic substrate, interactions occur, which cause rearrangements in the electronic structure of both components. These rearrangements can be described either as redistribution in electron density in real space or as a hybridization of molecular and metallic states. The former can be directly connected with adsorptioninduced work-function modifications while the latter reveals changes in orbital alignment and occupation [41]. In this section, we will investigate these effects for the monolayer being placed at the adsorption distances determined from x-ray standing wave (XSW) experiments [10]- [12]. This is achieved by performing a structural optimization of monolayer and slab while applying several constraints. The topmost metallic layer was allowed to relax while the remaining layers where kept fixed to their bulk values. Furthermore, the positions of the carbon atoms along the z-direction (which in this work is normal to the surface) was fixed to their experimental values while the xand y-components (components parallel to the surface) were free to relax. For the oxygen and hydrogen atoms no constraints were applied. Note that in XSW experiments, the distances to the surface are not determined with respect to the actual (relaxed) topmost layer but with respect to the position, which the topmost layer would have in the bulk. For this reason we adopt here the same definition.
The real-space charge redistribution, ρ(z), occurring due to monolayer-substrate interaction is calculated in the following way: where ρ met+mono (z) is the charge density of the combined system while ρ met (z) and ρ mono (z) represent the charge densities of the non-interacting metal and monolayer, respectively. The densities, which are in principle depending on x, y and z are integrated over the x-y-plane within the unit cell and divided by 2 as there are two molecules in the unit cell. ρ(z) is shown in figure 2 (bottom), where a positive value indicates a gain in electron density upon adsorption and a negative peak means a loss of electrons. Also shown in figure 2 (middle part) is Q(z) which is defined as z vac represents a point on the z-axis in the vacuum where ρ(z) has decayed to zero. For each z-position, Q is a measure of electron transfer between the region left to z to the region right to z where a negative (positive) value of Q(z) implies a left-to-right (right-to-left) electron transfer. We first focus on the Ag surface for which charge rearrangements have already been discussed [10,17,44]. These are of strongly oscillating nature, which is a feature typical for metal-organic systems (see e.g. [2,24,41,45]). On the metal side, one can see a pronounced loss of electron density right above the topmost metal layer accompanied by a weak gain on the topmost layer. We attribute this feature primarily to the push-back effect originating from Pauli repulsion [46]. On the molecular side, the plane of the carbon atoms is depleted of electron density, while the regions immediately above and below the plane gain electron density. As we will see later in the analysis of level alignment this feature is consistent with a gain of electrons in the π system counteracted by a loss of electrons from the σ system. In between the topmost metal layer and the carbon plane, Q reaches a pronounced minimum indicating a net metal-to-molecule electron transfer of 0.3 electrons per molecule. In the Cu case, the curves resemble the Ag case, which suggests a similar adsorption mechanism. However, as the molecule is located closer to the substrate, the push-back effect is stronger in comparison to the Ag case. From the minimum in Q, a smaller electron transfer of 0.14 electrons is determined.
In the Au case, a completely different situation is encountered. The rearrangements are much smaller and basically just the push-back feature can be recognized. This gives rise to an electron transfer of about 0.10 electrons from the interfacial region back to the metal. As a next step, we discuss the change in work function induced by monolayer adsorption. In figure 3, we show the potential energy of an electron in the system where PTCDA is adsorbed on Ag. The value of E pot in the vacuum region on the metal side of the slab (left side in the figure) is regarded as the work function of the bare substrate, φ 0 . On the molecular side, the work function of the PTCDA-covered surface, φ, is obtained (right side in the figure). Table 1 lists the work functions for the bare Ag(111), Cu(111) and Au(111) substrates. Au exhibits the highest work function followed by Cu and Ag in good agreement with UPS experiments [9,47]. For the covered surface, an increase in work function is obtained in the Ag case; for the Cu and Au cases (see table 1) the work function decreases with the effect being more pronounced for the latter. Hence, due to monolayer adsorption the work functions of the substrates have become more similar. In the UPS measurements the same trend is observed as the covered surfaces all exhibit identical work functions. The precise numerical value of the work-function change, φ, shows minor variations between theory and experiment (in the Ag case, also between available experimental data sets). However, the qualitative trend, φ most negative for Au followed by Cu and Ag, is observed in both cases (see table 1). In order to understand the origin for these differences we partition the work-function change φ into two contributions, one of purely geometrical nature, E vac , and the other, BD, directly associated with ρ. φ = E vac + BD.
The geometrical component is related to the carboxylic oxygen atoms (see figure 1) that are allowed to freely relax during the structural optimization. They bend down towards the surface (see table 1) indicative of the formation of a bond between the carboxylic oxygen atoms and neighboring metal atoms. This provides the otherwise flat and centrosymmetric molecule with a dipole moment in the direction perpendicular to the surface. Through the Helmholtz equation this gives rise to the potential step E vac . Note that this component reflects the change in work function for the idealized case where ρ = 0 is assumed to be zero and the molecule is already in its distorted geometry. In table 1, we see that E vac reduces the substrate work function in all three cases. The effect is more pronounced for the Ag and Cu cases and weaker for the Au case, which reflects the different reactivity of the surfaces. The bond dipole, BD, is connected to the change in electron density, ρ(z), discussed earlier. When solving the Poisson equation one obtains the change in the potential energy of an electron, E pot (z), induced by ρ(z) [2]. E pot (z) is shown in the top part of figure 2. In the Au case (rightmost plot), where the dominant feature is given by the push-back effect, E pot (z) continuously decreases when going from left to right. The overall decrease is equal to BD and amounts to −0.31 eV. In the Ag case, the  (111) surfaces, as obtained from DFT using PBE or LDA. Note that for the former DF the backbone was constrained to the XSW adsorption distances while for the latter this was not the case. Experimental values are also listed for comparison. They were taken from the following references: (a) [47], (b) [9], (c) [48], (d) [10], (e) [11] and (f) [12]. Quantities are given in eV or Å. decrease due to the push-back effect can also be recognized in the vicinity of the topmost metal layer. However, in the region between metal and molecule, E pot (z) becomes positive, and, therefore, BD is positive. In the Cu case, the situation is intermediate, as the two contributions are very similar and counterbalance each other. Hence, we can clarify the origin for the different work-function changes: while E vac decreases the work function in all cases, BD increases the work function in the Ag case and reduces the work function in the other cases. The dominant effect, determining the general trend, is BD. Next, we compare the alignment of the molecular orbitals with respect to the metal Fermi level. For this purpose projections of the DOS (PDOS) are shown in figure 4 (upper panels) for the three systems. The PDOS originating from the molecule (MDOS) is displayed as a white area. Also shown is the PDOS from the individual molecular orbitals (MODOS). Strictly speaking, as each molecular orbital forms two bands in the PTCDA monolayer (there are two inequivalent molecules in the surface unit cell) the MODOS gives the contributions to the MDOS arising from the two bands. In the lower panels, the peaks as obtained from the UPS spectra [9] are also indicated.  In agreement with previous work [10,17,18], the present calculations show that the LUMO is largely occupied on the Ag surface. The peaks belonging to this orbital are strongly broadened, and the peak maximum is located below the Fermi level. A similar observation is made for the Cu case, where the broadening is even more pronounced. In contrast, in the Au case the LUMO is located right above the Fermi level. The alignment and broadening of the LUMO is in very good agreement with experimental observations [9]. Also UPS locates the LUMO right below the Fermi level on Ag [9,30] and Cu [9], while in the Au case scanning tunneling spectroscopy measurements [49] have shown that the LUMO is located 1.33 eV above the Fermi level. Hence, DFT captures the most striking qualitative feature in the orbital alignment of the three substrates. Also the alignment of the highest occupied molecular orbital (HOMO) is in very satisfactory agreement to the UPS data on Ag and Cu cases; a slight shift upwards can be recognized in the calculations. This is typical for DFT in the GGA and is usually accounted for by aligning the calculated spectra on the first occupied peak and expanding the energy scale by a factor of 1.2-1.3 [41,50]. In the present case, such an alignment is not necessary as the LUMO is right at the Fermi level in both cases. On Au this overestimation is more pronounced; here the HOMO position is located about 1 eV too high in energy which has, however, a different origin: as the monolayer is weakly interacting with the substrate, the calculation correctly locates the LUMO above the Fermi level. Hence, the position of the HOMO cannot be located at lower energies as is allowed from the size of the DFT band gap. In the weakly interacting situation, DFT calculations underestimate the band gap [51]- [53] and, hence, locate the HOMO at higher energies.
The MODOS also allows for an orbital-resolved electron transfer analysis. For that purpose, the MODOS is integrated from minus infinity up to the Fermi level, which gives the occupation of each individual molecular orbital. The resulting orbital occupation diagram for PTCDA on Ag(111) is shown in figure 5. The LUMO is occupied by 67%, which corresponds to an electron transfer of about 1.4 electrons into the LUMO. As already pointed out (for the HOMO and HOMO-1) [48] this gain of charge in the LUMO is counteracted with a loss of charge from many formerly occupied orbitals, which do not quite reach the 100% level. The overall electron transfer, obtained from adding the occupation of all orbitals, is, therefore, considerably smaller, namely 0.47 e. Note that, due to the way the projections are carried out [42], this value corresponds exactly to the Mulliken charge.
This delicate balance between forward and backward donation closely resembles the one observed for F4-TCNQ on Cu(111) [41] and Ag(111) [54]. There, however, the electron transfer to the metal predominantly originates from four molecular orbitals, which are constructed from the lone pairs of the four cyano docking groups. The situation for PTCDA is more complex.
Here, the lone pair orbitals from the carboxylic oxygens contribute to about 20 orbitals. Hence, the loss of electrons is portioned over many molecular orbitals and does not give rise to a clear feature in the occupation diagram. For PTCDA on Cu, the occupation diagram is very similar and shows a LUMO occupied by 1.50 electrons and an overall electron transfer of 0.50 electrons. Hence, in contrast to the transfer obtained from Q (which was smaller for Cu compared to Ag), the Mulliken electron transfer indicates an equally strong electron transfer for Cu. 4 In the Au case, the orbital occupation analysis is qualitatively different and rather straightforward as all formerly unoccupied (occupied) orbitals are still almost fully unoccupied (occupied) on the surface. Here the overall electron transfer is with 0.05 electrons basically negligible in good agreement with the observation that the charge neutrality level is located in very close proximity to the Fermi level [51,52].

Metal-monolayer separation
After having established that the DFT electronic structure as calculated with PBE compares quite well with experiments when placing the monolayer at the measured distances, we investigate how quantities vary with the metal-monolayer separation distance. For that purpose, geometric relaxations as explained above were performed, now fixing the carbon atoms of the monolayer at distances d also different from experiment. We first focus on the position of the carboxylic oxygen atoms with respect to the carbon backbone ( figure 6(a)). The oxygen atoms are located below the carbon backbone in all cases and the bending of the oxygen atoms towards the surface builds up gradually when the molecule approaches the surface. On Ag and Cu this effect is strongest at the experimental positions and is reduced at 0.5 Å below, where also the oxygens start feeling repulsion. On Au, the effect is less pronounced at equivalent distances, which shows that the reduced downward-bending observed previously is not only a consequence of the higher adsorption distance. In XSW experiments, the downward bending has been observed on Ag [10]. The magnitude of the bending compares reasonably with the values found in the present calculations (table 1). On Au, the oxygen positions have not been determined experimentally, and in the Cu case the oxygen atoms were measured to be located further away from the surface than the molecular backbone (see table 1). This finding is not confirmed by the present calculations. Such a distortion does not occur, even for very small distances.
Also, the work-function modification displays a significant dependence on distance as shown in figure 6(b) 5 . On Au, φ is negative at all separations and gradually increases in magnitude as the monolayer approaches the surface. On Ag, φ is seen to be almost constant at larger distances. The decrease of the work function arising from the push back is counteracted by the charge transfer which, with the exception of the smallest considered d, results in a weak work-function increase. For the Cu substrate, the situation is intermediate between the Ag and Au cases as the plateau at larger distances can also be recognized here, but is not as pronounced. On Ag, the work function changes very significantly between 2.5 and 3.0 Å. Hence, a correct prediction of the increase or decrease requires a very precise description of the adsorption distance.
The distance dependence of the alignment of the LUMO with respect to the Fermi level, as derived from the MDOS, is shown in figure 6(c). It was determined from the centers of  the corresponding peaks in the MDOS. The position of the LUMO for the hypothetical noninteracting case, i.e. for the situation where the metal and molecule electronic states are not interacting but aligned to a common vacuum level is indicated by the horizontal lines. Hence, the graph also shows how the orbital shifts upwards or downwards due to the interactions. In general, a correlation between orbital shifts and work-function change can be recognized as in all cases the orbital moves downwards with decreasing distance. Also here, on Ag and Cu a plateau at larger distances can be seen. This behavior is not just observed for the LUMO but also for the HOMO ( figure 6(c)). It is attributed to the fact that the potential energy landscape at the position of the molecule is modified in accordance with the work-function change.
The downward shift of the HOMO is, however, not as strong as for the LUMO. As shown in figure 6(d) the molecular band-gap, E gap , reduces upon approaching the surface. This is reminiscent of the model of Newns-Anderson, which was discussed for the present system in [8]. According to this model a reduction of the band gap arises from dynamical image charge effects. However, this is an effect generally not captured by DFT [55]. Moreover, a proper description of the situation is hampered by the underestimation of DFT band gaps. In fact, GW calculations for the isolated PTCDA molecule have reported a band gap of 4.9 eV [53] and experiments for PTCDA multilayers have revealed about 4 eV [49]. Therefore, having the larger band gap as a starting point, the reduction upon adsorption should be much more pronounced. The small decrease, which is nevertheless observed in the present calculations, is attributed to the formation of a bond between the carboxylic oxygen and the metal surfaces, which renders the molecular backbone more quinoidal [17]. This destabilizes the HOMO and stabilizes the LUMO and, thus, reduces the gap.
A surprising behavior is observed on Ag(111) at large distances: rather than approaching zero, the change in work function converges to a constant value of +0.2 eV. This means that even at very large distances, there remains a finite charge transfer between molecule and metal, at least when assuming thermodynamic equilibrium. This is rooted in the fact that vacuum level alignment places the LUMO below the Fermi level in the Ag case (see figure 6(c), the vertical gray solid line lies below zero) and, therefore, the LUMO should be occupied, no matter how far the molecule is away from the surface. This full charge transfer can, however, not occur as it would result in a gigantic BD (at a separation of 10 Å BD would be 27 eV). Therefore, in the equilibrium situation a much smaller amount of charge is transferred, which shifts the molecular levels by the necessary amount such that the LUMO is unoccupied. GW calculations or calculations with hybrid functionals, that are known to yield larger band gaps [53] would not show this 0.2 V work-function change at large separations. We expect the situation to be similar to the Au and Cu cases where the LUMO is located above the Fermi level and, consequently, no change in work function and level alignment would be observed at large distances. Hence, the level alignment and work function on Ag at large separations are affected by artificial contributions, which result from the underestimation of the DFT band gap. Fortunately, as we will see below the error in the adsorption potential due to this artifact is rather small.
The adsorption potential, E ads (d), is calculated in the following way: where E met+mono (d) is the energy per unit cell of the metal-monolayer system. As explained in the previous sections, the monolayer is placed at specific distances d to the surface. E met is the energy of the relaxed isolated metallic slab, i.e. where the monolayer has been removed from the unit cell. E mono is the energy of the relaxed isolated monolayer. To estimate the error 14 due to the above-discussed artificial dipole on Ag, we calculated a system with the monolayer placed in the middle of the vacuum gap. This corresponds to a highly separated situation with d = 10.6 Å. The actual value of 0.012 eV is vanishingly small and will, therefore, not be considered in the following. As shown in figure 6(e), E ads is found to be largely repulsive for all three metals, where a very small binding contribution can be observed at large distances to the surface (about 4 Å). For the Ag substrate this has already been described in [14,15,19]. The E ads curves strongly contrast the XSW results, which show that the monolayer should be adsorbed much closer to the surface for all three metals (top part of figure 6(e)). Also the minimum in the bonding potential can be expected to be much more pronounced. Experimentally, the desorption energy of PTCDA is not accessible as the molecule dissociates prior to desorption (see e.g. [26]). However, for NTCDA, which has a naphthalene instead of perylene building block in the molecular center, thermal desorption spectroscopy at submonolayer coverage has reported desorption energies of 1.2 eV [56]. As NTCDA is about half the size of PTCDA, one can expect that the adsorption energies should be as large or even larger. Hence, although GGA calculations can reproduce the experimental trends related to the charge density when fixing the molecular backbone to the experimental positions, the adsorption potentials are highly inaccurate.

Impact of xc-functional
In this section, we will investigate what results are obtained using functionals other than PBE. We first focus on the PW91 xc-functional, another GGA parametrization, which has been used in [16] where a negative adsorption potential of −0.54 eV was reported together with a strong arching of the molecule. These qualitative differences in comparison to the PBE findings have been attributed to the different nature of the functional [17]. Furthermore, we also show results obtained with LDA as it has been found to bind PTCDA on Ag [17,18] at very reasonable binding distances.
As far as the DOS of the monolayer and the charge rearrangements for the three functionals are concerned, there is virtually no difference between the functionals as long as the same geometry is adopted. In figure 7, it is shown that the curves are almost on top of each other; just in the LDA case a closer look revels that the LUMO is slightly shifted towards higher energies and consequently the charge rearrangements on the molecule are also somewhat smaller. Also the adsorption potential as obtained with PW91 is very close to the PBE results ( figure 8). Hence, the differences between [16] (which reported significant binding using PW91) and [14,15] (which did not find binding using PBE) are not a consequence of the different GGA flavors but must have a different origin. We performed a structural relaxation of PTCDA on Ag employing PW91, starting from the molecule adsorbed flat on the surface (using four metallic layers) at the experimental position. We also found that the backbone is repelled by the metal while the oxygen positions remain constant. When using 0.02 eV Å −1 for the maximum final force as the criterion for the termination of the structural optimization (as in [16]) the reported arching could be reproduced. However, when reducing the convergence threshold to 0.01 eV Å −1 the molecule continues to move away from the surface. In fact, a geometry where the molecule is placed flat on the surface at the larger distance of 4.3 Å gives a lower value in the adsorption potential (by 0.08 eV) than the arched configuration. The reason for the discrepancy concerning the bonding potential, is attributed to the fact, that in [16] the molecule and not the monolayer is taken as the reference point in equation (5) [57]. As monolayer formation already lowers the energy by 0.5 eV per molecule the negative adsorption potential reported in [16] arises from this contribution. In conclusion, the PW91 functional gives results very similar to PBE, as expected. This is not the case for the LDA. In agreement with earlier works on Ag [17,18], a structural optimization of the monolayer where no constraints are imposed on the molecule, gives binding distances in reasonable agreement with XSW experiments and realistic negative values for the adsorption potential (−3.04 −3.20 and −1.80 eV for the Ag, Cu and Au surfaces, respectively). Although the calculations place the carbon backbone somewhat closer to the surface (see d in table 1) they correctly reproduce the relative ordering, d being smallest for Cu and largest for Au. The same holds for the level alignment as the LUMO is located at −0.49, −0.61 and 0.14 eV relative to the Fermi level for the Ag, Cu and Au cases. For the work functions, however, the experimental trends are not captured. A reduction is obtained in all cases and is strongest on Cu. The reasons for this behavior are rooted in the unconstrained structural optimization which, also here, leads to arching (see d carb O , d anhy O and d in table 1). As a result, E vac , which reduces the work function, is too pronounced 6 .
The calculations presented so far indicate that the following distinct mechanisms are operative when PTCDA adsorbs on the surfaces: (1) A bond forms between the carboxylic oxygen and neighboring metal atoms, which is reflected in the downward bend of the carboxylic oxygen atoms. This leads to arching of the molecule and a reduction of the work function. The effect was seen to be more pronounced in the Ag and Cu cases as compared to the Au case. Simultaneously, this bond formation leads to a change in bond length alternation of the whole carbon backbone driving it towards a more quinoidal structure [17] thereby decreasing the band gap. This effect can be expected to be energetically unfavorable for the molecule and, therefore, the backbone has a tendency to inhibit the formation of the oxygen-metal bonds. This is in contrast to the situation of other molecules, such as F4-TCNQ [41] or thiol-based selfassembled monolayers [2], for which GGA already gives sizable binding energies and the backbones become more aromatic due to the bonding through the docking group. (2) Pauli repulsion, also known as the pillow effect or push-back effect, arises from the exchange interaction and is a repelling contribution leading to a reduction of the work function when the monolayer approaches the surface. (3) Metal-to-molecule electron transfer arises from the downward shift of the molecular levels, which in the Ag and Cu cases leads to an occupation of the LUMO. In the Au case, where the molecule is located at higher distances to the substrate, this downshift is not as pronounced and the LUMO remains unoccupied. Based on electrostatic arguments, this charge transfer is expected to be attractive and leads to an increase in work function.
Neither of the functionals discussed so far includes dispersive interactions, which are operative between the polarizable perylene backbone and the metal. For systems which are bound just by dispersive interactions (such as rare gas or benzene dimers [20,22] or graphene sheets [21]), PBE and LDA calculations are known to give rather similar bonding potentials as in the present cases. The former gives small binding energies at too high adsorption distances while the latter strongly binds these systems. This comparison suggests that, as also found for thiophene adsorbed onto Cu(110) [24], for PTCDA on metals, the dominant bonding effect arises from the dispersion interactions.

Dispersion forces
As a last step, we will, therefore, investigate what is found when employing a recently developed functional, which has been designed to include vdW interactions from first principles. The binding energy, E nl c , from non-local correlations, predominantly arising from dispersion forces, is calculated in the following way [20,21]: where the kernel K ( r , r ) has been derived based on the adiabatic connection fluctuation dissipation theorem, the RPA, and using a plasmon pole model for the polarizability. The vdW-DF, E vdW xc [n], is, then, given by [20,21]: where all terms are calculated based on the self-consistent PBE charge density n. For the exchange part, the revised version of PBE [39] is chosen as this functional, in comparison to PBE, was found to resemble exact exchange calculations more closely [20]. Figure 9 shows the adsorption potential when using the revPBE xc-functional (top left, lines and symbols) and when using E revPBE [n] (top left, bold lines). The curves are strongly repulsive. In the top right part of figure 9, the adsorption potential obtained from E nl c alone is displayed. The curve is purely attractive for the three different substrates and rather similar for the three surfaces, albeit Au gives a slightly stronger attraction than Cu and Ag. Note that dispersion forces yield a very strong binding contribution of about 5-6 eV per molecule at the experimental distances.
The adsorption potential as obtained with E vdW xc is shown at the bottom of figure 9. It exhibits a pronounced minimum of about −2 eV, which is found for all three substrates at around 3.5 Å. In the Ag case, the result can be compared to [19] where E nl c as calculated within the RPA is combined with an exact exchange approach [19]. In this way, a similar adsorption potential was obtained giving a somewhat smaller binding distance (3.2 Å) and a comparable adsorption energy of 2.5 eV. Furthermore, a study based on an embedded atom model reported an adsorption energy of about 3.5 eV for PTCDA on Ag(111) [58].
Compared to the experimental situation the bonding distances are systematically overestimated (see figure 9, bottom) with the biggest deviations obtained for the Cu surface. Hence the vdW-DF corrects the failure of the GGA functionals insofar as it results in reasonable binding energies but it still shows a tendency to bind the monolayers at too high distances. Several reasons for these deviations could be listed. First, energies are calculated on the basis of the self-consistent PBE charge density, while a fully rigorous treatment would require a self-consistent calculation with the vdW-DF. As argued in [22] this can, however, be expected to yield only a small correction. Furthermore, the plasmon pole model used to describe the polarizability might be too approximative to capture the attraction from dispersion forces correctly. However, the good agreement with [19] suggests that the so-induced error is most likely not significant. The deviation could also be due to specific choices in the construction of the vdW-DF (equation (7)): revPBE is used instead of PBE, which was motivated by the fact that the former is closer to exact exchange calculations. Furthermore, the GGA part of the correlation is abandoned to avoid double counting. The latter could be a problem especially for strongly interacting electronic systems.

Conclusions
We investigated the adsorption of PTCDA on the (111) surface of the three coinage metals, Ag, Cu and Au, by means of DFT calculations. The focus was laid on understanding substratespecific trends in work-function modification, level alignment and adsorption distance. Several xc-correlation functionals, LDA, GGA in the PBE and PW91 flavor and the vdW-DF were used for that purpose.
We found that, when fixing the monolayer to the distances measured by XSW, the workfunction modification and level alignment as calculated with PBE are consistent with the results obtained from UPS measurements. The work function is reduced in the Au case as here the molecule is physisorbed, the LUMO remains unoccupied and only Pauli push-back is encountered. In the Ag and Cu cases, the LUMO is strongly occupied while several formerly occupied orbitals are slightly depleted. This leads to an overall metal-to-molecule electron transfer of about 0.5 electrons, which counteracts the push-back effect. In the Ag cases, this leads to an increase in work function while on Cu, where the push-back effect is stronger, the work function is decreased. This observation does not depend on the choice of DF. This shows that, with a precise knowledge of the adsorption geometry, a reliable prediction for workfunction changes and level alignments can be achieved with all functionals.
The calculated adsorption distances and binding energies, however, reveal a strong sensitivity to the xc-functional. GGA, which is not capable of describing dispersion forces leads to clearly wrong adsorption energies, while LDA, which suffers from the same short-coming, gives results rather close to the XSW data due to a fortuitous cancelation of errors. The attraction arising from dispersion forces amounts to several eV for all three substrates. Including this nonlocal effect, the GGA results can be improved giving adsorption energies of about 2 eV per molecule together with adsorption distances of about 3.5 Å for all substrates. The latter are not fully consistent with the XSW results revealing that further developments from the theoretical side are still required.