Non-collinear Magnetic Atomic Cluster Expansion for Iron

The Atomic Cluster Expansion (ACE) provides a formally complete basis for the local atomic environment. ACE is not limited to representing energies as a function of atomic positions and chemical species, but can be generalized to vectorial or tensorial properties and to incorporate further degrees of freedom (DOF). This is crucial for magnetic materials with potential energy surfaces that depend on atomic positions and atomic magnetic moments simultaneously. In this work, we employ the ACE formalism to develop a non-collinear magnetic ACE parametrization for the prototypical magnetic element Fe. The model is trained on a broad range of collinear and non-collinear magnetic structures calculated using spin density functional theory. We demonstrate that the non-collinear magnetic ACE is able to reproduce not only ground state properties of various magnetic phases of Fe but also the magnetic and lattice excitations that are essential for a correct description of the finite temperature behavior and properties of crystal defects.


I. INTRODUCTION
Recent advancements of data-driven methods and machine-learned (ML) interatomic potentials have led to dramatically improved descriptions of the potential energy surface (PES) for many material systems.However, the incorporation of spin degrees of freedom (DOF), which are crucial to capture finite temperature phenomena in magnetic materials, has remained a challenging endeavour.In spin density functional theory (SDFT), magnetizaton emerges from the competition of magnetic exchange and band energy contributions [61,62], where the energy required for reshuffling electrons in up and down spin channels depends on the local density of states (DOS).The bimodal DOS of iron in the body-centred crystal (bcc) structure affords large DOS values close to the Fermi level, leading to larger magnetic moments than in the face-centred cubic (fcc) structure with its more unimodal DOS that is lower at the Fermi level [63,64].This intricate interplay between magnetic and atomic structure implies that multi-atom multi-spin interactions are necessary for capturing different magnetic and atomic structures in a single model.
Unlike approaches that were derived from electronic structure theory and that seamlessly incorporate the complexity of magnetic interactions [65][66][67], classical interatomic potentials needed to be supplemented via suitable interaction terms that mimic the quantum exchange interactions.The simplest possibility was to employ a classical Heisenberg Hamiltonian [68], where the atomic spin operators are substituted by spin vectors and the exchange interactions are parameterized using firstprinciples calculations [69].Such strategies have been adopted also in most existing ML approaches for magnetic systems.
Nikolov et al. [70] augmented the spectral neighborhood analysis potential (SNAP) framework with a twospin bi-linear Heisenberg model with atomic magnetic moment magnitudes being fixed and independent of the environment.A similar approach, where a neural network was trained to describe contributions to the Heisenberg Hamiltonian based on the local magnetic environment, was developed by Yu et al. [71].However, this approach did not include information about the underlying lattice and treated the magnetic moments as unit vectors.Eckhoff et al. [72] extended the formalism based on Behler-Parrinello symmetry functions [73] in a framework that was limited to collinear configurations.Magnetic moments as additional DOF were incorporated by Novikov et al. [74] in the moment tensor potential framework [75].Even though the description was confined to collinear moments only, the magnetic moment tensor potential was able to reproduce a number of thermodynamic properties of bulk bcc Fe.Very recently, Domina et al. [76] extended the SNAP framework to deal with arbitrary vectorial fields and demonstrated its functionality by training to non-collinear spin configurations generated using a model Landau-Heisenberg Hamiltonian.Finally, aiming at large-scale spin-lattice dynamics simulations, Chapman et al. [77] added a neural network correction term to an embedded atom method potential augmented with a Heisenberg-Landau Hamiltonian.The model was successfully applied in finite temperature simulations of bulk Fe phases as well as complex defects.However, due to its simplicity, absolute errors were in some cases larger than a few tens of meV that are comparable to the fluctuations of exchange parameters with temperature.Thus, none of the existing magnetic ML approaches has so far succeeded in achieving a transferable and quantitatively accurate description of magnetic interactions suitable for modelling magnetism in different crystal structures.
We present an explicit treatment of non-collinear magnetic DOF within the atomic cluster expansion (ACE) [78,79].ACE provides a complete basis in the space of atomic environments [78,80] and accurate, transferable and numerically efficient parameterizations of ACE have been developed for diverse bonding environments including bulk metallic systems as well as covalent molecules [81][82][83].Thanks to ACE universality, additional scalar, vectorial or tensorial DOF can be incorporated seamlessly into ACE models [79].Specifically for magnetic systems, ACE provides a body-ordered decomposition of combined atomic and magnetic PES in terms of a complete set of basis functions that depend on atomic and magnetic DOF.The inclusion of magnetic DOF requires an extension of the ACE equivariant basis such that any transformation of the relevant translation and rotation symmetry group acting on both atomic and magnetic spaces leaves the energy invariant.Magnetic ACE can therefore be considered as a generalization of most existing magnetic ML models as well as the classical spin-cluster expansion (SCE) [84][85][86][87][88].
In this work, we develop a non-collinear magnetic ACE parameterization for the prototypical magnetic element Fe.The model is trained on a large dataset of both collinear and non-collinear DFT calculations and validated for a broad range of structural, thermodynamic and defect properties.The resulting interatomic potential is able to describe accurately complex potential energy landscapes of different magnetic and atomic phases of Fe as a function of both atomic positions and local magnetic moment vectors.

A. Reference DFT data
A comprehensive sampling of variations in both atomic positions and magnetic moments is crucial for the construction of any atomistic magnetic ML model.Sampling of the atomic DOF can be carried out following well established protocols employed in ML fitting of PES, commonly by choosing a set of structures and varying their geometries and atomic positions.In contrast, sampling of the magnetic DOF presents a significant difficulty, both from the computational and methodological point of view.Firstly, the number of required calculations increases drastically due to the additional spin degrees of freedom and, secondly, the local atomic magnetic moments need to be constrained to desired magnitudes and orientations [89].While it is in principle possible to fix both the direction and the magnitude of each atomic magnetic moment to a target vector [89], these calculations are computationally demanding.Furthermore, as atomic magnetic moments are computed by integrating over a sphere, different magnetization densities within the sphere may in principle lead to identical moments.
To generate the training dataset for magnetic ACE, we considered both conventional, unconstrained and collinear as well as constrained non-collinear spin po-FIG.1: Constant magnetic moment energy-volume curves for FM bcc and fcc phases computed using constrained DFT.The black curve marks the ground state configurations without any applied constrain.The two minima for fcc correspond to the high-and low-spin magnetic configurations.larized configurations.
These configurations ranged from various spin spirals in ideal bcc cells to supercells with random orientations of the moments and perturbed atomic positions.For bulk phases along the Bain transformation path, we sampled the magnitudes of the collinear magnetic moment over the whole physically reasonable range from 0 to ∼3 µ B /atom.The simultaneous sampling of both atomic and magnetic DOF enabled to generate a set of uniformly distributed configurations that are relevant for the properties of interest for a wide range of atomic densities as well as magnitudes and directions of the atomic magnetic moments.An example of data collected with this strategy is given in Fig. 1 for the bcc and fcc ferromagnetic (FM) phases.Each data point corresponds to the energy of either structure at a given volume and a constrained value of the magnetic moment.The ground state configurations are marked by the black curve.While the bcc phase has only one minimum, corresponding to the ground state FM bcc phase, the fcc phase exhibits two minima corresponding to high-spin and low-spin configurations.
The constrained magnetic calculations required convergence of the energy and forces with respect to a constraining penalty term [89].In some limited cases it was computationally prohibitive to achieve numerically small penalty contributions, mainly for configurations far from equilibrium such as highly distorted structures and defects (see SI for representative examples).Therefore, we excluded configurations for which the penalty energy was larger than ≈ 5 meV/atom as these would significantly increase the noise in the data and adversely affect the parameterization.
The resulting training dataset contained about 70 000 structures in total that can be divided into several categories, each associated with a particular property of interest.The categories are listed in Table I, where we specify the number of configurations and the range of volumes and magnetic moment magnitudes that we considered.The free atom data, obtained from calculations of a bcc unit cell with lattice parameter equal to 12 Å at different magnitudes of the magnetic moment, is used to fit the first order contribution of the expansion that characterizes the asymptotic large volume limit for each structure with a given magnetic moment.Detailed information about the free atom reference is provided in the supplementary material.

B. Training procedure
The fitting of the magnetic ACE potential for Fe was done following procedures that were established for the nonmagnetic ACE [81,90].A hierarchical basis extension was employed, starting from one-body contribution and adding gradually contributions with higher body orders.The first step of the fitting procedure consisted in the parameterization of the expansion coefficients for the first order contribution using the free magnetic atom data.This term can be reduced to a Ginzburg-Landau expansion n A n m 2n i , where a maximum number of three terms is commonly used [88,[91][92][93][94].In our parameterization, excellent agreement with the reference data (∼2 meV/atom error) could be obtained using four terms (n = 4).After the first order contribution was fixed, we fitted second order contributions.The ACE second order contributions are formally equivalent to a distance dependent Heisenberg Hamiltonian , and its bicubic term for l ′ max = 1, 2 and 3, respectively (see Methods).In addition, a third-order magnetic contribution, analogous to a screened three-spin interaction ijk K ijk (m i • m j )(m j • m k ), was also included.Angular contributions in higher order magnetic terms did not improve the fit significantly and were neglected, which reduced the number of basis functions significantly.Radial and angular indices for the atomic contributions were then incremented following a hierarchical basis expansion scheme [90], where contributions with increasing body order were gradually added.The cutoff distance of the ACE was set to 4.5 Å.However, it can be extended if necessary [85].Additional hyperparameters, relevant to magnetic DOF only, include the magnetic cutoff m cut set to 4 µ B , which defines the upper bound of the possible magnitude of atomic magnetic moments, and the upper bounds for magnetic radial functions and spherical harmonics n ′ max and l ′ max for each body order (see Methods for details).
The resulting model consists of 6519 parameters and its overall accuracy is equal to 8 meV/atom and 37 meV/ Å for energies and forces, respectively.The main limiting factor in reducing the error further was numerical noise in the reference data that originated from the magnetic moment confinement procedure.In addition, another parameterization was constructed with a particular focus on defect properties (see supplementary material).

C. Properties of the non-collinear magnetic ACE
We carried out a thorough validation of magnetic ACE against the reference DFT data and evaluated a broad range of properties of various bulk Fe phases that were not included explicitly in the training.The predicted volume-energy curves for the bcc and fcc magnetic and non-magnetic (NM) Fe phases are plotted in Fig. 2, where the corresponding cohesive energies are given with respect to the non-magnetic free atom.It is obvious that ACE predictions agree closely with the reference DFT data for all considered magnetic and non-magnetic phases, including the portion of the magnetic energy landscape where the NM to magnetic transitions take place.Moreover, our potential is able to distinguish subtly different magnetic states within one structure, such as the low-spin and high-spin states of the FM fcc structure.FIG.2: Volume-energy curves for both magnetic and non-magnetic structures of bcc (left) and fcc (right) with corresponding DFT data (small circles).
Variations of the magnetic energy as a function of magnetic moment are displayed for FM bcc in Fig. 3, where each curve corresponds to a constant volume.As expected, these dependencies are positive and monotonic for small volumes (dark blue curves), while above a certain critical volume their behavior qualitatively changes to include a minimum at finite value of the magnetic moment in analogy to a Landau expansion.In the limit of large volumes (dark red curves), the magnetic energy approaches the free atom value.Graphs for other bcc and fcc structures are given in the supplementary material.
Two contour plots of magnetic PES for FM bcc and fcc phases are shown in Fig. 4.These plots demonstrate that ACE can capture simultaneously PES of different phases over a broad range of volumes and magnetic moments (from zero up to ≈ 3.2 µ B ).In agreement with DFT, the bcc phase has a single global minimum at the corresponding equilibrium volume and magnetic moment, while the FM fcc phase exhibits two local minima corresponding to low-spin and high-spin states.
Equilibrium properties of the most important bulk Fe phases are listed in Table II.Further properties, such as phonon spectra and magnetic moment variations, are presented in the supplementary material.As one can see, TABLE I: A list of categories associated with a given target property.To each category we provide the total number of structures, the volume range in percentage of the equilibrium volume V 0 of the corresponding phase, the range of sampled magnetic moments, and the number of atoms in the simulation cell.In the case of both constrained and unconstrained supercell calculations, only the direction of the magnetic moments was fixed while their magnitude was self-consistently converged to the equilibrium value.Defects were calculated without any constrains regarding the atomic magnetic moments.the equilibrium lattice parameters, magnetic moments, and elastic constants of the magnetic phases are in good agreement with the reference DFT values.Larger discrepancies exist for the non-magnetic phases since only very few of these configurations were included in the training dataset.
The Bain transformation path is closely related to the bcc-fcc phase transformation.In the case of Fe, the energetics of this transformation depends sensitively on the magnetic state of both phases [95][96][97][98].Variations of energy along the Bain path, computed at the FM bcc equilibrium volume for different magnetic phases of Fe are shown in Fig. 5 for ACE and DFT.Unlike the ground state FM bcc phase, the AFM and NM bcc phases are not mechanically stable with respect to tetragonal distortion, as reflected by negative values of the magnetic states are almost identical, but both phases are unstable.The minimum energy AFM structure is a body-centered tetragonal phase with c/a ≈ 1.45.The excellent agreement between ACE and DFT for the Bain path is due incorporating the coupling between the magnetic and lattice DOF correctly, which is anomalously strong in Fe [97].
The energy barrier for spin rotations depends sensitively on angular interactions between atomic magnetic moments and changes in moment magnitudes.Here we demonstrate that ACE captures the energetics of spin rotation between FM and AFM bcc phases.In Fig. 6, we show the energetics associated with the rotation of one magnetic moment in a two-atom bcc cell.As the moment on the central atom is rotated, the magnetic configuration gradually transforms from FM to AFM.The contour plot in Fig. 6(b) depicts PES as a function of volume and rotation angle.The black arrow marks the minimum energy path between the FM and AFM phases.While the equilibrium volumes of both phases are not very different, the magnetic moment of the AFM phase is significantly lower than that of the FM phase (cf.also correctly reproduced by ACE, as shown in Fig. 6(c), where we plot the rotation energy barriers evaluated at constant magnetic moments.The minimum energy path (dashed grey curve), corresponding to a reduction of the absolute value of magnetic moment from 2.22 µ B in FM bcc to 1.25 µ B in AFM bcc, is in excellent agreement with the DFT reference (black points).
The energy of magnetic moment orientations that deviate only slightly from the collinear alignment can be described by lowest order contributions only, i.e., a bilinear Heisenberg model.From the distance dependent exchange interactions J ij in the bilinear Heisenberg model, the magnon spectrum can be obtained in adiabatic approximation as We determined the exchange interactions for different coordination shells following the real space approach by Liechtenstein et al. [99][100][101], where infinitesimal perturbations to the directions of two neighboring magnetic moments are applied.Calculating the energy δE ij for rotating two spins at atomic sites i and j by opposite infinitesimal angles ±θ/2 and comparing to the energy for rotating the two spins individually, δE i and δE j , results in The distance dependent exchange interactions are then obtained by fitting δE ij − (δE i + δE j ) with respect to the tilting angle for consecutive coordination shells in a large supercell.The resulting adiabatic magnon spectrum is shown in Fig. 7 with reference data obtained using the spin-polarized relativistic Korringa-Kohn-Rostoker (SPRKKR) [102] framework.Small discrepancies between the ACE and SPRKKR results visible for some high frequencies in the magnon spectrum can be at-tributed to the long range part of the magnetic interactions neglected in the present ACE parameterization.Nevertheless, the overall agreement is good, indicating the ability of our parameterization to describe spin spirals with different frequencies.

D. Finite temperature and defects
The magnetic ACE can be applied in large-scale finite temperature simulations to investigate properties that depend on both spin and lattice DOF.The ACE prediction of the FM to paramagnetic (PM) phase transition in bcc Fe is presented in Fig. 8.In a simulation with 3456 atoms, we employed coupled molecular dynamics (MD) -Monte Carlo (MC) sampling [103], where the atoms follow Langevin dynamics while MC is employed for updating the directions of the atomic magnetic moments at constant magnitudes.A direct simulation of the dynamics of the combined atomic and magnetic system is difficult due to the lack of numerically stable and efficient symplectic integrators [104].The predicted Curie temperature T C ≈ 950 K is somewhat underestimated in comparison with the experimental value of 1043 K.The variation of magnetization with temperature shown in Fig. 8 is consistent with previous theoretical studies [91,[105][106][107].
To demonstrate that ACE is able to capture properties of crystal defects, we also included several defect configurations in the DFT training data.However, as discussed in Sec.II A, it is often not possible to reach sufficiently small penalty energies in the constrained DFT calculations for such distorted configurations.Therefore, we needed to resort in many cases to unconstrained spinpolarized calculations only, which limited the sampling of the magnetic PES for defects.
Here we present results for three types of defects -a monovacancy, generalized stacking faults, and a screw dislocation.For most defects, the Heisenberg model is insufficient and it is necessary to include higher-order terms in the magnetic Hamiltonian [108].In addition, an accurate reproduction of defect properties can be achieved only if the coupling between spin and lattice excitations is taken into account.
The monovacancy formation and migration energies of 2.57 eV and 0.65 eV, respectively, agree well with the reference DFT data (equal to 2.17 and 0.67 eV, respectively).The generalized stacking fault energy surface, the γ-surface, for the {110} plane is shown in Fig. 9.The figure further shows cuts along the ⟨111⟩ direction on both the {110} and {211} planes that are related to atomic structures of 1  2 ⟨111⟩ screw dislocations.The ACE predictions for both cuts are in excellent agreement with the DFT reference.ACE also predicts the core structure of the 1 2 ⟨111⟩ screw dislocation, which governs the low temperature plasticity of Fe, in quantitative agreement with DFT reference [109,110] as demonstrated in Fig. 10.
To examine the ability of ACE to reproduce properties of other defects, we constructed a smaller training dataset containing also surfaces and interstitials and generated another ACE parameterization.As shown in the supplementary material, ACE can reproduce properties of these defects as well.

III. DISCUSSION
By incorporating magnetic DOF in the form of atomic magnetic moment vectors into ACE, we demonstrated that constrained non-collinear DFT reference data can be reproduced with excellent accuracy and transferability, exceeding those of existing magnetic ML interatomic potentials.We constructed a non-collinear ACE parametrization for Fe and validated it for a wide range of properties, including volume-energy curves, elastic moduli, phonon spectra, Bain transformation paths, spin rotations and magnon spectra, and point and extended defects.These tests showed that magnetic ACE is not only able to capture large structural and magnetic variations but also resolves subtle spin fluctuations that are crucial for a correct reproduction of phase transitions and thermodynamic properties.To this end it is necessary to include multi-atom multi-spin interactions that are missing in simple models with pairwise couplings between atoms and/or magnetic moments.In iron magnetic angular contributions of body order four and higher are numerically small and can be neglected.
The magnetic ACE was parameterized from DFT reference data that was generated by constraining both the magnitude and direction of the atomic magnetic moments.For configurations with defects or significant atomic displacements, it was often difficult to achieve self-consistency.Furthermore, as the constraints were implemented by integrating magnetization over spheres about atoms, various intra-atomic magnetization distributions could result in the same atomic magnetic moment.This non-uniqueness effectively led to noise in the DFT reference data that ultimately limited the accuracy of our parameterization.Therefore, there is a strong incentive to implement more advanced constraints in DFT that would help to increase the accuracy of magnetic ACE as well as other magnetic ML approaches.
The numerical efficiency of ACE enables to carry out large-scale molecular and spin dynamics simulations to study the dynamics of combined magnetic and structural phase transitions.Nevertheless, to predict the magnetic transition in bcc iron, we employed MD simulations for the atomic positions and MC sampling to vary the atomic magnetic moments.One of the reasons is that there is a lack of classical or semi-classical equations of motion and corresponding, numerically robust integrators applicable to combined atomic and spin dynamics in systems with multi-atom multi-spin interactions and including changes of magnitudes of magnetic moments [104].
The ACE for iron can be extend directly to multicomponent systems, such as technologically important magnetic alloys and carbides.While this is straightforward from a formal point of view, the generation of accurate and comprehensive DFT reference data for magnetic multicomponent materials is challenging.Here efficient sampling based on D-optimality active learning [111] extended to include magnetic DOF will help to reduce the number of required DFT reference calculations.

IV. METHODS
We provide a summary of the magnetic ACE formalism together with aspects of its implementation.Further explanations on implementation and workflow are available in the supplementary material.We also provide computational details of the DFT calculations and the combined MD-MC simulations that were employed for the calculation of the FM-PM transition.

A. Energies, forces and magnetic gradients
We define state variables σ ji of atom j neighboring atom i in terms of interatomic distances vectors r ji , chemical species µ j , , magnetic moments m j , etc. as with σ ii = (µ i , m i ).A neighbor density on atom i including atomic and magnetic contributions can then be written as Magnetic contributions also enter the single bond basis functions, where v = (knlmn ′ l ′ m ′ ) and the primed indices are used to label basis functions that depend on magnetic contributions, and The projection of the density in Eq. ( 4) on the corresponding single atom basis functions leads to the atomic basis A iv and A (0) iv and From the two atomic bases the tensor product basis is formed and symmetrized to ensure invariance with respect to rotation and inversion, leading to equivariant basis functions where C is a sparse matrix of products of the Clebsch-Gordan coefficients of the atomic and magnetic systems.
The coupling tree, used to form possible tuples v (see SI for an example), can be simplified assuming that spinorbit coupling can be neglected.This is typically an excellent approximation as the spin-orbit coupling energy is on the order of a few µeV for iron bulk systems.Then the atomic and magnetic systems can be completely decoupled and the total angular momenta of the atomic and magnetic channels couple to zero individually, leading to a significant reduction in the number of basis functions (see SI for a detailed explanation).A further reduction of the allowed combinations of atomic and magnetic indices can be obtained by requiring inversion invariance for both atomic and magnetic spaces by restricting the sum of the corresponding angular momenta to even numbers.
We represent the energy for atom i including atomic and magnetic contributions as a linear expansion where c is the vector of the expansion coefficients.The energy can be rewritten in terms of the c basis introduced in [79,81] as This expansion was used to fit the DFT energies and forces.In order that the expression reduces to the non-magnetic ACE when the magnetic moments are zero, the first order equivariant basis was taken as iµin ′ (m = 0)=1 by our choice of magnetic radial functions (see the following Sec.IV B).
Expressions for forces and magnetic gradients are obtained by taking the derivative of the energy with respect to atomic positions and magnetic moments, respectively, and are written in a compact notation as and The pairwise atomic forces f ki are given by and magnetic forces t k and t ki by and The calculation of the adjoints ω iµinlmn ′ l ′ m ′ and ω iµin ′ l ′ m ′ can be further decomposed to the evaluation of two distinct terms, where The adjoint ω iµin ′ l ′ m ′ does not contain the onsite basis contribution and is simply given by The summation over N in Eq. ( 20) starts from zero because even a single atom contributes to the total magnetic gradient.

B. Magnetic radial functions
The magnetic radial functions M µj µi n ′ l ′ used in this work exhibit a different functional form to their atomic counterparts that are given in terms of Chebyshev polynomials [78,90]).In particular, one has to ensure that the energy is invariant under time reversal symmetry, i.e., m i → −m i for every i.For these reasons, we chose a linear combination of Chebyshev polynomials T k as with The scaled distance x guarantees the invariance under time reversal symmetry where m cut is the cutoff for the magnetic moment magnitude.The expansion coefficients c µj µi n ′ l ′ k ′ for both magnetic and atomic radial functions are adjusted during the fitting procedure.
The constrained local moment approach [89] was employed to constrain either both size and direction or just the direction of the atomic magnetic moments.The exchange-correlation energy was represented using the Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation (GGA) method [117].We carried out carefully converged calculations with tight settings of the principal parameters in order to obtain accurate results for the energy, forces and magnetic moments.Specifically, the kinetic energy cutoff was set to 500 eV, the convergence threshold for the energy to 10 −5 eV and the k-mesh density to 0.18 Å−1 .The integration radius for the atomic magnetic moments (VASP parameter RWIGS) was kept constant at the value of the Fe PAW (1.302 Å).See SI for the convergence of magnetic moment magnitude with respect to the integration radius and for a discussion on the convergence of the penalty energy in the constrained local moment method.

D. MD-MC calculations
The MD-MC simulations of the FM-PM transition in bcc Fe consisted of alternating MD and MC steps.The MD simulations were performed using Langevin dynamics (from ASE [118] package) with a time step of 1 fs.The MC sampling included uniform spin rotations on a unit sphere.The simulation supercell had dimensions 12 × 12 × 12 of a bcc cell and contained 3456 atoms.The dimensions of the supercell were kept fixed at all temperatures so that the effect of thermal expansion was neglected.At each temperature, we carried out about 10 7 steps, with the initial 10% used for equilibration.

V. DATA AVAILABILITY
The authors declare that all data supporting the findings of this study are available from MR on reasonable request.

VI. CODE AVAILABILITY
The magnetic ACE presented in the manuscript was evaluated and parameterized with a magnetic version of the PACEmaker code [90].The magnetic ACE implementation is not yet distributed with PACEmaker, but we plan to make it available in a future release.

VII. AUTHOR CONTRIBUTIONS
MR implemented the magnetic ACE by extending nonmagnetic ACE implementations in Fortran and in Tensorflow.MR produced the DFT reference data with support from MM and fitted the magnetic ACE.MR carried out all calculations presented in the manuscript and drafted the manuscript.MM and RD supervised the work, advised on theory and edited the manuscript.AB, YL and RD supported the implementation of magnetic ACE and advised on numerical aspects of the implementation.

Calculation of the onsite atomic basis
using radial functions and spherical harmonics of the atom i; 3. Loop over the neighbors j of atom i to calculate the atomic basis A iµnlmn ′ l ′ m ′ ; 4. Loop over basis function indices to calculate products of the atomic basis for each body order (i.e.Eq. ( 9)) and the quantities dA iµnlmn ′ l ′ m ′ of Eqs. ( 19) and ( 21); 5.At this point the single atom energies can be evaluated; 6.The matrix Θ (N ) µiµnln ′ l ′ is determined from the expansion coefficients and then the adjoints of Eqs.(17) and (20) are calculated; 7. The atomic and magnetic forces can be finally assembled.
The final computational cost for the determination of the single atom energies ε i (and in general of any scalar) depends on the evaluation of the onsite A (0) basis, the atomic basis A and on the computation of the B basis.The scaling is ) respectively, where N c is the number of neighbors within the chosen cutoff radius and N max is the maximum order.

IV. TUNING OF THE INTEGRATION RADIUS
The parameter RWIGS in VASP can be adjusted for bulk phases to the corresponding Wigner-Seitz radius.However, this procedure is not easily applicable to distorted structures.Nevertheless, the atomic magnetic moment (and also the charge associated with an atom) depends on the value of the integration radius.In order to illustrate this point, we show in Fig. S2 the charge and the magnetic moment vs the integration radius for a single Fe atom in a large supercell.The results are obtained by projecting s, p and d atomic-like orbitals on the self-consistent charge density.Both the charge and the magnetic moment exhibit a considerable variation with respect to the integration radius, where the d electrons account for the largest contribution.For this orbital projection, the most appropriate integration radius for a free atom corresponds to about 2.5 Å, which results in a charge of almost 8 valence electrons and an atomic magnetic moment of almost 4 µ B .
Figure S3 shows similar dependencies for the FM bcc structure with lattice parameters equal to 2.85 and 4.0 Å, respectively.It is clear from these two plots that in the case of the bulk bcc phase, and generally for other bulk structures, the magnitude of the magnetic moment does not change much for moderate changes of the Wigner-Seitz radius, namely, between the default pseudopotential value of 1.302 Å and the actual Wigner-Seitz radius corresponding to the given volume.Therefore, to avoid any ambiguity, we decided to use the same (default) integration radius in all calculations to obtain consistent values of magnetic moments for arbitrary configuration.

V. ISSUES WITH THE CONVERGENCE OF THE PENALTY ENERGY TERM FOR CONSTRAINED DFT CALCULATIONS
A representative example showing a map of the penalty energy for antiferromagnetic (AFM) bcc as function of volume and constrained magnetic moment is shown in Fig. S4.One can see that the region of lowest penalty energy (∼1 meV/atom) is within a curved trench (dark red color) close to the lowest energy magnetic configurations as a function of volume (black curve).Away from this region the penalty energy increases up to to about 4.5 meV/atom for large-volume/small-moment configurations (the lower right plateau region) or more for configurations characterized by large magnetic moments and small volumes (the upper left blue region).Therefore, for our training, we restricted the dataset to regions characterized by low magnitudes of the penalty energy, i.e. by excluding both large and low volume data.

VI. SUPPLEMENTARY PLOTS
The learning curve for the magnetic ACE training is shown in Fig. S5 together with the correlation plot of predicted vs reference energies.
The equilibrium curves for the magnetic energy vs magnetic moment magnitude are displayed in Fig. S6, where the different magnetic phases of bcc and fcc show varying degrees of localization.For instance, in the case of bcc FM the moment at the equilibrium volume is sufficiently localized if compared to fcc FM.The latter shows indeed a more FIG.S2: Charge and magnetic moment per atom vs integration radius for a single iron atom in a large supercell.The atom has 8 valence electrons and a magnetic moment equal to 4µ B .The crosses mark the default integration radius RWIGS of 1.302 Å of the pseudopotential.pronounced itinerant character, leading to energetically favorable LSFs at finite temperatures and larger spread in the distribution of the moments [S119, S120].The related equilibrium moment magnitude vs volume curves are given in Fig. S7, where the equilibrium state for all the magnetic phases changes rapidly from the NM to the magnetic solution in a range of volumes from ≈ 6 Å3 /atom to 10 Å3 /atom and stabilizes around 12 Å3 /atom increasing then more gently up to the free atom limit.In the case of fcc FM the high-and low-spin states are also visible at the corresponding volumes.The non-smooth behavior of these curves is determined by the details of the integrated density of states (DOS) for the spin up and down channels, which are given by the competition between electronic kinetic and Coulomb energy.The FM phonon spectrum at the equilibrium magnetic moment is displayed, together with the corresponding DOS, in Fig. S8, where also an excellent agreement is obtained with the reference DFT spectrum.

VII. SPECIFIC MODEL FITTED ON DEFECTS ONLY
In this section we discuss the predictions obtained by fitting ACE on defects and bcc FM only.The chosen cutoff is equal to 6.0 Å and RMSE are 14 meV for energies and 153 meV/ Å for forces.A list of defect formation energies is shown in Tab.I together with the reference DFT data.The agreement is good for all the defects properties despite the lack of sampling of the magnetic DOF due to convergence problems of the constrained local moment method in this specific case.In Tab.II the predicted distribution of the magnetic moments around the monovacancy is shown with the reference DFT data.In the particular case of the (112) surface we also compared the values of the relaxed moments as a function of the layer depth, as shown in Fig. S9, to the reference.The agreement for this particular surface is within an acceptable range of ± 0.1 µ B /atom.The predictions for the 1/2 [111] screw dislocation core structure and γ-surface are identical to the results obtained with the potential presented in the main text.FIG.S3: Charge and magnetic moment per atom vs integration radius for the FM bcc Fe with lattice parameter equal to 2.85 Å (top panels) and 4.0 Å (bottom panels).The crosses mark the default integration radius RWIGS radius of 1.302 Å from the pseudopotential.The vertical dashed lines mark the Wigner-Seitz radius for the corresponding volume.

FIG. 3 :
FIG. 3: Magnetic energy vs magnetic moment magnitude at different volume for bcc FM.Dashed lines and black dots correspond to equilibrium volume ACE and reference DFT data, respectively.
FIG. 4: Contour plots for the bcc (left) and fcc (right) FM phases.
FIG. 5: Bain transformation paths between the FM, AFM and NM bcc and fcc phases (ACE: lines, DFT: circles).

FIG. 6 :
FIG. 6: Analysis of the FM to AFM transformation in the bcc phase via rotation of the spin on the central atom: a) A schematic picture of the transformation.b) A contour plot of PES as a function of volume vs rotation angle.c) FM to AFM spin rotation energy barriers at constant magnetic moment.The minimum energy path is marked by the grey dashed curve together with the DFT reference (black points).

FIG. 8 :
FIG. 8: Magnetization vs temperature.The vertical dashed line indicates the estimated Curie temperature T C .The experimental value of T C is 1043 K. Insets show snapshots of parts of the simulation cell for better visualization.
FIG. S4: Map of the penalty energy for bcc AFM.The equilibrium configurations are marked by the black curve.
FIG.S6: Magnetic energy vs magnetic moment magnitude curves for the different magnetic phases of bcc and fcc together with the DFT reference.

TABLE II :
Equilibrium properties of bcc and fcc phases of

TABLE I :
List of defect formation energies as predicted by ACE when both atoms and magnetic moments are relaxed in comparison with our reference DFT values.In the case of interstitials and surfaces the predicted formation energies are calculated by relaxing only the positions and keeping the moments fixed to the equilibrium bulk value.

TABLE II :
Values of the magnetic moments in µ B for the 1 st , 2 nd and 3 rd nearest-neighbor of the vacancy site.