High-Dimensional Neural Network Potentials for Magnetic Systems Using Spin-Dependent Atom-Centered Symmetry Functions

Machine learning potentials have emerged as a powerful tool to extend the time and length scales of first principles-quality simulations. Still, most machine learning potentials cannot distinguish different electronic spin orientations and thus are not applicable to materials in different magnetic states. Here, we propose spin-dependent atom-centered symmetry functions as a new type of descriptor taking the atomic spin degrees of freedom into account. When used as input for a high-dimensional neural network potential (HDNNP), accurate potential energy surfaces of multicomponent systems describing multiple magnetic states can be constructed. We demonstrate the performance of these magnetic HDNNPs for the case of manganese oxide, MnO. We show that the method predicts the magnetically distorted rhombohedral structure in excellent agreement with density functional theory and experiment. Its efficiency allows to determine the N\'{e}el temperature considering structural fluctuations, entropic effects, and defects. The method is general and is expected to be useful also for other types of systems like oligonuclear transition metal complexes.

In recent years, machine learning potentials (MLP), which allow to extend the time and length scales of first principles-quality atomistic simulations [1][2][3], have become a rapidly growing field of research. More and more complex systems have been investigated driving new developments and extending the applicability of MLPs. Many of the current MLPs can be classified into four generations [4,5]: The first generation of MLPs proposed already in 1995 [6] typically employs a single or a few feed-forward neural networks making them applicable to low-dimensional systems like small molecules in vacuum or diatomic molecules interacting with frozen surfaces of solids [7,8]. In 2007, second-generation MLPs have become available with the introduction of high-dimensional neural network potentials (HDNNP) [9][10][11][12]. These MLPs are applicable to systems containing thousands of atoms by making use of the locality of a large part of the atomic interactions. For this purpose, the potential energy is calculated as a sum of local environment-dependent atomic energies defined by a cutoff radius. Most modern MLPs belong to this second generation, for example, HDNNPs [9,13], Gaussian approximation potentials [14], moment tensor potentials [15], and many others [16][17][18]. The third generation of MLPs includes long-range interactions, mainly electrostatics but also dispersion, beyond the cutoff radius. The electrostatic interactions can be based on element-specific fixed charges [14,19] or local environment-dependent atomic charges expressed by machine learning [20][21][22]. Also message passing networks with explicit electrostatics have been proposed [23]. Finally, fourth-generation MLPs take non-local or even global dependencies in the electronic structure into account, and consequently the atomic charges can adapt to non-local charge transfer and even different global charge states. The first method of this generation has been the charge equilibration neural network technique (CENT) [24], and further methods like Becke population neural networks (BpopNN) [25] and fourth-generation (4G) HDNNPs [26] have been introduced recently.
In spite of these methodical advances, which extended the complexity of the systems that now can be studied and the physical phenomena that can be included, a remaining limitation of MLPs is the inability to take different spin orientations and thus magnetic interactions into account. The reason for this limitation is that typically MLPs are trained to represent the potential energy surface of one electronic state as a function of structural descriptors -most often but not necessarily the ground state of a system. With the exception of different global charge states in fourth-generation potentials, describing multiple electronic states usually requires to train separate MLPs for each state [27][28][29][30][31] or to use one MLP yielding a vector of output energies corresponding to the different excited states [32][33][34][35][36]. The unique structureenergy relation of a given state is a central component of nowadays MLPs, and energy changes resulting from different spin orientations give rise to the problem of contradictory information in the training process. Only very recently, magnetic moment tensor potentials (mMTP) [37] have been proposed as a first example of a MLP which is able to describe magnetic systems containing a single element like iron.
For studying magnetic materials, the description of one magnetic state only is not sufficient as the atomic magnetic moments fluctuate and change signs at finite temperatures. These atomic spin-flips often play a role already at ambient temperatures as the energy differences between different spin configurations are typically only in the order of a few meV atom −1 . For example, the Curie and Néel temperatures, at which a material loses its ferromagnetic or antiferromagnetic state, are typically below 1000 and 500 K, respectively [38]. Magnetic transitions often give rise to structural changes [39]. However, current implementations of MLPs are not able to capture these effects because the magnetic interaction is not included explicitly. For example, if an antiferromagnetic ground state is used for training the MLP, this state will also be the basis of atomistic simulations at higher temperatures irrespective of the true magnetic ground state under these conditions.
To calculate the energy contribution of a magnetic configuration in an atomistic simulation, models like the Ising model [40], the Heisenberg model [41], and the Hubbard model [42] are widely used. However, these models are based on lattices, and structural and spin changes at finite temperatures cannot be considered simultaneously. If both contributions are important, the only generally applicable option is the use of ab initio molecular dynamics, which is based on the explicit calculation of the electronic structure in each step and thus is able to distinguish different magnetic states. This approach, however, is inevitably associated with high computational costs restricting simulations to small systems and short time scales.
Beyond the field of potential-energy surfaces, several machine learning approaches have been developed to address the properties of magnetic compounds, for example for the prediction of magnetic moments [38] and ordering temperatures [43][44][45] from structural parameters. Moreover, machine learning methods are able to classify ferromagnetic and antiferromagnetic ground state materials [46] and to predict spin state splittings and metal-ligand bond distances in transition metal complexes [47,48]. Recently, a high-dimensional neural network approach for the prediction of oxidation and spin states has been developed [49].
In spite of these applications of machine learning to magnetic compounds, with the exception of mMTPs [37] to date there is no method enabling large-scale first principles-quality atomistic simulations of systems explicitly including magnetic interactions. The main reason for this lack of methods is the use of descriptors in current MLPs which exclusively depend on the atomic structure, like atom-centered symmetry functions (ACSF) [50], smooth overlap of atomic positions (SOAP) [51] and many others [52]. Exceptions are BpopNNs [25] and the recently introduced fourth generation of HDNNPs [26] in which the atomic charges are included as additional information besides the structural descriptors to provide qualitative information about the electronic structure.
Describing the spin configuration in a form suitable as input for MLPs is very challenging, as not only the absolute spins but also the relative orientations of the atomic spins and distances of the respective atoms are vital. Hence, to predict the energy and forces simultaneously as a function of the geometric and spin configuration, suitable spin-dependent descriptors are needed. Here, we propose such a descriptor to construct MLPs simultaneously applicable to multiple magnetic states based on a new type of spin-dependent atom-centered symmetry function (sACSF), which adds the description of the magnetic configuration for the case of collinear spin polarization. The sACSFs produce atomic energies as a function of the local geometric and spin environment and thus formally represent a second-generation potential that can also be combined with third and fourth generation MLPs to include additional physics like long-range electrostatic interactions. In this work we will benchmark sACSFs employing second-generation HDNNPs [9], extending these potentials by an explicit dependence on the full spin configuration space to describe magnetic interactions.
We choose manganese oxide, MnO, to assess the quality of the resulting magnetic HDNNPs (mHDNNP) that can be constructed using sACSFs, because of the well characterized antiferromagnetic ground state configuration [53,54], the Néel temperature of T exp N = 116 K [55,56], and the rhombohedral distortion of the antiferromagnetic phase with lattice constant a exp = 4.430Å and lattice angle α exp = 90.62 • at 8 K [57] making this system a very interesting and challenging benchmark case. The magnetic unit cell is a 2 × 2 × 2 supercell of the geometric unit cell and is built from (111) planes of ferromagnetically coupled Mn ions [54]. These planes couple antiferromagnetically to the neighboring planes. Spinpolarized density functional theory (DFT) calculations employing the hybrid functional PBE0 [58,59] yield the correct magnetic ground state, called AFM-II, with the rhombohedral lattice parameters a PBE0 = 4.40Å and α PBE0 = 90.88 • in good agreement with experiment [60]. The rhombohedral distortion is a consequence of the magnetic anisotropy, which arises from the magnetic dipole interactions [61].
Several questions have not yet been conclusively answered for this system because the required large simulation cells and long time scales are inaccessible by DFT and the determination of coupling constants for lattice models gets very complicated as soon as defects have to be considered. Using a mHDNNP we are able to perform high-throughput studies of the equilibrium geometries of different magnetic configurations as well as to simulate the antiferromagnetic to paramagnetic phase transition to determine the Néel temperature. Moreover, the mHDNNP enables the inclusion of defects, as Mn vacancies, to simulate their role in real materials.

Magnetic High-Dimensional Neural Network Potential
The reference data set of the mHDNNP consists of 3101 2 × 2 × 2 bulk supercells of MnO and Mn 0.969 O in various magnetic states and their corresponding HSE06 DFT energies and force components. The supercells include different displacements of the atomic positions from the ideal rock salt lattice as well as distortions of the lattice parameters. The construction of the reference data set is described in detail in the Supplementary Information [62].
Training this data set with conventional ACSFs yields an energy root mean squared error (RMSE) of about 11 meV atom −1 in the best potential we were able to construct. This RMSE is an order of magnitude higher than the usual HDNNP accuracy of 1 meV atom −1 because the ACSFs can only provide geometrical information to assign the energy. Thus, HDNNPs based on ACSFs only predict an averaged potential energy regardless of the magnetic configurations as the best compromise while, for example, the HSE06 DFT functional yields an energy difference of 45.9 meV atom −1 between the AFM-II order and the ferromagnetic (FM) order for the ideal rock salt MnO structure using the experimental lattice constant a exp .
By including sACSFs to distinguish the magnetic configurations the energy RMSE is reduced by one order of magnitude to about 1 meV atom −1 . Now, the energy difference between the AFM-II order and the FM order for the rock salt structure is resolved and is predicted to be 46.3 meV atom −1 in excellent agreement with the DFT reference data. Both results demonstrate the ability of the sACSFs to accurately describe the different magnetic configurations. The energy errors ∆E = E mHDNNP − E DFT and force component errors ∆F = F mHDNNP − F DFT are plotted as a function of the reference values in Figure 1 Table I for the training and test data sets. It can be clearly seen that even when including various magnetic states the accuracy of the mHDNNP is in the typical region of state-of-the-art MLPs of 1 meV atom −1 and 0.1 eV a −1 0 for energies and forces with a test energy RMSE of 1.11 meV atom −1 and a test force components RMSE of 0.066 eV a −1 0 . a 0 is the Bohr radius. In comparison, the RMSE of a recent magnetic moment tensor potential for defect-free body centered cubic iron restricted to a fixed lattice parameter is 2.0 meV atom −1 [37], i.e., in the same order of magnitude.

Magnetic Configurations
Besides the accurate description of the energies and forces the sACSFs enable to predict structural changes arising from the magnetic configuration. For example, for the rhombohedral MnO global minimum structure with AFM-II order (see Figure 2   . Also excited magnetic configurations are predicted in agreement with DFT. For example, the optimized lattice parameter of the resulting cubic lattice for the FM configuration (see Figure 2   In conclusion, the lattices of FM (cubic), AFM-I (tetragonal), and AFM-II orders (rhombohedral) as well as other magnetic configurations can be different due to the magnetic interaction. This cannot be described by a lattice model such as the Heisenberg spin Hamiltonian. Further, the mHDNNP can provide information about the influence of defects in the magnetic order. For example, one spin-flip in the AFM-II configuration per 2 × 2 × 2 supercell reduces the rhombohedral distortion by 0.09 • (DFT: 0.11 • ). The corresponding energy increase is 1.5 meV atom −1 (DFT: 1.5 meV atom −1 ).
The efficiency of the mHDNNP in combination with a basin-hopping Monte Carlo scheme [63] in which Monte Carlo (MC) spin-flips are employed instead of atomic displacements enables high throughput searches of the minima in spin configuration space which would be computationally too demanding employing DFT directly. The spin-flip basin hopping Monte Carlo (SFBHMC) simu-lations of 2 × 2 × 2 MnO supercells as well as molecular dynamics (MD) simulations including MC spin-flips (MDMC) of 6 × 6 × 6 MnO supercells confirmed the rhombohedral AFM-II magnetic order to be the global minimum in agreement with experiments [54]. However, if the lattice is restricted to be cubic, two degenerate global minima exist: The AFM-II magnetic configuration shown in Figure 2 (a) and another antiferromagnetic order shown in Figure 2 (b). A recalculation of the structure using DFT confirms this observation. Radial distribution functions (see Supplementary Information [62]) show that the same ferro-and antiferromagnetic interactions are present in the two configurations. This second cubic global minimum configuration stays cubic in an unconstrained optimization. The rhombohedral distortion of the AFM-II order can therefore be identified as the origin of the energetic preference. The energy gain by the rhombohedral distortion is 1.9 meV atom −1 (DFT: 2.1 meV atom −1 ).

Magnetic Interactions
For the ideal cubic structure a Heisenberg spin Hamiltonian can be constructed, which includes the magnetic interactions of atom i with its nearest neighbors n i and next nearest neighbors m i . The strengths of the magnetic coupling between the vector spin operators S i and S ni as well as S i and S mi are given by the exchange coupling constants J 1 and J 2 , respectively. k B is the Boltzmann constant. The exchange coupling constants can be determined employing energetic differences among the FM, AFM-I, and AFM-II configurations, and Equation (1) for these systems can be rearranged to yield with the spin S = 5 2 of the high-spin Mn II ions. Using mean field theory [64] the Néel temperature can then be calculated as Employing the HSE06 DFT energies using the experimental lattice parameter a exp = 4.430Å we obtain J 1 = −13.9 K, J 2 = −14.5 K, and T N = 255 K. The mHDNNP results match these values almost perfectly with J 1 = −14.0 K, J 2 = −14.6 K, and T N = 256 K. This agreement again underlines the performance of the sACSFs as well as of the mHDNNP method to describe the multiple potential energy surfaces of the MnO magnetic states. While the experimental Néel temperature of T exp N = 116 K [55,56] is lower, this overestimation of the DFT value compared to experiment was also found in previous studies, for example, the PBE0 functional yields T N = 240 K [60] and the HSE03 functional T N = 230 K [65].

Néel Temperature
Mean field theory misses the influence of the specific magnetic configurations of MnO on the determination of the Néel temperature. Moreover, the underlying Heisenberg spin Hamiltonian restricts the system to a fixed lattice. Employing the mHDNNP we can overcome both limitations step by step to reveal their influences on the Néel temperature. While conventional MC spin-flip simulations use a fixed lattice but allow to explore the specific states of MnO, by including N pT molecular dynamics trajectories in MDMC simulations both the atomic positions as well as the lattice parameters can adapt to the magnetic configurations, and accounting for thermal fluctuations can provide a more realistic description.
To determine the temperature of a phase transition the molar heat capacity can be employed because it shows a peak at the transition temperature for systems of finite size. Employing N pT MD simulations the heat capacity at constant pressure C p can be obtained from the fluctuations of the total energy, with N atoms atoms in the simulation cell, the mean simulation temperature T , the total energy per atom E tot (n) as a function of the MD time step n for the total number of simulation steps N steps , and the mean total energy per atom during the simulation E tot . N A is the Avogadro constant.
From MC spin-flip simulations the influence of the magnetic degrees of freedom on the heat capacity at constant volume C V can be calculated, The contribution of the atomic motions to C V is taken into account by the term 3k B N A , as the atomic motions are not considered in the conventional MC spin-flip simulations. The energy fluctuation is calculated from the potential energies per atom E(n) as a function of the MC step n for the total number of steps N steps compared to the mean potential energy per atom E.
To identify the AFM-II to paramagnetic (PM) transition, the temperature dependence of an order parameter C can be used, longs to the antiferromagnetic to paramagnetic transition. The Néel temperature is the same for the disordering process (MC AFM-II cub ) as well as the ordering process (MC cub ), i.e., no hysteresis is present. The AFM-II to PM (MC AFM-II cub ) as well as the second cubic global minimum to PM (MC 2 nd cub ) transition temperatures are identical because of the equal interactions of both cubic global minimum configurations. The ordering process can consequently also end in the second cubic global minimum configuration as happened in the MC cub simulations at 120, 160, and 200 K. However, below the Néel temperature the interconversion of both cubic global minima is kinetically hindered. In summary, the inclusion of specific state information increases the Néel temperature of MnO by about 50 K compared to mean field theory.
The restriction to a cubic lattice is an approximation because on the one hand the AFM-II configuration prefers a rhombohedral lattice. Employing this lattice (MC AFM-II min ) increases the Néel temperature above 480 K. This overestimation is a consequence of the decreased AFM-II energy and increased PM energy. On the other hand, an optimized lattice of the initial PM configuration (MC min ) leads to a Néel temperature of about 280 K. In conclusion, the choice of a specific fixed lattice can change the Néel temperature at least by 200 K. Adapting the lattice as well as the atomic positions to the magnetic order is consequently of major importance to obtain reliable results.
N pT MD steps enable to sample the thermodynamic equilibrium of the given magnetic configuration at pressure p and temperature T to get rid of simulation arti-facts due to a restricted geometric structure. The MDMC simulations are therefore a much better representation of realistic, experimental conditions. MDMC simulations predict a AFM-II to PM transition at 300 K (see MDMC AFM-II and MDMC in Figure 3 and 4). This temperature is similar to the result of the cubic lattice due to a compensation of two main influences: The optimized AFM-II lattice (at p = 1 bar) leads to an energy gain of 1.9 meV atom −1 compared to the cubic lattice with a exp whereby the mean energy gain of the PM configurations is 1.3 meV atom −1 . The data of the PM phase were calculated from optimization results of 1000 PM configurations, which were obtained during the 10 ns MDMC simulation at 400 K. In addition, thermal fluctuations of the atomic positions are included in MDMC simulations leading to thermal expansion of the MnO lattice with increasing temperature. To quantify the energy increase due to the increased lattice volume, optimizations at various pressures of the AFM-II and PM configurations were performed. If the optimized volume at p = 1 bar is increased to the volume given at 300 K in the MDMC simulations (see Figure 5), the energy of the AFM-II configuration is increased by 1.7 meV atom −1 and the mean energy of the PM configurations by 1.1 meV atom −1 (see Supplementary Information [62]). The cube root of the mean lattice volume, i.e., a hypothetical averaged cubic lattice constant, shows a discontinuous increase of 0.005Å at the Néel temperature ( Figure 5). This increase is similar to the difference of the optimized AFM-II and PM values which differ by 0.005Å. Experimentally an increase of about 0.004Å has been observed at the Néel temperature [66]. The mean lattice angle decreases from 90.77 to 90.00 • at the Néel temperature as shown in the Supplementary Information [62] matching the experimental PM lattice angle of 90.00 • [66]. This discontinuous change of the lattice volume and shape confirms the assignment to a first-order magnetic phase transition [67,68].
The linear thermal expansion coefficient of the paramagnetic phase has been measured to be α exp L = 12 · 10 −6 K −1 at 400 K [69]. From the MDMC simulations it can be calculated by the change of the lattice constant with temperature at constant pressure divided by the lattice constant, Employing the PM data a value of α mHDNNP L = 14 · 10 −6 K −1 at 400 K is obtained in good agreement with experiment.

Defects
Real materials cannot be considered as infinitely large periodic systems, because they contain surfaces and defects breaking the ideal AFM-II order. For example, Mn vacancies lead to imperfections in the magnetic order and increase the oxidation states of other Mn ions to compensate for the excess of O atoms ensuring overall charge neutrality. The impact of these defects can be predicted by the mHDNNP. In principle, it is also possible to study the role of surfaces as well as doping, but our current parameterization is based on Mn x O bulk data only, with x = 0.969 and 1, and thus the present mHDNNP is not applicable to surface structures or doped MnO. From the heat capacities shown in the Supplementary Information [62] and the order parameters in Figure 6 the Néel temperatures can be determined to be (298±1) K for MnO, (296±1) K for Mn 0.999 O, (275±1) K for Mn 0.991 O, and (235 ± 1) K for Mn 0.969 O. Consequently, the increasing Mn vacancy concentration decreases the Néel temperature, but the change per vacancy is reduced at higher concentrations. The difference between the theoretical and experimental Néel temperature is therefore not only a consequence of the overestimation by DFT underlying the mHDNNP. A model system containing defects like vacancies, which is more comparable to real, experimental conditions, reduces the difference as well and is therefore required for accurate predictions.

CONCLUSION
In this work we have introduced spin-dependent atomcentered symmetry functions (sACSF), which enable the construction of magnetic high-dimensional neural network potentials (mHDNNP) including the full magnetic and geometric configuration space of spin-polarized, multicomponent systems. Using MnO as model systems we show that energy errors as low as 1 meV atom −1 can be reached, which is an order of magnitude lower than employing structure-dependent descriptors only. Structural changes due to the magnetic order are accurately predicted for ground and magnetically excited states with errors of only about 0.001Å and 0.1 • compared to the hybrid density functional theory (DFT) reference calculations. Furthermore, the determination of exchange coupling constants, which agree within 0.1 K with the DFT reference, demonstrates the high quality in the description of the magnetic interactions. mHDNNPs combine the accuracy and generality of first principles methods with an efficiency close to spin lattice models. The Néel temperature of MnO calculated by mean field theory differs by about 50 K from the Monte Carlo result which explicitly samples the magnetic configurations. Including structural fluctuations in the prediction of magnetic transition temperatures is essential because fixing the lattice to the low-or hightemperature configuration can lead to differences of more than 200 K. mHDNNP-driven molecular dynamics simulations including Monte Carlo spin-flips reveal a small volume increase and the disappearance of the rhombohedral distortion at the Néel temperature of MnO, as the method is able to provide mean geometric and thermodynamic data of the paramagnetic phase. We find that Mn vacancies lead to a reduced Néel temperature, which is expected to be relevant for a comparison to experimental results. For instance, changing the stoichiometry from MnO to about Mn 0.99 O the mHDNNP predicts a reduction of the Néel temperature by about 25 K.
Beyond the present work, we expect the mHDNNP method to be a powerful tool for highly accurate, largescale atomistic simulations of systems involving different spin states, like a variety of magnetic bulk materials, surfaces, and interfaces as well as molecular transition metal complexes containing spin-polarized atoms. Theoretical predictions of the magnetic, geometric, and thermodynamical implications of surfaces, interfaces, defects, and doping can provide interesting control tactics of materials properties and finally for technological applications.

High-Dimensional Neural Network Potential
The MLP used in this work is a second-generation high-dimensional neural network potential (HDNNP) [9]. In this method the potential energy E is constructed as a sum of atomic energy contributions, for a system containing N elem elements and N m atoms atoms of element m. For each element an individual feedforward neural network is trained which can provide the atomic energy contribution as a function of the local chemical environment and which is evaluated as often as atoms of the respective element are present in the system. The structural descriptors G m n are vectors of many-body atom-centered symmetry functions (ACSF) [50], which fulfill the mandatory translational, rotational, and permutational invariances of the potential energy surface and serve as input vectors of the atomic neural networks. ACSFs describe the local chemical environment of a given central atom as a function of the positions of all neighboring atoms inside a cutoff sphere of radius R c . To include all energetically relevant interactions the cutoff radius has to be sufficiently large. Besides the positions of the atoms, only the elements have to be specified leading to a reactive potential being able to describe the making and breaking of bonds. The dimensionality of the ACSF vectors can be predefined for each element individually and does not depend on the specific atomic environments. This ensures that the number of input neurons of the atomic feed-forward neural networks remains constant during molecular dynamics simulations. After optimizing the parameters of the neural networks in a training process using the potential energies and atomic force components of reference systems obtained from DFT, the HDNNP can be applied in large-scale simulations at a small fraction of the computational costs. More details about HDNNPs, their construction and validation can be found in several recent reviews [5,[10][11][12].

Atom-Centered Symmetry Functions
Two types of ACSFs are most commonly used for the construction of HDNNPs: The radial symmetry functions and the angular symmetry functions with the cutoff function R ij is the distance between central atom i and neighboring atom j, θ is the angle j −i−k involving two neighbors j and k, and η, λ, and ζ are parameters defining the spatial shapes of the ACSFs. Consequently, the ACSF values only depend on the local geometric environment of the atoms. For multicomponent systems containing several elements ACSFs for all element combinations are explicitly included. A detailed discussion of the properties of conventional ACSFs and further functional forms can be found in Reference [50].

Spin-Dependent Atom-Centered Symmetry Functions
To describe the magnetic configuration, we now introduce an atomic spin coordinate with M S is the half-difference of the number of spin-up electrons n ↑ and spin-down electrons n ↓ of an atom i in a collinear spin-polarized calculation. Consequently, the atomic spin coordinate is equal to the net direction of the atomic spin, i.e., the sign of M S , unless the absolute atomic spin value is smaller than a threshold value M thres S , which we introduce to filter out noise in the atomic spin reference data arising from the ambiguity in assigning spins in electronic structure calculations. In this work M thres S is set to 0.25. The set of atomic spin coordinates can represent all possible collinear magnetic configurations of a system enabling to identify ferro-and antiferromagnetic spin arrangements as well as non-magnetic interactions. Still, we note that different atomic oxidation and spin states with the same spin orientations cannot be distinguished by the spin coordinate alone. Such situations can often be observed, for example, for high-and low-spin states of transition metal ions. However, often the resulting different orbital occupations give rise to structural changes in the local atomic environments like changes in bond lengths or Jahn-Teller distortions and can thus be described by the conventional spatial ACSFs as shown in our previous studies [49,70] as long as there is a unique relation between the geometric structure and the electronic configuration. In principle, also the inclusion of specific spin values beyond the relative sign might be of interest, but we leave this aspect to future work here.
To integrate the spin coordinates into the radial ACSFs, the radial spin-augmentation function (SAF) M x (s i , s j ) is employed, Different radial SAFs, with x = 0, +, −, are used to describe the interactions of same (ferromagnetic interactions) and opposite spin directions (antiferromagnetic interactions) respectively, Taking the absolute values in M + and M − ensures that the descriptor is only dependent on the relative spin direction. In this way, a simultaneous sign change of all atomic spins does not change the value of the sACSFs ensuring the invariance of the potential energy with respect to the absolute spin orientation. M 0 is used to describe the non-magnetic, purely geometry-dependent interactions between a s = 0 atom and a s = 0 atom or between two s = 0 atoms, which are not included in the other terms.
When constructing a mHDNNP, it is sufficient to use sACSFs only as input for the feed-forward neural networks of elements exhibiting atoms with non-zero spins. For the feed-forward neural networks of all other elements conventional ACSFs can be used. In the same way as ACSFs, sACSFs are constructed for individual element combinations. The choice, which SAFs, i.e., only M 0 or both M + and M − , are required for a given element combination, can be made before constructing the potential because in most systems the atoms of a given element are either all characterized by s = 0 or by s = 0. For instance, in MnO the manganese atoms exhibit s = 0 and the oxygen atoms correspond to s = 0. Still, the method is also applicable to other systems including partly magnetically active elements. These systems require the usage of M 0 * (see Supplementary Information [62]) instead of M 0 to explicitly separate the non-magnetic interactions from magnetic interactions as well as a careful choice of M thres S to assign physically meaningful spin coordinate values. For element combinations of partly magnetically active elements with (partly) magnetically active elements all radial SAFs M 0 * , M + , and M − are then required.
In a similar way angular sACSFs can be defined as containing the angular SAFs M xx (s i , s j , s k ). They allow to distinguish three different interactions of a central s = 0 atom i with two neighboring s = 0 atoms j and k: (1) s i = s j = s k , (2) s i = s j = s k , and (3) s i = s j = s k . The fourth possibility s i = s k = s j is equivalent to (3) since the sums over j and k in the angular sACSFs include the interactions j − i − k and k − i − j, as both j and k sum over all contributing neighbor atoms to exclude any dependence on the order of the atoms. An efficient separation of these interactions is given by the functions,  [62]) has to be used instead of M 00 to explicitly distinguish non-magnetic and magnetic interactions and the contributions of M ++ and M −− have to be further split as described in the Supplementary Information [62].
Employing sACSFs to calculate the atomic energy contributions of atoms of elements for which at least some atoms in the system are spin-polarized, the mHDNNP is able to distinguish different magnetic configurations and can predict the corresponding potential energies. For elements including exclusively atoms with M S = 0 the usual geometry-dependent ACSFs can be applied, which is equivalent to using only M 0 and M 00 in the sACSFs of these elements. Both, the training of only several selected magnetic configurations but also the training of the full magnetic configuration space, are possible in this way.

Computational Details
The collinear spin-polarized DFT reference calculations were performed employing the Fritz-Haber-Institute ab initio molecular simulations (FHI-aims) code (version 200112.2) [71,72]. The screened hybrid exchange-correlation functional HSE06 (ω = 0.11 a 0 ) [73][74][75] and the "intermediate" FHI-aims basis set of numeric atom-centered functions excluding the auxiliary 5g hydrogenic functions were employed. A Γ-centered kpoint grid of 2×2×2 was applied to calculate the 2×2×2 supercells of MnO (64 atoms without vacancies). The convergence criterion for the self-consistency cycle was set to 10 −6 eV for the total energies and 10 −4 eVÅ −1 for the forces. Hirshfeld spin moments were used to determine the atomic spins coordinates [76]. Further details are given in the Supplementary Information [62]. An extensive benchmark for manganese oxides employing hybrid DFT functionals can be found in our previous work [77]. The sACSFs were implemented in a modified version of the RuNNer code version 1.00 [11,12,78] to construct the mHDNNP. A cutoff radius of R c = 10.5 a 0 was used. A list of the employed parameters of the n m G sACSFs for each element m is given in Tables II and III. The feed-forward neural networks consist of n m G input neurons, three hidden layers with 20, 15, and 10 neurons, respectively, and one output neuron. The mHDNNP was trained using the cohesive energies, i.e., the total energy minus the sum of the free atom energies, and using atomic force components obtained from DFT calculations of reference structures in different magnetic states. 90% of these data were used for training, i.e., optimization of the mHDNNP's weight parameters. The remaining data were used as test set. Further details about the training can be found in the Supplementary Information [62].   [79,80] and the neural network potential package (n2p2) [81] as MD library for potentials generated with RuNNer. The n2p2 code was modified to enable the usage of sACSFs. The MD simulations with MC spin-flips (MDMC) employed 6 × 6 × 6 Mn x O supercells referring to the geometric unit cell, with x = 0.969, 0.991, 0.999, 1, i.e., 1701, 1720, 1727, and 1728 atoms. They were run in the isothermalisobaric (N pT ) ensemble at a pressure of p = 1 bar with a timestep of 1 fs applying the Nosé-Hoover thermostat and barostat with coupling constants of 0.1 ps and 1 ps, respectively [82,83]. MC spin-flips were performed after each time step. The spin-flip rate is not set to measured or calculated rates to study dynamic properties but to sample the thermodynamic equilibrium efficiently. In all MDMC simulations the system was equilibrated for 1 ns before the acquisition period of 10 ns. In all conventional MC spin-flip simulations, i.e., no MD steps in between the MC spin-flips, the equilibration was performed for 10 6 steps and the acquisition consisted of 10 7 steps.

DATA AVAILABILITY
The datasets generated and analyzed during the current study are available from the corresponding author on reasonable request.