Spirit: Multifunctional framework for atomistic spin simulations

The Spirit framework is designed for atomic-scale spin simulations of magnetic systems with arbitrary geometry and magnetic structure, providing a graphical user interface with powerful visualizations and an easy-to-use scripting interface. An extended Heisenberg-type spin-lattice Hamiltonian including competing exchange interactions between neighbors at arbitrary distances, higher-order exchange, Dzyaloshinskii-Moriya and dipole-dipole interactions is used to describe the energetics of a system of classical spins localized at atom positions. A variety of common simulation methods are implemented including Monte Carlo and various time evolution algorithms based on the Landau-Lifshitz-Gilbert (LLG) equation of motion. These methods can be used to determine static ground-state and metastable spin configurations, sample equilibrium and finite-temperature thermodynamical properties of magnetic materials and nanostructures, or calculate dynamical trajectories including spin torques induced by stochastic temperature or electric current. Methods for finding the mechanism and rate of thermally assisted transitions include the geodesic nudged elastic band method, which can be applied when both initial and final states are specified, and the minimum mode-following method when only the initial state is given. The lifetimes of magnetic states and rates of transitions can be evaluated within the harmonic approximation of transition-state theory. The framework offers performant central processing unit (CPU) and graphics processing unit (GPU) parallelizations. All methods are verified and applications to several systems, such as vortices, domain walls, skyrmions, and bobbers are described.

Figure 1

The general structure of the framework, which is separated into a core library with an application programming interface (API) layer and a set of user interfaces (UIs). The core library handles input-output and calculations, while the API layer provides an abstract way of interacting with the code through several programming languages. The UIs provide direct control of calculations, as well as real-time visualization and postprocessing features. The back end for numerical calculations can be used in single-threaded and central processing unit (CPU)- as well as graphics processing unit (GPU)-parallel calculations.

Figure 2

Iterations per second of a LLG simulation over side length $L$ of a simple cubic system for 1 thread, for 10 threads, and on a GPU. The CPU parallelization consistently increases performance by almost an order of magnitude. By using a GPU, another order of magnitude can be gained for large system sizes, while the GPU performance at small system sizes is limited by the overhead of cuda kernel launches. Calculations were performed on a linux system with Intel Core i9-7900X 3.30 GHz and NVIDIA GeForce GTX 1080 processors.

Figure 3

Helicity (3) of a ferromagnetic cube, composed of $50\times 50\times 50$ spins on a simple cubic lattice with constant $a=1\AA$ and nearest-neighbor exchange of $J=16.86$ meV. The stray-field-induced helicity $\nu$ (circles) and ${\nu }^2$ (triangles) are shown in dependence on the reduced external magnetic field $h$. Both $\nu$ and $h=B/\left({\mu }_0{M}_{\mathrm{s}}\right)$ are unitless parameters, with the saturation magnetization density $M_{\mathrm{s}}$. The fitted curves (solid lines) show that the dependence of ${\nu }^2$ close to the critical field is approximately linear and they give a critical field value of $h_{\mathrm{c}}=0.159$, which matches the expected value of $h_{\mathrm{c}}=0.158$, as shown in Ref. [34], within $1\%$. The two insets show illustrations of how the cube will be magnetized at $h=0$ (left) and $h=0.1$ (right).

Figure 4

A $30\times 30\times 30$ ferromagnet with $J=1$ meV, with an expected critical temperature of $T_{\mathrm{C}}\approx 16.71$ K. Normalized values of the total magnetization $M$, susceptibility $\chi$, specific heat $C_{\mathrm{V}}$, and fourth-order Binder cumulant $U_4$ are shown. The magnetization is fitted with $M\left(T\right)={\left(1-T/T_{\mathrm{c}}\right)}^b$. At each temperature, ${10}^4$ thermalization steps were made before taking ${10}^5$ samples. Monte Carlo calculations give $T_{\mathrm{c}}\approx 16.60$ K, an agreement with expectation within $1\%$. The exponent is fitted with $b\approx 0.33$. The inset shows the Binder cumulants for system sizes of $L=30$, $L=20$, $L=15$, and $L=10$, giving an intersection at $T_{\mathrm{I}}=16.5\pm 0.25$, which is an excellent agreement with the expected value of $T_{\mathrm{c}}$ within the temperature step of $0.5\mathrm{K}$.

Figure 5

LLG calculation (top) and numerical error (bottom) of a single spin in an external magnetic field of $B=1$ T with a damping of $\alpha =0.1$ and a time step of $\mathrm{d}t=10$ fs, using the Depondt method. The error is the difference between the numerically calculated value and the analytical solution (11). Note that the error may depend strongly on the time step and damping. While the Heun method matches well with results shown in Ref. [9], giving an error within ${10}^{-6}$, the Depondt method shows a lower error of around $3\times {10}^{-7}$ with respect to the analytical solution.

Figure 6

The average velocity of head-to-head domain wall (see top) for various values of the nonadiabatic parameter $\beta$. For $\beta =0.10$, the Walker breakdown occurs at approximately $u_W\approx 0.01$. For $\beta =0$, a critical current is at $u_c\approx 0.0414$. From this point, the relation $\langle v\rangle =\sqrt{u^2-u_c^2}/(1+{\alpha }^2)$ mentioned by Thiaville et al. [53] takes effect. The mentioned relation is fitted to the data for $\beta =0$. For $\beta =0.1$ and currents under the Walker breakdown and $\beta =0.02$, the dashed lines show linear fits. Open symbols denote rotation around the $x$ axis. The results from Ref. [48] are reproduced well.

Figure 7

An illustration of the GNEB method for a single-spin system (the Hamiltonian and corresponding parameters are given in Appendix pp7). The two-dimensional energy landscape is shown superimposed on a unit sphere. The initial guess (green), relaxed path (blue), and final path using climbing and falling images (red) are shown.

Figure 8

The skyrmion tube (SkT) is either cut in half by the nucleation of a pair of Bloch points in the center (red minimum energy path) or separated from the upper surface by nucleation of a single Bloch point (blue minimum energy path). At a field strength of $H=0.8H_{\mathrm{D}}$, both processes have almost equal energy barriers of $\mathrm{\Delta }{E}_{\mathrm{center}}=23.13J$ and $\mathrm{\Delta }{E}_{\mathrm{surface}}=22.81J$. A chiral bobber is formed (two when the skyrmion tube is cut in half), whose collapse has an energy barrier of $\mathrm{\Delta }{E}_{\mathrm{bobber}}=7.55J$. Note that the slight differences in the collapse of the chiral bobber between the two paths come from different initial paths.

Figure 9

Energy barriers for the nucleation of Bloch points (BP) at the surface (blue circles) and in the center (green square), as well as the nucleation (red triangles up) and collapse (red triangles down) of a chiral bobber for a cube of size $30\times 30\times 30$ over applied magnetic field $H$. Periodic boundary conditions are applied in the $xy$ plane. The BP nucleation at the surface and center represents collapse of a skyrmion tube, while the bobber nucleation represents the creation of a BP in an otherwise homogeneous sample.

Figure 10

Lifetime $\tau =1/{\mathrm{\Gamma }}^{\mathrm{HTST}}$ of an isolated skyrmion in a periodic two-dimensional system, with $J=1$ meV and $D=0.6$ meV, as a function of temperature $T$ and external magnetic field $H$. The lifetime is given on a logarithmic scale with isolines ranging from 1 ps up to 1 year. Because only a single transition mechanism is taken into account, the structure of the graph is simple.

Figure 11

A single spin under the exchange and DMI interaction with another spin. The energy landscape is two dimensional and is projected onto a sphere. (a) The gradient force field, pointing away from the maximum and toward the minima. (b) The effective force field, pointing toward the saddle point. The resulting paths for four different starting points are shown (black, gray, and white lines). See Appendix pp9 for a visualization of the corresponding minimum mode directions.

Figure 12

Fragment of hexagonal lattice of magnetic spins, which illustrates the definition of the topological charge on a discrete lattice as given in the main text. $A_l$ is the area of a spherical triangle defined by vectors ${\mathit{n}}_i$, ${\mathit{n}}_j$, and ${\mathit{n}}_k$ located at the vertices of a triangle of lattice points (indicated shaded).

Figure 13

A $30\times 30\times 30$ ferromagnet with $J=1$ meV, with an expected critical temperature of $T_{\mathrm{C}}\approx 16.71$ K. The energy per spin $E$ and normalized values of the total magnetization $M$, susceptibility $\chi$, specific heat $c_{\mathrm{V}}$, and fourth-order Binder cumulant $U_4$ are shown. The value obtained from the simulation is $T_c\approx 16.92$ K, an agreement with the expectation of $1.2\%$. The exponent is fitted with $b\approx 0.33$. At each temperature, $2\times {10}^5$ thermalization steps were made before taking ${10}^6$ samples.

Figure 14

Schematic visualization of the projection of the tangents for a single-spin system. After a tangent ${\tau }_{\nu }^{\mathrm{FD}}$ is determined by finite difference calculation, it needs to be projected onto the tangent plane to the spin configuration so that it correctly points along the path. This tangent is denoted ${\tau }_{\nu }^{\mathrm{proj}}$ and can be calculated, e.g., by removing the component in the direction of the image; see Eq. (F1). Note that the tangent vector ${\tau }_{\nu }$ needs to be normalized, which for a multispin system needs to be performed in $3N$ dimensions.

Figure 15

Field of minimum eigenmodes (white lines) of a single spin in anisotropy and the interaction field of a second, pinned spin, corresponding to the force fields shown in Fig. 11. The minimum mode-following paths are shown in gray colors. The dashed lines show the separation of the convex regions around the minima from the rest of the configuration space.