CC BY
17
УДК 538.931
Necessary Duration of Molecular Dynamics Simulation for Calculation of Self-Diffusion Coefficient during Migration of Different Point Defects in Nickel*
G.M. Poletaev1, V.V. Kovalenko2, N.M. Gurova1, M.A. Iliina3
'Polzunov Altai State Technical University (Barnaul, Russia) 2Siberian State Industrial University (Novokuznetsk, Russia) 3Financial University under the Government of the Russian Federation, Barnaul brunch (Barnaul, Russia)
Продолжительность молекулярно-динамического моделирования, необходимая для вычисления коэффициента самодиффузии при миграции различных точечных дефектов в никеле
Г.М. Полетаев1, В.В. Коваленко2, Н.М. Гурова1, М.А. Ильина3
'Алтайский государственный технический университет им. И.И. Ползунова (Барнаул, Россия)
2Сибирский государственный индустриальный университет (Новокузнецк, Россия)
3Финансовый университет при Правительстве РФ, Барнаульский филиал (Барнаул, Россия)
The evaluation of the necessary duration of molecular dynamics experiment for the calculation of the self-diffusion coefficient during migration of different point defects in Ni (vacancy, bivacancy, self-interstitial atom, hydrogen atom) is conducted in this paper. The mentioned defects have different mobility that results in different intensities of atoms displacements caused by migration of the defect. The accuracy of diffusion coefficient calculation is related to the accuracy of estimation of root-mean-square changes of atoms coordinates. Consequently, the accuracy increases with the increase of molecular-dynamic experiment duration t, the temperature T, and the mobility of the defect initiating the diffusion. To describe the interatomic interactions, the multi-particle Cleri-Rosato potential is used in the study. It is shown that the simulation duration of 100 ps is enough to calculate the diffusion coefficient when the temperature is higher than 0.6 of melting point. When calculating the diffusion coefficient of impurity in a metal crystal (for example, the hydrogen impurity), it is possible to decrease the root mean square error of displacement evaluation of impurity atoms by increasing the number of impurity atoms.
В работе проведена оценка продолжительности мо-лекулярно-динамического эксперимента, необходимой для расчета коэффициента самодиффузии при миграции различных точечных дефектов: вакансии, бивакан-сии, собственного межузельного атома, атома водорода. Перечисленные дефекты обладают различной подвижностью, в результате чего смещения атомов, возникающие вследствие миграции дефекта, имеют разную интенсивность. Погрешность определения коэффициента диффузии связана с точностью определения среднеквадра-тических изменений координат атомов, которая в свою очередь повышается с ростом продолжительности мо-лекулярно-динамического эксперимента ^ температуры Т и подвижности дефекта, инициирующего диффузию. Для описания межатомных взаимодействий в работе использовался многочастичный потенциал Клери — Розато. Показано, что для расчета коэффициента диффузии при миграции вакансии, бивакансии и межузель-ного атома при температурах выше 0,6 от температуры плавления достаточно моделирования в течение 100 пс. При расчете коэффициента диффузии примеси в кристалле металла, например водорода, одного примесного атома, как показано в настоящей работе, недостаточно. В данном случае снизить погрешность определения сред-неквадратических смещений примесных атомов можно путем введения большого числа атомов примеси.
*The reported study was partially supported by Ministry of Education and Science of the Russian Federation within the framework of the state task (project No. 3.4820.2017/8.9) and RFBR (project No. 16-48-190182-r_a).
Key words: molecular dynamics, diffusion, diffusion coefficient, point defect, vacancy, bivacancy, interstitial atom.
Ключевые слова: молекулярная динамика, диффузия, коэффициент диффузии, точечный дефект, вакансия, бивакансия, межузельный атом.
DOI 10.14258/izvasu(2018)1-06
1. Introduction
The molecular dynamics method is intensively used nowadays to model the atomic structure and processes occurring at the atomic and nano levels in various materials. This method allows studying a wide range of properties and obtaining many characteristics of materials, including kinetic and thermodynamic characteristics. The molecular dynamics method has been successfully used for thirty years to investigate diffusion involving various structural defects and to calculate diffusion characteristics, such as diffusion coefficient, diffusion activation energy, defect migration energy,
pre-exponential factor in the corresponding Arrhenius equation, using a computer model [1-4].
The diffusion coefficient in the molecular dynamics method is usually found using a well-known formula relating the diffusion coefficient and the diffusion path:
D = -
< Kff >
2t
(1)
where <x^ > — root mean square diffusion path of atoms, t — time.
Firstly, the average diffusion coefficients along the x, y and z axes are found:
1 N
N§(x» -
D =- '=1
2t
1N
N§(y* -У' D =N '=
here x , y and z0j — coordinates of initial position of the i-th atom; x, y. and z . — coordinates of the i-th atom at the moment of time t; N — the number of atoms in the calculation block. Then, the diffusion coefficient is calculated as the arithmetic mean of the coefficients D ,
D and D :
У г
D = —(D + D 3l x y
D
(3)
2t
D =
Z
1N
N § ^-2t
(2)
Once the molecular dynamics simulation is done, the calculation block should be cooled down before calculation of the diffusion coefficient. Otherwise, according to (2), thermal displacements of atoms relative to lattice sites will be taken into account during calculation of the mean square displacements of atoms.
The accuracy of diffusion coefficient calculation is related to the accuracy of evaluation of the mean square changes of atoms coordinates. So, the longer the molecular dynamics experiment duration t, the higher the temperature T and the smaller the migration energy of the defect that initiates the diffusion E , the more
m
the atomic displacements and the accuracy of diffusion coefficient calculation should be.
Typically, the computer experiment duration to calculate the diffusion coefficient is several hundred picoseconds (when the duration of one iteration is 1-10 fs) [3-9]. However, special studies on this topic are rare, especially, for different types of migrating defects.
The present study is devoted to estimation of molecular dynamics experiment duration that is necessary for diffusion coefficient calculation during migration of various point defects: a vacancy, a bivacancy, an intrinsic
interstitial atom, a hydrogen atom. For Ni which is used in this study, the migration energies of these defects are, respectively: 0.97 eV, 0.18 eV, 0.13 eV, 0.34 eV [4, 10].
2. Description of the model
The Ni calculation block of a cube form contains 8400 atoms. Periodic boundary conditions are imposed on the boundaries of the block. The interactions of nickel atoms with each other are described by the many-particle Clery-Rosato potential [11]. The experience of applying of this potential shows its helpfulness when various properties of metals should be taken into account [4, 10, 12]. When modeling of a hydrogen atom migration in a Ni crystal, interactions of Ni-H and H-H are described using Morse potentials found in [10] with the absorption energy, the activation energy of hydrogen diffusion, and the binding energy with a vacancy.
When a defect is introduced into the computational block before the main computer experiment, a preliminary dynamic relaxation of the structure near the defect is carried out. After relaxation, the computational block is cooled. A vacancy is introduced by removing one atom, a bivacancy - by removing two neighboring atoms. An interstitial atom, as is known, can be located in several positions in the crystal lattice [4]. However, when modeling its migration, its starting position does not affect the mechanism and intensity of migration. In this study, the intrinsic interstitial atom is introduced into the dumbbell position in the <100> direction, and the hydrogen atom — into the octahedral pore. According to [4, 10], these are the most energetically favorable positions for such defects.
The diffusion coefficient is calculated at different temperatures (sufficient for defect migration) with consideration of the molecular dynamics experiment duration. After a certain period, the calculation block is saved, then, for each case, the block is cooled, and the self-diffusion coefficient is calculated.
The temperature in the model is set through initial velocities of atoms according to the Maxwell-Boltzmann distribution. The total kinetic energy corresponds to the given temperature, and the total momentum of the calculation block is equal to zero. The time step in molecular dynamics experiments is 5 fs.
3. Results and discussion
Fig. 1 shows the dependence of the self-diffusion coefficient for the block of Ni crystal containing 8400 atoms and one vacancy on the duration of the molecular dynamics experiment at temperatures of 1000 K and 1500 K. The vacancy is the least mobile defect among the considered point defects, the energy of vacancy migration is greater than that of other defects.
An interstitial atom has the lowest migration energy in comparison with other considered defects, i.e., at the same temperature, diffusion with its participation
As it is shown in Fig. 1, the diffusion coefficient calculation error decreases when the duration of the molecular dynamics experiment increases. The error also decreases with the increase of temperature. As stated above, this is due to increase in atomic displacements, the number of which affects the mean square displacement calculation accuracy.
These graphs demonstrate that a computer experiment must have the duration of at least 100 ps. at temperatures above 1000 K for Ni (i.e., about 0.6 of the melting temperature).
A bivacancy is much more mobile than the vacancy, and it is shown, for example, by higher values of the diffusion coefficient (Fig. 2). According to the graphs in Fig. 2, simulation time of 100 ps is also sufficient to calculate the diffusion coefficient during bivacancy migration at temperatures above 1000 K. When the bivacancy migrates, there is a probability for the bivacancy to be divided into two separately migrating vacancies (that rises with temperature increasing). This peculiarity of bivacancy migration leads to an additional error of diffusion coefficient calculation.
proceeds comparatively more intensively. This is evident, for example, with comparatively higher values of the diffusion coefficient (Fig. 3).
a) b)
Fig. 1. Dependence of the self-diffusion coefficient for the block of Ni crystal containing 8400 atoms and one vacancy on the duration of the molecular dynamics experiment at temperatures of 1000 K (a) and 1500 K (b).
/« 1 c. _
^ 1 A
«2 14 O 12 o~ -
8 6 -
4 -
2 -
Q 50 S
C 25 50 75 100 125 150 175 200 225 2 t.F
20
fc 1H
CO
*7 1h
O
14
n 12
10
8
6
4
2
0
J 1
/ N
x /
X
25 50 75 100 125 150 175 200 225 250
t, ps
a) b)
Fig. 2. Dependence of the self-diffusion coefficient for the block of Ni crystal containing 8400 atoms and one bivacancy on the duration of the molecular dynamics experiment at temperatures of 1000 K (a) and 1500 K (b).
a) b)
Fig. 3. Dependence of the self-diffusion coefficient for the block of Ni crystal containing 8400 atoms and one interstitial atom on the duration of the molecular dynamics experiment at temperatures of 1000 K (a) and 1500 K (b).
A migration of an interstitial atom has its peculiarities. This defect, according to [4], consists of at least two migration mechanisms: the displacement and rotation of the <100> dumbbell and the crowdion mechanism. Apparently, this should introduce some error to the diffusion coefficient calculation, but this error is relatively small — as can be seen in Fig. 3. The diffusion coefficient for the interstitial atom migration is the most accurate in comparison with the values calculated for other considered point defects. It is sufficient for a molecular dynamics experiment to have its duration of 50-60 ps. to calculate the diffusion coefficient for Ni at temperatures
above 1000 K.
When calculating the diffusion coefficient of hydrogen atoms in a metal, only impurity atoms participate in the calculation of atomic displacements. With the introduction of one hydrogen atom, the calculation error for the displacement of impurity atoms (it is alone in this case) is extremely high, because the error for the diffusion coefficient is high (Fig. 4). In such cases, a greater number of defects (hydrogen atoms) are often introduced. As it is shown in Fig. 5, the calculation error for the diffusion coefficient of hydrogen is much lower for the case of 10 hydrogen atoms being introduced. The sufficient simulation duration to calculate the coefficient at the considered temperatures is also 100 ps.
a) b)
Fig. 4. Dependence of the diffusion coefficient of hydrogen in Ni at the introduction of one hydrogen atom in the calculation block containing 8400 atoms on the duration of the molecular dynamics experiment at temperatures of 1000 K (a) and 1500 K (b).
a) b)
Fig. 5. Dependence of the diffusion coefficient of hydrogen in Ni at the introduction of 10 hydrogen atoms in the calculation block containing 8400 atoms on the duration of the molecular dynamics experiment at temperatures of 1000 K (a) and 1500 K (b).
4. Summary
The evaluation of the necessary duration of a molecular dynamics experiment for the calculation of the diffusion coefficient during migration of different point defects in Ni (a vacancy, a bivacancy, a self-interstitial atom, a hydrogen atom) is conducted in this study. The listed defects have different mobility being the result of the displacements of atoms due to defect migrations with different intensities. The calculation error for the diffusion coefficient is related to the evaluation accuracy of the root-mean-square changes in the coordinates of atoms. This, in turn, increases with the duration of a molecular
dynamics experiment, the temperature and the mobility of the defect initiating the diffusion. It is shown that at the temperature higher than 0.6 of melting point, the simulation duration of 100 ps. is usually enough to calculate the diffusion coefficient during migration of a vacancy, a bivacancy, and an interstitial atom. Also, it is found out that one impurity atom is not enough to calculate the diffusion coefficient of an impurity, for example, a hydrogen atom. In this case, the calculation accuracy for the mean square displacements of impurity atoms can be increased by introducing a large number of impurity atoms.
References
1. Suzuki A., Mishin Y. Atomistic modeling of point defects and diffusion in copper grain boundary // Interface Science. — 2003. — №11. — P. 131-148.
2. Liu C.L., Plimpton S.J. Molecular-statics and molecular-dynamics study of diffusion along [001] tilt grain boundaries in Ag // Physical Review B. — 1995. — V. 51. — P. 4523-4529.
3. Frolov T., Mishin Y. Molecular dynamics modeling of self-diffusion along a triple junction // Physical Review B. — 2009. — V. 79. — 174110.
4. Poletaev G.M., Starostenkov M.D. Contributions of different mechanisms of self-diffusion in face-centered cubic metals under equilibrium conditions // Physics of the Solid State. — 2010. — V. 52, №6, P. 1146-1154.
5. Lipnitskii A.G., Nelasov I.V., Kolobov Yu.R. Self-Diffusion Parameters of Grain Boundaries and Triple Junctions in Nanocrystalline Materials // Defect and Diffusion Forum. — 2011. — V. 309-310. — P. 45-50.
6. Upmanyu M., Srolovitz D.J., Shvindlerman L.S., Gottstein G. Molecular dynamics simulation of triple junction migration // Acta Materialia. — 2002. — V. 50. — P. 1405-1420.
7. Mendelev M.I., Deng C., Schuh C.A., Srolovitz D.J. Comparison of molecular dynamics simulation methods
for the study of grain boundary migration // Modelling and Simulation in Materials Science and Engineering. — 2013. — V.21. — 045017.
8. Zhang H., Upmanyu M., Srolovitz D.J. Curvature driven grain boundary migration in aluminum: molecular dynamics simulations // Acta Materialia. — 2005. — V. 53. — P. 79-86.
9. Trautt Z.T., Mishin Y. Grain boundary migration and grain rotation studied by molecular dynamics // Acta Materialia. — 2012. — V. 60. — P. 2407-2424.
10. Poletaev G.M., Starostenkov M.D., Dmitriev S.V. Interatomic potentials in the systems Pd-H and Ni-H // Materials Physics and Mechanics. — 2016. — V. 27, №1. — P. 53-59.
11. Cleri F., Rosato V. Tight-binding potentials for transition metals and alloys // Physical Review B. — 1993. — V. 48., №1 — P. 22-33.
12. Poletaev G.M., Novoselova D.V., Kaygorodova V.M. The causes of formation of the triple junctions of grain boundaries containing excess free volume in fcc metals at crystallization // Solid State Phenomena. — 2016. — V. 249. — P. 3-8.
