diff --git a/doc/model/dpa4c.md b/doc/model/dpa4c.md index e75f286b54..df66b340d0 100644 --- a/doc/model/dpa4c.md +++ b/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 diff --git a/examples/spin/dpa4c/lmp/README.md b/examples/spin/dpa4c/lmp/README.md index aa4153f4f6..d13be52950 100644 --- a/examples/spin/dpa4c/lmp/README.md +++ b/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. diff --git a/examples/spin/dpa4c/lmp/in.lammps b/examples/spin/dpa4c/lmp/in.lammps index 050ce53903..d55fde1cdf 100644 --- a/examples/spin/dpa4c/lmp/in.lammps +++ b/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 diff --git a/source/install/test_cc_local.sh b/source/install/test_cc_local.sh index 17b4f2f593..8fc4fcc4d9 100755 --- a/source/install/test_cc_local.sh +++ b/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 diff --git a/source/lmp/pair_dpa4spin.cpp b/source/lmp/pair_dpa4spin.cpp index 0a8d6e1ba8..442a5f620b 100644 --- a/source/lmp/pair_dpa4spin.cpp +++ b/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 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; } diff --git a/source/lmp/pair_dpa4spin.h b/source/lmp/pair_dpa4spin.h index 7d0aba8977..63b3fb6e64 100644 --- a/source/lmp/pair_dpa4spin.h +++ b/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 [charge_spin - // ... ]`` and fix the charge state the keyword names. + // Load the artifact named by ``pair_style dpa4spin + // [charge_spin ... ] [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 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_; diff --git a/source/lmp/pair_dpa4spin_kokkos.cpp b/source/lmp/pair_dpa4spin_kokkos.cpp index f479b27559..50f1360eec 100644 --- a/source/lmp/pair_dpa4spin_kokkos.cpp +++ b/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::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::compute(int eflag, int vflag) { "dpa4spin/kk:scatter_f", Kokkos::RangePolicy(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) { diff --git a/source/lmp/plugin/deepmdplugin.cpp b/source/lmp/plugin/deepmdplugin.cpp index 9bc0b685d2..46922f4ff1 100644 --- a/source/lmp/plugin/deepmdplugin.cpp +++ b/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; diff --git a/source/lmp/tests/test_dpa4spin_magnetic_force_convention.py b/source/lmp/tests/test_dpa4spin_magnetic_force_convention.py new file mode 100644 index 0000000000..af1bbe38b1 --- /dev/null +++ b/source/lmp/tests/test_dpa4spin_magnetic_force_convention.py @@ -0,0 +1,506 @@ +# SPDX-License-Identifier: LGPL-3.0-or-later +"""Checks for the ``dpa4spin`` / ``dpa4spin/kk`` magnetic-force conventions. + +The LAMMPS adapter stores the model energy gradient with respect to the +magnetic moment in ``atom->fm`` under one of two explicit conventions, +selected by ``pair_style ... magnetic_force precession|cartesian`` and +reported by ``Pair::extract("fm_cartesian_force", dim)``: + +* precession (default, extract value 0): ``fm = scale * force_cvt * |m| / hbar + * force_mag``, matching the stock LAMMPS SPIN integrator; +* cartesian (extract value 1): ``fm = scale * ener_cvt * force_mag``, the raw + ``-dE/dS`` gradient with only the energy-per-model-spin unit conversion, + including at zero moment. + +NumPy helpers below document the conversion algebra. Behavioral tests drive +the registered ``dpa4spin`` pair style against +``deeppot_dpa4c_spin_lmp.pt2`` (no fparam/aparam/chg_spin) when PyLammps and +the style are available; they skip otherwise so CI without a dpa4spin build +stays green. ``dpa4spin/kk`` requires a GPU Kokkos backend and a compressed +canonical archive, so the Kokkos case is an explicit skip when unavailable. +""" + +from __future__ import ( + annotations, +) + +import ctypes +import os +import tempfile +from pathlib import ( + Path, +) + +import numpy as np +import pytest + +# Metal-units reduced Planck constant used by pair_dpa4spin{,_kokkos}.cpp. +HBAR_METAL = 6.5821191e-04 + +_MODEL_ENV = "DPA4SPIN_LMP_MODEL" +_DEFAULT_MODEL = ( + Path(__file__).resolve().parents[2] + / "tests" + / "infer" + / "deeppot_dpa4c_spin_lmp.pt2" +) + +# 6-atom NiO geometry shared with source/tests/infer/gen_dpa4c_spin.py +# (box edge 6 A = DPA4C ghost cutoff). Ni is spin-active; O is not. +_BOX = np.array([0.0, 6.0, 0.0, 6.0, 0.0, 6.0, 0.0, 0.0, 0.0]) +_COORD = np.array( + [ + [1.0, 1.0, 1.0], + [3.2, 1.4, 1.1], + [1.3, 1.8, 1.0], + [0.4, 1.2, 1.6], + [3.6, 2.0, 1.3], + [3.4, 0.7, 1.7], + ], + dtype=np.float64, +) +# Atom 0: spin-active Ni at zero moment (cartesian fm must stay nonzero). +# Atoms 1-2: nonzero Ni moments. Atoms 3-5: non-magnetic O. +_SPIN = np.array( + [ + [0.0, 0.0, 0.0], + [-0.07, 0.09, 0.03], + [0.02, -0.06, 0.08], + [0.0, 0.0, 0.0], + [0.0, 0.0, 0.0], + [0.0, 0.0, 0.0], + ], + dtype=np.float64, +) +_TYPE = np.array([1, 1, 1, 2, 2, 2], dtype=np.int32) + + +def store_magnetic_force( + force_mag: np.ndarray, + spin_norm: np.ndarray, + *, + cartesian: bool, + scale: float = 1.0, + ener_unit_cvt_factor: float = 1.0, + force_unit_cvt_factor: float = 1.0, +) -> np.ndarray: + """Apply the adapter's magnetic-force conversion (metal path by default). + + Parameters + ---------- + force_mag : np.ndarray + Model energy gradient with respect to the moment, shape (natoms, 3). + spin_norm : np.ndarray + Per-atom moment magnitude ``sp[:, 3]``, shape (natoms,). + cartesian : bool + True for ``magnetic_force cartesian``, False for precession. + scale : float + Pair ``scale[1][1]`` prefactor. + ener_unit_cvt_factor, force_unit_cvt_factor : float + LAMMPS unit conversions relative to the model eV / Angstrom system. + + Returns + ------- + np.ndarray + Values that the adapter accumulates into ``atom->fm``, shape (natoms, 3). + """ + force_mag = np.asarray(force_mag, dtype=np.float64) + spin_norm = np.asarray(spin_norm, dtype=np.float64) + if cartesian: + return scale * ener_unit_cvt_factor * force_mag + return (scale * force_unit_cvt_factor / HBAR_METAL) * spin_norm[:, None] * force_mag + + +@pytest.mark.parametrize( + "spin_norm", + [ + np.array([1.2737, 0.05, 0.0, 0.0]), + np.array([0.0, 0.0]), + np.array([2.5]), + ], +) +def test_cartesian_preserves_model_gradient_including_zero_moment( + spin_norm: np.ndarray, +) -> None: + """Cartesian mode must not multiply by |m| or divide by hbar.""" + rng = np.random.default_rng(0) + force_mag = rng.normal(size=(spin_norm.shape[0], 3)) + stored = store_magnetic_force(force_mag, spin_norm, cartesian=True) + np.testing.assert_allclose(stored, force_mag) + # Zero-moment atoms keep a non-zero Cartesian force when the model does. + zero = spin_norm == 0.0 + if np.any(zero): + assert np.any(np.abs(stored[zero]) > 0.0) + + +@pytest.mark.parametrize( + "spin_norm", + [ + np.array([1.2737, 0.05, 0.0, 0.0]), + np.array([0.0, 0.0]), + ], +) +def test_precession_scales_by_moment_over_hbar(spin_norm: np.ndarray) -> None: + """Precession mode multiplies by |m| / hbar; zero moment stays zero.""" + rng = np.random.default_rng(1) + force_mag = rng.normal(size=(spin_norm.shape[0], 3)) + stored = store_magnetic_force(force_mag, spin_norm, cartesian=False) + expected = force_mag * (spin_norm / HBAR_METAL)[:, None] + np.testing.assert_allclose(stored, expected) + zero = spin_norm == 0.0 + if np.any(zero): + np.testing.assert_array_equal(stored[zero], 0.0) + + +def test_precession_matches_cartesian_times_moment_over_hbar() -> None: + """On metal units the two conventions differ only by |m| / hbar.""" + spin_norm = np.array([1.1, 0.02, 0.0]) + force_mag = np.array( + [ + [0.3, -0.1, 0.2], + [1.0, 0.0, -0.5], + [4.0, 5.0, 6.0], + ], + dtype=np.float64, + ) + cartesian = store_magnetic_force(force_mag, spin_norm, cartesian=True) + precession = store_magnetic_force(force_mag, spin_norm, cartesian=False) + np.testing.assert_allclose( + precession, cartesian * (spin_norm / HBAR_METAL)[:, None] + ) + + +def test_adapter_declares_extract_key_and_keyword() -> None: + """The C++ adapter must keep the discoverable selection and extract API.""" + root = Path(__file__).resolve().parents[1] + pair_cpp = (root / "pair_dpa4spin.cpp").read_text(encoding="utf-8") + pair_h = (root / "pair_dpa4spin.h").read_text(encoding="utf-8") + pair_kk = (root / "pair_dpa4spin_kokkos.cpp").read_text(encoding="utf-8") + assert 'strcmp(str, "fm_cartesian_force")' in pair_cpp + assert "fm_cartesian_force" in pair_h + assert 'key == "magnetic_force"' in pair_cpp + assert 'value == "cartesian"' in pair_cpp + assert 'value == "precession"' in pair_cpp + assert "dim = 0" in pair_cpp + assert 'strcmp(update->unit_style, "metal")' in pair_cpp + # Primary issue target is dpa4spin/kk: scatter must keep both conventions. + assert "fm_cartesian_force" in pair_kk + assert "ener_unit_cvt_factor" in pair_kk + assert "cartesian ? 1.0 : sp(i, 3)" in pair_kk + + +def test_non_unit_energy_conversion_applies_only_to_cartesian() -> None: + """Cartesian uses ener_unit_cvt_factor; precession uses force_unit_cvt_factor.""" + spin_norm = np.array([1.0]) + force_mag = np.array([[2.0, 0.0, 0.0]]) + cartesian = store_magnetic_force( + force_mag, + spin_norm, + cartesian=True, + ener_unit_cvt_factor=23.060549, + force_unit_cvt_factor=99.0, + ) + precession = store_magnetic_force( + force_mag, + spin_norm, + cartesian=False, + ener_unit_cvt_factor=23.060549, + force_unit_cvt_factor=2.0, + ) + np.testing.assert_allclose(cartesian, [[46.121098, 0.0, 0.0]]) + np.testing.assert_allclose(precession, force_mag * (2.0 / HBAR_METAL)) + + +# --------------------------------------------------------------------------- +# Behavioral coverage through the real ``dpa4spin`` adapter +# --------------------------------------------------------------------------- + + +def _resolve_model() -> Path: + env = os.environ.get(_MODEL_ENV) + return Path(env) if env else _DEFAULT_MODEL + + +def _require_dpa4spin_host(): + """Import PyLammps and require a host ``dpa4spin`` registration.""" + try: + from lammps import ( + lammps, + ) + except ImportError: + pytest.skip("PyLammps is not available") + model = _resolve_model() + if not model.is_file(): + pytest.skip( + f"native-spin dpa4spin fixture missing ({model}); " + f"set {_MODEL_ENV} or run source/tests/infer/gen_dpa4c_spin_lmp.py" + ) + probe = lammps(cmdargs=["-log", "none", "-screen", "none"]) + try: + if not probe.has_style("pair", "dpa4spin"): + pytest.skip("pair style dpa4spin is not registered in this LAMMPS build") + finally: + probe.close() + return lammps, model + + +def _extract_fm_cartesian_force(lmp) -> int: + """Read the scalar int from ``Pair::extract("fm_cartesian_force", dim)``.""" + extract = lmp.lib.lammps_extract_pair + # The Python binding defaults restype to c_int, which truncates pointers. + previous = extract.restype + extract.restype = ctypes.c_void_p + try: + c_ptr = extract(lmp.lmp, b"fm_cartesian_force") + finally: + extract.restype = previous + if not c_ptr: + raise AssertionError("extract(fm_cartesian_force) returned null") + return int(ctypes.cast(c_ptr, ctypes.POINTER(ctypes.c_int)).contents.value) + + +def _open_dpa4spin(lammps_cls, data_file: Path, model: Path, pair_args: str): + """Build a single-rank metal/spin system with ``pair_style dpa4spin``.""" + lmp = lammps_cls(cmdargs=["-log", "none", "-screen", "none"]) + try: + for cmd in ( + "units metal", + "dimension 3", + "boundary p p p", + "atom_style spin", + "atom_modify map yes", + f"read_data {data_file.resolve()}", + "mass 1 58.6934", + "mass 2 15.9994", + f"pair_style dpa4spin {model.resolve()} {pair_args}".rstrip(), + "pair_coeff * * Ni O", + "neighbor 2.0 bin", + "neigh_modify every 1 delay 0 check yes", + "compute fmprop all property/atom fmx fmy fmz", + "compute spprop all property/atom spx spy spz sp", + "run 0", + ): + lmp.command(cmd) + except Exception: + lmp.close() + raise + return lmp + + +def _gather_fm_sp(lmp) -> tuple[np.ndarray, np.ndarray]: + natoms = int(lmp.get_natoms()) + fm = np.array(lmp.gather("c_fmprop", 1, 3), dtype=np.float64).reshape(natoms, 3) + sp = np.array(lmp.gather("c_spprop", 1, 4), dtype=np.float64).reshape(natoms, 4) + return fm, sp + + +@pytest.fixture(scope="module") +def dpa4spin_data(tmp_path_factory: pytest.TempPathFactory) -> Path: + """Write the mixed zero-/nonzero-moment NiO data file once per module.""" + from write_lmp_data import ( + write_lmp_data_spin, + ) + + path = tmp_path_factory.mktemp("dpa4spin_mf") / "nio_zero_and_finite.lmp" + write_lmp_data_spin(_BOX, _COORD, _SPIN, _TYPE, path) + return path + + +def test_dpa4spin_extract_reports_magnetic_force_selection( + dpa4spin_data: Path, +) -> None: + """``magnetic_force`` selection must be visible via Pair::extract.""" + lammps_cls, model = _require_dpa4spin_host() + for args, expect in ( + ("", 0), + ("magnetic_force precession", 0), + ("magnetic_force cartesian", 1), + ): + lmp = _open_dpa4spin(lammps_cls, dpa4spin_data, model, args) + try: + assert _extract_fm_cartesian_force(lmp) == expect + finally: + lmp.close() + + +def test_dpa4spin_cartesian_fm_matches_energy_fd_and_keeps_zero_moment( + dpa4spin_data: Path, +) -> None: + """Cartesian ``fm`` is ``-dE/dS`` and stays nonzero at zero moment on Ni.""" + lammps_cls, model = _require_dpa4spin_host() + lmp = _open_dpa4spin(lammps_cls, dpa4spin_data, model, "magnetic_force cartesian") + try: + fm, sp = _gather_fm_sp(lmp) + assert sp[0, 3] == pytest.approx(0.0) + assert np.linalg.norm(fm[0]) > 1e-8, ( + "spin-active Ni at zero moment must keep a nonzero Cartesian fm" + ) + assert sp[1, 3] > 0.0 + assert np.linalg.norm(fm[1]) > 1e-8 + + # Energy FD on the largest |fm| component of Ni atom 2 (1-based id 2). + # ``set atom ... spin`` takes magnitude then direction (mag > 0). + s0 = sp[1, :3] * sp[1, 3] + eps = 1e-5 + comp = int(np.argmax(np.abs(fm[1]))) + analytic0 = float(fm[1, comp]) + assert abs(analytic0) > 1e-8, ( + f"largest |fm[1]| component must be live for FD: fm[1]={fm[1]}" + ) + ehat = np.zeros(3, dtype=np.float64) + ehat[comp] = 1.0 + + def energy_at(moment: np.ndarray) -> float: + mag = float(np.linalg.norm(moment)) + assert mag > 0.0 + u = moment / mag + lmp.command( + f"set atom 2 spin {mag:.16e} {u[0]:.16e} {u[1]:.16e} {u[2]:.16e}" + ) + lmp.command("run 0") + return float(lmp.get_thermo("pe")) + + ep = energy_at(s0 + eps * ehat) + em = energy_at(s0 - eps * ehat) + energy_at(s0) + fm_restored, _ = _gather_fm_sp(lmp) + fd = -(ep - em) / (2.0 * eps) + analytic = float(fm_restored[1, comp]) + assert abs(analytic) > 1e-8, ( + f"analytic fm[{comp}] vanished after restore: analytic={analytic}" + ) + # Relative only: an absolute escape of ~1e-8 would let a scatter + # sign flip (+= -> -=) pass when |fm| is tiny; with |analytic| > 1e-8, + # |analytic - (-analytic)| = 2|analytic| fails both rel and a tight abs. + rel = abs(fd - analytic) / max(1e-12, abs(analytic) + abs(fd)) + assert rel < 5e-3, ( + f"cartesian fm[{comp}] vs energy FD: analytic={analytic}, fd={fd}, rel={rel}" + ) + finally: + lmp.close() + + +def test_dpa4spin_default_precession_scales_cartesian_by_moment_over_hbar( + dpa4spin_data: Path, +) -> None: + """Default precession ``fm`` equals cartesian ``fm`` times ``|m| / hbar``.""" + lammps_cls, model = _require_dpa4spin_host() + lmp_c = _open_dpa4spin(lammps_cls, dpa4spin_data, model, "magnetic_force cartesian") + try: + fm_c, sp = _gather_fm_sp(lmp_c) + finally: + lmp_c.close() + + lmp_p = _open_dpa4spin(lammps_cls, dpa4spin_data, model, "") + try: + assert _extract_fm_cartesian_force(lmp_p) == 0 + fm_p, _ = _gather_fm_sp(lmp_p) + finally: + lmp_p.close() + + # Zero-moment spin-active Ni: precession must vanish; cartesian must not. + np.testing.assert_allclose(fm_p[0], 0.0, atol=1e-14) + assert np.linalg.norm(fm_c[0]) > 1e-8 + + # Nonzero Ni moments: ratio must be |m|/hbar. + for i in (1, 2): + m = sp[i, 3] + expect = m / HBAR_METAL + for d in range(3): + if abs(fm_c[i, d]) < 1e-18: + assert abs(fm_p[i, d]) < 1e-14 + continue + ratio = fm_p[i, d] / fm_c[i, d] + np.testing.assert_allclose(ratio, expect, rtol=1e-8, atol=0.0) + + # Non-magnetic O: both conventions stay at zero. + np.testing.assert_allclose(fm_c[3:], 0.0, atol=1e-14) + np.testing.assert_allclose(fm_p[3:], 0.0, atol=1e-14) + + +def test_dpa4spin_kokkos_cartesian_skipped_without_gpu_backend() -> None: + """Exercise ``dpa4spin/kk`` when a GPU Kokkos build is present; else skip.""" + try: + from lammps import ( + lammps, + ) + except ImportError: + pytest.skip("PyLammps is not available") + + model = _resolve_model() + if not model.is_file(): + pytest.skip(f"fixture missing: {model}") + + probe = lammps(cmdargs=["-log", "none", "-screen", "none"]) + try: + has_kk = probe.has_style("pair", "dpa4spin/kk") + has_kokkos = probe.has_package("KOKKOS") + finally: + probe.close() + if not has_kk or not has_kokkos: + pytest.skip("dpa4spin/kk or KOKKOS package not available") + + # dpa4spin/kk rejects host Kokkos backends and graph-lower (uncompressed) + # archives; a GPU + compressed canonical model is required. + try: + lmp = lammps( + cmdargs=[ + "-k", + "on", + "g", + "1", + "-sf", + "kk", + "-log", + "none", + "-screen", + "none", + ] + ) + except Exception as exc: + pytest.skip(f"GPU Kokkos backend unavailable: {exc}") + + try: + if not lmp.has_style("pair", "dpa4spin/kk"): + pytest.skip("dpa4spin/kk not registered under Kokkos") + # Probe whether this archive is accepted; graph-lower must skip. + from write_lmp_data import ( + write_lmp_data_spin, + ) + + with tempfile.TemporaryDirectory() as tmp: + data = Path(tmp) / "nio.lmp" + write_lmp_data_spin(_BOX, _COORD, _SPIN, _TYPE, data) + try: + for cmd in ( + "units metal", + "dimension 3", + "boundary p p p", + "atom_style spin", + "atom_modify map yes", + f"read_data {data.resolve()}", + "mass 1 58.6934", + "mass 2 15.9994", + f"pair_style dpa4spin/kk {model.resolve()} " + "magnetic_force cartesian", + "pair_coeff * * Ni O", + "neighbor 2.0 bin", + "neigh_modify every 1 delay 0 check yes", + "compute fmprop all property/atom fmx fmy fmz", + "run 0", + ): + lmp.command(cmd) + except Exception as exc: + msg = str(exc) + if "GPU backend only" in msg or "compact canonical" in msg: + pytest.skip(f"dpa4spin/kk cannot run here: {msg}") + raise + assert _extract_fm_cartesian_force(lmp) == 1 + natoms = int(lmp.get_natoms()) + fm = np.array(lmp.gather("c_fmprop", 1, 3), dtype=np.float64).reshape( + natoms, 3 + ) + assert np.linalg.norm(fm[0]) > 1e-8 + assert np.linalg.norm(fm[1]) > 1e-8 + finally: + lmp.close() diff --git a/source/tests/infer/gen_dpa4c_spin_lmp.py b/source/tests/infer/gen_dpa4c_spin_lmp.py new file mode 100644 index 0000000000..e0b2fecf0d --- /dev/null +++ b/source/tests/infer/gen_dpa4c_spin_lmp.py @@ -0,0 +1,93 @@ +#!/usr/bin/env python3 +# SPDX-License-Identifier: LGPL-3.0-or-later +"""Generate a minimal native-spin DPA4C graph ``.pt2`` for LAMMPS ``dpa4spin``. + +Unlike ``gen_dpa4c_spin.py``, this archive carries no fparam, aparam, or +charge/spin condition so ``pair_style dpa4spin`` can load it without extra +keywords. Used by ``source/lmp/tests/test_dpa4spin_magnetic_force_convention.py``. +""" + +from __future__ import ( + annotations, +) + +import copy +import json +import os +import sys +import zipfile + +import numpy as np + +sys.path.insert(0, os.path.join(os.path.dirname(__file__), "..")) + +from dpa4_fixtures import ( + jitter_zero_arrays, +) +from gen_common import ( + ensure_inductor_compiler, + load_custom_ops, +) + +CFG = { + "type_map": ["Ni", "O"], + "descriptor": { + "type": "dpa4c", + "rcut": 4.0, + "channels": 8, + "lmax": 2, + "n_radial": 4, + "precision": "float64", + "seed": 7, + }, + "fitting_net": {"neuron": [8, 8], "precision": "float64", "seed": 7}, + "spin": {"use_spin": [True, False], "scheme": "native"}, +} + +_JITTER_SEED = 20260817 +OUT = os.path.join(os.path.dirname(__file__), "deeppot_dpa4c_spin_lmp.pt2") + + +def main() -> None: + from deepmd.dpmodel.model.model import ( + get_model, + ) + from deepmd.pt_expt.utils.serialization import ( + deserialize_to_file, + ) + + ensure_inductor_compiler() + load_custom_ops() + model = get_model(copy.deepcopy(CFG)) + md = jitter_zero_arrays(model.serialize(), np.random.default_rng(_JITTER_SEED)) + data = { + "model": md, + "model_def_script": CFG, + "backend": "dpmodel", + "software": "deepmd-kit", + "version": "3.0.0", + } + print(f"Exporting {OUT}") # noqa: T201 + deserialize_to_file(OUT, data, do_atomic_virial=True, lower_kind="graph") + with zipfile.ZipFile(OUT) as zf: + meta = json.loads(zf.read("model/extra/metadata.json").decode()) + print( # noqa: T201 + json.dumps( + { + k: meta.get(k) + for k in ( + "spin_scheme", + "has_comm_artifact", + "dim_fparam", + "dim_aparam", + "dim_chg_spin", + "lower_input_kind", + ) + }, + indent=2, + ) + ) + + +if __name__ == "__main__": + main()