diff --git a/source/api_cc/include/DataModifier.h b/source/api_cc/include/DataModifier.h
index 502d7fcf4e..1e611a3930 100644
--- a/source/api_cc/include/DataModifier.h
+++ b/source/api_cc/include/DataModifier.h
@@ -1,9 +1,92 @@
 // SPDX-License-Identifier: LGPL-3.0-or-later
 #pragma once
 
-#include "DeepPot.h"
+#include <memory>
+
+#include "common.h"
 
 namespace deepmd {
+/**
+ * @brief Dipole charge modifier. (Base class)
+ **/
+class DipoleChargeModifierBase {
+ public:
+  /**
+   * @brief Dipole charge modifier without initialization.
+   **/
+  DipoleChargeModifierBase(){};
+  /**
+   * @brief Dipole charge modifier without initialization.
+   * @param[in] model The name of the frozen model file.
+   * @param[in] gpu_rank The GPU rank. Default is 0.
+   * @param[in] name_scope The name scope.
+   **/
+  DipoleChargeModifierBase(const std::string& model,
+                           const int& gpu_rank = 0,
+                           const std::string& name_scope = "");
+  virtual ~DipoleChargeModifierBase(){};
+  /**
+   * @brief Initialize the dipole charge modifier.
+   * @param[in] model The name of the frozen model file.
+   * @param[in] gpu_rank The GPU rank. Default is 0.
+   * @param[in] name_scope The name scope.
+   **/
+  virtual void init(const std::string& model,
+                    const int& gpu_rank = 0,
+                    const std::string& name_scope = "") = 0;
+  /**
+   * @brief Evaluate the force and virial correction by using this dipole charge
+   *modifier.
+   * @param[out] dfcorr_ The force correction on each atom.
+   * @param[out] dvcorr_ The virial correction.
+   * @param[in] dcoord_ The coordinates of atoms. The array should be of size
+   *natoms x 3.
+   * @param[in] datype_ The atom types. The list should contain natoms ints.
+   * @param[in] dbox The cell of the region. The array should be of size 9.
+   * @param[in] pairs The pairs of atoms. The list should contain npairs pairs
+   *of ints.
+   * @param[in] delef_ The electric field on each atom. The array should be of
+   *size natoms x 3.
+   * @param[in] nghost The number of ghost atoms.
+   * @param[in] lmp_list The neighbor list.
+   @{
+   **/
+  virtual void computew(std::vector<double>& dfcorr_,
+                        std::vector<double>& dvcorr_,
+                        const std::vector<double>& dcoord_,
+                        const std::vector<int>& datype_,
+                        const std::vector<double>& dbox,
+                        const std::vector<std::pair<int, int>>& pairs,
+                        const std::vector<double>& delef_,
+                        const int nghost,
+                        const InputNlist& lmp_list) = 0;
+  virtual void computew(std::vector<float>& dfcorr_,
+                        std::vector<float>& dvcorr_,
+                        const std::vector<float>& dcoord_,
+                        const std::vector<int>& datype_,
+                        const std::vector<float>& dbox,
+                        const std::vector<std::pair<int, int>>& pairs,
+                        const std::vector<float>& delef_,
+                        const int nghost,
+                        const InputNlist& lmp_list) = 0;
+  /** @} */
+  /**
+   * @brief Get cutoff radius.
+   * @return double cutoff radius.
+   */
+  virtual double cutoff() const = 0;
+  /**
+   * @brief Get the number of atom types.
+   * @return int number of atom types.
+   */
+  virtual int numb_types() const = 0;
+  /**
+   * @brief Get the list of sel types.
+   * @return The list of sel types.
+   */
+  virtual std::vector<int> sel_types() const = 0;
+};
+
 /**
  * @brief Dipole charge modifier.
  **/
@@ -38,7 +121,6 @@ class DipoleChargeModifier {
    **/
   void print_summary(const std::string& pre) const;
 
- public:
   /**
    * @brief Evaluate the force and virial correction by using this dipole charge
    *modifier.
@@ -69,50 +151,20 @@ class DipoleChargeModifier {
    * @brief Get cutoff radius.
    * @return double cutoff radius.
    */
-  double cutoff() const {
-    assert(inited);
-    return rcut;
-  };
+  double cutoff() const;
   /**
    * @brief Get the number of atom types.
    * @return int number of atom types.
    */
-  int numb_types() const {
-    assert(inited);
-    return ntypes;
-  };
+  int numb_types() const;
   /**
    * @brief Get the list of sel types.
    * @return The list of sel types.
    */
-  std::vector<int> sel_types() const {
-    assert(inited);
-    return sel_type;
-  };
+  std::vector<int> sel_types() const;
 
  private:
-  tensorflow::Session* session;
-  std::string name_scope, name_prefix;
-  int num_intra_nthreads, num_inter_nthreads;
-  tensorflow::GraphDef* graph_def;
   bool inited;
-  double rcut;
-  int dtype;
-  double cell_size;
-  int ntypes;
-  std::string model_type;
-  std::vector<int> sel_type;
-  template <class VT>
-  VT get_scalar(const std::string& name) const;
-  template <class VT>
-  void get_vector(std::vector<VT>& vec, const std::string& name) const;
-  template <typename MODELTYPE, typename VALUETYPE>
-  void run_model(std::vector<VALUETYPE>& dforce,
-                 std::vector<VALUETYPE>& dvirial,
-                 tensorflow::Session* session,
-                 const std::vector<std::pair<std::string, tensorflow::Tensor>>&
-                     input_tensors,
-                 const AtomMap& atommap,
-                 const int nghost);
+  std::shared_ptr<deepmd::DipoleChargeModifierBase> dcm;
 };
 }  // namespace deepmd
diff --git a/source/api_cc/include/DataModifierTF.h b/source/api_cc/include/DataModifierTF.h
new file mode 100644
index 0000000000..5c44322c0c
--- /dev/null
+++ b/source/api_cc/include/DataModifierTF.h
@@ -0,0 +1,132 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+#pragma once
+
+#include "DataModifier.h"
+#include "common.h"
+
+namespace deepmd {
+/**
+ * @brief Dipole charge modifier.
+ **/
+class DipoleChargeModifierTF : public DipoleChargeModifierBase {
+ public:
+  /**
+   * @brief Dipole charge modifier without initialization.
+   **/
+  DipoleChargeModifierTF();
+  /**
+   * @brief Dipole charge modifier without initialization.
+   * @param[in] model The name of the frozen model file.
+   * @param[in] gpu_rank The GPU rank. Default is 0.
+   * @param[in] name_scope The name scope.
+   **/
+  DipoleChargeModifierTF(const std::string& model,
+                         const int& gpu_rank = 0,
+                         const std::string& name_scope = "");
+  ~DipoleChargeModifierTF();
+  /**
+   * @brief Initialize the dipole charge modifier.
+   * @param[in] model The name of the frozen model file.
+   * @param[in] gpu_rank The GPU rank. Default is 0.
+   * @param[in] name_scope The name scope.
+   **/
+  void init(const std::string& model,
+            const int& gpu_rank = 0,
+            const std::string& name_scope = "");
+
+ public:
+  /**
+   * @brief Evaluate the force and virial correction by using this dipole charge
+   *modifier.
+   * @param[out] dfcorr_ The force correction on each atom.
+   * @param[out] dvcorr_ The virial correction.
+   * @param[in] dcoord_ The coordinates of atoms. The array should be of size
+   *natoms x 3.
+   * @param[in] datype_ The atom types. The list should contain natoms ints.
+   * @param[in] dbox The cell of the region. The array should be of size 9.
+   * @param[in] pairs The pairs of atoms. The list should contain npairs pairs
+   *of ints.
+   * @param[in] delef_ The electric field on each atom. The array should be of
+   *size natoms x 3.
+   * @param[in] nghost The number of ghost atoms.
+   * @param[in] lmp_list The neighbor list.
+   **/
+  template <typename VALUETYPE>
+  void compute(std::vector<VALUETYPE>& dfcorr_,
+               std::vector<VALUETYPE>& dvcorr_,
+               const std::vector<VALUETYPE>& dcoord_,
+               const std::vector<int>& datype_,
+               const std::vector<VALUETYPE>& dbox,
+               const std::vector<std::pair<int, int>>& pairs,
+               const std::vector<VALUETYPE>& delef_,
+               const int nghost,
+               const InputNlist& lmp_list);
+  /**
+   * @brief Get cutoff radius.
+   * @return double cutoff radius.
+   */
+  double cutoff() const {
+    assert(inited);
+    return rcut;
+  };
+  /**
+   * @brief Get the number of atom types.
+   * @return int number of atom types.
+   */
+  int numb_types() const {
+    assert(inited);
+    return ntypes;
+  };
+  /**
+   * @brief Get the list of sel types.
+   * @return The list of sel types.
+   */
+  std::vector<int> sel_types() const {
+    assert(inited);
+    return sel_type;
+  };
+  void computew(std::vector<double>& dfcorr_,
+                std::vector<double>& dvcorr_,
+                const std::vector<double>& dcoord_,
+                const std::vector<int>& datype_,
+                const std::vector<double>& dbox,
+                const std::vector<std::pair<int, int>>& pairs,
+                const std::vector<double>& delef_,
+                const int nghost,
+                const InputNlist& lmp_list);
+  void computew(std::vector<float>& dfcorr_,
+                std::vector<float>& dvcorr_,
+                const std::vector<float>& dcoord_,
+                const std::vector<int>& datype_,
+                const std::vector<float>& dbox,
+                const std::vector<std::pair<int, int>>& pairs,
+                const std::vector<float>& delef_,
+                const int nghost,
+                const InputNlist& lmp_list);
+
+ private:
+  tensorflow::Session* session;
+  std::string name_scope, name_prefix;
+  int num_intra_nthreads, num_inter_nthreads;
+  tensorflow::GraphDef* graph_def;
+  bool inited;
+  double rcut;
+  int dtype;
+  double cell_size;
+  int ntypes;
+  std::string model_type;
+  std::vector<int> sel_type;
+  template <class VT>
+  VT get_scalar(const std::string& name) const;
+  template <class VT>
+  void get_vector(std::vector<VT>& vec, const std::string& name) const;
+  template <typename MODELTYPE, typename VALUETYPE>
+  void run_model(std::vector<VALUETYPE>& dforce,
+                 std::vector<VALUETYPE>& dvirial,
+                 tensorflow::Session* session,
+                 const std::vector<std::pair<std::string, tensorflow::Tensor>>&
+                     input_tensors,
+                 const AtomMap& atommap,
+                 const int nghost);
+};
+}  // namespace deepmd
diff --git a/source/api_cc/src/DataModifier.cc b/source/api_cc/src/DataModifier.cc
index c44bbceaa2..954c969c13 100644
--- a/source/api_cc/src/DataModifier.cc
+++ b/source/api_cc/src/DataModifier.cc
@@ -1,26 +1,21 @@
 // SPDX-License-Identifier: LGPL-3.0-or-later
 #include "DataModifier.h"
 
+#include "DataModifierTF.h"
+#include "common.h"
+
 using namespace deepmd;
-using namespace tensorflow;
 
-DipoleChargeModifier::DipoleChargeModifier()
-    : inited(false), graph_def(new GraphDef()) {}
+DipoleChargeModifier::DipoleChargeModifier() : inited(false) {}
 
 DipoleChargeModifier::DipoleChargeModifier(const std::string& model,
                                            const int& gpu_rank,
                                            const std::string& name_scope_)
-    : inited(false), name_scope(name_scope_), graph_def(new GraphDef()) {
-  try {
-    init(model, gpu_rank, name_scope_);
-  } catch (...) {
-    // Clean up and rethrow, as the destructor will not be called
-    delete graph_def;
-    throw;
-  }
+    : inited(false) {
+  init(model, gpu_rank, name_scope_);
 }
 
-DipoleChargeModifier::~DipoleChargeModifier() { delete graph_def; };
+DipoleChargeModifier::~DipoleChargeModifier(){};
 
 void DipoleChargeModifier::init(const std::string& model,
                                 const int& gpu_rank,
@@ -31,140 +26,22 @@ void DipoleChargeModifier::init(const std::string& model,
               << std::endl;
     return;
   }
-  name_scope = name_scope_;
-  SessionOptions options;
-  get_env_nthreads(num_intra_nthreads, num_inter_nthreads);
-  options.config.set_inter_op_parallelism_threads(num_inter_nthreads);
-  options.config.set_intra_op_parallelism_threads(num_intra_nthreads);
-  deepmd::load_op_library();
-  int gpu_num = -1;
-#if GOOGLE_CUDA || TENSORFLOW_USE_ROCM
-  DPGetDeviceCount(gpu_num);  // check current device environment
-  if (gpu_num > 0) {
-    options.config.set_allow_soft_placement(true);
-    options.config.mutable_gpu_options()->set_per_process_gpu_memory_fraction(
-        0.9);
-    options.config.mutable_gpu_options()->set_allow_growth(true);
-    DPErrcheck(DPSetDevice(gpu_rank % gpu_num));
-    std::string str = "/gpu:";
-    str += std::to_string(gpu_rank % gpu_num);
-    graph::SetDefaultDevice(str, graph_def);
-  }
-#endif  // GOOGLE_CUDA || TENSORFLOW_USE_ROCM
-  deepmd::check_status(NewSession(options, &session));
-  deepmd::check_status(ReadBinaryProto(Env::Default(), model, graph_def));
-  deepmd::check_status(session->Create(*graph_def));
-  // int nnodes = graph_def.node_size();
-  // for (int ii = 0; ii < nnodes; ++ii){
-  //   cout << ii << " \t " << graph_def.node(ii).name() << endl;
-  // }
-  dtype = session_get_dtype(session, "descrpt_attr/rcut");
-  if (dtype == tensorflow::DT_DOUBLE) {
-    rcut = get_scalar<double>("descrpt_attr/rcut");
+  // TODO: To implement detect_backend
+  DPBackend backend = deepmd::DPBackend::TensorFlow;
+  if (deepmd::DPBackend::TensorFlow == backend) {
+    // TODO: throw errors if TF backend is not built, without mentioning TF
+    dcm = std::make_shared<deepmd::DipoleChargeModifierTF>(model, gpu_rank,
+                                                           name_scope_);
+  } else if (deepmd::DPBackend::PyTorch == backend) {
+    throw deepmd::deepmd_exception("PyTorch backend is not supported yet");
+  } else if (deepmd::DPBackend::Paddle == backend) {
+    throw deepmd::deepmd_exception("PaddlePaddle backend is not supported yet");
   } else {
-    rcut = get_scalar<float>("descrpt_attr/rcut");
+    throw deepmd::deepmd_exception("Unknown file type");
   }
-  cell_size = rcut;
-  ntypes = get_scalar<int>("descrpt_attr/ntypes");
-  model_type = get_scalar<STRINGTYPE>("model_attr/model_type");
-  get_vector<int>(sel_type, "model_attr/sel_type");
-  sort(sel_type.begin(), sel_type.end());
   inited = true;
 }
 
-template <class VT>
-VT DipoleChargeModifier::get_scalar(const std::string& name) const {
-  return session_get_scalar<VT>(session, name, name_scope);
-}
-
-template <class VT>
-void DipoleChargeModifier::get_vector(std::vector<VT>& vec,
-                                      const std::string& name) const {
-  session_get_vector<VT>(vec, session, name, name_scope);
-}
-
-template <typename MODELTYPE, typename VALUETYPE>
-void DipoleChargeModifier::run_model(
-    std::vector<VALUETYPE>& dforce,
-    std::vector<VALUETYPE>& dvirial,
-    Session* session,
-    const std::vector<std::pair<std::string, Tensor>>& input_tensors,
-    const AtomMap& atommap,
-    const int nghost) {
-  unsigned nloc = atommap.get_type().size();
-  unsigned nall = nloc + nghost;
-  if (nloc == 0) {
-    dforce.clear();
-    dvirial.clear();
-    return;
-  }
-
-  std::vector<Tensor> output_tensors;
-  deepmd::check_status(session->Run(input_tensors,
-                                    {"o_dm_force", "o_dm_virial", "o_dm_av"},
-                                    {}, &output_tensors));
-  int cc = 0;
-  Tensor output_f = output_tensors[cc++];
-  Tensor output_v = output_tensors[cc++];
-  Tensor output_av = output_tensors[cc++];
-  assert(output_f.dims() == 2 && "dim of output tensor should be 2");
-  assert(output_v.dims() == 2 && "dim of output tensor should be 2");
-  assert(output_av.dims() == 2 && "dim of output tensor should be 2");
-  int nframes = output_f.dim_size(0);
-  int natoms = output_f.dim_size(1) / 3;
-  assert(output_f.dim_size(0) == 1 && "nframes should match");
-  assert(natoms == nall && "natoms should be nall");
-  assert(output_v.dim_size(0) == nframes && "nframes should match");
-  assert(output_v.dim_size(1) == 9 && "dof of virial should be 9");
-  assert(output_av.dim_size(0) == nframes && "nframes should match");
-  assert(output_av.dim_size(1) == natoms * 9 &&
-         "dof of atom virial should be 9 * natoms");
-
-  auto of = output_f.flat<MODELTYPE>();
-  auto ov = output_v.flat<MODELTYPE>();
-
-  dforce.resize(static_cast<size_t>(nall) * 3);
-  dvirial.resize(9);
-  for (int ii = 0; ii < nall * 3; ++ii) {
-    dforce[ii] = of(ii);
-  }
-  for (int ii = 0; ii < 9; ++ii) {
-    dvirial[ii] = ov(ii);
-  }
-}
-
-template void DipoleChargeModifier::run_model<double, double>(
-    std::vector<double>& dforce,
-    std::vector<double>& dvirial,
-    Session* session,
-    const std::vector<std::pair<std::string, Tensor>>& input_tensors,
-    const AtomMap& atommap,
-    const int nghost);
-
-template void DipoleChargeModifier::run_model<float, double>(
-    std::vector<double>& dforce,
-    std::vector<double>& dvirial,
-    Session* session,
-    const std::vector<std::pair<std::string, Tensor>>& input_tensors,
-    const AtomMap& atommap,
-    const int nghost);
-
-template void DipoleChargeModifier::run_model<double, float>(
-    std::vector<float>& dforce,
-    std::vector<float>& dvirial,
-    Session* session,
-    const std::vector<std::pair<std::string, Tensor>>& input_tensors,
-    const AtomMap& atommap,
-    const int nghost);
-
-template void DipoleChargeModifier::run_model<float, float>(
-    std::vector<float>& dforce,
-    std::vector<float>& dvirial,
-    Session* session,
-    const std::vector<std::pair<std::string, Tensor>>& input_tensors,
-    const AtomMap& atommap,
-    const int nghost);
-
 template <typename VALUETYPE>
 void DipoleChargeModifier::compute(
     std::vector<VALUETYPE>& dfcorr_,
@@ -176,148 +53,8 @@ void DipoleChargeModifier::compute(
     const std::vector<VALUETYPE>& delef_,
     const int nghost,
     const InputNlist& lmp_list) {
-  // firstly do selection
-  int nall = datype_.size();
-  int nloc = nall - nghost;
-  int nghost_real;
-  std::vector<int> real_fwd_map, real_bkw_map;
-  select_real_atoms(real_fwd_map, real_bkw_map, nghost_real, dcoord_, datype_,
-                    nghost, ntypes);
-  int nall_real = real_bkw_map.size();
-  int nloc_real = nall_real - nghost_real;
-  if (nloc_real == 0) {
-    dfcorr_.resize(static_cast<size_t>(nall) * 3);
-    dvcorr_.resize(9);
-    fill(dfcorr_.begin(), dfcorr_.end(), (VALUETYPE)0.0);
-    fill(dvcorr_.begin(), dvcorr_.end(), (VALUETYPE)0.0);
-    return;
-  }
-  // resize to nall_real
-  std::vector<VALUETYPE> dcoord_real;
-  std::vector<VALUETYPE> delef_real;
-  std::vector<int> datype_real;
-  dcoord_real.resize(static_cast<size_t>(nall_real) * 3);
-  delef_real.resize(static_cast<size_t>(nall_real) * 3);
-  datype_real.resize(nall_real);
-  // fwd map
-  select_map<VALUETYPE>(dcoord_real, dcoord_, real_fwd_map, 3);
-  select_map<VALUETYPE>(delef_real, delef_, real_fwd_map, 3);
-  select_map<int>(datype_real, datype_, real_fwd_map, 1);
-  // internal nlist
-  NeighborListData nlist_data;
-  nlist_data.copy_from_nlist(lmp_list);
-  nlist_data.shuffle_exclude_empty(real_fwd_map);
-  // sort atoms
-  AtomMap atommap(datype_real.begin(), datype_real.begin() + nloc_real);
-  assert(nloc_real == atommap.get_type().size());
-  const std::vector<int>& sort_fwd_map(atommap.get_fwd_map());
-  const std::vector<int>& sort_bkw_map(atommap.get_bkw_map());
-  // shuffle nlist
-  nlist_data.shuffle(atommap);
-  InputNlist nlist;
-  nlist_data.make_inlist(nlist);
-  // make input tensors
-  std::vector<std::pair<std::string, Tensor>> input_tensors;
-  int ret;
-  if (dtype == tensorflow::DT_DOUBLE) {
-    ret = session_input_tensors<double>(
-        input_tensors, dcoord_real, ntypes, datype_real, dbox, nlist,
-        std::vector<VALUETYPE>(), std::vector<VALUETYPE>(), atommap,
-        nghost_real, 0, name_scope);
-  } else {
-    ret = session_input_tensors<float>(
-        input_tensors, dcoord_real, ntypes, datype_real, dbox, nlist,
-        std::vector<VALUETYPE>(), std::vector<VALUETYPE>(), atommap,
-        nghost_real, 0, name_scope);
-  }
-  assert(nloc_real == ret);
-  // make bond idx map
-  std::vector<int> bd_idx(nall, -1);
-  for (int ii = 0; ii < pairs.size(); ++ii) {
-    bd_idx[pairs[ii].first] = pairs[ii].second;
-  }
-  // make extf by bond idx map
-  std::vector<int> dtype_sort_loc = atommap.get_type();
-  std::vector<VALUETYPE> dextf;
-  for (int ii = 0; ii < dtype_sort_loc.size(); ++ii) {
-    if (binary_search(sel_type.begin(), sel_type.end(), dtype_sort_loc[ii])) {
-      // selected atom
-      int first_idx = real_bkw_map[sort_bkw_map[ii]];
-      int second_idx = bd_idx[first_idx];
-      assert(second_idx >= 0);
-      dextf.push_back(delef_[second_idx * 3 + 0]);
-      dextf.push_back(delef_[second_idx * 3 + 1]);
-      dextf.push_back(delef_[second_idx * 3 + 2]);
-    }
-  }
-  // dextf should be loc and virtual
-  assert(dextf.size() == (nloc - nloc_real) * 3);
-  // make tensor for extf
-  int nframes = 1;
-  TensorShape extf_shape;
-  extf_shape.AddDim(nframes);
-  extf_shape.AddDim(dextf.size());
-  Tensor extf_tensor((tensorflow::DataType)dtype, extf_shape);
-  if (dtype == tensorflow::DT_DOUBLE) {
-    auto extf = extf_tensor.matrix<double>();
-    for (int ii = 0; ii < nframes; ++ii) {
-      for (int jj = 0; jj < extf.size(); ++jj) {
-        extf(ii, jj) = dextf[jj];
-      }
-    }
-  } else {
-    auto extf = extf_tensor.matrix<float>();
-    for (int ii = 0; ii < nframes; ++ii) {
-      for (int jj = 0; jj < extf.size(); ++jj) {
-        extf(ii, jj) = dextf[jj];
-      }
-    }
-  }
-  // append extf to input tensor
-  input_tensors.push_back({"t_ef", extf_tensor});
-  // run model
-  std::vector<VALUETYPE> dfcorr, dvcorr;
-  if (dtype == tensorflow::DT_DOUBLE) {
-    run_model<double>(dfcorr, dvcorr, session, input_tensors, atommap,
-                      nghost_real);
-  } else {
-    run_model<float>(dfcorr, dvcorr, session, input_tensors, atommap,
-                     nghost_real);
-  }
-  assert(dfcorr.size() == nall_real * 3);
-  // back map force
-  std::vector<VALUETYPE> dfcorr_1 = dfcorr;
-  atommap.backward<VALUETYPE>(dfcorr_1.begin(), dfcorr.begin(), 3);
-  assert(dfcorr_1.size() == nall_real * 3);
-  // resize to all and clear
-  std::vector<VALUETYPE> dfcorr_2(nall * 3);
-  fill(dfcorr_2.begin(), dfcorr_2.end(), (VALUETYPE)0.0);
-  // back map to original position
-  for (int ii = 0; ii < nall_real; ++ii) {
-    for (int dd = 0; dd < 3; ++dd) {
-      dfcorr_2[real_bkw_map[ii] * 3 + dd] += dfcorr_1[ii * 3 + dd];
-    }
-  }
-  // self correction of bonded force
-  for (int ii = 0; ii < pairs.size(); ++ii) {
-    for (int dd = 0; dd < 3; ++dd) {
-      dfcorr_2[pairs[ii].first * 3 + dd] += delef_[pairs[ii].second * 3 + dd];
-    }
-  }
-  // add ele contrinution
-  dfcorr_ = dfcorr_2;
-  // for (int ii = 0; ii < nloc; ++ii){
-  //   for (int dd = 0; dd < 3; ++dd){
-  //     dfcorr_[ii*3+dd] += delef_[ii*3+dd];
-  //   }
-  // }
-  for (int ii = 0; ii < nloc_real; ++ii) {
-    int oii = real_bkw_map[ii];
-    for (int dd = 0; dd < 3; ++dd) {
-      dfcorr_[oii * 3 + dd] += delef_[oii * 3 + dd];
-    }
-  }
-  dvcorr_ = dvcorr;
+  dcm->computew(dfcorr_, dvcorr_, dcoord_, datype_, dbox, pairs, delef_, nghost,
+                lmp_list);
 }
 
 template void DipoleChargeModifier::compute<double>(
@@ -345,3 +82,11 @@ template void DipoleChargeModifier::compute<float>(
 void DipoleChargeModifier::print_summary(const std::string& pre) const {
   deepmd::print_summary(pre);
 }
+
+double DipoleChargeModifier::cutoff() const { return dcm->cutoff(); }
+
+int DipoleChargeModifier::numb_types() const { return dcm->numb_types(); }
+
+std::vector<int> DipoleChargeModifier::sel_types() const {
+  return dcm->sel_types();
+}
diff --git a/source/api_cc/src/DataModifierTF.cc b/source/api_cc/src/DataModifierTF.cc
new file mode 100644
index 0000000000..219139cf89
--- /dev/null
+++ b/source/api_cc/src/DataModifierTF.cc
@@ -0,0 +1,363 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+#include "DataModifierTF.h"
+
+#include "common.h"
+
+using namespace deepmd;
+using namespace tensorflow;
+
+DipoleChargeModifierTF::DipoleChargeModifierTF()
+    : inited(false), graph_def(new GraphDef()) {}
+
+DipoleChargeModifierTF::DipoleChargeModifierTF(const std::string& model,
+                                               const int& gpu_rank,
+                                               const std::string& name_scope_)
+    : inited(false), name_scope(name_scope_), graph_def(new GraphDef()) {
+  try {
+    init(model, gpu_rank, name_scope_);
+  } catch (...) {
+    // Clean up and rethrow, as the destructor will not be called
+    delete graph_def;
+    throw;
+  }
+}
+
+DipoleChargeModifierTF::~DipoleChargeModifierTF() { delete graph_def; };
+
+void DipoleChargeModifierTF::init(const std::string& model,
+                                  const int& gpu_rank,
+                                  const std::string& name_scope_) {
+  if (inited) {
+    std::cerr << "WARNING: deepmd-kit should not be initialized twice, do "
+                 "nothing at the second call of initializer"
+              << std::endl;
+    return;
+  }
+  name_scope = name_scope_;
+  SessionOptions options;
+  get_env_nthreads(num_intra_nthreads, num_inter_nthreads);
+  options.config.set_inter_op_parallelism_threads(num_inter_nthreads);
+  options.config.set_intra_op_parallelism_threads(num_intra_nthreads);
+  deepmd::load_op_library();
+  int gpu_num = -1;
+#if GOOGLE_CUDA || TENSORFLOW_USE_ROCM
+  DPGetDeviceCount(gpu_num);  // check current device environment
+  if (gpu_num > 0) {
+    options.config.set_allow_soft_placement(true);
+    options.config.mutable_gpu_options()->set_per_process_gpu_memory_fraction(
+        0.9);
+    options.config.mutable_gpu_options()->set_allow_growth(true);
+    DPErrcheck(DPSetDevice(gpu_rank % gpu_num));
+    std::string str = "/gpu:";
+    str += std::to_string(gpu_rank % gpu_num);
+    graph::SetDefaultDevice(str, graph_def);
+  }
+#endif  // GOOGLE_CUDA || TENSORFLOW_USE_ROCM
+  deepmd::check_status(NewSession(options, &session));
+  deepmd::check_status(ReadBinaryProto(Env::Default(), model, graph_def));
+  deepmd::check_status(session->Create(*graph_def));
+  dtype = session_get_dtype(session, "descrpt_attr/rcut");
+  if (dtype == tensorflow::DT_DOUBLE) {
+    rcut = get_scalar<double>("descrpt_attr/rcut");
+  } else {
+    rcut = get_scalar<float>("descrpt_attr/rcut");
+  }
+  cell_size = rcut;
+  ntypes = get_scalar<int>("descrpt_attr/ntypes");
+  model_type = get_scalar<STRINGTYPE>("model_attr/model_type");
+  get_vector<int>(sel_type, "model_attr/sel_type");
+  sort(sel_type.begin(), sel_type.end());
+  inited = true;
+}
+
+template <class VT>
+VT DipoleChargeModifierTF::get_scalar(const std::string& name) const {
+  return session_get_scalar<VT>(session, name, name_scope);
+}
+
+template <class VT>
+void DipoleChargeModifierTF::get_vector(std::vector<VT>& vec,
+                                        const std::string& name) const {
+  session_get_vector<VT>(vec, session, name, name_scope);
+}
+
+template <typename MODELTYPE, typename VALUETYPE>
+void DipoleChargeModifierTF::run_model(
+    std::vector<VALUETYPE>& dforce,
+    std::vector<VALUETYPE>& dvirial,
+    Session* session,
+    const std::vector<std::pair<std::string, Tensor>>& input_tensors,
+    const AtomMap& atommap,
+    const int nghost) {
+  unsigned nloc = atommap.get_type().size();
+  unsigned nall = nloc + nghost;
+  if (nloc == 0) {
+    dforce.clear();
+    dvirial.clear();
+    return;
+  }
+
+  std::vector<Tensor> output_tensors;
+  deepmd::check_status(session->Run(input_tensors,
+                                    {"o_dm_force", "o_dm_virial", "o_dm_av"},
+                                    {}, &output_tensors));
+  int cc = 0;
+  Tensor output_f = output_tensors[cc++];
+  Tensor output_v = output_tensors[cc++];
+  Tensor output_av = output_tensors[cc++];
+  assert(output_f.dims() == 2 && "dim of output tensor should be 2");
+  assert(output_v.dims() == 2 && "dim of output tensor should be 2");
+  assert(output_av.dims() == 2 && "dim of output tensor should be 2");
+  int nframes = output_f.dim_size(0);
+  int natoms = output_f.dim_size(1) / 3;
+  assert(output_f.dim_size(0) == 1 && "nframes should match");
+  assert(natoms == nall && "natoms should be nall");
+  assert(output_v.dim_size(0) == nframes && "nframes should match");
+  assert(output_v.dim_size(1) == 9 && "dof of virial should be 9");
+  assert(output_av.dim_size(0) == nframes && "nframes should match");
+  assert(output_av.dim_size(1) == natoms * 9 &&
+         "dof of atom virial should be 9 * natoms");
+
+  auto of = output_f.flat<MODELTYPE>();
+  auto ov = output_v.flat<MODELTYPE>();
+
+  dforce.resize(nall * 3);
+  dvirial.resize(9);
+  for (int ii = 0; ii < nall * 3; ++ii) {
+    dforce[ii] = of(ii);
+  }
+  for (int ii = 0; ii < 9; ++ii) {
+    dvirial[ii] = ov(ii);
+  }
+}
+
+template void DipoleChargeModifierTF::run_model<double, double>(
+    std::vector<double>& dforce,
+    std::vector<double>& dvirial,
+    Session* session,
+    const std::vector<std::pair<std::string, Tensor>>& input_tensors,
+    const AtomMap& atommap,
+    const int nghost);
+
+template void DipoleChargeModifierTF::run_model<float, double>(
+    std::vector<double>& dforce,
+    std::vector<double>& dvirial,
+    Session* session,
+    const std::vector<std::pair<std::string, Tensor>>& input_tensors,
+    const AtomMap& atommap,
+    const int nghost);
+
+template void DipoleChargeModifierTF::run_model<double, float>(
+    std::vector<float>& dforce,
+    std::vector<float>& dvirial,
+    Session* session,
+    const std::vector<std::pair<std::string, Tensor>>& input_tensors,
+    const AtomMap& atommap,
+    const int nghost);
+
+template void DipoleChargeModifierTF::run_model<float, float>(
+    std::vector<float>& dforce,
+    std::vector<float>& dvirial,
+    Session* session,
+    const std::vector<std::pair<std::string, Tensor>>& input_tensors,
+    const AtomMap& atommap,
+    const int nghost);
+
+template <typename VALUETYPE>
+void DipoleChargeModifierTF::compute(
+    std::vector<VALUETYPE>& dfcorr_,
+    std::vector<VALUETYPE>& dvcorr_,
+    const std::vector<VALUETYPE>& dcoord_,
+    const std::vector<int>& datype_,
+    const std::vector<VALUETYPE>& dbox,
+    const std::vector<std::pair<int, int>>& pairs,
+    const std::vector<VALUETYPE>& delef_,
+    const int nghost,
+    const InputNlist& lmp_list) {
+  // firstly do selection
+  int nall = datype_.size();
+  int nloc = nall - nghost;
+  int nghost_real;
+  std::vector<int> real_fwd_map, real_bkw_map;
+  select_real_atoms(real_fwd_map, real_bkw_map, nghost_real, dcoord_, datype_,
+                    nghost, ntypes);
+  int nall_real = real_bkw_map.size();
+  int nloc_real = nall_real - nghost_real;
+  if (nloc_real == 0) {
+    dfcorr_.resize(nall * 3);
+    dvcorr_.resize(9);
+    fill(dfcorr_.begin(), dfcorr_.end(), (VALUETYPE)0.0);
+    fill(dvcorr_.begin(), dvcorr_.end(), (VALUETYPE)0.0);
+    return;
+  }
+  // resize to nall_real
+  std::vector<VALUETYPE> dcoord_real;
+  std::vector<VALUETYPE> delef_real;
+  std::vector<int> datype_real;
+  dcoord_real.resize(nall_real * 3);
+  delef_real.resize(nall_real * 3);
+  datype_real.resize(nall_real);
+  // fwd map
+  select_map<VALUETYPE>(dcoord_real, dcoord_, real_fwd_map, 3);
+  select_map<VALUETYPE>(delef_real, delef_, real_fwd_map, 3);
+  select_map<int>(datype_real, datype_, real_fwd_map, 1);
+  // internal nlist
+  NeighborListData nlist_data;
+  nlist_data.copy_from_nlist(lmp_list);
+  nlist_data.shuffle_exclude_empty(real_fwd_map);
+  // sort atoms
+  AtomMap atommap(datype_real.begin(), datype_real.begin() + nloc_real);
+  assert(nloc_real == atommap.get_type().size());
+  const std::vector<int>& sort_fwd_map(atommap.get_fwd_map());
+  const std::vector<int>& sort_bkw_map(atommap.get_bkw_map());
+  // shuffle nlist
+  nlist_data.shuffle(atommap);
+  InputNlist nlist;
+  nlist_data.make_inlist(nlist);
+  // make input tensors
+  std::vector<std::pair<std::string, Tensor>> input_tensors;
+  int ret;
+  if (dtype == tensorflow::DT_DOUBLE) {
+    ret = session_input_tensors<double>(
+        input_tensors, dcoord_real, ntypes, datype_real, dbox, nlist,
+        std::vector<VALUETYPE>(), std::vector<VALUETYPE>(), atommap,
+        nghost_real, 0, name_scope);
+  } else {
+    ret = session_input_tensors<float>(
+        input_tensors, dcoord_real, ntypes, datype_real, dbox, nlist,
+        std::vector<VALUETYPE>(), std::vector<VALUETYPE>(), atommap,
+        nghost_real, 0, name_scope);
+  }
+  assert(nloc_real == ret);
+  // make bond idx map
+  std::vector<int> bd_idx(nall, -1);
+  for (int ii = 0; ii < pairs.size(); ++ii) {
+    bd_idx[pairs[ii].first] = pairs[ii].second;
+  }
+  // make extf by bond idx map
+  std::vector<int> dtype_sort_loc = atommap.get_type();
+  std::vector<VALUETYPE> dextf;
+  for (int ii = 0; ii < dtype_sort_loc.size(); ++ii) {
+    if (binary_search(sel_type.begin(), sel_type.end(), dtype_sort_loc[ii])) {
+      // selected atom
+      int first_idx = real_bkw_map[sort_bkw_map[ii]];
+      int second_idx = bd_idx[first_idx];
+      assert(second_idx >= 0);
+      dextf.push_back(delef_[second_idx * 3 + 0]);
+      dextf.push_back(delef_[second_idx * 3 + 1]);
+      dextf.push_back(delef_[second_idx * 3 + 2]);
+    }
+  }
+  // dextf should be loc and virtual
+  assert(dextf.size() == (nloc - nloc_real) * 3);
+  // make tensor for extf
+  int nframes = 1;
+  TensorShape extf_shape;
+  extf_shape.AddDim(nframes);
+  extf_shape.AddDim(dextf.size());
+  Tensor extf_tensor((tensorflow::DataType)dtype, extf_shape);
+  if (dtype == tensorflow::DT_DOUBLE) {
+    auto extf = extf_tensor.matrix<double>();
+    for (int ii = 0; ii < nframes; ++ii) {
+      for (int jj = 0; jj < extf.size(); ++jj) {
+        extf(ii, jj) = dextf[jj];
+      }
+    }
+  } else {
+    auto extf = extf_tensor.matrix<float>();
+    for (int ii = 0; ii < nframes; ++ii) {
+      for (int jj = 0; jj < extf.size(); ++jj) {
+        extf(ii, jj) = dextf[jj];
+      }
+    }
+  }
+  // append extf to input tensor
+  input_tensors.push_back({"t_ef", extf_tensor});
+  // run model
+  std::vector<VALUETYPE> dfcorr, dvcorr;
+  if (dtype == tensorflow::DT_DOUBLE) {
+    run_model<double>(dfcorr, dvcorr, session, input_tensors, atommap,
+                      nghost_real);
+  } else {
+    run_model<float>(dfcorr, dvcorr, session, input_tensors, atommap,
+                     nghost_real);
+  }
+  assert(dfcorr.size() == nall_real * 3);
+  // back map force
+  std::vector<VALUETYPE> dfcorr_1 = dfcorr;
+  atommap.backward<VALUETYPE>(dfcorr_1.begin(), dfcorr.begin(), 3);
+  assert(dfcorr_1.size() == nall_real * 3);
+  // resize to all and clear
+  std::vector<VALUETYPE> dfcorr_2(nall * 3);
+  fill(dfcorr_2.begin(), dfcorr_2.end(), (VALUETYPE)0.0);
+  // back map to original position
+  for (int ii = 0; ii < nall_real; ++ii) {
+    for (int dd = 0; dd < 3; ++dd) {
+      dfcorr_2[real_bkw_map[ii] * 3 + dd] += dfcorr_1[ii * 3 + dd];
+    }
+  }
+  // self correction of bonded force
+  for (int ii = 0; ii < pairs.size(); ++ii) {
+    for (int dd = 0; dd < 3; ++dd) {
+      dfcorr_2[pairs[ii].first * 3 + dd] += delef_[pairs[ii].second * 3 + dd];
+    }
+  }
+  // add ele contrinution
+  dfcorr_ = dfcorr_2;
+  for (int ii = 0; ii < nloc_real; ++ii) {
+    int oii = real_bkw_map[ii];
+    for (int dd = 0; dd < 3; ++dd) {
+      dfcorr_[oii * 3 + dd] += delef_[oii * 3 + dd];
+    }
+  }
+  dvcorr_ = dvcorr;
+}
+
+template void DipoleChargeModifierTF::compute<double>(
+    std::vector<double>& dfcorr_,
+    std::vector<double>& dvcorr_,
+    const std::vector<double>& dcoord_,
+    const std::vector<int>& datype_,
+    const std::vector<double>& dbox,
+    const std::vector<std::pair<int, int>>& pairs,
+    const std::vector<double>& delef_,
+    const int nghost,
+    const InputNlist& lmp_list);
+
+template void DipoleChargeModifierTF::compute<float>(
+    std::vector<float>& dfcorr_,
+    std::vector<float>& dvcorr_,
+    const std::vector<float>& dcoord_,
+    const std::vector<int>& datype_,
+    const std::vector<float>& dbox,
+    const std::vector<std::pair<int, int>>& pairs,
+    const std::vector<float>& delef_,
+    const int nghost,
+    const InputNlist& lmp_list);
+
+void DipoleChargeModifierTF::computew(
+    std::vector<double>& dfcorr_,
+    std::vector<double>& dvcorr_,
+    const std::vector<double>& dcoord_,
+    const std::vector<int>& datype_,
+    const std::vector<double>& dbox,
+    const std::vector<std::pair<int, int>>& pairs,
+    const std::vector<double>& delef_,
+    const int nghost,
+    const InputNlist& lmp_list) {
+  compute(dfcorr_, dvcorr_, dcoord_, datype_, dbox, pairs, delef_, nghost,
+          lmp_list);
+}
+void DipoleChargeModifierTF::computew(
+    std::vector<float>& dfcorr_,
+    std::vector<float>& dvcorr_,
+    const std::vector<float>& dcoord_,
+    const std::vector<int>& datype_,
+    const std::vector<float>& dbox,
+    const std::vector<std::pair<int, int>>& pairs,
+    const std::vector<float>& delef_,
+    const int nghost,
+    const InputNlist& lmp_list) {
+  compute(dfcorr_, dvcorr_, dcoord_, datype_, dbox, pairs, delef_, nghost,
+          lmp_list);
+}