Automated Generation of Commensurate Magnetic Structures based on Spin Space Groups and Graph Theory

INTRODUCTION

A long-standing goal of computational magnetism is to describe magnetic ground states, collective excitations, and phase transitions in agreement with experiment. Effective Hamiltonians make this problem manageable by retaining only the degrees of freedom relevant to the phenomenon of interest. For magnetism, this leads to the development of the spin model that traces back from the Heisenberg model [1, 2, 3] to models with relativistic and higher-order interactions [4, 5, 6, 7, 8], and finally the recent general tensorial model [9]. Once their coefficients have been obtained from first-principles calculations, these models can be transferred into various post-processing calculations. Energy mapping is a central algorithm to extract spin interaction coefficients from the Schrödinger electronic structure. It uses the total energies of a set of magnetic configurations to fit into the energy landscape within a predefined framework of the spin model. The quality of the fit depends heavily on the choice of the configurations.

Symmetric magnetic structures are configurations that remain invariant under symmetry operations belonging to various groups. They are energy-degenerate points that are well-spread over the configuration space. Statistically, sampling symmetric structures can increase the linear stability of the fitting. They are also highly likely to be the local energy minimum, which is accessible via the self-consistent calculation. Even without the fitting procedure, the ground states of most magnets can be tractable via enumerating the symmetry structures and comparing their total energies. Therefore, it is of great significance to list symmetry magnetic structures, especially in an automated workflow so that density functional theory can be implemented in a high-throughput way. Previously, superpositions of spin spiral states were prepared by hand in AMATIS[9], catering only to the magnetic space groups [10]. This way of selection already reproduces a robust fitting of the spin interactions. Here, we aim to develop an automated workflow for generating symmetric structures in the expanded framework of the spin space group, which is important for diverse magnet data construction.

Current methods commonly involve two ways to generate symmetry magnetic structures. One imposes a magnetic space group directly and keeps the structures that obey it [11, 12], such as Bilbao Crystallographic Server and JANA2020 that construct commensurate and incommensurate structures using magnetic superspace groups [10, 13, 14, 15]. The other uses representation analysis within Landau theory [16, 17, 18, 19]. BASIREPS implements Bertaut’s representation analysis [20], and SARAh provides online representation-analysis tools [21]. More recently, spin-space-group symmetry has been incorporated into magnetic-structure generation. [22] constructed spin-symmetry-adapted structures from totally symmetric representations of oriented spin space groups, while [23] developed a symmetry-guided framework for systematically generating magnetic configurations and predicting magnetic ground states. Though these approaches provided important foundations for the recent magnet data structure constructions, there remain a few limitations for each of them. For instance, representation-analysis methods naturally incorporate prescribed propagation ($\mathbf{q}$) vectors, but their conventional formulations are not readily compatible with general spin-space-group symmetries. Conversely, existing spin-space-group-based methods can generate structures with prescribed spin symmetry, but do not directly provide systematic construction under specified single- or multi-$Q$ constraints. Explicit propagation-vector constraints are important because magnetic Bragg peaks may identify the relevant propagation vector(s) without uniquely determining the underlying magnetic structure [11, 24].

We therefore developed SpinGraph, an automated workflow that is compatible with both the aforementioned requirements of the spin space groups framework and the propagation-vector constraints. It uniquely formulates magnetic-structure generation as a graph-constrained problem, enabling an exact treatment of continuous spin-only components and a unified description of real-space, Fourier-space, supercell, and propagation-vector constraints within a single workflow. We demonstrate the method using CrI$_3$ as a representative example, where SpinGraph systematically generates symmetry-distinct magnetic structures under different supercell and $Q$ constraints and reproduces the expected symmetry relations. Using 396 generated configurations for energy mapping, we further recover the few dominant Heisenberg exchange, anisotropic interactions, and higher-order spin interactions consistent with the established magnetic behavior of CrI$_3$.

Results

SpinGraph framework

SpinGraph provides three schemes for calculating symmetry-constrained magnetic structures. The first constrains only the supercell size, the second constrains the propagation vectors of the magnetic structure, and the third constrains both the supercell size and the propagation vectors. The program then decomposes the continuous and discrete parts of the spin space group (SSG), combines them, and attempts to construct SSGs with different symmetry-breaking patterns. For each newly constructed SSG, the program identifies its generators and the infinitesimal generators of its Lie-group part, constructs the graph structure, and solves the linear null space. Finally, zero-total-magnetic-moment constraints or equal-magnetic-moment constraints are added to obtain the magnetic structures.

Spin space group action

An SSG is a symmetry group in which the spin and spatial parts are decoupled, and therefore has more degrees of freedom than a conventional magnetic space group (MSG)[25, 26, 27]. When the spin and spatial parts are locked together, the SSG reduces to an MSG. An operation $g \in G$ is written as

where $R_g$ and $\tau_g$ act on real-space crystallographic coordinates, $U_g \in SO(3)$ rotates the spin, and $\eta_g = -1$ marks an operation containing time reversal [27]. Their product $O_g \in O(3)$ is the operation that acts on the axial magnetic moment. Let $L$ be the set of discrete vectors with integer components, with the basis of the primitive Bravais lattice, and $L_{\mathrm{sc}} \subset L$ the set of supercell translation vectors. For the finite set $R_{\mathrm{sc}} = \{r_i\}$ of magnetic sites in that supercell, the action on a site position and its magnetic moment is

Here, $T_{ij}^{(g)} \in L_{\mathrm{sc}}$ is the supercell translation relating the transformed site to representative $r_j$, and $S_i$ is the axial magnetic moment at $r_i$. We use $q \cdot T$ for the crystallographic phase when the propagation vector $q$ is in reciprocal-lattice units and $T$ is a lattice translation, with the customary factor $2\pi$ included. For the complex Fourier amplitude $S_i(q)$ at site $i$ and propagation vector $q \in Q$, where $Q$ is the selected set of allowed propagation vectors. The reciprocal-space action is

Here, $L^*$ is the set of the reciprocal lattice which is one-to-one corresponds to $L$, $\rho_g$ denotes the reciprocal-space action of $g$, and $q' = \rho_g(q)$ is the transformed propagation vector. For a real magnetic structure, $q$ and $-q$ appear in pairs, and their Fourier components satisfy a complex-conjugation constraint. To treat this conjugation constraint uniformly within the graph structure, the complex Fourier components are realified into six-dimensional vectors. Supplementary Notes 1 and 4 derive the phase and compare the linear and anti-linear conventions for time reversal[28, 29].

Compatible symmetry

The three compatible groups used by the workflow are the supercell stabilizer $G_{\mathrm{sc}}$, the propagation-vector stabilizer $G_Q$, and their joint stabilizer $G_{\mathrm{sc},Q}$, defined by

$G_{\mathrm{sc}}$ denotes the set of operations that leave the crystal structure invariant when the supercell is taken as the reference unit. Similarly, $G_Q$ contains the operations that permute the selected propagation vectors. $G_{\mathrm{sc},Q}$ denotes the set of operations allowed under both the supercell constraint and the $Q$-set constraint. In the following, $G_{\mathrm{c}}$ denotes the symmetry group selected in the calculation, with pure supercell translations in $L_{\mathrm{sc}}$ understood to have been quotiented out. Supplementary Note 2 gives the subgroup proofs[28].

If $G_{\mathrm{c}}^{\circ}$ is the identity component, i.e., the connected subgroup containing the identity element and all symmetry operations continuously connected to it, the discrete symmetry is represented by the finite component group,

For a finite parent group, $G_{c}^{\circ}=\{e\}$ and this quotient is simply the original compatible group. We enumerate the subgroups of $\overline{G}_{c}$ and then restores the compatible spin-only constraints. The subgroups of the discrete component group are enumerated using its Cayley table, and details are provided in Supplementary Note 3 [28]. However, a continuous parent’s exact closed subgroups generally do not form a finite list [30, 31]. The algorithm therefore separates the finite and spin-only parts. For every enumerated subgroup $\overline{K}\leq\overline{G}_{c}$, it selects a spin-only subgroup $\mathcal{K}_{K}$ that is contained in the parent spin-only factor and is normalized by the finite subgroup,

Here and below, $h\in\overline{K}$ is understood through a chosen representative in $K$. For magnetic-structure construction, spin-only subgroups that impose the same constraints on the allowed spin directions are treated as equivalent and need not be solved repeatedly. The corresponding subgroup relations are discussed in Supplementary Note 3 [28].

Graph constraints and solution spaces

For a selected subgroup candidate $K\leq G_{c}$, let $\Gamma_{K}^{\mathrm{disc}}=\{g_{1},\ldots,g_{n_{K}}\}$ denotes the generators of its finite representative. We write $\mathbf{S}$ for a magnetic configuration and $\mathcal{D}(g)$ for the linear action of $g$. If

then the same relation holds for every product of these generators and their inverses. In particular,

If $K$ has a nontrivial identity component $K^{\circ}$, its connected spin-only action is imposed through the Lie algebra

where $\mathcal{L}_{a}$ is the infinitesimal action induced by $\mathbf{X}_{a}$. $SO(2)$ requires one independent generator, whereas $SO(3)$ requires three [32, 33, 34]. Coplanar constraints are imposed directly by projecting the spins onto the selected plane.

The discrete generator actions define either a real-space graph $G_{r}$ or a Fourier-space graph $G_{q}$. Let $\mathbf{x}_{i}$ be the variable at vertex $i$. Each directed edge $e$ stores a matrix $\mathbf{P}_{e}$ relating its source $s(e)$ to its target $t(e)$,

In real space, $\mathbf{x}_{i}=\mathbf{S}_{i}\in\mathbb{R}^{3}$. In Fourier space, it is the $\mathbb{R}^{6}$ form of a complex amplitude. In $G_{r}$, each site $\mathbf{r}_{i}\in\mathcal{R}_{\mathrm{sc}}$ is a vertex, and (2) supplies the edge matrix $\mathbf{O}_{g}$. In $G_{q}$, a vertex is a pair $(i,\mathbf{q})$, and the edge matrix also contains the phase in (3). Using the constructed graph and the transport relations, linear constraints can be obtained by identifying the fundamental cycles of the graph. A graph may contain multiple connected components, and independent degrees of freedom are not shared between different components, whereas vertices within the same component share the same set of independent degrees of freedom. For each connected component, a spanning tree transports a reference variable $\mathbf{x}_{0}$ through the component, $\mathbf{x}_{i}=\mathbf{M}_{i}\mathbf{x}_{0}$, while each remaining edge $e:i\to j$ imposes a cycle-consistency condition [35, 36].

Stacking these conditions together with the continuous and planar constraints gives

The third equation shows that the basis vectors at the other vertices within a connected component can be obtained through the transport relations and are not independent. Supplementary Note 5 gives the group constraints and fundamental-cycle details [28, 37, 38].

Fourier-space restriction and projection

Let $V_{\mathrm{sc}}$ be the space of all magnetic configurations in the selected supercell, and let $V_Q\subseteq V_{\mathrm{sc}}$ contain only the Fourier components in $Q$. In a real-space calculation restricted to $Q$, the projector $P_Q:V_{\mathrm{sc}}\to V_Q$ removes all other Fourier components. Because both the supercell constraint and the $Q$-set constraint must be satisfied, the corresponding symmetry group is $G_{\mathrm{sc},Q}$, and the resulting solution space is

Thus, projecting the real-space solution onto $Q$ gives the same solution space as solving directly in Fourier space. The equivalence is proved in Supplementary Notes 6 and 7 [28].

Although the real-space calculation is performed using $G_{\mathrm{sc},Q}$, the resulting magnetic structures may still contain propagation vectors outside $Q$. Therefore, the basis obtained from the real-space calculation must be projected onto the allowed subspace using the projection operator. Thus, the real-space basis is projected as

Construction of physical magnetic structures

After obtaining the symmetry-allowed basis, additional physical constraints are imposed directly on its global coefficients. A zero-total-spin condition further reduces the allowed coefficient space, and an equal-moment-length condition constrains the remaining coefficients. In contrast, normalizing each site independently may break these relations or introduce Fourier components outside $Q$. The detailed procedure is given in Supplementary Note 8 [28].

The symmetry used to construct a magnetic structure is not necessarily its full symmetry, because a particular solution may have additional operations. Therefore, after a structure is generated, we should check all allowed spatial and spin operations and determine its maximal stabilizer $K_{\max}$, after which the full symmetry operations are verified against the magnetic structure [26, 39, 40, 41].

SpinGraph workflow

The inputs are a crystal structure, the magnetic sites, the spin space group operations, and a supercell, a propagation-vector set $Q$, or both. The output contains the allowed magnetic basis together with the subgroup and graph information used to obtain it. Figure 1 summarizes the following steps.

1. Prepare the input data. SpinGraph reads the CIF, determines the crystallographic symmetry and Wyckoff positions with Spglib, and expands the selected magnetic sites into a supercell [42]. The parent magnetic or spin space group operations are then obtained from Spglib, FindSpinGroup, or a user-supplied SSG file [41]. For Fourier calculations, the program also completes the input propagation-vector set for real conditions.

2. Select the compatible operations. Each parent operation is tested using (4). Depending on the calculation mode, it gives $G_{\mathrm{sc}}$, $G_Q$, or $G_{\mathrm{sc},Q}$. The program verifies the identity, closure, and inverses of the group. For a continuous SSG, the finite component is used here, while the continuous spin constraints are retained separately.

3. Enumerate the subgroups. The finite subgroups are enumerated by cyclic extension of the Cayley table [43, 44, 45]. For a continuous parent, this step is applied to its finite component group, and each finite subgroup is then combined with the spin-only axial, plane-reflection, or trivial constraint allowed by (6). Then the corresponding groups are identified by Spglib or FindSpinGroup.

4. Build the constraint graph. The generators in $\Gamma^{\mathrm{disc}}_{K}$ are applied to all graph vertices. Real-space edges relate magnetic moments on different sites, while Fourier-space edges relate amplitudes at different sites and propagation vectors. Continuous spin symmetries are added as exact local linear constraints as self-loops.

5. Solve the graph constraints. The graph is divided into connected components [35, 36]. A spanning tree, generated by BFS [36], transports a reference vector through each component, and the remaining edges provide a complete set of consistency conditions. Solving these conditions gives the allowed real-space or Fourier-space basis and its free parameters.

6. Apply the additional restrictions. For a joint calculation $G_{\mathrm{sc},Q}$ in the real-space representation, SpinGraph projects the solved basis onto the requested propagation vectors and removes zero or dependent basis vectors. When magnetic structures are generated, optional zero-total-spin and equal-moment constraints are imposed on the complete global coefficient space. Only one common scaling factor is applied, while sitewise normalization is not used.

7. Generate and validate the output. The component bases are combined into the full magnetic basis, and generated moments are transformed to Cartesian coordinates. SpinGraph checks the complete subgroup, the continuous constraints, the prescribed Fourier support, and any enabled magnetic constraints. After equivalent structures are removed, the stabilizer of each magnetic configuration is recomputed.

Magnetic-structure generation and validation

Magnetic-structure generation using a CrI3 monolayer

We compare the three calculation modes using the same monolayer CrI3 structure. It belongs to $P\overline{3}1m$ (No. 162), and the two Cr atoms at the $2d$ Wyckoff position, $\mathbf{r}_1=(1/3,2/3,1/2)$ and $\mathbf{r}_2=(2/3,1/3,1/2)$, are chosen as magnetic sites. The parent symmetry is spin space group 162.157.8.1. For a direct comparison, we use the same embedded spin space group 10.6.1.3 in all three cases.

In $G_{\mathrm{sc}}$, we use a $(3 \times1 \times1)$ supercell and place no restriction on $Q$, and the Fourier decomposition may contain $\mathbf{q}=(0,0,0)$, $(1/3,0,0)$, and $(-1/3,0,0)$. For spin space group 10.6.1.3, the six Cr sites form three independent pairs. Each pair has an in-plane basis direction $\mathbf{u}_{\mathrm{sc}}=(\kappa,1,0)$, where $\kappa=(1+\sqrt{3}+\sqrt{6})/2$, and an out-of-plane direction $\mathbf{e}_z=(0,0,1)$. This gives six real free parameters, as listed in Table 1.

Site	Coordinate	Moment $S_{\mu\ell}$
Cr$_{10}$	$(1/3,2/3,1/2)$	$a_0\mathbf{u}_{\mathrm{sc}}+b_0\mathbf{e}_z$
Cr$_{11}$	$(4/3,2/3,1/2)$	$a_1\mathbf{u}_{\mathrm{sc}}+b_1\mathbf{e}_z$
Cr$_{12}$	$(7/3,2/3,1/2)$	$a_2\mathbf{u}_{\mathrm{sc}}+b_2\mathbf{e}_z$
Cr$_{20}$	$(2/3,1/3,1/2)$	$a_2\mathbf{u}_{\mathrm{sc}}+b_2\mathbf{e}_z$
Cr$_{21}$	$(5/3,1/3,1/2)$	$a_1\mathbf{u}_{\mathrm{sc}}+b_1\mathbf{e}_z$
Cr$_{22}$	$(8/3,1/3,1/2)$	$a_0\mathbf{u}_{\mathrm{sc}}+b_0\mathbf{e}_z$

AMATIS calculations for monolayer CrI3

In $G_Q$, we choose $\mathbf{q}=(1/3,0,0)$, so the effective set is $Q=\{\mathbf{q},-\mathbf{q}\}$. Spin space group 10.6.1.3 has one connected component and four real free parameters. Here, we define $\mathbf{A}=\alpha\mathbf{u}_{\mathrm{sc}}+\beta\mathbf{e}_z$, $\mathbf{B}=\gamma\mathbf{u}_{\mathrm{sc}}+\delta\mathbf{e}_z$, and $\omega=e^{2\pi i/3}$. The amplitudes at $\mathbf{q}$ are given in Table 2, while $S_\mu(-\mathbf{q})=S_\mu(\mathbf{q})^*$. The real-space spin is reconstructed as $S_{\mu\ell}=\sum_{\mathbf{q}\in Q}S_\mu(\mathbf{q})e^{-i\mathbf{q}\cdot\mathbf{T}_\ell}$, where $\mathbf{T}_\ell$ is the translation of cell $\ell$. The conjugate pair makes the reconstructed spin real and removes the $\mathbf{q}=(0,0,0)$ component.

Site	Coordinate	$S_\mu(\mathbf{q})$
Cr$_1$	$(1/3,2/3,1/2)$	$\mathbf{A}+i\mathbf{B}$
Cr$_2$	$(2/3,1/3,1/2)$	$\omega^2(\mathbf{A}-i\mathbf{B})$

In $G_{\mathrm{sc},Q}$, we use the same $(3 \times1 \times1)$ supercell, $Q=\{\mathbf{q},-\mathbf{q}\}$. The three Cr-site pairs found in $G_{\mathrm{sc}}$ remain, but the projection couples their coefficients and reduces the basis from six to four real parameters. Equivalently, the three moments on each Cr sublattice sum to zero, so the $\mathbf{q}=(0,0,0)$ component vanishes. Table 3 lists the resulting real-space basis. The $G_Q$ calculation imposes the restriction directly in Fourier space, while the $G_{\mathrm{sc},Q}$ calculation obtains the same restricted space by projecting the real-space basis, since $G_Q=G_{\mathrm{sc},Q}$ in this case.

Site	Coordinate	Moment $S_{\mu\ell}$
Cr$_{10}$	$(1/3, 2/3, 1/2)$	$-(a+c)\mathbf{u}_{\mathrm{sc}}-(b+d)\mathbf{e}_{z}$
Cr$_{11}$	$(4/3, 2/3, 1/2)$	$a\mathbf{u}_{\mathrm{sc}}+b\mathbf{e}_{z}$
Cr$_{12}$	$(7/3, 2/3, 1/2)$	$c\mathbf{u}_{\mathrm{sc}}+d\mathbf{e}_{z}$
Cr$_{20}$	$(2/3, 1/3, 1/2)$	$c\mathbf{u}_{\mathrm{sc}}+d\mathbf{e}_{z}$
Cr$_{21}$	$(5/3, 1/3, 1/2)$	$a\mathbf{u}_{\mathrm{sc}}+b\mathbf{e}_{z}$
Cr$_{22}$	$(8/3, 1/3, 1/2)$	$-(a+c)\mathbf{u}_{\mathrm{sc}}-(b+d)\mathbf{e}_{z}$

Multi-Wyckoff real-space validation using Fe$_4$GeTe$_2$

This example checks whether the algorithm keeps inequivalent magnetic Wyckoff orbits separate in one calculation. The Fe sublattice of Fe$_4$GeTe$_2$ in a $(1\times1\times1)$ cell has two orbits, Fe1 and Fe2, each containing six sites.

The chosen group $G_c$ has BNS label 166.101. Its real-space graph $G_r$ contains all 12 Fe sites but splits into two six-site components, $C_1$ and $C_2$, for Fe1 and Fe2. Symmetry relates the sites within each component, but no graph edge relates the amplitudes of the two inequivalent orbits. All Fe moments within one component are therefore parallel to the crystallographic $c$ axis, while $a_1$ and $a_2$ remain independent. Table 4 shows the result site by site.

Site	Fractional coordinate	Magnetic moment
Fe$_{11}$	$(0,0,0.88813)$	$(0,0,a_1)$
Fe$_{12}$	$(0,0,0.11187)$	$(0,0,a_1)$
Fe$_{13}$	$(2/3,1/3,0.22147)$	$(0,0,a_1)$
Fe$_{14}$	$(2/3,1/3,0.44520)$	$(0,0,a_1)$
Fe$_{15}$	$(1/3,2/3,0.55480)$	$(0,0,a_1)$
Fe$_{16}$	$(1/3,2/3,0.77853)$	$(0,0,a_1)$
Fe$_{21}$	$(0,0,0.81227)$	$(0,0,a_2)$
Fe$_{22}$	$(0,0,0.18773)$	$(0,0,a_2)$
Fe$_{23}$	$(2/3,1/3,0.14561)$	$(0,0,a_2)$
Fe$_{24}$	$(2/3,1/3,0.52106)$	$(0,0,a_2)$
Fe$_{25}$	$(1/3,2/3,0.47894)$	$(0,0,a_2)$
Fe$_{26}$	$(1/3,2/3,0.85439)$	$(0,0,a_2)$

However, Table 4 is not a prediction of the magnetic ground state. Symmetry fixes the allowed directions but not the values or relative sign of $a_1$ and $a_2$, and it does not compare the energies of the allowed states. That would require an energy calculation or experimental input. Here the result only shows that the program treats the two Fe orbits separately and combines their allowed moments correctly.

Spin interactions in monolayer CrI$_3$

The magnetic configurations generated by SpinGraph are used as input to AMATIS to extract the tensorial spin interactions. We consider monolayer CrI$_3$, a two-dimensional honeycomb insulator with space group $P\bar{3}1m$ (No. 162). The magnetic Cr atoms occupy the $2d$ Wyckoff position, with spin $S=3/2$.

For the rank-2 interactions, the cutoff radius is set to $r_c=13$ Å. This includes the single-ion anisotropy tensor $J_1^{2-1}$ and five two-site tensors up to $J_5^{2-2}$. For the rank-4 interactions, we include two-site tensors in which each site appears twice, up to the second-nearest-neighbor tensor $J_2^{4-2}$. No rank-4 single-ion term is included because $S < 2$. The spin model is fitted using the 396 magnetic configurations generated above. The smallest supercell commensurate with all these configurations is $4 \times2 \times1$, and the total energies of all configurations are therefore calculated using this supercell. The fitted rank-2 and rank-4 tensor components are provided in the Supplementary Note 10 [28], which is consistent with [9].

The extracted interactions are consistent with the experimental description of monolayer CrI$_3$ as an Ising-type ferromagnet[46]. Both the single-ion anisotropy, 0.52 meV per Cr, and the two-ion anisotropy, 0.43 meV per Cr, favor out-of-plane spin polarization. The first three nearest-neighbor Heisenberg interactions are dominant and give a Curie temperature of about 45 K. In addition to the two-ion anisotropy, the symmetric exchange $\Gamma^{ij}$ contains a bond-dependent Kitaev-like interaction. For the nearest-neighbor bond, its preferred plane is the Cr-I-Cr-I rhombus shown by the green plaquette in Figure 2 (a). Although inversion symmetry forbids a net Dzyaloshinskii–Moriya interaction, local inversion-symmetry breaking along second-nearest-neighbor Cr bonds allows antisymmetric exchange interactions of about 0.20 meV.

Discussion

SpinGraph provides a comprehensive workflow for generating commensurate magnetic structures, capable of producing a series of magnetic structures with different symmetry-breaking patterns. This method can simultaneously handle both spin space groups and magnetic space groups, and can properly treat the continuous components of these groups. Furthermore, the algorithm transforms symmetry constraints into graph constraints and constructs the corresponding linear null-space problem for solution, thereby avoiding the need to establish systems of linear equations containing a large number of redundant relations, while also offering good interpretability.

During the generation of magnetic structures, constraints such as zero total magnetic moment and equal moment magnitude can be further introduced to ensure that the resulting structures meet the corresponding computational requirements. The algorithm’s different computational modes share the same underlying graph structure: the real-space mode handles supercell periodicity, the Fourier mode handles a specified propagation-vector set $Q$, and the combined mode applies both types of constraints simultaneously. Consequently, this method can generate magnetic structures that simultaneously satisfy strict propagation-vector constraints and symmetry requirements, thereby covering cases that are difficult for traditional magnetic-structure generation programs to handle directly.

The final result is a space of candidate magnetic structures permitted under the given symmetry conditions, rather than the energetic ground state. The relative stability of the candidate structures still requires further evaluation through spin models, first-principles calculations, or experiments. In the future, spin superspace groups could be introduced to extend this framework to the generation of incommensurate magnetic structures.

METHODS

Symmetry identification and databases

Crystal symmetries and Wyckoff positions are identified using Spglib[42]. When the spin-space symmetry reduces to an ordinary magnetic space group, Spglib is also used to identify the corresponding magnetic subgroup and the magnetic symmetry of the final generated structure. For spin-space group calculations, the parent SSG is constructed from the spin-space group database implemented in FindSpinGroup[41]. The same database and identification routines are then used to identify the generated subgroups and the final spin-space group symmetry of each magnetic structure.

First-principles calculations

The ab initio calculations were performed using the Vienna ab initio Simulation Package (VASP)[47, 48] with the projector augmented-wave (PAW) method[49, 50]. The Cr $3p$, $3d$, and $4s$ orbitals and the I $5s$ and $5p$ orbitals were treated as valence states. A plane-wave cutoff energy of 520 eV and a $2 \times4 \times1$ $k$-point mesh were used for monolayer CrI$_3$. The electronic energy convergence criterion was set to $10^{-7}$ eV, and spin–orbit coupling and noncollinear magnetism were included.

The exchange-correlation functional was treated within the local spin-density approximation (LSDA)[51], together with a Hubbard $U = 0.5$ eV applied to the Cr $d$ orbitals, following Refs. [52, 53]. The experimental in-plane lattice constant $a = 6.72$ Å was used[46, 53]. The out-of-plane lattice parameter was set to $c = 30$ Å, providing a vacuum spacing exceeding 20 Å to avoid interactions between periodic layers. Atomic positions were relaxed until the residual forces were below 0.03 eV/Å.

A $4 \times2 \times1$ supercell was used for the total-energy calculations. The noncollinear magnetic configurations generated by SpinGraph were imposed using constrained magnetic moments[54], and the resulting total energies were used for the subsequent spin-interaction fitting.

AMATIS calculations for monolayer CrI$_3$

We extract the spin interactions from the generalized spin Hamiltonian using AMATIS (Automated Magnetism Analysis via Tensorial Interacting Spins)[9, 55]. The symmetry constraints of P$\bar{3}$1m (No. 162) and the spin power limit are switched on. The truncated length of the rank-2 tensorial interaction, including spin-orbit coupling (SOC), is set to 12Å, up to $J_2^{5/2}$; we also consider the non-relativistic part of $J_2^{6/2}$. This setting ensures that as many interactions as possible are taken into account for a supercell of $4 \times2 \times1$. For rank-4 tensors, the two-ion interactions with both site multiplicities of 2 are taken into account, up to the second-nearest-neighbor $J_{2}^{4-2}$. We use 396 SpinGraph-generated symmetric magnetic configurations to calculate the total energies from first principles. The choice of the truncated lengths as well as the number of configurations for the energy fitting is already validated by the two-sided convergence tests[9], yielding the genuine systematic error of 0.10 meV per unit cell.

Data Availability

The data generated and analysed during this study are available from the corresponding author upon reasonable request.

Code Availability

The source code, manual, and tutorial package of SpinGraph are released in Zenodo (10.5281/zenodo.21947648). The archived code will be made available to the editors and reviewers during peer review upon request. A public GitHub repository will be released after submission, and the Zenodo record will serve as the persistent archived version of the code associated with this work.

Acknowledgements

H.L. acknowledges financial support from the National Key R&D Programs of China (No. 2025YFA1411000 and 2024YFA1410200), the Beijing Natural Science Foundation (No. 2232055), and the National Natural Science Foundation of China (No. 12204027). S.B. acknowledges financial support from the European Research Council (ERC) grant 856538 (project "3D MAGIC") and from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through SFB 1238 (Project C1).

Author Contributions

H.L. conceived the idea and supervised the project. Y.W. developed the theoretical framework, implemented the SpinGraph code, analysed the results, and wrote the manuscript. B.D. performed the first-principles calculations and contributed the corresponding computational data. S.B. reviewed the manuscript and provided critical feedback. J.R. and W.Z. provided scientific guidance and contributed to discussions. All authors discussed the results, reviewed the manuscript, and approved the submitted version.

Competing Interests

The authors declare no competing interests.