A reference-free MEAM potential for α-Fe and γ-Fe

A reference-free modified embedded atom method (RF-MEAM) potential for iron has been constructed. The new potential is made to predict both bcc and fcc (α-Fe and γ-Fe) lattice properties, with a special interest in modelling in the 800–1300 K temperature range. This is the range in which transformations and key processes in steel occur. RF-MEAM potentials can be used directly in commonly used molecular dynamics simulation software (e.g. LAMMPS). The new potential is compared to several other (M)EAM potentials which are commonly used. It is demonstrated that the new potential combines good characteristics for point defect energies with free surface and stacking fault energies. Also the Nishiyama–Wassermann and Kurdjumov–Sachs orientation relation ratios and interface energies are reproduced, allowing for simulations of α-Fe and γ-Fe interphases.


Introduction
A key component for molecular dynamics (MD) simulations is an interatomic potential capable of predicting experimental material properties accurately. For iron several potentials have been constructed [1][2][3][4][5][6][7][8], several of these potentials were constructed using the Embedded Atom Method (EAM) [1,3,4,7] or modified embedded atom method (MEAM) [5,8] formalisms. More recently machine-learning potentials (MLPs) like Gaussian approximation potentials [9] and neural network potentials [10] have seen increased interest, however existing MLPs are not fit to experimental but rather density functional theory (DFT) data as MLPs require larger datasets to train than * Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. experiments can provide. An issue with DFT is that it finds different material properties for γ-Fe than what is experimentally observed. Additionally MLPs prove impractical for systems containing large numbers (∼10 6 ) of atoms as the MLP formalisms are computationally demanding. Furthermore current MLPs are made for highly specialised applications, e.g. dislocation dynamics [10], whereas the potential in this work is aimed at general simulation at finite temperatures.
The original EAM formalism has some shortcomings, particularly for systems where directional bonding plays a role, e.g. elements from the middle of the transition metal series especially bcc metals such as Fe, and semi-conductors [11]. The MEAM formalism has subsequently been made to address the issues concerning directional bonding. However MEAM potentials are tied to a particular reference structure and its associated zero-temperature Rose universal equation of state (EoS) to find its total energy [12]. Timonova and Thijsse [13] have developed an improved MEAM formalism which they have called reference free (RF)-MEAM. RF-MEAM removes the dependency of the potential on the Rose equation, thus removing the need for a particular reference structure, hence the name. The RF-MEAM formulation is used here to construct an improved interatomic potential for iron, applicable to both the ferrite (bcc) and austenite (fcc) phase.
Many EAM potentials have a similar analytical form as the RF-MEAM potential, e.g. [1,3,4,14,15]. It is therefore convenient to use these EAM potentials as a starting-point for fitting RF-MEAM potentials. The use of a prior EAM potential eases the finding of an acceptable local minimum from the many local minima in the parameter space. Existing EAM potentials can be enhanced by adding angular functions for the local electron density. In this work this method is illustrated for iron in section 2, but it may be applied to other systems for which an EAM potential is available, e.g. vanadium [14], or a system containing both iron and carbon [15].
The widely used EAM potential by Ackland et al [1] was chosen as the point of departure to construct a RF-MEAM potential as formulated by Duff et al [16]. The Ackland potential performs well for predicting several properties of α-Fe, notably various point defects. Simultaneously there is room for improvement in the description of the γ-Fe phase, as the elastic constants and lattice parameter are not well predicted by the potential for this phase. The new potential is fitted to experimental data where possible, as there is no direct experimental data available for the elastic constants of γ-Fe at 0 K an extrapolation of high-temperature neutron scattering data [17] is performed. Details for this extrapolation are provided in section 2.3.
The potential presented in this work performs well for both the α-Fe and γ-Fe lattice parameter, elastic constants, and several point-defects. This is not only an improvement on the EAM potentials, but also on some alternative MEAM potentials [5,8] that describe the formation energies of interstitials [5] and vacancies [8] less well. Results demonstrating the effectiveness of the new potential are presented in section 3. The ability to predict lattice properties for both α-Fe and γ-Fe lattices is particularly useful in modelling transformations and processes in steel, where processes may take place in both crystal structures.
In section 2 the experimentally measured data selected for the fit, in particular the elastic parameters for austenite, and the fitting procedure for the potential are described. In section 3 the results of the potential are compared to several other potentials and its performance is discussed. Finally, in appendix the potential parameters are given.

RF-MEAM formalism
The RF-MEAM formalism in this work follows the form used by Duff et al [16]. The total energy function for a system of N atoms of type ζ is given as Here ϕ ζi,ζj (r ij ) is the pair-potential between atoms i and j, separated by a distance r ij . E emb ζ (ρ i ) is the embedding energy function given as, where a quartic term was included to accommodate the starting potential from Ackland. ρ i is the fictitious local electron density at site i, The density term ρ (0) i is the spherical term present in EAM, the terms ρ (l) i where l = 1, 2, 3 are the angular terms which are added in the MEAM formalism. When t where f (l>0) ζi are partial electron density contributions, and P (l) (x) are Legendre polynomials, given by Rodrigues' formula P (l) (x) = 1 The (partial) electron density and pair-potentials are sums of cubic splines: where H (r) is the Heaviside function. In the fitting procedure a (n,l) i , a ζ , b ζ , c ζ , and d ζ are the parameters to be optimised. In contrast to the MEAM formalism the pairpotential is not tied to the reference structure via the Rose EoS. All fitting in this work was done using MEAMfit v2 [16]; the values of the parameters are given in section 4.
The pair potential as in equation (8) is valid above a cutoff distance (here 1.2 Å), below the cut-off distance the pair potential is divided into two parts. Below 0.9 Å the Biersack-Ziegler-Littmark [18] form is used. Between 0.9 Å and 1.2 Å an interpolation of the form e (B0+B1r+B2r 2 +B3r 3 ) is made by MEAMfit, with B i values such that the function and its first derivative are continuous at 0.9 Å and 1.2 Å. The resulting interpolation function is reconstructed within the LAMMPS implementation for RF-MEAM [19].

Methodology
The fitting of the RF-MEAM potential was performed using MEAMfit v2 [16]. The starting point in this work is the EAM potential by Ackland et al [1], this potential uses cubic splines as in MEAMfit. Since the Ackland potential is effectively an initial guess for the RF-MEAM potential it is important to establish where the prior potential performs well and also where room for improvement exists. The lattice constant for α-Fe is predicted well by the Ackland potential as are pointdefects; changes must not infringe upon the qualities of the existing potential. For γ-Fe several changes need to be made, most notably to the elastic constants and lattice parameter. The chosen fit-data is discussed in sections 2.3-2.6 as well as the detailed work-flow.
The strategy used in this work is as follows: (a) First the lattice parameter, cohesive energy, and bulk modulus are fitted for perfect α-Fe, γ-Fe, and ε-Fe lattices.
Here the angular terms do not contribute, so they are not included in the fit yet. (b) Next, the elastic parameters are fitted against data of strained lattices, at this point the angular terms become relevant and are included in the fitting procedure from this stage onward. Note the EAM portion of the potential are still included as optimisation parameters in the fitting procedure. (c) Once the lattice and elastic parameters are fitted, pointdefect energies are included, see section 2.5. At this stage inter-atomic forces are relevant, these forces ensure that the crystal remains stable when a point-defect is present. Reference data for forces are obtained from Vienna Ab-initio Simulation Package [20][21][22][23][24] (VASP 5.3.5) simulations. (d) Once the lattice parameters, cohesive energies, elastic constants, and point-defects are reproduced satisfactorily the potential can be used in LAMMPS [19] to verify the potential's performance in MD simulations.
The weighting of various fitting targets is chosen such that all the energies, forces, and stresses are equally weighted. This weighting is entered in the fit database used by MEAMfit v2 [16]. It is important to notice that each supercell returns one energy, six stress values, and for each atom in the supercell three forces. During the optimisation some weights are adjusted to correct for poorly fitted quantities, this is repeated until the resulting fit is satisfactory. MEAMfit v2 uses a Normalised Weighted RMSE for the optimisation function.
A notorious issue with iron, and other ferromagnetic systems, is the role of magnetisation. Magnetic ordering plays a key role in the bcc-fcc transition in iron [25]. An in depth study to incorporate magnetisation into an interatomic potential was conducted by Dudarev and Derlet [26]. The potential in that work was fitted to the magnetic ground states for both α-Fe and γ-Fe. Magnetic contributions at finite temperatures are not explicitly included in the potential fitted in this work.  (9) and (10) as found in this work.  (9) assuming c = 0. b Extrapolated to equation (10). c Extrapolated to equation (10), assuming B = 204 at 298 K, C ′ = 0 at 823.

Elastic constants for γ-Fe
The elastic constants for γ-Fe are experimentally determined at high-temperatures or in stabilizing alloys. However, the RF-MEAM potential fitted in this work is fitted to static data, so data at 0 K is needed. In table 1 it can be seen that DFT calculations at high temperature contradict experimental data. Therefore, instead of using DFT, another approach is chosen to approximate the required data using experimental data. In this work an extrapolation was performed using the experimental data of Neuhaus et al [17], as their data follows the expectations in [27] best where it is assumed that the elastic parameters for γ-Fe will not differ greatly from those of α-Fe. The data points were used for extrapolation using equation (9) which follows the experimental data for several metals, but most importantly for α-Fe as shown by Varshni [28] As there are just two datapoints equation (9) cannot be used as-is, hence it was assumed that ∂C ij /∂T = k, where k = ∆C ij /∆T is constant, i.e. c = 0. Such a linear fit is also performed for the bulk modulus and shear modulus in the work by Lindgren and Gyhlesten Back for pure iron, see [29]. A linear fit will produce a value for C ij that is likely too large. Results for extrapolating the elastic constants for γ-Fe from [17], as well as the data from [33][34][35]. The narrow solid line and dashed line (blue) and the + symbols are the linear, quadratic and literature values for C 11 respectively. The thick dash-dotted line and dotted line (red) and the 'x' symbols are the linear, quadratic and literature values for C 12 resp., and the narrow dotted line and dash-dotted line (black) are the quadratic extrapolation for C 11 and C 12 using that C ′ = 0 at 823 K and B = 204 GPa at 0 K respectively. The thick solid line and dashed line (purple) and the ' * ' symbols are the linear, quadratic and literature values for C 44 respectively.
Varshni also introduced an asymptotic approximation to equation (9) [28]: Effectively the extrapolations with equations (9) and (10) yield an upper and a lower value to the elastic constants. In table 1 the results as well as the experimental values from [17] are given. One expects dC ij /dT| T→0 = 0 as there is little expansion near absolute zero, therefore equation (10) provides a more physically consistent temperature dependence than the linear approximation.
As the value for C 12 is not explicitly given in [17] but rather C ′ = (C 11 − C 12 ) /2, an additional datapoint is entered, C ′ = 0 at 823 K [30,31], as this is roughly where the martensitic transformation would occur, where it is assumed that the martensitic transformation is caused by mechanical instability of the γ-Fe lattice, which is why C 11 < C 12 at low temperatures. Here it is used that the paramagnetic γ-Fe phase is unstable at low temperatures. Finally the assumption that B = (C 11 + 2C 12 ) /3 ≈ 204 GPa at 298 K (from [27]) is included, here C 11 and C 12 are solved simultaneously for equation (10) using a least-squares minimisation. Additionally the data is compared to a theoretical calculation based on calculated (ab initio) lattice parameters [32].
The extrapolation with equation (10) using the experimental values from [17,[33][34][35] performs best with the assumption that C ′ = 0 at 823 K, and the prediction of Ghosh and Olson [27] that B = 204 at 298 K. The results can be seen in figure 1. Furthermore, the temperature dependence of the elastic constants is similar to α-Fe as suggested in [27]. The potential is therefore fitted to B = 204 GPa, C 11 = 198 GPa, C 12 = 207 GPa, and C 44 = 114 GPa.

Crystal structure and elastic constants
First the non-angular contributions of the (fictitious) electrondensity are optimised, i.e. the EAM portion of the potential where the angular terms are set to 0, i.e. t (l) i = 0 for l = 1, 2, 3 as in equation (4), as well as the pair-potential terms. Fitting is performed on multiple defect free α-Fe lattices with varying lattice parameters. The lattice parameter is varied around the experimental value a = 2.8605 Å (±3%) [38,39], in [39] this value is found by extrapolation to 0 K using thermal expansion [40] with the data of [41]. The cohesive energy is set to E coh = −4.28 eV/atom [2,42]. For the lattice parameters near the equilibrium value the energy is calculated using VASP 5.3.5 with a Projector Augmented Wave Perdew Burke Ernzerhof potential (PAW-PBE) to establish the shape of the potential. Note that the PAW-PBE lattice parameter is significantly smaller than the experimental value, therefore the VASP lattice parameters were shifted as detailed in table 2.
After the perfect lattice data (for bcc, fcc, and hcp) is fitted, stresses are included. For α-Fe the elastic constants used are C 11 = 243 GPa, C 12 = 138 GPa, and C 44 = 122 GPa [41]. The elastic constants are determined at 4.2 K, here we assume that the temperature dependence near 0 K is negligible as changes are minimal at low temperatures [41]. The same approach is taken regarding the γ-Fe lattice. The target lattice parameter is a = 3.562 Å [38], this value was extrapolated from high temperature measurements and verified using FeNiMn alloys. The energy difference at 0 K isF ∆E bcc→fcc = 0.06 eV derived by thermodynamic assessment [43]. Ab initio studies suggest a value of 0.1 eV [44,45] for anti-ferromagnetic (AFM) γ-Fe at 0 K. Elastic parameters were extrapolated from high temperature measurements as part of this work as detailed in section 2.3. The potential is fitted to B = 204 GPa, C 11 = 198 GPa, C 12 = 207 GPa, and C 44 = 114 GPa.
An additional data point is a ε-Fe lattice, where a = 2.523 Å, c/a = 1.603 [46], the experimental data holds for room-temperature (296 ± 3 K). The values at zero pressure are extrapolations of an empirical function derived in [46]. For the cohesive energy ∆E bcc→hcp = 0.082 eV from ab initio calculations [44] was used as there is no direct experimental data available. The cohesive energy of ε-Fe is similar to that of γ-Fe or lower as found by some studies [2,44].

Point defects
Especially for (relaxed) point defects the RF-MEAM potential improves on the starting (EAM) potential by adding in angular dependencies for the local density function. Using the perfect lattices for α-Fe and γ-Fe using the target lattice parameters as in section 2.4. Several point defects are constructed. For α-Fe: a vacancy; a vacancy cluster of first-nearest neighbours (1NN), and second-nearest neighbours (2NN); common interstitials, <100>, <110>, <111> dumbbells, tetrahedral (TET) and octahedral (OCT). For defects in γ-Fe only a vacancy was included, as there are few experimentally known defect energies and DFT calculations for γ-Fe give results that differ from known experimental values. All the used target values, and their respective sources, are presented in tables 3-5.

Vacancies.
A notable change to most existing potentials is the vacancy formation energy in α-Fe. Following the work of Glensk et al [47] it is assumed that previous methods to estimate the vacancy formation energy have systematically given a value that is too low. Therefore values taken from DFT calculations were used. For the vacancy in γ-Fe a DFT value is used for AFM ordering as it is found that this is the most stable ordering [48]. For the bi-vacancy clusters in α-Fe the binding energies were taken from the detailed work of Domain et al [49]. In this work the indirect binding energy for cells with 128-2 atoms at constant volume was used, as the vacancy clusters are modelled for large lattices at equilibrium.

Interstitials.
In literature various values for interstitial energies, particularly for dumbbells, are presented. In this work the data by Marinica et al [50] were used, as it confirms earlier data [1,49,51,52]. Moreover, the various energies in the references show the same general order from most to least stable interstitial. In [50] the <100> dumbbell is unstable and relaxes to an octahedral site so no value was provided. Other studies [1,49,51,52] show that the formation energy is similar to, or as Marinica suggests higher than that of the octahedral site.
The pair-potential part of the Ackland potential uses a short range cut-off radius of 2.05 Å, this is too large for (unrelaxed) interstitials. Unrelaxed interstitials are included to provide data for interatomic forces, so the cut-off radius must be chosen accordingly. In order to accommodate for this the short range cut-off has been set to 1.2 Å.

Stacking faults and free surfaces (FSs)
Unstable stacking fault (USF) energies and FS energies are not used in the fit, however the performance of the potential presented in this work has been tested for several USFs and FSs. The reference data is mainly taken from DFT calculations. For α-Fe data is available from an experimental study [53], which has been included in table 3. The theoretical data presented for γ-Fe FSs and for α-Fe USFs is taken directly from DFT calculations.
Orientation relations (ORs) between α-Fe and γ-Fe like Nishiyama-Wassermann (NW) and Kurdjumov-Sachs (KS) are also tested, but not fitted to. For the OR's reference data is taken from DFT calculations [54,55].

Structural properties
Physical properties, as calculated with the RF-MEAM potential presented in this work using LAMMPS [19], are presented in tables 3-5. The new potential, given in table 6, is compared to several other potentials, notably the Ackland potential [1] on which it is built, labelled Ack04, the MEAM potential by Lee et al [5] (BJLee), a second MEAM potential by Lee et al [8] (TLee). Also added in the comparison are an analytical bond-order potential by Müller et al [2] (Mül07), another EAM potential by Mendelev et al [4] (Mend03), and finally a recent angular-dependent potential (ADP) by Starikov et al [56] (Stari21). Reference data is provided from various sources. Where possible, experimental data is provided otherwise a theoretical value is used.
The performance of the various listed potentials was checked using LAMMPS. The lattice parameter is determined by isotropic relaxation of the volume of a perfect lattice. For defects, and stacking faults the supercell was relaxed at constant volume. For the FSs the surface direction is allowed to expand to allow for relaxation. Except for the FSs periodic boundary conditions are used. The stopping criteria are set such that the energy change between iterations is 1 × 10 −8 eV, for forces the 2-norm of the global force vector must be less than 1 × 10 −8 eV Å −1 . All supercells were constructed using Atomsk [57]. Defect formation energies and surface energies are calculated from the excess energy using the cohesive energy for the relevant state. The total energy for a supercell containing N atoms in a perfect lattice is known, so the energy difference is caused by the introduced defect. All point defects in α-Fe were simulated in a 1024 atom cubic cell of 8 × 8 × 8 unit cells. For γ-Fe a cubic cell containing 864 atoms, i.e. 6 × 6 × 6 unit cells was used. The FSs have been simulated in a 6 × 6 × 6 unit cell system with periodic boundaries in two directions for both α-Fe and γ-Fe, leaving two 6 × 6 unit cell FSs. For the 6 × 6 × 6 unit cell the   separation of the two FSs is more than three times the maximum cut-off radius used in the potential. The FSs were also simulated for a smaller supercell, where the difference in FS energy with a 4 × 4 × 4 unit cell is in the order of 0.01 J m −2 .
Finally the USFs had an area of 6 × 6 unit cells, and the cell length perpendicular to the stacking fault is 16 unit cells, i.e. 8 unit cells between stacking faults, here periodic boundaries were used. The separation of the stacking fault interfaces is more than four times the maximum cut-off radius used in the potential. The lattice parameters and elastic constants found for the RF-MEAM potential reproduce the experimental values for both α-Fe at 0 K and our extrapolated experimental elastic constants for γ-Fe at 0 K. This in contrast to most other potentials as shown in table 3, except for the potential by Lee et al [8]. In general the γ-Fe lattice parameter is overestimated, except for Stari21 which was fitted to DFT data. For Stari21 the a fcc,0K : a bcc,0K ratio agrees with the experimental ratio. However the elastic parameters for Stari21 are slightly off for α-Fe. The RF-MEAM potential is among the potentials (BJLee, and Stari21) that find a high bulk modulus for γ-Fe, as was predicted by Ghosh and Olson [27].

Point defects
For point defects there are greater differences between the various potentials. TLee shows values for vacancy clusters in α-Fe, vacancy migration energies which are too low. Stari21 performs better, however the interstitial formation energies are relatively high. This makes the current RF-MEAM potential more suitable for modelling and simulating defects, most notably vacancy (migration) in α-Fe and γ-Fe than TLee. Note that the RF-MEAM potential performs best when combining vacancies and interstitials in a single system, as some other potentials perform well for specific types of point defects. The nearest-neighbour vacancy migration path is presented in figure 2, note the unphysical behaviour observed for the TLee potential.

Surfaces and stacking faults
The surface energies and USFs which were not part of the fitting procedure, however they are reproduced well for the RF-MEAM potential (table 5). Note, there exist relatively large differences between the various potentials and even some literature sources. Here it is seen that, except for BJLee and Stari21, other potentials underestimate surface and USF energies. Stari21 on the other hand finds relatively high energies.
Lastly the NW and KS-OR are evaluated with the new potential. The energies for both OR predicted by the RF-MEAM potential agree with literature [54,55]. The lattice ratios for α-Fe and γ-Fe follow the experimental data from [38]. This is contrasted by Mend03 and Stari21, which both find very high interface energies for the fcc-bcc ORs. These are valuable characteristics of the RF-MEAM potential when applied to interphase modelling. Hence the RF-MEAM potential is well-suited for static interphase modelling.
Finally for the common {110}-plane and {112}-plane general stacking faults (GSF) in bcc, when shifting in the <111>direction are simulated. The supercell used is 6 × 6 × 16 bcc unit cells in size with periodic boundary conditions. The stopping criteria are the same as for other simulations. The stacking fault is simulated by incrementally moving one half of the supercell with respect to the other along the <111>-direction.  Relaxation of atoms is only allowed perpendicular to the stacking fault plane, as in [56,69].
In figure 3 the RF-MEAM potential is compared to BJLee, TLee, and Stari21. The GSF energies for both planes are similar for the four tested potentials, only Stari21 produces higher values. The higher values correspond to DFT results [69], to which Stari21 was originally fitted.

Thermal expansion
For thermal expansion four perfect lattices (5 × 5 × 5 unit cells) with different lattice parameters are heated using the canonical (NVT) ensemble setting in LAMMPS. The cell is simulated for 2.5 ns, such that the equilibrium pressure can be determined. The equilibrium lattice parameter is determined by using that P = −dE/dV. A linear fit to the equilibrium pressure with respect to volume is made, for the equilibrium lattice parameter P (V) = 0.
At 300 K (room temperature) the calculated lattice parameter in α-Fe is 2.864 Å, and at 1100 K the α-Fe lattice parameter is predicted to be 2.886 Å, which is in agreement with the results from Acet [38]. Moreover the expansion in the temperature range 800-1300 K, shown in figure 4, follows the values presented in [38].
The thermal expansion of the RF-MEAM potential differs from other potentials. For α-Fe the growth is repressed particularly for temperatures below 700 K, and for γ-Fe the lattice is not stable at low temperatures. This instability needs to be  [66] is used, for fcc free surfaces [67], and for the ORs [55]. taken into account when using the potential in finite temperature simulations for γ-Fe. The instability makes it difficult to accurately determine the lattice parameter at 700 K, so there is a jump in figure 4.
Most other potentials show only linear growth for the γ-Fe lattice parameter, however this is not physical. At low temperatures the thermal expansion should be small and gradually increase to a linear growth as is seen for the literature data. The RF-MEAM potential and Stari21 display this for both α-Fe and γ-Fe, where TLee only has linear expansion. A peculiarity of Stari21 is that at room temperature the equilibrium lattice parameter has decreased which is unphysical behaviour.

Discussion
Most Fe-potentials are constructed with specific targets in mind; the TLee potential is made to model the bcc to fcc transition, but it fails when modelling vacancy structures. The ADP potential Stari21 is made with dislocations and point defects in mind, but is not well-suited for surfaces and stacking faults, or the fcc-bcc interface. The RF-MEAM potential presented here was fitted to simulate experimental data where possible. However, this is challenging as a wealth of fitting data, such as atomic forces, are available from DFT only. For the general structural data fitting to experimental data at 0 K was accomplished, however the elastic constants of γ-Fe at 0 K are fitted to extrapolated experimental data.
The new RF-MEAM potential can describe vacancy systems well, it also follows the FS energies, and fcc-bcc OR energies. Hence, the RF-MEAM potential is most suitable when simulating vacancy systems and surfaces or stacking faults, and particularly the interphase relation between α-Fe and γ-Fe. It must be noted that static simulations here work better than finite temperature simulations as the thermal expansion of both phases is slightly off. In figure 5, the specific strengths of the listed potentials in comparison to the RF-MEAM potential are visualised.

Conclusions
In this work a RF-MEAM potential for iron has been obtained for the purpose of simulating both α-Fe and γ-Fe. The presented potential reproduces many experimental material properties, such as defects, interfaces or surfaces, better than prior (M)EAM potentials. The current RF-MEAM potential is ready to for use in LAMMPS for calculations on large systems, for which MLP's are still impractical.
The potential presented in this work predicts accurate lattice properties and defect energies, most notably for vacancies, and self-diffusion. Additionally surface energies and stacking fault energies are predicted well. ORs between the bcc and fcc lattice are predicted well for both NW and KS. The thermal expansion in the 800-1300 K temperature range follows the experimental data well. Therefore the potential presented here can be used for simulations in α-Fe, and in γ-Fe, as well as in inter phase simulations.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).

Funding
This research was carried out under project number T17019m in the framework of the Research Program of the Materials Innovation Institute (M2i) supported by the Dutch government.

Appendix. The potential
In table 6 the parameters of the RF-MEAM potential presented in this work are given.
Fe,Fe , s  The potential parameters are labelled as in the equations (1)-(8) in section 2.1.