diff --git a/deepmd/pt_expt/utils/serialization.py b/deepmd/pt_expt/utils/serialization.py
index 88bd839454..9c03783574 100644
--- a/deepmd/pt_expt/utils/serialization.py
+++ b/deepmd/pt_expt/utils/serialization.py
@@ -460,6 +460,9 @@ def _collect_metadata(model: torch.nn.Module, is_spin: bool = False) -> dict:
         "nnei": sum(model.get_sel()),
         "dim_fparam": model.get_dim_fparam(),
         "dim_aparam": model.get_dim_aparam(),
+        "dim_chg_spin": (
+            model.get_dim_chg_spin() if hasattr(model, "get_dim_chg_spin") else 0
+        ),
         "mixed_types": model.mixed_types(),
         "has_default_fparam": model.has_default_fparam(),
         "default_fparam": model.get_default_fparam(),
diff --git a/doc/third-party/lammps-command.md b/doc/third-party/lammps-command.md
index 70df75d22d..3cc5d955cc 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 _aparam_from_compute_ or _charge_spin_ or _atomic_ or _relative_ or _relative_v_ or _aparam_ or _ttm_
 
 <pre>
     <i>out_file</i> value = filename
@@ -62,6 +62,8 @@ pair_style deepmd models ... keyword value ...
         id = compute id used to update the frame parameter.
     <i>aparam_from_compute</i> value = id
         id = compute id used to update the atom parameter.
+    <i>charge_spin</i> value = parameters
+        parameters = the per-frame charge/spin values (dim_chg_spin numbers) required by models trained with a charge/spin embedding (e.g. DPA-3 with add_chg_spin_ebd). If omitted, the model's stored default_chg_spin is used.
     <i>atomic</i> = no value is required.
         If this keyword is set, the force model deviation of each atom will be output.
     <i>relative</i> value = level
@@ -88,6 +90,8 @@ compute    TEMP all temp
 
 pair_style deepmd ener.pb aparam_from_compute 1
 compute    1 all ke/atom
+
+pair_style deepmd dpa3.pth charge_spin 1.0 2.0
 ```
 
 ### Description
@@ -112,6 +116,7 @@ If the keyword `fparam` is set, the given frame parameter(s) will be fed to the
 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 `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 `charge_spin` is set, the given per-frame charge/spin value(s) will be fed to models that were trained with a charge/spin embedding (e.g. DPA-3 with `add_chg_spin_ebd`). If the keyword is not set, the model's stored `default_chg_spin` (if any) is used.
 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.
 
 Only a single `pair_coeff` command is used with the deepmd style which specifies atom names. These are mapped to LAMMPS atom types (integers from 1 to Ntypes) by specifying Ntypes additional arguments after `* *` in the `pair_coeff` command.
diff --git a/source/api_c/include/c_api.h b/source/api_c/include/c_api.h
index 358480b0ad..31b9e0626c 100644
--- a/source/api_c/include/c_api.h
+++ b/source/api_c/include/c_api.h
@@ -12,7 +12,7 @@ extern "C" {
 /** C API version. Bumped whenever the API is changed.
  * @since API version 22
  */
-#define DP_C_API_VERSION 26
+#define DP_C_API_VERSION 27
 
 /**
  * @brief Neighbor list.
@@ -675,6 +675,198 @@ extern void DP_DeepPotComputeNListf2(DP_DeepPot* dp,
                                      float* atomic_energy,
                                      float* atomic_virial);
 
+/**
+ * @brief Evaluate the energy, force and virial by using a DP. (double version,
+ *with charge_spin)
+ * @version 3
+ * @param[in] dp The DP to use.
+ * @param[in] nframes The number of frames.
+ * @param[in] natoms The number of atoms.
+ * @param[in] coord The coordinates of atoms. The array should be of size natoms
+ *x 3.
+ * @param[in] atype The atom types. The array should contain natoms ints.
+ * @param[in] cell The cell of the region. The array should be of size 9. Pass
+ *NULL if pbc is not used.
+ * @param[in] fparam The frame parameters. The array can be of size nframes x
+ *dim_fparam.
+ * @param[in] aparam The atom parameters. The array can be of size nframes x
+ *natoms x dim_aparam.
+ * @param[in] charge_spin The per-frame charge/spin input. The array can be of
+ *size nframes x dim_chg_spin. Pass NULL to use the model's stored
+ *default_chg_spin.
+ * @param[out] energy Output energy.
+ * @param[out] force Output force. The array should be of size natoms x 3.
+ * @param[out] virial Output virial. The array should be of size 9.
+ * @param[out] atomic_energy Output atomic energy. The array should be of size
+ *natoms.
+ * @param[out] atomic_virial Output atomic virial. The array should be of size
+ *natoms x 9.
+ * @warning The output arrays should be allocated before calling this function.
+ *Pass NULL if not required.
+ * @since API version 27
+ **/
+extern void DP_DeepPotCompute3(DP_DeepPot* dp,
+                               const int nframes,
+                               const int natoms,
+                               const double* coord,
+                               const int* atype,
+                               const double* cell,
+                               const double* fparam,
+                               const double* aparam,
+                               const double* charge_spin,
+                               double* energy,
+                               double* force,
+                               double* virial,
+                               double* atomic_energy,
+                               double* atomic_virial);
+
+/**
+ * @brief Evaluate the energy, force and virial by using a DP. (float version,
+ *with charge_spin)
+ * @version 3
+ * @param[in] dp The DP to use.
+ * @param[in] nframes The number of frames.
+ * @param[in] natoms The number of atoms.
+ * @param[in] coord The coordinates of atoms. The array should be of size natoms
+ *x 3.
+ * @param[in] atype The atom types. The array should contain natoms ints.
+ * @param[in] cell The cell of the region. The array should be of size 9. Pass
+ *NULL if pbc is not used.
+ * @param[in] fparam The frame parameters. The array can be of size nframes x
+ *dim_fparam.
+ * @param[in] aparam The atom parameters. The array can be of size nframes x
+ *natoms x dim_aparam.
+ * @param[in] charge_spin The per-frame charge/spin input. The array can be of
+ *size nframes x dim_chg_spin. Pass NULL to use the model's stored
+ *default_chg_spin.
+ * @param[out] energy Output energy.
+ * @param[out] force Output force. The array should be of size natoms x 3.
+ * @param[out] virial Output virial. The array should be of size 9.
+ * @param[out] atomic_energy Output atomic energy. The array should be of size
+ *natoms.
+ * @param[out] atomic_virial Output atomic virial. The array should be of size
+ *natoms x 9.
+ * @warning The output arrays should be allocated before calling this function.
+ *Pass NULL if not required.
+ * @since API version 27
+ **/
+extern void DP_DeepPotComputef3(DP_DeepPot* dp,
+                                const int nframes,
+                                const int natoms,
+                                const float* coord,
+                                const int* atype,
+                                const float* cell,
+                                const float* fparam,
+                                const float* aparam,
+                                const float* charge_spin,
+                                double* energy,
+                                float* force,
+                                float* virial,
+                                float* atomic_energy,
+                                float* atomic_virial);
+
+/**
+ * @brief Evaluate the energy, force and virial by using a DP with the neighbor
+ *list. (double version, with charge_spin)
+ * @version 3
+ * @param[in] dp The DP to use.
+ * @param[in] nframes The number of frames.
+ * @param[in] natoms The number of atoms.
+ * @param[in] coord The coordinates of atoms. The array should be of size natoms
+ *x 3.
+ * @param[in] atype The atom types. The array should contain natoms ints.
+ * @param[in] cell The cell of the region. The array should be of size 9. Pass
+ *NULL if pbc is not used.
+ * @param[in] nghost The number of ghost atoms.
+ * @param[in] nlist The neighbor list.
+ * @param[in] ago Update the internal neighbour list if ago is 0.
+ * @param[in] fparam The frame parameters. The array can be of size nframes x
+ *dim_fparam.
+ * @param[in] aparam The atom parameters. The array can be of size nframes x
+ *natoms x dim_aparam.
+ * @param[in] charge_spin The per-frame charge/spin input. The array can be of
+ *size nframes x dim_chg_spin. Pass NULL to use the model's stored
+ *default_chg_spin.
+ * @param[out] energy Output energy.
+ * @param[out] force Output force. The array should be of size natoms x 3.
+ * @param[out] virial Output virial. The array should be of size 9.
+ * @param[out] atomic_energy Output atomic energy. The array should be of size
+ *natoms.
+ * @param[out] atomic_virial Output atomic virial. The array should be of size
+ *natoms x 9.
+ * @warning The output arrays should be allocated before calling this function.
+ *Pass NULL if not required.
+ * @since API version 27
+ **/
+extern void DP_DeepPotComputeNList3(DP_DeepPot* dp,
+                                    const int nframes,
+                                    const int natoms,
+                                    const double* coord,
+                                    const int* atype,
+                                    const double* cell,
+                                    const int nghost,
+                                    const DP_Nlist* nlist,
+                                    const int ago,
+                                    const double* fparam,
+                                    const double* aparam,
+                                    const double* charge_spin,
+                                    double* energy,
+                                    double* force,
+                                    double* virial,
+                                    double* atomic_energy,
+                                    double* atomic_virial);
+
+/**
+ * @brief Evaluate the energy, force and virial by using a DP with the neighbor
+ *list. (float version, with charge_spin)
+ * @version 3
+ * @param[in] dp The DP to use.
+ * @param[in] nframes The number of frames.
+ * @param[in] natoms The number of atoms.
+ * @param[in] coord The coordinates of atoms. The array should be of size natoms
+ *x 3.
+ * @param[in] atype The atom types. The array should contain natoms ints.
+ * @param[in] cell The cell of the region. The array should be of size 9. Pass
+ *NULL if pbc is not used.
+ * @param[in] nghost The number of ghost atoms.
+ * @param[in] nlist The neighbor list.
+ * @param[in] ago Update the internal neighbour list if ago is 0.
+ * @param[in] fparam The frame parameters. The array can be of size nframes x
+ *dim_fparam.
+ * @param[in] aparam The atom parameters. The array can be of size nframes x
+ *natoms x dim_aparam.
+ * @param[in] charge_spin The per-frame charge/spin input. The array can be of
+ *size nframes x dim_chg_spin. Pass NULL to use the model's stored
+ *default_chg_spin.
+ * @param[out] energy Output energy.
+ * @param[out] force Output force. The array should be of size natoms x 3.
+ * @param[out] virial Output virial. The array should be of size 9.
+ * @param[out] atomic_energy Output atomic energy. The array should be of size
+ *natoms.
+ * @param[out] atomic_virial Output atomic virial. The array should be of size
+ *natoms x 9.
+ * @warning The output arrays should be allocated before calling this function.
+ *Pass NULL if not required.
+ * @since API version 27
+ **/
+extern void DP_DeepPotComputeNListf3(DP_DeepPot* dp,
+                                     const int nframes,
+                                     const int natoms,
+                                     const float* coord,
+                                     const int* atype,
+                                     const float* cell,
+                                     const int nghost,
+                                     const DP_Nlist* nlist,
+                                     const int ago,
+                                     const float* fparam,
+                                     const float* aparam,
+                                     const float* charge_spin,
+                                     double* energy,
+                                     float* force,
+                                     float* virial,
+                                     float* atomic_energy,
+                                     float* atomic_virial);
+
 /**
  * @brief Evaluate the energy, force, magnetic force and virial by using a DP
  *spin model with the neighbor list. (float version)
@@ -1346,6 +1538,206 @@ void DP_DeepPotModelDeviComputeNListf2(DP_DeepPotModelDevi* dp,
                                        float* atomic_energy,
                                        float* atomic_virial);
 
+/**
+ * @brief Evaluate energy, force and virial with a DP model deviation.
+ *(double version, with charge_spin)
+ * @version 3
+ * @param[in] dp The DP model deviation to use.
+ * @param[in] nframes The number of frames. Only 1 is supported.
+ * @param[in] natoms The number of atoms.
+ * @param[in] coord The coordinates of atoms. The array should be of size natoms
+ *x 3.
+ * @param[in] atype The atom types. The array should contain natoms ints.
+ * @param[in] cell The cell of the region. The array should be of size 9. Pass
+ *NULL if pbc is not used.
+ * @param[in] fparam The frame parameters. The array can be of size nframes x
+ *dim_fparam.
+ * @param[in] aparam The atom parameters. The array can be of size nframes x
+ *natoms x dim_aparam.
+ * @param[in] charge_spin The per-frame charge/spin input. The array should be
+ *of size nframes x dim_chg_spin. Pass NULL to use the model's stored
+ *default_chg_spin.
+ * @param[out] energy Output energy of all models.
+ * @param[out] force Output force of all models. The array should be of size
+ *nmodels x natoms x 3.
+ * @param[out] virial Output virial of all models. The array should be of size
+ *nmodels x 9.
+ * @param[out] atomic_energy Output atomic energy of all models. The array
+ *should be of size nmodels x natoms.
+ * @param[out] atomic_virial Output atomic virial of all models. The array
+ *should be of size nmodels x natoms x 9.
+ * @warning The output arrays should be allocated before calling this function.
+ *Pass NULL if not required.
+ * @since API version 27
+ **/
+void DP_DeepPotModelDeviCompute3(DP_DeepPotModelDevi* dp,
+                                 const int nframes,
+                                 const int natoms,
+                                 const double* coord,
+                                 const int* atype,
+                                 const double* cell,
+                                 const double* fparam,
+                                 const double* aparam,
+                                 const double* charge_spin,
+                                 double* energy,
+                                 double* force,
+                                 double* virial,
+                                 double* atomic_energy,
+                                 double* atomic_virial);
+
+/**
+ * @brief Evaluate energy, force and virial with a DP model deviation.
+ *(float version, with charge_spin)
+ * @version 3
+ * @param[in] dp The DP model deviation to use.
+ * @param[in] nframes The number of frames. Only 1 is supported.
+ * @param[in] natoms The number of atoms.
+ * @param[in] coord The coordinates of atoms. The array should be of size natoms
+ *x 3.
+ * @param[in] atype The atom types. The array should contain natoms ints.
+ * @param[in] cell The cell of the region. The array should be of size 9. Pass
+ *NULL if pbc is not used.
+ * @param[in] fparam The frame parameters. The array can be of size nframes x
+ *dim_fparam.
+ * @param[in] aparam The atom parameters. The array can be of size nframes x
+ *natoms x dim_aparam.
+ * @param[in] charge_spin The per-frame charge/spin input. The array should be
+ *of size nframes x dim_chg_spin. Pass NULL to use the model's stored
+ *default_chg_spin.
+ * @param[out] energy Output energy of all models.
+ * @param[out] force Output force of all models. The array should be of size
+ *nmodels x natoms x 3.
+ * @param[out] virial Output virial of all models. The array should be of size
+ *nmodels x 9.
+ * @param[out] atomic_energy Output atomic energy of all models. The array
+ *should be of size nmodels x natoms.
+ * @param[out] atomic_virial Output atomic virial of all models. The array
+ *should be of size nmodels x natoms x 9.
+ * @warning The output arrays should be allocated before calling this function.
+ *Pass NULL if not required.
+ * @since API version 27
+ **/
+void DP_DeepPotModelDeviComputef3(DP_DeepPotModelDevi* dp,
+                                  const int nframes,
+                                  const int natoms,
+                                  const float* coord,
+                                  const int* atype,
+                                  const float* cell,
+                                  const float* fparam,
+                                  const float* aparam,
+                                  const float* charge_spin,
+                                  double* energy,
+                                  float* force,
+                                  float* virial,
+                                  float* atomic_energy,
+                                  float* atomic_virial);
+
+/**
+ * @brief Evaluate energy, force and virial with a DP model deviation and a
+ *neighbor list. (double version, with charge_spin)
+ * @version 3
+ * @param[in] dp The DP model deviation to use.
+ * @param[in] nframes The number of frames. Only 1 is supported.
+ * @param[in] natoms The number of atoms.
+ * @param[in] coord The coordinates of atoms. The array should be of size natoms
+ *x 3.
+ * @param[in] atype The atom types. The array should contain natoms ints.
+ * @param[in] cell The cell of the region. The array should be of size 9. Pass
+ *NULL if pbc is not used.
+ * @param[in] nghost The number of ghost atoms.
+ * @param[in] nlist The neighbor list.
+ * @param[in] ago Update the internal neighbour list if ago is 0.
+ * @param[in] fparam The frame parameters. The array can be of size nframes x
+ *dim_fparam.
+ * @param[in] aparam The atom parameters. The array can be of size nframes x
+ *natoms x dim_aparam.
+ * @param[in] charge_spin The per-frame charge/spin input. The array should be
+ *of size nframes x dim_chg_spin. Pass NULL to use the model's stored
+ *default_chg_spin.
+ * @param[out] energy Output energy of all models.
+ * @param[out] force Output force of all models. The array should be of size
+ *nmodels x natoms x 3.
+ * @param[out] virial Output virial of all models. The array should be of size
+ *nmodels x 9.
+ * @param[out] atomic_energy Output atomic energy of all models. The array
+ *should be of size nmodels x natoms.
+ * @param[out] atomic_virial Output atomic virial of all models. The array
+ *should be of size nmodels x natoms x 9.
+ * @warning The output arrays should be allocated before calling this function.
+ *Pass NULL if not required.
+ * @since API version 27
+ **/
+void DP_DeepPotModelDeviComputeNList3(DP_DeepPotModelDevi* dp,
+                                      const int nframes,
+                                      const int natoms,
+                                      const double* coord,
+                                      const int* atype,
+                                      const double* cell,
+                                      const int nghost,
+                                      const DP_Nlist* nlist,
+                                      const int ago,
+                                      const double* fparam,
+                                      const double* aparam,
+                                      const double* charge_spin,
+                                      double* energy,
+                                      double* force,
+                                      double* virial,
+                                      double* atomic_energy,
+                                      double* atomic_virial);
+
+/**
+ * @brief Evaluate energy, force and virial with a DP model deviation and a
+ *neighbor list. (float version, with charge_spin)
+ * @version 3
+ * @param[in] dp The DP model deviation to use.
+ * @param[in] nframes The number of frames. Only 1 is supported.
+ * @param[in] natoms The number of atoms.
+ * @param[in] coord The coordinates of atoms. The array should be of size natoms
+ *x 3.
+ * @param[in] atype The atom types. The array should contain natoms ints.
+ * @param[in] cell The cell of the region. The array should be of size 9. Pass
+ *NULL if pbc is not used.
+ * @param[in] nghost The number of ghost atoms.
+ * @param[in] nlist The neighbor list.
+ * @param[in] ago Update the internal neighbour list if ago is 0.
+ * @param[in] fparam The frame parameters. The array can be of size nframes x
+ *dim_fparam.
+ * @param[in] aparam The atom parameters. The array can be of size nframes x
+ *natoms x dim_aparam.
+ * @param[in] charge_spin The per-frame charge/spin input. The array should be
+ *of size nframes x dim_chg_spin. Pass NULL to use the model's stored
+ *default_chg_spin.
+ * @param[out] energy Output energy of all models.
+ * @param[out] force Output force of all models. The array should be of size
+ *nmodels x natoms x 3.
+ * @param[out] virial Output virial of all models. The array should be of size
+ *nmodels x 9.
+ * @param[out] atomic_energy Output atomic energy of all models. The array
+ *should be of size nmodels x natoms.
+ * @param[out] atomic_virial Output atomic virial of all models. The array
+ *should be of size nmodels x natoms x 9.
+ * @warning The output arrays should be allocated before calling this function.
+ *Pass NULL if not required.
+ * @since API version 27
+ **/
+void DP_DeepPotModelDeviComputeNListf3(DP_DeepPotModelDevi* dp,
+                                       const int nframes,
+                                       const int natoms,
+                                       const float* coord,
+                                       const int* atype,
+                                       const float* cell,
+                                       const int nghost,
+                                       const DP_Nlist* nlist,
+                                       const int ago,
+                                       const float* fparam,
+                                       const float* aparam,
+                                       const float* charge_spin,
+                                       double* energy,
+                                       float* force,
+                                       float* virial,
+                                       float* atomic_energy,
+                                       float* atomic_virial);
+
 /**
  * @brief Evaluate the energy, force, magnetic force and virial by using a DP
  *spin model deviation with neighbor list. (float version)
@@ -1583,6 +1975,15 @@ int DP_DeepPotGetDimFParam(DP_DeepPot* dp);
  */
 int DP_DeepPotGetDimAParam(DP_DeepPot* dp);
 
+/**
+ * @brief Get the dimension of the charge/spin input of a DP.
+ * @param[in] dp The DP to use.
+ * @return The dimension of the charge/spin input (0 if the model has no
+ * charge/spin embedding).
+ * @since API version 27
+ */
+int DP_DeepPotGetDimChgSpin(DP_DeepPot* dp);
+
 /**
  * @brief Check whether the atomic dimension of atomic parameters is nall
  * instead of nloc.
@@ -1630,6 +2031,14 @@ int DP_DeepPotModelDeviGetDimFParam(DP_DeepPotModelDevi* dp);
  */
 int DP_DeepPotModelDeviGetDimAParam(DP_DeepPotModelDevi* dp);
 
+/**
+ * @brief Get the dimension of the charge/spin input of a DP Model Deviation.
+ * @param[in] dp The DP Model Deviation to use.
+ * @return The dimension of the charge/spin input (0 if none).
+ * @since API version 27
+ */
+int DP_DeepPotModelDeviGetDimChgSpin(DP_DeepPotModelDevi* dp);
+
 /**
  * @brief Check whether the atomic dimension of atomic parameters is nall
  * instead of nloc.
diff --git a/source/api_c/include/deepmd.hpp b/source/api_c/include/deepmd.hpp
index c3ca40b75f..1a68d0e227 100644
--- a/source/api_c/include/deepmd.hpp
+++ b/source/api_c/include/deepmd.hpp
@@ -53,6 +53,7 @@ inline void _DP_DeepPotCompute(DP_DeepPot* dp,
                                const FPTYPE* cell,
                                const FPTYPE* fparam,
                                const FPTYPE* aparam,
+                               const FPTYPE* charge_spin,
                                double* energy,
                                FPTYPE* force,
                                FPTYPE* virial,
@@ -68,13 +69,22 @@ inline void _DP_DeepPotCompute<double>(DP_DeepPot* dp,
                                        const double* cell,
                                        const double* fparam,
                                        const double* aparam,
+                                       const double* charge_spin,
                                        double* energy,
                                        double* force,
                                        double* virial,
                                        double* atomic_energy,
                                        double* atomic_virial) {
-  DP_DeepPotCompute2(dp, nframes, natom, coord, atype, cell, fparam, aparam,
-                     energy, force, virial, atomic_energy, atomic_virial);
+  // charge_spin == nullptr keeps the version-2 entry point so models without a
+  // charge/spin embedding still work against an older libdeepmd_c.
+  if (charge_spin) {
+    DP_DeepPotCompute3(dp, nframes, natom, coord, atype, cell, fparam, aparam,
+                       charge_spin, energy, force, virial, atomic_energy,
+                       atomic_virial);
+  } else {
+    DP_DeepPotCompute2(dp, nframes, natom, coord, atype, cell, fparam, aparam,
+                       energy, force, virial, atomic_energy, atomic_virial);
+  }
 }
 
 template <>
@@ -86,13 +96,20 @@ inline void _DP_DeepPotCompute<float>(DP_DeepPot* dp,
                                       const float* cell,
                                       const float* fparam,
                                       const float* aparam,
+                                      const float* charge_spin,
                                       double* energy,
                                       float* force,
                                       float* virial,
                                       float* atomic_energy,
                                       float* atomic_virial) {
-  DP_DeepPotComputef2(dp, nframes, natom, coord, atype, cell, fparam, aparam,
-                      energy, force, virial, atomic_energy, atomic_virial);
+  if (charge_spin) {
+    DP_DeepPotComputef3(dp, nframes, natom, coord, atype, cell, fparam, aparam,
+                        charge_spin, energy, force, virial, atomic_energy,
+                        atomic_virial);
+  } else {
+    DP_DeepPotComputef2(dp, nframes, natom, coord, atype, cell, fparam, aparam,
+                        energy, force, virial, atomic_energy, atomic_virial);
+  }
 }
 
 // support spin
@@ -167,6 +184,7 @@ inline void _DP_DeepPotComputeNList(DP_DeepPot* dp,
                                     const int ago,
                                     const FPTYPE* fparam,
                                     const FPTYPE* aparam,
+                                    const FPTYPE* charge_spin,
                                     double* energy,
                                     FPTYPE* force,
                                     FPTYPE* virial,
@@ -185,14 +203,21 @@ inline void _DP_DeepPotComputeNList<double>(DP_DeepPot* dp,
                                             const int ago,
                                             const double* fparam,
                                             const double* aparam,
+                                            const double* charge_spin,
                                             double* energy,
                                             double* force,
                                             double* virial,
                                             double* atomic_energy,
                                             double* atomic_virial) {
-  DP_DeepPotComputeNList2(dp, nframes, natom, coord, atype, cell, nghost, nlist,
-                          ago, fparam, aparam, energy, force, virial,
-                          atomic_energy, atomic_virial);
+  if (charge_spin) {
+    DP_DeepPotComputeNList3(dp, nframes, natom, coord, atype, cell, nghost,
+                            nlist, ago, fparam, aparam, charge_spin, energy,
+                            force, virial, atomic_energy, atomic_virial);
+  } else {
+    DP_DeepPotComputeNList2(dp, nframes, natom, coord, atype, cell, nghost,
+                            nlist, ago, fparam, aparam, energy, force, virial,
+                            atomic_energy, atomic_virial);
+  }
 }
 
 template <>
@@ -207,14 +232,21 @@ inline void _DP_DeepPotComputeNList<float>(DP_DeepPot* dp,
                                            const int ago,
                                            const float* fparam,
                                            const float* aparam,
+                                           const float* charge_spin,
                                            double* energy,
                                            float* force,
                                            float* virial,
                                            float* atomic_energy,
                                            float* atomic_virial) {
-  DP_DeepPotComputeNListf2(dp, nframes, natom, coord, atype, cell, nghost,
-                           nlist, ago, fparam, aparam, energy, force, virial,
-                           atomic_energy, atomic_virial);
+  if (charge_spin) {
+    DP_DeepPotComputeNListf3(dp, nframes, natom, coord, atype, cell, nghost,
+                             nlist, ago, fparam, aparam, charge_spin, energy,
+                             force, virial, atomic_energy, atomic_virial);
+  } else {
+    DP_DeepPotComputeNListf2(dp, nframes, natom, coord, atype, cell, nghost,
+                             nlist, ago, fparam, aparam, energy, force, virial,
+                             atomic_energy, atomic_virial);
+  }
 }
 
 // support spin
@@ -347,6 +379,7 @@ inline void _DP_DeepPotModelDeviCompute(DP_DeepPotModelDevi* dp,
                                         const FPTYPE* cell,
                                         const FPTYPE* fparam,
                                         const FPTYPE* aparam,
+                                        const FPTYPE* charge_spin,
                                         double* energy,
                                         FPTYPE* force,
                                         FPTYPE* virial,
@@ -361,14 +394,21 @@ inline void _DP_DeepPotModelDeviCompute<double>(DP_DeepPotModelDevi* dp,
                                                 const double* cell,
                                                 const double* fparam,
                                                 const double* aparam,
+                                                const double* charge_spin,
                                                 double* energy,
                                                 double* force,
                                                 double* virial,
                                                 double* atomic_energy,
                                                 double* atomic_virial) {
-  DP_DeepPotModelDeviCompute2(dp, 1, natom, coord, atype, cell, fparam, aparam,
-                              energy, force, virial, atomic_energy,
-                              atomic_virial);
+  if (charge_spin) {
+    DP_DeepPotModelDeviCompute3(dp, 1, natom, coord, atype, cell, fparam,
+                                aparam, charge_spin, energy, force, virial,
+                                atomic_energy, atomic_virial);
+  } else {
+    DP_DeepPotModelDeviCompute2(dp, 1, natom, coord, atype, cell, fparam,
+                                aparam, energy, force, virial, atomic_energy,
+                                atomic_virial);
+  }
 }
 
 template <>
@@ -379,14 +419,21 @@ inline void _DP_DeepPotModelDeviCompute<float>(DP_DeepPotModelDevi* dp,
                                                const float* cell,
                                                const float* fparam,
                                                const float* aparam,
+                                               const float* charge_spin,
                                                double* energy,
                                                float* force,
                                                float* virial,
                                                float* atomic_energy,
                                                float* atomic_virial) {
-  DP_DeepPotModelDeviComputef2(dp, 1, natom, coord, atype, cell, fparam, aparam,
-                               energy, force, virial, atomic_energy,
-                               atomic_virial);
+  if (charge_spin) {
+    DP_DeepPotModelDeviComputef3(dp, 1, natom, coord, atype, cell, fparam,
+                                 aparam, charge_spin, energy, force, virial,
+                                 atomic_energy, atomic_virial);
+  } else {
+    DP_DeepPotModelDeviComputef2(dp, 1, natom, coord, atype, cell, fparam,
+                                 aparam, energy, force, virial, atomic_energy,
+                                 atomic_virial);
+  }
 }
 
 template <typename FPTYPE>
@@ -456,6 +503,7 @@ inline void _DP_DeepPotModelDeviComputeNList(DP_DeepPotModelDevi* dp,
                                              const int ago,
                                              const FPTYPE* fparam,
                                              const FPTYPE* aparam,
+                                             const FPTYPE* charge_spin,
                                              double* energy,
                                              FPTYPE* force,
                                              FPTYPE* virial,
@@ -473,14 +521,21 @@ inline void _DP_DeepPotModelDeviComputeNList<double>(DP_DeepPotModelDevi* dp,
                                                      const int ago,
                                                      const double* fparam,
                                                      const double* aparam,
+                                                     const double* charge_spin,
                                                      double* energy,
                                                      double* force,
                                                      double* virial,
                                                      double* atomic_energy,
                                                      double* atomic_virial) {
-  DP_DeepPotModelDeviComputeNList2(dp, 1, natom, coord, atype, cell, nghost,
-                                   nlist, ago, fparam, aparam, energy, force,
-                                   virial, atomic_energy, atomic_virial);
+  if (charge_spin) {
+    DP_DeepPotModelDeviComputeNList3(
+        dp, 1, natom, coord, atype, cell, nghost, nlist, ago, fparam, aparam,
+        charge_spin, energy, force, virial, atomic_energy, atomic_virial);
+  } else {
+    DP_DeepPotModelDeviComputeNList2(dp, 1, natom, coord, atype, cell, nghost,
+                                     nlist, ago, fparam, aparam, energy, force,
+                                     virial, atomic_energy, atomic_virial);
+  }
 }
 
 template <>
@@ -494,14 +549,21 @@ inline void _DP_DeepPotModelDeviComputeNList<float>(DP_DeepPotModelDevi* dp,
                                                     const int ago,
                                                     const float* fparam,
                                                     const float* aparam,
+                                                    const float* charge_spin,
                                                     double* energy,
                                                     float* force,
                                                     float* virial,
                                                     float* atomic_energy,
                                                     float* atomic_virial) {
-  DP_DeepPotModelDeviComputeNListf2(dp, 1, natom, coord, atype, cell, nghost,
-                                    nlist, ago, fparam, aparam, energy, force,
-                                    virial, atomic_energy, atomic_virial);
+  if (charge_spin) {
+    DP_DeepPotModelDeviComputeNListf3(
+        dp, 1, natom, coord, atype, cell, nghost, nlist, ago, fparam, aparam,
+        charge_spin, energy, force, virial, atomic_energy, atomic_virial);
+  } else {
+    DP_DeepPotModelDeviComputeNListf2(dp, 1, natom, coord, atype, cell, nghost,
+                                      nlist, ago, fparam, aparam, energy, force,
+                                      virial, atomic_energy, atomic_virial);
+  }
 }
 
 template <typename FPTYPE>
@@ -800,6 +862,34 @@ inline double* _DP_Get_Energy_Pointer(double& vec, const int nframes) {
 
 namespace deepmd {
 namespace hpp {
+/**
+ * @brief Validate a runtime charge_spin vector and return a raw pointer to it
+ * (or nullptr when empty, which falls back to the model's default_chg_spin).
+ * @param[in] charge_spin The per-frame charge/spin values.
+ * @param[in] dchgspin The model's charge/spin dimension (0 if unsupported).
+ * @param[in] nframes The number of frames.
+ * @return charge_spin.data() if non-empty and valid, otherwise nullptr.
+ */
+template <typename FPTYPE>
+inline const FPTYPE* validate_charge_spin(
+    const std::vector<FPTYPE>& charge_spin,
+    const int dchgspin,
+    const unsigned int nframes) {
+  if (charge_spin.empty()) {
+    return nullptr;
+  }
+  if (dchgspin == 0) {
+    throw deepmd::hpp::deepmd_exception(
+        "charge_spin was provided, but this model does not support "
+        "charge/spin conditioning");
+  }
+  if (charge_spin.size() != static_cast<size_t>(nframes) * dchgspin) {
+    throw deepmd::hpp::deepmd_exception(
+        "the dim of charge_spin provided is not consistent with what the "
+        "model uses");
+  }
+  return charge_spin.data();
+}
 /**
  * @brief Neighbor list.
  **/
@@ -1027,7 +1117,7 @@ class DeepPot : public DeepBaseModel {
   /**
    * @brief DP constructor without initialization.
    **/
-  DeepPot() : dp(nullptr) {};
+  DeepPot() : dp(nullptr), dchgspin(0) {};
   ~DeepPot() { DP_DeleteDeepPot(dp); };
   /**
    * @brief DP constructor with initialization.
@@ -1038,7 +1128,7 @@ class DeepPot : public DeepBaseModel {
   DeepPot(const std::string& model,
           const int& gpu_rank = 0,
           const std::string& file_content = "")
-      : dp(nullptr) {
+      : dp(nullptr), dchgspin(0) {
     try {
       init(model, gpu_rank, file_content);
     } catch (...) {
@@ -1069,11 +1159,22 @@ class DeepPot : public DeepBaseModel {
     DP_CHECK_OK(DP_DeepPotCheckOK, dp);
     dfparam = DP_DeepPotGetDimFParam(dp);
     daparam = DP_DeepPotGetDimAParam(dp);
+    dchgspin = DP_DeepPotGetDimChgSpin(dp);
     aparam_nall = DP_DeepPotIsAParamNAll(dp);
     has_default_fparam_ = DP_DeepPotHasDefaultFParam(dp);
     dpbase = (DP_DeepBaseModel*)dp;
   };
 
+  /**
+   * @brief Get the dimension of the charge/spin embedding input.
+   * @return The dimension of the charge/spin input (0 if the model has no
+   *charge/spin embedding).
+   **/
+  int dim_chg_spin() const {
+    assert(dp);
+    return dchgspin;
+  }
+
   /**
    * @brief Evaluate the energy, force and virial by using this DP.
    * @param[out] ener The system energy.
@@ -1092,6 +1193,10 @@ class DeepPot : public DeepBaseModel {
    * nframes x natoms x dim_aparam.
    * natoms x dim_aparam. Then all frames are assumed to be provided with the
    *same aparam.
+   * @param[in] charge_spin The per-frame charge/spin input. The array can be
+   *of size nframes x dim_chg_spin. Then all frames are assumed to be provided
+   *with the same charge_spin. Leave it empty to use the model's stored
+   *default_chg_spin.
    * @warning Natoms should not be zero when computing multiple frames.
    **/
   template <typename VALUETYPE, typename ENERGYVTYPE>
@@ -1103,7 +1208,8 @@ class DeepPot : public DeepBaseModel {
       const std::vector<int>& atype,
       const std::vector<VALUETYPE>& box,
       const std::vector<VALUETYPE>& fparam = std::vector<VALUETYPE>(),
-      const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>()) {
+      const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>(),
+      const std::vector<VALUETYPE>& charge_spin = std::vector<VALUETYPE>()) {
     unsigned int natoms = atype.size();
     unsigned int nframes = natoms > 0 ? coord.size() / natoms / 3 : 1;
     assert(nframes * natoms * 3 == coord.size());
@@ -1124,10 +1230,14 @@ class DeepPot : public DeepBaseModel {
     tile_fparam_aparam(aparam_, nframes, natoms * daparam, aparam);
     const VALUETYPE* fparam__ = !fparam_.empty() ? &fparam_[0] : nullptr;
     const VALUETYPE* aparam__ = !aparam_.empty() ? &aparam_[0] : nullptr;
+    // charge_spin routes to the version-3 C API; nullptr keeps version-2 so
+    // non-charge_spin models still work against an older libdeepmd_c.
+    const VALUETYPE* charge_spin__ =
+        validate_charge_spin(charge_spin, dchgspin, nframes);
 
     _DP_DeepPotCompute<VALUETYPE>(dp, nframes, natoms, coord_, atype_, box_,
-                                  fparam__, aparam__, ener_, force_, virial_,
-                                  nullptr, nullptr);
+                                  fparam__, aparam__, charge_spin__, ener_,
+                                  force_, virial_, nullptr, nullptr);
     DP_CHECK_OK(DP_DeepPotCheckOK, dp);
   };
   /**
@@ -1151,6 +1261,10 @@ class DeepPot : public DeepBaseModel {
    * nframes x natoms x dim_aparam.
    * natoms x dim_aparam. Then all frames are assumed to be provided with the
    *same aparam.
+   * @param[in] charge_spin The per-frame charge/spin input. The array can be
+   *of size nframes x dim_chg_spin. Then all frames are assumed to be provided
+   *with the same charge_spin. Leave it empty to use the model's stored
+   *default_chg_spin.
    * @warning Natoms should not be zero when computing multiple frames.
    **/
   template <typename VALUETYPE, typename ENERGYVTYPE>
@@ -1164,7 +1278,8 @@ class DeepPot : public DeepBaseModel {
       const std::vector<int>& atype,
       const std::vector<VALUETYPE>& box,
       const std::vector<VALUETYPE>& fparam = std::vector<VALUETYPE>(),
-      const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>()) {
+      const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>(),
+      const std::vector<VALUETYPE>& charge_spin = std::vector<VALUETYPE>()) {
     unsigned int natoms = atype.size();
     unsigned int nframes = natoms > 0 ? coord.size() / natoms / 3 : 1;
     assert(nframes * natoms * 3 == coord.size());
@@ -1190,10 +1305,12 @@ class DeepPot : public DeepBaseModel {
     tile_fparam_aparam(aparam_, nframes, natoms * daparam, aparam);
     const VALUETYPE* fparam__ = !fparam_.empty() ? &fparam_[0] : nullptr;
     const VALUETYPE* aparam__ = !aparam_.empty() ? &aparam_[0] : nullptr;
+    const VALUETYPE* charge_spin__ =
+        validate_charge_spin(charge_spin, dchgspin, nframes);
 
-    _DP_DeepPotCompute<VALUETYPE>(dp, nframes, natoms, coord_, atype_, box_,
-                                  fparam__, aparam__, ener_, force_, virial_,
-                                  atomic_ener_, atomic_virial_);
+    _DP_DeepPotCompute<VALUETYPE>(
+        dp, nframes, natoms, coord_, atype_, box_, fparam__, aparam__,
+        charge_spin__, ener_, force_, virial_, atomic_ener_, atomic_virial_);
     DP_CHECK_OK(DP_DeepPotCheckOK, dp);
   };
 
@@ -1219,6 +1336,10 @@ class DeepPot : public DeepBaseModel {
    * nframes x natoms x dim_aparam.
    * natoms x dim_aparam. Then all frames are assumed to be provided with the
    *same aparam.
+   * @param[in] charge_spin The per-frame charge/spin input. The array can be
+   *of size nframes x dim_chg_spin. Then all frames are assumed to be provided
+   *with the same charge_spin. Leave it empty to use the model's stored
+   *default_chg_spin.
    * @warning Natoms should not be zero when computing multiple frames.
    **/
   template <typename VALUETYPE, typename ENERGYVTYPE>
@@ -1233,7 +1354,8 @@ class DeepPot : public DeepBaseModel {
       const InputNlist& lmp_list,
       const int& ago,
       const std::vector<VALUETYPE>& fparam = std::vector<VALUETYPE>(),
-      const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>()) {
+      const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>(),
+      const std::vector<VALUETYPE>& charge_spin = std::vector<VALUETYPE>()) {
     unsigned int natoms = atype.size();
     unsigned int nframes = natoms > 0 ? coord.size() / natoms / 3 : 1;
     assert(nframes * natoms * 3 == coord.size());
@@ -1257,10 +1379,13 @@ class DeepPot : public DeepBaseModel {
                        aparam);
     const VALUETYPE* fparam__ = !fparam_.empty() ? &fparam_[0] : nullptr;
     const VALUETYPE* aparam__ = !aparam_.empty() ? &aparam_[0] : nullptr;
+    const VALUETYPE* charge_spin__ =
+        validate_charge_spin(charge_spin, dchgspin, nframes);
 
-    _DP_DeepPotComputeNList<VALUETYPE>(
-        dp, nframes, natoms, coord_, atype_, box_, nghost, lmp_list.nl, ago,
-        fparam__, aparam__, ener_, force_, virial_, nullptr, nullptr);
+    _DP_DeepPotComputeNList<VALUETYPE>(dp, nframes, natoms, coord_, atype_,
+                                       box_, nghost, lmp_list.nl, ago, fparam__,
+                                       aparam__, charge_spin__, ener_, force_,
+                                       virial_, nullptr, nullptr);
     DP_CHECK_OK(DP_DeepPotCheckOK, dp);
   };
   /**
@@ -1287,6 +1412,10 @@ class DeepPot : public DeepBaseModel {
    * nframes x natoms x dim_aparam.
    * natoms x dim_aparam. Then all frames are assumed to be provided with the
    *same aparam.
+   * @param[in] charge_spin The per-frame charge/spin input. The array can be
+   *of size nframes x dim_chg_spin. Then all frames are assumed to be provided
+   *with the same charge_spin. Leave it empty to use the model's stored
+   *default_chg_spin.
    * @warning Natoms should not be zero when computing multiple frames.
    **/
   template <typename VALUETYPE, typename ENERGYVTYPE>
@@ -1303,7 +1432,8 @@ class DeepPot : public DeepBaseModel {
       const InputNlist& lmp_list,
       const int& ago,
       const std::vector<VALUETYPE>& fparam = std::vector<VALUETYPE>(),
-      const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>()) {
+      const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>(),
+      const std::vector<VALUETYPE>& charge_spin = std::vector<VALUETYPE>()) {
     unsigned int natoms = atype.size();
     unsigned int nframes = natoms > 0 ? coord.size() / natoms / 3 : 1;
     assert(nframes * natoms * 3 == coord.size());
@@ -1332,11 +1462,13 @@ class DeepPot : public DeepBaseModel {
                        aparam);
     const VALUETYPE* fparam__ = !fparam_.empty() ? &fparam_[0] : nullptr;
     const VALUETYPE* aparam__ = !aparam_.empty() ? &aparam_[0] : nullptr;
+    const VALUETYPE* charge_spin__ =
+        validate_charge_spin(charge_spin, dchgspin, nframes);
 
     _DP_DeepPotComputeNList<VALUETYPE>(dp, nframes, natoms, coord_, atype_,
                                        box_, nghost, lmp_list.nl, ago, fparam__,
-                                       aparam__, ener_, force_, virial_,
-                                       atomic_ener_, atomic_virial_);
+                                       aparam__, charge_spin__, ener_, force_,
+                                       virial_, atomic_ener_, atomic_virial_);
     DP_CHECK_OK(DP_DeepPotCheckOK, dp);
   };
   /**
@@ -1465,6 +1597,7 @@ class DeepPot : public DeepBaseModel {
 
  private:
   DP_DeepPot* dp;
+  int dchgspin;
 };
 
 class DeepSpin : public DeepBaseModel {
@@ -1881,6 +2014,11 @@ class DeepBaseModelDevi {
     assert(dpbase);
     return has_default_fparam_;
   }
+  /**
+   * @brief Get the dimension of the charge/spin input.
+   * @return Always 0; charge_spin is not supported via the C-API path.
+   **/
+  int dim_chg_spin() const { return 0; }
   /**
    * @param[out] avg The average of vectors.
    * @param[in] xx The vectors of all models.
@@ -2048,13 +2186,14 @@ class DeepPotModelDevi : public DeepBaseModelDevi {
   /**
    * @brief DP model deviation constructor without initialization.
    **/
-  DeepPotModelDevi() : dp(nullptr) {};
+  DeepPotModelDevi() : dp(nullptr), dchgspin(0) {};
   ~DeepPotModelDevi() { DP_DeleteDeepPotModelDevi(dp); };
   /**
    * @brief DP model deviation constructor with initialization.
    * @param[in] models The names of the frozen model file.
    **/
-  DeepPotModelDevi(const std::vector<std::string>& models) : dp(nullptr) {
+  DeepPotModelDevi(const std::vector<std::string>& models)
+      : dp(nullptr), dchgspin(0) {
     try {
       init(models);
     } catch (...) {
@@ -2103,11 +2242,22 @@ class DeepPotModelDevi : public DeepBaseModelDevi {
     numb_models = models.size();
     dfparam = DP_DeepPotModelDeviGetDimFParam(dp);
     daparam = DP_DeepPotModelDeviGetDimAParam(dp);
+    dchgspin = DP_DeepPotModelDeviGetDimChgSpin(dp);
     aparam_nall = DP_DeepPotModelDeviIsAParamNAll(dp);
     has_default_fparam_ = DP_DeepPotModelDeviHasDefaultFParam(dp);
     dpbase = (DP_DeepBaseModelDevi*)dp;
   };
 
+  /**
+   * @brief Get the dimension of the charge/spin input.
+   * @return The dimension of the charge/spin input (0 if the models have no
+   *charge/spin embedding).
+   **/
+  int dim_chg_spin() const {
+    assert(dp);
+    return dchgspin;
+  }
+
   /**
    * @brief Evaluate the energy, force and virial by using this DP model
    *deviation.
@@ -2127,6 +2277,10 @@ class DeepPotModelDevi : public DeepBaseModelDevi {
    * nframes x natoms x dim_aparam.
    * natoms x dim_aparam. Then all frames are assumed to be provided with the
    *same aparam.
+   * @param[in] charge_spin The per-frame charge/spin input. The array can be
+   *of size nframes x dim_chg_spin. Then all frames are assumed to be provided
+   *with the same charge_spin. Leave it empty to use the model's stored
+   *default_chg_spin.
    **/
   template <typename VALUETYPE>
   void compute(
@@ -2137,7 +2291,9 @@ class DeepPotModelDevi : public DeepBaseModelDevi {
       const std::vector<int>& atype,
       const std::vector<VALUETYPE>& box,
       const std::vector<VALUETYPE>& fparam = std::vector<VALUETYPE>(),
-      const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>()) {
+      const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>(),
+      const std::vector<VALUETYPE>& charge_spin = std::vector<VALUETYPE>()) {
+    // charge_spin is not supported via the C-API model-deviation path.
     unsigned int natoms = atype.size();
     unsigned int nframes = 1;
     assert(natoms * 3 == coord.size());
@@ -2163,10 +2319,12 @@ class DeepPotModelDevi : public DeepBaseModelDevi {
     tile_fparam_aparam(aparam_, nframes, natoms * daparam, aparam);
     const VALUETYPE* fparam__ = !fparam_.empty() ? &fparam_[0] : nullptr;
     const VALUETYPE* aparam__ = !aparam_.empty() ? &aparam_[0] : nullptr;
+    const VALUETYPE* charge_spin__ =
+        validate_charge_spin(charge_spin, dchgspin, nframes);
 
-    _DP_DeepPotModelDeviCompute<VALUETYPE>(dp, natoms, coord_, atype_, box_,
-                                           fparam__, aparam__, ener_, force_,
-                                           virial_, nullptr, nullptr);
+    _DP_DeepPotModelDeviCompute<VALUETYPE>(
+        dp, natoms, coord_, atype_, box_, fparam__, aparam__, charge_spin__,
+        ener_, force_, virial_, nullptr, nullptr);
     DP_CHECK_OK(DP_DeepPotModelDeviCheckOK, dp);
 
     // reshape
@@ -2206,6 +2364,10 @@ class DeepPotModelDevi : public DeepBaseModelDevi {
    * nframes x natoms x dim_aparam.
    * natoms x dim_aparam. Then all frames are assumed to be provided with the
    *same aparam.
+   * @param[in] charge_spin The per-frame charge/spin input. The array can be
+   *of size nframes x dim_chg_spin. Then all frames are assumed to be provided
+   *with the same charge_spin. Leave it empty to use the model's stored
+   *default_chg_spin.
    **/
   template <typename VALUETYPE>
   void compute(
@@ -2218,7 +2380,9 @@ class DeepPotModelDevi : public DeepBaseModelDevi {
       const std::vector<int>& atype,
       const std::vector<VALUETYPE>& box,
       const std::vector<VALUETYPE>& fparam = std::vector<VALUETYPE>(),
-      const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>()) {
+      const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>(),
+      const std::vector<VALUETYPE>& charge_spin = std::vector<VALUETYPE>()) {
+    // charge_spin is not supported via the C-API model-deviation path.
     unsigned int natoms = atype.size();
     unsigned int nframes = 1;
     assert(natoms * 3 == coord.size());
@@ -2248,10 +2412,12 @@ class DeepPotModelDevi : public DeepBaseModelDevi {
     tile_fparam_aparam(aparam_, nframes, natoms * daparam, aparam);
     const VALUETYPE* fparam__ = !fparam_.empty() ? &fparam_[0] : nullptr;
     const VALUETYPE* aparam__ = !aparam_.empty() ? &aparam_[0] : nullptr;
+    const VALUETYPE* charge_spin__ =
+        validate_charge_spin(charge_spin, dchgspin, nframes);
 
     _DP_DeepPotModelDeviCompute<VALUETYPE>(
-        dp, natoms, coord_, atype_, box_, fparam__, aparam__, ener_, force_,
-        virial_, atomic_ener_, atomic_virial_);
+        dp, natoms, coord_, atype_, box_, fparam__, aparam__, charge_spin__,
+        ener_, force_, virial_, atomic_ener_, atomic_virial_);
     DP_CHECK_OK(DP_DeepPotModelDeviCheckOK, dp);
 
     // reshape
@@ -2303,6 +2469,10 @@ class DeepPotModelDevi : public DeepBaseModelDevi {
    * nframes x natoms x dim_aparam.
    * natoms x dim_aparam. Then all frames are assumed to be provided with the
    *same aparam.
+   * @param[in] charge_spin The per-frame charge/spin input. The array can be
+   *of size nframes x dim_chg_spin. Then all frames are assumed to be provided
+   *with the same charge_spin. Leave it empty to use the model's stored
+   *default_chg_spin.
    **/
   template <typename VALUETYPE>
   void compute(
@@ -2316,7 +2486,8 @@ class DeepPotModelDevi : public DeepBaseModelDevi {
       const InputNlist& lmp_list,
       const int& ago,
       const std::vector<VALUETYPE>& fparam = std::vector<VALUETYPE>(),
-      const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>()) {
+      const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>(),
+      const std::vector<VALUETYPE>& charge_spin = std::vector<VALUETYPE>()) {
     unsigned int natoms = atype.size();
     unsigned int nframes = 1;
     assert(natoms * 3 == coord.size());
@@ -2345,10 +2516,12 @@ class DeepPotModelDevi : public DeepBaseModelDevi {
                        aparam);
     const VALUETYPE* fparam__ = !fparam_.empty() ? &fparam_[0] : nullptr;
     const VALUETYPE* aparam__ = !aparam_.empty() ? &aparam_[0] : nullptr;
+    const VALUETYPE* charge_spin__ =
+        validate_charge_spin(charge_spin, dchgspin, nframes);
 
     _DP_DeepPotModelDeviComputeNList<VALUETYPE>(
         dp, natoms, coord_, atype_, box_, nghost, lmp_list.nl, ago, fparam__,
-        aparam__, ener_, force_, virial_, nullptr, nullptr);
+        aparam__, charge_spin__, ener_, force_, virial_, nullptr, nullptr);
     DP_CHECK_OK(DP_DeepPotModelDeviCheckOK, dp);
 
     // reshape
@@ -2391,6 +2564,10 @@ class DeepPotModelDevi : public DeepBaseModelDevi {
    * nframes x natoms x dim_aparam.
    * natoms x dim_aparam. Then all frames are assumed to be provided with the
    *same aparam.
+   * @param[in] charge_spin The per-frame charge/spin input. The array can be
+   *of size nframes x dim_chg_spin. Then all frames are assumed to be provided
+   *with the same charge_spin. Leave it empty to use the model's stored
+   *default_chg_spin.
    **/
   template <typename VALUETYPE>
   void compute(
@@ -2406,7 +2583,8 @@ class DeepPotModelDevi : public DeepBaseModelDevi {
       const InputNlist& lmp_list,
       const int& ago,
       const std::vector<VALUETYPE>& fparam = std::vector<VALUETYPE>(),
-      const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>()) {
+      const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>(),
+      const std::vector<VALUETYPE>& charge_spin = std::vector<VALUETYPE>()) {
     unsigned int natoms = atype.size();
     unsigned int nframes = 1;
     assert(natoms * 3 == coord.size());
@@ -2439,10 +2617,13 @@ class DeepPotModelDevi : public DeepBaseModelDevi {
                        aparam);
     const VALUETYPE* fparam__ = !fparam_.empty() ? &fparam_[0] : nullptr;
     const VALUETYPE* aparam__ = !aparam_.empty() ? &aparam_[0] : nullptr;
+    const VALUETYPE* charge_spin__ =
+        validate_charge_spin(charge_spin, dchgspin, nframes);
 
     _DP_DeepPotModelDeviComputeNList<VALUETYPE>(
         dp, natoms, coord_, atype_, box_, nghost, lmp_list.nl, ago, fparam__,
-        aparam__, ener_, force_, virial_, atomic_ener_, atomic_virial_);
+        aparam__, charge_spin__, ener_, force_, virial_, atomic_ener_,
+        atomic_virial_);
     DP_CHECK_OK(DP_DeepPotModelDeviCheckOK, dp);
 
     // reshape
@@ -2474,6 +2655,7 @@ class DeepPotModelDevi : public DeepBaseModelDevi {
 
  private:
   DP_DeepPotModelDevi* dp;
+  int dchgspin;
 };
 
 class DeepSpinModelDevi : public DeepBaseModelDevi {
diff --git a/source/api_c/src/c_api.cc b/source/api_c/src/c_api.cc
index b0e789648e..1346bfbf20 100644
--- a/source/api_c/src/c_api.cc
+++ b/source/api_c/src/c_api.cc
@@ -251,7 +251,8 @@ inline void DP_DeepPotCompute_variant(DP_DeepPot* dp,
                                       VALUETYPE* force,
                                       VALUETYPE* virial,
                                       VALUETYPE* atomic_energy,
-                                      VALUETYPE* atomic_virial) {
+                                      VALUETYPE* atomic_virial,
+                                      const VALUETYPE* charge_spin = nullptr) {
   // init C++ vectors from C arrays
   std::vector<VALUETYPE> coord_(coord, coord + nframes * natoms * 3);
   std::vector<int> atype_(atype, atype + natoms);
@@ -268,15 +269,22 @@ inline void DP_DeepPotCompute_variant(DP_DeepPot* dp,
   if (aparam) {
     aparam_.assign(aparam, aparam + nframes * natoms * dp->daparam);
   }
+  // charge_spin is converted to double for the api_cc::DeepPot interface
+  // (the compute backend stores/uses charge_spin as float64).
+  std::vector<double> charge_spin_;
+  if (charge_spin) {
+    charge_spin_.assign(charge_spin,
+                        charge_spin + nframes * dp->dp.dim_chg_spin());
+  }
   std::vector<double> e;
   std::vector<VALUETYPE> f, v, ae, av;
 
   if (atomic_energy || atomic_virial) {
     DP_REQUIRES_OK(dp, dp->dp.compute(e, f, v, ae, av, coord_, atype_, cell_,
-                                      fparam_, aparam_));
+                                      fparam_, aparam_, charge_spin_));
   } else {
-    DP_REQUIRES_OK(
-        dp, dp->dp.compute(e, f, v, coord_, atype_, cell_, fparam_, aparam_));
+    DP_REQUIRES_OK(dp, dp->dp.compute(e, f, v, coord_, atype_, cell_, fparam_,
+                                      aparam_, charge_spin_));
   }
   // copy from C++ vectors to C arrays, if not NULL pointer
   if (energy) {
@@ -308,7 +316,8 @@ template void DP_DeepPotCompute_variant<double>(DP_DeepPot* dp,
                                                 double* force,
                                                 double* virial,
                                                 double* atomic_energy,
-                                                double* atomic_virial);
+                                                double* atomic_virial,
+                                                const double* charge_spin);
 
 template void DP_DeepPotCompute_variant<float>(DP_DeepPot* dp,
                                                const int nframes,
@@ -322,7 +331,8 @@ template void DP_DeepPotCompute_variant<float>(DP_DeepPot* dp,
                                                float* force,
                                                float* virial,
                                                float* atomic_energy,
-                                               float* atomic_virial);
+                                               float* atomic_virial,
+                                               const float* charge_spin);
 // support spin
 template <typename VALUETYPE>
 inline void DP_DeepSpinCompute_variant(DP_DeepSpin* dp,
@@ -416,22 +426,24 @@ template void DP_DeepSpinCompute_variant<float>(DP_DeepSpin* dp,
                                                 float* atomic_virial);
 
 template <typename VALUETYPE>
-inline void DP_DeepPotComputeNList_variant(DP_DeepPot* dp,
-                                           const int nframes,
-                                           const int natoms,
-                                           const VALUETYPE* coord,
-                                           const int* atype,
-                                           const VALUETYPE* cell,
-                                           const int nghost,
-                                           const DP_Nlist* nlist,
-                                           const int ago,
-                                           const VALUETYPE* fparam,
-                                           const VALUETYPE* aparam,
-                                           double* energy,
-                                           VALUETYPE* force,
-                                           VALUETYPE* virial,
-                                           VALUETYPE* atomic_energy,
-                                           VALUETYPE* atomic_virial) {
+inline void DP_DeepPotComputeNList_variant(
+    DP_DeepPot* dp,
+    const int nframes,
+    const int natoms,
+    const VALUETYPE* coord,
+    const int* atype,
+    const VALUETYPE* cell,
+    const int nghost,
+    const DP_Nlist* nlist,
+    const int ago,
+    const VALUETYPE* fparam,
+    const VALUETYPE* aparam,
+    double* energy,
+    VALUETYPE* force,
+    VALUETYPE* virial,
+    VALUETYPE* atomic_energy,
+    VALUETYPE* atomic_virial,
+    const VALUETYPE* charge_spin = nullptr) {
   // init C++ vectors from C arrays
   std::vector<VALUETYPE> coord_(coord, coord + nframes * natoms * 3);
   std::vector<int> atype_(atype, atype + natoms);
@@ -451,16 +463,24 @@ inline void DP_DeepPotComputeNList_variant(DP_DeepPot* dp,
                                 (dp->aparam_nall ? natoms : (natoms - nghost)) *
                                 dp->daparam);
   }
+  // charge_spin is converted to double for the api_cc::DeepPot interface
+  // (the compute backend stores/uses charge_spin as float64).
+  std::vector<double> charge_spin_;
+  if (charge_spin) {
+    charge_spin_.assign(charge_spin,
+                        charge_spin + nframes * dp->dp.dim_chg_spin());
+  }
   std::vector<double> e;
   std::vector<VALUETYPE> f, v, ae, av;
 
   if (atomic_energy || atomic_virial) {
     DP_REQUIRES_OK(
         dp, dp->dp.compute(e, f, v, ae, av, coord_, atype_, cell_, nghost,
-                           nlist->nl, ago, fparam_, aparam_));
+                           nlist->nl, ago, fparam_, aparam_, charge_spin_));
   } else {
-    DP_REQUIRES_OK(dp, dp->dp.compute(e, f, v, coord_, atype_, cell_, nghost,
-                                      nlist->nl, ago, fparam_, aparam_));
+    DP_REQUIRES_OK(
+        dp, dp->dp.compute(e, f, v, coord_, atype_, cell_, nghost, nlist->nl,
+                           ago, fparam_, aparam_, charge_spin_));
   }
   // copy from C++ vectors to C arrays, if not NULL pointer
   if (energy) {
@@ -495,7 +515,8 @@ template void DP_DeepPotComputeNList_variant<double>(DP_DeepPot* dp,
                                                      double* force,
                                                      double* virial,
                                                      double* atomic_energy,
-                                                     double* atomic_virial);
+                                                     double* atomic_virial,
+                                                     const double* charge_spin);
 
 template void DP_DeepPotComputeNList_variant<float>(DP_DeepPot* dp,
                                                     const int nframes,
@@ -512,7 +533,8 @@ template void DP_DeepPotComputeNList_variant<float>(DP_DeepPot* dp,
                                                     float* force,
                                                     float* virial,
                                                     float* atomic_energy,
-                                                    float* atomic_virial);
+                                                    float* atomic_virial,
+                                                    const float* charge_spin);
 
 // support spin
 template <typename VALUETYPE>
@@ -708,19 +730,21 @@ inline void flatten_vector(std::vector<VALUETYPE>& onedv,
 }
 
 template <typename VALUETYPE>
-void DP_DeepPotModelDeviCompute_variant(DP_DeepPotModelDevi* dp,
-                                        const int nframes,
-                                        const int natoms,
-                                        const VALUETYPE* coord,
-                                        const int* atype,
-                                        const VALUETYPE* cell,
-                                        const VALUETYPE* fparam,
-                                        const VALUETYPE* aparam,
-                                        double* energy,
-                                        VALUETYPE* force,
-                                        VALUETYPE* virial,
-                                        VALUETYPE* atomic_energy,
-                                        VALUETYPE* atomic_virial) {
+void DP_DeepPotModelDeviCompute_variant(
+    DP_DeepPotModelDevi* dp,
+    const int nframes,
+    const int natoms,
+    const VALUETYPE* coord,
+    const int* atype,
+    const VALUETYPE* cell,
+    const VALUETYPE* fparam,
+    const VALUETYPE* aparam,
+    double* energy,
+    VALUETYPE* force,
+    VALUETYPE* virial,
+    VALUETYPE* atomic_energy,
+    VALUETYPE* atomic_virial,
+    const VALUETYPE* charge_spin = nullptr) {
   if (nframes > 1) {
     throw std::runtime_error("nframes > 1 not supported yet");
   }
@@ -740,16 +764,23 @@ void DP_DeepPotModelDeviCompute_variant(DP_DeepPotModelDevi* dp,
   if (aparam) {
     aparam_.assign(aparam, aparam + nframes * natoms * dp->daparam);
   }
+  // charge_spin is converted to double for the api_cc::DeepPot interface
+  // (the compute backend stores/uses charge_spin as float64).
+  std::vector<double> charge_spin_;
+  if (charge_spin) {
+    charge_spin_.assign(charge_spin,
+                        charge_spin + nframes * dp->dp.dim_chg_spin());
+  }
   // different from DeepPot
   std::vector<double> e;
   std::vector<std::vector<VALUETYPE>> f, v, ae, av;
 
   if (atomic_energy || atomic_virial) {
     DP_REQUIRES_OK(dp, dp->dp.compute(e, f, v, ae, av, coord_, atype_, cell_,
-                                      fparam_, aparam_));
+                                      fparam_, aparam_, charge_spin_));
   } else {
-    DP_REQUIRES_OK(
-        dp, dp->dp.compute(e, f, v, coord_, atype_, cell_, fparam_, aparam_));
+    DP_REQUIRES_OK(dp, dp->dp.compute(e, f, v, coord_, atype_, cell_, fparam_,
+                                      aparam_, charge_spin_));
   }
   // 2D vector to 2D array, flatten first
   if (energy) {
@@ -790,21 +821,24 @@ template void DP_DeepPotModelDeviCompute_variant<double>(
     double* force,
     double* virial,
     double* atomic_energy,
-    double* atomic_virial);
+    double* atomic_virial,
+    const double* charge_spin);
 
-template void DP_DeepPotModelDeviCompute_variant<float>(DP_DeepPotModelDevi* dp,
-                                                        const int nframes,
-                                                        const int natoms,
-                                                        const float* coord,
-                                                        const int* atype,
-                                                        const float* cell,
-                                                        const float* fparam,
-                                                        const float* aparam,
-                                                        double* energy,
-                                                        float* force,
-                                                        float* virial,
-                                                        float* atomic_energy,
-                                                        float* atomic_virial);
+template void DP_DeepPotModelDeviCompute_variant<float>(
+    DP_DeepPotModelDevi* dp,
+    const int nframes,
+    const int natoms,
+    const float* coord,
+    const int* atype,
+    const float* cell,
+    const float* fparam,
+    const float* aparam,
+    double* energy,
+    float* force,
+    float* virial,
+    float* atomic_energy,
+    float* atomic_virial,
+    const float* charge_spin);
 
 template <typename VALUETYPE>
 void DP_DeepSpinModelDeviCompute_variant(DP_DeepSpinModelDevi* dp,
@@ -919,22 +953,24 @@ template void DP_DeepSpinModelDeviCompute_variant<float>(
     float* atomic_virial);
 
 template <typename VALUETYPE>
-void DP_DeepPotModelDeviComputeNList_variant(DP_DeepPotModelDevi* dp,
-                                             const int nframes,
-                                             const int natoms,
-                                             const VALUETYPE* coord,
-                                             const int* atype,
-                                             const VALUETYPE* cell,
-                                             const int nghost,
-                                             const DP_Nlist* nlist,
-                                             const int ago,
-                                             const VALUETYPE* fparam,
-                                             const VALUETYPE* aparam,
-                                             double* energy,
-                                             VALUETYPE* force,
-                                             VALUETYPE* virial,
-                                             VALUETYPE* atomic_energy,
-                                             VALUETYPE* atomic_virial) {
+void DP_DeepPotModelDeviComputeNList_variant(
+    DP_DeepPotModelDevi* dp,
+    const int nframes,
+    const int natoms,
+    const VALUETYPE* coord,
+    const int* atype,
+    const VALUETYPE* cell,
+    const int nghost,
+    const DP_Nlist* nlist,
+    const int ago,
+    const VALUETYPE* fparam,
+    const VALUETYPE* aparam,
+    double* energy,
+    VALUETYPE* force,
+    VALUETYPE* virial,
+    VALUETYPE* atomic_energy,
+    VALUETYPE* atomic_virial,
+    const VALUETYPE* charge_spin = nullptr) {
   if (nframes > 1) {
     throw std::runtime_error("nframes > 1 not supported yet");
   }
@@ -956,6 +992,13 @@ void DP_DeepPotModelDeviComputeNList_variant(DP_DeepPotModelDevi* dp,
         aparam,
         aparam + (dp->aparam_nall ? natoms : (natoms - nghost)) * dp->daparam);
   }
+  // charge_spin is converted to double for the api_cc::DeepPot interface
+  // (the compute backend stores/uses charge_spin as float64).
+  std::vector<double> charge_spin_;
+  if (charge_spin) {
+    charge_spin_.assign(charge_spin,
+                        charge_spin + nframes * dp->dp.dim_chg_spin());
+  }
   // different from DeepPot
   std::vector<double> e;
   std::vector<std::vector<VALUETYPE>> f, v, ae, av;
@@ -963,10 +1006,11 @@ void DP_DeepPotModelDeviComputeNList_variant(DP_DeepPotModelDevi* dp,
   if (atomic_energy || atomic_virial) {
     DP_REQUIRES_OK(
         dp, dp->dp.compute(e, f, v, ae, av, coord_, atype_, cell_, nghost,
-                           nlist->nl, ago, fparam_, aparam_));
+                           nlist->nl, ago, fparam_, aparam_, charge_spin_));
   } else {
-    DP_REQUIRES_OK(dp, dp->dp.compute(e, f, v, coord_, atype_, cell_, nghost,
-                                      nlist->nl, ago, fparam_, aparam_));
+    DP_REQUIRES_OK(
+        dp, dp->dp.compute(e, f, v, coord_, atype_, cell_, nghost, nlist->nl,
+                           ago, fparam_, aparam_, charge_spin_));
   }
   // 2D vector to 2D array, flatten first
   if (energy) {
@@ -1010,7 +1054,8 @@ template void DP_DeepPotModelDeviComputeNList_variant<double>(
     double* force,
     double* virial,
     double* atomic_energy,
-    double* atomic_virial);
+    double* atomic_virial,
+    const double* charge_spin);
 
 template void DP_DeepPotModelDeviComputeNList_variant<float>(
     DP_DeepPotModelDevi* dp,
@@ -1028,7 +1073,8 @@ template void DP_DeepPotModelDeviComputeNList_variant<float>(
     float* force,
     float* virial,
     float* atomic_energy,
-    float* atomic_virial);
+    float* atomic_virial,
+    const float* charge_spin);
 
 // support spin multi model.
 template <typename VALUETYPE>
@@ -1702,6 +1748,90 @@ void DP_DeepPotComputeNListf2(DP_DeepPot* dp,
       aparam, energy, force, virial, atomic_energy, atomic_virial);
 }
 
+// charge_spin-aware variants (version 3): same as the version-2 functions
+// plus a per-frame charge_spin input (nframes x dim_chg_spin, always double).
+void DP_DeepPotCompute3(DP_DeepPot* dp,
+                        const int nframes,
+                        const int natoms,
+                        const double* coord,
+                        const int* atype,
+                        const double* cell,
+                        const double* fparam,
+                        const double* aparam,
+                        const double* charge_spin,
+                        double* energy,
+                        double* force,
+                        double* virial,
+                        double* atomic_energy,
+                        double* atomic_virial) {
+  DP_DeepPotCompute_variant<double>(dp, nframes, natoms, coord, atype, cell,
+                                    fparam, aparam, energy, force, virial,
+                                    atomic_energy, atomic_virial, charge_spin);
+}
+
+void DP_DeepPotComputef3(DP_DeepPot* dp,
+                         const int nframes,
+                         const int natoms,
+                         const float* coord,
+                         const int* atype,
+                         const float* cell,
+                         const float* fparam,
+                         const float* aparam,
+                         const float* charge_spin,
+                         double* energy,
+                         float* force,
+                         float* virial,
+                         float* atomic_energy,
+                         float* atomic_virial) {
+  DP_DeepPotCompute_variant<float>(dp, nframes, natoms, coord, atype, cell,
+                                   fparam, aparam, energy, force, virial,
+                                   atomic_energy, atomic_virial, charge_spin);
+}
+
+void DP_DeepPotComputeNList3(DP_DeepPot* dp,
+                             const int nframes,
+                             const int natoms,
+                             const double* coord,
+                             const int* atype,
+                             const double* cell,
+                             const int nghost,
+                             const DP_Nlist* nlist,
+                             const int ago,
+                             const double* fparam,
+                             const double* aparam,
+                             const double* charge_spin,
+                             double* energy,
+                             double* force,
+                             double* virial,
+                             double* atomic_energy,
+                             double* atomic_virial) {
+  DP_DeepPotComputeNList_variant<double>(
+      dp, nframes, natoms, coord, atype, cell, nghost, nlist, ago, fparam,
+      aparam, energy, force, virial, atomic_energy, atomic_virial, charge_spin);
+}
+
+void DP_DeepPotComputeNListf3(DP_DeepPot* dp,
+                              const int nframes,
+                              const int natoms,
+                              const float* coord,
+                              const int* atype,
+                              const float* cell,
+                              const int nghost,
+                              const DP_Nlist* nlist,
+                              const int ago,
+                              const float* fparam,
+                              const float* aparam,
+                              const float* charge_spin,
+                              double* energy,
+                              float* force,
+                              float* virial,
+                              float* atomic_energy,
+                              float* atomic_virial) {
+  DP_DeepPotComputeNList_variant<float>(
+      dp, nframes, natoms, coord, atype, cell, nghost, nlist, ago, fparam,
+      aparam, energy, force, virial, atomic_energy, atomic_virial, charge_spin);
+}
+
 void DP_DeepSpinComputeNListf2(DP_DeepSpin* dp,
                                const int nframes,
                                const int natoms,
@@ -1970,6 +2100,89 @@ void DP_DeepPotModelDeviComputeNListf2(DP_DeepPotModelDevi* dp,
       aparam, energy, force, virial, atomic_energy, atomic_virial);
 }
 
+// charge_spin-aware model-deviation variants (version 3).
+void DP_DeepPotModelDeviCompute3(DP_DeepPotModelDevi* dp,
+                                 const int nframes,
+                                 const int natoms,
+                                 const double* coord,
+                                 const int* atype,
+                                 const double* cell,
+                                 const double* fparam,
+                                 const double* aparam,
+                                 const double* charge_spin,
+                                 double* energy,
+                                 double* force,
+                                 double* virial,
+                                 double* atomic_energy,
+                                 double* atomic_virial) {
+  DP_DeepPotModelDeviCompute_variant<double>(
+      dp, nframes, natoms, coord, atype, cell, fparam, aparam, energy, force,
+      virial, atomic_energy, atomic_virial, charge_spin);
+}
+
+void DP_DeepPotModelDeviComputef3(DP_DeepPotModelDevi* dp,
+                                  const int nframes,
+                                  const int natoms,
+                                  const float* coord,
+                                  const int* atype,
+                                  const float* cell,
+                                  const float* fparam,
+                                  const float* aparam,
+                                  const float* charge_spin,
+                                  double* energy,
+                                  float* force,
+                                  float* virial,
+                                  float* atomic_energy,
+                                  float* atomic_virial) {
+  DP_DeepPotModelDeviCompute_variant<float>(
+      dp, nframes, natoms, coord, atype, cell, fparam, aparam, energy, force,
+      virial, atomic_energy, atomic_virial, charge_spin);
+}
+
+void DP_DeepPotModelDeviComputeNList3(DP_DeepPotModelDevi* dp,
+                                      const int nframes,
+                                      const int natoms,
+                                      const double* coord,
+                                      const int* atype,
+                                      const double* cell,
+                                      const int nghost,
+                                      const DP_Nlist* nlist,
+                                      const int ago,
+                                      const double* fparam,
+                                      const double* aparam,
+                                      const double* charge_spin,
+                                      double* energy,
+                                      double* force,
+                                      double* virial,
+                                      double* atomic_energy,
+                                      double* atomic_virial) {
+  DP_DeepPotModelDeviComputeNList_variant<double>(
+      dp, nframes, natoms, coord, atype, cell, nghost, nlist, ago, fparam,
+      aparam, energy, force, virial, atomic_energy, atomic_virial, charge_spin);
+}
+
+void DP_DeepPotModelDeviComputeNListf3(DP_DeepPotModelDevi* dp,
+                                       const int nframes,
+                                       const int natoms,
+                                       const float* coord,
+                                       const int* atype,
+                                       const float* cell,
+                                       const int nghost,
+                                       const DP_Nlist* nlist,
+                                       const int ago,
+                                       const float* fparam,
+                                       const float* aparam,
+                                       const float* charge_spin,
+                                       double* energy,
+                                       float* force,
+                                       float* virial,
+                                       float* atomic_energy,
+                                       float* atomic_virial) {
+  DP_DeepPotModelDeviComputeNList_variant<float>(
+      dp, nframes, natoms, coord, atype, cell, nghost, nlist, ago, fparam,
+      aparam, energy, force, virial, atomic_energy, atomic_virial, charge_spin);
+}
+
 void DP_DeepSpinModelDeviComputeNListf2(DP_DeepSpinModelDevi* dp,
                                         const int nframes,
                                         const int natoms,
@@ -2089,6 +2302,8 @@ int DP_DeepPotGetDimAParam(DP_DeepPot* dp) {
   return DP_DeepBaseModelGetDimAParam(static_cast<DP_DeepBaseModel*>(dp));
 }
 
+int DP_DeepPotGetDimChgSpin(DP_DeepPot* dp) { return dp->dp.dim_chg_spin(); }
+
 bool DP_DeepPotIsAParamNAll(DP_DeepPot* dp) {
   return DP_DeepBaseModelIsAParamNAll(static_cast<DP_DeepBaseModel*>(dp));
 }
@@ -2125,6 +2340,10 @@ int DP_DeepPotModelDeviGetDimAParam(DP_DeepPotModelDevi* dp) {
       static_cast<DP_DeepBaseModelDevi*>(dp));
 }
 
+int DP_DeepPotModelDeviGetDimChgSpin(DP_DeepPotModelDevi* dp) {
+  return dp->dp.dim_chg_spin();
+}
+
 bool DP_DeepPotModelDeviIsAParamNAll(DP_DeepPotModelDevi* dp) {
   return DP_DeepBaseModelDeviIsAParamNAll(
       static_cast<DP_DeepBaseModelDevi*>(dp));
diff --git a/source/api_cc/include/DeepPot.h b/source/api_cc/include/DeepPot.h
index 6a30eed7ea..177513ebb0 100644
--- a/source/api_cc/include/DeepPot.h
+++ b/source/api_cc/include/DeepPot.h
@@ -198,6 +198,115 @@ class DeepPotBackend : public DeepBaseModelBackend {
                                    const std::vector<float>& aparam,
                                    const bool atomic) = 0;
   /** @} */
+
+  /**
+   * @brief Get dimension of charge/spin condition inputs.
+   * Returns 0 for backends that do not support charge/spin conditioning.
+   **/
+  virtual int dim_chg_spin() const { return 0; }
+
+  // charge_spin-aware computew overloads.  Default implementations call the
+  // existing pure-virtual overloads (ignoring charge_spin) so that backends
+  // that do not support charge/spin do not need any changes.  DeepPotPTExpt
+  // overrides these to thread charge_spin through to the model.
+  virtual void computew(std::vector<double>& ener,
+                        std::vector<double>& force,
+                        std::vector<double>& virial,
+                        std::vector<double>& atom_energy,
+                        std::vector<double>& atom_virial,
+                        const std::vector<double>& coord,
+                        const std::vector<int>& atype,
+                        const std::vector<double>& box,
+                        const std::vector<double>& fparam,
+                        const std::vector<double>& aparam,
+                        const std::vector<double>& charge_spin,
+                        const bool atomic) {
+    computew(ener, force, virial, atom_energy, atom_virial, coord, atype, box,
+             fparam, aparam, atomic);
+  }
+  virtual void computew(std::vector<double>& ener,
+                        std::vector<float>& force,
+                        std::vector<float>& virial,
+                        std::vector<float>& atom_energy,
+                        std::vector<float>& atom_virial,
+                        const std::vector<float>& coord,
+                        const std::vector<int>& atype,
+                        const std::vector<float>& box,
+                        const std::vector<float>& fparam,
+                        const std::vector<float>& aparam,
+                        const std::vector<double>& charge_spin,
+                        const bool atomic) {
+    computew(ener, force, virial, atom_energy, atom_virial, coord, atype, box,
+             fparam, aparam, atomic);
+  }
+  virtual void computew(std::vector<double>& ener,
+                        std::vector<double>& force,
+                        std::vector<double>& virial,
+                        std::vector<double>& atom_energy,
+                        std::vector<double>& atom_virial,
+                        const std::vector<double>& coord,
+                        const std::vector<int>& atype,
+                        const std::vector<double>& box,
+                        const int nghost,
+                        const InputNlist& inlist,
+                        const int& ago,
+                        const std::vector<double>& fparam,
+                        const std::vector<double>& aparam,
+                        const std::vector<double>& charge_spin,
+                        const bool atomic) {
+    computew(ener, force, virial, atom_energy, atom_virial, coord, atype, box,
+             nghost, inlist, ago, fparam, aparam, atomic);
+  }
+  virtual void computew(std::vector<double>& ener,
+                        std::vector<float>& force,
+                        std::vector<float>& virial,
+                        std::vector<float>& atom_energy,
+                        std::vector<float>& atom_virial,
+                        const std::vector<float>& coord,
+                        const std::vector<int>& atype,
+                        const std::vector<float>& box,
+                        const int nghost,
+                        const InputNlist& inlist,
+                        const int& ago,
+                        const std::vector<float>& fparam,
+                        const std::vector<float>& aparam,
+                        const std::vector<double>& charge_spin,
+                        const bool atomic) {
+    computew(ener, force, virial, atom_energy, atom_virial, coord, atype, box,
+             nghost, inlist, ago, fparam, aparam, atomic);
+  }
+  virtual void computew_mixed_type(std::vector<double>& ener,
+                                   std::vector<double>& force,
+                                   std::vector<double>& virial,
+                                   std::vector<double>& atom_energy,
+                                   std::vector<double>& atom_virial,
+                                   const int& nframes,
+                                   const std::vector<double>& coord,
+                                   const std::vector<int>& atype,
+                                   const std::vector<double>& box,
+                                   const std::vector<double>& fparam,
+                                   const std::vector<double>& aparam,
+                                   const std::vector<double>& charge_spin,
+                                   const bool atomic) {
+    computew_mixed_type(ener, force, virial, atom_energy, atom_virial, nframes,
+                        coord, atype, box, fparam, aparam, atomic);
+  }
+  virtual void computew_mixed_type(std::vector<double>& ener,
+                                   std::vector<float>& force,
+                                   std::vector<float>& virial,
+                                   std::vector<float>& atom_energy,
+                                   std::vector<float>& atom_virial,
+                                   const int& nframes,
+                                   const std::vector<float>& coord,
+                                   const std::vector<int>& atype,
+                                   const std::vector<float>& box,
+                                   const std::vector<float>& fparam,
+                                   const std::vector<float>& aparam,
+                                   const std::vector<double>& charge_spin,
+                                   const bool atomic) {
+    computew_mixed_type(ener, force, virial, atom_energy, atom_virial, nframes,
+                        coord, atype, box, fparam, aparam, atomic);
+  }
 };
 
 /**
@@ -249,6 +358,10 @@ class DeepPot : public DeepBaseModel {
    * nframes x natoms x dim_aparam.
    * natoms x dim_aparam. Then all frames are assumed to be provided with the
    *same aparam.
+   * @param[in] charge_spin The charge/spin parameter. The array can be of size
+   *nframes x dim_chg_spin.
+   * dim_chg_spin. Then all frames are assumed to be provided with the same
+   *charge_spin. Leave it empty to use the model's stored default_chg_spin.
    * @{
    **/
   template <typename VALUETYPE>
@@ -259,7 +372,8 @@ class DeepPot : public DeepBaseModel {
                const std::vector<int>& atype,
                const std::vector<VALUETYPE>& box,
                const std::vector<VALUETYPE>& fparam = std::vector<VALUETYPE>(),
-               const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>());
+               const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>(),
+               const std::vector<double>& charge_spin = std::vector<double>());
   template <typename VALUETYPE>
   void compute(std::vector<ENERGYTYPE>& ener,
                std::vector<VALUETYPE>& force,
@@ -268,7 +382,8 @@ class DeepPot : public DeepBaseModel {
                const std::vector<int>& atype,
                const std::vector<VALUETYPE>& box,
                const std::vector<VALUETYPE>& fparam = std::vector<VALUETYPE>(),
-               const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>());
+               const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>(),
+               const std::vector<double>& charge_spin = std::vector<double>());
   /** @} */
   /**
    * @brief Evaluate the energy, force and virial by using this DP.
@@ -291,6 +406,10 @@ class DeepPot : public DeepBaseModel {
    * nframes x natoms x dim_aparam.
    * natoms x dim_aparam. Then all frames are assumed to be provided with the
    *same aparam.
+   * @param[in] charge_spin The charge/spin parameter. The array can be of size
+   *nframes x dim_chg_spin.
+   * dim_chg_spin. Then all frames are assumed to be provided with the same
+   *charge_spin. Leave it empty to use the model's stored default_chg_spin.
    * @{
    **/
   template <typename VALUETYPE>
@@ -304,7 +423,8 @@ class DeepPot : public DeepBaseModel {
                const InputNlist& inlist,
                const int& ago,
                const std::vector<VALUETYPE>& fparam = std::vector<VALUETYPE>(),
-               const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>());
+               const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>(),
+               const std::vector<double>& charge_spin = std::vector<double>());
   template <typename VALUETYPE>
   void compute(std::vector<ENERGYTYPE>& ener,
                std::vector<VALUETYPE>& force,
@@ -316,7 +436,8 @@ class DeepPot : public DeepBaseModel {
                const InputNlist& inlist,
                const int& ago,
                const std::vector<VALUETYPE>& fparam = std::vector<VALUETYPE>(),
-               const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>());
+               const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>(),
+               const std::vector<double>& charge_spin = std::vector<double>());
   /** @} */
   /**
    * @brief Evaluate the energy, force, virial, atomic energy, and atomic virial
@@ -339,6 +460,10 @@ class DeepPot : public DeepBaseModel {
    * nframes x natoms x dim_aparam.
    * natoms x dim_aparam. Then all frames are assumed to be provided with the
    *same aparam.
+   * @param[in] charge_spin The charge/spin parameter. The array can be of size
+   *nframes x dim_chg_spin.
+   * dim_chg_spin. Then all frames are assumed to be provided with the same
+   *charge_spin. Leave it empty to use the model's stored default_chg_spin.
    * @{
    **/
   template <typename VALUETYPE>
@@ -351,7 +476,8 @@ class DeepPot : public DeepBaseModel {
                const std::vector<int>& atype,
                const std::vector<VALUETYPE>& box,
                const std::vector<VALUETYPE>& fparam = std::vector<VALUETYPE>(),
-               const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>());
+               const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>(),
+               const std::vector<double>& charge_spin = std::vector<double>());
   template <typename VALUETYPE>
   void compute(std::vector<ENERGYTYPE>& ener,
                std::vector<VALUETYPE>& force,
@@ -362,7 +488,8 @@ class DeepPot : public DeepBaseModel {
                const std::vector<int>& atype,
                const std::vector<VALUETYPE>& box,
                const std::vector<VALUETYPE>& fparam = std::vector<VALUETYPE>(),
-               const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>());
+               const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>(),
+               const std::vector<double>& charge_spin = std::vector<double>());
   /** @} */
 
   /**
@@ -389,6 +516,10 @@ class DeepPot : public DeepBaseModel {
    * nframes x natoms x dim_aparam.
    * natoms x dim_aparam. Then all frames are assumed to be provided with the
    *same aparam.
+   * @param[in] charge_spin The charge/spin parameter. The array can be of size
+   *nframes x dim_chg_spin.
+   * dim_chg_spin. Then all frames are assumed to be provided with the same
+   *charge_spin. Leave it empty to use the model's stored default_chg_spin.
    * @{
    **/
   template <typename VALUETYPE>
@@ -404,7 +535,8 @@ class DeepPot : public DeepBaseModel {
                const InputNlist& lmp_list,
                const int& ago,
                const std::vector<VALUETYPE>& fparam = std::vector<VALUETYPE>(),
-               const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>());
+               const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>(),
+               const std::vector<double>& charge_spin = std::vector<double>());
   template <typename VALUETYPE>
   void compute(std::vector<ENERGYTYPE>& ener,
                std::vector<VALUETYPE>& force,
@@ -418,7 +550,8 @@ class DeepPot : public DeepBaseModel {
                const InputNlist& lmp_list,
                const int& ago,
                const std::vector<VALUETYPE>& fparam = std::vector<VALUETYPE>(),
-               const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>());
+               const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>(),
+               const std::vector<double>& charge_spin = std::vector<double>());
   /** @} */
   /**
    * @brief Evaluate the energy, force, and virial with the mixed type
@@ -441,6 +574,10 @@ class DeepPot : public DeepBaseModel {
    * nframes x natoms x dim_aparam.
    * natoms x dim_aparam. Then all frames are assumed to be provided with the
    *same aparam.
+   * @param[in] charge_spin The charge/spin parameter. The array can be of size
+   *nframes x dim_chg_spin.
+   * dim_chg_spin. Then all frames are assumed to be provided with the same
+   *charge_spin. Leave it empty to use the model's stored default_chg_spin.
    * @{
    **/
   template <typename VALUETYPE>
@@ -453,7 +590,8 @@ class DeepPot : public DeepBaseModel {
       const std::vector<int>& atype,
       const std::vector<VALUETYPE>& box,
       const std::vector<VALUETYPE>& fparam = std::vector<VALUETYPE>(),
-      const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>());
+      const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>(),
+      const std::vector<double>& charge_spin = std::vector<double>());
   template <typename VALUETYPE>
   void compute_mixed_type(
       std::vector<ENERGYTYPE>& ener,
@@ -464,7 +602,8 @@ class DeepPot : public DeepBaseModel {
       const std::vector<int>& atype,
       const std::vector<VALUETYPE>& box,
       const std::vector<VALUETYPE>& fparam = std::vector<VALUETYPE>(),
-      const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>());
+      const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>(),
+      const std::vector<double>& charge_spin = std::vector<double>());
   /** @} */
   /**
    * @brief Evaluate the energy, force, and virial with the mixed type
@@ -489,6 +628,10 @@ class DeepPot : public DeepBaseModel {
    * nframes x natoms x dim_aparam.
    * natoms x dim_aparam. Then all frames are assumed to be provided with the
    *same aparam.
+   * @param[in] charge_spin The charge/spin parameter. The array can be of size
+   *nframes x dim_chg_spin.
+   * dim_chg_spin. Then all frames are assumed to be provided with the same
+   *charge_spin. Leave it empty to use the model's stored default_chg_spin.
    * @{
    **/
   template <typename VALUETYPE>
@@ -503,7 +646,8 @@ class DeepPot : public DeepBaseModel {
       const std::vector<int>& atype,
       const std::vector<VALUETYPE>& box,
       const std::vector<VALUETYPE>& fparam = std::vector<VALUETYPE>(),
-      const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>());
+      const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>(),
+      const std::vector<double>& charge_spin = std::vector<double>());
   template <typename VALUETYPE>
   void compute_mixed_type(
       std::vector<ENERGYTYPE>& ener,
@@ -516,8 +660,12 @@ class DeepPot : public DeepBaseModel {
       const std::vector<int>& atype,
       const std::vector<VALUETYPE>& box,
       const std::vector<VALUETYPE>& fparam = std::vector<VALUETYPE>(),
-      const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>());
+      const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>(),
+      const std::vector<double>& charge_spin = std::vector<double>());
   /** @} */
+
+  int dim_chg_spin() const;
+
  protected:
   std::shared_ptr<deepmd::DeepPotBackend> dp;
 };
@@ -552,6 +700,16 @@ class DeepPotModelDevi : public DeepBaseModelDevi {
             const std::vector<std::string>& file_contents =
                 std::vector<std::string>());
 
+  /**
+   * @brief Get the dimension of the charge/spin input.
+   * @return The dimension of the charge/spin input (0 if the models have no
+   *charge/spin embedding). Taken from the first model; all models are assumed
+   *to share the same value.
+   **/
+  int dim_chg_spin() const {
+    return numb_models > 0 ? dps[0]->dim_chg_spin() : 0;
+  };
+
   /**
    * @brief Evaluate the energy, force and virial by using these DP models.
    * @param[out] all_ener The system energies of all models.
@@ -571,6 +729,10 @@ class DeepPotModelDevi : public DeepBaseModelDevi {
    * natoms x dim_aparam. Then all frames are assumed to be provided with the
    *same aparam. dim_aparam. Then all frames and atoms are provided with the
    *same aparam.
+   * @param[in] charge_spin The charge/spin parameter. The array can be of size
+   *nframes x dim_chg_spin.
+   * dim_chg_spin. Then all frames are assumed to be provided with the same
+   *charge_spin. Leave it empty to use the model's stored default_chg_spin.
    **/
   template <typename VALUETYPE>
   void compute(std::vector<ENERGYTYPE>& all_ener,
@@ -580,7 +742,8 @@ class DeepPotModelDevi : public DeepBaseModelDevi {
                const std::vector<int>& atype,
                const std::vector<VALUETYPE>& box,
                const std::vector<VALUETYPE>& fparam = std::vector<VALUETYPE>(),
-               const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>());
+               const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>(),
+               const std::vector<double>& charge_spin = std::vector<double>());
 
   /**
    * @brief Evaluate the energy, force, virial, atomic energy, and atomic virial
@@ -604,6 +767,10 @@ class DeepPotModelDevi : public DeepBaseModelDevi {
    * natoms x dim_aparam. Then all frames are assumed to be provided with the
    *same aparam. dim_aparam. Then all frames and atoms are provided with the
    *same aparam.
+   * @param[in] charge_spin The charge/spin parameter. The array can be of size
+   *nframes x dim_chg_spin.
+   * dim_chg_spin. Then all frames are assumed to be provided with the same
+   *charge_spin. Leave it empty to use the model's stored default_chg_spin.
    **/
   template <typename VALUETYPE>
   void compute(std::vector<ENERGYTYPE>& all_ener,
@@ -615,7 +782,8 @@ class DeepPotModelDevi : public DeepBaseModelDevi {
                const std::vector<int>& atype,
                const std::vector<VALUETYPE>& box,
                const std::vector<VALUETYPE>& fparam = std::vector<VALUETYPE>(),
-               const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>());
+               const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>(),
+               const std::vector<double>& charge_spin = std::vector<double>());
 
   /**
    * @brief Evaluate the energy, force and virial by using these DP models.
@@ -639,6 +807,10 @@ class DeepPotModelDevi : public DeepBaseModelDevi {
    * natoms x dim_aparam. Then all frames are assumed to be provided with the
    *same aparam. dim_aparam. Then all frames and atoms are provided with the
    *same aparam.
+   * @param[in] charge_spin The charge/spin parameter. The array can be of size
+   *nframes x dim_chg_spin.
+   * dim_chg_spin. Then all frames are assumed to be provided with the same
+   *charge_spin. Leave it empty to use the model's stored default_chg_spin.
    **/
   template <typename VALUETYPE>
   void compute(std::vector<ENERGYTYPE>& all_ener,
@@ -651,7 +823,8 @@ class DeepPotModelDevi : public DeepBaseModelDevi {
                const InputNlist& lmp_list,
                const int& ago,
                const std::vector<VALUETYPE>& fparam = std::vector<VALUETYPE>(),
-               const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>());
+               const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>(),
+               const std::vector<double>& charge_spin = std::vector<double>());
   /**
    * @brief Evaluate the energy, force, virial, atomic energy, and atomic virial
    *by using these DP models.
@@ -677,6 +850,10 @@ class DeepPotModelDevi : public DeepBaseModelDevi {
    * natoms x dim_aparam. Then all frames are assumed to be provided with the
    *same aparam. dim_aparam. Then all frames and atoms are provided with the
    *same aparam.
+   * @param[in] charge_spin The charge/spin parameter. The array can be of size
+   *nframes x dim_chg_spin.
+   * dim_chg_spin. Then all frames are assumed to be provided with the same
+   *charge_spin. Leave it empty to use the model's stored default_chg_spin.
    **/
   template <typename VALUETYPE>
   void compute(std::vector<ENERGYTYPE>& all_ener,
@@ -691,7 +868,8 @@ class DeepPotModelDevi : public DeepBaseModelDevi {
                const InputNlist& lmp_list,
                const int& ago,
                const std::vector<VALUETYPE>& fparam = std::vector<VALUETYPE>(),
-               const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>());
+               const std::vector<VALUETYPE>& aparam = std::vector<VALUETYPE>(),
+               const std::vector<double>& charge_spin = std::vector<double>());
 
  protected:
   std::vector<std::shared_ptr<deepmd::DeepPot>> dps;
diff --git a/source/api_cc/include/DeepPotPT.h b/source/api_cc/include/DeepPotPT.h
index 14b545dce4..deb487610f 100644
--- a/source/api_cc/include/DeepPotPT.h
+++ b/source/api_cc/include/DeepPotPT.h
@@ -73,6 +73,7 @@ class DeepPotPT : public DeepPotBackend {
                const std::vector<VALUETYPE>& box,
                const std::vector<VALUETYPE>& fparam,
                const std::vector<VALUETYPE>& aparam,
+               const std::vector<double>& charge_spin,
                const bool atomic);
   /**
    * @brief Evaluate the energy, force, virial, atomic energy, and atomic virial
@@ -114,6 +115,7 @@ class DeepPotPT : public DeepPotBackend {
                const int& ago,
                const std::vector<VALUETYPE>& fparam,
                const std::vector<VALUETYPE>& aparam,
+               const std::vector<double>& charge_spin,
                const bool atomic);
   /**
    * @brief Evaluate the energy, force, and virial with the mixed type
@@ -229,6 +231,15 @@ class DeepPotPT : public DeepPotBackend {
     assert(inited);
     return daparam;
   };
+  /**
+   * @brief Get the dimension of the charge/spin input.
+   * @return The dimension of the charge/spin input (0 if the model has no
+   *charge/spin embedding).
+   **/
+  int dim_chg_spin() const override {
+    assert(inited);
+    return dchgspin;
+  };
   /**
    * @brief Get the type map (element name of the atom types) of this model.
    * @param[out] type_map The type map of this model.
@@ -304,6 +315,11 @@ class DeepPotPT : public DeepPotBackend {
                 const std::vector<float>& fparam,
                 const std::vector<float>& aparam,
                 const bool atomic);
+  // Keep the base charge_spin-aware computew / computew_mixed_type overloads
+  // visible: the non-charge_spin declarations below would otherwise hide them
+  // for callers using the static type DeepPotPT.
+  using DeepPotBackend::computew;
+  using DeepPotBackend::computew_mixed_type;
   void computew_mixed_type(std::vector<double>& ener,
                            std::vector<double>& force,
                            std::vector<double>& virial,
@@ -329,6 +345,62 @@ class DeepPotPT : public DeepPotBackend {
                            const std::vector<float>& aparam,
                            const bool atomic);
 
+  // charge_spin overloads — pass runtime charge/spin per call
+  void computew(std::vector<double>& ener,
+                std::vector<double>& force,
+                std::vector<double>& virial,
+                std::vector<double>& atom_energy,
+                std::vector<double>& atom_virial,
+                const std::vector<double>& coord,
+                const std::vector<int>& atype,
+                const std::vector<double>& box,
+                const std::vector<double>& fparam,
+                const std::vector<double>& aparam,
+                const std::vector<double>& charge_spin,
+                const bool atomic) override;
+  void computew(std::vector<double>& ener,
+                std::vector<float>& force,
+                std::vector<float>& virial,
+                std::vector<float>& atom_energy,
+                std::vector<float>& atom_virial,
+                const std::vector<float>& coord,
+                const std::vector<int>& atype,
+                const std::vector<float>& box,
+                const std::vector<float>& fparam,
+                const std::vector<float>& aparam,
+                const std::vector<double>& charge_spin,
+                const bool atomic) override;
+  void computew(std::vector<double>& ener,
+                std::vector<double>& force,
+                std::vector<double>& virial,
+                std::vector<double>& atom_energy,
+                std::vector<double>& atom_virial,
+                const std::vector<double>& coord,
+                const std::vector<int>& atype,
+                const std::vector<double>& box,
+                const int nghost,
+                const InputNlist& inlist,
+                const int& ago,
+                const std::vector<double>& fparam,
+                const std::vector<double>& aparam,
+                const std::vector<double>& charge_spin,
+                const bool atomic) override;
+  void computew(std::vector<double>& ener,
+                std::vector<float>& force,
+                std::vector<float>& virial,
+                std::vector<float>& atom_energy,
+                std::vector<float>& atom_virial,
+                const std::vector<float>& coord,
+                const std::vector<int>& atype,
+                const std::vector<float>& box,
+                const int nghost,
+                const InputNlist& inlist,
+                const int& ago,
+                const std::vector<float>& fparam,
+                const std::vector<float>& aparam,
+                const std::vector<double>& charge_spin,
+                const bool atomic) override;
+
  private:
   int num_intra_nthreads, num_inter_nthreads;
   bool inited;
@@ -336,8 +408,10 @@ class DeepPotPT : public DeepPotBackend {
   int ntypes_spin;
   int dfparam;
   int daparam;
+  int dchgspin;
   bool aparam_nall;
   bool has_default_fparam_;
+  std::vector<double> default_chg_spin_;
   // copy neighbor list info from host
   torch::jit::script::Module module;
   double rcut;
diff --git a/source/api_cc/include/DeepPotPTExpt.h b/source/api_cc/include/DeepPotPTExpt.h
index 94631ab12a..461fc6f33c 100644
--- a/source/api_cc/include/DeepPotPTExpt.h
+++ b/source/api_cc/include/DeepPotPTExpt.h
@@ -77,6 +77,7 @@ class DeepPotPTExpt : public DeepPotBackend {
                const int& ago,
                const std::vector<VALUETYPE>& fparam,
                const std::vector<VALUETYPE>& aparam,
+               const std::vector<double>& charge_spin,
                const bool atomic);
   /**
    * @brief Evaluate without nlist (standalone — builds nlist, folds back).
@@ -92,6 +93,7 @@ class DeepPotPTExpt : public DeepPotBackend {
                const std::vector<VALUETYPE>& box,
                const std::vector<VALUETYPE>& fparam,
                const std::vector<VALUETYPE>& aparam,
+               const std::vector<double>& charge_spin,
                const bool atomic);
 
  public:
@@ -115,6 +117,10 @@ class DeepPotPTExpt : public DeepPotBackend {
     assert(inited);
     return daparam;
   };
+  int dim_chg_spin() const override {
+    assert(inited);
+    return dchgspin;
+  };
   void get_type_map(std::string& type_map);
   bool is_aparam_nall() const {
     assert(inited);
@@ -125,7 +131,8 @@ class DeepPotPTExpt : public DeepPotBackend {
     return has_default_fparam_;
   };
 
-  // forward to template class
+  // forward to template class (no charge_spin — uses default_chg_spin_
+  // fallback)
   void computew(std::vector<double>& ener,
                 std::vector<double>& force,
                 std::vector<double>& virial,
@@ -201,12 +208,94 @@ class DeepPotPTExpt : public DeepPotBackend {
                            const std::vector<float>& aparam,
                            const bool atomic);
 
+  // charge_spin overloads — pass runtime charge/spin per call
+  void computew(std::vector<double>& ener,
+                std::vector<double>& force,
+                std::vector<double>& virial,
+                std::vector<double>& atom_energy,
+                std::vector<double>& atom_virial,
+                const std::vector<double>& coord,
+                const std::vector<int>& atype,
+                const std::vector<double>& box,
+                const std::vector<double>& fparam,
+                const std::vector<double>& aparam,
+                const std::vector<double>& charge_spin,
+                const bool atomic) override;
+  void computew(std::vector<double>& ener,
+                std::vector<float>& force,
+                std::vector<float>& virial,
+                std::vector<float>& atom_energy,
+                std::vector<float>& atom_virial,
+                const std::vector<float>& coord,
+                const std::vector<int>& atype,
+                const std::vector<float>& box,
+                const std::vector<float>& fparam,
+                const std::vector<float>& aparam,
+                const std::vector<double>& charge_spin,
+                const bool atomic) override;
+  void computew(std::vector<double>& ener,
+                std::vector<double>& force,
+                std::vector<double>& virial,
+                std::vector<double>& atom_energy,
+                std::vector<double>& atom_virial,
+                const std::vector<double>& coord,
+                const std::vector<int>& atype,
+                const std::vector<double>& box,
+                const int nghost,
+                const InputNlist& inlist,
+                const int& ago,
+                const std::vector<double>& fparam,
+                const std::vector<double>& aparam,
+                const std::vector<double>& charge_spin,
+                const bool atomic) override;
+  void computew(std::vector<double>& ener,
+                std::vector<float>& force,
+                std::vector<float>& virial,
+                std::vector<float>& atom_energy,
+                std::vector<float>& atom_virial,
+                const std::vector<float>& coord,
+                const std::vector<int>& atype,
+                const std::vector<float>& box,
+                const int nghost,
+                const InputNlist& inlist,
+                const int& ago,
+                const std::vector<float>& fparam,
+                const std::vector<float>& aparam,
+                const std::vector<double>& charge_spin,
+                const bool atomic) override;
+  void computew_mixed_type(std::vector<double>& ener,
+                           std::vector<double>& force,
+                           std::vector<double>& virial,
+                           std::vector<double>& atom_energy,
+                           std::vector<double>& atom_virial,
+                           const int& nframes,
+                           const std::vector<double>& coord,
+                           const std::vector<int>& atype,
+                           const std::vector<double>& box,
+                           const std::vector<double>& fparam,
+                           const std::vector<double>& aparam,
+                           const std::vector<double>& charge_spin,
+                           const bool atomic) override;
+  void computew_mixed_type(std::vector<double>& ener,
+                           std::vector<float>& force,
+                           std::vector<float>& virial,
+                           std::vector<float>& atom_energy,
+                           std::vector<float>& atom_virial,
+                           const int& nframes,
+                           const std::vector<float>& coord,
+                           const std::vector<int>& atype,
+                           const std::vector<float>& box,
+                           const std::vector<float>& fparam,
+                           const std::vector<float>& aparam,
+                           const std::vector<double>& charge_spin,
+                           const bool atomic) override;
+
  private:
   bool inited;
   int ntypes;
   int dfparam;
   int daparam;
-  int dim_chg_spin;
+  int dchgspin;
   bool aparam_nall;
   bool has_default_fparam_;
   std::vector<double> default_fparam_;
@@ -255,6 +344,7 @@ class DeepPotPTExpt : public DeepPotBackend {
                        const std::vector<VALUETYPE>& box,
                        const std::vector<VALUETYPE>& fparam,
                        const std::vector<VALUETYPE>& aparam,
+                       const std::vector<double>& charge_spin,
                        const bool atomic);
 
   /**
@@ -272,6 +362,7 @@ class DeepPotPTExpt : public DeepPotBackend {
                                const std::vector<VALUETYPE>& box,
                                const std::vector<VALUETYPE>& fparam,
                                const std::vector<VALUETYPE>& aparam,
+                               const std::vector<double>& charge_spin,
                                const bool atomic);
 
   /**
@@ -282,6 +373,7 @@ class DeepPotPTExpt : public DeepPotBackend {
    * @param[in] mapping Mapping tensor.
    * @param[in] fparam Frame parameter tensor (or empty).
    * @param[in] aparam Atomic parameter tensor (or empty).
+   * @param[in] charge_spin Charge/spin tensor (or empty).
    * @return Vector of output tensors in sorted key order.
    */
   std::vector<torch::Tensor> run_model(const torch::Tensor& coord,
@@ -289,13 +381,15 @@ class DeepPotPTExpt : public DeepPotBackend {
                                        const torch::Tensor& nlist,
                                        const torch::Tensor& mapping,
                                        const torch::Tensor& fparam,
-                                       const torch::Tensor& aparam);
+                                       const torch::Tensor& aparam,
+                                       const torch::Tensor& charge_spin);
 
   /**
    * @brief Run the with-comm .pt2 artifact with comm tensors appended.
    *
    * @param[in] base 4-6 base inputs (coord, atype, nlist, mapping,
    *            fparam?, aparam?) — same as ``run_model``.
+   * @param[in] charge_spin Charge/spin tensor (or empty).
    * @param[in] comm_tensors 8 comm tensors in canonical positional
    *            order: send_list, send_proc, recv_proc, send_num,
    *            recv_num, communicator, nlocal, nghost.
@@ -307,6 +401,7 @@ class DeepPotPTExpt : public DeepPotBackend {
       const torch::Tensor& mapping,
       const torch::Tensor& fparam,
       const torch::Tensor& aparam,
+      const torch::Tensor& charge_spin,
       const std::vector<at::Tensor>& comm_tensors);
 
   /**
diff --git a/source/api_cc/src/DeepPot.cc b/source/api_cc/src/DeepPot.cc
index ae543e1e31..30c64ab83f 100644
--- a/source/api_cc/src/DeepPot.cc
+++ b/source/api_cc/src/DeepPot.cc
@@ -93,11 +93,12 @@ void DeepPot::compute(ENERGYTYPE& dener,
                       const std::vector<int>& datype_,
                       const std::vector<VALUETYPE>& dbox,
                       const std::vector<VALUETYPE>& fparam_,
-                      const std::vector<VALUETYPE>& aparam_) {
+                      const std::vector<VALUETYPE>& aparam_,
+                      const std::vector<double>& charge_spin) {
   std::vector<ENERGYTYPE> dener_;
   std::vector<VALUETYPE> datom_energy_, datom_virial_;
   dp->computew(dener_, dforce_, dvirial, datom_energy_, datom_virial_, dcoord_,
-               datype_, dbox, fparam_, aparam_, false);
+               datype_, dbox, fparam_, aparam_, charge_spin, false);
   dener = dener_[0];
 }
 
@@ -109,10 +110,11 @@ void DeepPot::compute(std::vector<ENERGYTYPE>& dener,
                       const std::vector<int>& datype_,
                       const std::vector<VALUETYPE>& dbox,
                       const std::vector<VALUETYPE>& fparam_,
-                      const std::vector<VALUETYPE>& aparam_) {
+                      const std::vector<VALUETYPE>& aparam_,
+                      const std::vector<double>& charge_spin) {
   std::vector<VALUETYPE> datom_energy_, datom_virial_;
   dp->computew(dener, dforce_, dvirial, datom_energy_, datom_virial_, dcoord_,
-               datype_, dbox, fparam_, aparam_, false);
+               datype_, dbox, fparam_, aparam_, charge_spin, false);
 }
 
 template void DeepPot::compute<double>(ENERGYTYPE& dener,
@@ -122,7 +124,8 @@ template void DeepPot::compute<double>(ENERGYTYPE& dener,
                                        const std::vector<int>& datype_,
                                        const std::vector<double>& dbox,
                                        const std::vector<double>& fparam,
-                                       const std::vector<double>& aparam);
+                                       const std::vector<double>& aparam,
+                                       const std::vector<double>& charge_spin);
 
 template void DeepPot::compute<float>(ENERGYTYPE& dener,
                                       std::vector<float>& dforce_,
@@ -131,7 +134,8 @@ template void DeepPot::compute<float>(ENERGYTYPE& dener,
                                       const std::vector<int>& datype_,
                                       const std::vector<float>& dbox,
                                       const std::vector<float>& fparam,
-                                      const std::vector<float>& aparam);
+                                      const std::vector<float>& aparam,
+                                      const std::vector<double>& charge_spin);
 
 template void DeepPot::compute<double>(std::vector<ENERGYTYPE>& dener,
                                        std::vector<double>& dforce_,
@@ -140,7 +144,8 @@ template void DeepPot::compute<double>(std::vector<ENERGYTYPE>& dener,
                                        const std::vector<int>& datype_,
                                        const std::vector<double>& dbox,
                                        const std::vector<double>& fparam,
-                                       const std::vector<double>& aparam);
+                                       const std::vector<double>& aparam,
+                                       const std::vector<double>& charge_spin);
 
 template void DeepPot::compute<float>(std::vector<ENERGYTYPE>& dener,
                                       std::vector<float>& dforce_,
@@ -149,7 +154,8 @@ template void DeepPot::compute<float>(std::vector<ENERGYTYPE>& dener,
                                       const std::vector<int>& datype_,
                                       const std::vector<float>& dbox,
                                       const std::vector<float>& fparam,
-                                      const std::vector<float>& aparam);
+                                      const std::vector<float>& aparam,
+                                      const std::vector<double>& charge_spin);
 
 template <typename VALUETYPE>
 void DeepPot::compute(ENERGYTYPE& dener,
@@ -162,11 +168,13 @@ void DeepPot::compute(ENERGYTYPE& dener,
                       const InputNlist& lmp_list,
                       const int& ago,
                       const std::vector<VALUETYPE>& fparam_,
-                      const std::vector<VALUETYPE>& aparam__) {
+                      const std::vector<VALUETYPE>& aparam__,
+                      const std::vector<double>& charge_spin) {
   std::vector<ENERGYTYPE> dener_;
   std::vector<VALUETYPE> datom_energy_, datom_virial_;
   dp->computew(dener_, dforce_, dvirial, datom_energy_, datom_virial_, dcoord_,
-               datype_, dbox, nghost, lmp_list, ago, fparam_, aparam__, false);
+               datype_, dbox, nghost, lmp_list, ago, fparam_, aparam__,
+               charge_spin, false);
   dener = dener_[0];
 }
 
@@ -181,10 +189,12 @@ void DeepPot::compute(std::vector<ENERGYTYPE>& dener,
                       const InputNlist& lmp_list,
                       const int& ago,
                       const std::vector<VALUETYPE>& fparam_,
-                      const std::vector<VALUETYPE>& aparam__) {
+                      const std::vector<VALUETYPE>& aparam__,
+                      const std::vector<double>& charge_spin) {
   std::vector<VALUETYPE> datom_energy_, datom_virial_;
   dp->computew(dener, dforce_, dvirial, datom_energy_, datom_virial_, dcoord_,
-               datype_, dbox, nghost, lmp_list, ago, fparam_, aparam__, false);
+               datype_, dbox, nghost, lmp_list, ago, fparam_, aparam__,
+               charge_spin, false);
 }
 
 template void DeepPot::compute<double>(ENERGYTYPE& dener,
@@ -197,7 +207,8 @@ template void DeepPot::compute<double>(ENERGYTYPE& dener,
                                        const InputNlist& lmp_list,
                                        const int& ago,
                                        const std::vector<double>& fparam,
-                                       const std::vector<double>& aparam_);
+                                       const std::vector<double>& aparam_,
+                                       const std::vector<double>& charge_spin);
 
 template void DeepPot::compute<float>(ENERGYTYPE& dener,
                                       std::vector<float>& dforce_,
@@ -209,7 +220,8 @@ template void DeepPot::compute<float>(ENERGYTYPE& dener,
                                       const InputNlist& lmp_list,
                                       const int& ago,
                                       const std::vector<float>& fparam,
-                                      const std::vector<float>& aparam_);
+                                      const std::vector<float>& aparam_,
+                                      const std::vector<double>& charge_spin);
 
 template void DeepPot::compute<double>(std::vector<ENERGYTYPE>& dener,
                                        std::vector<double>& dforce_,
@@ -221,7 +233,8 @@ template void DeepPot::compute<double>(std::vector<ENERGYTYPE>& dener,
                                        const InputNlist& lmp_list,
                                        const int& ago,
                                        const std::vector<double>& fparam,
-                                       const std::vector<double>& aparam_);
+                                       const std::vector<double>& aparam_,
+                                       const std::vector<double>& charge_spin);
 
 template void DeepPot::compute<float>(std::vector<ENERGYTYPE>& dener,
                                       std::vector<float>& dforce_,
@@ -233,7 +246,8 @@ template void DeepPot::compute<float>(std::vector<ENERGYTYPE>& dener,
                                       const InputNlist& lmp_list,
                                       const int& ago,
                                       const std::vector<float>& fparam,
-                                      const std::vector<float>& aparam_);
+                                      const std::vector<float>& aparam_,
+                                      const std::vector<double>& charge_spin);
 
 template <typename VALUETYPE>
 void DeepPot::compute(ENERGYTYPE& dener,
@@ -245,10 +259,11 @@ void DeepPot::compute(ENERGYTYPE& dener,
                       const std::vector<int>& datype_,
                       const std::vector<VALUETYPE>& dbox,
                       const std::vector<VALUETYPE>& fparam_,
-                      const std::vector<VALUETYPE>& aparam_) {
+                      const std::vector<VALUETYPE>& aparam_,
+                      const std::vector<double>& charge_spin) {
   std::vector<ENERGYTYPE> dener_;
   dp->computew(dener_, dforce_, dvirial, datom_energy_, datom_virial_, dcoord_,
-               datype_, dbox, fparam_, aparam_, true);
+               datype_, dbox, fparam_, aparam_, charge_spin, true);
   dener = dener_[0];
 }
 template <typename VALUETYPE>
@@ -261,9 +276,10 @@ void DeepPot::compute(std::vector<ENERGYTYPE>& dener,
                       const std::vector<int>& datype_,
                       const std::vector<VALUETYPE>& dbox,
                       const std::vector<VALUETYPE>& fparam_,
-                      const std::vector<VALUETYPE>& aparam_) {
+                      const std::vector<VALUETYPE>& aparam_,
+                      const std::vector<double>& charge_spin) {
   dp->computew(dener, dforce_, dvirial, datom_energy_, datom_virial_, dcoord_,
-               datype_, dbox, fparam_, aparam_, true);
+               datype_, dbox, fparam_, aparam_, charge_spin, true);
 }
 
 template void DeepPot::compute<double>(ENERGYTYPE& dener,
@@ -275,7 +291,8 @@ template void DeepPot::compute<double>(ENERGYTYPE& dener,
                                        const std::vector<int>& datype_,
                                        const std::vector<double>& dbox,
                                        const std::vector<double>& fparam,
-                                       const std::vector<double>& aparam);
+                                       const std::vector<double>& aparam,
+                                       const std::vector<double>& charge_spin);
 
 template void DeepPot::compute<float>(ENERGYTYPE& dener,
                                       std::vector<float>& dforce_,
@@ -286,7 +303,8 @@ template void DeepPot::compute<float>(ENERGYTYPE& dener,
                                       const std::vector<int>& datype_,
                                       const std::vector<float>& dbox,
                                       const std::vector<float>& fparam,
-                                      const std::vector<float>& aparam);
+                                      const std::vector<float>& aparam,
+                                      const std::vector<double>& charge_spin);
 
 template void DeepPot::compute<double>(std::vector<ENERGYTYPE>& dener,
                                        std::vector<double>& dforce_,
@@ -297,7 +315,8 @@ template void DeepPot::compute<double>(std::vector<ENERGYTYPE>& dener,
                                        const std::vector<int>& datype_,
                                        const std::vector<double>& dbox,
                                        const std::vector<double>& fparam,
-                                       const std::vector<double>& aparam);
+                                       const std::vector<double>& aparam,
+                                       const std::vector<double>& charge_spin);
 
 template void DeepPot::compute<float>(std::vector<ENERGYTYPE>& dener,
                                       std::vector<float>& dforce_,
@@ -308,7 +327,8 @@ template void DeepPot::compute<float>(std::vector<ENERGYTYPE>& dener,
                                       const std::vector<int>& datype_,
                                       const std::vector<float>& dbox,
                                       const std::vector<float>& fparam,
-                                      const std::vector<float>& aparam);
+                                      const std::vector<float>& aparam,
+                                      const std::vector<double>& charge_spin);
 
 template <typename VALUETYPE>
 void DeepPot::compute(ENERGYTYPE& dener,
@@ -323,10 +343,12 @@ void DeepPot::compute(ENERGYTYPE& dener,
                       const InputNlist& lmp_list,
                       const int& ago,
                       const std::vector<VALUETYPE>& fparam_,
-                      const std::vector<VALUETYPE>& aparam__) {
+                      const std::vector<VALUETYPE>& aparam__,
+                      const std::vector<double>& charge_spin) {
   std::vector<ENERGYTYPE> dener_;
   dp->computew(dener_, dforce_, dvirial, datom_energy_, datom_virial_, dcoord_,
-               datype_, dbox, nghost, lmp_list, ago, fparam_, aparam__, true);
+               datype_, dbox, nghost, lmp_list, ago, fparam_, aparam__,
+               charge_spin, true);
   dener = dener_[0];
 }
 template <typename VALUETYPE>
@@ -342,9 +364,11 @@ void DeepPot::compute(std::vector<ENERGYTYPE>& dener,
                       const InputNlist& lmp_list,
                       const int& ago,
                       const std::vector<VALUETYPE>& fparam_,
-                      const std::vector<VALUETYPE>& aparam__) {
+                      const std::vector<VALUETYPE>& aparam__,
+                      const std::vector<double>& charge_spin) {
   dp->computew(dener, dforce_, dvirial, datom_energy_, datom_virial_, dcoord_,
-               datype_, dbox, nghost, lmp_list, ago, fparam_, aparam__, true);
+               datype_, dbox, nghost, lmp_list, ago, fparam_, aparam__,
+               charge_spin, true);
 }
 
 template void DeepPot::compute<double>(ENERGYTYPE& dener,
@@ -359,7 +383,8 @@ template void DeepPot::compute<double>(ENERGYTYPE& dener,
                                        const InputNlist& lmp_list,
                                        const int& ago,
                                        const std::vector<double>& fparam,
-                                       const std::vector<double>& aparam_);
+                                       const std::vector<double>& aparam_,
+                                       const std::vector<double>& charge_spin);
 
 template void DeepPot::compute<float>(ENERGYTYPE& dener,
                                       std::vector<float>& dforce_,
@@ -373,7 +398,8 @@ template void DeepPot::compute<float>(ENERGYTYPE& dener,
                                       const InputNlist& lmp_list,
                                       const int& ago,
                                       const std::vector<float>& fparam,
-                                      const std::vector<float>& aparam_);
+                                      const std::vector<float>& aparam_,
+                                      const std::vector<double>& charge_spin);
 
 template void DeepPot::compute<double>(std::vector<ENERGYTYPE>& dener,
                                        std::vector<double>& dforce_,
@@ -387,7 +413,8 @@ template void DeepPot::compute<double>(std::vector<ENERGYTYPE>& dener,
                                        const InputNlist& lmp_list,
                                        const int& ago,
                                        const std::vector<double>& fparam,
-                                       const std::vector<double>& aparam_);
+                                       const std::vector<double>& aparam_,
+                                       const std::vector<double>& charge_spin);
 
 template void DeepPot::compute<float>(std::vector<ENERGYTYPE>& dener,
                                       std::vector<float>& dforce_,
@@ -401,7 +428,8 @@ template void DeepPot::compute<float>(std::vector<ENERGYTYPE>& dener,
                                       const InputNlist& lmp_list,
                                       const int& ago,
                                       const std::vector<float>& fparam,
-                                      const std::vector<float>& aparam_);
+                                      const std::vector<float>& aparam_,
+                                      const std::vector<double>& charge_spin);
 
 // mixed type
 template <typename VALUETYPE>
@@ -413,12 +441,13 @@ void DeepPot::compute_mixed_type(ENERGYTYPE& dener,
                                  const std::vector<int>& datype_,
                                  const std::vector<VALUETYPE>& dbox,
                                  const std::vector<VALUETYPE>& fparam_,
-                                 const std::vector<VALUETYPE>& aparam_) {
+                                 const std::vector<VALUETYPE>& aparam_,
+                                 const std::vector<double>& charge_spin) {
   std::vector<ENERGYTYPE> dener_;
   std::vector<VALUETYPE> datom_energy_, datom_virial_;
   dp->computew_mixed_type(dener_, dforce_, dvirial, datom_energy_,
                           datom_virial_, nframes, dcoord_, datype_, dbox,
-                          fparam_, aparam_, false);
+                          fparam_, aparam_, charge_spin, false);
   dener = dener_[0];
 }
 template <typename VALUETYPE>
@@ -430,11 +459,12 @@ void DeepPot::compute_mixed_type(std::vector<ENERGYTYPE>& dener,
                                  const std::vector<int>& datype_,
                                  const std::vector<VALUETYPE>& dbox,
                                  const std::vector<VALUETYPE>& fparam_,
-                                 const std::vector<VALUETYPE>& aparam_) {
+                                 const std::vector<VALUETYPE>& aparam_,
+                                 const std::vector<double>& charge_spin) {
   std::vector<VALUETYPE> datom_energy_, datom_virial_;
   dp->computew_mixed_type(dener, dforce_, dvirial, datom_energy_, datom_virial_,
                           nframes, dcoord_, datype_, dbox, fparam_, aparam_,
-                          false);
+                          charge_spin, false);
 }
 
 template void DeepPot::compute_mixed_type<double>(
@@ -446,7 +476,8 @@ template void DeepPot::compute_mixed_type<double>(
     const std::vector<int>& datype_,
     const std::vector<double>& dbox,
     const std::vector<double>& fparam,
-    const std::vector<double>& aparam);
+    const std::vector<double>& aparam,
+    const std::vector<double>& charge_spin);
 
 template void DeepPot::compute_mixed_type<float>(
     ENERGYTYPE& dener,
@@ -457,7 +488,8 @@ template void DeepPot::compute_mixed_type<float>(
     const std::vector<int>& datype_,
     const std::vector<float>& dbox,
     const std::vector<float>& fparam,
-    const std::vector<float>& aparam);
+    const std::vector<float>& aparam,
+    const std::vector<double>& charge_spin);
 
 template void DeepPot::compute_mixed_type<double>(
     std::vector<ENERGYTYPE>& dener,
@@ -468,7 +500,8 @@ template void DeepPot::compute_mixed_type<double>(
     const std::vector<int>& datype_,
     const std::vector<double>& dbox,
     const std::vector<double>& fparam,
-    const std::vector<double>& aparam);
+    const std::vector<double>& aparam,
+    const std::vector<double>& charge_spin);
 
 template void DeepPot::compute_mixed_type<float>(
     std::vector<ENERGYTYPE>& dener,
@@ -479,7 +512,8 @@ template void DeepPot::compute_mixed_type<float>(
     const std::vector<int>& datype_,
     const std::vector<float>& dbox,
     const std::vector<float>& fparam,
-    const std::vector<float>& aparam);
+    const std::vector<float>& aparam,
+    const std::vector<double>& charge_spin);
 
 template <typename VALUETYPE>
 void DeepPot::compute_mixed_type(ENERGYTYPE& dener,
@@ -492,11 +526,12 @@ void DeepPot::compute_mixed_type(ENERGYTYPE& dener,
                                  const std::vector<int>& datype_,
                                  const std::vector<VALUETYPE>& dbox,
                                  const std::vector<VALUETYPE>& fparam_,
-                                 const std::vector<VALUETYPE>& aparam_) {
+                                 const std::vector<VALUETYPE>& aparam_,
+                                 const std::vector<double>& charge_spin) {
   std::vector<ENERGYTYPE> dener_;
   dp->computew_mixed_type(dener_, dforce_, dvirial, datom_energy_,
                           datom_virial_, nframes, dcoord_, datype_, dbox,
-                          fparam_, aparam_, true);
+                          fparam_, aparam_, charge_spin, true);
   dener = dener_[0];
 }
 template <typename VALUETYPE>
@@ -510,10 +545,11 @@ void DeepPot::compute_mixed_type(std::vector<ENERGYTYPE>& dener,
                                  const std::vector<int>& datype_,
                                  const std::vector<VALUETYPE>& dbox,
                                  const std::vector<VALUETYPE>& fparam_,
-                                 const std::vector<VALUETYPE>& aparam_) {
+                                 const std::vector<VALUETYPE>& aparam_,
+                                 const std::vector<double>& charge_spin) {
   dp->computew_mixed_type(dener, dforce_, dvirial, datom_energy_, datom_virial_,
                           nframes, dcoord_, datype_, dbox, fparam_, aparam_,
-                          true);
+                          charge_spin, true);
 }
 
 template void DeepPot::compute_mixed_type<double>(
@@ -527,7 +563,8 @@ template void DeepPot::compute_mixed_type<double>(
     const std::vector<int>& datype_,
     const std::vector<double>& dbox,
     const std::vector<double>& fparam,
-    const std::vector<double>& aparam);
+    const std::vector<double>& aparam,
+    const std::vector<double>& charge_spin);
 
 template void DeepPot::compute_mixed_type<float>(
     ENERGYTYPE& dener,
@@ -540,7 +577,8 @@ template void DeepPot::compute_mixed_type<float>(
     const std::vector<int>& datype_,
     const std::vector<float>& dbox,
     const std::vector<float>& fparam,
-    const std::vector<float>& aparam);
+    const std::vector<float>& aparam,
+    const std::vector<double>& charge_spin);
 
 template void DeepPot::compute_mixed_type<double>(
     std::vector<ENERGYTYPE>& dener,
@@ -553,7 +591,8 @@ template void DeepPot::compute_mixed_type<double>(
     const std::vector<int>& datype_,
     const std::vector<double>& dbox,
     const std::vector<double>& fparam,
-    const std::vector<double>& aparam);
+    const std::vector<double>& aparam,
+    const std::vector<double>& charge_spin);
 
 template void DeepPot::compute_mixed_type<float>(
     std::vector<ENERGYTYPE>& dener,
@@ -566,7 +605,10 @@ template void DeepPot::compute_mixed_type<float>(
     const std::vector<int>& datype_,
     const std::vector<float>& dbox,
     const std::vector<float>& fparam,
-    const std::vector<float>& aparam);
+    const std::vector<float>& aparam,
+    const std::vector<double>& charge_spin);
+
+int DeepPot::dim_chg_spin() const { return dp->dim_chg_spin(); }
 
 DeepPotModelDevi::DeepPotModelDevi() {
   inited = false;
@@ -616,7 +658,8 @@ void DeepPotModelDevi::compute(std::vector<ENERGYTYPE>& all_energy,
                                const std::vector<int>& datype_,
                                const std::vector<VALUETYPE>& dbox,
                                const std::vector<VALUETYPE>& fparam,
-                               const std::vector<VALUETYPE>& aparam_) {
+                               const std::vector<VALUETYPE>& aparam_,
+                               const std::vector<double>& charge_spin) {
   // without nlist
   if (numb_models == 0) {
     return;
@@ -626,7 +669,7 @@ void DeepPotModelDevi::compute(std::vector<ENERGYTYPE>& all_energy,
   all_virial.resize(numb_models);
   for (unsigned ii = 0; ii < numb_models; ++ii) {
     dps[ii]->compute(all_energy[ii], all_force[ii], all_virial[ii], dcoord_,
-                     datype_, dbox, fparam, aparam_);
+                     datype_, dbox, fparam, aparam_, charge_spin);
   }
 }
 
@@ -638,7 +681,8 @@ template void DeepPotModelDevi::compute<double>(
     const std::vector<int>& datype_,
     const std::vector<double>& dbox,
     const std::vector<double>& fparam,
-    const std::vector<double>& aparam);
+    const std::vector<double>& aparam,
+    const std::vector<double>& charge_spin);
 
 template void DeepPotModelDevi::compute<float>(
     std::vector<ENERGYTYPE>& all_energy,
@@ -648,7 +692,8 @@ template void DeepPotModelDevi::compute<float>(
     const std::vector<int>& datype_,
     const std::vector<float>& dbox,
     const std::vector<float>& fparam,
-    const std::vector<float>& aparam);
+    const std::vector<float>& aparam,
+    const std::vector<double>& charge_spin);
 
 template <typename VALUETYPE>
 void DeepPotModelDevi::compute(
@@ -661,7 +706,8 @@ void DeepPotModelDevi::compute(
     const std::vector<int>& datype_,
     const std::vector<VALUETYPE>& dbox,
     const std::vector<VALUETYPE>& fparam,
-    const std::vector<VALUETYPE>& aparam_) {
+    const std::vector<VALUETYPE>& aparam_,
+    const std::vector<double>& charge_spin) {
   if (numb_models == 0) {
     return;
   }
@@ -673,7 +719,7 @@ void DeepPotModelDevi::compute(
   for (unsigned ii = 0; ii < numb_models; ++ii) {
     dps[ii]->compute(all_energy[ii], all_force[ii], all_virial[ii],
                      all_atom_energy[ii], all_atom_virial[ii], dcoord_, datype_,
-                     dbox, fparam, aparam_);
+                     dbox, fparam, aparam_, charge_spin);
   }
 }
 
@@ -687,7 +733,8 @@ template void DeepPotModelDevi::compute<double>(
     const std::vector<int>& datype_,
     const std::vector<double>& dbox,
     const std::vector<double>& fparam,
-    const std::vector<double>& aparam);
+    const std::vector<double>& aparam,
+    const std::vector<double>& charge_spin);
 
 template void DeepPotModelDevi::compute<float>(
     std::vector<ENERGYTYPE>& all_energy,
@@ -699,7 +746,8 @@ template void DeepPotModelDevi::compute<float>(
     const std::vector<int>& datype_,
     const std::vector<float>& dbox,
     const std::vector<float>& fparam,
-    const std::vector<float>& aparam);
+    const std::vector<float>& aparam,
+    const std::vector<double>& charge_spin);
 
 template <typename VALUETYPE>
 void DeepPotModelDevi::compute(std::vector<ENERGYTYPE>& all_energy,
@@ -712,7 +760,8 @@ void DeepPotModelDevi::compute(std::vector<ENERGYTYPE>& all_energy,
                                const InputNlist& lmp_list,
                                const int& ago,
                                const std::vector<VALUETYPE>& fparam,
-                               const std::vector<VALUETYPE>& aparam_) {
+                               const std::vector<VALUETYPE>& aparam_,
+                               const std::vector<double>& charge_spin) {
   if (numb_models == 0) {
     return;
   }
@@ -721,7 +770,8 @@ void DeepPotModelDevi::compute(std::vector<ENERGYTYPE>& all_energy,
   all_virial.resize(numb_models);
   for (unsigned ii = 0; ii < numb_models; ++ii) {
     dps[ii]->compute(all_energy[ii], all_force[ii], all_virial[ii], dcoord_,
-                     datype_, dbox, nghost, lmp_list, ago, fparam, aparam_);
+                     datype_, dbox, nghost, lmp_list, ago, fparam, aparam_,
+                     charge_spin);
   }
 }
 
@@ -736,7 +786,8 @@ template void DeepPotModelDevi::compute<double>(
     const InputNlist& lmp_list,
     const int& ago,
     const std::vector<double>& fparam,
-    const std::vector<double>& aparam);
+    const std::vector<double>& aparam,
+    const std::vector<double>& charge_spin);
 
 template void DeepPotModelDevi::compute<float>(
     std::vector<ENERGYTYPE>& all_energy,
@@ -749,7 +800,8 @@ template void DeepPotModelDevi::compute<float>(
     const InputNlist& lmp_list,
     const int& ago,
     const std::vector<float>& fparam,
-    const std::vector<float>& aparam);
+    const std::vector<float>& aparam,
+    const std::vector<double>& charge_spin);
 
 template <typename VALUETYPE>
 void DeepPotModelDevi::compute(
@@ -765,7 +817,8 @@ void DeepPotModelDevi::compute(
     const InputNlist& lmp_list,
     const int& ago,
     const std::vector<VALUETYPE>& fparam,
-    const std::vector<VALUETYPE>& aparam_) {
+    const std::vector<VALUETYPE>& aparam_,
+    const std::vector<double>& charge_spin) {
   if (numb_models == 0) {
     return;
   }
@@ -777,7 +830,7 @@ void DeepPotModelDevi::compute(
   for (unsigned ii = 0; ii < numb_models; ++ii) {
     dps[ii]->compute(all_energy[ii], all_force[ii], all_virial[ii],
                      all_atom_energy[ii], all_atom_virial[ii], dcoord_, datype_,
-                     dbox, nghost, lmp_list, ago, fparam, aparam_);
+                     dbox, nghost, lmp_list, ago, fparam, aparam_, charge_spin);
   }
 }
 
@@ -794,7 +847,8 @@ template void DeepPotModelDevi::compute<double>(
     const InputNlist& lmp_list,
     const int& ago,
     const std::vector<double>& fparam,
-    const std::vector<double>& aparam);
+    const std::vector<double>& aparam,
+    const std::vector<double>& charge_spin);
 
 template void DeepPotModelDevi::compute<float>(
     std::vector<ENERGYTYPE>& all_energy,
@@ -809,4 +863,5 @@ template void DeepPotModelDevi::compute<float>(
     const InputNlist& lmp_list,
     const int& ago,
     const std::vector<float>& fparam,
-    const std::vector<float>& aparam);
+    const std::vector<float>& aparam,
+    const std::vector<double>& charge_spin);
diff --git a/source/api_cc/src/DeepPotPT.cc b/source/api_cc/src/DeepPotPT.cc
index d69dbb8f82..d29783089c 100644
--- a/source/api_cc/src/DeepPotPT.cc
+++ b/source/api_cc/src/DeepPotPT.cc
@@ -135,6 +135,24 @@ void DeepPotPT::init(const std::string& model,
   } else {
     has_default_fparam_ = false;
   }
+  // Charge/spin embedding (e.g. DPA3 add_chg_spin_ebd).  Guarded with
+  // find_method so .pth models exported before charge_spin support
+  // (which lack these jit-exported methods) keep loading with dchgspin=0.
+  if (module.find_method("get_dim_chg_spin")) {
+    dchgspin = module.run_method("get_dim_chg_spin").toInt();
+  } else {
+    dchgspin = 0;
+  }
+  default_chg_spin_.clear();
+  if (dchgspin > 0 && module.find_method("get_default_chg_spin")) {
+    auto cs = module.run_method("get_default_chg_spin");
+    if (!cs.isNone()) {
+      torch::Tensor cs_t =
+          cs.toTensor().to(torch::kFloat64).to(torch::kCPU).view({-1});
+      default_chg_spin_.assign(cs_t.data_ptr<double>(),
+                               cs_t.data_ptr<double>() + cs_t.numel());
+    }
+  }
   inited = true;
 }
 
@@ -163,6 +181,7 @@ void DeepPotPT::compute(ENERGYVTYPE& ener,
                         const int& ago,
                         const std::vector<VALUETYPE>& fparam,
                         const std::vector<VALUETYPE>& aparam,
+                        const std::vector<double>& charge_spin,
                         const bool atomic) {
   torch::Device device(torch::kCUDA, gpu_id);
   if (!gpu_enabled) {
@@ -242,18 +261,76 @@ void DeepPotPT::compute(ENERGYVTYPE& ener,
             options)
             .to(device);
   }
-  auto outputs =
-      (do_message_passing)
-          ? module
-                .run_method("forward_lower", coord_wrapped_Tensor, atype_Tensor,
-                            firstneigh_tensor, mapping_tensor, fparam_tensor,
-                            aparam_tensor, do_atom_virial_tensor, comm_dict)
-                .toGenericDict()
-          : module
-                .run_method("forward_lower", coord_wrapped_Tensor, atype_Tensor,
-                            firstneigh_tensor, mapping_tensor, fparam_tensor,
-                            aparam_tensor, do_atom_virial_tensor)
-                .toGenericDict();
+  // Build charge_spin tensor (always float64): use the runtime value when
+  // provided, otherwise fall back to the model's stored default_chg_spin.
+  // Only threaded into forward_lower when the model has a charge/spin
+  // embedding (dchgspin > 0), so .pth models without it are unaffected.
+  c10::optional<torch::Tensor> charge_spin_tensor;
+  if (dchgspin > 0) {
+    auto dbl_options = torch::TensorOptions().dtype(torch::kFloat64);
+    if (!charge_spin.empty()) {
+      // Single-frame path: charge_spin must hold exactly dim_chg_spin values.
+      if (static_cast<int>(charge_spin.size()) != dchgspin) {
+        throw deepmd::deepmd_exception(
+            "charge_spin has " + std::to_string(charge_spin.size()) +
+            " values but the model expects dim_chg_spin=" +
+            std::to_string(dchgspin) + ".");
+      }
+      charge_spin_tensor =
+          torch::from_blob(const_cast<double*>(charge_spin.data()),
+                           {1, static_cast<std::int64_t>(charge_spin.size())},
+                           dbl_options)
+              .clone()
+              .to(device);
+    } else if (!default_chg_spin_.empty()) {
+      if (static_cast<int>(default_chg_spin_.size()) != dchgspin) {
+        throw deepmd::deepmd_exception(
+            "default_chg_spin has " + std::to_string(default_chg_spin_.size()) +
+            " values but the model expects dim_chg_spin=" +
+            std::to_string(dchgspin) + ".");
+      }
+      charge_spin_tensor =
+          torch::from_blob(const_cast<double*>(default_chg_spin_.data()),
+                           {1, dchgspin}, dbl_options)
+              .clone()
+              .to(device);
+    } else {
+      throw deepmd::deepmd_exception(
+          "charge_spin is empty and no default_chg_spin is available in the "
+          "model. Provide charge_spin explicitly or regenerate the model with "
+          "a default charge/spin value.");
+    }
+  }
+  c10::IValue outputs_ival;
+  if (dchgspin > 0) {
+    // charge_spin model. DPA3 (the only charge/spin descriptor) always uses
+    // message passing, so comm_dict is populated; for the non-message-passing
+    // edge case pass a real None for comm_dict (an empty Dict would wrongly
+    // flip the descriptor into parallel mode).
+    if (do_message_passing) {
+      outputs_ival = module.run_method(
+          "forward_lower", coord_wrapped_Tensor, atype_Tensor,
+          firstneigh_tensor, mapping_tensor, fparam_tensor, aparam_tensor,
+          do_atom_virial_tensor, comm_dict, charge_spin_tensor);
+    } else {
+      outputs_ival = module.run_method(
+          "forward_lower", coord_wrapped_Tensor, atype_Tensor,
+          firstneigh_tensor, mapping_tensor, fparam_tensor, aparam_tensor,
+          do_atom_virial_tensor, c10::IValue(), charge_spin_tensor);
+    }
+  } else {
+    outputs_ival =
+        (do_message_passing)
+            ? module.run_method("forward_lower", coord_wrapped_Tensor,
+                                atype_Tensor, firstneigh_tensor, mapping_tensor,
+                                fparam_tensor, aparam_tensor,
+                                do_atom_virial_tensor, comm_dict)
+            : module.run_method("forward_lower", coord_wrapped_Tensor,
+                                atype_Tensor, firstneigh_tensor, mapping_tensor,
+                                fparam_tensor, aparam_tensor,
+                                do_atom_virial_tensor);
+  }
+  auto outputs = outputs_ival.toGenericDict();
   c10::IValue energy_ = outputs.at("energy");
   c10::IValue force_ = outputs.at("extended_force");
   c10::IValue virial_ = outputs.at("virial");
@@ -313,6 +390,7 @@ template void DeepPotPT::compute<double, std::vector<ENERGYTYPE>>(
     const int& ago,
     const std::vector<double>& fparam,
     const std::vector<double>& aparam,
+    const std::vector<double>& charge_spin,
     const bool atomic);
 template void DeepPotPT::compute<float, std::vector<ENERGYTYPE>>(
     std::vector<ENERGYTYPE>& ener,
@@ -328,6 +406,7 @@ template void DeepPotPT::compute<float, std::vector<ENERGYTYPE>>(
     const int& ago,
     const std::vector<float>& fparam,
     const std::vector<float>& aparam,
+    const std::vector<double>& charge_spin,
     const bool atomic);
 template <typename VALUETYPE, typename ENERGYVTYPE>
 void DeepPotPT::compute(ENERGYVTYPE& ener,
@@ -340,6 +419,7 @@ void DeepPotPT::compute(ENERGYVTYPE& ener,
                         const std::vector<VALUETYPE>& box,
                         const std::vector<VALUETYPE>& fparam,
                         const std::vector<VALUETYPE>& aparam,
+                        const std::vector<double>& charge_spin,
                         const bool atomic) {
   torch::Device device(torch::kCUDA, gpu_id);
   if (!gpu_enabled) {
@@ -391,6 +471,47 @@ void DeepPotPT::compute(ENERGYVTYPE& ener,
   inputs.push_back(aparam_tensor);
   bool do_atom_virial_tensor = atomic;
   inputs.push_back(do_atom_virial_tensor);
+  // Append charge_spin (always float64) only for models with a charge/spin
+  // embedding, so models without it (whose forward() lacks the parameter)
+  // keep working. Uses the runtime value when provided, else the model's
+  // stored default_chg_spin.
+  if (dchgspin > 0) {
+    auto dbl_options = torch::TensorOptions().dtype(torch::kFloat64);
+    c10::optional<torch::Tensor> charge_spin_tensor;
+    if (!charge_spin.empty()) {
+      // Single-frame path: charge_spin must hold exactly dim_chg_spin values.
+      if (static_cast<int>(charge_spin.size()) != dchgspin) {
+        throw deepmd::deepmd_exception(
+            "charge_spin has " + std::to_string(charge_spin.size()) +
+            " values but the model expects dim_chg_spin=" +
+            std::to_string(dchgspin) + ".");
+      }
+      charge_spin_tensor =
+          torch::from_blob(const_cast<double*>(charge_spin.data()),
+                           {1, static_cast<std::int64_t>(charge_spin.size())},
+                           dbl_options)
+              .clone()
+              .to(device);
+    } else if (!default_chg_spin_.empty()) {
+      if (static_cast<int>(default_chg_spin_.size()) != dchgspin) {
+        throw deepmd::deepmd_exception(
+            "default_chg_spin has " + std::to_string(default_chg_spin_.size()) +
+            " values but the model expects dim_chg_spin=" +
+            std::to_string(dchgspin) + ".");
+      }
+      charge_spin_tensor =
+          torch::from_blob(const_cast<double*>(default_chg_spin_.data()),
+                           {1, dchgspin}, dbl_options)
+              .clone()
+              .to(device);
+    } else {
+      throw deepmd::deepmd_exception(
+          "charge_spin is empty and no default_chg_spin is available in the "
+          "model. Provide charge_spin explicitly or regenerate the model with "
+          "a default charge/spin value.");
+    }
+    inputs.push_back(charge_spin_tensor);
+  }
   c10::Dict<c10::IValue, c10::IValue> outputs =
       module.forward(inputs).toGenericDict();
   c10::IValue energy_ = outputs.at("energy");
@@ -437,6 +558,7 @@ template void DeepPotPT::compute<double, std::vector<ENERGYTYPE>>(
     const std::vector<double>& box,
     const std::vector<double>& fparam,
     const std::vector<double>& aparam,
+    const std::vector<double>& charge_spin,
     const bool atomic);
 template void DeepPotPT::compute<float, std::vector<ENERGYTYPE>>(
     std::vector<ENERGYTYPE>& ener,
@@ -449,6 +571,7 @@ template void DeepPotPT::compute<float, std::vector<ENERGYTYPE>>(
     const std::vector<float>& box,
     const std::vector<float>& fparam,
     const std::vector<float>& aparam,
+    const std::vector<double>& charge_spin,
     const bool atomic);
 void DeepPotPT::get_type_map(std::string& type_map) {
   auto ret = module.run_method("get_type_map").toList();
@@ -472,7 +595,79 @@ void DeepPotPT::computew(std::vector<double>& ener,
                          const bool atomic) {
   translate_error([&] {
     compute(ener, force, virial, atom_energy, atom_virial, coord, atype, box,
-            fparam, aparam, atomic);
+            fparam, aparam, {}, atomic);
+  });
+}
+void DeepPotPT::computew(std::vector<double>& ener,
+                         std::vector<float>& force,
+                         std::vector<float>& virial,
+                         std::vector<float>& atom_energy,
+                         std::vector<float>& atom_virial,
+                         const std::vector<float>& coord,
+                         const std::vector<int>& atype,
+                         const std::vector<float>& box,
+                         const std::vector<float>& fparam,
+                         const std::vector<float>& aparam,
+                         const bool atomic) {
+  translate_error([&] {
+    compute(ener, force, virial, atom_energy, atom_virial, coord, atype, box,
+            fparam, aparam, {}, atomic);
+  });
+}
+void DeepPotPT::computew(std::vector<double>& ener,
+                         std::vector<double>& force,
+                         std::vector<double>& virial,
+                         std::vector<double>& atom_energy,
+                         std::vector<double>& atom_virial,
+                         const std::vector<double>& coord,
+                         const std::vector<int>& atype,
+                         const std::vector<double>& box,
+                         const int nghost,
+                         const InputNlist& inlist,
+                         const int& ago,
+                         const std::vector<double>& fparam,
+                         const std::vector<double>& aparam,
+                         const bool atomic) {
+  translate_error([&] {
+    compute(ener, force, virial, atom_energy, atom_virial, coord, atype, box,
+            nghost, inlist, ago, fparam, aparam, {}, atomic);
+  });
+}
+void DeepPotPT::computew(std::vector<double>& ener,
+                         std::vector<float>& force,
+                         std::vector<float>& virial,
+                         std::vector<float>& atom_energy,
+                         std::vector<float>& atom_virial,
+                         const std::vector<float>& coord,
+                         const std::vector<int>& atype,
+                         const std::vector<float>& box,
+                         const int nghost,
+                         const InputNlist& inlist,
+                         const int& ago,
+                         const std::vector<float>& fparam,
+                         const std::vector<float>& aparam,
+                         const bool atomic) {
+  translate_error([&] {
+    compute(ener, force, virial, atom_energy, atom_virial, coord, atype, box,
+            nghost, inlist, ago, fparam, aparam, {}, atomic);
+  });
+}
+// charge_spin overloads — thread runtime charge/spin through to compute()
+void DeepPotPT::computew(std::vector<double>& ener,
+                         std::vector<double>& force,
+                         std::vector<double>& virial,
+                         std::vector<double>& atom_energy,
+                         std::vector<double>& atom_virial,
+                         const std::vector<double>& coord,
+                         const std::vector<int>& atype,
+                         const std::vector<double>& box,
+                         const std::vector<double>& fparam,
+                         const std::vector<double>& aparam,
+                         const std::vector<double>& charge_spin,
+                         const bool atomic) {
+  translate_error([&] {
+    compute(ener, force, virial, atom_energy, atom_virial, coord, atype, box,
+            fparam, aparam, charge_spin, atomic);
   });
 }
 void DeepPotPT::computew(std::vector<double>& ener,
@@ -485,10 +680,11 @@ void DeepPotPT::computew(std::vector<double>& ener,
                          const std::vector<float>& box,
                          const std::vector<float>& fparam,
                          const std::vector<float>& aparam,
+                         const std::vector<double>& charge_spin,
                          const bool atomic) {
   translate_error([&] {
     compute(ener, force, virial, atom_energy, atom_virial, coord, atype, box,
-            fparam, aparam, atomic);
+            fparam, aparam, charge_spin, atomic);
   });
 }
 void DeepPotPT::computew(std::vector<double>& ener,
@@ -504,10 +700,11 @@ void DeepPotPT::computew(std::vector<double>& ener,
                          const int& ago,
                          const std::vector<double>& fparam,
                          const std::vector<double>& aparam,
+                         const std::vector<double>& charge_spin,
                          const bool atomic) {
   translate_error([&] {
     compute(ener, force, virial, atom_energy, atom_virial, coord, atype, box,
-            nghost, inlist, ago, fparam, aparam, atomic);
+            nghost, inlist, ago, fparam, aparam, charge_spin, atomic);
   });
 }
 void DeepPotPT::computew(std::vector<double>& ener,
@@ -523,10 +720,11 @@ void DeepPotPT::computew(std::vector<double>& ener,
                          const int& ago,
                          const std::vector<float>& fparam,
                          const std::vector<float>& aparam,
+                         const std::vector<double>& charge_spin,
                          const bool atomic) {
   translate_error([&] {
     compute(ener, force, virial, atom_energy, atom_virial, coord, atype, box,
-            nghost, inlist, ago, fparam, aparam, atomic);
+            nghost, inlist, ago, fparam, aparam, charge_spin, atomic);
   });
 }
 void DeepPotPT::computew_mixed_type(std::vector<double>& ener,
diff --git a/source/api_cc/src/DeepPotPTExpt.cc b/source/api_cc/src/DeepPotPTExpt.cc
index c8c1bfcfad..caf73416c8 100644
--- a/source/api_cc/src/DeepPotPTExpt.cc
+++ b/source/api_cc/src/DeepPotPTExpt.cc
@@ -104,9 +104,9 @@ void DeepPotPTExpt::init(const std::string& model,
                : static_cast<int>(metadata["type_map"].as_array().size());
   dfparam = metadata["dim_fparam"].as_int();
   daparam = metadata["dim_aparam"].as_int();
-  dim_chg_spin = metadata.obj_val.count("dim_chg_spin")
-                     ? metadata["dim_chg_spin"].as_int()
-                     : 0;
+  dchgspin = metadata.obj_val.count("dim_chg_spin")
+                 ? metadata["dim_chg_spin"].as_int()
+                 : 0;
   aparam_nall = false;  // pt_expt models use nloc for aparam
   if (metadata.obj_val.count("has_default_fparam")) {
     has_default_fparam_ = metadata["has_default_fparam"].as_bool();
@@ -132,7 +132,7 @@ void DeepPotPTExpt::init(const std::string& model,
                 << std::endl;
     }
   }
-  default_chg_spin_ = read_default_chg_spin(metadata, dim_chg_spin);
+  default_chg_spin_ = read_default_chg_spin(metadata, dchgspin);
 
   if (metadata.obj_val.count("do_atomic_virial")) {
     do_atomic_virial = metadata["do_atomic_virial"].as_bool();
@@ -241,10 +241,10 @@ std::vector<torch::Tensor> DeepPotPTExpt::run_model(
     const torch::Tensor& nlist,
     const torch::Tensor& mapping,
     const torch::Tensor& fparam,
-    const torch::Tensor& aparam) {
-  // Only include fparam/aparam if the model was exported with them.
-  // When fparam/aparam are None at export time, AOTInductor compiles
-  // the model with fewer inputs (e.g. 4 instead of 6).
+    const torch::Tensor& aparam,
+    const torch::Tensor& charge_spin) {
+  // Only include fparam/aparam/charge_spin if the model was exported with them.
+  // When they are None at export time, AOTInductor compiles with fewer inputs.
   std::vector<torch::Tensor> inputs = {coord, atype, nlist, mapping};
   if (dfparam > 0) {
     inputs.push_back(fparam);
@@ -252,11 +252,7 @@ std::vector<torch::Tensor> DeepPotPTExpt::run_model(
   if (daparam > 0) {
     inputs.push_back(aparam);
   }
-  if (dim_chg_spin > 0) {
-    auto charge_spin = torch::tensor(default_chg_spin_, coord.options())
-                           .view({1, dim_chg_spin})
-                           .expand({coord.size(0), dim_chg_spin})
-                           .contiguous();
+  if (dchgspin > 0) {
     inputs.push_back(charge_spin);
   }
   return loader->run(inputs);
@@ -269,6 +265,7 @@ std::vector<torch::Tensor> DeepPotPTExpt::run_model_with_comm(
     const torch::Tensor& mapping,
     const torch::Tensor& fparam,
     const torch::Tensor& aparam,
+    const torch::Tensor& charge_spin,
     const std::vector<at::Tensor>& comm_tensors) {
   if (!with_comm_loader) {
     throw deepmd::deepmd_exception(
@@ -293,11 +290,7 @@ std::vector<torch::Tensor> DeepPotPTExpt::run_model_with_comm(
   if (daparam > 0) {
     inputs.push_back(aparam);
   }
-  if (dim_chg_spin > 0) {
-    auto charge_spin = torch::tensor(default_chg_spin_, coord.options())
-                           .view({1, dim_chg_spin})
-                           .expand({coord.size(0), dim_chg_spin})
-                           .contiguous();
+  if (dchgspin > 0) {
     inputs.push_back(charge_spin);
   }
   for (const auto& t : comm_tensors) {
@@ -334,6 +327,7 @@ void DeepPotPTExpt::compute(ENERGYVTYPE& ener,
                             const int& ago,
                             const std::vector<VALUETYPE>& fparam,
                             const std::vector<VALUETYPE>& aparam,
+                            const std::vector<double>& charge_spin,
                             const bool atomic) {
   // Fail fast before allocating any tensors: refuse to run if the caller
   // asked for atomic virial but the .pt2 was exported without it.
@@ -514,6 +508,39 @@ void DeepPotPTExpt::compute(ENERGYVTYPE& ener,
     aparam_tensor = torch::zeros({0}, options).to(device);
   }
 
+  // Build charge_spin tensor: use runtime value when provided, fall back to
+  // default_chg_spin_ stored in the .pt2 metadata.
+  at::Tensor charge_spin_tensor;
+  if (dchgspin > 0) {
+    auto dbl_options = torch::TensorOptions().dtype(torch::kFloat64);
+    if (!charge_spin.empty()) {
+      // Single-frame path: charge_spin must hold exactly dim_chg_spin values.
+      if (static_cast<int>(charge_spin.size()) != dchgspin) {
+        throw deepmd::deepmd_exception(
+            "charge_spin has " + std::to_string(charge_spin.size()) +
+            " values but the model expects dim_chg_spin=" +
+            std::to_string(dchgspin) + ".");
+      }
+      charge_spin_tensor =
+          torch::from_blob(const_cast<double*>(charge_spin.data()),
+                           {1, static_cast<std::int64_t>(charge_spin.size())},
+                           dbl_options)
+              .clone()
+              .to(device);
+    } else if (!default_chg_spin_.empty()) {
+      charge_spin_tensor =
+          torch::from_blob(const_cast<double*>(default_chg_spin_.data()),
+                           {1, dchgspin}, dbl_options)
+              .clone()
+              .to(device);
+    } else {
+      throw deepmd::deepmd_exception(
+          "charge_spin is empty and no default_chg_spin is available in the "
+          ".pt2 metadata. Provide charge_spin explicitly or regenerate the "
+          "model with a default charge/spin value.");
+    }
+  }
+
   // ``use_with_comm`` was computed earlier alongside the fail-fast
   // dispatch check.  Use the with-comm artifact for the multi-rank case
   // (the regular artifact uses the mapping tensor to gather ghost
@@ -545,10 +572,11 @@ void DeepPotPTExpt::compute(ENERGYVTYPE& ener,
     }
     flat_outputs = run_model_with_comm(
         coord_Tensor, atype_Tensor, firstneigh_tensor, mapping_tensor,
-        fparam_tensor, aparam_tensor, comm_tensors);
+        fparam_tensor, aparam_tensor, charge_spin_tensor, comm_tensors);
   } else {
-    flat_outputs = run_model(coord_Tensor, atype_Tensor, firstneigh_tensor,
-                             mapping_tensor, fparam_tensor, aparam_tensor);
+    flat_outputs =
+        run_model(coord_Tensor, atype_Tensor, firstneigh_tensor, mapping_tensor,
+                  fparam_tensor, aparam_tensor, charge_spin_tensor);
   }
 
   // Map flat outputs to internal keys
@@ -624,6 +652,7 @@ template void DeepPotPTExpt::compute<double, std::vector<ENERGYTYPE>>(
     const int& ago,
     const std::vector<double>& fparam,
     const std::vector<double>& aparam,
+    const std::vector<double>& charge_spin,
     const bool atomic);
 template void DeepPotPTExpt::compute<float, std::vector<ENERGYTYPE>>(
     std::vector<ENERGYTYPE>& ener,
@@ -639,6 +668,7 @@ template void DeepPotPTExpt::compute<float, std::vector<ENERGYTYPE>>(
     const int& ago,
     const std::vector<float>& fparam,
     const std::vector<float>& aparam,
+    const std::vector<double>& charge_spin,
     const bool atomic);
 
 template <typename VALUETYPE, typename ENERGYVTYPE>
@@ -652,6 +682,7 @@ void DeepPotPTExpt::compute(ENERGYVTYPE& ener,
                             const std::vector<VALUETYPE>& box,
                             const std::vector<VALUETYPE>& fparam,
                             const std::vector<VALUETYPE>& aparam,
+                            const std::vector<double>& charge_spin,
                             const bool atomic) {
   // Fail fast before allocating any tensors (same check as the nlist
   // overload — see its comment).
@@ -668,7 +699,7 @@ void DeepPotPTExpt::compute(ENERGYVTYPE& ener,
   if (nframes > 1) {
     // Multi-frame: loop over frames and concatenate
     compute_nframes(ener, force, virial, atom_energy, atom_virial, nframes,
-                    coord, atype, box, fparam, aparam, atomic);
+                    coord, atype, box, fparam, aparam, charge_spin, atomic);
     return;
   }
   // The .pt2 model only contains forward_common_lower, which requires
@@ -808,9 +839,43 @@ void DeepPotPTExpt::compute(ENERGYVTYPE& ener,
     aparam_tensor = torch::zeros({0}, options).to(device);
   }
 
+  // Build charge_spin tensor: use runtime value when provided, fall back to
+  // default_chg_spin_ stored in the .pt2 metadata.
+  at::Tensor charge_spin_tensor;
+  if (dchgspin > 0) {
+    auto dbl_options = torch::TensorOptions().dtype(torch::kFloat64);
+    if (!charge_spin.empty()) {
+      // Single-frame path: charge_spin must hold exactly dim_chg_spin values.
+      if (static_cast<int>(charge_spin.size()) != dchgspin) {
+        throw deepmd::deepmd_exception(
+            "charge_spin has " + std::to_string(charge_spin.size()) +
+            " values but the model expects dim_chg_spin=" +
+            std::to_string(dchgspin) + ".");
+      }
+      charge_spin_tensor =
+          torch::from_blob(const_cast<double*>(charge_spin.data()),
+                           {1, static_cast<std::int64_t>(charge_spin.size())},
+                           dbl_options)
+              .clone()
+              .to(device);
+    } else if (!default_chg_spin_.empty()) {
+      charge_spin_tensor =
+          torch::from_blob(const_cast<double*>(default_chg_spin_.data()),
+                           {1, dchgspin}, dbl_options)
+              .clone()
+              .to(device);
+    } else {
+      throw deepmd::deepmd_exception(
+          "charge_spin is empty and no default_chg_spin is available in the "
+          ".pt2 metadata. Provide charge_spin explicitly or regenerate the "
+          "model with a default charge/spin value.");
+    }
+  }
+
   // 5. Run the .pt2 model
-  auto flat_outputs = run_model(coord_Tensor, atype_Tensor, nlist_tensor,
-                                mapping_tensor, fparam_tensor, aparam_tensor);
+  auto flat_outputs =
+      run_model(coord_Tensor, atype_Tensor, nlist_tensor, mapping_tensor,
+                fparam_tensor, aparam_tensor, charge_spin_tensor);
 
   // 6. Map flat outputs to internal keys
   std::map<std::string, torch::Tensor> output_map;
@@ -874,10 +939,26 @@ void DeepPotPTExpt::compute_nframes(ENERGYVTYPE& ener,
                                     const std::vector<VALUETYPE>& box,
                                     const std::vector<VALUETYPE>& fparam,
                                     const std::vector<VALUETYPE>& aparam,
+                                    const std::vector<double>& charge_spin,
                                     const bool atomic) {
   int natoms = atype.size();
   int dap = aparam.empty() ? 0 : static_cast<int>(aparam.size()) / nframes;
   int dfp = fparam.empty() ? 0 : static_cast<int>(fparam.size()) / nframes;
+  // charge_spin may be empty (default fallback), a single dim_chg_spin vector
+  // (broadcast to all frames), or nframes * dim_chg_spin (per-frame). Reject
+  // anything else up-front to avoid out-of-range slicing in the loop.
+  if (!charge_spin.empty()) {
+    size_t s_dcsp = static_cast<size_t>(dchgspin);
+    if (charge_spin.size() != s_dcsp &&
+        charge_spin.size() != s_dcsp * static_cast<size_t>(nframes)) {
+      throw deepmd::deepmd_exception(
+          "charge_spin has " + std::to_string(charge_spin.size()) +
+          " values but the model expects dim_chg_spin=" +
+          std::to_string(dchgspin) + " (per frame) or " +
+          std::to_string(dchgspin * nframes) + " (for " +
+          std::to_string(nframes) + " frames).");
+    }
+  }
   ener.clear();
   force.clear();
   virial.clear();
@@ -907,10 +988,22 @@ void DeepPotPTExpt::compute_nframes(ENERGYVTYPE& ener,
       frame_aparam.assign(aparam.begin() + s_ff * s_dap,
                           aparam.begin() + (s_ff + 1) * s_dap);
     }
+    std::vector<double> frame_chg_spin;
+    if (!charge_spin.empty()) {
+      size_t s_dcsp = static_cast<size_t>(dchgspin);
+      if (charge_spin.size() == s_dcsp) {
+        // single charge/spin vector broadcast to every frame
+        frame_chg_spin = charge_spin;
+      } else {
+        frame_chg_spin.assign(charge_spin.begin() + s_ff * s_dcsp,
+                              charge_spin.begin() + (s_ff + 1) * s_dcsp);
+      }
+    }
     std::vector<ENERGYTYPE> frame_ener;
     std::vector<VALUETYPE> frame_force, frame_virial, frame_ae, frame_av;
     compute(frame_ener, frame_force, frame_virial, frame_ae, frame_av,
-            frame_coord, atype, frame_box, frame_fparam, frame_aparam, atomic);
+            frame_coord, atype, frame_box, frame_fparam, frame_aparam,
+            frame_chg_spin, atomic);
     ener.insert(ener.end(), frame_ener.begin(), frame_ener.end());
     force.insert(force.end(), frame_force.begin(), frame_force.end());
     virial.insert(virial.end(), frame_virial.begin(), frame_virial.end());
@@ -932,6 +1025,7 @@ template void DeepPotPTExpt::compute<double, std::vector<ENERGYTYPE>>(
     const std::vector<double>& box,
     const std::vector<double>& fparam,
     const std::vector<double>& aparam,
+    const std::vector<double>& charge_spin,
     const bool atomic);
 template void DeepPotPTExpt::compute<float, std::vector<ENERGYTYPE>>(
     std::vector<ENERGYTYPE>& ener,
@@ -944,6 +1038,7 @@ template void DeepPotPTExpt::compute<float, std::vector<ENERGYTYPE>>(
     const std::vector<float>& box,
     const std::vector<float>& fparam,
     const std::vector<float>& aparam,
+    const std::vector<double>& charge_spin,
     const bool atomic);
 
 void DeepPotPTExpt::get_type_map(std::string& type_map_str) {
@@ -953,7 +1048,7 @@ void DeepPotPTExpt::get_type_map(std::string& type_map_str) {
   }
 }
 
-// forward to template method
+// forward to template method (no runtime charge_spin — uses default_chg_spin_)
 void DeepPotPTExpt::computew(std::vector<double>& ener,
                              std::vector<double>& force,
                              std::vector<double>& virial,
@@ -967,7 +1062,7 @@ void DeepPotPTExpt::computew(std::vector<double>& ener,
                              const bool atomic) {
   translate_error([&] {
     compute(ener, force, virial, atom_energy, atom_virial, coord, atype, box,
-            fparam, aparam, atomic);
+            fparam, aparam, {}, atomic);
   });
 }
 void DeepPotPTExpt::computew(std::vector<double>& ener,
@@ -983,7 +1078,7 @@ void DeepPotPTExpt::computew(std::vector<double>& ener,
                              const bool atomic) {
   translate_error([&] {
     compute(ener, force, virial, atom_energy, atom_virial, coord, atype, box,
-            fparam, aparam, atomic);
+            fparam, aparam, {}, atomic);
   });
 }
 void DeepPotPTExpt::computew(std::vector<double>& ener,
@@ -1002,7 +1097,7 @@ void DeepPotPTExpt::computew(std::vector<double>& ener,
                              const bool atomic) {
   translate_error([&] {
     compute(ener, force, virial, atom_energy, atom_virial, coord, atype, box,
-            nghost, inlist, ago, fparam, aparam, atomic);
+            nghost, inlist, ago, fparam, aparam, {}, atomic);
   });
 }
 void DeepPotPTExpt::computew(std::vector<double>& ener,
@@ -1021,7 +1116,7 @@ void DeepPotPTExpt::computew(std::vector<double>& ener,
                              const bool atomic) {
   translate_error([&] {
     compute(ener, force, virial, atom_energy, atom_virial, coord, atype, box,
-            nghost, inlist, ago, fparam, aparam, atomic);
+            nghost, inlist, ago, fparam, aparam, {}, atomic);
   });
 }
 template <typename VALUETYPE>
@@ -1037,12 +1132,28 @@ void DeepPotPTExpt::compute_mixed_type_impl(
     const std::vector<VALUETYPE>& box,
     const std::vector<VALUETYPE>& fparam,
     const std::vector<VALUETYPE>& aparam,
+    const std::vector<double>& charge_spin,
     const bool atomic) {
   // Mixed-type: atype has nframes * natoms elements.
   // Loop over frames, each with its own atype slice.
   int natoms = static_cast<int>(atype.size()) / nframes;
   int dap = aparam.empty() ? 0 : static_cast<int>(aparam.size()) / nframes;
   int dfp = fparam.empty() ? 0 : static_cast<int>(fparam.size()) / nframes;
+  // charge_spin may be empty (default fallback), a single dim_chg_spin vector
+  // (broadcast to all frames), or nframes * dim_chg_spin (per-frame). Reject
+  // anything else up-front to avoid out-of-range slicing in the loop.
+  if (!charge_spin.empty()) {
+    size_t s_dcsp = static_cast<size_t>(dchgspin);
+    if (charge_spin.size() != s_dcsp &&
+        charge_spin.size() != s_dcsp * static_cast<size_t>(nframes)) {
+      throw deepmd::deepmd_exception(
+          "charge_spin has " + std::to_string(charge_spin.size()) +
+          " values but the model expects dim_chg_spin=" +
+          std::to_string(dchgspin) + " (per frame) or " +
+          std::to_string(dchgspin * nframes) + " (for " +
+          std::to_string(nframes) + " frames).");
+    }
+  }
   ener.clear();
   force.clear();
   virial.clear();
@@ -1074,11 +1185,22 @@ void DeepPotPTExpt::compute_mixed_type_impl(
       frame_aparam.assign(aparam.begin() + s_ff * s_dap,
                           aparam.begin() + (s_ff + 1) * s_dap);
     }
+    std::vector<double> frame_chg_spin;
+    if (!charge_spin.empty()) {
+      size_t s_dcsp = static_cast<size_t>(dchgspin);
+      if (charge_spin.size() == s_dcsp) {
+        // single charge/spin vector broadcast to every frame
+        frame_chg_spin = charge_spin;
+      } else {
+        frame_chg_spin.assign(charge_spin.begin() + s_ff * s_dcsp,
+                              charge_spin.begin() + (s_ff + 1) * s_dcsp);
+      }
+    }
     std::vector<ENERGYTYPE> frame_ener;
     std::vector<VALUETYPE> frame_force, frame_virial, frame_ae, frame_av;
     compute(frame_ener, frame_force, frame_virial, frame_ae, frame_av,
             frame_coord, frame_atype, frame_box, frame_fparam, frame_aparam,
-            atomic);
+            frame_chg_spin, atomic);
     ener.insert(ener.end(), frame_ener.begin(), frame_ener.end());
     force.insert(force.end(), frame_force.begin(), frame_force.end());
     virial.insert(virial.end(), frame_virial.begin(), frame_virial.end());
@@ -1103,7 +1225,121 @@ void DeepPotPTExpt::computew_mixed_type(std::vector<double>& ener,
                                         const bool atomic) {
   translate_error([&] {
     compute_mixed_type_impl(ener, force, virial, atom_energy, atom_virial,
-                            nframes, coord, atype, box, fparam, aparam, atomic);
+                            nframes, coord, atype, box, fparam, aparam, {},
+                            atomic);
+  });
+}
+void DeepPotPTExpt::computew_mixed_type(std::vector<double>& ener,
+                                        std::vector<float>& force,
+                                        std::vector<float>& virial,
+                                        std::vector<float>& atom_energy,
+                                        std::vector<float>& atom_virial,
+                                        const int& nframes,
+                                        const std::vector<float>& coord,
+                                        const std::vector<int>& atype,
+                                        const std::vector<float>& box,
+                                        const std::vector<float>& fparam,
+                                        const std::vector<float>& aparam,
+                                        const bool atomic) {
+  translate_error([&] {
+    compute_mixed_type_impl(ener, force, virial, atom_energy, atom_virial,
+                            nframes, coord, atype, box, fparam, aparam, {},
+                            atomic);
+  });
+}
+
+// charge_spin overloads — pass runtime charge/spin value through to compute()
+void DeepPotPTExpt::computew(std::vector<double>& ener,
+                             std::vector<double>& force,
+                             std::vector<double>& virial,
+                             std::vector<double>& atom_energy,
+                             std::vector<double>& atom_virial,
+                             const std::vector<double>& coord,
+                             const std::vector<int>& atype,
+                             const std::vector<double>& box,
+                             const std::vector<double>& fparam,
+                             const std::vector<double>& aparam,
+                             const std::vector<double>& charge_spin,
+                             const bool atomic) {
+  translate_error([&] {
+    compute(ener, force, virial, atom_energy, atom_virial, coord, atype, box,
+            fparam, aparam, charge_spin, atomic);
+  });
+}
+void DeepPotPTExpt::computew(std::vector<double>& ener,
+                             std::vector<float>& force,
+                             std::vector<float>& virial,
+                             std::vector<float>& atom_energy,
+                             std::vector<float>& atom_virial,
+                             const std::vector<float>& coord,
+                             const std::vector<int>& atype,
+                             const std::vector<float>& box,
+                             const std::vector<float>& fparam,
+                             const std::vector<float>& aparam,
+                             const std::vector<double>& charge_spin,
+                             const bool atomic) {
+  translate_error([&] {
+    compute(ener, force, virial, atom_energy, atom_virial, coord, atype, box,
+            fparam, aparam, charge_spin, atomic);
+  });
+}
+void DeepPotPTExpt::computew(std::vector<double>& ener,
+                             std::vector<double>& force,
+                             std::vector<double>& virial,
+                             std::vector<double>& atom_energy,
+                             std::vector<double>& atom_virial,
+                             const std::vector<double>& coord,
+                             const std::vector<int>& atype,
+                             const std::vector<double>& box,
+                             const int nghost,
+                             const InputNlist& inlist,
+                             const int& ago,
+                             const std::vector<double>& fparam,
+                             const std::vector<double>& aparam,
+                             const std::vector<double>& charge_spin,
+                             const bool atomic) {
+  translate_error([&] {
+    compute(ener, force, virial, atom_energy, atom_virial, coord, atype, box,
+            nghost, inlist, ago, fparam, aparam, charge_spin, atomic);
+  });
+}
+void DeepPotPTExpt::computew(std::vector<double>& ener,
+                             std::vector<float>& force,
+                             std::vector<float>& virial,
+                             std::vector<float>& atom_energy,
+                             std::vector<float>& atom_virial,
+                             const std::vector<float>& coord,
+                             const std::vector<int>& atype,
+                             const std::vector<float>& box,
+                             const int nghost,
+                             const InputNlist& inlist,
+                             const int& ago,
+                             const std::vector<float>& fparam,
+                             const std::vector<float>& aparam,
+                             const std::vector<double>& charge_spin,
+                             const bool atomic) {
+  translate_error([&] {
+    compute(ener, force, virial, atom_energy, atom_virial, coord, atype, box,
+            nghost, inlist, ago, fparam, aparam, charge_spin, atomic);
+  });
+}
+void DeepPotPTExpt::computew_mixed_type(std::vector<double>& ener,
+                                        std::vector<double>& force,
+                                        std::vector<double>& virial,
+                                        std::vector<double>& atom_energy,
+                                        std::vector<double>& atom_virial,
+                                        const int& nframes,
+                                        const std::vector<double>& coord,
+                                        const std::vector<int>& atype,
+                                        const std::vector<double>& box,
+                                        const std::vector<double>& fparam,
+                                        const std::vector<double>& aparam,
+                                        const std::vector<double>& charge_spin,
+                                        const bool atomic) {
+  translate_error([&] {
+    compute_mixed_type_impl(ener, force, virial, atom_energy, atom_virial,
+                            nframes, coord, atype, box, fparam, aparam,
+                            charge_spin, atomic);
   });
 }
 void DeepPotPTExpt::computew_mixed_type(std::vector<double>& ener,
@@ -1117,10 +1353,12 @@ void DeepPotPTExpt::computew_mixed_type(std::vector<double>& ener,
                                         const std::vector<float>& box,
                                         const std::vector<float>& fparam,
                                         const std::vector<float>& aparam,
+                                        const std::vector<double>& charge_spin,
                                         const bool atomic) {
   translate_error([&] {
     compute_mixed_type_impl(ener, force, virial, atom_energy, atom_virial,
-                            nframes, coord, atype, box, fparam, aparam, atomic);
+                            nframes, coord, atype, box, fparam, aparam,
+                            charge_spin, atomic);
   });
 }
 #endif
diff --git a/source/api_cc/tests/test_deeppot_chg_spin_pt.cc b/source/api_cc/tests/test_deeppot_chg_spin_pt.cc
new file mode 100644
index 0000000000..296ee3023c
--- /dev/null
+++ b/source/api_cc/tests/test_deeppot_chg_spin_pt.cc
@@ -0,0 +1,295 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Test C++ inference for the PT (.pth) backend with charge_spin input (DPA3).
+// Model + reference generated by source/tests/infer/gen_chg_spin.py.
+//
+// Mirrors the DPA3 .pth test layout (test_deeppot_dpa3_pt.cc): the PBC case
+// exercises the standalone path (forward), and the NoPbc case exercises the
+// LAMMPS neighbor-list path (forward_lower) without ghost atoms.  Both cases
+// cover explicit charge_spin and the default_chg_spin fallback.
+#include <gtest/gtest.h>
+
+#include <algorithm>
+#include <cmath>
+#include <fstream>
+#include <vector>
+
+#include "DeepPot.h"
+#include "expected_ref.h"
+#include "neighbor_list.h"
+#include "test_utils.h"
+
+#undef EPSILON
+#define EPSILON (std::is_same<VALUETYPE, double>::value ? 1e-10 : 1e-4)
+
+namespace {
+constexpr const char* kRefPath = "../../tests/infer/chg_spin.expected";
+constexpr const char* kModelPath = "../../tests/infer/chg_spin.pth";
+}  // namespace
+
+// ---------------------------------------------------------------------------
+// PBC: standalone path (DeepPot::compute that builds its own nlist -> forward)
+// ---------------------------------------------------------------------------
+template <class VALUETYPE>
+class TestInferDeepPotChgSpinPt : public ::testing::Test {
+ protected:
+  std::vector<VALUETYPE> coord = {12.83, 2.56, 2.18, 12.09, 2.87, 2.74,
+                                  00.25, 3.32, 1.68, 3.36,  3.00, 1.81,
+                                  3.51,  2.51, 2.60, 4.27,  3.22, 1.56};
+  std::vector<int> atype = {0, 1, 1, 0, 1, 1};
+  std::vector<VALUETYPE> box = {13., 0., 0., 0., 13., 0., 0., 0., 13.};
+  // charge_spin is always double regardless of VALUETYPE.
+  // [1.0, 2.0] maps to charge idx 101, spin idx 2 — both differ from the
+  // model's default_chg_spin [0.0, 1.0] (charge idx 100, spin idx 1).
+  std::vector<double> charge_spin_explicit = {1.0, 2.0};
+
+  std::vector<VALUETYPE> expected_f_explicit, expected_f_default;
+  int natoms;
+  double expected_tot_e_explicit, expected_tot_e_default;
+  std::vector<VALUETYPE> expected_tot_v_explicit, expected_tot_v_default;
+
+  deepmd::DeepPot dp;
+
+  static void accumulate(const std::vector<VALUETYPE>& e,
+                         const std::vector<VALUETYPE>& v,
+                         double& tot_e,
+                         std::vector<VALUETYPE>& tot_v) {
+    int n = e.size();
+    tot_e = 0.;
+    tot_v.assign(9, 0.);
+    for (int ii = 0; ii < n; ++ii) {
+      tot_e += e[ii];
+      for (int dd = 0; dd < 9; ++dd) {
+        tot_v[dd] += v[ii * 9 + dd];
+      }
+    }
+  }
+
+  void SetUp() override {
+#ifndef BUILD_PYTORCH
+    GTEST_SKIP() << "Skip because PyTorch support is not enabled.";
+#endif
+    deepmd_test::ExpectedRef ref;
+    ref.load(kRefPath);
+    std::vector<VALUETYPE> e_exp = ref.get<VALUETYPE>("explicit", "expected_e");
+    std::vector<VALUETYPE> e_def = ref.get<VALUETYPE>("default", "expected_e");
+    expected_f_explicit = ref.get<VALUETYPE>("explicit", "expected_f");
+    expected_f_default = ref.get<VALUETYPE>("default", "expected_f");
+    accumulate(e_exp, ref.get<VALUETYPE>("explicit", "expected_v"),
+               expected_tot_e_explicit, expected_tot_v_explicit);
+    accumulate(e_def, ref.get<VALUETYPE>("default", "expected_v"),
+               expected_tot_e_default, expected_tot_v_default);
+    natoms = e_exp.size();
+
+    {
+      std::ifstream model_file(kModelPath);
+      if (!model_file.good()) {
+        GTEST_SKIP() << "Skip because " << kModelPath
+                     << " was not generated (e.g. .pth export was skipped).";
+      }
+    }
+    dp.init(kModelPath);
+  }
+
+  void TearDown() override {}
+};
+
+TYPED_TEST_SUITE(TestInferDeepPotChgSpinPt, ValueTypes);
+
+TYPED_TEST(TestInferDeepPotChgSpinPt, dim_chg_spin) {
+  deepmd::DeepPot& dp = this->dp;
+  EXPECT_EQ(dp.dim_chg_spin(), 2);
+}
+
+TYPED_TEST(TestInferDeepPotChgSpinPt, cpu_build_nlist_explicit) {
+  using VALUETYPE = TypeParam;
+  std::vector<VALUETYPE>& coord = this->coord;
+  std::vector<int>& atype = this->atype;
+  std::vector<VALUETYPE>& box = this->box;
+  std::vector<double>& charge_spin = this->charge_spin_explicit;
+  std::vector<VALUETYPE>& expected_f = this->expected_f_explicit;
+  std::vector<VALUETYPE>& expected_tot_v = this->expected_tot_v_explicit;
+  int& natoms = this->natoms;
+  double& expected_tot_e = this->expected_tot_e_explicit;
+  deepmd::DeepPot& dp = this->dp;
+
+  double ener;
+  std::vector<VALUETYPE> force, virial;
+  dp.compute(ener, force, virial, coord, atype, box, {}, {}, charge_spin);
+
+  EXPECT_EQ(force.size(), static_cast<size_t>(natoms * 3));
+  EXPECT_EQ(virial.size(), 9u);
+  EXPECT_LT(fabs(ener - expected_tot_e), EPSILON);
+  for (int ii = 0; ii < natoms * 3; ++ii) {
+    EXPECT_LT(fabs(force[ii] - expected_f[ii]), EPSILON);
+  }
+  for (int ii = 0; ii < 9; ++ii) {
+    EXPECT_LT(fabs(virial[ii] - expected_tot_v[ii]), EPSILON);
+  }
+}
+
+TYPED_TEST(TestInferDeepPotChgSpinPt, cpu_build_nlist_default) {
+  using VALUETYPE = TypeParam;
+  std::vector<VALUETYPE>& coord = this->coord;
+  std::vector<int>& atype = this->atype;
+  std::vector<VALUETYPE>& box = this->box;
+  std::vector<VALUETYPE>& expected_f = this->expected_f_default;
+  std::vector<VALUETYPE>& expected_tot_v = this->expected_tot_v_default;
+  int& natoms = this->natoms;
+  double& expected_tot_e = this->expected_tot_e_default;
+  deepmd::DeepPot& dp = this->dp;
+
+  double ener;
+  std::vector<VALUETYPE> force, virial;
+  // Empty charge_spin triggers the default_chg_spin fallback ([0.0, 1.0]).
+  dp.compute(ener, force, virial, coord, atype, box);
+
+  EXPECT_EQ(force.size(), static_cast<size_t>(natoms * 3));
+  EXPECT_EQ(virial.size(), 9u);
+  EXPECT_LT(fabs(ener - expected_tot_e), EPSILON);
+  for (int ii = 0; ii < natoms * 3; ++ii) {
+    EXPECT_LT(fabs(force[ii] - expected_f[ii]), EPSILON);
+  }
+  for (int ii = 0; ii < 9; ++ii) {
+    EXPECT_LT(fabs(virial[ii] - expected_tot_v[ii]), EPSILON);
+  }
+}
+
+// ---------------------------------------------------------------------------
+// NoPbc: LAMMPS neighbor-list path (forward_lower), no ghost atoms.
+// Mirrors test_deeppot_dpa3_pt.cc::TestInferDeepPotDpa3PtNoPbc.cpu_lmp_nlist.
+// ---------------------------------------------------------------------------
+template <class VALUETYPE>
+class TestInferDeepPotChgSpinPtNoPbc : public ::testing::Test {
+ protected:
+  std::vector<VALUETYPE> coord = {12.83, 2.56, 2.18, 12.09, 2.87, 2.74,
+                                  00.25, 3.32, 1.68, 3.36,  3.00, 1.81,
+                                  3.51,  2.51, 2.60, 4.27,  3.22, 1.56};
+  std::vector<int> atype = {0, 1, 1, 0, 1, 1};
+  std::vector<VALUETYPE> box = {};
+  std::vector<double> charge_spin_explicit = {1.0, 2.0};
+
+  std::vector<VALUETYPE> expected_f_explicit, expected_f_default;
+  int natoms;
+  double expected_tot_e_explicit, expected_tot_e_default;
+  std::vector<VALUETYPE> expected_tot_v_explicit, expected_tot_v_default;
+
+  deepmd::DeepPot dp;
+
+  static void accumulate(const std::vector<VALUETYPE>& e,
+                         const std::vector<VALUETYPE>& v,
+                         double& tot_e,
+                         std::vector<VALUETYPE>& tot_v) {
+    int n = e.size();
+    tot_e = 0.;
+    tot_v.assign(9, 0.);
+    for (int ii = 0; ii < n; ++ii) {
+      tot_e += e[ii];
+      for (int dd = 0; dd < 9; ++dd) {
+        tot_v[dd] += v[ii * 9 + dd];
+      }
+    }
+  }
+
+  void SetUp() override {
+#ifndef BUILD_PYTORCH
+    GTEST_SKIP() << "Skip because PyTorch support is not enabled.";
+#endif
+    deepmd_test::ExpectedRef ref;
+    ref.load(kRefPath);
+    std::vector<VALUETYPE> e_exp =
+        ref.get<VALUETYPE>("nopbc_explicit", "expected_e");
+    std::vector<VALUETYPE> e_def =
+        ref.get<VALUETYPE>("nopbc_default", "expected_e");
+    expected_f_explicit = ref.get<VALUETYPE>("nopbc_explicit", "expected_f");
+    expected_f_default = ref.get<VALUETYPE>("nopbc_default", "expected_f");
+    accumulate(e_exp, ref.get<VALUETYPE>("nopbc_explicit", "expected_v"),
+               expected_tot_e_explicit, expected_tot_v_explicit);
+    accumulate(e_def, ref.get<VALUETYPE>("nopbc_default", "expected_v"),
+               expected_tot_e_default, expected_tot_v_default);
+    natoms = e_exp.size();
+
+    {
+      std::ifstream model_file(kModelPath);
+      if (!model_file.good()) {
+        GTEST_SKIP() << "Skip because " << kModelPath
+                     << " was not generated (e.g. .pth export was skipped).";
+      }
+    }
+    dp.init(kModelPath);
+  }
+
+  void TearDown() override {}
+};
+
+TYPED_TEST_SUITE(TestInferDeepPotChgSpinPtNoPbc, ValueTypes);
+
+TYPED_TEST(TestInferDeepPotChgSpinPtNoPbc, cpu_lmp_nlist_explicit) {
+  using VALUETYPE = TypeParam;
+  std::vector<VALUETYPE>& coord = this->coord;
+  std::vector<int>& atype = this->atype;
+  std::vector<VALUETYPE>& box = this->box;
+  std::vector<double>& charge_spin = this->charge_spin_explicit;
+  std::vector<VALUETYPE>& expected_f = this->expected_f_explicit;
+  std::vector<VALUETYPE>& expected_tot_v = this->expected_tot_v_explicit;
+  int& natoms = this->natoms;
+  double& expected_tot_e = this->expected_tot_e_explicit;
+  deepmd::DeepPot& dp = this->dp;
+
+  std::vector<std::vector<int>> nlist_data = {{1, 2, 3, 4, 5}, {0, 2, 3, 4, 5},
+                                              {0, 1, 3, 4, 5}, {0, 1, 2, 4, 5},
+                                              {0, 1, 2, 3, 5}, {0, 1, 2, 3, 4}};
+  std::vector<int> ilist(natoms), numneigh(natoms);
+  std::vector<int*> firstneigh(natoms);
+  deepmd::InputNlist inlist(natoms, &ilist[0], &numneigh[0], &firstneigh[0]);
+  convert_nlist(inlist, nlist_data);
+
+  double ener;
+  std::vector<VALUETYPE> force, virial;
+  dp.compute(ener, force, virial, coord, atype, box, 0, inlist, 0, {}, {},
+             charge_spin);
+
+  EXPECT_EQ(force.size(), static_cast<size_t>(natoms * 3));
+  EXPECT_EQ(virial.size(), 9u);
+  EXPECT_LT(fabs(ener - expected_tot_e), EPSILON);
+  for (int ii = 0; ii < natoms * 3; ++ii) {
+    EXPECT_LT(fabs(force[ii] - expected_f[ii]), EPSILON);
+  }
+  for (int ii = 0; ii < 9; ++ii) {
+    EXPECT_LT(fabs(virial[ii] - expected_tot_v[ii]), EPSILON);
+  }
+}
+
+TYPED_TEST(TestInferDeepPotChgSpinPtNoPbc, cpu_lmp_nlist_default) {
+  using VALUETYPE = TypeParam;
+  std::vector<VALUETYPE>& coord = this->coord;
+  std::vector<int>& atype = this->atype;
+  std::vector<VALUETYPE>& box = this->box;
+  std::vector<VALUETYPE>& expected_f = this->expected_f_default;
+  std::vector<VALUETYPE>& expected_tot_v = this->expected_tot_v_default;
+  int& natoms = this->natoms;
+  double& expected_tot_e = this->expected_tot_e_default;
+  deepmd::DeepPot& dp = this->dp;
+
+  std::vector<std::vector<int>> nlist_data = {{1, 2, 3, 4, 5}, {0, 2, 3, 4, 5},
+                                              {0, 1, 3, 4, 5}, {0, 1, 2, 4, 5},
+                                              {0, 1, 2, 3, 5}, {0, 1, 2, 3, 4}};
+  std::vector<int> ilist(natoms), numneigh(natoms);
+  std::vector<int*> firstneigh(natoms);
+  deepmd::InputNlist inlist(natoms, &ilist[0], &numneigh[0], &firstneigh[0]);
+  convert_nlist(inlist, nlist_data);
+
+  double ener;
+  std::vector<VALUETYPE> force, virial;
+  // Empty charge_spin triggers the default_chg_spin fallback ([0.0, 1.0]).
+  dp.compute(ener, force, virial, coord, atype, box, 0, inlist, 0);
+
+  EXPECT_EQ(force.size(), static_cast<size_t>(natoms * 3));
+  EXPECT_EQ(virial.size(), 9u);
+  EXPECT_LT(fabs(ener - expected_tot_e), EPSILON);
+  for (int ii = 0; ii < natoms * 3; ++ii) {
+    EXPECT_LT(fabs(force[ii] - expected_f[ii]), EPSILON);
+  }
+  for (int ii = 0; ii < 9; ++ii) {
+    EXPECT_LT(fabs(virial[ii] - expected_tot_v[ii]), EPSILON);
+  }
+}
diff --git a/source/api_cc/tests/test_deeppot_chg_spin_ptexpt.cc b/source/api_cc/tests/test_deeppot_chg_spin_ptexpt.cc
new file mode 100644
index 0000000000..f9db1a0f2b
--- /dev/null
+++ b/source/api_cc/tests/test_deeppot_chg_spin_ptexpt.cc
@@ -0,0 +1,263 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Test C++ inference for pt_expt (.pt2) backend with charge_spin input (DPA3).
+// Model generated by source/tests/infer/gen_chg_spin.py.
+#include <gtest/gtest.h>
+
+#include <algorithm>
+#include <cmath>
+#include <vector>
+
+#include "DeepPot.h"
+#include "DeepPotPTExpt.h"
+#include "expected_ref.h"
+#include "neighbor_list.h"
+#include "test_utils.h"
+
+#undef EPSILON
+#define EPSILON (std::is_same<VALUETYPE, double>::value ? 1e-7 : 1e-4)
+
+namespace {
+constexpr const char* kRefPath = "../../tests/infer/chg_spin.expected";
+constexpr const char* kModelPath = "../../tests/infer/chg_spin.pt2";
+}  // namespace
+
+template <class VALUETYPE>
+class TestInferDeepPotChgSpinPtExpt : public ::testing::Test {
+ protected:
+  std::vector<VALUETYPE> coord = {12.83, 2.56, 2.18, 12.09, 2.87, 2.74,
+                                  00.25, 3.32, 1.68, 3.36,  3.00, 1.81,
+                                  3.51,  2.51, 2.60, 4.27,  3.22, 1.56};
+  std::vector<int> atype = {0, 1, 1, 0, 1, 1};
+  std::vector<VALUETYPE> box = {13., 0., 0., 0., 13., 0., 0., 0., 13.};
+  // charge_spin is always double regardless of VALUETYPE.
+  // [1.0, 2.0] maps to charge idx 101, spin idx 2 — both differ from the
+  // model's default_chg_spin [0.0, 1.0] (charge idx 100, spin idx 1).
+  std::vector<double> charge_spin_explicit = {1.0, 2.0};
+
+  std::vector<VALUETYPE> expected_e_explicit;
+  std::vector<VALUETYPE> expected_f_explicit;
+  std::vector<VALUETYPE> expected_v_explicit;
+  std::vector<VALUETYPE> expected_e_default;
+  std::vector<VALUETYPE> expected_f_default;
+  std::vector<VALUETYPE> expected_v_default;
+  int natoms;
+  double expected_tot_e_explicit;
+  double expected_tot_e_default;
+  std::vector<VALUETYPE> expected_tot_v_explicit;
+  std::vector<VALUETYPE> expected_tot_v_default;
+
+  static deepmd::DeepPot dp;
+
+  static void SetUpTestSuite() {
+#if defined(BUILD_PYTORCH) && BUILD_PT_EXPT
+    dp.init(kModelPath);
+#endif
+  }
+
+  void SetUp() override {
+#if !defined(BUILD_PYTORCH) || !BUILD_PT_EXPT
+    GTEST_SKIP() << "Skip because PyTorch support is not enabled.";
+#endif
+    deepmd_test::ExpectedRef ref;
+    ref.load(kRefPath);
+    expected_e_explicit = ref.get<VALUETYPE>("explicit", "expected_e");
+    expected_f_explicit = ref.get<VALUETYPE>("explicit", "expected_f");
+    expected_v_explicit = ref.get<VALUETYPE>("explicit", "expected_v");
+    expected_e_default = ref.get<VALUETYPE>("default", "expected_e");
+    expected_f_default = ref.get<VALUETYPE>("default", "expected_f");
+    expected_v_default = ref.get<VALUETYPE>("default", "expected_v");
+
+    natoms = expected_e_explicit.size();
+    EXPECT_EQ(natoms * 3, static_cast<int>(expected_f_explicit.size()));
+    EXPECT_EQ(natoms * 9, static_cast<int>(expected_v_explicit.size()));
+    EXPECT_EQ(natoms, static_cast<int>(expected_e_default.size()));
+
+    expected_tot_e_explicit = 0.;
+    expected_tot_e_default = 0.;
+    expected_tot_v_explicit.assign(9, 0.);
+    expected_tot_v_default.assign(9, 0.);
+    for (int ii = 0; ii < natoms; ++ii) {
+      expected_tot_e_explicit += expected_e_explicit[ii];
+      expected_tot_e_default += expected_e_default[ii];
+    }
+    for (int ii = 0; ii < natoms; ++ii) {
+      for (int dd = 0; dd < 9; ++dd) {
+        expected_tot_v_explicit[dd] += expected_v_explicit[ii * 9 + dd];
+        expected_tot_v_default[dd] += expected_v_default[ii * 9 + dd];
+      }
+    }
+  }
+
+  void TearDown() override {}
+
+  static void TearDownTestSuite() { dp = deepmd::DeepPot(); }
+};
+
+template <class VALUETYPE>
+deepmd::DeepPot TestInferDeepPotChgSpinPtExpt<VALUETYPE>::dp;
+
+TYPED_TEST_SUITE(TestInferDeepPotChgSpinPtExpt, ValueTypes);
+
+TYPED_TEST(TestInferDeepPotChgSpinPtExpt, dim_chg_spin) {
+  deepmd::DeepPot& dp = this->dp;
+  EXPECT_EQ(dp.dim_chg_spin(), 2);
+}
+
+TYPED_TEST(TestInferDeepPotChgSpinPtExpt, cpu_build_nlist_explicit) {
+  using VALUETYPE = TypeParam;
+  std::vector<VALUETYPE>& coord = this->coord;
+  std::vector<int>& atype = this->atype;
+  std::vector<VALUETYPE>& box = this->box;
+  std::vector<double>& charge_spin = this->charge_spin_explicit;
+  std::vector<VALUETYPE>& expected_f = this->expected_f_explicit;
+  std::vector<VALUETYPE>& expected_tot_v = this->expected_tot_v_explicit;
+  int& natoms = this->natoms;
+  double& expected_tot_e = this->expected_tot_e_explicit;
+  deepmd::DeepPot& dp = this->dp;
+
+  double ener;
+  std::vector<VALUETYPE> force, virial;
+  dp.compute(ener, force, virial, coord, atype, box, {}, {}, charge_spin);
+
+  EXPECT_EQ(force.size(), static_cast<size_t>(natoms * 3));
+  EXPECT_EQ(virial.size(), 9u);
+  EXPECT_LT(fabs(ener - expected_tot_e), EPSILON);
+  for (int ii = 0; ii < natoms * 3; ++ii) {
+    EXPECT_LT(fabs(force[ii] - expected_f[ii]), EPSILON);
+  }
+  for (int ii = 0; ii < 9; ++ii) {
+    EXPECT_LT(fabs(virial[ii] - expected_tot_v[ii]), EPSILON);
+  }
+}
+
+TYPED_TEST(TestInferDeepPotChgSpinPtExpt, cpu_build_nlist_default) {
+  using VALUETYPE = TypeParam;
+  std::vector<VALUETYPE>& coord = this->coord;
+  std::vector<int>& atype = this->atype;
+  std::vector<VALUETYPE>& box = this->box;
+  std::vector<VALUETYPE>& expected_f = this->expected_f_default;
+  std::vector<VALUETYPE>& expected_tot_v = this->expected_tot_v_default;
+  int& natoms = this->natoms;
+  double& expected_tot_e = this->expected_tot_e_default;
+  deepmd::DeepPot& dp = this->dp;
+
+  double ener;
+  std::vector<VALUETYPE> force, virial;
+  // Empty charge_spin triggers default_chg_spin_ fallback ([0.0, 1.0])
+  dp.compute(ener, force, virial, coord, atype, box);
+
+  EXPECT_EQ(force.size(), static_cast<size_t>(natoms * 3));
+  EXPECT_EQ(virial.size(), 9u);
+  EXPECT_LT(fabs(ener - expected_tot_e), EPSILON);
+  for (int ii = 0; ii < natoms * 3; ++ii) {
+    EXPECT_LT(fabs(force[ii] - expected_f[ii]), EPSILON);
+  }
+  for (int ii = 0; ii < 9; ++ii) {
+    EXPECT_LT(fabs(virial[ii] - expected_tot_v[ii]), EPSILON);
+  }
+}
+
+TYPED_TEST(TestInferDeepPotChgSpinPtExpt, cpu_lmp_nlist_explicit) {
+  using VALUETYPE = TypeParam;
+  std::vector<VALUETYPE>& coord = this->coord;
+  std::vector<int>& atype = this->atype;
+  std::vector<VALUETYPE>& box = this->box;
+  std::vector<double>& charge_spin = this->charge_spin_explicit;
+  std::vector<VALUETYPE>& expected_f = this->expected_f_explicit;
+  std::vector<VALUETYPE>& expected_tot_v = this->expected_tot_v_explicit;
+  int& natoms = this->natoms;
+  double& expected_tot_e = this->expected_tot_e_explicit;
+  deepmd::DeepPot& dp = this->dp;
+
+  float rc = dp.cutoff();
+  int nloc = coord.size() / 3;
+  std::vector<VALUETYPE> coord_cpy;
+  std::vector<int> atype_cpy, mapping;
+  std::vector<std::vector<int>> nlist_data;
+  _build_nlist<VALUETYPE>(nlist_data, coord_cpy, atype_cpy, mapping, coord,
+                          atype, box, rc);
+  int nall = coord_cpy.size() / 3;
+  std::vector<int> ilist(nloc), numneigh(nloc);
+  std::vector<int*> firstneigh(nloc);
+  deepmd::InputNlist inlist(nloc, &ilist[0], &numneigh[0], &firstneigh[0]);
+  convert_nlist(inlist, nlist_data);
+  // DPA3 repflows uses mapping; pass it so the C++ backend can
+  // correctly map ghost atoms back to local atoms.
+  inlist.mapping = mapping.data();
+
+  double ener;
+  std::vector<VALUETYPE> force_, virial;
+  dp.compute(ener, force_, virial, coord_cpy, atype_cpy, box, nall - nloc,
+             inlist, 0, {}, {}, charge_spin);
+  std::vector<VALUETYPE> force;
+  _fold_back<VALUETYPE>(force, force_, mapping, nloc, nall, 3);
+
+  EXPECT_EQ(force.size(), static_cast<size_t>(natoms * 3));
+  EXPECT_EQ(virial.size(), 9u);
+  EXPECT_LT(fabs(ener - expected_tot_e), EPSILON);
+  for (int ii = 0; ii < natoms * 3; ++ii) {
+    EXPECT_LT(fabs(force[ii] - expected_f[ii]), EPSILON);
+  }
+  for (int ii = 0; ii < 9; ++ii) {
+    EXPECT_LT(fabs(virial[ii] - expected_tot_v[ii]), EPSILON);
+  }
+
+  // Verify ago=1 (cached nlist) gives same result
+  ener = 0.;
+  std::fill(force_.begin(), force_.end(), 0.0);
+  std::fill(virial.begin(), virial.end(), 0.0);
+  dp.compute(ener, force_, virial, coord_cpy, atype_cpy, box, nall - nloc,
+             inlist, 1, {}, {}, charge_spin);
+  _fold_back<VALUETYPE>(force, force_, mapping, nloc, nall, 3);
+
+  EXPECT_LT(fabs(ener - expected_tot_e), EPSILON);
+  for (int ii = 0; ii < natoms * 3; ++ii) {
+    EXPECT_LT(fabs(force[ii] - expected_f[ii]), EPSILON);
+  }
+}
+
+TYPED_TEST(TestInferDeepPotChgSpinPtExpt, cpu_lmp_nlist_default) {
+  using VALUETYPE = TypeParam;
+  std::vector<VALUETYPE>& coord = this->coord;
+  std::vector<int>& atype = this->atype;
+  std::vector<VALUETYPE>& box = this->box;
+  std::vector<VALUETYPE>& expected_f = this->expected_f_default;
+  std::vector<VALUETYPE>& expected_tot_v = this->expected_tot_v_default;
+  int& natoms = this->natoms;
+  double& expected_tot_e = this->expected_tot_e_default;
+  deepmd::DeepPot& dp = this->dp;
+
+  float rc = dp.cutoff();
+  int nloc = coord.size() / 3;
+  std::vector<VALUETYPE> coord_cpy;
+  std::vector<int> atype_cpy, mapping;
+  std::vector<std::vector<int>> nlist_data;
+  _build_nlist<VALUETYPE>(nlist_data, coord_cpy, atype_cpy, mapping, coord,
+                          atype, box, rc);
+  int nall = coord_cpy.size() / 3;
+  std::vector<int> ilist(nloc), numneigh(nloc);
+  std::vector<int*> firstneigh(nloc);
+  deepmd::InputNlist inlist(nloc, &ilist[0], &numneigh[0], &firstneigh[0]);
+  convert_nlist(inlist, nlist_data);
+  // DPA3 repflows uses mapping; pass it so the C++ backend can
+  // correctly map ghost atoms back to local atoms.
+  inlist.mapping = mapping.data();
+
+  double ener;
+  std::vector<VALUETYPE> force_, virial;
+  // Empty charge_spin triggers default_chg_spin_ fallback ([0.0, 1.0])
+  dp.compute(ener, force_, virial, coord_cpy, atype_cpy, box, nall - nloc,
+             inlist, 0);
+  std::vector<VALUETYPE> force;
+  _fold_back<VALUETYPE>(force, force_, mapping, nloc, nall, 3);
+
+  EXPECT_EQ(force.size(), static_cast<size_t>(natoms * 3));
+  EXPECT_EQ(virial.size(), 9u);
+  EXPECT_LT(fabs(ener - expected_tot_e), EPSILON);
+  for (int ii = 0; ii < natoms * 3; ++ii) {
+    EXPECT_LT(fabs(force[ii] - expected_f[ii]), EPSILON);
+  }
+  for (int ii = 0; ii < 9; ++ii) {
+    EXPECT_LT(fabs(virial[ii] - expected_tot_v[ii]), EPSILON);
+  }
+}
diff --git a/source/install/test_cc_local.sh b/source/install/test_cc_local.sh
index 5ddbf0eecc..d5b2f3afb0 100755
--- a/source/install/test_cc_local.sh
+++ b/source/install/test_cc_local.sh
@@ -87,9 +87,12 @@ else:
 	PID5=$!
 	env ${_GEN_ENV} python ${INFER_SCRIPT_PATH}/gen_model_devi.py &
 	PID6=$!
+	env ${_GEN_ENV} python ${INFER_SCRIPT_PATH}/gen_chg_spin.py &
+	PID9=$!
 	wait $PID4
 	wait $PID5
 	wait $PID6
+	wait $PID9
 
 	env ${_GEN_ENV} python ${INFER_SCRIPT_PATH}/gen_spin.py &
 	PID7=$!
diff --git a/source/lmp/pair_base.cpp b/source/lmp/pair_base.cpp
index 63e5462590..faadd086b5 100644
--- a/source/lmp/pair_base.cpp
+++ b/source/lmp/pair_base.cpp
@@ -337,6 +337,9 @@ PairDeepBaseModel::PairDeepBaseModel(
   eps_v = 0.;
   scale = NULL;
   do_ttm = false;
+  dim_fparam = 0;
+  dim_aparam = 0;
+  dim_chg_spin = 0;
   do_compute_fparam = false;
   do_compute_aparam = false;
   single_model = false;
diff --git a/source/lmp/pair_base.h b/source/lmp/pair_base.h
index 1dd4b84041..93c5cb7b39 100644
--- a/source/lmp/pair_base.h
+++ b/source/lmp/pair_base.h
@@ -65,6 +65,7 @@ class PairDeepBaseModel : public Pair {
   std::string out_file;
   int dim_fparam;
   int dim_aparam;
+  int dim_chg_spin;
   int out_each;
   int out_rel;
   int out_rel_v;
@@ -80,6 +81,7 @@ class PairDeepBaseModel : public Pair {
   std::map<int, int> old_idx_map;
   std::vector<double> fparam;
   std::vector<double> aparam;
+  std::vector<double> charge_spin;
   double eps;
   double eps_v;
 
diff --git a/source/lmp/pair_deepmd.cpp b/source/lmp/pair_deepmd.cpp
index fb27fa0ff6..c1c4c29635 100644
--- a/source/lmp/pair_deepmd.cpp
+++ b/source/lmp/pair_deepmd.cpp
@@ -249,7 +249,7 @@ void PairDeepMD::compute(int eflag, int vflag) {
       if (!(eflag_atom || cvflag_atom)) {
         try {
           deep_pot.compute(dener, dforce, dvirial, dcoord, dtype, dbox, nghost,
-                           lmp_list, ago, fparam, daparam);
+                           lmp_list, ago, fparam, daparam, charge_spin);
         } catch (deepmd_compat::deepmd_exception& e) {
           error->one(FLERR, e.what());
         }
@@ -260,7 +260,8 @@ void PairDeepMD::compute(int eflag, int vflag) {
         vector<double> dvatom(nall * 9, 0);
         try {
           deep_pot.compute(dener, dforce, dvirial, deatom, dvatom, dcoord,
-                           dtype, dbox, nghost, lmp_list, ago, fparam, daparam);
+                           dtype, dbox, nghost, lmp_list, ago, fparam, daparam,
+                           charge_spin);
         } catch (deepmd_compat::deepmd_exception& e) {
           error->one(FLERR, e.what());
         }
@@ -312,7 +313,7 @@ void PairDeepMD::compute(int eflag, int vflag) {
         try {
           deep_pot_model_devi.compute(all_energy, all_force, all_virial, dcoord,
                                       dtype, dbox, nghost, lmp_list, ago,
-                                      fparam, daparam);
+                                      fparam, daparam, charge_spin);
         } catch (deepmd_compat::deepmd_exception& e) {
           error->one(FLERR, e.what());
         }
@@ -321,7 +322,7 @@ void PairDeepMD::compute(int eflag, int vflag) {
           deep_pot_model_devi.compute(all_energy, all_force, all_virial,
                                       all_atom_energy, all_atom_virial, dcoord,
                                       dtype, dbox, nghost, lmp_list, ago,
-                                      fparam, daparam);
+                                      fparam, daparam, charge_spin);
         } catch (deepmd_compat::deepmd_exception& e) {
           error->one(FLERR, e.what());
         }
@@ -534,6 +535,7 @@ static bool is_key(const string& input) {
   keys.push_back("aparam");
   keys.push_back("fparam_from_compute");
   keys.push_back("aparam_from_compute");
+  keys.push_back("charge_spin");
   keys.push_back("ttm");
   keys.push_back("atomic");
   keys.push_back("relative");
@@ -577,6 +579,7 @@ void PairDeepMD::settings(int narg, char** arg) {
     numb_types_spin = deep_pot.numb_types_spin();
     dim_fparam = deep_pot.dim_fparam();
     dim_aparam = deep_pot.dim_aparam();
+    dim_chg_spin = deep_pot.dim_chg_spin();
   } else {
     try {
       deep_pot.init(arg[0], get_node_rank(), get_file_content(arg[0]));
@@ -590,11 +593,13 @@ void PairDeepMD::settings(int narg, char** arg) {
     numb_types_spin = deep_pot_model_devi.numb_types_spin();
     dim_fparam = deep_pot_model_devi.dim_fparam();
     dim_aparam = deep_pot_model_devi.dim_aparam();
+    dim_chg_spin = deep_pot_model_devi.dim_chg_spin();
     assert(cutoff == deep_pot.cutoff() * dist_unit_cvt_factor);
     assert(numb_types == deep_pot.numb_types());
     assert(numb_types_spin == deep_pot.numb_types_spin());
     assert(dim_fparam == deep_pot.dim_fparam());
     assert(dim_aparam == deep_pot.dim_aparam());
+    assert(dim_chg_spin == deep_pot.dim_chg_spin());
   }
 
   out_freq = 100;
@@ -604,6 +609,7 @@ void PairDeepMD::settings(int narg, char** arg) {
   eps = 0.;
   fparam.clear();
   aparam.clear();
+  charge_spin.clear();
   while (iarg < narg) {
     if (!is_key(arg[iarg])) {
       error->all(FLERR,
@@ -686,6 +692,17 @@ void PairDeepMD::settings(int narg, char** arg) {
       do_compute_aparam = true;
       compute_aparam_id = arg[iarg + 1];
       iarg += 1 + 1;
+    } else if (string(arg[iarg]) == string("charge_spin")) {
+      for (int ii = 0; ii < dim_chg_spin; ++ii) {
+        if (iarg + 1 + ii >= narg || is_key(arg[iarg + 1 + ii])) {
+          char tmp[1024];
+          sprintf(tmp, "Illegal charge_spin, the dimension should be %d",
+                  dim_chg_spin);
+          error->all(FLERR, tmp);
+        }
+        charge_spin.push_back(atof(arg[iarg + 1 + ii]));
+      }
+      iarg += 1 + dim_chg_spin;
     } else if (string(arg[iarg]) == string("atomic")) {
       out_each = 1;
       iarg += 1;
diff --git a/source/lmp/tests/test_lammps_chg_spin_pt2.py b/source/lmp/tests/test_lammps_chg_spin_pt2.py
new file mode 100644
index 0000000000..a995b25061
--- /dev/null
+++ b/source/lmp/tests/test_lammps_chg_spin_pt2.py
@@ -0,0 +1,129 @@
+# SPDX-License-Identifier: LGPL-3.0-or-later
+"""Test the LAMMPS ``charge_spin`` pair_style keyword (.pt2 DPA3 model).
+
+Exercises the parsing and wiring of the ``charge_spin`` keyword through to the
+backend: the model (chg_spin.pt2, DPA3 with add_chg_spin_ebd=True and
+default_chg_spin=[0.0, 1.0]) is generated by source/tests/infer/gen_chg_spin.py,
+together with the reference values in chg_spin.expected.
+"""
+
+import os
+from pathlib import (
+    Path,
+)
+
+import numpy as np
+import pytest
+from expected_ref import (
+    read_expected_ref,
+)
+from lammps import (
+    PyLammps,
+)
+from write_lmp_data import (
+    write_lmp_data,
+)
+
+pt2_file = Path(__file__).parent.parent.parent / "tests" / "infer" / "chg_spin.pt2"
+ref_file = Path(__file__).parent.parent.parent / "tests" / "infer" / "chg_spin.expected"
+data_file = Path(__file__).parent / "data_chg_spin_pt2.lmp"
+
+# Reference values written by source/tests/infer/gen_chg_spin.py (PBC case):
+#   "default"  -- no charge_spin keyword -> stored default_chg_spin [0.0, 1.0]
+#   "explicit" -- charge_spin 1.0 2.0
+# Guarded with try/except because gen_chg_spin.py only runs when PyTorch is
+# built; matrices that disable PyTorch still load this file at collection time.
+try:
+    _ref = read_expected_ref(ref_file)
+    expected_e_default = float(np.sum(_ref["default"]["expected_e"]))
+    expected_f_default = _ref["default"]["expected_f"].reshape(6, 3)
+    expected_e_explicit = float(np.sum(_ref["explicit"]["expected_e"]))
+    expected_f_explicit = _ref["explicit"]["expected_f"].reshape(6, 3)
+except FileNotFoundError:
+    expected_e_default = expected_f_default = None
+    expected_e_explicit = expected_f_explicit = None
+
+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],
+    ]
+)
+# atype [0, 1, 1, 0, 1, 1] (O, H, H, O, H, H) -> LAMMPS types [1, 2, 2, 1, 2, 2]
+type_OH = np.array([1, 2, 2, 1, 2, 2])
+
+
+def setup_module() -> None:
+    if os.environ.get("ENABLE_PYTORCH", "1") != "1":
+        pytest.skip("Skip test because PyTorch support is not enabled.")
+    write_lmp_data(box, coord, type_OH, data_file)
+
+
+def teardown_module() -> None:
+    if data_file.exists():
+        os.remove(data_file)
+
+
+def _lammps(data_file, units="metal") -> PyLammps:
+    lammps = PyLammps()
+    lammps.units(units)
+    lammps.boundary("p p p")
+    lammps.atom_style("atomic")
+    # DPA3 (message passing) resolves ghost atoms via the LAMMPS atom-map.
+    lammps.atom_modify("map yes")
+    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.mass("2 2")
+    lammps.timestep(0.0005)
+    lammps.fix("1 all nve")
+    return lammps
+
+
+@pytest.fixture
+def lammps():
+    if not pt2_file.is_file():
+        pytest.skip(f"Skip because {pt2_file} was not generated.")
+    lmp = _lammps(data_file=data_file)
+    yield lmp
+    lmp.close()
+
+
+def test_pair_deepmd_charge_spin_default(lammps) -> None:
+    """No charge_spin keyword -> the model's stored default_chg_spin is used."""
+    lammps.pair_style(f"deepmd {pt2_file.resolve()}")
+    lammps.pair_coeff("* *")
+    lammps.run(0)
+    assert lammps.eval("pe") == pytest.approx(expected_e_default)
+    for ii in range(6):
+        assert lammps.atoms[ii].force == pytest.approx(
+            expected_f_default[lammps.atoms[ii].id - 1]
+        )
+    lammps.run(1)
+
+
+def test_pair_deepmd_charge_spin_explicit(lammps) -> None:
+    """Explicit ``charge_spin`` keyword is parsed and threaded to the model."""
+    lammps.pair_style(f"deepmd {pt2_file.resolve()} charge_spin 1.0 2.0")
+    lammps.pair_coeff("* *")
+    lammps.run(0)
+    assert lammps.eval("pe") == pytest.approx(expected_e_explicit)
+    for ii in range(6):
+        assert lammps.atoms[ii].force == pytest.approx(
+            expected_f_explicit[lammps.atoms[ii].id - 1]
+        )
+    lammps.run(1)
+
+
+def test_charge_spin_changes_result(lammps) -> None:
+    """Different charge_spin must give a different energy (keyword takes effect)."""
+    lammps.pair_style(f"deepmd {pt2_file.resolve()} charge_spin 1.0 2.0")
+    lammps.pair_coeff("* *")
+    lammps.run(0)
+    assert lammps.eval("pe") != pytest.approx(expected_e_default)
diff --git a/source/tests/infer/gen_chg_spin.py b/source/tests/infer/gen_chg_spin.py
new file mode 100644
index 0000000000..327d967af4
--- /dev/null
+++ b/source/tests/infer/gen_chg_spin.py
@@ -0,0 +1,240 @@
+#!/usr/bin/env python3
+# SPDX-License-Identifier: LGPL-3.0-or-later
+"""Generate chg_spin.pt2 / chg_spin.pth test models and reference values.
+
+Creates a DPA3 model with add_chg_spin_ebd=True and default_chg_spin=[0.0, 1.0],
+exports it to both .pt2 (AOTInductor) and .pth (TorchScript) from the same
+weights, and writes chg_spin.expected with two sections:
+  [default]  -- eval with no charge_spin (uses stored default [0.0, 1.0])
+  [explicit] -- eval with charge_spin=[1.0, 2.0]
+The .pth is verified to match the .pt2 reference for both sections.
+"""
+
+import copy
+import os
+import sys
+
+import numpy as np
+
+sys.path.insert(0, os.path.join(os.path.dirname(__file__), "..", "..", ".."))
+
+from gen_common import (
+    ensure_inductor_compiler,
+    load_custom_ops,
+    write_expected_ref,
+)
+
+
+def main():
+    from deepmd.dpmodel.model.model import (
+        get_model,
+    )
+
+    ensure_inductor_compiler()
+
+    config = {
+        "type_map": ["O", "H"],
+        "descriptor": {
+            "type": "dpa3",
+            "repflow": {
+                "n_dim": 8,
+                "e_dim": 5,
+                "a_dim": 4,
+                "nlayers": 2,
+                "e_rcut": 6.0,
+                "e_rcut_smth": 2.0,
+                "e_sel": 30,
+                "a_rcut": 4.0,
+                "a_rcut_smth": 2.0,
+                "a_sel": 20,
+                "axis_neuron": 4,
+                "update_angle": True,
+                "update_style": "res_residual",
+                "update_residual_init": "const",
+                "smooth_edge_update": True,
+            },
+            "concat_output_tebd": True,
+            "precision": "float64",
+            "add_chg_spin_ebd": True,
+            "default_chg_spin": [0.0, 1.0],
+            "seed": 1,
+        },
+        "fitting_net": {"neuron": [5, 5, 5], "resnet_dt": True, "seed": 1},
+    }
+
+    model = get_model(copy.deepcopy(config))
+    model_dict = model.serialize()
+
+    data = {
+        "model": model_dict,
+        "model_def_script": config,
+        "backend": "dpmodel",
+        "software": "deepmd-kit",
+        "version": "3.0.0",
+    }
+
+    from deepmd.pt.utils.serialization import (
+        deserialize_to_file as pt_deserialize_to_file,
+    )
+    from deepmd.pt_expt.utils.serialization import (
+        deserialize_to_file as pt_expt_deserialize_to_file,
+    )
+
+    load_custom_ops()
+
+    base_dir = os.path.dirname(__file__)
+    pt2_path = os.path.join(base_dir, "chg_spin.pt2")
+    print(f"Exporting to {pt2_path} ...")  # noqa: T201
+    pt_expt_deserialize_to_file(pt2_path, copy.deepcopy(data), do_atomic_virial=True)
+
+    pth_path = os.path.join(base_dir, "chg_spin.pth")
+    # Remove any stale .pth first so a failed export below cannot leave an old
+    # artifact that the parity check would then validate against.
+    if os.path.exists(pth_path):
+        os.remove(pth_path)
+    pth_generated = False
+    print(f"Exporting to {pth_path} ...")  # noqa: T201
+    try:
+        pt_deserialize_to_file(pth_path, copy.deepcopy(data))
+        pth_generated = True
+    except RuntimeError as e:
+        # Custom ops may not be available in all build environments; .pth
+        # generation is not critical (the .pth test skips if the file is
+        # missing / PyTorch support is off).
+        print(f"WARNING: .pth export failed ({e}), skipping.")  # noqa: T201
+    print("Export done.")  # noqa: T201
+
+    from deepmd.infer import (
+        DeepPot,
+    )
+
+    dp = DeepPot(pt2_path)
+
+    dim = dp.deep_eval.get_dim_chg_spin()
+    assert dim == 2, f"Expected dim_chg_spin == 2, got {dim}"
+    print(f"dim_chg_spin = {dim}")  # noqa: T201
+
+    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,
+        ],
+        dtype=np.float64,
+    )
+    atype = [0, 1, 1, 0, 1, 1]
+    box = np.array([13.0, 0.0, 0.0, 0.0, 13.0, 0.0, 0.0, 0.0, 13.0], dtype=np.float64)
+
+    # Default charge_spin: no argument → uses stored [0.0, 1.0]
+    e_def, f_def, v_def, ae_def, av_def = dp.eval(coord, box, atype, atomic=True)
+    print(f"\n// Default charge_spin total energy: {e_def[0, 0]:.18e}")  # noqa: T201
+
+    # Explicit charge_spin = [1.0, 2.0] (charge idx 101, spin idx 2 — both
+    # differ from the default [0.0, 1.0] which maps to charge idx 100, spin 1)
+    charge_spin_exp = np.array([[1.0, 2.0]], dtype=np.float64)
+    e_exp, f_exp, v_exp, ae_exp, av_exp = dp.eval(
+        coord, box, atype, atomic=True, charge_spin=charge_spin_exp
+    )
+    print(f"\n// Explicit charge_spin [1.0, 2.0] total energy: {e_exp[0, 0]:.18e}")  # noqa: T201
+
+    # Sanity: different charge_spin should yield different outputs
+    assert not np.allclose(e_def, e_exp), (
+        "Default and explicit charge_spin gave identical energy — charge_spin may be ignored"
+    )
+
+    # NoPbc variants (box=None) — used by the .pth nlist test, which mirrors the
+    # established DPA3 .pth pattern (NoPbc lmp_nlist, nghost=0) and so needs its
+    # own reference values.
+    e_np_def, f_np_def, v_np_def, ae_np_def, av_np_def = dp.eval(
+        coord, None, atype, atomic=True
+    )
+    e_np_exp, f_np_exp, v_np_exp, ae_np_exp, av_np_exp = dp.eval(
+        coord, None, atype, atomic=True, charge_spin=charge_spin_exp
+    )
+    assert not np.allclose(e_np_def, e_np_exp), (
+        "NoPbc: default and explicit charge_spin gave identical energy"
+    )
+
+    ref_path = os.path.join(base_dir, "chg_spin.expected")
+    write_expected_ref(
+        ref_path,
+        sections={
+            "default": {
+                "expected_e": ae_def[0, :, 0],
+                "expected_f": f_def[0],
+                "expected_v": av_def[0],
+            },
+            "explicit": {
+                "expected_e": ae_exp[0, :, 0],
+                "expected_f": f_exp[0],
+                "expected_v": av_exp[0],
+            },
+            "nopbc_default": {
+                "expected_e": ae_np_def[0, :, 0],
+                "expected_f": f_np_def[0],
+                "expected_v": av_np_def[0],
+            },
+            "nopbc_explicit": {
+                "expected_e": ae_np_exp[0, :, 0],
+                "expected_f": f_np_exp[0],
+                "expected_v": av_np_exp[0],
+            },
+        },
+        source_script="source/tests/infer/gen_chg_spin.py",
+    )
+    print(f"Wrote {ref_path}")  # noqa: T201
+
+    # ---- Verify .pth reproduces the .pt2 reference (PBC + NoPbc) ----
+    if pth_generated:
+        dp_pth = DeepPot(pth_path)
+        # Note: the .pth (pt) DeepEval does not expose get_dim_chg_spin (only
+        # the .pt2 / pt_expt one does); the energy/force parity below is the
+        # real check that charge_spin is threaded for the .pth backend.
+        tol = 1e-10
+        e_def_p, f_def_p, v_def_p = dp_pth.eval(coord, box, atype)
+        e_exp_p, f_exp_p, v_exp_p = dp_pth.eval(
+            coord, box, atype, charge_spin=charge_spin_exp
+        )
+        np.testing.assert_allclose(e_def_p, e_def, atol=tol, err_msg="pth default e")
+        np.testing.assert_allclose(f_def_p, f_def, atol=tol, err_msg="pth default f")
+        np.testing.assert_allclose(e_exp_p, e_exp, atol=tol, err_msg="pth explicit e")
+        np.testing.assert_allclose(f_exp_p, f_exp, atol=tol, err_msg="pth explicit f")
+        # NoPbc parity (the .pth nlist test uses NoPbc)
+        e_np_def_p, f_np_def_p, _ = dp_pth.eval(coord, None, atype)
+        e_np_exp_p, f_np_exp_p, _ = dp_pth.eval(
+            coord, None, atype, charge_spin=charge_spin_exp
+        )
+        np.testing.assert_allclose(
+            e_np_def_p, e_np_def, atol=tol, err_msg="pth nopbc default e"
+        )
+        np.testing.assert_allclose(
+            e_np_exp_p, e_np_exp, atol=tol, err_msg="pth nopbc explicit e"
+        )
+        # different charge_spin must change the .pth output too
+        assert not np.allclose(e_def_p, e_exp_p), (
+            ".pth: default and explicit charge_spin gave identical energy"
+        )
+        print("// .pth matches .pt2 reference (PBC + NoPbc).")  # noqa: T201
+    else:
+        print("\n// Skipping .pth verification (file not generated).")  # noqa: T201
+
+    print("\nDone!")  # noqa: T201
+
+
+if __name__ == "__main__":
+    main()
diff --git a/source/tests/pt_expt/infer/test_deep_eval.py b/source/tests/pt_expt/infer/test_deep_eval.py
index a921a94a7f..f831b365f6 100644
--- a/source/tests/pt_expt/infer/test_deep_eval.py
+++ b/source/tests/pt_expt/infer/test_deep_eval.py
@@ -2297,5 +2297,119 @@ def test_eval_descriptor_vesin_matches_native(self) -> None:
             np.testing.assert_allclose(d_native, d_vesin, rtol=1e-10, atol=1e-10)
 
 
+class TestDeepEvalEnerChgSpinPt2(unittest.TestCase):
+    """Test .pt2 inference with charge_spin (DPA3 with add_chg_spin_ebd=True)."""
+
+    @classmethod
+    def setUpClass(cls) -> None:
+        import copy
+
+        from deepmd.dpmodel.model.model import get_model as dp_get_model
+
+        cls.nt = 2
+        cls.default_chg_spin = [0.5, 0.8]
+
+        config = {
+            "type_map": ["O", "H"],
+            "descriptor": {
+                "type": "dpa3",
+                "repflow": {
+                    "n_dim": 8,
+                    "e_dim": 5,
+                    "a_dim": 4,
+                    "nlayers": 2,
+                    "e_rcut": 6.0,
+                    "e_rcut_smth": 2.0,
+                    "e_sel": 30,
+                    "a_rcut": 4.0,
+                    "a_rcut_smth": 2.0,
+                    "a_sel": 20,
+                    "axis_neuron": 4,
+                    "update_angle": True,
+                    "update_style": "res_residual",
+                    "update_residual_init": "const",
+                    "smooth_edge_update": True,
+                },
+                "concat_output_tebd": True,
+                "precision": "float64",
+                "add_chg_spin_ebd": True,
+                "default_chg_spin": cls.default_chg_spin,
+                "seed": GLOBAL_SEED,
+            },
+            "fitting_net": {
+                "neuron": [5, 5, 5],
+                "resnet_dt": True,
+                "seed": GLOBAL_SEED,
+            },
+        }
+
+        model = dp_get_model(copy.deepcopy(config))
+        cls.model_data = {
+            "model": model.serialize(),
+            "model_def_script": config,
+            "backend": "dpmodel",
+            "software": "deepmd-kit",
+            "version": "3.0.0",
+        }
+
+        cls.tmpfile = tempfile.NamedTemporaryFile(suffix=".pt2", delete=False)
+        cls.tmpfile.close()
+        torch.set_default_device(None)
+        try:
+            deserialize_to_file(cls.tmpfile.name, cls.model_data, do_atomic_virial=True)
+        finally:
+            torch.set_default_device("cuda:9999999")
+
+        cls.dp = DeepPot(cls.tmpfile.name)
+
+    @classmethod
+    def tearDownClass(cls) -> None:
+        import os
+
+        os.unlink(cls.tmpfile.name)
+
+    def test_dim_chg_spin(self) -> None:
+        self.assertEqual(self.dp.deep_eval.get_dim_chg_spin(), 2)
+
+    def test_charge_spin_takes_effect(self) -> None:
+        """Different charge_spin values must produce different outputs."""
+        rng = np.random.default_rng(GLOBAL_SEED)
+        natoms = 5
+        coords = rng.random((1, natoms, 3)) * 8.0
+        cells = np.eye(3).reshape(1, 9) * 10.0
+        atom_types = np.array([i % self.nt for i in range(natoms)], dtype=np.int32)
+
+        e0, f0, v0 = self.dp.eval(
+            coords, cells, atom_types, charge_spin=np.array([[0.0, 0.0]])
+        )
+        e1, f1, v1 = self.dp.eval(
+            coords, cells, atom_types, charge_spin=np.array([[1.0, 2.0]])
+        )
+
+        assert not np.allclose(e0, e1), (
+            "Changing charge_spin did not change output — charge_spin may be ignored"
+        )
+
+    def test_default_matches_explicit_default(self) -> None:
+        """Eval without charge_spin should use stored default and match explicit."""
+        rng = np.random.default_rng(GLOBAL_SEED)
+        natoms = 5
+        coords = rng.random((1, natoms, 3)) * 8.0
+        cells = np.eye(3).reshape(1, 9) * 10.0
+        atom_types = np.array([i % self.nt for i in range(natoms)], dtype=np.int32)
+
+        e_no, f_no, v_no = self.dp.eval(coords, cells, atom_types)
+        e_ex, f_ex, v_ex = self.dp.eval(
+            coords,
+            cells,
+            atom_types,
+            charge_spin=np.array([self.default_chg_spin]),
+        )
+
+        np.testing.assert_allclose(e_no, e_ex, atol=1e-10)
+        np.testing.assert_allclose(f_no, f_ex, atol=1e-10)
+        np.testing.assert_allclose(v_no, v_ex, atol=1e-10)
+
+
 if __name__ == "__main__":
     unittest.main()