Interatomic Fe-Cr potential for modeling kinetics on Fe surfaces

To enable accurate molecular dynamics simulations of iron-chromium alloys with surfaces, we develop, based on density-functional-theory (DFT) calculations, a new interatomic Fe-Cr potential in the Tersoff formalism. Contrary to previous potential models, which have been designed for bulk Fe-Cr, we extend our potential fitting database to include not only conventional bulk properties but also surface-segregation energies of Cr in bcc Fe. In terms of reproducing our DFT results for the bulk properties, the new potential is found to be superior to the previously developed Tersoff potential and competitive with the concentration-dependent and two-band embedded-atom-method potentials. For Cr segregation toward the surface of an Fe-Cr alloy, only the new potential agrees with our DFT calculations in predicting preferential segregation of Cr to the topmost surface layer, instead of the second layer preferred by the other potentials. We expect this rectification to foster future research, e.g., on the mechanisms of corrosion resistance of stainless steels at the atomic level.


I. INTRODUCTION
Iron-chromium alloys are not only scientifically interesting due to their peculiar material properties, but also play an important technological role as the base component for stainless steels [1,2]. From a basic-science perspective, Fe-Cr alloys at varying relative compositions and structures exhibit intricate phenomena such as giant magnetoresistance [3], a spin-ice phase [4], and the "475°C embrittlement" effect [5]. From an application point of view, precision design of steels-beyond the traditional method of empirical trial and error-benefits from an atomic-level understanding of the ins and outs of Fe-Cr alloys. Of particular importance is to explore the mechanisms by which the Cr atoms render the open surfaces of Crcontaining steels resistant to corrosion at and above Cr concentrations of 9-10% [6].
The static properties of the iron-chromium system have been extensively studied at the atomic level. For example, surface [7] and interface [8] energies and surface segregation energies of chromium [9] have been calculated for Fe-Cr alloys using ab initio methods. Importantly, such calculations have shown that the aforementioned critical Cr concentration for the onset of the corrosion resistance in stainless steels closely coincides with an onset of anomalous surface segregation of Cr in Fe-Cr alloys [10,11]. This observation suggests that the segregation of Cr toward the surface plays a crucial role in facilitating the formation of the protective, self-healing layer of iron and chromium oxides that is known to be the reason for the corrosion resistivity of stainless steels [12][13][14][15][16].
The current ab initio modeling methods give reasonably accurate predictions of the static properties of Fe-Cr alloys. The drawback of these methods is that they are computationally heavy, to the extent that the length scales in the modeling can be no larger than a few nanometers and that the costs of studying time-dependent phenomena are still largely prohibitive. In an attempt to overcome these limitations, semi- * pekko.kuopanportti@gmail.com empirical potential models have been developed for the Fe-Cr system. The most recent of these are the concentrationdependent embedded-atom model (CDEAM) [17], the twoband embedded-atom model (2BEAM) [18], and the Tersoff potential [19]. These potentials have been fitted to various basic material properties such as cohesion energies, lattice constants, and elastic properties, and, in general, they describe the point-defect energetics and the solubility of chromium in iron at low concentrations fairly well. As such, the semiempirical Fe-Cr potential models are well suited for modeling bulk Fe-Cr alloys and have been employed, in conjunction with Monte Carlo methods, to study diffusion coefficients and precipitation kinetics [20], vacancy migration near grain boundaries [21], and equilibrium configurations of phase-separated Fe-Cr alloys [22] (using an earlier 2BEAM parametrization [23]).
The aforementioned semi-empirical potential models [17][18][19] turn out to be, however, much less successful in predicting properties of Fe-Cr surfaces [24]. Their disagreement with ab initio calculations is succinctly demonstrated by inspecting the segregation of a Cr atom from the bulk to the surface of a bcc Fe crystal: all three semi-empirical models predict the Cr atom to segregate to the second layer of Fe atoms, whereas according to ab initio calculations the first layer (i.e., the surface itself) should be preferred [24]. In addition, all three models fail to reproduce the ab initio results of Fe-Cr interface energies [8]. As a consequence of these shortcomings, no Monte Carlo simulations of atom kinetics near Fe-Cr surfaces have been performed, and the current atomic-level knowledge of surface physics in Fe-Cr alloys is far from satisfactory, despite the practical importance of the topic in relation to the corrosion of steel surfaces.
In this paper, we develop a new semi-empirical potential model that is suitable for investigating surface physics in Fe-Cr alloys at the atomic level, yet offers performance on par with the existing models in describing the bulk alloy. We choose to work in the Tersoff formalism and use previously developed potentials for the homonuclear Fe-Fe [25,26] and Cr-Cr [19] interactions while constructing a new one for the heteronuclear Fe-Cr interaction. This choice of homonuclear potentials makes the new potential model directly compatible with the previously developed Fe-C [19,27] and Cr-C [19] models [28], so that they can all be combined into an Fe-Cr-C potential model for stainless steel.
The remainder of this paper is organized as follows. In Sec. II, we introduce the Tersoff potential formalism, describe the fitting procedure by which we develop the new Fe-Cr potential, and outline the DFT methods that we use to generate the target data for the fitting procedure. Section III is devoted to presenting the new potential and benchmarking it against the pre-existing potential models and DFT calculations, with both bulk and surface properties considered.
Finally, we summarize our main results and discuss their implications in Sec. IV.

A. Potential formalism
The reactive Tersoff formalism [29][30][31] adopted in this work originates from Pauling's concept of bond order [32,33]; it can also be formally linked [31] to both the tightbinding scheme [34] and the embedded-atom method [35,36]. Since the same formalism has already been described extensively elsewhere [31,[37][38][39], we will give here only a brief overview. In the molecular dynamics code LAMMPS [40], this formalism is available as the potential style tersoff/zbl. Although the only atomic types considered in this work are the elements Fe and Cr, we present the potential formalism below for a general system with an arbitrary number of atomic types.
Let each atom in the system (regardless of its type) be assigned a unique ordinal number, which we will denote using the roman indices i, j, k ∈ N. Furthermore, let (ν i ) i be a finite sequence such that ν i gives the type of the ith atom, with i ∈ {1, . . . , N} and N denoting the total number of atoms in the system [41]. In the Tersoff formalism, the total potential energy E tot of the system can then be written as a sum of individual bond energies: where r i j is the distance between atoms i and j, V ZBL ν i ν j is the universal Ziegler-Biersack-Littmark (ZBL) potential [42] for elements ν i and ν j , V Ter ν i ν j is the pure Tersoff potential for these elements, and is a Fermi function used to join the short-range ZBL and longer-range Tersoff parts smoothly together. The values of the parameters b F and r F are chosen manually such that the potential is essentially the unmodified Tersoff potential at and past equilibrium bonding distances and that a smooth transition to the ZBL potential at short separations is obtained for all realistic coordination numbers.
Incorporation of the ZBL potential [as done in Eq. (1)] is needed to make the potential formalism suitable for modeling nonequilibrium phenomena such as melting or high-energy particle irradiation processes, which typically involve repulsive short-distance interactions originating mainly from the screened Coulomb repulsion between the positively charged nuclei. The ZBL potential is written as where Z ν i is the atomic number of element ν i , with a 0 denoting the Bohr radius, and φ is the universal screening function This screening function has been fitted to the interaction energy between ions, and its accuracy is ∼10% [42]. The Tersoff part V Ter is what chiefly determines the equilibrium properties of the system. It is written as where f c is a cutoff function for the pair interaction, V R is a repulsive and V A an attractive pair potential, and b is a bondorder term that describes three-body interactions and angularity. The pair potentials are of the Morse-like form where D 0 and r 0 are the bond energy and length of the dimer molecule, respectively, and S > 1 is a dimensionless parameter that adjusts the relative strengths of the repulsive and attractive terms. The parameter β is related to the ground-state vibrational frequency ω and the reduced mass µ of the dimer according to The bond-order term is given by where In Eq. (10), θ i jk is the angle between the vectors r i j = r j − r i and r ik , and g is the angular function where γ, c, d, and h are adjustable parameters.
The cutoff function f c appearing in Eqs. (6) and (10) is continuously differentiable and is defined piecewise as where R and D determine, respectively, the center and width of the cutoff interval. Typically, R is chosen to lie midway between the second-and third-nearest neighbors in the relevant equilibrium crystal.

B. Fitting procedure
In order to devise a well-performing Fe-Cr potential in the Tersoff formalism, we use the following approach: The parameters for the Cr-Cr interaction are all taken completely unchanged from Ref. [19], while the parameters for the Fe-Fe interaction are taken unchanged from Refs. [25,26] with the exception that the value of the cutoff distance R Fe Fe is increased to 3.5 Å from the original 3.15 Å to avoid an unphysical increase in Young's modulus of elasticity at elevated temperatures [43]. Furthermore, the two parameters b F, Fe Cr and r F, Fe Cr appearing in the Fermi function [Eq. (2)] of the Fe-Cr potential are chosen to be the same as in Ref. [19]. After making these initial choices, we are left with a total of 17 Fe-Cr potential parameters that we then determine by numerical optimization.
We implement the numerical optimization in MATLAB [44]. To this end, the fitting of the 17 non-predetermined potential parameters is formulated as the nonlinear constrained leastsquares minimization problem where ξ i ∈ [0, 1] is a normalized optimization variable that maps to the closed interval between the minimum and maximum values we allow for the ith potential parameter. The minimum (maximum) allowed values are chosen small (large) enough to have no significant impact on the optimal solution. The target function T is the weighted square sum of the differences between the target values t j and the potential model's predictions p j , where w j ≥ 0 and N data is the number of evaluated quantities in the fitting database. The constrained minimization problem (13) is solved using the trust-region-reflective algorithm [45][46][47] implemented in the MATLAB function lsqnonlin. The potential predictions p j are evaluated by calling the minimize command in LAMMPS [40]. The target values in our fitting database are determined beforehand by DFT calculations, with the database consisting of the following quantities (N data = 40): • FeCr dimer bond length r 0 (target value 2.1428 Å); • formation energies of nine Cr point defects in bcc Fe (see Table III for the specific defects and the target formation energies); • mixing energy of Fe 52 Cr 2 (7 different configurations);  [48]. For a given Cr concentration c Cr N Cr / (N Fe + N Cr ), the SQS cell is a single, computationally constructed 54-atom cell whose lattice sites are occupied by Fe and Cr atoms in such a way that the resulting periodic lattice mimics as closely as possible the first few, physically most relevant radial correlation functions of a perfectly random lattice with the same c Cr . We generate the SQS cells using the mcsqs code [49] of the Alloy Theoretic Automated Toolkit [50].

C. DFT calculations
Before solving the optimization problem (13), we carry out ab initio DFT calculations to determine the target values t j in Eq. (14). All our DFT calculations are performed using the GPAW code [51,52] (version 1.1.0) and the Atomic Simulation Environment (ASE) [53] (version 3.11). Valencecore interactions are modeled with the projector augmentedwave method (GPAW/PAW version 0.8). The generalizedgradient approximation is used in the form of the Perdew-Burke-Ernzerhof exchange-correlation functional [54].
All the bulk DFT calculations are carried out using a 3×3×3 cubic simulation cell of 54 atoms. Wave functions are represented on a real-space grid of 48 × 48 × 48 points, and Brillouin-zone integrations are performed on a Monkhorst-Pack grid [55] of 4 × 4 × 4 k points. All atomic coordinates are relaxed using the Broyden-Fletcher-Goldfarb-Shanno algorithm until the Hellmann-Feynman forces are less than 0.05 eV/Å.
The DFT calculations involving an Fe(100) surface are performed for a 2 × 2 × 5 slab geometry with a distance of 24 Å between the two surfaces. A real-space grid of 48 × 48 × 288 points and a Monkhorst-Pack grid of 4 × 4 × 1 k points are used. The two centermost atomic layers of the slab are fixed to their bulk positions, while all other atoms are relaxed using the FIRE algorithm [56] with a force tolerance of 0.05 eV/Å.

A. New Fe-Cr potential
The numerically optimized parameter values for the new Fe-Cr Tersoff potential are shown in Tables I and II. The corresponding potential input file for the LAMMPS potential style tersoff/zbl is available at [URL].
It is worth noting from Table II that the values of the six three-body parameters α pertaining to the heteronuclear Fe-Cr part have apparently not been numerically optimized in the old Tersoff parametrization of Ref. [19], on the grounds of them all having the same exact value of 1. This should be contrasted with the new parametrization, where all six heteronuclear α parameters have been fully incorporated into the optimization. The resulting increase in the dimensionality of the optimization space may in part explain why we have succeeded in significantly improving the Tersoff potential's agreement with the ab initio target data for both bulk and surface properties, as demonstrated below.
B. Comparison of the potential models in regard to bulk Fe-Cr Let us start with the properties of bulk Fe-Cr alloys and compare the predictions of the new Tersoff potential to those of our DFT calculations, the CDEAM potential, the 2BEAM potential, and the old Tersoff potential. In particular, Figs. 1 and 2 show, respectively, the mixing energy and the Cowley short-range order parameter of the bcc Fe-Cr alloy as a function of the chromium concentration, and Table III lists the formation energies of various isolated Cr defects in bcc iron. In Fig. 1, the alloy mixing energy is determined as the formation energy per atom averaged over different randomalloy configurations at a given chromium concentration c Cr N Cr / (N Fe + N Cr ). The formation energy of a given structure, in turn, is defined in this work as where E tot is the total potential energy of the computational cell, N ν is the number of atoms of element ν in the cell, and E coh (ν) is the cohesive energy of element ν. The cohesive energy E coh (ν) is determined as the total potential energy of a cell containing only atoms of element ν divided by the number Table I. Parameters for the new Tersoff potential of the Fe-Cr system. The Fe-Fe part is the same as originally introduced in Ref. [25] and subsequently augmented with the additional parameters b F and r F [Eq. (2)] in Ref. [26], except for a slightly larger cutoff R that we adopt to avoid an unphysical increase in Young's modulus of elasticity at elevated temperatures; the Cr-Cr part is from Ref. [19]; the Fe-Cr part given in the last column is derived in this work. Note that the three-body parameter α has only one value associated with each of the two homonuclear potentials (α FeFeFe and α CrCrCr , as shown here) but assumes six distinct values for the heteronuclear Fe-Cr part (as listed separately in Table II). All the two-body Fe-Cr parameters are symmetric with respect to interchange of the atomic types (i.e., γ Fe Cr = γ Cr Fe , and similarly for the others).  [57], ab initio calculations predict a negative (positive) mixing energy of Fe-Cr alloys for chromium concentrations below (above) ∼ 6%. As can be seen from Fig. 1(b), our DFT calculations corroborate this result. Although the zerocrossing behavior is qualitatively reproduced by all four potential models under consideration, there are quantitative differences: the new Tersoff potential yields the best fit to the ab initio mixing energies for Cr concentrations < 6%, closely followed by the CDEAM potential, while the 2BEAM and old Tersoff potentials predict the zero crossing to lie at a higher Cr concentration [ Fig. 1(b)]. For chromium-rich alloys, which are included in Fig. 1(a), the new Tersoff potential matches the DFT values significantly better than the other three potentials. For example, at 50% Cr concentration, where the alloy mixing energy has been calculated with the SQS cell [48], the mixing energy given by the new Tersoff potential is only 0.2% larger than the DFT value of 4.663 eV, whereas the CDEAM, 2BEAM and old Tersoff potentials differ, respectively, by 15%, −9.7%, and −65% from the DFT result. The formation energies of selected Cr point defects in bcc Fe are presented in Table III. In general, all four potential models are in fairly good agreement with the target DFT values, with percentage errors typically < 10%. If we take the nine Cr point defects listed in Table III and compute the symmetric mean absolute percentage error (SMAPE) of the predictions of each potential model [58], we obtain the following ranking (from best to worst): new Tersoff (SMAPE of 2.1%), old Tersoff (7.2%), CDEAM (13%), and 2BEAM (17%). On average, the new Tersoff potential is therefore in a better agreement with our ab initio point-defect energies than the other three potential models. We note in particular that the new Tersoff potential provides the best match for the negative DFT value of the Cr substitutional defect.

Interaction
Besides comparing formation energies, we may examine the degree of short-range ordering in a thermally equilibrated Fe-Cr alloy. To this end, the Cowley short-range order parameters α can be defined as [60] α (k) where Z (k) Fe and Z (k) Cr are the average numbers of Fe and Cr atoms in the kth neighbor shell. We further define the linear combination (2) Cr 14 (17) relevant to bcc lattices. If β < 0, a Cr atom prefers to have Fe atoms as its nearest neighbors; if β > 0, Cr prefers Cr neighbors; β = 0 corresponds to a random alloy.  simulations in the canonical ensemble, because this has been found to be more efficient in moving the atoms than the conventional Metropolis algorithm [24]. It should be noted that these Monte Carlo-MD calculations are pure equilibrium simulations with no kinetics involved. Figure 2 shows the order parameter β at 700 K as a function of the Cr concentration for the four potential models, along with three experimental data points determined by diffuseneutron-scattering measurements [59]. While none of the four potential models yields a particularly good fit to the experimental data, the CDEAM and 2BEAM potentials at least qualitatively reproduce the observed change of sign of β from negative to positive around c Cr = 0.1 [ Fig. 2(a)]. For both Tersoff potentials, β remains negative up to high Cr concentrations, eventually turning positive at c Cr ≈ 0.37 for the new Tersoff potential and at c Cr ≈ 0.76 for the old Tersoff potential [ Fig. 2(b)]. We thus conclude that the new Tersoff potential, although largely failing to match the experimental data, still performs noticeably better than the old Tersoff potential in describing the short-range ordering in Fe-Cr alloys.
We have also computed the migration energy barries related to the diffusion of a vacancy-Cr-substitutional pair in bulk bcc Fe using the nudged elastic band (NEB) method [61,62] implemented in LAMMPS. We consider both a process where a Cr subsitutional moves to a vacancy in the nearest-neighbor site and a process where the vacancy migrates from the nearest-neighbor to the second-nearest-neighbor site of the Cr atom (see Fig. 3). The obtained barrier energies for the four potential models are listed in Table IV along with DFT results by Messina, Nastar, Garnier, Domain, and Olsson [63]. By far the best agreement with the DFT values is given by the 2BEAM potential (with the largest relative difference being < 6%), followed by the CDEAM potential, which overestimates the Cr-vacancy migration energy by 53% but is otherwise close to the target values. The new Tersoff potential is in fairly good agreement with DFT for the Cr-vacancy migration energy E  denotes the migration energy barrier for the solute-vacancy jump, and E mig i j is for an iron atom moving from site j to a vacant site i. Table IV lists the values of these energies barriers within the different models under consideration. We now move to scenarios involving surfaces of Fe-Cr alloys and investigate whether we have achieved our goal of improving upon the performance of the existing potential models in predicting properties of such surfaces.
Let us first consider the segregation of Cr atoms to the (100) surface of bcc Fe. The segregation energy of Cr from a given reference region A to another region B is defined as the net change in energy when a Cr atom is transferred from A to B and an Fe atom from B to A, E Cr segr,A→B = E tot (Fe at A, Cr at B) − E tot (Fe at B, Cr at A), (18) where E tot is the total energy of the computational cell. In our case, B will be one of the top surface layers and A will be a bulk site. The energies are calculated using a slab geometry with periodic boundary conditions in y and z directions and the slab thickness L x large enough that further increases in L x cause no changes in the segregation energies. The reference bulk position A is taken from the middle of the slab. For nonzero background Cr concentration, the segregation energy values are averaged over 10 000 random alloy configurations.
The Cr segregation energies to the topmost surface layer (E L1 segr ) and the second topmost layer (E L2 segr ) are shown in Figs. 4(a) and 4(b), respectively, as a function of the Cr concentration c Cr for both the DFT and the four potential models. Due to computational limitations, the GPAW DFT calculations are limited to zero background Cr concentration. The DFT results for c Cr > 0 are obtained with a basis set of exact muffin-tin orbitals (EMTO) [64][65][66] in combination with the coherent potential approximation [67,68], which circumvents the need to average over a large number of alloy configurations as done for the potential models. Note that only the zero-concentration segregation energies of the new Tersoff potential have been explicitly fitted (with the target values being the GPAW ones); the data for c Cr > 0 should be regarded as predictions of the model.
The old Tersoff potential from Ref. [19] is observed to be far away from the ab initio segregation energies for both layers. The CDEAM and 2BEAM potentials give a fairly good fit to the EMTO results for E L1 segr but are far off from both the GPAW and EMTO results for E L2 segr . Importantly, as can be observed from Fig. 4(c), all three pre-existing models (CDEAM, 2BEAM, and old Tersoff) yield E L2 segr − E L1 segr < 0 at all Cr concentrations shown and thus predict segregation of Cr to the second atomic layer instead of the topmost one. This is in stark contrast to the ab initio results, for which E L2 segr − E L1 segr is positive at all investigated concentrations. Fortunately, this shortcoming is fixed by our new Tersoff potential, for which E L2 segr − E L1 segr > 0 in the entire range 0 ≤ c Cr ≤ 0.48 and the fit to the GPAW value of E L2 segr − E L1 segr = 0.3696 eV at c Cr = 0 is excellent, the relative error being 2.5%.
We have also investigated the migration of Cr in the vicinity of an Fe(100) surface. In Table V, we list the migration energy barriers for the process where a surface Cr atom moves to a vacancy at a nearest-neighbor site in the second layer. Although we do not have a DFT value for this migration energy to compare against, it is worthwhile to note that the 2BEAM potential and the new Tersoff potential both yield an energy barrier close to zero whereas the value given by the old Tersoff potential is much higher than the others.  The asterisk is the GPAW ab initio result at c Cr = 0 and the squares are the ab initio results from Ref. [24] obtained using the EMTO method. The lines are guides to the eye.

IV. DISCUSSION
In summary, we have presented a new interatomic Fe-Cr potential that is suitable for molecular dynamics simulations of surface physics in Fe-Cr alloys. The potential was for- mulated in the Tersoff formalism, allowing it to be combined with previously developed Fe-C and Cr-C potentials [19] into a model of the stainless-steel system Fe-Cr-C. The new potential parameters were optimized by fitting to a structural database consisting of ab initio results not only for bulk alloys but also for the segregation of Cr atoms to the (100) surface of an Fe-Cr alloy. We compared the performance of the new potential to that of three pre-existing Fe-Cr potentials with regard to the fitting database as well as to bulk-and surfacerelated quantities not involved in the fitting. For the bulk properties, the new Tersoff potential was found to perform at the same overall level as the pre-existing EAM potentials considered, namely, the CDEAM and the 2BEAM potentials. On the one hand, the new Tersoff potential provided the closest match of all the tested potentials to our DFTcalculated mixing energies and point-defect energies. On the other hand, the 2BEAM potential was the best performer when it came to the short-range ordering in random Fe-Cr alloys (qualitatively reproducing the experimentally observed sign change of the order parameter at c Cr ≈ 0.1) and to the vacancy-Cr diffusion in bcc Fe (yielding energy barriers within 6% of the DFT values). For the tested bulk data, the new Tersoff potential performed significantly better than the old Tersoff potential from Ref. [19], even though only the latter was fitted exclusively to bulk quantities.
The main objective we set at the beginning was for the new potential model to outperform the existing models by correctly predicting the surface-segregation behavior of Cr in Fe-Cr alloys. We achieved this goal in the sense that the new Tersoff potential yields the same ordering of Cr segregation energies as the DFT calculations, namely, E L1 segr < E L2 segr . This is in contrast to the other potential models, for which E L1 segr > E L2 segr and which thus predict Cr to segregate primarily to the second layer instead of the first. The agreement between the new Tersoff potential and the DFT calculations is far from perfect, however. For example, the new Tersoff potential predicts a sign change in the first-layer segregation energy at a significantly smaller Cr concentration (∼ 2.5%) than do the DFTbased EMTO calculations (∼ 8%) [10,24]. In light of the above, we conclude that the new Tersoff potential appears to be the best potential for scenarios where both surface and bulk properties of Fe-Cr alloys are of importance. In bulk systems without boundaries, however, the 2BEAM potential is perhaps the most optimal choice, especially if we take into account its lower computational cost compared to the Tersoff potentials. The old Tersoff potential performs worse than the new one in almost all of our tests, and hence it seems difficult to justify its future use.