Embedded Diatomics-in-Molecules Potential Energy Function for

6842

J. Phys. Chem. B 1998, 102, 6842-6860

Embedded Diatomics-in-Molecules Potential Energy Function for Methyl Radical and Methane on Nickel Surfaces Steven E. Wonchoba and Donald G. Truhlar* Department of Chemistry and Supercomputer Institute, UniVersity of Minnesota, Minneapolis, Minnesota 55455-0431 ReceiVed: March 27, 1998; In Final Form: May 22, 1998

We have developed a potential energy function for the interaction of H, CH3, and CH4 with crystalline Ni surfaces and the dissociation reaction of CH4 on such surfaces. The potential is based on the embedded diatomics-in-molecules (EDIM) formalism, and it involves mixing the semiempirical diatomics-in-molecules (DIM) valence bond method for the covalent part of the system (CH4) with the embedded atom method (EAM) for the metal. A new molecular modeling technique is developed to accomplish this, and it employs a physically motivated mapping of electron density into an effective interatomic distance. The EDIM potential energy surface is calibrated against experimental findings and ab initio electronic structure calculations for CH3 and CH4 adsorbed on Ni(111), and it reproduces the best available binding geometries, energies, and frequencies quite well. Based on the behavior of the methyl torsion mode for CH3 in the 3-fold hollow site, our results indicate that the minimum energy orientation for this site has the H atoms staggered with respect to the three Ni atoms that form the hollow site. This differs from the prediction of previous calculations, and we present a possible explanation for this discrepancy. Finally, we calculated the barrier height for dissociative chemisorption of methane to produce CH3 and H, and we obtained 13.8 kcal/mol, which is within the range of estimates obtained from experiment.

1. Introduction C-H bond activation by catalysts is a topic of great industrial interest and has significant implications for the global economy. Steam reforming, i.e., the reaction of CH4 and water to produce hydrogen and CO, is a key commercial process1-3 for the formation of hydrogen or synthesis gas. Metallic nickel is the catalyst for all commercial steam reforming processes. Furthermore, adsorbed hydrocarbon radicals are likely intermediates in a variety of reactions catalyzed by transition metals.3 As a result, there have been many studies, both experimental and theoretical, of the interaction of hydrocarbons with metal surfaces. Ceyer and co-workers have examined the dissociative chemisorption of methane4-6 and the hydrogenation of ethylene7 on Ni(111) and have proposed mechanisms for both. They have also analyzed CH3 adsorption on Ni(111) and, based on spectroscopic data, were able to determine the specific adsorption site.8 As a first step forward performing atomistic dynamical simulations of the proposed mechanisms for surfacecatalyzed dissociative adsorption and hydrogenation of hydrocarbons, in this paper we develop a potential energy function (PEF) that accurately describes the interaction of these hydrocarbons and their fragments with Ni. Previous theoretical work on this problem has been directed to understanding the geometries and energetics of various adsorption sites, and many calculations9-25 have been carried out for hydrocarbons interacting with Ni surfaces with the goal of attaining such understanding. A long-term goal is to understand the structure and energetics of transition states for the C-H bond activation on metal surfaces; this topic has also been studied, although less widely.16,20,24 Theoretical work on chemisorption of CH3 on metal surfaces has recently been reviewed.25 The goal of the present paper is to create a new PEF for studying the dynamics of methane and its first decomposition

product, methyl radical, on a Ni surface. The PEF presented in this paper incorporates three components. The first of these is an atom-metal potential function called EAM7, and it is used to model the dynamics of an H atom interacting with a metallic surface. EAM7 is a modification of the EAM6 potential function,26 which was presented previously and used for extensive dynamical calculations. The changes made to EAM6 to produce EAM7 are presented in this paper. The second component of the present work is an embedded-diatomics-inmolecules (EDIM) formalism,27 which is used here, with appropriate extensions, to model the dynamics of C-H bonds interacting with a metallic surface. The third component is a gas-phase CH5 potential called J128,29 that was developed to model the dynamics of the gas-phase reaction CH3 + H2 f CH4 + H. The J1 PEF is based on an earlier potential developed by Raff30 for the same reaction. The Raff potential contains several features (functional discontinuities, lack of an out-ofplane bend vibration for methyl) that prohibit its use in reactionpath calculations. The J1 potential features minimal modifications to the Raff function that allow it to be used for dynamics calculations, and in addition, it has a reasonable value for the transition-state barrier height. By using the EDIM formalism to combine features of the EAM7 and J1 potential functions, we are able to develop a PEF that models CH3 and CH4 interacting with the surface of metallic Ni. In this paper, we present the full potential function, and we examine selected adsorption configurations of CH3 and CH4 on the Ni surface. Wherever possible, we compare our calculated geometries, energies, and vibrational frequencies to experiment4,8,31 and previous theoretical results.9-12,17-20,23,24 A distinguishing feature of the PEF presented here is that it includes the full permutational symmetry of the four hydrogen atoms. This is particularly important for trajectory studies (or

S1089-5647(98)01666-6 CCC: $15.00 © 1998 American Chemical Society Published on Web 08/12/1998

Embedded Diatomics-in-Molecules Potential other molecular dynamics studies) of methane-surface collisions since it is not clear at the start of a trajectory which or how many hydrogen atoms will interact most directly with the surface. In addition, geometries in which more than one C-H bond is significantly distorted may be important. The EDIM method27 is based on the embedded atom method32-34 (EAM) combined with semiempirical valence bond theory,35-38 in particular with the method of diatomics in molecules.39-41 The embedded atom method, and hence the current potential function, assumes that the interaction energy of the molecule with the surface is a function of the density created by the Ni surface at the atoms in the adsorbed molecule. The hydrogen-nickel interaction parameters in our methanenickel potential functions were previously optimized so that a single set of parameters yields reasonable results on all three low-index faces: Ni(100), Ni(110), and Ni(111).26 An advantage of a PEF with a parametrization in which the facial characteristics emerge naturally is that it should be suitable for studying the interaction of methyl radical and methane with any facial exposure of Ni (including kinks and steps), but in the present paper we limit our attention to the low-index faces of Ni. In section 2 we discuss the geometries of the various chemisorbed and physisorbed configurations for CH3 and CH4 on Ni. In section 3 we present the new potential function and outline the procedures used for obtaining its parameters. In section 4 we present some calculations with the new potential function and compare them to available experimental measurements and other theoretical calculations, and in section 5 we summarize the conclusions and discuss future work.

J. Phys. Chem. B, Vol. 102, No. 35, 1998 6843

Figure 1. CH3 adsorbed at 3-fold sites on Ni(111). Only primary zone atoms are shown. The CE orientation has the H atoms eclipsed with respect to the three Ni atoms that form the 3-fold hollow site. The C30 orientation has the H atoms rotated 30° from the CE orientation. The C60 orientation has the H atoms staggered with respect to the three Ni atoms. It is rotated 60° from the CE orientation.

Figure 2. CH3 adsorbed at bridge sites on Ni(111). Only primary zone atoms are shown. The B0 orientation is formed by translating the CH3 molecule in the CE orientation to an adjacent bridge site. The B30 orientation is formed by rotating the CH3 molecule in the B0 orientation by 30°.

2. Adsorption Sites In this section, we present an adsorption-configuration nomenclature which is necessary for the discussion of the various configurations of interest in the context of the remainder of the paper. In previous work,26 we made use of an adsorptionsite nomenclature for labeling minimum energy sites and saddle points for a single adatom interacting with a surface. That nomenclature is not sufficient for the present study of the interaction of CHn with a surface because the adsorption of a molecule requires specification of the molecular orientation as well as the site. In general, stationary points for a CH3 molecule adsorbed on a Ni surface have pyramidal molecular geometries with the C atom closest to the surface plane. Each adsorption configuration can be described by the Miller indices of the exposed Ni crystal face on which the CH3 molecule is adsorbed [e.g., (111), (100), or (110)], by the surface site to which the C atom bonds [e.g., atop (when the C atom is directly above a Ni atom), bridge (when the C atom is directly above a bridge formed by two Ni atoms), or center (when the C atom is directly above a hollow center site formed by multiple Ni atoms)], and by the local geometry of the three H atoms with respect to the C atom and the surface (e.g., the equilibrium values of the torsional and bending coordinates). Previous theoretical modeling of CH3 adsorbed on Ni has been addressed to the (111) surface, and therefore we will emphasize this crystal face. Several configurations of CH3 adsorption on Ni(111) are shown in Figures 1-3. For each configuration, the figures show those Ni atoms, called the primary zone, that lie within a sphere (of radius specified below) centered at the C atom. The primary zone atoms are allowed to move in geometry optimizations, frequency calculations, and dynamics calculations. The radius of the sphere differs for each

Figure 3. CH3 adsorbed at atop sites on Ni(111). Only primary zone atoms are shown. The A0 orientation is formed by translating the CH3 molecule in the CE orientation to an adjacent atop site. The A30 orientation is formed by rotating the CH3 molecule in the A0 orientation by 30°.

system, but it is chosen to enclose the smallest possible number of atoms in the primary zone to result in converged adsorption geometries, energies, and vibrational frequencies. For each case, the primary zone is embedded into a surrounding lattice of Ni atoms, called the secondary zone, to preserve the bulk structure of the lattice. The metal atom-metal atom interaction potential of EAM7, like EAM6, has a finite cutoff distance, and the secondary zone is chosen to be large enough that it includes all lattice atoms within the potential cutoff of every primary zone atom. The total number of included atoms in both zones varies from 670 to 689 for calculations on the (111) face of Ni. The radius of the primary zone is chosen such that increasing it does not change the adsorption geometries, energies, and frequencies by more than 1%. For CH3 in a 3-fold site (see Figure 1), the primary zone radius is 4.1 Å, and it contains nine Ni atoms, consisting of six surface level Ni atoms and three second-level Ni atoms. For CH3 at a bridge site (see Figure 2), the primary zone radius is 4.4 Å, and it again contains nine atoms; but now these consist of eight surface level Ni atoms and one second-level Ni atom. For CH3 at an atop site (see Figure 3), the primary zone radius is 5.0 Å, and it contains 10 Ni atoms: seven surface level Ni atoms and three second-level Ni atoms. Further details of the lattices are given in the Appendix.

6844 J. Phys. Chem. B, Vol. 102, No. 35, 1998

Wonchoba and Truhlar

SCHEME 1

Our analysis of the adsorption energies of CH3 on a Ni(111) surface will require a careful specification of sites. To do this as clearly as possible, we have developed the labeling scheme of Scheme 1. In this grid, boldface capital letters represent surface Ni atoms, lower case letters represent subsurface Ni atoms in the second layer, and asterisks represent Ni atoms in the third layer (or vacancies in the second layer). Three-fold sites are formed by three surface atoms in a triangle, e.g., BGH, HMN, HIN, etc. Therefore, HMN is a 3-fold site directly above a second plane vacancy, and HIN is a 3-fold site directly above a second plane atom. We label the adsorption configurations in Figures 1-3 by the location of the carbon atom with respect to the surface. A hollow site in the center of three adjacent Ni atoms is labeled C, a bridge site between two adjacent Ni atoms is labeled B, and a site atop a single Ni atom is labeled A. The CE configuration has the H atoms eclipsing the three Ni atoms that create the hollow 3-fold site. The C30 configuration has the H atoms rotated 30° from the CE configuration, and the C60 configuration has the H atoms rotated 60° from the CE configuration and oriented toward the bridge sites formed by the three Ni atoms that create the hollow 3-fold site. The CE and C60 configurations are also called the eclipsed and staggered 3-fold configurations, respectively. There are two distinct 3-fold hollow sites. One is immediately above a Ni atom in the second atomic layer, and the other is immediately above a vacancy in the second atomic layer. Atoms and molecules adsorbed at these sites are energetically and structurally similar, so henceforth, a 3-fold site will refer to the site immediately above the second plane vacancy. In cases where CH3 and H are adsorbed at adjacent sites, the CH3 radical is situated above the second-layer vacancy and the H atom is situated above the second layer atom. Bridge sites are at MN in this labeling scheme. In a B0 configuration, the H atoms of CH3 are in the same relative orientation as in a CE configuration. That is, start with CH3 in a CE configuration, and translate it as a whole so that the C atom is above MN. In a B30 configuration, the H atoms are in the same orientation as in a C30 configuration. That is, start with a CE configuration, rotate it 30° clockwise, and translate the whole CH3 so that the C atom is above MN. Atop sites are located at Ni atom H in the grid. For the A0 configuration, the hydrogen atoms are in the same relative orientation with respect to carbon as at the CE configuration. The A30 configuration is created by rotating the H atoms in the A0 configuration by 30° clockwise. Similar adsorption configurations occur on the (100) and (110) crystal faces. Figure 4 shows three important examples on the (100) face, namely, CH3 adsorbed on a 1-fold atop site (Figure 4a), at a 2-fold bridge site (Figure 4b), and at a 4-fold hollow center site (Figure 4e). Figure 5 shows three CH3 adsorption configurations on the (110) face, namely, CH3 adsorbed on a 1-fold atop site (Figure 5a), at a 2-fold bridge site (Figure 5b), and at a pseudo-3-fold site (Figure 5c). The

Figure 4. CH3 adsorbed at 1-fold atop (a), 2-fold bridge (b), and 4-fold center (c) sites on Ni(100). Only primary zone atoms are shown.

Figure 5. CH3 adsorbed at 1-fold atop (a), 2-fold bridge (b), and pseudo-3-fold (c) sites on Ni(110). Only primary zone atoms are shown.

geometries, energetics, and frequencies of these adsorption configurations are presented in section 4. Since CH4 is a closed shell system, it does not chemisorb on the surface, but there is a weakly bound physisorbed configuration. Following Yang and Whitten,20 who examined the dissociation of CH4 on Ni(111), we examine several physisorbed configurations of CH4 on Ni, all with the C atom positioned directly above a Ni surface atom. For all three low-index faces [(111), (100), and (110)], we found the most energetically favorable physisorbed configurations to involve tetrahedral CH4

Embedded Diatomics-in-Molecules Potential


TABLE 1: Parameters for EAM7 H/Ni Interaction Which Differ from EAM6 RZ ) -5.957 Å-5 RH1 ) -4.803 eV Å3γH1 RH2 ) -74.556 eV Å3γH2 AH1 ) 1.841 × 1011 eV Å15 DH1 ) -1.304 × 1011 eV Å15 Fcut3 ) 0.009 Å-3 a

Parameters for Effective Charge on H Atoma βZ ) 9.804 Å-4 Parameters for Primary Embedding Energy of H Atom βH1 ) 34.46 Å-3 βH2 ) 7.2784 Å-3 BH1 ) 3.276 × 109 eV Å12 EH1 ) 2.376 × 109 eV Å12 ∆3 ) 0.007 Å-3

γZ ) -3.715 Å-3 b

γH1 ) 0.3855 γH2 ) 1.0154 CH1 ) 1.510 × 107 eV Å9 FH1 ) -1.195 × 107 eV Å9

See eq 2. b See eq 4.

with a single H atom facing the surface and three H atoms pointing away. Henceforth, when we refer to a physisorbed configuration of CH4 on a Ni surface, we are referring to this orientation.

3. Potential Energy Function EDIM models have been presented previously for two- and three-body hydrogenic systems interacting with a metal surface.27 For four- and five-body hydrocarbons (methyl radical and methane), the full valence bond equations become very unwieldy,37,38 and an EDIM formalism based on a full valencebond configuration list would likewise become complicated. The alternative approach employed here is to treat each bond in the adsorbed species as an individual two-body system interacting with the surface and sum the four resulting interactions.28-30 With this approach, the total potential energy, Vtot, is a sum of six terms:

Vtot ) Vmet + VDIM(C-H1-Ni) + VDIM(C-H2-Ni) + VDIM(C-H3-Ni) + VDIM(C-H4-Ni) + Vbend (1) The first term in eq 1, Vmet, is the energy of the metal atoms. Terms two through five in eq 1 are the diatomics-in-molecules (DIM) expressions of the four C-H bonds interacting with the Ni surface (i.e., they are the four three-body interactions of the C atom, the Ni surface, and one of the H atoms). The final term in eq 1, Vbend, is the bend potential within the CHn group which largely controls the H-C-Ni and H-C-H angles. Each term in eq 1 is itself composed of several functions. We begin with some background information in section 3.1. In section 3.2, we present EAM7, which was mentioned in the Introduction and which enters into terms 2-5 of eq 1. In particular, we outline two changes made to the EAM6 H/Ni potential function26 which correct minor flaws in that potential. In section 3.3, we discuss the Vmet and VDIM terms in eq 1; we first present the functions required to introduce the carbon atom into the EAM and EDIM formalisms, and then we use these functions to specify the Vmet and VDIM terms. In section 3.4, we present the Vbend term in eq 1. In presenting Vbend, we discuss the changes and additions made to the previous DIM CH5 potential function28,29 to allow permutational symmetry to be enforced in the new potential function as well as the method by which we introduce the symmetrized potential into the embedded atom framework. Finally, in section 3.5, we summarize the fitting procedures which we used to optimize the PEF. Throughout this section, in order to keep it short, we omit most of the motivations for the functional forms; this material may be found in Appendix 2 in the supporting information. 3.1. Strategy. We distinguish three types of geometries. One set of geometries includes configurations where CH4 or

CH3 interacts strongly with the metal. This set includes the chemisorption configurations, and we call this the short-range interaction region. The adsorption energies, adsorption geometries, and vibrational frequencies in this region are fit both to experimental results of Ceyer and co-workers8 and also to theoretical findings by Yang and Whitten18,19 and Schu¨le, Siegbahn, and Wahlgren.12 The second set of geometries corresponds to physisorption of CH4 on the surface; we call this the moderate-range interaction region. Physisorbed molecules have not been observed experimentally, but Yang and Whitten20 have made calculations for physisorbed CH4, and we attempt to make our PEF qualitatively fit their results. The third set of geometries corresponds to the long-range dispersion region. An effort is made to insure that the long-range dispersion forces42 predicted by the new PEF are qualitatively correct. Note that there are no specific boundaries that separate the three regions; the division into three regions only applies to the strategy for fitting the PEF. We sometimes combine functions that are appropriate for different regions of the PEF by using Shepard interpolation43,44 or quintic splines. 3.2. Modification of the H/Ni Potential: EAM7. Reference 26 presents a general potential energy function for H on and in Ni. That potential function, called EAM6, is accurate for H adsorbed on the Ni surface, for interstitial H, and for various transition states connecting subsurface and on-surface sites. However, there are two problems with EAM6 which, although they do not affect any results presented in ref 26, need to be corrected in order to proceed with the current study. We call the corrected version of the H/Ni potential function EAM7, and we present the corrections in this section. The first problem with EAM6 is that the equation for the effective charge of the H atom (labeled ZH) gives an unphysical value of zero at the nucleus of the H atom (the proper value at the nucleus is 1.0). A consequence of this is that EAM6 predicts a favorable energy when an H atom occupies the same space as a Ni atom. To correct this, the effective charge for the H atom was modified to level off at 1.0 at the nucleus of the H atom by using a switching function. The EAM7 equation for the H atom effective charge as a function of the distance R from the nucleus is given as

{ [( )

ZH(R) )

R e RSW 1.0 + RZR5 + βZR4 + γZR3 3 R cj -R (2) exp + d1 sech [d2(R - d3)] R > RSW aj j)1 bj

∑

( )]

where RSW is chosen to equal 1.09 Å so that the fluctuation of eq 2 is small. The equation for the region of R greater than RSW is identical to that for EAM6, and the parameters for this part of the expression are given in ref 26. The parameters RZ, βZ, and γZ, given in Table 1, are chosen to make eq 2 continuous through the second derivative. Since eq 2 only differs from EAM6 in the region of R less than RSW, and since the H atom



never gets closer than 1.25 Å to a Ni atom in the work presented in ref 26, none of the results in that paper are affected as a result of this improvement. The second problem with EAM6 is in the moderate- and longrange density regions of the potential function, for which the treatment of the H atom in EAM6 is arbitrary. Preliminary calculations indicated that the embedding energy of the H atom had to be altered in the moderate-range interaction region (corresponding to physisorbed CH4) to obtain accurate energetics. The H atom embedding energy, as defined in ref 26, can be written as

FH(FjH,Fˆ (1) ˆ (2) H ,F H )

)

FPH(FjH)

+

FNL ˆ (1) ˆ (2) H (F H ,F H )

(3)

where FPH (the primary embedding energy) is a function of the sum of the electron densities, FjH, of all metal atoms in the system at the point at which the H atom is being embedded, and FNL H (the nonlocal embedding energy) is a function of two “powerˆ (2) treated” electron densities,26 Fˆ (1) H and F H . The power treatment indirectly accounts for the coordination number of the adsorbed atom and simulates gradient-corrected density functional effects,45,46 but the specifics of FNL H are beyond the scope of this paper, so we refer readers to ref 26 for more details on the FNL H term. In the current work, only FPH needed to be altered; FNL H remains unchanged from EAM6. To obtain the proper energetic behavior in the moderate- and long-range density regions, FPH was redefined for EAM7 as follows:

{

FPH(FjH) ) 1 RH1FjγH j H) H exp(-βH1F

0 e FjH e Fcut3 - ∆3

AH1(∆F5) + BH1(∆F5) + CH1(∆F5) 5

4

3

+ DH1(∆F6)5 + EH1(∆F6)4 + FH1(∆F6)3 Fcut3 - ∆3 e FjH e Fcut3

2 RH2FjγH j H) H exp(-βH2F

Fcut3 e FjH

a

i

ni

ℵi

ζi, Å-1

1 2 3 4 5 6 7 8 9 10

1 1 2 2 2 2 2 2 2 2

-0.208 14 -0.010 71 0.080 99 0.750 45 0.335 49 -0.147 65 0.282 41 0.546 97 0.231 95 0.010 25

10.272 73 17.919 79 1.998 41 2.880 51 5.072 78 7.938 82 1.853 35 2.728 08 4.914 35 12.302 41

See eqs 5-7. Ns ) 2 and Np ) 2.

with the surface, the EAM formalism requires three functions. Specifically, we need expressions for the electronic density of a C atom, the energy to embed a C atom into the electronic density created by other atoms, and a function representing the effective charge of the C core at various distances from the nucleus. Then, to make EDIM calculations involving multiple atoms interacting with the metal, singlet and triplet energy interaction functions involving the C atom are needed. We begin by describing the three EAM functions in section 3.3.1. We then describe the additional functions and parameters required by EDIM in section 3.3.2. Finally, in sections 3.3.3 and 3.3.4, we present Vmet and VDIM from eq 1 as functions of the EAM and EDIM functions presented in sections 3.3.1 and 3.3.2. 3.3.1. EAM Functions for Carbon-Nickel Interaction. The first EAM function needed is the electron density of a C atom. The 1s density of a C atom is significant only at distances less than 0.5 Å from the C nucleus; such distances are well below the smallest interatomic distances involving C atoms encountered in this study, so the density, Fa(R), of a single C atom at a distance R away from its nucleus is taken as the sum of the densities of the 2s and 2p orbitals. The expression for the total density has the same functional form as used by Daw and Baskes,32 namely,

(4)

Fa(R) ) NsFas (R) + NpFap(R)

(4a)

where Ns and Np are the effective number of s and p electrons for the C atom. The individual spherically averaged atomic densities of each orbital, Fas (R) and Fap(R), are taken from the restricted Hartree-Fock (RHF) calculations of Clementi and Roetti,47 which yield

where

∆F5 ) FjH - Fcut3

TABLE 2: Parameters in the Electron Density of a C Atoma

and

(∑ ) 10

Fa(s,p)(R) ) ∆F6 ) FjH - (Fcut3 - ∆3)

3.3. Incorporating Carbon and Hydrocarbons into the EAM and EDIM Formalisms; Vmet and VDIM. In this section, we present the functions that are required to introduce the carbon atom into the calculations. First, to model a C atom interacting

2

ℵiRi(R) /4π

(6)

i)1

(4b)

The final values of parameters defined by eq 4 are given in Table 1. The main difference between EAM7 and EAM6 for FPH is for FjH less than Fcut3, which is well below the density into which single H atoms adsorbed are embedded on or in Ni. For FjH greater than Fcut3, the EAM7 FPH is also refit but it equals the EAM6 FPH from ref 26 within 0.1%, and therefore the current PEF is very similar to EAM6 for H atom adsorbed on a Ni surface.

(5)

where i varies from 1 to 6 for the 2s orbital and from 7 to 10 for the 2p orbital, and

Ri(R) )

(2ζi)(ni+1/2) [(2ni)!](1/2)

R(ni - 1) exp(-ζiR)

(7)

The values for ℵi and ζi are taken from Clementi and Roetti.47 These values, along with the values of Ns, Np, and ni are given in Table 2. The next function needed to incorporate a C atom into the embedded atom framework is the energy to embed it in the electron density of the metal. Let FjC denote the electron density of the metal evaluated at the center of the C atom. The



TABLE 3: Parameter Set for the Potential Energy Function Parameters for embedding energy of C Atoma βC1 ) -16.04 Å-3 γC1 ) 0.4155 βC2 ) -35.42 Å-3 γC2 ) 0.5224 βC3 ) 28.54 Å-3 γC3 ) 0.8279 BC1 ) 8.340 × 1012 eV Å12 CC1 ) 5.489 × 109 eV Å9 EC1 ) 8.124 × 1012 eV Å12 FC1 ) -5.501 × 109 eV Å9 BC2 ) 1.648 × 1010 eV Å12 CC2 ) 7.474 × 109 eV Å9 EC2 ) 1.628 × 1010 eV Å12 FC2 ) -7.642 × 109 eV Å9 ∆1 ) 0.001 Å-3 Fcut2 ) 0.022 Å-3

RC1 ) -6.470 eV Å3γC1 RC2 ) -11.236 eV Å3γC2 RC3 ) -119.242 eV Å3γC3 AC1 ) 3.356 × 1015 eV Å15 DC1 ) -3.224 × 1015 eV Å15 AC2 ) 9.559 × 1011 eV Å15 DC2 ) -9.270 × 1011 eV Å15 Fcut1 ) 0.003 Å-3

∆2) 0.007 Å-3

b

a1 ) 9.8039 [8.0346]

Parameters for Effective Charge of a C Atom a2 ) 3.3380 Å-1 [3.2721] b1 ) -133.2758 [-31.044]

DCH ) 4.8668 eV a ) 1.706 Å-1

Parameters for Singlet and Triplet C-H Interaction Energiesc ∆CH ) 0.270 Req CH ) 1.094 Å b ) 0.136 Å-1 c ) 20.79 Å-1

Ah ) 64.282 Å3

Parameters for Triplet C-Ni Interaction Energyd Bh ) 0.0069 Å-3 Ch ) 0.5382

R ) 1.464 Å-1

Parameters for Available Electronic Density of H Atomse 0 RC-H ) 1.774 Å-1 RNi-H ) 1.9 Å

b2 ) 14.0996 Å-1 [8.916]

0 RC-H ) 1.09 Å

Equationsf

Parameters for H-Ni Map FHeq1 ) 0.010 Å-3 FHeq2 ) 0.180 Å-3 B ) -0.1551 C ) 0.0675 E ) -8.707 × 1024 Å13 F ) 4.292 × 1018 Å10

K1 ) 8.0 K2 ) 1.0 A ) 1.3067 Å D ) 4.038 × 1030 Å16 Fcut4 ) 1.0 × 10-6 Å-3

G ) 11.0 Å

Parameters for C-Ni Map Equationsg ) 0.0 Å-3 ) 0.024 Å-3 ) 0.180 Å-3 B′ ) -0.4368 C′ ) -0.1816

K1 ) 1.0 K2 ) 4.0 K3 )1.0 A′ ) 0.3466 Å

FCeq1 FCeq2 FCeq3

Parameters for Reference Angle in Out-of-Plane Bend Potentialh Rcut ) 10.0 Å Bφ ) -2.670 × 10-4 Å-3 Rφ0 ) 1.094 Å CH3 R3 ) 0.142 Å-3 h∆ ) 1.192 eV rad-2

τ ) 109.47° Aφ ) 4.742 Å-1 3 -2 fCH ∆ ) 0.596 eV rad Rmid ) 3.175 Å

Rm ) 1.0 Å

β3 ) 0.307 Å

Parameters for Proximity Functioni

a

See eq 8. b See eq 9. The values in brackets are the original values of Strand and Bonham (given for comparison, not used in this work). c See eqs 11-13. d See eq 16. e See eq 19. f See eqs 27-30. g See eqs 31-35. h See eqs 43-47. i See eq 50.

{

embedding energy of the C atom, denoted FC(Fj), is given as FC(FjC) ) 1 RC1FjγC j C) C exp(-βC1F

0 e FjC e Fcut1 - ∆1

AC1(∆F1) + BC1(∆F1) + CC1(∆F1) 5

4

3

+ DC1(∆F2)5 + EC1(∆F2)4 + FC1(∆F2)3 Fcut1 - ∆1 e FjC e Fcut1

2 RC2FjγC j C) C exp(-βC2F

Fcut1 e FjC e Fcut2 - ∆2

AC2(∆F3) + BC2(∆F3) + CC2(∆F3) 5

4

3

+ DC2(∆F4)5 + EC2(∆F4)4 + FC2(∆F4)3 Fcut2 - ∆2 e FjC e Fcut2

RC3FjγCC3 exp(-βC3FjC)

FjC > Fcut2

(8)

where

∆F1 ) FjC - Fcut1

(8a)

∆F2 ) FjC - (Fcut1 - ∆1)

(8b)

∆F3 ) FjC - Fcut2

(8c)

∆F4 ) FjC - (Fcut2 - ∆2)

(8d)

and

Figure 6. Energy to embed a C atom into a surrounding electron density, F.

Equation 8 joins three primary functions with two six-term quintic splines with widths ∆1 and ∆2. For the C atom embedding energy, we did not need to exploit the nonlocal embedding energy techniques used for the H embedding energy in ref 26. We want parameters appropriate to the sp2 and sp3 valence states of C in CH3 or CH4, not for the s2p2 of pure C. Thus we optimized the embedding function to data for these species, and this parametrization must overcome the use of Ns ) Np ) 2 in eq 5. The optimized parameters for eq 8 are given in Table 3, the resulting function is plotted in Figure 6, and the parametrization is explained below. The third function needed to incorporate a C atom into the EAM calculations is the effective nuclear charge on the C atom,


Wonchoba and Truhlar is given as

RCH ) a + b({tanh[c(R h - Req CH)] + 1}/2)

Figure 7. Effective charge of the C atom as a function of the distance from the nucleus. The function used in the present EDIM model is compared to the eight-term function of standard Bonham, which (to plotting accuracy) is indistinguishable from the result obtained with the two longest-range terms over this domain of R.

which we represented by

ZC(R) ) a1 exp(-a2R) + b1R exp(-b2R)

(9)

Table 3 gives the values of a1, b1, a2, and b2 and compares them to those of Strand and Bonham.47 The effective nuclear charge for the C atom with the present values of the parameters is plotted in Figure 7, where it is compared to the fit of Strand and Bonham.47 The effective nuclear charge on the C atom is used in the Coulombic pair repulsion of the EAM energy. The general pairrepulsion energy between two atoms, i and j, is given as

φij(Rij) )

Zi(Rij)Zj(Rij) Rij

2 ESCH ) DCH{1 - exp[-RCH(RCH - Req CH)]} - DCH

(11)

and

ETCH )

(

where R h is the average of the four C-H bonds in methane. The Morse parameters, DCH and Req CH, and the parameters a, b, and c in the function RCH are taken from previous work.29,30 These values are given in Table 3. The remaining parameter in eq 13, ∆CH, is a Sato parameter,49 which is optimized for the current PEF and is given in Table 3. The singlet C-Ni interaction in the EAM formalism is given as S ) FC(FjC) + EC-Ni

)

1 1 - ∆CH 2 (D {1 + exp[-RCH(RCH - Req CH)]} 2 1 + ∆CH CH DCH) (12)

where RCH is a function that is calibrated such that the C-H stretching frequencies agree with experiment,48 in both the methyl and methane limits. As in previous work,29 this function

∑j φC-Ni (RC-Ni ) j

j

(14)

where FjC is defined above as the total density created by metal atoms at the C atom, and j runs over all the metal atoms. Note that for a single adatom interacting with a surface, the total potential energy reduces to the sum over the singlet interactions of the adatom with all the metal atoms.26,27,50,51 The triplet energy for the C-Ni interaction is given as T ) hC-Ni(FjC){-FC(FjC) + EC-Ni

∑j φC-Ni (RC-Ni )} j

j

(15)

where hC-Ni(FjC) is an empirical, geometry-dependent generalization of the Sato parameter and is given by

(10)

where Rij is the interatomic distance between atoms i and j, and Zi and Zj are the effective nuclear charges of the two atoms. The three carbon functions discussed above, inserted in the standard EAM formalism, would allow one to perform calculations for a single C atom adsorbed on a Ni surface. However, such calculations are not recommended since the parameters are not optimized for a system with that much open shell character. To perform calculations for a C atom and three or four H atoms, we need an extended formalism, which we describe next. 3.3.2. EDIM Functions for Carbon-Hydrogen-Nickel Interaction. In the EDIM formalism, there are four important interaction energy terms: singlet and triplet adatom-adatom interactions and singlet and triplet adatom-metal interactions. The singlet and triplet C-H interaction energies are given respectively by the following Morse and anti-Morse potentials:

(13)

hC-Ni(FjC) )

1 tanh[Ah(FjC - Bh)] + Ch 2

(16)

where Ah, Bh, and Ch are optimized for the current PEF and given in Table 3. Equipped with these functions, we can define Vmet and VDIM in eq 1. 3.3.3. Vmet. The energy due to the metals, Vmet, is given as

Vmet )

∑i Vi

(17)

where

Vi ) Fi(FFi ) +

1 2

φij(Rij) ∑ i*j

(18)

where i and j run over all metal atoms in the system, the embedding function Fi for a Ni atom, given in ref 26, is an attractive function of its argument, and FFi is the total density of all other atoms in the system (metal and nonmetal) that is aVailable to metal atom i. Note that in the original EDIM method27 the embedding energy of metal atom i is a function of the total electron density Fi of all other atoms (both metal and nonmetal) evaluated at the center of atom i, whereas the embedding energy of nonmetal atom k is a function of the electron density Fjk of all metal atoms evaluated at the center of k. (The overbar denotes that nonmetal atoms are not included.) We retain this distinction here except that Fk is replaced by FFk , which is still a sum of contributions from all metal and nonmetal atoms, butsas explained nextsis the “free” or “available” density, not just the sum of the raw densities.



The available density for a non-hydrogen atom is equal to the its actual density. The available density of hydrogen atom k at metal site i is defined to be

(those involving the surface are actually effective distances) are used to specify the interactions among the six entities:

FFki ) FH(Rki) × 0 exp[-RNi-H(RRki - RNi-H )] 0 0 exp[-RNi-H(RRk - RNi-H )] + exp[-RC-H(RC-Hi - RC-H )] i

(19) where FH(Rki) is the actual electron density of the H atom at a distance Rki from its nucleus, Rki is the distance between atoms k and i, RNi-H and RC-H are Morse range parameters for Ni0 0 and RC-H H50 and C-H30 interactions, respectively, and RNi-H are parameters that were empirically adjusted to ensure that unphysical geometries have unfavorable energies. The four parameters in eq 19 are given in Table 3. 3.3.4. VDIM. VDIM(C-H-Ni) is the lowest energy root of the secular equation

det|H - SVDIM| ) 0

(20)

where S, the generalized configurational overlap matrix for two atoms interacting with Ni, is approximated by neglecting orbital overlap and is given as27,35-37

S11 ) S22 ) 4; S12 ) -2

(21)

and H, the electronic Hamiltonian, is given as

H11 ) 4(QC-H + QC-Ni + QH-Ni + JC-H) 2(JC-Ni + JH-Ni) (22a) H22 ) 4(QC-H + QC-Ni + QH-Ni + JC-Ni + JH-Ni) 2JC-H (22b) H12 ) -2(QC-H + QC-Ni + QH-Ni + JC-H + JC-Ni + JH-Ni) (22c) where

1 QA-B ) (ESAB + ETAB) 2

(23)

1 JA-B ) (ESAB - ETAB) 2

(24)

and ESAB and ETAB are singlet and triplet potential curves introduced above. Each of terms two through five in eq 1 represents a diatomic C-H fragment interacting with the Ni surface. The lower-energy root, VDIM(C-H-Ni), of eq 20 reduces to the following expression:

R1 ) RC-E1

R6 ) RE1-E2

R11 ) RE1-E4

R2 ) RC-E3

R7 ) RE3-E2

R12 ) RE1-E5

R3 ) RC-E4

R8 ) RE4-E2

R13 ) RE3-E4

R4 ) RC-E5

R9 ) RE5-E2

R14 ) RE3-E5

R5 ) RC-E2

R10 ) RE1-E3

R15 ) RE4-E5

It is not known beforehand which entity is the surface, and those distances involving the surface entity are effective (not literal) interatomic distances; these effective distances are defined below. The introduction of these effective distances is critically connected to the way that we enforce permutational symmetry with respect to interchanging the hydrogens. As preparation for discussing permutational symmetry, section 3.4.1 presents the method by which the surface is treated equivalently to a pseudoatom which will be called M. In section 3.4.2 we introduce the method used to ensure permutational symmetry for the present PEF. Then in section 3.4.3 we present the full bend potential. 3.4.1. EffectiVe Atom-Surface Distances: The Map Equations. To treat all possible bonding interactions of the carbon atom on the same footing, we define an effective distance between the metal surface and each of the other entities. Since the surface, unlike the other entities, is not a localized point in space, the effective distance is defined by making an analogy between the bonding capacity of the density created by the surface at an atom and the bonding capacity of a pseudoatom M that is a specific interatomic distance from that atom. To do this, we map the metallic electron density FjC encountered by a C atom approaching the surface onto an effective C-M interatomic distance, and we map the metallic density FjH encountered by an H atom approaching the surface onto an effective H-M interatomic distance, where M denotes a pseudometal atom that “represents” the surface. The pseudoatom M should, in some sense, be analogous to a hydrogen atom to facilitate the treatment of permutational symmetry. The H-M map is given by the following set of expressions:

[

RH-M(FjH) ) N1(FjH)

2

∑ i)1

KiReff,Hi (FjH) (FjH - FHeqi)2

]

(27)

where

VDIM(C-H-Ni) ) QC-H + QC-Ni + QH-Ni -

x[(JC-Ni + JH-Ni)2 + JC-H2 - (JC-Ni + JH-Ni)JC-H]

(26)

2

3.4. Permutationally Invariant Mapping of Surface Density. In this section, we present Vbend, which is the final term in eq 1. For CH4 interacting with a metal surface, there are six entities that need to be considered. Specifically, they are the carbon atom, the four hydrogen atoms, and the surface. The carbon atom will be denoted C, and the five non-carbon entities (the four hydrogen atoms and the surface) will henceforth be denoted E1, E2, E3, E4, and E5. Fifteen distances

{KiReff,H (FjH)} ∑ i)1

N1(FjH) ) 1/

(25)

Reff,H1(FjH) )

{

i

(28)

DFj5H + EFj4H + Fj3H + G FjH e Fcut4 (29) FjH > Fcut4 A(FjBH - FjCH)

and

Reff,H2 (FjH) ) 0.7 Å

(30)



Values for all parameters in these equations are given in Table 3. The C-M map is given by the following set of expressions:

[

RC-M(FjC) ) N2

3

∑ i)1

KiReff,Ci (FjC) (FjC - FCeqi )2

]

TABLE 4: The 20 Combinations of Entities E1 and E2. H(1)-H(4) Represent the Four H Atoms, and S Represents the Surfacea combination

E1

E2

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

H(1) H(1) H(1) H(1) H(2) H(2) H(2) H(2) H(3) H(3) H(3) H(3) H(4) H(4) H(4) H(4) S S S S

H(2) H(3) H(4) S H(1) H(3) H(4) S H(1) H(2) H(4) S H(1) H(2) H(3) S H(1) H(2) H(3) H(4)

(31)

where 3

{KiReff,C (FjC)} ∑ i)1

N2(FjC) ) 1/

i

(32)

Reff,C1(Fj) ) Reff,H1(Fj)

(33)

Reff,C2(FjC) ) A′(FjB′ - FjC′)

(34)

Reff,C3(FjC) ) 1.1

(35)

Labels E3, E4, and E5 may be assigned arbitrarily to the remaining three entities.

Values for all parameters in these equations are given in Table 3. 3.4.2. Permutational Symmetry. The equation for Vbend results from two fundamental assumptions. The first assumption is that entity E2 is the furthest entity from the C atom (i.e., that R5 has the largest value among R1-R5). The second assumption is that entity E1 is the closest non-carbon entity to entity E2 (i.e., that R6 has the smallest value among R6-R9). The remaining three non-carbon entities, E3, E4, and E5, are treated equivalently to one another. For geometries close to the reaction path of the H + CH4 f H2 + CH3 reaction, there is no ambiguity as to which of the H atoms correspond to entities E1 and E2. In particular, the reactant H atom is the entity furthest from the C atom (and is therefore treated as entity E2), and the H atom being extracted from methane to form H2 is treated as entity E1. However, for more general geometries or for the case where one H is replaced by the metal surface, the identification of entities E1 and E2 is not as straightforward. To illustrate this difficulty, consider the approach of methane to the metallic surface. Since the orientation of the approaching molecule is not known beforehand, it is not known which entities should be identified as E1 and E2. Furthermore, the molecule may rotate with respect to the surface over the course of its approach to the surface, thus changing which entities are associated with E1 and E2. To make the PEF a continuous function of geometry in regions where such switches occur, we consider every possible labeling of the entities and weight each one appropriately. There are 20 combinations of entities E1 and E2 that need to be examined in the calculation of Vbend. These combinations are shown in Table 4. The total bend potential is a linear combination of the 20 energy calculations:

the labeling conventions for E1 and E2 listed above (ωi is discussed in detail below), and Wnorm is a normalization factor,

and

ωiVbend,i)/Wnorm ∑ i)1

20

Wnorm )

(36)

where Vbend,i is the unweighted bend potential energy of the system calculated with combination i, ωi is a weighting factor that is a measure of how closely combination i corresponds to

ωi ∑ i)1

(37)

This inclusion of all possible combinations of entities E1 and E2 allows for a smooth transition between conformations, and it allows the PEF to be applicable to a CH4 molecule approaching the surface in any orientation. As a consequence of using eq 36, the potential function does not preferentially single out one H atom as the reactive one; it has the correct permutational symmetry. The weighting factor, ωi, in eqs 36 and 37, is given as

ωi )

(RE1-E3,4,5)P(RC-E2)2P (RE1-E2)P

(38)

where P is a positive integer, and the average interatomic distance between entity E1 and entities E3, E4, and E5 for a given combination is defined by

RE1-E3,4,5 )

(RE1-E3 + RE1-E4 + RE1-E5) 3

(39)

A value of P ) 6 is large enough that the present PEF reproduces energies, frequencies, and rate constants for H abstraction from CH4 calculated by the original J1 potential function to within 0.3%. 3.4.3. Vbend. We now describe the calculation of Vbend, which is given as

Vbend )

20

Vbend ) (

a

1

3

kj[(θj - θ0j )2 + (θj - θ0j )6] + ∑ 2 j)1 1

6

kj[(Rj - R0j )2 + (Rj - R0j )6] + Vop ∑ 2 j)4

(40)

where θj is one of the three angles that are formed by E1, the carbon atom, and E3, E4, and E5 (i.e., angles E1-C-E3, E1-



C-E4, and E1-C-E5); Rj is one of the three angles that are formed by the carbon atom with two of E3, E4, and E5 (i.e., angles E3-C-E4, E3-C-E5, and E4-C-E5); and θ0j and R0j are equilibrium angles defined in ref 28. The out-of-plane bend potential, Vop, and the force constants, kj, in eq 40 are interpreted and parametrized slightly differently than in previous work28,29 in light of the permutational symmetry of the PEF discussed above. We define Vop and kj in the following paragraphs. In the present PEF, the out-of-plane bend potential is given as 3

Vop ) f∆(RC-E1)

3

∆2i + h∆(RC-E1)∑∆4i ∑ i)1 i)1

(41)

where f∆ and h∆ are force constants defined below, and the angles are defined by25

∆i ) cos-1

(

(r2 - r1) × (r3 - r1)‚ri

The force constants in Vbend are defined as

)

|(r2 - r1) × (r3 - r1)|ri

- φ0(RC-E1)

(42)

k1 ) k0(RC-E1)F1F2

(48a)

k2 ) k0(RC-E1)F1F3

(48b)

k3 ) k0(RC-E1)F1F4

(48c)

k4 ) k0(RC-E1)F2F3

(43d)

k5 ) k0(RC-E1)F2F4

(48e)

k6 ) k0(RC-E1)F3F4

(48f)

and

where k0(RC-E1) is defined in ref 29, and the attenuation terms, F1 through F4, are given by

F1 ) [1 - P(RC-E2) exp(-a1RE1-E22)] × e )2)] × exp{-[a2 + a3 P(RC-E2) exp(-a4(RE1-E2 - RH-H e )2} (49a) (RC-E2 - RC-H

where r1, r2, and r3 are the three vectors associated with RC-E3, RC-E4, and RC-E5, and φ0 is a reference angle that defines an equilibrium value for the three angles formed by the three vectors and the C-E1 vector. This reference angle is defined below. Note that each of the 20 expressions (see eq 36) for Vop, will have a different set of values for r1, r2, r3, and RC-E1 depending upon which entities are chosen as E1 and E2. The reference angle is the substantial feature to change in the definition of VOP from ref 29. In the present PEF, the equation for the reference angle is taken as

φ0(RC-E1) ) τ + (τ - 90.0)[Sφ(RC-E1) - 1.0] RC-E1 e Rcut (43) RC-E1 > Rcut 90.0

{

Sφ(RC-E1) ) 1.0 - tanh[Aφ(RC-E1 -

3 exp(BφRC-E1 )]

(44) The out-of-plane bend force constants, f∆ and h∆, are given as

f∆(RC-E1) ) [1 - S3(RC-E1)]

e )2)] × exp{-[a2 + a3 exp(-a4(RE3-E2 - RH-H e )2} (49b) (RC-E3 - RC-H

F3 ) [1 - exp(-a1RE4-E22)] × e )2)] × exp{-[a2 + a3 exp(-a4(RE4-E2 - RH-H e )2} (49c) (RC-E4 - RC-H

and

F4 ) [1 - exp(-a1RE5-E22)] × e )2)] × exp{-[a2 + a3 exp(-a4(RE5-E2 - RH-H e )2} (49d) (RC-E5 - RC-H

where τ is the tetrahedral angle of 109.47°, and

Rφ0 )

F2 ) [1 - exp(-a1RE3-E22)] ×

3 fCH ∆

(45)

3 h∆(RC-E1) ) [1 - S3(RC-E1)]hCH ∆

(46)

and

where

S3(RC-E1) ) 1.0 tanh[R3(RC-E1 ) - Rφ0 )(RC-E1 - β3)2] (47) The parameters in eqs 43 through 47 are given in Table 3.

e e where a1, a2, a3, a4, RH-H , and RC-H are given in ref 28. In eq 49a, the factor P(RC-E2) is a measure of the proximity of entity E2 to the carbon atom and is given by

1 P(RC-E2) ) {1 - tanh[(RC-E2 - Rmid/Rm]} 2

(50)

where Rmid and Rm are parameters. (See Supporting Information for justification.) Since distortion of methane has been postulated4 to play a major role in the CH4 chemisorption process, we performed a validation test of the bending potential of methane against ab initio electronic structure calculations. We carried out a series of calculations for C3V methane in which the H-1 and C atoms define the symmetry axis, called Z. Then we varied the angle made by the other C-H bonds with this axis, decreasing this angle from the tetrahedral value of 109.47° down to very small angles corresponding to squashing all four H atoms on the same side of C. The ab initio calculations were carried out by Møller-Plesset second-order perturbation theory with the 6-311G** basis set,51 using the GAUSSIAN program.52 The comparison is given in Table 5. The present potential tracks the ab initio one quite well up to very severely distorted geometries, which gives us some confidence that the individual



TABLE 5: Distortion Energies (kcal/mol) for C3W Methanea

a

θ (deg)

MP2/6-311G**

present

109.47 105 100 95 90 80 70 60 50

0.0 1.4 5.9 13.2 22.5 45.1 70.8 102.2 151.8

0.0 1.3 5.2 10.8 17.3 32.7 55.0 91.8 149.5

The C-H bond lengths are frozen at 2.0673 a0.

terms in our fit are reasonably physical and that the need to introduce free density (see section 3.3.3) is not an artifact of an insufficiently repulsive bend potential. 3.5. Optimization Procedures. The optimization of the PEF was performed in several stages. In the first step, we adjusted the parameters such that the PEF yields binding energies, vibrational frequencies, and equilibrium geometries that are as close as possible to the available experimentally measured4,8 and theoretically calculated12,18,19 values for those quantities in the short-range interaction region of the potential. Most of the available experimental vibrational frequencies are for this region, specifically for CH3 adsorbed on Ni(111). During this stage of the fitting process, we were not concerned in detail with the behavior of the potential in the medium- and long-range interaction regions, because we intended to fine-tune the PEF in these regions by using piecewise potentials and Shepard interpolations as discussed previously. Also, since the information about chemisorbed configurations is known to a much greater quantitative accuracy than is information about physisorbed states or long-range interactions, and since the results produced by the PEF are very sensitive to choice of parameters, we wanted to isolate the fit to the short-range chemisorption region in order to produce as accurate a potential as possible for this region without adding additional constraints for behavior in other regions. As in previous work,26 the initial fitting to quantities in the chemisorption region involved minimizing an error function. The error function we minimized was

S)

(

∑R WR

)

Rcalc - Rlit Rlit

2

(51)

where R is an observable quantity, Rcalc is the value calculated for that quantity, Rlit is the literature value for that quantity, and WR is a weighting function which accounts for the quantitative differences in the levels of uncertainty among the observables. The unitless weighting function is given as

WR )

( ) Rlit ∆Rlit

2

(52)

where ∆Rlit is nominally the reported error in the quantity. There are two exceptions to this, however. In cases where no error is reported, which is often the case for frequencies, then the ∆Rlit value is empirically chosen to yield reasonable results in the optimization. Also, when the minimization of eq 52 results in localizing the error substantially on one particular quantity rather than spreading out the error somewhat evenly, then the ∆Rlit value for the offending observable is reduced. Table 6 gives the values of Rlit and ∆Rlit used in the minimization of eq 52. We used a genetic algorithm routine53 to optimize the parameter sets for the new functions introduced into the EAM formalism for the C atom and both the C and H map functions.

TABLE 6: Quantities Used in the Minimization of the Error Function, eq 53, for the Parametrization of the Potential Energy Functiona R

Rlit

∆Rlit

WR

E (3-fold eclipsed)b θ (3-fold eclipsed)d Zc (3-fold eclipsed)e C-H asymmetric stretch (E) C-H symmetric stretch (A1) asymmetric deformation (HCH bending modes) (E) symmetric deformation (umbrella mode) (A1) frustrated rotation (rocking modes) (E) methyl torsion (A2) C-Ni stretch (A1) frustrated translation (E)

-45.6 112 1.86 2730 2655 1320

c c c f f f

7.0 2 0.03 100 100 100

42.44 3136.0 3844.0 745.29 704.90 174.24

1230

f

100

151.29

965

f

100

93.12

485 385 160

f f f

100 100 100

23.52 14.82 2.56

ref

a The first three rows refer to CH3 adsorbed at the 3-fold eclipsed site. The remaining rows are frequencies (in cm-1) for CH3 adsorbed at the 3-fold staggered site; symmetry of the vibrational modes are given in parentheses. b E is the adsorption energy of CH3 in kcal/mol with respect to CH3 infinitely separated from the surface. c Schu¨le et al., ref 12; Yang and Whitten, ref 18; and Yang and Whitten, ref 19. d θ is the average of the angles that the three C-H bonds make with the C-Ni surface normal. e Zc is the C atom equilibrium height above the surface in angstroms. f Yang et al., ref 8.

We used all recommended53 niching and population sizing parameters in the genetic algorithm routine. When the genetic algorithm was unable to further minimize eq 53 in a reasonable amount of time, we assumed that we were in the vicinity of a global minimum, and we used a simplex algorithm54 routine to fine-tune the parameter set and isolate the global minimum. After we obtained an accurate fit to available quantities for chemisorption configurations, we further adjusted the PEF for physisorption configurations and long-range dispersion interaction. In general, this involved ensuring that the shape of the potential energy curves for CH3 and CH4 approaching Ni were physically accurate. In particular, it can be shown42 that when a Lennard-Jones gas approaches a metal surface at distances above the surface (z) greater than about 8 times the radius of an atom in the solid, the dispersion energy should enter as approximately (-C/z3), where C can be thought of as a LennardJones dispersion constant. For distances above the surface less than about 4 times the diameter of an atom in the solid, the dispersion energy is better approximated42 as (-C/z4). Also, the potential energy curve for CH4 approaching Ni should have a weakly bound physisorbed state before entering the repulsion region. The precise magnitude and location of this physisorbed state is not known experimentally, but we can obtain approximate qualitatively accurate behavior for the physisorption region based on theoretical calculations.20 We adjusted the PEF manually in the moderate- and long-range physisorption and dispersion regions to obtain appropriate results. So as not to damage the fit obtained for the chemisorption region, this stage of the fitting procedure thus resulted in the multipart functional forms for the hydrogen and carbon embedding energies (eqs 4 and 8) and the map equations (eqs 27-35). The final parameter set used for the current potential function is given in Tables 1-3. 4. Results The energetic results obtained with the present potential energy functions are presented in Tables 7-13. Tables 7-10 are for CH3 interacting with a Ni(111) surface. In Tables 7 and 9, all CH3 radicals are adsorbed in the C60 orientation,



TABLE 7: Energetic and Geometric Results for CH3 and H Interacting with Ni(111)a moving latticeb

rigid lattice CH4(g) CH3(g) + H(g) CH3(g) + H@HMN CH3(g) + H@GHM CH3@HIN + H(g) CH3@HIN + H adsorbed far awayc CH3@HIN + H@BGH (3rd adjacent) CH3@HIN + H@GHM (2nd adjacent)d CH3@HIN + H@HMN (adjacent)d CH3@HMN + H(g) CH3@HMN + H adsorbed far awayc CH3@HMN + H@BCH (3rd adjacent) CH3@HMN + H@BGH (2nd adjacent)d CH3@HMN + H@GHM (1st adjacent)d

E

E - EM

ZC

E

E - EM

ZC

-112.21 0.00 -63.49 -63.51 -43.81 -107.30 -107.54 -107.55 -106.27 -43.81 -107.30 -107.54 -107.54 -106.52

0.00 112.21 48.72 48.70 68.40 4.91 4.67 4.66 5.94 68.40 4.91 4.67 4.67 5.69

∞ ∞ ∞ ∞ 1.77 1.77 1.76 1.76 1.75 1.77 1.77 1.76 1.76 1.76

-112.21 0.00 -63.47 -63.49 -43.81 -107.32 -107.55 -107.55 -106.52 -43.82 -107.32 -107.55 -107.54 -106.50

0.0 112.21 48.74 48.72 68.40 4.89 4.66 4.66 5.69 68.39 4.89 4.66 4.67 5.71

∞ ∞ ∞ ∞ 1.76 1.76 1.75 1.74 1.76 1.76 1.76 1.75 1.75 1.76

a Energies in kcal/mol and distances in angstroms. b The zero of energy is the relaxed moving lattice for these values. The energy decreases 0.16 kcal/mol when the surface relaxes. c The distant site is an HMN site; i.e., the adsorbate is over a second-plane vacancy. d For H adsorbed at an adjacent or second adjacent site, the HIN site is slightly favored (by 0.02 kcal/mol).

TABLE 8: Comparison of Rigid-Lattice Results from Other Theoretical Calculationsa To the Present Results for Adsorption Energy of CH3 on Ni(111)a

CE configuration

C30 configuration C60 configuration B0 configuration A0 configuration

E

E - ECE

-32c -45 to -50d -38.7e

0 0 0

-36.5e -34.6e -35.5e -43 to -48d -33.9e

ZC

θ

C-H stretchb

C-Ni stretch

Literature

2.2e 1d 4.1e 3.8f 3.2e 2d 4.8e 5.7f 19g

1.91d 1.86e 1.56f 1.86e

112d 109.5e 117f 109.5e

1.86e 1.94f 1.98e 1.96d 2.03e

109.5e

2935d 2966e

353d 386e

2974d

This Work, Fixed CH3h 0.00 1.74

109.5e

296e 3168d

109.5e

CE configuration

-40.05

C30 configuration

-41.08

-1.03

1.74

109.5

C60 configuration

-41.46

-1.41

1.74

109.5

B0 configuration

-41.01

-0.96

1.77

109.5

A0 configuration

-39.69

0.36

1.95

109.5

CE configurationk

-42.40

This Work, Optimized CH3j 0.00 1.75 114

C60 configurationk

-43.81

-1.41

1.77

113

B0 configuration

-42.83

-0.43

1.78

113

A0 configuration

-41.62

0.78

1.97

110

109.5

416e

2753 (2)i 2623 2742 (2) 2611 2753 (2) 2607 2746 (2) 2614 2767 (2) 2724 2669, 2663 2564 2645, 2645 2515 2672, 2658 2555 2642, 2621 2617

505 491 462 479 526

425 416 400 460

a Energies, E, are given in kcal/mol, carbon atom equilibrium heights above the surface, Z , are given in angstroms, H-C-Ni angles, θ, are C given in degrees, and frequencies are given in cm-1. b The literature does not specify which C-H stretch frequency (symmetric or asymmetric) is being reported. For the current results, all C-H stretch frequencies are given. The nearly degenerate asymmetric stretch frequencies are given first, followed by the nondegenerate symmetric stretch. c Shustorovich, ref 10. d Schu¨le, Siegbahn, and Wahlgren, ref 12. e Yang and Whitten, refs 18 and 19. f Kratzer, Hammer, and Nørskov, ref 24. g Gavin, Reutt, and Muetterties, ref 9. h This set of results is obtained by using the same optimization procedures as Yang and Whitten18,19 (see text for details). i For fixed CH3, the result tabulated is the average of the two highest-frequency modes. j This set of results is obtained using full optimization of the CH adsorbate at each configuration. The C30 configuration site is not a stationary 3 point for this potential function, and therefore it does not optimize. k HMN site.

whereas in Tables 8 and 10 the orientations are specified. Tables 11 and 12 are for CH3 interacting with the (100) and (110) surfaces, respectively. Table 13 is for physisorption of CH4 on all three faces. Table 14 is for transition states. Minimum-energy sites (MESs) are identified by the absence of

imaginary frequencies. The ability to characterize stationary points by frequency analysis is one of the advantages of having an analytic potential energy function with analytic gradients. The zero of energy for all energies involving three hydrogen atoms is taken as the energy of CH3 infinitely far from the



TABLE 9: Moving-Lattice Results for CH3 Adsorbed in the C60 Orientation (Staggered 3-fold) on Ni(111) As Calculated in the Current Study and Compared to Experiment quantitya

experimentb

H@∞c

H@distantd

H@BCHe

H@BGHe

H@GHMe

C-H asymmetric stretch (E) C-H symmetric stretch (A1) asymmetric deformation (HCH bending modes) (E) symmetric deformation (umbrella mode) (A1) frustrated rotation (rocking modes) (E) methyl torsion (A2) C-Ni stretch (A1) frustrated translation (E) E E - EHg θ ZC

2730 2655 1320

2653, 2636 2533 1410, 1408

2653, 2636 2533 1410, 1408

2623, 2622 2515 1390,1390

2616, 2614 2508 1385, 1384

2644, 2620 2366 1394, 1389

1230

1572

1572

1504

1492

1488

965

964, 961

485 385 160f

235 454 79, 79 -43.82 19.65 113 1.76

964, 961 235 454 79, 79 -107.32 -43.85 113 1.76

892, 892 216 440 82, 80 -107.55 -44.08 113 1.75

867, 866 212 436 86, 79 -107.54 -44.07 112 1.75

880, 861 189 495 198, 74 -106.50 -43.03 113 1.76

a See caption of Table 4 for description of terms. In all cases, CH is adsorbed at site HMN. b Ref 8. c Results for CH adsorbed on a clean 3 3 Ni(111) surface. d Results for CH3 and H adsorbed at distant 3-fold sites. e BCH, BGH, and GHM are respectively third adjacent, second adjacent, and first adjacent sites. f The frequency of the doubly degenerate frustrated translation was only very tentatively assigned in ref 8, because the HREELS features in this range were obscured by surface phonon modes, which are also in this frequency range. g EH is defined in section 4; it equals -63.47 kcal/mol.

TABLE 10: Comparison of Vibrational Frequencies, Adsorption Energies, and Geometries of the CE, C60, B0, and A0 configurations for Rigid and Moving Ni(111) Surfaces. Degree of Degeneracy for Frequencies Is Given in Parenthesesa CEb

C60b

B0c

A0d

quantitya

rigid

moving

rigid

moving

rigid

moving

rigid

moving

C-H asymmetric stretches

2669 2663 2564 1387 1381 1591 943 924 144i 425 45i 20i -42.40 114 1.75

2667 2665 2565 1385 1379 1592 946 928 152i 461 52 12 -42.42 114 1.73

2645 2645 2515 1408 1408 1576 961 961 224 416 82 82 -43.81 113 1.77

2653 2636 2533 1410 1408 1572 964 901 235 454 79 79 -43.82 113 1.76

2672 2658 2555 1408 1370 1562 912 910 69i 400 89 27i -42.83 113 1.78

2678 2669 2568 1392 1380 1580 957 895 133i 443 62 44i -42.41 113 1.75

2642 2621 2617 1446 1442 1465 1029 1015 55i 460 53i 49i -41.62 110 1.97

2643 2631 2729 1444 1443 1490 1017 1011 61i 478 23i 13i -42.43 110 1.95

C-H symmetric stretch asymmetric deformations symmetric deformation frustrated rotations methyl torsion C-Ni stretch frustrated translations E θ ZC

a See caption of Table 4 for description of terms. b Adsorption site is HMN. This is the minimum-energy site (MES). c Has two imaginaryfrequency modes: one for lateral motion to a 3-fold site and one for rotation to a lower-energy B30 orientation. d Has three imaginary-frequency modes: two for lateral motion and one for rotation.

TABLE 11: Vibrational Frequencies (cm-1), Binding Energies (E, kcal/mol), C Atom Equilibrium Heights above the Surface Plane (ZC, Å), and H-C-Ni Surface Plane Angles (θ, deg) for CH3 Adsorbed at Various Sites on Rigid and Moving Ni(100) atop (Figure 4a)

bridgea (Figure 4b)

4-fold (Figure 4c)

quantity

rigid

moving

rigid

moving

rigid

moving

C-H asymmetric stretch C-H symmetric stretch H-C-H bend umbrella mode rocking mode methyl torsion C-Ni stretch frustrated translation Eb ZC θc

2663, 2659 2623 1416, 1404 1425 835, 831 77i 392 39i, 32i -42.71 1.98 109

2701, 2669 2599 1400, 1387 1453 757, 752 l55i 541 22i, 20i -40.47 2.00 109

2675, 2666 2559 1408, 1377 1535 929, 821 114 446 53, 44 -42.94 1.72 114(2), 111

2698, 2660 2559 1410, 1370 1565 975, 865 103 456 51, 36 -42.90 1.72 115(2), 111

2676, 2655 2572 1388, 1355 1580 916, 896 126i 445 52, 15i -46.38 1.60 117(2), 115

2671, 2666 2579 1384, 1346 1572 908, 865 169i 409 81i, 75i -41.40 1.59 117(2), 115

a MES. b Rigid adsorption energies are given with respect to the energy of a CH molecule infinitely separated from a rigid surface. Moving 3 adsorption energies are given with respect to the energy of a CH3 molecule infinitely separated from a lattice after the primary zone atoms have relaxed. b In cases where the equilibrium geometry of the adsorbed CH3 molecule is skewed with respect to the Ni surface plane resulting in more than one H-C-Ni surface plane angle, all angles are given, and the number of occurrences is given in parentheses.

surface. The zero of energy for all energies involving four hydrogens is CH3 and H both infinitely far from the surface and also infinitely separated from each other. Table 7 also gives

some energies relative to EM, which is the energy of a methane molecule infinitely far from the surface. Table 8 also gives some energies relative to ECE, which is the energy of a CH3



TABLE 12: Same as Table 10, Except for the Ni(110) Surface atop (Figure 5a)

bridge (Figure 5b)

3-folda (Figure 5c)

quantity

rigid

moving

rigid

moving

rigid

moving

C-H asymmetric stretch C-H symmetric stretch H-C-H bend umbrella mode rocking mode methyl torsion C-Ni stretch frustrated translation Eb ZC θc

2625, 2605 2477 1413, 1389 1583 935, 867 60i 773 35i, 18 -39.72 1.98 109(2), 104

2629, 2607 2478 1410, 1388 1571 912, 848 60i 749 36i, 27 -39.78 1.92 109(2),104

2619, 2609 2027 1378, 1369 1591 1111, 871 110 336 84, 18i -42.64 1.71 115(2),104

2611, 2600 2026 1379, 1365 1581 1109, 873 121 337 94, 14i -42.65 1.68 115(2), 104

2660, 2652 2248 1404, 1370 1625 1110, 843 46 420 104, 34 -42.70 1.68 118(2), 99

2648, 2640 2125 1406, 1371 1623 1111, 850 68 387 107, 38 -42.73 1.65 117(2), 110

a MES. b Rigid adsorption energies are given with respect to the energy of a CH3 molecule infinitely separated from a rigid surface. Moving adsorption energies are given with respect to the energy of a CH3 molecule infinitely separated from a lattice after the primary zone atoms have relaxed. c In cases where the equilibrium geometry of the adsorbed CH3 molecule is skewed with respect to the Ni surface plane resulting in more than one H-C-Ni surface plane angle, all angles are given.

TABLE 13: Selected Frequencies (cm-1), Physisorption Energies (E, kcal/mol), and C Atom Equilibrium Heights above the Surface Plane (ZC, Å) for CH4 Physisorbed on the Atop Sites of Rigid and Moving Ni(111), Ni(100), and Ni(110). Energies Are Given with Respect to CH4 Infinitely Separated from the Surface CH4/Ni(111)a

a

CH4/Ni(100)b

quantity

rigid

moving

rigid

moving

rigid

moving

methane-Ni stretch rocking modes E ZC

42 110, 104 -1.18 3.26

44 107, 106 -1.18 3.25

41 109, 100 -1.16 3.21

41 109, 100 -1.16 3.19

39 107, 91 -1.10 3.14

39 110, 91 -1.10 3.10

C above atop site and H atoms in same orientation as C60 configuration. b C above atop site; negligible barrier to rotating methyl group.

TABLE 14: Saddle Point Calculations for CH4 Chemisorption on Rigid Ni(111)a reactant: CH4 at infinity transition state product: CH3 at HIN (C60), H distant a

CH4/Ni(110)b

E

E - EM

ZC

-112.21 -98.36 -107.30

0.00 13.85 4.91

∞ 2.01 1.77

E in kcal/mol, ZC in Å.

adsorbed at site CE. Table 9 also gives energies relative to EH, which is the energy of CH3 infinitely separated from the surface and H adsorbed at site HMN. For calculations involving rigid surfaces, the zero of energy and EM, EEC, and EH are computed with rigid surfaces, and for calculations in which Ni atoms are movable, these Ni atoms were also relaxed in computing the zero of energy and EM, EEC, and EH. One can put these results on the same scale, if desired, by using the following fact: if six surface atoms and three subsurface atoms are allowed to relax in an isolated Ni(111) surface, with all possible adsorbates infinitely removed, the energy decreases by 0.16 kcal/mol. Thus the difference in the zero of energy for rigid and nonrigid results is 0.16 kcal/mol. 5. Discussion 5.1. Chemisorbed Configurations: Short-Range Interaction. In this subsection, we examine the three low-index faces in turn. 5.1.1. Ni(111). Several theoretical groups have previously examined CH3 adsorption on the Ni(111) surface, with widely varying results concerning the minimum-energy adsorption site.9,12,17-19,23,24 Gavin et al.,9 Schu¨le et al.,12 Yang and Whitten,18,19 and Kratzer et al.24 all found the minimum-energy site to be the 3-fold adsorption site. Schu¨le et al. calculated a more favorable energy for the 3-fold site than for the atop site, but the uncertainty of their calculations exceeded the calculated difference between the two sites, so they were unable to draw

a conclusion for the most favorable site based solely on energy calculations. Based on comparison of the experimentally measured4 C-Ni stretch frequency to their calculated C-Ni stretch frequencies for the atop and 3-fold site, however, they concluded that the 3-fold site was the favored adsorption site. On the other hand, de Koster and van Santen17 and Burghgraef et al.23 found the atop site to be the preferred site of CH3 adsorption on Ni(111). Ceyer and co-workers8 have used high-resolution electron energy loss spectroscopy (HREELS) to assign most of the vibrational frequencies of CH3 adsorbed on Ni(111). They were also able to determine, from the number of dipole active modes and the magnitude of the C-Ni stretch frequency, that the CH3 adsorbate must be at a site with C3V symmetry; this indicates that the H atoms are either eclipsed or staggered with respect to the three closest surface Ni atoms. In terms of Figure 1, then, the experiments are consistent with structures CE or C60, but not C30. Since the experimental evidence points to a preference for a CE or C60 3-fold site rather than a C30 3-fold or bridge or atop site, the question remains which 3-fold orientation (CE or C60) is preferred by the CH3 molecule. The theroetical groups12,18,19,24 who have examined this question conclude that the eclipsed orientation is preferred. Schu¨le et al.12 state that the barrier to rotation for the CH3 molecule in the 3-fold site is about 1 kcal/mol, in favor of the CE eclipsed orientation. Kratzer et al.24 also found the eclipsed orientation to be favored, but by 3.8 kcal/mol. Yang and Whitten concluded that the CE orientation is preferred by 4.1 kcal/mol from the following analysis. Yang and Whitten assumed, based on their study of NH3 on Ni(111),14 that the 3-fold site with the eclipsed H atoms (orientation CE) is the minimum-energy site of adsorption for CH3 on Ni(111). For this orientation, they froze the angle between the C-H bonds and the surface normal at 109.5° and the C-H bond distances at 1.07 Å, and they optimized the height of the C atom above the surface. Next, they varied the

6856 J. Phys. Chem. B, Vol. 102, No. 35, 1998 H-C-Ni surface angle from 75° to 150° and found the energy to be minimized and essentially unchanged for angles from 109.5° to 115°. Then, using a frozen tetrahedral CH3 geometry, they rotated the CH3 molecule 30° and 60°, keeping the C height above the surface plane constant, and calculated the adsorption energies of all of these orientations. Following this procedure using many of the early prototype potential functions over the course of the optimization, we also find the eclipsed site to be lower in energy than the staggered site, usually by about 1 kcal/ mol. (We note that with the final potential function, which is geared toward optimizing frequencies for the staggered orientation, the staggered site is lower in energy even using the somewhat eclipsed-biased optimization discussed above.) However, if we allow the adsorbed CH3 molecule to fully optimize at both the CE and C60 orientations independently, we always found that the energies were much closer to each other, and usually the staggered site is lower in energy (more so in the final potential function than in the early prototypes). Allowing lattice atoms to move increases the stability of the staggered site with respect to the eclipsed site even further. We conclude that the optimization of both the adsorbed molecule and the lattice atoms is very important in making a definitive decision about the preferred orientation of the CH3 molecule in the 3-fold site, so the results of Yang and Whitten alone do not indicate conclusively which orientation is preferred. In an attempt to make a more definitive assignment of the minimum energy configuration of adsorbed CH3, we consider the methyl torsion frequency. Ceyer and coworkers8 measured this frequency to be 485 cm-1. They indicated that although this frequency appears somewhat high for the hindered rotation around the C-Ni surface normal, it could be an indicator of significant interaction between the H atoms and the surface and a possible high potential barrier to rotation. Over the course of our optimization, although we were never able to obtain a methyl torsion frequency as high as 485 cm-1, the frequency was in all cases much higher (i.e., closer to 485 cm-1) for the staggered orientation than for the eclipsed orientation. In fact, in most cases the methyl torsion at the eclipsed configuration was in the range from 10 to 40 cm-1, and it was often imaginary, indicating that the eclipsed configuration was a saddle point. Kratzer et al.24 calculated a methyl torsion frequency of 129 cm-1 for CH3 in the CE configuration, which is a higher frequency than we were able to obtain for this configuration but is still very low compared to the experiment.8 Unfortunately, Yang and Whitten18,19 and Schu¨le et al.12 did not calculate full Hessians at the stationary points, and therefore it is not known whether the eclipsed site is a saddle point or a local minimum in their work. The only frequencies reported by Yang and Whitten and Schu¨le et al. are the C-H stretch and the C-Ni stretch (the latter of which are in very good agreement with experiment in both cases, further verifying the 3-fold site as the observed site). Additionally, in most cases over the course of our optimization, the staggered site was the preferred site energetically. Therefore we suggest, based on our optimizations, for which the methyl torsion for the staggered orientation was always in far better agreement with experiment than it was for the eclipsed orientation, that the observed orientation for CH3 adsorbed in the 3-fold site is with the H atoms staggered relative to the three nearest Ni atoms (orientation C60). Upon allowing the system to optimize to this orientation, we were able to obtain much better agreement with all the experimentally observed frequencies than we were able to obtain when we tried to force the eclipsed configuration to be the preferred orientation.

Wonchoba and Truhlar Next we compare our highly parametrized semiempirical model to ab initio theory. The first set of results listed in Table 8 gives the theoretically calculated value of the quantity from the literature, the second set of results is the quantity calculated by the current potential function using the partial geometry optimization procedure used by Yang and Whitten,18,19 and the third set gives the quantity calculated by the current potential function using a full optimization of the CH3 molecule. All quantities in this table are calculated with no moving Ni atoms in order to provide a consistent comparison to literature values. Note that since the second set of results are forced into a particular geometry that is not a minimum, they are not stationary points. Therefore, we present these results not as our best estimates to the quantities, or even reasonably accurate estimates to the quantities, but rather as a measure of comparison to other theoretical work. These results in Table 8 are indeed in reasonably good agreement with previous theoretical work. The third set of results are, however, our best estimates because they involve full minimization of the CH3 adsorbate at each site. Another interesting topic is the behavior of the C-H stretching vibration. Yang and Whitten18,19 calculated the frequency of the C-H stretch at the 3-fold site to be about 300 cm-1 higher than the experimental value,4 and they performed additional calculations to examine this discrepancy further. Specifically, they examined the C-H stretch frequency for a tilted CH3 radical (i.e., with a geometry in which two of the H-C-Ni surface angles are more than 109.5, and the other is less than 109.5) adsorbed on Ni(111) slightly shifted from the 3-fold site so that one of the H atoms is directly above a Ni atom. For this geometry, they calculated a frequency that is much closer to the experimental value. However, the spectroscopic measurements of Ceyer and co-workers,8 which show that the adsorbed species is of C3V symmetry, eliminate the tilted and shifted orientation as a possible adsorption site. When we use the same optimization route as Yang and Whitten (see the second set of results in Table 8), we calculate the C-H stretch to be higher than when we fully optimize the CH3 radical (see the third set of results in Table 8 and also Table 9). Thus, it is possible that the geometry of the 3-fold adsorbed CH3 radical was not fully optimized in the calculations of Yang and Whitten, and a more complete optimization of the geometry might bring the calculated C-H stretch frequency more in line with experiment. Table 9 gives results calculated for CH3 adsorbed in the staggered 3-fold orientation on a moving Ni(111) lattice and compares the results to experiment8 where available and to other theoretical calculations. These results represent our best estimates of these quantities. Again our calculated frequencies are in good agreement with experiment. In an attempt to examine the effects of coadsorbed H on the surface, we calculated another set of results for an H atom and the CH3 molecule adsorbed at adjacent 3-fold sites. In most cases, the frequencies go down, most notably for the frustrated rotations and the symmetric deformation. The symmetric deformation mode is in much better agreement with experiment when H is coadsorbed. These results are given in Table 9 as well. As discussed in ref 23, because the Ni-Ni interaction parameters were very carefully optimized to reproduce the bulk lattice constant, we can test the effect of assuming a rigid lattice. Some such effects can be seen by comparing Tables 8 and 9. Table 10 gives a more complete comparison of structural and energetic quantities plus all frequencies for the CE, C60, B0, and A0 configurations. The results in Table 10 indicate that



Figure 9. Potential energy diagram of CH4 approaching atop sites on Ni(111), Ni(100), and Ni(1 10). Figure 8. CH4 physisorbed on Ni(111). The C atom is directly above an atop site. Only primary zone atoms are shown.

lattice motion does not greatly affect the structures, energies, and frequencies of the adsorbed CH3 configurations. A similar conclusion was drawn23 for H adsorbed on Ni(111). Yang et al.21 calculated the effect of a subsurface hydrogen atom on the chemisorption energy of CH3 at a 3-fold site and found, assuming a bulk geometry for the Ni atoms, that the chemisorption energy increases from 38 to 47 kcal/mol. In the present calculations, with nine Ni atoms optimized in the presence of the adsorbate (CH3) and absorbate (H), we find a much smaller effect; in particular the chemisorption energy decreases from 43.8 to 44.8 kcal/mol. 5.1.2. Ni(100) and Ni(110). Most previous work done with CH3 on Ni has been for the (111) surface. However, our use of effective H-Ni surface and C-Ni surface distances based on density allows us to make calculations for CH3 adsorbed on any interface. Here we examine CH3 adsorbed on the other low-index faces, Ni(100) and Ni(110). Tables 11 and 12 present binding energies, geometries, and vibrational frequencies for CH3 adsorbed at various sites on rigid and moving Ni(100) and Ni(110). For the moving lattice calculations, the atoms allowed to move are those shown in Figures 4 and 5 for each adsorption site. From our calculations, the minimum energy adsorption site for CH3 on Ni(100) (Figure 4b) has the C atom 1.72 Å above a bridge site between two Ni atoms. One H atom is pointing toward a 4-fold vacancy, and the other two H atoms are staggered with respect to the first, pointing toward 2-fold bridge sites adjacent to the one that the C atom is above. The CH3 molecule is tilted with respect to the surface plane such that the angle formed by the H atom which points toward the 4-fold vacancy, the C atom, and the Ni surface normal is 110°. The other two H-C-Ni surface normal angles are 115°. The minimum-energy adsorption site for CH3 on Ni(110) (Figure 5c) has the C atom bound 1.65 Å above the surface plane in a pseudo-3-fold site formed by two adjacent surface plane Ni atoms and one second plane Ni atom. The H atoms are eclipsed toward the three Ni atoms which form the adsorption site. The molecule is tilted with respect to the surface plane such that the angle formed by the H atom which is eclipsed with the second plane Ni atom, the C atom, and the Ni surface normal is 100°. The other two H-C-Ni surface normal angles are 117°. 5.2. Physisorbed Configurations and Long-Range Interactions. Although no experimental data are available for physisorbed configurations of CH4 on Ni, previous theoretical results on Ni(111) are available from Yang and Whitten.20 Yang and Whitten examined several approaches of methane towards an atop site on Ni(111). The most stable orientation that they found is shown in Figure 8, and it has three H atoms pointing toward the surface. They found that methane physisorbs with

Figure 10. Potential energy diagram of CH3 approaching an atop site on Ni(111) and a least-squares fit to the dispersion region (higher than 3 Å above the surface). The value of C is 1068.4 Å4 kcal/mol.

this orientation at an energy of -0.7 kcal/mol and with the C atom approximately 3.5 Å above the surface. Table 13 gives our results for methane physisorbed on an atop site of Ni(111) in the orientation shown in Figure 8. Results are given for a rigid lattice and for a lattice with 11 moving atoms. The results are in good qualitative agreement with those of Yang and Whitten. Table 13 also gives results for CH4 physisorbed on atop sites of Ni(100) and Ni(110) (again with three H atoms pointing down). The internal geometries and frequencies of the CH4 molecules are nearly identical to those of gas-phase methane. Therefore, Table 12 gives only the calculated results that are specific to the physisorbed molecule (i.e., proximity to the Ni surface, Ni-methane stretching frequencies, rocking vibrational frequencies, and physisorption energies). Comparison of the results calculated with rigid and moving Ni atoms in Table 13 indicates that surface motion does not substantially affect these quantities. Figure 9 gives potential energy curves for methane approaching the atop sites of all three low-index faces of Ni. All three curves are very similar, as are the results for each crystal face in Table 13, indicating that the physisorbed state is too far above the surface to distinguish the corrugation of each crystal face. Finally, we consider CH3 approaching Ni(111), but still far enough away that the interaction is weak. Figure 10 is a plot of the total potential energy of a CH3 radical approaching an atop site on Ni(111). We saw above from Figure 9 that when the molecule is far from the surface, the specific structure of the crystal face does not greatly affect the energy and geometry of the approaching molecule. As a result, in the dispersion region (when the C atom is greater than roughly 3.5-4.0 Å above the surface), the potential energy curve for CH3 approaching any other site would look very similar to that for the atop site approach shown in Figure 10. To test the qualitative accuracy of the potential energy curve in the dispersion region, we also plot in Figure 10 a least-squares fit to this curve of the functional form -C/z4, where z is the distance from the center of the C atom to the center of the Ni atom at the atop site toward


Figure 11. Top and side views of the transition state for dissociative chemisorption of CH4 on Ni(111). Only primary zone atoms are shown. The bond length of the breaking C-H bond is 2.593 ao (1.37 Å).

which the CH3 molecule is approaching. The fit, where C is equal to 1068 Å4 kcal/mol, closely matches the potential curve, showing that the potential varies approximately as -z-4, in accord with the theoretical treatment of Steele.42 5.3. Transition State for Dissociative Chemisorption. Since we did not use any transition-state data or kinetics information for species involving carbon to parametrize the present surface, such data provide a useful test. Most previous work has not involved full gradient-based geometry optimization, but this is possible with our present surface. We succeeded in optimizing the saddle point for dissociative chemisorption of CH4 to produce CH3 and H on rigid Ni(111), and we verified that this stationary point has only one imaginary frequency. For the present PEF, the dissociative chemisorption reaction (starting with CH4 in the gas phase) is 4.9 kcal/mol endoergic (zeropoint-energy-exclusive) to produce CH3 at HIN and H widely separated. The transition state, shown in Figure 11, has an energy of 13.8 kcal/mol with respect to CH4 infinitely far from the Ni(111) surface (i.e., the barrier for the back reaction is 8.9 kcal/mol starting from distant sites). We estimate a zero-pointenergy-inclusive forward barrier height of 10.5 kcal/mol. The classical and zero-point-energy-inclusive barrier heights of the forward transition state, 13.8 and 10.5 kcal/mol, respectively, bracket the values for the activation energy, 12-13 kcal/mol, that Beebe et al.31 estimated from the experiments of Ceyer et al.4 and from their own experiments. This is very encouraging

Wonchoba and Truhlar since our experience with these kinds of processes is that the phenomenological experimental activation energies often agree with one or another of these theoretical quantities within a couple of kcal/mol or better. For example, we previously found55,56 that the high-temperature activation energy for siteto-site hopping of H on Cu and Ni is very close (in both cases) to the corresponding classical barrier height. The present dissociative chemisorption reaction, however, may be more dominated by tunneling, and it is possible that the activation energy is considerably less than the barrier height. For example, Jansen and Burghgraef57a obtained a barrier of 22.3 kcal/mol by adjusting a three-dimensional wave packet simulation to agree with the experimental data of Ceyer and co-workers4 at a single energy, and by comparing to experiment at other energies they conclude that the true barrier is a few kcal/mol higher. They point out though that the 12 degrees of freedom of CH4 not included in their simulations will also have an effect on the estimation of the true barrier. They further point out that there is “a whole range of experimental4,31,57b,58 values for the barrier height,” from 7.9 to 31.1 kcal/mol. The geometry of the 13.8 kcal/mol transition state corresponds to the breaking of the C-H bond with the carbon over a bridge site (HN in Scheme 1) formed by two adjacent Ni atoms. This is a different transition state than the one reported by Kratzer et al.,24 who calculated a transition state corresponding to the breaking of a C-H bond with the carbon over an atop Ni site. For this state, those authors calculate a classical energy of approximately 26 kcal/mol with respect to CH4 infinitely separated from the surface. Using our PEF with a rigid Ni surface, we calculate a value of 26.0 kcal/mol for the precise transition-state geometry specified by those authors. This comparison for the atop site gives further confidence in the reasonableness of the present analytic potential energy function. Thus, we seem to have accomplished our goal of creating a PEF that is useful for studying dissociative chemisorption transition states as well as stable adsorbed states. In earlier calculations, Anderson and Mulroney11 and Yang and Whitten22 also estimated an activation energy for dissociative chemisorption over an atop site and obtained 15 and 17 kcal/mol, respectively, relative to CH4 at infinity. It is not clear why their values are lower than the values obtained by Kratzer et al. for an atop site. Perhaps our conclusions about the geometry of the transition state are more significant than the absolute barrier height. Yang and Whitten22 concluded that “the stability of H and CH3 at nearby sites on the Ni(111) surface strongly implicates the atop Ni site as the most suitable site for CH4 dissociation.”25 As a consequence they examined only the atop pathway. In the present work we used analytic gradients and Newton searches starting at a wide variety of initial guess geometries to look for the true lowest energy saddle point, and we found that dissociation over a bridge is the most favorable. Further discussion of the transition states is planned for a follow-up paper on the reaction dynamics of dissociative chemisorption. 6. Concluding Remarks We have presented an analytical potential energy function for methyl radical and methane interacting with a Ni surface. The new potential function presented here combines features of the J1 potential energy function29 for the reaction CH3 + H2 T CH4 + H, the EAM6 potential energy function26 for H on and in Ni, and the EDIM potential energy functions of previous work27 for H2 and H3 on Ni. The essential modification of

Embedded Diatomics-in-Molecules Potential the gas-phase reaction potential is to introduce an interaction of the covalent adsorbate with a metal by using an embedding function of the metal density. The essential modification of the embedded-atom model for the metal is to introduce nonmetallic atoms that exhibit valence saturation. The CH4/ Ni system is modeled by embedding the four CH subsystems of CH4 in the metal surface by converting the density of the surface at each adatom into an effective atom-surface distance, which has a length scale comparable to an H-H distance. The potential function is symmetrized with respect to permuting H atoms so any orientation of the adsorbed molecule on the surface is permitted, and there are no functional discontinuities between any two adsorption geometries. This permutational symmetry will be especially important in future work, in which the present potential can be used to examine the chemical dynamics of dissociative chemisorption and associative desorption. The analytical surface reproduces experimental8 measurements of the vibrational frequencies of the adsorbed methyl radical on clean Ni(111) with a mean unsigned deviation from experiment of 108 cm-1. When an H atom is adsorbed at an adjacent 3-fold site, the mean unsigned deviation increases to 134 cm-1, but the error in the symmetric deformation mode is reduced by 25%. Combining our calculations of the methyl torsion frequency for CH3 adsorbed in the 3-fold hollow site of Ni(111) in both the staggered and eclipsed orientations with the observation of a relatively high frequency for this mode (485 cm-1) in the experimental work indicates to us that the minimum energy orientation is the staggered orientation. This conclusion disagrees with the ab initio results of Yang and Whitten.18,19 A possible explanation for this disagreement is that in the calculations of Yang and Whitten, the local geometry of the CH3 molecule was partially optimized for the eclipsed formation, and the optimized geometry of the staggered orientation was taken as a simple rigid rotation of the CH3 molecule by 60° from the eclipsed orientation. Applying the same optimization procedures as used by Yang and Whitten to our analytical potential energy function yields results comparable to theirs; however, our calculations show that allowing the adsorbed molecule to fully optimize at each orientation changes the energy by a significant amount. Additionally, the calculations of Yang and Whitten were performed with a rigid metal lattice, and our calculations also indicate that relaxation of the lattice is important to fully optimize the geometry and energy of the adsorbed species. Acknowledgment. The authors are grateful to Gregory D. Hawkins and John C. Tully for helpful discussions. This work was supported in part by the National Science Foundation. Supporting Information Available: Appendix 2, containing extra discussion for Section 3 and justification of the functional forms used in the fit (11 pages). Ordering information is given on any current masthead page. Appendix 1. Primary Zones and Secondary Zones for All Sites The primary zone atoms are the atoms enclosed within a sphere centered at the approximate location of the C atom for a given orientation. These atoms are allowed to move during the calculations. The primary zone is taken to be large enough that adding additional atoms to the primary zone does not change the results. The secondary zone atoms are the atoms enclosed within a sphere centered at the approximate location of the C atom for

J. Phys. Chem. B, Vol. 102, No. 35, 1998 6859 TABLE 15: Lattices site

Nprim

Nsec

111, atop 111, bridge 111, center 100, atop 100, bridge 100, center 110, atop 110, bridge 110, center

9 9 10 9 8 9 12 10 10

677 680 660 579 596 599 576 584 582

a given orientation. The radius of the sphere enclosing the atoms in the secondary zone is 11 Å larger than the radius of the sphere enclosing the primary zone. Table 15 specifies the actual lattices used. The second column is the number of primary zone atoms for the movingsurface calculations. The third column is the number of secondary zone atoms that are not part of the primary zone. For the rigid-surface calculations, the primary zone atoms are moved into the secondary zone; that is, the radius of the secondary zone does not change. The secondary zone could be smaller for the rigid surface calculations, but the computational savings would be negligible. For the (111) surface the moving atoms are G, H, I, M, N, S, f, g, and k in the notation of Scheme 1. These are the six surface atoms and three second-plane atoms that are closest to site HMN. For methyl groups adsorbed at HIN in Table 7, we checked that the results do not change if, instead of these nine atoms, the six subsurface atoms and three subsurface atoms around the HIN site are allowed to move. References and Notes (1) (a) Probstein, R. F.; Hicks, R. E. Synthetic Fuels; McGraw-Hill: New York, 1982. (b) van Hook, J. P. Catal. ReV. Sci. Eng. 1980, 21, 1. (c) Rostrup-Nielsen, J. R. In Catalysis Science; Anderson, J. R., Boudart, M., Eds.; Springer: Berlin, 1984; Vol. 5. (2) (a) Satterfield, C. N. Heterogeneous Catalysis in Practice; McGrawHill: New York, 1980. (b) Bond, G. C. Heterogeneous Catalysis; Oxford University Press: Oxford, 1987. (3) Campbell, I. M. Catalysis at Surfaces; Chapman and Hill: London, 1988. (4) Lee, M. B.; Yang, Q. Y.; Ceyer, S. T. J. Chem. Phys. 1987, 87, 2724. (5) Ceyer, S. T. Langmuir 1990, 6, 82. (6) Ceyer, S. T. Science 1990, 249, 133. (7) Daley, S. P.; Utz, A. L.; Trautman, T. R; Ceyer, S. T. J. Am. Chem. Soc. 1994, 116, 6001. (8) Yang, Q. Y.; Maynard, K. J.; Johnson, A. D.; Ceyer, S. T. J. Chem. Phys. 1995, 102, 7734. (9) Gavin, R. M.; Reutt, J.; Muetterties, E. L. Proc. Natl. Acad. Sci U.S.A. 1981, 78, 3981. (10) Shustorovich, E. Surf. Sci. 1986, 176, L863. (11) Anderson, A. B.; Mulroney, J. J. J. Phys. Chem. 1988, 92, 809. (12) Schu¨le, J.; Siegbahn, P.; Wahlgren, U. J. Chem. Phys. 1988, 89, 6982. (13) Yang, H.; Whitten, J. L. J. Chem. Phys. 1989, 91, 126. (14) Chattopadhyay, A.; Yang, H.; Whitten, J. L. J. Phys. Chem. 1990, 94, 6379. (15) Bell, A. T.; Shustorovich, E. J. Catal. 1990, 121, 1. (16) Swang, O.; Faegri, K., Jr.; Gropen, O.; Wahlgren, U.; Siegbahn, P. Chem Phys. 1991, 156, 379. (17) de Koster, A.; van Santen, R. A. J. Catal. 1991, 127, 141. (18) Yang, H.; Whitten, J. L. J. Am. Chem. Soc. 1991, 113, 6442. (19) Yang, H.; Whitten, J. L. Surf. Sci. 1991, 255, 193. (20) Yang, H.; Whitten, J. L. J. Chem. Phys. 1992, 96, 5529. (21) Yang, H.; Whitten, J. L.; Thomas, R. E.; Rudder, R. A.; Markanus, R. J. Surf Sci. Lett. 1992, 277, L95. (22) Yang, H.; Whitten, J. L.; Markanus, R. J. Appl. Surf Sci. 1994, 75, 12. (23) Burghgraef, H.; Jansen, A. P. J.; van Santen, R. A. Faraday Discuss. 1993, 96, 337; J. Chem. Phys. 1994, 101, 11012; Surf. Sci. 1995, 324, 345. (24) Kratzer, P.; Hammer, B.; N0rskov, J. K. J. Chem. Phys. 1996, 105, 5595.

6860 J. Phys. Chem. B, Vol. 102, No. 35, 1998 (25) Whitten, J. L.; Yang, H. Surf. Sci. Reports 1996, 24, 55. (26) Wonchoba, S. E.; Truhlar, D. G. Phys. ReV. B 1996, 53, 11222. (27) (a) Truong, T. N.; Truhlar, D. G.; Garrett, B. C. J. Phys. Chem. 1989, 93, 8227. (b) Truong, T. N.; Truhlar, D. G. J. Phys. Chem. 1990, 94, 8262. (28) Steckler, R.; Dykema, K. J.; Brown, F. B.; Hancock, G. C.; Truhlar, D. G.; Valencich, T. J. Chem. Phys. 1987, 87, 7024. (29) Joseph, T.; Steckler, R.; Truhlar, D. G. J. Chem. Phys. 1987, 87, 7036. (30) Raff, L. M. J. Chem. Phys. 1974, 60, 2220. (31) Beebe, T. P., Jr.; Goodman, D. W.; Kay, B. D.; Yates, J. T., Jr. J. Chem. Phys. 1987, 87, 2305. (32) Daw, M. S.; Baskes, M. I. Phys. ReV. B 1984, 29, 6443. (33) Rice, B. M.; Garrett, B. C.; Koszykowski, M. L.; Foiles, S. M.; Daw, M. S. J. Chem. Phys. 1990, 92, 775; 1994, 100, 8556(E). (34) Daw, M. S.; Foiles, S. M.; Baskes, M. I. Mater. Sci. Rep. 1993, 9, 251. (35) London, F. Z. Elektrochem. Angew. Chem. 1929, 35, 552. (36) Eyring, H.; Polanyi, M. Z. Phys. Chem. 1931, B12, 279. (37) Eyring, H.; Walters, J.; Kimball, G. Quantum Chemistry; Wiley: New York, 1944. (38) Pola´k, R. Paidarova´, I.; Kuntz, P. J. Int. J. Quantum Chem. 1997, 62, 659. (39) Cashion, J. K.; Herschbach, J. Chem. Phys. 1964, 40, 2358. (40) Ellison, F. O. J. Am. Chem. Soc. 1963, 85, 3640. (41) Tully, J. C. In Semiempirical Methods of Electronic Structure Calculation, Part A; Segal, G. A., Ed.; Plenum: New York, 1977; p 173. (42) Steele, W. A. The Interaction of Gases with Solid Surfaces; Pergamon Press: New York, 1974. (43) Ischtwan, J.; Collins, M. J. Chem. Phys. 1994, 100, 8080. (44) Nguyen, K. A.; Rossi, I.; Truhlar, D. G. J. Chem. Phys. 1995, 103, 5522. (45) Kryachenko, E. S.; Ludena, E. V. Energy Density Functional Theory of Many-Electron Systems; Kluwer: Dordrecht, 1990.

Wonchoba and Truhlar (46) Perdew, J. P.; Burke, K.; Ernzerhof, M. ACS Symp. Ser. 1996, 629, 453. (47) Clementi, E.; Roetti, C. A. Data Nucl. Data Tables 1974, 14, 177. Strand, T. G.; Bonham, R. A. J. Chem. Phys. 1964, 40, 1686. (48) JANAF Thermochemical Tables, 2nd ed.; Natl. Stand. Ref. Data Ser. 37; Stull, D. R., Prophet, H., Eds.; National Bureau of Standards: Washington, DC, 1971. (49) Sato, S. Bull. Chem. Soc. Jpn. 1955, 28, 450; J. Chem. Phys. 1955, 23, 592. (50) Walch, S. P.; Bauschlischer, C. W., Jr. Chem. Phys. Lett. 1982, 86, 66. (51) Hehre, W. J.; Radom, L.; Schleyer, P. v. R.; Pople, J. A. Ab Initio Molecular Orbital Theory; Wiley: New York, 1986. (52) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Gill, P. M. W.; Johnson, B. G.; Robb, M. A.; Cheeseman, J. R.; Keith, T.; Petersson, G. A.; Montgomery, J. A.; Raghavachari, K.; Al-Laham, M. A.; Zakrzewski, V. G.; Ortiz, J. V.; Foresman, J. B.; Cioslowski, J.; Stefanov, B. B.; Nanayakkara, A.; Challacombe, M., Peng, C. Y.; Ayala, P. Y.; Chen, W.; Wong, M. W.; Andres, J. L.; Replogle, E. S.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Binkley, J. S.; Defrees, D. J.; Baker, J.; Stewart, J. P.; HeadGordon, M.; Gonzalez, C.; Pople, J. A. Gaussian 94, Revision D.4; Gaussian, Inc.: Pittsburgh, PA, 1995. (53) Carroll, D. L. AIAA J. 1996, 34, 338. (54) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes in FORTRAN, 2nd ed.; Cambridge University Press: New York, 1992. (55) Wonchoba, S. E.; Hu, W.-P.; Truhlar, D. G. Phys. ReV. B 1995, 51, 9985. (56) Wonchoba, S. E.; Truhlar, D. G. J. Chem. Phys. 1993, 99, 9637. (57) (a) Jansen, A. P. J.; Burghgraef, H. Surf. Sci. 1995, 344, 149. (b) Luntz, A. C. Private communications to the authors of ref 57a. (58) Hamza, A. V.; Madix, R. J. Surf. Sci. 1987, 179, 25.

Embedded Diatomics-in-Molecules Potential Energy Function for

Recommend Documents