[3/6] Add TSPIN and Spin lattice integration support

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,32 @@
+# Spin regression checks
+
+Build on the target GPU with its CUDA architecture, for example:
+
+```sh
+make -C src -j8 CUDA_ARCH=-arch=sm_70 gpumd nep
+```
+
+Use Python with NumPy and pytest. Set RESULTS to a new directory outside the
+checkout. Multi-GPU checks require four allocated GPUs visible through
+CUDA_VISIBLE_DEVICES. GPUMD_COMMAND selects the runtime binary; otherwise tests
+use src/gpumd. Keep the frozen model/oracle values independent of the candidate.
+
+```sh
+python -m pytest --noconftest -q tests_pytest/test_spin_model_format.py
+python tests_pytest/validate_main_nep_spin3_training.py
+python tests_pytest/validate_spin_training_lifecycle.py
+python tests_pytest/validate_spin_training_multigpu.py --output "$RESULTS/training"
+python tests_pytest/validate_nep_spin3_runtime.py
+python tests_pytest/validate_nep_spin3_runtime.py --block-format
+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 -m pytest --noconftest -q tests_pytest/test_spin_tspin.py tests_pytest/test_spin_lattice_restart.py
+```
+
+Training checks cover five representative layouts, actual energy errors,
+magnetic loss components, model IO and restart. Flexible ZBL is rejected;
+universal ZBL remains supported. The CPU model codec is checked with strict
+C++14 compilation. Build the focused probes with `-std=c++14` as well.
+Initialization checks share `restart_parity.cu`; torque and response edge
+cases share `spin_mforce_split.cu` and `response_fitness.cu`.
+
+Use `-m "not slow"` for routine pytest and `-m slow` for extended statistics.

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,58 @@
+.. _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 every force evaluation using the actual
+stored moments 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.

doc/gpumd/input_parameters/dump_restart.rst

@@ -28,3 +28,19 @@ Caveats
 -------
 This keyword is not propagating.
 That means, its effect will not be passed from one run to the next.
+
+Spin lattice checkpoints
+------------------------
+
+Spin NVT/NPT runs save the atomic thermostat/barostat state and original ramp
+phase as paired ``nhc_state``/``nhc_time`` or ``mttk_state``/``mttk_time`` fields.
+Copy ``restart.xyz`` to ``model.xyz`` and retain the potential, atom order,
+timestep, temperature/pressure endpoints and ensemble parameters. Do not
+regenerate velocities. The original ramp continues and holds its final target
+beyond its saved length. Mismatched configurations are rejected; selecting a
+different lattice driver discards pending state. ``lattice off`` has no lattice
+thermostat checkpoint. Ordinary NHC/MTTK ensembles retain upstream behavior.
+
+This preserves the lattice bath, not every extended variable or random stream
+of the spin integrator. TSPIN's spin thermostat state is not checkpointed.
+Other stateful modifiers and analysis histories are outside this contract.

doc/gpumd/input_parameters/dump_xyz.rst

@@ -29,7 +29,7 @@ If it is omitted, data for the whole system will be output.
 :attr:`single` (the default) writes nine significant digits, which is enough to recover a single-precision value exactly.
 :attr:`double` writes 17 significant digits, which is enough to recover a double-precision value exactly.
 
-* Then one can write the properties to be output, and the allowed properties include: :attr:`mass`, :attr:`velocity`, :attr:`force`, :attr:`potential`, :attr:`virial`, :attr:`charge`, :attr:`bec`, :attr:`group_labels`, and :attr:`unwrapped_position`.
+* Then one can write the properties to be output, and the allowed properties include: :attr:`mass`, :attr:`velocity`, :attr:`force`, :attr:`potential`, :attr:`virial`, :attr:`charge`, :attr:`bec`, :attr:`group_labels`, :attr:`spin`, :attr:`mforce`, :attr:`spin_velocity`, and :attr:`unwrapped_position`.
 
 The :attr:`group_labels` property writes one integer column per grouping method, holding the label of the group each atom belongs to.
 
@@ -63,3 +63,9 @@ Caveats
 * For qNEP models, the charge values dumped out are predicted by the qNEP models.
   For other models, the charge values are those specified in :attr:`model.xyz` via :attr:`charge:R:1`.
 * The Born effective charge (:term:`BEC`) is only meaningful for a qNEP model trained with target :term:`BEC`.
+* ``spin`` and ``mforce`` are three-component vectors available when
+  ``model.xyz`` contains ``spin:R:3``. ``spin_velocity`` additionally requires
+  an initialized ``spin_vel:R:3`` property.
+* The per-atom ``virial`` property contains all nine row-major components
+  ``xx xy xz yx yy yz zx zy zz``. NEP_Spin does not symmetrize this per-atom
+  tensor.

doc/gpumd/input_parameters/ensemble.rst

@@ -10,6 +10,7 @@ There are different categories of methods accessible via this keyword, which are
 
 * :ref:`"standard" ensembles <kw_ensemble_standard>`
 * :ref:`MTTK integrators <kw_ensemble_mttk>`
+* :ref:`TSPIN integrators for Spin potentials <kw_ensemble_tspin>`
 * :ref:`integrators for thermal conductivity simulations <kw_ensemble_heat>`
 * :ref:`integrators for path integral molecular dynamics simulations <kw_ensemble_pimd>`
 * :ref:`MSST integrator for simulating compressive shock wave <kw_ensemble_msst>`
@@ -22,7 +23,7 @@ There are different categories of methods accessible via this keyword, which are
 Units and suggested parameters
 ------------------------------
 
-The units of temperature and pressure for this keyword are K and GPa, respectively. 
+The units of temperature and pressure for this keyword are K and GPa, respectively.
 
 The temperature coupling constant :attr:`<T_coup>` means :math:`\tau_T/\Delta t`, where :math:`\tau_T` is the relaxation time of the thermostat and :math:`\Delta t` is the time step for integration.
 We require :math:`\tau_T/\Delta t \geq 1` and a good choice is :math:`\tau_T/\Delta t \approx 100`.