From 19c6a688ffc3ce9e5d7ffb482f925dbd59a05880 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Yifan=20Li=E6=9D=8E=E4=B8=80=E5=B8=86?=
 <yifanl0716@gmail.com>
Date: Sun, 17 Sep 2023 19:50:33 -0500
Subject: [PATCH] lmp: support other units for both pair deepmd and fix dplr
 (#2800)
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

This PR is intended to do 2 things:
1. support other unit styles like si and nano (except lj);
2. support the non-metal units for fix dplr and compute deeptensor/atom.

Unittests for both features mentioned above have also been added in this
pull request.

---------

Signed-off-by: Yifan Li李一帆 <yifanl0716@gmail.com>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Jinzhe Zeng <jinzhe.zeng@rutgers.edu>
---
 doc/third-party/lammps-command.md      |  10 ++
 source/lmp/compute_deeptensor_atom.cpp |  31 ++---
 source/lmp/compute_deeptensor_atom.h   |   1 +
 source/lmp/fix_dplr.cpp                |  59 +++++----
 source/lmp/fix_dplr.h                  |   1 +
 source/lmp/pair_deepmd.cpp             | 140 ++++++++++++----------
 source/lmp/pair_deepmd.h               |   2 +-
 source/lmp/tests/constants.py          |  17 +++
 source/lmp/tests/test_deeptensor.py    |  59 ++++++++-
 source/lmp/tests/test_dplr.py          |  80 +++++++++++--
 source/lmp/tests/test_lammps.py        | 160 ++++++++++++++++---------
 source/lmp/tests/write_lmp_data.py     |  22 ++--
 12 files changed, 406 insertions(+), 176 deletions(-)
 create mode 100644 source/lmp/tests/constants.py

diff --git a/doc/third-party/lammps-command.md b/doc/third-party/lammps-command.md
index e1d482381f..a9d849bc7c 100644
--- a/doc/third-party/lammps-command.md
+++ b/doc/third-party/lammps-command.md
@@ -1,5 +1,14 @@
 # LAMMPS commands
 
+## units
+All units in LAMMPS except `lj` are supported. `lj` is not supported.
+
+The most commonly used units are `metal`, since the internal units of distance, energy, force, and charge in DeePMD-kit are `\AA`, `eV`, `eV / \AA`, and `proton charge`, respectively. These units are consistent with the `metal` units in LAMMPS.
+
+If one wants to use other units like `real` or `si`, it is welcome to do so. There is no need to do the unit conversion mannualy. The unit conversion is done automatically by LAMMPS.
+
+The only thing that one needs to take care is the unit of the output of `compute deeptensor/atom`. Working with `metal` units for `compute deeptensor/atom` is totally fine, since there is no unit conversion. For other unit styles, we currently assume that the output of the `compute deeptensor/atom` command has the unit of distance and have applied the unit conversion factor of distance. If a user wants to infer quantities with units other than distance, the user is encouraged to open a GitHub feature request, so that the unit conversion factor can be added.
+
 ## Enable DeePMD-kit plugin (plugin mode)
 
 If you are using the plugin mode, enable DeePMD-kit package in LAMMPS with `plugin` command:
@@ -119,6 +128,7 @@ dump            1 all custom 100 water.dump id type c_dipole[1] c_dipole[2] c_di
 
 ### Restrictions
 - The `deeptensor/atom` compute is provided in the USER-DEEPMD package, which is compiled from the DeePMD-kit, visit the [DeePMD-kit website](https://github.com/deepmodeling/deepmd-kit) for more information.
+- For the issue of using a unit style for `compute deeptensor/atom`, refer to the discussions in [units](#units) of this page.
 
 
 ## Long-range interaction
diff --git a/source/lmp/compute_deeptensor_atom.cpp b/source/lmp/compute_deeptensor_atom.cpp
index 0de523e1bf..2f4486002e 100644
--- a/source/lmp/compute_deeptensor_atom.cpp
+++ b/source/lmp/compute_deeptensor_atom.cpp
@@ -26,12 +26,12 @@ using namespace LAMMPS_NS;
 
 ComputeDeeptensorAtom::ComputeDeeptensorAtom(LAMMPS *lmp, int narg, char **arg)
     : Compute(lmp, narg, arg), dp(lmp), tensor(nullptr) {
-  if (!(strcmp(update->unit_style, "metal") == 0 ||
-        strcmp(update->unit_style, "real") == 0)) {
-    error->all(
-        FLERR,
-        "Compute deeptensor/atom requires metal or real unit; please set it by "
-        "\"units metal\" or \"units real\"");
+  if (strcmp(update->unit_style, "lj") == 0) {
+    error->all(FLERR,
+               "Compute deeptensor/atom does not support unit style lj. Please "
+               "use other "
+               "unit styles like metal or real unit instead. You may set it by "
+               "\"units metal\" or \"units real\"");
   }
 
   if (narg < 4) {
@@ -57,6 +57,8 @@ ComputeDeeptensorAtom::ComputeDeeptensorAtom(LAMMPS *lmp, int narg, char **arg)
   timeflag = 1;
 
   nmax = 0;
+
+  dist_unit_cvt_factor = force->angstrom;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -116,16 +118,17 @@ void ComputeDeeptensorAtom::compute_peratom() {
     dtype[ii] = type[ii] - 1;
   }
   // get box
-  dbox[0] = domain->h[0];  // xx
-  dbox[4] = domain->h[1];  // yy
-  dbox[8] = domain->h[2];  // zz
-  dbox[7] = domain->h[3];  // zy
-  dbox[6] = domain->h[4];  // zx
-  dbox[3] = domain->h[5];  // yx
+  dbox[0] = domain->h[0] / dist_unit_cvt_factor;  // xx
+  dbox[4] = domain->h[1] / dist_unit_cvt_factor;  // yy
+  dbox[8] = domain->h[2] / dist_unit_cvt_factor;  // zz
+  dbox[7] = domain->h[3] / dist_unit_cvt_factor;  // zy
+  dbox[6] = domain->h[4] / dist_unit_cvt_factor;  // zx
+  dbox[3] = domain->h[5] / dist_unit_cvt_factor;  // yx
   // get coord
   for (int ii = 0; ii < nall; ++ii) {
     for (int dd = 0; dd < 3; ++dd) {
-      dcoord[ii * 3 + dd] = x[ii][dd] - domain->boxlo[dd];
+      dcoord[ii * 3 + dd] =
+          (x[ii][dd] - domain->boxlo[dd]) / dist_unit_cvt_factor;
     }
   }
 
@@ -155,7 +158,7 @@ void ComputeDeeptensorAtom::compute_peratom() {
     // record when selected and in group
     if (selected && ingroup) {
       for (int jj = 0; jj < size_peratom_cols; ++jj) {
-        tensor[ii][jj] = atensor[iter_tensor + jj];
+        tensor[ii][jj] = atensor[iter_tensor + jj] * dist_unit_cvt_factor;
       }
     }
     // if not selected or not in group set to 0.
diff --git a/source/lmp/compute_deeptensor_atom.h b/source/lmp/compute_deeptensor_atom.h
index ab25da4246..a90283aa9e 100644
--- a/source/lmp/compute_deeptensor_atom.h
+++ b/source/lmp/compute_deeptensor_atom.h
@@ -36,6 +36,7 @@ class ComputeDeeptensorAtom : public Compute {
   void compute_peratom() override;
   double memory_usage() override;
   void init_list(int, class NeighList *) override;
+  double dist_unit_cvt_factor;
 
  private:
   int nmax;
diff --git a/source/lmp/fix_dplr.cpp b/source/lmp/fix_dplr.cpp
index a0df7efd24..77bf0d56c0 100644
--- a/source/lmp/fix_dplr.cpp
+++ b/source/lmp/fix_dplr.cpp
@@ -61,10 +61,11 @@ FixDPLR::FixDPLR(LAMMPS *lmp, int narg, char **arg)
   qe2f = force->qe2f;
   xstyle = ystyle = zstyle = NONE;
 
-  if (strcmp(update->unit_style, "metal") != 0) {
-    error->all(
-        FLERR,
-        "Fix dplr requires metal unit, please set it by \"units metal\"");
+  if (strcmp(update->unit_style, "lj") == 0) {
+    error->all(FLERR,
+               "Fix dplr does not support unit style lj. Please use other "
+               "unit styles like metal or real unit instead. You may set it by "
+               "\"units metal\" or \"units real\"");
   }
 
   int iarg = 3;
@@ -142,6 +143,9 @@ FixDPLR::FixDPLR(LAMMPS *lmp, int narg, char **arg)
   if (!pair_deepmd) {
     error->all(FLERR, "pair_style deepmd should be set before this fix\n");
   }
+  ener_unit_cvt_factor = pair_deepmd->ener_unit_cvt_factor;
+  dist_unit_cvt_factor = pair_deepmd->dist_unit_cvt_factor;
+  force_unit_cvt_factor = pair_deepmd->force_unit_cvt_factor;
 
   int n = atom->ntypes;
   std::vector<std::string> type_names = pair_deepmd->type_names;
@@ -445,16 +449,17 @@ void FixDPLR::pre_force(int vflag) {
     dtype[ii] = type_idx_map[type[ii] - 1];
   }
   // get box
-  dbox[0] = domain->h[0];  // xx
-  dbox[4] = domain->h[1];  // yy
-  dbox[8] = domain->h[2];  // zz
-  dbox[7] = domain->h[3];  // zy
-  dbox[6] = domain->h[4];  // zx
-  dbox[3] = domain->h[5];  // yx
+  dbox[0] = domain->h[0] / dist_unit_cvt_factor;  // xx
+  dbox[4] = domain->h[1] / dist_unit_cvt_factor;  // yy
+  dbox[8] = domain->h[2] / dist_unit_cvt_factor;  // zz
+  dbox[7] = domain->h[3] / dist_unit_cvt_factor;  // zy
+  dbox[6] = domain->h[4] / dist_unit_cvt_factor;  // zx
+  dbox[3] = domain->h[5] / dist_unit_cvt_factor;  // yx
   // get coord
   for (int ii = 0; ii < nall; ++ii) {
     for (int dd = 0; dd < 3; ++dd) {
-      dcoord[ii * 3 + dd] = x[ii][dd] - domain->boxlo[dd];
+      dcoord[ii * 3 + dd] =
+          (x[ii][dd] - domain->boxlo[dd]) / dist_unit_cvt_factor;
     }
   }
   // get lammps nlist
@@ -523,9 +528,11 @@ void FixDPLR::pre_force(int vflag) {
     int res_idx = sel_fwd[idx0];
     // int ret_idx = dpl_bwd[res_idx];
     for (int dd = 0; dd < 3; ++dd) {
-      x[idx1][dd] = x[idx0][dd] + tensor[res_idx * 3 + dd];
+      x[idx1][dd] =
+          x[idx0][dd] + tensor[res_idx * 3 + dd] * dist_unit_cvt_factor;
       // res_buff[idx1 * odim + dd] = tensor[res_idx * odim + dd];
-      dipole_recd[idx0 * 3 + dd] = tensor[res_idx * 3 + dd];
+      dipole_recd[idx0 * 3 + dd] =
+          tensor[res_idx * 3 + dd] * dist_unit_cvt_factor;
     }
   }
   // cout << "-------------------- fix/dplr: pre force " << endl;
@@ -568,17 +575,18 @@ void FixDPLR::post_force(int vflag) {
     for (int ii = 0; ii < nall; ++ii) {
       dtype[ii] = type_idx_map[type[ii] - 1];
     }
-    dbox[0] = domain->h[0];  // xx
-    dbox[4] = domain->h[1];  // yy
-    dbox[8] = domain->h[2];  // zz
-    dbox[7] = domain->h[3];  // zy
-    dbox[6] = domain->h[4];  // zx
-    dbox[3] = domain->h[5];  // yx
+    dbox[0] = domain->h[0] / dist_unit_cvt_factor;  // xx
+    dbox[4] = domain->h[1] / dist_unit_cvt_factor;  // yy
+    dbox[8] = domain->h[2] / dist_unit_cvt_factor;  // zz
+    dbox[7] = domain->h[3] / dist_unit_cvt_factor;  // zy
+    dbox[6] = domain->h[4] / dist_unit_cvt_factor;  // zx
+    dbox[3] = domain->h[5] / dist_unit_cvt_factor;  // yx
     // get coord
     double **x = atom->x;
     for (int ii = 0; ii < nall; ++ii) {
       for (int dd = 0; dd < 3; ++dd) {
-        dcoord[ii * 3 + dd] = x[ii][dd] - domain->boxlo[dd];
+        dcoord[ii * 3 + dd] =
+            (x[ii][dd] - domain->boxlo[dd]) / dist_unit_cvt_factor;
       }
     }
     // revise force according to efield
@@ -599,7 +607,7 @@ void FixDPLR::post_force(int vflag) {
     for (int ii = 0; ii < nlocal; ++ii) {
       double tmpf[3];
       for (int dd = 0; dd < 3; ++dd) {
-        tmpf[dd] = q[ii] * efield[dd];
+        tmpf[dd] = q[ii] * efield[dd] * force->qe2f;
       }
       for (int dd = 0; dd < 3; ++dd) {
         dfele[ii * 3 + dd] += tmpf[dd];
@@ -632,8 +640,17 @@ void FixDPLR::post_force(int vflag) {
   vector<FLOAT_PREC> dfcorr, dvcorr;
   // compute
   try {
+    for (int ii = 0; ii < nlocal * 3; ++ii) {
+      dfele[ii] /= force_unit_cvt_factor;
+    }
     dtm.compute(dfcorr, dvcorr, dcoord, dtype, dbox, valid_pairs, dfele, nghost,
                 lmp_list);
+    for (int ii = 0; ii < nlocal * 3; ++ii) {
+      dfcorr[ii] *= force_unit_cvt_factor;
+    }
+    for (int ii = 0; ii < 9; ++ii) {
+      dvcorr[ii] *= ener_unit_cvt_factor;
+    }
   } catch (deepmd_compat::deepmd_exception &e) {
     error->one(FLERR, e.what());
   }
diff --git a/source/lmp/fix_dplr.h b/source/lmp/fix_dplr.h
index 23ae1c818d..a6822fe4fe 100644
--- a/source/lmp/fix_dplr.h
+++ b/source/lmp/fix_dplr.h
@@ -54,6 +54,7 @@ class FixDPLR : public Fix {
   void unpack_reverse_comm(int, int *, double *) override;
   double compute_scalar(void) override;
   double compute_vector(int) override;
+  double ener_unit_cvt_factor, dist_unit_cvt_factor, force_unit_cvt_factor;
 
  private:
   PairDeepMD *pair_deepmd;
diff --git a/source/lmp/pair_deepmd.cpp b/source/lmp/pair_deepmd.cpp
index 489c31ff19..3fa592bf58 100644
--- a/source/lmp/pair_deepmd.cpp
+++ b/source/lmp/pair_deepmd.cpp
@@ -344,18 +344,16 @@ PairDeepMD::PairDeepMD(LAMMPS *lmp)
   if (lmp->citeme) {
     lmp->citeme->add(cite_user_deepmd_package);
   }
-  int unit_convert;
-  if (strcmp(update->unit_style, "metal") == 0) {
-    unit_convert = utils::NOCONVERT;
-  } else if (strcmp(update->unit_style, "real") == 0) {
-    unit_convert = utils::METAL2REAL;
-  } else {
+  if (strcmp(update->unit_style, "lj") == 0) {
     error->all(FLERR,
-               "Pair deepmd requires metal or real unit, please set it by "
+               "Pair deepmd does not support unit style lj. Please use other "
+               "unit styles like metal or real unit instead. You may set it by "
                "\"units metal\" or \"units real\"");
   }
-  ener_unit_cvt_factor =
-      utils::get_conversion_factor(utils::ENERGY, unit_convert);
+  ener_unit_cvt_factor = force->boltz / 8.617343e-5;
+  dist_unit_cvt_factor = force->angstrom;
+  force_unit_cvt_factor = ener_unit_cvt_factor / dist_unit_cvt_factor;
+
   restartinfo = 1;
 #if LAMMPS_VERSION_NUMBER >= 20201130
   centroidstressflag =
@@ -368,7 +366,6 @@ PairDeepMD::PairDeepMD(LAMMPS *lmp)
   pppmflag = 1;
   respa_enable = 0;
   writedata = 0;
-  unit_convert_flag = utils::get_supported_conversions(utils::ENERGY);
 
   cutoff = 0.;
   numb_types = 0;
@@ -480,17 +477,18 @@ void PairDeepMD::compute(int eflag, int vflag) {
   vector<double> daparam;
 
   // get box
-  dbox[0] = domain->h[0];  // xx
-  dbox[4] = domain->h[1];  // yy
-  dbox[8] = domain->h[2];  // zz
-  dbox[7] = domain->h[3];  // zy
-  dbox[6] = domain->h[4];  // zx
-  dbox[3] = domain->h[5];  // yx
+  dbox[0] = domain->h[0] / dist_unit_cvt_factor;  // xx
+  dbox[4] = domain->h[1] / dist_unit_cvt_factor;  // yy
+  dbox[8] = domain->h[2] / dist_unit_cvt_factor;  // zz
+  dbox[7] = domain->h[3] / dist_unit_cvt_factor;  // zy
+  dbox[6] = domain->h[4] / dist_unit_cvt_factor;  // zx
+  dbox[3] = domain->h[5] / dist_unit_cvt_factor;  // yx
 
   // get coord
   for (int ii = 0; ii < nall; ++ii) {
     for (int dd = 0; dd < 3; ++dd) {
-      dcoord[ii * 3 + dd] = x[ii][dd] - domain->boxlo[dd];
+      dcoord[ii * 3 + dd] =
+          (x[ii][dd] - domain->boxlo[dd]) / dist_unit_cvt_factor;
     }
   }
 
@@ -586,7 +584,7 @@ void PairDeepMD::compute(int eflag, int vflag) {
         }
         if (eflag_atom) {
           for (int ii = 0; ii < nlocal; ++ii) {
-            eatom[ii] += scale[1][1] * deatom[ii];
+            eatom[ii] += scale[1][1] * deatom[ii] * ener_unit_cvt_factor;
           }
         }
         // Added by Davide Tisi 2020
@@ -600,15 +598,24 @@ void PairDeepMD::compute(int eflag, int vflag) {
             // vatom[ii][3] += 1.0 * dvatom[9*ii+3];
             // vatom[ii][4] += 1.0 * dvatom[9*ii+6];
             // vatom[ii][5] += 1.0 * dvatom[9*ii+7];
-            cvatom[ii][0] += scale[1][1] * dvatom[9 * ii + 0];  // xx
-            cvatom[ii][1] += scale[1][1] * dvatom[9 * ii + 4];  // yy
-            cvatom[ii][2] += scale[1][1] * dvatom[9 * ii + 8];  // zz
-            cvatom[ii][3] += scale[1][1] * dvatom[9 * ii + 3];  // xy
-            cvatom[ii][4] += scale[1][1] * dvatom[9 * ii + 6];  // xz
-            cvatom[ii][5] += scale[1][1] * dvatom[9 * ii + 7];  // yz
-            cvatom[ii][6] += scale[1][1] * dvatom[9 * ii + 1];  // yx
-            cvatom[ii][7] += scale[1][1] * dvatom[9 * ii + 2];  // zx
-            cvatom[ii][8] += scale[1][1] * dvatom[9 * ii + 5];  // zy
+            cvatom[ii][0] +=
+                scale[1][1] * dvatom[9 * ii + 0] * ener_unit_cvt_factor;  // xx
+            cvatom[ii][1] +=
+                scale[1][1] * dvatom[9 * ii + 4] * ener_unit_cvt_factor;  // yy
+            cvatom[ii][2] +=
+                scale[1][1] * dvatom[9 * ii + 8] * ener_unit_cvt_factor;  // zz
+            cvatom[ii][3] +=
+                scale[1][1] * dvatom[9 * ii + 3] * ener_unit_cvt_factor;  // xy
+            cvatom[ii][4] +=
+                scale[1][1] * dvatom[9 * ii + 6] * ener_unit_cvt_factor;  // xz
+            cvatom[ii][5] +=
+                scale[1][1] * dvatom[9 * ii + 7] * ener_unit_cvt_factor;  // yz
+            cvatom[ii][6] +=
+                scale[1][1] * dvatom[9 * ii + 1] * ener_unit_cvt_factor;  // yx
+            cvatom[ii][7] +=
+                scale[1][1] * dvatom[9 * ii + 2] * ener_unit_cvt_factor;  // zx
+            cvatom[ii][8] +=
+                scale[1][1] * dvatom[9 * ii + 5] * ener_unit_cvt_factor;  // zy
           }
         }
       }
@@ -638,7 +645,7 @@ void PairDeepMD::compute(int eflag, int vflag) {
       dvatom = all_atom_virial[0];
       if (eflag_atom) {
         for (int ii = 0; ii < nlocal; ++ii) {
-          eatom[ii] += scale[1][1] * deatom[ii];
+          eatom[ii] += scale[1][1] * deatom[ii] * ener_unit_cvt_factor;
         }
       }
       // Added by Davide Tisi 2020
@@ -652,15 +659,24 @@ void PairDeepMD::compute(int eflag, int vflag) {
           // vatom[ii][3] += 1.0 * dvatom[9*ii+3];
           // vatom[ii][4] += 1.0 * dvatom[9*ii+6];
           // vatom[ii][5] += 1.0 * dvatom[9*ii+7];
-          cvatom[ii][0] += scale[1][1] * dvatom[9 * ii + 0];  // xx
-          cvatom[ii][1] += scale[1][1] * dvatom[9 * ii + 4];  // yy
-          cvatom[ii][2] += scale[1][1] * dvatom[9 * ii + 8];  // zz
-          cvatom[ii][3] += scale[1][1] * dvatom[9 * ii + 3];  // xy
-          cvatom[ii][4] += scale[1][1] * dvatom[9 * ii + 6];  // xz
-          cvatom[ii][5] += scale[1][1] * dvatom[9 * ii + 7];  // yz
-          cvatom[ii][6] += scale[1][1] * dvatom[9 * ii + 1];  // yx
-          cvatom[ii][7] += scale[1][1] * dvatom[9 * ii + 2];  // zx
-          cvatom[ii][8] += scale[1][1] * dvatom[9 * ii + 5];  // zy
+          cvatom[ii][0] +=
+              scale[1][1] * dvatom[9 * ii + 0] * ener_unit_cvt_factor;  // xx
+          cvatom[ii][1] +=
+              scale[1][1] * dvatom[9 * ii + 4] * ener_unit_cvt_factor;  // yy
+          cvatom[ii][2] +=
+              scale[1][1] * dvatom[9 * ii + 8] * ener_unit_cvt_factor;  // zz
+          cvatom[ii][3] +=
+              scale[1][1] * dvatom[9 * ii + 3] * ener_unit_cvt_factor;  // xy
+          cvatom[ii][4] +=
+              scale[1][1] * dvatom[9 * ii + 6] * ener_unit_cvt_factor;  // xz
+          cvatom[ii][5] +=
+              scale[1][1] * dvatom[9 * ii + 7] * ener_unit_cvt_factor;  // yz
+          cvatom[ii][6] +=
+              scale[1][1] * dvatom[9 * ii + 1] * ener_unit_cvt_factor;  // yx
+          cvatom[ii][7] +=
+              scale[1][1] * dvatom[9 * ii + 2] * ener_unit_cvt_factor;  // zx
+          cvatom[ii][8] +=
+              scale[1][1] * dvatom[9 * ii + 5] * ener_unit_cvt_factor;  // zy
         }
       }
       if (out_freq > 0 && update->ntimestep % out_freq == 0) {
@@ -729,12 +745,12 @@ void PairDeepMD::compute(int eflag, int vflag) {
           all_v_avg = sqrt(all_v_avg / 9);
         }
         if (rank == 0) {
-          all_v_max *= scale[1][1];
-          all_v_min *= scale[1][1];
-          all_v_avg *= scale[1][1];
-          all_f_max *= scale[1][1];
-          all_f_min *= scale[1][1];
-          all_f_avg *= scale[1][1];
+          all_v_max *= scale[1][1] * ener_unit_cvt_factor;
+          all_v_min *= scale[1][1] * ener_unit_cvt_factor;
+          all_v_avg *= scale[1][1] * ener_unit_cvt_factor;
+          all_f_max *= scale[1][1] * force_unit_cvt_factor;
+          all_f_min *= scale[1][1] * force_unit_cvt_factor;
+          all_f_avg *= scale[1][1] * force_unit_cvt_factor;
           fp << setw(12) << update->ntimestep << " " << setw(18) << all_v_max
              << " " << setw(18) << all_v_min << " " << setw(18) << all_v_avg
              << " " << setw(18) << all_f_max << " " << setw(18) << all_f_min
@@ -760,7 +776,8 @@ void PairDeepMD::compute(int eflag, int vflag) {
                       displacements, MPI_DOUBLE, 0, world);
           if (rank == 0) {
             for (int dd = 0; dd < all_nlocal; ++dd) {
-              std_f_all[tagrecv[dd] - 1] = stdfrecv[dd] * scale[1][1];
+              std_f_all[tagrecv[dd] - 1] =
+                  stdfrecv[dd] * scale[1][1] * force_unit_cvt_factor;
             }
             for (int dd = 0; dd < all_nlocal; ++dd) {
               fp << " " << setw(18) << std_f_all[dd];
@@ -790,7 +807,7 @@ void PairDeepMD::compute(int eflag, int vflag) {
   if (!atom->sp_flag) {
     for (int ii = 0; ii < nall; ++ii) {
       for (int dd = 0; dd < 3; ++dd) {
-        f[ii][dd] += scale[1][1] * dforce[3 * ii + dd];
+        f[ii][dd] += scale[1][1] * dforce[3 * ii + dd] * force_unit_cvt_factor;
       }
     }
   } else {
@@ -799,13 +816,14 @@ void PairDeepMD::compute(int eflag, int vflag) {
     for (int ii = 0; ii < nall; ++ii) {
       for (int dd = 0; dd < 3; ++dd) {
         int new_idx = new_idx_map[ii];
-        f[ii][dd] += scale[1][1] * dforce[3 * new_idx + dd];
+        f[ii][dd] +=
+            scale[1][1] * dforce[3 * new_idx + dd] * force_unit_cvt_factor;
         if (dtype[ii] < numb_types_spin && ii < nlocal) {
           fm[ii][dd] += scale[1][1] * dforce[3 * (new_idx + nlocal) + dd] /
-                        (hbar / spin_norm[dtype[ii]]);
+                        (hbar / spin_norm[dtype[ii]]) * force_unit_cvt_factor;
         } else if (dtype[ii] < numb_types_spin) {
           fm[ii][dd] += scale[1][1] * dforce[3 * (new_idx + nghost) + dd] /
-                        (hbar / spin_norm[dtype[ii]]);
+                        (hbar / spin_norm[dtype[ii]]) * force_unit_cvt_factor;
         }
       }
     }
@@ -819,15 +837,15 @@ void PairDeepMD::compute(int eflag, int vflag) {
 
   // accumulate energy and virial
   if (eflag) {
-    eng_vdwl += scale[1][1] * dener;
+    eng_vdwl += scale[1][1] * dener * ener_unit_cvt_factor;
   }
   if (vflag) {
-    virial[0] += 1.0 * dvirial[0] * scale[1][1];
-    virial[1] += 1.0 * dvirial[4] * scale[1][1];
-    virial[2] += 1.0 * dvirial[8] * scale[1][1];
-    virial[3] += 1.0 * dvirial[3] * scale[1][1];
-    virial[4] += 1.0 * dvirial[6] * scale[1][1];
-    virial[5] += 1.0 * dvirial[7] * scale[1][1];
+    virial[0] += 1.0 * dvirial[0] * scale[1][1] * ener_unit_cvt_factor;
+    virial[1] += 1.0 * dvirial[4] * scale[1][1] * ener_unit_cvt_factor;
+    virial[2] += 1.0 * dvirial[8] * scale[1][1] * ener_unit_cvt_factor;
+    virial[3] += 1.0 * dvirial[3] * scale[1][1] * ener_unit_cvt_factor;
+    virial[4] += 1.0 * dvirial[6] * scale[1][1] * ener_unit_cvt_factor;
+    virial[5] += 1.0 * dvirial[7] * scale[1][1] * ener_unit_cvt_factor;
   }
 }
 
@@ -854,7 +872,7 @@ void PairDeepMD::allocate() {
         continue;
       }
       setflag[i][j] = 1;
-      scale[i][j] = 1.0 * ener_unit_cvt_factor;
+      scale[i][j] = 1.0;
     }
   }
 }
@@ -904,7 +922,7 @@ void PairDeepMD::settings(int narg, char **arg) {
     } catch (deepmd_compat::deepmd_exception &e) {
       error->one(FLERR, e.what());
     }
-    cutoff = deep_pot.cutoff();
+    cutoff = deep_pot.cutoff() * dist_unit_cvt_factor;
     numb_types = deep_pot.numb_types();
     numb_types_spin = deep_pot.numb_types_spin();
     dim_fparam = deep_pot.dim_fparam();
@@ -917,12 +935,12 @@ void PairDeepMD::settings(int narg, char **arg) {
     } catch (deepmd_compat::deepmd_exception &e) {
       error->one(FLERR, e.what());
     }
-    cutoff = deep_pot_model_devi.cutoff();
+    cutoff = deep_pot_model_devi.cutoff() * dist_unit_cvt_factor;
     numb_types = deep_pot_model_devi.numb_types();
     numb_types_spin = deep_pot_model_devi.numb_types_spin();
     dim_fparam = deep_pot_model_devi.dim_fparam();
     dim_aparam = deep_pot_model_devi.dim_aparam();
-    assert(cutoff == deep_pot.cutoff());
+    assert(cutoff == deep_pot.cutoff() * dist_unit_cvt_factor);
     assert(numb_types == deep_pot.numb_types());
     assert(numb_types_spin == deep_pot.numb_types_spin());
     assert(dim_fparam == deep_pot.dim_fparam());
@@ -1197,7 +1215,7 @@ void PairDeepMD::coeff(int narg, char **arg) {
   for (int i = ilo; i <= ihi; i++) {
     for (int j = MAX(jlo, i); j <= jhi; j++) {
       setflag[i][j] = 1;
-      scale[i][j] = 1.0 * ener_unit_cvt_factor;
+      scale[i][j] = 1.0;
       if (i > numb_types || j > numb_types) {
         char warning_msg[1024];
         sprintf(warning_msg,
@@ -1244,7 +1262,7 @@ double PairDeepMD::init_one(int i, int j) {
   }
 
   if (setflag[i][j] == 0) {
-    scale[i][j] = 1.0 * ener_unit_cvt_factor;
+    scale[i][j] = 1.0;
   }
   scale[j][i] = scale[i][j];
 
diff --git a/source/lmp/pair_deepmd.h b/source/lmp/pair_deepmd.h
index d04ea97916..bde7745d36 100644
--- a/source/lmp/pair_deepmd.h
+++ b/source/lmp/pair_deepmd.h
@@ -76,6 +76,7 @@ class PairDeepMD : public Pair {
   std::vector<std::string> get_file_content(
       const std::vector<std::string> &models);
   std::vector<std::string> type_names;
+  double ener_unit_cvt_factor, dist_unit_cvt_factor, force_unit_cvt_factor;
 
  protected:
   virtual void allocate();
@@ -132,7 +133,6 @@ class PairDeepMD : public Pair {
   tagint *tagsend, *tagrecv;
   double *stdfsend, *stdfrecv;
   std::vector<int> type_idx_map;
-  double ener_unit_cvt_factor;
 };
 
 }  // namespace LAMMPS_NS
diff --git a/source/lmp/tests/constants.py b/source/lmp/tests/constants.py
new file mode 100644
index 0000000000..016bf66288
--- /dev/null
+++ b/source/lmp/tests/constants.py
@@ -0,0 +1,17 @@
+# SPDX-License-Identifier: LGPL-3.0-or-later
+# https://github.com/lammps/lammps/blob/1e1311cf401c5fc2614b5d6d0ff3230642b76597/src/update.cpp#L193
+nktv2p = 1.6021765e6
+nktv2p_real = 68568.415
+metal2real = 23.060549
+
+dist_metal2real = 1.0
+ener_metal2real = 23.060549
+force_metal2real = ener_metal2real / dist_metal2real
+mass_metal2real = 1.0
+charge_metal2real = 1.0
+
+dist_metal2si = 1.0e-10
+ener_metal2si = 1.3806504e-23 / 8.617343e-5
+force_metal2si = ener_metal2si / dist_metal2si
+mass_metal2si = 1e-3 / 6.02214e23
+charge_metal2si = 1.6021765e-19
diff --git a/source/lmp/tests/test_deeptensor.py b/source/lmp/tests/test_deeptensor.py
index 8c7d7d9640..3e684b386e 100644
--- a/source/lmp/tests/test_deeptensor.py
+++ b/source/lmp/tests/test_deeptensor.py
@@ -6,6 +6,7 @@
     Path,
 )
 
+import constants
 import numpy as np
 import pytest
 from lammps import (
@@ -23,6 +24,7 @@
 pb_file2 = Path(__file__).parent / "deepdipole_new.pb"
 system_file = Path(__file__).parent.parent.parent / "tests"
 data_file = Path(__file__).parent / "data.lmp"
+data_file_si = Path(__file__).parent / "data.si"
 data_type_map_file = Path(__file__).parent / "data_type_map.lmp"
 
 # this is as the same as python and c++ tests, test_deepdipole.py
@@ -75,24 +77,49 @@ def setup_module():
     write_lmp_data(box, coord, type_OH, data_file)
     # TODO
     # write_lmp_data(box, coord, type_HO, data_type_map_file)
+    write_lmp_data(
+        box * constants.dist_metal2si,
+        coord * constants.dist_metal2si,
+        type_OH,
+        data_file_si,
+    )
 
 
 def teardown_module():
     os.remove(data_file)
     # os.remove(data_type_map_file)
+    os.remove(data_file_si)
 
 
-def _lammps(data_file) -> PyLammps:
+def _lammps(data_file, units="metal") -> PyLammps:
     lammps = PyLammps()
-    lammps.units("metal")
+    lammps.units(units)
     lammps.boundary("p p p")
     lammps.atom_style("atomic")
-    lammps.neighbor("2.0 bin")
+    if units == "metal" or units == "real":
+        lammps.neighbor("2.0 bin")
+    elif units == "si":
+        lammps.neighbor("2.0e-10 bin")
+    else:
+        raise ValueError("units should be metal, real, or si")
     lammps.neigh_modify("every 10 delay 0 check no")
     lammps.read_data(data_file.resolve())
-    lammps.mass("1 16")
-    lammps.mass("2 2")
-    lammps.timestep(0.0005)
+    if units == "metal" or units == "real":
+        lammps.mass("1 16")
+        lammps.mass("2 2")
+    elif units == "si":
+        lammps.mass("1 %.10e" % (16 * constants.mass_metal2si))
+        lammps.mass("2 %.10e" % (2 * constants.mass_metal2si))
+    else:
+        raise ValueError("units should be metal, real, or si")
+    if units == "metal":
+        lammps.timestep(0.0005)
+    elif units == "real":
+        lammps.timestep(0.5)
+    elif units == "si":
+        lammps.timestep(5e-16)
+    else:
+        raise ValueError("units should be metal, real, or si")
     lammps.fix("1 all nve")
     return lammps
 
@@ -109,6 +136,13 @@ def lammps():
 #    yield _lammps(data_file=data_type_map_file)
 
 
+@pytest.fixture
+def lammps_si():
+    lmp = _lammps(data_file=data_file_si, units="si")
+    yield lmp
+    lmp.close()
+
+
 def test_compute_deeptensor_atom(lammps):
     lammps.pair_style(f"deepmd {pb_file.resolve()}")
     lammps.pair_coeff("* *")
@@ -120,3 +154,16 @@ def test_compute_deeptensor_atom(lammps):
     assert np.array(lammps.variables["tensor"].value) == pytest.approx(
         expected_d[idx_map]
     )
+
+
+def test_compute_deeptensor_atom_si(lammps_si):
+    lammps_si.pair_style(f"deepmd {pb_file.resolve()}")
+    lammps_si.pair_coeff("* *")
+    lammps_si.compute(f"tensor all deeptensor/atom {pb_file2.resolve()}")
+    lammps_si.variable("tensor atom c_tensor[1]")
+    lammps_si.dump("1 all custom 1 dump id c_tensor[1]")
+    lammps_si.run(0)
+    idx_map = lammps_si.lmp.numpy.extract_atom("id") - 1
+    assert np.array(lammps_si.variables["tensor"].value) == pytest.approx(
+        expected_d[idx_map] * constants.dist_metal2si
+    )
diff --git a/source/lmp/tests/test_dplr.py b/source/lmp/tests/test_dplr.py
index ceebb71310..9c8f1c0d4f 100644
--- a/source/lmp/tests/test_dplr.py
+++ b/source/lmp/tests/test_dplr.py
@@ -6,6 +6,7 @@
     Path,
 )
 
+import constants
 import numpy as np
 import pytest
 from lammps import (
@@ -20,6 +21,7 @@
 dipole_pbtxt_file = Path(__file__).parent / "lrdipole.pbtxt"
 dipole_pb_file = Path(__file__).parent / "lrdipole.pb"
 data_file = Path(__file__).parent / "data.lmp"
+data_file_si = Path(__file__).parent / "data.si"
 data_type_map_file = Path(__file__).parent / "data_type_map.lmp"
 
 # this is as the same as python and c++ tests, test_deeppot_a.py
@@ -261,9 +263,6 @@
 mesh = 10
 
 
-# https://github.com/lammps/lammps/blob/1e1311cf401c5fc2614b5d6d0ff3230642b76597/src/update.cpp#L193
-nktv2p = 1.6021765e6
-
 sp.check_output(
     "{} -m deepmd convert-from pbtxt -i {} -o {}".format(
         sys.executable,
@@ -280,21 +279,45 @@ def setup_module():
     write_lmp_data_full(
         box, coord, mol_list, type_HO, charge, data_type_map_file, bond_list, mass_list
     )
+    write_lmp_data_full(
+        box * constants.dist_metal2si,
+        coord * constants.dist_metal2si,
+        mol_list,
+        type_OH,
+        charge * constants.charge_metal2si,
+        data_file_si,
+        bond_list,
+        mass_list * constants.mass_metal2si,
+    )
 
 
 def teardown_module():
     os.remove(data_file)
+    os.remove(data_type_map_file)
+    os.remove(data_file_si)
 
 
-def _lammps(data_file, exclude_type="1 3") -> PyLammps:
+def _lammps(data_file, exclude_type="1 3", units="metal") -> PyLammps:
     lammps = PyLammps()
-    lammps.units("metal")
+    lammps.units(units)
     lammps.boundary("p p p")
     lammps.atom_style("full")
-    lammps.neighbor("0.2 bin")
+    if units == "metal" or units == "real":
+        lammps.neighbor("0.2 bin")
+    elif units == "si":
+        lammps.neighbor("2.0e-11 bin")
+    else:
+        raise ValueError("units should be metal, real, or si")
     lammps.neigh_modify("every 1 delay 0 check no exclude type " + exclude_type)
     lammps.read_data(data_file.resolve())
-    lammps.timestep(0.0005)
+    if units == "metal":
+        lammps.timestep(0.0005)
+    elif units == "real":
+        lammps.timestep(0.5)
+    elif units == "si":
+        lammps.timestep(5e-16)
+    else:
+        raise ValueError("units should be metal, real, or si")
     lammps.fix("1 all nve")
     return lammps
 
@@ -313,6 +336,13 @@ def lammps_type_map():
     lmp.close()
 
 
+@pytest.fixture
+def lammps_si():
+    lmp = _lammps(data_file=data_file_si, units="si")
+    yield lmp
+    lmp.close()
+
+
 def test_pair_deepmd_sr(lammps):
     lammps.pair_style(f"deepmd {pb_file.resolve()}")
     lammps.pair_coeff("* *")
@@ -349,7 +379,7 @@ def test_pair_deepmd_sr_virial(lammps):
     for ii in range(9):
         assert np.array(lammps.variables[f"virial{ii}"].value)[
             idx_list
-        ] / nktv2p == pytest.approx(expected_v_sr[:, ii])
+        ] / constants.nktv2p == pytest.approx(expected_v_sr[:, ii])
     os.remove("dump")
 
 
@@ -502,3 +532,37 @@ def test_pair_deepmd_lr_type_map(lammps_type_map):
                 expected_f_lr[lammps_type_map.atoms[ii].id - 1]
             )
     lammps_type_map.run(1)
+
+
+def test_pair_deepmd_lr_si(lammps_si):
+    lammps_si.pair_style(f"deepmd {pb_file.resolve()}")
+    lammps_si.pair_coeff("* *")
+    lammps_si.bond_style("zero")
+    lammps_si.bond_coeff("*")
+    lammps_si.special_bonds("lj/coul 1 1 1 angle no")
+    lammps_si.kspace_style("pppm/dplr 1e-5")
+    lammps_si.kspace_modify(
+        f"gewald {beta / constants.dist_metal2si:.6e} diff ik mesh {mesh:d} {mesh:d} {mesh:d}"
+    )
+    lammps_si.fix(
+        f"0 all dplr model {pb_file.resolve()} type_associate 1 3 bond_type 1"
+    )
+    lammps_si.fix_modify("0 virial yes")
+    lammps_si.run(0)
+    for ii in range(8):
+        if lammps_si.atoms[ii].id > 6:
+            assert lammps_si.atoms[ii].position == pytest.approx(
+                expected_WC[lammps_si.atoms[ii].id - 7] * constants.dist_metal2si
+            )
+    assert lammps_si.eval("elong") == pytest.approx(
+        expected_e_kspace * constants.ener_metal2si
+    )
+    assert lammps_si.eval("pe") == pytest.approx(
+        expected_e_lr * constants.ener_metal2si
+    )
+    for ii in range(8):
+        if lammps_si.atoms[ii].id <= 6:
+            assert lammps_si.atoms[ii].force == pytest.approx(
+                expected_f_lr[lammps_si.atoms[ii].id - 1] * constants.force_metal2si
+            )
+    lammps_si.run(1)
diff --git a/source/lmp/tests/test_lammps.py b/source/lmp/tests/test_lammps.py
index 78eaf7ea4e..028b403abf 100644
--- a/source/lmp/tests/test_lammps.py
+++ b/source/lmp/tests/test_lammps.py
@@ -6,6 +6,7 @@
     Path,
 )
 
+import constants
 import numpy as np
 import pytest
 from lammps import (
@@ -23,6 +24,7 @@
 pb_file2 = Path(__file__).parent / "graph2.pb"
 system_file = Path(__file__).parent.parent.parent / "tests"
 data_file = Path(__file__).parent / "data.lmp"
+data_file_si = Path(__file__).parent / "data.si"
 data_type_map_file = Path(__file__).parent / "data_type_map.lmp"
 md_file = Path(__file__).parent / "md.out"
 
@@ -215,10 +217,6 @@
 type_OH = np.array([1, 2, 2, 1, 2, 2])
 type_HO = np.array([2, 1, 1, 2, 1, 1])
 
-# https://github.com/lammps/lammps/blob/1e1311cf401c5fc2614b5d6d0ff3230642b76597/src/update.cpp#L193
-nktv2p = 1.6021765e6
-nktv2p_real = 68568.415
-metal2real = 23.060549
 
 sp.check_output(
     "{} -m deepmd convert-from pbtxt -i {} -o {}".format(
@@ -239,6 +237,12 @@
 def setup_module():
     write_lmp_data(box, coord, type_OH, data_file)
     write_lmp_data(box, coord, type_HO, data_type_map_file)
+    write_lmp_data(
+        box * constants.dist_metal2si,
+        coord * constants.dist_metal2si,
+        type_OH,
+        data_file_si,
+    )
 
 
 def teardown_module():
@@ -251,17 +255,30 @@ def _lammps(data_file, units="metal") -> PyLammps:
     lammps.units(units)
     lammps.boundary("p p p")
     lammps.atom_style("atomic")
-    lammps.neighbor("2.0 bin")
+    if units == "metal" or units == "real":
+        lammps.neighbor("2.0 bin")
+    elif units == "si":
+        lammps.neighbor("2.0e-10 bin")
+    else:
+        raise ValueError("units should be metal, real, or si")
     lammps.neigh_modify("every 10 delay 0 check no")
     lammps.read_data(data_file.resolve())
-    lammps.mass("1 16")
-    lammps.mass("2 2")
+    if units == "metal" or units == "real":
+        lammps.mass("1 16")
+        lammps.mass("2 2")
+    elif units == "si":
+        lammps.mass("1 %.10e" % (16 * constants.mass_metal2si))
+        lammps.mass("2 %.10e" % (2 * constants.mass_metal2si))
+    else:
+        raise ValueError("units should be metal, real, or si")
     if units == "metal":
         lammps.timestep(0.0005)
     elif units == "real":
         lammps.timestep(0.5)
+    elif units == "si":
+        lammps.timestep(5e-16)
     else:
-        raise ValueError("units should be metal or real")
+        raise ValueError("units should be metal, real, or si")
     lammps.fix("1 all nve")
     return lammps
 
@@ -287,6 +304,13 @@ def lammps_real():
     lmp.close()
 
 
+@pytest.fixture
+def lammps_si():
+    lmp = _lammps(data_file=data_file_si, units="si")
+    yield lmp
+    lmp.close()
+
+
 def test_pair_deepmd(lammps):
     lammps.pair_style(f"deepmd {pb_file.resolve()}")
     lammps.pair_coeff("* *")
@@ -319,7 +343,7 @@ def test_pair_deepmd_virial(lammps):
     for ii in range(9):
         assert np.array(
             lammps.variables[f"virial{ii}"].value
-        ) / nktv2p == pytest.approx(expected_v[idx_map, ii])
+        ) / constants.nktv2p == pytest.approx(expected_v[idx_map, ii])
 
 
 def test_pair_deepmd_model_devi(lammps):
@@ -374,7 +398,7 @@ def test_pair_deepmd_model_devi_virial(lammps):
     for ii in range(9):
         assert np.array(
             lammps.variables[f"virial{ii}"].value
-        ) / nktv2p == pytest.approx(expected_v[idx_map, ii])
+        ) / constants.nktv2p == pytest.approx(expected_v[idx_map, ii])
     # load model devi
     md = np.loadtxt(md_file.resolve())
     expected_md_f = np.linalg.norm(np.std([expected_f, expected_f2], axis=0), axis=1)
@@ -472,10 +496,12 @@ def test_pair_deepmd_real(lammps_real):
     lammps_real.pair_style(f"deepmd {pb_file.resolve()}")
     lammps_real.pair_coeff("* *")
     lammps_real.run(0)
-    assert lammps_real.eval("pe") == pytest.approx(expected_e * metal2real)
+    assert lammps_real.eval("pe") == pytest.approx(
+        expected_e * constants.ener_metal2real
+    )
     for ii in range(6):
         assert lammps_real.atoms[ii].force == pytest.approx(
-            expected_f[lammps_real.atoms[ii].id - 1] * metal2real
+            expected_f[lammps_real.atoms[ii].id - 1] * constants.force_metal2real
         )
     lammps_real.run(1)
 
@@ -491,16 +517,20 @@ def test_pair_deepmd_virial_real(lammps_real):
         "1 all custom 1 dump id " + " ".join([f"v_virial{ii}" for ii in range(9)])
     )
     lammps_real.run(0)
-    assert lammps_real.eval("pe") == pytest.approx(expected_e * metal2real)
+    assert lammps_real.eval("pe") == pytest.approx(
+        expected_e * constants.ener_metal2real
+    )
     for ii in range(6):
         assert lammps_real.atoms[ii].force == pytest.approx(
-            expected_f[lammps_real.atoms[ii].id - 1] * metal2real
+            expected_f[lammps_real.atoms[ii].id - 1] * constants.force_metal2real
         )
     idx_map = lammps_real.lmp.numpy.extract_atom("id") - 1
     for ii in range(9):
         assert np.array(
             lammps_real.variables[f"virial{ii}"].value
-        ) / nktv2p_real == pytest.approx(expected_v[idx_map, ii] * metal2real)
+        ) / constants.nktv2p_real == pytest.approx(
+            expected_v[idx_map, ii] * constants.ener_metal2real
+        )
 
 
 def test_pair_deepmd_model_devi_real(lammps_real):
@@ -511,25 +541,27 @@ def test_pair_deepmd_model_devi_real(lammps_real):
     )
     lammps_real.pair_coeff("* *")
     lammps_real.run(0)
-    assert lammps_real.eval("pe") == pytest.approx(expected_e * metal2real)
+    assert lammps_real.eval("pe") == pytest.approx(
+        expected_e * constants.ener_metal2real
+    )
     for ii in range(6):
         assert lammps_real.atoms[ii].force == pytest.approx(
-            expected_f[lammps_real.atoms[ii].id - 1] * metal2real
+            expected_f[lammps_real.atoms[ii].id - 1] * constants.force_metal2real
         )
     # load model devi
     md = np.loadtxt(md_file.resolve())
     expected_md_f = np.linalg.norm(np.std([expected_f, expected_f2], axis=0), axis=1)
-    assert md[7:] == pytest.approx(expected_md_f * metal2real)
-    assert md[4] == pytest.approx(np.max(expected_md_f) * metal2real)
-    assert md[5] == pytest.approx(np.min(expected_md_f) * metal2real)
-    assert md[6] == pytest.approx(np.mean(expected_md_f) * metal2real)
+    assert md[7:] == pytest.approx(expected_md_f * constants.force_metal2real)
+    assert md[4] == pytest.approx(np.max(expected_md_f) * constants.force_metal2real)
+    assert md[5] == pytest.approx(np.min(expected_md_f) * constants.force_metal2real)
+    assert md[6] == pytest.approx(np.mean(expected_md_f) * constants.force_metal2real)
     expected_md_v = (
         np.std([np.sum(expected_v, axis=0), np.sum(expected_v2, axis=0)], axis=0) / 6
     )
-    assert md[1] == pytest.approx(np.max(expected_md_v) * metal2real)
-    assert md[2] == pytest.approx(np.min(expected_md_v) * metal2real)
+    assert md[1] == pytest.approx(np.max(expected_md_v) * constants.ener_metal2real)
+    assert md[2] == pytest.approx(np.min(expected_md_v) * constants.ener_metal2real)
     assert md[3] == pytest.approx(
-        np.sqrt(np.mean(np.square(expected_md_v))) * metal2real
+        np.sqrt(np.mean(np.square(expected_md_v))) * constants.ener_metal2real
     )
 
 
@@ -548,30 +580,34 @@ def test_pair_deepmd_model_devi_virial_real(lammps_real):
         "1 all custom 1 dump id " + " ".join([f"v_virial{ii}" for ii in range(9)])
     )
     lammps_real.run(0)
-    assert lammps_real.eval("pe") == pytest.approx(expected_e * metal2real)
+    assert lammps_real.eval("pe") == pytest.approx(
+        expected_e * constants.ener_metal2real
+    )
     for ii in range(6):
         assert lammps_real.atoms[ii].force == pytest.approx(
-            expected_f[lammps_real.atoms[ii].id - 1] * metal2real
+            expected_f[lammps_real.atoms[ii].id - 1] * constants.force_metal2real
         )
     idx_map = lammps_real.lmp.numpy.extract_atom("id") - 1
     for ii in range(9):
         assert np.array(
             lammps_real.variables[f"virial{ii}"].value
-        ) / nktv2p_real == pytest.approx(expected_v[idx_map, ii] * metal2real)
+        ) / constants.nktv2p_real == pytest.approx(
+            expected_v[idx_map, ii] * constants.ener_metal2real
+        )
     # load model devi
     md = np.loadtxt(md_file.resolve())
     expected_md_f = np.linalg.norm(np.std([expected_f, expected_f2], axis=0), axis=1)
-    assert md[7:] == pytest.approx(expected_md_f * metal2real)
-    assert md[4] == pytest.approx(np.max(expected_md_f) * metal2real)
-    assert md[5] == pytest.approx(np.min(expected_md_f) * metal2real)
-    assert md[6] == pytest.approx(np.mean(expected_md_f) * metal2real)
+    assert md[7:] == pytest.approx(expected_md_f * constants.force_metal2real)
+    assert md[4] == pytest.approx(np.max(expected_md_f) * constants.force_metal2real)
+    assert md[5] == pytest.approx(np.min(expected_md_f) * constants.force_metal2real)
+    assert md[6] == pytest.approx(np.mean(expected_md_f) * constants.force_metal2real)
     expected_md_v = (
         np.std([np.sum(expected_v, axis=0), np.sum(expected_v2, axis=0)], axis=0) / 6
     )
-    assert md[1] == pytest.approx(np.max(expected_md_v) * metal2real)
-    assert md[2] == pytest.approx(np.min(expected_md_v) * metal2real)
+    assert md[1] == pytest.approx(np.max(expected_md_v) * constants.ener_metal2real)
+    assert md[2] == pytest.approx(np.min(expected_md_v) * constants.ener_metal2real)
     assert md[3] == pytest.approx(
-        np.sqrt(np.mean(np.square(expected_md_v))) * metal2real
+        np.sqrt(np.mean(np.square(expected_md_v))) * constants.ener_metal2real
     )
 
 
@@ -582,32 +618,34 @@ def test_pair_deepmd_model_devi_atomic_relative_real(lammps_real):
             pb_file.resolve(),
             pb_file2.resolve(),
             md_file.resolve(),
-            relative * metal2real,
+            relative * constants.force_metal2real,
         )
     )
     lammps_real.pair_coeff("* *")
     lammps_real.run(0)
-    assert lammps_real.eval("pe") == pytest.approx(expected_e * metal2real)
+    assert lammps_real.eval("pe") == pytest.approx(
+        expected_e * constants.ener_metal2real
+    )
     for ii in range(6):
         assert lammps_real.atoms[ii].force == pytest.approx(
-            expected_f[lammps_real.atoms[ii].id - 1] * metal2real
+            expected_f[lammps_real.atoms[ii].id - 1] * constants.force_metal2real
         )
     # load model devi
     md = np.loadtxt(md_file.resolve())
     norm = np.linalg.norm(np.mean([expected_f, expected_f2], axis=0), axis=1)
     expected_md_f = np.linalg.norm(np.std([expected_f, expected_f2], axis=0), axis=1)
     expected_md_f /= norm + relative
-    assert md[7:] == pytest.approx(expected_md_f * metal2real)
-    assert md[4] == pytest.approx(np.max(expected_md_f) * metal2real)
-    assert md[5] == pytest.approx(np.min(expected_md_f) * metal2real)
-    assert md[6] == pytest.approx(np.mean(expected_md_f) * metal2real)
+    assert md[7:] == pytest.approx(expected_md_f * constants.force_metal2real)
+    assert md[4] == pytest.approx(np.max(expected_md_f) * constants.force_metal2real)
+    assert md[5] == pytest.approx(np.min(expected_md_f) * constants.force_metal2real)
+    assert md[6] == pytest.approx(np.mean(expected_md_f) * constants.force_metal2real)
     expected_md_v = (
         np.std([np.sum(expected_v, axis=0), np.sum(expected_v2, axis=0)], axis=0) / 6
     )
-    assert md[1] == pytest.approx(np.max(expected_md_v) * metal2real)
-    assert md[2] == pytest.approx(np.min(expected_md_v) * metal2real)
+    assert md[1] == pytest.approx(np.max(expected_md_v) * constants.ener_metal2real)
+    assert md[2] == pytest.approx(np.min(expected_md_v) * constants.ener_metal2real)
     assert md[3] == pytest.approx(
-        np.sqrt(np.mean(np.square(expected_md_v))) * metal2real
+        np.sqrt(np.mean(np.square(expected_md_v))) * constants.ener_metal2real
     )
 
 
@@ -618,22 +656,24 @@ def test_pair_deepmd_model_devi_atomic_relative_v_real(lammps_real):
             pb_file.resolve(),
             pb_file2.resolve(),
             md_file.resolve(),
-            relative * metal2real,
+            relative * constants.ener_metal2real,
         )
     )
     lammps_real.pair_coeff("* *")
     lammps_real.run(0)
-    assert lammps_real.eval("pe") == pytest.approx(expected_e * metal2real)
+    assert lammps_real.eval("pe") == pytest.approx(
+        expected_e * constants.ener_metal2real
+    )
     for ii in range(6):
         assert lammps_real.atoms[ii].force == pytest.approx(
-            expected_f[lammps_real.atoms[ii].id - 1] * metal2real
+            expected_f[lammps_real.atoms[ii].id - 1] * constants.force_metal2real
         )
     md = np.loadtxt(md_file.resolve())
     expected_md_f = np.linalg.norm(np.std([expected_f, expected_f2], axis=0), axis=1)
-    assert md[7:] == pytest.approx(expected_md_f * metal2real)
-    assert md[4] == pytest.approx(np.max(expected_md_f) * metal2real)
-    assert md[5] == pytest.approx(np.min(expected_md_f) * metal2real)
-    assert md[6] == pytest.approx(np.mean(expected_md_f) * metal2real)
+    assert md[7:] == pytest.approx(expected_md_f * constants.force_metal2real)
+    assert md[4] == pytest.approx(np.max(expected_md_f) * constants.force_metal2real)
+    assert md[5] == pytest.approx(np.min(expected_md_f) * constants.force_metal2real)
+    assert md[6] == pytest.approx(np.mean(expected_md_f) * constants.force_metal2real)
     expected_md_v = (
         np.std([np.sum(expected_v, axis=0), np.sum(expected_v2, axis=0)], axis=0) / 6
     )
@@ -644,8 +684,20 @@ def test_pair_deepmd_model_devi_atomic_relative_v_real(lammps_real):
         / 6
     )
     expected_md_v /= norm + relative
-    assert md[1] == pytest.approx(np.max(expected_md_v) * metal2real)
-    assert md[2] == pytest.approx(np.min(expected_md_v) * metal2real)
+    assert md[1] == pytest.approx(np.max(expected_md_v) * constants.ener_metal2real)
+    assert md[2] == pytest.approx(np.min(expected_md_v) * constants.ener_metal2real)
     assert md[3] == pytest.approx(
-        np.sqrt(np.mean(np.square(expected_md_v))) * metal2real
+        np.sqrt(np.mean(np.square(expected_md_v))) * constants.ener_metal2real
     )
+
+
+def test_pair_deepmd_si(lammps_si):
+    lammps_si.pair_style(f"deepmd {pb_file.resolve()}")
+    lammps_si.pair_coeff("* *")
+    lammps_si.run(0)
+    assert lammps_si.eval("pe") == pytest.approx(expected_e * constants.ener_metal2si)
+    for ii in range(6):
+        assert lammps_si.atoms[ii].force == pytest.approx(
+            expected_f[lammps_si.atoms[ii].id - 1] * constants.force_metal2si
+        )
+    lammps_si.run(1)
diff --git a/source/lmp/tests/write_lmp_data.py b/source/lmp/tests/write_lmp_data.py
index 05a2760f4c..12e91764f1 100644
--- a/source/lmp/tests/write_lmp_data.py
+++ b/source/lmp/tests/write_lmp_data.py
@@ -11,13 +11,13 @@ def write_lmp_data(box, coord, type_list, file_name):
         f.write(comment_lmp_data + "\n")
         f.write("%d atoms\n" % (natom))
         f.write("%d atom types\n" % (ntype))
-        f.write(f"{box[0]:.8f} {box[1]:.8f} xlo xhi\n")
-        f.write(f"{box[2]:.8f} {box[3]:.8f} ylo yhi\n")
-        f.write(f"{box[4]:.8f} {box[5]:.8f} zlo zhi\n")
-        f.write(f"{box[6]:.8f} {box[7]:.8f} {box[8]:.8f} xy xz yz\n\nAtoms\n\n")
+        f.write(f"{box[0]:.10e} {box[1]:.10e} xlo xhi\n")
+        f.write(f"{box[2]:.10e} {box[3]:.10e} ylo yhi\n")
+        f.write(f"{box[4]:.10e} {box[5]:.10e} zlo zhi\n")
+        f.write(f"{box[6]:.10e} {box[7]:.10e} {box[8]:.10e} xy xz yz\n\nAtoms\n\n")
         for i in range(natom):
             f.write(
-                "%d %d %.8f %.8f %.8f\n"
+                "%d %d %.10e %.10e %.10e\n"
                 % (i + 1, type_list[i], coord[i][0], coord[i][1], coord[i][2])
             )
         f.write("\n")
@@ -38,17 +38,17 @@ def write_lmp_data_full(
         f.write("%d atom types\n" % (ntype))
         f.write("%d bonds\n" % (nbond_list.sum()))
         f.write("%d bond types\n" % (nbond_type))
-        f.write(f"{box[0]:.8f} {box[1]:.8f} xlo xhi\n")
-        f.write(f"{box[2]:.8f} {box[3]:.8f} ylo yhi\n")
-        f.write(f"{box[4]:.8f} {box[5]:.8f} zlo zhi\n")
-        f.write(f"{box[6]:.8f} {box[7]:.8f} {box[8]:.8f} xy xz yz\n")
+        f.write(f"{box[0]:.10e} {box[1]:.10e} xlo xhi\n")
+        f.write(f"{box[2]:.10e} {box[3]:.10e} ylo yhi\n")
+        f.write(f"{box[4]:.10e} {box[5]:.10e} zlo zhi\n")
+        f.write(f"{box[6]:.10e} {box[7]:.10e} {box[8]:.10e} xy xz yz\n")
         f.write("\nMasses\n\n")
         for i in range(3):
-            f.write(f"{i+1:d} {mass_list[i]:.6f}\n")
+            f.write(f"{i+1:d} {mass_list[i]:.10e}\n")
         f.write("\nAtoms\n\n")
         for i in range(natom):
             f.write(
-                "%d %d %d %d %.8f %.8f %.8f\n"
+                "%d %d %d %.10e %.10e %.10e %.10e\n"
                 % (
                     i + 1,
                     mol_list[i],