Molecular Simulation of Dissociation Phenomena of Gas Molecule on Metal Surface

Dissociative adsorption phenomena often occur in various fields of engineering, such as oxidation-reduction reactions, cleaning, adhesion, plating, plasma etching, sputtering, and tribology. These phenomena that involve surface reactions have attracted much attention and are analyzed both experimentally and numerically. However, when the surface has structures on the molecular scale, and the scale is not small enough compared to the system, the characteristics of a surface reaction cannot be sufficiently expressed macroscopically, for example by the rate equation. It is extremely difficult to analyze the characteristics of a nanoscale system experimentally due to the scale. Therefore, analysis by numerical calculation, in which the system structure and its electronic state are treated comprehensively, is more effective.


Introduction
Dissociative adsorption phenomena often occur in various fields of engineering, such as oxidation-reduction reactions, cleaning, adhesion, plating, plasma etching, sputtering, and tribology.These phenomena that involve surface reactions have attracted much attention and are analyzed both experimentally and numerically.However, when the surface has structures on the molecular scale, and the scale is not small enough compared to the system, the characteristics of a surface reaction cannot be sufficiently expressed macroscopically, for example by the rate equation.It is extremely difficult to analyze the characteristics of a nanoscale system experimentally due to the scale.Therefore, analysis by numerical calculation, in which the system structure and its electronic state are treated comprehensively, is more effective.
To analyze the surface reaction of these systems accurately, it is necessary to solve the electronic state by the first principle calculation based on quantum mechanics and to then obtain the energy state.The Molecular Orbital (MO) method is most accurate one.However it takes much calculation time and it is impossible to analyze the dissociation phenomena of gas molecule on metal surface because metal has many electrons.Recently, density functional theory (DFT) is one of the most commonly used methods for this process (Parr & Yang, 1989;Satoko & Onishi,1994).Based on the theory that the state of a system is expressed by the functional of the density distribution of the electron, this method can calculate a system state faster than methods that calculate the wave functions of each electron like MO.In the process of surface reaction analysis, this method is applied in various situations, such as specifying the reaction paths from the potential energy surface obtained by the method and calculating the reaction probability at the surface by the value of the absorption/activation energy and the transition state theory (Steinfeld et al., 1989).However, the effects of the motion of gas molecules impinging on the surface and the motion of surface atoms on the surface reaction cannot be considered because this method is applied under the assumption that the temperature is 0 K (fixed atom)., In order to analyze the flow dynamics, including the surface reaction, accurately, a method that considers the interaction between the electronic state of the system and the motion of atoms or molecules for which the time/space scale is distant, must be used.
Molecular Dynamics (MD) is one of the most suitable method for simulating the motion of atoms.In this method, the force between atoms is first obtained, and then the positions or velocities of the atoms are simulated by a time marching scheme.The quantum molecular dynamics (QMD) method, in which the interaction between atoms of a system is obtained according to the quantum calculation mentioned above, is the most precise method.In particular, the Car-Parrinello method (Car & Parrinello, 1985), in which the potential force of a system is calculated using the DFT to calculate the electronic state, is applied to analyze the process of oxidation-reduction reaction at the platinum surface (Wang & Balbuena, 2004;Jinnouchi & Okazaki, 2003).In general, the Car-Parrinello method is known as a first principle quantum molecular dynamics method.However, despite its precision, the Car-Parrinello method is not practical for analyzing a flow phenomenon as a statistical behavior of numerous motions of atoms or molecules because of the enormous calculation load.For the analysis of flow phenomena with surface reaction, the application of the multiscale method, in which only smaller-scale systems are analyzed by precise quantum calculation and the characteristics that affect surface reaction are modeled, is more appropriate.Flow phenomena with larger-scale surface reactions can then be analyzed by the MD method using this model, rather than treating the entire system using quantum calculation.
The tight binding (TB) method, which greatly improves the computation speed, is contrived by simplifying the first principle quantum molecular dynamics method.Regarding this method, as mentioned later, it has been reported that the computation speed was improved 5,000 times beyond that of the first principle quantum molecular dynamics method by calculating the electronic state using the extended Hückel method (Yonezawa et al, 2001a(Yonezawa et al, , 2001b)).Since this method enables faster computation and still has the characteristics of calculating the electronic state of a system according to quantum theory, the tight binding method is often applied within the field of chemistry for analyzing the surface reaction dynamics of relatively large systems.Using the tight binding method, the process of the production reaction of water in a fuel cell at the solid-gas interface has been simulated (Ishimoto et al., 2006).In addition, a hybrid method, which combines the tight binding method and the classical molecular dynamics method, is also being developed.
The tight binding method is still not sufficient for large-scale calculations that deal with the "flow" of a system because the calculation procedure is complicated because the electronic state of a system is calculated according to quantum theory.For calculations that deal with the statistical quantity of atomic motion, another method that determines the potential function, which is used in the classical molecular dynamics method using the results of the density functional theory, is also applied frequently.For instance, the potential, which is a function of the position or the orientation of an impinging molecule, is contrived by fitting the potential energy surface of a diatomic molecule on a transient metal surface at various orientations obtained by the density functional theory using an analytic function.The potential is often used to analyze the dissociation phenomena of hydrogen on a Pt or Pd surface (Beutl et al., 1995;Olsen et al., 1999Olsen et al., , 2002)).In this method, however, the motion of surface atoms cannot be considered, and therefore the effect of thermal motion of surface atoms on dissociation phenomena cannot be analyzed using this method.Owing to this defect, it has often been reported that the dissociation probability obtained by this method cannot be used to reproduce experimentally obtained data (Vincent et al., 2004).
The embedded atom method (EAM) is a scheme that treats the interaction between gas atoms and a metal surface by considering the effects of the surface electrons (Daw & Baskes, 1983;1984).The method is based on DFT, and the potential energy of the system is expressed as a sum of the energy embedded an atom in the electron density of the surface and pair-interaction energy.In EAM, the electron density of the system can be reflected in the interaction potential, and therefore the motion of an atom on a metal surface can be simulated accurately.Moreover, the method has the advantage of a smaller computational load than quantum molecular dynamics (QMD), which is also based on DFT.EAM has often been used to analyze the motion of atoms on a transient metal surface (Baskes, 1992;Baskes et al., 2007).However, there has been little research of the application of EAM to dissociation phenomena.To apply EAM to the analysis of dissociation probability, incorporating the motion of the metal surface atoms, brings the simulations closer to the real system.Moreover, this method can be easily expanded to a more complicated surface.
In this chapter, the methods used for multi-scale analysis of flow phenomena including a surface reaction are described, and a typical example in which these methods are applied to the analysis of flow phenomena including a surface reaction is presented.In Section 2, the simulation method of Density Functional Theory and Embedded Atom Method are described.In Section 3, the analysis of dissociation phenomena by the EAM is discussed.Section 4 summarized the chapter.

Simulation methods
This section describes the simulation method which was used in this chapter.Especially Density Functional Theory (DFT) and Embedded Atom Method (EAM) are described.

Density Functional Theory (DFT)
The density functional theory (DFT) is the method in which various values obtained by wave functions of a system and operators are expressed by the functional of electron density of the system.In this method, the energy of the system is obtained exactly.In the analysis of the surface reaction, macroscopic values such as dissociation probability are obtained from the potential energy surface of the system obtained by this method.Moreover, this method is effective for obtaining the database necessary to determine the parameters for the tight binding method or the embedded atom method mentioned below.In this subsection, the outline of the theory of the method is explained below.For details, the reader should refer to some references (Parr & Yang, 1989;Satoko & Onishi,1994).
Let us consider the system that consists of M nuclei and N electrons.Thus, the nuclei are fixed and only the state of the electron is considered (Born-Oppenheimer approximation).The Schrödinger equation, which expresses the ground state of a system, is given by: where Ψ o and E o denote a wave function that expresses the ground state of a system without degeneracy and the energy of electron at this state, respectively, and H denotes an operator, called the Hamiltonian, that expresses the total energy of the system.The operator is expressed using the kinetic energy of the electron, K, the Coulomb interaction between electrons, V ee , and the interaction with the external force field, V ex , as ee ex The external force field, V ex , in Eq. ( 2) is expressed as the sum of a function of only the position of each electron, v(r).The formula of the function, for instance, is expressed in atomic units as assuming that v(r) is a Coulomb interaction from the nuclei of the system.In Eq. ( 3), r i and R k denotes the position of electron i and that of nucleus k reduced by the Bohr radius, respectively, and Z k denotes the charge of nucleus k reduced by elementary electric charge (e=1.602×10−19 C).The Coulomb interaction between electrons is expressed as the sum of the Coulomb interaction between electron i and j as ee 1 1 The solution of Eq. ( 1) satisfies the equation below: where <Φ|Φ>=1, Φ is an eigenfunction of H, Ψ o is a function that minimizes the total energy of the system among the eigenfunctions of H.As mentioned above, when the formula of the interaction from the external force field, v(r), is determined, the Hamiltonian is obtained by Eq. ( 2), and therefore the wave function Ψ o and the total energy E o of the system are obtained by Eq. ( 5).Namely, the state of the system is a functional of the interaction with the external force field, v(r).Moreover, the electron density at position r is obtained by The basic principle of the density functional theory is that an interaction with the external force field, v(r), corresponds to a certain electron density, ρ(r), and vice versa, at the ground state without degeneracy.That is, an interaction from the external force field, v(r), is obtained by a functional of electron density, ρ(r) (Parr & Yang, 1989).The electron density at the ground state, ρ(r), is obtained as described below, rather than by Eq. ( 5), so as to minimize the total energy of the system.
In density functional theory, the electron density is obtained by Eq. ( 7).The total energy of the system, E[ρ], is expressed as where the sum of the kinetic energy, K [ρ], and the interaction between electrons, V ee [ρ], is expressed by F [ρ].The electron density at the ground state, ρ(r), minimizes E[ρ] in Eq. ( 8) and satisfies the following Euler equation: while considering the conservation of the number of electron as where μ is Lagrange's undetermined multiplier.In other words, it is necessary to solve the electron density, ρ(r), that satisfies Eq. ( 9) in the density functional theory.However, it is difficult to solve Eq. ( 9) because the interaction between electrons is included in the functional F[ρ].In the density functional theory, the electron density of the system, ρ(r), is obtained by solving formula of a virtual system in which a certain effective potential, v eff , affects electrons, and an interaction between electrons is not considered in the self-consistent field method.
Let us consider a system that consists of N electrons without interaction between electrons.The Hamiltonian of the system is expressed as where v v (r) denotes an interaction from an external force field and is a function of only the electron position.The wave function of the system at the ground state is expressed by the Slater determinant of the one-electron wave function, φ i , as The one-electron wave function, φ i , corresponds to the function in which the eigenenergy is lowest among all of the wave functions satisfying Eq. ( 13).
The electron density of the system at the ground state, ρ(r), is obtained by Eq. ( 6) using the wave functions obtained by Eq. ( 12).According to the theory mentioned above, the total energy of the system is a functional of electron density, ρ(r), and is obtained by However, the electron density, ρ(r), obtained from the wave function in Eq. ( 12) is that which minimizes the energy of Eq. ( 14) and satisfies the following Euler equation: because the electron density is that at the ground state, where K v [ρ] in Eq. ( 14) or Eq. ( 15) shows the term of the kinetic energy of the virtual system.Namely, if Eq. ( 9), which includes the interaction between electrons in the functional F [ρ], can be changed to Eq. ( 15), the electron density of the system can be obtained by Eqs. ( 6), ( 12), and ( 13) using an interaction with the external force field of the virtual system, v v (r).
The functional in Eq. ( 9), F[ρ], is expressed as where J[ρ] denotes a classical Coulomb interaction between electrons and is expressed as and the functional, denotes the sum of the non-classical terms of interactions between electrons, (V ee [ρ]-J[ρ]) and the difference in kinetic energy between the real system and virtual system, (K , is called the exchange-correlation energy.Using Eqs.( 9) and ( 16), we obtain where denotes the Coulomb interaction of the electrons and denotes the exchange-correlation energy.Assuming the value in Eq. ( 19), v(r)+(r)+v xc (r), to be an interaction with the external force field in the virtual system, v v (r) in Eq. ( 15) (referred to as effective potential, v eff (r)), namely, Equation ( 19) is the same as Eq. ( 15), and the electron density of the real system can be obtained using Eqs.( 6), (12), and (13).In order to obtain the effective potential, v eff (r), in Eq. ( 22), the term of the exchange-correlation energy, v xc (r), must be determined.The values are given as a functional of electron density.Several types of functional have been proposed in previous studies.The local density approximation (LDA), which is the most well known of these functionals, is obtained by assuming that the electron is uniformly distributed at a small region of the system.Other functionals, such as the generalized gradient approximation (GGA), which includes the information of the gradient of the electron density, or the local spin density approximation (LSDA), which is applied to the case of spin polarization, are developed considering the compensation between the accuracy and the calculation time.After the exchange-correlation energy is determined, the electron density is obtained by Eqs. ( 6), ( 12), ( 13), ( 20) and ( 22).However, it is necessary to obtain the electron density by the convergence scheme because (r) and v xc (r) are functionals of electron density.A flowchart of this scheme is shown in Fig. 1.The advantage of this method is the dramatic decrease in the calculation time needed to directly calculate the energy from the electron density without solving the wave functions of the system.Namely, the exact solution can be obtained from the electron density of an independent N-electron system without solving the Schrödinger equations of the 3N space including the interaction between electrons by including the interaction between electrons in the empirical exchange-correlation energy.

Embedded Atom Method (EAM)
DFT can calculate the energy barrier accurately and the dissociation probability can be obtained from the information of energy barrier.However, it is based on static transition state theory (TST), while in real system atoms and molecules move.To consider the motion of atoms, the FPMD must be used.However it costs too much calculation load to analyze the statistical characteristics of dissociation phenomena.It is necessary to use a faster method to analyze flow phenomena obtained as the statistical characteristics of the motion of atoms of the system.In this subsection, an outline of the embedded atom method is presented.The theory of this method is based on density functional theory, and this method is often applied to the analysis of the motion of atoms on a metal surface (Daw and Baskes, 1984).
In the embedded atom method, the energy of a system is expressed as the sum of the energy necessary to embed an atom into the background electron density and interaction between nuclei.The former corresponds to the energy from the electrons of the system, and the latter corresponds to the energy from the nuclei of the system.The energy of a system is then expressed as where F i () denotes the energy necessary to embed atom i into the background electron density ρ, and  ij (R ij ) denotes the core-core pair repulsion between atoms i and j separated by the distance R ij .In addition, ρ i denotes the electron density at atom i due to the remaining atoms of the system.The electron density is expressed by the superposition of the electron density of each atom as where ρ j a is the electron density contributed by atom j.It is necessary to determine the functions, F i (ρ),  ij (R ij ), and ρ j a (R ij ) included in Eqs. ( 23) and ( 24) in order to use the embedded atom method.An example of how to determine these functions or parameters is described below.
The electron density can be obtained by a quantum chemical calculation, such as the DFT or molecular orbital calculation.McLean and Clementi calculated the wave functions of various atoms by the Roothaan-Hartree-Fock method and generated a data table (Clementi & Roetti, 1974;McLean & McLean, 1981).The reader can obtain the electron density by the data table and Eq. ( 6).With respect to the function, F i (ρ i ), the energy of various atoms embedded in a homogeneous electron density as a function of electron gas density have been obtained by the density functional theory (Puska et al, 1981;Norskov, 1982).The core-core pair-repulsion,  ij (R ij ), is often determined such that the crystal structure of the bulk system can be reproduced by the EAM using the determined values of ρ j a (R ij ) and F i (ρ i ).These functions are determined for every type of atom.
The force acting on atom i is obtained by differentiating Eq. ( 23), as follows:

Results and discussion
In this section the example of the analysis of dissociation phenomena of gas molecule on metal surface is introduced.We use the dissociation of H 2 molecule on Pt(111) surface, which is important not only for the fundamental system of dissociation phenomena but also for the reaction at anode side in Polymer Electrolyte Fuel Cell.

Potential energy surface
Analysis of the potential energy surface is important when discussing dissociation phenomena on a metal surface.From this potential energy surface, the values of the dissociation path or the dissociation barrier are calculated, and the dissociation probability is then obtained.In this subsection, a method for calculating the potential energy surface of a platinum-hydrogen system according to the density functional theory and the abovedescribed analysis will be shown.All the calculations were conducted using the commercial density functional theory (DFT) program, DMol 3 .LDA/VWN and GGA/PBE were used as exchange-correlation functionals for geometry optimization and energy calculations, respectively.The double-numerical plus polarization functions (DNP) were used as the basis sets for all atoms.All electrons were considered in the calculation.The convergence tolerance for the self-consistent field (SCF) calculation was set at 1.010 -6 Ha and the convergence tolerance for geometry was set at 1.010 -5 Ha.Smearing was set at 5.010 -3 Ha for both calculations to improve convergence of the calculations.Spin was unrestricted during calculations.
Figure 2 shows the calculation model.In the calculation of the typical density functional theory, the surface is often expressed by three molecular layers in order to have a reasonable calculation load.At each layer, a platinum atom is placed at a location to its Miller indices correspond.Figure 2 shows the arrangement of the (111) surface.Because the calculation load determines the number of atoms representing one layer of the surface, the calculation of the direction of the surface is generally performed under the assumption of the periodic boundary condition, expressing the surface with 22 = 4 atoms.In the same figure, the basic cell is denoted by the white line.In the density functional theory, calculation is performed with periodic boundary condition in all directions.Figure 2 shows the technique used to avoid interactions normal to the surface by setting a long vacuum region.
Sometimes the calculation model that passed through the process described above is applied as a metal surface, but this causes instability when using the structure of bulk as is for the surface.At the actual surface, the structure changes slightly, causing surface relaxation.In the case of calculating a more realistic surface, reconstruction of the surface atom arrangement to achieve the most stable state after geometry optimization is one option.The LDA method is usually used for geometry optimization.
The orientation of molecule toward the surface is determined by locating the center of the molecule on a particular site on the surface (location of the molecule on a surface) and then defining the Euler angle of the molecular axis (, θ) in the space coordinate system.For the determined orientation, the energy of a system is calculated by changing the distance between nuclei, r, and the distance between the first surface layer and the center of the molecule, Z.The GGA method is most commonly used for the energy calculation.Repeated energy calculation while changing the values of r and Z enables the contour of the energy (potential energy surface) to be drawn under the condition that the lateral axis is r and the longitudinal axis is Z.Examples are shown in Fig. 3.The left-hand side represents the PES at =0° and θ=90° of top site, and the right-hand side represents that at =0° and θ=90° of brg site.In this calculation, r is changed at 0.5 Å intervals in the range between 0.5 Å and 2.5 Å, and Z is changed at 0.5 Å intervals in the range between 0.5 Å and 2.5 Å. Fig. 2. Simulation system for the potential energy surface.From the PES, the dissociation path or the height of the dissociation barrier is obtained as described below.Assuming that the hydrogen molecule passes through the route of lowest energy from the state before dissociation (State A) to that after dissociation (State B) in Fig. 3, the route corresponds to a ridge of the potential energy surface.Therefore, the route of reaction is determined by obtaining the ridge from State A to State B. The white dashed lines in Fig. 3 show the dissociation path of each potential energy surface.The energy on the dissociation path is calculated by defining the energy of State A as 0. The potential energy increases as the molecule approaches the surface.The maximum energy on the dissociation path is called the "dissociation barrier".Molecules of less energy than the dissociation barrier cannot dissociate.
In the transient state theory, the dissociation probability is obtained from the height of the dissociation barrier.The probability that the kinetic energy of molecules impinging on the surface is in the range of [e tr , e tr +de tr ] is obtained by Boltzmann distribution as where k denotes the Boltzmann's constant.In this theory, all of molecules that have a more energy than the dissociation barrier are considered to dissociate.This probability is obtained by integrating Eq. ( 26) in the range of [E b , ∞], as follows:

Energy transfer between gas molecule and solid atoms
In the analysis of dissociation probability by the transient state theory and the density functional theory mentioned above, the molecules that are reflected from the surface without dissociation after passing over the dissociation barrier are not considered when obtaining the dissociation probability because in this theory all molecules that have an energy greater than the dissociation barrier are considered to dissociate.Moreover, the effect of the direction of the motion or the rotational energy of molecules on the dissociation probability is not considered because the dissociation probability is determined by only the magnitude of the translational energy of the molecule.However, the potential energy surface is likely to change due to the motion of surface atoms.It is necessary to consider these motions of atoms in order to make the analysis more real.In this subsection, the method used to obtain the dissociation probability while considering the motion of atoms is described.The embedded atom method is used as the simulation method.The method used to determine the embedded function or the core-core pair-repulsion function is also described.
In the EAM, the potential energy, E pot , of the system, which consists of a Pt surface of N atoms and an H 2 molecule, is obtained by where subscripts i and j denote the given number of Pt atoms, subscript k denotes the given number of H atoms, and R denotes the distance between atoms.The functions F Pt (ρ i ) and F H (ρ i ) denote the energy necessary to embed a Pt atom or an H atom, respectively, into the background electron density, ρ i , and  Pt-Pt (R ij ),  H-H (R 12 ) and  Pt-H (R ik ), denote the core-core pair-repulsion potential between Pt atoms, between H atoms, and between a Pt atom and an H atom, respectively.In addition, ρ i and ρ k denote the electron density at Pt atom i or H atom k, respectively, by the remaining atoms.By defining these functions such that the results obtained by the EAM reproduce those obtained by the DFT, calculation by the EAM is faster, and the accuracy in the EAM is sufficient.The details of the form of the function are omitted herein.For these details, the reader should refer to some references (Tokuamsu and Ito, 2007;2011).
Next, an example of simulating the dissociation phenomena of H 2 on the Pt surface by the EAM potential determined mentioned above is described.A schematic diagram is shown in Fig. 4. The Pt(111) surface consists of 10103 = 300 Pt atoms in the x, y, and z direction, respectively.A periodic boundary condition is imposed in the x and y directions.The initial position of atoms is the same as that of a bulk crystal structure of Pt.The lattice constant of Pt is set at 3.92 Å.
When the temperature of the surface is controlled, it is often used to control the velocity of the atoms of the system.However, this method affects the motion of dissociating molecules because the velocity of atoms that directly interact with the dissociating molecules is changed artificially.For this reason, in this example, the temperature of the surface is controlled by phantom molecules (Blomer and Beylich, 1996).The initial velocities of the Pt atoms and the phantom atoms are given at random, according to the Boltzmann distribution at the temperature of T [K], and the system is relaxed toward the equilibrium state when the temperature of the system is controlled.The spring constant for phantom atoms, k, was set at 46.8 N/m, and the coefficient of the dumper, , was set at 5.184×10 -12 kg/s.The cutoff distance of the interaction was set at 15 Ǻ.Time integration was calculated by the leap-frog method (Allen & Tildesley, 1986) with a time interval, Δt, of 1 fs.Before simulating the dissociation of H 2 molecules, the relaxed Pt(111) surface, whose temperature was controlled to a target temperature, had to be obtained.The simulation was performed until the system reached equilibrium at the target temperature.In this simulation, the target temperature was set at T=350 K.
H 2 molecules impinged upon the top, brg, or fcc sites of the relaxed Pt(111) surface from a height of 5 Ǻ.An initial translational energy, E tr , was given to the H 2 molecule normal to the surface as the impinging energy.Neither rotational nor vibrational energy was given to the H 2 molecule.The orientation of the impinging H 2 molecule was given at random.
A typical example of the difference from the initial energy of kinetic energy of Pt atoms, kinetic energy of the impinging H 2 molecule, potential energy obtained in Eq. ( 28) and total energy are shown in Fig. 5. Bold line shows the total energy of the system and square, circle and triangle denotes the E pot obtained in Eq. ( 28), kinetic energy of Pt atoms and that of the impinging H 2 molecule, respectively.As shown in this figure, the energy transfers between each degree of freedom but total energy of the system was well conserved relative to the difference from initial energy of each degree of freedom.Simulations were performed for 40,000 steps or 10,000 steps, depending upon whether the impinging energy was smaller or larger than 0. In Fig. 6 (a), the impinging energy was set at E tr =0.25 eV.In this case, the velocity of impinging molecule decreases during t=5060 fs and the molecule collide with the Pt surface, which implies that the molecule passes over a dissociation barrier.The distance between H atoms becomes longer while the H 2 molecule migrates on the surface.In the present simulations, H 2 molecules having a distance between H atoms larger than 3.5 Ǻ were considered to dissociate. Figure 6 (b) shows the case not dissociating after collision.The impinging energy was set at E tr =0.1 eV.As shown in this figure, the molecule can pass over the dissociation barrier and collides with the surface.However, the molecule reflects and departs again directly.In this case the vibrational energy is excited because of the strong interaction to the surface.

Dissociation probability
The simulations mentioned in Sec.3.2 were performed 640 times, fixing the impinging energy and the site on the Pt(111) surface, and changing the orientation of the H 2 molecule at random.The dissociation probability was obtained against impinging energy at top, brg, and fcc sites from the ratio of dissociated to all cases.The initial configuration of the Pt(111) surface was changed in every simulation.The impinging energy was varied from 0.01 to 0.5 eV.The "dynamic" dissociation probabilities, P d , against impinging energy are shown in Fig. 7 (a).
By the way, the "static" dissociation probabilities, P s , can be obtained from the interaction potential between an H 2 atoms and a Pt surface.The detailed method was described in a reference (Tokumasu and Ito, 2011).The static dissociation probability is shown in Fig. 7 (b).Fig. 7 (a) shows that at brg and fcc sites, the dynamic dissociation probability increases with increasing impinging energy of the H 2 molecule, and becomes nearly constant at higher impinging energy.The tendency is similar to the static dissociation probability shown in Fig. 7 (b), implying that dynamic effects are not important to the dissociation at brg or fcc sites.
The dissociation probability at top site, however, decreased with increasing impinging energy.This tendency is very different from the static dissociation probability shown in Fig. 7 (b).Moreover,Fig. 7 (b) shows that the dissociation probability when the impinging energy is near 0 is about 25 %, which means that 75 % of the impinging molecules cannot pass over the dissociation barrier.However, Fig. 7 (a) shows that 90 % of the impinging molecules dissociates when the impinging energy is near 0. This contradiction shows the "dynamic effects" on the dissociation probability (Tokumasu & Ito, 2011), and therefore the effects at the top site are larger than that at brg or fcc sites.This contradiction was analyzed in detail.shows that H 2 molecules with very low impinging energy can reach the chemisorption layer by changing their orientation via interaction with the Pt(111) surface.These impinging molecules have no dissociation barrier, even if they have an initial orientation at which the dissociation barrier is large.The change in orientation of an impinging molecule via interaction with the surface to reduce the dissociation barrier is called the "steering effect (Darling & Holoway, 1994;Darling et al, 1998;Gross et al. 1995).Previous research indicated that an impinging molecule is easier to dissociate than expected in a static manner due to this effect when the impinging energy and rotational energy of the impinging molecule is very small (Darling et al, 1998).The steering effect is also observed at brg and fcc sites.The impinging molecules with very low impinging energy, however, cannot pass over the dissociation barrier regardless of orientation at brg or fcc sites, because the minimum dissociation barrier is more than 0.1 eV at these sites.For this reason, the steering effect is unimportant to dissociation phenomena at brg or fcc sites.The steering effect becomes remarkable only at sites where a molecule can reach the chemisorption layer with a very small dissociation barrier (or no dissociation barrier) at specific orientations, such as at the top site.
The dissociation probability, however, decreases rapidly with increasing impinging energy, although all impinging molecules reach the chemisorption layer.The tendency is different from that of brg and fcc sites, where almost all molecules that reach the chemisorption layer dissociate.The distribution of orientation of impinging molecules at the top site, when the molecules pass over the dissociation barrier and when the molecules dissociate, were investigated at very low impinging energy.The results showed that the orientation of the impinging molecule at which the dissociation barrier is very low is normal to the surface, and the orientation at which the molecules dissociate is parallel to the surface at the top site.Almost all the molecules can dissociate at the top site when the impinging energy is very low, because the molecules can change their orientation by the steering effect to easily pass over the dissociation barrier to reach the chemisorption layer (by aligning normal to the surface).They can then readjust their orientation, again by the steering effect, such that they can easily dissociate (by aligning parallel to the surface).With increasing impinging energy, however, the molecules have sufficient energy to reach the chemisorption layer without the steering effect, and do not have time to change their orientation to one in which they can easily dissociate.Therefore, they depart from the surface without dissociation after collision.
In dissociation at the top site, two steering effects are important: once when the molecule reaches the chemisorption layer, and again when the molecule dissociates.Therefore, the dominant factor in dissociation at the top site is the motion of the H 2 molecule at the chemisorption layer by steering effect, not the probability that the H 2 molecule reaches the chemisorption layer.The dynamic effects on dissociation probability are very important when dissociation phenomena at the top site are considered.

Dissociation Probability [-]
Impinging Energy [eV] : P d : P c 8 Fig. 8. Dissociation probability at brg and fcc sites, P d , and the probability of reaching the chemisorption layer, P c (top site).Circles and squares show P d and P c , respectively.

Comparison with experimental results
To check the valildity of the simulation, the dissociatino probability against the impinging angle was simulated and the resuts were compared with experimental results (Luntz et al., 1990).In this study, D 2 molecules were used as impinging gas molecules because in the experimental study with which we compared our results, D 2 molecules were used.The EAM potential constructed for the H 2 molecule was used for the D 2 molecule since the state of electrons can be regarded as the same for both of these molecules.Simulation system were almost the same as that obtained in Sec.3.1, but the incident angle was given to the impinging molecule.The incident polar angle was varied from 0 o (normal to the surface) to 60 o and the incident azimuthal angle was varied in a uniform random manner.Moreover, the initial position of H 2 molecule is given at random.The rotational energies were given according to Boltzmann distribution.The detailed method was described in a reference (Koido et al., 2011).
The results of MD simulations are shown in Fig. 9.The dissociation probability, P d , for each incident angle against initial translational energy is shown.It shows that P d increases along with the initial translational energy at all level of polar angle of incidence.When the angle is normal to the surface (θ i = 0°), even at very low translational energy, P d has a certain value and it increases along with the initial translational energy.This shows that there is no obvious translational energy threshold for dissociation probability to rise.It means that there is no energy barrier for dissociation when a gas molecule approaches to the surface with the orientation that takes the minimum potential among its possible potential energy surfaces.At a larger initial translational energy than E tr = 0.25 eV, P d does not increase.
At larger angles of incidence, P d increases more slowly with initial translational energy similar to the normal incidence case and stops increasing at higher initial translational energy.It should be noted that when the polar angle of incidence is off normal to the surface, at the lowest initial translational energy calculated (E tr =0.0025 eV), P d is a little higher than the value at little higher translational energy.The rise at very low translational energy suggests that the "Steering Effect" at the top site (Tokumasu and Ito, 2011) is important.When translational energy is very low, the possibility of a molecule to have both low translational energy and low rotational energy is relatively high.Moreover, it allows good chance for a molecule to rotate to the orientation that experiences low energy barrier to dissociate near the surface.It can be said that even though the result is the average of collisions on all the surface area, the "Steering Effect" at a particular site cannot be neglected.
These MD results are compared with experimental data (Luntz, 1990) at each incident polar angle.Although the MD results and the experimental data do not agree very well, they both show the following trends: P d increases along with initial translational energy, P d decreases along with increasing polar angle of incidence, and no initial translational energy threshold for dissociation probability is observed when the molecules approach along the surface normal.
If we see the hydrogen dissociation on other metal surfaces, the trends may be quite different.It was reported that hydrogen dissociation on Cu(111) surface requires at least around 0.5 eV of large activation energy (Luntz, 2009;Gross, 1998).Moreover, it was also reported both theoretically and experimentally that the dissociation probability is high at 0 translational energy and decreases along with the translational energy and increases at higher energy again on Pd(100) surface (Luntz, 2009;Gross, 1998).From the above discussion, MD simulations using the constructed EAM potential are capable of simulating molecular beam experiments to a certain degree.
The rise in the dissociation probability observed at low translational energies in the MD results cannot be compared with the experiments due to lack of experimental data.

Summary
In order to analyze the dissociation phenomena of gas molecule on metal surface which includes chemical reactions, it is necessary to perform quantum mechanical calculations considering the electronic states of materials.However, since this is not practical due to the large calculation time and therefore the analysis based on molecular dynamics, which needs dramatically small calculation time than quantum calculation, is desired.Therefore, a quantum mechanical method is first used to analyze the (nanoscale) electronic state, which dominates the reaction phenomena, and, using the obtained information, the empirical parameter applied to molecular dynamics is determined.The behavior of the molecule or its statistical (micro-scale) quantity is then analyzed.This analysis is referred to as "multi-scale analysis".In multi-scale analysis, it is important to properly understand the nanoscale phenomena that dominate the micro-scale phenomena.

Acknowledgment
I am very grateful to Dr. Koido for his results of comparing MD simulation with experimental data.All the simulations of this chapter were performed by a supercomputer system of Institute of Fluid Science, Tohoku University.

Fig. 1 .
Fig. 1.Flowchart for the calculation of electron density by the density functional theory.

Fig. 3 .
Fig. 3. Contour of the potential energy surface.The contour is plotted at 0.2 eV intervals.The white dashed lines show the dissociation path.

Fig. 5 .
Fig. 5. Typical example of differential from initial energy of degree of freedom.Bold line denotes the difference of initial energy of total energy.Square, circle and triangle denote those of potential energy, kinetic energy of Pt atoms and kinetic energy of H 2 molecule, respectively.

Fig. 6 .
Fig. 6.Distance from surface of the center of mass of the H 2 and the distance between H atoms. (a) the impinging H 2 molecule dissociates.(b) the impinging H 2 molecule does not dissociate.

Fig. 7 .
Fig. 7.The dissociation probabilities at each site.Circles, squares, and triangles show the results at top, brg, and fcc sites, respectively.(a) Dynamic dissociation probability, (b) Static dissociation probability.Initially, the probability that an impinging H 2 molecule reaches the chemisorption layer, P c , was obtained for the top site.The result is shown in Fig.8.As shown in this figure, all impinging molecules reach the chemisorption layer, even at very low impinging energy, unlike the results shown in Fig.7 (b).The behavior of the impinging H 2 molecules clearly shows that H 2 molecules with very low impinging energy can reach the chemisorption layer by changing their orientation via interaction with the Pt(111) surface.These impinging molecules have no dissociation barrier, even if they have an initial orientation at which the dissociation barrier is large.The change in orientation of an impinging molecule via interaction with the surface to reduce the dissociation barrier is called the "steering effect(Darling & Holoway, 1994;Darling et al, 1998;Gross et al. 1995).Previous research indicated that an impinging molecule is easier to dissociate than expected in a static manner due to this effect when the impinging energy and rotational energy of the impinging molecule is very small(Darling et al, 1998).The steering effect is also observed at brg and fcc sites.The impinging molecules with very low impinging energy, however, cannot pass over the dissociation barrier regardless of orientation at brg or fcc sites, because the minimum dissociation barrier is more than 0.1 eV at these sites.For this reason, the steering effect is unimportant to dissociation phenomena at brg or fcc sites.The steering effect becomes remarkable only at sites where a molecule can reach the chemisorption layer with a very small dissociation barrier (or no dissociation barrier) at specific orientations, such as at the top site.

Fig. 9 .
Fig. 9. Dissociation probability predicted by MD simulations (black symbols) as a function of initial translational energy, E tr for surface temperature, T s = 295 K, for different incident polar angles ((a) θ i =0 o , (b) θ i =30 o , (c) θ i =45 o , (d) θ i =60 o ) are compared with experimental results by Luntz et al. (white symbols).