Added generalized granular option to fix wall/gran and fix wall/gran/region; some minor bug fixes for pair granular

doc/src/pair_granular.txt

@@ -533,14 +533,16 @@ the pairwise interaction force.  However, the single() function also
 calculates 10 extra pairwise quantities.  The first 3 are the
 components of the tangential force between particles I and J, acting
 on particle I.  The 4th is the magnitude of this tangential force.
-The next 3 (5-7) are the components of the relative velocity in the
+The next 3 (5-7) are the components of the rolling torque acting on 
-normal direction (along the line joining the 2 sphere centers).  The
+particle I. The next entry (8) is the magnitude of the rolling torque.
-last 3 (8-10) the components of the relative velocity in the
+The next entry (9) is the magnitude of the twisting torque acting
-tangential direction.
+about the vector connecting the two particle centers.
+The last 3 (10-12) are the components of the vector connecting
+the centers of the two particles (x_I - x_J).
 
 These extra quantities can be accessed by the "compute
 pair/local"_compute_pair_local.html command, as {p1}, {p2}, ...,
-{p10}.
+{p12}.
 
 :line
 
@@ -568,6 +570,7 @@ compute depend on atom velocities.  See the
 [Related commands:]
 
 "pair_coeff"_pair_coeff.html
+"pair gran/*"_pair_gran.html
 
 [Default:] 
 

src/GRANULAR/fix_wall_gran.h

@@ -54,12 +54,11 @@ class FixWallGran : public Fix {
   void hertz_history(double, double, double, double, double *, double,
       double *, double *, double *, double *, double, double,
       double *, double *);
-  void dmt_rolling(double, double, double, double, double *, double,
+  void granular(double, double, double, double, double *, double,
-      double *, double *, double *, double *, double, double,
+        double *, double *, double *, double *, double, double,
-      double *, double *);
+        double *, double *);
- // void jkr_rolling(double, double, double, double, double *, double,
+
- //     double *, double *, double *, double *, double, double,
+  double pulloff_distance(double);
- //     double *, double *);
 
  protected:
   int wallstyle,wiggle,wshear,axis;
@@ -67,23 +66,43 @@ class FixWallGran : public Fix {
   bigint time_origin;
   double kn,kt,gamman,gammat,xmu;
 
-  //For DMT/ROLLING
+  //For granular
-  int normaldamp, rollingdamp;
+  //Model choices
-  double Emod, Gmod, alpha, Ecoh, kR, muR, etaR;
+  int normal_model, damping_model;
+  int tangential_model, roll_model, twist_model;
 
+  int beyond_contact;
+
+  //History flags
+  int normal_history, tangential_history, roll_history, twist_history;
+
+  //Indices of history entries
+  int normal_history_index;
+  int tangential_history_index;
+  int roll_history_index;
+  int twist_history_index;
+
+  //Material coefficients
+  double Emod, poiss, Gmod;
+
+  //Contact model coefficients
+  double normal_coeffs[4];
+  double tangential_coeffs[3];
+  double roll_coeffs[3];
+  double twist_coeffs[3];
 
   double lo,hi,cylradius;
   double amplitude,period,omega,vshear;
   double dt;
   char *idregion;
 
-  int history;       // if particle/wall interaction stores history
+  int use_history;       // if particle/wall interaction stores history
-  int shearupdate;   // flag for whether shear history is updated
+  int history_update;   // flag for whether shear history is updated
-  int sheardim;      // # of shear history values per contact
+  int size_history;      // # of shear history values per contact
 
   // shear history for single contact per particle
 
-  double **shearone;
+  double **history_one;
 
   // rigid body masses for use in granular interactions
 

src/GRANULAR/fix_wall_gran_region.cpp

@@ -39,7 +39,8 @@ using namespace MathConst;
 
 // same as FixWallGran
 
-enum{HOOKE,HOOKE_HISTORY,HERTZ_HISTORY,JKR_ROLLING,DMT_ROLLING};
+enum{HOOKE,HOOKE_HISTORY,HERTZ_HISTORY,GRANULAR};
+enum {NORMAL_HOOKE, NORMAL_HERTZ, HERTZ_MATERIAL, DMT, JKR};
 
 #define BIG 1.0e20
 
@@ -47,7 +48,7 @@ enum{HOOKE,HOOKE_HISTORY,HERTZ_HISTORY,JKR_ROLLING,DMT_ROLLING};
 
 FixWallGranRegion::FixWallGranRegion(LAMMPS *lmp, int narg, char **arg) :
       FixWallGran(lmp, narg, arg), region(NULL), region_style(NULL), ncontact(NULL),
-      walls(NULL), shearmany(NULL), c2r(NULL)
+      walls(NULL), history_many(NULL), c2r(NULL)
 {
   restart_global = 1;
   motion_resetflag = 0;
@@ -66,17 +67,17 @@ FixWallGranRegion::FixWallGranRegion(LAMMPS *lmp, int narg, char **arg) :
   // re-allocate atom-based arrays with nshear
   // do not register with Atom class, since parent class did that
 
-  memory->destroy(shearone);
+  memory->destroy(history_one);
-  shearone = NULL;
+  history_one = NULL;
 
   ncontact = NULL;
   walls = NULL;
-  shearmany = NULL;
+  history_many = NULL;
   grow_arrays(atom->nmax);
 
   // initialize shear history as if particle is not touching region
 
-  if (history) {
+  if (use_history) {
     int nlocal = atom->nlocal;
     for (int i = 0; i < nlocal; i++)
       ncontact[i] = 0;
@@ -92,7 +93,7 @@ FixWallGranRegion::~FixWallGranRegion()
 
   memory->destroy(ncontact);
   memory->destroy(walls);
-  memory->destroy(shearmany);
+  memory->destroy(history_many);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -138,8 +139,8 @@ void FixWallGranRegion::post_force(int /*vflag*/)
 
   // do not update shear history during setup
 
-  shearupdate = 1;
+  history_update = 1;
-  if (update->setupflag) shearupdate = 0;
+  if (update->setupflag) history_update = 0;
 
   // if just reneighbored:
   // update rigid body masses for owned atoms if using FixRigid
@@ -188,7 +189,12 @@ void FixWallGranRegion::post_force(int /*vflag*/)
     if (mask[i] & groupbit) {
       if (!region->match(x[i][0],x[i][1],x[i][2])) continue;
 
-      nc = region->surface(x[i][0],x[i][1],x[i][2],radius[i]);
+      if (pairstyle == GRANULAR && normal_model == JKR){
+        nc = region->surface(x[i][0],x[i][1],x[i][2],radius[i]+pulloff_distance(radius[i]));
+      }
+      else{
+        nc = region->surface(x[i][0],x[i][1],x[i][2],radius[i]);
+      }
       if (nc > tmax)
         error->one(FLERR,"Too many wall/gran/region contacts for one particle");
 
@@ -198,7 +204,7 @@ void FixWallGranRegion::post_force(int /*vflag*/)
       // also set c2r[] = indices into region->contact[] for each of N contacts
       // process zero or one contact here, otherwise invoke update_contacts()
 
-      if (history) {
+      if (use_history) {
         if (nc == 0) {
           ncontact[i] = 0;
           continue;
@@ -209,8 +215,8 @@ void FixWallGranRegion::post_force(int /*vflag*/)
           if (ncontact[i] == 0) {
             ncontact[i] = 1;
             walls[i][0] = iwall;
-            for (m = 0; m < sheardim; m++)
+            for (m = 0; m < size_history; m++)
-              shearmany[i][0][m] = 0.0;
+              history_many[i][0][m] = 0.0;
           } else if (ncontact[i] > 1 || iwall != walls[i][0])
             update_contacts(i,nc);
         } else update_contacts(i,nc);
@@ -224,13 +230,20 @@ void FixWallGranRegion::post_force(int /*vflag*/)
 
         rsq = region->contact[ic].r*region->contact[ic].r;
 
+        if (pairstyle == GRANULAR && normal_model == JKR){
+          if (history_many[i][c2r[ic]][0] == 0.0 && rsq > radius[i]*radius[i]){
+            for (m = 0; m < size_history; m++)
+              history_many[i][0][m] = 0.0;
+            continue;
+          }
+        }
+
         dx = region->contact[ic].delx;
         dy = region->contact[ic].dely;
         dz = region->contact[ic].delz;
 
         if (regiondynamic) region->velocity_contact(vwall, x[i], ic);
 
-
         // meff = effective mass of sphere
         // if I is part of rigid body, use body mass
 
@@ -259,13 +272,13 @@ void FixWallGranRegion::post_force(int /*vflag*/)
         else if (pairstyle == HOOKE_HISTORY)
           hooke_history(rsq,dx,dy,dz,vwall,v[i],f[i],
               omega[i],torque[i],radius[i],meff,
-              shearmany[i][c2r[ic]], contact);
+              history_many[i][c2r[ic]], contact);
         else if (pairstyle == HERTZ_HISTORY)
           hertz_history(rsq,dx,dy,dz,vwall,region->contact[ic].radius,
               v[i],f[i],omega[i],torque[i],
-              radius[i],meff,shearmany[i][c2r[ic]], contact);
+              radius[i],meff,history_many[i][c2r[ic]], contact);
-        else if (pairstyle == DMT_ROLLING)
+        else if (pairstyle == GRANULAR)
-          dmt_rolling(rsq,dx,dy,dz,vwall,region->contact[ic].radius, v[i],f[i],omega[i],torque[i], radius[i],meff,shearmany[i][c2r[ic]], contact);
+          granular(rsq,dx,dy,dz,vwall,region->contact[ic].radius, v[i],f[i],omega[i],torque[i], radius[i],meff,history_many[i][c2r[ic]], contact);
 
       }
     }
@@ -294,8 +307,8 @@ void FixWallGranRegion::update_contacts(int i, int nc)
       if (region->contact[m].iwall == walls[i][iold]) break;
     if (m >= nc) {
       ilast = ncontact[i]-1;
-      for (j = 0; j < sheardim; j++)
+      for (j = 0; j < size_history; j++)
-        shearmany[i][iold][j] = shearmany[i][ilast][j];
+        history_many[i][iold][j] = history_many[i][ilast][j];
       walls[i][iold] = walls[i][ilast];
       ncontact[i]--;
     } else iold++;
@@ -317,8 +330,8 @@ void FixWallGranRegion::update_contacts(int i, int nc)
       iadd = ncontact[i];
 
       c2r[iadd] = inew;
-      for (j = 0; j < sheardim; j++)
+      for (j = 0; j < size_history; j++)
-        shearmany[i][iadd][j] = 0.0;
+        history_many[i][iadd][j] = 0.0;
       walls[i][iadd] = iwall;
       ncontact[i]++;
     }
@@ -333,10 +346,10 @@ double FixWallGranRegion::memory_usage()
 {
   int nmax = atom->nmax;
   double bytes = 0.0;
-  if (history) {                                   // shear history
+  if (use_history) {                                   // shear history
     bytes += nmax * sizeof(int);                   // ncontact
     bytes += nmax*tmax * sizeof(int);              // walls
-    bytes += nmax*tmax*sheardim * sizeof(double);  // shearmany
+    bytes += nmax*tmax*size_history * sizeof(double);  // history_many
   }
   if (fix_rigid) bytes += nmax * sizeof(int);      // mass_rigid
   return bytes;
@@ -348,10 +361,10 @@ double FixWallGranRegion::memory_usage()
 
 void FixWallGranRegion::grow_arrays(int nmax)
 {
-  if (history) {
+  if (use_history) {
     memory->grow(ncontact,nmax,"fix_wall_gran:ncontact");
     memory->grow(walls,nmax,tmax,"fix_wall_gran:walls");
-    memory->grow(shearmany,nmax,tmax,sheardim,"fix_wall_gran:shearmany");
+    memory->grow(history_many,nmax,tmax,size_history,"fix_wall_gran:history_many");
   }
   if (peratom_flag){
     memory->grow(array_atom,nmax,size_peratom_cols,"fix_wall_gran:array_atom");
@@ -366,12 +379,12 @@ void FixWallGranRegion::copy_arrays(int i, int j, int /*delflag*/)
 {
   int m,n,iwall;
 
-  if (history){
+  if (use_history){
     n = ncontact[i];
     for (iwall = 0; iwall < n; iwall++) {
       walls[j][iwall] = walls[i][iwall];
-      for (m = 0; m < sheardim; m++)
+      for (m = 0; m < size_history; m++)
-        shearmany[j][iwall][m] = shearmany[i][iwall][m];
+        history_many[j][iwall][m] = history_many[i][iwall][m];
     }
     ncontact[j] = ncontact[i];
   }
@@ -388,7 +401,7 @@ void FixWallGranRegion::copy_arrays(int i, int j, int /*delflag*/)
 
 void FixWallGranRegion::set_arrays(int i)
 {
-  if (history)
+  if (use_history)
     ncontact[i] = 0;
   if (peratom_flag){
     for (int m = 0; m < size_peratom_cols; m++)
@@ -405,13 +418,13 @@ int FixWallGranRegion::pack_exchange(int i, double *buf)
   int m;
 
   int n = 0;
-  if (history){
+  if (use_history){
     int count = ncontact[i];
     buf[n++] = ubuf(count).d;
     for (int iwall = 0; iwall < count; iwall++) {
       buf[n++] = ubuf(walls[i][iwall]).d;
-      for (m = 0; m < sheardim; m++)
+      for (m = 0; m < size_history; m++)
-        buf[n++] = shearmany[i][iwall][m];
+        buf[n++] = history_many[i][iwall][m];
     }
   }
   if (peratom_flag){
@@ -432,12 +445,12 @@ int FixWallGranRegion::unpack_exchange(int nlocal, double *buf)
 
 
   int n = 0;
-  if (history){
+  if (use_history){
     int count = ncontact[nlocal] = (int) ubuf(buf[n++]).i;
     for (int iwall = 0; iwall < count; iwall++) {
       walls[nlocal][iwall] = (int) ubuf(buf[n++]).i;
-      for (m = 0; m < sheardim; m++)
+      for (m = 0; m < size_history; m++)
-        shearmany[nlocal][iwall][m] = buf[n++];
+        history_many[nlocal][iwall][m] = buf[n++];
     }
   }
   if (peratom_flag){
@@ -456,7 +469,7 @@ int FixWallGranRegion::pack_restart(int i, double *buf)
 {
   int m;
 
-  if (!history) return 0;
+  if (!use_history) return 0;
 
   int n = 1;
   int count = ncontact[i];
@@ -464,8 +477,8 @@ int FixWallGranRegion::pack_restart(int i, double *buf)
   buf[n++] = ubuf(count).d;
   for (int iwall = 0; iwall < count; iwall++) {
     buf[n++] = ubuf(walls[i][iwall]).d;
-    for (m = 0; m < sheardim; m++)
+    for (m = 0; m < size_history; m++)
-      buf[n++] = shearmany[i][iwall][m];
+      buf[n++] = history_many[i][iwall][m];
   }
   buf[0] = n;
   return n;
@@ -479,7 +492,7 @@ void FixWallGranRegion::unpack_restart(int nlocal, int nth)
 {
   int k;
 
-  if (!history) return;
+  if (!use_history) return;
 
   double **extra = atom->extra;
 
@@ -492,8 +505,8 @@ void FixWallGranRegion::unpack_restart(int nlocal, int nth)
   int count = ncontact[nlocal] = (int) ubuf(extra[nlocal][m++]).i;
   for (int iwall = 0; iwall < count; iwall++) {
     walls[nlocal][iwall] = (int) ubuf(extra[nlocal][m++]).i;
-    for (k = 0; k < sheardim; k++)
+    for (k = 0; k < size_history; k++)
-      shearmany[nlocal][iwall][k] = extra[nlocal][m++];
+      history_many[nlocal][iwall][k] = extra[nlocal][m++];
   }
 }
 
@@ -503,8 +516,8 @@ void FixWallGranRegion::unpack_restart(int nlocal, int nth)
 
 int FixWallGranRegion::maxsize_restart()
 {
-  if (!history) return 0;
+  if (!use_history) return 0;
-  return 2 + tmax*(sheardim+1);
+  return 2 + tmax*(size_history+1);
 }
 
 /* ----------------------------------------------------------------------
@@ -513,8 +526,8 @@ int FixWallGranRegion::maxsize_restart()
 
 int FixWallGranRegion::size_restart(int nlocal)
 {
-  if (!history) return 0;
+  if (!use_history) return 0;
-  return 2 + ncontact[nlocal]*(sheardim+1);
+  return 2 + ncontact[nlocal]*(size_history+1);
 }
 
 /* ----------------------------------------------------------------------

src/GRANULAR/fix_wall_gran_region.h

@@ -54,7 +54,7 @@ class FixWallGranRegion : public FixWallGran {
   int tmax;              // max # of region walls one particle can touch
   int *ncontact;         // # of shear contacts per particle
   int **walls;           // which wall each contact is with
-  double ***shearmany;   // shear history per particle per contact
+  double ***history_many;   // history per particle per contact
   int *c2r;              // contact to region mapping
                          // c2r[i] = index of Ith contact in
                          //   region-contact[] list of contacts

src/GRANULAR/pair_granular.cpp

@@ -45,6 +45,7 @@ using namespace MathConst;
 #define SIXROOT6 14.69693845669906728801   // 6*sqrt(6)
 #define INVROOT6 0.40824829046386307274    // 1/sqrt(6)
 #define FOURTHIRDS 1.333333333333333       // 4/3
+#define THREEQUARTERS 0.75                 // 3/4
 #define TWOPI 6.28318530717959             // 2*PI
 
 #define EPSILON 1e-10
@@ -63,7 +64,7 @@ PairGranular::PairGranular(LAMMPS *lmp) : Pair(lmp)
   no_virial_fdotr_compute = 1;
   fix_history = NULL;
 
-  single_extra = 9;
+  single_extra = 12;
   svector = new double[single_extra];
 
   neighprev = 0;
@@ -238,10 +239,11 @@ void PairGranular::compute(int eflag, int vflag)
       radsum = radi + radj;
 
       E = normal_coeffs[itype][jtype][0];
-      Reff = radi*radj/(radi+radj);
+      Reff = radi*radj/radsum;
       touchflag = false;
 
       if (normal_model[itype][jtype] == JKR){
+        E *= THREEQUARTERS;
         if (touch[jj]){
           R2 = Reff*Reff;
           coh = normal_coeffs[itype][jtype][3];
@@ -319,17 +321,17 @@ void PairGranular::compute(int eflag, int vflag)
           sqrt3 = MAX(0, 4*dR - t5 + SIXROOT6*coh*M_PI*R2/(E*t6));
           a = INVROOT6*(t6 + sqrt(sqrt3));
           a2 = a*a;
-          knfac = FOURTHIRDS*E*a;
+          knfac = normal_coeffs[itype][jtype][0]*a;
           Fne = knfac*a2/Reff - TWOPI*a2*sqrt(4*coh*E/(M_PI*a));
         }
         else{
           knfac = E; //Hooke
+          Fne = knfac*delta;
           a = sqrt(dR);
           if (normal_model[itype][jtype] != HOOKE){
             Fne *= a;
             knfac *= a;
           }
-          Fne = knfac*delta;
           if (normal_model[itype][jtype] == DMT)
             Fne -= 4*MY_PI*normal_coeffs[itype][jtype][3]*Reff;
         }
@@ -711,9 +713,9 @@ void PairGranular::coeff(int narg, char **arg)
     }
     else if (strcmp(arg[iarg], "hertz/material") == 0){
       int num_coeffs = 3;
-      if (iarg + num_coeffs >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for Hertz option");
+      if (iarg + num_coeffs >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for Hertz/material option");
-      normal_model_one = HERTZ;
+      normal_model_one = HERTZ_MATERIAL;
-      normal_coeffs_one[0] = force->numeric(FLERR,arg[iarg+1])*FOURTHIRDS; //E
+      normal_coeffs_one[0] = force->numeric(FLERR,arg[iarg+1]); //E
       normal_coeffs_one[1] = force->numeric(FLERR,arg[iarg+2]); //damping
       normal_coeffs_one[2] = force->numeric(FLERR,arg[iarg+3]); //Poisson's ratio
       iarg += num_coeffs+1;
@@ -721,10 +723,10 @@ void PairGranular::coeff(int narg, char **arg)
     else if (strcmp(arg[iarg], "dmt") == 0){
       if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for Hertz option");
       normal_model_one = DMT;
-      normal_coeffs_one[0] = force->numeric(FLERR,arg[iarg+1])*FOURTHIRDS; //E
+      normal_coeffs_one[0] = force->numeric(FLERR,arg[iarg+1]); //E
       normal_coeffs_one[1] = force->numeric(FLERR,arg[iarg+2]); //damping
       normal_coeffs_one[2] = force->numeric(FLERR,arg[iarg+3]); //Poisson's ratio
-      normal_coeffs_one[3] = force->numeric(FLERR,arg[iarg+3]); //cohesion
+      normal_coeffs_one[3] = force->numeric(FLERR,arg[iarg+4]); //cohesion
       iarg += 5;
     }
     else if (strcmp(arg[iarg], "jkr") == 0){
@@ -755,21 +757,27 @@ void PairGranular::coeff(int narg, char **arg)
       iarg += 1;
     }
     else if (strcmp(arg[iarg], "tangential") == 0){
-      if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for tangential model");
+      if (iarg + 1 >= narg) error->all(FLERR,"Illegal pair_coeff command, must specify tangential model after 'tangential' keyword");
       if (strcmp(arg[iarg+1], "linear_nohistory") == 0){
+        if (iarg + 3 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for tangential model");
         tangential_model_one = TANGENTIAL_NOHISTORY;
+        tangential_coeffs_one[0] = 0;
+        tangential_coeffs_one[1] = force->numeric(FLERR,arg[iarg+2]); //gammat
+        tangential_coeffs_one[2] = force->numeric(FLERR,arg[iarg+3]); //friction coeff.
+        iarg += 4;
       }
       else if (strcmp(arg[iarg+1], "linear_history") == 0){
+        if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for tangential model");
         tangential_model_one = TANGENTIAL_HISTORY;
         tangential_history = 1;
+        tangential_coeffs_one[0] = force->numeric(FLERR,arg[iarg+2]); //kt
+        tangential_coeffs_one[1] = force->numeric(FLERR,arg[iarg+3]); //gammat
+        tangential_coeffs_one[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff.
+        iarg += 5;
       }
       else{
         error->all(FLERR, "Illegal pair_coeff command, tangential model not recognized");
       }
-      tangential_coeffs_one[0] = force->numeric(FLERR,arg[iarg+2]); //kt
-      tangential_coeffs_one[1] = force->numeric(FLERR,arg[iarg+3]); //gammat
-      tangential_coeffs_one[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff.
-      iarg += 5;
     }
     else if (strcmp(arg[iarg], "rolling") == 0){
       if (iarg + 1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters");
@@ -842,7 +850,7 @@ void PairGranular::coeff(int narg, char **arg)
       if (normal_model_one != HERTZ && normal_model_one != HOOKE){
         Emod[i][j] = Emod[j][i] = normal_coeffs_one[0];
         poiss[i][j] = poiss[j][i] = normal_coeffs_one[2];
-        normal_coeffs[i][j][0] = normal_coeffs[j][i][0] = mix_stiffnessE(Emod[i][j], Emod[i][j], poiss[i][j], poiss[i][j]);
+        normal_coeffs[i][j][0] = normal_coeffs[j][i][0] = FOURTHIRDS*mix_stiffnessE(Emod[i][j], Emod[i][j], poiss[i][j], poiss[i][j]);
       }
       else{
         normal_coeffs[i][j][0] = normal_coeffs[j][i][0] = normal_coeffs_one[0];
@@ -1068,14 +1076,13 @@ double PairGranular::init_one(int i, int j)
     if (((maxrad_dynamic[i] > 0.0) && (maxrad_dynamic[j] > 0.0)) ||
         ((maxrad_dynamic[i] > 0.0) &&  (maxrad_frozen[j] > 0.0)) ||
         ((maxrad_frozen[i] > 0.0)  && (maxrad_dynamic[j] > 0.0))) { // radius info about both i and j exist
+
       cutoff = maxrad_dynamic[i]+maxrad_dynamic[j];
+      pulloff = 0.0;
       if (normal_model[i][j] == JKR){
         pulloff = pulloff_distance(maxrad_dynamic[i], maxrad_dynamic[j], i, j);
         cutoff += pulloff;
       }
-      else{
-        pulloff = 0;
-      }
 
       if (normal_model[i][j] == JKR)
         pulloff = pulloff_distance(maxrad_frozen[i], maxrad_dynamic[j], i, j);
@@ -1224,10 +1231,12 @@ double PairGranular::single(int i, int j, int itype, int jtype,
   radi = radius[i];
   radj = radius[j];
   radsum = radi + radj;
-  Reff = radi*radj/(radi+radj);
+  Reff = radi*radj/radsum;
 
   bool touchflag;
+  E = normal_coeffs[itype][jtype][0];
   if (normal_model[itype][jtype] == JKR){
+    E *= THREEQUARTERS;
     R2 = Reff*Reff;
     coh = normal_coeffs[itype][jtype][3];
     a = cbrt(9.0*M_PI*coh*R2/(4*E));
@@ -1333,7 +1342,7 @@ double PairGranular::single(int i, int j, int itype, int jtype,
     sqrt3 = MAX(0, 4*dR - t5 + SIXROOT6*coh*M_PI*R2/(E*t6));
     a = INVROOT6*(t6 + sqrt(sqrt3));
     a2 = a*a;
-    knfac = FOURTHIRDS*E*a;
+    knfac = normal_coeffs[itype][jtype][0]*a;
     Fne = knfac*a2/Reff - TWOPI*a2*sqrt(4*coh*E/(M_PI*a));
   }
   else{
@@ -1348,7 +1357,6 @@ double PairGranular::single(int i, int j, int itype, int jtype,
       Fne -= 4*MY_PI*normal_coeffs[itype][jtype][3]*Reff;
   }
 
-  //Consider restricting Hooke to only have 'velocity' as an option for damping?
   if (damping_model[itype][jtype] == VELOCITY){
     damp_normal = normal_coeffs[itype][jtype][1];
   }
@@ -1536,6 +1544,9 @@ double PairGranular::single(int i, int j, int itype, int jtype,
   svector[6] = fr3;
   svector[7] = fr;
   svector[8] = magtortwist;
+  svector[9] = delx;
+  svector[10] = dely;
+  svector[11] = delz;
   return 0.0;
 }
 
@@ -1614,7 +1625,7 @@ double PairGranular::pulloff_distance(double radi, double radj, int itype, int j
   Reff = radi*radj/(radi+radj);
   if (Reff <= 0) return 0;
   coh = normal_coeffs[itype][itype][3];
-  E = normal_coeffs[itype][jtype][0];
+  E = normal_coeffs[itype][jtype][0]*THREEQUARTERS;
   a = cbrt(9*M_PI*coh*Reff/(4*E));
   return a*a/Reff - 2*sqrt(M_PI*coh*a/E);
 }

src/GRANULAR/pair_granular.h

@@ -65,9 +65,9 @@ public:
 private:
   int size_history;
 
-  //Models choices
+  //Model choices
   int **normal_model, **damping_model;
-  double **tangential_model, **roll_model, **twist_model;
+  int **tangential_model, **roll_model, **twist_model;
 
   //History flags
   int normal_history, tangential_history, roll_history, twist_history;