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

src/GRANULAR/pair_gran_hertz_history.cpp

@@ -276,3 +276,149 @@ void PairGranHertzHistory::settings(int narg, char **arg)
   kn /= force->nktv2p;
   kt /= force->nktv2p;
 }
+
+/* ---------------------------------------------------------------------- */
+
+double PairGranHertzHistory::single(int i, int j, int itype, int jtype,
+				    double rsq,
+				    double factor_coul, double factor_lj,
+				    double &fforce)
+{
+  double radi,radj,radsum;
+  double r,rinv,rsqinv,delx,dely,delz;
+  double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3,wr1,wr2,wr3;
+  double meff,damp,ccel,polyhertz;
+  double vtr1,vtr2,vtr3,vrel,shrmag,rsht;
+  double fs1,fs2,fs3,fs,fn;
+
+  double *radius = atom->radius;
+  radi = radius[i];
+  radj = radius[j];
+  radsum = radi + radj;
+
+  if (rsq >= radsum*radsum) {
+    fforce = 0.0;
+    svector[0] = svector[1] = svector[2] = svector[3] = 0.0;
+    return 0.0;
+  }
+
+  r = sqrt(rsq);
+  rinv = 1.0/r;
+  rsqinv = 1.0/rsq;
+
+  // relative translational velocity
+
+  double **v = atom->v;
+  vr1 = v[i][0] - v[j][0];
+  vr2 = v[i][1] - v[j][1];
+  vr3 = v[i][2] - v[j][2];
+
+  // normal component
+  
+  double **x = atom->x;
+  delx = x[i][0] - x[j][0];
+  dely = x[i][1] - x[j][1];
+  delz = x[i][2] - x[j][2];
+
+  vnnr = vr1*delx + vr2*dely + vr3*delz;
+  vn1 = delx*vnnr * rsqinv;
+  vn2 = dely*vnnr * rsqinv;
+  vn3 = delz*vnnr * rsqinv;
+  
+  // tangential component
+
+  vt1 = vr1 - vn1;
+  vt2 = vr2 - vn2;
+  vt3 = vr3 - vn3;
+
+  // relative rotational velocity
+
+  double **omega = atom->omega;
+  wr1 = (radi*omega[i][0] + radj*omega[j][0]) * rinv;
+  wr2 = (radi*omega[i][1] + radj*omega[j][1]) * rinv;
+  wr3 = (radi*omega[i][2] + radj*omega[j][2]) * rinv;
+  
+  // normal force = Hertzian contact + normal velocity damping
+
+  double *rmass = atom->rmass;
+  double *mass = atom->mass;
+  int *mask = atom->mask;
+
+  if (rmass) {
+    meff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]);
+    if (mask[i] & freeze_group_bit) meff = rmass[j];
+    if (mask[j] & freeze_group_bit) meff = rmass[i];
+  } else {
+    meff = mass[itype]*mass[jtype] / (mass[itype]+mass[jtype]);
+    if (mask[i] & freeze_group_bit) meff = mass[jtype];
+    if (mask[j] & freeze_group_bit) meff = mass[itype];
+  }
+  
+  damp = meff*gamman*vnnr*rsqinv;
+  ccel = kn*(radsum-r)*rinv - damp;
+  polyhertz = sqrt((radsum-r)*radi*radj / radsum);
+  ccel *= polyhertz;
+
+  // relative velocities
+
+  vtr1 = vt1 - (delz*wr2-dely*wr3);
+  vtr2 = vt2 - (delx*wr3-delz*wr1);
+  vtr3 = vt3 - (dely*wr1-delx*wr2);
+  vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3;
+  vrel = sqrt(vrel);
+  
+  // shear history effects
+  // neighprev = index of found neigh on previous call
+  // search entire jnum list of neighbors of I for neighbor J
+  // start from neighprev, since will typically be next neighbor
+  // reset neighprev to 0 as necessary
+
+  int *jlist = list->firstneigh[i];
+  int jnum = list->numneigh[i];
+  int *touch = list->listgranhistory->firstneigh[i];
+  double *allshear = list->listgranhistory->firstdouble[i];
+
+  for (int jj = 0; jj < jnum; jj++) {
+    neighprev++;
+    if (neighprev >= jnum) neighprev = 0;
+    if (touch[neighprev] == j) break;
+  }
+
+  double *shear = &allshear[3*neighprev];
+  shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + 
+		shear[2]*shear[2]);
+  
+  // rotate shear displacements
+  
+  rsht = shear[0]*delx + shear[1]*dely + shear[2]*delz;
+  rsht *= rsqinv;
+  
+  // tangential forces = shear + tangential velocity damping
+  
+  fs1 = -polyhertz * (kt*shear[0] + meff*gammat*vtr1);
+  fs2 = -polyhertz * (kt*shear[1] + meff*gammat*vtr2);
+  fs3 = -polyhertz * (kt*shear[2] + meff*gammat*vtr3);
+  
+  // rescale frictional displacements and forces if needed
+  
+  fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3);
+  fn = xmu * fabs(ccel*r);
+  
+  if (fs > fn) {
+    if (shrmag != 0.0) {
+      fs1 *= fn/fs;
+      fs2 *= fn/fs;
+      fs3 *= fn/fs;
+      fs *= fn/fs;
+    } else fs1 = fs2 = fs3 = fs = 0.0;
+  }
+  
+  // set all forces and return no energy
+  
+  fforce = ccel;
+  svector[0] = fs1;
+  svector[1] = fs2;
+  svector[2] = fs3;
+  svector[3] = fs;
+  return 0.0;
+}

src/GRANULAR/pair_gran_hertz_history.h

@@ -29,6 +29,7 @@ class PairGranHertzHistory : public PairGranHookeHistory {
   PairGranHertzHistory(class LAMMPS *);
   virtual void compute(int, int);
   void settings(int, char **);
+  double single(int, int, int, int, double, double, double, double &);
 };
 
 }

src/GRANULAR/pair_gran_hooke.cpp

@@ -191,3 +191,111 @@ void PairGranHooke::compute(int eflag, int vflag)
 
   if (vflag_fdotr) virial_fdotr_compute();
 }
+
+/* ---------------------------------------------------------------------- */
+
+double PairGranHooke::single(int i, int j, int itype, int jtype, double rsq,
+			     double factor_coul, double factor_lj,
+			     double &fforce)
+{
+  double radi,radj,radsum,r,rinv,rsqinv;
+  double delx,dely,delz;
+  double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3,wr1,wr2,wr3;
+  double vtr1,vtr2,vtr3,vrel;
+  double meff,damp,ccel;
+  double fn,fs,ft,fs1,fs2,fs3;
+
+  double *radius = atom->radius;
+  radi = radius[i];
+  radj = radius[j];
+  radsum = radi + radj;
+
+  // zero out forces if caller requests non-touching pair outside cutoff
+
+  if (rsq >= radsum*radsum) {
+    fforce = 0.0;
+    svector[0] = svector[1] = svector[2] = svector[3] = 0.0;
+    return 0.0;
+  }
+
+  r = sqrt(rsq);
+  rinv = 1.0/r;
+  rsqinv = 1.0/rsq;
+
+  // relative translational velocity
+
+  double **v = atom->v;
+  vr1 = v[i][0] - v[j][0];
+  vr2 = v[i][1] - v[j][1];
+  vr3 = v[i][2] - v[j][2];
+  
+  // normal component
+
+  double **x = atom->x;
+  delx = x[i][0] - x[j][0];
+  dely = x[i][1] - x[j][1];
+  delz = x[i][2] - x[j][2];
+
+  vnnr = vr1*delx + vr2*dely + vr3*delz;
+  vn1 = delx*vnnr * rsqinv;
+  vn2 = dely*vnnr * rsqinv;
+  vn3 = delz*vnnr * rsqinv;
+
+  // tangential component
+
+  vt1 = vr1 - vn1;
+  vt2 = vr2 - vn2;
+  vt3 = vr3 - vn3;
+
+  // relative rotational velocity
+
+  double **omega = atom->omega;
+  wr1 = (radi*omega[i][0] + radj*omega[j][0]) * rinv;
+  wr2 = (radi*omega[i][1] + radj*omega[j][1]) * rinv;
+  wr3 = (radi*omega[i][2] + radj*omega[j][2]) * rinv;
+
+  // normal forces = Hookian contact + normal velocity damping
+
+  double *rmass = atom->rmass;
+  double *mass = atom->mass;
+  int *mask = atom->mask;
+
+  if (rmass) {
+    meff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]);
+    if (mask[i] & freeze_group_bit) meff = rmass[j];
+    if (mask[j] & freeze_group_bit) meff = rmass[i];
+  } else {
+    meff = mass[itype]*mass[jtype] / (mass[itype]+mass[jtype]);
+    if (mask[i] & freeze_group_bit) meff = mass[jtype];
+    if (mask[j] & freeze_group_bit) meff = mass[itype];
+  }
+
+  damp = meff*gamman*vnnr*rsqinv;
+  ccel = kn*(radsum-r)*rinv - damp;
+
+  // relative velocities
+
+  vtr1 = vt1 - (delz*wr2-dely*wr3);
+  vtr2 = vt2 - (delx*wr3-delz*wr1);
+  vtr3 = vt3 - (dely*wr1-delx*wr2);
+  vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3;
+  vrel = sqrt(vrel);
+
+  // force normalization
+
+  fn = xmu * fabs(ccel*r);
+  fs = meff*gammat*vrel;
+  if (vrel != 0.0) ft = MIN(fn,fs) / vrel;
+  else ft = 0.0;
+  
+  // set all forces and return no energy
+  
+  fforce = ccel;
+  svector[0] = -ft*vtr1;
+  svector[1] = -ft*vtr2;
+  svector[2] = -ft*vtr3;
+  svector[3] = sqrt(svector[0]*svector[0] + 
+		    svector[1]*svector[1] + 
+		    svector[2]*svector[2]);
+  return 0.0;
+}

src/GRANULAR/pair_gran_hooke.h

@@ -28,6 +28,7 @@ class PairGranHooke : public PairGranHookeHistory {
  public:
   PairGranHooke(class LAMMPS *);
   virtual void compute(int, int);
+  double single(int, int, int, int, double, double, double, double &);
 };
 
 }

src/GRANULAR/pair_gran_hooke_history.cpp

@@ -42,19 +42,24 @@ using namespace LAMMPS_NS;
 
 PairGranHookeHistory::PairGranHookeHistory(LAMMPS *lmp) : Pair(lmp)
 {
-  single_enable = 0;
+  single_enable = 1;
   no_virial_fdotr_compute = 1;
   history = 1;
   fix_history = NULL;
   suffix = NULL;
 
+  single_extra = 4;
+  svector = new double[4];
+
   laststep = -1;
+  neighprev = 0;
 }
 
 /* ---------------------------------------------------------------------- */
 
 PairGranHookeHistory::~PairGranHookeHistory()
 {
+  delete [] svector;
   if (fix_history) modify->delete_fix("SHEAR_HISTORY");
   if (suffix) delete[] suffix;
 
@@ -545,3 +550,147 @@ void PairGranHookeHistory::reset_dt()
 {
   dt = update->dt;
 }
+
+/* ---------------------------------------------------------------------- */
+
+double PairGranHookeHistory::single(int i, int j, int itype, int jtype,
+				    double rsq,
+				    double factor_coul, double factor_lj,
+				    double &fforce)
+{
+  double radi,radj,radsum;
+  double r,rinv,rsqinv,delx,dely,delz;
+  double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3,wr1,wr2,wr3;
+  double meff,damp,ccel,polyhertz;
+  double vtr1,vtr2,vtr3,vrel,shrmag,rsht;
+  double fs1,fs2,fs3,fs,fn;
+
+  double *radius = atom->radius;
+  radi = radius[i];
+  radj = radius[j];
+  radsum = radi + radj;
+
+  if (rsq >= radsum*radsum) {
+    fforce = 0.0;
+    svector[0] = svector[1] = svector[2] = svector[3] = 0.0;
+    return 0.0;
+  }
+
+  r = sqrt(rsq);
+  rinv = 1.0/r;
+  rsqinv = 1.0/rsq;
+
+  // relative translational velocity
+
+  double **v = atom->v;
+  vr1 = v[i][0] - v[j][0];
+  vr2 = v[i][1] - v[j][1];
+  vr3 = v[i][2] - v[j][2];
+
+  // normal component
+  
+  double **x = atom->x;
+  delx = x[i][0] - x[j][0];
+  dely = x[i][1] - x[j][1];
+  delz = x[i][2] - x[j][2];
+
+  vnnr = vr1*delx + vr2*dely + vr3*delz;
+  vn1 = delx*vnnr * rsqinv;
+  vn2 = dely*vnnr * rsqinv;
+  vn3 = delz*vnnr * rsqinv;
+  
+  // tangential component
+
+  vt1 = vr1 - vn1;
+  vt2 = vr2 - vn2;
+  vt3 = vr3 - vn3;
+
+  // relative rotational velocity
+
+  double **omega = atom->omega;
+  wr1 = (radi*omega[i][0] + radj*omega[j][0]) * rinv;
+  wr2 = (radi*omega[i][1] + radj*omega[j][1]) * rinv;
+  wr3 = (radi*omega[i][2] + radj*omega[j][2]) * rinv;
+  
+  // normal force = Hertzian contact + normal velocity damping
+
+  double *rmass = atom->rmass;
+  double *mass = atom->mass;
+  int *mask = atom->mask;
+
+  if (rmass) {
+    meff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]);
+    if (mask[i] & freeze_group_bit) meff = rmass[j];
+    if (mask[j] & freeze_group_bit) meff = rmass[i];
+  } else {
+    meff = mass[itype]*mass[jtype] / (mass[itype]+mass[jtype]);
+    if (mask[i] & freeze_group_bit) meff = mass[jtype];
+    if (mask[j] & freeze_group_bit) meff = mass[itype];
+  }
+  
+  damp = meff*gamman*vnnr*rsqinv;
+  ccel = kn*(radsum-r)*rinv - damp;
+
+  // relative velocities
+
+  vtr1 = vt1 - (delz*wr2-dely*wr3);
+  vtr2 = vt2 - (delx*wr3-delz*wr1);
+  vtr3 = vt3 - (dely*wr1-delx*wr2);
+  vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3;
+  vrel = sqrt(vrel);
+  
+  // shear history effects
+  // neighprev = index of found neigh on previous call
+  // search entire jnum list of neighbors of I for neighbor J
+  // start from neighprev, since will typically be next neighbor
+  // reset neighprev to 0 as necessary
+
+  int *jlist = list->firstneigh[i];
+  int jnum = list->numneigh[i];
+  int *touch = list->listgranhistory->firstneigh[i];
+  double *allshear = list->listgranhistory->firstdouble[i];
+
+  for (int jj = 0; jj < jnum; jj++) {
+    neighprev++;
+    if (neighprev >= jnum) neighprev = 0;
+    if (touch[neighprev] == j) break;
+  }
+
+  double *shear = &allshear[3*neighprev];
+  shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + 
+		shear[2]*shear[2]);
+  
+  // rotate shear displacements
+  
+  rsht = shear[0]*delx + shear[1]*dely + shear[2]*delz;
+  rsht *= rsqinv;
+  
+  // tangential forces = shear + tangential velocity damping
+  
+  fs1 = - (kt*shear[0] + meff*gammat*vtr1);
+  fs2 = - (kt*shear[1] + meff*gammat*vtr2);
+  fs3 = - (kt*shear[2] + meff*gammat*vtr3);
+  
+  // rescale frictional displacements and forces if needed
+  
+  fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3);
+  fn = xmu * fabs(ccel*r);
+  
+  if (fs > fn) {
+    if (shrmag != 0.0) {
+      fs1 *= fn/fs;
+      fs2 *= fn/fs;
+      fs3 *= fn/fs;
+      fs *= fn/fs;
+    } else fs1 = fs2 = fs3 = fs = 0.0;
+  }
+  
+  // set all forces and return no energy
+  
+  fforce = ccel;
+  svector[0] = fs1;
+  svector[1] = fs2;
+  svector[2] = fs3;
+  svector[3] = fs;
+  return 0.0;
+}

src/GRANULAR/pair_gran_hooke_history.h

@@ -38,6 +38,7 @@ class PairGranHookeHistory : public Pair {
   void read_restart(FILE *);
   void write_restart_settings(FILE *);
   void read_restart_settings(FILE *);
+  double single(int, int, int, int, double, double, double, double &);
   void reset_dt();
 
  protected:
@@ -55,6 +56,8 @@ class PairGranHookeHistory : public Pair {
   double *onerad_dynamic,*onerad_frozen;
   double *maxrad_dynamic,*maxrad_frozen;
 
+  int neighprev;
+
   void allocate();
 };
 

src/compute_pair_local.cpp

@@ -13,6 +13,7 @@
 
 #include "math.h"
 #include "string.h"
+#include "stdlib.h"
 #include "compute_pair_local.h"
 #include "atom.h"
 #include "update.h"
@@ -29,6 +30,8 @@ using namespace LAMMPS_NS;
 
 #define DELTA 10000
 
+enum{DIST,ENG,FORCE,FX,FY,FZ,PN};
+
 /* ---------------------------------------------------------------------- */
 
 ComputePairLocal::ComputePairLocal(LAMMPS *lmp, int narg, char **arg) :
@@ -41,18 +44,32 @@ ComputePairLocal::ComputePairLocal(LAMMPS *lmp, int narg, char **arg) :
   if (nvalues == 1) size_local_cols = 0;
   else size_local_cols = nvalues;
 
-  dflag = eflag = fflag = -1;
+  pstyle = new int[nvalues];
-  nvalues = 0;
+  pindex = new int[nvalues];
 
-  int i;
+  nvalues = 0;
   for (int iarg = 3; iarg < narg; iarg++) {
-    i = iarg-3;
+    if (strcmp(arg[iarg],"dist") == 0) pstyle[nvalues++] = DIST;
-    if (strcmp(arg[iarg],"dist") == 0) dflag = nvalues++;
+    else if (strcmp(arg[iarg],"eng") == 0) pstyle[nvalues++] = ENG;
-    else if (strcmp(arg[iarg],"eng") == 0) eflag = nvalues++;
+    else if (strcmp(arg[iarg],"force") == 0) pstyle[nvalues++] = FORCE;
-    else if (strcmp(arg[iarg],"force") == 0) fflag = nvalues++;
+    else if (strcmp(arg[iarg],"fx") == 0) pstyle[nvalues++] = FX;
-    else error->all(FLERR,"Invalid keyword in compute pair/local command");
+    else if (strcmp(arg[iarg],"fy") == 0) pstyle[nvalues++] = FY;
+    else if (strcmp(arg[iarg],"fz") == 0) pstyle[nvalues++] = FZ;
+    else if (arg[iarg][0] == 'p') {
+      int n = atoi(&arg[iarg][1]);
+      if (n <= 0) error->all(FLERR,
+			     "Invalid keyword in compute pair/local command");
+      pstyle[nvalues] = PN;
+      pindex[nvalues++] = n-1;
+    } else error->all(FLERR,"Invalid keyword in compute pair/local command");
   }
 
+  // set singleflag if need to call pair->single()
+
+  singleflag = 0;
+  for (int i = 0; i < nvalues; i++)
+    if (pstyle[i] != DIST) singleflag = 1;
+
   nmax = 0;
   vector = NULL;
   array = NULL;
@@ -64,17 +81,24 @@ ComputePairLocal::~ComputePairLocal()
 {
   memory->destroy(vector);
   memory->destroy(array);
+  delete [] pstyle;
+  delete [] pindex;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputePairLocal::init()
 {
-  if (force->pair == NULL) 
+  if (singleflag && force->pair == NULL) 
     error->all(FLERR,"No pair style is defined for compute pair/local");
-  if (force->pair->single_enable == 0)
+  if (singleflag && force->pair->single_enable == 0)
     error->all(FLERR,"Pair style does not support compute pair/local");
 
+  for (int i = 0; i < nvalues; i++)
+    if (pstyle[i] == PN && pindex[i] >= force->pair->single_extra)
+      error->all(FLERR,"Pair style does not have single field"
+		 " requested by compute pair/local");
+
   // need an occasional half neighbor list
 
   int irequest = neighbor->request((void *) this);
@@ -101,7 +125,7 @@ void ComputePairLocal::compute_local()
   ncount = compute_pairs(0);
   if (ncount > nmax) reallocate(ncount);
   size_local_rows = ncount;
-  ncount = compute_pairs(1);
+  compute_pairs(1);
 }
 
 /* ----------------------------------------------------------------------
@@ -117,7 +141,7 @@ int ComputePairLocal::compute_pairs(int flag)
   double xtmp,ytmp,ztmp,delx,dely,delz;
   double rsq,eng,fpair,factor_coul,factor_lj;
   int *ilist,*jlist,*numneigh,**firstneigh;
-  double *dbuf,*ebuf,*fbuf;
+  double *ptr;
 
   double **x = atom->x;
   int *type = atom->type;
@@ -138,23 +162,13 @@ int ComputePairLocal::compute_pairs(int flag)
 
   // loop over neighbors of my atoms
   // skip if I or J are not in group
-
+  // for flag = 0, just count pair interactions within force cutoff
-  if (flag) {
+  // for flag = 1, calculate requested output fields
-    if (nvalues == 1) {
-      if (dflag >= 0) dbuf = vector;
-      if (eflag >= 0) ebuf = vector;
-      if (fflag >= 0) fbuf = vector;
-    } else {
-      if (dflag >= 0) dbuf = &array[0][dflag];
-      if (eflag >= 0) ebuf = &array[0][eflag];
-      if (fflag >= 0) fbuf = &array[0][fflag];
-    }
-  }
 
   Pair *pair = force->pair;
   double **cutsq = force->pair->cutsq;
 
-  m = n = 0;
+  m = 0;
   for (ii = 0; ii < inum; ii++) {
     i = ilist[ii];
     if (!(mask[i] & groupbit)) continue;
@@ -183,13 +197,37 @@ int ComputePairLocal::compute_pairs(int flag)
       if (rsq >= cutsq[itype][jtype]) continue;
 	
       if (flag) {
-	if (dflag >= 0) dbuf[n] = sqrt(rsq);
+	if (singleflag)
-	if (eflag >= 0 || fflag >= 0) {
 	  eng = pair->single(i,j,itype,jtype,rsq,factor_coul,factor_lj,fpair);
-	  if (eflag >= 0) ebuf[n] = eng;
+
-	  if (fflag >= 0) fbuf[n] = sqrt(rsq)*fpair;
+	if (nvalues == 1) ptr = &vector[m];
+	else ptr = array[m];
+
+	for (n = 0; n < nvalues; n++) {
+	  switch (pstyle[n]) {
+	  case DIST:
+	    ptr[n] = sqrt(rsq);
+	    break;
+	  case ENG:
+	    ptr[n] = eng;
+	    break;
+	  case FORCE:
+	    ptr[n] = sqrt(rsq)*fpair;
+	    break;
+	  case FX:
+	    ptr[n] = delx*fpair;
+	    break;
+	  case FY:
+	    ptr[n] = dely*fpair;
+	    break;
+	  case FZ:
+	    ptr[n] = delz*fpair;
+	    break;
+	  case PN:
+	    ptr[n] = pair->svector[pindex[n]];
+	    break;
+	  }
 	}
-	n += nvalues;
       }
 
       m++;

src/compute_pair_local.h

@@ -37,6 +37,10 @@ class ComputePairLocal : public Compute {
   int nvalues,dflag,eflag,fflag;
   int ncount;
 
+  int *pstyle;              // style of each requested output 
+  int *pindex;              // for pI, index of the output (0 to M-1)
+  int singleflag;
+
   int nmax;
   double *vector;
   double **array;

src/pair.cpp

@@ -58,6 +58,8 @@ Pair::Pair(LAMMPS *lmp) : Pointers(lmp)
 
   nextra = 0;
   pvector = NULL;
+  single_extra = 0;
+  svector = NULL;
 
   // pair_modify settings
 

src/pair.h

@@ -53,6 +53,9 @@ class Pair : protected Pointers {
   int nextra;                    // # of extra quantities pair style calculates
   double *pvector;               // vector of extra pair quantities
 
+  int single_extra;              // number of extra single values calculated
+  double *svector;               // vector of extra single quantities
+  
   class NeighList *list;         // standard neighbor list used by most pairs
   class NeighList *listhalf;     // half list used by some pairs
   class NeighList *listfull;     // full list used by some pairs

src/pair_hybrid.cpp

@@ -45,6 +45,8 @@ PairHybrid::~PairHybrid()
     delete [] keywords;
   }
 
+  delete [] svector;
+
   if (allocated) {
     memory->destroy(setflag);
     memory->destroy(cutsq);
@@ -265,19 +267,32 @@ void PairHybrid::settings(int narg, char **arg)
     if (styles[m]) comm_reverse = MAX(comm_reverse,styles[m]->comm_reverse);
   }
 
-  // single_enable = 0 if any sub-style = 0
+  // single_enable = 1 if any sub-style is set
   // respa_enable = 1 if any sub-style is set
   // no_virial_fdotr_compute = 1 if any sub-style is set
   // ghostneigh = 1 if any sub-style is set
 
+  single_enable = 0;
   for (m = 0; m < nstyles; m++)
-    if (styles[m]->single_enable == 0) single_enable = 0;
+    if (styles[m]->single_enable) single_enable = 1;
   for (m = 0; m < nstyles; m++)
     if (styles[m]->respa_enable) respa_enable = 1;
   for (m = 0; m < nstyles; m++)
     if (styles[m]->no_virial_fdotr_compute) no_virial_fdotr_compute = 1;
   for (m = 0; m < nstyles; m++)
     if (styles[m]->ghostneigh) ghostneigh = 1;
+
+  // single_extra = min of all sub-style single_extra
+  // allocate svector
+
+  single_extra = styles[0]->single_extra;
+  for (m = 1; m < nstyles; m++)
+    single_extra = MIN(single_extra,styles[m]->single_extra);
+
+  if (single_extra) {
+    delete [] svector;
+    svector = new double[single_extra];
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -613,10 +628,17 @@ double PairHybrid::single(int i, int j, int itype, int jtype,
   for (int m = 0; m < nmap[itype][jtype]; m++) {
     if (rsq < styles[map[itype][jtype][m]]->cutsq[itype][jtype]) {
       if (styles[map[itype][jtype][m]]->single_enable == 0)
-	error->all(FLERR,"Pair hybrid sub-style does not support single call");
+	error->one(FLERR,"Pair hybrid sub-style does not support single call");
+
       esum += styles[map[itype][jtype][m]]->
 	single(i,j,itype,jtype,rsq,factor_coul,factor_lj,fone);
       fforce += fone;
+
+      // copy substyle extra values into hybrid's svector
+
+      if (single_extra && styles[map[itype][jtype][m]]->single_extra)
+	for (m = 0; m < single_extra; m++)
+	  svector[m] = styles[map[itype][jtype][m]]->svector[m];
     }
   }
 
@@ -666,7 +688,8 @@ void *PairHybrid::extract(char *str, int &dim)
       double *p_newvalue = (double *) ptr;
       double newvalue = *p_newvalue;
       if (cutptr && newvalue != cutvalue)
-	error->all(FLERR,"Coulomb cutoffs of pair hybrid sub-styles do not match");
+	error->all(FLERR,
+		   "Coulomb cutoffs of pair hybrid sub-styles do not match");
       cutptr = ptr;
       cutvalue = newvalue;
     } else if (ptr) return ptr;

src/random_mars.cpp

@@ -27,7 +27,7 @@ RanMars::RanMars(LAMMPS *lmp, int seed) : Pointers(lmp)
   double s,t;
 
   if (seed <= 0 || seed > 900000000)
-    error->all(FLERR,"Invalid seed for Marsaglia random # generator");
+    error->one(FLERR,"Invalid seed for Marsaglia random # generator");
 
   save = 0;
   u = new double[97+1];

src/random_park.cpp

@@ -40,7 +40,7 @@ using namespace LAMMPS_NS;
 RanPark::RanPark(LAMMPS *lmp, int seed_init) : Pointers(lmp)
 {
   if (seed_init <= 0)
-    error->all(FLERR,"Invalid seed for Park random # generator");
+    error->one(FLERR,"Invalid seed for Park random # generator");
   seed = seed_init;
   save = 0;
 }