Molecule-surface interaction from van der Waals-corrected semilocal density functionals: the example of thiophene on transition-metal surfaces

TL;DR

This paper evaluates various semi-local density functionals, especially with vdW corrections, for accurately modeling thiophene adsorption on transition-metal surfaces, emphasizing the importance of balanced short- and long-range interactions.

Contribution

It demonstrates the effectiveness of vdW-corrected semi-local density functionals in predicting adsorption energies and geometries for thiophene on metal surfaces.

Findings

vdW corrections improve adsorption energy predictions

SCAN functional estimates intermediate-range vdW effects well

balanced short- and long-range correlation is crucial for accuracy

Abstract

Semi-local density functional approximations are widely used. None of them can capture the long-range van der Waals (vdW) attraction between separated subsystems, but they differ remarkably in the extent to which they capture intermediate-range vdW effects responsible for equilibrium bonds between neighboring small closed-shell subsystems. The local density approximation (LDA) often overestimates this effect, while the Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation (GGA) underestimates it. The strongly-constrained and appropriately normed (SCAN) meta-GGA often estimates it well. All of these semi-local functionals require an additive non-local correction such as the revised Vydrov-Van Voorhis 2010 (rVV10) to capture the long-range part. This work reports adsorption energies and the corresponding geometry of the aromatic thiophene (C$_4$H$_4$S) bound to transition metal surfaces. The adsorption process requires a genuine interplay of covalent and weak binding and requires a simultaneously accurate description of surface and adsorption energies with the correct prediction of the adsorption site. All these quantities must come from well balanced short and long-range correlation effects for a universally applicable method for weak interactions with chemical accuracy.

Tables

Table 1. Table 1: The lattice constants (in Å) of substrate metals, compared to the experimental zero-point phonon corrected lattice constants.

	Copper	Silver	Gold
Reference57	3.60	4.06	4.06
PBE	3.63	4.15	4.16
PBEsol	3.57	4.05	4.08
SCAN	3.56	4.08	4.09
revSCAN	3.57	4.11	4.11
PBE+rVV10	3.62	4.11	4.13
PBEsol+rVV10	3.55	4.02	4.06
SCAN + rVV10	3.54	4.06	4.07
revSCAN + rVV10	3.55	4.06	4.08

Table 2. Table 2: Distances (in Å) between the sulphur atom in thiophene and the nucleus of the nearest atom of the metal surface. The PBE+vdW-dZK method was presented in ref. 40 and 45.

	Cu		Ag		Au
	d(Cu-S)	d(C-S)	d(Cu-S)	d(C-S)	d(Cu-S)	d(C-S)
Expt55, 56, 59	2.62$\pm$0.03	1.71	-	1.71	-	1.71
PBE+vdW-dZK	2.57	1.71	3.16	1.71	3.23	1.71
SCAN	2.70	1.71	3.00	1.71	2.98	1.71
revSCAN	2.88	1.71	3.03	1.70	3.02	1.70
SCAN + rVV10	2.57	1.71	2.97	1.71	2.93	1.71
revSCAN + rVV10	2.56	1.71	2.93	1.71	2.91	1.71
PBE + rVV10	2.88	1.71	3.00	1.71	2.98	1.71
PBEsol + rVV10	2.19	1.71	2.68	1.71	2.59	1.71

Table 3. Table 3: Adsorption energy (in eV) of thiophene on Cu(111), Ag(111) and Au(111) surfaces with GGA and meta-GGA-based approximations with respect to the most stable fcc-45 adsorption site. The estimated (see the supplementary information for more details) adsorption energies ( ± plus-or-minus \pm 0.2 eV) 38 , 24 from the TPD temperature maxima (Expt.), 59 , 61 , 66 and the results from the PBE+vdW-dZK model 40 are also shown.

Cu	Ag	Au
Expt59, 61, 66	-0.66	-0.52	-0.64
PBE+vdW-dZK	-0.60	-0.50	-0.56
SCAN	-0.44	-0.41	-0.45
revSCAN	-0.38	-0.34	-0.35
SCAN + rVV10	-0.83	-0.74	-0.81
revSCAN + rVV10	-1.06	-0.88	-0.95
PBE + rVV10	-0.61	-0.55	-0.63
PBEsol + rVV10	-1.22	-0.93	-1.06

Table 4. Table 4: Computed surface energies (in J/m 2 ) of different transition metals

Metal	PBE	PBEsol	SCAN	revSCAN	PBE+rVV10	PBEsol+rVV10	SCAN+rVV10	revSCAN+rVV10	Reference76, 77, 75
Cu	1.41	1.76	1.71	1.73	1.74	2.11	1.92	2.07	1.79 $\pm$ 0.19
Au	0.82	1.13	1.06	1.06	1.17	1.52	1.29	1.43	1.51 $\pm$ 0.16
Ag	0.80	1.09	1.04	1.04	1.11	1.43	1.24	1.36	1.25 $\pm$ 0.13
ME	-0.51	-0.19	-0.25	-0.24	-0.18	0.17	-0.03	0.10
MAE	0.51	0.19	0.25	0.24	0.18	0.17	0.12	0.15

Table 5. Table S1: Adsorption energy (in eV) of thiophene over Cu(100), Ag(100), and Au(100) surfaces with the meta-GGA based approximations and their corresponding vdW corrected versions with respect to the most stable adsorption site. Hollow-45 (the number 45 represents the value of angle of rotation) is found to be the most stable site of adsorption in 100 surfaces irrespective of the nature of metals and exchange-correlation(XC) functionals used.

	Cu	Ag	Au
Reference62	-0.63	-	-
SCAN18	-0.89	-0.49	-0.64
revSCAN22	-0.76	-0.40	-0.45
SCAN + rVV1016	-1.30	-0.82	-1.01
revSCAN + rVV10	-1.47	-0.94	-1.07

Table 6. Table S2: Adsorption distances (in Å) between the sulfur atom in thiophene and the nucleus of the nearest atom of the metal surface (100 surface) and the C-S bond length (in Å). Measurements are made on the hollow-45 site, which is the most stable site of adsorption in 100 metal surfaces.

	Cu		Ag		Au
	d(Cu-S)	d(C-S)	d(Cu-S)	d(C-S)	d(Cu-S)	d(C-S)
Reference55, 72	-	1.71	-	1.71	-	1.71
SCAN	2.13	1.77	2.89	1.71	2.54	1.72
revSCAN	2.17	1.76	2.97	1.70	2.98	1.70
SCAN + rVV10	2.13	1.77	2.75	1.71	2.45	1.72
revSCAN + rVV10	2.14	1.77	2.80	1.71	2.72	1.71

Table 7. Table S3: Adsorption energy (in eV) of thiophene on Cu(110), Ag(110), and Au(110) surfaces with meta-GGA based approximations and their corresponding vdW corrected versions, with respect to the most stable site of adsorption. Shortbridge-45 is found to be the most stable site of adsorption irrespective of the nature of metal, except for revSCAN which predicts longbridge-0 as the most stable site of adsorption. We did’t found the results from experiments to compare to our results.

	Cu	Ag	Au
SCAN	-0.74	-0.49	-0.74
revSCAN	-0.53	-0.38	-0.49
SCAN + rVV10	-1.12	-0.77	-1.05
revSCAN + rVV10	-1.23	-0.87	-1.06

Table 8. Table S4: Adsorption distances (in Å) between the atom in thiophene and the nucleus of the nearest atom of the metal surface (110 surface) and the C-S bond length (in Å). Measurements are made on the most stable site of adsorption in 110 metal surfaces.

	Cu		Ag		Au
	d(Cu-S)	d(C-S)	d(Cu-S)	d(C-S)	d(Cu-S)	d(C-S)
Reference55	-	1.71	-	1.71	-	1.71
SCAN	2.24	1.71	2.36	1.71	2.73	1.71
revSCAN	2.43	1.71	2.53	1.7	2.82	1.71
SCAN + rVV10	2.24	1.71	2.35	1.71	2.67	1.72
revSCAN + rVV10	2.27	1.72	2.48	1.71	2.66	1.71

Table 9. Table S5: Surface energy(in J / m 2 𝐽 superscript 𝑚 2 J/m^{2} ) of different surfaces of copper, gold and silver. The average of the surface energies of three different surfaces are compared against the reference 76 , 77 values in Table 4 in the main text. rVV10 corrections on meta-GGA’s significantly improve the surface energies. Both SCAN+rVV10 and revSCAN+rVV10 overall perform better for surface energies.

Metal	surface	Functionals
		PBE	PBEsol	SCAN	revSCAN	PBE+rVV10	PBEsol+rVV10	SCAN+rVV10	revSCAN+rVV10	Reference76, 77
	111	1.28	1.60	1.56	1.58	1.61	1.95	1.77	1.91
Cu	110	1.51	1.87	1.83	1.84	1.84	2.23	2.04	2.19	1.79 $\pm$ 0.19
	100	1.44	1.79	1.75	1.75	1.77	2.15	1.97	2.10
	111	0.69	0.97	0.91	0.91	1.04	1.35	1.13	1.28
Au	110	0.90	1.24	1.15	1.15	1.26	1.63	1.39	1.54	1.51 $\pm$ 0.16
	100	0.87	1.18	1.11	1.10	1.22	1.57	1.34	1.47
	111	0.73	0.99	0.95	0.96	1.03	1.33	1.15	1.28
Ag	110	0.86	1.16	1.11	1.11	1.17	1.51	1.31	1.44	1.25 $\pm$ 0.13
	100	0.81	1.10	1.06	1.05	1.12	1.44	1.25	1.36

Table 10. Table S6: Calculation of adsorption energies of thiophene over Cu(111), Au(111) and Ag(111) using Redhead’s Analysis. With the information of peak maximum temperature ( T m subscript 𝑇 𝑚 T_{m} ) and the heating rate ( β 𝛽 \beta ) utilized in the experiment’s 59 , 61 , 66 such energies can be estimated.

	$T_m$	$\beta$	Adsorption energy
Surface	(K)	(K/s)	(eV)
Cu(111)	23459	0.5	-0.66
Au(111)	24061	3	-0.64
Ag(111)	19066	1	-0.52

Equations

$$f_x(\alpha )=e^{-c_{1x}\frac{\alpha }{1-\alpha }}.\theta (1-\alpha )-d_x{e}^{\frac{C_{2x}}{1-\alpha }}.\theta (\alpha -1)$$

$$f_c(\alpha )=e^{-c_{1c}\frac{\alpha }{1-\alpha }}.\theta (1-\alpha )-d_c{e}^{\frac{C_{2c}}{1-\alpha }}.\theta (\alpha -1)$$

$$E_c^{nl}=\int drn(r)[\frac{1}{2}\int d{r}^{\prime }\varphi (r,r^{\prime })n(r^{\prime })+\beta ]$$

$$E_{vdW}=[-\frac{C_3}{(z-z_0{)}^3}-\frac{C_5}{(z-z_0{)}^5}]f_d$$

$$f_d=\frac{x^5}{(1+g{x}^2+h{x}^4+x^{10}{)}^{\frac{1}{2}}}$$

$${\alpha }_{renormalizedatom}=\frac{{\alpha }_{(freeatom)}}{4{\alpha }_{(C)}+4{\alpha }_{(H)}+{\alpha }_{(S)}}{\alpha }_{(thiophene)}$$

$$E_{ads}=E_{surface+thiophene}-E_{surface}-E_{thiophene}$$

$$\sigma =\underset{n\to \infty }{\lim }\frac{E(n)-nN{ϵ}_{bulk}}{2A}$$

$$E_{ad}=-R{T}_m[ln(\frac{\nu T_m}{\beta })-3.46]$$

Full text

\phone

+1 267 909 1660

Temple University] Department of Physics, Temple University, Philadelphia, PA, 19122

Department of Inorganic and Analytical Chemistry] Department of Chemistry, Budapest University of Technology and Economics, Budapest, Hungary, H-1111

Temple University] Department of Physics, Temple University, Philadelphia, PA, 19122

1 Introduction

The development and assessment of various vdW methods has been an intensive area of research in the past decade. 1, 2, 3, 4, 5, 6 Now the scientific community possesses a broad range of approximations7, 8, 9 of useful but limited accuracy. vdW methods approximate the long-range correlation which arises from the physics of collective plasmon oscillations. Wavefunction-based approximations such as Coupled Cluster (CC) methods naturally include vdW interactions, but are practically beyond the reach of the condensed matter community at this time. Alternatively the Random Phase Approximation (RPA) is nearly exact for the long-range, and is regarded as a benchmark to assess vdW methods in condensed matter physics.10 Even with increased computational power and increasingly efficient implementations, RPA has a limited practicality for materials.11

Most current vdW methods have been developed within the density functional context, in which the construction of self-consistent orbitals is a parallelizable step. Semilocal density functional theory (DFT) can be accurate for the short-range correlation, but misses the long-range part. The long-range vdW component is captured by either pairwise vdW methods12, 13, 14 or by nonlocal correlation functionals.7, 15, 16 Each of these approximations is then added to an appropriately chosen semilocal exchange-correlation functional. The vdW-DF7 and VV1015 non-local correlation functionals are based on approximations for the polarizability, and VV10 has a fitted short-range attenuation parameter (b) that adapts to the semilocal term. Many of the current vdW models are reasonably accurate and efficiently applicable to geometries and binding energies.

The recent rVV1016 correlation functional has been combined17 with the SCAN18 meta-GGA and has been successfully applied to various systems including interlayer binding energies, adsorption energies and structural properties.19 One major advantage of this method is its computational efficiency. Although SCAN+rVV10 delivers a generally reasonable description for various properties, it gives a disappointing treatment for some others. Examples include the overstructured radial distribution functions in liquid water,20 inaccurate structural and mechanical properties in PPTA,21 and inaccurate prediction of ground state properties of MnO and CoO.23 While SCAN captures intermediate-range vdW interactions, it may capture too much. revSCAN,22 a revised version of SCAN was constructed to diminish the intermediate-range vdW interaction.

This work explores the accuracy and precision of the SCAN+rVV10 and the revSCAN+rVV10 approximations for thiophene adsorption on the surface of coinage metals. For comparison we also include several GGA-based semilocal exchange-correlation functionals with rVVV10 correction into this assessment.

The adsorption of molecular species on metal surfaces is a relevant problem for both computational simulations and industry. In general the adsorption of organic species on metal surfaces can be a synergy of chemi - and physisorption, however, recent works on the adsorption of benzene on the surface of coinage metals reveal the large role of weak interactions.24 A recent work reported accurate SCAN+rVV10 binding energies for the adsorption of benzene on transition metal surfaces.17 A universally accurate approximation can be expected to capture adsorption sites, surface and binding energies simultaneously in the adsorption process.

The thiophene molecule is the smallest aromatic sulfur-containing compound. It is a natural choice as a test case for simulations. Thiophene is also a good test to study reactions that follow the catalytic desulfurization on metal or semiconducting surfaces. The adsorption of the thiophene molecule on metal surfaces turns out different than the extensively studied case of benzene adsorption. Depending on the underlying metal, the adsorption of thiophene on the metal surface can show some chemisorption character25 whose description requires a very accurate balance of local and nonlocal correlation in the meta-GGA and its partner van der Waals approximation.26 A simultaneously correct description of the adsorption energies and sites is a challenge to density functional theory.27

 

**Meta-GGA density functional approximations

**Exchange-correlation approximations within DFT can be broadly categorized by the different ingredients of the 5 rungs of a ladder of increasing nonlocality.28 The accuracy of an approximation usually increases when more ingredients are included. The increased accuracy often comes at the price of deteriorated efficiency, especially when non-local information is included.

Among the most accurate density functional approximations are the meta-generalized approximations (meta-GGA’s). Meta-GGA’s constitute the third-rung of the ladder. Commonly used meta-GGAs include one more ingredient beyond the GGA level, the kinetic energy density $\tau(r)=\frac{1}{2}\sum_{i=0}^{occ}|\nabla\phi_{i}|^{2}$ where the $\phi_{i}$’s are Kohn-Sham orbitals. Different dimensionless variables built from $\tau(r)$ have been proposed for use in meta-GGAs. The most successful dimensionless variable so far29, 30 is $\alpha(r)=\frac{\tau(r)-\tau^{vW}(r)}{\tau^{UEG}(r)}$, where $\tau(r)^{vW}=\frac{|\nabla n(r)|^{2}}{8n(r)}$ is the von-Weiszäcker kinetic energy density of a single-orbital system and $\tau^{UEG}(r)=\frac{3}{10}{(3\pi^{2})}^{\frac{2}{3}}n(r)^{\frac{5}{3}}$ is the kinetic energy density of the uniform electron gas. The $\alpha(r)$ variable is an iso-orbital indicator, recognizing different types of orbital overlap environments and directly related to the electron localization function. Single-orbital regions are identified by $\alpha(r)=0$, highly-overlapped metallic regions with slowly-varying electron densities are revealed by $\alpha(r)\approx 1$, and regions of non-bonding density overlap by $\alpha(r)>>1$.

Madsen et al.31 showed that inclusion of the kinetic energy densities enables meta-GGA’s to distinguish between dispersive and covalent interactions. A family of nonempirical constructions27 led to the development of the SCAN18 meta-GGA. SCAN satisfies 17 exact constraints, while preserving the ability to capture intermediate-range weak interactions. With tests and assessments on diverse systems, the SCAN meta-GGA has been a success-story in the past three years.18, 19, 32, 33, 34

Through the $\alpha$ dependence of the interpolation functions for exchange and correlation energy, SCAN can capture intermediate-range dispersive interactions as shown in the following equations:

[TABLE]

[TABLE]

Many physical situations require the long-range part of the correlation, mathematically described by a double integral in the three-dimensional space, and not captured by any semilocal density functional. The vdW correlation functional by Vydrov and Van Voorhis (VV10)15 and rVV1016 are the examples for a long-range functional that allows the nonlocal correlation energy and its derivatives to be efficiently evaluated in a plane wave framework, as pioneered by Román-Pérez and Soler.35 The long-range correlation is a double integral:15, 16

[TABLE]

The VV10 and rVV10 corrections are designed to vanish for the uniform electron gas. This feature makes it possible to pair the nonlocal correlation energy with the semilocal exchange-correlation energy by utilizing a parameter to damp the intermediate and short range contribution of the former. A critical ingredient in the kernel is the local band gap, a quantity that accounts for density inhomogeneity and makes VV10 and rVV10 applicable for molecular systems. Like VV10, the rVV10 correction has two adjustable parameters C and b inside the kernel $\phi$ that allow it to adopt to any semilocal functional. The values of the C and b parameters for rVV10 for SCAN were 0.0093 and 15.7 respectively.

A universally applicable and accurate vdW approximation should benefit from the interplay of the nonlocal and semilocal functionals. Aside from the particular form of the vdW correlation functional, the choice of the exchange functional has received considerable attention within this work. The choice of semilocal exchange has already attracted interest in the context of GGA density functional. The revPBE-GGA exchange functional chosen for vdW-DF often leads to too-large intermolecular binding distances and inaccurate binding energies.36 To remedy these problems, alternative ”less repulsive” exchange functionals have been proposed, which when incorporated within the vdW-DF scheme lead to much improved accuracy.

Earlier attempts emphasized the improvement of vdW-DF by exploring and developing alternatives to the original revPBE exchange. These studies were limited to PBE-based GGA’s, and the underlying semilocal exchange was fitted to the vdW functional in an empirical manner.

 

**Benchmark binding energies for the adsorption of thiophene on metals

**To properly assess the limitations of our approximations, we need accurate benchmark adsorption energies. After appropriate calibration, temperature programmed desorption (TPD) or thermal desorption spectroscopy can be used to evaluate the activation energy of desorption. The binding energy might be estimated from the temperature of maximum desorption via Redhead’s analysis.37 However, the estimated binding energies might display an uncertainty larger than the chemical accuracy of 0.04 eV required for an accurate description of the adsorption.38 A considerably more accurate complete analysis method would lead to more accurate results,37, 38 but no such results are available for thiophene on coinage metal surfaces according to our knowledge. The nonlocal random phase approximation (RPA)10, 39 could be a reliable reference for long-range vdW interactions. RPA calculations are, however beyond reach for large supercells at this time. Here in this work we use an approximation that is robust enough and mimics the RPA binding energies almost perfectly for the interaction of graphene and metal surfaces. This approximation40, which we will call PBE+vdW-dZK from this point onwards, models the long-range van der Waals correction for physisorption of graphene on metals with the damped Zaremba-Kohn (ZK)41, 42 second-order perturbation theory. In this model, quadrupole-surface interactions and screening effects are included. This model relies on accurate static polarizabilities from higher-level calculations43, 44 and predicts the vdW interaction from the $C_{3}$ and $C_{5}$ coefficients and the distances between the particle and surface plane through an expression whose large -z asymptote is

[TABLE]

with z being the distance between the particle and the surface and $z_{0}$ the reference plane position. The damping factor for eqn.(4) is

[TABLE]

where, $x=\frac{z-z_{0}}{b}>0$, $g=\frac{2b^{2}C_{3}}{C_{5}}$ and $h=\frac{10b^{4}{C_{3}}^{2}}{{C_{5}}^{2}}$. The cutoff parameter ’b’ was choosen to be 3.3 bohr.45 Instead of taking the static dipole polarizability of the thiophene molecule, we base our $C_{3}$ coefficients on the ”renormalized atom” approach.45 The best polarizability for a particular atom (H, C, or S) in thiophene is then renormalized as

[TABLE]

With the static polarizabilities we can find separate $C_{3}$ and $C_{5}$ coefficients for each of the three elements in thiophene. The formula of renormalization that we are using is constructed for a ”particle” interacting with a metal surface. A molecule such as thiophene, even if not of a large size, cannot be treated as a particle. Here we treat it as a collection of renormalized atoms. Treating the whole thiophene as a particle would make its quadrupole polarizability grow roughly as the 5/3 power of its dipole polarizability, and would overestimate $C_{5}$ significantly.

 

**Parameterization of rVV10

**In the present work we are following the approach of parameterization of rVV10 for SCAN.17 Since the change of the value of C parameter does not significantly improve the binding curve,16, 17, 46 we keep the value of C fixed. But, we optimize the b value by fitting to the CCSD(T)47 binding energy curve of the Argon dimer. Notice that a recent empirical potential energy function for Ar dimer48 showed excellent agreement with CCSD(T) and CCSDT(Q) results.49 The use of such ab initio derived potential functions for reference can be justified.50 For all the calculations, the Argon dimer was placed in a cubic supercell of 25 Å. All the calculations were done using a single point gamma-centered k-mesh.

Figure 1 displays the binding energy curves of the Argon dimer from SCAN and revSCAN and the corresponding rVV10 corrected version with the CCSD(T) data.47 The revSCAN is more underbinding than SCAN in the intermediate range due to its construction, suggesting its need for a stronger van der Waals correction. We determine the b parameter for rVV10 for revSCAN by fitting to nine data points around the equilibrium distance with respect to the CCSD(T) binding curve. With 2.89% of mean absolute percentage error (MAPE), the b parameter was determined to be 9.4. This value is slightly smaller than the b = 9.8 suggested in the original revSCAN+VV10. This smaller value leads to stronger dispersion interaction. The SCAN+rVV10 has MAPE of 3.32% with the original b = 15.7.

For comparison we have combined the PBE51 and PBEsol52 GGA’s with the rVV10 correction shown in Figure 2. With MAPE of 1.76%, the b parameter for PBEsol+rVV10 is found to be 9.7 while 9.8 is the b parameter determined for PBE+rVV10. Surprisingly PBEsol gives less binding than PBE. One of the reasons for the larger b parameter value for PBE+rVV10 is its inability to give the minimum position correctly. While all other methods such as SCAN+rVV10, revSCAN+rVV10 and PBEsol+rVV10 give the minimum around 3.775 Å in agreement with CCSD(T), the PBE+rVV10 yields the minimum at around 4.0 Å. Our results show that PBEsol+rVV10 gives the best fit to Argon dimer followed by revSCAN+rVV10 and SCAN+rVV10 while PBE+rVV10 gives the relatively worst fit.

2 Computational details

All the DFT calculations were performed using the projector-augmented wave (PAW) formalism implemented in the Vienna ab initio simulation package (VASP) code. The lattice constants of silver, gold and copper were obtained by the geometry relaxations of their respective bulk structures using different XC functionals. (4x4) supercell of 111, 110 and 100 surfaces using optimized lattice constants were built in the atomic simulation environment (ASE).53, 54 The supercell has a five-atomic-layer thickness. In order to prevent the interactions due to the periodic images, a vacuum of 12 Å was added along the z-direction. To reduce the computational cost, the positions of the atoms on the bottom three layers were kept fixed, only allowing atoms on the two top layers to relax. The thiophene molecule was constructed using the reference C-S, C-C and C-H bond lengths55, 56 which was allowed to relax in the slab whose dimension was identical to that of the surface on which the adsorption occurs. Initially the thiophene was placed in a parallel orientation 3 Å above the top metal layer and was allowed to relax. The center of mass and the azimuthal angle of the thiophene was used for the classification of the geometry. High symmetry sites, namely top, hollow, bridge, shortbridge, longbridge, fcc and hcp,53, 54 were used as sites for adsorption, e.g., top-45 indicates the center of mass of the thiophene adsorbed on top of the metal atom with a symmetry axis rotated by 45*∘* from the direction of metal rows. The surface, the thiophene and the thiophene over the surface were all separately relaxed. PAW potentials as recommended in the VASP manual were used for all the calculations. A plane wave cutoff of 650 eV and a thermal-smearing temperature of $k_{B}T$ = 0.1 eV were chosen for both the bulk and surface calculations. The Brillouin zone was sampled using a 4x4x1 Monkhorst-Pack k-point set for surfaces while 20x20x20 were used in the case of bulk relaxations. The structure of thiophene was optimized before adsorbing over one side of the slab (i.e., coverage of 1/16). Since both energies and equilibrium distances are dependent on the site of the adsorption, all major high symmetry sites were chosen for relaxing thiophene over metallic surfaces. The adsorption energy was calculated by subtracting the energy of the combined system (surface + thiophene) from the energy of the surface alone and the energy of the thiophene alone.

[TABLE]

3 Results and discussion

3.1 Lattice constants of transition metals

In the adsorption process, the organic molecule binds to the metal surface. A full picture about the exchange-correlation (XC) approximations utilized must include the lattice constants, which impact the geometry during the adsorption process. We have assessed SCAN, revSCAN, PBE, PBEsol and their long-range van der Waals-corrected versions for the lattice constant of the three transition metals. Table 1 displays the lattice constants, which are compared to the zero-point-phonon corrected experimental lattice constants.57

The results in Table 1 show that PBEsol systematically shortens the lattice constants by roughly about 0.07 Å compared to the too-long PBE values making the PBEsol results quite accurate. Although SCAN is accurate for many bonding situations,19 it is not particularly accurate for the lattice constants of the transition metals studied here. The revSCAN results are slightly longer than the SCAN values making the revSCAN values slightly better for copper but worse for silver and gold. The application of rVV10 shortens the lattice constants. For copper it improves the too-long PBE results but worsens the already shorter PBEsol, SCAN and revSCAN results. Though this shortening effect helps to obtain excellent latice constants for Gold and Silver with SCAN+rVV10, revSCAN+rVV10, it is not enough to correct the too-long PBE values.

3.2 Adsorption energies and geometry of thiophene on Cu(111), Ag(111) and Au(111)

We have assessed the adsorption of the thiophene molecule on three crystal faces of copper, silver and gold, considering the adsorption energies, the adsorption geometry and the tilt angle between thiophene and the metal surface. Moving from Cu toward Au, the nature of the adsorption on these three metal surfaces changes. The adsorption on Cu(111) is a mixture of covalent and weak interactions, while the interaction on Au(111) is dominated by weak van der Waals interactions 58, 59, 60, 61. Though, the experiments do not give the precise value of adsorption energies59, they 62, 61 overall report a strong dependence on the coverage of the thiophene adsorption. The structural information of the adsorbed molecule on the metallic surface such as molecule-surface distance, the angle of the adsorbed molecule and surface, and the adsorption sites vary with increasing coverage.

Irrespective of the exchange-correlation functional and vdW correction applied, the adsorption of thiophene on (111) surfaces of the different metals displays some common features. Based on the given coverage, our calculations find that the fcc-45 is the most stable site of adsorption for all metals. The same adsorption site was found to be the most stable with our benchmark PBE+vdW-dZK approximation, and experiments support these results too. The predicted fcc-45 adsorption site for Cu(111) is close to the top adsorption site predicted by experiments for Cu(111)59, 63. The difference is just the result of the choice of the position of the reference point in thiophene.

According to the experimental results, an increased tilt angle is observed with increasing coverage, while lower coverage prefers a slightly lower tilt angle59 for thiophene adsorbed on Cu(111). Though a higher tilt angle of 55*∘* was found at a significantly higher coverage of thiophene on Au(111) than ours,64, 65 experiments indicate the preference of thiophene lying flat on both Ag(111) and Au(111) at low coverage.61, 66, 67

Our PBE+vdW-dZK method gives a tilt angle of 7*∘* - 17*∘* for thiophene adsorbed on Cu(111) surface, a tilt angle of 1*∘* - 2*∘* for thiophene adsorbed on Ag(111) and almost zero tilt angle on Au(111). Our other methods give a tilt angle of 1*∘* - 6*∘* for thiophene adsorbed on Cu(111), a tilt angle of 1*∘* - 2*∘* on Ag(111), and tilt angles of 1*∘* - 4*∘* on Au(111).

Experiments59 for adsorption of thiophene over copper suggest a range of tilt angle of 20*∘$\pm$3∘* at the coverage of 0.05 ML, and 25*∘$\pm$4∘* at the coverage of 0.1 ML. Our PBE+vdW-dZK results giving the tilt angle of 17*∘* at the most stable site at the coverage of 0.0625 ML, agree very well with the experiments.

The minimum Cu-S distance of 2.57 Å that we show in Table 2 from the PBE+vdW-dZK method based on LZK theory40 is in close agreement to the experiments.59 Though both SCAN and revSCAN predict slightly longer Cu-S distance, SCAN+rVV10 and revSCAN+rVV10 results are closer to the experiments.

We have not found precise experimental data for the distance between the adsorbed S atom of thiophene and the Ag(111) and Au(111) surface. Both SCAN and revSCAN and the corresponding rVV10 corrected methods yield almost identical adsorption distance irrespective of whether the surface is Ag(111) or Au(111). The PBE+vdW-dZK method gives a slightly larger distance of 3.23 Å for Au(111) and 3.16 Å for Ag(111), compared to other methods discussed here. However, the latter result gives a very good match with the previously studied68 PBE-vdWsurf method for the distance of Ag and S atoms for Ag(111). The adsorption distance predicted by SCAN+rVV10 is close to the results of Maurer et al.68 for thiophene adsorbed on Au(111).

The results in Table 3 show that our theoretical benchmark the PBE+vdW-dZK approximation adsorption energies are in a good agreement with the experimental adsorption energies estimated from TPD temperature maxima 59, 61, 66 using the Redhead’s model37 (see the supplementary information for more details). Our analysis shows that the adsorption energies estimated properly from Redhead’s model are considerably more precise for thiophene and benzene adsorption on coinage metal surfaces than suggested before38, 24 (cf. supplementary information). Notice the slightly different coverage given in most of the experiments.59, 58 Both SCAN and revSCAN underbind thiophene on Cu(111), Ag(111), and Au(111) surfaces compared to our theoretical reference. However, with the added rVV10 corrections they overbind. SCAN+rVV10 works less well for the thiophene adsorption on Cu(111) than for the adsorption of benzene.17 revSCAN by construction was designed to remove some of the intermediate range interactions of SCAN, so an underbinding of thiophene on coinage metal surfaces is expected, but surprisingly revSCAN+rVV10 is more overbinding than SCAN+rVV10. The reason behind the strong overbinding of revSCAN+rVV10 is the inclusion of relatively larger vdW correction through smaller b value. SCAN and in particular revSCAN are significantly underbinding for adsorption of thiophene on Ag(111) and Au(111) surfaces too, but adding the rVV10 corrections again leads to overbinding. The SCAN+rVV10 is overbinding compared to the earlier results from PBE+vdWsurf 68 and B86bPBE-XDM approximations38 too. Comparison with the relevant results of Christian et al.38 shows that B86bPBE-XDM results do not reflect the qualitative tendency that Cu and Au surfaces bind the thiophene about equally strongly and slightly stronger than Ag. This tendency is reproduced by all methods in Table 3, and especially well quantitatively reproduced by PBE+vdW-dZK, SCAN+rVV10, and PBE+rVV10.

For comparative purposes we show PBE+rVV10 and PBEsol+rVV10 results in Table 3. PBEsol is known to contain a certain amount of medium-range correlation compared to PBE. The combination of these approximations with rVV10 can serve as a simplified model of the more sophisticated meta-GGA’s with rVV10. Due to the lower-order gradient correction in the correlation energy and the decreased medium-range enhancement in its exchange the revSCAN meta-GGA resembles PBE, while SCAN exhibits more analogy with PBEsol. Being inspired by this analogy, we have computed the adsorption energy, distances and tilting of thiophene with PBE+rVV10 and PBEsol+rVV10.

Surprisingly PBE+rVV10 is more reliable than SCAN+rVV10 or revSCAN+rVV10. The adsorption energy on Cu(111) is overestimated by revSCAN+rVV10 and SCAN+rVV10, while PBE+rVV10 predicts less overbinding. This reliability of PBE+rVV10 is present for Ag(111) and Au(111). Adsorption energies on Ag(111) and Au(111) from PBE+rVV10 not only agree with the PBE+vdW-dZK results, but are very close to the estimated experimental values. The PBEsol+rVV10 approximation, although it is remarkably accurate for diverse properties including the binding energy of Xe on Cu(111) and Ag(111),69 turns out less successful in the case of adsorption of thiophene on Cu(111), Ag(111) and Au(111). PBEsol+rVV10 predicts too large binding energies and too short adsorption distances. The PBE+rVV10 method however is unable to yield the moderate tilting of thiophene over Cu(111), and is able to predict the almost parallel orientation over Ag(111) and Au(111) surfaces. The predicted adsorption distance from PBE+rVV10 is slightly longer than the value predicted by the experiments59 for thiophene over Cu(111). We are in lack of an accurate value of adsorption distance from experiments for thiophene over Ag(111) and Au(111) surfaces, but PBE+rVV10 values agree with earlier results from the $PBE+vdW^{Surf}$ method.68

It is more surprising that revSCAN+rVV10 significantly overbinds thiophene on all three transition metals, more than SCAN+rVV10 does. The naïve expectation from the combination of revSCAN and rVV10 is a more balanced description of weak interaction. The overall conclusion is that the combination of SCAN, revSCAN and rVV10 does not work accurately in general. A similar conclusion was drawn about the SCAN+MBD approximation.70 When the MBD method was combined with SCAN, the effective range of SCAN depended on system size.

Such long-range corrections need a different empirical cutoff parameter for each semi-local functional, in order to avoid a misrepresentation of intermediate-range vdW interaction. The pairing of semi-local and nonlocal terms can work well for some systems and fail for others.70 In particular, the pairing of SCAN with rVV10 works well for layered materials and for a benzene molecule adsorbed on coinage metals. But for liquid water, for some molecular crystals, and for the problem considered here, adsorption of thiophene on transition metals, this pairing overbinds. A familiar proposed solution would be to start from a semi-local functional that has little (PBE) or no (revSCAN) intermediate-range vdW interaction, and get the intermediate-range contribution from rVV10. This solution is consistent with the fact that semi-local functionals yield intermediate-range vdW attraction from their exchange energy terms. But we show here that this solution does not always work. While the PBE based vdW corrected methods performed better, revSCAN+rVV10 performed poorly. Both PBE+vdW-dZK and PBE+rVV10 methods gave overall better adsorption energies. PBE+vdW-dZK in particular was outstanding. It not only predicted the adsorption energies reasonably well, but also predicted the adsorption distance and tilting of the thiophene correctly.

Our results for the Cu, Ag and Au (100) and (110) surfaces can be found in the Supplementary Material. In general, the adsorption of thiophene on the (100) and (110) crystal faces is less explored experimentally. With all methods applied in our work, except SCAN we generally observe a stronger overbinding of the thiophene molecule on the Cu(100) than on Cu(111) surface.71, 72, 73, 74 Considering the hollow-45 site as the most stable adsorption site, the SCAN and revSCAN adsorption energy of -0.89 eV and -0.76 eV are reasonably accurate even without vdW correction, while SCAN+rVV10 and revSCAN+rVV10 both overbind. The revPBE-vdW and optB88-vdW approximation were found give the similar trend predicting -0.50 eV and -0.85 eV binding energy, respectively.73, 74 Among the approximations that we have applied, the revSCAN approaches these values the best, giving -0.76 eV binding energy at the hollow-45 site.

3.3 Surface energies of the three transition metals with GGA and meta-GGAs, and their vdW-corrected versions

Practically, due to cancellation in Eq. 7, the surface energy makes a minor contribution to the adsorption energy. Independent of the adsorption energies, the surface energies obtained by GGA and meta-GGA-based methods and their vdW corrections can still shed some more light on how the long-range vdW correlation works with the short-range exchange and correlation. Recently, Patra et al.75 demonstrated that the long-range vdW correction has a significant role in the surface energy. The surface energy is evaluated from the energies of the bulk and the slab as:

[TABLE]

with the information of the surface area A and n as the number of layers.

Our PBE, PBEsol, SCAN and SCAN+rVV10 results displayed in Table 4 agree with the surface energies found in ref.76, 77, 75 The results in Table 4 show that PBE seriously underestimates the surface energies in agreement with the results for jellium surface energies.52 This underestimation is well corrected by the rVV10 terms for Cu but not for Ag and Au. PBEsol, SCAN and revSCAN all yield excellent results for Cu, but not for Ag and Au. Addition of rVV10 corrections to PBEsol, SCAN and revSCAN worsen considerably the already excellent surface energies for Cu but yield better surface energies for Ag and Au. The best overall performer is SCAN+rVV10 followed by revSCAN+rVV10. These results suggest that, the Cu surface energy can be accurately calculated with a proper GGA or metaGGA without any rVV10 corrections, but accurate surface energies for Ag and Au require rVV10 corrections.

4 Conclusion

To describe the equilibrium adsorption of thiophene on coinage metal surfaces, we selected several promising density functionals such as SCAN, revSCAN, PBE and PBEsol. As GGA and meta-GGA functionals miss the long-range dispersion, we added rVV10 corrections as already suggested for SCAN. It is known that pairwise-interaction models such as rVV10 show incorrect asymptotic behavior, but as it was shown earlier this error is negligible around the equilibrium distance.

rVV10 correction requires one parameter b to adjust the long range part of interaction to the range of interaction present in the given parent functional. To obtain a fitted value we have chosen for reference the calculated CCSD(T) potential energy curve for Ar dimer. This curve is in reasonable agreement with the experiment and with higher level CCSDT(Q) curves around the well and the attractive branch of the curve. Our results for combining rVV10 to revSCAN, PBE, and PBEsol are b = 9.4, 9.8, and 9.7, respectively. With these parameters the revSCAN+rVV10 and PBEsol+rVV10 reproduce correctly the minimum energy around the correct 3.77 Å equilibrium distance, however the PBE+rVV10 yields minimum around 4 ∘. The order from best to worst fit is the following: PBEsol+rVV10, revSCAN+rVV10, SCAN+rVV10, and PBE+rVV10. In our calculations we apply the original b = 15.7 for SCAN+rVV10.

Our results for the lattice constants for Cu, Ag, and Au show that SCAN is not as accurate as PBEsol that yields the best lattice constants. The SCAN and revSCAN lattice constants agree well with each other, too short for Cu and too-long for Ag and Au. After the rVV10 correction the lattice constants shorten by 0.1-0.2 Å for Cu and by 0.2-0.5 Å for Ag and Au. The PBE+rVV10 lattice constants remain too-long and the PBEsol+rVV10 lattice constants become too short except for Au. The rVV10 corrected SCAN and revSCAN lattice constants are too short for Cu, excellent for Ag, and 0.1-0.2 Å too-long for Au.

We have shown that PBE seriously underestimates the surface energies of Cu, Ag, and Au, in agreement with the results obtained for jellium surface energies. This underestimation can be remedied by rVV10 correction for Cu but not for Ag and Au. PBEsol, SCAN and revSCAN yield excellent surface energies for Cu, but rVV10 corrections worsen the results.

Our results show that the lowest energy adsorption site on the coinage metal (111) surfaces is fcc-45 by all methods used in this paper, in agreement with experiment. For metal thiophene distances and adsorption energies we have chosen the PBE+vdW-dZK method for reference as it mimics the very accurate but computationally too demanding RPA binding energies perfectly for the interaction of graphene and metal surfaces. For the Cu-S distance revSCAN+rVV10, SCAN+rVV10 yield good agreement with our reference method and with the experiment.

According to our calculations SCAN and revSCAN underbind thiophene on Cu(111), Ag(111), and Au(111) surfaces by 0.1-0.3 eV. The rVV10 correction adds 0.3-0.5 eV to the binding energy making revSCAN+rVV10, SCAN+rVV10 overbinding by 0.2-0.4 eV. PBE+vdW-dZK and PBE+rVV10 show excellent agreement with estimated experimental results. PBEsol+rVV10 yields serious (0.4-0.6 eV) overbinding in accordance with the too short metal-S distance. Our calculations reflect the qualitative tendency that Cu and Au surfaces bind the thiophene about equally strongly and slightly stronger than Ag surface. This tendency is quantitatively reproduced by PBE+vdW-dZK, SCAN+rVV10, and PBE+rVV10.

We have demonstrated that good results of the rVV10 corrected density functionals for the well depth and the attractive region of the Ar dimer dissociation curve do not guarantee good results for thiophene adsorption on coinage metals. The order from best to worst fit is the following: PBEsol+rVV10, revSCAN+rVV10, SCAN+rVV10, and PBE+rVV10. However, the order of the performance for thiophene adsorption is the opposite with PBE+rVV10 being the best and PBEsol+rVV10 being the worst. We clearly show that the good fit to the Ar dimer curve does not guarantee good adsorption energies of polar molecules, e.g., thiophene on coinage metals.

The best method for thiophene adsorption is PBE+vdW-dZK, which is not only quantitatively correct for the adsorption energies but also predicts the ordering of adsorption energies for copper, gold and silver right along with the tilting angles and adsorption distances in good agreement with the experiment. {acknowledgement} We acknowledge the donors of the American Chemical Society Petroleum Research Fund for support of this research. We thank Center for Computational Design of Functional Layered Materials (CCDM) and High Performance Computing (HPC) at Temple University for providing the computational facilities. We thank John P. Perdew, Haowei Peng, Jefferson E. Bates, Hemanadhan Myneni, and Niraj K. Nepal for their valuable help during the calculations.

5 Supporting info:

**Redhead’s peak maximum method:

**Redhead’s peak maximum method (Readhead’s Analysis)37, 78 has been utilized to estimate the adsorption energy of the thiophene over Cu(111), Au(111) and Ag(111) surfaces from the thermal desorption spectra (TDS) . With the information of the peak maximum temperature ($T_{m}$), adsorption energy ($E_{ad}$) is calculated as:

[TABLE]

Where R is the ideal gas constant, $\nu$ is the pre exponential factor of desorption (commonly called prefactor) and $\beta$ is the heating rate. The value of $\nu$ is chosen 37 as $10^{13}s^{-1}$. The value of $T_{m}$ and $\beta$ are utilized from the TDS (TPD) experiments.59, 61, 66 Notice that the accuracy of Redhead’s method was questioned, but our analysis shows that despite its simplicity it gives reasonable adsorption energies for benzene and thiophene on the coinage metal surfaces. For example it was stated that ”Temperature-programed desorption (TPD) studies of bezene/Ag(111) yield a desorption temperature of 220 K. Depending on the value of the empirical prefactor in the Redhead equation utilized to convert this temperature into  EAds, values ranging from 0.43 eV to 0.80 eV have been reported in the literature [16–19].”24 Substitution $T_{m}$ = 220 K and $\beta$ = 1 K/s yields EAds = -0.60 eV (not -0.43 eV or -0.8 eV). Applying the more correct 230 K, the vanishing surface coverage limit, gives -0.63 eV. This value is within the error bar of the result of the complete analysis: -0.68 $\pm$ 0.05eV for Benzene on Ag(111) in the limit of vanishing surface coverage. It seems that $\nu$ = $10^{13}s^{-1}$ might yield reasonable values. The claim that ”Coverage dependence of the pre-exponential factor can introduce uncertainties of up to 0.2 eV for the molecules considered here.”38 is not supported by our results if the lowest published coverage leading to highest $T_{m}$ is considered and $\nu$ = $10^{13}s^{-1}$ in the Redhead’s adsorption energy estimation. The precision of the Redhead’s method is around or better than 0.1 eV for thiophene or benzene adsorption on coinage metals. Notice also that deviations of eq. (1) from the analytically correct expression are within 1.5% provided that $\nu$/$\beta$ falls between $10^{8}$ and $10^{13}K^{-1}$ for first order desorption.