- Email: [email protected]
Hybrid simulation between molecular dynamics and binary collision approximation codes for hydrogen injection into carbon materials
Journal of Nuclear Materials 415 (2011) S208–S211
Contents lists available at ScienceDirect
Journal of Nuclear Materials journal homepage: www.elsevier.com/locate/jnucmat
Hybrid simulation between molecular dynamics and binary collision approximation codes for hydrogen injection into carbon materials Seiki Saito a,⇑, Atsushi M. Ito b, Arimichi Takayama b, Takahiro Kenmotsu c, Hiroaki Nakamura a,b a
Department of Energy Engineering and Science, Graduate School of Engineering, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8603, Japan Fundamental Physics Simulation Research Division, National Institute for Fusion Science, Toki, Gifu 509-5292, Japan c Faculty of Life and Medical Sciences, Doshisha University, Kyotanabe 610-0394, Japan b
a r t i c l e
i n f o
Article history: Available online 30 December 2010
a b s t r a c t Molecular dynamics (MD) simulation with modified Brenner’s reactive empirical bond order (REBO) potential is a powerful tool to investigate plasma wall interaction on divertor plates in a nuclear fusion device. However, the size of MD simulation box is generally set less than several nm because of the limits of a computer performance. To extend the size of the MD simulation, we develop a hybrid simulation code between MD code using REBO potential and binary collision approximation (BCA) code. Using the BCA code instead of computing all particles with a high kinetic energy for every step in the MD simulation, considerable computation time is saved. By demonstrating a hydrogen atom injection into a graphite by the hybrid simulation code, it is found that the hybrid simulation code works efficiently in a large simulation box. 2010 Elsevier B.V. All rights reserved.
1. Introduction
2. AC"T-MD hybrid simulation
To understand a mechanism of chemical and physical interactions between hydrogen plasmas and divertor plates in a nuclear fusion device, it is necessary to study elementary processes of the reactions. In molecular dynamics (MD) simulations, the equations of motion of atoms are solved numerically. We investigated plasma–wall interaction on the divertor plates made of carbon materials by the MD simulation with modified Brenner’s reactive empirical bond order (REBO) potential [1,2] in the previous works [3–5]. A typical scale length of the MD simulation box in these works is several nm. In order to expand the simulation box to more realistic scale length, i.e., several lm, we develop a hybrid simulation code between the MD simulation code with the REBO potential and atomic collision in any structured target (AC"T) code [6]. We note here that AC"T is an extended code of atomic collision in amorphous target (ACAT) code [7]. While ACAT code treats only amorphous structures, AC"T code treats arbitrary structures including amorphous and crystalline structure. In the AC"T code, where binary collision approximation (BCA) is used, a two-body interaction is calculated instead of computing all particles for every step in the MD simulation. Therefore computation time is saved. As the partner of the MD simulation, we choose the AC"T code for its fast calculation of a location and a velocity of the particle.
To explain the hybrid simulation more concretely, we consider a hydrogen atom injection into a graphite (Fig. 1). The injection kinetic energy of the hydrogen atom is set to higher kinetic energy than Eth at which a simulation algorithm changes from the AC"T to the MD simulation. While the kinetic energy of the hydrogen atom is higher than Eth, trajectories of the hydrogen atom and the surrounding carbon atoms are calculated by the AC"T simulation. The hydrogen atom dissipates its kinetic energy by interacting with carbon atoms in graphite and then the kinetic energy becomes lower than Eth. At that moment, the MD simulation starts to calculate the motions of the hydrogen atom and the surrounding carbon atoms in a box which is called MD simulation box. The main problem of hybridization between the two simulation codes is how to define a threshold kinetic energy Eth. To determine the threshold kinetic energy in our simulation, incident angle dependence of reflection rate in the case of hydrogen injection into amorphous carbon materials was investigated by both ACAT and MD simulation [8]. Comparing the MD with ACAT simulation, two simulation results coincide in higher kinetic energy than 200 eV. Thus we set Eth = 200 eV. Different from AC"T simulation, our MD simulation code cannot treat inelastic interactions of projectile ions. However, inelastic interactions are negligible in the energy range of our simulation. In the collision process of the AC"T simulation, the motion of carbon atoms which collide with the incident hydrogen depends on their received kinetic energy Er. When Er is lower than the
⇑ Corresponding author. E-mail address: [email protected] (S. Saito). 0022-3115/$ - see front matter 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.jnucmat.2010.12.233
S209
S. Saito et al. / Journal of Nuclear Materials 415 (2011) S208–S211
hydrogen AC
initialization signal
T
initialization signal
BCA calculation
MD domain
1005
MD initial material data
MD
MD calculation
MD 34 347
Fig. 2. Data flow diagram between the processes of the AC"T-MD hybrid simulation code.
321 Fig. 1. Schematic diagram of the AC"T-MD hybrid simulation of hydrogen injection into a graphite.
2.2. Molecular dynamics simulation The equation of motion in the MD simulation is written by
binding energy of covalent bond between carbon atoms, the collision carbon atom cannot move. On the other hand, when Er is larger than the binding energy, the collided carbon atom recoils. The motion of recoil carbon atoms is solved in the same way that the motion of the incident hydrogen is solved in AC"T and MD simulation. In order to realize the above hybrid simulation, we use Multiple Program Multiple Data (MPMD) [9] which was developed to execute different programs parallely in one simulation with synchronizing their data. Using the MPMD, we control four programs, i.e., (i) master program, (ii) AC"T program, (iii) MD program and (iv) splitting program as shown in Fig. 2. The master process is installed between the AC"T and MD processes for data synchronization. After an initialization of each process (Fig. 2 (1) and (2)), the AC"T process (Fig. 2 (3)) is executed until the kinetic energy of the incident hydrogen and all recoil atoms become lower than Eth. When the kinetic energies of all atoms become lower than Eth, the AC"T process sends, to the master process, a list of the positions and velocities of the moving atoms as a result of the AC"T simulation. The master process (Fig. 2 (4)) informs the splitter process of the completion of the AC"T process with the list of moved atoms. The splitter process updates its material data with the list of moved atoms. Then, the MD simulation boxes are picked up from the entire target material and they are sent to the MD processes which are running parallely (Fig. 2 (5)). The MD simulation is started in each MD process (Fig. 2 (6)). We will explain the algorithm of the AC"T and the MD simulation in the following subsections.
r_ i ¼
pi ; mi
p_ i ¼
@U ; @ri
ð3Þ
where ri, pi and mi are the position, momentum and the mass of the ith atom. The function U is the total potential energy. In our MD simulation, the modified REBO potential is used for the potential U,
U
X
U ij ;
ðfrg; fhB g; fhDH gÞV A ðr Þ; U ij V R½ij ðrij Þ b ij ½ij ij
ð4Þ
i;j>i
where rij is the distance between the ith and the jth atoms. The functions V R½ij and V A½ij represent repulsion and attraction, respec generates a multi-body force. The bond angle tively. The function b ij hBjik is the angle between the vector from the ith atom to the jth atom and the vector from the ith atom to the kth atom. The dihedral angle hDH kijl is the angle between the plane passing through the kth, ith, and jth atoms and the plane passing through the ith, jth, and lth atoms. Moreover, the parameters {r}, {hB}, and {hDH} denote all sets of rij ; hBjik ; and hDH kijl , respectively (for details of the modified Brenner REBO potential, see Refs. [10,11]). The second-order symplectic integration [12] is used to execute the time integration of the equation of motion where the time step is set to 5 1018 s. Brenner’s REBO potential was developed for hydrocarbons that can model intramolecular chemical bonding in a variety of small hydrocarbon molecules [1], so that the REBO potential is accurate for low potential energy. To the contrary, Thomas-Fermi potential is accurate for high potential energy where strong repulsive force is dominant. In our hybrid simulation, the REBO and the Thomas-Fermi potentials are adopted for low and high potential energy, respectively. Therefore, our hybrid simulation has an advantage over the pure MD and BCA simulation.
2.1. AC"T simulation
CII
A binary collision between a projectile and a target atom is necessary in the AC"T simulation [6]. The final position and velocity of the projectile and the target atom are obtained analytically in a potential V(r), where r is the distance between the projectile and the target atom. In the AC"T code similarly to the ACAT code, the Moliere approximation is adopted to the ThomasFermi potential V(r) as follows:
Z 1 Z 2 e2 Uðr=aÞ; r 0:3x UðxÞ 0:35e þ 0:55e1:2x þ 0:10e6:0x ; VðrÞ ¼
CI
ð1Þ
14
ð2Þ
where a is a screening length, Z1 and Z2 are atomic numbers of the projectile and the target atom, respectively. Unlike the MD simulation, the AC"T simulation directly gives us the asymptotic trajectory of the projectile collided with the target.
34 Fig. 3. The structure of a MD simulation box. The MD simulation box is divided into two regions, CI and CII. The motion of atoms in CI is simulated. The bonding states of atoms in CII are necessary to obtain the precise forces acting on all atoms in CI.
S210
S. Saito et al. / Journal of Nuclear Materials 415 (2011) S208–S211
Unlike the AC"T simulation, the MD simulation cannot treat a large scale material. To save the simulation time, cubic simulation boxes whose side is 14 Å long are picked up in the MD simulation, and the simulation for the atoms in each box is executed in parallel. This parallelization reduces the simulation time. Fig. 3 shows the structure of an MD simulation box. To obtain precise forces acting on all atoms in CI, the bonding states of atoms in CII are necessary (see Eq. (4)). Namely, in the MD process, the equations of motions of atoms in CI are solved, and forces acting on the atoms are obtained by using surrounding atoms in CII. We comment here a special case that two or more cubic boxes (CII) overlap. This case corresponds to the case that the final positions of moving atoms in AC"T are close. These overlapping boxes are simulated in one CPU process.
3. Numerical demonstration We demonstrate an AC"T-MD hybrid simulation for a hydrogen atom injection into an ideal graphite which is 321 Å long, 347 Å wide and 1005 Å deep. The graphite consists of 300 graphene sheets. Periodic boundary conditions are imposed on the horizontal direction. The initial temperature of the graphite is set to 0 K.
(a)
The hydrogen atom with the kinetic energy 1 keV is injected vertically into the graphite. Thousand simulations with the same initial graphite and randomly changed injection position of hydrogen have been performed. To compare the trajectories, pure AC"T simulations have also been performed. Fig. 4 shows the histograms of the depth d\ and the horizontal distance dk from the injection point at the surface to the final positions in the result of both the AC"T-MD hybrid and pure AC"T simulation. The histogram of the depth d\ by AC"T-MD has lower average and standard deviation than those by pure AC"T. This behavior is caused by the difference of the interaction potentials between AC"T and MD. In AC"T, the twobody interaction potential is used to calculate the force between particles. On the other hand, the multi-body interaction is treated in the MD simulation. Therefore, the kinetic energy of the moving particle is dissipated to the surrounding particle more easily in the MD simulation than in the AC"T simulation. This fact derives the shorter average of depth d\ in AC"T-MD. A comparison between the averages of dk in both simulations shows an agreement except one difference, that is, the standard
(a)
1005
AC T-MD average.: 384.5 std. deviation: 140.8
0.05
347
relative frequency
0.06
0.04
AC T average: 394.9 std. deviation: 156.7
0.03 0.02
(b)
0.01 0.00
321
0
200
400
600
800
A
(b) AC T-MD average.: 206.5 std. deviation: 90.0
1005
depth d [ ]
347
long tail
relative frequency
0.06
321
0.05
(c)
0.04 AC T average: 206.9 std. deviation: 84.1
0.03 0.02 0.01 0.00
0
200
400
600
800
horizontal distance d|| [ ] Fig. 4. Histograms of the depth d\ and the horizontal distance dk of 902 injections which were retained in the simulation box as the result of both the pure AC"T and AC"T-MD hybrid simulation. Ninety eight trajectories went out of the simulation box.
Fig. 5. Thirty sample trajectories of incident hydrogen that move in graphite calculated by the (a) pure AC"T code and (b) AC"T-MD hybrid simulation, respectively. (c) Detailed image of the MD simulation part for trajectory A in (b).
S. Saito et al. / Journal of Nuclear Materials 415 (2011) S208–S211
deviation by AC"T-MD is larger than that by pure AC"T. This difference is driven by the fact that the histogram by AC"T-MD has a long tail. The long tail appears in the histogram of dk by the AC"T-MD simulation, because hydrogen atoms horizontally move for a long distance in the interlayer region of the graphite. Fig. 5a and b shows 30 samples chosen randomly in 1000 trajectories simulated by the AC"T and the AC"T-MD hybrid code, respectively. Fig. 5b shows the trajectories in the same initial condition calculated by the AC"T-MD hybrid simulation. The gray and black lines in Fig. 5b are obtained by the AC"T and the MD simulation, respectively. The portion of the trajectories of AC"T (gray lines) in Fig. 5b are exactly the same as the trajectories at kinetic energies above 200 eV in Fig. 5a. Comparing these two figures, it is clearly found that the trajectories obtained by the MD simulation in AC"T-MD tend to move horizontally for a long distance in the interlayer region of the graphite. In our previous research [5], it was shown that there is an adsorption site 1.1Å above each carbon atom of a graphene sheet. As shown in Fig. 5c, the hydrogen atom is finally trapped at the adsorption site. This phenomenon suggests that our hybrid simulation treats the chemical processes more realistically than pure AC"T simulations. Thus, it is concluded that our hybrid simulation is useful to investigate the dynamical processes where chemical reactions take place in lm scale (ex., chemical sputtering and hydrogen blister generation). 4. Summary The AC"T-MD hybrid simulation code has been developed to treat large materials with wide range of kinetic energy in plasma surface interaction. The motion of hydrogen and all recoil atoms are solved by BCA method using the AC"T code while its kinetic energy is higher than Eth. When the kinetic energy becomes lower than Eth, the MD simulation starts to calculate the motions of the hydrogen and all recoil atoms with their surrounding carbon atoms in a local box.
S211
The AC"T-MD hybrid simulation was performed and compared with pure AC"T simulations of 1000 hydrogen atom injections into a graphite. Their results suggested that the AC"T-MD hybrid simulation tends to derive a smaller value of depth than the AC"T simulation does. This fact is caused by the interaction between the surrounding atoms, which is ignored in the AC"T simulation. Acknowledgements This work is supported by the National Institutes of Natural Sciences’ undertaking Forming Bases for Interdisciplinary and International Research through Cooperation Across Fields of Study and Collaborative Research Programs (NIFS09KEIN0091 and NIFS10KEIN0161) and a Grants-in-Aid for Scientific Research (No. 19055005) from the Ministry of Education, Culture, Sports, Science and Technology, Japan. This work is performed with the support and under the auspices of the NIFS Collaboration Research Programs (NIFS09KTAN008 and NIFS10KTAS005). References [1] D.W. Brenner, O.A. Shenderova, J.A. Harrison, S.J. Stuart, B. Ni, S.B. Sinnott, J. Phys.: Condens. Matter 14 (2002) 783. [2] A. Ito, H. Nakamura, Commun. Comput. Phys. 4 (2008) 592. [3] A. Ito, H. Nakamura, J. Plasma Phys. 72 (2006) 805. [4] A. Ito, H. Nakamura, A. Takayama, J. Phys. Soc. Jpn. 77 (2008) 114602. [5] S. Saito, A. Ito and H. Nakamura, Plasma Fusion Res., special issue, 5 (2010) S2076. [6] A. Takayama, S. Saito, A.M. Ito, T. Kenmotsu and H. Nakamura, Jpn. J. Appl. Phys., in press. [arXiv:1005.3344]. [7] Y. Yamamura, Y. Mizuno, IPPJ-AM-40, Inst. Plasma Phys, Nagoya Univ., 1985. [8] A. Ito, T. Kenmotsu, Y. Kikuhara, K. Ohya, K. Inai, Y. Wang, S. Irle, K. Morokuma, H. Nakamura, The 4th IAEA-TM on the Theory of Plasma Instabilities, Kyoto University, Japan, 2009. May 18-19. [9] A. Takayama, K. Shimizu, Y. Tomita, T. Takizuka, J. Plasma Fusion Res. SERIES 9 (2010) 604. [10] A. Ito, H. Nakamura, Jpn. J. Appl. Phys. 47 (2008) 4715. [11] A. Ito, H. Nakamura, A. Takayama, J. Phys. Soc. Jpn. 77 (2008) 114602. [12] M. Suzuki, J. Math. Phys. 26 (1985) 601.
Contents lists available at ScienceDirect
Journal of Nuclear Materials journal homepage: www.elsevier.com/locate/jnucmat
Hybrid simulation between molecular dynamics and binary collision approximation codes for hydrogen injection into carbon materials Seiki Saito a,⇑, Atsushi M. Ito b, Arimichi Takayama b, Takahiro Kenmotsu c, Hiroaki Nakamura a,b a
Department of Energy Engineering and Science, Graduate School of Engineering, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8603, Japan Fundamental Physics Simulation Research Division, National Institute for Fusion Science, Toki, Gifu 509-5292, Japan c Faculty of Life and Medical Sciences, Doshisha University, Kyotanabe 610-0394, Japan b
a r t i c l e
i n f o
Article history: Available online 30 December 2010
a b s t r a c t Molecular dynamics (MD) simulation with modified Brenner’s reactive empirical bond order (REBO) potential is a powerful tool to investigate plasma wall interaction on divertor plates in a nuclear fusion device. However, the size of MD simulation box is generally set less than several nm because of the limits of a computer performance. To extend the size of the MD simulation, we develop a hybrid simulation code between MD code using REBO potential and binary collision approximation (BCA) code. Using the BCA code instead of computing all particles with a high kinetic energy for every step in the MD simulation, considerable computation time is saved. By demonstrating a hydrogen atom injection into a graphite by the hybrid simulation code, it is found that the hybrid simulation code works efficiently in a large simulation box. 2010 Elsevier B.V. All rights reserved.
1. Introduction
2. AC"T-MD hybrid simulation
To understand a mechanism of chemical and physical interactions between hydrogen plasmas and divertor plates in a nuclear fusion device, it is necessary to study elementary processes of the reactions. In molecular dynamics (MD) simulations, the equations of motion of atoms are solved numerically. We investigated plasma–wall interaction on the divertor plates made of carbon materials by the MD simulation with modified Brenner’s reactive empirical bond order (REBO) potential [1,2] in the previous works [3–5]. A typical scale length of the MD simulation box in these works is several nm. In order to expand the simulation box to more realistic scale length, i.e., several lm, we develop a hybrid simulation code between the MD simulation code with the REBO potential and atomic collision in any structured target (AC"T) code [6]. We note here that AC"T is an extended code of atomic collision in amorphous target (ACAT) code [7]. While ACAT code treats only amorphous structures, AC"T code treats arbitrary structures including amorphous and crystalline structure. In the AC"T code, where binary collision approximation (BCA) is used, a two-body interaction is calculated instead of computing all particles for every step in the MD simulation. Therefore computation time is saved. As the partner of the MD simulation, we choose the AC"T code for its fast calculation of a location and a velocity of the particle.
To explain the hybrid simulation more concretely, we consider a hydrogen atom injection into a graphite (Fig. 1). The injection kinetic energy of the hydrogen atom is set to higher kinetic energy than Eth at which a simulation algorithm changes from the AC"T to the MD simulation. While the kinetic energy of the hydrogen atom is higher than Eth, trajectories of the hydrogen atom and the surrounding carbon atoms are calculated by the AC"T simulation. The hydrogen atom dissipates its kinetic energy by interacting with carbon atoms in graphite and then the kinetic energy becomes lower than Eth. At that moment, the MD simulation starts to calculate the motions of the hydrogen atom and the surrounding carbon atoms in a box which is called MD simulation box. The main problem of hybridization between the two simulation codes is how to define a threshold kinetic energy Eth. To determine the threshold kinetic energy in our simulation, incident angle dependence of reflection rate in the case of hydrogen injection into amorphous carbon materials was investigated by both ACAT and MD simulation [8]. Comparing the MD with ACAT simulation, two simulation results coincide in higher kinetic energy than 200 eV. Thus we set Eth = 200 eV. Different from AC"T simulation, our MD simulation code cannot treat inelastic interactions of projectile ions. However, inelastic interactions are negligible in the energy range of our simulation. In the collision process of the AC"T simulation, the motion of carbon atoms which collide with the incident hydrogen depends on their received kinetic energy Er. When Er is lower than the
⇑ Corresponding author. E-mail address: [email protected] (S. Saito). 0022-3115/$ - see front matter 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.jnucmat.2010.12.233
S209
S. Saito et al. / Journal of Nuclear Materials 415 (2011) S208–S211
hydrogen AC
initialization signal
T
initialization signal
BCA calculation
MD domain
1005
MD initial material data
MD
MD calculation
MD 34 347
Fig. 2. Data flow diagram between the processes of the AC"T-MD hybrid simulation code.
321 Fig. 1. Schematic diagram of the AC"T-MD hybrid simulation of hydrogen injection into a graphite.
2.2. Molecular dynamics simulation The equation of motion in the MD simulation is written by
binding energy of covalent bond between carbon atoms, the collision carbon atom cannot move. On the other hand, when Er is larger than the binding energy, the collided carbon atom recoils. The motion of recoil carbon atoms is solved in the same way that the motion of the incident hydrogen is solved in AC"T and MD simulation. In order to realize the above hybrid simulation, we use Multiple Program Multiple Data (MPMD) [9] which was developed to execute different programs parallely in one simulation with synchronizing their data. Using the MPMD, we control four programs, i.e., (i) master program, (ii) AC"T program, (iii) MD program and (iv) splitting program as shown in Fig. 2. The master process is installed between the AC"T and MD processes for data synchronization. After an initialization of each process (Fig. 2 (1) and (2)), the AC"T process (Fig. 2 (3)) is executed until the kinetic energy of the incident hydrogen and all recoil atoms become lower than Eth. When the kinetic energies of all atoms become lower than Eth, the AC"T process sends, to the master process, a list of the positions and velocities of the moving atoms as a result of the AC"T simulation. The master process (Fig. 2 (4)) informs the splitter process of the completion of the AC"T process with the list of moved atoms. The splitter process updates its material data with the list of moved atoms. Then, the MD simulation boxes are picked up from the entire target material and they are sent to the MD processes which are running parallely (Fig. 2 (5)). The MD simulation is started in each MD process (Fig. 2 (6)). We will explain the algorithm of the AC"T and the MD simulation in the following subsections.
r_ i ¼
pi ; mi
p_ i ¼
@U ; @ri
ð3Þ
where ri, pi and mi are the position, momentum and the mass of the ith atom. The function U is the total potential energy. In our MD simulation, the modified REBO potential is used for the potential U,
U
X
U ij ;
ðfrg; fhB g; fhDH gÞV A ðr Þ; U ij V R½ij ðrij Þ b ij ½ij ij
ð4Þ
i;j>i
where rij is the distance between the ith and the jth atoms. The functions V R½ij and V A½ij represent repulsion and attraction, respec generates a multi-body force. The bond angle tively. The function b ij hBjik is the angle between the vector from the ith atom to the jth atom and the vector from the ith atom to the kth atom. The dihedral angle hDH kijl is the angle between the plane passing through the kth, ith, and jth atoms and the plane passing through the ith, jth, and lth atoms. Moreover, the parameters {r}, {hB}, and {hDH} denote all sets of rij ; hBjik ; and hDH kijl , respectively (for details of the modified Brenner REBO potential, see Refs. [10,11]). The second-order symplectic integration [12] is used to execute the time integration of the equation of motion where the time step is set to 5 1018 s. Brenner’s REBO potential was developed for hydrocarbons that can model intramolecular chemical bonding in a variety of small hydrocarbon molecules [1], so that the REBO potential is accurate for low potential energy. To the contrary, Thomas-Fermi potential is accurate for high potential energy where strong repulsive force is dominant. In our hybrid simulation, the REBO and the Thomas-Fermi potentials are adopted for low and high potential energy, respectively. Therefore, our hybrid simulation has an advantage over the pure MD and BCA simulation.
2.1. AC"T simulation
CII
A binary collision between a projectile and a target atom is necessary in the AC"T simulation [6]. The final position and velocity of the projectile and the target atom are obtained analytically in a potential V(r), where r is the distance between the projectile and the target atom. In the AC"T code similarly to the ACAT code, the Moliere approximation is adopted to the ThomasFermi potential V(r) as follows:
Z 1 Z 2 e2 Uðr=aÞ; r 0:3x UðxÞ 0:35e þ 0:55e1:2x þ 0:10e6:0x ; VðrÞ ¼
CI
ð1Þ
14
ð2Þ
where a is a screening length, Z1 and Z2 are atomic numbers of the projectile and the target atom, respectively. Unlike the MD simulation, the AC"T simulation directly gives us the asymptotic trajectory of the projectile collided with the target.
34 Fig. 3. The structure of a MD simulation box. The MD simulation box is divided into two regions, CI and CII. The motion of atoms in CI is simulated. The bonding states of atoms in CII are necessary to obtain the precise forces acting on all atoms in CI.
S210
S. Saito et al. / Journal of Nuclear Materials 415 (2011) S208–S211
Unlike the AC"T simulation, the MD simulation cannot treat a large scale material. To save the simulation time, cubic simulation boxes whose side is 14 Å long are picked up in the MD simulation, and the simulation for the atoms in each box is executed in parallel. This parallelization reduces the simulation time. Fig. 3 shows the structure of an MD simulation box. To obtain precise forces acting on all atoms in CI, the bonding states of atoms in CII are necessary (see Eq. (4)). Namely, in the MD process, the equations of motions of atoms in CI are solved, and forces acting on the atoms are obtained by using surrounding atoms in CII. We comment here a special case that two or more cubic boxes (CII) overlap. This case corresponds to the case that the final positions of moving atoms in AC"T are close. These overlapping boxes are simulated in one CPU process.
3. Numerical demonstration We demonstrate an AC"T-MD hybrid simulation for a hydrogen atom injection into an ideal graphite which is 321 Å long, 347 Å wide and 1005 Å deep. The graphite consists of 300 graphene sheets. Periodic boundary conditions are imposed on the horizontal direction. The initial temperature of the graphite is set to 0 K.
(a)
The hydrogen atom with the kinetic energy 1 keV is injected vertically into the graphite. Thousand simulations with the same initial graphite and randomly changed injection position of hydrogen have been performed. To compare the trajectories, pure AC"T simulations have also been performed. Fig. 4 shows the histograms of the depth d\ and the horizontal distance dk from the injection point at the surface to the final positions in the result of both the AC"T-MD hybrid and pure AC"T simulation. The histogram of the depth d\ by AC"T-MD has lower average and standard deviation than those by pure AC"T. This behavior is caused by the difference of the interaction potentials between AC"T and MD. In AC"T, the twobody interaction potential is used to calculate the force between particles. On the other hand, the multi-body interaction is treated in the MD simulation. Therefore, the kinetic energy of the moving particle is dissipated to the surrounding particle more easily in the MD simulation than in the AC"T simulation. This fact derives the shorter average of depth d\ in AC"T-MD. A comparison between the averages of dk in both simulations shows an agreement except one difference, that is, the standard
(a)
1005
AC T-MD average.: 384.5 std. deviation: 140.8
0.05
347
relative frequency
0.06
0.04
AC T average: 394.9 std. deviation: 156.7
0.03 0.02
(b)
0.01 0.00
321
0
200
400
600
800
A
(b) AC T-MD average.: 206.5 std. deviation: 90.0
1005
depth d [ ]
347
long tail
relative frequency
0.06
321
0.05
(c)
0.04 AC T average: 206.9 std. deviation: 84.1
0.03 0.02 0.01 0.00
0
200
400
600
800
horizontal distance d|| [ ] Fig. 4. Histograms of the depth d\ and the horizontal distance dk of 902 injections which were retained in the simulation box as the result of both the pure AC"T and AC"T-MD hybrid simulation. Ninety eight trajectories went out of the simulation box.
Fig. 5. Thirty sample trajectories of incident hydrogen that move in graphite calculated by the (a) pure AC"T code and (b) AC"T-MD hybrid simulation, respectively. (c) Detailed image of the MD simulation part for trajectory A in (b).
S. Saito et al. / Journal of Nuclear Materials 415 (2011) S208–S211
deviation by AC"T-MD is larger than that by pure AC"T. This difference is driven by the fact that the histogram by AC"T-MD has a long tail. The long tail appears in the histogram of dk by the AC"T-MD simulation, because hydrogen atoms horizontally move for a long distance in the interlayer region of the graphite. Fig. 5a and b shows 30 samples chosen randomly in 1000 trajectories simulated by the AC"T and the AC"T-MD hybrid code, respectively. Fig. 5b shows the trajectories in the same initial condition calculated by the AC"T-MD hybrid simulation. The gray and black lines in Fig. 5b are obtained by the AC"T and the MD simulation, respectively. The portion of the trajectories of AC"T (gray lines) in Fig. 5b are exactly the same as the trajectories at kinetic energies above 200 eV in Fig. 5a. Comparing these two figures, it is clearly found that the trajectories obtained by the MD simulation in AC"T-MD tend to move horizontally for a long distance in the interlayer region of the graphite. In our previous research [5], it was shown that there is an adsorption site 1.1Å above each carbon atom of a graphene sheet. As shown in Fig. 5c, the hydrogen atom is finally trapped at the adsorption site. This phenomenon suggests that our hybrid simulation treats the chemical processes more realistically than pure AC"T simulations. Thus, it is concluded that our hybrid simulation is useful to investigate the dynamical processes where chemical reactions take place in lm scale (ex., chemical sputtering and hydrogen blister generation). 4. Summary The AC"T-MD hybrid simulation code has been developed to treat large materials with wide range of kinetic energy in plasma surface interaction. The motion of hydrogen and all recoil atoms are solved by BCA method using the AC"T code while its kinetic energy is higher than Eth. When the kinetic energy becomes lower than Eth, the MD simulation starts to calculate the motions of the hydrogen and all recoil atoms with their surrounding carbon atoms in a local box.
S211
The AC"T-MD hybrid simulation was performed and compared with pure AC"T simulations of 1000 hydrogen atom injections into a graphite. Their results suggested that the AC"T-MD hybrid simulation tends to derive a smaller value of depth than the AC"T simulation does. This fact is caused by the interaction between the surrounding atoms, which is ignored in the AC"T simulation. Acknowledgements This work is supported by the National Institutes of Natural Sciences’ undertaking Forming Bases for Interdisciplinary and International Research through Cooperation Across Fields of Study and Collaborative Research Programs (NIFS09KEIN0091 and NIFS10KEIN0161) and a Grants-in-Aid for Scientific Research (No. 19055005) from the Ministry of Education, Culture, Sports, Science and Technology, Japan. This work is performed with the support and under the auspices of the NIFS Collaboration Research Programs (NIFS09KTAN008 and NIFS10KTAS005). References [1] D.W. Brenner, O.A. Shenderova, J.A. Harrison, S.J. Stuart, B. Ni, S.B. Sinnott, J. Phys.: Condens. Matter 14 (2002) 783. [2] A. Ito, H. Nakamura, Commun. Comput. Phys. 4 (2008) 592. [3] A. Ito, H. Nakamura, J. Plasma Phys. 72 (2006) 805. [4] A. Ito, H. Nakamura, A. Takayama, J. Phys. Soc. Jpn. 77 (2008) 114602. [5] S. Saito, A. Ito and H. Nakamura, Plasma Fusion Res., special issue, 5 (2010) S2076. [6] A. Takayama, S. Saito, A.M. Ito, T. Kenmotsu and H. Nakamura, Jpn. J. Appl. Phys., in press. [arXiv:1005.3344]. [7] Y. Yamamura, Y. Mizuno, IPPJ-AM-40, Inst. Plasma Phys, Nagoya Univ., 1985. [8] A. Ito, T. Kenmotsu, Y. Kikuhara, K. Ohya, K. Inai, Y. Wang, S. Irle, K. Morokuma, H. Nakamura, The 4th IAEA-TM on the Theory of Plasma Instabilities, Kyoto University, Japan, 2009. May 18-19. [9] A. Takayama, K. Shimizu, Y. Tomita, T. Takizuka, J. Plasma Fusion Res. SERIES 9 (2010) 604. [10] A. Ito, H. Nakamura, Jpn. J. Appl. Phys. 47 (2008) 4715. [11] A. Ito, H. Nakamura, A. Takayama, J. Phys. Soc. Jpn. 77 (2008) 114602. [12] M. Suzuki, J. Math. Phys. 26 (1985) 601.