Accurate Neural Network Description of Surface Phonons in Reactive Gas-Surface Dynamics: N2 + Ru(0001).

Ab initio molecular dynamics (AIMD) simulations enable the accurate description of reactive molecule-surface scattering especially if energy transfer involving surface phonons is important. However, presently, the computational expense of AIMD rules out its application to systems where reaction probabilities are smaller than about 1%. Here we show that this problem can be overcome by a high-dimensional neural network fit of the molecule-surface interaction potential, which also incorporates the dependence on phonons by taking into account all degrees of freedom of the surface explicitly. As shown for N2 + Ru(0001), which is a prototypical case for highly activated dissociative chemisorption, the method allows an accurate description of the coupling of molecular and surface atom motion and accurately accounts for vibrational properties of the employed slab model of Ru(0001). The neural network potential allows reaction probabilities as low as 10-5 to be computed, showing good agreement with experimental results.

The accurate simulation of atomic and molecular scattering from metal surfaces often requires a reliable description of surface phonon motion,[1−4] especially if the projectile is heavier than dihydrogen. Although it may often be possible to obtain reasonably accurate descriptions of the scattering dynamics with a reduced dimensionality representation of the surface vibrations,[5] for quantitative agreement, it is mandatory to use an ab initio molecular dynamics (AIMD) approach.[6,7] AIMD recently enabled a chemically accurate description of the reaction of CHD3 with Ni(111).[8] The high accuracy of AIMD stems from its ability to include the motion of the atoms in the upper surface layers by computing the interatomic forces on-the-fly through electronic structure methods, which for molecules scattering from metal surfaces usually rely on density functional theory (DFT).[9] However, AIMD is computationally expensive, presently making simulations of long-lasting or low-probability chemical processes infeasible. This poses a severe challenge: how can we enable AIMD-quality simulations of molecule–surface reaction processes associated with low probabilities at affordable computational costs?

Here we address this question for the N2 + Ru(0001) reaction, which is considered a paradigm for highly activated dissociative adsorption. Dissociation of N2 over Ru surfaces is also of practical interest as this reaction is the rate-determining step in ammonia synthesis, for which Ru is an efficient catalyst.[10,11] The N2 + Ru(0001) system has therefore received a lot of attention in a series of experimental studies.[12−18] These experiments have shown that the dissociation probability S0 of N2 is very small (for the collisional energy range of ∼1–2.75 eV, the measured values are S0 ≈ 10–6–10–3),[16−18] even though the minimum barrier height is only about 1.8 eV.[17,19,20] It has been suggested that the low reactivity is due to nonadiabatic effects, such as electron–hole pair excitation.[21] Subsequent quasi-classical trajectory (QCT) calculations carried out within the Born–Oppenheimer static surface (BOSS) approximation showed that the low reactivity could be well explained by the barrier representing an extremely narrow bottleneck to reaction.[19,20] However, discrepancies between theory and experiment remained, and it was suggested that surface phonons might play an important role.[19,20]

As simulating the low probability N2 + Ru(0001) reaction with AIMD would require running too many (∼5 × 105) expensive AIMD trajectories, instead we investigate the role of phonons by developing a DFT-quality high-dimensional (HD) neural network potential (NNP) for N2 + Ru(0001). The applicability of artificial neural networks to describing chemical reactions has already been proven for, for example, bulk materials,[22,23] organic reactions,[24,25] simulations of solid–liquid interfaces,[26] molecules interacting with static surfaces,[27−29] and, very recently, energy transfer to metal surfaces in scattering of HCl from Au(111).[30] However, in the latter molecule–surface scattering study, reaction probabilities and scattering probabilities were not yet calculated, although final rotational state distributions in scattering were presented. In this work, for the first time, we apply the HD-NNP approach to reactive molecule–surface scattering with inclusion of surface atom motion.

Here, in order to construct a HD-NNP for N2 + Ru(0001), we employ the method of Behler–Parrinello,[31] through which the energy is expressed as a sum of environmentally dependent atomic energies, and the atomic environments are described by many-body atom-centered symmetry functions.[32,33] The accuracy of the HD-NNP is assessed by comparison with the DFT reference data for N2 + Ru(0001) as well as with available experimental results regarding the phonons of Ru. After this validation step, the HD-NNP is used to generate reaction probabilities for comparison with experiments[17,18] and to quantitatively analyze the role of phonons and the surface temperature (Ts).

The HD-NNP is developed from a data set of 25 000 N2 + Ru(0001) structures, including the concomitant information on their energy and forces. To obtain the data, we carried out DFT calculations with the RPBE functional[34] as implemented in the Vienna Ab initio Simulation Package (VASP),[35−38] which was also used in earlier studies.[19,20] The DFT calculations used a 3 × 3 surface unit cell and seven atomic layers of Ru atoms, of which the bottom layer was kept stationary while the N2 molecule was placed on the top side of the slab. The calculations employed a plane wave cutoff of 550 eV and a 7 × 7 × 1 Γ-point centered mesh of k-points.

The HD-NNP was constructed using the RuNNer package.[32,39] The DFT data included 5000 configurations of two nitrogen atoms on an ideally frozen surface and 20 000 configurations on the thermally disordered Ru(0001) surface. For the latter, the thermal displacements of the atoms from their equilibrium positions were sampled randomly from a Gaussian distribution, making sure to also include data for surface temperatures (up to Ts = 900 K) higher than those used in the experiments simulated below (up to Ts = 600 K). For both the N and the Ru atoms, radial and angular symmetry functions were used to describe interactions with neighbor atoms within their chemical environment up to a cutoff radius of 5.82 Å.[31,32]

As a first test of the HD-NNP, we examine two-dimensional (2D) cuts through the potential energy surface (PES), showing it as a function of the molecular coordinates r and Z (see Figure a) for impact on the bridge site at two different orientations of the molecule (Figure b,c). The 2D cuts provided by the HD-NNP are compared with DFT data computed on a fine grid and then interpolated with splines. In this particular case, the surface atoms were kept frozen at their equilibrium positions to facilitate the comparison. The dissociation geometries studied are close to the minimum barrier geometry found in ref (40). Direct comparison shows that the HD-NNP accurately predicts the DFT data for the geometries and the static surface configuration considered. To check whether a similar accuracy is obtained also for other dissociation geometries and for distorted surface geometries, the 25 000 points on which the HD-NNP is based was divided in a test set of 5000 points (which were picked at random from the complete set) and a training set of 20 000 points. With the fit based on the training set, the total energy of each structure present in the test data set was evaluated and a root-mean-square (RMS) deviation of 38 meV (less than 1 kcal/mol) from the DFT data was obtained for the total energies. The same or higher accuracy is expected from the HD-NNP constructed with the aid of the complete set, which was used in the reactive scattering calculations reported below.

Figure 1

(a) Illustration of molecular coordinates in the center of mass system and of four dissociation geometries. (b,c) 2D elbow maps of the PES, for N2 on Ru(0001), as a function of the (r, Z) coordinates for θ = 80° and φ = 40° (b) and for θ = φ = 90° (c). The blue dashed and orange solid curves describe the HD-NNP and the spline interpolated raw DFT data, respectively. The mesh size used for the interpolation is specified as a green rectangle. The insets schematically depict the orientation of N2 with respect to the surface. The energy difference between successive contour lines is 0.4 eV.

The so-constructed HD-NNP also accurately describes vibrational properties of the Ru(0001) slab, as evidenced by the corresponding phonon band structure and density of states. Our phonon calculations are based on the so-called “direct method”[42] as implemented in an in-house version of the phonopy code,[43,44] interfacing both the VASP and LAMMPS[45] codes and using the same computational setups as described before. We have verified that 3 × 3 multiples of the surface unit cell yield converged phonon band structures and densities of states. The results for the phonon band structure are summarized in Figure a for the path given by the high-symmetry points Γ̅–K̅–M̅–Γ̅ through the surface Brillouin zone and compared with the results obtained directly from RPBE-DFT calculations and with experimental results (ref (41)). Note that for accurate extraction of the phonon dispersion relation, the RMS error in the interatomic forces, from which the required force constants are derived, should be less than 0.001 eV/Å. The acoustic phonon mode frequencies in Figure a obtained from the HD-NNP are in excellent agreement with both the RPBE-DFT and the experimental results.[41] Likewise, the frequencies of the highest optical phonon modes at point Γ̅, that is, ω(Γ̅) ≈ 31–32 meV, are consistent with the results obtained with high-resolution electron–electron-loss spectroscopy (30.4–32 meV).[46]

Figure 2

(a) Phonon dispersion of Ru(0001) along the path through the surface Brillouin zone given by the high-symmetry points Γ̅–K̅–M̅–Γ̅. The blue dots are extracted from HD-NNP, and the orange solid curves are from RPBE-DFT calculations. The green squares indicate the acoustic phonon modes obtained from experimental measurements (ref (41)). (b) Corresponding phonon density of states as a function of the frequency of phonon modes. (c) Variations of the reaction barrier height for the four representative dissociation geometries described in Figure a, with the vertical displacement of an individual Ru atom in the topmost surface layer. The plus markers indicate the reaction barrier heights calculated through DFT for the same barrier geometry as that found from HD-NNP. (d) Vertical shift of the barrier location along the molecular coordinate Z as a function of the lattice coordinate Q⊥.

The good agreement between the NN and DFT descriptions of the phonons is not restricted to the depicted Γ̅–K̅–M̅–Γ̅ path; there is also excellent agreement between the RPBE-DFT and the HD-NNP results for the phonon density of states g(ω) for the calculations where all phonon modes in the surface Brillouin zone are taken into account (Figure b). The density function g(ω) also enables us to evaluate the effective Debye temperature TD of Ru(0001). To do this, the density of states for acoustic phonon modes needs to be fit as gD(ω) ∝ ω2 within the Debye model. The resulting TD from the HD-NNP is 410 K, which agrees well with TD = 423 K as obtained directly from our RPBE-DFT calculations and also with the values of 400–410 K obtained in low-energy electron diffraction measurements.[47,48] Altogether, the good agreement noted between the HD-NNP and DFT suggests that the HD-NNP yields an accurate description of the lattice dynamics of the Ru(0001) surface.

To further validate the NN fit, we also examine the quality of the HD-NNP for the coupling between molecular and surface atom motion, at some stationary points. Starting from the ideal rigid Ru(0001) surface and using the HD-NNP, barrier geometries and heights are computed for the four (X, Y, φ) dissociation geometries indicated in Figure a, through three-dimensional Newton–Raphson optimizations of (r, Z, θ) geometries. In Figure c, the coupling of the reaction barrier height (Eb) to the motion of an individual Ru atom in the uppermost surface layer is investigated separately for the aforementioned dissociation geometries. Vertical displacements (Q⊥) of the surface Ru atom nearest to N2 are considered. An approximation to Eb is also calculated directly through DFT by computing the DFT energy for the molecular coordinates found from the HD-NNP and shown by “+” markers for reference. As seen, the barrier height exhibits roughly a linear dependence on Q⊥, which agrees with the findings of similar studies on other adsorbate–surface systems.[49,50] In a similar vein, the vertical shift of the barrier location Zb along the molecular coordinate Z also depends approximately linearly on Q⊥, as illustrated in Figure d (in the procedure to obtain Zb through DFT calculations, the r and θ coordinates were kept at the values obtained from the HD-NNP). In summary, and following the definitions of ref (49), the electronic coupling β = −ΔEb/Q⊥ and the mechanical coupling α = ΔZb/Q⊥, which can be taken as representative of the coupling between molecular and surface atom motion, are well described with the HD-NNP. This suggests that the HD-NNP may be used with confidence to assess the effects of surface temperature and of energy exchange with phonons on the dissociative chemisorption of N2 on Ru(0001).

As a final test of the NN fit, we also checked whether the HD-NNP accurately describes the minimum barrier geometry and barrier height. For the HD-NNP, the minimum barrier geometry and height were found by performing a Newton–Raphson search incorporating all molecular coordinates but keeping the surface atoms fixed in their ideal geometry. The search was started from a bridge-to-hollow adsorption geometry as it was previously shown[20,40] that this dissociation geometry is close to the minimum barrier geometry. For comparison, we also calculated the minimum barrier coordinates and energy directly with RPBE-DFT using the climbing-image nudged elastic band[51] (CI NEB) technique. The results are summarized in Table . The barrier geometry and height derived from the HD-NNP are in excellent agreement with the raw DFT results.

Table 1

Minimum Barrier Geometry and Its Height Eb As Obtained from the HD-NNP and RPBE-DFT Calculationsa

coordinates/height of barrier	HD-NNP	RPBE-DFT
r (Å)	1.738	1.741
Z (Å)	1.312	1.318
X	0.504	0.506
Y	0.504	0.506
θ	84.3°	84.0°
φ	30.0°	30.0°
Eb (eV)	1.83	1.84

For the RPBE-DFT calculations, the CI NEB method has been employed.

We now turn to the effect of phonons on the reaction dynamics (Figure a,b). Comparisons were made with three types of experiments,[17,18] which we refer to here as Exps. 1–3. In Exp. 1, a nozzle temperature (Tn) of 1100 K was used, and seeding with H2 was carried out. The average collision energy was varied by changing the seeding ratio. Exp. 2 concerns a single experiment performed for Tn = 1150 K using He as a seeding gas. Both experiments were done by the same group and used Ts = 575 K for the surface (ref (18)). Careful error analysis was performed, estimating the extracted reaction probabilities to be accurate to within a factor 5. In Exp. 3 (ref (17)), Ts = 600 K and Tn was varied between 350 and 1000 K to change the collision energy using seeding in H2. The results were not accompanied by error analysis. We do not consider the experiments performed by the same group and published in ref (16) as it is not clear how this group could have achieved collision energies as high as 4 eV with a hydrogen-seeded beam using Tn up to 700 K (assuming an average translational energy of 2.7kBTn for H2[52] and no velocity slip, the average collision energy could not have been higher than 2.3 eV in these experiments). When comparing our dynamics results to the experiments considered, we give the highest weight to Exps. 1 and 2 as they come with thorough error analysis.

Figure 3

Dissociative sticking probability of N2 on Ru(0001) as a function of its average translational energy. (a) Quasi-classical dynamics calculations are compared with two types of experiments[18] for N2 beams seeded by H2 (Exp. 1) and He (Exp. 2). The full squares indicate the sticking probability of N2 for Ts = 575 K including the effect of phonons (HD-NNP, Ts = 575). The probability S0 is also shown for Ts = 0 with the surface treated as mobile (HD-NNP, Ts = 0). The green triangles show the results obtained within the ideal rigid surface approximation (6D dynamics). (b) Same as that in (a), but the sticking probability is compared with experimental results of ref (17) (Exp. 3) for Ts = 600 K.

In performing the dynamics calculations using the mobile surface, the seven-layer slab of Ru(0001) was first equilibrated to the surface temperature used in the experiments. The initial quantum states of N2 were statistically populated in accordance with the value of Tn used in the experiments, and our calculations mimic the velocity distribution of the molecular beams in the experiments. In our quasi-classical dynamics calculations, dissociative chemisorption was operationally defined to occur once r > 2.5req, where req = 1.117 Å denotes our RPBE value for the equilibrium N–N distance in the gas phase. In this study, 5 × 105 trajectories for incident energies of Ek < 2 eV and 105 trajectories for Ek ≥ 2 eV were run, ensuring statistically converged results over the incidence energy range of 1.5–3.25 eV. The outcome of our quasi-classical simulations is summarized in Figure . Taking into account the uncertainties in the reaction probabilities measured in Exp. 1, the calculated sticking probabilities S0 computed with the RPBE functional are in good agreement with the experimental values (cf. the full squares and circles in Figure a), and the same holds true for Exp. 2 (cf. the blank circle and star marker). There are sizable discrepancies between the calculated results and the results of Exp. 3, especially for collision energies of >2.0 eV. However, as noted before, we attach greater weight to the comparison with the results of Exps. 1 and 2. Future calculations could assess the possible effects that electron–hole pair excitation might have on the calculated reaction probability. Most likely, their inclusion in the dynamical model should somewhat diminish the calculated reaction probabilities, although work in which the effect of their inclusion was investigated for similar systems suggests a much smaller role for electron–hole pair excitation than for the phonons for a similar N2 + metal system.[3] Also, the error analysis of ref (18) clearly suggests that it would be worthwhile to invest efforts in improving the accuracy of the measurements for this paradigmatic reactive molecule–surface system.

The HD-NNP also allows us to assess the influence of surface temperature Ts on the reactivity by comparing with Ts = 0, that is, the Ru surface is static at the beginning of each trajectory but additionally accounts for energy dissipation into phonons during the collision of N2 with the surface, as well as with 6D dynamics performed within the rigid surface approximation. As to the latter, we note that for the highest collision energy studied the rigid surface result overestimates the value obtained for Ts = 575 K by more than a factor 2 (green triangles vs red squares). Clearly, allowing energy transfer from the molecule to the surface yields a decreased reaction probability. Compared to the 0 K simulations, raising Ts to 575 K enhances the reaction probability by factors of 3–5 for incidence energies between 1.75 and 2.25 eV, emphasizing the need to model the effect of Ts for this important reaction. We emphasize again that the HD-NNP allows us to study these effects and take them into account in comparing with experiment even for low collision energies where reaction probabilities are on the order 10–5–10–4, which would at present have been impossible with AIMD simulations.

In summary, we have tested the applicability of a HD-NNP for describing the interaction of a molecule with a thermally distorted surface, taking the N2 + Ru(0001) system, which is a paradigm of activated dissociative chemisorption, as an example. The HD-NNP fitted to RPBE-DFT calculations on this system accurately describes the lattice dynamics within the employed slab model of Ru(0001), the electronic and mechanical coupling between molecular and surface atom motion, and the interaction of the molecule with the metal surface. Furthermore, the HD-NNP allows sticking probabilities as small as 10–5 to be computed accurately in QCT calculations simulating both molecular and surface atom motion. Sticking probabilities computed with the HD-NNP are in good agreement with experimental molecular beam sticking probabilities accompanied by a solid error analysis. The HD-NNP tested here allows accurate calculation of reaction probabilities for cases where these probabilities are too small for AIMD to be applicable.