diff --git a/source/api_c/include/c_api.h b/source/api_c/include/c_api.h
index 911813e428..cac6de377a 100644
--- a/source/api_c/include/c_api.h
+++ b/source/api_c/include/c_api.h
@@ -512,6 +512,147 @@ extern DP_DeepPotModelDevi* DP_NewDeepPotModelDeviWithParam(
  */
 extern void DP_DeleteDeepPotModelDevi(DP_DeepPotModelDevi* dp);
 
+/**
+ * @brief Evaluate the energy, force and virial by using a DP model deviation
+ *with neighbor list. (double version)
+ * @param[in] dp The DP model deviation to use.
+ * @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] box The cell of the region. The array should be of size 9. Pass
+ *NULL if pbc is not used.
+ * @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.
+ **/
+extern void DP_DeepPotModelDeviCompute(DP_DeepPotModelDevi* dp,
+                                       const int natom,
+                                       const double* coord,
+                                       const int* atype,
+                                       const double* cell,
+                                       double* energy,
+                                       double* force,
+                                       double* virial,
+                                       double* atomic_energy,
+                                       double* atomic_virial);
+
+/**
+ * @brief Evaluate the energy, force and virial by using a DP model deviation
+ *with neighbor list. (float version)
+ * @param[in] dp The DP model deviation to use.
+ * @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] box The cell of the region. The array should be of size 9. Pass
+ *NULL if pbc is not used.
+ * @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.
+ **/
+extern void DP_DeepPotModelDeviComputef(DP_DeepPotModelDevi* dp,
+                                        const int natom,
+                                        const float* coord,
+                                        const int* atype,
+                                        const float* cell,
+                                        double* energy,
+                                        float* force,
+                                        float* virial,
+                                        float* atomic_energy,
+                                        float* atomic_virial);
+
+/**
+ * @brief Evaluate the energy, force and virial by using a DP model deviation
+ *with neighbor list. (double version)
+ * @version 2
+ * @param[in] dp The DP model deviation to use.
+ * @param[in] nframes The number of frames. Only support 1 for now.
+ * @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] box 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[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.
+ **/
+void DP_DeepPotModelDeviCompute2(DP_DeepPotModelDevi* dp,
+                                 const int nframes,
+                                 const int natoms,
+                                 const double* coord,
+                                 const int* atype,
+                                 const double* cell,
+                                 const double* fparam,
+                                 const double* aparam,
+                                 double* energy,
+                                 double* force,
+                                 double* virial,
+                                 double* atomic_energy,
+                                 double* atomic_virial);
+/**
+ * @brief Evaluate the energy, force and virial by using a DP model deviation
+ *with neighbor list. (float version)
+ * @version 2
+ * @param[in] dp The DP model deviation to use.
+ * @param[in] nframes The number of frames. Only support 1 for now.
+ * @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] box 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[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.
+ **/
+void DP_DeepPotModelDeviComputef2(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);
+
 /**
  * @brief Evaluate the energy, force and virial by using a DP model deviation
  *with neighbor list. (double version)
diff --git a/source/api_c/include/deepmd.hpp b/source/api_c/include/deepmd.hpp
index 16b8f08cad..f348a3dbe9 100644
--- a/source/api_c/include/deepmd.hpp
+++ b/source/api_c/include/deepmd.hpp
@@ -210,6 +210,56 @@ inline void _DP_DeepPotComputeMixedType<float>(DP_DeepPot *dp,
                               atomic_virial);
 }
 
+template <typename FPTYPE>
+inline void _DP_DeepPotModelDeviCompute(DP_DeepPotModelDevi *dp,
+                                        const int natom,
+                                        const FPTYPE *coord,
+                                        const int *atype,
+                                        const FPTYPE *cell,
+                                        const FPTYPE *fparam,
+                                        const FPTYPE *aparam,
+                                        double *energy,
+                                        FPTYPE *force,
+                                        FPTYPE *virial,
+                                        FPTYPE *atomic_energy,
+                                        FPTYPE *atomic_virial);
+
+template <>
+inline void _DP_DeepPotModelDeviCompute<double>(DP_DeepPotModelDevi *dp,
+                                                const int natom,
+                                                const double *coord,
+                                                const int *atype,
+                                                const double *cell,
+                                                const double *fparam,
+                                                const double *aparam,
+                                                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);
+}
+
+template <>
+inline void _DP_DeepPotModelDeviCompute<float>(DP_DeepPotModelDevi *dp,
+                                               const int natom,
+                                               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) {
+  DP_DeepPotModelDeviComputef2(dp, 1, natom, coord, atype, cell, fparam, aparam,
+                               energy, force, virial, atomic_energy,
+                               atomic_virial);
+}
+
 template <typename FPTYPE>
 inline void _DP_DeepPotModelDeviComputeNList(DP_DeepPotModelDevi *dp,
                                              const int natom,
@@ -1179,6 +1229,189 @@ class DeepPotModelDevi {
    * @param[in] atype The atom types. The list should contain natoms ints.
    * @param[in] box The cell of the region. The array should be of size nframes
    *x 9 (PBC) or empty (no PBC).
+   * @param[in] fparam The frame parameter. The array can be of size :
+   * nframes x dim_fparam.
+   * dim_fparam. Then all frames are assumed to be provided with the same
+   *fparam.
+   * @param[in] aparam The atomic parameter The array can be of size :
+   * nframes x natoms x dim_aparam.
+   * natoms x dim_aparam. Then all frames are assumed to be provided with the
+   *same aparam.
+   **/
+  template <typename VALUETYPE>
+  void compute(
+      std::vector<double> &ener,
+      std::vector<std::vector<VALUETYPE>> &force,
+      std::vector<std::vector<VALUETYPE>> &virial,
+      const std::vector<VALUETYPE> &coord,
+      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>()) {
+    unsigned int natoms = atype.size();
+    unsigned int nframes = 1;
+    assert(natoms * 3 == coord.size());
+    if (!box.empty()) {
+      assert(box.size() == 9);
+    }
+    const VALUETYPE *coord_ = &coord[0];
+    const VALUETYPE *box_ = !box.empty() ? &box[0] : nullptr;
+    const int *atype_ = &atype[0];
+
+    // memory will be continous for std::vector but not std::vector<std::vector>
+    std::vector<double> energy_flat(numb_models);
+    std::vector<VALUETYPE> force_flat(static_cast<size_t>(numb_models) *
+                                      natoms * 3);
+    std::vector<VALUETYPE> virial_flat(numb_models * 9);
+    double *ener_ = &energy_flat[0];
+    VALUETYPE *force_ = &force_flat[0];
+    VALUETYPE *virial_ = &virial_flat[0];
+    std::vector<VALUETYPE> fparam_, aparam_;
+    validate_fparam_aparam(nframes, natoms, fparam, aparam);
+    tile_fparam_aparam(fparam_, nframes, dfparam, fparam);
+    tile_fparam_aparam(aparam_, nframes, natoms * daparam, aparam);
+    const VALUETYPE *fparam__ = !fparam_.empty() ? &fparam_[0] : nullptr;
+    const VALUETYPE *aparam__ = !aparam_.empty() ? &aparam_[0] : nullptr;
+
+    _DP_DeepPotModelDeviCompute<VALUETYPE>(dp, natoms, coord_, atype_, box_,
+                                           fparam__, aparam__, ener_, force_,
+                                           virial_, nullptr, nullptr);
+    DP_CHECK_OK(DP_DeepPotModelDeviCheckOK, dp);
+
+    // reshape
+    ener.resize(numb_models);
+    force.resize(numb_models);
+    virial.resize(numb_models);
+    for (int i = 0; i < numb_models; i++) {
+      ener[i] = energy_flat[i];
+      force[i].resize(static_cast<size_t>(natoms) * 3);
+      virial[i].resize(9);
+      for (int j = 0; j < natoms * 3; j++) {
+        force[i][j] = force_flat[i * natoms * 3 + j];
+      }
+      for (int j = 0; j < 9; j++) {
+        virial[i][j] = virial_flat[i * 9 + j];
+      }
+    }
+  };
+  /**
+   * @brief Evaluate the energy, force, virial, atomic energy, and atomic virial
+   *by using this DP model deviation.
+   * @param[out] ener The system energy.
+   * @param[out] force The force on each atom.
+   * @param[out] virial The virial.
+   * @param[out] atom_energy The atomic energy.
+   * @param[out] atom_virial The atomic virial.
+   * @param[in] coord The coordinates of atoms. The array should be of size
+   *nframes x natoms x 3.
+   * @param[in] atype The atom types. The list should contain natoms ints.
+   * @param[in] box The cell of the region. The array should be of size nframes
+   *x 9 (PBC) or empty (no PBC).
+   * @param[in] fparam The frame parameter. The array can be of size :
+   * nframes x dim_fparam.
+   * dim_fparam. Then all frames are assumed to be provided with the same
+   *fparam.
+   * @param[in] aparam The atomic parameter The array can be of size :
+   * nframes x natoms x dim_aparam.
+   * natoms x dim_aparam. Then all frames are assumed to be provided with the
+   *same aparam.
+   **/
+  template <typename VALUETYPE>
+  void compute(
+      std::vector<double> &ener,
+      std::vector<std::vector<VALUETYPE>> &force,
+      std::vector<std::vector<VALUETYPE>> &virial,
+      std::vector<std::vector<VALUETYPE>> &atom_energy,
+      std::vector<std::vector<VALUETYPE>> &atom_virial,
+      const std::vector<VALUETYPE> &coord,
+      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>()) {
+    unsigned int natoms = atype.size();
+    unsigned int nframes = 1;
+    assert(natoms * 3 == coord.size());
+    if (!box.empty()) {
+      assert(box.size() == 9);
+    }
+    const VALUETYPE *coord_ = &coord[0];
+    const VALUETYPE *box_ = !box.empty() ? &box[0] : nullptr;
+    const int *atype_ = &atype[0];
+
+    std::vector<double> energy_flat(numb_models);
+    std::vector<VALUETYPE> force_flat(static_cast<size_t>(numb_models) *
+                                      natoms * 3);
+    std::vector<VALUETYPE> virial_flat(numb_models * 9);
+    std::vector<VALUETYPE> atom_energy_flat(static_cast<size_t>(numb_models) *
+                                            natoms);
+    std::vector<VALUETYPE> atom_virial_flat(static_cast<size_t>(numb_models) *
+                                            natoms * 9);
+    double *ener_ = &energy_flat[0];
+    VALUETYPE *force_ = &force_flat[0];
+    VALUETYPE *virial_ = &virial_flat[0];
+    VALUETYPE *atomic_ener_ = &atom_energy_flat[0];
+    VALUETYPE *atomic_virial_ = &atom_virial_flat[0];
+    std::vector<VALUETYPE> fparam_, aparam_;
+    validate_fparam_aparam(nframes, natoms, fparam, aparam);
+    tile_fparam_aparam(fparam_, nframes, dfparam, fparam);
+    tile_fparam_aparam(aparam_, nframes, natoms * daparam, aparam);
+    const VALUETYPE *fparam__ = !fparam_.empty() ? &fparam_[0] : nullptr;
+    const VALUETYPE *aparam__ = !aparam_.empty() ? &aparam_[0] : nullptr;
+
+    _DP_DeepPotModelDeviCompute<VALUETYPE>(
+        dp, natoms, coord_, atype_, box_, fparam__, aparam__, ener_, force_,
+        virial_, atomic_ener_, atomic_virial_);
+    DP_CHECK_OK(DP_DeepPotModelDeviCheckOK, dp);
+
+    // reshape
+    ener.resize(numb_models);
+    force.resize(numb_models);
+    virial.resize(numb_models);
+    atom_energy.resize(numb_models);
+    atom_virial.resize(numb_models);
+    for (int i = 0; i < numb_models; i++) {
+      ener[i] = energy_flat[i];
+      force[i].resize(static_cast<size_t>(natoms) * 3);
+      virial[i].resize(9);
+      atom_energy[i].resize(natoms);
+      atom_virial[i].resize(static_cast<size_t>(natoms) * 9);
+      for (int j = 0; j < natoms * 3; j++) {
+        force[i][j] = force_flat[i * natoms * 3 + j];
+      }
+      for (int j = 0; j < 9; j++) {
+        virial[i][j] = virial_flat[i * 9 + j];
+      }
+      for (int j = 0; j < natoms; j++) {
+        atom_energy[i][j] = atom_energy_flat[i * natoms + j];
+      }
+      for (int j = 0; j < natoms * 9; j++) {
+        atom_virial[i][j] = atom_virial_flat[i * natoms * 9 + j];
+      }
+    }
+  };
+
+  /**
+   * @brief Evaluate the energy, force and virial by using this DP model
+   *deviation.
+   * @param[out] ener The system energy.
+   * @param[out] force The force on each atom.
+   * @param[out] virial The virial.
+   * @param[in] coord The coordinates of atoms. The array should be of size
+   *nframes x natoms x 3.
+   * @param[in] atype The atom types. The list should contain natoms ints.
+   * @param[in] box The cell of the region. The array should be of size nframes
+   *x 9 (PBC) or empty (no PBC).
+   * @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 parameter. The array can be of size :
+   * nframes x dim_fparam.
+   * dim_fparam. Then all frames are assumed to be provided with the same
+   *fparam.
+   * @param[in] aparam The atomic parameter The array can be of size :
+   * nframes x natoms x dim_aparam.
+   * natoms x dim_aparam. Then all frames are assumed to be provided with the
+   *same aparam.
    **/
   template <typename VALUETYPE>
   void compute(
@@ -1255,6 +1488,17 @@ class DeepPotModelDevi {
    * @param[in] atype The atom types. The list should contain natoms ints.
    * @param[in] box The cell of the region. The array should be of size nframes
    *x 9 (PBC) or empty (no PBC).
+   * @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 parameter. The array can be of size :
+   * nframes x dim_fparam.
+   * dim_fparam. Then all frames are assumed to be provided with the same
+   *fparam.
+   * @param[in] aparam The atomic parameter The array can be of size :
+   * nframes x natoms x dim_aparam.
+   * natoms x dim_aparam. Then all frames are assumed to be provided with the
+   *same aparam.
    **/
   template <typename VALUETYPE>
   void compute(
diff --git a/source/api_c/src/c_api.cc b/source/api_c/src/c_api.cc
index 79dc486e0d..e21cd48ffa 100644
--- a/source/api_c/src/c_api.cc
+++ b/source/api_c/src/c_api.cc
@@ -414,6 +414,100 @@ 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) {
+  if (nframes > 1) {
+    throw std::runtime_error("nframes > 1 not supported yet");
+  }
+  // init C++ vectors from C arrays
+  std::vector<VALUETYPE> coord_(coord, coord + natoms * 3);
+  std::vector<int> atype_(atype, atype + natoms);
+  std::vector<VALUETYPE> cell_;
+  if (cell) {
+    // pbc
+    cell_.assign(cell, cell + 9);
+  }
+  std::vector<VALUETYPE> fparam_;
+  if (fparam) {
+    fparam_.assign(fparam, fparam + dp->dfparam);
+  }
+  std::vector<VALUETYPE> aparam_;
+  if (aparam) {
+    aparam_.assign(aparam, aparam + nframes * natoms * dp->daparam);
+  }
+  // different from DeepPot
+  std::vector<double> e;
+  std::vector<std::vector<VALUETYPE>> f, v, ae, av;
+
+  DP_REQUIRES_OK(dp, dp->dp.compute(e, f, v, ae, av, coord_, atype_, cell_,
+                                    fparam_, aparam_));
+  // 2D vector to 2D array, flatten first
+  if (energy) {
+    std::copy(e.begin(), e.end(), energy);
+  }
+  if (force) {
+    std::vector<VALUETYPE> f_flat;
+    flatten_vector(f_flat, f);
+    std::copy(f_flat.begin(), f_flat.end(), force);
+  }
+  if (virial) {
+    std::vector<VALUETYPE> v_flat;
+    flatten_vector(v_flat, v);
+    std::copy(v_flat.begin(), v_flat.end(), virial);
+  }
+  if (atomic_energy) {
+    std::vector<VALUETYPE> ae_flat;
+    flatten_vector(ae_flat, ae);
+    std::copy(ae_flat.begin(), ae_flat.end(), atomic_energy);
+  }
+  if (atomic_virial) {
+    std::vector<VALUETYPE> av_flat;
+    flatten_vector(av_flat, av);
+    std::copy(av_flat.begin(), av_flat.end(), atomic_virial);
+  }
+}
+
+template void DP_DeepPotModelDeviCompute_variant<double>(
+    DP_DeepPotModelDevi* dp,
+    const int nframes,
+    const int natoms,
+    const double* coord,
+    const int* atype,
+    const double* cell,
+    const double* fparam,
+    const double* aparam,
+    double* energy,
+    double* force,
+    double* virial,
+    double* atomic_energy,
+    double* 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);
+
 template <typename VALUETYPE>
 void DP_DeepPotModelDeviComputeNList_variant(DP_DeepPotModelDevi* dp,
                                              const int nframes,
@@ -1054,6 +1148,72 @@ const char* DP_DeepPotCheckOK(DP_DeepPot* dp) {
   return string_to_char(dp->exception);
 }
 
+void DP_DeepPotModelDeviCompute(DP_DeepPotModelDevi* dp,
+                                const int natoms,
+                                const double* coord,
+                                const int* atype,
+                                const double* cell,
+                                double* energy,
+                                double* force,
+                                double* virial,
+                                double* atomic_energy,
+                                double* atomic_virial) {
+  DP_DeepPotModelDeviCompute_variant<double>(dp, 1, natoms, coord, atype, cell,
+                                             NULL, NULL, energy, force, virial,
+                                             atomic_energy, atomic_virial);
+}
+
+void DP_DeepPotModelDeviComputef(DP_DeepPotModelDevi* dp,
+                                 const int natoms,
+                                 const float* coord,
+                                 const int* atype,
+                                 const float* cell,
+                                 double* energy,
+                                 float* force,
+                                 float* virial,
+                                 float* atomic_energy,
+                                 float* atomic_virial) {
+  DP_DeepPotModelDeviCompute_variant<float>(dp, 1, natoms, coord, atype, cell,
+                                            NULL, NULL, energy, force, virial,
+                                            atomic_energy, atomic_virial);
+}
+
+void DP_DeepPotModelDeviCompute2(DP_DeepPotModelDevi* dp,
+                                 const int nframes,
+                                 const int natoms,
+                                 const double* coord,
+                                 const int* atype,
+                                 const double* cell,
+                                 const double* fparam,
+                                 const double* aparam,
+                                 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);
+}
+
+void DP_DeepPotModelDeviComputef2(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) {
+  DP_DeepPotModelDeviCompute_variant<float>(
+      dp, nframes, natoms, coord, atype, cell, fparam, aparam, energy, force,
+      virial, atomic_energy, atomic_virial);
+}
+
 void DP_DeepPotModelDeviComputeNList(DP_DeepPotModelDevi* dp,
                                      const int natoms,
                                      const double* coord,
diff --git a/source/api_c/tests/test_deeppot_model_devi_fparam_aparam.cc b/source/api_c/tests/test_deeppot_model_devi_fparam_aparam.cc
index 39bc1140d0..5e77529021 100644
--- a/source/api_c/tests/test_deeppot_model_devi_fparam_aparam.cc
+++ b/source/api_c/tests/test_deeppot_model_devi_fparam_aparam.cc
@@ -18,7 +18,7 @@ class TestInferDeepPotModeDeviFparamAparam : public ::testing::Test {
   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<int> atype = {0, 0, 0, 0, 0, 0};
   std::vector<VALUETYPE> box = {13., 0., 0., 0., 13., 0., 0., 0., 13.};
   std::vector<VALUETYPE> fparam = {0.25852028};
   std::vector<VALUETYPE> aparam = {0.25852028, 0.25852028, 0.25852028,
@@ -75,6 +75,96 @@ TYPED_TEST(TestInferDeepPotModeDeviFparamAparam, attrs) {
   EXPECT_EQ(dp1.dim_aparam(), dp_md.dim_aparam());
 }
 
+TYPED_TEST(TestInferDeepPotModeDeviFparamAparam, cpu_build_nlist) {
+  using VALUETYPE = TypeParam;
+  std::vector<VALUETYPE>& coord = this->coord;
+  std::vector<int>& atype = this->atype;
+  std::vector<VALUETYPE>& box = this->box;
+  std::vector<VALUETYPE>& fparam = this->fparam;
+  std::vector<VALUETYPE>& aparam = this->aparam;
+  int& natoms = this->natoms;
+  deepmd::hpp::DeepPot& dp0 = this->dp0;
+  deepmd::hpp::DeepPot& dp1 = this->dp1;
+  deepmd::hpp::DeepPotModelDevi& dp_md = this->dp_md;
+  float rc = dp_md.cutoff();
+  int nloc = coord.size() / 3;
+
+  int nmodel = 2;
+  std::vector<double> edir(nmodel), emd;
+  std::vector<std::vector<VALUETYPE> > fdir(nmodel), vdir(nmodel), fmd(nmodel),
+      vmd;
+  dp0.compute(edir[0], fdir[0], vdir[0], coord, atype, box, fparam, aparam);
+  dp1.compute(edir[1], fdir[1], vdir[1], coord, atype, box, fparam, aparam);
+  dp_md.compute(emd, fmd, vmd, coord, atype, box, fparam, aparam);
+
+  EXPECT_EQ(edir.size(), emd.size());
+  EXPECT_EQ(fdir.size(), fmd.size());
+  EXPECT_EQ(vdir.size(), vmd.size());
+  for (int kk = 0; kk < nmodel; ++kk) {
+    EXPECT_EQ(fdir[kk].size(), fmd[kk].size());
+    EXPECT_EQ(vdir[kk].size(), vmd[kk].size());
+  }
+  for (int kk = 0; kk < nmodel; ++kk) {
+    EXPECT_LT(fabs(edir[kk] - emd[kk]), EPSILON);
+    for (int ii = 0; ii < fdir[0].size(); ++ii) {
+      EXPECT_LT(fabs(fdir[kk][ii] - fmd[kk][ii]), EPSILON);
+    }
+    for (int ii = 0; ii < vdir[0].size(); ++ii) {
+      EXPECT_LT(fabs(vdir[kk][ii] - vmd[kk][ii]), EPSILON);
+    }
+  }
+}
+
+TYPED_TEST(TestInferDeepPotModeDeviFparamAparam, cpu_build_nlist_atomic) {
+  using VALUETYPE = TypeParam;
+  std::vector<VALUETYPE>& coord = this->coord;
+  std::vector<int>& atype = this->atype;
+  std::vector<VALUETYPE>& box = this->box;
+  std::vector<VALUETYPE>& fparam = this->fparam;
+  std::vector<VALUETYPE>& aparam = this->aparam;
+  int& natoms = this->natoms;
+  deepmd::hpp::DeepPot& dp0 = this->dp0;
+  deepmd::hpp::DeepPot& dp1 = this->dp1;
+  deepmd::hpp::DeepPotModelDevi& dp_md = this->dp_md;
+
+  int nmodel = 2;
+  std::vector<double> edir(nmodel), emd;
+  std::vector<std::vector<VALUETYPE> > fdir(nmodel), vdir(nmodel), fmd(nmodel),
+      vmd, aedir(nmodel), aemd, avdir(nmodel), avmd(nmodel);
+  dp0.compute(edir[0], fdir[0], vdir[0], aedir[0], avdir[0], coord, atype, box,
+              fparam, aparam);
+  dp1.compute(edir[1], fdir[1], vdir[1], aedir[1], avdir[1], coord, atype, box,
+              fparam, aparam);
+  dp_md.compute(emd, fmd, vmd, aemd, avmd, coord, atype, box, fparam, aparam);
+
+  EXPECT_EQ(edir.size(), emd.size());
+  EXPECT_EQ(fdir.size(), fmd.size());
+  EXPECT_EQ(vdir.size(), vmd.size());
+  EXPECT_EQ(aedir.size(), aemd.size());
+  EXPECT_EQ(avdir.size(), avmd.size());
+  for (int kk = 0; kk < nmodel; ++kk) {
+    EXPECT_EQ(fdir[kk].size(), fmd[kk].size());
+    EXPECT_EQ(vdir[kk].size(), vmd[kk].size());
+    EXPECT_EQ(aedir[kk].size(), aemd[kk].size());
+    EXPECT_EQ(avdir[kk].size(), avmd[kk].size());
+  }
+  for (int kk = 0; kk < nmodel; ++kk) {
+    EXPECT_LT(fabs(edir[kk] - emd[kk]), EPSILON);
+    for (int ii = 0; ii < fdir[0].size(); ++ii) {
+      EXPECT_LT(fabs(fdir[kk][ii] - fmd[kk][ii]), EPSILON);
+    }
+    for (int ii = 0; ii < vdir[0].size(); ++ii) {
+      EXPECT_LT(fabs(vdir[kk][ii] - vmd[kk][ii]), EPSILON);
+    }
+    for (int ii = 0; ii < aedir[0].size(); ++ii) {
+      EXPECT_LT(fabs(aedir[kk][ii] - aemd[kk][ii]), EPSILON);
+    }
+    for (int ii = 0; ii < avdir[0].size(); ++ii) {
+      EXPECT_LT(fabs(avdir[kk][ii] - avmd[kk][ii]), EPSILON);
+    }
+  }
+}
+
 TYPED_TEST(TestInferDeepPotModeDeviFparamAparam, cpu_lmp_list) {
   using VALUETYPE = TypeParam;
   std::vector<VALUETYPE>& coord = this->coord;
diff --git a/source/api_c/tests/test_deeppot_model_devi_hpp.cc b/source/api_c/tests/test_deeppot_model_devi_hpp.cc
index 75b0fab880..5751c4f69c 100644
--- a/source/api_c/tests/test_deeppot_model_devi_hpp.cc
+++ b/source/api_c/tests/test_deeppot_model_devi_hpp.cc
@@ -116,6 +116,90 @@ TYPED_TEST(TestInferDeepPotModeDevi, attrs) {
   // EXPECT_EQ(dp1.dim_aparam(), dp_md.dim_aparam());
 }
 
+TYPED_TEST(TestInferDeepPotModeDevi, cpu_build_nlist) {
+  using VALUETYPE = TypeParam;
+  std::vector<VALUETYPE>& coord = this->coord;
+  std::vector<int>& atype = this->atype;
+  std::vector<VALUETYPE>& box = this->box;
+  int& natoms = this->natoms;
+  deepmd::hpp::DeepPot& dp0 = this->dp0;
+  deepmd::hpp::DeepPot& dp1 = this->dp1;
+  deepmd::hpp::DeepPotModelDevi& dp_md = this->dp_md;
+  float rc = dp_md.cutoff();
+  int nloc = coord.size() / 3;
+
+  int nmodel = 2;
+  std::vector<double> edir(nmodel), emd;
+  std::vector<std::vector<VALUETYPE> > fdir(nmodel), vdir(nmodel), fmd(nmodel),
+      vmd;
+  dp0.compute(edir[0], fdir[0], vdir[0], coord, atype, box);
+  dp1.compute(edir[1], fdir[1], vdir[1], coord, atype, box);
+  dp_md.compute(emd, fmd, vmd, coord, atype, box);
+
+  EXPECT_EQ(edir.size(), emd.size());
+  EXPECT_EQ(fdir.size(), fmd.size());
+  EXPECT_EQ(vdir.size(), vmd.size());
+  for (int kk = 0; kk < nmodel; ++kk) {
+    EXPECT_EQ(fdir[kk].size(), fmd[kk].size());
+    EXPECT_EQ(vdir[kk].size(), vmd[kk].size());
+  }
+  for (int kk = 0; kk < nmodel; ++kk) {
+    EXPECT_LT(fabs(edir[kk] - emd[kk]), EPSILON);
+    for (int ii = 0; ii < fdir[0].size(); ++ii) {
+      EXPECT_LT(fabs(fdir[kk][ii] - fmd[kk][ii]), EPSILON);
+    }
+    for (int ii = 0; ii < vdir[0].size(); ++ii) {
+      EXPECT_LT(fabs(vdir[kk][ii] - vmd[kk][ii]), EPSILON);
+    }
+  }
+}
+
+TYPED_TEST(TestInferDeepPotModeDevi, cpu_build_nlist_atomic) {
+  using VALUETYPE = TypeParam;
+  std::vector<VALUETYPE>& coord = this->coord;
+  std::vector<int>& atype = this->atype;
+  std::vector<VALUETYPE>& box = this->box;
+  int& natoms = this->natoms;
+  deepmd::hpp::DeepPot& dp0 = this->dp0;
+  deepmd::hpp::DeepPot& dp1 = this->dp1;
+  deepmd::hpp::DeepPotModelDevi& dp_md = this->dp_md;
+
+  int nmodel = 2;
+  std::vector<double> edir(nmodel), emd;
+  std::vector<std::vector<VALUETYPE> > fdir(nmodel), vdir(nmodel), fmd(nmodel),
+      vmd, aedir(nmodel), aemd, avdir(nmodel), avmd(nmodel);
+  dp0.compute(edir[0], fdir[0], vdir[0], aedir[0], avdir[0], coord, atype, box);
+  dp1.compute(edir[1], fdir[1], vdir[1], aedir[1], avdir[1], coord, atype, box);
+  dp_md.compute(emd, fmd, vmd, aemd, avmd, coord, atype, box);
+
+  EXPECT_EQ(edir.size(), emd.size());
+  EXPECT_EQ(fdir.size(), fmd.size());
+  EXPECT_EQ(vdir.size(), vmd.size());
+  EXPECT_EQ(aedir.size(), aemd.size());
+  EXPECT_EQ(avdir.size(), avmd.size());
+  for (int kk = 0; kk < nmodel; ++kk) {
+    EXPECT_EQ(fdir[kk].size(), fmd[kk].size());
+    EXPECT_EQ(vdir[kk].size(), vmd[kk].size());
+    EXPECT_EQ(aedir[kk].size(), aemd[kk].size());
+    EXPECT_EQ(avdir[kk].size(), avmd[kk].size());
+  }
+  for (int kk = 0; kk < nmodel; ++kk) {
+    EXPECT_LT(fabs(edir[kk] - emd[kk]), EPSILON);
+    for (int ii = 0; ii < fdir[0].size(); ++ii) {
+      EXPECT_LT(fabs(fdir[kk][ii] - fmd[kk][ii]), EPSILON);
+    }
+    for (int ii = 0; ii < vdir[0].size(); ++ii) {
+      EXPECT_LT(fabs(vdir[kk][ii] - vmd[kk][ii]), EPSILON);
+    }
+    for (int ii = 0; ii < aedir[0].size(); ++ii) {
+      EXPECT_LT(fabs(aedir[kk][ii] - aemd[kk][ii]), EPSILON);
+    }
+    for (int ii = 0; ii < avdir[0].size(); ++ii) {
+      EXPECT_LT(fabs(avdir[kk][ii] - avmd[kk][ii]), EPSILON);
+    }
+  }
+}
+
 TYPED_TEST(TestInferDeepPotModeDevi, cpu_lmp_list) {
   using VALUETYPE = TypeParam;
   std::vector<VALUETYPE>& coord = this->coord;
diff --git a/source/api_cc/include/DeepPot.h b/source/api_cc/include/DeepPot.h
index a8aedde510..fa35d4b13a 100644
--- a/source/api_cc/include/DeepPot.h
+++ b/source/api_cc/include/DeepPot.h
@@ -629,6 +629,71 @@ class DeepPotModelDevi {
             const std::vector<std::string>& file_contents =
                 std::vector<std::string>());
 
+  /**
+   * @brief Evaluate the energy, force and virial by using these DP models.
+   * @param[out] all_ener The system energies of all models.
+   * @param[out] all_force The forces on each atom of all models.
+   * @param[out] all_virial The virials of all models.
+   * @param[in] coord The coordinates of atoms. The array should be of size
+   *nframes x natoms x 3.
+   * @param[in] atype The atom types. The list should contain natoms ints.
+   * @param[in] box The cell of the region. The array should be of size nframes
+   *x 9.
+   * @param[in] fparam The frame parameter. The array can be of size :
+   * nframes x dim_fparam.
+   * dim_fparam. Then all frames are assumed to be provided with the same
+   *fparam.
+   * @param[in] aparam The atomic parameter The array can be of size :
+   * nframes x natoms x dim_aparam.
+   * 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.
+   **/
+  template <typename VALUETYPE>
+  void compute(std::vector<ENERGYTYPE>& all_ener,
+               std::vector<std::vector<VALUETYPE> >& all_force,
+               std::vector<std::vector<VALUETYPE> >& all_virial,
+               const std::vector<VALUETYPE>& coord,
+               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>());
+
+  /**
+   * @brief Evaluate the energy, force, virial, atomic energy, and atomic virial
+   *by using these DP models.
+   * @param[out] all_ener The system energies of all models.
+   * @param[out] all_force The forces on each atom of all models.
+   * @param[out] all_virial The virials of all models.
+   * @param[out] all_atom_energy The atomic energies of all models.
+   * @param[out] all_atom_virial The atomic virials of all models.
+   * @param[in] coord The coordinates of atoms. The array should be of size
+   *nframes x natoms x 3.
+   * @param[in] atype The atom types. The list should contain natoms ints.
+   * @param[in] box The cell of the region. The array should be of size nframes
+   *x 9.
+   * @param[in] fparam The frame parameter. The array can be of size :
+   * nframes x dim_fparam.
+   * dim_fparam. Then all frames are assumed to be provided with the same
+   *fparam.
+   * @param[in] aparam The atomic parameter The array can be of size :
+   * nframes x natoms x dim_aparam.
+   * 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.
+   **/
+  template <typename VALUETYPE>
+  void compute(std::vector<ENERGYTYPE>& all_ener,
+               std::vector<std::vector<VALUETYPE> >& all_force,
+               std::vector<std::vector<VALUETYPE> >& all_virial,
+               std::vector<std::vector<VALUETYPE> >& all_atom_energy,
+               std::vector<std::vector<VALUETYPE> >& all_atom_virial,
+               const std::vector<VALUETYPE>& coord,
+               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>());
+
   /**
    * @brief Evaluate the energy, force and virial by using these DP models.
    * @param[out] all_ener The system energies of all models.
diff --git a/source/api_cc/src/DeepPot.cc b/source/api_cc/src/DeepPot.cc
index 498f35f46b..dd04e8e5e1 100644
--- a/source/api_cc/src/DeepPot.cc
+++ b/source/api_cc/src/DeepPot.cc
@@ -598,6 +598,99 @@ void DeepPotModelDevi::init(const std::vector<std::string>& models,
   inited = true;
 }
 
+template <typename VALUETYPE>
+void DeepPotModelDevi::compute(std::vector<ENERGYTYPE>& all_energy,
+                               std::vector<std::vector<VALUETYPE>>& all_force,
+                               std::vector<std::vector<VALUETYPE>>& all_virial,
+                               const std::vector<VALUETYPE>& dcoord_,
+                               const std::vector<int>& datype_,
+                               const std::vector<VALUETYPE>& dbox,
+                               const std::vector<VALUETYPE>& fparam,
+                               const std::vector<VALUETYPE>& aparam_) {
+  // without nlist
+  if (numb_models == 0) {
+    return;
+  }
+  all_energy.resize(numb_models);
+  all_force.resize(numb_models);
+  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_);
+  }
+}
+
+template void DeepPotModelDevi::compute<double>(
+    std::vector<ENERGYTYPE>& all_energy,
+    std::vector<std::vector<double>>& all_force,
+    std::vector<std::vector<double>>& all_virial,
+    const std::vector<double>& dcoord_,
+    const std::vector<int>& datype_,
+    const std::vector<double>& dbox,
+    const std::vector<double>& fparam,
+    const std::vector<double>& aparam);
+
+template void DeepPotModelDevi::compute<float>(
+    std::vector<ENERGYTYPE>& all_energy,
+    std::vector<std::vector<float>>& all_force,
+    std::vector<std::vector<float>>& all_virial,
+    const std::vector<float>& dcoord_,
+    const std::vector<int>& datype_,
+    const std::vector<float>& dbox,
+    const std::vector<float>& fparam,
+    const std::vector<float>& aparam);
+
+template <typename VALUETYPE>
+void DeepPotModelDevi::compute(
+    std::vector<ENERGYTYPE>& all_energy,
+    std::vector<std::vector<VALUETYPE>>& all_force,
+    std::vector<std::vector<VALUETYPE>>& all_virial,
+    std::vector<std::vector<VALUETYPE>>& all_atom_energy,
+    std::vector<std::vector<VALUETYPE>>& all_atom_virial,
+    const std::vector<VALUETYPE>& dcoord_,
+    const std::vector<int>& datype_,
+    const std::vector<VALUETYPE>& dbox,
+    const std::vector<VALUETYPE>& fparam,
+    const std::vector<VALUETYPE>& aparam_) {
+  if (numb_models == 0) {
+    return;
+  }
+  all_energy.resize(numb_models);
+  all_force.resize(numb_models);
+  all_virial.resize(numb_models);
+  all_atom_energy.resize(numb_models);
+  all_atom_virial.resize(numb_models);
+  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_);
+  }
+}
+
+template void DeepPotModelDevi::compute<double>(
+    std::vector<ENERGYTYPE>& all_energy,
+    std::vector<std::vector<double>>& all_force,
+    std::vector<std::vector<double>>& all_virial,
+    std::vector<std::vector<double>>& all_atom_energy,
+    std::vector<std::vector<double>>& all_atom_virial,
+    const std::vector<double>& dcoord_,
+    const std::vector<int>& datype_,
+    const std::vector<double>& dbox,
+    const std::vector<double>& fparam,
+    const std::vector<double>& aparam);
+
+template void DeepPotModelDevi::compute<float>(
+    std::vector<ENERGYTYPE>& all_energy,
+    std::vector<std::vector<float>>& all_force,
+    std::vector<std::vector<float>>& all_virial,
+    std::vector<std::vector<float>>& all_atom_energy,
+    std::vector<std::vector<float>>& all_atom_virial,
+    const std::vector<float>& dcoord_,
+    const std::vector<int>& datype_,
+    const std::vector<float>& dbox,
+    const std::vector<float>& fparam,
+    const std::vector<float>& aparam);
+
 template <typename VALUETYPE>
 void DeepPotModelDevi::compute(std::vector<ENERGYTYPE>& all_energy,
                                std::vector<std::vector<VALUETYPE>>& all_force,
diff --git a/source/api_cc/tests/test_deeppot_model_devi.cc b/source/api_cc/tests/test_deeppot_model_devi.cc
index 3ceba6a94e..5ce4b11c7f 100644
--- a/source/api_cc/tests/test_deeppot_model_devi.cc
+++ b/source/api_cc/tests/test_deeppot_model_devi.cc
@@ -120,6 +120,90 @@ TYPED_TEST(TestInferDeepPotModeDevi, attrs) {
   EXPECT_EQ(dp1.dim_aparam(), dp_md.dim_aparam());
 }
 
+TYPED_TEST(TestInferDeepPotModeDevi, cpu_build_nlist) {
+  using VALUETYPE = TypeParam;
+  std::vector<VALUETYPE>& coord = this->coord;
+  std::vector<int>& atype = this->atype;
+  std::vector<VALUETYPE>& box = this->box;
+  int& natoms = this->natoms;
+  deepmd::DeepPot& dp0 = this->dp0;
+  deepmd::DeepPot& dp1 = this->dp1;
+  deepmd::DeepPotModelDevi& dp_md = this->dp_md;
+  float rc = dp_md.cutoff();
+  int nloc = coord.size() / 3;
+
+  int nmodel = 2;
+  std::vector<double> edir(nmodel), emd;
+  std::vector<std::vector<VALUETYPE> > fdir(nmodel), vdir(nmodel), fmd(nmodel),
+      vmd;
+  dp0.compute(edir[0], fdir[0], vdir[0], coord, atype, box);
+  dp1.compute(edir[1], fdir[1], vdir[1], coord, atype, box);
+  dp_md.compute(emd, fmd, vmd, coord, atype, box);
+
+  EXPECT_EQ(edir.size(), emd.size());
+  EXPECT_EQ(fdir.size(), fmd.size());
+  EXPECT_EQ(vdir.size(), vmd.size());
+  for (int kk = 0; kk < nmodel; ++kk) {
+    EXPECT_EQ(fdir[kk].size(), fmd[kk].size());
+    EXPECT_EQ(vdir[kk].size(), vmd[kk].size());
+  }
+  for (int kk = 0; kk < nmodel; ++kk) {
+    EXPECT_LT(fabs(edir[kk] - emd[kk]), EPSILON);
+    for (int ii = 0; ii < fdir[0].size(); ++ii) {
+      EXPECT_LT(fabs(fdir[kk][ii] - fmd[kk][ii]), EPSILON);
+    }
+    for (int ii = 0; ii < vdir[0].size(); ++ii) {
+      EXPECT_LT(fabs(vdir[kk][ii] - vmd[kk][ii]), EPSILON);
+    }
+  }
+}
+
+TYPED_TEST(TestInferDeepPotModeDevi, cpu_build_nlist_atomic) {
+  using VALUETYPE = TypeParam;
+  std::vector<VALUETYPE>& coord = this->coord;
+  std::vector<int>& atype = this->atype;
+  std::vector<VALUETYPE>& box = this->box;
+  int& natoms = this->natoms;
+  deepmd::DeepPot& dp0 = this->dp0;
+  deepmd::DeepPot& dp1 = this->dp1;
+  deepmd::DeepPotModelDevi& dp_md = this->dp_md;
+
+  int nmodel = 2;
+  std::vector<double> edir(nmodel), emd;
+  std::vector<std::vector<VALUETYPE> > fdir(nmodel), vdir(nmodel), fmd(nmodel),
+      vmd, aedir(nmodel), aemd, avdir(nmodel), avmd(nmodel);
+  dp0.compute(edir[0], fdir[0], vdir[0], aedir[0], avdir[0], coord, atype, box);
+  dp1.compute(edir[1], fdir[1], vdir[1], aedir[1], avdir[1], coord, atype, box);
+  dp_md.compute(emd, fmd, vmd, aemd, avmd, coord, atype, box);
+
+  EXPECT_EQ(edir.size(), emd.size());
+  EXPECT_EQ(fdir.size(), fmd.size());
+  EXPECT_EQ(vdir.size(), vmd.size());
+  EXPECT_EQ(aedir.size(), aemd.size());
+  EXPECT_EQ(avdir.size(), avmd.size());
+  for (int kk = 0; kk < nmodel; ++kk) {
+    EXPECT_EQ(fdir[kk].size(), fmd[kk].size());
+    EXPECT_EQ(vdir[kk].size(), vmd[kk].size());
+    EXPECT_EQ(aedir[kk].size(), aemd[kk].size());
+    EXPECT_EQ(avdir[kk].size(), avmd[kk].size());
+  }
+  for (int kk = 0; kk < nmodel; ++kk) {
+    EXPECT_LT(fabs(edir[kk] - emd[kk]), EPSILON);
+    for (int ii = 0; ii < fdir[0].size(); ++ii) {
+      EXPECT_LT(fabs(fdir[kk][ii] - fmd[kk][ii]), EPSILON);
+    }
+    for (int ii = 0; ii < vdir[0].size(); ++ii) {
+      EXPECT_LT(fabs(vdir[kk][ii] - vmd[kk][ii]), EPSILON);
+    }
+    for (int ii = 0; ii < aedir[0].size(); ++ii) {
+      EXPECT_LT(fabs(aedir[kk][ii] - aemd[kk][ii]), EPSILON);
+    }
+    for (int ii = 0; ii < avdir[0].size(); ++ii) {
+      EXPECT_LT(fabs(avdir[kk][ii] - avmd[kk][ii]), EPSILON);
+    }
+  }
+}
+
 TYPED_TEST(TestInferDeepPotModeDevi, cpu_lmp_list) {
   using VALUETYPE = TypeParam;
   std::vector<VALUETYPE>& coord = this->coord;
diff --git a/source/api_cc/tests/test_deeppot_model_devi_fparam_aparam.cc b/source/api_cc/tests/test_deeppot_model_devi_fparam_aparam.cc
index 4647b30139..fb1797961d 100644
--- a/source/api_cc/tests/test_deeppot_model_devi_fparam_aparam.cc
+++ b/source/api_cc/tests/test_deeppot_model_devi_fparam_aparam.cc
@@ -19,7 +19,7 @@ class TestInferDeepPotModeDeviFparamAparam : public ::testing::Test {
   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<int> atype = {0, 0, 0, 0, 0, 0};
   std::vector<VALUETYPE> box = {13., 0., 0., 0., 13., 0., 0., 0., 13.};
   std::vector<VALUETYPE> fparam = {0.25852028};
   std::vector<VALUETYPE> aparam = {0.25852028, 0.25852028, 0.25852028,
@@ -76,6 +76,96 @@ TYPED_TEST(TestInferDeepPotModeDeviFparamAparam, attrs) {
   EXPECT_EQ(dp1.dim_aparam(), dp_md.dim_aparam());
 }
 
+TYPED_TEST(TestInferDeepPotModeDeviFparamAparam, cpu_build_nlist) {
+  using VALUETYPE = TypeParam;
+  std::vector<VALUETYPE>& coord = this->coord;
+  std::vector<int>& atype = this->atype;
+  std::vector<VALUETYPE>& box = this->box;
+  std::vector<VALUETYPE>& fparam = this->fparam;
+  std::vector<VALUETYPE>& aparam = this->aparam;
+  int& natoms = this->natoms;
+  deepmd::DeepPot& dp0 = this->dp0;
+  deepmd::DeepPot& dp1 = this->dp1;
+  deepmd::DeepPotModelDevi& dp_md = this->dp_md;
+  float rc = dp_md.cutoff();
+  int nloc = coord.size() / 3;
+
+  int nmodel = 2;
+  std::vector<double> edir(nmodel), emd;
+  std::vector<std::vector<VALUETYPE> > fdir(nmodel), vdir(nmodel), fmd(nmodel),
+      vmd;
+  dp0.compute(edir[0], fdir[0], vdir[0], coord, atype, box, fparam, aparam);
+  dp1.compute(edir[1], fdir[1], vdir[1], coord, atype, box, fparam, aparam);
+  dp_md.compute(emd, fmd, vmd, coord, atype, box, fparam, aparam);
+
+  EXPECT_EQ(edir.size(), emd.size());
+  EXPECT_EQ(fdir.size(), fmd.size());
+  EXPECT_EQ(vdir.size(), vmd.size());
+  for (int kk = 0; kk < nmodel; ++kk) {
+    EXPECT_EQ(fdir[kk].size(), fmd[kk].size());
+    EXPECT_EQ(vdir[kk].size(), vmd[kk].size());
+  }
+  for (int kk = 0; kk < nmodel; ++kk) {
+    EXPECT_LT(fabs(edir[kk] - emd[kk]), EPSILON);
+    for (int ii = 0; ii < fdir[0].size(); ++ii) {
+      EXPECT_LT(fabs(fdir[kk][ii] - fmd[kk][ii]), EPSILON);
+    }
+    for (int ii = 0; ii < vdir[0].size(); ++ii) {
+      EXPECT_LT(fabs(vdir[kk][ii] - vmd[kk][ii]), EPSILON);
+    }
+  }
+}
+
+TYPED_TEST(TestInferDeepPotModeDeviFparamAparam, cpu_build_nlist_atomic) {
+  using VALUETYPE = TypeParam;
+  std::vector<VALUETYPE>& coord = this->coord;
+  std::vector<int>& atype = this->atype;
+  std::vector<VALUETYPE>& box = this->box;
+  std::vector<VALUETYPE>& fparam = this->fparam;
+  std::vector<VALUETYPE>& aparam = this->aparam;
+  int& natoms = this->natoms;
+  deepmd::DeepPot& dp0 = this->dp0;
+  deepmd::DeepPot& dp1 = this->dp1;
+  deepmd::DeepPotModelDevi& dp_md = this->dp_md;
+
+  int nmodel = 2;
+  std::vector<double> edir(nmodel), emd;
+  std::vector<std::vector<VALUETYPE> > fdir(nmodel), vdir(nmodel), fmd(nmodel),
+      vmd, aedir(nmodel), aemd, avdir(nmodel), avmd(nmodel);
+  dp0.compute(edir[0], fdir[0], vdir[0], aedir[0], avdir[0], coord, atype, box,
+              fparam, aparam);
+  dp1.compute(edir[1], fdir[1], vdir[1], aedir[1], avdir[1], coord, atype, box,
+              fparam, aparam);
+  dp_md.compute(emd, fmd, vmd, aemd, avmd, coord, atype, box, fparam, aparam);
+
+  EXPECT_EQ(edir.size(), emd.size());
+  EXPECT_EQ(fdir.size(), fmd.size());
+  EXPECT_EQ(vdir.size(), vmd.size());
+  EXPECT_EQ(aedir.size(), aemd.size());
+  EXPECT_EQ(avdir.size(), avmd.size());
+  for (int kk = 0; kk < nmodel; ++kk) {
+    EXPECT_EQ(fdir[kk].size(), fmd[kk].size());
+    EXPECT_EQ(vdir[kk].size(), vmd[kk].size());
+    EXPECT_EQ(aedir[kk].size(), aemd[kk].size());
+    EXPECT_EQ(avdir[kk].size(), avmd[kk].size());
+  }
+  for (int kk = 0; kk < nmodel; ++kk) {
+    EXPECT_LT(fabs(edir[kk] - emd[kk]), EPSILON);
+    for (int ii = 0; ii < fdir[0].size(); ++ii) {
+      EXPECT_LT(fabs(fdir[kk][ii] - fmd[kk][ii]), EPSILON);
+    }
+    for (int ii = 0; ii < vdir[0].size(); ++ii) {
+      EXPECT_LT(fabs(vdir[kk][ii] - vmd[kk][ii]), EPSILON);
+    }
+    for (int ii = 0; ii < aedir[0].size(); ++ii) {
+      EXPECT_LT(fabs(aedir[kk][ii] - aemd[kk][ii]), EPSILON);
+    }
+    for (int ii = 0; ii < avdir[0].size(); ++ii) {
+      EXPECT_LT(fabs(avdir[kk][ii] - avmd[kk][ii]), EPSILON);
+    }
+  }
+}
+
 TYPED_TEST(TestInferDeepPotModeDeviFparamAparam, cpu_lmp_list) {
   using VALUETYPE = TypeParam;
   std::vector<VALUETYPE>& coord = this->coord;