Global Potentials for the Interaction between Rare Gases and

Apr 30, 2013 - Global potentials for the physisorption of rare-gas atoms on graphene and graphite, amenable for a variety of dynamics simulations, are...
1 downloads 15 Views 2MB Size
Article pubs.acs.org/JPCC

Global Potentials for the Interaction between Rare Gases and Graphene-Based Surfaces: An Atom−Bond Pairwise Additive Representation Massimiliano Bartolomei,*,† Estela Carmona-Novillo,† Marta I. Hernández,† José Campos-Martínez,† and Fernando Pirani‡ †

Instituto de Física Fundamental, Consejo Superior de Investigaciones Científicas (IFF-CSIC), Serrano 123, 28006 Madrid, Spain Dipartimento di Chimica, Università di Perugia, Perugia, Italia



ABSTRACT: Global potentials for the physisorption of raregas atoms on graphene and graphite, amenable for a variety of dynamics simulations, are reported. An atom−bond pairwise additive form of the potential is used, where the interaction pairs, represented by proper analytical functions, are constituted by the Rg atom (Rg = He, Ne, Ar, Kr) and the C−C bonds of the graphene sheet(s). The parameters of the atom−bond pair potential, derived from the polarizability of the interacting partners, are fine-tuned, exploiting calculations of the prototypical Rg−coronene system using high-level electronic structure methods and large basis sets. The atom−graphene/graphite potential is further expanded in a Fourier series, and it is found that for an accurate representation of the interaction only a small number of corrugation terms need to be added to the laterally averaged potential. Furthermore, this corrugation part of the potential is both identical for Rg−graphene and Rg− graphite; in other words, inner layers of graphite only play a role in the laterally averaged Rg−graphite potential. For all systems, the hollow at the center of the carbon ring is the preferred adsorption site, although diffusion barriers are low. The present results compare well with previous data regarding well depths and equilibrium distances at different adsorption sites and, for graphite, the long-range dispersion coefficient C3. In addition, binding energies (eigenvalues of the laterally averaged potentials) are in a fairly good agreement with experimental determinations, providing further support for the reliability of the potentials. SAPT).7 Along these lines, another example is that of ref 8, where CCSD(T) interaction energies of argon with benzene and coronene are used to correct DFT calculations of the raregas physisorption on graphene and graphite. In addition, empirical potentials, besides being backed by experiments, usually have the advantage of being represented by global analytical functions. For rare gas−graphite, they date back to the 1980s9−14 and are built from experiments such as atom-surface scattering, electron diffraction, or isosteric heat adsorption measurements. This body of work was reviewed by Vidali et al.,15 where data relevant for the assessment of theoretical calculations are provided. Finally, we must stand out the pairwise additive approach, where the total potential is obtained as a sum of two-body interactions between the adatom and the constituent atoms of the substrate,10,16 recently employed, for instance, in the description of the interaction between water molecules and graphite.17 Our general goal is to obtain potentials to represent the interactions between molecules and graphene-based surfaces, with the requirement of reliability as well as suitability for the

1. INTRODUCTION The interaction of atoms and molecules with graphene-based surfaces has attracted considerable attention because of several potential applications in nanotechnology and material sciences. For example, the unique physical properties of graphene allow the detection of a single atom (or molecule) adsorbed on its surface, providing a resolution that is unmatched by any other detection technique.1 A reliable description of noncovalent interactions with graphitic surfaces is by no means a simple task for contemporary computational chemistry. The use of accurate wave function-based methods, such as coupled-cluster with single, double, and perturbative triple excitations (CCSD(T)), is computationally prohibitive for such extended systems. The performance of density functional theory (DFT) methods depends to a large extent on the functional employed.2,3 Nowadays, the best choice for both the density functional and a posteriori dispersion correction scheme seems to be not univocal, although much progress has been recently done.4,5 An interesting alternative consists of the use of highly accurate methods on smaller prototypical systems, whose results are then extended or extrapolated to the complete graphene-like planes. That is the case of ref 6, where the interaction between water and a series of acenes of increasing size was determined using DFT-symmetry adapted perturbation theory (DFT© XXXX American Chemical Society

Received: February 15, 2013 Revised: April 24, 2013

A

dx.doi.org/10.1021/jp401635t | J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C

Article

use in a variety of dynamics simulations such as scattering,18−20 diffusion, 21 transmission of molecules through porous graphene,22 and so on. To this end, the potential should be global (covering the full range of positions of the atom over the surface) and it should be given by a simple and compact functional form with a limited number of parameters having a physical meaning. The atom−bond pairwise additive representation, proposed by one of the authors and collaborators,23,24 meets these requirements. In this approach, the van der Waals interaction energy between an atom and a molecule (or between an atom and a solid substrate) is obtained as a sum of atom−molecular bond pairwise contributions. It is based on the concept of bond polarizability and its additivity property25 and makes explicit use of the bond polarizability tensor components26 for the evaluation of the atom−bond contributions. The approach is conceptually similar to the well-known atom−atom pairwise schemes, with the extra benefit of naturally including high-order (three-body) effects. A more realistic description of the interaction anisotropy is then achieved because, with respect to atom−atom schemes, the present formulation accounts better for the electronic charge distribution along the molecular frame. So far the atom−bond approach has been successfully applied to the study of the interaction between a number of closed shell atoms/ions and hydrocarbons of different nature.24,27−37 In this work, we apply the atom−bond approach to study the physical interaction between a rare gas (Rg = He, Ne, Ar, Kr) atom and graphene/graphite. The strategy adopted shares features from some of the approaches outlined above and basically involves three steps. First, electronic structure calculations are performed for Rg−coronene, a smaller but affordable system, to obtain accurate interaction energies.6,8,38 In this regard, we have used the DFT-SAPT method as well as “coupled” supermolecular second-order Møller−Plesset perturbation theory (MP2C),39 with extended basis sets involved. In a second step, these results are used to test and fine-tune the parameters of the Rg−(C−C) pair potential, the basic ingredient of the atom−bond model. Finally, using these optimized parameters, Rg−graphene/graphite interactions are computed within the atom−bond approach by summing over all atom−bond pairs in the extended system until convergence. A global interaction potential is obtained in a fairly compact form as an expansion in a Fourier series involving a limited number of terms. The paper is organized as follows. In Section 2, quantumchemistry calculations for Rg−coronene are described and details of the atom−bond model and its application to Rg− graphene/graphite systems are provided. In Section 3, results are presented and compared with previous calculations as well as with measured binding energies. Finally, conclusions are given in Section 4.

Figure 1. (a) Investigated rare-gas−coronene complex. The hollow (H), top (T), and bridge (B) geometries correspond to three different perpendicular approaches of the rare-gas atom on the coronene plane. (b) Unit cell of graphene/graphite. Carbon atoms and C−C bonds of the outermost layer (graphene) are shown in dark green, while those of the next layer (α-graphite) are given in light green. Sites H, T, and B are also indicated.

approximation of a graphene plane,6,8,38 with the advantage that it can be treated with high-level theories. Calculations were carried out by using the DFT-SAPT7 and MP2C39 methods as implemented in the Molpro2009.1 package.40 The choice of these methodologies comes from an optimal compromise between accuracy and computational cost. In fact, for a proper estimation of these weak interactions, the use of large basis sets (including diffuse functions) is crucial to fully recover the elusive long-range dispersion contribution. This stringent requirement prevents the practical use of the paradigmatic CCSD(T) method for this system due to its higher computational demand. We therefore rely, on the one hand, on DFT-SAPT, a well-established method which combines a DFT description of the intramonomer electron correlation with a perturbative treatment of the intermolecular non covalent contribution, and whose accuracy is in general comparable with that of CCSD(T).7,41 On the other hand, we employ the more recent MP2C method, where the deficiency of MP2 in describing the uncoupled Hartree−Fock (HF) dispersion contribution is corrected with a term determined by using time-dependent DFT response theory.39 The accuracy of MP2C has been recently assessed by extensive calculations of Rg−fullerene interaction energies,42 and it was concluded that the MP2C results are in good agreement with DFT-SAPT calculations, with the advantage of a lower computational cost. An important issue has been the choice of the basis set needed for obtaining reliable interaction energies within a practical calculation time. After several tests on the lightest He−coronene system, we have concluded that at least the augcc-pVTZ43 basis set is required. Inclusion of diffuse functions is fundamental, as calculations with smaller (cc-pVTZ43) basis sets underestimate interaction energies by ∼50% (near the equilibrium position). However, problems of linear dependences appear when employing the aug-cc-pVTZ basis set for the

2. THEORETICAL METHODOLOGY 2.1. Rg−Coronene Calculations. The Rg−coronene system, depicted in Figure 1 (top panel), is chosen as the prototype for carrying out accurate quantum-chemical calculations. Coronene (C24H12) is a planar polycyclic aromatic hydrocarbon (PAH) that shares with graphene a repeating unit consisting of a benzene-type ring. Moreover, it is the smallest PAH for which all carbon atoms of the central ring are chemically bound to other carbon atoms (not hydrogen). Therefore, coronene can be considered as the smallest B

dx.doi.org/10.1021/jp401635t | J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C

Article

coronene molecule. To minimize these issues we have used a modified version of the aug-cc-pVTZ basis set, with the exponents of the most diffuse functions of each angular momentum type multiplied by 2.3, as suggested in ref 6. Test calculations on helium−anthracene using the standard aug-ccpVTZ and aug-cc-pVQZ basis sets show that this scaling of the exponents has a negligible effect on the net interaction energy and that the employed modified set is large enough for the considered hydrocarbon. Tests enlarging the basis set for the rare gas atom were also made. It was found, for instance, that interaction energies near the equilibrium position become deeper by ∼3% with the enlarged aug-cc-pV5Z(Rg)/aug-ccpVTZ(coronene) basis set (if compared with the aug-ccpVTZ(Rg)/aug-cc-pVTZ(coronene) results). All reported calculations have been carried out with the larger basis set, and MP2C interaction energies have been corrected for the basis set superposition error by the counterpoise method of Boys and Bernardi.44 To make both DFT-SAPT and MP2C calculations tractable, the density-fitting method45 has been applied to approximate the two-electron repulsion integrals. In DFT-SAPT, the exchange-correlation functional has been approximated by the Perdew−Burke−Ernzerhof (PBE) generalized gradient functional,46 whose asymptotic behavior has been corrected as described by Grüning et al.47 In this approach, the difference between the ionization and HOMO energies of each monomer is needed; for all monomers, we took the experimental ionization potentials of ref 48, while the HOMO energies have been calculated at the DFT level of theory with the corresponding basis set as used for DFT-SAPT computations. Finally, the δEHF contribution49 was added to the DFT-SAPT energies; this is a term obtained at the HF level of theory that accounts for induction and exchange-induction effects higher than second-order and which is generally small for neutral or highly symmetric systems even if not negligible. For coronene, the experimental C−C bond lengths and C− C−C angles appropriate for graphite were employed (1.420 Å and 120°, respectively50), whereas C−H bond lengths and C− C−H angles were chosen to be 1.09 Å and 120°, respectively. Interaction energies have been computed for three different geometries named hollow (H), top (T), and bridge (B), identifying different perpendicular approaches of the atom onto the coronene plane, as shown in Figure 1. A large grid of intermolecular distances (from 2.5 to 9 Å) has been used for the MP2C calculations, while a shorter one (ranging from 2.5 to 4.5 Å) has been considered for the DFT-SAPT calculations due to its high computational cost (especially for Ar and Kr). Results are reported in Figures 2−5 together with those obtained with the atom−bond approach, described in the next section. 2.2. Atom−Bond Pairwise Additive Approach. In the atom−bond approach,23,24 the Rg−coronene interaction is represented by considering the chemical bonds of the molecule as the interaction centers on the molecular frame. The total intermolecular energy Vtot is given as a sum of atom−bond pair contributions (i), each one represented by a potential term Vi: Vtot =

Figure 2. He−coronene intermolecular interaction for the three investigated geometries. The solid lines correspond to the atom−bond results for the optimized parameters of Table 1, while dashed lines refer to those obtained with the initial parameters (see text); open circles and full squares correspond to the MP2C and DFT-SAPT calculations, respectively.

bond is here assumed to be as an independent diatomic subunit having electronic charge distribution of cylindrical symmetry. Moreover, the reference point for each bond is set to coincide in the present application with the geometric bond center because the dispersion center and the bond center coincide for C−C and are almost coincident for C−H.23,24 The parametrization adopted for the pair potential Vi is of the Improved Lennard-Jones (ILJ) type:51 V i(R i , θi) = f (xi) εi(θi) 6⎤ ⎡ ⎛ 1 ⎞n(xi) n(xi) ⎛ 1 ⎞ ⎥ 6 =⎢ − ⎜ ⎟ ⎜ ⎟ ⎢ n(xi) − 6 ⎝ xi ⎠ n(xi) − 6 ⎝ xi ⎠ ⎥⎦ ⎣

where xi is a reduced distance Ri xi = R mi(θi)

(3)

and where εi(θi) and Rmi(θi) are the well depth and the location of the atom−bond interaction, respectively. They vary with the Jacobian angle as24

∑ Vi (1)

i

(2)

i

The atom−bond pair potential, V , depends on both the distance Ri between the atom and the bond reference point as well as on θi, the Jacobi angle describing the orientation of the atom relative to the bond axis. It is important to note that each C

εi(θi) = εi⊥sin 2(θi) + εi cos2(θi)

(4)

⊥ R mi(θi) = R mi sin 2(θi) + R mi cos2(θi)

(5)

dx.doi.org/10.1021/jp401635t | J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C

Article

Figure 5. As in Figure 2, for the Kr−coronene interaction.

Figure 3. As in Figure 2, for Ne−coronene.

so that ε⊥i , ε∥i , R⊥mi, and R∥mi are, respectively, well depths and equilibrium distances for the perpendicular and parallel approaches of the atom to the bond. In this way, f(xi), the reduced form of the atom−bond potentials,24 is taken to be the same for all relative orientations, as discussed in refs 52−55. Finally, in eq 2, the n exponent is given by the following expression:24

n(xi) = β + 4.0xi2

(6)

where β is a parameter defining the shape of the potential and depends on the nature and hardness of interacting particles. It is worth noting that the ILJ function used here is much more flexible than its classic Lennard-Jones(12,6) precursor and gives a more realistic representation of both the size repulsion and the long-range dispersion attraction, as discussed in detail in ref 51. In particular, at intermediate and large intermolecular distances such model properly accounts for the full dispersion attraction due to dipole−dipole and higher order multipole interaction and for possible damping effects as the distance decreases. This improvement in the potential formulation is crucial to correctly predict the long-range attraction features in complex polyatomic systems, exhibiting many interaction centers distributed on the molecular frame, as the present ones. For the Rg−coronene system, we count 30 Rg−(C−C) and 12 Rg−(C−H) contributions. In the present study, the β parameter has been fixed to 8.5 and 9.0 for the Rg−(C−C) and Rg−(C−H) interactions, respectively, which are typical values for van der Waals interactions in neutral−neutral systems. For the remaining parameters of the Rg−(C−C) pair potential, we have taken as initial values those reported in ref 51, obtained using known empirical correlation formulas23,56 that take into account rare-gas and C−C bond polarizabilities. As for the

Figure 4. As in Figure 2, for Ar−coronene interaction. Additionally, a calculation from ref 8 obtained at the CCSD(T) level of theory is reported in solid diamonds.

D

dx.doi.org/10.1021/jp401635t | J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C

Article

where τ1 and τ2, varying in the interval [−1/2,1/2], give the projection of the atom position onto the unit cell. The interaction potential is expanded in a Fourier series as

parameters of the Rg−(C−H) interaction pairs, we have considered as initial values those of ref 24 and applied them to several rare gas−hydrocarbon systems (see, for instance, refs 24 and 57). These values have been fine-tuned on the basis of a comparison with results of the electronic structure calculations of the previous subsection. The optimized parameters deviate from their initial values not more than a few percent for R⊥m,R∥m and 10% for ε⊥,ε∥, and they are presented in Table 1. The Rg− coronene atom−bond interaction energies are depicted in Figures 2−5, in comparison with the DFT-SAPT and MP2C calculations.

V (z , τ ) =

R⊥m

R∥m

ε⊥

ε∥

β

He−CC Ne−CC Ar−CC Kr−CC He−CH Ne−CH Ar−CH Kr−CH

3.589 3.643 3.851 3.974 3.234 3.316 3.641 3.782

3.936 3.995 4.147 4.249 3.584 3.655 3.851 3.987

0.721 1.297 3.553 4.250 1.364 2.544 4.814 5.667

0.866 1.809 5.207 5.906 0.924 1.782 3.981 4.781

8.5 8.5 8.5 8.5 9.0 9.0 9.0 9.0

(9)

g

where g is the reciprocal lattice vector and the coefficients are given by Wg(z) =

1 aS

∫unit cell exp[−i g·τ]V (z , τ) dτ

V (z , τ ) = V0(z) + V1(z)[cos(2πτ1) + cos(2πτ2) − cos(2π (τ1 + τ2))] + V2(z)[cos(2π (τ1 − τ2)) − cos(2π (2τ1 + τ2)) − cos(2π (τ1 + 2τ2))]

Perpendicular and parallel components of Rm and ε are in angstroms and in millielectronvolts, respectively. β is dimensionless. a

+ V3(z)[cos(4πτ1) + cos(4πτ2) + cos(4π (τ1 + τ2))] + ...

d [3i + 2 d a 2 = [3i − 2

⎡ ℏ2 d 2 ⎤ + V0(z) − ε⎥ψ (z) = 0 ⎢− 2 ⎣ 2m dz ⎦

(12)

in which m is the mass of the incident rare gas atom, V0(z) is the corresponding interaction potential, and ε is the binding energy. The equation has been solved numerically with a Numerov integrator58,59 using typically 10 000 to 20 000 points in the z coordinate depending on the (mass) system.

3. RESULTS AND DISCUSSION 3.1. Rg−Coronene. Results for the Rg−coronene interaction (Rg = He, Ne, Ar, Kr) are reported in Figures 2−5, respectively. In general, it can be noticed that the most attractive interaction corresponds to the hollow configuration, while the top and bridge sites are almost coincident, with the latter being slightly more attractive, as could be anticipated on the basis of steric considerations (see Figure 1). Regarding the electronic structure calculations, it can be seen that both DFTSAPT and MP2C methodologies give consistent interaction energies, the reference DFT-SAPT results being slightly less repulsive in the barrier region. For Ar−coronene, the accurate CCSD(T) calculation of ref 8 at the equilibrium distance is also shown in Figure 4, and it can be seen that present calculations

3 j] 3 j]

(11)

The newly defined coefficients Vi(z) (i = 0−3) are related to Wg(z) as follows: V0 = Re[Wg0], V1 = 2Re[Wg1], V2 = 2Re[Wg2], and V3 = 2Re[Wg3], where the absolute values of g0, g1, g2, and g3 are 0, 4π/3d, 4π/√3d, and 8π/3d, respectively. 2.4. Binding Energies of Rg−Graphene/Graphite. Exploiting the analytical formulation of the interaction, we have computed the binding energies corresponding to the socalled laterally averaged potential, V0, observables that can be compared with available experimental data in some cases or theoretical results in some others. To obtain them, we solve the Schrödinger type equation

2.3. Rg−Graphene/Graphite Interaction: Global Potentials. Rg−graphene interaction potentials have been obtained by generalizing the atom−bond pairwise additive approach, that is, by summing neighboring Rg−(C−C) pair potentials until convergence on the total interaction energy is reached (typically four or five significant digits). For Rg− graphite, the procedure is completed by adding as many graphene sheets as needed to reach convergence with the same criterion. In this work, a model of α-graphite was considered with an interlayer distance of 3.35 Å. It is important to stress that within the present model we are not taking into account possible effects due to the more extended π-electron conjugation in graphene with respect to that in coronene; moreover, in the case of graphite we are considering that each graphene sheet contributes independently to the total interaction and therefore possible screening effects due to the outermost layers are not considered. We expect that the abovementioned effects should give a minor contribution as discussed in next section. A top view of the unit cell of graphene/graphite is shown in Figure 1 (lower panel). The xy plane of the Cartesian coordinate system coincides with on the outermost (graphene) layer. The origin of this system is chosen at the bridge position in the center of the unit cell. The unit lattice vectors are a1 =

(10)

where in the previous equation the integration covers the unit cell and aS is its area. A numerical integration on a grid of ∼400 (τ1, τ2) points was carried out to actually compute such coefficients. It was found that the imaginary part of Wg(z) is negligible and, besides, that Wg(z) ≈ Wg′ when |g| = |g′|. As found here, the expansion of Rg−graphene/graphite involves a very small number of significant terms, and it can be written as

Table 1. Optimized Parameters for the Rg-(C−C) and Rg(C−H) Pair Potentials (see eqs 2−6)a atom−bond pair

∑ Wg(z) exp[ig·τ ]

(7)

where d = 1.42 Å is the C−C distance and i,j are unit vectors parallel to x,y axes. The position of the rare-gas atom is given by r = (z,τ), where z is the distance from the atom to the outermost layer of the substrate and τ = τ1a1 + τ2a 2 (8) E

dx.doi.org/10.1021/jp401635t | J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C

Article

agree quite well with this result. As for the atom−bond interaction model, it can be seen that the curves obtained with the optimized parameters of Table 1 reproduce quite well the DFT-SAPT and MP2C energy profiles. It must be stressed that the optimized atom−bond parameters come from a fine-tuning of the initial values of ref 51 to better reproduce the quantum chemistry calculations for the considered geometries. Nevertheless, it can be noticed that the atom−bond approach also works rather well even in the case of using nonoptimized parameters (see the dashed lines in Figures 2−5). This confirms the versatility of the atom−bond approach, which, starting from parameters originated by meaningful physical properties such as atom and bond polarizabilities, is capable of reproducing the main features of Rg−coronene interaction. 3.2. Rg−Graphene/Graphite Global Potentials. The atom−bond pairwise additive approach has been extrapolated to represent the Rg−graphene/graphite interaction, as described in Section 2.3. The parameters of the atom−bond pair potential are those of Table 1 which optimize the agreement with the DFT-SAPT and MP2C Rg-coronene energies. As an example, we present in Figure 6 a contour

Figure 7. Rare gas−graphene and rare gas−graphite interaction potentials as a function of the distance from the surface. The depicted curves refer to the laterally averaged interaction contribution V0(z), that is, the first term of the expansion given in eq 11.

going from the lighter to the heavier atoms, even if in all cases the increase in De is always ∼10%. Nevertheless, it is interesting to note that the equilibrium distances for graphite almost match the corresponding ones for graphene. Numerical values of the equilibrium features of V0(z) are given in Tables 2 and 3. In Figure 8, we report the coefficients Vi(z) (i = 1,2,3) of eq 11, responsible for the corrugation of the interaction. It has been checked that expansion up to i = 3 is sufficient for an accurate representation of the global potential V(z,τ). This result indicates that the corrugation of these surfaces is small and that dynamics calculations using these potentials should not be too computationally demanding. It can be seen that all anisotropic coefficients have a simple exponential-like behavior, with V1 being the dominant one. Moreover, a general trend for the different rare gases is seen, where the larger coefficients correspond to the heavier atoms. It is important to stress that the coefficients of Figure 8 are identical for both graphene and graphite. This result is a consequence of the short-range character of the corrugation part of the interaction in atom− graphene. Taking into account that within the present model each layer of graphite independently influences the total interaction, all layers, except the outermost one, are too distant from the rare-gas atom for contributing to the corrugation part of the interaction. In other words, the inner layers of graphite only affect the laterally averaged part of the interaction, making it more attractive if compared with graphene. (See Figure 7.) Table 2 collects our results concerning the main features of the Rg−graphite interactions and compares them with results of previous works. Starting with the laterally averaged interaction, present well depth De, equilibrium distance ze, and asymptotic dispersion coefficient C3 (directly provided by

Figure 6. Ar−graphene interaction potential as a function of the surface coordinates (eq 11) for z = 3.4 Å, the equilibrium distance of the laterally averaged interaction. Contour values are in millielectronvolts.

plot of the Ar−graphene potential for an atom-surface distance z = 3.4 Å, which nearly corresponds to the absolute minimum of the interaction. It can be noticed that the present procedure provides a quite smooth potential all over the unit cell. As expected, the minimum of the potential is seen at the hollow position. It is found that corrugation is quite small. Indeed, very low diffusion barriers are noticed with saddle points located at the bridge configuration and with slightly higher values at the top sites. A qualitatively similar behavior is found for the other rare gas atoms studied and in the interactions with graphite. Laterally averaged potentials (V0(z) in eq 11) are shown in Figure 7 as functions of the atom-surface distance, z. For both graphene and graphite, it can be seen that the well depth, De, progressively increases from the lighter to the heavier rare-gas atom. This feature can be interpreted considering the polarizability of the involved rare gas, which goes as 1.4, 2.7, 11.1, and 16.8 au for He, Ne, Ar and Kr, respectively.60 A similar behavior occurs for the equilibrium distance, ze, and this typical trend can be related again to the atomic polarizability which also quantifies the size of the involved rare-gas atom. As a consequence, the van der Waals radius increases when going from He to Kr.48 For a given atom, well depths and long-range attractive tails are larger for graphite than for graphene, as expected, and this difference becomes more important when F

dx.doi.org/10.1021/jp401635t | J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C

Article

Table 2. Well Depth, De (in meV), and Equilibrium Distances ze (in Å) for He, Ne, Ar, and Kr on Different Sites (Hollow, Top and Bridge) of the Graphite Substrate Together with the Values Predicted for the Laterally Averaged Potential V0(z)a De

graphite

ze

C3

Bridge He Ne Ar Kr

17.8 34.5 106.2 133.1

(11.3b)

He Ne Ar Kr

17.7 34.3 105.9 132.8

(11.0b)

He Ne Ar Kr

18.5 35.8 109.0 135.8

(106.4b,100.5d,111.0c,116.0e) (137.9c)

He Ne Ar Kr

17.9 34.8 106.9 133.8

(16.6 ± 0.4f) (32.6f) (96 ± 2f) (125 ± 5f)

(99.8b,106.6c) (133.7c)

(107.0b,106.0c) (133.2c) (13.0b)

3.18 3.22 3.40 3.50 Top 3.18 3.22 3.40 3.51 Hollow 3.13 3.17 3.36 3.48 Average Potential 3.17 3.21 3.39 3.50

(3.66b) (3.70b,3.37c) (3.48c) (3.71b) (3.80b,3.37c) (3.48c) (3.60b) (3.65b,3.32c,3.33e) (3.44c) (2.65f) (2.79f) (3.1 ± 0.1f) (3.2 ± 0.2f)

191.1 370.5 1317.3 1803.3

(180 ± 15f) (346f) (1210f) (1730f)

a

In the last case, the long range C3 dispersion coefficient (in meV·Å3) is also reported. In parentheses previous estimations are given for comparison. DFT/vdW-WF calculations of ref 2. cTheoretical study of ref 61. dDFT/CC calculations of ref 8. eDFT calculations of ref 3. fBest estimate values of ref 15. b

Table 3. Well-Depth Energies De (in meV) at Equilibrium Distances ze (in Å) for He, Ne, Ar, and Kr on Different Sites (Bridge, Top, and Hollow) of the Graphene Substrate Together with the Values Predicted for the Laterally Averaged Potential V0(z)a De

graphene He Ne Ar Kr

16.4 31.7 96.6 120.2

He Ne Ar Kr

16.3 31.6 96.3 119.9

He Ne Ar Kr

17.1 33.0 99.2 122.7

He Ne Ar Kr

16.5 32.0 97.2 120.8

Bridge (10.6b) 3.19 3.23 (78.0b) 3.41 3.52 Top (11.2b) 3.20 3.24 (78.6b) 3.42 3.52 Hollow (11.6b) 3.14 3.18 (78.0b,93.3c) 3.38 3.49 Average Potential 3.18 3.22 3.40 3.51

ze

C4

(3.71b) (3.81b)

(3.75b) (3.87b)

(3.61b) (3.79b,3.29c)

1690.2 3444.7 13053.8 18446.0

a

In the case of V0(z), the long range C4 dispersion coefficient (in meV·Å4) is also reported. Some previous theoretical estimations, in parentheses, are given for comparison. bDFT/vdW-WF calculations of ref 2. cDFT/CC calculations of ref 8.

Figure 8. Corrugation terms Vi(z) of the Fourier series expansion given in eq 11. It has been found that for a specific rare gas the graphene and graphite potential expansions provide identical terms.

the adopted potential model) are compared with the best estimations reported in the work of Vidali et al.15 It can be seen that present values are slightly but systematically larger than

those of ref 15. Note that the largest discrepancy concerns the equilibrium distance for He and Ne. It is worth pointing out that the values selected in ref 15 mostly correspond to the lower limit of the collected values (for instance, ze values for Ne G

dx.doi.org/10.1021/jp401635t | J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C

Article

Table 4. Low-Lying Six Energy Levels εn (in meV) of He, Ne, Ar, and Kr on Graphene and Graphite As Obtained from the Present V0(z) Potentialsa 3

He−graphene 3 He−graphite 4

He−graphene He−graphite

4

Ne−graphene Ne−graphite Ar−graphene Ar−graphite

Kr−graphene Kr−graphite

work

ε0

ε1

ε2

this work this work expb this work this work expc expd this work this work best estimatione this work this work best estimatione expf this work this work

12.08 13.35 11.62 12.63 13.92 12.42 12.27 29.49 32.21 30.1 94.23 103.8 99 ± 4 119 ± 2 118.6 131.4

5.63 6.63 5.38 6.68 7.73 6.59 6.56 24.82 27.40

2.04 2.73 1.78 2.97 3.76 3.05 3.08 20.64 23.08

88.40 97.8

82.82 92.0

114.2 126.9

109.9 122.4

ε3

ε4

ε5

0.50 0.89

0.05 0.20

0.03

1.03 1.56 1.24 1.20 16.93 19.22

0.24 0.52 0.47 0.49 13.68 15.82

0.02 0.13 0.16 0.13 10.87 12.85

77.46 86.4

72.29 81.1

67.33 75.9

105.7 118.1

101.5 113.9

97.48 109.7

a

Experimental determinations are given for comparison. bScattering experiment.9 cScattering experiment.13 dScattering experiment.63 eBest estimate values of ref 15. fIsosteric heat measurement.64

At this point it is worth mentioning the work of Carlos and Cole,10 where a global analytical interaction potential for He− graphite was reported and which has been extensively used to describe such a system so far. We were unable to find detailed values of the interaction features and therefore we have established only a qualitative comparison with Figure 2 of that work,10 where energy profiles of the three typical sites are depicted. The hollow site is the most stable one with a well depth of ∼18 meV, in very good agreement with our estimation, but the equilibrium distance is smaller (∼2.6 Å). Also in accordance with present work, the bridge site is slightly preferred over the top one but the corrugation is somewhat larger: while their difference between hollow and bridge well depths is ∼3 meV, ours is on the order of 1 meV. Table 3 collects present results for the Rg−graphene interaction and compares them with data available from previous calculations. For Ar−graphene at the hollow site, we achieve a good agreement with the DFT/CC approach of ref 8, where accurate CCSD(T) calculations on smaller reference systems are combined with a DFT scheme to provide a reliable description of the interaction with extended graphene-based surfaces. However, this approach has not been used yet for the assessment of the anisotropy or the dependence of the interaction energies on the approaching configuration. The DFT/vdW-WF calculations of ref 2 are more complete but tend to overestimate the values of the equilibrium distances and are not accurate enough to predict reliable values of the lateral variation of the well depths, as already commented for the case of Rg−graphite. Finally, in Table 3, we also report the equilibrium features of the interaction, as well as an estimation of the main dispersion coefficient, C4, of the laterally averaged potential, V0(z). Note that C4 values for graphene, as those of C3 for graphite, are directly provided by the asymptotic behavior of the adopted potential functions. In particular, according to the ILJ potential formulation extended to Rg− monolayer systems, C4 is defined as De·z4e , where De and ze are the equilibrium parameters of the laterally averaged potentials of Table 3. In the case of graphite, C3 is obtained, according to Vidali et al.,11 as De·l3, where De is the well depth of the laterally averaged potentials of Table 2 while l is a distance different

range from 2.78 to 2.99 Å therein). For He, in addition, one should clarify that the ze value reported by Vidali et al. actually corresponds to 2.85 Å in the original reference,62 as noted by Ambrosetti et al.2 Besides this, differences with the reference values are within