refactoring of neighbor history

src/GRANULAR/pair_gran_hertz_history.cpp

@@ -24,6 +24,7 @@
 #include "update.h"
 #include "force.h"
 #include "fix.h"
+#include "fix_neigh_history.h"
 #include "neighbor.h"
 #include "neigh_list.h"
 #include "comm.h"
@@ -95,8 +96,8 @@ void PairGranHertzHistory::compute(int eflag, int vflag)
   ilist = list->ilist;
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
-  firsttouch = list->listhistory->firstneigh;
+  firsttouch = fix_history->firstflag;
-  firstshear = list->listhistory->firstdouble;
+  firstshear = fix_history->firstvalue;
 
   // loop over neighbors of my atoms
 
@@ -407,7 +408,7 @@ double PairGranHertzHistory::single(int i, int j, int itype, int jtype,
 
   int jnum = list->numneigh[i];
   int *jlist = list->firstneigh[i];
-  double *allshear = list->listhistory->firstdouble[i];
+  double *allshear = fix_history->firstvalue[i];
 
   for (int jj = 0; jj < jnum; jj++) {
     neighprev++;

src/GRANULAR/pair_gran_hooke_history.cpp

@@ -27,7 +27,7 @@
 #include "update.h"
 #include "modify.h"
 #include "fix.h"
-#include "fix_shear_history.h"
+#include "fix_neigh_history.h"
 #include "comm.h"
 #include "neighbor.h"
 #include "neigh_list.h"
@@ -64,7 +64,7 @@ PairGranHookeHistory::PairGranHookeHistory(LAMMPS *lmp) : Pair(lmp)
 PairGranHookeHistory::~PairGranHookeHistory()
 {
   delete [] svector;
-  if (fix_history) modify->delete_fix("SHEAR_HISTORY");
+  if (fix_history) modify->delete_fix("NEIGH_HISTORY");
 
   if (allocated) {
     memory->destroy(setflag);
@@ -137,8 +137,8 @@ void PairGranHookeHistory::compute(int eflag, int vflag)
   ilist = list->ilist;
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
-  firsttouch = listhistory->firstneigh;
+  firsttouch = fix_history->firstflag;
-  firstshear = listhistory->firstdouble;
+  firstshear = fix_history->firstvalue;
 
   // loop over neighbors of my atoms
 
@@ -400,35 +400,28 @@ void PairGranHookeHistory::init_style()
   if (comm->ghost_velocity == 0)
     error->all(FLERR,"Pair granular requires ghost atoms store velocity");
 
-  // need a granular neigh list and optionally a granular history neigh list
+  // need a granular neigh list
 
   int irequest = neighbor->request(this,instance_me);
   neighbor->requests[irequest]->size = 1;
-  if (history) {
+  if (history) neighbor->requests[irequest]->history = 1;
-    irequest = neighbor->request(this,instance_me);
-    neighbor->requests[irequest]->id = 1;
-    neighbor->requests[irequest]->history = 1;
-    neighbor->requests[irequest]->dnum = 3;
-  }
 
   dt = update->dt;
 
-  // if shear history is stored:
   // if first init, create Fix needed for storing shear history
 
   if (history && fix_history == NULL) {
     char dnumstr[16];
     sprintf(dnumstr,"%d",3);
     char **fixarg = new char*[4];
-    fixarg[0] = (char *) "SHEAR_HISTORY";
+    fixarg[0] = (char *) "NEIGH_HISTORY";
     fixarg[1] = (char *) "all";
-    fixarg[2] = (char *) "SHEAR_HISTORY";
+    fixarg[2] = (char *) "NEIGH_HISTORY";
     fixarg[3] = dnumstr;
-    modify->add_fix(4,fixarg);
+    modify->add_fix(4,fixarg,1);
     delete [] fixarg;
-    fix_history = (FixShearHistory *) modify->fix[modify->nfix-1];
+    fix_history = (FixNeighHistory *) modify->fix[modify->nfix-1];
     fix_history->pair = this;
-    neighbor->requests[irequest]->fix_history = fix_history;
   }
 
   // check for FixFreeze and set freeze_group_bit
@@ -494,23 +487,12 @@ void PairGranHookeHistory::init_style()
   // set fix which stores history info
 
   if (history) {
-    int ifix = modify->find_fix("SHEAR_HISTORY");
+    int ifix = modify->find_fix("NEIGH_HISTORY");
-    if (ifix < 0) error->all(FLERR,"Could not find pair fix ID");
+    if (ifix < 0) error->all(FLERR,"Could not find pair fix neigh history ID");
-    fix_history = (FixShearHistory *) modify->fix[ifix];
+    fix_history = (FixNeighHistory *) modify->fix[ifix];
   }
 }
 
-/* ----------------------------------------------------------------------
-   neighbor callback to inform pair style of neighbor list to use
-   optional granular history list
-------------------------------------------------------------------------- */
-
-void PairGranHookeHistory::init_list(int id, NeighList *ptr)
-{
-  if (id == 0) list = ptr;
-  else if (id == 1) listhistory = ptr;
-}
-
 /* ----------------------------------------------------------------------
    init for one type pair i,j and corresponding j,i
 ------------------------------------------------------------------------- */
@@ -704,7 +686,7 @@ double PairGranHookeHistory::single(int i, int j, int itype, int jtype,
 
   int jnum = list->numneigh[i];
   int *jlist = list->firstneigh[i];
-  double *allshear = list->listhistory->firstdouble[i];
+  double *allshear = fix_history->firstvalue[i];
 
   for (int jj = 0; jj < jnum; jj++) {
     neighprev++;
@@ -797,14 +779,3 @@ double PairGranHookeHistory::memory_usage()
   double bytes = nmax * sizeof(double);
   return bytes;
 }
-
-/* ----------------------------------------------------------------------
-   return ptr to FixShearHistory class
-   called by Neighbor when setting up neighbor lists
-------------------------------------------------------------------------- */
-
-void *PairGranHookeHistory::extract(const char *str, int &dim)
-{
-  if (strcmp(str,"history") == 0) return (void *) fix_history;
-  return NULL;
-}

src/GRANULAR/pair_gran_hooke_history.h

@@ -32,7 +32,6 @@ class PairGranHookeHistory : public Pair {
   virtual void settings(int, char **);
   void coeff(int, char **);
   void init_style();
-  void init_list(int, class NeighList *);
   double init_one(int, int);
   void write_restart(FILE *);
   void read_restart(FILE *);
@@ -43,7 +42,6 @@ class PairGranHookeHistory : public Pair {
   int pack_forward_comm(int, int *, double *, int, int *);
   void unpack_forward_comm(int, int, double *);
   double memory_usage();
-  void *extract(const char *, int &);
 
  protected:
   double kn,kt,gamman,gammat,xmu;
@@ -56,7 +54,7 @@ class PairGranHookeHistory : public Pair {
   double *onerad_dynamic,*onerad_frozen;
   double *maxrad_dynamic,*maxrad_frozen;
 
-  class FixShearHistory *fix_history;
+  class FixNeighHistory *fix_history;
 
   // storage of rigid body masses for use in granular interactions
 

src/fix.h

@@ -113,6 +113,7 @@ class Fix : protected Pointers {
   virtual void setup(int) {}
   virtual void setup_pre_exchange() {}
   virtual void setup_pre_neighbor() {}
+  virtual void setup_post_neighbor() {}
   virtual void setup_pre_force(int) {}
   virtual void setup_pre_reverse(int, int) {}
   virtual void min_setup(int) {}
@@ -120,6 +121,7 @@ class Fix : protected Pointers {
   virtual void post_integrate() {}
   virtual void pre_exchange() {}
   virtual void pre_neighbor() {}
+  virtual void post_neighbor() {}
   virtual void pre_force(int) {}
   virtual void pre_reverse(int,int) {}
   virtual void post_force(int) {}
@@ -155,6 +157,7 @@ class Fix : protected Pointers {
 
   virtual void min_pre_exchange() {}
   virtual void min_pre_neighbor() {}
+  virtual void min_post_neighbor() {}
   virtual void min_pre_force(int) {}
   virtual void min_pre_reverse(int,int) {}
   virtual void min_post_force(int) {}
@@ -244,25 +247,27 @@ namespace FixConst {
   static const int POST_INTEGRATE =          1<<1;
   static const int PRE_EXCHANGE =            1<<2;
   static const int PRE_NEIGHBOR =            1<<3;
-  static const int PRE_FORCE =               1<<4;
+  static const int POST_NEIGHBOR =           1<<4;
-  static const int PRE_REVERSE =             1<<5;
+  static const int PRE_FORCE =               1<<5;
-  static const int POST_FORCE =              1<<6;
+  static const int PRE_REVERSE =             1<<6;
-  static const int FINAL_INTEGRATE =         1<<7;
+  static const int POST_FORCE =              1<<7;
-  static const int END_OF_STEP =             1<<8;
+  static const int FINAL_INTEGRATE =         1<<8;
-  static const int POST_RUN =                1<<9;
+  static const int END_OF_STEP =             1<<9;
-  static const int THERMO_ENERGY =           1<<10;
+  static const int POST_RUN =                1<<10;
-  static const int INITIAL_INTEGRATE_RESPA = 1<<11;
+  static const int THERMO_ENERGY =           1<<11;
-  static const int POST_INTEGRATE_RESPA =    1<<12;
+  static const int INITIAL_INTEGRATE_RESPA = 1<<12;
-  static const int PRE_FORCE_RESPA =         1<<13;
+  static const int POST_INTEGRATE_RESPA =    1<<13;
-  static const int POST_FORCE_RESPA =        1<<14;
+  static const int PRE_FORCE_RESPA =         1<<14;
-  static const int FINAL_INTEGRATE_RESPA =   1<<15;
+  static const int POST_FORCE_RESPA =        1<<15;
-  static const int MIN_PRE_EXCHANGE =        1<<16;
+  static const int FINAL_INTEGRATE_RESPA =   1<<16;
-  static const int MIN_PRE_NEIGHBOR =        1<<17;
+  static const int MIN_PRE_EXCHANGE =        1<<17;
-  static const int MIN_PRE_FORCE =           1<<18;
+  static const int MIN_PRE_NEIGHBOR =        1<<18;
-  static const int MIN_PRE_REVERSE =         1<<19;
+  static const int MIN_POST_NEIGHBOR =       1<<19;
-  static const int MIN_POST_FORCE =          1<<20;
+  static const int MIN_PRE_FORCE =           1<<20;
-  static const int MIN_ENERGY =              1<<21;
+  static const int MIN_PRE_REVERSE =         1<<21;
-  static const int FIX_CONST_LAST =          1<<22;
+  static const int MIN_POST_FORCE =          1<<22;
+  static const int MIN_ENERGY =              1<<23;
+  static const int FIX_CONST_LAST =          1<<24;
 }
 
 }

src/fix_neigh_history.cpp

@@ -14,7 +14,7 @@
 #include <mpi.h>
 #include <string.h>
 #include <stdio.h>
-#include "fix_shear_history.h"
+#include "fix_neigh_history.h"
 #include "atom.h"
 #include "comm.h"
 #include "neighbor.h"
@@ -33,12 +33,12 @@ enum{DEFAULT,NPARTNER,PERPARTNER};
 
 /* ---------------------------------------------------------------------- */
 
-FixShearHistory::FixShearHistory(LAMMPS *lmp, int narg, char **arg) :
+FixNeighHistory::FixNeighHistory(LAMMPS *lmp, int narg, char **arg) :
   Fix(lmp, narg, arg),
-  npartner(NULL), partner(NULL), shearpartner(NULL), pair(NULL), 
+  npartner(NULL), partner(NULL), valuepartner(NULL), pair(NULL), 
-  ipage(NULL), dpage(NULL)
+  ipage_atom(NULL), dpage_atom(NULL), ipage_neigh(NULL), dpage_neigh(NULL)
 {
-  if (narg != 4) error->all(FLERR,"Illegal fix SHEAR_HISTORY command");
+  if (narg != 4) error->all(FLERR,"Illegal fix NEIGH_HISTORY command");
 
   restart_peratom = 1;
   create_attribute = 1;
@@ -48,9 +48,12 @@ FixShearHistory::FixShearHistory(LAMMPS *lmp, int narg, char **arg) :
   dnum = force->inumeric(FLERR,arg[3]);
   dnumbytes = dnum * sizeof(double);
 
+  zeroes = new double[dnum];
+  for (int i = 0; i < dnum; i++) zeroes[i] = 0.0;
+
   onesided = 0;
-  if (strcmp(id,"LINE_SHEAR_HISTORY") == 0) onesided = 1;
+  if (strcmp(id,"LINE_NEIGH_HISTORY") == 0) onesided = 1;
-  if (strcmp(id,"TRI_SHEAR_HISTORY") == 0) onesided = 1;
+  if (strcmp(id,"TRI_NEIGH_HISTORY") == 0) onesided = 1;
 
   if (newton_pair) comm_reverse = 1;   // just for single npartner value
                                        // variable-size history communicated via
@@ -65,11 +68,24 @@ FixShearHistory::FixShearHistory(LAMMPS *lmp, int narg, char **arg) :
 
   pgsize = oneatom = 0;
 
+  // other per-atom vectors
+
+  firstflag = NULL;
+  firstvalue = NULL;
+  maxatom = 0;
+
+  // per-atom and per-neighbor data structs
+
+  ipage_atom = NULL;
+  dpage_atom = NULL;
+  ipage_neigh = NULL;
+  dpage_neigh = NULL;
+
   // initialize npartner to 0 so neighbor list creation is OK the 1st time
 
   int nlocal = atom->nlocal;
   for (int i = 0; i < nlocal; i++) npartner[i] = 0;
-  maxtouch = 0;
+  maxpartner = 0;
 
   nlocal_neigh = nall_neigh = 0;
   commflag = DEFAULT;
@@ -77,7 +93,7 @@ FixShearHistory::FixShearHistory(LAMMPS *lmp, int narg, char **arg) :
 
 /* ---------------------------------------------------------------------- */
 
-FixShearHistory::~FixShearHistory()
+FixNeighHistory::~FixNeighHistory()
 {
   // unregister this fix so atom class doesn't invoke it any more
 
@@ -86,86 +102,111 @@ FixShearHistory::~FixShearHistory()
 
   // delete locally stored arrays
 
+  delete [] zeroes;
+
+  memory->sfree(firstflag);
+  memory->sfree(firstvalue);
+
   memory->destroy(npartner);
   memory->sfree(partner);
-  memory->sfree(shearpartner);
+  memory->sfree(valuepartner);
+
+  delete [] ipage_atom;
+  delete [] dpage_atom;
+  delete [] ipage_neigh;
+  delete [] dpage_neigh;
 
   // to better detect use-after-delete errors
 
+  firstflag = NULL;
+  firstvalue = NULL;
+
   pair = NULL;
   npartner = NULL;
   partner = NULL;
-  shearpartner = NULL;
+  valuepartner = NULL;
-
-  delete [] ipage;
-  delete [] dpage;
 }
 
 /* ---------------------------------------------------------------------- */
 
-int FixShearHistory::setmask()
+int FixNeighHistory::setmask()
 {
   int mask = 0;
   mask |= PRE_EXCHANGE;
   mask |= MIN_PRE_EXCHANGE;
+  mask |= POST_NEIGHBOR;
+  mask |= MIN_POST_NEIGHBOR;
   return mask;
 }
 
 /* ---------------------------------------------------------------------- */
 
-void FixShearHistory::init()
+void FixNeighHistory::init()
 {
   if (atom->tag_enable == 0)
-    error->all(FLERR,"Granular shear history requires atoms have IDs");
+    error->all(FLERR,"Neighbor history requires atoms have IDs");
 
   allocate_pages();
 }
 
 /* ----------------------------------------------------------------------
    create pages if first time or if neighbor pgsize/oneatom has changed
-   note that latter could cause shear history info to be discarded
+   note that latter could cause neighbor history info to be discarded
 ------------------------------------------------------------------------- */
 
-void FixShearHistory::allocate_pages()
+void FixNeighHistory::allocate_pages()
 {
   int create = 0;
-  if (ipage == NULL) create = 1;
+  if (ipage_atom == NULL) create = 1;
   if (pgsize != neighbor->pgsize) create = 1;
   if (oneatom != neighbor->oneatom) create = 1;
 
   if (create) {
-    delete [] ipage;
+    delete [] ipage_atom;
-    delete [] dpage;
+    delete [] dpage_atom;
+    delete [] ipage_neigh;
+    delete [] dpage_neigh;
 
     pgsize = neighbor->pgsize;
     oneatom = neighbor->oneatom;
     int nmypage = comm->nthreads;
-    ipage = new MyPage<tagint>[nmypage];
+    ipage_atom = new MyPage<tagint>[nmypage];
-    dpage = new MyPage<double>[nmypage];
+    dpage_atom = new MyPage<double>[nmypage];
+    ipage_neigh = new MyPage<int>[nmypage];
+    dpage_neigh = new MyPage<double>[nmypage];
     for (int i = 0; i < nmypage; i++) {
-      ipage[i].init(oneatom,pgsize);
+      ipage_atom[i].init(oneatom,pgsize);
-      dpage[i].init(dnum*oneatom,dnum*pgsize);
+      dpage_atom[i].init(dnum*oneatom,dnum*pgsize);
+      ipage_neigh[i].init(oneatom,pgsize);
+      dpage_neigh[i].init(dnum*oneatom,dnum*pgsize);
     }
   }
 }
 
+/* ---------------------------------------------------------------------- */
+
+void FixNeighHistory::setup_post_neighbor()
+{
+  post_neighbor();
+}
+
 /* ----------------------------------------------------------------------
-   copy shear partner info from neighbor lists to atom arrays
+   copy partner info from neighbor data structs (NDS) to atom arrays
-   should be called whenever neighbor list stores current history info
+   should be called whenever NDS store current history info
-     and need to store the info with owned atoms
+     and need to transfer the info to owned atoms
-   e.g. so atoms can migrate to new procs or between runs
+   e.g. when atoms migrate to new procs, new neigh list built, or between runs
-     when atoms may be added or deleted (neighbor list becomes out-of-date)
+     when atoms may be added or deleted (NDS becomes out-of-date)
-   the next granular neigh list build will put this info back into neigh list
+   the next post_neighbor() will put this info back into new NDS
    called during run before atom exchanges, including for restart files
    called at end of run via post_run()
    do not call during setup of run (setup_pre_exchange)
-     b/c there is no guarantee of a current neigh list (even on continued run)
+     b/c there is no guarantee of a current NDS (even on continued run)
    if run command does a 2nd run with pre = no, then no neigh list
      will be built, but old neigh list will still have the info
    onesided and newton on and newton off versions
 ------------------------------------------------------------------------- */
 
-void FixShearHistory::pre_exchange()
+void FixNeighHistory::pre_exchange()
 {
   if (onesided) pre_exchange_onesided();
   else if (newton_pair) pre_exchange_newton();
@@ -178,60 +219,57 @@ void FixShearHistory::pre_exchange()
    only store history info with spheres
 ------------------------------------------------------------------------- */
 
-void FixShearHistory::pre_exchange_onesided()
+void FixNeighHistory::pre_exchange_onesided()
 {
   int i,j,ii,jj,m,n,inum,jnum;
   int *ilist,*jlist,*numneigh,**firstneigh;
-  int *touch,**firsttouch;
+  int *allflags;
-  double *shear,*allshear,**firstshear;
+  double *allvalues,*onevalues;
 
   // NOTE: all operations until very end are with nlocal_neigh <= current nlocal
   // b/c previous neigh list was built with nlocal_neigh
   // nlocal can be larger if other fixes added atoms at this pre_exchange()
 
-  // zero npartner for owned atoms
+  // clear two paged data structures
-  // clear 2 page data structures
 
-  for (i = 0; i < nlocal_neigh; i++) npartner[i] = 0;
+  ipage_atom->reset();
-
+  dpage_atom->reset();
-  ipage->reset();
-  dpage->reset();
 
   // 1st loop over neighbor list, I = sphere, J = tri
   // only calculate npartner for owned spheres
 
+  for (i = 0; i < nlocal_neigh; i++) npartner[i] = 0;
+
   tagint *tag = atom->tag;
   NeighList *list = pair->list;
   inum = list->inum;
   ilist = list->ilist;
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
-  firsttouch = list->listhistory->firstneigh;
-  firstshear = list->listhistory->firstdouble;
 
   for (ii = 0; ii < inum; ii++) {
     i = ilist[ii];
     jlist = firstneigh[i];
     jnum = numneigh[i];
-    touch = firsttouch[i];
+    allflags = firstflag[i];
 
     for (jj = 0; jj < jnum; jj++)
-      if (touch[jj]) npartner[i]++;
+      if (allflags[jj]) npartner[i]++;
   }
 
-  // get page chunks to store atom IDs and shear history for owned atoms
+  // get page chunks to store partner IDs and values for owned atoms
 
   for (ii = 0; ii < inum; ii++) {
     i = ilist[ii];
     n = npartner[i];
-    partner[i] = ipage->get(n);
+    partner[i] = ipage_atom->get(n);
-    shearpartner[i] = dpage->get(dnum*n);
+    valuepartner[i] = dpage_atom->get(dnum*n);
-    if (partner[i] == NULL || shearpartner[i] == NULL)
+    if (partner[i] == NULL || valuepartner[i] == NULL)
-      error->one(FLERR,"Shear history overflow, boost neigh_modify one");
+      error->one(FLERR,"Neighbor history overflow, boost neigh_modify one");
   }
 
   // 2nd loop over neighbor list, I = sphere, J = tri
-  // store atom IDs and shear history for owned spheres
+  // store partner IDs and values for owned+ghost atoms
   // re-zero npartner to use as counter
 
   for (i = 0; i < nlocal_neigh; i++) npartner[i] = 0;
@@ -239,28 +277,28 @@ void FixShearHistory::pre_exchange_onesided()
   for (ii = 0; ii < inum; ii++) {
     i = ilist[ii];
     jlist = firstneigh[i];
-    allshear = firstshear[i];
     jnum = numneigh[i];
-    touch = firsttouch[i];
+    allflags = firstflag[i];
+    allvalues = firstvalue[i];
 
     for (jj = 0; jj < jnum; jj++) {
-      if (touch[jj]) {
+      if (allflags[jj]) {
-        shear = &allshear[dnum*jj];
+        onevalues = &allvalues[dnum*jj];
         j = jlist[jj];
         j &= NEIGHMASK;
         m = npartner[i]++;
         partner[i][m] = tag[j];
-        memcpy(&shearpartner[i][dnum*m],shear,dnumbytes);
+        memcpy(&valuepartner[i][dnum*m],onevalues,dnumbytes);
       }
     }
   }
 
-  // set maxtouch = max # of partners of any owned atom
+  // set maxpartner = max # of partners of any owned atom
   // bump up comm->maxexchange_fix if necessary
   
-  maxtouch = 0;
+  maxpartner = 0;
-  for (i = 0; i < nlocal_neigh; i++) maxtouch = MAX(maxtouch,npartner[i]);
+  for (i = 0; i < nlocal_neigh; i++) maxpartner = MAX(maxpartner,npartner[i]);
-  comm->maxexchange_fix = MAX(comm->maxexchange_fix,(dnum+1)*maxtouch+1);
+  comm->maxexchange_fix = MAX(comm->maxexchange_fix,(dnum+1)*maxpartner+1);
 
   // zero npartner values from previous nlocal_neigh to current nlocal
 
@@ -269,50 +307,47 @@ void FixShearHistory::pre_exchange_onesided()
 }
 
 /* ----------------------------------------------------------------------
-   newton on version, for sphere/sphere contacts
+   newton ON version
-   performs reverse comm to acquire shear partner info from ghost atoms
+   performs reverse comm to acquire partner values from ghost atoms
 ------------------------------------------------------------------------- */
 
-void FixShearHistory::pre_exchange_newton()
+void FixNeighHistory::pre_exchange_newton()
 {
   int i,j,ii,jj,m,n,inum,jnum;
   int *ilist,*jlist,*numneigh,**firstneigh;
-  int *touch,**firsttouch;
+  int *allflags;
-  double *shear,*shearj,*allshear,**firstshear;
+  double *allvalues,*onevalues,*jvalues;
 
   // NOTE: all operations until very end are with 
   //   nlocal_neigh  <= current nlocal and nall_neigh
-  // b/c previous neigh list was built with nlocal_neigh,nghost_neigh
+  // b/c previous neigh list was built with nlocal_neigh & nghost_neigh
   // nlocal can be larger if other fixes added atoms at this pre_exchange()
 
-  // zero npartner for owned+ghost atoms
+  // clear two paged data structures
-  // clear 2 page data structures
 
-  for (i = 0; i < nall_neigh; i++) npartner[i] = 0;
+  ipage_atom->reset();
-
+  dpage_atom->reset();
-  ipage->reset();
-  dpage->reset();
 
   // 1st loop over neighbor list
   // calculate npartner for owned+ghost atoms
 
+  for (i = 0; i < nall_neigh; i++) npartner[i] = 0;
+
   tagint *tag = atom->tag;
   NeighList *list = pair->list;
   inum = list->inum;
   ilist = list->ilist;
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
-  firsttouch = list->listhistory->firstneigh;
-  firstshear = list->listhistory->firstdouble;
 
   for (ii = 0; ii < inum; ii++) {
     i = ilist[ii];
     jlist = firstneigh[i];
     jnum = numneigh[i];
-    touch = firsttouch[i];
+    allflags = firstflag[i];
 
     for (jj = 0; jj < jnum; jj++) {
-      if (touch[jj]) {
+      if (allflags[jj]) {
         npartner[i]++;
         j = jlist[jj];
         j &= NEIGHMASK;
@@ -326,29 +361,29 @@ void FixShearHistory::pre_exchange_newton()
   commflag = NPARTNER;
   comm->reverse_comm_fix(this,0);
 
-  // get page chunks to store atom IDs and shear history for owned+ghost atoms
+  // get page chunks to store partner IDs and values for owned+ghost atoms
 
   for (ii = 0; ii < inum; ii++) {
     i = ilist[ii];
     n = npartner[i];
-    partner[i] = ipage->get(n);
+    partner[i] = ipage_atom->get(n);
-    shearpartner[i] = dpage->get(dnum*n);
+    valuepartner[i] = dpage_atom->get(dnum*n);
-    if (partner[i] == NULL || shearpartner[i] == NULL) {
+    if (partner[i] == NULL || valuepartner[i] == NULL) {
-      error->one(FLERR,"Shear history overflow, boost neigh_modify one");
+      error->one(FLERR,"Neighbor history overflow, boost neigh_modify one");
     }
   }
 
   for (i = nlocal_neigh; i < nall_neigh; i++) {
     n = npartner[i];
-    partner[i] = ipage->get(n);
+    partner[i] = ipage_atom->get(n);
-    shearpartner[i] = dpage->get(dnum*n);
+    valuepartner[i] = dpage_atom->get(dnum*n);
-    if (partner[i] == NULL || shearpartner[i] == NULL) {
+    if (partner[i] == NULL || valuepartner[i] == NULL) {
-      error->one(FLERR,"Shear history overflow, boost neigh_modify one");
+      error->one(FLERR,"Neighbor history overflow, boost neigh_modify one");
     }
   }
 
   // 2nd loop over neighbor list
-  // store atom IDs and shear history for owned+ghost atoms
+  // store partner IDs and values for owned+ghost atoms
   // re-zero npartner to use as counter
 
   for (i = 0; i < nall_neigh; i++) npartner[i] = 0;
@@ -356,40 +391,40 @@ void FixShearHistory::pre_exchange_newton()
   for (ii = 0; ii < inum; ii++) {
     i = ilist[ii];
     jlist = firstneigh[i];
-    allshear = firstshear[i];
     jnum = numneigh[i];
-    touch = firsttouch[i];
+    allflags = firstflag[i];
+    allvalues = firstvalue[i];
 
     for (jj = 0; jj < jnum; jj++) {
-      if (touch[jj]) {
+      if (allflags[jj]) {
-        shear = &allshear[dnum*jj];
+        onevalues = &allvalues[dnum*jj];
         j = jlist[jj];
         j &= NEIGHMASK;
         m = npartner[i]++;
         partner[i][m] = tag[j];
-        memcpy(&shearpartner[i][dnum*m],shear,dnumbytes);
+        memcpy(&valuepartner[i][dnum*m],onevalues,dnumbytes);
         m = npartner[j]++;
         partner[j][m] = tag[i];
-        shearj = &shearpartner[j][dnum*m];
+        jvalues = &valuepartner[j][dnum*m];
-        for (n = 0; n < dnum; n++) shearj[n] = -shear[n];
+        for (n = 0; n < dnum; n++) jvalues[n] = -onevalues[n];
       }
     }
   }
 
   // perform reverse comm to augment
-  // owned atom partner/shearpartner with ghost info
+  // owned atom partner/valuepartner with ghost info
   // use variable variant b/c size of packed data can be arbitrarily large
   //   if many touching neighbors for large particle
 
   commflag = PERPARTNER;
   comm->reverse_comm_fix_variable(this);
 
-  // set maxtouch = max # of partners of any owned atom
+  // set maxpartner = max # of partners of any owned atom
   // bump up comm->maxexchange_fix if necessary
 
-  maxtouch = 0;
+  maxpartner = 0;
-  for (i = 0; i < nlocal_neigh; i++) maxtouch = MAX(maxtouch,npartner[i]);
+  for (i = 0; i < nlocal_neigh; i++) maxpartner = MAX(maxpartner,npartner[i]);
-  comm->maxexchange_fix = MAX(comm->maxexchange_fix,(dnum+1)*maxtouch+1);
+  comm->maxexchange_fix = MAX(comm->maxexchange_fix,4*maxpartner+1);
 
   // zero npartner values from previous nlocal_neigh to current nlocal
 
@@ -398,49 +433,47 @@ void FixShearHistory::pre_exchange_newton()
 }
 
 /* ----------------------------------------------------------------------
-   newton off version, for sphere/sphere contacts
+   newton OFF version
-   newton OFF works with smaller vectors that don't include ghost info
+   do not need partner values from ghost atoms
+   assume J values are negative of I values
 ------------------------------------------------------------------------- */
 
-void FixShearHistory::pre_exchange_no_newton()
+void FixNeighHistory::pre_exchange_no_newton()
 {
   int i,j,ii,jj,m,n,inum,jnum;
   int *ilist,*jlist,*numneigh,**firstneigh;
-  int *touch,**firsttouch;
+  int *allflags;
-  double *shear,*shearj,*allshear,**firstshear;
+  double *allvalues,*onevalues,*jvalues;
 
   // NOTE: all operations until very end are with nlocal_neigh <= current nlocal
   // b/c previous neigh list was built with nlocal_neigh
   // nlocal can be larger if other fixes added atoms at this pre_exchange()
 
-  // zero npartner for owned atoms
+  // clear two paged data structures
-  // clear 2 page data structures
 
-  for (i = 0; i < nlocal_neigh; i++) npartner[i] = 0;
+  ipage_atom->reset();
-
+  dpage_atom->reset();
-  ipage->reset();
-  dpage->reset();
 
   // 1st loop over neighbor list
   // calculate npartner for owned atoms
 
+  for (i = 0; i < nlocal_neigh; i++) npartner[i] = 0;
+
   tagint *tag = atom->tag;
   NeighList *list = pair->list;
   inum = list->inum;
   ilist = list->ilist;
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
-  firsttouch = list->listhistory->firstneigh;
-  firstshear = list->listhistory->firstdouble;
 
   for (ii = 0; ii < inum; ii++) {
     i = ilist[ii];
     jlist = firstneigh[i];
     jnum = numneigh[i];
-    touch = firsttouch[i];
+    allflags = firstflag[i];
 
     for (jj = 0; jj < jnum; jj++) {
-      if (touch[jj]) {
+      if (allflags[jj]) {
         npartner[i]++;
         j = jlist[jj];
         j &= NEIGHMASK;
@@ -449,19 +482,19 @@ void FixShearHistory::pre_exchange_no_newton()
     }
   }
 
-  // get page chunks to store atom IDs and shear history for owned atoms
+  // get page chunks to store partner IDs and values for owned atoms
 
   for (ii = 0; ii < inum; ii++) {
     i = ilist[ii];
     n = npartner[i];
-    partner[i] = ipage->get(n);
+    partner[i] = ipage_atom->get(n);
-    shearpartner[i] = dpage->get(dnum*n);
+    valuepartner[i] = dpage_atom->get(dnum*n);
-    if (partner[i] == NULL || shearpartner[i] == NULL)
+    if (partner[i] == NULL || valuepartner[i] == NULL)
-      error->one(FLERR,"Shear history overflow, boost neigh_modify one");
+      error->one(FLERR,"Neighbor history overflow, boost neigh_modify one");
   }
 
   // 2nd loop over neighbor list
-  // store atom IDs and shear history for owned atoms
+  // store partner IDs and values for owned+ghost atoms
   // re-zero npartner to use as counter
 
   for (i = 0; i < nlocal_neigh; i++) npartner[i] = 0;
@@ -469,34 +502,34 @@ void FixShearHistory::pre_exchange_no_newton()
   for (ii = 0; ii < inum; ii++) {
     i = ilist[ii];
     jlist = firstneigh[i];
-    allshear = firstshear[i];
     jnum = numneigh[i];
-    touch = firsttouch[i];
+    allflags = firstflag[i];
+    allvalues = firstvalue[i];
 
     for (jj = 0; jj < jnum; jj++) {
-      if (touch[jj]) {
+      if (allflags[jj]) {
-        shear = &allshear[dnum*jj];
+        onevalues = &allvalues[dnum*jj];
         j = jlist[jj];
         j &= NEIGHMASK;
         m = npartner[i]++;
         partner[i][m] = tag[j];
-        memcpy(&shearpartner[i][dnum*m],shear,dnumbytes);
+        memcpy(&valuepartner[i][dnum*m],onevalues,dnumbytes);
         if (j < nlocal_neigh) {
           m = npartner[j]++;
           partner[j][m] = tag[i];
-          shearj = &shearpartner[j][dnum*m];
+          jvalues = &valuepartner[j][dnum*m];
-          for (n = 0; n < dnum; n++) shearj[n] = -shear[n];
+          for (n = 0; n < dnum; n++) jvalues[n] = -onevalues[n];
         }
       }
     }
   }
 
-  // set maxtouch = max # of partners of any owned atom
+  // set maxpartner = max # of partners of any owned atom
   // bump up comm->maxexchange_fix if necessary
   
-  maxtouch = 0;
+  maxpartner = 0;
-  for (i = 0; i < nlocal_neigh; i++) maxtouch = MAX(maxtouch,npartner[i]);
+  for (i = 0; i < nlocal_neigh; i++) maxpartner = MAX(maxpartner,npartner[i]);
-  comm->maxexchange_fix = MAX(comm->maxexchange_fix,(dnum+1)*maxtouch+1);
+  comm->maxexchange_fix = MAX(comm->maxexchange_fix,(dnum+1)*maxpartner+1);
 
   // zero npartner values from previous nlocal_neigh to current nlocal
 
@@ -506,14 +539,109 @@ void FixShearHistory::pre_exchange_no_newton()
 
 /* ---------------------------------------------------------------------- */
 
-void FixShearHistory::min_pre_exchange()
+void FixNeighHistory::min_pre_exchange()
 {
   pre_exchange();
 }
 
+/* ----------------------------------------------------------------------
+   called after neighbor list is build
+   recover history info stored temporarily in per-atom partner lists
+     and store afresh in per-neighbor firstflag and firstvalue lists
+------------------------------------------------------------------------- */
+
+void FixNeighHistory::post_neighbor()
+{
+  int i,j,m,ii,jj,nn,np,inum,jnum,rflag;
+  tagint jtag;
+  double xtmp,ytmp,ztmp,delx,dely,delz;
+  double radi,rsq,radsum;
+  int *ilist,*jlist,*numneigh,**firstneigh;
+  int *allflags;
+  double *allvalues;
+
+  // store atom counts used for new neighbor list which was just built
+
+  int nlocal = atom->nlocal;
+  int nall = nlocal + atom->nghost;
+  nlocal_neigh = nlocal;
+  nall_neigh = nall;
+
+  // realloc firstflag and firstvalue if needed
+
+  if (maxatom < nlocal) {
+    memory->sfree(firstflag);
+    memory->sfree(firstvalue);
+    maxatom = nall;
+    firstflag = (int **) 
+      memory->smalloc(maxatom*sizeof(int *),"neighbor_history:firstflag");
+    firstvalue = (double **) 
+      memory->smalloc(maxatom*sizeof(double *),"neighbor_history:firstvalue");
+  }
+
+  // loop over newly built neighbor list
+  // repopulate entire per-neighbor data structs
+  //   whether with old-neigh partner info or zeroes
+
+  ipage_neigh->reset();
+  dpage_neigh->reset();
+
+  tagint *tag = atom->tag;
+  NeighList *list = pair->list;
+  inum = list->inum;
+  ilist = list->ilist;
+  numneigh = list->numneigh;
+  firstneigh = list->firstneigh;
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+    firstflag[i] = allflags = ipage_neigh->get(jnum);
+    firstvalue[i] = allvalues = dpage_neigh->get(jnum*dnum);
+    np = npartner[i];
+    nn = 0;
+
+    for (jj = 0; jj < jnum; jj++) {
+      j = jlist[jj];
+      rflag = sbmask(j);
+      j &= NEIGHMASK;
+      jlist[jj] = j;
+
+      // rflag = 1 if r < radsum in npair_size() method
+      // preserve neigh history info if tag[j] is in old-neigh partner list
+      // this test could be more geometrically precise for two sphere/line/tri
+
+      if (rflag) {
+        jtag = tag[j];
+        for (m = 0; m < np; m++)
+          if (partner[i][m] == jtag) break;
+        if (m < np) {
+          allflags[jj] = 1;
+          memcpy(&allvalues[nn],&valuepartner[i][dnum*m],dnumbytes);
+        } else {
+          allflags[jj] = 0;
+          memcpy(&allvalues[nn],zeroes,dnumbytes);
+        }
+      } else {
+        allflags[jj] = 0;
+        memcpy(&allvalues[nn],zeroes,dnumbytes);
+      }
+      nn += dnum;
+    }
+  }
+}
+
 /* ---------------------------------------------------------------------- */
 
-void FixShearHistory::post_run()
+void FixNeighHistory::min_post_neighbor()
+{
+  post_neighbor();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixNeighHistory::post_run()
 {
   pre_exchange();
 }
@@ -522,17 +650,21 @@ void FixShearHistory::post_run()
    memory usage of local atom-based arrays
 ------------------------------------------------------------------------- */
 
-double FixShearHistory::memory_usage()
+double FixNeighHistory::memory_usage()
 {
   int nmax = atom->nmax;
-  double bytes = nmax * sizeof(int);
+  double bytes = nmax * sizeof(int);    // npartner
-  bytes += nmax * sizeof(tagint *);
+  bytes += nmax * sizeof(tagint *);     // partner
-  bytes += nmax * sizeof(double *);
+  bytes += nmax * sizeof(double *);     // valuepartner
+  bytes += maxatom * sizeof(int *);     // firstflag
+  bytes += maxatom * sizeof(double *);  // firstvalue
 
   int nmypage = comm->nthreads;
   for (int i = 0; i < nmypage; i++) {
-    bytes += ipage[i].size();
+    bytes += ipage_atom[i].size();
-    bytes += dpage[i].size();
+    bytes += dpage_atom[i].size();
+    bytes += ipage_neigh[i].size();
+    bytes += dpage_neigh[i].size();
   }
 
   return bytes;
@@ -542,38 +674,38 @@ double FixShearHistory::memory_usage()
    allocate local atom-based arrays
 ------------------------------------------------------------------------- */
 
-void FixShearHistory::grow_arrays(int nmax)
+void FixNeighHistory::grow_arrays(int nmax)
 {
-  memory->grow(npartner,nmax,"shear_history:npartner");
+  memory->grow(npartner,nmax,"neighbor_history:npartner");
   partner = (tagint **) memory->srealloc(partner,nmax*sizeof(tagint *),
-                                         "shear_history:partner");
+                                         "neighbor_history:partner");
-  shearpartner = (double **) memory->srealloc(shearpartner,
+  valuepartner = (double **) memory->srealloc(valuepartner,
                                               nmax*sizeof(double *),
-                                              "shear_history:shearpartner");
+                                              "neighbor_history:valuepartner");
 }
 
 /* ----------------------------------------------------------------------
    copy values within local atom-based arrays
 ------------------------------------------------------------------------- */
 
-void FixShearHistory::copy_arrays(int i, int j, int delflag)
+void FixNeighHistory::copy_arrays(int i, int j, int delflag)
 {
-  // just copy pointers for partner and shearpartner
+  // just copy pointers for partner and valuepartner
-  // b/c can't overwrite chunk allocation inside ipage,dpage
+  // b/c can't overwrite chunk allocation inside ipage_atom,dpage_atom
   // incoming atoms in unpack_exchange just grab new chunks
   // so are orphaning chunks for migrating atoms
-  // OK, b/c will reset ipage,dpage on next reneighboring
+  // OK, b/c will reset ipage_atom,dpage_atom on next reneighboring
 
   npartner[j] = npartner[i];
   partner[j] = partner[i];
-  shearpartner[j] = shearpartner[i];
+  valuepartner[j] = valuepartner[i];
 }
 
 /* ----------------------------------------------------------------------
    initialize one atom's array values, called when atom is created
 ------------------------------------------------------------------------- */
 
-void FixShearHistory::set_arrays(int i)
+void FixNeighHistory::set_arrays(int i)
 {
   npartner[i] = 0;
 }
@@ -582,7 +714,7 @@ void FixShearHistory::set_arrays(int i)
    only called by Comm::reverse_comm_fix_variable for PERPARTNER mode
 ------------------------------------------------------------------------- */
 
-int FixShearHistory::pack_reverse_comm_size(int n, int first)
+int FixNeighHistory::pack_reverse_comm_size(int n, int first)
 {
   int i,last;
 
@@ -590,7 +722,7 @@ int FixShearHistory::pack_reverse_comm_size(int n, int first)
   last = first + n;
 
   for (i = first; i < last; i++)
-    m += 1 + (dnum+1)*npartner[i];
+    m += 1 + 4*npartner[i];
 
   return m;
 }
@@ -599,7 +731,7 @@ int FixShearHistory::pack_reverse_comm_size(int n, int first)
    two modes: NPARTNER and PERPARTNER
 ------------------------------------------------------------------------- */
 
-int FixShearHistory::pack_reverse_comm(int n, int first, double *buf)
+int FixNeighHistory::pack_reverse_comm(int n, int first, double *buf)
 {
   int i,k,last;
 
@@ -615,11 +747,11 @@ int FixShearHistory::pack_reverse_comm(int n, int first, double *buf)
       buf[m++] = npartner[i];
       for (k = 0; k < npartner[i]; k++) {
         buf[m++] = partner[i][k];
-        memcpy(&buf[m],&shearpartner[i][dnum*k],dnumbytes);
+        memcpy(&buf[m],&valuepartner[i][dnum*k],dnumbytes);
         m += dnum;
       }
     }
-  } else error->all(FLERR,"Unsupported comm mode in shear history");
+  } else error->all(FLERR,"Unsupported comm mode in neighbor history");
 
   return m;
 }
@@ -628,7 +760,7 @@ int FixShearHistory::pack_reverse_comm(int n, int first, double *buf)
    two modes: NPARTNER and PERPARTNER
 ------------------------------------------------------------------------- */
 
-void FixShearHistory::unpack_reverse_comm(int n, int *list, double *buf)
+void FixNeighHistory::unpack_reverse_comm(int n, int *list, double *buf)
 {
   int i,j,k,kk,ncount;
 
@@ -646,18 +778,18 @@ void FixShearHistory::unpack_reverse_comm(int n, int *list, double *buf)
       for (k = 0; k < ncount; k++) {
         kk = npartner[j]++;
         partner[j][kk] = static_cast<tagint> (buf[m++]);
-        memcpy(&shearpartner[j][dnum*kk],&buf[m],dnumbytes);
+        memcpy(&valuepartner[j][dnum*kk],&buf[m],dnumbytes);
         m += dnum;
       }
     }
-  } else error->all(FLERR,"Unsupported comm mode in shear history");
+  } else error->all(FLERR,"Unsupported comm mode in neighbor history");
 }
 
 /* ----------------------------------------------------------------------
    pack values in local atom-based arrays for exchange with another proc
 ------------------------------------------------------------------------- */
 
-int FixShearHistory::pack_exchange(int i, double *buf)
+int FixNeighHistory::pack_exchange(int i, double *buf)
 {
   // NOTE: how do I know comm buf is big enough if extreme # of touching neighs
   // Comm::BUFEXTRA may need to be increased
@@ -666,7 +798,7 @@ int FixShearHistory::pack_exchange(int i, double *buf)
   buf[m++] = npartner[i];
   for (int n = 0; n < npartner[i]; n++) {
     buf[m++] = partner[i][n];
-    memcpy(&buf[m],&shearpartner[i][dnum*n],dnumbytes);
+    memcpy(&buf[m],&valuepartner[i][dnum*n],dnumbytes);
     m += dnum;
   }
   return m;
@@ -676,18 +808,18 @@ int FixShearHistory::pack_exchange(int i, double *buf)
    unpack values in local atom-based arrays from exchange with another proc
 ------------------------------------------------------------------------- */
 
-int FixShearHistory::unpack_exchange(int nlocal, double *buf)
+int FixNeighHistory::unpack_exchange(int nlocal, double *buf)
 {
-  // allocate new chunks from ipage,dpage for incoming values
+  // allocate new chunks from ipage_atom,dpage_atom for incoming values
 
   int m = 0;
   npartner[nlocal] = static_cast<int> (buf[m++]);
-  maxtouch = MAX(maxtouch,npartner[nlocal]);
+  maxpartner = MAX(maxpartner,npartner[nlocal]);
-  partner[nlocal] = ipage->get(npartner[nlocal]);
+  partner[nlocal] = ipage_atom->get(npartner[nlocal]);
-  shearpartner[nlocal] = dpage->get(dnum*npartner[nlocal]);
+  valuepartner[nlocal] = dpage_atom->get(dnum*npartner[nlocal]);
   for (int n = 0; n < npartner[nlocal]; n++) {
     partner[nlocal][n] = static_cast<tagint> (buf[m++]);
-    memcpy(&shearpartner[nlocal][dnum*n],&buf[m],dnumbytes);
+    memcpy(&valuepartner[nlocal][dnum*n],&buf[m],dnumbytes);
     m += dnum;
   }
   return m;
@@ -697,13 +829,13 @@ int FixShearHistory::unpack_exchange(int nlocal, double *buf)
    pack values in local atom-based arrays for restart file
 ------------------------------------------------------------------------- */
 
-int FixShearHistory::pack_restart(int i, double *buf)
+int FixNeighHistory::pack_restart(int i, double *buf)
 {
   int m = 1;
   buf[m++] = npartner[i];
   for (int n = 0; n < npartner[i]; n++) {
     buf[m++] = partner[i][n];
-    memcpy(&buf[m],&shearpartner[i][dnum*n],dnumbytes);
+    memcpy(&buf[m],&valuepartner[i][dnum*n],dnumbytes);
     m += dnum;
   }
   buf[0] = m;
@@ -714,11 +846,11 @@ int FixShearHistory::pack_restart(int i, double *buf)
    unpack values from atom->extra array to restart the fix
 ------------------------------------------------------------------------- */
 
-void FixShearHistory::unpack_restart(int nlocal, int nth)
+void FixNeighHistory::unpack_restart(int nlocal, int nth)
 {
-  // ipage = NULL if being called from granular pair style init()
+  // ipage_atom = NULL if being called from granular pair style init()
 
-  if (ipage == NULL) allocate_pages();
+  if (ipage_atom == NULL) allocate_pages();
 
   // skip to Nth set of extra values
 
@@ -728,15 +860,15 @@ void FixShearHistory::unpack_restart(int nlocal, int nth)
   for (int i = 0; i < nth; i++) m += static_cast<int> (extra[nlocal][m]);
   m++;
 
-  // allocate new chunks from ipage,dpage for incoming values
+  // allocate new chunks from ipage_atom,dpage_atom for incoming values
 
   npartner[nlocal] = static_cast<int> (extra[nlocal][m++]);
-  maxtouch = MAX(maxtouch,npartner[nlocal]);
+  maxpartner = MAX(maxpartner,npartner[nlocal]);
-  partner[nlocal] = ipage->get(npartner[nlocal]);
+  partner[nlocal] = ipage_atom->get(npartner[nlocal]);
-  shearpartner[nlocal] = dpage->get(dnum*npartner[nlocal]);
+  valuepartner[nlocal] = dpage_atom->get(dnum*npartner[nlocal]);
   for (int n = 0; n < npartner[nlocal]; n++) {
     partner[nlocal][n] = static_cast<tagint> (extra[nlocal][m++]);
-    memcpy(&shearpartner[nlocal][dnum*n],&extra[nlocal][m],dnumbytes);
+    memcpy(&valuepartner[nlocal][dnum*n],&extra[nlocal][m],dnumbytes);
     m += dnum;
   }
 }
@@ -745,20 +877,20 @@ void FixShearHistory::unpack_restart(int nlocal, int nth)
    maxsize of any atom's restart data
 ------------------------------------------------------------------------- */
 
-int FixShearHistory::maxsize_restart()
+int FixNeighHistory::maxsize_restart()
 {
-  // maxtouch_all = max # of touching partners across all procs
+  // maxpartner_all = max # of touching partners across all procs
 
-  int maxtouch_all;
+  int maxpartner_all;
-  MPI_Allreduce(&maxtouch,&maxtouch_all,1,MPI_INT,MPI_MAX,world);
+  MPI_Allreduce(&maxpartner,&maxpartner_all,1,MPI_INT,MPI_MAX,world);
-  return (dnum+1)*maxtouch_all + 2;
+  return (dnum+1)*maxpartner_all + 2;
 }
 
 /* ----------------------------------------------------------------------
    size of atom nlocal's restart data
 ------------------------------------------------------------------------- */
 
-int FixShearHistory::size_restart(int nlocal)
+int FixNeighHistory::size_restart(int nlocal)
 {
   return (dnum+1)*npartner[nlocal] + 2;
 }

src/fix_neigh_history.h

@@ -13,38 +13,35 @@
 
 #ifdef FIX_CLASS
 
-FixStyle(SHEAR_HISTORY,FixShearHistory)
+FixStyle(NEIGH_HISTORY,FixNeighHistory)
 
 #else
 
-#ifndef LMP_FIX_SHEAR_HISTORY_H
+#ifndef LMP_FIX_NEIGH_HISTORY_H
-#define LMP_FIX_SHEAR_HISTORY_H
+#define LMP_FIX_NEIGH_HISTORY_H
 
 #include "fix.h"
 #include "my_page.h"
 
 namespace LAMMPS_NS {
 
-class FixShearHistory : public Fix {
+class FixNeighHistory : public Fix {
-  //friend class Neighbor;
-  //friend class PairGranHookeHistory;
-  friend class PairLineGranHookeHistory;
-  friend class PairTriGranHookeHistory;
-
  public:
   int nlocal_neigh;             // nlocal at last time neigh list was built
   int nall_neigh;               // ditto for nlocal+nghost
-  int *npartner;                // # of touching partners of each atom
+  int **firstflag;              // ptr to each atom's neighbor flsg
-  tagint **partner;             // global atom IDs for the partners
+  double **firstvalue;          // ptr to each atom's values
-  double **shearpartner;        // shear values with the partner
+  class Pair *pair;             // ptr to pair style that uses neighbor history
-  class Pair *pair;             // ptr to pair style that uses shear history
 
-  FixShearHistory(class LAMMPS *, int, char **);
+  FixNeighHistory(class LAMMPS *, int, char **);
-  ~FixShearHistory();
+  ~FixNeighHistory();
   int setmask();
   void init();
+  void setup_post_neighbor();
   virtual void pre_exchange();
   void min_pre_exchange();
+  void post_neighbor();
+  void min_post_neighbor();
   void post_run();
 
   double memory_usage();
@@ -64,20 +61,40 @@ class FixShearHistory : public Fix {
 
  protected:
   int newton_pair;              // same as force setting
-  int dnum,dnumbytes;           // dnum = # of shear history values
+  int dnum,dnumbytes;           // dnum = # of values per neighbor
   int onesided;                 // 1 for line/tri history, else 0
 
-  int maxtouch;                 // max # of touching partners for my atoms
+  int maxatom;                  // max size of firstflag and firstvalue
   int commflag;                 // mode of reverse comm to get ghost info
+  double *zeroes;
+
+  // per-atom data structures
+  // partners = flagged neighbors of an atom
+
+  int *npartner;                // # of partners of each atom
+  tagint **partner;             // global atom IDs for the partners
+  double **valuepartner;        // values for the partners
+  int maxpartner;               // max # of partners for any of my atoms
+
+  // per-atom data structs pointed to by partner & valuepartner
 
   int pgsize,oneatom;           // copy of settings in Neighbor
-  MyPage<tagint> *ipage;        // pages of partner atom IDs
+  MyPage<tagint> *ipage_atom;   // pages of partner atom IDs
-  MyPage<double> *dpage;        // pages of shear history with partners
+  MyPage<double> *dpage_atom;   // pages of partner values
+
+  // per-neighbor data structs pointed to by firstflag & firstvalue
+
+  MyPage<int> *ipage_neigh;     // pages of local atom indices
+  MyPage<double> *dpage_neigh;  // pages of partner values
 
   void pre_exchange_onesided();
   void pre_exchange_newton();
   void pre_exchange_no_newton();
   void allocate_pages();
+
+  inline int sbmask(int j) const {
+    return j >> SBBITS & 3;
+  }
 };
 
 }

src/modify.cpp

@@ -42,7 +42,7 @@ Modify::Modify(LAMMPS *lmp) : Pointers(lmp)
 {
   nfix = maxfix = 0;
   n_initial_integrate = n_post_integrate = 0;
-  n_pre_exchange = n_pre_neighbor = 0;
+  n_pre_exchange = n_pre_neighbor = n_post_neighbor = 0;
   n_pre_force = n_pre_reverse = n_post_force = 0;
   n_final_integrate = n_end_of_step = n_thermo_energy = 0;
   n_thermo_energy_atom = 0;
@@ -54,14 +54,14 @@ Modify::Modify(LAMMPS *lmp) : Pointers(lmp)
   fix = NULL;
   fmask = NULL;
   list_initial_integrate = list_post_integrate = NULL;
-  list_pre_exchange = list_pre_neighbor = NULL;
+  list_pre_exchange = list_pre_neighbor = list_post_neighbor = NULL;
   list_pre_force = list_pre_reverse = list_post_force = NULL;
   list_final_integrate = list_end_of_step = NULL;
   list_thermo_energy = list_thermo_energy_atom = NULL;
   list_initial_integrate_respa = list_post_integrate_respa = NULL;
   list_pre_force_respa = list_post_force_respa = NULL;
   list_final_integrate_respa = NULL;
-  list_min_pre_exchange = list_min_pre_neighbor = NULL;
+  list_min_pre_exchange = list_min_pre_neighbor = list_min_post_neighbor = NULL;
   list_min_pre_force = list_min_pre_reverse = list_min_post_force = NULL;
   list_min_energy = NULL;
 
@@ -123,6 +123,7 @@ Modify::~Modify()
   delete [] list_post_integrate;
   delete [] list_pre_exchange;
   delete [] list_pre_neighbor;
+  delete [] list_post_neighbor;
   delete [] list_pre_force;
   delete [] list_pre_reverse;
   delete [] list_post_force;
@@ -137,6 +138,7 @@ Modify::~Modify()
   delete [] list_final_integrate_respa;
   delete [] list_min_pre_exchange;
   delete [] list_min_pre_neighbor;
+  delete [] list_min_post_neighbor;
   delete [] list_min_pre_force;
   delete [] list_min_pre_reverse;
   delete [] list_min_post_force;
@@ -169,6 +171,7 @@ void Modify::init()
   list_init(POST_INTEGRATE,n_post_integrate,list_post_integrate);
   list_init(PRE_EXCHANGE,n_pre_exchange,list_pre_exchange);
   list_init(PRE_NEIGHBOR,n_pre_neighbor,list_pre_neighbor);
+  list_init(POST_NEIGHBOR,n_post_neighbor,list_post_neighbor);
   list_init(PRE_FORCE,n_pre_force,list_pre_force);
   list_init(PRE_REVERSE,n_pre_reverse,list_pre_reverse);
   list_init(POST_FORCE,n_post_force,list_post_force);
@@ -190,6 +193,7 @@ void Modify::init()
 
   list_init(MIN_PRE_EXCHANGE,n_min_pre_exchange,list_min_pre_exchange);
   list_init(MIN_PRE_NEIGHBOR,n_min_pre_neighbor,list_min_pre_neighbor);
+  list_init(MIN_POST_NEIGHBOR,n_min_post_neighbor,list_min_post_neighbor);
   list_init(MIN_PRE_FORCE,n_min_pre_force,list_min_pre_force);
   list_init(MIN_PRE_REVERSE,n_min_pre_reverse,list_min_pre_reverse);
   list_init(MIN_POST_FORCE,n_min_post_force,list_min_post_force);
@@ -329,6 +333,21 @@ void Modify::setup_pre_neighbor()
       fix[list_min_pre_neighbor[i]]->setup_pre_neighbor();
 }
 
+/* ----------------------------------------------------------------------
+   setup post_neighbor call, only for fixes that define post_neighbor
+   called from Verlet, RESPA
+------------------------------------------------------------------------- */
+
+void Modify::setup_post_neighbor()
+{
+  if (update->whichflag == 1)
+    for (int i = 0; i < n_post_neighbor; i++)
+      fix[list_post_neighbor[i]]->setup_post_neighbor();
+  else if (update->whichflag == 2)
+    for (int i = 0; i < n_min_post_neighbor; i++)
+      fix[list_min_post_neighbor[i]]->setup_post_neighbor();
+}
+
 /* ----------------------------------------------------------------------
    setup pre_force call, only for fixes that define pre_force
    called from Verlet, RESPA, Min
@@ -399,6 +418,16 @@ void Modify::pre_neighbor()
     fix[list_pre_neighbor[i]]->pre_neighbor();
 }
 
+/* ----------------------------------------------------------------------
+   post_neighbor call, only for relevant fixes
+------------------------------------------------------------------------- */
+
+void Modify::post_neighbor()
+{
+  for (int i = 0; i < n_post_neighbor; i++)
+    fix[list_post_neighbor[i]]->post_neighbor();
+}
+
 /* ----------------------------------------------------------------------
    pre_force call, only for relevant fixes
 ------------------------------------------------------------------------- */
@@ -589,6 +618,16 @@ void Modify::min_pre_neighbor()
     fix[list_min_pre_neighbor[i]]->min_pre_neighbor();
 }
 
+/* ----------------------------------------------------------------------
+   minimizer post-neighbor call, only for relevant fixes
+------------------------------------------------------------------------- */
+
+void Modify::min_post_neighbor()
+{
+  for (int i = 0; i < n_min_post_neighbor; i++)
+    fix[list_min_post_neighbor[i]]->min_post_neighbor();
+}
+
 /* ----------------------------------------------------------------------
    minimizer pre-force call, only for relevant fixes
 ------------------------------------------------------------------------- */

src/modify.h

@@ -29,12 +29,13 @@ class Modify : protected Pointers {
 
  public:
   int nfix,maxfix;
-  int n_initial_integrate,n_post_integrate,n_pre_exchange,n_pre_neighbor;
+  int n_initial_integrate,n_post_integrate,n_pre_exchange;
+  int n_pre_neighbor,n_post_neighbor;
   int n_pre_force,n_pre_reverse,n_post_force;
   int n_final_integrate,n_end_of_step,n_thermo_energy,n_thermo_energy_atom;
   int n_initial_integrate_respa,n_post_integrate_respa;
   int n_pre_force_respa,n_post_force_respa,n_final_integrate_respa;
-  int n_min_pre_exchange,n_min_pre_neighbor;
+  int n_min_pre_exchange,n_min_pre_neighbor,n_min_post_neighbor;
   int n_min_pre_force,n_min_pre_reverse,n_min_post_force,n_min_energy;
 
   int restart_pbc_any;       // 1 if any fix sets restart_pbc
@@ -53,12 +54,14 @@ class Modify : protected Pointers {
   virtual void setup(int);
   virtual void setup_pre_exchange();
   virtual void setup_pre_neighbor();
+  virtual void setup_post_neighbor();
   virtual void setup_pre_force(int);
   virtual void setup_pre_reverse(int, int);
   virtual void initial_integrate(int);
   virtual void post_integrate();
   virtual void pre_exchange();
   virtual void pre_neighbor();
+  virtual void post_neighbor();
   virtual void pre_force(int);
   virtual void pre_reverse(int,int);
   virtual void post_force(int);
@@ -78,6 +81,7 @@ class Modify : protected Pointers {
 
   virtual void min_pre_exchange();
   virtual void min_pre_neighbor();
+  virtual void min_post_neighbor();
   virtual void min_pre_force(int);
   virtual void min_pre_reverse(int,int);
   virtual void min_post_force(int);
@@ -123,14 +127,14 @@ class Modify : protected Pointers {
   // lists of fixes to apply at different stages of timestep
 
   int *list_initial_integrate,*list_post_integrate;
-  int *list_pre_exchange,*list_pre_neighbor;
+  int *list_pre_exchange,*list_pre_neighbor,*list_post_neighbor;
   int *list_pre_force,*list_pre_reverse,*list_post_force;
   int *list_final_integrate,*list_end_of_step,*list_thermo_energy;
   int *list_thermo_energy_atom;
   int *list_initial_integrate_respa,*list_post_integrate_respa;
   int *list_pre_force_respa,*list_post_force_respa;
   int *list_final_integrate_respa;
-  int *list_min_pre_exchange,*list_min_pre_neighbor;
+  int *list_min_pre_exchange,*list_min_pre_neighbor,*list_min_post_neighbor;
   int *list_min_pre_force,*list_min_pre_reverse,*list_min_post_force;
   int *list_min_energy;
 

src/neigh_list.cpp

@@ -40,16 +40,15 @@ NeighList::NeighList(LAMMPS *lmp) : Pointers(lmp)
   ilist = NULL;
   numneigh = NULL;
   firstneigh = NULL;
-  firstdouble = NULL;
 
   // defaults, but may be reset by post_constructor()
 
   occasional = 0;
   ghost = 0;
   ssa = 0;
+  history = 0;
   copy = 0;
   copymode = 0;
-  dnum = 0;
 
   // ptrs
 
@@ -60,9 +59,6 @@ NeighList::NeighList(LAMMPS *lmp) : Pointers(lmp)
   listskip = NULL;
   listfull = NULL;
 
-  listhistory = NULL;
-  fix_history = NULL;
-
   respamiddle = 0;
   listinner = NULL;
   listmiddle = NULL;
@@ -70,7 +66,6 @@ NeighList::NeighList(LAMMPS *lmp) : Pointers(lmp)
   fix_bond = NULL;
 
   ipage = NULL;
-  dpage = NULL;
 
   // Kokkos package
 
@@ -92,10 +87,7 @@ NeighList::~NeighList()
     memory->destroy(ilist);
     memory->destroy(numneigh);
     memory->sfree(firstneigh);
-    memory->sfree(firstdouble);
-
     delete [] ipage;
-    delete [] dpage;
   }
 
   delete [] iskip;
@@ -108,7 +100,6 @@ NeighList::~NeighList()
    copy -> set listcopy for list to copy from
    skip -> set listskip for list to skip from, create copy of itype,ijtype
    halffull -> set listfull for full list to derive from
-   history -> set LH and FH ptrs in partner list that uses the history info
    respaouter -> set listinner/listmiddle for other rRESPA lists
    bond -> set fix_bond to Fix that made the request
 ------------------------------------------------------------------------- */
@@ -120,8 +111,8 @@ void NeighList::post_constructor(NeighRequest *nq)
   occasional = nq->occasional;
   ghost = nq->ghost;
   ssa = nq->ssa;
+  history = nq->history;
   copy = nq->copy;
-  dnum = nq->dnum;
 
   if (nq->copy)
     listcopy = neighbor->lists[nq->copylist];
@@ -141,13 +132,6 @@ void NeighList::post_constructor(NeighRequest *nq)
   if (nq->halffull)
     listfull = neighbor->lists[nq->halffulllist];
 
-  if (nq->history) {
-    neighbor->lists[nq->historylist]->listhistory = this;
-    int tmp;
-    neighbor->lists[nq->historylist]->fix_history = 
-      (Fix *) ((Pair *) nq->requestor)->extract("history",tmp);
-  }
-  
   if (nq->respaouter) {
     if (nq->respamiddlelist < 0) {
       respamiddle = 0;
@@ -173,12 +157,6 @@ void NeighList::setup_pages(int pgsize_caller, int oneatom_caller)
   ipage = new MyPage<int>[nmypage];
   for (int i = 0; i < nmypage; i++)
     ipage[i].init(oneatom,pgsize,PGDELTA);
-
-  if (dnum) {
-    dpage = new MyPage<double>[nmypage];
-    for (int i = 0; i < nmypage; i++)
-      dpage[i].init(dnum*oneatom,dnum*pgsize,PGDELTA);
-  } else dpage = NULL;
 }
 
 /* ----------------------------------------------------------------------
@@ -195,11 +173,11 @@ void NeighList::grow(int nlocal, int nall)
 {
   // trigger grow() in children before possible return
 
-  if (listhistory) listhistory->grow(nlocal,nall);
   if (listinner) listinner->grow(nlocal,nall);
   if (listmiddle) listmiddle->grow(nlocal,nall);
 
   // skip if data structs are already big enough
+
   if (ssa) {
     if ((nlocal * 3) + nall <= maxatom) return;
   } else if (ghost) {
@@ -218,11 +196,6 @@ void NeighList::grow(int nlocal, int nall)
   memory->create(numneigh,maxatom,"neighlist:numneigh");
   firstneigh = (int **) memory->smalloc(maxatom*sizeof(int *),
                                         "neighlist:firstneigh");
-  if (dnum) {
-    memory->sfree(firstdouble);
-    firstdouble = (double **) memory->smalloc(maxatom*sizeof(double *),
-                                              "neighlist:firstdouble");
-  }
 }
 
 /* ----------------------------------------------------------------------
@@ -262,13 +235,11 @@ void NeighList::print_attributes()
   printf("  %d = kokkos host\n",rq->kokkos_host);
   printf("  %d = kokkos device\n",rq->kokkos_device);
   printf("  %d = ssa flag\n",ssa);
-  printf("  %d = dnum\n",dnum);
   printf("\n");
   printf("  %d = skip flag\n",rq->skip);
   printf("  %d = off2on\n",rq->off2on);
   printf("  %d = copy flag\n",rq->copy);
   printf("  %d = half/full\n",rq->halffull);
-  printf("  %d = history/partner\n",rq->history_partner);
   printf("\n");
 }
 
@@ -292,12 +263,5 @@ bigint NeighList::memory_usage()
       bytes += ipage[i].size();
   }
 
-  if (dnum && dpage) {
-    for (int i = 0; i < nmypage; i++) {
-      bytes += maxatom * sizeof(double *);
-      bytes += dpage[i].size();
-    }
-  }
-
   return bytes;
 }

src/neigh_list.h

@@ -34,9 +34,9 @@ class NeighList : protected Pointers {
   int occasional;                  // 0 if build every reneighbor, 1 if not
   int ghost;                       // 1 if list stores neighbors of ghosts
   int ssa;                         // 1 if list stores Shardlow data
-  int copy;                        // 1 if this list is (host) copied from another list
+  int history;                     // 1 if there is neigh history (FixNeighHist)
+  int copy;                        // 1 if this list is copied from another list
   int copymode;                    // 1 if this is a Kokkos on-device copy
-  int dnum;                        // # of doubles per neighbor, 0 if none
 
   // data structs to store neighbor pairs I,J and associated values
 
@@ -45,13 +45,11 @@ class NeighList : protected Pointers {
   int *ilist;                      // local indices of I atoms
   int *numneigh;                   // # of J neighbors for each I atom
   int **firstneigh;                // ptr to 1st J int value of each I atom
-  double **firstdouble;            // ptr to 1st J double value of each I atom
   int maxatom;                     // size of allocated per-atom arrays
 
   int pgsize;                      // size of each page
   int oneatom;                     // max size for one atom
   MyPage<int> *ipage;              // pages of neighbor indices
-  MyPage<double> *dpage;           // pages of neighbor doubles, if dnum > 0
 
   // atom types to skip when building list
   // copied info from corresponding request into realloced vec/array
@@ -65,9 +63,6 @@ class NeighList : protected Pointers {
   NeighList *listskip;          // me = skip list, point to list I skip from
   NeighList *listfull;          // me = half list, point to full I derive from
 
-  NeighList *listhistory;       // list storing neigh history
-  class Fix *fix_history;       // fix that stores history info
-
   int respamiddle;              // 1 if this respaouter has middle list
   NeighList *listinner;         // me = respaouter, point to respainner
   NeighList *listmiddle;        // me = respaouter, point to respamiddle

src/neigh_request.cpp

@@ -42,7 +42,7 @@ NeighRequest::NeighRequest(LAMMPS *lmp) : Pointers(lmp)
   // default is use newton_pair setting in force
   // default is no neighbors of ghosts
   // default is use cutoffs, not size of particles
-  // default is no additional neighbor history info
+  // default is no associated neighbor history info in FixNeighHistory
   // default is no one-sided sphere/surface interactions (when size = 1)
   // default is neighbors of atoms, not bonds
   // default is no multilevel rRESPA neighbors
@@ -68,8 +68,6 @@ NeighRequest::NeighRequest(LAMMPS *lmp) : Pointers(lmp)
   cut = 0;
   cutoff = 0.0;
 
-  dnum = 0;
-
   // skip info, default is no skipping
   
   skip = 0;
@@ -88,8 +86,6 @@ NeighRequest::NeighRequest(LAMMPS *lmp) : Pointers(lmp)
   copylist = -1;
   halffull = 0;
   halffulllist = -1;
-  history_partner = 0;
-  historylist = -1;
   respaouterlist = -1;
   respamiddlelist = -1;
   respainnerlist = -1;
@@ -158,8 +154,6 @@ int NeighRequest::identical(NeighRequest *other)
   if (copy != other->copy) same = 0;
   if (cutoff != other->cutoff) same = 0;
 
-  if (dnum != other->dnum) same = 0;
-
   if (skip != other->skip) same = 0;
   if (skip) same = same_skip(other);
 
@@ -226,8 +220,6 @@ void NeighRequest::copy_request(NeighRequest *other, int skipflag)
   cut = other->cut;
   cutoff = other->cutoff;
 
-  dnum = other->dnum;
-
   iskip = NULL;
   ijskip = NULL;
 

src/neigh_request.h

@@ -59,7 +59,7 @@ class NeighRequest : protected Pointers {
   
   int ghost;             // 1 if includes ghost atom neighbors
   int size;              // 1 if pair cutoff set by particle radius
-  int history;           // 1 if stores neighbor history info
+  int history;           // 1 if there is also neigh history info (FixNeighHist)
   int granonesided;      // 1 if one-sided granular list for 
                          //   sphere/surf interactions
   int respainner;        // 1 if a rRESPA inner list
@@ -74,8 +74,6 @@ class NeighRequest : protected Pointers {
   int cut;               // 1 if use a non-standard cutoff length
   double cutoff;         // special cutoff distance for this list
 
-  int dnum;              // # of extra floating point values stored in list
-
   // flags set by pair hybrid
 
   int skip;              // 1 if this list skips atom types from another list
@@ -100,10 +98,6 @@ class NeighRequest : protected Pointers {
   int halffull;          // 1 if half list computed from another full list
   int halffulllist;      // index of full list to derive half from
 
-  int history_partner;   // 1 if this list partners with a history list
-  int historylist;       // index of the associated history list
-                         // for history = 1, index of the non-history partner
-
   int respaouterlist;    // index of respaouter/middle/inner lists
   int respamiddlelist;   // which this rREPSA list is associated with
   int respainnerlist;    // each rRESPA style list points at the others

src/neighbor.cpp

@@ -133,8 +133,6 @@ pairclass(NULL), pairnames(NULL), pairmasks(NULL)
   old_pgsize = pgsize;
   old_oneatom = oneatom;
 
-  zeroes = NULL;
-
   binclass = NULL;
   binnames = NULL;
   binmasks = NULL;
@@ -208,8 +206,6 @@ Neighbor::~Neighbor()
     if (old_requests[i]) delete old_requests[i];
   memory->sfree(old_requests);
 
-  delete [] zeroes;
-
   delete [] binclass;
   delete [] binnames;
   delete [] binmasks;
@@ -666,14 +662,12 @@ int Neighbor::init_pair()
   // purpose is to avoid duplicate or inefficient builds
   // may add new requests if a needed request to derive from does not exist
   // methods:
-  //   (1) other = point history and rRESPA lists at their partner lists
+  //   (1) other = point rRESPA lists at their partner lists
   //   (2) skip = create any new non-skip lists needed by pair hybrid skip lists
   //   (3) granular = adjust parent and skip lists for granular onesided usage
   //   (4) h/f = pair up any matching half/full lists
   //   (5) copy = convert as many lists as possible to copy lists
   // order of morph methods matters:
-  //   (1) before (2), b/c (2) needs to know history partner pairings
-  //   (2) after (1), b/c (2) may also need to create new history lists
   //   (3) after (2), b/c it adjusts lists created by (2)
   //   (4) after (2) and (3), 
   //       b/c (2) may create new full lists, (3) may change them
@@ -827,23 +821,14 @@ int Neighbor::init_pair()
   }
 
   // allocate initial pages for each list, except if copy flag set
-  // allocate dnum vector of zeroes if set
   
-  int dnummax = 0;
   for (i = 0; i < nlist; i++) {
     if (lists[i]->copy) continue;
     lists[i]->setup_pages(pgsize,oneatom);
-    dnummax = MAX(dnummax,lists[i]->dnum);
-  }
-  
-  if (dnummax) {
-    delete [] zeroes;
-    zeroes = new double[dnummax];
-    for (i = 0; i < dnummax; i++) zeroes[i] = 0.0;
   }
 
   // first-time allocation of per-atom data for lists that are built and store
-  // lists that are not built: granhistory, respa inner/middle (no neigh_pair)
+  // lists that are not built: respa inner/middle (no neigh_pair)
   // lists that do not store: copy 
   // use atom->nmax for both grow() args
   //   i.e. grow first time to expanded size to avoid future reallocs
@@ -923,7 +908,7 @@ int Neighbor::init_pair()
 
 /* ----------------------------------------------------------------------
    scan NeighRequests to set additional flags
-   only for history, respaouter, custom cutoff lists
+   only for respaouter, custom cutoff lists
 ------------------------------------------------------------------------- */
 
 void Neighbor::morph_other()
@@ -933,14 +918,6 @@ void Neighbor::morph_other()
   for (int i = 0; i < nrequest; i++) {
     irq = requests[i];
 
-    // if history, point this list and partner list at each other
-
-    if (irq->history) {
-      irq->historylist = i-1;
-      requests[i-1]->history_partner = 1;
-      requests[i-1]->historylist = i;
-    }
-
     // if respaouter, point all associated rRESPA lists at each other
 
     if (irq->respaouter) {
@@ -987,7 +964,6 @@ void Neighbor::morph_skip()
     // halffull list and its full parent may both skip,
     //   but are checked to insure matching skip info
 
-    if (irq->history) continue;
     if (irq->respainner || irq->respamiddle) continue;
     if (irq->halffull) continue;
     if (irq->copy) continue;
@@ -1022,11 +998,11 @@ void Neighbor::morph_skip()
       //   or their data structures are different
       // this includes custom cutoff set by requestor
       // no need to check respaouter b/c it stores same pairs
-      // no need to check dnum b/c only set for history
       // NOTE: need check for 2 Kokkos flags?
 
       if (irq->ghost != jrq->ghost) continue;
       if (irq->size != jrq->size) continue;
+      if (irq->history != jrq->history) continue;
       if (irq->bond != jrq->bond) continue;
       if (irq->omp != jrq->omp) continue;
       if (irq->intel != jrq->intel) continue;
@@ -1045,8 +1021,8 @@ void Neighbor::morph_skip()
     // else create a new identical list except non-skip
     // for new list, set neigh = 1, skip = 0, no skip vec/array,
     //   copy unique flag (since copy_request() will not do it)
-    // note: parents of skip lists do not have associated history list
+    // note: parents of skip lists do not have associated history
-    //   b/c child skip lists store their own history info
+    //   b/c child skip lists have the associated history
 
     if (j < nrequest) irq->skiplist = j;
     else {
@@ -1107,7 +1083,6 @@ void Neighbor::morph_granular()
       if (onesided == 2) break;
     }
 
-
     // if onesided = 2, parent has children with both granonesided = 0/1
     // force parent newton off (newton = 2) to enable onesided skip by child
     // set parent granonesided = 0, so it stores all neighs in usual manner
@@ -1159,7 +1134,6 @@ void Neighbor::morph_halffull()
     // these lists are created other ways, no need for halffull
     // do want to process skip lists
 
-    if (irq->history) continue;
     if (irq->respainner || irq->respamiddle) continue;
     if (irq->copy) continue;
 
@@ -1180,10 +1154,10 @@ void Neighbor::morph_halffull()
       //   or their data structures are different
       // this includes custom cutoff set by requestor
       // no need to check respaouter b/c it stores same pairs
-      // no need to check dnum b/c only set for history
 
       if (irq->ghost != jrq->ghost) continue;
       if (irq->size != jrq->size) continue;
+      if (irq->history != jrq->history) continue;
       if (irq->bond != jrq->bond) continue;
       if (irq->omp != jrq->omp) continue;
       if (irq->intel != jrq->intel) continue;
@@ -1233,7 +1207,6 @@ void Neighbor::morph_copy()
     // these lists are created other ways, no need to copy
     // skip lists are eligible to become a copy list
 
-    if (irq->history) continue;
     if (irq->respainner || irq->respamiddle) continue;
     
     // check all other lists
@@ -1272,9 +1245,8 @@ void Neighbor::morph_copy()
 
       if (irq->ghost && !jrq->ghost) continue;
 
-      // do not copy from a history list or a respa middle/inner list
+      // do not copy from a respa middle/inner list
 
-      if (jrq->history) continue;
       if (jrq->respamiddle) continue;
       if (jrq->respainner) continue;
 
@@ -1284,10 +1256,10 @@ void Neighbor::morph_copy()
       // this includes custom cutoff set by requestor
       // no need to check respaouter b/c it stores same pairs
       // no need to check omp b/c it stores same pairs
-      // no need to check dnum b/c only set for history
       // NOTE: need check for 2 Kokkos flags?
 
       if (irq->size != jrq->size) continue;
+      if (irq->history != jrq->history) continue;
       if (irq->bond != jrq->bond) continue;
       if (irq->intel != jrq->intel) continue;
       if (irq->kokkos_host != jrq->kokkos_host) continue;
@@ -1535,9 +1507,7 @@ void Neighbor::print_pairwise_info()
 
         // order these to get single output of most relevant
 
-        if (rq->history) 
+        if (rq->copy)
-          fprintf(out,", history for (%d)",rq->historylist+1);
-        else if (rq->copy)
           fprintf(out,", copy from (%d)",rq->copylist+1);
         else if (rq->halffull) 
           fprintf(out,", half/full from (%d)",rq->halffulllist+1);
@@ -1659,7 +1629,6 @@ int Neighbor::choose_bin(NeighRequest *rq)
 
   if (style == NSQ) return 0;
   if (rq->skip || rq->copy || rq->halffull) return 0;
-  if (rq->history) return 0;
   if (rq->respainner || rq->respamiddle) return 0;
 
   // use request settings to match exactly one NBin class mask
@@ -1701,7 +1670,6 @@ int Neighbor::choose_stencil(NeighRequest *rq)
 
   if (style == NSQ) return 0;
   if (rq->skip || rq->copy || rq->halffull) return 0;
-  if (rq->history) return 0;
   if (rq->respainner || rq->respamiddle) return 0;
 
   // convert newton request to newtflag = on or off
@@ -1795,7 +1763,6 @@ int Neighbor::choose_pair(NeighRequest *rq)
 {
   // no neighbor list build performed
 
-  if (rq->history) return 0;
   if (rq->respainner || rq->respamiddle) return 0;
 
   // error check for includegroup with ghost neighbor request

src/neighbor.h

@@ -54,7 +54,6 @@ class Neighbor : protected Pointers {
 
   double *bboxlo,*bboxhi;          // ptrs to full domain bounding box
                                    // different for orthog vs triclinic
-  double *zeroes;                  // vector of zeroes for shear history init
 
   // exclusion info, used by NeighPair
 

src/npair.cpp

@@ -66,7 +66,6 @@ void NPair::copy_neighbor_info()
   cut_inner_sq = neighbor->cut_inner_sq;
   cut_middle_sq = neighbor->cut_middle_sq;
   cut_middle_inside_sq = neighbor->cut_middle_inside_sq;
-  zeroes = neighbor->zeroes;
   bboxlo = neighbor->bboxlo;
   bboxhi = neighbor->bboxhi;
 

src/npair.h

@@ -47,7 +47,6 @@ class NPair : protected Pointers {
   double cut_inner_sq;
   double cut_middle_sq;
   double cut_middle_inside_sq;
-  double *zeroes;
   double *bboxlo,*bboxhi;
 
   // exclusion data from Neighbor class

src/npair_copy.cpp

@@ -40,7 +40,5 @@ void NPairCopy::build(NeighList *list)
   list->ilist = listcopy->ilist;
   list->numneigh = listcopy->numneigh;
   list->firstneigh = listcopy->firstneigh;
-  list->firstdouble = listcopy->firstdouble;
   list->ipage = listcopy->ipage;
-  list->dpage = listcopy->dpage;
 }

src/npair_half_size_bin_newtoff.cpp

@@ -17,9 +17,6 @@
 #include "neigh_list.h"
 #include "atom.h"
 #include "atom_vec.h"
-#include "molecule.h"
-#include "domain.h"
-#include "fix_shear_history.h"
 #include "my_page.h"
 #include "error.h"
 
@@ -32,7 +29,6 @@ NPairHalfSizeBinNewtoff::NPairHalfSizeBinNewtoff(LAMMPS *lmp) : NPair(lmp) {}
 /* ----------------------------------------------------------------------
    size particles
    binned neighbor list construction with partial Newton's 3rd law
-   shear history must be accounted for when a neighbor pair is added
    each owned atom i checks own bin and surrounding bins in non-Newton stencil
    pair stored once if i,j are both owned and i < j
    pair stored by me if j is ghost (also stored by proc owning j)
@@ -40,20 +36,10 @@ NPairHalfSizeBinNewtoff::NPairHalfSizeBinNewtoff(LAMMPS *lmp) : NPair(lmp) {}
 
 void NPairHalfSizeBinNewtoff::build(NeighList *list)
 {
-  int i,j,k,m,n,nn,ibin,dnum,dnumbytes;
+  int i,j,k,m,n,nn,ibin;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   double radi,radsum,cutsq;
-  int *neighptr,*touchptr;
+  int *neighptr;
-  double *shearptr;
-
-  int *npartner;
-  tagint **partner;
-  double **shearpartner;
-  int **firsttouch;
-  double **firstshear;
-  MyPage<int> *ipage_touch;
-  MyPage<double> *dpage_shear;
-  NeighList *listhistory;
 
   double **x = atom->x;
   double *radius = atom->radius;
@@ -64,42 +50,20 @@ void NPairHalfSizeBinNewtoff::build(NeighList *list)
   int nlocal = atom->nlocal;
   if (includegroup) nlocal = atom->nfirst;
 
+  int history = list->history;
   int *ilist = list->ilist;
   int *numneigh = list->numneigh;
   int **firstneigh = list->firstneigh;
   MyPage<int> *ipage = list->ipage;
 
-  FixShearHistory *fix_history = (FixShearHistory *) list->fix_history;
+  int mask_history = 3 << SBBITS;
-  if (fix_history) {
-    fix_history->nlocal_neigh = nlocal;
-    fix_history->nall_neigh = nlocal + atom->nghost;
-    npartner = fix_history->npartner;
-    partner = fix_history->partner;
-    shearpartner = fix_history->shearpartner;
-    listhistory = list->listhistory;
-    firsttouch = listhistory->firstneigh;
-    firstshear = listhistory->firstdouble;
-    ipage_touch = listhistory->ipage;
-    dpage_shear = listhistory->dpage;
-    dnum = listhistory->dnum;
-    dnumbytes = dnum * sizeof(double);
-  }
 
   int inum = 0;
   ipage->reset();
-  if (fix_history) {
-    ipage_touch->reset();
-    dpage_shear->reset();
-  }
 
   for (i = 0; i < nlocal; i++) {
     n = 0;
     neighptr = ipage->vget();
-    if (fix_history) {
-      nn = 0;
-      touchptr = ipage_touch->vget();
-      shearptr = dpage_shear->vget();
-    }
 
     xtmp = x[i][0];
     ytmp = x[i][1];
@@ -125,29 +89,10 @@ void NPairHalfSizeBinNewtoff::build(NeighList *list)
         cutsq = (radsum+skin) * (radsum+skin);
         
         if (rsq <= cutsq) {
-          neighptr[n] = j;
+          if (history && rsq < radsum*radsum)
-
+            neighptr[n++] = j ^ mask_history;
-          if (fix_history) {
+          else
-            if (rsq < radsum*radsum) {
+            neighptr[n++] = j;
-              for (m = 0; m < npartner[i]; m++)
-                if (partner[i][m] == tag[j]) break;
-              if (m < npartner[i]) {
-                touchptr[n] = 1;
-                memcpy(&shearptr[nn],&shearpartner[i][dnum*m],dnumbytes);
-                nn += dnum;
-              } else {
-                touchptr[n] = 0;
-                memcpy(&shearptr[nn],zeroes,dnumbytes);
-                nn += dnum;
-              }
-            } else {
-              touchptr[n] = 0;
-              memcpy(&shearptr[nn],zeroes,dnumbytes);
-              nn += dnum;
-            }
-          }
-
-          n++;
         }
       }
     }
@@ -158,13 +103,6 @@ void NPairHalfSizeBinNewtoff::build(NeighList *list)
     ipage->vgot(n);
     if (ipage->status())
       error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
-
-    if (fix_history) {
-      firsttouch[i] = touchptr;
-      firstshear[i] = shearptr;
-      ipage_touch->vgot(n);
-      dpage_shear->vgot(nn);
-    }
   }
 
   list->inum = inum;

src/npair_half_size_bin_newton.cpp

@@ -17,9 +17,6 @@
 #include "neigh_list.h"
 #include "atom.h"
 #include "atom_vec.h"
-#include "molecule.h"
-#include "domain.h"
-#include "fix_shear_history.h"
 #include "my_page.h"
 #include "error.h"
 
@@ -32,27 +29,16 @@ NPairHalfSizeBinNewton::NPairHalfSizeBinNewton(LAMMPS *lmp) : NPair(lmp) {}
 /* ----------------------------------------------------------------------
    size particles
    binned neighbor list construction with full Newton's 3rd law
-   shear history must be accounted for when a neighbor pair is added
    each owned atom i checks its own bin and other bins in Newton stencil
    every pair stored exactly once by some processor
 ------------------------------------------------------------------------- */
 
 void NPairHalfSizeBinNewton::build(NeighList *list)
 {
-  int i,j,k,m,n,nn,ibin,dnum,dnumbytes;
+  int i,j,k,m,n,nn,ibin;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   double radi,radsum,cutsq;
-  int *neighptr,*touchptr;
+  int *neighptr;
-  double *shearptr;
-
-  int *npartner;
-  tagint **partner;
-  double **shearpartner;
-  int **firsttouch;
-  double **firstshear;
-  MyPage<int> *ipage_touch;
-  MyPage<double> *dpage_shear;
-  NeighList *listhistory;
 
   double **x = atom->x;
   double *radius = atom->radius;
@@ -63,42 +49,20 @@ void NPairHalfSizeBinNewton::build(NeighList *list)
   int nlocal = atom->nlocal;
   if (includegroup) nlocal = atom->nfirst;
 
+  int history = list->history;
   int *ilist = list->ilist;
   int *numneigh = list->numneigh;
   int **firstneigh = list->firstneigh;
   MyPage<int> *ipage = list->ipage;
 
-  FixShearHistory *fix_history = (FixShearHistory *) list->fix_history;
+  int mask_history = 3 << SBBITS;
-  if (fix_history) {
-    fix_history->nlocal_neigh = nlocal;
-    fix_history->nall_neigh = nlocal + atom->nghost;
-    npartner = fix_history->npartner;
-    partner = fix_history->partner;
-    shearpartner = fix_history->shearpartner;
-    listhistory = list->listhistory;
-    firsttouch = listhistory->firstneigh;
-    firstshear = listhistory->firstdouble;
-    ipage_touch = listhistory->ipage;
-    dpage_shear = listhistory->dpage;
-    dnum = listhistory->dnum;
-    dnumbytes = dnum * sizeof(double);
-  }
 
   int inum = 0;
   ipage->reset();
-  if (fix_history) {
-    ipage_touch->reset();
-    dpage_shear->reset();
-  }
 
   for (i = 0; i < nlocal; i++) {
     n = 0;
     neighptr = ipage->vget();
-    if (fix_history) {
-      nn = 0;
-      touchptr = ipage_touch->vget();
-      shearptr = dpage_shear->vget();
-    }
 
     xtmp = x[i][0];
     ytmp = x[i][1];
@@ -128,29 +92,10 @@ void NPairHalfSizeBinNewton::build(NeighList *list)
       cutsq = (radsum+skin) * (radsum+skin);
 
       if (rsq <= cutsq) {
-        neighptr[n] = j;
+        if (history && rsq < radsum*radsum) 
-
+          neighptr[n++] = j ^ mask_history;
-        if (fix_history) {
+        else 
-          if (rsq < radsum*radsum) {
+          neighptr[n++] = j;
-            for (m = 0; m < npartner[i]; m++)
-              if (partner[i][m] == tag[j]) break;
-            if (m < npartner[i]) {
-              touchptr[n] = 1;
-              memcpy(&shearptr[nn],&shearpartner[i][dnum*m],dnumbytes);
-              nn += dnum;
-            } else {
-              touchptr[n] = 0;
-              memcpy(&shearptr[nn],zeroes,dnumbytes);
-              nn += dnum;
-            }
-          } else {
-            touchptr[n] = 0;
-            memcpy(&shearptr[nn],zeroes,dnumbytes);
-            nn += dnum;
-          }
-        }
-
-        n++;
       }
     }
 
@@ -169,29 +114,10 @@ void NPairHalfSizeBinNewton::build(NeighList *list)
         cutsq = (radsum+skin) * (radsum+skin);
 
         if (rsq <= cutsq) {
-          neighptr[n] = j;
+          if (history && rsq < radsum*radsum) 
-
+            neighptr[n++] = j ^ mask_history;
-          if (fix_history) {
+          else
-            if (rsq < radsum*radsum) {
+            neighptr[n++] = j;
-              for (m = 0; m < npartner[i]; m++)
-                if (partner[i][m] == tag[j]) break;
-              if (m < npartner[i]) {
-                touchptr[n] = 1;
-                memcpy(&shearptr[nn],&shearpartner[i][dnum*m],dnumbytes);
-                nn += dnum;
-              } else {
-                touchptr[n] = 0;
-                memcpy(&shearptr[nn],zeroes,dnumbytes);
-                nn += dnum;
-              }
-            } else {
-              touchptr[n] = 0;
-              memcpy(&shearptr[nn],zeroes,dnumbytes);
-              nn += dnum;
-            }
-          }
-
-          n++;
         }
       }
     }
@@ -202,13 +128,6 @@ void NPairHalfSizeBinNewton::build(NeighList *list)
     ipage->vgot(n);
     if (ipage->status())
       error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
-
-    if (fix_history) {
-      firsttouch[i] = touchptr;
-      firstshear[i] = shearptr;
-      ipage_touch->vgot(n);
-      dpage_shear->vgot(nn);
-    }
   }
 
   list->inum = inum;

src/npair_half_size_bin_newton_tri.cpp

@@ -17,9 +17,6 @@
 #include "neigh_list.h"
 #include "atom.h"
 #include "atom_vec.h"
-#include "molecule.h"
-#include "domain.h"
-#include "fix_shear_history.h"
 #include "my_page.h"
 #include "error.h"
 
@@ -33,27 +30,16 @@ NPairHalfSizeBinNewtonTri::NPairHalfSizeBinNewtonTri(LAMMPS *lmp) :
 /* ----------------------------------------------------------------------
    size particles
    binned neighbor list construction with Newton's 3rd law for triclinic
-   shear history must be accounted for when a neighbor pair is added
    each owned atom i checks its own bin and other bins in triclinic stencil
    every pair stored exactly once by some processor
 ------------------------------------------------------------------------- */
 
 void NPairHalfSizeBinNewtonTri::build(NeighList *list)
 {
-  int i,j,k,m,n,nn,ibin,dnum,dnumbytes;
+  int i,j,k,m,n,nn,ibin;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   double radi,radsum,cutsq;
-  int *neighptr,*touchptr;
+  int *neighptr;
-  double *shearptr;
-
-  int *npartner;
-  tagint **partner;
-  double **shearpartner;
-  int **firsttouch;
-  double **firstshear;
-  MyPage<int> *ipage_touch;
-  MyPage<double> *dpage_shear;
-  NeighList *listhistory;
 
   double **x = atom->x;
   double *radius = atom->radius;
@@ -64,42 +50,20 @@ void NPairHalfSizeBinNewtonTri::build(NeighList *list)
   int nlocal = atom->nlocal;
   if (includegroup) nlocal = atom->nfirst;
 
+  int history = list->history;
   int *ilist = list->ilist;
   int *numneigh = list->numneigh;
   int **firstneigh = list->firstneigh;
   MyPage<int> *ipage = list->ipage;
 
-  FixShearHistory *fix_history = (FixShearHistory *) list->fix_history;
+  int mask_history = 3 << SBBITS;
-  if (fix_history) {
-    fix_history->nlocal_neigh = nlocal;
-    fix_history->nall_neigh = nlocal + atom->nghost;
-    npartner = fix_history->npartner;
-    partner = fix_history->partner;
-    shearpartner = fix_history->shearpartner;
-    listhistory = list->listhistory;
-    firsttouch = listhistory->firstneigh;
-    firstshear = listhistory->firstdouble;
-    ipage_touch = listhistory->ipage;
-    dpage_shear = listhistory->dpage;
-    dnum = listhistory->dnum;
-    dnumbytes = dnum * sizeof(double);
-  }
 
   int inum = 0;
   ipage->reset();
-  if (fix_history) {
-    ipage_touch->reset();
-    dpage_shear->reset();
-  }
 
   for (i = 0; i < nlocal; i++) {
     n = 0;
     neighptr = ipage->vget();
-    if (fix_history) {
-      nn = 0;
-      touchptr = ipage_touch->vget();
-      shearptr = dpage_shear->vget();
-    }
 
     xtmp = x[i][0];
     ytmp = x[i][1];
@@ -134,29 +98,10 @@ void NPairHalfSizeBinNewtonTri::build(NeighList *list)
         cutsq = (radsum+skin) * (radsum+skin);
 
         if (rsq <= cutsq) {
+          if (history && rsq < radsum*radsum)
+            neighptr[n++] = j ^ mask_history;
+          else
             neighptr[n++] = j;
-
-          if (fix_history) {
-            if (rsq < radsum*radsum) {
-              for (m = 0; m < npartner[i]; m++)
-                if (partner[i][m] == tag[j]) break;
-              if (m < npartner[i]) {
-                touchptr[n] = 1;
-                memcpy(&shearptr[nn],&shearpartner[i][dnum*m],dnumbytes);
-                nn += dnum;
-              } else {
-                touchptr[n] = 0;
-                memcpy(&shearptr[nn],zeroes,dnumbytes);
-                nn += dnum;
-              }
-            } else {
-              touchptr[n] = 0;
-              memcpy(&shearptr[nn],zeroes,dnumbytes);
-              nn += dnum;
-            }
-          }
-
-          n++;
         }
       }
     }
@@ -167,13 +112,6 @@ void NPairHalfSizeBinNewtonTri::build(NeighList *list)
     ipage->vgot(n);
     if (ipage->status())
       error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
-
-    if (fix_history) {
-      firsttouch[i] = touchptr;
-      firstshear[i] = shearptr;
-      ipage_touch->vgot(n);
-      dpage_shear->vgot(nn);
-    }
   }
 
   list->inum = inum;

src/npair_half_size_nsq_newtoff.cpp

@@ -18,9 +18,6 @@
 #include "atom.h"
 #include "atom_vec.h"
 #include "group.h"
-#include "molecule.h"
-#include "domain.h"
-#include "fix_shear_history.h"
 #include "my_page.h"
 #include "error.h"
 
@@ -33,27 +30,16 @@ NPairHalfSizeNsqNewtoff::NPairHalfSizeNsqNewtoff(LAMMPS *lmp) : NPair(lmp) {}
 /* ----------------------------------------------------------------------
    size particles
    N^2 / 2 search for neighbor pairs with partial Newton's 3rd law
-   shear history must be accounted for when a neighbor pair is added
    pair added to list if atoms i and j are both owned and i < j
    pair added if j is ghost (also stored by proc owning j)
 ------------------------------------------------------------------------- */
 
 void NPairHalfSizeNsqNewtoff::build(NeighList *list)
 {
-  int i,j,m,n,nn,bitmask,dnum,dnumbytes;
+  int i,j,m,n,nn,bitmask;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   double radi,radsum,cutsq;
-  int *neighptr,*touchptr;
+  int *neighptr;
-  double *shearptr;
-
-  int *npartner;
-  tagint **partner;
-  double **shearpartner;
-  int **firsttouch;
-  double **firstshear;
-  MyPage<int> *ipage_touch;
-  MyPage<double> *dpage_shear;
-  NeighList *listhistory;
 
   double **x = atom->x;
   double *radius = atom->radius;
@@ -68,42 +54,20 @@ void NPairHalfSizeNsqNewtoff::build(NeighList *list)
     bitmask = group->bitmask[includegroup];
   }
 
+  int history = list->history;
   int *ilist = list->ilist;
   int *numneigh = list->numneigh;
   int **firstneigh = list->firstneigh;
   MyPage<int> *ipage = list->ipage;
 
-  FixShearHistory *fix_history = (FixShearHistory *) list->fix_history;
+  int mask_history = 3 << SBBITS;
-  if (fix_history) {
-    fix_history->nlocal_neigh = nlocal;
-    fix_history->nall_neigh = nall;
-    npartner = fix_history->npartner;
-    partner = fix_history->partner;
-    shearpartner = fix_history->shearpartner;
-    listhistory = list->listhistory;
-    firsttouch = listhistory->firstneigh;
-    firstshear = listhistory->firstdouble;
-    ipage_touch = listhistory->ipage;
-    dpage_shear = listhistory->dpage;
-    dnum = listhistory->dnum;
-    dnumbytes = dnum * sizeof(double);
-  }
 
   int inum = 0;
   ipage->reset();
-  if (fix_history) {
-    ipage_touch->reset();
-    dpage_shear->reset();
-  }
 
   for (i = 0; i < nlocal; i++) {
     n = 0;
     neighptr = ipage->vget();
-    if (fix_history) {
-      nn = 0;
-      touchptr = ipage_touch->vget();
-      shearptr = dpage_shear->vget();
-    }
 
     xtmp = x[i][0];
     ytmp = x[i][1];
@@ -124,29 +88,10 @@ void NPairHalfSizeNsqNewtoff::build(NeighList *list)
       cutsq = (radsum+skin) * (radsum+skin);
 
       if (rsq <= cutsq) {
-        neighptr[n] = j;
+        if (history && rsq < radsum*radsum)
-
+          neighptr[n++] = j ^ mask_history;
-        if (fix_history) {
+        else
-          if (rsq < radsum*radsum) {
+          neighptr[n++] = j;
-            for (m = 0; m < npartner[i]; m++)
-              if (partner[i][m] == tag[j]) break;
-            if (m < npartner[i]) {
-              touchptr[n] = 1;
-              memcpy(&shearptr[nn],&shearpartner[i][dnum*m],dnumbytes);
-              nn += dnum;
-            } else {
-              touchptr[n] = 0;
-              memcpy(&shearptr[nn],zeroes,dnumbytes);
-              nn += dnum;
-            }
-          } else {
-            touchptr[n] = 0;
-            memcpy(&shearptr[nn],zeroes,dnumbytes);
-            nn += dnum;
-          }
-        }
-
-        n++;
       }
     }
 
@@ -156,13 +101,6 @@ void NPairHalfSizeNsqNewtoff::build(NeighList *list)
     ipage->vgot(n);
     if (ipage->status())
       error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
-
-    if (fix_history) {
-      firsttouch[i] = touchptr;
-      firstshear[i] = shearptr;
-      ipage_touch->vgot(n);
-      dpage_shear->vgot(nn);
-    }
   }
 
   list->inum = inum;

src/npair_half_size_nsq_newton.cpp

@@ -18,9 +18,6 @@
 #include "atom.h"
 #include "atom_vec.h"
 #include "group.h"
-#include "molecule.h"
-#include "domain.h"
-#include "fix_shear_history.h"
 #include "my_page.h"
 #include "error.h"
 
@@ -33,7 +30,6 @@ NPairHalfSizeNsqNewton::NPairHalfSizeNsqNewton(LAMMPS *lmp) : NPair(lmp) {}
 /* ----------------------------------------------------------------------
    size particles
    N^2 / 2 search for neighbor pairs with full Newton's 3rd law
-   shear history must be accounted for when a neighbor pair is added
    pair added to list if atoms i and j are both owned and i < j
    if j is ghost only me or other proc adds pair
    decision based on itag,jtag tests
@@ -41,20 +37,10 @@ NPairHalfSizeNsqNewton::NPairHalfSizeNsqNewton(LAMMPS *lmp) : NPair(lmp) {}
 
 void NPairHalfSizeNsqNewton::build(NeighList *list)
 {
-  int i,j,m,n,nn,itag,jtag,bitmask,dnum,dnumbytes;
+  int i,j,m,n,nn,itag,jtag,bitmask;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   double radi,radsum,cutsq;
-  int *neighptr,*touchptr;
+  int *neighptr;
-  double *shearptr;
-
-  int *npartner;
-  tagint **partner;
-  double **shearpartner;
-  int **firsttouch;
-  double **firstshear;
-  MyPage<int> *ipage_touch;
-  MyPage<double> *dpage_shear;
-  NeighList *listhistory;
 
   double **x = atom->x;
   double *radius = atom->radius;
@@ -69,42 +55,20 @@ void NPairHalfSizeNsqNewton::build(NeighList *list)
     bitmask = group->bitmask[includegroup];
   }
 
+  int history = list->history;
   int *ilist = list->ilist;
   int *numneigh = list->numneigh;
   int **firstneigh = list->firstneigh;
   MyPage<int> *ipage = list->ipage;
 
-  FixShearHistory *fix_history = (FixShearHistory *) list->fix_history;
+  int mask_history = 3 << SBBITS;
-  if (fix_history) {
-    fix_history->nlocal_neigh = nlocal;
-    fix_history->nall_neigh = nall;
-    npartner = fix_history->npartner;
-    partner = fix_history->partner;
-    shearpartner = fix_history->shearpartner;
-    listhistory = list->listhistory;
-    firsttouch = listhistory->firstneigh;
-    firstshear = listhistory->firstdouble;
-    ipage_touch = listhistory->ipage;
-    dpage_shear = listhistory->dpage;
-    dnum = listhistory->dnum;
-    dnumbytes = dnum * sizeof(double);
-  }
 
   int inum = 0;
   ipage->reset();
-  if (fix_history) {
-    ipage_touch->reset();
-    dpage_shear->reset();
-  }
 
   for (i = 0; i < nlocal; i++) {
     n = 0;
     neighptr = ipage->vget();
-    if (fix_history) {
-      nn = 0;
-      touchptr = ipage_touch->vget();
-      shearptr = dpage_shear->vget();
-    }
 
     itag = tag[i];
     xtmp = x[i][0];
@@ -142,29 +106,10 @@ void NPairHalfSizeNsqNewton::build(NeighList *list)
       cutsq = (radsum+skin) * (radsum+skin);
 
       if (rsq <= cutsq) {
-        neighptr[n] = j;
+        if (history && rsq < radsum*radsum)
-
+          neighptr[n++] = j ^ mask_history;
-        if (fix_history) {
+        else
-          if (rsq < radsum*radsum) {
+          neighptr[n++] = j;
-            for (m = 0; m < npartner[i]; m++)
-              if (partner[i][m] == tag[j]) break;
-            if (m < npartner[i]) {
-              touchptr[n] = 1;
-              memcpy(&shearptr[nn],&shearpartner[i][dnum*m],dnumbytes);
-              nn += dnum;
-            } else {
-              touchptr[n] = 0;
-              memcpy(&shearptr[nn],zeroes,dnumbytes);
-              nn += dnum;
-            }
-          } else {
-            touchptr[n] = 0;
-            memcpy(&shearptr[nn],zeroes,dnumbytes);
-            nn += dnum;
-          }
-        }
-
-        n++;
       }
     }
 
@@ -174,13 +119,6 @@ void NPairHalfSizeNsqNewton::build(NeighList *list)
     ipage->vgot(n);
     if (ipage->status())
       error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
-
-    if (fix_history) {
-      firsttouch[i] = touchptr;
-      firstshear[i] = shearptr;
-      ipage_touch->vgot(n);
-      dpage_shear->vgot(nn);
-    }
   }
 
   list->inum = inum;

src/npair_skip_size.cpp

@@ -17,9 +17,6 @@
 #include "neigh_list.h"
 #include "atom.h"
 #include "atom_vec.h"
-#include "molecule.h"
-#include "domain.h"
-#include "fix_shear_history.h"
 #include "my_page.h"
 #include "error.h"
 
@@ -32,24 +29,13 @@ NPairSkipSize::NPairSkipSize(LAMMPS *lmp) : NPair(lmp) {}
 /* ----------------------------------------------------------------------
    build skip list for subset of types from parent list
    iskip and ijskip flag which atom types and type pairs to skip
-   if list requests it, preserve shear history via fix shear/history
 ------------------------------------------------------------------------- */
 
 void NPairSkipSize::build(NeighList *list)
 {
   int i,j,ii,jj,m,n,nn,itype,jnum,joriginal,dnum,dnumbytes;
   tagint jtag;
-  int *neighptr,*jlist,*touchptr;
+  int *neighptr,*jlist;
-  double *shearptr;
-
-  int *npartner;
-  tagint **partner;
-  double **shearpartner;
-  int **firsttouch;
-  double **firstshear;
-  MyPage<int> *ipage_touch;
-  MyPage<double> *dpage_shear;
-  NeighList *listhistory;
 
   tagint *tag = atom->tag;
   int *type = atom->type;
@@ -68,28 +54,8 @@ void NPairSkipSize::build(NeighList *list)
   int *iskip = list->iskip;
   int **ijskip = list->ijskip;
 
-  FixShearHistory *fix_history = (FixShearHistory *) list->fix_history;
-  if (fix_history) {
-    fix_history->nlocal_neigh = nlocal;
-    fix_history->nall_neigh = nlocal + atom->nghost;
-    npartner = fix_history->npartner;
-    partner = fix_history->partner;
-    shearpartner = fix_history->shearpartner;
-    listhistory = list->listhistory;
-    firsttouch = listhistory->firstneigh;
-    firstshear = listhistory->firstdouble;
-    ipage_touch = listhistory->ipage;
-    dpage_shear = listhistory->dpage;
-    dnum = listhistory->dnum;
-    dnumbytes = dnum * sizeof(double);
-  }
-
   int inum = 0;
   ipage->reset();
-  if (fix_history) {
-    ipage_touch->reset();
-    dpage_shear->reset();
-  }
 
   // loop over atoms in other list
   // skip I atom entirely if iskip is set for type[I]
@@ -102,13 +68,8 @@ void NPairSkipSize::build(NeighList *list)
 
     n = 0;
     neighptr = ipage->vget();
-    if (fix_history) {
-      nn = 0;
-      touchptr = ipage_touch->vget();
-      shearptr = dpage_shear->vget();
-    }
 
-    // loop over parent non-skip size list and optionally its history info
+    // loop over parent non-skip size list
 
     jlist = firstneigh_skip[i];
     jnum = numneigh_skip[i];
@@ -117,29 +78,7 @@ void NPairSkipSize::build(NeighList *list)
       joriginal = jlist[jj];
       j = joriginal & NEIGHMASK;
       if (ijskip[itype][type[j]]) continue;
-      neighptr[n] = joriginal;
+      neighptr[n++] = joriginal;
-
-      // no numeric test for current touch
-      // just use FSH partner list to infer it
-      // would require distance calculation for spheres
-      // more complex calculation for surfs
-
-      if (fix_history) {
-        jtag = tag[j];
-        for (m = 0; m < npartner[i]; m++)
-          if (partner[i][m] == jtag) break;
-        if (m < npartner[i]) {
-          touchptr[n] = 1;
-          memcpy(&shearptr[nn],&shearpartner[i][dnum*m],dnumbytes);
-          nn += dnum;
-        } else {
-          touchptr[n] = 0;
-          memcpy(&shearptr[nn],zeroes,dnumbytes);
-          nn += dnum;
-        }
-      }
-
-      n++;
     }
 
     ilist[inum++] = i;
@@ -148,13 +87,6 @@ void NPairSkipSize::build(NeighList *list)
     ipage->vgot(n);
     if (ipage->status())
       error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
-
-    if (fix_history) {
-      firsttouch[i] = touchptr;
-      firstshear[i] = shearptr;
-      ipage_touch->vgot(n);
-      dpage_shear->vgot(nn);
-    }
   }
 
   list->inum = inum;

src/npair_skip_size_off2on.cpp

@@ -17,9 +17,6 @@
 #include "neigh_list.h"
 #include "atom.h"
 #include "atom_vec.h"
-#include "molecule.h"
-#include "domain.h"
-#include "fix_shear_history.h"
 #include "my_page.h"
 #include "error.h"
 
@@ -33,24 +30,13 @@ NPairSkipSizeOff2on::NPairSkipSizeOff2on(LAMMPS *lmp) : NPair(lmp) {}
    build skip list for subset of types from parent list
    iskip and ijskip flag which atom types and type pairs to skip
    parent non-skip list used newton off, this skip list is newton on
-   if list requests it, preserve shear history via fix shear/history
 ------------------------------------------------------------------------- */
 
 void NPairSkipSizeOff2on::build(NeighList *list)
 {
   int i,j,ii,jj,m,n,nn,itype,jnum,joriginal,dnum,dnumbytes;
   tagint itag,jtag;
-  int *neighptr,*jlist,*touchptr;
+  int *neighptr,*jlist;
-  double *shearptr;
-
-  int *npartner;
-  tagint **partner;
-  double **shearpartner;
-  int **firsttouch;
-  double **firstshear;
-  MyPage<int> *ipage_touch;
-  MyPage<double> *dpage_shear;
-  NeighList *listhistory;
 
   tagint *tag = atom->tag;
   int *type = atom->type;
@@ -69,28 +55,8 @@ void NPairSkipSizeOff2on::build(NeighList *list)
   int *iskip = list->iskip;
   int **ijskip = list->ijskip;
 
-  FixShearHistory *fix_history = (FixShearHistory *) list->fix_history;
-  if (fix_history) {
-    fix_history->nlocal_neigh = nlocal;
-    fix_history->nall_neigh = nlocal + atom->nghost;
-    npartner = fix_history->npartner;
-    partner = fix_history->partner;
-    shearpartner = fix_history->shearpartner;
-    listhistory = list->listhistory;
-    firsttouch = listhistory->firstneigh;
-    firstshear = listhistory->firstdouble;
-    ipage_touch = listhistory->ipage;
-    dpage_shear = listhistory->dpage;
-    dnum = listhistory->dnum;
-    dnumbytes = dnum * sizeof(double);
-  }
-
   int inum = 0;
   ipage->reset();
-  if (fix_history) {
-    ipage_touch->reset();
-    dpage_shear->reset();
-  }
 
   // loop over atoms in other list
   // skip I atom entirely if iskip is set for type[I]
@@ -104,11 +70,6 @@ void NPairSkipSizeOff2on::build(NeighList *list)
 
     n = 0;
     neighptr = ipage->vget();
-    if (fix_history) {
-      nn = 0;
-      touchptr = ipage_touch->vget();
-      shearptr = dpage_shear->vget();
-    }
 
     // loop over parent non-skip size list and optionally its history info
 
@@ -125,28 +86,7 @@ void NPairSkipSizeOff2on::build(NeighList *list)
       jtag = tag[j];
       if (j >= nlocal && jtag < itag) continue;
 
-      neighptr[n] = joriginal;
+      neighptr[n++] = joriginal;
-
-      // no numeric test for current touch
-      // just use FSH partner list to infer it
-      // would require distance calculation for spheres
-      // more complex calculation for surfs
-
-      if (fix_history) {
-        for (m = 0; m < npartner[i]; m++)
-          if (partner[i][m] == jtag) break;
-        if (m < npartner[i]) {
-          touchptr[n] = 1;
-          memcpy(&shearptr[nn],&shearpartner[i][dnum*m],dnumbytes);
-          nn += dnum;
-        } else {
-          touchptr[n] = 0;
-          memcpy(&shearptr[nn],zeroes,dnumbytes);
-          nn += dnum;
-        }
-      }
-
-      n++;
     }
 
     ilist[inum++] = i;
@@ -155,13 +95,6 @@ void NPairSkipSizeOff2on::build(NeighList *list)
     ipage->vgot(n);
     if (ipage->status())
       error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
-
-    if (fix_history) {
-      firsttouch[i] = touchptr;
-      firstshear[i] = shearptr;
-      ipage_touch->vgot(n);
-      dpage_shear->vgot(nn);
-    }
   }
 
   list->inum = inum;

src/npair_skip_size_off2on_oneside.cpp

@@ -17,9 +17,7 @@
 #include "neigh_list.h"
 #include "atom.h"
 #include "atom_vec.h"
-#include "molecule.h"
 #include "domain.h"
-#include "fix_shear_history.h"
 #include "my_page.h"
 #include "error.h"
 
@@ -35,7 +33,6 @@ NPairSkipSizeOff2onOneside::NPairSkipSizeOff2onOneside(LAMMPS *lmp) :
    iskip and ijskip flag which atom types and type pairs to skip
    parent non-skip list used newton off and was not onesided,
      this skip list is newton on and onesided
-   if list requests it, preserve shear history via fix shear/history
 ------------------------------------------------------------------------- */
 
 void NPairSkipSizeOff2onOneside::build(NeighList *list)
@@ -44,15 +41,6 @@ void NPairSkipSizeOff2onOneside::build(NeighList *list)
   tagint jtag;
   int *surf,*jlist;
 
-  int *npartner;
-  tagint **partner;
-  double **shearpartner;
-  int **firsttouch;
-  double **firstshear;
-  MyPage<int> *ipage_touch;
-  MyPage<double> *dpage_shear;
-  NeighList *listhistory;
-
   tagint *tag = atom->tag;
   int *type = atom->type;
   int nlocal = atom->nlocal;
@@ -73,28 +61,8 @@ void NPairSkipSizeOff2onOneside::build(NeighList *list)
   if (domain->dimension == 2) surf = atom->line;
   else surf = atom->tri;
 
-  FixShearHistory *fix_history = (FixShearHistory *) list->fix_history;
-  if (fix_history) {
-    fix_history->nlocal_neigh = nlocal;
-    fix_history->nall_neigh = nlocal + atom->nghost;
-    npartner = fix_history->npartner;
-    partner = fix_history->partner;
-    shearpartner = fix_history->shearpartner;
-    listhistory = list->listhistory;
-    firsttouch = listhistory->firstneigh;
-    firstshear = listhistory->firstdouble;
-    ipage_touch = listhistory->ipage;
-    dpage_shear = listhistory->dpage;
-    dnum = listhistory->dnum;
-    dnumbytes = dnum * sizeof(double);
-  }
-
   int inum = 0;
   ipage->reset();
-  if (fix_history) {
-    ipage_touch->reset();
-    dpage_shear->reset();
-  }
 
   // two loops over parent list required, one to count, one to store
   // because onesided constraint means pair I,J may be stored with I or J
@@ -139,7 +107,7 @@ void NPairSkipSizeOff2onOneside::build(NeighList *list)
     }
   }
 
-  // allocate all per-atom neigh list chunks, including history
+  // allocate all per-atom neigh list chunks
 
   for (i = 0; i < nlocal; i++) {
     if (numneigh[i] == 0) continue;
@@ -147,10 +115,6 @@ void NPairSkipSizeOff2onOneside::build(NeighList *list)
     firstneigh[i] = ipage->get(n);
     if (ipage->status())
       error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
-    if (fix_history) {
-      firsttouch[i] = ipage_touch->get(n);
-      firstshear[i] = dpage_shear->get(dnum*n);
-    }
   }
 
   // second loop over atoms in other list to store neighbors
@@ -189,32 +153,11 @@ void NPairSkipSizeOff2onOneside::build(NeighList *list)
       // OK, b/c there is no special list flagging for surfs
 
       firstneigh[i][numneigh[i]] = j;
-
-      // no numeric test for current touch
-      // just use FSH partner list to infer it
-      // would require complex calculation for surfs
-
-      if (fix_history) {
-        jtag = tag[j];
-        n = numneigh[i];
-        nn = dnum*n;
-        for (m = 0; m < npartner[i]; m++)
-          if (partner[i][m] == jtag) break;
-        if (m < npartner[i]) {
-          firsttouch[i][n] = 1;
-          memcpy(&firstshear[i][nn],&shearpartner[i][dnum*m],dnumbytes);
-        } else {
-          firsttouch[i][n] = 0;
-          memcpy(&firstshear[i][nn],zeroes,dnumbytes);
-        }
-      }
-
       numneigh[i]++;
       if (flip) i = j;
     }
 
     // only add atom I to ilist if it has neighbors
-    // fix shear/history allows for this in pre_exchange_onesided()
 
     if (numneigh[i]) ilist[inum++] = i;
   }

src/verlet.cpp

@@ -90,10 +90,9 @@ void Verlet::setup(int flag)
   if (comm->me == 0 && screen) {
     fprintf(screen,"Setting up Verlet run ...\n");
     if (flag) {
-      fprintf(screen,"  Unit style    : %s\n", update->unit_style);
+      fprintf(screen,"  Unit style    : %s\n",update->unit_style);
-      fprintf(screen,"  Current step  : " BIGINT_FORMAT "\n",
+      fprintf(screen,"  Current step  : " BIGINT_FORMAT "\n",update->ntimestep);
-              update->ntimestep);
+      fprintf(screen,"  Time step     : %g\n",update->dt);
-      fprintf(screen,"  Time step     : %g\n", update->dt);
       timer->print_timeout(screen);
     }
   }
@@ -122,6 +121,7 @@ void Verlet::setup(int flag)
   domain->box_too_small_check();
   modify->setup_pre_neighbor();
   neighbor->build();
+  modify->setup_post_neighbor();
   neighbor->ncalls = 0;
 
   // compute all forces
@@ -183,6 +183,7 @@ void Verlet::setup_minimal(int flag)
     domain->box_too_small_check();
     modify->setup_pre_neighbor();
     neighbor->build();
+    modify->setup_post_neighbor();
     neighbor->ncalls = 0;
   }
 
@@ -227,6 +228,7 @@ void Verlet::run(int n)
   int n_post_integrate = modify->n_post_integrate;
   int n_pre_exchange = modify->n_pre_exchange;
   int n_pre_neighbor = modify->n_pre_neighbor;
+  int n_post_neighbor = modify->n_post_neighbor;
   int n_pre_force = modify->n_pre_force;
   int n_pre_reverse = modify->n_pre_reverse;
   int n_post_force = modify->n_post_force;
@@ -284,6 +286,10 @@ void Verlet::run(int n)
       }
       neighbor->build();
       timer->stamp(Timer::NEIGH);
+      if (n_post_neighbor) {
+        modify->post_neighbor();
+        timer->stamp(Timer::MODIFY);
+      }
     }
 
     // force computations