Repository: uf3/uf3 (language: python)

Issue:
Keys can only be integers in the current uf3 demo
Currently the demo Jupyter notebook can only process keys of integer type, which might need some modifications so that it can be used with json files with predefined data split (like mlearn data)


Discussion / hints:

Fix patch:
diff --git a/lammps_plugin/ML-UF3/uf3_pair_bspline.h b/lammps_plugin/ML-UF3/uf3_pair_bspline.h
index 046c07d7..aa3f1e8c 100644
--- a/lammps_plugin/ML-UF3/uf3_pair_bspline.h
+++ b/lammps_plugin/ML-UF3/uf3_pair_bspline.h
@@ -33,10 +33,9 @@ class uf3_pair_bspline {
   std::vector<double> coeff_vect, dncoeff_vect;
   std::vector<uf3_bspline_basis3> bspline_bases;
   std::vector<uf3_bspline_basis2> dnbspline_bases;
-  int knot_spacing_type;
   int get_starting_index_uniform(double), get_starting_index_nonuniform(double);
   int (uf3_pair_bspline::*get_starting_index)(double);
-  double knot_spacing=0; 
+  //double knot_spacing=0; 
   LAMMPS *lmp;
 
  public:
@@ -46,6 +45,8 @@ class uf3_pair_bspline {
                    const std::vector<double> &ucoeff_vect,
                    const int &uknot_spacing_type);
   ~uf3_pair_bspline();
+  int knot_spacing_type;
+  double knot_spacing=0;
   double ret_val[2];
   double *eval(double value_rij);
   double memory_usage();
diff --git a/lammps_plugin/ML-UF3/uf3_triplet_bspline.cpp b/lammps_plugin/ML-UF3/uf3_triplet_bspline.cpp
index 642bcce3..6c5a5a19 100644
--- a/lammps_plugin/ML-UF3/uf3_triplet_bspline.cpp
+++ b/lammps_plugin/ML-UF3/uf3_triplet_bspline.cpp
@@ -181,11 +181,17 @@ double *uf3_triplet_bspline::eval(double value_rij, double value_rik, double val
   ret_val[3] = 0;
 
   for (int i = 0; i < 4; i++) {
+    const double basis_iji = basis_ij[i]; // prevent repeated access of same memory location
     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];
-      }
+      const double factor = basis_iji * basis_ik[j]; // prevent repeated access of same memory location
+      const double* slice = &coeff_matrix[i + iknot_ij][j + iknot_ik][iknot_jk]; // declare a contigues 1D slice of memory
+      double tmp[4]; // declare tmp array that holds the 4 tmp values so the can be computed simultaniously in 4 separate registeres.
+      tmp[0] = slice[0] * basis_jk[0];
+      tmp[1] = slice[1] * basis_jk[1];
+      tmp[2] = slice[2] * basis_jk[2];
+      tmp[3] = slice[3] * basis_jk[3];
+      double sum = tmp[0] + tmp[1] + tmp[2] + tmp[3];
+      ret_val[0] += factor * sum; // use 1 fused multiply-add (FMA)
     }
   }
 
@@ -210,29 +216,46 @@ double *uf3_triplet_bspline::eval(double value_rij, double value_rik, double val
   dnbasis_jk[3] = 0;
 
   for (int i = 0; i < 3; i++) {
+    const double dnbasis_iji = dnbasis_ij[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];
-      }
+      const double factor = dnbasis_iji * basis_ik[j];
+      const double* slice = &dncoeff_matrix_ij[iknot_ij + i][iknot_ik + j][iknot_jk];
+      double tmp[4];
+      tmp[0] = slice[0] * basis_jk[0];
+      tmp[1] = slice[1] * basis_jk[1];
+      tmp[2] = slice[2] * basis_jk[2];
+      tmp[3] = slice[3] * basis_jk[3];
+      double sum = tmp[0] + tmp[1] + tmp[2] + tmp[3];
+      ret_val[1] += factor * sum;
     }
   }
 
   for (int i = 0; i < 4; i++) {
+    const double basis_iji = basis_ij[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];
-      }
+      const double factor = basis_iji * dnbasis_ik[j];
+      const double* slice = &dncoeff_matrix_ik[iknot_ij + i][iknot_ik + j][iknot_jk];
+      double tmp[4];
+      tmp[0] = slice[0] * basis_jk[0];
+      tmp[1] = slice[1] * basis_jk[1];
+      tmp[2] = slice[2] * basis_jk[2];
+      tmp[3] = slice[3] * basis_jk[3];
+      double sum = tmp[0] + tmp[1] + tmp[2] + tmp[3];
+      ret_val[2] += factor * sum;
     }
   }
 
   for (int i = 0; i < 4; i++) {
+    const double basis_iji = basis_ij[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];
-      }
+      const double factor = basis_iji * basis_ik[j];
+      const double* slice = &dncoeff_matrix_jk[iknot_ij + i][iknot_ik + j][iknot_jk];
+      double tmp[3];
+      tmp[0] = slice[0] * dnbasis_jk[0];
+      tmp[1] = slice[1] * dnbasis_jk[1];
+      tmp[2] = slice[2] * dnbasis_jk[2];
+      double sum = tmp[0] + tmp[1] + tmp[2];
+      ret_val[3] += factor * sum;
     }
   }
 
diff --git a/lammps_plugin/ML-UF3/uf3_triplet_bspline.h b/lammps_plugin/ML-UF3/uf3_triplet_bspline.h
index bdf30f82..f5cf6a1b 100644
--- a/lammps_plugin/ML-UF3/uf3_triplet_bspline.h
+++ b/lammps_plugin/ML-UF3/uf3_triplet_bspline.h
@@ -30,11 +30,11 @@ class uf3_triplet_bspline {
   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;
-  int knot_spacing_type;
   int get_starting_index_uniform(const std::vector<double>, int, double, double);
   int get_starting_index_nonuniform(const std::vector<double>, int, double, double);
   int (uf3_triplet_bspline::*get_starting_index)(const std::vector<double>, int, double, double);
-  double knot_spacing_ij=0,knot_spacing_ik=0,knot_spacing_jk=0; 
+  //double knot_spacing_ij=0,knot_spacing_ik=0,knot_spacing_jk=0; 
+  //double _alignvar(, 8) ret_val[4]; // Force memory alignment on 8 byte boundaries
   double ret_val[4];
 
   int starting_knot(const std::vector<double>, int, double);
@@ -46,6 +46,8 @@ class uf3_triplet_bspline {
                       const std::vector<std::vector<std::vector<double>>> &ucoeff_matrix,
                       const int &uknot_spacing_type);
   ~uf3_triplet_bspline();
+  int knot_spacing_type;
+  double knot_spacing_ij=0,knot_spacing_ik=0,knot_spacing_jk=0;
   double *eval(double value_rij, double value_rik, double value_rjk);
 
   double memory_usage();
diff --git a/lammps_plugin/ML-UF3_KOKKOS/pair_uf3_kokkos.cpp b/lammps_plugin/ML-UF3_KOKKOS/pair_uf3_kokkos.cpp
new file mode 100644
index 00000000..dbcdd22e
--- /dev/null
+++ b/lammps_plugin/ML-UF3_KOKKOS/pair_uf3_kokkos.cpp
@@ -0,0 +1,1722 @@
+/* ----------------------------------------------------------------------
+   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_kokkos.h"
+
+#include "atom_kokkos.h"
+#include "atom_masks.h"
+#include "comm.h"
+#include "error.h"
+#include "force.h"
+#include "kokkos.h"
+#include "kokkos_type.h"
+#include "math_const.h"
+#include "memory.h"
+#include "memory_kokkos.h"
+#include "neigh_list_kokkos.h"
+#include "neigh_request.h"
+#include "neighbor.h"
+#include "pair_kokkos.h"
+#include "text_file_reader.h"
+#include <algorithm>
+#include <cmath>
+#include <utility>
+
+#include <cmath>
+
+using namespace LAMMPS_NS;
+using namespace MathConst;
+
+template <class DeviceType> PairUF3Kokkos<DeviceType>::PairUF3Kokkos(LAMMPS *lmp) : PairUF3(lmp)
+{
+  respa_enable = 0;
+
+  //kokkosable = 1;
+  atomKK = (AtomKokkos *) atom;
+  execution_space = ExecutionSpaceFromDevice<DeviceType>::space;
+  datamask_read = X_MASK | F_MASK | TAG_MASK | TYPE_MASK | ENERGY_MASK | VIRIAL_MASK;
+  datamask_modify = F_MASK | ENERGY_MASK | VIRIAL_MASK;
+}
+
+template <class DeviceType> PairUF3Kokkos<DeviceType>::~PairUF3Kokkos()
+{
+  if (!copymode) {
+    memoryKK->destroy_kokkos(k_eatom, eatom); //destory eatom from host, set it to nullptr
+                                              //Also set k_eatom to empty View
+    memoryKK->destroy_kokkos(k_vatom, vatom);
+    memoryKK->destroy_kokkos(k_cutsq,cutsq);
+    destroy_3d(k_cut_3b,cut_3b);
+    destroy_4d(k_min_cut_3b,min_cut_3b);
+    eatom = NULL;
+    vatom = NULL;
+    cvatom = NULL;
+  }
+}
+
+template <class DeviceType>
+template <typename TYPE> 
+void PairUF3Kokkos<DeviceType>::destroy_3d(TYPE data, typename TYPE::value_type*** &array)
+{
+  if (array == nullptr) return;
+  data = TYPE();
+  memory->sfree(array);
+  array = nullptr;
+}
+
+template <class DeviceType>
+template <typename TYPE> 
+void PairUF3Kokkos<DeviceType>::destroy_4d(TYPE data, typename TYPE::value_type**** &array)
+{
+  if (array == nullptr) return;
+  data = TYPE();
+  memory->sfree(array);
+  array = nullptr;
+}
+
+/* ----------------------------------------------------------------------
+ *     global settings
+ * ---------------------------------------------------------------------- */
+
+template <class DeviceType> void PairUF3Kokkos<DeviceType>::settings(int narg, char **arg)
+{
+  PairUF3::settings(narg, arg);
+  //1. Determines whether the simulation is 2-body or 2 and 3-body
+  //2. Set nbody_flag, num_of_elements, pot_3b
+}
+
+/* ----------------------------------------------------------------------
+ *    set coeffs for one or more type pairs
+ * ---------------------------------------------------------------------- */
+template <class DeviceType> void PairUF3Kokkos<DeviceType>::coeff(int narg, char **arg)
+{
+  PairUF3::coeff(narg,arg);
+  //Also calls allocate internally
+  //Grows arrays to the right dimensions --> setflag, cutsq, cut, knot_spacing_type_2b,
+  //n2b_knot, n2b_coeff, UFBS2b, n2b_knot[i], n2b_coeff[i], UFBS2b[i], setflag_3b,
+  //cut_3b, cut_3b_list, min_cut_3b, knot_spacing_type_3b, cut_3b_list, n3b_knot_matrix,
+  //UFBS3b, neighshort
+  //
+  //Also reads the pot_files_internally
+
+
+/*  copy_2d(d_cutsq, cutsq, num_of_elements + 1, num_of_elements + 1); //copy cutsq from
+  //host to device memory
+  if (pot_3b) {
+    copy_3d(d_cut_3b, cut_3b, num_of_elements + 1, num_of_elements + 1, num_of_elements + 1);
+    copy_2d(d_cut_3b_list, cut_3b_list, num_of_elements + 1, num_of_elements + 1);
+  } else {
+    //why are we allocating space to d_cut_3b, d_cut_3b_list if the pot is 2-body only?
+    Kokkos::realloc(d_cut_3b, num_of_elements + 1, num_of_elements + 1, num_of_elements + 1);
+    Kokkos::realloc(d_cut_3b_list, num_of_elements + 1, num_of_elements + 1);
+  }
+  //No allocation for device equivalent of --> setflag, cut, knot_spacing_type_2b,
+  //n2b_knot, n2b_coeff, n2b_knot[i], n2b_coeff[i], setflag_3b, cut_3b, 
+  //cut_3b_list, min_cut_3b, knot_spacing_type_3b, cut_3b_list, n3b_knot_matrix,
+  //neighshort
+  
+  //UFBS2b and UFBS3b are array of objects. Bad idea to use kokkos view(array)
+  //for it
+  create_2b_coefficients();
+  if (pot_3b) create_3b_coefficients();*/
+}
+
+template<class DeviceType>
+void PairUF3Kokkos<DeviceType>::allocate()
+{
+  if (!allocated) PairUF3::allocate();
+  int n = atom->ntypes;
+  memory->destroy(cutsq); //Why are we destroying cutsq? cutsq is allocated when
+  //PairUF3::coeff or PairUF3::allocate is called; in the next step when k_cutsq
+  //is created cut_3b is set to point to the host array of k_cutsq
+  //memory->destroy(cut_3b);
+
+  memoryKK->create_kokkos(k_cutsq,cutsq,n+1,n+1,"pair:cutsq"); 
+  d_cutsq = k_cutsq.template view<DeviceType>(); //assignment; get the device
+  //view of k_cutsq and assign it to d_cutsq; in the header file we just
+  //decleared d_cutsq's type
+  memoryKK->create_kokkos(k_cut_3b,n+1,n+1,n+1,"threebody:cut");
+  memoryKK->create_kokkos(k_min_cut_3b,n+1,n+1,n+1,3,"threebody:cut");
+  d_cut_3b = k_cut_3b.template view<DeviceType>();
+  d_min_cut_3b = k_min_cut_3b.template view<DeviceType>();
+}
+
+
+/* ----------------------------------------------------------------------
+   init specific to this pair style
+------------------------------------------------------------------------- */
+
+template <class DeviceType> void PairUF3Kokkos<DeviceType>::init_style()
+{
+
+  PairUF3::init_style();
+
+  neighflag = lmp->kokkos->neighflag;
+
+  auto request = neighbor->find_request(this);
+  request->set_kokkos_host(std::is_same<DeviceType, LMPHostType>::value &&
+                           !std::is_same<DeviceType, LMPDeviceType>::value);
+  request->set_kokkos_device(std::is_same<DeviceType, LMPDeviceType>::value);
+
+  request->enable_full();
+  request->enable_ghost();
+}
+
+/* ----------------------------------------------------------------------
+   init list sets the pointer to full neighbour list requested in previous function
+------------------------------------------------------------------------- */
+
+template <class DeviceType>
+void PairUF3Kokkos<DeviceType>::init_list(int /*id*/, class NeighList *ptr)
+{
+  list = ptr;
+}
+
+/* ----------------------------------------------------------------------
+   init for one type pair i,j and corresponding j,i
+------------------------------------------------------------------------- */
+template <class DeviceType> double PairUF3Kokkos<DeviceType>::init_one(int i, int j)
+{
+  double cutone = PairUF3::init_one(i, j);
+
+  if (!coefficients_created) create_coefficients();
+
+  k_cutsq.h_view(i,j) = k_cutsq.h_view(j,i) = cutone*cutone; //Update the k_cutsq's
+  //host memory
+  k_cutsq.template modify<LMPHostType>(); //Record that k_cutsq's host memory has
+  //been modified
+
+  return cutone;
+}
+
+template <class DeviceType> void PairUF3Kokkos<DeviceType>::create_coefficients()
+{
+  coefficients_created = 1;
+  
+  /*for (int i = 1; i < num_of_elements + 1; i++) {
+    for (int j = 1; j < num_of_elements + 1; j++) {
+      //Check for knot_spacing type
+      //Currently only uniform is supported
+      if (UFBS2b[i][j].knot_spacing_type != 0)
+        error->all(FLERR,"UF3Kokkos: Currently only uniform knot-spacing is suupoted");
+    }
+  }*/
+
+  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 (UFBS3b[i][j][k].knot_spacing_type != 0)
+            error->all(FLERR,"UF3Kokkos: Currently only uniform knot-spacing is suupoted");*/
+          k_cut_3b.h_view(i,j,k) = cut_3b[i][j][k];
+
+          // Notice the order of min_cut_3b[i][j][k]
+          //In min_cut_3b[i][j][k],
+          //min_cut_3b[i][j][k][0] is the knot_vector along jk,
+          //min_cut_3b[i][j][k][1] is the knot_vector along ik,
+          //min_cut_3b[i][j][k][2] is the knot_vector along ij,
+          //see pair_uf3.cpp for more details
+          k_min_cut_3b.h_view(i,j,k,0) = min_cut_3b[i][j][k][0];
+          k_min_cut_3b.h_view(i,j,k,1) = min_cut_3b[i][j][k][1];
+          k_min_cut_3b.h_view(i,j,k,2) = min_cut_3b[i][j][k][2];
+          //k_cut_3b.h_view(i,k,j) = cut_3b[i][k][j];
+        }
+      }
+    }
+    k_cut_3b.template modify<LMPHostType>();
+    k_min_cut_3b.template modify<LMPHostType>();
+  }
+  //copy_2d(d_cutsq, cutsq, num_of_elements + 1, num_of_elements + 1); //copy cutsq from
+  //host to device memory
+  /*if (pot_3b) {
+    copy_3d(d_cut_3b, cut_3b, num_of_elements + 1, num_of_elements + 1, num_of_elements + 1);
+    copy_2d(d_cut_3b_list, cut_3b_list, num_of_elements + 1, num_of_elements + 1);
+  } else {
+    //why are we allocating space to d_cut_3b, d_cut_3b_list if the pot is 2-body only?
+    Kokkos::realloc(d_cut_3b, num_of_elements + 1, num_of_elements + 1, num_of_elements + 1);
+    Kokkos::realloc(d_cut_3b_list, num_of_elements + 1, num_of_elements + 1);
+  }*/
+  //No allocation for device equivalent of --> setflag, cut, knot_spacing_type_2b,
+  //n2b_knot, n2b_coeff, n2b_knot[i], n2b_coeff[i], setflag_3b, cut_3b, 
+  //cut_3b_list, min_cut_3b, knot_spacing_type_3b, cut_3b_list, n3b_knot_matrix,
+  //neighshort
+  
+  //UFBS2b and UFBS3b are array of objects. Bad idea to use kokkos view(array)
+  //for it
+  create_2b_coefficients();
+  if (pot_3b) create_3b_coefficients();
+
+}
+
+template <class DeviceType> void PairUF3Kokkos<DeviceType>::create_2b_coefficients()
+{
+
+  // Setup interaction pair map
+  //TODO: Instead of using map2b and map3b use simple indexing
+  Kokkos::realloc(map2b, num_of_elements + 1, num_of_elements + 1);
+  auto map2b_view = Kokkos::create_mirror(map2b);
+
+  int interaction_count = 0;
+  for (int i = 1; i < num_of_elements + 1; i++) {
+    for (int j = i; j < num_of_elements + 1; j++) {
+      map2b_view(i, j) = interaction_count;
+      map2b_view(j, i) = interaction_count++;
+    }
+  }
+  Kokkos::deep_copy(map2b, map2b_view);
+
+  // Count max knots for array size
+
+  int max_knots = 0;
+  for (int i = 1; i < n2b_knot.size(); i++)
+    for (int j = i; j < n2b_knot[i].size(); j++) max_knots = max(max_knots, n2b_knot[i][j].size());
+
+  // Copy coefficients to view
+
+  Kokkos::realloc(d_coefficients_2b, interaction_count, max_knots - 4);
+  auto d_coefficients_2b_view = Kokkos::create_mirror(d_coefficients_2b);
+
+  for (int i = 1; i < num_of_elements + 1; i++) {
+    for (int j = i; j < num_of_elements + 1; j++) {
+      for (int k = 0; k < n2b_coeff[i][j].size(); k++) {
+        d_coefficients_2b_view(map2b_view(i, j), k) = n2b_coeff[i][j][k];
+      }
+    }
+  }
+  Kokkos::deep_copy(d_coefficients_2b, d_coefficients_2b_view);
+
+  // Copy knots from array to view
+
+  Kokkos::realloc(d_n2b_knot, interaction_count, max_knots);
+  Kokkos::realloc(d_n2b_knot_spacings, interaction_count);
+  auto d_n2b_knot_view = Kokkos::create_mirror(d_n2b_knot);
+  auto d_n2b_knot_spacings_view = Kokkos::create_mirror(d_n2b_knot_spacings);
+
+  for (int i = 1; i < num_of_elements + 1; i++) {
+    for (int j = i; j < num_of_elements + 1; j++) {
+      for (int k = 0; k < n2b_knot[i][j].size(); k++) {
+        d_n2b_knot_view(map2b_view(i, j), k) = n2b_knot[i][j][k];
+      }
+      d_n2b_knot_spacings_view(map2b_view(i, j)) = UFBS2b[i][j].knot_spacing;
+    }
+  }
+
+  Kokkos::deep_copy(d_n2b_knot, d_n2b_knot_view);
+  Kokkos::deep_copy(d_n2b_knot_spacings, d_n2b_knot_spacings_view);
+  // Set spline constants
+
+  Kokkos::realloc(constants_2b, interaction_count, max_knots - 4);
+  auto constants_2b_view = Kokkos::create_mirror(constants_2b);
+
+  for (int i = 1; i < num_of_elements + 1; i++) {
+    for (int j = i; j < num_of_elements + 1; j++) {
+      for (int l = 0; l < n2b_knot[i][j].size() - 4; l++) {
+        auto c = get_constants(&n2b_knot[i][j][l], n2b_coeff[i][j][l]);
+        for (int k = 0; k < 16; k++)
+          constants_2b_view(map2b_view(i, j), l, k) = (std::isinf(c[k]) || std::isnan(c[k])) ? 0
+                                                                                             : c[k];
+      }
+    }
+  }
+  Kokkos::deep_copy(constants_2b, constants_2b_view);
+
+  Kokkos::realloc(dnconstants_2b, interaction_count, max_knots - 5);
+  auto dnconstants_2b_view = Kokkos::create_mirror(dnconstants_2b);
+
+  for (int i = 1; i < num_of_elements + 1; i++) {
+    for (int j = i; j < num_of_elements + 1; j++) {
+      for (int l = 0; l < n2b_knot[i][j].size() - 5; l++) {
+        double dntemp4 = 3 / (n2b_knot[i][j][l + 4] - n2b_knot[i][j][l + 1]);
+        double coeff = (n2b_coeff[i][j][l + 1] - n2b_coeff[i][j][l]) * dntemp4;
+        auto c = get_dnconstants(&n2b_knot[i][j][l + 1], coeff);
+        for (int k = 0; k < 9; k++)
+          dnconstants_2b_view(map2b_view(i, j), l, k) =
+              (std::isinf(c[k]) || std::isnan(c[k])) ? 0 : c[k];
+      }
+    }
+  }
+  Kokkos::deep_copy(dnconstants_2b, dnconstants_2b_view);
+}
+
+template <class DeviceType> void PairUF3Kokkos<DeviceType>::create_3b_coefficients()
+{
+  // Init interaction map for 3B
+
+  Kokkos::realloc(map3b, num_of_elements + 1, num_of_elements + 1, num_of_elements + 1);
+  auto map3b_view = Kokkos::create_mirror(map3b);
+
+  int interaction_count = 0;
+  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++) {
+        map3b_view(i, j, k) = interaction_count;
+        interaction_count++;
+       // map3b_view(i, k, j) = interaction_count++;
+      }
+    }
+  }
+  Kokkos::deep_copy(map3b, map3b_view);
+
+  // Count max knots for view
+
+  int max_knots = 0;
+  //In n3b_knot_matrix[i][j][k],
+  //n3b_knot_matrix[i][j][k][0] is the knot_vector along jk,
+  //n3b_knot_matrix[i][j][k][1] is the knot_vector along ik,
+  //n3b_knot_matrix[i][j][k][2] is the knot_vector along ij,
+  //see pair_uf3.cpp for more details
+
+  for (int i = 1; i < n3b_knot_matrix.size(); i++)
+    for (int j = 1; j < n3b_knot_matrix[i].size(); j++)
+      for (int k = 1; k < n3b_knot_matrix[i][j].size(); k++)
+        max_knots =
+            max(max_knots,
+                max(n3b_knot_matrix[i][j][k][0].size(),
+                    max(n3b_knot_matrix[i][j][k][1].size(), n3b_knot_matrix[i][j][k][2].size())));
+
+  // Init knot matrix view
+
+  Kokkos::realloc(d_n3b_knot_matrix, interaction_count, 3, max_knots);
+  Kokkos::realloc(d_n3b_knot_matrix_spacings, interaction_count, 3);
+  auto d_n3b_knot_matrix_view = Kokkos::create_mirror(d_n3b_knot_matrix);
+  auto d_n3b_knot_matrix_spacings_view = Kokkos::create_mirror(d_n3b_knot_matrix_spacings);
+
+  for (int i = 1; i < n3b_knot_matrix.size(); i++)
+    for (int j = 1; j < n3b_knot_matrix[i].size(); j++)
+      for (int k = 1; k < n3b_knot_matrix[i][j].size(); k++) {
+        for (int m = 0; m < n3b_knot_matrix[i][j][k][0].size(); m++)
+          d_n3b_knot_matrix_view(map3b_view(i, j, k), 0, m) = n3b_knot_matrix[i][j][k][0][m];
+        for (int m = 0; m < n3b_knot_matrix[i][j][k][1].size(); m++)
+          d_n3b_knot_matrix_view(map3b_view(i, j, k), 1, m) = n3b_knot_matrix[i][j][k][1][m];
+        for (int m = 0; m < n3b_knot_matrix[i][j][k][2].size(); m++)
+          d_n3b_knot_matrix_view(map3b_view(i, j, k), 2, m) = n3b_knot_matrix[i][j][k][2][m];
+
+        d_n3b_knot_matrix_spacings_view(map3b_view(i, j, k),2) = UFBS3b[i][j][k].knot_spacing_ij;
+        d_n3b_knot_matrix_spacings_view(map3b_view(i, j, k),1) = UFBS3b[i][j][k].knot_spacing_ik;
+        d_n3b_knot_matrix_spacings_view(map3b_view(i, j, k),0) = UFBS3b[i][j][k].knot_spacing_jk;
+      }
+  Kokkos::deep_copy(d_n3b_knot_matrix, d_n3b_knot_matrix_view);
+  Kokkos::deep_copy(d_n3b_knot_matrix_spacings, d_n3b_knot_matrix_spacings_view);
+
+  // Set knots spacings
+
+  Kokkos::realloc(d_n3b_knot_spacings, interaction_count, 3);
+  auto d_n3b_knot_spacings_view = Kokkos::create_mirror(d_n3b_knot_spacings);
+
+  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++) {
+        d_n3b_knot_spacings_view(map3b_view(i, j, k), 0) =
+            1 / (n3b_knot_matrix[i][j][k][0][5] - n3b_knot_matrix[i][j][k][0][4]);
+        d_n3b_knot_spacings_view(map3b_view(i, j, k), 1) =
+            1 / (n3b_knot_matrix[i][j][k][1][5] - n3b_knot_matrix[i][j][k][1][4]);
+        d_n3b_knot_spacings_view(map3b_view(i, j, k), 2) =
+            1 / (n3b_knot_matrix[i][j][k][2][5] - n3b_knot_matrix[i][j][k][2][4]);
+      }
+    }
+  }
+  Kokkos::deep_copy(d_n3b_knot_spacings, d_n3b_knot_spacings_view);
+
+  // Copy coefficients
+
+  Kokkos::realloc(d_coefficients_3b, interaction_count, max_knots - 4, max_knots - 4,
+                  max_knots - 4);
+  auto d_coefficients_3b_view = Kokkos::create_mirror(d_coefficients_3b);
+
+  for (int n = 1; n < num_of_elements + 1; n++) {
+    for (int m = 1; m < num_of_elements + 1; m++) {
+      for (int o = 1; o < num_of_elements + 1; o++) {
+        std::string key = std::to_string(n) + std::to_string(m) + std::to_string(o);
+        for (int i = 0; i < n3b_coeff_matrix[key].size(); i++) {
+          for (int j = 0; j < n3b_coeff_matrix[key][i].size(); j++) {
+            for (int k = 0; k < n3b_coeff_matrix[key][i][j].size() - 1; k++) {
+              d_coefficients_3b_view(map3b_view(n, m, o), i, j, k) = n3b_coeff_matrix[key][i][j][k];
+            }
+          }
+        }
+      }
+    }
+  }
+  Kokkos::deep_copy(d_coefficients_3b, d_coefficients_3b_view);
+  //
+  // Create derivative coefficients
+
+  // TODO: Shrink size
+  Kokkos::realloc(d_dncoefficients_3b, interaction_count, 3, max_knots - 4, max_knots - 4,
+                  max_knots - 4);
+  auto d_dncoefficients_3b_view = Kokkos::create_mirror(d_dncoefficients_3b);
+
+  //Notice the order for d_dncoefficients_3b_view(map3b_view(n, m, o), X, i, j, k)
+  //d_dncoefficients_3b_view(map3b_view(n, m, o), 2, i, j, k) --> coeff for rjk
+  //d_dncoefficients_3b_view(map3b_view(n, m, o), 1, i, j, k) --> coeff for rik
+  //d_dncoefficients_3b_view(map3b_view(n, m, o), 0, i, j, k) --> coeff for rij
+  //
+  //This is because-
+  //In n3b_knot_matrix[i][j][k],
+  //n3b_knot_matrix[i][j][k][0] is the knot_vector along jk,
+  //n3b_knot_matrix[i][j][k][1] is the knot_vector along ik,
+  //n3b_knot_matrix[i][j][k][2] is the knot_vector along ij,
+  //see pair_uf3.cpp for more details
+
+  for (int n = 1; n < num_of_elements + 1; n++) {
+    for (int m = 1; m < num_of_elements + 1; m++) {
+      for (int o = 1; o < num_of_elements + 1; o++) {
+        std::string key = std::to_string(n) + std::to_string(m) + std::to_string(o);
+        for (int i = 0; i < n3b_coeff_matrix[key].size(); i++) {
+          for (int j = 0; j < n3b_coeff_matrix[key][i].size(); j++) {
+            for (int k = 0; k < n3b_coeff_matrix[key][i][j].size() - 1; k++) {
+              F_FLOAT dntemp4 =
+                  3 / (n3b_knot_matrix[n][m][o][0][k + 4] - n3b_knot_matrix[n][m][o][0][k + 1]);
+              d_dncoefficients_3b_view(map3b_view(n, m, o), 2, i, j, k) =
+                  (n3b_coeff_matrix[key][i][j][k + 1] - n3b_coeff_matrix[key][i][j][k]) * dntemp4;
+            }
+          }
+        }
+
+        for (int i = 0; i < n3b_coeff_matrix[key].size(); i++) {
+          std::vector<std::vector<F_FLOAT>> dncoeff_vect2;
+          for (int j = 0; j < n3b_coeff_matrix[key][i].size() - 1; j++) {
+            F_FLOAT dntemp4 =
+                3 / (n3b_knot_matrix[n][m][o][1][j + 4] - n3b_knot_matrix[n][m][o][1][j + 1]);
+            std::vector<F_FLOAT> dncoeff_vect;
+            for (int k = 0; k < n3b_coeff_matrix[key][i][j].size(); k++) {
+              d_dncoefficients_3b_view(map3b_view(n, m, o), 1, i, j, k) =
+                  (n3b_coeff_matrix[key][i][j + 1][k] - n3b_coeff_matrix[key][i][j][k]) * dntemp4;
+            }
+          }
+        }
+
+        for (int i = 0; i < n3b_coeff_matrix[key].size() - 1; i++) {
+          F_FLOAT dntemp4 =
+              3 / (n3b_knot_matrix[n][m][o][2][i + 4] - n3b_knot_matrix[n][m][o][2][i + 1]);
+          for (int j = 0; j < n3b_coeff_matrix[key][i].size(); j++) {
+            for (int k = 0; k < n3b_coeff_matrix[key][i][j].size(); k++) {
+              d_dncoefficients_3b_view(map3b_view(n, m, o), 0, i, j, k) =
+                  (n3b_coeff_matrix[key][i + 1][j][k] - n3b_coeff_matrix[key][i][j][k]) * dntemp4;
+            }
+          }
+        }
+      }
+    }
+  }
+  Kokkos::deep_copy(d_dncoefficients_3b, d_dncoefficients_3b_view);
+
+  // Set spline constants
+
+  Kokkos::realloc(constants_3b, interaction_count, 3, max_knots - 4);
+  auto constants_3b_view = Kokkos::create_mirror(constants_3b);
+
+  //In n3b_knot_matrix[i][j][k],
+  //n3b_knot_matrix[i][j][k][0] is the knot_vector along jk,
+  //n3b_knot_matrix[i][j][k][1] is the knot_vector along ik,
+  //n3b_knot_matrix[i][j][k][2] is the knot_vector along ij,
+  //see pair_uf3.cpp for more details
+  for (int n = 1; n < num_of_elements + 1; n++) {
+    for (int m = 1; m < num_of_elements + 1; m++) {
+      for (int o = 1; o < num_of_elements + 1; o++) {
+        for (int l = 0; l < n3b_knot_matrix[n][m][o][2].size() - 4; l++) {
+          auto c = get_constants(&n3b_knot_matrix[n][m][o][2][l], 1);
+          for (int k = 0; k < 16; k++)
+            constants_3b_view(map3b_view(n, m, o), 0, l, k) =
+                (std::isinf(c[k]) || std::isnan(c[k])) ? 0 : c[k];
+        }
+        for (int l = 0; l < n3b_knot_matrix[n][m][o][1].size() - 4; l++) {
+          auto c = get_constants(&n3b_knot_matrix[n][m][o][1][l], 1);
+          for (int k = 0; k < 16; k++)
+            constants_3b_view(map3b_view(n, m, o), 1, l, k) =
+                (std::isinf(c[k]) || std::isnan(c[k])) ? 0 : c[k];
+        }
+        for (int l = 0; l < n3b_knot_matrix[n][m][o][0].size() - 4; l++) {
+          auto c = get_constants(&n3b_knot_matrix[n][m][o][0][l], 1);
+          for (int k = 0; k < 16; k++)
+            constants_3b_view(map3b_view(n, m, o), 2, l, k) =
+                (std::isinf(c[k]) || std::isnan(c[k])) ? 0 : c[k];
+        }
+      }
+    }
+  }
+  Kokkos::deep_copy(constants_3b, constants_3b_view);
+
+  Kokkos::realloc(dnconstants_3b, interaction_count, 3, max_knots - 6);
+  auto dnconstants_3b_view = Kokkos::create_mirror(dnconstants_3b);
+
+  for (int n = 1; n < num_of_elements + 1; n++) {
+    for (int m = 1; m < num_of_elements + 1; m++) {
+      for (int o = 1; o < num_of_elements + 1; o++) {
+        for (int l = 1; l < n3b_knot_matrix[n][m][o][2].size() - 5; l++) {
+          auto c = get_dnconstants(&n3b_knot_matrix[n][m][o][2][l], 1);
+          for (int k = 0; k < 9; k++)
+            dnconstants_3b_view(map3b_view(n, m, o), 0, l - 1, k) =
+                (std::isinf(c[k]) || std::isnan(c[k])) ? 0 : c[k];
+        }
+        for (int l = 1; l < n3b_knot_matrix[n][m][o][1].size() - 5; l++) {
+          auto c = get_dnconstants(&n3b_knot_matrix[n][m][o][1][l], 1);
+          for (int k = 0; k < 9; k++)
+            dnconstants_3b_view(map3b_view(n, m, o), 1, l - 1, k) =
+                (std::isinf(c[k]) || std::isnan(c[k])) ? 0 : c[k];
+        }
+        for (int l = 1; l < n3b_knot_matrix[n][m][o][0].size() - 5; l++) {
+          auto c = get_dnconstants(&n3b_knot_matrix[n][m][o][0][l], 1);
+          for (int k = 0; k < 9; k++)
+            dnconstants_3b_view(map3b_view(n, m, o), 2, l - 1, k) =
+                (std::isinf(c[k]) || std::isnan(c[k])) ? 0 : c[k];
+        }
+      }
+    }
+  }
+  Kokkos::deep_copy(dnconstants_3b, dnconstants_3b_view);
+}
+
+template <class DeviceType>
+template <int EVFLAG>
+KOKKOS_INLINE_FUNCTION void PairUF3Kokkos<DeviceType>::twobody(const int itype, const int jtype,
+                                                               const F_FLOAT r, F_FLOAT &evdwl,
+                                                               F_FLOAT &fpair) const
+{
+
+  // Find knot starting position
+  int interaction_id = map2b(itype, jtype);
+  int start_index = 3;
+  while (r > d_n2b_knot(interaction_id, start_index + 1)) start_index++;
+  //int start_index = 3+(int)((r-d_n2b_knot(interaction_id,0))/d_n2b_knot_spacings(interaction_id));
+
+  F_FLOAT r_values[4];
+  r_values[0] = 1;
+  r_values[1] = r;
+  r_values[2] = r_values[1] * r_values[1];
+
+  if (EVFLAG) {
+    r_values[3] = r_values[2] * r_values[1];
+    // Calculate energy
+    evdwl = constants_2b(interaction_id, start_index, 0);
+    evdwl += r_values[1] * constants_2b(interaction_id, start_index, 1);
+    evdwl += r_values[2] * constants_2b(interaction_id, start_index, 2);
+    evdwl += r_values[3] * constants_2b(interaction_id, start_index, 3);
+    evdwl += constants_2b(interaction_id, start_index - 1, 4);
+    evdwl += r_values[1] * constants_2b(interaction_id, start_index - 1, 5);
+    evdwl += r_values[2] * constants_2b(interaction_id, start_index - 1, 6);
+    evdwl += r_values[3] * constants_2b(interaction_id, start_index - 1, 7);
+    evdwl += constants_2b(interaction_id, start_index - 2, 8);
+    evdwl += r_values[1] * constants_2b(interaction_id, start_index - 2, 9);
+    evdwl += r_values[2] * constants_2b(interaction_id, start_index - 2, 10);
+    evdwl += r_values[3] * constants_2b(interaction_id, start_index - 2, 11);
+    evdwl += constants_2b(interaction_id, start_index - 3, 12);
+    evdwl += r_values[1] * constants_2b(interaction_id, start_index - 3, 13);
+    evdwl += r_values[2] * constants_2b(interaction_id, start_index - 3, 14);
+    evdwl += r_values[3] * constants_2b(interaction_id, start_index - 3, 15);
+  }
+
+  // Calculate force
+  fpair = dnconstants_2b(interaction_id, start_index - 1, 0);
+  fpair += r_values[1] * dnconstants_2b(interaction_id, start_index - 1, 1);
+  fpair += r_values[2] * dnconstants_2b(interaction_id, start_index - 1, 2);
+  fpair += dnconstants_2b(interaction_id, start_index - 2, 3);
+  fpair += r_values[1] * dnconstants_2b(interaction_id, start_index - 2, 4);
+  fpair += r_values[2] * dnconstants_2b(interaction_id, start_index - 2, 5);
+  fpair += dnconstants_2b(interaction_id, start_index - 3, 6);
+  fpair += r_values[1] * dnconstants_2b(interaction_id, start_index - 3, 7);
+  fpair += r_values[2] * dnconstants_2b(interaction_id, start_index - 3, 8);
+}
+
+template <class DeviceType>
+template <int EVFLAG>
+KOKKOS_INLINE_FUNCTION void PairUF3Kokkos<DeviceType>::threebody(
+    const int itype, const int jtype, const int ktype, const F_FLOAT value_rij,
+    const F_FLOAT value_rik, const F_FLOAT value_rjk, F_FLOAT &evdwl, F_FLOAT (&fforce)[3]) const
+{
+  evdwl = 0;
+  fforce[0] = 0;
+  fforce[1] = 0;
+  fforce[2] = 0;
+
+  F_FLOAT evals[3][4];
+  F_FLOAT dnevals[3][4];
+  int start_indices[3];
+  F_FLOAT r[3] = {value_rij, value_rik, value_rjk};
+  int interaction_id = map3b(itype, jtype, ktype);
+
+  auto coefficients =
+      Kokkos::subview(d_coefficients_3b, interaction_id, Kokkos::ALL, Kokkos::ALL, Kokkos::ALL);
+  auto dncoefficients = Kokkos::subview(d_dncoefficients_3b, interaction_id, Kokkos::ALL,
+                                        Kokkos::ALL, Kokkos::ALL, Kokkos::ALL);
+  //Notice the 2-d in d_n3b_knot_matrix
+  //
+  //In d_n3b_knot_matrix[i][j][k],
+  //d_n3b_knot_matrix[i][j][k][0] is the knot_vector along jk,
+  //d_n3b_knot_matrix[i][j][k][1] is the knot_vector along ik,
+  //d_n3b_knot_matrix[i][j][k][2] is the knot_vector along ij,
+  //
+  //and r[0] = rij, r[1] = rik and r[2] = rjk
+  //see n3b_knot_matrix and pair_uf3.cpp for more details
+  for (int d = 0; d < 3; d++) {
+    start_indices[d] = 3;
+    while (r[d] > d_n3b_knot_matrix(interaction_id, 2-d, start_indices[d] + 1)) start_indices[d]++;
+    //start_indices[d] = 3+(int)((r[d]-d_n3b_knot_matrix(interaction_id, 2-d, 0))/d_n3b_knot_matrix_spacings(interaction_id, 2-d));
+
+    F_FLOAT r_values[4];
+    r_values[0] = 1;
+    r_values[1] = r[d];
+    r_values[2] = r_values[1] * r_values[1];
+
+    r_values[3] = r_values[2] * r_values[1];
+
+    // Calculate energy
+    evals[d][0] = constants_3b(interaction_id, d, start_indices[d], 0);
+    evals[d][0] += r_values[1] * constants_3b(interaction_id, d, start_indices[d], 1);
+    evals[d][0] += r_values[2] * constants_3b(interaction_id, d, start_indices[d], 2);
+    evals[d][0] += r_values[3] * constants_3b(interaction_id, d, start_indices[d], 3);
+    evals[d][1] = constants_3b(interaction_id, d, start_indices[d] - 1, 4);
+    evals[d][1] += r_values[1] * constants_3b(interaction_id, d, start_indices[d] - 1, 5);
+    evals[d][1] += r_values[2] * constants_3b(interaction_id, d, start_indices[d] - 1, 6);
+    evals[d][1] += r_values[3] * constants_3b(interaction_id, d, start_indices[d] - 1, 7);
+    evals[d][2] = constants_3b(interaction_id, d, start_indices[d] - 2, 8);
+    evals[d][2] += r_values[1] * constants_3b(interaction_id, d, start_indices[d] - 2, 9);
+    evals[d][2] += r_values[2] * constants_3b(interaction_id, d, start_indices[d] - 2, 10);
+    evals[d][2] += r_values[3] * constants_3b(interaction_id, d, start_indices[d] - 2, 11);
+    evals[d][3] = constants_3b(interaction_id, d, start_indices[d] - 3, 12);
+    evals[d][3] += r_values[1] * constants_3b(interaction_id, d, start_indices[d] - 3, 13);
+    evals[d][3] += r_values[2] * constants_3b(interaction_id, d, start_indices[d] - 3, 14);
+    evals[d][3] += r_values[3] * constants_3b(interaction_id, d, start_indices[d] - 3, 15);
+
+    dnevals[d][0] = dnconstants_3b(interaction_id, d, start_indices[d] - 1, 0);
+    dnevals[d][0] += r_values[1] * dnconstants_3b(interaction_id, d, start_indices[d] - 1, 1);
+    dnevals[d][0] += r_values[2] * dnconstants_3b(interaction_id, d, start_indices[d] - 1, 2);
+    dnevals[d][1] = dnconstants_3b(interaction_id, d, start_indices[d] - 2, 3);
+    dnevals[d][1] += r_values[1] * dnconstants_3b(interaction_id, d, start_indices[d] - 2, 4);
+    dnevals[d][1] += r_values[2] * dnconstants_3b(interaction_id, d, start_indices[d] - 2, 5);
+    dnevals[d][2] = dnconstants_3b(interaction_id, d, start_indices[d] - 3, 6);
+    dnevals[d][2] += r_values[1] * dnconstants_3b(interaction_id, d, start_indices[d] - 3, 7);
+    dnevals[d][2] += r_values[2] * dnconstants_3b(interaction_id, d, start_indices[d] - 3, 8);
+    dnevals[d][3] = 0;
+  }
+
+  if (EVFLAG) {
+    for (int i = 0; i < 4; i++) {
+      for (int j = 0; j < 4; j++) {
+        for (int k = 0; k < 4; k++) {
+          evdwl += coefficients(start_indices[0] - i, start_indices[1] - j, start_indices[2] - k) *
+              evals[0][i] * evals[1][j] * evals[2][k];
+        }
+      }
+    }
+  }
+
+  for (int i = 0; i < 3; i++) {
+    for (int j = 0; j < 4; j++) {
+      for (int k = 0; k < 4; k++) {
+        fforce[0] += dncoefficients(0, start_indices[0] - 3 + i, start_indices[1] - 3 + j,
+                                    start_indices[2] - 3 + k) *
+            dnevals[0][2 - i] * evals[1][3 - j] * evals[2][3 - k];
+      }
+    }
+  }
+
+  for (int i = 0; i < 4; i++) {
+    for (int j = 0; j < 3; j++) {
+      for (int k = 0; k < 4; k++) {
+        fforce[1] += dncoefficients(1, start_indices[0] - 3 + i, start_indices[1] - 3 + j,
+                                    start_indices[2] - 3 + k) *
+            evals[0][3 - i] * dnevals[1][2 - j] * evals[2][3 - k];
+      }
+    }
+  }
+
+  for (int i = 0; i < 4; i++) {
+    for (int j = 0; j < 4; j++) {
+      for (int k = 0; k < 3; k++) {
+        fforce[2] += dncoefficients(2, start_indices[0] - 3 + i, start_indices[1] - 3 + j,
+                                    start_indices[2] - 3 + k) *
+            evals[0][3 - i] * evals[1][3 - j] * dnevals[2][2 - k];
+      }
+    }
+  }
+}
+
+/* ----------------------------------------------------------------------
+   init specific to this pair style
+------------------------------------------------------------------------- */
+
+/*template <class DeviceType> void PairUF3Kokkos<DeviceType>::init_style()
+{
+
+  PairUF3::init_style();
+
+  neighflag = lmp->kokkos->neighflag;
+
+  auto request = neighbor->find_request(this);
+  request->set_kokkos_host(std::is_same<DeviceType, LMPHostType>::value &&
+                           !std::is_same<DeviceType, LMPDeviceType>::value);
+  request->set_kokkos_device(std::is_same<DeviceType, LMPDeviceType>::value);
+
+  request->enable_full();
+  // request->enable_ghost();
+}*/
+
+/* ----------------------------------------------------------------------
+   init list sets the pointer to full neighbour list requested in previous function
+------------------------------------------------------------------------- */
+
+//template <class DeviceType>
+//void PairUF3Kokkos<DeviceType>::init_list(int /*id*/, class NeighList *ptr)
+//{
+//  list = ptr;
+//}
+
+template <class DeviceType> void PairUF3Kokkos<DeviceType>::compute(int eflag_in, int vflag_in)
+{
+  eflag = eflag_in;
+  vflag = vflag_in;
+
+  if (neighflag == FULL) no_virial_fdotr_compute = 1;
+
+  ev_init(eflag, vflag, 0);
+
+  // reallocate per-atom arrays if necessary
+
+  if (eflag_atom) {
+    memoryKK->destroy_kokkos(k_eatom, eatom);
+    memoryKK->create_kokkos(k_eatom, eatom, maxeatom, "pair:eatom");
+    d_eatom = k_eatom.view<DeviceType>();
+  }
+  if (vflag_atom) {
+    memoryKK->destroy_kokkos(k_vatom, vatom);
+    memoryKK->create_kokkos(k_vatom, vatom, maxvatom, "pair:vatom");
+    d_vatom = k_vatom.view<DeviceType>();
+  }
+
+  atomKK->sync(execution_space, datamask_read);
+  if (eflag || vflag) atomKK->modified(execution_space,datamask_modify);
+  else atomKK->modified(execution_space,F_MASK);
+
+  x = atomKK->k_x.template view<DeviceType>();
+  f = atomKK->k_f.template view<DeviceType>();
+  tag = atomKK->k_tag.template view<DeviceType>();
+  type = atomKK->k_type.template view<DeviceType>();
+  nlocal = atom->nlocal;
+  newton_pair = force->newton_pair;
+  nall = atom->nlocal + atom->nghost;
+  k_cutsq.template sync<DeviceType>(); //Sync the device memory of k_cutsq with
+  //the array from the host memory; this updates d_cutsq also
+  k_cut_3b.template sync<DeviceType>();
+  k_min_cut_3b.template sync<DeviceType>();
+  
+  inum = list->inum;
+  const int ignum = inum + list->gnum;
+  NeighListKokkos<DeviceType> *k_list = static_cast<NeighListKokkos<DeviceType> *>(list);
+  d_ilist = k_list->d_ilist;
+  d_numneigh = k_list->d_numneigh;
+  d_neighbors = k_list->d_neighbors;
+
+  copymode = 1;
+
+  escatter = ScatterEType(d_eatom);
+  fscatter = ScatterFType(f);
+  vscatter = ScatterVType(d_vatom);
+  //cvscatter = ScatterCVType(d_cvatom);
+
+  EV_FLOAT ev;
+  EV_FLOAT ev_all;
+
+  // build short neighbor list
+
+  int max_neighs = d_neighbors.extent(1);
+
+  if (((int)d_neighbors_short.extent(1) != max_neighs) || 
+          ((int)d_neighbors_short.extent(0) != ignum)) {
+    d_neighbors_short = Kokkos::View<int **, DeviceType>("UF3::neighbors_short", ignum, max_neighs);
+  }
+  if (d_numneigh_short.extent(0) != ignum)
+    d_numneigh_short = Kokkos::View<int *, DeviceType>("UF3::numneighs_short", ignum);
+  Kokkos::parallel_for(
+      Kokkos::RangePolicy<DeviceType, TagPairUF3ComputeShortNeigh>(0, ignum), *this);
+
+  // loop over neighbor list of my atoms
+
+  if (evflag){
+    Kokkos::parallel_reduce(
+        Kokkos::RangePolicy<DeviceType, TagPairUF3ComputeFullA<FULL, 1>>(0, inum), *this, ev);
+  }
+  else{
+    Kokkos::parallel_for(
+        Kokkos::RangePolicy<DeviceType, TagPairUF3ComputeFullA<FULL, 0>>(0, inum), *this);
+  }
+  ev_all += ev;
+
+  Kokkos::Experimental::contribute(d_eatom, escatter);
+  Kokkos::Experimental::contribute(d_vatom, vscatter);
+  //Kokkos::Experimental::contribute(d_cvatom, cvscatter);
+  Kokkos::Experimental::contribute(f, fscatter);
+
+  if (eflag_global) eng_vdwl += ev_all.evdwl;
+  if (vflag_global) {
+    virial[0] += ev_all.v[0];
+    virial[1] += ev_all.v[1];
+    virial[2] += ev_all.v[2];
+    virial[3] += ev_all.v[3];
+    virial[4] += ev_all.v[4];
+    virial[5] += ev_all.v[5];
+  }
+
+  if (eflag_atom) {
+    k_eatom.template modify<DeviceType>();
+    k_eatom.template sync<LMPHostType>();
+  }
+
+  if (vflag_atom) {
+    k_vatom.template modify<DeviceType>();
+    k_vatom.template sync<LMPHostType>();
+  }
+
+  if (cvflag_atom) {
+    //k_cvatom.template modify<DeviceType>();
+    //k_cvatom.template sync<LMPHostType>();
+  }
+
+  if (vflag_fdotr) pair_virial_fdotr_compute(this);
+
+  copymode = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+template <class DeviceType>
+KOKKOS_INLINE_FUNCTION void PairUF3Kokkos<DeviceType>::operator()(TagPairUF3ComputeShortNeigh,
+                                                                  const int &ii) const
+{
+  const int i = d_ilist[ii];
+  const X_FLOAT xtmp = x(i, 0);
+  const X_FLOAT ytmp = x(i, 1);
+  const X_FLOAT ztmp = x(i, 2);
+
+  const int jnum = d_numneigh[i];
+  int inside = 0;
+  for (int jj = 0; jj < jnum; jj++) {
+    int j = d_neighbors(i, jj);
+    j &= NEIGHMASK;
+
+    const X_FLOAT delx = xtmp - x(j, 0);
+    const X_FLOAT dely = ytmp - x(j, 1);
+    const X_FLOAT delz = ztmp - x(j, 2);
+    const F_FLOAT rsq = delx * delx + dely * dely + delz * delz;
+
+    const int itype = type[i];
+    const int jtype = type[j];
+
+    if (rsq <= d_cutsq(itype, jtype)) {
+      d_neighbors_short(i, inside) = j;
+      inside++;
+    }
+  }
+  d_numneigh_short(i) = inside;
+}
+
+/* ---------------------------------------------------------------------- */
+
+template <class DeviceType>
+template <int NEIGHFLAG, int EVFLAG>
+KOKKOS_INLINE_FUNCTION void
+PairUF3Kokkos<DeviceType>::operator()(TagPairUF3ComputeFullA<NEIGHFLAG, EVFLAG>, const int &ii,
+                                      EV_FLOAT &ev) const
+{
+  // The f array is duplicated for OpenMP, atomic for CUDA, and neither for Serial
+
+  auto v_f = vscatter.access();
+  auto a_f = fscatter.access();
+  auto a_cvatom = cvscatter.access();
+
+  F_FLOAT del_rji[3], del_rki[3], del_rkj[3], triangle_eval[3];
+  F_FLOAT fij[3], fik[3], fjk[3];
+  F_FLOAT fji[3], fki[3], fkj[3];
+  F_FLOAT Fj[3], Fk[3];
+  F_FLOAT evdwl = 0, evdwl3 = 0;
+  F_FLOAT fpair = 0;
+
+  const int i = d_ilist[ii];
+
+  const tagint itag = tag[i];
+  const int itype = type[i];
+  const X_FLOAT xtmp = x(i, 0);
+  const X_FLOAT ytmp = x(i, 1);
+  const X_FLOAT ztmp = x(i, 2);
+
+  // two-body interactions
+
+  const int jnum = d_numneigh_short[i];
+
+  F_FLOAT fxtmpi = 0.0;
+  F_FLOAT fytmpi = 0.0;
+  F_FLOAT fztmpi = 0.0;
+
+  for (int jj = 0; jj < jnum; jj++) {
+    int j = d_neighbors_short(i, jj);
+    j &= NEIGHMASK;
+    const tagint jtag = tag[j];
+
+    const int jtype = type[j];
+    
+    const X_FLOAT delx = xtmp - x(j, 0);
+    const X_FLOAT dely = ytmp - x(j, 1);
+    const X_FLOAT delz = ztmp - x(j, 2);
+    const F_FLOAT rsq = delx * delx + dely * dely + delz * delz;
+
+    if (rsq >= d_cutsq(itype, jtype)) continue;
+
+    const F_FLOAT rij = sqrt(rsq);
+    this->template twobody<EVFLAG>(itype, jtype, rij, evdwl, fpair);
+    
+    fpair = -fpair / rij;
+
+    fxtmpi += delx * fpair;
+    fytmpi += dely * fpair;
+    fztmpi += delz * fpair;
+    a_f(j, 0) -= delx * fpair;
+    a_f(j, 1) -= dely * fpair;
+    a_f(j, 2) -= delz * fpair;
+
+    if (EVFLAG) {
+      if (eflag) ev.evdwl += evdwl;
+      if (vflag_either || eflag_atom)
+        this->template ev_tally<NEIGHFLAG>(ev, i, j, evdwl, fpair, delx, dely, delz);
+    }
+  }
+
+  // 3-body interaction
+  // jth atom
+  const int jnumm1 = jnum - 1;
+  for (int jj = 0; jj < jnumm1; jj++) {
+    int j = d_neighbors_short(i, jj);
+    j &= NEIGHMASK;
+    const int jtype = type[j];
+    del_rji[0] = x(j, 0) - xtmp;
+    del_rji[1] = x(j, 1) - ytmp;
+    del_rji[2] = x(j, 2) - ztmp;
+    F_FLOAT rij = sqrt(del_rji[0] * del_rji[0] + del_rji[1] * del_rji[1] + del_rji[2] * del_rji[2]);
+
+    F_FLOAT fxtmpj = 0.0;
+    F_FLOAT fytmpj = 0.0;
+    F_FLOAT fztmpj = 0.0;
+
+    for (int kk = jj + 1; kk < jnum; kk++) {
+      int k = d_neighbors_short(i, kk);
+      k &= NEIGHMASK;
+      const int ktype = type[k];
+
+      // Notice the order of d_min_cut_3b[i][j][k]
+      //In d_min_cut_3b[i][j][k],
+      //d_min_cut_3b[i][j][k][0] is the knot_vector along jk,
+      //d_min_cut_3b[i][j][k][1] is the knot_vector along ik,
+      //d_min_cut_3b[i][j][k][2] is the knot_vector along ij,
+      //see pair_uf3.cpp for more details
+      if (rij < d_min_cut_3b(itype, jtype, ktype, 2)) continue;
+      if (rij > d_cut_3b(itype, jtype, ktype)) continue;
+
+      del_rki[0] = x(k, 0) - xtmp;
+      del_rki[1] = x(k, 1) - ytmp;
+      del_rki[2] = x(k, 2) - ztmp;
+      F_FLOAT rik =
+          sqrt(del_rki[0] * del_rki[0] + del_rki[1] * del_rki[1] + del_rki[2] * del_rki[2]);
+
+      if (rik < d_min_cut_3b(itype, jtype, ktype, 1)) continue;
+      if (rik > d_cut_3b(itype, ktype, jtype)) continue;
+
+      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);
+      F_FLOAT rjk =
+          sqrt(del_rkj[0] * del_rkj[0] + del_rkj[1] * del_rkj[1] + del_rkj[2] * del_rkj[2]);
+      if (rjk < d_min_cut_3b(itype, jtype, ktype, 0)) continue;
+      this->template threebody<EVFLAG>(itype, jtype, ktype, rij, rik, rjk, evdwl3, triangle_eval);
+
+      fij[0] = *(triangle_eval + 0) * (del_rji[0] / rij);
+      fji[0] = -fij[0];
+      fik[0] = *(triangle_eval + 1) * (del_rki[0] / rik);
+      fki[0] = -fik[0];
+      fjk[0] = *(triangle_eval + 2) * (del_rkj[0] / rjk);
+      fkj[0] = -fjk[0];
+
+      fij[1] = *(triangle_eval + 0) * (del_rji[1] / rij);
+      fji[1] = -fij[1];
+      fik[1] = *(triangle_eval + 1) * (del_rki[1] / rik);
+      fki[1] = -fik[1];
+      fjk[1] = *(triangle_eval + 2) * (del_rkj[1] / rjk);
+      fkj[1] = -fjk[1];
+
+      fij[2] = *(triangle_eval + 0) * (del_rji[2] / rij);
+      fji[2] = -fij[2];
+      fik[2] = *(triangle_eval + 1) * (del_rki[2] / rik);
+      fki[2] = -fik[2];
+      fjk[2] = *(triangle_eval + 2) * (del_rkj[2] / rjk);
+      fkj[2] = -fjk[2];
+
+      Fj[0] = fji[0] + fjk[0];
+      Fj[1] = fji[1] + fjk[1];
+      Fj[2] = fji[2] + fjk[2];
+
+      Fk[0] = fki[0] + fkj[0];
+      Fk[1] = fki[1] + fkj[1];
+      Fk[2] = fki[2] + fkj[2];
+
+      fxtmpi += (fij[0] + fik[0]);
+      fytmpi += (fij[1] + fik[1]);
+      fztmpi += (fij[2] + fik[2]);
+      fxtmpj += Fj[0];
+      fytmpj += Fj[1];
+      fztmpj += Fj[2];
+      a_f(k, 0) += Fk[0];
+      a_f(k, 1) += Fk[1];
+      a_f(k, 2) += Fk[2];
+
+      if (EVFLAG) {
+        if (eflag) { ev.evdwl += evdwl3; }
+        if (vflag_either || eflag_atom) {
+          this->template ev_tally3<NEIGHFLAG>(ev, i, j, k, evdwl3, 0.0, Fj, Fk, del_rji, del_rki);
+          if (cvflag_atom) {
+
+            F_FLOAT 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]);
+            a_cvatom(i, 0) += ric[0] * (-Fj[0] - Fk[0]);
+            a_cvatom(i, 1) += ric[1] * (-Fj[1] - Fk[1]);
+            a_cvatom(i, 2) += ric[2] * (-Fj[2] - Fk[2]);
+            a_cvatom(i, 3) += ric[0] * (-Fj[1] - Fk[1]);
+            a_cvatom(i, 4) += ric[0] * (-Fj[2] - Fk[2]);
+            a_cvatom(i, 5) += ric[1] * (-Fj[2] - Fk[2]);
+            a_cvatom(i, 6) += ric[1] * (-Fj[0] - Fk[0]);
+            a_cvatom(i, 7) += ric[2] * (-Fj[0] - Fk[0]);
+            a_cvatom(i, 8) += ric[2] * (-Fj[1] - Fk[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]);
+
+            a_cvatom(j, 0) += rjc[0] * Fj[0];
+            a_cvatom(j, 1) += rjc[1] * Fj[1];
+            a_cvatom(j, 2) += rjc[2] * Fj[2];
+            a_cvatom(j, 3) += rjc[0] * Fj[1];
+            a_cvatom(j, 4) += rjc[0] * Fj[2];
+            a_cvatom(j, 5) += rjc[1] * Fj[2];
+            a_cvatom(j, 6) += rjc[1] * Fj[0];
+            a_cvatom(j, 7) += rjc[2] * Fj[0];
+            a_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]);
+
+            a_cvatom(k, 0) += rkc[0] * Fk[0];
+            a_cvatom(k, 1) += rkc[1] * Fk[1];
+            a_cvatom(k, 2) += rkc[2] * Fk[2];
+            a_cvatom(k, 3) += rkc[0] * Fk[1];
+            a_cvatom(k, 4) += rkc[0] * Fk[2];
+            a_cvatom(k, 5) += rkc[1] * Fk[2];
+            a_cvatom(k, 6) += rkc[1] * Fk[0];
+            a_cvatom(k, 7) += rkc[2] * Fk[0];
+            a_cvatom(k, 8) += rkc[2] * Fk[1];
+          }
+        }
+      }
+    }
+    a_f(j, 0) += fxtmpj;
+    a_f(j, 1) += fytmpj;
+    a_f(j, 2) += fztmpj;
+  }
+
+  a_f(i, 0) += fxtmpi;
+  a_f(i, 1) += fytmpi;
+  a_f(i, 2) += fztmpi;
+}
+
+template <class DeviceType>
+template <int NEIGHFLAG, int EVFLAG>
+KOKKOS_INLINE_FUNCTION void
+PairUF3Kokkos<DeviceType>::operator()(TagPairUF3ComputeFullA<NEIGHFLAG, EVFLAG>,
+                                      const int &ii) const
+{
+  EV_FLOAT ev;
+  this->template operator()<NEIGHFLAG, EVFLAG>(TagPairUF3ComputeFullA<NEIGHFLAG, EVFLAG>(), ii, ev);
+}
+
+/* ---------------------------------------------------------------------- */
+
+template <class DeviceType>
+template <int NEIGHFLAG>
+KOKKOS_INLINE_FUNCTION void
+PairUF3Kokkos<DeviceType>::ev_tally(EV_FLOAT &ev, const int &i, const int &j, const F_FLOAT &epair,
+                                    const F_FLOAT &fpair, const F_FLOAT &delx, const F_FLOAT &dely,
+                                    const F_FLOAT &delz) const
+{
+
+  // The eatom and vatom arrays are duplicated for OpenMP, atomic for CUDA,
+  // and neither for Serial
+
+  auto a_eatom = escatter.access();
+  auto a_vatom = vscatter.access();
+  auto a_cvatom = cvscatter.access();
+
+  if (eflag_atom) {
+    const E_FLOAT epairhalf = 0.5 * epair;
+    a_eatom[i] += epairhalf;
+    a_eatom[j] += epairhalf;
+  }
+
+  if (vflag_either) {
+    const E_FLOAT v0 = delx * delx * fpair;
+    const E_FLOAT v1 = dely * dely * fpair;
+    const E_FLOAT v2 = delz * delz * fpair;
+    const E_FLOAT v3 = delx * dely * fpair;
+    const E_FLOAT v4 = delx * delz * fpair;
+    const E_FLOAT v5 = dely * delz * fpair;
+
+    if (vflag_global) {
+      ev.v[0] += v0;
+      ev.v[1] += v1;
+      ev.v[2] += v2;
+      ev.v[3] += v3;
+      ev.v[4] += v4;
+      ev.v[5] += v5;
+    }
+
+    if (vflag_atom) {
+      a_vatom(i, 0) += 0.5 * v0;
+      a_vatom(i, 1) += 0.5 * v1;
+      a_vatom(i, 2) += 0.5 * v2;
+      a_vatom(i, 3) += 0.5 * v3;
+      a_vatom(i, 4) += 0.5 * v4;
+      a_vatom(i, 5) += 0.5 * v5;
+
+      a_vatom(j, 0) += 0.5 * v0;
+      a_vatom(j, 1) += 0.5 * v1;
+      a_vatom(j, 2) += 0.5 * v2;
+      a_vatom(j, 3) += 0.5 * v3;
+      a_vatom(j, 4) += 0.5 * v4;
+      a_vatom(j, 5) += 0.5 * v5;
+    }
+
+    if (cvflag_atom) {
+      a_cvatom(i, 0) += 0.5 * v0;
+      a_cvatom(i, 1) += 0.5 * v1;
+      a_cvatom(i, 2) += 0.5 * v2;
+      a_cvatom(i, 3) += 0.5 * v3;
+      a_cvatom(i, 4) += 0.5 * v4;
+      a_cvatom(i, 5) += 0.5 * v5;
+      a_cvatom(i, 6) += 0.5 * v3;
+      a_cvatom(i, 7) += 0.5 * v4;
+      a_cvatom(i, 8) += 0.5 * v5;
+      a_cvatom(j, 0) += 0.5 * v0;
+      a_cvatom(j, 1) += 0.5 * v1;
+      a_cvatom(j, 2) += 0.5 * v2;
+      a_cvatom(j, 3) += 0.5 * v3;
+      a_cvatom(j, 4) += 0.5 * v4;
+      a_cvatom(j, 5) += 0.5 * v5;
+      a_cvatom(j, 6) += 0.5 * v3;
+      a_cvatom(j, 7) += 0.5 * v4;
+      a_cvatom(j, 8) += 0.5 * v5;
+    }
+  }
+}
+
+/* ----------------------------------------------------------------------
+   tally eng_vdwl and virial into global and per-atom accumulators
+   called by SW and hbond potentials, newton_pair is always on
+   virial = riFi + rjFj + rkFk = (rj-ri) Fj + (rk-ri) Fk = drji*fj + drki*fk
+ ------------------------------------------------------------------------- */
+
+template <class DeviceType>
+template <int NEIGHFLAG>
+KOKKOS_INLINE_FUNCTION void
+PairUF3Kokkos<DeviceType>::ev_tally3(EV_FLOAT &ev, const int &i, const int &j, int &k,
+                                     const F_FLOAT &evdwl, const F_FLOAT &ecoul, F_FLOAT *fj,
+                                     F_FLOAT *fk, F_FLOAT *drji, F_FLOAT *drki) const
+{
+  F_FLOAT epairthird, v[6];
+
+  // The eatom and vatom arrays are duplicated for OpenMP, atomic for CUDA,
+  // and neither for Serial
+
+  auto a_eatom = escatter.access();
+  auto a_vatom = vscatter.access();
+
+  if (eflag_atom) {
+    epairthird = THIRD * (evdwl + ecoul);
+    a_eatom[i] += epairthird;
+    a_eatom[j] += epairthird;
+    a_eatom[k] += epairthird;
+  }
+
+  if (vflag_either) {
+    v[0] = drji[0] * fj[0] + drki[0] * fk[0];
+    v[1] = drji[1] * fj[1] + drki[1] * fk[1];
+    v[2] = drji[2] * fj[2] + drki[2] * fk[2];
+    v[3] = drji[0] * fj[1] + drki[0] * fk[1];
+    v[4] = drji[0] * fj[2] + drki[0] * fk[2];
+    v[5] = drji[1] * fj[2] + drki[1] * fk[2];
+
+    if (vflag_global) {
+      ev.v[0] += v[0];
+      ev.v[1] += v[1];
+      ev.v[2] += v[2];
+      ev.v[3] += v[3];
+      ev.v[4] += v[4];
+      ev.v[5] += v[5];
+    }
+
+    if (vflag_atom) {
+      a_vatom(i, 0) += THIRD * v[0];
+      a_vatom(i, 1) += THIRD * v[1];
+      a_vatom(i, 2) += THIRD * v[2];
+      a_vatom(i, 3) += THIRD * v[3];
+      a_vatom(i, 4) += THIRD * v[4];
+      a_vatom(i, 5) += THIRD * v[5];
+
+      a_vatom(j, 0) += THIRD * v[0];
+      a_vatom(j, 1) += THIRD * v[1];
+      a_vatom(j, 2) += THIRD * v[2];
+      a_vatom(j, 3) += THIRD * v[3];
+      a_vatom(j, 4) += THIRD * v[4];
+      a_vatom(j, 5) += THIRD * v[5];
+
+      a_vatom(k, 0) += THIRD * v[0];
+      a_vatom(k, 1) += THIRD * v[1];
+      a_vatom(k, 2) += THIRD * v[2];
+      a_vatom(k, 3) += THIRD * v[3];
+      a_vatom(k, 4) += THIRD * v[4];
+      a_vatom(k, 5) += THIRD * v[5];
+    }
+  }
+}
+
+/* ----------------------------------------------------------------------
+   tally eng_vdwl and virial into global and per-atom accumulators
+   called by SW and hbond potentials, newton_pair is always on
+   virial = riFi + rjFj + rkFk = (rj-ri) Fj + (rk-ri) Fk = drji*fj + drki*fk
+ ------------------------------------------------------------------------- */
+
+template <class DeviceType>
+template <typename T, typename V>
+void PairUF3Kokkos<DeviceType>::copy_2d(V &d, T **h, int m, int n)
+{
+  Kokkos::View<T **> tmp("pair::tmp", m, n); //Create tmp view(array) on
+  //device memory
+
+  //auto h_view = Kokkos::create_mirror_view(tmp);
+  auto h_view = Kokkos::create_mirror(tmp); //Create a mirror of the device
+  //view(array) tmp, as deep_copy is only possible for mirror views
+
+  for (int i = 0; i < m; i++) {
+    for (int j = 0; j < n; j++) {
+      h_view(i, j) = h[i][j]; //fill mirror
+    }
+    //views with data from normal array 'h' which always lives on host memory
+  }
+
+  Kokkos::deep_copy(tmp, h_view); //Deepcopy data from h_view(host) to tmp(device)
+
+  d = tmp;
+}
+
+template <class DeviceType>
+template <typename T, typename V>
+void PairUF3Kokkos<DeviceType>::copy_3d(V &d, T ***h, int m, int n, int o)
+{
+  Kokkos::View<T ***> tmp("pair::tmp", m, n, o); //Create tmp view(array) on
+  //device memory
+
+  //auto h_view = Kokkos::create_mirror_view(tmp); //create_mirror always copies
+  //the data. create_mirror_view only copies data if the host cannot access the 
+  //data
+  auto h_view = Kokkos::create_mirror(tmp); //Create a mirror of the device
+  //view(array) tmp, as deep_copy is only possible for mirror views
+  for (int i = 0; i < m; i++) {
+    for (int j = 0; j < n; j++) {
+      for (int k = 0; k < o; k++) { h_view(i, j, k) = h[i][j][k]; } //fill mirror
+      //views with data from normal array 'h' which always lives on host memory
+    }
+  }
+
+  Kokkos::deep_copy(tmp, h_view); //Deepcopy data from h_view(host) to tmp(device)
+
+  d = tmp;
+}
+
+template <class DeviceType>
+std::vector<F_FLOAT> PairUF3Kokkos<DeviceType>::get_constants(double *knots, double coefficient)
+{
+
+  std::vector<F_FLOAT> constants(16);
+
+  constants[0] = 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]));
+  constants[1] = 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]));
+  constants[2] = 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]));
+  constants[3] = 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[4] = 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)));
+  constants[5] = 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)));
+  constants[6] = 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)));
+  constants[7] = 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[8] = 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)));
+  constants[9] = 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)));
+  constants[10] = 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)));
+  constants[11] = 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[12] = 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)));
+  constants[13] = 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)));
+  constants[14] = 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)));
+  constants[15] = 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)));
+
+  return constants;
+}
+
+template <class DeviceType>
+std::vector<F_FLOAT> PairUF3Kokkos<DeviceType>::get_dnconstants(double *knots, double coefficient)
+{
+  std::vector<F_FLOAT> constants(9);
+
+  constants[0] = coefficient *
+      (pow(knots[0], 2) /
+       (pow(knots[0], 2) - knots[0] * knots[1] - knots[0] * knots[2] + knots[1] * knots[2]));
+  constants[1] = coefficient *
+      (-2 * knots[0] /
+       (pow(knots[0], 2) - knots[0] * knots[1] - knots[0] * knots[2] + knots[1] * knots[2]));
+  constants[2] = coefficient *
+      (1 / (pow(knots[0], 2) - knots[0] * knots[1] - knots[0] * knots[2] + knots[1] * knots[2]));
+  constants[3] = 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)));
+  constants[4] = 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)));
+  constants[5] = 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[6] = coefficient *
+      (pow(knots[3], 2) /
+       (knots[1] * knots[2] - knots[1] * knots[3] - knots[2] * knots[3] + pow(knots[3], 2)));
+  constants[7] = coefficient *
+      (-2 * knots[3] /
+       (knots[1] * knots[2] - knots[1] * knots[3] - knots[2] * knots[3] + pow(knots[3], 2)));
+  constants[8] = coefficient *
+      (1 / (knots[1] * knots[2] - knots[1] * knots[3] - knots[2] * knots[3] + pow(knots[3], 2)));
+
+  return constants;
+}
+
+template <class DeviceType>
+double PairUF3Kokkos<DeviceType>::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);
+  int interaction_id = map2b(itype, jtype);
+  int start_index = 3;
+  while (r > d_n2b_knot(interaction_id, start_index + 1)) start_index++;
+
+  if (r < d_cutsq(itype, jtype)) {
+    F_FLOAT r_values[4];
+    r_values[0] = 1;
+    r_values[1] = r;
+    r_values[2] = r_values[1] * r_values[1];
+    r_values[3] = r_values[2] * r_values[1];
+
+    // Calculate energy
+    value = constants_2b(interaction_id, start_index, 0);
+    value += r_values[1] * constants_2b(interaction_id, start_index, 1);
+    value += r_values[2] * constants_2b(interaction_id, start_index, 2);
+    value += r_values[3] * constants_2b(interaction_id, start_index, 3);
+    value += constants_2b(interaction_id, start_index - 1, 4);
+    value += r_values[1] * constants_2b(interaction_id, start_index - 1, 5);
+    value += r_values[2] * constants_2b(interaction_id, start_index - 1, 6);
+    value += r_values[3] * constants_2b(interaction_id, start_index - 1, 7);
+    value += constants_2b(interaction_id, start_index - 2, 8);
+    value += r_values[1] * constants_2b(interaction_id, start_index - 2, 9);
+    value += r_values[2] * constants_2b(interaction_id, start_index - 2, 10);
+    value += r_values[3] * constants_2b(interaction_id, start_index - 2, 11);
+    value += constants_2b(interaction_id, start_index - 3, 12);
+    value += r_values[1] * constants_2b(interaction_id, start_index - 3, 13);
+    value += r_values[2] * constants_2b(interaction_id, start_index - 3, 14);
+    value += r_values[3] * constants_2b(interaction_id, start_index - 3, 15);
+
+    // Calculate force
+    fforce = dnconstants_2b(interaction_id, start_index - 1, 0);
+    fforce += r_values[1] * dnconstants_2b(interaction_id, start_index - 1, 1);
+    fforce += r_values[2] * dnconstants_2b(interaction_id, start_index - 1, 2);
+    fforce += dnconstants_2b(interaction_id, start_index - 2, 3);
+    fforce += r_values[1] * dnconstants_2b(interaction_id, start_index - 2, 4);
+    fforce += r_values[2] * dnconstants_2b(interaction_id, start_index - 2, 5);
+    fforce += dnconstants_2b(interaction_id, start_index - 3, 6);
+    fforce += r_values[1] * dnconstants_2b(interaction_id, start_index - 3, 7);
+    fforce += r_values[2] * dnconstants_2b(interaction_id, start_index - 3, 8);
+  }
+
+  return factor_lj * value;
+}
+
+namespace LAMMPS_NS {
+template class PairUF3Kokkos<LMPDeviceType>;
+#ifdef KOKKOS_ENABLE_CUDA
+template class PairUF3Kokkos<LMPHostType>;
+#endif
+}    // namespace LAMMPS_NS
diff --git a/lammps_plugin/ML-UF3_KOKKOS/pair_uf3_kokkos.h b/lammps_plugin/ML-UF3_KOKKOS/pair_uf3_kokkos.h
new file mode 100644
index 00000000..193c65d6
--- /dev/null
+++ b/lammps_plugin/ML-UF3_KOKKOS/pair_uf3_kokkos.h
@@ -0,0 +1,199 @@
+/* -*- 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/kk,PairUF3Kokkos<LMPDeviceType>)
+PairStyle(uf3/kk/device,PairUF3Kokkos<LMPDeviceType>)
+// clang-format on
+#else
+
+#ifndef LMP_PAIR_UF3_KOKKOS_H
+#define LMP_PAIR_UF3_KOKKOS_H
+
+#include "kokkos.h"
+#include "pair_kokkos.h"
+#include "pair_uf3.h"
+
+template <int NEIGHFLAG, int EVFLAG> struct TagPairUF3ComputeFullA {};
+struct TagPairUF3ComputeShortNeigh {};
+
+namespace LAMMPS_NS {
+
+template <class DeviceType> class PairUF3Kokkos : public PairUF3 {
+ public:
+  PairUF3Kokkos(class LAMMPS *);
+  ~PairUF3Kokkos() override;
+  void compute(int, int) override;
+  void settings(int, char **) override;
+  void coeff(int, char **) override;
+  void allocate();
+  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;
+
+  template <typename T, typename V> void copy_2d(V &d, T **h, int m, int n);
+  template <typename T, typename V> void copy_3d(V &d, T ***h, int m, int n, int o);
+
+  template <int NEIGHFLAG, int EVFLAG>
+  KOKKOS_INLINE_FUNCTION void operator()(TagPairUF3ComputeFullA<NEIGHFLAG, EVFLAG>, const int &,
+                                         EV_FLOAT &) const;
+
+  template <int NEIGHFLAG, int EVFLAG>
+  KOKKOS_INLINE_FUNCTION void operator()(TagPairUF3ComputeFullA<NEIGHFLAG, EVFLAG>,
+                                         const int &) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void operator()(TagPairUF3ComputeShortNeigh, const int &) const;
+
+  enum { EnabledNeighFlags = FULL };
+  enum { COUL_FLAG = 0 };
+  typedef DeviceType device_type;
+  typedef ArrayTypes<DeviceType> AT;
+  typedef EV_FLOAT value_type;
+
+ protected:
+  typename AT::tdual_ffloat_2d k_cutsq;//Create a DualView, defination of tdual_ffloat_2d in kokkos_type.h
+  typename AT::t_ffloat_2d d_cutsq; //t_ffloat_2d = t_dev ==> Creates a new View d_cutsq
+  //the type of d_cutsq is decided by the Device(not host) type for the DualView k_cutsq
+  //Meaning the memory location of d_cutsq is the same as the Device(not host) memory location of
+  //k_cutsq
+  typedef Kokkos::DualView<F_FLOAT***, Kokkos::LayoutRight, DeviceType> tdual_ffloat_3d;
+  typedef Kokkos::DualView<F_FLOAT****, Kokkos::LayoutRight, DeviceType> tdual_ffloat_4d;
+  tdual_ffloat_3d k_cut_3b;
+  tdual_ffloat_4d k_min_cut_3b;
+  typename tdual_ffloat_3d::t_dev d_cut_3b;
+  typename tdual_ffloat_4d::t_dev d_min_cut_3b;
+  template <typename TYPE> void destroy_3d(TYPE data, typename TYPE::value_type*** &array);
+  template <typename TYPE> void destroy_4d(TYPE data, typename TYPE::value_type**** &array);
+  Kokkos::View<F_FLOAT **, LMPDeviceType::array_layout, LMPDeviceType> /*d_cutsq,*/ d_cut_3b_list;
+  //Kokkos::View<F_FLOAT ***, LMPDeviceType::array_layout, LMPDeviceType> d_cut_3b;
+
+  Kokkos::View<F_FLOAT **, LMPDeviceType::array_layout, LMPDeviceType> d_coefficients_2b;
+  Kokkos::View<F_FLOAT **, LMPDeviceType::array_layout, LMPDeviceType> d_dncoefficients_2b;
+  Kokkos::View<F_FLOAT **, LMPDeviceType::array_layout, LMPDeviceType> d_n2b_knot;
+  Kokkos::View<F_FLOAT *, LMPDeviceType::array_layout, LMPDeviceType> d_n2b_knot_spacings;
+  Kokkos::View<int **, LMPDeviceType::array_layout, LMPDeviceType> map2b;
+  Kokkos::View<F_FLOAT[4][4], LMPDeviceType::array_layout, LMPDeviceType> constants;
+  Kokkos::View<F_FLOAT[3][3], LMPDeviceType::array_layout, LMPDeviceType> dnconstants;
+  Kokkos::View<F_FLOAT ***, LMPDeviceType::array_layout, LMPDeviceType> d_n3b_knot_matrix;
+  Kokkos::View<F_FLOAT ****, LMPDeviceType::array_layout, LMPDeviceType> d_coefficients_3b;
+  Kokkos::View<F_FLOAT *****, LMPDeviceType::array_layout, LMPDeviceType> d_dncoefficients_3b;
+  Kokkos::View<F_FLOAT **, LMPDeviceType::array_layout, LMPDeviceType> d_n3b_knot_spacings;
+  Kokkos::View<F_FLOAT **, LMPDeviceType::array_layout, LMPDeviceType> d_n3b_knot_matrix_spacings;
+  Kokkos::View<int ***, LMPDeviceType::array_layout, LMPDeviceType> map3b;
+
+  Kokkos::View<F_FLOAT **[16], LMPDeviceType::array_layout, LMPDeviceType> constants_2b;
+  Kokkos::View<F_FLOAT **[9], LMPDeviceType::array_layout, LMPDeviceType> dnconstants_2b;
+  Kokkos::View<F_FLOAT ***[16], LMPDeviceType::array_layout, LMPDeviceType> constants_3b;
+  Kokkos::View<F_FLOAT ***[9], LMPDeviceType::array_layout, LMPDeviceType> dnconstants_3b;
+
+  std::vector<F_FLOAT> get_constants(double *knots, double coefficient);
+  std::vector<F_FLOAT> get_dnconstants(double *knots, double coefficient);
+
+  int coefficients_created = 0;
+  void create_coefficients();
+  void create_3b_coefficients();
+  void create_2b_coefficients();
+  std::vector<F_FLOAT> get_coefficients(const double *knots, const double coefficient) const;
+  std::vector<F_FLOAT> get_dncoefficients(const double *knots, const double coefficient) const;
+
+  template <int EVFLAG>
+  void twobody(const int itype, const int jtype, const F_FLOAT r, F_FLOAT &evdwl,
+               F_FLOAT &fpair) const;
+  template <int EVFLAG>
+  void threebody(const int itype, const int jtype, const int ktype, const F_FLOAT value_rij,
+                 const F_FLOAT value_rik, const F_FLOAT value_rjk, F_FLOAT &evdwl3,
+                 F_FLOAT (&fforce)[3]) const;
+
+  template <int NEIGHFLAG>
+  KOKKOS_INLINE_FUNCTION void
+  ev_tally(EV_FLOAT &ev, const int &i, const int &j, const F_FLOAT &epair, const F_FLOAT &fpair,
+           const F_FLOAT &delx, const F_FLOAT &dely, const F_FLOAT &delz) const;
+
+  template <int NEIGHFLAG>
+  KOKKOS_INLINE_FUNCTION void ev_tally3(EV_FLOAT &ev, const int &i, const int &j, int &k,
+                                        const F_FLOAT &evdwl, const F_FLOAT &ecoul, F_FLOAT *fj,
+                                        F_FLOAT *fk, F_FLOAT *drji, F_FLOAT *drki) const;
+
+  typename AT::t_x_array_randomread x;
+  typename AT::t_f_array f;
+  typename AT::t_tagint_1d tag;
+  typename AT::t_int_1d_randomread type;
+
+  DAT::tdual_efloat_1d k_eatom;
+  DAT::tdual_virial_array k_vatom;
+  typename AT::t_efloat_1d d_eatom;
+  typename AT::t_virial_array d_vatom;
+
+  using ScatterFType = Kokkos::Experimental::ScatterView<F_FLOAT *[3], Kokkos::LayoutRight,
+                                                         typename DeviceType::memory_space>;
+  ScatterFType fscatter;
+  using ScatterVType = Kokkos::Experimental::ScatterView<F_FLOAT *[6], Kokkos::LayoutRight,
+                                                         typename DeviceType::memory_space>;
+  ScatterVType vscatter;
+  using ScatterCVType = Kokkos::Experimental::ScatterView<F_FLOAT *[9], Kokkos::LayoutRight,
+                                                          typename DeviceType::memory_space>;
+  ScatterCVType cvscatter;
+  using ScatterEType = Kokkos::Experimental::ScatterView<E_FLOAT *, LMPDeviceType::array_layout,
+                                                         typename DeviceType::memory_space>;
+  ScatterEType escatter;
+
+  typename AT::t_neighbors_2d d_neighbors;
+  typename AT::t_int_1d_randomread d_ilist;
+  typename AT::t_int_1d_randomread d_numneigh;
+
+  int neighflag, newton_pair;
+  int nlocal, nall, eflag, vflag;
+
+  int inum;
+  Kokkos::View<int **, DeviceType> d_neighbors_short;
+  Kokkos::View<int *, DeviceType> d_numneigh_short;
+
+  friend void pair_virial_fdotr_compute<PairUF3Kokkos>(PairUF3Kokkos *);
+};
+
+KOKKOS_INLINE_FUNCTION int min(int i, int j)
+{
+  return i < j ? i : j;
+}
+KOKKOS_INLINE_FUNCTION int max(int i, int j)
+{
+  return i > j ? i : j;
+}
+
+}    // 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/README.rst b/lammps_plugin/README.rst
index 26fb1d69..9b7e88ed 100644
--- a/lammps_plugin/README.rst
+++ b/lammps_plugin/README.rst
@@ -38,9 +38,21 @@ If you want to compile with MPI support enabled, you additionally have to set th
 
 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` found here_.
+The source files for kokkos accelerator variants of uf3 (i.e. :code:`pair_style uf3/kk` with GPU support) are available in the :code:`ML-UF3_KOKKOS` directory. To compile lammps with support for the kokkos accelerator variants of uf3 do-
 
-.. _here: https://github.com/monk-04/uf3/tree/lammps_implementation_v2/lammps_plugin#kokkos
+.. code:: bash
+
+ cp ML-UF3_KOKKOS/pair_uf3_kokkos.cpp ML-UF3_KOKKOS/pair_uf3_kokkos.h LAMMPS_BASE_DIR/src/KOKKOS/.
+ mkdir LAMMPS_BASE_DIR/build_kokkos
+ cmake -D LAMMPS_EXCEPTIONS=ON -D BUILD_SHARED_LIBS=ON -D PKG_KOKKOS=yes KOKKOS_FLAGS -D PKG_ML-UF3=ON -D CMAKE_POSITION_INDEPENDENT_CODE=ON -D CMAKE_EXE_FLAGS="-dynamic" ../cmake
+
+where :code:`KOKKOS_FLAGS` are machine specific flags. For example when compiling for Nvidia A100 GPUs-
+
+.. code:: bash
+
+ cmake -D LAMMPS_EXCEPTIONS=ON -D BUILD_SHARED_LIBS=ON -D PKG_KOKKOS=yes -D Kokkos_ARCH_AMPERE80=ON -D Kokkos_ENABLE_CUDA=yes -D PKG_ML-UF3=ON -D CMAKE_POSITION_INDEPENDENT_CODE=ON -D CMAKE_EXE_FLAGS="-dynamic" ../cmake
+
+Refer to `lammps documentation <https://docs.lammps.org/Speed_kokkos.html>`_ for more details.
 
 =====
 Running lammps with UF3 potential
diff --git a/lammps_plugin/documentation_lammps_format/pair_uf3.rst b/lammps_plugin/documentation_lammps_format/pair_uf3.rst
index b2413382..f5821fb6 100644
--- a/lammps_plugin/documentation_lammps_format/pair_uf3.rst
+++ b/lammps_plugin/documentation_lammps_format/pair_uf3.rst
@@ -37,10 +37,15 @@ Examples
     pair_coeff 2 2 Sn_Sn
 
     pair_style uf3 3 2
-    pair_coeff 1 * Nb_Nb
-    pair_coeff 2 * Sn_Sn
-    pair_style 3b 1 * * Nb_Sn_Sn
-    pair_style 3b 2 * * Sn_Sn_Sn
+    pair_coeff 1 1 Nb_Nb
+    pair_coeff 1 2 Nb_Sn
+    pair_coeff 2 2 Sn_Sn
+    pair_style 3b 1 1 1 Nb_Nb_Nb
+    pair_style 3b 1 1 2 Nb_Nb_Sn
+    pair_style 3b 1 2 2 Nb_Sn_Sn
+    pair_style 3b 2 1 1 Sn_Nb_Nb
+    pair_style 3b 2 1 2 Sn_Nb_Sn
+    pair_style 3b 2 2 2 Sn_Sn_Sn
 
 Description
 """""""""""
@@ -51,13 +56,13 @@ The *uf3* style computes the :ref:`Ultra-Fast Force Fields (UF3) <Xie23>` potent
 
     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_2(r_{ij}) & = \sum_{n=0}^N 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})
+    V_3 (r_{ij},r_{ik},r_{jk}) & = \sum_{l=0}^N_l \sum_{m=0}^N_m \sum_{n=0}^N_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.
+where :math:`V_2(r_{ij})` and :math:`V_3 (r_{ij},r_{ik},r_{jk})` are the two- and three-body interactions, respectively. 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:`N`, :math:`N_l`, :math:`N_m`, and :math:`N_n` denote the number of basis functions per spline or tensor spline dimension.
 
-The UF3 LAMMPS potential files are provided using multiple pair_coeff commands. 
+The UF3 LAMMPS potential files are provided using multiple pair_coeff commands. A single UF3 LAMMPS potential file contains information about one particular interaction only.
 
 .. note::
 
@@ -78,6 +83,8 @@ As an example, if a LAMMPS simulation contains 2 atom types (elements 'A' and 'B
     pair_coeff 3b 2 1 2 B_A_B
     pair_coeff 3b 2 2 2 B_B_B
 
+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/uf3/uf3/tree/master/lammps_plugin/scripts>`_) for generating the UF3 LAMMPS potential files from the UF3 JSON potentials.
+
 LAMMPS wild-card character "*" can also be used to specify a single UF3 LAMMPS potential file for multiple interaction. For example- 
 
 .. code-block:: LAMMPS
@@ -87,9 +94,10 @@ LAMMPS wild-card character "*" can also be used to specify a single UF3 LAMMPS p
     pair_coeff 3b 1 * * A_A_A
     pair_coeff 3b 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/uf3/uf3/tree/master/lammps_plugin/scripts>`_) for generating the UF3 LAMMPS potential files from the UF3 JSON potentials.
+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. Note, using a single interaction file for all types of interactions is **not** the recommended way of using :code:`pair_style uf3` and will often lead to **incorrect results**.
+
 
-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-
+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 end with :code:`#` characters. Following shows the format of a generic 2-body UF3 LAMMPS potential file-
 
 .. code-block:: LAMMPS
 
@@ -174,4 +182,4 @@ none
 
 .. _Xie23:
 
-**(Xie23)** S. R. Xie, M. Rupp, and R. G. Hennig, "Ultra-fast interpretable machine-learning potentials", preprint arXiv:2110.00624v2 (2023)
+**(Xie23)** Xie, S.R., Rupp, M. & Hennig, R.G. Ultra-fast interpretable machine-learning potentials. npj Comput Mater 9, 162 (2023). https://doi.org/10.1038/s41524-023-01092-7
diff --git a/lammps_plugin/scripts/README.rst b/lammps_plugin/scripts/README.rst
index 5c0a0d50..ca57387d 100644
--- a/lammps_plugin/scripts/README.rst
+++ b/lammps_plugin/scripts/README.rst
@@ -1,4 +1,4 @@
-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).
+This document describes how to generate UF3 lammps potential from an UF3 model file(JSON).
 
 Usuage-
 
@@ -6,9 +6,7 @@ Usuage-
 
    python generate_uf3_lammps_pots.py 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. 
+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`. 
 
 
 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/timing_info/README.rst b/lammps_plugin/timing_info/README.rst
index e8d1cd77..3e1c7935 100644
--- a/lammps_plugin/timing_info/README.rst
+++ b/lammps_plugin/timing_info/README.rst
@@ -1,3 +1,5 @@
+*This documented is outdated and should not be referred to*
+
 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.
diff --git a/lammps_plugin/tungsten_example/W_nve.in b/lammps_plugin/tungsten_example/W_nve.in
index 32071a00..5689e297 100644
--- a/lammps_plugin/tungsten_example/W_nve.in
+++ b/lammps_plugin/tungsten_example/W_nve.in
@@ -15,7 +15,7 @@ neigh_modify    every 5 delay 0
 pair_style      zero 10
 pair_style      uf3 3 1
 pair_coeff      1 1 W_W
-pair_coeff      1 1 W_W_W
+pair_coeff      3b 1 1 W_W_W
 
 # ============= give all atoms velocity
 timestep        0.001
diff --git a/lammps_plugin/tungsten_example/relax_lammps_3b.in b/lammps_plugin/tungsten_example/relax_lammps_3b.in
index 1be90759..f72d0e42 100644
--- a/lammps_plugin/tungsten_example/relax_lammps_3b.in
+++ b/lammps_plugin/tungsten_example/relax_lammps_3b.in
@@ -11,7 +11,7 @@ neighbor	0.0 bin
 pair_style  zero 10.0
 pair_style  uf3 3 1
 pair_coeff  1 1 W_W
-pair_coeff  1 1 1 W_W_W
+pair_coeff  3b 1 1 1 W_W_W
 
 #----------------------Run Minimization-------------------------
 reset_timestep	0


Test patch:
diff --git a/lammps_plugin/unittest/A_A.uf3_pot b/lammps_plugin/unittest/A_A.uf3_pot
index cb9371d8..a50b48ba 100644
--- a/lammps_plugin/unittest/A_A.uf3_pot
+++ b/lammps_plugin/unittest/A_A.uf3_pot
@@ -1,7 +1,7 @@
 #UF3 POT
-2B 0 3
+2B 0 3 nk
 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
index 1f089fda..18f79c3b 100644
--- a/lammps_plugin/unittest/A_A_A.uf3_pot
+++ b/lammps_plugin/unittest/A_A_A.uf3_pot
@@ -1,5 +1,5 @@
 #UF3 POT
-3B 0 3
+3B 0 3 nk
 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
@@ -126,4 +126,4 @@
 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
 0.0 0.0 0.0 0.0 0.0 0.0 0.0 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/manybody-pair-uf3.yaml b/lammps_plugin/unittest/manybody-pair-uf3.yaml
index 259f66e6..8ba4a24e 100644
--- a/lammps_plugin/unittest/manybody-pair-uf3.yaml
+++ b/lammps_plugin/unittest/manybody-pair-uf3.yaml
@@ -1,9 +1,9 @@
 ---
 lammps_version: 24 Mar 2022
 tags: generated
-date_generated: Thu Jun  8 23:16:30 2023
-epsilon: 1e-9
-skip_tests: ""
+date_generated: Tue Nov 28 14:51:03 2023
+epsilon: 1e-09
+skip_tests: 
 prerequisites: ! |
   pair uf3
 pre_commands: ! |
@@ -15,7 +15,7 @@ input_file: in.manybody
 pair_style: uf3 3 8
 pair_coeff: ! |
   * * A_A.uf3_pot
-  * * * A_A_A.uf3_pot
+  3b * * * A_A_A.uf3_pot
 extract: ! ""
 natoms: 64
 init_vdwl: -76.14388662099438