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

Svetoslav Nikolov, Mitchell A. Wood, Atilla Canggi, Jean-Bernanrd Maillet, Mihai Cosmin Marinica, Aidan P. Thompson, Michael P. Desjarlais, Julien Tranchida

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

► To cite this version:

Svetoslav Nikolov, Mitchell A. Wood, Atilla Canggi, Jean-Bernanrd Maillet, Mihai Cosmin Marinica, et al.. Quantum-accurate magneto-elastic predictions with classical spin-lattice dynamics. npj Computational Materials, 2022, 7, pp.153. 10.1038/s41524-021-00617-2 . cea-03651978

HAL Id: cea-03651978

https://cea.hal.science/cea-03651978

Submitted on 26 Apr 2022

**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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Quantum-accurate magneto-elastic predictions with classical spin-lattice dynamics

Svetoslav Nikolov¹, Mitchell A. Wood¹, Attila Cangi², Jean-Bernard Maillet³,⁴,
Mihai-Cosmin Marinica⁵, Aidan P. Thompson¹, Michael P. Desjarlais⁶, and Julien
Tranchida¹,*

¹Computational Multiscale Department, Sandia National Laboratories, P.O. Box 5800, MS 1322, Albuquerque, NM 87185

²Center for Advanced Systems Understanding (CASUS), Helmholtz-Zentrum Dresden-Rossendorf, 02826 Görlitz, Germany

³CEA - DAM, DIF, Arpajon Cedex F-91297, France

⁴Université Paris-Saclay, CEA, LMCE, 91680 Bruyères-le-Châtel, France

⁵Université Paris-Saclay, CEA, Service de Recherches de Métallurgie Physique, Gif-sur-Yvette 91191, France

⁶Sandia National Laboratories, P.O. Box 5800, MS 1322, Albuquerque, NM 87185

*jtranch@sandia.gov

## ABSTRACT

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 α-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¹,², steels³, high-entropy alloys⁴,⁵, nuclear fuels such as uranium dioxide⁶, magnetic oxides⁷,⁸, and numerous other classes of functional materials⁹. 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¹⁰,¹¹ or exposed to strong electromagnetic fields such as in coherent lights sources¹²,¹³, pulsed power and high magnetic fields facilities¹⁴,¹⁵. 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¹⁸. Without a scalable coupled spin-lattice dynamics simulation environment, that properly accounts for thermal expansion and magnetic contribution

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Constant pressure heat capacity of α-iron versus temperature. The black triangles denote experimental measurements¹⁶,¹⁷, the red squares our simulation results, and black dashed line indicates the experimental Curie transition temperature. This illustrates the well-known ferromagnetic-paramagnetic phase transition, where the heat capacity diverges at the Curie temperature.</figcaption></figure>

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-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-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:

$$E = \sum_{i=1}^{N} \epsilon (\{\mathbf{r}_{ij}, \mathbf{s}_i\}), \quad (1)$$

where $\mathbf{r}_{ij} = \mathbf{r}_i - \mathbf{r}_j$ denotes the relative position between atoms i and j, $\mathbf{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 distance-dependent Heisenberg Hamiltonian with an embedded-atom-method (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 magneto-elastic 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 $\alpha$ 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 $\alpha \to \gamma$ and $\alpha \to \epsilon$ 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 magneto-dynamic 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 $\alpha$-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><img src="image_4.png" /><figcaption>**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.</figcaption></figure>

similar k-point densities) as detailed in the "Methods" sec-
tion. 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 parametriza-
tion steps. After parameterizing our atomic spin Hamiltonian
by leveraging DFT spin-spiral results, its energy, forces, and
stress contributions are subtracted from each atomic config-
uration 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 first-
principles data. Optimization is handled by the DAKOTA
software package^44 in both fitting steps. For the SNAP com-
ponent 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 data-
driven 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 magneto-
elastic ML-IAP can directly be used to run spin-lattice calcu-
lations in LAMMPS<sup>23,40,45</sup>.

Magneto-Static Accuracy

We first assess the quantitative agreement of our magneto-
elastic ML-IAP by comparing with DFT results where mag-
netic 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 corre-
sponding pressure can be evaluated from both DFT and clas-
sical magneto-elastic potential calculations. Details about
definition and computation of spin-spirals can be found in
the "Methods" section. Equation-of-state calculations (en-
ergy 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 rep-
resent 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) repre-
sents 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
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Plots of the equation of state data from *first-principles* calculations (VASP computations) and our magneto-elastic ML-IAP (LAMMPS computations) for seven different spin-spirals: a) Γ point b) vectors along the ΓH high-symmetry line, and c) vectors along the ΓP high-symmetry line. Visualizations of the corresponding spin-spiral supercells and associated q-vectors are shown to the right of and above each plot, respectively.</figcaption></figure>

q-vectors on each high-symmetry line ($q = 0.01$ and $q = 0.1$ along $\Gamma H$, $q = 0.07$ and $q = 0.14$ along $\Gamma P$) in the pressure range relevant for the α-phase of iron (up to 13 GPa which corresponds to the $\alpha \rightarrow \epsilon$ 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³⁹. 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
<figure><img src="image_6.png" /></figure>

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

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Plots a-f show magnetoelastic data obtained with our magneto-elastic ML-IAP. The green (○), blue (♦), and red (□) markers indicate the choice of equilibration conditions: "fixed-volume conditions" (FVC), "pressure-controlled conditions" (PCC) and "pressure-controlled and magnetization-controlled conditions" (PCMCC), respectively. In all plots, experimental data (extracted from five different references$^{16,17,46-48}$) is denoted by the filled triangles (▲), and the dotted black lines (···) represent the experimental Curie temperature. The plots in a-b) show magnetization and specific heat comparisons between different ensembles and experiments. The light blue region in (b) indicates the low temperature regime $T \lesssim 250$ K where quantum effects reduce the experimental heat capacity below the classical Dulong-Petit limiting value of $3R^{49}$. The data in plot c) illustrates how the lattice expands with temperature. An inherent offset exists between our model (trained to match the DFT data at 0 K) and experimental measurements. Plots d-f show (d) bulk modulus, (e) $(c_{11} - c_{12})/2$ shear constant, and (f) $c_{44}$ shear constant for the three aforementioned sets of conditions.</figcaption></figure>

to the pressure. By alternating thermalization (coupled spin-lattice 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 equilibra-
tion conditions in Figure 4.a-f (FVC : ○, PCC : ◇, PCMCC
: □). The agreement of the following magneto-elastic prop-
erties with experimental results is assessed: magnetization
(Figure 4.a), heat-capacity C<sub>p</sub> (Figure 4.b), thermal expansion
(cell volume on Figure 4.c), bulk modulus (Figure 4.d), and
two shear constants, (c<sub>11</sub> − c<sub>12</sub>)/2 and c<sub>44</sub> (Figure 4.e-f). The
"Spin-Lattice Dynamics" subsection of the "Methods" sec-
tion details the computation of those temperature-dependent
elastic constants.

We first work under the FVC (○), keeping a constant vol-
ume 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 com-
puting 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<sup>-1</sup>K<sup>-1</sup>, in good agreement
with the Dulong-Petit value of 3R<sup>49</sup>. The magnetic contri-
bution 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 semi-
transparent 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<sub>p</sub> 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. Inter-
estingly, at the Curie temperature (here 716K), the increasing
pressure exhibits an inflection point, confirming the impor-
tance 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 ma-
terial. Under these conditions a nearly temperature-invariant
c<sub>44</sub> response is predicted, which is in strong contrast to trends
in experiment. Despite these shortcomings, the FVC calcula-
tions actually match the experimental data for shear constant
(c<sub>11</sub> − c<sub>12</sub>)/2 relatively well throughout the entire tempera-
ture 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 un-
der 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 in-
creases, 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ₚ measurement. However,
as in the FVC, the low agreement between the experimental
and computed magnetization evolution leads to an offset in
the initial Cₚ and does not match the Dulong-Petit value at low
temperature. The PCC fares better in reproducing the experi-
mental bulk modulus up to the Curie transition (no hardening
observed). PCC also does better in terms of the shear constant
c₄₄, as it is able to reproduce the thermal softening seen in
experiments. However, for shear constant (c₁₁ − c₁₂)/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 exper-
iment. Most prominently, the resulting Cₚ agrees well with
experiments (Figure 4.b). The Dulong-Petit value is recovered
at low temperature, and the Cₚ discontinuity at the Curie tran-
sition is well captured. The thermal expansion trend is also in
much better agreement with experiments, with very compara-
ble slopes between approximately 200 and 750K (Figure 4.c).
Up to approximately 600 K, PCMCC agree very well with
the experimental values for (c₁₁ − c₁₂)/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 in-
cluding additional, finite-temperature objective functions in
the fitting procedure. The PCMCC prediction of the shear
constant c44 closely matches the PCC data. This tends to indi-
cate that this shear constant c44 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 gener-
ation of magneto-elastic ML-IAPs which enable large-scale
spin-lattice dynamics simulations for any magnetic material
in LAMMPS. This framework was demonstrated by gener-
ating a robust magneto-elastic ML-IAP for α-iron. First we
investigated the magneto-static accuracy (energy and pres-
sure) 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. Sub-
sequently, 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 ferromagnetic-
paramagnetic phase transition from spin-lattice dynamics sim-
ulations 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 dy-
namics 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 spin-
lattice 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⁴⁹,⁵⁰ (the
difference with our model is highlighted by the blue area on
Figure 4.b). However, computing Cₚ across the Curie tran-
sition 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⁵¹. An empiri-
cal 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ₚ
curve of α-iron up to the Curie temperature⁵². 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 ther-
momechanical 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 ad-
justed by parameterizing the spin potential on a smaller range
of the high-symmetry lines (see Figure 5), or by adding an ob-
jective 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 quan-
tities being its derivatives (C<sub>p</sub>, elastic constants, ...). This is
reducing the agreement versus experiments of the magneto-
dynamic accuracy measurements displayed on Figure 4 at low
temperature, and can be seen as a limitation of our classical
approach<sup>53</sup>.

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<sup>54,55</sup>. A straight-
forward extension of this work could combine recently devel-
oped extended spin Hamiltonians with *first-principles* studies,
and apply our formalism to extend our α-iron magneto-elastic
ML-IAP to account for defect configurations<sup>56,57</sup>, Cr cluster-
ing<sup>58,59</sup>, and magneto-structural phase-transitions<sup>11,60</sup>.

Enhanced magnetic thermostats have also been proposed in
order to better match the experimental magnetic transition ver-
sus temperature⁶¹,⁶². 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⁶³. 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 im-
proving 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 complemen-
tary approach could be leveraged to improve the quantum-
accuracy of our magneto-elastic ML-IAPs.

In summary, we have presented a new computational frame-
work for near quantum-accuracy simulations of magneto-
elastic 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.

## Methods

### 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 self-consistency 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 \times 10 \times 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:

$$ \mathbf{s}_j = \sin \theta \cos(\mathbf{q} \cdot \mathbf{R}_{0j}) \hat{\mathbf{x}} + \sin \theta \sin(\mathbf{q} \cdot \mathbf{R}_{0j}) \hat{\mathbf{y}} + \cos \theta \hat{\mathbf{z}}, \quad (2) $$

where $\mathbf{q}$ is the spin-spiral vector, $\mathbf{R}_{0j}$ is the position of atom j relative to a central atom 0, $\mathbf{s}_j$ is the spin on atom j, and $\theta$ is a constant angle between the spins and the spin-spiral vector (often referred to as "cone angle")$^{68}$. $\hat{\mathbf{x}}, \hat{\mathbf{y}},$ and $\hat{\mathbf{z}}$ are the unit vectors along [100], [010], and [001], respectively. Our calculations are restricted to $\theta = \pi/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 \times 10 \times 10$ k-point 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 $\mathbf{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:

$$ \mathcal{H}_{\text{mag}} = -\sum_{i \neq j}^{N} J(r_{ij}) [\mathbf{s}_i \cdot \mathbf{s}_j - 1] \\ -\sum_{i \neq j}^{N} K(r_{ij}) [(\mathbf{s}_i \cdot \mathbf{s}_j)^2 - 1], \quad (3) $$

where $\mathbf{s}_i$ and $\mathbf{s}_j$ are classical atomic spins of unit length located on atoms i and j, $J(r_{ij})$ and $K(r_{ij})$ (in eV) are magnetic exchange functions, and $r_{ij}$ 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_{ij})$ and $K(r_{ij})$, the interatomic dependence is taken into account through the following function based on an approximation of the Bethe-Slater curve$^{75,76}$:

$$ f(r) = 4\alpha \left(\frac{r}{\delta}\right)^2 \left(1 - \gamma \left(\frac{r}{\delta}\right)^2\right) \exp\left(-\left(\frac{r}{\delta}\right)^2\right) \Theta(R_c - r), \quad (4) $$

where $\alpha$ denotes the interaction energy, $\delta$ the interaction decay length, $\gamma$ a dimensionless curvature parameter, r a distance, and $\Theta(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 ($\mathbf{\omega}_i$), mechanical forces ($\mathbf{F}_i$), and their corresponding virial components ($W(\mathbf{r}^N)$) are derived:

$$ \mathbf{\omega}_i = \frac{1}{\hbar} \sum_j^{N_i} J(r_{ij}) \mathbf{s}_j + K(r_{ij}) (\mathbf{s}_i \cdot \mathbf{s}_j) \mathbf{s}_j, \quad (5) $$

$$ \mathbf{F}_i = \sum_j^{N_i} \frac{dJ(r_{ij})}{d_{rij}} [\mathbf{s}_i \cdot \mathbf{s}_j - 1] \mathbf{e}_{ij} \\ + \frac{dK(r_{ij})}{d_{rij}} [(\mathbf{s}_i \cdot \mathbf{s}_j)^2 - 1] \mathbf{e}_{ij}, \quad (6) $$

$$ W(\mathbf{r}^N) = \sum_{i=1}^{N} \mathbf{r}_i \cdot \mathbf{F}_i, \quad (7) $$

where $\mathbf{r}^N$ denotes a 3N size vector of all atomic positions and $\mathbf{r}_i$ the position vector of atom i. We note that the virial
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Comparison of spin-spiral results along sections of the $\Gamma$H and $\Gamma$P high-symmetry lines. The upper plot displays the per-atom energy, the middle one the atomic moment fluctuations (in Bohr magneton per atom), and on the bottom the evolution of the pressure. The energy and pressure fluctuations are plotted with respect to the magnetic ground state at the $\Gamma$ point. The green and red dots represent experimental measurements obtained by Loong *et al.*$^{79}$ and Lynn$^{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 $\mu_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$-vectors on which our spin Hamiltonian was parametrized.</figcaption></figure>

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, $\Gamma$H and $\Gamma$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–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 $\Gamma$ point. We thus decided to parameterize our model only on spin-spirals corresponding to $q$-vectors for which the spin norm deviates from the ferromagnetic value ($\approx 2.2 \mu_B$/atom at the $\Gamma$ point) by less than 5%. The red dashed lines in Figure 5 delimit this $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_{ij})$ and $K(r_{ij})$ 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}$:

$$H_{\text{cubic}} = -\sum_{i=1}^{N} K_1 [(\mathbf{s}_i \cdot \hat{\mathbf{x}})^2 (\mathbf{s}_i \cdot \hat{\mathbf{y}})^2 + (\mathbf{s}_i \cdot \hat{\mathbf{y}})^2 (\mathbf{s}_i \cdot \hat{\mathbf{z}})^2 + \dots \\ (\mathbf{s}_i \cdot \hat{\mathbf{x}})^2 (\mathbf{s}_i \cdot \hat{\mathbf{z}})^2] + K_2^{(c)} (\mathbf{s}_i \cdot \hat{\mathbf{x}})^2 (\mathbf{s}_i \cdot \hat{\mathbf{y}})^2 (\mathbf{s}_i \cdot \hat{\mathbf{z}})^2, \quad (8)$$

with $K_1 = 0.001$ eV and $K_2^{(c)} = 0.0005$ eV the intensity coefficients corresponding to $\alpha$-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 $\alpha$- (BCC), $\varepsilon$- (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⁴⁵ to isolate the non-magnetic 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⁸⁷ or electrostatic interaction term⁸⁸ 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⁴⁴. 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⁹⁰ 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²³,⁴⁰, and set by the spin-lattice Hamiltonian below:

$$ \mathcal{H}_{sl}(\mathbf{r}, \mathbf{p}, \mathbf{s}) = \mathcal{H}_{mag}(\mathbf{r}, \mathbf{s}) + \sum_{i=1}^{N} \frac{|\mathbf{p}|^2}{2m_i} + \sum_{i,j=1}^{N} V_{SNAP}(r_{ij}) \quad (9) $$

where $\mathcal{H}_{mag}$ is the spin Hamiltonian defined by the combination of Eq. (3) and Eq. (8). The term $V_{SNAP}(r_{ij})$ 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 $\mathbf{p}$ and the mass of particle $i$ is $m_i$. Based on this spin-lattice Hamiltonian and leveraging the generalized spin-lattice Poisson bracket as defined by Yang et al.⁷⁸, the equations of motion can be defined as:

$$ \frac{d\mathbf{r}_i}{dt} = \frac{\mathbf{p}_i}{m_i} \quad (10) $$

$$ \frac{d\mathbf{p}_i}{dt} = \sum_{j, i \neq j}^{N} \left[ -\frac{dV_{SNAP}(r_{ij})}{dr_{ij}} + \frac{dJ(r_{ij})}{dr_{ij}}(\mathbf{s}_i \cdot \mathbf{s}_j) + \dots \right] \left( \mathbf{e}_{ij} - \frac{\gamma_L}{m_i} \mathbf{p}_i + \mathbf{f}(t) \right) \quad (11) $$

$$ \frac{d\mathbf{s}_i}{dt} = \frac{1}{(1 + \lambda^2)} \left[ (\mathbf{\omega}_i + \mathbf{\eta}(t)) \times \mathbf{s}_i + \dots \right. \\ \left. \lambda \mathbf{s}_i \times (\mathbf{\omega}_i \times \mathbf{s}_i) \right] \quad (12) $$

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 $\gamma_L$ is the Langevin damping constant for the lattice and $\mathbf{f}$ is a fluctuating force following Gaussian statistics given below⁴⁰.

$$ \langle \mathbf{f}(t) \rangle = 0 \quad (13) $$

$$ \langle f_{\alpha}(t) f_{\beta}(t') \rangle = 2k_{B}T_{l}\gamma_{L}\delta_{\alpha\beta}\delta(t-t') \quad (14) $$

The fluctuating force $\mathbf{f}$ is coupled to $\gamma_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 $\alpha$ and $\beta$ 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), $\lambda$ is the transverse damping constant and $\mathbf{\omega}_i$ is a spin force analogue as shown in Eq. (5). The variable $\mathbf{\eta}(t)$ is a random vector whose components are drawn from a Gaussian probability distribution given below:

$$ \langle \mathbf{\eta}(t) \rangle = 0 \quad (15) $$

$$ \langle \eta_{\alpha}(t) \eta_{\beta}(t') \rangle = D_{S} \delta_{\alpha\beta} \delta(t-t') \quad (16) $$

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\pi\lambda k_B T_s / \hbar^3$.

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
<table><caption>**Table 1.** Training set for linear SNAP model adapted from Ref. [89] to include explicit spin degrees of freedom. Regression of SNAP coefficients takes into account both configuration energies and forces from DFT, optimization of group weights is applied to either term independently. Weighted linear regression is carried out via reported optimal fit weights, values have already been scaled by the number of training points each group contributes. The two last columns report the obtained mean-absolute errors (MAEs) in eV per atom.</caption><thead><tr><th></th><th># of Config.</th><th># of Forces</th><th>Target Property</th><th>Energy Fit Weight</th><th>Forces Fit Weight</th><th>Energy MAE (eV)</th><th>Forces MAE (eV·Å<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Eq. of State</td><td>403</td><td>65286</td><td>Volumetric Deform</td><td>4.2·10<sup>3</sup></td><td>2.0·10<sup>5</sup></td><td>1.6·10<sup>-2</sup></td><td>2.4·10<sup>-1</sup></td></tr><tr><td>DFT-MD, 300K</td><td>40</td><td>15360</td><td>Bulk phonons</td><td>2.9·10<sup>5</sup></td><td>1.1·10<sup>5</sup></td><td>5.2·10<sup>-4</sup></td><td>2.4·10<sup>-1</sup></td></tr><tr><td>Liquid w/ Spins</td><td>10</td><td>3000</td><td>Magnetic Disorder</td><td>5.5·10<sup>1</sup></td><td>1.9·10<sup>4</sup></td><td>2.0·10<sup>-1</sup></td><td>5.9·10<sup>-1</sup></td></tr><tr><td>Liquid w/o Spins</td><td>52</td><td>15300</td><td>Structural Disorder</td><td>3.3·10<sup>3</sup></td><td>2.0·10<sup>4</sup></td><td>2.2·10<sup>-1</sup></td><td>8.0·10<sup>-1</sup></td></tr><tr><td>Point Defects</td><td>10</td><td>3096</td><td>Defect Energetics</td><td>1.4·10<sup>2</sup></td><td>3.5·10<sup>4</sup></td><td>2.8·10<sup>-2</sup></td><td>1.1·10<sup>-1</sup></td></tr><tr><td>Martensitic Transform</td><td>168</td><td>1008</td><td>α → ε</td><td>4.0·10<sup>2</sup></td><td>2.3·10<sup>3</sup></td><td>9.2·10<sup>-2</sup></td><td>2.3·10<sup>-1</sup></td></tr></tbody></table>

<table><caption>**Table 2.** Objective functions of the DAKOTA optimization with ground truth values taken from the present DFT calculations (at zero Kelvin) or experiments⁹⁰. Percent error is used as the objective function to avoid artificial importance scaling based on units of the target property.</caption><thead><tr><th></th><th>SNAP</th><th>Exp/DFT</th><th>Units</th><th>Error %</th></tr></thead><tbody><tr><td>c<sub>11</sub></td><td>243.25</td><td>239.55</td><td>GPa</td><td>1.54%</td></tr><tr><td>c<sub>12</sub></td><td>135.65</td><td>138.1</td><td>GPa</td><td>1.77%</td></tr><tr><td>c<sub>44</sub></td><td>118.73</td><td>120.75</td><td>GPa</td><td>1.67%</td></tr><tr><td>Bulk modulus</td><td>171.52</td><td>169.55</td><td>GPa</td><td>1.16%</td></tr><tr><td>0.5(c<sub>11</sub>-c<sub>12</sub>)</td><td>53.8</td><td>51.9</td><td>GPa</td><td>3.66%</td></tr><tr><td>Poisson ratio</td><td>0.358</td><td>0.36</td><td>-</td><td>1.10%</td></tr><tr><td>bcc energy</td><td>-8.25</td><td>-8.26</td><td>eV</td><td>0.02%</td></tr><tr><td>bcc lat. const.</td><td>2.838</td><td>2.83</td><td>Å</td><td>0.30%</td></tr></tbody></table>

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).

<figure><img src="image_10.png" /><figcaption>**Figure 6.** Symbols (○) and (□) associated to the left axis represent the pressure evolution in the FVC and PCMCC (similar to the PCC), as a function of lattice temperature. Dotted lines associated to the right axis represent the spin thermostat temperature for the PCMCC and FVC (similar to PCC) as a function of the lattice temperature.</figcaption></figure>

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-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}$:
$$T_{s,post}(T_l) = T_l - a_1 \qquad (17)$$

$$T_{s,pre}(T_l) = a_2 \exp \left[ - \left( \frac{T_l - a_3}{a_4} \right)^2 \right] - a_5 \qquad (18)$$

$$f_{sw}(T_l) = \frac{1}{2} \left[ 1 + \tanh \left( \frac{T_l - a_6}{a_7} \right) \right] \qquad (19)$$

$$T_s(T_l) = f_{sw} T_{s,post} + (1 - f_{sw}) T_{s,pre} \qquad (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.

## References

1. Tatsumoto, E. & Okamoto, T. Temperature dependence of the magnetostriction constants in iron and silicon iron. *J. Phys. Soc. Jpn.* **14**, 1588–1594 (1959).

2. Bahl, C. R. H. & Nielsen, K. K. The effect of demagnetization on the magnetocaloric properties of gadolinium. *J. Appl. Phys.* **105**, 013916 (2009).

3. Tavares, S., Fruchart, D., Miraglia, S. & Laborie, D. Magnetic properties of an aisi 420 martensitic stainless steel. *J. alloys compounds* **312**, 307–314 (2000).

4. Huang, S., Holmström, E., Eriksson, O. & Vitos, L. Mapping the magnetic transition temperatures for medium- and high-entropy alloys. *Intermetallics* **95**, 80–84 (2018).

5. Rao, Z. *et al*. Unveiling the mechanism of abnormal magnetic behavior of fenicomcu high-entropy alloys through a joint experimental-theoretical study. *Phys. Rev. Mater.* **4**, 014402 (2020).

6. Jaime, M. *et al*. Piezomagnetism and magnetoelastic memory in uranium dioxide. *Nat. communications* **8**, 1–7 (2017).

7. Nussle, T., Thibaudeau, P. & Nicolis, S. Dynamic magnetostriction for antiferromagnets. *Phys. Rev. B* **100**, 214428 (2019).

8. Lejman, M. *et al*. Magnetoelastic and magnetoelectric couplings across the antiferromagnetic transition in multiferroic bifeo 3. *Phys. Rev. B* **99**, 104103 (2019).

9. Patrick, C. E., Marchant, G. A. & Staunton, J. B. Spin orientation and magnetostriction of tb 1- x dy x fe 2 from first principles. *Phys. Rev. Appl.* **14**, 014091 (2020).

10. Graham, R., Morosin, B., Venturini, E. & Carr, M. Materials Modification and Synthesis Under High Pressure Shock Compression. *Annu. Rev. Mater. Sci.* **16**, 315–341, DOI: 10.1146/annurev.ms.16.080186.001531 (1986).

11. Surh, M. P., Benedict, L. X. & Sadigh, B. Magnetostructural transition kinetics in shocked iron. *Phys. review letters* **117**, 085701 (2016).

12. Moses, E. I., Boyd, R. N., Remington, B. A., Keane, C. J. & Al-Ayat, R. The national ignition facility: Ushering in a new age for high energy density science. *Phys. Plasmas* **16**, 041006, DOI: 10.1063/1.3116505 (2009).

13. Tschentscher, T. *et al*. Photon beam transport and scientific instruments at the european xfel. *Appl. Sci.* **7**, DOI: 10.3390/app7060592 (2017).

14. Tan, X., Chan, S., Han, K. & Xu, H. Combined effects of magnetic interaction and domain wall pinning on the coercivity in a bulk nd 60 fe 30 al 10 ferromagnet. *Sci. reports* **4**, 6805 (2014).

15. Gràcia-Condal, A. *et al*. Multicaloric effects in metamagnetic heusler ni-mn-in under uniaxial stress and magnetic field. *Appl. Phys. Rev.* **7**, 041406 (2020).

16. Wallace, D. C., Sidles, P. & Danielson, G. Specific heat of high purity iron by a pulse heating method. *J. applied physics* **31**, 168–176 (1960).

17. Touloukian, Y. & Buyco, E. Thermophysical properties of matter, vol. 4, specific heat. *IFI/Plenum, New York* (1970).

18. Chandler, D. *Introduction to modern statistical mechanics* (1987).

19. Horstemeyer, M. F. The near Future: ICME for the Creation of New Materials and Structures. In *Integrated Computational Materials Engineering (ICME) for Metals*, chap. 10, 410–423, DOI: 10.1002/9781118342664.ch10 (John Wiley & Sons, Ltd, 2012).

20. van der Giessen, E. *et al*. Roadmap on multiscale materials modeling. *Model. Simul. Mater. Sci. Eng.* **28**, 043001 (2020).

21. Alder, B. J. & Wainwright, T. E. Studies in Molecular Dynamics. I. General Method. *The J. Chem. Phys.* **31**, 459–466, DOI: 10.1063/1.1730376 (1959).

22. Rapaport, D. C. *The art of molecular dynamics simulation* (Cambridge university press, 2004).

23. Plimpton, S. Fast parallel algorithms for short-range molecular dynamics. *J. computational physics* **117**, 1–19 (1995).
24. Voter, A. F., Montalenti, F. & Germann, T. C. Extending the time scale in atomistic simulation of materials. *Annu. Rev. Mater. Res.* **32**, 321–346 (2002).

25. Zepeda-Ruiz, L. A., Stukowski, A., Oppelstrup, T. & Bulatov, V. V. Probing the limits of metal plasticity with molecular dynamics simulations. *Nature* **550**, 492–495 (2017).

26. Huan, T. D. et al. A universal strategy for the creation of machine learning-based atomistic force fields. *npj Comput. Mater.* **3**, DOI: 10.1038/s41524-017-0042-y (2017).

27. Smith, J. S., Isayev, O. & Roitberg, A. E. Ani-1: an extensible neural network potential with dft accuracy at force field computational cost. *Chem. Sci.* **8**, 3192–3203, DOI: 10.1039/C6SC05720A (2017).

28. Zhang, L., Han, J., Wang, H., Car, R. & E. W. Deep Potential Molecular Dynamics: A Scalable Model with the Accuracy of Quantum Mechanics. *Phys. Rev. Lett.* **120**, 143001, DOI: 10.1103/PhysRevLett.120.143001 (2018).

29. Bartók, A. P., Payne, M. C., Kondor, R. & Csányi, G. Gaussian Approximation Potentials: The Accuracy of Quantum Mechanics, without the Electrons. *Phys. Rev. Lett.* **104**, 136403, DOI: 10.1103/PhysRevLett.104.136403 (2010).

30. Jaramillo-Botero, A., Naserifar, S. & Goddard, W. A. General Multiobjective Force Field Optimization Framework, with Application to Reactive Force Fields for Silicon Carbide. *J. Chem. Theory Comput.* **10**, 1426–1439, DOI: 10.1021/ct5001044 (2014).

31. Lubbers, N., Smith, J. S. & Barros, K. Hierarchical modeling of molecular energies using a deep neural network. *The J. Chem. Phys.* **148**, 241715, DOI: 10.1063/1.5011181 (2018).

32. Thompson, A. P., Swiler, L. P., Trott, C. R., Foiles, S. M. & Tucker, G. J. Spectral neighbor analysis method for automated generation of quantum-accurate interatomic potentials. *J. Comput. Phys.* **285**, 316–330, DOI: 10.1016/j.jcp.2014.12.018 (2015).

33. Kohn, W. & Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. *Phys. Rev.* **140**, A1133–A1138, DOI: 10.1103/PhysRev.140.A1133 (1965).

34. Li, X.-G., Chen, C., Zheng, H., Zuo, Y. & Ong, S. P. Complex strengthening mechanisms in the nbmotaw multiprincipal element alloy. *npj Comput. Mater.* **6**, 1–10 (2020).

35. Cusentino, M., Wood, M. & Thompson, A. Suppression of helium bubble nucleation in beryllium exposed tungsten surfaces. *Nucl. Fusion* **60**, 126018 (2020).

36. Dragoni, D., Daff, T. D., Csányi, G. & Marzari, N. Achieving dft accuracy with a machine-learning interatomic potential: Thermomechanics and defects in bcc ferromagnetic iron. *Phys. Rev. Mater.* **2**, 013808 (2018).

37. Ma, P.-W., Woo, C. & Dudarev, S. Large-scale simulation of the spin-lattice dynamics in ferromagnetic iron. *Phys. review B* **78**, 024434 (2008).

38. Ma, P.-W., Dudarev, S. & Woo, C. Spilady: A parallel cpu and gpu code for spin-lattice magnetic molecular dynamics simulations. *Comput. Phys. Commun.* **207**, 350–361 (2016).

39. Ma, P.-W. & Dudarev, S. Atomistic spin-lattice dynamics. *Handb. Mater. Model. Methods: Theory Model.* 1017–1035 (2020).

40. Tranchida, J., Plimpton, S., Thibaudeau, P. & Thompson, A. P. Massively parallel symplectic algorithm for coupled magnetic spin dynamics and molecular dynamics. *J. Comput. Phys.* **372**, 406–425 (2018).

41. Dos Santos, G. et al. Size-and temperature-dependent magnetization of iron nanoclusters. *Phys. Rev. B* **102**, 184426 (2020).

42. Zhou, Y., Tranchida, J., Ge, Y., Murthy, J. & Fisher, T. S. Atomistic simulation of phonon and magnon thermal transport across the ferromagnetic-paramagnetic transition. *Phys. Rev. B* **101**, 224303 (2020).

43. Ma, P.-W., Dudarev, S. & Wróbel, J. S. Dynamic simulation of structural phase transitions in magnetic iron. *Phys. Rev. B* **96**, 094418 (2017).

44. Eldred, M. S. et al. Dakota, a multilevel parallel object-oriented framework for design optimization, parameter estimation, uncertainty quantification, and sensitivity analysis. Tech. Rep., Citeseer (2006).

45. Thompson, A. P., Swiler, L. P., Trott, C. R., Foiles, S. M. & Tucker, G. J. Spectral neighbor analysis method for automated generation of quantum-accurate interatomic potentials. *J. Comput. Phys.* **285**, 316–330 (2015).

46. Crangle, J. & Goodman, G. The magnetization of pure iron and nickel. *Proc. Royal Soc. London. A. Math. Phys. Sci.* **321**, 477–491 (1971).

47. Seki, I. & Nagata, K. Lattice constant of iron and austenite including its supersaturation phase of carbon. *ISIJ international* **45**, 1789–1794 (2005).

48. Basinski, Z. S., Hume-Rothery, W. & Sutton, A. The lattice expansion of iron. *Proc. Royal Soc. London. Ser. A. Math. Phys. Sci.* **229**, 459–467 (1955).

49. Ashcroft, N. W., Mermin, N. D. & Wei, D. Solid State Physics (Cengage Learning Asia Pte Limited, 2016).

50. Dragoni, D., Ceresoli, D. & Marzari, N. Thermoelastic properties of α-iron from first-principles. *Phys. Rev. B* **91**, 104105 (2015).

51. Dragoni, D., Ceresoli, D. & Marzari, N. Vibrational and thermoelastic properties of bcc iron from selected eam potentials. *Comput. Mater. Sci.* **152**, 99–106 (2018).
52. Körmann, F. et al. Free energy of bcc iron: Integrated ab initio derivation of vibrational, electronic, and magnetic contributions. *Phys. Rev. B* **78**, 033102 (2008).

53. Arakawa, K. et al. Quantum de-trapping and transport of heavy defects in tungsten. *Nat. materials* **19**, 508–511 (2020).

54. Drautz, R. & Fähnle, M. Spin-cluster expansion: Parametrization of the general adiabatic magnetic energy surface with ab initio accuracy. *Phys. Rev. B* **69**, 104404 (2004).

55. Drautz, R. Atomic cluster expansion of scalar, vectorial, and tensorial properties including magnetism and charge transfer. *Phys. Rev. B* **102**, 024104 (2020).

56. Marinica, M.-C., Willaime, F. & Crocombette, J.-P. Irradiation-induced formation of nanocrystallites with c 15 laves phase structure in bcc iron. *Phys. review letters* **108**, 025501 (2012).

57. Chapman, J. B., Ma, P.-W. & Dudarev, S. L. Effect of non-heisenberg magnetic interactions on defects in ferromagnetic iron. *Phys. Rev. B* **102**, 224106 (2020).

58. Chapman, J. B., Ma, P.-W. & Dudarev, S. L. Dynamics of magnetism in fe-cr alloys with cr clustering. *Phys. Rev. B* **99**, 184413 (2019).

59. Klaver, T., Drautz, R. & Finnis, M. Magnetism and thermodynamics of defect-free fe-cr alloys. *Phys. Rev. B* **74**, 094435 (2006).

60. Kalantar, D. et al. Direct observation of the α-ε transition in shock-compressed iron via nanosecond x-ray diffraction. *Phys. review letters* **95**, 075502 (2005).

61. Woo, C., Wen, H., Semenov, A., Dudarev, S. & Ma, P.-W. Quantum heat bath for spin-lattice dynamics. *Phys. Rev. B* **91**, 104306 (2015).

62. Bergqvist, L. & Bergman, A. Realistic finite temperature simulations of magnetic systems using quantum statistics. *Phys. Rev. Mater.* **2**, 013802 (2018).

63. Novikov, I., Grabowski, B., Kormann, F. & Shapeev, A. Machine-learning interatomic potentials reproduce vibrational and magnetic degrees of freedom. *arXiv preprint arXiv:2012.12763* (2020).

64. Kresse, G. & Furthmüller, J. Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. *Comput. Mater. Sci.* **6**, 15–50, DOI: https://doi.org/10.1016/0927-0256(96)00008-0 (1996).

65. Kresse, G. & Joubert, D. From ultrasoft pseudopotentials to the projector augmented-wave method. *Phys. Rev. B* **59**, 1758–1775, DOI: 10.1103/PhysRevB.59.1758 (1999).

66. Perdew, J. P., Burke, K. & Ernzerhof, M. Generalized gradient approximation made simple. *Phys. review letters* **77**, 3865 (1996).

67. Blöchl, P. E. Projector augmented-wave method. *Phys. review B* **50**, 17953 (1994).

68. Zimmermann, B. et al. Comparison of first-principles methods to extract magnetic parameters in ultrathin films: Co/pt (111). *Phys. Rev. B* **99**, 214426 (2019).

69. Halilov, S., Perlov, A., Oppeneer, P. & Eschrig, H. Magnon spectrum and related finite-temperature magnetic properties: A first-principle approach. *EPL (Europhysics Lett.* **39**, 91 (1997).

70. Kurz, P., Förster, F., Nordström, L., Bihlmayer, G. & Blügel, S. Ab initio treatment of noncollinear magnets with the full-potential linearized augmented plane wave method. *Phys. Rev. B* **69**, 024415 (2004).

71. Sandratskii, L. Noncollinear magnetism in itinerant-electron systems: theory and applications. *Adv. Phys.* **47**, 91–160 (1998).

72. Marsman, M. & Hafner, J. Broken symmetries in the crystalline and magnetic structures of γ-iron. *Phys. Rev. B* **66**, 224409 (2002).

73. Rosengaard, N. & Johansson, B. Finite-temperature study of itinerant ferromagnetism in fe, co, and ni. *Phys. Rev. B* **55**, 14975 (1997).

74. Szilva, A. et al. Interatomic exchange interactions for finite-temperature magnetism and nonequilibrium spin dynamics. *Phys. review letters* **111**, 127204 (2013).

75. Kaneyoshi, T. *Introduction to amorphous magnets* (World Scientific Publishing Company, 1992).

76. Yosida, K., Mattis, D. C. & Yosida, K. *THEORY OF MAGNETISM.: Edition en anglais*, vol. 122 (Springer Science & Business Media, 1996).

77. Pajda, M., Kudrnovský, J., Turek, I., Drchal, V. & Bruno, P. Ab initio calculations of exchange interactions, spin-wave stiffness constants, and curie temperatures of fe, co, and ni. *Phys. Rev. B* **64**, 174402 (2001).

78. Yang, K.-H. & Hirschfelder, J. O. Generalizations of classical poisson brackets to include spin. *Phys. Rev. A* **22**, 1814 (1980).

79. Loong, C.-K., Carpenter, J., Lynn, J., Robinson, R. & Mook, H. Neutron scattering study of the magnetic excitations in ferromagnetic iron at high energy transfers. *J. applied physics* **55**, 1895–1897 (1984).

80. Lynn, J. Temperature dependence of the magnetic excitations in iron. *Phys. Rev. B* **11**, 2624 (1975).

81. Leger, J., Loriers-Susse, C. & Vodar, B. Pressure effect on the curie temperatures of transition metals and alloys. *Phys. Rev. B* **6**, 4250 (1972).

82. Morán, S., Ederer, C. & Fähnle, M. Ab initio electron theory for magnetism in fe: Pressure dependence of spin-wave energies, exchange parameters, and curie temperature. *Phys. Rev. B* **67**, 012407 (2003).

83. Körmann, F., Dick, A., Hickel, T. & Neugebauer, J. Pressure dependence of the curie temperature in bcc iron
studied by ab initio simulations. *Phys. Rev. B* **79**, 184406 (2009).

84. Skomski, R. et al. *Simple models of magnetism* (Oxford University Press on Demand, 2008).

85. Wood, M. A. & Thompson, A. P. Extending the accuracy of the snap interatomic potential form. *The J. Chem. Phys.* **148**, 241721 (2018).

86. Zuo, Y. et al. Performance and cost assessment of machine learning interatomic potentials. *The J. Phys. Chem. A* **124**, 731–745 (2020).

87. Wood, M. A., Cusentino, M. A., Wirth, B. D. & Thompson, A. P. Data-driven material models for atomistic simulation. *Phys. Rev. B* **99**, 184305 (2019).

88. Deng, Z., Chen, C., Li, X.-G. & Ong, S. P. An electrostatic spectral neighbor analysis potential for lithium nitride. *npj Comput. Mater.* **5**, 1–8 (2019).

89. Goryaeva, A. M., Maillet, J.-B. & Marinica, M.-C. Towards better efficiency of interatomic linear machine learning potentials. *Comput. Mater. Sci.* **166**, 200–209 (2019).

90. Adams, J. J., Agosta, D., Leisure, R. & Ledbetter, H. Elastic constants of monocrystal iron from 3 to 500 k. *J. applied physics* **100**, 113530 (2006).

## Acknowledgements

All authors thank Mark Wilson for his detailed review and edits. Sandia National Laboratories is a multimission laboratory managed and operated by National Technology & Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International Inc., for the U.S. Department of Energy's National Nuclear Security Administration under contract DE-NA0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government. AC acknowledges funding from the Center for Advanced Systems Understanding (CASUS) which is financed by the German Federal Ministry of Education and Research (BMBF) and by the Saxon State Ministry for Science, Art, and Tourism (SMWK) with tax funds on the basis of the budget approved by the Saxon State Parliament.

## Author contributions statement

AC, MAW, MPD and JT performed the DFT calculations. JBM, MCM, JT and MAW generated the Database of configurations. JT implemented the extended spin Hamiltonian and the magnetic pressure computation in LAMMPS, and parametrized it on *first-principles* calculations. MAW and SN trained the SNAP potential. JT and SN performed the magneto-static calculations. SN, APT, MAW and JT performed the magneto-dynamics calculations. All authors participated in conceiving the research and writing the manuscript.

## Competing interests

The authors declare no competing interests.