Article pubs.acs.org/JPCA

Argon Interaction with Gold Surfaces: Ab Initio-Assisted Determination of Pair Ar−Au Potentials for Molecular Dynamics Simulations Romain Grenier,† Quy-Dong To,† María Pilar de Lara-Castells,‡ and Céline Léonard*,† †

Université Paris-Est, Laboratoire Modélisation et Simulation Multi Echelle, MSME UMR 8208 CNRS, 5 bd Descartes, 77454 Marne-la-Vallée, France ‡ Instituto de Física Fundamental (CSIC), Serrano 123, E-28006 Madrid, Spain S Supporting Information *

ABSTRACT: Global potentials for the interaction between the Ar atom and gold surfaces are investigated and Ar−Au pair potentials suitable for molecular dynamics simulations are derived. Using a periodic plane-wave representation of the electronic wave function, the nonlocal van-der-Waals vdW-DF2 and vdW-OptB86 approaches have been proved to describe better the interaction. These global interaction potentials have been decomposed to produce pair potentials. Then, the pair potentials have been compared with those derived by combining the dispersionless density functional dlDF for the repulsive part with an effective pairwise dispersion interaction. These repulsive potentials have been obtained from the decomposition of the repulsive interaction between the Ar atom and the Au2 and Au4 clusters and the dispersion coefficients have been evaluated by means of ab initio calculations on the Ar+Au2 complex using symmetry adapted perturbation theory. The pair potentials agree very well with those evaluated through periodic vdW-DF2 calculations. For benchmarking purposes, CCSD(T) calculations have also been performed for the ArAu and Ar+Au2 systems using large basis sets and extrapolations to the complete basis set limit. This work highlights that ab initio calculations using very small surface clusters can be used either as an independent cross-check to compare the performance of state-of-the-art vdW-corrected periodic DFT approaches or, directly, to calculate the pair potentials necessary in further molecular dynamics calculations.

1. INTRODUCTION In the context of nanotechnologies, the multiscale approaches based on molecular dynamics (MD) simulations have become essential. For instance, the huge surface/volume ratio of microporous media implies that surface effects must be taken into account when fluid flows are simulated in microporous media, such as slip velocity, temperature jump, adsorption, etc. In rarefied regimes, the continuum models are not suitable anymore and statistical molecular simulations must be employed. However, they are not suitable at scales larger than the nanometer. The high dependence of the process at small scale on the chemical nature of both the fluid and the surface is expected to provoke noticeable effects at larger scales. During the past few years, a new bottom-up multiscale multistep methodology has been developed.1−3 As a first step, the interaction potentials between the different atomic species composing the system must be determined or chosen from the literature. The second step involves MD simulations to extract accommodation parameters between the gas and the surface. These parameters are then transferred at a higher scale modeling step. Next, a diffuse-specular scattering kernel is used to model the interaction of the gas with the surface such that MD calculations can be carried out to simulate the complete flow. The velocities and density profiles can be © XXXX American Chemical Society

analyzed and compared with the extended Navier−Stokes equation solutions. The MD simulations are also able to provide the mean free path variations, which can be used to capture the velocity behavior in the zone near the wall (Knudsen layer) within the framework of the extended Navier− Stokes (NS) equations.4 As a solid surface material, it is well-known that gold has numerous technological relevant properties (it is easy to cover any surface with a thin layer, it has good conductivities, ...). As the interacting gas, rare gases are monatomic, almost nonreactive, and often used in nanotechnologies as transport gas. In particular, argon is a moderate weighted rare gas atom that has no the complexities related to lighter atoms such as helium, attaining the weakest van-der-Waals (vdW) interactions, or those of heavier atoms such as krypton, for which relativistic effects can be significant. Argon is also the third main constituent of the Earth’s atmosphere. Hence, the present study is dedicated to the interaction of an Ar atom with gold surfaces. The dimer Ar−Ar interaction potential is well-known with high accuracy.5 Numerous gold solid potentials exist in the Received: April 20, 2015 Revised: June 4, 2015

A

DOI: 10.1021/acs.jpca.5b03769 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A literature, as pairwise Lennard-Jones-type potentials,6 or effective many body quantum Sutton−Chen,7 EAM,8,9 and MEAM10 potentials, for which the local electron density surrounding the metal atoms is taken into account. Most of them are empirical and have been fitted to reproduce both bulk and surfaces properties (surface energy, melting temperature, phonon vibrations, ...). As often occurs with many other gassurface interacting systems, the Ar−Au(111) potential has rarely been investigated. To the best of our knowledge, the only reported works are those by Chizmeshya and Zaremba11,12 and the refitted potential of Tang and Toennies.13 These theoretical predictions for the adsorption energy of the Ar atom on the Au(111) surface are based on a scattering theory formulation for He-surface metal interactions, using the atomic Hartree− Fock (HF) approximation to model the repulsive interaction. The attractive part was determined by Zaremba and Kohn from second-order perturbation theory.14 These approaches have been validated for the interaction of a He atom with metallic surfaces for which scattering-based experimental data exist. However, a global interaction potential of a gas atom with a surface is not adapted for MD simulations in which the temperature-induced atomic displacements has to be accounted for. An Au−Ar pair potential is instead needed. The question is how to obtain the pair potential for an Ar atom interacting with a gold surface from ab initio calculations, because this pair potential is certainly different from the Ar−Au diatomic potential. In fact, the interaction of the Ar atom with a metallic surface can be associated with the Ar interaction with its image charge in the metal.15 The attractive part can be approximate by V (Z ) = −

⎛1 ⎞ + o⎜ 5 ⎟ ⎝Z ⎠ (Z − Z0) C3

3

or

−

on graphene,18 the determination of the He-TiO2(110) potential energy surface to model soft, 4He droplet-assisted, deposition processes,19 and the interaction of the silver dimer with graphene.20 In the present work, this approach has been applied to the Ar atom interaction with small gold clusters (Au2 and Au4). The application of the two approaches will help to validate the results because they are grounded differently. It must be noticed that the ab initio schemes developed by de Lara-Castells et al.17−20 were proposed to calculate global interaction adsorbate−surface potentials either than a direct determination of pair potentials to be used in further MD calculations, which is the present purpose. The paper is organized as follows. In the first section, the Ar−Au diatomic potential is computed using highly accurate ab initio methods and the resulting benchmark potential is compared with those computed with state-of-the-art periodic DFT approaches. The second section is devoted to the test of the considered density functionals on the determination of the Ar/gold surfaces interaction, and the third one to the derivation of pair potentials using the two approaches described previously. These pair potentials have then been used to reconstruct a global interaction potential. In this work, the distances and the energies are given in angström (1 Å = 10−10 m) and electronvolt (1 eV = 1.602176565(35) × 10−19 m2 kg s−2) units, respectively.

2. THE ISOLATED AR + AU DIATOMIC SYSTEM As a first test to choose the adequate vdW-corrected DFT treatment for the extended system, we have applied different plane-wave-based DFT approaches to calculate potential energy curves for the Ar + Au interaction, comparing them with that computed through high-level ab initio calculations. 2.1. Computational Details. First, electronic calculations on the Au + Ar diatomic system have been performed through coupled-cluster singles, doubles, and noniterative triples [CCSD(T)] calculations by means of the MOLPRO code.21,22 The spin-restricted open-shell coupled cluster, RCCSD(T),23,24 method has been employed using the augcc-pwCVnZ and aug-cc-pwCVnZ-PP basis sets, for Ar and Au respectively, with n = 3, 4, 5.25,26 The inner electrons of the Au atoms have been replaced by a (60-electron) relativistic pseudopotential, ECP60MDF,27 whereas the remaining 19 electrons, namely, the electrons occupying the 5s, 5p, 5d, and 6s atomic orbitals, have been described by the associated valence basis sets. The 1s, 2s and 2p orbitals of the Ar atom have been frozen in a core space. Then, the remaining 27 electrons have been correlated in the 5s, 5p, 5d, 6s orbitals of the Au atom, as well as the Ar 3s and 3p valence orbitals. All interaction energies were counterpoise (CP)-corrected,28 which is crucial for weakly bound complexes such as vdW complexes. The complete basis set (CBS) limit of the interaction energies was estimated by applying the scheme of Helgaker and coworkers for the CCSD(T) correlation energies29,30

⎛ 1 ⎞ + o⎜ 4 ⎟ ⎝Z ⎠ Z

C3

3

(1)

where Z is the distance of the Ar atom to the surface plane. It has been proved that π C3 ≈ nsC6 (2) 6 where ns denotes the homogeneous substrate atom density16 and C6 is the vdW parameter such that the attractive part of the interaction of Ar with a single Au surface atom can be approximated as v (r ) = −

C6 r

6

−

C8 r

8

⎛ 1 ⎞ + o⎜ 10 ⎟ ⎝r ⎠

(3)

where r is the distance between Ar and Au atoms. The present work focuses on the ab initio determination of the interaction potential between an Ar atom with a gold surface to extract the Ar−Au pair potential suitable for MD simulations. Two different computational strategies have been followed and confronted: (i) the application of periodic DFT approaches based on the use of plane waves for the description of the electronic wave function. This treatment has the advantage of accounting for the electron delocalization on the metal surface. Because the choice of both the DFT functional and the van-der-Waals correction scheme are crucial, several vdW-corrected DFT approaches have been tested: (ii) the use of dispersionless density functional theory coupled with an ab initio-based determination of dispersion coefficients. The second approach has been successfully applied to the interaction of a He atom with clusters modeling the TiO2(110) surface,17 the 4He nanodroplet collision dynamics

E(n) = ECBS + A n−3

(4)

with n = 3, 4, 5, using aug-cc-pwCVnZ and aug-cc-pwCVnZ-PP basis sets. The HF interaction energies were fixed to those obtained with the largest basis. The resulting potential energy curve is plotted in Figure 1. DFT calculations have been performed using the VASP code31,32 with the projector augmented-wave (PAW) method33,34 for the electron−ion interactions. The generalized B

DOI: 10.1021/acs.jpca.5b03769 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

estimates the lattice constant, Klimeš et al.40 have modified the nonlocal vdW-DF treatment by replacing the revPBE exchange functional approach with the reoptimized version of the B86b functional, leading to the vdW-OptB86b approach. To perform the plane-wave-based calculations for the diatomic, a cube with a side length of 46 Å was employed to avoid interdiatomic interactions. Only the Γ point was sampled. An energy cutoff of 266 eV and a Gaussian smearing of width σ = 0.001 eV have been employed in all calculations, which were spin-polarized. Using these numerical parameters, the interaction energy was found to be converged to better than 1 meV. Interaction energies have been computed for atomic distances between 3.0 and 10.0 Å and have been fitted to the following analytical expression: vdiat(r ) = v0e−αr −

Figure 1. Potential energy curves of the ArAu diatomic in the X2Σ+ electronic ground state.

r

f (β r ) −

6 6

C8 r8

f8 (βr )

(6)

where r is the Ar−Au internuclear distance and ⎛ (βr )n ⎞ fn (r ) = 1 − e−βr ⎜1 + r + ... + ⎟ ⎝ n! ⎠

gradient approximation (GGA) has been used for the exchange−correlation functional with the parametrization of Perdew, Burke, and Ernzerhof (PBE)35 as well as the revised versions revPBE36 and PBEsol37 of the same functional, with the latter optimized for the study of solids. Being gradientdependent, these semilocal density functionals are designed to include short- and medium-range correlation effects only. The long-range dispersion corrections have been taken into account by two approaches: (i) the DFT-D3 approach of Grimme38,39 and (ii) the nonlocal van der Waals density functionals (vdWDF) vdW-OptB86b40 and vdW-DF2.41 In contrast to the empirical underpinnings of the DFT-D3 approach, the vdW-DF method includes long-range vdW corrections in a nonempirical way Exc = ExGGA + EcLDA + Ecnl

C6

(7)

are damping functions. The analytical function of eq 6 is a popular Buckingham-type potential, with an exponential Born− Mayer component.42 To reduce the divergent character of the attractive terms at short interatomic distances, these terms are multiplied by the Tang−Toennies damping functions,43 f n(r), which are close to 1 for large r but continuously go to zero when r decreases. The resulting parameters have been compiled in Table 1 and compared with experimental data.44,45 2.2. Results. Previous experimental studies on the ArAu diatomic indicate that the potential well depth is about D0 = 18.5 ± 2 meV. This value including the zero-point energy, De, should equal 20.1 ± 2 meV if the harmonic vibrational energy of 3.26 meV, from RCCSD(T) CBS calculations, is considered. This value is in good agreement with both the RCCSD(T) CBS value and that reported by Gardner et al. (22.91 meV).46 Table 1 and Figure 1 show that the DFT-based results are mainly out of the range from the RCCSD(T) CBS counterpart. The revPBE and PBE well depths are more than 8 meV less attractive whereas De(PBEsol) is too large by more than 7 meV, with the equilibrium distance re being ∼0.5 Å too short. The well depths of the vdW-corrected PBE-D3 and vdW-DF2 potentials are more than 10 meV too attractive and slightly shifted to shorter equilibrium distances. Only the vdWOptB86b approach produces De and re values in close agreement with the RCCSD(T) CBS benchmarking. However, the medium- and long-range attractive tail of the vdW-OptB86b

(5)

Enlc ,

where the nonlocal correlation energy functional, is based on electron densities interacting via a model dielectric function and approximately accounts for dispersion interaction. Lee et al.41 have demonstrated that the performance of the vdW-DF exchange−correlation functional can depend sensitively upon the choice of the exchange functional. They have replaced the revPBE exchange functional of the original vdW-DF treatment by a revised version of the PW86 functional, giving rise to the vdW-DF2 functional. The long-range part of the nonlocal correlation functional, Enlc , has also been improved by applying a large-N asymptote gradient correction in determining the vdW kernel. On the contrary, as the vdW-DF2 treatment over-

Table 1. Parameters for the ArAu (X2Σ+) Interaction Potential Calculated by Different Methods and Fitted to Eq 6

a

method

De (meV)

req (Å)

v0 (×106 meV)

α (Å−1)

C6 (meV Å6)

C8 (×106 meV Å8)

β (Å−1)

RMS (meV)

RCCSD(T) CBS revPBE PBE PBEsol PBE-D3 vdW-OptB86b vdW-DF2 VMI exp44 REMPI exp45 RCCSD(Val,T)/d-awCV∞Z46

22.639 6.875 13.945 29.840 43.582 23.232 34.413 18.5 ± 2a 16.12 ± 2a 22.91

3.711 4.501 3.781 3.249 3.601 3.781 3.650

1.7493 0.5870 1.5826 2.2964 3.1075 0.3628 3.1075

2.815 2.426 2.840 3.224 2.769 2.360 2.850

67008 99731 20213 9803 42320 204816 38099

1.7374 0.9286 1.7429 1.2044 4.8569 0.6549 3.6044

4.660 5.625 4.690 3.958 4.656 4.733 4.566

0.008 0.459 0.210 0.342 0.458 0.391 0.368

3.711

D0 values C

DOI: 10.1021/acs.jpca.5b03769 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Table 2. Optimized Cubic Cell and Potential Parameters for the Ar + Au(111) Global Interaction Potential Calculated by Means of Different DFT Methods Using a 2 × 2 3-Layer Supercell Model, Considering the Atop Adsorption Site, and Fitted by the Expression of Eq 8

a

method

a (Å)

De (meV)

Zeq (Å)

V0 (×106meV)

γ (Å−1)

C3 (meV Å3)

Z0 (Å)

RMS (meV)

PBE PBEsol PBE-D3 vdW-OptB86b vdW-DF2 theory13 theory11 theory12 theory14 theory47

4.18 4.10 4.11 4.13 4.35a

10.24 13.61 135.68 126.87 100.36 90.9 100 85.1

4.312 3.800 3.518 3.549 3.685 2.19 1.91 2.08

1.176 13.16 3.541 0.864 5.001

2.84 3.67 2.87 2.57 2.92

422 142 2309 4598 1659

1.325 2.021 1.501 0.794 1.679

0.50 0.13 0.63 0.90 1.35

1811 1768 1847

0.145 0.189

The value of the cell parameter used in the present interaction potential calculations is 4.08 Å.48

test has been performed on the Au(111) surface. This surface has been modeled using a repeated slab of 3 atomic layers. A 2 × 2 supercell and vacuum spacing of 25 Å have been used to minimize the interactions between the periodic images along the three directions of the slab. An energy cutoff of 266 eV has been set for all the calculations, employing a 8 × 8 × 1 Monkhorst−Pack integration grid. Using the Methfessel− Paxton smearing of first order with a width of 0.01 eV for the Au atoms, the interaction energy was found to be converged to within 0.3 meV. All calculations were spinpolarized and the bulk cell parameter was first optimized for each density functional. Next, a single Ar atom approaching at the atop adsorption surface site was considered. Structural relaxation effects of the gold atom positions were not included after testing their small influence. Ar−surface potential energy curves have been calculated for vertical height values above the surface (referred to as Z) between 2.8 and 12.0 Å, having been fitted by the following analytical expression:

potential decays much more slowly than that of the RCCSD(T) CBS potential (Figure 1), as reflected in the large differences between the C6 and C8 values (Table 1). As can be seen in Table 1, the values of the PBE-D3 and vdW-DF2 C 6 coefficients are the closest ones to the corresponding RCCSD(T) CBS value. As clearly apparent in Figure 1, the long-range vdW corrections of vdW-OptB86b and vdW-DF2 approaches are rather different: the vdW-DF2 correction is shorter ranged. As mentioned in ref 20, this might be due to the fact that the gradient correction factor for the nonlocal Enlc functional is more than twice as large, following atomic scaling laws. As will be discussed later, this is also found for the global Au + surface interaction potential as well as the pair potential to be used in MD simulations. At this point, it should also be stressed that the validation of DFT-based approaches on the Au + Ar diatomic system should be taken carefully. The long-ranged correction is mainly determined by dispersion C6 and C8 coefficients, which are expected to have reasonable transferable properties between molecular models of the extended system. Hence, the different long-range attractive tails of the vdWOptB86b and vdW-DF2 potentials might be reflected on the extended Ar + surface system. However, this is not necessarily true at the regime where the electronic densities of the two monomers start to overlap. This is clearly evident in Figure 1 by the failure of the PBEsol approach to provide a good potential energy curve. In fact, the parameters of the PBEsol functional were optimized to provide accurate values of solid lattice constants, at the slowing varying density regime, compromising then the performance of the density functional on molecular systems. The revPBE, PBE, and PBEsol sequence in Figure 1 emphasizes the enormous dependence on the exchange functional at the repulsive potential region, which becomes smoothed out when vdW-DF approaches are used. This is another clear argument to retain only vdW-DF approaches for the determination of the interaction and pair potentials. Another reason limiting the value of the Au + Ar complex to provide the pair potential is its open-shell character. As discussed below, the closed-shell Au2 + Ar complex is better used for that purpose.

Vsurf (Z) = V0e−γZ −

C3 (Z − Z0)3

(8)

where Z is the distance between the Ar atom and the Au(111) atop atom position. Table 2 collects the resulting parameters, and the interaction potentials are shown in Figure 2.

3. THE AR/GOLD SURFACE 3.1. Structural Models and Computational Details. Aimed to evaluate the performance of the considered DFTbased treatments on the determination of the interaction potential of an Ar atom with gold surfaces, a first comparative

Figure 2. Ar + Au(111) global interaction potentials calculated by different DFT methods using a 2 × 2 3-layer supercell, and considering the atop adsorption site. D

DOI: 10.1021/acs.jpca.5b03769 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A 3.2. Results. The results of this work have been compared with those of pioneering studies in which the C3 vdW coefficients were determined in terms of the dynamic polarizability of the adatom and the dynamic linear dielectric function of the substrate. In the work of Bruch,47 C3 parameters were evaluated for several adatom−substrate systems by using experimental response values. On the contrary, Chizmeshya, Zaremba, and Kohn11,12,14 developed a theory to determine physisorption potentials from experimental data. Specifically, applying the jelium model and perturbation theory, they proposed a model to determine the C3 parameter within a term accounting for the position of a reference plane [i.e., Z0 of eq 1]. The most elaborate version of the model rare gas−metal interaction potentials used scattering theory for the description of the repulsive potential regions. Next, Tang and Toennies13 improved the physisorption potential model proposed by Chizmeshya and Zaremba using damping functions to better connect the repulsive and the attractive potential regions. The improved model was capable of predicting a potential minimum for the Xe atom on gold surfaces. From Figure 2 and Table 2, it is clear that the vdWuncorrected PBE and PBEsol functionals are not suitable for describing the Ar−Au(111) interaction, as foreshadowed by the results for the diatomic. The PBE-D3 and vdW-DF-based treatments produce parameters in good agreement with the previous results of Chizmeshya and Zaremba,11,12 Zaremba and Kohn,14 Bruch,47 and Tang and Toennies13 with the exception of the geometrical parameters Zeq and Z0, which are more than 1 Å larger. However, notice that the values of the Zeq−Z0 differences are comparable. The vdW-DF2 approach provides the closest parameters to the previous data, especially for the well depth and the dispersion energy coefficient C3, but it fails to accurately reproduce the gold bulk lattice constant. Hence, the experimental value of a = 4.0786 Å from X-ray diffraction has been used48 when the vdW-DF2 approach is applied to the Ar−surface system. The comparative analysis of Chen et al.49 showed that the vdW-DF2 functional provides equilibrium distances, adsorption energies, and vibrational energies in good agreement with the experimental data for several rare gas + metal systems, without including the Ar + gold surface pair. The vdW-OptB86b treatment leads to a more reasonable a value, but the agreement of the interaction potential parameters with the reported values is less good. The PBE-D3 approach leads to the deepest potential minima. However, the C3 dispersion coefficient and geometrical parameters obtained with the PBE-D3 treatment are rather close to the vdW-DF2 ones. Jordan and collaborators have also demonstrated the good performance of state-of-the-art vdW-corrected DFTbased methods, such as the nonlocal optPBE-vdW approach, to treat dispersion effects on surfaces.50 On the basis of this analysis, the vdW-DF2 and vdW-OptB86b approaches have been chosen for the determination of the MD pair potentials. The corrugation of the interaction potential has also been examined for the Au(111) and Au(100) surfaces using the vdW-DF2 and vdW-OptB86b approaches. To this end, a larger 3 × 3 4-layer supercell model has been used. Four and three adsorption sites have been considered for the Ar atom on the Au(111) [referred to as hcp, fcc, bridge, and atop] and the Au(100) surface [denoted as hollow, bridge, and atop], respectively (Figure S1, Supporting Information). Using the vdW-OptB86b and vdW-DF2 approaches, the adsorption site dependence of the interaction potential well depths, De, is illustrated in Figures 3 and 4 for the Au(111) and Au(100)

Figure 3. Corrugation of the interaction potential upon the adsorption of a single Ar atom on the Au(111) surface (considering a 3 × 3 supercell with 4 layers). The reference, for each curve, is the potential well depth corresponding to the hcp adsorption site, i.e., De = 130.455 meV (vdW-OptB86b) and De = 101.318 meV (vdW-DF2), respectively. For the values resulting from the recomposition procedure, the reference is the adsorption energy on the hcp site of the Au(111) surface, as given in Table 4.

Figure 4. Corrugation of the interaction potential upon the adsorption of a single Ar atom on the Au(100) surface (using a 3 × 3 supercell model with 4 layers). The reference, for each curve, is the potential well depth corresponding to the hcp adsorption site on the Au(111) surface, i.e., De = 130.455 meV (vdW-OptB86b) and De = 101.318 meV (vdW-DF2), respectively. For the values resulting from the recomposition procedure, the reference is the adsorption energy on the hcp site of the Au(111) surface, as given in Table 4.

surfaces, respectively. These figures show the adsorption energy differences with respect to that attained on the (most stable) hcp site of the Au(111) surface. Using both nonlocal vdW-DF approaches, the adsorption energy differences on the Au(111) and the Au(100) surfaces are smaller than 2 meV. As opposed to previous theoretical studies characterizing the atop position as the most stable site for rare gases on metal surfaces,49,51 the hollow site is found to be more attractive for the Ar atom adsorbed on both Au(111) and Au(100) surfaces. As shown in ref 49, however, the employment of nonlocal vdW functionals reduces strongly the interaction potential corrugation as compared with the LDA or DFT-D2 approaches. As can be seen in Figure 3, the potential corrugation on the Au(111) surface is smaller than 2 meV whereas that on the Au(100) surface is larger, as expected from the more pronounced roughness of the Au(100) surface as compared with the Au(111) one. The difference between the atop and hollow Au(100) adsorption energies reaches 7 meV using the vdWOptB86b approach. However, this difference corresponds to just ≈5% of the total adsorption energy. The application of the E

DOI: 10.1021/acs.jpca.5b03769 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A method of increments52 at CCSD(T) level on the Xe +Mg(0001) system highlighted that the perturbation of the approaching Xe atom is more effectively screened when adsorbed on the atop position. Again, the differences between atop and hollow adsorption energies were found to be very small (5 meV and ∼6% of the total interaction energy). These differences were reflected in the appearance of attractive intramonomer correlation contributions for the atop position, somehow compensating the Pauli repulsion, rather than in favored intermonomer (dispersion-like) correlation effects.

with Z defined as in eq 10. Using this functional form, only the vdW-OptB86b and vdW-DF2 global interaction potentials associated with the Ar atom approaching the hcp surface site have been decomposed. We have tested that the inclusion of additional C10/r10 terms in eq 12 is irrelevant, reducing the root-mean-square of the fitting by less than 0.1 meV. In fact, the values of the optimized C10 coefficients turned out to be negligible (about 0.1 meV Å10), modifying the well depth of the fitted pair potentials by less than 0.01 meV. 4.1.2. From Ar + Au2 and Ar + Au4 Clusters. The pair potentials have also been obtained by following an ab initiobased strategy developed by de Lara-Castells and collaborators to calculate accurate vdW-dominated adsorbate−surface global interactions.17−20 It extends the dispersionless density functional dlDF approach developed by Szalewicz and collaborators53 to adsorbate−surface interactions by including periodic boundary conditions.19 Szalewicz and collaborators also developed a general-purpose (effective interatomic pairwise) dispersion functional named Das to be used along with the dlDF approach.53 The general-purpose Das parameters were extracted by Pernal et al.53 from the dispersion energies calculated on weakly bound complexes by means of symmetry adapted perturbation theory with a DFT-based description of the monomers (i.e., the SAPT(DFT) approach). As proposed by de Lara-Castells et al.,17−20 the dlDF+Das approach can also be applied by determining the dispersion Das parameters on surface cluster models up to convergence. The proposed strategy also introduces a scheme to parametrize the dispersion interaction (so-called incremental Das*) by calculating two- and three-body dispersion terms at CCSD(T) level via the method of increments.19,20 The efficiency of the periodic dlDF + incremental D*as approach relies on the transferability properties of the dispersion coefficients upon augmenting the surface cluster model as well as the short-ranged nature of the dispersionless correlation contribution. Previous works have considered TiO2(110), graphene, and graphite substrates whereas a metal surface is considered here for the first time. Because pair potentials are necessary instead of global interaction potentials, the procedure can be simplified as compared with that described in previous studies.17−20 The pair potential can be decomposed as a sum of repulsive plus attractive terms. As already mentioned in the Introduction, the attractive part of the Ar + Au potential interaction is mainly due to the dispersion whereas the dispersionless interaction is repulsive at short range. As extensively discussed in previous works,17−19 the dispersionless correlation contribution is repulsive at short-range due to so-called correlation (virtual orbital) space truncation effect exerted by the noble gas atom on the localized orbitals centered on the surface atoms. The benchmarking via CCSD(T) calculations has previously shown that the dlDF approach53 provides accurate estimations of this correlation energy contribution19 and, more generally, all the interaction energy components different than the dispersion. The dlDF functional differs from the Minnesota M05-2X functional54 in the number and values of DFT parameters that were optimized to reproduce benchmark dispersionless interaction energies of weakly bound dimers.53 These benchmark dispersionless interaction energies were obtained by removing from total CCSD(T) interaction energies the dispersion contributions, as evaluated with the SAPT(DFT) approach for a training set of molecular systems.55,56 The dlDF calculations have been performed for two complexes, Au2 + Ar (Figure S2, Supporting Information)

4. DETERMINATION OF THE AR−AU PAIR POTENTIAL 4.1. Computational Details. 4.1.1. From Ar/Gold Surface Global Interaction Potential. Once the global interaction potentials have been computed and validated by means of a comparative analysis with previous experimental works, the pair potentials can be extracted for MD simulations. By following the strategy described in Introduction, we decomposed the interaction potential, Vsurf, into a sum of pair potentials vpair, natcluster

Vsurf ({ri};Z) ≈

natcluster

∑

vpair[ri(Z)] =

∑

i

v0e−αri −

i

C6 ri 6

−

C8 ri 8 (9)

in which Z corresponds to the vertical height of the Ar atom above the basal plane of the gold substrate. This plane contains the gold atoms located at the upper layer of the slab. The term ri is defined as the distance between the Ar atom and the ith Au atom with a parametric dependence on Z: ri(Z) =

xi 2 + yi 2 + (Z + zi)2

(10)

with {xi, yi, zi} as the Cartesian coordinates of the same atom. In what follows, the Z parametric dependence of the ri distances will be omitted. The v0, α, and C6 parameters have been adjusted so that the total Vsurf values fit the vdW-OptB86b or vdW-DF2 interaction energies. We have shown in the previous section that the potential corrugation on the Au(111) and Au(100) surfaces can be ignored and that these two surfaces bear comparable adsorption energies. Previous studies have revealed that the Au(111) surface is the most stable.6−10 On the basis of these considerations, the vpair suitable for MD simulations has been determined as an average of the pair potentials associated with the main adsorption sites of the Au(111) surface. Specifically, the average potential was obtained from these curves using geometrical (position density counting) weighting in the unit cell, vave =

vatop + v hcp + vfcc + 2v bridge (11)

5

A 147 atom 3-layer Au(111) cluster has been chosen for the decomposition, after ensuring that it provides converged vpair potential energy curves to within 0.1 meV, i.e., less than 1% of the vpair well depths. We have also tested the performance of the pair potential function used to fit the diatomic potential (eq 6): nat cluster

Vsurf ({ri};Z) ≈

∑

vpair(ri)

i nat cluster

=

∑ i

v0e−αri −

C6 ri 6

f6 (βri) −

C8 ri 8

f8 (βri) (12) F

DOI: 10.1021/acs.jpca.5b03769 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

providing the dlDF(Au2) + Das and dlDF(Au4) + Das pair potentials. The Das parametrization has also been used to estimate dispersion energies for the Au4 + Ar system. These estimated dispersion energies agreed to within 2% with those directly calculated using the SAPT-PBE0 approach. Using the Das parameters for the Ar atom reported by Pernal et al.,53 we obtained C6(Au) = 192.73 eV Å6, C8(Au) = 1619.4 eV Å8, and β(Au) = 2.681 Å −1. The C6 value compares very well with that calculated by Hatz et al.59 (215.53(33) eV Å6) from the isotropic part of the Au2−Au2 interaction. It is also consistent with that extracted by Tonigold and Groß60 (220.2 eV Å6) from the adsorption energies of small aromatic molecules on the Au(111) surface, using density functionals aided with semiempirical corrections for dispersion effects. However, our calculated C6 value clearly disagrees with that reported by Grimme38 (615 eV Å6) for the D2 parametrization. For benchmarking purposes, CCSD(T)/AV5Z calculations have also been performed for the Au2 + Ar system, with the resulting interaction potentials having been decomposed using eq 12. 4.2. Results. The pair potentials are shown in Figure 5, and the resulting fitting pair parameters are collected in Table 3.

and Au4 + Ar (Figure S3, Supporting Information), using the aug-cc-pVnZ (Ar) and aug-cc-pVnZ-PP (Au) basis sets with n = 3 and 5, respectively.26,27,57 The aug-cc-pV5Z (Ar) and aug-ccpV5Z-PP (Au) basis sets will be collectively referred to as AV5Z. For the Au2 + Ar complex, the Au2 bond length has been kept fixed to the equilibrium geometry and the Ar−Au2 Tshaped configuration was considered, varying the distance, Z, between the Ar atom and the Au2 bond midpoint. For the Au4 + Ar complex, the gold cluster was modeled as a tetrahedron and different distances, Z, from the Ar atom to the midpoint of its triangular base were considered (Figure S3, Supporting Information). For both clusters, the dlDF energies have been decomposed and fitted to pair repulsive potentials, as described above (eq 9): natcluster

Vcluster({ri};Z) ≈

∑

natcluster

vpair(ri) =

∑

i

v0e−αri

i

(13)

with Z defined as in eq 10. Next, the dispersion has been included in the pair potentials via the effective interatomic pairwise functional form Das as follows, nat cluster

∑

Das =

i

−

−

C6ArC6Au ri

6

C8ArC8Au f8 ( ri 8

f6 ( β Ar β Au ri)

β Ar β Au ri)

(14)

where ri is the internuclear distance between the Ar atom and the ith Au surface atom, as defined above. This functional form has been shown to be successful in fitting the benchmark intermonomer correlation (dispersion-like) contributions calculated at CCSD(T) level via the application of the method of increments.17−20 The effective reduction in the two-body dispersion terms of eq 14 due to the inclusion of three-body incremental terms was also accounted for,18,20 with the damping and C8 parameters playing a role. Screening effects were also analyzed for the Ag2/graphene interaction (see Table 2 of ref 20), showing that they give rise to an effective reduction of two-body incremental terms of about 1%. For the He/ TiO2(110) interaction,17 it was found that additional C10/r10 are necessary to account for the interaction from sublayer surface atoms if the damping factors are omitted. These studies reflect that damping functions might be useful in effectively accounting for screening effects. Also, the total (two- plus three-body) intermonomer correlation contributions were found to agree well with those using DFT-based SAPT.18,20 The Das parameters, (C6ArC6Au )1/2 = 76785 meV Å6, Au 1/2 (CAr = 1.0669 × 106 meV Å8, and (βArβAu)1/2 = 3.051 8 C8 ) −1 (2) Å , were obtained by fitting SAPT E(2) disp + Ex‑disp energies. The SAPT calculations were performed on the Ar + Au2 system using the PBE0 functional58 for the monomers and the AV5Z basis set. The SAPT-PBE0 approach was chosen to extract the dispersion energies on the basis of previous benchmarking on the He/TiO2(110) interaction,17 revealing a very good agreement with those determined at CCSD(T) level. We have checked that the extrapolation to the CBS limit does not improve the values of the present AV5Z Das parameters. For example, the estimated Das C6 coefficient at the CBS limit differs by less than 0.1% from the AV5Z value. The calculated Das parameters were combined with the dlDF repulsive pair potentials from both the Au2 + Ar and Au4 + Ar interactions,

Figure 5. Pair potentials obtained from dlDF+Das calculations for the Au2 + Ar and Au4 + Ar complexes, from CCSD(T) calculations for the Au2 + Ar cluster, and from vdW-OptB86b and vdW-DF2 calculations on the Au(111) + Ar extended system (see also Table 3). The dispersion energies determined with the SAPT-PBE0 approach are also shown for comparison purposes.

The main conclusions drawn from the comparison of the pair potentials are (i) the vdW-DF2 and dlDF+Das approaches provide very similar pair potentials despite the very different backgrounds of the two treatments, which is the most important finding of the present work. In this way, the C6 and C8 parameters of the Das functional lie in the range of the nonlocal vdW-DF counterparts. In particular, the value of the C6 Das coefficient is between the vdW-DF2 and vdW-OptB86b C6 values. Actually, when comparing the C6 and C8 coefficients obtained from the Das parametrization and those evaluated from the fitting of total interaction energies (Table 3), one should take into account that the former ones encompass dispersion effects only, whereas the latter ones might reflect other attractive interaction energy contributions that, within the dlDF approach, could be embodied as an effective reduction of the repulsive interaction terms. (ii) The convergence of the repulsive pair parameters from dlDF calculations is almost G

DOI: 10.1021/acs.jpca.5b03769 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A Table 3. Pair Potential Parameters Corresponding to the Analytical Form of Eq 6

a

method

De (meV)

req (Å)

v0 (×106 meV)

α (Å−1)

C6 (meV Å6)

C8 (×106 meV Å8)

β (Å−1)

hcp vdW-OptB86b av vdW-OptB86b hcp vdW-DF2 av vdW-DF2 CCSD(T)/AV5Z PBE0 SAPT/AV5Z dlDF(Au2)+Das dlDF(Au4)+Das

15.435 13.792 12.317 12.605 13.606 18.513 11.360 10.335

4.0911 4.2190 4.2237 4.2010 4.1014 3.9030 4.2848 4.3014

1.3239

2.569

99343

1.8447

4.444

1.651

3.5925

2.864

44999

2.4813

4.864

1.745

2.8611 2.9715 16.230a 5.7663a

2.867 2.933 3.356a 3.090a

58927 62185 76785 76785

1.8879 1.7786 1.0669 1.0669

5.820 4.478 3.051 3.051

0.011 (0.7)b 0.067 (1.0)b 0.969a (0.8)b 0.830a

RMS (meV)

Only the repulsive part has been fitted. bEstimation of the basis set incompleteness effect at the minimum geometry in meV.

already reached for the Au2 cluster, with the dlDF calculations on the Ar + Au4 complex providing very similar results. The pair potential extracted from the Ar + Au4 calculations is slightly more repulsive at the minimum (Figure 5). Notice, however, that the energy difference (about 2 meV) is twice as large as the root-mean-square (RMS) of the fitting (about 1 meV, Table 3). The energy differences between vdW-DF2 and dlDF+Das interaction energies are also of the same order of magnitude than the RMS of the fitting from the vdW-DF2 interaction energies (about 2 meV, Table 3). The errors introduced by the basis set incompleteness should be also considered (to within 1 meV at the minimum, Table 3) as well as those of the VASP calculations (less than 1 meV, see above). (iii) As compared with both the CCSD(T)/AV5Z and the SAPT-PBE0 benchmarking, the long-range tail provided by the vdW-OptB86b treatment decays too slowly. The energies of the nuclear bound states supported by the different pair potentials can also guide the comparative analysis because they are determined by a balance between the repulsive and attractive potential regions. The vibrational energy levels of the pair potentials have been computed up to v = 9 using the Cooley algorithm by means of the NUMEROV code.61 Figure 6 shows the nuclear bound state energies as a function of the

vibrational energies are the closest ones, confirming that the corresponding pair potentials are in very good agreement, as already noticed from the results presented in Table 3 and Figure 5. Aimed to examine the quality of the calculated pair potentials, the global gas/surface interaction potentials have been reconstructed by considering that the Ar atom approaches the Au(111) surface in front of the hcp adsorption site. The cluster model of the semi-infinite surface was composed of 3 layers and 147 atoms, having checked that the addition of a fourth layer leaves the recomposed potential almost unmodified. The resulting potential energy curves have been plotted in Figure 7. The interaction energies directly calculated using the

Figure 7. Recomposed Ar−Au(111) global interaction potentials on the hcp site, using a 147 atom 3-layer cluster and considering the dlDF +Das and vdW-DF pair potentials (Table 4). The markers correspond to VASP calculation results.

vdW functionals through VASP calculations have also been represented to allow a comparison with the reconstructed potential curves. To facilitate the comparative analysis, the interaction energies have been fitted by eq 8. The fitting parameters have been collected in Table 4 such that they can be directly compared with the results presented in Table 2. The Vsurf parameters obtained before the decomposition and after the recomposition have been compared, as well as those resulting from the recomposition potential calculated from the average pair potentials. The results presented in Table 4 and Figure 7 indicate that the qualities of the dlDF+Das and vdWDF2 pair potentials are good. Remarkably, the recomposed curves obtained from both dlDF+Das pair potentials are very similar and close to the vdW-DF2 recomposed potential energy curves. The fitting parameters from the energies directly

Figure 6. Vibrational energy levels (meV) supported by the calculated pair potentials (see also Table 3).

vibrational quantum number v. The vibrational energies are in millielectronvolts, and the energy reference corresponds to the dissociation energy of the potential. Notice that the two vdWDF2-based approaches provide almost identical vibrational energies. All the vibrational energy profiles bear a parallel behavior as a function of v up to v = 3, with the exception of the dlDF(Au4)+Das profile. The vdW-DF2 and dlDF(Au2)+Das H

DOI: 10.1021/acs.jpca.5b03769 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Table 4. Recomposed Interaction Potentials Considering the hcp Adsorption Site and a Gold Cluster Model of the Au(111) Surface (Composed of 3 Layers and 147 Au Atoms)a method

De (meV)

Zeq (Å)

V0 (×106 meV)

γ (Å−1)

C3 (meV Å3)

Z0 (Å)

RMS (meV)

points hcp vdW-OptB86b recompo hcp vdW-OptB86b av hcp vdW-OptB86bb points hcp vdW-DF2 recompo hcp vdW-DF2 av hcp vdW-DF2b recompo hcp dlDF(Au2)+Das recompo hcp dlDF(Au4)+Das

130.455 131.117 132.400 101.318 101.751 103.839 95.070 91.104

3.4987 3.4793 3.5045 3.6321 3.6070 3.5865 3.7278 3.7143

1.2377 1.2211 1.0037 3.6747 4.3032 4.1687 6.9416 3.8100

2.5687 2.5994 2.5082 2.8627 2.9385 2.9349 3.0610 2.9023

2845 2946 3173 1695 1645 1652 1811 1867

1.3460 1.2758 1.2722 1.6369 1.6179 1.6154 1.5359 1.4933

2.488 1.822 1.696 2.626 1.376 1.293 0.606 0.671

a

The pair potentials are built with the parameters collected in Table 3 and the function expressed in eq 6. The experimental cubic cell parameter, a = 4.0786 Å has been used,48 with the exception of the vdW-OptB86b potential for which the calculated value of a has been employed (Table 2). b Evaluated without including damping functions.

approximation is only about 10% for the atop adsorption site and, hence, its impact should not be significant in MD simulations on the Au(111) surface. The trends found for the effective pair potentials can now be compared with those discussed for the case of the isolated Ar− Au diatomic. The results on the Ar−Au diatomic already showed that the long-range vdW-OptB86b and vdW-DF2 dispersion corrections are rather different. As already mentioned, the gradient correction factor within the nonlocal vdW internal functional of the vdW-DF2 approach is larger than that included in the vdW-OptB86b approach and was determined from atomic scaling laws. This difference qualitatively explains why the vdW-DF2 long-range tail decays more rapidly. As a validation system for the external exchange functional, the Ar−Au diatomic is less useful, as expected from the open-shell character of the Au atom. As highlighted by Halonen and collaborators in theoretical studies of metal clusters,59,63,64 diatomic metal fragments are better used as the building blocks to design accurate models of the interaction. Previous studies on the He/TiO2(110) interaction have shown that vdW-corrected DFT-based approaches failing on molecular models of the surface17 can give a reasonable description when its extended nature is accounted for.19 In contrast, a good behavior of the dlDF+Das approach on small cluster models17 has been found to be reflected when using periodic models.19 Besides the transferability properties of the Das parametrization, this can be understood by considering that the dispersionless correlation effects provided by the dlDF approach are short-ranged, and the partial single-determinant (exact) treatment of the exchange contribution. By application of the dlDF+Das approach, the convergence of both (repulsive) dispersionless and (attractive) dispersion-accounting contributions upon increasing the number of Au2 subunits modeling the metal substrate has allowed to ensure the adequacy of the treatment.

calculated compare very well with those obtained from the pair potentials after the recomposition, proving that the determination of pair potentials is consistent. The interaction potentials resulting of the recomposition from the average pair potentials lead to parameters in less good agreement with those determined from the energies calculated on the hcp site, as could be expected. Still, the agreement is good enough to conclude that these averaged potentials could be used in cases where the adsorption site of the rare gas atom is unknown as, for instance, in MD simulations. As can be seen in Table 4, the vdW-DF2 C3 parameters are very similar to the dlDF+Das ones, corresponding to superimposed curves for Z(Ar−Au) > 4 Å. The dlDF+D as recomposed potentials are also almost identical (Figure 7), showing that the small differences in the repulsive dlDF part of the corresponding pair potentials are smoothed out at the recomposition step. Moreover, the C3 parameters extracted from either vdW-DF2 calculations or dlDF+Das pair potentials are in good agreement with previous reported results,11−13 as well as the well depth values. Notice how vdW-OptB86b C3 parameters are modified in going from Table 2 to Table 4, showing that the supercell size must be large enough to provide a well converged potential when this functional is used. The potential corrugation resulting from using the pair potentials is also shown in Figures 3 and 4 and compared with the VASP results. The application of the different pair potentials leads to a very similar potential corrugation on the Au(111) surface, whereas only small differences are observed on the Au(100) surface. Although the main trends of the direct VASP calculations are found, the exact values are not reproduced. In particular, the potential corrugation at the atop adsorption site is less well reproduced. The fact that all pair potentials lead to very similar potential corrugations indicates that the lack of accuracy is related to the use of pair potentials for the recomposition of the global interaction potentials. In particular, many-body interaction terms are certainly missed when the global interaction is recomposed from effective 2-body pair interactions. Previous benchmarking at the CCSD(T) level on the Xe + Mg(0001) system by Voloshina62 highlighted that surface many-body correlation effects play a relevant role for the adsorption of Xe onto the atop site, at a variance with the hollow site. These effects are associated with attractive dispersionless correlation contributions which cannot be captured through dlDF calculations using small clusters and a periodic extension, such as that implemented by de Lara-Castells et al.,18 is probably necessary. In any case, the error introduced with the pairwise

5. CONCLUSIONS This work has been mainly aimed at determining the pair potentials which are necessary in molecular dynamics simulations of the Ar atom motion on gold surfaces. Two extremely different approaches have been tested and compared. The first one consisted of the decomposition of the global interaction potential calculated by means of plane-wave-based periodic calculations with state-of-the-art vdW-corrected approaches and the VASP code. Two treatments involving nonlocal vdW functionals have been examined: the vdW-DF2 approach of Lee et al.41 and the vdW-OptB86b approach of I

DOI: 10.1021/acs.jpca.5b03769 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A Klimes̃ et al.40 The second strategy has been successfully applied to the calculation of global interaction potential for van der Waals dominated adsorbate−surface interactions,17−20 having been adapted in this work to direct determination of the pair potentials. In this way, the dispersionless DFT developed by Pernal et al.53 has been first applied for the determination of the repulsive pair potential region, and then coupled with ab initio calculations for the attractive dispersion contribution. Specifically, the attractive part was calculated using the D as functional proposed by Szalewicz and collaborators,53 with the dispersion coefficients and damping factors extracted through SAPT-PBE0 calculations. The pair potentials parameters were obtained from the interaction of an Ar atom with small gold clusters, Au2 and Au4. For benchmarking purposes, CCSD(T) calculations were also accomplished on Ar + Au and Ar + Au2 molecular systems. In concluding, our analysis indicates that pair potentials, to be used in further MD simulations of the Ar atom on gold surfaces, are better determined from the application of the dlDF +Das and vdW-DF2 approaches. Although the combined DFT/ post-HF dlDF+Das approach has been applied using a very small cluster model of the extended system, the vdW-corrected (vdW-DF2) DFT-based computations were realized using a periodic surface representation. Despite being very different approaches, the two sets of pair potentials produced extremely similar recomposed (global) interaction potentials for the Ar atom approaching the gold surface. Aimed to determine scattering data and accommodation coefficients through MD simulations of the Ar atom collision dynamics on the target surface, work is in progress to use the calculated pair potentials. Further studies considering different metal surfaces are necessary to analyze how widely applicable and general the strategy applied in this work is for MD pair potential determinations. In fact, the emergence of a rapid growing number of vdW-corrected DFT-based treatments has made convenient the parallel development and application of costefficient treatments mainly grounded on post-HF ab initio theory for cross-validation purposes.

■

FIS2011-2861-C02-01 from DGI, Spain (FEDER), and the Supercomputer Center Cesga. The aid of the COST Action CM1405 “Molecules in Motion” (MOLIM) is also much appreciated.

■

(1) Léonard, C.; Brites, V.; Tung Pham, T.; To, Q.-D.; Lauriat, G. Influence of the Pairwise Potential on the Tangential Momentum Accommodation Coefficient: a Multi-Scale Study Applied to the Argon on Pt(111) System. Eur. Phys. J. B 2013, 86, 164. (2) Tung Pham, T.; To, Q.-D.; Lauriat, G.; Léonard, C.; Vo, V. H. Effects of Surface Morphology and Anisotropy on the Tangential Momentum Accommodation Coefficient Between Pt(100) and Ar. Phys. Rev. E 2012, 86, 051201. (3) To, Q.-D.; Tung Pham, T.; Brites, V.; Léonard, C.; Lauriat, G. Multi-Scale Study of Gas Slip Flows in Nano-Channels. J. Heat Transfer 2015, 137, 091002. (4) To, Q.-D.; Léonard, C.; Lauriat, G. Free-Path Distribution and Knudsen-Layer Modeling for Gaseous Flows in the Transition Regime. Phys. Rev. E 2015, 91, 023015. (5) Patkowski, K. On the Accuracy of Explicitly Correlated CoupledCluster Interaction Energies - Have Orbital Results Been Beaten Yet? J. Chem. Phys. 2012, 137, 034103. (6) Heinz, H.; Vaia, R. A.; Farmer, B. L.; Naik, R. R. Accurate Simulation of Surfaces and Interfaces of Face-Centered Cubic Metals Using 12−6 and 9−6 Lennard-Jones Potentials. J. Phys. Chem. C 2008, 112, 17281−17290. (7) Kimura, Y.; Qi, Y.; Cagin, T.; Goddard, W. A., III. The Quantum Sutton-Chen Many-Body Potential for Properties of FCC Metals. MRS Symp. Ser. 1999, 554, 43−48. (8) Daw, M. S.; Foiles, S. M.; Baskes, M. I. The Embedded-Atom Method: a Review of Theory and Applications. Mater. Sci. Rep. 1993, 9, 251−310. (9) Grochola, G.; Russo, S. P.; Snook, I. K. On Fitting a Gold Embedded Atom Method Potential Using the Force Matching Method. J. Chem. Phys. 2005, 123, 204719. (10) Lee, B.-J.; Shim, J.-H.; Baskes, M. I. Semiempirical Atomic Potentials for the FCC Metals Cu, Ag, Au, Ni, Pd, Pt, Al, and Pb Based on First and Second Nearest-Neighbor Modified Embedded Atom Method. Phys. Rev. B: Condens. Matter Mater. Phys. 2003, 68, 144112. (11) Chizmeshya, A.; Zaremba, E. The Interaction of Rare Gas Atoms with Metal Surfaces: a Scattering Theory Approach. Surf. Sci. 1992, 268, 432−456. (12) Chizmeshya, A.; Zaremba, E. Interaction of Rare Gas Atoms with Metal Surface: A Pseudo-Potential Approach. Surf. Sci. 1989, 220, 443−470. (13) Tang, K. T.; Toennies, J. P. Recalculation of Physisorption Potentials of Rare Gases on Noble Metals. Surf. Sci. Lett. 1992, 279, L203−L206. (14) Zaremba, E.; Kohn, W. Van der Waals Interaction Between an Atom and a Solid Surface. Phys. Rev. B 1976, 13, 2270−2285. (15) Vidali, G.; Ihm, G.; Kim, H. Y.; Cole, M. W. Potentials of Physical Adsorption. Surf. Sci. Rep. 1991, 12, 133−181. (16) Liebrecht, M.; Cole, M. W. Many-Body Effects in Physical Adsorption. J. Low Temp. Phys. 2012, 169, 316−323. (17) de Lara-Castells, M. P.; Stoll, H.; Mitrushchenkov, A. O. Assessing the Performance of Dispersionless and DispersionAccounting Methods: Helium Interaction with Cluster Models of the TiO2(110) Surface. J. Phys. Chem. A 2014, 118, 6367−6384. (18) de Lara-Castells, M. P.; Stoll, H.; Civalleri, B.; Causá, M.; Voloshina, E.; Mitrsuchenkov, A. O.; Pí, M. A Combined Periodic Density Functional and Incremental Wave-Function-based Approach for the Time-Resolved Dynamics of 4He Nanodroplets on Surfaces: 4 He/Graphene. J. Chem. Phys. (Communication) 2014, 141, 151102. (19) de Lara-Castells, M. P.; Aguirre, N. F.; Stoll, H.; Mitrsuchenkov, A. O.; Mateo, D.; Pí, M. Unraveling the 4He Droplet-Mediated SoftLanding from Ab Initio-Assisted and Time Resolved Density

ASSOCIATED CONTENT

S Supporting Information *

Figure S1: Adsorption sites on the Au(111) and Au(100) surfaces. Figure S2: Structural model considered for the Ar atom interaction with the Au2 cluster. Figure S3: Structural model considered for the Ar atom interaction with the Au4 cluster. The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/ acs.jpca.5b03769.

■

REFERENCES

AUTHOR INFORMATION

Corresponding Author

*C. Léonard. Phone: +33 (0)1 60 95 73 18. Fax: +33 (0)1 60 95 73 20. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors thank Doctor Alexander O. Mitrushchenkov, Université Paris-Est Marne-la-Vallée, for technical assistance. M. P. de. L. C acknowledges the Laboratoire Modélisation et Simulation Multi Echelle (MSME, UPEMLV) and CNRS for an invited professor fellowship as well as the Grant No. J

DOI: 10.1021/acs.jpca.5b03769 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A Functional Simulations: Au@4He300/TiO2(110). J. Chem. Phys. (Communication) 2015, 142, 131101. (20) de Lara-Castells, M. P.; Mitrushchenkov, A. O.; Stoll, H. Combining Density Functional and Incremental Post-Hartree-Fock Approaches for van-der-Waals Dominated Adsorbate-Surface Interactions: Ag2/Graphene. J. Chem. Phys. 2015, 143, 102804. (21) Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M. Molpro: a General-Purpose Quantum Chemistry Program Package. WIREs Comput. Mol. Sci. 2012, 2, 242−253. (22) Version 2012.1, MOLPRO is a package of ab initio programs written by Werner, H.-J., Knowles, P. J., Knizia, G., Manby, F. R., Schütz, M., Celani, P., Korona, T., Lindh, R., Mitrushenkov, A., Rauhut, G. et al., 2012 See http://www.molpro.net for more details. (23) Knowles, P. J.; Hampel, C.; Werner, H.-J. Coupled Cluster Theory for High Spin, Open Shell Reference Wave Functions. J. Chem. Phys. 1993, 99, 5219−5227; J. Chem. Phys. 2000, 112, 3106−3107. (24) Watts, J. D.; Gauss, J.; Bartlett, R. J. Coupled-Cluster Methods with Noniterative Triple Excitations for Restricted Open-Shell Hartree-Fock and Other General Single Determinant Reference Functions. Energies and Analytical Gradients. J. Chem. Phys. 1993, 98, 8718−8733. (25) Peterson, K. A.; Dunning, T. H., Jr. Accurate Correlation Consistent Basis Sets for Molecular Core-Valence Correlation Effects. The Second Row Atoms Al - Ar, and the First Row Atoms B - Ne Revisited. J. Chem. Phys. 2002, 117, 10548−10560. (26) Peterson, K. A.; C. Puzzarini, C. Systematically Convergent Basis Sets for Transition Metals. II. Pseudopotential-Based Correlation Consistent Basis Sets for the Group 11 (Cu, Ag, Au) and 12 (Zn, Cd, Hg) Elements. Theor. Chem. Acc. 2005, 114, 283−296. (27) Figgen, D.; Rauhut, G. M.; Dolg, G.; Stoll, H. Energy-Consistent Pseudopotentials for Group 11 and 12 Atoms: Adjustment to MultiConfiguration Dirac-Hartree-Fock Data. Chem. Phys. 2005, 311, 227− 244. (28) Boys, S. F.; Bernardi, F. The Calculation of Small Molecular Interactions by the Differences of Separate Total Energies. Some Procedures With Reduced Errors. Mol. Phys. 1970, 19, 553−566. (29) Halkier, A.; Helgaker, T.; Jørgensen, P.; Klopper, W.; Koch, H.; Olsen, J.; Wilson, A. K. Basis-Set Convergence in Correlated Calculations on Ne, N2, and H2O. Chem. Phys. Lett. 1998, 286, 243−252. (30) Helgaker, T.; Klopper, W.; Koch, H.; Noga, J. Basis-Set Convergence of Correlated Calculations on Water. J. Chem. Phys. 1997, 106, 9639−9646. (31) Kresse, G.; Furthmüller, J. Efficiency of Ab Initio Total Energy Calculations for Metals and Semiconductors using a Plane-Wave Basis Set. Comput. Mater. Sci. 1996, 6, 15−50. (32) Kresse, G.; Furthmüller, J. Efficient Iterative Schemes for Ab Initio Total-Energy Calculations Using a Plane-Wave Basis Set. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 11169−11186. (33) Blöch, P. E. Projector Augmented-Wave Method. Phys. Rev. B: Condens. Matter Mater. Phys. 1994, 50, 17953. (34) Kresse, G.; Joubert, D. From Ultrasoft Pseudopotentials to the Projector Augmented-Wave Method. Phys. Rev. B: Condens. Matter Mater. Phys. 1999, 59, 1758. (35) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865; Phys. Rev. Lett. 1997, 78, 1396. (36) Zhang, Y.; Yang, W. Comment on -Generalized Gradient Approximation Made Simple-. Phys. Rev. Lett. 1998, 80, 890. (37) Perdew, J. P.; Ruzsinszky, A.; Csonka, G. I.; Vydrov, O. A.; Scuseria, G. E.; Constantin, L. A.; Zhou, X.; Burke, K. Restoring the Density-Gradient Expansion for Exchange in Solids and Surfaces. Phys. Rev. Lett. 2008, 100, 136406; Phys. Rev. Lett. 2009, 102, 039902. (38) Grimme, S. Semi-Empirical GGA-type Density Functional Constructed With a Long-Range Dispersion Correction. J. Comput. Chem. 2006, 27, 1787−1799. (39) Grimme, S.; Ehrlich, S.; Goerigki, L. Effect of the Damping Function in Dispersion Corrected Density Functional Theory. J. Comput. Chem. 2011, 32, 1456−1465.

(40) Klimeš, J.; Bowler, R.; Michaelides, A. Van der Waals Density Functionals Applied to olids. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 83, 195131. (41) Lee, K.; Murray, E. D.; Kong, L.; Lundqvist, B. I.; Langreth, D. C. Higher-Accuracy van der Waals Density Functional. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 82, 081101(R). (42) Buckingham, R. A. The Classical Equation of State of Gaseous Helium, Neon and Argon. Proc. R. Soc. London, Ser. A 1938, 168, 264− 283. (43) Tang, K. T.; Toennies, J. P. An Improved Simple Model for the van der Waals Potential Based on Universal Damping Functions for the Dispersion Coefficients. J. Chem. Phys. 1984, 80, 3726−3741. (44) Hopkins, W. S.; Woodham, A. P.; Plowright, R. J.; Wright, T. G.; Mackenzie, S. R. A Velocity Map Imaging Study of Gold-Rare Gas Complexes: Au-Ar, Au-Kr, and Au-Xe. J. Chem. Phys. 2010, 132, 214303. (45) Knight, A. M.; Stangassinger, A.; Duncan, M. A. Photoionization Spectroscopy of Au-Ar. Chem. Phys. Lett. 1997, 273, 265−271. (46) Gardner, A. M.; Plowright, R. J.; Watkins, M. J.; Wright, T. G.; Breckenridge, W. H. Theoretical Study of the X2Σ+ States of the Neutral CM-RG Complexes (CM = coinage metal, Cu, Ag, and Au and RG = rare gas, He-Rn). J. Chem. Phys. 2010, 132, 184301. (47) Bruch, L. W. Theory of Physisorption Interactions. Surf. Sci. 1983, 125, 194−217. (48) JCPDS-International Centre for Diffraction Data. (49) Chen, D.-L.; Al-Saidi, W. A.; Johnson, J. K. The Role of van der Waals Interactions in the Adsorption of Noble Gases on Metal Surfaces. J. Phys.: Condens. Matter 2012, 24, 424211. (50) Sorescu, D. C.; Byrd, E. F. C.; Rice, B. M.; Jordan, K. D. Assessing the Performances of Dispersion-Corrected Density Functional Methods for Predicting the Crystallographic Properties of High Nitrogen Energetic Salts. J. Chem. Theory Comput. 2014, 10, 4982− 4994. (51) Diehl, R. D.; Seyller, Th.; Caragiu, M.; Leatherman, G. S.; Ferralis, N.; Pussi, K.; Kaukasoina, P.; Lindroos, M. The Adsorption Sites of Rare Gases on Metallix Surfaces: a Review. J. Phys.: Condens. Matter 2004, 16, S2839−S2862. (52) Stoll, H. On the Correlation-Energy of Graphite. J. Chem. Phys. 1992, 97, 8448−8454. (53) Pernal, K.; Podeszwa, R.; Patkowski, K.; Szalewicz, K. Dispersionless Density Functional Theory. Phys. Rev. Lett. 2009, 103, 263201. (54) Zhao, Y.; Schultz, N. E.; Truhlar, D. G. Design of Density Functionals by Combining the Method of Constraint Satisfaction with Parametrization for Thermochemistry, Thermochemical Kinetics, and Noncovalent Interactions. J. Chem. Theory Comput. 2006, 2, 364−382. (55) Jeziorski, B.; Moszynski, R.; Szalewicz, K. Perturbation Theory Approach to Intermolecular Potential Energy Surfaces of van der Waals Complexes. Chem. Rev. 1994, 94, 1887−1930. (56) Misquitta, A. J.; Podeszwa, R.; Jeziorski, B.; Szalewicz, K. Intermolecular Potentials Based on Symmetry-Adapted Perturbation Theory with Dispersion Energies from Time-Dependent DensityFunctional Calculations. J. Chem. Phys. 2005, 123, 214103−214114. (57) Woon, D. E.; Dunning, T. H., Jr. Gaussian Basis Sets for Use in Correlated Molecular Calculations. III. The Atoms Aluminum Through Argon. J. Chem. Phys. 1993, 98, 1358−1371. (58) Adamo, C.; Barone, V. Toward Reliable Density Functional Methods Without Adjustable Parameters: The PBE0 Model. J. Chem. Phys. 1999, 110, 6158−6170. (59) Hatz, R.; Korpinen, M.; Hänninen, V.; Halonen, L. Characterization of the Dispersion Interactions and an Ab Initio Study of van der Waals Potential Energy Parameters for Coinage Metal Clusters. J. Phys. Chem. A 2012, 116, 11685−11693. (60) Tonigold, K.; Groß, A. Adsorption of Small Aromatic Molecules on the (111) Surfaces of Noble Metals: A Density Functional Theory Study with Semiempirical Corrections for Dispersion Effects. J. Chem. Phys. 2010, 132, 224701. K

DOI: 10.1021/acs.jpca.5b03769 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A (61) The NUMEROV method, written by Senekowitsch, J., Rosmus, P., Johan Wolfgang Goethe Universit-t, Frankfurt-am-Main, Germany, 1985. (62) Voloshina, E. Local Correlation Method for Metals: Benchmarks for Surface and Adsorption Energies. Phys. Rev. B: Condens. Matter Mater. Phys. 2012, 85, 045444. (63) Hatz, R.; Hänninen, V.; Halonen, L. Dispersion Interactions in Small Zinc, Cadmium, and Mercury Clusters. J. Phys. Chem. A 2014, 118, 5734−5740. (64) Hatz, R.; Hänninen, V.; Halonen, L. Pair-Potential Approach to Accurate Dispersion Energies between Group 12 (Zn, Cd, Hg) Clusters. J. Phys. Chem. A 2014, 118, 12274−12279.

L

DOI: 10.1021/acs.jpca.5b03769 J. Phys. Chem. A XXXX, XXX, XXX−XXX