The conventional theoretical approach to calculating rate coefficients of an elementary surface reaction is based on transition-state theory (TST), 1,2 often with the harmonic approximation (HA). HA-TST estimates the enthalpic and entropic contributions to the rate coefficient using energies and harmonic vibrational frequencies computed using a first principles method, such as density functional theory (DFT).
There is ample evidence indicating that HA could lead to significant errors in the calculated entropy, due to anharmonicity in soft modes, and this is compounded by the uncertainties introduced by the generalized gradient approximation (GGA) in DFT. 3 Significant efforts have recently been devoted to more accurate molecular dynamics (MD) based free-energy calculations, using either on-the-fly DFT forces or machine learned potentials. [4][5][6][7][8] Although these MD calculations allow the extensive sampling of the thermally accessible phase space and thus avoid HA, the classical nature of the sampling method neglects nuclear quantum effects, such as zero-point energy and tunneling.
Many surface processes involving hydrogen may exhibit strong nuclear quantum effects. Accurate theoretical predictions of thermal rate coefficients influenced by quantum effects are of great interest in the field of heterogeneous catalysis. The rate theory based on ring-polymer molecular dynamics (RPMD) 9,10 leverages the isomorphism between a quantum statistical property and that of a fictitious classical necklace of harmonically connected beads. 11 This Feynman path-integral based method is thus capable of accounting for nuclear quantum effects such as zero-point energy and tunneling. 12 Anharmonicity is also included because the rate coefficient is computed on a realistic (non-harmonic) potential energy surface (PES). Applications to gas-phase bimolecular reactions have demonstrated its accuracy when compared with benchmark quantum results. 13 Recently, the RPMD rate theory has been extended to surface processes, 8,[14][15][16][17] for which high-dimensional PESs can now be routinely constructed from first principles data. 18 Because a single bead RPMD is equivalent to the classical limit, it is ideally suited for understanding the origin of nuclear quantum effects.
The recombinative desorption (RD) of surface atomic hydrogens (2H* → H2) from metal surfaces is not only an ideal test system for theory due to its simplicity, but also of great importance in heterogeneous catalysis. Until very recently, however, experimentally measured rates of this elementary surface reaction exhibited large uncertainties, which made it meaningless to compare with theoretical studies. The situation has greatly improved by the latest progress in experimental rate measurements using Velocity Resolved Kinetics 19 (VRK), which yielded highly accurate rate coefficients for this process, 20 thus providing highly desirable experimental benchmarks for theoretical calculations.
Because of the light mass of hydrogen, nuclear quantum effects could play an important role in such reactions. Indeed, the H2 dissociative adsorption (DA) on Pt(111), which is the microscopic reverse of RD, is shown in our recent RPMD rate calculations on a DFT PES to possess significant nuclear quantum effects. 16 More recently, we have implemented RPMD to compute thermal rates of activated H2/D2 DA and RD on Cu(111), 17 where excellent agreement with experimental data was also obtained.
Furthermore, significant tunneling through the activation barrier was shown at high temperatures up to 400 K.
In this work, we use the RPMD method to compute the RD rate coefficients of hydrogen and its isotopes on Pt(111) using a DFT PES, as well as its modification calibrated with experimental data. The computed RPMD rate coefficients on the modified PES are in very good (within a factor of two) agreement with the measurements. Our analysis further reveals that tunnelling is relatively unimportant, whereas the dominant nuclear quantum effect is related to the reactant zero-point energy.
Here, we focus on RD of hydrogen and related isotopomers from an infinite surface:
2H * (ad)/2D * (ad) H2(g)/D2(g) or H * (ad) + D * (ad) HD(g)
where * indicates the adsorbed state. The differential rate equations can be expressed in terms of either the molecular concentration in the gas phase or the coverage on the surface:
In experiments, it is more straightforward to measure the surface coverage, which is in units of 1/cm 2 . This results in rate coefficients 𝑘𝑘(𝑇𝑇) in units of cm 2 /s. 20 The details of the general RPMD rate theory have been extensively discussed elsewhere, 21 so we focus here on the RD process involving adsorbed atomic hydrogens and deuteriums. 17 The quantum mechanical thermal rate coefficient can be written as, 22 𝑘𝑘(𝑇𝑇) =
where 𝑄𝑄 𝑟𝑟 (𝑇𝑇) is the reactant partition function per unit volume and 𝑐𝑐 𝑓𝑓𝑓𝑓 (𝑡𝑡) the Kubotransformed flux-side correlation function 9,10 defined in the N-particle phase space {𝒑𝒑, 𝒒𝒒}:
Here, the ring-polymer Hamiltonian 𝐻𝐻 𝑁𝑁 (𝒑𝒑, 𝒒𝒒) is given as,
where each atom is replaced by a necklace of n beads that are interconnected via harmonic springs with the spring frequency of 𝜔𝜔 𝑁𝑁 =
is the n-bead reciprocal temperature, 𝒑𝒑 𝑖𝑖 (𝑗𝑗) and 𝒒𝒒 𝑖𝑖 (𝑗𝑗) are the momentum and position vectors of the j th bead of the i th atom, respectively, mi is the atomic mass of the i th atom, V is the PES of the N-atom system.
The dividing surface is defined by 𝛿𝛿[𝜉𝜉(𝒒𝒒)] along the reaction coordinate 𝜉𝜉(𝒒𝒒)
and ℎ[𝜉𝜉(𝒒𝒒 𝑑𝑑 )] is a Heaviside function, which counts the fraction of trajectories (𝒑𝒑 𝑑𝑑 , 𝒒𝒒 𝑑𝑑 ) passing through the dividing surface to the product side at time t, and 𝑣𝑣 𝜉𝜉 (𝒑𝒑, 𝒒𝒒)
is the initial velocity along the reaction coordinate 𝜉𝜉(𝒒𝒒),
In our treatment of the RD process, 17 two dividing surfaces are defined based on the ring polymer centroid. The first dividing surface is placed in the reactant region where two hydrogen atoms are adsorbed at two adjacent fcc sites,
where r is the distance between two surface hydrogen atoms and 𝑟𝑟 𝑎𝑎𝑑𝑑𝑓𝑓 is the corresponding distance at the adsorbed state of the PES. This represents an approximation in which the diffusion barrier for an adsorbed atomic hydrogen is ignored and these atoms are effectively treated as a two-dimensional free gas. This approximation is reasonable for H on Pt at the zero-coverage limit, as the diffusion barrier is low (~0.04 eV). The second dividing surface is placed in the transition-state region and can be defined as,
where Z is the vertical distance between the center of mass (COM) of the incipient molecule and the surface, 𝑍𝑍 ‡ and 𝑟𝑟 ‡ are the corresponding distances at the saddle point of the PES.
In analogy to the treatment of gas-phase bimolecular reactions, 23 the reaction coordinate 𝜉𝜉(𝒒𝒒̄) that connects the two dividing surfaces is given by, 17
where
With this definition, 𝜉𝜉 → 0 as 𝑠𝑠 0 → 0 and 𝜉𝜉 → 1 as 𝑠𝑠 1 → 0.
In practical calculations, the RPMD rate coefficient is presented in the Bennett-Chandler factorized form, 25,26 𝑘𝑘 RPMD (𝑇𝑇) = 𝑘𝑘 QTST (𝑇𝑇; 𝜉𝜉 ‡ )𝜅𝜅�𝑡𝑡 → 𝑡𝑡 𝑝𝑝 ; 𝜉𝜉 ‡ �.
Here, the static contribution of the rate constant is 𝑘𝑘 QTST (𝑇𝑇; 𝜉𝜉 ‡ ), the centroid-density quantum transition-state theory (QTST) rate coefficient. 27 It can be evaluated at the peak (𝜉𝜉 ‡ ) of the potential of mean force (PMF) along the reaction coordinate. For a RD process, we have 17
where 𝜇𝜇 𝑟𝑟 is the molecular reduced mass, and 𝑊𝑊(𝜉𝜉) is the PMF calculated by using umbrella integration: 21, 28, 29
with the probability distribution
Here Nwindows is the number of umbrella windows placed along the reaction coordinate, with each window assigned a specific value 𝜉𝜉 𝑖𝑖 . Ni represents the total number of steps sampled for window i. 𝜉𝜉 ̄𝑖𝑖 and 𝜎𝜎 𝑖𝑖 2 are the mean value and the variance calculated for the ith window from the trajectory generated by the modified ring polymer Hamiltonian,
where ki is the force constant which defines the strength of the harmonic bias in window i and
On the other hand, 𝜅𝜅�𝑡𝑡 → 𝑡𝑡 𝑝𝑝 ; 𝜉𝜉 ‡ � in Eq. ( 12) is the ring polymer transmission coefficient in the long-time limit, which is a dynamical correction accounting for the recrossing of the dividing surface 𝜉𝜉(𝒒𝒒̄) = 𝜉𝜉 ‡ with tp as a ''plateau'' time. This ensures that the resulting RPMD rate coefficient k RPMD (T) is independent of the choice of the dividing surface. 23 Specifically, it is given by 13
where the angular brackets denote a canonical average and the subscripts on the brackets indicate that the averages are taken over the constrained ensemble at 𝜉𝜉(𝒒𝒒̄) = 𝜉𝜉 ‡ and the factor of 𝑓𝑓 𝜉𝜉 (𝒒𝒒̄) -1 is a metric tensor correction for the effect of the constraint. 30 The original PES employed in this work was the one recently developed by us, 16 using the permutation invariant polynomial-neural network (PIP-NN) method 31,32 based on 1457 DFT points. A contour plot for the PES at the top site of Pt(111) is shown in Fig. 1 as a function of Z and r with the angular coordinates optimized. For more details of the DFT calculations and PES fitting, the reader is referred to our recent publication. 16 Note that the PES is based on a rigid surface approximation, which has been widely used due to the large mass disparity between hydrogen and metal atoms. 14 To improve the agreement with experimental data, this PES was modified as described below.
All calculations were performed with a home-made FORTRAN code in Cartesian atomic coordinates with the ring polymer Hamiltonian in Eq. ( 6), with no restriction on the overall rotational or translational motion of the H2 moiety on the surface. For the calculation of the PMF, the reaction coordinate was divided into 141 equally spaced windows over the 𝜉𝜉 range from -0.05 to 1.35. In each sampling window, 80 trajectories were propagated for 120 ps with a quadratic bias potential with a force constant of 2.72 (T/K) eV, after an initial equilibration period of 20 ps with the Andersen thermostat 33 for thermalization. The free-energy profile was obtained by splicing the data in each window, after the removal of the harmonic bias potential, as discussed in the previous section.
Having obtained the peak (𝜉𝜉 ‡ ) of the PMF at each temperature, the transmission coefficient was computed. A long "parent" trajectory of 2 ns was propagated with the ring-polymer centroid constrained at the peak of the PMF via the SHAKE algorithm 34 after a thermalization period of 20 ps. Configurations were saved successively in every 2 ps. For each of these configurations, 150 separate unconstrained "child" trajectories were spawned with different initial momenta sampled from a Boltzmann distribution.
These trajectories were then propagated for 150 fs to make sure that the transmission coefficient reaches a plateau. The time step is selected to be 0.1 fs in all RPMD calculations. The number of beads is successfully increased until convergence. Purely classical rate coefficients were also calculated by setting the number of beads to one.
For H2 RD on Pt(111), the thermal rate coefficients were calculated on the original DFT PES at a number of temperatures ranging from 300 to 1000 K. The parameter rads in Eq. ( 8) was 2.84 Å for all of the temperatures considered and the s1 dividing surface temperatures. In Figs. 2(b) and 2(c), the classical and RPMD PMFs are compared at
two representative temperatures (500 and 953 K). Interestingly, the RPMD free-energy barrier is slightly higher than the classical counterpart, a point to which we will later return.
The time dependence of the RPMD transmission coefficient, 𝜅𝜅�𝑡𝑡 → 𝑡𝑡 𝑝𝑝 ; 𝜉𝜉 ‡ �, the dynamical factor, is shown in Fig. 3 for several temperatures. Clearly, all the transmission coefficients quickly reach plateau values after initial drops from unity. The value is approximately 0.4 for all temperatures, suggesting significant recrossing.
The converged RPMD rate coefficients on the original DFT PES for H2 RD on Pt(111) and related calculation parameters are listed in Table I. They are also compared with the available experiment 20 in Fig. 4 and Table II. The RPMD rate coefficients are in remarkably good agreement with the experimental counterparts, differing by a maximum factor of 6 over a broad temperature range. It is important to note that the theoretical rate coefficients have been divided by a factor of 2 to account for the exchange symmetry of the two indistinguishable H particles. 17 Besides, as discussed in ref. 20, only one out of four electron spin combinations of two hydrogen atoms yields a stable gas phase molecule. Therefore, we account for this electronic quantum effect by dividing the RPMD rate coefficients by an additional factor of 4.
The relatively small discrepancy between theory and experimental can presumably be attributed to the following factors. First, the rigid surface approximation in the model used here ignores the contribution of the surface lattice motion. Despite the low mass of the H atoms relative to metal atoms, the motion of the surface atoms may still influence the kinetics, particularly at high temperatures. In principle, the impact of surface motion can be considered in RPMD calculations, but it requires a highdimensional PES that includes all the surface degrees of freedom. Second, the approximate treatment of the adsorbed H atoms as 2D free translators and the neglect of the diffusion barrier might also lead to significant errors. Removal of this approximation requires a better characterization of the reactant state in the RPMD formulation. Third, the RPMD rate theory is an approximate quantum theory, which may itself also contribute to the uncertainties. It would be desirable to be validated by exact quantum dynamical calculations for this surface process, which is possible using for example a wave packet method. Finally, errors might arise from the inadequacy of the GGA functional to compute accurate potential energies.
It is worth mentioning that the activation energy derived from an Arrhenius fit based on the RPMD rate coefficients in Table I is 0.66 eV, which is 0.09 eV lower than the experimentally estimated barrier of 0.75 eV in the zero-coverage limit. 20 To correct this discrepancy, the following empirical correction potential (in eV) was added in the original PES,
where a = -0.09 eV. This scaling does not alter the PES landscape and the vibrational frequencies, thus the impact on the entropy change is small. A comparison of the on the modified PES and experimental rate coefficients is much better than that on the original DFT PES over the temperature range, suggesting the error in the original PES is largely contained in the enthalpic contribution. The slightly larger experiment-theory differences at high temperatures are likely due to the rigid surface approximation in our model.
In Fig. 4, the RPMD rate coefficients are also compared with the classical values (1 bead RPMD). Interestingly, the rate coefficients from converged RPMD simulations are lower than those from classical limit at all temperatures except 300 K, mainly due to the higher RPMD free-energy barrier shown in Fig. 2
To make a comprehensive comparison with the experiment, we also calculated the RPMD rate coefficients for HD and D2 RD on Pt(111) based on the modified PES.
When calculating D2, a factor of 2 was also introduced to account for the exchange symmetry of the two indistinguishable D particles. However, for the calculation of HD, since the two particles (H and D) are not the same, the factor is not needed. The PMFs and transmission coefficients for these two isotopologues are given in SI, along with other results. As shown in the upper panel of Fig. 5, the calculated RPMD rate coefficients for HD are the largest followed by H2 and D2, due apparently to distinguishability of the two particles. Once we have gathered all the rate data for HD, H2 and D2, we can use them to simulate the isotopic branching fractions, and subsequently compare them to the experimental results. 20 The lower panel of Fig. 5 reveals that the relative trend of the isotopic branching fractions derived from RPMD rate coefficients are also consistent with the experimental data derived from VRK measurements. Note that the isotopic branching fraction is defined as the ratio of the yield (Y) of each isotopologue to the total yield of all species, where Y is the time integral of the product (HD/H2/D2) flux. 20
The main advantage of the RPMD approach is its ability to capture quantum effects such as zero-point energy and tunneling. In our recent work, 16 we identified small but significant quantum effects in H2 DA on a Pt(111) surface. Here, we found that the reverse RD process also has a small nuclear quantum effect, as evidenced by the relatively minor difference between the RPMD and classical rate coefficients shown in Fig. 4. This is easy to understand based on the fact that the RD is a highly endothermic process with a relatively low exit channel barrier. Tunneling through such a highly asymmetric barrier is known to be ineffective. 35 The larger difference between the classical and quantum rates in our earlier work 20 is attributable to a different benchmark in defining the classical limit.
We further argue that the zero-point energy effect plays a more prominent role in the hydrogen RD on Pt(111) than tunneling, at least within the experimental temperature range. This conclusion is based on several observations. First, tunneling often manifests in an enhanced reactivity at low temperatures in the form of a curved deviation from the Arrhenius behavior. However, it can be seen from Fig. 4 that there is only a rather slight curvature in the temperature dependence of the RPMD rate coefficients near the lowest temperature investigated here. Second, the classical rate coefficients are slightly higher than RPMD counterparts at all temperatures except 300 K. Especially, when the temperature is above 400 K, the difference between the logs of the classical and RPMD rate coefficients remains nearly constant. If tunneling were to contribute to the reactivity, the quantum rate coefficients would be larger than classical counterparts, and the difference between the classical and RPMD rate coefficients should increase with decreasing temperature. 16,17 The behaviors exhibited by the RPMD results in comparison with the classical counterparts strongly suggest that the zero-point energy plays a larger role than tunneling. However, Fig. 4 indicates that at the low temperature of 300 K, the RPMD rate coefficient slightly exceeds the classical result, suggesting a growing contribution of tunneling as the temperature decreases.
In order to more explicitly show that the zero-point energy plays a major role in this process, we calculated the zero-point energies of H2 and D2 at the transition and adsorbed states, respectively, using the harmonic vibrational frequencies from the PES.
As indicated in Table S-II in SI, the zero-point energy of the H2 transition state is 0.039 eV larger than that of the adsorbed state, consistent with the higher RPMD free-energy barrier than classical barrier as shown in Fig. 2(b)-(c). Furthermore, in Table S-III of SI, we calculated the free energy at the transition and adsorbed states using the harmonic vibrational frequencies for both H2 and D2. As shown, the free-energy difference of H2 between the transition and adsorbed states is 0.011 eV higher than that of D2. On the other hand, when comparing the RPMD free-energy barrier of H2 and D2 in Fig. S-4 of SI, we observed that the RPMD free-energy barrier of H2 is approximately 0.010 eV higher than that of D2, which is also consistent with the results obtained from the harmonic vibrational frequencies. Although the HA-TST theory is not quantitatively accurate, the trend is consistent with our RPMD results and provides further evidence for the role of the zero-point energy in the hydrogen RD on Pt(111).
The less prominent role of tunneling in this process is readily understood because the minimum energy path (MEP) at the top site for Pt(111) has no potential barrier. As a result, the PMF is dominated by the MEP. Clearly, the RPMD free-energy barrier is higher than the classical counterpart in both cases.
In this work, we combine a newly adapted RPMD rate theory with a firstprinciples PES to compute the RD rate coefficients of atomic hydrogens and its isotopes from a Pt(111) surface. The RPMD rate coefficients of H2 RD exhibit very good agreement with experiment (with a factor 2), after scaling the PES to match the experimental activation energy. Additionally, we also simulate the isotopic branching fractions of H2/D2/HD measured in previous VRK experiments. The results show that the experimental isotope effect is well captured by our calculated RD rates. More importantly, we demonstrate that the nuclear quantum effect is relatively small and primary source of the quantum effects in hydrogen RD on Pt(111) stems from the reactant zero-point energy.
The success in this prototypical system suggests the great promise of the RPMD rate theory in characterizing the kinetics of surface reactions. Given the linear scaling of the RPMD approach, and the increasingly available high-dimensional PESs, we expect its wide applications to more complex surface reactions, particularly for polyatomic species with the surface motion. Such progress, coupled with experimental advances in measuring the rate coefficients for surface processes, should significantly improve our quantitative understanding of the reaction kinetics of these interfacial reactions.
The Supporting Information is available free of charge at xxx. Comparison of the PESs, additional RPMD results different isotopomers, and HA-
Table I. Summary of centroid-density QTST rate coefficients, transmission coefficients and RPMD rate coefficients for H2 RD on the original DFT PES in temperature interval between 300 and 1000 K. The centroid-density QTST and RPMD rate coefficients are given in cm 2 /s and the numbers in parentheses denote powers of 10. Note that the RPMD rate coefficients have included the symmetry number and electron spin factor. T (K) Nbead k QTST κ k RPMD 300 32 4.75(-14) 0.40 2.37(-15) 400 32 3.95(-11) 0.38 1.88(-12) 500 32 2.18(-09) 0.36 9.81(-11) 653 16 9.14(-08) 0.36 4.11(-09) 753 16 4.59(-07) 0.36 2.07(-08) 853 953 16 16 1.56(-06) 4.17(-06) 0.37 0.36 7.21(-08) 1.88(-07) Table II. RPMD rate coefficients for H2 RD on the original and modified PESs in temperature interval between 653 and 953 K. The experimental rate coefficients are obtained from Ref. 20. The rate coefficients are given in cm 2 /s and the numbers in parentheses denote powers of 10. Note that the RPMD rate coefficients have included the symmetry number and electron spin factor. T (K) k RPMD Original PES k RPMD Modified PES Expt. 653 753 853 953 4.11(-09) 2.07(-08) 7.21(-08) 1.88(-07) 8.14(-10) 4.99(-09) 2.11(-08) 6.38(-08) 6.61(-10) 5.27(-09) 1.69(-08) 3.01(-08) fixed on the top site of Pt(111) and the angular coordinates optimized, which contains two submerged saddle points. comparison of the classical (dashed blue line) and RPMD (solid red line) PMFs at two temperatures (b and c) for H2 RD on the original DFT PES. RD on the original DFT PES. RD on the original DFT PES. The RPMD results on the modified PES are also shown for comparison. modified PES. Lower panel: Isotopic branching fraction from experiments (symbols) and RPMD calculations on the modified PES (lines). Note that some symbols have been shifted by ±5 K for clarity.