# Collinear-spin machine learned interatomic potential for Fe7Cr2Ni alloy

Lakshmi Shenoy, Christopher D. Woodgate, Julie B. Staunton, Albert P. Bartók, Charlotte Becquart, Christophe Domain, James R. Kermode

► To cite this version:

Lakshmi Shenoy, Christopher D. Woodgate, Julie B. Staunton, Albert P. Bartók, Charlotte Becquart, et al., Collinear-spin machine learned interatomic potential for Fe7Cr2Ni alloy. Physical Review Materials, 2024, Physical Review Materials, 8 (3), 10.1103/physrevmaterials.8.033804. [arXiv:2403.07286]

HAL Id: hal-04522806

https://hal.univ-lille.fr/hal-04522806

Submitted on 27 Mar 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Collinear-spin machine learned interatomic potential for Fe7Cr2Ni alloy

Lakshmi Shenoy, Christopher D. Woodgate, Julie B. Staunton, Albert P. Bartók, Charlotte Becquart, Christophe Domain, James R. Kermode

► To cite this version:

Lakshmi Shenoy, Christopher D. Woodgate, Julie B. Staunton, Albert P. Bartók, Charlotte Becquart, et al., Collinear-spin machine learned interatomic potential for Fe7Cr2Ni alloy. Physical Review Materials, 2024, Physical Review Materials, 8 (3), 10.1103/physrevmaterials.8.033804. [arXiv:2403.07286] [hal-04522806]

HAL Id: hal-04522806

https://lilloa.hal.science/hal-04522806v1

Submitted on 27 Mar 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Collinear-spin machine learned interatomic potential for Fe<sub>7</sub>Cr<sub>2</sub>Ni alloy

Lakshmi Shenoy<sup>1,1,*</sup> Christopher D. Woodgate<sup>2</sup> Julie B. Staunton<sup>2</sup> Albert P. Bartók<sup>1,2</sup> Charlotte S. Becquart<sup>3</sup>
Christophe Domain<sup>4</sup> and James R. Kermode<sup>1</sup>

<sup>1</sup>Warwick Centre for Predictive Modelling, School of Engineering, University of Warwick, Coventry CV4 7AL, United Kingdom

<sup>2</sup>Department of Physics, University of Warwick, Coventry CV4 7AL, United Kingdom

<sup>3</sup>Univ. Lille, CNRS, INRAE, Centrale Lille, UMR 8207 - UMET - Unité Matériaux et Transformations, F-59000 Lille, France

<sup>4</sup>Electricite de France, EDF Recherche et Developpement, Departement Materiaux et Mecanique des Composants,
Les Renardieres, F-77250 Moret sur Loing, France

<figure><img src="image_2.png" /></figure>

(Received 4 October 2023; accepted 27 February 2024; published 22 March 2024)

We have developed a machine learned interatomic potential for the prototypical austenitic steel Fe<sub>7</sub>Cr<sub>2</sub>Ni, using the Gaussian approximation potential (GAP) framework. This GAP can model the alloy's properties with close to density functional theory (DFT) accuracy, while at the same time allowing us to access larger length and time scales than expensive first-principles methods. We also extended the GAP input descriptors to approximate the effects of collinear spins (spin GAP), and demonstrate how this extended model successfully predicts structural distortions due to antiferromagnetic and paramagnetic spin states. We demonstrate the application of the spin GAP model for bulk properties and vacancies and validate against DFT. These results are a step towards modeling the atomistic origins of ageing in austenitic steels with higher accuracy.

DOI: 10.1103/PhysRevMaterials.8.033804

I. INTRODUCTION

Austenitic stainless steels are key structural materials with
high mechanical strength and corrosion resistance. They have
many applications, from everyday tools such as household
appliances and cars, to structural components of buildings,
bridges, and industrial machines [1]. Their desirable proper-
ties come from alloying Fe with a concentration of 16–20%
Cr for rust resistance, and about 10% Ni to maintain a
crack-resistant face-centred-cubic crystal structure at all tem-
peratures where it remains solid (< 1800K) [2]. Specific
properties of the steel can be fine tuned by adding small
concentrations (< 1%) of other solutes such as Mo, Mn, P, C
etc., depending on the application of interest [3]. However, the
core properties come from the Fe-Cr-Ni mix, so in this paper
we focus on the prototypical austenitic steel with composition
Fe<sub>7</sub>Cr<sub>2</sub>Ni.

A key application of austenitic steels is in nuclear power
plants. The austenitic steel grades 304 and 316, which have
a composition very close to the model alloy Fe<sub>7</sub>Cr<sub>2</sub>Ni, are
used in the primary containment barriers of light-water nu-
clear reactor pressure vessels (RPV) [4]. The RPV interiors
are exposed to constant neutron bombardment and a wide
range of thermal fluctuations from cryogenic temperatures
during maintenance to operating temperature of ~600 K. Over

time, these harsh conditions trigger the formation of radiation-
induced point defects in the RPV interiors, which eventually
cluster to form larger defects such as voids and dislocation
loops. Macroscopically, these manifest as creep, swelling,
stress corrosion, and other ageing effects, which are safety
concerns [5]. To improve RPV designs and facilitate mitiga-
tion efforts, there is interest in understanding the mechanisms
behind these ageing phenomena.

Experimental data for aging phenomena are limited due
to the challenges in conducting controlled experiments in
extreme conditions, and also challenges in measuring atomic
quantities experimentally [6–8]. We therefore require compu-
tational methods to simulate these phenomena across length
scales, to fill in the gaps between the available experimental
data, and to explain the observed aging effects.

There are various tiers of first-principles methods available for materials modeling, with higher accuracy being accompanied by (prohibitively) higher computational costs. While more complex phenomena such as nuclear cascades or spectroscopic excitations require higher levels of theory, the first principles method of choice for modeling mechanical properties as in our case is density functional theory (DFT), in which the electronic structure of materials is calculated from approximations of quantum mechanics [9]. There have been efforts to study the Fe<sub>7</sub>Cr<sub>2</sub>Ni alloy using DFT [5,10–12]. However, as the computational cost of DFT calculations scales cubically with system size, these studies were restricted to supercells containing up to 256 atoms. These supercells are too small to model extended defects such as voids or dislocations. Also, as DFT calculations are computationally expensive, these studies were restricted to limited statistics, which is not ideal for modeling the variations due to chemical complexity of concentrated alloys.

<footnote>*lakshmi.shenoy@warwick.ac.uk</footnote>
Collinear-spin machine learned interatomic potential for Fe<sub>7</sub>Cr<sub>2</sub>Ni alloy

Lakshmi Shenoy<sup>1,1,*</sup> Christopher D. Woodgate<sup>2</sup> Julie B. Staunton<sup>2</sup> Albert P. Bartók<sup>1,2</sup> Charlotte S. Becquart<sup>3</sup>
Christophe Domain<sup>4</sup> and James R. Kermode<sup>1</sup>

<sup>1</sup>Warwick Centre for Predictive Modelling, School of Engineering, University of Warwick, Coventry CV4 7AL, United Kingdom

<sup>2</sup>Department of Physics, University of Warwick, Coventry CV4 7AL, United Kingdom

<sup>3</sup>Univ. Lille, CNRS, INRAE, Centrale Lille, UMR 8207 - UMET - Unité Matériaux et Transformations, F-59000 Lille, France

<sup>4</sup>Electricite de France, EDF Recherche et Developpement, Departement Materiaux et Mecanique des Composants,
Les Renardieres, F-77250 Moret sur Loing, France

<figure><img src="image_2.png" /></figure>

(Received 4 October 2023; accepted 27 February 2024; published 22 March 2024)

We have developed a machine learned interatomic potential for the prototypical austenitic steel Fe<sub>7</sub>Cr<sub>2</sub>Ni, using the Gaussian approximation potential (GAP) framework. This GAP can model the alloy's properties with close to density functional theory (DFT) accuracy, while at the same time allowing us to access larger length and time scales than expensive first-principles methods. We also extended the GAP input descriptors to approximate the effects of collinear spins (spin GAP), and demonstrate how this extended model successfully predicts structural distortions due to antiferromagnetic and paramagnetic spin states. We demonstrate the application of the spin GAP model for bulk properties and vacancies and validate against DFT. These results are a step towards modeling the atomistic origins of ageing in austenitic steels with higher accuracy.

DOI: 10.1103/PhysRevMaterials.8.033804

I. INTRODUCTION

Austenitic stainless steels are key structural materials with
high mechanical strength and corrosion resistance. They have
many applications, from everyday tools such as household
appliances and cars, to structural components of buildings,
bridges, and industrial machines [1]. Their desirable proper-
ties come from alloying Fe with a concentration of 16–20%
Cr for rust resistance, and about 10% Ni to maintain a
crack-resistant face-centred-cubic crystal structure at all tem-
peratures where it remains solid (< 1800K) [2]. Specific
properties of the steel can be fine tuned by adding small
concentrations (< 1%) of other solutes such as Mo, Mn, P, C
etc., depending on the application of interest [3]. However, the
core properties come from the Fe-Cr-Ni mix, so in this paper
we focus on the prototypical austenitic steel with composition
Fe<sub>7</sub>Cr<sub>2</sub>Ni.

A key application of austenitic steels is in nuclear power
plants. The austenitic steel grades 304 and 316, which have
a composition very close to the model alloy Fe<sub>7</sub>Cr<sub>2</sub>Ni, are
used in the primary containment barriers of light-water nu-
clear reactor pressure vessels (RPV) [4]. The RPV interiors
are exposed to constant neutron bombardment and a wide
range of thermal fluctuations from cryogenic temperatures
during maintenance to operating temperature of ~600 K. Over

time, these harsh conditions trigger the formation of radiation-
induced point defects in the RPV interiors, which eventually
cluster to form larger defects such as voids and dislocation
loops. Macroscopically, these manifest as creep, swelling,
stress corrosion, and other ageing effects, which are safety
concerns [5]. To improve RPV designs and facilitate mitiga-
tion efforts, there is interest in understanding the mechanisms
behind these ageing phenomena.

Experimental data for aging phenomena are limited due
to the challenges in conducting controlled experiments in
extreme conditions, and also challenges in measuring atomic
quantities experimentally [6–8]. We therefore require compu-
tational methods to simulate these phenomena across length
scales, to fill in the gaps between the available experimental
data, and to explain the observed aging effects.

There are various tiers of first-principles methods available for materials modeling, with higher accuracy being accompanied by (prohibitively) higher computational costs. While more complex phenomena such as nuclear cascades or spectroscopic excitations require higher levels of theory, the first principles method of choice for modeling mechanical properties as in our case is density functional theory (DFT), in which the electronic structure of materials is calculated from approximations of quantum mechanics [9]. There have been efforts to study the Fe<sub>7</sub>Cr<sub>2</sub>Ni alloy using DFT [5,10–12]. However, as the computational cost of DFT calculations scales cubically with system size, these studies were restricted to supercells containing up to 256 atoms. These supercells are too small to model extended defects such as voids or dislocations. Also, as DFT calculations are computationally expensive, these studies were restricted to limited statistics, which is not ideal for modeling the variations due to chemical complexity of concentrated alloys.

<footnote>*lakshmi.shenoy@warwick.ac.uk</footnote>
Computationally more affordable alternatives for atomistic modeling of materials that allow us to simulate larger length and time scales are classical interatomic potentials. These models have a fixed mathematical form and are fit to experimental and DFT data for the material properties of interest. For metallic systems, the main classical models used are the embedded atom models (EAM) [13], the modified EAM [14], and the Finnis Sinclair model [15]. While there have been many efforts to study alloys using classical interatomic potentials [16–20], these studies are restricted to just a few properties owing to the fixed mathematical form of empirical models. For instance, Bonny *et al.* developed a FeCrNi EAM for plasticity [18] and another FeCrNi EAM for point defects [19] (referred to as EAM-13 throughout this paper), but neither can simulate both types of defects as they were trained for their respective target properties. We would ideally like to have a model that can simulate multiple defects and their interactions.

A successful emerging approach that can address the shortcoming of both DFT and EAM approaches is machine learned interatomic potentials (MLIPs) [21]. MLIPs are trained on a database of first-principles data, generally from DFT. Their highly flexible mathematical forms allow them to interpolate the training data smoothly and so fit the potential energy surface with near DFT accuracy, while at the same time being much cheaper to use for predictions than first-principles methods directly. A major advantage of the MLIPs over classical potentials like EAM is that they can be trained on multiple defects, by including representative configurations in the training database. Moreover, the accuracy and scope of an MLIP can be improved iteratively by appending new training configurations to its training database, chosen either manually using knowledge about the system of interest or using active learning schemes [22].

A range of MLIP approaches have enjoyed considerable success, including kernel methods [23] linear and nonlinear expansions of the potential in a polynomial basis such as the atomic cluster expansion (ACE) [24] or moment tensor potential (MTP) [25] frameworks, and most recently message-passing neural networks [26–28]. In this study, we use the Gaussian approximation potential (GAP) framework, a kernel based method [23]. The GAP model has physically motivated kernel choices and hyperparameters, making it intuitive to understand. It uses Gaussian process regression to fit the potential energy surface; the Bayesian interpretation of this procedure allows us to derive error estimates for GAP predictions, which is a useful feature of this choice of MLIP. The GAP model has been fitted and validated for many elemental materials such as Si [29], W [30], Fe [31], and C [32], demonstrating the effectiveness of this model. More recently, it has also been used for binary alloys [33] and the ternary alloy Ge-Sb-Te [34]. In this paper, we fit and validate a GAP model for Fe<sub>7</sub>Cr<sub>2</sub>Ni. We also extend the GAP descriptors to incorporate the effects of the spin-polarized electronic structure of the alloy in terms of collinear spins associated with the Fe atoms, and demonstrate how this spin GAP model improves predictions compared to the standard GAP model.

This paper is organized as follows: Section II outlines how the training database was assembled, and includes details of

<table><caption>TABLE I. Number of training and testing configurations for each subgroup DBx, comprising either 256-atom bulk or 255-atom vacancy configurations.</caption><thead><tr><th></th><th>Config. type</th><th>Temperature</th><th>Train</th><th>Test</th></tr></thead><tbody><tr><td>DB1</td><td>relaxed bulk</td><td>0 K</td><td>111</td><td>12</td></tr><tr><td>DB2</td><td>sheared bulk</td><td>~0 K</td><td>190</td><td>75</td></tr><tr><td>DB3</td><td>relaxed vacancy</td><td>0 K</td><td>36</td><td>8</td></tr><tr><td>DB4</td><td>vacancy optimisation</td><td>~0 K</td><td>72</td><td>20</td></tr><tr><td>DB5</td><td>bulk MD</td><td>1000/1500/2000 K</td><td>60</td><td>13</td></tr><tr><td>DB6</td><td>vacancy MD</td><td>1000/1500/2000 K</td><td>71</td><td>20</td></tr><tr><td>DB7</td><td>bulk MD quench</td><td>1000 K → 0 K</td><td>38</td><td>4</td></tr><tr><td>DB8</td><td>vacancy MD quench</td><td>1000 K → 0 K</td><td>44</td><td>28</td></tr><tr><td></td><td>element substitutions</td><td>0K</td><td>0</td><td>42</td></tr></tbody></table>

the standard GAP and spin GAP fitting procedures. Section III discusses the main results for bulk and point defect properties of Fe<sub>7</sub>Cr<sub>2</sub>Ni using the spin GAP. The results are validated against DFT and experiments, and compared with EAM and standard GAP predictions.

## II. METHODOLOGY

### A. Training Database

The first step to training a MLIP is to assemble a training database of high accuracy data relevant to properties of interest. Our training database comprises configurations relevant to the bulk, finite temperature, and point defect properties of the alloy. It can be categorized into eight subdatabases (DBx) based on configuration type, as summarized in Table I. Figure 1 shows the kernel principal component analysis (KPCA) of the training database configurations using the SOAP kernel [35], colored by subdatabases. The horizontal axis KPCA component can be interpreted as being roughly correlated with the temperature of the configurations, while the vertical axis KPCA component is correlated with strain.

Perfect-lattice bulk and monovacancy supercells (containing 256 and 255 atoms respectively) with different element

<figure><img src="image_3.png" /><figcaption>FIG. 1. Kernel principal component analysis (KPCA) of training data in the three-species SOAP descriptor space. This approach allows us to identify gaps in the training configuration space; e.g., DB7-8 were included to sample the gap between the DB1-4 and DB5-6 clusters.</figcaption></figure>
Computationally more affordable alternatives for atomistic modeling of materials that allow us to simulate larger length and time scales are classical interatomic potentials. These models have a fixed mathematical form and are fit to experimental and DFT data for the material properties of interest. For metallic systems, the main classical models used are the embedded atom models (EAM) [13], the modified EAM [14], and the Finnis Sinclair model [15]. While there have been many efforts to study alloys using classical interatomic potentials [16–20], these studies are restricted to just a few properties owing to the fixed mathematical form of empirical models. For instance, Bonny *et al.* developed a FeCrNi EAM for plasticity [18] and another FeCrNi EAM for point defects [19] (referred to as EAM-13 throughout this paper), but neither can simulate both types of defects as they were trained for their respective target properties. We would ideally like to have a model that can simulate multiple defects and their interactions.

A successful emerging approach that can address the shortcoming of both DFT and EAM approaches is machine learned interatomic potentials (MLIPs) [21]. MLIPs are trained on a database of first-principles data, generally from DFT. Their highly flexible mathematical forms allow them to interpolate the training data smoothly and so fit the potential energy surface with near DFT accuracy, while at the same time being much cheaper to use for predictions than first-principles methods directly. A major advantage of the MLIPs over classical potentials like EAM is that they can be trained on multiple defects, by including representative configurations in the training database. Moreover, the accuracy and scope of an MLIP can be improved iteratively by appending new training configurations to its training database, chosen either manually using knowledge about the system of interest or using active learning schemes [22].

A range of MLIP approaches have enjoyed considerable success, including kernel methods [23] linear and nonlinear expansions of the potential in a polynomial basis such as the atomic cluster expansion (ACE) [24] or moment tensor potential (MTP) [25] frameworks, and most recently message-passing neural networks [26–28]. In this study, we use the Gaussian approximation potential (GAP) framework, a kernel based method [23]. The GAP model has physically motivated kernel choices and hyperparameters, making it intuitive to understand. It uses Gaussian process regression to fit the potential energy surface; the Bayesian interpretation of this procedure allows us to derive error estimates for GAP predictions, which is a useful feature of this choice of MLIP. The GAP model has been fitted and validated for many elemental materials such as Si [29], W [30], Fe [31], and C [32], demonstrating the effectiveness of this model. More recently, it has also been used for binary alloys [33] and the ternary alloy Ge-Sb-Te [34]. In this paper, we fit and validate a GAP model for Fe<sub>7</sub>Cr<sub>2</sub>Ni. We also extend the GAP descriptors to incorporate the effects of the spin-polarized electronic structure of the alloy in terms of collinear spins associated with the Fe atoms, and demonstrate how this spin GAP model improves predictions compared to the standard GAP model.

This paper is organized as follows: Section II outlines how the training database was assembled, and includes details of

<table><caption>TABLE I. Number of training and testing configurations for each subgroup DBx, comprising either 256-atom bulk or 255-atom vacancy configurations.</caption><thead><tr><th></th><th>Config. type</th><th>Temperature</th><th>Train</th><th>Test</th></tr></thead><tbody><tr><td>DB1</td><td>relaxed bulk</td><td>0 K</td><td>111</td><td>12</td></tr><tr><td>DB2</td><td>sheared bulk</td><td>~0 K</td><td>190</td><td>75</td></tr><tr><td>DB3</td><td>relaxed vacancy</td><td>0 K</td><td>36</td><td>8</td></tr><tr><td>DB4</td><td>vacancy optimisation</td><td>~0 K</td><td>72</td><td>20</td></tr><tr><td>DB5</td><td>bulk MD</td><td>1000/1500/2000 K</td><td>60</td><td>13</td></tr><tr><td>DB6</td><td>vacancy MD</td><td>1000/1500/2000 K</td><td>71</td><td>20</td></tr><tr><td>DB7</td><td>bulk MD quench</td><td>1000 K → 0 K</td><td>38</td><td>4</td></tr><tr><td>DB8</td><td>vacancy MD quench</td><td>1000 K → 0 K</td><td>44</td><td>28</td></tr><tr><td></td><td>element substitutions</td><td>0K</td><td>0</td><td>42</td></tr></tbody></table>

the standard GAP and spin GAP fitting procedures. Section III discusses the main results for bulk and point defect properties of Fe<sub>7</sub>Cr<sub>2</sub>Ni using the spin GAP. The results are validated against DFT and experiments, and compared with EAM and standard GAP predictions.

## II. METHODOLOGY

### A. Training Database

The first step to training a MLIP is to assemble a training database of high accuracy data relevant to properties of interest. Our training database comprises configurations relevant to the bulk, finite temperature, and point defect properties of the alloy. It can be categorized into eight subdatabases (DBx) based on configuration type, as summarized in Table I. Figure 1 shows the kernel principal component analysis (KPCA) of the training database configurations using the SOAP kernel [35], colored by subdatabases. The horizontal axis KPCA component can be interpreted as being roughly correlated with the temperature of the configurations, while the vertical axis KPCA component is correlated with strain.

Perfect-lattice bulk and monovacancy supercells (containing 256 and 255 atoms respectively) with different element

<figure><img src="image_3.png" /><figcaption>FIG. 1. Kernel principal component analysis (KPCA) of training data in the three-species SOAP descriptor space. This approach allows us to identify gaps in the training configuration space; e.g., DB7-8 were included to sample the gap between the DB1-4 and DB5-6 clusters.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIG. 2. Protocol for generating bulk and vacancy training structures that sample different elemental configurations, with realistic short-range ordering derived from effective medium (EM) theory. In this case, the DFT relaxed off-lattice displacements shown by the red arrows of step 3 are for AFM layering along the $\hat{x}$ axis (AFM-X).</figcaption></figure>

distributions that capture relevant short-range ordering were generated using an effective medium theory (EM) as described in Sec. II B. These were fully relaxed as shown in the schematic of Fig. 2, using DFT calculations as described in Sec. II C, to generate the training subdatabases DB1 and DB3 respectively. Random strains were applied to a subset of DB1 samples, to generate DB2, aimed at providing reference material, which can inform an accurate description of the alloy's elastic behavior. Snapshots from the DFT geometry optimisation trajectory of the EM monovacancy structures were compiled to form DB4, to provide details of the potential energy surface close to 0 K.

The EM structures were generated to capture short-range ordering at temperatures of 600 K, 1000 K, 1500 K, and 2000 K, respectively. We ran ab initio molecular dynamics (MD) for the 1000–2000 K configurations from DB1 and DB3 at the temperatures they were generated for, and uncorrelated snapshots from these MD trajectories were compiled to form DB5 and DB6, respectively. In Fig. 1, we can see that the high-temperature DB5-6 samples sit quite far apart from the low-temperature DB1-4 samples along the first KPCA component. To aid the GAP model with better interpolation of training data, we need to populate such gaps between training data clusters along the principal axes of the kernel space. To this end, we extend the database by including snapshots from the geometry optimisation of a few bulk and monovacancy MD structures, comprising DB7 and DB8, respectively. These configurations are seen to sample the space between low- and high-temperature clusters in the KPCA plot Fig. 1, and their inclusion decreased the errors of our models. The GAP potentials (both with and without the Fe spin extensions) were trained on this combined database DB1-8 as described in Sec. II D.

## B. Effective Medium Theory

Past DFT studies on FeCrNi have modelled the alloy using special quasirandom structures (SQS) [5,10,11]. A limitation of using SQS for magnetic materials is that in aiming to produce local environments reflective of the random alloy, they neglect atomic short-range order. It would be more effective to train a model using configurations in which temperature-dependent short-range order is present, especially in the case of point defects where the energetics are correlated to the concentration of species decorating the defect [5].

To obtain atomic configurations with physically motivated atomic short-range order, we draw samples from equilibrated, lattice-based Monte Carlo (MC) simulations, as depicted by the first step in Fig. 2. A configuration is described by a set of lattice site occupancies {$\xi_{i\alpha}$}, where $\xi_{i\alpha} = 1$ if site $i$ is occupied by an atom of chemical species $\alpha$, and 0 otherwise. The internal energy is described by a simple pairwise Bragg-Williams Hamiltonian [36,37],

$$H = \sum_{i\alpha; j\alpha'} V_{i\alpha; j\alpha'} \xi_{i\alpha} \xi_{j\alpha'}, \quad (1)$$

where the atom-atom interchange parameter $V_{i\alpha; j\alpha'}$ describes the energy associated with an atom of species $i$ on site $\alpha$ interacting with an atom of species $j$ on site $\alpha'$. We assume these interchange parameters are isotropic, homogeneous, and have finite range, which simplifies Eq. (1).

The $V_{i\alpha; j\alpha'}$ are obtained using the $S^{(2)}$ theory for multicomponent alloys [38], the details and implementation of which have been discussed extensively in earlier papers [39–42]. The theory uses the Korringa-Kohn-Rostoker (KKR) formulation of DFT, with disorder described via the coherent potential approximation (CPA) [43], producing an effective medium representing the electronic structure of the disordered alloy. Our method captures the effects of an ordered magnetic state, in our case the single-layer antiferromagnetic (AFM) state, on the short-range-order state [41]. It is also possible to treat vacancies in this formalism at low concentrations by including a chemical species with no associated electrons, and no nuclear charge.

In this paper, we used the HUTSEPOT code [44] to construct the self-consistent potentials of DFT within the KKR-CPA formalism. The interchange parameters are fitted to the first four coordination shells, and they are tabulated in Appendix A for reference. The MC simulations used a cells of 256 atoms with periodic boundary conditions applied in all three directions. Simulations were equilibrated at the desired temperature and then decorrelated samples were drawn 25 600 MC steps apart, i.e., 100 MC steps per lattice site. The samples drawn from high-temperature simulations will have little to no short-range order, while those drawn at lower temperatures will have a degree of short-range order present. These lattice-based EM configurations are then relaxed with DFT using setting described in the next section, to get the cell distortion and off-lattice atomic displacements as depicted in the last step of Fig. 2.
<figure><img src="image_4.png" /><figcaption>FIG. 2. Protocol for generating bulk and vacancy training structures that sample different elemental configurations, with realistic short-range ordering derived from effective medium (EM) theory. In this case, the DFT relaxed off-lattice displacements shown by the red arrows of step 3 are for AFM layering along the $\hat{x}$ axis (AFM-X).</figcaption></figure>

distributions that capture relevant short-range ordering were generated using an effective medium theory (EM) as described in Sec. II B. These were fully relaxed as shown in the schematic of Fig. 2, using DFT calculations as described in Sec. II C, to generate the training subdatabases DB1 and DB3 respectively. Random strains were applied to a subset of DB1 samples, to generate DB2, aimed at providing reference material, which can inform an accurate description of the alloy's elastic behavior. Snapshots from the DFT geometry optimisation trajectory of the EM monovacancy structures were compiled to form DB4, to provide details of the potential energy surface close to 0 K.

The EM structures were generated to capture short-range ordering at temperatures of 600 K, 1000 K, 1500 K, and 2000 K, respectively. We ran ab initio molecular dynamics (MD) for the 1000–2000 K configurations from DB1 and DB3 at the temperatures they were generated for, and uncorrelated snapshots from these MD trajectories were compiled to form DB5 and DB6, respectively. In Fig. 1, we can see that the high-temperature DB5-6 samples sit quite far apart from the low-temperature DB1-4 samples along the first KPCA component. To aid the GAP model with better interpolation of training data, we need to populate such gaps between training data clusters along the principal axes of the kernel space. To this end, we extend the database by including snapshots from the geometry optimisation of a few bulk and monovacancy MD structures, comprising DB7 and DB8, respectively. These configurations are seen to sample the space between low- and high-temperature clusters in the KPCA plot Fig. 1, and their inclusion decreased the errors of our models. The GAP potentials (both with and without the Fe spin extensions) were trained on this combined database DB1-8 as described in Sec. II D.

## B. Effective Medium Theory

Past DFT studies on FeCrNi have modelled the alloy using special quasirandom structures (SQS) [5,10,11]. A limitation of using SQS for magnetic materials is that in aiming to produce local environments reflective of the random alloy, they neglect atomic short-range order. It would be more effective to train a model using configurations in which temperature-dependent short-range order is present, especially in the case of point defects where the energetics are correlated to the concentration of species decorating the defect [5].

To obtain atomic configurations with physically motivated atomic short-range order, we draw samples from equilibrated, lattice-based Monte Carlo (MC) simulations, as depicted by the first step in Fig. 2. A configuration is described by a set of lattice site occupancies {$\xi_{i\alpha}$}, where $\xi_{i\alpha} = 1$ if site $i$ is occupied by an atom of chemical species $\alpha$, and 0 otherwise. The internal energy is described by a simple pairwise Bragg-Williams Hamiltonian [36,37],

$$H = \sum_{i\alpha; j\alpha'} V_{i\alpha; j\alpha'} \xi_{i\alpha} \xi_{j\alpha'}, \quad (1)$$

where the atom-atom interchange parameter $V_{i\alpha; j\alpha'}$ describes the energy associated with an atom of species $i$ on site $\alpha$ interacting with an atom of species $j$ on site $\alpha'$. We assume these interchange parameters are isotropic, homogeneous, and have finite range, which simplifies Eq. (1).

The $V_{i\alpha; j\alpha'}$ are obtained using the $S^{(2)}$ theory for multicomponent alloys [38], the details and implementation of which have been discussed extensively in earlier papers [39–42]. The theory uses the Korringa-Kohn-Rostoker (KKR) formulation of DFT, with disorder described via the coherent potential approximation (CPA) [43], producing an effective medium representing the electronic structure of the disordered alloy. Our method captures the effects of an ordered magnetic state, in our case the single-layer antiferromagnetic (AFM) state, on the short-range-order state [41]. It is also possible to treat vacancies in this formalism at low concentrations by including a chemical species with no associated electrons, and no nuclear charge.

In this paper, we used the HUTSEPOT code [44] to construct the self-consistent potentials of DFT within the KKR-CPA formalism. The interchange parameters are fitted to the first four coordination shells, and they are tabulated in Appendix A for reference. The MC simulations used a cells of 256 atoms with periodic boundary conditions applied in all three directions. Simulations were equilibrated at the desired temperature and then decorrelated samples were drawn 25 600 MC steps apart, i.e., 100 MC steps per lattice site. The samples drawn from high-temperature simulations will have little to no short-range order, while those drawn at lower temperatures will have a degree of short-range order present. These lattice-based EM configurations are then relaxed with DFT using setting described in the next section, to get the cell distortion and off-lattice atomic displacements as depicted in the last step of Fig. 2.
### C. DFT Settings

The Vienna Ab initio Software Package (VASP) [45–47] was used to compute ground-state DFT structure, energy, atomic forces, and cell stresses for the training database. MLIPs can only be as accurate as the data they are trained on, so the choice of DFT settings cap the achievable accuracy. A key choice when setting up DFT calculations is which exchange-correlation (XC) functional to use. For structural properties of transition metals, the Perdew–Burke–Ernzerhof [48] flavor of generalized gradient approximation (GGA) functionals is well established in the literature as a good choice, and so this is the XC functional we use. However, for more complex phenomena involving electronic excitations, we would require beyond-GGA approaches such as hybrid exchange, or even beyond-DFT methods such as the GW approach, and so we cannot hope to model such phenomena using a MLIP trained on GGA-based DFT.

Depending on the range of atomic interactions that we hope to model, the choice of pseudopotential would also affect the DFT accuracy. For example, modeling nuclear cascades would require very short-range interactions that benefit from all-electron DFT codes or semicore pseudopotentials [49]. However, since we are targeting mechanical properties, the key players are the outer shell electrons of the alloying elements, so the accuracy we require is not compromised by using softer pseudopotentials for the core electrons. We use the standard VASP pseudopotentials based on the projector augmented wave (PAW) method [50], with eight, six, and ten valence electrons for Fe, Cr, and Ni respectively.

As Fe and Ni are magnetic materials, and Cr has shown a spin-density wave at low temperatures [51,52], it is essential to turn on the spin setting for more accurate DFT calculations of this alloy. Piochaud et al. reported that their tests with noncollinear spins generally relaxed to collinear arrangements [5], so we constrain our systems to collinear spin. Also, while austenitic steel is paramagnetic (PM) at operation temperatures of nuclear reactors, single-layer antiferromagnetic (AFM) ordering was found to be the most stable ground state [5]. To ensure we get ground-state properties right, we impose an AFM ordering for DFT relaxations of our training samples. Initial magnetic moment magnitudes are overestimated to allow for relaxation, and are set at 3.0, 2.0, and 1.0 $\mu_B$ for Fe, Cr, and Ni respectively. These atomic magnetic moments evolve during DFT electronic minimisation, and are computed in VASP by integrating the spin density within a sphere centered around the atom of interest.

A Monkhorst-Pack k-point mesh [53] of size 3 × 3 × 3 was used to sample the Brillouin zone for all the training samples. Training data was computed with a plane-wave energy cutoff of 600 eV and an electronic self-consistency criterion of 10⁻⁷ eV. Preliminary steps, i.e., initial geometry optimisations and ab initio MD trajectories, were done using a lower self-consistency criterion of 10⁻⁴ eV. Electronic minimisation were performed using the preconditioned RMM-DIIS algorithm, while ionic position optimisations were done using the conjugate gradient method. The ab initio MD trajectories were run for the NVT ensemble, using the Nose-Hoover thermostat [54,55] and timesteps of either 5 fs (for T = 1000 K, 1500 K) or 2 fs (for T = 2000 K). The initial ~200 steps (exact number

depending on when the total energy stabilizes) of the trajectories were discarded, and uncorrelated samples were taken from the remaining 1000–1500 steps of each trajectory. These were then evaluated with the higher precision DFT settings before being included in the training database.

### D. Potential Fitting

This section outlines how the GAP fits the Born-Oppenheimer potential energy surface (PES) by mapping the DFT training structures to their DFT energies, forces, and stresses. As with other interatomic potentials, the GAP assumes that the energy $E_{N_c}$ of a configuration with $N_c$ atoms can be written as

$$E_{N_c} = \sum_{i=1}^{N_c} \epsilon(\hat{\mathbf{q}}_i) \qquad (2)$$

where $\epsilon(\hat{\mathbf{q}}_i)$ is the energy contribution from atom $i$ in the configuration, and $\hat{\mathbf{q}}_i$ is a descriptor vector that captures the unique features of the atom $i$ environment. The atomic energies are modelled as a linear combination of kernels [56] given by

$$\epsilon(\hat{\mathbf{q}}_i) = \sum_{j}^{M} \alpha_j K(\hat{\mathbf{q}}_j, \hat{\mathbf{q}}_i) = \mathbf{k}_i^T \alpha \qquad (3)$$

where the kernel basis functions $K(\hat{\mathbf{q}}_j, \hat{\mathbf{q}}_i)$ capture the similarity between environment $i$ and the $M$ environments of the basis set. As GAP training databases are generally large (e.g., $N \sim 160 000$ atomic environments in this paper), the usual strategy is to use only a small portion $M \ll N$ of the training set for the basis set. This makes the fit computationally tractable, while still retaining accuracy as the training database contains many groups of similar configurations that can be sparsified without much loss of information. These $M$ sparse points are selected by CUR decomposition, which searches for maximally dissimilar environments [29], ensuring that the chosen basis is representative of the variety across the full training database.

Moving to the left-hand side of Eq. (3), DFT does not give direct access to the atomic energies $\epsilon$, but rather, it gives us total energies $E_{N_c}$ of the training configurations, along with atomic forces and virial stresses. As forces and stress are properties that can be derived from atomic energies (e.g., $\partial/\partial x$ for an atomic force along the Cartesian $x$ direction), we can back-calculate the unknown $N$ component vector $\epsilon$ of atomic energies, from the known $D$ component vector $y$ of DFT data, by the relation $y = L^T \epsilon$, where $L$ is a linear differential operator $L$ [30]. This allows us to derive the covariance kernel of the unknown target atomic energies ($K_{NN}$) from the covariance kernel of the known DFT training data ($K_{DD}$) using the relation $K_{DD} = L^T K_{NN} L$. The sparse ridge-regression solution to Eq. (3) regularized by a diagonal tolerance matrix $\Lambda = \sigma_v^2 \mathbf{I}$ is

$$\alpha = [K_{MM} + K_{MN} L \Lambda^{-1} L K_{NM}]^{-1} K_{MN} L \Lambda^{-1} y. \qquad (4)$$

With these optimized coefficients $\alpha$, we can make predictions for new atomic environments $\hat{\mathbf{q}}_*$ using Eq. (3). In this paper, we use the smooth overlap of atomic positions (SOAP)
### C. DFT Settings

The Vienna Ab initio Software Package (VASP) [45–47] was used to compute ground-state DFT structure, energy, atomic forces, and cell stresses for the training database. MLIPs can only be as accurate as the data they are trained on, so the choice of DFT settings cap the achievable accuracy. A key choice when setting up DFT calculations is which exchange-correlation (XC) functional to use. For structural properties of transition metals, the Perdew–Burke–Ernzerhof [48] flavor of generalized gradient approximation (GGA) functionals is well established in the literature as a good choice, and so this is the XC functional we use. However, for more complex phenomena involving electronic excitations, we would require beyond-GGA approaches such as hybrid exchange, or even beyond-DFT methods such as the GW approach, and so we cannot hope to model such phenomena using a MLIP trained on GGA-based DFT.

Depending on the range of atomic interactions that we hope to model, the choice of pseudopotential would also affect the DFT accuracy. For example, modeling nuclear cascades would require very short-range interactions that benefit from all-electron DFT codes or semicore pseudopotentials [49]. However, since we are targeting mechanical properties, the key players are the outer shell electrons of the alloying elements, so the accuracy we require is not compromised by using softer pseudopotentials for the core electrons. We use the standard VASP pseudopotentials based on the projector augmented wave (PAW) method [50], with eight, six, and ten valence electrons for Fe, Cr, and Ni respectively.

As Fe and Ni are magnetic materials, and Cr has shown a spin-density wave at low temperatures [51,52], it is essential to turn on the spin setting for more accurate DFT calculations of this alloy. Piochaud et al. reported that their tests with noncollinear spins generally relaxed to collinear arrangements [5], so we constrain our systems to collinear spin. Also, while austenitic steel is paramagnetic (PM) at operation temperatures of nuclear reactors, single-layer antiferromagnetic (AFM) ordering was found to be the most stable ground state [5]. To ensure we get ground-state properties right, we impose an AFM ordering for DFT relaxations of our training samples. Initial magnetic moment magnitudes are overestimated to allow for relaxation, and are set at 3.0, 2.0, and 1.0 $\mu_B$ for Fe, Cr, and Ni respectively. These atomic magnetic moments evolve during DFT electronic minimisation, and are computed in VASP by integrating the spin density within a sphere centered around the atom of interest.

A Monkhorst-Pack k-point mesh [53] of size 3 × 3 × 3 was used to sample the Brillouin zone for all the training samples. Training data was computed with a plane-wave energy cutoff of 600 eV and an electronic self-consistency criterion of 10⁻⁷ eV. Preliminary steps, i.e., initial geometry optimisations and ab initio MD trajectories, were done using a lower self-consistency criterion of 10⁻⁴ eV. Electronic minimisation were performed using the preconditioned RMM-DIIS algorithm, while ionic position optimisations were done using the conjugate gradient method. The ab initio MD trajectories were run for the NVT ensemble, using the Nose-Hoover thermostat [54,55] and timesteps of either 5 fs (for T = 1000 K, 1500 K) or 2 fs (for T = 2000 K). The initial ~200 steps (exact number

depending on when the total energy stabilizes) of the trajectories were discarded, and uncorrelated samples were taken from the remaining 1000–1500 steps of each trajectory. These were then evaluated with the higher precision DFT settings before being included in the training database.

### D. Potential Fitting

This section outlines how the GAP fits the Born-Oppenheimer potential energy surface (PES) by mapping the DFT training structures to their DFT energies, forces, and stresses. As with other interatomic potentials, the GAP assumes that the energy $E_{N_c}$ of a configuration with $N_c$ atoms can be written as

$$E_{N_c} = \sum_{i=1}^{N_c} \epsilon(\hat{\mathbf{q}}_i) \qquad (2)$$

where $\epsilon(\hat{\mathbf{q}}_i)$ is the energy contribution from atom $i$ in the configuration, and $\hat{\mathbf{q}}_i$ is a descriptor vector that captures the unique features of the atom $i$ environment. The atomic energies are modelled as a linear combination of kernels [56] given by

$$\epsilon(\hat{\mathbf{q}}_i) = \sum_{j}^{M} \alpha_j K(\hat{\mathbf{q}}_j, \hat{\mathbf{q}}_i) = \mathbf{k}_i^T \alpha \qquad (3)$$

where the kernel basis functions $K(\hat{\mathbf{q}}_j, \hat{\mathbf{q}}_i)$ capture the similarity between environment $i$ and the $M$ environments of the basis set. As GAP training databases are generally large (e.g., $N \sim 160 000$ atomic environments in this paper), the usual strategy is to use only a small portion $M \ll N$ of the training set for the basis set. This makes the fit computationally tractable, while still retaining accuracy as the training database contains many groups of similar configurations that can be sparsified without much loss of information. These $M$ sparse points are selected by CUR decomposition, which searches for maximally dissimilar environments [29], ensuring that the chosen basis is representative of the variety across the full training database.

Moving to the left-hand side of Eq. (3), DFT does not give direct access to the atomic energies $\epsilon$, but rather, it gives us total energies $E_{N_c}$ of the training configurations, along with atomic forces and virial stresses. As forces and stress are properties that can be derived from atomic energies (e.g., $\partial/\partial x$ for an atomic force along the Cartesian $x$ direction), we can back-calculate the unknown $N$ component vector $\epsilon$ of atomic energies, from the known $D$ component vector $y$ of DFT data, by the relation $y = L^T \epsilon$, where $L$ is a linear differential operator $L$ [30]. This allows us to derive the covariance kernel of the unknown target atomic energies ($K_{NN}$) from the covariance kernel of the known DFT training data ($K_{DD}$) using the relation $K_{DD} = L^T K_{NN} L$. The sparse ridge-regression solution to Eq. (3) regularized by a diagonal tolerance matrix $\Lambda = \sigma_v^2 \mathbf{I}$ is

$$\alpha = [K_{MM} + K_{MN} L \Lambda^{-1} L K_{NM}]^{-1} K_{MN} L \Lambda^{-1} y. \qquad (4)$$

With these optimized coefficients $\alpha$, we can make predictions for new atomic environments $\hat{\mathbf{q}}_*$ using Eq. (3). In this paper, we use the smooth overlap of atomic positions (SOAP)
<table><caption>TABLE II. Parameters for the collinear-spin GAP.</caption><thead><tr><th>Parameter</th><th>Two-body</th><th>Many-body</th></tr></thead><tbody><tr><td>Kernel choice</td><td>Gaussian</td><td>SOAP</td></tr><tr><td>Radial cutoff</td><td>6.0 Å</td><td>6.0 Å</td></tr><tr><td>Cutoff transition width</td><td>1.0 Å</td><td>1.0 Å</td></tr><tr><td>δ</td><td>1.0 eV</td><td>0.1 eV</td></tr><tr><td>Sparse points</td><td>100</td><td>800</td></tr><tr><td>θ</td><td>1.0</td><td></td></tr><tr><td>n<sub>max</sub>, l<sub>max</sub></td><td></td><td>10,10</td></tr><tr><td>ζ</td><td></td><td>4</td></tr><tr><td>σ<sub>atom</sub></td><td></td><td>1.0 Å</td></tr></tbody></table>

descriptor

$$ \mathbf{q}_i = \sum_{m=-l}^{l} (c_{nlm}^j)^* c_{nlm}^j ; \quad \hat{\mathbf{q}}_i = \frac{\mathbf{q}_i}{|\mathbf{q}_i|} \quad (5) $$

where $c_{nlm}$ are expansion coefficients for the atomic neighbor density power spectrum [35]. For computational tractability, this sum is truncated at radial-basis functions up to order $n < n_{max}$ and spherical harmonics up to degree $l < l_{max}$. The SOAP descriptor is suitable for describing atomic environments since it is translationally, rotationally, and permutationally invariant [35]. The corresponding covariance kernel is

$$ K(\hat{\mathbf{q}}_i, \hat{\mathbf{q}}_j) = \delta^2|\hat{\mathbf{q}}_i \cdot \hat{\mathbf{q}}_j|^\zeta \quad (6) $$

where $\delta$ is a factor that corresponds to the spread in energy of the training dataset, and $\zeta$ is tuned to amplify variations in the dot product, to pick out fine features.

While the SOAP kernel is good at capturing smoothly varying details of the PES, it is not as effective in capturing sharp features as well, such as exchange repulsion when atoms get close to each other, or phase transitions [30]. It has been found to be useful to use multiple kernels targeting different scales, so that each kernel does not have to simultaneously fit vastly different features [29]. We hence use simple two-body kernels as well, to roughly capture the main trends of the PES, so that the SOAP kernels can focus on finer details. The two-body descriptor is just the pairwise separation between atoms, and we use the squared-exponential form

$$ K(\hat{q}_i, \hat{q}_j) = \delta^2 \exp \left( -\frac{1}{2} \left( \frac{\hat{q}_i - \hat{q}_j}{\theta_j} \right)^2 \right) \quad (7) $$

for the two-body kernel, where $\theta_j$ is a lengthscale hyperparameter. For a three-species alloy, we require six two-body kernels (one for each pair of species), and three SOAP kernels (one for each species). The kernel hyperparameters of Eqs. (5), (6), and (7) that we use are summarized in Table II. In addition these, we $\sigma_v^{\text{energy}} = 0.001 \text{ eV/atom}$, $\sigma_v^{\text{force}} = 0.1 \text{ eV/Å}$, and $\sigma_v^{\text{virial}} = 0.05 \text{ eV/atom}$ for the regularization described in Eq. (4), values chosen based on the accuracy of our DFT training data. The MPI parallel gap_fit code from the QUIP software package [57,58] was used to fit the GAP on the full dataset of Table I.

The cumulative distribution functions of the training and test energy errors of the resulting nonmagnetic GAP are

<figure><img src="image_5.png" /><figcaption>FIG. 3. Cumulative distribution functions of energy training (solid lines) and test (dashed lines) errors for the standard SOAP + GAP model (blue), a GAP model using compressed SOAP descriptors (orange) and a GAP model with an additional species included to differentiate between spin-up and spin-down Fe sites (spin GAP; green). Compression does not significantly reduce accuracy (orange vs blue), and the spin GAP model (green) performs significantly better.</figcaption></figure>

shown in Fig. 3. Similar plots for the force and stress errors have been included in Appendix B. The GAP has energy training and test error within 2 meV/atom and 6 meV/atom respectively, which are comparable to the differences in DFT energies due to choice of AFM layering direction. The train and test force errors are below 200 meV/Å for low-temperature configurations and 400 meV/Å for high-temperature configurations. The relatively large force errors for low-temperature configurations are because the GAP cannot fit the directional spin-induced face-centered tetragonal (FCT) ground-state structure, as the GAP does not have spin degrees of freedom. The stress training and test errors are within 2 GPa and 4 GPa respectively, which are small as they are less than 2% of the typical elastic constant values that are about ~200GPa.

## E. Spin Model

As the SOAP and two-body descriptors do not have any information on spin, the standard GAP is spin-agnostic and cannot capture the tetrahedral deformation of AFM layering, and this limits the accuracy of its predictions. To address this, we introduce a fictitious fourth species to the GAP model, such that we have different descriptors for up spin (Fe ↑) and down spin (Fe ↓). This expanded descriptor set gives the model the flexibility to learn the tetrahedral deformation of AFM layering, and could be extended to the PM phase too. We henceforth refer to this 4-species potential as the spin GAP.

As the SOAP descriptor length scales quadratically with number of species, slowing down training and evaluation times prohibitively, we needed to use descriptor compression methods to make the training computationally tractable, such as those introduced in [59]. We use combined radial and species compression with tensor product coupling across 300
<table><caption>TABLE II. Parameters for the collinear-spin GAP.</caption><thead><tr><th>Parameter</th><th>Two-body</th><th>Many-body</th></tr></thead><tbody><tr><td>Kernel choice</td><td>Gaussian</td><td>SOAP</td></tr><tr><td>Radial cutoff</td><td>6.0 Å</td><td>6.0 Å</td></tr><tr><td>Cutoff transition width</td><td>1.0 Å</td><td>1.0 Å</td></tr><tr><td>δ</td><td>1.0 eV</td><td>0.1 eV</td></tr><tr><td>Sparse points</td><td>100</td><td>800</td></tr><tr><td>θ</td><td>1.0</td><td></td></tr><tr><td>n<sub>max</sub>, l<sub>max</sub></td><td></td><td>10,10</td></tr><tr><td>ζ</td><td></td><td>4</td></tr><tr><td>σ<sub>atom</sub></td><td></td><td>1.0 Å</td></tr></tbody></table>

descriptor

$$ \mathbf{q}_i = \sum_{m=-l}^{l} (c_{nlm}^j)^* c_{nlm}^j ; \quad \hat{\mathbf{q}}_i = \frac{\mathbf{q}_i}{|\mathbf{q}_i|} \quad (5) $$

where $c_{nlm}$ are expansion coefficients for the atomic neighbor density power spectrum [35]. For computational tractability, this sum is truncated at radial-basis functions up to order $n < n_{max}$ and spherical harmonics up to degree $l < l_{max}$. The SOAP descriptor is suitable for describing atomic environments since it is translationally, rotationally, and permutationally invariant [35]. The corresponding covariance kernel is

$$ K(\hat{\mathbf{q}}_i, \hat{\mathbf{q}}_j) = \delta^2|\hat{\mathbf{q}}_i \cdot \hat{\mathbf{q}}_j|^\zeta \quad (6) $$

where $\delta$ is a factor that corresponds to the spread in energy of the training dataset, and $\zeta$ is tuned to amplify variations in the dot product, to pick out fine features.

While the SOAP kernel is good at capturing smoothly varying details of the PES, it is not as effective in capturing sharp features as well, such as exchange repulsion when atoms get close to each other, or phase transitions [30]. It has been found to be useful to use multiple kernels targeting different scales, so that each kernel does not have to simultaneously fit vastly different features [29]. We hence use simple two-body kernels as well, to roughly capture the main trends of the PES, so that the SOAP kernels can focus on finer details. The two-body descriptor is just the pairwise separation between atoms, and we use the squared-exponential form

$$ K(\hat{q}_i, \hat{q}_j) = \delta^2 \exp \left( -\frac{1}{2} \left( \frac{\hat{q}_i - \hat{q}_j}{\theta_j} \right)^2 \right) \quad (7) $$

for the two-body kernel, where $\theta_j$ is a lengthscale hyperparameter. For a three-species alloy, we require six two-body kernels (one for each pair of species), and three SOAP kernels (one for each species). The kernel hyperparameters of Eqs. (5), (6), and (7) that we use are summarized in Table II. In addition these, we $\sigma_v^{\text{energy}} = 0.001 \text{ eV/atom}$, $\sigma_v^{\text{force}} = 0.1 \text{ eV/Å}$, and $\sigma_v^{\text{virial}} = 0.05 \text{ eV/atom}$ for the regularization described in Eq. (4), values chosen based on the accuracy of our DFT training data. The MPI parallel gap_fit code from the QUIP software package [57,58] was used to fit the GAP on the full dataset of Table I.

The cumulative distribution functions of the training and test energy errors of the resulting nonmagnetic GAP are

<figure><img src="image_5.png" /><figcaption>FIG. 3. Cumulative distribution functions of energy training (solid lines) and test (dashed lines) errors for the standard SOAP + GAP model (blue), a GAP model using compressed SOAP descriptors (orange) and a GAP model with an additional species included to differentiate between spin-up and spin-down Fe sites (spin GAP; green). Compression does not significantly reduce accuracy (orange vs blue), and the spin GAP model (green) performs significantly better.</figcaption></figure>

shown in Fig. 3. Similar plots for the force and stress errors have been included in Appendix B. The GAP has energy training and test error within 2 meV/atom and 6 meV/atom respectively, which are comparable to the differences in DFT energies due to choice of AFM layering direction. The train and test force errors are below 200 meV/Å for low-temperature configurations and 400 meV/Å for high-temperature configurations. The relatively large force errors for low-temperature configurations are because the GAP cannot fit the directional spin-induced face-centered tetragonal (FCT) ground-state structure, as the GAP does not have spin degrees of freedom. The stress training and test errors are within 2 GPa and 4 GPa respectively, which are small as they are less than 2% of the typical elastic constant values that are about ~200GPa.

## E. Spin Model

As the SOAP and two-body descriptors do not have any information on spin, the standard GAP is spin-agnostic and cannot capture the tetrahedral deformation of AFM layering, and this limits the accuracy of its predictions. To address this, we introduce a fictitious fourth species to the GAP model, such that we have different descriptors for up spin (Fe ↑) and down spin (Fe ↓). This expanded descriptor set gives the model the flexibility to learn the tetrahedral deformation of AFM layering, and could be extended to the PM phase too. We henceforth refer to this 4-species potential as the spin GAP.

As the SOAP descriptor length scales quadratically with number of species, slowing down training and evaluation times prohibitively, we needed to use descriptor compression methods to make the training computationally tractable, such as those introduced in [59]. We use combined radial and species compression with tensor product coupling across 300
<table><caption>TABLE III. Lattice parameters and elastic constants of the FeCrNi alloy computed using DFT, GAP, spin GAP, and EAM-13 at 0 K. For comparison, DFT results for AFM Fe at 0 K and experimental values for finite temperatures have also been included. Dashes indicate that the value is equal to the previous value by symmetry.</caption><thead><tr><th></th><th>DFT-Fe</th><th>DFT</th><th>GAP</th><th>Spin GAP</th><th>EAM13</th><th>Expt (&gt; 0 K)</th></tr></thead><tbody><tr><td>Phase</td><td>FCT</td><td>FCT</td><td>FCO</td><td>FCT</td><td>FCC</td><td>~FCC <sup>a</sup></td></tr><tr><td>a (Å)</td><td>3.42</td><td>3.50</td><td>3.49</td><td>3.50</td><td>3.55</td><td>3.52–3.59 <sup>a</sup></td></tr><tr><td>b (Å)</td><td>–</td><td>–</td><td>3.53</td><td>–</td><td>–</td><td>–</td></tr><tr><td>c (Å)</td><td>3.68</td><td>3.58</td><td>3.56</td><td>3.58</td><td>–</td><td>–</td></tr><tr><td>C<sub>11</sub>, C<sub>22</sub>, C<sub>33</sub> (GPa)</td><td>321, –, 254</td><td>259, –, 249</td><td>299, 243, 232</td><td>272, –, 268</td><td>333, –</td><td>204–226<sup>b</sup></td></tr><tr><td>C<sub>12</sub>, C<sub>13</sub>, C<sub>23</sub> (GPa)</td><td>246, 96, –</td><td>164, 124, –</td><td>175, 163, 151</td><td>181, 126, –</td><td>196, –, –</td><td>132–134<sup>b</sup></td></tr><tr><td>C<sub>44</sub>, C<sub>55</sub>, C<sub>66</sub> (GPa)</td><td>166, –, 262</td><td>147, –, 191</td><td>155, 171, 186</td><td>148, –, 188</td><td>161, –, –</td><td>111–122<sup>b</sup></td></tr></tbody></table>

mixed radial and species channels. To check the effects of compression, a compressed version of the three-species GAP was also fit. Figure 3 shows that the compressed GAP has comparable errors to the standard GAP, verifying that there is no significant loss of accuracy when compression is used for this system. The spin GAP, which also uses compressed descriptors, has significantly improved training and test errors over those of the standard GAP. Typical energy errors decrease by 1 meV/atom, force errors by about 50 meV/Å, and stress training errors decrease by about 0.2 GPa.

### III. RESULTS

#### A. Bulk properties

The ground-state bulk properties of the FeCrNi alloy predicted by the GAP and the spin GAP are summarized in Table III, along with the corresponding DFT, EAM, and experimental values for comparison, as well as DFT for pure Fe in an AFM spin state. The DFT results for the ground-state structure of both AFM-Fe and AFM-FeCrNi are FCT, with the longer lattice constant corresponding to the AFM layering direction. This tetragonal distortion is more significant in pure AFM Fe (0.26 Å) compared to the alloy (0.08 Å), as the high concentration of Cr and Ni in the alloy lattice reduce the magnetism-induced geometry effects. Relaxation of EM structures with EAM-13 gives a lattice constant within the DFT and experimental range; however, EAM-13 fails to predict the lattice distortion, leading to a symmetric face centered cubic (FCC) lattice instead, as seen from the first set of bars in Fig. 4(a).

In comparison, the standard GAP predicts that a distorted structure is more energetically favorable than a symmetric FCC structure, as it is trained on relaxed DFT structures that have FCT lattices with zero atomic forces. However, since the standard GAP has no spin degrees of freedom, it does not have the framework to process the spin inputs required to elongate one direction over the other two. Instead, the GAP training correlates the distortion to minor differences in the alloy composition along the different axes in the relaxed training set geometries, and so predicts a face-centred orthorhombic (FCO) relaxed lattice based on composition along the different axes of test structures. As seen from the second set of bars in Fig. 4(a), the three lattice constants of the GAP relaxed FCO structure are equally spaced in a similar range to that

of the two DFT relaxed FCT lattice constants. The spin GAP, on the other hand, has a more flexible model form that can capture the FCT distortion associated with alternating Fe ↑ and Fe ↓ layers. The third set of bars in Fig. 4(a) show that the spin GAP successfully predicts lattice constants in excellent agreement with DFT, and the experimental values in Table III.

The cumulative distribution functions of the nonzero elastic matrix components, computed using the matscipy PYTHON package [63], of EAM-13, spin GAP, and DFT are shown in Fig. 4(b). As the EAM predicts an FCC structure, it has only three distinct nonzero elastic constants {$C_{11}$, $C_{12}$, $C_{44}$}. The DFT and spin GAP predict an FCT structure, and so have six distinct nonzero elastic constants {$C_{11}$, $C_{33}$, $C_{12}$, $C_{13}$, $C_{44}$, $C_{66}$}. The standard GAP on the other hand has nine distinct nonzero elastic constants as it predicts an FCO structure, the mean values of which are in Table III.

As seen in the first two sets of bars in Fig. 4(b), the spin GAP slightly overestimates $C_{11}$ and $C_{33}$ by ~8% each, but gets their relative magnitude right. This is a significant improvement from EAM-13, which overestimates the DFT $C_{11}$ by 30%. The EAM-13 $C_{11}$ is in fact closer to that of pure AFM Fe reported in Table III, indicating that EAM-13 does not capture the softening of this elastic mode due to alloying elements Cr and Ni, while the spin GAP does capture this softening. Moving to the off-diagonal elastic components, the spin GAP $C_{13}$ agrees very well with DFT, and is an improvement over the EAM-13 that overestimates it by 20%. The spin GAP performs less well with $C_{12}$, overestimating it by 10%, but again this is a significant improvement over EAM-13, which overestimates this component by 58%. Lastly, the spin GAP $C_{44}$ and $C_{66}$ match very well with DFT as seen in the last two sets of bars in Fig. 4(b), while EAM-13 predicts an averaged value between the DFT {$C_{44}$, $C_{66}$} range. Overall, the GAP elastic constants agree well with DFT. The 8–10% overestimation of modes {$C_{11}$, $C_{33}$, $C_{13}$} indicate that there are still errors in getting the behavior of the long vs short axis of the spin-induced tetrahedral deformation. This is likely because our spin model is still quite a simple approximation to the full physics of the alloy—it excludes contributions from Cr and Ni spins and ignores spin magnitudes.

Given that the spin GAP predicts the different elastic constants well, it unsurprisingly also predicts the equation of state reasonably well. Figure 5(a) compares the three spin GAP E-V curves (AFM spin layering along $\hat{X}$, $\hat{Y}$, and $\hat{Z}$ directions) of a 256-atom test configuration with respect to its DFT EV

<footnote><sup>a</sup>Temperature range 750–900 K. Ref. [60].</footnote>

<footnote><sup>b</sup>Temperature range 4–295 K. Ref. [61,62].</footnote>
<table><caption>TABLE III. Lattice parameters and elastic constants of the FeCrNi alloy computed using DFT, GAP, spin GAP, and EAM-13 at 0 K. For comparison, DFT results for AFM Fe at 0 K and experimental values for finite temperatures have also been included. Dashes indicate that the value is equal to the previous value by symmetry.</caption><thead><tr><th></th><th>DFT-Fe</th><th>DFT</th><th>GAP</th><th>Spin GAP</th><th>EAM13</th><th>Expt (&gt; 0 K)</th></tr></thead><tbody><tr><td>Phase</td><td>FCT</td><td>FCT</td><td>FCO</td><td>FCT</td><td>FCC</td><td>~FCC <sup>a</sup></td></tr><tr><td>a (Å)</td><td>3.42</td><td>3.50</td><td>3.49</td><td>3.50</td><td>3.55</td><td>3.52–3.59 <sup>a</sup></td></tr><tr><td>b (Å)</td><td>–</td><td>–</td><td>3.53</td><td>–</td><td>–</td><td>–</td></tr><tr><td>c (Å)</td><td>3.68</td><td>3.58</td><td>3.56</td><td>3.58</td><td>–</td><td>–</td></tr><tr><td>C<sub>11</sub>, C<sub>22</sub>, C<sub>33</sub> (GPa)</td><td>321, –, 254</td><td>259, –, 249</td><td>299, 243, 232</td><td>272, –, 268</td><td>333, –</td><td>204–226<sup>b</sup></td></tr><tr><td>C<sub>12</sub>, C<sub>13</sub>, C<sub>23</sub> (GPa)</td><td>246, 96, –</td><td>164, 124, –</td><td>175, 163, 151</td><td>181, 126, –</td><td>196, –, –</td><td>132–134<sup>b</sup></td></tr><tr><td>C<sub>44</sub>, C<sub>55</sub>, C<sub>66</sub> (GPa)</td><td>166, –, 262</td><td>147, –, 191</td><td>155, 171, 186</td><td>148, –, 188</td><td>161, –, –</td><td>111–122<sup>b</sup></td></tr></tbody></table>

mixed radial and species channels. To check the effects of compression, a compressed version of the three-species GAP was also fit. Figure 3 shows that the compressed GAP has comparable errors to the standard GAP, verifying that there is no significant loss of accuracy when compression is used for this system. The spin GAP, which also uses compressed descriptors, has significantly improved training and test errors over those of the standard GAP. Typical energy errors decrease by 1 meV/atom, force errors by about 50 meV/Å, and stress training errors decrease by about 0.2 GPa.

### III. RESULTS

#### A. Bulk properties

The ground-state bulk properties of the FeCrNi alloy predicted by the GAP and the spin GAP are summarized in Table III, along with the corresponding DFT, EAM, and experimental values for comparison, as well as DFT for pure Fe in an AFM spin state. The DFT results for the ground-state structure of both AFM-Fe and AFM-FeCrNi are FCT, with the longer lattice constant corresponding to the AFM layering direction. This tetragonal distortion is more significant in pure AFM Fe (0.26 Å) compared to the alloy (0.08 Å), as the high concentration of Cr and Ni in the alloy lattice reduce the magnetism-induced geometry effects. Relaxation of EM structures with EAM-13 gives a lattice constant within the DFT and experimental range; however, EAM-13 fails to predict the lattice distortion, leading to a symmetric face centered cubic (FCC) lattice instead, as seen from the first set of bars in Fig. 4(a).

In comparison, the standard GAP predicts that a distorted structure is more energetically favorable than a symmetric FCC structure, as it is trained on relaxed DFT structures that have FCT lattices with zero atomic forces. However, since the standard GAP has no spin degrees of freedom, it does not have the framework to process the spin inputs required to elongate one direction over the other two. Instead, the GAP training correlates the distortion to minor differences in the alloy composition along the different axes in the relaxed training set geometries, and so predicts a face-centred orthorhombic (FCO) relaxed lattice based on composition along the different axes of test structures. As seen from the second set of bars in Fig. 4(a), the three lattice constants of the GAP relaxed FCO structure are equally spaced in a similar range to that

of the two DFT relaxed FCT lattice constants. The spin GAP, on the other hand, has a more flexible model form that can capture the FCT distortion associated with alternating Fe ↑ and Fe ↓ layers. The third set of bars in Fig. 4(a) show that the spin GAP successfully predicts lattice constants in excellent agreement with DFT, and the experimental values in Table III.

The cumulative distribution functions of the nonzero elastic matrix components, computed using the matscipy PYTHON package [63], of EAM-13, spin GAP, and DFT are shown in Fig. 4(b). As the EAM predicts an FCC structure, it has only three distinct nonzero elastic constants {$C_{11}$, $C_{12}$, $C_{44}$}. The DFT and spin GAP predict an FCT structure, and so have six distinct nonzero elastic constants {$C_{11}$, $C_{33}$, $C_{12}$, $C_{13}$, $C_{44}$, $C_{66}$}. The standard GAP on the other hand has nine distinct nonzero elastic constants as it predicts an FCO structure, the mean values of which are in Table III.

As seen in the first two sets of bars in Fig. 4(b), the spin GAP slightly overestimates $C_{11}$ and $C_{33}$ by ~8% each, but gets their relative magnitude right. This is a significant improvement from EAM-13, which overestimates the DFT $C_{11}$ by 30%. The EAM-13 $C_{11}$ is in fact closer to that of pure AFM Fe reported in Table III, indicating that EAM-13 does not capture the softening of this elastic mode due to alloying elements Cr and Ni, while the spin GAP does capture this softening. Moving to the off-diagonal elastic components, the spin GAP $C_{13}$ agrees very well with DFT, and is an improvement over the EAM-13 that overestimates it by 20%. The spin GAP performs less well with $C_{12}$, overestimating it by 10%, but again this is a significant improvement over EAM-13, which overestimates this component by 58%. Lastly, the spin GAP $C_{44}$ and $C_{66}$ match very well with DFT as seen in the last two sets of bars in Fig. 4(b), while EAM-13 predicts an averaged value between the DFT [$C_{44}$, $C_{66}$] range. Overall, the GAP elastic constants agree well with DFT. The 8–10% overestimation of modes {$C_{11}$, $C_{33}$, $C_{13}$} indicate that there are still errors in getting the behavior of the long vs short axis of the spin-induced tetrahedral deformation. This is likely because our spin model is still quite a simple approximation to the full physics of the alloy—it excludes contributions from Cr and Ni spins and ignores spin magnitudes.

Given that the spin GAP predicts the different elastic constants well, it unsurprisingly also predicts the equation of state reasonably well. Figure 5(a) compares the three spin GAP E-V curves (AFM spin layering along $\hat{X}$, $\hat{Y}$, and $\hat{Z}$ directions) of a 256-atom test configuration with respect to its DFT EV

<footnote><sup>a</sup>Temperature range 750–900 K. Ref. [60].</footnote>

<footnote><sup>b</sup>Temperature range 4–295 K. Ref. [61,62].</footnote>
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>FIG. 4. (a) Lattice constants (a), (b), (c) averaged over 400 configurations predicted by EAM, GAP, and spin GAP for the AFM state, with respect to the DFT AFM lattice constants averaged over 111 configurations. (b) The six independent elastic constants of FCT structures predicted by spin GAP (green) and EAM (pink) from the average over 100 configurations each, compared to the DFT reference values computed for one configuration.</figcaption></figure>

curves. We can see than the spin GAP EV curves fall within the range of the DFT EV curves, even if they do not match exactly. The spin GAP also gets the relaxed volume (E-V curve minimum) right, whereas the EAM-13 overestimates the relaxed volume as it is restricted to the FCC lattice. The bulk modulii computed from fitting an equation of state to these spin GAP E-V curves are in very good agreement with DFT, overestimating it by just 4 GPa (~2%), as seen in Fig. 5(b) (green distribution vs black dashed DFT reference), whereas the EAM-13 (pink distribution) overestimates the bulk modulus by 21%.

As the spin GAP is orders of magnitude faster than DFT, it allows for much more sampling of different random configurations, or larger supercells than the 256-atom ones used for DFT. Figure 6 shows the mean and spread (one standard deviation) of the spin GAP relaxed energy distributions at different supercell sizes, computed over 40 different configurations for each size. We see that the spin GAP relaxed energy converges for supercells with more than 200 atoms, and its converged energy is in good agreement with the DFT reference values (black crosses) within the converged size range. For smaller sizes, as is the case when applying DFT to random alloys, finite-size effects are significant leading to unconverged energy predictions. Also, the distributions at smaller sizes are broader, because each configuration's energy is computed over a smaller set of random local environments, leading to a larger variance between different configurations.

For the cell size used for training (256-atoms), the full distribution of spin GAP relaxed energies is in Fig. 10(a), and we can see that its peak agrees very well with the DFT mean. The spin-agnostic GAP on the other hand underestimates the ground-state energy by about 2 meV/atom. We attribute this underestimation to the fact that, while the GAP correlates the tendency to distort to minor changes in composition, it associates FCT directionality to noise as it does not have spin degrees of freedom. The GAP train-

ing regularizes the forces of the low energy DFT training structures, and extrapolates to predicting FCO as a lower-energy structure than FCT. The spin GAP overcomes this model form error by associating the tetragonal distortion with the AFM layering.

## B. Chemical potentials

An essential ingredient for accurately computing energetics of point defects and their clusters is the chemical potential $\mu$ of each element in the alloy. This is the energy required to move an atom of the element from vacuum (infinitely far away from any interactions) to its bulk lattice site in the alloy supercell. The energy $E^{A \to B}$ required to substitute an element A with another element B at a lattice site in an alloy is equal to

$$E^{A \to B} = \mu^{B} - \mu^{A}. \qquad (8)$$

This gives three equations for $A, B \in \{\text{Fe, Cr, Ni}\}$, of which only two are independent. To solve for the three unknown chemical potentials, we require three independent equations, and we can get this third equation from the sum of chemical potentials of the N atoms in a configuration, which gives the bulk energy $E_N$ of that configuration. For $\text{Fe}_7\text{Cr}_2\text{Ni}$, this is expressed by the equation

$$E_N = (0.7\mu^{\text{Fe}} + 0.2\mu^{\text{Cr}} + 0.1\mu^{\text{Ni}})N. \qquad (9)$$

The chemical potentials $\mu^{\text{Fe}}, \mu^{\text{Cr}}$, and $\mu^{\text{Ni}}$ can be determined for any lattice site in a cell by computing the two relevant substitution energies, and then solving the system of equations defined by Eqs. (8) and (9). For the spin GAP, we solve this three-species system of equations twice: first for $\{\text{Fe} \uparrow, \text{Cr}, \text{Ni}\}$, and then for $\{\text{Fe} \downarrow, \text{Cr}, \text{Ni}\}$, and we then take the average of the two resulting chemical potentials per species to get $\mu^{\text{Fe}}, \mu^{\text{Cr}}$, and $\mu^{\text{Ni}}$.
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>FIG. 4. (a) Lattice constants (a), (b), (c) averaged over 400 configurations predicted by EAM, GAP, and spin GAP for the AFM state, with respect to the DFT AFM lattice constants averaged over 111 configurations. (b) The six independent elastic constants of FCT structures predicted by spin GAP (green) and EAM (pink) from the average over 100 configurations each, compared to the DFT reference values computed for one configuration.</figcaption></figure>

curves. We can see than the spin GAP EV curves fall within the range of the DFT EV curves, even if they do not match exactly. The spin GAP also gets the relaxed volume (E-V curve minimum) right, whereas the EAM-13 overestimates the relaxed volume as it is restricted to the FCC lattice. The bulk modulii computed from fitting an equation of state to these spin GAP E-V curves are in very good agreement with DFT, overestimating it by just 4 GPa (~2%), as seen in Fig. 5(b) (green distribution vs black dashed DFT reference), whereas the EAM-13 (pink distribution) overestimates the bulk modulus by 21%.

As the spin GAP is orders of magnitude faster than DFT, it allows for much more sampling of different random configurations, or larger supercells than the 256-atom ones used for DFT. Figure 6 shows the mean and spread (one standard deviation) of the spin GAP relaxed energy distributions at different supercell sizes, computed over 40 different configurations for each size. We see that the spin GAP relaxed energy converges for supercells with more than 200 atoms, and its converged energy is in good agreement with the DFT reference values (black crosses) within the converged size range. For smaller sizes, as is the case when applying DFT to random alloys, finite-size effects are significant leading to unconverged energy predictions. Also, the distributions at smaller sizes are broader, because each configuration's energy is computed over a smaller set of random local environments, leading to a larger variance between different configurations.

For the cell size used for training (256-atoms), the full distribution of spin GAP relaxed energies is in Fig. 10(a), and we can see that its peak agrees very well with the DFT mean. The spin-agnostic GAP on the other hand underestimates the ground-state energy by about 2 meV/atom. We attribute this underestimation to the fact that, while the GAP correlates the tendency to distort to minor changes in composition, it associates FCT directionality to noise as it does not have spin degrees of freedom. The GAP train-

ing regularizes the forces of the low energy DFT training structures, and extrapolates to predicting FCO as a lower-energy structure than FCT. The spin GAP overcomes this model form error by associating the tetragonal distortion with the AFM layering.

## B. Chemical potentials

An essential ingredient for accurately computing energetics of point defects and their clusters is the chemical potential $\mu$ of each element in the alloy. This is the energy required to move an atom of the element from vacuum (infinitely far away from any interactions) to its bulk lattice site in the alloy supercell. The energy $E^{A \to B}$ required to substitute an element A with another element B at a lattice site in an alloy is equal to

$$E^{A \to B} = \mu^{B} - \mu^{A}. \qquad (8)$$

This gives three equations for $A, B \in \{\text{Fe}, \text{Cr}, \text{Ni}\}$, of which only two are independent. To solve for the three unknown chemical potentials, we require three independent equations, and we can get this third equation from the sum of chemical potentials of the N atoms in a configuration, which gives the bulk energy $E_N$ of that configuration. For $\text{Fe}_7\text{Cr}_2\text{Ni}$, this is expressed by the equation

$$E_N = (0.7\mu^{\text{Fe}} + 0.2\mu^{\text{Cr}} + 0.1\mu^{\text{Ni}})N. \qquad (9)$$

The chemical potentials $\mu^{\text{Fe}}, \mu^{\text{Cr}},$ and $\mu^{\text{Ni}}$ can be determined for any lattice site in a cell by computing the two relevant substitution energies, and then solving the system of equations defined by Eqs. (8) and (9). For the spin GAP, we solve this three-species system of equations twice: first for $\{\text{Fe} \uparrow, \text{Cr}, \text{Ni}\}$, and then for $\{\text{Fe} \downarrow, \text{Cr}, \text{Ni}\}$, and we then take the average of the two resulting chemical potentials per species to get $\mu^{\text{Fe}}, \mu^{\text{Cr}},$ and $\mu^{\text{Ni}}$.
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 5. (a) Energy volume (E-V) curves predicted by DFT (crosses) and spin GAP (solid lines) with AFM layering in ˆx, ˆy, and ˆz directions, and the corresponding spin-agnostic EAM-13 (dotted-pink line). (b) Bulk modulus computed from a Birch-Murnaghan equation of state fit to the E-V curves.</figcaption></figure>

This procedure was carried out for randomly chosen lattice sites of 800 different configurations using the spin GAP. The distributions of substitution energies are shown in Fig. 7(a), with reference DFT minimum values $E_{\text{subs}}^{\text{A} \to \text{B}}$ from [5]. We see that the minimum substitution energies from literature (using the same DFT settings as this paper) agree well with the minima of the corresponding spin GAP substitution energy distributions. The resulting chemical potential distributions are shown in the violin plots of Fig. 7(b). DFT chemical potentials averaged over five configurations are also labeled using black crosses. We can see that the spin GAP $\langle \mu^{\text{Fe}} \rangle$ agrees very well with DFT. The spin GAP $\langle \mu^{\text{Ni}} \rangle$ and $\langle \mu^{\text{Cr}} \rangle$ are slightly offset from the corresponding DFT averages by +2% and -2% respectively; this is probably a result of not accounting fully for the magnetic nature of Ni and Cr in this alloy.

<figure><img src="image_10.png" /><figcaption>FIG. 6. Convergence of spin GAP ground-state energy with respect to supercell size. The green dots and bars show the mean and standard deviation respectively, of the spin GAP energy distribution of 40 configurations at each cell size. The spin GAP energy converges to 1 meV/atom for cells with > 200 atoms, and the converged energy agrees well with DFT (black crosses).</figcaption></figure>

### C. Monovacancies

To compute vacancy formation energies, we start with a relaxed bulk configuration of N atoms and energy $E_N$, remove one atom of a given species “el”, and then carry out a fixed-cell relaxation of the atomic positions to get the energy $E_{N-1}^{\text{el}}$ of the vacancy configuration. As the vacancy configuration has one atom less than the corresponding bulk system, the comparison of the two configurations requires the chemical potential $\mu^{\text{el}}$ of the atom removed to form the vacancy, giving the formula

$$E_{\text{vac}}^{\text{el}} = \mu^{\text{el}} + E_{N-1}^{\text{el}} - E_N \quad (10)$$

where quantities $(\mu^{\text{el}} + E_{N-1}^{\text{el}})$ and $E_N$ both have N atoms each. We first compute the intermediate step $(E_{N-1}^{\text{el}} - E_N)$ for 800 different configurations using the spin GAP and the results are shown in Fig. 7(c). In this case, the spin GAP result for Cr agrees well with DFT, whereas those for Fe and Ni are underestimated. These offsets are consistent with the spin GAP underestimating $E^{\text{Fe} \to \text{Ni}}$ slightly, as seen previously in Fig. 7(a).

In fact, the small species-dependent errors in $\mu^{\text{el}}$ and $E_{N-1}^{\text{el}} - E_N$ cancel each other out to give a consistent behavior in vacancy energy $E_{\text{vac}}$, irrespective of which element was removed to form the vacancy. This is evident from the consistent mean and spread of the $E_{\text{vac}}$ energy distributions denoted by the bars in Fig. 7(d). This cancellation of errors makes sense as $\mu^{\text{el}}$ involves introducing an atom into the bulk whereas $E_{N-1}^{\text{el}} - E_N$ involves the opposite procedure of removing an atom from the bulk structure. The resulting prediction of $E_{\text{vac}}^{\text{SpinGAP}} = 1.80$ eV is in good agreement with the DFT value of $E_{\text{vac}}^{\text{DFT}} = 1.96$ eV from the literature [5]. This prediction for $E_{\text{vac}}$ is significantly more accurate than the standard spin-agnostic GAP trained on the same database, which predicts $E_{\text{vac}}^{\text{GAP}} = 1.6$ eV. The spin GAP also performs marginally better than EAM-13, which gives a mean prediction of $E_{\text{vac}} = 1.77$ eV as seen in Fig. 7(d). We attribute the
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 5. (a) Energy volume (E-V) curves predicted by DFT (crosses) and spin GAP (solid lines) with AFM layering in ˆx, ˆy, and ˆz directions, and the corresponding spin-agnostic EAM-13 (dotted-pink line). (b) Bulk modulus computed from a Birch-Murnaghan equation of state fit to the E-V curves.</figcaption></figure>

This procedure was carried out for randomly chosen lattice sites of 800 different configurations using the spin GAP. The distributions of substitution energies are shown in Fig. 7(a), with reference DFT minimum values $E_{\text{subs}}^{\text{A} \to \text{B}}$ from [5]. We see that the minimum substitution energies from literature (using the same DFT settings as this paper) agree well with the minima of the corresponding spin GAP substitution energy distributions. The resulting chemical potential distributions are shown in the violin plots of Fig. 7(b). DFT chemical potentials averaged over five configurations are also labeled using black crosses. We can see that the spin GAP $\langle \mu^{\text{Fe}} \rangle$ agrees very well with DFT. The spin GAP $\langle \mu^{\text{Ni}} \rangle$ and $\langle \mu^{\text{Cr}} \rangle$ are slightly offset from the corresponding DFT averages by +2% and -2% respectively; this is probably a result of not accounting fully for the magnetic nature of Ni and Cr in this alloy.

<figure><img src="image_10.png" /><figcaption>FIG. 6. Convergence of spin GAP ground-state energy with respect to supercell size. The green dots and bars show the mean and standard deviation respectively, of the spin GAP energy distribution of 40 configurations at each cell size. The spin GAP energy converges to 1 meV/atom for cells with > 200 atoms, and the converged energy agrees well with DFT (black crosses).</figcaption></figure>

### C. Monovacancies

To compute vacancy formation energies, we start with a relaxed bulk configuration of N atoms and energy $E_N$, remove one atom of a given species “el”, and then carry out a fixed-cell relaxation of the atomic positions to get the energy $E_{N-1}^{\text{el}}$ of the vacancy configuration. As the vacancy configuration has one atom less than the corresponding bulk system, the comparison of the two configurations requires the chemical potential $\mu^{\text{el}}$ of the atom removed to form the vacancy, giving the formula

$$E_{\text{vac}}^{\text{el}} = \mu^{\text{el}} + E_{N-1}^{\text{el}} - E_N \quad (10)$$

where quantities $(\mu^{\text{el}} + E_{N-1}^{\text{el}})$ and $E_N$ both have N atoms each. We first compute the intermediate step $(E_{N-1}^{\text{el}} - E_N)$ for 800 different configurations using the spin GAP and the results are shown in Fig. 7(c). In this case, the spin GAP result for Cr agrees well with DFT, whereas those for Fe and Ni are underestimated. These offsets are consistent with the spin GAP underestimating $E^{\text{Fe} \to \text{Ni}}$ slightly, as seen previously in Fig. 7(a).

In fact, the small species-dependent errors in $\mu^{\text{el}}$ and $E_{N-1}^{\text{el}} - E_N$ cancel each other out to give a consistent behavior in vacancy energy $E_{\text{vac}}$, irrespective of which element was removed to form the vacancy. This is evident from the consistent mean and spread of the $E_{\text{vac}}$ energy distributions denoted by the bars in Fig. 7(d). This cancellation of errors makes sense as $\mu^{\text{el}}$ involves introducing an atom into the bulk whereas $E_{N-1}^{\text{el}} - E_N$ involves the opposite procedure of removing an atom from the bulk structure. The resulting prediction of $E_{\text{vac}}^{\text{SpinGAP}} = 1.80$ eV is in good agreement with the DFT value of $E_{\text{vac}}^{\text{DFT}} = 1.96$ eV from the literature [5]. This prediction for $E_{\text{vac}}$ is significantly more accurate than the standard spin-agnostic GAP trained on the same database, which predicts $E_{\text{vac}}^{\text{GAP}} = 1.6$ eV. The spin GAP also performs marginally better than EAM-13, which gives a mean prediction of $E_{\text{vac}} = 1.77$ eV as seen in Fig. 7(d). We attribute the
<figure><img src="image_11.png" /><figcaption>FIG. 7. (a) Distribution of substitution energies for elemental swaps predicted by spin GAP (100 configurations each) and corresponding DFT minima (dashed lines) computed from 40 DFT calculations performed by Piochaud *et al.* [5]. Distribution of (b) chemical potentials $\mu$ and (c) energy differences $E_{N-1} - E_N$ predicted by the spin GAP from 800 configurations (256-atom) with respect to DFT means from five configurations (256-atom). (d) Mean and twice the standard deviations of spin GAP and EAM-13 vacancy formation energies $E_{vac}$ of 200 different EM configurations respectively, split by vacancy element type. The DFT mean and range of $E_{vac}$ are taken from [5].</figcaption></figure>

remaining error of –0.2 eV in the spin GAP prediction to the incomplete magnetic model, as we do not account for spin polarization of the electron density of Cr and Ni or for variations in Fe magnetic-moment magnitudes near to defect sites. A DFT study on this alloy finds that the DFT vacancy formation energy drops from 1.98 eV down to 1.80 eV when they switch from imposing paramagnetic spins to ferrimagnetic spins on the same chemical composition, showing that changes in the magnetic state alone can indeed lead to an order of 0.2 eV changes in $E_{vac}$ [12], further suggesting that the spin GAP’s underestimation of $E_{vac}$ is likely due to limitations of our spin model.

On the other hand, a recently reported fully noncollinear magnetic atomic cluster expansion (ACE) model for iron also makes similar magnitudes errors in vacancy formation

energies (errors of 0.4 eV and 0.2 eV when trained without and with extended defect sub-databases respectively), despite being a much more rigorous magnetic model [64]. In their case, the errors probably come from challenges in sampling their multiplicatively larger spin space (as spin magnitudes and angles are included). Whereas in our case, our simpler spin model comes with a relatively constrained spin space that is much easier to sample well. It is a promising result that we get comparable performance, suggesting that our spin model captures the main magnetic contributions despite being relatively simple. It would be interesting to try relaxing a few of the assumptions of our simple spin model while still keeping the spin space manageable, but that is beyond the scope of this paper and will be investigated in future works.
<figure><img src="image_11.png" /><figcaption>FIG. 7. (a) Distribution of substitution energies for elemental swaps predicted by spin GAP (100 configurations each) and corresponding DFT minima (dashed lines) computed from 40 DFT calculations performed by Piochaud *et al.* [5]. Distribution of (b) chemical potentials $\mu$ and (c) energy differences $E_{N-1} - E_N$ predicted by the spin GAP from 800 configurations (256-atom) with respect to DFT means from five configurations (256-atom). (d) Mean and twice the standard deviations of spin GAP and EAM-13 vacancy formation energies $E_{vac}$ of 200 different EM configurations respectively, split by vacancy element type. The DFT mean and range of $E_{vac}$ are taken from [5].</figcaption></figure>

remaining error of –0.2 eV in the spin GAP prediction to the incomplete magnetic model, as we do not account for spin polarization of the electron density of Cr and Ni or for variations in Fe magnetic-moment magnitudes near to defect sites. A DFT study on this alloy finds that the DFT vacancy formation energy drops from 1.98 eV down to 1.80 eV when they switch from imposing paramagnetic spins to ferrimagnetic spins on the same chemical composition, showing that changes in the magnetic state alone can indeed lead to an order of 0.2 eV changes in $E_{vac}$ [12], further suggesting that the spin GAP’s underestimation of $E_{vac}$ is likely due to limitations of our spin model.

On the other hand, a recently reported fully noncollinear magnetic atomic cluster expansion (ACE) model for iron also makes similar magnitudes errors in vacancy formation

energies (errors of 0.4 eV and 0.2 eV when trained without and with extended defect sub-databases respectively), despite being a much more rigorous magnetic model [64]. In their case, the errors probably come from challenges in sampling their multiplicatively larger spin space (as spin magnitudes and angles are included). Whereas in our case, our simpler spin model comes with a relatively constrained spin space that is much easier to sample well. It is a promising result that we get comparable performance, suggesting that our spin model captures the main magnetic contributions despite being relatively simple. It would be interesting to try relaxing a few of the assumptions of our simple spin model while still keeping the spin space manageable, but that is beyond the scope of this paper and will be investigated in future works.
<table><caption>TABLE IV. Linear correlation of vacancy formation energy with respect to number of first nearest neighbours (#1nn) of type Cr and Ni respectively, from the literature [5,10] (columns 2 and 3) and our paper (column 5), and for ($E_{N-1}^{Fe} - E_N$) from our paper (column 4).</caption><thead><tr><th>#1nn</th><th>E<sub>vac</sub> ([5])</th><th>E<sub>vac</sub> ([10])</th><th>(E<sub>N-1</sub> - E<sub>N</sub>)</th><th>E<sub>vac</sub></th></tr></thead><tbody><tr><td>Cr</td><td>0.05</td><td>0.02</td><td>0.02</td><td>0.01</td></tr><tr><td>Ni</td><td>-0.04</td><td>-0.04</td><td>-0.01</td><td>-0.004</td></tr></tbody></table>

### D. Atomic Short Range Order

DFT studies [5,10] conclude that vacancies in the Fe<sub>7</sub>Cr<sub>2</sub>Ni are more stable in Ni-rich environments and less stable in Cr-rich environments. The gradients reported in the literature for the linear correlation between $E_{vac}$ and number of Cr and Ni neighbours respectively have been included in Table IV for reference. The two studies agree well on the correlation with Ni neighbors, whereas Manzoor *et al.* find a weaker correlation with respect to Cr neighbors compared to Piochaud *et al.*, which they attribute to the fact that they have four times the amount of sampling [10]. In both cases, however, they use constant precomputed values for the chemical potentials, so trends in their vacancy formation energy Eq. (10) are effectively proportional to trends in just $E_{N-1}^{el} - E_N$.

The computational expense of calculating chemical potentials, i.e., two fixed-cell relaxations of element substitutions Eq. (8) for each vacancy, makes it impractical to compute it for every vacancy using DFT, but this is now possible to run cheaply with the spin GAP. Figure 8 from left to right shows the trends in energy difference $E_{N-1}^{el} - E_N$, the corresponding chemical potentials, and the vacancy formation energies with respect to number of Cr neighbors as predicted by the spin GAP. These results show that the chemical potential is also

correlated to the number of Cr neighbors, with the opposite slope to the energy difference $E_{N-1}^{el} - E_N$. Computing vacancy formation energy involves adding up these two quantities [see Eq. (10)] and so their opposing slopes offset each other, reducing the overall magnitude of positive correlation between $E_{vac}$ and number of Cr neighbors. Therefore, while we get consistent gradient of 0.02 eV/Cr for our $E_{N-1}^{el} - E_N$ and the $E_{vac}$ from the literature (see values in Table IV), the inclusion of variations in chemical potentials decreases the overall correlation. This makes sense as the chemical potential accounts for majority of the heterogeneity in the vacancies neighborhood. The remnant positive correlation is likely because of the volume of the vacancy, which has previously been seen to affect stability [10]. As Cr atoms are the smallest species in the alloy, they relax less into the vacancy, leading to larger vacancies that are less stable and so have a higher formation energy.

The inverse trends (but to a weaker extent) of the three subplots of Fig. 8 are seen with respect to number of Ni neighbors, the gradients of which are reported in Table IV. The spin GAP predicts a smaller negative correlation of -0.01 eV/Ni for $E_{N-1}^{el} - E_N$ compared to the -0.04 eV/Ni values in the literature. Similar to the discrepancy between the two DFT studies for Cr trends, our lower correlation could be due to more variety in sampling. Piochaud *et al.* generate all their vacancy structure by removing different atoms from one 256-atom SQS structure, and Manzoor *et al.* do the same over two 256-atom SQS structures, leading to a large overlap in the neighborhood of the vacancies they consider. On the other hand, our spin GAP results were computed for vacancies introduced in 200 different starting structures from effective medium theory, and so sample a larger variety in long range vacancy environments. Despite the weaker correlations, we do still get a small negative relation of -0.004 eV/Ni for $E_{vac}$ and number of Ni neighbors. Similar to the previous reasoning on

<figure><img src="image_12.png" /><figcaption>FIG. 8. Linear correlation between (a) energy difference of vacancy vs bulk configurations $E_{N-1}^{Fe} - E_N$ for vacancies generated by removing an iron atom each from the 200 different bulk configurations, (b) corresponding iron chemical potentials $\mu_{Fe}$ for those vacancy sites, and (c) corresponding vacancy formation energy $E_{vac}$ computed by adding the quantities of (a) and (b), all with respect to number of Cr atoms amongst the 12 nearest neighbors of the vacancy.</figcaption></figure>
<table><caption>TABLE IV. Linear correlation of vacancy formation energy with respect to number of first nearest neighbours (#1nn) of type Cr and Ni respectively, from the literature [5,10] (columns 2 and 3) and our paper (column 5), and for ($E_{N-1}^{\text{Fe}} - E_N$) from our paper (column 4).</caption><thead><tr><th>#1nn</th><th><i>E</i><sub>vac</sub> ([5])</th><th><i>E</i><sub>vac</sub> ([10])</th><th>(<i>E</i><sub>N-1</sub> - <i>E</i><sub>N</sub>)</th><th><i>E</i><sub>vac</sub></th></tr></thead><tbody><tr><td>Cr</td><td>0.05</td><td>0.02</td><td>0.02</td><td>0.01</td></tr><tr><td>Ni</td><td>-0.04</td><td>-0.04</td><td>-0.01</td><td>-0.004</td></tr></tbody></table>

### D. Atomic Short Range Order

DFT studies [5,10] conclude that vacancies in the Fe<sub>7</sub>Cr<sub>2</sub>Ni are more stable in Ni-rich environments and less stable in Cr-rich environments. The gradients reported in the literature for the linear correlation between $E_{vac}$ and number of Cr and Ni neighbours respectively have been included in Table IV for reference. The two studies agree well on the correlation with Ni neighbors, whereas Manzoor *et al.* find a weaker correlation with respect to Cr neighbors compared to Piochaud *et al.*, which they attribute to the fact that they have four times the amount of sampling [10]. In both cases, however, they use constant precomputed values for the chemical potentials, so trends in their vacancy formation energy Eq. (10) are effectively proportional to trends in just $E_{N-1}^{\text{el}} - E_N$.

The computational expense of calculating chemical potentials, i.e., two fixed-cell relaxations of element substitutions Eq. (8) for each vacancy, makes it impractical to compute it for every vacancy using DFT, but this is now possible to run cheaply with the spin GAP. Figure 8 from left to right shows the trends in energy difference $E_{N-1}^{\text{el}} - E_N$, the corresponding chemical potentials, and the vacancy formation energies with respect to number of Cr neighbors as predicted by the spin GAP. These results show that the chemical potential is also

correlated to the number of Cr neighbors, with the opposite slope to the energy difference $E_{N-1}^{\text{Fe}} - E_N$. Computing vacancy formation energy involves adding up these two quantities [see Eq. (10)] and so their opposing slopes offset each other, reducing the overall magnitude of positive correlation between $E_{vac}$ and number of Cr neighbors. Therefore, while we get consistent gradient of 0.02 eV/Cr for our $E_{N-1}^{\text{Fe}} - E_N$ and the $E_{vac}$ from the literature (see values in Table IV), the inclusion of variations in chemical potentials decreases the overall correlation. This makes sense as the chemical potential accounts for majority of the heterogeneity in the vacancies neighborhood. The remnant positive correlation is likely because of the volume of the vacancy, which has previously been seen to affect stability [10]. As Cr atoms are the smallest species in the alloy, they relax less into the vacancy, leading to larger vacancies that are less stable and so have a higher formation energy.

The inverse trends (but to a weaker extent) of the three subplots of Fig. 8 are seen with respect to number of Ni neighbors, the gradients of which are reported in Table IV. The spin GAP predicts a smaller negative correlation of -0.01 eV/Ni for $E_{N-1}^{\text{Fe}} - E_N$ compared to the -0.04 eV/Ni values in the literature. Similar to the discrepancy between the two DFT studies for Cr trends, our lower correlation could be due to more variety in sampling. Piochaud *et al.* generate all their vacancy structure by removing different atoms from one 256-atom SQS structure, and Manzoor *et al.* do the same over two 256-atom SQS structures, leading to a large overlap in the neighborhood of the vacancies they consider. On the other hand, our spin GAP results were computed for vacancies introduced in 200 different starting structures from effective medium theory, and so sample a larger variety in long range vacancy environments. Despite the weaker correlations, we do still get a small negative relation of -0.004 eV/Ni for $E_{vac}$ and number of Ni neighbors. Similar to the previous reasoning on

<figure><img src="image_12.png" /><figcaption>FIG. 8. Linear correlation between (a) energy difference of vacancy vs bulk configurations $E_{N-1}^{\text{Fe}} - E_N$ for vacancies generated by removing an iron atom each from the 200 different bulk configurations, (b) corresponding iron chemical potentials $\mu_{\text{Fe}}$ for those vacancy sites, and (c) corresponding vacancy formation energy $E_{vac}$ computed by adding the quantities of (a) and (b), all with respect to number of Cr atoms amongst the 12 nearest neighbors of the vacancy.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>FIG. 9. Variation in ground-state energy per atom of bulk alloy with variation of Cr → Fe concentrations (blue) and Ni → Fe concentrations (orange), where 0% refers to the target 70:20:10 ratio of Fe:Cr:Ni. We see good agreement with DFT values (black triangles) for this range of ±4% variation in composition.</figcaption></figure>

stability of vacancy volumes, Ni atoms are the largest species in the alloy, so having more Ni neighbors leads to smaller vacancy volumes, which are marginally more stable.

## E. Alloy Composition

In the 304 and 316 austenitic stainless steel grades, the typical Cr and Ni concentrations can vary up to ±4% from the model 70:20:10 ratio of Fe, Cr, and Ni in our training set. It would hence be useful to check if our spin GAP can extrapolate to small variations in the alloying ratios. Figure 9 shows the results for bulk relaxation for configurations within ±4% variations in Cr (blue) and Ni (orange) concentrations. We see very good agreement with the DFT reference values in this range, validating that the spin GAP can indeed be used to extrapolate to compositions about the training ratio. This is because the GAP framework breaks down each configuration into local atomic environments (LAEs) as in Eq. (2), and the training set is likely to have a few Cr or Ni rich/poor LAEs since we sampled a large number of EM configurations, and these LAEs inform the potential for variations in alloy composition about the target ratio.

The gradients of the linear trends in Fig. 9 are 0.0132 eV/%Cr and 0.0280 eV/%Ni. These convert to 1.32 eV and 2.80 eV for single atom Cr → Fe and Ni → Fe swaps respectively, which are consistent with the respective mean substitution energies of Fig. 7(a). This linear trend indicates that for small variations in alloying concentrations, the substitution energies do not change much. The chemical potential, however, would vary, since their computation involves the bulk energy as well via Eq. (9), and these bulk energies vary linearly as seen in Fig. 9. Since we have accurate bulk and substitution energies for the ±4% concentration range, we can expect all derivative properties computed in the previous results sections to follow through reliably.

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>FIG. 10. (a) DFT mean values (dashed lines) and spin GAP distributions for ground-state energies for AFM (green) vs PM states using uniform random 50% up-down iron spins (unshaded orange) and after Monte Carlo optimisation with spin flip moves (shaded orange), comprising 400, 50, and 50 configurations respectively. (b) Correlation plot of spin GAP relaxed energies for 13 PM reference DFT configurations with respect to their DFT relaxed energies.</figcaption></figure>

## F. Paramagnetism

Although the spin GAP was trained only on AFM DFT data, in this section we assess how well it extrapolates to the paramagnetic (PM) state. Comparing the magnetic ground-state energies, Fig. 10(a) shows that the spin GAP correctly predicts the PM ground state of this alloy to be of higher energy, and hence less stable at 0 K, than the AFM ground state. The unshaded orange distribution comprises energies of the starting configurations initialized with perfectly random 50% up and down spin iron atoms, whereas the shaded orange distribution comprises their equilibrated energies after using Metropolis Monte Carlo with spin flip moves to optimize their paramagnetic spin states. We see that the equilibrated PM energy distribution agrees very well with the mean PM energy from DFT marked by the orange-dashed line in Fig. 10(a),
<figure><img src="image_13.png" /><figcaption>FIG. 9. Variation in ground-state energy per atom of bulk alloy with variation of Cr → Fe concentrations (blue) and Ni → Fe concentrations (orange), where 0% refers to the target 70:20:10 ratio of Fe:Cr:Ni. We see good agreement with DFT values (black triangles) for this range of ±4% variation in composition.</figcaption></figure>

stability of vacancy volumes, Ni atoms are the largest species in the alloy, so having more Ni neighbors leads to smaller vacancy volumes, which are marginally more stable.

### E. Alloy Composition

In the 304 and 316 austenitic stainless steel grades, the typical Cr and Ni concentrations can vary up to ±4% from the model 70:20:10 ratio of Fe, Cr, and Ni in our training set. It would hence be useful to check if our spin GAP can extrapolate to small variations in the alloying ratios. Figure 9 shows the results for bulk relaxation for configurations within ±4% variations in Cr (blue) and Ni (orange) concentrations. We see very good agreement with the DFT reference values in this range, validating that the spin GAP can indeed be used to extrapolate to compositions about the training ratio. This is because the GAP framework breaks down each configuration into local atomic environments (LAEs) as in Eq. (2), and the training set is likely to have a few Cr or Ni rich/poor LAEs since we sampled a large number of EM configurations, and these LAEs inform the potential for variations in alloy composition about the target ratio.

The gradients of the linear trends in Fig. 9 are 0.0132 eV/%Cr and 0.0280 eV/%Ni. These convert to 1.32 eV and 2.80 eV for single atom Cr → Fe and Ni → Fe swaps respectively, which are consistent with the respective mean substitution energies of Fig. 7(a). This linear trend indicates that for small variations in alloying concentrations, the substitution energies do not change much. The chemical potential, however, would vary, since their computation involves the bulk energy as well via Eq. (9), and these bulk energies vary linearly as seen in Fig. 9. Since we have accurate bulk and substitution energies for the ±4% concentration range, we can expect all derivative properties computed in the previous results sections to follow through reliably.

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>FIG. 10. (a) DFT mean values (dashed lines) and spin GAP distributions for ground-state energies for AFM (green) vs PM states using uniform random 50% up-down iron spins (unshaded orange) and after Monte Carlo optimisation with spin flip moves (shaded orange), comprising 400, 50, and 50 configurations respectively. (b) Correlation plot of spin GAP relaxed energies for 13 PM reference DFT configurations with respect to their DFT relaxed energies.</figcaption></figure>

### F. Paramagnetism

Although the spin GAP was trained only on AFM DFT data, in this section we assess how well it extrapolates to the paramagnetic (PM) state. Comparing the magnetic ground-state energies, Fig. 10(a) shows that the spin GAP correctly predicts the PM ground state of this alloy to be of higher energy, and hence less stable at 0 K, than the AFM ground state. The unshaded orange distribution comprises energies of the starting configurations initialized with perfectly random 50% up and down spin iron atoms, whereas the shaded orange distribution comprises their equilibrated energies after using Metropolis Monte Carlo with spin flip moves to optimize their paramagnetic spin states. We see that the equilibrated PM energy distribution agrees very well with the mean PM energy from DFT marked by the orange-dashed line in Fig. 10(a),
<figure><img src="image_16.png" /><figcaption>FIG. 11. Paramagnetic lattice constants (orange bars) and elastic constants (green bars) of the spin GAP averaged over 50 configurations, compared to their respective DFT reference values (black bars). The DFT reference for lattice constants are an average over our 13 DFT PM test configurations, while those for the elastic constants are from the study by Antillon *et al.* [12].</figcaption></figure>

showing that the spin GAP can be used to generate reliable PM spin state configurations despite not being specifically trained for it. Figure 10(b) shows the direct correlation plot for energies of the 13 paramagnetic DFT test structures, where spin GAP was used to evaluate these test structure with the same spin state as that from DFT. We see that the spin GAP reproduces the DFT energies of these paramagnetic configurations very well, as all the points lie along the diagonal with the maximum deviation being 6 meV/atom.

For the PM ground-state structure, DFT predicts a FCC structure with an average lattice constant of 3.53 Å. The spin GAP accurately predicts the paramagnetic FCC structure, with lattice constants that agree remarkably well with DFT, as seen from the orange bars of Fig. 11. The lattice constants of all samples were sorted in ascending order and assigned to *a*, *b*, and *c*, respectively, and Fig. 11 shows their mean values. The small ~0.005 Å difference between the three PM lattice constants is probably due to minor variations in ratio of the three species along the three axes, and this too is consistent between the spin GAP and DFT. As seen from the green bars in Fig. 11, the spin GAP agrees reasonably well with the paramagnetic elastic constants as well, agreeing within 5–8% to their respective DFT reference values from the study by Antillon *et al.* [12]. The fact that the spin GAP reproduces these properties is a promising result, as the paramagnetic state is relevant to modeling the temperature range of most applications of Fe<sub>7</sub>Cr<sub>2</sub>Ni based austenitic steels.

## IV. CONCLUSIONS

In summary, we have developed two machine learning interatomic potentials—a standard GAP and a collinear spin GAP—for a model Fe<sub>7</sub>Cr<sub>2</sub>Ni alloy representative of austenitic steels. They are each trained on the same DFT database containing 159 k atomic environments, and validated against as independent test set. We demonstrate the shortcomings

of the standard GAP in not being able to capture the ground-state tetragonal distortion due to AFM layering observed in this alloy, as the standard GAP formalism cannot account for magnetic spins. We then proposed an extended model incorporating spin to correct for this issue, and verified that the spin GAP predicts the ground-state structure, energy, elastic properties, and vacancies in very good agreement with DFT.

Atomistic modeling using this spin GAP could be used to compute high accuracy inputs to larger scale microstructural models. Our spin GAP is seen to extrapolate well to small variations in alloy composition making is useful for atomistic studies of more variety in austenitic steel grades. The spin model is also seen to describe the paramagnetic state well, which is relevant to modeling the alloy at operating temperatures of the steel components in nuclear and industrial machinery. Beyond the properties validated in this paper, it could also be extended via iterative re-training to model more phenomena of interest such as diffusion, phase transformations or responses to radiation. For instance, by adding and testing the description of grain boundary structures, the spin GAP could be used to study segregation of Cr at grain boundaries under radiation. Another avenue for future work would be to extend this spin GAP to study hydrogen embrittlement in austenitic steel, relevant for applications in future fuel pipelines. This spin GAP is a good starting point for efforts at modeling austenitic steel with close to DFT accuracy.

The DFT training data, interatomic potential and property predictions supporting this work are freely available from Zenodo [65].

<table><caption>TABLE V. Fitted pairwise atom-atom interchange parameters obtained from the effective medium theory for the composition Fe<sub>7</sub>Cr<sub>2</sub>Ni, which are used in lattice-based Monte Carlo simulations to sample configuration space. All energies are in units of meV.</caption><thead><tr><th></th><th>V<sup>(1)</sup><sub>αα'</sub></th><th>Fe</th><th>Cr</th><th>Ni</th></tr></thead><tbody><tr><th scope="row">Fe</th><td></td><td>0.85</td><td>-5.28</td><td>4.58</td></tr><tr><th scope="row">Cr</th><td></td><td>-5.28</td><td>23.44</td><td>-9.95</td></tr><tr><th scope="row">Ni</th><td></td><td>4.58</td><td>-9.95</td><td>-12.18</td></tr></tbody><thead><tr><th></th><th>V<sup>(2)</sup><sub>αα'</sub></th><th>Fe</th><th>Cr</th><th>Ni</th></tr></thead><tbody><tr><th scope="row">Fe</th><td></td><td>0.12</td><td>-0.28</td><td>-0.29</td></tr><tr><th scope="row">Cr</th><td></td><td>-0.28</td><td>2.80</td><td>-3.64</td></tr><tr><th scope="row">Ni</th><td></td><td>-0.29</td><td>-3.64</td><td>9.33</td></tr></tbody><thead><tr><th></th><th>V<sup>(3)</sup><sub>αα'</sub></th><th>Fe</th><th>Cr</th><th>Ni</th></tr></thead><tbody><tr><th scope="row">Fe</th><td></td><td>-0.03</td><td>-0.16</td><td>0.53</td></tr><tr><th scope="row">Cr</th><td></td><td>-0.16</td><td>1.41</td><td>-1.70</td></tr><tr><th scope="row">Ni</th><td></td><td>0.53</td><td>-1.70</td><td>-0.30</td></tr></tbody><thead><tr><th></th><th>V<sup>(4)</sup><sub>αα'</sub></th><th>Fe</th><th>Cr</th><th>Ni</th></tr></thead><tbody><tr><th scope="row">Fe</th><td></td><td>-0.06</td><td>0.51</td><td>-0.59</td></tr><tr><th scope="row">Cr</th><td></td><td>0.51</td><td>-1.03</td><td>-1.50</td></tr><tr><th scope="row">Ni</th><td></td><td>-0.59</td><td>-1.50</td><td>7.15</td></tr></tbody></table>
<figure><img src="image_16.png" /><figcaption>FIG. 11. Paramagnetic lattice constants (orange bars) and elastic constants (green bars) of the spin GAP averaged over 50 configurations, compared to their respective DFT reference values (black bars). The DFT reference for lattice constants are an average over our 13 DFT PM test configurations, while those for the elastic constants are from the study by Antillon *et al.* [12].</figcaption></figure>

showing that the spin GAP can be used to generate reliable PM spin state configurations despite not being specifically trained for it. Figure 10(b) shows the direct correlation plot for energies of the 13 paramagnetic DFT test structures, where spin GAP was used to evaluate these test structure with the same spin state as that from DFT. We see that the spin GAP reproduces the DFT energies of these paramagnetic configurations very well, as all the points lie along the diagonal with the maximum deviation being 6 meV/atom.

For the PM ground-state structure, DFT predicts a FCC structure with an average lattice constant of 3.53 Å. The spin GAP accurately predicts the paramagnetic FCC structure, with lattice constants that agree remarkably well with DFT, as seen from the orange bars of Fig. 11. The lattice constants of all samples were sorted in ascending order and assigned to *a*, *b*, and *c*, respectively, and Fig. 11 shows their mean values. The small ~0.005 Å difference between the three PM lattice constants is probably due to minor variations in ratio of the three species along the three axes, and this too is consistent between the spin GAP and DFT. As seen from the green bars in Fig. 11, the spin GAP agrees reasonably well with the paramagnetic elastic constants as well, agreeing within 5–8% to their respective DFT reference values from the study by Antillon *et al.* [12]. The fact that the spin GAP reproduces these properties is a promising result, as the paramagnetic state is relevant to modeling the temperature range of most applications of Fe<sub>7</sub>Cr<sub>2</sub>Ni based austenitic steels.

## IV. CONCLUSIONS

In summary, we have developed two machine learning interatomic potentials—a standard GAP and a collinear spin GAP—for a model Fe<sub>7</sub>Cr<sub>2</sub>Ni alloy representative of austenitic steels. They are each trained on the same DFT database containing 159 k atomic environments, and validated against as independent test set. We demonstrate the shortcomings

of the standard GAP in not being able to capture the ground-state tetragonal distortion due to AFM layering observed in this alloy, as the standard GAP formalism cannot account for magnetic spins. We then proposed an extended model incorporating spin to correct for this issue, and verified that the spin GAP predicts the ground-state structure, energy, elastic properties, and vacancies in very good agreement with DFT.

Atomistic modeling using this spin GAP could be used to compute high accuracy inputs to larger scale microstructural models. Our spin GAP is seen to extrapolate well to small variations in alloy composition making is useful for atomistic studies of more variety in austenitic steel grades. The spin model is also seen to describe the paramagnetic state well, which is relevant to modeling the alloy at operating temperatures of the steel components in nuclear and industrial machinery. Beyond the properties validated in this paper, it could also be extended via iterative re-training to model more phenomena of interest such as diffusion, phase transformations or responses to radiation. For instance, by adding and testing the description of grain boundary structures, the spin GAP could be used to study segregation of Cr at grain boundaries under radiation. Another avenue for future work would be to extend this spin GAP to study hydrogen embrittlement in austenitic steel, relevant for applications in future fuel pipelines. This spin GAP is a good starting point for efforts at modeling austenitic steel with close to DFT accuracy.

The DFT training data, interatomic potential and property predictions supporting this work are freely available from Zenodo [65].

<table><caption>TABLE V. Fitted pairwise atom-atom interchange parameters obtained from the effective medium theory for the composition Fe<sub>7</sub>Cr<sub>2</sub>Ni, which are used in lattice-based Monte Carlo simulations to sample configuration space. All energies are in units of meV.</caption><thead><tr><th></th><th>V<sup>(1)</sup><sub>αα'</sub></th><th>Fe</th><th>Cr</th><th>Ni</th></tr></thead><tbody><tr><th scope="row">Fe</th><td></td><td>0.85</td><td>-5.28</td><td>4.58</td></tr><tr><th scope="row">Cr</th><td></td><td>-5.28</td><td>23.44</td><td>-9.95</td></tr><tr><th scope="row">Ni</th><td></td><td>4.58</td><td>-9.95</td><td>-12.18</td></tr></tbody><thead><tr><th></th><th>V<sup>(2)</sup><sub>αα'</sub></th><th>Fe</th><th>Cr</th><th>Ni</th></tr></thead><tbody><tr><th scope="row">Fe</th><td></td><td>0.12</td><td>-0.28</td><td>-0.29</td></tr><tr><th scope="row">Cr</th><td></td><td>-0.28</td><td>2.80</td><td>-3.64</td></tr><tr><th scope="row">Ni</th><td></td><td>-0.29</td><td>-3.64</td><td>9.33</td></tr></tbody><thead><tr><th></th><th>V<sup>(3)</sup><sub>αα'</sub></th><th>Fe</th><th>Cr</th><th>Ni</th></tr></thead><tbody><tr><th scope="row">Fe</th><td></td><td>-0.03</td><td>-0.16</td><td>0.53</td></tr><tr><th scope="row">Cr</th><td></td><td>-0.16</td><td>1.41</td><td>-1.70</td></tr><tr><th scope="row">Ni</th><td></td><td>0.53</td><td>-1.70</td><td>-0.30</td></tr></tbody><thead><tr><th></th><th>V<sup>(4)</sup><sub>αα'</sub></th><th>Fe</th><th>Cr</th><th>Ni</th></tr></thead><tbody><tr><th scope="row">Fe</th><td></td><td>-0.06</td><td>0.51</td><td>-0.59</td></tr><tr><th scope="row">Cr</th><td></td><td>0.51</td><td>-1.03</td><td>-1.50</td></tr><tr><th scope="row">Ni</th><td></td><td>-0.59</td><td>-1.50</td><td>7.15</td></tr></tbody></table>
<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /><figcaption>FIG. 12. Cumulative distributions of (a) training and (b) test errors in forces, and in (c) stresses, for the standard SOAP GAP (blue), GAP with compressed SOAP descriptors (orange), and spin GAP (green).</figcaption></figure>

## ACKNOWLEDGMENTS

L.S. and C.D.W. are supported by studentships within the Engineering and Physical Sciences Research Council

supported Centre for Doctoral Training in Modelling of Heterogeneous Systems, Grant No. EP/S022848/1. J.R.K. acknowledges funding from the Leverhulme Trust under Grant No. RPG-2017-191. A.P.B. acknowledges support from the CASTEP-USER project, funded by the Engineering and Physical Sciences Research Council under the Grant Agreement EP/W030438/1. L.S., C.S.B., C.D., and J.R.K. acknowledge funding from the ENTENTE consortium funded by the European Commission under Grant Agreement No. 900018. L.S., C.B., C.D., and J.R.K. acknowledge PRACE support from IRENE at TGCC under Grant No. 2020235585. L.S. and J.R.K. acknowledge usage of the ARCHER2 facility for which access was obtained via the UKCP consortium and funded by EPSRC Grants No. EP/P022561/1 and No. EP/X035891/1. Calculations were performed using the Sulis Tier 2 HPC platform hosted by the Scientific Computing Research Technology Platform at the University of Warwick. Sulis is funded by EPSRC Grant EP/T022108/1 and the HPC Midlands+ consortium. Further computing facilities were provided by the Scientific Computing Research Technology Platform of the University of Warwick.

J.R.K. designed the research, with input from all authors, and supervised L.S. (with support from A.P.B.). With support from J.B.S., C.D.W. performed the KKR-CPA calculations and Monte Carlo simulations to generate the EM supercells. L.S. assembled the database of training configurations, carried out the majority of the DFT calculations, fitted the GAP model (with support from A.P.B. and J.R.K.) and carried out the validation tests. C.D. carried out further DFT calculations, and together with C.S.B. supported on alloy and defect modeling. L.S. drafted the manuscript with support from J.R.K. All authors revised the paper and approved its final version.

## APPENDIX A: EFFECTIVE MEDIUM THEORY PARAMETERS

Tabulated in Table V are the pairwise atom-atom interchange parameters of Eq. (1) obtained for the Fe<sub>7</sub>Cr<sub>2</sub>Ni alloy considered in this paper. The lattice-based model assumes interactions are isotropic and homogeneous, and we write $V_{\alpha\alpha'}^{(n)}$ to denote the energy associated with an atom of species $\alpha$ interacting with an atom of species $\alpha'$ at the $n$th nearest-neighbor distance. We find that it is sufficient to fit interactions to the first four-neighbor distances to accurately capture the data produced using the $S^{(2)}$ theory.

## APPENDIX B: POTENTIAL TRAINING & TESTING

Figure 12 shows the training and testing errors for forces and stresses for the GAP, compressed GAP, and spin GAP. In all cases, the spin GAP is seen to perform the best (lower errors). The errors at which these trends plateau were discussed previously in Sec. II E. The composition of the database is reported in Table I.
<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /><figcaption>FIG. 12. Cumulative distributions of (a) training and (b) test errors in forces, and in (c) stresses, for the standard SOAP GAP (blue), GAP with compressed SOAP descriptors (orange), and spin GAP (green).</figcaption></figure>

## ACKNOWLEDGMENTS

L.S. and C.D.W. are supported by studentships within the Engineering and Physical Sciences Research Council

supported Centre for Doctoral Training in Modelling of Heterogeneous Systems, Grant No. EP/S022848/1. J.R.K. acknowledges funding from the Leverhulme Trust under Grant No. RPG-2017-191. A.P.B. acknowledges support from the CASTEP-USER project, funded by the Engineering and Physical Sciences Research Council under the Grant Agreement EP/W030438/1. L.S., C.S.B., C.D., and J.R.K. acknowledge funding from the ENTENTE consortium funded by the European Commission under Grant Agreement No. 900018. L.S., C.B., C.D., and J.R.K. acknowledge PRACE support from IRENE at TGCC under Grant No. 2020235585. L.S. and J.R.K. acknowledge usage of the ARCHER2 facility for which access was obtained via the UKCP consortium and funded by EPSRC Grants No. EP/P022561/1 and No. EP/X035891/1. Calculations were performed using the Sulis Tier 2 HPC platform hosted by the Scientific Computing Research Technology Platform at the University of Warwick. Sulis is funded by EPSRC Grant EP/T022108/1 and the HPC Midlands+ consortium. Further computing facilities were provided by the Scientific Computing Research Technology Platform of the University of Warwick.

J.R.K. designed the research, with input from all authors, and supervised L.S. (with support from A.P.B.). With support from J.B.S., C.D.W. performed the KKR-CPA calculations and Monte Carlo simulations to generate the EM supercells. L.S. assembled the database of training configurations, carried out the majority of the DFT calculations, fitted the GAP model (with support from A.P.B. and J.R.K.) and carried out the validation tests. C.D. carried out further DFT calculations, and together with C.S.B. supported on alloy and defect modeling. L.S. drafted the manuscript with support from J.R.K. All authors revised the paper and approved its final version.

## APPENDIX A: EFFECTIVE MEDIUM THEORY PARAMETERS

Tabulated in Table V are the pairwise atom-atom interchange parameters of Eq. (1) obtained for the Fe<sub>7</sub>Cr<sub>2</sub>Ni alloy considered in this paper. The lattice-based model assumes interactions are isotropic and homogeneous, and we write $V_{\alpha\alpha'}^{(n)}$ to denote the energy associated with an atom of species $\alpha$ interacting with an atom of species $\alpha'$ at the $n$th nearest-neighbor distance. We find that it is sufficient to fit interactions to the first four-neighbor distances to accurately capture the data produced using the $S^{(2)}$ theory.

## APPENDIX B: POTENTIAL TRAINING & TESTING

Figure 12 shows the training and testing errors for forces and stresses for the GAP, compressed GAP, and spin GAP. In all cases, the spin GAP is seen to perform the best (lower errors). The errors at which these trends plateau were discussed previously in Sec. II E. The composition of the database is reported in Table I.
[1] J. G. Parr, A. Hanson, and R. Lula, *Stainless Steel* (American Society for Metals, Metals Park, OH, 1993).

[2] P. Marshall, *Austenitic Stainless Steels: Microstructure and Mechanical Properties* (Elsevier Applied Science Publishing, London, 1984).

[3] H. Ledbetter and M. Austin, *Mater. Sci. Eng.* **70**, 143 (1985).

[4] G. S. Was, *Fundamentals of Radiation Materials Science* (Springer, Berlin, 2007).

[5] J. B. Piochaud, T. P. C. Klaver, G. Adjanor, P. Olsson, C. Domain, and C. S. Becquart, *Phys. Rev. B* **89**, 024101 (2014).

[6] O. Dimitrov and C. Dimitrov, *J. Nucl. Mater.* **105**, 39 (1982).

[7] C. Dimitrov and O. Dimitrov, *J. Phys. F* **14**, 793 (1984).

[8] M. Benkaddour, D. Bahia, Y. Belkadi, J. Barkani, and A. Benkaddour, *Mater. Today: Proceedings* **45**, 7359 (2021).

[9] W. Kohn and L. J. Sham, *Phys. Rev.* **140**, A1133 (1965).

[10] A. Manzoor, Y. Zhang, and D. S. Aidhy, *Comput. Mater. Sci.* **198**, 110669 (2021).

[11] Y. Zhang, A. Manzoor, C. Jiang, D. Aidhy, and D. Schwen, *Comput. Mater. Sci.* **190**, 110308 (2021).

[12] E. Antillon, N. Bernstein, and M. D. Johannes, *Phys. Rev. Mater.* **6**, 053801 (2022).

[13] R. A. Johnson, *Phys. Rev. B* **39**, 12554 (1989).

[14] B. Jelinek, S. Groh, M. F. Horstemeyer, J. Houze, S. G. Kim, G. J. Wagner, A. Moitra, and M. I. Baskes, *Phys. Rev. B* **85**, 245102 (2012).

[15] H. Rafii-Tabar and A. P. Sulton, *Philos. Mag. Lett.* **63**, 217 (1991).

[16] R. W. Smith and G. S. Was, *Phys. Rev. B* **40**, 10322 (1989).

[17] M. Grujicic and X. Zhou, *Calphad* **17**, 383 (1993).

[18] G. Bonny, D. Terentyev, R. C. Pasianot, S. Poncé, and A. Bakaev, *Modell. Simul. Mater. Sci. Eng.* **19**, 085008 (2011).

[19] G. Bonny, N. Castin, and D. Terentyev, *Modell. Simul. Mater. Sci. Eng.* **21**, 085004 (2013).

[20] A. Daramola, G. Bonny, G. Adjanor, C. Domain, G. Monnet, and A. Fraczkiewicz, *Comput. Mater. Sci.* **203**, 111165 (2022).

[21] V. L. Deringer, M. A. Caro, and G. Csányi, *Adv. Mater.* **31**, 1902765 (2019).

[22] C. van der Oord, M. Sachs, D. P. Kovács, C. Ortner, and G. Csányi, *npj Comput. Mater.* **9**, 168 (2023).

[23] A. P. Bartók, M. C. Payne, R. Kondor, and G. Csányi, *Phys. Rev. Lett.* **104**, 136403 (2010).

[24] R. Drautz, *Phys. Rev. B* **99**, 014104 (2019).

[25] A. V. Shapeev, *Multiscale Model. Simul.* **14**, 1153 (2016).

[26] I. Batatia, D. P. Kovács, G. N. C. Simm, C. Ortner, and G. Csányi, *Adv. Neural Inf. Process. Syst.* **35**, 11423 (2022).

[27] S. Batzner, A. Musaelian, L. Sun, M. Geiger, J. P. Mailoa, M. Kornbluth, N. Molinari, T. E. Smidt, and B. Kozinsky, *Nat. Commun.* **13**, 2453 (2022).

[28] R. Lot, F. Pellegrini, Y. Shaidu, and E. Küçükbenli, *Comput. Phys. Commun.* **256**, 107402 (2020).

[29] A. P. Bartók, J. Kermode, N. Bernstein, and G. Csányi, *Phys. Rev. X* **8**, 041048 (2018).

[30] W. J. Szlachta, A. P. Bartók, and G. Csányi, *Phys. Rev. B* **90**, 104108 (2014).

[31] D. Dragoni, T. D. Daff, G. Csányi, and N. Marzari, *Phys. Rev. Mater.* **2**, 013808 (2018).

[32] V. L. Deringer and G. Csányi, *Phys. Rev. B* **95**, 094203 (2017).

[33] C. W. Rosenbrock, K. Gubaev, A. V. Shapeev, L. B. Pártay, N. Bernstein, G. Csányi, and G. L. W. Hart, *npj Comput. Mater.* **7**, 24 (2021).

[34] F. C. Mocanu, K. Konstantinou, T. H. Lee, N. Bernstein, V. L. Deringer, G. Csányi, and S. R. Elliott, *J. Phys. Chem. B* **122**, 8998 (2018).

[35] A. P. Bartók, R. Kondor, and G. Csányi, *Phys. Rev. B* **87**, 184115 (2013).

[36] W. L. Bragg and E. J. Williams, *Proc. R. Soc. London, Ser. A* **145**, 699 (1934).

[37] W. L. Bragg and E. J. Williams, *Proc. R. Soc. London, Ser. A* **151**, 540 (1935).

[38] S. N. Khan, J. B. Staunton, and G. M. Stocks, *Phys. Rev. B* **93**, 054206 (2016).

[39] C. D. Woodgate and J. B. Staunton, *Phys. Rev. B* **105**, 115124 (2022).

[40] C. D. Woodgate and J. B. Staunton, *Phys. Rev. Mater.* **7**, 013801 (2023).

[41] C. D. Woodgate, D. Hedlund, L. H. Lewis, and J. B. Staunton, *Phys. Rev. Mater.* **7**, 053801 (2023).

[42] C. D. Woodgate, L. H. Lewis, and J. B. Staunton, *arXiv:2401.02809*.

[43] J. S. Faulkner, G. M. Stocks, and Y. Wang, *Multiple Scattering Theory Electronic Structure of Solids*, 1st ed. (IOP Publishing, Bristol, 2018).

[44] M. Hoffmann, A. Ernst, W. Hergert, V. N. Antonov, W. A. Adeagbo, R. M. Geilhufe, and H. Ben Hamed, *Physica Status Solidi (b)* **257**, 1900671 (2020).

[45] G. Kresse and J. Hafner, *Phys. Rev. B* **49**, 14251 (1994).

[46] G. Kresse and J. Furthmüller, *Comput. Mater. Sci.* **6**, 15 (1996).

[47] G. Kresse and J. Furthmüller, *Phys. Rev. B* **54**, 11169 (1996).

[48] J. P. Perdew, K. Burke, and M. Ernzerhof, *Phys. Rev. Lett.* **77**, 3865 (1996).

[49] P. Olsson, C. Becquart, and C. Domain, *Mater. Res. Lett.* **4**, 219 (2016).

[50] G. Kresse and D. Joubert, *Phys. Rev. B* **59**, 1758 (1999).

[51] S. Cottenier, B. D. Vries, J. Meersschaut, and M. Rots, *J. Phys.: Condens. Matter* **14**, 3275 (2002).

[52] R. Hafner, D. Spišák, R. Lorenz, and J. Hafner, *Phys. Rev. B* **65**, 184432 (2002).

[53] H. J. Monkhorst and J. D. Pack, *Phys. Rev. B* **13**, 5188 (1976).

[54] S. Nosé, *J. Chem. Phys.* **81**, 511 (1984).

[55] W. G. Hoover, *Phys. Rev. A* **31**, 1695 (1985).

[56] C. E. Rasmussen and C. K. I. Williams, *Gaussian Processes for Machine Learning* (MIT Press, Cambridge, MA, 2006).

[57] S. Klawohn, J. P. Darby, J. R. Kermode, G. Csányi, M. A. Caro, and A. P. Bartók, *J. Chem. Phys.* **159**, 174108 (2023).

[58] S. Klawohn, J. R. Kermode, and A. P. Bartók, *Mach. Learn.: Sci. Technol.* **4**, 015020 (2023).

[59] J. P. Darby, D. P. Kovács, I. Batatia, M. A. Caro, G. L. W. Hart, C. Ortner, and G. Csányi, *Phys. Rev. Lett.* **131**, 028001 (2023).

[60] A. Marucco, *Mater. Sci. Eng.*, **A** **189**, 267 (1994).

[61] H. M. Ledbetter, *J. Appl. Phys.* **52**, 1587 (1981).

[62] V. Rajevac, M. Hoelzel, S. A. Danilkin, A. Hoser, and H. Fuess, *J. Phys.: Condens. Matter* **16**, 2609 (2004).

[63] P. Grigorev et al., *J. Open Source Software* **9**, 93 (2024).

[64] M. Rinaldi, M. Mrovec, A. Bochkarev, Y. Lysogorskiy, and R. Drautz, *npj Comput. Mater.* **10**, 12 (2024).

[65] https://doi.org/10.5281/zenodo.10577535.
[1] J. G. Parr, A. Hanson, and R. Lula, *Stainless Steel* (American Society for Metals, Metals Park, OH, 1993).

[2] P. Marshall, *Austenitic Stainless Steels: Microstructure and Mechanical Properties* (Elsevier Applied Science Publishing, London, 1984).

[3] H. Ledbetter and M. Austin, *Mater. Sci. Eng.* **70**, 143 (1985).

[4] G. S. Was, *Fundamentals of Radiation Materials Science* (Springer, Berlin, 2007).

[5] J. B. Piochaud, T. P. C. Klaver, G. Adjanor, P. Olsson, C. Domain, and C. S. Becquart, *Phys. Rev. B* **89**, 024101 (2014).

[6] O. Dimitrov and C. Dimitrov, *J. Nucl. Mater.* **105**, 39 (1982).

[7] C. Dimitrov and O. Dimitrov, *J. Phys. F* **14**, 793 (1984).

[8] M. Benkaddour, D. Bahia, Y. Belkadi, J. Barkani, and A. Benkaddour, *Mater. Today: Proceedings* **45**, 7359 (2021).

[9] W. Kohn and L. J. Sham, *Phys. Rev.* **140**, A1133 (1965).

[10] A. Manzoor, Y. Zhang, and D. S. Aidhy, *Comput. Mater. Sci.* **198**, 110669 (2021).

[11] Y. Zhang, A. Manzoor, C. Jiang, D. Aidhy, and D. Schwen, *Comput. Mater. Sci.* **190**, 110308 (2021).

[12] E. Antillon, N. Bernstein, and M. D. Johannes, *Phys. Rev. Mater.* **6**, 053801 (2022).

[13] R. A. Johnson, *Phys. Rev. B* **39**, 12554 (1989).

[14] B. Jelinek, S. Groh, M. F. Horstemeyer, J. Houze, S. G. Kim, G. J. Wagner, A. Moitra, and M. I. Baskes, *Phys. Rev. B* **85**, 245102 (2012).

[15] H. Rafii-Tabar and A. P. Sulton, *Philos. Mag. Lett.* **63**, 217 (1991).

[16] R. W. Smith and G. S. Was, *Phys. Rev. B* **40**, 10322 (1989).

[17] M. Grujicic and X. Zhou, *Calphad* **17**, 383 (1993).

[18] G. Bonny, D. Terentyev, R. C. Pasianot, S. Poncé, and A. Bakaev, *Modell. Simul. Mater. Sci. Eng.* **19**, 085008 (2011).

[19] G. Bonny, N. Castin, and D. Terentyev, *Modell. Simul. Mater. Sci. Eng.* **21**, 085004 (2013).

[20] A. Daramola, G. Bonny, G. Adjanor, C. Domain, G. Monnet, and A. Fraczkiewicz, *Comput. Mater. Sci.* **203**, 111165 (2022).

[21] V. L. Deringer, M. A. Caro, and G. Csányi, *Adv. Mater.* **31**, 1902765 (2019).

[22] C. van der Oord, M. Sachs, D. P. Kovács, C. Ortner, and G. Csányi, *npj Comput. Mater.* **9**, 168 (2023).

[23] A. P. Bartók, M. C. Payne, R. Kondor, and G. Csányi, *Phys. Rev. Lett.* **104**, 136403 (2010).

[24] R. Drautz, *Phys. Rev. B* **99**, 014104 (2019).

[25] A. V. Shapeev, *Multiscale Model. Simul.* **14**, 1153 (2016).

[26] I. Batatia, D. P. Kovács, G. N. C. Simm, C. Ortner, and G. Csányi, *Adv. Neural Inf. Process. Syst.* **35**, 11423 (2022).

[27] S. Batzner, A. Musaelian, L. Sun, M. Geiger, J. P. Mailoa, M. Kornbluth, N. Molinari, T. E. Smidt, and B. Kozinsky, *Nat. Commun.* **13**, 2453 (2022).

[28] R. Lot, F. Pellegrini, Y. Shaidu, and E. Küçükbenli, *Comput. Phys. Commun.* **256**, 107402 (2020).

[29] A. P. Bartók, J. Kermode, N. Bernstein, and G. Csányi, *Phys. Rev. X* **8**, 041048 (2018).

[30] W. J. Szlachta, A. P. Bartók, and G. Csányi, *Phys. Rev. B* **90**, 104108 (2014).

[31] D. Dragoni, T. D. Daff, G. Csányi, and N. Marzari, *Phys. Rev. Mater.* **2**, 013808 (2018).

[32] V. L. Deringer and G. Csányi, *Phys. Rev. B* **95**, 094203 (2017).

[33] C. W. Rosenbrock, K. Gubaev, A. V. Shapeev, L. B. Pártay, N. Bernstein, G. Csányi, and G. L. W. Hart, *npj Comput. Mater.* **7**, 24 (2021).

[34] F. C. Mocanu, K. Konstantinou, T. H. Lee, N. Bernstein, V. L. Deringer, G. Csányi, and S. R. Elliott, *J. Phys. Chem. B* **122**, 8998 (2018).

[35] A. P. Bartók, R. Kondor, and G. Csányi, *Phys. Rev. B* **87**, 184115 (2013).

[36] W. L. Bragg and E. J. Williams, *Proc. R. Soc. London, Ser. A* **145**, 699 (1934).

[37] W. L. Bragg and E. J. Williams, *Proc. R. Soc. London, Ser. A* **151**, 540 (1935).

[38] S. N. Khan, J. B. Staunton, and G. M. Stocks, *Phys. Rev. B* **93**, 054206 (2016).

[39] C. D. Woodgate and J. B. Staunton, *Phys. Rev. B* **105**, 115124 (2022).

[40] C. D. Woodgate and J. B. Staunton, *Phys. Rev. Mater.* **7**, 013801 (2023).

[41] C. D. Woodgate, D. Hedlund, L. H. Lewis, and J. B. Staunton, *Phys. Rev. Mater.* **7**, 053801 (2023).

[42] C. D. Woodgate, L. H. Lewis, and J. B. Staunton, *arXiv:2401.02809*.

[43] J. S. Faulkner, G. M. Stocks, and Y. Wang, *Multiple Scattering Theory Electronic Structure of Solids*, 1st ed. (IOP Publishing, Bristol, 2018).

[44] M. Hoffmann, A. Ernst, W. Hergert, V. N. Antonov, W. A. Adeagbo, R. M. Geilhufe, and H. Ben Hamed, *Physica Status Solidi (b)* **257**, 1900671 (2020).

[45] G. Kresse and J. Hafner, *Phys. Rev. B* **49**, 14251 (1994).

[46] G. Kresse and J. Furthmüller, *Comput. Mater. Sci.* **6**, 15 (1996).

[47] G. Kresse and J. Furthmüller, *Phys. Rev. B* **54**, 11169 (1996).

[48] J. P. Perdew, K. Burke, and M. Ernzerhof, *Phys. Rev. Lett.* **77**, 3865 (1996).

[49] P. Olsson, C. Becquart, and C. Domain, *Mater. Res. Lett.* **4**, 219 (2016).

[50] G. Kresse and D. Joubert, *Phys. Rev. B* **59**, 1758 (1999).

[51] S. Cottenier, B. D. Vries, J. Meersschaut, and M. Rots, *J. Phys.: Condens. Matter* **14**, 3275 (2002).

[52] R. Hafner, D. Spišák, R. Lorenz, and J. Hafner, *Phys. Rev. B* **65**, 184432 (2002).

[53] H. J. Monkhorst and J. D. Pack, *Phys. Rev. B* **13**, 5188 (1976).

[54] S. Nosé, *J. Chem. Phys.* **81**, 511 (1984).

[55] W. G. Hoover, *Phys. Rev. A* **31**, 1695 (1985).

[56] C. E. Rasmussen and C. K. I. Williams, *Gaussian Processes for Machine Learning* (MIT Press, Cambridge, MA, 2006).

[57] S. Klawohn, J. P. Darby, J. R. Kermode, G. Csányi, M. A. Caro, and A. P. Bartók, *J. Chem. Phys.* **159**, 174108 (2023).

[58] S. Klawohn, J. R. Kermode, and A. P. Bartók, *Mach. Learn.: Sci. Technol.* **4**, 015020 (2023).

[59] J. P. Darby, D. P. Kovács, I. Batatia, M. A. Caro, G. L. W. Hart, C. Ortner, and G. Csányi, *Phys. Rev. Lett.* **131**, 028001 (2023).

[60] A. Marucco, *Mater. Sci. Eng.*, **A** **189**, 267 (1994).

[61] H. M. Ledbetter, *J. Appl. Phys.* **52**, 1587 (1981).

[62] V. Rajevac, M. Hoelzel, S. A. Danilkin, A. Hoser, and H. Fuess, *J. Phys.: Condens. Matter* **16**, 2609 (2004).

[63] P. Grigorev *et al.*, *J. Open Source Software* **9**, 93 (2024).

[64] M. Rinaldi, M. Mrovec, A. Bochkarev, Y. Lysogorskiy, and R. Drautz, *npj Comput. Mater.* **10**, 12 (2024).

[65] https://doi.org/10.5281/zenodo.10577535.