diff --git a/.gitignore b/.gitignore
index 9f63a65219..69f82afbb5 100644
--- a/.gitignore
+++ b/.gitignore
@@ -28,7 +28,7 @@ venv*
 .vscode/**
 _build
 _templates
-API_CC
+doc/API_CC/
 doc/api_py/
 doc/api_core/
 doc/api_c/
diff --git a/deepmd/pt/model/task/fitting.py b/deepmd/pt/model/task/fitting.py
index 22bbf6165b..490d739f59 100644
--- a/deepmd/pt/model/task/fitting.py
+++ b/deepmd/pt/model/task/fitting.py
@@ -657,3 +657,8 @@ def _forward_common(
         outs = torch.where(mask[:, :, None], outs, 0.0)
         results.update({self.var_name: outs})
         return results
+
+    @torch.jit.export
+    def get_task_dim(self) -> int:
+        """Get the output dimension of the fitting net."""
+        return self._net_out_dim()
diff --git a/source/api_cc/include/DeepTensorPT.h b/source/api_cc/include/DeepTensorPT.h
new file mode 100644
index 0000000000..c602fc53e0
--- /dev/null
+++ b/source/api_cc/include/DeepTensorPT.h
@@ -0,0 +1,247 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+#pragma once
+
+#include <torch/script.h>
+#include <torch/torch.h>
+
+#include "DeepTensor.h"
+
+namespace deepmd {
+/**
+ * @brief PyTorch implementation for Deep Tensor.
+ **/
+class DeepTensorPT : public DeepTensorBase {
+ public:
+  /**
+   * @brief Deep Tensor constructor without initialization.
+   **/
+  DeepTensorPT();
+  virtual ~DeepTensorPT();
+  /**
+   * @brief Deep Tensor constructor with 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 Name scopes of operations.
+   **/
+  DeepTensorPT(const std::string& model,
+               const int& gpu_rank = 0,
+               const std::string& name_scope = "");
+  /**
+   * @brief Initialize the Deep Tensor.
+   * @param[in] model The name of the frozen model file.
+   * @param[in] gpu_rank The GPU rank. Default is 0.
+   * @param[in] name_scope Name scopes of operations.
+   **/
+  void init(const std::string& model,
+            const int& gpu_rank = 0,
+            const std::string& name_scope = "");
+
+ private:
+  /**
+   * @brief Evaluate the global tensor and component-wise force and virial.
+   * @param[out] global_tensor The global tensor to evaluate.
+   * @param[out] force The component-wise force of the global tensor, size odim
+   *x natoms x 3.
+   * @param[out] virial The component-wise virial of the global tensor, size
+   *odim x 9.
+   * @param[out] atom_tensor The atomic tensor value of the model, size natoms x
+   *odim.
+   * @param[out] atom_virial The component-wise atomic virial of the global
+   *tensor, size odim x natoms x 9.
+   * @param[in] coord The coordinates of atoms. The array should be of size
+   *natoms x 3.
+   * @param[in] atype The atom types. The list should contain natoms ints.
+   * @param[in] box The cell of the region. The array should be of size 9.
+   * @param[in] request_deriv Whether to request the derivative of the global
+   * tensor, including force and virial.
+   **/
+  template <typename VALUETYPE>
+  void compute(std::vector<VALUETYPE>& global_tensor,
+               std::vector<VALUETYPE>& force,
+               std::vector<VALUETYPE>& virial,
+               std::vector<VALUETYPE>& atom_tensor,
+               std::vector<VALUETYPE>& atom_virial,
+               const std::vector<VALUETYPE>& coord,
+               const std::vector<int>& atype,
+               const std::vector<VALUETYPE>& box,
+               const bool request_deriv);
+  /**
+   * @brief Evaluate the global tensor and component-wise force and virial.
+   * @param[out] global_tensor The global tensor to evaluate.
+   * @param[out] force The component-wise force of the global tensor, size odim
+   *x natoms x 3.
+   * @param[out] virial The component-wise virial of the global tensor, size
+   *odim x 9.
+   * @param[out] atom_tensor The atomic tensor value of the model, size natoms x
+   *odim.
+   * @param[out] atom_virial The component-wise atomic virial of the global
+   *tensor, size odim x natoms x 9.
+   * @param[in] coord The coordinates of atoms. The array should be of size
+   *natoms x 3.
+   * @param[in] atype The atom types. The list should contain natoms ints.
+   * @param[in] box The cell of the region. The array should be of size 9.
+   * @param[in] nghost The number of ghost atoms.
+   * @param[in] inlist The input neighbour list.
+   * @param[in] request_deriv Whether to request the derivative of the global
+   * tensor, including force and virial.
+   **/
+  template <typename VALUETYPE>
+  void compute(std::vector<VALUETYPE>& global_tensor,
+               std::vector<VALUETYPE>& force,
+               std::vector<VALUETYPE>& virial,
+               std::vector<VALUETYPE>& atom_tensor,
+               std::vector<VALUETYPE>& atom_virial,
+               const std::vector<VALUETYPE>& coord,
+               const std::vector<int>& atype,
+               const std::vector<VALUETYPE>& box,
+               const int nghost,
+               const InputNlist& inlist,
+               const bool request_deriv);
+
+ public:
+  /**
+   * @brief Get the cutoff radius.
+   * @return The cutoff radius.
+   **/
+  double cutoff() const {
+    assert(inited);
+    return rcut;
+  };
+  /**
+   * @brief Get the number of types.
+   * @return The number of types.
+   **/
+  int numb_types() const {
+    assert(inited);
+    return ntypes;
+  };
+  /**
+   * @brief Get the output dimension.
+   * @return The output dimension.
+   **/
+  int output_dim() const {
+    assert(inited);
+    return odim;
+  };
+  /**
+   * @brief Get the list of sel types.
+   * @return The list of sel types.
+   */
+  const std::vector<int>& sel_types() const {
+    assert(inited);
+    return sel_type;
+  };
+  /**
+   * @brief Get the type map (element name of the atom types) of this model.
+   * @param[out] type_map The type map of this model.
+   **/
+  void get_type_map(std::string& type_map);
+
+  /**
+   * @brief Evaluate the global tensor and component-wise force and virial.
+   * @param[out] global_tensor The global tensor to evaluate.
+   * @param[out] force The component-wise force of the global tensor, size odim
+   *x natoms x 3.
+   * @param[out] virial The component-wise virial of the global tensor, size
+   *odim x 9.
+   * @param[out] atom_tensor The atomic tensor value of the model, size natoms x
+   *odim.
+   * @param[out] atom_virial The component-wise atomic virial of the global
+   *tensor, size odim x natoms x 9.
+   * @param[in] coord The coordinates of atoms. The array should be of size
+   *natoms x 3.
+   * @param[in] atype The atom types. The list should contain natoms ints.
+   * @param[in] box The cell of the region. The array should be of size 9.
+   * @param[in] request_deriv Whether to request the derivative of the global
+   * tensor, including force and virial.
+   * @{
+   **/
+  void computew(std::vector<double>& global_tensor,
+                std::vector<double>& force,
+                std::vector<double>& virial,
+                std::vector<double>& atom_tensor,
+                std::vector<double>& atom_virial,
+                const std::vector<double>& coord,
+                const std::vector<int>& atype,
+                const std::vector<double>& box,
+                const bool request_deriv);
+  void computew(std::vector<float>& global_tensor,
+                std::vector<float>& force,
+                std::vector<float>& virial,
+                std::vector<float>& atom_tensor,
+                std::vector<float>& atom_virial,
+                const std::vector<float>& coord,
+                const std::vector<int>& atype,
+                const std::vector<float>& box,
+                const bool request_deriv);
+  /** @} */
+  /**
+   * @brief Evaluate the global tensor and component-wise force and virial.
+   * @param[out] global_tensor The global tensor to evaluate.
+   * @param[out] force The component-wise force of the global tensor, size odim
+   *x natoms x 3.
+   * @param[out] virial The component-wise virial of the global tensor, size
+   *odim x 9.
+   * @param[out] atom_tensor The atomic tensor value of the model, size natoms x
+   *odim.
+   * @param[out] atom_virial The component-wise atomic virial of the global
+   *tensor, size odim x natoms x 9.
+   * @param[in] coord The coordinates of atoms. The array should be of size
+   *natoms x 3.
+   * @param[in] atype The atom types. The list should contain natoms ints.
+   * @param[in] box The cell of the region. The array should be of size 9.
+   * @param[in] nghost The number of ghost atoms.
+   * @param[in] inlist The input neighbour list.
+   * @param[in] request_deriv Whether to request the derivative of the global
+   * tensor, including force and virial.
+   * @{
+   **/
+  void computew(std::vector<double>& global_tensor,
+                std::vector<double>& force,
+                std::vector<double>& virial,
+                std::vector<double>& atom_tensor,
+                std::vector<double>& atom_virial,
+                const std::vector<double>& coord,
+                const std::vector<int>& atype,
+                const std::vector<double>& box,
+                const int nghost,
+                const InputNlist& inlist,
+                const bool request_deriv);
+  void computew(std::vector<float>& global_tensor,
+                std::vector<float>& force,
+                std::vector<float>& virial,
+                std::vector<float>& atom_tensor,
+                std::vector<float>& atom_virial,
+                const std::vector<float>& coord,
+                const std::vector<int>& atype,
+                const std::vector<float>& box,
+                const int nghost,
+                const InputNlist& inlist,
+                const bool request_deriv);
+  /** @} */
+
+ private:
+  int num_intra_nthreads, num_inter_nthreads;
+  bool inited;
+  double rcut;
+  int ntypes;
+  mutable int odim;
+  std::vector<int> sel_type;
+  std::string name_scope;
+  // PyTorch module and device management
+  mutable torch::jit::script::Module module;
+  int gpu_id;
+  bool gpu_enabled;
+  NeighborListData nlist_data;
+  // Neighbor list tensors for efficient computation
+  at::Tensor firstneigh_tensor;
+
+  /**
+   * @brief Translate PyTorch exceptions to the DeePMD-kit exception.
+   * @param[in] f The function to run.
+   * @example translate_error([&](){...});
+   */
+  void translate_error(std::function<void()> f);
+};
+
+}  // namespace deepmd
diff --git a/source/api_cc/src/DeepTensor.cc b/source/api_cc/src/DeepTensor.cc
index a9031472e6..a40a16b79c 100644
--- a/source/api_cc/src/DeepTensor.cc
+++ b/source/api_cc/src/DeepTensor.cc
@@ -6,6 +6,9 @@
 #ifdef BUILD_TENSORFLOW
 #include "DeepTensorTF.h"
 #endif
+#ifdef BUILD_PYTORCH
+#include "DeepTensorPT.h"
+#endif
 #include "common.h"
 
 using namespace deepmd;
@@ -38,7 +41,11 @@ void DeepTensor::init(const std::string &model,
     throw deepmd::deepmd_exception("TensorFlow backend is not built.");
 #endif
   } else if (deepmd::DPBackend::PyTorch == backend) {
-    throw deepmd::deepmd_exception("PyTorch backend is not supported yet");
+#ifdef BUILD_PYTORCH
+    dt = std::make_shared<deepmd::DeepTensorPT>(model, gpu_rank, name_scope_);
+#else
+    throw deepmd::deepmd_exception("PyTorch backend is not built.");
+#endif
   } else if (deepmd::DPBackend::Paddle == backend) {
     throw deepmd::deepmd_exception("PaddlePaddle backend is not supported yet");
   } else {
diff --git a/source/api_cc/src/DeepTensorPT.cc b/source/api_cc/src/DeepTensorPT.cc
new file mode 100644
index 0000000000..1636f3af95
--- /dev/null
+++ b/source/api_cc/src/DeepTensorPT.cc
@@ -0,0 +1,482 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+#ifdef BUILD_PYTORCH
+#include "DeepTensorPT.h"
+
+#include <torch/csrc/jit/runtime/jit_exception.h>
+
+#include <cstdint>
+#include <numeric>  // for std::iota
+#include <sstream>
+
+#include "common.h"
+#include "device.h"
+#include "errors.h"
+
+using namespace deepmd;
+
+static torch::Tensor createNlistTensor(
+    const std::vector<std::vector<int>>& data) {
+  size_t total_size = 0;
+  for (const auto& row : data) {
+    total_size += row.size();
+  }
+  std::vector<int> flat_data;
+  flat_data.reserve(total_size);
+  for (const auto& row : data) {
+    flat_data.insert(flat_data.end(), row.begin(), row.end());
+  }
+
+  torch::Tensor flat_tensor = torch::tensor(flat_data, torch::kInt32);
+  int nloc = data.size();
+  int nnei = nloc > 0 ? total_size / nloc : 0;
+  return flat_tensor.view({1, nloc, nnei});
+}
+
+void DeepTensorPT::translate_error(std::function<void()> f) {
+  try {
+    f();
+    // it seems that libtorch may throw different types of exceptions which are
+    // inherbited from different base classes
+    // https://github.com/pytorch/pytorch/blob/13316a8d4642454012d34da0d742f1ba93fc0667/torch/csrc/jit/runtime/interpreter.cpp#L924-L939
+  } catch (const c10::Error& e) {
+    throw deepmd::deepmd_exception("DeePMD-kit PyTorch backend error: " +
+                                   std::string(e.what()));
+  } catch (const torch::jit::JITException& e) {
+    throw deepmd::deepmd_exception("DeePMD-kit PyTorch backend JIT error: " +
+                                   std::string(e.what()));
+  } catch (const std::runtime_error& e) {
+    throw deepmd::deepmd_exception("DeePMD-kit PyTorch backend error: " +
+                                   std::string(e.what()));
+  }
+}
+
+DeepTensorPT::DeepTensorPT() : inited(false) {}
+
+DeepTensorPT::DeepTensorPT(const std::string& model,
+                           const int& gpu_rank,
+                           const std::string& name_scope_)
+    : inited(false), name_scope(name_scope_) {
+  try {
+    translate_error([&] { init(model, gpu_rank, name_scope_); });
+  } catch (...) {
+    // Clean up and rethrow, as the destructor will not be called
+    throw;
+  }
+}
+
+void DeepTensorPT::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_;
+  deepmd::load_op_library();
+  int gpu_num = torch::cuda::device_count();
+  if (gpu_num > 0) {
+    gpu_id = gpu_rank % gpu_num;
+  } else {
+    gpu_id = 0;
+  }
+  torch::Device device(torch::kCUDA, gpu_id);
+  gpu_enabled = torch::cuda::is_available();
+  if (!gpu_enabled) {
+    device = torch::Device(torch::kCPU);
+    std::cout << "load model from: " << model << " to cpu " << std::endl;
+  } else {
+#if GOOGLE_CUDA || TENSORFLOW_USE_ROCM
+    DPErrcheck(DPSetDevice(gpu_id));
+#endif  // GOOGLE_CUDA || TENSORFLOW_USE_ROCM
+    std::cout << "load model from: " << model << " to gpu " << gpu_id
+              << std::endl;
+  }
+  std::unordered_map<std::string, std::string> metadata = {{"type", ""}};
+  module = torch::jit::load(model, device, metadata);
+  module.eval();
+
+  get_env_nthreads(num_intra_nthreads, num_inter_nthreads);
+  if (num_inter_nthreads) {
+    try {
+      at::set_num_interop_threads(num_inter_nthreads);
+    } catch (...) {
+    }
+  }
+  if (num_intra_nthreads) {
+    try {
+      at::set_num_threads(num_intra_nthreads);
+    } catch (...) {
+    }
+  }
+
+  // Get model properties using run_method for C++ interface
+  auto rcut_result = module.run_method("get_rcut");
+  rcut = rcut_result.toDouble();
+
+  auto ntypes_result = module.run_method("get_ntypes");
+  ntypes = ntypes_result.toInt();
+
+  // Get task dimension from model method
+  auto task_dim_result = module.run_method("get_task_dim");
+  odim = task_dim_result.toInt();
+
+  // Get type map and set up sel_type
+  auto type_map_result = module.run_method("get_type_map");
+  auto type_map_list = type_map_result.toList();
+  sel_type.clear();
+
+  // For PyTorch models, all types are included (the backend handles exclusions
+  // internally) The model always outputs all types, but some results may be
+  // zero
+  for (size_t i = 0; i < type_map_list.size(); ++i) {
+    sel_type.push_back(i);
+  }
+  inited = true;
+}
+
+DeepTensorPT::~DeepTensorPT() {}
+
+void DeepTensorPT::get_type_map(std::string& type_map) {
+  auto type_map_result = module.run_method("get_type_map");
+  auto type_map_list = type_map_result.toList();
+  type_map.clear();
+  for (const torch::IValue& element : type_map_list) {
+    if (!type_map.empty()) {
+      type_map += " ";
+    }
+    type_map += torch::str(element);
+  }
+}
+
+template <typename VALUETYPE>
+void DeepTensorPT::compute(std::vector<VALUETYPE>& global_tensor,
+                           std::vector<VALUETYPE>& force,
+                           std::vector<VALUETYPE>& virial,
+                           std::vector<VALUETYPE>& atom_tensor,
+                           std::vector<VALUETYPE>& atom_virial,
+                           const std::vector<VALUETYPE>& coord,
+                           const std::vector<int>& atype,
+                           const std::vector<VALUETYPE>& box,
+                           const bool request_deriv) {
+  torch::Device device(torch::kCUDA, gpu_id);
+  if (!gpu_enabled) {
+    device = torch::Device(torch::kCPU);
+  }
+
+  int natoms = atype.size();
+  auto options = torch::TensorOptions().dtype(torch::kFloat64);
+  torch::ScalarType floatType = torch::kFloat64;
+  if (std::is_same<VALUETYPE, float>::value) {
+    options = torch::TensorOptions().dtype(torch::kFloat32);
+    floatType = torch::kFloat32;
+  }
+  auto int_options = torch::TensorOptions().dtype(torch::kInt64);
+
+  // Convert inputs to tensors
+  std::vector<VALUETYPE> coord_wrapped = coord;
+  at::Tensor coord_tensor =
+      torch::from_blob(coord_wrapped.data(), {1, natoms, 3}, options)
+          .to(device);
+
+  std::vector<std::int64_t> atype_64(atype.begin(), atype.end());
+  at::Tensor atype_tensor =
+      torch::from_blob(atype_64.data(), {1, natoms}, int_options).to(device);
+
+  c10::optional<torch::Tensor> box_tensor;
+  if (!box.empty()) {
+    box_tensor =
+        torch::from_blob(const_cast<VALUETYPE*>(box.data()), {1, 9}, options)
+            .to(device);
+  }
+
+  // Create input vector
+  std::vector<torch::jit::IValue> inputs;
+  inputs.push_back(coord_tensor);
+  inputs.push_back(atype_tensor);
+  inputs.push_back(box_tensor);
+
+  // Add None for fparam and aparam (not used by tensor models)
+  inputs.push_back(torch::jit::IValue());  // fparam = None
+  inputs.push_back(torch::jit::IValue());  // aparam = None
+  inputs.push_back(request_deriv);         // do_atomic_virial
+
+  // Forward pass through model
+  c10::Dict<c10::IValue, c10::IValue> outputs =
+      module.forward(inputs).toGenericDict();
+
+  // Extract global dipole/polar results
+  c10::IValue global_out;
+  if (outputs.contains("global_dipole")) {
+    global_out = outputs.at("global_dipole");
+  } else if (outputs.contains("global_polar")) {
+    global_out = outputs.at("global_polar");
+  } else {
+    throw deepmd::deepmd_exception(
+        "Cannot find global tensor output in model results");
+  }
+  torch::Tensor flat_global_ = global_out.toTensor().view({-1}).to(floatType);
+  torch::Tensor cpu_global_ = flat_global_.to(torch::kCPU);
+  global_tensor.assign(cpu_global_.data_ptr<VALUETYPE>(),
+                       cpu_global_.data_ptr<VALUETYPE>() + cpu_global_.numel());
+
+  // Extract atomic dipole/polar results
+  c10::IValue atom_out;
+  if (outputs.contains("dipole")) {
+    atom_out = outputs.at("dipole");
+  } else if (outputs.contains("polar")) {
+    atom_out = outputs.at("polar");
+  } else {
+    throw deepmd::deepmd_exception(
+        "Cannot find atomic tensor output in model results");
+  }
+  torch::Tensor flat_atom_ = atom_out.toTensor().view({-1}).to(floatType);
+  torch::Tensor cpu_atom_ = flat_atom_.to(torch::kCPU);
+  atom_tensor.assign(cpu_atom_.data_ptr<VALUETYPE>(),
+                     cpu_atom_.data_ptr<VALUETYPE>() + cpu_atom_.numel());
+
+  // Extract force results
+  c10::IValue force_ = outputs.at("force");
+  torch::Tensor flat_force_ = force_.toTensor().view({-1}).to(floatType);
+  torch::Tensor cpu_force_ = flat_force_.to(torch::kCPU);
+  force.assign(cpu_force_.data_ptr<VALUETYPE>(),
+               cpu_force_.data_ptr<VALUETYPE>() + cpu_force_.numel());
+
+  // Extract virial results
+  c10::IValue virial_ = outputs.at("virial");
+  torch::Tensor flat_virial_ = virial_.toTensor().view({-1}).to(floatType);
+  torch::Tensor cpu_virial_ = flat_virial_.to(torch::kCPU);
+  virial.assign(cpu_virial_.data_ptr<VALUETYPE>(),
+                cpu_virial_.data_ptr<VALUETYPE>() + cpu_virial_.numel());
+  // Extract atomic virial results if requested
+  if (request_deriv) {
+    c10::IValue atom_virial_ = outputs.at("atom_virial");
+    torch::Tensor flat_atom_virial_ =
+        atom_virial_.toTensor().view({-1}).to(floatType);
+    torch::Tensor cpu_atom_virial_ = flat_atom_virial_.to(torch::kCPU);
+    atom_virial.assign(
+        cpu_atom_virial_.data_ptr<VALUETYPE>(),
+        cpu_atom_virial_.data_ptr<VALUETYPE>() + cpu_atom_virial_.numel());
+  } else {
+    atom_virial.clear();
+  }
+}
+
+template <typename VALUETYPE>
+void DeepTensorPT::compute(std::vector<VALUETYPE>& global_tensor,
+                           std::vector<VALUETYPE>& force,
+                           std::vector<VALUETYPE>& virial,
+                           std::vector<VALUETYPE>& atom_tensor,
+                           std::vector<VALUETYPE>& atom_virial,
+                           const std::vector<VALUETYPE>& coord,
+                           const std::vector<int>& atype,
+                           const std::vector<VALUETYPE>& box,
+                           const int nghost,
+                           const InputNlist& lmp_list,
+                           const bool request_deriv) {
+  torch::Device device(torch::kCUDA, gpu_id);
+  if (!gpu_enabled) {
+    device = torch::Device(torch::kCPU);
+  }
+
+  int natoms = atype.size();
+  auto options = torch::TensorOptions().dtype(torch::kFloat64);
+  torch::ScalarType floatType = torch::kFloat64;
+  if (std::is_same<VALUETYPE, float>::value) {
+    options = torch::TensorOptions().dtype(torch::kFloat32);
+    floatType = torch::kFloat32;
+  }
+  auto int32_option =
+      torch::TensorOptions().device(torch::kCPU).dtype(torch::kInt32);
+  auto int_option =
+      torch::TensorOptions().device(torch::kCPU).dtype(torch::kInt64);
+
+  // Select real atoms following DeepPotPT pattern
+  std::vector<VALUETYPE> dcoord, aparam_;
+  std::vector<int> datype, fwd_map, bkw_map;
+  int nghost_real, nall_real, nloc_real;
+  int nall = natoms;
+  int nframes = 1;
+  std::vector<VALUETYPE> aparam;  // Empty for tensor models
+  select_real_atoms_coord(dcoord, datype, aparam_, nghost_real, fwd_map,
+                          bkw_map, nall_real, nloc_real, coord, atype, aparam,
+                          nghost, ntypes, nframes, 0, nall, false);
+
+  std::vector<VALUETYPE> coord_wrapped = dcoord;
+  at::Tensor coord_wrapped_Tensor =
+      torch::from_blob(coord_wrapped.data(), {1, nall_real, 3}, options)
+          .to(device);
+  std::vector<std::int64_t> atype_64(datype.begin(), datype.end());
+  at::Tensor atype_Tensor =
+      torch::from_blob(atype_64.data(), {1, nall_real}, int_option).to(device);
+
+  // Process neighbor list following DeepPotPT pattern
+  nlist_data.copy_from_nlist(lmp_list, nall - nghost);
+  nlist_data.shuffle_exclude_empty(fwd_map);
+  nlist_data.padding();
+
+  at::Tensor firstneigh = createNlistTensor(nlist_data.jlist);
+  firstneigh_tensor = firstneigh.to(torch::kInt64).to(device);
+
+  bool do_atom_virial_tensor = request_deriv;
+  c10::optional<torch::Tensor> fparam_tensor;
+  c10::optional<torch::Tensor> aparam_tensor;
+  c10::optional<torch::Tensor> mapping_tensor;
+
+  // Use forward_lower method following DeepPotPT pattern
+  c10::Dict<c10::IValue, c10::IValue> outputs =
+      module
+          .run_method("forward_lower", coord_wrapped_Tensor, atype_Tensor,
+                      firstneigh_tensor, mapping_tensor, fparam_tensor,
+                      aparam_tensor, do_atom_virial_tensor)
+          .toGenericDict();
+
+  // Extract outputs following DeepPotPT pattern
+  c10::IValue global_dipole_;
+  if (outputs.contains("global_dipole")) {
+    global_dipole_ = outputs.at("global_dipole");
+  } else if (outputs.contains("global_polar")) {
+    global_dipole_ = outputs.at("global_polar");
+  } else {
+    throw deepmd::deepmd_exception(
+        "Cannot find global tensor output in model results");
+  }
+  // in Python, here used double; however, in TF C++, float is used
+  // for consistency, we use float
+  torch::Tensor flat_global_ =
+      global_dipole_.toTensor().view({-1}).to(floatType);
+  torch::Tensor cpu_global_ = flat_global_.to(torch::kCPU);
+  global_tensor.assign(cpu_global_.data_ptr<VALUETYPE>(),
+                       cpu_global_.data_ptr<VALUETYPE>() + cpu_global_.numel());
+
+  c10::IValue force_ = outputs.at("extended_force");
+  torch::Tensor flat_force_ = force_.toTensor().view({-1}).to(floatType);
+  torch::Tensor cpu_force_ = flat_force_.to(torch::kCPU);
+  std::vector<VALUETYPE> dforce;
+  dforce.assign(cpu_force_.data_ptr<VALUETYPE>(),
+                cpu_force_.data_ptr<VALUETYPE>() + cpu_force_.numel());
+
+  c10::IValue virial_ = outputs.at("virial");
+  torch::Tensor flat_virial_ = virial_.toTensor().view({-1}).to(floatType);
+  torch::Tensor cpu_virial_ = flat_virial_.to(torch::kCPU);
+  virial.assign(cpu_virial_.data_ptr<VALUETYPE>(),
+                cpu_virial_.data_ptr<VALUETYPE>() + cpu_virial_.numel());
+
+  // bkw map for forces
+  force.resize(static_cast<size_t>(nframes) * odim * fwd_map.size() * 3);
+  for (int kk = 0; kk < odim; ++kk) {
+    select_map<VALUETYPE>(force.begin() + kk * fwd_map.size() * 3,
+                          dforce.begin() + kk * bkw_map.size() * 3, bkw_map, 3);
+  }
+
+  // Extract atomic dipoles/polars if available
+  c10::IValue atom_tensor_output;
+  int task_dim;
+  if (outputs.contains("dipole")) {
+    atom_tensor_output = outputs.at("dipole");
+    task_dim = 3;  // dipole has 3 components
+  } else if (outputs.contains("polar")) {
+    atom_tensor_output = outputs.at("polar");
+    task_dim = 9;  // polarizability has 9 components typically
+  } else {
+    throw deepmd::deepmd_exception(
+        "Cannot find atomic tensor output in model results");
+  }
+
+  torch::Tensor flat_atom_tensor_ =
+      atom_tensor_output.toTensor().view({-1}).to(floatType);
+  torch::Tensor cpu_atom_tensor_ = flat_atom_tensor_.to(torch::kCPU);
+  std::vector<VALUETYPE> datom_tensor;
+  datom_tensor.assign(
+      cpu_atom_tensor_.data_ptr<VALUETYPE>(),
+      cpu_atom_tensor_.data_ptr<VALUETYPE>() + cpu_atom_tensor_.numel());
+  atom_tensor.resize(static_cast<size_t>(nframes) * fwd_map.size() * task_dim);
+  select_map<VALUETYPE>(atom_tensor, datom_tensor, bkw_map, task_dim, nframes,
+                        fwd_map.size(), nall_real);
+
+  if (request_deriv) {
+    c10::IValue atom_virial_ = outputs.at("extended_virial");
+    torch::Tensor flat_atom_virial_ =
+        atom_virial_.toTensor().view({-1}).to(floatType);
+    torch::Tensor cpu_atom_virial_ = flat_atom_virial_.to(torch::kCPU);
+    std::vector<VALUETYPE> datom_virial;
+    datom_virial.assign(
+        cpu_atom_virial_.data_ptr<VALUETYPE>(),
+        cpu_atom_virial_.data_ptr<VALUETYPE>() + cpu_atom_virial_.numel());
+    atom_virial.resize(static_cast<size_t>(nframes) * odim * fwd_map.size() *
+                       9);
+    for (int kk = 0; kk < odim; ++kk) {
+      select_map<VALUETYPE>(atom_virial.begin() + kk * fwd_map.size() * 9,
+                            datom_virial.begin() + kk * bkw_map.size() * 9,
+                            bkw_map, 9);
+    }
+  }
+}
+
+// Public wrapper functions
+void DeepTensorPT::computew(std::vector<double>& global_tensor,
+                            std::vector<double>& force,
+                            std::vector<double>& virial,
+                            std::vector<double>& atom_tensor,
+                            std::vector<double>& atom_virial,
+                            const std::vector<double>& coord,
+                            const std::vector<int>& atype,
+                            const std::vector<double>& box,
+                            const bool request_deriv) {
+  translate_error([&] {
+    compute(global_tensor, force, virial, atom_tensor, atom_virial, coord,
+            atype, box, request_deriv);
+  });
+}
+
+void DeepTensorPT::computew(std::vector<float>& global_tensor,
+                            std::vector<float>& force,
+                            std::vector<float>& virial,
+                            std::vector<float>& atom_tensor,
+                            std::vector<float>& atom_virial,
+                            const std::vector<float>& coord,
+                            const std::vector<int>& atype,
+                            const std::vector<float>& box,
+                            const bool request_deriv) {
+  translate_error([&] {
+    compute(global_tensor, force, virial, atom_tensor, atom_virial, coord,
+            atype, box, request_deriv);
+  });
+}
+
+void DeepTensorPT::computew(std::vector<double>& global_tensor,
+                            std::vector<double>& force,
+                            std::vector<double>& virial,
+                            std::vector<double>& atom_tensor,
+                            std::vector<double>& atom_virial,
+                            const std::vector<double>& coord,
+                            const std::vector<int>& atype,
+                            const std::vector<double>& box,
+                            const int nghost,
+                            const InputNlist& inlist,
+                            const bool request_deriv) {
+  translate_error([&] {
+    compute(global_tensor, force, virial, atom_tensor, atom_virial, coord,
+            atype, box, nghost, inlist, request_deriv);
+  });
+}
+
+void DeepTensorPT::computew(std::vector<float>& global_tensor,
+                            std::vector<float>& force,
+                            std::vector<float>& virial,
+                            std::vector<float>& atom_tensor,
+                            std::vector<float>& atom_virial,
+                            const std::vector<float>& coord,
+                            const std::vector<int>& atype,
+                            const std::vector<float>& box,
+                            const int nghost,
+                            const InputNlist& inlist,
+                            const bool request_deriv) {
+  translate_error([&] {
+    compute(global_tensor, force, virial, atom_tensor, atom_virial, coord,
+            atype, box, nghost, inlist, request_deriv);
+  });
+}
+
+#endif  // BUILD_PYTORCH
diff --git a/source/api_cc/tests/test_deepdipole_pt.cc b/source/api_cc/tests/test_deepdipole_pt.cc
new file mode 100644
index 0000000000..70e46dd9e9
--- /dev/null
+++ b/source/api_cc/tests/test_deepdipole_pt.cc
@@ -0,0 +1,173 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+#include <fcntl.h>
+#include <gtest/gtest.h>
+#include <sys/stat.h>
+#include <sys/types.h>
+
+#include <algorithm>
+#include <cmath>
+#include <fstream>
+#include <sstream>
+#include <vector>
+
+#include "DeepTensor.h"
+#include "neighbor_list.h"
+#include "test_utils.h"
+
+template <class VALUETYPE>
+class TestInferDeepTensorPt : public ::testing::Test {
+ protected:
+  std::vector<VALUETYPE> coord = {12.83, 2.56, 2.18, 12.09, 2.87, 2.74,
+                                  00.25, 3.32, 1.68, 3.36,  3.00, 1.81,
+                                  3.51,  2.51, 2.60, 4.27,  3.22, 1.56};
+  std::vector<int> atype = {0, 1, 1, 0, 1, 1};
+  std::vector<VALUETYPE> box = {13., 0., 0., 0., 13., 0., 0., 0., 13.};
+
+  // Expected global tensor values from Python inference
+  std::vector<VALUETYPE> expected_global_tensor = {0.2338104, 0.23701073,
+                                                   0.2334505};
+
+  // Expected atomic tensor values from Python inference (flattened)
+  std::vector<VALUETYPE> expected_atom_tensor = {-0.1808925408386811,
+                                                 0.3190798607195795,
+                                                 0.04760079958216837,
+                                                 -0.0,
+                                                 -0.0,
+                                                 0.0,
+                                                 0.0,
+                                                 0.0,
+                                                 -0.0,
+                                                 0.4147029447879755,
+                                                 -0.08206913353381971,
+                                                 0.1858497008385067,
+                                                 0.0,
+                                                 -0.0,
+                                                 0.0,
+                                                 0.0,
+                                                 0.0,
+                                                 -0.0};
+
+  int natoms = 6;
+  int output_dim = 3;
+
+  deepmd::DeepTensor dt;
+
+  void SetUp() override {
+    std::string file_name = "../../tests/infer/deepdipole_pt.pth";
+    dt.init(file_name);
+  };
+
+  void TearDown() override {};
+};
+
+TYPED_TEST_SUITE(TestInferDeepTensorPt, ValueTypes);
+
+TYPED_TEST(TestInferDeepTensorPt, cpu_build_nlist) {
+  using VALUETYPE = TypeParam;
+  std::vector<VALUETYPE>& coord = this->coord;
+  std::vector<int>& atype = this->atype;
+  std::vector<VALUETYPE>& box = this->box;
+  std::vector<VALUETYPE>& expected_global_tensor = this->expected_global_tensor;
+  std::vector<VALUETYPE>& expected_atom_tensor = this->expected_atom_tensor;
+  int& natoms = this->natoms;
+  int& output_dim = this->output_dim;
+  deepmd::DeepTensor& dt = this->dt;
+  // Use reasonable tolerance for minimal trained model
+  double tensor_tol = 1e-6;
+
+  std::vector<VALUETYPE> global_tensor, force, virial, atom_tensor, atom_virial;
+
+  dt.compute(global_tensor, force, virial, atom_tensor, atom_virial, coord,
+             atype, box);
+
+  EXPECT_EQ(global_tensor.size(), output_dim);
+  EXPECT_EQ(atom_tensor.size(), natoms * output_dim);
+  EXPECT_EQ(force.size(), natoms * output_dim * 3);
+  EXPECT_EQ(virial.size(), output_dim * 9);
+  EXPECT_EQ(atom_virial.size(), natoms * output_dim * 9);
+
+  for (int ii = 0; ii < output_dim; ++ii) {
+    EXPECT_LT(fabs(global_tensor[ii] - expected_global_tensor[ii]), tensor_tol);
+  }
+
+  for (int ii = 0; ii < natoms * output_dim; ++ii) {
+    EXPECT_LT(fabs(atom_tensor[ii] - expected_atom_tensor[ii]), tensor_tol);
+  }
+}
+
+TYPED_TEST(TestInferDeepTensorPt, cpu_lmp_nlist) {
+  using VALUETYPE = TypeParam;
+  std::vector<VALUETYPE>& coord = this->coord;
+  std::vector<int>& atype = this->atype;
+  std::vector<VALUETYPE>& box = this->box;
+  std::vector<VALUETYPE>& expected_global_tensor = this->expected_global_tensor;
+  std::vector<VALUETYPE>& expected_atom_tensor = this->expected_atom_tensor;
+  int& natoms = this->natoms;
+  int& output_dim = this->output_dim;
+  deepmd::DeepTensor& dt = this->dt;
+  double ener_tol = 1e-6;
+
+  float rc = dt.cutoff();
+  int nloc = coord.size() / 3;
+  std::vector<VALUETYPE> coord_cpy;
+  std::vector<int> atype_cpy, mapping;
+  std::vector<std::vector<int> > nlist_data;
+  _build_nlist<VALUETYPE>(nlist_data, coord_cpy, atype_cpy, mapping, coord,
+                          atype, box, rc);
+  int nall = coord_cpy.size() / 3;
+  std::vector<int> ilist(nloc), numneigh(nloc);
+  std::vector<int*> firstneigh(nloc);
+  deepmd::InputNlist inlist(nloc, &ilist[0], &numneigh[0], &firstneigh[0]);
+  convert_nlist(inlist, nlist_data);
+
+  std::vector<VALUETYPE> global_tensor, force, virial, atom_tensor, atom_virial;
+
+  dt.compute(global_tensor, force, virial, atom_tensor, atom_virial, coord_cpy,
+             atype_cpy, box, nall - nloc, inlist);
+
+  EXPECT_EQ(global_tensor.size(), output_dim);
+  EXPECT_EQ(atom_tensor.size(), nall * output_dim);
+
+  for (int ii = 0; ii < output_dim; ++ii) {
+    EXPECT_LT(fabs(global_tensor[ii] - expected_global_tensor[ii]), ener_tol);
+  }
+
+  for (int ii = 0; ii < natoms * output_dim; ++ii) {
+    EXPECT_LT(fabs(atom_tensor[ii] - expected_atom_tensor[ii]), ener_tol);
+  }
+}
+
+TYPED_TEST(TestInferDeepTensorPt, print_summary) {
+  deepmd::DeepTensor& dt = this->dt;
+  dt.print_summary("");
+}
+
+TYPED_TEST(TestInferDeepTensorPt, get_type_map) {
+  deepmd::DeepTensor& dt = this->dt;
+  std::string type_map_str;
+  dt.get_type_map(type_map_str);
+  // Parse the type map string manually
+  std::vector<std::string> type_map;
+  std::istringstream iss(type_map_str);
+  std::string token;
+  while (iss >> token) {
+    type_map.push_back(token);
+  }
+  EXPECT_EQ(type_map.size(), 2);
+  EXPECT_EQ(type_map[0], "O");
+  EXPECT_EQ(type_map[1], "H");
+}
+
+TYPED_TEST(TestInferDeepTensorPt, get_properties) {
+  deepmd::DeepTensor& dt = this->dt;
+
+  EXPECT_EQ(dt.numb_types(), 2);
+  EXPECT_EQ(dt.output_dim(), 3);
+  EXPECT_DOUBLE_EQ(dt.cutoff(), 4.0);
+
+  std::vector<int> sel_types = dt.sel_types();
+  EXPECT_EQ(sel_types.size(), 2);  // PyTorch models always return all types
+  EXPECT_EQ(sel_types[0], 0);      // Type 0 (O)
+  EXPECT_EQ(sel_types[1],
+            1);  // Type 1 (H) - included but may have zero results
+}
diff --git a/source/tests/infer/deepdipole_pt.pth b/source/tests/infer/deepdipole_pt.pth
new file mode 100644
index 0000000000..4c93a1b864
Binary files /dev/null and b/source/tests/infer/deepdipole_pt.pth differ