initial implementation of automated unit conversion.

this includes a tester program and implementation into pair style tersoff

src/MANYBODY/pair_tersoff.cpp

@@ -49,6 +49,7 @@ PairTersoff::PairTersoff(LAMMPS *lmp) : Pair(lmp)
   restartinfo = 0;
   one_coeff = 1;
   manybody_flag = 1;
+  unit_convert_flag = utils::get_supported_conversions(utils::ENERGY);
 
   nelements = 0;
   elements = NULL;
@@ -400,8 +401,16 @@ void PairTersoff::read_file(char *file)
   // open file on proc 0
 
   if (comm->me == 0) {
-    PotentialFileReader reader(lmp, file, "Tersoff");
+    PotentialFileReader reader(lmp, file, "Tersoff", unit_convert_flag);
-    char * line;
+    char *line;
+
+    // transparently convert units for supported conversions
+
+    int unit_convert = reader.get_unit_convert();
+    double conversion_factor = utils::get_conversion_factor(utils::ENERGY,
+                                                            unit_convert);
+    printf("unit_conver=%d\n",unit_convert);
+    printf("conversion_factor=%g\n",conversion_factor);
 
     while((line = reader.next_line(NPARAMS_PER_LINE))) {
       try {
@@ -453,6 +462,11 @@ void PairTersoff::read_file(char *file)
         params[nparams].lam1      = values.next_double();
         params[nparams].biga      = values.next_double();
         params[nparams].powermint = int(params[nparams].powerm);
+
+        if (unit_convert) {
+          params[nparams].biga *= conversion_factor;
+          params[nparams].bigb *= conversion_factor;
+        }
       } catch (TokenizerException & e) {
         error->one(FLERR, e.what());
       }

src/pair.cpp

@@ -66,6 +66,7 @@ Pair::Pair(LAMMPS *lmp) : Pointers(lmp)
   no_virial_fdotr_compute = 0;
   writedata = 0;
   ghostneigh = 0;
+  unit_convert_flag = utils::NOCONVERT;
 
   nextra = 0;
   pvector = NULL;

src/pair.h

@@ -53,6 +53,7 @@ class Pair : protected Pointers {
   int respa_enable;              // 1 if inner/middle/outer rRESPA routines
   int one_coeff;                 // 1 if allows only one coeff * * call
   int manybody_flag;             // 1 if a manybody potential
+  int unit_convert_flag;         // value != 0 indicates support for unit conversion.
   int no_virial_fdotr_compute;   // 1 if does not invoke virial_fdotr_compute()
   int writedata;                 // 1 if writes coeffs to data file
   int ghostneigh;                // 1 if pair style needs neighbors of ghosts

src/potential_file_reader.cpp

@@ -31,11 +31,13 @@ using namespace LAMMPS_NS;
 
 PotentialFileReader::PotentialFileReader(LAMMPS *lmp,
                                          const std::string &filename,
-                                         const std::string &potential_name) :
+                                         const std::string &potential_name,
+                                         const int auto_convert) :
   Pointers(lmp),
   reader(nullptr),
   filename(filename),
-  filetype(potential_name + " potential")
+  filetype(potential_name + " potential"),
+  unit_convert(auto_convert)
 {
   if (comm->me != 0) {
     error->one(FLERR, "FileReader should only be called by proc 0!");
@@ -155,9 +157,21 @@ TextFileReader * PotentialFileReader::open_potential(const std::string& path) {
       utils::logmesg(lmp, fmt::format("Reading potential file {} with DATE: {}\n", filename, date));
     }
 
-    if (!units.empty() && (units != unit_style)) {
+    if (units.empty()) {
-      lmp->error->one(FLERR, fmt::format("Potential file {} requires {} units "
+      unit_convert == utils::NOCONVERT;
-                                         "but {} units are in use",filename, units, unit_style));
+    } else {
+      if (units == unit_style) {
+        unit_convert = utils::NOCONVERT;
+      } else {
+        if ((units == "metal") && (unit_style == "real") && (unit_convert | utils::METAL2REAL)) {
+          unit_convert = utils::METAL2REAL;
+        } else if ((units == "real") && (unit_style == "metal") && (unit_convert | utils::REAL2METAL)) {
+          unit_convert = utils::REAL2METAL;
+        } else {
+          lmp->error->one(FLERR, fmt::format("Potential file {} requires {} units "
+                                             "but {} units are in use",filename, units, unit_style));
+        }
+      }
     }
 
     return new TextFileReader(filepath, filetype);

src/potential_file_reader.h

@@ -28,22 +28,25 @@ namespace LAMMPS_NS
 {
   class PotentialFileReader : protected Pointers {
   protected:
-    TextFileReader * reader;
+    TextFileReader *reader;
     std::string filename;
     std::string filetype;
+    int unit_convert;
 
-    TextFileReader * open_potential(const std::string& path);
+    TextFileReader *open_potential(const std::string& path);
 
   public:
-    PotentialFileReader(class LAMMPS *lmp, const std::string &filename, const std::string &potential_name);
+    PotentialFileReader(class LAMMPS *lmp, const std::string &filename,
+                        const std::string &potential_name,
+                        const int auto_convert = 0);
     virtual ~PotentialFileReader();
 
     void ignore_comments(bool value);
 
     void skip_line();
-    char * next_line(int nparams = 0);
+    char *next_line(int nparams = 0);
-    void next_dvector(double * list, int n);
+    void next_dvector(double *list, int n);
-    ValueTokenizer next_values(int nparams, const std::string & separators = TOKENIZER_DEFAULT_SEPARATORS);
+    ValueTokenizer next_values(int nparams, const std::string &separators = TOKENIZER_DEFAULT_SEPARATORS);
 
     // convenience functions
     double next_double();
@@ -51,6 +54,9 @@ namespace LAMMPS_NS
     tagint next_tagint();
     bigint next_bigint();
     std::string next_string();
+
+    // unit conversion info
+    int get_unit_convert() const { return unit_convert; }
   };
 
 } // namespace LAMMPS_NS

src/utils.cpp

@@ -583,6 +583,36 @@ std::string utils::get_potential_units(const std::string & path, const std::stri
   return "";
 }
 
+/* ----------------------------------------------------------------------
+   return bitmask of supported conversions for a given property
+------------------------------------------------------------------------- */
+int utils::get_supported_conversions(const int property)
+{
+  if (property == ENERGY) {
+    return METAL2REAL | REAL2METAL;
+  }
+  return NOCONVERT;
+}
+
+/* ----------------------------------------------------------------------
+   return conversion factor for a given property and conversion setting
+   return 0.0 if unknown.
+------------------------------------------------------------------------- */
+
+double utils::get_conversion_factor(const int property, const int conversion)
+{
+  if (property == ENERGY) {
+    if (conversion == NOCONVERT) {
+      return 1.0;
+    } else if (conversion == METAL2REAL) {
+      return 23.060549;
+    } else if (conversion == REAL2METAL) {
+      return 0.04336410204284381954;
+    }
+  }
+  return 0.0;
+}
+
 /* ------------------------------------------------------------------ */
 
 extern "C" {

src/utils.h

@@ -240,6 +240,24 @@ namespace LAMMPS_NS {
      * \return UNITS field if present
      */
     std::string get_potential_units(const std::string & path, const std::string & potential_name);
+
+    enum { NOCONVERT = 0, METAL2REAL = 1, REAL2METAL = 1<<1 };
+    enum { UNKNOWN = 0, ENERGY };
+
+    /**
+     * \brief Return bitmask of available conversion factors for a given propert
+     * \param property property to be converted
+     * \return bitmask indicating available conversions
+     */
+    int get_supported_conversions(const int property);
+
+    /**
+     * \brief Return unit conversion factor for given property and selected from/to units
+     * \param property property to be converted
+     * \param conversion constant indicating the conversion
+     * \return conversion factor
+     */
+    double get_conversion_factor(const int property, const int conversion);
   }
 }
 

unittest/formats/CMakeLists.txt

@@ -3,6 +3,11 @@ add_executable(test_atom_styles test_atom_styles.cpp)
 target_link_libraries(test_atom_styles PRIVATE lammps GTest::GMock GTest::GTest)
 add_test(NAME AtomStyles COMMAND test_atom_styles WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
 
+add_executable(test_pair_unit_convert test_pair_unit_convert.cpp)
+target_link_libraries(test_pair_unit_convert PRIVATE lammps GTest::GMock GTest::GTest)
+add_test(NAME PairUnitConvert COMMAND test_pair_unit_convert WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
+set_tests_properties(PairUnitConvert PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR}")
+
 add_executable(test_potential_file_reader test_potential_file_reader.cpp)
 target_link_libraries(test_potential_file_reader PRIVATE lammps GTest::GMock GTest::GTest)
 add_test(NAME PotentialFileReader COMMAND test_potential_file_reader WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})

unittest/formats/test_pair_unit_convert.cpp

@@ -0,0 +1,206 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "atom.h"
+#include "force.h"
+#include "info.h"
+#include "input.h"
+#include "lammps.h"
+#include "pair.h"
+#include "utils.h"
+#include "gmock/gmock.h"
+#include "gtest/gtest.h"
+
+#include <cstdio>
+#include <cstring>
+#include <mpi.h>
+
+// whether to print verbose output (i.e. not capturing LAMMPS screen output).
+bool verbose = false;
+
+namespace LAMMPS_NS {
+using ::testing::Eq;
+
+// eV to kcal/mol conversion constant (CODATA 2018)
+const double ev_convert = utils::get_conversion_factor(utils::ENERGY, utils::METAL2REAL);
+const double rel_error  = 1.0e-13;
+
+class PairUnitConvertTest : public ::testing::Test {
+protected:
+    LAMMPS *lmp;
+    Info *info;
+    double fold[4][3];
+
+    void SetUp() override
+    {
+        const char *args[] = {"PairUnitConvertTest", "-log", "none", "-echo", "screen", "-nocite"};
+        char **argv        = (char **)args;
+        int argc           = sizeof(args) / sizeof(char *);
+        if (!verbose) ::testing::internal::CaptureStdout();
+        lmp = new LAMMPS(argc, argv, MPI_COMM_WORLD);
+        if (!verbose) ::testing::internal::GetCapturedStdout();
+        ASSERT_NE(lmp, nullptr);
+        if (!verbose) ::testing::internal::CaptureStdout();
+        info = new Info(lmp);
+        lmp->input->one("units metal");
+        lmp->input->one("dimension 3");
+        lmp->input->one("region box block -4 4 -4 4 -4 4");
+        lmp->input->one("create_box 2 box");
+        lmp->input->one("create_atoms 1 single -1.1  1.2  0.0 units box");
+        lmp->input->one("create_atoms 1 single -1.2 -1.1  0.0 units box");
+        lmp->input->one("create_atoms 2 single  0.9  1.0  0.0 units box");
+        lmp->input->one("create_atoms 2 single  1.0 -0.9  0.0 units box");
+        lmp->input->one("pair_style zero 4.0");
+        lmp->input->one("pair_coeff * *");
+        lmp->input->one("mass * 1.0");
+        lmp->input->one("write_data test_pair_unit_convert.data nocoeff");
+        lmp->input->one("clear");
+        if (!verbose) ::testing::internal::GetCapturedStdout();
+    }
+
+    void TearDown() override
+    {
+        if (!verbose) ::testing::internal::CaptureStdout();
+        delete info;
+        delete lmp;
+        if (!verbose) ::testing::internal::GetCapturedStdout();
+        remove("test_pair_unit_convert.data");
+    }
+};
+
+TEST_F(PairUnitConvertTest, zero)
+{
+    // check if the prerequisite pair style is available
+    if (!info->has_style("pair", "zero")) GTEST_SKIP();
+
+    if (!verbose) ::testing::internal::CaptureStdout();
+    lmp->input->one("units metal");
+    lmp->input->one("read_data test_pair_unit_convert.data");
+    lmp->input->one("pair_style zero 6.0");
+    lmp->input->one("pair_coeff * *");
+    lmp->input->one("run 0 post no");
+    if (!verbose) ::testing::internal::GetCapturedStdout();
+
+    // copy energy and force from first step
+    double eold = lmp->force->pair->eng_vdwl + lmp->force->pair->eng_coul;
+    double **f  = lmp->atom->f;
+    for (int i = 0; i < 4; ++i)
+        for (int j = 0; j < 3; ++j)
+            fold[i][j] = f[i][j];
+
+    if (!verbose) ::testing::internal::CaptureStdout();
+    lmp->input->one("clear");
+    lmp->input->one("units real");
+    lmp->input->one("read_data test_pair_unit_convert.data");
+    lmp->input->one("pair_style zero 6.0");
+    lmp->input->one("pair_coeff * *");
+    lmp->input->one("run 0 post no");
+    if (!verbose) ::testing::internal::GetCapturedStdout();
+
+    double enew = lmp->force->pair->eng_vdwl + lmp->force->pair->eng_coul;
+    EXPECT_NEAR(ev_convert * eold, enew, fabs(enew * rel_error));
+
+    f = lmp->atom->f;
+    for (int i = 0; i < 4; ++i)
+        for (int j = 0; j < 3; ++j)
+            EXPECT_NEAR(ev_convert * fold[i][j], f[i][j], fabs(f[i][j] * rel_error));
+}
+
+TEST_F(PairUnitConvertTest, lj_cut)
+{
+    // check if the prerequisite pair style is available
+    if (!info->has_style("pair", "lj/cut")) GTEST_SKIP();
+
+    if (!verbose) ::testing::internal::CaptureStdout();
+    lmp->input->one("units metal");
+    lmp->input->one("read_data test_pair_unit_convert.data");
+    lmp->input->one("pair_style lj/cut 6.0");
+    lmp->input->one("pair_coeff * * 0.01014286346782117 3.504");
+    lmp->input->one("run 0 post no");
+    if (!verbose) ::testing::internal::GetCapturedStdout();
+
+    // copy energy and force from first step
+    double eold = lmp->force->pair->eng_vdwl + lmp->force->pair->eng_coul;
+    double **f  = lmp->atom->f;
+    for (int i = 0; i < 4; ++i)
+        for (int j = 0; j < 3; ++j)
+            fold[i][j] = f[i][j];
+
+    if (!verbose) ::testing::internal::CaptureStdout();
+    lmp->input->one("clear");
+    lmp->input->one("units real");
+    lmp->input->one("read_data test_pair_unit_convert.data");
+    lmp->input->one("pair_style lj/cut 6.0");
+    lmp->input->one("pair_coeff * * 0.2339 3.504");
+    lmp->input->one("run 0 post no");
+    if (!verbose) ::testing::internal::GetCapturedStdout();
+
+    double enew = lmp->force->pair->eng_vdwl + lmp->force->pair->eng_coul;
+    EXPECT_NEAR(ev_convert * eold, enew, fabs(enew * rel_error));
+
+    f = lmp->atom->f;
+    for (int i = 0; i < 4; ++i)
+        for (int j = 0; j < 3; ++j)
+            EXPECT_NEAR(ev_convert * fold[i][j], f[i][j], fabs(f[i][j] * rel_error));
+}
+
+TEST_F(PairUnitConvertTest, tersoff)
+{
+    // check if the prerequisite pair style is available
+    if (!info->has_style("pair", "tersoff")) GTEST_SKIP();
+
+    if (!verbose) ::testing::internal::CaptureStdout();
+    lmp->input->one("units metal");
+    lmp->input->one("read_data test_pair_unit_convert.data");
+    lmp->input->one("pair_style tersoff");
+    lmp->input->one("pair_coeff * * SiC.tersoff Si C");
+    lmp->input->one("run 0 post no");
+    if (!verbose) ::testing::internal::GetCapturedStdout();
+
+    // copy energy and force from first step
+    double eold = lmp->force->pair->eng_vdwl + lmp->force->pair->eng_coul;
+    double **f  = lmp->atom->f;
+    for (int i = 0; i < 4; ++i)
+        for (int j = 0; j < 3; ++j)
+            fold[i][j] = f[i][j];
+
+    if (!verbose) ::testing::internal::CaptureStdout();
+    lmp->input->one("clear");
+    lmp->input->one("units real");
+    lmp->input->one("read_data test_pair_unit_convert.data");
+    lmp->input->one("pair_style tersoff");
+    lmp->input->one("pair_coeff * * SiC.tersoff Si C");
+    lmp->input->one("run 0 post no");
+    if (!verbose) ::testing::internal::GetCapturedStdout();
+
+    double enew = lmp->force->pair->eng_vdwl + lmp->force->pair->eng_coul;
+    EXPECT_NEAR(ev_convert * eold, enew, fabs(enew * rel_error));
+
+    f = lmp->atom->f;
+    for (int i = 0; i < 4; ++i)
+        for (int j = 0; j < 3; ++j)
+            EXPECT_NEAR(ev_convert * fold[i][j], f[i][j], fabs(f[i][j] * rel_error));
+}
+
+} // namespace LAMMPS_NS
+
+int main(int argc, char **argv)
+{
+    MPI_Init(&argc, &argv);
+    ::testing::InitGoogleMock(&argc, argv);
+    if ((argc > 1) && (strcmp(argv[1], "-v") == 0)) verbose = true;
+
+    int rv = RUN_ALL_TESTS();
+    MPI_Finalize();
+    return rv;
+}