Repository: uf3/uf3 (language: python)

Issue:
Two-and-three body support for UFCalculator and UFLammps
- [WIP] Three-body contributions for UFCalculator for multicomponent systems are experiencing some minor bugs
- Three-body contributions for UFLammps will require a tensor-product spline evaluator in C++ for LAMMPS


Discussion / hints:
Note that training and testing with two-and-three body interactions is fully functional for both single-component and multi-component systems.

Fix patch:
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 00000000..c7d5a2b5
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,58 @@
+
+uf3.egg-info/dependency_links.txt
+
+uf3.egg-info/PKG-INFO
+
+uf3.egg-info/requires.txt
+
+uf3.egg-info/SOURCES.txt
+
+uf3.egg-info/top_level.txt
+
+uf3/__pycache__/__init__.cpython-37.pyc
+
+uf3/data/__pycache__/__init__.cpython-37.pyc
+
+uf3/data/__pycache__/composition.cpython-37.pyc
+
+uf3/data/__pycache__/geometry.cpython-37.pyc
+
+uf3/data/__pycache__/io.cpython-37.pyc
+
+uf3/forcefield/__pycache__/__init__.cpython-37.pyc
+
+uf3/forcefield/__pycache__/calculator.cpython-37.pyc
+
+uf3/forcefield/__pycache__/lammps.cpython-37.pyc
+
+uf3/forcefield/properties/__pycache__/elastic.cpython-37.pyc
+
+uf3/forcefield/properties/__pycache__/phonon.cpython-37.pyc
+
+uf3/regression/__pycache__/__init__.cpython-37.pyc
+
+uf3/regression/__pycache__/least_squares.cpython-37.pyc
+
+uf3/regression/__pycache__/regularize.cpython-37.pyc
+
+uf3/representation/__pycache__/__init__.cpython-37.pyc
+
+uf3/representation/__pycache__/angles.cpython-37.pyc
+
+uf3/representation/__pycache__/bspline.cpython-37.pyc
+
+uf3/representation/__pycache__/distances.cpython-37.pyc
+
+uf3/representation/__pycache__/process.cpython-37.pyc
+
+uf3/util/__pycache__/__init__.cpython-37.pyc
+
+uf3/util/__pycache__/json_io.cpython-37.pyc
+
+uf3/util/__pycache__/parallel.cpython-37.pyc
+
+uf3/util/__pycache__/plotting.cpython-37.pyc
+
+uf3/util/__pycache__/subsample.cpython-37.pyc
+
+lammps_plugin/ML-UF3/.vscode/
diff --git a/lammps_plugin/ML-UF3/.clang-format b/lammps_plugin/ML-UF3/.clang-format
new file mode 100644
index 00000000..d5bb8c37
--- /dev/null
+++ b/lammps_plugin/ML-UF3/.clang-format
@@ -0,0 +1,32 @@
+---
+Language:                             Cpp
+BasedOnStyle:                         LLVM
+AccessModifierOffset:                 -1
+AlignConsecutiveAssignments:          false
+AlignConsecutiveDeclarations:         false
+AlignEscapedNewlines:                 Left
+AlignOperands:                        false
+AllowShortBlocksOnASingleLine:        true
+AllowShortFunctionsOnASingleLine:     Inline
+AllowShortIfStatementsOnASingleLine:  WithoutElse
+AllowShortLambdasOnASingleLine:       None
+AllowShortLoopsOnASingleLine:         true
+BraceWrapping:
+  AfterFunction:                      true
+BreakBeforeBraces:                    Custom
+BreakConstructorInitializers:         AfterColon
+BreakInheritanceList:                 AfterColon
+ColumnLimit:                          100
+IndentCaseLabels:                     true
+IndentWidth:                          2
+NamespaceIndentation:                 Inner
+ObjCBlockIndentWidth:                 2
+PenaltyBreakAssignment:               4
+ReflowComments:                       false
+SpaceAfterCStyleCast:                 true
+SpacesBeforeTrailingComments:         4
+SpacesInContainerLiterals:            false
+Standard:                             Cpp11
+TabWidth:                             2
+UseTab:                               Never
+...
diff --git a/lammps_plugin/ML-UF3/pair_uf3.cpp b/lammps_plugin/ML-UF3/pair_uf3.cpp
new file mode 100644
index 00000000..f016d42b
--- /dev/null
+++ b/lammps_plugin/ML-UF3/pair_uf3.cpp
@@ -0,0 +1,1198 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, 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.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ *    Contributing authors: Ajinkya Hire (U of Florida), 
+ *                          Hendrik Kraß (U of Constance),
+ *                          Richard Hennig (U of Florida)
+ * ---------------------------------------------------------------------- */
+
+#include "pair_uf3.h"
+#include "uf3_pair_bspline.h"
+#include "uf3_triplet_bspline.h"
+
+#include "atom.h"
+#include "comm.h"
+#include "error.h"
+#include "force.h"
+#include "math_const.h"
+#include "memory.h"
+#include "neigh_list.h"
+#include "neighbor.h"
+#include "text_file_reader.h"
+
+#include <cmath>
+
+using namespace LAMMPS_NS;
+using namespace MathConst;
+
+PairUF3::PairUF3(LAMMPS *lmp) : Pair(lmp)
+{
+  single_enable = 1;    // 1 if single() routine exists
+  restartinfo = 0;      // 1 if pair style writes restart info
+  maxshort = 10;
+  neighshort = nullptr;
+  centroidstressflag = CENTROID_AVAIL;
+  manybody_flag = 1;
+  one_coeff = 0; //if 1 then allow only one coeff call of form 'pair_coeff * *'
+                 //by setting it to 0 we will allow multiple 'pair_coeff' calls
+  bsplines_created = 0;
+}
+
+PairUF3::~PairUF3()
+{
+  if (copymode) return;
+  if (allocated) {
+    memory->destroy(setflag);
+    memory->destroy(cutsq);
+    memory->destroy(cut);
+
+    if (pot_3b) {
+      memory->destroy(setflag_3b);
+      memory->destroy(cut_3b);
+      memory->destroy(cut_3b_list);
+      memory->destroy(min_cut_3b);
+      memory->destroy(neighshort);
+    }
+  }
+}
+
+/* ----------------------------------------------------------------------
+ *     global settings
+ * ---------------------------------------------------------------------- */
+
+void PairUF3::settings(int narg, char **arg)
+{
+
+  if (narg != 2)
+    error->all(FLERR, "UF3: Invalid number of argument in pair settings\n\
+            Are you running 2-body or 2 & 3-body UF potential\n\
+            Also how many elements?");
+  nbody_flag = utils::numeric(FLERR, arg[0], true, lmp);
+  num_of_elements = utils::numeric(FLERR, arg[1], true, lmp);    // atom->ntypes;
+  if (num_of_elements != atom->ntypes) {
+    if (comm->me == 0)
+      utils::logmesg(lmp, "\nUF3: Number of elements provided in the input file and \
+number of elements detected by lammps in the structure are not same\n\
+     proceed with caution\n");
+  }
+  if (nbody_flag == 2) {
+    pot_3b = false;
+    n2body_pot_files = num_of_elements * (num_of_elements + 1) / 2;
+    tot_pot_files = n2body_pot_files;
+  } else if (nbody_flag == 3) {
+    pot_3b = true;
+    n2body_pot_files = num_of_elements * (num_of_elements + 1) / 2;
+    n3body_pot_files = num_of_elements * (num_of_elements * (num_of_elements + 1) / 2);
+    tot_pot_files = n2body_pot_files + n3body_pot_files;
+  } else
+    error->all(FLERR, "UF3: UF3 not yet implemented for {}-body", nbody_flag);
+}
+
+/* ----------------------------------------------------------------------
+ *    set coeffs for one or more type pairs
+ * ---------------------------------------------------------------------- */
+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\
+             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");
+
+  if (narg == 3 || narg == 4){
+    int ilo, ihi, jlo, jhi, klo, khi;
+    utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error);
+    utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error);
+    if (narg == 4)
+      utils::bounds(FLERR, arg[2], 1, atom->ntypes, klo, khi, error);
+ 
+
+    if (narg == 3){
+      for (int i = ilo; i <= ihi; i++) {
+        for (int j = MAX(jlo, i); j <= jhi; j++) {
+          if (comm->me == 0)
+            utils::logmesg(lmp, "\nUF3: Opening {} file\n", arg[2]);
+          uf3_read_pot_file(i,j,arg[2]);
+        }
+      }
+    }
+
+    if (narg == 4){
+      for (int i = ilo; i <= ihi; i++) {
+        for (int j = jlo; j <= jhi; j++) {
+          for (int k = MAX(klo, jlo); k <= khi; k++) {
+            if (comm->me == 0)
+              utils::logmesg(lmp, "\nUF3: Opening {} file\n", arg[3]);
+            uf3_read_pot_file(i,j,k,arg[3]);
+          }
+        }
+      }
+    }
+  }
+
+  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");
+
+    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");
+
+    // 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++) {
+        if (setflag[i][j] != 1)
+          error->all(FLERR,"UF3: Not all 2-body UF potentials are set, \n\
+                missing potential file for {}-{} interaction",i, j);
+      }
+    }
+
+    if (pot_3b) {
+      for (int i = 1; i < num_of_elements + 1; i++) {
+        for (int j = 1; j < num_of_elements + 1; j++) {
+          for (int k = 1; k < num_of_elements + 1; k++) {
+            if (setflag_3b[i][j][k] != 1)
+              error->all(FLERR,"UF3: Not all 3-body UF potentials are set, \n\
+                    missing potential file for {}-{}-{} interaction", i, j, k);
+          }
+        }
+      }
+    }
+  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[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]);
+          UFBS3b[i][k][j] = UFBS3b[i][j][k];
+        }
+      }
+    }
+  }*/
+  }
+}
+
+void PairUF3::allocate()
+{
+  allocated = 1;
+
+  // Contains info about wether UF potential were found for type i and j
+  memory->create(setflag, num_of_elements + 1, num_of_elements + 1, "pair:setflag");
+
+  // Contains info about 2-body cutoff distance for type i and j
+  // cutsq is the global variable
+  // Even though we are making cutsq don't manually change the default values
+  // Lammps take care of setting the value
+  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 knot_vect of 2-body potential for type i and j
+  n2b_knot.resize(num_of_elements + 1);
+  n2b_coeff.resize(num_of_elements + 1);
+  UFBS2b.resize(num_of_elements + 1);
+  for (int i = 1; i < num_of_elements + 1; i++) {
+    n2b_knot[i].resize(num_of_elements + 1);
+    n2b_coeff[i].resize(num_of_elements + 1);
+    UFBS2b[i].resize(num_of_elements + 1);
+  }
+  if (pot_3b) {
+    // Contains info about wether UF potential were found for type i, j and k
+    memory->create(setflag_3b, num_of_elements + 1, num_of_elements + 1, num_of_elements + 1,
+                   "pair:setflag_3b");
+    // Contains info about 3-body cutoff distance for type i, j and k
+    memory->create(cut_3b, num_of_elements + 1, num_of_elements + 1, num_of_elements + 1,
+                   "pair:cut_3b");
+    // Contains info about 3-body cutoff distance for type i, j and k
+    // for constructing 3-body list
+    memory->create(cut_3b_list, num_of_elements + 1, num_of_elements + 1, "pair:cut_3b_list");
+    // 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");
+
+    // setting cut_3b and setflag = 0
+    for (int i = 1; i < num_of_elements + 1; i++) {
+      for (int j = 1; j < num_of_elements + 1; j++) {
+        cut_3b_list[i][j] = 0;
+        for (int k = 1; k < num_of_elements + 1; k++) {
+          cut_3b[i][j][k] = 0;
+          min_cut_3b[i][j][k][0] = 0;
+          min_cut_3b[i][j][k][1] = 0;
+          min_cut_3b[i][j][k][2] = 0;
+        }
+      }
+    }
+    n3b_knot_matrix.resize(num_of_elements + 1);
+    UFBS3b.resize(num_of_elements + 1);
+    for (int i = 1; i < num_of_elements + 1; i++) {
+      n3b_knot_matrix[i].resize(num_of_elements + 1);
+      UFBS3b[i].resize(num_of_elements + 1);
+      for (int j = 1; j < num_of_elements + 1; j++) {
+        n3b_knot_matrix[i][j].resize(num_of_elements + 1);
+        UFBS3b[i][j].resize(num_of_elements + 1);
+      }
+    }
+    memory->create(neighshort, maxshort, "pair:neighshort");
+  }
+}
+
+void PairUF3::uf3_read_pot_file(int itype, int jtype, char *potf_name)
+{
+  utils::logmesg(lmp, "UF3: {} file should contain UF3 potential for {} {}\n", \
+          potf_name, itype, jtype);
+
+  FILE *fp;
+  fp = utils::open_potential(potf_name, lmp, nullptr);
+
+  TextFileReader txtfilereader(fp, "UF3:POTFP");
+  txtfilereader.ignore_comments = false;
+
+  std::string temp_line = txtfilereader.next_line(1);
+  Tokenizer file_header(temp_line);
+  
+  if (file_header.count() != 2)
+    error->all(FLERR, "UF3: Expected only two words on 1st line of {} but found \n\
+            {} word/s",potf_name,file_header.count());
+
+  if (file_header.contains("#UF3 POT") == 0)
+    error->all(FLERR, "UF3: {} file is not UF3 POT type, 1st line of UF3 POT \n\
+            files contain '#UF3 POT'. Found {} in the header",potf_name,temp_line);
+
+  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\
+            Found {}",temp_line);
+   
+  std::string nbody_on_file = fp2nd_line.next_string();
+  if (utils::strmatch(nbody_on_file,"2B"))
+    utils::logmesg(lmp, "UF3: File {} contains 2-body UF3 potential\n",potf_name);
+  else
+    error->all(FLERR, "UF3: Expected a 2B UF3 file but found {}",
+            nbody_on_file);
+
+  int leading_trim = fp2nd_line.next_int();
+  int trailing_trim = fp2nd_line.next_int();
+  if (leading_trim != 0)
+    error->all(FLERR, "UF3: Current implementation is throughly tested only for\n\
+            leading_trim=0\n");
+  if (trailing_trim != 3)
+    error->all(FLERR, "UF3: Current implementation is throughly tested only for\n\
+            trailing_trim=3\n");
+
+  temp_line = txtfilereader.next_line(1);
+  ValueTokenizer fp3rd_line(temp_line);
+  if (fp3rd_line.count() != 2)
+    error->all(FLERR, "UF3: Expected only 2 numbers on 3rd line =>\n\
+            Rij_CUTOFF NUM_OF_KNOTS\n\
+            Found {} number/s",fp3rd_line.count());
+  
+  //cut is used in init_one which is called by pair.cpp at line 267 where the return of init_one is squared
+  cut[itype][jtype] = fp3rd_line.next_double();
+  cut[jtype][itype] = cut[itype][jtype];
+
+  int num_knots_2b = fp3rd_line.next_int();
+
+  temp_line = txtfilereader.next_line(num_knots_2b);
+  ValueTokenizer fp4th_line(temp_line);
+  
+  if (fp4th_line.count() != num_knots_2b)
+      error->all(FLERR, "UF3: Expected {} numbers on 4th line but found {} numbers",
+              num_knots_2b,fp4th_line.count());
+
+  n2b_knot[itype][jtype].resize(num_knots_2b);
+  n2b_knot[jtype][itype].resize(num_knots_2b);
+  for (int k = 0; k < num_knots_2b; k++) {
+    n2b_knot[itype][jtype][k] = fp4th_line.next_double();
+    n2b_knot[jtype][itype][k] = n2b_knot[itype][jtype][k];
+  }
+
+  temp_line = txtfilereader.next_line(1);
+  ValueTokenizer fp5th_line(temp_line);
+  int num_of_coeff_2b = fp5th_line.next_int();
+
+  temp_line = txtfilereader.next_line(num_of_coeff_2b);
+  ValueTokenizer fp6th_line(temp_line);
+
+  if (fp6th_line.count() != num_of_coeff_2b)
+    error->all(FLERR, "UF3: Expected {} numbers on 6th line but found {} numbers",
+            num_of_coeff_2b, fp6th_line.count());
+
+  n2b_coeff[itype][jtype].resize(num_of_coeff_2b);
+  n2b_coeff[jtype][itype].resize(num_of_coeff_2b);
+  for (int k = 0; k < num_of_coeff_2b; k++) {
+    n2b_coeff[itype][jtype][k] = fp6th_line.next_double();
+    n2b_coeff[jtype][itype][k] = n2b_coeff[itype][jtype][k];
+  }
+  
+  if (n2b_knot[itype][jtype].size() != n2b_coeff[itype][jtype].size() + 4) {
+    error->all(FLERR, "UF3: {} has incorrect knot and coeff data nknots!=ncoeffs + 3 +1",
+               potf_name);
+  }
+  setflag[itype][jtype] = 1;
+  setflag[jtype][itype] = 1;
+}
+
+
+void PairUF3::uf3_read_pot_file(int itype, int jtype, int ktype, char *potf_name)
+{
+  utils::logmesg(lmp, "UF3: {} file should contain UF3 potential for {} {} {}\n",
+          potf_name, itype, jtype, ktype);
+
+  FILE *fp;
+  fp = utils::open_potential(potf_name, lmp, nullptr);
+
+  TextFileReader txtfilereader(fp, "UF3:POTFP");
+  txtfilereader.ignore_comments = false;
+
+  std::string temp_line = txtfilereader.next_line(1);
+  Tokenizer file_header(temp_line);
+  
+  if (file_header.count() != 2)
+    error->all(FLERR, "UF3: Expected only two words on 1st line of {} but found \n\
+            {} word/s",potf_name,file_header.count());
+
+  if (file_header.contains("#UF3 POT") == 0)
+    error->all(FLERR, "UF3: {} file is not UF3 POT type, 1st line of UF3 POT \n\
+            files contain '#UF3 POT'. Found {} in the header",potf_name,temp_line);
+
+  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\
+            Found {}",temp_line);
+   
+  std::string nbody_on_file = fp2nd_line.next_string();
+
+  if (utils::strmatch(nbody_on_file,"3B"))
+    utils::logmesg(lmp, "UF3: File {} contains 3-body UF3 potential\n",potf_name);
+  else
+    error->all(FLERR, "UF3: Expected a 3B UF3 file but found {}",
+            nbody_on_file);
+
+  int leading_trim = fp2nd_line.next_int();
+  int trailing_trim = fp2nd_line.next_int();
+  if (leading_trim != 0)
+    error->all(FLERR, "UF3: Current implementation is throughly tested only for\n\
+            leading_trim=0\n");
+  if (trailing_trim != 3)
+    error->all(FLERR, "UF3: Current implementation is throughly tested only for\n\
+            trailing_trim=3\n");
+
+  temp_line = txtfilereader.next_line(6);
+  ValueTokenizer fp3rd_line(temp_line);
+
+  if (fp3rd_line.count() != 6)
+    error->all(FLERR, "UF3: Expected only 6 numbers on 3rd line =>\n\
+            Rjk_CUTOFF Rik_CUTOFF Rij_CUTOFF NUM_OF_KNOTS_JK NUM_OF_KNOTS_IK NUM_OF_KNOTS_IJ\n\
+            Found {} number/s",fp3rd_line.count());
+  
+  double cut3b_rjk = fp3rd_line.next_double();
+  double cut3b_rij = fp3rd_line.next_double();
+  double cut3b_rik = fp3rd_line.next_double();
+
+  if (cut3b_rij != cut3b_rik) {
+     error->all(FLERR, "UF3: rij!=rik, Current implementation only works for rij=rik");
+  }
+
+  if (2 * cut3b_rik != cut3b_rjk) {
+    error->all(FLERR, "UF3: 2rij=2rik!=rik, Current implementation only works \n\
+            for 2rij=2rik!=rik");
+  }
+  
+  cut_3b_list[itype][jtype] = std::max(cut3b_rij, cut_3b_list[itype][jtype]);
+  cut_3b_list[itype][ktype] = std::max(cut_3b_list[itype][ktype], cut3b_rik);
+  
+  cut_3b[itype][jtype][ktype] = cut3b_rij;
+  cut_3b[itype][ktype][jtype] = cut3b_rik;
+
+  int num_knots_3b_jk = fp3rd_line.next_int();
+  temp_line = txtfilereader.next_line(num_knots_3b_jk);
+  ValueTokenizer fp4th_line(temp_line);
+
+  if (fp4th_line.count() != num_knots_3b_jk)
+    error->all(FLERR, "UF3: Expected {} numbers on 4th line but found {} numbers",
+            num_knots_3b_jk, fp4th_line.count());
+ 
+
+  n3b_knot_matrix[itype][jtype][ktype].resize(3);
+  n3b_knot_matrix[itype][ktype][jtype].resize(3);
+
+  n3b_knot_matrix[itype][jtype][ktype][0].resize(num_knots_3b_jk);
+  n3b_knot_matrix[itype][ktype][jtype][0].resize(num_knots_3b_jk);
+
+  for (int i = 0; i < num_knots_3b_jk; i++) {
+    n3b_knot_matrix[itype][jtype][ktype][0][i] = fp4th_line.next_double();
+    n3b_knot_matrix[itype][ktype][jtype][0][i] =
+        n3b_knot_matrix[itype][jtype][ktype][0][i];
+  }
+
+  min_cut_3b[itype][jtype][ktype][0] = n3b_knot_matrix[itype][jtype][ktype][0][0];
+  min_cut_3b[itype][ktype][jtype][0] = min_cut_3b[itype][ktype][jtype][0];
+  if (comm->me == 0)
+      utils::logmesg(lmp, "UF3: 3b min cutoff 0 {} {}\n", potf_name,
+                     min_cut_3b[itype][jtype][ktype][0], 
+                     min_cut_3b[itype][ktype][jtype][0]);
+
+  int num_knots_3b_ik = fp3rd_line.next_int();
+  temp_line = txtfilereader.next_line(num_knots_3b_ik);
+  ValueTokenizer fp5th_line(temp_line);
+
+  if (fp5th_line.count() != num_knots_3b_ik)
+    error->all(FLERR, "UF3: Expected {} numbers on 5th line but found {} numbers",
+            num_knots_3b_ik, fp5th_line.count());
+
+  n3b_knot_matrix[itype][jtype][ktype][1].resize(num_knots_3b_ik);
+  n3b_knot_matrix[itype][ktype][jtype][1].resize(num_knots_3b_ik);
+  for (int i = 0; i < num_knots_3b_ik; i++) {
+      n3b_knot_matrix[itype][jtype][ktype][1][i] = fp5th_line.next_double();
+      n3b_knot_matrix[itype][ktype][jtype][1][i] =
+          n3b_knot_matrix[itype][jtype][ktype][1][i];
+    }
+
+  min_cut_3b[itype][jtype][ktype][1] = n3b_knot_matrix[itype][jtype][ktype][1][0];
+  min_cut_3b[itype][ktype][jtype][1] = min_cut_3b[itype][ktype][jtype][1];
+  if (comm->me == 0)
+    utils::logmesg(lmp, "UF3: 3b min cutoff 1 {} {}\n", potf_name,
+            min_cut_3b[itype][jtype][ktype][1], min_cut_3b[itype][ktype][jtype][1]);
+
+  int num_knots_3b_ij = fp3rd_line.next_int();
+  temp_line = txtfilereader.next_line(num_knots_3b_ij);
+  ValueTokenizer fp6th_line(temp_line);
+
+  if (fp6th_line.count() != num_knots_3b_ij)
+    error->all(FLERR, "UF3: Expected {} numbers on 6th line but found {} numbers",
+            num_knots_3b_ij, fp5th_line.count());
+
+  n3b_knot_matrix[itype][jtype][ktype][2].resize(num_knots_3b_ij);
+  n3b_knot_matrix[itype][ktype][jtype][2].resize(num_knots_3b_ij);
+  for (int i = 0; i < num_knots_3b_ij; i++) {
+    n3b_knot_matrix[itype][jtype][ktype][2][i] = fp6th_line.next_double();
+    n3b_knot_matrix[itype][ktype][jtype][2][i] = 
+        n3b_knot_matrix[itype][jtype][ktype][2][i];
+    }
+
+  min_cut_3b[itype][jtype][ktype][2] = n3b_knot_matrix[itype][jtype][ktype][2][0];
+  min_cut_3b[itype][ktype][jtype][2] = min_cut_3b[itype][ktype][jtype][2];
+  if (comm->me == 0)
+    utils::logmesg(lmp, "UF3: 3b min cutoff 2 {} {}\n", potf_name,
+            min_cut_3b[itype][jtype][ktype][2], min_cut_3b[itype][ktype][jtype][2]);
+
+  temp_line = txtfilereader.next_line(3);
+  ValueTokenizer fp7th_line(temp_line);
+
+  if (fp7th_line.count() != 3)
+    error->all(FLERR, "UF3: Expected 3 numbers on 7th line =>\n\
+           SHAPE_OF_COEFF_MATRIX[I][J][K] \n\
+           found {} numbers", fp7th_line.count());
+
+  coeff_matrix_dim1 = fp7th_line.next_int();
+  coeff_matrix_dim2 = fp7th_line.next_int();
+  coeff_matrix_dim3 = fp7th_line.next_int();
+
+  if (n3b_knot_matrix[itype][jtype][ktype][0].size() != coeff_matrix_dim3 + 3 + 1)
+    error->all(FLERR, "UF3: {} has incorrect knot (NUM_OF_KNOTS_JK) and \n\
+            coeff (coeff_matrix_dim3) data nknots!=ncoeffs + 3 +1", potf_name);
+
+  if (n3b_knot_matrix[itype][jtype][ktype][1].size() != coeff_matrix_dim2 + 3 + 1) 
+    error->all(FLERR, "UF3: {} has incorrect knot (NUM_OF_KNOTS_IK) and \n\
+            coeff (coeff_matrix_dim2) data nknots!=ncoeffs + 3 +1",potf_name);
+
+  if (n3b_knot_matrix[itype][jtype][ktype][2].size() != coeff_matrix_dim1 + 3 + 1) 
+    error->all(FLERR, "UF3: {} has incorrect knot (NUM_OF_KNOTS_IJ) and \n\
+            coeff ()coeff_matrix_dim1 data nknots!=ncoeffs + 3 +1",potf_name);
+
+  coeff_matrix_elements_len = coeff_matrix_dim3;
+  
+  std::string key = std::to_string(itype) + std::to_string(jtype) + std::to_string(ktype);
+  n3b_coeff_matrix[key].resize(coeff_matrix_dim1);
+
+  int line_count = 0;
+  for (int i = 0; i < coeff_matrix_dim1; i++) {
+    n3b_coeff_matrix[key][i].resize(coeff_matrix_dim2);
+    for (int j = 0; j < coeff_matrix_dim2; j++) {
+      temp_line = txtfilereader.next_line(coeff_matrix_elements_len);
+      ValueTokenizer coeff_line(temp_line);
+      n3b_coeff_matrix[key][i][j].resize(coeff_matrix_dim3);
+
+      if (coeff_line.count() != coeff_matrix_elements_len)
+        error->all(FLERR, "UF3: Expected {} numbers on {}th line but found \n\
+                {} numbers",coeff_matrix_elements_len, line_count+8, coeff_line.count());
+      for (int k = 0; k < coeff_matrix_dim3; k++) {
+        n3b_coeff_matrix[key][i][j][k] = coeff_line.next_double();
+      }
+      line_count += 1;
+    }
+  }
+
+  key = std::to_string(itype) + std::to_string(ktype) + std::to_string(jtype);
+  n3b_coeff_matrix[key] = n3b_coeff_matrix[std::to_string(itype) +
+      std::to_string(jtype) + std::to_string(ktype)];
+
+  setflag_3b[itype][jtype][ktype] = 1;
+  setflag_3b[itype][ktype][jtype] = 1;
+
+}
+
+void PairUF3::uf3_read_pot_file(char *potf_name)
+{
+  if (comm->me == 0) utils::logmesg(lmp, "\nUF3: Opening {} file\n", potf_name);
+
+  FILE *fp;
+  fp = utils::open_potential(potf_name, lmp, nullptr);
+  // if (fp) error->all(FLERR,"UF3: {} file not found",potf_name);
+
+  TextFileReader txtfilereader(fp, "UF3:POTFP");
+  txtfilereader.ignore_comments = false;
+
+  std::string temp_line = txtfilereader.next_line(2);
+  Tokenizer fp1st_line(temp_line);
+
+  if (fp1st_line.contains("#UF3 POT") == 0)
+    error->all(FLERR, "UF3: {} file is not UF3 POT type, found type {} {} on the file", potf_name,
+               fp1st_line.next(), fp1st_line.next());
+
+  if (comm->me == 0)
+    utils::logmesg(lmp, "UF3: {} file is of type {} {}\n", potf_name, fp1st_line.next(),
+                   fp1st_line.next());
+
+  temp_line = txtfilereader.next_line(1);
+  Tokenizer fp2nd_line(temp_line);
+  if (fp2nd_line.contains("2B") == 1) {
+    temp_line = txtfilereader.next_line(4);
+    ValueTokenizer fp3rd_line(temp_line);
+    int temp_type1 = fp3rd_line.next_int();
+    int temp_type2 = fp3rd_line.next_int();
+    if (comm->me == 0)
+      utils::logmesg(lmp, "UF3: {} file contains 2-body UF3 potential for {} {}\n", potf_name,
+                     temp_type1, temp_type2);
+
+    //cut is used in init_one which is called by pair.cpp at line 267 where the return of init_one is squared
+    cut[temp_type1][temp_type2] = fp3rd_line.next_double();
+    // if(comm->me==0) utils::logmesg(lmp,"UF3: Cutoff {}\n",cutsq[temp_type1][temp_type2]);
+    cut[temp_type2][temp_type1] = cut[temp_type1][temp_type2];
+
+    int temp_line_len = fp3rd_line.next_int();
+
+    temp_line = txtfilereader.next_line(temp_line_len);
+    ValueTokenizer fp4th_line(temp_line);
+
+    n2b_knot[temp_type1][temp_type2].resize(temp_line_len);
+    n2b_knot[temp_type2][temp_type1].resize(temp_line_len);
+    for (int k = 0; k < temp_line_len; k++) {
+      n2b_knot[temp_type1][temp_type2][k] = fp4th_line.next_double();
+      n2b_knot[temp_type2][temp_type1][k] = n2b_knot[temp_type1][temp_type2][k];
+    }
+
+    temp_line = txtfilereader.next_line(1);
+    ValueTokenizer fp5th_line(temp_line);
+
+    temp_line_len = fp5th_line.next_int();
+
+    temp_line = txtfilereader.next_line(temp_line_len);
+    // utils::logmesg(lmp,"UF3:11 {}",temp_line);
+    ValueTokenizer fp6th_line(temp_line);
+    // if(comm->me==0) utils::logmesg(lmp,"UF3: {}\n",temp_line_len);
+    n2b_coeff[temp_type1][temp_type2].resize(temp_line_len);
+    n2b_coeff[temp_type2][temp_type1].resize(temp_line_len);
+
+    for (int k = 0; k < temp_line_len; k++) {
+      n2b_coeff[temp_type1][temp_type2][k] = fp6th_line.next_double();
+      n2b_coeff[temp_type2][temp_type1][k] = n2b_coeff[temp_type1][temp_type2][k];
+      // if(comm->me==0) utils::logmesg(lmp,"UF3: {}\n",n2b_coeff[temp_type1][temp_type2][k]);
+    }
+    // for(int i=0;i<n2b_coeff[temp_type1][temp_type2].size();i++) if(comm->me==0) utils::logmesg(lmp,"UF3: {}\n",n2b_coeff[temp_type1][temp_type2][i]);
+    if (n2b_knot[temp_type1][temp_type2].size() != n2b_coeff[temp_type1][temp_type2].size() + 4) {
+      error->all(FLERR, "UF3: {} has incorrect knot and coeff data nknots!=ncoeffs + 3 +1",
+                 potf_name);
+    }
+    setflag[temp_type1][temp_type2] = 1;
+    setflag[temp_type2][temp_type1] = 1;
+  } else if (fp2nd_line.contains("3B") == 1) {
+    temp_line = txtfilereader.next_line(9);
+    ValueTokenizer fp3rd_line(temp_line);
+    int temp_type1 = fp3rd_line.next_int();
+    int temp_type2 = fp3rd_line.next_int();
+    int temp_type3 = fp3rd_line.next_int();
+    if (comm->me == 0)
+      utils::logmesg(lmp, "UF3: {} file contains 3-body UF3 potential for {} {} {}\n", potf_name,
+                     temp_type1, temp_type2, temp_type3);
+
+    double cut3b_rjk = fp3rd_line.next_double();
+    double cut3b_rij = fp3rd_line.next_double();
+    // cut_3b[temp_type1][temp_type2] = std::max(cut3b_rij,
+    // cut_3b[temp_type1][temp_type2]);
+    cut_3b_list[temp_type1][temp_type2] = std::max(cut3b_rij, cut_3b_list[temp_type1][temp_type2]);
+    double cut3b_rik = fp3rd_line.next_double();
+    if (cut3b_rij != cut3b_rik) {
+      error->all(FLERR, "UF3: rij!=rik, Current implementation only works for rij=rik");
+    }
+    if (2 * cut3b_rik != cut3b_rjk) {
+      error->all(FLERR,
+                 "UF3: 2rij=2rik!=rik, Current implementation only works for 2rij=2rik!=rik");
+    }
+    // cut_3b[temp_type1][temp_type3] = std::max(cut_3b[temp_type1][temp_type3],cut3b_rik);
+    cut_3b_list[temp_type1][temp_type3] = std::max(cut_3b_list[temp_type1][temp_type3], cut3b_rik);
+    cut_3b[temp_type1][temp_type2][temp_type3] = cut3b_rij;
+    cut_3b[temp_type1][temp_type3][temp_type2] = cut3b_rik;
+
+    int temp_line_len = fp3rd_line.next_int();
+    temp_line = txtfilereader.next_line(temp_line_len);
+    ValueTokenizer fp4th_line(temp_line);
+
+    n3b_knot_matrix[temp_type1][temp_type2][temp_type3].resize(3);
+    n3b_knot_matrix[temp_type1][temp_type3][temp_type2].resize(3);
+
+    n3b_knot_matrix[temp_type1][temp_type2][temp_type3][0].resize(temp_line_len);
+    n3b_knot_matrix[temp_type1][temp_type3][temp_type2][0].resize(temp_line_len);
+    for (int i = 0; i < temp_line_len; i++) {
+      n3b_knot_matrix[temp_type1][temp_type2][temp_type3][0][i] = fp4th_line.next_double();
+      n3b_knot_matrix[temp_type1][temp_type3][temp_type2][0][i] =
+          n3b_knot_matrix[temp_type1][temp_type2][temp_type3][0][i];
+    }
+
+    min_cut_3b[temp_type1][temp_type2][temp_type3][0] = n3b_knot_matrix[temp_type1][temp_type2][temp_type3][0][0];
+    min_cut_3b[temp_type1][temp_type3][temp_type2][0] = min_cut_3b[temp_type1][temp_type3][temp_type2][0];
+    if (comm->me == 0)
+      utils::logmesg(lmp, "UF3: 3b min cutoff 0 {} {}\n", potf_name,
+                     min_cut_3b[temp_type1][temp_type2][temp_type3][0], min_cut_3b[temp_type1][temp_type3][temp_type2][0]);
+
+    temp_line_len = fp3rd_line.next_int();
+    temp_line = txtfilereader.next_line(temp_line_len);
+    ValueTokenizer fp5th_line(temp_line);
+    n3b_knot_matrix[temp_type1][temp_type2][temp_type3][1].resize(temp_line_len);
+    n3b_knot_matrix[temp_type1][temp_type3][temp_type2][1].resize(temp_line_len);
+    for (int i = 0; i < temp_line_len; i++) {
+      n3b_knot_matrix[temp_type1][temp_type2][temp_type3][1][i] = fp5th_line.next_double();
+      n3b_knot_matrix[temp_type1][temp_type3][temp_type2][1][i] =
+          n3b_knot_matrix[temp_type1][temp_type2][temp_type3][1][i];
+    }
+
+    min_cut_3b[temp_type1][temp_type2][temp_type3][1] = n3b_knot_matrix[temp_type1][temp_type2][temp_type3][1][0];
+    min_cut_3b[temp_type1][temp_type3][temp_type2][1] = min_cut_3b[temp_type1][temp_type3][temp_type2][1];
+    if (comm->me == 0)
+      utils::logmesg(lmp, "UF3: 3b min cutoff 1 {} {}\n", potf_name,
+                     min_cut_3b[temp_type1][temp_type2][temp_type3][1], min_cut_3b[temp_type1][temp_type3][temp_type2][1]);
+
+    temp_line_len = fp3rd_line.next_int();
+    temp_line = txtfilereader.next_line(temp_line_len);
+    ValueTokenizer fp6th_line(temp_line);
+    n3b_knot_matrix[temp_type1][temp_type2][temp_type3][2].resize(temp_line_len);
+    n3b_knot_matrix[temp_type1][temp_type3][temp_type2][2].resize(temp_line_len);
+    for (int i = 0; i < temp_line_len; i++) {
+      n3b_knot_matrix[temp_type1][temp_type2][temp_type3][2][i] = fp6th_line.next_double();
+      n3b_knot_matrix[temp_type1][temp_type3][temp_type2][2][i] =
+          n3b_knot_matrix[temp_type1][temp_type2][temp_type3][2][i];
+    }
+
+    min_cut_3b[temp_type1][temp_type2][temp_type3][2] = n3b_knot_matrix[temp_type1][temp_type2][temp_type3][2][0];
+    min_cut_3b[temp_type1][temp_type3][temp_type2][2] = min_cut_3b[temp_type1][temp_type3][temp_type2][2];
+    if (comm->me == 0)
+      utils::logmesg(lmp, "UF3: 3b min cutoff 2 {} {}\n", potf_name,
+                     min_cut_3b[temp_type1][temp_type2][temp_type3][2], min_cut_3b[temp_type1][temp_type3][temp_type2][2]);
+
+    temp_line = txtfilereader.next_line(3);
+    ValueTokenizer fp7th_line(temp_line);
+
+    coeff_matrix_dim1 = fp7th_line.next_int();
+    coeff_matrix_dim2 = fp7th_line.next_int();
+    coeff_matrix_dim3 = fp7th_line.next_int();
+    if (n3b_knot_matrix[temp_type1][temp_type2][temp_type3][0].size() !=
+        coeff_matrix_dim3 + 3 + 1) {
+      error->all(FLERR, "UF3: {} has incorrect knot and coeff data nknots!=ncoeffs + 3 +1",
+                 potf_name);
+    }
+    if (n3b_knot_matrix[temp_type1][temp_type2][temp_type3][1].size() !=
+        coeff_matrix_dim2 + 3 + 1) {
+      error->all(FLERR, "UF3: {} has incorrect knot and coeff data nknots!=ncoeffs + 3 +1",
+                 potf_name);
+    }
+    if (n3b_knot_matrix[temp_type1][temp_type2][temp_type3][2].size() !=
+        coeff_matrix_dim1 + 3 + 1) {
+      error->all(FLERR, "UF3: {} has incorrect knot and coeff data nknots!=ncoeffs + 3 +1",
+                 potf_name);
+    }
+
+    coeff_matrix_elements_len = coeff_matrix_dim3;
+
+    std::string key =
+        std::to_string(temp_type1) + std::to_string(temp_type2) + std::to_string(temp_type3);
+    n3b_coeff_matrix[key].resize(coeff_matrix_dim1);
+    for (int i = 0; i < coeff_matrix_dim1; i++) {
+      n3b_coeff_matrix[key][i].resize(coeff_matrix_dim2);
+      for (int j = 0; j < coeff_matrix_dim2; j++) {
+        temp_line = txtfilereader.next_line(coeff_matrix_elements_len);
+        ValueTokenizer coeff_line(temp_line);
+        n3b_coeff_matrix[key][i][j].resize(coeff_matrix_dim3);
+        for (int k = 0; k < coeff_matrix_dim3; k++) {
+          n3b_coeff_matrix[key][i][j][k] = coeff_line.next_double();
+        }
+      }
+    }
+
+    key = std::to_string(temp_type1) + std::to_string(temp_type3) + std::to_string(temp_type2);
+    n3b_coeff_matrix[key] =
+        n3b_coeff_matrix[std::to_string(temp_type1) + std::to_string(temp_type2) +
+                         std::to_string(temp_type3)];
+    setflag_3b[temp_type1][temp_type2][temp_type3] = 1;
+    setflag_3b[temp_type1][temp_type3][temp_type2] = 1;
+  } else
+    error->all(
+        FLERR,
+        "UF3: {} file does not contain right words indicating whether it is 2 or 3 body potential",
+        potf_name);
+}
+
+/* ----------------------------------------------------------------------
+   init specific to this pair style
+------------------------------------------------------------------------- */
+void PairUF3::init_style()
+{
+  if (force->newton_pair == 0) error->all(FLERR, "UF3: Pair style requires newton pair on");
+  // request a default neighbor list
+  neighbor->add_request(this, NeighConst::REQ_FULL);
+}
+
+/* ----------------------------------------------------------------------
+   init list sets the pointer to full neighbour list requested in previous function
+------------------------------------------------------------------------- */
+
+void PairUF3::init_list(int /*id*/, class NeighList *ptr)
+{
+  list = ptr;
+}
+
+/* ----------------------------------------------------------------------
+   init for one type pair i,j and corresponding j,i
+------------------------------------------------------------------------- */
+double PairUF3::init_one(int i /*i*/, int /*j*/ j)
+{
+
+  if (!bsplines_created) create_bsplines();
+
+  //init_one is called by pair.cpp at line 267 where it is squred
+  //at line 268
+  return cut[i][j];
+}
+
+void PairUF3::create_bsplines()
+{
+  bsplines_created = 1;
+  for (int i = 1; i < num_of_elements + 1; i++) {
+    for (int j = 1; j < num_of_elements + 1; j++) {
+      if (setflag[i][j] != 1)
+        error->all(FLERR,"UF3: Not all 2-body UF potentials are set, \n\
+                missing potential file for {}-{} interaction",i, j);
+    }
+  }
+  if (pot_3b) {
+    for (int i = 1; i < num_of_elements + 1; i++) {
+      for (int j = 1; j < num_of_elements + 1; j++) {
+        for (int k = 1; k < num_of_elements + 1; k++) {
+          if (setflag_3b[i][j][k] != 1)
+            error->all(FLERR,"UF3: Not all 3-body UF potentials are set, \n\
+                    missing potential file for {}-{}-{} interaction", i, j, k);
+        }
+      }
+    }
+  }
+
+  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[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]);
+          UFBS3b[i][k][j] = UFBS3b[i][j][k];
+        }
+      }
+    }
+  }
+}
+
+void PairUF3::compute(int eflag, int vflag)
+{
+  int i, j, k, ii, jj, kk, inum, jnum, itype, jtype, ktype;
+  double xtmp, ytmp, ztmp, delx, dely, delz, evdwl, fpair, fx, fy, fz;
+  double del_rji[3], del_rki[3], del_rkj[3];
+  double fij[3], fik[3], fjk[3];
+  double fji[3], fki[3], fkj[3];
+  double Fi[3], Fj[3], Fk[3];
+  double rsq, rij, rik, rjk;
+  int *ilist, *jlist, *numneigh, **firstneigh;
+
+  ev_init(eflag, vflag);
+
+  double **x = atom->x;
+  double **f = atom->f;
+  int *type = atom->type;
+  int nlocal = atom->nlocal;
+  int newton_pair = force->newton_pair;
+
+  inum = list->inum;
+  ilist = list->ilist;
+  numneigh = list->numneigh;
+  firstneigh = list->firstneigh;
+  // loop over neighbors of my atoms
+  for (ii = 0; ii < inum; ii++) {
+    evdwl = 0;
+    i = ilist[ii];
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+    itype = type[i];
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+    int numshort = 0;
+    for (jj = 0; jj < jnum; jj++) {
+      fx = 0;
+      fy = 0;
+      fz = 0;
+      j = jlist[jj];
+      j &= NEIGHMASK;
+
+      delx = xtmp - x[j][0];
+      dely = ytmp - x[j][1];
+      delz = ztmp - x[j][2];
+
+      rsq = delx * delx + dely * dely + delz * delz;
+      jtype = type[j];
+      if (rsq < cutsq[itype][jtype]) {
+        rij = sqrt(rsq);
+
+        if (pot_3b) {
+          if (rij <= cut_3b_list[itype][jtype]) {
+            neighshort[numshort] = j;
+            if (numshort >= maxshort - 1) {
+              maxshort += maxshort / 2;
+              memory->grow(neighshort, maxshort, "pair:neighshort");
+            }
+            numshort = numshort + 1;
+          }
+        }
+
+        double *pair_eval = UFBS2b[itype][jtype].eval(rij);
+
+        fpair = -1 * pair_eval[1] / rij;
+
+        fx = delx * fpair;
+        fy = dely * fpair;
+        fz = delz * fpair;
+
+        f[i][0] += fx;
+        f[i][1] += fy;
+        f[i][2] += fz;
+        f[j][0] -= fx;
+        f[j][1] -= fy;
+        f[j][2] -= fz;
+
+        if (eflag) evdwl = pair_eval[0];
+
+        if (evflag) {
+          ev_tally_xyz(i, j, nlocal, newton_pair, evdwl, 0.0, fx, fy, fz, delx, dely, delz);
+
+          // Centroid Stress
+          if (vflag_either && cvflag_atom) {
+            double v[6];
+
+            v[0] = delx * fx;
+            v[1] = dely * fy;
+            v[2] = delz * fz;
+            v[3] = delx * fy;
+            v[4] = delx * fz;
+            v[5] = dely * fz;
+
+            cvatom[i][0] += 0.5 * v[0];
+            cvatom[i][1] += 0.5 * v[1];
+            cvatom[i][2] += 0.5 * v[2];
+            cvatom[i][3] += 0.5 * v[3];
+            cvatom[i][4] += 0.5 * v[4];
+            cvatom[i][5] += 0.5 * v[5];
+            cvatom[i][6] += 0.5 * v[3];
+            cvatom[i][7] += 0.5 * v[4];
+            cvatom[i][8] += 0.5 * v[5];
+
+            cvatom[j][0] += 0.5 * v[0];
+            cvatom[j][1] += 0.5 * v[1];
+            cvatom[j][2] += 0.5 * v[2];
+            cvatom[j][3] += 0.5 * v[3];
+            cvatom[j][4] += 0.5 * v[4];
+            cvatom[j][5] += 0.5 * v[5];
+            cvatom[j][6] += 0.5 * v[3];
+            cvatom[j][7] += 0.5 * v[4];
+            cvatom[j][8] += 0.5 * v[5];
+          }
+        }
+      }
+    }
+
+    // 3-body interaction
+    // jth atom
+    jnum = numshort - 1;
+    for (jj = 0; jj < jnum; jj++) {
+      fij[0] = fji[0] = 0;
+      fij[1] = fji[1] = 0;
+      fij[2] = fji[2] = 0;
+      j = neighshort[jj];
+      jtype = type[j];
+      del_rji[0] = x[j][0] - xtmp;
+      del_rji[1] = x[j][1] - ytmp;
+      del_rji[2] = x[j][2] - ztmp;
+      rij =
+          sqrt(((del_rji[0] * del_rji[0]) + (del_rji[1] * del_rji[1]) + (del_rji[2] * del_rji[2])));
+
+      // kth atom
+      for (kk = jj + 1; kk < numshort; kk++) {
+
+        fik[0] = fki[0] = 0;
+        fik[1] = fki[1] = 0;
+        fik[2] = fki[2] = 0;
+
+        fjk[0] = fkj[0] = 0;
+        fjk[1] = fkj[1] = 0;
+        fjk[2] = fkj[2] = 0;
+
+        k = neighshort[kk];
+        ktype = type[k];
+        del_rki[0] = x[k][0] - xtmp;
+        del_rki[1] = x[k][1] - ytmp;
+        del_rki[2] = x[k][2] - ztmp;
+        rik = sqrt(
+            ((del_rki[0] * del_rki[0]) + (del_rki[1] * del_rki[1]) + (del_rki[2] * del_rki[2])));
+
+        if ((rij <= cut_3b[itype][jtype][ktype]) && (rik <= cut_3b[itype][ktype][jtype]) &&
+                (rij >= min_cut_3b[itype][jtype][ktype][0]) && 
+                (rik >= min_cut_3b[itype][jtype][ktype][1])) {
+
+          del_rkj[0] = x[k][0] - x[j][0];
+          del_rkj[1] = x[k][1] - x[j][1];
+          del_rkj[2] = x[k][2] - x[j][2];
+          rjk = sqrt(
+              ((del_rkj[0] * del_rkj[0]) + (del_rkj[1] * del_rkj[1]) + (del_rkj[2] * del_rkj[2])));
+
+          if (rjk >= min_cut_3b[itype][jtype][ktype][2]){
+            //utils::logmesg(lmp, "UF3: {} {} {}",itype,jtype,ktype);
+            double *triangle_eval = UFBS3b[itype][jtype][ktype].eval(rij, rik, rjk);
+
+            fij[0] = *(triangle_eval + 1) * (del_rji[0] / rij);
+            fji[0] = -fij[0];
+            fik[0] = *(triangle_eval + 2) * (del_rki[0] / rik);
+            fki[0] = -fik[0];
+            fjk[0] = *(triangle_eval + 3) * (del_rkj[0] / rjk);
+            fkj[0] = -fjk[0];
+
+            fij[1] = *(triangle_eval + 1) * (del_rji[1] / rij);
+            fji[1] = -fij[1];
+            fik[1] = *(triangle_eval + 2) * (del_rki[1] / rik);
+            fki[1] = -fik[1];
+            fjk[1] = *(triangle_eval + 3) * (del_rkj[1] / rjk);
+            fkj[1] = -fjk[1];
+
+            fij[2] = *(triangle_eval + 1) * (del_rji[2] / rij);
+            fji[2] = -fij[2];
+            fik[2] = *(triangle_eval + 2) * (del_rki[2] / rik);
+            fki[2] = -fik[2];
+            fjk[2] = *(triangle_eval + 3) * (del_rkj[2] / rjk);
+            fkj[2] = -fjk[2];
+
+            Fi[0] = fij[0] + fik[0];
+            Fi[1] = fij[1] + fik[1];
+            Fi[2] = fij[2] + fik[2];
+            f[i][0] += Fi[0];
+            f[i][1] += Fi[1];
+            f[i][2] += Fi[2];
+
+            Fj[0] = fji[0] + fjk[0];
+            Fj[1] = fji[1] + fjk[1];
+            Fj[2] = fji[2] + fjk[2];
+            f[j][0] += Fj[0];
+            f[j][1] += Fj[1];
+            f[j][2] += Fj[2];
+
+            Fk[0] = fki[0] + fkj[0];
+            Fk[1] = fki[1] + fkj[1];
+            Fk[2] = fki[2] + fkj[2];
+            f[k][0] += Fk[0];
+            f[k][1] += Fk[1];
+            f[k][2] += Fk[2];
+
+            if (eflag) evdwl = *triangle_eval;
+
+            if (evflag) { ev_tally3(i, j, k, evdwl, 0, Fj, Fk, del_rji, del_rki);
+              // Centroid stress 3-body term
+              if (vflag_either && cvflag_atom) {
+                double ric[3];
+                ric[0] = THIRD * (-del_rji[0] - del_rki[0]);
+                ric[1] = THIRD * (-del_rji[1] - del_rki[1]);
+                ric[2] = THIRD * (-del_rji[2] - del_rki[2]);
+
+                cvatom[i][0] += ric[0] * Fi[0];
+                cvatom[i][1] += ric[1] * Fi[1];
+                cvatom[i][2] += ric[2] * Fi[2];
+                cvatom[i][3] += ric[0] * Fi[1];
+                cvatom[i][4] += ric[0] * Fi[2];
+                cvatom[i][5] += ric[1] * Fi[2];
+                cvatom[i][6] += ric[1] * Fi[0];
+                cvatom[i][7] += ric[2] * Fi[0];
+                cvatom[i][8] += ric[2] * Fi[1];
+
+                double rjc[3];
+                rjc[0] = THIRD * (del_rji[0] - del_rkj[0]);
+                rjc[1] = THIRD * (del_rji[1] - del_rkj[1]);
+                rjc[2] = THIRD * (del_rji[2] - del_rkj[2]);
+
+                cvatom[j][0] += rjc[0] * Fj[0];
+                cvatom[j][1] += rjc[1] * Fj[1];
+                cvatom[j][2] += rjc[2] * Fj[2];
+                cvatom[j][3] += rjc[0] * Fj[1];
+                cvatom[j][4] += rjc[0] * Fj[2];
+                cvatom[j][5] += rjc[1] * Fj[2];
+                cvatom[j][6] += rjc[1] * Fj[0];
+                cvatom[j][7] += rjc[2] * Fj[0];
+                cvatom[j][8] += rjc[2] * Fj[1];
+
+                double rkc[3];
+                rkc[0] = THIRD * (del_rki[0] + del_rkj[0]);
+                rkc[1] = THIRD * (del_rki[1] + del_rkj[1]);
+                rkc[2] = THIRD * (del_rki[2] + del_rkj[2]);
+
+                cvatom[k][0] += rkc[0] * Fk[0];
+                cvatom[k][1] += rkc[1] * Fk[1];
+                cvatom[k][2] += rkc[2] * Fk[2];
+                cvatom[k][3] += rkc[0] * Fk[1];
+                cvatom[k][4] += rkc[0] * Fk[2];
+                cvatom[k][5] += rkc[1] * Fk[2];
+                cvatom[k][6] += rkc[1] * Fk[0];
+                cvatom[k][7] += rkc[2] * Fk[0];
+                cvatom[k][8] += rkc[2] * Fk[1];
+              }
+            }
+          }
+        }
+      }
+    }
+  }
+  if (vflag_fdotr) virial_fdotr_compute();
+}
+
+double PairUF3::single(int /*i*/, int /*j*/, int itype, int jtype, double rsq,
+                       double /*factor_coul*/, double factor_lj, double &fforce)
+{
+  double value = 0.0;
+  double r = sqrt(rsq);
+
+  if (r < cut[itype][jtype]) {
+    double *pair_eval = UFBS2b[itype][jtype].eval(r);
+    value = pair_eval[0];
+    fforce = factor_lj * pair_eval[1];
+  }
+
+  return factor_lj * value;
+}
+
+double PairUF3::memory_usage()
+{
+  utils::logmesg(lmp, "\nUF3: memory_usage called");
+  double bytes = Pair::memory_usage();
+
+  bytes = 0;
+
+  bytes += (double)5*sizeof(double);     //num_of_elements, nbody_flag, 
+                                        //n2body_pot_files, n3body_pot_files, 
+                                        //tot_pot_files;
+
+  bytes += (double)5*sizeof(double);    //bsplines_created, coeff_matrix_dim1,
+                                        //coeff_matrix_dim2, coeff_matrix_dim3, 
+                                        //coeff_matrix_elements_len
+  bytes += (double)(num_of_elements+1)*(num_of_elements+1)*\
+           (num_of_elements+1)*sizeof(double);      //***setflag_3b
+
+  bytes += (double)(num_of_elements+1)*(num_of_elements+1)*sizeof(double); //cut 
+
+  bytes += (double)(num_of_elements+1)*(num_of_elements+1)*\
+           (num_of_elements+1)*sizeof(double);      //***cut_3b
+
+  bytes += (double)(num_of_elements+1)*(num_of_elements+1)*sizeof(double); //cut_3b_list
+
+  bytes += (double)(num_of_elements+1)*(num_of_elements+1)*\
+           (num_of_elements+1)*3*sizeof(double);    //min_cut_3b
+
+  for (int i=1; i < num_of_elements+1; i++){
+    for (int j=i; j < num_of_elements+1; j++){
+      bytes += (double)2*n2b_knot[i][j].size()*sizeof(double);      //n2b_knot
+      bytes += (double)2*n2b_coeff[i][j].size()*sizeof(double);     //n2b_coeff
+    }
+    if (pot_3b){
+      for (int j = 1; j < num_of_elements + 1; j++) {
+        for (int k = j; k < num_of_elements + 1; k++) {
+          bytes += (double)2*n3b_knot_matrix[i][j][k][0].size()*sizeof(double);
+          bytes += (double)2*n3b_knot_matrix[i][j][k][1].size()*sizeof(double);
+          bytes += (double)2*n3b_knot_matrix[i][j][k][2].size()*sizeof(double);
+
+          std::string key = std::to_string(i) + std::to_string(j) + std::to_string(k);
+
+          for (int l=0; l < n3b_coeff_matrix[key].size(); l++){
+            for (int m=0; m < n3b_coeff_matrix[key][l].size(); m++){
+              bytes += (double)2*n3b_coeff_matrix[key][l][m].size()*sizeof(double);
+              //key = ijk
+              //key = ikj
+            }
+          }
+        }
+      }
+    }
+  }
+
+  for (int i = 1; i < num_of_elements + 1; i++) {
+    for (int j = i; j < num_of_elements + 1; j++){
+        bytes += (double)2*UFBS2b[i][j].memory_usage(); //UFBS2b[i][j] UFBS2b[j][1]
+    }
+    if (pot_3b) {
+      for (int j = 1; j < num_of_elements + 1; j++) {
+        for (int k = j; k < num_of_elements + 1; k++) {
+          bytes += (double)2*UFBS3b[i][j][k].memory_usage();
+        }
+      }
+    }
+  }
+  
+  bytes += (double)(maxshort+1)*sizeof(int);            //neighshort, maxshort
+
+  return bytes;
+}
+
diff --git a/lammps_plugin/ML-UF3/pair_uf3.h b/lammps_plugin/ML-UF3/pair_uf3.h
new file mode 100644
index 00000000..c8ff7882
--- /dev/null
+++ b/lammps_plugin/ML-UF3/pair_uf3.h
@@ -0,0 +1,88 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, 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.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ *    Contributing authors: Ajinkya Hire(U of Florida), 
+ *                          Hendrik Kraß (U of Constance),
+ *                          Richard Hennig (U of Florida)
+ * ---------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+// clang-format off
+PairStyle(uf3,PairUF3);
+// clang-format on
+#else
+
+#ifndef LMP_PAIR_UF3_H
+#define LMP_PAIR_UF3_H
+
+#include "uf3_pair_bspline.h"
+#include "uf3_triplet_bspline.h"
+
+#include "pair.h"
+
+#include <unordered_map>
+namespace LAMMPS_NS {
+
+class PairUF3 : public Pair {
+ public:
+  PairUF3(class LAMMPS *);
+  ~PairUF3() override;
+  void compute(int, int) override;
+  void settings(int, char **) override;
+  void coeff(int, char **) override;
+  void init_style() override;
+  void init_list(int, class NeighList *) override;    // needed for ptr to full neigh list
+  double init_one(int, int) override;                 // needed for cutoff radius for neighbour list
+  double single(int, int, int, int, double, double, double, double &) override;
+
+  double memory_usage() override;
+
+ protected:
+  void uf3_read_pot_file(char *potf_name);
+  void uf3_read_pot_file(int i, int j, char *potf_name);
+  void uf3_read_pot_file(int i, int j, int k, char *potf_name);
+  int num_of_elements, nbody_flag, n2body_pot_files, n3body_pot_files, tot_pot_files;
+  int bsplines_created;
+  int coeff_matrix_dim1, coeff_matrix_dim2, coeff_matrix_dim3, coeff_matrix_elements_len;
+  bool pot_3b;
+  int ***setflag_3b;
+  double **cut, ***cut_3b, **cut_3b_list, ****min_cut_3b;
+  virtual void allocate();
+  void create_bsplines();
+  std::vector<std::vector<std::vector<double>>> n2b_knot, n2b_coeff;
+  std::vector<std::vector<std::vector<std::vector<std::vector<double>>>>> n3b_knot_matrix;
+  std::unordered_map<std::string, std::vector<std::vector<std::vector<double>>>> n3b_coeff_matrix;
+  std::vector<std::vector<uf3_pair_bspline>> UFBS2b;
+  std::vector<std::vector<std::vector<uf3_triplet_bspline>>> UFBS3b;
+  int *neighshort, maxshort;    // short neighbor list array for 3body interaction
+};
+
+}    // namespace LAMMPS_NS
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Incorrect args for pair coefficients
+
+Self-explanatory.  Check the input script or data file.
+
+*/
diff --git a/lammps_plugin/ML-UF3/uf3_bspline_basis2.cpp b/lammps_plugin/ML-UF3/uf3_bspline_basis2.cpp
new file mode 100644
index 00000000..8ae1991c
--- /dev/null
+++ b/lammps_plugin/ML-UF3/uf3_bspline_basis2.cpp
@@ -0,0 +1,88 @@
+#include "uf3_bspline_basis2.h"
+
+#include "utils.h"
+#include <vector>
+
+using namespace LAMMPS_NS;
+
+// Constructor
+// Initializes coefficients and knots
+// Requires [knots] to have length 4
+uf3_bspline_basis2::uf3_bspline_basis2(LAMMPS *ulmp, const double *knots, double coefficient)
+{
+  lmp = ulmp;
+
+  double c0, c1, c2;
+
+  c0 = coefficient *
+      (pow(knots[0], 2) /
+       (pow(knots[0], 2) - knots[0] * knots[1] - knots[0] * knots[2] + knots[1] * knots[2]));
+  c1 = coefficient *
+      (-2 * knots[0] /
+       (pow(knots[0], 2) - knots[0] * knots[1] - knots[0] * knots[2] + knots[1] * knots[2]));
+  c2 = coefficient *
+      (1 / (pow(knots[0], 2) - knots[0] * knots[1] - knots[0] * knots[2] + knots[1] * knots[2]));
+  constants.push_back(c0);
+  constants.push_back(c1);
+  constants.push_back(c2);
+  c0 = coefficient *
+      (-knots[1] * knots[3] /
+           (pow(knots[1], 2) - knots[1] * knots[2] - knots[1] * knots[3] + knots[2] * knots[3]) -
+       knots[0] * knots[2] /
+           (knots[0] * knots[1] - knots[0] * knots[2] - knots[1] * knots[2] + pow(knots[2], 2)));
+  c1 = coefficient *
+      (knots[1] /
+           (pow(knots[1], 2) - knots[1] * knots[2] - knots[1] * knots[3] + knots[2] * knots[3]) +
+       knots[3] /
+           (pow(knots[1], 2) - knots[1] * knots[2] - knots[1] * knots[3] + knots[2] * knots[3]) +
+       knots[0] /
+           (knots[0] * knots[1] - knots[0] * knots[2] - knots[1] * knots[2] + pow(knots[2], 2)) +
+       knots[2] /
+           (knots[0] * knots[1] - knots[0] * knots[2] - knots[1] * knots[2] + pow(knots[2], 2)));
+  c2 = coefficient *
+      (-1 / (pow(knots[1], 2) - knots[1] * knots[2] - knots[1] * knots[3] + knots[2] * knots[3]) -
+       1 / (knots[0] * knots[1] - knots[0] * knots[2] - knots[1] * knots[2] + pow(knots[2], 2)));
+  constants.push_back(c0);
+  constants.push_back(c1);
+  constants.push_back(c2);
+  c0 = coefficient *
+      (pow(knots[3], 2) /
+       (knots[1] * knots[2] - knots[1] * knots[3] - knots[2] * knots[3] + pow(knots[3], 2)));
+  c1 = coefficient *
+      (-2 * knots[3] /
+       (knots[1] * knots[2] - knots[1] * knots[3] - knots[2] * knots[3] + pow(knots[3], 2)));
+  c2 = coefficient *
+      (1 / (knots[1] * knots[2] - knots[1] * knots[3] - knots[2] * knots[3] + pow(knots[3], 2)));
+  constants.push_back(c0);
+  constants.push_back(c1);
+  constants.push_back(c2);
+}
+
+uf3_bspline_basis2::~uf3_bspline_basis2() {}
+
+// Evaluate outer-left part of spline
+double uf3_bspline_basis2::eval0(double rsq, double r)
+{
+  return rsq * constants[2] + r * constants[1] + constants[0];
+}
+
+// Evaluate center-left part of spline
+double uf3_bspline_basis2::eval1(double rsq, double r)
+{
+  return rsq * constants[5] + r * constants[4] + constants[3];
+}
+
+// Evaluate center-right part of spline
+double uf3_bspline_basis2::eval2(double rsq, double r)
+{
+  return rsq * constants[8] + r * constants[7] + constants[6];
+}
+
+double uf3_bspline_basis2::memory_usage()
+{
+  double bytes = 0;
+
+  bytes += (double)constants.size()*sizeof(double);
+
+  return bytes;
+}
diff --git a/lammps_plugin/ML-UF3/uf3_bspline_basis2.h b/lammps_plugin/ML-UF3/uf3_bspline_basis2.h
new file mode 100644
index 00000000..8551b097
--- /dev/null
+++ b/lammps_plugin/ML-UF3/uf3_bspline_basis2.h
@@ -0,0 +1,32 @@
+//De Boor's algorithm @
+//https://pages.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/B-spline/de-Boor.html
+//For values outside the domain,
+//extrapoltaes the left(right) hand side piece of the curve
+//Only works for bspline degree upto 3 becuase of definiation of P
+//
+#include "pointers.h"
+
+#include <vector>
+
+#ifndef UF3_BSPLINE_BASIS2_H
+#define UF3_BSPLINE_BASIS2_H
+
+namespace LAMMPS_NS {
+
+class uf3_bspline_basis2 {
+ private:
+  LAMMPS *lmp;
+  std::vector<double> constants;
+
+ public:
+  uf3_bspline_basis2(LAMMPS *ulmp, const double *knots, double coefficient);
+  ~uf3_bspline_basis2();
+  double eval0(double, double);
+  double eval1(double, double);
+  double eval2(double, double);
+
+  double memory_usage();
+};
+
+}    // namespace LAMMPS_NS
+#endif
diff --git a/lammps_plugin/ML-UF3/uf3_bspline_basis3.cpp b/lammps_plugin/ML-UF3/uf3_bspline_basis3.cpp
new file mode 100644
index 00000000..f66ac0d1
--- /dev/null
+++ b/lammps_plugin/ML-UF3/uf3_bspline_basis3.cpp
@@ -0,0 +1,324 @@
+#include "uf3_bspline_basis3.h"
+
+#include "utils.h"
+#include <vector>
+
+using namespace LAMMPS_NS;
+
+// Constructor
+// Initializes coefficients and knots
+// [knots] needs to have length 4
+uf3_bspline_basis3::uf3_bspline_basis3(LAMMPS *ulmp, const double *knots, double coefficient)
+{
+  lmp = ulmp;
+
+  double c0, c1, c2, c3;
+
+  c0 = coefficient *
+      (-pow(knots[0], 3) /
+       (-pow(knots[0], 3) + pow(knots[0], 2) * knots[1] + pow(knots[0], 2) * knots[2] +
+        pow(knots[0], 2) * knots[3] - knots[0] * knots[1] * knots[2] -
+        knots[0] * knots[1] * knots[3] - knots[0] * knots[2] * knots[3] +
+        knots[1] * knots[2] * knots[3]));
+  c1 = coefficient *
+      (3 * pow(knots[0], 2) /
+       (-pow(knots[0], 3) + pow(knots[0], 2) * knots[1] + pow(knots[0], 2) * knots[2] +
+        pow(knots[0], 2) * knots[3] - knots[0] * knots[1] * knots[2] -
+        knots[0] * knots[1] * knots[3] - knots[0] * knots[2] * knots[3] +
+        knots[1] * knots[2] * knots[3]));
+  c2 = coefficient *
+      (-3 * knots[0] /
+       (-pow(knots[0], 3) + pow(knots[0], 2) * knots[1] + pow(knots[0], 2) * knots[2] +
+        pow(knots[0], 2) * knots[3] - knots[0] * knots[1] * knots[2] -
+        knots[0] * knots[1] * knots[3] - knots[0] * knots[2] * knots[3] +
+        knots[1] * knots[2] * knots[3]));
+  c3 = coefficient *
+      (1 /
+       (-pow(knots[0], 3) + pow(knots[0], 2) * knots[1] + pow(knots[0], 2) * knots[2] +
+        pow(knots[0], 2) * knots[3] - knots[0] * knots[1] * knots[2] -
+        knots[0] * knots[1] * knots[3] - knots[0] * knots[2] * knots[3] +
+        knots[1] * knots[2] * knots[3]));
+  constants.push_back(c0);
+  constants.push_back(c1);
+  constants.push_back(c2);
+  constants.push_back(c3);
+  c0 = coefficient *
+      (pow(knots[1], 2) * knots[4] /
+           (-pow(knots[1], 3) + pow(knots[1], 2) * knots[2] + pow(knots[1], 2) * knots[3] +
+            pow(knots[1], 2) * knots[4] - knots[1] * knots[2] * knots[3] -
+            knots[1] * knots[2] * knots[4] - knots[1] * knots[3] * knots[4] +
+            knots[2] * knots[3] * knots[4]) +
+       pow(knots[0], 2) * knots[2] /
+           (-pow(knots[0], 2) * knots[1] + pow(knots[0], 2) * knots[2] +
+            knots[0] * knots[1] * knots[2] + knots[0] * knots[1] * knots[3] -
+            knots[0] * pow(knots[2], 2) - knots[0] * knots[2] * knots[3] -
+            knots[1] * knots[2] * knots[3] + pow(knots[2], 2) * knots[3]) +
+       knots[0] * knots[1] * knots[3] /
+           (-knots[0] * pow(knots[1], 2) + knots[0] * knots[1] * knots[2] +
+            knots[0] * knots[1] * knots[3] - knots[0] * knots[2] * knots[3] +
+            pow(knots[1], 2) * knots[3] - knots[1] * knots[2] * knots[3] -
+            knots[1] * pow(knots[3], 2) + knots[2] * pow(knots[3], 2)));
+  c1 = coefficient *
+      (-pow(knots[1], 2) /
+           (-pow(knots[1], 3) + pow(knots[1], 2) * knots[2] + pow(knots[1], 2) * knots[3] +
+            pow(knots[1], 2) * knots[4] - knots[1] * knots[2] * knots[3] -
+            knots[1] * knots[2] * knots[4] - knots[1] * knots[3] * knots[4] +
+            knots[2] * knots[3] * knots[4]) -
+       2 * knots[1] * knots[4] /
+           (-pow(knots[1], 3) + pow(knots[1], 2) * knots[2] + pow(knots[1], 2) * knots[3] +
+            pow(knots[1], 2) * knots[4] - knots[1] * knots[2] * knots[3] -
+            knots[1] * knots[2] * knots[4] - knots[1] * knots[3] * knots[4] +
+            knots[2] * knots[3] * knots[4]) -
+       pow(knots[0], 2) /
+           (-pow(knots[0], 2) * knots[1] + pow(knots[0], 2) * knots[2] +
+            knots[0] * knots[1] * knots[2] + knots[0] * knots[1] * knots[3] -
+            knots[0] * pow(knots[2], 2) - knots[0] * knots[2] * knots[3] -
+            knots[1] * knots[2] * knots[3] + pow(knots[2], 2) * knots[3]) -
+       2 * knots[0] * knots[2] /
+           (-pow(knots[0], 2) * knots[1] + pow(knots[0], 2) * knots[2] +
+            knots[0] * knots[1] * knots[2] + knots[0] * knots[1] * knots[3] -
+            knots[0] * pow(knots[2], 2) - knots[0] * knots[2] * knots[3] -
+            knots[1] * knots[2] * knots[3] + pow(knots[2], 2) * knots[3]) -
+       knots[0] * knots[1] /
+           (-knots[0] * pow(knots[1], 2) + knots[0] * knots[1] * knots[2] +
+            knots[0] * knots[1] * knots[3] - knots[0] * knots[2] * knots[3] +
+            pow(knots[1], 2) * knots[3] - knots[1] * knots[2] * knots[3] -
+            knots[1] * pow(knots[3], 2) + knots[2] * pow(knots[3], 2)) -
+       knots[0] * knots[3] /
+           (-knots[0] * pow(knots[1], 2) + knots[0] * knots[1] * knots[2] +
+            knots[0] * knots[1] * knots[3] - knots[0] * knots[2] * knots[3] +
+            pow(knots[1], 2) * knots[3] - knots[1] * knots[2] * knots[3] -
+            knots[1] * pow(knots[3], 2) + knots[2] * pow(knots[3], 2)) -
+       knots[1] * knots[3] /
+           (-knots[0] * pow(knots[1], 2) + knots[0] * knots[1] * knots[2] +
+            knots[0] * knots[1] * knots[3] - knots[0] * knots[2] * knots[3] +
+            pow(knots[1], 2) * knots[3] - knots[1] * knots[2] * knots[3] -
+            knots[1] * pow(knots[3], 2) + knots[2] * pow(knots[3], 2)));
+  c2 = coefficient *
+      (2 * knots[1] /
+           (-pow(knots[1], 3) + pow(knots[1], 2) * knots[2] + pow(knots[1], 2) * knots[3] +
+            pow(knots[1], 2) * knots[4] - knots[1] * knots[2] * knots[3] -
+            knots[1] * knots[2] * knots[4] - knots[1] * knots[3] * knots[4] +
+            knots[2] * knots[3] * knots[4]) +
+       knots[4] /
+           (-pow(knots[1], 3) + pow(knots[1], 2) * knots[2] + pow(knots[1], 2) * knots[3] +
+            pow(knots[1], 2) * knots[4] - knots[1] * knots[2] * knots[3] -
+            knots[1] * knots[2] * knots[4] - knots[1] * knots[3] * knots[4] +
+            knots[2] * knots[3] * knots[4]) +
+       2 * knots[0] /
+           (-pow(knots[0], 2) * knots[1] + pow(knots[0], 2) * knots[2] +
+            knots[0] * knots[1] * knots[2] + knots[0] * knots[1] * knots[3] -
+            knots[0] * pow(knots[2], 2) - knots[0] * knots[2] * knots[3] -
+            knots[1] * knots[2] * knots[3] + pow(knots[2], 2) * knots[3]) +
+       knots[2] /
+           (-pow(knots[0], 2) * knots[1] + pow(knots[0], 2) * knots[2] +
+            knots[0] * knots[1] * knots[2] + knots[0] * knots[1] * knots[3] -
+            knots[0] * pow(knots[2], 2) - knots[0] * knots[2] * knots[3] -
+            knots[1] * knots[2] * knots[3] + pow(knots[2], 2) * knots[3]) +
+       knots[0] /
+           (-knots[0] * pow(knots[1], 2) + knots[0] * knots[1] * knots[2] +
+            knots[0] * knots[1] * knots[3] - knots[0] * knots[2] * knots[3] +
+            pow(knots[1], 2) * knots[3] - knots[1] * knots[2] * knots[3] -
+            knots[1] * pow(knots[3], 2) + knots[2] * pow(knots[3], 2)) +
+       knots[1] /
+           (-knots[0] * pow(knots[1], 2) + knots[0] * knots[1] * knots[2] +
+            knots[0] * knots[1] * knots[3] - knots[0] * knots[2] * knots[3] +
+            pow(knots[1], 2) * knots[3] - knots[1] * knots[2] * knots[3] -
+            knots[1] * pow(knots[3], 2) + knots[2] * pow(knots[3], 2)) +
+       knots[3] /
+           (-knots[0] * pow(knots[1], 2) + knots[0] * knots[1] * knots[2] +
+            knots[0] * knots[1] * knots[3] - knots[0] * knots[2] * knots[3] +
+            pow(knots[1], 2) * knots[3] - knots[1] * knots[2] * knots[3] -
+            knots[1] * pow(knots[3], 2) + knots[2] * pow(knots[3], 2)));
+  c3 = coefficient *
+      (-1 /
+           (-pow(knots[1], 3) + pow(knots[1], 2) * knots[2] + pow(knots[1], 2) * knots[3] +
+            pow(knots[1], 2) * knots[4] - knots[1] * knots[2] * knots[3] -
+            knots[1] * knots[2] * knots[4] - knots[1] * knots[3] * knots[4] +
+            knots[2] * knots[3] * knots[4]) -
+       1 /
+           (-pow(knots[0], 2) * knots[1] + pow(knots[0], 2) * knots[2] +
+            knots[0] * knots[1] * knots[2] + knots[0] * knots[1] * knots[3] -
+            knots[0] * pow(knots[2], 2) - knots[0] * knots[2] * knots[3] -
+            knots[1] * knots[2] * knots[3] + pow(knots[2], 2) * knots[3]) -
+       1 /
+           (-knots[0] * pow(knots[1], 2) + knots[0] * knots[1] * knots[2] +
+            knots[0] * knots[1] * knots[3] - knots[0] * knots[2] * knots[3] +
+            pow(knots[1], 2) * knots[3] - knots[1] * knots[2] * knots[3] -
+            knots[1] * pow(knots[3], 2) + knots[2] * pow(knots[3], 2)));
+  constants.push_back(c0);
+  constants.push_back(c1);
+  constants.push_back(c2);
+  constants.push_back(c3);
+  c0 = coefficient *
+      (-knots[0] * pow(knots[3], 2) /
+           (-knots[0] * knots[1] * knots[2] + knots[0] * knots[1] * knots[3] +
+            knots[0] * knots[2] * knots[3] - knots[0] * pow(knots[3], 2) +
+            knots[1] * knots[2] * knots[3] - knots[1] * pow(knots[3], 2) -
+            knots[2] * pow(knots[3], 2) + pow(knots[3], 3)) -
+       knots[1] * knots[3] * knots[4] /
+           (-pow(knots[1], 2) * knots[2] + pow(knots[1], 2) * knots[3] +
+            knots[1] * knots[2] * knots[3] + knots[1] * knots[2] * knots[4] -
+            knots[1] * pow(knots[3], 2) - knots[1] * knots[3] * knots[4] -
+            knots[2] * knots[3] * knots[4] + pow(knots[3], 2) * knots[4]) -
+       knots[2] * pow(knots[4], 2) /
+           (-knots[1] * pow(knots[2], 2) + knots[1] * knots[2] * knots[3] +
+            knots[1] * knots[2] * knots[4] - knots[1] * knots[3] * knots[4] +
+            pow(knots[2], 2) * knots[4] - knots[2] * knots[3] * knots[4] -
+            knots[2] * pow(knots[4], 2) + knots[3] * pow(knots[4], 2)));
+  c1 = coefficient *
+      (2 * knots[0] * knots[3] /
+           (-knots[0] * knots[1] * knots[2] + knots[0] * knots[1] * knots[3] +
+            knots[0] * knots[2] * knots[3] - knots[0] * pow(knots[3], 2) +
+            knots[1] * knots[2] * knots[3] - knots[1] * pow(knots[3], 2) -
+            knots[2] * pow(knots[3], 2) + pow(knots[3], 3)) +
+       pow(knots[3], 2) /
+           (-knots[0] * knots[1] * knots[2] + knots[0] * knots[1] * knots[3] +
+            knots[0] * knots[2] * knots[3] - knots[0] * pow(knots[3], 2) +
+            knots[1] * knots[2] * knots[3] - knots[1] * pow(knots[3], 2) -
+            knots[2] * pow(knots[3], 2) + pow(knots[3], 3)) +
+       knots[1] * knots[3] /
+           (-pow(knots[1], 2) * knots[2] + pow(knots[1], 2) * knots[3] +
+            knots[1] * knots[2] * knots[3] + knots[1] * knots[2] * knots[4] -
+            knots[1] * pow(knots[3], 2) - knots[1] * knots[3] * knots[4] -
+            knots[2] * knots[3] * knots[4] + pow(knots[3], 2) * knots[4]) +
+       knots[1] * knots[4] /
+           (-pow(knots[1], 2) * knots[2] + pow(knots[1], 2) * knots[3] +
+            knots[1] * knots[2] * knots[3] + knots[1] * knots[2] * knots[4] -
+            knots[1] * pow(knots[3], 2) - knots[1] * knots[3] * knots[4] -
+            knots[2] * knots[3] * knots[4] + pow(knots[3], 2) * knots[4]) +
+       knots[3] * knots[4] /
+           (-pow(knots[1], 2) * knots[2] + pow(knots[1], 2) * knots[3] +
+            knots[1] * knots[2] * knots[3] + knots[1] * knots[2] * knots[4] -
+            knots[1] * pow(knots[3], 2) - knots[1] * knots[3] * knots[4] -
+            knots[2] * knots[3] * knots[4] + pow(knots[3], 2) * knots[4]) +
+       2 * knots[2] * knots[4] /
+           (-knots[1] * pow(knots[2], 2) + knots[1] * knots[2] * knots[3] +
+            knots[1] * knots[2] * knots[4] - knots[1] * knots[3] * knots[4] +
+            pow(knots[2], 2) * knots[4] - knots[2] * knots[3] * knots[4] -
+            knots[2] * pow(knots[4], 2) + knots[3] * pow(knots[4], 2)) +
+       pow(knots[4], 2) /
+           (-knots[1] * pow(knots[2], 2) + knots[1] * knots[2] * knots[3] +
+            knots[1] * knots[2] * knots[4] - knots[1] * knots[3] * knots[4] +
+            pow(knots[2], 2) * knots[4] - knots[2] * knots[3] * knots[4] -
+            knots[2] * pow(knots[4], 2) + knots[3] * pow(knots[4], 2)));
+  c2 = coefficient *
+      (-knots[0] /
+           (-knots[0] * knots[1] * knots[2] + knots[0] * knots[1] * knots[3] +
+            knots[0] * knots[2] * knots[3] - knots[0] * pow(knots[3], 2) +
+            knots[1] * knots[2] * knots[3] - knots[1] * pow(knots[3], 2) -
+            knots[2] * pow(knots[3], 2) + pow(knots[3], 3)) -
+       2 * knots[3] /
+           (-knots[0] * knots[1] * knots[2] + knots[0] * knots[1] * knots[3] +
+            knots[0] * knots[2] * knots[3] - knots[0] * pow(knots[3], 2) +
+            knots[1] * knots[2] * knots[3] - knots[1] * pow(knots[3], 2) -
+            knots[2] * pow(knots[3], 2) + pow(knots[3], 3)) -
+       knots[1] /
+           (-pow(knots[1], 2) * knots[2] + pow(knots[1], 2) * knots[3] +
+            knots[1] * knots[2] * knots[3] + knots[1] * knots[2] * knots[4] -
+            knots[1] * pow(knots[3], 2) - knots[1] * knots[3] * knots[4] -
+            knots[2] * knots[3] * knots[4] + pow(knots[3], 2) * knots[4]) -
+       knots[3] /
+           (-pow(knots[1], 2) * knots[2] + pow(knots[1], 2) * knots[3] +
+            knots[1] * knots[2] * knots[3] + knots[1] * knots[2] * knots[4] -
+            knots[1] * pow(knots[3], 2) - knots[1] * knots[3] * knots[4] -
+            knots[2] * knots[3] * knots[4] + pow(knots[3], 2) * knots[4]) -
+       knots[4] /
+           (-pow(knots[1], 2) * knots[2] + pow(knots[1], 2) * knots[3] +
+            knots[1] * knots[2] * knots[3] + knots[1] * knots[2] * knots[4] -
+            knots[1] * pow(knots[3], 2) - knots[1] * knots[3] * knots[4] -
+            knots[2] * knots[3] * knots[4] + pow(knots[3], 2) * knots[4]) -
+       knots[2] /
+           (-knots[1] * pow(knots[2], 2) + knots[1] * knots[2] * knots[3] +
+            knots[1] * knots[2] * knots[4] - knots[1] * knots[3] * knots[4] +
+            pow(knots[2], 2) * knots[4] - knots[2] * knots[3] * knots[4] -
+            knots[2] * pow(knots[4], 2) + knots[3] * pow(knots[4], 2)) -
+       2 * knots[4] /
+           (-knots[1] * pow(knots[2], 2) + knots[1] * knots[2] * knots[3] +
+            knots[1] * knots[2] * knots[4] - knots[1] * knots[3] * knots[4] +
+            pow(knots[2], 2) * knots[4] - knots[2] * knots[3] * knots[4] -
+            knots[2] * pow(knots[4], 2) + knots[3] * pow(knots[4], 2)));
+  c3 = coefficient *
+      (1 /
+           (-knots[0] * knots[1] * knots[2] + knots[0] * knots[1] * knots[3] +
+            knots[0] * knots[2] * knots[3] - knots[0] * pow(knots[3], 2) +
+            knots[1] * knots[2] * knots[3] - knots[1] * pow(knots[3], 2) -
+            knots[2] * pow(knots[3], 2) + pow(knots[3], 3)) +
+       1 /
+           (-pow(knots[1], 2) * knots[2] + pow(knots[1], 2) * knots[3] +
+            knots[1] * knots[2] * knots[3] + knots[1] * knots[2] * knots[4] -
+            knots[1] * pow(knots[3], 2) - knots[1] * knots[3] * knots[4] -
+            knots[2] * knots[3] * knots[4] + pow(knots[3], 2) * knots[4]) +
+       1 /
+           (-knots[1] * pow(knots[2], 2) + knots[1] * knots[2] * knots[3] +
+            knots[1] * knots[2] * knots[4] - knots[1] * knots[3] * knots[4] +
+            pow(knots[2], 2) * knots[4] - knots[2] * knots[3] * knots[4] -
+            knots[2] * pow(knots[4], 2) + knots[3] * pow(knots[4], 2)));
+  constants.push_back(c0);
+  constants.push_back(c1);
+  constants.push_back(c2);
+  constants.push_back(c3);
+  c0 = coefficient *
+      (pow(knots[4], 3) /
+       (-knots[1] * knots[2] * knots[3] + knots[1] * knots[2] * knots[4] +
+        knots[1] * knots[3] * knots[4] - knots[1] * pow(knots[4], 2) +
+        knots[2] * knots[3] * knots[4] - knots[2] * pow(knots[4], 2) - knots[3] * pow(knots[4], 2) +
+        pow(knots[4], 3)));
+  c1 = coefficient *
+      (-3 * pow(knots[4], 2) /
+       (-knots[1] * knots[2] * knots[3] + knots[1] * knots[2] * knots[4] +
+        knots[1] * knots[3] * knots[4] - knots[1] * pow(knots[4], 2) +
+        knots[2] * knots[3] * knots[4] - knots[2] * pow(knots[4], 2) - knots[3] * pow(knots[4], 2) +
+        pow(knots[4], 3)));
+  c2 = coefficient *
+      (3 * knots[4] /
+       (-knots[1] * knots[2] * knots[3] + knots[1] * knots[2] * knots[4] +
+        knots[1] * knots[3] * knots[4] - knots[1] * pow(knots[4], 2) +
+        knots[2] * knots[3] * knots[4] - knots[2] * pow(knots[4], 2) - knots[3] * pow(knots[4], 2) +
+        pow(knots[4], 3)));
+  c3 = coefficient *
+      (-1 /
+       (-knots[1] * knots[2] * knots[3] + knots[1] * knots[2] * knots[4] +
+        knots[1] * knots[3] * knots[4] - knots[1] * pow(knots[4], 2) +
+        knots[2] * knots[3] * knots[4] - knots[2] * pow(knots[4], 2) - knots[3] * pow(knots[4], 2) +
+        pow(knots[4], 3)));
+  constants.push_back(c0);
+  constants.push_back(c1);
+  constants.push_back(c2);
+  constants.push_back(c3);
+}
+
+uf3_bspline_basis3::~uf3_bspline_basis3() {}
+
+// Evaluate outer-left part of spline
+double uf3_bspline_basis3::eval0(double rth, double rsq, double r)
+{
+  return rth * constants[3] + rsq * constants[2] + r * constants[1] + constants[0];
+}
+
+// Evaluate center-left part of spline
+double uf3_bspline_basis3::eval1(double rth, double rsq, double r)
+{
+  return rth * constants[7] + rsq * constants[6] + r * constants[5] + constants[4];
+}
+
+// Evaluate center-right part of spline
+double uf3_bspline_basis3::eval2(double rth, double rsq, double r)
+{
+  return rth * constants[11] + rsq * constants[10] + r * constants[9] + constants[8];
+}
+
+// Evaluate outer-right part of spline
+double uf3_bspline_basis3::eval3(double rth, double rsq, double r)
+{
+  return rth * constants[15] + rsq * constants[14] + r * constants[13] + constants[12];
+}
+
+double uf3_bspline_basis3::memory_usage()
+{
+  double bytes = 0;
+
+  bytes += (double)constants.size()*sizeof(double);
+
+  return bytes;
+}
diff --git a/lammps_plugin/ML-UF3/uf3_bspline_basis3.h b/lammps_plugin/ML-UF3/uf3_bspline_basis3.h
new file mode 100644
index 00000000..d29d9b08
--- /dev/null
+++ b/lammps_plugin/ML-UF3/uf3_bspline_basis3.h
@@ -0,0 +1,33 @@
+//De Boor's algorithm @
+//https://pages.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/B-spline/de-Boor.html
+//For values outside the domain,
+//extrapoltaes the left(right) hand side piece of the curve
+//Only works for bspline degree upto 3 becuase of definiation of P
+//
+#include "pointers.h"
+
+#include <vector>
+
+#ifndef UF3_BSPLINE_BASIS3_H
+#define UF3_BSPLINE_BASIS3_H
+
+namespace LAMMPS_NS {
+
+class uf3_bspline_basis3 {
+ private:
+  LAMMPS *lmp;
+  std::vector<double> constants;
+
+ public:
+  uf3_bspline_basis3(LAMMPS *ulmp, const double *knots, double coefficient);
+  ~uf3_bspline_basis3();
+  double eval0(double, double, double);
+  double eval1(double, double, double);
+  double eval2(double, double, double);
+  double eval3(double, double, double);
+
+  double memory_usage();
+};
+
+}    // namespace LAMMPS_NS
+#endif
diff --git a/lammps_plugin/ML-UF3/uf3_pair_bspline.cpp b/lammps_plugin/ML-UF3/uf3_pair_bspline.cpp
new file mode 100644
index 00000000..8db46f5f
--- /dev/null
+++ b/lammps_plugin/ML-UF3/uf3_pair_bspline.cpp
@@ -0,0 +1,109 @@
+#include "uf3_pair_bspline.h"
+
+#include "uf3_bspline_basis2.h"
+#include "uf3_bspline_basis3.h"
+
+#include "utils.h"
+#include <vector>
+
+using namespace LAMMPS_NS;
+
+// Dummy constructor
+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)
+{
+  lmp = ulmp;
+  knot_vect = uknot_vect;
+  coeff_vect = ucoeff_vect;
+
+  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]));
+
+  // Initialize Coefficients and Knots for Derivatives
+  // The last coefficient needs to be droped
+  for (int i = 0; i < coeff_vect_size - 1; i++) {
+    double dntemp4 = 3 / (knot_vect[i + 4] - knot_vect[i + 1]);
+    dncoeff_vect.push_back((coeff_vect[i + 1] - coeff_vect[i]) * dntemp4);
+  }
+  //What we have is a clamped bspline -->i.e value of the bspline curve at the 
+  //knots with multiplicity equal to the degree of bspline is equal to the coefficient
+  //
+  //Therefore for the derivative bspline the very first and last knot needs to be droped
+  //to change their multiplicity from 4 (necessary condition for clamped cubic bspline)
+  //to 3 (necessary condition for clamped quadratic bspline)
+  //
+  //Also if the coeff vector size of decreases by 1 for the derivative bspline
+  //knots size needs to go down by 2 as ==> knots = coefficient + degree + 1
+  for (int i = 1; i < knot_vect_size - 1; i++) dnknot_vect.push_back(knot_vect[i]);
+
+  // Initialize B-Spline Derivative Basis Functions
+  for (int i = 0; i < dnknot_vect.size() - 3; i++)
+    dnbspline_bases.push_back(uf3_bspline_basis2(lmp, &dnknot_vect[i], dncoeff_vect[i]));
+}
+
+uf3_pair_bspline::~uf3_pair_bspline() {}
+
+double *uf3_pair_bspline::eval(double r)
+{
+
+  // Find knot starting position
+
+  int start_index;
+  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;
+
+  double rsq = r * r;
+  double rth = rsq * r;
+
+  // Calculate energy
+
+  ret_val[0] = bspline_bases[knot_affect_start + 3].eval0(rth, rsq, r);
+  ret_val[0] += bspline_bases[knot_affect_start + 2].eval1(rth, rsq, r);
+  ret_val[0] += bspline_bases[knot_affect_start + 1].eval2(rth, rsq, r);
+  ret_val[0] += bspline_bases[knot_affect_start].eval3(rth, rsq, r);
+
+  // Calculate force
+
+  ret_val[1] = dnbspline_bases[knot_affect_start + 2].eval0(rsq, r);
+  ret_val[1] += dnbspline_bases[knot_affect_start + 1].eval1(rsq, r);
+  ret_val[1] += dnbspline_bases[knot_affect_start].eval2(rsq, r);
+
+  return ret_val;
+}
+
+double uf3_pair_bspline::memory_usage()
+{
+  double bytes = 0;
+
+  bytes += (double)2*sizeof(int);                           //knot_vect_size,
+                                                            //coeff_vect_size
+  bytes += (double)knot_vect.size()*sizeof(double);         //knot_vect
+  bytes += (double)dnknot_vect.size()*sizeof(double);       //dnknot_vect
+  bytes += (double)coeff_vect.size()*sizeof(double);        //coeff_vect
+  bytes += (double)dncoeff_vect.size()*sizeof(double);      //dncoeff_vect
+
+  for (int i = 0; i < knot_vect.size() - 4; i++)
+    bytes += (double)bspline_bases[i].memory_usage();       //bspline_basis3
+
+  for (int i = 0; i < dnknot_vect.size() - 3; i++)
+    bytes += (double)dnbspline_bases[i].memory_usage();     //bspline_basis2
+
+  bytes += (double)2*sizeof(double);    //ret_val
+}
diff --git a/lammps_plugin/ML-UF3/uf3_pair_bspline.h b/lammps_plugin/ML-UF3/uf3_pair_bspline.h
new file mode 100644
index 00000000..a6d537cc
--- /dev/null
+++ b/lammps_plugin/ML-UF3/uf3_pair_bspline.h
@@ -0,0 +1,50 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/ 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.
+------------------------------------------------------------------------- */
+
+// De Boor's algorithm @
+// https://pages.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/B-spline/de-Boor.html
+// For values outside the domain, it exhibits undefined behavior.
+// Uses fixed B-Spline degree 3.
+
+#include "pointers.h"
+
+#include "uf3_bspline_basis2.h"
+#include "uf3_bspline_basis3.h"
+
+#include <vector>
+
+#ifndef UF3_PAIR_BSPLINE_H
+#define UF3_PAIR_BSPLINE_H
+
+namespace LAMMPS_NS {
+
+class uf3_pair_bspline {
+ private:
+  int knot_vect_size, coeff_vect_size;
+  std::vector<double> knot_vect, dnknot_vect;
+  std::vector<double> coeff_vect, dncoeff_vect;
+  std::vector<uf3_bspline_basis3> bspline_bases;
+  std::vector<uf3_bspline_basis2> dnbspline_bases;
+  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);
+  ~uf3_pair_bspline();
+  double ret_val[2];
+  double *eval(double value_rij);
+
+  double memory_usage();
+};
+}    // namespace LAMMPS_NS
+#endif
diff --git a/lammps_plugin/ML-UF3/uf3_triplet_bspline.cpp b/lammps_plugin/ML-UF3/uf3_triplet_bspline.cpp
new file mode 100644
index 00000000..403c6969
--- /dev/null
+++ b/lammps_plugin/ML-UF3/uf3_triplet_bspline.cpp
@@ -0,0 +1,282 @@
+#include "uf3_triplet_bspline.h"
+#include <iostream>
+#include <vector>
+
+using namespace LAMMPS_NS;
+
+// Dummy constructor
+uf3_triplet_bspline::uf3_triplet_bspline(){};
+
+// Construct a new 3D B-Spline
+uf3_triplet_bspline::uf3_triplet_bspline(
+    LAMMPS *ulmp, const std::vector<std::vector<double>> &uknot_matrix,
+    const std::vector<std::vector<std::vector<double>>> &ucoeff_matrix)
+{
+  lmp = ulmp;
+  knot_matrix = uknot_matrix;
+  coeff_matrix = ucoeff_matrix;
+
+  knot_vect_size_ij = knot_matrix[2].size();
+  knot_vect_size_ik = knot_matrix[1].size();
+  knot_vect_size_jk = knot_matrix[0].size();
+
+  int resolution_ij = knot_vect_size_ij - 4;
+  int resolution_ik = knot_vect_size_ik - 4;
+  int resolution_jk = knot_vect_size_jk - 4;
+
+  // Cache Spline Basis Functions
+  for (int l = 0; l < resolution_ij; l++) {
+    bsplines_ij.push_back(uf3_bspline_basis3(lmp, &knot_matrix[2][l], 1));
+  }
+
+  for (int l = 0; l < resolution_ik; l++) {
+    // Reuse jk Basis if Knots match
+    if (knot_matrix[1][l] == knot_matrix[2][l] && knot_matrix[1][l + 1] == knot_matrix[2][l + 1] &&
+        knot_matrix[1][l + 2] == knot_matrix[2][l + 2] &&
+        knot_matrix[1][l + 3] == knot_matrix[2][l + 3])
+      bsplines_ik.push_back(bsplines_ij[l]);
+    else
+      bsplines_ik.push_back(uf3_bspline_basis3(lmp, &knot_matrix[1][l], 1));
+  }
+
+  for (int l = 0; l < resolution_jk; l++) {
+    bsplines_jk.push_back(uf3_bspline_basis3(lmp, &knot_matrix[0][l], 1));
+  }
+
+  // Initialize Coefficients for Derivatives
+  for (int i = 0; i < coeff_matrix.size(); i++) {
+    std::vector<std::vector<double>> dncoeff_vect2;
+    for (int j = 0; j < coeff_matrix[0].size(); j++) {
+      std::vector<double> dncoeff_vect;
+      for (int k = 0; k < coeff_matrix[0][0].size() - 1; k++) {
+        double dntemp4 = 3 / (knot_matrix[0][k + 4] - knot_matrix[0][k + 1]);
+        dncoeff_vect.push_back((coeff_matrix[i][j][k + 1] - coeff_matrix[i][j][k]) * dntemp4);
+      }
+      dncoeff_vect2.push_back(dncoeff_vect);
+    }
+    dncoeff_matrix_jk.push_back(dncoeff_vect2);
+  }
+
+  for (int i = 0; i < coeff_matrix.size(); i++) {
+    std::vector<std::vector<double>> dncoeff_vect2;
+    for (int j = 0; j < coeff_matrix[0].size() - 1; j++) {
+      double dntemp4 = 3 / (knot_matrix[1][j + 4] - knot_matrix[1][j + 1]);
+      std::vector<double> dncoeff_vect;
+      for (int k = 0; k < coeff_matrix[0][0].size(); k++) {
+        dncoeff_vect.push_back((coeff_matrix[i][j + 1][k] - coeff_matrix[i][j][k]) * dntemp4);
+      }
+      dncoeff_vect2.push_back(dncoeff_vect);
+    }
+    dncoeff_matrix_ik.push_back(dncoeff_vect2);
+  }
+
+  for (int i = 0; i < coeff_matrix.size() - 1; i++) {
+    std::vector<std::vector<double>> dncoeff_vect2;
+    double dntemp4 = 3 / (knot_matrix[2][i + 4] - knot_matrix[2][i + 1]);
+    for (int j = 0; j < coeff_matrix[0].size(); j++) {
+      std::vector<double> dncoeff_vect;
+      for (int k = 0; k < coeff_matrix[0][0].size(); k++) {
+        dncoeff_vect.push_back((coeff_matrix[i + 1][j][k] - coeff_matrix[i][j][k]) * dntemp4);
+      }
+      dncoeff_vect2.push_back(dncoeff_vect);
+    }
+    dncoeff_matrix_ij.push_back(dncoeff_vect2);
+  }
+
+  std::vector<std::vector<double>> dnknot_matrix;
+  for (int i = 0; i < knot_matrix.size(); i++) {
+    std::vector<double> dnknot_vect;
+    for (int j = 1; j < knot_matrix[0].size() - 1; j++) {
+      dnknot_vect.push_back(knot_matrix[i][j]);
+    }
+    dnknot_matrix.push_back(dnknot_vect);
+  }
+
+  // Cache Derivative Spline Basis Functions
+  for (int l = 0; l < resolution_ij - 1; l++) {
+    dnbsplines_ij.push_back(uf3_bspline_basis2(lmp, &dnknot_matrix[2][l], 1));
+  }
+
+  for (int l = 0; l < resolution_ik - 1; l++) {
+    // Reuse jk Basis if Knots match
+    if (dnknot_matrix[1][l] == dnknot_matrix[2][l] &&
+        dnknot_matrix[1][l + 1] == dnknot_matrix[2][l + 1] &&
+        dnknot_matrix[1][l + 2] == dnknot_matrix[2][l + 2])
+      dnbsplines_ik.push_back(dnbsplines_ij[l]);
+    else
+      dnbsplines_ik.push_back(uf3_bspline_basis2(lmp, &dnknot_matrix[1][l], 1));
+  }
+
+  for (int l = 0; l < resolution_jk - 1; l++) {
+    dnbsplines_jk.push_back(uf3_bspline_basis2(lmp, &dnknot_matrix[0][l], 1));
+  }
+}
+
+// Destructor
+uf3_triplet_bspline::~uf3_triplet_bspline() {}
+
+// Evaluate 3D B-Spline value
+double *uf3_triplet_bspline::eval(double value_rij, double value_rik, double value_rjk)
+{
+
+  // Find starting knots
+
+  int iknot_ij = starting_knot(knot_matrix[2], knot_vect_size_ij, value_rij) - 3;
+  int iknot_ik = starting_knot(knot_matrix[1], knot_vect_size_ik, value_rik) - 3;
+  int iknot_jk = starting_knot(knot_matrix[0], knot_vect_size_jk, value_rjk) - 3;
+
+  double rsq_ij = value_rij * value_rij;
+  double rsq_ik = value_rik * value_rik;
+  double rsq_jk = value_rjk * value_rjk;
+  double rth_ij = rsq_ij * value_rij;
+  double rth_ik = rsq_ik * value_rik;
+  double rth_jk = rsq_jk * value_rjk;
+
+  // Calculate energies
+
+  double basis_ij[4];
+  basis_ij[0] = bsplines_ij[iknot_ij].eval3(rth_ij, rsq_ij, value_rij);
+  basis_ij[1] = bsplines_ij[iknot_ij + 1].eval2(rth_ij, rsq_ij, value_rij);
+  basis_ij[2] = bsplines_ij[iknot_ij + 2].eval1(rth_ij, rsq_ij, value_rij);
+  basis_ij[3] = bsplines_ij[iknot_ij + 3].eval0(rth_ij, rsq_ij, value_rij);
+
+  double basis_ik[4];
+  basis_ik[0] = bsplines_ik[iknot_ik].eval3(rth_ik, rsq_ik, value_rik);
+  basis_ik[1] = bsplines_ik[iknot_ik + 1].eval2(rth_ik, rsq_ik, value_rik);
+  basis_ik[2] = bsplines_ik[iknot_ik + 2].eval1(rth_ik, rsq_ik, value_rik);
+  basis_ik[3] = bsplines_ik[iknot_ik + 3].eval0(rth_ik, rsq_ik, value_rik);
+
+  double basis_jk[4];
+  basis_jk[0] = bsplines_jk[iknot_jk].eval3(rth_jk, rsq_jk, value_rjk);
+  basis_jk[1] = bsplines_jk[iknot_jk + 1].eval2(rth_jk, rsq_jk, value_rjk);
+  basis_jk[2] = bsplines_jk[iknot_jk + 2].eval1(rth_jk, rsq_jk, value_rjk);
+  basis_jk[3] = bsplines_jk[iknot_jk + 3].eval0(rth_jk, rsq_jk, value_rjk);
+
+  ret_val[0] = 0;
+  ret_val[1] = 0;
+  ret_val[2] = 0;
+  ret_val[3] = 0;
+
+  for (int i = 0; i < 4; i++) {
+    for (int j = 0; j < 4; j++) {
+      for (int k = 0; k < 4; k++) {
+        ret_val[0] += coeff_matrix[i + iknot_ij][j + iknot_ik][k + iknot_jk] * basis_ij[i] *
+            basis_ik[j] * basis_jk[k];
+      }
+    }
+  }
+
+  // Calculate forces
+
+  double dnbasis_ij[4];
+  dnbasis_ij[0] = dnbsplines_ij[iknot_ij].eval2(rsq_ij, value_rij);
+  dnbasis_ij[1] = dnbsplines_ij[iknot_ij + 1].eval1(rsq_ij, value_rij);
+  dnbasis_ij[2] = dnbsplines_ij[iknot_ij + 2].eval0(rsq_ij, value_rij);
+  dnbasis_ij[3] = 0;
+
+  double dnbasis_ik[4];
+  dnbasis_ik[0] = dnbsplines_ik[iknot_ik].eval2(rsq_ik, value_rik);
+  dnbasis_ik[1] = dnbsplines_ik[iknot_ik + 1].eval1(rsq_ik, value_rik);
+  dnbasis_ik[2] = dnbsplines_ik[iknot_ik + 2].eval0(rsq_ik, value_rik);
+  dnbasis_ik[3] = 0;
+
+  double dnbasis_jk[4];
+  dnbasis_jk[0] = dnbsplines_jk[iknot_jk].eval2(rsq_jk, value_rjk);
+  dnbasis_jk[1] = dnbsplines_jk[iknot_jk + 1].eval1(rsq_jk, value_rjk);
+  dnbasis_jk[2] = dnbsplines_jk[iknot_jk + 2].eval0(rsq_jk, value_rjk);
+  dnbasis_jk[3] = 0;
+
+  for (int i = 0; i < 3; i++) {
+    for (int j = 0; j < 4; j++) {
+      for (int k = 0; k < 4; k++) {
+        ret_val[1] += dncoeff_matrix_ij[iknot_ij + i][iknot_ik + j][iknot_jk + k] * dnbasis_ij[i] *
+            basis_ik[j] * basis_jk[k];
+      }
+    }
+  }
+
+  for (int i = 0; i < 4; i++) {
+    for (int j = 0; j < 3; j++) {
+      for (int k = 0; k < 4; k++) {
+        ret_val[2] += dncoeff_matrix_ik[iknot_ij + i][iknot_ik + j][iknot_jk + k] * basis_ij[i] *
+            dnbasis_ik[j] * basis_jk[k];
+      }
+    }
+  }
+
+  for (int i = 0; i < 4; i++) {
+    for (int j = 0; j < 4; j++) {
+      for (int k = 0; k < 3; k++) {
+        ret_val[3] += dncoeff_matrix_jk[iknot_ij + i][iknot_ik + j][iknot_jk + k] * basis_ij[i] *
+            basis_ik[j] * dnbasis_jk[k];
+      }
+    }
+  }
+
+  return ret_val;
+}
+
+// Find starting knot for spline evaluation
+
+int uf3_triplet_bspline::starting_knot(const std::vector<double> knot_vect, int knot_vect_size,
+                                       double r)
+{
+  if (knot_vect.front() <= r && r < knot_vect.back()) {
+    for (int i = 3; i < knot_vect_size - 1; i++) {
+      if (knot_vect[i] <= r && r < knot_vect[i + 1]) return i;
+    }
+  }
+
+  return 0;
+}
+
+double uf3_triplet_bspline::memory_usage()
+{
+  double bytes = 0;
+
+  bytes += (double) 3*sizeof(int);          //knot_vect_size_ij,
+                                            //knot_vect_size_ik,
+                                            //knot_vect_size_jk;
+
+  for (int i=0; i<coeff_matrix.size(); i++)
+    for (int j=0; j<coeff_matrix[i].size(); j++)
+      bytes += (double)coeff_matrix[i][j].size()*sizeof(double);
+
+
+  for (int i=0; i<dncoeff_matrix_ij.size(); i++)
+    for (int j=0; j<dncoeff_matrix_ij[i].size(); j++)
+      bytes += (double)dncoeff_matrix_ij[i][j].size()*sizeof(double);
+
+  for (int i=0; i<dncoeff_matrix_ik.size(); i++)
+    for (int j=0; j<dncoeff_matrix_ik[i].size(); j++)
+      bytes += (double)dncoeff_matrix_ik[i][j].size()*sizeof(double);
+
+  for (int i=0; i<dncoeff_matrix_jk.size(); i++)
+    for (int j=0; j<dncoeff_matrix_jk[i].size(); j++)
+      bytes += (double)dncoeff_matrix_jk[i][j].size()*sizeof(double);
+
+  bytes += (double)knot_matrix[0].size()*sizeof(double);
+  bytes += (double)knot_matrix[1].size()*sizeof(double);
+  bytes += (double)knot_matrix[2].size()*sizeof(double);
+
+  for (int i=0; i < bsplines_ij.size(); i++)
+    bytes += (double)bsplines_ij[i].memory_usage();
+
+  for (int i=0; i < bsplines_ik.size(); i++)
+    bytes += (double)bsplines_ik[i].memory_usage();
+
+  for (int i=0; i < bsplines_jk.size(); i++)
+    bytes += (double)bsplines_jk[i].memory_usage();
+
+  for (int i=0; i < dnbsplines_ij.size(); i++)
+    bytes += (double)dnbsplines_ij[i].memory_usage();
+
+  for (int i=0; i < dnbsplines_ik.size(); i++)
+    bytes += (double)dnbsplines_ik[i].memory_usage();
+
+  for (int i=0; i < dnbsplines_jk.size(); i++)
+    bytes += (double)dnbsplines_jk[i].memory_usage();
+
+  bytes += (double)4*sizeof(double);                        //ret_val
+  return bytes;
+}
diff --git a/lammps_plugin/ML-UF3/uf3_triplet_bspline.h b/lammps_plugin/ML-UF3/uf3_triplet_bspline.h
new file mode 100644
index 00000000..dbb7f2c8
--- /dev/null
+++ b/lammps_plugin/ML-UF3/uf3_triplet_bspline.h
@@ -0,0 +1,48 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/ 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 "pointers.h"
+
+#include "uf3_bspline_basis2.h"
+#include "uf3_bspline_basis3.h"
+#include "uf3_pair_bspline.h"
+
+#include <vector>
+
+#ifndef UF3_TRIPLET_BSPLINE_H
+#define UF3_TRIPLET_BSPLINE_H
+
+namespace LAMMPS_NS {
+class uf3_triplet_bspline {
+ private:
+  LAMMPS *lmp;
+  int knot_vect_size_ij, knot_vect_size_ik, knot_vect_size_jk;
+  std::vector<std::vector<std::vector<double>>> coeff_matrix, dncoeff_matrix_ij, dncoeff_matrix_ik,
+      dncoeff_matrix_jk;
+  std::vector<std::vector<double>> knot_matrix;
+  std::vector<uf3_bspline_basis3> bsplines_ij, bsplines_ik, bsplines_jk;
+  std::vector<uf3_bspline_basis2> dnbsplines_ij, dnbsplines_ik, dnbsplines_jk;
+  double ret_val[4];
+
+  int starting_knot(const std::vector<double>, int, double);
+
+ public:
+  //Dummy Constructor
+  uf3_triplet_bspline();
+  uf3_triplet_bspline(LAMMPS *ulmp, const std::vector<std::vector<double>> &uknot_matrix,
+                      const std::vector<std::vector<std::vector<double>>> &ucoeff_matrix);
+  ~uf3_triplet_bspline();
+  double *eval(double value_rij, double value_rik, double value_rjk);
+
+  double memory_usage();
+};
+}    // namespace LAMMPS_NS
+#endif
diff --git a/lammps_plugin/README.rst b/lammps_plugin/README.rst
new file mode 100644
index 00000000..36ff281b
--- /dev/null
+++ b/lammps_plugin/README.rst
@@ -0,0 +1,164 @@
+This document describes how to use UF3 potentials in `lammps <https://www.lammps.org/>`_ MD code. See the `tungsten_example <https://github.com/monk-04/uf3/tree/lammps_implementation/lammps_plugin/tungsten_example>`_ directory for an example of lammps input file and UF3 lammps potential files.
+
+.. contents:: Contents
+	:depth: 1
+	:local: 
+
+=====
+Compiling lammps with UF3 library
+=====
+
+Before running lammps with UF3 potentials, lammps must be re-compiled with the :code:`pair_uf3` and other supporting libraries contained in :code:`ML-UF3` directory.
+
+Copy the :code:`ML-UF3` directory to the appropriate lammps source code directory.
+
+.. code:: bash
+
+   mv ML-UF3 LAMMPS_BASE_DIR/src/.
+
+If compiling lammps with :code:`CMake`, add the :code:`ML-UF3` keyword to :code:`set(STANDARD_PACKAGES)` list in :code:`LAMMPS_BASE_DIR/cmake/CMakeLists.txt` file (search for :code:`ML-SNAP` and add :code:`ML-UF3` in the very next line). Go to the :code:`build` directory in :code:`LAMMPS_BASE_DIR` (make one if it doesn't exist) and compile lammps.
+
+.. code:: bash
+
+ cd LAMMPS_BASE_DIR/build
+ cmake ../cmake/ -D PKG_ML-UF3=yes
+ cmake --build .
+
+This will compile lammps with support for UF3 potentials. For more information on how to build lammps see this link_.
+
+.. _link: https://docs.lammps.org/Build.html
+
+If you want to compile with MPI support enabled, you additionally have to set the BUILD_MPI flag and make sure that a suitable MPI runtime can be found by the compiler. To make the executables distinguishable, the `LAMMPS_MACHINE` option sets a suffix for the output file.
+
+.. code:: bash
+
+ cd LAMMPS_BASE_DIR/build
+ cmake ../cmake/ -D PKG_ML-UF3=yes -D BUILD_MPI=yes -D LAMMPS_MACHINE=mpi # Compiles to lmp_mpi
+ cmake --build .
+
+Kokkos
+=====
+The kokkos accelerator variants of uf3 (i.e. :code:`pair_style uf3/kk` with GPU support) is under active development. For an early access refer to the :code:`README` at this link_.
+
+.. _link: https://github.com/monk-04/uf3/tree/lammps_implementation_v2/lammps_plugin
+
+=====
+Running lammps with UF3 potential
+=====
+
+**The old style of listing all UF3 LAMMPS potential files after pair_coeff * * in single line still works but a warning is printed**
+
+To use UF3 potentials in lammps just add the following tags to the lammps input file-
+
+.. code:: bash
+
+    pair_style uf3 3 1
+    pair_coeff 1 1 W_W 
+    pair_coeff 1 1 1 W_W_W
+
+The 'uf3' keyword in :code:`pair_style` invokes the UF3 potentials in lammps. The number next to the :code:`uf3` keyword tells lammps whether the user wants to run the MD code with just 2-body or 2 and 3-body UF3 potentials. The last number of this line specifies the number of elemnts in the system. So in the above example, the user wants to run MD simulation with UF3 potentials containing both 2-body and 3-body interactions on a system containing only 1 element.
+
+The :code:`pair_coeff` tag is used to read in the user-provided UF3 lammps potential files. These files can be generated directly from the :code:`json` potential files of UF3. We recommend using the :code:`generate_uf3_lammps_pots.py` script (`found here <https://github.com/monk-04/uf3/tree/lammps_implementation/lammps_plugin/scripts>`_) for generating the UF3 lammps potential files. It will also additionally print lines that should be added to the lammps input file for using UF3 lammps potential files.
+
+After :code:`pair_coeff` specify the interactions (two numbers for 2-body, three numbers for 3-body) followed by the name of the potential file. The user can also use asterisks:code:`*` for wild-card characters. In this case the behaviour is similar to other LAMMPS :code:`pair_style` for example LJ. The user can also specify. Make sure these files are present in the current run directory or in directories where lammps can find them.
+
+As an example for a multicomponet system containing elements 'A' and 'B' the above lines can be-
+
+.. code:: bash
+
+   pair_style uf3 3 2
+   pair_coeff 1 1 A_A
+   pair_coeff 1 2 A_B
+   pair_coeff 2 2 B_B
+   pair_coeff 1 1 1 A_A_A
+   pair_coeff 1 1 2 A_A_B
+   pair_coeff 1 2 2 A_B_B
+   pair_coeff 2 1 1 B_A_A
+   pair_coeff 2 1 2 B_A_B
+   pair_coeff 2 2 2 B_B_B
+
+Following format is also a valid for system containing elements 'A' and 'B'
+
+.. code:: bash
+
+   pair_style uf3 3 2
+   pair_coeff * * A_A
+   pair_coeff 1 * * A_A_A
+   pair_coeff 2 * * B_B_B
+
+   
+Alternatively, if the user wishes to use only the 2-body interactions from a model containing both two and three body interaction simply change the number next to :code:`uf3` to :code:`2` and don't list the three body interaction files in the :code:`pair_coeff` line. Beware! Using only the 2-body interaction from a model containing both 2 and 3-body is not recommended and will give wrong results!
+
+.. code:: bash
+  pair_style uf3 2 2
+  pair_coeff 1 1 A_A
+  pair_coeff 1 2 A_B
+  pair_coeff 2 2 B_B
+  
+=====
+Structure of UF3 lammps potential file
+=====
+
+This section describes the format of the UF3 lammps potential file. Not following the format can lead to unexpected error in the MD simulation and sometimes unexplained core dumps.
+
+
+2-body potential
+====
+
+**The old UF3 LAMMPS potential files can still be used but a warning is printed**
+
+The 2-body UF3 lammps potential file should have the following format-
+
+.. code:: bash
+
+    #UF3 POT
+    2B LEADING_TRIM TRAILING_TRIM
+    Rij_CUTOFF NUM_OF_KNOTS
+    BSPLINE_KNOTS
+    NUM_OF_COEFF
+    COEFF
+    #
+
+The first line of all UF3 lammps potential files should start with :code:`#UF3 POT` characters. The next line indicates whether the file contains UF3 lammps potential data for 2-body (:code:`2B`) or 3-body (:code:`3B`) interaction. This is followed by :code:`LEADING_TRIM` and :code:`TRAILING_TRIM` number. The current implementation is only tested for :code:`LEADING_TRIM=0` and :code:`TRAILING_TRIM=3`. If other values are used LAMMPS will stop with an error message.
+
+The :code:`Rij_CUTOFF` sets the 2-body cutoff for the interaction described by the potential file. :code:`NUM_OF_KNOTS` is the number of knots (or the length of the knot vector) present on the very next line. The :code:`BSPLINE_KNOTS` line should contain all the knots in ascending order. :code:`NUM_OF_COEFF` is the number of coefficients in the :code:`COEFF` line. All the numbers in the BSPLINE_KNOTS and COEFF line should be space-separated. 
+
+3-body potential
+====
+
+**The old UF3 LAMMPS potential files can still be used but a warning is printed**
+
+The 3-body UF3 lammps potential file has a format similar to the 2-body potential file-
+
+.. code:: bash
+
+    #UF3 POT
+    3B LEADING_TRIM TRAILING_TRIM
+    Rjk_CUTOFF Rik_CUTOFF Rij_CUTOFF NUM_OF_KNOTS_JK NUM_OF_KNOTS_IK NUM_OF_KNOTS_IJ
+    BSPLINE_KNOTS_FOR_JK
+    BSPLINE_KNOTS_FOR_IK
+    BSPLINE_KNOTS_FOR_IJ
+    SHAPE_OF_COEFF_MATRIX[I][J][K]
+    COEFF_MATRIX[0][0][K]
+    COEFF_MATRIX[0][1][K]
+    COEFF_MATRIX[0][2][K]
+    .
+    .
+    .
+    COEFF_MATRIX[1][0][K]
+    COEFF_MATRIX[1][1][K]
+    COEFF_MATRIX[1][2][K]
+    .
+    .
+    .
+    #
+
+
+The first line is similar to the 2-body potential file and the second line has :code:`3B` characters indicating that this file describes 3-body interaction. The :code:`3B` is followed by :code:`LEADING_TRIM` and :code:`TRAILING_TRIM` number. The current implementation is only tested for :code:`LEADING_TRIM=0` and :code:`TRAILING_TRIM=3`. If other values are used LAMMPS will stop with an error message.
+
+Similar to the 2-body potential file, the third line sets the cutoffs and length of the knots. The cutoff distance between atom-type 1 and 2 is :code:`Rij_CUTOFF`, atom-type 1 and 3 is :code:`Rik_CUTOFF` and between 2 and 3 is :code:`Rjk_CUTOFF`. **Note the current implementation works only for UF3 potentials with cutoff distances for 3-body interactions that follows** :code:`2Rij_CUTOFF=2Rik_CUTOFF=Rjk_CUTOFF` **relation.**
+
+The :code:`BSPLINE_KNOTS_FOR_JK`, :code:`BSPLINE_KNOTS_FOR_IK`, and :code:`BSPLINE_KNOTS_FOR_IJ` lines (note the order) contain the knots in increasing order for atoms J and K, I and K, and atoms I and J respectively. The number of knots is defined by the :code:`NUM_OF_KNOTS_*` characters in the previous line.
+The shape of the coefficient matrix is defined on the :code:`SHAPE_OF_COEFF_MATRIX[I][J][K]` line followed by the columns of the coefficient matrix, one per line, as shown above. For example, if the coefficient matrix has the shape of 8x8x13, then :code:`SHAPE_OF_COEFF_MATRIX[I][J][K]` will be :code:`8 8 13` followed by 64 (8x8) lines each containing 13 coefficients seperated by space.
+
+All the UF3 lammps potential files end with :code:`#` character.
diff --git a/lammps_plugin/documentation_lammps_format/pair_uf3.rst b/lammps_plugin/documentation_lammps_format/pair_uf3.rst
new file mode 100644
index 00000000..4e57d08b
--- /dev/null
+++ b/lammps_plugin/documentation_lammps_format/pair_uf3.rst
@@ -0,0 +1,177 @@
+.. index:: pair_style uf3
+.. index:: pair_style uf3/kk
+
+pair_style uf3 command
+======================
+
+Accelerator Variants: *uf3/kk*
+
+Syntax
+""""""
+
+.. code-block:: LAMMPS
+
+    pair_style style BodyFlag NumAtomType
+
+ * style = *uf3* or *uf3/kk*
+ 
+   .. parsed-literal::
+
+        BodyFlag = Indicates whether to calculate only 2-body or 2 and 3-body interactions. Possible values- 2 or 3
+        NumAtomType = Number of atoms types in the simulation
+
+
+
+Examples
+""""""""
+
+.. code-block:: LAMMPS
+
+    pair_style uf3 3 1
+    pair_coeff 1 1 Nb_Nb
+    pair_coeff 1 1 1 Nb_Nb_Nb
+
+    pair_style uf3 2 2
+    pair_coeff 1 1 Nb_Nb
+    pair_coeff 1 2 Nb_Sn
+    pair_coeff 2 2 Sn_Sn
+
+    pair_style uf3 3 2
+    pair_coeff 1 * Nb_Nb
+    pair_coeff 2 * Sn_Sn
+    pair_style 1 * * Nb_Sn_Sn
+    pair_style 2 * * Sn_Sn_Sn
+
+Description
+"""""""""""
+
+The *uf3* style computes the :ref:`Ultra-Fast Force Fields (UF3) <Xie23>` potential, a machine-learning interatomic potential. In UF3, the total energy of the system is defined via two- and three-body interactions:
+
+.. math::
+
+    E & = \sum_{i,j} V_2(r_{ij}) + \sum_{i,j,k} V_3 (r_{ij},r_{ik},r_{jk})
+
+    V_2(r_{ij}) & = \sum_{n=0}^K c_n B_n(r_{ij})
+
+    V_3 (r_{ij},r_{ik},r_{jk}) & = \sum_{l=0}^K_l \sum_{m=0}^K_m \sum_{n=0}^K_n c_{l,m,n} B_l(r_{ij}) B_m(r_{ik}) B_n(r_{jk})
+
+where :math:`V_2(r_{ij})` and :math:`V_3 (r_{ij},r_{ik},r_{jk})` are the two- and three-body interactions. For the two-body the summation is over all neighbours J and for the three-body the summation is over all neighbors J and K of atom I within a cutoff distance determined from the potential files. :math:`B_n(r_{ij})` are the cubic bspline basis, :math:`c_n` and :math:`c_{l,m,n}` are the machine-learned interaction parameters and :math:`K`, :math:`K_l`, :math:`K_m`, and :math:`K_n` denote the number of basis functions per spline or tensor spline dimension.
+
+The UF3 LAMMPS potential files are provided using single or multiple pair_coeff commands. 
+
+.. note::
+
+    Unlike other MANYBODY and ML potentials in LAMMPS, the atom type for which the specified potential file should be used for is not determined from the potential file. But is rather determined from the user provided atom type numbers after pair_coeff.
+
+As an example, if a LAMMPS simulation contains 2 atom types (elements 'A' and 'B'), the pair_coeff command will be-
+
+.. code-block:: LAMMPS
+
+    pair_style uf3 3 2
+    pair_coeff 1 1 A_A
+    pair_coeff 1 2 A_B
+    pair_coeff 2 2 B_B
+    pair_coeff 1 1 1 A_A_A
+    pair_coeff 1 1 2 A_A_B
+    pair_coeff 1 2 2 A_B_B
+    pair_coeff 2 1 1 B_A_A
+    pair_coeff 2 1 2 B_A_B
+    pair_coeff 2 2 2 B_B_B
+
+LAMMPS wild-card character "*" can also be used to specify a single UF3 LAMMPS potential file for multiple interaction. For example- 
+
+.. code-block:: LAMMPS
+
+    pair_style uf3 3 2
+    pair_coeff * * A_A
+    pair_coeff 1 * * A_A_A
+    pair_coeff 2 * * B_B_B
+
+The file A_A will be used for 2-body interaction between atom types 1-1, 1-2 and 2-2; file A_A_A will be used 3-body interaction for atom types 1-1-1, 1-1-2, 1-2-2; and so on. If a value of "2" is specified in the :code:`pair_style uf3` command, only the two-body potential files are needed. For 3-body interaction the first atom type is the central atom. We recommend using the :code:`generate_uf3_lammps_pots.py` script (found `here <https://github.com/monk-04/uf3/tree/lammps_impl    ementation/lammps_plugin/scripts>`_) for generating the UF3 LAMMPS potential files from the UF3 JSON potentials.
+
+UF3 LAMMPS potential files in the *potentials* directory of the LAMMPS distribution have a ".uf3" suffix. All UF3 LAMMPS potential files should start with :code:`#UF3 POT` and with :code:`#` characters. Following shows the format of a generic 2-body UF3 LAMMPS potential file-
+
+.. code-block:: LAMMPS
+
+    #UF3 POT
+    2B LEADING_TRIM TRAILING_TRIM
+    Rij_CUTOFF NUM_OF_KNOTS
+    BSPLINE_KNOTS
+    NUM_OF_COEFF
+    COEFF
+    #
+
+The second line indicates whether the potential file contains data for 2-body (:code:`2B`) or 3-body (:code:`3B`) interaction. This is followed by :code:`LEADING_TRIM` and :code:`TRAILING_TRIM` number on the same line. The current implementation is only tested for :code:`LEADING_TRIM=0` and :code:`TRAILING_TRIM=3`. If other values are used LAMMPS is terminated after issuing an error message. The :code:`Rij_CUTOFF` sets the 2-body cutoff for the interaction described by the potential file. :code:`NUM_OF_KNOTS` is the number of knots (or the length of the knot vector) present on the very next line. The :code:`BSPLINE_KNOTS` line should contain all the knots in ascending order. :code:`NUM_OF_COEFF` is the number of coefficients in the :code:`COEFF` line. All the numbers in the BSPLINE_KNOTS and COEFF line should be space-separated.
+
+The format of a generic 3-body UF3 LAMMPS potential file is as follow-
+
+.. code-block:: LAMMPS
+    #UF3 POT
+    3B LEADING_TRIM TRAILING_TRIM
+    Rjk_CUTOFF Rik_CUTOFF Rij_CUTOFF NUM_OF_KNOTS_JK NUM_OF_KNOTS_IK NUM_OF_KNOTS_IJ
+    BSPLINE_KNOTS_FOR_JK
+    BSPLINE_KNOTS_FOR_IK
+    BSPLINE_KNOTS_FOR_IJ
+    SHAPE_OF_COEFF_MATRIX[I][J][K]
+    COEFF_MATRIX[0][0][K]
+    COEFF_MATRIX[0][1][K]
+    COEFF_MATRIX[0][2][K]
+    .
+    .
+    .
+    COEFF_MATRIX[1][0][K]
+    COEFF_MATRIX[1][1][K]
+    COEFF_MATRIX[1][2][K]
+    .
+    .
+    .
+    #
+
+Similar to the 2-body potential file, the third line sets the cutoffs and length of the knots. The cutoff distance between atom-type I and J is :code:`Rij_CUTOFF`, atom-type I and K is :code:`Rik_CUTOFF` and between J and K is :code:`Rjk_CUTOFF`.
+
+.. note::
+
+    The current implementation only works for UF3 potentials with cutoff distances for 3-body interactions that follows :code:`2Rij_CUTOFF=2Rik_CUTOFF=Rjk_CUTOFF` relation.
+
+The :code:`BSPLINE_KNOTS_FOR_JK`, :code:`BSPLINE_KNOTS_FOR_IK`, and :code:`BSPLINE_KNOTS_FOR_IJ` lines (note the order) contain the knots in increasing order for atoms J and K, I and K, and atoms I and J respectively. The number of knots is defined by the :code:`NUM_OF_KNOTS_*` characters in the previous line.
+The shape of the coefficient matrix is defined on the :code:`SHAPE_OF_COEFF_MATRIX[I][J][K]` line followed by the columns of the coefficient matrix, one per line, as shown above. For example, if the coefficient matrix has the shape of 8x8x13, then :code:`SHAPE_OF_COEFF_MATRIX[I][J][K]` will be :code:`8 8 13` followed by 64 (8x8) lines each containing 13 coefficients seperated by space.
+
+
+Mixing, shift, table, tail correction, restart, rRESPA info
+"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+For atom type pairs I,J and I != J, where types I and J correspond to two different element types, mixing is performed by LAMMPS as described above from values in the potential file.
+
+This pair style does not support the :doc:`pair_modify <pair_modify>` shift, table, and tail options.
+
+This pair style does not write its information to :doc:`binary restart files <restart>`, since it is stored in potential files.
+
+This pair style can only be used via the *pair* keyword of the :doc:`run_style respa <run_style>` command.  It does not support the *inner*, *middle*, *outer* keywords.
+
+The single() function of 'uf3' pair style only return the 2-body interaction energy.
+
+Restrictions
+""""""""""""
+
+The 'uf3' pair style is part of the ML-UF3 package. It is only enabled if LAMMPS was built with that package. See the :doc:`Build package <Build_package>` page for more info.
+
+This pair style requires the :doc:`newton <newton>` setting to be "on".
+
+The UF3 LAMMPS potential files provided with LAMMPS (see the potentials directory) are parameterized for metal :doc:`units <units>`.
+
+
+Related commands
+""""""""""""""""
+
+:doc:`pair_coeff <pair_coeff>`
+
+Default
+"""""""
+
+none
+
+----------
+
+.. _Xie23:
+
+**(Xie23)** S. R. Xie, M. Rupp, and R. G. Hennig, "Ultra-fast interpretable machine-learning potentials", preprint arXiv:2110.00624v2 (2023)
diff --git a/lammps_plugin/scripts/README.rst b/lammps_plugin/scripts/README.rst
new file mode 100644
index 00000000..83975ff4
--- /dev/null
+++ b/lammps_plugin/scripts/README.rst
@@ -0,0 +1,15 @@
+This document describes how to generate UF3 lammps potential from an input structure file(`POSCAR <https://www.vasp.at/wiki/index.php/POSCAR>`_), and UF3 model file(JSON).
+
+Usuage-
+
+.. code:: bash
+
+   python generate_uf3_lammps_pots.py NAME_of_POSCAR NAME_of_UF3_model_File DIR_PATH_or_NAME 
+
+Will create :code:`DIR_PATH_or_NAME` if it does not exists and will write UF3 lammps potential files to :code:`DIR_PATH_or_NAME`. Make sure NAME_of_POSCAR and NAME_of_UF3_model_File are present in the same directory from which :code:`generate_uf3_lammps_pots.py` is executed.
+
+A lammps structure file with the name :code:`lammps.struct` will also be generated from the provided POSCAR file. This structure file should be used for running the lammps simulation. The potential files produced can also be used with other structures. But care must be taken to ensure that the elemental species map onto the right species number in the lammps structure file. 
+
+
+It will also additionally print lines that should be added to the lammps input file for using UF3 lammps potential files.
+ 
diff --git a/lammps_plugin/scripts/generate_uf3_lammps_pots.py b/lammps_plugin/scripts/generate_uf3_lammps_pots.py
new file mode 100644
index 00000000..fccfe012
--- /dev/null
+++ b/lammps_plugin/scripts/generate_uf3_lammps_pots.py
@@ -0,0 +1,219 @@
+from uf3.regression import least_squares
+from pymatgen.core.structure import Structure
+from pymatgen.core import Element
+from uf3.data import composition
+import numpy as np
+import os, sys
+
+if len(sys.argv) != 3:
+    #raise ValueError("Invalid number of arguments. Enter name of the structure\n\
+    #        file readable by pymatgen, UF3 model file, and name \n\
+    #        directory to write the UF3 potential files")
+    raise ValueError("Invalid number of arguments. Enter UF3 model file and name\n\
+            of directory to write the UF3 potential files")
+
+#struct = Structure.from_file(sys.argv[1])
+model = least_squares.WeightedLinearModel.from_json(sys.argv[1])
+
+#struct_elements = set(struct.symbol_set)
+model_elements = set(model.bspline_config.chemical_system.element_list)
+pot_dir = sys.argv[2]
+
+def create_element_map_for_lammps(uf3_chem_sys):
+    """Returns dict
+
+    Creates a mapping between element symbol and lammps element species,
+    to be used while creating UF3 lammps file. Takes UF3 composition object as
+    the input
+    """
+    lemap = {}
+    species_count = 1
+    for i in uf3_chem_sys.interactions_map[2]:
+        if i[0]==i[1]:
+            if i[0] not in lemap:
+                lemap[i[0]] = species_count
+                species_count += 1
+            else:
+                pass
+        else:
+            for j in i:
+                if j not in lemap:
+                    lemap[j] = species_count
+                    species_count += 1
+                else:
+                    pass
+    return lemap
+"""
+def write_lammps_ip_struct(struct_obj,element_map):
+    """#Returns str
+
+    #Converts a POSCAR to lammps structure format and writes 'lammps.struct' file
+    #Takes pymatgen structure object as the input
+"""
+    struct_dict = struct_obj.as_dict()
+
+    w_filename = 'lammps.struct'
+
+    lx = struct_obj.lattice.a
+    xy = struct_obj.lattice.b * np.cos(struct_obj.lattice.angles[2]*np.pi/180)
+    xz = struct_obj.lattice.c * np.cos(struct_obj.lattice.angles[1]*np.pi/180)
+
+    ly = np.sqrt(np.square(struct_obj.lattice.b) - np.square(xy))
+
+    yz = ((struct_obj.lattice.b*struct_obj.lattice.c*
+           np.cos(struct_obj.lattice.angles[0]*np.pi/180))-(xy*xz))/(ly)
+
+    lz = np.sqrt(np.square(struct_obj.lattice.c)-np.square(xz)-np.square(yz))
+
+    new_H = np.array([[lx,0,0],
+         [xy,ly,0],
+          [xz,yz,lz]])
+
+    fp = open(w_filename,'w')
+    fp.write("# Converted by ACH the great\n\n")
+    fp.write("%i atoms\n"%struct.num_sites)
+    fp.write("%i atom types\n\n"%len(struct.symbol_set))
+
+    fp.write("0.000000%.6f   xlo xhi\n"%lx)
+    fp.write("0.000000%.6f   ylo yhi\n"%ly)
+    fp.write("0.000000%.6f   zlo zhi\n\n"%lz)
+    fp.write("  %.6f%.6f%.6f   xy xz yz\n\n"%(xy,xz,yz))
+
+    fp.write("Masses\n\n")
+
+    for key in element_map:
+        at_mass = Element(key)
+        fp.write("  %i %.5f\n"%(element_map[key],at_mass.atomic_mass))
+        fp.write("\n")
+
+        fp.write("Atoms\n\n")
+    for i in range(0,struct.num_sites):
+        key = struct_dict['sites'][i]['species'][0]['element']
+        temp_cord = np.matmul(struct_dict['sites'][i]['abc'],new_H)
+        fp.write("   %i  %i %.6f %.6f %.6f\n"%(i+1,element_map[key],temp_cord[0],temp_cord[1],temp_cord[2]))
+    fp.close()
+    return w_filename
+"""
+def write_uf3_lammps_pot_files(chemical_sys,model,pot_dir):
+    """Returns list
+
+    Creates and writes UF3 lammps potential files. Takes UF3 composition object,
+    UF3 model and name of potential directory as input. Will overwrite the files
+    if files with the same exists
+    """
+    overwrite = True
+    if not os.path.exists(pot_dir):
+        os.mkdir(pot_dir)
+    files = {}
+
+    for interaction in chemical_sys.interactions_map[2]:
+        key = '_'.join(interaction)
+        files[key] = "#UF3 POT\n"
+        files[key] += "2B %i %i\n"%(model.bspline_config.leading_trim,model.bspline_config.trailing_trim)
+
+        files[key] += str(model.bspline_config.r_max_map[interaction]) + " " + \
+                str(len(model.bspline_config.knots_map[interaction]))+"\n"
+
+        files[key] += " ".join(['{:.17g}'.format(v) for v in \
+                model.bspline_config.knots_map[interaction]]) + "\n"
+
+        files[key] += str(model.bspline_config.get_interaction_partitions()[0][interaction]) \
+                + "\n"
+
+        start_index = model.bspline_config.get_interaction_partitions()[1][interaction]
+        length = model.bspline_config.get_interaction_partitions()[0][interaction]
+        files[key] += " ".join(['{:.17g}'.format(v) for v in \
+                model.coefficients[start_index:start_index + length]]) + "\n"
+
+        files[key] += "#"
+
+    if 3 in model.bspline_config.interactions_map:
+        for interaction in model.bspline_config.interactions_map[3]:
+            key = '_'.join(interaction)
+            files[key] = "#UF3 POT\n"
+            files[key] += "3B %i %i\n"%(model.bspline_config.leading_trim,model.bspline_config.trailing_trim)
+
+            files[key] += str(model.bspline_config.r_max_map[interaction][2]) \
+                    + " " + str(model.bspline_config.r_max_map[interaction][1]) \
+                    + " " + str(model.bspline_config.r_max_map[interaction][0]) + " "
+
+            files[key] += str(len(model.bspline_config.knots_map[interaction][2])) \
+                    + " " + str(len(model.bspline_config.knots_map[interaction][1])) + " " + str(len(model.bspline_config.knots_map[interaction][0])) + "\n"
+            files[key] += " ".join(['{:.17g}'.format(v) for v in \
+                    model.bspline_config.knots_map[interaction][2]]) + "\n"
+
+            files[key] += " ".join(['{:.17g}'.format(v) for v in \
+                    model.bspline_config.knots_map[interaction][1]]) + "\n"
+
+            files[key] += " ".join(['{:.17g}'.format(v) for v in \
+                    model.bspline_config.knots_map[interaction][0]]) + "\n"
+
+            solutions = least_squares.arrange_coefficients(model.coefficients, \
+                    model.bspline_config)
+
+            decompressed = model.bspline_config.decompress_3B( \
+                    solutions[(interaction[0], interaction[1],interaction[2])], \
+                    (interaction[0], interaction[1],interaction[2]))
+
+            files[key] += str(decompressed.shape[0]) + " " \
+                    + str(decompressed.shape[1]) + " " \
+                    + str(decompressed.shape[2]) + "\n"
+
+            for i in range(decompressed.shape[0]):
+                for j in range(decompressed.shape[1]):
+                    files[key] += ' '.join(map(str, decompressed[i,j]))
+                    files[key] += "\n"
+                    
+            files[key] += "#"
+
+    for k, v in files.items():
+        if not overwrite and os.path.exists(pot_dir + k):
+            continue
+        with open(pot_dir +"/"+ k, "w") as f:
+            f.write(v)
+    return files.keys()
+
+"""
+if len(model_elements.intersection(struct_elements))==len(struct_elements):
+    chemical_sys = composition.ChemicalSystem(element_list=list(struct_elements),\
+            degree=model.bspline_config.degree)
+    element_map = create_element_map_for_lammps(chemical_sys)
+    print(element_map)
+    write_lammps_ip_struct(struct_obj=struct,element_map=element_map)
+    pot_files = write_uf3_lammps_pot_files(chemical_sys=chemical_sys,model=model,pot_dir=pot_dir)
+    pot_files = list(pot_files)
+    lines = "pair_style uf3 %i %i\n"%(model.bspline_config.degree,len(struct_elements))
+
+    for interaction in chemical_sys.interactions_map[2]:
+        lines += "pair_coeff %i %i %s/%s\n"%(element_map[interaction[0]],
+                element_map[interaction[1]],pot_dir,'_'.join(interaction))
+
+    if 3 in model.bspline_config.interactions_map:
+        for interaction in model.bspline_config.interactions_map[3]:
+            lines += "pair_coeff %i %i %i %s/%s\n"%(element_map[interaction[0]],
+                    element_map[interaction[1]],element_map[interaction[2]],
+                    pot_dir,'_'.join(interaction))
+
+    print("***Add the following lines to the lammps input script***")
+    print(lines)
+else:
+    raise RuntimeError("Elements in the provided structure file does not match the elements in ")
+"""
+
+chemical_sys = model.bspline_config.chemical_system
+
+pot_files = write_uf3_lammps_pot_files(chemical_sys=chemical_sys,model=model,pot_dir=pot_dir)
+pot_files = list(pot_files)
+lines = "pair_style uf3 %i %i"%(model.bspline_config.degree,len(chemical_sys.element_list))
+
+"""
+for interaction in chemical_sys.interactions_map[2]:
+    lines += " %s/%s"%(pot_dir,'_'.join(interaction))
+
+if 3 in model.bspline_config.interactions_map:
+    for interaction in model.bspline_config.interactions_map[3]:
+        lines += " %s/%s"%(pot_dir,'_'.join(interaction))
+"""
+print("\n\n***Add the following line to the lammps input script followed by the 'pair_coeff' line/s***\n\n")
+print(lines)
+print("\n\n")
diff --git a/lammps_plugin/timing_info/README.rst b/lammps_plugin/timing_info/README.rst
new file mode 100644
index 00000000..e8d1cd77
--- /dev/null
+++ b/lammps_plugin/timing_info/README.rst
@@ -0,0 +1,47 @@
+This document describes the timing performance of lammps implementation of UF3.
+
+The tests were performed using a not yet published UF3 potential for the Nb-Sn system. The simulation cell consists of an Nb3Sn structure downloaded from the `Materials Project <materialsproject.org/>`_. The lammps simulations were carried out on the `Hipegator machine <https://help.rc.ufl.edu/doc/Available_Node_Features>`_ using a single AMD EPYC 7702 64-Core processor.
+
+For benchmarking the implementation, the perfomance is compared to simulations performed with LJ potentials. For this comparison a LJ potnetial was fitted to 2-body part of the UF3 potential.
+
+All the simulations utilized perodic boundary conditions and the `NVE fix <https://docs.lammps.org/fix_nve.html>`_. The atoms were first given initial velocities corresponding to 3000K (Gaussian distribution). All the simulations were terminated after 10000 steps.
+
+In the following tables performance of the current implementation is given.
+
+.. list-table:: Nb3Sn 2x2x2 supercell- 64 atoms, 1fs timestep, 1 core
+    :header-rows: 1
+
+    * - Potential
+      - Speed (timestep/s)
+    * - LJ
+      - 6513.272
+    * - UF2
+      - 1163.371
+    * - UF3
+      - 226.488
+
+.. list-table:: Nb3Sn 4x4x4 supercell- 512 atoms, 1fs timestep, 1 core
+    :header-rows: 1
+
+    * - Potential
+      - Speed (timestep/s)
+    * - LJ
+      - 960.561
+    * - UF2
+      - 115.658
+    * - UF3
+      - 28.214
+
+
+.. list-table:: Performance of UF3 with system size, 1fs timestep, 1 core
+    :header-rows: 1
+
+    * - Supercell size
+      - Speed (timestep/s)
+    * - 2x2x2
+      - 226.488
+    * - 3x3x3
+      - 63.358
+    * - 4x4x4
+      - 28.214
+
diff --git a/lammps_plugin/tungsten_example/W_W b/lammps_plugin/tungsten_example/W_W
new file mode 100644
index 00000000..01d4b971
--- /dev/null
+++ b/lammps_plugin/tungsten_example/W_W
@@ -0,0 +1,7 @@
+#UF3 POT
+2B 0 3
+5.5 22
+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
+18
+85.256465889606673 85.256465889606673 60.118514370108457 39.372994322346095 23.019905746319598 11.095292589050336 2.5770289203919474 -0.21591409154943711 -0.41582603147952274 -0.29361205978352245 -0.23467230507950282 -0.13083777191249607 -0.031024404500506326 -0.059144784016542738 -0.037145211195670137 0 0 0
+#
\ No newline at end of file
diff --git a/lammps_plugin/tungsten_example/W_W_W b/lammps_plugin/tungsten_example/W_W_W
new file mode 100644
index 00000000..3e71f65d
--- /dev/null
+++ b/lammps_plugin/tungsten_example/W_W_W
@@ -0,0 +1,89 @@
+#UF3 POT
+3B 0 3
+7.0 3.5 3.5 19 13 13
+1.5 1.5 1.5 1.5 1.9583333333333333 2.4166666666666665 2.875 3.333333333333333 3.7916666666666665 4.25 4.708333333333333 5.1666666666666661 5.625 6.083333333333333 6.5416666666666661 7 7 7 7
+1.5 1.5 1.5 1.5 1.8333333333333333 2.1666666666666665 2.5 2.833333333333333 3.1666666666666665 3.5 3.5 3.5 3.5
+1.5 1.5 1.5 1.5 1.8333333333333333 2.1666666666666665 2.5 2.833333333333333 3.1666666666666665 3.5 3.5 3.5 3.5
+9 9 15
+-1.0025267383014886e-06 1.627666892216261e-10 8.48457178601592e-09 1.0410120456668556e-08 1.1330088302187643e-06 1.4224457260740182e-05 3.749471162873949e-05 3.695233634472074e-05 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+4.922778833696604e-05 -1.0222423928319826e-06 -4.8075020979943344e-05 -3.8518708184701e-05 -9.794643241404816e-06 -5.109554485113297e-07 -4.104591531909916e-10 5.419385825688335e-09 2.3261619203718816e-06 0.0 0.0 0.0 0.0 0.0 0.0
+2.829068261701631e-05 3.3565662462635585e-05 -4.4282774329572976e-05 5.3019338472701156e-05 1.656981200092865e-05 5.764723793849365e-06 8.90212763080714e-06 -1.266718085517006e-06 2.431925564724538e-06 -3.3416793274446366e-07 0.0 0.0 0.0 0.0 0.0
+9.661425891180097e-09 1.1180525909838732e-07 7.372128540386131e-06 9.39083445715028e-06 1.0838653682096786e-05 9.806075321172179e-05 7.23071321840456e-05 0.00011393368280907902 7.532201192649856e-05 1.4915718900782886e-05 0.0 0.0 0.0 0.0 0.0
+4.674109331563019e-07 -5.406616138544008e-09 2.0217516565542138e-08 5.815480232140697e-06 2.377543599749233e-05 4.214830983861595e-05 -7.397778428926203e-05 7.061700585974997e-05 3.398811616466058e-05 3.0440209444922143e-05 0.00013514584901612643 0.0 0.0 0.0 0.0
+-2.18344999249825e-05 -0.00011697904895508256 -5.553581955175406e-06 2.658838454117769e-08 1.2146153008512734e-05 1.918298817375017e-05 5.40154814815751e-05 1.2944187401636162e-05 3.710536690671586e-05 1.6795050072301748e-05 -7.229556489685409e-06 8.543467581763247e-06 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+4.922778833696604e-05 -1.0222423928319826e-06 -4.8075020979943344e-05 -3.8518708184701e-05 -9.794643241404816e-06 -5.109554485113297e-07 -4.104591531909916e-10 5.419385825688335e-09 2.3261619203718816e-06 0.0 0.0 0.0 0.0 0.0 0.0
+4.7777593982974866e-05 2.1830415880551608e-05 -8.74161436707018e-06 7.718518436280754e-09 3.872620009304499e-06 2.2808275890118977e-05 4.9693195966817054e-05 -5.792733170014431e-06 4.574471594251372e-05 -2.840624725005251e-05 0.0 0.0 0.0 0.0 0.0
+-2.983543176137801e-05 4.9033674735689904e-05 4.195810226183578e-05 0.00010403045298571696 0.0003809288520722337 3.373314734861335e-05 1.0875278462846613e-05 2.680082116524925e-05 2.6236988936225364e-05 3.6338704048890435e-06 0.0 0.0 0.0 0.0 0.0
+2.1399512390776296e-05 2.7931318963336216e-05 1.601611374550643e-05 0.0011462578563508437 0.0020096052304490785 0.0010339536769749804 0.001014125101157229 0.00021102390605096824 -9.581571365103867e-06 -2.280019726143644e-05 3.57351467975471e-07 0.0 0.0 0.0 0.0
+6.2467013661098965e-06 2.658299496065299e-05 0.0007087760404542255 -0.002278488280232273 -0.007044271575650293 0.0038614975616579586 0.009260853457130498 -0.0010473278276317726 -8.516349321648054e-05 0.00020262630262751473 1.3414729868362445e-06 7.953262907187799e-08 0.0 0.0 0.0
+2.3603382351924323e-06 7.146280072503116e-05 0.0012361202483760594 -0.005734465078354871 -0.017135562936976125 -0.0016664758722632735 0.010195949308015286 -0.005412690557758046 -0.0029107216230642853 -0.0007612908343076258 -0.00013449226099999524 2.1468761472671134e-05 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+2.829068261701631e-05 3.3565662462635585e-05 -4.4282774329572976e-05 5.3019338472701156e-05 1.656981200092865e-05 5.764723793849365e-06 8.90212763080714e-06 -1.266718085517006e-06 2.431925564724538e-06 -3.3416793274446366e-07 0.0 0.0 0.0 0.0 0.0
+-2.983543176137801e-05 4.9033674735689904e-05 4.195810226183578e-05 0.00010403045298571696 0.0003809288520722337 3.373314734861335e-05 1.0875278462846613e-05 2.680082116524925e-05 2.6236988936225364e-05 3.6338704048890435e-06 0.0 0.0 0.0 0.0 0.0
+4.805582975622838e-05 3.5904852202535306e-05 0.0001861581391418004 0.011742061168666256 0.024776619721352023 0.009905830782079064 0.001165275555230131 0.0004071734306161239 0.00037033988296356686 4.258576293935207e-05 -1.0503156955159468e-06 0.0 0.0 0.0 0.0
+-5.680777858210279e-07 7.084011901144194e-05 0.0050124193537796876 0.1029510036679177 0.24007141064120838 0.17886074568759117 0.03778564255583381 -0.01538608134214139 -0.01488601265585909 -0.005013137591006775 1.477967512985198e-07 2.827576677896298e-05 0.0 0.0 0.0
+3.4701391129480595e-05 0.005680635772167662 0.07108357438929552 0.18745843365735923 0.18473391598372302 0.2461411417071668 -0.005111666129467331 -0.24038059036230178 -0.19258370027930355 -0.048958604017684836 -0.0025501911244683846 -4.982008439539256e-06 0.0 0.0 0.0
+1.8054960278101902e-05 0.011611823828717463 0.14139350699353365 0.1678085431015684 -0.22152577092526715 0.26507903411379746 0.08279760718466518 0.02879361559120084 0.04448600527937973 -0.006339073225905362 0.012182164968354765 -0.0016647733469859114 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+9.661425891180097e-09 1.1180525909838732e-07 7.372128540386131e-06 9.39083445715028e-06 1.0838653682096786e-05 9.806075321172179e-05 7.23071321840456e-05 0.00011393368280907902 7.532201192649856e-05 1.4915718900782886e-05 0.0 0.0 0.0 0.0 0.0
+2.1399512390776296e-05 2.7931318963336216e-05 1.601611374550643e-05 0.0011462578563508437 0.0020096052304490785 0.0010339536769749804 0.001014125101157229 0.00021102390605096824 -9.581571365103867e-06 -2.280019726143644e-05 3.57351467975471e-07 0.0 0.0 0.0 0.0
+-5.680777858210279e-07 7.084011901144194e-05 0.0050124193537796876 0.1029510036679177 0.24007141064120838 0.17886074568759117 0.03778564255583381 -0.01538608134214139 -0.01488601265585909 -0.005013137591006775 1.477967512985198e-07 2.827576677896298e-05 0.0 0.0 0.0
+3.224364873093737e-05 0.0005926077146340155 0.03682711928362166 0.046898795979650645 0.054125869772576704 0.49020499621758107 0.3614418794239106 0.569703720925391 0.37654073726442117 0.07453613145477757 0.002320760953344357 -4.1320810202123745e-05 0.0 0.0 0.0
+7.181321030642065e-05 0.029059269505937943 0.11884075783303522 0.21076878762972284 -0.36996025259611676 0.35312700844142586 0.16990097184555128 0.15210775765373416 0.6757140775652752 -0.10922320044264704 -0.5849697059462919 -0.027759887396949002 0.0 0.0 0.0
+9.734779590447442e-05 0.06072685888500135 0.09586373042400446 0.2700889438647233 0.06469057221256519 0.1854847883254132 0.08392153868089386 -0.03626027544891736 0.04232832520121794 0.23849583983886427 0.10911730115168922 -0.04380214407208333 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+4.674109331563019e-07 -5.406616138544008e-09 2.0217516565542138e-08 5.815480232140697e-06 2.377543599749233e-05 4.214830983861595e-05 -7.397778428926203e-05 7.061700585974997e-05 3.398811616466058e-05 3.0440209444922143e-05 0.00013514584901612643 0.0 0.0 0.0 0.0
+6.2467013661098965e-06 2.658299496065299e-05 0.0007087760404542255 -0.002278488280232273 -0.007044271575650293 0.0038614975616579586 0.009260853457130498 -0.0010473278276317726 -8.516349321648054e-05 0.00020262630262751473 1.3414729868362445e-06 7.953262907187799e-08 0.0 0.0 0.0
+3.4701391129480595e-05 0.005680635772167662 0.07108357438929552 0.18745843365735923 0.18473391598372302 0.2461411417071668 -0.005111666129467331 -0.24038059036230178 -0.19258370027930355 -0.048958604017684836 -0.0025501911244683846 -4.982008439539256e-06 0.0 0.0 0.0
+7.181321030642065e-05 0.029059269505937943 0.11884075783303522 0.21076878762972284 -0.36996025259611676 0.35312700844142586 0.16990097184555128 0.15210775765373416 0.6757140775652752 -0.10922320044264704 -0.5849697059462919 -0.027759887396949002 0.0 0.0 0.0
+3.8792902483003914e-05 0.01940944659869733 0.11401272918676479 0.24845559032094336 -0.029040137137018487 0.22760022061135138 -0.14408030764137164 -0.1502303415258775 0.24414290246285578 0.20048840049956973 -0.3098919852331769 0.38455675655513366 0.0 0.0 0.0
+0.0001449154156948679 0.05388868537335966 0.1340445791122114 0.1304729995514306 0.02041924003827766 0.11401975265164792 0.13440972558465583 -0.056806908474571434 0.03945347719095876 0.02504259931146615 0.04957547514580214 0.01854515348232565 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+-2.18344999249825e-05 -0.00011697904895508256 -5.553581955175406e-06 2.658838454117769e-08 1.2146153008512734e-05 1.918298817375017e-05 5.40154814815751e-05 1.2944187401636162e-05 3.710536690671586e-05 1.6795050072301748e-05 -7.229556489685409e-06 8.543467581763247e-06 0.0 0.0 0.0
+2.3603382351924323e-06 7.146280072503116e-05 0.0012361202483760594 -0.005734465078354871 -0.017135562936976125 -0.0016664758722632735 0.010195949308015286 -0.005412690557758046 -0.0029107216230642853 -0.0007612908343076258 -0.00013449226099999524 2.1468761472671134e-05 0.0 0.0 0.0
+1.8054960278101902e-05 0.011611823828717463 0.14139350699353365 0.1678085431015684 -0.22152577092526715 0.26507903411379746 0.08279760718466518 0.02879361559120084 0.04448600527937973 -0.006339073225905362 0.012182164968354765 -0.0016647733469859114 0.0 0.0 0.0
+9.734779590447442e-05 0.06072685888500135 0.09586373042400446 0.2700889438647233 0.06469057221256519 0.1854847883254132 0.08392153868089386 -0.03626027544891736 0.04232832520121794 0.23849583983886427 0.10911730115168922 -0.04380214407208333 0.0 0.0 0.0
+0.0001449154156948679 0.05388868537335966 0.1340445791122114 0.1304729995514306 0.02041924003827766 0.11401975265164792 0.13440972558465583 -0.056806908474571434 0.03945347719095876 0.02504259931146615 0.04957547514580214 0.01854515348232565 0.0 0.0 0.0
+5.588250079169304e-05 -0.012017079745903034 -0.0008295411205340048 0.00043274236043781584 0.036933605134660964 -0.017759680455812197 0.15248826375477179 -0.022974467887332546 0.005116007779072725 -0.01249692329646853 -0.12288470938617652 0.022731129576009593 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+#
\ No newline at end of file
diff --git a/lammps_plugin/tungsten_example/W_nve.in b/lammps_plugin/tungsten_example/W_nve.in
new file mode 100644
index 00000000..32071a00
--- /dev/null
+++ b/lammps_plugin/tungsten_example/W_nve.in
@@ -0,0 +1,31 @@
+# # ============= basic parameters
+units           metal
+boundary        p p p
+atom_style      atomic
+
+# # Define Variables
+variable        T equal 3000
+
+# # ============= import structure file
+read_data       struct.lammps
+mass            1 183.84
+neighbor        0.0 bin
+neigh_modify    every 5 delay 0
+# ============= force field
+pair_style      zero 10
+pair_style      uf3 3 1
+pair_coeff      1 1 W_W
+pair_coeff      1 1 W_W_W
+
+# ============= give all atoms velocity
+timestep        0.001
+velocity        all create $T 4928459 dist gaussian
+
+# ============= Equilibrate the entire model
+fix             1  all nve
+run             10000
+
+#######################################
+#
+## SIMULATION DONE
+print "All done"
diff --git a/lammps_plugin/tungsten_example/relax_lammps_2b.in b/lammps_plugin/tungsten_example/relax_lammps_2b.in
new file mode 100644
index 00000000..c68b9cd9
--- /dev/null
+++ b/lammps_plugin/tungsten_example/relax_lammps_2b.in
@@ -0,0 +1,26 @@
+# # ============= basic parameters
+units 		metal
+boundary 	p p p
+atom_style 	atomic
+
+# # ============= import structure file
+read_data	struct.lammps
+mass		1 183.84
+neighbor	0.0 bin
+# ============= force field
+pair_style  zero 10.0
+pair_style  uf3 2 1
+pair_coeff  1 1 W_W
+
+#----------------------Run Minimization-------------------------
+reset_timestep	0
+fix 1 all box/relax iso 0
+thermo 10
+thermo_style custom step pe lx ly lz press pxx pyy pzz  
+min_style cg
+minimize 1e-15 1e-15 5000 5000
+#
+#######################################
+#
+## SIMULATION DONE
+print "All done"
diff --git a/lammps_plugin/tungsten_example/relax_lammps_3b.in b/lammps_plugin/tungsten_example/relax_lammps_3b.in
new file mode 100644
index 00000000..1be90759
--- /dev/null
+++ b/lammps_plugin/tungsten_example/relax_lammps_3b.in
@@ -0,0 +1,27 @@
+# # ============= basic parameters
+units 		metal
+boundary 	p p p
+atom_style 	atomic
+
+# # ============= import structure file
+read_data	struct.lammps
+mass		1 183.84
+neighbor	0.0 bin
+# ============= force field
+pair_style  zero 10.0
+pair_style  uf3 3 1
+pair_coeff  1 1 W_W
+pair_coeff  1 1 1 W_W_W
+
+#----------------------Run Minimization-------------------------
+reset_timestep	0
+fix 1 all box/relax iso 0
+thermo 10
+thermo_style custom step pe lx ly lz press pxx pyy pzz  
+min_style cg
+minimize 1e-15 1e-15 5000 5000
+#
+#######################################
+#
+## SIMULATION DONE
+print "All done"
diff --git a/lammps_plugin/tungsten_example/struct.lammps b/lammps_plugin/tungsten_example/struct.lammps
new file mode 100644
index 00000000..c72ef718
--- /dev/null
+++ b/lammps_plugin/tungsten_example/struct.lammps
@@ -0,0 +1,67 @@
+# Converted from POSCAR to lammps format
+
+54 atoms
+1 atom types
+
+0.000000    9.510947   xlo xhi
+0.000000    9.510947   ylo yhi
+0.000000    9.510947   zlo zhi
+
+  0.000000    0.000000    0.000000   xy xz yz
+
+Atoms
+
+   1  1     0.000000 0.000000 0.000000
+   2  1     0.000000 0.000000 3.170312
+   3  1     0.000000 0.000000 6.340635
+   4  1     9.510947 3.170312 0.000000
+   5  1     9.510947 3.170312 3.170312
+   6  1     9.510947 3.170312 6.340635
+   7  1     9.510947 6.340635 0.000000
+   8  1     9.510947 6.340635 3.170312
+   9  1     9.510947 6.340635 6.340635
+  10  1     3.170312 0.000000 0.000000
+  11  1     3.170312 0.000000 3.170312
+  12  1     3.170312 0.000000 6.340635
+  13  1     3.170312 3.170312 0.000000
+  14  1     3.170312 3.170312 3.170312
+  15  1     3.170312 3.170312 6.340635
+  16  1     3.170312 6.340635 0.000000
+  17  1     3.170312 6.340635 3.170312
+  18  1     3.170312 6.340635 6.340635
+  19  1     6.340635 0.000000 0.000000
+  20  1     6.340635 0.000000 3.170312
+  21  1     6.340635 0.000000 6.340635
+  22  1     6.340635 3.170312 0.000000
+  23  1     6.340635 3.170312 3.170312
+  24  1     6.340635 3.170312 6.340635
+  25  1     6.340635 6.340635 0.000000
+  26  1     6.340635 6.340635 3.170312
+  27  1     6.340635 6.340635 6.340635
+  28  1     1.585161 1.585161 1.585161
+  29  1     1.585161 1.585161 4.755473
+  30  1     1.585161 1.585161 7.925786
+  31  1     1.585161 4.755474 1.585161
+  32  1     1.585161 4.755474 4.755473
+  33  1     1.585161 4.755474 7.925786
+  34  1     1.585161 7.925786 1.585161
+  35  1     1.585161 7.925786 4.755473
+  36  1     1.585161 7.925786 7.925786
+  37  1     4.755474 1.585161 1.585161
+  38  1     4.755474 1.585161 4.755473
+  39  1     4.755474 1.585161 7.925786
+  40  1     4.755474 4.755474 1.585161
+  41  1     4.755474 4.755474 4.755473
+  42  1     4.755474 4.755474 7.925786
+  43  1     4.755474 7.925786 1.585161
+  44  1     4.755474 7.925786 4.755473
+  45  1     4.755474 7.925786 7.925786
+  46  1     7.925786 1.585161 1.585161
+  47  1     7.925786 1.585161 4.755473
+  48  1     7.925786 1.585161 7.925786
+  49  1     7.925786 4.755474 1.585161
+  50  1     7.925786 4.755474 4.755473
+  51  1     7.925786 4.755474 7.925786
+  52  1     7.925786 7.925786 1.585161
+  53  1     7.925786 7.925786 4.755473
+  54  1     7.925786 7.925786 7.925786


Test patch:
diff --git a/lammps_plugin/unittest/A_A.uf3_pot b/lammps_plugin/unittest/A_A.uf3_pot
new file mode 100644
index 00000000..cb9371d8
--- /dev/null
+++ b/lammps_plugin/unittest/A_A.uf3_pot
@@ -0,0 +1,7 @@
+#UF3 POT
+2B 0 3
+8.0 31
+0.001 0.001 0.001 0.001 0.33429166666666665 0.66758333333333331 1.000875 1.3341666666666665 1.6674583333333333 2.00075 2.3340416666666663 2.6673333333333331 3.0006249999999999 3.3339166666666666 3.667208333333333 4.0004999999999997 4.3337916666666665 4.6670833333333333 5.000375 5.3336666666666668 5.6669583333333335 6.0002500000000003 6.3335416666666671 6.6668333333333338 7.0001249999999997 7.3334166666666665 7.6667083333333332 8 8 8 8
+27
+79.140244588519465 79.140244588519465 55.85833391113556 36.597903318706138 21.358952811231141 12.290000872768841 1.9593931914091953 -0.65697974623243804 -0.85177956270573463 -0.68929688239869991 -0.46787243412973262 -0.27624655899523165 -0.11912921944351409 -0.056302369393035338 -0.0049812809608429064 0.0085637634684603507 0.0034716161454604712 -0.0058751075573311978 -0.005453415412748467 -0.0015123194244718201 0.0011577919587182201 0.001583772506713282 -0.00049823976100720228 -0.0013902809146717273 0 0 0
+#
\ No newline at end of file
diff --git a/lammps_plugin/unittest/A_A_A.uf3_pot b/lammps_plugin/unittest/A_A_A.uf3_pot
new file mode 100644
index 00000000..1f089fda
--- /dev/null
+++ b/lammps_plugin/unittest/A_A_A.uf3_pot
@@ -0,0 +1,129 @@
+#UF3 POT
+3B 0 3
+8.0 4.0 4.0 23 15 15
+0.001 0.001 0.001 0.001 0.50093749999999992 1.000875 1.5008124999999999 2.00075 2.5006874999999997 3.0006249999999999 3.5005624999999996 4.0004999999999997 4.5004375000000003 5.000375 5.5003124999999997 6.0002500000000003 6.5001875 7.0001249999999997 7.5000625000000003 8 8 8 8
+0.001 0.001 0.001 0.001 0.50087499999999996 1.00075 1.5006249999999999 2.0005000000000002 2.500375 3.0002499999999999 3.5001250000000002 4 4 4 4
+0.001 0.001 0.001 0.001 0.50087499999999996 1.00075 1.5006249999999999 2.0005000000000002 2.500375 3.0002499999999999 3.5001250000000002 4 4 4 4
+11 11 19
+-1.1790416072105636e-06 6.589114265858035e-08 2.1094970065385374e-06 4.7014910818419987e-07 8.288423734406254e-06 0.000186151370764668 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+3.2366723157260956e-05 2.6208783380066457e-05 1.0239981836366566e-05 1.8487993936404763e-05 2.1943710009352506e-05 2.6899947783571087e-07 5.002786118380638e-06 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+4.0696915445106504e-05 6.184559535738335e-05 5.870203846064511e-05 -1.112085789465658e-05 -4.7600144351359347e-07 -5.861397594145908e-08 1.4524208784805573e-08 2.300649782987421e-06 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+5.282737863089323e-05 6.328946161646202e-06 1.8329325276370316e-05 1.6423572667388823e-05 1.2653184610977003e-06 7.181714140248046e-06 3.491501462345434e-06 -7.285463619241614e-06 3.2609159022388403e-06 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+6.9413198850914024e-06 2.107626397843018e-05 1.8155172114721186e-05 2.0928626557075606e-06 1.5632037328512312e-06 -2.7335717313450097e-07 -7.2126792356200426e-09 9.213093725547886e-09 4.186629643010996e-08 8.198811769753182e-08 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+1.4231393739579515e-06 3.630746449160232e-07 9.318604659023228e-07 4.92311430374376e-07 -3.701479331898353e-09 2.1280257031614452e-07 1.2240989510544568e-06 5.3432540178806065e-06 2.043230389835189e-06 3.2740024159475547e-07 6.717304982644579e-07 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 5.196157077391817e-07 6.938124100654148e-08 7.597564197383545e-08 1.863740632660483e-07 4.437837629589167e-07 5.453941063185757e-07 1.5602917821833568e-06 3.404289212094662e-07 9.967447994956849e-07 5.8845599651090215e-06 1.5052240335012455e-05 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 3.196534127251379e-06 -1.1872677866681306e-06 -2.5678892066098854e-08 5.139117071436217e-09 1.1142431390092631e-06 2.0605776537608227e-06 5.297265009242829e-06 6.713907186544732e-06 2.7028644452395994e-06 1.149242737988068e-06 2.2449682976927855e-06 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+3.2366723157260956e-05 2.6208783380066457e-05 1.0239981836366566e-05 1.8487993936404763e-05 2.1943710009352506e-05 2.6899947783571087e-07 5.002786118380638e-06 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+9.36844864368356e-06 1.9223644807676324e-05 1.9979026246524356e-05 3.627062812549574e-05 9.775578281629195e-06 -5.894357549683858e-06 6.470814473155067e-07 2.31805322174729e-06 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+2.2218614171918013e-06 5.325319655352672e-06 7.766746363269582e-06 9.361315506075464e-06 5.0417710282874456e-05 9.822946186678772e-05 0.00026400568406806884 0.00033610865151919737 0.00013239814531221768 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+4.4850777249143735e-05 7.094600012126306e-05 0.00030581781354430576 0.00044661036994300023 0.00016699596636619577 1.5860625743775105e-05 9.74250537001798e-07 5.385650613476577e-06 8.091278451728344e-06 1.2460869401480828e-05 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0001028691918645833 0.0002737040057685444 0.0003861446001781946 0.0004042287651515365 0.0017229200225725174 0.003198296698131205 0.008774096120579751 0.011237818178923189 0.004334800036723805 0.0007344916552783145 -0.0001506915192259342 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+-1.2482511756156149e-05 1.365124801275985e-05 1.3894049203809568e-05 2.3985465221727954e-05 3.3458449092465795e-05 0.00028172299406359233 0.00040056109827889085 0.0004621959325200118 0.0034637215474633033 0.009153352872912168 0.012804683731760212 0.010674833967812809 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.005346497059990333 -0.002010348201210142 -0.0010943235863089423 0.0001661513182702165 0.00012025969610516196 2.4949866002221845e-05 6.627236360802077e-06 3.003757825105864e-06 3.997348910159012e-05 0.000427961841918743 0.0007451357800599296 0.0011219432594133996 0.009685550613014016 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.013431011504370738 0.007883426617122005 0.007935899204760883 0.005880150773602205 0.0009832099103910489 0.005414528729313218 0.0015950126575825377 0.00024127039666882992 -1.5674461809944553e-05 1.0711548076574028e-05 7.830483572860064e-06 0.00011012649333888752 0.0005497452692208139 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+4.0696915445106504e-05 6.184559535738335e-05 5.870203846064511e-05 -1.112085789465658e-05 -4.7600144351359347e-07 -5.861397594145908e-08 1.4524208784805573e-08 2.300649782987421e-06 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+2.2218614171918013e-06 5.325319655352672e-06 7.766746363269582e-06 9.361315506075464e-06 5.0417710282874456e-05 9.822946186678772e-05 0.00026400568406806884 0.00033610865151919737 0.00013239814531221768 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.00013804169495254125 0.0012238504051117233 0.01461439973735456 0.010197731078827295 0.009003775355755566 0.030381456320656558 0.024785731678029766 0.004637019267552505 -0.005938106654005813 -0.002605150959220643 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+-0.0002216303488927365 0.0001541319391627563 4.83626397765333e-05 2.1041902272582753e-05 0.00026610797279588076 0.00045665788403242036 0.00017325291338578903 0.0035336618936866277 0.018540440861910777 0.003501320637152642 0.002219074201926699 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.01652594711887213 0.011725805945708163 0.011173124362203699 0.0024178633067081135 -0.01796173780303683 -0.011618890946870497 -0.0008528234196397706 -9.706084806556783e-05 -2.754043401157181e-06 0.0001324948483342069 5.482811058752758e-05 0.0003886706609323921 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.008855607599180339 0.015278371659152929 0.0023879884417463693 0.000935102721182137 0.0034651523786104546 0.005282714096158778 0.0012741283037854573 0.006268847413803995 0.004013755514818873 0.010363477891131097 0.007109323912817858 -0.0002600093944564617 -0.00034633976332068713 0.0 0.0 0.0 0.0 0.0 0.0
+-2.698386662730078e-05 2.754973422173369e-07 0.0001320534807487939 0.004372176148977807 0.0015642026255259442 0.004248632573013906 0.00040885420395593786 0.001088966135412402 0.004766232525411325 0.002205157762668968 0.007459939889093756 0.005587608653898612 -0.004720162133268877 -0.006593174803103767 0.0 0.0 0.0 0.0 0.0
+-7.940872165606751e-05 2.276762148612182e-06 0.00010635762128769112 2.049233578255131e-05 0.00042145671490654473 0.01249692872936893 0.0020370917425772224 0.0017316344055948985 0.0006574004028558345 0.0012283310563930355 0.0014629455315045585 -4.492954039177435e-06 0.00029547002108771967 0.002683727758662211 0.004096710661285439 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+5.282737863089323e-05 6.328946161646202e-06 1.8329325276370316e-05 1.6423572667388823e-05 1.2653184610977003e-06 7.181714140248046e-06 3.491501462345434e-06 -7.285463619241614e-06 3.2609159022388403e-06 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+4.4850777249143735e-05 7.094600012126306e-05 0.00030581781354430576 0.00044661036994300023 0.00016699596636619577 1.5860625743775105e-05 9.74250537001798e-07 5.385650613476577e-06 8.091278451728344e-06 1.2460869401480828e-05 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+-0.0002216303488927365 0.0001541319391627563 4.83626397765333e-05 2.1041902272582753e-05 0.00026610797279588076 0.00045665788403242036 0.00017325291338578903 0.0035336618936866277 0.018540440861910777 0.003501320637152642 0.002219074201926699 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.003917873880791907 -0.0007874629498581528 -2.4595030318112164e-05 -2.8508297646329816e-06 1.8504666071760445e-07 0.0001190800388356091 0.0035373487148805376 0.00037674157183609377 0.0012087894330956167 0.0009651695201594091 4.769364472898923e-05 0.0004612591073953361 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.00022640297379380707 -0.0005111351623843819 0.0002024994190007784 0.00044921539785371963 0.001393999756415734 0.0011754659552919043 7.761533188706794e-05 5.8170561410888746e-05 -1.0918989217761552e-05 -2.9455190099531973e-07 2.0845708467284646e-07 1.0072947935068441e-06 2.449241542240889e-06 0.0 0.0 0.0 0.0 0.0 0.0
+5.3159155089314414e-05 1.2012438398909825e-05 3.274115317951001e-05 1.689353008824745e-05 -4.232110203859359e-07 7.656160171407207e-06 5.754938338062256e-05 0.0002687156220968384 9.768544317740195e-05 1.0129360414562531e-05 2.293272526112481e-05 1.8260890221186993e-05 1.7288534885724222e-06 1.8856324749638164e-06 0.0 0.0 0.0 0.0 0.0
+4.542360076931743e-06 1.2329326209575631e-05 1.7173803033436737e-05 5.0102066463061734e-05 1.09067765324765e-05 4.930240898900306e-05 0.00028721835291257015 0.0007503332386451459 0.00015238128535605624 -6.391111549761724e-05 -4.677072820313549e-06 -5.572104125200205e-06 2.707663268609677e-05 5.022846595129856e-05 0.00012970076559689836 0.0 0.0 0.0 0.0
+0.0001650881267658455 6.548635015912796e-05 3.282822556024051e-05 7.205996516588195e-05 0.00030710620226134084 0.0007279352507540159 0.000909918529220897 0.0017257383928761386 0.00047543969768972346 -0.00030462814537952123 2.7321207199326783e-05 0.00010130812246147248 5.4817489360932934e-05 0.00012640103175376577 0.00018995742264027741 0.00023929212126678798 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+6.9413198850914024e-06 2.107626397843018e-05 1.8155172114721186e-05 2.0928626557075606e-06 1.5632037328512312e-06 -2.7335717313450097e-07 -7.2126792356200426e-09 9.213093725547886e-09 4.186629643010996e-08 8.198811769753182e-08 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0001028691918645833 0.0002737040057685444 0.0003861446001781946 0.0004042287651515365 0.0017229200225725174 0.003198296698131205 0.008774096120579751 0.011237818178923189 0.004334800036723805 0.0007344916552783145 -0.0001506915192259342 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.01652594711887213 0.011725805945708163 0.011173124362203699 0.0024178633067081135 -0.01796173780303683 -0.011618890946870497 -0.0008528234196397706 -9.706084806556783e-05 -2.754043401157181e-06 0.0001324948483342069 5.482811058752758e-05 0.0003886706609323921 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.00022640297379380707 -0.0005111351623843819 0.0002024994190007784 0.00044921539785371963 0.001393999756415734 0.0011754659552919043 7.761533188706794e-05 5.8170561410888746e-05 -1.0918989217761552e-05 -2.9455190099531973e-07 2.0845708467284646e-07 1.0072947935068441e-06 2.449241542240889e-06 0.0 0.0 0.0 0.0 0.0 0.0
+0.0016181071043329578 0.003179470547568356 0.008661622548635572 0.011071785334468471 0.004294892778359652 0.0017845979744737465 0.0034643761195723064 0.015112039067322293 0.022192108732694595 0.008134230944897397 0.0007595380961610584 1.6727218309602107e-05 0.00012823915020345735 0.0001971442066043176 0.0 0.0 0.0 0.0 0.0
+0.00033854327480422193 0.0032940012133255356 0.008932075729876752 0.012661798131960687 0.013919159699477152 0.08208818801401566 0.15644219608737447 0.4341503084393359 0.557964838826116 0.21173866865770563 0.03478604116524652 -0.007972916161324952 -0.0007132020154210059 0.0004420932693293155 0.0003860741867263207 0.0 0.0 0.0 0.0
+0.0005943091729493132 0.0009672866630600067 0.009096489467732383 0.01304005169719466 0.01734885651081947 0.1652472792755658 0.44880261475702005 0.6329167707872334 0.5244926606398645 0.26270394857828266 -0.10487776273847933 -0.05845317239353218 0.0032135333670803676 0.003777214235598332 0.0007681993725802362 0.00011959437074006901 0.0 0.0 0.0
+6.489304732024981e-05 0.0011748014116558024 0.014162497237899634 0.024470256017495366 0.046651858009005745 0.46415330969555396 0.6536160620847673 0.38295932353650225 0.39318393240673155 0.2877811024545165 0.03747496490739291 0.2654106688863148 0.074634565935104 0.012330295300167044 -0.0004996457463809098 0.00023089219471653216 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+1.4231393739579515e-06 3.630746449160232e-07 9.318604659023228e-07 4.92311430374376e-07 -3.701479331898353e-09 2.1280257031614452e-07 1.2240989510544568e-06 5.3432540178806065e-06 2.043230389835189e-06 3.2740024159475547e-07 6.717304982644579e-07 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+-1.2482511756156149e-05 1.365124801275985e-05 1.3894049203809568e-05 2.3985465221727954e-05 3.3458449092465795e-05 0.00028172299406359233 0.00040056109827889085 0.0004621959325200118 0.0034637215474633033 0.009153352872912168 0.012804683731760212 0.010674833967812809 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.008855607599180339 0.015278371659152929 0.0023879884417463693 0.000935102721182137 0.0034651523786104546 0.005282714096158778 0.0012741283037854573 0.006268847413803995 0.004013755514818873 0.010363477891131097 0.007109323912817858 -0.0002600093944564617 -0.00034633976332068713 0.0 0.0 0.0 0.0 0.0 0.0
+5.3159155089314414e-05 1.2012438398909825e-05 3.274115317951001e-05 1.689353008824745e-05 -4.232110203859359e-07 7.656160171407207e-06 5.754938338062256e-05 0.0002687156220968384 9.768544317740195e-05 1.0129360414562531e-05 2.293272526112481e-05 1.8260890221186993e-05 1.7288534885724222e-06 1.8856324749638164e-06 0.0 0.0 0.0 0.0 0.0
+0.00033854327480422193 0.0032940012133255356 0.008932075729876752 0.012661798131960687 0.013919159699477152 0.08208818801401566 0.15644219608737447 0.4341503084393359 0.557964838826116 0.21173866865770563 0.03478604116524652 -0.007972916161324952 -0.0007132020154210059 0.0004420932693293155 0.0003860741867263207 0.0 0.0 0.0 0.0
+0.00018639122271027446 0.0034767153815636618 0.018231622622978436 0.0042988446466234575 0.05300133577632108 0.7250033316881788 0.4935009383219143 0.42606157195551264 1.5039118559972142 1.2303936880370434 0.22548363428435172 -0.30298333788301807 -0.13073656424500055 -0.010403074319359695 0.005059645288601829 0.0015181780063355998 0.0 0.0 0.0
+0.0005854549915295117 0.008715343377777027 0.015086646356161998 0.00551449475752511 0.16771008674532747 0.9262472909137329 0.16287454314633398 0.10026082406557575 0.8138824998965698 0.5734788502649438 0.5447903546528722 0.10618098607332271 -0.9103493511412221 -0.576660316279193 -0.03402023384277208 -0.0031478924709104684 0.0 0.0 0.0
+-0.00011090752856021898 0.0043805342288213535 0.0016927729283659975 0.012954305343473369 0.4368797400066579 0.7614408833500013 0.10944765261392181 0.0412051224385603 0.15805652281823218 0.2585893276179897 0.05949650977132904 0.30356693149425945 0.19493804255113664 0.5122323478146567 0.35279041052227494 -0.007488446744915854 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 5.196157077391817e-07 6.938124100654148e-08 7.597564197383545e-08 1.863740632660483e-07 4.437837629589167e-07 5.453941063185757e-07 1.5602917821833568e-06 3.404289212094662e-07 9.967447994956849e-07 5.8845599651090215e-06 1.5052240335012455e-05 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.005346497059990333 -0.002010348201210142 -0.0010943235863089423 0.0001661513182702165 0.00012025969610516196 2.4949866002221845e-05 6.627236360802077e-06 3.003757825105864e-06 3.997348910159012e-05 0.000427961841918743 0.0007451357800599296 0.0011219432594133996 0.009685550613014016 0.0 0.0 0.0 0.0 0.0 0.0
+-2.698386662730078e-05 2.754973422173369e-07 0.0001320534807487939 0.004372176148977807 0.0015642026255259442 0.004248632573013906 0.00040885420395593786 0.001088966135412402 0.004766232525411325 0.002205157762668968 0.007459939889093756 0.005587608653898612 -0.004720162133268877 -0.006593174803103767 0.0 0.0 0.0 0.0 0.0
+4.542360076931743e-06 1.2329326209575631e-05 1.7173803033436737e-05 5.0102066463061734e-05 1.09067765324765e-05 4.930240898900306e-05 0.00028721835291257015 0.0007503332386451459 0.00015238128535605624 -6.391111549761724e-05 -4.677072820313549e-06 -5.572104125200205e-06 2.707663268609677e-05 5.022846595129856e-05 0.00012970076559689836 0.0 0.0 0.0 0.0
+0.0005943091729493132 0.0009672866630600067 0.009096489467732383 0.01304005169719466 0.01734885651081947 0.1652472792755658 0.44880261475702005 0.6329167707872334 0.5244926606398645 0.26270394857828266 -0.10487776273847933 -0.05845317239353218 0.0032135333670803676 0.003777214235598332 0.0007681993725802362 0.00011959437074006901 0.0 0.0 0.0
+0.0005854549915295117 0.008715343377777027 0.015086646356161998 0.00551449475752511 0.16771008674532747 0.9262472909137329 0.16287454314633398 0.10026082406557575 0.8138824998965698 0.5734788502649438 0.5447903546528722 0.10618098607332271 -0.9103493511412221 -0.576660316279193 -0.03402023384277208 -0.0031478924709104684 0.0 0.0 0.0
+-0.011510449622067839 -0.0007610595804959427 1.0172131902385016e-05 0.00448473230635448 0.2199020425072735 0.07466247888373397 0.21339517449435372 0.013531084627798973 0.04473458040783101 0.2378714243611067 0.10478901497777 0.37112701147924365 0.27817806337533985 -0.2385293501359752 -0.33578482057600856 -0.004009555340792723 0.0 0.0 0.0
+4.1934884152796484e-05 0.0035902058234846823 0.0005824017329515572 0.01441762880686721 0.6293209695433368 0.09572421881003676 0.08586090788885907 0.031601835145355577 0.06121858995657405 0.07315778164682188 -0.0013308261229131521 0.012882908917780034 0.13273058077548822 0.2035107778991338 0.1961396131043037 -0.04170252478795313 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 3.196534127251379e-06 -1.1872677866681306e-06 -2.5678892066098854e-08 5.139117071436217e-09 1.1142431390092631e-06 2.0605776537608227e-06 5.297265009242829e-06 6.713907186544732e-06 2.7028644452395994e-06 1.149242737988068e-06 2.2449682976927855e-06 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.013431011504370738 0.007883426617122005 0.007935899204760883 0.005880150773602205 0.0009832099103910489 0.005414528729313218 0.0015950126575825377 0.00024127039666882992 -1.5674461809944553e-05 1.0711548076574028e-05 7.830483572860064e-06 0.00011012649333888752 0.0005497452692208139 0.0 0.0 0.0 0.0 0.0
+-7.940872165606751e-05 2.276762148612182e-06 0.00010635762128769112 2.049233578255131e-05 0.00042145671490654473 0.01249692872936893 0.0020370917425772224 0.0017316344055948985 0.0006574004028558345 0.0012283310563930355 0.0014629455315045585 -4.492954039177435e-06 0.00029547002108771967 0.002683727758662211 0.004096710661285439 0.0 0.0 0.0 0.0
+0.0001650881267658455 6.548635015912796e-05 3.282822556024051e-05 7.205996516588195e-05 0.00030710620226134084 0.0007279352507540159 0.000909918529220897 0.0017257383928761386 0.00047543969768972346 -0.00030462814537952123 2.7321207199326783e-05 0.00010130812246147248 5.4817489360932934e-05 0.00012640103175376577 0.00018995742264027741 0.00023929212126678798 0.0 0.0 0.0
+6.489304732024981e-05 0.0011748014116558024 0.014162497237899634 0.024470256017495366 0.046651858009005745 0.46415330969555396 0.6536160620847673 0.38295932353650225 0.39318393240673155 0.2877811024545165 0.03747496490739291 0.2654106688863148 0.074634565935104 0.012330295300167044 -0.0004996457463809098 0.00023089219471653216 0.0 0.0 0.0
+-0.00011090752856021898 0.0043805342288213535 0.0016927729283659975 0.012954305343473369 0.4368797400066579 0.7614408833500013 0.10944765261392181 0.0412051224385603 0.15805652281823218 0.2585893276179897 0.05949650977132904 0.30356693149425945 0.19493804255113664 0.5122323478146567 0.35279041052227494 -0.007488446744915854 0.0 0.0 0.0
+4.1934884152796484e-05 0.0035902058234846823 0.0005824017329515572 0.01441762880686721 0.6293209695433368 0.09572421881003676 0.08586090788885907 0.031601835145355577 0.06121858995657405 0.07315778164682188 -0.0013308261229131521 0.012882908917780034 0.13273058077548822 0.2035107778991338 0.1961396131043037 -0.04170252478795313 0.0 0.0 0.0
+-0.0008465335016788498 -9.208409992139663e-05 -1.1210629044433908e-05 0.004300763141885697 0.1800503541691201 0.016819879476467067 0.060737350690215776 0.042330087306607714 0.0007006095666399378 0.02251352126872946 0.011185142175963118 -0.02716513333374441 0.00952941322650142 0.02210443704516212 0.07011148238258197 0.057840211345517194 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+#
\ No newline at end of file
diff --git a/lammps_plugin/unittest/README.rst b/lammps_plugin/unittest/README.rst
new file mode 100644
index 00000000..89aa0c46
--- /dev/null
+++ b/lammps_plugin/unittest/README.rst
@@ -0,0 +1,20 @@
+This document describes how to add ML-UF3 unittest to LAMMPS. Refer to `this link <https://docs.lammps.org/Developer_unittest.html>`_ for more information on unittest in LAMMPS.
+
+To run the ML-UF3 unittest, copy the contents of :code:`unittest` directory (i.e :code:`A_A.uf3_pot`, :code:`A_A_A.uf3_pot` and :code:`manybody-pair-uf3.yaml`) to the unittest directory of LAMMPS and re-build LAMMPS-
+
+.. code:: bash
+
+    cp unittest/* LAMMPS_BASE_DIR/unittest/force-styles/tests/
+    cd LAMMPS_BASE_DIR/build
+    cmake ../cmake/ -D PKG_ML-UF3=yes -D ENABLE_TESTING=on
+    cmake --build .
+
+After the build is sucessful go to the LAMMPS unittest directory for pair style and run-
+
+.. code:: bash
+
+    cd LAMMPS_BASE_DIR/unittest/force-styles/tests/
+    LAMMPS_BASE_DIR/build/test_pair_style manybody-pair-uf3.yaml
+
+
+The current implementation will only pass 2 out of 8 tests, remaining 6 tests will be skipped.
diff --git a/lammps_plugin/unittest/manybody-pair-uf3.yaml b/lammps_plugin/unittest/manybody-pair-uf3.yaml
new file mode 100644
index 00000000..259f66e6
--- /dev/null
+++ b/lammps_plugin/unittest/manybody-pair-uf3.yaml
@@ -0,0 +1,159 @@
+---
+lammps_version: 24 Mar 2022
+tags: generated
+date_generated: Thu Jun  8 23:16:30 2023
+epsilon: 1e-9
+skip_tests: ""
+prerequisites: ! |
+  pair uf3
+pre_commands: ! |
+  variable newton_pair delete
+  variable newton_pair index on
+  variable newton_bond delete
+post_commands: ! ""
+input_file: in.manybody
+pair_style: uf3 3 8
+pair_coeff: ! |
+  * * A_A.uf3_pot
+  * * * A_A_A.uf3_pot
+extract: ! ""
+natoms: 64
+init_vdwl: -76.14388662099438
+init_coul: 0
+init_stress: ! |2-
+   3.1223073343802071e+02  3.1503555484293474e+02  3.2087032195384182e+02 -5.2677023646012433e+00  4.1046361968856566e+01 -2.2705704820012654e-01
+init_forces: ! |2
+    1 -1.0963106297354930e+00  1.9921565797217811e+00  2.0176595423650685e+00
+    2 -2.5744974244934786e+00  8.6065313692841872e-01 -1.3343920771683084e+00
+    3 -7.5762202587571881e-01 -5.1086473746213934e-01  1.7774798100697495e+00
+    4 -1.2651178900120015e+00  2.5481168050091734e+00  1.0332353551246649e+00
+    5 -3.3301075059618213e-02 -8.6936885426915711e-01 -8.0361144939346540e-01
+    6  5.9576544655966956e-02 -5.7569733665007693e-02 -2.5260577270195245e-01
+    7 -7.8223973805485159e-01 -1.5872724248886485e+00 -4.0690678808175756e-01
+    8  1.2837377243355602e-01  6.8616887866365453e-02  3.9279992110159728e-02
+    9  1.3013607224018784e+00 -2.2234802020121042e-01 -2.5820065882172409e+00
+   10 -3.5389541898719123e-01  1.0418734269769769e+00 -6.7534263859128518e-01
+   11  1.5048713773196754e+00 -1.1607474819622305e+00 -4.1343086960946002e-01
+   12 -3.8382036374205457e+00 -1.8953194768009614e+00 -1.5975045274049304e+00
+   13 -5.1243974635480005e-01  2.8005048190722506e+00 -8.5188627697957486e-01
+   14 -4.3260615266171099e-01  1.7104546503175775e+00  4.8896255947620082e-01
+   15 -2.1729171109783061e+00  1.9610148525080886e+00 -1.7550881080125238e+00
+   16  1.3063694982664014e+00  4.0801741093678640e-01  3.3816303090167859e+00
+   17  1.0103076935514768e+00  1.3062727533124225e+00 -3.0591856171743643e+00
+   18  6.6364797168941592e-01  1.3269858177825717e+00  1.0859880651465830e+00
+   19 -4.3237772669504843e-01  1.5179601478654310e+00 -8.9399081937433889e-01
+   20 -4.1109776137187977e+00  4.8942123015814726e-01  3.0324705210921541e-02
+   21  1.0443613152888414e+00 -2.7611218991110009e-01 -3.2334333893456133e+00
+   22 -3.8840635954358733e+00  1.8543888965122868e+00 -2.7044170995178298e+00
+   23  6.4824470795616951e-01  1.1930964693491897e+00  2.3472683895454201e+00
+   24  5.6518885203578750e-01  8.9024666222906623e-01  9.9558410495963645e-01
+   25 -2.3884920507811258e-01 -1.2236748552119361e-01  9.6166740424011798e-04
+   26 -9.4060753459907698e-01 -2.9688412181022056e-01  6.7007032584453752e-01
+   27  5.0390685311588923e-01 -5.8477364114704944e-01  1.6320689076693460e+00
+   28 -8.1564347181637331e-01  1.9951358496458818e-01  3.2879220126777038e+00
+   29 -1.6141037582891811e+00  3.7621933923780954e-01  5.6865013193370151e-01
+   30 -1.2051389982059610e+00  6.5471000907223187e-02  4.6548238063076404e-02
+   31 -8.3799333021355227e-01 -4.3387478264068147e-01 -3.4488785440063413e-01
+   32 -1.8712481645378094e+00  3.1431436055998407e-02  1.1920833582466677e+00
+   33  2.1132311819974756e+00 -7.8762333052322075e-01  2.9319369786200289e+00
+   34 -7.4212669325944880e-01  1.5089695276247311e-01  6.4837671979385259e-01
+   35  4.7876606586762549e-01 -5.2894027298845681e-01  6.5312176822087242e-01
+   36  9.4335406181387005e-01 -7.7270400013223828e-01  4.8506498341757304e-01
+   37 -6.9776523536821422e-01  1.5814045923629079e+00 -1.7145687025150753e+00
+   38  8.3390581678419395e-01 -2.3460369438656256e-01 -8.1845978792256724e-01
+   39  7.9082785819764490e-01 -1.4014821253885934e+00 -3.6171284136791626e+00
+   40  1.5528260681499937e+00 -1.9423308463414859e+00  2.7454733155675826e-01
+   41  1.6230449781222470e-01  4.0254440068393893e-01 -1.0728426614941826e+00
+   42 -2.5760797297848943e+00 -2.3727612877133377e+00 -7.3694736943877159e-01
+   43 -1.4658113294445060e-01 -3.7059450064886161e+00 -2.2611910919567893e-01
+   44  2.6125365540590240e+00 -2.0393457867642488e+00 -1.7353429519549574e-01
+   45  9.1508218547652620e-01  3.2090468646350390e-01 -1.3730810235354041e+00
+   46 -1.2861926535120600e+00 -1.8325412123528377e+00  2.7409156132103112e+00
+   47  5.2214882788544981e-01  9.0702150750152088e-02  1.3758849361839385e+00
+   48  8.7060691872545093e-01  1.0333012026994193e+00 -8.8450736609033931e-01
+   49  1.6921455867723978e-01 -3.2865843167979367e+00 -2.3941507623279072e+00
+   50  1.3029435213640246e+00  1.1566980491369294e+00 -7.8373321422495534e-02
+   51  1.7782167191801962e-01  3.5869618077998595e+00  2.1417753790319543e+00
+   52  3.4939214706481048e+00 -2.4127970289820255e+00  4.1443505260596725e+00
+   53  4.5955988135622799e-01  1.4913218496577223e+00  1.3076728090591363e+00
+   54  7.0400726037068106e-01 -9.4898269328552198e-01 -7.0340401843204670e-01
+   55  1.1319661424097816e+00 -2.3442414589969114e+00 -7.7166661047173946e-01
+   56  7.7292681497946214e-01  5.4202239892193216e-01 -1.0429033367200278e+00
+   57  1.1664627895682855e+00  1.9743121270468009e-01 -3.8302192241786348e-01
+   58  1.3342985230821185e+00 -2.6808564460978351e-01 -9.3875599645237040e-01
+   59  2.1547446695381884e+00  2.5922697594917221e-01  5.5883175680123842e-01
+   60 -1.1202730918333845e+00 -4.0836109063032069e+00 -3.1463158656990915e+00
+   61  7.5258775499505959e-01  1.4742795022217277e+00 -2.2056849646259416e-01
+   62 -2.1194607090574338e+00  9.6304617778841872e-01 -5.8648933450219842e-01
+   63  2.2116847362243819e+00 -6.9485816680348522e-01 -1.2888780585377166e+00
+   64  2.0946943533672595e+00  1.7817828615230797e+00  5.2222100516662051e+00
+run_vdwl: -76.1335425447406
+run_coul: 0
+run_stress: ! |2-
+   3.1227357278065733e+02  3.1510436259931976e+02  3.2097655273455166e+02 -5.2932374901106582e+00  4.0956468170617640e+01  1.5018802669860862e-01
+run_forces: ! |2
+    1 -1.1136068191144672e+00  1.9868928531451695e+00  2.0284715885889444e+00
+    2 -2.5718148688418596e+00  8.3532761170683545e-01 -1.3404690568364115e+00
+    3 -7.6676194475345183e-01 -5.0716795928532021e-01  1.7997537336742999e+00
+    4 -1.2720923579033645e+00  2.5419855255447907e+00  1.0258862987566395e+00
+    5 -4.3524258688024453e-02 -8.8004954612209696e-01 -7.9382095260610019e-01
+    6  8.9678254474529928e-02 -5.0742940694548030e-02 -2.4348008365156723e-01
+    7 -7.7364435925734953e-01 -1.5558296160053275e+00 -4.2051790182009818e-01
+    8  1.1143851887986839e-01  7.2497057799814191e-02  3.6624183488660304e-02
+    9  1.2839141188983776e+00 -2.4381166390795905e-01 -2.5926851536420066e+00
+   10 -3.6606370803571597e-01  1.0472513106007628e+00 -6.5579353351390024e-01
+   11  1.4796344435746969e+00 -1.1469564140842350e+00 -4.2766370678100751e-01
+   12 -3.8341194520315707e+00 -1.8877013457658414e+00 -1.5901839974079914e+00
+   13 -4.6650968579193675e-01  2.7962961910932282e+00 -8.2237329615475063e-01
+   14 -4.5362402854231521e-01  1.7048830533652388e+00  4.9165017205992945e-01
+   15 -2.1790730324228798e+00  1.9552241266001966e+00 -1.7453842965710837e+00
+   16  1.3028733875052598e+00  4.1039519382071737e-01  3.3783851223727419e+00
+   17  1.0049330607671125e+00  1.3021263787131616e+00 -3.0544446010601831e+00
+   18  6.3052176771903379e-01  1.3324251901771336e+00  1.0977786571785106e+00
+   19 -4.5964705284123741e-01  1.4850585747490403e+00 -8.9955925151919691e-01
+   20 -4.0843338981113764e+00  4.6843741575982639e-01  4.6947857407657889e-03
+   21  1.0581909871649149e+00 -2.9084634536392007e-01 -3.2540160500533157e+00
+   22 -3.8927799316015426e+00  1.8587085622785937e+00 -2.7125312751977178e+00
+   23  6.6350903350923252e-01  1.1973095378961089e+00  2.3524290511459478e+00
+   24  5.7143336751982488e-01  8.7715571174154605e-01  9.9432218349388102e-01
+   25 -2.4703347062300032e-01 -1.2710986429222224e-01 -1.5749033311971992e-02
+   26 -9.3451713960198701e-01 -2.7778133649821901e-01  6.6713867661147919e-01
+   27  4.9808843860149610e-01 -5.8536464204160887e-01  1.6301914265685915e+00
+   28 -8.2771277475903027e-01  2.1570720428840873e-01  3.2880151472957344e+00
+   29 -1.6048620343875755e+00  3.8486480806372897e-01  5.6967398899129529e-01
+   30 -1.1973961350199096e+00  5.1593120685595580e-02  4.1293831832017062e-02
+   31 -8.4115448475027121e-01 -4.4444260925279200e-01 -3.3858426829044269e-01
+   32 -1.8597577591090164e+00  1.2810085646854485e-02  1.1797889462030640e+00
+   33  2.1159519472471811e+00 -7.8729199670032701e-01  2.9290939088097181e+00
+   34 -7.4188497270023746e-01  1.5645497560825464e-01  6.6033973957472436e-01
+   35  4.8203360935099837e-01 -5.2533146218590032e-01  6.5589135580856639e-01
+   36  9.4037223416537397e-01 -7.6094150822319195e-01  4.5734538655438040e-01
+   37 -6.9906373360474205e-01  1.5746536313295925e+00 -1.7263112595330459e+00
+   38  8.5615988284237132e-01 -2.2431722964535125e-01 -8.0332888328255958e-01
+   39  7.8065137836247200e-01 -1.3983715176027893e+00 -3.6344006191116991e+00
+   40  1.5756754145077525e+00 -1.9548680562354619e+00  2.4171639489710248e-01
+   41  1.7946227297557377e-01  4.1626928569421628e-01 -1.0542445180049111e+00
+   42 -2.5855819286639807e+00 -2.3809591696792194e+00 -7.5444018878363772e-01
+   43 -1.5651979567151933e-01 -3.7133196766462597e+00 -2.3417783002479109e-01
+   44  2.6088788324017109e+00 -2.0370826629544276e+00 -1.7748087222007777e-01
+   45  9.0750981150111454e-01  3.1320571808181447e-01 -1.3570279945304848e+00
+   46 -1.2844641842483453e+00 -1.8271514736854049e+00  2.7652384797305016e+00
+   47  5.2355177969578193e-01  9.5756521123463834e-02  1.3725765177085687e+00
+   48  8.7731662768264451e-01  1.0414860079651591e+00 -9.0032888888892382e-01
+   49  1.4644557239036082e-01 -3.3233388609839873e+00 -2.4250592659007468e+00
+   50  1.3159886335597091e+00  1.1534831728413786e+00 -7.9023436269654135e-02
+   51  1.5810961991408728e-01  3.6227053406802825e+00  2.1622129324293375e+00
+   52  3.5238822669377128e+00 -2.4112486219526210e+00  4.1691651240037055e+00
+   53  4.6279178075715538e-01  1.5349821135997805e+00  1.3624008261786760e+00
+   54  7.0473247201702627e-01 -9.3593778743269240e-01 -7.0089892988315006e-01
+   55  1.1161312397801737e+00 -2.3385526086116117e+00 -7.5151515336312336e-01
+   56  7.6996392118239054e-01  5.3225925152027631e-01 -1.0367667800909899e+00
+   57  1.1705705890031866e+00  1.9234764366836088e-01 -3.8394944778612211e-01
+   58  1.3292217712869172e+00 -2.6420177775466325e-01 -9.5699172084102668e-01
+   59  2.1842122534627268e+00  2.6220494593377625e-01  5.4393643560229021e-01
+   60 -1.1457317620961742e+00 -4.1194027577964150e+00 -3.1920012226137322e+00
+   61  7.5501968427331045e-01  1.4824037577796831e+00 -2.0845605296052469e-01
+   62 -2.1181153169676303e+00  9.6226604332888710e-01 -5.7042120154066434e-01
+   63  2.2161829432347631e+00 -7.0548049758710407e-01 -1.3062673143062038e+00
+   64  2.1263589729936818e+00  1.8281740261598369e+00  5.2543331430537403e+00
+...