Add spin-dependent NEP training and spin–lattice dynamics

.gitignore

@@ -32,5 +32,7 @@ __pycache__/
 .hypothesis/
 README.gnep.md
 push-to-both.sh
+/tmp/qscaler-ab/
+/tmp/spin-local-policy/
 
 src/.wpe_obj/

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,155 @@
+# Spin regression checks
+
+Run these commands from the repository root. Use Python with NumPy and pytest,
+a C++ compiler for the format tests, and CUDA with a supported NVIDIA GPU for
+the runtime tests. Load your site's CUDA environment before building; the tests
+inherit the environment and do not configure library paths.
+
+## Model format and runtime
+
+The format tests compile the production codec on the CPU:
+
+```sh
+python -m pytest --noconftest -q tests_pytest/test_spin_model_format.py
+```
+
+Build both executables on the target machine. Set `CUDA_ARCH` for its GPU; this
+example targets V100:
+
+```sh
+make -C src -j8 CUDA_ARCH=-arch=sm_70 gpumd nep
+python tests_pytest/validate_nep_spin3_runtime.py
+python tests_pytest/validate_nep_spin3_runtime.py --block-format
+python -m pytest --noconftest -q tests_pytest/test_spin_fire.py \
+  tests_pytest/test_spin_sib.py tests_pytest/test_spin_tspin.py
+python tests_pytest/validate_spin_tspin_runtime.py
+```
+
+These tests default to `src/gpumd`; `GPUMD_COMMAND` selects another executable.
+`--noconftest` runs the self-contained Spin tests without importing the ordinary
+NEP suite's ASE/Calorine fixtures. Runtime tests require a working GPU executable.
+
+`tests_pytest/fixtures/nep_spin3` contains two legacy models (C=2 and C=3), one
+structure, and three frozen oracle files. The C=2 typewise case reuses the C=2
+uniform parameters and changes only `spin_cutoff 6.0` to `spin_cutoff 4.0 6.0`.
+Each oracle records the SHA256 of its complete model input; the loader verifies
+that hash, including the generated typewise header. The versioned `model nep_spin`
+format and legacy files are checked against the same frozen energy, force,
+magnetic force, and virial values. Lifecycle checks cover prediction, MD,
+restart, and compatibility with existing files and inputs. Unsupported versions
+and conflicting selectors must fail. Training examples use `spin_version 1`.
+Do not regenerate expected values with the candidate implementation when fixing
+a test.
+
+## Integration temperature contract
+
+After building `gpumd`, link this driver probe against its objects, excluding
+`main_gpumd/main.o`. It checks that the initial temperature cache stays fixed,
+while the per-step lattice and Spin targets follow the ramp. It also checks a
+non-temperature-driving ensemble with an unusable base target.
+
+```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"
+```
+
+## Training and multiple GPUs
+
+The following self-contained training checks generate their small datasets in
+temporary directories:
+
+```sh
+python tests_pytest/validate_main_nep_spin3_training.py
+python tests_pytest/validate_spin_training_lifecycle.py
+```
+
+For GPU scaling checks, expose four allocated GPUs through `CUDA_VISIBLE_DEVICES`.
+Choose an output directory outside the checkout:
+
+```sh
+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"
+```
+
+The training comparison uses one, two, and four GPUs, multiple structures per
+batch, and independently recomputed energy errors. The runtime comparison covers
+domain partitioning, boundary conditions, energy, forces, and magnetic forces.
+
+Additional probes must be linked against the current NEP objects, excluding its
+`main` object. On Linux/CUDA, after building `nep`, for example:
+
+```sh
+for probe in spin_mforce_split restart_parity response_fitness; do
+  objects=$(find src/main_nep src/utilities -name '*.cu' ! -path 'src/main_nep/main.cu' \
+    | sed 's/\.cu$/.o/')
+  nvcc -std=c++17 -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
+```
+
+Set `RESULTS` to an existing absolute directory before using these commands.
+The restart check covers both ordinary NEP and Spin, preserves the SNES search
+mean separately from the exported elite, and checks scaler import. The split
+probe checks parallel/perpendicular magnetic losses against an independent
+projection calculation, including zero spins and parser rejection cases.
+
+## GLSD and lattice-driver restart
+
+The standalone GLSD probe checks deterministic convergence and additive-noise
+statistics against independent analytic expectations. It needs no trained model:
+
+```sh
+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"
+GPUMD_COMMAND=src/gpumd python tests_pytest/validate_spin_sib_runtime.py
+```
+
+Full runtime and restart checks require four allocated GPUs, a compatible Fe
+Spin model, and independently preserved reference executables. Use absolute
+paths and fresh output directories outside the checkout:
+
+```sh
+python tests_pytest/validate_spin_glsd_runtime.py \
+  --binary /path/to/current/gpumd --reference /path/to/pre-glsd/gpumd \
+  --mobility-reference /path/to/pre-damping/gpumd \
+  --model /path/to/Fe/nep.txt --output "$RESULTS/glsd"
+python tests_pytest/validate_mttk_restart.py \
+  --binary /path/to/current/gpumd --reference /path/to/pre-isolation/gpumd \
+  --model /path/to/Fe/nep.txt --output "$RESULTS/mttk"
+python tests_pytest/validate_nhc_restart.py \
+  --binary /path/to/current/gpumd --reference /path/to/pre-isolation/gpumd \
+  --model /path/to/Fe/nep.txt --output "$RESULTS/nhc"
+GPUMD_COMMAND=/path/to/current/gpumd python -m pytest --noconftest -q \
+  tests_pytest/test_spin_lattice_restart.py
+```
+
+These reference executables and the Fe model are external prerequisites, not
+bundled fixtures. Record their source revisions and hashes; never substitute
+the candidate for an independent reference. The GLSD mobility reference must
+support the former dimensional input. Restart checks compare continuous and
+resumed trajectories step by step, exercise invalid inputs, and use missing-state
+controls to detect incomplete restoration. The lattice restart reference must
+support the existing endpoint format: its checkpoints are resumed by the
+candidate to check compatibility. Ordinary NHC/MTTK cases must emit no Spin
+lattice metadata and must ignore old bath state when starting an ordinary run.
+See the
+[restart documentation](../doc/gpumd/input_parameters/dump_restart.rst) for the
+supported continuation contracts and the
+[GLSD documentation](../doc/gpumd/input_parameters/ensemble_glsd.rst) for the
+physical equation and timestep limits.
+
+Keep run logs, machine details and numerical reports outside the source tree.
+Passing these regressions establishes the covered software contracts; it does
+not establish magnetic phase behavior or extrapolative model accuracy.

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.