Quantum-accurate magneto-elastic predictions with classical spin-lattice dynamics

A data-driven framework is presented for building magneto-elastic machine-learning interatomic potentials (ML-IAPs) for large-scale spin-lattice dynamics simulations. The magneto-elastic ML-IAPs are constructed by coupling a collective atomic spin model with an ML-IAP. Together they represent a potential energy surface from which the mechanical forces on the atoms and the precession dynamics of the atomic spins are computed. Both the atomic spin model and the ML-IAP are parametrized on data from first-principles calculations. We demonstrate the efficacy of our data-driven framework across magneto-structural phase transitions by generating a magneto-elastic ML-IAP for {\alpha}-iron. The combined potential energy surface yields excellent agreement with first-principles magneto-elastic calculations and quantitative predictions of diverse materials properties including bulk modulus, magnetization, and specific heat across the ferromagnetic-paramagnetic phase transition.


Introduction
Magnetism strongly influences thermomechanical properties in a large variety of materials, such as single-element magnetic metals 1, 2 , steels 3 , high-entropy alloys 4,5 , nuclear fuels such as uranium dioxide 6 , magnetic oxides 7,8 , and numerous other classes of functional materials 9 . Despite the critical role of magnetism in the aforementioned materials classes, modeling efforts to study the interplay between structural and magnetic properties have been notably lacking. Furthermore, there are unanswered scientific questions regarding the significance of magnetism in matter that is shock-compressed 10,11 or exposed to strong electromagnetic fields such as in coherent lights sources 12,13 , pulsed power and high magnetic fields facilities 14,15 . Properties of interest include phase transitions, thermal stability of magnetic defects, magneto-mechanical couplings, but many of these subjects are challenging or prohibited by state of the art computational tools.
A prime but simple example of the computational advance made herein is the heat-capacity of α-iron displayed in Figure 1. The experimental measurement of the heat capacity C p diverges at the magnetic Curie transition, characteristic of a second-order phase transition 18 . Without a scalable coupled spin-lattice dynamics simulation environment, that properly accounts for thermal expansion and magnetic contribution  16,17 , the red squares our simulation results, and black dashed line indicates the experimental Curie transition temperature. This illustrates the well-known ferromagneticparamagnetic phase transition, where the heat capacity diverges at the Curie temperature.
to the pressure, reproducing the divergence of C p (and of other thermomechanical properties) at the critical point is not possible. Accurate numerical simulations are critical for enabling technological advances, as they shape our fundamental understanding of the underlying solid state physics that dictates material behavior. Developing high fidelity models however is challenging, because it necessitates capturing physical phenomena that occur across several length and time scales. This can only be achieved with sufficiently accurate multiscale simulation tools 19,20 , which is the focus of this work. Classical molecular dynamics (MD) simulations 21 provide a useful framework for multiscale modeling by leveraging interatomic potentials (IAPs) to represent the dynamics of atoms on a Born-Oppenheimer potential energy surface (PES) 22 . By utilizing massively parallel algorithms 23 and long time-scale methodologies 24 , MD enables bridging first-principles with continuum-scale simulations 25 .
The absorption of machine learning (ML) techniques into the creation of interatomic potentials has lead to classical MD simulations that approach the accuracy of first-principles methods. A large number of these highly accurate ML-IAPs [26][27][28][29][30][31][32] have been developed. In general, they are parameterized on training data (configuration energy, atomic forces) from first-principles methods like density functional theory (DFT) 33 and utilize different flavors of ML model forms to construct the PES. While they have proven to be useful for large-scale simulations of materials properties 34,35 , further progress in multiscale modeling is hampered by the limitation of ML-IAPs to non-magnetic materials phenomena. Even with highly accurate ML-IAPs, state-of-the-art MD simulations cannot reproduce the divergent behavior of C p near the critical point ( Figure 1) because they fail to account for the magnetic degrees of freedom 36 .
Coupling atomic spin dynamics with classical MD has been pioneered by Ma et al. [37][38][39] . Herein, a classical magnetic spin is assigned to each atom in addition to its position leading to a 6N-dimensional PES (5N if the magnetic spin norms are fixed), instead of the common 3N-dimensional PES in classical MD: where r r r i j = r r r i −r r r j denotes the relative position between atoms i and j, s s s i the classical spin assigned to atom i, and N the number of atoms in the system. In most classical spin-lattice calculations, the 6N-dimensional PES is constructed by introducing an atomic spin model on top of a mechanical IAP 37 . For example, a common approach is to combine a distancedependent Heisenberg Hamiltonian with an embedded-atommethod (EAM) potential 39,40 . While these prior approaches recover experimental properties on a qualitative level 41,42 , their combined representation of phononic and magnetic degrees of freedom is not sufficiently consistent for providing quantitative predictions at the level of first-principles results. More recently, Ma et al. developed a magneto-elastic IAP for magnetic iron based on data from first-principles calculations 43 . However, this remained an isolated attempt as there is no general methodology for generating a magneto-elastic PES in a classical context that enables large-scale spin-lattice dynamics simulations for any magnetic material.
In this work, we overcome this methodological obstacle by providing a data-driven framework for generating magnetoelastic ML-IAPs that (1) provide a consistent representation of both mechanical and magnetic degrees of freedom and (2) achieve near first-principles accuracy. We refer to our new class of IAPs as "magneto-elastic ML-IAPs" as they generate a consistent PES accurately representing the magnetic degrees of freedom and the interplay between magnetic and elastic phenomena. Our framework couples an atomic spin model (Heisenberg Hamiltonian) with an ML-IAP and provides a unified magneto-elastic PES which yields the correct mechanical forces on the atoms in the MD framework. The Heisenberg Hamiltonian is parameterized with data from DFT spin-spiral calculations at different degrees of lattice compression. In constructing the ML-IAP, we leverage the flexible and data-driven spectral neighbor analysis potential (SNAP) methodology 32 which is trained on a database of magnetic configurations generated using DFT calculations.
We apply our framework to generate a magneto-elastic ML-IAP for the α phase of iron. We demonstrate that our potential is transferable to an extended area of the phase diagram, corresponding to a temperature and pressure range of 0 to 1200 K and 0 to 13 GPa (up to the α → γ and α → ε transitions, respectively). The Curie temperature, which experimentally occurs at approximately 1045 K, lies within this parameter space. After presenting our training workflow, the "Results" section will probe the "quantum-accuracy" of our magneto-elastic ML-IAP by performing magneto-static comparisons to first-principles measurements. We then stress that our generated magneto-elastic ML-IAP can also be directly used in the LAMMPS package 23 to perform magnetodynamic simulations that take into account both the thermal expansion of the lattice and magnetic pressure due to spin disorder. This enables us to maintain a constant ambient pressure throughout all calculations of thermomechanical properties, consistent with conditions prevalent in experiments. As illustrated in Figure 1, our framework allows us to perform the first pressure-controlled quantitative prediction of the critical behavior across a second-order phase transition within a classical spin-lattice dynamics simulation.

Results
In this section we outline our advancements in magnetic materials modeling. We first present our training workflow and subsequently assess our results by comparing both static and dynamic properties in α-iron against first-principles calculations and experiments. Figure 2 displays our training workflow. Further details to each box in this diagram are presented as a subsection in the "Methods" section. All atomic configurations in the training set result from first-principles calculations performed with the same DFT setup (same pseudo-potential and energy cutoff, Figure 2. Magneto-elastic ML-IAP training workflow. A training set of DFT calculations is partitioned into those that train the SNAP interatomic potential and those that train the spin Hamiltonian, respectively. A non-magnetic interatomic potential is fit to configuration energies and atomic forces after the spin Hamiltonian contribution is subtracted and is validated against magneto-elastic properties computed in LAMMPS. Optimization of the spin Hamiltonian and interatomic potential parameters is handled by DAKOTA. similar k-point densities) as detailed in the "Methods" section. In contrast to traditional force-matching approaches in the development of classical IAPs, we treat the magnetic and phononic degrees of freedom in the PES in a consistent and unified manner, as indicated by the exchange of information between spin Hamiltonian and SNAP potential parametrization steps. After parameterizing our atomic spin Hamiltonian by leveraging DFT spin-spiral results, its energy, forces, and stress contributions are subtracted from each atomic configuration in the first-principles training set. The ML-IAP is then trained to reproduce the non-magnetic component of the first-principles data. Finally, both components of the magneto-elastic PES are recombined to construct a unified magneto-elastic ML-IAP that is consistently trained on firstprinciples data. Optimization is handled by the DAKOTA software package 44 in both fitting steps. For the SNAP component of the potential, DAKOTA optimizes the radial cutoff of the interaction along with the weights of each training data set (energy and force weights) to generate different candidate potentials. Those candidates potentials are then recombined with the spin Hamiltonian and tested against selected objective functions (mean-absolute errors (MAEs) in lattice constants, cohesive energies, elastic constants, forces and total energies). Table 1 summarizes the different groups of training data, the optimal weights obtained for each of those groups, and the corresponding energy and force MAEs. The target values for the objective functions are based on both experimental and DFT data, as outlined in Table 2. Objective function evaluations are done within LAMMPS 23 .
Herein, the critical innovation that enables a leap forward in predictive simulations of magnetic materials is this datadriven workflow. Magnetic and phononic contributions to the PES are taken into account explicitly and any miscounting is avoided (for example, no double counting of the magnetic energy or contribution to the pressure). The obtained magnetoelastic ML-IAP can directly be used to run spin-lattice calculations in LAMMPS 23,40,45 .

Magneto-Static Accuracy
We first assess the quantitative agreement of our magnetoelastic ML-IAP by comparing with DFT results where magnetic order and elastic deformations are coupled. This is done by leveraging a particular subset of spin configurations referred to as spin-spirals, for which the energy and corresponding pressure can be evaluated from both DFT and classical magneto-elastic potential calculations. Details about definition and computation of spin-spirals can be found in the "Methods" section. Equation-of-state calculations (energy and pressure versus volume) are performed at the Γ point (corresponding to the purely ferromagnetic state) and for spin-spirals corresponding to q-vectors along the ΓH and ΓP high-symmetry lines. The calculations at the Γ point represent the magnetic ground state and, hence, serve as a point of reference for the spin spiral calculations. The geometric orientation of the various computed spin spirals is visualized in Figure 3. The first set (q = 0.01 along ΓH and q = 0.07 along ΓP) represents "long" spirals, close to the Γ point, the second set (q = 0.1 along ΓH and q = 0.14 in ΓP) represents spirals with intermediate periodicity, and the last set (q = 0.2 along ΓH and q = 0.21 along ΓP) is chosen close to the borders of the magnetic training set (see red demarcation lines in Figure 5 in the "Methods" section). The DFT results are obtained by leveraging the generalized Bloch theorem, whereas our classical spin-lattice calculations were performed by generating the corresponding supercells (details given in the "Methods" section).
Excellent agreement between our classical spin-lattice model and DFT is achieved at the Γ point and for the two first q-vectors on each high-symmetry line (q = 0.01 and q = 0.1 along ΓH, q = 0.07 and q = 0.14 along ΓP) in the pressure range relevant for the α-phase of iron (up to 13 GPa which corresponds to the α → ε transition). At higher q-vector values, the energy and pressure predictions of our atomic spin-lattice model still agree reasonably well with the DFT calculations. The observed deviation from the DFT results can be explained by the limitations of our atomic spin-lattice model: as both the pressure and the relative angle between neighboring spins increase, fluctuations of the atomic spin norms become more important. As discussed in the "Methods" and "Discussion" sections, these are not included in the Hamiltonian of our atomic spin-lattice model.

Magneto-Dynamic Accuracy
Turning now to spin-lattice dynamics calculations based on our magneto-elastic ML-IAP (as detailed in the "Methods" section), we assess the quantitative accuracy with respect to experimental measurements of changes in magnetic and thermoelastic properties as the material is heated. In making this comparison, it is necessary to choose which thermodynamic state variables will be held fixed and which will be allowed to vary with temperature. Spin-lattice dynamics algorithms have been developed for simulations in a canonical ensemble (CE) which preserves the number of particles, the volume, and the temperature in the system 39 . Our first set of simulation conditions, referred to as "fixed-volume conditions" (FVC), hold the volume fixed while running dynamics in the CE at specified values of the lattice and spin temperatures. A disadvantage of this choice is that the pressure steadily increases as heat is added to the material, in contradiction to the experimental observations, which are conducted at constant pressure. To this date, an isobaric spin-lattice algorithm has not been developed (preserving the system's pressure rather than its volume). However, our methodology as implemented in LAMMPS enables us to compute the magnetic contribution to the pressure. By alternating thermalization (coupled spinlattice dynamics in a CE) and pressure equilibration (frozen spin configuration in an isobaric ensemble) steps, it is possible to control the pressure of our spin-lattice system. Hence, we refer to calculations performed in this pressure-controlled CE as "pressure-controlled conditions" (PCC). In both conditions, the temperature of the spin and lattice subsystems is set using two separate Langevin thermostats (one acting on the spins, the other on the lattice) 39 . Finally, this enables us to define a third set of conditions: in addition to controlling the pressure, the spin thermostat can be set to match a given magnetization value (i.e., the experimental magnetization) rather than a temperature. We refer to this as "pressure-controlled and magnetization-controlled conditions" (PCMCC). Figure 6 in the "Methods" section displays the different definitions of the spin temperature and the evolution of the pressure for those three different conditions.
In practice, FVC, PCC and PCMCC only differ in their equilibration conditions (control of pressure and / or magnetization), as each of the corresponding simulations are performed in a canonical ensemble. We illustrate the predictive capability of our magneto-elastic ML-IAP in α-iron for these equilibration conditions in Figure 4.a-f (FVC : , PCC : , PCMCC : ). The agreement of the following magneto-elastic properties with experimental results is assessed: magnetization (Figure 4.a), heat-capacity C p (Figure 4.b), thermal expansion (cell volume on Figure 4.c), bulk modulus (Figure 4.d), and two shear constants, (c 11 − c 12 )/2 and c 44 (Figure 4.e-f). The "Spin-Lattice Dynamics" subsection of the "Methods" section details the computation of those temperature-dependent elastic constants.
We first work under the FVC ( ), keeping a constant volume and equal spin and lattice temperatures (Figure 4.c and Figure 6). At constant volume, our model predicts a Curie temperature of approximately 716K (Figure 4.a). Specific heat calculations shown in Figure 4.b were performed by computing the derivative of the internal energy, taking both the lattice and magnetic contributions into account. The SNAP contribution (lattice only) was first isolated and determined to be a constant value of 26.4 Jmol −1 K −1 , in good agreement with the Dulong-Petit value of 3R 49 . The magnetic contribution offsets the total specific heat at low temperature, as the magnetization steadily decreases (thus steadily increasing the magnetic energy). Also at low temperature, deviation between simulations and experiment (highlighted by the semitransparent blue region in Figure 4.b) occurs due to quantum effects which reduce the experimental heat capacity below the 3R value. The FVC heat-capacity is determined at constant volume, although we use the symbol C p on the axis label because the enhanced simulations described below are indeed conducted at constant pressure conditions. In those constant volume conditions, the pressure evolution with temperature increase is substantial (up to 12 GPa, almost corresponding to the α → ε transition, as can be seen on Figure 6), which has a strong impact on the underlying elastic properties. Interestingly, at the Curie temperature (here 716K), the increasing pressure exhibits an inflection point, confirming the importance of spin fluctuations on the thermoelastic properties. The temperature dependence of three elastic constants is shown in Figure 4.d-f. For the bulk modulus, FVC does not agree well with experimental data, especially at higher temperatures. The FVC results tend to overestimate the stiffness, which most likely arises from the build-up of thermal stresses in the material. Under these conditions a nearly temperature-invariant c 44 response is predicted, which is in strong contrast to trends in experiment. Despite these shortcomings, the FVC calculations actually match the experimental data for shear constant (c 11 − c 12 )/2 relatively well throughout the entire temperature range. In general, the fixed volume assumption made under FVC fails to account for thermal expansion, leading to incorrect elastic predictions.
We correct this shortcoming of the model by working under PCC ( ) which allows for thermal expansion. As can be seen on Figure 4.c, the cell volumes are relaxed at each finite temperature, until the pressure in the system drops to approximately 0 GPa. As shown in Figure 4.a, the thermal expansion incorrectly moves the onset of Curie transition to approximately 536K. As the average interatomic distance increases, the strength of the exchange interaction is lowered, thus decreasing the transition temperature. The computed heat-capacity (Figure 4.b) now corresponds to the derivative of the free energy, and to an actual C p measurement. However, as in the FVC, the low agreement between the experimental and computed magnetization evolution leads to an offset in the initial C p and does not match the Dulong-Petit value at low temperature. The PCC fares better in reproducing the experimental bulk modulus up to the Curie transition (no hardening observed). PCC also does better in terms of the shear constant c 44 , as it is able to reproduce the thermal softening seen in experiments. However, for shear constant (c 11 − c 12 )/2, PCC underestimates the extent of the thermal softening. Overall, PCC does better than FVC in terms of elastic properties, but deviates more in terms of magnetic predictions compared to experiment. By shifting the Curie transition towards lower temperatures, it reduces the range of validity of our elastic calculations.
In order to improve the magnetic predictions of α-iron, we finally consider the PCMCC scheme ( ). In addition to allowing for thermal expansion similarly to the PCC, we also set the spin thermostat temperature in order to reproduce the experimental magnetization. Below the Curie transition, the spin temperature increases more slowly than the lattice temperature, while above the Curie transition, it increases at the same rate as the lattice temperature (see Figure 6 in the "Methods section"). Figure 4.a shows that the obtained magnetization under PCMCC closely matches that of experiment. Most prominently, the resulting C p agrees well with experiments (Figure 4.b). The Dulong-Petit value is recovered at low temperature, and the C p discontinuity at the Curie transition is well captured. The thermal expansion trend is also in much better agreement with experiments, with very comparable slopes between approximately 200 and 750K (Figure 4.c). Up to approximately 600 K, PCMCC agree very well with the experimental values for (c 11 − c 12 )/2 (Figure 4.e) but at 800-1000K a slight hardening is observed, which contradicts experimental data. For the bulk modulus, PCMCC correctly predicts the nearly linear trend up to the Curie temperature.
We note that in all three sets of conditions, a rapid increase of about 25-30 GPa in the bulk modulus is observed as we move across the critical point. This jump was found to be strongly impacted by the underlying mechanical potential. The prediction accuracy could possibly be improved by including additional, finite-temperature objective functions in the fitting procedure. The PCMCC prediction of the shear constant c 44 closely matches the PCC data. This tends to indicate that this shear constant c 44 is not impacted significantly by the spin dynamics. For both pressure controlled conditions (PCC and PCMCC) the maximum deviation from experiments occurs near 700K and is approximately 14%.

Discussion
We presented a data-driven framework for automated generation of magneto-elastic ML-IAPs which enable large-scale spin-lattice dynamics simulations for any magnetic material in LAMMPS. This framework was demonstrated by generating a robust magneto-elastic ML-IAP for α-iron. First we investigated the magneto-static accuracy (energy and pressure) with respect to equivalent first-principles calculations. It was demonstrated that the generated magneto-elastic ML-IAP (which represents the corresponding 5-N dimensional PES) is in close agreement with first-principles magneto-elastic calculations. This was achieved by properly partitioning the PES into magnetic and mechanical degrees of freedom. Subsequently, we investigated the magneto-dynamic accuracy by comparing predicted finite temperature magneto-elastic properties (magnetization, heat-capacity, thermal-expansion, bulk modulus, and shear constants) across the ferromagneticparamagnetic phase transition from spin-lattice dynamics simulations against data from experiments. In the course of this, we analyzed the choice of simulation conditions (control of pressure and magnetization) and highlighted the importance of thermal and magnetic pressure contributions. This is an important advance over traditional classical magnetization dynamics methods, where contributions from thermal expansion or spin pressure due to disorder are negated. We demonstrated that spin-lattice dynamics simulations of controlled pressure and constrained magnetization yields qualitative agreement with the measured magneto-elastic properties.
Our framework enables predictions of critical properties across the second-order phase transition within classical spinlattice dynamics simulations, such as the divergent behavior of the heat capacity around the Curie temperature (Figure 1  and 4.b). We provide a more comprehensive perspective on our results by comparing them within the context of other first-principles and classical methods. At low temperature, first-principles methods can capture the electronic component of the heat-capacity, up to the Dulong-Petit value 49,50 (the difference with our model is highlighted by the blue area on Figure 4.b). However, computing C p across the Curie transition requires a dynamic treatment of large spin-ensembles whose calculation is computationally expensive in terms of first-principles methods. Classical IAPs do not explicitly treat magnetic degrees of freedom and, thus, cannot reproduce the effects of this magnetic phase-transition 51 . An empirical model which is based on first-principles calculations and accounts for electronic, phononic and magnetic degrees of freedom gave excellent agreement with the experimental C p curve of α-iron up to the Curie temperature 52 . However, this model does not extend above the Curie temperature, does not account for the pressure generated by the corresponding spin configurations, and cannot be easily extended to other thermomechanical properties. Thus, for a range of temperature from about 250 K to 1200 K, our model provides with a set of very good predictions, obtained for the computational cost of classical MD calculations only.
We conclude the discussion of our results by pointing out limitations of the present method and future prospects. First, note that the agreement to the experimental Curie transition (T c ≈ 716K in a fixed volume calculation) could have been adjusted by parameterizing the spin potential on a smaller range of the high-symmetry lines (see Figure 5), or by adding an objective function aimed at matching the experimental value in the spin-potential fitting procedure. However, this additional constraint would have worsened the agreement of our model with the DFT energy and pressure results (as displayed on Figure 3) and would contradict the overall objective of this work.
For temperatures below approximately 250 K, our classical framework cannot access the quantized free energy, and is thus enable to accurately reproduce the trends of all the quantities being its derivatives (C p , elastic constants, ...). This is reducing the agreement versus experiments of the magnetodynamic accuracy measurements displayed on Figure 4 at low temperature, and can be seen as a limitation of our classical approach 53 .
Another limitation of our work lies in the simplicity of the spin Hamiltonian model used. Extended spin Hamiltonians, such as spin-cluster expansions, might be a promising route to improving the accuracy of the magnetic component of the PES by both accounting for the fluctuation of the magnetic moment magnitudes and many-body spin interactions 54,55 . A straightforward extension of this work could combine recently developed extended spin Hamiltonians with first-principles studies, and apply our formalism to extend our α-iron magneto-elastic ML-IAP to account for defect configurations 56, 57 , Cr clustering 58,59 , and magneto-structural phase-transitions 11,60 .
Enhanced magnetic thermostats have also been proposed in order to better match the experimental magnetic transition versus temperature 61,62 . Such thermostats could be implemented in LAMMPS and used to replace the magnetization-controlled conditions defined in the "Results" section. This could extend the range validity of our framework to areas of phase-diagrams where the magnetization distribution is not well measured (for example in the ε phase of iron).
A recent study added a magnetic contribution to the set of descriptors used in a moment-tensor ML-IAP 63 . Although this approach does not explicitly simulate the magnetization dynamics (and its effects on thermomechanical properties), the authors demonstrated remarkable improvement in terms of error convergence. At this stage of our work, we believe improving the modeling of the magnetic component of the PES remains our first priority (and thus implementing and fitting improved spin Hamiltonians, as discussed above). However, depending on the success of this first effort, this complementary approach could be leveraged to improve the quantumaccuracy of our magneto-elastic ML-IAPs.
In summary, we have presented a new computational framework for near quantum-accuracy simulations of magnetoelastic materials properties. By leveraging the flexibility of ML-IAPs, our data-driven workflow enables to model the in-terplay between magnetic and phononic dynamics for a large class of magnetic materials. Furthermore, our straightforward connection to the LAMMPS package makes it possible to perform large-scale quantitative magneto-elastic predictions over controlled pressure and temperature spaces, hitherto study unexplored magneto-dynamics properties of materials.

Density functional theory calculations
Parameterizing both the ML-IAP and the magnetic Heisenberg Hamiltonian relies on data computed using spin-dependent DFT calculations. They were performed using VASP 64,65 . In all calculations the PBE 66 exchange-correlation functional was employed. We used PAW pseudopotentials 67 with 8 valence electrons and a core radius of r c = 2.3 a B . The plane wave cutoff was set to 320 eV and the convergence in each selfconsistency cycle was set to 10 −8 . The Fermi-Dirac smearing scheme with a width of 0.026 eV was used. The Brillouin zone was sampled on a 10 × 10 × 10 grid of k-points. The number of bands used was 224 per atom.

Spin-Spiral Calculations
Spin-spirals define a subset of non-collinear magnetic states. In this work, we leverage spin-spirals as a convenient tool to perform one-to-one comparisons between first-principles and classical magneto-elastic calculations. They can be defined as follows: s s s j = sin θ cos(q q q · R R R 0 j )x x x + sin θ sin(q q q · R R R 0 j )ŷ y y + cos θẑ z z , (2) where q q q is the spin-spiral vector, R R R 0 j is the position of atom j relative to a central atom 0, s s s j is the spin on atom j, and θ is a constant angle between the spins and the spin-spiral vector (often referred to as "cone angle") 68 .x x x,ŷ y y, andẑ z z are the unit vectors along [100], [010], and [001], respectively. Our calculations are restricted to θ = π/2, corresponding to flat spin-spirals in the (001) plane.
First-principles calculations of the per-atom energy and the pressure corresponding to spin-spiral states are performed using DFT by leveraging the frozen-magnon approach 69, 70 and the generalized Bloch theorem 71 as implemented in VASP 72 . We consider a primitive cell of one atom. A 10 × 10 × 10 kpoint grid, an energy cutoff of 320 eV, and 224 bands proved sufficient to reach the level of accuracy expected in our model (as can be seen in Figure 5).
Classical calculations are performed by using Eq. (2) to generate supercells accommodating the spin-spirals corresponding to the q q q-vectors used in the DFT calculations. Based on a given supercell and a spin Hamiltonian, the per-atom energy and pressure are computed using the SPIN package of LAMMPS 23,40 .

Spin Hamiltonian
A spin Hamiltonian is used to model the energy, mechanical forces, and pressure contributions of magnetic configurations.
Rosengaard and Johansson 73 and Szilva et al. 74 showed that adding a biquadratic term to the classical Heisenberg Hamiltonian improves the accuracy of magnetic excitations in 3-d transition ferromagnets. We adopted their Hamiltonian form: where s s s i and s s s j are classical atomic spins of unit length located on atoms i and j, J (r i j ) and K (r i j ) (in eV) are magnetic exchange functions, and r i j is the interatomic distance between magnetic atoms i and j. The two terms in Eq. 3 are offset by subtracting the spin ground state (corresponding to a purely ferromagnetic situation), as detailed in Ma et al. 37 .
Although this offset of the exchange energy does not affect the precession dynamics of the spins, it allows to offset the corresponding mechanical forces. Without this additional term, the magnetic contribution to the forces and the pressure are not zero at the energy ground state. For the exchange interaction terms J (r i j ) and K (r i j ), the interatomic dependence is taken into account through the following function based on an approximation of the Bethe-Slater curve 75,76 : where α denotes the interaction energy, δ the interaction decay length, γ a dimensionless curvature parameter, r a distance, and Θ (R c − r) a Heaviside step function for the radial cutoff R c . This assumes that the interaction decays rapidly with the interatomic distance, consistent with former calculations 74,77 . We set R c = 5Å to include five neighbor shells, as Pajda et al. 77 showed that the exchange interaction decays slower along the [111] direction in α-iron. Using Eq. (3) and leveraging the generalized spin-lattice Poisson bracket as defined by Yang et al. 78 , the magnetic precession vectors (ω ω ω i ), mechanical forces (F F F i ), and their corresponding virial components (W r r r N ) are derived: where r r r N denotes a 3N size vector of all atomic positions and r r r i the position vector of atom i. We note that the virial  80 . In all three plots, the dashed lines correspond to the DFT results, and the continuous lines to our classical model results, whereas the line color (black or blue) corresponds to the lattice compression (0 or 2%, respectively). In the middle plot, the green dashed horizontal line represent the experimental equilibrium value (2.2 µ B per atom), which is the constant value chosen in our model. In all three plots, the red vertical dashed lines are delimiting the q q q-vectors on which our spin Hamiltonian was parametrized.
components enable computing the spin contribution to the pressure.
The spin Hamiltonian is used to reproduce spin-spiral energy and pressure reference results obtained from DFT. They are sampled along two high-symmetry lines, ΓH and ΓP, and for two different lattice constant values (corresponding to the equilibrium bulk value and to a lattice compression of 2%). This allows us to encapsulate in the model the influence of lattice compression on the spin stiffness and the Curie temperature, which was experimentally and theoretically predicted to be small [81][82][83] . Figure 5 displays the excellent agreement obtained between our first-principles spin-spiral energies and experimental measurements.
Our current spin Hamiltonian does not account for fluctuations of the magnetic moment magnitudes, i.e. the norm of atomic spins remains constant in our calculations. As can be seen in Figure 5, this is not the case for our DFT results, as those fluctuations can become important when departing from the Γ point. We thus decided to parameterize our model only on spin-spirals corresponding to q q q-vectors for which the spin norm deviates from the ferromagnetic value (≈ 2.2 µ B /atom at the Γ point) by less than 5%. The red dashed lines in Figure 5 delimit this q q q-vector range.
Finally, we used the single objective genetic algorithm within the DAKOTA software package 44 to optimize the six coefficients of J (r i j ) and K (r i j ) in order to obtain the best possible agreement between our reference DFT spin-spiral energy and pressure results and our spin model. Figure 5 displays the obtained result. As can be seen in Figure 4, for a fixed-volume calculation, our spin Hamiltonian predicts a Curie temperature of 716K. Note that a better match of the DFT spin-spiral energies would yield a larger spin-stiffness, and thus a better agreement for the Curie temperature. However, this would worsen the pressure agreement.
Spin-orbit coupling effects were included by accounting for an iron-type cubic anisotropy 84 : with K 1 = 0.001 eV and K (c) 2 = 0.0005 eV the intensity coefficients corresponding to α-iron. The cubic anisotropy was only included to run calculations, but ignored in the fitting procedure as its intensity is below the range of accuracy of our ML-IAP.
In all our classical spin-lattice dynamics calculations, our system size remained small compared to the typical magnetic domain-wall width in iron 84 . Thus, long-range dipole-dipole interactions could safely be neglected.

SNAP Potential
For this work, an interatomic potential for iron was developed that is specifically parameterized for use in coupled spin and molecular dynamics simulations. Training data for a Spectral Neighborhood Analysis Potential (SNAP) 45,85,86 was collected to constrain the fit to the pressure and temperature phase space of < 20GPa and < 2000K. The set of non-colinear, spin-polarized VASP calculations includes α-(BCC), ε-(HCP) and liquid-iron, Table 1 displays the quantity of each training type and target properties that are captured therein. Optimization of a SNAP potential necessitates that the generated training database be broken into these groups (rows in Table 1) such that the weighted linear regression can (de-)emphasize different parts in search of a global minima in objective function errors. Each training group is assigned a unique weight for its' associated energies and atomic forces for each candidate potential, optimization of these weights is controlled by DAKOTA. Regression is carried out using singular value decomposition with a squared loss function (L2 norm). In order to avoid double counting, and properly simulate the magnetic properties of iron in classical MD, we have adapted the SNAP fitting protocol 45 to isolate the nonmagnetic energy and forces from the generated training data. To do so, the fitted biquadratic spin Hamiltonian is evaluated for every atom in the training set, and its' contribution to the total energy and per-atom forces is subtracted. This is akin to previous uses of an ion core repulsion 87 or electrostatic interaction term 88 as a reference potential while fitting SNAP models.
Optimization of the SNAP potential was achieved using a single objective genetic algorithm within the DAKOTA software package 44 . Radial cutoff distance, training group weights and number of bispectrum descriptors were varied to minimize a set of objective functions, as percent error to available DFT or experimental 90 data, that encapsulate the desired mechanical properties of Fe. These objective functions specific to α-iron are listed in Table 2, and the RMSE energy and force regression errors are included in optimization as well. In all objectives, our linear SNAP model with 31 bispectrum descriptors achieves accuracy in all mechanical properties within a few percent of experiment/DFT. Additionally, lattice constants and cohesive energies of γ-(FCC) and ε-iron (HCP) phases were fit, but given far less priority with respect to the α-iron mechanical properties resulting in ∼ 6 − 7% errors with respect to DFT. Importantly, each of the objective functions were evaluated including the magnetic spin contributions to avoid unforeseen changes in property predictions. A full breakdown of the optimal training group weights and mean absolute energy/force errors are given in Table 1. Group weights listed have been adjusted by the number of configurations or forces they are applied to, therefore allowing for larger group weights to be (cautiously) interpreted more valuable at meeting the set of targeted objective functions. This optimized Fe-SNAP interatomic potential is contained as Supplemental Material along with LAMMPS input scripts used in the following section.

Spin-Lattice Dynamics
Calculations are performed following the spin-lattice dynamics approach as implemented in the SPIN package of LAMMPS 23,40 , and set by the spin-lattice Hamiltonian below: H sl (r r r, p p p, s s s) = H mag (r r r, s s s) where H mag is the spin Hamiltonian defined by the combination of Eq. (3) and Eq. (8). The term V SNAP (r i j ) is our SNAP ML-IAP. The second term on the right in Eq. (9), represents the kinetic energy, where the particle momentum is given as p p p and the mass of particle i is m i . Based on this spinlattice Hamiltonian and leveraging the generalized spin-lattice Poisson bracket as defined by Yang et al. 78 , the equations of motion can be defined as: Particle positions are advanced according to Eq. (10). The derivative of the momentum, given in Eq. (11), is dependent not only on the mechanical potential but the magnetic exchange functions as well. Here γ L is the Langevin damping constant for the lattice and f f f is a fluctuating force following Gaussian statistics given below 40 .
The fluctuating force f f f is coupled to γ L via the fluctuation dissipation theorem as shown in Eq. (14). Here k B is the Boltzmann constant, T l is the lattice temperature, and α and β are coordinates. Shown in Eq. (12) is the stochastic Landau-Lifshitz-Gilbert equation which describes the precessional motion of spins under the influence of thermal noise. In Eq. (12), λ is the transverse damping constant and ω ω ω i is a spin force analogue as shown in Eq. (5). The variable η η η(t) is a random vector whose components are drawn from a Gaussian probability distribution given below: where the amplitude of the noise D S can be related to the temperature of the external spin bath T s according to D S = 2πλ k B T s /h 39 .
SD-MD calculations are carried out using a 20x20x20 BCC cell. The BCC lattice is oriented along each of the coordinate directions. The MD timestep in all cases is set to 0.1 femtoseconds. The damping constants are set to 0.1 (Gilbert damping, no units) for the spin thermostat, and to 0.1 picoseconds for the lattice thermostat. Initially all spins start out aligned in the  90 . Percent error is used as the objective function to avoid artificial importance scaling based on units of the target property.
z-direction. To measure the magnetic properties for the canonical ensemble we initially thermalize the system under NVT dynamics at the target spin/lattice temperatures for 40 picoseconds and then sample the target properties for 10 picoseconds. For pressure-controlled simulations (see PCC and MCPCC in the "Results" section), after the initial 40 picosecond of temperature equilibration we freeze the spin configuration and run isobaric-isothermal NPT dynamics in order to allow the system to thermally expand (still accounting for the effect of the "magnetic" pressure, generated by the spin Hamiltonian). The pressure damping parameter is set to 10 picoseconds. The pressure equilibration run is terminated once the system pressure drops below 0.05 GPa. After this, the spin configuration is unfrozen and another equilibration run is carried out under NVT dynamics for 20 picoseconds. Unfreezing the spin configurations causes a small jump in the pressure, typically within the range of +/-2 GPa. To reduce this pressure fluctuation, a series of uniform isotropic box deformations are performed under the NVE ensemble. During this procedure the box is deformed in 0.02% increments every 2 picoseconds until the magnitude of the pressure is reduced to negligible values (< 10 MPa). Figure 6 displays the pressure profiles obtained within the FVC and PCMCC (similar to the PCC). For the magnetization-controlled conditions (PCMCC in the "Results" section), the spin temperature is adjusted to match the experimental magnetization values. Spin temperature adjustments are made based on the magnetization curve obtained in the pressure-controlled conditions (PCC in the "Results" section). The corresponding spin-lattice temperature relationship is shown in Eqs. (17)(18)(19)(20). Here the fitting coefficients are given as a 1 = 471.6, a 2 = 6362, a 3 = 2774, a 4 = 1119, a 5 = 13.6, a 6 = 1043.3, and a 7 = 0.1, respectively. The functions T s,pre and T s,post prescribe how the spin temperature varies before and after the critical point. At the critical point we use a switching function f sw to smoothly transition from T s,pre to T s,post :

11/15
T s,post (T l ) = T l − a 1 (17) T s,pre (T l ) = a 2 exp − T l − a 3 a 4 2 − a 5 (18) f sw (T l ) = 1 2 1 + tanh T l − a 6 a 7 (19) T s (T l ) = f sw T s,post + (1 − f sw )T s,pre (20) Figure 6 displays the spin temperature profiles for the FVC (and, similarly, the PCC), and the PCMCC. After the magnetic measurements we compute elastic constants by performing both uniaxial and shear deformations along each of the coordinate directions and planes. The magnitude of these deformations in all cases is 2% of the box length. Following each deformation the box is relaxed for 3 picoseconds. After this relaxation the stresses are sampled for 2 picoseconds.

Data Availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Code Availability
The code which was used to train the SNAP potential is available from: https://github.com/FitSNAP/ FitSNAP.