Chinese Journal of Chemical Physics  2017, Vol. 30 Issue (1): 71-76

The article information

Yue-mei Sun, Xiang-jian Shen, Xiao-hong Yan
孙月梅, 沈祥建, 颜晓红
Molecular Dynamics Study of Hydrogen Dissociation on Pd Surfaces using Reactive Force Fields
基于反应力场的分子动力学方法研究氢在钯表面的分解
Chinese Journal of Chemical Physics, 2017, 30(1): 71-76
化学物理学报, 2017, 30(1): 71-76
http://dx.doi.org/10.1063/1674-0068/30/cjcp1605096

Article history

Received on: May 6, 2016
Accepted on: November 16, 2016
Molecular Dynamics Study of Hydrogen Dissociation on Pd Surfaces using Reactive Force Fields
Yue-mei Suna, Xiang-jian Shenb, Xiao-hong Yana     
Dated: Received on May 6, 2016; Accepted on November 16, 2016
a. College of Science, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China;
b. Research Center of Heterogeneous Catalysis and Engineering Sciences, School of Chemical Engineering and Energy, Zhengzhou University, Zhengzhou 450001, China
*Author to whom correspondence should be addressed. Yue-mei Sun,E-mail:sunyuemei@nuaa.edu.cn
Abstract: Developing a widely-used reactive force field is meaningful to explore the fundamental reaction mechanism on gas-surface chemical reaction dynamics due to its very high computational efficiency. We here present a study of hydrogen and its deuterated molecules dissociation on Pd surfaces based on a full-dimensional potential energy surface (PES) constructed by using a simple second moment approximation reactive force field (SMA RFF). Although the descriptions of the adsorbate-substrate interaction contain only the dissociation reaction of H2/Pd(111) system, a good transferability of SMA potential energy surface (PES) is shown to investigate the hydrogen dissociation on Pd(100). Our simulation results show that, the dissociation probabilities of H2 and its deuterated molecules on Pd(111) and Pd(100) surfaces keep non-monotonous variations with respect to the incident energy Ei, which is in good agreement with the previous ab initio molecular dynamics. Furthermore, for the oriented molecules, the dissociation probabilities of the oriented H2 (D2 and T2) molecule have the same orientation dependence behavior as those oriented HD (HT and DT) molecules.
Key words: Surface reaction dynamics    Hydrogen dissociation    Reactive force fields    Isotope effect    
I. INTRODUCTION

Molecule-surface interactions have long attracted much attention due to a lot of important surface chemical reactions in heterogeneous catalysis [1-3]. The dissociative adsorption of hydrogen (H2) molecule on metal and alloy surfaces is one benchmark system [4-5]. With a large number of supersonic molecular beams and state-resolved experiments [6-15], much attention to H2 dissociation on bare Pd surfaces has historically been paid because of this typical non-activated gas-surface reaction. From the view of dynamics, many important dynamical studies [16-22] using classical molecular dynamics, ab initio molecular dynamics (AIMD) and quantum dynamics approaches have been devoted to the fundamental understanding of novel dissociation mechanisms and characters, e.g., the dissociation probability (S0) as a function of the initial incident energy presents a non-monotonous dependence due to the significant trapping molecular 'precursor'; H2 initially rotating in cartwheel-like states has the lower S0 than that in helicopter-like states. Recently, much attention is paid to the partially deuterated HD molecules since it is an ideal molecule for experiments on rotational excitation. There have been some diffractive, rotationally elastic and inelastic scattering experiments on HD molecules interacting with Pt [6], Ni [8], Cu [13], and Au [7] surfaces as well as Pd [15]. The large difference between various metal surfaces was observed in the amount of rotational excitation. One significant factor is attributed to the fact that molecular center of mass of the HD molecule is displaced from the corresponding molecular center of charge which effectively affects angular momentum during collision process. Kingma and coworker [23] first presented a quantum dynamical study of dissociative adsorption of HD and H2 molecules interacting with an activated Pt(111) and reported slightly smaller dissociative adsorption probabilities of HD molecules. Another dynamical study of the HD and H2 molecules interacting with Cu(110) and Pd(111) surfaces was done by Ramirez et al. [24] by classical trajectory method based on a six dimensional potential energy surface (PES) using the corrugation reducing procedure (CRP) [25-26]. Very recently, using stimulated Raman pump and linearly polarized laser, the well-defined, non-statistical aligned or oriented HD and H2 molecules can be experimentally prepared in gas phase [14, 27]. One interesting alignment effect of H2 interacting with Pd(111) has been experimentally investigated by Sitz group [14]. However, few theoretical attempts are made by straightforwardly emphasizing on molecular orientation dependence of the partially deuterated HD molecule interacting with a non-activated metal surface.

In this work, we present the good transferability of a simple second moment approximation reactive force field through some simulations on hydrogen and its deuterated molecules dissociation on Pd(111) and Pd(100) surfaces.

II. METHOD A. Second moment approximation reactive force field (SMA RFF)

The interaction potential between the incident molecules and the substrate (Pd(111) and Pd(100) surfaces) is described by a full-dimensional potential energy surface (PES) based on SMA RFF. We use exactly the same expressions for the SMA RFF as those in Ref.[28]. In order to determine the parameters in the RFF, we use a database containing only results obtained from DFT calculations with the help of the Vienna ab initio simulation package. The DFT calculations are performed with the generalized gradient approximation of Perdew and Wang (PW91). Plane-waves are used for expanding the wave functions with a cutoff energy equal to 200 eV. We use ultra soft pseudo potentials for valence electrons and a 4×4×1 k-points grid for the k-space sampling. The super cell approach is adopted which includes a slab of 5 Pd layers with a (3×3) Pd(111) surface cell and a vacuum space corresponding to five Pd layers. We take the energy of H2 at a point far from the surface (Z=7.0 Å) as the zero of the potential energy of the system. DFT calculated configurations are exactly the same as those in Ref.[28].

We use the widely spread least square scheme to carry out the fitting by minimizing the chi-square function defined in Ref.[28]. The fitting is carried out with only data related to the total energy. All the data in the database do not have the same importance. Therefore, there is a weight to each data according to its importance. In Ref.[28], the weight coefficient is determined by using the following expression for SMA RFF,

(1)

In Ref.[29], Shen et al. used the following formula to determine the weight coefficient instead of Eq.(1):

(2)

The results showed that SMA RFF parameterized with Eq.(2) gives a much more reliable PES than SMA RFF with Eq.(1). According to Ref.[29], the SMA RFF parameterized with Eq.(2) leading to qualitatively correct results for the dynamics (i.e., nonmonotonous variation of S0 with respect to Ei for H2/Pd(111)). Moreover, the quantitative difference between various results (e.g. CRP-MD and AIMD) for S0 being of the order of 30% (the same order as that found for the experimental results obtained by different groups). The results indicated that the cutoff energy used in the chi-square function for fitting has an important influence on the accuracy of the parameterized RFF. The details about the accuracy and fitting procedure of such SMA(PES) can refer to our previous work [28-30].

In our present work, the constructed PES using SMARFF is the same as those performed by Shen and coworkers in Ref.[29]. Here we simply mention that this developed SMA(PES) only contains the information about H2/Pd(111) system [29, 30]. In Fig. 1, the variation of the total energy with the molecule-surface distance, Z, is shown for fcc-fcc pathway on clean Pd(111) (Fig. 1(a)) and hollow-top-hollow pathway on clean Pd(100) (Fig. 1(b)) with the bond of H2 being fixed to 0.75 Å. We see that the curve given by our SMA(PES) is qualitatively similar to that given by DFT calculation for both surfaces. It is quite remarkable that the SMA RFF constructed with a database containing only information on H2/Pd(111) allows for obtaining accurate results for H2/Pd(100). The main reason for the success of transference is due to the similarity of their electronic structures.

FIG. 1 Variation of the potential energy with respect to the distance of H2 to the surface with the bond of H2 being fixed to 0.75 Å and parallel to the surface along fcc-fcc pathway on clean Pd(111) (a) and hollow-top-hollow pathway on clean Pd(100). (b) DFT results (solid squares) and SMA RFF results (open circles).
B. Dynamics calculations

Our simulations on the dissociative chemisorption of various H2 and its deuterated molecules (such as D2, T2, HD, HT, and DT) with different incident energies on Pd(111) and Pd(100) surfaces have been carried out by using classical molecular dynamics method [31]. Verleta lgorithm is used for integrating the equation of motion with a time step of 0.3 fs, which ensures a conservation of total energy within 0.5 meV in a micro canonical ensemble NVE during a period of 5.0 ps. The starting height of incident molecules is 6.0 Å above the surface and only normal incidence is considered in the present work. The normal incident energy varies from 10 meV to 400 meV for each molecule. The initial impact position of incident molecules is chosen randomly. For a random molecule, the initial orientation is chosen randomly, while for an oriented molecule, it is prepared with a given orientation angle θ, i.e., the angle between molecular bond and the surface normal. The given orientation angle varies from 0° to 180° with an interval angle of 30°. We performed only classical trajectory simulations, i.e., without the zero point vibrational energy being included in the initial state of all the molecules. We consider the molecule as dissociated once its bond length is larger than 2.25 Å and two H (D or T) atoms have the opposite velocities along the bond. A molecule is considered as reflected when it returns to the starting height with a velocity pointing to the gas phase. The dissociation probability (S0) is the ratio of the number of reactive trajectories (Ndis) to the number of total trajectories (Ntot), i.e., S0=Ndis/Ntot. It is calculated with 1000 trajectories and the statistical errors are less than 0.02. All the dynamics simulations are carried out on non-rigid surfaces, i.e., simulations involving the motion of some substrate atoms, three bottom layers of the slab are fixed while the uppermost two layers are allowed to move.

III. RESULTS AND DISCUSSION A. Transfer ability of SMA(PES)

First of all, we demonstrate the S0 of H2 dissociation on a non-rigid Pd(111) and a non-rigid Pd(100) surface with the previous theoretical and experimental studies [9, 11, 28, 32]. In Fig. 2, it is clear to note that the S0 obtained from our SMA RFF has the similar tendency with these theoretical results such as AIMD methods. Importantly, these results indicate that our SMA RFF have a good transferability for investigating hydrogen dissociation on Pd(100) surface, although the descriptions of the adsorbate-substrate interaction contain only the dissociation reaction of H2/Pd(111) system.

FIG. 2 Dissociation probability S0 of random H2 molecule as a function of incident energy Ei calculated from REBO (solid squares) [28], SMA RFF (solid circles), AIMD (solid triangles) [32] and from experiment (open squares and circles)[9, 11] (a) on a non-rigid Pd(111) and (b) on a non-rigid Pd(100).
B. S0 of random H2, D2, T2, HD, HT, and DT molecules

In order to understand molecular orientation dependence behavior on the dissociation probability, we carry out some various classical molecular dynamics simulations on random H2, D2, T2, HD, HT, and DT molecules dissociation on Pd(111) and Pd(100) surfaces respectively. For the non-activated hydrogen/Pd gas-surface reactions, hydrogen molecule on clean Pd surfaces with low incident energies combines the direct dissociation mechanism with the dynamic trapping dissociation mechanism. Consequently, it yields the minimum of S0 at low incident energy. Figure 3 presents the final S0 of the random H2, D2, T2, HD, HT, and DT molecules dissociation on non-rigid Pd(111) and non-rigid Pd(100) surfaces as a function of incident energy Ei. Clearly, all the sticking curves of these random molecules have the similar dissociation behavior, i.e., the S0 keeps the non-monotonous variations with respect to incident energy Ei. On the non-rigid Pd (111) substrate, in the region of high incident energy, there are negligible discrepancies in the S0 between the random H2, D2, T2, HD, HT and DT molecules, while some apparent differences appear in the region of low incident energy. The random D2 and HD molecules have lower minimum S0, than that random H2 molecule. And the random T2, HT, and DT molecules have the lowest S0. Such small differences might be attributed to the mass effect since no zero-point vibrational energy difference is given between the random H2, D2, T2, HD, HT, and DT molecules in our present simulations. Moreover, the dissociation probability of the HD molecule we calculated is given as about 0.7 at incident energy Ei=74 meV, which is consistent with that experimental measurement from the clean Pd(111) surface [15]. The observed S0 differences between the random H2, D2, and HD molecules on the non-rigid Pd(111) surface are in good agreement with those obtained by Ramirez et al. [24]. Similarly, on the non-rigid Pd(100) substrate, the S0 of these molecules are very close in the whole region of incident energy. In this section, we discussed the negligible isotope effect on the S0 of the random H2, D2, T2, HD, HT, and DT molecules, and next to study the dissociation probabilities of the initially oriented H2, D2, T2, HD, HT, and DT molecules and address the important molecular orientation dependence behavior.

FIG. 3 Dissociation probability S0 for random H2, D2, T2, HD, HT, and DT molecule as a function of incident energy Ei (a) on a non-rigid Pd(111) surface and (b) on a non-rigid Pd(100) surface.
C. S0 of oriented H2, D2 , T2, HD, HT, and DT molecules

As illustrated above, the oriented H2, D2, T2, HD, HT, and DT molecules are generated with the given initial orientation angle θ. In our simulations, six given orientation angles for these oriented molecules are under consideration, which are θ=30°, 60°, 90°, 120°, 150° and 180°. Figures 4 and 5 respectively present the final S0 of six different initially oriented H2, D2, T2, HD, HT, and DT molecules dissociation on the non-rigid Pd(100) and the non-rigid Pd(111) surfaces as a function of incident energy Ei. For all these oriented H2, D2, T2, HD, HT and DT molecules, their sticking curves still keep the non-monotonous behavior with respect to incident energy Ei. Interestingly, the oriented H2 (D2 and T2) molecules with a given θ (shown in Fig. 4 on Pd(100) surface and Fig. 5 on Pd(111) surface) have the almost same S0 as that of the given orientation angle (180°-θ). It is attributed to the homonuclear isotopmer H2 (D2 and T2) without distinguishing. Moreover, these non-monotonous sticking curves of the oriented H2 (D2 and T2) molecules yield different low minimum S0 values in the region of very low incident energy. The S0 of the oriented H2 (D2 and T2) with the given θ=180° is apparently much lower than that of other given orientation angles θ=30°, 60° and 90°. The oriented H2 (D2 and T2) molecules with the orientation angle of 90° has the highest S0 and can be totally dissociated from the incident energy of 100 meV on Pd(111) and 50 meV on Pd(100). Such 100% dissociation behavior indicates that the most favorable molecular orientation for H2 (D2 and T2) dissociation is H-H axis parallel to the metal surface. It is in good agreement with the previous most of dynamics studies [17]. Our dynamical simulation results provide the evidence for understanding diatomic H2 (D2 and T2) molecule dissociation mechanism. Furthermore, the S0 of the oriented H2 (D2 and T2) molecules increases with respect to the orientation angle θ from 30° to 90°, and then decreases with respect to the orientation angle θ from 90° to 180°.

FIG. 4 Dissociation probability S0 for the initially oriented molecules dissociation on a non-rigid Pd(100) surface with different given orientation angles θ (30°, 60°, 90°, 120°, 150°, and 180°) as a function of incident energy Ei. (a) H2, (b) HD, (c) D2, (d) HT, (e) T2, (f) DT.
FIG. 5 Dissociation probability S0 for the initially oriented molecules dissociation on a non-rigid Pd(111) surface with different given orientation angles θ (30°, 60°, 90°, 120°, 150°, and 180°) as a function of incident energy Ei. (a) H2, (b) HD, (c) D2, (d) HT, (e) T2, (f) DT.

However, for the initially oriented heteronuclear HD (HT and DT molecules) shown in Fig. 4 (b) (d), (f) on Pd(100) surface and Fig. 5(b), (d), (f), on Pd(111) surface, it is obviously distinct from those oriented homonuclear molecules. It is to note that the dissociation probability of the oriented heteronuclear molecule with a given orientation angle θ does not have the similar tendency to that of the orientation angle (180°-θ), although the S0 of the oriented molecule varies with respect to the increase of orientation angle θ and has the lowest dissociation probability with a given θ=180°. Moreover, for the initially oriented HD, HT, and DT on Pd(100) surface (see Fig. 4 (b), (d), and (f)), the dissociation probability firstly increases with the given θ from 30° to 120° and then decreases with the given θ from 120° to 180°. Such interesting variations in the dissociation behavior indicate the fact that the effective favorable molecular orientation for the oriented HD, HT, and DT molecules dissociation on Pd(100) surface does not keep the H-D axis parallel to the metal surface. The highest dissociation probability S0 we obtained is for those near the given θ=120°, instead of the given θ=90° shown for the homonuclear oriented molecules. This indicates that the effective dissociation orientation for molecule with asymmetric mass distribution is distinct from that with symmetric mass distribution in the dynamics simulations.

IV. CONCLUSION

In summary, we have carried out some various classical molecular dynamics simulations on H2 molecules dissociation on Pd(111) and Pd(100) surfaces as well as the D2, T2, HD, HT, and DT molecules based on the SMA(PES) we developed. Our simulation results show that, the SMA(PES) has a good transferability for investigating hydrogen dissociation on Pd(100) surface, although the descriptions of the adsorbate-substrate interaction contain only the dissociation reaction of H2/Pd(111) system. The dissociation probabilities S0 of the random and oriented H2, D2, T2, HD, HT and DT molecules keep the non-monotonous behavior with respect to the incident energy Ei. It is in good agreement with those previous dynamics studies. The initially oriented homonuclear H2, D2, and T2 molecules with a given θ has the almost same dissociation probability as that the orientation angle of a given (180°-θ). While the initially oriented heteronuclear HD, HT, and DT molecules with a given θ has different dissociation probability from that orientation angle given (180°-θ). Furthermore, the effective favorable orientation angle of the oriented H2, D2 and T2 on both Pd(100) and Pd(111) is about 90°, i.e., that H-H axis parallel to the metal surface. It shifts toward the higher orientation angle for the oriented heteronuclear HD, HT and DT molecules dissociation on Pd(111) and Pd(100) surfaces. With a good transferability of reactive force field and the important dependence of molecular orientation observed in the present work, it opens very attractive perspectives for investigating various aspects of complex polyatomic dissociation on metal surfaces and more generally provides an accurate and tractable strategy for studying reaction dynamics of more and more complex surface reactions.

V. ACKNOWLEDGMENTS

This work was supported by the National Natural Science Foundation of China (No.21506053) and Doctoral Scientific Research Foundation Project (KYY15023).

Reference
[1] H. A. Michelsen, C. T. Rettner, and D. J. Auerbach, Surface Reactions, Springer Series in Surface Science Vol. 34, Berlin:Springer-Verlag, 185(1994).
[2] Jiang and H. Guo B., Phys. Rev. Lett. 114 , 166101 (2015). DOI:10.1103/PhysRevLett.114.166101
[3] X. J. Shen, A. Lozano, W. Dong, H. F. Busnengo, and X. H. Yan, Phys. Rev. Lett. 112 , 046101 (2014). DOI:10.1103/PhysRevLett.112.046101
[4] R. Darling and S. Holloway G., Rep. Prog. Phys. 58 , 1595 (1995). DOI:10.1088/0034-4885/58/12/001
[5] G. J. Kroes, A. Gross, E. J. Baerends, M. Scheffler, and A. McCormack D., Acc. Chem. Res. 35 , 193 (2002). DOI:10.1021/ar010104u
[6] J. P. Cowin, C. F. Yu, and L. Wharton, Surf. Sci. 161 , 221 (1985). DOI:10.1016/0039-6028(85)90738-1
[7] U. Harten, J. P. Toennies, and C. Woell, J. Chem. Phys. 85 , 2249 (1986). DOI:10.1063/1.451121
[8] R. Berndt, J. P. Toennies, and C. Woell, J. Chem. Phys. 92 , 1468 (1990). DOI:10.1063/1.458105
[9] M. Beutl, M. Riedler, and K. D. Rendulic, Chem. Phys. Lett. 247 , 249 (1995). DOI:10.1016/0009-2614(95)01208-5
[10] M. Beutl, K. D. Rendulic, and G. R. Castro, J. Chem. Soc. Faraday Trans. 91 , 3639 (1995). DOI:10.1039/ft9959103639
[11] Gostein and G. O. Sitz M., J. Chem. Phys. 106 , 7378 (1997). DOI:10.1063/1.473699
[12] H. Hou, S. J. Gulding, C. T. Rettner, A. M. Wodtke, and D. J. Auerbach, Science 277 , 80 (1997). DOI:10.1126/science.277.5322.80
[13] L. V. Goncharova, J. Braun, A. V. Ermakov, G. G. Bishop, D. M. Smilgies, and B. J. Hinch, J. Chem. Phys. 115 , 7713 (2001). DOI:10.1063/1.1403001
[14] M. Isakson, PhD Dissertation, Austin:The University of Texas, (2002).
[15] C. Shackman and G. O. Sitz L., J. Chem. Phys. 122 , 114702 (2005). DOI:10.1063/1.1861884
[16] H. F. Busnengo, W. Dong, P. Sautet, and A. Salin, Phys. Rev. Lett. 87 , 127601 (2001). DOI:10.1103/PhysRevLett.87.127601
[17] H. F. Busnengo, E. Pijper, G. J. Kroes, and A. Salin, J. Chem. Phys. 119 , 12553 (2003). DOI:10.1063/1.1626535
[18] C. Diaz, H. F. Busnengo, F. Martin, and A. Salin, J. Chem. Phys. 118 , 2886 (2003). DOI:10.1063/1.1524160
[19] A. Di Cesare M., H. F. Busnengo, W. Dong, and A. Salin, J. Chem. Phys. 118 , 11226 (2003). DOI:10.1063/1.1575208
[20] H. F. Busnengo, W. Dong, and A. Salin, Phys. Rev. Lett. 93 , 236103 (2004). DOI:10.1103/PhysRevLett.93.236103
[21] Dianat and A. Gross A., J. Chem. Phys. 120 , 5339 (2004). DOI:10.1063/1.1647519
[22] C. Diaz, M. F. Somers, G. J. Kroes, H. F. Busnengo, A. Salin, and F. Martin, Phys. Rev B72 , 035401 (2005).
[23] S. M. Kingma, M. F. Somers, E. Pijper, G. J. Kroes, R. A. Olsen, and E. J. Baerends, J. Chem. Phys. 118 , 4190 (2003). DOI:10.1063/1.1540981
[24] Ramirez and H. Busnengo C., Surf. Sci. 603 , 3171 (2009). DOI:10.1016/j.susc.2009.07.017
[25] H. F. Busnengo, C. Crespos, W. Dong, J. C. Rayez, and A. Salin, J. Chem. Phys. 116 , 9005 (2002). DOI:10.1063/1.1471248
[26] A. Salin, J. Chem. Phys. 124 , 104704 (2006). DOI:10.1063/1.2178357
[27] N. C. Bartlett, D. J. Miller, R. N. Zare, D. Sofikitis, T. P. Rakitzis, and A. J. Alexander, J. Chem. Phys. 129 , 084312 (2008). DOI:10.1063/1.2973628
[28] Y. Xiao, W. Dong, and H. Busnengo, J. Chem. Phys. 132 , 014704 (2010). DOI:10.1063/1.3265854
[29] X. J. Shen, W. Dong, Y. Xiao, and X. H. Yan, J. Chem. Phys. 135 , 167101 (2011). DOI:10.1063/1.3652026
[30] X. J. Shen, Y. Xiao, W. Dong, X. H. Yan, and H. Busnengo, Comput. Theor. Chem. 990 , 152 (2012). DOI:10.1016/j.comptc.2012.03.012
[31] M. Allen and D. Tildsley, Computer Simulation of liquids, Oxford:Clarendon Press, (1991).
[32] Groß and A. Dianat A., Phys. Rev. Lett. 98 , 206107 (2007). DOI:10.1103/PhysRevLett.98.206107
基于反应力场的分子动力学方法研究氢在钯表面的分解
孙月梅a, 沈祥建b, 颜晓红a     
a. 南京航空航天大学理学院, 南京 210016;
b. 郑州大学化工与能源学院, 多相催化与工程科学研究中心, 郑州 450001
摘要: 基于二阶矩近似反应力场方法构建的全维度势能面研究了氢分子及其同位素分子在钯表面的分解过程.在构建势能面的过程中数据库中只包含了氢分子与钯(111)表面相互作用的相关信息,该势能面在研究氢分子在钯(100)表面上的分解过程中表现出了非常好的可转移性.结果表明,氢分子及其同位素分子在钯(111)与钯(100)表面上的分解系数S0均随着入射能量的增加呈现非单调变化,并且通过固定分子取向的方法发现同核分子(H2、D2和T2)最有利分解取向角为90°,而异核分子(HD、HT和DT)受质心偏移的影响,其最有利分解取向角向大角度偏移.
关键词: 表面反应动力学    氢分解    反应力场    同位素效应