reaxc: Fixes to allow multiple object instances

Move LR into system structure
Make fix_reaxc have unique id for each instance

closes #1432

src/KOKKOS/pair_reaxc_kokkos.cpp

@@ -143,7 +143,7 @@ template<class DeviceType>
 void PairReaxCKokkos<DeviceType>::init_style()
 {
   PairReaxC::init_style();
-  if (fix_reax) modify->delete_fix("REAXC"); // not needed in the Kokkos version
+  if (fix_reax) modify->delete_fix(fix_id.c_str()); // not needed in the Kokkos version
   fix_reax = NULL;
 
   // irequest = neigh request made by parent class
@@ -340,6 +340,7 @@ void PairReaxCKokkos<DeviceType>::init_md()
 {
   // init_taper()
   F_FLOAT d1, d7, swa, swa2, swa3, swb, swb2, swb3;
+  LR_lookup_table ** & LR = system->LR;
 
   swa = control->nonb_low;
   swb = control->nonb_cut;
@@ -437,6 +438,7 @@ int PairReaxCKokkos<DeviceType>::Init_Lookup_Tables()
   double dr;
   double *h, *fh, *fvdw, *fele, *fCEvd, *fCEclmb;
   double v0_vdw, v0_ele, vlast_vdw, vlast_ele;
+  LR_lookup_table ** & LR = system->LR;
 
   /* initializations */
   v0_vdw = 0;
@@ -541,6 +543,7 @@ void PairReaxCKokkos<DeviceType>::Deallocate_Lookup_Tables()
 {
   int i, j;
   int ntypes;
+  LR_lookup_table ** & LR = system->LR;
 
   ntypes = atom->ntypes;
 

src/USER-OMP/pair_reaxc_omp.cpp

@@ -33,6 +33,8 @@
    High Performance Computing Applications, to appear.
  ------------------------------------------------------------------------- */
 
+#include <sstream>
+
 #include "pair_reaxc_omp.h"
 #include "atom.h"
 #include "update.h"
@@ -97,6 +99,10 @@ PairReaxCOMP::PairReaxCOMP(LAMMPS *lmp) : PairReaxC(lmp), ThrOMP(lmp, THR_PAIR)
 
   num_nbrs_offset = NULL;
 
+  std::stringstream ss;
+  ss << "REAXC_COMP_" << std::dec << Pair::instance_me;
+  fix_id = ss.str();
+
 #ifdef OMP_TIMING
   for (int i=0;i<LASTTIMINGINDEX;i++) {
     ompTimingData[i] = 0;
@@ -373,7 +379,7 @@ void PairReaxCOMP::init_style( )
 
   if (fix_reax == NULL) {
     char **fixarg = new char*[3];
-    fixarg[0] = (char *) "REAXC";
+    fixarg[0] = (char *) fix_id.c_str();
     fixarg[1] = (char *) "all";
     fixarg[2] = (char *) "REAXC";
     modify->add_fix(3,fixarg);
@@ -643,4 +649,3 @@ void PairReaxCOMP::FindBond()
     }
   }
 }
-

src/USER-OMP/pair_reaxc_omp.h

@@ -76,6 +76,7 @@ class PairReaxCOMP : public PairReaxC, public ThrOMP {
   }
 
  protected:
+  std::string fix_id;
   virtual void setup();
   virtual void write_reax_atoms();
   virtual int estimate_reax_lists();

src/USER-OMP/reaxc_nonbonded_omp.cpp

@@ -281,6 +281,7 @@ void Tabulated_vdW_Coulomb_Energy_OMP(reax_system *system,control_params *contro
   int tid = 0;
 #endif
   long froffset = (system->N * tid);
+  LR_lookup_table ** & LR = system->LR;
 
   class PairReaxCOMP *pair_reax_ptr;
   pair_reax_ptr = static_cast<class PairReaxCOMP*>(system->pair_ptr);

src/USER-REAXC/pair_reaxc.cpp

@@ -20,6 +20,8 @@
    Hybrid and hybrid/overlay compatibility added by Ray Shan (Sandia)
 ------------------------------------------------------------------------- */
 
+#include <sstream>
+
 #include "pair_reaxc.h"
 #include "atom.h"
 #include "update.h"
@@ -74,6 +76,10 @@ PairReaxC::PairReaxC(LAMMPS *lmp) : Pair(lmp)
   manybody_flag = 1;
   ghostneigh = 1;
 
+  std::stringstream ss;
+  ss << "REAXC_" << std::dec << Pair::instance_me;
+  fix_id = ss.str();
+
   system = (reax_system *)
     memory->smalloc(sizeof(reax_system),"reax:system");
   memset(system,0,sizeof(reax_system));
@@ -135,7 +141,7 @@ PairReaxC::~PairReaxC()
 {
   if (copymode) return;
 
-  if (fix_reax) modify->delete_fix("REAXC");
+  if (fix_reax) modify->delete_fix(fix_id.c_str());
 
   if (setup_flag) {
     Close_Output_Files( system, control, out_control, mpi_data );
@@ -412,7 +418,7 @@ void PairReaxC::init_style( )
 
   if (fix_reax == NULL) {
     char **fixarg = new char*[3];
-    fixarg[0] = (char *) "REAXC";
+    fixarg[0] = (char *) fix_id.c_str();
     fixarg[1] = (char *) "all";
     fixarg[2] = (char *) "REAXC";
     modify->add_fix(3,fixarg);

src/USER-REAXC/pair_reaxc.h

@@ -60,6 +60,7 @@ class PairReaxC : public Pair {
   bigint ngroup;
 
  protected:
+  std::string fix_id;
   double cutmax;
   int nelements;                // # of unique elements
   char **elements;              // names of unique elements

src/USER-REAXC/reaxc_lookup.cpp

@@ -29,7 +29,6 @@
 #include "reaxc_nonbonded.h"
 #include "reaxc_tool_box.h"
 
-LR_lookup_table **LR;
 
 void Tridiagonal_Solve( const double *a, const double *b,
                         double *c, double *d, double *x, unsigned int n){
@@ -157,6 +156,7 @@ int Init_Lookup_Tables( reax_system *system, control_params *control,
   double dr;
   double *h, *fh, *fvdw, *fele, *fCEvd, *fCEclmb;
   double v0_vdw, v0_ele, vlast_vdw, vlast_ele;
+  LR_lookup_table ** & LR = system->LR;
 
   /* initializations */
   v0_vdw = 0;
@@ -273,6 +273,7 @@ void Deallocate_Lookup_Tables( reax_system *system )
 {
   int i, j;
   int ntypes;
+  LR_lookup_table ** & LR = system->LR;
 
   ntypes = system->reax_param.num_atom_types;
 

src/USER-REAXC/reaxc_nonbonded.cpp

@@ -220,6 +220,7 @@ void Tabulated_vdW_Coulomb_Energy( reax_system *system,control_params *control,
   far_neighbor_data *nbr_pj;
   reax_list *far_nbrs;
   LR_lookup_table *t;
+  LR_lookup_table ** & LR = system->LR;
 
   natoms = system->n;
   far_nbrs = (*lists) + FAR_NBRS;

src/USER-REAXC/reaxc_types.h

@@ -395,6 +395,8 @@ typedef struct
   double ghost_cutoff;
 } boundary_cutoff;
 
+
+struct _LR_lookup_table;  // forward declaration
 struct _reax_system
 {
   reax_interaction reax_param;
@@ -418,6 +420,8 @@ struct _reax_system
   int mincap;
   double safezone, saferzone;
 
+  _LR_lookup_table **LR;
+
   int omp_active;
 };
 typedef _reax_system reax_system;
@@ -883,7 +887,7 @@ struct cubic_spline_coef
 
 
 
-typedef struct
+typedef struct _LR_lookup_table
 {
   double xmin, xmax;
   int n;
@@ -897,7 +901,6 @@ typedef struct
   cubic_spline_coef *vdW, *CEvd;
   cubic_spline_coef *ele, *CEclmb;
 } LR_lookup_table;
-extern LR_lookup_table **LR;
 
 /* function pointer defs */
 typedef void (*evolve_function)(reax_system*, control_params*,