[2/6] Add the NEP-spin potential to GPUMD

CMakeLists.txt

@@ -30,7 +30,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

@@ -55,8 +55,8 @@ The main source directories are:
 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`.
 
 Training in `src/main_nep/nep_spin/` uses the shared magnetic descriptor and
-force headers in `src/force/nep_spin/`. GPUMD Spin force dispatch is added
-separately; the headers do not enable a runtime potential by themselves.
+force headers in `src/force/nep_spin/`. The GPUMD potential, fixed-spin runtime support and spatial decomposition
+also live in `src/force/nep_spin/`. See [runtime validation](spin_runtime_validation.md).
 
 ## 3. Main `gpumd` architecture
 

developer/spin_runtime_validation.md

@@ -0,0 +1,53 @@
+# Fixed-spin runtime validation
+
+Build `gpumd` and `nep` on a CUDA node, choosing the device architecture:
+
+```sh
+make -C src -j8 CUDA_ARCH=-arch=sm_70 gpumd nep
+export GPUMD_COMMAND="$PWD/src/gpumd"
+export NEP_COMMAND="$PWD/src/nep"
+```
+
+Use Python with NumPy and pytest. Set `RESULTS` to a new output directory outside
+this checkout. The first two GPU scripts need one GPU; run them with one
+allocated device visible. The third requires four allocated GPUs in one process.
+
+```sh
+python -m pytest --noconftest -q tests_pytest/test_spin_model_format.py tests_pytest/test_spin_validation.py
+python tests_pytest/validate_nep_spin3_runtime.py
+python tests_pytest/validate_mfield_runtime.py
+python tests_pytest/validate_spin_multigpu.py --binary "$GPUMD_COMMAND" --model tests_pytest/fixtures/nep_spin3/nep_spin_o3c2_uniform.nep --output "$RESULTS/runtime"
+```
+
+The runtime oracle checks cover energy, atomic force, magnetic force and all
+nine virial components, plus energy-reference invariance, excluded-spin
+masking, replication and fixed-spin restart. Model-selection checks include
+uppercase fields, attached comments and a model record after numeric blocks.
+Small-cell checks cover capacity growth and a 10,000-atom thin cell compared
+with an equivalent expanded cell. Split XYZ output is checked after restart. The frozen reference provenance
+is documented in [fixtures](../tests_pytest/fixtures/nep_spin3/README.md).
+The external-field script checks analytic Zeeman energy and magnetic force,
+group and region selection, boundary crossing, excluded spins and invalid inputs.
+Multi-GPU comparisons use componentwise tolerances and reject nonfinite values.
+They cover small cells, nonuniform and triclinic systems, cutoff boundaries,
+partition directions, motion across partition boundaries and energy references.
+An additional supported `n_max 4 8`, `l_max 8 2 1` model needs 51,840 bytes
+in the direct angular-force path and therefore exercises its 48 KiB fallback
+on one, two and four GPUs with nonzero forces.
+
+Run the [training validation](spin_validation.md) against the same build to
+check the shared descriptor implementation and native/JIT agreement. Keep its
+model, seed, budget and GPU settings fixed when comparing revisions.
+
+For ordinary GPUMD and NEP behavior, build an unmodified upstream checkout and
+run the full regression suite. `BASELINE_GPUMD` and `BASELINE_NEP` must identify
+those independently built upstream executables:
+
+```sh
+python tests_regression/run_regression.py --repo-root "$PWD" \
+  --baseline "$BASELINE_GPUMD" --candidate "$GPUMD_COMMAND" \
+  --baseline-nep "$BASELINE_NEP" --candidate-nep "$NEP_COMMAND" --suite full
+```
+
+These checks validate implementation contracts, not fit quality or performance.
+Spin integrators and their checkpoint state are outside this suite.

developer/spin_validation.md

@@ -90,3 +90,5 @@ implementation to make a failing check pass. Models use `model nep_spin`,
 
 Keep build logs, source/binary hashes, timing and numerical reports outside the
 source tree.
+
+For GPUMD fixed-spin inference and restart checks, see [Spin runtime validation](spin_runtime_validation.md).

doc/gpumd/input_files/model_xyz.rst

@@ -51,6 +51,7 @@ 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*)
   * ``group:I:number_of_grouping_methods`` grouping methods (*Optional*)
 
 
@@ -67,6 +68,16 @@ 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).
+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.
+Active spins are not normalized; training and runtime data must use the same
+spin convention. Elements excluded by the model's ``spin_type`` have zero
+effective spin and zero magnetic force, while their geometry is retained.
+The ``spin`` property 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,54 @@
+.. _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). Coordinates are evaluated after periodic boundary wrapping.
+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
+only on selected atoms whose elements belong to the model's ``spin_type``.
+Other elements retain zero effective spin and zero magnetic force.
+
+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.
+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,5 @@ Caveats
 -------
 This keyword is not propagating.
 That means, its effect will not be passed from one run to the next.
+
+NEP-spin models also save the per-atom Cartesian ``spin:R:3`` values.

doc/gpumd/input_parameters/dump_xyz.rst

@@ -34,7 +34,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`, 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.
 
@@ -72,3 +72,8 @@ 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``.
+* 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/index.rst

@@ -38,6 +38,7 @@ Below you can find a listing of keywords for the ``run.in`` input file.
    ensemble_nphug
    add_force
    add_efield
+   add_mfield
    add_spring
    deposit
    electron_stop

doc/gpumd/input_parameters/potential.rst

@@ -30,7 +30,7 @@ relative or absolute path) of the potential file to be used. Some potential
 types require additional parameters, as described in their corresponding
 sections.
 
-For standard NEP and temperature-dependent NEP models, a multi-GPU run can
+For standard NEP, temperature-dependent NEP, and ``nep_spin`` models, a multi-GPU run can
 optionally specify the spatial partition direction::
 
   potential <potential_filename> [partition_direction]
@@ -41,7 +41,7 @@ triclinic box. If the direction is omitted, GPUMD selects it automatically,
 normally along the thickest direction of the box.
 
 The optional partition direction applies only to the multi-GPU implementation
-of standard NEP and temperature-dependent NEP models. Charge NEP models
+of standard NEP, temperature-dependent NEP, and ``nep_spin`` models. Charge NEP models
 currently use the single-GPU implementation in ``gpumd``.
 
 Example
@@ -55,3 +55,17 @@ For a multi-GPU NEP run, one can force the decomposition along a specific
 direction, for example::
 
    potential Si_NEP.txt x
+
+NEP-spin
+--------
+
+A :ref:`NEP-spin model <nep_spin>` uses all GPUs visible to the process to
+compute the same atomic system. The optional partition direction also applies
+to NEP-spin. Small systems may not benefit from multiple GPUs; global atom
+arrays remain on GPU 0.
+
+Provide ``spin:R:3`` in :ref:`model.xyz <model_xyz>`. Fully periodic orthogonal
+and triclinic cells are supported. Spins remain fixed under ``nve`` and
+``nvt_nhc``; other ensembles and multiple-potential runs are not supported.
+See :ref:`the Spin runtime instructions <nep_spin_runtime>` for output,
+restart and external-field usage.

doc/installation.rst

@@ -97,7 +97,7 @@ use ``CUDA_VISIBLE_DEVICES`` to select a subset, for example::
   CUDA_VISIBLE_DEVICES=0,1 ./nep
 
 The ``gpumd`` executable also uses all visible GPUs when running standard NEP
-or temperature-dependent NEP models that support the multi-GPU implementation.
+or temperature-dependent NEP models, and ``nep_spin`` spin models.
 For example::
 
   CUDA_VISIBLE_DEVICES=0,1 ./gpumd

doc/potentials/nep_spin.rst

@@ -609,3 +609,45 @@ selected Spin layout. The magnetic kernels share their formulas with the
 native path. If runtime compilation is unavailable, native kernels remain
 available. Compilation adds startup time, so very short runs may be slower;
 ``nep_compile 0`` disables specialization. Prediction uses the native path.
+
+.. _nep_spin_runtime:
+
+Fixed-spin molecular dynamics
+-----------------------------
+
+Load the exported model with ``potential nep.txt`` and provide Cartesian
+magnetic moments through ``spin:R:3`` in :ref:`model.xyz <model_xyz>`.
+The ``nve`` and ``nvt_nhc`` ensembles evolve positions and velocities while
+keeping the magnetic moments fixed. Only one potential is allowed, and all
+three cell directions must be periodic.
+
+Atoms whose elements are absent from ``spin_type`` have zero effective spin.
+Their input spin values are ignored; they still contribute geometric
+environments. Their magnetic forces remain zero, including in an applied
+external field. The reader does not normalize the active spins.
+
+For example::
+
+   potential nep.txt
+   velocity 300
+   ensemble nve
+   time_step 1
+   dump_xyz 100 trajectory.xyz spin mforce
+   dump_restart 1000
+   run 1000
+
+Use ``dump_xyz`` with ``spin`` and ``mforce`` to write the moments and magnetic
+forces. ``dump_restart`` preserves the fixed spins in ``restart.xyz``;
+copy this file to ``model.xyz`` and use the same potential to resume.
+``replicate`` copies the spin vectors with their atoms.
+
+The :ref:`add_mfield keyword <kw_add_mfield>` applies a magnetic induction in
+tesla to a group or fractional-coordinate region. Magnetic force is
+:math:`-\partial E/\partial \mathbf{s}` in eV/Bohr magneton; it is evaluated
+although spins are fixed. Spin dynamics and spin velocities are not part of
+this runtime interface.
+
+Multi-GPU execution uses every visible GPU in a single process. Use
+``potential nep.txt x`` (or ``y``/``z``) to select the spatial partition
+direction. Small cells may require replicated geometry and offer little
+speedup. GPU 0 retains the global arrays.

src/force/force.cu

@@ -46,6 +46,7 @@ The driver class calculating force and related quantities.
 #include "utilities/gpu_macro.cuh"
 #include "utilities/nep_model_file.cuh"
 #include "utilities/nep_model_input.cuh"
+#include "utilities/nep_spin_model_file.cuh"
 #include "utilities/read_file.cuh"
 #include "utilities/run_input.cuh"
 #include <cstring>
@@ -114,6 +115,13 @@ std::unique_ptr<Potential> Force::create_potential(
              strcmp(potential_name, "qnep2_vdw") == 0) {
     potential.reset(new NEP_Charge_VDW(tokens[1].c_str(), number_of_atoms, run_input));
     is_nep = true;
+  } else if (strcmp(potential_name, "nep_spin") == 0) {
+    if (!potentials.empty()) {
+      PRINT_INPUT_ERROR("NEP_Spin must be the only potential in a run.");
+    }
+    potential = spin_runtime_.create_potential(tokens, number_of_atoms);
+    is_nep = true;
+    has_spin_potential_ = true;
   } else if (
     strcmp(potential_name, "nep4_charge1") == 0 ||
     strcmp(potential_name, "nep4_charge2") == 0 ||
@@ -234,7 +242,9 @@ void Force::parse_potential(
     PRINT_INPUT_ERROR("reading error for potential file.");
   }
   int num_types = get_number_of_types(fid_potential);
-  if (NEP_Model_Input::is_new_format(potential_name)) {
+  if (NEP_Model_Input::is_new_format(potential_name) && NEP_Spin_Model_File::is_spin_model(tokens[1])) {
+    std::strcpy(potential_name, "nep_spin");
+  } else if (NEP_Model_Input::is_new_format(potential_name)) {
     NEP_Model_File model;
     std::ifstream input(tokens[1]);
     std::string error;
@@ -276,6 +286,9 @@ void Force::parse_potential(
   // Move the pointer into the list of potentials
   potentials.push_back(std::move(potential));
   // Check if a non-NEP potential has previously been defined
+  if (has_spin_potential_ && potentials.size() != 1) {
+    PRINT_INPUT_ERROR("NEP_Spin must be the only potential in a run.");
+  }
   has_non_nep = has_non_nep || !is_nep;
   if (potentials.size() > 1 && has_non_nep) {
     PRINT_INPUT_ERROR("Multiple potentials may only be used with NEP potentials.\n");
@@ -613,3 +626,44 @@ void Force::compute(
     virial_per_atom);
 
 }
+
+void Force::compute_spin(
+  Box& box,
+  GPU_Vector<double>& position_per_atom,
+  GPU_Vector<int>& type,
+  std::vector<Group>& group,
+  GPU_Vector<double>& potential_per_atom,
+  GPU_Vector<double>& force_per_atom,
+  GPU_Vector<double>& virial_per_atom,
+  GPU_Vector<double>& spin_per_atom,
+  GPU_Vector<double>& mforce_per_atom,
+  int* position_image)
+{
+  if (!has_spin_potential_ || potentials.size() != 1) {
+    PRINT_INPUT_ERROR("Spin-aware force dispatch requires one NEP_Spin potential.");
+  }
+
+  const int number_of_atoms = type.size();
+  if (
+    spin_per_atom.size() != static_cast<std::size_t>(number_of_atoms) * 3 ||
+    mforce_per_atom.size() != static_cast<std::size_t>(number_of_atoms) * 3) {
+    PRINT_INPUT_ERROR("NEP_Spin requires three spin and mforce components per atom.");
+  }
+  prepare_compute(
+    number_of_atoms, box, position_per_atom, potential_per_atom,
+    force_per_atom, virial_per_atom, position_image);
+  mforce_per_atom.fill(0.0);
+  GPU_CHECK_KERNEL
+
+  potentials[0]->compute(
+    box,
+    type,
+    position_per_atom,
+    spin_per_atom,
+    potential_per_atom,
+    force_per_atom,
+    virial_per_atom,
+    mforce_per_atom);
+  spin_runtime_.apply_field(
+    box, position_per_atom, group, type, spin_per_atom, potential_per_atom, mforce_per_atom);
+}