Interatomic potential for Fe–Cr–Ni–N system based on the second nearest-neighbor modified embedded-atom method

Abstract The interatomic potential for Fe–Cr–Ni–N system based on the second nearest-neighbour modified embedded-atom method has been developed in this work. The potential is based on those for the corresponding lower order systems. The potential parameters for the binary systems, Cr–N, Ni–N, Ni–Fe and Ni–Cr, were determined by fitting the lattice constants, elastic properties, heat of solution and defect binding energies. The potential parameters for the ternary systems were calculated based on the corresponding binary systems. Then, all of them were applied to the quaternary system Fe–Cr–Ni–N to confirm their validity by a simulation of the lattice constants of AISI 316 austenitic stainless steel with a range of nitrogen content. The results were in good agreement with the previous observations and calculations.


Introduction
The strengthening of Fe-Cr-Ni austenitic stainless steel by nitriding has long been recognised. [1] Above the solubility limit of nitrogen, the formation of nitrides can play a pivotal role in improvement of the strength and hardness. Below the solubility limit of nitrogen, the improvement is responsible for the presence of N atoms in interstitial sites, which builds strong interactions with the lattice defects present in the steel and has a significant influence on the steel properties.
Atomistic level computer simulation is an ideal tool to understand the effects of N on the atomic structure and the strengthening. However, it requires the knowledge of interaction forces among Fe, Ni, Cr and N atoms. Several (semi-)empirical interatomic potentials have been devoted to the Fe-Cr-Ni-N system or its subsystems. [2][3][4][5][6] Particularly, the embedded-atom method (EAM) potential for the Fe-Cr-Ni-N quaternary system has been developed by Grujicic and Zhou [6], since it is widely accepted that the empirical interatomic potential based on EAM is suitable for metals. However, it is difficult to describe correctly some properties of the quaternary system using the EAM potential. For example, face-centred cubic (fcc) Fe is more stable than body-centred cubic (bcc) Fe, and the calculated heat of solution of Cr in Ni is opposite to the relevant experiments using this potential. [6] The second nearest-neighbour modified EAM (2NN MEAM) can resolve the above-mentioned problem of EAM by modifying the EAM formalism, because it considers not only the first nearest-neighbour interactions but also partially the second nearest-neighbour interactions. [7] Various interatomic potentials have been developed using this potential formalism. In addition to the unary systems Fe, Cr, Ni and N, [7][8][9][10] the 2NN MEAM has been applied to describe the binary systems Fe-Cr [11] and Fe-N [8]. However, the parameters for the quaternary system Fe-Cr-Ni-N remain lacking.
The objective of the present work is to develop the potential parameters for Cr-N, Ni-N, Ni-Fe and Ni-Cr binary systems, and the interatomic potential parameters for the Fe-Cr-Ni-N quaternary system based on the corresponding lower order systems. The reliability of these potential parameters is then examined by computing the physical properties and comparing them with the available observation and calculation results.

Formalism
The details of the 2NN MEAM formalism have been reported elsewhere [7,9], and only the central formulas are summarised here. In the 2NN MEAM, the total energy of a system E tot can be calculated by where E i is the total energy of an atom, F i is the embedding energy of atom i embedded in a background electron density ρ i , ϕ(r ij ) is the repulsive pair wise potential and S ij is the screen function between atoms i and j with a distance r ij . In Equations (4)-(6), d is an adjustable parameter, E c is the cohesive energy, r e is the first nearest-neighbour distance in the equilibrium reference structure, B is the bulk modulus and Ω is the equilibrium atomic volume of the reference structure. The parameters E c , r e (or Ω), α (or B) and d can be determined by fitting the properties of the reference structures.

Determination of potential parameters
As mentioned above, the potential parameters for Fe, Cr, Ni and N unary systems, as well as Fe-N and Fe-Cr binary systems, are available from Lee et al. [7][8][9][10][11]. In the present work, these potential parameters were taken without any modification, as listed in Tables 1 and 2. The potential parameters for Cr-N, Ni-N, Ni-Fe and Ni-Cr binary systems were required before determination of the parameters for higher order systems.
From the 2NN MEAM formulas, [7] the main potential parameters are E c , r e , α, ρ 0 , d, C min and C max . The first three are given experimentally for the reference structures are real phases. For the fictitious NaCl-type NiN and L1 2 -type Ni 3 Cr, the first three parameters have been determined by first-principle calculations, [6,[13][14][15][16][17][18] except the α of Ni 3 Cr was determined by fitting the elastic constants of Ni-Cr alloy in this work. ρ 0 is the ratio between the electron density scaling factors of the corresponding unary elements, and d is an average value of those for pure elements.
C min (i−k−j) and C max (i−k−j) are parameters to estimate the screening factor of atom k to the interaction between two neighbouring atoms i and j. In the case of binary system, there are four cases that one of the interacting atoms and/or the screening atoms can be different types, i.e. i−k−j = A-B-A, B-A-B, A-A-B, A-B-B, where A and B denote the components of the binary system. As discussed by Lee and Baskes [7], the key difference of the 2NN MEAM from the 1NN MEAM is that the second nearest-neighbour interactions are partially considered during the procedure for the determination of ϕ(r) values. Moreover, the third nearest-neighbour interactions (even more) have to be completely screened. Otherwise, it will cause an inconsistency and error in energy calculations, because only the interactions between the first and the second nearest neighbours are considered in the definition of the pair interaction in Equation (3). The range of C min and C max depends on the character of the reference structure. In the present work, the range of C min for L1 2 -type and NaCl-type reference structure was from 0 to 1, and the effect of such third nearest-neighbour interactions was minimised using a radial cut-off function with a cut-off distance between second In order to describe binary systems, the pair interaction ϕ(r ij ) between different elements is determined by estimating the total energy per atom of a reference structure using the zerotemperature universal equation of state. [12] Then, the value of the pair interaction ϕ(r ij ) can be evaluated from the known values of the total energy per atom and the embedding energy, as a function of the nearest-neighbour distance. In the 2NN MEAM, a perfectly ordered binary intermetallic compound is considered as a reference structure. In this compound, the first nearest-neighbour atoms with respect to any atom are the same type but different from the atom, while the second nearest-neighbour atoms and the atom are the same type. For the Cr-N and Ni-N systems, NaCl-type CrN and NiN-ordered structures are chosen as the reference structures. For the Ni-Fe and Ni-Cr systems, L1 2 -type Ni 3 Fe and Ni 3 Cr-ordered structures are chosen as the reference structures. The Ni 3 Fe and CrN are real compounds, while the fictitious NaCl-type NiN and L1 2 -type Ni 3 Cr are constructed for the fitting purpose. In the L1 2 -type Ni 3 Fe structure, for example, the total energy per atom (for 3/4 Ni atom + 1/4 Fe atom) is given by where Z 1 and Z 2 are the numbers of first nearest-neighbour and second nearest-neighbour atoms in the L1 2 -type Ni 3 Fe structure, respectively. In this case, Z 1 and Z 2 are 12 and 6, respectively. a′ is the ratio between the second and first nearest-neighbour distances. S Ni and S Fe are the screening factors for second nearest-neighbour interactions between Ni atoms and between Fe atoms, respectively. They can be calculated by the screening parameters, C min and C max . Therefore, the pair interaction between Ni and Fe can be deduced in the following form: The value of E u Ni3Fe (r) can be obtained from the universal equation of state [12]: where and (2) Table 1. Potential parameters for Fe, cr, ni, and n. the units of the cohesive energy E c , and the equilibrium nearest-neighbor distance r e are eV and Å, respectively. a ref. [8],, b ref. [9],, c ref. [10],, d ref. [7]. and third nearest-neighbour distances. Setting the C max to 2.8 could ensure the first nearest neighbours of L1 2 -type and NaCltype structure were completely unscreened for reasonably large thermal vibration. The potential parameters for the ternary systems are obtained by a combination of all the unary and binary parameters. Here, only the C min (i−k−j) and C max (i−k−j) values are required to be determined. There are three additional cases, i.e. i−k−j = A-B-C, B-C-A and C-A-B, where A, B and C denote the components of the ternary systems. When no information on physical or thermodynamic properties of the ternary systems is available for parameter optimisation, the assumption raised by Lee [19] can be used: the degree of screening by a third atom to the interaction between the other two different type atoms is averaged by where C means C min or C max . The final 2NN MEAM parameters for the ternary systems are listed in Table 3.
As for the quaternary system, there is no specific parameter, since only three atoms are involved in the screening consideration Table 2. Potential parameters for the binary systems. the units of the cohesive energy E c , and the equilibrium nearest-neighbor distance r e are eV and Å, respectively. a ref. [11], b ref. [7].  Table 3. Potential parameters of C min and C max for ni-Fe-cr, Fe-cr-n, ni-Fe-n and ni-cr-n ternary systems.   [43], c ref. [44], d ref. [45], e ref. [13], f ref. [46], g ref. [47], h ref. [23], i ref. [24], j ref. [14], k ref. [38], l ref. [6], m ref. [15], n ref. [16], o ref. [17], p ref. [18], q ref. [20], r ref. [48], s ref. [49], t ref. [25], u ref. [21], v ref. [4], w ref. [50].  [20] was employed in all calculations. It is noted that the meanings of several parameters in the 2NN MEAM potentials, such as C(i−k−j), d, etc., are different from those in LAMMPS, and the potentials can work correctly only when special options are used. An example of LAMMPS input script for the lattice constant calculation of the Fe-Cr-Ni-N alloy system is provided in the supplementary material. All calculations were performed at 0 K allowing full relaxations of individual atoms. Periodic boundary conditions were applied in all three spatial directions. For the binary alloys, the properties calculated in the present study can be divided into two groups. The first group is the properties which are used for parameter optimisation. The lattice constants, bulk modulus, and cohesive energy of the reference structures and the dilute heat of solution, belong to this group. The comparisons between the current calculation results and the previous results of experiments and/or calculations show the quality of fitting. The second group is the properties from the previous experiments and/or calculations, but those are not used for parameter optimisation. The lattice constants and cohesive energy of random solutions, the elastic constants and the binding energies of pairs of defects, which are available from previous literatures, belong to this group. Through these comparisons with previous results, the transferability and the applicability of the present potentials are shown. In the calculations for binary alloys, the supercell sizes were 10 × 10 × 10 unit cells for solid solutions in the 2NN MEAM. In other words, the simulation of the quaternary system is performed by a combination of the parameters for the lower order systems.

Calculation of physical properties
The potential parameters determined by the above procedure were then used to compute the properties in order to evaluate   [24], c ref. [26], d ref. [27], e ref. [25], f ref. [28], g ref. [29], h ref. [30], i ref. [31], j ref. [32], k ref. [33], l ref. [34], m ref. [35], n ref. [36], o ref. [37]. with bcc, fcc and hcp structures. For NaCl-type and L1 2 -type reference structures, the supercell sizes were 5 × 5 × 5 unit cells and 10 × 10 × 10 unit cells, respectively. It was confirmed that all calculation results were independent of the selected supercell sizes. To simulate the solid solution, Fe, Cr or Ni atoms were randomly distributed on the lattice sites. The radial cut-off distance here was chosen 3.6 Å for bcc structures, 4.5 Å for fcc structures and 4 Å for hcp structures. As shown in Table 4, the properties of the first group are correctly reproduced. The errors of most results are within 0.5%, except those of dilute heat of solution. In spite of this, the most values of dilute heat of solution are more accurate than the previous EAM results. [6] The lattice constants and the cohesive energies of random solutions for Ni-Cr and the cohesive energy for Ni-Fe alloys, are then calculated, which are shown in Figures 1 and 2. For Ni-Cr alloy (Figure 1), comparing with the experimental results, [4,21] the lattice constants are underestimated less than 0.42%, and the cohesive energies are underestimated less than 0.19%. For Table 6. calculated solute-solute and vacancy-solute binding energies (E b in eV) of ni-Fe system in first (1nn) and second (2nn) nearest neighbor using the present 2nn MeaM parameters, in comparison with the previous calculations.  Fe, Cr or Ni atoms on the lattice sites, N atoms were randomly distributed in the octahedral interstices. The radial cut-off distance here was chosen 5 Å due to the lattice expansion with the increase in the N content. The calculation results are compared with the previous experimental data of gaseous nitriding of AISI 316 [40,41] in Figure 4. It depicts a reasonable agreement between the present calculation and literature information on the lattice constants of AISI 316 austenitic stainless steel at various nitrogen contents.

Conclusions
The 2NN MEAM interatomic potential for Fe-Cr-Ni-N system has been developed in this work. The potential is based on those for the corresponding lower order systems. It has been shown that (i) the 2NN MEAM potentials for the Cr-N, Ni-N, Ni-Fe and Ni-Cr binary systems can reproduce various fundamental physical properties; (ii) when these potentials are applied to the Fe-Cr-Ni-N system, they can reproduce the lattice constants of AISI 316 austenitic stainless steel when the N content varies from 0 to 40 at. %. Therefore, the present 2NN MEAM potential parameters have relatively high transferability and applicability for the Fe-Cr-Ni-N system and the lower order systems.

Disclosure statement
No potential conflict of interest was reported by the authors. the Ni-Fe alloy (Figure 2), both the cohesive energies of the bcc and fcc phases are calculated. The bcc-fcc transition with this potential occurs at about 24 at. % Ni, which is consistent with the previous results. [22] The calculated properties of the second group are compared with the previous experimental data and first-principle calculations [23][24][25][26][27][28][29][30][31][32][33][34][35][36][37] in Table 5. These properties are correctly reproduced using the 2NN MEAM potential, and the errors are within 2%, except the errors of C 12 for CrN, C 44 for Ni-19.6 at% Cr and B for Cr 2 N are relatively large.

Funding
The binding energies of defect pairs are compared with the previous EAM calculations, angle-dependent potential (ADP) calculations and first-principle calculations [2,38] in Table 6. These EAM potentials were developed by Meyer [39] and Bonny [2], and ADP potential was developed by Mishin [28]. Some obvious deviation appeared in part of the results using these previous potentials. For the 2NN MEAM potential, almost all of the results are located in the range of the first-principle data [2], except the vacancy-solute binding energies of Ni-Fe system in second nearest-neighbour (E V−Ni b (2nn)). For the ternary Fe-Cr-Ni system, the cohesive energies of both the bcc and fcc phases as a function of Ni content are calculated and shown in Figure 3. The supercell size and the radial cut-off were the same as those for the binary alloy systems. The Fe content was fixed at 10 at.%, because the phase diagram at this composition merely shows a bcc-fcc transition. It can be seen that the transition occurs at about 15 at. % Ni, which is consistent with the Fe-Cr-Ni phase diagram. The elastic constants and bulk modulus of several fcc Fe-based alloys are calculated and listed in Table 7. The C 11 is correctly reproduced with this potential, and the error is within 10%. The maximum error of B is about 20%, the maximum error of C 12 is about 30%. However, the C 44 is underestimated by about 50%, which should be kept in mind in future applications of this potential.
For quaternary system, only the lattice constants of Fe-18Cr-12Ni alloy system with a range of N content are calculated as a means to confirm the applicability of the present potential, since there are few of available experimental data for the quaternary system. The contents of Cr and Ni in this alloy system are the same as those in AISI 316. In the calculation, the supercell size was 30 × 30 × 30 unit cells. Besides the randomly distributed