# Machine Learning Exchange Fields for Ab-initio Spin Dynamics

Yuqiang Gao  School of Physics and Electronic Information, Anhui Province Key
Laboratory for Control and Applications of Optoelectronic Information
Materials, Anhui Normal University, Wuhu 241000, PR China  University of
Twente, Faculty of Science and Technology and MESA+ Institute for
Nanotechnology, P.O. Box 217, 7500 AE Enschede, The Netherlands  Menno Bokdam
University of Twente, Faculty of Science and Technology and MESA+ Institute
for Nanotechnology, P.O. Box 217, 7500 AE Enschede, The Netherlands  Paul J.
Kelly  University of Twente, Faculty of Science and Technology and MESA+
Institute for Nanotechnology, P.O. Box 217, 7500 AE Enschede, The Netherlands

###### Abstract

We add the magnetic degrees of freedom to the widely used Gaussian
Approximation Potential of machine learning (ML) and present a model that
describes the potential energy surface of a crystal based on the atomic
coordinates as well as their noncollinear magnetic moments. Assuming an
adiabatic approximation for the spin directions and magnitudes, the ML model
depends solely on spin coordinates and orientation, resulting in computational
efficiency and enabling ab initio spin dynamics. Leveraging rotational
symmetries of magnetic interactions, the ML model can incorporate various
magnetic interactions, expanding into two-body, three-body terms, etc.,
following the spirit of cluster expansion. For simplicity, we implement the ML
model with a two-body form for the exchange interaction. Comparing total
energies and local fields predicted by the model for noncollinear spin
arrangements with explicit results of constrained noncollinear density
functional calculations for bcc Fe yields excellent results, within 1 meV/spin
for the total energy. Further optimization, including three-body and other
terms, is expected to encompass diverse magnetic interactions and enhance the
model’s accuracy. This will extend the model’s applicability to a wide range
of materials and facilitate the machine learning ab initio spin dynamics.

###### pacs:

75.70.Ak, 73.22.-f, 75.30.Hx, 75.50.Pp

##  I  Introduction

At finite temperatures, the magnetic moments of magnetic materials vary with
time in magnitude and direction influencing the potential energy surface (PES)
that describes the motion of the ions. The purpose of this publication is to
extend to itinerant magnetic materials recently developed methods for
performing first-principles (FP) molecular dynamics (MD) simulations. This
requires taking the spin degrees of freedom into consideration by extending
the concept of “machine-learning force fields” (ML-FF) to “machine-learning
exchange fields” (ML-EF) and will ultimately entail performing coupled ionic
and spin dynamics. We begin with a brief summary of the recent relevant
progress in first-principles (or ab-initio) molecular dynamics, FP-MD (or
AIMD).

In the Car-Parrinello ab-initio molecular dynamics (AIMD) scheme [1], the
interactions between the atoms comprising molecules and solids are described
in terms of full quantum mechanical (QM) solutions for the constituent
electrons. AIMD founders ultimately on the need to describe the electron and
ion dynamics on the same short time scale necessary to keep the electrons in
their ground state [2]. For metals with no gap between occupied and unoccupied
states, it turns out to be more efficient to combine long molecular-dynamics
(MD) time steps for the ions with full self-consistent field (SCF) solutions
of the Kohn-Sham equations of Density Functional Theory (DFT) [3, 4] for the
electrons. Every such step yields the total electronic energy that plays the
role of a potential energy for the ionic motion, the forces acting on all the
atoms, and the stress tensor [5]. The resulting Born-Oppenheimer MD (BO-MD)
procedure is limited by the DFT calculation to time scales of tens of
picoseconds (ps) for thousands of atoms [2] (hundreds of ps for hundreds of
atoms [6]). For comparison, typical optical vibration frequencies are in the
range $10^{12}-10^{13}$ Hertz. In these simulations, many of the
computationally expensive DFT calculations are unnecessary because at
temperatures at which bonds are only seldom broken, much of the simulation
time is spent repeatedly exploring a tiny portion of the 3$N$ dimensional
coordinate space. The effectiveness of modelling a huge diversity of systems
using periodic boundary conditions is such that $N$ is determined by the size
of the unit cell required to model the system of interest and far smaller than
the $N\approx 10^{23}$ atoms composing a typical solid. Typically $N$ should
be of order $10^{3}$ to avoid artifacts of the artificial periodicity but such
large systems are usually prohibitively expensive. For metallic or nonpolar
materials, the range of the interatomic force constants is quite limited, to
of order $\sim 30$ atoms (the long range electrostatic interactions in polar
solids can be taken care of classically without significantly increasing the
computational cost). The challenge is to find an efficient parameterization of
the DFT BO potential energy surface (PES) in this enormously reduced but still
very large coordinate space.

Behler and Parrinello decomposed the total potential energy $U_{\rm tot}$ of a
system of interacting atoms into a sum of local atomic contributions $U_{i}$

$U_{\rm tot}=\sum_{i=1}^{N}U_{i}(\\{{\bf r}_{j}\\}-{\bf r}_{i}),$ (1)

and replaced the Cartesian coordinates $\\{{\bf r}_{j}\\}$ on which $U_{i}$
depends with coordinates (“descriptors”) that reflect the translational,
rotational, reflectional and permutational symmetry of atom $i$, so-called
“atom-centered symmetry functions” (ACSF) [2]. This innovation received a
great deal of attention and has stimulated much subsequent activity. Observing
that the PES is a relatively smooth function of the coordinates on which it
depends and using the same atomic decomposition of the potential energy as
[2], Bartok et al. introduced the density $\rho_{i}({\bf r})$ at position
${\bf r}$ of the atomic neighbours of each atom $i$

$\rho_{i}({\bf r})=\sum_{j}^{r_{ij}<r_{\rm cut}}f_{\rm
cut}(r_{ij})\,\,\delta({\bf r}-{\bf r}_{ij}),$ (2)

in terms of $\delta$ functions at each atomic position. Here the index $j$
runs over the neighbors of atom $i$ at ${\bf r}_{i}$ within some radius
$r_{\rm cut}$, ${\bf r}_{ij}={\bf r}_{j}-{\bf r}_{i}$ and $r=|{\bf r}|$. The
smooth cutoff function $f_{\rm cut}(r)$ removes information about the
structure beyond $r_{\rm cut}$ and a commonly made choice is $f_{\rm
cut}(r)=[1+\cos(\pi r/r_{\rm cut})]/2$ [2, 7, 8]. $\rho({\bf r})$ should
characterize an atomic environment uniquely. Bartok et al. described the
radial dependence in terms of an angle, expanded $\rho$ in “hyperspherical”
harmonics and reduced the determination of the BO-PES to interpolating atomic
energies $U_{i}$ in this space. To do this, they introduced a nonparametric
method called “Gaussian process regression” and called the resulting PES the
“Gaussian Approximation Potential” (GAP) [7]. Because of the delta functions
used to describe the atomic positions, the expression (2) for $\rho_{i}({\bf
r})$ is not smooth and this leads to difficulties characterizing the
similarity of different atomic environments. Replacing the $\delta$ functions
with Gaussians led to a smooth measure of similarity that was termed “Smooth
Overlap of Atomic Positions” (SOAP) [9].

In spite of the great reduction in the size of the coordinate space made
possible by the approximations outlined above, a huge number of calculations
is still needed to determine the BO-PES to perform AIMD, even for simple
materials. A considerable improvement was described by Jinnouchi et al.
whereby new DFT calculations are only performed for an uncharted volume of
coordinate space when the MD simulation explores that volume of coordinate
space, so-called “on-the-fly machine learning (ML)”. The criterion for
deciding to perform a new DFT calculation is in essence based upon the
distance of the new region of coordinate space from regions already explored
and the estimate of the error in the potential energy and forces based upon a
Bayesian regression analysis that lends itself to full automation [10, 8].

In the “machine-learning force-field (ML-FF)” approach just sketched, the
electronic degrees of freedom have been effectively integrated out and
Newton’s equations of motion are solved for the ionic degrees of freedom using
an effective force field that is determined quantum mechanically and
interpolated in the coordinate space. As a result a much longer time step
appropriate to the ionic dynamics can be used. For magnetic materials we will
develop an analogous procedure whereby the spatial distribution of the spin
density will be integrated out to yield atomic moments and effective
(exchange) field which enter Landau-Lifshitz-like (LL) equations that will be
solved to describe the dynamics of atomic magnetic moments.

The LL equations describe the time variation of a magnetization ${\bf M}({\bf
r},t)$ in an effective magnetic field ${\bf H}_{\rm eff}(t)$ as

$\frac{\partial{\bf M}({\bf r},t)}{\partial t}=-\gamma{\bf M}({\bf
r},t)\times{\bf H}_{\rm eff}({\bf r},t)$ (3)

where $\gamma$ is the gyromagnetic ratio $g\mu_{B}/\hbar$,
$\mu_{B}=e\hbar/2m_{e}$ is the Bohr magneton and $g\sim 2$ and

${\bf H}_{\rm eff}({\bf r},t)=\frac{\partial F[M]}{\partial{\bf M}}$ (4)

[11, 12]. The effective field can usually be decomposed into contributions
from applied, dipolar demagnetization, crystal-anisotropy and exchange fields.
In this manuscript we will be focussing on bulk itinerant ferromagnets like Fe
so that the applied, demagnetization and anisotropy fields can be neglected by
comparison with the exchange field. In (3) and (4), ${\bf M}({\bf r},t)$ and
${\bf H}_{\rm eff}({\bf r},t)$ are vector fields and solution of the LL
equation forms the subject of “micromagnetism” [11]. Since we are interested
in a first-principles description of magnetization dynamics, we will replace
the continuum vector fields with discrete atomic quantities as in “atomistic
spin dynamics (ASD)” [13, 14]. We will present a novel descriptor for
noncollinear spins which takes into account both the spin orientation and
atomic coordinates within the GAP framework [7, 15]. Our ML model effectively
captures the DFT total energies and effective magnetic field of non-collinear
spin structures at different temperatures. This achievement demonstrates the
broad applicability and accuracy of our approach with the simplified ML model.

The paper is organized as follows. Section II introduces the methodology used
to build the descriptor for the non-collinear spin structures. Section II.1
presents the power spectrum of the SOAP descriptor. In Section II.2, we expand
the spin orientationvector in spherical harmonics and introduce the SOAP
analogue for spins, the “smooth overlap of spin orientation” (SOSO). In
Sections II.3 and II.4, we derive the descriptor and kernels for magnetic
exchange and dipole-dipole interactions, respectively. In Section III, the
performance of the ML model is tested. Because the energy of magnetic dipole-
dipole, magnetic anisotropy (and other) interactions are so small compared to
the exchange interaction in cubic Fe, we only consider the exchange
interaction here. In Section III.1, we conduct constrained self-consistent
first principle calculations for non-collinear spin systems to train and test
the ML model. In Section III.2, the ML model is tested on a Heisenberg
Hamiltonian with varying magnitude of the magnetic moment and atomic
positions. We compare our model to similar recent work in Section IV. The
adiabatic approximation for the spin simplifies the ML model, thereby
manifesting a notable enhancement in computational efficiency and enabling the
machine learning ab initio spin dynamics, as summarized in Section V.

##  II METHODOLOGY

The main challenge in developing an effective ML model for magnetic materials
is to efficiently represent the local environment seen by a spin ${\bf s}$
characterized by its magnitude $s$, its orientation ${\bf\hat{e}}={\bf s}/s$
and its position ${\bf r}$. At finite temperatures it is subject to
longitudinal and transverse fluctuations [16] that lead to a huge phase space
for the non-collinear spin system. The typical interatomic exchange parameters
$J_{ij}$ are small (e.g. $<30$ meV for Fe, Co, Ni [17]) compared to the
characteristic electronic energies such as intraatomic exchange, interatomic
hopping, etc. Thus the spin orientation corresponds to the slow degrees of
freedom ($\sim\,$ps) while the evolution of the spin magnitudes is determined
by the change of the electronic wave functions (fast degrees of freedom
$\sim\,$fs). By analogy with the Born-Oppenheimer approximation in molecular
dynamics [18], the adiabatic approximation between the rotation of the spin
orientation and fluctuations of the magnitude of the spin can be made because
of their different time scales [19, *Antropov:prb96]. Therefore, the ML model
can be constructed solely based on the spin orientation and its location
(which is just the corresponding atomic position). This will greatly contract
the phase space for the non-collinear system and improve the efficiency of the
ML model. Variation of the spin magnitude can be included at the expense of
expanding the dimension of the descriptor space to include $s$ explicitly.

A key to the success of the ML-FFs was the incorporation of the translational
and rotational invariance in the formulation of the atomic interactions. To
design the magnetic PES, mapping the local spin configurations onto a
symmetry-adapted descriptor can reduce the phase space of noncollinear spin
systems enormously. The target property should be invariant under
translational, rotational, reflectional and permutational symmetry operations.
These symmetries can be found from a magnetic Hamiltonian [21] such as

$\displaystyle H$ $\displaystyle=-J\sum_{i\neq j}{\bf s}_{i}\cdot{\bf
s}_{j}-P\sum_{i\neq j}\frac{3({\bf s}_{i}\cdot{\bf r}_{ij})({\bf
s}_{j}\cdot{\bf r}_{ij})-{\bf s}_{i}\cdot{\bf s}_{j}}{|r_{ij}|^{3}}-$ (5)
$\displaystyle-K\sum_{i}({\bf s}_{i}\cdot{\bf e}_{K})^{2}-\sum_{i\neq
j}D_{ij}\cdot[{\bf s}_{i}\times{\bf s}_{j}]$

where on the right-hand side the four most important magnetic interactions in
magnetic materials are listed. From left to right, these are the interatomic
exchange, magnetic dipole-dipole, magnetocrystalline anisotropy and
Dzyaloshinskii–Moriya (DMI) interactions. ${\bf s}_{i}$ and ${\bf s}_{j}$
represent the local spin on sites $i$ and $j$, respectively; ${\bf e}_{K}$ is
the direction of the anisotropy axis. By making use of the rotational symmetry
of the magnetic interactions, the phase space can be contracted. In this way
the ML PES would accurately capture isotropic magnetic interactions. The
energy of magnetic dipole-dipole interaction is normally of order 0.1 meV per
atom and is comparable to the size of the DMI interaction, Table 1. For the
anisotropic part that depends on the crystal geometry, including
magnetocrystalline anisotropy (that is considerably smaller in energy scale,
of the order of $\mu$eV/atom compared to isotropic exchange interactions that
is of the order of 100 meV/atom) as indicated in Table 1, an additional
correction should be incorporated. However, in this work, we will focus on the
isotropic magnetic interaction, in particular the exchange interactions that
dominate the total energy of the spin systems.

Table 1: Calculated magnetic energies of bulk bcc Fe, hcp Co and fcc Ni (meV/atom). MAE: magnetocrystalline anisotropy ($\mu$eV/atom), Co from Ref. [22], Fe and Ni from Ref.[23]; DMI: Dzyaloshinskii-Moriya interaction (meV) in thin film or multilayer interface [24, 25, 26].  | Fe | Co | Ni
---|---|---|---
$E_{\rm FM}-E_{\rm NM}$ | 561 | 208 | 58
$E_{\rm FM}-E_{\rm AFM}$ | 459 | 203 | 56
$E_{\rm MA}$ | 1.4 | 65 | 2.7
DMI | 0.15 | 0.18 | 0.12

In the spirit of the cluster expansion [27], the total energy can be expanded
in a sum of $n$-body (cluster) interaction terms ($n=1,2,...,\infty$). For
example, the simple Coulomb interaction was expanded in two-body and three-
body terms in the ML force field in [8]. Although there have been discussions
about its completeness, it would appear that the accuracy of present ML force
fields considering up to three-body terms is already sufficiently high; it has
been argued that decomposing the total energy into a sum of atom-centered
contributions mitigates the influence of fundamental deficiencies of this
approach [28]. Therefore, in this work on magnetic interactions, we will only
consider two and three-body terms. Depending on the rotational symmetry of
magnetic interactions, the spin and coordinates vectors are coupled in
different way as shown in Table 2. The exchange interaction depends on the two
spin vectors, ${\bf s}_{i}$ and ${\bf s}_{j}$, along with their separation,
$r_{ij}$. For an isotropic exchange interaction, the Hamiltonian must remain
invariant when both ${\bf s}_{i}$ and ${\bf s}_{j}$ are rotated by the same
angle. The three-body term, involving higher-order exchange interactions,
depends on the three spin vectors, ${\bf s}_{i},{\bf s}_{j}$ and ${\bf s}_{k}$
as well as the separations $d_{ij}$ and $d_{ik}$ from ${\bf s}_{i}$, where the
Hamiltonian maintains its invariance under simultaneous rotations of ${\bf
s}_{i},{\bf s}_{j}$ and ${\bf s}_{k}$. In contrast, the classical magnetic
dipole-dipole interaction in (5) depends on ${\bf s}_{i},{\bf s}_{j}$, and
their relative position vector, ${\bf r}_{ij}$. The corresponding two-body
term retains rotational invariance when ${\bf s}_{i},{\bf s}_{j},$ and ${\bf
r}_{ij}$ are rotated by the same angle. In this work, we will develop
descriptors for exchange interactions and dipole-dipole interactions with
rotational symmetries extending up to three-body terms.

Table 2: The coupled vectors that are rotationally invariant for different descriptors. “2b” and “3b” denote two-body and three-body, respectively.  | ${\bf r}_{ij}$ | ${\bf r}_{ik}$ | ${\bf\hat{s}}_{i}$ | ${\bf\hat{s}}_{j}$ | ${\bf\hat{s}}_{k}$
---|---|---|---|---|---
Exchange-2b | - | - | coupled | coupled | -
Exchange-3b | - | - | coupled | coupled | coupled
Dipole-2b | coupled | - | coupled | coupled | -
Dipole-3b | coupled | coupled | coupled | coupled | coupled

In contrast to conventional machine learning force fields, which primarily
address the Coulomb interactions based on the local atomic environments,
machine learning force fields for magnetic systems require an expansion of the
phase space to include not only atomic positions but also spin orientations.
To compare the similarity of two spin systems, we will extend the SOAP method,
originally designed to distinguish two local structures with the smooth
overlap of atomic positions, to calculate the similarity of two local spin
structures with what we will call “Smooth Overlap of Spin Orientations
(SOSO)”. Considering that magnetic interactions are influenced by atomic
separations, the SOSO framework will include the smooth overlap of atomic
positions as well.

###  II.1 Smooth overlap of atomic positions (SOAP)

In the GAP framework, the potential energy $U$ of a structure with $N_{a}$
atoms is approximated as a sum of local energies $U_{i}$ as

$U=\sum_{i=1}^{N_{a}}U_{i}.$ (6)

Each local energy $U_{i}$ is assumed to be fully determined by the local
environment of atom $i$ described by the smooth atomic density distribution

$\rho_{i}(r)=\sum_{j=1}^{N_{a}}f_{\rm cut}({\bf r}_{ij})g({\bf r}-{\bf
r}_{ij}),$ (7)

in which $g(r)$ is the normalized Gaussian function

$g({\bf r}-{\bf r}_{ij})=\frac{1}{(\sqrt{2\sigma_{\rm
atom}}\pi)^{3}}\exp\Big{(}-\frac{({\bf r}-{\bf r}_{ij})^{2}}{2\sigma_{\rm
atom}^{2}}\Big{)}.$ (8)

By expanding it in plane waves [29], the Gaussian function can be expressed in
spherical harmonics

$\displaystyle\rho_{i}({\bf r})=\sum_{j=1}^{N_{a}}f_{\rm cut}({\bf
r}_{ij})\frac{1}{(\sqrt{2\sigma_{\rm
atom}}\pi)^{3}}\exp\Big{(}-\frac{r^{2}+r_{ij}^{2}}{2\sigma_{\rm
atom}^{2}}\Big{)}$ (9) $\displaystyle 4\pi\sum_{l=0}^{L_{\rm
max}}\sum_{m=-l}^{m=l}j_{l}\Big{(}\frac{rr_{ij}}{\sigma_{\rm
atom}^{2}}\Big{)}Y_{lm}^{*}(\hat{\bf r}_{ij})Y_{lm}(\hat{\bf r}).$

where we use spherical coordinates ${\bf r}=(r,\theta,\phi)$ and the notation
${\bf\hat{r}}=(\theta,\phi)$ for the unit vector. $j_{l}$ is the spherical
Bessel function. On expanding the radial part in a set of radial basis
functions $\chi(r)$ [30], the atomic density can be rewritten as

$\rho_{i}({\bf
r})=\sum_{n=1}^{N_{R}^{l}}\sum_{l=1}^{L_{max}}\sum_{m=-l}^{l}c_{nlm}^{i}\chi_{nl}(r)Y_{lm}(\hat{\bf
r}),$ (10)

where the coefficient $c_{nlm}^{i}$ has the form

$c_{nlm}^{i}=\sum_{j=1}^{N_{a}}h_{nl}(r_{ij})Y_{lm}^{*}(\hat{\bf r}_{ij}),$
(11)

with which the rotationally invariant three-body term [8] that is equivalent
to the power spectrum [9]

$p^{i}_{nn^{\prime}l}=\sum_{m}c^{i}_{nlm}c^{i*}_{n^{\prime}lm}.$ (12)

can be constructed. The three-body term decomposes the local atomic
environment of atom $i$ into triplets involving atoms $i$, $j$ and $k$
represented by the power spectrum. Then the similarity kernel of two different
atomic configurations can be calculated by the overlap of these sets in two
different configurations as in (13).

In (11), the radial coefficient $c^{i}_{nlm}$ is expressed as a sum of
spherical harmonics describing the positions of atoms $j$ with respect to the
central atom $i$. Multiplication of two $c^{i}_{nlm}$ generates all possible
combinations of $r_{ij}$ and $r_{ik}$, which corresponds to the three-body
term. The sum over $m$ for the coupling of two spherical harmonics with the
same angular momentum number $l$ in (12) indicates the rotational invariance
according to the addition theorem of spherical harmonics whereby two angular
momenta are coupled to zero. The kernel that measures the similarity of two
local structures ${\bf X}$ and ${\bf X}^{\prime}$ can be indicated from the
smooth overlap of the atomic position (SOAP) of two structures and written as

$K({\bf X},{\bf
X}^{\prime})=\sum_{nn^{\prime}l}p_{nn^{\prime}l}p^{\prime}_{nn^{\prime}l}.$
(13)

where each vector ${\bf X}_{i}$ collects all coefficients $c^{i}_{n}$ and
$p^{i}_{nn^{\prime}l}$ for a specific local configuration $\rho_{i}({\bf r})$
.

###  II.2 Smooth overlap of spin orientation (SOSO)

By analogy with SOAP, we use a Gaussian distribution function to describe the
spin orientation ${\bf\hat{e}}$ and the overlap between two different spin
structures. Because of the “adiabatic approximation” for spin according to
which the spin orientation varies slowly in time whereas its magnitude varies
more rapidly [31], we only consider the direction of the spin on atom $i$
expressed in terms of the unit vector ${\bf\hat{e}}_{i}$

$g({\bf\hat{e}}-{\bf\hat{e}}_{i})\sim\exp\Big{(}-\frac{({\bf\hat{e}}-{\bf\hat{e}}_{i})^{2}}{2\sigma^{2}}\Big{)}$
(14)

where $\sigma$ determines the width of the Gaussian. To obtain the
normalization factor, we integrate the Gaussian function over the surface of
the unit sphere with ${\bf\hat{e}}_{i}=(0,0,1)$

$\displaystyle\int
g({\bf\hat{e}}-{\bf\hat{e}}_{i})d{\bf\hat{e}}\sim\int\limits_{0}^{\pi}\exp\Big{(}-\frac{1}{\sigma^{2}}\Big{)}\exp\Big{(}\frac{\cos\theta}{\sigma^{2}}\Big{)}d\theta$
(15) $\displaystyle=\exp\Big{(}-\frac{1}{\sigma^{2}}\Big{)}\pi
I_{0}\Big{(}\frac{1}{\sigma^{2}}\Big{)},$

where $I_{0}(x)$ is the modified Bessel function of the first kind. The
normalised Gaussian function can be rewritten as

$\displaystyle
g({\bf\hat{e}}-{\bf\hat{e}}_{i})=\frac{\exp(\frac{1}{\sigma^{2}})}{\pi
I_{0}(\frac{1}{\sigma^{2}})}\exp\Big{(}-\frac{({\bf\hat{e}}-{\bf\hat{e}}_{i})^{2}}{2\sigma^{2}}\Big{)},$
(16)

which represents the spin orientation distribution (SOD) on the unit sphere at
site $i$. Expanding in spherical harmonics, it can be rewritten as

$\displaystyle g({\bf\hat{e}}-{\bf\hat{e}}_{i})$
$\displaystyle=\frac{\exp(\frac{1}{\sigma^{2}})}{\pi
I_{0}(\frac{1}{\sigma^{2}})}4\pi\exp\Big{(}-\frac{1}{\sigma^{2}}\Big{)}$ (17)
$\displaystyle\sum_{l=0}^{L_{\rm
max}}\sum_{m=-l}^{l}j_{l}\Big{(}\frac{1}{\sigma^{2}}\Big{)}Y_{lm}^{*}({\bf\hat{e}}_{i})Y_{lm}({\bf\hat{e}})$
$\displaystyle=\frac{4}{I_{0}\Big{(}\cfrac{1}{\sigma^{2}}\Big{)}}$
$\displaystyle\sum_{l=0}^{L_{\rm
max}}\sum_{m=-l}^{l}j_{l}\Big{(}\frac{1}{\sigma^{2}}\Big{)}Y_{lm}^{*}({\bf\hat{e}}_{i})Y_{lm}({\bf\hat{e}}).$

Figure 1:  Schematic of exchange interaction between two neighbouring spins.
${\bf e}_{i}$, ${\bf e}_{j}$ represent the unit direction vectors of spins $i$
(${\bf s}_{i}$) and $j$ (${\bf s}_{j}$), respectively. ${\bf r}_{ij}$ denotes
the relative position vector between two spins.

###  II.3 Descriptor for interatomic exchange interaction

First we consider the descriptor for the dominant interatomic exchange
interaction in magnetic materials. It can be constructed by substituting the
SODs for sites $i$ and $j$ in (7). Introducing the spin degree of freedom into
(7), which adds the SODs defined in Eq.(17) for spins at sites $i$ and $j$.

$\rho_{i}({\bf r},{\bf\hat{e}}_{1},{\bf\hat{e}}_{2})=\sum_{j=1}^{N_{a}}f_{\rm
cut}(r_{ij})g({\bf r}-{\bf
r}_{ij})g({\bf\hat{e}}_{1}-{\bf\hat{e}}_{i})g({\bf\hat{e}}_{2}-{\bf\hat{e}}_{j}),$
(18)

in which ${\bf\hat{e}}_{1}$ and ${\bf\hat{e}}_{2}$ are independent unit spin
variable. Because of the short range of the exchange interaction, a cutoff
function $f_{\rm cut}(r_{ij})$ will be used to only consider the spins within
the cutoff distance.

In the classical Heisenberg Hamiltonian (5), the interatomic exchange energy
$E_{\rm ex}=\sum_{i\neq j}J_{ij}(R){\bf s}_{i}\cdot{\bf s}_{j}$. Because
$J_{ij}(R)$ is distance dependent and $E_{\rm ex}$ is invariant under
simultaneous rotations of ${\bf s}_{i}$ and ${\bf s}_{j}$, the spin density
distribution descriptor can be obtained by integrating the descriptor in (18)
over the unit sphere ${\bf\hat{r}}$ because the exchange interaction between
spins depends on their separation $r$, not their orientation ${\bf\hat{r}}$.

$\rho_{i}(r,{\bf e}_{1},{\bf e}_{2})=\int
d{\bm{\hat{r}}}\sum_{j=1}^{N_{a}}f_{\rm cut}(r_{ij})g({\bf r}-{\bf
r}_{ij})g({\bf e}_{1}-{\bf e}_{i})g({\bf e}_{2}-{\bf e}_{j}).$ (19)

Consider the rotational invariance of ${\bf s}_{i}.{\bf s}_{j}$. Under an
arbitrary rotation $\widehat{R}$ the overlap $S$ of two local spin structures
becomes

$S(\widehat{R})=\iiint
r^{2}drd{\bf\hat{e}}_{1}d{\bf\hat{e}}_{2}\rho(r,{\bf\hat{e}}_{1},{\bf\hat{e}}_{2})\rho^{\prime}(r,\widehat{R}{\bf\hat{e}}_{1},\widehat{R}{\bf\hat{e}}_{2})=\sum_{\begin{subarray}{c}m_{1},m_{2}\\\
m^{\prime}_{1},m^{\prime}_{2}\\\
l_{1},l_{2}\end{subarray}}C_{l_{1}l_{2}m_{1}m_{2}}^{*}C^{\prime}_{l_{1}l_{2}m^{\prime}_{1}m^{\prime}_{2}}D_{m_{1}m^{\prime}_{1}}^{(l_{1})}(R)D_{m_{2}m^{\prime}_{2}}^{(l_{2})}(R)$
(20)

where we use the relation

$Y_{lm}^{*}(\widehat{R}{\bf\hat{e}})=\sum_{m^{\prime}=-l}^{l}Y_{lm^{\prime}}^{*}({\bf\hat{e}})D_{m^{\prime}m}^{(l)}(\widehat{R})$
(21)

and $D_{m^{\prime}m}^{(l)}(\widehat{R})$ are the Wigner rotation matrices for
the rotation $\widehat{R}$ [9, 32].

The similarity kernel can be derived as

$\displaystyle K(\rho,\rho^{\prime})=\int
S(\widehat{R})dR=\sum_{\begin{subarray}{c}m_{1},m_{2}\\\
m^{\prime}_{1},m^{\prime}_{2}\\\
l_{1},l_{2}\end{subarray}}C_{l_{1}l_{2}m_{1}m_{2}}^{*}C^{\prime}_{l_{1}l_{2}m^{\prime}_{1}m^{\prime}_{2}}\int
D_{m_{1}m^{\prime}_{1}}^{(l_{1})}(\widehat{R})D_{m_{2}m^{\prime}_{2}}^{(l_{2})}(\widehat{R})d\widehat{R}$
(22) $\displaystyle=\sum_{\begin{subarray}{c}m_{1},m_{2}\\\
m^{\prime}_{1},m^{\prime}_{2}\\\
l_{1},l_{2}\end{subarray}}C_{l_{1}l_{2}m_{1}m_{2}}^{*}C^{\prime}_{l_{1}l_{2}m^{\prime}_{1}m^{\prime}_{2}}(-1)^{m_{1}-m^{\prime}_{1}}\int
D_{-m_{1},-m^{\prime}_{1}}^{(l_{1})}(\widehat{R})^{*}D_{m_{2}m^{\prime}_{2}}^{(l_{2})}(\widehat{R})d\widehat{R}$
$\displaystyle=\sum_{\begin{subarray}{c}m_{1},m^{\prime}_{1},l\end{subarray}}\frac{8\pi^{2}(-1)^{m_{1}+m^{\prime}_{1}}}{2l+1}C_{lm_{1},-m_{1}}^{*}C^{\prime}_{lm^{\prime}_{1},-m^{\prime}_{1}}$

where $d\widehat{R}$ indicates integration over all possible rotations [9] and

$Y^{*}_{lm}(r)=(-1)^{m}Y_{l,-m}(r).$ (23)

Here the kernel does not integrate the square of $S(R)$ as in the case of
SOAP. Because exchange interaction is in principle a two-body interaction, the
kernel defined here only includes the two-body term of the exchange
interaction. For three-body or higher order terms, it can be constructed from
the descriptor defined in (44) (see Appendix A).

Now we have given the full expression for the kernel of two-body term of
exchange interaction. The exchange interaction does not depend on the
direction of ${\bf r}_{ij}$, so only the scalar radial distribution part
$r_{ij}$ left. For computational efficiency, we expend the radial part of spin
density into orthogonal radial basis functions $\chi_{nl}(r)$. Then the radial
distribution density of spin $\rho_{i}(r)$ can be rewritten as

$\displaystyle\rho_{i}(r)=\frac{1}{\sqrt{4\pi}}\sum_{n=1}^{N_{R}}C_{n00}^{i}\chi_{nl}(r).$
(24)

Therefore, by summing the coefficients in (22) over $m$, we can get

$\displaystyle
C_{nl}^{i}=\sum_{j}^{N_{a}}C_{n00}^{i}\Big{(}\frac{4}{I_{0}(\sigma^{-2})}\Big{)}^{2}\sum_{m=-l}^{l}j_{l}^{2}(\sigma^{-2})Y_{lm}^{*}({\bf\hat{e}}_{i})Y_{lm}({\bf\hat{e}}_{j})$
(25)
$\displaystyle=\sum_{j}^{N_{a}}C_{n00}^{i}\Big{(}\frac{4}{I_{0}(\sigma^{-2})}\Big{)}^{2}j_{l}^{2}(\sigma^{-2})\frac{2l+1}{4\pi}P_{l}(\cos\theta_{ij})$

Adding the prefactor $\sqrt{\frac{8\pi^{2}}{2l+1}}$,

$C_{nl}^{i}=\sum_{j}^{N_{a}}C_{n00}^{i}\left(\frac{4}{I_{0}(\cfrac{1}{\sigma^{2}})}\right)^{2}j_{l}^{2}(\frac{1}{\sigma^{2}})\sqrt{\frac{2l+1}{8\pi}}P_{l}(\cos\theta_{ij}).$
(26)

The similarity kernel can be rewritten as

$K(\rho_{i},\rho^{\prime}_{i})=\sum_{n,l}C_{nl}^{i*}C_{nl}^{{}^{\prime}i}.$
(27)

###  II.4 Descriptor for magnetic dipole-dipole like interaction

Here, we also give the descriptor for magnetic dipole-dipole like interactions
that maintain the rotational symmetry of ${\bf s}_{i},{\bf s}_{j},$ and ${\bf
r}_{ij}$.

$\rho_{i}({\bf
r},{\bf\hat{e}}_{1},{\bf\hat{e}}_{2})=\sum_{j=1}^{N_{a}}f_{cut}(r_{ij})g({\bf
r}-{\bf
r}_{ij})g({\bf\hat{e}}_{1}-{\bf\hat{e}}_{i})g({\bf\hat{e}}_{2}-{\bf\hat{e}}_{j}).$
(28)

The overlap of two local spin structures

$\displaystyle S(\widehat{R})$ $\displaystyle=\int\\!\\!\\!\int\\!\\!\\!\int
r^{2}d{\bf\hat{r}}d{\bf\hat{e}}_{1}d{\bf\hat{e}}_{2}\rho({\bf
r},{\bf\hat{e}}_{1},{\bf\hat{e}}_{2})\rho^{\prime}(\widehat{R}{\bf
r},\widehat{R}{\bf\hat{e}}_{1},\widehat{R}{\bf\hat{e}}_{2})$ (29)
$\displaystyle=\sum_{\begin{subarray}{c}m,m_{1},m_{2}\\\
m^{\prime},m^{\prime}_{1},m^{\prime}_{2}\\\
l,l_{1},l_{2}\end{subarray}}C_{ll_{1}l_{2};mm_{1}m_{2}}^{*}C^{\prime}_{ll_{1}l_{2};m^{\prime}m^{\prime}_{1}m^{\prime}_{2}}$
$\displaystyle
D_{mm^{\prime}}^{l}({\bf\hat{e}})D_{m_{1}m^{\prime}_{1}}^{l_{1}}({\bf\hat{e}})D_{m_{2}m^{\prime}_{2}}^{l_{2}}({\bf\hat{e}}).$

To get the similarity kernel, we need to transform the Wigner D-matrix to be
complex conjugate of another D-matrix by

$D_{mm^{\prime}}^{l}(R)=(-1)^{m-m^{\prime}}D_{-m,-m^{\prime}}^{l}(R)^{*}.$
(30)

Then, the integrate of three D-matrix can be expressed as

$\displaystyle\int
D_{mm^{\prime}}^{l}(R)^{*}D_{m_{1}m^{\prime}_{1}}^{l_{1}}(R)D_{m_{2}m^{\prime}_{2}}^{l_{2}}(R)dR$
(31)
$\displaystyle=\frac{8\pi^{2}}{2l+1}<lm|l_{1}m_{1};l_{2}m_{2}>\times<lm^{\prime}|l_{1}m^{\prime}_{1};l_{2}m^{\prime}_{2}>,$

in which the angle brackets are Clebsch–Gordan coefficients.

Finally, the similarity kernel can be written as

$\displaystyle K(\rho,\rho^{\prime})$ $\displaystyle=\int
S(R)dR=\sum_{\begin{subarray}{c}m,m_{1},m_{2}\\\
m^{\prime},m^{\prime}_{1},m^{\prime}_{2}\\\
l,l_{1},l_{2}\end{subarray}}C_{ll_{1}l_{2}mm_{1}m_{2}}^{*}C^{\prime}_{ll_{1}l_{2}m^{\prime}m^{\prime}_{1}m^{\prime}_{2}}\int
D_{mm^{\prime}}^{l}(R)D_{m_{1}m^{\prime}_{1}}^{l_{1}}(R)D_{m_{2}m^{\prime}_{2}}^{l_{2}}(R)dR$
(32) $\displaystyle=\sum_{\begin{subarray}{c}m,m_{1},m_{2}\\\
m^{\prime},m^{\prime}_{1},m^{\prime}_{2}\\\
l,l_{1},l_{2}\end{subarray}}C_{ll_{1}l_{2}mm_{1}m_{2}}^{*}C^{\prime}_{ll_{1}l_{2}m^{\prime}m^{\prime}_{1}m^{\prime}_{2}}(-1)^{m-m^{\prime}}\int
D_{-m,-m^{\prime}}^{l}(R)^{*}D_{m_{1}m^{\prime}_{1}}^{l_{1}}(R)D_{m_{2}m^{\prime}_{2}}^{l_{2}}(R)dR$
$\displaystyle=\sum_{\begin{subarray}{c}m,m_{1},m_{2}\\\
m^{\prime},m^{\prime}_{1},m^{\prime}_{2}\\\
l,l_{1},l_{2}\end{subarray}}C_{ll_{1}l_{2}mm_{1}m_{2}}^{*}C^{\prime}_{ll_{1}l_{2}m^{\prime}m^{\prime}_{1}m^{\prime}_{2}}\frac{8\pi^{2}(-1)^{m+m^{\prime}}}{2l+1}<l,-m|l_{1}m_{1};l_{2}m_{2}>\times<l,-m^{\prime}|l_{1}m^{\prime}_{1};l_{2}m^{\prime}_{2}>.$

The Clebsch–Gordan coefficients will only be nonzero when

$\displaystyle|l_{1}-l_{2}|<l<l_{1}+l_{2}$ (33a) $\displaystyle
m_{1}+m_{2}=-m$ (33b) $\displaystyle
m^{\prime}_{1}+m^{\prime}_{2}=-m^{\prime}$ (33c)

Therefore, the coefficient in (32) can be rewritten as radial part
$C_{nlm_{1}+m_{2}}^{i}$ and spin vector part $C_{l_{1}l_{2},m_{1}m_{2}}^{i}$

$C_{nlm_{1}+m_{2}}^{i}=\sum_{j=1}^{N_{a}}h_{nl}(r_{ij})Y_{lm_{1}+m_{2}}({\bf\hat{r}}_{ij})$
(34)

$\displaystyle
C_{l_{1}l_{2},m_{1}m_{2}}^{i}=\sum_{j}^{N_{a}}\Bigg{(}\frac{4}{I_{0}(\cfrac{1}{\sigma^{2}})}\Bigg{)}^{2}j_{l_{1}}\Big{(}\frac{1}{\sigma^{2}}\Big{)}j_{l_{2}}\Big{(}\frac{1}{\sigma^{2}}\Big{)}$
(35) $\displaystyle
Y_{l_{1}m_{1}}^{*}({\bf\hat{e}}_{j})Y_{l_{2}m_{2}}^{*}({\bf\hat{e}}_{k})$
$\displaystyle\sqrt{\frac{8\pi^{2}}{2l+1}}<l,m_{1}+m_{2}|l_{1}m_{1};l_{2}m_{2}>.$

Now we have got the similarity kernel for two-body term of magnetic dipole-
dipole interaction. For the three-body term, we can not use the traditional
power spectrum, because the spin direction at spin $i$ is shared by the spin
$j$ and $k$. Here, we can define the three body spin distribution as

$\displaystyle\rho_{i}(r_{1},r_{2},e_{1},e_{2},e_{3})$
$\displaystyle=\sum_{j\neq k}^{N_{a}}f_{cut}({\bf r}_{ij})g({\bf r}_{1}-{\bf
r}_{ij})g({\bf r}_{2}-{\bf r}_{ik})$ (36)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!g({\bf e}_{1}-{\bf
e}_{i})g({\bf e}_{2}-{\bf e}_{j})g({\bf e}_{3}-{\bf e}_{k}).$

Figure 2:  Comparison of total DFT energies calculated for non-collinear
configurations of 16 atom supercells of Fe with those predicted using the ML
model at (a) 300K and (d) 1000K, respectively. All energies are given with
respect to the energy of the corresponding collinear system and RMSE is the
root mean square error. The ML exchange field model is trained with 25
different noncollinear spin configurations. (b,c) and (e,f) Comparison of the
effective magnetic fields (arrows) on the Fe atoms from explicit DFT
calculations (b) and (e) with the ML predictions (c) and (f) at 300K (b) and
(c) and at 1000K (e) and (f), respectively. The unit of the effective magnetic
field is eV/$\mu_{B}$ and is scaled by a factor of 4 for visualization
purposes.

##  III Training and testing

Because the interatomic exchange interaction constitutes the primary
contribution to the total energy in ferromagnets as indicated in Table 1, a
GAP-SOSO ML model is constructed analogous to the GAP-SOAP framework,
incorporating a two-body term of the exchange interaction, see (26). To assess
the model’s performance we carried out tests for noncollinear spin structures
of bcc Fe. We began by generating noncollinear training sets for $2\times
2\times 2$ supercells of a simple cubic description of bcc Fe containing 16
atoms with temperatures ranging from 0K to 1000K. To do this, we used the
Uppsala atomic spin dynamics (UppASD) package [13, 14]. Once equilibrium was
reached at each temperature, 50 spin structures were selected at random and
divided into sets of 25 for training with the remaining 25 being used for
testing.

###  III.1 Tests on DFT constrained calculations

To calculate the total energy of the non-collinear spin systems accurately, we
performed constrained density functional theory calculations with the Vienna
Ab-initio Simulation Package (vasp). The spin configurations were constrained
so their orientations matched those of the training set output from UppASD and
the perpendicular part of the internal effective field, about which the spins
precess, was determined as the negative of the constraining field [33]
Exchange and correlation effects were described in the generalized gradient
approximation as given by Perdew, Burke, and Ernzerhof [34]. An energy cut-off
of 500 eV with energy converging to $10^{-6}\,$eV was used. An equivalent
36$\times$36$\times$36 k-point sampling for 1$\times 1\times$ was used. The
DFT energies of the noncollinear spin systems used for training were defined
with respect to the corresponding collinear ferromagnetic system as

$E_{\rm DFT}=E_{\rm FM}^{\rm nc}-E_{\rm FM}^{\rm c},$ (37)

where $E_{\rm FM}^{\rm nc}$ and $E_{\rm FM}^{\rm c}$ are the total energy for
noncolllinear and collinear systems, respectively. Because of the small size
of the magnetocrystalline anisotropy energy (MAE) for cubic systems [35], spin
orbit coupling was not included. Its effect will be studied at some later date
for lower symmetry systems with larger MAEs [36, 37].

In the ML model, a cut-off radius of 7 Å was used to account for the exchange
interaction between Fe atoms. The maximum values of $n$ and $l$ for the
spherical harmonics were set to 12 and 6, respectively.Fig. 2 compares total
DFT energies and effective magnetic fields calculated for noncollinear
configurations with the corresponding quantities predicted by the ML model at
300 K and 1000 K. The results demonstrate that the ML magnetic PES, trained
with a very small set of 25 spin structures, accurately predicts the total
energies and effective exchange fields of noncollinear spin structures. The
root-mean-square error (RMSE) of the total energy is less than 1 meV/spin at
300 K, increasing to 1.44 meV/spin at 1000K. The larger RMSE at 1000 K
compared to that at 300K suggests the need for a more extensive training set
to represent the phase space of spin structures at higher temperatures more
faithfully. To obtain the complete magnetic PES of bcc Fe may require training
the model with sets that sample a larger part of the entire phase spaces
ranging from 0K to the phase transition temperature.

The effective magnetic field shown in Fig. 2 is the transverse component of
the total effective magnetic field for each spin which generates the torque
that describes the spin precession. This component is obtained using the
expression

$\displaystyle{\bf B}_{\perp}^{i}=-\frac{\partial{E_{i}({\bf\hat{e}}_{i},{\bf
r}_{i})}}{\partial{{\bf\hat{e}}_{i}}}.$ (38)

One notable point is that the total energy learned by our model does not
depend on the spin magnitude $|s_{i}|$, the ${\bf B}_{\perp}^{i}$ can be
obtained directly by the partial derivative of the total energy with respect
to the spin unit vector ${\bf\hat{e}}_{i}$. This is confirmed by the results
shown in the right-hand panels of Fig. 2. While ${\bf B}_{\perp}^{i}$ gives
rise to a torque on the spin, the longitudinal part of the effective field
primarily governs the magnitude of the spin, which is not included in the
present model; the noncollinear spin configurations for training and
prediction only involve spin directions.

In principle, the exchange field which constrains the spin to the preferred
direction in DFT calculations is equivalent to ${\bf B}_{\perp}^{i}$ so we can
compare the predicted ${\bf B}_{\perp}^{i}$ and the constraining magnetic
field obtained from DFT calculations. As seen in Fig. 2, the ML model
accurately predicts the perpendicular component of the effective magnetic
field which can be leveraged in spin dynamics simulations. While the spin
direction from the DFT and predictions align closely (maximum $5^{\circ}$
difference), we observe that the magnitude of the effective exchange field
from DFT calculations exceeds the prediction by a factor ranging from 1.0 to
1.2 on a few atomic sites, see Fig. 2 (e,f). This discrepancy can be
attributed to the non-zero constraining penalty term in the constrained DFT
calculations [33], where the spin direction approaches the preferred direction
as the penalty term goes to zero. The constraining magnetic field obtained
from constrained DFT calculations (as implemented in vasp) is given by

$\displaystyle{\bf B}_{\rm
Con}^{i}=-\frac{\partial{E_{P}}}{\partial{\bm{s}_{i}}}=-\lambda_{i}({\bf\hat{e}}_{i}-{\bf\hat{e}}_{i}^{F})$
(39)

in which $E_{P}$ denotes the penalty energy, $\lambda_{i}$ represents the
parameter used in the constrained DFT calculations and ${\bf\hat{e}}_{i}$ and
${\bf\hat{e}}_{i}^{F}$ represent the unit spin direction and the preferred
spin direction of atom $i$ [33]. From (39), we can see that the constraining
magnetic field obtained from DFT is either perpendicular to the spin direction
or the preferred spin direction. Estimating the perpendicular exchange field
using (38) generates an exchange field perpendicular to ${\bf\hat{e}}_{i}$. As
a consequence, the magnitude and direction of the predicted exchange field
maydeviate from he DFT value but this can be improved as the penalty term
converges to zero. However, this disparity does not compromise the accuracy of
the effective magnetic field prediction with our ML model. It is worth noting
that we train the ML model using spin configuration {${\bf\hat{e}}_{i}$}
obtained from constrained DFT calculations, rather than the preferred spin
configuration {${\bf\hat{e}}_{i}^{F}$}.

A very useful feature of the MLFF is its portability allowing the model
trained on small supercells to accurately predict the properties of
significantly larger supercells. To assess the portability of our model, we
employed the ML model trained with bcc Fe in a 16 atom supercell at 300K to
predict the total energy of bcc Fe in $3\times 3\times 3$ cubic supercell
containing 54 Fe atoms also at the same temperature. Because of the
computational expense of constrained noncollinear calculations in larger
supercells, only 12 samples were selected for testing. As demonstrated in Fig.
3, the trained model effectively predicts the total energy of the larger
supercell, with the RMSE remaining consistent compared to that in the $2\times
2\times 2$ supercell.

Figure 3:  Comparison of total DFT energies calculated for non-collinear
configurations of 54 atom supercells of Fe with those predicted using the ML
model. The ML model are trained with 2 $\times$ 2 $\times$ 2 supercell of bcc
Fe.

###  III.2 Tests on Heisenberg Hamiltonian

The simplified Hamiltonian based on Heisenberg model can be constructed as

$\displaystyle E=\sum_{i\neq j}J_{ij}(R){\bf s}_{i}\cdot{\bf s}_{j}.$ (40)

where only the exchange interaction is included and the exchange parameters
$J_{ij}(R)$ were fitted from DFT calculations. 2550 random noncollinear spin
configurations were produced in a $2\times 2\times 2$ supercell of bcc Fe,
with the average angle between spins varing from $0$ to $\pi$. Among these
configurations, 50 were allocated for training the model, while the remaining
2500 were reserved for testing. The Heisenberg Hamiltonian was utilized to
compute the total energy of each spin configuration. As depicted in Fig. 4
(a), the nearly perfect prediction achieved by the ML model, with just 50
training sets, underscores its exceptional performance. The error
distribution, as shown in the inset, predominantly centers around 0 with an
RMSE of 0.07 meV/atom. Notably, in this Heisenberg Hamiltonian, the magnetic
moment remains fixed at 2.23 $\mu_{B}$ and the exchange interaction parameter
soley reliant on the distance between two spins.

In a real system, the variations in magnetic moment and exchange interactions
with temperature, particularly in noncollinear spin configurations, render
accurate predictions challenging with a simplistic model that only considers
spin directions. As illustrated in Fig. 7, the magnetic moment of bcc Fe
shifts from 2.23 $\mu_{B}$ in FM state to 1.60 $\mu_{B}$ with antiparallel
spins at the equilibrium lattice constant of GGA. According to the ”adiabatic
approximation” of spin, the magnitude of the moment should be slaved to the
spin configurations. Our ML model that solely includes the spin direction
should be enough to predict the total energy. However, to check the impact of
the variation of the local magnetic moment on the performance of our model, we
constructed a Heisenberg Hamiltonian with the magnetic moment of spin $i$ set
by

$\displaystyle M_{i}=2.2-\sum_{i\neq j}^{j}\frac{\arccos({\bf e}_{i}\cdot{\bf
e}_{j})}{n\pi}\times 0.6,$ (41)

in which n represents the total number of neighbours of spin $i$. As indicated
in Fig. 4 (b), the spin configuration dependent magnetic moment elevates the
RMSE of the prediction to 0.84 meV/atom, still remaining small compared to the
total energy. As shown in Fig. 2, the RMSE increases at 1000 K compared to 300
K, attributed not only to the larger phase space at 1000K but also to the
variation of the magnetic moment due to increased spin disorder at higher
temperatures. The prediction accuracy can be enhanced with larger training
sets. Although the changes in magnetic moment are not incorporated in our ML
model, it proficiently predicts the total energy of the Heisenberg model. The
energy scale ($\sim$ 4-10 eV) of the Heisenberg model are much larger than
that in DFT calculations, primarily because the angels between spins in the
Heisenberg model are distributed in the range [0,$\pi$], whereas the angles in
DFT calculations obtained from Uppsala ASD simulations are much smaller ($\sim
0.4\pi$).

Figure 4:  The comparison between total energy from the prediction and
Heisenberg Hamiltonian with (a) fixed magnetic moment and spin configuration
dependent magnetic moment in 2 $\times$ 2 $\times$ 2 supercell of bcc Fe. The
inset shows the error distribution in the unit of meV/atom.

It is worth noting that the noncollinear spin systems in both Fig. 2 and Fig.
4 are situated at fixed atomic positions. In our ML model, the descriptor
relies on both the spin direction and atomic coordinates. To evaluate the
model’s performance on systems with varying spin configurations and atomic
distributions, we introduce a Gaussian distribution of lattice displacement in
the noncollinear spin systems from their equilibrium positions. The Gaussian
probability distribution is expressed as

$P(\Delta r)=\frac{1}{\sqrt{2\pi}\sigma}\exp\Big{(}-\frac{(\Delta
r)^{2}}{2\sigma^{2}}\Big{)},$ (42)

where $\sigma$ represents the broadening of the displacement and $\Delta r$
denotes the displacement length. The Curie temperature of bcc Fe is around
1043K, corresponding to a thermal energy of 82 meV. To mimic the spin
configuration and lattice displacement in actual systems, the average angle
between spin at 1000 K can be estimated to $0.4\pi$, considering the energy
needed to flip the spin shown in Table 1. As illustrated in Fig. 7, the
lattice expands slightly with temperature, e.g. 1.5% at Curie temperature.
Here we choose the $\sigma$ to be 0.07 Åand 0.25 Å, which correspond to 1.2 %
and 4.4 % of lattice expansion, respectively.

As proposed in , the distance-dependent exchange interaction $J(r)$ can be
parameterized using the double exchange interaction form

$J(r)=\left\\{\begin{array}[]{lr}\frac{c}{r^{3}}e^{-r/r_{0}},&{\rm if}\quad
r\leq r_{c}\vspace{1ex}\\\ 0,&{\rm otherwise}.\end{array}\right.$ (43)

in which $r$ denotes the distance between Fe atoms, and $r_{c}$ represents the
cutoff length of the exchange interaction. In the Heisenberg model, $c$ and
$r_{0}$ are set to be 0.94 eV/Å3 and 2.0 Å by fitting the DFT calculations,
respectively. Although it can not exactly reproduce the exchange interaction
of bcc Fe, this model serves only to test the model’s performance on the spin-
lattice coupling. Because the exchange interaction between the nearest and
second nearest neighbours contribute significantly to bcc Fe, $r_{c}$ is set
to 1.4 times the lattice constant.

Figure 5:  The comparison between total energy of the prediction and
Heisenberg Hamiltonian with atomic displacements in 2 $\times$ 2 $\times$ 2
supercell of bcc Fe. The atomic displacements are sampled in Gaussian
distribution with (a) $\sigma=0.07$ Å and (b) $0.25$ Å.

Fig. 5 demonstrates the performance of the ML model on the spin-lattice
coupling. The model can well predict the total energy for systems with diverse
spin and atomic structures. The RMSE marginally increases as the atomic
displacement broadens from 0.07 Å to 0.25 Å. Furthermore, it should be
acknowledged that the atomic displacement is expected to influence the
magnetization density and consequently the magnitude of the moment, as
indicated in Appendix B. However, this is beyond the scope of this paper.

##  IV Discussion

In this paper, we introduced a novel descriptor for noncollinear spin systems
within the GAP-SOAP framework. Under the ”spin adiabatic approximation”, the
descriptor relies solely on the spin directions and the corresponding atomic
coordinates where the spins are situated. By evaluating the smooth overlap of
the spin directions in two local spin configurations, the similarity betweeen
them can be accurately assessed. The formulation of the descriptor and
similarity kernel for both the exchange interaction and magnetic dipole-dipole
interaction has been elucidated, and a descriptor incorporating a two-body
term of the exchange interaction has been implemented and tested. The
remarkable performance in predicting the total energy and internal effective
magnetic field of noncollinear spin systems with a relative small training set
underscores the robust capability of the proposed model.

###  IV.1 ML model vs the magnetic cluster expansion

The equivalence between the GAP-SOAP power spectrum and the atomic cluster
expansion with three-body term in spherical harmonics [38] suggests that the
cluster expansion including spin and atomic coordinates will be analogous to
the descriptor for spin and coordinates within the framework of GAP-SOAP. In
the atomic cluster expansion, the angular momentum coupled to zero (L=0)
signifies the rotation invariance. Recently, Matteo Rinaldi proposed a
complete set for noncollinear spin systems based on the magnetic cluster
expansion [39]. In comparison to their model, we have introduced an ”adiabatic
approximation” for the spin magnitude and spin direction, simplifying the
model and enhancing its computational efficiency, while the training sets in
atomic cluster expansion are in the order of $10^{4}$ to achieve proper
accuracy [39]. The estimated variation in the moment’s magnitude of bcc Fe, as
the spin configurations transition from FM to AFM, is approximately 26 %, as
deduced from Appendix B. At 1000K, the average spin angles is approximately
$0.4\pi$, obtained from spin dynamics simulations, suggesting that the
variation of the magnetic moment below the Curie temperature of bcc Fe is
restricted to roughly 10 %. As denoted in LABEL:eq:perp_B, the perpendicular
component of the effective magnetic field relies solely on the spin unit
vector rather than the spin magnitude. Consequently, the internal magnetic
field can be accurately predicted without including the spin magnitude.

Our work have demonstrated that a two-body term of exchange interaction alone
can sufficiently capture the exchange field. Another significant distinction
lies in our incorporation of rotational invariance for various magnetic
interactions. Instead of encompassing the two-spin, three-spin, and many-spin
interaction term as in the atomic cluster expansion, we have constructed two-
body, three-body, and many-body terms for different magnetic interactions,
each possessing unique rotational symmetries. This approach makes our ML model
more efficient and precise to include the magnetic interactions in interest.

###  IV.2 ML model vs Spectral neighbor representation

ML models based on the spectral neighbor representation for noncollinear spin
have been proposed [40]. This model introduces the rotational invariance of
the neighbour vector field. However, the system is invariant under the
rotation of the entire system and does not account for the rotational
symmetries for different kinds of magnetic interactions. For instance, the
exchange interaction should remian invariant under the rotation of $s_{i}$ and
$s_{j}$. Moreover, the vector spin on the rotation origin will destroy this
kind of rotational invariance, although specific techniques can be applied to
solve this issue [40].

Suzuki [41] extended the spectral neighbour representation to the GAP-SOAP
framework, which and is quite similar to our model. Their model does not rely
on the rotation origin but the rotational invariance for different types of
magnetic interactions are not specified in the model. They found that a
higher-order partial spectrum (trispectrum) is necessary to accurately
distinguish magnetic structures with different magnetic anisotropy. However,
the trispectrum is equivalent to the five-body term descriptor for the
magnetic dipole-dipole interaction in our model, which demands an extensive
training set. By utlizing the spin adiabatic approximation and incorporating
rational invariance for specific magnetic interactions, such as exchange
interaction, we have demonstrated that only a two-body term of exchange
interaction in our model can accurately describe the total magnetic energy and
internal effective magnetic field through training on DFT constrained
calculations with a relative small training set.

##  V Conclusion

In this paper, we have developed a highly accurate and efficient descriptor
for noncollinear spin systems within the GAP-SOAP framework. The descriptors
are designed with specific rotational invariant symmetries, tailored to
different types of magnetic interactions. Leveraging the spin ”adiabatic
approximation” in the ML model has enhanced its efficiency, requiring only
small training sets to achieve robust performance. Through constrained DFT
calculations, we trained the ML model that only incorporates a simple two-body
term of the exchange interaction with 25 different noncollinear spin
configurations at each temperature. Despite the potential numerical error
induced by the non-zero constraining penalty term in the DFT calculations, the
ML model accurately forecasts the total energy and internal effective magnetic
field with exceptional precision (around 1 meV/spin for the total energy and a
maximum 5-degree difference for the magnetic field) relative to the DFT
outcomes. In conclusion, an efficient ML model for the noncollinear spin
systems proposed in this paper can predict the total energy and effective
magnetic field in the accuracy of DFT, at a significantly lower computational
cost. Furthermore, this advancement paves the way for machine learning spin
dynamics and a combination of ML spin dynamics and molecular dynamics can be
expected, enabling investigation into various phenomena stemming from spin-
lattice coupling effects, such as magneto-elastic effects, and magnon-phonon
interactions.

###### Acknowledgements.

This work was financially supported by the “Nederlandse Organisatie voor
Wetenschappelijk Onderzoek” (NWO) through the research programme of the former
“Stichting voor Fundamenteel Onderzoek der Materie,” (NWO-I, formerly FOM) and
through the use of supercomputer facilities of NWO “Exacte Wetenschappen”
(Physical Sciences). Y.G. thanks AHNU start-up grant for financial support.

##  Appendix A Three body term for magnetic exchange like interaction

Three body term descriptor

$\rho_{i}(r_{1},r_{2},{\bf e}_{1},{\bf e}_{2},{\bf e}_{3})=\iint
d{\bf\hat{r}}_{1}d{\bf\hat{r}}_{2}\sum_{j\neq k}^{N_{a}}f_{\rm
cut}(r_{ij})g({\bf r}_{1}-{\bf r}_{ij})g({\bf r}_{2}-{\bf r}_{ik})g({\bf
e}_{1}-{\bf e}_{i})g({\bf e}_{2}-{\bf e}_{j})g({\bf e}_{3}-{\bf e}_{k})$ (44)

The overlap $S$ of the three-body term descriptor is

$\displaystyle S(R)=\iiint
r_{1}^{2}dr_{1}r_{2}^{2}dr_{2}de_{1}de_{2}de_{3}\rho(r_{1},r_{2},e_{1},e_{2},e_{3})\rho^{\prime}(r_{1},r_{2},Re_{1},Re_{2},Re_{3})$
(45) $\displaystyle=\sum_{\begin{subarray}{c}m,m_{1},m_{2}\\\
m^{\prime},m^{\prime}_{1},m^{\prime}_{2}\\\
l,l_{1},l_{2}\end{subarray}}C_{ll_{1}l_{2}m_{1}m_{2}m_{3}}^{*}C^{\prime}_{ll_{1}l_{2}m^{\prime}m^{\prime}_{1}m^{\prime}_{2}}D_{mm^{\prime}}^{l}(R)D_{m_{1}m^{\prime}_{1}}^{l_{1}}(R)D_{m_{2}m^{\prime}_{2}}^{l_{2}}(R)$

The kernel for the three-body term descriptor

$\displaystyle K(\rho,\rho^{\prime})$ $\displaystyle=\int
S(R)dR=\sum_{\begin{subarray}{c}m,m_{1},m_{2}\\\
m^{\prime},m^{\prime}_{1},m^{\prime}_{2}\\\
l,l_{1},l_{2}\end{subarray}}C_{ll_{1}l_{2}mm_{1}m_{2}}^{*}C^{\prime}_{ll_{1}l_{2}m^{\prime}m^{\prime}_{1}m^{\prime}_{2}}\int
D_{mm^{\prime}}^{l}(R)D_{m_{1}m^{\prime}_{1}}^{l_{1}}(R)D_{m_{2}m^{\prime}_{2}}^{l_{2}}(R)dR$
(46) $\displaystyle=\sum_{\begin{subarray}{c}m,m_{1},m_{2}\\\
m^{\prime},m^{\prime}_{1},m^{\prime}_{2}\\\
l,l_{1},l_{2}\end{subarray}}C_{ll_{1}l_{2}mm_{1}m_{2}}^{*}C^{\prime}_{ll_{1}l_{2}m^{\prime}m^{\prime}_{1}m^{\prime}_{2}}(-1)^{m-m^{\prime}}\int
D_{-m,-m^{\prime}}^{l}(R)^{*}D_{m_{1}m^{\prime}_{1}}^{l_{1}}(R)D_{m_{2}m^{\prime}_{2}}^{l_{2}}(R)dR$
$\displaystyle=\sum_{\begin{subarray}{c}m,m_{1},m_{2}\\\
m^{\prime},m^{\prime}_{1},m^{\prime}_{2}\\\
l,l_{1},l_{2}\end{subarray}}C_{ll_{1}l_{2}mm_{1}m_{2}}^{*}C^{\prime}_{ll_{1}l_{2}m^{\prime}m^{\prime}_{1}m^{\prime}_{2}}\frac{8\pi^{2}(-1)^{m+m^{\prime}}}{2l+1}<l,-m|l_{1}m_{1};l_{2}m_{2}>\times<l,-m^{\prime}|l_{1}m^{\prime}_{1};l_{2}m^{\prime}_{2}>$

The atomic densities are written as

$\displaystyle\rho_{i}(r_{1})=\frac{1}{\sqrt{4\pi}}\sum_{n=1}^{N_{R}}C_{n}^{i}\chi_{nl}(r_{1})\;\;\;\;;\;\;\;\;\rho_{i}(r_{2})=\frac{1}{\sqrt{4\pi}}\sum_{n=1}^{N_{R}}C_{n}^{i}\chi_{nl}(r_{2})$
(47) $\displaystyle
C_{ll_{1}l_{2},m_{1}m_{2}}^{i}=\sum_{j}^{N_{a}}\Bigg{(}\frac{4}{I_{0}(\cfrac{1}{\sigma^{2}})}\Bigg{)}^{3}j_{l}\Big{(}\frac{1}{\sigma^{2}}\Big{)}j_{l_{1}}\Big{(}\frac{1}{\sigma^{2}}\Big{)}j_{l_{2}}\Big{(}\frac{1}{\sigma^{2}}\Big{)}Y_{lm_{1}+m_{2}}({\bf\hat{e}}_{i})Y_{l_{1}m_{1}}^{*}({\bf\hat{e}}_{j})Y_{l_{2}m_{2}}^{*}({\bf\hat{e}}_{k})$
$\displaystyle\times\sqrt{\frac{8\pi^{2}}{2l+1}}<l,m_{1}+m_{2}|l_{1}m_{1};l_{2}m_{2}>$

Finally, the coefficients for the kernel become

$\displaystyle
C_{nn^{\prime}ll_{1}l_{2},m_{1}m_{2}}^{i}=\frac{1}{4\pi}C_{n00}^{i}C_{n^{\prime}00}^{i}\sum_{j,k}^{N_{a}}\Bigg{(}\frac{4}{I_{0}(\cfrac{1}{\sigma^{2}})}\Bigg{)}^{3}j_{l}\Big{(}\frac{1}{\sigma^{2}}\Big{)}j_{l_{1}}\Big{(}\frac{1}{\sigma^{2}}\Big{)}j_{l_{2}}\Big{(}\frac{1}{\sigma^{2}}\Big{)}Y_{lm_{1}+m_{2}}({\bf\hat{e}}_{i})Y_{l_{1}m_{1}}^{*}({\bf\hat{e}}_{j})Y_{l_{2}m_{2}}^{*}({\bf\hat{e}}_{k})$
(48)
$\displaystyle\times\sqrt{\frac{8\pi^{2}}{2l+1}}<l,m_{1}+m_{2}|l_{1}m_{1};l_{2}m_{2}>$

The kernel

$K(\rho_{i},\rho^{\prime}_{i})=\sum_{nn^{\prime}ll_{1}l_{2},m_{1}m_{2}}C_{nn^{\prime}ll_{1}l_{2},m_{1}m_{2}}^{i*}C_{nn^{\prime}ll_{1}l_{2},m^{\prime}_{1}m^{\prime}_{2}}^{{}^{\prime}i}$
(49)

Another kind of ”three body term” with a bit more strict constraints on the
relative atomic position can be defined by the independent rotation invariant
symmetry of ($r_{ij}$,$r_{ik}$) and ($e_{i}$,$e_{j}$, $e_{k}$), which lies in-
between the three body terms of magnetic dipole-dipole interaction and
exchange interaction. This three body term should be more efficient in
describing the non-collinear effect which originates from the interaction
between neighbouring spins.

The overlap of the new “three body term” descriptor

$\displaystyle S(R_{1},R_{2})$ $\displaystyle=\iiint
r_{1}^{2}dr_{1}r_{2}^{2}dr_{2}d{\bf\hat{e}}_{1}d{\bf\hat{e}}_{2}d{\bf\hat{e}}_{3}\rho(r_{1},r_{2},{\bf\hat{e}}_{1},{\bf\hat{e}}_{2},{\bf\hat{e}}_{3})\rho^{\prime}(R_{1}r_{1},R_{1}r_{2},R_{2}{\bf\hat{e}}_{1},R_{2}{\bf\hat{e}}_{2},R_{2}{\bf\hat{e}}_{3})$
(50)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!=\sum_{\begin{subarray}{c}m,m_{1},m_{2}\\\
m^{\prime},m^{\prime}_{1},m^{\prime}_{2}\\\
l,l_{1},l_{2}\end{subarray}}\\!\\!\\!\\!\\!\\!C_{LL_{1}MM_{1}}^{*}C^{\prime}_{LL_{1}M^{\prime}M^{\prime}_{1}}C_{l{l_{1}}{l_{2}}m_{1}m_{2}m_{3}}^{*}C^{\prime}_{ll_{1}l_{2}m^{\prime}m^{\prime}_{1}m^{\prime}_{2}}D_{MM^{\prime}}^{L}(R_{1})D_{M_{1}M^{\prime}_{1}}^{L_{1}}(R_{1})D_{mm^{\prime}}^{l}(R_{2})D_{m_{1}m^{\prime}_{1}}^{l_{1}}(R_{2})D_{m_{2}m^{\prime}_{2}}^{l_{2}}(R_{2})$

The kernel for the three-body term descriptor

$\displaystyle K(\rho,\rho^{\prime})$ $\displaystyle=\iint
S(R_{1},R_{2})dR_{1}dR_{2}$ (51)
$\displaystyle=\sum_{\begin{subarray}{c}m,m_{1},m_{2}\\\
m^{\prime},m^{\prime}_{1},m^{\prime}_{2}\\\
l,l_{1},l_{2}\end{subarray}}C_{LL_{1}MM_{1}}^{*}C^{\prime}_{LL_{1}M^{\prime}M^{\prime}_{1}}C_{ll_{1}l_{2}mm_{1}m_{2}}^{*}C^{\prime}_{ll_{1}l_{2}m^{\prime}m^{\prime}_{1}m^{\prime}_{2}}$
$\displaystyle\times$ $\displaystyle\int
D_{MM^{\prime}}^{L}(R_{1})D_{M_{1}M^{\prime}_{1}}^{L_{1}}(R_{1})dR_{1}\int
D_{mm^{\prime}}^{l}(R)D_{m_{1}m^{\prime}_{1}}^{l_{1}}(R)D_{m_{2}m^{\prime}_{2}}^{l_{2}}(R)dR_{2}$
$\displaystyle=\sum_{\begin{subarray}{c}m_{1},m_{2}\\\
m^{\prime}_{1},m^{\prime}_{2}\\\
l_{1},l_{2}\end{subarray}}C_{l_{1}l_{2}m_{1}m_{2}}^{*}C^{\prime}_{l_{1}l_{2}m^{\prime}_{1}m^{\prime}_{2}}(-1)^{m_{1}-m^{\prime}_{1}}\int
D_{-m_{1},-m^{\prime}_{1}}^{l_{1}}(R)^{*}(R)D_{m_{2}m^{\prime}_{2}}^{l_{2}}(R)dR$
$\displaystyle\times$ $\displaystyle\sum_{\begin{subarray}{c}m,m_{1},m_{2}\\\
m^{\prime},m^{\prime}_{1},m^{\prime}_{2}\\\
l,l_{1},l_{2}\end{subarray}}C_{ll_{1}l_{2}mm_{1}m_{2}}^{*}C^{\prime}_{ll_{1}l_{2}m^{\prime}m^{\prime}_{1}m^{\prime}_{2}}(-1)^{m-m^{\prime}}\int
D_{-m,-m^{\prime}}^{l}(R)^{*}D_{m_{1}m^{\prime}_{1}}^{l_{1}}(R)D_{m_{2}m^{\prime}_{2}}^{l_{2}}(R)dR$
$\displaystyle=\sum_{\begin{subarray}{c}M_{1},M^{\prime}_{1},L\end{subarray}}\frac{8\pi^{2}(-1)^{M_{1}+M^{\prime}_{1}}}{2L+1}C_{LM_{1},-M_{1}}^{*}C^{\prime}_{LM^{\prime}_{1},-M^{\prime}_{1}}\sum_{\begin{subarray}{c}m,m_{1},m_{2}\\\
m^{\prime},m^{\prime}_{1},m^{\prime}_{2}\\\
l,l_{1},l_{2}\end{subarray}}C_{ll_{1}l_{2}mm_{1}m_{2}}^{*}C^{\prime}_{ll_{1}l_{2}m^{\prime}m^{\prime}_{1}m^{\prime}_{2}}$
$\displaystyle\times\frac{8\pi^{2}(-1)^{m+m^{\prime}}}{2l+1}<l,-m|l_{1}m_{1};l_{2}m_{2}><l,-m^{\prime}|l_{1}m^{\prime}_{1};l_{2}m^{\prime}_{2}>$

$\displaystyle\rho_{i}(r_{1})=\frac{1}{\sqrt{4\pi}}\sum_{n=1}^{N_{R}}C_{n}^{i}\chi_{nl}(r_{1})\;\;\;\;;\;\;\;\;\rho_{i}(r_{2})=\frac{1}{\sqrt{4\pi}}\sum_{n=1}^{N_{R}}C_{n}^{i}\chi_{nl}(r_{2})$
(52) $\displaystyle
C_{ll_{1}l_{2},m_{1}m_{2}}^{i}=\sum_{j}^{N_{a}}\Bigg{(}\frac{4}{I_{0}(\cfrac{1}{\sigma^{2}})}\Bigg{)}^{3}j_{l}\Big{(}\frac{1}{\sigma^{2}}\Big{)}j_{l_{1}}\Big{(}\frac{1}{\sigma^{2}}\Big{)}j_{l_{2}}\Big{(}\frac{1}{\sigma^{2}}\Big{)}Y_{lm_{1}+m_{2}}({\bf\hat{e}}_{i})Y_{l_{1}m_{1}}^{*}({\bf\hat{e}}_{j})Y_{l_{2}m_{2}}^{*}({\bf\hat{e}}_{k})$
$\displaystyle\times\sqrt{\frac{8\pi^{2}}{2l+1}}<l,m_{1}+m_{2}|l_{1}m_{1};l_{2}m_{2}>$

Finally, the coefficients for the kernel become

$\displaystyle
C_{nn^{\prime}ll_{1}l_{2},m_{1}m_{2}}^{i}=\frac{1}{4\pi}C_{n00}^{i}C_{n^{\prime}00}^{i}\sum_{j,k}^{N_{a}}\Bigg{(}\frac{4}{I_{0}(\cfrac{1}{\sigma^{2}})}\Bigg{)}^{3}j_{l}\Big{(}\frac{1}{\sigma^{2}}\Big{)}j_{l_{1}}\Big{(}\frac{1}{\sigma^{2}}\Big{)}j_{l_{2}}\Big{(}\frac{1}{\sigma^{2}}\Big{)}Y_{lm_{1}+m_{2}}({\bf\hat{e}}_{i})Y_{l_{1}m_{1}}^{*}({\bf\hat{e}}_{j})Y_{l_{2}m_{2}}^{*}({\bf\hat{e}}_{k})$
(53)
$\displaystyle\times\sqrt{\frac{8\pi^{2}}{2l+1}}<l,m_{1}+m_{2}|l_{1}m_{1};l_{2}m_{2}>$

where the kernel

$K(\rho_{i},\rho^{\prime}_{i})=\sum_{nn^{\prime}ll_{1}l_{2},m_{1}m_{2}}C_{nn^{\prime}ll_{1}l_{2},m_{1}m_{2}}^{i*}C_{nn^{\prime}ll_{1}l_{2},m^{\prime}_{1}m^{\prime}_{2}}^{{}^{\prime}i}$
(54)

##  Appendix B The lattice constant dependence of total energy and local
magnetic moment of bcc Fe

Figure 6:  The total energy and local magnetic moment of bcc Fe as a function
of lattice constant for FM, AFM and NM states calculated with LDA.  Figure 7:
The total energy and local magnetic moment of bcc Fe as a function of lattice
constant for FM, AFM and NM states calculated with GGA.

## References

  * Car and Parrinello [1985] R. Car and M. Parrinello, Unified approach for molecular dynamics and density-functional theory, Phys. Rev. Lett. 55, 2471 (1985).
  * Behler and Parrinello [2007] J. Behler and M. Parrinello, Generalized neural-network representation of high-dimensional potential-energy surfaces, Phys. Rev. Lett. 98, 146401 (2007).
  * Hohenberg and Kohn [1964] P. Hohenberg and W. Kohn, Inhomogeneous electron gas, Phys. Rev. 136, B864 (1964).
  * Kohn and Sham [1965] W. Kohn and L. J. Sham, Self-consistent equations including exchange and correlation effects, Phys. Rev. 140, A1133 (1965).
  * Other and Els [1978] A. N. Other and S. W. Els, We need a reference here, Adv. Phys. 27, 799 (1978).
  * Behler [2016] J. Behler, Perspective: Machine learning potentials for atomistic simulations, J. Chem. Phys. 145, 170901 (2016).
  * Bartók  _et al._ [2010] A. P. Bartók, M. C. Payne, R. Kondor, and G. Csányi, Gaussian Approximation Potentials: The Accuracy of Quantum Mechanics, without the Electrons, Phys. Rev. Lett. 104, 136403 (2010).
  * Jinnouchi  _et al._ [2019a] R. Jinnouchi, F. Karsai, and G. Kresse, On-the-fly machine learning force field generation: Application to melting points, Phys. Rev. B 100, 014105 (2019a).
  * Bartók  _et al._ [2013] A. P. Bartók, R. Kondor, and G. Csányi, On representing chemical environments, Phys. Rev. B 87, 184115 (2013).
  * Jinnouchi  _et al._ [2019b] R. Jinnouchi, J. Lahnsteiner, F. Karsai, G. Kresse, and M. Bokdam, Phase Transitions of Hybrid Perovskites Simulated by Machine-Learning Force Fields Trained on the Fly with Bayesian Inference, Phys. Rev. Lett. 122, 225701 (2019b).
  * Miltat  _et al._ [2002] J. Miltat, G. Albuquerque, and A. Thiaville, An Introduction to Micromagnetics in the Dynamic Regime (Springer, Berlin, 2002) Chap. 1, pp. 1–34.
  * Tserkovnyak  _et al._ [2005] Y. Tserkovnyak, A. Brataas, G. E. W. Bauer, and B. I. Halperin, Nonlocal magnetization dynamics in ferromagnetic nanostructures, Rev. Mod. Phys. 77, 1375 (2005).
  * Skubic  _et al._ [2008] B. Skubic, J. Hellsvik, L. Nordström, and O. Eriksson, A method for atomistic spin dynamics simulations: implementation and examples, J. Phys.: Condens. Matter 20, 315203 (2008).
  * Eriksson  _et al._ [2017] O. Eriksson, A. Bergman, L. Bergqvist, and J. Hellsvik,  _Atomistic Spin Dynamics: Foundations and Applications_ (Oxford University Press, Oxford, 2017).
  * Bartók and Csányi [2015] A. P. Bartók and G. Csányi, Gaussian approximation potentials: A brief tutorial introduction, International Journal of Quantum Chemistry 115, 1051 (2015).
  * Heine  _et al._ [2021] M. Heine, O. Hellman, and D. Broido, Temperature-dependent renormalization of magnetic interactions by thermal, magnetic, and lattice disorder from first principles, Phys. Rev. B 103, 184409 (2021).
  * Pajda  _et al._ [2001] M. Pajda, J. Kudrnovskỳ, I. Turek, V. Drchal, and P. Bruno, Ab initio calculations of exchange interactions, spin-wave stiffness constants, and Curie temperatures of Fe, Co, and Ni, Phys. Rev. B 64, 174402 (2001).
  * Born and Oppenheimer [1927] M. Born and R. J. Oppenheimer, Zur Quantentheorie der Molekeln (On the Quantum Theory of Molecules), Annalen der Physik 84, 457 (1927).
  * Antropov  _et al._ [1995] V. P. Antropov, M. I. Katsnelson, M. van Schilfgaarde, and B. N. Harmon, Ab Initio spin dynamics in magnets, Phys. Rev. Lett. 75, 729 (1995).
  * Antropov  _et al._ [1996] V. P. Antropov, M. I. Katsnelson, B. N. Harmon, M. van Schilfgaarde, and D. Kusnezov, Spin dynamics in magnets: Equation of motion and finite temperature effects, Phys. Rev. B 54, 1019 (1996).
  * White [1983] R. M. White,  _Quantum theory of magnetism_, 2nd ed., edited by M. Carona, P. Fulde, K. von Klitzing, and H.-J. Queisser, Springer Series in Solid-State Sciences, Vol. 32 (Springer, Berlin, 1983).
  * Paige  _et al._ [1984] D. M. Paige, B. Szpunar, and B. K. Tanner, The magnetocrystalline anisotropy of cobalt, J. Magn. Magn. Mater. 44, 239 (1984).
  * Escudier [1975] P. Escudier, L’anisotropie de l’aimantation: un paramètre important de l’étude de l’anisotropie magnétocristalline (magnetization anisotropy - important parameter for study of magnetocrystalline anisotropy) (1975) pp. 125–173.
  * Chen  _et al._ [2013a] G. Chen, J. Zhu, A. Quesada, J. Li, A. T. N’Diaye, Y. Huo, T. P. Ma, Y. Chen, H. Y. Kwon, C. Won, Z. Q. Qiu, A. K. Schmid, and Y. Z. Wu, Novel chiral magnetic domain wall structure in Fe/Ni/Cu (001) films, Phys. Rev. Lett. 110, 177204 (2013a).
  * Yang  _et al._ [2018] H. Yang, G. Chen, A. A. C. Cotta, A. T. N’Diaye, S. A. Nikolaev, E. A. Soares, W. A. A. Macedo, K. Liu, A. K. Schmid, A. Fert, and M. Chshiev, Significant Dzyaloshinskii-Moriya interaction at graphene-ferromagnet interfaces due to the Rashba effect, Nat. Mater. 17, 605 (2018).
  * Chen  _et al._ [2013b] G. Chen, T. Ma, A. T. N’Diaye, H. Kwon, C. Won, Y. Wu, and A. K. Schmid, Tailoring the chirality of magnetic domain walls by interface engineering, Nat. Commun. 4, 2671 (2013b).
  * Kikuchi [1951] R. Kikuchi, A theory of cooperative phenomena, Phys. Rev. 81, 988 (1951).
  * Pozdnyakov  _et al._ [2020] S. N. Pozdnyakov, M. J. Willatt, A. P. Bartók, C. Ortner, G. Csányi, and M. Ceriotti, Incompleteness of Atomic Structure Representations, Phys. Rev. Lett. 125, 166001 (2020).
  * Kaufmann and Baumeister [1989] K. Kaufmann and W. Baumeister, Single-centre expansion of Gaussian basis functions and the angular decomposition of their overlap integrals, J. Phys. B: Atomic, Molecular and Optical Physics 22, 1 (1989).
  * foo [a] (a), In the QUIP code from Cambridge University, the radial parts are expanded in a set of equally spaced Gaussian functions while in the vasp ML package the radial basis functions are normalised spherical Bessel functions. The Gaussian basis functions need to be orthogonalised. The normalised spherical Bessel functions are mutually orthogonal. .
  * foo [b] (b), The Born-Oppenheimer approximation makes use of the different time scales of electronic and nuclear motion determined by the electron and nuclear masses to integrate out the electronic degrees of freedom and calculate the PES that determines nuclear motion [18]. Because there is no moment of inertia associated with the electron spin, we must use the characteristic frequencies of different spin interactions to establish an adiabatic approximation [19, *Antropov:prb96]. .
  * Wigner [1959] E. P. Wigner, Group theory and its application to the quantum mechanics of atomic spectra, in  _Group Theory and its application to the quantum mechanics of atomic spectra_ (Academic Press, 1959).
  * Ma and Dudarev [2015] P.-W. Ma and S. L. Dudarev, Constrained density functional for noncollinear magnetism, Phys. Rev. B 91, 054420 (2015).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Daalderop  _et al._ [1990] G. H. O. Daalderop, P. J. Kelly, and M. F. H. Schuurmans, First-principles calculation of the magnetocrystalline anisotropy energy of iron, cobalt and nickel, Phys. Rev. B 41, 11919 (1990).
  * Daalderop  _et al._ [1992] G. H. O. Daalderop, P. J. Kelly, and F. J. A. den Broeder, Prediction and confirmation of perpendicular magnetic anisotropy in Co/Ni multilayers, Phys. Rev. Lett. 68, 682 (1992).
  * Daalderop  _et al._ [1994] G. H. O. Daalderop, P. J. Kelly, and M. F. H. Schuurmans, Magnetic anisotropy of a free-standing Co monolayer and of multilayers which contain Co monolayers, Phys. Rev. B 50, 9989 (1994).
  * Drautz [2019] R. Drautz, Atomic cluster expansion for accurate and transferable interatomic potentials, Phys. Rev. B 99, 014104 (2019).
  * Rinaldi  _et al._ [2023] M. Rinaldi, M. Mrovec, A. Bochkarev, Y. Lysogorskiy, and R. Drautz, Non-collinear magnetic atomic cluster expansion for iron, arXiv:2305.15137 10.48550/arXiv.2305.15137 (2023).
  * Domina  _et al._ [2022] M. Domina, M. Cobelli, and S. Sanvito, Spectral neighbor representation for vector fields: Machine learning potentials including spin, Phys. Rev. B 105, 214439 (2022).
  * Suzuki  _et al._ [2023] M.-T. Suzuki, T. Nomoto, E. V. Morooka, Y. Yanagi, and H. Kusunose, High-performance descriptor for magnetic materials: Accurate discrimination of magnetic structure, Phys. Rev. B 108, 014403 (2023).