feat(lammps): Cartesian magnetic-force mode for dpa4spin (#6077)

doc/model/dpa4c.md

@@ -398,6 +398,20 @@ already performs, and the magnetic force is reduced back onto owning atoms
 alongside the conservative force, so domain decomposition needs no additional
 exchange.
 
+By default the pair styles store the precession magnetic force expected by the
+LAMMPS SPIN package (the model energy gradient scaled by the moment magnitude
+over `hbar`). A Newtonian Cartesian integrator instead needs the raw energy
+gradient `-dE/dS`. Select it with `magnetic_force cartesian` (requires
+`units metal`); the default is `magnetic_force precession`:
+
+```lammps
+pair_style        dpa4spin/kk compressed_model.pt2 magnetic_force cartesian
+```
+
+The selected convention is exposed as the scalar integer
+`Pair::extract("fm_cartesian_force", dim)` (zero for precession, one for
+Cartesian) so a consuming integrator can reject an incompatible configuration.
+
 `min_style spin` reads the magnetic force from the pair style and relaxes the
 moment directions. Spin dynamics through `fix nve/spin` requires a LAMMPS build
 whose fix recognizes this pair style, because the stock fix accumulates the

examples/spin/dpa4c/lmp/README.md

@@ -67,3 +67,12 @@ they stay fixed through that stage.
 Spin dynamics through `fix nve/spin` requires a LAMMPS build whose fix
 recognizes this pair style, because the stock fix accumulates the magnetic
 force only from pair styles matching its own name pattern.
+
+The default magnetic-force convention is precession (`|m| / hbar` scaling of
+the model energy gradient). A separately maintained Newtonian Cartesian
+integrator can request `-dE/dS` with only the energy-per-model-spin unit
+conversion via `magnetic_force cartesian` on the `pair_style` line
+(`units metal`) and can confirm the choice through
+`Pair::extract("fm_cartesian_force", dim)`. On `dpa4spin/kk`, pair that keyword
+with the compressed compact-canonical archive (not the graph-lower freeze
+output); host `dpa4spin` may use either archive.

examples/spin/dpa4c/lmp/in.lammps

@@ -11,6 +11,15 @@
 # graph and the moment in device memory:
 #
 #   lmp -k on g 1 -sf kk -in in.lammps
+#
+# By default atom->fm holds the precession magnetic force. For a Newtonian
+# Cartesian integrator that needs -dE/dS, append `magnetic_force cartesian`
+# (units metal) and read Pair::extract("fm_cartesian_force", dim). Host
+# `dpa4spin` accepts the graph-lower archive from freeze; `dpa4spin/kk`
+# requires the compact canonical archive from `dp --pt-expt compress`
+# (see README.md):
+#
+#   pair_style        dpa4spin/kk compressed_model.pt2 magnetic_force cartesian
 
 units             metal
 dimension         3

source/install/test_cc_local.sh

@@ -134,6 +134,11 @@ else:
 	# exports NO with-comm artifact (single-rank only by construction).
 	env ${_GEN_ENV} python ${INFER_SCRIPT_PATH}/gen_dpa4_spin_zbl.py &
 	PID14=$!
+	# Minimal native-spin DPA4C for LAMMPS pair_style dpa4spin (no
+	# fparam/aparam/chg_spin). CPU-safe; used by the magnetic-force
+	# convention behavioral tests.
+	env ${_GEN_ENV} python ${INFER_SCRIPT_PATH}/gen_dpa4c_spin_lmp.py &
+	PID14B=$!
 	# Native-spin DPA4C: the archives above all ship a with-comm artifact
 	# and are therefore served by DeepSpinPTExpt, which leaves
 	# NativeSpinPTExpt -- the class the dispatch selects for a native-spin
@@ -149,6 +154,7 @@ else:
 	wait $PID12
 	wait $PID13
 	wait $PID14
+	wait $PID14B
 	if [ "${DP_VARIANT}" = "cuda" ]; then
 		wait $PID15
 	fi

source/lmp/pair_dpa4spin.cpp

@@ -107,9 +107,12 @@ static const char cite_user_deepmd_package[] =
     "}\n\n";
 
 namespace {
-// Reduced Planck constant in eV.ps. The model reports the magnetic force as
-// the energy gradient with respect to the magnetic moment, while LAMMPS stores
-// the precession force; the two differ by the factor hbar / |m|.
+// Reduced Planck constant in eV.ps (metal units). The model reports the
+// magnetic force as the energy gradient with respect to the magnetic moment.
+// Under magnetic_force precession (the default), that gradient is stored in
+// atom->fm after multiplication by |m| / hbar. Under magnetic_force cartesian
+// the gradient is stored with only the energy-per-model-spin unit conversion,
+// and this constant is unused.
 constexpr double kHBar = 6.5821191e-04;
 // Positions of the LAMMPS global virial components (xx, yy, zz, xy, xz, yz)
 // within the nine-component tensor the model reports.
@@ -120,7 +123,11 @@ constexpr int kCentroidVirialMap[9] = {0, 4, 8, 3, 6, 7, 1, 2, 5};
 }  // namespace
 
 PairDPA4Spin::PairDPA4Spin(LAMMPS* lmp)
-    : Pair(lmp), scale(nullptr), cutoff(0.0), commdata_(nullptr) {
+    : Pair(lmp),
+      scale(nullptr),
+      cutoff(0.0),
+      fm_cartesian_force(0),
+      commdata_(nullptr) {
   if (lmp->citeme) {
     lmp->citeme->add(cite_user_deepmd_package);
   }
@@ -237,46 +244,85 @@ void PairDPA4Spin::settings(int narg, char** arg) {
   }
   cutoff = deep_spin.cutoff() * dist_unit_cvt_factor;
 
-  // How many values name a charge state is a property of the artifact, so the
-  // keyword can only be read once the model is loaded.
+  // Optional keywords follow the artifact path. How many values name a charge
+  // state is a property of the artifact, so charge_spin can only be read once
+  // the model is loaded. magnetic_force selects the convention stored in
+  // atom->fm and is independent of whether the model uses native spin.
   const int dim_chg_spin = deep_spin.dim_chg_spin();
   std::vector<double> charge_spin;
+  bool magnetic_force_seen = false;
+  fm_cartesian_force = 0;
   int iarg = 1;
   while (iarg < narg) {
-    if (std::string(arg[iarg]) != "charge_spin") {
+    const std::string key(arg[iarg]);
+    if (key == "charge_spin") {
+      // One charge state holds for the whole run, so a second occurrence names
+      // no state the style could serve. Values already read mark the keyword
+      // as seen: a successful read always contributes at least one value.
+      if (!charge_spin.empty()) {
+        error->all(FLERR,
+                   "Illegal pair_style command: keyword charge_spin names the "
+                   "single charge state of the whole run, so it may appear "
+                   "only once.");
+      }
+      if (dim_chg_spin == 0) {
+        error->all(FLERR,
+                   "Illegal pair_style command: the artifact served by pair "
+                   "style " +
+                       style +
+                       " carries no charge/spin condition, so keyword "
+                       "charge_spin names nothing it can serve.");
+      }
+      if (iarg + dim_chg_spin >= narg) {
+        error->all(FLERR,
+                   "Illegal pair_style command: keyword charge_spin names a "
+                   "charge state with " +
+                       std::to_string(dim_chg_spin) + " value(s).");
+      }
+      for (int ii = 0; ii < dim_chg_spin; ++ii) {
+        charge_spin.push_back(
+            utils::numeric(FLERR, arg[iarg + 1 + ii], false, lmp));
+      }
+      iarg += 1 + dim_chg_spin;
+    } else if (key == "magnetic_force") {
+      if (magnetic_force_seen) {
+        error->all(FLERR,
+                   "Illegal pair_style command: keyword magnetic_force may "
+                   "appear only once.");
+      }
+      if (iarg + 1 >= narg) {
+        error->all(FLERR,
+                   "Illegal pair_style command: keyword magnetic_force takes "
+                   "one of precession or cartesian.");
+      }
+      const std::string value(arg[iarg + 1]);
+      if (value == "precession") {
+        fm_cartesian_force = 0;
+      } else if (value == "cartesian") {
+        fm_cartesian_force = 1;
+      } else {
+        error->all(FLERR,
+                   "Illegal pair_style command: keyword magnetic_force takes "
+                   "one of precession or cartesian, not '" +
+                       value + "'.");
+      }
+      magnetic_force_seen = true;
+      iarg += 2;
+    } else {
       error->all(FLERR, "Illegal pair_style command: pair style " + style +
                             " takes the artifact path followed by the optional "
-                            "keyword charge_spin, not '" +
-                            std::string(arg[iarg]) + "'.");
-    }
-    // One charge state holds for the whole run, so a second occurrence names
-    // no state the style could serve.  Values already read mark the keyword
-    // as seen: a successful read always contributes at least one value.
-    if (!charge_spin.empty()) {
-      error->all(FLERR,
-                 "Illegal pair_style command: keyword charge_spin names the "
-                 "single charge state of the whole run, so it may appear "
-                 "only once.");
-    }
-    if (dim_chg_spin == 0) {
-      error->all(FLERR,
-                 "Illegal pair_style command: the artifact served by pair "
-                 "style " +
-                     style +
-                     " carries no charge/spin condition, so keyword "
-                     "charge_spin names nothing it can serve.");
-    }
-    if (iarg + dim_chg_spin >= narg) {
-      error->all(FLERR,
-                 "Illegal pair_style command: keyword charge_spin names a "
-                 "charge state with " +
-                     std::to_string(dim_chg_spin) + " value(s).");
+                            "keywords charge_spin and magnetic_force, not '" +
+                            key + "'.");
     }
-    for (int ii = 0; ii < dim_chg_spin; ++ii) {
-      charge_spin.push_back(
-          utils::numeric(FLERR, arg[iarg + 1 + ii], false, lmp));
-    }
-    iarg += 1 + dim_chg_spin;
+  }
+
+  // The Cartesian convention stores an energy-per-spin gradient. Its only
+  // validated unit system is metal, where the model eV and the LAMMPS energy
+  // unit coincide and the precession path's metal hbar remains consistent.
+  if (fm_cartesian_force && strcmp(update->unit_style, "metal") != 0) {
+    error->all(FLERR,
+               "pair style " + style +
+                   " keyword magnetic_force cartesian requires units metal.");
   }
 
   // A charge/spin condition named on the pair_style line holds for the whole
@@ -305,6 +351,8 @@ void PairDPA4Spin::settings(int narg, char** arg) {
     }
     utils::logmesg(lmp, "  using charge_spin:  {}\n", values);
   }
+  utils::logmesg(lmp, "  magnetic_force:     {}\n",
+                 fm_cartesian_force ? "cartesian" : "precession");
 }
 
 /* ----------------------------------------------------------------------
@@ -536,16 +584,23 @@ void PairDPA4Spin::compute(int eflag, int vflag) {
   // === Step 3. Accumulate the model outputs ===
   // The model reports a force and a magnetic force for every atom the neighbor
   // list covers, ghosts included; the spin atom style folds the ghost rows
-  // onto their owners in its reverse communication. Dividing the magnetic
-  // force by hbar / |m| turns the energy gradient with respect to the moment
-  // into the precession force LAMMPS stores. Multiplication by |m| preserves
-  // that relation while making a zero moment well-defined.
+  // onto their owners in its reverse communication. In the default precession
+  // mode, dividing the magnetic force by hbar / |m| turns the energy gradient
+  // with respect to the moment into the precession force LAMMPS stores;
+  // multiplication by |m| preserves that relation while making a zero moment
+  // well-defined. In Cartesian mode the energy gradient is stored directly,
+  // with only the energy-per-model-spin unit conversion, including at zero
+  // moment.
   const double force_scale = scale[1][1] * force_unit_cvt_factor;
-  const double magnetic_force_scale = force_scale / kHBar;
+  const double magnetic_force_scale = fm_cartesian_force
+                                          ? scale[1][1] * ener_unit_cvt_factor
+                                          : force_scale / kHBar;
   for (int ii = 0; ii < nall; ++ii) {
+    const double moment_factor = fm_cartesian_force ? 1.0 : sp[ii][3];
     for (int dd = 0; dd < 3; ++dd) {
       f[ii][dd] += dforce[3 * ii + dd] * force_scale;
-      fm[ii][dd] += dforce_mag[3 * ii + dd] * sp[ii][3] * magnetic_force_scale;
+      fm[ii][dd] +=
+          dforce_mag[3 * ii + dd] * moment_factor * magnetic_force_scale;
     }
   }
 
@@ -579,5 +634,9 @@ void* PairDPA4Spin::extract(const char* str, int& dim) {
     dim = 2;
     return (void*)scale;
   }
+  if (strcmp(str, "fm_cartesian_force") == 0) {
+    dim = 0;
+    return (void*)&fm_cartesian_force;
+  }
   return nullptr;
 }

source/lmp/pair_dpa4spin.h

@@ -66,8 +66,10 @@ class PairDPA4Spin : public Pair {
   PairDPA4Spin(class LAMMPS*);
   ~PairDPA4Spin() override;
 
-  // Load the artifact named by ``pair_style dpa4spin <model> [charge_spin
-  // <v1> ... <vN>]`` and fix the charge state the keyword names.
+  // Load the artifact named by ``pair_style dpa4spin <model>
+  // [charge_spin <v1> ... <vN>] [magnetic_force precession|cartesian]`` and
+  // fix the optional charge state and magnetic-force convention the keywords
+  // name. ``magnetic_force`` defaults to precession.
   void settings(int, char**) override;
   // Resolve the LAMMPS atom types onto the element list of the model.
   void coeff(int, char**) override;
@@ -98,6 +100,14 @@ class PairDPA4Spin : public Pair {
   std::vector<int> type_idx_map;
   // LAMMPS unit system relative to the eV / Angstrom system of the model.
   double ener_unit_cvt_factor, dist_unit_cvt_factor, force_unit_cvt_factor;
+  // Magnetic-force convention stored in ``atom->fm``. Zero selects the
+  // precession force (energy gradient scaled by the moment magnitude over
+  // hbar), which is the default and what the stock LAMMPS SPIN integrator
+  // expects. One selects the Cartesian energy gradient ``-dE/dS`` with only
+  // the energy-per-model-spin unit conversion, which a Newtonian spin
+  // integrator consumes. Exposed to callers through
+  // ``extract("fm_cartesian_force", dim)``.
+  int fm_cartesian_force;
 
  private:
   CommBrickDPA4Spin* commdata_;

source/lmp/pair_dpa4spin_kokkos.cpp

@@ -18,9 +18,12 @@
 using namespace LAMMPS_NS;
 
 namespace {
-// Reduced Planck constant in eV.ps. The model reports the magnetic force as
-// the energy gradient with respect to the magnetic moment, while LAMMPS stores
-// the precession force; the two differ by the moment magnitude over hbar.
+// Reduced Planck constant in eV.ps (metal units). The model reports the
+// magnetic force as the energy gradient with respect to the magnetic moment.
+// Under magnetic_force precession (the default), that gradient is stored in
+// atom->fm after multiplication by |m| / hbar. Under magnetic_force cartesian
+// the gradient is stored with only the energy-per-model-spin unit conversion,
+// and this constant is unused.
 constexpr double kHBar = 6.5821191e-04;
 }  // namespace
 
@@ -333,7 +336,11 @@ void PairDPA4SpinKokkos<DeviceType>::compute(int eflag, int vflag) {
   // communication stages.
   auto model2loc = compact_graph.d_model2loc;
   const double fscale = scale[1][1] * force_unit_cvt_factor;
-  const double fmscale = scale[1][1] * force_unit_cvt_factor / kHBar;
+  // Precession: energy gradient times |m| / hbar. Cartesian: energy gradient
+  // with only the energy-per-model-spin unit conversion (no |m|, no hbar).
+  const int cartesian = fm_cartesian_force;
+  const double fmscale =
+      cartesian ? scale[1][1] * ener_unit_cvt_factor : fscale / kHBar;
   reverse_virial = false;
   reverse_used_host = false;
   // The KOKKOS package runs 'newton off', so the integrator's force_clear only
@@ -364,15 +371,14 @@ void PairDPA4SpinKokkos<DeviceType>::compute(int eflag, int vflag) {
       "dpa4spin/kk:scatter_f", Kokkos::RangePolicy<DeviceType>(0, nnode_m),
       KOKKOS_LAMBDA(const int m) {
         const int i = model2loc(m);
-        // The precession force carried by ``fm`` is the energy gradient with
-        // respect to the moment, rescaled by the moment magnitude over hbar.
-        const double moment = fmscale * sp(i, 3);
+        // Precession multiplies by the moment magnitude; Cartesian does not.
+        const double moment_factor = cartesian ? 1.0 : sp(i, 3);
         f(i, 0) += fscale * out_force(3 * m + 0);
         f(i, 1) += fscale * out_force(3 * m + 1);
         f(i, 2) += fscale * out_force(3 * m + 2);
-        fm(i, 0) += moment * out_force_mag(3 * m + 0);
-        fm(i, 1) += moment * out_force_mag(3 * m + 1);
-        fm(i, 2) += moment * out_force_mag(3 * m + 2);
+        fm(i, 0) += fmscale * moment_factor * out_force_mag(3 * m + 0);
+        fm(i, 1) += fmscale * moment_factor * out_force_mag(3 * m + 1);
+        fm(i, 2) += fmscale * moment_factor * out_force_mag(3 * m + 2);
       });
   atomKK->modified(execution_space, F_MASK | FM_MASK);
   if (multi_rank) {

source/lmp/plugin/deepmdplugin.cpp

@@ -13,6 +13,7 @@
 #endif
 #include "pair_deepmd.h"
 #include "pair_deepspin.h"
+#include "pair_dpa4spin.h"
 #include "version.h"
 #if LAMMPS_VERSION_NUMBER >= 20220328
 #include "pppm_dplr.h"
@@ -22,6 +23,7 @@ using namespace LAMMPS_NS;
 
 static Pair* pairdeepmd(LAMMPS* lmp) { return new PairDeepMD(lmp); }
 static Pair* pairdeepspin(LAMMPS* lmp) { return new PairDeepSpin(lmp); }
+static Pair* pairdpa4spin(LAMMPS* lmp) { return new PairDPA4Spin(lmp); }
 
 #ifdef LMP_KOKKOS
 // Runtime plugins do not consume the PairStyle declarations used by LAMMPS'
@@ -109,6 +111,15 @@ extern "C" void lammpsplugin_init(void* lmp, void* handle, void* regfunc) {
   plugin.handle = handle;
   (*register_plugin)(&plugin, lmp);
 
+  plugin.version = LAMMPS_VERSION;
+  plugin.style = "pair";
+  plugin.name = "dpa4spin";
+  plugin.info = "dpa4spin pair style " STR_GIT_SUMM;
+  plugin.author = "Tiancheng Li";
+  plugin.creator.v1 = (lammpsplugin_factory1*)&pairdpa4spin;
+  plugin.handle = handle;
+  (*register_plugin)(&plugin, lmp);
+
   plugin.style = "compute";
   plugin.name = "deeptensor/atom";
   plugin.info = "compute deeptensor/atom " STR_GIT_SUMM;