Uncovering multi-channel magnetic hopfion annihilation via a single-node, billion-spin-scale atomistic framework Qichen Xu*1, 2 and Anna Delin1, 2, 3 1
arXiv:2605.26016v1 [physics.comp-ph] 25 May 2026
Department of Applied Physics, School of Engineering Sciences, KTH Royal Institute of Technology, AlbaNova University Center, SE-10691 Stockholm, Sweden 2 Swedish e-Science Research Center, KTH Royal Institute of Technology, SE-10044 Stockholm, Sweden∗ 3 Wallenberg Initiative Materials Science for Sustainability (WISE), KTH Royal Institute of Technology, SE-10044 Stockholm, Sweden (Dated: May 26, 2026) Modern atomistic spin simulations combine long stochastic trajectories, thermodynamic sampling, static optimization and multi-image transition-path workflows, all of which rely on repeated evaluation of spin Hamiltonians and become computationally prohibitive on the large lattices required for three-dimensional magnetic textures. We introduce SpinX, a GPU-native atomistic spin simulation framework built around a unified Hamiltonian interface and multiple user-selectable computational backends. Its core is a crystallographic sublattice decomposition that reformulates translationally invariant spin interactions as multi-channel tensor convolutions, enabling dense, sparse and FFT-based convolution backends, while irregular systems are handled by pair-list evaluation and long-range dipolar fields by reciprocal-space FFT. Implemented in JAX, SpinX supports deterministic and stochastic Landau–Lifshitz–Gilbert dynamics, Monte Carlo sampling, static optimization, dynamical spectroscopy and string and geodesic nudged elastic band transition-path calculations on heterogeneous accelerator platforms. A validated mixed-precision mode combines fp32 field evaluation with fp64 spin-state propagation. We validate SpinX against analytical single-spin dynamics, finite-size thermodynamics of bcc Fe and transverse dynamic structure factors. Performance benchmarks show peak throughput exceeding 10 billion spin-site operations per second on a single accelerator and aggregate single-node workloads of over 1 billion atomic spins. Applying this framework to an exchange-stabilized magnetic hopfion, we uncover two competing annihilation channels on a million-spin atomistic lattice: a previously reported axial-collapse pathway and a distinct lateral-rupture pathway with a different transition morphology and activation barrier. By combining accelerator-native throughput with large-scale transition-state workflows, SpinX establishes atomistic spin simulation as a practical route to studying three-dimensional magnetic textures and their energy landscapes.
MAIN
Many phenomena in modern magnetism — including lattice-scale frustrated exchange, non-collinear spin textures, thermal fluctuations and topological singularities such as Bloch points — couple atomic-scale interactions to mesoscale magnetic textures [1–3]. Characterizing the stability and energy landscapes of such textures requires integrated simulation workflows that combine long stochastic trajectories, thermodynamic sampling, static minimization and multi-image transition-path searches [4–8]. This simultaneous demand for large spatial volumes, atomistic lattice resolution and repeated path evaluations defines a computational regime that is not well served neither by small atomistic cells nor by continuum descriptions that smooth over lattice-scale singular events [9]. Three-dimensional topological solitons such as magnetic hopfions are a particularly demanding example [10–13]. Their stability is determined not only by the relaxed metastable texture, but also by the topology and energetics of the annihilation pathways connecting them to the trivial ferromagnetic state [4, 5]. Resolving
such pathways in practice can require million-spin-scale lattices to separate the texture from spurious boundary effects while preserving atomistic resolution [4, 5]. Several mature simulation packages have shaped the current landscape of atomistic and continuum spin modelling. UppASD and VAMPIRE provide established platforms for atomistic spin dynamics, thermodynamic sampling and finite-temperature magnetic phenomena, while Spirit integrates atomistic dynamics with Monte Carlo sampling, transition-path methods and visualization [14– 17]. On the continuum side, mumax3 has demonstrated the transformative impact of GPU acceleration for highthroughput micromagnetic simulations. These developments have made many classes of spin simulations routine. However, emerging workflows for three-dimensional magnetic textures increasingly require a different combination of capabilities, including atomistic lattice resolution, accelerator-native execution across heterogeneous GPU architectures, multiple mathematically equivalent Hamiltonian backends, a unified treatment of regular and irregular interaction graphs, automatic-differentiation compatibility and multi-image transition-path calculations at the million-spin-per-image scale. Existing atomistic codes typically rely on solver-specific data structures and explicit interaction lists, which are general but do not naturally expose regular crystal Hamiltonians as
2
FIG. 1. Modular architecture and backend dispatch in SpinX. (a) Overview of the SpinX software stack. User inputs, including TOML configuration files, tabulated Hamiltonian parameters, and spin-state files, are parsed and checked in the simulation setup layer. The physics functionality is organized as decoupled Python modules for spin dynamics, thermodynamic sampling, static optimization, analysis, minimum-energy-path calculations, and data storage. These modules call a unified Hamiltonian backend, which dispatches effective-field evaluations to specialized kernels, including sublattice dense convolution, sublattice sparse convolution, FFT-based convolution, pairlist evaluation, and dipolar FFT. JAX/XLA compilation provides hardwareportable execution on CPUs, NVIDIA GPUs, and AMD GPUs, while trajectories, checkpoints, path data, and computed observables are stored in HDF5 format. (b) Backend selection for effective-field evaluation. For ideal crystals with translational symmetry, SpinX decomposes the magnetic lattice into sublattice channels (analogous to colour channels in image processing) and maps pairwise bilinear interactions onto dense, sparse, or FFT-based tensor convolutions. The FFT backend precomputes the interaction kernel in reciprocal space and evaluates the field by multiplying the spin spectrum followed by an inverse FFT. Long-range dipolar interactions use the same reciprocal-space convolution principle through a dedicated dipole-FFT backend. For non-ideal crystals, including chemical disorder, structural defects, and explicit site-to-site connectivity, SpinX uses a generalized pairlist representation that stores local Hamiltonian parameters such as the exchange Jij and Dzyaloshinskii– Moriya vectors Dij .
convolutional or FFT-convolutional workloads. A framework that combines these capabilities within a single compiler-driven Hamiltonian engine would make largescale atomistic workflows more composable, portable and directly compatible with differentiable programming and machine-learning interfaces. Here we introduce SpinX, a GPU-native atomistic spin simulation framework designed for high-throughput and large-scale calculations on modern accelerators. The central design principle of SpinX is to express spin Hamiltonians through a unified interface that covers scalar and tensorial bilinear interactions as well as higherorder terms such as biquadratic, scalar-chirality and quadruplet interactions, and to evaluate the corresponding effective fields through user-selectable computational
backends. For ideal crystals with translational symmetry, crystallographic sublattice decomposition reformulates translationally invariant pair interactions as multichannel tensor convolutions: each crystallographic sublattice plays the role of a colour channel in image processing, and interactions between sublattices are encoded as cross-channel kernels. This mapping enables dense, sparse and FFT-based convolution backends to be implemented using Accelerated Linear Algebra (XLA) convolution, FFT and fused array primitives, while preserving a common Hamiltonian interface[18]. For non-ideal crystals with chemical disorder, structural defects or explicit site-resolved couplings, the same Hamiltonian representation is evaluated through a generalized pair-list backend. Long-range dipolar interactions are treated by
3
FIG. 2. Numerical validation of SpinX using analytical, thermodynamic, and dynamical benchmarks. (a) Single-spin deterministic LLG benchmark. The trajectory root-mean-square error is computed over a 500 ps trajectory sampled every 1 ps by comparing Heun integration with the analytical precession solution for time steps ∆t = 0.1-100 fs. In mixed-precision mode, effective-field evaluations are performed in single precision, while the spin state is stored and advanced in double precision. This mixed-precision scheme closely follows the fp64 baseline over the tested range, whereas pure fp32 exhibits a round-off error floor at small time steps. (b) Thermodynamic validation by Monte Carlo sampling for a classical bcc Fe Heisenberg model. The magnetization M , Binder cumulant U4 , magnetic susceptibility χ, and specific heat Cv are shown as functions of temperature, M , χ, and Cv are scaled for visual comparison. The shaded region marks the estimated critical-temperature region. The inset shows Binder-cumulant crossings for system sizes L = 12-36, used for finite-size estimation of the ferromagnetic-toparamagnetic transition temperature. (c) Dynamical validation from the transverse dynamic structure factor. The spectrum S(q, ω) = Sxx (q, ω) + Syy (q, ω) is evaluated along the high-symmetry Γ-X-M -Γ path, with intensity shown in arbitrary units on a logarithmic scale. The white dashed curve denotes the adiabatic magnon spectrum (AMS) calculated from linear spin-wave theory, which follows the dominant simulated spectral weight across the Brillouin-zone path.
a reciprocal-space FFT convolution backend. This backend structure allows SpinX to exploit lattice regularity when available, while retaining a general representation for irregular lattice settings. SpinX is implemented in JAX and compiled through XLA, enabling hardware-portable execution across heterogeneous accelerator platforms, including both NVIDIA and AMD GPUs[19]. To our knowledge, SpinX is the first atomistic spin-dynamics framework to systematically employ a mixed-precision execution mode in which effective-field evaluation is performed in fp32 while spin-state propagation is retained in fp64. This separates the most bandwidth-intensive part of the calculation from the accuracy-sensitive integration state, reducing memory traffic and exploiting the lower-precision arithmetic that modern accelerators natively accelerate, while preserving long-time integration accuracy. The framework follows a modular functional design: deterministic and stochastic Landau–Lifshitz–Gilbert (LLG) dynamics, heat-bath and Metropolis Monte Carlo sampling, static optimization, spectral analysis and minimumenergy-path calculations all call the same shared Hamiltonian engine. This avoids duplicating physics logic across solvers and allows independent modules and backends to be cross-validated against one another. Because the framework is written in JAX, registered Hamiltonian terms compose naturally with automatic differentiation, batching and just-in-time compilation, providing a direct interface to gradient-based optimization, sensitivity analysis and machine-learning workflows. We validate SpinX using analytical, thermodynamic and dynamical benchmarks. Single-spin deterministic dynamics are benchmarked against the analytical LLG
precession solution across multiple time-step and precision settings [1, 2]. Thermodynamic sampling is validated through heat-bath Monte Carlo simulations of a classical body-centred cubic (bcc) Fe Heisenberg model, including a Binder cumulant finite-size scaling analysis of the ferromagnetic transition [20]. Dynamical correlations are tested by computing the transverse dynamic structure factor and comparing the dominant simulated spectral weight with the adiabatic magnon spectrum obtained from linear spin-wave theory [20]. Together, these tests verify both the physical accuracy of the solvers and the consistency of the shared Hamiltonian engine. Performance benchmarks show that SpinX achieves peak throughput exceeding 10 billion spin-site operations per second on a single accelerator and supports aggregate single-node workloads of over 1 billion atomic spins. This scale is particularly vital for path-based calculations, where a multi-image path-ensemble multiplies the effective number of spin degrees of freedom. Rather than exposing a single universally fastest kernel, the framework provides a set of mathematically equivalent backends with distinct performance profiles. Dense convolution, sparse convolution, FFT convolution, pair-list evaluation and dipolar FFT each dominate in different regions of interaction range, lattice regularity, floating-point precision and memory footprint. This regime-based architecture is essential because different physical workloads–ranging from short-range dynamics and extended frustrated exchange to long-range dipolar interactions, thermal sampling and static minimization–each present unique computational bottlenecks. As a demonstration of the scientific capability enabled by this framework, we investigate annihilation pathways
4
FIG. 3. GPU performance and scaling benchmarks of SpinX. Benchmarks were performed on NVIDIA GH200 GPUs (Dardel, KTH, GPU portion only) and AMD MI250X accelerators (LUMI, CSC). For MI250X, throughput is reported per Graphics Compute Die (GCD), each physical card contains two GCDs. Unless otherwise stated, timings exclude JIT compilation and include device synchronization. Error bars denote the standard deviation over repeated executions.(a) Single-device backend throughput for a nearest-neighbour exchange model. Benchmarks are executed on a 200 × 200 × 200 spin system. Dense convolution, sparse convolution, FFT convolution, and pairlist backends are compared on NVIDIA and AMD hardware in fp64, mixed-precision, and fp32 modes. Solid bars show Heun integration throughput, whereas hatched bars show isolated effectivefield evaluation throughput, the two quantities should be interpreted independently rather than stacked. Missing bars indicate unsupported backend/precision combinations.(b) Dependence of throughput on Hamiltonian complexity. Left, mixed-precision throughput as the exchange-interaction shell count increases from one to five. Right, throughput for nearest-neighbour models containing exchange only (J), exchange plus Dzyaloshinskii–Moriya interaction (J +D), and exchange plus DMI plus long-range dipolar interactions (J + D + dipolar). Solid and hatched bars denote Heun integration and isolated effective-field evaluation, respectively.(c) Single-node multi-GPU scaling. Spatial domain decomposition is benchmarked using up to four devices for a 204 × 200 × 200 nearest-neighbour model, with the first dimension chosen to be divisible by the number of devices. Ensemble scaling distributes independent trajectories across devices. Translucent background bars, read against the right axis, indicate the maximum number of spins that fit in memory for the mixed-precision dense-convolution workflow.(d) Monte Carlo sampling throughput as a function of exchange shell count. Heat-bath and Metropolis updates are compared on NVIDIA and AMD hardware, with throughput reported as billion attempted spin updates per second.(e) Static-minimization throughput as a function of exchange shell count. Velocity projection (VP), Barzilai–Borwein (BB), and L-BFGS are benchmarked over a fixed number of iterations, with throughput reported as billion spin-site updates per second.
of an exchange-stabilized magnetic hopfion in a minimal frustrated-exchange model [4, 5]. The goal is not to provide material-specific predictions, but to resolve competing topological collapse mechanisms in a controlled atomistic setting. Previous transition-state studies of exchange-stabilized hopfions have focused on an axialcollapse mechanism, in which Bloch-point-like defects nu-
cleate near the hopfion axis and mediate collapse of the toroidal texture [4, 5]. Using string-method initialization followed by geodesic nudged elastic band refinement on a (100 × 100 × 100) atomistic lattice, corresponding to (106 ) spins per image, under open boundary conditions, we identify a second annihilation channel in the same Hamiltonian: a lateral-rupture pathway. In this
5
FIG. 4. Two competing annihilation pathways of magnetic hopfions. (a) Refined transition paths connecting the initial hopfion state to the ferromagnetic (FM) state. Both pathways are obtained by string-method initialization followed by GNEB refinement on an atomistic lattice, with open boundary conditions and the hopfion initialized near the center of the simulation box, the large surrounding ferromagnetic region reduces spurious boundary effects on the annihilation mechanism. The blue path corresponds to a lateral-rupture channel, whereas the red path corresponds to an axial-collapse channel. Reaction coordinates are normalized independently for each path. The lateral-rupture saddle region was further refined using a denser image distribution and stringent CI-GNEB relaxation to improve the resolution of the local rupture process. Energies are shown in units of the nearest-neighbour exchange constant J1 . Open stars mark the highest-energy images along the two paths. Insets show magnified views of the barrier regions, with activation barriers measured relative to the local metastable minimum on the hopfion side. The lower inset highlights the saddle-region energy offset between the two channels. (b) Representative configurations along the lateral-rupture path. The 3D spin textures are visualized using cross-sectional spin-phase tubes derived from the mz = 0 isosurface, with color representing the local spin phase. The hopfion develops a side-localized rupture, where Bloch-point-like singular regions mediate the topological opening of the toroidal texture before relaxation toward the FM state. Numbered labels correspond to the sequential images marked in panel (a), with the star denoting the highest-energy image. (c) Representative configurations along the axial-collapse path. The hopfion core shrinks along the axial direction and collapses toward the FM state through a more symmetric channel. The visualization style matches panel (b), and the numbered labels correspond to the red markers in panel (a).
pathway, the hopfion tube develops a side-localized constriction and is severed by localized Bloch-point-like singular regions away from the central axis before relaxing toward the ferromagnetic state. This process is distinct from boundary escape, as the rupture occurs on the hopfion tube rather than through the sample boundary [5]. To our knowledge, such a lateral-rupture mechanism has not been reported for exchange-stabilized hopfions in frustrated-exchange atomistic models. The existence of two competing transition states with distinct activation barriers expands the picture of hopfion decay and illustrates how integrated atomistic workflows at the millionspin-per-image scale can uncover topological transition
channels that are difficult to resolve without a unified large-scale simulation workflow.
A unified Hamiltonian interface with multiple computational backends
SpinX is organized around a shared Hamiltonian interface rather than solver-specific implementations. Energy, effective-field and torque evaluations form the common computational core used by deterministic and stochastic Landau–Lifshitz–Gilbert (LLG) dynamics, Monte Carlo sampling, static minimization, spectral analysis and
6 minimum-energy-path calculations. As shown in Figure 1a, user-provided configurations, interaction tables and spin-state files are first processed by the simulation setup layer, which validates the model and constructs the requested workflow. The resulting spin Hamiltonian is then accessed by decoupled Python modules for dynamics, sampling, optimization, analysis and minimumenergy-path calculations, while storage modules record trajectories, checkpoints and observables. The computational modules differ in how they manipulate energies, fields, torques or trajectories, but they call the same shared Hamiltonian engine where physical model evaluations are required. The backend layer provides distinct numerical realizations of the same physical model. SpinX supports scalar and tensorial bilinear interactions, together with higherorder terms such as scalar biquadratic, scalar-chirality and quadruplet interactions in supported workflows. For pairwise bilinear terms, a unified 3 × 3 interaction tensor Jijαβ provides a common representation for isotropic exchange, symmetric anisotropic exchange and antisymmetric Dzyaloshinskii–Moriya interactions. This tensor representation maps directly onto local parameters such as the exchange scalar Jij and the Dzyaloshinskii–Moriya vector Dij (Figure 1b), allowing the same Hamiltonian to be evaluated by different backends without modifying the high-level physics modules that call it. For ideal crystals with translational symmetry, SpinX uses crystallographic sublattice decomposition (Figure 1b). The magnetic lattice is decomposed into discrete sublattice channels, conceptually analogous to colour channels in image processing, and translationally invariant pair interactions are encoded as cross-channel convolution kernels. In this representation, effective-field evaluation becomes a multi-channel tensor-convolution problem. The dense convolution backend handles compact regular kernels. The sparse convolution backend exploits structured sparse-shell interaction patterns. The FFT convolution backend treats extended translationally invariant interactions in reciprocal space. For the FFT backend, the interaction kernel is precomputed in q-space, multiplied by the transformed spin spectrum and transformed back to real space to obtain the effective field. Long-range dipolar interactions use the same reciprocal-space convolution principle through a dedicated dipole-FFT path in supported dynamics workflows. For non-ideal crystals, translational symmetry can be broken by chemical disorder, structural defects, vacancies or explicitly site-resolved couplings. In such heterogeneous settings, tensor convolution is no longer always the natural data structure. SpinX therefore provides a generalized pairlist backend, in which interactions are evaluated from explicitly indexed site-to-site couplings and their corresponding Hamiltonian entries. This representation preserves the same abstract Hamiltonian interface used by the convolutional backends while providing the structural flexibility required for site-dependent interactions and explicitly indexed interaction patterns.
The choice of computational backend is usercontrolled. SpinX does not assume a single universally optimal kernel. Instead, it exposes several numerically consistent realizations of the same Hamiltonian core within their overlapping capability domains. Dense convolution, sparse convolution, FFT convolution, dipole FFT and pairlist evaluation are suited to different combinations of lattice geometry, interaction range, floatingpoint precision and hardware platform. The performance benchmarks in Figure 3 provide practical guidance for selecting a backend for a given simulation regime.
Validation across dynamics, thermodynamics and spectroscopy
We next validate SpinX at the level of time integration, statistical sampling and dynamical correlations. These tests probe different components of the framework, including the deterministic LLG integrator, the mixedprecision execution mode, the Monte Carlo sampling routines and the spectral post-processing pipeline. Because all three benchmarks use the same model parsing layer and shared Hamiltonian engine, their agreement provides an end-to-end validation of the simulation stack. We first consider the simplest deterministic setting with an analytical reference: a single spin precessing under a uniform, constant magnetic field. Figure 2a shows the trajectory root-mean-square error (RMSE) of the Heun integration scheme over a 500 ps trajectory sampled every 1 ps, plotted as a function of the integration time step ∆t. In fp64, the error follows the expected time-discretization trend as the time step is reduced. The mixed-precision mode, in which effective-field evaluations are performed in fp32 while the spin state is stored and advanced in fp64, closely follows the fp64 baseline over the tested range. In contrast, the pure fp32 calculation develops a round-off error floor at small time steps, where further reducing ∆t no longer improves the trajectory error. This benchmark verifies the deterministic integrator and supports the use of mixed precision in large-scale production runs. Thermodynamic sampling is validated using a classical bcc Fe Heisenberg model with a multi-shell exchange profile. Figure 2b shows the temperature dependence of the magnetization M , Binder cumulant U4 , magnetic susceptibility χ and specific heat Cv . The magnetization decreases across the transition region, while χ and Cv show enhanced fluctuations in the same temperature window. The Binder-cumulant analysis provides a finitesize estimate of the ferromagnetic-to-paramagnetic transition: the inset shows the crossing behaviour of U4 for system sizes L = 12, 20, 28 and 36. This calculation tests the heat-bath Monte Carlo implementation, local energy and field updates, and the ability of SpinX to reproduce standard finite-temperature behaviour in an atomistic Heisenberg model. We further evaluate the dynamical correlation pipeline
7 by computing the transverse dynamic structure factor, S(q, ω) = Sxx (q, ω) + Syy (q, ω), at T = 5 K along the high-symmetry Γ–X–M –Γ path. The low temperature suppresses strong thermal broadening, while the transverse projection removes the static longitudinal magnetization component, making this benchmark suitable for direct comparison with the zero-temperature adiabatic magnon spectrum. As shown in Figure 2c, the dominant simulated spectral weight follows the adiabatic magnon spectrum (AMS) obtained from linear spin-wave theory. This agreement validates the generation of spindynamics trajectories, the space–time Fourier analysis and the reciprocal-space path construction used in the spectroscopy module. It also checks the sign conventions, coordinate transformations and energy-to-frequency conversion used to compare finite-temperature spin dynamics with the zero-temperature spin-wave reference. Together, these analytical, thermodynamic and spectroscopic benchmarks validate the main physical solvers in SpinX and show that the shared Hamiltonian engine gives consistent results across deterministic dynamics, thermal sampling and dynamical response calculations. Backend-to-backend consistency of energies and effective fields across the dense, sparse, FFT and pairlist backends is documented in the Supplementary Information.
GPU performance regimes and scaling
We then benchmark the computational performance of SpinX across backend choices, precision modes, hardware platforms and simulation workflows. All timings exclude JIT compilation and include device synchronization after execution. Throughput is reported as the number of spin sites processed per second, with the precise operation defined by each benchmark: field evaluations, Heun spin steps, Monte Carlo attempted spin updates or optimization iterations. For AMD MI250X, results are reported per Graphics Compute Die (GCD), whereas NVIDIA results use the GPU portion of the NVIDIA GH200 Grace Hopper Superchip. Figure 3a compares the single-device throughput of the four main Hamiltonian backends on a 200 × 200 × 200 nearest-neighbour exchange model. The solid bars report full Heun integration throughput, while the hatched bars report isolated effective-field evaluation throughput. These two measurements should be interpreted independently. The field-only benchmark measures the raw throughput of the Hamiltonian backend, whereas the full Heun benchmark additionally includes the predictorcorrector update, torque evaluation and spin normalization. Several backend and precision combinations exceed 109 spin sites per second for this short-range model, with the sparse convolution backend giving the highest throughput in the nearest-neighbour regime. The effect of precision is backend and platform dependent, while mixed precision provides a practical trade-off between speed, memory use and the accuracy behaviour estab-
lished in Figure 2a. The relative performance of the backends changes as the Hamiltonian becomes more complex. Figure 3b, left, shows the mixed-precision throughput as the number of exchange shells is increased from one to five. This benchmark illustrates why SpinX exposes multiple backend realizations instead of a single default kernel. Backends designed for compact nearest-neighbour interactions are not necessarily optimal for extended frustrated-exchange models, and the preferred backend changes with interaction range and memory-access pattern. Figure 3b, right, compares models with exchange only, exchange plus Dzyaloshinskii–Moriya interaction, and exchange plus DMI plus long-range dipolar fields. Adding DMI increases the local tensor structure of the field evaluation, while the dipolar term introduces a reciprocal-space contribution. The resulting overhead depends on the backend because different implementations reuse intermediate data, memory layouts and Fourier-space representations differently. We then evaluate single-node multi-GPU scaling for the mixed-precision dense-convolution workflow (Figure 3c). Two scaling modes are considered. In spatial domain decomposition, the lattice is partitioned across devices and each device advances a subdomain of the same simulation. In ensemble scaling, independent trajectories are distributed over devices. Ensemble scaling provides an upper-bound reference because it requires little inter-device communication, whereas spatial decomposition tests distributed execution of a single large simulation. Both modes scale across multiple devices on NVIDIA and AMD hardware. The translucent background bars (right axis) indicate the maximum system size that fits in memory, demonstrating that multi-GPU execution increases not only throughput but also the accessible spin capacity. Monte Carlo sampling has a different performance profile from deterministic time integration, since each update touches only a single spin and its neighbours rather than the entire lattice. Figure 3d reports fp64 heatbath and Metropolis throughput as a function of exchange shell count. One Monte Carlo sweep corresponds to one attempted update per spin. The cost of these routines is controlled not only by the number of interacting shells, but also by random-number generation, proposal construction, local field or energy evaluation and masked spin replacement. As a result, the shell-count dependence differs from the deterministic field-evaluation benchmarks. Finally, Figure 3e benchmarks fp64 static minimization under fixed-iteration settings. Velocity projection, Barzilai–Borwein and L-BFGS are compared as implementation workloads rather than as globally tuned optimizer rankings. The goal is to quantify the cost of executing common minimization updates on large spin systems with the shared Hamiltonian backend. Convergence speed depends on optimizer parameters, stopping criteria and the energy landscape, whereas the fixed-iteration
8 throughput reported here isolates the computational cost of each implementation. Together, these benchmarks show that SpinX does not rely on a single universally optimal backend. Its performance comes from exposing multiple numerically consistent realizations of the same Hamiltonian and allowing the user to select the backend appropriate for the lattice regularity, interaction range, precision mode, hardware platform and simulation workflow. This multi-backend design is essential for covering short-range dynamics, extended frustrated exchange, dipolar fields, Monte Carlo sampling, static optimization and multi-image path calculations within one framework.
Two competing annihilation channels of an exchange-stabilized hopfion
We finally use SpinX to study annihilation pathways of an exchange-stabilized hopfion in a minimal frustratedexchange model. The calculation is designed as a controlled atomistic test case rather than a material-specific prediction. The initial state is a metastable hopfion placed near the centre of a 100 × 100 × 100 spin lattice, and the final state is the ferromagnetic state. Open boundary conditions are used, with a surrounding ferromagnetic background separating the hopfion from the sample boundary. This geometry helps distinguish annihilation of the hopfion tube from simple escape through the outer boundary. Figure 4a shows two refined paths connecting the same initial and final states under the same Hamiltonian. The paths are constructed from string or NEB initializations and refined using GNEB. The blue path corresponds to a lateral-rupture channel, while the red path corresponds to an axial-collapse channel. The reaction coordinate is normalized independently for each path, and the energy is reported in units of the nearest-neighbour exchange J1 . The open stars mark the highest-energy images along the two paths. For the lateral-rupture channel, the highenergy region was further resolved by increasing the image density around the local maximum and applying local CI-GNEB refinement to better localize the saddle region. The magnified insets show the activation barriers of each channel, measured from the local metastable minimum on the hopfion side. The lower inset shows the saddleenergy difference |∆Es | ≈ 2.9 J1 between the two channels. The lateral-rupture pathway is shown in Figure 4b. Along this channel, the hopfion tube first develops a sidelocalized constriction. Near the highest-energy image, localized Bloch-point-like singular regions appear away from the central axis and rupture the toroidal texture from the side. The path then descends toward the ferromagnetic basin without an additional barrier comparable to the main saddle. This mechanism is distinct from boundary escape, since the rupture occurs on the hopfion tube itself rather than through motion of the entire
texture out of the sample. The axial-collapse pathway in Figure 4c provides a reference channel in the same model. In this case, the texture evolves more symmetrically along the path: the hopfion core shrinks along the axial direction, and Blochpoint-like singular regions appear near the central axis before the texture collapses toward the ferromagnetic state. This pathway resembles previously reported axial or core-collapse mechanisms of exchange-stabilized hopfions. In contrast, the lateral-rupture channel involves a side-localized rupture of the toroidal tube and leads to a distinct saddle morphology and activation barrier. The coexistence of these two paths shows that hopfion annihilation in the same frustrated-exchange model is not exhausted by a single pathway. The energy landscape contains both an axial-collapse route and a lateral-rupture route, each involving a different topologychanging event. This result demonstrates the utility of large-scale atomistic path calculations in SpinX, where many-image transition-path workflows can be carried out on million-spin configurations and used to resolve secondary annihilation channels that are difficult to identify from relaxation or direct dynamics alone.
METHODS Spin Hamiltonian and effective fields
SpinX represents each atomic magnetic moment as µi = µi µB Si , where µi is the moment magnitude in units of the Bohr magneton µB , and Si is a unit vector, |Si | = 1. Throughout this work, spin Hamiltonian parameters are stored internally in meV. The Bohr magneton is therefore written as meV/T
µB
≈ 5.7884 × 10−2 meV/T,
when converting energy gradients to magnetic fields in tesla. The total Hamiltonian is written as a sum of enabled contributions, H = Hbil + HZ + Hani + Hdip + Hho , where Hbil contains pairwise bilinear interactions, HZ is the Zeeman coupling, Hani denotes onsite anisotropy, Hdip is the dipolar interaction, and Hho collects supported higher-order terms. Terms that are not enabled in a given calculation contribute zero[21]. For an unordered pair representation, the pairwise bilinear Hamiltonian is X X αβ Hbil = − Jij Siα Sjβ , i>j α,β
9 αβ is a 3 × 3 interaction tensor in meV. Equivawhere Jij αβ , lently, writing Jij for the matrix with elements Jij X Hbil = − ST i Jij Sj . i>j
α,β
The factor 1/2 appears only in this ordered-pair representation and avoids double counting of reciprocal directed bonds. For scalar exchange and the Dzyaloshinskii-Moriya interaction, the tensor action is defined by X αβ β α Jij Sj = Jij Sjα + (Sj × Dij ) . β
Equivalently, y z 0 Dij −Dij z x 0 Dij Jij = Jij I + −Dij . y x Dij −Dij 0
With this convention, the energy of an unordered interacting pair is Hij = −Jij Si · Sj − Dij · (Si × Sj ), with Dji = −Dij . Symmetric anisotropic exchange is αβ included through the symmetric part of Jij . The Zeeman contribution is X meV/T ext HZ = − µi µB Bi · Si . i
The onsite anisotropy used by SpinX has the form XX 2 Hani = − Kiℓ (Si · n̂iℓ ) , ℓ
where Kiℓ is stored in meV and n̂iℓ is the anisotropy axis. Additional supported terms include dipolar interactions and selected higher-order interactions, including scalar biquadratic, scalar-chirality and quadruplet terms in workflows where they are enabled. Not all higher-order terms are available in every computational backend. The effective field [2]is defined in the atomistic spindynamics convention as Beff i =−
1
j,β
and then converts it to a physical field through
In the implementation, unordered pair interactions are expanded into reciprocal directed bonds, and the equivalent ordered-pair form is 1 X X αβ α β βα αβ Hbil = − Jij Si Sj , Jji = Jij . 2 i,j
i
SpinX first evaluates the reduced bilinear field in meV, X αβ β bbil Jij Sj , i,α =
∂H
meV/T ∂S i µi µB
,
where Beff i is in tesla when H is in meV. For the bilinear tensor term this gives X αβ β 1 bil Jij Sj . Bi,α = meV/T µi µB j,β
Bbil i =
bbil i meV/T
.
µi µB
The same effective-field convention is used by the dynamics, static optimization, spectral-analysis and minimumenergy-path modules. Monte Carlo sampling uses the same Hamiltonian parameters through local fields or local energy differences, depending on the update scheme. Backend representations
For translationally invariant crystals, SpinX decomposes the magnetic lattice into crystallographic sublattice channels. A spin on sublattice a in unit cell R is denoted Sa (R). The reduced bilinear field on sublattice a can be written as a multi-channel convolution[22], X X X αβ bα Kab (δ)Sbβ (R + δ), a (R) = b
δ
β
αβ where Kab (δ) stores the interaction tensor between sublattices a and b at displacement δ. Dense convolution, sparse convolution and FFT convolution are different numerical realizations of this same expression. The dense backend stores compact regular kernels and evaluates the reduced field using convolution primitives. The sparse backend stores only non-zero displacement entries and evaluates structured finite-range kernels without forming a dense stencil. The FFT backend is used for fully periodic translationally invariant interactions by transforming the spin field to reciprocal space, multiplying by the precomputed kernel spectrum and transforming back, X X αβ b̂α K̂ab (q)Ŝbβ (q). a (q) = b
β
Long-range dipolar fields are evaluated using the same reciprocal-space convolution principle in workflows where the dipolar FFT path is enabled[23, 24]. Finite regular supercells with open boundaries can still be handled by the regular sublattice backends through the appropriate padding or truncation rule. When translational regularity is broken by chemical disorder, structural defects or explicitly site-resolved couplings, SpinX uses a generalized pairlist backend. Each pairlist entry stores the source site, target site and corresponding Hamiltonian parameters. The reduced bilinear field is evaluated as X X bα Jpαβ Sjβp , i = p:ip =i β
10 where p indexes pairlist entries, ip is the target site and jp is the source neighbour associated with entry p. This backend preserves the same Hamiltonian interface while removing the requirement of translational symmetry.
Atomistic spin dynamics
Deterministic dynamics are integrated using the Landau-Lifshitz-Gilbert equation[25, 26], dSi γ eff Si × Beff , =− i + αi Si × Si × Bi 2 dt 1 + αi where γ is the gyromagnetic ratio and αi is the Gilbert damping parameter. The main time integrator used in this work is a normalized Heun scheme. Defining the right-hand side of the LLG equation as fi (S), the predictor step is S̃i = N [Sni + ∆t fi (Sn )] , followed by the corrected step ∆t n+1 n n Si = N Si + fi (S ) + fi (S̃) . 2 Here N [x] =
x |x|
denotes projection back to the unit sphere.
Stochastic Landau-Lifshitz-Gilbert simulations
For stochastic LLG simulations, a Gaussian thermal field is added to the effective field[27, 28], eff th Beff i → Bi + Bi .
In discrete time, the thermal field is written as Bth,n = σi ηin , i where the Cartesian components of ηin are independent standard normal random variables. For the LLG convention used here, s 2αi kB T σi = . meV/T |γ|µi µB ∆t Here kB T is evaluated in meV. This noise amplitude corresponds to the LLG convention used in the equation above. Random fields are generated reproducibly from explicit JAX random keys. In mixed-precision mode, effective-field evaluation is performed in fp32, while the spin state and time integration are kept in fp64.
Monte Carlo sampling
SpinX implements heat-bath and Metropolis Monte Carlo updates for supported spin Hamiltonians[29, 30]. The heat-bath update is used when the local energy is linear in the updated spin, as in bilinear exchange models with optional Zeeman terms. The Monte Carlo benchmarks in this work use this class of Hamiltonians. In a local update, the local field acting on site i is held fixed while a new spin direction is sampled or proposed. For the heat-bath update, the spin is sampled from P (Si ) ∝ exp βbloc i · Si , where bloc is the reduced local field in meV and β = i (kB T )−1 , with kB T in meV. Equivalently, using the physical local field Bloc i , meV/T
bloc i = µi µB
Bloc i .
In the Metropolis update, a trial spin S′i is accepted with probability Pacc = min [1, exp(−β∆Hi )] , where ∆Hi is the local energy change in meV associated with the trial move. For compatible regular lattices, SpinX uses colourbased site-major updates so that spins updated simultaneously do not interact within the selected interaction range. The number of colours is chosen according to the neighbour range of the active Hamiltonian, enabling parallel updates while preserving the local dependency structure of the Monte Carlo move. Thermodynamic and dynamical observables
SpinX stores the total magnetic moment as X Mtot = µi Si , i
in units of µB . For equal magnetic moments, the plotted reduced magnetization is 1 X m= Si , M = |m|. N i For the temperature scans shown in Figure 2, the plotted magnetization, susceptibility and specific heat are scaled for visual comparison. The Binder cumulant is evaluated from the sampled magnetization magnitude as U4 = 1 −
⟨M 4 ⟩ . 3⟨M 2 ⟩2
For the plotted finite-temperature curves, SpinX uses the following fluctuation-based estimators before rescaling for visual comparison, N χrep = ⟨M 2 ⟩ − ⟨M ⟩2 , T
11 and Cv,rep =
N ⟨e2 ⟩ − ⟨e⟩2 , T2
e=
H . N
These reported quantities are used consistently for the scaled curves in Figure 2. For dynamical spectra, the time average is first subtracted from each spin component,
For |ui | = 0, the exponential-map expression is evaluated in its limiting form. The minimum-energy-path calculations in this work use the exponential-map retraction unless otherwise stated. Minimum-energy paths are computed using string and geodesic nudged elastic band methods[36–38]. An image n is a full spin configuration (n)
Yn = {Si }N i=1 .
δSiα (t) = Siα (t) − ⟨Siα ⟩t . The spatial Fourier amplitude is then computed as 1 X Mα (q, t) = √ µi δSiα (t) exp(−iq · ri ). N i A Hann window is applied to Mα (q, t) along t before the time Fourier transform to reduce spectral leakage. The spectral tensor is formed as Sαβ (q, ω) = Mα∗ (q, ω)Mβ (q, ω). Overall normalization factors associated with the time window, window function and trajectory length are omitted because the spectra in Figure 2 are plotted in arbitrary units. For the ferromagnetic spectra shown in Figure 2, the equilibrium magnetization is along the z axis and the plotted transverse intensity is S⊥ (q, ω) = Sxx (q, ω) + Syy (q, ω). More generally, SpinX can also form a q-transverse projection using Pαβ = δαβ − q̂α q̂β . The numerical Fourier transform returns ordinary frequencies f , and the plotted energy is E = hf = ℏω. Static optimization and minimum-energy paths
Static minimization is performed on the product manifold of unit spin spheres[31]. Given the Euclidean gradient ∇i H =
∂H , ∂Si
the corresponding tangent gradient is
Optimization algorithms such as velocity projection, Barzilai-Borwein, L-BFGS and Adam update tangent directions and retract the spin back to the unit sphere[32– 35]. SpinX supports both normalized and exponentialmap retractions. The normalized retraction is Si + ui , |Si + ui |
whereas the exponential-map retraction is Rexp Si (ui ) = cos |ui |Si + sin |ui |
ui , |ui |
d(Yn , Ym ) =
" X
arccos2
(n) (m) Si · Si
# 1/2
.
i
In the string method, each image is first displaced by the projected physical force on the spin manifold. The path is then reparameterized along the geodesic reaction coordinate so that the image distribution follows the chosen arc-length spacing. This separates the relaxation of the path shape from the maintenance of image spacing and avoids introducing an explicit spring force. For image n, the GNEB physical force is the negative energy gradient projected onto the tangent space of the spin manifold. The GNEB force separates this physical force into a component perpendicular to the path tangent and a spring force parallel to the path, spring FGNEB = Fphys . n n,⊥ + Fn,∥
For climbing-image refinement, the highest-energy image is driven toward the saddle region by reversing the tangential component of the physical force, phys FCI − 2 Fphys · τ̂n τ̂n . n = Fn n The vector τ̂n is the normalized path tangent in image space, and all forces are projected onto the spin tangent space. Endpoints can be pinned, and image redistribution is applied when needed to maintain path resolution.
Performance measurements
Gi = ∇i H − (∇i H · Si )Si .
Rnorm (ui ) = Si
Distances between images are evaluated using the geodesic metric on the spin manifold,
ui · Si = 0.
Performance benchmarks report throughput in spin sites processed per second. For effective-field benchmarks, one spin-site operation corresponds to one evaluated effective field for one spin. For Heun dynamics, one spin-site operation corresponds to one full Heun time step for one spin, including the predictor-corrector update. For Monte Carlo, throughput is reported as attempted spin updates per second. For static optimization, throughput is reported as spin-site updates per fixed optimization iteration. All timings exclude JIT compilation and include device synchronization after execution. AMD MI250X results are reported per Graphics Compute Die (GCD), whereas
12 NVIDIA GH200 results are reported per GPU. The hardware and software environment used for the benchmarks is summarized in Table I. Benchmark system sizes, precision modes, backend choices and workload definitions are specified in the corresponding figure captions. TABLE I. Hardware and software environment for performance benchmarks. Platform Accelerator Software stack NVIDIA GH200 GPU JAX 0.9.2, CUDA 12 AMD MI250X GCD JAX 0.9.2, ROCm 7
Hopfion calculation
The exchange-stabilized hopfion in Figure 4 is computed using a minimal four-shell frustrated-exchange model on a simple-cubic lattice[39]. The model is used as a controlled atomistic realization of a frustratedexchange Hamiltonian that stabilizes Hopf solitons, rather than as a material-specific parameter set. The Hamiltonian used for this calculation is X H=− Jij Si · Sj , i>j
with non-zero couplings on the first four neighbour shells. We use J1 ≃ 12.4 meV
and a lateral-rupture channel. The paths shown in Figure 4 are assembled from GNEB-refined segments with pinned endpoints. For the lateral-rupture channel, the image density near the highest-energy region is increased and a local CI-GNEB refinement is applied to the saddle region. In this local refinement, the maximum perpendicular force on the climbing image is reduced below 5 × 10−6 J1 per spin. Outside the locally refined CI region, the maximum perpendicular force along the path is kept below ∼ 10−3 J1 per spin. Here forces are measured as energy gradients with respect to dimensionless spin coordinates. Bloch-point-like regions are identified from localized peaks in a discretized topological-charge diagnostic on the magnetic lattice. The spin textures in Figure 4 are rendered as streamtubes traced on the mz = 0 isosurface, with colours representing the local in-plane spin phase. These markers are used to visualize topology-changing regions along the transition path[9, 40].
I.
CODE AND DATA AVAILABILITY
The source code of SpinX v0.1.0, corresponding to the version used for the calculations in this work, is publicly available at: https://github.com/QichenXuResearch/SpinX. This repository is released as an initial review snapshot and will continue to evolve during and after the review process. The repository includes the core SpinX package, command-line interface, validation tests and minimal executable examples.
as a representative energy scale, and J2 = 0.188, J1
J3 = −0.274, J1
J4 = −0.161. J1
All energies in Figure 4 are reported in units of J1 . The initial hopfion configuration with Hopf index QH = 1 is prepared from a toroidal Hopf-map ansatz, S(r) = (cos Φ(r) sin Θ(r), sin Φ(r) sin Θ(r), cos Θ(r)) , where Θ(r) and Φ(r) encode the linked-fibre structure of the Hopf map. The ansatz is placed near the centre of a 100 × 100 × 100 simple-cubic lattice with open boundary conditions and a surrounding ferromagnetic region large enough to reduce spurious boundary effects on the annihilation mechanism. The ansatz is relaxed using static minimization with the velocity-projection method to obtain the metastable hopfion endpoint. The final endpoint is the uniform ferromagnetic state under the same Hamiltonian and boundary conditions. Transition paths between the relaxed hopfion and the ferromagnetic state are generated from multiple trial paths and refined using string and geodesic nudged elastic band calculations. Two reproducible annihilation channels are retained for analysis: an axial-collapse channel
II.
ACKNOWLEDGMENTS
The authors thank Filipp N. Rybakov (Uppsala University), Olle Eriksson (Uppsala University), Manuel Pereiro(Uppsala University), Pavel Bessarab (Linnaeus University) and Liuzhen Yang for many fruitful discussions. We also thank Johan Hellsvik (KTH, Dardel), Olle Eriksson (CSC, LUMI) for support with GPU resources. The authors used AI-assisted tools to improve the language of the manuscript. Financial support from the Swedish Research Council (Vetenskapsrådet, VR; Grant No. 2024-04986), and the Knut and Alice Wallenberg Foundation (KAW; Grant No. 2022.0108), is acknowledged. The Wallenberg Initiative Materials Science for Sustainability (WISE) funded by the Knut and Alice Wallenberg Foundation and The Swedish e-Science Research Centre (SeRC) are also acknowledged. The computations/data handling were enabled by resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS), partially funded by the Swedish Research Council through grant agreement no. 2022-06725.
13 III.
COMPETING INTERESTS
All authors declare that they have no conflicts of interest.
[1] V. P. Antropov, M. I. Katsnelson, M. van Schilfgaarde, and B. N. Harmon, Physical Review Letters 75, 729 (1995). [2] V. P. Antropov, M. I. Katsnelson, B. N. Harmon, M. van Schilfgaarde, and D. Kusnezov, Physical Review B 54, 1019 (1996). [3] Q. Xu, Z. Shen, A. Edström, Z. Lu, A. Bergman, D. Thonig, W. Yin, O. Eriksson, and A. Delin, npj Computational Materials 11, 56 (2025). [4] M. Sallermann, H. Jónsson, and S. Blügel, Physical Review B 107, 104404 (2023). [5] I. S. Lobanov and V. M. Uzdin, Physical Review B 107, 104405 (2023). [6] Q. Xu and A. Delin, npj Computational Materials (2026). [7] Q. Xu, I. P. Miranda, M. Pereiro, F. N. Rybakov, D. Thonig, E. Sjöqvist, P. F. Bessarab, A. Bergman, O. Eriksson, P. Herman, et al., Physical Review Research 5, 043199 (2023). [8] Q. Xu, Z. Shen, M. Pereiro, E. Sjöqvist, P. Herman, O. Eriksson, and A. Delin, Communications Physics 6, 239 (2023). [9] V. M. Kuchkin, A. Haller, Š. Liščák, M. P. Adams, V. Rai, E. P. Sinaga, A. Michels, and T. L. Schmidt, Physical Review Research 7, 013195 (2025). [10] Y. Liu, R. K. Lake, and J. Zang, Physical Review B 98, 174437 (2018). [11] X. S. Wang, A. Qaiumzadeh, and A. Brataas, Physical Review Letters 123, 147203 (2019). [12] N. Kent, N. Reynolds, D. Raftrey, I. T. G. Campbell, S. Virasawmy, S. Dhuey, R. V. Chopdekar, A. HierroRodriguez, A. Sorrentino, E. Pereiro, S. Ferrer, F. Hellman, P. Sutcliffe, and P. Fischer, Nature Communications 12, 1562 (2021). [13] F. Zheng, N. S. Kiselev, F. N. Rybakov, L. Yang, W. Shi, S. Blügel, and R. E. Dunin-Borkowski, Nature 623, 718 (2023). [14] O. Eriksson, A. Bergman, L. Bergqvist, and J. Hellsvik, Atomistic spin dynamics: foundations and applications (Oxford university press, 2017). [15] Q. Xu, O. Eriksson, and A. Delin, arXiv preprint arXiv:2309.17367 (2023). [16] R. F. Evans, W. J. Fan, P. Chureemart, T. A. Ostler, M. O. Ellis, and R. W. Chantrell, Journal of Physics: Condensed Matter 26, 103202 (2014). [17] G. P. Müller, M. Hoffmann, C. Dißelkamp, D. Schürhoff, S. Mavros, M. Sallermann, N. S. Kiselev, H. Jónsson, and S. Blügel, Physical Review B 99, 224414 (2019). [18] A. Sabne, “Xla : Compiling machine learning for peak performance,” (2020).
[19] J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, Y. Katariya, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Milne, and Q. Zhang, “JAX: composable transformations of Python+NumPy programs,” (2018). [20] M. Pajda, J. Kudrnovský, I. Turek, V. Drchal, and P. Bruno, Physical Review B 64, 174402 (2001). [21] W. Heisenberg, Zeitschrift für Physik 49, 619 (1928). [22] J. W. Cooley and J. W. Tukey, Mathematics of Computation 19, 297 (1965). [23] D. V. Berkov and N. L. Gorn, Journal of Magnetism and Magnetic Materials 122, 204 (1993). [24] C. J. Garcı́a-Cervera, Journal of Computational Physics 190, 1 (2003). [25] L. D. Landau and E. M. Lifshitz, Physikalische Zeitschrift der Sowjetunion 8, 153 (1935). [26] T. L. Gilbert, IEEE Transactions on Magnetics 40, 3443 (2004). [27] W. F. Brown, Physical Review 130, 1677 (1963). [28] R. F. L. Evans, W. J. Fan, P. Chureemart, T. A. Ostler, M. O. A. Ellis, and R. W. Chantrell, Journal of Physics: Condensed Matter 26, 103202 (2014). [29] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, The Journal of Chemical Physics 21, 1087 (1953). [30] M. Creutz, Physical Review D 36, 515 (1987). [31] P.-A. Absil, R. Mahony, and R. Sepulchre, Optimization Algorithms on Matrix Manifolds (Princeton University Press, Princeton, 2008). [32] E. Bitzek, P. Koskinen, F. Gähler, M. Moseler, and P. Gumbsch, Physical Review Letters 97, 170201 (2006). [33] J. Barzilai and J. M. Borwein, IMA Journal of Numerical Analysis 8, 141 (1988). [34] D. C. Liu and J. Nocedal, Mathematical Programming 45, 503 (1989). [35] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” (2015), arXiv:1412.6980 [cs.LG]. [36] W. E, W. Ren, and E. Vanden-Eijnden, Physical Review B 66, 052301 (2002). [37] G. Henkelman, B. P. Uberuaga, and H. Jónsson, The Journal of Chemical Physics 113, 9901 (2000). [38] P. F. Bessarab, V. M. Uzdin, and H. Jónsson, Computer Physics Communications 196, 335 (2015). [39] F. N. Rybakov, N. S. Kiselev, A. B. Borisov, L. Döring, C. Melcher, and S. Blügel, APL materials 10 (2022). [40] B. Berg and M. Lüscher, Nuclear Physics B 190, 412 (1981).