Fixing few bugs with histories, removing indices from granular model

src/GRANULAR/fix_wall_gran.cpp

@@ -492,7 +492,6 @@ void FixWallGran::post_force(int /*vflag*/)
     }
 
     // Reset model and copy initial geometric data
-    model->i = i;
     model->dx[0] = dx;
     model->dx[1] = dy;
     model->dx[2] = dz;

src/GRANULAR/fix_wall_gran_region.cpp

@@ -229,7 +229,6 @@ void FixWallGranRegion::post_force(int /*vflag*/)
     for (int ic = 0; ic < nc; ic++) {
 
       // Reset model and copy initial geometric data
-      model->i = i;
       model->dx[0] = region->contact[ic].delx;
       model->dx[1] = region->contact[ic].dely;
       model->dx[2] = region->contact[ic].delz;

src/GRANULAR/granular_model.h

@@ -87,7 +87,6 @@ class GranularModel : protected Pointers {
   double radi, radj, meff, dt, Ti, Tj, area;
   double Fntot, magtortwist;
 
-  int i, j;
   double *xi, *xj, *vi, *vj, *omegai, *omegaj;
   double fs[3], fr[3], ft[3];
 

src/GRANULAR/gsm_tangential.cpp

@@ -44,7 +44,7 @@ GSMTangentialNone::GSMTangentialNone(GranularModel *gm, LAMMPS *lmp) : GSMTangen
 GSMTangentialLinearNoHistory::GSMTangentialLinearNoHistory(GranularModel *gm, LAMMPS *lmp) : GSMTangential(gm, lmp)
 {
   num_coeffs = 2;
-  size_history = 3;
+  size_history = 0;
 }
 
 /* ---------------------------------------------------------------------- */

src/GRANULAR/pair_granular.cpp

@@ -122,6 +122,9 @@ void PairGranular::compute(int eflag, int vflag)
   bool touchflag = false;
   const bool history_update = update->setupflag == 0;
 
+  for (int i = 0; i < vec_models.size(); i++)
+    vec_models[i].history_update = history_update;
+
   ev_init(eflag,vflag);
 
   // update rigid body info for owned & ghost atoms if using FixRigid masses
@@ -190,13 +193,12 @@ void PairGranular::compute(int eflag, int vflag)
       jtype = type[j];
 
       // Reset model and copy initial geometric data
-      models[itype][jtype]->i = i;
-      models[itype][jtype]->j = j;
       models[itype][jtype]->xi = x[i];
       models[itype][jtype]->xj = x[j];
       models[itype][jtype]->radi = radius[i];
       models[itype][jtype]->radj = radius[j];
-      models[itype][jtype]->touch = touch[jj];
+      if (use_history || models[itype][jtype]->beyond_contact)
+        models[itype][jtype]->touch = touch[jj];
 
       touchflag = models[itype][jtype]->check_contact();
 
@@ -234,12 +236,11 @@ void PairGranular::compute(int eflag, int vflag)
       models[itype][jtype]->vj = v[j];
       models[itype][jtype]->omegai = omega[i];
       models[itype][jtype]->omegaj = omega[j];
-      models[itype][jtype]->history_update = history_update;
       if (use_history) {
         history = &allhistory[size_history*jj];
         models[itype][jtype]->history = history;
       }
-;
+
       if (heat_flag) {
         models[itype][jtype]->Ti = temperature[i];
         models[itype][jtype]->Tj = temperature[j];
@@ -411,7 +412,10 @@ void PairGranular::init_style()
 
   for (auto &model : vec_models) {
     model.init();
-    if (model.beyond_contact) beyond_contact = 1;
+    if (model.beyond_contact) {
+      beyond_contact = 1;
+      use_history = 1; // Need to track if in contact
+    }
     if (model.size_history != 0) use_history = 1;
 
     for (i = 0; i < NSUBMODELS; i++)
@@ -700,14 +704,28 @@ double PairGranular::single(int i, int j, int itype, int jtype,
   double *radius = atom->radius;
 
   // Reset model and copy initial geometric data
-  models[itype][jtype]->i = i;
-  models[itype][jtype]->j = j;
   models[itype][jtype]->xi = x[i];
   models[itype][jtype]->xj = x[j];
   models[itype][jtype]->radi = radius[i];
   models[itype][jtype]->radj = radius[j];
   models[i][j]->history_update = 0; // Don't update history
 
+  // If history is needed
+  jnum = list->numneigh[i];
+  jlist = list->firstneigh[i];
+  if (use_history) {
+    if ((fix_history == nullptr) || (fix_history->firstvalue == nullptr))
+      error->one(FLERR,"Pair granular single computation needs history");
+    allhistory = fix_history->firstvalue[i];
+    for (int jj = 0; jj < jnum; jj++) {
+      neighprev++;
+      if (neighprev >= jnum) neighprev = 0;
+      if (jlist[neighprev] == j) break;
+    }
+    history = &allhistory[size_history*neighprev];
+    models[itype][jtype]->touch = fix_history->firstflag[i][neighprev];
+  }
+
   touchflag = models[itype][jtype]->check_contact();
 
   if (!touchflag) {
@@ -733,22 +751,6 @@ double PairGranular::single(int i, int j, int itype, int jtype,
   if (mask[i] & freeze_group_bit) meff = mj;
   if (mask[j] & freeze_group_bit) meff = mi;
 
-  // if any history is needed
-  jnum = list->numneigh[i];
-  jlist = list->firstneigh[i];
-
-  if (use_history) {
-    if ((fix_history == nullptr) || (fix_history->firstvalue == nullptr))
-      error->one(FLERR,"Pair granular single computation needs history");
-    allhistory = fix_history->firstvalue[i];
-    for (int jj = 0; jj < jnum; jj++) {
-      neighprev++;
-      if (neighprev >= jnum) neighprev = 0;
-      if (jlist[neighprev] == j) break;
-    }
-    history = &allhistory[size_history*neighprev];
-  }
-
   // Copy additional information and calculate forces
   double **v = atom->v;
   double **omega = atom->omega;