git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@5881 f3b2605a-c512-4ea7-a41b-209d697bcdaa

src/GRANULAR/fix_wall_gran.cpp

@@ -182,6 +182,7 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
     shear[i][0] = shear[i][1] = shear[i][2] = 0.0;
 
   time_origin = update->ntimestep;
+  laststep = -1;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -282,6 +283,9 @@ void FixWallGran::post_force(int vflag)
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
+  if (update->ntimestep > laststep) shearupdate = 1;
+  else shearupdate = 0;
+
   for (int i = 0; i < nlocal; i++) {
     if (mask[i] & groupbit) {
 
@@ -338,6 +342,8 @@ void FixWallGran::post_force(int vflag)
       }
     }
   }
+
+  laststep = update->ntimestep;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -489,18 +495,22 @@ void FixWallGran::hooke_history(double rsq, double dx, double dy, double dz,
 
   // shear history effects
 
-  shear[0] += vtr1*dt;
-  shear[1] += vtr2*dt;
-  shear[2] += vtr3*dt;
+  if (shearupdate) {
+    shear[0] += vtr1*dt;
+    shear[1] += vtr2*dt;
+    shear[2] += vtr3*dt;
+  }
   shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + shear[2]*shear[2]);
 
   // rotate shear displacements
 
   rsht = shear[0]*dx + shear[1]*dy + shear[2]*dz;
   rsht = rsht*rsqinv;
-  shear[0] -= rsht*dx;
-  shear[1] -= rsht*dy;
-  shear[2] -= rsht*dz;
+  if (shearupdate) {
+    shear[0] -= rsht*dx;
+    shear[1] -= rsht*dy;
+    shear[2] -= rsht*dz;
+  }
 
   // tangential forces = shear + tangential velocity damping
 

src/GRANULAR/fix_wall_gran.h

@@ -55,8 +55,10 @@ class FixWallGran : public Fix {
   int nlevels_respa;
   int time_origin;
 
+  bigint laststep;
   int *touch;
   double **shear;
+  int shearupdate;
 
   void hooke(double, double, double, double, double *,
 	     double *, double *, double *, double *, double, double);

src/GRANULAR/pair_gran_hertz_history.cpp

@@ -21,6 +21,7 @@
 #include "string.h"
 #include "pair_gran_hertz_history.h"
 #include "atom.h"
+#include "update.h"
 #include "force.h"
 #include "neigh_list.h"
 #include "error.h"
@@ -55,6 +56,9 @@ void PairGranHertzHistory::compute(int eflag, int vflag)
   if (eflag || vflag) ev_setup(eflag,vflag);
   else evflag = vflag_fdotr = 0;
 
+  if (eflag || vflag) ev_setup(eflag,vflag);
+  else evflag = vflag_fdotr = 0;
+
   double **x = atom->x;
   double **v = atom->v;
   double **f = atom->f;
@@ -168,9 +172,11 @@ void PairGranHertzHistory::compute(int eflag, int vflag)
 
 	touch[jj] = 1;
 	shear = &allshear[3*jj];
-        shear[0] += vtr1*dt;
-        shear[1] += vtr2*dt;
-        shear[2] += vtr3*dt;
+	if (shearupdate) {
+	  shear[0] += vtr1*dt;
+	  shear[1] += vtr2*dt;
+	  shear[2] += vtr3*dt;
+	}
         shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + 
 		      shear[2]*shear[2]);
 
@@ -178,9 +184,11 @@ void PairGranHertzHistory::compute(int eflag, int vflag)
 
 	rsht = shear[0]*delx + shear[1]*dely + shear[2]*delz;
 	rsht *= rsqinv;
-        shear[0] -= rsht*delx;
-        shear[1] -= rsht*dely;
-        shear[2] -= rsht*delz;
+	if (shearupdate) {
+	  shear[0] -= rsht*delx;
+	  shear[1] -= rsht*dely;
+	  shear[2] -= rsht*delz;
+	}
 
 	// tangential forces = shear + tangential velocity damping
 
@@ -237,6 +245,8 @@ void PairGranHertzHistory::compute(int eflag, int vflag)
       }
     }
   }
+
+  laststep = update->ntimestep;
 }
 
 /* ----------------------------------------------------------------------

src/GRANULAR/pair_gran_hooke_history.cpp

@@ -49,6 +49,8 @@ PairGranHookeHistory::PairGranHookeHistory(LAMMPS *lmp) : Pair(lmp)
   no_virial_compute = 1;
   history = 1;
   fix_history = NULL;
+
+  laststep = -1;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -88,6 +90,9 @@ void PairGranHookeHistory::compute(int eflag, int vflag)
   if (eflag || vflag) ev_setup(eflag,vflag);
   else evflag = vflag_fdotr = 0;
 
+  int shearupdate = 0;
+  if (update->ntimestep > laststep) shearupdate = 1;
+
   double **x = atom->x;
   double **v = atom->v;
   double **f = atom->f;
@@ -199,9 +204,12 @@ void PairGranHookeHistory::compute(int eflag, int vflag)
 
 	touch[jj] = 1;
 	shear = &allshear[3*jj];
-        shear[0] += vtr1*dt;
-        shear[1] += vtr2*dt;
-        shear[2] += vtr3*dt;
+
+	if (shearupdate) {
+	  shear[0] += vtr1*dt;
+	  shear[1] += vtr2*dt;
+	  shear[2] += vtr3*dt;
+	}
         shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] +
 		      shear[2]*shear[2]);
 
@@ -209,9 +217,11 @@ void PairGranHookeHistory::compute(int eflag, int vflag)
 
 	rsht = shear[0]*delx + shear[1]*dely + shear[2]*delz;
 	rsht *= rsqinv;
-        shear[0] -= rsht*delx;
-        shear[1] -= rsht*dely;
-        shear[2] -= rsht*delz;
+	if (shearupdate) {
+	  shear[0] -= rsht*delx;
+	  shear[1] -= rsht*dely;
+	  shear[2] -= rsht*delz;
+	}
 
 	// tangential forces = shear + tangential velocity damping
 
@@ -268,6 +278,8 @@ void PairGranHookeHistory::compute(int eflag, int vflag)
       }
     }
   }
+
+  laststep = update->ntimestep;
 }
 
 /* ----------------------------------------------------------------------

src/GRANULAR/pair_gran_hooke_history.h

@@ -46,7 +46,10 @@ class PairGranHookeHistory : public Pair {
   double dt;
   int freeze_group_bit;
   int history;
+
+  bigint laststep;
   class FixShearHistory *fix_history;
+  int shearupdate;
 
   double *onerad_dynamic,*onerad_frozen;
   double *maxrad_dynamic,*maxrad_frozen;

src/fix.h

@@ -92,6 +92,7 @@ class Fix : protected Pointers {
   virtual void init() {}
   virtual void init_list(int, class NeighList *) {}
   virtual void setup(int) {}
+  virtual void setup_pre_exchange() {}
   virtual void setup_pre_force(int) {}
   virtual void min_setup(int) {}
   virtual void initial_integrate(int) {}

src/fix_shear_history.cpp

@@ -33,8 +33,12 @@ using namespace LAMMPS_NS;
 FixShearHistory::FixShearHistory(LAMMPS *lmp, int narg, char **arg) :
   Fix(lmp, narg, arg)
 {
+  // set time_depend so that history will be preserved correctly
+  // across multiple runs via laststep setting in granular pair styles
+
   restart_peratom = 1;
   create_attribute = 1;
+  time_depend = 1;
 
   // perform initial allocation of atom-based arrays
   // register with atom class
@@ -85,6 +89,13 @@ void FixShearHistory::init()
     error->all("Pair style granular with history requires atoms have IDs");
 }
 
+/* ---------------------------------------------------------------------- */
+
+void FixShearHistory::setup_pre_exchange()
+{
+  pre_exchange();
+}
+
 /* ----------------------------------------------------------------------
    copy shear partner info from neighbor lists to atom arrays
    so can be exchanged with atoms

src/fix_shear_history.h

@@ -34,6 +34,7 @@ class FixShearHistory : public Fix {
   ~FixShearHistory();
   int setmask();
   void init();
+  void setup_pre_exchange();
   void pre_exchange();
 
   double memory_usage();

src/modify.cpp

@@ -255,7 +255,17 @@ void Modify::setup(int vflag)
 }
 
 /* ----------------------------------------------------------------------
-   setup pre_force call, only for relevant fixes
+   setup pre_exchange call, only for fixes that define pre_exchange
+------------------------------------------------------------------------- */
+
+void Modify::setup_pre_exchange()
+{
+  for (int i = 0; i < n_pre_exchange; i++)
+    fix[list_pre_exchange[i]]->setup_pre_exchange();
+}
+
+/* ----------------------------------------------------------------------
+   setup pre_force call, only for fixes that define pre_force
 ------------------------------------------------------------------------- */
 
 void Modify::setup_pre_force(int vflag)

src/modify.h

@@ -43,6 +43,7 @@ class Modify : protected Pointers {
   ~Modify();
   void init();
   void setup(int);
+  void setup_pre_exchange();
   void setup_pre_force(int);
   void initial_integrate(int);
   void post_integrate();

src/respa.cpp

@@ -327,6 +327,7 @@ void Respa::setup()
   // build neighbor lists
 
   atom->setup();
+  modify->setup_pre_exchange();
   if (triclinic) domain->x2lamda(atom->nlocal);
   domain->pbc();
   domain->reset_box();

src/verlet.cpp

@@ -90,6 +90,7 @@ void Verlet::setup()
   // build neighbor lists
 
   atom->setup();
+  modify->setup_pre_exchange();
   if (triclinic) domain->x2lamda(atom->nlocal);
   domain->pbc();
   domain->reset_box();