diff --git a/source/api_cc/include/DeepTensor.h b/source/api_cc/include/DeepTensor.h
index 96cf04fdbf..4e1e47307b 100644
--- a/source/api_cc/include/DeepTensor.h
+++ b/source/api_cc/include/DeepTensor.h
@@ -177,14 +177,6 @@ class DeepTensor
 		  const AtomMap<VALUETYPE> &		atommap, 
 		  const std::vector<int> &		sel_fwd,
 		  const int				nghost = 0);
-  void run_model (std::vector<VALUETYPE> &		dglobal_tensor_,
-		  std::vector<VALUETYPE> &	dforce_,
-		  std::vector<VALUETYPE> &	dvirial_,
-		  tensorflow::Session *			session, 
-		  const std::vector<std::pair<std::string, tensorflow::Tensor>> & input_tensors,
-		  const AtomMap<VALUETYPE> &		atommap, 
-		  const std::vector<int> &		sel_fwd,
-		  const int				nghost = 0);
   void run_model (std::vector<VALUETYPE> &		dglobal_tensor_,
 		  std::vector<VALUETYPE> &	dforce_,
 		  std::vector<VALUETYPE> &	dvirial_,
@@ -205,20 +197,6 @@ class DeepTensor
 		      const std::vector<VALUETYPE> &	box, 
 		      const int				nghost,
 		      const InputNlist&			inlist);
-  void compute_inner (std::vector<VALUETYPE> &		global_tensor,
-		      std::vector<VALUETYPE> &	force,
-		      std::vector<VALUETYPE> &	virial,
-		      const std::vector<VALUETYPE> &	coord,
-		      const std::vector<int> &		atype,
-		      const std::vector<VALUETYPE> &	box);
-  void compute_inner (std::vector<VALUETYPE> &		global_tensor,
-		      std::vector<VALUETYPE> &	force,
-		      std::vector<VALUETYPE> &	virial,
-		      const std::vector<VALUETYPE> &	coord,
-		      const std::vector<int> &		atype,
-		      const std::vector<VALUETYPE> &	box, 
-		      const int				nghost,
-		      const InputNlist&			inlist);
   void compute_inner (std::vector<VALUETYPE> &		global_tensor,
 		      std::vector<VALUETYPE> &	force,
 		      std::vector<VALUETYPE> &	virial,
diff --git a/source/api_cc/include/common.h b/source/api_cc/include/common.h
index 8d9e53c20d..7c769c4b1e 100644
--- a/source/api_cc/include/common.h
+++ b/source/api_cc/include/common.h
@@ -92,6 +92,13 @@ select_map_inv(std::vector<VT> & out,
 	   const std::vector<int > & fwd_map, 
 	   const int & stride);
 
+template<typename VT>
+void 
+select_map_inv(typename std::vector<VT >::iterator out,
+	   const typename std::vector<VT >::const_iterator in, 
+	   const std::vector<int > & fwd_map, 
+	   const int & stride);
+
 /**
 * @brief Get the number of threads from the environment variable.
 * @param[out] num_intra_nthreads The number of intra threads. Read from TF_INTRA_OP_PARALLELISM_THREADS.
diff --git a/source/api_cc/src/DeepTensor.cc b/source/api_cc/src/DeepTensor.cc
index c9d39e42c6..e6ec684f28 100644
--- a/source/api_cc/src/DeepTensor.cc
+++ b/source/api_cc/src/DeepTensor.cc
@@ -130,77 +130,6 @@ run_model (std::vector<VALUETYPE> &	d_tensor_,
   select_map<VALUETYPE>(d_tensor_, d_tensor, sel_srt, odim);
 }
 
-void
-DeepTensor::
-run_model (std::vector<VALUETYPE> &		dglobal_tensor_,
-		  std::vector<VALUETYPE> &	dforce_,
-		  std::vector<VALUETYPE> &	dvirial_,
-		  tensorflow::Session *			session, 
-		  const std::vector<std::pair<std::string, tensorflow::Tensor>> & input_tensors,
-		  const AtomMap<VALUETYPE> &		atommap, 
-		  const std::vector<int> &		sel_fwd,
-		  const int				nghost)
-{
-  unsigned nloc = atommap.get_type().size();
-  unsigned nall = nloc + nghost;
-  if (nloc == 0) {
-    // return empty
-    dglobal_tensor_.clear();
-    dforce_.clear();
-    dvirial_.clear();
-    return;
-  }
-
-  std::vector<Tensor> output_tensors;
-  deepmd::check_status (session->Run(input_tensors, 
-			    {name_prefix(name_scope) + "o_global_" + model_type, 
-			     name_prefix(name_scope) + "o_force", 
-			     name_prefix(name_scope) + "o_virial"},
-			    {}, 
-			    &output_tensors));
-
-  Tensor output_gt = output_tensors[0];
-  Tensor output_f = output_tensors[1];
-  Tensor output_v = output_tensors[2];
-  // this is the new model, output has to be rank 2 tensor
-  assert (output_gt.dims() == 2), "dim of output tensor should be 2";
-  assert (output_f.dims() == 2), "dim of output tensor should be 2";
-  assert (output_v.dims() == 2), "dim of output tensor should be 2";
-  // also check the tensor shapes
-  assert (output_gt.dim_size(0) == 1), "nframes should match";
-  assert (output_gt.dim_size(1) == odim), "dof of global tensor should be odim";  
-  assert (output_f.dim_size(0) == 1), "nframes should match";
-  assert (output_f.dim_size(1) == odim * nall * 3), "dof of force natoms should be odim * nall * 3";
-  assert (output_v.dim_size(0) == 1), "nframes should match";
-  assert (output_v.dim_size(1) == odim * 9), "dof of virial should be odim * 9";
-
-  auto ogt = output_gt.flat <ENERGYTYPE> ();
-  auto of = output_f.flat <VALUETYPE> ();
-  auto ov = output_v.flat <VALUETYPE> ();
-
-  // global tensor
-  dglobal_tensor_.resize(odim);
-  for (unsigned ii = 0; ii < odim; ++ii){
-    dglobal_tensor_[ii] = ogt(ii);
-  }
-
-  // component-wise force
-  std::vector<VALUETYPE> dforce (3 * nall * odim);
-  for (unsigned ii = 0; ii < odim * nall * 3; ++ii){
-    dforce[ii] = of(ii);
-  }
-  dforce_ = dforce;
-  for (unsigned dd = 0; dd < odim; ++dd){
-    atommap.backward (dforce_.begin() + (dd * nall * 3), dforce.begin() + (dd * nall * 3), 3);
-  }
-
-  // component-wise virial
-  dvirial_.resize(odim * 9);
-  for (unsigned ii = 0; ii < odim * 9; ++ii){
-    dvirial_[ii] = ov(ii);
-  }
-}
-
 void
 DeepTensor::
 run_model (std::vector<VALUETYPE> &		dglobal_tensor_,
@@ -367,23 +296,8 @@ compute (std::vector<VALUETYPE> &	dglobal_tensor_,
 	 const std::vector<int> &	datype_,
 	 const std::vector<VALUETYPE> &	dbox)
 {
-  std::vector<VALUETYPE> dcoord, dforce;
-  std::vector<int> datype, fwd_map, bkw_map;
-  int nghost_real;
-  select_real_atoms(fwd_map, bkw_map, nghost_real, dcoord_, datype_, 0, ntypes);
-  assert(nghost_real == 0);
-  // resize to nall_real
-  dcoord.resize(bkw_map.size() * 3);
-  datype.resize(bkw_map.size());
-  // fwd map
-  select_map<VALUETYPE>(dcoord, dcoord_, fwd_map, 3);
-  select_map<int>(datype, datype_, fwd_map, 1);
-  compute_inner(dglobal_tensor_, dforce, dvirial_, dcoord, datype, dbox);
-  // bkw map
-  dforce_.resize(odim * fwd_map.size() * 3);
-  for(int kk = 0; kk < odim; ++kk){
-    select_map<VALUETYPE>(dforce_.begin() + kk * fwd_map.size() * 3, dforce.begin() + kk * bkw_map.size() * 3, bkw_map, 3);
-  }
+  std::vector<VALUETYPE> tmp_at_, tmp_av_;
+  compute(dglobal_tensor_, dforce_, dvirial_, tmp_at_, tmp_av_, dcoord_, datype_, dbox);
 }
 
 void
@@ -397,28 +311,8 @@ compute (std::vector<VALUETYPE> &	dglobal_tensor_,
 	 const int			nghost,
 	 const InputNlist &	lmp_list)
 {
-  std::vector<VALUETYPE> dcoord, dforce;
-  std::vector<int> datype, fwd_map, bkw_map;
-  int nghost_real;
-  select_real_atoms(fwd_map, bkw_map, nghost_real, dcoord_, datype_, nghost, ntypes);
-  // resize to nall_real
-  dcoord.resize(bkw_map.size() * 3);
-  datype.resize(bkw_map.size());
-  // fwd map
-  select_map<VALUETYPE>(dcoord, dcoord_, fwd_map, 3);
-  select_map<int>(datype, datype_, fwd_map, 1);
-  // internal nlist
-  NeighborListData nlist_data;
-  nlist_data.copy_from_nlist(lmp_list);
-  nlist_data.shuffle_exclude_empty(fwd_map);  
-  InputNlist nlist;
-  nlist_data.make_inlist(nlist);
-  compute_inner(dglobal_tensor_, dforce, dvirial_, dcoord, datype, dbox, nghost_real, nlist);
-  // bkw map
-  dforce_.resize(odim * fwd_map.size() * 3);
-  for(int kk = 0; kk < odim; ++kk){
-    select_map<VALUETYPE>(dforce_.begin() + kk * fwd_map.size() * 3, dforce.begin() + kk * bkw_map.size() * 3, bkw_map, 3);
-  }
+  std::vector<VALUETYPE> tmp_at_, tmp_av_;
+  compute(dglobal_tensor_, dforce_, dvirial_, tmp_at_, tmp_av_, dcoord_, datype_, dbox, nghost, lmp_list);
 }
 
 void
@@ -554,67 +448,6 @@ compute_inner (std::vector<VALUETYPE> &		dtensor_,
   run_model (dtensor_, session, input_tensors, atommap, sel_fwd, nghost);
 }
 
-void
-DeepTensor::
-compute_inner (std::vector<VALUETYPE> &		dglobal_tensor_,
-	       std::vector<VALUETYPE> &	dforce_,
-	       std::vector<VALUETYPE> &	dvirial_,
-	       const std::vector<VALUETYPE> &	dcoord_,
-	       const std::vector<int> &		datype_,
-	       const std::vector<VALUETYPE> &	dbox)
-{
-  int nall = dcoord_.size() / 3;
-  int nloc = nall;
-  AtomMap<VALUETYPE> atommap (datype_.begin(), datype_.begin() + nloc);
-  assert (nloc == atommap.get_type().size());
-  
-  std::vector<int> sel_fwd, sel_bkw;
-  int nghost_sel;
-  // this gives the raw selection map, will pass to run model
-  select_by_type(sel_fwd, sel_bkw, nghost_sel, dcoord_, datype_, 0, sel_type);
-
-  std::vector<std::pair<std::string, Tensor>> input_tensors;
-  int ret = session_input_tensors (input_tensors, dcoord_, ntypes, datype_, dbox, cell_size, std::vector<VALUETYPE>(), std::vector<VALUETYPE>(), atommap, name_scope);
-  assert (ret == nloc);
-
-  run_model (dglobal_tensor_, dforce_, dvirial_, session, input_tensors, atommap, sel_fwd);
-}
-
-void
-DeepTensor::
-compute_inner (std::vector<VALUETYPE> &		dglobal_tensor_,
-	       std::vector<VALUETYPE> &	dforce_,
-	       std::vector<VALUETYPE> &	dvirial_,
-	       const std::vector<VALUETYPE> &	dcoord_,
-	       const std::vector<int> &		datype_,
-	       const std::vector<VALUETYPE> &	dbox, 
-	       const int			nghost,
-	       const InputNlist &	nlist_)
-{
-  int nall = dcoord_.size() / 3;
-  int nloc = nall - nghost;
-  AtomMap<VALUETYPE> atommap (datype_.begin(), datype_.begin() + nloc);
-  assert (nloc == atommap.get_type().size());
-
-  std::vector<int> sel_fwd, sel_bkw;
-  int nghost_sel;
-  // this gives the raw selection map, will pass to run model
-  select_by_type(sel_fwd, sel_bkw, nghost_sel, dcoord_, datype_, nghost, sel_type);
-  sel_fwd.resize(nloc);
-
-  NeighborListData nlist_data;
-  nlist_data.copy_from_nlist(nlist_);
-  nlist_data.shuffle(atommap);
-  InputNlist nlist;
-  nlist_data.make_inlist(nlist);
-
-  std::vector<std::pair<std::string, Tensor>> input_tensors;
-  int ret = session_input_tensors (input_tensors, dcoord_, ntypes, datype_, dbox, nlist, std::vector<VALUETYPE>(), std::vector<VALUETYPE>(), atommap, nghost, 0, name_scope);
-  assert (nloc == ret);
-
-  run_model (dglobal_tensor_, dforce_, dvirial_, session, input_tensors, atommap, sel_fwd, nghost);
-}
-
 void
 DeepTensor::
 compute_inner (std::vector<VALUETYPE> &		dglobal_tensor_,
diff --git a/source/api_cc/src/common.cc b/source/api_cc/src/common.cc
index 85eea3683f..92883f04fe 100644
--- a/source/api_cc/src/common.cc
+++ b/source/api_cc/src/common.cc
@@ -618,6 +618,24 @@ select_map_inv(std::vector<VT> & out,
   }
 }
 
+template<typename VT>
+void 
+deepmd::
+select_map_inv(typename std::vector<VT >::iterator out,
+	   const typename std::vector<VT >::const_iterator in, 
+	   const std::vector<int > & idx_map, 
+	   const int & stride)
+{
+  for (int ii = 0; ii < idx_map.size(); ++ii){
+    if (idx_map[ii] >= 0) {
+      int from_ii = idx_map[ii];
+      for (int dd = 0; dd < stride; ++dd){
+	*(out + ii * stride + dd) = *(in + from_ii * stride + dd);
+      }
+    }
+  }
+}
+
 
 template
 int
@@ -641,19 +659,29 @@ select_map<int>(
 template
 void 
 deepmd::
-select_map<int>(typename std::vector<int >::iterator out,
-	   const typename std::vector<int >::const_iterator in, 
-	   const std::vector<int > & idx_map, 
-	   const int & stride);
+select_map<int>(
+    typename std::vector<int >::iterator out,
+    const typename std::vector<int >::const_iterator in, 
+    const std::vector<int > & idx_map, 
+    const int & stride);
 
 template
 void 
 deepmd::
 select_map_inv<int>(
-	   std::vector<int> & out,
-	   const std::vector<int > & in,
-	   const std::vector<int > & idx_map, 
-	   const int & stride);
+    std::vector<int> & out,
+    const std::vector<int > & in,
+    const std::vector<int > & idx_map, 
+    const int & stride);
+
+template
+void 
+deepmd::
+select_map_inv<int>(
+    typename std::vector<int >::iterator out,
+    const typename std::vector<int >::const_iterator in, 
+    const std::vector<int > & idx_map, 
+    const int & stride);
 
 
 template
@@ -678,7 +706,8 @@ select_map<float>(
 template
 void 
 deepmd::
-select_map<float>(typename std::vector<float >::iterator out,
+select_map<float>(
+    typename std::vector<float >::iterator out,
     const typename std::vector<float >::const_iterator in, 
     const std::vector<int > & idx_map, 
     const int & stride);
@@ -692,6 +721,15 @@ select_map_inv<float>(
     const std::vector<int > & idx_map, 
     const int & stride);
 
+template
+void 
+deepmd::
+select_map_inv<float>(
+    typename std::vector<float >::iterator out,
+    const typename std::vector<float >::const_iterator in, 
+    const std::vector<int > & idx_map, 
+    const int & stride);
+
 
 template
 double
@@ -715,7 +753,8 @@ select_map<double>(
 template
 void 
 deepmd::
-select_map<double >(typename std::vector<double >::iterator out,
+select_map<double >(
+    typename std::vector<double >::iterator out,
     const typename std::vector<double >::const_iterator in, 
     const std::vector<int > & idx_map, 
     const int & stride);
@@ -729,6 +768,15 @@ select_map_inv<double>(
     const std::vector<int > & idx_map, 
     const int & stride);
 
+template
+void 
+deepmd::
+select_map_inv<double >(
+    typename std::vector<double >::iterator out,
+    const typename std::vector<double >::const_iterator in, 
+    const std::vector<int > & idx_map, 
+    const int & stride);
+
 
 template
 deepmd::STRINGTYPE
@@ -749,6 +797,15 @@ select_map<deepmd::STRINGTYPE>(
     const std::vector<int > & idx_map, 
     const int & stride);
 
+template
+void 
+deepmd::
+select_map<deepmd::STRINGTYPE >(
+    typename std::vector<deepmd::STRINGTYPE >::iterator out,
+    const typename std::vector<deepmd::STRINGTYPE >::const_iterator in, 
+    const std::vector<int > & idx_map, 
+    const int & stride);
+
 template
 void 
 deepmd::
@@ -756,4 +813,13 @@ select_map_inv<deepmd::STRINGTYPE>(
     std::vector<deepmd::STRINGTYPE> & out,
     const std::vector<deepmd::STRINGTYPE > & in,
     const std::vector<int > & idx_map, 
-    const int & stride);
\ No newline at end of file
+    const int & stride);
+
+template
+void 
+deepmd::
+select_map_inv<deepmd::STRINGTYPE >(
+    typename std::vector<deepmd::STRINGTYPE >::iterator out,
+    const typename std::vector<deepmd::STRINGTYPE >::const_iterator in, 
+    const std::vector<int > & idx_map, 
+    const int & stride);
diff --git a/source/api_cc/tests/test_utils.h b/source/api_cc/tests/test_utils.h
index c4c9338014..d2d9d6c261 100644
--- a/source/api_cc/tests/test_utils.h
+++ b/source/api_cc/tests/test_utils.h
@@ -3,42 +3,35 @@
 
 inline void 
 _fold_back(
-    std::vector<double > &out,
-    const std::vector<double > &in,
+    typename std::vector<double >::iterator out,
+    const typename std::vector<double >::const_iterator in, 
     const std::vector<int> &mapping,
     const int nloc,
     const int nall,
     const int ndim)
 {
-  out.resize(nloc*ndim);
-  std::copy(in.begin(), in.begin() + nloc*ndim, out.begin());
+  // out.resize(nloc*ndim);
+  std::copy(in, in + nloc*ndim, out);
   for(int ii = nloc; ii < nall; ++ii){
     int in_idx = ii;
     int out_idx = mapping[in_idx];
     for(int dd = 0; dd < ndim; ++dd){
-      out[out_idx * ndim + dd] += in[in_idx * ndim + dd];
+      *(out + out_idx * ndim + dd) += *(in + in_idx * ndim + dd);
     }
   }
 }
 
 inline void 
 _fold_back(
-    typename std::vector<double >::iterator out,
-	  const typename std::vector<double >::const_iterator in, 
+    std::vector<double > &out,
+    const std::vector<double > &in,
     const std::vector<int> &mapping,
     const int nloc,
     const int nall,
     const int ndim)
 {
-  // out.resize(nloc*ndim);
-  std::copy(in, in + nloc*ndim, out);
-  for(int ii = nloc; ii < nall; ++ii){
-    int in_idx = ii;
-    int out_idx = mapping[in_idx];
-    for(int dd = 0; dd < ndim; ++dd){
-      *(out + out_idx * ndim + dd) += *(in + in_idx * ndim + dd);
-    }
-  }
+  out.resize(nloc*ndim);
+  _fold_back(out.begin(), in.begin(), mapping, nloc, nall, ndim);
 }
 
 inline void
diff --git a/source/tests/common.py b/source/tests/common.py
index 971d2b0708..5a63666825 100644
--- a/source/tests/common.py
+++ b/source/tests/common.py
@@ -368,8 +368,9 @@ def finite_difference(f, x, delta=1e-6):
     for idx in np.ndindex(*in_shape):
         diff = np.zeros(in_shape)
         diff[idx] += delta
-        y1 = f(x+diff)
-        res[(Ellipsis, *idx)] = (y1-y0) / delta
+        y1p = f(x+diff)
+        y1n = f(x-diff)
+        res[(Ellipsis, *idx)] = (y1p - y1n) / (2 * delta)
     return res