Density functional theory with modified dispersion correction for metals applied to self-assembled monolayers of thiols on Au(111)

Using sound physical principles we modify the DFT-D2 atom pairwise semiempirical dispersion correction to density functional theory to work for metallic systems and in particular self-assembled monolayers of thiols on gold surfaces. We test our approximation for two functionals PBE-D and revPBE-D for lattice parameters and cohesive energies for Ni, Pd, Pt, Cu, Ag, and Au, adsorption energies of CO on (111) surfaces of Pd, Pt, Cu, Ag, and Au, and adsorption energy of benzene on Ag(111) and Au(111). Agreement with experimental data is substantially improved. We apply the method to self-assembled monolayers of alkanethiols on Au(111) and find reasonable agreement for PBE-D and revPBE-D for both physisorption of n-alkanethiols as well as dissociative chemisorption of dimethyl disulfide as an Au-adatom-dithiolate complex. By modifying the C 6 coefficient for Au, we obtain quantitative agreement for physisorption and chemisorption for both PBE-D and revPBE-D using the same set of parameters. Our results confirm that inclusion of dispersion forces is crucial for any quantitative analysis of the thiol and thiolate bonds to the gold surface using quantum chemical calculations.


Introduction
Density functional theory (DFT) is the method of choice for first-principles calculations in condensed matter systems and has contributed greatly to our understanding of metallic systems such as heterogeneous catalysis of ammonia synthesis [1].
However, conventional DFT functionals do not take into account van der Waals interactions, that is, London dispersion. These interactions are crucial for many systems such as interlayer bonding in graphite [2] and biological systems. Research in the last decade has led to dispersion being included in DFT in many ways [3]. Some methods that have the correct asymptotic 1/ 6 behavior are nonlocal dispersiondensity functional (vdW-DF) [4,5] and semiempirical atom pairwise dispersion [6][7][8][9]. Some highly parameterized meta-GGA functionals also include short-range dispersion effects, like the M06 family of functionals [10], but do not have the correct long-range asymptotic behavior [3]. The vdW-DF functional takes into account electronic effects such as electron transfer from first-principles. Its accuracy for normal thermochemistry is however not well established yet. It is furthermore not well defined for spin-polarized systems, such as Fe, Ni, Co, and their alloys.
Dispersion effects included via semiempirical atom pairwise interactions using the DFT-D2 or DFT-D3 methods by Grimme et al. have been shown to give quite accurate thermochemistry for both covalently bonded systems and systems dominated by dispersion forces [9]. They also have advantages in that they have negligible computational time and add valuable insight, because the contribution from dispersion is easily separated out from the standard DFT energy. Another advantage is that when using DFT-D, the difference between various DFT functionals tends to decrease [9], as the global scaling factor 6 is larger for more local functionals and smaller for functionals with stronger longrange interactions. For PBE 6 = 0.75, and for revPBE the original parameters, for example, see [2,[11][12][13][14], and for parameters specifically modified for ionic surfaces [15]. The DFT-D3 method is a refinement of the DFT-D2 and has been implemented and evaluated in a plane-wave code recently [16]. The authors found that physisorption on Ag(111) was too strong, but this could be alleviated through inclusion of threebody terms. Other studies of adsorption on metal surfaces using dispersion corrected DFT have found that adsorption energies are overestimated [17], unless 6 parameters are decreased by fitting to ab initio data for metal clusters [18] or by including screening effects [19,20].
vdW-DF has been used as a nonlocal correction to adsorption for a few metallic systems. Examples include benzene on noble metals [21], benzene and carbon monoxide on Au stepped surfaces [22], azobenzene on Ag(111) [23], and carbon monoxide on Pt(111) [24].
The focus of our paper is self-assembled monolayers (SAMs) of alkanethiols on gold [25]. Experimental data show that gold adatoms are part of the bond between the SAM and the gold surface [26][27][28][29][30]. Rather subtle changes to the thiol chain length and chemistry can also change the observed structure of the SAM [26]. Pure DFT for thiol chemistry on gold is far from being able to quantitatively reproduce the adsorption behavior [31], and results depend strongly on the chosen functional [32]. This could in fact primarily be a result of neglecting the dispersion contribution. Dispersion forces determine interchain interactions, and it is reasonable to assume that dispersion interaction is important for interactions between sulphur and large metal atoms such as gold (the 6 coefficient for S is more than three times larger than that for C). The overly attractive PBE functional includes dispersion to a small degree and predicts weak binding in noble gas dimers, where more repulsive functionals like BLYP or revPBE do not bind at all [3].
In this paper we explain the main reasons why semiempirical dispersion treatment based on atom pairwise potentials is inadequate for metallic systems and propose a simple method for how to build a physically sound yet simple model for including dispersion forces for metals based on the DFT-D method. The method is not as rigorous as the explicit method in [19] but captures many of the same properties and has the advantage that it can be used together with the common method DFT-D2, which is available in most planewave codes. We apply it to physisorption of some linear alkanethiols as well as chemisorption of methylthiolate on Au(111) and show that dispersion forces account for about half of the adsorption energy of methylthiolate.

Computational Details
We start by noting that in the DFT-D atom pairwise method no screening of the dispersion interaction is taken into account. Most of the net dispersion energy comes from frequencies of the polarizability that lie in the visible part of the electromagnetic spectrum ( Figure 6.3 [38]). The plasma frequency of metals is also in the visible region [39], and therefore metals will quite efficiently screen dispersion interactions, making them even more short-ranged than non-metallic systems. Although this screening depends on the whole spectrum of dynamic polarizabilities in the system, we simplify the general behavior into a model that is much easier to implement in any plane-wave code with the DFT-D2 method available. We treat the valence electrons as being inactive for dispersion interactions and assume that only the bound core electrons contribute to the atom pairwise interaction in DFT-D. The valence electrons constitute the metallic states and will screen the dispersion interaction between the ionic cores. We introduce a hard cutoff for London dispersion of 12 bohr to simulate the screening of the valence electrons. We also replace atomic 6 coefficients and dispersion radii for any transition metal by the noble gas in the row above that metal in the periodic table. Because the DFT-D2 method only includes 6 parameters, no higher terms in the dispersion expansion have been considered. This might be necessary for similar treatment of the DFT-D3 method. The screening distance is chosen, rather arbitrarily, as the distance where the calculated lattice parameters for PBE-D and revPBE-D agree for gold, with the goal in mind to minimize differences between density functionals. The screening represented in our work by a hard cutoff has been included more rigorously using the Lifshitz-Zaremba-Kohn (LZK) theory for dispersion interaction between an atom and a solid surface [19]. The many-body effect of the metallic screening can be rewritten in terms of pairwise potentials with reduced atomic 6 parameters. This suggests that our simple alternative approach for a pairwise method is reasonable. Our approach treats all transition metals in the same row equally, while the LZK theory shows that the screening varies slightly depending on the metal.
Our method includes the major physics, makes use of existing implementations of current DFT packages, and provides a quick set of consistent parameters for all transition metals for use as is or as a starting point for further parameter refinement. Although in principle a hard cutoff makes the potential energy surface discontinuous, no problems were encountered during any geometry optimizations in this work. In a few occasions the BFGS optimization procedure failed, but using damped dynamics worked for all cases. Because of the cutoff introduced to handle screening from metallic states any long-range asymptotic interactions are missing for the parts of the simulation cell, which have no metal atoms, such as interchain interactions in SAMs. Both of these effects could be improved by making the cutoff via a smooth function that also depends on the position in the cell (i.e., is only active inside the metal).

DFT Calculations
All DFT calculations were performed using plane waves and ultrasoft pseudopotentials using Quantum Espresso (QE) version 4.2.1 [41]. The kinetic energy cutoff was 25 Ry and the density cutoff 250 Ry for all adsorption studies and 35/350 Ry for bulk calculations. Adsorption energies were converged to within 0.02 eV using these parameters. All pseudopotentials were taken from the QE pseudopotential library [42]. Both the PBE [43] and revPBE [44] functionals have been used in the paper. All surface slab calculations are made at the equilibrium lattice parameter for the corresponding functional and included four metal layers and at least 10Å vacuum between the slabs. A Gaussian smearing of 0.01 Ry was used for the Brillouin zone integration. The modules/mm dispersion.f90 file of QE was modified before the program was compiled: the 6 coefficients and dispersion radii were taken from the DFT-D2 data in the dftd3.f file downloadable from the DFT-D3 web page [45]. The parameters for all transition metals were replaced by the noble gas in the row above. We also had to increase the parameter × by a factor of 16 to handle the short cutoff values we use for London dispersion. The modified mm dispersion.f90 files are available upon request from the author.
Bulk calculations were made for the six metals Cu, Ag, Au, Ni, Pd, and Pt, and we calculated lattice parameters and cohesive energies using a 6 × 6 × 6 Monkhorst Pack k-point mesh [46]. Calculations for Ni were spin polarized.
We have studied CO adsorption on the close packed (111) surfaces of the fcc metals Cu, Ag, Au, Pd, and Pt. The choice of metals covers the range from very weak adsorption to strong chemisorption and experimentally different adsorption sites. The CO coverage was 0.25 monolayers, and a 2 × 2 surface unit cell was used. The Monkhorst-Pack grid for k-points was 4 × 4 × 1. All CO adsorption energies have been corrected using an empirical method based on the vibrational frequency related to the surface coordination [33,47]. The correction to the adsorption energy is 1.8 − 0.0008 * ] CO , where ] CO is the internal stretch frequency of the adsorbed CO molecule. We used the vibrational frequencies from [33]. The adsorption energy is as a result corrected by (on average) +0.11 eV for top adsorption, +0.25 eV for bridge adsorption, and +0.32 eV for hollow adsorption.
Benzene adsorption on Ag(111) and Au(111) serves as an example of physisorbed systems for which dispersion is the main mode of binding and adsorption on Pt(111) as a system where benzene adsorbs quite strongly. These calculations were performed with a 4 × 4 unit cell and a 2 × 2 × 1 Monkhorst Pack k-point mesh.
We have also studied methylthiolate adsorption on Au(111) in several adsorption geometries: on the (111) surface, adsorbed on top of an Au adatom on (111), and as a dithiol-Au adatom complex on (111). All adsorption energies presented for adatom systems include the formation energy of the adatom from a reservoir of bulk gold atoms [31]. The most stable adsorption geometry is the cis conformation of the dithiol-Au-adatom complex, consistent with a previous study [31]. The complex is shown in Figure 2. Adsorption energies ads were calculated in a standard fashion: This means that the more negative the adsorption energy, the stronger the molecule adsorbs on the surface.

Bulk Metals.
The inclusion of dispersion forces in our modified DFT-D2 method improves lattice parameters and cohesive energies for both PBE and revPBE ( Table 1). The mean deviation, mean absolute deviation, and maximum deviation of the lattice parameters are all improved, and the results for the two dispersion corrected functionals are similar. The calculated cohesive energies behave similarly. This indicates that bulk properties of the late transition metals improve using the dispersion correction, with no major shortcomings. Furthermore, the difference between the revPBE and PBE functionals decreases significantly when dispersion is included. Unmodified DFT-D2 improves lattice constants and worsens cohesive energies compared to nondispersion corrected functionals, while our modified version improves both lattice parameters and cohesive energies.

Chemisorption of CO and Physisorption of Benzene.
Adsorption energies for CO on transition metals are significantly improved for revPBE-D as compared to revPBE, whereas results on average are unchanged for PBE-D as compared to PBE (Table 2(a)). PBE gives stronger chemisorption than revPBE, consistent with a previous comparison between similar functionals PW91 and RPBE [33]. Weak adsorption is improved for PBE-D whereas strong chemisorption is overestimated. Again the difference between the two functionals is greatly reduced. The dispersion correction does not change the predicted adsorption site on either Pd or Pt, where multiple sites were investigated. We performed a calculation for CO on Pt(111) using the vdW-DF2 functional as well, which gives an adsorption energy of +0. 23   whereas our results indicate that a self-consistent calculation gives quite different results. We reoptimized the geometries for adsorption of CO on Pt(111) with PBE (no dispersion), which gave adsorption energies −1.43 eV (top) and −1.42 eV (hollow). This shows that the adsorption geometries for chemisorbed CO are not changed significantly from the PBE-D geometries, in particular considering that Pt has a strong dispersion interaction. We also performed a set of calculations using PBE-D (unmodified Grimme) and found adsorption energies of −1.92 eV (top) and −1.84 eV (hollow). The unmodified dispersion correction thus strongly overbinds this system and is less suitable for adsorption onto metal surfaces if used as is.
For adsorption of benzene on Ag(111) both functionals give good results, overestimating the experimental adsorption energy (−0.42 eV [34,35]) by about 0.1 eV (Table 2(b)). On Au(111), the error is slightly larger; PBE-D overestimates the experimental adsorption energy (−0.62 eV [36]) by 0.13 eV and revPBE-D by 0.28 eV. As expected, the entire binding energy is a result of dispersion. The calculations for strongly bound benzene on Pt(111) show that both functionals match the experimental adsorption energy −1.70 eV very well (Table 2(b)). Our approach gives equally good or better results for benzene adsorption on Au(111) and Ag(111) than Møller-Plesset 2nd order perturbation theory (MP2) [35], at a much smaller computational cost. We have also performed calculations for benzene on Au(111) with a cutoff for London dispersion of 200 Bohr, and the adsorption energy is then overestimated by ∼0.3 eV for PBE-D and ∼0.55 eV for revPBE-D, showing the importance of keeping the dispersion forces short-ranged via the hard cutoff.

Self-Assembled Monolayers of Thiols on Au(111).
We start our discussion of alkanethiols on Au(111) with the physisorption of methanethiol, CH 3 SH, up to butane thiol, CH 3 (CH 2 ) 3 SH. The experimental physisorption energy is linearly dependent on the chain length for nonbranched thiols [40], and the slope of the line is thus related to the strength of the CH 2 -Au interaction, whereas the -intercept is related to the SH-Au interaction [40]. Our calculated physisorption energies for the thiols on gold are too strong (Figure 3(a) and Table 3(a)), as was physisorption of benzene on gold. Again revPBE-D overestimates more than PBE-D. Both the slope and the -intercept deviate from experiment, which indicates that both Au-C and Au-S dispersion interactions are too strong. No adatom geometries were investigated in this paper for physisorbed thiols.
We have also calculated the dissociative adsorption energy of dimethyl disulfide adsorbed on Au(111) as RS-, RS-Au-adatom and as RS-Au-SR adatom-dithiol moiety (Table 3(b)). The RS-Au-SR has previously been found to be the most stable adsorption geometry by DFT calculations [31], and the dissociative chemisorption energy of disulfides is independent of chain length [40], which results in a straight line for chemisorption in Figure 3. We also find the cis-RS-Au-RS adsorption geometry to be the most stable one. Table 3: (a) Physisorption energies of ethanethiol and butanethiol on Au(111). All energies are in eV. The physisorbed R-SH species are referenced against corresponding gas phase thiol. 6 coefficients are in Ry * Bohr 6 . (b) Dissociative adsorption energy of dimethyl disulfide as a dithiol-adatom surface complex on Au(111) including the formation of adatom species from bulk Au atoms (Figure 2). The PBE and revPBE results are obtained from the PBE-D and revPBE-D results by subtracting the contribution from dispersion forces. 6 coefficients are in Ry * Bohr 6 . The adsorption energy is overestimated by 0.35 eV (Figure 3(a)),but agreement with experiments is better than for nondispersion corrected functionals [31]. Our results qualitatively reproduce TPD results of small thiol and disulfide measurements. Both physisorption and dissociative chemisorption are too strong, but the balance between the two is reasonable. If we extrapolate our results, the physisorption and the dissociative chemisorption lines cross at a chain length of 9 carbon atoms for revPBE-D and 11 carbon atoms for PBE-D. This is to be compared to experimental data for which the lines cross for linear alkanethiols with 14 carbon atoms [40].
Physisorption of benzene and thiols on Au(111) as well as chemisorption of thiolate on Au(111) are overestimated. This coupled with the fact that revPBE-D overestimates more than PBE-D implies that our simple approximation leads to an overestimation of surface-adsorbate dispersion interaction between any molecule and gold. This in turn implies that the systematic error stems from the dispersion parameters for Au and that the 6 parameter is too high.
As a simple refinement and incentive to other groups, we try a 6 coefficient for Au that is half as high. The adsorption energy of benzene on Au(111) becomes −0.51 eV, which is 0.1 eV higher than the experimental value, but on the other hand both the PBE-D and revPBE-D functionals give the same result, which is a strong indication that the 6 coefficient is more reasonable after modification. More importantly, the result for thiol SAM's on gold is highly encouraging (Figure 3(b)).
Both functionals now obtain the right balance between physisorption of alkanethiols and chemisorption of dialkanedithiols, with a crossing occurring for chain length = 14 where physisorption is as strong as chemisorption, in agreement with experiments. The revPBE-D results agree with experiments, while PBE-D only slightly overestimates both physisorption (0.1 eV) and chemisorption (0.2 eV).
Quite importantly, our results indicate that the chemisorption energy of thiols binding as a thiolate via a dithiolate-Au-adatom surface complex is strongly influenced by dispersion interactions (65% of the total adsorption energy for revPBE-D and 40% for PBE-D using the modified DFT-D2 parameters of Table 4). Probably the only reason why PBE is moderately successful for SAM calculations is because it contains some long-range effects in its construction [3], whereas revPBE does not. It also strongly suggests that any accurate theoretical treatment of thiols adsorbed on gold needs to take into account dispersion. Our parameters that do so simultaneously for PBE-D and revPBE-D for the thiol/gold system are presented in Table 4. They should work for any molecule adsorbing onto gold surfaces.

Conclusions
We have modified the DFT-D method of Grimme [7] for use in metallic systems by introducing a hard cutoff of 12 Bohr to the dispersion interaction and replacing the 6 coefficients for metal atoms by those of the noble gas in the row above. Our model takes into account screening of the dispersion forces by the conducting valence electrons in the metal, includes the main physics, and reproduces a wide variety of experimental data for both bulk metallic systems as well as adsorption onto metal surfaces. It is particularly useful in systems where both chemisorption and dispersion interactions compete and need to be taken into account simultaneously, such as self-assembled thiol monolayers on gold.
Our calculated adsorption energies are in agreement with experimental data for both physisorbed and dissociatively chemisorbed thiols on gold. Furthermore, our calculations show that including dispersion is crucial in order to obtain that level of accuracy for this kind of system.