diff --git a/src/GRANULAR/fix_pour.cpp b/src/GRANULAR/fix_pour.cpp
index 2255f64eb2..b54920f7b2 100644
--- a/src/GRANULAR/fix_pour.cpp
+++ b/src/GRANULAR/fix_pour.cpp
@@ -40,6 +40,7 @@ using namespace MathConst;
 
 enum{ATOM,MOLECULE};
 enum{ONE,RANGE,POLY};
+enum{CONSTANT,EQUAL};    // same as FixGravity
 
 #define EPSILON 0.001
 #define SMALL 1.0e-10
@@ -318,6 +319,10 @@ void FixPour::init()
   if (ifix == -1)
     error->all(FLERR,"No fix gravity defined for fix pour");
 
+  int varflag = ((FixGravity *) modify->fix[ifix])->varflag;
+  if (varflag != CONSTANT) 
+    error->all(FLERR,"Fix gravity for fix pour must be constant");
+
   double xgrav = ((FixGravity *) modify->fix[ifix])->xgrav;
   double ygrav = ((FixGravity *) modify->fix[ifix])->ygrav;
   double zgrav = ((FixGravity *) modify->fix[ifix])->zgrav;
diff --git a/src/MOLECULE/atom_vec_bond.cpp b/src/MOLECULE/atom_vec_bond.cpp
index 30ad7d83e1..00e8e3260d 100644
--- a/src/MOLECULE/atom_vec_bond.cpp
+++ b/src/MOLECULE/atom_vec_bond.cpp
@@ -60,6 +60,20 @@ AtomVecBond::~AtomVecBond()
   delete [] bond_negative;
 }
 
+/* ----------------------------------------------------------------------
+   grow atom arrays
+   must set local copy of body ptr
+   needed in replicate when 2 atom classes exist and pack_restart() is called
+------------------------------------------------------------------------- */
+
+void AtomVecBond::grow(int n)
+{
+  AtomVec::grow(n);
+  num_bond = atom->num_bond;
+  bond_type = atom->bond_type;
+}
+
+
 /* ----------------------------------------------------------------------
    modify values for AtomVec::pack_restart() to pack
 ------------------------------------------------------------------------- */
@@ -76,8 +90,8 @@ void AtomVecBond::pack_restart_pre(int ilocal)
 
   // flip any negative types to positive and flag which ones
 
-  int *num_bond = atom->num_bond;
-  int **bond_type = atom->bond_type;
+  //int *num_bond = atom->num_bond;
+  //int **bond_type = atom->bond_type;
 
   any_bond_negative = 0;
   for (int m = 0; m < num_bond[ilocal]; m++) {
@@ -98,8 +112,8 @@ void AtomVecBond::pack_restart_post(int ilocal)
   // restore the flagged types to their negative values
 
   if (any_bond_negative) {
-    int *num_bond = atom->num_bond;
-    int **bond_type = atom->bond_type;
+    //int *num_bond = atom->num_bond;
+    //int **bond_type = atom->bond_type;
     for (int m = 0; m < num_bond[ilocal]; m++)
       if (bond_negative[m]) bond_type[ilocal][m] = -bond_type[ilocal][m];
   }
diff --git a/src/MOLECULE/atom_vec_bond.h b/src/MOLECULE/atom_vec_bond.h
index c2bf7d6680..90c4b1f217 100644
--- a/src/MOLECULE/atom_vec_bond.h
+++ b/src/MOLECULE/atom_vec_bond.h
@@ -28,6 +28,7 @@ class AtomVecBond : public AtomVec {
  public:
   AtomVecBond(class LAMMPS *);
   ~AtomVecBond();
+  void grow(int);
   void pack_restart_pre(int);
   void pack_restart_post(int);
   void unpack_restart_init(int);
@@ -37,6 +38,9 @@ class AtomVecBond : public AtomVec {
   int any_bond_negative;
   int bond_per_atom;
   int *bond_negative;
+
+  int *num_bond;
+  int **bond_type;
 };
 
 }
diff --git a/src/MOLECULE/atom_vec_full.cpp b/src/MOLECULE/atom_vec_full.cpp
index 9ab0a296e0..be7c92ce93 100644
--- a/src/MOLECULE/atom_vec_full.cpp
+++ b/src/MOLECULE/atom_vec_full.cpp
@@ -128,7 +128,7 @@ void AtomVecFull::pack_restart_pre(int ilocal)
   int *num_improper = atom->num_improper;
   int **improper_type = atom->improper_type;
 
-  int any_bond_negative = 0;
+  any_bond_negative = 0;
   for (int m = 0; m < num_bond[ilocal]; m++) {
     if (bond_type[ilocal][m] < 0) {
       bond_negative[m] = 1;
@@ -137,7 +137,7 @@ void AtomVecFull::pack_restart_pre(int ilocal)
     } else bond_negative[m] = 0;
   }
 
-  int any_angle_negative = 0;
+  any_angle_negative = 0;
   for (int m = 0; m < num_angle[ilocal]; m++) {
     if (angle_type[ilocal][m] < 0) {
       angle_negative[m] = 1;
@@ -146,7 +146,7 @@ void AtomVecFull::pack_restart_pre(int ilocal)
     } else angle_negative[m] = 0;
   }
 
-  int any_dihedral_negative = 0;
+  any_dihedral_negative = 0;
   for (int m = 0; m < num_dihedral[ilocal]; m++) {
     if (dihedral_type[ilocal][m] < 0) {
       dihedral_negative[m] = 1;
@@ -155,7 +155,7 @@ void AtomVecFull::pack_restart_pre(int ilocal)
     } else dihedral_negative[m] = 0;
   }
 
-  int any_improper_negative = 0;
+  any_improper_negative = 0;
   for (int m = 0; m < num_improper[ilocal]; m++) {
     if (improper_type[ilocal][m] < 0) {
       improper_negative[m] = 1;
diff --git a/src/MOLECULE/atom_vec_molecular.cpp b/src/MOLECULE/atom_vec_molecular.cpp
index 52947ceb71..ea15216aee 100644
--- a/src/MOLECULE/atom_vec_molecular.cpp
+++ b/src/MOLECULE/atom_vec_molecular.cpp
@@ -128,7 +128,7 @@ void AtomVecMolecular::pack_restart_pre(int ilocal)
   int *num_improper = atom->num_improper;
   int **improper_type = atom->improper_type;
 
-  int any_bond_negative = 0;
+  any_bond_negative = 0;
   for (int m = 0; m < num_bond[ilocal]; m++) {
     if (bond_type[ilocal][m] < 0) {
       bond_negative[m] = 1;
@@ -137,7 +137,7 @@ void AtomVecMolecular::pack_restart_pre(int ilocal)
     } else bond_negative[m] = 0;
   }
 
-  int any_angle_negative = 0;
+  any_angle_negative = 0;
   for (int m = 0; m < num_angle[ilocal]; m++) {
     if (angle_type[ilocal][m] < 0) {
       angle_negative[m] = 1;
@@ -146,7 +146,7 @@ void AtomVecMolecular::pack_restart_pre(int ilocal)
     } else angle_negative[m] = 0;
   }
 
-  int any_dihedral_negative = 0;
+  any_dihedral_negative = 0;
   for (int m = 0; m < num_dihedral[ilocal]; m++) {
     if (dihedral_type[ilocal][m] < 0) {
       dihedral_negative[m] = 1;
@@ -155,7 +155,7 @@ void AtomVecMolecular::pack_restart_pre(int ilocal)
     } else dihedral_negative[m] = 0;
   }
 
-  int any_improper_negative = 0;
+  any_improper_negative = 0;
   for (int m = 0; m < num_improper[ilocal]; m++) {
     if (improper_type[ilocal][m] < 0) {
       improper_negative[m] = 1;
diff --git a/src/USER-AWPMD/atom_vec_wavepacket.cpp b/src/USER-AWPMD/atom_vec_wavepacket.cpp
index bce334a7b3..d643ae8e0a 100644
--- a/src/USER-AWPMD/atom_vec_wavepacket.cpp
+++ b/src/USER-AWPMD/atom_vec_wavepacket.cpp
@@ -18,13 +18,7 @@
 #include "atom_vec_wavepacket.h"
 #include <cstring>
 #include "atom.h"
-#include "comm.h"
-#include "domain.h"
-#include "modify.h"
-#include "fix.h"
-#include "memory.h"
 #include "error.h"
-#include "utils.h"
 
 using namespace LAMMPS_NS;
 
@@ -32,1079 +26,69 @@ using namespace LAMMPS_NS;
 
 AtomVecWavepacket::AtomVecWavepacket(LAMMPS *lmp) : AtomVec(lmp)
 {
-  comm_x_only = comm_f_only = 0;
-
   mass_type = 1;
   molecular = 0;
   forceclearflag = 1;
 
-  size_forward = 4; // coords[3]+radius[1]
-  size_reverse = 10; // force[3]+erforce[1]+ervelforce[1]+vforce[3]+csforce[2]
-  size_border = 10; // coords[3]+tag[1]+type[1]+mask[1]+q[1]+spin[1]+eradius[1]+etag[1]
-  size_velocity = 6; // +velocities[3]+ ervel[1]+cs[2]
-  size_data_atom = 11; // for input file: 1-tag 2-type 3-q 4-spin 5-eradius 6-etag 7-cs_re 8-cs_im 9-x 10-y 11-z
-  size_data_vel = 5; // for input file: vx vy vz ervel <??>
-  xcol_data = 9; // starting column for x data
-
   atom->wavepacket_flag = 1;
-  atom->electron_flag = 1; // compatible with eff
+
+  atom->electron_flag = 1;    // compatible with eff
   atom->q_flag = atom->spin_flag = atom->eradius_flag =
     atom->ervel_flag = atom->erforce_flag = 1;
+  atom->cs_flag = atom->csforce_flag = 
+    atom->vforce_flag = atom->ervelforce_flag = atom->etag_flag = 1;
 
-  atom->cs_flag = atom->csforce_flag = atom->vforce_flag = atom->ervelforce_flag = atom->etag_flag = 1;
+  // strings with peratom variables to include in each AtomVec method
+  // strings cannot contain fields in corresponding AtomVec default strings
+  // order of fields in a string does not matter
+  // except: fields_data_atom & fields_data_vel must match data file
+
+  fields_grow = (char *) 
+    "q spin eradius ervel erforce cs csforce "
+    "vforce ervelforce etag";
+  fields_copy = (char *) "q spin eradius ervel cs etag";
+  fields_comm = (char *) "eradius";
+  fields_comm_vel = (char *) "eradius ervel cs";
+  fields_reverse = (char *) "erforce ervelforce vforce csforce";
+  fields_border = (char *) "q spin eradius etag";
+  fields_border_vel = (char *) "q spin eradius etag ervel cs";
+  fields_exchange = (char *) "q spin eradius ervel etag cs";
+  fields_restart = (char *) "q spin eradius ervel etag cs";
+  fields_create = (char *) "q spin eradius ervel etag cs";
+  fields_data_atom = (char *) "id type q spin eradius etag cs x";
+  fields_data_vel = (char *) "id v ervel";
+
+  setup_fields();
 }
 
 /* ----------------------------------------------------------------------
-   grow atom-electron arrays
-   n = 0 grows arrays by a chunk
-   n > 0 allocates arrays to size n
+   clear extra forces starting at atom N
+   nbytes = # of bytes to clear for a per-atom vector
 ------------------------------------------------------------------------- */
 
-void AtomVecWavepacket::grow(int n)
-{
-  if (n == 0) grow_nmax();
-  else nmax = n;
-  atom->nmax = nmax;
-
-  tag = memory->grow(atom->tag,nmax,"atom:tag");
-  type = memory->grow(atom->type,nmax,"atom:type");
-  mask = memory->grow(atom->mask,nmax,"atom:mask");
-  image = memory->grow(atom->image,nmax,"atom:image");
-  x = memory->grow(atom->x,nmax,3,"atom:x");
-  v = memory->grow(atom->v,nmax,3,"atom:v");
-  f = memory->grow(atom->f,nmax*comm->nthreads,3,"atom:f");
-
-  q = memory->grow(atom->q,nmax,"atom:q");
-  spin = memory->grow(atom->spin,nmax,"atom:spin");
-  eradius = memory->grow(atom->eradius,nmax,"atom:eradius");
-  ervel = memory->grow(atom->ervel,nmax,"atom:ervel");
-  erforce = memory->grow(atom->erforce,nmax*comm->nthreads,"atom:erforce");
-
-  cs = memory->grow(atom->cs,2*nmax,"atom:cs");
-  csforce = memory->grow(atom->csforce,2*nmax,"atom:csforce");
-  vforce = memory->grow(atom->vforce,3*nmax,"atom:vforce");
-  ervelforce = memory->grow(atom->ervelforce,nmax,"atom:ervelforce");
-  etag = memory->grow(atom->etag,nmax,"atom:etag");
-
-  if (atom->nextra_grow)
-    for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
-      modify->fix[atom->extra_grow[iextra]]->grow_arrays(nmax);
-}
-
-/* ----------------------------------------------------------------------
-   reset local array ptrs
-------------------------------------------------------------------------- */
-
-void AtomVecWavepacket::grow_reset()
-{
-  tag = atom->tag; type = atom->type;
-  mask = atom->mask; image = atom->image;
-  x = atom->x; v = atom->v; f = atom->f;
-  q = atom->q;
-  eradius = atom->eradius; ervel = atom->ervel; erforce = atom->erforce;
-
-  cs = atom->cs;
-  csforce = atom->csforce;
-  vforce = atom->vforce;
-  ervelforce = atom->ervelforce;
-  etag = atom->etag;
-}
-
-/* ----------------------------------------------------------------------
-   copy atom I info to atom J
-------------------------------------------------------------------------- */
-
-void AtomVecWavepacket::copy(int i, int j, int delflag)
-{
-  tag[j] = tag[i];
-  type[j] = type[i];
-  mask[j] = mask[i];
-  image[j] = image[i];
-  x[j][0] = x[i][0];
-  x[j][1] = x[i][1];
-  x[j][2] = x[i][2];
-  v[j][0] = v[i][0];
-  v[j][1] = v[i][1];
-  v[j][2] = v[i][2];
-
-  q[j] = q[i];
-  spin[j] = spin[i];
-  eradius[j] = eradius[i];
-  ervel[j] = ervel[i];
-
-  cs[2*j] = cs[2*i];
-  cs[2*j+1] = cs[2*i+1];
-  etag[j] = etag[i];
-
-
-  if (atom->nextra_grow)
-    for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
-      modify->fix[atom->extra_grow[iextra]]->copy_arrays(i,j,delflag);
-}
-
-/* ---------------------------------------------------------------------- */
-
 void AtomVecWavepacket::force_clear(int n, size_t nbytes)
 {
-  memset(&erforce[n],0,nbytes);
-}
-
-/* ---------------------------------------------------------------------- */
-// this will be used as partial pack for unsplit Hartree packets (v, ervel not regarded as separate variables)
-
-int AtomVecWavepacket::pack_comm(int n, int *list, double *buf,
-                               int pbc_flag, int *pbc)
-{
-  int i,j,m;
-  double dx,dy,dz;
-
-  m = 0;
-  if (pbc_flag == 0) {
-    for (i = 0; i < n; i++) {
-      j = list[i];
-      buf[m++] = x[j][0];
-      buf[m++] = x[j][1];
-      buf[m++] = x[j][2];
-      buf[m++] = eradius[j];
-    }
-  } else {
-    if (domain->triclinic == 0) {
-      dx = pbc[0]*domain->xprd;
-      dy = pbc[1]*domain->yprd;
-      dz = pbc[2]*domain->zprd;
-    } else {
-      dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz;
-      dy = pbc[1]*domain->yprd + pbc[3]*domain->yz;
-      dz = pbc[2]*domain->zprd;
-    }
-    for (i = 0; i < n; i++) {
-      j = list[i];
-      buf[m++] = x[j][0] + dx;
-      buf[m++] = x[j][1] + dy;
-      buf[m++] = x[j][2] + dz;
-      buf[m++] = eradius[j];
-    }
-  }
-  return m;
-}
-
-/* ---------------------------------------------------------------------- */
-// this is a complete pack of all 'position' variables of AWPMD
-
-int AtomVecWavepacket::pack_comm_vel(int n, int *list, double *buf,
-                                   int pbc_flag, int *pbc)
-{
-  int i,j,m;
-  double dx,dy,dz,dvx,dvy,dvz;
-
-  m = 0;
-  if (pbc_flag == 0) {
-    for (i = 0; i < n; i++) {
-      j = list[i];
-      buf[m++] = x[j][0];
-      buf[m++] = x[j][1];
-      buf[m++] = x[j][2];
-      buf[m++] = eradius[j];
-      buf[m++] = v[j][0];
-      buf[m++] = v[j][1];
-      buf[m++] = v[j][2];
-
-      buf[m++] = ervel[j];
-      buf[m++] = cs[2*j];
-      buf[m++] = cs[2*j+1];
-    }
-  } else {
-    if (domain->triclinic == 0) {
-      dx = pbc[0]*domain->xprd;
-      dy = pbc[1]*domain->yprd;
-      dz = pbc[2]*domain->zprd;
-    } else {
-      dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz;
-      dy = pbc[1]*domain->yprd + pbc[3]*domain->yz;
-      dz = pbc[2]*domain->zprd;
-    }
-    if (!deform_vremap) {
-      for (i = 0; i < n; i++) {
-        j = list[i];
-        buf[m++] = x[j][0] + dx;
-        buf[m++] = x[j][1] + dy;
-        buf[m++] = x[j][2] + dz;
-        buf[m++] = eradius[j];
-        buf[m++] = v[j][0];
-        buf[m++] = v[j][1];
-        buf[m++] = v[j][2];
-
-        buf[m++] = ervel[j];
-        buf[m++] = cs[2*j];
-        buf[m++] = cs[2*j+1];
-      }
-    } else {
-      dvx = pbc[0]*h_rate[0] + pbc[5]*h_rate[5] + pbc[4]*h_rate[4];
-      dvy = pbc[1]*h_rate[1] + pbc[3]*h_rate[3];
-      dvz = pbc[2]*h_rate[2];
-      for (i = 0; i < n; i++) {
-        j = list[i];
-        buf[m++] = x[j][0] + dx;
-        buf[m++] = x[j][1] + dy;
-        buf[m++] = x[j][2] + dz;
-        buf[m++] = eradius[j];
-        if (mask[i] & deform_groupbit) {
-          buf[m++] = v[j][0] + dvx;
-          buf[m++] = v[j][1] + dvy;
-          buf[m++] = v[j][2] + dvz;
-        } else {
-          buf[m++] = v[j][0];
-          buf[m++] = v[j][1];
-          buf[m++] = v[j][2];
-        }
-        buf[m++] = ervel[j];
-        buf[m++] = cs[2*j];
-        buf[m++] = cs[2*j+1];
-      }
-    }
-  }
-  return m;
-}
-
-/* ---------------------------------------------------------------------- */
-
-int AtomVecWavepacket::pack_comm_hybrid(int n, int *list, double *buf)
-{
-  int i,j,m;
-
-  m = 0;
-  for (i = 0; i < n; i++) {
-    j = list[i];
-    buf[m++] = eradius[j];
-    buf[m++] = ervel[j];
-    buf[m++] = cs[2*j];
-    buf[m++] = cs[2*j+1];
-  }
-  return m;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void AtomVecWavepacket::unpack_comm(int n, int first, double *buf)
-{
-  int i,m,last;
-
-  m = 0;
-  last = first + n;
-  for (i = first; i < last; i++) {
-    x[i][0] = buf[m++];
-    x[i][1] = buf[m++];
-    x[i][2] = buf[m++];
-    eradius[i] = buf[m++];
-  }
-}
-
-/* ---------------------------------------------------------------------- */
-
-void AtomVecWavepacket::unpack_comm_vel(int n, int first, double *buf)
-{
-  int i,m,last;
-
-  m = 0;
-  last = first + n;
-  for (i = first; i < last; i++) {
-    x[i][0] = buf[m++];
-    x[i][1] = buf[m++];
-    x[i][2] = buf[m++];
-    eradius[i] = buf[m++];
-    v[i][0] = buf[m++];
-    v[i][1] = buf[m++];
-    v[i][2] = buf[m++];
-
-    ervel[i] = buf[m++];
-    cs[2*i] =  buf[m++];
-    cs[2*i+1] = buf[m++];
-  }
-}
-
-/* ---------------------------------------------------------------------- */
-
-int AtomVecWavepacket::unpack_comm_hybrid(int n, int first, double *buf)
-{
-  int i,m,last;
-
-  m = 0;
-  last = first + n;
-  for (i = first; i < last; i++){
-    eradius[i] = buf[m++];
-    ervel[i] = buf[m++];
-    cs[2*i] =  buf[m++];
-    cs[2*i+1] = buf[m++];
-  }
-  return m;
-}
-
-/* ---------------------------------------------------------------------- */
-
-int AtomVecWavepacket::pack_reverse(int n, int first, double *buf)
-{
-  int i,m,last;
-
-  m = 0;
-  last = first + n;
-  for (i = first; i < last; i++) { //10
-    buf[m++] = f[i][0];
-    buf[m++] = f[i][1];
-    buf[m++] = f[i][2];
-    buf[m++] = erforce[i];
-
-    buf[m++] = ervelforce[i];
-    buf[m++] = vforce[3*i];
-    buf[m++] = vforce[3*i+1];
-    buf[m++] = vforce[3*i+2];
-    buf[m++] = csforce[2*i];
-    buf[m++] = csforce[2*i+1];
-  }
-  return m;
-}
-
-/* ---------------------------------------------------------------------- */
-
-int AtomVecWavepacket::pack_reverse_hybrid(int n, int first, double *buf)
-{
-  int i,m,last;
-
-  m = 0;
-  last = first + n;
-  for (i = first; i < last; i++){
-    buf[m++] = erforce[i];
-
-    buf[m++] = ervelforce[i];
-    buf[m++] = vforce[3*i];
-    buf[m++] = vforce[3*i+1];
-    buf[m++] = vforce[3*i+2];
-    buf[m++] = csforce[2*i];
-    buf[m++] = csforce[2*i+1];
-  }
-  return m;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void AtomVecWavepacket::unpack_reverse(int n, int *list, double *buf)
-{
-  int i,j,m;
-
-  m = 0;
-  for (i = 0; i < n; i++) {
-    j = list[i];
-    f[j][0] += buf[m++];
-    f[j][1] += buf[m++];
-    f[j][2] += buf[m++];
-    erforce[j] += buf[m++];
-
-    ervelforce[j] += buf[m++];
-    vforce[3*j] += buf[m++];
-    vforce[3*j+1] += buf[m++];
-    vforce[3*j+2] += buf[m++];
-    csforce[2*j] += buf[m++];
-    csforce[2*j+1] += buf[m++];
-  }
-}
-
-/* ---------------------------------------------------------------------- */
-
-int AtomVecWavepacket::unpack_reverse_hybrid(int n, int *list, double *buf)
-{
-  int i,j,m;
-
-  m = 0;
-  for (i = 0; i < n; i++) {
-    j = list[i];
-    erforce[j] += buf[m++];
-
-    ervelforce[j] += buf[m++];
-    vforce[3*j] += buf[m++];
-    vforce[3*j+1] += buf[m++];
-    vforce[3*j+2] += buf[m++];
-    csforce[2*j] += buf[m++];
-    csforce[2*j+1] += buf[m++];
-  }
-  return m;
-}
-
-/* ---------------------------------------------------------------------- */
-// will be used for Hartree unsplit version (the etag is added however)
-int AtomVecWavepacket::pack_border(int n, int *list, double *buf,
-                                 int pbc_flag, int *pbc)
-{
-  int i,j,m;
-  double dx,dy,dz;
-
-  m = 0;
-  if (pbc_flag == 0) {
-    for (i = 0; i < n; i++) {
-      j = list[i];
-      buf[m++] = x[j][0];
-      buf[m++] = x[j][1];
-      buf[m++] = x[j][2];
-      buf[m++] = ubuf(tag[j]).d;
-      buf[m++] = ubuf(type[j]).d;
-      buf[m++] = ubuf(mask[j]).d;
-      buf[m++] = q[j];
-      buf[m++] = ubuf(spin[j]).d;
-      buf[m++] = eradius[j];
-      buf[m++] = ubuf(etag[j]).d;
-    }
-  } else {
-    if (domain->triclinic == 0) {
-      dx = pbc[0]*domain->xprd;
-      dy = pbc[1]*domain->yprd;
-      dz = pbc[2]*domain->zprd;
-    } else {
-      dx = pbc[0];
-      dy = pbc[1];
-      dz = pbc[2];
-    }
-    for (i = 0; i < n; i++) {
-      j = list[i];
-      buf[m++] = x[j][0] + dx;
-      buf[m++] = x[j][1] + dy;
-      buf[m++] = x[j][2] + dz;
-      buf[m++] = ubuf(tag[j]).d;
-      buf[m++] = ubuf(type[j]).d;
-      buf[m++] = ubuf(mask[j]).d;
-      buf[m++] = q[j];
-      buf[m++] = ubuf(spin[j]).d;
-      buf[m++] = eradius[j];
-      buf[m++] = ubuf(etag[j]).d;
-    }
-  }
-
-  if (atom->nextra_border)
-    for (int iextra = 0; iextra < atom->nextra_border; iextra++)
-      m += modify->fix[atom->extra_border[iextra]]->pack_border(n,list,&buf[m]);
-
-  return m;
-}
-
-/* ---------------------------------------------------------------------- */
-
-int AtomVecWavepacket::pack_border_vel(int n, int *list, double *buf,
-                                     int pbc_flag, int *pbc)
-{
-  int i,j,m;
-  double dx,dy,dz,dvx,dvy,dvz;
-
-  m = 0;
-  if (pbc_flag == 0) {
-    for (i = 0; i < n; i++) {
-      j = list[i];
-      buf[m++] = x[j][0];
-      buf[m++] = x[j][1];
-      buf[m++] = x[j][2];
-      buf[m++] = ubuf(tag[j]).d;
-      buf[m++] = ubuf(type[j]).d;
-      buf[m++] = ubuf(mask[j]).d;
-      buf[m++] = q[j];
-      buf[m++] = ubuf(spin[j]).d;
-      buf[m++] = eradius[j];
-      buf[m++] = ubuf(etag[j]).d;
-
-      buf[m++] = v[j][0];
-      buf[m++] = v[j][1];
-      buf[m++] = v[j][2];
-
-      buf[m++] = ervel[j];
-      buf[m++] = cs[2*j];
-      buf[m++] = cs[2*j+1];
-    }
-  } else {
-    if (domain->triclinic == 0) {
-      dx = pbc[0]*domain->xprd;
-      dy = pbc[1]*domain->yprd;
-      dz = pbc[2]*domain->zprd;
-    } else {
-      dx = pbc[0];
-      dy = pbc[1];
-      dz = pbc[2];
-    }
-    if (domain->triclinic == 0) {
-      for (i = 0; i < n; i++) {
-        j = list[i];
-        buf[m++] = x[j][0] + dx;
-        buf[m++] = x[j][1] + dy;
-        buf[m++] = x[j][2] + dz;
-        buf[m++] = ubuf(tag[j]).d;
-        buf[m++] = ubuf(type[j]).d;
-        buf[m++] = ubuf(mask[j]).d;
-        buf[m++] = q[j];
-        buf[m++] = ubuf(spin[j]).d;
-        buf[m++] = eradius[j];
-        buf[m++] = ubuf(etag[j]).d;
-
-        buf[m++] = v[j][0];
-        buf[m++] = v[j][1];
-        buf[m++] = v[j][2];
-
-
-        buf[m++] = ervel[j];
-        buf[m++] = cs[2*j];
-        buf[m++] = cs[2*j+1];
-      }
-    } else {
-      dvx = pbc[0]*h_rate[0] + pbc[5]*h_rate[5] + pbc[4]*h_rate[4];
-      dvy = pbc[1]*h_rate[1] + pbc[3]*h_rate[3];
-      dvz = pbc[2]*h_rate[2];
-      for (i = 0; i < n; i++) {
-        j = list[i];
-        buf[m++] = x[j][0] + dx;
-        buf[m++] = x[j][1] + dy;
-        buf[m++] = x[j][2] + dz;
-        buf[m++] = ubuf(tag[j]).d;
-        buf[m++] = ubuf(type[j]).d;
-        buf[m++] = ubuf(mask[j]).d;
-        buf[m++] = q[j];
-        buf[m++] = ubuf(spin[j]).d;
-        buf[m++] = eradius[j];
-        buf[m++] = ubuf(etag[j]).d;
-
-        if (mask[i] & deform_groupbit) {
-          buf[m++] = v[j][0] + dvx;
-          buf[m++] = v[j][1] + dvy;
-          buf[m++] = v[j][2] + dvz;
-        } else {
-          buf[m++] = v[j][0];
-          buf[m++] = v[j][1];
-          buf[m++] = v[j][2];
-        }
-
-        buf[m++] = ervel[j];
-        buf[m++] = cs[2*j];
-        buf[m++] = cs[2*j+1];
-      }
-    }
-  }
-
-  if (atom->nextra_border)
-    for (int iextra = 0; iextra < atom->nextra_border; iextra++)
-      m += modify->fix[atom->extra_border[iextra]]->pack_border(n,list,&buf[m]);
-
-  return m;
-}
-
-/* ---------------------------------------------------------------------- */
-
-int AtomVecWavepacket::pack_border_hybrid(int n, int *list, double *buf)
-{
-  int i,j,m;
-
-  m = 0;
-  for (i = 0; i < n; i++) {
-    j = list[i];
-    buf[m++] = q[j];
-    buf[m++] = ubuf(spin[j]).d;
-    buf[m++] = eradius[j];
-
-    buf[m++] = ubuf(etag[j]).d;
-    buf[m++] = ervel[j];
-    buf[m++] = cs[2*j];
-    buf[m++] = cs[2*j+1];
-  }
-  return m;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void AtomVecWavepacket::unpack_border(int n, int first, double *buf)
-{
-  int i,m,last;
-
-  m = 0;
-  last = first + n;
-  for (i = first; i < last; i++) {
-    if (i == nmax) grow(0);
-    x[i][0] = buf[m++];
-    x[i][1] = buf[m++];
-    x[i][2] = buf[m++];
-    tag[i] = (tagint) ubuf(buf[m++]).i;
-    type[i] = (int) ubuf(buf[m++]).i;
-    mask[i] = (int) ubuf(buf[m++]).i;
-    q[i] = buf[m++];
-    spin[i] = (int) ubuf(buf[m++]).i;
-    eradius[i] = buf[m++];
-    etag[i] = (int) ubuf(buf[m++]).i;
-  }
-
-  if (atom->nextra_border)
-    for (int iextra = 0; iextra < atom->nextra_border; iextra++)
-      m += modify->fix[atom->extra_border[iextra]]->
-        unpack_border(n,first,&buf[m]);
-}
-
-/* ---------------------------------------------------------------------- */
-
-void AtomVecWavepacket::unpack_border_vel(int n, int first, double *buf)
-{
-  int i,m,last;
-
-  m = 0;
-  last = first + n;
-  for (i = first; i < last; i++) {
-    if (i == nmax) grow(0);
-    x[i][0] = buf[m++];
-    x[i][1] = buf[m++];
-    x[i][2] = buf[m++];
-    tag[i] = (tagint) ubuf(buf[m++]).i;
-    type[i] = (int) ubuf(buf[m++]).i;
-    mask[i] = (int) ubuf(buf[m++]).i;
-    q[i] = buf[m++];
-    spin[i] = (int) ubuf(buf[m++]).i;
-    eradius[i] = buf[m++];
-    etag[i] = (int) ubuf(buf[m++]).i;
-    v[i][0] = buf[m++];
-    v[i][1] = buf[m++];
-    v[i][2] = buf[m++];
-    ervel[i] = buf[m++];
-    cs[2*i] = buf[m++];
-    cs[2*i+1] = buf[m++];
-  }
-
-  if (atom->nextra_border)
-    for (int iextra = 0; iextra < atom->nextra_border; iextra++)
-      m += modify->fix[atom->extra_border[iextra]]->
-        unpack_border(n,first,&buf[m]);
-}
-
-/* ---------------------------------------------------------------------- */
-
-int AtomVecWavepacket::unpack_border_hybrid(int n, int first, double *buf)
-{
-  int i,m,last;
-
-  m = 0;
-  last = first + n;
-  for (i = first; i < last; i++) {
-    q[i] = buf[m++];
-    spin[i] = (int) ubuf(buf[m++]).i;
-    eradius[i] = buf[m++];
-    etag[i] = (int) ubuf(buf[m++]).i;
-    ervel[i] = buf[m++];
-    cs[2*i] = buf[m++];
-    cs[2*i+1] = buf[m++];
-  }
-  return m;
+  memset(&atom->erforce[n],0,nbytes);
 }
 
 /* ----------------------------------------------------------------------
-   pack data for atom I for sending to another proc
-   xyz must be 1st 3 values, so comm::exchange() can test on them
+   initialize non-zero atom quantities
+   make each atom a proton
 ------------------------------------------------------------------------- */
 
-int AtomVecWavepacket::pack_exchange(int i, double *buf)
+void AtomVecWavepacket::create_atom_post(int ilocal)
 {
-  int m = 1;
-  buf[m++] = x[i][0];
-  buf[m++] = x[i][1];
-  buf[m++] = x[i][2];
-  buf[m++] = v[i][0];
-  buf[m++] = v[i][1];
-  buf[m++] = v[i][2];
-  buf[m++] = ubuf(tag[i]).d;
-  buf[m++] = ubuf(type[i]).d;
-  buf[m++] = ubuf(mask[i]).d;
-  buf[m++] = ubuf(image[i]).d;
-
-  buf[m++] = q[i];
-  buf[m++] = ubuf(spin[i]).d;
-  buf[m++] = eradius[i];
-  buf[m++] = ervel[i];
-
-  buf[m++] = ubuf(etag[i]).d;
-  buf[m++] = cs[2*i];
-  buf[m++] = cs[2*i+1];
-
-  if (atom->nextra_grow)
-    for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
-      m += modify->fix[atom->extra_grow[iextra]]->pack_exchange(i,&buf[m]);
-
-  buf[0] = m;
-  return m;
-}
-
-/* ---------------------------------------------------------------------- */
-
-int AtomVecWavepacket::unpack_exchange(double *buf)
-{
-  int nlocal = atom->nlocal;
-  if (nlocal == nmax) grow(0);
-
-  int m = 1;
-  x[nlocal][0] = buf[m++];
-  x[nlocal][1] = buf[m++];
-  x[nlocal][2] = buf[m++];
-  v[nlocal][0] = buf[m++];
-  v[nlocal][1] = buf[m++];
-  v[nlocal][2] = buf[m++];
-  tag[nlocal] = (tagint) ubuf(buf[m++]).i;
-  type[nlocal] = (int) ubuf(buf[m++]).i;
-  mask[nlocal] = (int) ubuf(buf[m++]).i;
-  image[nlocal] = (imageint) ubuf(buf[m++]).i;
-
-  q[nlocal] = buf[m++];
-  spin[nlocal] = (int) ubuf(buf[m++]).i;
-  eradius[nlocal] = buf[m++];
-  ervel[nlocal] = buf[m++];
-
-  etag[nlocal] = (int) ubuf(buf[m++]).i;
-  cs[2*nlocal] = buf[m++];
-  cs[2*nlocal+1] = buf[m++];
-
-  if (atom->nextra_grow)
-    for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
-      m += modify->fix[atom->extra_grow[iextra]]->
-        unpack_exchange(nlocal,&buf[m]);
-
-  atom->nlocal++;
-  return m;
+  atom->q[ilocal] = 1.0;
 }
 
 /* ----------------------------------------------------------------------
-   size of restart data for all atoms owned by this proc
-   include extra data stored by fixes
+   modify what AtomVec::data_atom() just unpacked
+   or initialize other atom quantities
 ------------------------------------------------------------------------- */
 
-int AtomVecWavepacket::size_restart()
+void AtomVecWavepacket::data_atom_post(int ilocal)
 {
-  int i;
-
-  int nlocal = atom->nlocal;
-  int n = 18 * nlocal;        // Associated with pack_restart
-
-  if (atom->nextra_restart)
-    for (int iextra = 0; iextra < atom->nextra_restart; iextra++)
-      for (i = 0; i < nlocal; i++)
-        n += modify->fix[atom->extra_restart[iextra]]->size_restart(i);
-
-  return n;
-}
-
-/* ----------------------------------------------------------------------
-   pack atom I's data for restart file including extra quantities
-   xyz must be 1st 3 values, so that read_restart can test on them
-   molecular types may be negative, but write as positive
-------------------------------------------------------------------------- */
-
-int AtomVecWavepacket::pack_restart(int i, double *buf)
-{
-  int m = 1;
-  buf[m++] = x[i][0];
-  buf[m++] = x[i][1];
-  buf[m++] = x[i][2];
-  buf[m++] = ubuf(tag[i]).d;
-  buf[m++] = ubuf(type[i]).d;
-  buf[m++] = ubuf(mask[i]).d;
-  buf[m++] = ubuf(image[i]).d;
-  buf[m++] = v[i][0];
-  buf[m++] = v[i][1];
-  buf[m++] = v[i][2];
-
-  buf[m++] = q[i];
-  buf[m++] = ubuf(spin[i]).d;
-  buf[m++] = eradius[i];
-  buf[m++] = ervel[i];
-
-  buf[m++] = ubuf(etag[i]).d;
-  buf[m++] = cs[2*i];
-  buf[m++] = cs[2*i+1];
-
-  if (atom->nextra_restart)
-    for (int iextra = 0; iextra < atom->nextra_restart; iextra++)
-      m += modify->fix[atom->extra_restart[iextra]]->pack_restart(i,&buf[m]);
-
-  buf[0] = m;
-  return m;
-}
-
-/* ----------------------------------------------------------------------
-   unpack data for one atom from restart file including extra quantities
-------------------------------------------------------------------------- */
-
-int AtomVecWavepacket::unpack_restart(double *buf)
-{
-  int nlocal = atom->nlocal;
-  if (nlocal == nmax) {
-    grow(0);
-    if (atom->nextra_store)
-      memory->grow(atom->extra,nmax,atom->nextra_store,"atom:extra");
-  }
-
-  int m = 1;
-  x[nlocal][0] = buf[m++];
-  x[nlocal][1] = buf[m++];
-  x[nlocal][2] = buf[m++];
-  tag[nlocal] = (tagint) ubuf(buf[m++]).i;
-  type[nlocal] = (int) ubuf(buf[m++]).i;
-  mask[nlocal] = (int) ubuf(buf[m++]).i;
-  image[nlocal] = (imageint) ubuf(buf[m++]).i;
-  v[nlocal][0] = buf[m++];
-  v[nlocal][1] = buf[m++];
-  v[nlocal][2] = buf[m++];
-
-  q[nlocal] = buf[m++];
-  spin[nlocal] = (int) ubuf(buf[m++]).i;
-  eradius[nlocal] = buf[m++];
-  ervel[nlocal] = buf[m++];
-
-  etag[nlocal] = (int) ubuf(buf[m++]).i;
-  cs[2*nlocal] = buf[m++];
-  cs[2*nlocal+1] = buf[m++];
-
-  double **extra = atom->extra;
-  if (atom->nextra_store) {
-    int size = static_cast<int> (buf[0]) - m;
-    for (int i = 0; i < size; i++) extra[nlocal][i] = buf[m++];
-  }
-
-  atom->nlocal++;
-  return m;
-}
-
-/* ----------------------------------------------------------------------
-   create one atom of itype at coord
-   set other values to defaults
-   AWPMD: creates a proton
-------------------------------------------------------------------------- */
-
-void AtomVecWavepacket::create_atom(int itype, double *coord)
-{
-  int nlocal = atom->nlocal;
-  if (nlocal == nmax) grow(0);
-
-  tag[nlocal] = 0;
-  type[nlocal] = itype;
-  x[nlocal][0] = coord[0];
-  x[nlocal][1] = coord[1];
-  x[nlocal][2] = coord[2];
-  mask[nlocal] = 1;
-  image[nlocal] = ((imageint) IMGMAX << IMG2BITS) |
-    ((imageint) IMGMAX << IMGBITS) | IMGMAX;
-  v[nlocal][0] = 0.0;
-  v[nlocal][1] = 0.0;
-  v[nlocal][2] = 0.0;
-
-  q[nlocal] = 1.;
-  spin[nlocal] = 0;
-  eradius[nlocal] = 0.0;
-  ervel[nlocal] = 0.0;
-
-  etag[nlocal] = 0;
-  cs[2*nlocal] = 0.;
-  cs[2*nlocal+1] = 0.;
-
-  atom->nlocal++;
-}
-
-/* ----------------------------------------------------------------------
-   unpack one line from Atoms section of data file
-   initialize other atom quantities
-   AWPMD: 0-tag 1-type 2-q 3-spin 4-eradius 5-etag 6-cs_re 7-cs_im
-------------------------------------------------------------------------- */
-
-void AtomVecWavepacket::data_atom(double *coord, imageint imagetmp,
-                                  char **values)
-{
-  int nlocal = atom->nlocal;
-
-  if (nlocal == nmax) grow(0);
-
-  tag[nlocal] = utils::tnumeric(FLERR,values[0],true,lmp);
-  type[nlocal] = utils::inumeric(FLERR,values[1],true,lmp);
-  if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes)
-    error->one(FLERR,"Invalid atom type in Atoms section of data file");
-
-  q[nlocal] = utils::numeric(FLERR,values[2],true,lmp);
-  spin[nlocal] = utils::inumeric(FLERR,values[3],true,lmp);
-  eradius[nlocal] = utils::numeric(FLERR,values[4],true,lmp);
-  if (eradius[nlocal] < 0.0)
-    error->one(FLERR,"Invalid eradius in Atoms section of data file");
-
-  etag[nlocal] = utils::inumeric(FLERR,values[5],true,lmp);
-  cs[2*nlocal] = utils::numeric(FLERR,values[6],true,lmp);
-  cs[2*nlocal+1] = utils::numeric(FLERR,values[7],true,lmp);
-
-  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;
-  ervel[nlocal] = 0.0;
-
-  atom->nlocal++;
-}
-
-/* ----------------------------------------------------------------------
-   unpack hybrid quantities from one line in Atoms section of data file
-   initialize other atom quantities for this sub-style
-------------------------------------------------------------------------- */
-
-int AtomVecWavepacket::data_atom_hybrid(int nlocal, char **values)
-{
-  q[nlocal] = utils::numeric(FLERR,values[0],true,lmp);
-  spin[nlocal] = utils::inumeric(FLERR,values[1],true,lmp);
-  eradius[nlocal] = utils::numeric(FLERR,values[2],true,lmp);
-  if (eradius[nlocal] < 0.0)
-    error->one(FLERR,"Invalid eradius in Atoms section of data file");
-
-  etag[nlocal] = utils::inumeric(FLERR,values[3],true,lmp);
-  cs[2*nlocal] = utils::inumeric(FLERR,values[4],true,lmp);
-  cs[2*nlocal+1] = utils::numeric(FLERR,values[5],true,lmp);
-
-  v[nlocal][0] = 0.0;
-  v[nlocal][1] = 0.0;
-  v[nlocal][2] = 0.0;
-  ervel[nlocal] = 0.0;
-
-  return 3;
-}
-
-/* ----------------------------------------------------------------------
-   unpack one line from Velocities section of data file
-------------------------------------------------------------------------- */
-
-void AtomVecWavepacket::data_vel(int m, char **values)
-{
-  v[m][0] = utils::numeric(FLERR,values[0],true,lmp);
-  v[m][1] = utils::numeric(FLERR,values[1],true,lmp);
-  v[m][2] = utils::numeric(FLERR,values[2],true,lmp);
-  ervel[m] = utils::numeric(FLERR,values[3],true,lmp);
-}
-
-/* ----------------------------------------------------------------------
-   unpack hybrid quantities from one line in Velocities section of data file
-------------------------------------------------------------------------- */
-
-int AtomVecWavepacket::data_vel_hybrid(int m, char **values)
-{
-  ervel[m] = utils::numeric(FLERR,values[0],true,lmp);
-  return 1;
-}
-
-/* ----------------------------------------------------------------------
-   pack atom info for data file including 3 image flags
-------------------------------------------------------------------------- */
-
-void AtomVecWavepacket::pack_data(double **buf)
-{
-  int nlocal = atom->nlocal;
-  for (int i = 0; i < nlocal; i++) {
-    buf[i][0] = ubuf(tag[i]).d;
-    buf[i][1] = ubuf(type[i]).d;
-    buf[i][2] = q[i];
-    buf[i][3] = ubuf(spin[i]).d;
-    buf[i][4] = eradius[i];
-    buf[i][5] = ubuf(etag[i]).d;
-    buf[i][6] = cs[2*i];
-    buf[i][7] = cs[2*i+1];
-    buf[i][8] = x[i][0];
-    buf[i][9] = x[i][1];
-    buf[i][10] = x[i][2];
-    buf[i][11] = ubuf((image[i] & IMGMASK) - IMGMAX).d;
-    buf[i][12] = ubuf((image[i] >> IMGBITS & IMGMASK) - IMGMAX).d;
-    buf[i][13] = ubuf((image[i] >> IMG2BITS) - IMGMAX).d;
-  }
-}
-
-/* ----------------------------------------------------------------------
-   pack hybrid atom info for data file
-------------------------------------------------------------------------- */
-
-int AtomVecWavepacket::pack_data_hybrid(int i, double *buf)
-{
-  buf[0] = q[i];
-  buf[1] = ubuf(spin[i]).d;
-  buf[2] = eradius[i];
-  buf[3] = ubuf(etag[i]).d;
-  buf[4] = cs[2*i];
-  buf[5] = cs[2*i+1];
-  return 6;
-}
-
-/* ----------------------------------------------------------------------
-   write atom info to data file including 3 image flags
-------------------------------------------------------------------------- */
-
-void AtomVecWavepacket::write_data(FILE *fp, int n, double **buf)
-{
-  for (int i = 0; i < n; i++)
-    fprintf(fp,TAGINT_FORMAT
-            " %d %-1.16e %d %-1.16e %d %-1.16e %-1.16e %-1.16e "
-            "%-1.16e %-1.16e %d %d %d\n",
-            (tagint) ubuf(buf[i][0]).i,(int) ubuf(buf[i][1]).i,
-            buf[i][2],(int) ubuf(buf[i][3]).i,buf[i][4],
-            (int) ubuf(buf[i][5]).i,buf[i][6],buf[i][8],
-            buf[i][8],buf[i][9],buf[i][10],
-            (int) ubuf(buf[i][11]).i,(int) ubuf(buf[i][12]).i,
-            (int) ubuf(buf[i][13]).i);
-}
-
-/* ----------------------------------------------------------------------
-   write hybrid atom info to data file
-------------------------------------------------------------------------- */
-
-int AtomVecWavepacket::write_data_hybrid(FILE *fp, double *buf)
-{
-  fprintf(fp," %-1.16e %d %-1.16e %d %-1.16e %-1.16e",
-          buf[0],(int) ubuf(buf[1]).i,buf[2],(int) ubuf(buf[3]).i,
-          buf[4],buf[5]);
-  return 6;
-}
-
-/* ----------------------------------------------------------------------
-   pack velocity info for data file
-------------------------------------------------------------------------- */
-
-void AtomVecWavepacket::pack_vel(double **buf)
-{
-  int nlocal = atom->nlocal;
-  for (int i = 0; i < nlocal; i++) {
-    buf[i][0] = ubuf(tag[i]).d;
-    buf[i][1] = v[i][0];
-    buf[i][2] = v[i][1];
-    buf[i][3] = v[i][2];
-    buf[i][4] = ervel[i];
-  }
-}
-
-/* ----------------------------------------------------------------------
-   pack hybrid velocity info for data file
-------------------------------------------------------------------------- */
-
-int AtomVecWavepacket::pack_vel_hybrid(int i, double *buf)
-{
-  buf[0] = ervel[i];
-  return 1;
-}
-
-/* ----------------------------------------------------------------------
-   write velocity info to data file
-------------------------------------------------------------------------- */
-
-void AtomVecWavepacket::write_vel(FILE *fp, int n, double **buf)
-{
-  for (int i = 0; i < n; i++)
-    fprintf(fp,TAGINT_FORMAT " %-1.16e %-1.16e %-1.16e %-1.16e\n",
-            (tagint) ubuf(buf[i][0]).i,buf[i][1],buf[i][2],buf[i][3],buf[i][4]);
-}
-
-/* ----------------------------------------------------------------------
-   write hybrid velocity info to data file
-------------------------------------------------------------------------- */
-
-int AtomVecWavepacket::write_vel_hybrid(FILE *fp, double *buf)
-{
-  fprintf(fp," %-1.16e",buf[0]);
-  return 1;
+  atom->ervel[ilocal] = 0.0;
 }
 
 /* ----------------------------------------------------------------------
@@ -1131,27 +115,31 @@ void AtomVecWavepacket::pack_property_atom(int index, double *buf,
 {
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
-  int n = 0;
 
+  int n = 0;
   if (index == 0) {
+    int *spin = atom->spin;
     for (int i = 0; i < nlocal; i++) {
       if (mask[i] & groupbit) buf[n] = spin[i];
       else buf[n] = 0.0;
       n += nvalues;
     }
   } else if (index == 1) {
+    double *eradius = atom->eradius;
     for (int i = 0; i < nlocal; i++) {
       if (mask[i] & groupbit) buf[n] = eradius[i];
       else buf[n] = 0.0;
       n += nvalues;
     }
   } else if (index == 2) {
+    double *ervel = atom->ervel;
     for (int i = 0; i < nlocal; i++) {
       if (mask[i] & groupbit) buf[n] = ervel[i];
       else buf[n] = 0.0;
       n += nvalues;
     }
   } else if (index == 3) {
+    double *erforce = atom->erforce;
     for (int i = 0; i < nlocal; i++) {
       if (mask[i] & groupbit) buf[n] = erforce[i];
       else buf[n] = 0.0;
@@ -1159,35 +147,3 @@ void AtomVecWavepacket::pack_property_atom(int index, double *buf,
     }
   }
 }
-
-/* ----------------------------------------------------------------------
-   return # of bytes of allocated memory
-------------------------------------------------------------------------- */
-
-bigint AtomVecWavepacket::memory_usage()
-{
-  bigint bytes = 0;
-
-  if (atom->memcheck("tag")) bytes += memory->usage(tag,nmax);
-  if (atom->memcheck("type")) bytes += memory->usage(type,nmax);
-  if (atom->memcheck("mask")) bytes += memory->usage(mask,nmax);
-  if (atom->memcheck("image")) bytes += memory->usage(image,nmax);
-  if (atom->memcheck("x")) bytes += memory->usage(x,nmax,3);
-  if (atom->memcheck("v")) bytes += memory->usage(v,nmax,3);
-  if (atom->memcheck("f")) bytes += memory->usage(f,nmax*comm->nthreads,3);
-
-  if (atom->memcheck("q")) bytes += memory->usage(q,nmax);
-  if (atom->memcheck("spin")) bytes += memory->usage(spin,nmax);
-  if (atom->memcheck("eradius")) bytes += memory->usage(eradius,nmax);
-  if (atom->memcheck("ervel")) bytes += memory->usage(ervel,nmax);
-  if (atom->memcheck("erforce"))
-    bytes += memory->usage(erforce,nmax*comm->nthreads);
-
-  if (atom->memcheck("ervelforce")) bytes += memory->usage(ervelforce,nmax);
-  if (atom->memcheck("cs")) bytes += memory->usage(cs,2*nmax);
-  if (atom->memcheck("csforce")) bytes += memory->usage(csforce,2*nmax);
-  if (atom->memcheck("vforce")) bytes += memory->usage(vforce,3*nmax);
-  if (atom->memcheck("etag")) bytes += memory->usage(etag,nmax);
-
-  return bytes;
-}
diff --git a/src/USER-AWPMD/atom_vec_wavepacket.h b/src/USER-AWPMD/atom_vec_wavepacket.h
index d1a0c7c7f2..e7db15db14 100644
--- a/src/USER-AWPMD/atom_vec_wavepacket.h
+++ b/src/USER-AWPMD/atom_vec_wavepacket.h
@@ -11,11 +11,6 @@
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
-/* ----------------------------------------------------------------------
-   Contributing author: Ilya Valuev (JIHT RAS)
-------------------------------------------------------------------------- */
-
-
 #ifdef ATOM_CLASS
 
 AtomStyle(wavepacket,AtomVecWavepacket)
@@ -32,77 +27,11 @@ namespace LAMMPS_NS {
 class AtomVecWavepacket : public AtomVec {
 public:
   AtomVecWavepacket(class LAMMPS *);
-  ~AtomVecWavepacket() {}
-  void grow(int);
-  void grow_reset();
-  void copy(int, int, int);
   void force_clear(int, size_t);
-  int pack_comm(int, int *, double *, int, int *);
-  int pack_comm_vel(int, int *, double *, int, int *);
-  int pack_comm_hybrid(int, int *, double *);
-  void unpack_comm(int, int, double *);
-  void unpack_comm_vel(int, int, double *);
-  int unpack_comm_hybrid(int, int, double *);
-  int pack_reverse(int, int, double *);
-  int pack_reverse_hybrid(int, int, double *);
-  void unpack_reverse(int, int *, double *);
-  int unpack_reverse_hybrid(int, int *, double *);
-  int pack_border(int, int *, double *, int, int *);
-  int pack_border_vel(int, int *, double *, int, int *);
-  int pack_border_hybrid(int, int *, double *);
-  void unpack_border(int, int, double *);
-  void unpack_border_vel(int, int, double *);
-  int unpack_border_hybrid(int, int, double *);
-  int pack_exchange(int, double *);
-  int unpack_exchange(double *);
-  int size_restart();
-  int pack_restart(int, double *);
-  int unpack_restart(double *);
-  void create_atom(int, double *);
-  void data_atom(double *, imageint, char **);
-  int data_atom_hybrid(int, char **);
-  void data_vel(int, char **);
-  int data_vel_hybrid(int, char **);
-  void pack_data(double **);
-  int pack_data_hybrid(int, double *);
-  void write_data(FILE *, int, double **);
-  int write_data_hybrid(FILE *, double *);
-  void pack_vel(double **);
-  int pack_vel_hybrid(int, double *);
-  void write_vel(FILE *, int, double **);
-  int write_vel_hybrid(FILE *, double *);
+  void create_atom_post(int);
+  void data_atom_post(int);
   int property_atom(char *);
   void pack_property_atom(int, double *, int, int);
-  bigint memory_usage();
-
-private:
-  tagint *tag;
-  int *type,*mask;
-  imageint *image;
-  double **x,**v,**f;
-
-  ///\en spin: -1 or 1 for electron, 0 for ion (compatible with eff)
-  int *spin;
-  ///\en charge: must be specified in the corresponding units (-1 for electron in real units, eff compatible)
-  double *q;
-  ///\en width of the wavepacket (compatible with eff)
-  double *eradius;
-  ///\en width velocity for the wavepacket (compatible with eff)
-  double *ervel;
-  ///\en (generalized) force on width  (compatible with eff)
-  double *erforce;
-
-  // AWPMD- specific:
-  ///\en electron tag: must be the same for the WPs belonging to the same electron
-  int *etag;
-  ///\en wavepacket split coefficients: cre, cim, size is 2*N
-  double *cs;
-  ///\en force on wavepacket split coefficients: re, im, size is 2*N
-  double *csforce;
-  ///\en (generalized) force on velocity, size is 3*N
-  double *vforce;
-   ///\en (generalized) force on radius velocity, size is N
-  double *ervelforce;
 };
 
 }
diff --git a/src/USER-AWPMD/fix_nve_awpmd.cpp b/src/USER-AWPMD/fix_nve_awpmd.cpp
index b4a1cbf72a..2aa2e7680b 100644
--- a/src/USER-AWPMD/fix_nve_awpmd.cpp
+++ b/src/USER-AWPMD/fix_nve_awpmd.cpp
@@ -74,8 +74,6 @@ void FixNVEAwpmd::init()
 
 void FixNVEAwpmd::initial_integrate(int /* vflag */)
 {
-
-
   // update v,vr and x,radius of atoms in group
 
   double **x = atom->x;
@@ -84,7 +82,7 @@ void FixNVEAwpmd::initial_integrate(int /* vflag */)
   double *ervel = atom->ervel;
   double **f = atom->f;
   double *erforce = atom->erforce;
-  double *vforce=atom->vforce;
+  double **vforce=atom->vforce;
   double *ervelforce=atom->ervelforce;
 
   double *mass = atom->mass;
@@ -101,7 +99,7 @@ void FixNVEAwpmd::initial_integrate(int /* vflag */)
       double dtfm = dtf / mass[type[i]];
       double dtfmr=dtfm;
       for(int j=0;j<3;j++){
-        x[i][j] += dtv*vforce[3*i+j];
+        x[i][j] += dtv*vforce[i][j];
         v[i][j] += dtfm*f[i][j];
       }
       eradius[i]+= dtv*ervelforce[i];
diff --git a/src/USER-DPD/atom_vec_dpd.cpp b/src/USER-DPD/atom_vec_dpd.cpp
index ce35178d8d..124a081191 100644
--- a/src/USER-DPD/atom_vec_dpd.cpp
+++ b/src/USER-DPD/atom_vec_dpd.cpp
@@ -47,7 +47,7 @@ AtomVecDPD::AtomVecDPD(LAMMPS *lmp) : AtomVec(lmp)
   fields_restart = (char *) "dpdTheta uCond uMech uChem";
   fields_create = (char *) "rho dpdTheta uCond uMech uChem uCG uCGnew duChem";
   fields_data_atom = (char *) "id type dpdTheta x";
-  fields_data_vel = (char *) "id v omega";
+  fields_data_vel = (char *) "id v";
 
   setup_fields();
 }
diff --git a/src/USER-EFF/atom_vec_electron.cpp b/src/USER-EFF/atom_vec_electron.cpp
index 3688f7f582..0ac56c48a3 100644
--- a/src/USER-EFF/atom_vec_electron.cpp
+++ b/src/USER-EFF/atom_vec_electron.cpp
@@ -44,8 +44,6 @@ AtomVecElectron::AtomVecElectron(LAMMPS *lmp) : AtomVec(lmp)
   molecular = 0;
   forceclearflag = 1;
 
-  atom->ecp_flag = 0;
-
   atom->electron_flag = 1;
   atom->q_flag = atom->spin_flag = atom->eradius_flag =
     atom->ervel_flag = atom->erforce_flag = 1;
@@ -99,7 +97,6 @@ void AtomVecElectron::create_atom_post(int ilocal)
 void AtomVecElectron::data_atom_post(int ilocal)
 {
   atom->ervel[ilocal] = 0.0;
-  if (atom->spin[ilocal] == 3) atom->ecp_flag = 1;
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/USER-EFF/pair_eff_cut.cpp b/src/USER-EFF/pair_eff_cut.cpp
index e7aed14030..f9333f4bec 100644
--- a/src/USER-EFF/pair_eff_cut.cpp
+++ b/src/USER-EFF/pair_eff_cut.cpp
@@ -801,7 +801,7 @@ void PairEffCut::settings(int narg, char **arg)
 
   int atype;
   int iarg = 1;
-  int ecp_found = 0;
+  ecp_found = 0;
 
   while (iarg < narg) {
     if (strcmp(arg[iarg],"limit/eradius") == 0) {
@@ -821,17 +821,15 @@ void PairEffCut::settings(int narg, char **arg)
         else if (strcmp(arg[iarg+1],"O") == 0) ecp_type[atype] = 8;
         else if (strcmp(arg[iarg+1],"Al") == 0) ecp_type[atype] = 13;
         else if (strcmp(arg[iarg+1],"Si") == 0) ecp_type[atype] = 14;
-        else error->all(FLERR, "Note: there are no default parameters for this atom ECP\n");
+        else error->all(FLERR, "No default parameters for this atom ECP\n");
         iarg += 2;
         ecp_found = 1;
       }
-    }
+    } else error->all(FLERR,"Illegal pair style command");
   }
 
-  if (!ecp_found && atom->ecp_flag)
-    error->all(FLERR,"Need to specify ECP type on pair_style command");
-
   // Need to introduce 2 new constants w/out changing update.cpp
+
   if (force->qqr2e==332.06371) {        // i.e. Real units chosen
     h2e = 627.509;                      // hartree->kcal/mol
     hhmss2e = 175.72044219620075;       // hartree->kcal/mol * (Bohr->Angstrom)^2
@@ -872,9 +870,24 @@ void PairEffCut::init_style()
 
   if (update->whichflag == 1) {
     if (force->qqr2e == 332.06371 && update->dt == 1.0)
-      error->all(FLERR,"You must lower the default real units timestep for pEFF ");
+      error->all(FLERR,"Must lower the default real units timestep for pEFF ");
   }
 
+  // check if any atom's spin = 3 and ECP type was not set
+
+  int *spin = atom->spin;
+  int nlocal = atom->nlocal;
+
+  int flag = 0;
+  for (int i = 0; i < nlocal; i++)
+    if (spin[i] == 3) flag = 1;
+
+  int flagall;
+  MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
+
+  if (flagall && !ecp_found)
+    error->all(FLERR,"Need to specify ECP type on pair_style command");
+
   // need a half neigh list and optionally a granular history neigh list
 
   neighbor->request(this,instance_me);
diff --git a/src/USER-EFF/pair_eff_cut.h b/src/USER-EFF/pair_eff_cut.h
index 63dabe5db8..bd04344373 100644
--- a/src/USER-EFF/pair_eff_cut.h
+++ b/src/USER-EFF/pair_eff_cut.h
@@ -46,10 +46,12 @@ class PairEffCut : public Pair {
 
  private:
   int limit_eradius_flag, pressure_with_evirials_flag;
+  int ecp_found;
   double cut_global;
   double **cut;
   int ecp_type[100];
-  double PAULI_CORE_A[100], PAULI_CORE_B[100], PAULI_CORE_C[100], PAULI_CORE_D[100], PAULI_CORE_E[100];
+  double PAULI_CORE_A[100],PAULI_CORE_B[100],PAULI_CORE_C[100];
+  double PAULI_CORE_D[100],PAULI_CORE_E[100];
   double hhmss2e, h2e;
 
   int nmax;
diff --git a/src/USER-MESO/atom_vec_edpd.cpp b/src/USER-MESO/atom_vec_edpd.cpp
index 7a031c256b..e06ac633ec 100644
--- a/src/USER-MESO/atom_vec_edpd.cpp
+++ b/src/USER-MESO/atom_vec_edpd.cpp
@@ -41,16 +41,16 @@ AtomVecEDPD::AtomVecEDPD(LAMMPS *lmp) : AtomVec(lmp)
   // order of fields in a string does not matter
   // except: fields_data_atom & fields_data_vel must match data file
 
-  fields_grow = (char *) "edpd_cv edpd_temp edpd_flux vest";
-  fields_copy = (char *) "edpd_cv edpd_temp edpd_flux vest";
-  fields_comm = (char *) "edpd_temp vest";
-  fields_comm_vel = (char *) "edpd_temp vest";
+  fields_grow = (char *) "edpd_cv edpd_temp edpd_flux vest vest_temp";
+  fields_copy = (char *) "edpd_cv edpd_temp edpd_flux vest vest_temp";
+  fields_comm = (char *) "edpd_temp vest vest_temp";
+  fields_comm_vel = (char *) "edpd_temp vest vest_temp";
   fields_reverse = (char *) "edpd_flux";
-  fields_border = (char *) "edpd_cv edpd_temp vest";
-  fields_border_vel = (char *) "edpd_cv edpd_temp vest";
-  fields_exchange = (char *) "edpd_cv edpd_temp vest";
-  fields_restart = (char * ) "edpd_cv edpd_temp vest";
-  fields_create = (char *) "edpd_cv edpd_temp edpd_flux vest";
+  fields_border = (char *) "edpd_cv edpd_temp vest vest_temp";
+  fields_border_vel = (char *) "edpd_cv edpd_temp vest vest_temp";
+  fields_exchange = (char *) "edpd_cv edpd_temp vest vest_temp";
+  fields_restart = (char * ) "edpd_cv edpd_temp vest vest_temp";
+  fields_create = (char *) "edpd_cv edpd_temp edpd_flux vest vest_temp";
   fields_data_atom = (char *) "id type edpd_temp edpd_cv x";
   fields_data_vel = (char *) "id v";
 
@@ -85,7 +85,7 @@ void AtomVecEDPD::create_atom_post(int ilocal)
 {
   atom->edpd_temp[ilocal] = 1.0;
   atom->edpd_cv[ilocal]= 1.0e5;
-  atom->vest[ilocal][3] = atom->edpd_temp[ilocal];
+  atom->vest_temp[ilocal] = atom->edpd_temp[ilocal];
 }
 
 /* ----------------------------------------------------------------------
@@ -99,5 +99,5 @@ void AtomVecEDPD::data_atom_post(int ilocal)
   atom->vest[ilocal][0] = 0.0;
   atom->vest[ilocal][1] = 0.0;
   atom->vest[ilocal][2] = 0.0;
-  atom->vest[ilocal][3] = atom->edpd_temp[ilocal];
+  atom->vest_temp[ilocal] = atom->edpd_temp[ilocal];
 }
diff --git a/src/USER-MESO/atom_vec_tdpd.cpp b/src/USER-MESO/atom_vec_tdpd.cpp
index 55284f69a2..5734fcf9ad 100644
--- a/src/USER-MESO/atom_vec_tdpd.cpp
+++ b/src/USER-MESO/atom_vec_tdpd.cpp
@@ -63,7 +63,7 @@ void AtomVecTDPD::process_args(int narg, char **arg)
   cc_species = atom->cc_species;
 
   atom->add_peratom_change_columns("cc",cc_species);
-  atom->add_peratom_change_columns("cc_species",cc_species);
+  atom->add_peratom_change_columns("cc_flux",cc_species);
 
   // delay setting up of fields until now
 
diff --git a/src/USER-MESO/fix_mvv_edpd.cpp b/src/USER-MESO/fix_mvv_edpd.cpp
index bd9cd9cc2a..3294d8d682 100644
--- a/src/USER-MESO/fix_mvv_edpd.cpp
+++ b/src/USER-MESO/fix_mvv_edpd.cpp
@@ -88,6 +88,7 @@ void FixMvvEDPD::initial_integrate(int /*vflag*/)
   double *edpd_flux = atom->edpd_flux;
   double *edpd_cv = atom->edpd_cv;
   double **vest = atom->vest;
+  double *vest_temp = atom->vest_temp;
   double *rmass = atom->rmass;
   double *mass = atom->mass;
   int *type = atom->type;
@@ -105,7 +106,7 @@ void FixMvvEDPD::initial_integrate(int /*vflag*/)
      vest[i][0] = v[i][0] + dtfm * f[i][0];
      vest[i][1] = v[i][1] + dtfm * f[i][1];
      vest[i][2] = v[i][2] + dtfm * f[i][2];
-     vest[i][3] = edpd_temp[i] + dtT * edpd_flux[i];
+     vest_temp[i] = edpd_temp[i] + dtT * edpd_flux[i];
 
      x[i][0] += dtv * vest[i][0];
      x[i][1] += dtv * vest[i][1];
@@ -131,6 +132,7 @@ void FixMvvEDPD::final_integrate()
   double *edpd_flux = atom->edpd_flux;
   double *edpd_cv = atom->edpd_cv;
   double **vest = atom->vest;
+  double *vest_temp = atom->vest_temp;
   double *rmass = atom->rmass;
   double *mass = atom->mass;
   int *type = atom->type;
@@ -148,7 +150,7 @@ void FixMvvEDPD::final_integrate()
      v[i][0] = vest[i][0] + dtfm * f[i][0];
      v[i][1] = vest[i][1] + dtfm * f[i][1];
      v[i][2] = vest[i][2] + dtfm * f[i][2];
-     edpd_temp[i] = vest[i][3] + dtT * edpd_flux[i];
+     edpd_temp[i] = vest_temp[i] + dtT * edpd_flux[i];
   }
 }
 
diff --git a/src/USER-SMD/atom_vec_smd.cpp b/src/USER-SMD/atom_vec_smd.cpp
index 0fada3dc2a..9f0ed04a09 100644
--- a/src/USER-SMD/atom_vec_smd.cpp
+++ b/src/USER-SMD/atom_vec_smd.cpp
@@ -63,7 +63,7 @@ AtomVecSMD::AtomVecSMD(LAMMPS *lmp) : AtomVec(lmp)
 
   fields_grow = (char *) 
     "de vfrac rmass x0 radius contact_radius molecule "
-    "smd_data_9 e vest tlsph_stress "
+    "smd_data_9 e vest smd_stress "
     "eff_plastic_strain eff_plastic_strain_rate damage";
   fields_copy = (char *) 
     "vfrac rmass x0 radius contact_radius molecule e "
@@ -91,7 +91,7 @@ AtomVecSMD::AtomVecSMD(LAMMPS *lmp) : AtomVec(lmp)
     "eff_plastic_strain eff_plastic_strain_rate smd_data_9 smd_stress damage";
   fields_data_atom = (char *) 
     "id type molecule vfrac rmass radius contact_radius x";
-  fields_data_vel = (char *) "id v vest";
+  fields_data_vel = (char *) "id v";
 
   // set these array sizes based on defines
 
diff --git a/src/atom.cpp b/src/atom.cpp
index 8fd859478f..0fbee8f583 100644
--- a/src/atom.cpp
+++ b/src/atom.cpp
@@ -566,6 +566,7 @@ void Atom::peratom_create()
 
   add_peratom("edpd_cv",&edpd_cv,DOUBLE,0);
   add_peratom("edpd_temp",&edpd_temp,DOUBLE,0);
+  add_peratom("vest_temp",&vest_temp,DOUBLE,0);
   add_peratom("edpd_flux",&edpd_flux,DOUBLE,0,1);     // set per-thread flag
   add_peratom("cc",&cc,DOUBLE,1);
   add_peratom("cc_flux",&cc_flux,DOUBLE,1,1);         // set per-thread flag
@@ -625,8 +626,11 @@ void Atom::add_peratom(const char *name, void *address,
 
 void Atom::add_peratom_change_columns(const char *name, int cols)
 {
+  int i;
   for (int i = 0; i < nperatom; i++)
     if (strcmp(name,peratom[i].name) == 0) peratom[i].cols = cols;
+  if (i == nperatom) 
+    error->all(FLERR,"Could not find name of peratom array for column change");
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/atom.h b/src/atom.h
index 007d364a7b..bc1eb1a7d7 100644
--- a/src/atom.h
+++ b/src/atom.h
@@ -124,6 +124,7 @@ class Atom : protected Pointers {
 
   double **cc,**cc_flux;         // cc = chemical concentration
   double *edpd_temp,*edpd_flux;  // temperature and heat flux
+  double *vest_temp;
   double *edpd_cv;               // heat capacity
   int cc_species;
 
@@ -162,10 +163,6 @@ class Atom : protected Pointers {
 
   int sp_flag;
 
-  // USER-EFF package
-
-  int ecp_flag;
-
   // USER-SMD package
 
   int smd_flag;
diff --git a/src/atom_vec.cpp b/src/atom_vec.cpp
index bbfccadb89..e31c235760 100644
--- a/src/atom_vec.cpp
+++ b/src/atom_vec.cpp
@@ -45,6 +45,7 @@ AtomVec::AtomVec(LAMMPS *lmp) : Pointers(lmp)
   forceclearflag = 0;
   maxexchange = 0;
   bonus_flag = 0;
+  size_forward_bonus = size_border_bonus = 0;
 
   kokkosable = 0;
 
@@ -176,7 +177,7 @@ int AtomVec::grow_nmax_bonus(int nmax_bonus)
 
 void AtomVec::grow(int n)
 {
-  int i,datatype,cols,maxcols;
+  int datatype,cols,maxcols;
   void *pdata,*plength;
 
   if (n == 0) grow_nmax();
@@ -193,7 +194,7 @@ void AtomVec::grow(int n)
   v = memory->grow(atom->v,nmax,3,"atom:v");
   f = memory->grow(atom->f,nmax*nthreads,3,"atom:f");
 
-  for (i = 0; i < ngrow; i++) {
+  for (int i = 0; i < ngrow; i++) {
     pdata = mgrow.pdata[i];
     datatype = mgrow.datatype[i];
     cols = mgrow.cols[i];
@@ -415,7 +416,7 @@ int AtomVec::pack_comm(int n, int *list, double *buf,
 /* ---------------------------------------------------------------------- */
 
 int AtomVec::pack_comm_vel(int n, int *list, double *buf,
-                                 int pbc_flag, int *pbc)
+                           int pbc_flag, int *pbc)
 {
   int i,j,m,mm,nn,datatype,cols;
   double dx,dy,dz,dvx,dvy,dvz;
@@ -1375,15 +1376,15 @@ int AtomVec::size_restart()
 
   if (nrestart) {
     for (nn = 0; nn < nrestart; nn++) {
-      cols = mrestart.cols[i];
+      cols = mrestart.cols[nn];
       if (cols == 0) n += nlocal;
       else if (cols > 0) n += cols*nlocal;
       else {
         collength = mrestart.collength[nn];
         plength = mrestart.plength[nn];
         for (i = 0; i < nlocal; i++) {
-          if (collength) ncols = (*((int ***) plength))[nlocal][collength-1];
-          else ncols = (*((int **) plength))[nlocal];
+          if (collength) ncols = (*((int ***) plength))[i][collength-1];
+          else ncols = (*((int **) plength))[i];
           n += ncols;
         }
       }
@@ -1435,7 +1436,7 @@ int AtomVec::pack_restart(int i, double *buf)
       if (cols == 0) {
         double *vec = *((double **) pdata);
         buf[m++] = vec[i];
-      } else if (ncols > 0) {
+      } else if (cols > 0) {
         double **array = *((double ***) pdata);
         for (mm = 0; mm < cols; mm++)
           buf[m++] = array[i][mm];
@@ -1469,7 +1470,7 @@ int AtomVec::pack_restart(int i, double *buf)
       if (cols == 0) {
         bigint *vec = *((bigint **) pdata);
         buf[m++] = ubuf(vec[i]).d;
-      } else if (ncols > 0) {
+      } else if (cols > 0) {
         bigint **array = *((bigint ***) pdata);
         for (mm = 0; mm < cols; mm++)
           buf[m++] = ubuf(array[i][mm]).d;
@@ -1705,7 +1706,7 @@ void AtomVec::data_atom(double *coord, imageint imagetmp, char **values)
         vec[nlocal] = utils::numeric(FLERR,values[ivalue++],true,lmp);
       } else {
         double **array = *((double ***) pdata);
-        if (array == atom->x) {      // already set by coord arg
+        if (array == atom->x) {      // x was already set by coord arg
           ivalue += cols;
           continue;
         }
@@ -1878,9 +1879,9 @@ void AtomVec::data_vel(int ilocal, char **values)
   v[ilocal][1] = utils::numeric(FLERR,values[1],true,lmp);
   v[ilocal][2] = utils::numeric(FLERR,values[2],true,lmp);
 
-  if (ndata_vel) {
+  if (ndata_vel > 2) {
     int ivalue = 3;
-    for (n = 0; n < ndata_vel; n++) {
+    for (n = 2; n < ndata_vel; n++) {
       pdata = mdata_vel.pdata[n];
       datatype = mdata_vel.datatype[n];
       cols = mdata_vel.cols[n];
@@ -2434,6 +2435,7 @@ void AtomVec::setup_fields()
     else size_data_atom += cols;
   }
 
+
   size_data_vel = 0;
   for (n = 0; n < ndata_vel; n++) {
     cols = mdata_vel.cols[n];
diff --git a/src/atom_vec.h b/src/atom_vec.h
index 69ea1d51ab..5911c8414b 100644
--- a/src/atom_vec.h
+++ b/src/atom_vec.h
@@ -74,7 +74,7 @@ class AtomVec : protected Pointers {
 
   virtual void force_clear(int, size_t) {}
 
-  void grow(int);
+  virtual void grow(int);
   void copy(int, int, int);
 
   virtual void copy_bonus(int, int, int) {}
@@ -122,13 +122,11 @@ class AtomVec : protected Pointers {
 
   void data_atom(double *, imageint, char **);
   virtual void data_atom_post(int) {}
-
-  void data_atom_bonus(int, char **) {}
-  void data_body(int, int, int, int *, double *) {}
+  virtual void data_atom_bonus(int, char **) {}
+  virtual void data_body(int, int, int, int *, double *) {}
 
   void pack_data(double **);
   void write_data(FILE *, int, double **);
-
   virtual void pack_data_pre(int) {}
   virtual void pack_data_post(int) {}
 
@@ -145,19 +143,19 @@ class AtomVec : protected Pointers {
   int pack_improper(tagint **);
   void write_improper(FILE *, int, tagint **, int);
 
-  int property_atom(char *) {return -1;}
-  void pack_property_atom(int, double *, int, int) {}
+  virtual int property_atom(char *) {return -1;}
+  virtual void pack_property_atom(int, double *, int, int) {}
 
   bigint memory_usage();
   virtual bigint memory_usage_bonus() {}
 
  protected:
-  int nmax;                             // local copy of atom->nmax
-  int deform_vremap;                    // local copy of domain properties
+  int nmax;                    // local copy of atom->nmax
+  int deform_vremap;           // local copy of domain properties
   int deform_groupbit;
   double *h_rate;
 
-  tagint *tag;                          // peratom fields common to all styles
+  tagint *tag;                 // peratom fields common to all styles
   int *type,*mask;
   imageint *image;
   double **x,**v,**f;
diff --git a/src/atom_vec_body.cpp b/src/atom_vec_body.cpp
index 19ac69cef3..a04f23a47c 100644
--- a/src/atom_vec_body.cpp
+++ b/src/atom_vec_body.cpp
@@ -40,8 +40,8 @@ AtomVecBody::AtomVecBody(LAMMPS *lmp) : AtomVec(lmp)
   // size_data_bonus is not used by Atom class for body style
 
   size_forward_bonus = 4;
-  size_border_bonus = 9;
-  size_restart_bonus_one = 9;
+  size_border_bonus = 10;
+  size_restart_bonus_one = 10;
   size_data_bonus = 0;
 
   atom->body_flag = 1;
@@ -72,11 +72,9 @@ AtomVecBody::AtomVecBody(LAMMPS *lmp) : AtomVec(lmp)
   fields_border_vel = (char *) "radius rmass angmom";
   fields_exchange = (char *) "radius rmass angmom";
   fields_restart = (char *) "radius rmass angmom";
-  fields_create = (char *) "radius rmass angmom tri";
+  fields_create = (char *) "radius rmass angmom body";
   fields_data_atom = (char *) "id type body rmass x";
   fields_data_vel = (char *) "id v angmom";
-
-  setup_fields();
 }
 
 /* ---------------------------------------------------------------------- */
@@ -127,6 +125,20 @@ void AtomVecBody::process_args(int narg, char **arg)
 
   size_forward_bonus += bptr->size_forward;
   size_border_bonus += bptr->size_border;
+
+  setup_fields();
+}
+
+/* ----------------------------------------------------------------------
+   grow atom arrays
+   must set local copy of body ptr
+   needed in replicate when 2 atom classes exist and pack_restart() is called
+------------------------------------------------------------------------- */
+
+void AtomVecBody::grow(int n)
+{
+  AtomVec::grow(n);
+  body = atom->body;
 }
 
 /* ----------------------------------------------------------------------
@@ -150,8 +162,6 @@ void AtomVecBody::grow_bonus()
 
 void AtomVecBody::copy_bonus(int i, int j, int delflag)
 {
-  int *body = atom->body;
-
   // if deleting atom J via delflag and J has bonus data, then delete it
 
   if (delflag && body[j] >= 0) {
@@ -206,8 +216,6 @@ int AtomVecBody::pack_comm_bonus(int n, int *list, double *buf)
   int i,j,m;
   double *quat;
 
-  int *body = atom->body;
-
   m = 0;
   for (i = 0; i < n; i++) {
     j = list[i];
@@ -231,8 +239,6 @@ void AtomVecBody::unpack_comm_bonus(int n, int first, double *buf)
   int i,m,last;
   double *quat;
 
-  int *body = atom->body;
-
   m = 0;
   last = first + n;
   for (i = first; i < last; i++) {
@@ -254,8 +260,6 @@ int AtomVecBody::pack_border_bonus(int n, int *list, double *buf)
   int i,j,m;
   double *quat,*inertia;
 
-  int *body = atom->body;
-
   m = 0;
   for (i = 0; i < n; i++) {
     j = list[i];
@@ -287,8 +291,6 @@ int AtomVecBody::unpack_border_bonus(int n, int first, double *buf)
   int i,j,m,last;
   double *quat,*inertia;
 
-  int *body = atom->body;
-
   m = 0;
   last = first + n;
   for (i = first; i < last; i++) {
@@ -330,8 +332,6 @@ int AtomVecBody::pack_exchange_bonus(int i, double *buf)
 {
   int m = 0;
 
-  int *body = atom->body;
-
   if (body[i] < 0) buf[m++] = ubuf(0).d;
   else {
     buf[m++] = ubuf(1).d;
@@ -363,8 +363,6 @@ int AtomVecBody::unpack_exchange_bonus(int ilocal, double *buf)
 {
   int m = 0;
 
-  int *body = atom->body;
-
   body[ilocal] = (int) ubuf(buf[m++]).i;
   if (body[ilocal] == 0) body[ilocal] = -1;
   else {
@@ -408,8 +406,6 @@ int AtomVecBody::size_restart_bonus()
 {
   int i;
 
-  int *body = atom->body;
-
   int n = 0;
   int nlocal = atom->nlocal;
   for (i = 0; i < nlocal; i++) {
@@ -418,8 +414,7 @@ int AtomVecBody::size_restart_bonus()
       if (intdoubleratio == 1) n += bonus[body[i]].ninteger;
       else n += (bonus[body[i]].ninteger+1)/2;
       n += bonus[body[i]].ndouble;
-    }
-    n++;
+    } else n++;
   }
 
   return n;
@@ -435,8 +430,6 @@ int AtomVecBody::pack_restart_bonus(int i, double *buf)
 {
   int m = 0;
 
-  int *body = atom->body;
-
   if (body[i] < 0) buf[m++] = ubuf(0).d;
   else {
     buf[m++] = ubuf(1).d;
@@ -470,8 +463,6 @@ int AtomVecBody::unpack_restart_bonus(int ilocal, double *buf)
 {
   int m = 0;
 
-  int *body = atom->body;
-
   body[ilocal] = (int) ubuf(buf[m++]).i;
   if (body[ilocal] == 0) body[ilocal] = -1;
   else {
diff --git a/src/atom_vec_body.h b/src/atom_vec_body.h
index 939b01878d..3f32c8223c 100644
--- a/src/atom_vec_body.h
+++ b/src/atom_vec_body.h
@@ -43,6 +43,7 @@ class AtomVecBody : public AtomVec {
   ~AtomVecBody();
   void process_args(int, char **);
 
+  void grow(int);
   void copy_bonus(int, int, int);
   void clear_bonus();
   int pack_comm_bonus(int, int *, double *);
@@ -74,6 +75,8 @@ class AtomVecBody : public AtomVec {
   int intdoubleratio;            // sizeof(double) / sizeof(int)
   int body_flag;
 
+  int *body;
+
   MyPoolChunk<int> *icp;
   MyPoolChunk<double> *dcp;
 
diff --git a/src/atom_vec_ellipsoid.cpp b/src/atom_vec_ellipsoid.cpp
index 9166293384..e8ab6d5613 100644
--- a/src/atom_vec_ellipsoid.cpp
+++ b/src/atom_vec_ellipsoid.cpp
@@ -38,7 +38,7 @@ AtomVecEllipsoid::AtomVecEllipsoid(LAMMPS *lmp) : AtomVec(lmp)
 
   size_forward_bonus = 4;
   size_border_bonus = 8;
-  size_restart_bonus_one = 7;
+  size_restart_bonus_one = 8;
   size_data_bonus = 8;
 
   atom->ellipsoid_flag = 1;
@@ -320,7 +320,7 @@ int AtomVecEllipsoid::size_restart_bonus()
   int nlocal = atom->nlocal;
   for (i = 0; i < nlocal; i++) {
     if (ellipsoid[i] >= 0) n += size_restart_bonus_one;
-    n++;
+    else n++;
   }
 
   return n;
diff --git a/src/atom_vec_line.cpp b/src/atom_vec_line.cpp
index ded2f88c2f..7ab697b349 100644
--- a/src/atom_vec_line.cpp
+++ b/src/atom_vec_line.cpp
@@ -37,7 +37,7 @@ AtomVecLine::AtomVecLine(LAMMPS *lmp) : AtomVec(lmp)
 
   size_forward_bonus = 1;
   size_border_bonus = 3;
-  size_restart_bonus_one = 2;
+  size_restart_bonus_one = 3;
   size_data_bonus = 5;
 
   atom->line_flag = 1;
@@ -286,7 +286,7 @@ int AtomVecLine::size_restart_bonus()
   int nlocal = atom->nlocal;
   for (i = 0; i < nlocal; i++) {
     if (line[i] >= 0) n += size_restart_bonus_one;
-    n++;
+    else n++;
   }
 
   return n;
diff --git a/src/atom_vec_sphere.cpp b/src/atom_vec_sphere.cpp
index 81f1b02fdb..b601fc8f7c 100644
--- a/src/atom_vec_sphere.cpp
+++ b/src/atom_vec_sphere.cpp
@@ -61,18 +61,22 @@ AtomVecSphere::AtomVecSphere(LAMMPS *lmp) : AtomVec(lmp)
 
 void AtomVecSphere::process_args(int narg, char **arg)
 {
-  if (narg == 0) return;
-  if (narg != 1) error->all(FLERR,"Illegal atom_style sphere command");
-
-  radvary = utils::numeric(FLERR,arg[0],true,lmp);
-  if (radvary < 0 || radvary > 1)
+  if (narg != 0 && narg != 1) 
     error->all(FLERR,"Illegal atom_style sphere command");
-  if (radvary == 0) return;
+
+  radvary = 0;
+  if (narg == 1) {
+    radvary = utils::numeric(FLERR,arg[0],true,lmp);
+    if (radvary < 0 || radvary > 1)
+      error->all(FLERR,"Illegal atom_style sphere command");
+  }
 
   // dynamic particle radius and mass must be communicated every step
 
-  fields_comm = (char *) "radius rmass";
-  fields_comm_vel = (char *) "radius rmass omega";
+  if (radvary) {
+    fields_comm = (char *) "radius rmass";
+    fields_comm_vel = (char *) "radius rmass omega";
+  }
 
   // delay setting up of fields until now
 
diff --git a/src/atom_vec_tri.cpp b/src/atom_vec_tri.cpp
index 32d75c16f3..7345646a1e 100644
--- a/src/atom_vec_tri.cpp
+++ b/src/atom_vec_tri.cpp
@@ -38,7 +38,7 @@ AtomVecTri::AtomVecTri(LAMMPS *lmp) : AtomVec(lmp)
 
   size_forward_bonus = 4;
   size_border_bonus = 17;
-  size_restart_bonus_one = 16;
+  size_restart_bonus_one = 17;
   size_data_bonus = 10;
 
   atom->tri_flag = 1;
@@ -64,7 +64,7 @@ AtomVecTri::AtomVecTri(LAMMPS *lmp) : AtomVec(lmp)
   fields_border_vel = (char *) "molecule radius rmass omega";
   fields_exchange = (char *) "molecule radius rmass omega angmom";
   fields_restart = (char *) "molecule radius rmass omega angmom";
-  fields_create = (char *) "molecule radius rmass omega angmom line";
+  fields_create = (char *) "molecule radius rmass omega angmom tri";
   fields_data_atom = (char *) "id molecule type tri rmass x";
   fields_data_vel = (char *) "id v omega angmom";
 
@@ -381,7 +381,7 @@ int AtomVecTri::size_restart_bonus()
   int nlocal = atom->nlocal;
   for (i = 0; i < nlocal; i++) {
     if (tri[i] >= 0) n += size_restart_bonus_one;
-    n++;
+    else n++;
   }
 
   return n;
diff --git a/src/fix_gravity.cpp b/src/fix_gravity.cpp
index 14ba913c01..e00bbb4e17 100644
--- a/src/fix_gravity.cpp
+++ b/src/fix_gravity.cpp
@@ -31,7 +31,7 @@ using namespace FixConst;
 using namespace MathConst;
 
 enum{CHUTE,SPHERICAL,VECTOR};
-enum{CONSTANT,EQUAL};
+enum{CONSTANT,EQUAL};          // same as FixPour
 
 /* ---------------------------------------------------------------------- */
 
@@ -134,6 +134,16 @@ FixGravity::FixGravity(LAMMPS *lmp, int narg, char **arg) :
 
   eflag = 0;
   egrav = 0.0;
+
+  // set gravity components once and for all if CONSTANT
+
+  varflag = CONSTANT;
+  if (mstyle != CONSTANT || vstyle != CONSTANT || pstyle != CONSTANT ||
+      tstyle != CONSTANT || xstyle != CONSTANT || ystyle != CONSTANT ||
+      zstyle != CONSTANT) varflag = EQUAL;
+
+  if (varflag == CONSTANT) set_acceleration();
+
 }
 
 /* ---------------------------------------------------------------------- */
@@ -222,15 +232,6 @@ void FixGravity::init()
     if (!input->variable->equalstyle(zvar))
       error->all(FLERR,"Variable for fix gravity is invalid style");
   }
-
-  varflag = CONSTANT;
-  if (mstyle != CONSTANT || vstyle != CONSTANT || pstyle != CONSTANT ||
-      tstyle != CONSTANT || xstyle != CONSTANT || ystyle != CONSTANT ||
-      zstyle != CONSTANT) varflag = EQUAL;
-
-  // set gravity components once and for all
-
-  if (varflag == CONSTANT) set_acceleration();
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/replicate.cpp b/src/replicate.cpp
index 1617ab0313..3659f7cf7a 100644
--- a/src/replicate.cpp
+++ b/src/replicate.cpp
@@ -331,35 +331,43 @@ void Replicate::command(int narg, char **arg)
 
     if (me == 0 && screen) {
       fprintf(screen,"  bounding box image = (%i %i %i) to (%i %i %i)\n",
-              _imagelo[0],_imagelo[1],_imagelo[2],_imagehi[0],_imagehi[1],_imagehi[2]);
+              _imagelo[0],_imagelo[1],_imagelo[2],
+              _imagehi[0],_imagehi[1],_imagehi[2]);
       fprintf(screen,"  bounding box extra memory = %.2f MB\n",
               (double)size_buf_all*sizeof(double)/1024/1024);
     }
 
     // rnk offsets
 
-    int * disp_buf_rnk;
+    int *disp_buf_rnk;
     memory->create(disp_buf_rnk, nprocs, "replicate:disp_buf_rnk");
     disp_buf_rnk[0] = 0;
-    for (i=1; i<nprocs; ++i) disp_buf_rnk[i] = disp_buf_rnk[i-1] + size_buf_rnk[i-1];
+    for (i = 1; i < nprocs; i++) 
+      disp_buf_rnk[i] = disp_buf_rnk[i-1] + size_buf_rnk[i-1];
 
     // allgather buf_all
 
     double * buf_all;
     memory->create(buf_all, size_buf_all, "replicate:buf_all");
 
-    MPI_Allgatherv(buf, n, MPI_DOUBLE, buf_all, size_buf_rnk, disp_buf_rnk, MPI_DOUBLE, world);
+    MPI_Allgatherv(buf,n,MPI_DOUBLE,buf_all,size_buf_rnk,disp_buf_rnk,
+                   MPI_DOUBLE,world);
 
     // bounding box of original unwrapped system
 
     double _orig_lo[3], _orig_hi[3];
     if (triclinic) {
-      _orig_lo[0] = domain->boxlo[0] + _imagelo[0] * old_xprd + _imagelo[1] * old_xy + _imagelo[2] * old_xz;
-      _orig_lo[1] = domain->boxlo[1] + _imagelo[1] * old_yprd + _imagelo[2] * old_yz;
+      _orig_lo[0] = domain->boxlo[0] + 
+        _imagelo[0] * old_xprd + _imagelo[1] * old_xy + _imagelo[2] * old_xz;
+      _orig_lo[1] = domain->boxlo[1] + 
+        _imagelo[1] * old_yprd + _imagelo[2] * old_yz;
       _orig_lo[2] = domain->boxlo[2] + _imagelo[2] * old_zprd;
 
-      _orig_hi[0] = domain->boxlo[0] + (_imagehi[0]+1) * old_xprd + (_imagehi[1]+1) * old_xy + (_imagehi[2]+1) * old_xz;
-      _orig_hi[1] = domain->boxlo[1] + (_imagehi[1]+1) * old_yprd + (_imagehi[2]+1) * old_yz;
+      _orig_hi[0] = domain->boxlo[0] + 
+        (_imagehi[0]+1) * old_xprd + 
+        (_imagehi[1]+1) * old_xy + (_imagehi[2]+1) * old_xz;
+      _orig_hi[1] = domain->boxlo[1] + 
+        (_imagehi[1]+1) * old_yprd + (_imagehi[2]+1) * old_yz;
       _orig_hi[2] = domain->boxlo[2] + (_imagehi[2]+1) * old_zprd;
     } else {
       _orig_lo[0] = domain->boxlo[0] + _imagelo[0] * old_xprd;
@@ -605,7 +613,8 @@ void Replicate::command(int narg, char **arg)
     MPI_Reduce(&num_replicas_added, &sum, 1, MPI_INT, MPI_SUM, 0, world);
     double avg = (double) sum / nprocs;
     if (me == 0 && screen)
-      fprintf(screen,"  average # of replicas added to proc = %.2f out of %i (%.2f %%)\n",
+      fprintf(screen,"  average # of replicas added to proc = %.2f "
+              "out of %i (%.2f %%)\n",
               avg,nx*ny*nz,avg/(nx*ny*nz)*100.0);
 
   } else {
diff --git a/src/verlet.cpp b/src/verlet.cpp
index 8cd6fe940d..fe9645618a 100644
--- a/src/verlet.cpp
+++ b/src/verlet.cpp
@@ -109,7 +109,7 @@ void Verlet::setup(int flag)
   domain->pbc();
   domain->reset_box();
   comm->setup();
-  if (neighbor->style) neighbor->setup_bins();
+  if (neighbor->style) neighbor->setup_bins(); 
   comm->exchange();
   if (atom->sortfreq > 0) atom->sort();
   comm->borders();