更新SOGKSpace以支持CubeS₂的cubes2_phi_max参数，重构核参数计算，添加自能量修正功能

Author: Johnsonjz

deepmd/pt/model/model/sog_model.py

@@ -113,20 +113,31 @@ def _build_sog_lib_direct_kernel(
             dtype=real_dtype,
         )
 
-        kernel = sog_lib.Sog(
-            sog_arguments={
-                "use_atomwise": False,
-                "n_dl": float(fitting.n_dl),
-                "amp": amp_internal_runtime,
+        nlayers = getattr(self.atomic_model.descriptor, "nlayers", 1) if hasattr(self.atomic_model, "descriptor") else 1
+        # Build sog_arguments dict — prefer cubes2_phi_max, fall back to n_dl
+        sog_args: dict = {
+            "use_atomwise": False,
+            "amp": amp_internal_runtime,
                 "bandwidth": bw2_runtime,
                 "kernel_param_mode": "internal",
                 "kernel_tensor_mode": "external",
                 "remove_self_interaction": bool(fitting.remove_self_interaction),
                 "nufft": False,
                 "use_nufft": False,
+                "use_cubes2_fft": True,
+                "nlayers": nlayers,
                 "norm_factor": E2_PER_ANGSTROM_TO_EV,
                 "trainable_kernel": False,
-            },
+            }
+        # Prefer cubes2_phi_max (new API), fall back to n_dl (legacy)
+        if getattr(fitting, "cubes2_phi_max", None) is not None:
+            sog_args["cubes2_phi_max"] = float(fitting.cubes2_phi_max)
+        elif getattr(fitting, "n_dl", None) is not None:
+            sog_args["n_dl"] = float(fitting.n_dl)
+        # else: auto-default from SOG lib's Table III
+
+        kernel = sog_lib.Sog(
+            sog_arguments=sog_args,
             r_cut=float(self.get_rcut()),
         )
         return kernel

deepmd/pt/model/task/sog_energy_fitting.py

@@ -101,8 +101,11 @@ class SOGEnergyFittingNet(LRFittingNet):
         Base bandwidth used by SOG parameterization.
     M : int
         Number of geometric bandwidth levels.
-    n_dl : float
-        NUFFT long-range grid density control factor.
+    n_dl : float, optional (deprecated)
+        Legacy grid density control. Use `cubes2_phi_max` instead.
+    cubes2_phi_max : float, optional
+        φ = Δ/r_c grid control factor. Auto-defaults from Predescu 2020 Table III
+        when not specified. See sog lib documentation for recommended values.
     remove_self_interaction : bool
         If True, remove self interaction term in long-range correction.
     external_kspace : bool
@@ -141,7 +144,8 @@ def __init__(
         b: float | torch.Tensor | None = None,
         sigma: float | torch.Tensor | None = None,
         M: int | None = None,
-        n_dl: float | int = 1.0,
+        n_dl: float | int | None = None,
+        cubes2_phi_max: float | None = None,
         remove_self_interaction: bool = False,
         external_kspace: bool = False,
         **kwargs: Any,
@@ -174,27 +178,25 @@ def __init__(
             **kwargs,
         )
         if b is None:
-            b_tensor = torch.as_tensor(SOG_DEFAULT_B, dtype=dtype, device=device)
+            b_value = SOG_DEFAULT_B  # sog lib default (b=2)
         else:
             b_tensor = torch.as_tensor(b, dtype=dtype, device=device)
-        if b_tensor.numel() == 0:
-            b_tensor = torch.as_tensor(SOG_DEFAULT_B, dtype=dtype, device=device)
-        b_value = float(b_tensor.reshape(-1)[0].item())
+            b_value = float(b_tensor.reshape(-1)[0].item())
         if b_value <= 0.0:
             raise ValueError("`b` should be positive.")
 
         if sigma is None:
-            sigma_tensor = torch.as_tensor(SOG_DEFAULT_SIGMA, dtype=dtype, device=device)
+            sigma_value = SOG_DEFAULT_SIGMA  # will be overridden by sog lib via rcut
         else:
             sigma_tensor = torch.as_tensor(sigma, dtype=dtype, device=device)
-        if sigma_tensor.numel() == 0:
-            sigma_tensor = torch.as_tensor(SOG_DEFAULT_SIGMA, dtype=dtype, device=device)
-        sigma_value = float(sigma_tensor.reshape(-1)[0].item())
+            sigma_value = float(sigma_tensor.reshape(-1)[0].item())
         if sigma_value <= 0.0:
             raise ValueError("`sigma` should be positive.")
 
-        m_value = SOG_DEFAULT_M if M is None else int(M)
-        m_value = max(1, m_value)
+        if M is None:
+            m_value = SOG_DEFAULT_M  # sog lib default (M=12)
+        else:
+            m_value = max(1, int(M))
 
         if bandwidth is None:
             b_base = torch.tensor(b_value, dtype=dtype, device=device)
@@ -233,11 +235,18 @@ def __init__(
         # Store amp as sog-lib internal amplitude (already includes bw^2 factor).
         amp_tensor *= bandwidth_tensor
 
-        n_dl_value = float(n_dl)
-        if (not np.isfinite(n_dl_value)) or n_dl_value <= 0.0:
-            raise ValueError("`n_dl` should be a positive finite number.")
+        # Grid control: prefer cubes2_phi_max, fall back to n_dl (deprecated)
+        if n_dl is not None:
+            n_dl_value = float(n_dl)
+            if (not np.isfinite(n_dl_value)) or n_dl_value <= 0.0:
+                raise ValueError("`n_dl` should be a positive finite number.")
+        if cubes2_phi_max is not None:
+            phi_val = float(cubes2_phi_max)
+            if (not np.isfinite(phi_val)) or phi_val <= 0.0:
+                raise ValueError("`cubes2_phi_max` should be positive finite.")
 
-        self.n_dl = n_dl_value
+        self.n_dl = float(n_dl) if n_dl is not None else None
+        self.cubes2_phi_max = float(cubes2_phi_max) if cubes2_phi_max is not None else None
         self.amp = torch.nn.Parameter(
             amp_tensor,
             requires_grad=bool(self.trainable),
@@ -282,7 +291,10 @@ def serialize(self) -> dict:
         data["b"] = float(self.b)
         data["sigma"] = float(self.sigma)
         data["M"] = int(self.M)
-        data["n_dl"] = self.n_dl
+        if self.cubes2_phi_max is not None:
+            data["cubes2_phi_max"] = self.cubes2_phi_max
+        if self.n_dl is not None:
+            data["n_dl"] = self.n_dl  # legacy
         data["remove_self_interaction"] = bool(self.remove_self_interaction)
         data["external_kspace"] = bool(self.external_kspace)
         return data

deepmd/utils/argcheck.py

@@ -2276,11 +2276,13 @@ def fitting_sog_energy() -> list[Argument]:
             doc=doc_only_pt_supported + doc_M,
         ),
         Argument(
-            "n_dl",
+            "cubes2_phi_max",
             [float, int],
             optional=True,
-            default=1.0,
-            doc=doc_only_pt_supported + doc_n_dl,
+            default=None,
+            doc=doc_only_pt_supported
+            + "φ = Δ/r_c grid control for CubeS₂ FFT solver. "
+            + "Auto-defaults from Predescu 2020 Table III when not set.",
         ),
         Argument(
             "remove_self_interaction",

source/build_pt_gcc13/log.lammps

@@ -1,4 +1,96 @@
-LAMMPS (30 Mar 2026 - Development - patch_30Mar2026-1074-g7f680de296)
+LAMMPS (30 Mar 2026 - Development - 7f680de296)
 OMP_NUM_THREADS environment is not set. Defaulting to 1 thread.
   using 1 OpenMP thread(s) per MPI task
+# Minimal SOG virial test: orthogonal box, single step.
+
+units           metal
+atom_style      charge
+boundary        p p p
+
+plugin load     /data/home/public/jiangzhen/dp/deepmd-kit/build/sog_lmp_develop_plugin/lmp/plugin/libdeepmd_lmp.so
+Loading plugin: deepmd pair style  by Han Wang
+Loading plugin: deepspin pair style  by Duo Zhang
+Loading plugin: compute deeptensor/atom  by Han Wang
+Loading plugin: fix dplr  by Han Wang
+Loading plugin: kspace pppm/dplr  by Han Wang
+Loading plugin: kspace les  by DeepMD contributors
+Loading plugin: kspace sog  by DeepMD contributors
+
+lattice         fcc 4.0
+Lattice spacing in x,y,z = 4 4 4
+region          box block 0 3 0 3 0 3
+create_box      1 box
+Created orthogonal box = (0 0 0) to (12 12 12)
+  1 by 1 by 1 MPI processor grid
+create_atoms    1 random 256 12345 NULL
+Created 256 atoms
+  using lattice units in orthogonal box = (0 0 0) to (12 12 12)
+  create_atoms CPU = 0.000 seconds
+
+mass            1 1.0
+set             type 1 charge 0.5
+Setting atom values ...
+  256 settings made for charge
+
+pair_style      coul/long 10.0
+pair_coeff      1 1
+kspace_style    sog 1e-6 spline cubes2_4 b 2.0 m 12 n_dl 2.0 remove_self_interaction yes use_finufft no
+
+neighbor        2.0 bin
+neigh_modify    every 1 delay 0 check yes
+
+timestep        0.001
+thermo          1
+
+run             1
+WARNING: No fixes with time integration, atoms won't move
+For more information see https://docs.lammps.org/err0028 (src/verlet.cpp:60)
+WARNING: System is not charge neutral, net charge = 128
+For more information see https://docs.lammps.org/err0029 (src/kspace.cpp:328)
+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 = 12
+  ghost atom cutoff = 12
+  binsize = 6, bins = 2 2 2
+  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
+SOG virial check (analytic vs r_i-F_i): xx: ana=-0.0630387 rF=-0.35395 diff=0.291 (82.19%) yy: ana=-0.0544381 rF=0.0280625 diff=-0.0825 (293.99%) zz: ana=0.0213709 rF=0.065452 diff=-0.0441 (67.35%) xy: ana=0.000874211 rF=0.128377 diff=-0.128 (99.32%) xz: ana=-0.000350028 rF=0.0197886 diff=-0.0201 (101.77%) yz: ana=0.00274303 rF=-0.0172367 diff=0.02 (115.91%)
+Per MPI rank memory allocation (min/avg/max) = 4.053 | 4.053 | 4.053 Mbytes
+   Step          Temp          E_pair         E_mol          TotEng         Press     
+         0   0              42582.875      0              42582.875      13200403     
+SOG virial check (analytic vs r_i-F_i): xx: ana=-0.0630387 rF=-0.35395 diff=0.291 (82.19%) yy: ana=-0.0544381 rF=0.0280625 diff=-0.0825 (293.99%) zz: ana=0.0213709 rF=0.065452 diff=-0.0441 (67.35%) xy: ana=0.000874211 rF=0.128377 diff=-0.128 (99.32%) xz: ana=-0.000350028 rF=0.0197886 diff=-0.0201 (101.77%) yz: ana=0.00274303 rF=-0.0172367 diff=0.02 (115.91%)
+         1   0              42582.875      0              42582.875      13200403     
+Loop time of 0.00210485 on 1 procs for 1 steps with 256 atoms
+
+Performance: 41.048 ns/day, 0.585 hours/ns, 475.093 timesteps/s, 121.624 katom-step/s
+99.2% CPU use with 1 MPI tasks x 1 OpenMP threads
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 0.0018666  | 0.0018666  | 0.0018666  |   0.0 | 88.68
+Kspace  | 0.00020703 | 0.00020703 | 0.00020703 |   0.0 |  9.84
+Neigh   | 0          | 0          | 0          |   0.0 |  0.00
+Comm    | 2.0878e-05 | 2.0878e-05 | 2.0878e-05 |   0.0 |  0.99
+Output  | 4.647e-06  | 4.647e-06  | 4.647e-06  |   0.0 |  0.22
+Modify  | 5.88e-07   | 5.88e-07   | 5.88e-07   |   0.0 |  0.03
+Other   |            | 5.057e-06  |            |       |  0.24
+
+Nlocal:            256 ave         256 max         256 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Nghost:           6656 ave        6656 max        6656 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Neighs:         137723 ave      137723 max      137723 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+
+Total # of neighbors = 137723
+Ave neighs/atom = 537.98047
+Neighbor list builds = 0
+Dangerous builds = 0
 Total wall time: 0:00:00