[5/6] Add variable-magnitude GLSD spin dynamics

CMakeLists.txt

@@ -22,6 +22,7 @@ file(GLOB GPUMD_SOURCES
     src/minimize/*.cu
     src/phonon/*.cu
     src/integrate/*.cu
+    src/integrate/spin/*.cu
     src/mc/*.cu
     src/force/*.cu
     src/force/*/*.cu
@@ -30,7 +31,7 @@ file(GLOB GPUMD_SOURCES
     src/utilities/*.cu)
 add_executable(gpumd ${GPUMD_SOURCES})
 target_include_directories(gpumd PRIVATE src)
-target_link_libraries(gpumd PRIVATE CUDA::cublas CUDA::cusolver CUDA::cufft)
+target_link_libraries(gpumd PRIVATE CUDA::cublas CUDA::cusolver CUDA::cufft Threads::Threads)
 
 # --- nep ---
 file(GLOB NEP_SOURCES

developer/README.md

@@ -52,7 +52,12 @@ The main source directories are:
 | `src/utilities/` | Shared utilities, GPU wrappers, input helpers, and containers |
 | `src/main_nep/` | `nep` executable |
 
-Files added under an existing source directory are normally collected automatically by the makefiles. For `nep`, the CUDA and HIP makefiles also collect `.cu` files one directory below `src/main_nep/`.
+Files added under an existing source directory are normally collected automatically by the makefiles. The CUDA and HIP makefiles also collect `.cu` files one directory below `src/force/` for `gpumd`, and `src/main_nep/` for `nep`.
+
+Spin inference, multi-GPU evaluation, external magnetic fields, and Spin-specific
+small-box kernels live in `src/force/nep_spin/`. Training in
+`src/main_nep/nep_spin/` includes the same magnetic descriptor and force headers;
+ordinary NEP keeps its own small-box kernels in `src/force/nep_small_box.cuh`.
 
 ## 3. Main `gpumd` architecture
 
@@ -199,6 +204,7 @@ The `nep` executable is implemented under `src/main_nep/` and shares utilities w
 | Top level | Entry point, model selection, optimizer, scheduling, and shared data structures |
 | `common/` | Reusable NEP parameter, data, loss, I/O, and runtime-compilation helpers |
 | `nep/` | Ordinary and temperature-dependent NEP |
+| `nep_spin/` | Spin3 parameters, magnetic labels/neighbors, response supervision, and model I/O |
 | `nep_charge/` | NEP charge modes 1 and 2 |
 | `nep_vdw/` | NEP with long-range dispersion |
 | `nep_charge_vdw/` | NEP with charge and dispersion |
@@ -233,6 +239,16 @@ Shared structure and dataset code handles common data. Model-specific labels, GP
 
 `evaluate()` receives contiguous candidate parameter vectors on the host. The model owns intermediate losses; `accumulate_fitness()` adds their contributions after common regularization. Keep model-specific reduction and transfer decisions behind these interfaces, so later GPU reduction or population parallelism does not require model branches in `Fitness` or `SNES`.
 
+Spin3 keeps `Structure` and `Dataset` unchanged and stores magnetic labels and device
+buffers in `Spin_Data`, indexed by `NEP_Batch_Data`. Labels are matched through
+`Structure::index_in_file` after batch sorting or energy-difference grouping.
+Its reader supplies the largest structural/magnetic cutoff to the shared geometry
+reader. `response.xyz` has separate device data and contributes only response loss;
+it does not fit the energy reference or scalers. Spin3 uses the shared
+`energy_reference` and block model/restart formats while preserving its
+fixed-reference scaler and the common generation-zero lifecycle. Legacy Spin3
+model and two-column restart files remain readable.
+
 ### Adding a training model
 
 1. Implement the `Potential` interface in the model's subdirectory, with its parameter parsing, labels, workspaces, losses, and output handling.
@@ -261,3 +277,17 @@ The basic internal units are:
 - charge: elementary charge
 
 User-facing changes to keywords, file formats, units, or physical definitions should be reflected in the documentation under `doc/`.
+
+Spin moments are in muB and magnetic forces in eV/muB. See the
+[model input units](../doc/gpumd/input_files/model_xyz.rst) and
+[external magnetic field](../doc/gpumd/input_parameters/add_mfield.rst)
+documentation for the input and Zeeman conventions.
+
+See [Spin regression checks](spin_validation.md) for portable format, training,
+and GPU validation commands.
+
+Spin lattice coupling and restart handling live in `src/integrate/spin/`.
+`Spin_Lattice_NHC` and `Spin_Lattice_MTTK` reuse the upstream mechanical
+integrators and add endpoint persistence and resumed temperature/pressure
+protocols only for Spin ensembles. Ordinary NHC/MTTK classes remain unchanged.
+The XYZ reader and restart writer delegate Spin metadata to this directory.

developer/spin_validation.md

@@ -0,0 +1,115 @@
+# Spin regression checks
+
+Run from the repository root with Python, NumPy, pytest, a C++ compiler and,
+for runtime checks, CUDA and an allocated NVIDIA GPU. Build `src/gpumd` and
+`src/nep` for that GPU; `GPUMD_COMMAND` and `NEP_COMMAND` can select other builds.
+Set `RESULTS` to an existing absolute directory outside the checkout.
+
+## Routine checks
+
+```sh
+python -m pytest --noconftest -q -m 'not slow' tests_pytest/test_spin_*.py
+python tests_pytest/validate_nep_spin3_runtime.py --block-format
+python tests_pytest/validate_nep_spin3_runtime.py --fixtures o3c2_uniform --cases chiral_soc
+python tests_pytest/validate_mfield_runtime.py
+python tests_pytest/validate_main_nep_spin3_training.py
+python tests_pytest/validate_spin_training_lifecycle.py
+```
+
+`--noconftest` avoids the ordinary NEP suite's ASE/Calorine fixtures. The CPU codec
+checks new and legacy model formats. Runtime checks cover independent energy,
+force, magnetic force and virial oracles, integration formulas, noise scaling,
+Spin lattice restart continuity. TSPIN
+and SIB each have one pytest implementation; the external-field SIB oracle is
+part of that suite. Training checks cover actual errors, defaults, masks,
+normalization, model export and restart. The main parity matrix has five cases:
+C2 with selective spins, C3 with all spins and nontrivial restart layout, C1 with
+typewise cutoff and torque supervision, ZBL, and a small periodic box with a wider
+magnetic cutoff. Lifecycle checks retain both full and subset spin environments.
+Detailed malformed model headers belong to the CPU codec suite; the CLI keeps a
+representative rejection check. Spin rejects flexible ZBL in training and in
+both legacy and block models; universal ZBL remains in the main parity matrix. The ordinary `tests_regression` suite
+remains a separate regression gate.
+
+The six files in `fixtures/nep_spin3` are two models, one structure and three
+frozen numerical oracles. The typewise-cutoff case derives its header from the
+C=2 model. SHA256 checks bind each oracle to its model. Never regenerate expected
+values with the candidate implementation to make a failing test pass.
+
+## Extended GPU checks
+
+Run these when changing the corresponding numerical kernels, decomposition,
+JIT or checkpoint contracts, and before merging those changes. The multi-GPU
+scripts require four allocated GPUs visible through `CUDA_VISIBLE_DEVICES`.
+
+```sh
+python -m pytest --noconftest -q -m slow tests_pytest/test_spin_*.py
+python tests_pytest/validate_main_nep_spin3_training.py --check-jit
+python tests_pytest/validate_spin_training_multigpu.py --output "$RESULTS/training"
+python tests_pytest/validate_spin_multigpu.py --binary src/gpumd \
+  --model tests_pytest/fixtures/nep_spin3/nep4_spin3_o3c2_uniform.nep \
+  --block-format --output "$RESULTS/runtime"
+python tests_pytest/validate_spin_glsd_runtime.py --binary src/gpumd \
+  --output "$RESULTS/glsd"
+python -m pytest --noconftest -q tests_pytest/test_spin_lattice_restart.py
+```
+
+The slow pytest cases check TSPIN temperature statistics and restart trajectories
+with the bundled NEP model; routine restart checks use a zero-force analytic
+model. Run lattice restart pytest on one and four GPUs to cover both execution
+paths. The GLSD runtime suite uses the bundled model by default and also checks
+analytic trajectories, timestep convergence, diffusion and the restart stream.
+No private Fe model or pre-refactor executable is required. Training multi-GPU
+checks use multiple structures per batch and independently recomputed energy
+errors. Optional historical binary comparisons are development evidence, not
+prerequisites for these regressions.
+
+## Focused compiled probes
+
+Build against current sources and objects. These commands assume V100; change
+`-arch=sm_70` for the allocated GPU. Rebuild the executable objects after source
+changes. The small probes check contracts not exposed by the command-line output.
+
+```sh
+c++ -std=c++14 -Isrc tests_pytest/rotation_response.cpp -o "$RESULTS/rotation_response"
+"$RESULTS/rotation_response"
+objects=$(find src/main_nep src/utilities -name '*.cu' ! -path 'src/main_nep/main.cu' \
+  | sed 's/\.cu$/.o/')
+for probe in spin_mforce_split restart_parity response_fitness; do
+  nvcc -std=c++14 -O2 -arch=sm_70 -Isrc "tests_pytest/$probe.cu" $objects \
+    -lcublas -lcusolver -lcufft -ldl -lpthread -o "$RESULTS/$probe"
+done
+python tests_pytest/validate_spin_mforce_split.py \
+  --oracle "$RESULTS/spin_mforce_split" --output "$RESULTS/split"
+RESTART_PROBE="$RESULTS/restart_parity" python tests_pytest/validate_restart_scaler.py
+RESPONSE_FITNESS_PROBE="$RESULTS/response_fitness" \
+  python tests_pytest/validate_response_file_runtime.py
+nvcc -std=c++14 -O3 -arch=sm_70 -Isrc tests_pytest/spin_glsd_analytic.cu \
+  src/integrate/spin_glsd.cu -o "$RESULTS/spin_glsd_analytic"
+python tests_pytest/validate_spin_glsd_analytic.py "$RESULTS/spin_glsd_analytic"
+```
+
+The magnetic-loss probe checks independent projections, including zero spins.
+The restart probe distinguishes the SNES mean from the exported elite. The
+response probes check rotation geometry, loss normalization and data isolation.
+Initialization and variable-type assignment are checked by the restart probe;
+torque and response-loss edge cases share the magnetic and response probes.
+GLSD's harmonic probe supplies independent statistical expectations
+that a bounded ANN potential cannot supply.
+
+For integration temperature changes, link `integrate_temperature.cu` against the
+GPUMD objects excluding `main_gpumd/main.o`; it checks the fixed initial cache and
+per-step temperature ramps:
+
+```sh
+objects=$(find src/main_gpumd src/minimize src/phonon src/integrate src/mc \
+  src/force src/measure src/model src/utilities -name '*.cu' \
+  ! -path 'src/main_gpumd/main.cu' | sed 's/\.cu$/.o/')
+nvcc -std=c++14 -O2 -arch=sm_70 -Isrc tests_pytest/integrate_temperature.cu \
+  $objects -lcublas -lcusolver -lcufft -lpthread -o "$RESULTS/integrate_temperature"
+"$RESULTS/integrate_temperature"
+```
+
+Keep timing, source/binary hashes and run logs outside the source tree. These
+checks establish software contracts, not magnetic phase behavior or predictive
+accuracy on unseen structures.

doc/gpumd/input_files/model_xyz.rst

@@ -51,6 +51,8 @@ Essentially any keyword is allowed, but we only read the following ones:
   * ``mass:R:1`` mass (*Optional*: default mass values will be used when this is missing)
   * ``charge:R:1`` charge (*Optional*: default charge is zero)  
   * ``vel:R:3`` velocity vector (*Optional*)
+  * ``spin:R:3`` magnetic-moment vector in Bohr magnetons (muB) (*Required for NEP_Spin*)
+  * ``spin_vel:R:3`` spin velocity vector (*Optional; requires ``spin:R:3``*)
   * ``group:I:number_of_grouping_methods`` grouping methods (*Optional*)
 
 
@@ -67,6 +69,15 @@ The mass should be given in units of the unified atomic mass unit (amu).
 The charge is in units of e (proton charge).
 The cell dimensions and atom coordinates should be given in units of Ångstrom. 
 Velocities need to be specified in units of Å/fs.
+Spin is specified in Bohr magnetons (muB). Spin velocity is in muB/fs.
+Magnetic force is the negative energy derivative with respect to spin, in eV/muB.
+The spin torque, spin crossed with magnetic force, is in eV.
+External magnetic induction B is specified in tesla (T).
+Spin vectors store Cartesian magnetic moments, not normalized directions.
+Readers do not normalize them or rescale model parameters; training and
+runtime data must use the same spin convention.
+Property names may not be repeated. NEP_Spin fails if ``spin:R:3`` is absent;
+missing spin values are never filled with zeros.
 
 Example
 -------

doc/gpumd/input_files/run_in.rst

@@ -101,6 +101,10 @@ Simulation setup
      - Apply an electric field to selected atoms
      - No
      - No
+   * - :ref:`add_mfield <kw_add_mfield>`
+     - Apply a magnetic field to an atom group or fractional region
+     - No
+     - Yes
    * - :ref:`add_spring <kw_add_spring>`
      - Add spring interactions to selected atom groups
      - No

doc/gpumd/input_parameters/add_mfield.rst

@@ -0,0 +1,59 @@
+.. _kw_add_mfield:
+
+:attr:`add_mfield`
+==================
+
+Apply a constant external magnetic induction to an atom group or a fractional-coordinate region::
+
+  add_mfield <group_method> <group_id> <Bx> <By> <Bz>
+  add_mfield region <amin> <amax> <bmin> <bmax> <cmin> <cmax> <Bx> <By> <Bz>
+
+Grouping method and group ID select an existing model.xyz group, as in
+add_efield. The three finite field components are in tesla (T). There is no
+separate magnitude or direction normalization. For a 10 T field along z::
+
+  add_mfield 0 0 0 0 10
+
+Group selection follows atom identities even when atoms move. Region selection
+instead follows current fractional coordinates along cell vectors a, b, and c,
+using the same coordinate conversion and membership test as heat_lan.
+Bounds must be finite, in [0, 1], and satisfy min < max. Intervals are half-open:
+[min, max). The region scales and tilts with the current cell, including under
+NPT. Coordinates are evaluated after normal periodic boundary wrapping; there
+is no extra wrapping in nonperiodic directions. No group labels are needed for
+region mode.
+
+For a 10 T field along z in the first half of the cell along a::
+
+  add_mfield region 0 0.5 0 1 0 1 0 0 10
+
+Membership is recomputed on every force evaluation: atoms entering the region
+receive the field, and atoms leaving it no longer do. The sharp boundary gives
+a discontinuous Zeeman energy. No mechanical boundary impulse is implemented;
+energy conservation is not guaranteed for atoms crossing this boundary. This
+is a local magnetic drive, not a smooth conservative gradient-field model.
+
+A Spin potential and spin:R:3 are required for a nonzero field. The field acts
+on all selected magnetic moments, independently of the model spin_dof_type
+response mask. It does not impose an integration constraint.
+
+Each command replaces the previous field and group/region selection. The setting
+persists across run commands. To turn it off::
+
+  add_mfield 0 0 0 0 0
+
+Spin is in muB, magnetic force is in eV/muB, and energy is in eV. With
+c_B = 5.7883818060e-5, the added magnetic force is c_B times the field vector;
+the added per-atom Zeeman energy is minus spin dotted with this force.
+Both contributions are recomputed at initial, SIB midpoint, GLSD predictor,
+and endpoint states using the actual stored moments, including the
+unnormalized SIB chord midpoint, and do not accumulate across evaluations.
+Output energies and magnetic forces include these contributions.
+No additional g factor is applied. There is no direct mechanical force or
+virial contribution. Time-dependent fields and linear gradients are not supported.
+
+The field is magnetic induction B, not H in A/m. The integrators consume
+magnetic force in eV/muB, so tesla-valued inputs must undergo this conversion
+before integration; the conversion does not change their gamma convention.
+For example, a 2 muB moment aligned with 1 T has Zeeman energy approximately
+-1.1576e-4 eV and magnetic force +5.7884e-5 eV/muB along the field.