diff --git a/doc/third-party/lammps-command.md b/doc/third-party/lammps-command.md index 70df75d22d..dfb88b5fa4 100644 --- a/doc/third-party/lammps-command.md +++ b/doc/third-party/lammps-command.md @@ -49,7 +49,7 @@ pair_style deepmd models ... keyword value ... - models = frozen model(s) to compute the interaction. If multiple models are provided, then only the first model serves to provide energy and force prediction for each timestep of molecular dynamics, and the model deviation will be computed among all models every `out_freq` timesteps. -- keyword = _out_file_ or _out_freq_ or _fparam_ or _fparam_from_compute_ or _aparam_from_compute_ or _atomic_ or _relative_ or _relative_v_ or _aparam_ or _ttm_ +- keyword = _out_file_ or _out_freq_ or _fparam_ or _fparam_from_compute_ or _fparam_from_fix_ or _aparam_from_compute_ or _atomic_ or _relative_ or _relative_v_ or _aparam_ or _ttm_ 
     out_file value = filename
@@ -60,6 +60,9 @@ pair_style deepmd models ... keyword value ...
         parameters = one or more frame parameters required for model evaluation.
     fparam_from_compute value = id
         id = compute id used to update the frame parameter.
+    fparam_from_fix value = id [index]
+        id = fix ID used to update the frame parameter.
+        index = optional 1-based starting index for a global fix vector.
     aparam_from_compute value = id
         id = compute id used to update the atom parameter.
     atomic = no value is required.
@@ -110,6 +113,7 @@ $$E_{v_i}=\frac{\left|D_{v_i}\right|}{\left|v_i\right|+l}$$
 
 If the keyword `fparam` is set, the given frame parameter(s) will be fed to the model.
 If the keyword `fparam_from_compute` is set, the global parameter(s) from compute command (e.g., temperature from [compute temp command](https://docs.lammps.org/compute_temp.html)) will be fed to the model as the frame parameter(s).
+If the keyword `fparam_from_fix` is set, the global parameter(s) from fix command will be fed to the model as the frame parameter(s). This is intended for generalized coordinates that are naturally carried by a fix, for example the potentiostat variable in `fix uvt`. See [compute `deepmd/fparam/dedn`](#compute-deepmdfparamdedn) for evaluating the corresponding energy derivative.
 If the keyword `aparam_from_compute` is set, the atomic parameter(s) from compute command (e.g., per-atom translational kinetic energy from [compute ke/atom command](https://docs.lammps.org/compute_ke_atom.html)) will be fed to the model as the atom parameter(s).
 If the keyword `aparam` is set, the given atomic parameter(s) will be fed to the model, where each atom is assumed to have the same atomic parameter(s).
 If the keyword `ttm` is set, electronic temperatures from [fix ttm command](https://docs.lammps.org/fix_ttm.html) will be fed to the model as the atomic parameters.
@@ -198,6 +202,23 @@ Please note that the virial and atomic virial are not currently supported in spi
 
 - The `deepspin` pair style is provided in the USER-DEEPMD package, which is compiled from the DeePMD-kit, visit the [DeePMD-kit website](https://github.com/deepmodeling/deepmd-kit) for more information.
 
+## compute `deepmd/fparam/dedn`
+
+The `deepmd/fparam/dedn` compute evaluates a finite-difference derivative of the model energy with respect to a chosen frame parameter source.
+
+```lammps
+compute ID group-ID deepmd/fparam/dedn source [delta]
+```
+
+- `source` can be a global variable (`v_name`), a compute (`c_ID`) or a fix (`f_ID`).
+- `source[index]` may be used for vector-valued computes or fixes, with 1-based indexing.
+- `delta` is the perturbation used in the central difference formula. If omitted, a small default perturbation is used.
+- The compute performs two additional model-energy evaluations, at `source + delta`
+  and `source - delta`. It does not consume a direct derivative tensor from the model.
+- This path currently requires `pair_style deepmd` with one model and one frame-parameter
+  dimension. It therefore works with existing supported backends and models that accept
+  a scalar frame parameter; no `o_dE_dN` or other derivative output is required.
+
 ## Compute tensorial properties
 
 The DeePMD-kit package provides the compute `deeptensor/atom` for computing atomic tensorial properties.
diff --git a/source/lmp/compute_deepmd_fparam_dedn.cpp b/source/lmp/compute_deepmd_fparam_dedn.cpp
new file mode 100644
index 0000000000..88ee8825b8
--- /dev/null
+++ b/source/lmp/compute_deepmd_fparam_dedn.cpp
@@ -0,0 +1,202 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+#include "compute_deepmd_fparam_dedn.h"
+
+#include 
+#include 
+#include 
+#include 
+#include 
+
+#include "comm.h"
+#include "compute.h"
+#include "error.h"
+#include "fix.h"
+#include "force.h"
+#include "input.h"
+#include "modify.h"
+#include "update.h"
+#include "variable.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+ComputeDeepmdFparamDedn::ComputeDeepmdFparamDedn(LAMMPS* lmp,
+                                                 int narg,
+                                                 char** arg)
+    : Compute(lmp, narg, arg), source_index(-1), delta(1.0e-6), pair(nullptr) {
+  if (narg < 4) {
+    error->all(FLERR, "Illegal compute deepmd/fparam/dedn command");
+  }
+
+  std::string token = arg[3];
+  auto lb = token.find('[');
+  auto rb = token.find(']');
+  if (lb != std::string::npos || rb != std::string::npos) {
+    if (lb == std::string::npos || rb == std::string::npos || rb <= lb + 1 ||
+        rb != token.size() - 1) {
+      error->all(FLERR, "Illegal source specification in compute command");
+    }
+    std::string idx = token.substr(lb + 1, rb - lb - 1);
+    char* endptr = nullptr;
+    errno = 0;
+    long one_based = std::strtol(idx.c_str(), &endptr, 10);
+    if (endptr == idx.c_str() || *endptr != '\0' || errno == ERANGE ||
+        one_based < 1 ||
+        one_based > static_cast(std::numeric_limits::max()) + 1L) {
+      error->all(FLERR, "Source index must be a positive 1-based integer");
+    }
+    source_index = static_cast(one_based - 1);
+    token = token.substr(0, lb);
+  }
+
+  if (token.rfind("v_", 0) == 0) {
+    source_type = SRC_VAR;
+    source_id = token.substr(2);
+    if (source_index >= 0) {
+      error->all(FLERR,
+                 "Variable source for compute deepmd/fparam/dedn must be "
+                 "scalar");
+    }
+  } else if (token.rfind("c_", 0) == 0) {
+    source_type = SRC_COMPUTE;
+    source_id = token.substr(2);
+  } else if (token.rfind("f_", 0) == 0) {
+    source_type = SRC_FIX;
+    source_id = token.substr(2);
+  } else {
+    error->all(FLERR, "Source must be a variable, compute, or fix reference");
+  }
+  if (source_id.empty()) {
+    error->all(FLERR, "Source ID must not be empty");
+  }
+
+  int iarg = 4;
+  if (iarg < narg) {
+    char* endptr = nullptr;
+    errno = 0;
+    delta = std::strtod(arg[iarg], &endptr);
+    if (endptr == arg[iarg] || *endptr != '\0' || errno == ERANGE ||
+        !std::isfinite(delta)) {
+      error->all(FLERR,
+                 "delta must be a finite number in compute deepmd/fparam/dedn");
+    }
+    ++iarg;
+  }
+  if (iarg != narg) {
+    error->all(FLERR, "Illegal compute deepmd/fparam/dedn command");
+  }
+  if (delta <= 0.0) {
+    error->all(FLERR, "delta must be > 0 in compute deepmd/fparam/dedn");
+  }
+
+  scalar_flag = 1;
+  extscalar = 1;
+  timeflag = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeDeepmdFparamDedn::~ComputeDeepmdFparamDedn() = default;
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeDeepmdFparamDedn::init() {
+  if (!force->pair) {
+    error->all(FLERR,
+               "compute deepmd/fparam/dedn requires an active pair style");
+  }
+  pair = dynamic_cast(force->pair);
+  if (!pair) {
+    error->all(
+        FLERR,
+        "compute deepmd/fparam/dedn currently requires pair_style deepmd");
+  }
+  if (pair->get_dim_fparam() != 1) {
+    error->all(FLERR,
+               "compute deepmd/fparam/dedn currently supports a single "
+               "frame-parameter dimension only");
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+double ComputeDeepmdFparamDedn::get_source_value() {
+  if (source_type == SRC_VAR) {
+    int ivar = input->variable->find(source_id.c_str());
+    if (ivar < 0) {
+      error->all(FLERR, "Variable source not found: " + source_id);
+    }
+    return input->variable->compute_equal(ivar);
+  }
+
+  if (source_type == SRC_COMPUTE) {
+    int icompute = modify->find_compute(source_id);
+    if (icompute < 0) {
+      error->all(FLERR, "Compute source not found: " + source_id);
+    }
+    Compute* compute = modify->compute[icompute];
+    if (!compute) {
+      error->all(FLERR, "Compute source not found: " + source_id);
+    }
+    if (source_index < 0) {
+      if (!compute->scalar_flag) {
+        error->all(FLERR, "Compute source is not scalar: " + source_id);
+      }
+      if (!(compute->invoked_flag & Compute::INVOKED_SCALAR)) {
+        compute->compute_scalar();
+        compute->invoked_flag |= Compute::INVOKED_SCALAR;
+      }
+      return compute->scalar;
+    }
+    if (!compute->vector_flag) {
+      error->all(FLERR, "Compute source is not vector-valued: " + source_id);
+    }
+    if (!(compute->invoked_flag & Compute::INVOKED_VECTOR)) {
+      compute->compute_vector();
+      compute->invoked_flag |= Compute::INVOKED_VECTOR;
+    }
+    if (source_index >= compute->size_vector) {
+      error->all(FLERR, "Compute source index is out of range: " + source_id);
+    }
+    return compute->vector[source_index];
+  }
+
+  int ifix = modify->find_fix(source_id);
+  if (ifix < 0) {
+    error->all(FLERR, "Fix source not found: " + source_id);
+  }
+  Fix* fix = modify->fix[ifix];
+  if (!fix) {
+    error->all(FLERR, "Fix source not found: " + source_id);
+  }
+  if (source_index < 0) {
+    if (!fix->scalar_flag) {
+      error->all(FLERR, "Fix source is not scalar: " + source_id);
+    }
+    return fix->compute_scalar();
+  }
+  if (!fix->vector_flag) {
+    error->all(FLERR, "Fix source is not vector-valued: " + source_id);
+  }
+  if (fix->size_vector <= source_index) {
+    error->all(FLERR, "Fix source index is out of range: " + source_id);
+  }
+  return fix->compute_vector(source_index);
+}
+
+/* ---------------------------------------------------------------------- */
+
+double ComputeDeepmdFparamDedn::compute_scalar() {
+  invoked_scalar = update->ntimestep;
+  double fparam0 = get_source_value();
+  std::vector fparam_plus(1, fparam0 + delta);
+  std::vector fparam_minus(1, fparam0 - delta);
+
+  double one = (pair->eval_energy_with_fparam(fparam_plus) -
+                pair->eval_energy_with_fparam(fparam_minus)) /
+               (2.0 * delta);
+
+  MPI_Allreduce(&one, &scalar, 1, MPI_DOUBLE, MPI_SUM, world);
+  return scalar;
+}
diff --git a/source/lmp/compute_deepmd_fparam_dedn.h b/source/lmp/compute_deepmd_fparam_dedn.h
new file mode 100644
index 0000000000..66f4808031
--- /dev/null
+++ b/source/lmp/compute_deepmd_fparam_dedn.h
@@ -0,0 +1,38 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+#ifdef COMPUTE_CLASS
+// clang-format off
+ComputeStyle(deepmd/fparam/dedn, ComputeDeepmdFparamDedn)
+// clang-format on
+#else
+
+#ifndef LMP_COMPUTE_DEEPMD_FPARAM_DEDN_H
+#define LMP_COMPUTE_DEEPMD_FPARAM_DEDN_H
+
+#include "compute.h"
+#include "pair_deepmd.h"
+
+namespace LAMMPS_NS {
+
+class ComputeDeepmdFparamDedn : public Compute {
+ public:
+  ComputeDeepmdFparamDedn(class LAMMPS*, int, char**);
+  ~ComputeDeepmdFparamDedn() override;
+  void init() override;
+  double compute_scalar() override;
+
+ private:
+  enum SourceType { SRC_VAR, SRC_COMPUTE, SRC_FIX };
+
+  SourceType source_type;
+  std::string source_id;
+  int source_index;
+  double delta;
+  PairDeepMD* pair;
+
+  double get_source_value();
+};
+
+}  // namespace LAMMPS_NS
+
+#endif
+#endif
diff --git a/source/lmp/pair_base.cpp b/source/lmp/pair_base.cpp
index 63e5462590..2639bc74ba 100644
--- a/source/lmp/pair_base.cpp
+++ b/source/lmp/pair_base.cpp
@@ -164,6 +164,52 @@ void PairDeepBaseModel::make_fparam_from_compute(vector& fparam) {
   }
 }
 
+void PairDeepBaseModel::make_fparam_from_fix(vector& fparam) {
+  assert(do_fix_fparam);
+
+  int ifix = modify->find_fix(fix_fparam_id);
+  if (ifix < 0) {
+    error->all(FLERR, "fix id is not found: " + fix_fparam_id);
+  }
+  Fix* fix = modify->fix[ifix];
+
+  if (!fix) {
+    error->all(FLERR, "fix id is not found: " + fix_fparam_id);
+  }
+  fparam.resize(dim_fparam);
+
+  if (fix_fparam_index < 0) {
+    if (dim_fparam == 1) {
+      if (!fix->scalar_flag) {
+        error->all(FLERR, "fix " + fix_fparam_id +
+                              " does not provide a scalar for fparam");
+      }
+      fparam[0] = fix->compute_scalar();
+    } else if (dim_fparam > 1) {
+      if (!fix->scalar_flag) {
+        error->all(FLERR, "fix " + fix_fparam_id +
+                              " does not provide a scalar for fparam");
+      }
+      double value = fix->compute_scalar();
+      for (int jj = 0; jj < dim_fparam; ++jj) {
+        fparam[jj] = value;
+      }
+    }
+  } else {
+    if (!fix->vector_flag) {
+      error->all(FLERR, "fix " + fix_fparam_id +
+                            " does not provide a vector for fparam");
+    }
+    if (fix_fparam_index > fix->size_vector - dim_fparam) {
+      error->all(FLERR, "fix " + fix_fparam_id +
+                            " vector is shorter than fparam dimension");
+    }
+    for (int jj = 0; jj < dim_fparam; ++jj) {
+      fparam[jj] = fix->compute_vector(fix_fparam_index + jj);
+    }
+  }
+}
+
 void PairDeepBaseModel::make_aparam_from_compute(vector& aparam) {
   assert(do_compute_aparam);
 
@@ -338,6 +384,8 @@ PairDeepBaseModel::PairDeepBaseModel(
   scale = NULL;
   do_ttm = false;
   do_compute_fparam = false;
+  do_fix_fparam = false;
+  fix_fparam_index = -1;
   do_compute_aparam = false;
   single_model = false;
   multi_models_mod_devi = false;
diff --git a/source/lmp/pair_base.h b/source/lmp/pair_base.h
index 1dd4b84041..f7c3b34808 100644
--- a/source/lmp/pair_base.h
+++ b/source/lmp/pair_base.h
@@ -43,6 +43,7 @@ class PairDeepBaseModel : public Pair {
   void print_summary(const std::string pre) const;
   int get_node_rank();
   void cum_sum(std::map&, std::map&);
+  int get_dim_fparam() const { return dim_fparam; }
 
   std::string get_file_content(const std::string& model);
   std::vector get_file_content(
@@ -86,6 +87,10 @@ class PairDeepBaseModel : public Pair {
   void make_fparam_from_compute(std::vector& fparam);
   bool do_compute_fparam;
   std::string compute_fparam_id;
+  void make_fparam_from_fix(std::vector& fparam);
+  bool do_fix_fparam;
+  std::string fix_fparam_id;
+  int fix_fparam_index;
   void make_aparam_from_compute(std::vector& aparam);
   bool do_compute_aparam;
   std::string compute_aparam_id;
diff --git a/source/lmp/pair_deepmd.cpp b/source/lmp/pair_deepmd.cpp
index fb27fa0ff6..d60388814f 100644
--- a/source/lmp/pair_deepmd.cpp
+++ b/source/lmp/pair_deepmd.cpp
@@ -2,6 +2,8 @@
 #include 
 
 #include 
+#include 
+#include 
 #include 
 #include 
 #include 
@@ -119,7 +121,8 @@ static const char cite_user_deepmd_package[] =
 
 PairDeepMD::PairDeepMD(LAMMPS* lmp)
     : PairDeepBaseModel(
-          lmp, cite_user_deepmd_package, deep_pot, deep_pot_model_devi) {
+          lmp, cite_user_deepmd_package, deep_pot, deep_pot_model_devi),
+      commdata_(nullptr) {
   // Constructor body can be empty
 }
 
@@ -127,6 +130,112 @@ PairDeepMD::~PairDeepMD() {
   // Ensure base class destructor is called
 }
 
+double PairDeepMD::eval_energy_with_fparam(
+    const std::vector& fparam_override) {
+  if (numb_models != 1) {
+    error->all(FLERR,
+               "deepmd/fparam/dedn currently supports single-model pair_style "
+               "only");
+  }
+  if (atom->sp_flag) {
+    error->all(FLERR,
+               "Pair style 'deepmd' does not support spin atoms, please use "
+               "pair style 'deepspin' instead.");
+  }
+
+  bool do_ghost = true;
+  commdata_ = (CommBrickDeepMD*)comm;
+  double** x = atom->x;
+  int* type = atom->type;
+  int nlocal = atom->nlocal;
+  int nghost = 0;
+  if (do_ghost) {
+    nghost = atom->nghost;
+  }
+  int nall = nlocal + nghost;
+
+  std::vector dtype(nall);
+  for (int ii = 0; ii < nall; ++ii) {
+    dtype[ii] = type_idx_map[type[ii] - 1];
+  }
+
+  double dener(0);
+  std::vector dforce(nall * 3);
+  std::vector dvirial(9, 0);
+  std::vector dcoord(nall * 3, 0.);
+  std::vector dbox(9, 0);
+  std::vector daparam;
+
+  if (fparam_override.size() != static_cast(dim_fparam)) {
+    error->all(FLERR, "fparam override has the wrong dimension");
+  }
+
+  // get box
+  dbox[0] = domain->h[0] / dist_unit_cvt_factor;  // xx
+  dbox[4] = domain->h[1] / dist_unit_cvt_factor;  // yy
+  dbox[8] = domain->h[2] / dist_unit_cvt_factor;  // zz
+  dbox[7] = domain->h[3] / dist_unit_cvt_factor;  // zy
+  dbox[6] = domain->h[4] / dist_unit_cvt_factor;  // zx
+  dbox[3] = domain->h[5] / dist_unit_cvt_factor;  // yx
+
+  // get coord
+  for (int ii = 0; ii < nall; ++ii) {
+    for (int dd = 0; dd < 3; ++dd) {
+      dcoord[ii * 3 + dd] =
+          (x[ii][dd] - domain->boxlo[dd]) / dist_unit_cvt_factor;
+    }
+  }
+
+  // mapping (for DPA-2/3 .pt2 GNN models that gather ghost features via
+  // the LAMMPS atom-map; harmless for other models).
+  std::vector mapping_vec(nall, -1);
+  if (comm->nprocs == 1 && atom->map_style != Atom::MAP_NONE) {
+    for (size_t ii = 0; ii < nall; ++ii) {
+      mapping_vec[ii] = atom->map(atom->tag[ii]);
+    }
+  }
+
+  if (do_compute_aparam) {
+    make_aparam_from_compute(daparam);
+  } else if (aparam.size() > 0) {
+    make_uniform_aparam(daparam, aparam, nlocal);
+  } else if (do_ttm) {
+#ifdef USE_TTM
+    if (dim_aparam > 0) {
+      make_ttm_aparam(daparam);
+    }
+#endif
+  }
+  int ago = neighbor->ago;
+
+  if (do_ghost) {
+    if (!list) {
+      error->all(FLERR,
+                 "deepmd/fparam/dedn requires an available pair neighbor list");
+    }
+    deepmd_compat::InputNlist lmp_list(
+        list->inum, list->ilist, list->numneigh, list->firstneigh,
+        commdata_->nswap, commdata_->sendnum, commdata_->recvnum,
+        commdata_->firstrecv, commdata_->sendlist, commdata_->sendproc,
+        commdata_->recvproc, &world, comm->nprocs);
+    lmp_list.set_mask(NEIGHMASK);
+    if (comm->nprocs == 1 && atom->map_style != Atom::MAP_NONE) {
+      lmp_list.set_mapping(mapping_vec.data());
+    }
+
+    try {
+      deep_pot.compute(dener, dforce, dvirial, dcoord, dtype, dbox, nghost,
+                       lmp_list, ago, fparam_override, daparam);
+    } catch (deepmd_compat::deepmd_exception& e) {
+      error->one(FLERR, e.what());
+    }
+  } else {
+    error->all(FLERR, "unknown computational branch");
+  }
+
+  return scale[1][1] * dener * ener_unit_cvt_factor;
+}
+
 void PairDeepMD::compute(int eflag, int vflag) {
   if (numb_models == 0) {
     return;
@@ -214,6 +323,8 @@ void PairDeepMD::compute(int eflag, int vflag) {
 
   if (do_compute_fparam) {
     make_fparam_from_compute(fparam);
+  } else if (do_fix_fparam) {
+    make_fparam_from_fix(fparam);
   }
 
   // int ago = numb_models > 1 ? 0 : neighbor->ago;
@@ -533,6 +644,7 @@ static bool is_key(const string& input) {
   keys.push_back("fparam");
   keys.push_back("aparam");
   keys.push_back("fparam_from_compute");
+  keys.push_back("fparam_from_fix");
   keys.push_back("aparam_from_compute");
   keys.push_back("ttm");
   keys.push_back("atomic");
@@ -675,6 +787,31 @@ void PairDeepMD::settings(int narg, char** arg) {
       do_compute_fparam = true;
       compute_fparam_id = arg[iarg + 1];
       iarg += 1 + 1;
+    } else if (string(arg[iarg]) == string("fparam_from_fix")) {
+      if (iarg + 1 >= narg || is_key(arg[iarg + 1])) {
+        error->all(FLERR,
+                   "invalid fparam_from_fix key: should be "
+                   "fparam_from_fix fix_fparam_id(str) [fix_vector_index]");
+      }
+      do_fix_fparam = true;
+      fix_fparam_id = arg[iarg + 1];
+      fix_fparam_index = -1;
+      if (iarg + 2 < narg && !is_key(arg[iarg + 2])) {
+        char* endptr = nullptr;
+        errno = 0;
+        long one_based = std::strtol(arg[iarg + 2], &endptr, 10);
+        if (endptr == arg[iarg + 2] || *endptr != '\0' || errno == ERANGE ||
+            one_based < 1 ||
+            one_based > static_cast(std::numeric_limits::max())) {
+          error->all(FLERR,
+                     "invalid fparam_from_fix key: vector index must be a "
+                     "positive 1-based integer");
+        }
+        fix_fparam_index = static_cast(one_based - 1);
+        iarg += 3;
+      } else {
+        iarg += 2;
+      }
     } else if (string(arg[iarg]) == string("aparam_from_compute")) {
       for (int ii = 0; ii < 1; ++ii) {
         if (iarg + 1 + ii >= narg || is_key(arg[iarg + 1 + ii])) {
@@ -725,6 +862,15 @@ void PairDeepMD::settings(int narg, char** arg) {
         FLERR,
         "fparam and fparam_from_compute should NOT be set simultaneously");
   }
+  if (do_fix_fparam && fparam.size() > 0) {
+    error->all(FLERR,
+               "fparam and fparam_from_fix should NOT be set simultaneously");
+  }
+  if (do_fix_fparam && do_compute_fparam) {
+    error->all(FLERR,
+               "fparam_from_compute and fparam_from_fix should NOT be set "
+               "simultaneously");
+  }
 
   if (comm->me == 0) {
     if (numb_models > 1 && out_freq > 0) {
@@ -769,6 +915,14 @@ void PairDeepMD::settings(int narg, char** arg) {
       cout << pre << "using compute id (fparam):      ";
       cout << compute_fparam_id << "  " << endl;
     }
+    if (do_fix_fparam) {
+      cout << pre << "using fix id (fparam):          ";
+      cout << fix_fparam_id;
+      if (fix_fparam_index >= 0) {
+        cout << "[" << fix_fparam_index + 1 << "]";
+      }
+      cout << "  " << endl;
+    }
     if (do_compute_aparam) {
       cout << pre << "using compute id (aparam):      ";
       cout << compute_aparam_id << "  " << endl;
diff --git a/source/lmp/pair_deepmd.h b/source/lmp/pair_deepmd.h
index 6d54a69fe6..86b2017e93 100644
--- a/source/lmp/pair_deepmd.h
+++ b/source/lmp/pair_deepmd.h
@@ -49,6 +49,7 @@ class PairDeepMD : public PairDeepBaseModel {
   void compute(int, int) override;
   int pack_reverse_comm(int, int, double*) override;
   void unpack_reverse_comm(int, int*, double*) override;
+  double eval_energy_with_fparam(const std::vector& fparam_override);
 
  protected:
   deepmd_compat::DeepPot deep_pot;
diff --git a/source/lmp/plugin/deepmdplugin.cpp b/source/lmp/plugin/deepmdplugin.cpp
index d3b54f8e41..b49b034c41 100644
--- a/source/lmp/plugin/deepmdplugin.cpp
+++ b/source/lmp/plugin/deepmdplugin.cpp
@@ -2,6 +2,7 @@
 /**
  * See https://docs.lammps.org/Developer_plugins.html
  */
+#include "compute_deepmd_fparam_dedn.h"
 #include "compute_deeptensor_atom.h"
 #include "deepmd_version.h"
 #include "fix_dplr.h"
@@ -22,6 +23,10 @@ static Compute* computedeepmdtensoratom(LAMMPS* lmp, int narg, char** arg) {
   return new ComputeDeeptensorAtom(lmp, narg, arg);
 }
 
+static Compute* computedeepmdfparamdedn(LAMMPS* lmp, int narg, char** arg) {
+  return new ComputeDeepmdFparamDedn(lmp, narg, arg);
+}
+
 static Fix* fixdplr(LAMMPS* lmp, int narg, char** arg) {
   return new FixDPLR(lmp, narg, arg);
 }
@@ -59,6 +64,13 @@ extern "C" void lammpsplugin_init(void* lmp, void* handle, void* regfunc) {
   plugin.creator.v2 = (lammpsplugin_factory2*)&computedeepmdtensoratom;
   (*register_plugin)(&plugin, lmp);
 
+  plugin.style = "compute";
+  plugin.name = "deepmd/fparam/dedn";
+  plugin.info = "compute deepmd/fparam/dedn " STR_GIT_SUMM;
+  plugin.author = "Li Fu";
+  plugin.creator.v2 = (lammpsplugin_factory2*)&computedeepmdfparamdedn;
+  (*register_plugin)(&plugin, lmp);
+
   plugin.style = "fix";
   plugin.name = "dplr";
   plugin.info = "fix dplr " STR_GIT_SUMM;
diff --git a/source/lmp/tests/test_lammps_fparam_from_fix_dedn.py b/source/lmp/tests/test_lammps_fparam_from_fix_dedn.py
new file mode 100644
index 0000000000..c929380051
--- /dev/null
+++ b/source/lmp/tests/test_lammps_fparam_from_fix_dedn.py
@@ -0,0 +1,134 @@
+# SPDX-License-Identifier: LGPL-3.0-or-later
+"""Test DeepMD fparam from fix and numeric dE/dfparam from the current state."""
+
+import os
+from pathlib import (
+    Path,
+)
+
+import numpy as np
+import pytest
+from lammps import (
+    PyLammps,
+)
+from model_convert import (
+    ensure_converted_pb,
+)
+from write_lmp_data import (
+    write_lmp_data,
+)
+
+pbtxt_file = (
+    Path(__file__).parent.parent.parent / "tests" / "infer" / "fparam_aparam.pbtxt"
+)
+pb_file = Path(__file__).parent / "fparam_aparam.pb"
+data_file = Path(__file__).parent / "data.lmp"
+
+box = np.array([0, 13, 0, 13, 0, 13, 0, 0, 0])
+coord = np.array(
+    [
+        [12.83, 2.56, 2.18],
+        [12.09, 2.87, 2.74],
+        [0.25, 3.32, 1.68],
+        [3.36, 3.00, 1.81],
+        [3.51, 2.51, 2.60],
+        [4.27, 3.22, 1.56],
+    ]
+)
+type_OH = np.array([1, 1, 1, 1, 1, 1])
+
+
+def setup_module() -> None:
+    """Create the converted model and data file needed by this test module."""
+    if os.environ.get("ENABLE_TENSORFLOW", "1") != "1":
+        pytest.skip("Skip test because TensorFlow support is not enabled.")
+    ensure_converted_pb(pbtxt_file, pb_file)
+    write_lmp_data(box, coord, type_OH, data_file)
+
+
+def teardown_module() -> None:
+    """Remove the temporary LAMMPS data file after the tests finish."""
+    if data_file.exists():
+        os.remove(data_file)
+
+
+def _lammps(fp_value, units="metal") -> PyLammps:
+    """Build a LAMMPS instance configured for frame-parameter derivative tests."""
+    lammps = PyLammps()
+    lammps.units(units)
+    lammps.boundary("p p p")
+    lammps.atom_style("atomic")
+    lammps.neighbor("2.0 bin")
+    lammps.neigh_modify("every 10 delay 0 check no")
+    lammps.read_data(data_file.resolve())
+    lammps.mass("1 16")
+    lammps.timestep(0.0005)
+    lammps.fix("1 all nve")
+    lammps.variable("fp equal " + str(fp_value))
+    lammps.variable("dummy equal 0.0")
+    lammps.fix("fpfix all ave/time 1 1 1 v_dummy v_fp")
+    lammps.pair_style(
+        f"deepmd {pb_file.resolve()} fparam_from_fix fpfix 2 aparam 0.25852028"
+    )
+    lammps.pair_coeff("* *")
+    lammps.compute("dedn all deepmd/fparam/dedn f_fpfix[2]")
+    return lammps
+
+
+@pytest.fixture
+def lammps():
+    """Provide a ready-to-run LAMMPS instance for the default frame parameter."""
+    lmp = _lammps(fp_value=0.25852028)
+    yield lmp
+    lmp.close()
+
+
+def _energy_at_fp(fp_value):
+    """Evaluate the potential energy at a chosen frame-parameter value."""
+    lmp = _lammps(fp_value=fp_value)
+    try:
+        lmp.run(1)
+        return lmp.eval("pe")
+    finally:
+        lmp.close()
+
+
+def _energy_with_direct_fparam(fp_value):
+    """Evaluate the potential energy using a direct fparam setting."""
+    lmp = PyLammps()
+    try:
+        lmp.units("metal")
+        lmp.boundary("p p p")
+        lmp.atom_style("atomic")
+        lmp.neighbor("2.0 bin")
+        lmp.neigh_modify("every 10 delay 0 check no")
+        lmp.read_data(data_file.resolve())
+        lmp.mass("1 16")
+        lmp.timestep(0.0005)
+        lmp.fix("1 all nve")
+        lmp.pair_style(
+            f"deepmd {pb_file.resolve()} fparam {fp_value} aparam 0.25852028"
+        )
+        lmp.pair_coeff("* *")
+        lmp.run(1)
+        return lmp.eval("pe")
+    finally:
+        lmp.close()
+
+
+def test_pair_fparam_from_fix(lammps) -> None:
+    """Check that fparam_from_fix matches the direct fparam path."""
+    lammps.run(1)
+    assert lammps.eval("pe") == pytest.approx(_energy_with_direct_fparam(0.25852028))
+
+
+def test_compute_deepmd_fparam_dedn_without_direct_model_output(lammps) -> None:
+    """Check finite differences with a model that has no direct dE/dN output."""
+    eps = 1.0e-6
+    assert "o_dE_dN" not in pbtxt_file.read_text()
+    lammps.run(1)
+    dedn = lammps.eval("c_dedn")
+    ref = (_energy_at_fp(0.25852028 + eps) - _energy_at_fp(0.25852028 - eps)) / (
+        2.0 * eps
+    )
+    assert dedn == pytest.approx(ref, rel=1.0e-4, abs=1.0e-4)