# Thermally Averaged Magnetic Anisotropy Tensors via Machine Learning Based on
Gaussian Moments

Viktor Zaverkin  Julia Netz  Fabian Zills  Andreas Köhn  Johannes Kästner
<EMAIL_ADDRESS>Institute for Theoretical Chemistry,
University of Stuttgart, Pfaffenwaldring 55, 70569 Stuttgart, Germany

###### Abstract

We propose a machine learning method to model molecular tensorial quantities,
namely the magnetic anisotropy tensor, based on the Gaussian-moment neural-
network approach. We demonstrate that the proposed methodology can achieve an
accuracy of 0.3–0.4 $\text{cm}^{-1}$ and has excellent generalization
capability for out-of-sample configurations. Moreover, in combination with
machine-learned interatomic potential energies based on Gaussian moments, our
approach can be applied to study the dynamic behavior of magnetic anisotropy
tensors and provide a unique insight into spin-phonon relaxation.

###### keywords:

American Chemical Society, LaTeX

IR,NMR,UV

##  1  Introduction

There is an ongoing interest in the investigation of magnetic properties of
transition metal complexes and their possible applications as single-molecule
magnets (SMMs), molecular quantum bits, and spintronic devices.1, 2, 3, 4 One
of the most important factors influencing the magnetic properties is the
magnetic anisotropy in the ground spin state.5, 2 Tailoring promising
complexes in a way to achieve a large barrier for magnetic relaxation and
minimizing other decay pathways6, 7, 8, 9 requires a detailed understanding of
the underlying principles that determine the structure-property relationships.

While sampling based on first-principles methods can provide some insight, it
is usually limited by its high computational cost. Approximate methods allow
for simulations even taking into account the large size of conformational and
chemical space of interest. Machine learning approaches with their ability to
learn any complex non-linear relationship between a structure and its related
property and with their generalization capability are perfect candidates for
the specific task of modeling the magnetic anisotropy. During the last
decades, machine learning methods have been gaining importance in several
fields in computational chemistry allowing for, e.g., the construction of
machine-learned force fields and prediction of vibrational spectra.10, 11, 12,
13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30 The
ability of neural networks (NNs) to interpolate any non-linear functional
relationship 31 promoted their broad application in computational chemistry
and materials science. NNs were initially applied to represent potential
energy surfaces (PESs) of small atomistic systems10, 11 and were later
extended to high-dimensional systems.13 Once trained, the computational cost
of machine-learned potentials (MLPs) based on NNs is independent of the number
of data points used for training.

Case studies of the application of machine learning algorithms to study the
magnetic properties of molecules and materials are somewhat rare. Recently, an
approach for the construction of machine-learned interatomic potentials which
are capable of reproducing both vibrational and magnetic degrees of freedom
was proposed32 and applied to bcc iron. Closer to the present work are the
investigations presented in Refs. 33, 34, where tensorial properties such as
the zero-field splitting (ZFS) tensor $\mathbf{D}$ and the Zeeman-splitting
(ZS) tensor $\mathbf{g}$ were modeled by machine learning approaches.

In this work, we build upon the methodology developed in our group for
constructing efficient and accurate interatomic potentials, referred to as
Gaussian moment neural network (GM-NN). 24, 30 The GM-NN uses neural networks
(NNs) to map novel symmetry-preserving local atomic descriptors, Gaussian
moments (GMs), to auxiliary atomic quantities and includes both the geometric
and alchemical information about the atomic species of both the central and
neighbor atoms. For all atomic contributions, only a single NN has to be
trained, in contrast to using an individual NN for each species as frequently
done in the literature.13, 15 To allow for efficient training on properties
such as $\mathbf{D}$ tensors we introduce a novel approach for the encoding of
relevant invariances into the output of the respective machine learning model.
Moreover, we propose a neural network architecture for the task of modeling
the magnetic anisotropy of SMMs.

To assess the quality of the surrogate models obtained employing the proposed
approach we thoroughly benchmark the predictive accuracy on three promising
SMMs: [Co(N2S2O4C8H10)2]2-, 6 [Fe(tpa)Ph]-, 7, 35 and [Ni(HIM2-py)2NO3]+
complexes. 8 Finally, we use the [Co(N2S2O4C8H10)2]2- complex as a case study
to investigate the applicability of the proposed approach to the dynamics of
the $\mathbf{D}$ tensor and examine its ability to model spin-phonon coupling
effects.

In this work, we employ the open-source package Gaussian moment neural network
(GM-NN) for constructing MLPs and applying them in atomistic simulations. The
GM-NN source code is available free-of-charge from gitlab.com/zaverkin_v/gmnn.

##  2  Machine-Learning Models

This section first introduces the representation of an atomistic system used
to encode molecular and solid-state structures suitable for machine learning
(ML) models. Second, we describe the approach to machine learning modeling the
magnetic anisotropy tensors including the neural network architecture and its
training, based on our previous work on Gaussian moments. 24, 30 Finally, we
outline the construction of the machine-learned interatomic potential energy
used in Section 4.2 for molecular dynamics simulations of
[Co(N2S2O4C8H10)2]2-.

###  2.1  Molecular Representation

A molecular or solid-state system is defined by its Cartesian coordinates
$\mathbf{r}_{i}\in\mathbb{R}^{3}$ and its atomic numbers $Z_{i}$ which, in the
following discussion, are combined to
$\mathcal{S}=\\{\mathbf{r}_{i},Z_{i}\\}_{i=1}^{N_{\mathrm{at}}}$ for
simplicity. In the community of neural network (NN) model chemistry, it is
usual to approximate properties by a sum of their atomic contributions. For
example, the total energy of a system can be approximated13 by a sum of atomic
energies $\mathcal{E}_{i}$

$\mathcal{E}\left(\mathcal{S},\boldsymbol{\theta}\right)=\sum_{i=1}^{N_{\mathrm{at}}}\mathcal{E}_{i}\left(\mathbf{G}_{i},\boldsymbol{\theta}\right).$
(1)

Total charge and the dipole moments can be defined via the atomic point
charges $q_{i}$

$\begin{split}&Q_{\mathrm{tot}}\left(\mathcal{S},\boldsymbol{\theta}\right)=\sum_{i=1}^{N_{\mathrm{at}}}q_{i}\left(\mathbf{G}_{i},\boldsymbol{\theta}\right),\\\
&\boldsymbol{\mu}\left(\mathcal{S},\boldsymbol{\theta}\right)=\sum_{i=1}^{N_{\mathrm{at}}}q_{i}\left(\mathbf{G}_{i},\boldsymbol{\theta}\right)\mathbf{r}_{i},\end{split}$
(2)

where each $\mathcal{E}_{i}$ and $q_{i}$ are given as the output of an
artificial NN, and all relevant invariances are encoded in the local
representation $\mathbf{G}_{i}(\mathcal{S},\boldsymbol{\beta})$ with
$\boldsymbol{\beta}$ being trainable parameters. Thus, the main purpose of
machine learning is to find such parameters $\boldsymbol{\theta}$ that the
mapping $f:\mathcal{S}\mapsto\\{E,Q_{\mathrm{tot}},\boldsymbol{\mu},\dots\\}$
is as close as possible to the reference values,
$\\{E^{\mathrm{ref}},Q_{\mathrm{tot}}^{\mathrm{ref}},\boldsymbol{\mu}^{\mathrm{ref}},\dots\\}$.

The most challenging aspect of a machine learning model applied to a
physicochemical problem is the definition of an appropriate representation of
an atomistic system. In this work, we employ Gaussian moment representation
which defines the local atomic environment as24, 30

$\begin{split}\boldsymbol{\Psi}_{i,L,s}=\sum_{j\neq
i}R_{Z_{i},Z_{j},s}\left(r_{ij},\boldsymbol{\beta}\right)\hat{\mathbf{r}}_{ij}^{\otimes
L},\end{split}$ (3)

where we use the radial and angular components of the atomic distance vector
$\mathbf{r}_{ij}=\mathbf{r}_{i}-\mathbf{r}_{j}$, i.e.
$r_{ij}=\lVert\mathbf{r}_{ij}\rVert_{2}$ and
$\hat{\mathbf{r}}_{ij}=\mathbf{r}_{ij}/r_{ij}$, as inputs to nonlinear radial
functions $R_{Z_{i},Z_{j},s}\left(r_{ij},\boldsymbol{\beta}\right)$ with
trainable parameters $\boldsymbol{\beta}$ and the $L$-fold tensor product
$\hat{\mathbf{r}}_{ij}^{\otimes
L}=\hat{\mathbf{r}}_{ij}\otimes\cdots\otimes\hat{\mathbf{r}}_{ij}$,
respectively. $Z_{i}$ and $Z_{j}$ correspond to the nuclear charges of the
central atom $i$ and its atomic neighbors $j$. As nonlinear radial functions
we employ a weighted sum of Gaussian functions, rescaled by the cosine cutoff
function, 13 see Ref. 24, 30. The parameters $\boldsymbol{\beta}$ are
optimized during training similar to other parameters of an atomistic NN by
minimization of a loss function. A rotationally invariant representation is
obtained by computing full tensor contractions as described in Ref. 24, 30.

As a final remark on the representation of an atomistic system, we want to
discuss the possibility of learning different physicochemical properties by a
single ML model. In general, it should be possible, given that the
corresponding properties are available for the same atomistic system.
Unfortunately, in this work, the energy and atomic forces are calculated for
the full periodic system while the magnetic anisotropy tensors are computed
for isolated molecules cut from them, see Fig. 1. That implies that both
atomistic systems provide different sets of input features and, therefore,
cannot directly be combined. For this reason, in the following sections we
describe the construction of two separate ML models for learning $\mathbf{D}$
tensors and potential energy surfaces.

###  2.2  The $\mathbf{D}$ Tensor

Mediated by spin-orbit coupling, the components of spin multiplets (with spin
quantum numbers $S\geq 1$) split characteristically even in the absence of an
external magnetic field (zero-field splitting, ZFS). This effect is usually
described by a phenomenological spin Hamiltonian36

$H_{\mathrm{ZFS}}=\hat{\mathbf{S}}\cdot\mathbf{D}\cdot\hat{\mathbf{S}}$ (4)

where $\hat{\mathbf{S}}$ is a (pseudo) spin operator and D is a $3\times 3$
symmetric, traceless tensor, usually called ZFS tensor or D tensor. The
anisotropy parameters $D$ and $E$ (often simply referred to as $D$ and $E$
values) can be then obtained from the eigenvalues of $\mathbf{D}$ ($X$, $Y$,
and $Z$) as

$\begin{split}&D=\frac{3}{2}Z,\\\ &E=\frac{1}{2}(X-Y),\end{split}$ (5)

where $Z$ corresponds to the eigenvalue of $\mathbf{D}$ with largest absolute
value. The normalized eigenvectors of D are the anisotropy axes of the system.

####  2.2.1 Incorporation of Symmetries into the ML Framework

The aim of this work is to predict $\mathbf{D}$ tensors using artificial
neural networks (NNs). The $\mathbf{D}$ tensor is a (1) traceless (2)
symmetric tensor, which implies that
$\operatorname{Tr}\mathbf{D}=\sum_{i}D_{ii}=0$ and $D_{ij}=D_{ji}$.
Additionally, it (3) transforms under rotation as
$\tilde{\mathbf{D}}=\mathbf{R}\mathbf{D}\mathbf{R}^{T}$, where $\mathbf{R}$ is
an orthogonal matrix, but is (4) invariant to translations. These properties
have to be encoded into the respective machine learning model to allow for
efficient training. Moreover, unlike the energy, the $\mathbf{D}$ tensor
cannot be reduced to a natural scalar atomic property which makes its
prediction more difficult using the atomistic neural networks.

To the best of our knowledge, there are only two studies in the literature of
modeling tensorial properties such as $\mathbf{g}$ and $\mathbf{D}$
tensors.33, 34 In Ref. 33, where ridge regression was used to fit $\mathbf{D}$
tensors, the rotational equivariance was imposed by requiring the regression
coefficients to transform as spherical tensors, i.e. employing the sum over
the Wigner matrices corresponding to rigid rotations. In Ref. 34, where
artificial neural networks were used to fit $\mathbf{g}$ tensors, a rotation
operator $\Lambda$ was defined by, e.g., the Kabsh algorithm to impose the
desired symmetries. In this work, we present a more rigorous and efficient way
of encoding the symmetries (1)–(4) into machine learning algorithms.

We approach the modeling of zero-field splitting tensors by introducing a
fictitious atomic quantity $m_{i}$ assigned to each atom $i$. A tensor which
satisfies (1)–(3) can be defined using the tensor product of Cartesian
vectors, similar to the quadrupole moment, employing the aforementioned atomic
quantity $m_{i}$

$\mathbf{D}\left(\mathcal{S},\boldsymbol{\theta}\right)=\sum_{i=1}^{N_{\mathrm{at}}}m_{i}\left(\mathbf{G}_{i},\boldsymbol{\theta}\right)\cdot\left(3\mathbf{r}_{i}\otimes\mathbf{r}_{i}-\|\mathbf{r}_{i}\|_{2}^{2}\mathbb{I}_{3}\right)~{},$
(6)

where $\otimes$ is the tensor product, $\|\mathbf{r}_{i}\|_{2}$ is the length
of the respective Cartesian vector, and $\mathbb{I}_{3}$ is a $3\times 3$
identity matrix.

Unfortunately, the above expression violates translational invariance (4) of
the respective tensorial property. This issue can be resolved by shifting the
coordinate system by, e.g.,
$\bar{\mathbf{r}}=1/N_{\mathrm{at}}\sum_{i}\mathbf{r}_{i}$, which is system
dependent. Note that the shift can be defined by any other procedure since its
aim is only to impose the translational invariance (4). However, we should
emphasize that we have found that using the pairwise distances
$\mathbf{r}_{ij}=\mathbf{r}_{i}-\mathbf{r}_{j}$, where $i$ can be selected to
be the metal atom, could be disadvantageous in terms of predictive accuracy
due to missing contributions from the central atom.

Using the respective shift vector $\bar{\mathbf{r}}$ results in a new position
vector $\mathbf{r}_{i}\rightarrow\mathbf{r}_{i}-\bar{\mathbf{r}}$.
Additionally, we have found that it is advantageous to use the normalized
vectors $\hat{\mathbf{r}}_{i}=\mathbf{r}_{i}/\|\mathbf{r}_{i}\|_{2}$. In
total, the resulting expression reads

$\mathbf{D}\left(\mathcal{S},\boldsymbol{\theta}\right)=\sum_{i=1}^{N_{\mathrm{at}}}m_{i}\left(\mathbf{G}_{i},\boldsymbol{\theta}\right)\cdot\left(3\hat{\mathbf{r}}_{i}\otimes\hat{\mathbf{r}}_{i}-\mathbb{I}_{3}\right)~{},$
(7)

where $m_{i}$ is predicted by an atomistic NN.

As a final remark we want to point out that any symmetry of a tensorial
property $\mathbf{P}$ can be modelled by the procedure proposed above if one
defines $\mathbf{P}=\sum_{i}m_{i}\mathbf{A}_{i}$, where $\mathbf{A}_{i}$ is a
tensor satisfying the symmetry of $\mathbf{P}$ and $m_{i}$ is a machine-
learned scalar value. In our example, we define
$\mathbf{A}_{i}=3\hat{\mathbf{r}}_{i}\otimes\hat{\mathbf{r}}_{i}-\mathbb{I}_{3}$
to impose (1)–(4), but if the respective property is, e.g, not traceless one
could use $\mathbf{A}_{i}=\hat{\mathbf{r}}_{i}\otimes\hat{\mathbf{r}}_{i}$.

####  2.2.2 Network Architecture and Training

We use a fully-connected feed-forward neural network consisting of two hidden
layers of the following functional form, similar to our previous work, 30

$\begin{split}y_{i}\left(\mathbf{G}_{i},\boldsymbol{\theta}\right)&=0.1\cdot\mathbf{b}^{(3)}+\frac{1}{\sqrt{d_{2}}}\mathbf{W}^{(3)}\phi\left(0.1\cdot\mathbf{b}^{(2)}+\right.\\\
&\quad\left.\frac{1}{\sqrt{d_{1}}}\mathbf{W}^{(2)}\phi\left(0.1\cdot\mathbf{b}^{(1)}+\frac{1}{\sqrt{d_{0}}}\mathbf{W}^{(1)}\mathbf{G}_{i}\right)\right)~{},\end{split}$
(8)

where $\mathbf{W}^{(l)}$ and $\mathbf{b}^{(l)}$ are weights and biases of the
respective layer $l$. As an input to the neural network, $\mathbf{G}_{i}$, we
use the recently proposed trainable local invariant representation based on
Gaussian moments.24, 30 The parameters $0.1$ and $1/\sqrt{d_{l}}$ correspond
to the so-called NTK parameterisation.37 We initialize weights of the fully-
connected part by drawing the respective entries from a normal distribution
with zero mean and unit variance. The trainable bias vectors are initialized
to zero. As an activation function we use the Swish/SiLU activation
function38, 39, 40 $\phi(x)=\alpha x/\left(1+\exp\left(-x\right)\right)$
multiplied by a scalar $\alpha$. We choose $\alpha\approx 1.6765$ such that
$\mathbb{E}_{x\sim\mathcal{N}\left(0,1\right)}\phi(x)^{2}=1$, i.e., the
activation function preserves the second moment if the input is standard
Gaussian. 41, 42, 43

In order to aid the training process, the output of the neural network can be
scaled and shifted by the standard deviation $\sigma$ and the mean $\mu$ of
the reference $\mathbf{D}^{\mathrm{ref}}$ tensor values, similar to the
scaling and shifting the atomic energy output. 24, 30 Note that we used for
the computation of $\sigma$ and $\mu$ only those elements of $D_{ij}$ which
satisfy $i\geq j$ to avoid double counting and excluded one diagonal element
since the respective tensor is traceless. The convergence of the model can be
improved even further by making these parameters trainable as well as
dependent on the atomic species, i.e. $\sigma_{Z_{i}}$ and $\mu_{Z_{i}}$. The
final output of the network reads

$m_{i}\left(\mathbf{G}_{i},\boldsymbol{\theta}\right)=y_{i}\left(\mathbf{G}_{i},\boldsymbol{\theta}\right)\sigma_{Z_{i}}+\mu_{Z_{i}}~{}.$
(9)

To train the neural network on reference values for
$\mathbf{D}^{\mathrm{ref}}$ tensors, we minimize the following loss function

$\begin{split}\mathcal{L}_{\mathbf{D}}\left(\boldsymbol{\theta}\right)=\sum_{k=1}^{N_{\mathrm{Train}}}&\lVert\mathbf{D}_{k}^{\mathrm{ref}}-\mathbf{D}(\mathcal{S}_{k},\boldsymbol{\theta})\rVert_{2}^{2}~{},\end{split}$
(10)

to optimize the respective parameters of the trainable representation, fully-
connected neural network part as well as the parameters which scale and shift
the output of the neural network. Note that we train only on those elements of
$\left(\mathbf{D}\right)_{ij}$ tensor which satisfy $i\geq j$ and we define
$\mathbf{D}^{\mathrm{ref}}$ and $\mathbf{D}$ as

$\begin{split}&\mathbf{D}^{\mathrm{ref}}=\Big{(}D_{11}^{\mathrm{ref}},D_{12}^{\mathrm{ref}},D_{13}^{\mathrm{ref}},D_{22}^{\mathrm{ref}},D_{23}^{\mathrm{ref}},D_{33}^{\mathrm{ref}}\Big{)}^{T}\\\
&\mathbf{D}=\Big{(}D_{11},D_{12},D_{13},D_{22},D_{23},-\left(D_{11}+D_{22}\right)\Big{)}^{T}~{}.\end{split}$
(11)

To minimize the loss function in Equation (10) the Adam optimizer 44 with
hyper-parameters $\beta_{1}=0.9$, $\beta_{2}=0.999$, $\epsilon=10^{-7}$, and a
mini-batch of 32 structures is employed. Moreover, we allow for layer-wise
learning rates which decay linearly to zero by multiplying them with $(1-r)$
where $r=\mathrm{step}/\mathrm{max\\_step}$. Overall, we use an initial
learning rate of $0.02$ for the parameters of the fully connected layers,
$0.02$ for the parameters of the trainable GM representation, $0.025$ and
$0.025$ for the shift and scale parameters.

Throughout this work we used an architecture with an input dimension of
$d_{0}=360$, two hidden layers with $d_{1}=d_{2}=512$ hidden neurons, and an
output layer which has a single $d_{3}=1$ output neuron. Each model in Section
4 is trained for 1000 epochs. Overfitting was prevented by using the early
stopping technique. 45 After each epoch, the mean absolute errors (MAE) of the
$\mathbf{D}$ tensor was evaluated on the validation set. After training, the
model with the minimal MAE on the validation set was selected for further
applications. While the selected hyper-parameters worked reasonably well on
the selected systems for us, we want to emphasize that other trade-offs
between the number of training epochs and the initial learning rates can be
achieved.

Note that to compute machine-learned $\mathbf{D}$ tensors during an MD
simulation we interfaced our approach with the ASE package (v. 3.21.0). 46 For
tracking the accuracy we employed the query-by-committee (QbC) approach 47
during MD simulations. For this purpose, we trained a committee of 3 models on
the same split of the data set but using randomly initialized parameters and
reported the obtained uncertainty

$\sigma_{\mathrm{ens}}\left(\mathcal{S}\right)=\sqrt{\frac{1}{N_{\mathrm{ens}}}\sum_{i=1}^{N_{\mathrm{ens}}}\left(y_{i}\left(\mathcal{S}\right)-\bar{y}\left(\mathcal{S}\right)\right)^{2}}~{},$
(12)

where $N_{\mathrm{ens}}$ is the number of models in the committee, i.e.
$N_{\mathrm{ens}}=3$.
$\bar{y}\left(\mathcal{S}\right)=1/N_{\mathrm{ens}}\sum_{i=1}^{N_{\mathrm{ens}}}y_{i}\left(\mathcal{S}\right)$
is the mean of the property prediction (energy, atomic force element, or
$\mathbf{D}$-tensor element, respectively) over the committee. All models were
trained within the Tensorflow framework 48 on an NVIDIA Tesla V100-SXM-32GB
GPU. The training of an ensemble of 3 models for 1000 epochs took from 4 min
(100 [Co(N2S2O4C8H10)2]2- structures) to 3 h (2900 [Ni(HIM2-py)2NO3]+
structures).

###  2.3  Potential Energy

Thermal averaging of magnetic anisotropy tensors requires an interatomic
potential that fulfills two main premises. The underlying model has to produce
sufficiently accurate atomic forces for molecular dynamics (MD) simulations,
i.e., comparable to the level of theory employed for the generation of the
training data, and allow for efficient computations, i.e., comparable to
empirical force fields. Machine-learning algorithms have found a broad
application in computational chemistry since they satisfy both conditions
provided a suitable molecular descriptor that encodes structural and
alchemical information.

As for the $\mathbf{D}$ tensors, we employed the Gaussian moment neural
network (GM-NN) approach24, 30 to construct the PES. Again, we use a neural
network with an input dimension of $d_{0}=360$, two hidden layers with
$d_{1}=d_{2}=512$ hidden neurons, and an output layer which has a single
$d_{3}=1$ output neuron.

To train the GM-NN model the combined loss function

$\begin{split}\mathcal{L}_{\mathcal{E},\mathbf{F}}\left(\boldsymbol{\theta}\right)=\sum_{k=1}^{N_{\mathrm{Train}}}&\Bigg{[}\lambda_{\mathcal{E}}\lVert\mathcal{E}_{k}^{\mathrm{ref}}-\mathcal{E}(\mathcal{S}_{k},\boldsymbol{\theta})\rVert_{2}^{2}+\\\
&\quad\frac{\lambda_{F}}{3N_{\mathrm{at}}^{(k)}}\sum_{i=1}^{N_{\mathrm{at}}^{(k)}}\lVert\mathbf{F}_{i,k}^{\mathrm{ref}}-\mathbf{F}_{i}\left(\mathcal{S}_{k},\boldsymbol{\theta}\right)\rVert_{2}^{2}\Bigg{]}~{},\end{split}$
(13)

is minimized. Here, $N_{\mathrm{at}}^{(k)}$ is the number of atoms in the
respective structure. The reference values for the energy and atomic force are
denoted by $\mathcal{E}_{k}^{\mathrm{ref}}$ and
$\mathbf{F}_{i,k}^{\mathrm{ref}}$, respectively. The parameters
$\lambda_{\mathcal{E}}$ and $\lambda_{F}$ were set to $1~{}\text{au}$ and
$12N_{\mathrm{at}}^{(k)}$ au Å2, respectively. The network was trained using
the Adam optimizer 44 with $32$ molecules per mini-batch. The layer-wise
learning rate was set to $0.03$ for the parameters of the fully connected
layers, $0.02$ for the trainable representation, 0.05 and 0.001 for the shift
and scale parameters of atomic energies respectively. All learning rates we
allowed to decay linearly to zero. For more information, see Ref. 30.

To run MD simulations with the machine-learned potentials (MLPs) we interfaced
the GM-NN approach with the ASE package (v. 3.21.0). 46 For tracking the
accuracy of MLPs we employed the query-by-committee (QbC) approach 47 during
MD simulations. For this purpose, we trained a committee of 3 models on the
same split of the data set but using randomly initialized parameters and
reported on the attained uncertainty, see Equation (12). All models were
trained within the Tensorflow framework 48 on an NVIDIA Tesla V100-SXM-32GB
GPU. The training of an ensemble of 3 models for 1000 epochs took about 40 h
(3100 [Co(N2S2O4C8H10)2]2- periodic structures).

##  3  Test Systems and Computational Details

In this section, we describe the generation of the data used to construct
machine-learned interatomic potentials as well as machine learning models for
$\mathbf{D}$ anisotropy tensors and an overview of selected test systems.

###  3.1  Test System Description

We selected the following three promising candidates for SMMs (in fact even
single-ion magnets, SIMs): [Co(N2S2O4C8H10)2]2-, 6 [Fe(tpa)Ph]-, 7, 35 and
[Ni(HIM2-py)2NO3]+. 8 Fig. 1 illustrates the unit cells of the systems and the
corresponding complexes used in the cluster models.

Figure 1: Illustration of the structure of (left) the [Co(N2S2O4C8H10)2]2-
unit cell, (middle) the [Fe(tpa)Ph]- unit cell, and (right) the
[Ni(HIM2-py)2NO3]+ unit cell. The color code is: Co green, Fe orange, Ni
turquoise, Na brown, S yellow, C grey, O red, N blue, and H white.

#####  [Co(N2S2O4C8H10)2]2-.

One of the most promising high-anisotropy SIMs has been reported by Rechkemmer
et al. 6 It has a single cobalt ion center which is bound to two doubly
deprotonated 1,2-bis(methanesulfonamido)benzene ligands resulting in a
distorted tetrahedral coordination sphere. In the crystal, the charge is
compensated by two NHEt3+ cations. The zero-field splitting parameter
$D^{\mathrm{exp}}=-115\pm 20~{}\text{cm}^{-1}$ has been experimentally
determined by fitting to AC and DC susceptibility and magnetic hysteresis
measurements. 6

#####  [Fe(tpa)Ph]-.

This is one example out of a big family of similar complexes which all feature
an iron center ion in trigonal-pyramidal surrounding of a pyrrolide ligand
tpa${}^{\text{R}}$. 7 The complex studied here has phenyl attached to the main
pyrrolide ligand (R=Ph). The counter ion in the crystal is Na+(H3COC2H4OCH3)3.
Experimentally the magnetic parameters have been determined to be
$D^{\mathrm{exp}}=-26\pm 2~{}\text{cm}^{-1}$,
$E^{\mathrm{exp}}=5~{}\text{cm}^{-1}$ by AC and DC magnetic susceptibility
measurements and fitting procedures. 7 This family of complexes was also
theoretically studied by Atanasov et al. on the level of complete active space
self-consistent field (CASSCF) wave functions in conjunction with N-electron
valence perturbation theory (NEVPT2) and quasi-degenerate perturbation theory
(QDPT). 35

#####  [Ni(HIM2-py)2NO3]+.

One example of a nickel complex with a high magnetic anisotropy has been
synthesized by Rogez et al. 8. The nickel ion is influenced by a highly
distorted octahedral coordination sphere consisting of two bidentate ligands
HIM2-py as well as an O,O’-chelating nitrate ligand. In the crystal, the
positive charge is compensated by a second NO3- ion without direct contact to
the nickel ion. Measurements of the magnetization versus field, HF-HFEPR, and
FDMRS resulted in fitted magnetic parameters $D^{\mathrm{exp}}=-10.1\pm
0.1~{}\text{cm}^{-1}$ and $E^{\mathrm{exp}}=0.202\pm 0.01~{}\text{cm}^{-1}$.8

###  3.2 Data Generated via AIMD

For all training data sets an ab initio molecular dynamics (AIMD) calculation
has been carried out using the PAW method 49, 50 with the PBE functional, 51
as implemented in the VASP program package.52, 53, 54 A Hubbard $U$ correction
term has been used with the simplified (rotationally invariant) approach
introduced by Dudarev _et al._55 The used values of $U$ are 3.3 eV for Co, 4.0
eV for Fe and 6.4 eV for Ni.56 Dispersion corrections were applied by the zero
damping DFT-D3 method.57

As a starting point, the crystal structure was first optimized using an energy
cutoff for the plane-wave basis of 600 eV and the ’accurate precision’
settings in VASP, the projection operators were evaluated in real-space. The
Brillouin zone was sampled by a Monkhorst-Pack grid (Co system:
2$\times$2$\times$2; Fe system: 3$\times$3$\times$2; Ni system:
4$\times$4$\times$2). The stress tensor was calculated and all degrees of
freedom were allowed to change in relaxation.

After that, we performed an AIMD simulation with VASP, using the same settings
as before but with a cutoff energy of 400 eV and the ’normal precision’
settings. Only the $\Gamma$ point was considered. For each system shown in
Fig. 1 we calculated 5000 time steps of 1 fs for the temperatures 100K, 300K,
400K, 450K, and 500K each using a Nosé–Hoover thermostat as implemented in
VASP.

We chose more structures from the MD runs at the higher temperatures, i.e.
200, 600, 800, 900, and 1000 sample structures randomly chosen from the MD
runs at 100K, 300K, 400K, 450K, and 500K, respectively. This was done because
more diverse conformations are visited at higher temperature and therefore
dynamics performed at a higher temperature contains more information necessary
for the construction of reliable interatomic potentials.

###  3.3 D-Tensor Data

To generate reference values for the $\mathbf{D}$ tensor on which models were
subsequently trained, we cut the corresponding molecules from the periodic
structures of the AIMD simulation, as illustrated in Fig. 1. Effects of the
neighboring molecules and the crystal field on the $\mathbf{D}$ tensor have
thus been neglected for the present study. In total, we obtained 3500
configurations for each test system, for which magnetic properties were
calculated using the Molpro program package.58 The orbitals were optimized by
the configuration-averaged Hartree-Fock (CAHF) procedure59, 60, 61, 62 using
the Karlsruhe def2-SVP basis sets.63 The active space consisted of the five d
orbitals in each of the systems (see below). Based on these orbitals, a
complete-active space configuration interaction (CASCI) calculation was
performed to obtain the spin-free states. These were afterwards used in a
spin-orbit configuration interaction (SO-CI) calculation, using a mean-field
spin-orbit operator64, 65, 66, based on the CAHF average density for the mean
field.

The following setup is used for each of the test systems: For
[Co(N2S2O4C8H10)2]2- we chose a CAS(7,5) and the SO-CI calculation was carried
out using 40 doublet and 10 quartet CASCI states; for [Fe(tpa)Ph]- we chose a
CAS(6,5) and the SO-CI calculation was based on 50 singlet, 45 triplet, and 5
quintet CASCI states; for [Ni(HIM2-py)2NO3]+ we finally chose a CAS(8,5) and
the SO-CI calculation used 15 singlet and 10 triplet CASCI states.

From the SO-CI calculations, the $\mathbf{D}$ tensor was extracted using the
pseudo-spin procedure described by Chibotaru and Ungur. 67

##  4 Results and Discussion

Here we apply the proposed approach for machine learning symmetric traceless
tensors, namely $\mathbf{D}$ tensors, to systems described in Section 3. For
the example of [Co(N2S2O4C8H10)2]2- we demonstrate the applicability of our
approach to studying relevant magnetic properties of SMMs.

###  4.1  Machine Learning of the D Tensor

To study the performance of the proposed approach for machine learning
magnetic anisotropy tensors, or specifically $\mathbf{D}$ tensors, we trained
our model on [Co(N2S2O4C8H10)2]2-, [Fe(tpa)Ph]-, and [Ni(HIM2-py)2NO3]+ data.
Each data set contained structures and $\mathbf{D}$ tensor elements of 3500
configurations. In Fig. 2 we report the mean absolute error (MAE) and the
root-mean-square error (RMSE) in $\text{cm}^{-1}$. The cutoff radius employed
in the definition of invariant atomic representation was set to $8.0$ Å. All
other hyper parameters are described in Section 2.2.2.

Figure 2: Learning curves for the [Co(N2S2O4C8H10)2]2-, [Fe(tpa)Ph]-, and
[Ni(HIM2-py)2NO3]+ data sets. The root-mean-squared errors (RMSE) and the mean
absolute errors (MAE) of the $\mathbf{D}$-tensor components are plotted
against the training set size. Linear fits are displayed for clarity and
shaded areas denote the 95 % confidence intervals for linear regression.

Fig. 2 reports the MAE and RMSE of $\mathbf{D}$ tensor predictions as a
function of the number of training samples. For all training set sizes, we
have randomly drawn 300 additional structures as validation data to track
overfitting during training. Since the validation data indirectly influence
the selected set of trainable parameters, all values presented in Fig. 2 are
obtained for the test data that have not been seen during training. For
example, for 2900 training data we have used the remaining 300 structures to
test the model, while for 1500 training data 1700 structures were used for the
same purpose.

In general, we notice that the proposed approach leads to models that learn
quite efficiently on the reference data. For example, the RMSE for the
[Co(N2S2O4C8H10)2]2- data set is reduced by a factor of two when doubling the
training data set size. It can be observed that the RMSE for
[Ni(HIM2-py)2NO3]+ is slightly higher compared to other systems, while
[Fe(tpa)Ph]- shows lower RMSE values for smaller training data set sizes. The
former observation can be explained by the higher flexibility of
[Ni(HIM2-py)2NO3]+. This leads to a larger conformational space sampled during
ab initio molecular dynamics (AIMD) in Section 3.2 and as a result to a
broader range of $\mathbf{D}$-tensor elements, see Fig. S2. For [Fe(tpa)Ph]-
the situation is different since structural differences lead only to a slight
variation in $\mathbf{D}$-tensor elements, see Fig. S1, compared to other
systems, Fig. 3 and Fig. S2.

A direct comparison to previous models dealing with machine learning magnetic
anisotropy tensors, e.g., Refs. 33, 34, is impossible since no standard
benchmark is available in the literature, in contrast to, e.g., QM9 68, 69 and
MD17 70, 71, 72 data sets for testing machine-learned interatomic potentials.
Therefore, we use the presented data sets, available free-of-charge from Ref.
73, to benchmark the models developed for predicting $\mathbf{D}$ tensors.

In previous work,33 for [Co(N2S2O4C8H10)2]2- an RMSE value of 2.2
$\text{cm}^{-1}$ was obtained when training on 900 structures, while we obtain
1.5 $\text{cm}^{-1}$. In Ref. 33 the training data set was generated starting
from an optimized structure in a vacuum and then displacing atoms by a maximum
of $\pm 0.05$ Å (500 structures), $\pm 0.1$ Å (500 structures), and $\pm 0.2$
Å (500 structures), while we extracted the conformations from an AIMD
simulation. Thus, we expect that our data sets cover a broader conformational
space, making the training more difficult (which would lead to higher MAE/RMSE
values) but facilitating the prediction of properties for out-of-sample
configurations.

In Fig. 3, we show the correlation of the reference and machine-learned
results of each $D_{ij}$ element in the $\mathbf{D}$ tensor for the
[Co(N2S2O4C8H10)2]2- system. We again use only those structures the machine
learning model has not seen during training. From Fig. 3 we clearly see a
perfect correlation between both values in accordance with the low MAE (RMSE)
value of 0.30 $\text{cm}^{-1}$ (0.58 $\text{cm}^{-1}$), obtained by training
the model on 2900 reference structures. Similar results were obtained for
[Fe(tpa)Ph]- and [Ni(HIM2-py)2NO3]+, see Fig. S1 and Fig. S2. In total, taking
into account the broadly sampled conformational space and an excellent out-of-
sample predictive accuracy, our approach should be applicable to large-scale
molecular dynamics simulations.

Figure 3: Correlation of the machine-learned symmetric elements ($i\geq j$) of
the zero-field splitting tensor ($D_{ij}$) with the corresponding reference
values ($D_{ij}^{\mathrm{ref}}$) for all structures in the test
[Co(N2S2O4C8H10)2]2- data. For all predictions, the model trained on 2900
reference structures was used. The respective coordinate system of the
periodic box as well as an example substructure for which $\mathbf{D}$ tensor
was computed are shown as an inset.

###  4.2  Dynamics of the D Tensor

We assess the quality of machine-learned $\mathbf{D}$ tensors by applying the
machine learning methodology to the [Co(N2S2O4C8H10)2]2- system to study its
magnetic/physicochemical properties. While a detailed analysis is out of the
scope of this paper, we provide merely a qualitative overview to demonstrate
the broad applicability of our approach.

To analyze the properties of [Co(N2S2O4C8H10)2]2- we ran molecular dynamics
(MD) simulations in the canonical (NVT) statistical ensemble, carried out
within the ASE simulation package46 using a Langevin thermostat at the
temperatures of 25, 50, 75, and 100 K. All MD runs were performed over 2.5 ns
using a time step of 0.5 fs. The atomic velocities were initialized with a
Maxwell–Boltzmann distribution for the temperatures of 25, 50, 75, and 100 K,
respectively.

Forces for molecular dynamics were generated by an ensemble of three machine-
learned interatomic potentials, see Section 2.3. To train the machine-learned
interatomic potentials a cutoff radius of 6.5 Å was employed. The ensembling
technique provides us with an error estimate of the potential during
simulation. We have found the machine-learned potential to be very accurate
with the uncertainty between models ranging from $0.10\pm 0.06$ kcal/mol/Å to
$0.15\pm 0.10$ kcal/mol/Å, see Fig. S3. For $\mathbf{D}$ tensor predictions we
also use an ensemble of 3 models with an uncertainty ranging from $0.11\pm
0.06~{}\text{cm}^{-1}$ to $0.20\pm 0.11~{}\text{cm}^{-1}$, see Fig S3. These
values allow us to claim that our machine-learned models are suitable for the
following analysis.

####  4.2.1 Dependence of the D Tensor on the Structure

To study the structural dependence of the $\mathbf{D}$ tensor we evaluated the
average Co-N distance $d_{\ce{CoN}}$ and the tetrahedral order parameter
$q_{\widehat{\ce{NCoN}}}$ as 74

$\begin{split}&d_{\ce{CoN}}=\frac{1}{4}\sum_{i=1}^{4}d_{i}~{},\\\
&q_{\widehat{\ce{NCoN}}}=1-\frac{3}{8}\sum_{i=1}^{3}\sum_{j=i+1}^{4}\left(\cos\left(\theta_{ij}\right)+\frac{1}{3}\right)^{2}~{},\end{split}$
(14)

where $d_{i}$ and $\theta_{ij}$ are the distance and angle formed by the metal
center and its neighboring nitrogen atoms. Note that in the literature usually
the average angle between nitrogen atoms belonging to the same ligand is used
as order parameter, 33 or its deviation from a perfect tetrahedral angle
$T_{d}=109.5^{\circ}$ is discussed. 75 We have employed the tetrahedral order
parameter $q_{\widehat{\ce{NCoN}}}$ often used when studying the structure of
liquid water. 76 This parameter contains the information about angular
distribution but is rescaled in such a way that its value varies between 0 (if
the arrangement of all atoms is random) and 1 (in a perfect tetrahedral
network). Moreover, it has a value of 0.5 for a perfect quadratic planar
configuration, which allows us to look into both possible configurations of
the first coordination sphere of [Co(N2S2O4C8H10)2]2-, tetrahedral and
quadratic planar.

Fig. 4 shows the correlation of the magnetic axial anisotropy
$D=\tfrac{3}{2}Z$, where $Z$ is the eigenvalue of the $\mathbf{D}$ tensor with
the largest absolute value, with the order parameters presented above. Data
obtained from MD at 100 K are shown since they provide the broadest range of
conformations and the broadest range of $D$ values. From Fig. 4 we see that
the average Co-N distance $d_{\ce{CoN}}$ correlates only marginally with the
magnetic axial anisotropy $D$ for which we obtained a linear correlation
coefficient of $-0.22$. This is in agreement with results found in Ref. 33
taking into consideration that a minimal average Co-N distance of about 1.98 Å
could be sampled at 100 K. Note that in Ref. 33 a stronger correlation was
found for distances less than 2.0 Å. However, from our MD simulations we see
that such configurations are too rare at 100 K to be taken into account.

Figure 4: Correlation of the magnetic axial anisotropy $D=\tfrac{3}{2}Z$,
where $Z$ is the largest eigenvalue of the $\mathbf{D}$ tensor, with (left)
the average Co-N distance $d_{\ce{CoN}}$ and (right) tetrahedral order
parameter $q_{\widehat{\ce{NCoN}}}$. The latter yields a value of 1 for the
perfect tetrahedral structure, 0.5 for a perfect quadratic planar structure,
and 0 for a random mutual arrangement of central and neighboring atoms. The
tetrahedral order parameter $q_{\widehat{\ce{NCoN}}}$ yields a linear
correlation coefficient of 0.88, while the average Co-N distance does not
correlate with $D$ and yields a linear correlation coefficient of $-0.22$.

Fig. 4 clearly shows a strong dependence of the magnetic axial anisotropy $D$
on the tetrahedral order parameter $q_{\widehat{\ce{NCoN}}}$ for which we
obtained a linear correlation coefficient of 0.88. When changing the structure
from tetrahedral coordination to a quadratic planar one, a strong increase of
the magnetic anisotropy is observed. Our results are in perfect agreement with
recent results that suggest that the $\widehat{\ce{NCoN}}$ angle represents
the main path for improving Co2+ single-ion magnets. 33, 6, 77, 78, 79, 80

####  4.2.2 Thermal Distribution of the D Tensor

Besides the correlation of the magnetic axial anisotropy $D$ with structural
order parameters $d_{\ce{CoN}}$ and $d_{\widehat{\ce{NCoN}}}$ we study the
variation in magnitude and orientation due to thermal fluctuations in the
$\mathbf{D}$ anisotropy tensor. Fig. 5 (left) shows the distribution of the
individual $D_{ij}$ components of the $\mathbf{D}$ tensor sampled over 2.5 ns.
It can be seen that all elements $D_{ij}$ are approximately normally
distributed with broader distributions for higher temperatures. The elements
are symmetric with $D_{ij}=D_{ji}$, hence only $D_{ij}$ for $i\geq j$ are
shown. It should be noted that we have found that the mean of the distribution
of each element $D_{ij}$ remains almost unchanged and the maximal deviation
equals 0.1 $\text{cm}^{-1}$. However, the distribution of each element
$D_{ij}$ is broadened by approximately a factor of 2 when increasing the
temperature from 25 K to 100 K.

In order to relate the dynamics of the magnetic axial anisotropy
$D=\tfrac{3}{2}Z$ to the structural dynamics of the [Co(N2S2O4C8H10)2]2-
complex we display the respective distribution in Fig. 5 (right) along with
the distribution for $E=\left|X-Y\right|/2$. Note that $X$, $Y$, and $Z$ are
the eigenvalues of $\mathbf{D}$ with $Z$ being the one with the largest
absolute value. We have found the values for $D$ to range from $-105.7\pm
3.1~{}\text{cm}^{-1}$ to $-106.0\pm 6.2~{}\text{cm}^{-1}$, i.e. the mean is
almost temperature-independent while the standard deviation is doubled when
increasing the temperature by a factor of 4 in accordance with results for
$D_{ij}$. Note that our values estimated from MD simulations with machine
learning are very close to the experimental one of $-115\pm
20~{}\text{cm}^{-1}$. 6 The values for $E$ range from $1.2\pm
0.2~{}\text{cm}^{-1}$ to $1.3\pm 0.5~{}\text{cm}^{-1}$.

Figure 5: (Left) Thermal broadening of the elements $D_{ij}$ with $i\geq j$ of
the zero-field splitting tensor $\mathbf{D}$ obtained from sampling
configurations over 2.5 ns long molecular dynamics. All values were computed
employing the machine learning model described in Section 2.2. For all
predictions, the model trained on 2900 reference structures was used. (Right)
Temperature dependence of $D=\tfrac{3}{2}Z$ and $E=\left|X-Y\right|/2$, where
$X$, $Y$, and $Z$ are the eigenvalues of the $\mathbf{D}$ tensor, with $Z$
being its largest absolute eigenvalue.

####  4.2.3 Time Correlation Functions and Spin-Phonon Coupling

Finally, and mainly as an outlook to future applications of our approach, we
study the coupling of the $\mathbf{D}$ tensor of the [Co(N2S2O4C8H10)2]2-
complex to the dynamics of the periodic atomic structure (velocity vectors
$\mathbf{v}=\dot{\mathbf{x}}$). For this purpose, we define time-dependent
$\mathbf{D}$ tensor and velocity autocorrelation functions (ACFs) as

$\begin{split}&C_{\mathbf{D}\mathbf{D}}\left(t\right)=\frac{1}{6}\sum_{i=1}^{3}\sum_{j\geq
i}^{3}\frac{\langle D_{ij}\left(0\right)D_{ij}\left(t\right)\rangle}{\langle
D_{ij}\left(0\right)D_{ij}\left(0\right)\rangle}~{},\\\
&C_{\mathbf{v}\mathbf{v}}\left(t\right)=\frac{1}{3N_{\mathrm{at}}}\sum_{i=1}^{N_{\mathrm{at}}}\sum_{j=1}^{3}\frac{\langle
v_{ij}\left(0\right)v_{ij}\left(t\right)\rangle}{\langle
v_{ij}\left(0\right)v_{ij}\left(0\right)\rangle}~{},\end{split}$ (15)

respectively. Using these ACFs it is possible to compute the corresponding
spectra by performing a Fourier transform or by employing the maximum-entropy
approach. 81, 82, 83 In this work we compute the power spectral density
function $S_{\omega}$ of interest employing the memspectrum package. 84

Fig. 6 shows the velocity and $\mathbf{D}$ tensor ACFs (left) as well as the
respective power spectral density functions
$S_{\omega}[C_{\mathbf{v}\mathbf{v}}\left(t\right)]$ and
$S_{\omega}[C_{\mathbf{D}\mathbf{D}}\left(t\right)]$ (right). From the latter,
in principle, the spin-phonon coupling coefficients $\partial D/\partial
q_{\alpha}$ beyond the harmonic approximation can be calculated, which provide
the interaction strength between the spin and the atomic movements.85 The
present study is restricted to the demonstration of the capability of the
proposed approach, a detailed analysis will be presented in future work. For
now, we computed the spectra for the $\Gamma$-point only, while the inclusion
of the full Brillouin zone ($k$-dependence) may be important to estimate the
spin lifetime. Even with our preliminary spectra, we are able to deduce
similar conclusions as previous work85 found for [Fe(tpa)Ph]- in that the most
important vibrations in the spin-phonon relaxation process are the low-energy
ones. They are predominantly populated under typical experimental conditions.

Figure 6: (Left) Velocity and $\mathbf{D}$ tensor autocorrelation functions,
$C_{\mathbf{v}\mathbf{v}}\left(t\right)$ and
$C_{\mathbf{D}\mathbf{D}}\left(t\right)$ respectively, obtained by averaging
over 5 independent MD trajectories at 50 K of a length of 500 ps. (Right)
Power spectral density function $S_{\omega}$ of the velocity and $\mathbf{D}$
tensor autocorrelation functions. The gray vertical lines at the bottom of the
graph correspond to harmonic frequencies.

##  5  Conclusion

In this work, we have presented a machine learning approach based on Gaussian
moments 24, 30 for tensorial properties on the example of the zero-field
splitting ($\mathbf{D}$) tensor. It enables an efficient prediction and
modeling of magnetic properties of single-molecule magnets. The presented
approach was extensively tested on three systems [Co(N2S2O4C8H10)2]2-,
[Fe(tpa)Ph]-, and [Ni(HIM2-py)2NO3]+ in terms of its predictive accuracy and
its reliability during real-time simulations.

Training the proposed model on the respective reference $\mathbf{D}$ values
for the [Co(N2S2O4C8H10)2]2- we have observed an improved accuracy compared to
previous studies, 33 especially taking into account the larger sampled
configurational space. In total, for all systems tested in this work we could
achieve an MAE of 0.3–0.4 $\text{cm}^{-1}$ and an RMSE of 0.6–0.7
$\text{cm}^{-1}$ for the models trained on 2900 reference structures.
Moreover, we have shown that the proposed approach, once trained on a
sufficiently large configurational space, is able to predict $\mathbf{D}$
tensor values for millions of conformations not seen before with a negligibly
small uncertainty of 0.11–0.20 $\text{cm}^{-1}$ obtained by employing the
query-by-committee approach. This demonstrates the excellent generalization
capability of our approach.

In combination with machine-learned interatomic potentials we were able to run
2.5 ns long molecular dynamics simulations at temperatures of 25, 50, 75, and
100 K. Using the respective trajectories we could analyze several properties
of the [Co(N2S2O4C8H10)2]2- complex taken as an example. Analysing the
dependence of the magnetic axial anisotropy on average Co-N distance and an
$\widehat{\ce{NCoN}}$ angle-dependent order parameter
$q_{\widehat{\ce{NCoN}}}$ we have found that the $\widehat{\ce{NCoN}}$ angle
represents the main path for improving Co2+ single-ion magnets, in accordance
with recent results. 33, 6, 77, 78, 79, 80 Moreover, we could estimate the
thermal average of $D$ and $E$ values. For the former, a very good agreement
with the experiment has been found, while for the latter no experimental data
are available.

Besides the structure, we investigated the dynamical behavior of the
$\mathbf{D}$ tensor of the [Co(N2S2O4C8H10)2]2- complex. Even in such
preliminary work, we observed the expected behavior that the low-energy
vibrations are important for the spin-phonon relaxation process.

In summary, our developments aim to provide an alternative way for the
efficient modeling of magnetic properties of single molecular magnets via
machine learning. While the current setup is based on a relatively simple
complete-active space configuration interaction treatment of the molecular
properties, our approach can be easily extended to more elaborate methods,
e.g. by using transfer learning. Future work will furthermore deal with an
application of the proposed methodology to allow for a detailed analysis of
spin-phonon relaxation processes.

We thank the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)
for supporting this work by funding EXC 2075 - 390740016 under Germany’s
Excellence Strategy and through grant no INST 40/575-1 FUGG (JUSTUS 2
cluster). We acknowledge the support by the Stuttgart Center for Simulation
Science (SimTech). We also like to acknowledge the support by the state of
Baden-Württemberg through the bwHPC consortium for providing computer time. V.
Z. acknowledges financial support received in the form of a PhD scholarship
from the Studienstiftung des Deutschen Volkes (German National Academic
Foundation).

## References

  * Gatteschi et al. 2006 Gatteschi, D.; Sessoli, R.; Villain, J. _Molecular Nanomagnets_ ; Oxford University Press, 2006
  * Craig and Murrie 2015 Craig, G. A.; Murrie, M. 3d single-ion magnets. _Chem. Soc. Rev._ 2015, _44_ , 2135–2147
  * Gaita-Ariño et al. 2019 Gaita-Ariño, A.; Luis, F.; Hill, S.; Coronado, E. Molecular spins for quantum computation. _Nat. Chem._ 2019, _11_ , 301–309
  * Sanvito 2011 Sanvito, S. Molecular spintronics. _Chem. Soc. Rev._ 2011, _40_ , 3336–3355
  * Gatteschi and Sessoli 2003 Gatteschi, D.; Sessoli, R. Quantum Tunneling of Magnetization and Related Phenomena in Molecular Materials. _Angew. Chem. Int. Ed._ 2003, _42_ , 268–297
  * Rechkemmer et al. 2016 Rechkemmer, Y.; Breitgoff, F. D.; van der Meer, M.; Atanasov, M.; Hakl, M.; Orlita, M.; Neugebauer, P.; Neese, F.; Sarkar, B.; van Slageren, J. A four-coordinate cobalt(II) single-ion magnet with coercivity and a very high energy barrier. _Nat. Commun._ 2016, _7_ , 10467
  * Harman et al. 2010 Harman, W. H.; Harris, T. D.; Freedman, D. E.; Fong, H.; Chang, A.; Rinehart, J. D.; Ozarowski, A.; Sougrati, M. T.; Grandjean, F.; Long, G. J.; Long, J. R.; Chang, C. J. Slow Magnetic Relaxation in a Family of Trigonal Pyramidal Iron(II) Pyrrolide Complexes. _J. Am. Chem. Soc._ 2010, _132_ , 18115–18126
  * Rogez et al. 2005 Rogez, G.; Rebilly, J.-N.; Barra, A.-L.; Sorace, L.; Blondin, G.; Kirchner, N.; Duran, M.; van Slageren, J.; Parsons, S.; Ricard, L.; Marvilliers, A.; Mallah, T. Very Large Ising-Type Magnetic Anisotropy in a Mononuclear Ni(II) Complex. _Angew. Chem. Int. Ed._ 2005, _44_ , 1876 –1879
  * Bamberger et al. 2021 Bamberger, H.; Albold, U.; Midlíková, J. D.; Su, C.-Y.; Deibel, N.; Hunger, D.; Hallmen, P. P.; Neugebauer, P.; Beerhues, J.; Demeshko, S.; Meyer, F.; Sarkar, B.; van Slageren, J. Iron(II), Cobalt(II), and Nickel(II) Complexes of Bis(sulfonamido)benzenes: Redox Properties, Large Zero-Field Splittings, and Single-Ion Magnets. _Inorg. Chem._ 2021, _60_ , 2953–2963
  * Blank et al. 1995 Blank, T. B.; Brown, S. D.; Calhoun, A. W.; Doren, D. J. Neural network models of potential energy surfaces. _J. Chem. Phys._ 1995, _103_ , 4129–4137
  * Lorenz et al. 2004 Lorenz, S.; Groß, A.; Scheffler, M. Representing high-dimensional potential-energy surfaces for reactions at surfaces by neural networks. _Chem. Phys. Lett._ 2004, _395_ , 210–215
  * Lorenz et al. 2006 Lorenz, S.; Scheffler, M.; Gross, A. Descriptions of surface chemical reactions using a neural network representation of the potential-energy surface. _Phys. Rev. B_ 2006, _73_ , 115431
  * Behler and Parrinello 2007 Behler, J.; Parrinello, M. Generalized Neural-Network Representation of High-Dimensional Potential-Energy Surfaces. _Phys. Rev. Lett._ 2007, _98_ , 146401
  * Behler 2011 Behler, J. Neural network potential-energy surfaces in chemistry: a tool for large-scale simulations. _Phys. Chem. Chem. Phys._ 2011, _13_ , 17930–17955
  * Behler 2016 Behler, J. Perspective: Machine learning potentials for atomistic simulations. _J. Chem. Phys._ 2016, _145_ , 170901
  * Shapeev 2016 Shapeev, A. V. Moment Tensor Potentials: A Class of Systematically Improvable Interatomic Potentials. _Multiscale Model. Simul._ 2016, _14_ , 1153–1173
  * Gastegger et al. 2017 Gastegger, M.; Behler, J.; Marquetand, P. Machine learning molecular dynamics for the simulation of infrared spectra. _Chem. Sci._ 2017, _8_ , 6924–6935
  * Gubaev et al. 2018 Gubaev, K.; Podryabinkin, E. V.; Shapeev, A. V. Machine learning of molecular properties: Locality and active learning. _J. Chem. Phys._ 2018, _148_ , 241727
  * Dral et al. 2018 Dral, P. O.; Barbatti, M.; Thiel, W. Nonadiabatic Excited-State Dynamics with Machine Learning. _J. Phys. Chem. Lett._ 2018, _9_ , 5660–5663
  * Zhang et al. 2018 Zhang, L.; Han, J.; Wang, H.; Car, R.; E, W. Deep Potential Molecular Dynamics: A Scalable Model with the Accuracy of Quantum Mechanics. _Phys. Rev. Lett._ 2018, _120_ , 143001
  * Yao et al. 2018 Yao, K.; Herr, J. E.; Toth, D. W.; Mckintyre, R.; Parkhill, J. The TensorMol-0.1 model chemistry: a neural network augmented with long-range physics. _Chem. Sci._ 2018, _9_ , 2261–2269
  * Westermayr et al. 2019 Westermayr, J.; Gastegger, M.; Menger, M. F. S. J.; Mai, S.; González, L.; Marquetand, P. Machine learning enables long time scale molecular photodynamics simulations. _Chem. Sci._ 2019, _10_ , 8100–8107
  * Unke and Meuwly 2019 Unke, O. T.; Meuwly, M. PhysNet: A Neural Network for Predicting Energies, Forces, Dipole Moments, and Partial Charges. _J. Chem. Theory Comput._ 2019, _15_ , 3678–3693
  * Zaverkin and Kästner 2020 Zaverkin, V.; Kästner, J. Gaussian Moments as Physically Inspired Molecular Descriptors for Accurate and Scalable Machine Learning Potentials. _J. Chem. Theory Comput._ 2020, _16_ , 5410–5421
  * Molpeceres et al. 2020 Molpeceres, G.; Zaverkin, V.; Kästner, J. Neural-network assisted study of nitrogen atom dynamics on amorphous solid water – I. adsorption and desorption. _Mon. Not. R. Astron. Soc._ 2020, _499_ , 1373–1384
  * Dral 2020 Dral, P. O. Quantum Chemistry in the Age of Machine Learning. _J. Phys. Chem. Lett._ 2020, _11_ , 2336–2347
  * Molpeceres et al. 2021 Molpeceres, G.; Zaverkin, V.; Watanabe, N.; Kästner, J. Binding energies and sticking coefficients of H2 on crystalline and amorphous CO ice. _A &A_ 2021, _648_ , A84
  * Christensen et al. 2020 Christensen, A. S.; Bratholm, L. A.; Faber, F. A.; von Lilienfeld, O. FCHL revisited: Faster and more accurate quantum machine learning. _J. Chem. Phys._ 2020, _152_ , 44107
  * Faber et al. 2018 Faber, F. A.; Christensen, A. S.; Huang, B.; von Lilienfeld, O. A. Alchemical and structural distribution based representation for universal quantum machine learning. _J. Chem. Phys._ 2018, _148_ , 241717
  * Zaverkin et al. 2021 Zaverkin, V.; Holzmüller, D.; Steinwart, I.; Kästner, J. Fast and Sample-Efficient Interatomic Neural Network Potentials for Molecules and Materials Based on Gaussian Moments. _J. Chem. Theory Comput._ 2021, _17_ , 6658–6670
  * Hornik 1991 Hornik, K. Approximation capabilities of multilayer feedforward networks. _Neural Netw._ 1991, _4_ , 251 – 257
  * Novikov et al. 2021 Novikov, I.; Grabowski, B.; Kormann, F.; Shapeev, A. Machine-learning interatomic potentials reproduce vibrational and magnetic degrees of freedom. _ArXiv_ 2021, _abs/2012.12763_
  * Lunghi and Sanvito 2020 Lunghi, A.; Sanvito, S. Surfing Multiple Conformation-Property Landscapes via Machine Learning: Designing Single-Ion Magnetic Anisotropy. _J. Phys. Chem. C_ 2020, _124_ , 5802–5806
  * Lunghi 2020 Lunghi, A. Insights into the Spin-Lattice Dynamics of Organic Radicals Beyond Molecular Tumbling: A Combined Molecular Dynamics and Machine-Learning Approach. _Appl. Magn. Reson._ 2020, _51_ , 1343–1356
  * Atanasov et al. 2011 Atanasov, M.; Ganyushin, D.; Pantazis, D. A.; Sivalingam, K.; Neese, F. Detailed Ab Initio First-Principles Study of the Magnetic Anisotropy in a Family of Trigonal Pyramidal Iron(II) Pyrrolide Complexes. _Inorg. Chem._ 2011, _50_ , 7460–7477
  * Abragam and Bleaney 1970 Abragam, A.; Bleaney, B. _Electron Paramagnetic Resonance of Transition Ions_ ; Clarendon, 1970
  * Jacot et al. 2018 Jacot, A.; Gabriel, F.; Hongler, C. Neural Tangent Kernel: Convergence and Generalization in Neural Networks. NeurIPS. 2018
  * Hendrycks and Gimpel 2016 Hendrycks, D.; Gimpel, K. Gaussian Error Linear Units (GELUs). _arXiv_ 2016, _abs/1606.08415_
  * Elfwing et al. 2018 Elfwing, S.; Uchibe, E.; Doya, K. Sigmoid-Weighted Linear Units for Neural Network Function Approximation in Reinforcement Learning. _Neural Netw._ 2018, _107_ , 3–11
  * Ramachandran et al. 2018 Ramachandran, P.; Zoph, B.; Le, Q. V. Searching for Activation Functions. _ArXiv_ 2018, _abs/1710.05941_
  * Klambauer et al. 2017 Klambauer, G.; Unterthiner, T.; Mayr, A.; Hochreiter, S. Self-Normalizing Neural Networks. NeurIPS. 2017
  * Arora et al. 2019 Arora, S.; Du, S.; Hu, W.; Li, Z.; Salakhutdinov, R.; Wang, R. On Exact Computation with an Infinitely Wide Neural Net. NeurIPS. 2019
  * Lu et al. 2020 Lu, Y.; Gould, S.; Ajanthan, T. Bidirectional Self-Normalizing Neural Networks. _ArXiv_ 2020, _abs/2006.12169_
  * Kingma and Ba 2015 Kingma, D. P.; Ba, J. Adam: A Method for Stochastic Optimization. _CoRR_ 2015, _abs/1412.6980_
  * Prechelt 2012 Prechelt, L. In _Neural Networks: Tricks of the Trade: Second Edition_ ; Montavon, G., Orr, G. B., Müller, K.-R., Eds.; Springer Berlin Heidelberg: Berlin, Heidelberg, 2012; pp 53–67
  * Larsen et al. 2017 Larsen, A. H.; Mortensen, J. J.; Blomqvist, J.; Castelli, I. E.; Christensen, R.; Dułak, M.; Friis, J.; Groves, M. N.; Hammer, B.; Hargus, C.; Hermes, E. D.; Jennings, P. C.; Jensen, P. B.; Kermode, J.; Kitchin, J. R.; Kolsbjerg, E. L.; Kubal, J.; Kaasbjerg, K.; Lysgaard, S.; Maronsson, J. B.; Maxson, T.; Olsen, T.; Pastewka, L.; Peterson, A.; Rostgaard, C.; Schiøtz, J.; Schütt, O.; Strange, M.; Thygesen, K. S.; Vegge, T.; Vilhelmsen, L.; Walter, M.; Zeng, Z.; Jacobsen, K. W. The atomic simulation environment—a Python library for working with atoms. _J. Phys. Condens. Matter_ 2017, _29_ , 273002
  * Settles 2009 Settles, B. _Active Learning Literature Survey_ ; Computer Sciences Technical Report 1648, 2009
  * Abadi et al. 2015 Abadi, M.; Agarwal, A.; Barham, P.; Brevdo, E.; Chen, Z.; Citro, C.; Corrado, G. S.; Davis, A.; Dean, J.; Devin, M.; Ghemawat, S.; Goodfellow, I.; Harp, A.; Irving, G.; Isard, M.; Jia, Y.; Jozefowicz, R.; Kaiser, L.; Kudlur, M.; Levenberg, J.; Mané, D.; Monga, R.; Moore, S.; Murray, D.; Olah, C.; Schuster, M.; Shlens, J.; Steiner, B.; Sutskever, I.; Talwar, K.; Tucker, P.; Vanhoucke, V.; Vasudevan, V.; Viégas, F.; Vinyals, O.; Warden, P.; Wattenberg, M.; Wicke, M.; Yu, Y.; Zheng, X. TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems. 2015; https://www.tensorflow.org/, Software available from tensorflow.org
  * Blöchl 1994 Blöchl, P. E. Projector augmented-wave method. _Phys. Rev. B_ 1994, _50_ , 17953
  * Kresse and Joubert 1999 Kresse, G.; Joubert, D. From ultrasoft pseudopotentials to the projector augmented-wave method. _Phys. Rev. B_ 1999, _59_ , 1758–1775
  * Perdew et al. 1997 Perdew, J.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. _Phys. Rev. Lett._ 1997, _78_ , 1396
  * Kresse and Hafner 1993 Kresse, G.; Hafner, J. Ab initio molecular dynamics for liquid metals. _Phys. Rev. B_ 1993, _47_ , 558
  * Kresse and Furthmüller 1996 Kresse, G.; Furthmüller, J. Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. _Comput. Mater. Sci._ 1996, _6_ , 15 – 50
  * Kresse and Furthmüller 1996 Kresse, G.; Furthmüller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. _Phys. Rev. B_ 1996, _54_ , 11169
  * Dudarev et al. 1998 Dudarev, S. L.; Botton, G. A.; Savrasov, S. Y.; Humphreys, C. J.; Sutton, A. P. Electron-energy-loss spectra and the structural stability of nickel oxide: An LSDA+U study. _Phys. Rev. B_ 1998, _57_ , 1505
  * Mann et al. 2016 Mann, G. W.; Lee, K.; Cococcioni, M.; Smit, B.; Neaton, J. B. First-principles Hubbard U approach for small molecule binding in metal-organic frameworks. _J. Chem. Phys._ 2016, _144_ , 174104
  * Grimme et al. 2010 Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. _J. Chem. Phys._ 2010, _132_ , 154104
  * Werner et al. 2020 Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M.; others MOLPRO, version 2020.0, a package of ab initio programs. 2020; see https://www.molpro.net.
  * McWeeny 1996 McWeeny, R. _Methods of Molecular Quantum Mechanics_ ; Academic Press, 1996
  * McWeeny 1974 McWeeny, R. SCF theory for excited states. _Mol. Phys._ 1974, _28_ , 1273–1282
  * Hallmen et al. 2017 Hallmen, P. P.; Köppl, C.; Rauhut, G.; Stoll, H.; Van Slageren, J. Fast and reliable ab initio calculation of crystal field splittings in lanthanide complexes. _J. Chem. Phys_ 2017, _147_ , 164101
  * Calvello et al. 2018 Calvello, S.; Piccardo, M.; Rao, S. V.; Soncini, A. CERES: An ab initio code dedicated to the calculation of the electronic structure and magnetic properties of lanthanide complexes . _J. Comput. Chem._ 2018, _39_ , 328–337
  * Weigend and Ahlrichs 2005 Weigend, F.; Ahlrichs, R. Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy. _Phys. Chem. Chem. Phys._ 2005, _7_ , 3297–3305
  * Marian and Wahlgren 1996 Marian, C. M.; Wahlgren, U. A new mean-field and ECP-based spin-orbit method. Applications to Pt and PtH. _Chem. Phys. Lett._ 1996, _251_ , 357–364
  * Hess et al. 1996 Hess, B. A.; Marian, C. M.; Wahlgren, U.; Gropen, O. A mean-field spin-orbit method applicable to correlated wavefunctions. _Chem. Phys. Lett._ 1996, _251_ , 365–371
  * Berning et al. 2000 Berning, A.; Schweizer, M.; Werner, H.-J.; Knowles, P. J.; Palmieri, P. Spin-orbit matrix elements for internally contracted multireference configuration interaction wavefunctions. _Mol. Phys._ 2000, _98_ , 1823–1833
  * Chibotaru and Ungur 2012 Chibotaru, L. F.; Ungur, L. Ab initio calculation of anisotropic magnetic properties of complexes. I. Unique definition of pseudospin Hamiltonians and their derivation. _J. Chem. Phys._ 2012, _137_ , 064112
  * Ruddigkeit et al. 2012 Ruddigkeit, L.; van Deursen, R.; Blum, L. C.; Reymond, J.-L. Enumeration of 166 Billion Organic Small Molecules in the Chemical Universe Database GDB-17. _J. Chem. Inf. Model._ 2012, _52_ , 2864–2875
  * Ramakrishnan et al. 2014 Ramakrishnan, R.; Dral, P. O.; Rupp, M.; von Lilienfeld, O. A. Quantum chemistry structures and properties of 134 kilo molecules. _Sci. Data_ 2014, _1_ , 140022
  * Schütt et al. 2017 Schütt, K. T.; Arbabzadah, F.; Chmiela, S.; Müller, K. R.; Tkatchenko, A. Quantum-chemical insights from deep tensor neural networks. _Nat. Commun._ 2017, _8_ , 13890
  * Chmiela et al. 2017 Chmiela, S.; Tkatchenko, A.; Sauceda, H. E.; Poltavsky, I.; Schütt, K. T.; Müller, K.-R. Machine learning of accurate energy-conserving molecular force fields. _Sci. Adv._ 2017, _3_
  * Chmiela et al. 2018 Chmiela, S.; Sauceda, H. E.; Müller, K.-R.; Tkatchenko, A. Towards exact molecular dynamics simulations with machine-learned force fields. _Nat. Commun._ 2018, _9_ , 3887
  * Zaverkin et al. 2021 Zaverkin, V.; Netz, J.; Zills, F.; Köhn, A.; Kästner, J. Thermally Averaged Magnetic Anisotropy Tensors via Machine Learning Based on Gaussian Moments (Version v1) [Data set]. 2021; https://doi.org/10.5281/zenodo.5172156
  * Errington and Debenedetti 2001 Errington, J.; Debenedetti, P. Relationship between structural order and the anomalies of liquid water. _Nature_ 2001, _409_ , 318–321
  * Titiš et al. 2013 Titiš, J.; Miklovič, J.; Boča, R. Magnetostructural study of tetracoordinate cobalt(II) complexes. _Inorg. Chem. Comm._ 2013, _35_ , 72–75
  * Liu et al. 2018 Liu, J.; He, X.; Zhang, J. Z. H.; Qi, L.-W. Hydrogen-bond structure dynamics in bulk water: insights from ab initio simulations with coupled cluster theory. _Chem. Sci._ 2018, _9_ , 2065–2073
  * Fataftah et al. 2014 Fataftah, M. S.; Zadrozny, J. M.; Rogers, D. M.; Freedman, D. E. A Mononuclear Transition Metal Single-Molecule Magnet in a Nuclear Spin-Free Ligand Environment. _Inorg. Chem._ 2014, _53_ , 10716–10721
  * Carl et al. 2015 Carl, E.; Demeshko, S.; Meyer, F.; Stalke, D. Triimidosulfonates as Acute Bite-Angle Chelates: Slow Relaxation of the Magnetization in Zero Field and Hysteresis Loop of a CoII Complex. _Chem. Eur. J._ 2015, _21_ , 10109–10115
  * Suturina et al. 2017 Suturina, E. A.; Nehrkorn, J.; Zadrozny, J. M.; Liu, J.; Atanasov, M.; Weyhermüller, T.; Maganas, D.; Hill, S.; Schnegg, A.; Bill, E.; Long, J. R.; Neese, F. Magneto-Structural Correlations in Pseudotetrahedral Forms of the [Co(SPh)4]2– Complex Probed by Magnetometry, MCD Spectroscopy, Advanced EPR Techniques, and ab Initio Electronic Structure Calculations. _Inorg. Chem._ 2017, _56_ , 3102–3118
  * Wu et al. 2019 Wu, T.; Zhai, Y.-Q.; Deng, Y.-F.; Chen, W.-P.; Zhang, T.; Zheng, Y.-Z. Correlating magnetic anisotropy with the subtle coordination geometry variation of a series of cobalt(ii)-sulfonamide complexes. _Dalton Trans._ 2019, _48_ , 15419–15426
  * Burg 1975 Burg, J. P. _Maximum Entropy Spectral Analysis_ ; Stanford University, 1975
  * Vos 2013 Vos, K. A Fast Implementation of Burg’s Method. 2013, Available: www.opuscodec.org/docs/vos_fastburg.pdf
  * Martini et al. 2021 Martini, A.; Schmidt, S.; del Pozzo, W. Maximum Entropy Spectral Analysis: a case study. _ArXiv_ 2021, _abs/2106.09499_
  * Martini et al. 2021 Martini, A.; Schmidt, S.; del Pozzo, W. Maximum Entropy Spectrum. https://github.com/martini-alessandro/Maximum-Entropy-Spectrum, 2021
  * Lunghi et al. 2017 Lunghi, A.; Totti, F.; Sanvito, S.; Sessoli, R. Intra-molecular origin of the spin-phonon coupling in slow-relaxing molecular magnets. _Chem. Sci._ 2017, _8_ , 6051–6059