A seven-dimensional quantum dynamics study of the dissociative chemisorption of H2O on Cu(111): effects of azimuthal angles and azimuthal angle-averaging

A seven-dimensional quantum dynamics study for the dissociative chemisorption of H2O on Cu(111) is reported, using the time-dependent wave-packet approach.

The interaction of molecular species with metal surfaces plays a pivotal role in understanding heterogeneous catalysis. Dissociative chemisorption, which features the surface assisted dissociation of a gas-phase molecule into chemisorbed products, is oen the initial and rate-limiting step in many industrial heterogeneous processes such as steam reforming, ammonia synthesis and the water-gas shi reaction. 1 Due to its signicance, tremendous progress has been made for the dissociative chemisorption of molecules in the past two decades. [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16] A large number of molecular beam experiments have revealed the nonstatistical dynamics including mode specicity and bond selectivity in dissociative chemisorption processes. [3][4][5][6]9,[12][13][14][15][16] However, theoretically it is still challenging to investigate those processes at a full-dimensional quantum mechanical level, in particular for the polyatomic molecules dissociating on surfaces, mainly owing to the difficulties in constructing reliable high-dimensional potential energy surfaces (PESs) and developing quantum mechanical methodologies. 17 Such full-dimensional quantum dynamics calculations were limited to diatomic molecules dissociating on metal surfaces, [18][19][20][21][22][23][24][25][26][27][28][29][30][31][32] such as state-of-the-art dynamics calculations for prototypical surface reactions involving H 2 , although many problems are not yet resolved, e.g., the effect of phonons. 18,19,22,27,32 Recently, much more attention has been focused on polyatomic molecules, such as water and methane. The ab initio molecular dynamics approach, which computes the potential energies and forces along the trajectory on the y, has been employed to study dissociative chemisorption dynamics. [33][34][35] Nonetheless, such calculations are computationally too expensive, and they also neglect quantum effects. Although quasiclassical trajectory (QCT) calculations based on global PESs are considerably less expensive and shed valuable light on the dynamics, 36 the accurate full quantum characterization of dissociative chemisorption of polyatomic molecules is highly desirable in view of potentially important quantum effects such as tunneling, zero-point energy and resonances.
The dissociative chemisorption of water on transition-metal surfaces is an essential fundamental step of steam reforming and the water-gas shi reaction. 37 The dynamics of this gas-surface reaction is quantum mechanical in nature due to large zero-point energies and high barriers, where tunneling is important. Since a total of 9 degrees (9D) of freedom on the rigid surface should be considered, a fully coupled 9D quantum treatment of this gas-surface reaction is still formidable. Tiwari and co-workers employed the pseudo-diatomic molecule model to investigate the dynamics for the H 2 O/Cu(111) system, which was based on a three dimensional LEPS PES. 38 Guo and coworkers carried out six-dimensional (6D) quantum dynamics calculations to investigate the mode specicity and bond selectivity for the systems of H 2 O/HOD on a rigid at Cu(111) surface, based on their 6D PESs tted by permutationally invariant polynomials. [39][40][41] The lateral coordinates (X and Y) and the azimuthal angle are xed at the transition state geometry, without considering the effect of impact sites and surface corrugation. The rst molecular beam experiment on the dissociative chemisorption of D 2 O on Ni(111) revealed a strong enhancement in reactivity upon excitation of the antisymmetric stretch of D 2 O. 16 The theoretical results obtained using the sixdimensional quantum approach by Guo and co-workers discussed above semiquantitatively reproduced the observed mode specicity. 16 Recently, Farjamnia and Jackson employed a fully quantum model to investigate the dissociative chemisorption of water on Ni(111), 42 based on the reaction path Hamiltonian (RPH). 43 Although the RPH approach is full-dimensional, it is approximate and only accurate near the minimum energy path. A quantitative description of a gas-surface reaction can only be achieved with a fully coupled quantum mechanical approach based on an accurate global PES.
Very recently, Jiang and Guo performed quasiclassical trajectory and (6 + 1)D quantum dynamics calculations to explore the dynamics of the dissociative chemisorption of D 2 O on rigid Ni(111), using a new nine-dimensional PES developed by a permutation invariant polynomial-neural network (PIP-NN) method. 44,45 In their (6 + 1)D quantum model, only six coordinates were treated explicitly, but the azimuthal angle was xed at the value of the corresponding saddle point of specic sites, thus the at surface approximation was invoked. The inuence of impact sites and incident angles was examined, indicating that the reactivity depends on the site of impact, which is determined by the topography of the PES, but not the static barrier height alone. This suggests that the simple energy shiing site-averaging model employed in many previous studies, 35,46 based on the assumption that dissociation probabilities at different impact sites have the same energy dependence, but varying with the barrier height, should not be quite reliable. A more rigorous site-averaging treatment is to average reduced-dimensional xed-site dissociation probabilities over multiple symmetric sites, as was reported recently by our group for HCl/DCl dissociating on Au(111) 28,29 and H 2 dissociating on Cu(111). 27 The validity of this site-averaging approximation should generally hold in many molecule-surface interactions such as the dissociative chemisorption of water on metal surfaces. Therefore, one can eventually obtain the fully coupled quantum dynamical probabilities by averaging the sevendimensional (7D) xed-site results, instead of directly doing 9D calculations, which is currently computationally more formidable.
In this article, we report the rst 7D quantum dynamics calculations for the dissociative chemisorption of H 2 O on a corrugated, rigid Cu(111) surface, which are represented in molecular coordinates (shown in Fig. 1) by a 7D time-dependent wave packet method. The translational coordinates of the H 2 O moiety in the XY plane are xed at a specic site of impact. The rst globally accurate PES in full 9 dimensions for water dissociative chemisorption on Cu(111) is developed by neural network (NN) tting to roughly 80 000 DFT energy points, employing the C 3v symmetry (the P3m1 plane group symmetry 47 ) of the Cu(111) surface (shown in Fig. S1 †). This t results in an overall very small root mean square error (RMSE) of only 9.0 meV, but is signicantly smaller (6.0 meV) for energy points below 2.0 eV relative to the H 2 O + Cu(111) asymptote, representing the unprecedented tting accuracy for PESs of polyatomic-surface reactions. The 6D site-specic calculations were carried out for various azimuthal angles as well, in an attempt to explore whether there exists the similar averaging approach for the azimuthal angle, as the site-averaging approximation for lateral coordinates (X and Y) that has been veried in gassurface reactions by our group (details of the PES and timedependent wave packet calculations are given in the ESI †).
Contour plots of the PES as a function of the vertical coordinate Z and the distance between the dissociating H atom and the center of mass of nondissociated OH r 1 are illustrated in Fig. 2(a)-(d) for dissociations xed at the transition state (TS), top, bridge and hcp sites, respectively, with other coordinates xed at the corresponding saddle-point geometries, which are also shown in the inset. Overall, these xed-site L-shape reaction paths are quite smooth. The optimized transition-state (TS) geometry by DFT calculations for the dissociation moves towards the hcp site, with one of the O-H bonds stretching out to 1.51Å and the other keeping around the equilibrium value of 0.982Å, which are in accord with results obtained on the PES. As shown, the locations of the saddle points on the four impact sites indicate that the dissociation of H 2 O on Cu(111) is a typical late barrier reaction. The static barrier heights for the xed TS, top, bridge and hcp sites are 1.08 eV, 1.17 eV, 1.13 eV and 1.09 eV, respectively, suggesting that dissociation on the xed TS site and hcp site is more favorable. Because the overall behavior of the contour plots for the xed fcc site is nearly the same as that for the hcp site despite the slightly higher (0.01 eV) barrier height of the fcc site, we do not present the results for the fcc site here.
Based on the newly constructed 9D PES, we carried out the rst 7D time-dependent wave packet study for the dissociative chemisorption of H 2 O on Cu(111), with the translational coordinates of the H 2 O moiety (X and Y) xed at a specic site of impact. These calculations give an excellent opportunity to quantitatively investigate the inuence of the azimuthal angle on the dissociation probabilities and the validity of the sixdimensional quantum model, which neglects two lateral surface coordinates and the azimuthal angle as was done previously. [39][40][41] Fig . 3 shows the seven-dimensional dissociation probabilities for H 2 O initially in the ground rovibrational state for the xed TS, top, bridge and hcp sites, together with the sixdimensional results with their azimuthal angles a xed at the corresponding saddle-point geometries. The 7D dissociation probability for the top site is about 10 times smaller than those for the other sites, although the static barrier height for the top site is just slightly higher than those for the other sites, indicating the signicance of the dynamical investigations as has been found in diatomic molecule-surface reactions. [27][28][29] On the whole, the 7D and 6D dissociation probabilities increase steadily with increasing kinetic energies; however, signicant differences are seen between the 7D and 6D results. The 6D dissociation probabilities with xed azimuthal angles are all considerably larger than the 7D results in the whole energy region, in particular for the TS site and hcp site. As seen, the 7D probabilities for the TS site and hcp site are roughly 2.5 times and twice as small as the corresponding 6D probabilities at a kinetic energy of 2.0 eV, but these factors become much larger in low energy regions. Similar trends are also found for the other two xed sites, with the 6D probabilities larger than the 7D probabilities by roughly 30% and 70% at a kinetic energy of 2.0 eV, respectively, for the top and bridge sites.  To show more clearly the comparison of 6D and 7D dissociation probabilities in very low energy regions below the classical barrier height, we depicted the corresponding results in a logarithmic scale down to 10 À6 in Fig. 4. We can see obviously the differences between the 6D and 7D results in very low energy regions become more signicant, with the 7D dissociation probabilities generally smaller than the corresponding 6D results by a factor of around 9. As a result, the six-dimensional model neglecting the azimuthal angle in quantum dynamical calculations may not reach the quantitative accuracy as compared with the seven-dimensional results.
To better understand the inuence of azimuthal angles on the dynamical results for the title molecule-surface reaction, we present in Fig. 5(a)-(d) the static barrier heights for different azimuthal angles on the PES at xed TS, top, bridge and hcp sites, respectively. The site-specic barrier heights for the azimuthal angle a varying from 0 to 360 with an interval of 10 , with other coordinates xed at the corresponding saddlepoint geometries are given. Overall, the barrier heights as a function of azimuthal angle oscillate with some periodic properties due to the symmetry of the PES. The barrier heights for the TS, bridge and hcp sites considerably change with varying azimuthal angles, with the smallest magnitude just above 1.0 eV and largest one around 2.5 eV. In contrast, as shown in Fig. 5(b), the barrier height for the top site uctuates very slightly with the azimuthal angle and nearly maintains a magnitude of around 1.2 eV in the entire angle region. It is interesting then to see the site-specic dissociation probabilities resulting from different azimuthal angles. Fig. 6 shows the six-dimensional dissociation probabilities for several xed azimuthal angles as well as the seven-dimensional results for H 2 30 ], respectively, for the TS, top, bridge and hcp sites, which correspond to the order of decreasing barrier heights, as shown in Fig. 5. As expected, we can see large differences between the 6D dissociation probabilities with various xed azimuthal angles and 7D dissociation probabilities. As seen from Fig. 6(a) for the TS site, the 6D dissociation probability increases gradually and the threshold shis to lower energy as the azimuthal angle a increases from 230 to 280 , which corresponds to the barrier height decreasing from roughly 2.5 eV to 1.1 eV in Fig. 5(a). In addition, the magnitude of dissociation probability varies from roughly 0.04 to 0.76 at a kinetic energy of 2.0 eV, differing by a factor of 18 times, which implies that the differences among 6D dissociation probabilities with different azimuthal angles can be very signicant, probably due to the large differences of barrier heights. The results for the bridge and hcp sites shown in Fig. 6(c) and (d) demonstrate similar features, indicating that the 6D results depend strongly on the azimuthal angle. However, the 6D dissociation probabilities for the top site in Fig. 6(b) do not change that much as compared with the other three sites, with the magnitude varying from 0.016 to 0.026 at a kinetic energy of 2.0 eV. This is quite understandable as the Fig. 4 Comparisons of 7D dissociation probabilities and 6D dissociation probabilities down to 10 À6 in a logarithmic scale.  differences of barrier heights with different azimuthal angles for the top site are very small, as shown in Fig. 5(b). Interestingly, the magnitudes of the 7D dissociation probabilities for the four xed sites are close to their corresponding medium magnitudes in the 6D results. Then an intriguing question is whether we can nd a similar averaging approach for the azimuthal angle in the H 2 O/Cu(111) system, i.e., accurately reproducing the 7D dissociation probability by averaging over the 6D azimuthal angle-xed results, as has been veried in the site-averaging approach for the surface lateral coordinates (X and Y) in diatomic molecule-surface reactions. [27][28][29] The 6D azimuthal angle-averaged dissociation probabilities are compared with the 7D dissociation probabilities for H 2 O initially in the ground rovibrational state at the xed TS, top, bridge and hcp sites in Fig. 7. As shown, the 18 angle-and 36angle averaged results are basically identical, indicating that the 6D angle-averaged results have already converged with the 18 angles included in the averaging, ranging from 0 to 360 . It is interesting that the overall behavior of the angle-averaged dissociation probabilities resemble those of the 7D dissociation probabilities. Compared with the 6D dissociation probabilities with the azimuthal angle xed at the corresponding saddlepoint geometries, also shown in the blue curves in Fig. 7, the agreement between the angle-averaged and the 7D dissociation probabilities is signicantly improved.
In particular for the top site, as shown in Fig. 7(b), the agreement between the angle-averaged dissociation probabilities and the 7D results is excellent and impressive in the entire kinetic energy region of [1.3, 2.0] eV. For the TS and hcp sites, the angle-averaged results can accurately reproduce the 7D results as well, except for a kinetic energy higher than 1.6 eV. The angle-averaged results become slightly smaller than the 7D results when the kinetic energy exceeds 1.6 eV, but the agreement between them is still good. Although this agreement for the bridge site is less satisfactory, with the angle-averaged probabilities smaller than the 7D results by a factor of roughly 18%, the former can reproduce the latter in the threshold region. The excellent agreement achieved between the 6D angleaveraged and 7D probabilities for the top site results from the relatively weak dependence of the 6D probabilities on azimuthal angles, as discussed above. Overall, the 7D dissociation probabilities for xed sites can be well reproduced by averaging the 6D results over 18 azimuthal angles, in particular in low kinetic energy regions.
To summarize, in this study we carried out the rst sevendimensional quantum dynamics calculations for the dissociative chemisorption of H 2 O on a rigid Cu(111) surface with H 2 O xed at TS, top, bridge and hcp sites, based on an accurate fulldimensional (9D) PES newly developed by neural network tting to roughly 80 000 DFT points. The current 7D quantum dynamical model with one more degree of freedom coupled, i.e., the azimuthal angle, represents much more computational effort compared with the previous 6D quantum model which neglects the azimuthal angle and surface lateral coordinates. It gives a good chance to quantitatively investigate the dependence of azimuthal angles on the dissociation probabilities, as well as to identify the quantitative relationship between the 7D and 6D dissociation probabilities. The calculated 7D dissociation probabilities at the TS, top, bridge and hcp sites are quite different from the corresponding 6D probabilities with various xed azimuthal angles, indicating the great importance of the 7D quantum dynamical investigations. A new azimuthal angleaveraging approach is proposed that one can accurately reproduce the 7D probability by averaging the 6D azimuthal angle xed probabilities over 18 or 36 angles, ranging from 0 to 360 , in particular in low kinetic energy regions, because very satisfactory agreement has been achieved between the 7D probabilities and the 18/36 azimuthal angle-averaged 6D probabilities for H 2 O initially in the ground rovibrational state. The validity of this new approach for excited rovibrational states of H 2 O will be veried in the near future. Furthermore, in principle we are capable of calculating 7D dissociation probabilities for many sites of impact, which can be employed in the site-averaging approach to accurately approximate the full-dimensional (9D) dissociation probability for the dissociative chemisorption of H 2 O on Cu(111). Fig. 7 The angle-averaged dissociation probabilities by averaging the 6D results over 18 and 36 azimuthal angles, 6D dissociation probabilities with the azimuthal angle fixed at the saddle-point geometry, and 7D dissociation probabilities for the fixed TS (a), top (b), bridge (c) and hcp (d) sites.