Extreme-Scale Atomistic Simulation of Real-Temperature Magnetic Skyrmion Dynamics by Coupled Spin-Lattice Modeling Pin Chen1 , Cheng-bing Chen2 , Hai Liu1 , Yuewen Huang1 , Kangyou Zhong1 , Hai-Jun Zhao3 , Liu-Liu Han4,5 , Guixin Guo1 , Jiang Li1 , Dan Huang1 , Ben Xu2,* , Yutong Lu1,6,*
arXiv:2606.14073v1 [cs.DC] 12 Jun 2026
1
Sun Yat-sen University, National Supercomputer Center in Guangzhou, Guangzhou, China 2 Graduate School of China Academy of Engineering Physics, Beijing 100089, China 3 Key Laboratory of Quantum Materials and Devices (MOE), School of Physics, Southeast University, Nanjing 211189, China 4 Suzhou Laboratory, Suzhou 215000, China 5 State Key Laboratory of Powder Metallurgy, Central South University, Changsha 410083, China 6 Sun Yat-sen University, National Supercomputer Center in Shenzhen, Shenzhen, China *
Corresponding authors: [email protected], [email protected]
Abstract—Real-temperature topological magnetic dynamics in functional materials is governed by coupled lattice and spin evolution, yet remains inaccessible to predictive simulation at device-relevant scales. As a flagship example, thermally driven helix-to-skyrmion transformation in FeGe requires atomistic resolution, explicit lattice motion, and micrometer-scale domains to resolve device-scale topological texture formation. We combine a spin-constrained density-functional-theory-trained neuroevolution potential with a structure-preserving spin–lattice integrator within one machine-learned framework. Architecture– specific optimizations, kernel fusion, SVE2 vectorization, and NUMA-aware data layout deliver a seven orders-of-magnitude speedup over prior spin-aware methods. Deployed on LineShine exascale supercomputer, the full application scales to 12.45 million CPU cores with 89.7% weak-scaling efficiency, enabling simulations of 1.34 trillion atoms and an equal number of spins while reaching 48.5 PFLOPS in double precision. The simulations directly resolve real-temperature skyrmion nucleation and reorganization at previously inaccessible scales, establishing a new regime for predictive simulation of coupled spin–lattice topological magnetic dynamics. Index Terms—magnetic skyrmions, machine-learning interatomic potential, spin–lattice dynamics, extreme-scale simulation
1. J USTIFICATION FOR THE ACM G ORDON B ELL P RIZE We performed the largest atomistic simulation of a magnetic material to date: a 1.34-trillion spins and 1.34-trillion atoms coupled spin–lattice system executed on 12.45 million ARM CPU cores, achieving first-principles accuracy, reaching a sustained performance of 48.5 PFLOPS in double precision, and sustaining 89.7% weak-scaling efficiency at full system scale. This unprecedented scale enables the first direct atomistic observation of skyrmion dynamics at device-relevant length scales. 2. P ERFORMANCE ATTRIBUTES
Performance attribute
This submission
Category of achievement Type of method used Results reported on basis of Precision reported System scale Measurement mechanism
Peak performance, scalability Explicit Whole application Double precision Full-scale machine (20,480 nodes) Timers and FLOP counting
3. OVERVIEW OF THE P ROBLEM Magnetic skyrmions have emerged as one of the most compelling topological spin textures for future low-power information technologies. This is because they can combine nanoscale or mesoscale manipulability with nontrivial topology, therefore enabling rich dynamical functionality. More broadly, Nobel Prize–recognized advances in spin-dependent transport and topological phases have established spin functionality and topology as two defining themes of modern condensed-matter physics [1], [2]. Among the known skyrmion-hosting material classes, noncentrosymmetric chiral magnets provide one of the most important platforms because the competition between Heisenberg exchange(J) and the Dzyaloshinskii-Moriya interaction(D) naturally sets the characteristic magnetic length scale in the experimentally useful tens-of-nanometers range [4], [5]. Additionally, the three-dimensional nature of bulk DMI materials enables hosting of skyrmion strings and more complex topological particles like hopfions, offering richer physics and greater computational complexity than 2D systems. FeGe is a prototypical example: its characteristic helix/skyrmion size is about 70–80 nm [6], while its bulk critical helimagnetic ordering temperature lies close to room temperature, around 278 K, with a more precise single-crystal determination giving 279.1 K [7]. These two features make FeGe especially attractive for studying real-temperature skyrmion formation and dynamics.
Fig. 1: Schematic of a magnetic helix (a) and skyrmion (b) in a chiral magnet. Color coding represents the z component of magnetization. Arrows are plotted for every 200 atoms; the typical helix period and radius of skyrmion in FeGe is ∼70–80 nm [3].
linear scaling with system size. (2) a symplectic spin-lattice dynamics integrator designed for magnetic machine learning library of interatomic potential (MMLIP) that propagates the coupled lattice and spin degrees of freedom with long-time numerical stability, and including the longitudinal fluctuation of magnetic moment and (3) an extreme-scale implementation of spin-lattice dynamics (termed magnetic molecular dynamics) within LAMMPS optimized for a many-core ARM supercomputer, enabling simulations on 20,480 nodes (12.45 million CPU cores). Our method enables atomistic simulations of up to 1.34 × 1012 magnetic atoms, corresponding to a physical domain of 3.02 µm × 2.41 µm × 2.41 µm and reaching experimentally relevant device dimensions. In the present model, each atom carries one local magnetic moment, so the number of spins is equal to the number of atoms. Throughout the paper, system size is reported by atom count unless otherwise noted.
However, these same features also make predictive simulation exceptionally demanding. First, a physically meaningful calculation must extend far beyond a single texture unit and instead contain many competing periods, boundaries, and fluctuation–mediated nucleation channels, which immediately push the system size toward the mesoscopic regime. Second, because the transformation occurs near experimentally relevant temperatures, temperature cannot be treated as a purely phenomenological parameter: thermal fluctuation, lattice motion, and spin–lattice energy exchange must be resolved explicitly. Third, the characteristic length scale is controlled by the delicate balance between exchange and DMI (J/D), so even modest errors in the effective J/D ratio can shift the helix pitch, skyrmion size, and phase competition. Finally, the relevant evolution unfolds on nanosecond time scales, which makes brute-force high-fidelity simulation prohibitively expensive unless both the potential evaluation and the coupled time integration can be executed at extreme throughput. Therefore, the evolution from a low-temperature helix-dominated state toward thermally stabilized skyrmionic textures in FeGe represents a central but still incompletely resolved problem [5], [8]. The central challenge is therefore not simply to simulate a larger or hotter magnetic system, but to resolve, within one physically consistent framework, how real-temperature lattice motion and spin dynamics jointly drive the transformation from helical order to skyrmionic textures and their subsequent reorganization on experimentally relevant scales. We address this challenge through three key components: (1) a magnetic machine-learning interatomic potential based on an extended spin-included neuro-evolution potential (NEPSPIN), developed from the standard NEP [9], which learns the spin-lattice coupled potential energy surface from constrained density functional theory data of magnetic excite configurations [10], and replaces repeated electronic-structure evaluations with local atomic inference, reducing the computational cost from cubic-scaling first-principles calculations to overall
4. C URRENT S TATE OF THE A RT State-of-the-art simulation of thermally driven topological magnetic phenomena requires a unified spin–lattice framework that treats atomic positions and local magnetic moments on an equal footing, enabling consistent computation of forces and magnetic torques for time integration. This requirement is particularly critical for modeling skyrmion formation, where the underlying physics spans multiple scales: from atom-resolved lattice distortions and texture reorganization at 10–100 nm, to collective dynamics across mesoscopic or even micrometer-sized domains over extended timescales. In the context of Gordon Bell-class computing, these scientific demands translate into stringent system requirements: the online evaluation of spin–lattice dynamics must be local, massively parallel, communication-efficient, and acceleratorfriendly, or production-scale simulation of real-temperature topological nucleation remains out of reach. Table I summarizes representative approaches across this landscape, including atom-only ML-MD baselines from the 2020 Gordon Bell Prize (DeepMD [12]) and the 2025 Gordon Bell finalist (XS-
2
TABLE I: Performance comparison of representative atomistic and spin–lattice simulation frameworks at scale. ↑ = higher is better; ↓ = lower is better. Abbreviations: D.O.F. = degrees of freedom; TtS = time-to-solution; A = atom-only; S+A = coupled spin–lattice; Norm. TtS = TtS normalized by model parameter count, i.e. s/(atom·parameter·step). ∗ Atom-only method without magnetic degrees of freedom; included as Gordon Bell performance baselines. † Projected from 1/16 machine (9,936 nodes) to full machine; see [11] for details. ‡ Mixed precision (FP64/FP32/BF16); not a pure FP64 result. All other sustained FLOPS entries are double precision. Work
Year
Method
System
D.O.F.
# atoms/ spins
# CPU cores
# GPUs
Machine
Sustained [FLOPS] ↑
TtS ↓ [s/step/atom]
Norm. TtS ↓ [s/(atom·param ·step)]
Baseline [12]∗ (double) Baseline [12]∗ (double) Guo et al. [11]∗ (double) Guo et al. [11]∗ (double) Guo et al. [11]∗† (double) Guo et al. [11]∗† (double) Razakh et al. [13]∗ (mixed)
2020 2020 2022 2022 2022 2022 2025
DeepMD DeepMD DeepMD DeepMD DeepMD DeepMD XS-NNQMD
Cu H2 O Cu H2 O Cu H2 O PbTiO3
A A A A A A A
127 M 679 M 3.4 B 3.9 B 17.3 B 24.9 B 1.23 T
27.3 K 27.3 K 27.3 K 27.3 K 7.63 M 7.63 M 1.04 M
27.3 K 27.3 K 27.3 K 27.3 K — — 60 K
Summit Summit Summit Summit Fugaku Fugaku Aurora
91 P 80 P 43.7 P 46.3 P 119 P 124.8 P 1.87 E‡
8.1 × 10−10 8.1 × 10−10 1.1 × 10−10 — 4.1 × 10−11 — —
— — — — — — 1.876 × 10−15
Dudarev et al. [14] Yu et al. [15] Yang et al. [16]
2008 2024 2024
Spilady spinGNN++ DeepSPIN
Fe — —
S+A S+A S+A
128 K 4,608 11,328
— — 5.6k
— 16 —
— — —
— — —
— — 4.52 × 10−3
— — —
This work
2026
NEP + SPIN
FeGe
S+A
168 B 1.34 T
12.45 M
—
LineShine
48.5 P 43.3 P
1.79 × 10−11 1.82 × 10−11
3.97 × 10−15 4.05 × 10−15
The remaining difficulty, however, is now less conceptual than computational. At inference of spin and lattice configuration and energy, the practical cost depends strongly on descriptor locality, kernel fusion, memory-access efficiency, and hardware suitability. In this regard, models such as DeepSPIN [16] represent an important step toward explicit spin-aware learning, but implementations based on generalpurpose deep-learning frameworks may introduce substantial inference overhead in long-time, large-system simulations. By contrast, NEP-style local descriptors [29] are much closer to the throughput and scalability requirements of production HPC workflows, but in their current form they do not explicitly include spin degrees of freedom or generate magnetic torques.
NNQMD [13]), existing spin–lattice methods, and the present work. In this perspective, existing methods still trade physical fidelity with computational scale. Continuum micromagnetics can reach large domains, but temperature usually enters only through effective noise or fitted parameters, so lattice heating and lattice-to-spin energy transfer are not treated explicitly [17]. Atomistic spin dynamics resolves individual moments, yet the lattice is often frozen or replaced by a thermal bath [18]. A more complete route is on-the-fly first-principles spin– lattice dynamics. In practice, however, its cost is set not only by electronic-structure evaluations, but also by the need to follow nonadiabatic transitions among multiple electronic or spin states, often in a surface-hopping-type framework [19], [20]. More scalable alternatives include DFT-parameterized spin Hamiltonians, spin- or magnetic-cluster expansions, and classical spin–lattice dynamics [14], [21]–[24]. These methods are efficient and interpretable, but their magnetic couplings are usually fixed in advance rather than updated self-consistently with temperature-driven atomic motion. This limits transferability when the structure changes strongly, the pathway is complex, or real-temperature evolution is essential.
The central state-of-the-art gap is therefore not simply one of accuracy or scale in isolation, but of closing physical completeness and exascale deployability within the same online potential-evaluation pipeline. Among the currently available spin-aware ML approaches, SpinGNN [26] explicitly reports million-atom GPU-parallel spin–lattice simulations, whereas DeepSPIN and SpinGNN++ [15], [16] place stronger emphasis on physical fidelity and model expressiveness. Publicly available papers still rarely report Gordon-Bell-style metrics such as peak FLOPS or seconds per step per atom, which makes their exascale deployability difficult to assess directly; as shown in Table I, most spin–lattice methods lack reported performance metrics entirely.
Recent spin-aware machine-learned models close this gap by describing energy, force, and magnetic torque within a unified representation, through explicit incorporation of magnetic environments into the descriptor or neural architecture. This category includes spin–lattice hybrid ML potentials [25], magnetic atomic cluster expansion (ACE) variants [22], graphbased magnetic networks [15], [26], and other torque-aware models [16], [27], [28] developed for non-collinear and finitetemperature settings. Their promise is substantial: they offer a route toward combining anharmonic lattice dynamics, complex magnetic interactions, and local spin–lattice feedback within a single differentiable surrogate.
This work targets precisely this gap. By extending the NEP descriptor to incorporate non-collinear spin degrees of freedom, coupling it with a spin–lattice integrator designed to preserve the key geometric structure of the coupled dynamics and maintain long-time numerical stability, and co-designing the entire pipeline for a many-core ARM supercomputer, we achieve coupled spin–lattice simulations of up to 1.34 T atoms at a time-to-solution of 1.79 × 10−11 s/step/atom in double
3
precision. The following sections describe the algorithmic and implementation innovations that make this possible.
multiple force/field reevaluations within a single time step, the spin update must be scheduled last among time-integration operations.
5. I NNOVATIONS R EALIZED
B. Implementation Innovations
Porting NEPSPIN-based spin–lattice dynamics to the target many-core ARM platform exposes several challenges rooted in the mismatch between the model’s computational patterns and the hardware micro-architecture: the deep 16-domain NUMA hierarchy demands careful affinity management; the absence of a shared L3 cache penalizes data reuse across memory-bound descriptor kernels; and the dominance of GEMV-like inference operations underutilizes the SME matrix units designed for outer-product GEMM workloads. Addressing these challenges motivates the algorithmic and implementation innovations described in the following and illustrated in Fig. 2.
1) Unified spin–lattice execution in LAMMPS: The proposed framework is realized in LAMMPS by coupling the NEPSPIN-based force-field evaluation with the spin–lattice time integrator equipping the above self-consistent updating scheme. At each MD step, the NEPSPIN module performs descriptor construction and inference on the unified energy surface E(R, S), where R and S denote the atomic positions and spin configurations respectively, producing both lattice and spin driving terms for the subsequent coupled update. The spin–lattice integrator then advances the state (R, S) with atomic forces and magnetic torques in a single update procedure, followed by halo exchange for the next step. This implementation eliminates the need for separate lattice and magnetic solvers, and preserves consistency between force evaluation, spin evolution, and parallel data movement throughout the simulation. 2) Fused multi-physics force kernel: In the original NEPSPIN implementation, forces and magnetic torques are evaluated by three independent functions, each performing a full neighbor-list traversal with its own distance computation, cutoff test, and Chebyshev basis recurrence. Because all three share the same radial cutoff, we fuse them into a single kernel that walks the neighbor list once and evaluates all force and torque contributions per pair in a single pass. This eliminates two redundant neighbor traversals (∼300 instructions per neighbor × ∼110 neighbors per atom) and two redundant Chebyshev recurrences (244 instructions per neighbor each), reducing the fused-kernel runtime by 23.9%. 3) SVE2 predicated vectorization with pre-staging: The irregular neighbor list—where each atom has a different number of valid neighbors after cutoff filtering—prevents direct vectorization of the inner loop. We address this with a twophase pre-staging strategy. In Phase A, a scalar pass filters valid neighbors through cutoff tests and packs their coordinates, distances, spins, and element types into a contiguous, alignas(64) SoA buffer. In Phase B, the packed buffer is processed in 8-wide SVE2 batches (matching the 512-bit SVL for FP64). Three techniques are critical within the SVE2 kernel: (i) online Chebyshev recurrence—basis functions Tk+1 = 2xTk − Tk−1 are evaluated inside the vector register file via a running pair of SVE vectors, avoiding the sizelesstype array limitation of ARM SVE; (ii) predicated multitype dispatch—svsel_f64 selects per-lane coefficients for different element types (e.g. Fe vs. Ge) in a single instruction, eliminating branch divergence and type-sorted gather/scatter; (iii) gather-load coefficient inner products—strided ANN coefficients are loaded via svld1_gather and contracted with contiguous basis vectors in one fused multiply–horizontalreduce sequence. Together with a column-major-to-row-major transposition of the descriptor-derivative array g_Fp (placing each atom’s 66 Fp parameters in ∼8 contiguous cache lines),
A. Algorithmic Innovations 1) A spin-aware local descriptor pipeline: We extended the existing NEP inference path [9] to incorporate magnetic information while preserving the same locality, cutoff-based neighbor traversal, and regular evaluation pattern as in the structural descriptor pipeline. In practice, the local feature vector augments the original structural channels with three groups of magnetic channels that correspond to onsite, pairwise, and angular spin-lattice information. This organization allows the magnetic extension to reuse the same radial-basis and angularaccumulation infrastructure as the structural part, keeping the descriptor evaluation local and implementation-friendly. 2) Implementation-oriented descriptor organization: The magnetic channels are organized in three groups. The first group collects onsite information from the local spin state. The second group evaluates pairwise spin-bond couplings over the neighbor list using the same radial carrier as the structural interaction channels. The third group evaluates spin-weighted angular information by accumulating directional channels over neighbors and contracting them into rotationally invariant quantities. Additional mixed channels are formed through local contractions of these accumulated quantities. From an implementation perspective, all magnetic channels follow the same basic computational pattern: local neighbor traversal, channel-wise accumulation, and small dense contractions. As a result, the magnetic extension increases arithmetic work without introducing new global data dependencies or irregular communication paths. 3) Self-consistent midpoint spin update: For cases with strong feedback between the spin state and the effective field, we augment the explicit predictor-corrector update with an optional self-consistent midpoint iteration. Starting from the beginning-of-step spin state, the algorithm repeatedly forms a midpoint configuration, reevaluates the force and effective field at that midpoint, and reapplies the existing one-step update until either a convergence criterion or an iteration cap is reached. An accelerated fixed-point variant with regularization can be enabled for strongly nonlinear cases. This preserves the original operator structure while improving robustness under strong state dependence. Because the procedure may trigger
4
(a) NEP Atomic Environment Representation (a1) Spin-Lattice Representation
(b) Key Optimization
(a2) NEP inference pipeline
(b2) Data Pre-arrangement Coefficient Repacking
Description construction Lattice q𝑟
Lattice q𝑎
Before:
C0
After:
C0
C1
C2
C2
C4
C1
type 0
C3
type 1 Contiguous for FMOPA column load
Row-major Fp[N × D]
Column-major Fp[D × N]
𝒒𝒂
𝒒𝒔𝟏
𝒒𝒔𝟐
𝒒𝒔𝟑
d0
d1
d2
d3
Per-type ANN
a0
0
1
2
3
h = tanh (𝑊0 (t) ∙ q – 𝑏0 (t))
a1
4
5
6
a2
8
9
10
b2
∑
d0
d1
d2
d3
a0
0
1
2
3
7
a1
4
5
6
7
11
a2
8
9
10
11
transform
atom 0 → indices 0,1,2,3 contiguous → cache-friendly
atom 0 → indices 0,3,6,9 stride = N between dimensions
𝐸𝑖 = 𝑤1 (𝑡)𝑇 ∙ h –𝑏1
E F ω σ
(b4) Kernel Fusion + SME pipeline SME pipeline
Kernel Fusion 𝐸𝑡𝑜𝑡𝑎𝑙 = ∑𝐸𝑖 F = - ∂E/ ∂R
b4
8 neighbors/ batch
SMSTART
Atomic force
8 ZA titles
Basis
Magnetic torque
LLG Equation Coeff.
Newton Equation
b3
ω = - ∂E/ ∂S
…
C5
(b3) 𝐅𝐩 Array Layout
𝑞𝑠𝑐𝑎𝑙𝑒𝑑 = q ∙ 𝑞𝑠𝑐𝑎𝑙𝑒𝑑
𝒒𝒓
ZA Tile
Fused Force Kernel Streaming-in
Z0 Z1 Z2 Z3
type 0
type 1
ZA0
ZA2
ZA1
ZA3
ZA4
ZA6
ZA5
ZA7
8x8 FP64
(b1) SVE2 Pre-staging + Batch Vectorization Pre-Staging
Scattered
repack
b1
q = [q𝑟 | q𝑎 | q𝑠 ]
…
C3
Spin qs
SMSTORT Force /Virial/ Spin Torque
SVE2 register (512-bit)
Fig. 2: Fine-grained overview of the DFT-trained NEPSPIN framework for spin-lattice molecular dynamics. (a) Atomicenvironment representation and NEP inference workflow for spin-lattice interactions. The workflow in panel (a2) is partitioned into four highlighted stages, denoted b1–b4. (b) Performance-oriented implementation optimizations corresponding one-to-one to these four stages in (a2), indicating precisely where data pre-arrangement, array-layout transformation, SVE2 pre-staging and batch vectorization, and kernel fusion with SME execution are applied along the pipeline.
these optimizations deliver a combined 29.6% reduction in loop time over the fused kernel alone (21.74 s→15.30 s). 4) SME three-stage pipeline with predicate-driven type disambiguation: The most architecture-specific innovation is the reformulation of the dominant coefficient inner products as outer-product GEMM operations executed on the ARM SME matrix engine. The NEPSPIN force kernel is restructured into a three-phase pipeline operating on batches of 8 neighbors. In the preparation phase, scalar distance computation, cutoff filtering, and Chebyshev basis recurrence produce perneighbor basis vectors fn12 [k] and dfn12 [k], written into an SoA buffer with a [basis][batch] layout so that the batch dimension is contiguous for FMOPA row operands. Shared Chebyshev polynomials Tk , Uk are computed once per neighbor pair and reused across all four magnetic sub-terms (Exchange, DMI, ANI, SIA), eliminating three redundant recurrences (∼120 FLOP per pair). Precomputed products fp·dC and fp·Cv are factored out of the sub-term loops, reducing perk cost from 6 to 2 FLOP. The SME GEMM phase performs the coefficient–basis inner products that dominate the per-neighbor cost. A mixed-type neighbor batch poses a data-layout challenge: in a naı̈ve implementation, Fe and Ge neighbors must be separated by gather/scatter operations before their respective coefficient vectors can be contracted, adding 35% overhead. We eliminate
this overhead through predicate-driven type disambiguation. An __arm_locally_streaming function maps the eight available ZA Tiles (8 × 8 FP64 each) into four logical groups: ZA0–ZA1 and ZA2–ZA3 accumulate radial coefficient–basis products (Crad × dfn/fn) for Fe and Ge respectively, while ZA4–ZA5 and ZA6–ZA7 accumulate the corresponding spin products (Cspin × dfn/fn). Within each FMOPA instruction, the column predicate is set according to the element type of each lane, masking inactive lanes to zero so that Fe and Ge contributions accumulate into their designated tiles without data reshuffling. After the basis loop completes, the two pertype tiles in each group are reduced by element-wise addition. For a basis size of 8 (9 recurrence iterations), the entire GEMM phase issues 126 instructions per 8-neighbor batch, equivalent to 1,008 FMA operations. In the post-processing phase, each neighbor’s force and torque contributions are assembled from the GEMM results using the precomputed tables fp·dC/fp·Cv and the cached geometric/spin quantities, with results accumulated into threadprivate force buffers followed by a parallel reduction. The three-stage pipeline delivers an additional 11.2% looptime reduction (13.63 s→12.11 s) on top of the SVE2 optimizations. Across all five architecture-specific optimizations, the cumulative single-node speedup is 2.36× over the OpenMP baseline and 70.9× over the serial code (Fig. 5).
5
TABLE II: Per-node hardware specification.
LX2 Processor 0 NUMA Domain
Compute Die 0 NUMA Domain
NUMA Domain DDR Ctrl
DDR Ctrl
SDMA Engine
SDMA Engine
DDR Ctrl
DDR Ctrl
Compute Die 1 NUMA Domain
Attribute
Value
Processor ISA Extensions Cores per processor Physical cores per node HBM per processor DDR per processor NUMA domains per node DP peak per processor Network bandwidth per node
Armv9-based LX2 (×2 per node) Armv9-A SVE2, SME 304 (2 dies × 152) 608 32 GB (8 stacks, 4 TB/s) 256 GB (4 NUMA domains per die) 16 60.3 TFLOP/s (SME) 1.6 Tb/s (LingQi fat-tree)
(12.45 million cores), corresponding to full-scale machine allocation. B. Benchmark application and setup To measure performance, we used the real-temperature helix-to-skyrmion transformation in FeGe as a wholeapplication benchmark for coupled spin–lattice dynamics. This benchmark is particularly demanding because every stage of the simulation pipeline remains active throughout the transition, including communication, descriptor evaluation, force-and-torque inference, and coupled time integration. Our performance-optimized spin–lattice dynamics evolves coupled atomic and spin degrees of freedom for up to 1.34 × 1012 atoms and an equal number of spins in a FeGe system of size 3018×2414×2414 Å (3.02×2.41×2.41 µm). We first prepare a low-temperature helical state under experimentally relevant geometry and boundary conditions, and then drive the system through the temperature regime in which skyrmion seeds emerge, proliferate, and reorganize into extended topological textures. The measured workload therefore includes neighbor construction, halo exchange, spin-aware descriptor evaluation, unified force-and-torque inference, and long-time coupled time integration.
Fig. 3: Internal architecture of an LX2 processor, comprising two compute dies, on-package HBM stacks, SDMA engines, and four DDR-attached NUMA domains. Combined with MPI spatial decomposition across nodes and OpenMP thread parallelism within each rank, this three-level execution model enables the coupled spin–lattice application to scale to 20,480 nodes and 12.45 million CPU cores. 6. H OW P ERFORMANCE WAS M EASURED A. System architecture and hardware platform All measurements were performed on the LineShine supercomputer, a newly deployed exascale system unveiled in April 2026 by the National Supercomputing Center in Shenzhen (NSCC-SZ), China. Each node contains two Armv9-based LX2 processors; the internal architecture of a single processor is illustrated in Fig. 3. The LX2 integrates two compute dies (304 cores total) and eight on-package HBM stacks (32 GB, 4 TB/s aggregate bandwidth) in a single package. Each compute die contains 152 cores and is paired with 128 GB of offpackage DDR memory organized into four NUMA domains, for a total of 256 GB DDR per processor. Within each NUMA domain, cores share HBM; each compute die includes a dedicated SDMA engine for data movement between DDR and HBM. At the node level, the two LX2 processors provide 608 physical cores, 64 GB HBM, and 512 GB DDR, arranged in 16 NUMA domains. The memory subsystem has no shared last-level cache, making performance highly sensitive to data locality—a constraint addressed by the layout and pre-staging optimizations described in Section 5. The complete per-node hardware specification is summarized in Table II. The LX2 supports FP64, FP32, FP16, and INT8 through SME and SVE units, delivering up to 60.3 TFLOP/s FP64 per processor (120.6 TFLOP/s per node). Large-scale compute nodes are interconnected via the LingQi high-speed network with a dual-plane, multi-rail fat-tree topology, providing 1.6 Tb/s bandwidth per node. The entire system delivers over 2.5 EFLOP/s peak performance in FP64. Our experiments use up to 20,480 nodes
Fig. 4: Schematic illustration of the helical magnetic structure in FeGe. Arrows represent local spin orientations with a helical period of λ = 57.3 nm. Colors indicate the x-component of the spin Sx , ranging from blue (Sx = −1) to red (Sx = +1). Capturing this process requires substantially more than identifying equilibrium phases: the helix pitch must first be reproduced quantitatively, because it sets the characteristic magnetic length scale through the competition among exchange, Dzyaloshinskii–Moriya interaction, anisotropy, and
6
Definition
Loop time (s) atom-step/s/node Parallel eff. η Speedup Sustained FLOPS
wall-clock time of run loop Natom × Nstep /Loop time / Nnode T1 /TN (weak) or S/Sideal (strong) Tbase /TN modeled FLOPs per step × Nstep /Loop time
Cumul. speedup 3.0
20 10
0 ine ine ion ion out ure asel ce Fus torizat jor Lay struct E Pipel B r e P OM +Fo E2 Vec ow-Ma ular R +SM V g R S n + +Fp +A
2.5 2.0 1.5 1.0
Fig. 5: Ablation study of single-node optimizations. Each bar shows the loop time after cumulatively applying the listed optimization. Annotations indicate the incremental reduction from the previous step. Gray bar: OMP thread-parallel baseline; blue bars: architecture-specific optimizations.
TABLE III: Performance metrics used in this work. All are derived from the MD Loop time; Natom denotes the total atom count, Nstep the number of time steps, and subscripts 1 and N refer to single-node and N -node measurements. Metric
30
Cumulative speedup
Loop time (s)
Zeeman energy. As shown in Fig. 4, our model reproduces the zero-field helix pitch of FeGe over the temperature range of interest. Performance is characterized by five complementary metrics, summarized in Table III. All metrics are derived from a common timing primitive, the MD Loop time, defined as the wall-clock time spent in the time-stepping run loop, excluding one-time setup costs. The application-level floatingpoint throughput (FLOPS) is estimated from the dominant NEPSPIN workload in each MD step, including radial and angular pair interactions, ANN evaluation, spin-dependent terms, and optional short-range interactions; the modeled FLOPs per step are summed across all MPI ranks and divided by the Loop time. Different subsections of Section 7 select the metric most suited to each analysis: loop time for the ablation study, atom-step/s for cross-framework comparison, and parallel efficiency, speedup, and sustained FLOPS for scalability characterization.
walks and two redundant Chebyshev basis evaluations per atom; this reduces loop time by 23.9% (28.57 s→21.74 s). SVE2 vectorization (Step 2) restructures the fused kernel into a batch processing pipeline that collects valid neighbors in contiguous SoA buffers and processes them in 8-wide SVE2 batches with predicated type-aware coefficient selection, delivering an additional 23.3% reduction (21.74 s→16.67 s). The remaining three optimizations target data layout and instruction-level efficiency: Fp row-major layout (Step 3) transposes the descriptor-derivative array so that each atom’s 66 Fp parameters occupy contiguous cache lines, improving L1 hit rate (−8.2%); angular descriptor restructure (Step 4) hoists shared Chebyshev polynomials and precomputes fp·dC / fp · Cv products across Exchange/DMI/ANI/SIA sub-terms, cutting redundant FLOPs (−11.0%); The SME three-stage pipeline (Step 5) reformulates the coefficient inner products as outer-product GEMM operations using all eight ZA Tiles with predicate-driven type disambiguation, achieving a further 11.2% reduction. Cumulatively, the five architecture-specific optimizations achieve a 2.36× speedup over the OMP baseline (28.57 s→12.11 s), and the total speedup from serial is 70.9×. 2) Comparison with DeePMD: To contextualize NEPSPIN, we compare it against DeePMD [12], the ML-potential framework featured in a prior Gordon Bell Prize, in both accuracy and computational throughput. Accuracy: Table IV compares the prediction errors of NEPSPIN and DeePMD on the same FeGe spin–lattice validation set derived from spin-constrained DFT calculations. NEPSPIN achieves an energy RMSE of 1.85 meV/atom (vs. 1.69 for DeePMD), a force RMSE of 45.67 meV/Å (vs. 46.34), and a torque RMSE of 11.16 meV/µB (vs. 12.58). Both models reach comparable accuracy across all three quantities, confirming that the lightweight NEP architecture does not sacrifice physical fidelity relative to the deep neural-network potential. Throughput scaling with system size: To ensure a controlled comparison with the DeePMD-based Gordon Bell campaigns [11], [12], which all adopted liquid water as their
7. P ERFORMANCE R ESULTS A. Single-node performance Before presenting multi-node scaling results, we establish single-node efficiency through two complementary analyses: an ablation study validating each architecture-specific optimization, and a direct comparison with DeePMD [12] on identical hardware. DeePMD is chosen as the reference because it is the ML interatomic-potential framework featured in the 2020 Gordon Bell Prize and remains the most widely adopted baseline for large-scale ML-MD performance evaluation. We compare the two frameworks in both prediction accuracy and single-node computational throughput to demonstrate that NEPSPIN achieves competitive physical fidelity while delivering substantially higher efficiency on the target ARM architecture. 1) Optimization ablation study: Starting from a serial baseline of 858.04 s per loop (single-threaded NEPSPIN force evaluation at one node), OpenMP thread parallelism across 608 cores reduces the loop time to 28.57 s (30.0×). Then, five architecture-specific optimizations are applied cumulatively. Fig. 5 reports the measured loop time after each step. The first two optimizations yield the largest individual gains. Spin-radial force fusion (Step 1) merges three independent force kernels –radial, spin-lattice, and magnetic torque—into a single neighbor traversal, eliminating two redundant neighbor
7
TABLE IV: Accuracy comparison between NEPSPIN and DeePMD on the FeGe spin–lattice validation set. Values are placeholders pending final validation. Property Energy RMSE Force RMSE Magnetic Torque RMSE
Unit
NEPSPIN
DeePMD
meV/atom meV/Å meV/µB
1.85 45.67 11.16
1.69 46.34 12.58
B. Weak scaling In the weak-scaling study, the per-node workload is held constant while the number of nodes increases from 1 to 20,480. Two representative configurations are evaluated: a small case with 8.19 × 106 atoms per node and a large case with 6.55 × 107 atoms per node (8× larger). At full scale (20,480 nodes, 12.45 M cores), the small case evolves 1.68 × 1011 atoms (167.77 billion) and the large case evolves 1.34 × 1012 atoms (1.34 trillion). Figure 7(a) reports the parallel efficiency. At full machine scale (20,480 nodes), the small case retains 89.7% efficiency and the large case 85.3%. The small case is consistently more efficient, owing to better cache reuse within the percore L1/L2 hierarchy on the LX2 architecture. The primary source of efficiency loss at extreme scale is the increasing surface-to-volume ratio of each MPI sub-domain, which raises the ghost-atom fraction and inter-node communication volume. The high arithmetic intensity of the NEPSPIN angular kernels—O(Nneigh × L2max ) FLOPs per atom ver2/3 sus O(Nlocal ) halo exchange—amortizes this overhead and maintains high efficiency across the full scaling range. Figure 7(b) shows the sustained application-level FLOPS. Both configurations scale nearly linearly with node count, reaching 48.5 PFLOPS (small case) and 43.3 PFLOPS (large case) at 20,480 nodes. The single-node baseline is 2.64 TFLOPS (small) and 2.48 TFLOPS (large), both in double precision.
Speed per node (atom-step/s per node)
Scaling to larger simulation systems
107
1.00in7eT) h (LineS
679M0) 1.56sB.) 3.92B) 24.9oj.B) ea it'2 it'2 pr (Summ(Fugaku m (Summ (Fugaku
106
DeepMD, 2020 Summit (GPU, V100×6) DeepMD, 2022 Summit (GPU, V100×6) DeepMD, 2022 Fugaku measured DeepMD, 2022 Fugaku projected NEP, This work (LineShine)
105 108
109
1010
Number of atoms
1011
1012
Fig. 6: Per-node throughput (atom-step/s) of NEPSPIN on LineShine (1–20,480 nodes, water system) compared with DeepMD on Summit [12] and Fugaku [11]. Colored markers denote the maximum atom count of each campaign.
C. Strong scaling Strong-scaling behaviour is characterised with two fixedsize systems—67 B and 268 B atoms—spanning 2,048 to 20,480 nodes (Figure 8 and Table V). Both cases exhibit nearideal speedup at moderate concurrency: the first doubling of node count delivers 1.99× and 1.98× speedup, respectively. At full machine scale the 268 B-atom case retains 96.0% parallel efficiency (4.80× over 4,096 nodes), while the 67 Batom case reaches 89.6% (8.96× over 2,048 nodes) despite a 10× node-count range that reduces the per-node workload to only ∼3.3 M atoms. The efficiency gap between the two cases is consistent with the scaling characteristics of shortrange force models: because the NEPSPIN force evaluation scales as O(Nlocal · Nneigh ) while ghost-atom exchange grows 2/3 as O(Nlocal ), smaller sub-domains expose a larger communication fraction. Nevertheless, the architecture-aware optimizations described in Section 5 keep efficiency above 89% across the entire scaling range.
benchmark system, we reproduce the same water simulation on LineShine using our NEPSPIN engine with 49,152,000 atoms per node. This choice isolates differences in framework efficiency and hardware platform from those in model complexity or physical system properties. The total throughput scales from 4.38 × 106 atom-step/s on a single node to 6.77 × 1010 atomstep/s on 20,480 nodes, representing a 1.55×104 -fold increase and corresponding to approximately 1.0 × 1012 atoms simulated simultaneously. Compared to the DeepMD results on Fugaku [11], which achieved a maximum measured throughput of 3.29 × 108 atom-step/s at 1/16 of the machine (∼9,874 nodes, 1.56 × 109 atoms), our implementation on LineShine surpasses this by more than two orders of magnitude in absolute throughput while extending the accessible system size from 1.6 × 109 to over 1.0 × 1012 atoms. As illustrated by the milestone markers in Figure 6, each successive Gordon Bell campaign has pushed the atom-count frontier—from 679 M (Summit 2020) to 3.9 B (Summit 2022) and 1.56 B measured on Fugaku 2022—while LineShine now reaches 1.0 T atoms. On a per-node basis, LineShine delivers ∼130× the throughput of Fugaku and ∼10× that of Summit 2022, despite providing only ∼2.8× the FP64 peak per node. This disproportionate advantage reflects the lower arithmetic intensity of the NEP descriptor relative to DeepMD, combined with the LX2 node’s efficient utilization of SME/SVE vector units and on-package HBM for memorybound neighbor-list operations.
D. Sustained floating-point performance At the full scale of 20,480 nodes (12.45 M cores), the application sustains 48.5 PFLOPS (small case) and 43.3 PFLOPS (large case) in double precision, as reported in Figure 7(b). The near-linear growth of sustained FLOPS with node count confirms that the computation-to-communication ratio remains favorable across the entire scaling range. At single-node level, the small case achieves 2.64 TFLOPS and the large case 2.48 TFLOPS per node; at 20,480 nodes, these per node throughputs decrease to 2.37 TFLOPS and 2.11 TFLOPS
8
0.8
100.0 P 1.00 1.00 1.000.990.990.990.980.980.970.97 0.94 0.95 0.940.90 0.90 0.91 1.00 0.99 0.99 0.99 1.00 0.99 0.99 0.94 0.91 0.95 0.92 0.900.880.87 0.86 0.85
0.6
10.0 P
0.4 0.2
100.0 T 10.0 T
Ideal Small case (8M atoms/node) Large case (65M atoms/node)
1.0 T
1 2 4 8 16 32 64 12 8 25 6 5 1,012 2,024 4,0486 8,196 16 92 ,3 20 84 ,48 0
0.0
1.0 P
(b)
39.3 P 48.5 P 43.3 P 19.5 P 34.8 P 10.1 P 17.7 P 5.1 P 8.9 P 2.6 P 4.6 P 1.3 P 2.3 P 651.8 T 1.2 P 330.4 T 595.7 T 165.4 T 301.2 T 83.2 T 156.4 T 41.8 T 78.8 T 21.0 T 39.4 T 10.5 T 19.7 T 5.3 T 9.8 T 2.6 T 4.9 T 2.5 T
Ideal Small case (8M atoms/node) Large case (65M atoms/node)
1 2 4 8 16 32 64 12 8 25 6 5 1,012 2,024 4,0486 8, 96 16 192 ,3 20 84 ,48 0
Efficiency
1.0
(a)
FLOPS
1.2
Number of nodes
Number of nodes
Fig. 7: Weak-scaling results from 1 to 20,480 nodes. (a) Parallel efficiency for the small case (8 M atoms/node) and large case (65 M atoms/node). (b) Sustained application-level FLOPS for both configurations. Dashed lines indicate ideal scaling.
128
117.04 s
64 51.10 s
59.18 s 25.67 s
16
13.21 s
2,048
8,192
Ideal Small case (67B atoms fixed) Large case (268B atoms fixed)
8
Number of nodes
Case
Nodes
Atoms/node
67 B
2,048 4,096 8,192 16,384 20,480
268 B
4,096 8,192 16,384 20,480
30.42 s 24.38 s 6.80 s
5.71 s
16,38 4 20,48 0
32
4,096
Loop time (s)
TABLE V: Strong-scaling results for two fixed-size systems. Speedup and efficiency are relative to the smallest node count used for each case (2,048 for the small case; 4,096 for the large case).
Fig. 8: Strong-scaling loop time for two fixed-size systems. Green triangles: 6.71 × 1010 atoms (67.1 billion), scaled from 2,048 to 20,480 nodes (10× range). Blue squares: 2.68×1011 atoms (268.4 billion), scaled from 4,096 to 20,480 nodes (5× range). Dashed lines indicate ideal linear speedup from each respective baseline node count.
Loop time (s)
Speedup
Eff. (%)
32.77 M 16.38 M 8.19 M 4.10 M 3.28 M
51.10 25.67 13.21 6.80 5.71
1.00× 1.99× 3.87× 7.52× 8.96×
100.0 99.5 96.7 93.9 89.6
65.54 M 32.77 M 16.38 M 13.11 M
117.04 59.18 30.42 24.38
1.00× 1.98× 3.85× 4.80×
100.0 98.9 96.2 96.0
120.6 TFLOP/s × 20,480 ≈ 2.47 EFLOP/s), the sustained 48.5 PFLOPS corresponds to approximately 2.0% of the FP64 peak. Low peak utilization is, in fact, characteristic of particle-based and neighbor-list-driven applications on modern hardware. Fedeli et al. [30] report that electromagnetic particle-in-cell (PIC) codes—whose dominant kernels (current deposition and field gathering) share the same gather/scatter memory-access pattern as MLIP neighbor-list operations— sustain only 1–3% of the DP peak on Fugaku A64FX CPUs, with the HPCG benchmark proposed as a more representative roofline reference than HPL. They note that even the highly optimized VPIC code achieved only 13.2% utilization on Roadrunner (2008), and that 2–10% is the typical range for production PIC codes across architectures. Our ∼2% utilization on the ARM-based LineShine is consistent with this picture: the NEPSPIN kernels are dominated by neighborlist traversal, descriptor accumulation via irregular gather/scatter, and shallow-network GEMV, all of which are memorybandwidth-bound and cannot saturate the SME matrix engine designed for dense outer-product GEMM.
respectively, consistent with the 89.7% and 85.3% parallel efficiencies reported above. Performance is limited primarily by the gradual increase in communication overhead rather than by load imbalance or memory bandwidth saturation, indicating that architecture-aware kernel design and SoA data layout effectively exploit the SVE2 and SME capabilities of the underlying hardware. All results are obtained in double precision because accurate spin–lattice dynamics—particularly the spinorbit-mediated torque terms—requires full FP64 fidelity throughout the force and torque evaluation pipeline. At the deployed scale of 20,480 nodes (theoretical peak
9
Fig. 9: Representative subregion extracted from the full-machine extreme-scale spin–lattice simulation, illustrating the spatial evolution from helical order to skyrmionic textures in a magnetic stripe. The displayed region has a size of 1.141 µm × 169.7 nm × 1.8 nm and is taken from the full production run under the same simulation protocol. The system is initialized from a random spin distribution and evolved at T = 160 K under a magnetic-field gradient along the x direction, with field values of 0, 0.05, 0.1, and 0.15 T. The upper panel (a) and (b) shows the full texture distribution at T = 160 K and afterrelaxation, (b) The lower-left panel enlarges a helical region, and the side view beneath it illustrates the spin rotation within the helix, with the white dashed lines tracing the helical winding. The lower-center panel (d) enlarges a skyrmion region, where three skyrmions exhibit a lattice-like arrangement. The lower-right panel (e) enlarges a area where a locally disordered spin and atom configuration develops as the helical stripe breaks. The figure highlights the ability of the present framework to capture both global texture reorganization and local topological spin structure at real temperature. By comparison, DeepMD-kit reports ∼22% utilization on both Fugaku and Summit [11], owing to its deep embedding networks whose inference is dominated by dense matrix multiplications with high arithmetic intensity. However, the Fugaku figure of 119 PFLOPS is a projected estimate extrapolated from runs on 1/16 of the machine, not a full-scale measurement. NEP, by contrast, adopts a compact Chebyshev polynomial descriptor with O(103 ) trainable parameters— deliberately trading arithmetic density for per-atom efficiency. This design delivers ∼10–130× higher per-node throughput than DeepMD (Section 7-A2) at comparable accuracy, making time-to-solution the more relevant metric for this class of applications.
texture landscape in the stripe (Fig. 9), ranging from multihelix textures to isolated skyrmions and skyrmion-latticelike regions. The enlarged local snapshot further shows that skyrmion nucleation proceeds through a localized rupture of the helical texture, accompanied by the formation of a strongly twisted and partially disordered intermediate region, which subsequently evolves into a topological seed. Importantly, this helix-breaking event is observed only when finite temperature is applied to both the lattice and spin degrees of freedom together with the external magnetic field. In contrast, under the same field but without thermal activation, the helical texture remains intact within the simulation time window and no skyrmion nucleation is observed. This indicates that the magnetic field alone is insufficient to overcome the topological and energetic barrier associated with helix breaking, whereas thermal fluctuations in the coupled spin–lattice system are essential for activating the transition pathway. More broadly, this work demonstrates that extreme-scale atomistic spin–lattice simulation can serve as a predictive tool for linking atomic structure, magnetic order, microstructure, and functional response in magnetic materials. This capability is directly relevant to topological-spin-texture-based spintronic devices, to microstructure- and defect-controlled magnetic performance, and to the broader class of advanced functional materials governed by coupled order parameters.
8. I MPLICATIONS A flagship application enabled by the present framework is the real-temperature transformation from helical order to skyrmion textures in FeGe. At the same time, the micrometer-scale simulation domain accommodates many helical and skyrmionic periods within a single atomistically resolved calculation, providing access to competing boundaries, fluctuation-mediated nucleation channels, and texture reorganization pathways. Under the thermal and field protocol described in the text, the system develops a spatially heterogeneous topological
10
ACKNOWLEDGMENT
[15] H. Yu, B. Liu, Y. Zhong, L. Hong, J. Ji, C. Xu, X. Gong, and H. Xiang, “Physics-informed time-reversal equivariant neural network potential for magnetic materials,” Physical Review B, vol. 110, p. 104427, 2024. [Online]. Available: https://doi.org/10.1103/PhysRevB.110.104427 [16] T. Yang, Z. Cai, Z. Huang, W. Tang, R. Shi, A. Godfrey, H. Liu, Y. Lin, C.-W. Nan, M. Ye, L. Zhang, K. Wang, H. Wang, and B. Xu, “Deep learning illuminates spin and lattice interaction in magnetic materials,” Physical Review B, vol. 110, p. 064427, 2024. [Online]. Available: https://doi.org/10.1103/PhysRevB.110.064427 [17] O. Chubykalo-Fesenko, U. Nowak, R. W. Chantrell, and D. Garanin, “Dynamic approach for micromagnetics close to the curie temperature,” Physical Review B, vol. 74, no. 9, p. 094436, 2006. [Online]. Available: https://doi.org/10.1103/PhysRevB.74.094436 [18] R. F. L. Evans, W. J. Fan, P. Chureemart, T. A. Ostler, M. O. A. Ellis, and R. W. Chantrell, “Atomistic spin model simulations of magnetic nanomaterials,” Journal of Physics: Condensed Matter, vol. 26, no. 10, p. 103202, 2014. [Online]. Available: https://doi.org/10.1088/09538984/26/10/103202 [19] J. C. Tully, “Molecular dynamics with electronic transitions,” The Journal of Chemical Physics, vol. 93, no. 2, pp. 1061–1071, 1990. [Online]. Available: https://doi.org/10.1063/1.459170 [20] R. Crespo-Otero and M. Barbatti, “Recent advances and perspectives on nonadiabatic mixed quantum-classical dynamics,” Chemical Reviews, vol. 118, no. 15, pp. 7026–7068, 2018. [Online]. Available: https://doi.org/10.1021/acs.chemrev.7b00577 [21] X. He, N. Helbig, M. J. Verstraete, and E. Bousquet, “Tb2j: A python package for computing magnetic interaction parameters,” Computer Physics Communications, vol. 264, p. 107938, 2021. [Online]. Available: https://doi.org/10.1016/j.cpc.2021.107938 [22] M. Rinaldi, M. Mrovec, A. Bochkarev, Y. Lysogorskiy, and R. Drautz, “Non-collinear magnetic atomic cluster expansion for iron,” npj Computational Materials, vol. 10, p. 12, 2024. [Online]. Available: https://doi.org/10.1038/s41524-024-01196-8 [23] M. Y. Lavrentiev, D. Nguyen-Manh, and S. L. Dudarev, “Magnetic cluster expansion model for bcc-fcc transitions in fe and fe-cr alloys,” Physical Review B, vol. 81, no. 18, p. 184202, 2010. [Online]. Available: https://doi.org/10.1103/PhysRevB.81.184202 [24] J. Tranchida, S. J. Plimpton, P. Thibaudeau, and A. P. Thompson, “Massively parallel symplectic algorithm for coupled magnetic spin dynamics and molecular dynamics,” Journal of Computational Physics, vol. 372, pp. 406–425, 2018. [Online]. Available: https://doi.org/10.1016/j.jcp.2018.06.042 [25] J. B. J. Chapman and P.-W. Ma, “A machine-learned spin-lattice potential for dynamic simulations of defective magnetic iron,” Scientific Reports, vol. 12, p. 22451, 2022. [Online]. Available: https://doi.org/10.1038/s41598-022-25682-5 [26] H. Yu, Y. Zhong, L. Hong, C. Xu, W. Ren, X. Gong, and H. Xiang, “Spin-dependent graph neural network potential for magnetic materials,” Physical Review B, vol. 109, p. 144426, 2024. [Online]. Available: https://doi.org/10.1103/PhysRevB.109.144426 [27] I. Novikov, B. Grabowski, F. Körmann, and A. Shapeev, “Magnetic moment tensor potentials for collinear spin-polarized materials reproduce different magnetic states of bcc fe,” npj Computational Materials, vol. 8, p. 13, 2022. [Online]. Available: https://doi.org/10.1038/s41524-022-00696-9 [28] A. S. Kotykhov, K. Gubaev, M. Hodapp, C. Tantardini, A. V. Shapeev, and I. S. Novikov, “Constrained dft-based magnetic machinelearning potentials for magnetic alloys: a case study of fe–al,” Scientific Reports, vol. 13, p. 19728, 2023. [Online]. Available: https://doi.org/10.1038/s41598-023-46951-x [29] Z. Fan, Z. Zeng, C. Zhang, Y. Wang, K. Song, H. Dong, Y. Chen, and T. Ala-Nissila, “Neuroevolution machine learning potentials: Combining high accuracy and low cost in atomistic simulations and application to heat transport,” Physical Review B, vol. 104, p. 104309, 2021. [Online]. Available: https://doi.org/10.1103/PhysRevB.104.104309 [30] L. Fedeli, A. Huebl, F. Boillod-Cerneux, T. Clark, K. Gott, C. Hillairet, S. Jaure, A. Leblanc, R. Lehe, A. Myers, C. Piber, O. Shapoval, R. Sandberg, M. Zaim, J.-L. Vay, and H. Vincenti, “Pushing the frontier in the design of laser-based electron accelerators with groundbreaking mesh-refined particle-incell simulations on exascale-class supercomputers,” in Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis (SC ’22), 2022. [Online]. Available: https://doi.org/10.1109/SC41404.2022.00008
This work was supported in part by the National Key R&D Program of China under Grant No. 2025YFB3003603. R EFERENCES [1] A. Fert, “Nobel lecture: Origin, development, and future of spintronics,” Reviews of modern physics, vol. 80, no. 4, pp. 1517–1530, 2008. [Online]. Available: https://doi.org/10.1103/RevModPhys.80.1517 [2] M. Z. Hasan and C. L. Kane, “Colloquium: topological insulators,” Reviews of modern physics, vol. 82, no. 4, pp. 3045–3067, 2010. [Online]. Available: https://doi.org/10.1103/RevModPhys.82.3045 [3] X. Yu, N. Kanazawa, Y. Onose, K. Kimoto, W. Zhang, S. Ishiwata, Y. Matsui, and Y. Tokura, “Near room-temperature formation of a skyrmion crystal in thin-films of the helimagnet FeGe,” Nature materials, vol. 10, no. 2, pp. 106–109, 2011. [Online]. Available: https://doi.org/10.1038/nmat2916 [4] P. Bak and M. H. Jensen, “Theory of helical magnetic structures and phase transitions in MnSi and FeGe,” Journal of Physics C: Solid State Physics, vol. 13, no. 31, pp. L881–L885, 1980. [Online]. Available: https://doi.org/10.1088/0022-3719/13/31/002 [5] N. Nagaosa and Y. Tokura, “Topological properties and dynamics of magnetic skyrmions,” Nature nanotechnology, vol. 8, no. 12, pp. 899– 911, 2013. [Online]. Available: https://doi.org/10.1038/nnano.2013.243 [6] X. Z. Yu, N. Kanazawa, Y. Onose, K. Kimoto, W. Z. Zhang, S. Ishiwata, Y. Matsui, and Y. Tokura, “Near room-temperature formation of a skyrmion crystal in thin-films of the helimagnet FeGe,” Nature Materials, vol. 10, no. 2, pp. 106–109, 2011. [Online]. Available: https://doi.org/10.1038/nmat2916 [7] L. Xu, H. Han, J. Fan, D. Shi, D. Hu, H. Du, L. Zhang, Y. Zhang, and H. Yang, “Magnetic entropy change and accurate determination of Curie temperature in single-crystalline helimagnet FeGe,” EPL (Europhysics Letters), vol. 117, no. 4, p. 47004, 2017. [Online]. Available: https://doi.org/10.1209/0295-5075/117/47004 [8] X. Yu, Y. Onose, N. Kanazawa, J. H. Park, J. Han, Y. Matsui, N. Nagaosa, and Y. Tokura, “Real-space observation of a twodimensional skyrmion crystal,” Nature, vol. 465, no. 7300, pp. 901–904, 2010. [Online]. Available: https://doi.org/10.1038/nature09124 [9] Z. Fan, Z. Zeng, C. Zhang, Y. Wang, K. Song, H. Dong, Y. Chen, and T. Ala-Nissila, “Neuroevolution machine learning potentials: Combining high accuracy and low cost in atomistic simulations and application to heat transport,” Physical Review B, vol. 104, no. 10, p. 104309, 2021. [Online]. Available: https://doi.org/10.1103/PhysRevB.104.104309 [10] D. Zheng, X. Peng, Y. Huang, Y. Wang, D. Zhang, Z. Huang, Z. Cai, L. Zhang, M. Chen, B. Xu et al., “Integrating deep-learningbased magnetic model and non-collinear spin-constrained method: methodology, implementation and application,” npj Computational Materials, 2026. [Online]. Available: https://doi.org/10.1038/s41524025-01923-9 [11] Z. Guo, D. Lu, Y. Yan, S. Hu, R. Liu, G. Tan, N. Sun, W. Jiang, L. Liu, Y. Chen et al., “Extending the limit of molecular dynamics with ab initio accuracy to 10 billion atoms,” in Proceedings of the 27th ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming, 2022, pp. 205–218. [Online]. Available: https://doi.org/10.1145/3503221.3508425 [12] W. Jia, H. Wang, M. Chen, D. Lu, L. Lin, R. Car, E. Weinan, and L. Zhang, “Pushing the limit of molecular dynamics with ab initio accuracy to 100 million atoms with machine learning,” in SC20: International conference for high performance computing, networking, storage and analysis. IEEE, 2020, pp. 1–14. [Online]. Available: https://doi.org/10.1109/SC41405.2020.00009 [13] T. M. Razakh, T. Linker, Y. Luo, N. Piroozan, J. Pennycook, N. Kumar, A. Musaelian, A. Johansson, B. Kozinsky, R. K. Kalia, P. Vashishta, F. Shimojo, S. Hattori, K.-i. Nomura, and A. Nakano, “Multiscale light-matter dynamics in quantum materials: from electrons to topological superlattices,” in Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis (SC ’25), 2025, pp. 36–47. [Online]. Available: https://doi.org/10.1145/3712285.3771785 [14] P.-W. Ma, C. H. Woo, and S. L. Dudarev, “Large-scale simulation of the spin-lattice dynamics in ferromagnetic iron,” Physical Review B, vol. 78, no. 2, p. 024434, 2008. [Online]. Available: https://doi.org/10.1103/PhysRevB.78.024434
11