The dissociative chemisorption of methane on Ni(100) and Ni(111): classical and quantum studies based on the reaction path Hamiltonian.

The dissociative chemisorption of methane on Ni(100) and Ni(111): Classical and quantum studies based on the reaction path Hamiltonian Michael Mastromatteo and Bret Jackson Citation: The Journal of Chemical Physics 139, 194701 (2013); doi: 10.1063/1.4829678 View online: http://dx.doi.org/10.1063/1.4829678 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/139/19?ver=pdfcov Published by the AIP Publishing Articles you may be interested in The dissociative chemisorption of methane on Ni(111): The effects of molecular vibration and lattice motion J. Chem. Phys. 138, 174705 (2013); 10.1063/1.4802008 The dissociative chemisorption of methane on Ni(100): Reaction path description of mode-selective chemistry J. Chem. Phys. 135, 114701 (2011); 10.1063/1.3634073 Quantum dynamics of dissociative chemisorption of CH 4 on Ni(111): Influence of the bending vibration J. Chem. Phys. 133, 144308 (2010); 10.1063/1.3491031 The temperature dependence of methane dissociation on Ni(111) and Pt(111): Mixed quantum-classical studies of the lattice response J. Chem. Phys. 132, 134702 (2010); 10.1063/1.3357415 Dissociative chemisorption of methane on Ni(100): Threshold energy from CH 4 (2ν 3 ) eigenstate-resolved sticking measurements J. Chem. Phys. 119, 6407 (2003); 10.1063/1.1613935

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 131.170.6.51 On: Tue, 11 Aug 2015 22:27:59

THE JOURNAL OF CHEMICAL PHYSICS 139, 194701 (2013)

The dissociative chemisorption of methane on Ni(100) and Ni(111): Classical and quantum studies based on the reaction path Hamiltonian Michael Mastromatteo and Bret Jacksona) Department of Chemistry, University of Massachusetts, Amherst, Massachusetts 01003, USA

(Received 19 September 2013; accepted 27 October 2013; published online 15 November 2013) Electronic structure methods based on density functional theory are used to construct a reaction path Hamiltonian for CH4 dissociation on the Ni(100) and Ni(111) surfaces. Both quantum and quasiclassical trajectory approaches are used to compute dissociative sticking probabilities, including all molecular degrees of freedom and the effects of lattice motion. Both approaches show a large enhancement in sticking when the incident molecule is vibrationally excited, and both can reproduce the mode specificity observed in experiments. However, the quasi-classical calculations significantly overestimate the ground state dissociative sticking at all energies, and the magnitude of the enhancement in sticking with vibrational excitation is much smaller than that computed using the quantum approach or observed in the experiments. The origin of this behavior is an unphysical flow of zero point energy from the nine normal vibrational modes into the reaction coordinate, giving large values for reaction at energies below the activation energy. Perturbative assumptions made in the quantum studies are shown to be accurate at all energies studied. © 2013 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4829678] I. INTRODUCTION

The rate-limiting step in steam reforming, the major industrial source for syngas and molecular hydrogen, is the dissociative chemisorption of methane on a Ni catalyst. As a result, there have been numerous molecular beam studies of this reaction, primarily on Ni and Pt surfaces,1–3 and it has become the prototypical polyatomic-surface reaction. The barriers to methane dissociation on smooth Ni and Pt surfaces are large,4 on the order of an eV, and a single C–H bond breaks as the molecule collides with the surface, leaving chemisorbed H and CH3 fragments. Molecular beam studies have shown that the zero-coverage sticking probability, S0 , increases with both the translational energy of the molecule normal to the surface and its vibrational energy, but that this behavior is nonstatistical.2, 3 For dissociation on Ni(100), for example, exciting the symmetric stretch of CH4 leads to a much larger increase in S0 than if the same amount of energy were put into translational motion.5 On the other hand, using this energy to excite the antisymmetric stretch leads to a slightly smaller increase in S0 than if it was used to increase translational motion.6 The ν 4 bending mode has been observed to be much less efficacious for promoting dissociation than either stretch.7 Similar mode-selective behavior has been observed on Ni(111) and other metal surfaces, though the details vary from surface-to-surface.2, 3 Bond-selective chemistry has also been observed. For example, in reactions of CD3 H on Ni(111), excitation of the v1 stretch preferentially breaks the C–H bond relative to a C–D bond by a 30:1 ratio.8 The dissociative chemisorption of methane on metals is also activated by the thermal energy of the lattice, and S0 has been a) Author to whom correspondence should be addressed. Electronic mail:

[email protected]

0021-9606/2013/139(19)/194701/9/$30.00

observed to increase substantially with increasing substrate temperature.9–16 A full theoretical understanding of these behaviors has been hampered by the complexity of the problem; it is difficult to construct an accurate potential energy surface (PES) that includes all or most of the 15 molecular degrees of freedom (DOFs), and it is impossible to implement a quantum reactive scattering calculation that treats all of these DOFs exactly. In addition, it is necessary to include the effects of lattice motion and substrate temperature for any direct comparison with experimental data. Until very recently, quantum studies of this problem have been limited to relatively lowdimensional models and/or approximate potentials.10, 11, 17–31 One way around the problem of constructing a full PES is to do ab initio molecular dynamics (AIMD). Unfortunately, a direct calculation of sticking in this manner may require a great many trajectories, as measured values for S0 are typically very small,2 below 10−2 . Also, this approach ignores quantum effects such as tunneling, though tunneling contributions to S0 may be less important at higher substrate temperatures.32 Sacchi et al. have used AIMD to study the recombinative desorption of methane from Pt(110)-(1×2)33, 34 and Ni(100),35 by starting trajectories at the transition state and integrating them away from the surface. While they cannot compute S0 in this manner, a vibrational analysis of the desorbed molecules combined with detailed balance arguments found behavior consistent with observed mode specificity. Recently, two high-dimensional quantum approaches have been used to explore the dynamics of methane dissociation. The first is based on the Reaction Path Hamiltonian36, 37 (RPH), in which the PES is assumed to be harmonic with respect to displacements away from the reaction (minimum energy) path. Given this approximation, it is easier to treat problems with many DOFs, and all of the terms defining the

139, 194701-1

© 2013 AIP Publishing LLC


194701-2

M. Mastromatteo and B. Jackson

Hamiltonian can be computed in a relatively straightforward manner using density functional theory (DFT).25 In studies on Ni(100),32 Ni(111),38 and Pt(110)-(1×2),39 quantum dissociation probabilities were computed using a 13-DOF model which assumed rotational adiabaticity, a rigid surface, and impact at a single surface site (the top site). Sudden models were then used to average over the remaining two molecular DOFs, corresponding to impact sites on the surface unit cell,32 as well as to introduce the effects of lattice motion.24, 26 Dissociative sticking probabilities were computed as a function of collision energy, molecular vibrational state, and substrate temperature, entirely from first principles. Comparison with experiment was good, both with regard to the magnitude of the reactivity, its variation with collision energy and substrate temperature, and the enhancement of this reactivity with vibrational excitation. More recently, Guo and coworkers40 developed a 12-DOF PES for methane dissociation on Ni(111), fitting to a very large number of DFT total energy calculations, again for the case of a fixed surface impact site (top site) and a rigid surface, but also assuming that the surface is flat. Using this PES, dissociation probabilities were computed quantum mechanically for a reduced 8-DOF model system, after making additional assumptions about the symmetry of the non-reacting methyl group. Surface impact site averaging and thermal lattice effects were then introduced using the same sudden approaches used in the RPH studies. The agreement with experiment was also very good, reproducing the vibrational efficacies observed in the experiments.40 This PES was then used in a 12-DOF quasi-classical trajectory (QCT) study of the dissociation of CH4 , CHD3 , and CH2 D2 on Ni(111).41 While impact site and lattice motion effects were not included, reaction probabilities at above threshold translational energies were compared for different initial vibrational states, and the trends were consistent with experimental data regarding vibrational efficacy. This behavior was further elucidated using the recently developed Sudden Vector Projection model.42 In this paper, we use both quantum and QCT methods with our reaction path PES to examine the dissociative chemisorption of CH4 on Ni(100) and Ni(111), with several goals in mind. First, we want to test the accuracy of the QCT approach for this reaction. It is now possible to fit force fields of very high dimensionality to thousands of DFT calculations, making full-dimensional QCT studies of this reaction likely in the near future. AIMD is another very powerful method for studying these complex systems, and there will be increased application of AIMD to this and similar hydrocarbon reactions on metals. It is essential to understand the limitations of QCT, and in particular, there may be problems with “incorrect” IVR (intramolecular vibrational redistribution), where the large amount of zero point energy (ZPE) in each of the nine vibrational modes of CH4 can improperly flow into other DOF. In this paper, both QCT and quantum RPH approaches will be applied to our rigid surface 13-DOF model systems. Well-tested sudden approximations will then be used to average these results over surface impact sites and include lattice motion effects, so that we can compute sticking probabilities over a wide range of energies, at experimental temperatures. The first question is to what extent the large amount of zero

J. Chem. Phys. 139, 194701 (2013)

point vibrational energy in the molecule will artificially enhance reactivity, particularly at low incident energies. Another question is how well the quantum and QCT results agree with respect to the variation in reaction probability with translational energy and vibrational excitation. Can QCT methods accurately reproduce the mode-selectivity observed in the experiments? Results from the Guo group41 suggest that the vibrational efficacies are lower in the QCT studies than observed in the experiments, but they did not compute the sticking, S0 , or make a direct QCT-quantum comparison. Finally, are these errors, as well as the neglect of tunneling, even important once thermal averaging over the lattice vibrations are included? It is hoped that these classical trajectory studies, in combination with our quantum approach, may shed additional light on the nature of mode-selective chemistry. In addition to addressing these questions, we have improved our quantum model, and present our new results for dissociation on both Ni(100) and (111) here for the first time. Specifically, we have increased our vibrational basis to include all 2-quanta states, and have correspondingly expanded our Hamiltonian to include terms that are higher order in the vibrationally non-adiabatic coupling. In this paper, we use QCT methods to test the accuracy of these perturbative assumptions. In Sec. II, we briefly describe our DFT calculations and present an overview of our RPH model and sudden approaches. In Sec. III, we present the results for both rigidsurface single-site reaction probabilities and fully averaged dissociative sticking probabilities, for CH4 on Ni(100) and Ni(111). QCT and fully quantum results are compared with each other and with existing experimental data. We close with some conclusions in Sec. IV.

II. THEORY A. Constructing the RPH

To construct our RPH, we use the Climbing ImageNudged Elastic Band method43, 44 to locate the reaction path, the minimum energy path from the reactant to the product configuration, passing over the transition state. Total energies along this path are computed using the DFT-based Vienna ab initio simulation package (VASP), developed at the Institut für Materialphysik of the Universität Wien.45–49 A supercell with periodic boundary conditions is used to represent our system as an infinite slab, with a large vacuum space above the slab to separate it from its repeated images. The interactions between the ionic cores and the electrons are described by fully nonlocal optimized projector augmentedwave (PAW) potentials,49, 50 and exchange-correlation effects are treated using the Perdew-Burke-Ernzerhof (PBE) functional.51, 52 Full details can be found in our earlier studies on Ni(100)4, 25, 32 and Ni(111).4, 38 The distance along the 15 (dxi )2 , and the xi are reaction path is s, where (ds)2 = i=1

the mass-weighted Cartesian coordinates of the CH4 nuclei. At numerous points along this path we compute the total energy, V0 (s), and perform a normal mode analysis. This gives the normal vibrational coordinates {Qk } and corresponding


194701-3


J. Chem. Phys. 139, 194701 (2013)

frequencies {ωk (s)}, k = 1-14, describing displacements orthogonal to the reaction path at point s, in the harmonic approximation. To construct the RPH, we ignore higher order (anharmonic) terms in the PES. The eigenvectors, Li,k (s), define the transformation between {xi } and our RPH coordinates s and {Qk } xi = ai (s) +

14

Li,k (s) Qk ,

(1)

k=1

where ai (s) is the configuration of the molecule at s, in terms of our mass-weighted Cartesian coordinates. Changing variables, the RPH has the form37 H = Hvib + V0 (s) +

1 (ps − πs )2 / (1 + bss )2 , 2

(2)

where Hvib =

14 1

1 Pk2 + ωk2 (s) Q2k , 2 2

k=1

πs =

14 14

Qk Pj Bk,j (s) ,

(3)

(4)

k=1 j =1

and bss =

14

Qk Bk,15 (s) .

(5)

k=1

The momenta conjugate to s and {Qk } are ps and {Pk }, respectively, and the vibrationally non-adiabatic couplings in the kinetic energy operator are given by Bk,j (s) =

15 dLi,k i=1

ds

Li,j (s) .

FIG. 1. Energies of the 14 normal modes of CH4 , along the reaction path s, for dissociation on Ni(100). The 8 vibrations with A symmetry are labeled 1 –8 , and the 6 modes with A symmetry are labeled 1 –6 . The lines are fits to a function used in the trajectory studies.

(6)

On both Ni(100) and Ni(111), the carbon is almost directly over the top site at the transition state, with the dissociating C– H bond angled towards the surface. The vibrational frequencies ¯ωk (s) describing motion orthogonal to the reaction path are plotted in Fig. 1 for CH4 on Ni(100). The points correspond to our VASP-based normal mode calculations, and the lines are fits to a 10-parameter function used to define the PES in the QCT studies. The fit is quite good. When the molecule is far above the surface (large negative s), 9 of the frequencies are non-zero: the triply degenerate antisymmetric stretch, ν 3 , the symmetric stretch, ν 1 , the doubly degenerate bend, ν 2 , and the triply degenerate bend, ν 4 . The interaction with the surface removes these degeneracies and softens the symmetric stretch and some of the bends near the transition state, s = 0. We find similar behavior on all of the surfaces we have explored on Ni and Pt.4, 32, 38, 39 Because of this softening, ZPE corrections lower the barrier heights by about 0.1 eV. The remaining 5 frequencies correspond to molecular rotation and translation parallel to the surface. These are zero, asymptotically, but as the molecule approaches the surface they become hindered rotations and translations. Our reaction path is symmetric with respect to reflection through a plane that lies perpendicular to the surface and along the dissociating C–H bond, and the normal modes are either symmetric (A ) or antisymmetric (A )

with respect to reflection through this plane. We label the A modes 1 –8 and the A modes 1 –6 . This symmetry is important, as the vibrationally non-adiabatic couplings, Bij , only couple modes of the same symmetry. In Fig. 2, we plot the curvature couplings, Bk,15 , for the A modes (Bk,15 = 0 for the A vibrations). Again, the points correspond to the calculations and the lines are fits to a 10-parameter function used in the QCT PES. While these couplings are particularly large for the bending vibrations, it is important to note that the efficacies for promoting vibration are typically larger for the stretching vibrations. The frequencies and couplings for CH4 dissociation on Ni(111) are very similar to those plotted in Figs. 1 and 2.

FIG. 2. Vibrationally non-adiabatic couplings, Bi,15 , along the reaction path s, for CH4 dissociation on Ni(100). The lines are fits to a function used in the trajectory studies.


194701-4


J. Chem. Phys. 139, 194701 (2013)

B. Quantum equations of motion

Our close-coupled wave packet RPH approach to reactive scattering is described elsewhere,25, 32 though we provide some details here. We write our total molecular wave function (t) = ψn (s, {Qk }u ; t)n ({Qk }b ; s), (7) n

where {Qk }b are the normal coordinates of the 9 modes that are bound at large negative s, and {Qk }u are the 5 that are unbound (ωk = 0). The n are eigenfunctions of that part of Hvib corresponding to the bound modes (1 –6 and 1 –3 in Fig. 1), bound with eigenvalues ¯ωk (s)[nk + 12 ] and quantum numbers k

n = {nk }b . We expand the quantum form of the RPH37 in a power series in bss . The operators bss and π s describe energy flow between the vibrationally adiabatic states; bss links states that differ by one vibrational quantum and π s couples states that differ by two quanta. The time-dependent Schrödinger equation then leads to coupled equations of motion32 for the ψ n (s, {Qk }u ; t). Two of these {Qk }u correspond to X and Y, the motion of the molecular center of mass parallel to the surface. It is reasonable to assume that this motion is slow on collision timescales, given the relatively large total molecular mass, the large collision energies, and our normal incidence conditions. Also, experiment suggests that this reaction is relatively insensitive to parallel motion of the molecule.53 We thus implement our scattering calculation for fixed values X0 and Y0 , and average the resulting reaction probabilities over all impact sites on the surface unit cell. The remaining {Qk }u describe molecular rotation. The rotational temperature in the molecular beams is about 10 K, and we assume that the molecule is initially in the ground rotational state. As the molecule approaches the metal, there are always one or more H atoms pointing towards the surface, and only minor angular reorientation is required to enter the transition state. As the moment of inertia is relatively small, and there is likely to be a reasonable amount of rotational steering, we assume rotational adiabaticity, writing ψn (s, {Qk }u ; t) = χn (s; X0, Y0 ; t)R0 ,

(8)

where R0 is the ground state eigenfunction of that part of the Hamiltonian containing the 3 remaining {Qk }u describing molecular rotation. Given the spherical shape of CH4 , one expects that the collision will not lead to significant rotational excitation, and experiment suggests that moderate rotational excitation of the incident molecule does little to modify the reactivity.54 At t = 0, consists of a gaussian wave packet χ n0 (s; X0 , Y0 ; 0) in some initial vibrational channel n0 , far from the surface. At later times, the interaction with the surface generates wave packets on the other vibrationally adiabatic channels n. These wave packets evolve on the effective potentials 1 , ¯ωk (s) nk + Veff,n (s) = V0 (s) + 2 k=1 14

(9)

and the parametric dependence of n on s leads to derivative coupling terms, so that curve crossing becomes increas-

ingly likely at higher velocities as well as for larger Bij . Standard techniques32, 38 are used to propagate the coupled wave packets, and the reactive flux at large positive s is Fourier transformed in time on each channel n, giving both stateresolved and energy-resolved reaction probabilities for all incident energies Ei included in the initial wave packet.55, 56 In practice, we implement this rigid-surface, rotationally adiabatic, fixed-site calculation for X0 and Y0 corresponding to the top site, where the barrier is lowest. This gives us a reaction probability P0 (Ei , n0 ). We then average over impact sites, approximating the reaction probability at other sites by shifting P0 (Ei , n0 ) along the energy axis by an amount equal to the change in barrier height.32 This has worked well in earlier studies,32, 38, 39 including the recent work of Guo and coworkers.40 We then introduce the effects of lattice motion using well-tested models.24, 26 DFT-based methods are used to examine how the location and height of the barrier changes with the motion of the lattice atom over which the methane dissociates. On all Ni and Pt surfaces we have examined, the barrier height can decrease significantly as this metal atom moves out of the surface plane,4 leading to a strong variation in dissociative sticking with substrate temperature. Sudden models are used to average over all displacements of the lattice atom, for a substrate temperature Ts . The final result is the dissociative sticking probability, S0 (Ei , n0 , Ts ). Finally, we note that in our earlier studies, the sum over n in Eq. (7) included only the ground state and the 9 singly excited vibrational states. Consistent with this, we truncated our Hamiltonian to first order in bss and π s . In this paper, we now also include in all 45 doubly excited states of the molecule. Correspondingly, we expand H to include all 2-quanta terms: second order in bss and first order in π s . While these improvements allow us to describe 2-quanta initial states, we find that they do not appreciably modify our earlier results for the ground or 1-quantum initial states.

C. QCT methods

The QCT studies are implemented in standard fashion. We use the same rigid-lattice, rotationally adiabatic, fixedsite Hamiltonian used in our quantum studies. Thus, as for the quantum studies, the initial molecules are assumed to be in the rotational ground state. To compute P0 (Ei , n0 ), we typically run 10 000 trajectories for each Ei and n0 . For reaction probabilities below 0.005, we run 100 000 trajectories. We randomly sample the initial phases of the normal modes, with quantized energies, including ZPE. These results are then averaged over surface impact sites and lattice motion, using the same methods as for the quantum results, to get S0 (Ei , n0 , Ts ). We note that the RPH, Eq. (2), has a singularity at bss = −1. In practice, this is not a problem. For CH4 dissociation on Ni(111) and Ni(100), the most trajectories we have had to abort for a particular Ei is 1.6% and 4.7%, respectively. As this happens only at higher incident energies, where P0 (Ei , n0 ) is large, this does not add appreciable error to our calculations. In order to test the perturbative approximations made in our quantum formalism, we also consider the case where we


194701-5


J. Chem. Phys. 139, 194701 (2013)

expand the denominator in Eq. (2) to low order. In particular, we write the last term in Eq. (2) as 1 (ps − πs )2 (1 − bss )2 , 2

(10)

where we have approximated (1 + bss )−1 ≈ (1 − bss ). We note that this is not the same result that a Taylor series expansion of Eq. (2) would give. However, the Eq. (10) form ensures that this kinetic energy operator remains ≥0, while the Taylor series form does not. Consequently, this form always leads to stable trajectories while the Taylor series form does not. In any case, as we show in Sec. III, this perturbative form gives very similar results to those using the full Hamiltonian. These observations are consistent with our quantum studies, 2 term does not significantly modify where addition of the bss the results. III. RESULTS AND DISCUSSION

In Fig. 3, we plot QCT results for the single-impact-site rigid-lattice reaction probability P0 (Ei , n0 ), for CH4 incident on Ni(111). The molecule is initially in the ground vibrational state or one of the three vibrationally excited states indicated. Results are also plotted for the classical perturbative approximation discussed in Sec. II C. First, we see that vibrational excitation enhances reactivity for all three modes, as expected for a late barrier system. While the 1ν 1 and 1ν 3 states are nearly degenerate, the 1ν 1 excitation is clearly more effective at enhancing dissociation. While no experimental data exist for molecules in the 1ν 1 state, this increased efficacy relative to 1ν 3 is observed on Ni(100).5, 6 Second, we see that the perturbative approximation produces results very similar to the non-perturbative case, with the largest error occurring at higher energies, as one might expect. A concern with QCT and AIMD approaches is that classical mechanics tends to overestimate the amount of IVR from

FIG. 3. Classical single-site static lattice reaction probabilities P0 as a function of incident energy, for CH4 dissociation on Ni(111). Results are plotted for methane initially in the ground (black), 1ν 1 (blue), 1ν 3 (red), and 1ν 4 (brown) vibrational states. Results are shown using both the full Hamiltonian (solid lines) and the perturbative form discussed in the text (dashed).

an initially excited vibration onto the other DOFs. In addition, there is a significant amount of ZPE in CH4 , 1.2 eV, and classical mechanics allows this energy to flow unphysically into the other molecular DOFs. This can be difficult to monitor during the reaction, as the Hamiltonian is non-separable, but we can compare with quantum results, and also look for reactions that clearly violate quantum ZPE restrictions. In Fig. 3, we note the activation energy, Eact , equal to the barrier height with all ZPE corrections. For molecules in the vibrational and rotational ground state, reaction for Ei < Eact should only occur via tunneling, and our classical P0 would optimally be zero at these energies. Instead, we see that the ground state P0 is relatively large in this region, and equal to 0.34 for both the full and perturbative calculations at Ei = Eact . In Fig. 4, we replot the ground and 1ν 3 data on a semi-log plot, to better enhance this low-probability region, and also to better compare with the quantum RPH results, which are also plotted. Results are also shown for the 2ν 3 state, an average of P0 over the six 2-quanta combinations of the 1 , 2 , and 1 modes in Fig. 1. Below threshold, the quantum reaction probabilities drop exponentially, while our QCT reaction probabilities remain large. This is most problematic for the ground state, where P0 can be sizable in the tunneling region. Once the translational energy plus any initial vibrational excitation is greater than or equal to Eact , the results are better, but for ground state CH4 this is very high energy. However, for the initially excited states, where there actually is a significant amount of vibrational energy transfer into the reaction coordinate, the quantum and classical results are in reasonable agreement over a broad range of energy. If we average the results of Fig. 4 over surface impact sites and introduce the effects of lattice motion, we get the dissociative sticking probabilities, S0 , plotted in Fig. 5. There are several important results. First, our quantum RPH

FIG. 4. Single-site static lattice reaction probabilities P0 as a function of incident energy, for CH4 dissociation on Ni(111). Results are plotted for methane initially in the ground (black), 1ν 3 (red), and 2ν 3 (green) vibrational states. QCT results are shown for both the full Hamiltonian (thin solid lines) and the perturbative form discussed in the text (dashed). Quantum results are plotted as thick solid lines.


194701-6


FIG. 5. Dissociative sticking probabilities S0 as a function of incident energy, for CH4 dissociation on Ni(111) at 475 K. Results are plotted for methane initially in the ground (black), 1ν 3 (red), and 2ν 3 (green) vibrational states. QCT results are shown for both the full Hamiltonian (thin solid lines) and the perturbative form discussed in the text (dashed). Quantum results are plotted as thick solid lines. The symbols are experimental data from the groups of Utz57 (A) and Beck (R).58

approach, coupled with our sudden model for treating lattice motion, is in good agreement with experimental data from the groups of Beck58 and Utz,57 particularly for the excited states. Our methods clearly reproduce the strong enhancement in S0 with vibrational excitation, as well as with surface temperature. As in earlier studies the method tends to overestimate the ground state reactivity near saturation, possibly due to the assumption of rotational adiabaticity. At higher energies, the molecules are less likely to properly reorient, and a better treatment would lead to lower values of P0 and S0 . Another factor may be our use of the PBE functional, which often gives dissociation barriers that are too low. Second, the QCT results for vibrationally excited states, at least for this model system, are in good agreement with the quantum results, and thus also the experimental data. The QCT results are a bit large for the 1ν 3 state, but overall they are reasonable. However, the QCT approach seriously overestimates the ground state S0 , due to ZPE problems. As a result, the quantum calculation gives larger efficacies than the QCT study. Guo and Jiang41 note similar behavior, comparing efficacies from their QCT studies with those from their quantum studies, though we note it was not possible for their rigid-surface fixed-site QCT probabilities to be compared directly with the quantum sticking probabilities. We see here that this lower efficacy is the result of having a ground state reactivity that is too large, due to improper IVR from the large ZPE. We note that this may be a general problem for QCT and AIMD studies involving molecules with several H atoms. A third observation is that the small differences between the full and perturbative QCT results become even more minor after impact site and lattice averaging. This is consistent with our quantum studies, where adding higher order terms in bss , for a given vibrationally adiabatic basis set, did not significantly change the reaction probabilities.

J. Chem. Phys. 139, 194701 (2013)

FIG. 6. Single-site static lattice reaction probabilities P0 as a function of incident energy, for CH4 dissociation on Ni(100). Results are plotted for methane initially in the ground (black), 1ν 3 (red), and 1ν 1 (blue) vibrational states. QCT results are shown for both the full Hamiltonian (thin solid lines) and the perturbative form discussed in the text (dashed). Quantum results are plotted as thick solid lines.

In Fig. 6, we plot P0 for CH4 dissociation on the Ni(100) surface. While the activation energy is lower on this surface relative to the (111) surface by about 0.16 eV, consistent with experimental measurements,5, 6, 57, 58 the results are otherwise similar. For the molecule initially in the ground state, the QCT results are again too large for energies below Eact . For the vibrationally excited states, the results are a bit better, as classical (over-the-barrier) processes are allowed for incident energies down to about 0.4 eV. In Fig. 7, we plot the corresponding fully averaged sticking probabilities and compare with experiment. Again, the QCT efficacies are smaller than the quantum ones, but mostly because the ground state QCT sticking

FIG. 7. Dissociative sticking probabilities S0 as a function of incident energy, for CH4 dissociation on Ni(100) at 475 K. Results are plotted for methane initially in the ground (black), 1ν 3 (red), and 1ν 1 (blue) vibrational states. QCT results are shown for both the full Hamiltonian (thin solid lines) and the perturbative form discussed in the text (dashed). Quantum results are plotted as thick solid lines. The symbols are experimental data from the groups of Utz6 (A) and Beck5 (R).


194701-7


probability is too large, by over an order of magnitude at lower energies relative to the quantum result. QCT and quantum sticking probabilities are again similar for the excited states, and the results are again consistent with available data. Note that both quantum and QCT approaches using the RPH are capable of reproducing the larger vibrational efficacy of the 1ν 1 state relative to the nearly degenerate 1ν 3 state. In the RPH formulation, reactions are vibrationally enhanced when the non-adiabatic couplings lead to transitions to other vibrationally adiabatic states with lower energy. The excess energy is converted into motion along the reaction path, which at the transition state corresponds to bond breaking. The most effective transition is thus to the vibrationally adiabatic ground state. Our earlier work32 has shown that part of the reason that the 1ν 3 state is less efficacious is that one of the three components of the 1ν 3 state, the 1 mode, cannot couple to the ground state, since B1 ,F = 0 by symmetry. The 1ν 1 efficacy is also large because mode softening lowers the effective barrier of the 3 (1ν 1 ) state. This also increases the reactivity of the 2 state, which strongly couples to the 3 . Thus, overall, the efficacies of both states are large, but the 1ν 1 efficacy is larger. As for Ni(111), the perturbative and full QCT results for P0 on Ni(100) are very similar, and nearly indistinguishable when averaged to compute S0 . Note also that for Ni(111) the perturbative results for P0 are a bit too large at lower energies and too small at higher energies. This is apparently not typical for the perturbative approach, since for Ni(100) this situation is reversed. As noted, as a Ni atom moves out of the plane of the surface, away from the bulk, the barrier for dissociation over this atom decreases dramatically. In an earlier study,32 we showed that at high substrate temperatures, sticking is dominated by reactions at lattice sites where the Ni atom is sufficiently puckered to allow for an over-the-barrier dissociation. Thus, tunneling does not make a significant contribution to S0 at these high temperatures, at least according to our models. Put another way, at these high temperatures the sticking probabilities in Figs. 5 and 7 are relatively insensitive to values of the static surface probabilities, P0 , below a few percent. Thus, the inability of the QCT approach to include tunneling is not an issue here. Of course, the QCT tends to overestimate P0 at these energies in any case, for reasons discussed. Unfortunately, the unphysical flow of ZPE within the molecule leads to errors in the QCT calculations at much larger values of P0 . Another way to measure the extent of incorrect IVR is to examine the vibrational and translational energy distributions of molecules backscattered into the asymptotic region, where H is separable. In Fig. 8, we plot the probability distributions for both total center-of-mass translational energy and total vibrational energy, for unreacted CH4 molecules scattered from Ni(111), using our QCT approach. The incident translational energy, 0.9 eV, is just below the activation energy, 0.945 eV, and the methane is initially in the ground vibrational state. The initial distributions are represented by the vertical lines. An exact quantum calculation would give distributions for the scattered molecules similar to the initial distribution, with additional small delta-function peaks corresponding to vibrational excitation, which is weak for this system. We see that both distributions are broadened, and that

J. Chem. Phys. 139, 194701 (2013)

FIG. 8. Probability distributions for the center-of-mass translational energy and total vibrational energy, for CH4 scattered from Ni(111) at an incident energy of 0.9 eV. The molecule is initially in the ground vibrational state. QCT results are shown for the initial state and for the scattered unreacted molecules.

there is a net transfer of energy from translation to vibration, on average, of about 0.03 eV. This is expected, as there is significantly more energy in translation than in any one vibrational degree of freedom. The broadening in the translational energy distribution is roughly consistent with the observation that 30% of the trajectories cross the barrier (unphysically) at this energy, keeping in mind that the V-to-T transfer from the ZPE might be even larger for the trajectories that went on to react. However, the half-widths at half-max of these distributions are only about 10% of the total energy, suggesting that the QCT are reasonably well behaved. We show that this is not the case in Fig. 9, where we plot individual probability distributions for the 9 normal modes, for the same conditions as Fig. 8. The vertical lines again correspond to the initial energy distributions (the 3 initial distribution covers that for

FIG. 9. Probability distributions for vibrational energy in the 1 , 3 , 6 , 1 , and 3 normal modes, for CH4 scattered from Ni(111) at an incident energy of 0.9 eV. The molecule is initially in the ground vibrational state. QCT results are shown for the initial state and for the scattered unreacted molecules.


194701-8


the 6 , as these modes are degenerate). In the absence of incorrect IVR, the final distributions would be the same, with additional delta function peaks for any vibrational excitations. There is significant broadening of the 1 state. The 2 distribution (not shown) is very similar to that of the 3 mode, and these are even broader. Distributions for the 4 and 5 modes (not shown) are similar to the 6 mode, and all three peak near zero energy, indicating a significant amount of energy transfer into other DOFs. The relative amount of broadening and the heights of the peaks near zero energy scale roughly with the magnitude of the curvature coupling, Bq,F , in Fig. 2. For the A states this coupling is zero, and any broadening and bad IVR comes from the generally smaller coriolis term, π s . As a result, the distributions for the A modes are less broad. In spite of this broadening and the shifts in the peaks towards zero energy, the average energy in each mode changes by only a few meV from the initial values. The exception is the 3 symmetric stretch, which increases by 25 meV. As shown in Fig. 1, this is the mode most strongly perturbed by the interaction of the molecule with the surface. We observe similar behavior in the scattered distributions for CH4 on Ni(100).

IV. CONCLUSIONS

We have used DFT-based electronic structure methods to construct a reaction path Hamiltonian for CH4 dissociation on the Ni(100) and Ni(111) surfaces. Both quantum and quasiclassical trajectory approaches are used to compute reaction probabilities for a 13-DOF system, corresponding to collision with a rigid surface at the top site, where the barrier is lowest. Impact site averaging and the effects of lattice motion are then introduced using sudden models, to get the dissociative sticking probability. All parameters used in our models are computed from DFT. We have improved our quantum model by expanding the vibrationally adiabatic basis set to include all two-quanta excitations, and have correspondingly expanded our Hamiltonian to include terms that are second order in the vibrationally non-adiabatic coupling bss . New quantum results for methane dissociation on Ni(100) and (111) are presented, and while we can now treat molecules initially in twoquanta excited states, results for single-quantum states are not significantly different from our old results. We have used QCT to show that our perturbative expansion of the RPH in the curvature coupling term, bss , is very accurate at all energies. We find that both the quantum and classical approaches, when applied to the RPH, show a large enhancement in dissociative sticking when the incident molecule is vibrationally excited, as has been observed in numerous experiments. In addition, both methods can reproduce the observed mode specificity. In particular, for CH4 dissociation on Ni(100), both the quantum and classical methods give a much larger sticking probability for molecules excited to the 1ν 1 state (symmetric stretch) than for those excited to the nearly degenerate 1ν 3 state (antisymmetric stretch), consistent with experimental observation.5, 6 However, the magnitude of this enhancement in sticking with vibrational excitation, when computed using QCT, is much smaller than that computed using our quantum approach or observed in the experiments.

J. Chem. Phys. 139, 194701 (2013)

When compared with experiment, the quantum-RPH method does a very good job at reproducing the dissociative sticking probabilities for molecules in vibrationally excited states. However, for molecules initially in the ground state, while the quantum sticking probabilities are reasonable at low incident energies they are almost an order of magnitude larger than experimental values near saturation. This may be due to our assumption of rotational adiabaticity, as molecules would have less time to reorient at higher collision energies, lowering S0 . Another factor may be our use of the PBE functional, which often gives dissociation barriers that are too low. On the other hand, given this Hamiltonian, QCT methods significantly overestimate the ground state dissociative sticking at all energies. At the highest energies, the QCT results are only a bit larger than the quantum, but at lower energies the QCT sticking can be over two orders of magnitude larger than in the experiments. We find that the origin of this behavior is an unphysical flow of zero point energy from the nine normal vibrational modes into the reaction coordinate, giving large values for P0 at energies below the activation energy. For any particular trajectory, there can be a significant redistribution of the vibrational ZPE, or any additional excitation energy, among the various molecular DOFs. Thus, in addition to significantly overestimating the ground state reactivity, the relative enhancement in S0 with vibrational excitation is too small when computed via QCT. Thus, some care must be taken when using molecular dynamics or AIMD approaches to study these reactions. ACKNOWLEDGMENTS

B. Jackson gratefully acknowledges support from the Division of Chemical Sciences, Office of Basic Energy Sciences, Office of Energy Research, (U.S.) Department of Energy (DOE), under Grant No. DE-FG02-87ER13744. 1 J.

H. Larsen and I. Chorkendorff, Surf. Sci. Rep. 35, 163 (1999). B. F. Juurlink, D. R. Killelea, and A. L. Utz, Prog. Surf. Sci. 84, 69 (2009). 3 A. L. Utz, Curr. Opin. Solid State Mater. Sci. 13, 4 (2009). 4 S. Nave, A. K. Tiwari, and B. Jackson, J. Chem. Phys. 132, 054705 (2010). 5 P. Maroni, D. C. Papageorgopoulos, M. Sacchi, T. T. Dang, R. D. Beck, and T. R. Rizzo, Phys. Rev. Lett. 94, 246104 (2005). 6 L. B. F. Juurlink, P. R. McCabe, R. R. Smith, C. L. DiCologero, and A. L. Utz, Phys. Rev. Lett. 83, 868 (1999). 7 L. B. F. Juurlink, R. R. Smith, D. R. Killelea, and A. L. Utz, Phys. Rev. Lett. 94, 208303 (2005). 8 D. R. Killelea, V. L. Campbell, N. S. Shuman, and A. L. Utz, Science 319, 790 (2008). 9 A. C. Luntz and D. S. Bethune, J. Chem. Phys. 90, 1274 (1989). 10 J. Harris, J. Simon, A. C. Luntz, C. B. Mullins, and C. T. Rettner, Phys. Rev. Lett. 67, 652 (1991). 11 A. C. Luntz and J. Harris, Surf. Sci. 258, 397 (1991). 12 D. J. Oakes, M. R. S. Mccoustra, and M. A. Chesters, Faraday Discuss. 96, 325 (1993). 13 P. M. Holmblad, J. Wambach, and I. Chorkendorff, J. Chem. Phys. 102, 8255 (1995). 14 R. C. Egeberg, S. Ullmann, I. Alstrup, C. B. Mullins, and I. Chorkendorff, Surf. Sci. 497, 183 (2002). 15 D. R. Killelea, V. L. Campbell, N. S. Shuman, R. R. Smith, and A. L. Utz, J. Phys. Chem. C 113, 20618 (2009). 16 R. Bisson, M. Sacchi, and R. D. Beck, J. Chem. Phys. 132, 094702 (2010). 17 A. C. Luntz, J. Chem. Phys. 102, 8264 (1995). 18 M. N. Carré and B. Jackson, J. Chem. Phys. 108, 3722 (1998). 2 L.


194701-9


J. Chem. Phys. 139, 194701 (2013)

19 Y. Xiang, J. Z. H. Zhang, and D. Y. Wang, J. Chem. Phys. 117, 7698 (2002).

39 D.

20 Y.

40 B.

Xiang and J. Z. H. Zhang, J. Chem. Phys. 118, 8954 (2003). 21 S. Nave and B. Jackson, Phys. Rev. Lett. 98, 173003 (2007). 22 S. Nave and B. Jackson, J. Chem. Phys. 127, 224702 (2007). 23 S. Nave and B. Jackson, J. Chem. Phys. 130, 054701 (2009). 24 A. K. Tiwari, S. Nave, and B. Jackson, Phys. Rev. Lett. 103, 253201 (2009). 25 S. Nave and B. Jackson, Phys. Rev. B 81, 233408 (2010). 26 A. K. Tiwari, S. Nave, and B. Jackson, J. Chem. Phys. 132, 134702 (2010). 27 M. M. Teixidor and F. Huarte-Larranaga, Chem. Phys. 399, 264 (2012). 28 A. P. J. Jansen and H. Burghgraef, Surf. Sci. 344, 149 (1995). 29 R. Milot and A. P. J. Jansen, Phys. Rev. B 61, 15657 (2000). 30 G. P. Krishnamohan, R. A. Olsen, G. J. Kroes, F. Gatti, and S. Woittequand, J. Chem. Phys. 133, 144308 (2010). 31 K. G. Prasanna, R. A. Olsen, A. Valdes, and G. J. Kroes, Phys. Chem. Chem. Phys. 12, 7654 (2010). 32 B. Jackson and S. Nave, J. Chem. Phys. 135, 114701 (2011). 33 M. Sacchi, D. J. Wales, and S. J. Jenkins, J. Phys. Chem. C 115, 21832 (2011). 34 M. Sacchi, D. J. Wales, and S. J. Jenkins, Comput. Theor. Chem. 990, 144 (2012). 35 M. Sacchi, D. J. Wales, and S. J. Jenkins, Phys. Chem. Chem. Phys. 14, 15879 (2012). 36 R. A. Marcus, J. Chem. Phys. 45, 4493 (1966). 37 W. H. Miller, N. C. Handy, and J. E. Adams, J. Chem. Phys. 72, 99 (1980). 38 B. Jackson and S. Nave, J. Chem. Phys. 138, 174705 (2013).

Han, S. Nave, and B. Jackson, J. Phys. Chem. A 117, 8651 (2013). Jiang, R. Liu, J. Li, D. Q. Xie, M. H. Yang, and H. Guo, Chem. Sci. 4, 3249 (2013). 41 B. Jiang and H. Guo, J. Phys. Chem. C 117, 16127 (2013). 42 B. Jiang and H. Guo, J. Chem. Phys. 138, 234104 (2013). 43 G. Henkelman, B. P. Uberuaga, and H. Jónsson, J. Chem. Phys. 113, 9901 (2000). 44 G. Henkelman and H. Jónsson, J. Chem. Phys. 113, 9978 (2000). 45 G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993). 46 G. Kresse and J. Hafner, Phys. Rev. B 49, 14251 (1994). 47 G. Kresse and J. Furthmuller, Phys. Rev. B 54, 11169 (1996). 48 G. Kresse and J. Furthmuller, Comput. Mater. Sci. 6, 15 (1996). 49 G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999). 50 P. E. Blöchl, Phys. Rev. B 50, 17953 (1994). 51 J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996). 52 J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 78, 1396 (1997). 53 M. B. Lee, Q. Y. Yang, and S. T. Ceyer, J. Chem. Phys. 87, 2724 (1987). 54 L. B. F. Juurlink, R. R. Smith, and A. L. Utz, Faraday Discuss. 117, 147 (2000). 55 J. Q. Dai and J. Z. H. Zhang, J. Phys. Chem. 100, 6898 (1996). 56 D. H. Zhang, Q. Wu, and J. Z. H. Zhang, J. Chem. Phys. 102, 124 (1995). 57 R. R. Smith, D. R. Killelea, D. F. DelSesto, and A. L. Utz, Science 304, 992 (2004). 58 R. Bisson, M. Sacchi, T. T. Dang, B. Yoder, P. Maroni, and R. D. Beck, J. Phys. Chem. A 111, 12679 (2007).


First-principles quantum dynamical theory for the dissociative chemisorption of H2O on rigid Cu(111).

Methane dissociative chemisorption and detailed balance on Pt(111): dynamical constraints and the modest influence of tunneling.

Dissociative chemisorption of methane on Ni and Pt surfaces: mode-specific chemistry and the effects of lattice motion.

The dissociative chemisorption of water on Ni(111): Mode- and bond-selective chemistry on metal surfaces.

Quantum modelling of hydrogen chemisorption on graphene and graphite.

classical limit based on selective use of the quantum potential.

Validity of the site-averaging approximation for modeling the dissociative chemisorption of H2 on Cu(111) surface: a quantum dynamics study on two potential energy surfaces.

Dynamics of water dissociative chemisorption on Ni(111): effects of impact sites and incident angles.

The classical and quantum dynamics of molecular spins on graphene.

Rotational and steric effects in water dissociative chemisorption on Ni(111).

Six-dimensional quantum dynamics for dissociative chemisorption of H2 and D2 on Ag(111) on a permutation invariant potential energy surface.

Another Look at the Mechanisms of Hydride Transfer Enzymes with Quantum and Classical Transition Path Sampling.

The mechanism of chemisorption of hydrogen atom on graphene: insights from the reaction force and reaction electronic flux.

A transfer hamiltonian model for devices based on quantum dot arrays.

Dissociative electron attachment studies on acetone.

Molecular machines operating on the nanoscale: from classical to quantum.

reaction system for in situ and time-resolved DRIFTS studies of catalytic reactions on solid surfaces.

Realization of a quantum Hamiltonian Boolean logic gate on the Si(001):H surface.

The quantum-optics Hamiltonian in the Multipolar gauge.

Understanding the reaction mechanism and intermediate stabilization in mammalian serine racemase using multiscale quantum-classical simulations.

Effects of the interplay between initial state and Hamiltonian on the thermalization of isolated quantum many-body systems.

The effect of classical and quantum dynamics on vibrational frequency shifts of H2 in clathrate hydrates.

Six-dimensional quantum dynamics study for the dissociative adsorption of HCl on Au(111) surface.

Six-dimensional quantum dynamics study for the dissociative adsorption of DCl on Au(111) surface.