All spin-dynamics simulations shown throughout the paper were performed with the Spirit code (Müller et al., 2019; Müller et al., 2021). The Spirit code provides a framework for atomic-scale spin simulations and combines both a graphical user interface as well as an easy accessible Python API. All simulations performed with Spirit are based on an extended Heisenberg Hamiltonian describing the interaction of spins S ⃗ i = M ⃗ i / μ i ( μ i = | M ⃗ i | ) sitting at lattice sites i. It can be written in its most general form as H = − ∑ ⟨ i j ⟩ J i j S ⃗ i ⋅ S ⃗ j − ∑ ⟨ i j ⟩ D ⃗ i j ⋅ S ⃗ i × S ⃗ j − ∑ i K ⊥ ( S ⃗ i ⋅ e ⃗ ̂ z ) 2 − ∑ i μ i B ⃗ ⋅ S ⃗ i − μ 0 8 π ∑ i ≠ j 3 S ⃗ i ⋅ r ⃗ ̂ i j S ⃗ j ⋅ r ⃗ ̂ i j − S ⃗ i ⋅ S ⃗ j r i j 3 − ∑ ⟨ i j k l ⟩ K i j k l S ⃗ i ⋅ S ⃗ j S ⃗ k ⋅ S ⃗ l . (1)

Here, the first line contains the isotropic and antisymmetric exchange interactions, the later also referred to as Dzyaloshinskii-Moriya interaction. The second and third line describe the on-site anisotropy, the Zeeman energy due to an external magnetic field B ⃗ , as well as the dipolar contribution. The last term allows to include higher-order exchange interactions (Hoffmann and Blügel, 2020) such as the conventional four-spin or four-spin-four-site interaction (Heinze et al., 2011), the four-spin-three-site interaction (Krönlein et al., 2018), as well as the biquadratic interaction (Szilva et al., 2013). The list of pairs ⟨ij⟩ and quadruplets ⟨ijkl⟩ as well as their respective parameters J _ij , D ⃗ i j , and K _ijkl can be defined by the user based on the desired use case. Furthermore, the system geometry such as the lattice symmetry and lattice size can be chosen arbitrarily and Spirit allows to introduce defects such as vacancies or atoms of different types. To obtain ground state as well as thermal properties of the investigated system, either the Monte Carlo method based on a Metropolis algorithm or Landau-Lifshitz-Gilbert dynamics can be used.

A more detailed description of the Spirit framework as well as its further functionalities, such as the possibility to calculate lifetimes of magnetic textures based on the combination of geodesic nudged elastic band and harmonic transition state theory calculations, can be found in Ref. (Müller et al., 2019).

AiiDA’s plugin system allows to combine various simulation codes and methods (to date more than 60 plugins exist already (The AiiDA team, 2021)) on the same footing while augmenting the calculation done through the AiiDA infrastructure with the stored data provenance. Albeit their significance in research on magnetic materials, spin-dynamics calculations have not been at the center of the AiiDA community so far. To the best of our knowledge, besides the AiiDA-Spirit plugin presented here, only a first version of the AiiDA-UppASD plugin (Xu et al., 2021) exists for the UppASD code (Skubic et al., 2008) to combine AiiDA with a spin dynamics simulation engine.

In the context of AiiDA, a calculation plugin needs to be able to generate typical input files that are required to run a calculation through a bash script that will be generated when a calculation is submitted to a computer or as a job on a supercomputer. At the heart of the AiiDA-Spirit plugin lies the SpiritCalculation that connects the Spirit code via the Spirit Python API to AiiDA. The Layout of the SpiritCalculation is shown in Figure 1. To run a Spirit calculation a number of input nodes are required:

• a structure node describing the lattice of spins (i.e., their positions in the unit cell),

Additionally, input modes that trigger special features of the Spirit code such as disorder and defects in the structure or pinning of spins to certain directions can be controlled with the corresponding optional input nodes. The SpiritCalculation then implements the functionality to translate this information into the appropriate input files and runs the calculation using the Spirit Python API. The AiiDA daemon automatically takes care of creating a suitable job script, copying necessary input files, and of submitting and monitoring the calculation run. Once the calculation job finishes, important output files are copied back to the retrieved folder in the AiiDA file repository associated to the AiiDA database. Then, the SpiritParser extracts useful information that should be stored in the database. For the example of a LLG calculation this entails settings such as the number of LLG steps until convergence, the used wall clock time on the computer where the calcualtion ran, an array of the energies (i.e. exchange energy per spin), and the initial and final directions of the spins in the magnetization array.

Apart from the SpiritCalculation and SpiritParser, AiiDA-Spirit comes with some tools that can be used in the typical jupyter notebook environment that is often used in the context of AiiDA. In particular we mention the show_spins tool of AiiDA-Spirit which provides the Spirit visualization capabilities in a simple Python API. This consists of a WebAssembly and WebGL version of the VFRendering library (Vfrendering, 2021) in combination with a JavaScript interface that can be used to visualize the directions of the spins from the web-browser based environment natural to jupyter notebooks.

The density functional theory (DFT) results of this work were produced within the generalized gradient approximation (GGA-PBE) (Perdew et al., 1996) using the full-potential scalar-relativistic Korringa-Kohn-Rostoker Green’s function method (KKR) (Ebert et al., 2011) as implemented in the JuKKR code package (The JuKKR developers, 2021). We use an ℓ _max = 3 cutoff in the angular momentum expansion with an exact description of the atomic cells (Stefanou et al., 1990; Stefanou and Zeller, 1991). After the self-consistent DFT calculations, the method of infinitesimal rotations (Liechtenstein et al., 1987) was used to compute the exchange interaction parameters J _ij . The series of DFT calculations in this study are orchestrated using the AiiDA-KKR (Rüßmann et al., 2021a) plugins to the AiiDA infrastructure (Huber et al., 2020). The complete dataset that includes the full provenance of the calculations is made publicly available in the materials cloud repository (Talirz et al., 2020; Rüßmann et al., 2021).

To illustrate the usage of AiiDA-Spirit we first consider a toy model consisting of a single layer of spins in a simple-cubic lattice. The complete example is part of the dataset that accompanies this publication (Rüßmann et al., 2021). We assume only nearest neighbor interactions with isotropic exchange J ₁ = 10 meV and Dzyaloshinskii-Moriya interactions with a strength of D ₁ = 6 meV. This choice of parameters does not reflect any concrete physical system but is chosen for illustration purposes because it is known to produce skyrmions with small radii. The generation of the corresponding input node for the SpiritCalculation where including the directions of the DMI vectors can be seen in the following code snippet.

Here i, j index the lattice site in the unit cell (situated at r ⃗ i , r ⃗ j ) and da, db, dc describe offsets into unit cells further away such that R ⃗ i j = r ⃗ j − r ⃗ i + d a a ⃗ + d b b ⃗ + d c c ⃗ describes the distance between two spins ( a ⃗ , b ⃗ , c ⃗ are the Bravais vectors of the lattice), Jij and Dx, Dy, Dz denote the exchange interaction and the three components of the DMI vector. For this example we consider an isolated layer of spins with periodic boundary conditions in the plane. We choose a supercell for the SpiritCalculation of 50 × 50 × 1 spins. Furthermore we apply an external field of various strength (in the code snippet we show the input parameters for a value of 25 T) in the direction perpendicular to the film in the following code snippet.

Starting from random orientations of the spins we then perform a time evolution using the LLG method with the Depondt solver (Depondt and Mertens, 2009). The parameters for the LLG calculations are summarized in Table 1.

To harness the high-throughput capabilities of the AiiDA-Spirit plugin we perform a series of SpiritCalculations to screen a range of external fields and temperatures. We change the temperature from 0 to 75 K in 2.5 K steps and vary the external field from − 50 T to + 50 T in steps of 2.5 T. The calculations for each parameters set are repeated 5 times starting from different random orientations of the spins for statistical averaging. This amounts to 31 × 41 × 5 = 6,355 individual SpiritCalculations that were submitted to an in-house compute cluster. We stress that the AiiDA daemon (Uhrin et al., 2021) conveniently takes care of creating submission scripts and automatically retrieves and parses the outcome of the calculations without the need for any user interaction. A visualization of the dataset and the provenance graph for this application is shown in Figure 2.

In order to analyze the outcome of the SpiritCalculations we chose to investigate the topological charge in the simulation cell at the end of the LLG simulation. For a continuous vector field m ⃗ it is defined as ρ T = 1 4 π ∫ m ⃗ ⋅ ∂ m ⃗ ∂ x × ∂ m ⃗ ∂ y d x d y . (2)

We added a custom post-processing step to the SpiritCalculation which uses the get_topological_charge function of the spirit Python API. This function calculates the topological charge from the discretized form of Eq. 2 as a summation over all contributions of triangles formed by neighboring spins in the simulation cell (Müller et al., 2019).

Figure 3 shows the outcome of these simulations where the topological charge is shown for all 1,271 pairs (T, B _z ) together with selected spin configurations of representative calculations marked by the symbols (Figures 3B–G). The real-space spin configuration at the end of the LLG calculations were visualized using the show_spins tool of the AiiDA-Spirit plugin. It can be seen that a small external field leads to the appearance of skyrmions which in this case have a topological charge of ± 1, depending whether they form in a ferromagnetic background of spins pointing in − z (Figure 3C) or + z (Figure 3E) direction. In general, the topological charge counts the difference between the amount of skyrmions in up-domains (ρ _T > 0) and skyrmions in down-domains (ρ _T < 0) as seen for vanishing external field in (Figure 3D) where several skyrmions with opposite topological charges lead to a near cancellation of the total topological charge.

At very large fields, the Zeeman exchange coupling term becomes larger than the DMI energy and a homogeneous ferromagnet forms (Figures 3B,G). Temperature fluctuations tend to deform the skyrmions (Figure 3F) and destabilize them. Thus, at elevated temperatures smaller magnitudes of the external magnetic field leads to a vanishing topological charge. As highlighted by the black arrow in (Figure 3A), this is however only true up to a certain critical temperature. For T > 40 K the topological charge increases again which can be explained by the energy barrier of skyrmion formation and destruction. At these elevated temperatures the fluctuations of the spin directions are larger than 40 K ⋅ k _B ≈ 3.45 meV which we conjecture is the energy barrier for skyrmion formation. While the energy barrier can in principle be calculated by performing geodesic nudged elastic band calculations, this is beyond the scope of this paper and therefore will be omitted. The larger temperature fluctuations also prohibit reaching the force convergence criterion set in the LLG calculation which means that the LLG simulation runs until the maximal simulation time of 100 ps is reached. During this simulation time, skyrmions can spontaneously form and disappear which results in a finite topological charge measured at the end of the run. In the future the real time dynamics of skyrmion creation and collapse may be the focus of the investigation. However, this approach may become unfeasible for situations where the skyrmion lifetime is very long compared to the typical time step in LLG calculations. Finally, we highlight that with increasing temperature fluctuations we also find a larger variance in the number of skyrmions when averaging over the five different starting configuration for each pair (T, B _z ). This supports our interpretation that skyrmions are spontaneously created and annihilated by temperature fluctuations.

The Monte Carlo (MC) method is a well established tool in physics which, when applied to spin systems, allows to estimate the critical temperature of the magnetic ordering (Curie temperature) (Binder and Heermann, 1997). The Spirit code (Müller et al., 2019) implements a Metropolis algorithm which can be used from AiiDA-Spirit by choosing the mc simulation method (instead of the previously used LLG method). We demonstrate the MC at the example of a simple-cubic ferromagnet with only nearest neighbor interactions J ₁ = 1 meV. We perform calculations for varying supercell sizes between 10 × 10 × 10 and 40 × 40 × 40 with the MC parameters given in Table 2. The results of the calculation are shown in Figure 4 where, together with the total magnetization M, the isothermal susceptibility χ = 1 k B T ( 〈 M 2 〉 − 〈 M 〉 2 ) (3) with M = | 1 N ∑ i S ⃗ i | the average magnetization of the sample is shown. We see that M(T) converges with increasing supercell size indicating that boundary effects become less important. The corresponding susceptibilities show a diverging behavior at T _c . Our calculation results agree well with the expected value of T _c = 1.44J ₁/k _B = 16.71 K. We stress that these calculations require a series of steps consisting of, for example, thermalization and decorrelation steps at each temperature value in the scanning interval. Within AiiDA-Spirit this complexity is conveniently absorbed in the SpiritCalculation which greatly facilitates the application of MC calculations.

We now demonstrate how the integration of the Spirit code into the AiiDA framework through the AiiDA-Spirit plugin can facilitate multi-scale modeling for magnetic materials. In this example we first calculate the exchange interaction parameters for γ–Fe using density functional theory which are then passed to the AiiDA-Spirit plugin for LLG simulations.

The γ phase of Fe is a metastable high-temperature phase where the atoms crystallize in the fcc lattice (Knöpfle et al., 2000; Sjöstedt and Nordström, 2002). This has a drastic consequence on the exchange interactions where, in contrast to the ferromagnetic bcc Fe, frustrated exchange interactions can lead to the formation of spin-spirals. Experimentally this structure of Fe can be realized in a Cu matrix (Tsunoda, 1989; Tsunoda et al., 1993). It is known that a variation of the lattice constant of γ–Fe can have drastic consequences for the magnetic ordering (Sjöstedt and Nordström, 2002). Here, we investigate bulk crystals of γ–Fe for varying lattice constants between a lat = 3.2 A ̊ and a lat = 4.0 A ̊ around the lattice constant of Cu ( a lat = 3.6 A ̊ ) .

Figure 5 summarizes the results of the DFT calculations that were done with the AiiDA-KKR plugin (see methods section for numerical details). The total energy as a function of the lattice constant (shown in panel Figure 5A) reveals a phase transition from the low-spin state (for a lat < 3.6 A ̊ ) to the high-spin state ( a lat ≥ 3.6 A ̊ ) of γ–Fe as seen in the jump of the spin moment to μ _s > 2.5 μ _B . This coincides with a smaller exchange splitting seen in the density of states (Figure 5B) and consequently a smaller value of the spin moment (Figure 5C). For lattice constants below 3.37 A ̊ we find that the magnetic moment vanishes. Panel (Figure 5D) shows the calculated exchange interactions J _ij as a function of distance between two Fe atoms. Clearly, the sign of the nearest neighbor interaction shows the most drastic change with the transition from high-spin to low-spin state at smaller lattice constant of γ–Fe. While in the high-spin state the first and second nearest neighbor interaction are both ferromagnetic (J _ij > 0), for the low-spin state the nearest neighbor interaction changes from being weakly ferromagnetic to antiferromagnetic (J _ij < 0).

In the following, the consequences of this change for the magnetic ordering are investigated based on a series of LLG calculations using the AiiDA-Spirit plugin. In the DFT calculation we use the primitive cell which contains a single atom in the unit cell. For the spirit calculations we map the calculated exchange interactions onto the conventional unit cell consisting of four atoms. The parameters of the LLG simulations are summarized in Table 3. We study the magnetic ordering in a 40 × 40 × 40 × 4 = 256,000 spins supercell without external magnetic fields and at temperature T = 0 K. Here we focus on the ground state that forms and therefore neglect effects of temperature fluctuations and external fields which can, for example if T is high enough, overcome the energy barrier between different (metastable) magnetic orderings. We further neglect the influence of anisotropy (K _⊥ = 0) in this work and we also do not include higher order exchange terms (K _i,j,k,l = 0). We choose open boundary conditions in order to not bias the eventually resulting spin-spiral wavelength by the periodicity of the supercell. Table 4 summarizes the DFT calculated values for magnetic moments and exchange coupling constants for varying lattice constants that were used in the respective SpiritCalculations. Note that the exchange coupling constants are only shown up to the seventh shell of neighbors but the calculations included pairs up to the 15th shell that are however much smaller than the values reported in Table 4.

We start the discussion of the LLG calculations with the results for γ–Fe in the lattice constant of Cu ( a lat = 3.61 A ̊ ) . Figure 6A shows the resulting spin texture at the end of the LLG calculation for the central layer of spins in the yz-plane (shown in the inset Figure 6B). We can see that a spin-spiral forms in z-direction with ferromagnetically ordered spins in y-direction. At the open boundaries we see that the missing neighbors on one side influence the magnetic ordering that deviates from the spin-spiral in the center for a distance of about five lattice constants. In order to quantify the spin-spiral we pick the two cardinal directions in this plane (indicated by blue and orange lines in Figure 5A) and extract the z-component of the spin S _z . We combine the projections onto the yz-plane from two adjacent layers of spins (indicated by the two grey planes in Figure 5B) to not restrict our analysis to a single sub-lattice only. This allows to describe also antiferromagnetic structures in the sub-lattices with ferromagnetic ordering within one sub-lattice, which will be important later. Figure 6C shows that, except for boundary effects, S _z stays constants when following the y-direction. Along the z-axis we see a complex beating pattern with the site index that can be decomposed into two π/2-shifted spin-spirals in the two different sub-lattices. The corresponding Fourier transformation F j ( q ) = 1 2 π ∫ e − i q r j S z ( r j ) d r j (4) with j = y, z computed with the fast Fourier transform algorithm (FFT) is shown in (Figure 6D). As expected, the FFT of the predominantly ferromagnetically ordered spins along the y-direction | F y | mainly shows a signal at q = 0 whereas | F z | shows the appearance of four peaks at q ≈ ± 0.2 2 π a lat and q ≈ ± 0.8 2 π a lat = ( 1 − 0.2 ) 2 π / a lat which are attributed to the two π/2-shifted oscillations in the two sub-lattices. We point out that the reflection symmetry around q = 0 is a consequence of the real-valued input to the FFT and is therefore not discussed further. As seen from the FFT in y-direction and from the corresponding spin texture in (Figure 6A) the spins in the direction perpendicular to the propagation direction of the spin spiral (i.e., the z-direction in this example) are ordered ferromagnetically. Therefore the spin-spiral wavevector is q ⃗ = ( 0,0,0.2 ) 2 π / a lat . After having characterized the spin-spiral ground state of γ–Fe we continue with a discussion of the magnetic ordering depending on the changing exchange coupling parameters with changing lattice constant.

Figure 7 summarizes the LLG calculations for γ–Fe for varying lattice constants using the respective set of exchange parameters shown in Figure 5D. The lines in Figure 7A show the energy at the end of the LLG calculation starting either from a random spin configuration (E, dashed orange) or from the ferromagnetic (E _FM, solid blue) state. We find that for lattice constants a lat ≥ 3.65 A ̊ the ferromagnetic state minimizes the energy (E − E _FM = 0). We attribute this to the increasing ferromagnetic interaction for nearest neighbor spins in the high-spin state which was discussed with Figure 5D. At smaller lattice constants ( a lat ≤ 3.6 A ̊ ) , the ferromagnet (FM) is not the ground state anymore. Here we find either a spin-spiral (SS) ground state or an antiferromagnetic (AFM) phase. Figure 7C shows three representative images of the SS, AFM and FM states.

In the SS state the magnetization rotates from left to right (i.e. along the y-axis) and shows antiparallel alignment of the rows in z-direction. Along the x-axis (direction perpendicular to the drawn plane) the spins are aligned ferromagnetically, except for boundary effects at the open ends of the simulation cell (seen in the direction of the first layer of spins). In z-direction, adjacent layers are antiferromagnetically ordered. Thus the spin-spiral wavevector for these lattice constants has the form q ⃗ = ( 0 , q , 1 ) 2 π / a lat which due to the cubic symmetry of the crystal is equivalent to q ⃗ = ( q , 0,1 ) 2 π / a lat . Note that q ⃗ AF = ( 0,0,1 ) 2 π / a lat is the antiferromagnet because the distance between two layers in the (0, 0, 1) direction of the fcc lattice is a _lat/2.

In the AFM phase the direction of the spins separate into four sub-lattices that correspond to the four atoms in the conventional fcc unit cell. Within each sub-lattice the spins are aligned parallel and form a right angle with their neighboring spins from different sub-lattices. This is highlighted with a red box in the middle panel of (Figure 7C). In the FM phase (lower panel) all spins point in the same direction. Note that in all these calculations the spins can collectively rotate since we neglected contributions from single-ion anisotropies and do not apply an external field.

The summed magnitude of the Fourier transform along the three cardinal axes | F ( q ) | = ∑ j = x , y , z | F j ( q ) | (5) is shown in Figure 7B for different lattice constants. Note that we have summed here over the symmetry-equivalent directions along the x- y- and z-directions because of the rotational invariance of the complete spin-structure. Starting from the smallest lattice constant of a lat = 3.39 A ̊ we see a peak in | F ( q ) | at q ⃗ = ( ± 0.195 , 0,1 ) 2 π / a lat which corresponds to a wavelength of the spin-spiral of 5.13 a _lat. Here we focus our discussion on the peak at smaller q values as the same arguments hold for the second peak at 2π/a _lat − q as discussed above. With increasing lattice constant the spin-spiral wavelength increases to 7.82 a _lat/2 ( q ⃗ = ( ± 0.128 , 0,1 ) 2 π / a lat ) at a lattice constant of a lat = 3.5 A ̊ . For 3.5 A ̊ < a lat < 3.6 A ̊ the AFM state is found, which in the Fourier transform is characterized by the dominating peak at q = ±2π/a _lat. Note that we still get a considerable signal at q = 0 because we sum over all three cardinal directions and there is ferromagnetic ordering along one direction (see Figure 7C). As discussed above, for larger lattice constants the spin-spiral state briefly shows up again at a lat = 3.61 A ̊ with a wavelength of 5.59 a _lat which however has ferromagnetically ordered spins in both directions perpendicular to the direction of spin-spiral propagation ( q ⃗ = ( ± 0.179 , 0,0 ) 2 π / a lat ) until for a lat ≥ 3.65 A ̊ the FM state is found which only shows a significant Fourier amplitude at q = 0.

The appearance of the AFM phase for 3.5 A ̊ < a lat < 3.6 A ̊ can be attributed to the sign change of the nearest neighbor interaction from ferro- to antiferromagnetic. To verify this hypothesis we employ a series of LLG calculations through AiiDA-Spirit. We chose to start from the setup of the calculation for a lat = 3.61 A ̊ , which was found to reproduce the spin-spiral phase. We then modify the strength of the nearest neighbor interaction J ₁ ranging from − 5 meV to + 15 meV and run LLG calculations starting from the FM state, from the AFM state, from the SS phase and random spin orientations. For the AFM state we construct the row-wise AFM orientation of the spins. All other parameters for the LLG calculation are kept constant. In total this is another set of 84 SpiritCalculations where we find that starting from a random spin configuration coincides with starting from the SS. The random starting point is therefore omitted in the following discussion. Figures 8A,B show the dependence of the energy per spin at the end of the LLG calculations. We point out that the LLG calculations for the FM and AFM states are converged in the very first iteration which indicates that the FM and AFM phases are local minima in the energy landscape. We find that the SS state is lowest in energy with a maximal energy gain of ∼ 7 meV/spin in the transition region where E _FM − E _AFM changes sign. From the final spin structure of the spin-spiral solution we proceed with an analysis of the Fourier components as introduced in equations 4 and 5. This is shown in Figure 8C. As highlighted by the grey line, we see an increase in the spin-spiral wavevector with increasing J ₁ up to the point where E _FM and E _AFM cross around J ₁ = 7.3 meV. As in the previous analysis for changing lattice constant we find that the spin-spiral state is characterized by two wavevectors at q and 2π/a _lat − q (grey dashed line). Furthermore, the Fourier transform for all states show significant signals at q = 0 (indicating parallel spins) and q = ±2π/a _lat (indicating antiparallel spins).

Overall we can conclude that the resulting spin-texture in the 256,000 spin unit cell with open boundary conditions is a result of the complex competition of distance-dependent exchange couplings that favor ferromagnetic or antiferromagnetic alignments of spins or can compete and stabilize spin-spiral ground states.