Merge pull request #10 from monk-04/lammps_implementation

Syncing to lammps_implementation

examples/tungsten_extxyz/model_2.json

@@ -1,7 +1,7 @@
 {
     "coefficients": {
-        "W": -7.7318674079835654,
-        "W-W": [68.784447407502043, 68.784447407502043, 48.652658262777322, 32.029843994880935, 18.916004603812866, 9.1783267118526144, 2.5038442647278982, 0.19953086430922376, -0.24027843128574466, -0.30607305653387484, -0.2169403990617684, 0.0020888316560132309, 0.011224689615986055, -0.0045318560134850359, 0.019103809780459731, 0, 0, 0]
+        "W": -9.3238674348400146,
+        "W-W": [72.09603726779028, 72.09603726779028, 50.837630561204165, 33.330632287281254, 19.575042446021545, 9.1648391313336894, 2.6141295842339058, 0.30838268122808682, -0.13315123440110913, -0.20440641753947206, -0.13045211036857915, 0.062789656733124477, 0.045162501429063438, 0.016332261041836296, 0.027809874694213737, 0, 0, 0]
     },
     "knots": {
         "W-W": [0.001, 0.001, 0.001, 0.001, 0.36759999999999998, 0.73419999999999996, 1.1007999999999998, 1.4673999999999998, 1.8339999999999999, 2.2005999999999997, 2.5671999999999997, 2.9337999999999997, 3.3003999999999998, 3.6669999999999998, 4.0335999999999999, 4.4001999999999999, 4.7667999999999999, 5.1334, 5.5, 5.5, 5.5, 5.5]

lammps_plugin/ML-UF3/pair_uf3.cpp

@@ -106,14 +106,20 @@ void PairUF3::coeff(int narg, char **arg)
 {
   if (!allocated) allocate();
 
-  if (narg != 3 && narg != 4)
-     error->warning(FLERR, "\nUF3: WARNING!! It seems that you are using the \n\
+  if (narg != 3 && narg != 4){
+     /*error->warning(FLERR, "\nUF3: WARNING!! It seems that you are using the \n\
              older style of specifying UF3 POT files. This style of listing \n\
              all the potential files on a single line will be depcrecated in \n\
-             the next version of ML-UF3");
-    //error->all(FLERR, "UF3: Invalid number of argument in pair coeff; \n\
-    //        Provide the species number followed by the LAMMPS POT file");
-
+             the next version of ML-UF3");*/
+    if (narg == tot_pot_files + 2)
+      error->all(FLERR, "UF3 The old style of listing all the potential files\n\
+              on a single line is depcrecated");
+    else
+      error->all(FLERR, "UF3: Invalid number of argument in pair coeff;\n\
+              Provide the species number followed by the LAMMPS POT file\n\
+              Eg. 'pair_coeff 1 1 POT_FILE' for 2-body and \n\
+              'pair_coeff 1 2 2 POT_FILE' for 3-body.");
+  }
   if (narg == 3 || narg == 4){
     int ilo, ihi, jlo, jhi, klo, khi;
     utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error);
@@ -145,7 +151,7 @@ void PairUF3::coeff(int narg, char **arg)
     }
   }
 
-  else{
+  /*else{
     if (narg != tot_pot_files + 2)
       error->all(FLERR,"UF3: Invalid number of argument in pair coeff; \n\
               Number of potential files provided is not correct");
@@ -158,7 +164,7 @@ void PairUF3::coeff(int narg, char **arg)
     // open UF3 potential file on all proc
     for (int i = 2; i < narg; i++) { uf3_read_pot_file(arg[i]); }
     if (!bsplines_created) create_bsplines();
-    /*
+    
     // setflag check needed here
     for (int i = 1; i < num_of_elements + 1; i++) {
       for (int j = 1; j < num_of_elements + 1; j++) {
@@ -194,8 +200,8 @@ void PairUF3::coeff(int narg, char **arg)
         }
       }
     }
-  }*/
   }
+  }*/
 }
 
 void PairUF3::allocate()
@@ -212,6 +218,9 @@ void PairUF3::allocate()
   memory->create(cutsq, num_of_elements + 1, num_of_elements + 1, "pair:cutsq");
   // cut is specific to this pair style. We will set the values in cut
   memory->create(cut, num_of_elements + 1, num_of_elements + 1, "pair:cut");
+  //Contains info about type of knot_spacing--> 0 = uniform knot spacing (default)
+  //1 = non-uniform knot spacing
+  memory->create(knot_spacing_type_2b, num_of_elements + 1, num_of_elements + 1, "pair:knot_spacing_2b");
 
   // Contains knot_vect of 2-body potential for type i and j
   n2b_knot.resize(num_of_elements + 1);
@@ -235,6 +244,11 @@ void PairUF3::allocate()
     // Contains info about minimum 3-body cutoff distance for type i, j and k
     memory->create(min_cut_3b, num_of_elements + 1, num_of_elements + 1, num_of_elements + 1, 3,
                     "pair:min_cut_3b");
+  //Contains info about type of knot_spacing--> 0 = uniform knot spacing (default)
+  //1 = non-uniform knot spacing
+  memory->create(knot_spacing_type_3b, num_of_elements + 1, num_of_elements + 1,
+          num_of_elements + 1, "pair:knot_spacing_3b");
+
 
     // setting cut_3b and setflag = 0
     for (int i = 1; i < num_of_elements + 1; i++) {
@@ -287,9 +301,9 @@ void PairUF3::uf3_read_pot_file(int itype, int jtype, char *potf_name)
   temp_line = txtfilereader.next_line(1);
   ValueTokenizer fp2nd_line(temp_line);
 
-  if (fp2nd_line.count() != 3)
-    error->all(FLERR, "UF3: Expected 3 words on 2nd line =>\n\
-            nBody leading_trim trailing_trim\n\
+  if (fp2nd_line.count() != 4)
+    error->all(FLERR, "UF3: Expected 4 words on 2nd line =>\n\
+            nBody leading_trim trailing_trim type_of_knot_spacing\n\
             Found {}",temp_line);
 
   std::string nbody_on_file = fp2nd_line.next_string();
@@ -307,6 +321,25 @@ void PairUF3::uf3_read_pot_file(int itype, int jtype, char *potf_name)
   if (trailing_trim != 3)
     error->all(FLERR, "UF3: Current implementation is throughly tested only for\n\
             trailing_trim=3\n");
+
+  std::string knot_type = fp2nd_line.next_string();
+  if (utils::strmatch(knot_type,"uk")){
+    utils::logmesg(lmp, "UF3: File {} contains 2-body UF3 potential with uniform\n\
+              knot spacing\n",potf_name);
+    knot_spacing_type_2b[itype][jtype] = 0;
+    knot_spacing_type_2b[jtype][itype] = 0;
+  }   
+  else if (utils::strmatch(knot_type,"nk")){
+    utils::logmesg(lmp, "UF3: File {} contains 2-body UF3 potential with non-uniform\n\
+            knot spacing\n",potf_name);
+    knot_spacing_type_2b[itype][jtype] = 1;
+    knot_spacing_type_2b[jtype][itype] = 1;
+    /*error->all(FLERR, "UF3: Current implementation only works with uniform\n\
+            knot spacing");*/
+  }
+  else
+    error->all(FLERR, "UF3: Expected either 'uk'(uniform-knots) or 'nk'(non-uniform knots)\n\
+            Found {} on the 2nd line of {} pot file",knot_type,potf_name);
 
   temp_line = txtfilereader.next_line(1);
   ValueTokenizer fp3rd_line(temp_line);
@@ -387,11 +420,11 @@ void PairUF3::uf3_read_pot_file(int itype, int jtype, int ktype, char *potf_name
   temp_line = txtfilereader.next_line(1);
   ValueTokenizer fp2nd_line(temp_line);
 
-  if (fp2nd_line.count() != 3)
+  if (fp2nd_line.count() != 4)
     error->all(FLERR, "UF3: Expected 3 words on 2nd line =>\n\
-            nBody leading_trim trailing_trim\n\
+            nBody leading_trim trailing_trim type_of_knot_spacing\n\
             Found {}",temp_line);
-   
+
   std::string nbody_on_file = fp2nd_line.next_string();
 
   if (utils::strmatch(nbody_on_file,"3B"))
@@ -408,7 +441,26 @@ void PairUF3::uf3_read_pot_file(int itype, int jtype, int ktype, char *potf_name
   if (trailing_trim != 3)
     error->all(FLERR, "UF3: Current implementation is throughly tested only for\n\
             trailing_trim=3\n");
-
+
+  std::string knot_type = fp2nd_line.next_string();
+  if (utils::strmatch(knot_type,"uk")){
+    utils::logmesg(lmp, "UF3: File {} contains 3-body UF3 potential with uniform\n\
+            knot spacing\n",potf_name);
+    knot_spacing_type_3b[itype][jtype][ktype] = 0;
+    knot_spacing_type_3b[itype][ktype][jtype] = 0;
+  }
+  else if (utils::strmatch(knot_type,"nk")){
+    utils::logmesg(lmp, "UF3: File {} contains 3-body UF3 potential with non-uniform\n\
+            knot spacing\n",potf_name);
+    knot_spacing_type_3b[itype][jtype][ktype] = 1;
+    knot_spacing_type_3b[itype][ktype][jtype] = 1;
+    /*error->all(FLERR, "UF3: Current implementation only works with uniform\n\
+            knot spacing");*/
+  }
+  else
+    error->all(FLERR, "UF3: Expected either 'uk'(uniform-knots) or 'nk'(non-uniform knots)\n\
+            Found {} on the 2nd line of {} pot file",knot_type,potf_name);
+
   temp_line = txtfilereader.next_line(6);
   ValueTokenizer fp3rd_line(temp_line);
 
@@ -848,18 +900,21 @@ void PairUF3::create_bsplines()
 
   for (int i = 1; i < num_of_elements + 1; i++) {
     for (int j = i; j < num_of_elements + 1; j++) {
-      UFBS2b[i][j] = uf3_pair_bspline(lmp, n2b_knot[i][j], n2b_coeff[i][j]);
+      UFBS2b[i][j] = uf3_pair_bspline(lmp, n2b_knot[i][j], n2b_coeff[i][j],
+              knot_spacing_type_2b[i][j]);
       UFBS2b[j][i] = UFBS2b[i][j];
     }
     if (pot_3b) {
       for (int j = 1; j < num_of_elements + 1; j++) {
         for (int k = j; k < num_of_elements + 1; k++) {
           std::string key = std::to_string(i) + std::to_string(j) + std::to_string(k);
           UFBS3b[i][j][k] = 
-              uf3_triplet_bspline(lmp, n3b_knot_matrix[i][j][k], n3b_coeff_matrix[key]);
+              uf3_triplet_bspline(lmp, n3b_knot_matrix[i][j][k], n3b_coeff_matrix[key],
+                      knot_spacing_type_3b[i][j][k]);
           std::string key2 = std::to_string(i) + std::to_string(k) + std::to_string(j);
           UFBS3b[i][k][j] =
-              uf3_triplet_bspline(lmp, n3b_knot_matrix[i][k][j], n3b_coeff_matrix[key2]);
+              uf3_triplet_bspline(lmp, n3b_knot_matrix[i][k][j], n3b_coeff_matrix[key2],
+                      knot_spacing_type_3b[i][k][j]);
         }
       }
     }

lammps_plugin/ML-UF3/pair_uf3.h

@@ -56,7 +56,7 @@ class PairUF3 : public Pair {
   int bsplines_created;
   int coeff_matrix_dim1, coeff_matrix_dim2, coeff_matrix_dim3, coeff_matrix_elements_len;
   bool pot_3b;
-  int ***setflag_3b;
+  int ***setflag_3b, **knot_spacing_type_2b, ***knot_spacing_type_3b;
   double **cut, ***cut_3b, **cut_3b_list, ****min_cut_3b;
   virtual void allocate();
   void create_bsplines();

lammps_plugin/ML-UF3/uf3_pair_bspline.cpp

@@ -4,6 +4,7 @@
 #include "uf3_bspline_basis3.h"
 
 #include "utils.h"
+#include "error.h"
 #include <vector>
 
 using namespace LAMMPS_NS;
@@ -14,15 +15,30 @@ uf3_pair_bspline::uf3_pair_bspline() {}
 // Constructor
 // Passing vectors by reference
 uf3_pair_bspline::uf3_pair_bspline(LAMMPS *ulmp, const std::vector<double> &uknot_vect,
-                                   const std::vector<double> &ucoeff_vect)
+                                   const std::vector<double> &ucoeff_vect, 
+                                   const int &uknot_spacing_type)
 {
   lmp = ulmp;
   knot_vect = uknot_vect;
   coeff_vect = ucoeff_vect;
+
+  knot_spacing_type = uknot_spacing_type;
+  if (knot_spacing_type==0){
+    knot_spacing = knot_vect[4]-knot_vect[3];
+    get_starting_index=&uf3_pair_bspline::get_starting_index_uniform;
+  }
+  else if (knot_spacing_type==1){
+    knot_spacing = 0;
+    get_starting_index=&uf3_pair_bspline::get_starting_index_nonuniform;
+  }
 
+  else
+    lmp->error->all(FLERR, "UF3: Expected either '0'(uniform-knots) or \n\
+            '1'(non-uniform knots)");
+
   knot_vect_size = uknot_vect.size();
   coeff_vect_size = ucoeff_vect.size();
-
+  
   // Initialize B-Spline Basis Functions
   for (int i = 0; i < knot_vect.size() - 4; i++)
     bspline_bases.push_back(uf3_bspline_basis3(lmp, &knot_vect[i], coeff_vect[i]));
@@ -51,21 +67,38 @@ uf3_pair_bspline::uf3_pair_bspline(LAMMPS *ulmp, const std::vector<double> &ukno
 
 uf3_pair_bspline::~uf3_pair_bspline() {}
 
+int uf3_pair_bspline::get_starting_index_uniform(double r)
+{
+  return 3+(int)((r-knot_vect[0])/knot_spacing);
+}
+
+int uf3_pair_bspline::get_starting_index_nonuniform(double r)
+{
+  if (knot_vect.front() <= r && r < knot_vect.back()) {
+    //Determine the interval for value_rij
+    for (int i = 3; i < knot_vect_size - 1; ++i) {
+      if (knot_vect[i] <= r && r < knot_vect[i + 1]) {
+        return i;
+      }
+    }
+  }
+}
+
 double *uf3_pair_bspline::eval(double r)
 {
 
   // Find knot starting position
 
-  int start_index;
-  if (knot_vect.front() <= r && r < knot_vect.back()) {
+  int start_index=(this->*get_starting_index)(r);
+  /*if (knot_vect.front() <= r && r < knot_vect.back()) {
     //Determine the interval for value_rij
     for (int i = 3; i < knot_vect_size - 1; ++i) {
       if (knot_vect[i] <= r && r < knot_vect[i + 1]) {
         start_index = i;
         break;
       }
     }
-  }
+  }*/
 
   int knot_affect_start = start_index - 3;
 

lammps_plugin/ML-UF3/uf3_pair_bspline.h

@@ -33,17 +33,21 @@ class uf3_pair_bspline {
   std::vector<double> coeff_vect, dncoeff_vect;
   std::vector<uf3_bspline_basis3> bspline_bases;
   std::vector<uf3_bspline_basis2> dnbspline_bases;
+  int knot_spacing_type;
+  int get_starting_index_uniform(double), get_starting_index_nonuniform(double);
+  int (uf3_pair_bspline::*get_starting_index)(double);
+  double knot_spacing=0; 
   LAMMPS *lmp;
 
  public:
   // dummy constructor
   uf3_pair_bspline();
   uf3_pair_bspline(LAMMPS *ulmp, const std::vector<double> &uknot_vect,
-                   const std::vector<double> &ucoeff_vect);
+                   const std::vector<double> &ucoeff_vect,
+                   const int &uknot_spacing_type);
   ~uf3_pair_bspline();
   double ret_val[2];
   double *eval(double value_rij);
-
   double memory_usage();
 };
 }    // namespace LAMMPS_NS