Enhance SOG KSpace Implementation and Documentation

- Updated sog.h to include new methods and variables for improved functionality, including grid settings and spectral kernel calculations.
- Added a comprehensive developer workflow document for building, validating, and debugging the SOG and LES kspace styles.
- Introduced a detailed derivation document for the SOG discrete Green function on the internal mesh-FFT path, outlining the mathematical framework and implementation specifics.

Author: Johnsonjz

log.lammps

@@ -1,65 +1,67 @@
 LAMMPS (22 Jul 2025 - Update 4)
-  using 1 OpenMP thread(s) per MPI task
-Loaded 0 plugins from /data/zyjin/dp_pt/dp_devel/deepmd-kit-devel/build/sog_lmp_stable/lmp/plugin
-units metal
-atom_style charge
-boundary p p p
-plugin load /data/zyjin/dp_pt/dp_devel/deepmd-kit-devel/build/sog_lmp_stable/lmp/plugin/libdeepmd_lmp.so
+units           metal
+plugin          load /data/zyjin/dp_pt/dp_devel/deepmd-kit-devel/build/sog_lmp_stable/lmp/plugin/libdeepmd_lmp.so
 Loading plugin: deepmd pair style v3.1.3-83-g3cba3c7a-dirty by Han Wang
 Loading plugin: deepspin pair style v3.1.3-83-g3cba3c7a-dirty by Duo Zhang
 Loading plugin: compute deeptensor/atom v3.1.3-83-g3cba3c7a-dirty by Han Wang
 Loading plugin: fix dplr v3.1.3-83-g3cba3c7a-dirty by Han Wang
 Loading plugin: kspace pppm/dplr v3.1.3-83-g3cba3c7a-dirty by Han Wang
 Loading plugin: kspace les v3.1.3-83-g3cba3c7a-dirty by DeepMD contributors
 Loading plugin: kspace sog v3.1.3-83-g3cba3c7a-dirty by DeepMD contributors
-region box block 0 5 0 5 0 5
-create_box 1 box
-Created orthogonal box = (0 0 0) to (5 5 5)
+atom_style      charge
+region          box block 0 10 0 10 0 10
+create_box      1 box
+Created orthogonal box = (0 0 0) to (10 10 10)
   1 by 1 by 1 MPI processor grid
-create_atoms 1 single 1 1 1
+create_atoms    1 single 5.0 5.0 5.0
 Created 1 atoms
-  using lattice units in orthogonal box = (0 0 0) to (5 5 5)
+  using lattice units in orthogonal box = (0 0 0) to (10 10 10)
   create_atoms CPU = 0.000 seconds
-set atom 1 charge 0.1
+mass            1 1.0
+set             atom 1 charge 0.1
 Setting atom values ...
   1 settings made for charge
-mass 1 1.0
-pair_style coul/long 2.0
-pair_coeff * *
-kspace_style les 1e-6 n_dl 1.0 sigma 1.0 remove_self_interaction yes
-run 0
+
+pair_style      coul/long 6.0
+pair_coeff      * *
+
+kspace_style    sog 1.0e-6 amp 0.1 bandwidth 0.5
+thermo_style    custom step pe
+thermo          1
+run             0
 WARNING: No fixes with time integration, atoms won't move
 For more information see https://docs.lammps.org/err0028 (src/verlet.cpp:60)
 PPPM initialization ...
 WARNING: System is not charge neutral, net charge = 0.1
-For more information see https://docs.lammps.org/err0029 (src/kspace.cpp:330)
+For more information see https://docs.lammps.org/err0029 (src/kspace.cpp:335)
   using 12-bit tables for long-range coulomb
-  G vector (1/distance) = 0.70710678
-  grid = 8 8 8
+WARNING: For better accuracy use 'pair_modify table 0' (src/kspace.cpp:373)
+  G vector (1/distance) = 1
+  grid = 30 30 30
   stencil order = 5
-  estimated absolute RMS force accuracy = 0.0024650382
-  estimated relative force accuracy = 0.00017118743
+  estimated absolute RMS force accuracy = 3.4929588e-07
+  estimated relative force accuracy = 2.4257256e-08
   using double precision KISS FFT
-  3d grid and FFT values/proc = 4913 512
+  3d grid and FFT values/proc = 68921 27000
 Generated 0 of 0 mixed pair_coeff terms from geometric mixing rule
 Neighbor list info ...
   update: every = 1 steps, delay = 0 steps, check = yes
   max neighbors/atom: 2000, page size: 100000
-  master list distance cutoff = 4
-  ghost atom cutoff = 4
-  binsize = 2, bins = 3 3 3
+  master list distance cutoff = 8
+  ghost atom cutoff = 8
+  binsize = 4, bins = 3 3 3
   1 neighbor lists, perpetual/occasional/extra = 1 0 0
   (1) pair coul/long, perpetual
       attributes: half, newton on
       pair build: half/bin/atomonly/newton
       stencil: half/bin/3d
       bin: standard
-Per MPI rank memory allocation (min/avg/max) = 3.446 | 3.446 | 3.446 Mbytes
-   Step          Temp          E_pair         E_mol          TotEng         Press     
-         0   0             -0.040864926    0             -0.040864926   -127.91624    
-Loop time of 1.292e-06 on 1 procs for 0 steps with 1 atoms
+Per MPI rank memory allocation (min/avg/max) = 8.426 | 8.426 | 8.426 Mbytes
+   Step         PotEng    
+         0   0.001285803  
+Loop time of 1.302e-06 on 1 procs for 0 steps with 1 atoms
 
-154.8% CPU use with 1 MPI tasks x 1 OpenMP threads
+153.6% CPU use with 1 MPI tasks x no OpenMP threads
 
 MPI task timing breakdown:
 Section |  min time  |  avg time  |  max time  |%varavg| %total
@@ -70,7 +72,7 @@ Neigh   | 0          | 0          | 0          |   0.0 |  0.00
 Comm    | 0          | 0          | 0          |   0.0 |  0.00
 Output  | 0          | 0          | 0          |   0.0 |  0.00
 Modify  | 0          | 0          | 0          |   0.0 |  0.00
-Other   |            | 1.292e-06  |            |       |100.00
+Other   |            | 1.302e-06  |            |       |100.00
 
 Nlocal:              1 ave           1 max           1 min
 Histogram: 1 0 0 0 0 0 0 0 0 0

source/lmp/pair_deepmd.cpp

@@ -2,10 +2,14 @@
 #include <string.h>
 
 #include <cassert>
+#include <cctype>
+#include <cmath>
+#include <cstdlib>
 #include <iomanip>
 #include <iostream>
 #include <limits>
 #include <map>
+#include <regex>
 #include <sstream>
 
 #include "atom.h"
@@ -35,6 +39,21 @@
 using namespace LAMMPS_NS;
 using namespace std;
 
+namespace {
+
+bool has_json_suffix(const std::string& path) {
+  if (path.size() < 5) {
+    return false;
+  }
+  std::string suffix = path.substr(path.size() - 5);
+  for (char& c : suffix) {
+    c = static_cast<char>(std::tolower(static_cast<unsigned char>(c)));
+  }
+  return suffix == ".json";
+}
+
+}  // namespace
+
 static const char cite_user_deepmd_package[] =
     "USER-DEEPMD package:\n\n"
     "@article{Wang_ComputPhysCommun_2018_v228_p178,\n"
@@ -121,12 +140,90 @@ PairDeepMD::PairDeepMD(LAMMPS* lmp)
     : PairDeepBaseModel(
           lmp, cite_user_deepmd_package, deep_pot, deep_pot_model_devi) {
   latent_charge_to_q = false;
+  model_n_dl = 0.0;
+  has_model_n_dl = 0;
 }
 
 PairDeepMD::~PairDeepMD() {
   // Ensure base class destructor is called
 }
 
+bool PairDeepMD::parse_model_n_dl_from_json(const std::string& model_path,
+                                            double& n_dl_out) const {
+  if (!has_json_suffix(model_path)) {
+    return false;
+  }
+
+  std::string content;
+  try {
+    deepmd_compat::read_file_to_string(model_path, content);
+  } catch (...) {
+    return false;
+  }
+
+  static const std::regex n_dl_pattern(
+      R"("n_dl"\s*:\s*([-+]?(?:\d+\.?\d*|\.\d+)(?:[eE][-+]?\d+)?))");
+
+  std::smatch match;
+  if (!std::regex_search(content, match, n_dl_pattern)) {
+    return false;
+  }
+
+  const double candidate = std::strtod(match[1].str().c_str(), nullptr);
+  if (!(std::isfinite(candidate) && candidate > 0.0)) {
+    return false;
+  }
+
+  n_dl_out = candidate;
+  return true;
+}
+
+void PairDeepMD::detect_model_n_dl(const std::vector<std::string>& models) {
+  has_model_n_dl = 0;
+  model_n_dl = 0.0;
+
+  double first_n_dl = 0.0;
+  bool found = false;
+
+  for (const auto& model_path : models) {
+    double parsed = 0.0;
+    if (!parse_model_n_dl_from_json(model_path, parsed)) {
+      continue;
+    }
+    if (!found) {
+      first_n_dl = parsed;
+      found = true;
+      continue;
+    }
+    if (std::fabs(parsed - first_n_dl) > 1e-12 && comm->me == 0) {
+      error->warning(
+          FLERR,
+          "pair_style deepmd found inconsistent n_dl across model JSON files; "
+          "using the first model value");
+    }
+  }
+
+  if (found) {
+    has_model_n_dl = 1;
+    model_n_dl = first_n_dl;
+  }
+}
+
+void* PairDeepMD::extract(const char* str, int& dim) {
+  if (strcmp(str, "deepmd_model_n_dl") == 0) {
+    dim = 0;
+    if (has_model_n_dl) {
+      return (void*)&model_n_dl;
+    }
+    return nullptr;
+  }
+  if (strcmp(str, "deepmd_model_n_dl_available") == 0) {
+    dim = 0;
+    return (void*)&has_model_n_dl;
+  }
+  return PairDeepBaseModel::extract(str, dim);
+}
+
 void PairDeepMD::compute(int eflag, int vflag) {
   if (numb_models == 0) {
     return;
@@ -615,6 +712,7 @@ void PairDeepMD::settings(int narg, char** arg) {
   for (int ii = 0; ii < iarg; ++ii) {
     models.push_back(arg[ii]);
   }
+  detect_model_n_dl(models);
   numb_models = models.size();
   if (numb_models == 1) {
     try {

source/lmp/pair_deepmd.h

@@ -44,6 +44,7 @@ class PairDeepMD : public PairDeepBaseModel {
  public:
   PairDeepMD(class LAMMPS*);
   ~PairDeepMD() override;
+  void* extract(const char*, int&) override;
   void settings(int, char**) override;
   void coeff(int, char**) override;
   void compute(int, int) override;
@@ -55,8 +56,14 @@ class PairDeepMD : public PairDeepBaseModel {
   deepmd_compat::DeepPotModelDevi deep_pot_model_devi;
 
  private:
+  bool parse_model_n_dl_from_json(const std::string& model_path,
+                                  double& n_dl_out) const;
+  void detect_model_n_dl(const std::vector<std::string>& models);
+
   CommBrickDeepMD* commdata_;
   bool latent_charge_to_q;
+  double model_n_dl;
+  int has_model_n_dl;
 };
 
 }  // namespace LAMMPS_NS