A Molecular Dynamics Study for the Dissociation Phenomena of Gas Molecule on Metal Surface ∗

Dissociation phenomena of gas molecule on a metal surface, especially, the eﬀect of motion of a gas molecule on dissociation probability of the molecule on a metal surface was analyzed by Molecular Dynamics (MD) method. Platinum (111) surface and hydrogen were chosen as the metal surface and the gas molecule, respectively. Embedded Atom Method (EAM) was used to express the interaction potential as the functional of the electron density. Parameters or functions of the EAM potential were determined so that the characteristics of the interaction potential obtained by the EAM method were consistent with those obtained by Density Functional Theory (DFT). Dissociation probabilities at speciﬁc sites (top, brg and fcc) were obtained by MD method against impinging energy. On the other hand, the dissociation probabilities were determined from dissociation barriers of the sites and orientations of hydrogen molecule. These results were compared with each other and the eﬀect of motion of atoms or molecules on dissociation probability was analyzed. [DOI: 10.1380/ejssnt.2010.211]


I. INTRODUCTION
The mechanism of dissociative adsorption phenomena of a gas molecule on a metal surface becomes very important and this problem has been investigated using various techniques. Recently, molecular simulation has become a more effective technique. At present, analyses of dissociation path or dissociation barrier are mainly performed by Density Functional Theory (DFT) [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19]. Olsen et al. obtained the Potential Energy Surface (PES) of various orientations of H 2 molecule on various transition metal surfaces by the DFT, and constructed a potential in which the total energy of the system is expressed by analytical functions of relative position, orientation, and the bond length of the molecule by interpolating the DFT data [12,13]. In DFT, however, the temperature of the system is assumed to be 0 K, that is, all the atoms of the surface are fixed, and therefore the effect of motion of the atoms of the surface on the dissociation phenomena cannot be considered. It has been reported that the dissociation probability obtained by DFT differs from that obtained experimentally [20]. It is necessary to analyze the phenomena by Molecular Dynamics (MD) method in order to consider the effect of the motion of the atoms on the dissociation phenomena. Embedded Atom Method (EAM) is a method to consider the interaction between gas and metal surface [21][22][23][24][25][26]. This method is based on the DFT, and the total energy of a system is expressed * This paper was presented at 10th International Conference on Atomically Controlled Surfaces, Interfaces and Nanostructures (ACSIN-10), Granada Conference Centre, Spain, 21-25 September, 2009. † Corresponding author: tokumasu@ifs.tohoku.ac.jp as the sum of the energy required to embed an atom into the background electron density and that of the two-body interaction potential between atoms. The EAM method is often used to analyze the motion of surface atoms on a transition metal surface and has been applied to an fcc metal surface or its alloys [26]. However, there have been little analyses of the dissociation phenomena with moving atoms by the EAM, although Xie et al. used the EAM method to obtain the PES to calculate the dissociation probability of Cu(111)/H 2 with a lattice detect [27]. Analysis of the dissociation phenomena of gas molecule by the EAM method considering the motion of atoms can obtain more realistic information about the phenomena.
We have analyzed the effect of motion of surface atoms on dissociation phenomena by MD method [28] and concluded that the transition region of dissociation probability is wider when the surface atoms move. However, the parameters of the EAM potential were obtained based on the data of only the brg site and therefore the EAM potential does not have the accuracy enough to analyze the difference of the characteristics of dissociation phenomena between each specific site (top brg and fcc). In the present paper, the parameter of the EAM potential was reconstructed to reproduce the characteristics of dissociation barrier at specific sites, and the difference of the characteristics of dissociation phenomena at the sites was analyzed. A Pt(111) surface was chosen as the metal surface, and hydrogen was chosen as the gas molecule of the system. Simulations in which an H 2 molecule impinges to a Pt(111) surface were performed several times by changing the orientation of the impinging H 2 molecule and dissociation probability at specific site (top, brg and fcc) of the surface was obtained. From the results, the effect of the motion of the atoms, especially gas molecules, on the dissociation probability at specific sites was analyzed.

A. Shapes of the functions
In EAM, the potential energy, E pot , of a system, which consists N Pt atoms and an H 2 molecule, is given by where subscripts i and j denote Pt atoms, k denotes H atoms, and R denotes distance between atoms. The function, F (ρ), is an energy required to embed an atom into background electron density, ρ, and expresses interaction between the atom and the electron of the system. The functions, F Pt and F H , are energies necessary to embed a Pt atom and an H atom, respectively. The function, φ(R), is a two-body interaction potential between atoms at a distance of R and φ Pt−Pt (R), φ Pt−H (R) and φ H−H (R) are the two-body interaction potentials between two Pt atoms, a Pt atom and an H atom, and two H atoms, respectively. The functions, ρ i and ρ k , are electron densities at the position of Pt atom i and H atom k, respectively, contributed by the surrounding atoms. In the EAM method, the electron density is assumed to be the superposition of the electron density contributed by each surrounding atom mentioned below [25].
where ρ a Pt (R) is the electron density at the distance of R contributed by a Pt atom, and ρ a H (R) is that contributed by an H atom.
As shown in Eqs.
have to be determined to define the EAM potential. These functions and parameters were determined so that the energy of Pt bulk or H 2 molecule, and the electron density at specific sites (top, brg and fcc) of the Pt(111) surface are consistent with that obtained by DFT. These reference data were obtained by the commercial DFT program, DMol 3 . GGA/PBE was used as an exchange-correlation functional. A double-numerical plus polarization functions (DNP) was used as basis sets. Only valence electrons were considered. The convergence tolerance for self-consistent field (SCF) calculations was set at 1.0×10 −6 Ha. Spin was unrestricted during calculation. These determined functions were adjusted so that the dissociation barriers of the sites are consistent with those obtained by DFT calculation [12]. The validity of the potential was mentioned in the next section.

B. Verification of the EAM potential
In this subsection, the EAM potential mentioned above was verified by comparing the characteristics of the Pt(111) surface obtained by the EAM potential with those obtained by DFT. The schematic diagram of the system is shown in Fig. 1.
Firstly, electron densities at specific sites (top, brg and fcc) of the Pt(111) surface were calculated by the EAM potential and the results were compared with those obtained by DFT. The results are shown in Fig. 2. As shown in this figure, it is confirmed that the EAM potential can reproduce the electron density at each site obtained by DFT entirely. This figure also shows that there is a difference between the results of EAM and those of DFT at brg site near the surface. The difference, however, doesn't affect the dissociation barrier of the site because the position of the dissociation barrier is more than 1.5Å above the surface.
Next, the PESs of H 2 /Pt(111) system were obtained by the EAM potential and the dissociation barriers of the sites were obtained from the PESs. The distance between Pt atoms were set at the same distance of Pt bulk system and the effect of surface relaxation was not considered. An H 2 molecule was placed so that the center of mass of the molecule was set on a site of the surface and the potential energy was calculated by changing the bond length of the H 2 molecule, r, and the height of the center of mass of the molecule from the surface, Z. The intervals of ∆r and ∆Z were set at 0.03Å and the range was set at   Table  I (a). Those obtained by DFT in a previous research [12] are also shown in Table I (b). In this table, the orientation of a molecule was set so that the axis of the molecule is parallel to the Pt(111) surface and φ denotes the angle between the axis of the H 2 molecule and x axis of the coordinate system as shown in Fig. 1. As shown in Table I, the EAM potential can reproduce the dissociation barrier at each site of Pt(111) surface obtained by DFT. With regards to the position of the dissociation barrier, moreover, the difference between the results of EAM and those of DFT is around 10 %. Consequently, it is verified that the EAM potential constructed in this paper can reproduce the dissociation phenomena of H 2 /Pt(111) system accurately.

A. Simulation method
Simulations in which an H 2 molecule impinges to a Pt(111) surface were performed by the EAM potential obtained in Section II, and the motion of the H 2 molecule during collision was observed. The simulation system is shown in Fig. 3 The Pt surface was constructed of 500 Pt atoms aligned in five layers in the z direction, each of which consists of 10 Pt atoms (27.72Å) in the x direction and 10 Pt atoms (24.00Å) in the y direction. Periodic boundary conditions were imposed in the x and y directions. The initial positions of the Pt atoms were the same as the configuration of a Pt bulk system. The lattice constant of Pt was set at 3.92Å. For the temperature control to conserve the Pt surface at the target temperature, T [K], the method by which the velocity of the Pt atoms is scaled is often used. In this method, however, the velocity of the Pt atoms that directly interact with the impinging H 2 molecule is changed artificially, which greatly affects the motion of the H 2 molecules during collision. In the present paper, the temperature of the system was controlled by phantom atoms [30]. In the simulations in which the temperature of the surface was controlled, the velocity of the Pt and phantom atoms were given according to the Boltzmann distribution of the temperature of T [K] and the Pt surface was relaxed to reach an equilibrium state at T [K] before an H 2 molecule impinges on the Pt surface. The potential energy of this simulation was expressed as the first and fifth terms of Eq. (1). The Velocity Verlet algorithm [29] was used to integrate the equation of motion, where the time step was set at ∆t = 10 fs.
An H 2 molecule impinged at the specific sites (top, brg and fcc) of the Pt surface from a height of 5Å. The impinging energy, E tr was given to the H 2 molecule. The rotational energy, E rot was set at 0. The direction of the velocity of the molecule was given as normal to the Pt surface, and the orientation of the molecule was given randomly. The Velocity Verlet algorithm [29] was also used in the same manner as the simulation of the relaxation of the Pt surface. The time step, however, must be smaller than that of the simulation of relaxation of the Pt surface because the specific time of the rotational motion of H 2 molecule is much smaller than that of the vibrational motion of Pt atoms. In this paper, the time step was set at ∆t = 0.5 fs, and simulation was performed 40,000 time steps. It is confirmed that the total energy of the system was well conserved with sufficient accuracy.

B. Dissociation probability obtained by MD simulation
Firstly, an H 2 molecule impinged at the brg site of the Pt(111) surface at T = 350 K, and the motion of the molecule was analyzed. The impinging energy was given as E tr = 0.25 eV. Examples of the collision of the H 2 molecule with the Pt(111) surface are shown in Fig. 4.
As shown in Fig. 4 (a), the H 2 molecule collides with the Pt(111) surface and migrates. In this case, the length between H-H bond is larger, that is, the H 2 molecule dissociates during the migration on the Pt(111) surface. From the results, an H 2 molecule was determined to dissociate when the length of the H-H bond is more than 3.5Å on the standpoint of the force acting on H atom. On the other hand, in Fig. 4 (b), the H 2 molecule detaches from the surface after collision without dissociation though it passes over the dissociation barrier once. This case is determined to "dissociate" in the static analysis, for instance, transition state theory, because the molecule passes over the dissociation barrier. In this case, the classification of dissociation is different between static and dynamic analysis. The simulations mentioned above were performed 640 times by changing the orientation of the H 2 molecule at random at each site or impinging energy. The dissociation probabilities, P d , at specific sites against the impinging energy were obtained as the ratio of dissociative cases to all cases. The results are shown in Fig. 5.
As shown in this figure, the dissociation probabilities at brg and fcc site increase with the increase of impinging energy of the H 2 molecule. The result is reasonable because the impinging molecule having higher impinging energy can pass over the higher dissociation barrier. The dissociation probability at top site, however, decreases with the increase of impinging energy. This tendency is very different from the static analysis and the effect of motion of atoms on the dissociation probability is larger at top site than that at brg or fcc site. In this paper, the cause of the tendency was analyzed in detail. Firstly, the distributions of dissociation barrier at each site were obtained by changing the orientation of the H 2 molecule at random. The results are shown in Fig. 6. As shown in this figure, the distribution of dissociation barrier at top site has a peak near 0 eV. The peak corresponds to "no barrier dissociation", in which an impinging molecule can reach the strong interaction region with the Pt(111) surface, where the H 2 molecule can dissociate, without passing over the dissociation barrier. The amount of the peak is 25 % of all the probability. This means that 25 % of all the molecule impinging to a top site can reach the strong interaction region without passing over the dissociation barrier even if the molecules have only a very small impinging energy. The probability that an H 2 molecule can reach the strong interaction region, P c , at the top site was also obtained from the results of MD simulation. The results are shown in Fig. 7. The dissociation probability, P d , at top site shown in Fig. 5 are shown in this figure again.
As shown in this figure, almost all the molecules can reach the strong interaction region even at very small impinging energy. The result is not consistent with that of Fig. 6, which describes that only the 25 % of all the molecules having very small impinging energy can reach the strong interaction region at top site. The dynamics of the process was analyzed in detail. Consequently, it was found out that an H 2 molecule whose impinging energy is very small can reach the strong interaction region without dissociation barrier by changing their orientation due to the interaction from Pt(111) surface during the impingement even though they have the initial orientation of the molecule whose dissociation barrier is higher. The effect that a molecule changes its orientation by the interaction with the surface was called as "steering effect" and a previous research has implied that an H 2 molecule whose impinging and rotational energies are very small is easy to dissociate due to this effect [20]. The steering effect is also observed at brg or fcc site. At these sites, however, an H 2 molecule whose impinging energy is less than 0.1 eV cannot reach the strong interaction region because even the lowest dissociation barrier at these sites is more than 0.1 eV. At this condition, the steering effect does not have the effect on the dissociation probability. From these results, it can be said that the "steering effect" has much effect on the dissociation probability at the sites such as top site, where the orientation of no dissociation barrier exists. At the top site, however, the dissociation probability decreases with the increase of impinging energy though all the molecules can reach the strong interaction region. The distribution of orientation in the case that the molecules pass over the dissociation barrier and that the molecules dissociate were analyzed. The results show that the dissociation barrier becomes lower when the axis of molecule is normal to the surface, while the molecule can dissociate more easily when the axis of molecule is parallel to the surface. The molecule can change its orientation normal to the surface firstly so that it can get low dissociation barrier, and then change the orientation parallel to the surface so that the molecule can dissociate easily by steering effect when the impinging energy is very small. The molecule having higher impinging energy, however, does not have enough time to change their orientation and they detach from the surface after collision, as shown in the case (b) of Fig. 4. As shown in these results, it can be said that the dissociation probability must be considered including the motion of atoms or molecules, especially in the case that the dissociation barrier is very small.

IV. CONCLUSION
In this paper, the EAM potential was reconstructed to reproduce the characteristics of dissociation barrier at specific sites and the difference of the characteristics of dissociation phenomena at the sites were analyzed. A Pt(111) surface was chosen as the metal surface, and hydrogen was chosen as the gas molecule of the system. Functions or parameters of the EAM potential were determined so that the characteristics of the interaction potential obtained by the EAM were consistent with those obtained by DFT. It was confirmed that the EAM potential can reproduce the distribution of electron density or the dissociation barrier at specific sites (top, brg and fcc) obtained by DFT. Using this potential, simulations that an H 2 molecule impinged to a Pt(111) surface were performed and the motion of the molecule during the collision and dissociation was observed. The simulations were performed 640 times, and the dissociation probability of an H 2 molecule on a Pt(111) surface was obtained. The dissociation probability at brg and fcc site increases with the increase of impinging energy, which is consistent with the static analysis such as transition state theory. At top site, however, the dissociation probability decreases with the increase of impinging energy. It was confirmed that the "Steering effect" affects the dissociation probability at top site and it can be said that the dissociation probability must be considered including the dynamic effects, especially in the case that the dissociation barrier is very small.