From cc44a57275762f812cfa3f61b8360a3f36f57cf4 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Fri, 25 Jul 2008 14:41:48 +0000
Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@1967
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/GRANULAR/atom_vec_granular.cpp |  2 +-
 src/PERI/atom_vec_peri.cpp         | 95 +++++++++++++++---------------
 src/PERI/atom_vec_peri.h           |  2 +-
 src/PERI/compute_damage_atom.cpp   |  2 +-
 src/PERI/fix_peri_neigh.cpp        | 29 ++++++---
 src/PERI/fix_peri_neigh.h          |  3 +-
 src/PERI/pair_peri_pmb.cpp         | 57 +++++++++---------
 src/atom.cpp                       |  8 +--
 src/atom.h                         |  5 +-
 src/set.cpp                        | 37 ++++++++----
 10 files changed, 133 insertions(+), 107 deletions(-)

diff --git a/src/GRANULAR/atom_vec_granular.cpp b/src/GRANULAR/atom_vec_granular.cpp
index b81541daf4..f7b8d97956 100644
--- a/src/GRANULAR/atom_vec_granular.cpp
+++ b/src/GRANULAR/atom_vec_granular.cpp
@@ -460,7 +460,7 @@ int AtomVecGranular::size_restart()
   int i;
 
   int nlocal = atom->nlocal;
-  int n = 19 * nlocal;
+  int n = 16 * nlocal;
 
   if (atom->nextra_restart)
     for (int iextra = 0; iextra < atom->nextra_restart; iextra++) 
diff --git a/src/PERI/atom_vec_peri.cpp b/src/PERI/atom_vec_peri.cpp
index dbc87887e1..1d06de95e3 100644
--- a/src/PERI/atom_vec_peri.cpp
+++ b/src/PERI/atom_vec_peri.cpp
@@ -37,13 +37,12 @@ AtomVec(lmp, narg, arg)
   comm_x_only = 0;
   size_comm = 4;
   size_reverse = 3;
-  size_border = 13;
+  size_border = 11;
   size_data_atom = 7;
   size_data_vel = 4;
-  xcol_data = 3;
+  xcol_data = 5;
 
-  atom->rmass_flag = atom->vfrac_flag = atom->s0_flag = 1;
-  atom->vinter_flag = 1;
+  atom->vfrac_flag = atom->density_flag = atom->rmass_flag = 1;
 }
 
 /* ----------------------------------------------------------------------
@@ -72,13 +71,13 @@ void AtomVecPeri::grow(int n)
 
   vfrac = atom->vfrac = (double *) 
     memory->srealloc(atom->vfrac,nmax*sizeof(double),"atom:vfrac");
+  density = atom->density = (double *)
+    memory->srealloc(atom->density,nmax*sizeof(double),"atom:density");
   rmass = atom->rmass = (double *) 
     memory->srealloc(atom->rmass,nmax*sizeof(double),"atom:rmass");
   s0 = atom->s0 = (double *) 
     memory->srealloc(atom->s0,nmax*sizeof(double),"atom:s0");
   x0 = atom->x0 = memory->grow_2d_double_array(atom->x0,nmax,3,"atom:x0");
-  vinter = atom->vinter = (double *) 
-    memory->srealloc(atom->vinter,nmax*sizeof(double),"atom:vinter");
  
   if (atom->nextra_grow)
     for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
@@ -101,12 +100,12 @@ void AtomVecPeri::copy(int i, int j)
   v[j][2] = v[i][2];
 
   vfrac[j] = vfrac[i];
+  density[j] = density[i];
   rmass[j] = rmass[i];
-  s0[j]    = s0[i];
+  s0[j] = s0[i];
   x0[j][0] = x0[i][0];
   x0[j][1] = x0[i][1];
   x0[j][2] = x0[i][2];
-  vinter[j] = vinter[i];
 
   if (atom->nextra_grow)
     for (int iextra = 0; iextra < atom->nextra_grow; iextra++) 
@@ -172,7 +171,7 @@ void AtomVecPeri::unpack_comm(int n, int first, double *buf)
     x[i][0] = buf[m++];
     x[i][1] = buf[m++];
     x[i][2] = buf[m++];
-    s0[i]   = buf[m++];
+    s0[i] = buf[m++];
   }
 }
 
@@ -234,12 +233,10 @@ int AtomVecPeri::pack_border(int n, int *list, double *buf,
       buf[m++] = type[j];
       buf[m++] = mask[j];
       buf[m++] = vfrac[j];
-      buf[m++] = rmass[j];
       buf[m++] = s0[j];
       buf[m++] = x0[j][0];
       buf[m++] = x0[j][1];
       buf[m++] = x0[j][2];
-      buf[m++] = vinter[j];
     }
   } else {
     if (domain->triclinic == 0) {
@@ -260,12 +257,10 @@ int AtomVecPeri::pack_border(int n, int *list, double *buf,
       buf[m++] = type[j];
       buf[m++] = mask[j];
       buf[m++] = vfrac[j];
-      buf[m++] = rmass[j];
       buf[m++] = s0[j];
       buf[m++] = x0[j][0];
       buf[m++] = x0[j][1];
       buf[m++] = x0[j][2];
-      buf[m++] = vinter[j];
     }
   }
   return m;
@@ -276,13 +271,11 @@ int AtomVecPeri::pack_border(int n, int *list, double *buf,
 int AtomVecPeri::pack_border_one(int i, double *buf)
 {
   buf[0] = vfrac[i];
-  buf[1] = rmass[i];
-  buf[2] = s0[i];
-  buf[3] = x0[i][0];
-  buf[4] = x0[i][1];
-  buf[5] = x0[i][2];
-  buf[6] = vinter[i];
-  return 7;
+  buf[1] = s0[i];
+  buf[2] = x0[i][0];
+  buf[3] = x0[i][1];
+  buf[4] = x0[i][2];
+  return 5;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -302,12 +295,10 @@ void AtomVecPeri::unpack_border(int n, int first, double *buf)
     type[i] = static_cast<int> (buf[m++]);
     mask[i] = static_cast<int> (buf[m++]);
     vfrac[i] = buf[m++];
-    rmass[i] = buf[m++];
     s0[i] = buf[m++];
     x0[i][0] = buf[m++];
     x0[i][1] = buf[m++];
     x0[i][2] = buf[m++];
-    vinter[i] = buf[m++];
   }
 }
 
@@ -316,13 +307,11 @@ void AtomVecPeri::unpack_border(int n, int first, double *buf)
 int AtomVecPeri::unpack_border_one(int i, double *buf)
 {
   vfrac[i] = buf[0];
-  rmass[i] = buf[1];
-  s0[i]    = buf[2];
-  x0[i][0] = buf[3];
-  x0[i][1] = buf[4];
-  x0[i][2] = buf[5];
-  vinter[i]= buf[6];
-  return 7;
+  s0[i] = buf[1];
+  x0[i][0] = buf[2];
+  x0[i][1] = buf[3];
+  x0[i][2] = buf[4];
+  return 5;
 }
 
 /* ----------------------------------------------------------------------
@@ -345,12 +334,12 @@ int AtomVecPeri::pack_exchange(int i, double *buf)
   buf[m++] = image[i];
 
   buf[m++] = vfrac[i];
+  buf[m++] = density[i];
   buf[m++] = rmass[i];
   buf[m++] = s0[i];
   buf[m++] = x0[i][0];
   buf[m++] = x0[i][1];
   buf[m++] = x0[i][2];
-  buf[m++] = vinter[i];
 
   if (atom->nextra_grow)
     for (int iextra = 0; iextra < atom->nextra_grow; iextra++) 
@@ -380,12 +369,12 @@ int AtomVecPeri::unpack_exchange(double *buf)
   image[nlocal] = static_cast<int> (buf[m++]);
 
   vfrac[nlocal] = buf[m++];
+  density[nlocal] = buf[m++];
   rmass[nlocal] = buf[m++];
   s0[nlocal] = buf[m++];
   x0[nlocal][0] = buf[m++];
   x0[nlocal][1] = buf[m++];
   x0[nlocal][2] = buf[m++];
-  vinter[nlocal] = buf[m++];
 
   if (atom->nextra_grow)
     for (int iextra = 0; iextra < atom->nextra_grow; iextra++) 
@@ -437,13 +426,13 @@ int AtomVecPeri::pack_restart(int i, double *buf)
   buf[m++] = v[i][1];
   buf[m++] = v[i][2];
  
-  buf[m++] = rmass[i];
   buf[m++] = vfrac[i];
+  buf[m++] = density[i];
+  buf[m++] = rmass[i];
   buf[m++] = s0[i];
   buf[m++] = x0[i][0];
   buf[m++] = x0[i][1];
   buf[m++] = x0[i][2];
-  buf[m++] = vinter[i];
  
   if (atom->nextra_restart)
     for (int iextra = 0; iextra < atom->nextra_restart; iextra++)
@@ -480,13 +469,13 @@ int AtomVecPeri::unpack_restart(double *buf)
   v[nlocal][1] = buf[m++];
   v[nlocal][2] = buf[m++];
  
-  rmass[nlocal] = buf[m++];
   vfrac[nlocal] = buf[m++];
+  density[nlocal] = buf[m++];
+  rmass[nlocal] = buf[m++];
   s0[nlocal] = buf[m++];
   x0[nlocal][0] = buf[m++];
   x0[nlocal][1] = buf[m++];
   x0[nlocal][2] = buf[m++];
-  vinter[nlocal] = buf[m++];
  
   double **extra = atom->extra;
   if (atom->nextra_store) {
@@ -520,12 +509,12 @@ void AtomVecPeri::create_atom(int itype, double *coord)
   v[nlocal][2] = 0.0;
  
   vfrac[nlocal] = 1.0;
-  rmass[nlocal] = 1.0;
+  density[nlocal] = 1.0;
+  rmass[nlocal] = density[nlocal];
   s0[nlocal] = MAXDOUBLE;
-  x0[nlocal][0] = x[nlocal][0];
-  x0[nlocal][1] = x[nlocal][1];
-  x0[nlocal][2] = x[nlocal][2];
-  vinter[nlocal] = 0.0;
+  x0[nlocal][0] = coord[0];
+  x0[nlocal][1] = coord[1];
+  x0[nlocal][2] = coord[2];
  
   atom->nlocal++;
 }
@@ -548,21 +537,25 @@ void AtomVecPeri::data_atom(double *coord, int imagetmp, char **values)
   if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes)
     error->one("Invalid atom type in Atoms section of data file");
 
-  vfrac[nlocal] = atof(values[5]);
-  rmass[nlocal] = atof(values[6]);
-  s0[nlocal] = MAXDOUBLE;
-  vinter[nlocal] = 0.0;
+  vfrac[nlocal] = atof(values[3]);
+  density[nlocal] = atof(values[4]);
+  rmass[nlocal] = density[nlocal];
 
   x[nlocal][0] = coord[0];
   x[nlocal][1] = coord[1];
   x[nlocal][2] = coord[2];
- 
+
   image[nlocal] = imagetmp;
  
   mask[nlocal] = 1;
   v[nlocal][0] = 0.0;
   v[nlocal][1] = 0.0;
   v[nlocal][2] = 0.0;
+
+  s0[nlocal] = MAXDOUBLE;
+  x0[nlocal][0] = coord[0];
+  x0[nlocal][1] = coord[1];
+  x0[nlocal][2] = coord[2];
  
   atom->nlocal++;
 }
@@ -576,13 +569,17 @@ void AtomVecPeri::data_atom(double *coord, int imagetmp, char **values)
 int AtomVecPeri::data_atom_hybrid(int nlocal, char **values)
 {
   vfrac[nlocal] = atof(values[0]);
-  rmass[nlocal] = atof(values[1]);
+  density[nlocal] = atof(values[1]);
+  rmass[nlocal] = density[nlocal];
 
-  s0[nlocal] = MAXDOUBLE;
-  vinter[nlocal] = 0.0;
   v[nlocal][0] = 0.0;
   v[nlocal][1] = 0.0;
   v[nlocal][2] = 0.0;
+
+  s0[nlocal] = MAXDOUBLE;
+  x0[nlocal][0] = x[nlocal][0];
+  x0[nlocal][1] = x[nlocal][1];
+  x0[nlocal][2] = x[nlocal][2];
  
   return 2;
 }
@@ -604,10 +601,10 @@ double AtomVecPeri::memory_usage()
   if (atom->memcheck("f")) bytes += nmax*3 * sizeof(double);
  
   if (atom->memcheck("vfrac")) bytes += nmax * sizeof(double);
+  if (atom->memcheck("density")) bytes += nmax * sizeof(double);
   if (atom->memcheck("rmass")) bytes += nmax * sizeof(double);
   if (atom->memcheck("s0")) bytes += nmax * sizeof(double);
   if (atom->memcheck("x0")) bytes += nmax*3 * sizeof(double);
-  if (atom->memcheck("vinter")) bytes += nmax * sizeof(double);
  
   return bytes;
 }
diff --git a/src/PERI/atom_vec_peri.h b/src/PERI/atom_vec_peri.h
index 0c785c2400..a115d90121 100755
--- a/src/PERI/atom_vec_peri.h
+++ b/src/PERI/atom_vec_peri.h
@@ -46,7 +46,7 @@ class AtomVecPeri : public AtomVec {
  private:
   int *tag,*type,*mask,*image;
   double **x,**v,**f;
-  double *rmass,*vfrac,*s0,**x0,*vinter;
+  double *vfrac,*density,*rmass,*s0,**x0;
 };
 
 }
diff --git a/src/PERI/compute_damage_atom.cpp b/src/PERI/compute_damage_atom.cpp
index b3549031dc..9a002439e9 100644
--- a/src/PERI/compute_damage_atom.cpp
+++ b/src/PERI/compute_damage_atom.cpp
@@ -87,8 +87,8 @@ void ComputeDamageAtom::compute_peratom()
   int nlocal = atom->nlocal;
   int *mask = atom->mask;
   double **x = atom->x;
-  double *vinter = atom->vinter;
   double *vfrac = atom->vfrac;
+  double *vinter = ((FixPeriNeigh *) modify->fix[ifix_peri])->vinter;
   int **partner = ((FixPeriNeigh *) modify->fix[ifix_peri])->partner;
   int *npartner = ((FixPeriNeigh *) modify->fix[ifix_peri])->npartner;
   int i,j,jj,jnum;
diff --git a/src/PERI/fix_peri_neigh.cpp b/src/PERI/fix_peri_neigh.cpp
index 509dfab3eb..624e346355 100644
--- a/src/PERI/fix_peri_neigh.cpp
+++ b/src/PERI/fix_peri_neigh.cpp
@@ -50,6 +50,8 @@ FixPeriNeigh::FixPeriNeigh(LAMMPS *lmp,int narg, char **arg) :
   npartner = NULL;
   partner = NULL;
   r0 = NULL;
+  vinter = NULL;
+
   grow_arrays(atom->nmax);
   atom->add_callback(0);
   atom->add_callback(1);
@@ -75,6 +77,7 @@ FixPeriNeigh::~FixPeriNeigh()
   memory->sfree(npartner);
   memory->destroy_2d_int_array(partner);
   memory->destroy_2d_double_array(r0);
+  memory->sfree(vinter);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -120,7 +123,6 @@ void FixPeriNeigh::setup(int vflag)
 
   double **x = atom->x;
   double *vfrac = atom->vfrac;
-  double *vinter = atom->vinter;
   int *type = atom->type;
   int *tag = atom->tag;
   int nlocal = atom->nlocal;
@@ -144,8 +146,6 @@ void FixPeriNeigh::setup(int vflag)
   Pair *anypair = force->pair_match("peri");
   double **cutsq = anypair->cutsq;
 
-  double *s0 = atom->s0;
-
   for (ii = 0; ii < inum; ii++) {
     i = ilist[ii];
     xtmp = x[i][0];
@@ -173,9 +173,8 @@ void FixPeriNeigh::setup(int vflag)
   int maxall;
   MPI_Allreduce(&maxpartner,&maxall,1,MPI_INT,MPI_MAX,world);
 
-  // realloc arrays stored by fix with correct value for maxpartner
+  // realloc arrays with correct value for maxpartner
 
-  memory->sfree(npartner);
   memory->destroy_2d_int_array(partner);
   memory->destroy_2d_double_array(r0);
   npartner = NULL;
@@ -184,8 +183,12 @@ void FixPeriNeigh::setup(int vflag)
   grow_arrays(atom->nmax);
 
   // create partner list and r0 values from neighbor list
+  // compute vinter for each atom
 
-  for (i = 0; i < nlocal; i++) npartner[i] = 0;
+  for (i = 0; i < nlocal; i++) {
+    npartner[i] = 0;
+    vinter[i] = 0.0;
+  }
 
   for (ii = 0; ii < inum; ii++) {
     i = ilist[ii];
@@ -195,7 +198,6 @@ void FixPeriNeigh::setup(int vflag)
     itype = type[i];
     jlist = firstneigh[i];
     jnum = numneigh[i];
-    vinter[i] = 0.0;
 
     for (jj = 0; jj < jnum; jj++) {
       j = jlist[jj];
@@ -246,6 +248,7 @@ double FixPeriNeigh::memory_usage()
   int bytes = nmax * sizeof(int);
   bytes += nmax*maxpartner * sizeof(int);
   bytes += nmax*maxpartner * sizeof(double);
+  bytes += nmax * sizeof(double);
   return bytes;
 }
 
@@ -260,6 +263,8 @@ void FixPeriNeigh::grow_arrays(int nmax)
   partner = memory->grow_2d_int_array(partner,nmax,maxpartner,
 				      "peri_neigh:partner");
   r0 = memory->grow_2d_double_array(r0,nmax,maxpartner,"peri_neigh:r0");
+  vinter = (double *) memory->srealloc(vinter,nmax*sizeof(double),
+				       "peri_neigh:vinter");
 }
 
 /* ----------------------------------------------------------------------
@@ -273,6 +278,7 @@ void FixPeriNeigh::copy_arrays(int i, int j)
     partner[j][m] = partner[i][m];
     r0[j][m] = r0[i][m];
   }
+  vinter[j] = vinter[i];
 }
 
 /* ----------------------------------------------------------------------
@@ -290,7 +296,9 @@ int FixPeriNeigh::pack_exchange(int i, double *buf)
     buf[m++] = partner[i][n];
     buf[m++] = r0[i][n];
   }
+
   buf[0] = m/2;
+  buf[m++] = vinter[i];
   return m;
 }
 
@@ -306,6 +314,7 @@ int FixPeriNeigh::unpack_exchange(int nlocal, double *buf)
     partner[nlocal][n] = static_cast<int> (buf[m++]);
     r0[nlocal][n] = buf[m++];
   }
+  vinter[nlocal] = buf[m++];
   return m;
 }
 
@@ -322,6 +331,7 @@ int FixPeriNeigh::pack_restart(int i, double *buf)
     buf[m++] = partner[i][n];
     buf[m++] = r0[i][n];
   }
+  buf[m++] = vinter[i];
   return m;
 }
 
@@ -344,6 +354,7 @@ void FixPeriNeigh::unpack_restart(int nlocal, int nth)
     partner[nlocal][n] = static_cast<int> (extra[nlocal][m++]);
     r0[nlocal][n] = extra[nlocal][m++];
   }
+  vinter[nlocal] = extra[nlocal][m++];
 }
 
 /* ----------------------------------------------------------------------
@@ -352,7 +363,7 @@ void FixPeriNeigh::unpack_restart(int nlocal, int nth)
 
 int FixPeriNeigh::maxsize_restart()
 {
-  return 2*maxpartner + 2;
+  return 2*maxpartner + 3;
 }
 
 /* ----------------------------------------------------------------------
@@ -361,5 +372,5 @@ int FixPeriNeigh::maxsize_restart()
 
 int FixPeriNeigh::size_restart(int nlocal)
 {
-  return 2*npartner[nlocal] + 2;
+  return 2*npartner[nlocal] + 3;
 }
diff --git a/src/PERI/fix_peri_neigh.h b/src/PERI/fix_peri_neigh.h
index 85f50e27b2..933153b8c4 100644
--- a/src/PERI/fix_peri_neigh.h
+++ b/src/PERI/fix_peri_neigh.h
@@ -41,11 +41,12 @@ class FixPeriNeigh : public Fix {
   int maxsize_restart();
 
  private:
-  int first;
+  int first;                 // flag for first time initialization
   int maxpartner;            // max # of peridynamic neighs for any atom
   int *npartner;             // # of neighbors for each atom
   int **partner;             // neighs for each atom, stored as global IDs
   double **r0;               // initial distance to partners
+  double *vinter;            // sum of vfrac for bonded neighbors
 
   class NeighList *list;
 };
diff --git a/src/PERI/pair_peri_pmb.cpp b/src/PERI/pair_peri_pmb.cpp
index 3d0b38f049..24dfb298f6 100644
--- a/src/PERI/pair_peri_pmb.cpp
+++ b/src/PERI/pair_peri_pmb.cpp
@@ -96,7 +96,6 @@ void PairPeriPMB::compute(int eflag, int vflag)
   double *vfrac = atom->vfrac;
   double *s0 = atom->s0;
   double **x0 = atom->x0;
-
   double **r0   = ((FixPeriNeigh *) modify->fix[ifix_peri])->r0;
   int **partner = ((FixPeriNeigh *) modify->fix[ifix_peri])->partner;
   int *npartner = ((FixPeriNeigh *) modify->fix[ifix_peri])->npartner;
@@ -112,13 +111,15 @@ void PairPeriPMB::compute(int eflag, int vflag)
 
   int nall = atom->nlocal + atom->nghost;
   int newton_pair = force->newton_pair;
+  int nonperiodic = domain->nonperiodic;
 
   inum = list->inum;
   ilist = list->ilist;
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
 
- // loop over neighbors of my atoms
+  // loop over neighbors of my atoms
+  // need minimg() for x0 difference since not ghosted
 
   for (ii = 0; ii < inum; ii++) {
     i = ilist[ii];
@@ -139,24 +140,23 @@ void PairPeriPMB::compute(int eflag, int vflag)
       delx = xtmp - x[j][0];
       dely = ytmp - x[j][1];
       delz = ztmp - x[j][2];
-      domain->minimum_image(delx,dely,delz);
       rsq = delx*delx + dely*dely + delz*delz;
       delx0 = xtmp0 - x0[j][0];
       dely0 = ytmp0 - x0[j][1];
       delz0 = ztmp0 - x0[j][2];
-      domain->minimum_image(delx0,dely0,delz0);
+      if (nonperiodic == 0) domain->minimum_image(delx0,dely0,delz0);
       rsq0 = delx0*delx0 + dely0*dely0 + delz0*delz0;
       jtype = type[j];
  
       r = sqrt(rsq);
 
-      // shortrange interaction distance based on initial particle position
+      // short-range interaction distance based on initial particle position
       // 0.9 and 1.35 are constants
 
       d_ij = MIN(0.9*sqrt(rsq0),1.35*lc);
 
-      // apply short-range contact forces
-      // 15 is a constant taken from the EMU Theory Manual
+      // short-range contact forces
+      // 15 is constant taken from the EMU Theory Manual
       // Silling, 12 May 2005, p 18
 
       if (r < d_ij) {
@@ -191,13 +191,13 @@ void PairPeriPMB::compute(int eflag, int vflag)
     s0_new = (double *) memory->smalloc(nmax*sizeof(double),"pair:s0_new");
   }
 
-  // first = flag indicating if this is first neighbor of particle i
+  // loop over my particles and their partners
+  // partner list contains all bond partners, so I-J appears twice
+  // if bond already broken, skip this partner
+  // first = true if this is first neighbor of particle i
 
   bool first;
 
-  // loop over my particles and their partners
-  // if bond already broken, skip this partner
-
   for (i = 0; i < nlocal; i++) {
     xtmp = x[i][0];
     ytmp = x[i][1];
@@ -211,22 +211,22 @@ void PairPeriPMB::compute(int eflag, int vflag)
       if (partner[i][jj] == 0) continue;
       j = atom->map(partner[i][jj]);
 
-      // check if we've lost a partner without first breaking bond
+      // check if lost a partner without first breaking bond
 
       if (j < 0) {
         partner[i][jj] = 0;
         continue;
       }
 
-      // compute force density and add to PD equation of motion
+      // compute force density, add to PD equation of motion
 
       delx = xtmp - x[j][0];
       dely = ytmp - x[j][1];
       delz = ztmp - x[j][2];
-      domain->minimum_image(delx,dely,delz);
+      if (nonperiodic == 0) domain->minimum_image(delx,dely,delz);
       rsq = delx*delx + dely*dely + delz*delz;
       jtype = type[j];
-      delta = sqrt(cutsq[itype][jtype]); // the horizon
+      delta = sqrt(cutsq[itype][jtype]);
       r = sqrt(rsq);
       dr = r - r0[i][jj];
 
@@ -234,7 +234,7 @@ void PairPeriPMB::compute(int eflag, int vflag)
 
       if (fabs(dr) < 2.2204e-016) dr = 0.0;
 
-      // apply scaling for vfrac[j] if particle j near the horizon
+      // scale vfrac[j] if particle j near the horizon
 
       if ((fabs(r0[i][jj] - delta)) <= half_lc)
         vfrac_scale = (-1.0/(2*half_lc))*(r0[i][jj]) + 
@@ -250,11 +250,13 @@ void PairPeriPMB::compute(int eflag, int vflag)
       f[i][1] += dely*fbond;
       f[i][2] += delz*fbond;
 
-      if (eflag) evdwl = rk*dr;
-      if (evflag) ev_tally(i,j,nlocal,newton_pair,
-			   evdwl,0.0,fbond,delx,dely,delz);
+      // hardwire newton flag off and use 1/2 since I-J is double counted
 
-      // find stretch in bond i-j and break if necessary
+      if (eflag) evdwl = rk*dr;
+      //if (evflag) ev_tally(i,j,nlocal,0,
+      //		   0.5*evdwl,0.0,0.5*fbond,delx,dely,delz);
+
+      // find stretch in bond I-J and break if necessary
       // use s0 from previous timestep
 
       if (stretch > MIN(s0[i],s0[j])) partner[i][jj] = 0;
@@ -266,14 +268,13 @@ void PairPeriPMB::compute(int eflag, int vflag)
       else
          s0_new[i] = MAX(s0_new[i],
 			 s00[itype][jtype] - (alpha[itype][jtype] * stretch));
-
-      // first neighbor of particle i has now been examined
-
       first = false;
     }
   }
 
-  // update with newly computed s0
+  if (vflag_fdotr) virial_compute();
+
+  // store new s0
 
   for (i = 0; i < nlocal; i++) s0[i] = s0_new[i]; 
 }
@@ -353,6 +354,7 @@ double PairPeriPMB::init_one(int i, int j)
   cutsq[j][i] = cutsq[i][j];
 
   // set other j,i parameters
+
   kspring[j][i] = kspring[i][j];
   alpha[j][i] = alpha[i][j];
   s00[j][i] = s00[i][j];
@@ -460,18 +462,17 @@ double PairPeriPMB::single(int i, int j, int itype, int jtype, double rsq,
 
   double *vfrac = atom->vfrac;
   double **x0 = atom->x0;
-  double lc = domain->lattice->xlattice;
-
-  double half_lc = 0.5*lc;
   double **r0   = ((FixPeriNeigh *) modify->fix[ifix_peri])->r0;
   int **partner = ((FixPeriNeigh *) modify->fix[ifix_peri])->partner;
   int *npartner = ((FixPeriNeigh *) modify->fix[ifix_peri])->npartner;
 
+  double lc = domain->lattice->xlattice;
+  double half_lc = 0.5*lc;
 
   delx0 = x0[i][0] - x0[j][0];
   dely0 = x0[i][1] - x0[j][1];
   delz0 = x0[i][2] - x0[j][2];
-  domain->minimum_image(delx0,dely0,delz0);
+  if (domain->nonperiodic == 0) domain->minimum_image(delx0,dely0,delz0);
   rsq0 = delx0*delx0 + dely0*dely0 + delz0*delz0;
 
   d_ij = MIN(0.9*sqrt(rsq0),1.35*lc);
diff --git a/src/atom.cpp b/src/atom.cpp
index f13b2cbc0e..8cfb2fd65b 100644
--- a/src/atom.cpp
+++ b/src/atom.cpp
@@ -64,7 +64,8 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
   q = NULL;
   mu = NULL;
   quat = omega = angmom = torque = NULL;
-  radius = density = rmass = vfrac = s0 = vinter = NULL;
+  radius = density = rmass = NULL;
+  vfrac = s0 = NULL;
   x0 = NULL;
 
   maxspecial = 1;
@@ -91,7 +92,7 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
   molecule_flag = 0;
   q_flag = mu_flag = 0;
   quat_flag = omega_flag = angmom_flag = torque_flag = 0;
-  radius_flag = density_flag = rmass_flag = vfrac_flag = s0_flag = vinter_flag = 0;
+  radius_flag = density_flag = rmass_flag = vfrac_flag = 0;
 
   // ntype-length arrays
 
@@ -164,7 +165,6 @@ Atom::~Atom()
   memory->sfree(rmass);
   memory->sfree(vfrac);
   memory->sfree(s0);
-  memory->sfree(vinter);
   memory->destroy_2d_double_array(x0);
 
   memory->sfree(molecule);
@@ -293,7 +293,7 @@ void Atom::init()
 }
 
 /* ----------------------------------------------------------------------
-   return 1 if style matches atom style hybrid sub-style
+   return 1 if style matches atom style or hybrid sub-style
    else return 0
 ------------------------------------------------------------------------- */
 
diff --git a/src/atom.h b/src/atom.h
index 4b9905b933..3690678841 100644
--- a/src/atom.h
+++ b/src/atom.h
@@ -48,7 +48,8 @@ class Atom : protected Pointers {
   int *molecule;
   double *q,**mu;
   double **quat,**omega,**angmom,**torque;
-  double *radius,*density,*rmass,*vfrac,*s0,**x0,*vinter;
+  double *radius,*density,*rmass,*vfrac,*s0;
+  double **x0;
 
   int maxspecial;
   int **nspecial,**special;
@@ -75,7 +76,7 @@ class Atom : protected Pointers {
   int molecule_flag;
   int q_flag,mu_flag;
   int quat_flag,omega_flag,angmom_flag,torque_flag;
-  int radius_flag,density_flag,rmass_flag,vfrac_flag,s0_flag,vinter_flag;
+  int radius_flag,density_flag,rmass_flag,vfrac_flag;
 
   // extra peratom info in restart file destined for fix & diag 
 
diff --git a/src/set.cpp b/src/set.cpp
index 12f9557c10..caf28a8e5f 100644
--- a/src/set.cpp
+++ b/src/set.cpp
@@ -185,10 +185,8 @@ void Set::command(int narg, char **arg)
     } else if (strcmp(arg[iarg],"density") == 0) {
       if (iarg+2 > narg) error->all("Illegal set command");
       dvalue = atof(arg[iarg+1]);
-      /*
       if (!atom->density_flag)
         error->all("Cannot set this attribute for this atom style");
-      */
       set(DENSITY);
       iarg += 2;
     } else if (strcmp(arg[iarg],"volume") == 0) {
@@ -302,6 +300,13 @@ void Set::selection(int n)
 
 void Set::set(int keyword)
 {
+  // set granular and peridynamic flags
+
+  int granflag = 0;
+  int periflag = 0;
+  if (atom->style_match("granular")) granflag = 1;
+  if (atom->style_match("peri")) periflag = 1;
+  
   if (keyword == DIPOLE) atom->check_dipole();
 
   selection(atom->nlocal);
@@ -318,6 +323,9 @@ void Set::set(int keyword)
       else if (keyword == VY) atom->v[i][1] = dvalue;
       else if (keyword == VZ) atom->v[i][2] = dvalue;
       else if (keyword == CHARGE) atom->q[i] = dvalue;
+
+      // diameter setting triggers setting of rmass
+
       else if (keyword == DIAMETER) {
 	atom->radius[i] = 0.5 * dvalue;
 	if (domain->dimension == 3)
@@ -326,17 +334,24 @@ void Set::set(int keyword)
 	else
 	  atom->rmass[i] = PI * 
 	    atom->radius[i]*atom->radius[i] * atom->density[i];
+
+      // density setting triggers setting of rmass
+      // for granular, rmass is function of diameter and density
+      // for peri, rmass stores density directly
+	
       } else if (keyword == DENSITY) {
-	atom->rmass[i] = dvalue;
-	/*
 	atom->density[i] = dvalue;
-	if (domain->dimension == 3)
-	  atom->rmass[i] = 4.0*PI/3.0 * 
-	    atom->radius[i]*atom->radius[i]*atom->radius[i] * atom->density[i];
-	else
-	  atom->rmass[i] = PI * 
-	    atom->radius[i]*atom->radius[i] * atom->density[i];
-	*/
+	if (granflag) {
+	  if (domain->dimension == 3)
+	    atom->rmass[i] = 4.0*PI/3.0 * 
+	      atom->radius[i]*atom->radius[i]*atom->radius[i] * 
+	      atom->density[i];
+	  else
+	    atom->rmass[i] = PI * 
+	      atom->radius[i]*atom->radius[i] * atom->density[i];
+	}
+	if (periflag) atom->rmass[i] = dvalue;
+
       } else if (keyword == VOLUME) atom->vfrac[i] = dvalue;
 
       else if (keyword == DIPOLE) {