Dynamical stability of radiation-induced C15 clusters in iron

Density functional theory predicts clusters in the form of the C15 Laves phase to be the most stable cluster of selfinterstitials in iron at small sizes. The C15 clusters can form as a result of irradiation, but their prevalence and survival in harsh irradiation conditions have not been thoroughly studied. Using a new bond-order potential optimised for molecular dynamics simulations of radiation damage, we explore the dynamical stability of the C15 clusters in iron under irradiation conditions. We find that small C15 clusters make up 5–20% of the interstitial clusters formed directly in cascades. In continuous irradiation, C15 clusters are frequently formed, after which they remain highly stable and grow by absorbing nearby single interstitial atoms. Growth of C15 clusters ultimately leads to collapse into dislocation loops, most frequently into 1/2〈1 1 1〉 loops and only rarely collapsing into 〈1 0 0〉 loops at low temperatures. The population, size, and collapse of C15 clusters during continuous irradiation correlates well with their formation energies relative to dislocation loops calculated at zero Kelvin.


Introduction
Irradiation of iron and iron-based alloys results in accumulation of defects and defect clusters. Self-interstitial atoms (SIAs) and vacancies cluster together as dislocation loops, voids, and other defect clusters. Unique to iron, small SIA clusters are most stable in the form of agglomerates of C15 Laves crystals, coherently embedded in the bcc structure as predicted by density functional theory [1]. Molecular dynamics (MD) simulations show that C15 clusters can form in radiation-induced collision cascades [1,2] or by direct clustering of migrating SIAs [3], similar to the conditions during electron irradiation. The presence of these highly stable and stationary three-dimensional clusters sets iron apart from non-magnetic bcc metals, like tungsten, where rapidly migrating 1/2 1 1 1 loops are the dominant radiation-induced clusters [4][5][6][7]. At larger cluster sizes, however, dislocation loops eventually become the most stable configuration in iron too [8]. The small size of the C15 clusters make them difficult to observe experimentally, and the first experimental measurements are yet to be published. Hence, the presence and effects of C15 clusters on the microstructural evolution of iron under irradiation still remains unclear.
Atomistic simulations can be used to provide useful predictions of the formation, stability, and evolution of defect clusters, such as the C15 clusters. Even though the relative stability at zero temperature of various defect clus-ters can be estimated from density functional theory, observing the formation and evolution of C15 clusters requires larger-scale classical MD simulations. One critical restriction is, however, the lack of interatomic potentials that correctly predict the C15 clusters as the most stable interstitial-rich clusters in iron. To date, only the embedded atom method potentials by Marinica et al. [1,9] correctly stabilise the C15 cluster in relation to parallel 1 1 1 interstitial configurations. Nevertheless, these potentials also incorrectly predicts 1 0 0 dislocation loops to become more stable than 1/2 1 1 1 loops at sizes larger than around 80 SIAs [9]. This makes it difficult to assess the reliability of simulated cluster transformations, for example transformations from C15 to dislocation loops due to cascade overlap [7], or collapse of growing C15 clusters into a dislocation loop of either type [3,10].
The aim of this article is to study the formation, stability, growth, and collapse of radiation-induced C15 clusters in iron. To achieve this, we first develop a new Tersofftype analytical bond-order potential that, unlike all previous potentials for radiation damage, predicts the correct relative stability between C15 clusters and 1/2 1 1 1 and 1 0 0 dislocation loops across all sizes. We use the new potential to investigate C15 clusters in iron and discuss our results in relation to previous observations and predictions by other interatomic potentials.
Müller et al. [11] with the short-range addition by Björkas and Nordlund [12]. The goal was to adjust the parameters in order to improve the energetics of single interstitials and vacancies, including migration energies, and to correctly stabilise small C15 clusters relative to dislocation loops in accordance with DFT predictions [1,8]. To reproduce the stability of C15 clusters, we found it helpful to monitor the difference in cohesive energy and lattice mismatch between bulk C15 and bcc iron [13]. The shortrange connection to the ZBL potential was also optimised to better reproduce the many-body repulsion at intermediate interatomic distances. A more detailed description of the fitting strategy along with extensive benchmarking of the optimised ABOP is provided in the Supplementary material available online.

Molecular statics and dynamics simulations
Static energy minimisations and phonon calculations were carried out using lammps [14] within the Atomic Simulation Environment (ASE) [15]. For molecular dynamics simulations, we used the parcas code [16,17].
The formation energies of SIA clusters were calculated by minimsing the positions and pressure of a bcc system with around 54 000 atoms. The collision cascade simulations were carried under identical conditions as outlined in [2] (from which the simulation data obtained with the AM04 [18] and M07-B [1,2,9] embedded atom method potentials were taken and further analysed in this work). All cells were quenched to zero Kelvin and an overall zero pressure over a period of 5 ps, to remove thermal displacements before the detailed analysis of the defect clusters. Wigner-Seitz analysis was used to isolate the interstitials and vacancies. Following that, all C15-like clusters were automatically identified by analysing the bond angles and planes of all interstitial clusters (after isolating the interstitial dumbbells using the Wigner-Seitz analysis, C15 clusters contain connected dumbbells with 60-and 120degree angles on 1 1 1 planes). This method of identifying C15-like clusters was found to be sufficiently robust when analysing highly damaged systems [2] as well as large amounts of overlapping cascade data [7].
In the C15 growth simulations, an initial C15 cluster containing 17 SIAs was coherently inserted in a bcc system of 128 000 atoms and relaxed at zero pressure at 500 K. Dumbbell self-interstitials in a randomly sampled 1 1 0 direction were then inserted at intervals of 200 ps. The SIAs were placed at a random lattice site roughly 5-10 A from the C15-bcc interface. After a new interstitial was added, the system was allowed to freely evolve in the N V E ensemble. The simulation approach is similar to the C15 growth simulations carried out in Ref. [10]. During the 200 ps, the interstitial was likely to migrate towards the C15 cluster and eventually attach itself to it, resulting in a slowly growing C15 cluster. The simulations continued until the system contained 150 SIAs (corresponding to roughly 25 ns), before which the C15 cluster in almost  all cases had collapsed or partially transformed into a dislocation loop. A total of 30 simulations were carried out in the ABOP to gather statistics. In addition, 20 growth simulations were performed with the M07-B potential for comparison.

Validation of the new interatomic potential
Before using the new interatomic potential to study the formation and growth of C15 clusters, we here shortly demonstrate its applicability by calculating properties relevant for simulations of radiation-induced defect clusters. In addition to the results presented below, the ABOP describes well the elastic and vibrational properties of bcc iron and the threshold displacement energies, making it suitable for collision cascade simulations. These and additional validation results are presented in the Supplementary material. We compare the results with two embedded atom method (EAM) potentials (AM04 [18] and M07-B [1,2,9]).  During the refitting of the ABOP, we made sure that the basic properties of iron in the bulk C15 phase given by DFT are well reproduced. Tab. 1 shows the cohesive energies, lattice constants, and bulk moduli of the bcc ground state compared to bulk C15 iron. Additionally, we calculated the phonon dispersion of bulk C15 and compared with DFT results from Ref. [13], as shown in Fig. 1. The fairly complex phonon dispersion with up to 18 branches is reproduced by the ABOP in overall good agreement with DFT, providing confidence that a realistic dynamical behaviour of large C15 clusters in α-iron can be expected. In contrast, we found that bulk C15 in the EAM potentials is dynamically more unstable, and completely unstable in the AM04 potential (as evidenced by the appearance of acoustic phonon branches with imaginary frequencies).
For self-interstitial clusters in bcc iron, DFT predicts C15 clusters to energetically favoured over dislocation loops at small cluster sizes. At some critical size, 1/2 1 1 1 dislocation loops become lower in energy. At large cluster sizes, 1/2 1 1 1 are lowest in energy, followed by the 1 0 0 loops and C15 clusters. This stability trend is qualitatively reproduced by the ABOP, as seen in Fig. 2 where formation energies of the three cluster types are plotted as functions of cluster size. However, like all EAM potentials, the crossovers in stability between the cluster types are underestimated. This is illustrated in Fig. 3, where the difference in formation energy between C15 and 1/2 1 1 1 loops, and between 1 0 0 and 1/2 1 1 1 loops are plotted for the ABOP, the two EAM potentials, and the DFT-based scaling laws from Ref. [8]. The ABOP and EAM curves are plotted based on fits to the same scaling laws as the DFT data. As is clear from Fig. 3, the ABOP presents a critical improvement over the EAM potentials in terms of the relative stability of the two dislocation loops. 1 0 0 loops are energetically very close to 1/2 1 1 1 loops in the EAM potentials, with the M07-B potential even predicting a crossover at around 80 SIAs, above which the 1 0 0 loops are incorrectly lower in energy. Only the ABOP reproduces a clear difference in energy between the two dislocation loops. Nevertheless, it is worth noting that since theoretical work has shown that the increasingly anisotropic elasticity at higher temperatures eventually leads to 1 0 0 loops being favoured [20], the EAM potentials can provide a qualitatively more realistic description at high temperatures, while the ABOP can be seen as more accurate at low temperatures.

C15 cluster statistics in single cascades
Having validated the revised ABOP for properties relevant for radiation damage, we use the potential to investigate the formation of C15 clusters in single collision cascades. In total, 150 cascades with the PKA energies 10, 20, and 50 keV were carried out (50 simulations per energy). Results from the EAM potentials were taken from Ref. [2] and further analysed. 10-18% of the interstitial clusters were identified as C15-like clusters in the ABOP, 5-10% in the M07-B potential, and 5-8% in the AM04 potential, with no clear dependence on the PKA energy. Only clusters containing three or more interstitials were analysed (note however that the smallest C15 clusters contains only a net amount of two interstitials, but shows up as six interstitials and four vacancies after the Wigner-Seitz analysis, and were therefore also considered in the analysis). Previous simulations using the M07 potential reported that about 5% of the clusters were C15-like [1], which is in line with our results with the M07-B potential (which only differ from the M07 potential at short interatomic distances [2]). Similar fractions can be expected at higher PKA energies, as cascade-splitting then becomes Potential   increasingly more likely and individual sub-cascades lead to similar cluster statistics.

Growth and collapse of C15 clusters
Previous studies have suggested that once C15 clusters are formed, they can grow by absorbing migrating interstitials [1,3,10]. Zhang et al. used the M07 EAM potential to study the growth of C15 clusters by inserting single SIAs close to the cluster [10]. They observed that the growing C15 cluster eventually collapses into dislocation loops at sizes in the range 56-147 SIAs, with the most likely product being 1 0 0 loops. Similarly, Chartier et al. recently observed, using the same interatomic potential, that the product of a collapsing C15 cluster depends on its size. A larger C15 cluster was seen to nucleate both 1 0 0 and 1/2 1 1 1 loops [3]. As noted previously, however, the M07 potential incorrectly predicts the 1 0 0 loops to be lower in energy at sizes above roughly 80 SIAs. Hence, it remains unclear whether collapse into 1 0 0 loops is the most probable transformation path of large C15 clusters, or if these observations were due to an artefact of the used interatomic potential. We therefore carried out similar growth simulations using the ABOP, which gives the correct relative stabilities of the dislocation loops at all sizes.
In line with the results by Zhang et al. [10], we found that the C15 clusters grow by capturing the nearby SIAs, and eventually become dynamically unstable and collapse into dislocation loops. Table. 2 summarises the statistics of the resulting dislocation loops. The vast majority of the collapsing C15 clusters transform into 1/2 1 1 1 loops in the ABOP. In contrast, the M07-B potential produces both loop types with roughly equal probability. The size at which the transformation into dislocation loops occurs, is highly stochastic. In the ABOP, in a few cases the C15 cluster collapsed into a loop already at 40-50 SIAs, while in some cases dislocation segments did not appear until cluster sizes above 120 SIAs. A similar size range was also observed for the M07-B potential, although the stochastic difference was smaller, with the transformation taking place at around 70-80 SIAs in the majority of the cases. In this size range, the 1 0 0 and 1/2 1 1 1 loops have almost identical formation energies in the M07-B potential (Fig. 3), which corresponds well to the observation that the probability of the C15 cluster transforming into either loop type is close to 50%. In contrast, the ABOP favours 1/2 1 1 1 loops at all sizes, which leads to a very low probability of forming 1 0 0 loops from collapsing C15 clusters. Hence, the most likely outcome of collapsing C15 clusters directly correlates with the zero-kelvin formation energies of the two dislocation loops. This indicates that collapsing C15 clusters does not follow a transformation path that favours either loop type, and that instead it will in most cases simply restructure into the lowest-energy configuration. In a few cases, the transformation into dislocation loops did not complete, and the cluster remained as a dislocation segment terminating in a small stable C15 cluster even after growth to cluster sizes of 150 SIAs (these are labelled "mixed" in Tab. 2). Fig. 4 shows snapshots of the two possible transformation paths into dislocation loops. Interestingly, in the cases that resulted in 1/2 1 1 1 loops, the C15 cluster rarely transformed entirely into a loop. Instead, only part of the C15 cluster collapsed into a 1/2 1 1 1 loop, which is then eventually detached from the cluster while leaving a small C15 cluster behind as the loop rapidly migrates away. The remaining small C15 cluster could then continue growing, and eventually reach an unstable size again. In other words, collapse of growing C15 clusters does not necessarily mean that the C15 clusters will vanish. Instead, they can act as stationary nucleation sites and emitters of 1/2 1 1 1 loops.
In the ABOP, only one case out of 30 resulted in a 1 0 0 loop. Our results therefore confirm that C15 clusters can collapse into 1 0 0 loops, but that it is not very likely at low temperatures. However, as mentioned previously, it is known from both theoretical work [20] and experimental observations [21] that 1 0 0 loops become more stable as the temperature increases. It is, therefore, reasonable to assume that the frequent formation of 1 0 0 loops from unstable C15 clusters in the M07-B potential is more representative of high-temperature (> 600 K [20]) irradiation (regardless of the temperature used in the simulations), while the results from the ABOP can be assumed to be more valid for lower temperatures.
The low probability of collapse into 1 0 0 loops predicted by the ABOP can be justified based on recent Object kinetic Monte Carlo (OKMC) simulations compared with experimental results. Balbuena et al. studied the population of 1 0 0 loops in a thin foil using OKMC, with the assumption that all 1 0 0 loops formed from formation and collapse of C15 clusters [22]. They assumed that a certain fraction of small clusters transform into C15, and a fraction of these subsequently collapsed into 1 0 0 loops. Good agreement with experiments was obtained when assuming that the total probability for this 1 0 0 nucleation sequence is 0.1%. The corresponding 1 0 0 nucleation probabilities can be estimated based on our MD simulations. The ABOP predicts 10-18% of clusters formed in cascades to be C15 clusters, and that upon growth 3.33% (1/30) of them collapse into 1 0 0 loops, yielding a nucleation probability in the range 0.33-0.6%. Considering the poor statistics of the growth simulations (1/30), the uncertainty of this probability is large, and is therefore roughly consistent with the 0.1% empirical estimate by Balbuena et al. [22]. In contrast, the corresponding nucleation probability in the M07-B potential is 2.5-5%, indicating an overestimation of 1 0 0 loops compared to experiments.

Accumulation of C15 clusters under continuous irradiation
In the previous sections we have confirmed that C15 clusters can form directly as a result of a collision cascade, and that they are energetically and dynamically stable up to a critical size, after which they collapse into dislocation loops. In this section, we observe the entire life cycle of C15 clusters in continuous irradiation, from formation to collapse or absorption by a larger cluster. This is similar to the recent work by Chartier et al., where electron irradiation was mimicked by inserting Frenkel pairs at regular intervals [3]. Our work deals with overlapping cascades induced by recoils from neutron or ion irradiation at low temperatures where thermal migration is negligible (due to the time scale limitation of MD). We analyse the population of C15 clusters in simulations of 2000 cumulative 5 keV collision cascades. New simulations are carried out with the ABOP, and the results are compared with further analysis of our previous simulation data from Ref. [2] using the EAM potentials. Fig. 5a shows the fraction of clusters that are identified as C15-like as a function of dose in the three different interatomic potentials. The data are averaged over three separate simulation runs. The AM04 potential shows by far the lowest fractions of C15 clusters, as is expected since the stability of C15 clusters compared to dislocation loops in this potential is strongly underestimated (Fig. 3). Both the ABOP and the M07-B potential show a saturation of around 30% C15 clusters out of all interstitial clusters  containing more than two interstitial atoms. The identified cluster fractions are, however, slightly misleading for two reasons. First, dislocation loops are frequently attached to C15 clusters, in which case the entire cluster is here counted as one C15-like cluster. Second, as the dose increases, the mobile 1/2 1 1 1 loops are rapidly growing to large sizes as they absorb small clusters, resulting in a rapid decrease in the total number of clusters. Therefore, in Fig. 5b we show the fraction of interstitial atoms identified as belonging to a C15-like geometry. These fractions show that after reaching a maximum C15 population at around 0.02-0.05 dpa, the fraction of C15-like interstitials starts decreasing until they start saturating at around 10% in the ABOP and M07-B potentials. Both the ABOP and the M07-B potentials show the same trend.
Figs. 5c-d show the mean and maximum sizes of the C15 clusters as a function of dose. The majority of C15 clusters are small at all doses, with the average size between 5 and 10 interstitials in all potentials. The maximum size of any C15 cluster reaches about 30 interstitials in the ABOP and M07-B potentials, and around 10 interstitials in the AM04 potential. This corresponds to diameters in the 1-2 nm range. The maximum sizes are very close to the crossovers in energy between C15 and 1/2 1 1 1 loops (Fig. 3). The maximum sizes are much smaller than observed in the controlled growth simulations, which is due to the continuous irradiation and presence of larger mobile dislocation loops resulting in frequent cluster transformations due to cascade overlap, and absorption of smaller clusters. These effects are, however, severely overestimated due to the rapid dose rate required in MD simulations, and in experimental conditions the C15 clusters would likely grow more undisturbed towards larger sizes. Additionally, since all potentials underestimate the cross-over in energy between C15 clusters and loops, the size range of C15 clusters observed here are likely underestimated accordingly. We also note that overall the ABOP and M07-B potentials show remarkably similar results in Fig. 5, owing to the similar predicted stability of C15 clusters.
The evolution of the C15 cluster statistics shown in Fig. 5 can be understood by visually analysing the defect structure at different doses. Fig. 6 shows snapshots of the interstitial atoms at different doses in the ABOP, with C15-like clusters coloured orange and other interstitials black. C15 clusters are frequently forming already at very low doses (Fig. 6a). The C15 population is largely controlled by the presence and growth of mobile 1/2 1 1 1 loops. At doses below 0.05 dpa (Fig. 6b), the clusters are still relatively small, and in the range where C15 clusters are energetically favoured. As the clusters absorb more interstitials, 1/2 1 1 1 loops eventually start dominating and absorbing other small clusters, including C15 clusters, as they migrate rapidly through the simulation box. This is visible in Fig. 5b as a stable decrease in interstitial fraction and is illustrated in Fig. 6c. At higher doses, the C15 population has reached an equilibrium where new clusters are formed roughly at the same rate as they are absorbed by the large 1/2 1 1 1 loops, corresponding to Fig. 6d.

Summary and conclusions
Using molecular dynamics simulations and a new interatomic bond-order potential that reproduces the correct relative stability of C15 clusters, 1 0 0 loops and 1/2 1 1 1 loops, we investigated the formation, growth, collapse, and population of self-interstitial C15 Laves phase clusters in α-iron. We found that about 5-20% of the interstitial clusters formed in collision cascades are C15-like clusters. Furthermore, we observed that when growing C15 clusters eventually collapse, 1/2 1 1 1 loops is the predominant product at least at low temperatures. During collapse of C15 into 1/2 1 1 1 loops, the loop is often detached and escapes before a full transformation is complete, leaving behind a small stable C15 cluster. Hence, C15 clusters can act as long-lived and stationary emitters of 1/2 1 1 1 loops. Our results also show that C15 clusters can in rare occasions collapse into 1 0 0 loops, in line with earlier studies, and that this is likely more common at higher temperatures. Under rapid continuous irradiation, up to 30% of the self-interstitial clusters or 10-20% of the total number of self-interstitial atoms were identified as C15-like. The majority of C15 clusters remain small in size, in the range of 1-2 nm. We conclude that the frequent formation and large populations of the small and immobile C15 clusters can have a significant impact on the microstructural evolution of iron-based materials.

Acknowledgements
This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 and 2019-2020 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. Grants of computer capacity from CSC -IT Center for Science, Finland, as well as from the Finnish Grid and Cloud Infrastructure (persistent identifier urn:nbn:fi:research-infras-2016072533) are gratefully acknowledged.

Data Availability
The raw data required to reproduce these findings are available upon request from the authors. The processed data required to reproduce these findings are available upon request from the authors.

Potential function
The equations defining the ABOP are listed below. For a more complete description of the parameters involved, see Refs. [1,2].
The ABOP is joined with the the universal repulsive Ziegler-Biersack-Littmark potential [3], V ZBL (r ij ), according to where .  2. Notes on the potential refitting As mentioned in the main article, our starting point when developing the potential was the existing ABOP parametrisation by Müller et al. [4] with the short-range addition by Björkas and Nordlund [5] (which we will call "ABOP-MEABN", and the revised version simply "ABOP"). The goal was to adjust the parameters in order to improve the energetics of single interstitials and vacancies, including migration energies, and to (more importantly) correctly stabilise small C15 clusters relative to dislocation loops in accordance with DFT predictions. To this end, the parameters affecting these properties were iteratively and manually adjusted until a satisfactory agreement with the target data was achieved. The stability of C15 clusters was achieved by changing the optimal interatomic angle defined by the parameter h, which was adjusted to reproduce the cohesive energy of bulk C15 and formation energies of small di-and tetra-SIA C15 clusters in bcc iron. Reproducing the relative stability of C15 clusters and dislocation loops in quantitative agreement with DFT for all sizes was not possible.
Decreasing h further can increase the stability of large C15 clusters towards DFT accuracy, but consequently leads to poor elastic constants of bulk bcc and too-stable smaller C15 clusters. Furthermore, we found that the energies of single and small SIA configurations as well as the vacancy formation energy could be improved by increasing the value of the α parameter, which initially had been set to zero by Müller et al. A non-zero α enables a dependence on the relative bond lengths in the three-body contribution. Müller et al. noted that using a non-zero α resulted in a reduced thermal stability of the fcc phase. However, since we are here mainly interested in radiation damage in bcc iron, sacrificing the stability of the fcc phase of iron is acceptable in favour of significantly improving the energetics of small defect configurations.
The cutoff range of the original parametrisation by Müller, and especially following the short-range modification in [5], is very close to the second-nearest Supplementary material 3 neighbour in the bcc crystal, leading to peculiar temperature effects (see Fig. 2). We therefore extended the cutoff range slightly in order to achieve a reasonable thermal expansion and better elastic properties at finite temperatures. The new cutoff radius was chosen to still be between the second and third nearest neighbour atom in bcc, while also being in-between neighbour shells in e.g. fcc and bulk C15 to avoid unphysical effects. In addition, the cutoff range was carefully chosen to remove spurious energy barriers in the migration paths of vacancies and self-interstitials. Sudden changes in the total coordination of the migrating atom due to atoms entering and exiting the cutoff range were seen to lead to unphysical peaks in the migration barriers.
Finally, to account for the above parameter adjustments, we slightly modified the values of the dimer energy and bond length (D 0 and r 0 ) to reproduce good cohesive energies and lattice constants of bcc iron and other phases. All unmentioned parameters were kept unchanged, including the parameters of the angular function (except the parameter h), and the parameters controlling the stiffness and the dependence on coordination (β, S, and γ). We note that the adjustments we have made eliminates the spontaneous bcc-fcc phase transition, which is correctly reproduced by the original parametrisation.
The short-range part of the original Fe ABOP introduced in [5] was also changed by modifying the parameters of the Fermi function used for the smooth transition from the ZBL potential to the Tersoff-like near-equilibrium part. While the previous Fermi parameters lead to reasonable threshold displacement energies, the many-body interaction curves when moving an atom in a rigid lattice along a given crystal direction (also called "quasi-static drag" or "sudden relaxation" curves) are not in agreement with DFT results. The new parameters of the Fermi function were chosen so that these repulsive many-body interaction curves are accurately reproduced (Fig. 4). Previously, when modifying the short-range part of an EAM potential [6], we used these manybody interaction curves to validate the potential after having obtained reasonable threshold displacement energies. Here, we found that a reverse approach is much more efficient. That is, tuning the short-range potential to reproduce the repulsive many-body energies from DFT automatically resulted in good threshold displacement energies (as long as other relevant properties, such as the formation energy and volume of the lowest-energy self-interstitial configuration, are well reproduced).
The original ABOP overestimates the melting temperature by about 500 K. Our revised version produces an even higher melting point, about 2450 K compared to the experimental value 1811 K. The significant overestimation of the melting point may affect the results of collision cascades by limiting the mixing efficiency during the heat spike. However, lowering the melting point is not possible without affecting more important properties.

Benchmarking the revised interatomic potential
In what follows, we show the results obtained from benchmarking tests of the revised ABOP. Most results are compared with the previous Fe ABOP (ABOP-MEABN) and two EAM potentials (AM04 [7] and M07-B [8,6]).
The formation and binding energies of single and small clusters of self-interstitial atoms (SIAs) and vacancies were calculated by relaxing non-cubic bcc systems containing the defects and around 2 000 atoms. All defect configurations were relaxed to a force convergence of 10 −8 eV/Å at zero pressure.

Supplementary material
4 Threshold displacement energies (TDEs) were simulated using the same methods as described in details in [9,6]. A total of 5 000 random directions were sampled to obtain full angular map of the TDEs.        [14]. Since the potentials were fitted to slightly different experimental data, the elastic constants are normalised by the value at 0 K. The temperatures are also given as fractions of the melting temperature, due to different melting points in the potentials.

135
Distance along given crystal direction (Å) Potential energy difference (eV) Figure 4: Energy difference for step-wise static movement of one atom along different crystal directions in a rigid bcc lattice. DFT data are from Ref. [17]. Exp. DFT-psd M07-B AM04 ABOP Figure 6: Threshold displacement energies. The DFT-psd data are from Ref. [17]. The arrow at the experimental value at [1 1 0] indicate the estimated value > 30 eV [18].  Figure 9: Peierls barrier of screw dislocation movement. Interestingly, despite reproducing the shape and saddle point energy of the [1 1 1] string displacement (Fig. 8) and the stacking fault energy profiles (Fig. 7), the ABOP fails completely in reproducing the corresponding shape and saddle point energy of the Peierls barrier. While all potentials shown here are unable to reproduce the DFT profile, they all, however, predict the correct compact non-degenerate core structure.