Benchmark for Ab Initio Prediction of Magnetic Structures based on Cluster Multipole Theory

The cluster multipole (CMP) expansion for magnetic structures provides a scheme to systematically generate candidate magnetic structures specifically including noncollinear magnetic configurations adapted to the crystal symmetry of a given material. A comparison with the experimental data collected on MAGNDATA shows that the most stable magnetic configurations in nature are linear combinations of only few CMPs. Furthermore, a high-throughput generalized local spin-density approximation (LSDA) calculation, in which each candidate magnetic structure was considered as an initial guess, was performed. We benchmark the predictive power of CMP+LSDA with $2\,935$ calculations, which show that (i) the CMP expansion administers an exhaustive list of candidate magnetic structures, (ii) CMP+LSDA can narrow down the possible magnetic configurations to a handful of computed configurations, and (iii) LSDA reproduces the experimental magnetic configurations with an accuracy of $\pm0.5\,\mu_{\text{B}}$.


I. INTRODUCTION
The grand challenge in first-principles calculation for magnetic materials is whether we can predict the experimental magnetic structure for a given material. Among a variety of possible functional materials, antiferromagnetic (AFM) compounds are a fascinating playground for materials design as they facilitate a wide range of fundamental phenomena and possible applications.
For example, in the context of AFM spintronics [1] there is a particular interest in noncollinear antiferromagnetism sparked by (i) its robustness against perturbations due to magnetic fields, (ii) a quasi-absence of magnetic stray fields disturbing for instance nearby electronic devices, and (iii) ultrafast dynamics of AFM domainwalls [2], as well as (iv) its ability to generate large magnetotransport effects [3]. Hence, the optimization of AFM materials would open the door for applications such as seamless and low-maintenance energy generation, ultrafast spintronics and robust data retention, as well as be a guide towards advancing fundamental understanding of magnetotransport.
However, first-principles calculations with the local spin density approximation (LSDA) in the framework of density functional theory (DFT) for magnetic materials have a problem: It is still an open question how accurately LSDA can reproduce the experimental magnetic ground state. While LSDA has been widely used in studies on various magnets [4], there has been no systematic benchmark calculation for noncollinear AFM materials. Previous attempts have been restricted to collinear magnetism [5] or even stricter symmetry constrains [6][7][8]. In regard to noncollinear AFM materials, high-throughput calculations have been limited to setting the experimentally determined magnetic configuration as an initial guess [9]. The lack of a systematic benchmark calculation is a consequence of the fact that it is a highly non-trivial task to investigate all the local minima in the LSDA energy landscape. Indeed, to search for all the (meta-)stable states, we need an exhaustive list of physically reasonable magnetic configurations for which LSDA calculations can be performed.
To this end, we devise the so-called cluster multipole (CMP) expansion, which enables the expansion of an arbitrary magnetic configuration in terms of an orthogonal basis set of magnetic multipole configurations. By means of the CMP expansion, a list of initial magnetic structures for self-consistent LSDA calculations is efficiently and systematically generated. With this at hand, a systematic high-throughput calculation with 2935 calculations has been performed.
The structure of the paper is as follows: In Section II, we explain the basic idea of CMP (Section II A) and setup for LSDA calculation (Section II B). In Section III, the magnetic configuration of 131 materials is predicted using a combination of the CMP expansion and LSDA (CMP+LSDA). A comparison to the experimental data shows that the magnetic ground state can be narrowed down to be among a handful of computed configurations and LSDA reproduces the experimental on-site magnetic moment with an accuracy of approximately ±0.5 µ B . This benchmark, which is summarized in Section IV, thus provides a solid foundation for the ab initio predictions of various magnetic properties.

II. METHODS
In this section we shortly discuss the employed methods, namely the CMP expansion and generalized LSDA. As the CMP expansion is a rather novel approach [10][11][12], it shall be motivated and set out in some detail. However for more background and details of the algorithm we refer the reader to Ref. [13]. LSDA on the other hand is a well arXiv:2008.13669v1 [cond-mat.mtrl-sci] 31 Aug 2020 established method [14]. It is available as part of many ab initio packages [15][16][17][18][19][20] in its generalized version [21,22], which is applicable to noncollinear AFM configurations. Here, we chose to use VASP [15,23,24] and hence we merely elaborate on the setup details employed in this study.

A. Cluster Multipole expansion
The cluster multipole (CMP) expansion for magnetic structures [10,13] provides an orthogonal basis set of magnetic configurations, which are symmetrized based on the crystallographic point group. In order to motivate the expansion, let us consider the vector Poisson equation: where j(r) = c∇×M (r) is the current density and M (r) is the magnetization density. Here, the Coulomb gauge ∇ · A(r) = 0 is invoked and the potential outside of the magnetization density is considered. The rotational invariance of ∇ 2 allows the vector gauge potential A(r) to be written in terms of spherical harmonics Y l m (θ, φ) as Here, l = −i (r × ∇) is the dimensionless orbital angular momentum with quantum number l ≥ 1, and magnetic quantum number −l ≤ m ≤ l [25,26]. Following M.-T. Suzuki et al. in Ref. [13] the CMPs for a magnetic configuration on a point form |m = (m 1 , m 2 , ..., m N ) T read m i is a local magnetic moment on the magnetic site i at position r i . For a given point group the point form is a set of all symmetrically equivalent points and can be classified into Wyckoff positions [27] in analogy to the Wyckoff positions of space groups. Here, N is the multiplicity of the Wyckoff position of the point form, that constitutes the magnetic configuration. As introduced by Ref. [13] a point form carrying a magnetic configuration is referred to as (magnetic) cluster in the context of the CMP expansion for magnetic structures. In contrast to Ref. [13], here we do not introduce toroidal moments. Symmetrization according to irreducible representations of the crystallographic point group allows for a physically meaningful expansion w.r.t. symmetrized harmonics where γ indicates the irreducible representation including the existing components of it. Here, the tabulated coefficients [28] c γ lm are chosen to be real valued. By mapping lγ → n through a Gram-Schmidt orthonormalization scheme the CMP basis is computed. This procedure is visualized in the Supplemental Material [29].
The CMP basis can be written as where e (n) i is a unit vector of a local magnetic moment on the magnetic site i. By convention n = 1, 2, 3 corresponds to ferromagnetism, while n ≥ 4 corresponds to more complicated higher order magnetic configurations including noncollinear magnetism. The definition of |n coincides with e µ lγ in Ref. [13] up to the choice of nor- In case that the period of the magnetic order coincides with that of the crystal structure, the propagation vector of the magnetic order q is zero. The magnetic structure is said to exhibit q = 0 magnetism. Note that 3 continuous degrees of freedom of rotation of the magnetic moment per magnetic site for a total of N magnetic sites yields 3N linearly independent magnetic configurations and thus n = 1, ..., 3N . In this work, the configuration space of q = 0 magnetic structures is explored.
The CMP basis defined in Equation (5) is complete and obeys the orthogonality relation Finally, the symmetry-adapted CMP coefficient reads In case of more than one inequivalent site exhibiting a magnetic moment, the space of all possible magnetic configurations is spanned by where d is the number of clusters. Based on the above, an arbitrary magnetic configuration can be expanded as Any two magnetic configurations on the same magnetic sites can be compared by an overlap, which we define as Lastly, notice that each CMP carries a definite order and irreducible representation (irrep). Additionally, CMPs of same order and irrep can be enumerated by a label y. This is a convention to write for instance |n(6 T 2u ; y) , where the CMP labeled n is the y-th CMP of 6th order and irrep T 2u . We recall that the 6th order multipole is called 64-pole in the 2 l -nomenclature.

B. Setup for LSDA
The LSDA calculations are performed by the Vienna Ab initio Simulation Package (VASP) in version 5.4 [15,23,24] and the flags are set appropriate to noncollinear LSDA calculation including spin-orbit coupling as described in the Supplemental Material [29]. The pseudopotentials were chosen such that d-electrons in transition metals and f -electrons in Lanthanoids and Actinoids are treated as valence electrons. The default exchange correlation functional, i.e. generalized gradient approximation [31], is used.
The VASP input is created by the aid of the Python Materials Genomics (pymatgen) package [32]. In particular, we used subroutines based on spglib [33]. The magnetic configurations are created by a code authored by M.-T. Suzuki, which employs the TSPACE library [34].

III. RESULTS AND DISCUSSION
In this Section, we want to explore the following two main aspects: (i) Is the CMP expansion a physically meaningful description of magnetic configurations? Namely, here the premise for a physically meaningful description constitutes that naturally occurring magnetic configurations can be characterized by one or few symmetrically related CMPs. It can be understood in the same sense as atomic orbitals are a meaningful basis to describe electrons bound to a free atom, i.e. the probability distribution of one electron is described by one or few degenerate atomic orbitals. In fact, this analogy extents to molecular orbitals in a complex, where the underlying spherical harmonics are symmetrized according to site symmetry.
(ii) Can LSDA predict the most stable magnetic configuration by the aid of an exhaustive list of candidate magnetic configurations for a given crystal? In fact, the predictive power of the combination of the CMP expansion and LSDA (CMP+LSDA) ought to be seen as a composition of the following issues: (a) Is there evidence to assume that the list of candidate magnetic configurations generated by the CMP basis is exhaustive?
(b) Can the experimentally determined magnetic configuration be found among all LSDA results? Note that the similarity between two magnetic configurations is expressed by the overlap defined in Equation (12). In addition, we compare the magnetic space group, which crucially influences physical properties.
(c) Can LSDA correctly assign the lowest total free energy to the experimental magnetic configuration?

A. The investigated materials and workflow
After preluding these questions, let us start by focusing on the experimental data found on MAGNDATA [35]. This commendable collection of meticulously gathered neutron diffraction measurements and other measurements, e.g. optomagnetic response, is still growing and by no means complete. The MAGNDATA entries used in this study were personally double-checked with the experimental references  and the specific compounds are listed in the Supplemental Material [29].
These materials explicitly contain transition metals, Lanthanoids and Actinoids with on-site magnetic moments and most data entries are fully AFM or show only weak ferromagnetism. The magnetic configurations considered here possess zero propagation vector, which limits the available data to about 400 entries in MAGN-DATA. Moreover, entries corresponding to duplicates in respect to higher temperature, pressure or external magnetic field phases are excluded from this study. Finally, some large unit cells are omitted for efficiency reasons.
The still evolving nature of this database inclined us to take a differentiated perspective on each entry: For some materials the size of the magnetic moment is well-determined, while the magnetic order could not be uniquely identified. And conversely, some materials have a well-known symmetry, despite the lack of an exactly determined size of the magnetic moment. Therefore, in this study a total of 131 materials are analyzed, albeit they are distinguished in 122 entries with known magnetic order and 116 entries with known on-site magnetic moment. Figure 1 (a) presents the number of CMPs needed to describe a magnetic cluster featured in the experimental magnetic configuration over the total number of degrees of freedom in the corresponding magnetic cluster. Here, a non-zero CMP component is a so-called active CMP in analogy to the terminology used w.r.t. irreducible representations. The number of degrees of freedom per cluster is naturally equivalent to the order of the CMP basis.
The data shown in Figure 1 (a) comprises 162 magnetic clusters in 122 materials, among which 69 are classified to be collinear, 53 are noncollinear. In particular, 10 are coplanar and 43 are noncoplanar, as indicated by the color of the circles. Meanwhile, the size of the circle indicates the rate of occurrence.
A well-chosen basis is able to express a configuration in terms of few non-zero components. In this regard, remarkably 48.77% of all clusters are characterized by a single active CMP. And only 6 clusters, i.e. 3.70%, of the clusters in the experimental configurations are linear combinations of more than three CMPs.
The construction of the CMP basis, as visualized in the Supplimental Material [29], might intuitively wake the expectation that the number of active CMPs per cluster for a collinear magnetic structure is equal or less than three. Nevertheless, that could not have been generally expected for the noncollinear case. This intuition is empirically confirmed in Figure 1 (a), where all collinear circles are as expected reported below three active CMPs. In the case of noncollinear magnetic configurations, on the other hand, ≤ 3 contributing CMPs per cluster strongly suggests that the basis is particularly well-chosen. Thus, the CMP expansion of experimental configurations in Figure 1 (a) establishes the CMP basis to be a particularly suitable basis. The pie charts in Figure 1 give an overview of the composition of all 131 materials. In particular, Figure 1 (b) shows the orbital character of the valence electrons on the magnetic site. The majority of the materials features transition metals with emerging d-orbital magnetism, while a minority of 25% observes f -orbital magnetism. Secondly, the pie chart in Figure 1 (c) presents the underlying Bravais lattice and fortifies a balanced mixture comprising of all lattice types.
After we have discussed the known experimental properties, let us move on to setting up a predictive scheme. In Figure 2 the computational workflow is organized in four steps: input, setup, calculation, and analysis. The input is taken in form of (magnetic) CIF files [157] from the database MAGNDATA.
Step 2 in Figure 2, the setup, includes reading the magnetic CIF files, creating the list of candidate magnetic configurations and writing the input files for VASP by the aid of pymatgen. Crucially, in this step the CMP basis is obtained as described in Section II A, which does not require the experimental magnetic configuration as an input, but merely the choice of magnetic clusters.
We presume the following heuristic rule holds: The magnetic ground state favors either pure CMPs or linear combinations of CMPs that combine equally weighted CMPs of the same order and same irrep.
This heuristic rule prompts us to extend the list of inital candidate magnetic configurations by linear combinations of same order and same irrep. Neglecting linear combinations of pairs yields (Y − 1)Y additional guesses, for Y being the number of CMPs with same order and same irrep.
In the case of more than one magnetic cluster, d ≥ 2, this would lead to too many additional guesses. For the 73 materials in concern, where d ≥ 2, we chose to combine only the exact same multipole projected onto a different magnetic cluster. In other words, the linear combination of CMPs with same order, same irrep and same y is taken, c.f. the last paragraph of Section II A. Now this similarly leads to (Y − 1)Y additional guesses, but Y is the number of CMPs, which are distinct only w.r.t. c j .
Step 3 in Figure 2, the VASP calculation, is performed as described in Section II B. The total number of LSDA calculations necessary is equal to the number of candidates. The list of candidates is composed of in total d j 3N (cj ) CMP basis magnetic configurations and accordingly many times (Y − 1)Y additional guesses. This amounts to a total of 2935 calculations including all 131 materials in this study.
Step 4 in Figure 2, the analysis, involves determining characteristic quantities of each calculation. First, all possible domains of the converged magnetic configuration are computed. To that end each space group operation combined with time reversal operations ±1 is applied, which leads to either (a) covering the magnetic configuration and thus the operation is element of the magnetic space group, or (b) a new magnetic domain. Considering the set of operations that leave the magnetic configuration invariant, we determine the magnetic space group devising the IDENTIFY MAGNETIC GROUP application on the Bilbao Crystallographic Server [158].
All calculations of a given material and their domains are cross-checked with each other in order to filter how many distinct magnetic configurations and, thus, distinct local minima in the LSDA total free energy landscape have been identified. Quantities such as the total free energy and the size of the magnetic moment per site are averaged over all calculation corresponding to the same local minimum. The calculation with the lowest total free energy among all LSDA calculations of a given material is the CMP+LSDA global minimum.
Note that the list of candidates created as discussed in step 3 is not free of duplicates corresponding to different domains of the same magnetic configuration. In the Supplemental Material [29] the CMP basis for YMnO 3 is constructed. Then, linear combinations of CMPs and magnetic domains are eluded by hands of that example. The candidates corresponding to different domains could be excluded to avoid unnecessary numerical cost. This amounts to a total of 2313 unique calculations for all 131 materials in this study, which comprise of 35.75% additional guesses.
To conclude step 4 in Figure 2, all possible domains are considered when computing the overlaps of (i) the experimental and the initial candidate's magnetic configuration, O exp,init , (ii) the experimental and the converged LSDA calculation's final magnetic configuration, O exp,fin , and (iii) the initial candidate's and the converged LSDA calculation's final magnetic configuration, O fin,init , as defined in Equation (12).
In total this study identifies 2005 CMP+LSDA local minima starting from 2313 unique candidates. As mentioned, we performed 2935 including some redundant candidates in this study. Instead of excluding these redundant candidates that correspond to different domains of the same magnetic configuration, we used them to statistically analyze the reproducibility. In a nutshell, the reproducibility is the probability to converge to the same local minimum, when repeating the LSDA calculation. More details are described in the Supplemental Material [29]. In this study the reproducibility reaches 0.79 on a scale from 0 to 1, where 1 refers to perfect reproducibility.

B. The performance of candidate magnetic configurations
The high computational cost is justified, only if the list of candidates can be expected to be exhaustive. Let us recall that the CMP basis defined in Equation (5) spans the space of all possible magnetic configurations. Each CMP is characterized by its order and irrep. First, we want to argue that the candidate's irrep is likely to prevail throughout the LSDA calculation. As the CMP basis is complete and, thus, any irrep that could be active in a given system explicitly appears in the CMP basis, the former corroborates that the CMP basis is a good starting point.  Figure 3 (a) shows a histogram of the overlap of the final magnetic configuration and the initial candidate, O fin,init . In particular, O fin,init ≈ 1 corresponds to the candidate's magnetic configuration remaining almost identical during the iterations. In that case, the candidate appears to be in close vicinity to a local minimum in the total free energy landscape of LSDA. We see that the uppermost bin, with 46.54% of all calculations, accounts for more calculations than any other bin.
On the other hand, if the candidate does not correspond to a minimum in the total free energy, the calculation is expected to yield a small overlap: O fin,init 1. If the system converges to a magnetic configuration, which is a linear combination of the inital candidate and another magnetic configuration, a finite O fin,init occurs.
There is a related scenario in which the system con-verges to a magnetic configuration that is of the same irrep, but does not including the CMP of the initial candidate. That case can be characterized by O fin,init ≈ 0 and σ irrep = 0, warranted the definition of the variance of the irrep reads Here, σ irrep is defined such that, if the same irreps appear with the same weight in the candidate's CMP expansion and in the CMP expansion of the converged calculation, then σ irrep = 0. In a nutshell, |M n | indicates to what percentage the n-th CMP contributes to the expansion and B nn is a boolean giving zero weight to equal irreps. The colorbar in Figure 3 (a) corresponds to σ irrep defined in Equation (13). The variance of the irrep is less than 10%, σ irrep < 0.1, in 78.16% of all LSDA calculations. In other words, the inital irrep is highly likely to be active in the final magnetic configuration.
The inset of Figure 3  As a more general statement, we have shown that the candidate's irrep is statistically likely to prevail throughout the LSDA calculation. Conversely, the most stable magnetic configuration is less likely to be found, if the irrep is not among the list of candidates. Hence, creating a list of candidates building upon the CMP basis is an efficient solution to assure all possible irreps are among the candidates.
From this point of view, it seems unnecessary to introduce additional guesses as candidates that are equally weighted CMPs of same order and same irrep. However, using the experimental data as a guide once more, the advantages of including additional guesses into the list of candidates becomes clear. Figure 3 (b) presents the maximum overlap of the CMP basis and the experiment, max all O init,exp . The histogram shows a probability density strongly peaked close to one. Additionally, there are side peaks at 1/3 and 1/2. This bias towards 1/3 and 1/2 can be appreciated when considering the aforementioned heuristic rule once again.
Namely, the magnetic ground state favors either pure CMPs or linear combinations of CMPs that combine equally weighted CMPs of the same order and same irrep. at each CMP order CMPs basis configurations occur in sets of 1, 2 or 3 configurations in the expansion. Hence, favored linear combinations projected onto a CMP basis configuration are prone to yield overlap of 1, 1/2 or 1/3.
In comparison, Figure 3 (c), displays the maximum overlap of initial candidate and the experiment, max all O init,exp , w.r.t. the complete list of candidates, which contains the CMP basis configurations as well as additional guesses. The introduction of additional guesses, following the heuristic rule, can effectively avoid side peaks at 1/3 and 1/2 and thus takes into account linear combinations common in materials existing in nature.
As Figure 3 (a) showed, most magnetic configurations remain close to the initial magnetic configuration. Therefore it is paramount to start from an exhaustive list of magnetic configurations.
The dark blue and light blue colors in Figure 3 (b) and (c) indicate, that the magnetic space group found experimentally is identical to the magnetic space group of the candidate or not, respectively. Considering all candidates, as in Figure 3 (c), 117 of 122 magnetic space groups agree. This is an improved agreement rate compared to considering only the CMP basis, as in Figure 3 (b), where 110 magnetic space groups agree. It is noteworthy that some magnetic space groups only enter the list of candidates through the additional guesses.
A final argument in favor of introducing additional guesses is that in total we find 655 of 2005, hence 32.67%, of the local minima in the LSDA energy landscape only thanks to the additional guesses. Even among the CMP+LSDA minima with the minimum total free energy 23 are thanks to the additional guesses, as well as 15 of the (local) minima most similar to the experiment. Therefore, with the collection of arguments mentioned above, we have justified expectation that the list of candidate magnetic configurations is exhaustive. In the following, let us investigate whether the experimentally determined magnetic configuration is present among all LSDA results and how we might predict the likely experimental magnetic configuration for an unknown material.

C. Analysis of CMP+LSDA local minima
Following the workflow in Figure 2 all final LSDA results are scrutinized for their similarity. Some LSDA results correspond to the same local minimum in the LSDA free energy landscape and as such they are grouped in CMP+LSDA local minima.
The overlap of each CMP+LSDA local minimum with the experimental magnetic configuration is computed according to Equation (12). The CMP+LSDA minimum that yields the maximum overlap with the experiment max {loc min} O fin,exp (MaxOExp) is termed to be the most similar CMP+LSDA local minimum to the experiment. A worthwhile run should yield max {loc min} O fin,exp ≈ 1, entailing that MaxOExp is indeed very similar to the experiment. Additionally, the magnetic space group (mspg) should agree with the experimentally detected symmetry.   [144]. While the parent spg is R3c (167) the small tilting results in P1 (2.4) for the mspg. In the CPM expansion, the experimental configuration is described by two CMP basis configurations of order 5. However, they do not observe the same irreducible representation. In particular, the main contribution is A 1g and the tilting is due to contributions of E u . In CMP+LSDA the most stable configuration is pure A 1g without any tilting. So that, although the overlap max {loc min} O exp,vasp = 0.9658, the mspg predicted by CMP+LSDA is R3c (167.103) not P1 (2.4) as found experimentally.
In total 84.43% among MaxOExp yield the correct mspg. This is to say that neither the overlap nor the mspg alone are a sufficient criterion whether the experimental configuration is correctly predicted or not.
In comparison, only 16.17% of all CMP+LSDA minima yield the experimental mspg. However, for 90.16% of the materials at least one CMP+LSDA minima yields the experimental mspg. As mentioned among MaxOExp 84.43% yield the experimental mspg.
Another characteristic CMP+LSDA minimum is the CMP+LSDA global minimum, which observes the minimum total free energy in LSDA. Among all CMP+LSDA global minima only 37.70% yield the experimental mspg. This shows that the mspg of the CMP+LSDA global minima is more likely to agree with the experimental mspg than a random CMP+LSDA minimum, but the CMP+LSDA global minima is not adequately predicting the mspg.
Let us continue by analyzing the LSDA total free energy of the CMP+LSDA minima in more detail. Each CMP+LSDA minimum is attributed one or more LSDA results, as multiple candidates might converge to the same minimum. An average over these attributed LSDA results leads to the material dependent and magnetic configuration dependent total free energy of a specific CMP+LSDA minimum F lm . The CMP+LSDA global minimum observes the minimum total free energy F min .
In order to compare the total free energy across materials, we take a normalized relative total free energy that reads Here, N is the total number of degrees of freedom, i.e. the sum of the oder of basis over all clusters that observe a magnetic moment in LSDA in that material. Figure 4 (b.1) and (b.2) present the distribution of CMP+LSDA minima over the normalized relative total free energy of materials featuring d-orbital magnetism and f -orbital magnetism, respectively. The energy scale is logarithmic in units of meV. And the lowermost bin, representing the CMP+LSDA global minima, would theoretically lie precisely at zero. However, for the obvious practical reasons, namely that log(0) → −∞, it is added at the lower edge. The remaining bins represent the distribution of CMP+LSDA local minima ρ d/f,tot . A key question is, whether MaxOExp tend to be close to the total free energy minimum.
In Figure 4 (b.1) and (b.2) the color intensity classifies all CMP+LSDA minima according to agreement/disagreement with the experimental mspg labeled by o/x, respectively. Additionally, the minima are classified according to being MaxOExp or not. Overall the total free energy distributions ρ d/f,tot span across many orders of magnitude. Albeit, ρ f,tot is more concentrated in the energy range 1 meV up to 1000 meV.
The data shows that in total 43 of 122 (35.25%) of the CMP+LSDA global minima coincide with MaxOExp. Hence, the magnetic configuration with the minimum total free energy in this study does not, at this point, identify the expected experimental configuration. Nevertheless, MaxOExp might tend towards smaller total free energy. In order to gain more insight, we ask if MaxOExp data points follow the same distribution as an arbitrary local minimum in ρ d/f,tot .
Two distributions can be compared in terms of a Q-Q plot [159], where the x-axis represents the quantile of the reference distribution and the y-axis represents the quantile of the sample distribution. Let us define the quantile, q s/r for a sample/reference distribution of local minima {lm k }, where k = 0, .., K − 1 and the local minima (lm) are ordered by F lm k ≤ F lm k+1 . The k/(K − 1) quantile q k is given by Hence, the 0.5 quantile is simply the median value and the 0.1 quantile is the point that divides the distribution such that 90% of the local minima have greater total free energy. The Inset of Figure 4 (b.1) shows the Q-Q plot comparing quantiles of ρ d,MaxOExp , as the sample distribution, with ρ d,tot , as the reference distribution. For each data point in the smaller sample distribution the quantile is computed, as explained above. Subsequently, q d,MaxOExp is juxtaposed against q d,tot .
If the two datasets are sampled from the same underlying distribution ρ d,MaxOExp = ρ d,tot , all points align on the median. The quantile is defined on the same axis as the original distribution, i.e. q d,MaxOExp and q d,tot are defined on (F lm k − F min )/N .
The Q-Q plot in the inset of Figure 4 (b.1) shows significant deviation from the median. Indeed, the slow incline up to approximately 10 meV reveals an accumulation of MaxOExp towards lower free energy. For dorbital magnetism we find 77.66% of MaxOExp below 1 meV. On average each material has 4.45 CMP+LSDA local minima below 1 meV. In particular, in this dataset the material with the maximum number of CMP+LSDA local minima has 18 minima below 1 meV. This shows that CMP+LSDA successfully narrows down the possible magnetic configurations for a new material featuring d-orbital magnetism to a handful of CMP+LSDA local minima, that are highly likely to be close to the experimental observation.
The inset of Figure 4 (b.2) shows the analogous Q-Q plot for f -orbital magnetism. Here, the quantiles basically align on the median suggesting that ρ f,MaxOExp = ρ f,tot . Moreover, for f -orbital magnetism we find only 32.43% of MaxOExp below 1 meV. Although in case of f -orbital magnetism the consideration of the total free energy seems to fail in narrowing down the number of possible magnetic configurations, at least the CMP+LSDA run itself proposes a set of 10 − 15 possible magnetic configurations.
The presented data opens a gateway to identifying a handful of magnetic configurations as CMP+LSDA local minima for a given material among which the experimentally stable magnetic space group and exact configuration is highly likely to be found. Yet it has not been possible to uniquely identify the ground state based on the LSDA total free energy. Although CMP+LSDA yield local minima with the experimental mspg and local minima with large overlap with the experimental magnetic configuration, LSDA fails to assign a low total free energy compared to other local minima.

D. The magnetic moment per site
Besides the magnetic configuration, the size of the onsite magnetic moment crucially influences the magnetic properties of a material. Hence, it is interesting to ask, if the magnetic moment estimated by LSDA is close to the experimentally determined magnetic moment per site. In the literature [160] it is well-known that complexes containing first row transition metals with open 3d orbitals are dominated by crystal field splitting. This is referred to as strong field regime. Further, the ground state of complexes containing Lantanides with open 4f orbitals are dominated by spin-orbit coupling. Complementary, this is referred to as weak field regime. Let us explore the implications by looking closer at the element-dependence of the on-site magnetic moment. Figure 5 presents the on-site magnetic moment averaged over sites within one magnetic cluster as a function of elements sorted by increasing no. of electrons. In particular, the average magnetic moment per site reads and, thus, the average is taken within each magnetic cluster c j , only. The columns show the case of 3d-orbital magnetism and 4f -orbital magnetism, respectively. Figure 5 (a.1) gives an overview of the experimental results µ exp for 3d-orbital magnetism. We see that within compounds featuring the same magnetic element vastly different on-site magnetic moments are reported. This is referred to as compound dependence in the following discussion. Overall, the maximum on-site magnetic moment per element frames a dome shape with a clear maximum at Mn closely followed by Fe. In comparison, Figure 5 (a.2) shows the on-site magnetic moment µ th predicted by CMP+LSDA. Here, µ th is taken to be the magnetic moment of the magnetic configuration with MaxOExp, which has the most similar magnetic order compared to the experiment. We can see very good agreement in the overall tendency between experiment and CMP+LSDA.
A strong crystal field represents a real and time reversal invariant perturbation that forces a real-valued ground state which effectively quenches the orbital angular moment operator (L ≡ 0) as discussed in many text books, see e.g. Ref. [161]. Therefore the spin contribution alone is expected to constitute the on-site magnetic moment. Fortunately, in contrast to the experiment the numeric calculation grants direct access to the spin contribution µ s,th and the angular momentum contribution µ l,th to the on-site magnetic moment µ th = µ s,th + µ l,th .
(18) Figure 5 (a.3) presents the absolute values µ s,th and µ l,th . The data clearly confirms that the angular momentum is almost entirely quenched in LSDA. Only for the heavier elements, where spin-orbit coupling becomes more relevant [162], a small contribution is given by µ l,th . In other words, LSDA supports that for compounds with more than half-filled 3d bands the angular momentum is only partially quenched.
The dominant µ s,th can be directly compared to the spin-only magnetic moment in the ionic limit. It is computed within the Russel-Saunders (or L-S) coupling scheme and is given by with spin quantum number s for the total spin operator S. The total spin S of the electronic configuration 3d n with n electrons is essentially constructed by following Hund's first rules. Albeit in real complexes the electron configuration can be in the high spin (hs) or the low spin (ls) configuration depending on the crystal field strength compared to the intra-orbital Coulomb repulsion. This yields different spin-only magnetic moments µ (3d) s,ion in the ionic limit for electronic configurations of the form 3d n (hs/ls).
In Figure 5 s,ion is displayed as a reference for various possible electronic configurations. Here, we assumed octahedral complexes for the crystal field splitting. the The maximum magnetic moment is consistent with the experiment and CMP+LSDA calculation realized for Mn 2+ or Fe 3+ in the ionic limit. Additionally, the ionic limit already hints towards possible reasons for the observed compound dependence. Namely, we expect the formal oxidation state and the crystal field strength to introduce compound dependence. Further compound  dependence arises due to the exact symmetry including small distortions as introduced by the Jahn-Teller effect and the choice of ligands via the nephelauxetic effect, which describes the delocalization of metal electrons through covalent bonds with the ligands.
Let us now move on to the case of compounds featuring Lanthanides shown in the right column of Figure 5. As mentioned, in the weak field regime spin-orbit coupling is strong compared to the crystal field effect. Therefore the orbital angular momentum operator L cannot be neglected and the magnetic moment is computed in the j-j coupling scheme in terms of the total angular momentum J . In the ionic limit, the electronic ground state can be determined following all three Hund's rules [163] for a given shell configuration 4f n with n electrons. The magnetic moment in terms of the total angular momentum quantum number j then reads with the Landé g-factor (g j ). Representative, we compute the magnetic moment µ (4f ) III-ion for all 3+ ions. Note that in fact, Eu 2+ for instance is expected to resemble Ga 3+ because both have a 4f 7 electronic configuration. Figure 5 (b.1) shows the experimental results µ (4f ) exp for 4f -orbital magnetism in comparison to µ (4f ) III-ion . Similar to the 3d-orbital magnetism, different compounds featuring the same magnetic element observe vastly different µ (4f ) exp , however the origin must be different as we will see. A comparison to the CMP+LSDA results presented in Figure 5 (b.2) shows good agreement of the overall characteristic behaviour. In both, experiment and CMP+LSDA, the magnetic moment is just below the ionic limit and a small (large) dome forms in the less (more) than half-filled region.
Noticeably, the compound dependence in the CMP+LSDA results is reduced compared to the experiment. By a more detailed analysis of the experimental data, the compound dependence in 4f -orbital magnetism is revealed to arise when long-range order cannot be established very well experimentally. LSDA naturally assumes a well-established long-range order by design as it is a zero temperature method. Specific cases are considered in the discussion of Figure 6. l,III-ion in the ionic limit for 3+ ions: Prominently, the destructively (constructively) coupling for less (more) than half-filling is confirmed and visualized. Further, the spin contribution µ (4f ) s,th very closely aligns with the ionic limit. This can be expected as 4f electrons barely delocalize by covalently bonding with the surrounding ligands. The orbital contribution µ (4f ) l,th shows a clearly reduced value compared to µ (4f ) l,III-ion . This might be interpreted as partial quenching of L in LSDA, which is supported by the observation that the reduction of µ (4f ) l,III-ion is stronger for lighter elements. So far it has become clear that there is no systematic overestimation of the on-site magnetic moment by CMP+LSDA. However naively one might anyways expect a general underestimation due to the lack of treatment of strong electronic correlation effects in LSDA, albeit strong electronic correlation is expected in particular in 3d and 4f -bands. As we will see in the following, the data defies this general expectation of an underestimated on-site moment. To this end, let us compare µ th and µ exp compound-wise, or rather cluster-wise for all compounds. Figure 6 juxtaposes the average magnetic moment per site µ th of the magnetic configuration with MaxOExp and the experimentally measured magnetic moment per site µ exp . If for a magnetic cluster µ th ≈ µ exp , the data point is in close vicinity to the median and the size of the magnetic moment per site is well-estimated. In Figure 6 (a), each cluster c j is represented by a star, whose color indicates the orbital character of the magnetic site and the number of points indicates which magnetic element forms the cluster. For instance, the 5-pointed dark red star corresponds to a Mn-cluster, since Mn atom has five 3d electrons. At first sight, there is no general over-or underestimation seen in the scatter plot.
Moreover, the data suggests that the uncertainty of LSDA is reflected in the absolute deviation of |µ th −µ exp |, rather than some relative deviation of the magnetic moment |µ th −µ exp |/|µ th +µ exp |. Indeed, 51.90% of the magnetic moments are within ±0.5 µ B , and beyond 77.22% obey |µ th − µ exp | ≤ 1 µ B . Concomitantly, in the small magnetic moment regime, that is approximately µ 2 µ B , no reliable prediction is possible. In the mid to high magnetic moment regime, on the other hand, a mostly accurate prediction is made.
There is an accumulation of 3d data points within 2 µ B < µ th < 5 µ B , whose center of mass closely aligns with the median. However, an apparent lack of precision leads to a wide spread around the median. Despite another accumulation of 4f data points in the range of 6 µ B < µ th < 10 µ B showing similarly a high accuracy with a center of mass near the median, we also see many outliers with 4f -orbital character across the entire range of the on-site magnetic moments.
The specific group of three outliers at 4 µ B < µ exp < 5 µ B and 7 µ B < µ th < 8 µ B correspond to Er-clusters in Er 2 Sn 2 O 7 , Er 2 Ru 2 O 7 and Er 2 Pt 2 O 7 , listed from left to right. In the ionic limit, the ground state electronic configuration of Er 3+ is 4 I 15/2 with the Landé g-factor (g j ) of 6/5. Therefore, µ (Er) III-ion is estimated to be 9.58 µ B using Equation (20). We see that µ th of the three outliers are considerably less than µ III-ion . In fact, the three outliers are known candidates for realizing a spin liquid phase due to the presence of magnetic frustration, as described in Ref. [118], [71] and [119] and hence present highly non-trivial cases.
The two outliers with µ th > 9 µ B correspond to Hoclusters. Both data points are contributed by the same material HoMnO 3 , which contains two inequivalent Hosites on top of a Mn-cluster. The latter orders at T = 78.5 K and is well-estimated by CMP+LSDA with µ (Mn) exp = 3.32 µ B and µ (Mn) th = 3.47 µ B . On the other hand, experimentally ordering of the two Ho-clusters is subject to controversy [130,[164][165][166][167]. It seems unclear from an experimental perspective whether one or both Ho-sites order even down to approximately 2 K. Generally, the long range ordering of magnetic moments on Ho-sites is suggested to occur at much lower temperature compared to Mn-sites. As mentioned above, a strict comparison of the LSDA result to µ exp is inappropriate in the case that proper long-range ordering cannot be established experimentally. Nevertheless, LSDA can be compared to the ionic limit, similar to the discussion on the three materials containing Er. The ground state electronic configuration of Ho 3+ is 5 I 8 , which yields µ exp is inappropriate. In Figure 6 (b), again µ th and µ exp are compared, but additionally the color indicates whether or not the compound is expected to be frustrated. Here, the expectation of frustration is based on whether nearest neighbors form rings of odd number of magnetic sites. Assuming AFM coupling this geometrically leads to magnetic frustration. Hence, we take advantage of the database being specifically focused on antiferromagnets. Furthermore, rings of even number of magnetic sites could potentially also yield a magnetically frustrated system, if the AFM coupling is anisotropic, such as in the Kitaev model. We hence note, that the definition of expected frustration used here is imprecise and only suitable for a quick superficial classification. Figure 6 (b) shows that indeed the well-estimated 4fclusters in the large magnetic moment regime are not expected to feature magnetic frustration. The discussed group of three outliers on the other hand are expected to be frustrated. Data points with 4f -orbital character in the small magnetic moment regime µ th < 2 µ B are likewise expected to be magnetically frustrated and are not particularly well-estimated. Although, we expect that µ th is overestimated when the system is frustrated, many clusters that are expected to be magnetically frustrated are not necessarily overestimated. And some outliers are-at least in the approximate definition employed here-not expected to be frustrated. However, as we have seen for HoMnO 3 there might be other non-trivial phenomena preventing a proper long-range order. Hence, the geometrically expected magnetic frustration is not a sufficient indicator for overestimation of the magnetic moment. Figure 6 (c) displays the filling on a colormap from 0 to 1, where 0.5 correspond to half-filling. Here, the filling is defined as the ratio between the number of d or f electrons in each magnetic atom and the number of orbitals. For the number of electrons, we consider the charge neutral state, i.e., the ionized state is not taken account. Less (more) than half-filled 4f and 5f -clusters appear in the underestimated (overestimated) region. Figure 6 (d) addresses the number of magnetic clusters present in a specific compound. The data points corresponding to single cluster (red), and multiple clusters (blue) appear to be evenly distributed. Let us divert the attention towards data points with µ th ≈ 0. It should be noted that these are not paramagnetic solutions. Two scenarios can yield µ th ≈ 0: Either another cluster bears most of the on-site magnetic moment, or the spin contribution to the magnetic moment µ s is canceled by the orbital contribution to the magnetic moment µ l . Figure 6 (e) shows the normalized orbital contribution in LSDA to the total magnetic moment µ th = |µ s + µ l |. Below the median in the small magnetic moment regime, indeed many clusters with less than half-filled orbitals observe µ l /(|µ l | + |µ s |) ≈ 0.5. In these instances, µ s and µ l adopt opposing signs and thus the contributions in fact cancel. Clusters of heavier Lanthanides are wellestimated solely as a result of including µ l . Considering, once more Figure 5

IV. CONCLUSION AND OUTLOOK
This study is a benchmark of an ab initio prediction of the magnetic ground state using a novel approach termed CMP+LSDA. This scheme devises a combination of the cluster multipole (CMP) expansion and the local spin-density approximation (LSDA) for noncollinear magnetism. We find that materials existent in nature are well-described in terms of only few CMPs and infer the CMP expansion basis to be a suitable basis for magnetic configurations. Additionally, the experimental data suggests that the magnetic ground state favors either pure CMPs or linear combinations of CMPs having the same expansion order and same irreducible representation. Guided by this heuristic rule an exhaustive list of initial candidate magnetic configurations for LSDA calculations is created.
A high-throughput calculation of 2935 LSDA calculations using VASP led to a handful of CMP+LSDA local minima corresponding to different possible magnetic configurations for each material. 90.16% of materials yield the experimental magnetic space group for at least one of the CMP+LSDA local minima. Furthermore, the maximum overlap between the experimental magnetic con-figuration and the CMP+LSDA local minima exceeds 0.75-with 1 corresponding to equivalence-in 70.99% of all materials.
An ab initio prediction of the most stable magnetic configuration in the experiment is guided by a comparison of the total free energy in LSDA of the the possible magnetic configurations for each material. In particular, the local minimum with the larges overlap with the experiment (MaxOExp) is expected to yield the lowest total free energy. Indeed, for materials featuring magnetic sites with d-orbital magnetism, MaxOExp is in great majority of the cases less than 1 meV above the so-called CMP+LSDA global minimum. On the other hand, the same could not be confirmed for f -orbital magnetism. In fact, MaxOExp for f -orbital magnetism shows no tendency towards lower total free energy. The implementation of LSDA used in this study did not necessarily assign the lowest total free energy to the local minimum with the larges overlap with the experiment. Nevertheless, we want to emphasize that CMP+LSDA succeeded to narrow down the number of possible magnetic ground states. This is achieved in parts thanks to a list of candidate magnetic configurations that is tailored to account for details of the symmetry of the crystallographic unit cell. With this CMP enables LSDA to identify a feasible number of local minima, that put data screening and AFM material design within reach.
In addition, this study showed that the on-site magnetic moment could be estimated surprisingly well by LSDA. The precision of the predicted magnetic moment is estimated to be roughly ±0.5 µ B . Some outliers arise from a lack of long-range order in the experiment. This can be due to extremely low transition temperatures but also due to magnetic frustration. Despite some explainable outliers, the prediction shows no major systematic over-or underestimation of the on-site magnetic moment in LSDA. In contrast to the experiment, the LSDA calculation grants additional insight into the balance of spin contribution and orbital angular momentum contribution to the total magnetic moment. The first row transition metals prove to be well-described by Russel-Saunders coupling applicable within the strong field regime. The orbital angular momentum is quenched and the spin-only ionic limit can be used as a reference. The case of Lanthanides, on the other hand, is representative for systems in the weak field regime. The on-site magnetic moment is well-described in the j-j coupling scheme. In the end of III, we speculate that LSDA might have slightly overestimates the crystal field effects compared to the strength of spin-orbit coupling. This could explain why materials governed by crystal field splitting-such as the compounds with d-orbital magnetism-are assigned appropriate total free energy by LSDA. Yet, materials governed by spin-orbit coupling-such as Lanthanidesthe experimental magnetic configuration is not assigned the lowest total free energy by LSDA. The balance between spin-orbit coupling and crystal field splitting becomes particularly crucial for lighter 4f -elements and heavier 3d-elements, where the orbital angular momentum is partially quenched. Actually there is a third energy scale, the intra-orbital Coulomb repulsion that might play a crucial role. The approximate inclusion of such an intra-orbital Coulomb repulsion led to the development of LSDA+U [168]. It is widely believed that LSDA+U might improve the theoretical treatment, although it is not well-tested as no large-scale benchmark has been performed.
The starting point of this study is the experimental database MAGNDATA [35]. It conveniently facilitated testing and benchmarking of our ab initio scheme. Generally, experimental databases [169][170][171][172][173][174][175][176][177][178] not only facilitate testing and benchmarking of theoretical methods, but also data mining in the experimentally explored chemical space. Indeed, for some nonmagnetic functional materials an informed search and optimization has led to promising discoveries [179][180][181][182][183][184][185][186][187][188][189][190][191]. However so far, apart from few pioneering works [5][6][7][8][9] that are constrained to specific cases, these breakthroughs in material design have not yet been matched by similar advances with respect to AFM materials. Certainly one of the major obstacles is that compared to databases of crystal structures with more than 200 000 entries, MAGNDATA has to date a modest amount of about 1 000 entries. This is because the experimental determination of the magnetic configuration is much more involved than that of the crystal structure. Given this situation, it is an urgent challenge to construct a large-scale computational database of AFM materials. The presented benchmark provides a crucial step in laying a solid foundation for the construction of such a computational database of AFM materials.

V. ACKNOWLEDGMENTS
We are thankful for useful discussions with Stephan Huebsch, Takashi Koretsune and Saeed Bahramy.
Moreover, we gratefully acknowledge the Center for Computational Materials Science, Institute for Materials Research, Tohoku University for the use of MASAMUNE-IMR (MAterials science Supercomputing system for Advanced MUlti-scale simulations towards NExt-generation -Institute for Materials Research). (Project No. 19S0005) This work was supported by a Grant-in-Aid for Scientific Research (No. 19H05825, and No. 16H06345) from the Ministry of Education, Culture, Sports, Science and Technology, and CREST (JPMJCR18T3) from the Japan Science and Technology Agency. correspond to components of e µ lγ in Ref. [13] except that here and magnetic toroidal (MT) multipoles, respectively. As we do not expand in terms of MT multipoles, those components appear as higher order magnetic multipoles here, and thus completeness is still ensured.