Theoretical study of the adsorption of benzene on coinage metals

Summary The adsorption of benzene on the M(111), M(100) and M(110) surfaces of the coinage metals copper (M = Cu), silver (M = Ag) and gold (M = Au) is studied on the basis of density functional theory (DFT) calculations with an empirical dispersion correction (D3). Variants of the Perdew–Burke–Ernzerhof functionals (PBE, RPBE and RevPBE) in combination with different versions of the dispersion correction (D3 and D3(BJ)) are compared. PBE-D3, PBE-D3(BJ) and RPBE-D3 give similar results which exhibit a good agreement with experimental data. RevPBE-D3 and RevPBE-D3(BJ) tend to overestimate adsorption energies. The inclusion of three-center terms (PBE-D3(ABC)) leads to a slightly better agreement with the experiment in most cases. Vertical adsorbate–substrate distances are calculated and compared to previous theoretical results. The observed trends for the surfaces and metals are consistent with the calculated adsorption energies.


Introduction
The adsorption of organic molecules on metals is of great interest since the formation of thin films and self-assembled monolayers opens the way toward a functionalization of surfaces [1][2][3][4][5][6][7][8]. The adsorbed molecules often contain an aromatic framework that can be substituted with functional groups. The bonding between the surface and the adsorbate is an interplay between electrostatic interaction, including charge transfer (CT) to the surface, and covalent contributions [9][10][11]. In addition, it was found that dispersion interaction plays a crucial role for the adsorption of large aromatic compounds on metal surfaces [9][10][11]. This holds in particular for the adsorption on the coinage metals copper, silver and gold. Therefore, a theoretical treatment of this process requires methods that provide an accurate description of these weak interactions. Density functional theory (DFT) is established as a standard method for quantum-chemical solid-state calculations [12]. However, DFT has the well known shortcoming that it fails to describe dispersion effects. Consequently, standard DFT methods are not suitable for the calculation of the adsorption of aromatic compounds. In the last years much effort has been directed to the development of DFT methods that eliminate this shortage [13][14][15][16][17][18][19][20][21][22][23][24]. One of them is a damped empirical correction called DFT-D3 which was proposed by Grimme et al. for molecular systems [13]. The D3-dispersion correction to the DFT energy is calculated by summation over pair potentials. Nonadditive effects of dispersion interaction can be treated on the basis of three-body terms D3(ABC) [13]. The most recent DFT-D3(BJ) method [16] differs from the original DFT-D3 essentially only in the damping function for short range interaction. Due to the computational efficiency of the D3 correction schemes it is possible to perform DFT-D3 calculations with nearly the same computational effort as standard DFT calculations. Only the calculation of three-body terms can become expensive for large systems, and is therefore usually carried out only for the final energy estimation. For molecular systems one can obtain with DFT-D3 results that are close to coupled cluster singles doubles with perturbative triples (CCSD(T)) results at the cost of GGA-DFT calculations. In recent years the DFT-D3 method has been extended to periodic systems [25]. In a few cases [26,27] it is observed that the metal C 6 dispersion coefficients for bulk systems can be largely reduced compared to the values of free atoms. Indeed this does not apply to the coinage metals Cu, Ag, and Au. The C 6 parameters for these metals are already converged [28]. In addition, it is known that metal substrates show significant dispersion screening effects that can modify the polarizabilities and C 6 coefficients of adsorbed molecules [29,30]. In principle these effects should be included in the coordination number dependent C 6 coefficients of the D3 correction. We checked this by calculating the C 3 coefficients for the benzene adsorption on the Au(111) surface.
Previous theoretical studies of the adsorption of organic compounds on silver and gold surfaces resulted in a good agreement with experimental results [9][10][11]26]. However, a systematic comparison of the different DFT-D3 approaches is still missing. One aim of this work is the comparison of different DFT-D3 methods for the description of the adsorption of aromatic compounds on this surfaces. We therefore present a theoretical study of the adsorption of benzene on the M(111), M(100) and M(110) surfaces of the coinage metals copper, silver and gold. The benzene molecule was selected because it is a building block of many organic compounds that are used for surface functionalization. As mentioned above the binding between the metal and the adsorbate is often dominated by dispersion interaction to the aromatic framework. Therefore the study of the bonding between benzene and the metal surfaces is of great interest. In addition, although the adsorption of benzene on some of these surfaces has been subject of previous theoretical studies [29][30][31][32][33][34][35][36][37][38][39], this is the first work in which the adsorp-tion on the most important coinage metal surfaces is systematically studied with the same method. Different from our previous study on benzene/Ag(111) [25], we apply a variety of DFT methods and dispersion corrections, and investigate all lowindex surfaces.

Computational methodologies
We used the plane-wave code VASP [40][41][42] in combination with the projector-augmented wave method to account for the core electrons [43] for all calculations. We applied our recent implementation [25] of Grimme's dispersion correction (DFT-D3) [13,16]. The dispersion corrected DFT-D3 energy E DFT−D is calculated by adding an empirical correction energy E disp to the DFT energy E DFT , see Equation 1. (1) In this work we used the gradient-corrected PBE [44], RPBE [45], and RevPBE [46] functionals in combination with the original D3 [13] as well as with the newer D3(BJ) dispersion correction [16]. These methods were chosen since they represent a selection of standard GGA functionals, which are available in most of the software for quantum chemical solid state studies. The computational effort of hybrid functionals is too large for these systems and cheaper methods like DFTB-D3 [47] lack suitable parameters, e.g., for gold. We also checked the impact of the three-body terms (D3(ABC)). These threecenter terms D3(ABC) were introduced in the D3 correction scheme since the long-range part of the interaction between three ground-state atoms is not exactly equal to the pairwise interaction energies [13]. From third-order perturbation theory one gets the Axilrod-Teller-Muto dispersion term for the consideration of the non-additivity of dispersion interaction. D3(ABC) is calculated according to Equation 2 (2) The coefficients are the geometric mean of the C 6 -coefficients, θ a , θ b and θ c are the angles of the triangle, which is formed by the three atoms, and f damp is the damping function. The Axilrod-Teller-Muto dispersion terms can be neglected for molecular systems [13]. However, recent studies indicate that their impact is larger for periodic systems that are more densely packed.
We performed calculations with the PBE-D3, PBE-D3(BJ), PBE-D3(ABC), RPBE-D3, RevPBE-D3, and RevPBE-D3(BJ) approaches. We used a cutoff energy of 400 eV for the planewave valence basis and a 3 × 3 × 1 k-point mesh for reciprocal-space integration. Preliminary convergence studies showed that these values are an optimal compromise between accuracy and computational efficiency. The calculated structure parameters and adsorption energies changed by several mÅ and kJ/mol, respectively, when larger cutoff energies 600 eV and 900 eV were used. However, the reported trends did not change. Geometry optimizations were performed with the PBE-D3 method, which was already used for the study of PTCDA on the Ag(111), Ag(100) and Ag(110) surfaces [9]. For these calculations we chose an energy convergence criterion of 10 −6 eV for the SCF energy, and of 5·10 −3 eV/Å for the ionic relaxation (forces are converged if smaller than 5·10 −3 eV/Å). The first three atomic layers of M(100) and M(111) and the first five atomic layers of M(110) were relaxed while the atoms of the lowermost layers were kept at their bulk-like positions. We give two different values for the adsorption distance. d 1 is the vertical distance between the adsorbate and the topmost layer of the surface. The adsorption of benzene effects a slight relaxation of the surface. Therefore we calculate d 1 with respect to the mean value of the z-coordinates. d 2 is the distance between the adsorbate and the hypothetical topmost surface layer for an unrelaxed surface. The latter is given since this enables a comparison to data which have been derived from normal incidence X-ray standing waves (NIX-SW) spectroscopy. This may be useful for a comparison with future experimental results, similar to our previous studies [9][10][11]. Both distances are determined by calculating the difference between the averaged z-coordinate of the carbon atoms and the surface atoms.
Potential curves are obtained on the basis of single-point calculations with the PBE-D3, PBE-D3(BJ), PBE-D3(ABC), RPBE-D3, RevPBE-D3(BJ) and RPBE-D3 approaches. We stepwise altered the distance between the benzene molecule and the unrelaxed surface from 2.5 Å to 5.0 Å. Since d 1 and d 2 do not differ for the potential curve we give only one adsorption distance d. Adsorption energies E ads are calculated according to the supramolecular approach: [48,49] was performed in order to study the net charge transfer between the surfaces and the adsorbed molecule.

Results and Discussion
First we identified the most stable adsorption sites on the basis of PBE-D3 optimizations by placing the benzene molecule on different adsorption sites and calculating the adsorption energies. The following discussion is limited to the most stable adsorption sites for each surface. Other adsorption sites are higher in energy by 2 to 15 kJ/mol. This points to a relatively high mobility of benzene on the surfaces assuming that the corresponding activation barriers are of similar magnitude. We did not calculate the transition states because this was not the scope of the present study. We also checked that the orientation of the flat lying molecules with respect to the underlying surfaces has no significant influence to the adsorption energies and geometries: A stepwise rotation of the molecule around the z-axis changes E ads by less than 1 kJ/mol in each case.
The results of the optimizations are summarized in Table 1. It was found that all metals have the same preferred adsorption sites for the M(111) and M(100) surfaces, the threefold hollow site for M(111) and the fourfold hollow site for M(100), see Figure 1. As a general trend we observe that structures with shorter C-metal distances are more stable than the others. For the M(110) surface it was found that Cu prefers benzene adsorption on a fourfold hollow position, whereas Ag and Au prefer a bridge position, see Figure 1. The notation of the adsorption places refers to the center of the benzene molecule.
The binding between benzene and copper is the largest of all investigated surfaces. E ads is calculated to E ads = −117 kJ/mol for Cu(110), E ads = −114 kJ/mol for Cu(100) and E ads = −97 kJ/mol for Cu(111). For gold the adsorption energies are E ads = −85 kJ/mol for the Au(110), E ads = −87 kJ/mol for the Au(100), and E ads = −85 for the Au(111) surface. The lowest adsorption energies are calculated for the silver surfaces, E ads = −76 kJ/mol for Ag(110), E ads = −75 kJ/mol for Ag(100) and E ads = −72 kJ/mol for Ag(111). It is worth to mention that the adsorption energies on a given metal are almost independent from the surface type. The sole exception is the Cu(110) surface the adsorption energy of which is about −20 kJ/mol smaller than E ads on Cu(100) and Cu(110). For all other systems we find that the adsorption energies on the selected surface planes are within the range of 3 kJ/mol. The results are summarized in Figure 2.     (110). This trend is consistent with the calculated adsorption energies, which are smallest for silver and largest for copper. Unfortunately, to the best of our knowledge no experimentally determined distances are available for these systems. However, previous experimental and theoretical studies of PTCDA on the Ag(111), Ag(100) and Ag(110) surfaces [9,10] indicate that the PBE-D3 and PBE-D3(BJ) approaches give accurate adsorption distances.
A closer look at the optimized structures reveals that the benzene molecule and the underlying surfaces are only slightly affected by adsorption. Therefore, the following comparison of different DFT-D approaches has been performed on the basis of potential curves with fixed structures of benzene and surface. The results of these calculations are summarized in Table 2. The potential curves for the silver surfaces are shown in Figure 3. The potential curves for the copper and gold surfaces are included in Supporting Information File 1.
In all cases the vertical distances obtained with the potential curves for PBE-D3 are in good agreement with the results of the full geometry optimization. This confirms the validity of the simplified approach.
It was found, that PBE-D3 and PBE-D3(BJ) give similar results for all surfaces. The PBE-D3(BJ) adsorption energies tend to be about 4 kJ/mol larger in absolute value than the PBE-D3 energies. This behavior is to be expected: Both methods essentially only differ in the empirical damping function for short-range  Table 2 we give the results of PBE-D3(ABC) calculations. The three-body correction to dispersion is repulsive in this case. This is in line with a previous study of the influence of the three-body terms to periodic systems. The PBE-D3(ABC) adsorption energies are 11 to 20 kJ/mol less negative than the PBE-D3 energies. As expected it was found that none of the pure DFT functionals is able to give a correct description of the adsorption. Most potential curves (not shown) are repulsive over the entire distance range. Only the potential curves obtained with the PBE functional exhibit some very flat minima in the range of −10 kJ/mol at larger distances. Accordingly the calculated adsorption energies can be almost solely ascribed to the dispersion correction. The contribution of E disp to the adsorption energy is larger than 90% for all systems. This confirms that the surface-adsorbate interaction is dominated by dispersion interaction.
Nevertheless we also calculated the charge transfer between benzene and the metal surfaces on the basis of a Bader analysis in order to investigate electrostatic contributions to the interactions. For the Cu(111), Ag(111) and all gold surfaces we observe only a small charge transfer, between 0.02 a.u. (for Ag(111)) and 0.12 a.u. (for Au(100)) from benzene to the metal surface. For the other systems we calculate a small charge transfer from the surface to the benzene molecule in the range of 0.03 (for Ag(100)) to 0.07 a.u. (for Cu(110)). As expected this charge transfer is much smaller than the one between functionalized aromatic compounds and coinage metal surfaces [9,11]. The different direction of the charge transfer may be In Table 3  Therefore we conclude that RevPBE-D3 and RevPBE-D3(BJ) methods are not suited for the calculation of aromatic organic compounds on these metal surfaces. However, the application of the other DFT-D methods treated in this work can be recommended.
The differences in adsorption distances for the recommended methods are in the range from 0.03 to 0.18 Å. As expected, PBE-D3(ABC) gives the largest distances for all systems due to the repulsive nature of the three-center terms. However, the deviation from the PBE-D3 distance is less than 0.07 Å. Therefore it is possible to neglect these contributions in structure optimizations without significant loss of accuracy, which is advantageous since the calculation of the three-center terms is rather expensive for large systems.    [29] In Table 4 we give the C 3 coefficients for the benzene adsorption on the Au(111) surface. The data are compared to the C 3 coefficient for the PBE+vdW surf functional, which is constructed to reproduce the exact values [29]. The C 3 coefficients are fitted from the pure dispersion interaction term at large distances according to the method described in [29].
It was found that PBE-D3, PBE-D3(BJ) RPBE-D3, RevPBE-D3 and RevPBE-D3(BJ) overestimate (by 30%) the C 3 coefficients compared to PBE+vdW surf whereas PBE-D3(ABC) underestimates them (by 50%). It has to be mentioned, however, that the values of the C 3 coefficients strongly depend on the functional. For example, in [29] values between 4 and 9 eV·Å 3 have been reported. Accordingly, we think that our deviations are within a reliable range.

Conclusion
The The adsorption distances are almost the same on these two surfaces due to the similar van der Waals radii of 172 pm for Ag [55] and 166 pm for Au [55]. Hence, the larger C 6 coefficients of 317.2 a.u. for the gold atoms [56] (compared to 268.6 a.u. for the silver atoms [56]) result in a larger dispersion interaction between the surface and the substrate. Copper has the smallest polarizability of the three metals. The C 6 coefficients for the surface atoms is 175.0 a.u. within the D3-correction [56]. However, the benzene molecules come closer to the copper surfaces due to the smaller van der Waals radius of 140 pm [55]. From there the dispersion interaction to the copper surface is largest although Cu has the smallest C 6 coefficients. PBE-D3, PBE-D3(BJ) and RPBE-D3 tend to slightly overestimate the adsorption energies in comparison to experiment, in particular for the Cu(111) surface. This effect is reduced when the three-body correction to dispersion is considered. The PBE-D3(ABC) adsorption energies are smaller in absolute value by 10 to 20 kJ/mol compared to the standard PBE-D3 values. This leads in most cases to a slightly better agreement with the available experimental results. As a result of this work we recommend DFT-D methods like PBE-D3, PBE-D3(BJ) or RPBE-D3 for the theoretical study of the adsorption of aromatic compounds on metal surfaces. Due to the high computational cost of the evaluation of three-center terms, we suggest to perform geometry optimizations with PBE-D3 followed by single-point calculations with PBE-D3(ABC) for adsorption energies. Surprisingly, we realize that the RevPBE-D3 and RevPBE-D3(BJ) methods, which yield a more realistic description of the adsorption of small molecules on ionic surfaces [26], seem to be not suitable for the present systems. Since the pure RevPBE potential curves are repulsive one can ascribe the observed overestimation of the adsorption energies solely to the dispersion correction.
We give adsorption distances with respect to the topmost layer of the relaxed surface, and with respect to the hypothetical topmost surface layer for an unrelaxed surface. The latter allows for a comparison to data from NIX-SW spectroscopy. It was found, that vertical distances are smallest for the M(110) surfaces and largest for the M(111) surfaces, in accordance with the trends of the adsorption energies. The distances on the copper surfaces are in the range from 2.35 to 2.86 Å, the distances on the silver surfaces in the range from 2.78 to 3.17 Å, and the distances on the gold surfaces are in the range from 2.84 to 3.18 Å. It will be interesting to compare these data with future experimental data.

Supporting Information File 1
Potential curves for adsorption of benzene on copper and gold surfaces.