Accurate Simulations of the Reaction of H2 on a Curved Pt Crystal through Machine Learning.

Theoretical studies on molecule-metal surface reactions have so far been limited to small surface unit cells due to computational costs. Here, for the first time molecular dynamics simulations on very large surface unit cells at the level of density functional theory are performed, allowing a direct comparison to experiments performed on a curved crystal. Specifically, the reaction of D2 on a curved Pt crystal is investigated with a neural network potential (NNP). The developed NNP is also accurate for surface unit cells considerably larger than those that have been included in the training data, allowing dynamical simulations on very large surface unit cells that otherwise would have been intractable. Important and complex aspects of the reaction mechanism are discovered such as diffusion and a shadow effect of the step. Furthermore, conclusions from simulations on smaller surface unit cells cannot always be transfered to larger surface unit cells, limiting the applicability of theoretical studies of smaller surface unit cells to heterogeneous catalysts with small defect densities.

Heterogeneous catalysis is vitally important to many industrial processes. To improve these processes, fundamental insights can be gained by performing molecule–metal surface reaction (MMSR) experiments and simulations. For instance, the shape of the metal surface is important for the overall reactivity of an MMSR as different surface facets can yield different elementary reaction rates. Fortunately, the overall reactivity tends to be dominated by one or only a few rate-controlling states such as the dissociative chemisorption transition state (TS) at specific metal surface facets.[1−4] Especially defects such as steps and kinks are often dominant locations for reactions in industrial processes.[5] However, research involving MMSRs tends to focus on surfaces with small unit cells,[6] which are not always the same as the surfaces that are relevant to heterogeneous catalysis. Computational research is particularly limited by the unit cell size of the investigated surface because the computational cost scales rapidly with the cell size and concomitant number of atoms. Experiments, on the contrary, can more easily investigate MMSRs on large unit cells, e.g., by employing well-defined stepped surfaces or curved crystals.[7−12]

Fortunately, advances in machine learning have enabled previously intractable computational studies. For example, high-dimensional neural network potentials (HDNNPs)[13] can be employed to accurately describe reactive scattering of diatomic molecules from flat metal surfaces, while explicitly modeling surface atom motion.[14−17] Furthermore, HDNNPs allow the use of more accurate but also computationally more demanding density functionals (DFs) such as meta-generalized gradient approximation (MGGA) DFs,[17] which are currently intractable for density functional molecular dynamics (DFMD) studies of MMSRs. Furthermore, polyatomic molecules reacting on metal surfaces can also be accurately described, e.g., CO2 + Ni(100)[18] and CHD3 + Cu(111).[19] With respect to the metal surface, HDNNPs are accurate for describing not only flat surfaces with small unit cells but also nanoparticles[20] and clusters.[21]

However, so far, reactive scattering of molecules from metal surfaces with large unit cells has not been computationally investigated at the level of density functional theory (DFT).[22] The large amount of atoms in such surfaces prohibits the use of DFMD. Furthermore, developing an HDNNP for a large surface unit cell is not straightforward because a large DFT training set is required, which, although cheaper than DFMD, is still a computationally expensive endeavor. Also, such an HDNNP would be able to describe only a single surface facet, precluding a direct comparison with experiments on curved crystals. Recently, it has been shown that the embedded atom neutral network (EANN) approach[23] can be transferred between different surface facets by describing reactive scattering of H2 + Cu(100), Cu(111), Cu(110), and Cu(211)[24] and of CH4 + Ir(111) and Ir(332).[25] Unfortunately, its applicability and accuracy with respect to surface unit cells considerably larger than those included in the training data remain unclear, which is also true for other approaches such as kernel-based regression,[26−28] reactive force fields,[29−33] and deep neural networks.[34,35] Another complicating factor is that neural networks are generally terrible at extrapolating, where the extrapolation is a consequence of evaluating structures in a phase space that is not included in the training data of the HDNNP, although it should also be noted that an HDNNP has been found to be able to extrapolate in the case of protonated water clusters, but still with an accuracy one magnitude lower than in the interpolation regime, i.e., for structures in the phase space that has been included in the development of the HDNNP.[36] For a further overview of machine learning approaches for MMSRs, see the reviews in refs (37−47). In this work, I will show that it is possible to develop an accurate HDNNP that is trained on smaller surface unit cells and that can be transferred to considerably larger surface unit cells.

For the reaction of low-energy H2 on stepped Pt surfaces, there has been a long-standing lack of clarity with regard to the reaction mechanism: Does H2 react directly on the step without major trapping, or is the reaction mediated by trapping?[48−52] Recently, Lent et al. concluded from reactive scattering of H2 from a curved Pt crystal that low-energy H2 always reacts on the step, either directly on impact or via dynamical trapping in the cusp of the step toward the top step edge, but without long-range diffusion.[53] In subsequent experimental work by Jansen and Juurlink, the step sticking cross section has been determined.[54] Interestingly, Jansen and Juurlink found that at a low incidence energy the step sticking cross section was dependent on the step density. The reason for this step density dependence was unknown, and it was suggested that theoretical dynamical studies are required to understand the dependence. Therefore, in this work, the HDNNP is developed specifically for the reaction of H2 on a large variety of Pt surface facets, including large surface unit cells. Because high-fidelity experimental sticking probability data are available for the dissociation of H2 on a curved Pt crystal and molecular beam experiments remain the gold standard for benchmarking theory involving MMSRs to experiment,[6] this reaction will also serve as an excellent benchmark for the approach presented here.

The HDNNP is trained on Pt surfaces containing (111) terraces and (100) steps, e.g., Pt(211). Specifically, the training data include the (111), (211), (533), (322), (755), (433), and (977) surfaces, adding up to ∼88000 structures (see also section S1 of the Supporting Information). Note that the stepped surfaces are similar, where the step is always the same [i.e., a (100)-like step] but where the (111) terrace length varies. Furthermore, one important detail of the construction of the HDNNP is that surface atom motion is omitted from the training data and concomitantly in the molecular dynamics (MD) simulations, but it is also expected that this omission does not influence the results presented in this work (see section S4). As I will show, the HDNNP is accurate not only for the surface unit cells included in the training but also, more importantly, for quasi-classical trajectory (QCT)[55] MD simulations performed for larger surface unit cells not included in the training. This makes (a large number of) MD simulations on very large surface unit cells at the accuracy of DFMD tractable. In this work, QCT-MD simulations are performed up to Pt(171515), which is twice as large as the largest surface unit cell included in the training data [i.e., Pt(977)] and contains 240 “unique” surface atoms in the MD simulations.

The HDNNP reproduces DFT calculations excellently (Figure a), where the HDNNP reproduces 99.3% of the DFT energies within chemical accuracy (i.e., 4.2 kJ/mol) and 99.98% within 2-fold chemical accuracy, yielding an RMSE of 0.024 kJ mol–1 atom–1 for the energies. Furthermore, in panels b and c of Figure , QCT-MD simulations for H2 + Pt(111) and Pt(211) with the HDNNP yield sticking probabilities in good agreement with MD simulations employing corrugation reducing procedure (CRP)[58,59] potential energy surfaces (PESs).[56,57] Also, the outstanding fit quality of the HDNNP is underlined by the excellent agreement with DFMD calculations for H2 + Pt(211) using the same computational setup (Figure c). Note that the DFT calculations on which the CRP PESs are based employed computational setups slightly different than those of the DFT calculations in this work on which the HDNNP is based. Nevertheless, the same DF is used (i.e., PBEα57-vdW-DF2[56,60,61]) and, therefore, PESs based on DFT calculations from either this work or that of Ghassemi and co-workers[56,57] should yield similar sticking probabilities. Moreover, the PBEα57-vdW-DF2 DF has been shown to be chemically accurate for H2 + Pt(111) and Pt(211).[56,57] Because these two surfaces can be considered the two extremes of the probed surfaces in this work in terms of the step density, it is likely that the employed DF is chemically accurate for all surfaces considered in this work, which, as I will show below, is found to be the case. For technical details regarding the MD simulations, see section S2.

Figure 1

(a) Absolute error of the energies in the training and test data predicted by the HD-NNP compared to DFT calculations. The dashed line indicates chemical accuracy, i.e., 4.2 kJ/mol. (b and c) Sticking probability of H2 on Pt(111) and Pt(211), respectively, as a function of incidence energy. The orange circles and green diamonds indicate results obtained from the HDNNP and DFMD, respectively. The blue squares indicate results obtained from the CRP PESs of refs (56) and (57). Note that the DFT calculations used for the CRP PESs employ computational setups slightly different than those in this work, but the same DF, i.e., PBEα57-vdW-DF2.

Figure a compares the computed and measured[54] sticking probability of D2 on stepped Pt surfaces as a function of step density. Throughout this work, theoretical results correspond to a specific Miller index and concomitant step density, whereas experimental results correspond to an ensemble of facets, yielding an average step density. Fortunately, the experimental variance is sufficiently low to compare experimental results obtained from an ensemble of facets with theoretical results obtained from specific facets.[8] At intermediate (2.9 kJ/mol) and high (13.9 kJ/mol) incidence energies, the HDNNP yields sticking probabilities in good agreement with experiment. At a low incidence energy (0.7 kJ/mol), the HDNNP seems to overestimate the sticking probability considerably, but the agreement is actually chemically accurate for a wide variety of step densities and incidence energies (Figure S3). Because the sticking probability decreases rapidly at a low incidence energy, even a small shift along the incidence energy axis can yield large differences in the sticking probability, as is the case here. It should also be emphasized that this shift is considerably smaller than 4.2 kJ/mol; i.e., the experimental sticking probability is reproduced well within chemical accuracy (Figure S3). Furthermore, previous computational work on H2/D2 + Pt(211)[57] also yielded sticking probabilities considerably higher than those from experiment at a low incidence energy. Possibly the employed DF or the QCT approach is the cause of the overestimation, but additional future work is required to understand this discrepancy. Nevertheless, the discrepancy shown here between experiment and theory is not an intrinsic part of the HDNNP, nor is it large. Therefore, I can consider the HDNNP to be not only accurate in reproducing DFT calculations but also accurate in reproducing experimental results. This fact is important because even though step densities of <0.53 nm–1 have not been included in the training of the HDNNP, the HDNNP still yields accurate results for low step densities that would have been intractable to simulate with DFMD and difficult to include in the training data of the HDNNP.

Figure 2

(a) Sticking probability of D2 on stepped Pt surfaces as a function of step density for several incidence energies. The filled (empty) circles with solid (dashed) lines indicate theory (experiment). The incidence energies shown are 0.7 kJ/mol (blue), 2.9 kJ/mol (orange), and 13.9 kJ/mol (green). Note that surfaces with a step density of <0.53 nm–1 are not included in the training data of the HD-NNP. (b) Same as panel a but showing the step sticking cross section (computed with eq ) instead of the sticking probability.

Jansen and Juurlink have also determined the step sticking cross section of the (100)-like step in a curved Pt crystal.[54] At a low incidence energy, the sticking probability of D2 on Pt presumably originates solely from the steps and only from molecules impacting the surface on or close to the step[53] and rapidly decays with the incidence energy. Therefore, site specific sticking probabilities (i.e., S0step and S0terrace) can be determined from the overall sticking probability S0 as a function of incidence energy Ei as follows:[54]where a–d are parameters fitted to the sticking probability and fstep and fterrace are the fractions of the surface unit cell covered in steps and terraces, respectively. Then, the step sticking cross section (AstepS0step) can be determined as follows:[54]where wstep is the width of the unit cell along the step edge and dstep is the step density. Note that the site specific sticking probabilities obtained from eq are employed only when calculating the step sticking cross section with eq (e.g., in Figure b). In all other cases, the location at the moment of reaction (i.e., r = 1.5 Å) is taken from the MD simulations and assigned to the step (terrace) as indicated by the (non)shaded areas in Figure .

Figure 3

(a–f) Initial location (i.e., t = 0, left panels) and the location at the moment of reaction (defined as r = 1.5 Å, right panels) of D2 reacting on Pt(433) for several incidence energies. The colors indicate the probability density, where integration over the entire unit cell yields unity. The circles indicate the top layer atoms, where the black circles indicate step atoms and the red and pink circles indicate terrace atoms. Furthermore, the pink (red) circles indicate that the top site (does not) lies in the shadow of the step top sites. The shaded area associated with the step is indicated by the black lines, whereas the nonshaded area is associated with the terrace. (g) Side and (h) top views of Pt(433). The black, red, and pink atoms correspond to the circles in panels a–f. Moreover, the gray atoms correspond to sublayer terrace atoms.

A comparison between the theoretical and experimental step sticking cross section (Figure b) yields results similar to those for the sticking probability. At low incidence energies, the agreement is qualitative, and at higher incidence energies, the agreement is improved. In general, the agreement for the step sticking cross section is lower than for the sticking probability, which is caused by the slight overestimation of the sticking probability at low incidence energies. Because the model assumes that at a low incidence energy D2 sticks only near the step, an overestimation of the sticking probability at a low incidence energy causes an overall overestimation of the step sticking cross section, even at higher incidence energies (see also section S5). Fortunately, again this is not an intrinsic problem of the HDNNP (vide supra) and the HDNNP can reproduce experimental trends.

Interestingly, Jansen and Juurlink observed a dependence of the step sticking cross section on the step density at low incidence energies, for which the reason was unclear.[54] They mentioned that, among others, short-range diffusion of D2 on the surface might play a role in the step density dependence but that they had no way of confirming this hypothesis. With the simulations presented in this work, the reaction dynamics can be investigated to test this hypothesis, which I will do now. Panels a–f of Figure show the initial location (beginning of the simulation, i.e., t = 0; panels a, c, and e) and the location at the moment of reaction (defined as r = 1.5 Å; panels b, d, and f) of reacting D2 on Pt(433), and panels g and h of Figure show side and top views, respectively, of the Pt(433) surface. Furthermore, Figure shows the step/terrace ratio of the initial and reaction locations as a function of the step density. At the lowest incidence energy (0.7 kJ/mol), the initial location is slightly more concentrated at the step, but most of the distribution is delocalized (Figures a and 4a). When the incidence energy increases, the initial location is increasingly more concentrated at either the step or near a top layer surface atom (Figures c,e and 4b,c). In contrast, the location at the moment of reaction is always strongly concentrated at the top step edge and top layer surface atoms (Figures b,d,f and 4); i.e., the reactive sites are the top step edge and the terrace top sites. Moreover, it is clear that dynamical trapping occurs not only through the cusp near the step, as was previously thought,[50,53] but also on the entire surface [which is also confirmed visually for a few trajectories on Pt(877)]. In short, at a low incidence energy a considerable amount of diffusion of D2 can occur prior to the dissociation, whereas at a high incidence energy only a small amount of diffusion is observed, independent of the step density.

Figure 4

Ratio with which D2 reacts on the step instead of the terrace as a function of the step density for several incidence energies. The blue circles indicate the initial location (i.e., t = 0), whereas the orange circles indicate the location at the moment of reaction (i.e., r = 1.5 Å). The black dashed line indicates the statistical ratio of the unit cell. The error bars represent 68% confidence intervals.

The observations regarding the diffusion of D2 and the reactive sites also offer an explanation for why at a low incidence energy the step sticking cross section depends on the step density. Because the molecule can travel across the terrace and find a favorable location with respect to dissociation (i.e., the top step edge or a terrace top site), the terrace contributes to the overall reactivity, even at a low incidence energy, and effectively increases the step sticking cross section in Figure b computed with eq when the step density decreases. Furthermore, even at the lowest incidence energy, the reaction partly takes place on the terrace instead of the step (Figure ), increasing the step sticking cross section even further (see also section S5). Also, the reactivity of the terrace top sites scales inversely with the step density, again increasing the step sticking cross section for lower step densities. In contrast, the minimum barrier height and geometry on the step, which is an indicator of the reactivity of the step in the absence of diffusion, are hardly affected by the terrace length. For instance, the minimum barrier height on Pt(877) is only 2.2 kJ/mol lower than on Pt(533) (see also section S3). As such, the approximate model overestimates the step sticking cross section (especially for larger surface unit cells) because it relies on the assumption that at a low incidence energy all reaction occurs near the step and that diffusion does not play a significant role. In fact, the step sticking cross section at the moment of reaction even seems visually to be independent of the step density (Figure S5).

As stated earlier, Lent et al. concluded that at a low incidence energy D2 always reacts directly on the step or via dynamical trapping in the cusp of the step toward the top step edge, but without long-range diffusion.[53] The significance of the diffusion of D2 across the entire surface in the reaction mechanism found here might lead one to think that both are possible: D2 reacts directly or via long-range diffusion (i.e., trapping-mediated) on the step. However, because the diffusion occurs within a rather short time (<20 ps) and length (<100 Å) scale before the molecule either reacts or desorbs (see also Figure S2), the reaction is not trapping-mediated because this would involve considerably longer time and length scales. Furthermore, Lent et al. also concluded that the reaction at a low incidence energy occurs only on the step. However, reaction takes place at the step and on the terrace. Therefore, I conclude that the reaction of low-energy D2 occurs both on the step and the terrace (although still predominantly on the step), either directly on impact or via short-range diffusion across the terrace.

Another intriguing aspect of the reaction mechanism is that the terrace reactivity at a low incidence energy is not uniform due to a shadow effect of the step. Panels a–c of Figure show that the terrace top sites that lie on the line perpendicular to the step edge “in the shadow” of the step top sites are less reactive than the terrace top sites that do not lie in the shadow of the step top sites [see Figure for a visualization of which terrace top sites are in the shadow of the step top sites; here, the pink (red) terrace circles/atoms (do not) lie in the shadow of the step top sites, i.e., the black circles/atoms]. Although still present for larger surface unit cells, this effect is considerably less pronounced when the terrace length is increased. Also, at higher incidence energies this shadow effect is still noticeable for large step densities (Figure b,c). The shadow effect is caused by a considerable difference between the top site barriers that lie in or out of the shadow of the step top sites on smaller stepped surface unit cells, whereas the difference in barrier heights and concomitantly the shadow effect largely disappear on larger stepped surface unit cells (Figure d). Furthermore, as mentioned above, in general the terrace top site barrier heights decrease with an increase in terrace length (Figure d). As such, one can argue that extrapolation of conclusions from reaction dynamics on small stepped surface unit cells might not always be valid in describing catalysts with low defect densities (i.e., low step density), even though doing so is the present norm, requiring simulations on large unit cells instead.

Figure 5

(a–c) Ratio of D2 reacting near terrace top sites that lie in the shadow of the top layer step atoms compared to the total sticking probability at the terrace. Blue, orange, and green indicate incidence energies of (a) 0.7, (b) 2.9, and (c) 13.9 kJ/mol, respectively. The black dashed line indicates the statistical ratio of the unit cell. The error bars represent 68% confidence intervals. (d) Terrace top site barriers of the dissociation of H2 on several stepped Pt surfaces. Several step densities are shown (see the legend; the corresponding Miller index is indicated in parentheses). The empty circles are top sites that lie in the shadow of the step top sites, whereas the solid circles do not. The first and last terrace top sites correspond to the sites next to the bottom and top step edge, respectively (see also Figure g,h).

In summary, in this work I perform for the first time a large number of MD simulations at the level of DFT that can be directly compared to a molecular beam experiment performed on a curved crystal, i.e., the reaction of D2 on a curved Pt crystal. The MD simulations are made tractable by developing an HDNNP with DFT data consisting of H2 interacting with a flat (111) surface and with several smaller stepped surface unit cells with (100) steps and (111) terraces [from (211) to (977)]. The developed HDNNP can accurately describe not only the reaction of H2 on surface unit cells that have been included in the training data but also on considerably larger surface unit cells that have not been included in the training data, allowing dynamical simulations on very large surface unit cells that otherwise would have been intractable. Moreover, our understanding of the reaction mechanisms of H2 on stepped Pt surfaces is improved. One of the reaction pathways was thought to be dynamical trapping in the cusp of the step before the molecule is attracted by the step and can dissociate, but here it is observed that the dynamical trapping occurs on the entire terrace before the molecule dissociates on the step or, more rarely, on a terrace top site. Furthermore, previous experimental work determined the step sticking cross section, where a dependence of the step sticking cross section on the step density of the surface was observed. With the help of the MD simulations performed in this work, the underlying cause for the step density dependence is identified. At a low incidence energy, H2 can easily travel large distances across the surface prior to the dissociation of the molecule, even at stepped surfaces with very large terrace lengths. Due to assumptions made in the approximate model employed to obtain the step sticking cross section, the computed step sticking cross section increases substantially with a decrease in step density, even though the actual cross section observed in the MD simulations seems to be independent of the step density. Furthermore, the step changes the reactivity across the surface in a complex fashion, of which the effect is dependent on both the incidence energy of the molecule and the step density of the surface. Therefore, reaction dynamics obtained on small stepped surface unit cells are not always the same as those on large stepped surface unit cells. Similarly, fundamental understanding of catalysts with a low defect density might not always be possible through simulations or experiments performed on small stepped surface unit cells, especially when diffusion plays an important role, but requires studies on large stepped surface unit cells instead.

Computational Methods

Here, I provide a short summary of the computational details for the DFT calculations used to construct the HDNNP with RuNNer.[38,40,62] All DFT calculations have been performed with a user-modified version of the Vienna Ab-initio Simulation Package (VASP version 5.3.5)[63−68] to allow the use of the PBEα57-vdW-DF2 DF.[56,60,61] A plane wave kinetic energy cutoff of 400 eV has been used. Furthermore, a Γ-centered k-point grid is employed, ranging from 6 × 6 × 1 [Pt(111)] to 4× 6× 1 [Pt(977)] to ensure that a reasonably consistent basis set is employed between the different facets. For the (111) facet, a 3 × 3 supercell is employed, whereas for the stepped surfaces, a 1 × 3 supercell is employed. A vacuum distance of 15 Å and five layers have been used, where the top three layers have been relaxed in all directions.