Thermomechanical response of NiTi shape-memory nanoprecipitates in TiV alloys

We study the properties of NiTi shape-memory nanoparticles coherently embedded in TiV matrices using three-dimensional atomistic simulations based on the modiﬁed embedded-atom method. To this end, we develop and present a suitable NiTiV potential for our simulations. Employing this potential, we identify the conditions under which the martensitic phase transformation of such a nanoparticle is triggered—speciﬁcally, how these conditions can be tuned by modifying the size of the particle, the composition of the surrounding matrix, or the temperature and strain state of the system. Using these insights, we establish how the transformation temperature of such particles can be inﬂuenced and discuss the practical implications in the context of shape-memory strengthened alloys.


I. INTRODUCTION
Near-stoichiometric NiTi alloys have been known to display shape-memory properties since the 1960s and the alloy itself has since become known as Nitinol [1,2]. Over the decades, it has been employed in a wide array of different applications ranging from stents to actuators. However, it remains impractical as a structural material due to its low stiffness and strength in both austenite and martensite phases [3,4].
An illustration of this mechanical inadequacy is the Nitinol paperclip, which can be deformed to the point where it is simply a piece of bent wire. Without fail, it recovers its original shape upon contact with boiling water, yet remains useless as an actual paperclip as it lacks sufficient strength to hold even a few sheets of paper together.
A potential route to mechanically strengthen NiTi alloys could be using their intrinsic shape-memory property by dispersing coherent austenitic NiTi nanoprecipitates within metallic matrices of appropriate lattice constants. In general, the properties of nanostructured NiTi alloys may differ significantly from those of bulk Nitinol, yet they retain their characteristic phase transformation and the resulting shapememory effect [5,6]. It has been speculated that this could lead to a phase-transformation-based strengthening effect [7]. In particular, hopes are that it is possible to keep nanoprecipitates of supercooled austenite embedded in a stiff matrix below the critical temperature. In this case, the stress field of cracks propagating in the vicinity of such an unstable nanoprecipitate could trigger the martensitic phase transformation, which should lead to a stress relief, stopping crack, or dislocation propagation. Such a material design approach is considered promising to confer self-healing features to metallic alloys [7]. Alternative concepts involve acoustic or magnetic detection of the martensitic transformation of such embedded shapememory particles to facilitate damage detection [8,9].
In this work, we study whether the ternary NiTiV system is a promising candidate material to realize such a scenario. Specifically, we show by means of atomistic molecular-dynamics simulations that a sufficiently stiff TiV matrix can indeed suppress the martensitic phase transformation of such nanoparticles. We report the transformation temperature as a function of matrix composition and study how readily the phase transformation can be triggered by heat treatments and externally applied stress and strain. The results are discussed from both a theoretical and an application-oriented perspective.

II. METHODOLOGY
We use three-dimensional atomistic molecular dynamics (MD) simulations based on the modified embedded-atom method (MEAM) to study the behavior of NiTi nanoprecipitates embedded in TiV matrices (Fig. 1) as a function of temperature as well as external stress and strain.
For this purpose, we have developed an interatomic potential for the ternary NiTiV system based on the 2NN MEAM formalism [10][11][12]. This approach has already been successfully applied for various metals and covalent materials overcoming several drawbacks of the original MEAM model [13]. Because the 2NN MEAM description of a ternary system is based on constituent unary and binary descriptions, potentials for the pure V, Ni, and Ti systems and the binary VNi, VTi, and NiTi systems are required as prerequisites for the construction of the ternary potential. 2NN MEAM potentials for pure Ni, pure Ti, and the binary NiTi system are provided by our previous work [6]. For pure V, a potential was previously reported by Lee et al. [11], but this potential reproduces physical properties at finite temperatures less accurately-in particular the properties related to melting. In the present study, a potential for the pure V system was developed first and potentials for binary VNi and VTi systems were subsequently developed based on the unary potentials. The ternary potential for the NiTiV system was then completed by extending the constituent unary and binary potentials. Details of the fitting procedure are given in Appendix A and a comprehensive benchmark in the Supplemental Material [41]. The potential itself is also made available online [42].
A fitting database of atomic energies and forces was constructed by performing density-functional theory (DFT) calculations for various atomic configurations. In order to ensure transferability of the potential to a variety of applications, configurations resulting from various temperature, strain, and defect conditions were included in the database. Then, the optimization of potential parameters was performed by minimizing errors between the DFT database and corresponding potential results. The required DFT calculations have been performed using the VASP software [14,15]. The computational details of these calculations can be found in Appendix B.
Cubic supercells with edge lengths l = 225Å, 273Å, and 290Å containing 70 × 70 × 70, 85 × 85 × 85, and 90 × 90 × 90 body-centered cubic (bcc) unit cells, respectively, were constructed. A spherical region with a diameter of d = 10-17 nm was cut from the center of each cell and replaced with a B2based NiTi nanoprecipitate such that the interface remained fully coherent. Various amounts of Ti atoms in the remaining matrix region were randomly replaced by V atoms, resulting in cells with matrix V content of 27.5, 35, 42.5, and 50 %. The V content x V is always given as the molar fraction of the matrix, not with respect to the entire cell. These concentrations have been chosen such that the two-phase region in the NiTiV phase diagram as suggested in experimental literature [16,17] should be adequately sampled. We believe that this two-phase region offers an excellent chance for practical realization of nanodispersed NiTi precipitates in a stiff V-rich matrix.
After construction, periodic boundary conditions were applied to these cells in all three directions. Equilibration was performed at T = 550 K in terms of both internal coordinates and cell shape degrees of freedom. Then, heating and cooling were simulated in an NpT ensemble with predefined temperature schedules changing the temperature by δT = 15 K per 10 000 steps starting at T = 485 K down to 170 K, then back up again (cyclic heat treatment). The results of these simulations are reported and discussed in Sec. III.
Tensile tests were simulated in the same supercells by applying strain along one of the three crystallographically equivalent 100 axes (from now on referred to as strain axis). This has been done at several fixed temperatures in the range 260-300 K. Generally the strain rate was chosen as 5 × 10 7 s −1 , but other strain rates were also tested to investigate the impact of the strain rate. The respective specimens were strained under tensile load until 3.5% tensile strain, then successively unloaded at a strain rate with the same absolute value and opposite sign. Corresponding results are discussed in Sec. IV.
All MEAM simulations were carried out in LAMMPS [18] employing a time step of δt = 2 fs.

III. COMPOSITE UNDER HEAT TREATMENT
We want to mimic an experimental situation where the nanoparticles are well separated. To ensure that our simulation cells are sufficiently large for that purpose and to achieve convergence with respect to the edge length l (see geometry in Fig. 1), careful tests have been performed. Too small values of l could lead to (short-range) interaction between two periodic images or to a situation where the matrix is no longer large enough to hold its lattice constant reliably against the eigenstrain of the embedded precipitate.
To this end, cubic supercells containing 70 × 70 × 70, 85 × 85 × 85, and 90 × 90 × 90 (cubic) bcc unit cells and matrix V contents of 27.5%, 35%, 42.5%, and 50% have been constructed with a spherical nanoprecipitate according to the approach outlined in Sec. II. These cells have been equilibrated, then subjected to heat treatment. Every 10 000 MD steps (corresponding to 20 ps of simulation time), a snapshot of the real-space cell was analyzed and the phase fractions of B19 -type martensite and B2-type austenite were obtained using an adaptive cutoff common-neighbor analysis (AC-CNA) as implemented in the OVITO software [19]. These phase fractions were then studied as a function of temperature T for the various precipitate and cell sizes.
Results for the largest particle diameter d = 17 nm and the three different cell sizes l are shown in Fig. 2. We define the critical transformation temperature M s to be the point of largest slope, which is easier to specify accurately than the actual onset of martensite formation. It is evident from Fig. 2 that the resulting transformation temperature is converged to within ±10 K and the saturation fraction to within a few % already for the supercell with 70 3 bcc unit cells, corresponding to an edge length of l ≈ 225Å. Similar convergence tests have been performed for the other matrix V contents with the same finding (not shown here). In the following, all presented results have been obtained from cells containing 70 3 (cubic) unit cells.
Some selected snapshots of 100 cuts through these MD cells are shown in Fig. 3. The top row shows 100 cuts for a matrix V content of x V = 35%, the middle row for x V = 42.5%, and the bottom row for x V = 50%. Columns pertain to different temperatures, above the transformation temperature (left), at the start of the transformation (middle), and below (right). Atoms have been color coded according to the result of the AC-CNA: Atoms embedded in a mostly bcc-like environment are depicted in blue, fcc-like in green, B19 -like in red. Atoms not readily associated with either of these three are shown in gray.
For the 35% V cell at 390 K [ Fig. 3(a)], most of the cell is blue (bcc), which suggests that the matrix is stable in its β state and the nanoprecipitate is purely austenitic. Some small anomalies accumulate near the interface between nanoprecipitate and matrix. These anomalies provide strain relief, as the lattice constant of Ti 65 V 35 Particle-size dependence of the martensitic transformation temperature. Triggering the transformation requires significantly lower temperatures at small particle sizes. Comparison to nanocrystalline NiTi is based on a computational study employing MEAM [20].
according to MEAM at T = 485 K). This results in a distortion of the lattice, which the AC-CNA attributes partly to fcc (green) and partly to an unknown crystal structure (shown in gray). The anomalies become smaller as the V matrix content is increased [Figs. 3(d) and 3(g)], which is likely a combination of the interfacial lattice mismatch decreasing and the matrix simultaneously becoming stiffer with increasing V content.
The snapshots in the middle column [Figs. 3(b), 3(e) 3(h)] were taken exactly at the onset of martensite formation, which we determined to be at 260 K for the lower V contents and 250 K for 50% V. These figures clearly show that martensite nucleates in the nanoparticle's interior. This nucleation behavior is similar to what has been reported for nanocrystalline NiTi [20]. Finally, Figs. 3(c), 3(f) and 3(i) show the low-temperature state at 180 K, where the nanoprecipitates are mostly transformed into the martensite. As can be seen, several domains of martensite form, which are separated by domain boundaries (gray areas). Again, a similar behavior has been observed previously for nanocrystalline NiTi [20].
Studying systematically a number of cells with various matrix V content x V and precipitate sizes d yields information on the influence of precipitate size and composition on the martensitic transformation temperature M s . The results are shown in Figs. 4 and 5 in two different representations. We can infer from Fig. 4 that for particle sizes below 16 nm the martensitic temperature is dependent on the particle size. In this region, the variation of M s with d is strong, making it possible to reduce M s by more than 50 K when changing the diameter from 16 nm to 12 nm.
For particles larger than ≈ 16 nm, the transformation temperature M s almost corresponds to that of B2, which has a martensitic start temperature of M s ≈ 270 K (MEAM value [6], experiment [21]: M s = 339 K). Figure 4 suggests that M s does not quite approach the MEAM bulk value of 270 K even for large d, which may have one of two causes. First, due to the coherent interface, even the largest particle sizes do not correspond to bulk B2 due to the lattice mismatch. Second, our definition of the transformation temperature (largest slope in phase fraction) is slightly different from the one used by Ko et al. [6] (onset of change in volume), which leads to a shift of about one temperature step (cf. Fig. 2), corresponding to about 10 K.
The above described qualitative trend of the martensitic temperature, i.e., decreasing with decreasing particle size, appears to be generic and has indeed also been reported for nanocrystalline NiTi [20] (black dashed line in Fig. 4).
From the perspective of realizing applications, which require austenite to be present at operating temperatures, Figs. 4 and 5 suggest that this particular choice of matrix is not ideally suited to suppress the martensitic transformation at room temperature. However, before discarding this specific choice of material, it has to be noted that this may be a result of the MEAM potential intrinsically underestimating the experimental transformation temperature M s by 60-70 K (even for bulk Nitinol [6]). Hence, it is more advisable to plot the change δT of the transformation temperature M s instead. δT quantifies how M s changes with respect to bulk Nitinol as a function of particle size d and matrix V content x V . Figure 6 shows a contour map of this shift. The contour map was obtained from a bilinear fit δT = a 0 + a 1 x V + a 2 d + a 3 x V d, with fitting coefficients a i , i = 0-3. The equipotential lines indicate lines of constant temperature shift δT with respect to M s in perfect bulk Nitinol, which is 270 K (MEAM value). Figure 6 visualizes how the transformation temperature can be tuned by the vanadium content of the surrounding bulk material and the size of the nanoprecipitates. In particular, smaller particles and higher x V lead to a lower transformation temperature M s .
For the purpose of designing an alloy, this is an important finding. In order to benefit from a phase-transformation-based strengthening effect, it is crucial that the alloy contains austenitic nanoparticles in an unstable, supercooled state, where the transformation can readily be triggered by the stress field of a moving dislocation or a propagating crack. Figure 6 033610-4 provides a guideline where such a supercooled state can be found, depending on the expected operating temperature of the material.
In fact, Fig. 6 suggests that sufficiently large vanadium concentrations x V can suppress the martensitic transformation to values significantly lower than the martensite formation temperature. This can be rationalized by recalling that there is a substantial change in lattice constants associated with the martensitic transformation, which changes the crystal structure from the B2 austenite at high temperatures to the low-temperature B19 structure [22,23].
The B2 austenite obeys cubic symmetry constraints. According to results of Bührer et al. [23] (which we found to be in excellent agreement with both our DFT and MEAM calculations), it exhibits NiTi nearest-neighbor distances of 2.61Å and second-nearest-neighbor distances of 3.02Å for Ni-Ni or Ti-Ti. The B19 martensite on the other hand features nonequivalent NiTi nearest-neighbor pairs ranging from 2.47Å to 3.33Å and second-nearest-neighbor distances from 2.58 to 3.21Å.
Therefore, some crystallographic directions will be affected by compressive strains, others by tensile strains. The elastic work required to deform the matrix to accommodate the particle's expansion and compression suppresses the martensitic transformation. Changing the bulk V content x V changes the stiffness of the matrix. Figure 7 shows the 0 K elastic constants of a disordered 6750 atom TiV cell (without nanoparticle) as computed from the MEAM potential. Adding more and more V to the material increases the stiffness, which in turn leads to more resistance against deformation associated with a nanoparticle's transformation. As mentioned already, the results indicate that it may be possible to manufacture a situation where such nanoparticles are kept in their austenitic state at operating temperatures of a structural material.

IV. COMPOSITE UNDER TENSILE STRESS
Since our goal is to understand how the nanoparticles impact the mechanical response of the nanocomposite Elastic stiffness coefficients C ij of TiV. The coefficients shown correspond to the C 11 (red), C 12 (blue), and C 44 (pink) stiffnesses. A general trend of increasing C ij stiffnesses with increasing V content is observed. C 12 is fairly constant over a large concentration region, which is related to the elastic instability of pure bcc titanium (for details, see main text). material, we study how the composite and its nanoparticles behave under tensile load. Therefore, the same supercells as discussed in the previous section were subjected to simulated tensile experiments. Such simulations have been performed systematically at T ∈ {260 K, 280 K, 300 K} with matrix V contents of x V ∈ {27.5%, 35%, 42.5%, 50%}, using straincontrolled tensile loading along the y direction and allowing for relaxation along x and z directions.
Most experimental techniques investigate the stress-strain response at strain rates of ≈ 10 −3 s −1 , though some setups for high-strain rate experiments exist [24]. Due to the femtosecond dynamics of atomistic molecular dynamics, strain has to be applied at unphysically high strain rates, several orders of magnitude higher than in experiment. We have analyzed the strain dependence of our results in a practically accessible regime of 10 7 s −1 to 10 8 s −1 . We find that neither the transformation temperature nor the saturation phase fractions at high strains change significantly in this regime. We therefore applied strain in the following simulations at a constant (engineering) strain rate of 5 × 10 7 s −1 . Figure 8 shows the resulting microstructures of such a strain-controlled tensile loading cycle, performed on specimens with x V = 35% and at 280 K. The atoms have again been color coded according to the results of an AC-CNA. Figure 8(a) corresponds to the unstrained nanoparticle, before applying any strain. At a tensile strain of ε ≈ 1.8% (not shown), the nanoprecipitate starts to partially transform into its martensitic state. Figure 8(b) shows the nanoparticle at its peak load ε = 3.5%, after the transformation has occurred. Not all of the austenite has transformed and a saturation effect with respect to the amount of martensite within the precipitate is observed.
After reaching a strain of 3.5%, the simulation cell was unloaded with a strain rate of the same absolute value, but opposite sign. Figure 8(c) shows the fully unloaded nanoparticle, which retains a significant portion of ≈30% martensite. From similar simulations for the different matrix V contents x V specified above, Fig. 9 can be constructed. In these plots, the microstructures shown in Fig. 8 correspond to the 35% panel. The diagram suggests that increasing the vanadium concentration leads to an increase in the critical strain ε c required to initiate the transformation. Thus, very high matrix V contents (x V = 50%) require nearly 3% of external tensile strain to trigger the martensitic phase transformation.
Similarly to the decrease of M s with increasing x V discussed in the previous chapter, we also attribute this shift of ε c to a stiffening of the matrix with increasing V content (cf. Fig. 7). Since the formation of martensite is associated with a change in lattice parameters, a stiff matrix leads to a more effective suppression of the transformation. The willingness of these particles to trigger can thus be tuned via the matrix composition, which is a major finding on the way to utilize the shapememory effect for a controlled active strain-relief mechanism.
The same argument also explains why the phase fraction is very much a function of the bulk stiffness and hence the matrix V content x V . Also, Fig. 9 suggests that for very high x V > 50% the transformation is indeed reverted by the eigenstresses of the matrix acting on the precipitate itself. Figure 10 shows the phase fractions of martensite and austenite within the nanoprecipitate as a function of applied strain for fixed particle size d = 17 nm and fixed x V = 35% for several temperatures. Figure 10(a) illustrates that at high temperatures precipitates do not transform even under significant tensile strain. Note that the two phase fractions occasionally do not add up to unity. This is a result of the AC-CNA sometimes not being able to identify an atom to be located in an explicitly martensitic or austenitic environment. These are visible as occasional gray and green atoms in Fig. 8, but are few in number and thus do not affect the analysis above. Figures 8-10 all suggest that the martensitic transformation is not fully reversed upon unloading the sample for matrix V contents below 50%. A significant amount of martensite remains present in the particles even after fully unloading it. In a final step, we simulated the effect of temperature on such a previously strained specimen to study whether austenite can be restored by heat treatment. To this end, we have heated the simulation cell resulting from a strain-controlled loading cycle as discussed in the previous section from T = 280 K to 485 K. The heating rate was chosen to be 1 K/ps. After each temperature step, a complete AC-CNA was performed.

V. HEAT TREATMENT OF STRAINED NANOPARTICLES
The temperature-dependent martensite content within such a nanoparticle is shown in Fig. 11. As can be seen, most of the martensite can be transformed back to its austenite B2 state by applying such a heat treatment. Figure 8(d) substantiates that the resulting microstructure of the nanoparticle at 385 K after heat treatment corresponds to the original austenite state (colored blue according to the AC-CNA).

VI. SUMMARY
We have presented a ternary 2NN MEAM potential for the NiTiV system and applied it to various atomistic molecular dynamics simulations of austenitic B2 nanoprecipitates embedded in a TiV matrix. Our findings allow three major conclusions.
First, these nanoprecipitates can be kept in a supercooled austenitic state well below the critical martensitic temperature of bulk Ni 50 Ti 50 .
Second, the required temperatures and stresses to trigger the martensitic transformation can be tuned using the composition of the matrix. The derived phase diagram shown in Fig. 6 provides a practical route to choose a desired transformation temperature by selecting a suitable nanoparticle size and matrix composition.
Third, while the transformation is not reversible by unloading the specimen except at high temperatures (T ≈ 300 K) or high matrix V content (x V ≈ 50), in all studied cases it was possible to restore the austenite by heat-treating the alloy after full unloading.
This suggests that further investigations of the NiTiV system are promising. Related systems, such as the NiTiNb alloy, which also offers a similar two-phase region [25] between a B2 phase and a solid solution, also seem promising candidates. These alloys might offer a window of opportunity for a microstructure with a custom-tailored martensitic phase transformation that triggers at just the right stresses and temperatures  11. Transformation of retained martensite into austenite during heat treatment. The phase fraction is given relative to the volume of the nanoparticle. Heating the nanoparticle reverts the martensitic transformation almost completely even at temperatures as low as ≈310 K. The microstructure after such a heat treatment is shown in Fig. 8(d).
to, e.g., delay crack growth or pin dislocation movement at the expected operating temperatures. To the best of the authors' knowledge, similar NiTiV-based compounds have so far only been studied for their hydrogen-related properties [26] at much larger precipitate sizes. Experimental work on synthesizing the compound with nanodispersed precipitates such as discussed in the present paper is ongoing, but promising preliminary results have been recently shown at an international conference on self-healing in Friedrichshafen, Germany [27].

ACKNOWLEDGMENTS
Funding by the Deutsche Forschungsgemeinschaft (SPP 1568) and by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (Grant Agreement No. 639211) is gratefully acknowledged.

APPENDIX A: OPTIMIZATION OF THE NiTiV 2NN MEAM POTENTIAL
In the 2NN MEAM formalism, 14 independent parameters are required to describe a pure element. Four of these parameters, the cohesive energy E c , the equilibrium nearest-neighbor distance r e , the bulk modulus B of the reference structure, and the adjustable parameter d are related to the universal equation of state. Seven further parameters, the decay lengths β 0 , β 1 , β 2 , β 3 and three weighting factors t 1 , t 2 , and t 3 are related to the electron density. The parameter A defines the embedding function, and the parameters C min and C max are responsible for the many-body screening. To describe a binary alloy, 13 independent parameters are necessary in addition to the 14 unary parameters for each species. Again, four parameters are related to the universal equation of state, but now the cohesive energy is replaced by the enthalpy of formation of the reference structure E f . Eight further coefficients (four independent C min and four independent C max ) are required to describe the many-body screening. The last parameter is an electron density scaling factor ρ 0 . Finally, to generalize the formalism to a ternary alloy, six parameters are required in addition to the constituent unary and binary parameters. These are all related to the many-body screening (three C min and three C max ). More details on the meaning of the 2NN MEAM parameters are given in the literature [10][11][12]28,29].
For pure Ni and pure Ti, we used the optimized unary parameters from our previous work [6], as reproduced in Table I. The optimization of the pure V potential was accomplished by partitioning DFT data into two groups. The first group of DFT results comprised structural energies and forces between atoms in perfect and defected (with vacancies and with thermal vibrations from ab initio MD) configurations of bcc, face-centered cubic (fcc), and hexagonal-close-packed (hcp) phases that were used in the parameter optimization. The second group of DFT results comprised phonons and thermodynamic properties which were used for the evaluation of the developed potential (see Supplemental Material [41]). The optimization was conducted by minimizing errors between energies and forces predicted by the DFT calculations and those predicted by temporary potential parameters. The bcc structure was selected as a reference structure and the experimentally reported [30] cohesive energy of bcc V 5.30 eV was selected as the E c value. The resulting unary parameters for the V potential are given in Table I. For NiTi, we used the optimized binary parameters from our previous work [6], as reproduced in Table II. In that work [6], we found that the alloy parameters introduced in the 2NN MEAM model do not provide enough flexibility for the force-matching process in the binary NiTi system, because the pair interaction description is determined not from a specific functional form but from the universal equation of state. The solution to this problem was a pragmatic approach that focused on optimizing only a particular set of properties related to the shape-memory transformation. In the present work, we followed a similar pragmatic approach when fitting potentials for NiV, TiV, and NiTiV. In particular, the optimization was performed using the energies (but not forces) of atomic configurations obtained by DFT calculations. TABLE I. 2NN MEAM potential parameter sets for the pure Ni, Ti, and V systems. The following properties are dimensionful: the cohesive energy E c (eV/atom), the equilibrium nearest-neighbor distance r e (Å), and the bulk modulus B (GPa). See Appendix A for the definition of the other parameters. The reference structures are bcc Ti, fcc Ni, and bcc V. The parameters for Ni and Ti are reproduced from Ref. [6]. Further, the optimization was performed by directly fitting to the following properties related to configurations of solid solutions: the dilute heat of the solution, the vacancy-solute binding energy, the solute-solute binding energy, and the solute migration energy. These properties were calculated by considering atomic relaxations for each configuration. Such a fitting procedure implies that the resulting binary and ternary potentials are mainly intended to properly describe the properties related to solid solutions. We especially focus on these properties because the present application of the developed potential is based on the analysis of the interaction between a solid solution phase and a NiTi shape-memory alloy particle. The optimized binary and ternary parameters are compiled in Table II and Table III, respectively. We confirmed that a radial cutoff distance of 5.0Å is sufficiently large to reproduce various physical properties as well as the martensitic phase transformation of NiTi shapememory alloys with a sufficient computational efficiency. All simulations in the present study were performed using this radial cutoff distance.
For the computational details of the DFT calculations used to fit the potential, see Appendix B. Extensive transferability and accuracy tests of the described potential are available in the Supplemental Material [41]. The potential itself is available online [42].

APPENDIX B: COMPUTATIONAL DETAILS OF THE FIRST-PRINCIPLES CALCULATIONS
The DFT calculations were conducted using the VASP code [14,15] and the projector-augmented wave method [31,32] within the Perdew-Burke-Ernzerhof generalized-gradient approximation [33] for the exchange-correlation functional. A cutoff energy of 350 eV for the plane-wave basis set and the Methfessel-Paxton smearing method [34] with a width of 0.1 eV were used. A k-point mesh of 19 × 19 × 19 k points was selected for the bcc primitive unit cell, and the corresponding k-point density was employed for other structures and supercells. Magnetism was included by considering spinpolarized calculations for Ni-rich supercells and intermetallic compounds. The lattice constants and bulk modulus were calculated by employing the Birch-Murnaghan equation of state [35,36]. For the calculation of defect properties, atomic positions were relaxed at a constant volume and cell shape with the convergence criteria for energy and forces set to 10 −6 eV and 10 −2 eV/Å, respectively.
The vacancy migration energy with a suitable saddle-point configuration was calculated using the nudged elastic band method [37,38]. For the calculation of the surface energy, rectangular cells with 20Å to 25Å thick slabs and a vacuum region of 10Å were employed. Phonon calculations were performed based on the direct force constant approach using the Phonopy code [39]. A body-centered cubic supercell of 250 atoms was used for the phonon calculation with the convergence criteria for energy and forces set to 10 −8 eV and 10 −4 eV/Å, respectively.
Atomic configurations at finite temperatures were obtained by conducting two-step DFT calculations following the concept of the UP-TILD method [40]. First, ab initio MD simulations were performed based on relatively low DFT convergence parameters for a total of 1000 steps with a time step of 1.5 fs. Then, well-converged calculations with a higher cutoff energy and denser k-point mesh were performed using snapshots from ab initio MD simulations to determine accurate energies and forces of each configuration at finite temperatures. 033610-9