SCIENCE CHINA Chemistry

. ARTICLES . . SPECIAL ISSUE. Chemical Methodology

January 2014 Vol. 57 No. 1: 147–155 doi: 10.1007/s11426-013-5005-7

Six-dimensional potential energy surface of the dissociative chemisorption of HCl on Au(111) using neural networks LIU TianHui, FU BiNa & ZHANG Dong H.∗ State Key Laboratory of Molecular Reaction Dynamics and Center for Theoretical Computational Chemistry; Dalian Institute of Chemical Physics, Chinese Academy of Sciences, Dalian 116023, China Received July 30, 2013; accepted August 19, 2013; published online October 28, 2013

We constructed a six-dimensional potential energy surface (PES) for the dissociative chemisorption of HCl on Au(111) using the neural networks method based on roughly 70000 energies obtained from extensive density functional theory (DFT) calculations. The resulting PES is accurate and smooth, based on the small fitting errors and good agreement between the fitted PES and the direct DFT calculations. Time-dependent wave packet calculations show that the potential energy surface is very well converged with respect to the number of DFT data points, as well as to the fitting process. The dissociation probabilities of HCl initially in the ground rovibrational state from six-dimensional quantum dynamical calculations are quite diﬀerent from the four-dimensional fixed-site calculations, indicating it is essential to perform full-dimensional quantum dynamical studies for the title molecule-surface interaction system. potential energy surface, DFT calculations, neural networks, wave packet, six-dimensional

1 Introduction The dissociative chemisorption of molecular species on transition-metal surfaces is one of the most important reactions in heterogeneous catalytic processes. It is often the rate-limiting step in most of chemical reactions, and numerous eﬀorts have been devoted to studying this dissociation process both experimentally and theoretically in the last several decades [1–15]. Due to the diﬃculties in constructing reliable potential energy surfaces (PESs) and developing quantum dynamics methodologies, very limited quantum dynamical calculations were carried out for dissociative adsorption reactions, mainly with the chemisorption of the simple H2 molecule on metal surfaces [16–26]. A full-dimensional quantum dynamical description of polyatomic dissociative chemisorption is still extremely challenging. Several quantum dynamical studies of CH4 on Ni surfaces resorted to employ low-dimensional models [27–31], as 15 degrees of freedom should be considered on a rigid surface, rendering it *Corresponding author (email: [email protected])

c Science China Press and Springer-Verlag Berlin Heidelberg 2013

diﬃcult to develop an accurate, global PES and to carry out quantum dynamical calculations. Recently, Guo and co-workers [32–34] investigated the dissociative chemisorption of H2 O on Cu(111) using the reduced-dimensional (sixdimensional) quantum dynamical calculations. An accurate six-dimensional PES was fit by the permutationally invariant polynomial approach of Bowman and co-workers [35,36] based on energy points from density functional theory (DFT) calculations. In this article we report an accurate, six-dimensional PES of the dissociative chemisorption of HCl on Au(111), constructed using neural network fitting to DFT energy points. In early 1990s, Lykke et al. [37] carried out a molecule beamsurface scattering experiment of HCl scattered from Au(111), which indicated a direct inelastic scattering mechanism. Recently, Wodtke and co-workers [38–41] investigated the energy transfer between gas and solid interfaces for this system, and a transition from an electronically adiabatic mechanical mechanism to an electronically nonadiabatic mechanism involving excited electron-hole pairs as the surface temperature increases was reported. To the best of our knowledge, none chem.scichina.com

link.springer.com

148

Liu TH, et al.

Sci China Chem

of theoretical investigations was reported for the title reaction. As a result, we have done this work and present the details of the PES constructed using the neural network approach in this paper. As the motion of surface atoms and excited electron-hole pairs were ignored in the current adiabatic PES, and all the calculations were based on the rigid surface approximation, the current theory is not able to directly describe the experiment. Neural networks (NN) are general fitting methods that can be used, in principle, to fit any function. Over the past two decades, they have been applied to fit many high-dimensional PESs for isolated molecules (up to 30 dimensions for Nmethylacetamide based on DFT (B3LYP) data points) and for molecule-surface interactions (up to 12 dimensions for H2 at Si(100)-(2 × 1) surface based on DFT (LDA) calculations) [42–53]. Scheﬄer and co-workers [54–56] presented the neural network scheme to fit the PESs for H2 dissociation on Pd(100) surface and O2 dissociation on Al(111) surface based on DFT energies, where symmetry-adapted coordinates were employed as the input to the NN fit, representing the symmetry of the surface. However, that approach was complicated by technical diﬃculties in achieving a proper consideration of the symmetry of the surface, and therefore to overcome these limitations they introduced a symmetry-adapted NN representation of PESs for moleculesurface interactions, which was based on a new type of symmetry functions that takes the symmetry of the surface potential exactly into account [57]. The accuracy and eﬃciency of such symmetry-adapted NNs was illustrated by the application to a six-dimensional PES describing the O2 /Al(111) system [57, 58]. Very recently, our group has successfully applied the neural networks to calculate the full-dimensional (six-dimensional) PESs of the prototype four-atom reactions, i.e., OH+H2 and OH+CO reactions [59,60]. The current PES of HCl molecule dissociation on Au (111) using NNs is well converged with respect to the number of DFT energy points, as well as to the fitting process, examined by time-dependent wave packet calculations. The rest of the paper is organized as follows: In Section 2, we give the details for DFT calculations of energy points, and NNs used in our fitting and the fitting procedure. Some properties of the PES are presented and comparisons of dissociation probabilities from six-dimensional and fourdimensional fixed-site quantum dynamical calculations using the time-dependent wave packet (TDWP) method are made in Section 3. The final section contains a brief discussion.

2 Methods

January (2014) Vol. 57 No. 1

were expanded in a planewave basis set [64]. The electron exchange correlation eﬀects were treated within the generalized gradient approximation (GGA) [65], using the Perdew-Wang (PW91) functional [66]. The Au(111) substrates consist of four layers with a 2 × 2 surface unit cell(1/4-ML coverage). A vacuum region between two repeated cells is 16 Å, the Monkhorst-Pack k-points grid mesh is 5 × 5 × 1 [67] and the planewave expansion is truncated at the kinetic energy of 400 eV. The optimized lattice constants for bulk Au is 4.1766 Å in this work, which agrees well with the experimental value (4.0783 Å) [68]. To validate our slab model, we calculated the relative barrier energy Eb and the relative product energy Ep with diﬀerent computation details. Eb and Ep are defined as Eb = b - ∞ and p - ∞ , where b , p , ∞ are the absolute barrier, product and asymptotic energies, respectively. The detailed parameters used in the DFT calculations are shown in Table 1, with the cutoﬀ kinetic energies ranging from 300 eV to 450 eV, the number of layers ranging from four to five, and the size of the supercell ranging from 2 × 2 to 3 × 3. We can see that the convergence of these parameters is good, indicating the reliability of the parameters currently used in this study. Due to the large number of the energy points we needed to obtain an accurate PES for this system, the moderate slab model was found to be the best compromise between eﬃciency and accuracy and was used for the computation of the PES. The climbing-image nudged elastic band (CI-NEB) method was used to determine the transition state (TS) and minimum energy path [69, 70]. Our calculations are considered to be converged when all forces are smaller than 0.01 eV/Å. To fit a PES from a given set of molecular configurations, we employed the feed forward NN with two hidden layers connecting the input layer and output layer, denoted as I − J − K − 1 NN. It has I nodes in the input layer, which equals to the number of degrees of freedom or input cartesian coordinates employed here (see below) for a molecular configuration, and one nodes in the output layer corresponding to the potential energy value of the input configuration. The two hidden layers have J and K neurons, respectively. The output of jth neuron in the first hidden layer is I (w1j,i × xi ) , y1j = f 1 b1j +

j = 1, 2, · · · , J

(1)

i=1

and the output of kth neuron in the second hidden layer is J (w2k, j × y1j ) , y2k = f 2 b2k +

k = 1, 2, · · · , K

(2)

j=1

All the planewave DFT calculations were carried out with the Vienna ab initio simulation package (VASP) [61,62]. The interaction between ionic cores and electrons was described by fully nonlocal optimized projector augmented-wave (PAW) potentials [63], and the Kohn-Sham valence electronic states

and consequently the final output is given by y = b31 +

K k=1

(w31,k × y2k ),

(3)

Liu TH, et al.

Table 1

Sci China Chem

Comparisons of calculated barrier height (Eb ) and relative product

energies (Ep ) of HCl dissociation on Au(111) surface, keeping the 5 × 5 × 1 Monkhorst-Pack k-points grid mesh Computation details

Eb

Ep

Four layers, 1/4 ML, 300 eV

0.639

0.256

Four layers, 1/4 ML, 350 eV

0.655

0.279

Four layers, 1/4 ML, 400 eV

0.661

0.275

Four layers, 1/4 ML, 450 eV

0.664

0.276

Five layers, 1/4 ML, 400 eV

0.593

0.276

Four layers, 1/9 ML, 400 eV

0.631

0.247

where xi (i = 1, · · · , 6) are the Cartesian coordinates of HCl in this HCl/Au(111) system when the surface is fixed, the weights wlj,i connect the ith neuron of (l − 1)th layer and the jth neuron of lth layer, and the biases blj determine the threshold of the jth neuron of lth layer, f 1 and f 2 are transfer functions taken as hyperbolic tangent functions. The full six-dimensional coordinates of the HCl/Au(111) system, namely X, Y, Z, r, θ and φ, are shown in Figure 1 with the surface fixed in the XY plane to describe the molecular configuration, where Z is the distance between the center of mass (COM) of HCl and the surface, r is the bond length of HCl, and θ and φ are two molecular orientation angles. This set of molecular coordinates is also used in the six-dimensional quantum dynamics calculations using the TDWP approach, for which the detailed methodologies and results are expected to be reported in the following work. During the PES fitting, a series of NNs with diﬀerent number of neurons for these two hidden layers was constructed and trained to get an optimal architecture, because the architecture (structure) of NN influences strongly the quality of fit. For a given NN architecture, the weights and biases in Eq. (1)–(3) can be updated through proper NN training using

149

January (2014) Vol. 57 No. 1

the Levenberg-Marquardt algorithm [71]. The root mean square error (RMSE) is used to measure the performance of NN. n 1 RMSE = (Efit − Eab initio )2 n i=1 The Levenberg-Marquardt algorithm used for our NN fitting is an iteratively fitting procedure with high computational cost. The iteration time grows rapidly with respect to the size of fitting samples and network parameters. To accelerate the NN fitting procedure, we divided the data points into several parts according to their locations (entrance valley, interaction region, product valley), trained and tested these parts separately to obtain desired performances. Finally these segmental parts can be connected to yield a global PES. The division of the PES into several parts can substantially reduce the coordinate space for each part, making it feasible to reach desired RMSE with less number of neurons. It not only speeds up the training procedure, but also gives us the flexibility to put more eﬀorts and/or more DFT data points in a region of especial importance to dynamics. It is worthwhile to point out that the potential energies are periodical with respect to the three coordinates, X, Y and φ in this HCl/Au system, especially the X and Y which correspond to the lateral symmetry. Exploiting the C3v symmetry of the Au(111) surface, it is suﬃcient to map the PES only for the molecular configurations inside the irreducible triangle of the surface unit cell which is spanned by the top, hollow and fcc sites, as shown in Figure 2. This approach significantly reduces the number of required electronic structure calculations. We applied the symmetry functions proposed by Bohler et al. [57,72] to the NN fit of the current HCl/Au(111) system, but found that the overall fitting error in term of RMSE increased considerably on the same set of DFT data points for the present system. Instead of this approach, we

ez

Cl

hcp r

ey

Z

θ

top

H

φ

bridge

ex fcc

Ω =120° Au

Figure 1

Molecular coordinates for the dissociative chemisorption of HCl

on Au(111), where (X,Y,Z) are the COM coordinates of HCl, r is the bond

Figure 2

length of HCl, θ is the polar angle, and φ is the azimuthal angle. Here Ω is

red triangle), the atoms in the first (second) layer are shown in white (gray)

The irreducible triangle of Au(111) surface unit cell (shown in

the skew angle and equals 120◦ .

spheres.

Liu TH, et al.

Sci China Chem

included the small range of data points near the boundary outside the triangle into the training set, to get a better fit that resolves the problem of the discontinuity in potential derivatives at the boundary originating from the switch to another site. The energies of these data points are not necessary to be calculated again and can be achieved by copying from the energies of points being equivalent for symmetry reasons inside the irreducible triangle. This procedure will increase 30% of fitting data size, but it is appropriate to generate a continuous PES that is suﬃciently accurate for the dynamics simulations.

January (2014) Vol. 57 No. 1

8000 76393 ab initio points 6000 Number of points

150

4000

2000

0 −0.4

3 Results We use the Cartesian coordinates (X1 ,Y1 , Z1 ) of the H atom and (X2 , Y2 , Z2 ) of the Cl atom as the input coordinates for the NN fit. The Cartesian coordinates are simply transformed from the molecular coordinates (X, Y, Z, r, θ , φ) as follows. X1 = X + r

m2 sin(Ω − φ) sin(θ), m1 + m2 sin(Ω)

m1 sin(Ω − φ) sin(θ), X2 = X − r m1 + m2 sin(Ω) Y1 = Y + r

m2 sin(φ) sin(θ), m1 + m2 sin(Ω)

m1 sin(φ) Y2 = Y − r sin(θ), m1 + m2 sin(Ω) Z1 = Z + r and

m2 cos(θ), m1 + m2

m1 cos(θ) m1 + m2 where the skewing angle Ω = 2π/3. A total of 76393 DFT energy points were calculated. Considering the continuity of potentials at the boundary of the irreducible triangle, the data size was enlarged to 99201 energy points using the method we mentioned above and was included in fitting the PES of the title molecule-surface interaction system. Figure 3 shows the potential energy distribution of the whole data set used in the fitting procedure. One can see that the data points are nearly equally distributed in the energy region of [−0.4, 4.0] eV, with a large number of energy points (roughly 65% of the data set) fall in the energy region of 0.0–2.0 eV above the HCl + Au(111) asymptote, which is the region of most importance to dynamics. To accelerate the NN fitting procedure as was described in detail [59, 60], we divided the data points into four parts: entrance part (I), interaction parts (II and III), and export part (IV), with some overlaps between I and II parts, II and III parts, III and IV parts as shown in Table 2. Roughly 70% of data points were calculated in the interaction parts II and III, and rest of them were distributed in the entrance part and export part where HCl is far away from the surface or it has Z2 = Z − r

Figure 3

0.6

1.6 2.6 Potential energy (eV)

3.6

The distribution of electronic energies of all the DFT data points

relative to the HCl + Au(111) asymptote used for the PES fit. Table 2

The range of the four parts divided and the number of energy

points calculated in each part (Z is the vertical distance between the COM of HCl and the surface, r is the bond length of HCl, and ZH is the vertical distance between the H atom and the surface). The distances are in Å Region

Number of data points

Part I (a)

Z 3.9 and ZH 3.0

15299

Part II (b)

(Z 4.1 or ZH 3.0) and r 1.5

22323

Part III (c) (Z 4.1 or ZH 3.0) and r 2.5 and r 1.4 Part IV (d)

(Z 4.1 or ZH 3.0) and r 2.4

46617 23819

been dissociated. Less neurons are required to fit these parts separately to reach a desired RMSE, leading to a much more eﬃcient fitting without any loss of accuracy. Convergence properties of these parts are tested independently with more eﬀort paid to the important regions of a PES. For each part we obtained three fits with small RMSEs using diﬀerent NN architectures or the same architecture with diﬀeren weights and biases, as shown in Table 3. We can see that the diﬀerences of the RMSEs between current three fits in each part are quite small (< 0.1 meV), with the maximal diﬀerence of 0.4 meV for the RMSEs of the three fits in part III, ranging from 3.94 to 4.31 meV. The combination of the 12 segmental fits shown in Table 3 results in a total of 81 global PESs, one of which is denoted as the NN1 PES here with the NN structures of 6-50-40-1, 6-60-45-1, 6-80-50-1 and 6-65-60-1, and RMSEs of 2.25 meV, 1.18 meV, 4.31 meV and 6.63 meV for part I, II, III and IV, respectively. The dissociation probabilities for HCl in the ground rovibrational state (v = 0, j = 0) were calculated using the six-dimensional TDWP approach on the three global PESs which were randomly selected from resulting 80 PESs except NN1. These probability curves are compared with that calculated on NN1 in Figure 4, in the kinetic energy region of 0.4 eV and 1.6 eV. As displayed in Figure 4, the dissociation probability is a monotonically increasing function as the kinetic energy. The probabilities obtained on the three PESs randomly selected are essentially identical to

Liu TH, et al.

Sci China Chem

151

January (2014) Vol. 57 No. 1

Table 3 Neural network structure parameters and fitting errors (meV) for

Table 4 Neural network structure parameters and fitting errors for NN1

the four segmental parts (shown in parentheses). In each part we used diﬀer-

and NN2 PESs

ent weights and biases to give three fits Structure

Fitting points

Fit 1

Fit 2

Fit 3

NN1 part I

6-50-40-1

15299

2.25

Part I (a)

6-50-40-1 (2.25)

6-50-40-1 (2.28)

6-50-45-1 (2.19)

NN1 part II

6-60-45-1

22323

1.18

Part II (b)

6-60-45-1 (1.18)

6-60-45-1 (1.23)

6-60-45-1 (1.27)

NN1 part III

6-80-50-1

46617

4.31

Part III (c)

6-80-50-1 (4.31)

6-80-50-1 (4.26)

6-90-40-1 (3.94)

NN1 part III

6-65-50-1

23879

6.63

Part IV (d)

6-65-60-1 (6.63)

6-65-60-1 (6.66)

6-65-55-1 (6.69)

NN1 overall 6-50-40-1

12239

2.30

NN2 part I

Dissociation probability

0.5

0.4

a1+b1+c1+d1 NN1 PES a2+b2+c2+d2 a3+b3+c1+d2 a2+b3+c3+d3

4.5

NN2 part II

6-60-45-1

17858

2.84

NN2 part III

6-80-40-1

37293

7.02

NN2 part IV

6-65-60-1

19103

NN2 overall

6.78 5.81

0.3 0.2 0.2 0.1

0 0.4

0.8 1.2 Kinetic energy (eV)

Efit − Eab initio (eV)

0.1

Figure 4

RMSE (meV)

1.6

0

The six-dimensional dissociation probabilities for the dissocia-

−0.1

tive chemisorption of HCl (v = 0, j = 0) calculated on four PESs we selected randomly in the total 81 PESs (a, b, c, and d denote part I, II, III, and IV, −0.2 −0.5

respectively).

Figure 5

2.5 1.5 Potential energy (eV)

3.5

The fitting errors for all the data points in the NN1 PES, as a

function of their corresponding DFT energies relative to the HCl + Au(111) asymptote.

0.5 NN2 NN1 Dissociation probability

that on the NN1 PES, indicating that the PES is well converged with the fitting procedure. The NN1 PES fits well to the DFT energies as can be seen from the fitting errors shown in Table 4, as well as from Figure 5 that presents the fitting errors for all the data points in the NN1 PES as a function of their corresponding energies. It is seen that most of data points (> 97.0%) have the fitting errors that are smaller than 10 meV, but only about 80 energy points are shown to be outliers for which the fitting errors are larger than 50 meV. These energy points with relatively large fitting errors are found to fall in the quite high energy region (> 2.0 eV). The overall RMSE of NN1 is 4.5 meV, but significantly smaller (2.95 meV) for energy points below 2.0 eV relative to HCl + Au(111) asymptote. Furthermore, to check the convergency, 20% of data points were removed randomly in each part of the NN1 PES and another PES, denoted as NN2, was fitted based on the remaining data points. Using the same procedure mentioned for NN1 above, we selected one fit in each of the four segmental parts with a relatively small RMSE, and combine the four fits to generate the global NN2 PES. The structure parameters and fitting errors of the NN2 PES we obtained are shown in Table 4. The dissociation probabilities for the HCl molecule on the ground rovibrational state are calculated on NN1 and NN2 PESs and illustrated in Figure 6. As can be seen from Figure

0.5

0.4

0.3

0.2

0.1

0 0.4

Figure 6

0.8 1.2 Kinetic energy (eV)

1.6

The six-dimensional dissociation probabilities for the dissocia-

tive chemisorption of HCl (v=0, j=0) on Au(111) on the NN1 and NN2 PESs.

6, the probability curve obtained on the NN2 PES, which was fitted on 80% randomly selected points of NN1, agrees very well to that on the NN1 PES, indicating both PESs are well

Liu TH, et al.

Sci China Chem

converged with respect to the number of data points as well as the fitting procedure. The minimum energy path of the dissociative chemisorption of HCl on Au(111) obtained by the CI-NEB method are shown in Figure 7, with the energies relative to the reactants HCl + Au(111) asymptote (zero potential energy). This path was discretized by eight images between the reactant and product. The reactant HCl is 6.0 Å above the bridge site of the Au(111) surface, and is almost parallel to the surface. There exists a van der waals well of about −0.04 eV, where HCl is around the equilibrium distance of 1.30 and 3.55 Å above the surface, but it moves slightly away from the bridge site. The HCl molecule at the optimized transition state (TS) moves towards the fcc site, with the bond length stretched to 1.917 Å. The structure parameters optimized at the transition state is X = 1.018 Å, Y = 0.576 Å, Z = 2.394 Å, r = 1.917 Å, θ = 48.3 (degree), φ = 25.1(degree), and the barrier height is 0.661 eV, which agrees well with the transition state geometry obtained from the fitted PES of X = 1.049 Å, Y = 0.575 Å, Z = 2.388 Å, r = 1.92 Å, θ = 48.4(degree), φ = 27.4 (degree) with the barrier height of 0.658 eV. The product configuration corresponds to the position where H is adsorbed on a fcc site, and Cl is adsorbed on a nearest-neighbor fcc site. The minimum energy path of the fitted PES and the DFT energy points calculated at the same geometries are also compared and shown in Figure 8. As we can see, the very good agreement is achieved. Therefore, the current PES constructed using NN for the dissociative adsorption of HCl on Au(111) is suﬃciently accurate to be used in investigating the quantum reaction dynamics. The contour plots of the NN1 PES as a function of Z and r for HCl fixed at bridge, hollow and top sites are illustrated in Figure 9(a), (b), and (c), respectively, with other coordinates (θ and φ) optimized. Overall, the L-shape reaction paths for the three fixed sites are quite smooth. The equilibrium distance of the free HCl molecule is 1.28 Å on the fitted NN1 PES, in good agreement with the DFT result. As seen, the

0.7

Energy (eV)

0.5

0.3

0.1

−0.1

Figure 7

Reaction coordinate

Minimum energy pathway for the dissociative chemisorption of

HCl on Au(111) obtained using CINEB method with eight images.

January (2014) Vol. 57 No. 1

0.7

0.5 Energy (eV)

152

DFT FIT

0.3

0.1

−0.1 Reaction coordinate

Figure 8 Comparison of the fitted NN1 PES and corresponding DFT energies along the minimum energy pathway.

reaction barriers are located around r = 3.76 bohr and Z = 4.39 bohr, r = 3.76 bohr and Z = 4.47 bohr, and r = 3.62 bohr and Z = 4.76 bohr for the bridge, hollow and top sites, respectively, indicative of the late barriers for the fixed sites of the title molecule-surface interaction system. Hence, the vibrational excitation of HCl is more eﬃcient in enhancing the reactivity of the reaction than the translational motion. The reaction barrier height for the bridge site is lowest (0.64 eV), compared to those for the hollow site of 0.78 eV and the top site of 0.83 eV, indicating the dissociation on the fixed bridge site is much more favored. Because the overall behavior of the contour plot for the fixed fcc site is similar to that of the hollow site despite of the slightly higher (0.02 eV) barrier height of the fcc site, we do not present the plot for the fcc site here. Figure 10 shows the dissociation probability for HCl initially in the ground rovibrational state (v = 0, j = 0) calculated using the six-dimensional (6D) TDWP method, together with the four-dimensional (4D) results for which HCl is fixed at the bridge, fcc, hollow and top cites, in the kinetic energy region of 0.4 eV and 1.6 eV. Overall, all the dissociation probabilities are smooth functions of the kinetic energy, and rise steadily with increase of the kinetic energy. The 4D dissociation probability for the fixed bridge site is much larger than the rest of three 4D probabilities and 6D probability in the entire energy region, presenting the lowest threshold (roughly 0.56 eV) compared to other results. This is quite understandable that the bare barrier height for the fixed bridge site is lowest as seen from Figure 9. The overall behavior of 4D results for the fixed fcc site is very similar to that of the hollow site, with the probability curves rise from the threshold of around 0.75 eV. The 4D dissociation probability for the hollow site is slightly larger than that for the fcc site just above the threshold resulting from the slightly lower barrier for the hollow cite as mentioned above; however, when the kinetic energy increases further, the two results become comparable in the medium energy region, and the dissociation probability for the fcc cite turns to be larger than the hollow site at the kinetic energy above 1.3 eV. It is interesting that the

Liu TH, et al.

Sci China Chem

6.0

0.8 (a)

bridge site

5.5

Dissociation probability

0.7 * TS

0.7

r (bohr)

0.6

4.0 3.5

0.6

3.0 2.5 2.0 3.5

4D bridge site 4D fcc site 4D hollow site 4D top site ×100 6D

Eb = 0.64 eV

5.0 4.5

153

January (2014) Vol. 57 No. 1

4.0

4.5

5.0 5.5 Z (bohr)

6.0

6.5

0.6

0.4

0.2

7.0 0 0.4

6.0 (b)

hollow site

5.5

Eb = 0.78 eV

5.0

Figure 10

0.8 1.2 Kinetic energy (eV)

1.6

Comparison of dissociation probabilities for the scattering of

HCl (v = 0, j = 0) from Au(111) between the six-dimensional calculation and four-dimensional fixed-site calculations.

lity is smaller than that of the bridge site in the entire energy region.

0.8

0.7 * TS 0.7

4.0

0.8

r (bohr)

4.5

3.5 3.0

4 Conclusions

2.5 2.0 3.5

4.0

4.5

5.0 5.5 Z (bohr)

6.0

6.5

7.0

6.0 (c)

top site

5.5

Eb = 0.83 eV

5.0

3.5

0.8 * TS 0.8

0.9

4.0

0.9

r (bohr)

4.5

3.0 2.5 2.0 4.0

Figure 9

4.5

5.0

5.5 Z (bohr)

6.0

6.5

7.0

(a) Contour plot of the NN1 PES as a function of Z and r for the

bridge site, with other coordinates optimized. The contours are relative to the HCl + Au(111) asymptote with an interval of 0.1 eV. The reaction barriers are also indicated by red star. (b) Same as (a) except for the hollow site. (c) Same as (a) except for the top site.

dissociation probability for the fixed top site is about 100–200 times smaller than the other results, though the bare barrier height for the top site is only 0.05 eV higher than the hollow site and 0.19 eV higher than the bridge site. Again as seem from Figure 10, the 6D probability is quite diﬀerent from the results calculated by reduced dimensional (4D) fixed-site approach, indicating the importance of full-dimensional calculations. The threshold of the 6D dissociation probability is about 0.02 eV higher than that of the 4D bridge site, in accord with the diﬀerence of the bare barrier heights of them shown in Figures 8 and 9, and the magnitude of the 6D probabi-

To summarize, we present a six-dimensional PES for the dissociative chemisorption of HCl on Au(111) based on NN fitting of more than 70000 DFT energy points, which was calculated using VASP with the moderate slab model. The resulting PES is accurate, with a very small RMSE of about 4.5 meV up to the potential energy of 4.0 eV relative to the HCl + Au(111) asymptote, and very good agreement between the PES and direct DFT calculations is achieved for the minimum energy path. Furthermore, the PES is well converged with respect to the number of data points as well as the fitting procedure, from extensive quantum dynamical calculations using the TDWP method. The four-dimensional fixedsite calculations for HCl initially in the ground rovibrational state show that the dissociation process is much more favored on the fixed bridge site, due to the lowest reaction barrier compared to hollow, fcc, and top sites. The dissociation probability for the top site is several hundred times smaller than those for the other sites, though the reaction barrier is not significantly (only 0.19 eV) higher than the bridge site. The comparison of the six-dimensional and four-dimensional fixed-site dissociation probabilities shows large diﬀerences, indicating it is very important to do full-dimensional calculations for the title molecule-surface interaction system. This is the first global PES for the dissociative chemisorption of HCl on Au(111), which is suﬃciently accurate to be used for the quantum dynamical investigations, and the following work for the detailed quantum dynamical calculations on this PES is underway and expected to be reported later. This work was supported by the National Natural Science Foundation of China (91221301 and 90921014), the Chinese Academy of Sciences, and Ministry of Science and Technology of China.

154

Liu TH, et al.

Sci China Chem

1 Beck RD, Maroni P, Papageorgopoulos DC, Dang TT, Schmid MP, Rizzo TR. Vibrational mode-specific reaction of methane on a nickel surface. S cience, 2003, 302: 98–100 2 Gates GA, Darling GR, Holloway S. A theoretical study of the vibrational excitation of NO/Ag(111). J Chem Phys, 1994, 101: 6281–6288 3 Juurlink LBF. Comparative study of C-H stretch and bend vibrations in methane activation on Ni(100) and Ni(111). Phys Rev Lett, 2005, 94: 208303(1–4) 4 Juurlink LBF, Killelea DR, Utz AL. State-resolved probes of methane dissociation dynamics. Prog Surf Sci, 2009, 84: 69–134 5 Killelea DR, Campbell VL, Shuman NS, Utz AL. Bond-selective control of a heterogeneously catalyzed reaction. S cience, 2008, 319: 790– 793 6 Kimman J, Rettner C, Auerbach D, Barker J, Tully J. Correlation between kinetic-energy transfer to rotation and to phonons in gas-surface collisions of NO with Ag(111). Phys Rev Lett, 1986, 57: 2053–2056 7 Kurten T, Biczysko M, Rajam¨aki T, Laasonen K, Halonen L. Computational study of the adsorption energetics and vibrational wavenumbers of NH3 adsorbed on the Ni(111) surface. J Phys Chem B, 2005, 109: 8954–8960 8 Lee MB, Yang QY, Ceyer ST. Dynamics of the activated dissociative chemisorptions of CH4 and implication for the pressure gap in catalysis: A molecular beam-high resolution electron energy loss study. J Chem Phys, 1987, 87: 2724–2741 9 Luntz AC. CH4 dissociation on Ni(100): Comparison of a direct dynamical model to molecular beam experiments. J Chem Phys, 1995, 102: 8264–8269 10 Luntz AC. Chemistry. Toward vibrational mode control in catalysis. Science, 2003, 302: 70–71 11 Maroni P, Papageorgopoulos DC, Sacchi M, Dang TT, Beck RD, Rizzo TR. State-resolved gas-surface reactivity of methane in the symmetric C-H stretch vibration on Ni(100). Phys Rev Lett, 2005, 94: 246104 (1-4) 12 Michelsen HA, Rettner CT, Auerbach DJ, Zare RN. Eﬀect of rotation on the translational and vibrational energy dependence of the dissociative adsorption of D2 on Cu(111). J Chem Phys, 993, 98: 8294–8307 13 Rettner CT, Michelsen HA, and Auerbach DJ. Quantum-state-specific dynamics of the issociative adsorption and associative desorption of H2 at a Cu(111) surface. J Chem Phys, 1995, 102: 4625–4641 14 Schmid MP, Maroni P, Beck RD, Rizzo TR. Surface reactivity of highly vibrationally excited molecules prepared by pulsed laser excitation: CH4 (2ν3 ) on Ni(100). J Chem Phys, 2002, 117: 8603–8606 15 Smith RR, Killelea DR, DelSesto DF, Utz AL. Preference for vibrational over translational energy in a gas-surface reaction. Science, 2004, 304: 992–995 16 Dai J, Light JC. Six dimensional quantum dynamics study for dissociative adsorption of H2 on Cu(111) surface. J Chem Phys, 1997, 107: 1676–1679 17 Dai J, Zhang JZH. Quantum adsorption dynamics of a diatomic molecule on surface: Four-dimensional fixed-site model for H2 on Cu(111). J Chem Phys, 1995, 102: 6280–6289 18 Dai J, Sheng J, Zhang JZH. Symmetry and rotational orientation eﬀects in dissociative adsorption of diatomic molecules on metals: H2 and HD on Cu(111). J Chem Phys, 1994, 101: 1555–1563 19 Gr¨uneich A, Cruz AJ, Jackson B. The dynamics of H2 dissociation on Cu and Ni surfaces. Mixed quantum-classical studies with all degrees of freedom. J Chem Phys, 1993, 98: 5800–5808 20 Kroes GJ, Wiesenekker G, Baerends EJ, Mowrey RC, Neuhauser D. Dissociative chemisorption of H2 on Cu(100): A four-dimensional study of the eﬀect of parallel translational motion on the reaction dynamics. J Chem Phys, 1996, 105: 5979–5998 21 Saalfrank P, Miller WH. Time-independent quantum dynamics for diatom-surface scattering. J Chem Phys, 1993, 98: 9040–9052

January (2014) Vol. 57 No. 1

22 Sheng J, Zhang JZH. Quantum dynamics studies of adsorption and desorption of hydrogen at a Cu(111) surface. J Chem Phys, 1993, 99: 1373–1381 23 Vincent JK, Olsen RA, Kroes GJ, Luppi M, Baerends EJ. Sixdimensional quantum dynamics of dissociative chemisorption of H2 on Ru(0001). J Chem Phys, 2005, 122: 44701(1–8) 24 Pijper E, Kroes GJ, Olsen RA, Baerends EJ. Reactive and diﬀractive scattering of H2 from Pt(111) studied using a six-dimensional wave packet method. J Chem Phys, 2002, 117: 5885 25 Pijper E, Kroes GJ, Olsen RA, Baerends EJ. Dissociative and diﬀractive scattering of H2 from Pt(111): A four-dimensional quantum dynamics study. J Chem Phys, 2002, 116: 9435–9448 26 Eichler A, Hafner J, Groß A, Scheﬄer M. Trends in the chemical reactivity of surfaces studied by ab initio quantum-dynamics calculations. Phys Rev B, 1999, 59: 13297-13300 27 Krishnamohan GP, Olsen RA, Kroes GJ, Gatti F, Woittequand S. Quantum dynamics of dissociative chemisorption of CH4 on Ni(111): Influence of the bending vibration. J Chem Phys, 2010, 133: 144308(1–15) 28 Nave S, Jackson B. Methane dissociation on Ni(111) and Pt(111): Energetic and dynamical studies. J Chem Phys, 2009, 130: 054701(1–14) 29 Nave S, Tiwari AK, Jackson B. Methane dissociation and adsorption on Ni(111), Pt(111), Ni(100), Pt(100), and Pt(110)-(1×2): Energetic study. J Chem Phys, 2010, 132: 054705(1–12) 30 Jackson B, Nave S. The dissociative chemisorption of methane on Ni(100): Reaction path description of mode-selective chemistry. J Chem Phys, 2011, 135: 114701(1–12) 31 Jiang B, Liu R, Li J, Xie D, Yang M, Guo H. Mode selectivity in methane dissociative chemisorption on Ni(111). Chem Sci, 2013, 4: 3249–3254 32 Jiang B, Xie D, Guo H. Vibrationally mediated bond selective dissociative chemisorptions of HOD on Cu(111). Chem Sci, 2013, 4: 503–508 33 Jiang B, Li J, Xie D, Guo H. Eﬀects of reactant internal excitation and orientation on dissociative chemisorption of H2 O on Cu(111): Quasiseven-dimensional quantum dynamics on a refined potential energy surface. J Chem Phys, 2013, 138: 044704(1–9) 34 Jiang B, Ren X, Xie D, Guo H. Enhancing dissociative chemisorption of H2 O on Cu(111) via vibrational excitation. Proc Natl Acad Sci USA, 2012, 109: 10224–10227 35 Bowman JM, Czak´o G, Fu B. High-dimensional ab initio potential energy surfaces for reaction dynamics calculations. Phys Chem Phys Chem, 2011, 13: 8094–8111 36 Braams BJ, Bowman JM. Permutationally invariant potential energy surfaces in high dimensionality. Int Rev Phys Chem, 2009, 28: 577– 606 37 Lykke KR, Kay BD. Rotationally inelastic gas-surface scattering: HCl from Au(111). J Chem Phys, 1990, 92: 2614–2623 38 Cooper R, Rahinov I, Yuan C, Yang X, Auerbach DJ, Wodtke AM. Efficient translational excitation of a solid metal surface: State-to-state translational energy distributions of vibrational ground state HCl scattering from Au(111). J Vac Sci Technol A, 2009, 27: 907–912 39 Rahinov I, Cooper R, Yuan C, Yang X, Auerbach DJ, Wodtke AM. Eﬃcient vibrational and translational excitations of a solid metal surface: State-to-state time-of-flight measurements of HCl(v = 2, J = 1) scattering from Au(111). J Chem Phys, 2008, 129: 214708(1–16) 40 Ran Q, Matsiev D, Auerbach DJ, Wodtke AM. Observation of a change of vibrational excitation mechanism with surface temperature: HCl collisions with Au(111). Phys Rev Lett, 2007, 98: 237601(1–4) 41 Ran Q, Matsiev D, Auerbach DJ, Wodtke AM. Direct translation-tovibrational energy transfer of HCl on gold: Measurement of absolute vibrational excitation probabilities. Nucl Instrum Meth B, 2007, 258: 1–6 42 Behler J, Parrinello M. Generalized neural-network representation of high-dimensional potential-energy surfaces. Phys Rev Lett, 2007, 98:

Liu TH, et al.

Sci China Chem

146401(1–4) 43 Blank TB, Brown SD, Calhoun AW, Doren DJ. Neural network models of potential energy surfaces. J Chem Phys, 1995, 103: 4129–4137 44 Brown DFR, Gibbs MN, Clary DC. Combining ab initio computations, neural networks, and diﬀusion Monte Carlo: An eﬃcient method to treat weakly bound molecules. J Chem Phys, 1996, 105: 7597–7604 45 Darley MG, Handley CM, Popelier PLA. Beyond point charges: Dynamic polarization from neural net predicted multipole moments. J Chem T heory Comput, 2008, 4: 1435–1448 46 Jose KV, Artrith N, Behler J. Construction of high-dimensional neural network potentials using environment-dependent atom pairs. J Chem Phys, 2012, 136: 194111(1–15) 47 Le HM, Huynh S, Raﬀ LM. Molecular dissociation of hydrogen peroxide (HOOH) on a neural network ab initio potential surface with a new configuration sampling method involving gradient fitting. J Chem Phys, 2009, 131: 014107(1–10) 48 Le HM, Raﬀ LM. Cis→trans, trans→cis isomerizations and N-O bond dissociation of nitrous acid (HONO) on an ab initio potential surface obtained by novelty sampling and feed-forward neural network fitting. J Chem Phys, 2008, 128: 194310(1–11) 49 Le HM, Raﬀ LM. Molecular dynamics investigation of the bimolecular reaction BeH + H2 →BeH2 + H on an ab initio potential-energy surface obtained using neural network methods with both potential and gradient accuracy determination. J Phys Chem A, 2009, 114: 45–53 50 Malshe M, Raﬀ LM, Rockley MG, Hagan M, Agrawal PM, Komanduri R. Theoretical investigation of the dissociation dynamics of vibrationally excited vinyl bromide on an ab initio potential-energy surface obtained using modified novelty sampling and feedforward neural networks. II. Numerical application of the method. J Chem Phys, 2007, 127: 134105(1–7) 51 Manzhos S, Wang X, Dawes R, Carrington T. A nested moleculeindependent neural network approach for high-quality potential fits. J Phys Chem A, 2005, 110: 5295–5304 52 Pukrittayakamee A, Malshe M, Hagan M, Raﬀ LM, Narulkar R, Bukkapatnum S, Komanduri R. Simultaneous fitting of a potentialenergy surface and its corresponding force fields using feedforward neural networks. J Chem Phys, 2009, 130: 134101(1–10) 53 Raﬀ LM, Malshe M, Hagan M, Doughan DI, Rockley MG, Komanduri R. Ab initio potential-energy surfaces for complex, multichannel systems using modified novelty sampling and feedforward neural networks. J Chem Phys, 2005, 122: 84104(1–16) 54 Behler J, Delley B, Lorenz S, Reuter K, Scheﬄer M. Dissociation of O2 at Al(111): The role of spin selection rules. Phys Rev Lett, 2005, 94: 036104(1–4) 55 Lorenz S, Groß A, Scheﬄer M. Representing high-dimensional potential-energy surfaces for reactions at surfaces by neural networks. Chem Phys Lett, 2004, 395: 210–215

January (2014) Vol. 57 No. 1

155

56 Lorenz S, Scheﬄer M, Gross A. Descriptions of surface chemical reactions using a neural network representation of the potential-energy surface. Phys Rev B, 2006, 73: 115431(1–13) 57 Behler J, Lorenz S, Reuter K. Representing molecule-surface interactions with symmetry-adapted neural networks. J Chem Phys, 2007, 127: 014705(1–11) 58 Behler J, Reuter K, Scheﬄer M. Nonadiabatic eﬀects in the dissociation of oxygen molecules at the Al(111) surface. Phys Rev B, 2008, 77: 115421(1–16) 59 Chen J, Xu X, Xu X, Zhang DH. Communication: An accurate global potential energy surface for the OH + CO → H + CO2 reaction using neural networks. J Chem Phys, 2013, 138: 221104(1–4) 60 Chen J, Xu X, Xu X, Zhang DH. A global potential energy surface for the H2 + OH ↔ H2 O + H reaction using neural networks. J Chem Phys, 2013, 138: 154301(1–8) 61 Kresse G, Furthm¨uller J. Eﬃciency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. Comp Mater Sci, 1996, 6: 15–50 62 Kresse G and Furthm¨uller J. Eﬃcient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys Rev B, 1996, 54: 11169–11186 63 Bl¨ochl PE. Projector augmented-wave method. Phys Rev B, 1994, 50: 17953–17979 64 Kresse G, Joubert D. From ultrasoft pseudopotentials to the projector augmented-wave method. Phys Rev B, 1999, 59: 1758–1775 65 Perdew JP, Burke K, Ernzerhof M. Generalized gradient approximation made simple. Phys Rev Lett, 1996, 77: 3865–3868 66 Perdew JP, Chevary JA, Vosko SH, Jackson KA, Pederson MR, Singh DJ, Fiolhais C. Atoms, molecules, solids, and surfaces: Applications of the generalized gradient approximation for exchange and correlation. Phys Rev B, 1992, 46: 6671–6687 67 Monkhorst HJ, Pack JD. Special points for Brillouin-zone integrations. Phys Rev B, 1976, 13: 5188–5192 68 Haynes WM. CRC Handbook o f Chemistry and Physics. 93rd edition (Internet version 2013). Boca Raton, Fl: CRC Press/Taylor and Francis 69 Henkelman G, Uberuaga BP, J´onsson H. A climbing image nudged elastic band method for finding saddle points and minimum energy paths. J Chem Phys, 2000, 113: 9901–9904 70 Jonsson H, Mills G, Jacobsen KW. Nudged elastic band method f or f inding minimum energy paths o f transitions. Classical and Quantum Dynamics in Condensed Phase S imulations. Eds. Berne BJ, Ciccotti G, Coker DF. Singapore: World Scientific, 1998 71 Hagan MT, Menhaj MB. Training feedforward networks with the Marquardt algorithm. IEEE T rans Neural Netw, 1994, 5: 989–993 72 Behler J. Neural network potential-energy surfaces in chemistry: A tool for large-scale simulations. Phys Chem Chem Phys, 2011, 13: 17930– 17955

. ARTICLES . . SPECIAL ISSUE. Chemical Methodology

January 2014 Vol. 57 No. 1: 147–155 doi: 10.1007/s11426-013-5005-7

Six-dimensional potential energy surface of the dissociative chemisorption of HCl on Au(111) using neural networks LIU TianHui, FU BiNa & ZHANG Dong H.∗ State Key Laboratory of Molecular Reaction Dynamics and Center for Theoretical Computational Chemistry; Dalian Institute of Chemical Physics, Chinese Academy of Sciences, Dalian 116023, China Received July 30, 2013; accepted August 19, 2013; published online October 28, 2013

We constructed a six-dimensional potential energy surface (PES) for the dissociative chemisorption of HCl on Au(111) using the neural networks method based on roughly 70000 energies obtained from extensive density functional theory (DFT) calculations. The resulting PES is accurate and smooth, based on the small fitting errors and good agreement between the fitted PES and the direct DFT calculations. Time-dependent wave packet calculations show that the potential energy surface is very well converged with respect to the number of DFT data points, as well as to the fitting process. The dissociation probabilities of HCl initially in the ground rovibrational state from six-dimensional quantum dynamical calculations are quite diﬀerent from the four-dimensional fixed-site calculations, indicating it is essential to perform full-dimensional quantum dynamical studies for the title molecule-surface interaction system. potential energy surface, DFT calculations, neural networks, wave packet, six-dimensional

1 Introduction The dissociative chemisorption of molecular species on transition-metal surfaces is one of the most important reactions in heterogeneous catalytic processes. It is often the rate-limiting step in most of chemical reactions, and numerous eﬀorts have been devoted to studying this dissociation process both experimentally and theoretically in the last several decades [1–15]. Due to the diﬃculties in constructing reliable potential energy surfaces (PESs) and developing quantum dynamics methodologies, very limited quantum dynamical calculations were carried out for dissociative adsorption reactions, mainly with the chemisorption of the simple H2 molecule on metal surfaces [16–26]. A full-dimensional quantum dynamical description of polyatomic dissociative chemisorption is still extremely challenging. Several quantum dynamical studies of CH4 on Ni surfaces resorted to employ low-dimensional models [27–31], as 15 degrees of freedom should be considered on a rigid surface, rendering it *Corresponding author (email: [email protected])

c Science China Press and Springer-Verlag Berlin Heidelberg 2013

diﬃcult to develop an accurate, global PES and to carry out quantum dynamical calculations. Recently, Guo and co-workers [32–34] investigated the dissociative chemisorption of H2 O on Cu(111) using the reduced-dimensional (sixdimensional) quantum dynamical calculations. An accurate six-dimensional PES was fit by the permutationally invariant polynomial approach of Bowman and co-workers [35,36] based on energy points from density functional theory (DFT) calculations. In this article we report an accurate, six-dimensional PES of the dissociative chemisorption of HCl on Au(111), constructed using neural network fitting to DFT energy points. In early 1990s, Lykke et al. [37] carried out a molecule beamsurface scattering experiment of HCl scattered from Au(111), which indicated a direct inelastic scattering mechanism. Recently, Wodtke and co-workers [38–41] investigated the energy transfer between gas and solid interfaces for this system, and a transition from an electronically adiabatic mechanical mechanism to an electronically nonadiabatic mechanism involving excited electron-hole pairs as the surface temperature increases was reported. To the best of our knowledge, none chem.scichina.com

link.springer.com

148

Liu TH, et al.

Sci China Chem

of theoretical investigations was reported for the title reaction. As a result, we have done this work and present the details of the PES constructed using the neural network approach in this paper. As the motion of surface atoms and excited electron-hole pairs were ignored in the current adiabatic PES, and all the calculations were based on the rigid surface approximation, the current theory is not able to directly describe the experiment. Neural networks (NN) are general fitting methods that can be used, in principle, to fit any function. Over the past two decades, they have been applied to fit many high-dimensional PESs for isolated molecules (up to 30 dimensions for Nmethylacetamide based on DFT (B3LYP) data points) and for molecule-surface interactions (up to 12 dimensions for H2 at Si(100)-(2 × 1) surface based on DFT (LDA) calculations) [42–53]. Scheﬄer and co-workers [54–56] presented the neural network scheme to fit the PESs for H2 dissociation on Pd(100) surface and O2 dissociation on Al(111) surface based on DFT energies, where symmetry-adapted coordinates were employed as the input to the NN fit, representing the symmetry of the surface. However, that approach was complicated by technical diﬃculties in achieving a proper consideration of the symmetry of the surface, and therefore to overcome these limitations they introduced a symmetry-adapted NN representation of PESs for moleculesurface interactions, which was based on a new type of symmetry functions that takes the symmetry of the surface potential exactly into account [57]. The accuracy and eﬃciency of such symmetry-adapted NNs was illustrated by the application to a six-dimensional PES describing the O2 /Al(111) system [57, 58]. Very recently, our group has successfully applied the neural networks to calculate the full-dimensional (six-dimensional) PESs of the prototype four-atom reactions, i.e., OH+H2 and OH+CO reactions [59,60]. The current PES of HCl molecule dissociation on Au (111) using NNs is well converged with respect to the number of DFT energy points, as well as to the fitting process, examined by time-dependent wave packet calculations. The rest of the paper is organized as follows: In Section 2, we give the details for DFT calculations of energy points, and NNs used in our fitting and the fitting procedure. Some properties of the PES are presented and comparisons of dissociation probabilities from six-dimensional and fourdimensional fixed-site quantum dynamical calculations using the time-dependent wave packet (TDWP) method are made in Section 3. The final section contains a brief discussion.

2 Methods

January (2014) Vol. 57 No. 1

were expanded in a planewave basis set [64]. The electron exchange correlation eﬀects were treated within the generalized gradient approximation (GGA) [65], using the Perdew-Wang (PW91) functional [66]. The Au(111) substrates consist of four layers with a 2 × 2 surface unit cell(1/4-ML coverage). A vacuum region between two repeated cells is 16 Å, the Monkhorst-Pack k-points grid mesh is 5 × 5 × 1 [67] and the planewave expansion is truncated at the kinetic energy of 400 eV. The optimized lattice constants for bulk Au is 4.1766 Å in this work, which agrees well with the experimental value (4.0783 Å) [68]. To validate our slab model, we calculated the relative barrier energy Eb and the relative product energy Ep with diﬀerent computation details. Eb and Ep are defined as Eb = b - ∞ and p - ∞ , where b , p , ∞ are the absolute barrier, product and asymptotic energies, respectively. The detailed parameters used in the DFT calculations are shown in Table 1, with the cutoﬀ kinetic energies ranging from 300 eV to 450 eV, the number of layers ranging from four to five, and the size of the supercell ranging from 2 × 2 to 3 × 3. We can see that the convergence of these parameters is good, indicating the reliability of the parameters currently used in this study. Due to the large number of the energy points we needed to obtain an accurate PES for this system, the moderate slab model was found to be the best compromise between eﬃciency and accuracy and was used for the computation of the PES. The climbing-image nudged elastic band (CI-NEB) method was used to determine the transition state (TS) and minimum energy path [69, 70]. Our calculations are considered to be converged when all forces are smaller than 0.01 eV/Å. To fit a PES from a given set of molecular configurations, we employed the feed forward NN with two hidden layers connecting the input layer and output layer, denoted as I − J − K − 1 NN. It has I nodes in the input layer, which equals to the number of degrees of freedom or input cartesian coordinates employed here (see below) for a molecular configuration, and one nodes in the output layer corresponding to the potential energy value of the input configuration. The two hidden layers have J and K neurons, respectively. The output of jth neuron in the first hidden layer is I (w1j,i × xi ) , y1j = f 1 b1j +

j = 1, 2, · · · , J

(1)

i=1

and the output of kth neuron in the second hidden layer is J (w2k, j × y1j ) , y2k = f 2 b2k +

k = 1, 2, · · · , K

(2)

j=1

All the planewave DFT calculations were carried out with the Vienna ab initio simulation package (VASP) [61,62]. The interaction between ionic cores and electrons was described by fully nonlocal optimized projector augmented-wave (PAW) potentials [63], and the Kohn-Sham valence electronic states

and consequently the final output is given by y = b31 +

K k=1

(w31,k × y2k ),

(3)

Liu TH, et al.

Table 1

Sci China Chem

Comparisons of calculated barrier height (Eb ) and relative product

energies (Ep ) of HCl dissociation on Au(111) surface, keeping the 5 × 5 × 1 Monkhorst-Pack k-points grid mesh Computation details

Eb

Ep

Four layers, 1/4 ML, 300 eV

0.639

0.256

Four layers, 1/4 ML, 350 eV

0.655

0.279

Four layers, 1/4 ML, 400 eV

0.661

0.275

Four layers, 1/4 ML, 450 eV

0.664

0.276

Five layers, 1/4 ML, 400 eV

0.593

0.276

Four layers, 1/9 ML, 400 eV

0.631

0.247

where xi (i = 1, · · · , 6) are the Cartesian coordinates of HCl in this HCl/Au(111) system when the surface is fixed, the weights wlj,i connect the ith neuron of (l − 1)th layer and the jth neuron of lth layer, and the biases blj determine the threshold of the jth neuron of lth layer, f 1 and f 2 are transfer functions taken as hyperbolic tangent functions. The full six-dimensional coordinates of the HCl/Au(111) system, namely X, Y, Z, r, θ and φ, are shown in Figure 1 with the surface fixed in the XY plane to describe the molecular configuration, where Z is the distance between the center of mass (COM) of HCl and the surface, r is the bond length of HCl, and θ and φ are two molecular orientation angles. This set of molecular coordinates is also used in the six-dimensional quantum dynamics calculations using the TDWP approach, for which the detailed methodologies and results are expected to be reported in the following work. During the PES fitting, a series of NNs with diﬀerent number of neurons for these two hidden layers was constructed and trained to get an optimal architecture, because the architecture (structure) of NN influences strongly the quality of fit. For a given NN architecture, the weights and biases in Eq. (1)–(3) can be updated through proper NN training using

149

January (2014) Vol. 57 No. 1

the Levenberg-Marquardt algorithm [71]. The root mean square error (RMSE) is used to measure the performance of NN. n 1 RMSE = (Efit − Eab initio )2 n i=1 The Levenberg-Marquardt algorithm used for our NN fitting is an iteratively fitting procedure with high computational cost. The iteration time grows rapidly with respect to the size of fitting samples and network parameters. To accelerate the NN fitting procedure, we divided the data points into several parts according to their locations (entrance valley, interaction region, product valley), trained and tested these parts separately to obtain desired performances. Finally these segmental parts can be connected to yield a global PES. The division of the PES into several parts can substantially reduce the coordinate space for each part, making it feasible to reach desired RMSE with less number of neurons. It not only speeds up the training procedure, but also gives us the flexibility to put more eﬀorts and/or more DFT data points in a region of especial importance to dynamics. It is worthwhile to point out that the potential energies are periodical with respect to the three coordinates, X, Y and φ in this HCl/Au system, especially the X and Y which correspond to the lateral symmetry. Exploiting the C3v symmetry of the Au(111) surface, it is suﬃcient to map the PES only for the molecular configurations inside the irreducible triangle of the surface unit cell which is spanned by the top, hollow and fcc sites, as shown in Figure 2. This approach significantly reduces the number of required electronic structure calculations. We applied the symmetry functions proposed by Bohler et al. [57,72] to the NN fit of the current HCl/Au(111) system, but found that the overall fitting error in term of RMSE increased considerably on the same set of DFT data points for the present system. Instead of this approach, we

ez

Cl

hcp r

ey

Z

θ

top

H

φ

bridge

ex fcc

Ω =120° Au

Figure 1

Molecular coordinates for the dissociative chemisorption of HCl

on Au(111), where (X,Y,Z) are the COM coordinates of HCl, r is the bond

Figure 2

length of HCl, θ is the polar angle, and φ is the azimuthal angle. Here Ω is

red triangle), the atoms in the first (second) layer are shown in white (gray)

The irreducible triangle of Au(111) surface unit cell (shown in

the skew angle and equals 120◦ .

spheres.

Liu TH, et al.

Sci China Chem

included the small range of data points near the boundary outside the triangle into the training set, to get a better fit that resolves the problem of the discontinuity in potential derivatives at the boundary originating from the switch to another site. The energies of these data points are not necessary to be calculated again and can be achieved by copying from the energies of points being equivalent for symmetry reasons inside the irreducible triangle. This procedure will increase 30% of fitting data size, but it is appropriate to generate a continuous PES that is suﬃciently accurate for the dynamics simulations.

January (2014) Vol. 57 No. 1

8000 76393 ab initio points 6000 Number of points

150

4000

2000

0 −0.4

3 Results We use the Cartesian coordinates (X1 ,Y1 , Z1 ) of the H atom and (X2 , Y2 , Z2 ) of the Cl atom as the input coordinates for the NN fit. The Cartesian coordinates are simply transformed from the molecular coordinates (X, Y, Z, r, θ , φ) as follows. X1 = X + r

m2 sin(Ω − φ) sin(θ), m1 + m2 sin(Ω)

m1 sin(Ω − φ) sin(θ), X2 = X − r m1 + m2 sin(Ω) Y1 = Y + r

m2 sin(φ) sin(θ), m1 + m2 sin(Ω)

m1 sin(φ) Y2 = Y − r sin(θ), m1 + m2 sin(Ω) Z1 = Z + r and

m2 cos(θ), m1 + m2

m1 cos(θ) m1 + m2 where the skewing angle Ω = 2π/3. A total of 76393 DFT energy points were calculated. Considering the continuity of potentials at the boundary of the irreducible triangle, the data size was enlarged to 99201 energy points using the method we mentioned above and was included in fitting the PES of the title molecule-surface interaction system. Figure 3 shows the potential energy distribution of the whole data set used in the fitting procedure. One can see that the data points are nearly equally distributed in the energy region of [−0.4, 4.0] eV, with a large number of energy points (roughly 65% of the data set) fall in the energy region of 0.0–2.0 eV above the HCl + Au(111) asymptote, which is the region of most importance to dynamics. To accelerate the NN fitting procedure as was described in detail [59, 60], we divided the data points into four parts: entrance part (I), interaction parts (II and III), and export part (IV), with some overlaps between I and II parts, II and III parts, III and IV parts as shown in Table 2. Roughly 70% of data points were calculated in the interaction parts II and III, and rest of them were distributed in the entrance part and export part where HCl is far away from the surface or it has Z2 = Z − r

Figure 3

0.6

1.6 2.6 Potential energy (eV)

3.6

The distribution of electronic energies of all the DFT data points

relative to the HCl + Au(111) asymptote used for the PES fit. Table 2

The range of the four parts divided and the number of energy

points calculated in each part (Z is the vertical distance between the COM of HCl and the surface, r is the bond length of HCl, and ZH is the vertical distance between the H atom and the surface). The distances are in Å Region

Number of data points

Part I (a)

Z 3.9 and ZH 3.0

15299

Part II (b)

(Z 4.1 or ZH 3.0) and r 1.5

22323

Part III (c) (Z 4.1 or ZH 3.0) and r 2.5 and r 1.4 Part IV (d)

(Z 4.1 or ZH 3.0) and r 2.4

46617 23819

been dissociated. Less neurons are required to fit these parts separately to reach a desired RMSE, leading to a much more eﬃcient fitting without any loss of accuracy. Convergence properties of these parts are tested independently with more eﬀort paid to the important regions of a PES. For each part we obtained three fits with small RMSEs using diﬀerent NN architectures or the same architecture with diﬀeren weights and biases, as shown in Table 3. We can see that the diﬀerences of the RMSEs between current three fits in each part are quite small (< 0.1 meV), with the maximal diﬀerence of 0.4 meV for the RMSEs of the three fits in part III, ranging from 3.94 to 4.31 meV. The combination of the 12 segmental fits shown in Table 3 results in a total of 81 global PESs, one of which is denoted as the NN1 PES here with the NN structures of 6-50-40-1, 6-60-45-1, 6-80-50-1 and 6-65-60-1, and RMSEs of 2.25 meV, 1.18 meV, 4.31 meV and 6.63 meV for part I, II, III and IV, respectively. The dissociation probabilities for HCl in the ground rovibrational state (v = 0, j = 0) were calculated using the six-dimensional TDWP approach on the three global PESs which were randomly selected from resulting 80 PESs except NN1. These probability curves are compared with that calculated on NN1 in Figure 4, in the kinetic energy region of 0.4 eV and 1.6 eV. As displayed in Figure 4, the dissociation probability is a monotonically increasing function as the kinetic energy. The probabilities obtained on the three PESs randomly selected are essentially identical to

Liu TH, et al.

Sci China Chem

151

January (2014) Vol. 57 No. 1

Table 3 Neural network structure parameters and fitting errors (meV) for

Table 4 Neural network structure parameters and fitting errors for NN1

the four segmental parts (shown in parentheses). In each part we used diﬀer-

and NN2 PESs

ent weights and biases to give three fits Structure

Fitting points

Fit 1

Fit 2

Fit 3

NN1 part I

6-50-40-1

15299

2.25

Part I (a)

6-50-40-1 (2.25)

6-50-40-1 (2.28)

6-50-45-1 (2.19)

NN1 part II

6-60-45-1

22323

1.18

Part II (b)

6-60-45-1 (1.18)

6-60-45-1 (1.23)

6-60-45-1 (1.27)

NN1 part III

6-80-50-1

46617

4.31

Part III (c)

6-80-50-1 (4.31)

6-80-50-1 (4.26)

6-90-40-1 (3.94)

NN1 part III

6-65-50-1

23879

6.63

Part IV (d)

6-65-60-1 (6.63)

6-65-60-1 (6.66)

6-65-55-1 (6.69)

NN1 overall 6-50-40-1

12239

2.30

NN2 part I

Dissociation probability

0.5

0.4

a1+b1+c1+d1 NN1 PES a2+b2+c2+d2 a3+b3+c1+d2 a2+b3+c3+d3

4.5

NN2 part II

6-60-45-1

17858

2.84

NN2 part III

6-80-40-1

37293

7.02

NN2 part IV

6-65-60-1

19103

NN2 overall

6.78 5.81

0.3 0.2 0.2 0.1

0 0.4

0.8 1.2 Kinetic energy (eV)

Efit − Eab initio (eV)

0.1

Figure 4

RMSE (meV)

1.6

0

The six-dimensional dissociation probabilities for the dissocia-

−0.1

tive chemisorption of HCl (v = 0, j = 0) calculated on four PESs we selected randomly in the total 81 PESs (a, b, c, and d denote part I, II, III, and IV, −0.2 −0.5

respectively).

Figure 5

2.5 1.5 Potential energy (eV)

3.5

The fitting errors for all the data points in the NN1 PES, as a

function of their corresponding DFT energies relative to the HCl + Au(111) asymptote.

0.5 NN2 NN1 Dissociation probability

that on the NN1 PES, indicating that the PES is well converged with the fitting procedure. The NN1 PES fits well to the DFT energies as can be seen from the fitting errors shown in Table 4, as well as from Figure 5 that presents the fitting errors for all the data points in the NN1 PES as a function of their corresponding energies. It is seen that most of data points (> 97.0%) have the fitting errors that are smaller than 10 meV, but only about 80 energy points are shown to be outliers for which the fitting errors are larger than 50 meV. These energy points with relatively large fitting errors are found to fall in the quite high energy region (> 2.0 eV). The overall RMSE of NN1 is 4.5 meV, but significantly smaller (2.95 meV) for energy points below 2.0 eV relative to HCl + Au(111) asymptote. Furthermore, to check the convergency, 20% of data points were removed randomly in each part of the NN1 PES and another PES, denoted as NN2, was fitted based on the remaining data points. Using the same procedure mentioned for NN1 above, we selected one fit in each of the four segmental parts with a relatively small RMSE, and combine the four fits to generate the global NN2 PES. The structure parameters and fitting errors of the NN2 PES we obtained are shown in Table 4. The dissociation probabilities for the HCl molecule on the ground rovibrational state are calculated on NN1 and NN2 PESs and illustrated in Figure 6. As can be seen from Figure

0.5

0.4

0.3

0.2

0.1

0 0.4

Figure 6

0.8 1.2 Kinetic energy (eV)

1.6

The six-dimensional dissociation probabilities for the dissocia-

tive chemisorption of HCl (v=0, j=0) on Au(111) on the NN1 and NN2 PESs.

6, the probability curve obtained on the NN2 PES, which was fitted on 80% randomly selected points of NN1, agrees very well to that on the NN1 PES, indicating both PESs are well

Liu TH, et al.

Sci China Chem

converged with respect to the number of data points as well as the fitting procedure. The minimum energy path of the dissociative chemisorption of HCl on Au(111) obtained by the CI-NEB method are shown in Figure 7, with the energies relative to the reactants HCl + Au(111) asymptote (zero potential energy). This path was discretized by eight images between the reactant and product. The reactant HCl is 6.0 Å above the bridge site of the Au(111) surface, and is almost parallel to the surface. There exists a van der waals well of about −0.04 eV, where HCl is around the equilibrium distance of 1.30 and 3.55 Å above the surface, but it moves slightly away from the bridge site. The HCl molecule at the optimized transition state (TS) moves towards the fcc site, with the bond length stretched to 1.917 Å. The structure parameters optimized at the transition state is X = 1.018 Å, Y = 0.576 Å, Z = 2.394 Å, r = 1.917 Å, θ = 48.3 (degree), φ = 25.1(degree), and the barrier height is 0.661 eV, which agrees well with the transition state geometry obtained from the fitted PES of X = 1.049 Å, Y = 0.575 Å, Z = 2.388 Å, r = 1.92 Å, θ = 48.4(degree), φ = 27.4 (degree) with the barrier height of 0.658 eV. The product configuration corresponds to the position where H is adsorbed on a fcc site, and Cl is adsorbed on a nearest-neighbor fcc site. The minimum energy path of the fitted PES and the DFT energy points calculated at the same geometries are also compared and shown in Figure 8. As we can see, the very good agreement is achieved. Therefore, the current PES constructed using NN for the dissociative adsorption of HCl on Au(111) is suﬃciently accurate to be used in investigating the quantum reaction dynamics. The contour plots of the NN1 PES as a function of Z and r for HCl fixed at bridge, hollow and top sites are illustrated in Figure 9(a), (b), and (c), respectively, with other coordinates (θ and φ) optimized. Overall, the L-shape reaction paths for the three fixed sites are quite smooth. The equilibrium distance of the free HCl molecule is 1.28 Å on the fitted NN1 PES, in good agreement with the DFT result. As seen, the

0.7

Energy (eV)

0.5

0.3

0.1

−0.1

Figure 7

Reaction coordinate

Minimum energy pathway for the dissociative chemisorption of

HCl on Au(111) obtained using CINEB method with eight images.

January (2014) Vol. 57 No. 1

0.7

0.5 Energy (eV)

152

DFT FIT

0.3

0.1

−0.1 Reaction coordinate

Figure 8 Comparison of the fitted NN1 PES and corresponding DFT energies along the minimum energy pathway.

reaction barriers are located around r = 3.76 bohr and Z = 4.39 bohr, r = 3.76 bohr and Z = 4.47 bohr, and r = 3.62 bohr and Z = 4.76 bohr for the bridge, hollow and top sites, respectively, indicative of the late barriers for the fixed sites of the title molecule-surface interaction system. Hence, the vibrational excitation of HCl is more eﬃcient in enhancing the reactivity of the reaction than the translational motion. The reaction barrier height for the bridge site is lowest (0.64 eV), compared to those for the hollow site of 0.78 eV and the top site of 0.83 eV, indicating the dissociation on the fixed bridge site is much more favored. Because the overall behavior of the contour plot for the fixed fcc site is similar to that of the hollow site despite of the slightly higher (0.02 eV) barrier height of the fcc site, we do not present the plot for the fcc site here. Figure 10 shows the dissociation probability for HCl initially in the ground rovibrational state (v = 0, j = 0) calculated using the six-dimensional (6D) TDWP method, together with the four-dimensional (4D) results for which HCl is fixed at the bridge, fcc, hollow and top cites, in the kinetic energy region of 0.4 eV and 1.6 eV. Overall, all the dissociation probabilities are smooth functions of the kinetic energy, and rise steadily with increase of the kinetic energy. The 4D dissociation probability for the fixed bridge site is much larger than the rest of three 4D probabilities and 6D probability in the entire energy region, presenting the lowest threshold (roughly 0.56 eV) compared to other results. This is quite understandable that the bare barrier height for the fixed bridge site is lowest as seen from Figure 9. The overall behavior of 4D results for the fixed fcc site is very similar to that of the hollow site, with the probability curves rise from the threshold of around 0.75 eV. The 4D dissociation probability for the hollow site is slightly larger than that for the fcc site just above the threshold resulting from the slightly lower barrier for the hollow cite as mentioned above; however, when the kinetic energy increases further, the two results become comparable in the medium energy region, and the dissociation probability for the fcc cite turns to be larger than the hollow site at the kinetic energy above 1.3 eV. It is interesting that the

Liu TH, et al.

Sci China Chem

6.0

0.8 (a)

bridge site

5.5

Dissociation probability

0.7 * TS

0.7

r (bohr)

0.6

4.0 3.5

0.6

3.0 2.5 2.0 3.5

4D bridge site 4D fcc site 4D hollow site 4D top site ×100 6D

Eb = 0.64 eV

5.0 4.5

153

January (2014) Vol. 57 No. 1

4.0

4.5

5.0 5.5 Z (bohr)

6.0

6.5

0.6

0.4

0.2

7.0 0 0.4

6.0 (b)

hollow site

5.5

Eb = 0.78 eV

5.0

Figure 10

0.8 1.2 Kinetic energy (eV)

1.6

Comparison of dissociation probabilities for the scattering of

HCl (v = 0, j = 0) from Au(111) between the six-dimensional calculation and four-dimensional fixed-site calculations.

lity is smaller than that of the bridge site in the entire energy region.

0.8

0.7 * TS 0.7

4.0

0.8

r (bohr)

4.5

3.5 3.0

4 Conclusions

2.5 2.0 3.5

4.0

4.5

5.0 5.5 Z (bohr)

6.0

6.5

7.0

6.0 (c)

top site

5.5

Eb = 0.83 eV

5.0

3.5

0.8 * TS 0.8

0.9

4.0

0.9

r (bohr)

4.5

3.0 2.5 2.0 4.0

Figure 9

4.5

5.0

5.5 Z (bohr)

6.0

6.5

7.0

(a) Contour plot of the NN1 PES as a function of Z and r for the

bridge site, with other coordinates optimized. The contours are relative to the HCl + Au(111) asymptote with an interval of 0.1 eV. The reaction barriers are also indicated by red star. (b) Same as (a) except for the hollow site. (c) Same as (a) except for the top site.

dissociation probability for the fixed top site is about 100–200 times smaller than the other results, though the bare barrier height for the top site is only 0.05 eV higher than the hollow site and 0.19 eV higher than the bridge site. Again as seem from Figure 10, the 6D probability is quite diﬀerent from the results calculated by reduced dimensional (4D) fixed-site approach, indicating the importance of full-dimensional calculations. The threshold of the 6D dissociation probability is about 0.02 eV higher than that of the 4D bridge site, in accord with the diﬀerence of the bare barrier heights of them shown in Figures 8 and 9, and the magnitude of the 6D probabi-

To summarize, we present a six-dimensional PES for the dissociative chemisorption of HCl on Au(111) based on NN fitting of more than 70000 DFT energy points, which was calculated using VASP with the moderate slab model. The resulting PES is accurate, with a very small RMSE of about 4.5 meV up to the potential energy of 4.0 eV relative to the HCl + Au(111) asymptote, and very good agreement between the PES and direct DFT calculations is achieved for the minimum energy path. Furthermore, the PES is well converged with respect to the number of data points as well as the fitting procedure, from extensive quantum dynamical calculations using the TDWP method. The four-dimensional fixedsite calculations for HCl initially in the ground rovibrational state show that the dissociation process is much more favored on the fixed bridge site, due to the lowest reaction barrier compared to hollow, fcc, and top sites. The dissociation probability for the top site is several hundred times smaller than those for the other sites, though the reaction barrier is not significantly (only 0.19 eV) higher than the bridge site. The comparison of the six-dimensional and four-dimensional fixed-site dissociation probabilities shows large diﬀerences, indicating it is very important to do full-dimensional calculations for the title molecule-surface interaction system. This is the first global PES for the dissociative chemisorption of HCl on Au(111), which is suﬃciently accurate to be used for the quantum dynamical investigations, and the following work for the detailed quantum dynamical calculations on this PES is underway and expected to be reported later. This work was supported by the National Natural Science Foundation of China (91221301 and 90921014), the Chinese Academy of Sciences, and Ministry of Science and Technology of China.

154

Liu TH, et al.

Sci China Chem

1 Beck RD, Maroni P, Papageorgopoulos DC, Dang TT, Schmid MP, Rizzo TR. Vibrational mode-specific reaction of methane on a nickel surface. S cience, 2003, 302: 98–100 2 Gates GA, Darling GR, Holloway S. A theoretical study of the vibrational excitation of NO/Ag(111). J Chem Phys, 1994, 101: 6281–6288 3 Juurlink LBF. Comparative study of C-H stretch and bend vibrations in methane activation on Ni(100) and Ni(111). Phys Rev Lett, 2005, 94: 208303(1–4) 4 Juurlink LBF, Killelea DR, Utz AL. State-resolved probes of methane dissociation dynamics. Prog Surf Sci, 2009, 84: 69–134 5 Killelea DR, Campbell VL, Shuman NS, Utz AL. Bond-selective control of a heterogeneously catalyzed reaction. S cience, 2008, 319: 790– 793 6 Kimman J, Rettner C, Auerbach D, Barker J, Tully J. Correlation between kinetic-energy transfer to rotation and to phonons in gas-surface collisions of NO with Ag(111). Phys Rev Lett, 1986, 57: 2053–2056 7 Kurten T, Biczysko M, Rajam¨aki T, Laasonen K, Halonen L. Computational study of the adsorption energetics and vibrational wavenumbers of NH3 adsorbed on the Ni(111) surface. J Phys Chem B, 2005, 109: 8954–8960 8 Lee MB, Yang QY, Ceyer ST. Dynamics of the activated dissociative chemisorptions of CH4 and implication for the pressure gap in catalysis: A molecular beam-high resolution electron energy loss study. J Chem Phys, 1987, 87: 2724–2741 9 Luntz AC. CH4 dissociation on Ni(100): Comparison of a direct dynamical model to molecular beam experiments. J Chem Phys, 1995, 102: 8264–8269 10 Luntz AC. Chemistry. Toward vibrational mode control in catalysis. Science, 2003, 302: 70–71 11 Maroni P, Papageorgopoulos DC, Sacchi M, Dang TT, Beck RD, Rizzo TR. State-resolved gas-surface reactivity of methane in the symmetric C-H stretch vibration on Ni(100). Phys Rev Lett, 2005, 94: 246104 (1-4) 12 Michelsen HA, Rettner CT, Auerbach DJ, Zare RN. Eﬀect of rotation on the translational and vibrational energy dependence of the dissociative adsorption of D2 on Cu(111). J Chem Phys, 993, 98: 8294–8307 13 Rettner CT, Michelsen HA, and Auerbach DJ. Quantum-state-specific dynamics of the issociative adsorption and associative desorption of H2 at a Cu(111) surface. J Chem Phys, 1995, 102: 4625–4641 14 Schmid MP, Maroni P, Beck RD, Rizzo TR. Surface reactivity of highly vibrationally excited molecules prepared by pulsed laser excitation: CH4 (2ν3 ) on Ni(100). J Chem Phys, 2002, 117: 8603–8606 15 Smith RR, Killelea DR, DelSesto DF, Utz AL. Preference for vibrational over translational energy in a gas-surface reaction. Science, 2004, 304: 992–995 16 Dai J, Light JC. Six dimensional quantum dynamics study for dissociative adsorption of H2 on Cu(111) surface. J Chem Phys, 1997, 107: 1676–1679 17 Dai J, Zhang JZH. Quantum adsorption dynamics of a diatomic molecule on surface: Four-dimensional fixed-site model for H2 on Cu(111). J Chem Phys, 1995, 102: 6280–6289 18 Dai J, Sheng J, Zhang JZH. Symmetry and rotational orientation eﬀects in dissociative adsorption of diatomic molecules on metals: H2 and HD on Cu(111). J Chem Phys, 1994, 101: 1555–1563 19 Gr¨uneich A, Cruz AJ, Jackson B. The dynamics of H2 dissociation on Cu and Ni surfaces. Mixed quantum-classical studies with all degrees of freedom. J Chem Phys, 1993, 98: 5800–5808 20 Kroes GJ, Wiesenekker G, Baerends EJ, Mowrey RC, Neuhauser D. Dissociative chemisorption of H2 on Cu(100): A four-dimensional study of the eﬀect of parallel translational motion on the reaction dynamics. J Chem Phys, 1996, 105: 5979–5998 21 Saalfrank P, Miller WH. Time-independent quantum dynamics for diatom-surface scattering. J Chem Phys, 1993, 98: 9040–9052

January (2014) Vol. 57 No. 1

22 Sheng J, Zhang JZH. Quantum dynamics studies of adsorption and desorption of hydrogen at a Cu(111) surface. J Chem Phys, 1993, 99: 1373–1381 23 Vincent JK, Olsen RA, Kroes GJ, Luppi M, Baerends EJ. Sixdimensional quantum dynamics of dissociative chemisorption of H2 on Ru(0001). J Chem Phys, 2005, 122: 44701(1–8) 24 Pijper E, Kroes GJ, Olsen RA, Baerends EJ. Reactive and diﬀractive scattering of H2 from Pt(111) studied using a six-dimensional wave packet method. J Chem Phys, 2002, 117: 5885 25 Pijper E, Kroes GJ, Olsen RA, Baerends EJ. Dissociative and diﬀractive scattering of H2 from Pt(111): A four-dimensional quantum dynamics study. J Chem Phys, 2002, 116: 9435–9448 26 Eichler A, Hafner J, Groß A, Scheﬄer M. Trends in the chemical reactivity of surfaces studied by ab initio quantum-dynamics calculations. Phys Rev B, 1999, 59: 13297-13300 27 Krishnamohan GP, Olsen RA, Kroes GJ, Gatti F, Woittequand S. Quantum dynamics of dissociative chemisorption of CH4 on Ni(111): Influence of the bending vibration. J Chem Phys, 2010, 133: 144308(1–15) 28 Nave S, Jackson B. Methane dissociation on Ni(111) and Pt(111): Energetic and dynamical studies. J Chem Phys, 2009, 130: 054701(1–14) 29 Nave S, Tiwari AK, Jackson B. Methane dissociation and adsorption on Ni(111), Pt(111), Ni(100), Pt(100), and Pt(110)-(1×2): Energetic study. J Chem Phys, 2010, 132: 054705(1–12) 30 Jackson B, Nave S. The dissociative chemisorption of methane on Ni(100): Reaction path description of mode-selective chemistry. J Chem Phys, 2011, 135: 114701(1–12) 31 Jiang B, Liu R, Li J, Xie D, Yang M, Guo H. Mode selectivity in methane dissociative chemisorption on Ni(111). Chem Sci, 2013, 4: 3249–3254 32 Jiang B, Xie D, Guo H. Vibrationally mediated bond selective dissociative chemisorptions of HOD on Cu(111). Chem Sci, 2013, 4: 503–508 33 Jiang B, Li J, Xie D, Guo H. Eﬀects of reactant internal excitation and orientation on dissociative chemisorption of H2 O on Cu(111): Quasiseven-dimensional quantum dynamics on a refined potential energy surface. J Chem Phys, 2013, 138: 044704(1–9) 34 Jiang B, Ren X, Xie D, Guo H. Enhancing dissociative chemisorption of H2 O on Cu(111) via vibrational excitation. Proc Natl Acad Sci USA, 2012, 109: 10224–10227 35 Bowman JM, Czak´o G, Fu B. High-dimensional ab initio potential energy surfaces for reaction dynamics calculations. Phys Chem Phys Chem, 2011, 13: 8094–8111 36 Braams BJ, Bowman JM. Permutationally invariant potential energy surfaces in high dimensionality. Int Rev Phys Chem, 2009, 28: 577– 606 37 Lykke KR, Kay BD. Rotationally inelastic gas-surface scattering: HCl from Au(111). J Chem Phys, 1990, 92: 2614–2623 38 Cooper R, Rahinov I, Yuan C, Yang X, Auerbach DJ, Wodtke AM. Efficient translational excitation of a solid metal surface: State-to-state translational energy distributions of vibrational ground state HCl scattering from Au(111). J Vac Sci Technol A, 2009, 27: 907–912 39 Rahinov I, Cooper R, Yuan C, Yang X, Auerbach DJ, Wodtke AM. Eﬃcient vibrational and translational excitations of a solid metal surface: State-to-state time-of-flight measurements of HCl(v = 2, J = 1) scattering from Au(111). J Chem Phys, 2008, 129: 214708(1–16) 40 Ran Q, Matsiev D, Auerbach DJ, Wodtke AM. Observation of a change of vibrational excitation mechanism with surface temperature: HCl collisions with Au(111). Phys Rev Lett, 2007, 98: 237601(1–4) 41 Ran Q, Matsiev D, Auerbach DJ, Wodtke AM. Direct translation-tovibrational energy transfer of HCl on gold: Measurement of absolute vibrational excitation probabilities. Nucl Instrum Meth B, 2007, 258: 1–6 42 Behler J, Parrinello M. Generalized neural-network representation of high-dimensional potential-energy surfaces. Phys Rev Lett, 2007, 98:

Liu TH, et al.

Sci China Chem

146401(1–4) 43 Blank TB, Brown SD, Calhoun AW, Doren DJ. Neural network models of potential energy surfaces. J Chem Phys, 1995, 103: 4129–4137 44 Brown DFR, Gibbs MN, Clary DC. Combining ab initio computations, neural networks, and diﬀusion Monte Carlo: An eﬃcient method to treat weakly bound molecules. J Chem Phys, 1996, 105: 7597–7604 45 Darley MG, Handley CM, Popelier PLA. Beyond point charges: Dynamic polarization from neural net predicted multipole moments. J Chem T heory Comput, 2008, 4: 1435–1448 46 Jose KV, Artrith N, Behler J. Construction of high-dimensional neural network potentials using environment-dependent atom pairs. J Chem Phys, 2012, 136: 194111(1–15) 47 Le HM, Huynh S, Raﬀ LM. Molecular dissociation of hydrogen peroxide (HOOH) on a neural network ab initio potential surface with a new configuration sampling method involving gradient fitting. J Chem Phys, 2009, 131: 014107(1–10) 48 Le HM, Raﬀ LM. Cis→trans, trans→cis isomerizations and N-O bond dissociation of nitrous acid (HONO) on an ab initio potential surface obtained by novelty sampling and feed-forward neural network fitting. J Chem Phys, 2008, 128: 194310(1–11) 49 Le HM, Raﬀ LM. Molecular dynamics investigation of the bimolecular reaction BeH + H2 →BeH2 + H on an ab initio potential-energy surface obtained using neural network methods with both potential and gradient accuracy determination. J Phys Chem A, 2009, 114: 45–53 50 Malshe M, Raﬀ LM, Rockley MG, Hagan M, Agrawal PM, Komanduri R. Theoretical investigation of the dissociation dynamics of vibrationally excited vinyl bromide on an ab initio potential-energy surface obtained using modified novelty sampling and feedforward neural networks. II. Numerical application of the method. J Chem Phys, 2007, 127: 134105(1–7) 51 Manzhos S, Wang X, Dawes R, Carrington T. A nested moleculeindependent neural network approach for high-quality potential fits. J Phys Chem A, 2005, 110: 5295–5304 52 Pukrittayakamee A, Malshe M, Hagan M, Raﬀ LM, Narulkar R, Bukkapatnum S, Komanduri R. Simultaneous fitting of a potentialenergy surface and its corresponding force fields using feedforward neural networks. J Chem Phys, 2009, 130: 134101(1–10) 53 Raﬀ LM, Malshe M, Hagan M, Doughan DI, Rockley MG, Komanduri R. Ab initio potential-energy surfaces for complex, multichannel systems using modified novelty sampling and feedforward neural networks. J Chem Phys, 2005, 122: 84104(1–16) 54 Behler J, Delley B, Lorenz S, Reuter K, Scheﬄer M. Dissociation of O2 at Al(111): The role of spin selection rules. Phys Rev Lett, 2005, 94: 036104(1–4) 55 Lorenz S, Groß A, Scheﬄer M. Representing high-dimensional potential-energy surfaces for reactions at surfaces by neural networks. Chem Phys Lett, 2004, 395: 210–215

January (2014) Vol. 57 No. 1

155

56 Lorenz S, Scheﬄer M, Gross A. Descriptions of surface chemical reactions using a neural network representation of the potential-energy surface. Phys Rev B, 2006, 73: 115431(1–13) 57 Behler J, Lorenz S, Reuter K. Representing molecule-surface interactions with symmetry-adapted neural networks. J Chem Phys, 2007, 127: 014705(1–11) 58 Behler J, Reuter K, Scheﬄer M. Nonadiabatic eﬀects in the dissociation of oxygen molecules at the Al(111) surface. Phys Rev B, 2008, 77: 115421(1–16) 59 Chen J, Xu X, Xu X, Zhang DH. Communication: An accurate global potential energy surface for the OH + CO → H + CO2 reaction using neural networks. J Chem Phys, 2013, 138: 221104(1–4) 60 Chen J, Xu X, Xu X, Zhang DH. A global potential energy surface for the H2 + OH ↔ H2 O + H reaction using neural networks. J Chem Phys, 2013, 138: 154301(1–8) 61 Kresse G, Furthm¨uller J. Eﬃciency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. Comp Mater Sci, 1996, 6: 15–50 62 Kresse G and Furthm¨uller J. Eﬃcient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys Rev B, 1996, 54: 11169–11186 63 Bl¨ochl PE. Projector augmented-wave method. Phys Rev B, 1994, 50: 17953–17979 64 Kresse G, Joubert D. From ultrasoft pseudopotentials to the projector augmented-wave method. Phys Rev B, 1999, 59: 1758–1775 65 Perdew JP, Burke K, Ernzerhof M. Generalized gradient approximation made simple. Phys Rev Lett, 1996, 77: 3865–3868 66 Perdew JP, Chevary JA, Vosko SH, Jackson KA, Pederson MR, Singh DJ, Fiolhais C. Atoms, molecules, solids, and surfaces: Applications of the generalized gradient approximation for exchange and correlation. Phys Rev B, 1992, 46: 6671–6687 67 Monkhorst HJ, Pack JD. Special points for Brillouin-zone integrations. Phys Rev B, 1976, 13: 5188–5192 68 Haynes WM. CRC Handbook o f Chemistry and Physics. 93rd edition (Internet version 2013). Boca Raton, Fl: CRC Press/Taylor and Francis 69 Henkelman G, Uberuaga BP, J´onsson H. A climbing image nudged elastic band method for finding saddle points and minimum energy paths. J Chem Phys, 2000, 113: 9901–9904 70 Jonsson H, Mills G, Jacobsen KW. Nudged elastic band method f or f inding minimum energy paths o f transitions. Classical and Quantum Dynamics in Condensed Phase S imulations. Eds. Berne BJ, Ciccotti G, Coker DF. Singapore: World Scientific, 1998 71 Hagan MT, Menhaj MB. Training feedforward networks with the Marquardt algorithm. IEEE T rans Neural Netw, 1994, 5: 989–993 72 Behler J. Neural network potential-energy surfaces in chemistry: A tool for large-scale simulations. Phys Chem Chem Phys, 2011, 13: 17930– 17955