diff --git a/source/api_cc/tests/.gitignore b/source/api_cc/tests/.gitignore
new file mode 100644
index 0000000000..2c32bfb8cf
--- /dev/null
+++ b/source/api_cc/tests/.gitignore
@@ -0,0 +1,5 @@
+CTestTestfile.cmake
+Makefile
+cmake_install.cmake
+runUnitTests
+version.h
diff --git a/source/api_cc/tests/CMakeLists.txt b/source/api_cc/tests/CMakeLists.txt
new file mode 100644
index 0000000000..3cdddd0d56
--- /dev/null
+++ b/source/api_cc/tests/CMakeLists.txt
@@ -0,0 +1,108 @@
+cmake_minimum_required(VERSION 3.9)
+project(deepmd_api_test)
+set(CMAKE_LINK_WHAT_YOU_USE TRUE)
+
+if (NOT DEFINED BUILD_CPP_IF) 
+  set(BUILD_CPP_IF TRUE)
+endif (NOT DEFINED BUILD_CPP_IF)
+add_definitions ("-DHIGH_PREC")
+
+enable_testing()
+
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11")
+
+# model version
+file(READ ${PROJECT_SOURCE_DIR}/../../config/MODEL_VER MODEL_VERSION)
+string(REPLACE "\n" " " MODEL_VERSION ${MODEL_VERSION})
+message(STATUS "Supported model version: ${MODEL_VERSION}")
+
+set(libname "deepmd")
+set(LIB_BASE_DIR ${CMAKE_SOURCE_DIR}/../../lib)
+include_directories(${LIB_BASE_DIR}/include)
+file(GLOB LIB_SRC ${LIB_BASE_DIR}/src/*.cc ${LIB_BASE_DIR}/src/*.cpp)
+add_library(${libname} SHARED ${LIB_SRC})
+
+set(apiname "deepmd_api")
+set(API_BASE_DIR ${CMAKE_SOURCE_DIR}/../)
+include_directories(${API_BASE_DIR}/include)
+include_directories(${CMAKE_SOURCE_DIR})
+file(GLOB API_SRC ${API_BASE_DIR}/src/*.cc ${API_BASE_DIR}/src/*.cpp)
+add_library(${apiname} SHARED ${API_SRC})
+configure_file(
+  ${API_BASE_DIR}/include/version.h.in
+  ${CMAKE_SOURCE_DIR}/version.h
+  @ONLY
+)
+
+set(opname "deepmd_op")
+set(OP_BASE_DIR ${CMAKE_SOURCE_DIR}/../../op)
+# file(GLOB OP_SRC ${OP_BASE_DIR}/*.cc)
+file(GLOB OP_SRC ${OP_BASE_DIR}/prod_force.cc ${OP_BASE_DIR}/prod_virial.cc ${OP_BASE_DIR}/descrpt.cc ${OP_BASE_DIR}/descrpt_se_a_ef.cc ${OP_BASE_DIR}/descrpt_se_a_ef.cc ${OP_BASE_DIR}/descrpt_se_a_ef_para.cc ${OP_BASE_DIR}/descrpt_se_a_ef_vert.cc ${OP_BASE_DIR}/pair_tab.cc ${OP_BASE_DIR}/prod_force_multi_device.cc ${OP_BASE_DIR}/prod_virial_multi_device.cc ${OP_BASE_DIR}/soft_min.cc ${OP_BASE_DIR}/soft_min_force.cc ${OP_BASE_DIR}/soft_min_virial.cc ${OP_BASE_DIR}/ewald_recp.cc ${OP_BASE_DIR}/gelu_multi_device.cc ${OP_BASE_DIR}/map_aparam.cc ${OP_BASE_DIR}/neighbor_stat.cc ${OP_BASE_DIR}/unaggregated_grad.cc ${OP_BASE_DIR}/tabulate_multi_device.cc ${OP_BASE_DIR}/prod_env_mat_multi_device.cc)
+add_library(${opname} SHARED ${OP_SRC})
+
+list (APPEND CMAKE_MODULE_PATH ${CMAKE_SOURCE_DIR}/../../cmake/)
+find_package(tensorflow REQUIRED)
+include_directories(${TensorFlow_INCLUDE_DIRS})
+
+find_package(Threads)
+# find openmp
+find_package(OpenMP)
+if (OPENMP_FOUND)
+    set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}")
+    set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
+endif()
+
+# define USE_CUDA_TOOLKIT
+if (DEFINED USE_CUDA_TOOLKIT)
+  if (USE_CUDA_TOOLKIT)
+    find_package(CUDA REQUIRED)
+  else()
+    message(STATUS "Will not build nv GPU support")
+  endif()
+else()
+  find_package(CUDA QUIET)  
+  if (CUDA_FOUND)
+    set(USE_CUDA_TOOLKIT TRUE)
+    message(STATUS "Found CUDA in ${CUDA_TOOLKIT_ROOT_DIR}, build nv GPU support")
+  else()
+    set(USE_CUDA_TOOLKIT FALSE)
+    message(STATUS "No cuda support found, will not build nv GPU support")
+  endif()
+endif()
+
+if (USE_CUDA_TOOLKIT)
+  add_definitions("-D GOOGLE_CUDA")
+  include_directories(${CUDA_INCLUDE_DIRS})
+  add_subdirectory(${LIB_BASE_DIR}/src/cuda cuda_binary_dir)
+endif()
+
+file(GLOB TEST_SRC test_*.cc)
+add_executable( runUnitTests ${TEST_SRC} )
+
+if (USE_CUDA_TOOLKIT)
+  target_link_libraries(runUnitTests gtest gtest_main ${libname} pthread deepmd_op_cuda)
+else()
+  target_link_libraries(runUnitTests gtest gtest_main ${libname} ${apiname} ${opname} pthread ${TensorFlow_LIBRARY})
+endif()
+add_test( runUnitTests runUnitTests )
+
+find_package(GTest)
+if(NOT GTEST_LIBRARY)
+  configure_file(../../cmake/googletest.cmake.in googletest-download/CMakeLists.txt)
+  execute_process(COMMAND ${CMAKE_COMMAND} -G "${CMAKE_GENERATOR}" .
+    RESULT_VARIABLE result
+    WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/googletest-download )
+  if(result)
+    message(FATAL_ERROR "CMake step for googletest failed: ${result}")
+  endif()
+  execute_process(COMMAND ${CMAKE_COMMAND} --build .
+    RESULT_VARIABLE result
+    WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/googletest-download )
+  if(result)
+    message(FATAL_ERROR "Build step for googletest failed: ${result}")
+  endif()
+  set(gtest_force_shared_crt ON CACHE BOOL "" FORCE)
+  add_subdirectory(${CMAKE_CURRENT_BINARY_DIR}/googletest-src ${CMAKE_CURRENT_BINARY_DIR}/googletest-build EXCLUDE_FROM_ALL)
+else ()
+  include_directories(${GTEST_INCLUDE_DIRS})
+endif ()
diff --git a/source/api_cc/tests/test_deepdipole.cc b/source/api_cc/tests/test_deepdipole.cc
new file mode 100644
index 0000000000..fbc7388f7a
--- /dev/null
+++ b/source/api_cc/tests/test_deepdipole.cc
@@ -0,0 +1,105 @@
+#include <gtest/gtest.h>
+#include <cmath>
+#include <algorithm>
+#include <fstream>
+#include <vector>
+#include "DeepTensor.h"
+#include "neighbor_list.h"
+#include "test_utils.h"
+
+#include "google/protobuf/text_format.h"
+#include "google/protobuf/io/zero_copy_stream_impl.h"
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>  
+
+class TestInferDeepDipole : public ::testing::Test
+{  
+protected:  
+  std::vector<double> 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<double> box = {
+    13., 0., 0., 0., 13., 0., 0., 0., 13.
+  };
+  std::vector<double> expected_d = {
+    -9.274180565967479195e-01,2.698028341272042496e+00,2.521268387140979117e-01,2.927260638453461628e+00,-8.571926301526779923e-01,1.667785136187720063e+00
+  };
+  int natoms;
+
+  deepmd::DeepTensor dp;
+
+  void SetUp() override {
+    std::string file_name = "../../tests/infer/deepdipole.pbtxt";
+    int fd = open(file_name.c_str(), O_RDONLY);
+    tensorflow::protobuf::io::ZeroCopyInputStream* input = new tensorflow::protobuf::io::FileInputStream(fd);
+    tensorflow::GraphDef graph_def;
+    tensorflow::protobuf::TextFormat::Parse(input, &graph_def);
+    delete input;
+    std::fstream output("deepdipole.pb", std::ios::out | std::ios::trunc | std::ios::binary);
+    graph_def.SerializeToOstream(&output);
+    // check the string by the following commands
+    // string txt;
+    // tensorflow::protobuf::TextFormat::PrintToString(graph_def, &txt);
+
+    dp.init("deepdipole.pb");
+
+    natoms = expected_d.size();
+  };
+
+  void TearDown() override {
+    remove( "deepdipole.pb" ) ;
+  };
+};
+
+
+TEST_F(TestInferDeepDipole, cpu_build_nlist)
+{
+  EXPECT_EQ(dp.cutoff(), 4.);
+  EXPECT_EQ(dp.numb_types(), 2);
+  EXPECT_EQ(dp.output_dim(), 3);
+  std::vector<int> sel_types = dp.sel_types();
+  EXPECT_EQ(sel_types.size(), 1);
+  EXPECT_EQ(sel_types[0], 0);
+
+  std::vector<double> value;
+  dp.compute(value, coord, atype, box);
+
+  EXPECT_EQ(value.size(), expected_d.size());
+  for(int ii = 0; ii < expected_d.size(); ++ii){
+    EXPECT_LT(fabs(value[ii] - expected_d[ii]), 1e-10);
+  }
+}
+
+TEST_F(TestInferDeepDipole, cpu_lmp_nlist)
+{
+  float rc = dp.cutoff();
+  int nloc = coord.size() / 3;  
+  std::vector<double> coord_cpy;
+  std::vector<int> atype_cpy, mapping;  
+  std::vector<int> ilist(nloc), numneigh(nloc);
+  std::vector<int*> firstneigh(nloc);
+  std::vector<std::vector<int > > nlist_data;
+  deepmd::InputNlist inlist(nloc, &ilist[0], &numneigh[0], &firstneigh[0]);
+  _build_nlist(nlist_data, coord_cpy, atype_cpy, mapping,
+	       coord, atype, box, rc);
+  int nall = coord_cpy.size() / 3;
+  convert_nlist(inlist, nlist_data);  
+
+  std::vector<double> value;
+  dp.compute(value, coord_cpy, atype_cpy, box, nall-nloc, inlist);
+
+  EXPECT_EQ(value.size(), expected_d.size());
+  for(int ii = 0; ii < expected_d.size(); ++ii){
+    EXPECT_LT(fabs(value[ii] - expected_d[ii]), 1e-10);
+  }
+}
+
diff --git a/source/api_cc/tests/test_deeppolar.cc b/source/api_cc/tests/test_deeppolar.cc
new file mode 100644
index 0000000000..a76e84179a
--- /dev/null
+++ b/source/api_cc/tests/test_deeppolar.cc
@@ -0,0 +1,105 @@
+#include <gtest/gtest.h>
+#include <cmath>
+#include <algorithm>
+#include <fstream>
+#include <vector>
+#include "DeepTensor.h"
+#include "neighbor_list.h"
+#include "test_utils.h"
+
+#include "google/protobuf/text_format.h"
+#include "google/protobuf/io/zero_copy_stream_impl.h"
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>  
+
+class TestInferDeepPolar : public ::testing::Test
+{  
+protected:  
+  std::vector<double> 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<double> box = {
+    13., 0., 0., 0., 13., 0., 0., 0., 13.
+  };
+  std::vector<double> expected_d = {
+    1.061407927405987051e-01,-3.569013342133873778e-01,-2.862108976089940138e-02,-3.569013342133875444e-01,1.304367268874677244e+00,1.037647501453442256e-01,-2.862108976089940138e-02,1.037647501453441284e-01,8.100521520762453409e-03,1.236797829492216616e+00,-3.717307430531632262e-01,7.371515676976750919e-01,-3.717307430531630041e-01,1.127222682121889058e-01,-2.239181552775717510e-01,7.371515676976746478e-01,-2.239181552775717787e-01,4.448255365635306879e-01
+  };
+  int natoms;
+
+  deepmd::DeepTensor dp;
+
+  void SetUp() override {
+    std::string file_name = "../../tests/infer/deeppolar.pbtxt";
+    int fd = open(file_name.c_str(), O_RDONLY);
+    tensorflow::protobuf::io::ZeroCopyInputStream* input = new tensorflow::protobuf::io::FileInputStream(fd);
+    tensorflow::GraphDef graph_def;
+    tensorflow::protobuf::TextFormat::Parse(input, &graph_def);
+    delete input;
+    std::fstream output("deeppolar.pb", std::ios::out | std::ios::trunc | std::ios::binary);
+    graph_def.SerializeToOstream(&output);
+    // check the string by the following commands
+    // string txt;
+    // tensorflow::protobuf::TextFormat::PrintToString(graph_def, &txt);
+
+    dp.init("deeppolar.pb");
+
+    natoms = expected_d.size();
+  };
+
+  void TearDown() override {
+    remove( "deeppolar.pb" ) ;
+  };
+};
+
+
+TEST_F(TestInferDeepPolar, cpu_build_nlist)
+{
+  EXPECT_EQ(dp.cutoff(), 6.);
+  EXPECT_EQ(dp.numb_types(), 2);
+  EXPECT_EQ(dp.output_dim(), 9);
+  std::vector<int> sel_types = dp.sel_types();
+  EXPECT_EQ(sel_types.size(), 1);
+  EXPECT_EQ(sel_types[0], 0);
+
+  std::vector<double> value;
+  dp.compute(value, coord, atype, box);
+
+  EXPECT_EQ(value.size(), expected_d.size());
+  for(int ii = 0; ii < expected_d.size(); ++ii){
+    EXPECT_LT(fabs(value[ii] - expected_d[ii]), 1e-10);
+  }
+}
+
+TEST_F(TestInferDeepPolar, cpu_lmp_nlist)
+{
+  float rc = dp.cutoff();
+  int nloc = coord.size() / 3;  
+  std::vector<double> coord_cpy;
+  std::vector<int> atype_cpy, mapping;  
+  std::vector<int> ilist(nloc), numneigh(nloc);
+  std::vector<int*> firstneigh(nloc);
+  std::vector<std::vector<int > > nlist_data;
+  deepmd::InputNlist inlist(nloc, &ilist[0], &numneigh[0], &firstneigh[0]);
+  _build_nlist(nlist_data, coord_cpy, atype_cpy, mapping,
+	       coord, atype, box, rc);
+  int nall = coord_cpy.size() / 3;
+  convert_nlist(inlist, nlist_data);  
+
+  std::vector<double> value;
+  dp.compute(value, coord_cpy, atype_cpy, box, nall-nloc, inlist);
+
+  EXPECT_EQ(value.size(), expected_d.size());
+  for(int ii = 0; ii < expected_d.size(); ++ii){
+    EXPECT_LT(fabs(value[ii] - expected_d[ii]), 1e-10);
+  }
+}
+
diff --git a/source/api_cc/tests/test_deeppot_a.cc b/source/api_cc/tests/test_deeppot_a.cc
new file mode 100644
index 0000000000..61336c4e33
--- /dev/null
+++ b/source/api_cc/tests/test_deeppot_a.cc
@@ -0,0 +1,480 @@
+#include <gtest/gtest.h>
+#include <cmath>
+#include <algorithm>
+#include <fstream>
+#include <vector>
+#include "DeepPot.h"
+#include "neighbor_list.h"
+#include "test_utils.h"
+
+#include "google/protobuf/text_format.h"
+#include "google/protobuf/io/zero_copy_stream_impl.h"
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>  
+
+class TestInferDeepPotA : public ::testing::Test
+{  
+protected:  
+  std::vector<double> 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<double> box = {
+    13., 0., 0., 0., 13., 0., 0., 0., 13.
+  };
+  std::vector<double> expected_e = {
+    -9.275780747115504710e+01,-1.863501786584258468e+02,-1.863392472863538103e+02,-9.279281325486221021e+01,-1.863671545232153903e+02,-1.863619822847602165e+02
+  };
+  std::vector<double> expected_f = {
+    -3.034045420701179663e-01,8.405844663871177014e-01,7.696947487118485642e-02,7.662001266663505117e-01,-1.880601391333554251e-01,-6.183333871091722944e-01,-5.036172391059643427e-01,-6.529525836149027151e-01,5.432962643022043459e-01,6.382357912332115024e-01,-1.748518296794561167e-01,3.457363524891907125e-01,1.286482986991941552e-03,3.757251165286925043e-01,-5.972588700887541124e-01,-5.987006197104716154e-01,-2.004450304880958100e-01,2.495901655353461868e-01
+  };
+  std::vector<double> expected_v = {
+    -2.912234126853306959e-01,-3.800610846612756388e-02,2.776624987489437202e-01,-5.053761003913598976e-02,-3.152373041953385746e-01,1.060894290092162379e-01,2.826389131596073745e-01,1.039129970665329250e-01,-2.584378792325942586e-01,-3.121722367954994914e-01,8.483275876786681990e-02,2.524662342344257682e-01,4.142176771106586414e-02,-3.820285230785245428e-02,-2.727311173065460545e-02,2.668859789777112135e-01,-6.448243569420382404e-02,-2.121731470426218846e-01,-8.624335220278558922e-02,-1.809695356746038597e-01,1.529875294531883312e-01,-1.283658185172031341e-01,-1.992682279795223999e-01,1.409924999632362341e-01,1.398322735274434292e-01,1.804318474574856390e-01,-1.470309318999652726e-01,-2.593983661598450730e-01,-4.236536279233147489e-02,3.386387920184946720e-02,-4.174017537818433543e-02,-1.003500282164128260e-01,1.525690815194478966e-01,3.398976109910181037e-02,1.522253908435125536e-01,-2.349125581341701963e-01,9.515545977581392825e-04,-1.643218849228543846e-02,1.993234765412972564e-02,6.027265332209678569e-04,-9.563256398907417355e-02,1.510815124001868293e-01,-7.738094816888557714e-03,1.502832772532304295e-01,-2.380965783745832010e-01,-2.309456719810296654e-01,-6.666961081213038098e-02,7.955566551234216632e-02,-8.099093777937517447e-02,-3.386641099800401927e-02,4.447884755740908608e-02,1.008593228579038742e-01,4.556718179228393811e-02,-6.078081273849572641e-02
+  };
+  int natoms;
+  double expected_tot_e;
+  std::vector<double>expected_tot_v;
+
+  deepmd::DeepPot dp;
+
+  void SetUp() override {
+    std::string file_name = "../../tests/infer/deeppot.pbtxt";
+    int fd = open(file_name.c_str(), O_RDONLY);
+    tensorflow::protobuf::io::ZeroCopyInputStream* input = new tensorflow::protobuf::io::FileInputStream(fd);
+    tensorflow::GraphDef graph_def;
+    tensorflow::protobuf::TextFormat::Parse(input, &graph_def);
+    delete input;
+    std::fstream output("deeppot.pb", std::ios::out | std::ios::trunc | std::ios::binary);
+    graph_def.SerializeToOstream(&output);
+    // check the string by the following commands
+    // string txt;
+    // tensorflow::protobuf::TextFormat::PrintToString(graph_def, &txt);
+
+    dp.init("deeppot.pb");
+
+    natoms = expected_e.size();
+    EXPECT_EQ(natoms * 3, expected_f.size());
+    EXPECT_EQ(natoms * 9, expected_v.size());
+    expected_tot_e = 0.;
+    expected_tot_v.resize(9);
+    std::fill(expected_tot_v.begin(), expected_tot_v.end(), 0.);
+    for(int ii = 0; ii < natoms; ++ii){
+      expected_tot_e += expected_e[ii];
+    }
+    for(int ii = 0; ii < natoms; ++ii){
+      for(int dd = 0; dd < 9; ++dd){
+	expected_tot_v[dd] += expected_v[ii*9+dd];
+      }
+    }
+  };
+
+  void TearDown() override {
+    remove( "deeppot.pb" ) ;
+  };
+};
+
+TEST_F(TestInferDeepPotA, cpu_build_nlist)
+{
+  double ener;
+  std::vector<double> force, virial;
+  dp.compute(ener, force, virial, coord, atype, box);
+
+  EXPECT_EQ(force.size(), natoms*3);
+  EXPECT_EQ(virial.size(), 9);
+
+  EXPECT_LT(fabs(ener - expected_tot_e), 1e-10);
+  for(int ii = 0; ii < natoms*3; ++ii){
+    EXPECT_LT(fabs(force[ii] - expected_f[ii]), 1e-10);    
+  }
+  for(int ii = 0; ii < 3*3; ++ii){
+    EXPECT_LT(fabs(virial[ii] - expected_tot_v[ii]), 1e-10);
+  }
+}
+
+TEST_F(TestInferDeepPotA, cpu_build_nlist_numfv)
+{
+  class MyModel : public EnergyModelTest<double>
+  {
+    deepmd::DeepPot & mydp;
+    const std::vector<int > atype;
+public:
+    MyModel(
+	deepmd::DeepPot & dp_,
+	const std::vector<int> & atype_
+	) : mydp(dp_), atype(atype_) {};
+    virtual void compute (
+	double & ener,
+	std::vector<double> &	force,
+	std::vector<double> &	virial,
+	const std::vector<double> & coord,
+	const std::vector<double> & box) {
+      mydp.compute(ener, force, virial, coord, atype, box);
+    }
+  };
+  MyModel model(dp, atype);
+  model.test_f(coord, box);
+  model.test_v(coord, box);
+  std::vector<double> box_(box);
+  box_[1] -= 0.4;
+  model.test_f(coord, box_);
+  model.test_v(coord, box_);
+  box_[2] += 0.5;
+  model.test_f(coord, box_);
+  model.test_v(coord, box_);
+  box_[4] += 0.2;
+  model.test_f(coord, box_);
+  model.test_v(coord, box_);
+  box_[3] -= 0.3;
+  model.test_f(coord, box_);
+  model.test_v(coord, box_);
+  box_[6] -= 0.7;
+  model.test_f(coord, box_);
+  model.test_v(coord, box_);
+  box_[7] += 0.6;
+  model.test_f(coord, box_);
+  model.test_v(coord, box_);
+}
+
+
+TEST_F(TestInferDeepPotA, cpu_build_nlist_atomic)
+{
+  double ener;
+  std::vector<double> force, virial, atom_ener, atom_vir;
+  dp.compute(ener, force, virial, atom_ener, atom_vir, coord, atype, box);
+  
+  EXPECT_EQ(force.size(), natoms*3);
+  EXPECT_EQ(virial.size(), 9);
+  EXPECT_EQ(atom_ener.size(), natoms);
+  EXPECT_EQ(atom_vir.size(), natoms*9);
+
+  EXPECT_LT(fabs(ener - expected_tot_e), 1e-10);
+  for(int ii = 0; ii < natoms*3; ++ii){
+    EXPECT_LT(fabs(force[ii] - expected_f[ii]), 1e-10);    
+  }
+  for(int ii = 0; ii < 3*3; ++ii){
+    EXPECT_LT(fabs(virial[ii] - expected_tot_v[ii]), 1e-10);
+  }
+  for(int ii = 0; ii < natoms; ++ii){
+    EXPECT_LT(fabs(atom_ener[ii] - expected_e[ii]), 1e-10);
+  }
+  for(int ii = 0; ii < natoms*9; ++ii){
+    EXPECT_LT(fabs(atom_vir[ii] - expected_v[ii]), 1e-10);
+  }
+}
+
+
+TEST_F(TestInferDeepPotA, cpu_lmp_nlist)
+{
+  float rc = dp.cutoff();
+  int nloc = coord.size() / 3;  
+  std::vector<double> coord_cpy;
+  std::vector<int> atype_cpy, mapping;  
+  std::vector<std::vector<int > > nlist_data;
+  _build_nlist(nlist_data, coord_cpy, atype_cpy, mapping,
+	       coord, atype, box, rc);
+  int nall = coord_cpy.size() / 3;
+  std::vector<int> ilist(nloc), numneigh(nloc);
+  std::vector<int*> firstneigh(nloc);
+  deepmd::InputNlist inlist(nloc, &ilist[0], &numneigh[0], &firstneigh[0]);
+  convert_nlist(inlist, nlist_data);  
+  
+  double ener;
+  std::vector<double> force_, virial;
+  dp.compute(ener, force_, virial, coord_cpy, atype_cpy, box, nall-nloc, inlist, 0);
+  std::vector<double> force;
+  _fold_back(force, force_, mapping, nloc, nall, 3);
+
+  EXPECT_EQ(force.size(), natoms*3);
+  EXPECT_EQ(virial.size(), 9);
+
+  EXPECT_LT(fabs(ener - expected_tot_e), 1e-10);
+  for(int ii = 0; ii < natoms*3; ++ii){
+    EXPECT_LT(fabs(force[ii] - expected_f[ii]), 1e-10);    
+  }
+  for(int ii = 0; ii < 3*3; ++ii){
+    EXPECT_LT(fabs(virial[ii] - expected_tot_v[ii]), 1e-10);
+  }
+
+  ener = 0.;
+  std::fill(force_.begin(), force_.end(), 0.0);
+  std::fill(virial.begin(), virial.end(), 0.0);
+  dp.compute(ener, force_, virial, coord_cpy, atype_cpy, box, nall-nloc, inlist, 1);
+  _fold_back(force, force_, mapping, nloc, nall, 3);
+
+  EXPECT_EQ(force.size(), natoms*3);
+  EXPECT_EQ(virial.size(), 9);
+
+  EXPECT_LT(fabs(ener - expected_tot_e), 1e-10);
+  for(int ii = 0; ii < natoms*3; ++ii){
+    EXPECT_LT(fabs(force[ii] - expected_f[ii]), 1e-10);    
+  }
+  for(int ii = 0; ii < 3*3; ++ii){
+    EXPECT_LT(fabs(virial[ii] - expected_tot_v[ii]), 1e-10);
+  }
+}
+
+
+TEST_F(TestInferDeepPotA, cpu_lmp_nlist_atomic)
+{
+  float rc = dp.cutoff();
+  int nloc = coord.size() / 3;  
+  std::vector<double> coord_cpy;
+  std::vector<int> atype_cpy, mapping;  
+  std::vector<std::vector<int > > nlist_data;
+  _build_nlist(nlist_data, coord_cpy, atype_cpy, mapping,
+	       coord, atype, box, rc);
+  int nall = coord_cpy.size() / 3;
+  std::vector<int> ilist(nloc), numneigh(nloc);
+  std::vector<int*> firstneigh(nloc);
+  deepmd::InputNlist inlist(nloc, &ilist[0], &numneigh[0], &firstneigh[0]);
+  convert_nlist(inlist, nlist_data);  
+  
+  double ener;
+  std::vector<double> force_, atom_ener_, atom_vir_, virial;
+  std::vector<double> force, atom_ener, atom_vir;
+  dp.compute(ener, force_, virial, atom_ener_, atom_vir_, coord_cpy, atype_cpy, box, nall-nloc, inlist, 0);
+  _fold_back(force, force_, mapping, nloc, nall, 3);
+  _fold_back(atom_ener, atom_ener_, mapping, nloc, nall, 1);
+  _fold_back(atom_vir, atom_vir_, mapping, nloc, nall, 9);
+
+  EXPECT_EQ(force.size(), natoms*3);
+  EXPECT_EQ(virial.size(), 9);
+  EXPECT_EQ(atom_ener.size(), natoms);
+  EXPECT_EQ(atom_vir.size(), natoms*9);
+
+  EXPECT_LT(fabs(ener - expected_tot_e), 1e-10);
+  for(int ii = 0; ii < natoms*3; ++ii){
+    EXPECT_LT(fabs(force[ii] - expected_f[ii]), 1e-10);    
+  }
+  for(int ii = 0; ii < 3*3; ++ii){
+    EXPECT_LT(fabs(virial[ii] - expected_tot_v[ii]), 1e-10);
+  }
+  for(int ii = 0; ii < natoms; ++ii){
+    EXPECT_LT(fabs(atom_ener[ii] - expected_e[ii]), 1e-10);
+  }
+  for(int ii = 0; ii < natoms*9; ++ii){
+    EXPECT_LT(fabs(atom_vir[ii] - expected_v[ii]), 1e-10);
+  }
+
+  ener = 0.;
+  std::fill(force_.begin(), force_.end(), 0.0);
+  std::fill(virial.begin(), virial.end(), 0.0);
+  std::fill(atom_ener_.begin(), atom_ener_.end(), 0.0);
+  std::fill(atom_vir_.begin(), atom_vir_.end(), 0.0);  
+  dp.compute(ener, force_, virial, atom_ener_, atom_vir_, coord_cpy, atype_cpy, box, nall-nloc, inlist, 1);
+  _fold_back(force, force_, mapping, nloc, nall, 3);
+  _fold_back(atom_ener, atom_ener_, mapping, nloc, nall, 1);
+  _fold_back(atom_vir, atom_vir_, mapping, nloc, nall, 9);
+
+  EXPECT_EQ(force.size(), natoms*3);
+  EXPECT_EQ(virial.size(), 9);
+  EXPECT_EQ(atom_ener.size(), natoms);
+  EXPECT_EQ(atom_vir.size(), natoms*9);
+
+  EXPECT_LT(fabs(ener - expected_tot_e), 1e-10);
+  for(int ii = 0; ii < natoms*3; ++ii){
+    EXPECT_LT(fabs(force[ii] - expected_f[ii]), 1e-10);    
+  }
+  for(int ii = 0; ii < 3*3; ++ii){
+    EXPECT_LT(fabs(virial[ii] - expected_tot_v[ii]), 1e-10);
+  }
+  for(int ii = 0; ii < natoms; ++ii){
+    EXPECT_LT(fabs(atom_ener[ii] - expected_e[ii]), 1e-10);
+  }
+  for(int ii = 0; ii < natoms*9; ++ii){
+    EXPECT_LT(fabs(atom_vir[ii] - expected_v[ii]), 1e-10);
+  }
+}
+
+
+TEST_F(TestInferDeepPotA, cpu_lmp_nlist_2rc)
+{
+  float rc = dp.cutoff();
+  int nloc = coord.size() / 3;  
+  std::vector<double> coord_cpy;
+  std::vector<int> atype_cpy, mapping;  
+  std::vector<std::vector<int > > nlist_data;
+  _build_nlist(nlist_data, coord_cpy, atype_cpy, mapping,
+	       coord, atype, box, rc*2);
+  int nall = coord_cpy.size() / 3;
+  std::vector<int> ilist(nloc), numneigh(nloc);
+  std::vector<int*> firstneigh(nloc);
+  deepmd::InputNlist inlist(nloc, &ilist[0], &numneigh[0], &firstneigh[0]);
+  convert_nlist(inlist, nlist_data);  
+  
+  double ener;
+  std::vector<double> force_(nall*3, 0.0), virial(9, 0.0);
+  dp.compute(ener, force_, virial, coord_cpy, atype_cpy, box, nall-nloc, inlist, 0);
+  std::vector<double> force;
+  _fold_back(force, force_, mapping, nloc, nall, 3);
+
+  EXPECT_EQ(force.size(), natoms*3);
+  EXPECT_EQ(virial.size(), 9);
+
+  EXPECT_LT(fabs(ener - expected_tot_e), 1e-10);
+  for(int ii = 0; ii < natoms*3; ++ii){
+    EXPECT_LT(fabs(force[ii] - expected_f[ii]), 1e-10);    
+  }
+  for(int ii = 0; ii < 3*3; ++ii){
+    EXPECT_LT(fabs(virial[ii] - expected_tot_v[ii]), 1e-10);
+  }
+
+  ener = 0.;
+  std::fill(force_.begin(), force_.end(), 0.0);
+  std::fill(virial.begin(), virial.end(), 0.0);
+  dp.compute(ener, force_, virial, coord_cpy, atype_cpy, box, nall-nloc, inlist, 1);
+  _fold_back(force, force_, mapping, nloc, nall, 3);
+
+  EXPECT_EQ(force.size(), natoms*3);
+  EXPECT_EQ(virial.size(), 9);
+
+  EXPECT_LT(fabs(ener - expected_tot_e), 1e-10);
+  for(int ii = 0; ii < natoms*3; ++ii){
+    EXPECT_LT(fabs(force[ii] - expected_f[ii]), 1e-10);    
+  }
+  for(int ii = 0; ii < 3*3; ++ii){
+    EXPECT_LT(fabs(virial[ii] - expected_tot_v[ii]), 1e-10);
+  }
+}
+
+
+TEST_F(TestInferDeepPotA, cpu_lmp_nlist_type_sel)
+{
+  float rc = dp.cutoff();
+
+  // add vir atoms
+  int nvir = 2;
+  std::vector<double> coord_vir(nvir*3);
+  std::vector<int> atype_vir(nvir, 2);
+  for(int ii = 0; ii < nvir; ++ii){
+    coord_vir[ii] = coord[ii];
+  }  
+  coord.insert(coord.begin(), coord_vir.begin(), coord_vir.end());
+  atype.insert(atype.begin(), atype_vir.begin(), atype_vir.end());
+  natoms += nvir;
+  std::vector<double> expected_f_vir(nvir*3, 0.0);
+  expected_f.insert(expected_f.begin(), expected_f_vir.begin(), expected_f_vir.end());
+
+  // build nlist
+  int nloc = coord.size() / 3;  
+  std::vector<double> coord_cpy;
+  std::vector<int> atype_cpy, mapping;  
+  std::vector<std::vector<int > > nlist_data;
+  _build_nlist(nlist_data, coord_cpy, atype_cpy, mapping,
+	       coord, atype, box, rc);
+  int nall = coord_cpy.size() / 3;
+  std::vector<int> ilist(nloc), numneigh(nloc);
+  std::vector<int*> firstneigh(nloc);
+  deepmd::InputNlist inlist(nloc, &ilist[0], &numneigh[0], &firstneigh[0]);
+  convert_nlist(inlist, nlist_data);  
+
+  // dp compute
+  double ener;
+  std::vector<double> force_(nall*3, 0.0), virial(9, 0.0);
+  dp.compute(ener, force_, virial, coord_cpy, atype_cpy, box, nall-nloc, inlist, 0);
+  // fold back
+  std::vector<double> force;
+  _fold_back(force, force_, mapping, nloc, nall, 3);
+
+  EXPECT_EQ(force.size(), natoms*3);
+  EXPECT_EQ(virial.size(), 9);
+
+  EXPECT_LT(fabs(ener - expected_tot_e), 1e-10);
+  for(int ii = 0; ii < natoms*3; ++ii){
+    EXPECT_LT(fabs(force[ii] - expected_f[ii]), 1e-10);    
+  }
+  for(int ii = 0; ii < 3*3; ++ii){
+    EXPECT_LT(fabs(virial[ii] - expected_tot_v[ii]), 1e-10);
+  }
+}
+
+
+
+class TestInferDeepPotANoPbc : public ::testing::Test
+{  
+protected:  
+  std::vector<double> 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<double> box = {};
+  std::vector<double> expected_e = {
+    -9.255934839310273787e+01,-1.863253376736990106e+02,-1.857237299341402945e+02,-9.279308539717486326e+01,-1.863708105823244239e+02,-1.863635196514972563e+02
+  };
+  std::vector<double> expected_f = {
+    -2.161037360255332107e+00,9.052994347015581589e-01,1.635379623977007979e+00,2.161037360255332107e+00,-9.052994347015581589e-01,-1.635379623977007979e+00,-1.167128117249453811e-02,1.371975700096064992e-03,-1.575265180249604477e-03,6.226508593971802341e-01,-1.816734122009256991e-01,3.561766019664774907e-01,-1.406075393906316626e-02,3.789140061530929526e-01,-6.018777878642909140e-01,-5.969188242856223736e-01,-1.986125696522633155e-01,2.472764510780630642e-01    
+  };
+  std::vector<double> expected_v = {
+    -7.042445481792056761e-01,2.950213647777754078e-01,5.329418202437231633e-01,2.950213647777752968e-01,-1.235900311906896754e-01,-2.232594111831812944e-01,5.329418202437232743e-01,-2.232594111831813499e-01,-4.033073234276823849e-01,-8.949230984097404917e-01,3.749002169013777030e-01,6.772391014992630298e-01,3.749002169013777586e-01,-1.570527935667933583e-01,-2.837082722496912512e-01,6.772391014992631408e-01,-2.837082722496912512e-01,-5.125052659994422388e-01,4.858210330291591605e-02,-6.902596153269104431e-03,6.682612642430500391e-03,-5.612247004554610057e-03,9.767795567660207592e-04,-9.773758942738038254e-04,5.638322117219018645e-03,-9.483806049779926932e-04,8.493873281881353637e-04,-2.941738570564985666e-01,-4.482529909499673171e-02,4.091569840186781021e-02,-4.509020615859140463e-02,-1.013919988807244071e-01,1.551440772665269030e-01,4.181857726606644232e-02,1.547200233064863484e-01,-2.398213304685777592e-01,-3.218625798524068354e-02,-1.012438450438508421e-02,1.271639330380921855e-02,3.072814938490859779e-03,-9.556241797915024372e-02,1.512251983492413077e-01,-8.277872384009607454e-03,1.505412040827929787e-01,-2.386150620881526407e-01,-2.312295470054945568e-01,-6.631490213524345034e-02,7.932427266386249398e-02,-8.053754366323923053e-02,-3.294595881137418747e-02,4.342495071150231922e-02,1.004599500126941436e-01,4.450400364869536163e-02,-5.951077548033092968e-02
+  };
+  int natoms;
+  double expected_tot_e;
+  std::vector<double>expected_tot_v;
+
+  deepmd::DeepPot dp;
+
+  void SetUp() override {
+    std::string file_name = "../../tests/infer/deeppot.pbtxt";
+    int fd = open(file_name.c_str(), O_RDONLY);
+    tensorflow::protobuf::io::ZeroCopyInputStream* input = new tensorflow::protobuf::io::FileInputStream(fd);
+    tensorflow::GraphDef graph_def;
+    tensorflow::protobuf::TextFormat::Parse(input, &graph_def);
+    delete input;
+    std::fstream output("deeppot.pb", std::ios::out | std::ios::trunc | std::ios::binary);
+    graph_def.SerializeToOstream(&output);
+
+    dp.init("deeppot.pb");
+
+    natoms = expected_e.size();
+    EXPECT_EQ(natoms * 3, expected_f.size());
+    EXPECT_EQ(natoms * 9, expected_v.size());
+    expected_tot_e = 0.;
+    expected_tot_v.resize(9);
+    std::fill(expected_tot_v.begin(), expected_tot_v.end(), 0.);
+    for(int ii = 0; ii < natoms; ++ii){
+      expected_tot_e += expected_e[ii];
+    }
+    for(int ii = 0; ii < natoms; ++ii){
+      for(int dd = 0; dd < 9; ++dd){
+	expected_tot_v[dd] += expected_v[ii*9+dd];
+      }
+    }
+  };
+
+  void TearDown() override {
+    remove( "deeppot.pb" ) ;
+  };
+};
+
+TEST_F(TestInferDeepPotANoPbc, cpu_build_nlist)
+{
+  double ener;
+  std::vector<double> force, virial;
+  dp.compute(ener, force, virial, coord, atype, box);
+
+  EXPECT_EQ(force.size(), natoms*3);
+  EXPECT_EQ(virial.size(), 9);
+
+  EXPECT_LT(fabs(ener - expected_tot_e), 1e-10);
+  for(int ii = 0; ii < natoms*3; ++ii){
+    EXPECT_LT(fabs(force[ii] - expected_f[ii]), 1e-10);    
+  }
+  for(int ii = 0; ii < 3*3; ++ii){
+    EXPECT_LT(fabs(virial[ii] - expected_tot_v[ii]), 1e-10);
+  }
+}
diff --git a/source/api_cc/tests/test_deeppot_model_devi.cc b/source/api_cc/tests/test_deeppot_model_devi.cc
new file mode 100644
index 0000000000..678b30a2ea
--- /dev/null
+++ b/source/api_cc/tests/test_deeppot_model_devi.cc
@@ -0,0 +1,216 @@
+#include <gtest/gtest.h>
+#include <cmath>
+#include <algorithm>
+#include <fstream>
+#include <vector>
+#include "DeepPot.h"
+#include "neighbor_list.h"
+#include "test_utils.h"
+
+#include "google/protobuf/text_format.h"
+#include "google/protobuf/io/zero_copy_stream_impl.h"
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>  
+
+class TestInferDeepPotModeDevi : public ::testing::Test
+{  
+protected:  
+  std::vector<double> 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<double> box = {
+    13., 0., 0., 0., 13., 0., 0., 0., 13.
+  };
+  int natoms;
+
+  deepmd::DeepPot dp0;
+  deepmd::DeepPot dp1;
+  deepmd::DeepPotModelDevi dp_md;
+
+  void SetUp() override {
+    {
+      std::string file_name = "../../tests/infer/deeppot.pbtxt";
+      int fd = open(file_name.c_str(), O_RDONLY);
+      tensorflow::protobuf::io::ZeroCopyInputStream* input = new tensorflow::protobuf::io::FileInputStream(fd);
+      tensorflow::GraphDef graph_def;
+      tensorflow::protobuf::TextFormat::Parse(input, &graph_def);
+      delete input;
+      std::fstream output("deeppot.pb", std::ios::out | std::ios::trunc | std::ios::binary);
+      graph_def.SerializeToOstream(&output);
+      dp0.init("deeppot.pb");
+    }
+    {
+      std::string file_name = "../../tests/infer/deeppot-1.pbtxt";
+      int fd = open(file_name.c_str(), O_RDONLY);
+      tensorflow::protobuf::io::ZeroCopyInputStream* input = new tensorflow::protobuf::io::FileInputStream(fd);
+      tensorflow::GraphDef graph_def;
+      tensorflow::protobuf::TextFormat::Parse(input, &graph_def);
+      delete input;
+      std::fstream output("deeppot-1.pb", std::ios::out | std::ios::trunc | std::ios::binary);
+      graph_def.SerializeToOstream(&output);
+      dp1.init("deeppot-1.pb");
+    }
+    dp_md.init(std::vector<std::string>({"deeppot.pb", "deeppot-1.pb"}));
+  };
+
+  void TearDown() override {
+    remove( "deeppot.pb" ) ;
+    remove( "deeppot-1.pb" ) ;
+  };
+};
+
+TEST_F(TestInferDeepPotModeDevi, attrs)
+{
+  EXPECT_EQ(dp0.cutoff(), dp_md.cutoff());
+  EXPECT_EQ(dp0.numb_types(), dp_md.numb_types());
+  EXPECT_EQ(dp0.dim_fparam(), dp_md.dim_fparam());
+  EXPECT_EQ(dp0.dim_aparam(), dp_md.dim_aparam());
+  EXPECT_EQ(dp1.cutoff(), dp_md.cutoff());
+  EXPECT_EQ(dp1.numb_types(), dp_md.numb_types());
+  EXPECT_EQ(dp1.dim_fparam(), dp_md.dim_fparam());
+  EXPECT_EQ(dp1.dim_aparam(), dp_md.dim_aparam());
+}
+
+
+TEST_F(TestInferDeepPotModeDevi, cpu_lmp_list)
+{
+  float rc = dp_md.cutoff();
+  int nloc = coord.size() / 3;  
+  std::vector<double> coord_cpy;
+  std::vector<int> atype_cpy, mapping;  
+  std::vector<std::vector<int > > nlist_data;
+  _build_nlist(nlist_data, coord_cpy, atype_cpy, mapping,
+	       coord, atype, box, rc);
+  int nall = coord_cpy.size() / 3;
+  std::vector<int> ilist(nloc), numneigh(nloc);
+  std::vector<int*> firstneigh(nloc);
+  deepmd::InputNlist inlist(nloc, &ilist[0], &numneigh[0], &firstneigh[0]);
+  convert_nlist(inlist, nlist_data);  
+
+  int nmodel = 2;
+  std::vector<double > edir(nmodel), emd;
+  std::vector<std::vector<double> > fdir_(nmodel), fdir(nmodel), vdir(nmodel), fmd_, fmd(nmodel), vmd;
+  dp0.compute(edir[0], fdir_[0], vdir[0], coord_cpy, atype_cpy, box, nall-nloc, inlist, 0);
+  dp1.compute(edir[1], fdir_[1], vdir[1], coord_cpy, atype_cpy, box, nall-nloc, inlist, 0);
+  dp_md.compute(emd, fmd_, vmd, coord_cpy, atype_cpy, box, nall-nloc, inlist, 0);
+  for(int kk = 0; kk < nmodel; ++kk){
+    _fold_back(fdir[kk], fdir_[kk], mapping, nloc, nall, 3);
+    _fold_back(fmd[kk], fmd_[kk], mapping, nloc, nall, 3);
+  }
+  
+
+  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]), 1e-10);
+    for(int ii = 0; ii < fdir[0].size(); ++ii){
+      EXPECT_LT(fabs(fdir[kk][ii] - fmd[kk][ii]), 1e-10);
+    }
+    for(int ii = 0; ii < vdir[0].size(); ++ii){
+      EXPECT_LT(fabs(vdir[kk][ii] - vmd[kk][ii]), 1e-10);
+    }
+  }
+
+  std::vector<double > avg_f, std_f;
+  dp_md.compute_avg(avg_f, fmd);
+  dp_md.compute_std_f(std_f, avg_f, fmd);
+  
+  // manual compute std f
+  std::vector<double > manual_std_f(nloc);
+  EXPECT_EQ(fmd[0].size(), nloc * 3);
+  for(int ii = 0; ii < nloc; ++ii){
+    std::vector<double > avg_f(3, 0.0);
+    for(int dd = 0; dd < 3; ++dd){
+      for(int kk = 0; kk < nmodel; ++kk){
+	avg_f[dd] += fmd[kk][ii*3+dd];
+      }
+      avg_f[dd] /= (nmodel) * 1.0;
+    }
+    double std = 0.;
+    for(int kk = 0; kk < nmodel; ++kk){
+      for(int dd = 0; dd < 3; ++dd){
+	double tmp = fmd[kk][ii*3+dd] - avg_f[dd];
+	std += tmp * tmp;
+      }
+    }
+    std /= nmodel * 1.0;
+    manual_std_f[ii] = sqrt(std);
+  }
+
+  EXPECT_EQ(manual_std_f.size(), std_f.size());
+  for(int ii = 0; ii < std_f.size(); ++ii){
+    EXPECT_LT(fabs(manual_std_f[ii] - std_f[ii]), 1e-10);
+  }
+}
+
+
+TEST_F(TestInferDeepPotModeDevi, cpu_lmp_list_atomic)
+{
+  float rc = dp_md.cutoff();
+  int nloc = coord.size() / 3;  
+  std::vector<double> coord_cpy;
+  std::vector<int> atype_cpy, mapping;  
+  std::vector<std::vector<int > > nlist_data;
+  _build_nlist(nlist_data, coord_cpy, atype_cpy, mapping,
+	       coord, atype, box, rc);
+  int nall = coord_cpy.size() / 3;
+  std::vector<int> ilist(nloc), numneigh(nloc);
+  std::vector<int*> firstneigh(nloc);
+  deepmd::InputNlist inlist(nloc, &ilist[0], &numneigh[0], &firstneigh[0]);
+  convert_nlist(inlist, nlist_data);  
+
+  int nmodel = 2;
+  std::vector<double > edir(nmodel), emd;
+  std::vector<std::vector<double> > fdir_(nmodel), fdir(nmodel), vdir(nmodel), fmd_, fmd(nmodel), vmd, aedir(nmodel), aemd, avdir(nmodel), avdir_(nmodel), avmd(nmodel), avmd_;
+  dp0.compute(edir[0], fdir_[0], vdir[0], aedir[0], avdir_[0], coord_cpy, atype_cpy, box, nall-nloc, inlist, 0);
+  dp1.compute(edir[1], fdir_[1], vdir[1], aedir[1], avdir_[1], coord_cpy, atype_cpy, box, nall-nloc, inlist, 0);
+  dp_md.compute(emd, fmd_, vmd, aemd, avmd_, coord_cpy, atype_cpy, box, nall-nloc, inlist, 0);
+  for(int kk = 0; kk < nmodel; ++kk){
+    _fold_back(fdir[kk], fdir_[kk], mapping, nloc, nall, 3);
+    _fold_back(fmd[kk], fmd_[kk], mapping, nloc, nall, 3);
+    _fold_back(avdir[kk], avdir_[kk], mapping, nloc, nall, 9);
+    _fold_back(avmd[kk], avmd_[kk], mapping, nloc, nall, 9);
+  }  
+
+  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]), 1e-10);
+    for(int ii = 0; ii < fdir[0].size(); ++ii){
+      EXPECT_LT(fabs(fdir[kk][ii] - fmd[kk][ii]), 1e-10);
+    }
+    for(int ii = 0; ii < vdir[0].size(); ++ii){
+      EXPECT_LT(fabs(vdir[kk][ii] - vmd[kk][ii]), 1e-10);
+    }
+    for(int ii = 0; ii < aedir[0].size(); ++ii){
+      EXPECT_LT(fabs(aedir[kk][ii] - aemd[kk][ii]), 1e-10);
+    }
+    for(int ii = 0; ii < avdir[0].size(); ++ii){
+      EXPECT_LT(fabs(avdir[kk][ii] - avmd[kk][ii]), 1e-10);
+    }
+  }
+}
+
diff --git a/source/api_cc/tests/test_deeppot_r.cc b/source/api_cc/tests/test_deeppot_r.cc
new file mode 100644
index 0000000000..2d6af9f6ae
--- /dev/null
+++ b/source/api_cc/tests/test_deeppot_r.cc
@@ -0,0 +1,480 @@
+#include <gtest/gtest.h>
+#include <cmath>
+#include <algorithm>
+#include <fstream>
+#include <vector>
+#include "DeepPot.h"
+#include "neighbor_list.h"
+#include "test_utils.h"
+
+#include "google/protobuf/text_format.h"
+#include "google/protobuf/io/zero_copy_stream_impl.h"
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>  
+
+class TestInferDeepPotR : public ::testing::Test
+{  
+protected:  
+  std::vector<double> 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<double> box = {
+    13., 0., 0., 0., 13., 0., 0., 0., 13.
+  };
+  std::vector<double> expected_e = {
+    -9.320909762801588272e+01,-1.868020345400987878e+02,-1.868011172371355997e+02,-9.320868430396934912e+01,-1.868010398844378415e+02,-1.868016706555875999e+02
+  };
+  std::vector<double> expected_f = {
+    6.385312846474267391e-04,-6.460452911141417731e-03,-5.652405655332678417e-04,-7.516468794343579736e-03,1.128804614240160216e-03,5.531937784564192051e-03,1.914138124904981664e-03,5.601819906021693503e-03,-5.131359585752605541e-03,-4.847104424804288617e-03,1.992071550328819614e-03,-4.028159855157302516e-03,1.236340684486603517e-03,-5.373955841338794344e-03,8.312829460571366513e-03,8.574563125108854156e-03,3.111712681889538742e-03,-4.120007238692381148e-03
+  };
+  std::vector<double> expected_v = {
+    5.844056241889131371e-03,4.663973497239899614e-04,-2.268382127762904633e-03,4.663973497239897988e-04,2.349338784202595950e-03,-6.908546513234039253e-04,-2.268382127762904633e-03,-6.908546513234039253e-04,2.040499248150800561e-03,4.238130266437327605e-03,-1.539867187443782223e-04,-2.393101333240631613e-03,-1.539867187443782223e-04,4.410341945447907377e-04,9.544239698119633068e-06,-2.393101333240631613e-03,9.544239698119578858e-06,1.877785959095269654e-03,5.798992562057291543e-03,6.943392552230453693e-04,-1.180376879311998773e-03,6.943392552230453693e-04,1.686725132156275536e-03,-1.461632060145726542e-03,-1.180376879311998556e-03,-1.461632060145726325e-03,1.749543733794208444e-03,7.173915604192910439e-03,3.903218041111061569e-04,-5.747400467123527524e-04,3.903218041111061569e-04,1.208289706621179949e-03,-1.826828914132010932e-03,-5.747400467123527524e-04,-1.826828914132011148e-03,2.856960586657185906e-03,4.067553030177322240e-03,-3.267469855253819430e-05,-6.980667859103454904e-05,-3.267469855253830272e-05,1.387653029234650918e-03,-2.096820720698671855e-03,-6.980667859103444062e-05,-2.096820720698671855e-03,3.218305506720191278e-03,4.753992590355240674e-03,1.224911338353675992e-03,-1.683421934571502484e-03,1.224911338353676209e-03,7.332113564901583539e-04,-1.025577052190138451e-03,-1.683421934571502484e-03,-1.025577052190138234e-03,1.456681925652047018e-03
+  };
+  int natoms;
+  double expected_tot_e;
+  std::vector<double>expected_tot_v;
+
+  deepmd::DeepPot dp;
+
+  void SetUp() override {
+    std::string file_name = "../../tests/infer/deeppot-r.pbtxt";
+    int fd = open(file_name.c_str(), O_RDONLY);
+    tensorflow::protobuf::io::ZeroCopyInputStream* input = new tensorflow::protobuf::io::FileInputStream(fd);
+    tensorflow::GraphDef graph_def;
+    tensorflow::protobuf::TextFormat::Parse(input, &graph_def);
+    delete input;
+    std::fstream output("deeppot.pb", std::ios::out | std::ios::trunc | std::ios::binary);
+    graph_def.SerializeToOstream(&output);
+    // check the string by the following commands
+    // string txt;
+    // tensorflow::protobuf::TextFormat::PrintToString(graph_def, &txt);
+
+    dp.init("deeppot.pb");
+
+    natoms = expected_e.size();
+    EXPECT_EQ(natoms * 3, expected_f.size());
+    EXPECT_EQ(natoms * 9, expected_v.size());
+    expected_tot_e = 0.;
+    expected_tot_v.resize(9);
+    std::fill(expected_tot_v.begin(), expected_tot_v.end(), 0.);
+    for(int ii = 0; ii < natoms; ++ii){
+      expected_tot_e += expected_e[ii];
+    }
+    for(int ii = 0; ii < natoms; ++ii){
+      for(int dd = 0; dd < 9; ++dd){
+	expected_tot_v[dd] += expected_v[ii*9+dd];
+      }
+    }
+  };
+
+  void TearDown() override {
+    remove( "deeppot.pb" ) ;
+  };
+};
+
+
+TEST_F(TestInferDeepPotR, cpu_build_nlist)
+{
+  double ener;
+  std::vector<double> force, virial;
+  dp.compute(ener, force, virial, coord, atype, box);
+
+  EXPECT_EQ(force.size(), natoms*3);
+  EXPECT_EQ(virial.size(), 9);
+
+  EXPECT_LT(fabs(ener - expected_tot_e), 1e-10);
+  for(int ii = 0; ii < natoms*3; ++ii){
+    EXPECT_LT(fabs(force[ii] - expected_f[ii]), 1e-10);    
+  }
+  for(int ii = 0; ii < 3*3; ++ii){
+    EXPECT_LT(fabs(virial[ii] - expected_tot_v[ii]), 1e-10);
+  }
+}
+
+TEST_F(TestInferDeepPotR, cpu_build_nlist_numfv)
+{
+  class MyModel : public EnergyModelTest<double>
+  {
+    deepmd::DeepPot & mydp;
+    const std::vector<int > & atype;
+public:
+    MyModel(
+	deepmd::DeepPot & dp_,
+	const std::vector<int> & atype_
+	) : mydp(dp_), atype(atype_) {};
+    virtual void compute (
+	double & ener,
+	std::vector<double> &	force,
+	std::vector<double> &	virial,
+	const std::vector<double> & coord,
+	const std::vector<double> & box) {
+      mydp.compute(ener, force, virial, coord, atype, box);
+    }
+  };
+  MyModel model(dp, atype);
+  model.test_f(coord, box);
+  model.test_v(coord, box);
+  std::vector<double> box_(box);
+  box_[1] -= 0.4;
+  model.test_f(coord, box_);
+  model.test_v(coord, box_);
+  box_[2] += 0.5;
+  model.test_f(coord, box_);
+  model.test_v(coord, box_);
+  box_[4] += 0.2;
+  model.test_f(coord, box_);
+  model.test_v(coord, box_);
+  box_[3] -= 0.3;
+  model.test_f(coord, box_);
+  model.test_v(coord, box_);
+  box_[6] -= 0.7;
+  model.test_f(coord, box_);
+  model.test_v(coord, box_);
+  box_[7] += 0.6;
+  model.test_f(coord, box_);
+  model.test_v(coord, box_);
+}
+
+TEST_F(TestInferDeepPotR, cpu_build_nlist_atomic)
+{
+  double ener;
+  std::vector<double> force, virial, atom_ener, atom_vir;
+  dp.compute(ener, force, virial, atom_ener, atom_vir, coord, atype, box);
+  
+  EXPECT_EQ(force.size(), natoms*3);
+  EXPECT_EQ(virial.size(), 9);
+  EXPECT_EQ(atom_ener.size(), natoms);
+  EXPECT_EQ(atom_vir.size(), natoms*9);
+
+  EXPECT_LT(fabs(ener - expected_tot_e), 1e-10);
+  for(int ii = 0; ii < natoms*3; ++ii){
+    EXPECT_LT(fabs(force[ii] - expected_f[ii]), 1e-10);    
+  }
+  for(int ii = 0; ii < 3*3; ++ii){
+    EXPECT_LT(fabs(virial[ii] - expected_tot_v[ii]), 1e-10);
+  }
+  for(int ii = 0; ii < natoms; ++ii){
+    EXPECT_LT(fabs(atom_ener[ii] - expected_e[ii]), 1e-10);
+  }
+  for(int ii = 0; ii < natoms*9; ++ii){
+    EXPECT_LT(fabs(atom_vir[ii] - expected_v[ii]), 1e-10);
+  }
+}
+
+
+TEST_F(TestInferDeepPotR, cpu_lmp_nlist)
+{
+  float rc = dp.cutoff();
+  int nloc = coord.size() / 3;  
+  std::vector<double> coord_cpy;
+  std::vector<int> atype_cpy, mapping;  
+  std::vector<std::vector<int > > nlist_data;
+  _build_nlist(nlist_data, coord_cpy, atype_cpy, mapping,
+	       coord, atype, box, rc);
+  int nall = coord_cpy.size() / 3;
+  std::vector<int> ilist(nloc), numneigh(nloc);
+  std::vector<int*> firstneigh(nloc);
+  deepmd::InputNlist inlist(nloc, &ilist[0], &numneigh[0], &firstneigh[0]);
+  convert_nlist(inlist, nlist_data);  
+  
+  double ener;
+  std::vector<double> force_, virial;
+  dp.compute(ener, force_, virial, coord_cpy, atype_cpy, box, nall-nloc, inlist, 0);
+  std::vector<double> force;
+  _fold_back(force, force_, mapping, nloc, nall, 3);
+
+  EXPECT_EQ(force.size(), natoms*3);
+  EXPECT_EQ(virial.size(), 9);
+
+  EXPECT_LT(fabs(ener - expected_tot_e), 1e-10);
+  for(int ii = 0; ii < natoms*3; ++ii){
+    EXPECT_LT(fabs(force[ii] - expected_f[ii]), 1e-10);    
+  }
+  for(int ii = 0; ii < 3*3; ++ii){
+    EXPECT_LT(fabs(virial[ii] - expected_tot_v[ii]), 1e-10);
+  }
+
+  ener = 0.;
+  std::fill(force_.begin(), force_.end(), 0.0);
+  std::fill(virial.begin(), virial.end(), 0.0);
+  dp.compute(ener, force_, virial, coord_cpy, atype_cpy, box, nall-nloc, inlist, 1);
+  _fold_back(force, force_, mapping, nloc, nall, 3);
+
+  EXPECT_EQ(force.size(), natoms*3);
+  EXPECT_EQ(virial.size(), 9);
+
+  EXPECT_LT(fabs(ener - expected_tot_e), 1e-10);
+  for(int ii = 0; ii < natoms*3; ++ii){
+    EXPECT_LT(fabs(force[ii] - expected_f[ii]), 1e-10);    
+  }
+  for(int ii = 0; ii < 3*3; ++ii){
+    EXPECT_LT(fabs(virial[ii] - expected_tot_v[ii]), 1e-10);
+  }
+}
+
+
+TEST_F(TestInferDeepPotR, cpu_lmp_nlist_atomic)
+{
+  float rc = dp.cutoff();
+  int nloc = coord.size() / 3;  
+  std::vector<double> coord_cpy;
+  std::vector<int> atype_cpy, mapping;  
+  std::vector<std::vector<int > > nlist_data;
+  _build_nlist(nlist_data, coord_cpy, atype_cpy, mapping,
+	       coord, atype, box, rc);
+  int nall = coord_cpy.size() / 3;
+  std::vector<int> ilist(nloc), numneigh(nloc);
+  std::vector<int*> firstneigh(nloc);
+  deepmd::InputNlist inlist(nloc, &ilist[0], &numneigh[0], &firstneigh[0]);
+  convert_nlist(inlist, nlist_data);  
+  
+  double ener;
+  std::vector<double> force_, atom_ener_, atom_vir_, virial;
+  std::vector<double> force, atom_ener, atom_vir;
+  dp.compute(ener, force_, virial, atom_ener_, atom_vir_, coord_cpy, atype_cpy, box, nall-nloc, inlist, 0);
+  _fold_back(force, force_, mapping, nloc, nall, 3);
+  _fold_back(atom_ener, atom_ener_, mapping, nloc, nall, 1);
+  _fold_back(atom_vir, atom_vir_, mapping, nloc, nall, 9);
+
+  EXPECT_EQ(force.size(), natoms*3);
+  EXPECT_EQ(virial.size(), 9);
+  EXPECT_EQ(atom_ener.size(), natoms);
+  EXPECT_EQ(atom_vir.size(), natoms*9);
+
+  EXPECT_LT(fabs(ener - expected_tot_e), 1e-10);
+  for(int ii = 0; ii < natoms*3; ++ii){
+    EXPECT_LT(fabs(force[ii] - expected_f[ii]), 1e-10);    
+  }
+  for(int ii = 0; ii < 3*3; ++ii){
+    EXPECT_LT(fabs(virial[ii] - expected_tot_v[ii]), 1e-10);
+  }
+  for(int ii = 0; ii < natoms; ++ii){
+    EXPECT_LT(fabs(atom_ener[ii] - expected_e[ii]), 1e-10);
+  }
+  for(int ii = 0; ii < natoms*9; ++ii){
+    EXPECT_LT(fabs(atom_vir[ii] - expected_v[ii]), 1e-10);
+  }
+
+  ener = 0.;
+  std::fill(force_.begin(), force_.end(), 0.0);
+  std::fill(virial.begin(), virial.end(), 0.0);
+  std::fill(atom_ener_.begin(), atom_ener_.end(), 0.0);
+  std::fill(atom_vir_.begin(), atom_vir_.end(), 0.0);  
+  dp.compute(ener, force_, virial, atom_ener_, atom_vir_, coord_cpy, atype_cpy, box, nall-nloc, inlist, 1);
+  _fold_back(force, force_, mapping, nloc, nall, 3);
+  _fold_back(atom_ener, atom_ener_, mapping, nloc, nall, 1);
+  _fold_back(atom_vir, atom_vir_, mapping, nloc, nall, 9);
+
+  EXPECT_EQ(force.size(), natoms*3);
+  EXPECT_EQ(virial.size(), 9);
+  EXPECT_EQ(atom_ener.size(), natoms);
+  EXPECT_EQ(atom_vir.size(), natoms*9);
+
+  EXPECT_LT(fabs(ener - expected_tot_e), 1e-10);
+  for(int ii = 0; ii < natoms*3; ++ii){
+    EXPECT_LT(fabs(force[ii] - expected_f[ii]), 1e-10);    
+  }
+  for(int ii = 0; ii < 3*3; ++ii){
+    EXPECT_LT(fabs(virial[ii] - expected_tot_v[ii]), 1e-10);
+  }
+  for(int ii = 0; ii < natoms; ++ii){
+    EXPECT_LT(fabs(atom_ener[ii] - expected_e[ii]), 1e-10);
+  }
+  for(int ii = 0; ii < natoms*9; ++ii){
+    EXPECT_LT(fabs(atom_vir[ii] - expected_v[ii]), 1e-10);
+  }
+}
+
+
+TEST_F(TestInferDeepPotR, cpu_lmp_nlist_2rc)
+{
+  float rc = dp.cutoff();
+  int nloc = coord.size() / 3;  
+  std::vector<double> coord_cpy;
+  std::vector<int> atype_cpy, mapping;  
+  std::vector<std::vector<int > > nlist_data;
+  _build_nlist(nlist_data, coord_cpy, atype_cpy, mapping,
+	       coord, atype, box, rc*2);
+  int nall = coord_cpy.size() / 3;
+  std::vector<int> ilist(nloc), numneigh(nloc);
+  std::vector<int*> firstneigh(nloc);
+  deepmd::InputNlist inlist(nloc, &ilist[0], &numneigh[0], &firstneigh[0]);
+  convert_nlist(inlist, nlist_data);  
+  
+  double ener;
+  std::vector<double> force_(nall*3, 0.0), virial(9, 0.0);
+  dp.compute(ener, force_, virial, coord_cpy, atype_cpy, box, nall-nloc, inlist, 0);
+  std::vector<double> force;
+  _fold_back(force, force_, mapping, nloc, nall, 3);
+
+  EXPECT_EQ(force.size(), natoms*3);
+  EXPECT_EQ(virial.size(), 9);
+
+  EXPECT_LT(fabs(ener - expected_tot_e), 1e-10);
+  for(int ii = 0; ii < natoms*3; ++ii){
+    EXPECT_LT(fabs(force[ii] - expected_f[ii]), 1e-10);    
+  }
+  for(int ii = 0; ii < 3*3; ++ii){
+    EXPECT_LT(fabs(virial[ii] - expected_tot_v[ii]), 1e-10);
+  }
+
+  ener = 0.;
+  std::fill(force_.begin(), force_.end(), 0.0);
+  std::fill(virial.begin(), virial.end(), 0.0);
+  dp.compute(ener, force_, virial, coord_cpy, atype_cpy, box, nall-nloc, inlist, 1);
+  _fold_back(force, force_, mapping, nloc, nall, 3);
+
+  EXPECT_EQ(force.size(), natoms*3);
+  EXPECT_EQ(virial.size(), 9);
+
+  EXPECT_LT(fabs(ener - expected_tot_e), 1e-10);
+  for(int ii = 0; ii < natoms*3; ++ii){
+    EXPECT_LT(fabs(force[ii] - expected_f[ii]), 1e-10);    
+  }
+  for(int ii = 0; ii < 3*3; ++ii){
+    EXPECT_LT(fabs(virial[ii] - expected_tot_v[ii]), 1e-10);
+  }
+}
+
+
+TEST_F(TestInferDeepPotR, cpu_lmp_nlist_type_sel)
+{
+  float rc = dp.cutoff();
+
+  // add vir atoms
+  int nvir = 2;
+  std::vector<double> coord_vir(nvir*3);
+  std::vector<int> atype_vir(nvir, 2);
+  for(int ii = 0; ii < nvir; ++ii){
+    coord_vir[ii] = coord[ii];
+  }  
+  coord.insert(coord.begin(), coord_vir.begin(), coord_vir.end());
+  atype.insert(atype.begin(), atype_vir.begin(), atype_vir.end());
+  natoms += nvir;
+  std::vector<double> expected_f_vir(nvir*3, 0.0);
+  expected_f.insert(expected_f.begin(), expected_f_vir.begin(), expected_f_vir.end());
+
+  // build nlist
+  int nloc = coord.size() / 3;  
+  std::vector<double> coord_cpy;
+  std::vector<int> atype_cpy, mapping;  
+  std::vector<std::vector<int > > nlist_data;
+  _build_nlist(nlist_data, coord_cpy, atype_cpy, mapping,
+	       coord, atype, box, rc);
+  int nall = coord_cpy.size() / 3;
+  std::vector<int> ilist(nloc), numneigh(nloc);
+  std::vector<int*> firstneigh(nloc);
+  deepmd::InputNlist inlist(nloc, &ilist[0], &numneigh[0], &firstneigh[0]);
+  convert_nlist(inlist, nlist_data);  
+
+  // dp compute
+  double ener;
+  std::vector<double> force_(nall*3, 0.0), virial(9, 0.0);
+  dp.compute(ener, force_, virial, coord_cpy, atype_cpy, box, nall-nloc, inlist, 0);
+  // fold back
+  std::vector<double> force;
+  _fold_back(force, force_, mapping, nloc, nall, 3);
+
+  EXPECT_EQ(force.size(), natoms*3);
+  EXPECT_EQ(virial.size(), 9);
+
+  EXPECT_LT(fabs(ener - expected_tot_e), 1e-10);
+  for(int ii = 0; ii < natoms*3; ++ii){
+    EXPECT_LT(fabs(force[ii] - expected_f[ii]), 1e-10);    
+  }
+  for(int ii = 0; ii < 3*3; ++ii){
+    EXPECT_LT(fabs(virial[ii] - expected_tot_v[ii]), 1e-10);
+  }
+}
+
+
+
+class TestInferDeepPotRNoPbc : public ::testing::Test
+{  
+protected:  
+  std::vector<double> 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<double> box = {};
+  std::vector<double> expected_e = {
+    -9.321213823508108476e+01,-1.868044102481340758e+02,-1.868067983858651075e+02,-9.320899631301440991e+01,-1.868014559732615112e+02,-1.868017660713088617e+02
+  };
+  std::vector<double> expected_f = {
+    4.578151103701261042e-03,-1.917874111009987628e-03,-3.464546781179331930e-03,-4.578151103701261042e-03,1.917874111009987628e-03,3.464546781179331930e-03,-2.624402581721222913e-03,3.566275128489623933e-04,-2.859315986763691776e-04,-5.767787273464367384e-03,1.907053583551196647e-03,-3.889064429673861831e-03,1.786820066350549132e-04,-5.327197473636275694e-03,8.236236182834734409e-03,8.213507848550535492e-03,3.063516377236116545e-03,-4.061240154484504865e-03
+  };
+  std::vector<double> expected_v = {
+            1.984979026299632174e-03,-8.315452677741701822e-04,-1.502146290172694243e-03,-8.315452677741700738e-04,3.483500446080982317e-04,6.292774999372096039e-04,-1.502146290172694243e-03,6.292774999372097123e-04,1.136759354725281907e-03,1.402852790439301908e-03,-5.876815743732210226e-04,-1.061618327900012114e-03,-5.876815743732211311e-04,2.461909298049979960e-04,4.447320022283834766e-04,-1.061618327900012331e-03,4.447320022283834766e-04,8.033868427351443728e-04,4.143606961846296385e-03,-5.511382161123719835e-04,4.465413399437045397e-04,-5.511382161123719835e-04,1.082271054025323839e-04,-1.097918001262628728e-04,4.465413399437046481e-04,-1.097918001262628728e-04,1.220966982358671871e-04,5.263952004497593831e-03,2.395243710938091842e-04,-2.830378939414603329e-04,2.395243710938094010e-04,1.189969706598244898e-03,-1.805627331015851201e-03,-2.830378939414602245e-04,-1.805627331015851635e-03,2.801996513751836820e-03,2.208413501170402270e-03,5.331756287635716889e-05,-1.664423506603235218e-04,5.331756287635695205e-05,1.379626072862918072e-03,-2.094132943741625064e-03,-1.664423506603234133e-04,-2.094132943741625064e-03,3.199787996743366607e-03,4.047014004814953811e-03,1.137904999421357000e-03,-1.568106936614101698e-03,1.137904999421357217e-03,7.205982843216952307e-04,-1.011174600268313238e-03,-1.568106936614101698e-03,-1.011174600268313238e-03,1.435226522157425754e-03
+  };
+  int natoms;
+  double expected_tot_e;
+  std::vector<double>expected_tot_v;
+
+  deepmd::DeepPot dp;
+
+  void SetUp() override {
+    std::string file_name = "../../tests/infer/deeppot-r.pbtxt";
+    int fd = open(file_name.c_str(), O_RDONLY);
+    tensorflow::protobuf::io::ZeroCopyInputStream* input = new tensorflow::protobuf::io::FileInputStream(fd);
+    tensorflow::GraphDef graph_def;
+    tensorflow::protobuf::TextFormat::Parse(input, &graph_def);
+    delete input;
+    std::fstream output("deeppot.pb", std::ios::out | std::ios::trunc | std::ios::binary);
+    graph_def.SerializeToOstream(&output);
+
+    dp.init("deeppot.pb");
+
+    natoms = expected_e.size();
+    EXPECT_EQ(natoms * 3, expected_f.size());
+    EXPECT_EQ(natoms * 9, expected_v.size());
+    expected_tot_e = 0.;
+    expected_tot_v.resize(9);
+    std::fill(expected_tot_v.begin(), expected_tot_v.end(), 0.);
+    for(int ii = 0; ii < natoms; ++ii){
+      expected_tot_e += expected_e[ii];
+    }
+    for(int ii = 0; ii < natoms; ++ii){
+      for(int dd = 0; dd < 9; ++dd){
+	expected_tot_v[dd] += expected_v[ii*9+dd];
+      }
+    }
+  };
+
+  void TearDown() override {
+    remove( "deeppot.pb" ) ;
+  };
+};
+
+TEST_F(TestInferDeepPotRNoPbc, cpu_build_nlist)
+{
+  double ener;
+  std::vector<double> force, virial;
+  dp.compute(ener, force, virial, coord, atype, box);
+
+  EXPECT_EQ(force.size(), natoms*3);
+  EXPECT_EQ(virial.size(), 9);
+
+  EXPECT_LT(fabs(ener - expected_tot_e), 1e-10);
+  for(int ii = 0; ii < natoms*3; ++ii){
+    EXPECT_LT(fabs(force[ii] - expected_f[ii]), 1e-10);    
+  }
+  for(int ii = 0; ii < 3*3; ++ii){
+    EXPECT_LT(fabs(virial[ii] - expected_tot_v[ii]), 1e-10);
+  }
+}
diff --git a/source/api_cc/tests/test_dipolecharge.cc b/source/api_cc/tests/test_dipolecharge.cc
new file mode 100644
index 0000000000..b13c3c7d7a
--- /dev/null
+++ b/source/api_cc/tests/test_dipolecharge.cc
@@ -0,0 +1,240 @@
+#include <gtest/gtest.h>
+#include <cmath>
+#include <algorithm>
+#include <fstream>
+#include <vector>
+#include "DeepTensor.h"
+#include "DataModifier.h"
+#include "SimulationRegion.h"
+#include "ewald.h"
+#include "neighbor_list.h"
+#include "test_utils.h"
+
+#include "google/protobuf/text_format.h"
+#include "google/protobuf/io/zero_copy_stream_impl.h"
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>  
+
+class TestDipoleCharge : public ::testing::Test
+{  
+protected:  
+  std::vector<double> coord = {
+    4.6067455554,    8.8719311819,    6.3886531197,
+    4.0044515745,    4.2449530507,    7.7902855220,
+    2.6453069446,    0.8772647726,    1.2804446790,
+    1.1445332290,    0.0067366438,    1.8606485070,
+    7.1002867706,    5.0325506787,    3.1805888348,
+    4.5352891138,    7.7389683929,    9.4260970128,
+    2.1833238914,    9.0916071034,    7.2299906064,
+    4.1040157820,    1.0496745045,    5.4748315591,
+  };
+  std::vector<int> atype = {
+    0,3,2,1,3,4,1,4
+  };
+  std::vector<double> box = {
+    10., 0., 0., 0., 10., 0., 0., 0., 10.
+  };
+  std::vector<double> expected_e = {
+    3.671081837126222158e+00
+  };
+  std::vector<double> expected_f = {
+    8.786854427753210128e-01,-1.590752486903602159e-01,-2.709225006303785932e-01,-4.449513960033193438e-01,-1.564291540964127813e-01,2.139031741772115178e-02,1.219699614140521193e+00,-5.580358618499958734e-02,-3.878662478349682585e-01,-1.286685244990778854e+00,1.886475802950296488e-01,3.904450515493615437e-01,1.605017382138404849e-02,2.138016869742287995e-01,-2.617514921203008965e-02,2.877081057057793712e-01,-3.846449683844421763e-01,3.048855616906603894e-02,-9.075632811311897807e-01,-6.509653472431625731e-03,2.302010972126376787e-01,2.370565856822822726e-01,3.600133435593881881e-01,1.243887532859055609e-02
+  };
+  std::vector<double> expected_v = {
+    3.714071471995848417e-01,6.957130186032146613e-01,-1.158289779017217302e+00,6.957130186032139951e-01,-1.400130091653774933e+01,-3.631620234653316626e-01,-1.158289779017217302e+00,-3.631620234653316626e-01,3.805077486043773050e+00
+  };
+  std::vector<double> charge_map = {
+    1., 1., 1., 1., 1., -1., -3.
+  };
+  int natoms;
+  int ntypes;
+  std::vector<int> type_asso;
+  double expected_tot_e;
+  std::vector<double>expected_tot_v;
+
+  deepmd::DeepTensor dp;
+  deepmd::DipoleChargeModifier dm;
+
+  void SetUp() override {
+    std::string file_name = "../../tests/infer/dipolecharge_e.pbtxt";
+    int fd = open(file_name.c_str(), O_RDONLY);
+    tensorflow::protobuf::io::ZeroCopyInputStream* input = new tensorflow::protobuf::io::FileInputStream(fd);
+    tensorflow::GraphDef graph_def;
+    tensorflow::protobuf::TextFormat::Parse(input, &graph_def);
+    delete input;
+    std::string model = "dipolecharge_e.pb";
+    std::fstream output(model.c_str(), std::ios::out | std::ios::trunc | std::ios::binary);
+    graph_def.SerializeToOstream(&output);
+    // check the string by the following commands
+    // string txt;
+    // tensorflow::protobuf::TextFormat::PrintToString(graph_def, &txt);
+
+    // dp.init("dipolecharge_d.pb");
+    // dm.init("dipolecharge_d.pb");
+    dp.init(model, 0, "dipole_charge");
+    dm.init(model, 0, "dipole_charge");
+
+    natoms = atype.size();
+    ntypes = 5;
+    type_asso.resize(ntypes, -1);
+    type_asso[1] = 5;
+    type_asso[3] = 6;
+
+    EXPECT_EQ(natoms * 3, expected_f.size());
+    EXPECT_EQ(9, expected_v.size());
+  };
+
+  void TearDown() override {
+    remove( "dipolecharge_e.pb" ) ;
+  };
+};
+
+static bool
+_in_vec(const int & value,
+	const std::vector<int> & vec)
+{
+  // naive impl.
+  for(int ii = 0; ii < vec.size(); ++ii){
+    if(value == vec[ii]) return true;
+  }
+  return false;
+}
+
+TEST_F(TestDipoleCharge, cpu_lmp_nlist)
+{
+  // build nlist
+  // float rc = dp.cutoff();
+  float rc = 4.0;
+  int nloc = coord.size() / 3;  
+  std::vector<double> coord_cpy;
+  std::vector<int> atype_cpy, mapping;  
+  std::vector<std::vector<int > > nlist_data;
+  _build_nlist(nlist_data, coord_cpy, atype_cpy, mapping,
+  	       coord, atype, box, rc);
+  int nall = coord_cpy.size() / 3;
+  int nghost = nall - nloc;
+  std::vector<int> ilist(nloc), numneigh(nloc);
+  std::vector<int*> firstneigh(nloc);
+  deepmd::InputNlist inlist(nloc, &ilist[0], &numneigh[0], &firstneigh[0]);
+  convert_nlist(inlist, nlist_data);  
+
+  // evaluate dipole
+  std::vector<double> dipole, dipole_recd(nloc*3, 0.0);
+  dp.compute(dipole, coord_cpy, atype_cpy, box, nall-nloc, inlist);
+
+  // add virtual atoms to the system
+  // // a lot of mappings
+  std::vector<int> sel_types = dp.sel_types();
+  std::vector<int> sel_fwd, sel_bwd;
+  int sel_nghost;
+  deepmd::select_by_type(sel_fwd, sel_bwd, sel_nghost, coord_cpy, atype_cpy, nghost, sel_types);
+  int sel_nall = sel_bwd.size();
+  int sel_nloc = sel_nall - sel_nghost;
+  std::vector<int> sel_atype(sel_bwd.size());
+  deepmd::select_map<int>(sel_atype, atype, sel_fwd, 1);
+  deepmd::AtomMap<double> nnp_map(sel_atype.begin(), sel_atype.begin() + sel_nloc);
+  const std::vector<int> & sort_fwd_map(nnp_map.get_fwd_map());
+
+  // // add coords
+  std::vector<double > add_coord;
+  std::vector<int > add_atype;
+  std::vector<std::pair<int,int>> pairs;
+  for(int ii = 0; ii < nloc; ++ii){
+    if(_in_vec(atype[ii], sel_types)){
+      int res_idx = sort_fwd_map[sel_fwd[ii]];
+      std::vector<double > tmp_coord(3);
+      for(int dd = 0; dd < 3; ++dd){
+	tmp_coord[dd] = coord[ii*3+dd] + dipole[res_idx*3+dd];
+	dipole_recd[ii*3+dd] = dipole[res_idx*3+dd];
+      }
+      pairs.push_back(std::pair<int,int>(ii, add_atype.size()+atype.size()));
+      // std::cout << ii <<  " " 
+      // 		<< atype[ii] << " " 
+      // 		<< res_idx << " " 
+      // 		<< type_asso[atype[ii]] << " " 
+      // 		<< " pair "  
+      // 		<< pairs.back().first << " " << pairs.back().second << " "
+      // 		<< std::endl;
+      add_coord.insert(add_coord.end(), tmp_coord.begin(), tmp_coord.end());
+      add_atype.push_back(type_asso[atype[ii]]);      
+    }
+  }
+  coord.insert(coord.end(), add_coord.begin(), add_coord.end());
+  atype.insert(atype.end(), add_atype.begin(), add_atype.end());
+  nloc = atype.size();
+  EXPECT_EQ(atype.size()*3, coord.size());
+
+  // get charge value
+  std::vector<double> charge(nloc);
+  for(int ii = 0; ii < nloc; ++ii){
+    charge[ii] = charge_map[atype[ii]];
+  }
+  
+  // compute the recp part of the ele interaction
+  double eener;
+  std::vector<double> eforce, evirial;
+  deepmd::Region<double> region;
+  init_region_cpu(region, &box[0]);
+  deepmd::EwaldParameters<double> eparam;
+  eparam.beta = 0.2;
+  eparam.spacing = 4;
+  ewald_recp(eener, eforce, evirial, coord, charge, region, eparam);
+  
+  EXPECT_LT(fabs(eener - expected_e[0]), 1e-6);
+  EXPECT_EQ(eforce.size(), coord.size());
+  EXPECT_EQ(evirial.size(), 9);  
+
+  // extend the system with virtual atoms, and build nlist
+  _build_nlist(nlist_data, coord_cpy, atype_cpy, mapping,
+  	       coord, atype, box, rc);
+  nall = coord_cpy.size() / 3;
+  nghost = nall - nloc;
+  ilist.resize(nloc);
+  numneigh.resize(nloc);
+  firstneigh.resize(nloc);
+  inlist.inum = nloc;
+  inlist.ilist = &ilist[0];
+  inlist.numneigh = &numneigh[0];
+  inlist.firstneigh = &firstneigh[0];
+  convert_nlist(inlist, nlist_data);
+
+  // compute force and virial
+  std::vector<double > force_, force, virial;
+  dm.compute(force_, virial, coord_cpy, atype_cpy, box, pairs, eforce, nghost, inlist);
+  // for(int ii = 0; ii < force_.size(); ++ii){
+  //   std::cout << force_[ii] << " " ;
+  // }
+  // std::cout << std::endl;
+  _fold_back(force, force_, mapping, nloc, nall, 3);
+
+  // compare force
+  EXPECT_EQ(force.size(), nloc*3);
+  // note nloc > expected_f.size(), because nloc contains virtual atoms.
+  for(int ii = 0; ii < expected_f.size(); ++ii){
+    EXPECT_LT(fabs(force[ii] - expected_f[ii]), 1e-6);
+  }
+
+  // add recp virial and viral corr to virial
+  // virial = virial_recp + virial_dipolecharge + virial_corr
+  for (int dd0 = 0; dd0 < 3; ++dd0){
+    for (int dd1 = 0; dd1 < 3; ++dd1){
+      virial[dd0*3+dd1] += evirial[dd0*3+dd1];
+    }
+  }    
+  for(int ii = 0; ii < pairs.size(); ++ii){
+    int idx0 = pairs[ii].first;
+    int idx1 = pairs[ii].second;
+    for (int dd0 = 0; dd0 < 3; ++dd0){
+      for (int dd1 = 0; dd1 < 3; ++dd1){
+	virial[dd0*3+dd1] -= eforce[idx1*3+dd0] * dipole_recd[idx0*3+dd1];
+      }
+    }    
+  }
+  // compare virial
+  EXPECT_EQ(virial.size(), 3*3);
+  for(int ii = 0; ii < expected_v.size(); ++ii){
+    EXPECT_LT(fabs(virial[ii] - expected_v[ii]), 1e-5);
+  }
+}
+
diff --git a/source/api_cc/tests/test_ewald.cc b/source/api_cc/tests/test_ewald.cc
new file mode 100644
index 0000000000..e27b1087e3
--- /dev/null
+++ b/source/api_cc/tests/test_ewald.cc
@@ -0,0 +1,78 @@
+#include <gtest/gtest.h>
+#include <cmath>
+#include <algorithm>
+#include <fstream>
+#include <vector>
+#include "neighbor_list.h"
+#include "test_utils.h"
+#include "ewald.h"
+
+class TestInferEwald : public ::testing::Test
+{  
+protected:  
+  std::vector<double> 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<double> charge = {
+    -2, 1, 1, -2, 1, 1
+  };
+  std::vector<double> box = {
+    13., 0., 0., 0., 13., 0., 0., 0., 13.
+  };
+  void SetUp() override {
+  };
+  void TearDown() override {
+  };
+};
+
+TEST_F(TestInferEwald, cpu_numfv)
+{
+  class MyModel : public EnergyModelTest<double>
+  {
+    const std::vector<double > & charge;
+    deepmd::EwaldParameters<double> eparam;    
+public:
+    MyModel(
+	const std::vector<double> & charge_
+	) : charge(charge_) {
+      eparam.beta = 0.4;
+    };
+    virtual void compute (
+	double & ener,
+	std::vector<double> &	force,
+	std::vector<double> &	virial,
+	const std::vector<double> & coord,
+	const std::vector<double> & box) {
+      deepmd::Region<double> region;
+      init_region_cpu(region, &box[0]);
+      ewald_recp(ener, force, virial, coord, charge, region, eparam);
+    }
+  };
+  MyModel model(charge);
+  model.test_f(coord, box);
+  model.test_v(coord, box);
+  std::vector<double> box_(box);
+  box_[1] -= 0.2;
+  model.test_f(coord, box_);
+  model.test_v(coord, box_);
+  box_[2] += 0.5;
+  model.test_f(coord, box_);
+  model.test_v(coord, box_);
+  box_[4] += 0.2;
+  model.test_f(coord, box_);
+  model.test_v(coord, box_);
+  box_[3] -= 0.3;
+  model.test_f(coord, box_);
+  model.test_v(coord, box_);
+  box_[6] -= 0.7;
+  model.test_f(coord, box_);
+  model.test_v(coord, box_);
+  box_[7] += 0.6;
+  model.test_f(coord, box_);
+  model.test_v(coord, box_);
+}
diff --git a/source/api_cc/tests/test_utils.h b/source/api_cc/tests/test_utils.h
new file mode 100644
index 0000000000..e6a496e374
--- /dev/null
+++ b/source/api_cc/tests/test_utils.h
@@ -0,0 +1,135 @@
+#pragma once
+#include <cmath>
+
+inline void 
+_fold_back(
+    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.begin(), in.begin() + nloc*ndim, out.begin());
+  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];
+    }
+  }
+}
+
+inline void
+_build_nlist(
+    std::vector<std::vector<int>> &nlist_data,
+    std::vector<double > & coord_cpy,
+    std::vector<int > & atype_cpy,
+    std::vector<int > & mapping,
+    const std::vector<double > & coord,
+    const std::vector<int > & atype,
+    const std::vector<double > & box,
+    const float & rc)
+{
+  SimulationRegion<double > region;
+  region.reinitBox(&box[0]);
+  std::vector<int> ncell, ngcell;
+  copy_coord(coord_cpy, atype_cpy, mapping, ncell, ngcell, coord, atype, rc, region);
+  std::vector<int> nat_stt, ext_stt, ext_end;
+  nat_stt.resize(3);
+  ext_stt.resize(3);
+  ext_end.resize(3);
+  for (int dd = 0; dd < 3; ++dd){
+    ext_stt[dd] = -ngcell[dd];
+    ext_end[dd] = ncell[dd] + ngcell[dd];
+  }
+  int nloc = coord.size() / 3;
+  int nall = coord_cpy.size() / 3;
+  std::vector<std::vector<int>> nlist_r_cpy;
+  build_nlist(nlist_data, nlist_r_cpy, coord_cpy, nloc, rc, rc, nat_stt, ncell, ext_stt, ext_end, region, ncell);
+}
+
+template<typename VALUETYPE>
+class EnergyModelTest
+{
+  double hh = 1e-5;
+  double level = 1e-6;
+public:
+  virtual void compute (
+      VALUETYPE & ener,
+      std::vector<VALUETYPE> &	force,
+      std::vector<VALUETYPE> &	virial,
+      const std::vector<VALUETYPE> & coord,
+      const std::vector<VALUETYPE> & box
+      ) = 0;
+  void test_f (
+      const std::vector<VALUETYPE> & coord,
+      const std::vector<VALUETYPE> & box) {
+    int ndof = coord.size();
+    VALUETYPE ener;
+    std::vector<VALUETYPE> force, virial;
+    compute(ener, force, virial, coord, box);
+    for(int ii = 0; ii < ndof; ++ii){
+      std::vector<VALUETYPE> coord0(coord), coord1(coord);
+      VALUETYPE ener0, ener1;
+      std::vector<VALUETYPE> forcet, virialt;
+      coord0[ii] += hh;
+      coord1[ii] -= hh;
+      compute(ener0, forcet, virialt, coord0, box);
+      compute(ener1, forcet, virialt, coord1, box);
+      VALUETYPE num = - (ener0 - ener1) / (2.*hh);
+      VALUETYPE ana = force[ii];
+      EXPECT_LT(fabs(num - ana), level);
+    }
+  }  
+  void test_v(
+      const std::vector<VALUETYPE> & coord,
+      const std::vector<VALUETYPE> & box) {
+    std::vector<VALUETYPE> num_diff(9);
+    VALUETYPE ener;
+    std::vector<VALUETYPE> force, virial;
+    compute(ener, force, virial, coord, box);
+    deepmd::Region<VALUETYPE> region;
+    init_region_cpu(region, &box[0]);
+    for(int ii = 0; ii < 9; ++ii){
+      std::vector<VALUETYPE> box0(box), box1(box);
+      box0[ii] += hh;
+      box1[ii] -= hh;
+      deepmd::Region<VALUETYPE> region0, region1;
+      init_region_cpu(region0, &box0[0]);
+      init_region_cpu(region1, &box1[0]);
+      std::vector<VALUETYPE> coord0(coord), coord1(coord);
+      int natoms = coord.size() / 3;
+      for(int ii = 0; ii < natoms; ++ii){
+	VALUETYPE pi[3];
+	convert_to_inter_cpu(pi, region, &coord[ii*3]);
+	convert_to_phys_cpu(&coord0[ii*3], region0, pi);
+      }
+      for(int ii = 0; ii < natoms; ++ii){
+	VALUETYPE pi[3];
+	convert_to_inter_cpu(pi, region, &coord[ii*3]);
+	convert_to_phys_cpu(&coord1[ii*3], region1, pi);
+      }
+      VALUETYPE ener0, ener1;
+      std::vector<VALUETYPE> forcet, virialt;
+      compute(ener0, forcet, virialt, coord0, box0);
+      compute(ener1, forcet, virialt, coord1, box1);
+      num_diff[ii] = - (ener0 - ener1) / (2.*hh);
+    }
+    std::vector<VALUETYPE> num_virial(9, 0);
+    for(int dd0 = 0; dd0 < 3; ++dd0){
+      for(int dd1 = 0; dd1 < 3; ++dd1){
+	for(int dd = 0; dd < 3; ++dd){
+	  num_virial[dd0*3+dd1] += num_diff[dd*3+dd0] * box[dd*3+dd1];
+	  // num_virial[dd0*3+dd1] += num_diff[dd0*3+dd] * box[dd1*3+dd];
+	}
+      }
+    }
+    for(int ii = 0; ii < 9; ++ii){
+      EXPECT_LT(fabs(num_virial[ii] - virial[ii]), level);
+    }
+  }
+};
+
+