add unit conversion support to pair style tersoff/table

src/USER-MISC/pair_tersoff_table.cpp

@@ -33,6 +33,8 @@
 #include "comm.h"
 #include "memory.h"
 #include "utils.h"
+#include "tokenizer.h"
+#include "potential_file_reader.h"
 
 #include "error.h"
 
@@ -59,6 +61,7 @@ PairTersoffTable::PairTersoffTable(LAMMPS *lmp) : Pair(lmp)
   restartinfo = 0;
   one_coeff = 1;
   manybody_flag = 1;
+  unit_convert_flag = utils::get_supported_conversions(utils::ENERGY);
 
   nelements = 0;
   elements = nullptr;
@@ -840,92 +843,45 @@ double PairTersoffTable::init_one(int i, int j)
 
 void PairTersoffTable::read_file(char *file)
 {
-  int params_per_line = 17;
-  char **words = new char*[params_per_line+1];
-
   memory->sfree(params);
-  params = NULL;
+  params = nullptr;
   nparams = maxparam = 0;
 
   // open file on proc 0
 
-  FILE *fp;
   if (comm->me == 0) {
-    fp = force->open_potential(file);
+    PotentialFileReader reader(lmp, file, "TersoffTable", unit_convert_flag);
-    if (fp == NULL) {
+    char *line;
-      char str[128];
-      snprintf(str,128,"Cannot open Tersoff potential file %s",file);
-      error->one(FLERR,str);
-    }
-  }
     
-  // read each set of params from potential file
+    // transparently convert units for supported conversions
-  // one set of params can span multiple lines
-  // store params if all 3 element tags are in element list
 
-  int n,nwords,ielement,jelement,kelement;
+    int unit_convert = reader.get_unit_convert();
-  char line[MAXLINE],*ptr;
+    double conversion_factor = utils::get_conversion_factor(utils::ENERGY,
-  int eof = 0;
+                                                            unit_convert);
 
-  while (1) {
+    while((line = reader.next_line(NPARAMS_PER_LINE))) {
-    if (comm->me == 0) {
-      ptr = fgets(line,MAXLINE,fp);
-      if (ptr == NULL) {
-        eof = 1;
-        fclose(fp);
-      } else n = strlen(line) + 1;
-    }
-    MPI_Bcast(&eof,1,MPI_INT,0,world);
-    if (eof) break;
-    MPI_Bcast(&n,1,MPI_INT,0,world);
-    MPI_Bcast(line,n,MPI_CHAR,0,world);
 
-    // strip comment, skip line if blank
+      try {
+        ValueTokenizer values(line);
 
-    if ((ptr = strchr(line,'#'))) *ptr = '\0';
+        std::string iname = values.next_string();
-    nwords = utils::count_words(line);
+        std::string jname = values.next_string();
-    if (nwords == 0) continue;
+        std::string kname = values.next_string();
-
-    // concatenate additional lines until have params_per_line words
-
-    while (nwords < params_per_line) {
-      n = strlen(line);
-      if (comm->me == 0) {
-        ptr = fgets(&line[n],MAXLINE-n,fp);
-        if (ptr == NULL) {
-          eof = 1;
-          fclose(fp);
-        } else n = strlen(line) + 1;
-      }
-      MPI_Bcast(&eof,1,MPI_INT,0,world);
-      if (eof) break;
-      MPI_Bcast(&n,1,MPI_INT,0,world);
-      MPI_Bcast(line,n,MPI_CHAR,0,world);
-      if ((ptr = strchr(line,'#'))) *ptr = '\0';
-      nwords = utils::count_words(line);
-    }
-
-    if (nwords != params_per_line)
-      error->all(FLERR,"Incorrect format in Tersoff potential file");
-
-    // words = ptrs to all words in line
-
-    nwords = 0;
-    words[nwords++] = strtok(line," \t\n\r\f");
-    while ((words[nwords++] = strtok(NULL," \t\n\r\f"))) continue;
 
         // ielement,jelement,kelement = 1st args
         // if all 3 args are in element list, then parse this line
         // else skip to next entry in file
 
+        int ielement, jelement, kelement;
+
         for (ielement = 0; ielement < nelements; ielement++)
-      if (strcmp(words[0],elements[ielement]) == 0) break;
+          if (iname == elements[ielement]) break;
         if (ielement == nelements) continue;
         for (jelement = 0; jelement < nelements; jelement++)
-      if (strcmp(words[1],elements[jelement]) == 0) break;
+          if (jname == elements[jelement]) break;
         if (jelement == nelements) continue;
         for (kelement = 0; kelement < nelements; kelement++)
-      if (strcmp(words[2],elements[kelement]) == 0) break;
+          if (kname == elements[kelement]) break;
         if (kelement == nelements) continue;
 
         // load up parameter settings and error check their values
@@ -936,29 +892,35 @@ void PairTersoffTable::read_file(char *file)
                                               "pair:params");
         }
 
+        // some parameters are not used since only Tersoff_2 is implemented
+
         params[nparams].ielement  = ielement;
         params[nparams].jelement  = jelement;
         params[nparams].kelement  = kelement;
-    params[nparams].powerm = atof(words[3]); // not used (only tersoff_2 is implemented)
+        params[nparams].powerm    = values.next_double(); // not used
-    params[nparams].gamma = atof(words[4]); // not used (only tersoff_2 is implemented)
+        params[nparams].gamma     = values.next_double(); // not used
-    params[nparams].lam3 = atof(words[5]); // not used (only tersoff_2 is implemented)
+        params[nparams].lam3      = values.next_double(); // not used
-    params[nparams].c = atof(words[6]);
+        params[nparams].c         = values.next_double();
-    params[nparams].d = atof(words[7]);
+        params[nparams].d         = values.next_double();
-    params[nparams].h = atof(words[8]);
+        params[nparams].h         = values.next_double();
-    params[nparams].powern = atof(words[9]);
+        params[nparams].powern    = values.next_double();
-    params[nparams].beta = atof(words[10]);
+        params[nparams].beta      = values.next_double();
-    params[nparams].lam2 = atof(words[11]);
+        params[nparams].lam2      = values.next_double();
-    params[nparams].bigb = atof(words[12]);
+        params[nparams].bigb      = values.next_double();
-
+        double bigr               = values.next_double();
-    // current implementation is based on functional form
+        double bigd               = values.next_double();
-    // of tersoff_2 as reported in the reference paper
-
-    double bigr = atof(words[13]);
-    double bigd = atof(words[14]);
         params[nparams].cutoffR = bigr - bigd;
         params[nparams].cutoffS = bigr + bigd;
-    params[nparams].lam1 = atof(words[15]);
+        params[nparams].lam1      = values.next_double();
-    params[nparams].biga = atof(words[16]);
+        params[nparams].biga      = values.next_double();
+
+        if (unit_convert) {
+          params[nparams].biga *= conversion_factor;
+          params[nparams].bigb *= conversion_factor;
+        }
+      } catch (TokenizerException & e) {
+        error->one(FLERR, e.what());
+      }
 
       if (params[nparams].c < 0.0 ||
           params[nparams].d < 0.0 ||
@@ -971,15 +933,25 @@ void PairTersoffTable::read_file(char *file)
           params[nparams].cutoffR > params[nparams].cutoffS ||
           params[nparams].lam1 < 0.0 ||
           params[nparams].biga < 0.0
-    ) error->all(FLERR,"Illegal Tersoff parameter");
+        ) error->one(FLERR,"Illegal Tersoff parameter");
 
       // only tersoff_2 parametrization is implemented
+
       if (params[nparams].gamma != 1.0 || params[nparams].lam3 != 0.0)
-      error->all(FLERR,"Current tersoff/table pair_style implements only tersoff_2 parametrization");
+        error->one(FLERR,"Currently the tersoff/table pair_style only "
+                   "implements the Tersoff_2 parametrization");
       nparams++;
     }
+  }
 
-  delete [] words;
+  MPI_Bcast(&nparams, 1, MPI_INT, 0, world);
+  MPI_Bcast(&maxparam, 1, MPI_INT, 0, world);
+
+  if(comm->me != 0) {
+    params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), "pair:params");
+  }
+
+  MPI_Bcast(params, maxparam*sizeof(Param), MPI_BYTE, 0, world);
 }
 
 /* ---------------------------------------------------------------------- */

src/USER-MISC/pair_tersoff_table.h

@@ -40,8 +40,10 @@ class PairTersoffTable : public Pair {
   virtual void compute(int, int);
   void settings(int, char **);
   void coeff(int, char **);
-  double init_one(int, int);
   void init_style();
+  double init_one(int, int);
+
+  static const int NPARAMS_PER_LINE = 17;
 
  protected:
   struct Param {

unittest/formats/test_pair_unit_convert.cpp

@@ -214,6 +214,49 @@ TEST_F(PairUnitConvertTest, tersoff)
             EXPECT_NEAR(ev_convert * fold[i][j], f[i][j], fabs(f[i][j] * rel_error));
 }
 
+TEST_F(PairUnitConvertTest, tersoff_table)
+{
+    // check if the prerequisite pair style is available
+    if (!info->has_style("pair", "tersoff/table")) 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/table");
+    lmp->input->one("pair_coeff * * SiC.tersoff Si C");
+    lmp->input->one("run 0 post no");
+    if (!verbose) ::testing::internal::GetCapturedStdout();
+
+    // copy pressure, energy, and force from first step
+    double pold;
+    lmp->output->thermo->evaluate_keyword("press", &pold);
+    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/table");
+    lmp->input->one("pair_coeff * * SiC.tersoff Si C");
+    lmp->input->one("run 0 post no");
+    if (!verbose) ::testing::internal::GetCapturedStdout();
+
+    double pnew;
+    lmp->output->thermo->evaluate_keyword("press", &pnew);
+    EXPECT_NEAR(pold, p_convert * pnew, fabs(pnew * rel_error));
+    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, sw)
 {
     // check if the prerequisite pair style is available

unittest/formats/test_potential_file_reader.cpp

@@ -23,6 +23,7 @@
 #include "MANYBODY/pair_tersoff_mod_c.h"
 #include "MANYBODY/pair_tersoff_zbl.h"
 #include "MANYBODY/pair_vashishta.h"
+#include "USER-MISC/pair_tersoff_table.h"
 #include "input.h"
 #include "lammps.h"
 #include "potential_file_reader.h"
@@ -50,6 +51,7 @@ const int LAMMPS_NS::PairGW::NPARAMS_PER_LINE;
 const int LAMMPS_NS::PairGWZBL::NPARAMS_PER_LINE;
 const int LAMMPS_NS::PairNb3bHarmonic::NPARAMS_PER_LINE;
 const int LAMMPS_NS::PairVashishta::NPARAMS_PER_LINE;
+const int LAMMPS_NS::PairTersoffTable::NPARAMS_PER_LINE;
 
 class PotentialFileReaderTest : public ::testing::Test {
 protected:
@@ -140,6 +142,17 @@ TEST_F(PotentialFileReaderTest, TersoffModC)
     ASSERT_EQ(utils::count_words(line), PairTersoffMODC::NPARAMS_PER_LINE);
 }
 
+TEST_F(PotentialFileReaderTest, TersoffTable)
+{
+    if (!verbose) ::testing::internal::CaptureStdout();
+    lmp->input->one("units metal");
+    PotentialFileReader reader(lmp, "Si.tersoff", "TersoffTable");
+    if (!verbose) ::testing::internal::GetCapturedStdout();
+
+    auto line = reader.next_line(PairTersoffTable::NPARAMS_PER_LINE);
+    ASSERT_EQ(utils::count_words(line), PairTersoffTable::NPARAMS_PER_LINE);
+}
+
 TEST_F(PotentialFileReaderTest, TersoffZBL)
 {
     if (!verbose) ::testing::internal::CaptureStdout();