First-principles surface reaction rates by ring polymer molecular dynamics and neural network potential: role of anharmonicity and lattice motion

Abstract

Elementary gas–surface processes are essential steps in heterogeneous catalysis. A predictive understanding of catalytic mechanisms remains challenging due largely to difficulties in accurately characterizing the kinetics of such steps. Experimentally, thermal rates for elementary surface reactions can now be measured using a novel velocity imaging technique, providing a stringent testing ground for ab initio rate theories. Here, we propose to combine ring polymer molecular dynamics (RPMD) rate theory with state-of-the-art first-principles-determined neural network potential to calculate surface reaction rates. Taking NO desorption from Pd(111) as an example, we show that the harmonic approximation and the neglect of lattice motion in the commonly-used transition state theory overestimates and underestimates the entropy change during the desorption process, respectively, leading to opposite errors in rate coefficient predictions and artificial error cancellations. Including anharmonicity and lattice motion, our results reveal a generally neglected surface entropy change due to significant local structural change during desorption and obtain the right answer for the right reasons. Although quantum effects are found to be less important in this system, the proposed approach establishes a more reliable theoretical benchmark for accurately predicting the kinetics of elementary gas–surface processes.

---

1. Introduction

Heterogeneous catalysis plays an irreplaceable role in large-scale chemical production and energy conversion. But achieving a predictive understanding of heterogeneous catalysis remains a grand challenge, as each catalytic reaction may involve a chain of interfacial elementary processes.¹ Accurate kinetics of elementary reactions at catalyst surfaces are thus crucial for reliable microkinetic modelling of whole catalytic processes across a wide range of reaction conditions that can guide rational catalysis design.² Experimental methods for measuring surface reaction rate coefficients were limited until very recently. Wodtke, Kitsopoulos and coworkers have accurately measured the thermal rates of a number of representative elementary reactions at metal surfaces by detecting velocity-resolved kinetic traces via ion imaging techniques.^3–11 Such experiments not only shed valuable light onto site-dependent surface reaction mechanisms,⁵ inspiring further ab initio molecular dynamics (AIMD) investigations,^12,13 but also reveal quantum effects in the thermal recombination of H₂ at metal surfaces.¹⁰ These experimental data are of great value for benchmarking surface reaction rate theories.⁶

On the other hand, theoretical predictions of surface reaction rate coefficients, especially in heterogeneous catalysis, rely largely on traditional transition state theory (TST)^14,15 using energetic information of stationary points that can be routinely optimized by density functional theory (DFT). However, the traditional TST method for surface reactions¹⁶ often invokes several over-simplified assumptions. Firstly, it considers only the transition state (TS) with the lowest energy and neglects multiple site-dependence. Secondly, it uses a harmonic approximation to estimate partition functions. Finally, it assumes that the partition function of the surface site is unchanged by adsorption/desorption and neglects coupling between molecular and surface degrees of freedom (DOFs), not to mention the intrinsic absence of recrossing and quantum tunneling in the TST. With such simplifications, TST-predicted surface reaction rate coefficients sometimes deviate heavily from experimental data.^5,10 Improvements have been made to the TST for surface reactions by using, for example, the hindered translator and hindered rotor models^17–19 or complete potential energy sampling.²⁰ Lattice vibrational modes (i.e. phonons) and dynamical effects (i.e. recrossing) were, however, rarely included in these improved TST models. One exception was the recent work by Jackson²¹ who calculated the TST rate constant for CH₄ dissociation on Ir(111), adding 81 lattice DOFs to the partition functions.

There are several theoretical studies on surface reaction kinetics which try to go beyond this TST framework and include lattice motion. Gillan and Alfè derived the formulae for calculating the thermal desorption rate of the adsorbate(s) with varying coverage on a surface based on chemical potentials,²² which can be compared with a temperature programmed desorption (TPD) experiment. They applied this approach to the desorption of H₂O/D₂O on MgO(001) using an empirical potential^23,24 or AIMD at the DFT level.^22,25 Wang and Zhao calculated the diffusion rate coefficients of H on Ni(100) and the dissociation/recombination rate coefficients of CH₄ on Ni(111) using the quantum instanton approach based on embedded diatomics-in-molecules force fields.^26–28 Recently, Chen et al. computed thermal diffusion rates of CO on Ru(0001) with AIMD using a novel enhanced sampling scheme.²⁹ Galparsoro et al. directly calculated the reactive flux for thermal recombinative desorption of hydrogen from Cu(111) via classical AIMD and compared it with a thermal desorption experiment based on a high-temperature permeation method.³⁰ These studies, however, suffer from either the low accuracy of the empirical potential or the high cost of the AIMD simulation.

A recently proposed ring polymer molecular dynamics (RPMD) approach offers a promising way to calculate rate coefficients with high accuracy.^31–33 Exploiting the isomorphism between the statistical properties of a quantum system and those of a classical system in which each particle is represented by the centroid of a ring polymer,³⁴ the path integral-based RPMD method has proven to yield correct quantum mechanical rate coefficients in the harmonic potential, high-temperature, and short-time limits.^31,35 The RPMD rate theory has been widely applied in calculating rate coefficients of bimolecular reactions in the gas phase based on highly-accurate ab initio potential energy surfaces (PESs) and has been shown to be able to properly describe quantum effects, e.g., zero-point energy and tunneling, representing a superior alternative to TST-based methods.^36–50 However, RPMD has rarely been applied in surface reactions, where many more DOFs are involved. In an early study, Suleimanov applied RPMD rate theory to calculate the diffusion constants of hydrogen on Ni(100) with an embedded atom method interatomic potential.⁵¹ More recently, non-equilibrium RPMD simulations for gas–surface scattering was carried out to calculate dissociative sticking probabilities^52–55 and diffusion constants of hydrogen on a single-atom catalytic surface.⁵⁶

On the other hand, recently developed machine learning methods have greatly facilitated the construction of high-dimensional gas–surface PESs^57–61 with first-principles accuracy, which offer a new opportunity for the efficient prediction of surface reaction rates. In this work, we combine RPMD rate theory with a state-of-the-art neural network (NN) based PES including all relevant DOFs. This combination allows us to compute thermal rate coefficients of elementary gas–surface processes efficiently and accurately. We chose to study the desorption of NO on Pd(111) as the first example, for which accurate velocity-resolved kinetic data between 620 K and 800 K have been measured very recently.⁸ Our work differs from recent work by Chen et al., who calculated the diffusion constants of CO on Ru(0001) with a neural network potential on a rigid surface.⁶² As we will show below, surface relaxation can significantly affect the entropy change and thus the desorption rate coefficients.

The rest of this article is organized as follows. Section 2 presents the potential energy surface of the system of interest and an overview of the general RPMD rate theory adapted to gas–surface reactions. Section 3 discusses how the harmonic approximation and the neglect of surface DOFs affect the desorption rate coefficients of NO from Pd(111) where entropy plays a critical role. We provide conclusions in Section 4.

2. Computational details

2.1. Potential energy surface

A global PES including all relevant molecular and surface DOFs is necessary for computing rate coefficients using MD-based approaches. To achieve high efficiency, an analytical PES for the NO + Pd(111) system was constructed in this work by the recently developed embedded atom neural network (EANN) method, which has been detailed elsewhere.^63–65 In the EANN approach, the total energy is expressed as the sum of atomic energies and the complexity of the PES thus scales linearly with system size. Each atomic energy is dependent on its atomic environment, which is represented by an array of embedded atom density (EAD) descriptors.⁶⁶ Such EAD descriptors are obtained from the square of the linear combination of Gaussian-type orbitals (GTOs) located at neighbor atoms inside a cutoff sphere with a radius of 6 Å, serving as the input vector of each atomic NN to output the atomic energy. The GTO is expressed as,

φ(r_ij) = x^l_xy^l_yz^l_z exp(−α|r − r_s|²), (1)

where r_ij and r are the Cartesian coordinates of the embedded atom i relative to atom j and its norm respectively; the parameters α and r_s determine its radial distribution; while orbital angular momenta (l_x, l_y, l_z) specify its spatial distribution. In practice, α and r_s were set as α = 0.2 Å⁻² and r_s = 0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0 Å, and the total orbital angular momentum (L = l_x + l_y + l_z) was set to 0, 1, or 2. This setup led finally to 33 EAD descriptors in the input layer of each atomic NN consisting of two hidden layers with 40 and 60 neurons.

The training dataset containing 2503 DFT points with both energies and forces was generated by a trajectory-free active learning strategy.^67,68 All DFT calculations were performed by the Vienna ab initio simulation package (VASP)^69,70 with spin-polarization using the revised Perdew–Burke–Ernzerhof (RPBE) functional.⁷¹ The flat surface of Pd(111) was modeled by a four-layer slab, which has a 3 × 3 surface unit cell with the top two layers relaxed and a vacuum region of 15 Å in the Z direction. The kinetic energy in the plane wave basis set was truncated at 400 eV and the Brillouin zone was sampled by a 5 × 5 × 1 Monkhorst–Pack k-point mesh. The adsorption energy was found to converge within a few dozens of meV with respect to the number of layers and k-points.

2.2. Ring polymer molecular dynamics

The RPMD rate theory has been well established elsewhere,⁷² so here we give only a brief summary. The quantum mechanical expression of the thermal rate coefficient was first formulated in terms of the Kubo-transformed flux–side correlation function,^73–75where Q_r(T) is the reactant partition function and c_fs(t) the Kubo-transformed flux–side correlation function. In the RPMD framework, Craig and Manolopoulos derived a path-integral version of the Kubo-transformed flux–side correlation function,^31,33 namely,Here H_n(p,q) is the ring-polymer Hamiltonian of the N-atom system,where each atom is replaced by a ring polymer consisting of n beads that are interconnected via harmonic springs with a spring frequency of ω_n = 1/(β_nℏ); β_n = 1/nk_BT is the n-bead-reduced reciprocal temperature; p_i^(j) and q_i^(j) are the momentum and position vectors of the j^th bead of the i^th atom; m_i is the atomic mass of the i^th atom; and V is the original PES of the N-atom system. In particular, s̄(q) in eqn (3) is the dividing surface defined in terms of reaction coordinate ξ(q̄), namely,

s̄(q) = ξ(q̄₁, …, q̄_N) − ξ^‡ (5)

and v̄_s(p,q) is its conjugated speed,where ξ^‡ is the reaction coordinate at the position of the dividing surface; p̄_i and q̄_i are the centroid momentum and position vectors of the i^th atom. h[s̄(q_t)] in eqn (3) is a Heaviside function of the dividing surface, which counts the fraction of trajectories passing through the dividing surface to the product side at time t. It is worthwhile pointing out that the RPMD rate coefficient is independent of the choice of the dividing surface and naturally reduces to the classical MD result when using a single bead.³³

To describe molecular desorption from the surface, we defined the reaction coordinate ξ(q̄) as the difference between the height of the molecular center of mass (COM) and the average height of the top layer of the surface as

where M_ad (N_ad) and M_s (N_s) are the mass (number of atoms) of the adsorbate and the top layer surface; z̄_i and z̄_i′ are the z coordinates of the i^th atomic centroid of the adsorbate and the i′^th atomic centroid of the top layer surface. Since the NO desorption is barrierless, the TS is assumed to be the free molecule plus the clean surface in the asymptote; the position of the dividing surface ξ^‡ is thus placed in the asymptotic area (here, ξ^‡ = 13.4 bohr above the uppermost surface layer).

In practical calculations, it is advantageous to rewrite the RPMD rate coefficient in the Bennett–Chandler form,^76,77

k(T) = κ(t → ∞)k_QTST(T). (8)

Here, κ(t → ∞) = c_fs(t → ∞)/c_fs(t → 0⁺) is the transmission coefficient in the long-time limit accounting for the recrossing effect at the dividing surface, which can be calculated by sampling trajectories constrained at the dividing surface (ξ^‡) with constant temperature and propagating them to a sufficient long time with constant energy to estimate the ratio of the two flux–side correlation functions,

On the other hand, k_QTST(T) = c_fs(t → 0⁺)/Q_r(T) is the centroid-density quantum transition state theory (QTST) rate coefficient in the short-time limit.³³ The QTST rate coefficient can be conveniently evaluated by the centroid potential of mean force (PMF) along the reaction coordinate,

where the PMF, W(ξ), can be obtained via some standard enhanced sampling approaches like umbrella sampling⁷⁸ and umbrella integration.⁷⁹ All RPMD calculations have been performed with a modified version of the RPMDrate program.⁴¹

For comparison, we also computed rate coefficients with a harmonic TST (HTST) model,

where the transmission coefficient κ is simply assumed to be one; E₀ is the zero-point energy (ZPE)-corrected binding energy; is the partition function of the free NO molecule in a (1 × 1) unit cell excluding its translation in the z direction (i.e. the reaction coordinate); Q_s is the clean surface component and Q_ad is the adsorbate–surface system. Note that in common TST implementations, the partition function of the surface component is generally assumed to be unchanged during desorption, so that Q_s will be cancelled out and Q_ad will depend only on the DOFs of the adsorbate.

In both RPMD and HTST calculations, since the electronic DOFs were not explicitly taken into account, as suggested in the experimental work,⁸ we introduced approximate electronic partition functions, namely Q^el_g = 2 + 2 exp(−ΔE_1/2→3/2/k_BT) where ΔE_1/2→3/2 = 120 cm⁻¹ is the energy separation between the two lowest spin–orbital states of NO⁸⁰ and Q^el_ad = 2, for the free and adsorbed NO molecules, respectively.

3. Results and discussion

3.1. Potential energy surface

Consistent with previous theoretical studies,^81–83 the NO molecule is found to adsorb perpendicularly on the fcc site with the N-terminal facing down with an adsorption energy of 1.844 eV on the relaxed Pd(111) surface. After the ZPE correction, the calculated binding energy E₀ is 1.773 eV, very close to the experimental value (1.766 ± 0.024 eV) derived by fitting the measured data to a hindered-translator TST model.⁸ The NO adsorbate can diffuse on Pd(111) from the fcc site to the hcp site (∼0.045 eV higher in energy) through the bridge site and the corresponding barrier is 0.320 eV, in good agreement with the experimental estimate (0.29 ± 0.11 eV).⁸ It should be noted that the experimentally fitted diffusion barrier was not super accurate because the C_3v symmetry of the Pd(111) surface was not held well by simply using two orthogonal one-dimensional potentials.⁸ On the other hand, the dissociation barrier of NO on Pd(111) was reported as 2.7 eV using the RPBE functional.⁸¹ Experimentally, NO dissociation was also never detected at surface temperatures from 620 K to 800 K.⁸ This implies that NO dissociation is much more difficult than its desorption on Pd(111) and would not affect the computation of desorption rates. Consequently, the NO + Pd(111) EANN PES in this work covers only the region where the N–O distance (r) is less than 1.8 Å. A two-dimensional cut of the PES is illustrated in Fig. 1, demonstrating the energy barrierless feature of molecular desorption and the smoothness of the PES that is beneficial to long-time MD simulations.

Fig. 1 Two-dimensional contour plot of the EANN PES as a function of Z (the distance from the COM of NO to the Pd(111) surface) and r at the fcc site of a frozen Pd(111) surface with other coordinates optimized.

The prediction root mean squared errors (RMSEs) of the PES with respect to the total energy and atomic forces are 37.7 meV (per cell; that is, ∼1.9 meV per mobile atom) and 37.3 meV Å⁻¹, respectively. The stationary points of the PES compare well with the DFT data, as shown in Table 1, where the geometries and vibrational harmonic frequencies of NO are listed. The good agreement indicates that the EANN PES faithfully reproduces the DFT description along the reaction coordinate for NO desorption from Pd(111).

Table 1 Comparison of properties of stationary points on the EANN PES and calculated by DFT^a

Species	Property	EANN PES	DFT
a θ refers to the angle between the vector pointing from N to O and the surface normal.
Isolated NO	r (Å)	1.174	1.176
	ν (cm^−1)	1848.0	1897.2
Adsorbed NO	Z (Å)	1.923	1.931
	θ (°)	0	0
	r (Å)	1.217	1.216
	ν 1 (cm^−1)	1553.8	1547.5
	ν 2 (cm^−1)	432.4	433.6
	ν 3 (cm^−1)	432.3	433.1
	ν 4 (cm^−1)	320.6	324.6
	ν 5 (cm^−1)	182.0	184.7
	ν 6 (cm^−1)	182.0	183.0

---

3.2. Convergence tests on PMF

In the RPMD rate theory, the rate coefficient is obtained by calculating the transmission coefficient κ(t → ∞) and QTST rate coefficient k_QTST(T), separately. Specifically, k_QTST(T) depends on the PMF along the reaction coordinate. Here, the PMF is calculated by umbrella sampling⁷⁸ and umbrella integration.⁷⁹ In umbrella sampling, one defines a number of windows along the reaction coordinate ξ(q̄) and a biased ring-polymer Hamiltonian, H̃_n = (p,q), is constructed at each window,where k_l is the strength of the umbrella potential at the l^th window and ξ_l is the reaction coordinate at the position of the l^th window. Note that here, these windows are evenly separated by a fixed interval (Δξ_l) with an identical strength of the umbrella potential. After sampling the biased probability density distributions of ξ(q̄) in each window, the PMF profile is calculated by umbrella integration along the reaction coordinate.

As the first RPMD application for a surface desorption process, we performed extensive convergence tests in calculating the PMF profile. For simplicity, we first tested in the frozen surface calculations with one bead (RPMD-1, or classical MD) the influence of some parameters, including the strength of the umbrella potential (k_l), the window interval (Δξ_l) as well as the single trajectory time. In the umbrella integration, the probability density of ξ(q̄) in each window follows a normal distribution centered at ξ_l. It is advisable to have a sufficient overlap between the probability density distributions of neighboring windows. The ratio between the crossing point of two adjacent normal distributions and their average height was defined as an overlap factor (λ), which correlates k_l with Δξ_l, namely,

and was used in our tests in practice. As shown in Fig. 2, the PMF profile is relatively insensitive to these parameters over a reasonable range. Specifically, 64 windows (or equivalently Δξ = 0.2 bohr) and 100 consecutive trajectories each with a sampling time of 125 ps can well converge the PMF profile to less than 1 kJ mol⁻¹. However, including surface DOFs significantly enlarges the configuration space, thus requiring a much longer sampling time to cover the space. Consequently, the number of consecutive trajectories increases to 400 at 600 K and 800 K and even 600 at 700 K so as to converge the PMF profile on a relaxed Pd(111) surface with the top two mobile layers, reaching a maximum sampling time of 4.8 μs. For example, Fig. 2d shows that the PMF profile at 800 K converges well with λ = 0.75, Δξ_l = 0.2 bohr, and 400 trajectories with 125 ps of each single trajectory.

Fig. 2 Convergence tests of RPMD-1 (MD) calculations. (a) PMF profiles for λ = 0.25, 0.5 and 0.75 obtained by evolving 30 trajectories for 50 ps with Δξ_l equal to 0.2 bohr. (b) PMF profiles at Δξ_l of 0.2, 0.4 and 0.8 bohr sampled over 30 trajectories evolved for 50 ps with λ equal to 0.75. (c) PMF profiles sampled by evolving 30 trajectories for 50, 75, 100, 125, 150 ps with λ and Δξ_l set to 0.75 and 0.2 bohr. (d) PMF profiles at different numbers of trajectories for the relaxed surface at 800 K with λ = 0.75, Δξ_l = 0.2 bohr, and 125 ps for a single trajectory.

The converged PMF results with one bead on the frozen and relaxed surfaces at different temperatures are compared in Fig. 3. As expected, the increasing temperature gradually lowers the free energy of desorption, promoting the desorption of NO on Pd(111). Interestingly, the surface motion modifies the shape of the PMF profile, leaving the profile less harmonic along the reaction coordinate from the adsorbate (minimum) to the asymptote (TS). This is likely to be due to the involvement of surface DOFs allowing for anharmonic coupling between the molecular and surface vibrations, which may soften the frequency along the reaction coordinate. However, the overall free energy change is less affected by the lattice motion, the reasons for which will be discussed below. Moreover, to check the effect of NO diffusion on PMF convergence, snapshots of two representative trajectories at 700 K in the window near the PMF well are taken every 0.5 ps and are shown in Fig. 4. Clearly, NO can jump to neighboring sites and diffuse across unit cells. The diffusion occurs faster on the frozen surface than on the relaxed surface, because the kinetic energy of NO cannot dissipate to surface phonons in the former case. This explains why convergence of the PMF profiles for the relaxed surface is more difficult than for the frozen surface. In practice, the sufficiently long-time umbrella sampling by hundreds of consecutive trajectories adequately describes the diffusional behavior and guarantees the molecular traversal on the surface in determining the desorption PMF profiles.

Fig. 3 Converged PMF profiles for the frozen (solid lines) and relaxed (dotted lines) surfaces at 600 K (black), 700 K (blue) and 800 K (red), obtained by RPMD-1 (MD) calculations where one bead is used.

Fig. 4 Exemplary RPMD-1 (MD) trajectories showing NO diffusion on the relaxed (red) and frozen (blue) Pd(111) surfaces during the umbrella sampling in the adsorption well for 125 ps after thermostatting of 20 ps at 700 K. The top Pd atoms of the relaxed and frozen surfaces are represented by cyan and gray circles, respectively. Both trajectories are initiated from the fcc site.

3.3. Quantum effects

The RPMD rate theory is known to approximately capture quantum effects with multiple beads. While this has been routinely validated in polyatomic reactions in the gas phase,⁸⁴ multi-bead RPMD simulations in gas–surface systems that involve many more DOFs remain very time-consuming. Fortunately, since all atoms in this NO + Pd(111) system are heavy, it is unlikely that the quantum tunneling effect will play any significant role. The ZPE change from the adsorbed NO to the free NO molecule is the most probable source of quantum effects. To check this possibility, we have calculated PMFs for NO desorption on a frozen Pd(111) at 700 K with different numbers of beads and the results are compared in Fig. 5. It is clear that the PMF profile barely changes with the number of beads increasing from one to eight, indicating a negligible quantum effect in this system. Unlike hydrogen atom diffusion^51,85 (at low temperatures) and recombination¹⁰ on the metal surface where remarkable nuclear quantum effects were identified, our results suggest that quantum effects are less important in the case of molecular desorption of a heavy adsorbate, at least at high temperatures. Indeed, the diffusion of a hydrogen atom was also found to be not greatly influenced by quantum effects at moderate temperatures (e.g., 300 K).⁵⁶ This is understandable since the molecular desorption (rather than recombination) process involves neither the breaking of a strong chemical bond nor a big change in ZPE. In addition, the binding energy is so large that the desorption has to occur at a very high temperature. It is thus encouraging that a classical description (one bead) is sufficient to converge the PMF of this process, which is adopted for all subsequent calculations. Hereafter, the corresponding single-bead results with and without consideration of the transmission coefficient are referred to as MD and classical TST (cTST), respectively.

Fig. 5 PMF profiles calculated by RPMD calculations with 1 bead (black solid line), 4 beads (red dashed line), and 8 beads (blue dotted line) at 700 K where Pd(111) is frozen.

3.4. The effect of harmonic approximation

One of the major advantages of the RPMD method over the traditional TST is that the PMF is accurately determined from MD trajectories rather than invoking the harmonic approximation at static stationary points. To analyze the influence of anharmonicity alone, the surface is temporarily fixed at its equilibrium and the transmission coefficient is temporarily neglected. In this case, the MD rate coefficients reduce to the cTST ones, which are compared with those obtained by the HTST based on DFT-optimized energies and harmonic frequencies of stationary points on the frozen surface at various temperatures in Fig. 6a. Not surprisingly, the cTST rate coefficients are about 3–4 times lower than the HTST ones. To clearly analyze this difference, the changes in thermodynamic variables in the desorption process are evaluated. In the cTST case, the internal energy along the reaction coordinate U(ξ) was calculated from centroid virial theorem⁸⁶ with the system constrained at the reaction coordinate ξ via the umbrella potential,and the associated entropy S(ξ) was given by TS(ξ) = U(ξ) − W(ξ). In the HTST case, the thermodynamic variables were computed from their relationships with partition functions. The difference in the rate coefficients can be ascribed primarily to the difference in entropy changes from the adsorbate to the free molecule (at the dividing surface) calculated by the two methods, as shown in Fig. 6b. It is clear that the internal energy changes (ΔU) are quite close from the two methods, as they are mainly determined by the binding energy of NO. However, the harmonic approximation seems to more strongly underestimate the entropy of the adsorbate (with more low-frequency modes) than the free molecule (with a single high-frequency vibrational mode), leading to higher TΔS values in the HTST model, and thus lower free energy barriers than those in cTST. Our results suggest the importance of anharmonicity in all adsorbate modes.

Fig. 6 (a) Comparison between cTST rate coefficients (blue solid line with circles) and HTST ones (red dotted line) on the frozen Pd(111). (b) Changes in the free energy (black), internal energy (red) and the entropy (blue) estimated by umbrella sampling in cTST (solid lines) calculations and by analytical partition functions (dashed lines) where the adsorbed NO on the frozen Pd(111) is described with six harmonic oscillators and the x and y translational partition functions of the asymptotic NO molecule are calculated in a (1 × 1) unit cell.

3.5. The lattice effect

Another important advantage of the current RPMD method is that the couplings of DOFs belonging to the molecule and the substrate are fully included. In comparison, the traditional TST method is typically implemented by including vibrational frequencies associated with the adsorbate motion alone, while most of the time neglecting partition functions of surface DOFs, although the optimized stationary structures are obtained with surface relaxation. If the interactions between the adsorbate and the substrate in the reactant and barrier regions are very different, lattice motion could play a significant role. To check this possibility, the calculated cTST rate coefficients with frozen and relaxed surfaces are compared in Fig. 7a. At first glance, it is surprising that the desorption rate coefficients in the two cases are almost identical. This phenomenon can be rationalized after investigating the changes in internal energy, entropy and free energy separately in the desorption process, as shown in Fig. 7b. It turns out that the internal energy change on the non-frozen surface is higher than that on the frozen surface by about 71 meV, due largely to the increased binding energy via surface relaxation. In contrast, the entropy change on the non-frozen surface also becomes higher than that on the frozen surface by about 103 meV. The increase in entropy change in the former case can be understood as follows. Specifically, the molecule–surface interaction in the favorable adsorption site is so strong that the local surface configuration has been largely rearranged and the three Pd atoms closest to the NO molecule deviate most from their equilibrium positions, as shown in Fig. 8. This interaction confines the surface structure and leads naturally to a lower entropy than that for the unconstrained equilibrium surface structure after the molecule desorbs. This part of entropy change is absent in the frozen surface approximation. A similar entropy change of molecules adsorbed at different sites has actually been discussed before in the desorption of Pb from Mo(100) by Starr and Campbell, who found a large difference in the entropy of the adsorbates between the strongly bound step site and the weakly bound terrace site.⁸⁷ What is different is that the change in surface structure is essential for the entropy change here.

Fig. 7 (a) Comparison between the cTST rate coefficients for the frozen (blue solid line with circles) and the relaxed Pd(111) surfaces (red dashed line with triangles). (b) Changes in the free energy (black), internal energy (red) and entropy (blue) estimated by umbrella sampling for the frozen (solid lines) and relaxed (dashed lines) surfaces in cTST calculations.

Fig. 8 Rearrangement of local Pd(111) surface structures upon NO adsorption where the largest displacements (in Å) among Pd atoms before and after adsorption are marked and the displacement directions are indicated by black arrows. The O, N and Pd atoms with and without the adsorbate are represented by red, blue, cyan and gray spheres, respectively.

Overall, the effects of the changes in internal energy and entropy due to the lattice motion cancel out, leaving similar free energies and rate coefficients with or without surface DOFs. Similarly, in the RPMD calculations of H diffusion on Ni(100) with either one bead (classical) or multiple beads (quantum), a very minor lattice effect was also found on the diffusion coefficient, especially above the quantum-classical crossover temperature (70 K).⁵¹ For the dissociation of CH₄ on Ir(111), Jackson also found that the lattice vibrational modes only lead to a slight change to the HTST rate constant.²¹ This is possibly because the entropy change in the dissociation of CH₄ is mainly contributed by the breaking of the C–H bond rather than a change in the lattice structure.

3.6. The transmission coefficient

While the transmission coefficient (κ) is often neglected in the traditional TST, the RPMD rate theory accounts for this factor in a quantitative way, by the ratio of the flux–side correlation function in the long-time limit to that at time zero. In Fig. 9a, the time-dependent transmission coefficients for the frozen and relaxed surfaces at 700 K are compared. Interestingly, the transmission coefficient exhibits completely different behaviors in the two cases. For the frozen surface, κ experiences a plateau at the very beginning (the first ∼1 ps) and then decreases quite rapidly. Its value decays to ∼0.1 at 100 ps and appears to further decay to zero when the time is long enough. This behavior is uncommon in gas-phase reactions, where κ typically converges within tens to hundreds of fs.^45,46 For the relaxed surface, however, the transmission coefficient converges to 0.92 very quickly at ∼4 ps. This indicates that the recrossing is practically insignificant in this process, which is consistent with the experimental observation.⁸ The large difference in κ in the two cases may have to be understood in terms of energy dissipation. Fig. 10 illustrates the evolution of the kinetic energies of the NO molecule moving towards the surface and of Pd atoms in the top two layers, and the associated reaction coordinate (ξ) as a function of time during two representative trajectories for computing κ on the frozen and relaxed surfaces. In any case, NO molecules slide from the TS to the adsorption well and are accelerated by the attraction. When the surface is frozen, the kinetic energy of NO has no way to dissipate and only exchanges with the potential energy. Thanks to conservation of energy, the total energy of NO on the frozen surface is always larger than the adsorption energy. As a result, the highly energetic NO adsorbate will scatter back after sufficient time and recross the dividing surface with a velocity towards the vacuum (Fig. 10a), which makes a negative contribution to eqn (9). In this way, κ decays to zero in the long-time limit. By contrast, on the relaxed surface, the molecular kinetic energy can largely transfer to surface Pd atoms, which stabilizes the adsorbate and prevents recrossing (Fig. 10b).

Fig. 9 (a) Time-dependent MD transmission coefficients for the desorption of NO from the relaxed (red) and frozen (blue) Pd(111) surfaces at 700 K. (b) Time-dependent MD transmission coefficients for the relaxed Pd(111) surface at 600 (black), 700 (red) and 800 (green) K.

Fig. 10 The reaction coordinates and kinetic energies of the NO molecule and the Pd atoms in the upper two layers as a function of time for two exemplary trajectories for the frozen (a) and relaxed (b) surfaces at 700 K.

In addition, Fig. 9b shows that κ on the relaxed surface decreases slightly with increasing temperature, possibly because the recrossing is enhanced at a higher temperature as the adsorbate has a higher thermal energy. These results underscore the importance of considering energy dissipation in the adsorption/desorption rate calculation. Therefore, the transmission coefficients are included in the following comparison of MD rate coefficients with experimental ones. Our observation is somewhat different from that in the RPMD calculations for the site-to-site hopping of hydrogen on the Ni(100) surface, where the recrossing was found not to be affected significantly by the lattice motion.⁵¹ It is likely that the diffusion of the adsorbed hydrogen is weakly coupled with the surface motion.

3.7. Comparison between the HTST and experiment

In Fig. 11, the MD rate coefficients of NO desorption on the movable Pd(111) surface are compared with the measured velocity-resolved kinetic data.⁸ The MD results give a good prediction of the shape of the log k ∼ 1/T curve and the fitted activation energy (E_a = 1.731 eV) is very close to the experimental one (E_a = 1.71 ± 0.03 eV), as evidenced by the good agreement with the experimental binding energy E₀ shown above. Whereas absolute rate coefficients are 3–4 times larger than experimental ones, because the calculated pre-factor (A = 10^16.37 s⁻¹) relying on the multi-dimensional topography of the PES deviates from experimental data (A = 10^15.65±0.20 s⁻¹). This discrepancy may have multiple sources of errors. For example, the diffusion barrier is probably too high so that the adsorbate entropy is underestimated, leading to desorption rate coefficients that are too large. In addition, the 3 × 3 supercell may not be large enough to eliminate the adsorbate–adsorbate interactions and/or perfectly reflect the surface structure change closest to the adsorbate so that the surface entropy change is underestimated. Another possibility is the neglect of the interaction of the adsorbate and metallic electron–hole pairs which is expected to be much stronger for the adsorbed NO than for the asymptotic NO, as observed for the NO + Au(111) case.⁸⁸ This nonadiabatic effect may slow down the molecular motion so as to lower the desorption rate coefficients. In fact, Gu et al. have found that H* diffusion is slowed by electronic friction due to interaction with surface electron–hole pairs, leading to better agreement with experiment.⁵⁶ More investigations in this direction will be interesting.

Fig. 11 Comparison between two sets of HTST rate coefficients (blue solid line and green dashed line), MD ones (red solid line with triangles) for the relaxed surface including transmission coefficients and the experimental data⁸ (black ×). The HTST (Relaxed-6D) (blue solid line) model only considers partition functions of the NO molecule based on the optimized DFT parameters on the relaxed Pd(111), while the HTST (Relaxed-FD) (green dashed line) model additionally includes the surface partition functions via the harmonic approximation and DFT-determined frequencies of surface phonons.

Also compared are two sets of HTST results assuming κ = 1. It is important to note that the effect of surface entropy was not well realized in common TST applications. Instead, the effect of surface relaxation on the binding energy was taken into account by including surface DOFs when optimizing stationary structures, but only the partition functions associated with the adsorbate DOFs were usually considered. According to our analysis above, the absence of entropy change contributed by lattice motion will overestimate the free energy change. On the other hand, the harmonic approximation will underestimate the free energy change. The combination of both approximations will thus result in an artificial error cancellation in the commonly used HTST model, rendering a spurious coincidence with the MD results, as displayed as HTST (Relaxed-6D) in Fig. 11. In comparison, the full-dimensional HTST model (HTST (Relaxed-FD)) considering all surface vibrational modes yields much higher rate coefficients than the MD and experimental results. Note that experimental data were used to fit a hindered-translator TST model in ref. 8, assuming a square symmetry of the Pd(111) surface to determine the diffusion barrier and parallel translational frequencies, while all parameters here are taken from DFT calculations. Our results to some extent justify the use of the traditional TST model with the assumption of neglecting the surface entropy, at least in the molecular desorption process. More investigations for other surface reactions will be interesting.

4. Conclusions

In this work, we combine RPMD rate theory with a newly-constructed machine-learned PES from first principles to calculate desorption rate coefficients of NO from Pd(111) including all relevant DOFs. This treatment moves beyond the traditional TST for surface processes, including anharmonicity, adsorbate–substrate coupling, dynamical recrossing, and quantum effects. In particular, we find that the harmonic approximation underestimates the entropy of the adsorbate, resulting in an overestimation of desorption rate coefficients. When the surface is relaxed, the local surface structural change increases the surface entropy change during the desorption process, which offsets the increase in binding energy and keeps the free energy change close to that in the frozen surface approximation. Our results indicate that the commonly-used harmonic TST approach relying on stationary points and energies with surface relaxation but only considering partition functions of the molecule DOFs actually suffers from error cancellation. Moreover, we find that the transmission coefficient is physically sound only when the energy dissipation between the adsorbate and substrate is considered with a movable surface. Overall, while the binding energy is well captured, the theoretical desorption rate coefficients remain 3–4 times higher than the experimental ones, suggesting the need for further improvement in the current PES or the adiabatic dynamics based model. The proposed computational scheme will be useful for including quantum effects that are expected to be more important in surface reactions with light atoms^10,51 from first principles instead of using parametric models and be applicable to more general elementary surface reaction rates towards accurate microkinetic modeling of heterogeneous catalysis.

Data availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Author contributions

B. J. conceived the idea behind this work and led the supervision of this research. C. L. did all calculations with the initial supervision of Y. L. All authors contributed to the analysis and discussion of the calculated data. C. L. and B. J. wrote the first draft of the manuscript and all authors edited its subsequent versions.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by Strategic Priority Research Program of the Chinese Academy of Sciences (XDB0450101 to B. J.), CAS Project for Young Scientists in Basic Research (YSBR-005 to B. J.), National Natural Science Foundation of China (No. 22221003, 22073089, and 22033007 to B. J., No. 22173057 to Y. L.), Science and Technology Commission of Shanghai Municipality (No. 21JC1402700 to Y. L.). Calculations were done on the Supercomputing Center of USTC and Hefei Advanced Computing Center. We thank Dr Dmitriy Borodin and Prof. Alec M. Wodtke for providing us the experimental data and thoughtful discussion, and Prof. Hua Guo for helpful suggestions.

References

1. A. Hellman, E. J. Baerends, M. Biczysko, T. Bligaard, C. H. Christensen, D. C. Clary, S. Dahl, R. van Harrevelt, K. Honkala, H. Jonsson, G. J. Kroes, M. Luppi, U. Manthe, J. K. Nørskov, R. A. Olsen, J. Rossmeisl, E. Skúlason, C. S. Tautermann, A. J. C. Varandas and J. K. Vincent, J. Phys. Chem. B, 2006, 110, 17719–17735.
2. Z. Chen, H. Wang, N. Q. Su, S. Duan, T. Shen and X. Xu, ACS Catal., 2018, 8, 5816–5826.
3. K. Golibrzuch, P. R. Shirhatti, J. Geweke, J. Werdecker, A. Kandratsenka, D. J. Auerbach, A. M. Wodtke and C. Bartels, J. Am. Chem. Soc., 2015, 137, 1465–1475.
4. D. J. Harding, J. Neugebohren, H. Hahn, D. J. Auerbach, T. N. Kitsopoulos and A. M. Wodtke, J. Chem. Phys., 2017, 147, 013939.
5. J. Neugebohren, D. Borodin, H. W. Hahn, J. Altschäffel, A. Kandratsenka, D. J. Auerbach, C. T. Campbell, D. Schwarzer, D. J. Harding, A. M. Wodtke and T. N. Kitsopoulos, Nature, 2018, 558, 280–283.
6. G. B. Park, T. N. Kitsopoulos, D. Borodin, K. Golibrzuch, J. Neugebohren, D. J. Auerbach, C. T. Campbell and A. M. Wodtke, Nat. Rev. Chem., 2019, 3, 723–732.
7. D. Borodin, K. Golibrzuch, M. Schwarzer, J. Fingerhut, G. Skoulatakis, D. Schwarzer, T. Seelemann, T. Kitsopoulos and A. M. Wodtke, ACS Catal., 2020, 10, 14056–14066.
8. D. Borodin, I. Rahinov, J. Fingerhut, M. Schwarzer, S. Hörandl, G. Skoulatakis, D. Schwarzer, T. N. Kitsopoulos and A. M. Wodtke, J. Phys. Chem. C, 2021, 125, 11773–11781.
9. D. Borodin, I. Rahinov, O. Galparsoro, J. Fingerhut, M. Schwarzer, K. Golibrzuch, G. Skoulatakis, D. J. Auerbach, A. Kandratsenka, D. Schwarzer, T. N. Kitsopoulos and A. M. Wodtke, J. Am. Chem. Soc., 2021, 143, 18305–18316.
10. D. Borodin, N. Hertl, G. B. Park, M. Schwarzer, J. Fingerhut, Y. Wang, J. Zuo, F. Nitz, G. Skoulatakis, A. Kandratsenka, D. J. Auerbach, D. Schwarzer, H. Guo, T. N. Kitsopoulos and A. M. Wodtke, Science, 2022, 377, 394–398.
11. M. Schwarzer, N. Hertl, F. Nitz, D. Borodin, J. Fingerhut, T. N. Kitsopoulos and A. M. Wodtke, J. Phys. Chem. C, 2022, 126, 14500–14508.
12. L. Zhou, A. Kandratsenka, C. T. Campbell, A. M. Wodtke and H. Guo, Angew. Chem., Int. Ed., 2019, 58, 6916–6920.
13. M. del Cueto, X. Zhou, L. Zhou, Y. Zhang, B. Jiang and H. Guo, J. Phys. Chem. C, 2020, 124, 5174–5181.
14. H. Eyring, Chem. Rev., 1935, 17, 65–77.
15. M. G. Evans and M. Polanyi, Trans. Faraday Soc., 1935, 31, 875–894.
16. P. L. Houston, Chemical Kinetics and Reaction Dynamics, New York, 2001, vol. 78, pp. 1466.
17. C. T. Campbell, L. H. Sprowl and L. Árnadóttir, J. Phys. Chem. C, 2016, 120, 10283–10297.
18. L. H. Sprowl, C. T. Campbell and L. Árnadóttir, J. Phys. Chem. C, 2016, 120, 9719–9731.
19. R. B. McClurg, R. C. Flagan and W. A. Goddard III, J. Chem. Phys., 1997, 106, 6675–6680.
20. M. Jørgensen and H. Grönbeck, J. Phys. Chem. C, 2017, 121, 7199–7207.
21. B. Jackson, J. Chem. Phys., 2020, 153, 034704.
22. D. Alfè and M. J. Gillan, J. Phys.: Condens. Matter, 2006, 18, L451–L457.
23. H. Fox, M. J. Gillan and A. P. Horsfield, Surf. Sci., 2007, 601, 5016–5025.
24. H. Fox, M. J. Gillan and A. P. Horsfield, Surf. Sci., 2009, 603, 2171–2178.
25. D. Alfè and M. J. Gillan, J. Chem. Phys., 2007, 127, 114709.
26. W. Wang and Y. Zhao, J. Phys. Chem. A, 2015, 119, 12953–12961.
27. W. Wang and Y. Zhao, J. Chem. Phys., 2017, 147, 044703.
28. W. Wang and Y. Zhao, J. Chem. Phys., 2009, 130, 114708.
29. Z.-N. Chen, L. Shen, M. Yang, G. Fu and H. Hu, J. Phys. Chem. C, 2015, 119, 26422–26428.
30. O. Galparsoro, S. Kaufmann, D. J. Auerbach, A. Kandratsenka and A. M. Wodtke, Phys. Chem. Chem. Phys., 2020, 22, 17532–17539.
31. I. R. Craig and D. E. Manolopoulos, J. Chem. Phys., 2004, 121, 3368–3373.
32. I. R. Craig and D. E. Manolopoulos, J. Chem. Phys., 2005, 122, 084106.
33. I. R. Craig and D. E. Manolopoulos, J. Chem. Phys., 2005, 123, 034102.
34. D. Chandler and P. G. Wolynes, J. Chem. Phys., 1981, 74, 4078–4095.
35. S. Habershon, D. E. Manolopoulos, T. E. Markland and T. F. Miller III, Annu. Rev. Phys. Chem., 2013, 64, 387–413.
36. R. Collepardo-Guevara, Y. V. Suleimanov and D. E. Manolopoulos, J. Chem. Phys., 2009, 130, 174713.
37. Y. V. Suleimanov, R. Collepardo-Guevara and D. E. Manolopoulos, J. Chem. Phys., 2011, 134, 044131.
38. R. Pérez de Tudela, F. J. Aoiz, Y. V. Suleimanov and D. E. Manolopoulos, J. Phys. Chem. Lett., 2012, 3, 493–497.
39. Y. V. Suleimanov, R. Pérez de Tudela, P. G. Jambrina, J. F. Castillo, V. Sáez-Rábanos, D. E. Manolopoulos and F. J. Aoiz, Phys. Chem. Chem. Phys., 2013, 15, 3655–3665.
40. J. W. Allen, W. H. Green, Y. Li, H. Guo and Y. V. Suleimanov, J. Chem. Phys., 2013, 138, 221103.
41. Y. V. Suleimanov, J. W. Allen and W. H. Green, Comput. Phys. Commun., 2013, 184, 833–840.
42. R. Pérez de Tudela, Y. V. Suleimanov, J. O. Richardson, V. Sáez Rábanos, W. H. Green and F. J. Aoiz, J. Phys. Chem. Lett., 2014, 5, 4219–4224.
43. J. Espinosa-Garcia, M. Garcia-Chamorro, J. C. Corchado, S. Bhowmick and Y. V. Suleimanov, Phys. Chem. Chem. Phys., 2020, 22, 13790–13801.
44. Y. Li, Y. V. Suleimanov, J. Li, W. H. Green and H. Guo, J. Chem. Phys., 2013, 138, 094307.
45. Y. Li, Y. V. Suleimanov, M. Yang, W. H. Green and H. Guo, J. Phys. Chem. Lett., 2013, 4, 48–52.
46. Y. Li, Y. V. Suleimanov, W. H. Green and H. Guo, J. Phys. Chem. A, 2014, 118, 1989–1996.
47. Y. Li, Y. V. Suleimanov and H. Guo, J. Phys. Chem. Lett., 2014, 5, 700–705.
48. J. Zuo, Y. Li, H. Guo and D. Xie, J. Phys. Chem. A, 2016, 120, 3433–3440.
49. J. Zuo, Y. Li, H. Guo and D. Xie, J. Phys. Chem. A, 2017, 121, 5067.
50. H. Wang, J. Fang, H. Yang, J. Song and Y. Li, Chem. Phys. Lett., 2019, 730, 227–233.
51. Y. V. Suleimanov, J. Phys. Chem. C, 2012, 116, 11141–11153.
52. H. Jiang, M. Kammler, F. Ding, Y. Dorenkamp, F. R. Manby, A. M. Wodtke, T. F. Miller, A. Kandratsenka and O. Bünermann, Science, 2019, 364, 379.
53. H. Jiang, X. Tao, M. Kammler, F. Ding, A. M. Wodtke, A. Kandratsenka, T. F. Miller and O. Bünermann, J. Phys. Chem. Lett., 2021, 12, 1991–1996.
54. Q. Liu, L. Zhang, Y. Li and B. Jiang, J. Phys. Chem. Lett., 2019, 10, 7475–7481.
55. C. Li, Q. Liu, L. Zhang, Y. Li and B. Jiang, Mol. Phys., 2021, 119, e1941367.
56. K. Gu, C. Li, B. Jiang, S. Lin and H. Guo, J. Phys. Chem. C, 2022, 126, 17093–17101.
57. B. Jiang, J. Li and H. Guo, J. Phys. Chem. Lett., 2020, 11, 5120–5131.
58. K. Shakouri, J. Behler, J. Meyer and G.-J. Kroes, J. Phys. Chem. Lett., 2017, 8, 2131–2136.
59. R. Moiraghi, A. Lozano, E. Peterson, A. Utz, W. Dong and H. F. Busnengo, J. Phys. Chem. Lett., 2020, 11, 2211–2218.
60. B. Kolb, X. Luo, X. Zhou, B. Jiang and H. Guo, J. Phys. Chem. Lett., 2017, 8, 666–672.
61. A. Rivero Santamaría, M. Ramos, M. Alducin, H. F. Busnengo, R. Díez Muiño and J. I. Juaristi, J. Phys. Chem. A, 2021, 125, 2588–2600.
62. J. Chen, T. Jin, Y. Jiang, T. Shen, M. Yang and Z.-N. Chen, Chin. Chem. Lett., 2022, 33, 4936–4942.
63. Y. Zhang, C. Hu and B. Jiang, J. Phys. Chem. Lett., 2019, 10, 4962–4967.
64. Y. Zhang, C. Hu and B. Jiang, Phys. Chem. Chem. Phys., 2021, 23, 1815–1821.
65. Y. Zhang, J. Xia and B. Jiang, Phys. Rev. Lett., 2021, 127, 156002.
66. J. Xia, Y. Zhang and B. Jiang, Chin. J. Chem. Phys., 2021, 34, 695–703.
67. Q. Lin, Y. Zhang, B. Zhao and B. Jiang, J. Chem. Phys., 2020, 152, 154104.
68. Q. Lin, L. Zhang, Y. Zhang and B. Jiang, J. Chem. Theory Comput., 2021, 17, 2691–2701.
69. G. Kresse and J. Furthmuller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
70. G. Kresse and J. Furthmuller, Comput. Mater. Sci., 1996, 6, 15–50.
71. B. Hammer, L. B. Hansen and J. K. Nørskov, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 7413–7421.
72. R. Collepardo-Guevara, I. R. Craig and D. E. Manolopoulos, J. Chem. Phys., 2008, 128, 144502.
73. T. Yamamoto, J. Chem. Phys., 1960, 33, 281.
74. W. H. Miller, J. Chem. Phys., 1974, 61, 1823–1834.
75. W. H. Miller, S. D. Schwartz and J. W. Tromp, J. Chem. Phys., 1983, 79, 4889–4899.
76. C. H. Bennett, in Algorithms for Chemical Computations, ACS Symposium Series, ed. R. E. Christofferson, ACS, 1977, vol. 46.
77. D. Chandler, J. Chem. Phys., 1978, 68, 2959–2970.
78. G. M. Torrie and J. P. Valleau, J. Comput. Phys., 1977, 23, 187–199.
79. J. Kästner and W. Thiel, J. Chem. Phys., 2005, 123, 144104.
80. J. M. Brown, A. R. H. Cole and F. R. Honey, Mol. Phys., 1972, 23, 287–295.
81. B. Hammer, J. Catal., 2001, 199, 171–176.
82. M. Gajdoš, J. Hafner and A. Eichler, J. Phys.: Condens. Matter, 2005, 18, 13–40.
83. J. A. Herron, S. Tonelli and M. Mavrikakis, Surf. Sci., 2012, 606, 1670–1679.
84. Y. V. Suleimanov, F. J. Aoiz and H. Guo, J. Phys. Chem. A, 2016, 120, 8488–8502.
85. W. Fang, J. O. Richardson, J. Chen, X.-Z. Li and A. Michaelides, Phys. Rev. Lett., 2017, 119, 126001.
86. M. F. Herman, E. J. Bruskin and B. J. Berne, J. Chem. Phys., 1982, 76, 5150–5155.
87. D. E. Starr and C. T. Campbell, J. Am. Chem. Soc., 2008, 130, 7321–7327.
88. C. L. Box, Y. Zhang, R. Yin, B. Jiang and R. J. Maurer, JACS Au, 2021, 1, 164–173.

---