**2. Simulation methods**

308 Molecular Dynamics – Theoretical Developments and Applications in Nanotechnology and Energy

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

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, 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

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., 1999, 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

method and the classical molecular dynamics method, is also being developed.

cannot be used to reproduce experimentally obtained data (Vincent et al., 2004).

calculation.

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

#### **2.1 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:

$$H\Psi\_o = E\_o\Psi\_o \tag{1}$$

where Ψ*o* and *Eo* 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

$$H = K + V\_{\text{ee}} + V\_{\text{ex}} \tag{2}$$

Molecular Simulation of Dissociation Phenomena of Gas Molecule on Metal Surface 311

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)

> *F*

 

**r** (9)

*d N* **r r** (10)

*v*

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

Let us consider a system that consists of *N* electrons without interaction between electrons.

*H KV v*

<sup>1</sup> det ! *o N N*

ex v 1 1

*i i*

 

The one-electron wave function, *φi*, corresponds to the function in which the eigenenergy is

 <sup>2</sup> v

<sup>2</sup> *i i i ii <sup>v</sup>* 

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

*E K vd* vv v

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:

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

1 2

 

**<sup>r</sup>** (13)

**r rr** (14)

1 2 *N N*

<sup>2</sup>

**<sup>r</sup>** (11)

(12)

*i i*

**r**

while considering the conservation of the number of electron as

and satisfies the following Euler equation:

The Hamiltonian of the system is expressed as

Slater determinant of the one-electron wave function, *φi*, as

lowest among all of the wave functions satisfying Eq. (13).

1

energy of the system is a functional of electron density, *ρ*(*r*), and is obtained by

field method.

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

$$\text{cov}\left(\mathbf{r}\_{i}\right) = \sum\_{k=1}^{M} \frac{-Z\_{k}}{\left|\mathbf{r}\_{i} - \mathbf{R}\_{k}\right|}\tag{3}$$

assuming that *v*(*r*) is a Coulomb interaction from the nuclei of the system. In Eq. (3), *ri* and *Rk* denotes the position of electron *i* and that of nucleus *k* reduced by the Bohr radius, respectively, and *Zk* 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

$$V\_{\rm ee} = \sum\_{i=1}^{N} \sum\_{j>i}^{N} \frac{1}{|\mathbf{r}\_i - \mathbf{r}\_j|} \tag{4}$$

The solution of Eq. (1) satisfies the equation below:

$$E\_o[\Psi\_o] = \underset{\Phi \to \Psi\_o}{\text{Min}} \{ \Phi | H | \Phi \}\tag{5}$$

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 *Eo* 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

$$\rho(\mathbf{r}) = \left| \Psi\_o \right|^2(\mathbf{r}) \tag{6}$$

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.

$$E[\rho] = \underset{\rho^\circ \to \rho}{\text{Min}} \left( K[\rho^\circ] + V\_{\text{ee}}[\rho^\circ] + V\_{\text{ex}}[\rho^\circ] \right) \tag{7}$$

In density functional theory, the electron density is obtained by Eq. (7). The total energy of the system, *E*[*ρ*], is expressed as

$$E[\rho] = F[\rho] + \int \rho(\mathbf{r}) v(\mathbf{r}) d\mathbf{r} \tag{8}$$

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:

$$
\mu = \upsilon(\mathbf{r}) + \frac{\delta F[\rho]}{\delta \rho(\mathbf{r})} \tag{9}
$$

while considering the conservation of the number of electron as

310 Molecular Dynamics – Theoretical Developments and Applications in Nanotechnology and Energy

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

1

assuming that *v*(*r*) is a Coulomb interaction from the nuclei of the system. In Eq. (3), *ri* and *Rk* denotes the position of electron *i* and that of nucleus *k* reduced by the Bohr radius, respectively, and *Zk* 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

*k i k Z*

*k*

*M*

ee

*V*

1

Min

external force field, *v*(*r*). Moreover, the electron density at position *r* is obtained by

'

 

 *E H o o* 

 

1 *N N i ji i j*

*o*

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 *Eo* of the system are obtained by Eq. (5). Namely, the state of the system is a functional of the interaction with the

<sup>2</sup>

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

ee ex

In density functional theory, the electron density is obtained by Eq. (7). The total energy of

*E F vd*

 

*E KV V* Min '''

**r r**

*v*

Coulomb interaction between electron *i* and *j* as

The solution of Eq. (1) satisfies the equation below:

minimize the total energy of the system.

the system, *E*[*ρ*], is expressed as

*i*

atomic units as

*H KV V* ee ex (2)

**<sup>r</sup> r R** (3)

(5)

*<sup>o</sup>* **r r** (6)

(7)

**r rr** (8)

(4)

$$\int \rho(\mathbf{r})d\mathbf{r} = N \tag{10}$$

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

$$H = K + V\_{\text{ex}} = \sum\_{i=1}^{N} \left( -\frac{1}{2} \nabla\_i^2 \right) + \sum\_{i=1}^{N} v\_{\text{v}} \left( \mathbf{r}\_i \right) \tag{11}$$

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

$$\Psi\_o = \frac{1}{\sqrt{N!}} \det \left[ \varphi\_1 \varphi\_2 \cdots \varphi\_N \right] \tag{12}$$

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).

$$\left[-\frac{1}{2}\nabla\_i^2 + v\_\mathbf{v}\left(\mathbf{r}\_i\right)\right]\rho\_i = \varepsilon\_i\rho\_i\tag{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

$$E\_{\mathbf{v}}[\rho] = K\_{\mathbf{v}}[\rho] + \int \rho(\mathbf{r}) v\_{\mathbf{v}}(\mathbf{r}) d\mathbf{r} \tag{14}$$

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:

$$\mu = v\_{\rm v}(\mathbf{r}) + \frac{\delta \mathcal{K}\_{\rm v}[\rho]}{\delta \rho(\mathbf{r})} \tag{15}$$

Molecular Simulation of Dissociation Phenomena of Gas Molecule on Metal Surface 313

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

in in

(5-22)

(5-13)

 outin

(5-6)

No

(5-20)

(*r*) and *v*xc(*r*) are functionals of electron

density by the convergence scheme because

density. A flowchart of this scheme is shown in Fig. 1.

in

the empirical exchange-correlation energy.

**2.2 Embedded Atom Method (EAM)** 

 initin

out *<sup>o</sup>*

 outin ?

> out

'

*<sup>r</sup> rr r rrr*

> *<sup>r</sup>* <sup>2</sup> 2 1

 *iiiivi v* 

 

 ' '

*d vvv v eff <sup>v</sup> xc*

> *<sup>o</sup> <sup>N</sup> <sup>N</sup>* det <sup>21</sup> !

> > 2

Yes

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

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 3*N* space including the interaction between electrons by including the interaction between electrons in

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

<sup>1</sup> (5-12)

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

$$F[\rho] = K\_{\rm v}[\rho] + f[\rho] + E\_{\rm xc}[\rho] \tag{16}$$

where *J*[*ρ*] denotes a classical Coulomb interaction between electrons and is expressed as

$$J[\rho] = \frac{1}{2} \int \int \frac{\rho(\mathbf{r})\rho(\mathbf{r}')}{|\mathbf{r} - \mathbf{r}'|} d\mathbf{r} d\mathbf{r}' \tag{17}$$

and the functional,

$$E\_{\infty}\left[\boldsymbol{\rho}\right] = F\left[\boldsymbol{\rho}\right] - K\_{\text{v}}\left[\boldsymbol{\rho}\right] - J\left[\boldsymbol{\rho}\right] = \left(K\left[\boldsymbol{\rho}\right] - K\_{\text{v}}\left[\boldsymbol{\rho}\right]\right) + \left(V\_{\text{ee}}\left[\boldsymbol{\rho}\right] - J\left[\boldsymbol{\rho}\right]\right) \tag{18}$$

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*[*ρ*]-*K*v[*ρ*]). The functional, *E*xc[*ρ*], is called the exchange-correlation energy. Using Eqs.(9) and (16), we obtain

$$\boldsymbol{\mu} = \boldsymbol{\upsilon}(\mathbf{r}) + \frac{\delta \mathcal{K}\_{\mathbf{v}}[\boldsymbol{\rho}]}{\delta \boldsymbol{\rho}(\mathbf{r})} + \frac{\delta \mathcal{J}[\boldsymbol{\rho}]}{\delta \boldsymbol{\rho}(\mathbf{r})} + \frac{\delta \mathcal{E}\_{\mathbf{xc}}[\boldsymbol{\rho}]}{\delta \boldsymbol{\rho}(\mathbf{r})} = \boldsymbol{\upsilon}(\mathbf{r}) + \boldsymbol{\phi}(\mathbf{r}) + \boldsymbol{\upsilon}\_{\mathbf{xc}}(\mathbf{r}) + \frac{\delta \mathcal{K}\_{\mathbf{v}}[\boldsymbol{\rho}]}{\delta \boldsymbol{\rho}(\mathbf{r})} \tag{19}$$

where

$$\phi(\mathbf{r}) = \int \frac{\rho(\mathbf{r}')}{|\mathbf{r} - \mathbf{r}'|} d\mathbf{r}' \tag{20}$$

denotes the Coulomb interaction of the electrons and

$$\upsilon\_{\text{xc}}\left(\mathbf{r}\right) = \frac{\delta E\_{\text{xc}}\left[\rho\right]}{\delta\rho\left(\mathbf{r}\right)}\tag{21}$$

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,

$$\upsilon\_{\rm eff} \left( \mathbf{r} \right) = \upsilon\_{\rm v} \left( \mathbf{r} \right) = \upsilon \left( \mathbf{r} \right) + \varphi \left( \mathbf{r} \right) + \upsilon\_{\rm uc} \left( \mathbf{r} \right) \tag{22}$$

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

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

v

**r**

*FK JE*

where *J*[*ρ*] denotes a classical Coulomb interaction between electrons and is expressed as

*E FK J KK V J* xc

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*[*ρ*]-*K*v[*ρ*]). The functional, *E*xc[*ρ*], is called the exchange-correlation energy. Using Eqs.(9) and (16), we obtain

 

' ' '

 xc

**<sup>r</sup> <sup>r</sup>**

be an interaction with the external force field in the virtual system, *v*v(*r*) in Eq. (15) (referred

*v vv v* eff v **r rrr r**

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

*E*

**r r r r r**

v xc v

**rrr r**

*d*

*K JE K*

 <sup>1</sup> ' ' 2 ' *J d <sup>d</sup>* 

**r r**

**r r**

 

 v x 

*v*

 

 

xc

denotes the exchange-correlation energy. Assuming the value in Eq. (19), *v*(*r*)+

*v*

 

*v v v*

 

**r r r r**

 

with the external force field of the virtual system, *v*v(*r*).

The functional in Eq. (9), *F*[*ρ*], is expressed as

 

denotes the Coulomb interaction of the electrons and

to as effective potential, *v*eff(*r*)), namely,

 

and the functional,

where

 v

*K*

 

**r r**

 

 

<sup>v</sup> <sup>v</sup> ee (18)

xc

<sup>c</sup> (16)

 

(17)

(19)

 

(20)

(21)

xc (22)

(*r*)+*v*xc(*r*), to

(15)

**r**

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.

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

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 3*N* space including the interaction between electrons by including the interaction between electrons in the empirical exchange-correlation energy.

#### **2.2 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

$$E = \sum\_{i=1}^{N} F\_i \left(\rho\_i\right) + \frac{1}{2} \sum\_{i=1}^{N} \sum\_{j \neq i} \phi\_{ij} \left(R\_{ij}\right) \tag{23}$$

Molecular Simulation of Dissociation Phenomena of Gas Molecule on Metal Surface 315

In this section the example of the analysis of dissociation phenomena of gas molecule on metal surface is introduced. We use the dissociation of H2 molecule on Pt(111) surface, which is important not only for the fundamental system of dissociation phenomena but also

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, DMol3. 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.010-6 Ha and the convergence tolerance for geometry was set at 1.010-5 Ha. Smearing was set at 5.010-3 Ha for both calculations to improve convergence of the calculations. Spin was unrestricted

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 22 = 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

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

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

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

, *θ*) in the space coordinate system. For the

avoid interactions normal to the surface by setting a long vacuum region.

LDA method is usually used for geometry optimization.

defining the Euler angle of the molecular axis (

for the reaction at anode side in Polymer Electrolyte Fuel Cell.

**3. Results and discussion** 

**3.1 Potential energy surface** 

during calculations.

where *Fi*() denotes the energy necessary to embed atom *i* into the background electron density *ρ*, and *ij*(*Rij*) denotes the core-core pair repulsion between atoms *i* and *j* separated by the distance *Rij*. In addition, *ρ<sup>i</sup>* 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

$$\rho\_i = \sum\_{j \neq i}^{N} \rho\_j^a \left( R\_{ij} \right) \tag{24}$$

where *ρ<sup>j</sup> <sup>a</sup>* is the electron density contributed by atom *j*. It is necessary to determine the functions, *Fi*(*ρ*), *ij*(*Rij*), and *ρ<sup>j</sup> <sup>a</sup>* (*Rij*) 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, *Fi*(*ρ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*(*Rij*), is often determined such that the crystal structure of the bulk system can be reproduced by the EAM using the determined values of *ρ<sup>j</sup> <sup>a</sup>* (*Rij*) and *Fi*(*ρi*). These functions are determined for every type of atom.

The force acting on atom *i* is obtained by differentiating Eq. (23), as follows:

$$F = -\frac{\partial E}{\partial \mathbf{R}\_i} = -\sum\_{i=1}^{N} \frac{dF\_i(\rho\_i)}{d\rho\_i} \frac{\partial \rho\_i}{\partial \mathbf{R}\_i} - \frac{1}{2} \sum\_{i=1}^{N} \sum\_{j \neq i} \frac{\partial \phi\_{ij}\{R\_{ij}\}}{\partial \mathbf{R}\_{ij}} \frac{\partial \mathbf{R}\_{ij}}{\partial \mathbf{R}\_i} \tag{25}$$

### **3. Results and discussion**

314 Molecular Dynamics – Theoretical Developments and Applications in Nanotechnology and Energy

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

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

 *R*

) denotes the energy necessary to embed atom *i* into the background electron

*ij*(*Rij*) denotes the core-core pair repulsion between atoms *i* and *j* separated

(23)

(24)

*ij*(*Rij*), is often determined such that

*<sup>a</sup>* (*Rij*) included in Eqs. (23) and (24) in order to use the

*i i ij ij*

1 2

by the distance *Rij*. In addition, *ρ<sup>i</sup>* 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

> *<sup>N</sup> <sup>a</sup> i j ij j i*

*<sup>a</sup>* is the electron density contributed by atom *j*. It is necessary to determine the

1 2 *N N ij ij ij i i <sup>i</sup>*

*d R*

*i ii i i ji ij i E dF R R*

*<sup>a</sup>* (*Rij*) and *Fi*(*ρi*). These functions are determined for every type

**RR R** (25)

 *R*

embedded atom method. An example of how to determine these functions or parameters is

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, *Fi*(*ρ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;

the crystal structure of the bulk system can be reproduced by the EAM using the

1 1

The force acting on atom *i* is obtained by differentiating Eq. (23), as follows:

1 1

*i i j i*

*N N*

*E F*

and Baskes, 1984).

expressed as

where *Fi*(

where *ρ<sup>j</sup>*

of atom.

functions, *Fi*(*ρ*),

described below.

determined values of *ρ<sup>j</sup>*

density *ρ*, and

density of each atom as

*ij*(*Rij*), and *ρ<sup>j</sup>*

Norskov, 1982). The core-core pair-repulsion,

*F*

In this section the example of the analysis of dissociation phenomena of gas molecule on metal surface is introduced. We use the dissociation of H2 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.

### **3.1 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, DMol3. 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.010-6 Ha and the convergence tolerance for geometry was set at 1.010-5 Ha. Smearing was set at 5.010-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 22 = 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

Molecular Simulation of Dissociation Phenomena of Gas Molecule on Metal Surface 317

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

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

> 2 tr tr tr tr tr <sup>1</sup> 2 exp *<sup>e</sup> <sup>f</sup> e de <sup>e</sup> de*

where *k* denotes the Boltzmann's constant. In this theory, all of the molecules that have a more energy than the dissociation barrier are considered to dissociate. This probability is

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

In the EAM, the potential energy, *E*pot, of the system, which consists of a Pt surface of *N* 

*i k ij ik*

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

2 2 pot Pt H Pt Pt Pt H H H 12

1 1 1 1 1

*N N N N*

*i k i ji i k*

(28)

 

*R RR*

 

2

<sup>1</sup> <sup>2</sup> exp *Eb <sup>e</sup> <sup>P</sup> <sup>e</sup> de kT kT*

*kT kT*

tr tr tr

(26)

(27)

surface is in the range of [*e*tr, *e*tr+*de*tr] is obtained by Boltzmann distribution as

obtained by integrating Eq. (26) in the range of [*Eb*, ∞], as follows:

**3.2 Energy transfer between gas molecule and solid atoms** 

the core-core pair-repulsion function is also described.

atoms and an H2 molecule, is obtained by

*EF F* 

barrier cannot dissociate.

(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.

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.

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

(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

site. In this calculation, *r* is changed at 0.5 Å intervals in the range between 0.5 Å and 2.5 Å,

1st layer

*x*

0.60 0.20

1.0

B

Fig. 3. Contour of the potential energy surface. The contour is plotted at 0.2 eV intervals. The

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

1.5

2.0

Z Distance [Å]

2.5

3.0

0.20 0.60



B

0.5 1.0 1.5 2.0 2.5

H-H Distance [Å]

2nd layer

3rd layer

=0° and *θ*=90° of brg

=0° and *θ*=90° of top site, and the right-hand side represents that at

and *Z* is changed at 0.5 Å intervals in the range between 0.5 Å and 2.5 Å.

brg site

Calculation domain for DFT

*y*

top site

Fig. 2. Simulation system for the potential energy surface.

H2

0.60 0.20

white dashed lines show the dissociation path.

1.0

1.5

2.0

Z Distance [Å]

2.5

3.0

0.20

0.60



0.5 1.0 1.5 2.0 2.5

H-H Distance [Å]



A A

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

$$f\left(e\_{\rm tr}\right)de\_{\rm tr} = 2\left(\frac{1}{kT}\right)^2 e\_{\rm tr} \exp\left[-\frac{e\_{\rm tr}}{kT}\right] de\_{\rm tr} \tag{26}$$

where *k* denotes the Boltzmann's constant. In this theory, all of the 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 [*Eb*, ∞], as follows:

$$P = 2\left(\frac{1}{kT}\right)^2 \int\_{E\_b}^{v\_0} e\_{\rm tr} \exp\left[-\frac{e\_{\rm tr}}{kT}\right] de\_{\rm tr} \tag{27}$$

#### **3.2 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 H2 molecule, is obtained by

$$E\_{\rm pot} = \sum\_{i=1}^{N} F\_{\rm Pt} \left( \rho\_i \right) + \sum\_{k=1}^{2} F\_{\rm H} \left( \rho\_k \right) + \sum\_{i=1}^{N} \sum\_{j>i}^{N} \phi\_{\rm Pt-Pt} \left( R\_{ij} \right) + \sum\_{i=1}^{N} \sum\_{k=1}^{2} \phi\_{\rm Pt-H} \left( R\_{ik} \right) + \phi\_{\rm H-H} \left( R\_{12} \right) \tag{28}$$

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

Molecular Simulation of Dissociation Phenomena of Gas Molecule on Metal Surface 319

performed until the system reached equilibrium at the target temperature. In this

H2 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 H2 molecule normal to the surface as the impinging energy. Neither rotational nor vibrational energy was given to the

A typical example of the difference from the initial energy of kinetic energy of Pt atoms, kinetic energy of the impinging H2 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 H2 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.1 eV, respectively, because an H2 molecule with a smaller impinging energy takes more time to reach the surface. Typical examples of the behavior of an H2 molecule on a Pt(111) surface upon collision are shown in Fig. 6. Fig. 6 (a) and (b) show the dissociation case and no dissociation case, respectively. Both the distance of the center of mass of the H2 molecule from the surface and the distance between H atoms in the H2 molecule are shown.

> Time [fs] 0 50 100 150 200 250 300

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 H2 molecule,

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*=5060 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 H2 molecule migrates on the surface. In the present simulations, H2 molecules having a distance between H atoms larger than 3.5 Ǻ

H2 molecule. The orientation of the impinging H2 molecule was given at random.

simulation, the target temperature was set at *T*=350 K.

Difference from initial energy [-]


respectively.



0

1

2

3

*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(*Rij*), H-H(*R12*) and Pt-H(*Rik*), 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 H2 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 10103 = 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.

Fig. 4. Schematic diagram of simulation system.

Before simulating the dissociation of H2 molecules, the relaxed Pt(111) surface, whose temperature was controlled to a target temperature, had to be obtained. The simulation was

*F*H(*ρi*) denote the energy necessary to embed a Pt atom or an H atom, respectively, into the

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

Next, an example of simulating the dissociation phenomena of H2 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 10103 = 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

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.

*z*

atoms)

*y*

rot *e*

tr *e*

27.72 [Å] (10 atoms)

Before simulating the dissociation of H2 molecules, the relaxed Pt(111) surface, whose temperature was controlled to a target temperature, had to be obtained. The simulation was

*x*

H2 Molecule

Pt surface

24.00 [Å] (10

Fig. 4. Schematic diagram of simulation system.

H-H(*R12*) and

Pt-H(*Rik*), denote the core-core

Pt-Pt(*Rij*),

background electron density, *ρi*, and

Ito, 2007; 2011).

Pt is set at 3.92 Å.

performed until the system reached equilibrium at the target temperature. In this simulation, the target temperature was set at *T*=350 K.

H2 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 H2 molecule normal to the surface as the impinging energy. Neither rotational nor vibrational energy was given to the H2 molecule. The orientation of the impinging H2 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 H2 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 H2 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.1 eV, respectively, because an H2 molecule with a smaller impinging energy takes more time to reach the surface. Typical examples of the behavior of an H2 molecule on a Pt(111) surface upon collision are shown in Fig. 6. Fig. 6 (a) and (b) show the dissociation case and no dissociation case, respectively. Both the distance of the center of mass of the H2 molecule from the surface and the distance between H atoms in the H2 molecule are shown.

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 H2 molecule, respectively.

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*=5060 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 H2 molecule migrates on the surface. In the present simulations, H2 molecules having a distance between H atoms larger than 3.5 Ǻ

Molecular Simulation of Dissociation Phenomena of Gas Molecule on Metal Surface 321

"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

> 0.0 0.2 0.4 0.6 0.8 1.0

Dissociation probability [-]

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

Initially, the probability that an impinging H2 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 H2 molecules clearly shows that H2 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

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

0.0 0.1 0.2 0.3 0.4 0.5

: brg : top

: fcc

Impinging energy [eV]

0.0 0.1 0.2 0.3 0.4 0.5

(a) (b)

Impinging energy [eV]

: brg : top

: fcc

in detail.

0.0 0.2 0.4 0.6 0.8 1.0

at the top site.

dissociation probability.

Dissociation probability [-]

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.

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

#### **3.3 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 H2 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 H2 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 H2 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

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

0

2

4

Distance [Å]

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

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 H2 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

By the way, the "static" dissociation probabilities, *P*s, can be obtained from the interaction potential between an H2 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 H2 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

6

8

0 50 100 150

Time [fs]

: Distance between atoms : Distance from surface

(b)

strong interaction to the surface.

0

dissociate.

Fig. 7 (a).

**3.3 Dissociation probability** 

2

4

Distance [Å]

6

8

0 50 100 150

Time [fs]

: Distance between atoms : Distance from surface (a) "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.

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 H2 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 H2 molecules clearly shows that H2 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

Molecular Simulation of Dissociation Phenomena of Gas Molecule on Metal Surface 323

according to Boltzmann distribution. The detailed method was described in a reference

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

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*<sup>d</sup> 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

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

The rise in the dissociation probability observed at low translational energies in the MD

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

results cannot be compared with the experiments due to lack of experimental data.

(Koido et al., 2011).

be neglected.

surface normal.

**4. Summary** 

molecular beam experiments to a certain degree.

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 H2 molecule at the chemisorption layer by steering effect, not the probability that the H2 molecule reaches the chemisorption layer. The dynamic effects on dissociation probability are very important when dissociation phenomena at the top site are considered.

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.

#### **3.4 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, D2 molecules were used as impinging gas molecules because in the experimental study with which we compared our results, D2 molecules were used. The EAM potential constructed for the H2 molecule was used for the D2 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 0o (normal to the surface) to 60o and the incident azimuthal angle was varied in a uniform random manner. Moreover, the initial position of H2 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*<sup>d</sup> 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.

#### **4. Summary**

322 Molecular Dynamics – Theoretical Developments and Applications in Nanotechnology and Energy

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 H2 molecule at the chemisorption layer by steering effect, not the probability that the H2 molecule reaches the chemisorption layer. The dynamic effects on dissociation probability are very important

Impinging Energy [eV]

0.0 0.1 0.2 0.3 0.4 0.5

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.

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, D2 molecules were used as impinging gas molecules because in the experimental study with which we compared our results, D2 molecules were used. The EAM potential constructed for the H2 molecule was used for the D2 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 0o (normal to the surface) to 60o and the incident azimuthal angle was varied in a uniform random manner. Moreover, the initial position of H2 molecule is given at random. The rotational energies were given

: *P*<sup>d</sup> : *P*<sup>c</sup>

when dissociation phenomena at the top site are considered.

Dissociation Probability [-]

0.0

**3.4 Comparison with experimental results** 

0.2

0.4

0.6

0.8

1.0

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

Molecular Simulation of Dissociation Phenomena of Gas Molecule on Metal Surface 325

Allen, M. P. and Tildesley, D. J. (1986). Computer Simulation of Liquids, Oxford University

Baskes, M. I. (1992). Modified embedded-atom potentials for cubic materials and impurities,

Baskes, M. I. ; Srinivasan, S. G. ; Valone, S. M. & Hoagland, R. G. (2007), Multistate modified

Beutl, M.; Riedler M. & Rendulic K. D. (1995). Strong Rotational Effects in the Adsorption

Blömer, J and Beylich, A. E. (1997). MD-Simulation of Inelastic Molecular Collisions with

Car, R. & Parrinello, M. (1985). Unified Approach for Molecular Dynamics and Density-

Clementi, E. & Roetti, C. (1974). Roothaan-Hartree-Fock Atomic Wavefunctions Basis

Darling, G. R. & Holloway, S. (1994). Rotational motion and the dissociation of H2 on Cu(111), Journal of Chemical Physics, Vol. 101, No. 4, pp.3268—3281 Darling, G. R.; Key, M. & Holloway, S. (1998). The steering of molecules in simple

Daw, M. S. & Baskes, M. I. (1983). Semiempirical, Quantum Mechanical Calculation of Hydrogen Embrittlement in Metals, *Physical Review Letters*, Vol. 50, pp.1285-1288. Daw, M. S. & Baskes, M.I. (1984). Embedded--atom method: Derivation and application to

Gross, A. ; Wilke, s. & Scheffler, M. (1995). Six-Dimensional Quantum Dynamics of

Gross, A. (1998). Hydrogen Dissociation on Metal Surface - A Model System for Reactions

Ishimoto, R.; Jung. C. H.; Tsuboi, H.; Koyama, M.; Endou, A.; Kubo, M.; Del Carpio, C. A. &

Luntz, A. C.; Brown, J. K. & Williams, M. D. (1990), Molecular beam studies of H2 and D2

impurities, surfaces, and other defects in metals, *Physical Review B*, Vol. 29, pp.6443-

Adsorption and Desorption of H2 at Pd(100): Steering and Steric Effects, Physical

Miyamoto, A. (2006). Periodic Density Functional and Tight-Binding Quantum Chemical Molecular Dynamics Study of Catalytic Properties on Gamma-Al2O3 Supported Pt Catalysts, *Applied Catalysis A-General*, Vol. 305, (2006), pp.64-69. Jinnouchi, R. & Okazaki, K. (2003), New Insight into Microscale Transport Phenomena in PEFC by Quantum MD, *Microscale Thermophysical Engineering*, Vol. 7, pp.15-31. Koido, T., Tomarikawa, K., Yonemura, S. & Tokumasu, T. (2011), Molecular Dynamics

Study of the Effects of Translational Energy and Incident Angle on Dissociation Probability of the H2/D2-Pt(111) System, Journal of Applied Physics, Vol. 110, No.

dissociative chemisorption on Pt(111), Journal of Chemical Physics, Vol. 93, No. 7,

Funtional Theory. *Physical Review Letters,* Vol. 55, pp.2471-2474.

dissociation reactions, *Surface Science*, Vol. 400, pp.314--328

on Surface, *Applied Physics A*, Vol. 67, No. 6, pp.627-635.

Dynamics of H2/Pd(111): Evidence for Dynamical Steering, *Chemical Physics Letters*,

Condensed Matter Surfaces, *Proceedings of Rarefied Gas Dynamics*, pp.392-397,

Functions and Their Coefficients for Ground and Certain Excited States of Neutral and Ionization Atoms, Z < 54, *Atomic Data and Nuclear Data Tables*, Vol. 14, pp.177-

embedded atom method, *Physical Review B*, Vol. 75, 094113.

**6. References** 

478.

6453.

2, pp. 024301

pp.5240—5246

Press, (1989), pp.80-81.

Vol. 247, pp.249-252.

*Physical Review B*, Vol. 46, pp.2727-2742.

Beijing, China, August 19-23, 1996.

Review Letters, Vol. 75, pp.2718—2721

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.

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=0o, (b) *θ*i=30o, (c) *θ*i=45o, (d) *θ*i=60o) are compared with experimental results by Luntz et al. (white symbols).

#### **5. 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.

#### **6. References**

324 Molecular Dynamics – Theoretical Developments and Applications in Nanotechnology and Energy

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.

0.0 0.1 0.2 0.3 0.4 0.5

Translational Energy, *E*tr

0.0 0.1 0.2 0.3 0.4 0.5

Translational Energy, *E*tr

system of Institute of Fluid Science, Tohoku University.

(a) *θ*i = 0o (b) *θ*i = 30o

[eV]

[eV]

(c) *θ*i = 45o (d) *θ*i = 60o

0.0

0.0

0.2 0.4

Dissociation Probability,

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=0o, (b) *θ*i=30o, (c) *θ*i=45o, (d) *θ*i=60o) are compared with experimental

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

*P*d [-]

0.6

0.8 1.0

 MD Experiment

0.2

0.4

Dissociation Probability,

*P*d [-]

0.6

0.8

 MD Experiment

1.0

0.0 0.1 0.2 0.3 0.4 0.5

Translational Energy, *E*tr

0.0 0.1 0.2 0.3 0.4 0.5

Translational Energy, *E*tr

[eV]

[eV]

0.0

0.0

**5. Acknowledgment**

results by Luntz et al. (white symbols).

0.2

0.4

Dissociation Probability,

*P*d [-]

0.6

0.8

 MD Experiment

1.0

0.2

0.4

Dissociation Probability,

*P*d [-]

0.6

0.8

 MD Experiment

1.0


**1. Introduction**

It has attracted attention for decades already.

In three dimensions, polymer dynamics exhibits a rich and complex behavior which depends on the solvent conditions and polymer concentration (1; 2). That the dynamics of polymer chains at and near solid interfaces differs profoundly from that in the bulk is intuitively expected. Polymer adsorption on the surface is of technological and scientific importance in the field of colloids and biomolecules. Examples include the two-dimensional (2-D) diffusion of DNA oligonucleotides confined to biological interfaces such as cell membranes (3; 4). The diffusion of confined polymers at surface is always a fundamental, yet problematical topic in polymer physics (1; 2; 5–8). The behavior of polymers at the liquid-solid interface is crucial to technologies involving molecular surface placement (9; 10). Polymers adsorbed onto a surface to form thin films is an emerging topic of modern materials science (11; 12). They can be applied, for example, in the fields of biosensors, light-emitting diodes, nonlinear optical devices, and permeation-selective gas membranes (13–18). The fabrication of the thin films always takes place in vacuum or dilute solutions. The adsorption of the polymer chains can be controlled by varying multiple parameters such as the polymer-surface interaction, the solvent quality, surface nano-roughness, temperature, polymer chain length and so on. The knowledge about the adsorption dynamics and the thermodynamics of the equilibrium adsorption is crucial to understand and furthermore improve the property of the final product. Experimentally it is difficult to control well the above influences separately, thus, clear information on the exact effect of a specific environmental parameter is hard to obtain. Due to the fast-growing computation power nowadays, it is possible to utilize the computer as an "experimental apparatus" to solve the problem, with the aid of various simulation techniques.

**A Study of the Adsorption and Diffusion** 

**Chain on a Silicon Surface by Molecular** 

*1College of Chemistry Chemical Engineering and Materials Science,* 

*2Photoelectric Engineering College, Zaozhuang University, Shandong* 

**Dynamics Simulation** 

Dan Mu1 and Jian-Quan Li2

*Zaozhuang University, Shandong* 

*China* 

**16**

**Behavior of a Single Polydimethylsiloxane** 

Granick and co-workers studied poly(ethylene glycol) molecules adsorbed on solid surface by means of fluorescence microscopy (19–22). They found that the diffusion coefficient (*D*) of such chains scales with the degree of polymerization (*N*) as *N*−3/2, which is characteristic for

