Initial files

src/fix_pair_tracker.cpp

@@ -0,0 +1,444 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include <string.h>
+#include "fix_pair_tracking.h"
+#include "pair_tracker.h"
+#include "atom.h"
+#include "atom_vec.h"
+#include "update.h"
+#include "force.h"
+#include "pair.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "memory.h"
+#include "modify.h"
+#include "error.h"
+#include "group.h"
+#include "domain.h"
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+
+#define DELTA 100
+
+/* ---------------------------------------------------------------------- */
+
+FixPairTracking::FixPairTracking(LAMMPS *lmp, int narg, char **arg) :
+  Fix(lmp, narg, arg),
+  nvalues(0),
+  array(NULL), vector(NULL), pack_choice(NULL)
+{
+  if (narg < 3) error->all(FLERR,"Illegal fix pair/tracker command");
+  store_flag = 0;
+  local_flag = 1;
+  nvalues = narg - 4; 
+  
+  nevery = utils::inumeric(FLERR,arg[3],false,lmp);
+  if (nevery <= 0) error->all(FLERR,"Illegal fix pair/tracker command");  
+    
+  if (nvalues == 1) size_local_cols = 0;
+  else size_local_cols = nvalues;
+
+  pack_choice = new FnPtrPack[nvalues];
+  
+  int i;
+  for (int iarg = 4; iarg < narg; iarg++) {
+    i = iarg - 4;
+    
+    if (strcmp(arg[iarg],"id1") == 0) {
+      pack_choice[i] = &FixPairTracking::pack_id1;
+    } else if (strcmp(arg[iarg],"id2") == 0) {
+      pack_choice[i] = &FixPairTracking::pack_id2;
+           
+    } else if (strcmp(arg[iarg],"time/created") == 0) {
+      pack_choice[i] = &FixPairTracking::pack_time_created;
+    } else if (strcmp(arg[iarg],"time/broken") == 0) {
+      pack_choice[i] = &FixPairTracking::pack_time_broken;
+    } else if (strcmp(arg[iarg],"time/total") == 0) {
+      pack_choice[i] = &FixPairTracking::pack_time_total;
+      
+    } else if (strcmp(arg[iarg],"x") == 0) {
+      pack_choice[i] = &FixPairTracking::pack_x;    
+    } else if (strcmp(arg[iarg],"y") == 0) {
+      pack_choice[i] = &FixPairTracking::pack_y;    
+    } else if (strcmp(arg[iarg],"z") == 0) {
+      pack_choice[i] = &FixPairTracking::pack_z;    
+      
+    } else if (strcmp(arg[iarg],"xstore") == 0) {
+      pack_choice[i] = &FixPairTracking::pack_xstore;    
+      store_flag = 1;
+    } else if (strcmp(arg[iarg],"ystore") == 0) {
+      pack_choice[i] = &FixPairTracking::pack_ystore;    
+      store_flag = 1;
+    } else if (strcmp(arg[iarg],"zstore") == 0) {
+      pack_choice[i] = &FixPairTracking::pack_zstore;      
+      store_flag = 1;      
+   
+    } else if (strcmp(arg[iarg],"rmin") == 0) {
+      pack_choice[i] = &FixPairTracking::pack_rmin;    
+    } else if (strcmp(arg[iarg],"rave") == 0) {
+      pack_choice[i] = &FixPairTracking::pack_rave;     
+   
+    } else error->all(FLERR, "Invalid keyword in fix pair/tracker command");
+  }
+
+  nmax = 0;
+  ncount = 0;
+  vector = NULL;
+  array = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixPairTracking::~FixPairTracking()
+{
+  if (modify->nfix & store_flag == 1) {
+    modify->delete_fix(id_fix);
+    delete [] id_fix;
+  }
+  
+  delete [] pack_choice;
+
+  memory->destroy(vector);  
+  memory->destroy(array);
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixPairTracking::setmask()
+{
+  int mask = 0;
+  mask |= POST_FORCE;
+  return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixPairTracking::post_constructor()
+{
+  //If use stored x,y,z values, use store property (can transfer to ghost atoms) to store positions
+  
+  if(store_flag == 1){
+      
+    lx = domain->xprd;
+    ly = domain->yprd;
+    lz = domain->zprd;        
+
+    int nn = strlen(id) + strlen("_FIX_PROP_ATOM") + 1;
+    id_property_fix = new char[nn];
+    strcpy(id_property_fix,id);
+    strcat(id_property_fix,"_FIX_PROP_ATOM");
+    
+    int ifix = modify->find_fix(id_property_fix);
+    if (ifix < 0) {
+    
+      int n_x = strlen(id) + 4;
+      
+      char * lab1 = new char[n_x];
+      strcpy(lab1, "d_x");
+      strcat(lab1, id);
+      char * lab2 = new char[n_x];
+      strcpy(lab2, "d_y");
+      strcat(lab2, id);
+      char * lab3 = new char[n_x];
+      strcpy(lab3, "d_z");
+      strcat(lab3, id);
+        
+      char **newarg = new char*[8];
+      newarg[0] = id_property_fix;
+      newarg[1] = group->names[igroup];
+      newarg[2] = (char *) "property/atom"; 
+      newarg[3] = (char *) lab1;         
+      newarg[4] = (char *) lab2;         
+      newarg[5] = (char *) lab3; 
+      newarg[6] = (char *) "ghost";
+      newarg[7] = (char *) "yes";
+
+      modify->add_fix(8,newarg); 
+      //Needs ghost atoms to calculate CoM
+
+      int type_flag;
+      int col_flag;
+
+      strcpy(lab1, "x");
+      strcat(lab1, id);
+      strcpy(lab2, "y");
+      strcat(lab2, id);
+      strcpy(lab3, "z");
+      strcat(lab3, id);
+
+      index_x = atom->find_custom(lab1, type_flag, col_flag);
+      index_y = atom->find_custom(lab2, type_flag, col_flag);
+      index_z = atom->find_custom(lab3, type_flag, col_flag);
+      delete [] newarg;    
+      delete [] lab1;
+      delete [] lab2;
+      delete [] lab3;
+    } 
+    
+    ifix = modify->find_fix(id_fix);
+    if (ifix < 0) error->all(FLERR,"Could not find fix ID for fix broken/bond");
+    if (modify->fix[ifix]->restart_reset) {
+        modify->fix[ifix]->restart_reset = 0;
+    } else {
+
+      double *xi = atom->dvector[index_x];
+      double *yi = atom->dvector[index_y];
+      double *zi = atom->dvector[index_z];
+      
+      double **xs = atom->x;
+      int nlocal = atom->nlocal;
+      
+      for (int i = 0; i < nlocal; i++) {
+        xi[i] = xs[i][0];
+        yi[i] = xs[i][1];
+        zi[i] = xs[i][2];   
+      }
+    }    
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixPairTracking::init()
+{
+  // Set size of array/vector
+  ncount = 0;
+  
+  if (ncount > nmax) {
+      reallocate(ncount);}
+  size_local_rows = ncount;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixPairTracking::lost_contact(int i, int j, double n, double rs, double rm)
+{    
+  if (ncount == nmax) reallocate(ncount);
+  
+  index_i = i;
+  index_j = j;
+  rmin = rm;
+  rave = ra;
+  ntimestep = n;
+  
+  // fill vector or array with local values
+  if (nvalues == 1) {
+    (this->*pack_choice[0])(0);
+  } else {
+    for (int n = 0; n < nvalues; n++)
+      (this->*pack_choice[n])(n); 
+  }  
+  
+  ncount += 1;  
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixPairTracking::post_force(int /*vflag*/) 
+{    
+  if (update->ntimestep % nevery == 0) {
+    size_local_rows = ncount;
+    ncount = 0;  
+  }
+}
+
+
+/* ---------------------------------------------------------------------- */
+
+void FixPairTracking::reallocate(int n)
+{
+  // grow vector or array 
+  while (nmax <= n) nmax += DELTA;
+  
+  if (nvalues == 1) {
+    memory->grow(vector,nmax,"fix_broken_bonds:vector");
+    vector_local = vector;
+  } else {
+    memory->grow(array,nmax,nvalues,"fix_broken_bonds:array");
+    array_local = array;
+  }
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of local data
+------------------------------------------------------------------------- */
+
+double FixPairTracking::memory_usage()
+{
+  double bytes = nmax*nvalues * sizeof(double);
+  bytes += nmax*2 * sizeof(int);
+  return bytes;
+}
+
+/* ----------------------------------------------------------------------
+   one method for every keyword fix property/local can output
+   the atom property is packed into buf starting at n with stride nvalues
+   customize a new keyword by adding a method
+------------------------------------------------------------------------- */
+
+/* ---------------------------------------------------------------------- */
+
+void FixPairTracking::pack_time_created(int n) 
+{
+  if (nvalues == 1)
+    vector[ncount] = ntimestep;
+  else
+    array[ncount][n] = ntimestep;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixPairTracking::pack_time_broken(int n) 
+{
+  if (nvalues == 1)
+    vector[ncount] = update->ntimestep;
+  else
+    array[ncount][n] = update->ntimestep;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixPairTracking::pack_time_total(int n) 
+{
+  if (nvalues == 1)
+    vector[ncount] = update->ntimestep-ntimestep;
+  else
+    array[ncount][n] = update->ntimestep-ntimestep;
+}
+
+
+/* ---------------------------------------------------------------------- */
+
+void FixPairTracking::pack_id1(int n) 
+{
+  tagint *tag = atom->tag;
+  
+  if (nvalues == 1)
+    vector[ncount] = tag[index_i];
+  else
+    array[ncount][n] = tag[index_i];
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixPairTracking::pack_id2(int n)
+{
+  tagint *tag = atom->tag;
+
+  if (nvalues == 1)
+    vector[ncount] = tag[index_j];
+  else
+    array[ncount][n] = tag[index_j];
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixPairTracking::pack_x(int n)
+{
+  double lx_new = domain->xprd;
+  double **x = atom->x; 
+    
+  if (nvalues == 1)
+    vector[ncount] = (x[index_i][0] + x[index_j][0])/2;
+  else
+    array[ncount][n] = (x[index_i][0] + x[index_j][0])/2;
+}
+
+
+/* ---------------------------------------------------------------------- */
+
+void FixPairTracking::pack_y(int n)
+{
+  double lx_new = domain->yprd;
+  double **x = atom->x; 
+  
+  if (nvalues == 1)
+    vector[ncount] = (x[index_i][1] + x[index_j][1])/2;
+  else
+    array[ncount][n] = (x[index_i][1] + x[index_j][1])/2;
+}
+
+
+/* ---------------------------------------------------------------------- */
+
+void FixPairTracking::pack_z(int n)
+{
+  double lx_new = domain->zprd;
+  double **x = atom->x; 
+    
+  if (nvalues == 1)
+    vector[ncount] = (x[index_i][2] + x[index_j][2])/2;
+  else
+    array[ncount][n] = (x[index_i][2] + x[index_j][2])/2;
+}
+
+
+/* ---------------------------------------------------------------------- */
+
+void FixPairTracking::pack_xstore(int n)
+{
+  double *x = atom->dvector[index_x];
+        
+  if (nvalues == 1)
+    vector[ncount] = (x[index_i] + x[index_j])/2;
+  else
+    array[ncount][n] = (x[index_i] + x[index_j])/2;
+}
+
+
+/* ---------------------------------------------------------------------- */
+
+void FixPairTracking::pack_ystore(int n)
+{
+  double *y = atom->dvector[index_y];
+        
+  if (nvalues == 1)
+    vector[ncount] = (y[index_i] + y[index_j])/2;
+  else
+    array[ncount][n] = (y[index_i] + y[index_j])/2;
+}
+
+
+/* ---------------------------------------------------------------------- */
+
+void FixPairTracking::pack_zstore(int n)
+{
+  double *z = atom->dvector[index_z];
+        
+  if (nvalues == 1)
+    vector[ncount] = (z[index_i] + z[index_j])/2;
+  else
+    array[ncount][n] = (z[index_i] + z[index_j])/2;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixPairTracking::pack_rmin(int n)
+{
+  if (nvalues == 1)
+    vector[ncount] = rmin;
+  else
+    array[ncount][n] = rmin;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixPairTracking::pack_rave(int n)
+{
+  if (nvalues == 1)
+    vector[ncount] = rsum/(update->ntimestep-ntimestep);
+  else
+    array[ncount][n] = rsum/(update->ntimestep-ntimestep);
+}

src/fix_pair_tracker.h

@@ -0,0 +1,119 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+
+FixStyle(pair/tracking,FixPairTracking)
+
+#else
+
+#ifndef LMP_FIX_PAIR_TRACKING_H
+#define LMP_FIX_PAIR_TRACKING_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixPairTracking : public Fix {
+ public: 
+  FixPairTracking(class LAMMPS *, int, char **);
+  ~FixPairTracking();
+  int setmask();  
+  void post_constructor();
+  void init();
+  void post_force(int);
+  double memory_usage();
+  void lost_contact(int, int, double, double, double);
+    
+ private:
+  int nvalues, never;
+  int nmax;
+  int store_flag;
+  int index_i, index_j;
+  
+  double *vector;
+  double **array;
+  
+  double lx;
+  double ly;
+  double lz;    
+
+  int ncount;
+
+  void reallocate(int);
+
+  typedef void (FixPairTracker::*FnPtrPack)(int);
+  FnPtrPack *pack_choice;              // ptrs to pack functions
+
+  void pack_id1(int);
+  void pack_id2(int);
+  
+  void pack_time_created(int);  
+  void pack_time_broken(int);  
+  void pack_time_total(int);  
+  
+  void pack_x(int);  
+  void pack_y(int);  
+  void pack_z(int);  
+  
+  void pack_xstore(int);  
+  void pack_ystore(int);  
+  void pack_zstore(int);   
+  
+  void pack_rmin(int);
+  void pack_rave(int);
+  
+  char *id_fix;
+  int index_x, index_y, index_z;  
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Fix property/local cannot use these inputs together
+
+Only inputs that generate the same number of datums can be used
+togther.  E.g. bond and angle quantities cannot be mixed.
+
+E: Invalid keyword in compute property/local command
+
+Self-explanatory.
+
+E: Fix property/local does not (yet) work with atom_style template
+
+Self-explanatory.
+
+E: Fix property/local for property that isn't allocated
+
+Self-explanatory.
+
+E: No pair style is defined for compute property/local
+
+Self-explanatory.
+
+E: Pair style does not support compute property/local
+
+The pair style does not have a single() function, so it can
+not be invoked by fix bond/swap.
+
+*/

src/pair_tracker.cpp

@@ -0,0 +1,469 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://lammps.sandia.gov/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing authors: Leo Silbert (SNL), Gary Grest (SNL)
+------------------------------------------------------------------------- */
+
+#include "pair_tracker.h"
+
+#include <cmath>
+#include <cstring>
+
+#include "atom.h"
+#include "force.h"
+#include "update.h"
+#include "modify.h"
+#include "fix.h"
+#include "fix_dummy.h"
+#include "fix_neigh_history.h"
+#include "comm.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+PairTracker::PairTracker(LAMMPS *lmp) : Pair(lmp)
+{
+  single_enable = 1;
+  no_virial_fdotr_compute = 1;
+  history = 1;
+  size_history = 3;
+
+  neighprev = 0;
+  nondefault_history_transfer = 1;
+  
+  finitecutflag = 0;
+
+  // create dummy fix as placeholder for FixNeigTRACKistory
+  // this is so final order of Modify:fix will conform to input script
+
+  fix_history = nullptr;
+  modify->add_fix("NEIGH_HISTORY_TRACK_DUMMY all DUMMY");
+  fix_dummy = (FixDummy *) modify->fix[modify->nfix-1];
+}
+
+/* ---------------------------------------------------------------------- */
+
+PairTracker::~PairTracker()
+{
+  if (!fix_history) modify->delete_fix("NEIGH_HISTORY_TRACK_DUMMY");
+  else modify->delete_fix("NEIGH_HISTORY_TRACK");
+
+  if (allocated) {
+    memory->destroy(setflag);
+    memory->destroy(cutsq);
+    memory->destroy(cut);
+
+    delete [] onerad_dynamic;
+    delete [] onerad_frozen;
+    delete [] maxrad_dynamic;
+    delete [] maxrad_frozen;    
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairTracker::compute(int eflag, int vflag)
+{
+  int i,j,ii,jj,inum,jnum;
+  double xtmp,ytmp,ztmp,delx,dely,delz;
+  double radi,radj,radsum,rsq,r;
+  int *ilist,*jlist,*numneigh,**firstneigh;
+  int *touch,**firsttouch;
+  double *data,*alldata,**firstdata;
+
+  int updateflag = 1;
+  if (update->setupflag) updateflag = 0;
+
+
+  double **x = atom->x;
+  double *radius = atom->radius;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+  int newton_pair = force->newton_pair;
+
+  inum = list->inum;
+  ilist = list->ilist;
+  numneigh = list->numneigh;
+  firstneigh = list->firstneigh;
+  firsttouch = fix_history->firstflag;
+  firstdata = fix_history->firstvalue;
+
+  // loop over neighbors of my atoms
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+    if (finitecutflag) radi = radius[i];
+    touch = firsttouch[i];
+    alldata = firstdata[i];
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+
+    for (jj = 0; jj < jnum; jj++) {
+      j = jlist[jj];
+      j &= NEIGHMASK;
+
+      delx = xtmp - x[j][0];
+      dely = ytmp - x[j][1];
+      delz = ztmp - x[j][2];
+      rsq = delx*delx + dely*dely + delz*delz;
+      if (finitecutflag) {
+        radj = radius[j];
+        radsum = radi + radj;
+        
+        if (rsq >= radsum*radsum) {
+          // unset non-touching neighbors
+          if(touch[jj] == 1) {
+            fix_pair_tracker->lost_contact(i, j, data[0], data[1], data[2]);
+          }
+          touch[jj] = 0;
+          data = &alldata[size_history*jj];
+          data[0] = 0.0; // initial time
+          data[1] = 0.0; // sum of r
+          data[2] = 0.0; // min of r
+                
+        } else {        
+
+          data = &alldata[size_history*jj];
+          if (touch[jj] == 0) {
+            data[0] = update->ntimestep;
+            data[1] = 0.0;
+            data[2] = 0.0;
+          }
+          touch[jj] = 1;
+
+          if (updateflag) {
+            r = sqrt(rsq);
+            data[1] += r;
+            if(data[2] > r) data[2] = r;
+          }
+        }
+      } else {
+          
+      }
+
+      if (rsq >= radsum*radsum) {
+
+
+        f[i][0] += fx;
+        f[i][1] += fy;
+        f[i][2] += fz;
+
+
+        if (newton_pair || j < nlocal) {
+          f[j][0] -= fx;
+          f[j][1] -= fy;
+          f[j][2] -= fz;
+          torque[j][0] -= radj*tor1;
+          torque[j][1] -= radj*tor2;
+          torque[j][2] -= radj*tor3;
+        }
+
+      }
+    }
+  }
+}
+
+/* ----------------------------------------------------------------------
+   allocate all arrays
+------------------------------------------------------------------------- */
+
+void PairTracker::allocate()
+{
+  allocated = 1;
+  int n = atom->ntypes;
+
+  memory->create(setflag,n+1,n+1,"pair:setflag");
+  for (int i = 1; i <= n; i++)
+    for (int j = i; j <= n; j++)
+      setflag[i][j] = 0;
+
+  memory->create(cutsq,n+1,n+1,"pair:cutsq");
+  memory->create(cut,n+1,n+1,"pair:cut");
+
+  onerad_dynamic = new double[n+1];
+  onerad_frozen = new double[n+1];
+  maxrad_dynamic = new double[n+1];
+  maxrad_frozen = new double[n+1];  
+}
+
+/* ----------------------------------------------------------------------
+   global settings
+------------------------------------------------------------------------- */
+
+void PairTracker::settings(int narg, char **arg)
+{
+  if (narg != 0 || narg != 1) error->all(FLERR,"Illegal pair_style command");
+
+  if (narg == 1) {
+    if (strcmp(arg[0], "radius")) finitecutflag = 1;
+    else error->all(FLERR,"Illegal pair_style command");
+  }
+}
+
+/* ----------------------------------------------------------------------
+   set coeffs for one or more type pairs
+------------------------------------------------------------------------- */
+
+void PairTracker::coeff(int narg, char **arg)
+{
+  if (narg > 2 && finitecutflag) error->all(FLERR,"Incorrect args for pair coefficients");
+  if (narg != 3 && !finitecutflag) error->all(FLERR,"Incorrect args for pair coefficients");
+  if (!allocated) allocate();
+
+  int ilo,ihi,jlo,jhi;
+  utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error);
+  utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error);
+
+  double cut_one = 0.0;
+  if (!finitecutflag) cut_one = utils::numeric(FLERR,arg[3],false,lmp); 
+
+  int count = 0;
+  for (int i = ilo; i <= ihi; i++) {
+    for (int j = MAX(jlo,i); j <= jhi; j++) {
+      setflag[i][j] = 1;
+      cut[i][j] = cut_one;
+      count++;
+    }
+  }
+
+  if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
+}
+
+/* ----------------------------------------------------------------------
+   init specific to this pair style
+------------------------------------------------------------------------- */
+
+void PairTracker::init_style()
+{
+  int i;
+
+  // error and warning checks
+
+  if (!atom->radius_flag && finitecutflag)
+    error->all(FLERR,"Pair tracker requires atom attribute radius for finite cutoffs");
+
+  // need a history neigh list
+
+  int irequest = neighbor->request(this,instance_me);
+  if (finitecutflag) neighbor->requests[irequest]->size = 1;
+  neighbor->requests[irequest]->history = 1;
+
+  // if history is stored and first init, create Fix to store history
+  // it replaces FixDummy, created in the constructor
+  // this is so its order in the fix list is preserved
+
+  if (fix_history == nullptr) {
+    char dnumstr[16];
+    sprintf(dnumstr,"%d",size_history);
+    char **fixarg = new char*[4];
+    fixarg[0] = (char *) "NEIGH_HISTORY_TRACK";
+    fixarg[1] = (char *) "all";
+    fixarg[2] = (char *) "NEIGH_HISTORY";
+    fixarg[3] = dnumstr;
+    modify->replace_fix("NEIGH_HISTORY_TRACK_DUMMY",4,fixarg,1);
+    delete [] fixarg;
+    int ifix = modify->find_fix("NEIGH_HISTORY_TRACK");
+    fix_history = (FixNeigTRACKistory *) modify->fix[ifix];
+    fix_history->pair = this;
+  }
+
+  if (finitecutflag) {
+    // set maxrad_dynamic and maxrad_frozen for each type
+    // include future FixPour and FixDeposit particles as dynamic
+    
+    int itype;
+    for (i = 1; i <= atom->ntypes; i++) {
+      onerad_dynamic[i] = onerad_frozen[i] = 0.0;
+      if (ipour >= 0) {
+        itype = i;
+        onerad_dynamic[i] =
+          *((double *) modify->fix[ipour]->extract("radius",itype));
+      }
+      if (idep >= 0) {
+        itype = i;
+        onerad_dynamic[i] =
+          *((double *) modify->fix[idep]->extract("radius",itype));
+      }
+    }
+    
+    double *radius = atom->radius;
+    int *mask = atom->mask;
+    int *type = atom->type;
+    int nlocal = atom->nlocal;
+    
+    for (i = 0; i < nlocal; i++)
+      if (mask[i] & freeze_group_bit)
+        onerad_frozen[type[i]] = MAX(onerad_frozen[type[i]],radius[i]);
+      else
+        onerad_dynamic[type[i]] = MAX(onerad_dynamic[type[i]],radius[i]);
+    
+    MPI_Allreduce(&onerad_dynamic[1],&maxrad_dynamic[1],atom->ntypes,
+                  MPI_DOUBLE,MPI_MAX,world);
+    MPI_Allreduce(&onerad_frozen[1],&maxrad_frozen[1],atom->ntypes,
+                  MPI_DOUBLE,MPI_MAX,world);
+  }
+
+  int ifix = modify->find_fix("NEIGH_HISTORY_TRACK");
+  if (ifix < 0) error->all(FLERR,"Could not find pair fix neigh history ID");
+  fix_history = (FixNeigTRACKistory *) modify->fix[ifix];
+   
+  ifix = modify->find_fix_by_style("pair/tracker");
+  if(ifix < 0) error->all(FLERR,"Cannot use pair tracker without fix pair/tracker");
+  fix_pair_tracker = (FixPairTracker *) modify->fix[ifix];  
+}
+
+/* ----------------------------------------------------------------------
+   init for one type pair i,j and corresponding j,i
+------------------------------------------------------------------------- */
+
+double PairTracker::init_one(int i, int j)
+{
+  if (!allocated) allocate();
+
+  // always mix prefactors geometrically
+
+  if (setflag[i][j] == 0) {
+    cut[i][j] = mix_distance(cut[i][i],cut[j][j]);
+  }
+
+  cut[j][i] = cut[i][j];
+
+  double cutoff; 
+  
+  // if finite, cutoff = sum of max I,J radii for
+  // dynamic/dynamic & dynamic/frozen interactions, but not frozen/frozen
+  if (finitecutflag) {
+    cutoff = maxrad_dynamic[i]+maxrad_dynamic[j];
+    cutoff = MAX(cutoff,maxrad_frozen[i]+maxrad_dynamic[j]);
+    cutoff = MAX(cutoff,maxrad_dynamic[i]+maxrad_frozen[j]);
+  } else { 
+    cutoff = cut[i][j];
+  }
+  return cutoff;
+}
+
+/* ----------------------------------------------------------------------
+  proc 0 writes to restart file
+------------------------------------------------------------------------- */
+
+void PairTracker::write_restart(FILE *fp)
+{
+  write_restart_settings(fp);
+
+  int i,j;
+  for (i = 1; i <= atom->ntypes; i++)
+    for (j = i; j <= atom->ntypes; j++) {
+      fwrite(&setflag[i][j],sizeof(int),1,fp);
+      if (setflag[i][j]) {
+        fwrite(&cut[i][j],sizeof(double),1,fp);
+      }        
+    }
+}
+
+/* ----------------------------------------------------------------------
+  proc 0 reads from restart file, bcasts
+------------------------------------------------------------------------- */
+
+void PairTracker::read_restart(FILE *fp)
+{
+  read_restart_settings(fp);
+  allocate();
+
+  int i,j;
+  int me = comm->me;
+  for (i = 1; i <= atom->ntypes; i++)
+    for (j = i; j <= atom->ntypes; j++) {
+      if (me == 0) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,nullptr,error);
+      MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
+      if (setflag[i][j]) {
+        if (me == 0) {
+          utils::sfread(FLERR,&cut[i][j],sizeof(double),1,fp,nullptr,error);
+        }
+        MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
+      }      
+    }
+}
+
+/* ----------------------------------------------------------------------
+  proc 0 writes to restart file
+------------------------------------------------------------------------- */
+
+void PairTracker::write_restart_settings(FILE *fp)
+{
+  fwrite(&mix_flag,sizeof(int),1,fp);
+}
+
+/* ----------------------------------------------------------------------
+  proc 0 reads from restart file, bcasts
+------------------------------------------------------------------------- */
+
+void PairTracker::read_restart_settings(FILE *fp)
+{
+  if (comm->me == 0) {
+    utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,nullptr,error);
+  }
+  MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairTracker::single(int i, int j, int /*itype*/, int /*jtype*/,
+                                    double rsq,
+                                    double /*factor_coul*/, double /*factor_lj*/,
+                                    double &fforce)
+{
+  return 0.0;
+}
+
+/* ----------------------------------------------------------------------
+   transfer history during fix/neigh/history exchange
+   only needed if any history entries i-j are not just negative of j-i entries
+------------------------------------------------------------------------- */
+
+void PairTracker::transfer_history(double* source, double* target)
+{
+  for (int i = 0; i < size_history; i++)
+    target[i] = source[i];
+}
+
+
+/* ----------------------------------------------------------------------
+   self-interaction range of particle if finite particles
+------------------------------------------------------------------------- */
+
+double PairTracker::atom2cut(int i)
+{
+  double cut = atom->radius[i]*2;
+  return cut;
+}
+
+/* ----------------------------------------------------------------------
+   maximum interaction range for two finite particles
+------------------------------------------------------------------------- */
+
+double PairTracker::radii2cut(double r1, double r2)
+{
+  double cut = r1+r2;
+  return cut;
+}
+

src/pair_tracker.h

@@ -0,0 +1,101 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://lammps.sandia.gov/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+
+PairStyle(tracker,PairTracker)
+
+#else
+
+#ifndef LMP_PAIR_TRACKER_H
+#define LMP_PAIR_TRACKER_H
+
+#include "pair.h"
+
+namespace LAMMPS_NS {
+
+class PairTracker : public Pair {
+ public:
+  PairTracker(class LAMMPS *);
+  virtual ~PairTracker();
+  virtual void compute(int, int);
+  virtual void settings(int, char **);
+  void coeff(int, char **);
+  void init_style();
+  double init_one(int, int);
+  void write_restart(FILE *);
+  void read_restart(FILE *);
+  void write_restart_settings(FILE *);
+  void read_restart_settings(FILE *);
+  void reset_dt();
+  virtual double single(int, int, int, int, double, double, double, double &);
+  double atom2cut(int);
+  double radii2cut(double,double);
+
+ protected:
+  int sizeflag;
+  int history;
+  int size_history;
+  int neighprev;
+  double **cut;
+  double *onerad_dynamic,*onerad_frozen;
+  double *maxrad_dynamic,*maxrad_frozen;
+
+  class FixDummy *fix_dummy;
+  class FixNeighHistory *fix_history;
+  class FixPairTracker *fix_pair_tracker;
+
+  void transfer_history(double*, double*);
+  void allocate();
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Incorrect args for pair coefficients
+
+Self-explanatory.  Check the input script or data file.
+
+E: Pair granular requires atom attributes radius, rmass
+
+The atom style defined does not have these attributes.
+
+E: Pair granular requires ghost atoms store velocity
+
+Use the comm_modify vel yes command to enable this.
+
+E: Could not find pair fix neigh history ID
+
+UNDOCUMENTED
+
+U: Pair granular with shear history requires newton pair off
+
+This is a current restriction of the implementation of pair
+granular styles with history.
+
+U: Could not find pair fix ID
+
+A fix is created internally by the pair style to store shear
+history information.  You cannot delete it.
+
+*/