- Email: [email protected]
Empirical interatomic potential for Fe-N binary system based on Finnis–Sinclair potential
Computational Materials Science 174 (2020) 109500
Contents lists available at ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
Editor’s Choice
Empirical interatomic potential for Fe-N binary system based on Finnis–Sinclair potential
T
⁎
Katsutoshi Hyodoa, , Shinji Munetohb, Toshihiro Tsuchiyamab,c,d, Setsuo Takakic a
Graduate School of Engineering, Kyushu University, 744, Motooka, Nishi-ku, Fukuoka 819-0395, Japan Department of Materials Science and Engineering, Kyushu University, 744, Motooka, Nishi-ku, Fukuoka 819-0395, Japan c Research Center for Steel, Kyushu University, 744, Motooka, Nishi-ku, Fukuoka 819-0395, Japan d International Institute for Carbon-Neutral Energy Research (WPI-I2CNER), Kyushu University, 744, Motooka, Nishi-ku, Fukuoka 819-0395, Japan b
A R T I C LE I N FO
A B S T R A C T
Keywords: Molecular dynamics (MD) Density functional theory (DFT) Iron nitrides Iron Nitrogen Steels
A new empirical potential for the Fe-N binary system based on the Finnis–Sinclair (FS) potential was developed in order to analyze the strengthening mechanism of N addition in ferritic and austenitic steels. The potential parameters for the Fe-N and N-N interactions were determined to reproduce the physical properties of the first principles calculation data of four iron nitrides (ZnS type-FeN, ε - Fe3 N , γ ′ - Fe4 N and α′ ′ - Fe16 N2 ) obtained in this study. Furthermore, we investigated the reproducibility of physical properties of the other nitrides (NaCl type-FeN, ε - Fe2 N and ζ - Fe2 N ) and the N behavior in body-centered cubic (bcc) and face-centered cubic (fcc) iron. The results showed good agreement with reference data.
1. Introduction N is an interstitial alloying element of steel, which increases the strength of ferrite and austenite by solid-solution strengthening. In austenitic stainless steels, N is often used to improve the yield strength and tensile strength, while maintaining ductility [1]. The strengthening effect of N is greater than that of C in stainless steels owing to the atomic interaction between N and Cr, which forms interstitial-substitutional (I-S) pairs or short range order (SRO) [2,3]. However, this behavior of light element N has hardly been observed experimentally. Therefore, a molecular dynamics (MD) simulation, as an atomic simulation based on an empirical atomic potential, is a powerful tool for the analysis of the N behavior in steels. So far, four atomic potentials for the Fe-N interaction have been developed. The first potential is the Johnson potential for calculating the migration characteristics and energy in bcc iron [4]. The second potential by Grujicic and Zhou uses the embedded atom method (EAM) [5], where the potential reproduces the interstitial solid solution behavior of N in fcc iron [6]. The third potential by Lee et al. is based on the modified EAM (MEAM) potential [7], and describes the physical properties of an interstitial N in bcc iron, fcc iron, and iron nitrides [8]. The forth potential by You and Yan is based on the EAM potential for the Fe-N and N-N interactions and describes the N behavior in bcc iron mainly and some physical properties in fcc iron and nitrides [9]. However, these potentials only consider ε - Fe2 N and γ ′ - Fe4 N though numerous other nitrides have been observed by ⁎
experimental methods (ZnS type-FeN, NaCl type-FeN, ζ - Fe2 N , ε - Fe3 N and α′ ′ - Fe16 N2 ) [10]. In addition, the mechanical properties (bulk modulus, Young’s modulus, etc.) may not be reproducible. This study aimed to develop new potentials for the Fe-N binary system to solve these problems. For the development of potentials, the most suitable potential function must be selected according to its purpose, because the function influences the calculation time. Recently, large-scale simulation has become possible owing to the development of the supercomputer. For example, in the grain growth simulation by Okita et al. [11], a MD calculation using the FS potential with more than 100,000,000 Fe atoms was performed [12]. The MEAM is commonly used in N-body potentials because it contains the angle dependent term; however, the potential requires a longer time to calculate the energy, and thus, it is unsuitable for large-scale simulation. In this study, Fe-N and N-N potentials based on the FS potential which requires a shorter time for the energy calculation were developed, where the parameters were decided to reproduce the physical properties (lattice parameters, bulk modulus and cohesive energy) of ZnS type-FeN, ε - Fe3 N , γ ′ - Fe4 N and α′ ′ - Fe16 N2 . In addition, the reproducibility of those of the other nitrides (NaCl type-FeN, ε - Fe2 N and ζ - Fe2 N ) and the N behavior in bcc and fcc iron were confirmed. 2. Interatomic potential The FS potential consists of repulsive core-core interactions and the
Corresponding author. E-mail address: [email protected] (K. Hyodo).
https://doi.org/10.1016/j.commatsci.2019.109500 Received 15 October 2019; Received in revised form 20 December 2019; Accepted 20 December 2019 0927-0256/ © 2019 Elsevier B.V. All rights reserved.
Computational Materials Science 174 (2020) 109500
K. Hyodo, et al.
Table 1 Potential parameters for the Fe-Fe, Fe-N, and N-N interactions.
Fe-Fe [12] Fe-N N-N
β
d
A
c
c0
c1
c2
1.8 0.1762207 1.4977494
3.569745 3.0149734 3.2659692
1.828905 2.2332911 3.2944090
3.4 2.7206212 3.2568513
1.2371147 1.8251739 2.7698948
−0.3592185 0.2901673 −1.0586444
−0.0385607 −0.5639885 0.2608420
band energy given by second moment tight-binding model. The potential has been used for transition metals, noble metals and binary alloys [12–14]. In this study, the original FS potential function by Finnis and Sinclair is adopted as the potential function for Fe-N and N-N interactions. The potential function is given as follows:
Utot = −A ∑ f (ρi ) + i
1 2
B=V
(7)
The bulk modulus of seven nitrides by new potentials was also obtained by the same method. 4. Results and discussion
∑ V (rij) i≠j
d 2E dV 2
(1) 4.1. Nitride
f (ρi ) = ρi =
ρi ,
(2) The potential parameters for the Fe-N binary system are listed in Table 1. The parameters for Fe-Fe interaction were proposed by Finnis and Sinclair [12], while those for Fe-N and N-N interactions are the results obtained in this study. By using the potential parameters, the structural optimization of the seven nitrides (ZnS type-FeN, NaCl typeFeN, ε - Fe2 N , ζ - Fe2 N , ε - Fe3 N , γ ′ - Fe4 N and α′ ′ - Fe16 N2 ) was performed without fixing internal coordinates. As a result, structures of seven nitrides were reasonably reproduced, and the physical properties of each nitride agreed well with the experimental and DFT calculation results of previous studies listed in Table 2 [12,15–27]. In the DFT calculation of NaCl-FeN, spin-polarized (SP) and non-spin-polarized (NSP) cases were described. In addition, the physical properties of α-Fe, γ-Fe, and ε-Fe are also listed in the table to compare the properties of each nitride with those of the pure iron with the same structure as a sublattice of the nitride. Each nitride expands more than that of pure iron with the same structure as the nitride sublattice. Jack reported that the structure of Fe2N nitride slightly changes from ε - Fe2 N to ζ - Fe2 N when the N content exceeds the stoichiometric composition (Fe:N/ atom = 2:1) [28]. Because the cohesive energies in Table 2 demonstrate that ζ - Fe2 N possesses a lower cohesive energy than ε - Fe2 N , Jack’s report was reasonably supported by the calculations in this study.
∑ ϕ (rij) j
(3)
(r − d )2 + β (r − d )3/ d, r ⩽ d ϕ (r ) = ⎧ ⎨ 0, r > d ⎩
(4)
(r − c )2 (c0 + c1 r + c1 r 2), r ⩽ c V (r ) = ⎧ , 0, r > c ⎨ ⎩
(5)
where, i and j are atomic labels, and rij is the bond length between i and j. The potential parameters A, β, d, c, c0, c1 and c2 are determined to reproduce the physical properties of four nitrides. As fitting targets, the cohesive energies of ZnS type-FeN, ε − Fe3 N , γ ′ - Fe4 N and α′ ′ − Fe16 N2 as a function of the volume per atom were calculated by the first principles calculation. The results of this study were compared with the lattice parameters and bulk modulus of the four nitrides reported in previous studies [15–22]. The potential parameters were varied by using the Monte-Carlo fitting to decrease the residual sum of squares between the target energy-volume curves and the curves obtained by using the new parameters, after setting seed potential parameters. The parameter fitting was performed until the errors of the lattice parameters, bulk modulus and cohesive energy became less than 10%.
4.2. N in iron 3. Calculation of the physical properties of nitrides The octahedral site (O site) and tetrahedral site (T site) are well known as interstitial sites in transition metals, and the solution energy at the O site is generally lower than that at the T site in bcc and fcc iron [29]. The interstitial sites are illustrated in Fig. 1. To confirm the validity of the potential parameters, the difference in the solution energy of a N atom between at the O and T sites (ΔE[(T) − (O)]) was calculated in bcc and fcc iron. For the calculation, the N atom was inserted at the O site or T site in a pure bcc iron cell containing 2904 Fe atoms (11 × 11 × 12 bcc unit cells) or a fcc cell containing 2916 atoms (9 × 9 × 9 fcc unit cells). The O and T positions were set near the center positions of the bcc and fcc cells. Then, structural relaxation was performed with fixed internal coordinates of the N and Fe atoms around the periodic boundary. The following structural relaxations were also performed under the same condition. The values of ΔE[(T) − (O)] were calculated at 0.35 eV in bcc iron and 2.32 eV in fcc iron. The solution energy at the O site was lower in bcc and fcc iron, and this result agreed with the general theory. The migration barriers of N atoms in ferrite and austenite have been reported as 0.76 eV and 1.75 eV, respectively [30,31]. And then, Lee et al. reported 0.78 eV in bcc iron and 1.36 eV in fcc iron by the MEAM potential [8]. The energy barriers were defined as the energy difference between O and T sites in bcc iron and the energy barrier from O site to O site along the 〈1 1 0〉 direction in fcc iron. In this study, the diffusion paths in bcc and fcc iron were defined as from O to T to O and from O to O, respectively (shown in Fig. 1). The energy
The cohesive energy-volume relation was calculated by using the first principles calculation code Advance/PHASE. Total energy calculations and structural optimizations were performed by using the Perdew-Burke-Ernzerhof generalized gradient approximation (GGAPBE) of the density functional theory (DFT) calculation and ultrasoft pseudopotential. 340 eV was selected as the cutoff energy, and the k point meshes were set at 8 × 8 × 8, 6 × 6 × 6, 8 × 8 × 8 and 4 × 4 × 4 for ZnS type-FeN, ε - Fe3 N , γ ′ - Fe4 N and α′ ′ - Fe16 N2 , respectively. Then, the internal coordinates of the N atoms were fixed in structural optimizations of ε - Fe3 N and α′ ′ - Fe16 N2 . For the DFT calculations of the nitrides with ferro magnetic characteristics which have been confirmed by Mössbauer spectroscopy [10], the calculations were performed considering the spin. To obtain the cohesive energy, the total energy calculation of one Fe or N atom was performed in a 20 Å × 20 Å × 20 Å vacuum. In addition, the k point mesh was 6 × 6 × 6. The cohesive energy of nitrides is given as follows:
Ecoh (Fe x Ny) = Etot (Fe x Ny ) − xEtot (Fe) − yEtot (N)
(6)
In addition, to obtain the bulk modulus of the nitrides, the change in the cohesive energy corresponding to that in the volume per atom was calculated, and then, the energy points around the energy bottom were fitted to the quadratic curves (E) as a function of volume (V). The bulk modulus (B) was calculated by the following equation: 2
Computational Materials Science 174 (2020) 109500
K. Hyodo, et al.
Table 2 Structural parameters, bulk modulus and cohesive energies of pure irons and iron nitrides using the potential parameters for the Fe-Fe, Fe-N and N-N interactions.
α-Fe
γ-Fe
ε-Fe
Method
Structural parameters (Å)
Bulk modulus B (GPa)
Cohesive energy Ecoh (eV/atom)
Other work [12] Exp. [12] Other DFT [23]
a = 2.8665
173.1
−4.28
a = 2.8665 a = 2.84
173.1 174
−4.28
This work Exp. [24] Other DFT [23]
a = 3.6938 a = 3.65 a = 3.64
156
−4.23
This work
a = 2.6119 c = 4.2652 c/a = 1.633
Other DFT [25]
α″-Fe16N2
This work
This DFT
Exp. [15]
Other DFT [16]
ζ-Fe2N
ε-Fe3N
156
This work
a = 4.3444
210
−5.42
216
−5.40
−4.23 Other DFT [27]
√3a/c = 1.585 2a/b = 1.615 a = 4.411 b = 5.509 c = 4.771 √3a/c = 1.60 2a/b = 1.60 a = 4.3406
a = 2.525
b = 5.4480
c = 4.375 c/a = 1.73
c = 4.7544 √3a/c = 1.581 2a/b = 1.593
a = 5.6120 c = 6.2329 c/a = 1.111 a = 5.682 c = 6.255 c/a = 1.10 a = 5.72 c = 6.31 c/a = 1.1 a = 5.70
165
−4.67 ε-Fe2N
187
−5.04
Other DFT [27]
ZnS typeFeN
a = 3.79 a = 3.792
196
This work
a = 4.6096
184
201
a = 4.7323 c = 4.3148 √3c/a = 1.579 a = 4.794 c = 4.418 √3c/a = 1.596 a = 4.7634 c = 4.3273 √3c/a = 1.573
180
Exp. [17] Other DFT [18]
c = 4.3458 √3c/a = 1.633 a = 4.654 c = 4.324 √3c/a = 1.61 a = 4.6982 c = 4.3789 √3c/a = 1.614 a = 4.69140
This work
Exp. [26]
182 200
Other DFT [20]
Cohesive energy Ecoh (eV/atom)
Exp. [26]
a = 3.7507 a = 3.794
Exp. [19]
Bulk modulus B (GPa)
167
This work This DFT
This DFT
Structural parameters (Å)
b = 5.3792 c = 4.7483
c = 6.16 c/a = 1.08
γ′-Fe4N
Method
−4.99 −5.05
NaCl type-FeN −5.19
This work
a = 4.2115
246
−4.87
This DFT Exp. [21] Other DFT [22]
a = 4.236 a = 4.33 a = 4.215
277
−5.09
324
This work
a = 4.0450
384
Exp. [21] Other DFT [22]
a = 4.50 a = 4.015 (SP)
288 (SP)
a = 3.980 (NSP)
369 (NSP)
−5.61
−5.07
170.9
c = 4.34551 √3c/a = 1.6043
Δd denote the bond length before structural relaxation and the difference of the bond length between before and after structural relaxation, respectively. Although the relaxation of 2NN at the T site was larger than the first principles calculation result, the other values were reasonable. The tendency of the contraction or expansion of the Fe-N bonds obtained by these calculations were in good agreement with reference data.
barriers of these paths are shown in Fig. 2. This figure was obtained by moving the N atom from O to O in steps of 1/500 and performing the structural relaxation at every step. The migration barrier values were 0.35 eV in bcc iron and 1.48 eV in fcc iron, and these values were similar to the experimental results. In addition, the changes in the distance from the N atom to the first nearest neighbor (1NN) and second nearest neighbor (2NN) were estimated with structural relaxation after inserting a N atom in the O site, T site and substitutional site (S site) in the bcc cell. When inserting the N atom in the S site, an iron atom near the center position in the bcc cell was replaced with the N atom. Table 3 lists the calculation results of the fraction of relaxation (in percent of Δd/d0) by inserting the N atom in bcc iron. Here, d0 and
4.3. Phonon density of states of γ ′ - Fe4 N Above results, the static behaviors of N were discussed. Furthermore, we researched the phonon density of states of γ ′ - Fe4 N as 3
Computational Materials Science 174 (2020) 109500
K. Hyodo, et al.
Fig. 1. Illustrations of interstitial sites and diffusion paths (a) in bcc iron and (b) in fcc iron. Closed circles, triangles and allows indicate O sites, T sites and diffusion paths, respectively.
Fig. 2. Diffusion path and migration barrier of an interstitial N atom (a) in bcc iron and (b) in fcc iron.
powder containing a small amount of α-Fe (less than 8 vol%) [34]. The phonon density of states of a stainless steel (310ss) measured by the nuclear resonant inelastic scattering method was used as a substitute of the experimental data of γ-Fe, because the pure γ-Fe is unstable under the measurable condition [35]. As same as reference data, below 40 meV of MD simulations of γ ′ - Fe4 N and γ-Fe, the Fe-Fe phonon mode was confirmed. Above 50 meV, although the energy peak of Fe-N phonon mode in this study is located at low energy side, the peak was reproduced.
Table 3 Relaxation (in percent of Δd/d0) by a N atom at the O site, T site, or S site in bcc iron. Method
This work Exp. [32] Other DFT [33]
O site
T site
S site
1NN
2NN
1NN
2NN
1NN
2NN
+17 +35 +23.4
−4.1 −4.2 −1.8
+8.5
−7.9
−11
+3.6
+12.0
−0.6
−4.4
+2.2
5. Conclusions a dynamic behavior of N by using MD simulations based on the new potentials. MD simulations were performed under the constant NVT condition by using Langevin equations. A 5000 atoms γ ′ - Fe4 N cell (10 × 10 × 10 unit cells) and a γ-Fe cell containing 4000 Fe atoms as comparison data were prepared. These cells were kept at 300 K for 1.0 ps, and the phonon density of states was calculated by the Fourier transformation of the velocity auto-correlation function (Z(t) ) among the simulations. The function is defined by below equation:
Z (t ) =
〈v (t )·v (0) 〉 〈v (0)·v (0) 〉
To confirm the reliability of the new empirical potentials for the FeN and N-N interactions based on the FS potential, various examinations were performed in this study. The following results were obtained: (1) The potentials reproduce the crystal structure, lattice parameters, bulk modulus of all nitrides observed by experimental methods and calculated by first principles calculations. (2) The difference in the cohesive energy of ε - Fe2 N and ζ - Fe2 N by the potential corresponds to Jack’s report, where the structure of Fe2N changes from ε to ζ when the stoichiometric composition of Fe and N is 2:1. (3) The N diffusion behavior in bcc and fcc iron and the changes in the distance from the N atom to 1NN and 2NN in BCC iron were in good agreement with the experimental and first principles calculations results.
(8)
Here, v (t ) is the velocity vector of an atom at time t, and the bracket means calculating the average of all atoms every MD step. Fig. 3 shows the phonon densities of states of γ ′ - Fe4 N and γ-Fe by the MD simulations and the reference data. The experimental data of γ ′ - Fe4 N was obtained by the inelastic neutron scattering measurement of γ ′ - Fe4 N 4
Computational Materials Science 174 (2020) 109500
K. Hyodo, et al.
interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgement This work was supported by JSPS KAKENHI Grant Number JP19J20306 and JP15H05768. References [1] T. Tsuchiyama, H. Ito, K. Kataoka, S. Takaki, Metall. Mater. Trans. 34A (2003) 2591–2599. [2] M.L.G. Byrnes, M. Grujicic, W.S. Owen, Acta Metall. 35 (1987) 1853–1862. [3] M. Grujicic, W.S. Owen, Acta Metall. Mater. 43 (1995) 4201–4211. [4] R.A. Johnson, G.J. Dienes, A.C. Damask, Acta Metall. 12 (1964) 1215–1224. [5] M.S. Daw, M.I. Baskes, Phys. Rev. B 29 (1984) 6443–6453. [6] M. Grujicic, X.W. Zhou, Calphad 17 (1993) 383–413. [7] M.I. Baskes, Phys. Rev. B 46 (1992) 2727–2742. [8] B.-J. Lee, T.-H. Lee, S.-J. Kim, Acta Mater. 54 (2006) 4597–4607. [9] Y. You, M.F. Yan, Phil. Mag. Lett. 92 (2012) 656–667. [10] S. Nasu, T. Hinomura, Materia Japan 36 (1997) 35–39. [11] S. Okita, E. Miyoshi, S. Sakane, T. Takaki, M. Ohno, Y. Shibuta, Acta Mater. 153 (2018) 108–116. [12] M.W. Finnis, J.E. Sinclair, Phil. Mag. A 50 (1984) 45–55. [13] G.J. Ackland, G. Tichy, V. Vitek, M.W. Finnis, Phil. Mag. A 56 (1987) 735–756. [14] G.J. Ackland, D.J. Bacon, A.F. Calder, T. Harry, Phil. Mag. A 75 (1997) 713–732. [15] H. Tanaka, S. Nagakura, Y. Nakamura, Y. Hirotsu, Acta Mater. 45 (1997) 1401–1410. [16] Y.J. Shi, Y.L. Du, G. Chen, Physica B 407 (2012) 3423–3426. [17] J.L. Costa-Krämer, D.M. Borsa, J.M. García-Martín, M.S. Martín-González, D.O. Boerma, F. Briones, Phys. Rev. B 69 (2004) 144402. [18] T. Takahashi, J. Burghaus, D. Music, R. Dronskowski, J.M. Schneider, Acta Mater. 60 (2012) 2054–2060. [19] D. Rechenbach, H. Jacobs, J. Alloys Compd. 235 (1996) 15–22. [20] W.H. Zhang, Z.Q. Lv, Z.P. Shi, S.H. Sun, Z.H. Wang, W.T. Fu, J. Magn. Magn. Mater. 324 (2012) 2271–2276. [21] H. Nakagawa, S. Nasu, H. Fujii, M. Takahashi, F. Kanamaru, Hyperfine Interactions 69 (1991) 455–458. [22] P. Lukashev, W.R.L. Lambrecht, Phys. Rev. B 70 (2004) 245205. [23] D.E. Jiang, E.A. Carter, Phys. Rev. B 67 (2003) 214103. [24] M. Acet, H. Zähres, E.F. Wassermann, Phys. Rev. B 49 (1994) 6012–6017. [25] C.M. Fang, M.A. van Huis, Heliyon 3 (2017) e00408. [26] S. Nagakura, K. Tanehashi, J. Phys. Soc. Jpn. 25 (1968) 840–846. [27] C.M. Fang, M.A. van Huis, H.W. Zandbergen, Scripta Mater. 64 (2011) 296–299. [28] K.H. Jack, Acta Cryst. 5 (1952) 404–411. [29] E. Fujita, Materia Japan 6 (1967) 647–656. [30] A.D. Le Claire, Diffusion of C, N, and O in metals, in: H. Mehrer (Eds.), Diffusion in festen Metallen und Legierungen, Springer, Berlin, 1990, pp. 481–503. [31] Diffusion in metals, in: E.A. Brandes, G.B. Brook (Eds.), Smithells Metals Reference Book, seventh ed., Butterworth-Heinemann, Oxford, 1992, pp. 13. [32] K.H. Jack, Proc. Roy. Soc. A 208 (1951) 200–215. [33] C. Domain, C.S. Becquart, J. Foct, Phys. Rev. B 69 (2004) 144112. [34] K.-J. Kim, K. Sumiyama, K. Shibata, K. Suzuki, Mater. Sci. Eng. A 181–182 (1994) 1272–1276. [35] Y. Tsunoda, Y. Kurimoto, M. Seto, S. Kitao, Y. Yoda, Phys. Rev. B 66 (2002) 214304.
Fig. 3. Phonon densities of states of γ′-Fe4N and γ-Fe. Open and filled circles indicate experimental results of γ′-Fe4N and γ-Fe, respectively.
(4) The phonon density of states of γ ′ - Fe4 N was reproduced by MD simulations using the new potentials. CRediT authorship contribution statement Katsutoshi Hyodo: Software, Investigation, Writing - original draft, Visualization, Funding acquisition. Shinji Munetoh: Methodology, Software, Resources, Writing - review & editing. Toshihiro Tsuchiyama: Conceptualization, Resources, Writing - review & editing, Project administration. Setsuo Takaki: Supervision, Funding acquisition. Declaration of Competing Interest The authors declare that they have no known competing financial
5
Contents lists available at ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
Editor’s Choice
Empirical interatomic potential for Fe-N binary system based on Finnis–Sinclair potential
T
⁎
Katsutoshi Hyodoa, , Shinji Munetohb, Toshihiro Tsuchiyamab,c,d, Setsuo Takakic a
Graduate School of Engineering, Kyushu University, 744, Motooka, Nishi-ku, Fukuoka 819-0395, Japan Department of Materials Science and Engineering, Kyushu University, 744, Motooka, Nishi-ku, Fukuoka 819-0395, Japan c Research Center for Steel, Kyushu University, 744, Motooka, Nishi-ku, Fukuoka 819-0395, Japan d International Institute for Carbon-Neutral Energy Research (WPI-I2CNER), Kyushu University, 744, Motooka, Nishi-ku, Fukuoka 819-0395, Japan b
A R T I C LE I N FO
A B S T R A C T
Keywords: Molecular dynamics (MD) Density functional theory (DFT) Iron nitrides Iron Nitrogen Steels
A new empirical potential for the Fe-N binary system based on the Finnis–Sinclair (FS) potential was developed in order to analyze the strengthening mechanism of N addition in ferritic and austenitic steels. The potential parameters for the Fe-N and N-N interactions were determined to reproduce the physical properties of the first principles calculation data of four iron nitrides (ZnS type-FeN, ε - Fe3 N , γ ′ - Fe4 N and α′ ′ - Fe16 N2 ) obtained in this study. Furthermore, we investigated the reproducibility of physical properties of the other nitrides (NaCl type-FeN, ε - Fe2 N and ζ - Fe2 N ) and the N behavior in body-centered cubic (bcc) and face-centered cubic (fcc) iron. The results showed good agreement with reference data.
1. Introduction N is an interstitial alloying element of steel, which increases the strength of ferrite and austenite by solid-solution strengthening. In austenitic stainless steels, N is often used to improve the yield strength and tensile strength, while maintaining ductility [1]. The strengthening effect of N is greater than that of C in stainless steels owing to the atomic interaction between N and Cr, which forms interstitial-substitutional (I-S) pairs or short range order (SRO) [2,3]. However, this behavior of light element N has hardly been observed experimentally. Therefore, a molecular dynamics (MD) simulation, as an atomic simulation based on an empirical atomic potential, is a powerful tool for the analysis of the N behavior in steels. So far, four atomic potentials for the Fe-N interaction have been developed. The first potential is the Johnson potential for calculating the migration characteristics and energy in bcc iron [4]. The second potential by Grujicic and Zhou uses the embedded atom method (EAM) [5], where the potential reproduces the interstitial solid solution behavior of N in fcc iron [6]. The third potential by Lee et al. is based on the modified EAM (MEAM) potential [7], and describes the physical properties of an interstitial N in bcc iron, fcc iron, and iron nitrides [8]. The forth potential by You and Yan is based on the EAM potential for the Fe-N and N-N interactions and describes the N behavior in bcc iron mainly and some physical properties in fcc iron and nitrides [9]. However, these potentials only consider ε - Fe2 N and γ ′ - Fe4 N though numerous other nitrides have been observed by ⁎
experimental methods (ZnS type-FeN, NaCl type-FeN, ζ - Fe2 N , ε - Fe3 N and α′ ′ - Fe16 N2 ) [10]. In addition, the mechanical properties (bulk modulus, Young’s modulus, etc.) may not be reproducible. This study aimed to develop new potentials for the Fe-N binary system to solve these problems. For the development of potentials, the most suitable potential function must be selected according to its purpose, because the function influences the calculation time. Recently, large-scale simulation has become possible owing to the development of the supercomputer. For example, in the grain growth simulation by Okita et al. [11], a MD calculation using the FS potential with more than 100,000,000 Fe atoms was performed [12]. The MEAM is commonly used in N-body potentials because it contains the angle dependent term; however, the potential requires a longer time to calculate the energy, and thus, it is unsuitable for large-scale simulation. In this study, Fe-N and N-N potentials based on the FS potential which requires a shorter time for the energy calculation were developed, where the parameters were decided to reproduce the physical properties (lattice parameters, bulk modulus and cohesive energy) of ZnS type-FeN, ε - Fe3 N , γ ′ - Fe4 N and α′ ′ - Fe16 N2 . In addition, the reproducibility of those of the other nitrides (NaCl type-FeN, ε - Fe2 N and ζ - Fe2 N ) and the N behavior in bcc and fcc iron were confirmed. 2. Interatomic potential The FS potential consists of repulsive core-core interactions and the
Corresponding author. E-mail address: [email protected] (K. Hyodo).
https://doi.org/10.1016/j.commatsci.2019.109500 Received 15 October 2019; Received in revised form 20 December 2019; Accepted 20 December 2019 0927-0256/ © 2019 Elsevier B.V. All rights reserved.
Computational Materials Science 174 (2020) 109500
K. Hyodo, et al.
Table 1 Potential parameters for the Fe-Fe, Fe-N, and N-N interactions.
Fe-Fe [12] Fe-N N-N
β
d
A
c
c0
c1
c2
1.8 0.1762207 1.4977494
3.569745 3.0149734 3.2659692
1.828905 2.2332911 3.2944090
3.4 2.7206212 3.2568513
1.2371147 1.8251739 2.7698948
−0.3592185 0.2901673 −1.0586444
−0.0385607 −0.5639885 0.2608420
band energy given by second moment tight-binding model. The potential has been used for transition metals, noble metals and binary alloys [12–14]. In this study, the original FS potential function by Finnis and Sinclair is adopted as the potential function for Fe-N and N-N interactions. The potential function is given as follows:
Utot = −A ∑ f (ρi ) + i
1 2
B=V
(7)
The bulk modulus of seven nitrides by new potentials was also obtained by the same method. 4. Results and discussion
∑ V (rij) i≠j
d 2E dV 2
(1) 4.1. Nitride
f (ρi ) = ρi =
ρi ,
(2) The potential parameters for the Fe-N binary system are listed in Table 1. The parameters for Fe-Fe interaction were proposed by Finnis and Sinclair [12], while those for Fe-N and N-N interactions are the results obtained in this study. By using the potential parameters, the structural optimization of the seven nitrides (ZnS type-FeN, NaCl typeFeN, ε - Fe2 N , ζ - Fe2 N , ε - Fe3 N , γ ′ - Fe4 N and α′ ′ - Fe16 N2 ) was performed without fixing internal coordinates. As a result, structures of seven nitrides were reasonably reproduced, and the physical properties of each nitride agreed well with the experimental and DFT calculation results of previous studies listed in Table 2 [12,15–27]. In the DFT calculation of NaCl-FeN, spin-polarized (SP) and non-spin-polarized (NSP) cases were described. In addition, the physical properties of α-Fe, γ-Fe, and ε-Fe are also listed in the table to compare the properties of each nitride with those of the pure iron with the same structure as a sublattice of the nitride. Each nitride expands more than that of pure iron with the same structure as the nitride sublattice. Jack reported that the structure of Fe2N nitride slightly changes from ε - Fe2 N to ζ - Fe2 N when the N content exceeds the stoichiometric composition (Fe:N/ atom = 2:1) [28]. Because the cohesive energies in Table 2 demonstrate that ζ - Fe2 N possesses a lower cohesive energy than ε - Fe2 N , Jack’s report was reasonably supported by the calculations in this study.
∑ ϕ (rij) j
(3)
(r − d )2 + β (r − d )3/ d, r ⩽ d ϕ (r ) = ⎧ ⎨ 0, r > d ⎩
(4)
(r − c )2 (c0 + c1 r + c1 r 2), r ⩽ c V (r ) = ⎧ , 0, r > c ⎨ ⎩
(5)
where, i and j are atomic labels, and rij is the bond length between i and j. The potential parameters A, β, d, c, c0, c1 and c2 are determined to reproduce the physical properties of four nitrides. As fitting targets, the cohesive energies of ZnS type-FeN, ε − Fe3 N , γ ′ - Fe4 N and α′ ′ − Fe16 N2 as a function of the volume per atom were calculated by the first principles calculation. The results of this study were compared with the lattice parameters and bulk modulus of the four nitrides reported in previous studies [15–22]. The potential parameters were varied by using the Monte-Carlo fitting to decrease the residual sum of squares between the target energy-volume curves and the curves obtained by using the new parameters, after setting seed potential parameters. The parameter fitting was performed until the errors of the lattice parameters, bulk modulus and cohesive energy became less than 10%.
4.2. N in iron 3. Calculation of the physical properties of nitrides The octahedral site (O site) and tetrahedral site (T site) are well known as interstitial sites in transition metals, and the solution energy at the O site is generally lower than that at the T site in bcc and fcc iron [29]. The interstitial sites are illustrated in Fig. 1. To confirm the validity of the potential parameters, the difference in the solution energy of a N atom between at the O and T sites (ΔE[(T) − (O)]) was calculated in bcc and fcc iron. For the calculation, the N atom was inserted at the O site or T site in a pure bcc iron cell containing 2904 Fe atoms (11 × 11 × 12 bcc unit cells) or a fcc cell containing 2916 atoms (9 × 9 × 9 fcc unit cells). The O and T positions were set near the center positions of the bcc and fcc cells. Then, structural relaxation was performed with fixed internal coordinates of the N and Fe atoms around the periodic boundary. The following structural relaxations were also performed under the same condition. The values of ΔE[(T) − (O)] were calculated at 0.35 eV in bcc iron and 2.32 eV in fcc iron. The solution energy at the O site was lower in bcc and fcc iron, and this result agreed with the general theory. The migration barriers of N atoms in ferrite and austenite have been reported as 0.76 eV and 1.75 eV, respectively [30,31]. And then, Lee et al. reported 0.78 eV in bcc iron and 1.36 eV in fcc iron by the MEAM potential [8]. The energy barriers were defined as the energy difference between O and T sites in bcc iron and the energy barrier from O site to O site along the 〈1 1 0〉 direction in fcc iron. In this study, the diffusion paths in bcc and fcc iron were defined as from O to T to O and from O to O, respectively (shown in Fig. 1). The energy
The cohesive energy-volume relation was calculated by using the first principles calculation code Advance/PHASE. Total energy calculations and structural optimizations were performed by using the Perdew-Burke-Ernzerhof generalized gradient approximation (GGAPBE) of the density functional theory (DFT) calculation and ultrasoft pseudopotential. 340 eV was selected as the cutoff energy, and the k point meshes were set at 8 × 8 × 8, 6 × 6 × 6, 8 × 8 × 8 and 4 × 4 × 4 for ZnS type-FeN, ε - Fe3 N , γ ′ - Fe4 N and α′ ′ - Fe16 N2 , respectively. Then, the internal coordinates of the N atoms were fixed in structural optimizations of ε - Fe3 N and α′ ′ - Fe16 N2 . For the DFT calculations of the nitrides with ferro magnetic characteristics which have been confirmed by Mössbauer spectroscopy [10], the calculations were performed considering the spin. To obtain the cohesive energy, the total energy calculation of one Fe or N atom was performed in a 20 Å × 20 Å × 20 Å vacuum. In addition, the k point mesh was 6 × 6 × 6. The cohesive energy of nitrides is given as follows:
Ecoh (Fe x Ny) = Etot (Fe x Ny ) − xEtot (Fe) − yEtot (N)
(6)
In addition, to obtain the bulk modulus of the nitrides, the change in the cohesive energy corresponding to that in the volume per atom was calculated, and then, the energy points around the energy bottom were fitted to the quadratic curves (E) as a function of volume (V). The bulk modulus (B) was calculated by the following equation: 2
Computational Materials Science 174 (2020) 109500
K. Hyodo, et al.
Table 2 Structural parameters, bulk modulus and cohesive energies of pure irons and iron nitrides using the potential parameters for the Fe-Fe, Fe-N and N-N interactions.
α-Fe
γ-Fe
ε-Fe
Method
Structural parameters (Å)
Bulk modulus B (GPa)
Cohesive energy Ecoh (eV/atom)
Other work [12] Exp. [12] Other DFT [23]
a = 2.8665
173.1
−4.28
a = 2.8665 a = 2.84
173.1 174
−4.28
This work Exp. [24] Other DFT [23]
a = 3.6938 a = 3.65 a = 3.64
156
−4.23
This work
a = 2.6119 c = 4.2652 c/a = 1.633
Other DFT [25]
α″-Fe16N2
This work
This DFT
Exp. [15]
Other DFT [16]
ζ-Fe2N
ε-Fe3N
156
This work
a = 4.3444
210
−5.42
216
−5.40
−4.23 Other DFT [27]
√3a/c = 1.585 2a/b = 1.615 a = 4.411 b = 5.509 c = 4.771 √3a/c = 1.60 2a/b = 1.60 a = 4.3406
a = 2.525
b = 5.4480
c = 4.375 c/a = 1.73
c = 4.7544 √3a/c = 1.581 2a/b = 1.593
a = 5.6120 c = 6.2329 c/a = 1.111 a = 5.682 c = 6.255 c/a = 1.10 a = 5.72 c = 6.31 c/a = 1.1 a = 5.70
165
−4.67 ε-Fe2N
187
−5.04
Other DFT [27]
ZnS typeFeN
a = 3.79 a = 3.792
196
This work
a = 4.6096
184
201
a = 4.7323 c = 4.3148 √3c/a = 1.579 a = 4.794 c = 4.418 √3c/a = 1.596 a = 4.7634 c = 4.3273 √3c/a = 1.573
180
Exp. [17] Other DFT [18]
c = 4.3458 √3c/a = 1.633 a = 4.654 c = 4.324 √3c/a = 1.61 a = 4.6982 c = 4.3789 √3c/a = 1.614 a = 4.69140
This work
Exp. [26]
182 200
Other DFT [20]
Cohesive energy Ecoh (eV/atom)
Exp. [26]
a = 3.7507 a = 3.794
Exp. [19]
Bulk modulus B (GPa)
167
This work This DFT
This DFT
Structural parameters (Å)
b = 5.3792 c = 4.7483
c = 6.16 c/a = 1.08
γ′-Fe4N
Method
−4.99 −5.05
NaCl type-FeN −5.19
This work
a = 4.2115
246
−4.87
This DFT Exp. [21] Other DFT [22]
a = 4.236 a = 4.33 a = 4.215
277
−5.09
324
This work
a = 4.0450
384
Exp. [21] Other DFT [22]
a = 4.50 a = 4.015 (SP)
288 (SP)
a = 3.980 (NSP)
369 (NSP)
−5.61
−5.07
170.9
c = 4.34551 √3c/a = 1.6043
Δd denote the bond length before structural relaxation and the difference of the bond length between before and after structural relaxation, respectively. Although the relaxation of 2NN at the T site was larger than the first principles calculation result, the other values were reasonable. The tendency of the contraction or expansion of the Fe-N bonds obtained by these calculations were in good agreement with reference data.
barriers of these paths are shown in Fig. 2. This figure was obtained by moving the N atom from O to O in steps of 1/500 and performing the structural relaxation at every step. The migration barrier values were 0.35 eV in bcc iron and 1.48 eV in fcc iron, and these values were similar to the experimental results. In addition, the changes in the distance from the N atom to the first nearest neighbor (1NN) and second nearest neighbor (2NN) were estimated with structural relaxation after inserting a N atom in the O site, T site and substitutional site (S site) in the bcc cell. When inserting the N atom in the S site, an iron atom near the center position in the bcc cell was replaced with the N atom. Table 3 lists the calculation results of the fraction of relaxation (in percent of Δd/d0) by inserting the N atom in bcc iron. Here, d0 and
4.3. Phonon density of states of γ ′ - Fe4 N Above results, the static behaviors of N were discussed. Furthermore, we researched the phonon density of states of γ ′ - Fe4 N as 3
Computational Materials Science 174 (2020) 109500
K. Hyodo, et al.
Fig. 1. Illustrations of interstitial sites and diffusion paths (a) in bcc iron and (b) in fcc iron. Closed circles, triangles and allows indicate O sites, T sites and diffusion paths, respectively.
Fig. 2. Diffusion path and migration barrier of an interstitial N atom (a) in bcc iron and (b) in fcc iron.
powder containing a small amount of α-Fe (less than 8 vol%) [34]. The phonon density of states of a stainless steel (310ss) measured by the nuclear resonant inelastic scattering method was used as a substitute of the experimental data of γ-Fe, because the pure γ-Fe is unstable under the measurable condition [35]. As same as reference data, below 40 meV of MD simulations of γ ′ - Fe4 N and γ-Fe, the Fe-Fe phonon mode was confirmed. Above 50 meV, although the energy peak of Fe-N phonon mode in this study is located at low energy side, the peak was reproduced.
Table 3 Relaxation (in percent of Δd/d0) by a N atom at the O site, T site, or S site in bcc iron. Method
This work Exp. [32] Other DFT [33]
O site
T site
S site
1NN
2NN
1NN
2NN
1NN
2NN
+17 +35 +23.4
−4.1 −4.2 −1.8
+8.5
−7.9
−11
+3.6
+12.0
−0.6
−4.4
+2.2
5. Conclusions a dynamic behavior of N by using MD simulations based on the new potentials. MD simulations were performed under the constant NVT condition by using Langevin equations. A 5000 atoms γ ′ - Fe4 N cell (10 × 10 × 10 unit cells) and a γ-Fe cell containing 4000 Fe atoms as comparison data were prepared. These cells were kept at 300 K for 1.0 ps, and the phonon density of states was calculated by the Fourier transformation of the velocity auto-correlation function (Z(t) ) among the simulations. The function is defined by below equation:
Z (t ) =
〈v (t )·v (0) 〉 〈v (0)·v (0) 〉
To confirm the reliability of the new empirical potentials for the FeN and N-N interactions based on the FS potential, various examinations were performed in this study. The following results were obtained: (1) The potentials reproduce the crystal structure, lattice parameters, bulk modulus of all nitrides observed by experimental methods and calculated by first principles calculations. (2) The difference in the cohesive energy of ε - Fe2 N and ζ - Fe2 N by the potential corresponds to Jack’s report, where the structure of Fe2N changes from ε to ζ when the stoichiometric composition of Fe and N is 2:1. (3) The N diffusion behavior in bcc and fcc iron and the changes in the distance from the N atom to 1NN and 2NN in BCC iron were in good agreement with the experimental and first principles calculations results.
(8)
Here, v (t ) is the velocity vector of an atom at time t, and the bracket means calculating the average of all atoms every MD step. Fig. 3 shows the phonon densities of states of γ ′ - Fe4 N and γ-Fe by the MD simulations and the reference data. The experimental data of γ ′ - Fe4 N was obtained by the inelastic neutron scattering measurement of γ ′ - Fe4 N 4
Computational Materials Science 174 (2020) 109500
K. Hyodo, et al.
interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgement This work was supported by JSPS KAKENHI Grant Number JP19J20306 and JP15H05768. References [1] T. Tsuchiyama, H. Ito, K. Kataoka, S. Takaki, Metall. Mater. Trans. 34A (2003) 2591–2599. [2] M.L.G. Byrnes, M. Grujicic, W.S. Owen, Acta Metall. 35 (1987) 1853–1862. [3] M. Grujicic, W.S. Owen, Acta Metall. Mater. 43 (1995) 4201–4211. [4] R.A. Johnson, G.J. Dienes, A.C. Damask, Acta Metall. 12 (1964) 1215–1224. [5] M.S. Daw, M.I. Baskes, Phys. Rev. B 29 (1984) 6443–6453. [6] M. Grujicic, X.W. Zhou, Calphad 17 (1993) 383–413. [7] M.I. Baskes, Phys. Rev. B 46 (1992) 2727–2742. [8] B.-J. Lee, T.-H. Lee, S.-J. Kim, Acta Mater. 54 (2006) 4597–4607. [9] Y. You, M.F. Yan, Phil. Mag. Lett. 92 (2012) 656–667. [10] S. Nasu, T. Hinomura, Materia Japan 36 (1997) 35–39. [11] S. Okita, E. Miyoshi, S. Sakane, T. Takaki, M. Ohno, Y. Shibuta, Acta Mater. 153 (2018) 108–116. [12] M.W. Finnis, J.E. Sinclair, Phil. Mag. A 50 (1984) 45–55. [13] G.J. Ackland, G. Tichy, V. Vitek, M.W. Finnis, Phil. Mag. A 56 (1987) 735–756. [14] G.J. Ackland, D.J. Bacon, A.F. Calder, T. Harry, Phil. Mag. A 75 (1997) 713–732. [15] H. Tanaka, S. Nagakura, Y. Nakamura, Y. Hirotsu, Acta Mater. 45 (1997) 1401–1410. [16] Y.J. Shi, Y.L. Du, G. Chen, Physica B 407 (2012) 3423–3426. [17] J.L. Costa-Krämer, D.M. Borsa, J.M. García-Martín, M.S. Martín-González, D.O. Boerma, F. Briones, Phys. Rev. B 69 (2004) 144402. [18] T. Takahashi, J. Burghaus, D. Music, R. Dronskowski, J.M. Schneider, Acta Mater. 60 (2012) 2054–2060. [19] D. Rechenbach, H. Jacobs, J. Alloys Compd. 235 (1996) 15–22. [20] W.H. Zhang, Z.Q. Lv, Z.P. Shi, S.H. Sun, Z.H. Wang, W.T. Fu, J. Magn. Magn. Mater. 324 (2012) 2271–2276. [21] H. Nakagawa, S. Nasu, H. Fujii, M. Takahashi, F. Kanamaru, Hyperfine Interactions 69 (1991) 455–458. [22] P. Lukashev, W.R.L. Lambrecht, Phys. Rev. B 70 (2004) 245205. [23] D.E. Jiang, E.A. Carter, Phys. Rev. B 67 (2003) 214103. [24] M. Acet, H. Zähres, E.F. Wassermann, Phys. Rev. B 49 (1994) 6012–6017. [25] C.M. Fang, M.A. van Huis, Heliyon 3 (2017) e00408. [26] S. Nagakura, K. Tanehashi, J. Phys. Soc. Jpn. 25 (1968) 840–846. [27] C.M. Fang, M.A. van Huis, H.W. Zandbergen, Scripta Mater. 64 (2011) 296–299. [28] K.H. Jack, Acta Cryst. 5 (1952) 404–411. [29] E. Fujita, Materia Japan 6 (1967) 647–656. [30] A.D. Le Claire, Diffusion of C, N, and O in metals, in: H. Mehrer (Eds.), Diffusion in festen Metallen und Legierungen, Springer, Berlin, 1990, pp. 481–503. [31] Diffusion in metals, in: E.A. Brandes, G.B. Brook (Eds.), Smithells Metals Reference Book, seventh ed., Butterworth-Heinemann, Oxford, 1992, pp. 13. [32] K.H. Jack, Proc. Roy. Soc. A 208 (1951) 200–215. [33] C. Domain, C.S. Becquart, J. Foct, Phys. Rev. B 69 (2004) 144112. [34] K.-J. Kim, K. Sumiyama, K. Shibata, K. Suzuki, Mater. Sci. Eng. A 181–182 (1994) 1272–1276. [35] Y. Tsunoda, Y. Kurimoto, M. Seto, S. Kitao, Y. Yoda, Phys. Rev. B 66 (2002) 214304.
Fig. 3. Phonon densities of states of γ′-Fe4N and γ-Fe. Open and filled circles indicate experimental results of γ′-Fe4N and γ-Fe, respectively.
(4) The phonon density of states of γ ′ - Fe4 N was reproduced by MD simulations using the new potentials. CRediT authorship contribution statement Katsutoshi Hyodo: Software, Investigation, Writing - original draft, Visualization, Funding acquisition. Shinji Munetoh: Methodology, Software, Resources, Writing - review & editing. Toshihiro Tsuchiyama: Conceptualization, Resources, Writing - review & editing, Project administration. Setsuo Takaki: Supervision, Funding acquisition. Declaration of Competing Interest The authors declare that they have no known competing financial
5