diff --git a/src/USER-AWPMD/atom_vec_wavepacket.cpp b/src/USER-AWPMD/atom_vec_wavepacket.cpp
index 155ad88ca..213a770ea 100644
--- a/src/USER-AWPMD/atom_vec_wavepacket.cpp
+++ b/src/USER-AWPMD/atom_vec_wavepacket.cpp
@@ -1,1134 +1,1142 @@
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
 
    Copyright (2003) Sandia Corporation.  Under the terms of Contract
    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
    certain rights in this software.  This software is distributed under
    the GNU General Public License.
 
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
 /* ----------------------------------------------------------------------
    Contributing author: Ilya Valuev (JIHT, Moscow, Russia)
 ------------------------------------------------------------------------- */
 
 #include "math.h"
 #include "stdlib.h"
+#include "string.h"
 #include "atom_vec_wavepacket.h"
 #include "atom.h"
 #include "comm.h"
 #include "domain.h"
 #include "modify.h"
 #include "force.h"
 #include "fix.h"
 #include "memory.h"
 #include "error.h"
 
 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->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;
 }
 
 /* ----------------------------------------------------------------------
    grow atom-electron arrays
    n = 0 grows arrays by a chunk
    n > 0 allocates arrays to size n
 ------------------------------------------------------------------------- */
 
 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;
 }
 
 /* ----------------------------------------------------------------------
    pack data for atom I for sending to another proc
    xyz must be 1st 3 values, so comm::exchange() can test on them
 ------------------------------------------------------------------------- */
 
 int AtomVecWavepacket::pack_exchange(int i, double *buf)
 {
   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;
 }
 
 /* ----------------------------------------------------------------------
    size of restart data for all atoms owned by this proc
    include extra data stored by fixes
 ------------------------------------------------------------------------- */
 
 int AtomVecWavepacket::size_restart()
 {
   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] = ATOTAGINT(values[0]);
   type[nlocal] = atoi(values[1]);
   if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes)
     error->one(FLERR,"Invalid atom type in Atoms section of data file");
 
   q[nlocal] = atof(values[2]);
   spin[nlocal] = atoi(values[3]);
   eradius[nlocal] = atof(values[4]);
   if (eradius[nlocal] < 0.0)
     error->one(FLERR,"Invalid eradius in Atoms section of data file");
 
   etag[nlocal] = atoi(values[5]);
   cs[2*nlocal] = atoi(values[6]);
   cs[2*nlocal+1] = atof(values[7]);
 
   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] = atof(values[0]);
   spin[nlocal] = atoi(values[1]);
   eradius[nlocal] = atof(values[2]);
   if (eradius[nlocal] < 0.0)
     error->one(FLERR,"Invalid eradius in Atoms section of data file");
 
   etag[nlocal] = atoi(values[3]);
   cs[2*nlocal] = atoi(values[4]);
   cs[2*nlocal+1] = atof(values[5]);
 
   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] = atof(values[0]);
   v[m][1] = atof(values[1]);
   v[m][2] = atof(values[2]);
   ervel[m] = atof(values[3]);
 }
 
 /* ----------------------------------------------------------------------
    unpack hybrid quantities from one line in Velocities section of data file
 ------------------------------------------------------------------------- */
 
 int AtomVecWavepacket::data_vel_hybrid(int m, char **values)
 {
   ervel[m] = atof(values[0]);
   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;
 }
 
 /* ----------------------------------------------------------------------
    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 8265e817e..2e13687ad 100644
--- a/src/USER-AWPMD/atom_vec_wavepacket.h
+++ b/src/USER-AWPMD/atom_vec_wavepacket.h
@@ -1,108 +1,109 @@
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
 
    Copyright (2003) Sandia Corporation.  Under the terms of Contract
    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
    certain rights in this software.  This software is distributed under
    the GNU General Public License.
 
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
 /* ----------------------------------------------------------------------
    Contributing author: Ilya Valuev (JIHT RAS)
 ------------------------------------------------------------------------- */
 
 
 #ifdef ATOM_CLASS
 
 AtomStyle(wavepacket,AtomVecWavepacket)
 
 #else
 
 #ifndef LMP_ATOM_VEC_WAVEPACKET_H
 #define LMP_ATOM_VEC_WAVEPACKET_H
 
 #include "atom_vec.h"
 
 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 *);
   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 coeffcients: cre, cim, size is 2*N
   double *cs;
   ///\en force on wavepacket split coeffcients: 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;
 };
 
 }
 
 #endif
 #endif
diff --git a/src/USER-EFF/atom_vec_electron.cpp b/src/USER-EFF/atom_vec_electron.cpp
index 41837f126..7b658f5c2 100644
--- a/src/USER-EFF/atom_vec_electron.cpp
+++ b/src/USER-EFF/atom_vec_electron.cpp
@@ -1,991 +1,1000 @@
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
 
    Copyright (2003) Sandia Corporation.  Under the terms of Contract
    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
    certain rights in this software.  This software is distributed under
    the GNU General Public License.
 
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
 /* ----------------------------------------------------------------------
    Contributing author: Andres Jaramillo-Botero (Caltech)
 ------------------------------------------------------------------------- */
 
 #include "math.h"
 #include "stdlib.h"
+#include "string.h"
 #include "atom_vec_electron.h"
 #include "atom.h"
 #include "comm.h"
 #include "domain.h"
 #include "modify.h"
 #include "force.h"
 #include "fix.h"
 #include "citeme.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
 
 static const char cite_user_eff_package[] =
   "USER-EFF package:\n\n"
   "@Article{Jaramillo-Botero11,\n"
   " author = {A. Jaramillo-Botero, J. Su, A. Qi, W. A. Goddard III},\n"
   " title = {Large-Scale, Long-Term Nonadiabatic Electron Molecular Dynamics for Describing Material Properties and Phenomena in Extreme Environments},\n"
   " journal = {J.~Comp.~Chem.},\n"
   " year =    2011,\n"
   " volume =  32,\n"
   " pages =   {497--512}\n"
   "}\n\n";
 
 /* ---------------------------------------------------------------------- */
 
 AtomVecElectron::AtomVecElectron(LAMMPS *lmp) : AtomVec(lmp)
 {
   if (lmp->citeme) lmp->citeme->add(cite_user_eff_package);
 
   comm_x_only = comm_f_only = 0;
 
   mass_type = 1;
   molecular = 0;
+  forceclearflag = 1;
 
   size_forward = 4;
   size_reverse = 4;
   size_border = 9;
   size_velocity = 3;
   size_data_atom = 8;
   size_data_vel = 5;
   xcol_data = 6;
 
   atom->ecp_flag = 0;
 
   atom->electron_flag = 1;
   atom->q_flag = atom->spin_flag = atom->eradius_flag =
     atom->ervel_flag = atom->erforce_flag = 1;
 }
 
 /* ----------------------------------------------------------------------
    grow atom-electron arrays
    n = 0 grows arrays by a chunk
    n > 0 allocates arrays to size n
 ------------------------------------------------------------------------- */
 
 void AtomVecElectron::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");
 
   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 AtomVecElectron::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;
 }
 
 /* ----------------------------------------------------------------------
    copy atom I info to atom J
 ------------------------------------------------------------------------- */
 
 void AtomVecElectron::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];
 
   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 AtomVecElectron::force_clear(int n, size_t nbytes)
+{
+  memset(&erforce[n],0,nbytes);
+}
+
+/* ---------------------------------------------------------------------- */
+
 int AtomVecElectron::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;
 }
 
 /* ---------------------------------------------------------------------- */
 
 int AtomVecElectron::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];
     }
   } 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];
       }
     } 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];
         }
       }
     }
   }
   return m;
 }
 
 /* ---------------------------------------------------------------------- */
 
 int AtomVecElectron::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];
   }
   return m;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void AtomVecElectron::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 AtomVecElectron::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++];
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 int AtomVecElectron::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++];
   return m;
 }
 
 /* ---------------------------------------------------------------------- */
 
 int AtomVecElectron::pack_reverse(int n, int first, double *buf)
 {
   int i,m,last;
 
   m = 0;
   last = first + n;
   for (i = first; i < last; i++) {
     buf[m++] = f[i][0];
     buf[m++] = f[i][1];
     buf[m++] = f[i][2];
     buf[m++] = erforce[i];
   }
   return m;
 }
 
 /* ---------------------------------------------------------------------- */
 
 int AtomVecElectron::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];
   return m;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void AtomVecElectron::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++];
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 int AtomVecElectron::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++];
   }
   return m;
 }
 
 /* ---------------------------------------------------------------------- */
 
 int AtomVecElectron::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];
     }
   } 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];
     }
   }
 
   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 AtomVecElectron::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++] = v[j][0];
       buf[m++] = v[j][1];
       buf[m++] = v[j][2];
     }
   } 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++] = v[j][0];
         buf[m++] = v[j][1];
         buf[m++] = v[j][2];
       }
     } 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];
         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];
         }
       }
     }
   }
 
   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 AtomVecElectron::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];
   }
   return m;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void AtomVecElectron::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++];
   }
 
   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 AtomVecElectron::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++];
     v[i][0] = buf[m++];
     v[i][1] = buf[m++];
     v[i][2] = 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 AtomVecElectron::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++];
   }
   return m;
 }
 
 /* ----------------------------------------------------------------------
    pack data for atom I for sending to another proc
    xyz must be 1st 3 values, so comm::exchange() can test on them
 ------------------------------------------------------------------------- */
 
 int AtomVecElectron::pack_exchange(int i, double *buf)
 {
   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];
 
   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 AtomVecElectron::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++];
 
   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;
 }
 
 /* ----------------------------------------------------------------------
    size of restart data for all atoms owned by this proc
    include extra data stored by fixes
 ------------------------------------------------------------------------- */
 
 int AtomVecElectron::size_restart()
 {
   int i;
 
   int nlocal = atom->nlocal;
   int n = 15 * 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 AtomVecElectron::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];
 
   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 AtomVecElectron::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++];
 
   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
 ------------------------------------------------------------------------- */
 
 void AtomVecElectron::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] = 0.0;
   spin[nlocal] = 1;
   eradius[nlocal] = 1.0;
   ervel[nlocal] = 0.0;
 
   atom->nlocal++;
 }
 
 /* ----------------------------------------------------------------------
    unpack one line from Atoms section of data file
    initialize other atom quantities
 ------------------------------------------------------------------------- */
 
 void AtomVecElectron::data_atom(double *coord, imageint imagetmp, char **values)
 {
   int nlocal = atom->nlocal;
 
   if (nlocal == nmax) grow(0);
 
   tag[nlocal] = ATOTAGINT(values[0]);
   type[nlocal] = atoi(values[1]);
   if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes)
     error->one(FLERR,"Invalid atom type in Atoms section of data file");
 
   q[nlocal] = atof(values[2]);
   spin[nlocal] = atoi(values[3]);
   if (spin[nlocal] == 3) atom->ecp_flag = 1;
 
   eradius[nlocal] = atof(values[4]);
 
   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 AtomVecElectron::data_atom_hybrid(int nlocal, char **values)
 {
   q[nlocal] = atof(values[0]);
   spin[nlocal] = atoi(values[1]);
   eradius[nlocal] = atof(values[2]);
   if (eradius[nlocal] < 0.0)
     error->one(FLERR,"Invalid eradius in Atoms section of data file");
 
   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 AtomVecElectron::data_vel(int m, char **values)
 {
   v[m][0] = atof(values[0]);
   v[m][1] = atof(values[1]);
   v[m][2] = atof(values[2]);
   ervel[m] = atof(values[3]);
 }
 
 /* ----------------------------------------------------------------------
    unpack hybrid quantities from one line in Velocities section of data file
 ------------------------------------------------------------------------- */
 
 int AtomVecElectron::data_vel_hybrid(int m, char **values)
 {
   ervel[m] = atof(values[0]);
   return 1;
 }
 
 /* ----------------------------------------------------------------------
    pack atom info for data file including 3 image flags
 ------------------------------------------------------------------------- */
 
 void AtomVecElectron::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] = x[i][0];
     buf[i][6] = x[i][1];
     buf[i][7] = x[i][2];
     buf[i][8] = ubuf((image[i] & IMGMASK) - IMGMAX).d;
     buf[i][9] = ubuf((image[i] >> IMGBITS & IMGMASK) - IMGMAX).d;
     buf[i][10] = ubuf((image[i] >> IMG2BITS) - IMGMAX).d;
   }
 }
 
 /* ----------------------------------------------------------------------
    pack hybrid atom info for data file
 ------------------------------------------------------------------------- */
 
 int AtomVecElectron::pack_data_hybrid(int i, double *buf)
 {
   buf[0] = q[i];
   buf[1] = ubuf(spin[i]).d;
   buf[2] = eradius[i];
   return 3;
 }
 
 /* ----------------------------------------------------------------------
    write atom info to data file including 3 image flags
 ------------------------------------------------------------------------- */
 
 void AtomVecElectron::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 %-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],buf[i][5],buf[i][6],buf[i][7],
             (int) ubuf(buf[i][8]).i,(int) ubuf(buf[i][9]).i,
             (int) ubuf(buf[i][10]).i);
 }
 
 /* ----------------------------------------------------------------------
    write hybrid atom info to data file
 ------------------------------------------------------------------------- */
 
 int AtomVecElectron::write_data_hybrid(FILE *fp, double *buf)
 {
   fprintf(fp," %-1.16e %d %-1.16e",buf[0],(int) ubuf(buf[1]).i,buf[2]);
   return 3;
 }
 
 /* ----------------------------------------------------------------------
    pack velocity info for data file
 ------------------------------------------------------------------------- */
 
 void AtomVecElectron::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 velocity info for data file
 ------------------------------------------------------------------------- */
 
 int AtomVecElectron::pack_vel_hybrid(int i, double *buf)
 {
   buf[0] = ervel[i];
   return 1;
 }
 
 /* ----------------------------------------------------------------------
    write hybrid velocity info to data file
 ------------------------------------------------------------------------- */
 
 void AtomVecElectron::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 AtomVecElectron::write_vel_hybrid(FILE *fp, double *buf)
 {
   fprintf(fp," %-1.16e",buf[0]);
   return 1;
 }
 
 /* ----------------------------------------------------------------------
    return # of bytes of allocated memory
 ------------------------------------------------------------------------- */
 
 bigint AtomVecElectron::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);
 
   return bytes;
 }
diff --git a/src/USER-EFF/atom_vec_electron.h b/src/USER-EFF/atom_vec_electron.h
index f60461359..6a202c47f 100644
--- a/src/USER-EFF/atom_vec_electron.h
+++ b/src/USER-EFF/atom_vec_electron.h
@@ -1,82 +1,83 @@
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
 
    Copyright (2003) Sandia Corporation.  Under the terms of Contract
    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
    certain rights in this software.  This software is distributed under
    the GNU General Public License.
 
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
 #ifdef ATOM_CLASS
 
 AtomStyle(electron,AtomVecElectron)
 
 #else
 
 #ifndef LMP_ATOM_VEC_ELECTRON_H
 #define LMP_ATOM_VEC_ELECTRON_H
 
 #include "atom_vec.h"
 
 namespace LAMMPS_NS {
 
 class AtomVecElectron : public AtomVec {
  public:
   AtomVecElectron(class LAMMPS *);
   ~AtomVecElectron() {}
   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 *);
   bigint memory_usage();
 
  private:
   tagint *tag;
   int *type,*mask;
   imageint *image;
   double **x,**v,**f;
   int *spin;
   double *q,*eradius,*ervel,*erforce;
 };
 
 }
 
 #endif
 #endif
diff --git a/src/USER-SPH/atom_vec_meso.cpp b/src/USER-SPH/atom_vec_meso.cpp
index 3a57776d9..a3e1b8788 100644
--- a/src/USER-SPH/atom_vec_meso.cpp
+++ b/src/USER-SPH/atom_vec_meso.cpp
@@ -1,971 +1,981 @@
 /* ----------------------------------------------------------------------
  LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
  http://lammps.sandia.gov, Sandia National Laboratories
  Steve Plimpton, sjplimp@sandia.gov
 
  Copyright (2003) Sandia Corporation.  Under the terms of Contract
  DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
  certain rights in this software.  This software is distributed under
  the GNU General Public License.
 
  See the README file in the top-level LAMMPS directory.
  ------------------------------------------------------------------------- */
 
+#include "string.h"
 #include "stdlib.h"
 #include "atom_vec_meso.h"
 #include "atom.h"
 #include "comm.h"
 #include "domain.h"
 #include "modify.h"
 #include "fix.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
 AtomVecMeso::AtomVecMeso(LAMMPS *lmp) : AtomVec(lmp)
 {
   molecular = 0;
   mass_type = 1;
+  forceclearflag = 1;
 
   comm_x_only = 0; // we communicate not only x forward but also vest ...
   comm_f_only = 0; // we also communicate de and drho in reverse direction
   size_forward = 8; // 3 + rho + e + vest[3], that means we may only communicate 5 in hybrid
   size_reverse = 5; // 3 + drho + de
   size_border = 12; // 6 + rho + e + vest[3] + cv
   size_velocity = 3;
   size_data_atom = 8;
   size_data_vel = 4;
   xcol_data = 6;
 
   atom->e_flag = 1;
   atom->rho_flag = 1;
   atom->cv_flag = 1;
   atom->vest_flag = 1;
 }
 
 /* ----------------------------------------------------------------------
    grow atom arrays
    n = 0 grows arrays by a chunk
    n > 0 allocates arrays to size n
    ------------------------------------------------------------------------- */
 
 void AtomVecMeso::grow(int n)
 {
   if (n == 0) grow_nmax();
   else nmax = n;
   atom->nmax = nmax;
   if (nmax < 0 || nmax > MAXSMALLINT)
     error->one(FLERR,"Per-processor system is too big");
 
   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");
 
   rho = memory->grow(atom->rho, nmax, "atom:rho");
   drho = memory->grow(atom->drho, nmax*comm->nthreads, "atom:drho");
   e = memory->grow(atom->e, nmax, "atom:e");
   de = memory->grow(atom->de, nmax*comm->nthreads, "atom:de");
   vest = memory->grow(atom->vest, nmax, 3, "atom:vest");
   cv = memory->grow(atom->cv, nmax, "atom:cv");
 
   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 AtomVecMeso::grow_reset() {
   tag = atom->tag;
   type = atom->type;
   mask = atom->mask;
   image = atom->image;
   x = atom->x;
   v = atom->v;
   f = atom->f;
   rho = atom->rho;
   drho = atom->drho;
   e = atom->e;
   de = atom->de;
   vest = atom->vest;
   cv = atom->cv;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void AtomVecMeso::copy(int i, int j, int delflag) {
   //printf("in AtomVecMeso::copy\n");
   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];
 
   rho[j] = rho[i];
   drho[j] = drho[i];
   e[j] = e[i];
   de[j] = de[i];
   cv[j] = cv[i];
   vest[j][0] = vest[i][0];
   vest[j][1] = vest[i][1];
   vest[j][2] = vest[i][2];
 
   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 AtomVecMeso::force_clear(int n, size_t nbytes)
+{
+  memset(&de[n],0,nbytes);
+  memset(&drho[n],0,nbytes);
+}
+
+/* ---------------------------------------------------------------------- */
+
 int AtomVecMeso::pack_comm_hybrid(int n, int *list, double *buf) {
   //printf("in AtomVecMeso::pack_comm_hybrid\n");
   int i, j, m;
 
   m = 0;
   for (i = 0; i < n; i++) {
     j = list[i];
     buf[m++] = rho[j];
     buf[m++] = e[j];
     buf[m++] = vest[j][0];
     buf[m++] = vest[j][1];
     buf[m++] = vest[j][2];
   }
   return m;
 }
 
 /* ---------------------------------------------------------------------- */
 
 int AtomVecMeso::unpack_comm_hybrid(int n, int first, double *buf) {
   //printf("in AtomVecMeso::unpack_comm_hybrid\n");
   int i, m, last;
 
   m = 0;
   last = first + n;
   for (i = first; i < last; i++) {
     rho[i] = buf[m++];
     e[i] = buf[m++];
     vest[i][0] = buf[m++];
     vest[i][1] = buf[m++];
     vest[i][2] = buf[m++];
   }
   return m;
 }
 
 /* ---------------------------------------------------------------------- */
 
 int AtomVecMeso::pack_border_hybrid(int n, int *list, double *buf) {
   //printf("in AtomVecMeso::pack_border_hybrid\n");
   int i, j, m;
 
   m = 0;
   for (i = 0; i < n; i++) {
     j = list[i];
     buf[m++] = rho[j];
     buf[m++] = e[j];
     buf[m++] = vest[j][0];
     buf[m++] = vest[j][1];
     buf[m++] = vest[j][2];
   }
   return m;
 }
 
 /* ---------------------------------------------------------------------- */
 
 int AtomVecMeso::unpack_border_hybrid(int n, int first, double *buf) {
   //printf("in AtomVecMeso::unpack_border_hybrid\n");
   int i, m, last;
 
   m = 0;
   last = first + n;
   for (i = first; i < last; i++) {
     rho[i] = buf[m++];
     e[i] = buf[m++];
     vest[i][0] = buf[m++];
     vest[i][1] = buf[m++];
     vest[i][2] = buf[m++];
   }
   return m;
 }
 
 /* ---------------------------------------------------------------------- */
 
 int AtomVecMeso::pack_reverse_hybrid(int n, int first, double *buf) {
   //printf("in AtomVecMeso::pack_reverse_hybrid\n");
   int i, m, last;
 
   m = 0;
   last = first + n;
   for (i = first; i < last; i++) {
     buf[m++] = drho[i];
     buf[m++] = de[i];
   }
   return m;
 }
 
 /* ---------------------------------------------------------------------- */
 
 int AtomVecMeso::unpack_reverse_hybrid(int n, int *list, double *buf) {
   //printf("in AtomVecMeso::unpack_reverse_hybrid\n");
   int i, j, m;
 
   m = 0;
   for (i = 0; i < n; i++) {
     j = list[i];
     drho[j] += buf[m++];
     de[j] += buf[m++];
   }
   return m;
 }
 
 /* ---------------------------------------------------------------------- */
 
 int AtomVecMeso::pack_comm(int n, int *list, double *buf, int pbc_flag,
                            int *pbc) {
   //printf("in AtomVecMeso::pack_comm\n");
   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++] = rho[j];
       buf[m++] = e[j];
       buf[m++] = vest[j][0];
       buf[m++] = vest[j][1];
       buf[m++] = vest[j][2];
     }
   } 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++] = rho[j];
       buf[m++] = e[j];
       buf[m++] = vest[j][0];
       buf[m++] = vest[j][1];
       buf[m++] = vest[j][2];
     }
   }
   return m;
 }
 
 /* ---------------------------------------------------------------------- */
 
 int AtomVecMeso::pack_comm_vel(int n, int *list, double *buf, int pbc_flag,
                                int *pbc) {
   //printf("in AtomVecMeso::pack_comm_vel\n");
   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++] = v[j][0];
       buf[m++] = v[j][1];
       buf[m++] = v[j][2];
       buf[m++] = rho[j];
       buf[m++] = e[j];
       buf[m++] = vest[j][0];
       buf[m++] = vest[j][1];
       buf[m++] = vest[j][2];
     }
   } 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++] = v[j][0];
       buf[m++] = v[j][1];
       buf[m++] = v[j][2];
       buf[m++] = rho[j];
       buf[m++] = e[j];
       buf[m++] = vest[j][0];
       buf[m++] = vest[j][1];
       buf[m++] = vest[j][2];
     }
   }
   return m;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void AtomVecMeso::unpack_comm(int n, int first, double *buf) {
   //printf("in AtomVecMeso::unpack_comm\n");
   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++];
     rho[i] = buf[m++];
     e[i] = buf[m++];
     vest[i][0] = buf[m++];
     vest[i][1] = buf[m++];
     vest[i][2] = buf[m++];
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void AtomVecMeso::unpack_comm_vel(int n, int first, double *buf) {
   //printf("in AtomVecMeso::unpack_comm_vel\n");
   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++];
     v[i][0] = buf[m++];
     v[i][1] = buf[m++];
     v[i][2] = buf[m++];
     rho[i] = buf[m++];
     e[i] = buf[m++];
     vest[i][0] = buf[m++];
     vest[i][1] = buf[m++];
     vest[i][2] = buf[m++];
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 int AtomVecMeso::pack_reverse(int n, int first, double *buf) {
   //printf("in AtomVecMeso::pack_reverse\n");
   int i, m, last;
 
   m = 0;
   last = first + n;
   for (i = first; i < last; i++) {
     buf[m++] = f[i][0];
     buf[m++] = f[i][1];
     buf[m++] = f[i][2];
     buf[m++] = drho[i];
     buf[m++] = de[i];
   }
   return m;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void AtomVecMeso::unpack_reverse(int n, int *list, double *buf) {
   //printf("in AtomVecMeso::unpack_reverse\n");
   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++];
     drho[j] += buf[m++];
     de[j] += buf[m++];
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 int AtomVecMeso::pack_border(int n, int *list, double *buf, int pbc_flag,
                              int *pbc) {
   //printf("in AtomVecMeso::pack_border\n");
   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++] = rho[j];
       buf[m++] = e[j];
       buf[m++] = cv[j];
       buf[m++] = vest[j][0];
       buf[m++] = vest[j][1];
       buf[m++] = vest[j][2];
     }
   } 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++] = rho[j];
       buf[m++] = e[j];
       buf[m++] = cv[j];
       buf[m++] = vest[j][0];
       buf[m++] = vest[j][1];
       buf[m++] = vest[j][2];
     }
   }
 
   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 AtomVecMeso::pack_border_vel(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++] = v[j][0];
       buf[m++] = v[j][1];
       buf[m++] = v[j][2];
       buf[m++] = rho[j];
       buf[m++] = e[j];
       buf[m++] = cv[j];
       buf[m++] = vest[j][0];
       buf[m++] = vest[j][1];
       buf[m++] = vest[j][2];
     }
   } 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 (!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++] = ubuf(tag[j]).d;
         buf[m++] = ubuf(type[j]).d;
         buf[m++] = ubuf(mask[j]).d;
         buf[m++] = v[j][0];
         buf[m++] = v[j][1];
         buf[m++] = v[j][2];
         buf[m++] = rho[j];
         buf[m++] = e[j];
         buf[m++] = cv[j];
         buf[m++] = vest[j][0];
         buf[m++] = vest[j][1];
         buf[m++] = vest[j][2];
       }
     } else {
       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;
         if (mask[i] & deform_groupbit) {
           buf[m++] = v[j][0];
           buf[m++] = v[j][1];
           buf[m++] = v[j][2];
         } else {
           buf[m++] = v[j][0];
           buf[m++] = v[j][1];
           buf[m++] = v[j][2];
         }
         buf[m++] = rho[j];
         buf[m++] = e[j];
         buf[m++] = cv[j];
         buf[m++] = vest[j][0];
         buf[m++] = vest[j][1];
         buf[m++] = vest[j][2];
       }
     }
   }
 
   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;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void AtomVecMeso::unpack_border(int n, int first, double *buf) {
   //printf("in AtomVecMeso::unpack_border\n");
   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;
     rho[i] = buf[m++];
     e[i] = buf[m++];
     cv[i] = buf[m++];
     vest[i][0] = buf[m++];
     vest[i][1] = buf[m++];
     vest[i][2] = 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]);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void AtomVecMeso::unpack_border_vel(int n, int first, double *buf) {
   //printf("in AtomVecMeso::unpack_border_vel\n");
   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;
     v[i][0] = buf[m++];
     v[i][1] = buf[m++];
     v[i][2] = buf[m++];
     rho[i] = buf[m++];
     e[i] = buf[m++];
     cv[i] = buf[m++];
     vest[i][0] = buf[m++];
     vest[i][1] = buf[m++];
     vest[i][2] = 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]);
 }
 
 /* ----------------------------------------------------------------------
    pack data for atom I for sending to another proc
    xyz must be 1st 3 values, so comm::exchange() can test on them
    ------------------------------------------------------------------------- */
 
 int AtomVecMeso::pack_exchange(int i, double *buf) {
   //printf("in AtomVecMeso::pack_exchange\n");
   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++] = rho[i];
   buf[m++] = e[i];
   buf[m++] = cv[i];
   buf[m++] = vest[i][0];
   buf[m++] = vest[i][1];
   buf[m++] = vest[i][2];
 
   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 AtomVecMeso::unpack_exchange(double *buf) {
   //printf("in AtomVecMeso::unpack_exchange\n");
   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;
   rho[nlocal] = buf[m++];
   e[nlocal] = buf[m++];
   cv[nlocal] = buf[m++];
   vest[nlocal][0] = buf[m++];
   vest[nlocal][1] = buf[m++];
   vest[nlocal][2] = 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;
 }
 
 /* ----------------------------------------------------------------------
    size of restart data for all atoms owned by this proc
    include extra data stored by fixes
    ------------------------------------------------------------------------- */
 
 int AtomVecMeso::size_restart() {
   int i;
 
   int nlocal = atom->nlocal;
   int n = 17 * nlocal; // 11 + rho + e + cv + vest[3]
 
   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 AtomVecMeso::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++] = rho[i];
   buf[m++] = e[i];
   buf[m++] = cv[i];
   buf[m++] = vest[i][0];
   buf[m++] = vest[i][1];
   buf[m++] = vest[i][2];
 
   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 AtomVecMeso::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++];
   rho[nlocal] = buf[m++];
   e[nlocal] = buf[m++];
   cv[nlocal] = buf[m++];
   vest[nlocal][0] = buf[m++];
   vest[nlocal][1] = buf[m++];
   vest[nlocal][2] = 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
    ------------------------------------------------------------------------- */
 
 void AtomVecMeso::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;
   rho[nlocal] = 0.0;
   e[nlocal] = 0.0;
   cv[nlocal] = 1.0;
   vest[nlocal][0] = 0.0;
   vest[nlocal][1] = 0.0;
   vest[nlocal][2] = 0.0;
   de[nlocal] = 0.0;
   drho[nlocal] = 0.0;
 
   atom->nlocal++;
 }
 
 /* ----------------------------------------------------------------------
    unpack one line from Atoms section of data file
    initialize other atom quantities
    ------------------------------------------------------------------------- */
 
 void AtomVecMeso::data_atom(double *coord, imageint imagetmp, char **values) {
   int nlocal = atom->nlocal;
   if (nlocal == nmax) grow(0);
 
   tag[nlocal] = ATOTAGINT(values[0]);
   type[nlocal] = atoi(values[1]);
   if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes)
     error->one(FLERR,"Invalid atom type in Atoms section of data file");
 
   rho[nlocal] = atof(values[2]);
   e[nlocal] = atof(values[3]);
   cv[nlocal] = atof(values[4]);
 
   x[nlocal][0] = coord[0];
   x[nlocal][1] = coord[1];
   x[nlocal][2] = coord[2];
 
   //printf("rho=%f, e=%f, cv=%f, x=%f\n", rho[nlocal], e[nlocal], cv[nlocal], x[nlocal][0]);
 
   image[nlocal] = imagetmp;
 
   mask[nlocal] = 1;
   v[nlocal][0] = 0.0;
   v[nlocal][1] = 0.0;
   v[nlocal][2] = 0.0;
 
   vest[nlocal][0] = 0.0;
   vest[nlocal][1] = 0.0;
   vest[nlocal][2] = 0.0;
 
   de[nlocal] = 0.0;
   drho[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 AtomVecMeso::data_atom_hybrid(int nlocal, char **values) {
 
   rho[nlocal] = atof(values[0]);
   e[nlocal] = atof(values[1]);
   cv[nlocal] = atof(values[2]);
 
   return 3;
 }
 
 /* ----------------------------------------------------------------------
    pack atom info for data file including 3 image flags
 ------------------------------------------------------------------------- */
 
 void AtomVecMeso::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] = rho[i];
     buf[i][3] = e[i];
     buf[i][4] = cv[i];
     buf[i][5] = x[i][0];
     buf[i][6] = x[i][1];
     buf[i][7] = x[i][2];
     buf[i][8] = ubuf((image[i] & IMGMASK) - IMGMAX).d;
     buf[i][9] = ubuf((image[i] >> IMGBITS & IMGMASK) - IMGMAX).d;
     buf[i][10] = ubuf((image[i] >> IMG2BITS) - IMGMAX).d;
   }
 }
 
 /* ----------------------------------------------------------------------
    pack hybrid atom info for data file
 ------------------------------------------------------------------------- */
 
 int AtomVecMeso::pack_data_hybrid(int i, double *buf)
 {
   buf[0] = rho[i];
   buf[1] = e[i];
   buf[2] = cv[i];
   return 3;
 }
 
 /* ----------------------------------------------------------------------
    write atom info to data file including 3 image flags
 ------------------------------------------------------------------------- */
 
 void AtomVecMeso::write_data(FILE *fp, int n, double **buf)
 {
   for (int i = 0; i < n; i++)
     fprintf(fp,TAGINT_FORMAT 
             " %d %-1.16e %-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],buf[i][3],buf[i][4],
             buf[i][5],buf[i][6],buf[i][7],
             (int) ubuf(buf[i][8]).i,(int) ubuf(buf[i][9]).i,
             (int) ubuf(buf[i][10]).i);
 }
 
 /* ----------------------------------------------------------------------
    write hybrid atom info to data file
 ------------------------------------------------------------------------- */
 
 int AtomVecMeso::write_data_hybrid(FILE *fp, double *buf)
 {
   fprintf(fp," %-1.16e %-1.16e %-1.16e",buf[0],buf[1],buf[2]);
   return 3;
 }
 
 /* ----------------------------------------------------------------------
    return # of bytes of allocated memory
    ------------------------------------------------------------------------- */
 
 bigint AtomVecMeso::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("rho"))
     bytes += memory->usage(rho, nmax);
   if (atom->memcheck("drho"))
     bytes += memory->usage(drho, nmax*comm->nthreads);
   if (atom->memcheck("e"))
     bytes += memory->usage(e, nmax);
   if (atom->memcheck("de"))
     bytes += memory->usage(de, nmax*comm->nthreads);
   if (atom->memcheck("cv"))
     bytes += memory->usage(cv, nmax);
   if (atom->memcheck("vest"))
     bytes += memory->usage(vest, nmax);
 
   return bytes;
 }
diff --git a/src/USER-SPH/atom_vec_meso.h b/src/USER-SPH/atom_vec_meso.h
index aff8b7672..6f50b44c2 100644
--- a/src/USER-SPH/atom_vec_meso.h
+++ b/src/USER-SPH/atom_vec_meso.h
@@ -1,76 +1,77 @@
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
 
    Copyright (2003) Sandia Corporation.  Under the terms of Contract
    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
    certain rights in this software.  This software is distributed under
    the GNU General Public License.
 
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
 #ifdef ATOM_CLASS
 
 AtomStyle(meso,AtomVecMeso)
 
 #else
 
 #ifndef LMP_ATOM_VEC_MESO_H
 #define LMP_ATOM_VEC_MESO_H
 
 #include "atom_vec.h"
 
 namespace LAMMPS_NS {
 
 class AtomVecMeso : public AtomVec {
  public:
   AtomVecMeso(class LAMMPS *);
   ~AtomVecMeso() {}
   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 *);
   void unpack_comm(int, int, double *);
   void unpack_comm_vel(int, int, double *);
   int pack_reverse(int, int, double *);
   void unpack_reverse(int, int *, double *);
   int pack_comm_hybrid(int, int *, double *);
   int unpack_comm_hybrid(int, int, double *);
   int pack_border_hybrid(int, int *, double *);
   int unpack_border_hybrid(int, int, double *);
   int pack_reverse_hybrid(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 *);
   void unpack_border(int, int, double *);
   void unpack_border_vel(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 pack_data(double **);
   int pack_data_hybrid(int, double *);
   void write_data(FILE *, int, double **);
   int write_data_hybrid(FILE *, double *);
   bigint memory_usage();
 
  private:
   tagint *tag;
   int *type,*mask;
   imageint *image;
   double **x,**v,**f;
   double *rho, *drho, *e, *de, *cv;
   double **vest; // estimated velocity during force computation
 };
 
 }
 
 #endif
 #endif
diff --git a/src/atom.cpp b/src/atom.cpp
index cb394e44e..1d47b1114 100644
--- a/src/atom.cpp
+++ b/src/atom.cpp
@@ -1,1935 +1,1985 @@
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
 
    Copyright (2003) Sandia Corporation.  Under the terms of Contract
    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
    certain rights in this software.  This software is distributed under
    the GNU General Public License.
 
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
 #include "mpi.h"
 #include "math.h"
 #include "stdio.h"
 #include "stdlib.h"
 #include "string.h"
 #include "limits.h"
 #include "atom.h"
 #include "style_atom.h"
 #include "atom_vec.h"
 #include "atom_vec_ellipsoid.h"
 #include "comm.h"
 #include "neighbor.h"
 #include "force.h"
 #include "modify.h"
 #include "fix.h"
 #include "output.h"
 #include "thermo.h"
 #include "update.h"
 #include "domain.h"
 #include "group.h"
 #include "molecule.h"
 #include "accelerator_cuda.h"
 #include "atom_masks.h"
 #include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
 using namespace MathConst;
 
 #define DELTA 1
 #define DELTA_MEMSTR 1024
 #define EPSILON 1.0e-6
 #define CUDA_CHUNK 3000
 #define MAXBODY 20       // max # of lines in one body, also in ReadData class
 
 /* ---------------------------------------------------------------------- */
 
 Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
 {
   natoms = 0;
   nlocal = nghost = nmax = 0;
   ntypes = 0;
   nbondtypes = nangletypes = ndihedraltypes = nimpropertypes = 0;
   nbonds = nangles = ndihedrals = nimpropers = 0;
 
   firstgroupname = NULL;
   sortfreq = 1000;
   nextsort = 0;
   userbinsize = 0.0;
   maxbin = maxnext = 0;
   binhead = NULL;
   next = permute = NULL;
 
   // initialize atom arrays
   // customize by adding new array
 
   tag = NULL;
   type = mask = NULL;
   image = NULL;
   x = v = f = NULL;
 
   molecule = NULL;
   molindex = molatom = NULL;
   q = NULL;
   mu = NULL;
   omega = angmom = torque = NULL;
   radius = rmass = NULL;
+  ellipsoid = line = tri = body = NULL;
+
   vfrac = s0 = NULL;
   x0 = NULL;
-  ellipsoid = line = tri = body = NULL;
+
   spin = NULL;
   eradius = ervel = erforce = NULL;
   cs = csforce = vforce = ervelforce = NULL;
   etag = NULL;
-  rho = drho = NULL;
-  e = de = NULL;
-  cv = NULL;
+
+  rho = drho = e = de = cv = NULL;
   vest = NULL;
 
   bond_per_atom =  extra_bond_per_atom = 0;
   num_bond = NULL;
   bond_type = NULL;
   bond_atom = NULL;
 
   angle_per_atom = extra_angle_per_atom = 0;
   num_angle = NULL;
   angle_type = NULL;
   angle_atom1 = angle_atom2 = angle_atom3 = NULL;
 
   dihedral_per_atom = extra_dihedral_per_atom = 0;
   num_dihedral = NULL;
   dihedral_type = NULL;
   dihedral_atom1 = dihedral_atom2 = dihedral_atom3 = dihedral_atom4 = NULL;
 
   improper_per_atom = extra_improper_per_atom = 0;
   num_improper = NULL;
   improper_type = NULL;
   improper_atom1 = improper_atom2 = improper_atom3 = improper_atom4 = NULL;
 
   maxspecial = 1;
   nspecial = NULL;
   special = NULL;
 
   // user-defined molecules
 
   nmolecule = 0;
   molecules = NULL;
 
   // custom atom arrays
 
   nivector = ndvector = 0;
   ivector = NULL;
   dvector = NULL;
   iname = dname = NULL;
 
   // initialize atom style and array existence flags
   // customize by adding new flag
 
-  sphere_flag = ellipsoid_flag = line_flag = tri_flag = body_flag = 0;
-  peri_flag = electron_flag = 0;
+  sphere_flag = peri_flag = electron_flag = 0;
   wavepacket_flag = sph_flag = 0;
 
-  molecule_flag = q_flag = mu_flag = 0;
-  rmass_flag = radius_flag = omega_flag = torque_flag = angmom_flag = 0;
-  vfrac_flag = spin_flag = eradius_flag = ervel_flag = erforce_flag = 0;
-  cs_flag = csforce_flag = vforce_flag = ervelforce_flag= etag_flag = 0;
+  molecule_flag = 0;
+  q_flag = mu_flag = 0;
+  omega_flag = torque_flag = angmom_flag = 0;
+  radius_flag = rmass_flag = 0;
+  ellipsoid_flag = line_flag = tri_flag = body_flag = 0;
+
+  vfrac_flag = 0;
+  spin_flag = eradius_flag = ervel_flag = erforce_flag = ervelforce_flag = 0;
+  cs_flag = csforce_flag = vforce_flag = etag_flag = 0;
+
   rho_flag = e_flag = cv_flag = vest_flag = 0;
 
   // Peridynamic scale factor
 
   pdscale = 1.0;
 
   // ntype-length arrays
 
   mass = NULL;
   mass_setflag = NULL;
 
   // callback lists & extra restart info
 
   nextra_grow = nextra_restart = nextra_border = 0;
   extra_grow = extra_restart = extra_border = NULL;
   nextra_grow_max = nextra_restart_max = nextra_border_max = 0;
   nextra_store = 0;
   extra = NULL;
 
   // default atom ID and mapping values
 
   tag_enable = 1;
   map_style = map_user = 0;
   map_tag_max = -1;
   map_maxarray = map_nhash = -1;
 
   max_same = 0;
   sametag = NULL;
   map_array = NULL;
   map_bucket = NULL;
   map_hash = NULL;
 
   atom_style = NULL;
   avec = NULL;
 
   datamask = ALL_MASK;
   datamask_ext = ALL_MASK;
 }
 
 /* ---------------------------------------------------------------------- */
 
 Atom::~Atom()
 {
   delete [] atom_style;
   delete avec;
 
   delete [] firstgroupname;
   memory->destroy(binhead);
   memory->destroy(next);
   memory->destroy(permute);
 
   // delete atom arrays
   // customize by adding new array
 
   memory->destroy(tag);
   memory->destroy(type);
   memory->destroy(mask);
   memory->destroy(image);
   memory->destroy(x);
   memory->destroy(v);
   memory->destroy(f);
 
+  memory->destroy(molecule);
+  memory->destroy(molindex);
+  memory->destroy(molatom);
+
   memory->destroy(q);
   memory->destroy(mu);
   memory->destroy(omega);
   memory->destroy(angmom);
   memory->destroy(torque);
   memory->destroy(radius);
   memory->destroy(rmass);
-  memory->destroy(vfrac);
-  memory->destroy(s0);
-  memory->destroy(x0);
   memory->destroy(ellipsoid);
   memory->destroy(line);
   memory->destroy(tri);
   memory->destroy(body);
+
+  memory->destroy(vfrac);
+  memory->destroy(s0);
+  memory->destroy(x0);
+
   memory->destroy(spin);
   memory->destroy(eradius);
   memory->destroy(ervel);
   memory->destroy(erforce);
-
-  memory->destroy(molecule);
-  memory->destroy(molindex);
-  memory->destroy(molatom);
+  memory->destroy(ervelforce);
+  memory->destroy(cs);
+  memory->destroy(csforce);
+  memory->destroy(vforce);
+  memory->destroy(etag);
+
+  memory->destroy(rho);
+  memory->destroy(drho);
+  memory->destroy(e);
+  memory->destroy(de);
+  memory->destroy(cv);
+  memory->destroy(vest);
 
   memory->destroy(nspecial);
   memory->destroy(special);
 
   memory->destroy(num_bond);
   memory->destroy(bond_type);
   memory->destroy(bond_atom);
 
   memory->destroy(num_angle);
   memory->destroy(angle_type);
   memory->destroy(angle_atom1);
   memory->destroy(angle_atom2);
   memory->destroy(angle_atom3);
 
   memory->destroy(num_dihedral);
   memory->destroy(dihedral_type);
   memory->destroy(dihedral_atom1);
   memory->destroy(dihedral_atom2);
   memory->destroy(dihedral_atom3);
   memory->destroy(dihedral_atom4);
 
   memory->destroy(num_improper);
   memory->destroy(improper_type);
   memory->destroy(improper_atom1);
   memory->destroy(improper_atom2);
   memory->destroy(improper_atom3);
   memory->destroy(improper_atom4);
 
-  // delete user-defined molecules
-
-  for (int i = 0; i < nmolecule; i++) delete molecules[i];
-  memory->sfree(molecules);
-
   // delete custom atom arrays
 
   for (int i = 0; i < nivector; i++) {
     delete [] iname[i];
     memory->destroy(ivector[i]);
   }
   for (int i = 0; i < ndvector; i++) {
     delete [] dname[i];
     memory->destroy(dvector[i]);
   }
 
   memory->sfree(iname);
   memory->sfree(dname);
   memory->sfree(ivector);
   memory->sfree(dvector);
 
+  // delete user-defined molecules
+
+  for (int i = 0; i < nmolecule; i++) delete molecules[i];
+  memory->sfree(molecules);
+
   // delete per-type arrays
 
   delete [] mass;
   delete [] mass_setflag;
 
   // delete extra arrays
 
   memory->destroy(extra_grow);
   memory->destroy(extra_restart);
   memory->destroy(extra_border);
   memory->destroy(extra);
 
   // delete mapping data structures
 
   map_delete();
 }
 
 /* ----------------------------------------------------------------------
    copy modify settings from old Atom class to current Atom class
 ------------------------------------------------------------------------- */
 
 void Atom::settings(Atom *old)
 {
   tag_enable = old->tag_enable;
   map_user = old->map_user;
   map_style = old->map_style;
   sortfreq = old->sortfreq;
   userbinsize = old->userbinsize;
   if (old->firstgroupname) {
     int n = strlen(old->firstgroupname) + 1;
     firstgroupname = new char[n];
     strcpy(firstgroupname,old->firstgroupname);
   }
 }
 
 /* ----------------------------------------------------------------------
    create an AtomVec style
    called from lammps.cpp, input script, restart file, replicate
 ------------------------------------------------------------------------- */
 
 void Atom::create_avec(const char *style, int narg, char **arg, char *suffix)
 {
   delete [] atom_style;
   if (avec) delete avec;
 
   // unset atom style and array existence flags
   // may have been set by old avec
   // customize by adding new flag
 
-  sphere_flag = ellipsoid_flag = line_flag = tri_flag = 0;
-  peri_flag = electron_flag = 0;
+  sphere_flag = peri_flag = electron_flag = 0;
+  wavepacket_flag = sph_flag = 0;
 
-  molecule_flag = q_flag = mu_flag = 0;
-  rmass_flag = radius_flag = omega_flag = torque_flag = angmom_flag = 0;
-  vfrac_flag = spin_flag = eradius_flag = ervel_flag = erforce_flag = 0;
+  molecule_flag = 0;
+  q_flag = mu_flag = 0;
+  omega_flag = torque_flag = angmom_flag = 0;
+  radius_flag = rmass_flag = 0;
+  ellipsoid_flag = line_flag = tri_flag = body_flag = 0;
+
+  vfrac_flag = 0;
+  spin_flag = eradius_flag = ervel_flag = erforce_flag = ervelforce_flag = 0;
+  cs_flag = csforce_flag = vforce_flag = etag_flag = 0;
+
+  rho_flag = e_flag = cv_flag = vest_flag = 0;
 
   // create instance of AtomVec
   // use grow() to initialize atom-based arrays to length 1
   //   so that x[0][0] can always be referenced even if proc has no atoms
 
   int sflag;
   avec = new_avec(style,suffix,sflag);
   avec->store_args(narg,arg);
   avec->process_args(narg,arg);
   avec->grow(1);
 
   if (sflag) {
     char estyle[256];
     sprintf(estyle,"%s/%s",style,suffix);
     int n = strlen(estyle) + 1;
     atom_style = new char[n];
     strcpy(atom_style,estyle);
   } else {
     int n = strlen(style) + 1;
     atom_style = new char[n];
     strcpy(atom_style,style);
   }
 
   // if molecular system:
   // atom IDs must be defined
   // force atom map to be created
   // map style may be reset by map_init() and its call to map_style_set()
 
   molecular = avec->molecular;
   if (molecular && tag_enable == 0)
     error->all(FLERR,"Atom IDs must be used for molecular systems");
   if (molecular) map_style = 1;
 }
 
 /* ----------------------------------------------------------------------
    generate an AtomVec class, first with suffix appended
 ------------------------------------------------------------------------- */
 
 AtomVec *Atom::new_avec(const char *style, char *suffix, int &sflag)
 {
   if (suffix && lmp->suffix_enable) {
     sflag = 1;
     char estyle[256];
     sprintf(estyle,"%s/%s",style,suffix);
 
     if (0) return NULL;
 
 #define ATOM_CLASS
 #define AtomStyle(key,Class) \
     else if (strcmp(estyle,#key) == 0) return new Class(lmp);
 #include "style_atom.h"
 #undef AtomStyle
 #undef ATOM_CLASS
 
   }
 
   sflag = 0;
 
   if (0) return NULL;
 
 #define ATOM_CLASS
 #define AtomStyle(key,Class) \
   else if (strcmp(style,#key) == 0) return new Class(lmp);
 #include "style_atom.h"
 #undef ATOM_CLASS
 
   else error->all(FLERR,"Invalid atom style");
 
   return NULL;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Atom::init()
 {
   // delete extra array since it doesn't persist past first run
 
   if (nextra_store) {
     memory->destroy(extra);
     extra = NULL;
     nextra_store = 0;
   }
 
   // check arrays that are atom type in length
 
   check_mass();
 
   // setup of firstgroup
 
   if (firstgroupname) {
     firstgroup = group->find(firstgroupname);
     if (firstgroup < 0)
       error->all(FLERR,"Could not find atom_modify first group ID");
   } else firstgroup = -1;
 
   // init AtomVec
 
   avec->init();
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Atom::setup()
 {
   // setup bins for sorting
   // cannot do this in init() because uses neighbor cutoff
 
   if (sortfreq > 0) setup_sort_bins();
 }
 
 /* ----------------------------------------------------------------------
    return ptr to AtomVec class if matches style or to matching hybrid sub-class
    return NULL if no match
 ------------------------------------------------------------------------- */
 
 AtomVec *Atom::style_match(const char *style)
 {
   if (strcmp(atom_style,style) == 0) return avec;
   else if (strcmp(atom_style,"hybrid") == 0) {
     AtomVecHybrid *avec_hybrid = (AtomVecHybrid *) avec;
     for (int i = 0; i < avec_hybrid->nstyles; i++)
       if (strcmp(avec_hybrid->keywords[i],style) == 0)
         return avec_hybrid->styles[i];
   }
   return NULL;
 }
 
 /* ----------------------------------------------------------------------
    modify parameters of the atom style
    some options can only be invoked before simulation box is defined
    first and sort options cannot be used together
 ------------------------------------------------------------------------- */
 
 void Atom::modify_params(int narg, char **arg)
 {
   if (narg == 0) error->all(FLERR,"Illegal atom_modify command");
 
   int iarg = 0;
   while (iarg < narg) {
     if (strcmp(arg[iarg],"id") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal atom_modify command");
       if (domain->box_exist)
         error->all(FLERR,
                    "Atom_modify id command after simulation box is defined");
       if (strcmp(arg[iarg+1],"yes") == 0) tag_enable = 1;
       else if (strcmp(arg[iarg+1],"no") == 0) tag_enable = 2;
       else error->all(FLERR,"Illegal atom_modify command");
       iarg += 2;
     } if (strcmp(arg[iarg],"map") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal atom_modify command");
       if (domain->box_exist)
         error->all(FLERR,
                    "Atom_modify map command after simulation box is defined");
       if (strcmp(arg[iarg+1],"array") == 0) map_user = 1;
       else if (strcmp(arg[iarg+1],"hash") == 0) map_user = 2;
       else error->all(FLERR,"Illegal atom_modify command");
       map_style = map_user;
       iarg += 2;
     } else if (strcmp(arg[iarg],"first") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal atom_modify command");
       if (strcmp(arg[iarg+1],"all") == 0) {
         delete [] firstgroupname;
         firstgroupname = NULL;
       } else {
         int n = strlen(arg[iarg+1]) + 1;
         firstgroupname = new char[n];
         strcpy(firstgroupname,arg[iarg+1]);
         sortfreq = 0;
       }
       iarg += 2;
     } else if (strcmp(arg[iarg],"sort") == 0) {
       if (iarg+3 > narg) error->all(FLERR,"Illegal atom_modify command");
       sortfreq = force->inumeric(FLERR,arg[iarg+1]);
       userbinsize = force->numeric(FLERR,arg[iarg+2]);
       if (sortfreq < 0 || userbinsize < 0.0)
         error->all(FLERR,"Illegal atom_modify command");
       if (sortfreq >= 0 && firstgroupname)
         error->all(FLERR,"Atom_modify sort and first options "
                    "cannot be used together");
       iarg += 3;
     } else error->all(FLERR,"Illegal atom_modify command");
   }
 }
 
 /* ----------------------------------------------------------------------
    check that atom IDs are valid
    error if any atom ID < 0 or atom ID = MAXTAGINT
    if any atom ID > 0, error if any atom ID == 0
    if all atom IDs = 0, tag_enable must be 0
    OK if atom IDs > natoms
    NOTE: not checking that atom IDs are unique
 ------------------------------------------------------------------------- */
 
 void Atom::tag_check()
 {
   int nlocal = atom->nlocal;
   tagint *tag = atom->tag;
 
   tagint min = MAXTAGINT;
   tagint max = 0;
 
   for (int i = 0; i < nlocal; i++) {
     min = MIN(min,tag[i]);
     max = MAX(max,tag[i]);
   }
 
   tagint minall,maxall;
   MPI_Allreduce(&min,&minall,1,MPI_LMP_TAGINT,MPI_MIN,world);
   MPI_Allreduce(&max,&maxall,1,MPI_LMP_TAGINT,MPI_MAX,world);
 
   if (minall < 0) error->all(FLERR,"Atom ID is negative");
   if (maxall >= MAXTAGINT) error->all(FLERR,"Atom ID is too big");
   if (maxall > 0 && minall == 0) error->all(FLERR,"Atom ID is zero");
   if (maxall == 0 && tag_enable && natoms) 
     error->all(FLERR,"Not all atom IDs are 0");
 }
 
 /* ----------------------------------------------------------------------
    add unique tags to any atoms with tag = 0
    new tags are grouped by proc and start after max current tag
    called after creating new atoms
    error if new tags will exceed MAXTAGINT
 ------------------------------------------------------------------------- */
 
 void Atom::tag_extend()
 {
   // maxtag_all = max tag for all atoms
 
   tagint maxtag = 0;
   for (int i = 0; i < nlocal; i++) maxtag = MAX(maxtag,tag[i]);
   tagint maxtag_all;
   MPI_Allreduce(&maxtag,&maxtag_all,1,MPI_LMP_TAGINT,MPI_MAX,world);
 
   // DEBUG: useful for generating 64-bit IDs even for small systems
   // use only when LAMMPS is compiled with BIGBIG
 
   //maxtag_all += 1000000000000;
 
   // notag = # of atoms I own with no tag (tag = 0)
   // notag_sum = # of total atoms on procs <= me with no tag
 
   bigint notag = 0;
   for (int i = 0; i < nlocal; i++) if (tag[i] == 0) notag++;
 
   bigint notag_total;
   MPI_Allreduce(&notag,&notag_total,1,MPI_LMP_BIGINT,MPI_SUM,world);
   if (notag_total >= MAXTAGINT)
     error->all(FLERR,"New atom IDs exceed maximum allowed ID");
 
   bigint notag_sum;
   MPI_Scan(&notag,&notag_sum,1,MPI_LMP_BIGINT,MPI_SUM,world);
 
   // itag = 1st new tag that my untagged atoms should use
 
   tagint itag = maxtag_all + notag_sum - notag + 1;
   for (int i = 0; i < nlocal; i++) if (tag[i] == 0) tag[i] = itag++;
 }
 
 /* ----------------------------------------------------------------------
    check that atom IDs span range from 1 to Natoms inclusive
    return 0 if mintag != 1 or maxtag != Natoms
    return 1 if OK
    doesn't actually check if all tag values are used
 ------------------------------------------------------------------------- */
 
 int Atom::tag_consecutive()
 {
   tagint idmin = MAXTAGINT;
   tagint idmax = 0;
 
   for (int i = 0; i < nlocal; i++) {
     idmin = MIN(idmin,tag[i]);
     idmax = MAX(idmax,tag[i]);
   }
   tagint idminall,idmaxall;
   MPI_Allreduce(&idmin,&idminall,1,MPI_LMP_TAGINT,MPI_MIN,world);
   MPI_Allreduce(&idmax,&idmaxall,1,MPI_LMP_TAGINT,MPI_MAX,world);
 
   if (idminall != 1 || idmaxall != natoms) return 0;
   return 1;
 }
 
 /* ----------------------------------------------------------------------
    count and return words in a single line
    make copy of line before using strtok so as not to change line
    trim anything from '#' onward
 ------------------------------------------------------------------------- */
 
 int Atom::count_words(const char *line)
 {
   int n = strlen(line) + 1;
   char *copy;
   memory->create(copy,n,"atom:copy");
   strcpy(copy,line);
 
   char *ptr;
   if ((ptr = strchr(copy,'#'))) *ptr = '\0';
 
   if (strtok(copy," \t\n\r\f") == NULL) {
     memory->destroy(copy);
     return 0;
   }
   n = 1;
   while (strtok(NULL," \t\n\r\f")) n++;
 
   memory->destroy(copy);
   return n;
 }
 
 /* ----------------------------------------------------------------------
    deallocate molecular topology arrays
    done before realloc with (possibly) new 2nd dimension set to
      correctly initialized per-atom values, e.g. bond_per_atom
    needs to be called whenever 2nd dimensions are changed
      and these arrays are already pre-allocated,
      e.g. due to grow(1) in create_avec()
 ------------------------------------------------------------------------- */
 
 void Atom::deallocate_topology()
 {
   memory->destroy(atom->bond_type);
   memory->destroy(atom->bond_atom);
   atom->bond_type = NULL;
   atom->bond_atom = NULL;
 
   memory->destroy(atom->angle_type);
   memory->destroy(atom->angle_atom1);
   memory->destroy(atom->angle_atom2);
   memory->destroy(atom->angle_atom3);
   atom->angle_type = NULL;
   atom->angle_atom1 = atom->angle_atom2 = atom->angle_atom3 = NULL;
   
   memory->destroy(atom->dihedral_type);
   memory->destroy(atom->dihedral_atom1);
   memory->destroy(atom->dihedral_atom2);
   memory->destroy(atom->dihedral_atom3);
   memory->destroy(atom->dihedral_atom4);
   atom->dihedral_type = NULL;
   atom->dihedral_atom1 = atom->dihedral_atom2 = 
     atom->dihedral_atom3 = atom->dihedral_atom4 = NULL;
   
   memory->destroy(atom->improper_type);
   memory->destroy(atom->improper_atom1);
   memory->destroy(atom->improper_atom2);
   memory->destroy(atom->improper_atom3);
   memory->destroy(atom->improper_atom4);
   atom->improper_type = NULL;
   atom->improper_atom1 = atom->improper_atom2 = 
     atom->improper_atom3 = atom->improper_atom4 = NULL;
 }
 
 /* ----------------------------------------------------------------------
    unpack n lines from Atom section of data file
    call style-specific routine to parse line
 ------------------------------------------------------------------------- */
 
 void Atom::data_atoms(int n, char *buf)
 {
   int m,xptr,iptr;
   imageint imagedata;
   double xdata[3],lamda[3];
   double *coord;
   char *next;
 
   next = strchr(buf,'\n');
   *next = '\0';
   int nwords = count_words(buf);
   *next = '\n';
 
   if (nwords != avec->size_data_atom && nwords != avec->size_data_atom + 3)
     error->all(FLERR,"Incorrect atom format in data file");
 
   char **values = new char*[nwords];
 
   // set bounds for my proc
   // if periodic and I am lo/hi proc, adjust bounds by EPSILON
   // insures all data atoms will be owned even with round-off
 
   int triclinic = domain->triclinic;
 
   double epsilon[3];
   if (triclinic) epsilon[0] = epsilon[1] = epsilon[2] = EPSILON;
   else {
     epsilon[0] = domain->prd[0] * EPSILON;
     epsilon[1] = domain->prd[1] * EPSILON;
     epsilon[2] = domain->prd[2] * EPSILON;
   }
 
   double sublo[3],subhi[3];
   if (triclinic == 0) {
     sublo[0] = domain->sublo[0]; subhi[0] = domain->subhi[0];
     sublo[1] = domain->sublo[1]; subhi[1] = domain->subhi[1];
     sublo[2] = domain->sublo[2]; subhi[2] = domain->subhi[2];
   } else {
     sublo[0] = domain->sublo_lamda[0]; subhi[0] = domain->subhi_lamda[0];
     sublo[1] = domain->sublo_lamda[1]; subhi[1] = domain->subhi_lamda[1];
     sublo[2] = domain->sublo_lamda[2]; subhi[2] = domain->subhi_lamda[2];
   }
 
   if (domain->xperiodic) {
     if (comm->myloc[0] == 0) sublo[0] -= epsilon[0];
     if (comm->myloc[0] == comm->procgrid[0]-1) subhi[0] += epsilon[0];
   }
   if (domain->yperiodic) {
     if (comm->myloc[1] == 0) sublo[1] -= epsilon[1];
     if (comm->myloc[1] == comm->procgrid[1]-1) subhi[1] += epsilon[1];
   }
   if (domain->zperiodic) {
     if (comm->myloc[2] == 0) sublo[2] -= epsilon[2];
     if (comm->myloc[2] == comm->procgrid[2]-1) subhi[2] += epsilon[2];
   }
 
   // xptr = which word in line starts xyz coords
   // iptr = which word in line starts ix,iy,iz image flags
 
   xptr = avec->xcol_data - 1;
   int imageflag = 0;
   if (nwords > avec->size_data_atom) imageflag = 1;
   if (imageflag) iptr = nwords - 3;
 
   // loop over lines of atom data
   // tokenize the line into values
   // extract xyz coords and image flags
   // remap atom into simulation box
   // if atom is in my sub-domain, unpack its values
 
   for (int i = 0; i < n; i++) {
     next = strchr(buf,'\n');
 
     values[0] = strtok(buf," \t\n\r\f");
     if (values[0] == NULL)
       error->all(FLERR,"Incorrect atom format in data file");
     for (m = 1; m < nwords; m++) {
       values[m] = strtok(NULL," \t\n\r\f");
       if (values[m] == NULL)
         error->all(FLERR,"Incorrect atom format in data file");
     }
 
     if (imageflag)
       imagedata = ((imageint) (atoi(values[iptr]) + IMGMAX) & IMGMASK) |
         (((imageint) (atoi(values[iptr+1]) + IMGMAX) & IMGMASK) << IMGBITS) |
         (((imageint) (atoi(values[iptr+2]) + IMGMAX) & IMGMASK) << IMG2BITS);
     else imagedata = ((imageint) IMGMAX << IMG2BITS) |
            ((imageint) IMGMAX << IMGBITS) | IMGMAX;
     
     xdata[0] = atof(values[xptr]);
     xdata[1] = atof(values[xptr+1]);
     xdata[2] = atof(values[xptr+2]);
     domain->remap(xdata,imagedata);
     if (triclinic) {
       domain->x2lamda(xdata,lamda);
       coord = lamda;
     } else coord = xdata;
 
     if (coord[0] >= sublo[0] && coord[0] < subhi[0] &&
         coord[1] >= sublo[1] && coord[1] < subhi[1] &&
         coord[2] >= sublo[2] && coord[2] < subhi[2])
       avec->data_atom(xdata,imagedata,values);
 
     buf = next + 1;
   }
 
   delete [] values;
 }
 
 /* ----------------------------------------------------------------------
    unpack n lines from Velocity section of data file
    check that atom IDs are > 0 and <= map_tag_max
    call style-specific routine to parse line
 ------------------------------------------------------------------------- */
 
 void Atom::data_vels(int n, char *buf)
 {
   int j,m;
   tagint tagdata;
   char *next;
 
   next = strchr(buf,'\n');
   *next = '\0';
   int nwords = count_words(buf);
   *next = '\n';
 
   if (nwords != avec->size_data_vel)
     error->all(FLERR,"Incorrect velocity format in data file");
 
   char **values = new char*[nwords];
 
   // loop over lines of atom velocities
   // tokenize the line into values
   // if I own atom tag, unpack its values
 
   for (int i = 0; i < n; i++) {
     next = strchr(buf,'\n');
 
     values[0] = strtok(buf," \t\n\r\f");
     for (j = 1; j < nwords; j++)
       values[j] = strtok(NULL," \t\n\r\f");
 
     tagdata = ATOTAGINT(values[0]);
     if (tagdata <= 0 || tagdata > map_tag_max)
       error->one(FLERR,"Invalid atom ID in Velocities section of data file");
     if ((m = map(tagdata)) >= 0) avec->data_vel(m,&values[1]);
 
     buf = next + 1;
   }
 
   delete [] values;
 }
 
 /* ----------------------------------------------------------------------
    process N bonds read into buf from data files
    if count is non-NULL, just count bonds per atom
    else store them with atoms
    check that atom IDs are > 0 and <= map_tag_max
 ------------------------------------------------------------------------- */
 
 void Atom::data_bonds(int n, char *buf, int *count)
 {
   int m,tmp,itype;
   tagint atom1,atom2;
   char *next;
   int newton_bond = force->newton_bond;
 
   for (int i = 0; i < n; i++) {
     next = strchr(buf,'\n');
     *next = '\0';
     sscanf(buf,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT,
            &tmp,&itype,&atom1,&atom2);
     if (atom1 <= 0 || atom1 > map_tag_max ||
         atom2 <= 0 || atom2 > map_tag_max)
       error->one(FLERR,"Invalid atom ID in Bonds section of data file");
     if (itype <= 0 || itype > nbondtypes)
       error->one(FLERR,"Invalid bond type in Bonds section of data file");
     if ((m = map(atom1)) >= 0) {
       if (count) count[m]++;
       else {
         bond_type[m][num_bond[m]] = itype;
         bond_atom[m][num_bond[m]] = atom2;
         num_bond[m]++;
       }
     }
     if (newton_bond == 0) {
       if ((m = map(atom2)) >= 0) {
         if (count) count[m]++;
         else {
           bond_type[m][num_bond[m]] = itype;
           bond_atom[m][num_bond[m]] = atom1;
           num_bond[m]++;
         }
       }
     }
     buf = next + 1;
   }
 }
 
 /* ----------------------------------------------------------------------
    process N angles read into buf from data files
    if count is non-NULL, just count angles per atom
    else store them with atoms
    check that atom IDs are > 0 and <= map_tag_max
 ------------------------------------------------------------------------- */
 
 void Atom::data_angles(int n, char *buf, int *count)
 {
   int m,tmp,itype;
   tagint atom1,atom2,atom3;
   char *next;
   int newton_bond = force->newton_bond;
 
   for (int i = 0; i < n; i++) {
     next = strchr(buf,'\n');
     *next = '\0';
     sscanf(buf,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT,
            &tmp,&itype,&atom1,&atom2,&atom3);
     if (atom1 <= 0 || atom1 > map_tag_max ||
         atom2 <= 0 || atom2 > map_tag_max ||
         atom3 <= 0 || atom3 > map_tag_max)
       error->one(FLERR,"Invalid atom ID in Angles section of data file");
     if (itype <= 0 || itype > nangletypes)
       error->one(FLERR,"Invalid angle type in Angles section of data file");
     if ((m = map(atom2)) >= 0) {
       if (count) count[m]++;
       else {
         angle_type[m][num_angle[m]] = itype;
         angle_atom1[m][num_angle[m]] = atom1;
         angle_atom2[m][num_angle[m]] = atom2;
         angle_atom3[m][num_angle[m]] = atom3;
         num_angle[m]++;
       }
     }
     if (newton_bond == 0) {
       if ((m = map(atom1)) >= 0) {
         if (count) count[m]++;
         else {
           angle_type[m][num_angle[m]] = itype;
           angle_atom1[m][num_angle[m]] = atom1;
           angle_atom2[m][num_angle[m]] = atom2;
           angle_atom3[m][num_angle[m]] = atom3;
           num_angle[m]++;
         }
       }
       if ((m = map(atom3)) >= 0) {
         if (count) count[m]++;
         else {
           angle_type[m][num_angle[m]] = itype;
           angle_atom1[m][num_angle[m]] = atom1;
           angle_atom2[m][num_angle[m]] = atom2;
           angle_atom3[m][num_angle[m]] = atom3;
           num_angle[m]++;
         }
       }
     }
     buf = next + 1;
   }
 }
 
 /* ----------------------------------------------------------------------
    process N dihedrals read into buf from data files
    if count is non-NULL, just count diihedrals per atom
    else store them with atoms
    check that atom IDs are > 0 and <= map_tag_max
 ------------------------------------------------------------------------- */
 
 void Atom::data_dihedrals(int n, char *buf, int *count)
 {
   int m,tmp,itype;
   tagint atom1,atom2,atom3,atom4;
   char *next;
   int newton_bond = force->newton_bond;
 
   for (int i = 0; i < n; i++) {
     next = strchr(buf,'\n');
     *next = '\0';
     sscanf(buf,"%d %d " 
            TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT,
            &tmp,&itype,&atom1,&atom2,&atom3,&atom4);
     if (atom1 <= 0 || atom1 > map_tag_max ||
         atom2 <= 0 || atom2 > map_tag_max ||
         atom3 <= 0 || atom3 > map_tag_max ||
         atom4 <= 0 || atom4 > map_tag_max)
       error->one(FLERR,"Invalid atom ID in Dihedrals section of data file");
     if (itype <= 0 || itype > ndihedraltypes)
       error->one(FLERR,
                  "Invalid dihedral type in Dihedrals section of data file");
     if ((m = map(atom2)) >= 0) {
       if (count) count[m]++;
       else {
         dihedral_type[m][num_dihedral[m]] = itype;
         dihedral_atom1[m][num_dihedral[m]] = atom1;
         dihedral_atom2[m][num_dihedral[m]] = atom2;
         dihedral_atom3[m][num_dihedral[m]] = atom3;
         dihedral_atom4[m][num_dihedral[m]] = atom4;
         num_dihedral[m]++;
       }
     }
     if (newton_bond == 0) {
       if ((m = map(atom1)) >= 0) {
         if (count) count[m]++;
         else {
           dihedral_type[m][num_dihedral[m]] = itype;
           dihedral_atom1[m][num_dihedral[m]] = atom1;
           dihedral_atom2[m][num_dihedral[m]] = atom2;
           dihedral_atom3[m][num_dihedral[m]] = atom3;
           dihedral_atom4[m][num_dihedral[m]] = atom4;
           num_dihedral[m]++;
         }
       }
       if ((m = map(atom3)) >= 0) {
         if (count) count[m]++;
         else {
           dihedral_type[m][num_dihedral[m]] = itype;
           dihedral_atom1[m][num_dihedral[m]] = atom1;
           dihedral_atom2[m][num_dihedral[m]] = atom2;
           dihedral_atom3[m][num_dihedral[m]] = atom3;
           dihedral_atom4[m][num_dihedral[m]] = atom4;
           num_dihedral[m]++;
         }
       }
       if ((m = map(atom4)) >= 0) {
         if (count) count[m]++;
         else {
           dihedral_type[m][num_dihedral[m]] = itype;
           dihedral_atom1[m][num_dihedral[m]] = atom1;
           dihedral_atom2[m][num_dihedral[m]] = atom2;
           dihedral_atom3[m][num_dihedral[m]] = atom3;
           dihedral_atom4[m][num_dihedral[m]] = atom4;
           num_dihedral[m]++;
         }
       }
     }
     buf = next + 1;
   }
 }
 
 /* ----------------------------------------------------------------------
    process N impropers read into buf from data files
    if count is non-NULL, just count impropers per atom
    else store them with atoms
    check that atom IDs are > 0 and <= map_tag_max
 ------------------------------------------------------------------------- */
 
 void Atom::data_impropers(int n, char *buf, int *count)
 {
   int m,tmp,itype;
   tagint atom1,atom2,atom3,atom4;
   char *next;
   int newton_bond = force->newton_bond;
 
   for (int i = 0; i < n; i++) {
     next = strchr(buf,'\n');
     *next = '\0';
     sscanf(buf,"%d %d " 
            TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT,
            &tmp,&itype,&atom1,&atom2,&atom3,&atom4);
     if (atom1 <= 0 || atom1 > map_tag_max ||
         atom2 <= 0 || atom2 > map_tag_max ||
         atom3 <= 0 || atom3 > map_tag_max ||
         atom4 <= 0 || atom4 > map_tag_max)
       error->one(FLERR,"Invalid atom ID in Impropers section of data file");
     if (itype <= 0 || itype > nimpropertypes)
       error->one(FLERR,
                  "Invalid improper type in Impropers section of data file");
     if ((m = map(atom2)) >= 0) {
       if (count) count[m]++;
       else {
         improper_type[m][num_improper[m]] = itype;
         improper_atom1[m][num_improper[m]] = atom1;
         improper_atom2[m][num_improper[m]] = atom2;
         improper_atom3[m][num_improper[m]] = atom3;
         improper_atom4[m][num_improper[m]] = atom4;
         num_improper[m]++;
       }
     }
     if (newton_bond == 0) {
       if ((m = map(atom1)) >= 0) {
         if (count) count[m]++;
         else {
           improper_type[m][num_improper[m]] = itype;
           improper_atom1[m][num_improper[m]] = atom1;
           improper_atom2[m][num_improper[m]] = atom2;
           improper_atom3[m][num_improper[m]] = atom3;
           improper_atom4[m][num_improper[m]] = atom4;
           num_improper[m]++;
         }
       }
       if ((m = map(atom3)) >= 0) {
         if (count) count[m]++;
         else {
           improper_type[m][num_improper[m]] = itype;
           improper_atom1[m][num_improper[m]] = atom1;
           improper_atom2[m][num_improper[m]] = atom2;
           improper_atom3[m][num_improper[m]] = atom3;
           improper_atom4[m][num_improper[m]] = atom4;
           num_improper[m]++;
         }
       }
       if ((m = map(atom4)) >= 0) {
         if (count) count[m]++;
         else {
           improper_type[m][num_improper[m]] = itype;
           improper_atom1[m][num_improper[m]] = atom1;
           improper_atom2[m][num_improper[m]] = atom2;
           improper_atom3[m][num_improper[m]] = atom3;
           improper_atom4[m][num_improper[m]] = atom4;
           num_improper[m]++;
         }
       }
     }
     buf = next + 1;
   }
 }
 
 /* ----------------------------------------------------------------------
    unpack n lines from atom-style specific section of data file
    check that atom IDs are > 0 and <= map_tag_max
    call style-specific routine to parse line
 ------------------------------------------------------------------------- */
 
 void Atom::data_bonus(int n, char *buf, AtomVec *avec_bonus)
 {
   int j,m,tagdata;
   char *next;
 
   next = strchr(buf,'\n');
   *next = '\0';
   int nwords = count_words(buf);
   *next = '\n';
 
   if (nwords != avec_bonus->size_data_bonus)
     error->all(FLERR,"Incorrect bonus data format in data file");
 
   char **values = new char*[nwords];
 
   // loop over lines of bonus atom data
   // tokenize the line into values
   // if I own atom tag, unpack its values
 
   for (int i = 0; i < n; i++) {
     next = strchr(buf,'\n');
 
     values[0] = strtok(buf," \t\n\r\f");
     for (j = 1; j < nwords; j++)
       values[j] = strtok(NULL," \t\n\r\f");
 
     tagdata = ATOTAGINT(values[0]);
     if (tagdata <= 0 || tagdata > map_tag_max)
       error->one(FLERR,"Invalid atom ID in Bonus section of data file");
 
     // ok to call child's data_atom_bonus() method thru parent avec_bonus,
     // since data_bonus() was called with child ptr, and method is virtual
 
     if ((m = map(tagdata)) >= 0) avec_bonus->data_atom_bonus(m,&values[1]);
 
     buf = next + 1;
   }
 
   delete [] values;
 }
 
 /* ----------------------------------------------------------------------
    unpack n lines from atom-style specific section of data file
    check that atom IDs are > 0 and <= map_tag_max
    call style-specific routine to parse line
 ------------------------------------------------------------------------- */
 
 void Atom::data_bodies(int n, char *buf, AtomVecBody *avec_body)
 {
   int j,m,tagdata,ninteger,ndouble;
 
   char **ivalues = new char*[10*MAXBODY];
   char **dvalues = new char*[10*MAXBODY];
 
   // loop over lines of body data
   // tokenize the lines into ivalues and dvalues
   // if I own atom tag, unpack its values
 
   for (int i = 0; i < n; i++) {
     if (i == 0) tagdata = ATOTAGINT(strtok(buf," \t\n\r\f"));
     else tagdata = ATOTAGINT(strtok(NULL," \t\n\r\f"));
     ninteger = atoi(strtok(NULL," \t\n\r\f"));
     ndouble = atoi(strtok(NULL," \t\n\r\f"));
 
     for (j = 0; j < ninteger; j++)
       ivalues[j] = strtok(NULL," \t\n\r\f");
     for (j = 0; j < ndouble; j++)
       dvalues[j] = strtok(NULL," \t\n\r\f");
 
     if (tagdata <= 0 || tagdata > map_tag_max)
       error->one(FLERR,"Invalid atom ID in Bodies section of data file");
 
     if ((m = map(tagdata)) >= 0)
       avec_body->data_body(m,ninteger,ndouble,ivalues,dvalues);
   }
 
   delete [] ivalues;
   delete [] dvalues;
 }
 
 /* ----------------------------------------------------------------------
    allocate arrays of length ntypes
    only done after ntypes is set
 ------------------------------------------------------------------------- */
 
 void Atom::allocate_type_arrays()
 {
   if (avec->mass_type) {
     mass = new double[ntypes+1];
     mass_setflag = new int[ntypes+1];
     for (int itype = 1; itype <= ntypes; itype++) mass_setflag[itype] = 0;
   }
 }
 
 /* ----------------------------------------------------------------------
    set a mass and flag it as set
    called from reading of data file
 ------------------------------------------------------------------------- */
 
 void Atom::set_mass(const char *str)
 {
   if (mass == NULL) error->all(FLERR,"Cannot set mass for this atom style");
 
   int itype;
   double mass_one;
   int n = sscanf(str,"%d %lg",&itype,&mass_one);
   if (n != 2) error->all(FLERR,"Invalid mass line in data file");
 
   if (itype < 1 || itype > ntypes)
     error->all(FLERR,"Invalid type for mass set");
 
   mass[itype] = mass_one;
   mass_setflag[itype] = 1;
 
   if (mass[itype] <= 0.0) error->all(FLERR,"Invalid mass value");
 }
 
 /* ----------------------------------------------------------------------
    set a mass and flag it as set
    called from EAM pair routine
 ------------------------------------------------------------------------- */
 
 void Atom::set_mass(int itype, double value)
 {
   if (mass == NULL) error->all(FLERR,"Cannot set mass for this atom style");
   if (itype < 1 || itype > ntypes)
     error->all(FLERR,"Invalid type for mass set");
 
   mass[itype] = value;
   mass_setflag[itype] = 1;
 
   if (mass[itype] <= 0.0) error->all(FLERR,"Invalid mass value");
 }
 
 /* ----------------------------------------------------------------------
    set one or more masses and flag them as set
    called from reading of input script
 ------------------------------------------------------------------------- */
 
 void Atom::set_mass(int narg, char **arg)
 {
   if (mass == NULL) error->all(FLERR,"Cannot set mass for this atom style");
 
   int lo,hi;
   force->bounds(arg[0],ntypes,lo,hi);
   if (lo < 1 || hi > ntypes) error->all(FLERR,"Invalid type for mass set");
 
   for (int itype = lo; itype <= hi; itype++) {
     mass[itype] = atof(arg[1]);
     mass_setflag[itype] = 1;
 
     if (mass[itype] <= 0.0) error->all(FLERR,"Invalid mass value");
   }
 }
 
 /* ----------------------------------------------------------------------
    set all masses as read in from restart file
 ------------------------------------------------------------------------- */
 
 void Atom::set_mass(double *values)
 {
   for (int itype = 1; itype <= ntypes; itype++) {
     mass[itype] = values[itype];
     mass_setflag[itype] = 1;
   }
 }
 
 /* ----------------------------------------------------------------------
    check that all masses have been set
 ------------------------------------------------------------------------- */
 
 void Atom::check_mass()
 {
   if (mass == NULL) return;
   for (int itype = 1; itype <= ntypes; itype++)
     if (mass_setflag[itype] == 0) error->all(FLERR,"All masses are not set");
 }
 
 /* ----------------------------------------------------------------------
    check that radii of all particles of itype are the same
    return 1 if true, else return 0
    also return the radius value for that type
 ------------------------------------------------------------------------- */
 
 int Atom::radius_consistency(int itype, double &rad)
 {
   double value = -1.0;
   int flag = 0;
   for (int i = 0; i < nlocal; i++) {
     if (type[i] != itype) continue;
     if (value < 0.0) value = radius[i];
     else if (value != radius[i]) flag = 1;
   }
 
   int flagall;
   MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
   if (flagall) return 0;
 
   MPI_Allreduce(&value,&rad,1,MPI_DOUBLE,MPI_MAX,world);
   return 1;
 }
 
 /* ----------------------------------------------------------------------
    check that shape of all particles of itype are the same
    return 1 if true, else return 0
    also return the 3 shape params for itype
 ------------------------------------------------------------------------- */
 
 int Atom::shape_consistency(int itype,
                             double &shapex, double &shapey, double &shapez)
 {
   double zero[3] = {0.0, 0.0, 0.0};
   double one[3] = {-1.0, -1.0, -1.0};
   double *shape;
 
   AtomVecEllipsoid *avec_ellipsoid =
     (AtomVecEllipsoid *) style_match("ellipsoid");
   AtomVecEllipsoid::Bonus *bonus = avec_ellipsoid->bonus;
 
   int flag = 0;
   for (int i = 0; i < nlocal; i++) {
     if (type[i] != itype) continue;
     if (ellipsoid[i] < 0) shape = zero;
     else shape = bonus[ellipsoid[i]].shape;
 
     if (one[0] < 0.0) {
       one[0] = shape[0];
       one[1] = shape[1];
       one[2] = shape[2];
     } else if (one[0] != shape[0] || one[1] != shape[1] || one[2] != shape[2])
       flag = 1;
   }
 
   int flagall;
   MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
   if (flagall) return 0;
 
   double oneall[3];
   MPI_Allreduce(one,oneall,3,MPI_DOUBLE,MPI_MAX,world);
   shapex = oneall[0];
   shapey = oneall[1];
   shapez = oneall[2];
   return 1;
 }
 
 /* ----------------------------------------------------------------------
    add a new molecule template = set of molecules
 ------------------------------------------------------------------------- */
 
 void Atom::add_molecule(int narg, char **arg)
 {
   if (narg < 2) error->all(FLERR,"Illegal molecule command");
   if (find_molecule(arg[0]) >= 0) 
     error->all(FLERR,"Reuse of molecule template ID");
 
   int nprevious = nmolecule;
   nmolecule += narg-1;
   molecules = (Molecule **)
     memory->srealloc(molecules,nmolecule*sizeof(Molecule *),"atom::molecules");
 
   for (int i = 1; i < narg; i++) {
     molecules[nprevious] = new Molecule(lmp,arg[0],arg[i]);
     if (i == 1) molecules[nprevious]->nset = narg-1;
     else molecules[nprevious]->nset = 0;
     nprevious++;
   }
 }
 
 /* ----------------------------------------------------------------------
    find first molecule in set with template ID
    return -1 if does not exist
 ------------------------------------------------------------------------- */
 
 int Atom::find_molecule(char *id)
 {
   int imol;
   for (imol = 0; imol < nmolecule; imol++)
     if (strcmp(id,molecules[imol]->id) == 0) return imol;
   return -1;
 }
 
 /* ----------------------------------------------------------------------
    add info to current atom ilocal from molecule template onemol and its iatom
    offset = atom ID preceeding IDs of atoms in this molecule
    called by fixes and commands that add molecules
 ------------------------------------------------------------------------- */
 
 void Atom::add_molecule_atom(Molecule *onemol, int iatom,
                              int ilocal, tagint offset)
 {
   if (onemol->qflag && q_flag) q[ilocal] = onemol->q[iatom];
   if (onemol->radiusflag && radius_flag) radius[ilocal] = onemol->radius[iatom];
   if (onemol->rmassflag && rmass_flag) rmass[ilocal] = onemol->rmass[iatom];
   else if (rmass_flag) 
     rmass[ilocal] = 4.0*MY_PI/3.0 *
       radius[ilocal]*radius[ilocal]*radius[ilocal];
 
   if (molecular != 1) return;
 
   // add bond topology info
   // for molecular atom styles, but not atom style template
 
   if (avec->bonds_allow) {
     num_bond[ilocal] = onemol->num_bond[iatom];
     for (int i = 0; i < num_bond[ilocal]; i++) {
       bond_type[ilocal][i] = onemol->bond_type[iatom][i];
       bond_atom[ilocal][i] = onemol->bond_atom[iatom][i] + offset;
     }
   }
 
   if (avec->angles_allow) {
     num_angle[ilocal] = onemol->num_angle[iatom];
     for (int i = 0; i < num_angle[ilocal]; i++) {
       angle_type[ilocal][i] = onemol->angle_type[iatom][i];
       angle_atom1[ilocal][i] = onemol->angle_atom1[iatom][i] + offset;
       angle_atom2[ilocal][i] = onemol->angle_atom2[iatom][i] + offset;
       angle_atom3[ilocal][i] = onemol->angle_atom3[iatom][i] + offset;
     }
   }
 
   if (avec->dihedrals_allow) {
     num_dihedral[ilocal] = onemol->num_dihedral[iatom];
     for (int i = 0; i < num_dihedral[ilocal]; i++) {
       dihedral_type[ilocal][i] = onemol->dihedral_type[iatom][i];
       dihedral_atom1[ilocal][i] = onemol->dihedral_atom1[iatom][i] + offset;
       dihedral_atom2[ilocal][i] = onemol->dihedral_atom2[iatom][i] + offset;
       dihedral_atom3[ilocal][i] = onemol->dihedral_atom3[iatom][i] + offset;
       dihedral_atom4[ilocal][i] = onemol->dihedral_atom4[iatom][i] + offset;
     }
   }
 
   if (avec->impropers_allow) {
     num_improper[ilocal] = onemol->num_improper[iatom];
     for (int i = 0; i < num_improper[ilocal]; i++) {
       improper_type[ilocal][i] = onemol->improper_type[iatom][i];
       improper_atom1[ilocal][i] = onemol->improper_atom1[iatom][i] + offset;
       improper_atom2[ilocal][i] = onemol->improper_atom2[iatom][i] + offset;
       improper_atom3[ilocal][i] = onemol->improper_atom3[iatom][i] + offset;
       improper_atom4[ilocal][i] = onemol->improper_atom4[iatom][i] + offset;
     }
   }
 
   if (onemol->specialflag) {
     nspecial[ilocal][0] = onemol->nspecial[iatom][0];
     nspecial[ilocal][1] = onemol->nspecial[iatom][1];
     int n = nspecial[ilocal][2] = onemol->nspecial[iatom][2];
     for (int i = 0; i < n; i++)
       special[ilocal][i] = onemol->special[iatom][i] + offset;
   }
 }
 
 /* ----------------------------------------------------------------------
    reorder owned atoms so those in firstgroup appear first
    called by comm->exchange() if atom_modify first group is set
    only owned atoms exist at this point, no ghost atoms
 ------------------------------------------------------------------------- */
 
 void Atom::first_reorder()
 {
   // insure there is one extra atom location at end of arrays for swaps
 
   if (nlocal == nmax) avec->grow(0);
 
   // loop over owned atoms
   // nfirst = index of first atom not in firstgroup
   // when find firstgroup atom out of place, swap it with atom nfirst
 
   int bitmask = group->bitmask[firstgroup];
   nfirst = 0;
   while (nfirst < nlocal && mask[nfirst] & bitmask) nfirst++;
 
   for (int i = 0; i < nlocal; i++) {
     if (mask[i] & bitmask && i > nfirst) {
       avec->copy(i,nlocal,0);
       avec->copy(nfirst,i,0);
       avec->copy(nlocal,nfirst,0);
       while (nfirst < nlocal && mask[nfirst] & bitmask) nfirst++;
     }
   }
 }
 
 /* ----------------------------------------------------------------------
    perform spatial sort of atoms within my sub-domain
    always called between comm->exchange() and comm->borders()
    don't have to worry about clearing/setting atom->map since done in comm
 ------------------------------------------------------------------------- */
 
 void Atom::sort()
 {
   int i,m,n,ix,iy,iz,ibin,empty;
 
   // set next timestep for sorting to take place
 
   nextsort = (update->ntimestep/sortfreq)*sortfreq + sortfreq;
 
   // download data from GPU if necessary
 
   if (lmp->cuda && !lmp->cuda->oncpu) lmp->cuda->downloadAll();
 
   // re-setup sort bins if needed
 
   if (domain->box_change) setup_sort_bins();
   if (nbins == 1) return;
 
   // reallocate per-atom vectors if needed
 
   if (nlocal > maxnext) {
     memory->destroy(next);
     memory->destroy(permute);
     maxnext = atom->nmax;
     memory->create(next,maxnext,"atom:next");
     memory->create(permute,maxnext,"atom:permute");
   }
 
   // insure there is one extra atom location at end of arrays for swaps
 
   if (nlocal == nmax) avec->grow(0);
 
   // bin atoms in reverse order so linked list will be in forward order
 
   for (i = 0; i < nbins; i++) binhead[i] = -1;
 
   for (i = nlocal-1; i >= 0; i--) {
     ix = static_cast<int> ((x[i][0]-bboxlo[0])*bininvx);
     iy = static_cast<int> ((x[i][1]-bboxlo[1])*bininvy);
     iz = static_cast<int> ((x[i][2]-bboxlo[2])*bininvz);
     ix = MAX(ix,0);
     iy = MAX(iy,0);
     iz = MAX(iz,0);
     ix = MIN(ix,nbinx-1);
     iy = MIN(iy,nbiny-1);
     iz = MIN(iz,nbinz-1);
     ibin = iz*nbiny*nbinx + iy*nbinx + ix;
     next[i] = binhead[ibin];
     binhead[ibin] = i;
   }
 
   // permute = desired permutation of atoms
   // permute[I] = J means Ith new atom will be Jth old atom
 
   n = 0;
   for (m = 0; m < nbins; m++) {
     i = binhead[m];
     while (i >= 0) {
       permute[n++] = i;
       i = next[i];
     }
   }
 
   // current = current permutation, just reuse next vector
   // current[I] = J means Ith current atom is Jth old atom
 
   int *current = next;
   for (i = 0; i < nlocal; i++) current[i] = i;
 
   // reorder local atom list, when done, current = permute
   // perform "in place" using copy() to extra atom location at end of list
   // inner while loop processes one cycle of the permutation
   // copy before inner-loop moves an atom to end of atom list
   // copy after inner-loop moves atom at end of list back into list
   // empty = location in atom list that is currently empty
 
   for (i = 0; i < nlocal; i++) {
     if (current[i] == permute[i]) continue;
     avec->copy(i,nlocal,0);
     empty = i;
     while (permute[empty] != i) {
       avec->copy(permute[empty],empty,0);
       empty = current[empty] = permute[empty];
     }
     avec->copy(nlocal,empty,0);
     current[empty] = permute[empty];
   }
 
   // upload data back to GPU if necessary
 
   if (lmp->cuda && !lmp->cuda->oncpu) lmp->cuda->uploadAll();
 
   // sanity check that current = permute
 
   //int flag = 0;
   //for (i = 0; i < nlocal; i++)
   //  if (current[i] != permute[i]) flag = 1;
   //int flagall;
   //MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
   //if (flagall) error->all(FLERR,"Atom sort did not operate correctly");
 }
 
 /* ----------------------------------------------------------------------
    setup bins for spatial sorting of atoms
 ------------------------------------------------------------------------- */
 
 void Atom::setup_sort_bins()
 {
   // binsize:
   // user setting if explicitly set
   // 1/2 of neighbor cutoff for non-CUDA
   // CUDA_CHUNK atoms/proc for CUDA
   // check if neighbor cutoff = 0.0
 
   double binsize;
   if (userbinsize > 0.0) binsize = userbinsize;
   else if (!lmp->cuda) binsize = 0.5 * neighbor->cutneighmax;
   else {
     if (domain->dimension == 3) {
       double vol = (domain->boxhi[0]-domain->boxlo[0]) *
         (domain->boxhi[1]-domain->boxlo[1]) *
         (domain->boxhi[2]-domain->boxlo[2]);
       binsize = pow(1.0*CUDA_CHUNK/natoms*vol,1.0/3.0);
     } else {
       double area = (domain->boxhi[0]-domain->boxlo[0]) *
         (domain->boxhi[1]-domain->boxlo[1]);
       binsize = pow(1.0*CUDA_CHUNK/natoms*area,1.0/2.0);
     }
   }
   if (binsize == 0.0) error->all(FLERR,"Atom sorting has bin size = 0.0");
 
   double bininv = 1.0/binsize;
 
   // nbin xyz = local bins
   // bbox lo/hi = bounding box of my sub-domain
 
   if (domain->triclinic)
     domain->bbox(domain->sublo_lamda,domain->subhi_lamda,bboxlo,bboxhi);
   else {
     bboxlo[0] = domain->sublo[0];
     bboxlo[1] = domain->sublo[1];
     bboxlo[2] = domain->sublo[2];
     bboxhi[0] = domain->subhi[0];
     bboxhi[1] = domain->subhi[1];
     bboxhi[2] = domain->subhi[2];
   }
 
   nbinx = static_cast<int> ((bboxhi[0]-bboxlo[0]) * bininv);
   nbiny = static_cast<int> ((bboxhi[1]-bboxlo[1]) * bininv);
   nbinz = static_cast<int> ((bboxhi[2]-bboxlo[2]) * bininv);
   if (domain->dimension == 2) nbinz = 1;
 
   if (nbinx == 0) nbinx = 1;
   if (nbiny == 0) nbiny = 1;
   if (nbinz == 0) nbinz = 1;
 
   bininvx = nbinx / (bboxhi[0]-bboxlo[0]);
   bininvy = nbiny / (bboxhi[1]-bboxlo[1]);
   bininvz = nbinz / (bboxhi[2]-bboxlo[2]);
 
   if (1.0*nbinx*nbiny*nbinz > INT_MAX)
     error->one(FLERR,"Too many atom sorting bins");
 
   nbins = nbinx*nbiny*nbinz;
 
   // reallocate per-bin memory if needed
 
   if (nbins > maxbin) {
     memory->destroy(binhead);
     maxbin = nbins;
     memory->create(binhead,maxbin,"atom:binhead");
   }
 }
 
 /* ----------------------------------------------------------------------
    register a callback to a fix so it can manage atom-based arrays
    happens when fix is created
    flag = 0 for grow, 1 for restart, 2 for border comm
 ------------------------------------------------------------------------- */
 
 void Atom::add_callback(int flag)
 {
   int ifix;
 
   // find the fix
   // if find NULL ptr:
   //   it's this one, since it is being replaced and has just been deleted
   //   at this point in re-creation
   // if don't find NULL ptr:
   //   i is set to nfix = new one currently being added at end of list
 
   for (ifix = 0; ifix < modify->nfix; ifix++)
     if (modify->fix[ifix] == NULL) break;
 
   // add callback to lists, reallocating if necessary
 
   if (flag == 0) {
     if (nextra_grow == nextra_grow_max) {
       nextra_grow_max += DELTA;
       memory->grow(extra_grow,nextra_grow_max,"atom:extra_grow");
     }
     extra_grow[nextra_grow] = ifix;
     nextra_grow++;
   } else if (flag == 1) {
     if (nextra_restart == nextra_restart_max) {
       nextra_restart_max += DELTA;
       memory->grow(extra_restart,nextra_restart_max,"atom:extra_restart");
     }
     extra_restart[nextra_restart] = ifix;
     nextra_restart++;
   } else if (flag == 2) {
     if (nextra_border == nextra_border_max) {
       nextra_border_max += DELTA;
       memory->grow(extra_border,nextra_border_max,"atom:extra_border");
     }
     extra_border[nextra_border] = ifix;
     nextra_border++;
   }
 }
 
 /* ----------------------------------------------------------------------
    unregister a callback to a fix
    happens when fix is deleted, called by its destructor
    flag = 0 for grow, 1 for restart
 ------------------------------------------------------------------------- */
 
 void Atom::delete_callback(const char *id, int flag)
 {
   int ifix;
   for (ifix = 0; ifix < modify->nfix; ifix++)
     if (strcmp(id,modify->fix[ifix]->id) == 0) break;
 
   // compact the list of callbacks
 
   if (flag == 0) {
     int match;
     for (match = 0; match < nextra_grow; match++)
       if (extra_grow[match] == ifix) break;
     for (int i = match; i < nextra_grow-1; i++)
       extra_grow[i] = extra_grow[i+1];
     nextra_grow--;
 
   } else if (flag == 1) {
     int match;
     for (match = 0; match < nextra_restart; match++)
       if (extra_restart[match] == ifix) break;
     for (int i = match; i < nextra_restart-1; i++)
       extra_restart[i] = extra_restart[i+1];
     nextra_restart--;
 
   } else if (flag == 2) {
     int match;
     for (match = 0; match < nextra_border; match++)
       if (extra_border[match] == ifix) break;
     for (int i = match; i < nextra_border-1; i++)
       extra_border[i] = extra_border[i+1];
     nextra_border--;
   }
 }
 
 /* ----------------------------------------------------------------------
    decrement ptrs in callback lists to fixes beyond the deleted ifix
    happens after fix is deleted
 ------------------------------------------------------------------------- */
 
 void Atom::update_callback(int ifix)
 {
   for (int i = 0; i < nextra_grow; i++)
     if (extra_grow[i] > ifix) extra_grow[i]--;
   for (int i = 0; i < nextra_restart; i++)
     if (extra_restart[i] > ifix) extra_restart[i]--;
   for (int i = 0; i < nextra_border; i++)
     if (extra_border[i] > ifix) extra_border[i]--;
 }
 
 /* ----------------------------------------------------------------------
    find custom per-atom vector with name
    return index if found, and flag = 0/1 for int/double
    return -1 if not found
 ------------------------------------------------------------------------- */
 
 int Atom::find_custom(char *name, int &flag)
 {
   for (int i = 0; i < nivector; i++)
     if (iname[i] && strcmp(iname[i],name) == 0) {
       flag = 0;
       return i;
     }
 
   for (int i = 0; i < ndvector; i++)
     if (dname[i] && strcmp(dname[i],name) == 0) {
       flag = 1;
       return i;
     }
 
   return -1;
 }
 
 /* ----------------------------------------------------------------------
    add a custom variable with name of type flag = 0/1 for int/double
    assumes name does not already exist
    return index in ivector or dvector of its location
 ------------------------------------------------------------------------- */
 
 int Atom::add_custom(char *name, int flag)
 {
   int index;
 
   if (flag == 0) {
     index = nivector;
     nivector++;
     iname = (char **) memory->srealloc(iname,nivector*sizeof(char *),
                                        "atom:iname");
     int n = strlen(name) + 1;
     iname[index] = new char[n];
     strcpy(iname[index],name);
     ivector = (int **) memory->srealloc(ivector,nivector*sizeof(int *),
                                         "atom:ivector");
     memory->create(ivector[index],nmax,"atom:ivector");
   } else {
     index = ndvector;
     ndvector++;
     dname = (char **) memory->srealloc(dname,ndvector*sizeof(char *),
                                        "atom:dname");
     int n = strlen(name) + 1;
     dname[index] = new char[n];
     strcpy(dname[index],name);
     dvector = (double **) memory->srealloc(dvector,ndvector*sizeof(double *),
                                            "atom:dvector");
     memory->create(dvector[index],nmax,"atom:dvector");
   }
 
   return index;
 }
 
 /* ----------------------------------------------------------------------
    remove a custom variable of type flag = 0/1 for int/double at index
    free memory for vector and name and set ptrs to NULL
    ivector/dvector and iname/dname lists never shrink
 ------------------------------------------------------------------------- */
 
 void Atom::remove_custom(int flag, int index)
 {
   if (flag == 0) {
     memory->destroy(ivector[index]);
     ivector[index] = NULL;
     delete [] iname[index];
     iname[index] = NULL;
   } else {
     memory->destroy(dvector[index]);
     dvector[index] = NULL;
     delete [] dname[index];
     dname[index] = NULL;
   }
 }
 
 /* ----------------------------------------------------------------------
    return a pointer to a named internal variable
    if don't recognize name, return NULL
    customize by adding names
 ------------------------------------------------------------------------- */
 
 void *Atom::extract(char *name)
 {
   if (strcmp(name,"mass") == 0) return (void *) mass;
 
   if (strcmp(name,"id") == 0) return (void *) tag;
   if (strcmp(name,"type") == 0) return (void *) type;
   if (strcmp(name,"mask") == 0) return (void *) mask;
   if (strcmp(name,"image") == 0) return (void *) image;
   if (strcmp(name,"x") == 0) return (void *) x;
   if (strcmp(name,"v") == 0) return (void *) v;
   if (strcmp(name,"f") == 0) return (void *) f;
   if (strcmp(name,"molecule") == 0) return (void *) molecule;
   if (strcmp(name,"q") == 0) return (void *) q;
   if (strcmp(name,"mu") == 0) return (void *) mu;
   if (strcmp(name,"omega") == 0) return (void *) omega;
-  if (strcmp(name,"amgmom") == 0) return (void *) angmom;
+  if (strcmp(name,"angmom") == 0) return (void *) angmom;
   if (strcmp(name,"torque") == 0) return (void *) torque;
   if (strcmp(name,"radius") == 0) return (void *) radius;
   if (strcmp(name,"rmass") == 0) return (void *) rmass;
+  if (strcmp(name,"ellipsoid") == 0) return (void *) ellipsoid;
+  if (strcmp(name,"line") == 0) return (void *) line;
+  if (strcmp(name,"tri") == 0) return (void *) tri;
+
   if (strcmp(name,"vfrac") == 0) return (void *) vfrac;
   if (strcmp(name,"s0") == 0) return (void *) s0;
+  if (strcmp(name,"x0") == 0) return (void *) x0;
+
+  if (strcmp(name,"spin") == 0) return (void *) spin;
+  if (strcmp(name,"eradius") == 0) return (void *) eradius;
+  if (strcmp(name,"ervel") == 0) return (void *) ervel;
+  if (strcmp(name,"erforce") == 0) return (void *) erforce;
+  if (strcmp(name,"ervelforce") == 0) return (void *) ervelforce;
+  if (strcmp(name,"cs") == 0) return (void *) cs;
+  if (strcmp(name,"csforce") == 0) return (void *) csforce;
+  if (strcmp(name,"vforce") == 0) return (void *) vforce;
+  if (strcmp(name,"etag") == 0) return (void *) etag;
+
+  if (strcmp(name,"rho") == 0) return (void *) rho;
+  if (strcmp(name,"drho") == 0) return (void *) drho;
+  if (strcmp(name,"e") == 0) return (void *) e;
+  if (strcmp(name,"de") == 0) return (void *) de;
+  if (strcmp(name,"cv") == 0) return (void *) cv;
+  if (strcmp(name,"vest") == 0) return (void *) vest;
 
   return NULL;
 }
 
 /* ----------------------------------------------------------------------
    return # of bytes of allocated memory
    call to avec tallies per-atom vectors
    add in global to local mapping storage
 ------------------------------------------------------------------------- */
 
 bigint Atom::memory_usage()
 {
   memlength = DELTA_MEMSTR;
   memory->create(memstr,memlength,"atom:memstr");
   memstr[0] = '\0';
   bigint bytes = avec->memory_usage();
   memory->destroy(memstr);
 
   bytes += max_same*sizeof(int);
   if (map_style == 1)
     bytes += memory->usage(map_array,map_maxarray);
   else if (map_style == 2) {
     bytes += map_nbucket*sizeof(int);
     bytes += map_nhash*sizeof(HashElem);
   }
   if (maxnext) {
     bytes += memory->usage(next,maxnext);
     bytes += memory->usage(permute,maxnext);
   }
 
   return bytes;
 }
 
 /* ----------------------------------------------------------------------
    accumulate per-atom vec names in memstr, padded by spaces
    return 1 if padded str is not already in memlist, else 0
 ------------------------------------------------------------------------- */
 
 int Atom::memcheck(const char *str)
 {
   int n = strlen(str) + 3;
   char *padded = new char[n];
   strcpy(padded," ");
   strcat(padded,str);
   strcat(padded," ");
 
   if (strstr(memstr,padded)) {
     delete [] padded;
     return 0;
   }
 
   if (strlen(memstr) + n >= memlength) {
     memlength += DELTA_MEMSTR;
     memory->grow(memstr,memlength,"atom:memstr");
   }
 
   strcat(memstr,padded);
   delete [] padded;
   return 1;
 }
diff --git a/src/atom.h b/src/atom.h
index 5887b17c3..c6bebe88a 100644
--- a/src/atom.h
+++ b/src/atom.h
@@ -1,459 +1,468 @@
 /* -*- c++ -*- ----------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
 
    Copyright (2003) Sandia Corporation.  Under the terms of Contract
    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
    certain rights in this software.  This software is distributed under
    the GNU General Public License.
 
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
 #ifndef LMP_ATOM_H
 #define LMP_ATOM_H
 
 #include "pointers.h"
 
 namespace LAMMPS_NS {
 
 class Atom : protected Pointers {
  public:
   char *atom_style;
   class AtomVec *avec;
 
   // atom counts
 
   bigint natoms;                // total # of atoms in system, could be 0
                                 // natoms may not be current if atoms lost
   int nlocal,nghost;            // # of owned and ghost atoms on this proc
   int nmax;                     // max # of owned+ghost in arrays on this proc
   int tag_enable;               // 0/1 if atom ID tags are defined
   int molecular;                // 0 = atomic, 1 = standard molecular system,
                                 // 2 = molecule template system
 
   bigint nbonds,nangles,ndihedrals,nimpropers;
   int ntypes,nbondtypes,nangletypes,ndihedraltypes,nimpropertypes;
   int bond_per_atom,angle_per_atom,dihedral_per_atom,improper_per_atom;
   int extra_bond_per_atom,extra_angle_per_atom;
   int extra_dihedral_per_atom,extra_improper_per_atom;
 
   int firstgroup;               // store atoms in this group first, -1 if unset
   int nfirst;                   // # of atoms in first group on this proc
   char *firstgroupname;         // group-ID to store first, NULL if unset
 
   // per-atom arrays
   // customize by adding new array
 
   tagint *tag;
   int *type,*mask;
   imageint *image;
   double **x,**v,**f;
 
   tagint *molecule;
   int *molindex,*molatom;
+
   double *q,**mu;
   double **omega,**angmom,**torque;
-  double *radius,*rmass,*vfrac,*s0;
-  double **x0;
+  double *radius,*rmass;
   int *ellipsoid,*line,*tri,*body;
+
+  // PERI package
+
+  double *vfrac,*s0;
+  double **x0;
+
+  // USER-EFF and USER-AWPMD packages
+
   int *spin;
   double *eradius,*ervel,*erforce,*ervelforce;
   double *cs,*csforce,*vforce;
   int *etag;
-  double *rho,*drho;
-  double *e,*de;
+
+  // USER-SPH package
+
+  double *rho,*drho,*e,*de,*cv;
   double **vest;
-  double *cv;
 
   int **nspecial;               // 0,1,2 = cummulative # of 1-2,1-3,1-4 neighs
   tagint **special;             // IDs of 1-2,1-3,1-4 neighs of each atom
   int maxspecial;               // special[nlocal][maxspecial]
 
   int *num_bond;
   int **bond_type;
   tagint **bond_atom;
 
   int *num_angle;
   int **angle_type;
   tagint **angle_atom1,**angle_atom2,**angle_atom3;
 
   int *num_dihedral;
   int **dihedral_type;
   tagint **dihedral_atom1,**dihedral_atom2,**dihedral_atom3,**dihedral_atom4;
 
   int *num_improper;
   int **improper_type;
   tagint **improper_atom1,**improper_atom2,**improper_atom3,**improper_atom4;
 
   // custom arrays used by fix property/atom
 
   int **ivector;
   double **dvector;
   char **iname,**dname;
   int nivector,ndvector;
 
   // used by USER-CUDA to flag used per-atom arrays
 
   unsigned int datamask;
   unsigned int datamask_ext;
 
   // atom style and per-atom array existence flags
   // customize by adding new flag
 
   int sphere_flag,ellipsoid_flag,line_flag,tri_flag,body_flag;
   int peri_flag,electron_flag;
   int ecp_flag;
   int wavepacket_flag,sph_flag;
 
   int molecule_flag,molindex_flag,molatom_flag;
   int q_flag,mu_flag;
   int rmass_flag,radius_flag,omega_flag,torque_flag,angmom_flag;
   int vfrac_flag,spin_flag,eradius_flag,ervel_flag,erforce_flag;
   int cs_flag,csforce_flag,vforce_flag,ervelforce_flag,etag_flag;
   int rho_flag,e_flag,cv_flag,vest_flag;
 
   // Peridynamics scale factor, used by dump cfg
 
   double pdscale;
 
   // molecule templates
   // each template can be a set of consecutive molecules
   // each with same ID (stored in molecules)
   // 1st molecule in template stores nset = # in set
 
   int nmolecule;
   class Molecule **molecules;
 
   // extra peratom info in restart file destined for fix & diag
 
   double **extra;
 
   // per-type arrays
 
   double *mass;
   int *mass_setflag;
 
   // callback ptrs for atom arrays managed by fix classes
 
   int nextra_grow,nextra_restart,nextra_border;  // # of callbacks of each type
   int *extra_grow,*extra_restart,*extra_border;  // index of fix to callback to
   int nextra_grow_max,nextra_restart_max;        // size of callback lists
   int nextra_border_max;
   int nextra_store;
 
   int map_style;                  // style of atom map: 0=none, 1=array, 2=hash
   int map_user;                   // user selected style = same 0,1,2
   tagint map_tag_max;             // max atom ID that map() is setup for
 
   // spatial sorting of atoms
 
   int sortfreq;             // sort atoms every this many steps, 0 = off
   bigint nextsort;          // next timestep to sort on
 
   // indices of atoms with same ID
 
   int *sametag;      // sametag[I] = next atom with same ID, -1 if no more
 
   // functions
 
   Atom(class LAMMPS *);
   ~Atom();
 
   void settings(class Atom *);
   void create_avec(const char *, int, char **, char *suffix = NULL);
   class AtomVec *new_avec(const char *, char *, int &);
   void init();
   void setup();
 
   class AtomVec *style_match(const char *);
   void modify_params(int, char **);
   void tag_check();
   void tag_extend();
   int tag_consecutive();
 
   int parse_data(const char *);
   int count_words(const char *);
 
   void deallocate_topology();
 
   void data_atoms(int, char *);
   void data_vels(int, char *);
 
   void data_bonds(int, char *, int *);
   void data_angles(int, char *, int *);
   void data_dihedrals(int, char *, int *);
   void data_impropers(int, char *, int *);
 
   void data_bonus(int, char *, class AtomVec *);
   void data_bodies(int, char *, class AtomVecBody *);
 
   virtual void allocate_type_arrays();
   void set_mass(const char *);
   void set_mass(int, double);
   void set_mass(int, char **);
   void set_mass(double *);
   void check_mass();
 
   int radius_consistency(int, double &);
   int shape_consistency(int, double &, double &, double &);
 
   void add_molecule(int, char **);
   int find_molecule(char *);
   void add_molecule_atom(class Molecule *, int, int, tagint);
 
   void first_reorder();
   virtual void sort();
 
   void add_callback(int);
   void delete_callback(const char *, int);
   void update_callback(int);
 
   int find_custom(char *, int &);
   int add_custom(char *, int);
   void remove_custom(int, int);
 
   void *extract(char *);
 
   inline int* get_map_array() {return map_array;};
   inline int get_map_size() {return map_tag_max+1;};
 
   bigint memory_usage();
   int memcheck(const char *);
 
   // functions for global to local ID mapping
   // map lookup function inlined for efficiency
   // return -1 if no map defined
 
   inline int map(tagint global) {
     if (map_style == 1) return map_array[global];
     else if (map_style == 2) return map_find_hash(global);
     else return -1;
   };
 
   void map_init(int check = 1);
   void map_clear();
   void map_set();
   void map_one(tagint, int);
   int map_style_set();
   void map_delete();
   int map_find_hash(tagint);
 
  protected:
 
   // global to local ID mapping
 
   int *map_array;       // direct map via array that holds map_tag_max
   int map_maxarray;     // allocated size of map_array (1 larger than this)
 
   struct HashElem {     // hashed map
     tagint global;      // key to search on = global ID
     int local;          // value associated with key = local index
     int next;           // next entry in this bucket, -1 if last
   };
   int map_nhash;        // # of entries hash table can hold
   int map_nused;        // # of actual entries in hash table
   int map_free;         // ptr to 1st unused entry in hash table
   int map_nbucket;      // # of hash buckets
   int *map_bucket;      // ptr to 1st entry in each bucket
   HashElem *map_hash;   // hash table
 
   int max_same;         // allocated size of sametag
 
   // spatial sorting of atoms
 
   int nbins;                      // # of sorting bins
   int nbinx,nbiny,nbinz;          // bins in each dimension
   int maxbin;                     // max # of bins
   int maxnext;                    // max size of next,permute
   int *binhead;                   // 1st atom in each bin
   int *next;                      // next atom in bin
   int *permute;                   // permutation vector
   double userbinsize;             // requested sort bin size
   double bininvx,bininvy,bininvz; // inverse actual bin sizes
   double bboxlo[3],bboxhi[3];     // bounding box of my sub-domain
 
   int memlength;                  // allocated size of memstr
   char *memstr;                   // string of array names already counted
 
   void setup_sort_bins();
   int next_prime(int);
 };
 
 }
 
 #endif
 
 /* ERROR/WARNING messages:
 
 E: Atom IDs must be used for molecular systems
 
 Atom IDs are used to identify and find partner atoms in bonds.
 
 E: Invalid atom style
 
 The choice of atom style is unknown.
 
 E: Could not find atom_modify first group ID
 
 Self-explanatory.
 
 E: Illegal ... command
 
 Self-explanatory.  Check the input script syntax and compare to the
 documentation for the command.  You can use -echo screen as a
 command-line option when running LAMMPS to see the offending line.
 
 E: Atom_modify id command after simulation box is defined
 
 The atom_modify id command cannot be used after a read_data,
 read_restart, or create_box command.
 
 E: Atom_modify map command after simulation box is defined
 
 The atom_modify map command cannot be used after a read_data,
 read_restart, or create_box command.
 
 E: Atom_modify sort and first options cannot be used together
 
 Self-explanatory.
 
 E: Atom ID is negative
 
 Self-explanatory.
 
 E: Atom ID is too big
 
 The limit on atom IDs is set by the SMALLBIG, BIGBIG, SMALLSMALL
 setting in your Makefile.  See Section_start 2.2 of the manual for
 more details.
 
 E: Atom ID is zero
 
 Either all atoms IDs must be zero or none of them.
 
 E: Not all atom IDs are 0
 
 Either all atoms IDs must be zero or none of them.
 
 E: New atom IDs exceed maximum allowed ID
 
 See the setting for tagint in the src/lmptype.h file.
 
 E: Incorrect atom format in data file
 
 Number of values per atom line in the data file is not consistent with
 the atom style.
 
 E: Incorrect velocity format in data file
 
 Each atom style defines a format for the Velocity section
 of the data file.  The read-in lines do not match.
 
 E: Invalid atom ID in Velocities section of data file
 
 Atom IDs must be positive integers and within range of defined
 atoms.
 
 E: Invalid atom ID in Bonds section of data file
 
 Atom IDs must be positive integers and within range of defined
 atoms.
 
 E: Invalid bond type in Bonds section of data file
 
 Bond type must be positive integer and within range of specified bond
 types.
 
 E: Invalid atom ID in Angles section of data file
 
 Atom IDs must be positive integers and within range of defined
 atoms.
 
 E: Invalid angle type in Angles section of data file
 
 Angle type must be positive integer and within range of specified angle
 types.
 
 E: Invalid atom ID in Dihedrals section of data file
 
 Atom IDs must be positive integers and within range of defined
 atoms.
 
 E: Invalid dihedral type in Dihedrals section of data file
 
 Dihedral type must be positive integer and within range of specified
 dihedral types.
 
 E: Invalid atom ID in Impropers section of data file
 
 Atom IDs must be positive integers and within range of defined
 atoms.
 
 E: Invalid improper type in Impropers section of data file
 
 Improper type must be positive integer and within range of specified
 improper types.
 
 E: Incorrect bonus data format in data file
 
 See the read_data doc page for a description of how various kinds of
 bonus data must be formatted for certain atom styles.
 
 E: Invalid atom ID in Bonus section of data file
 
 Atom IDs must be positive integers and within range of defined
 atoms.
 
 E: Invalid atom ID in Bodies section of data file
 
 Atom IDs must be positive integers and within range of defined
 atoms.
 
 E: Cannot set mass for this atom style
 
 This atom style does not support mass settings for each atom type.
 Instead they are defined on a per-atom basis in the data file.
 
 E: Invalid mass line in data file
 
 Self-explanatory.
 
 E: Invalid type for mass set
 
 Mass command must set a type from 1-N where N is the number of atom
 types.
 
 E: Invalid mass value
 
 Self-explanatory.
 
 E: All masses are not set
 
 For atom styles that define masses for each atom type, all masses must
 be set in the data file or by the mass command before running a
 simulation.  They must also be set before using the velocity
 command.
 
 E: Reuse of molecule template ID
 
 The template IDs must be unique.
 
 E: Atom sort did not operate correctly
 
 This is an internal LAMMPS error.  Please report it to the
 developers.
 
 E: Atom sorting has bin size = 0.0
 
 The neighbor cutoff is being used as the bin size, but it is zero.
 Thus you must explicitly list a bin size in the atom_modify sort
 command or turn off sorting.
 
 E: Too many atom sorting bins
 
 This is likely due to an immense simulation box that has blown up
 to a large size.
 
 */
diff --git a/src/atom_vec.h b/src/atom_vec.h
index 8e9bb0f51..a58e7c6e9 100644
--- a/src/atom_vec.h
+++ b/src/atom_vec.h
@@ -1,160 +1,162 @@
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
 
    Copyright (2003) Sandia Corporation.  Under the terms of Contract
    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
    certain rights in this software.  This software is distributed under
    the GNU General Public License.
 
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
 #ifndef LMP_ATOM_VEC_H
 #define LMP_ATOM_VEC_H
 
 #include "stdio.h"
 #include "pointers.h"
 
 namespace LAMMPS_NS {
 
 class AtomVec : protected Pointers {
  public:
   int molecular;                       // 0 = atomic, 1 = molecular system
   int bonds_allow,angles_allow;        // 1 if bonds, angles are used
   int dihedrals_allow,impropers_allow; // 1 if dihedrals, impropers used
   int mass_type;                       // 1 if per-type masses
   int dipole_type;                     // 1 if per-type dipole moments
 
   int comm_x_only;                     // 1 if only exchange x in forward comm
   int comm_f_only;                     // 1 if only exchange f in reverse comm
 
   int size_forward;                    // # of values per atom in comm
   int size_reverse;                    // # in reverse comm
   int size_border;                     // # in border comm
   int size_velocity;                   // # of velocity based quantities
   int size_data_atom;                  // number of values in Atom line
   int size_data_vel;                   // number of values in Velocity line
   int size_data_bonus;                 // number of values in Bonus line
   int xcol_data;                       // column (1-N) where x is in Atom line
 
+  int forceclearflag;                  // 1 if has forceclear() method
   class Molecule **onemols;            // list of molecules for style template
   int nset;                            // # of molecules in list
 
   int cudable;                         // 1 if atom style is CUDA-enabled
   int kokkosable;                      // 1 if atom style is KOKKOS-enabled
   int *maxsend;                        // CUDA-specific variable
 
   int nargcopy;          // copy of command-line args for atom_style command
   char **argcopy;        // used when AtomVec is realloced (restart,replicate)
 
   AtomVec(class LAMMPS *);
   virtual ~AtomVec();
   void store_args(int, char **);
   virtual void process_args(int, char **);
   virtual void init();
 
   virtual void grow(int) = 0;
   virtual void grow_reset() = 0;
   virtual void copy(int, int, int) = 0;
   virtual void clear_bonus() {}
+  virtual void force_clear(int, size_t) {}
 
   virtual int pack_comm(int, int *, double *, int, int *) = 0;
   virtual int pack_comm_vel(int, int *, double *, int, int *) = 0;
   virtual int pack_comm_hybrid(int, int *, double *) {return 0;}
   virtual void unpack_comm(int, int, double *) = 0;
   virtual void unpack_comm_vel(int, int, double *) = 0;
   virtual int unpack_comm_hybrid(int, int, double *) {return 0;}
 
   virtual int pack_reverse(int, int, double *) = 0;
   virtual int pack_reverse_hybrid(int, int, double *) {return 0;}
   virtual void unpack_reverse(int, int *, double *) = 0;
   virtual int unpack_reverse_hybrid(int, int *, double *) {return 0;}
 
   virtual int pack_border(int, int *, double *, int, int *) = 0;
   virtual int pack_border_vel(int, int *, double *, int, int *) = 0;
   virtual int pack_border_hybrid(int, int *, double *) {return 0;}
   virtual void unpack_border(int, int, double *) = 0;
   virtual void unpack_border_vel(int, int, double *) = 0;
   virtual int unpack_border_hybrid(int, int, double *) {return 0;}
 
   virtual int pack_exchange(int, double *) = 0;
   virtual int unpack_exchange(double *) = 0;
 
   virtual int size_restart() = 0;
   virtual int pack_restart(int, double *) = 0;
   virtual int unpack_restart(double *) = 0;
 
   virtual void create_atom(int, double *) = 0;
 
   virtual void data_atom(double *, imageint, char **) = 0;
   virtual void data_atom_bonus(int, char **) {}
   virtual int data_atom_hybrid(int, char **) {return 0;}
   virtual void data_vel(int, char **);
   virtual int data_vel_hybrid(int, char **) {return 0;}
 
   virtual void pack_data(double **) = 0;
   virtual int pack_data_hybrid(int, double *) {return 0;}
   virtual void write_data(FILE *, int, double **) = 0;
   virtual int write_data_hybrid(FILE *, double *) {return 0;}
   virtual void pack_vel(double **);
   virtual int pack_vel_hybrid(int, double *) {return 0;}
   virtual void write_vel(FILE *, int, double **);
   virtual int write_vel_hybrid(FILE *, double *) {return 0;}
 
   int pack_bond(tagint **);
   void write_bond(FILE *, int, tagint **, int);
   int pack_angle(tagint **);
   void write_angle(FILE *, int, tagint **, int);
   void pack_dihedral(tagint **);
   void write_dihedral(FILE *, int, tagint **, int);
   void pack_improper(tagint **);
   void write_improper(FILE *, int, tagint **, int);
 
   virtual bigint memory_usage() = 0;
 
  protected:
   int nmax;                             // local copy of atom->nmax
   int deform_vremap;                    // local copy of domain properties
   int deform_groupbit;
   double *h_rate;
 
   // union data struct for packing 32-bit and 64-bit ints into double bufs
   // this avoids aliasing issues by having 2 pointers (double,int)
   //   to same buf memory
   // constructor for 32-bit int prevents compiler
   //   from possibly calling the double constructor when passed an int
   // copy to a double *buf:
   //   buf[m++] = ubuf(foo).d, where foo is a 32-bit or 64-bit int
   // copy from a double *buf:
   //   foo = (int) ubuf(buf[m++]).i;, where (int) or (tagint) match foo
   //   the cast prevents compiler warnings about possible truncation
 
   union ubuf {
     double d;
     int64_t i;
     ubuf(double arg) : d(arg) {}
     ubuf(int64_t arg) : i(arg) {}
     ubuf(int arg) : i(arg) {}
   };
 
   void grow_nmax();
   int grow_nmax_bonus(int);
 };
 
 }
 
 #endif
 
 /* ERROR/WARNING messages:
 
 E: Invalid atom_style command
 
 Self-explanatory.
 
 E: USER-CUDA package requires a cuda enabled atom_style
 
 Self-explanatory.
 
 */
diff --git a/src/min.cpp b/src/min.cpp
index 7093af4e4..116173b2d 100644
--- a/src/min.cpp
+++ b/src/min.cpp
@@ -1,788 +1,778 @@
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
 
    Copyright (2003) Sandia Corporation.  Under the terms of Contract
    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
    certain rights in this software.  This software is distributed under
    the GNU General Public License.
 
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
 /* ----------------------------------------------------------------------
    Contributing author: Aidan Thompson (SNL)
                         improved CG and backtrack ls, added quadratic ls
    Sources: Numerical Recipes frprmn routine
             "Conjugate Gradient Method Without the Agonizing Pain" by
             JR Shewchuk, http://www-2.cs.cmu.edu/~jrs/jrspapers.html#cg
 ------------------------------------------------------------------------- */
 
 #include "math.h"
 #include "stdlib.h"
 #include "string.h"
 #include "min.h"
 #include "atom.h"
+#include "atom_vec.h"
 #include "domain.h"
 #include "comm.h"
 #include "update.h"
 #include "modify.h"
 #include "fix_minimize.h"
 #include "compute.h"
 #include "neighbor.h"
 #include "force.h"
 #include "pair.h"
 #include "bond.h"
 #include "angle.h"
 #include "dihedral.h"
 #include "improper.h"
 #include "kspace.h"
 #include "output.h"
 #include "thermo.h"
 #include "timer.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
 Min::Min(LAMMPS *lmp) : Pointers(lmp)
 {
   dmax = 0.1;
   searchflag = 0;
   linestyle = 0;
 
   elist_global = elist_atom = NULL;
   vlist_global = vlist_atom = NULL;
 
   nextra_global = 0;
   fextra = NULL;
 
   nextra_atom = 0;
   xextra_atom = fextra_atom = NULL;
   extra_peratom = extra_nlen = NULL;
   extra_max = NULL;
   requestor = NULL;
 
   external_force_clear = 0;
 }
 
 /* ---------------------------------------------------------------------- */
 
 Min::~Min()
 {
   delete [] elist_global;
   delete [] elist_atom;
   delete [] vlist_global;
   delete [] vlist_atom;
 
   delete [] fextra;
 
   memory->sfree(xextra_atom);
   memory->sfree(fextra_atom);
   memory->destroy(extra_peratom);
   memory->destroy(extra_nlen);
   memory->destroy(extra_max);
   memory->sfree(requestor);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Min::init()
 {
   // create fix needed for storing atom-based quantities
   // will delete it at end of run
 
   char **fixarg = new char*[3];
   fixarg[0] = (char *) "MINIMIZE";
   fixarg[1] = (char *) "all";
   fixarg[2] = (char *) "MINIMIZE";
   modify->add_fix(3,fixarg);
   delete [] fixarg;
   fix_minimize = (FixMinimize *) modify->fix[modify->nfix-1];
 
   // clear out extra global and per-atom dof
   // will receive requests for new per-atom dof during pair init()
   // can then add vectors to fix_minimize in setup()
 
   nextra_global = 0;
   delete [] fextra;
   fextra = NULL;
 
   nextra_atom = 0;
   memory->sfree(xextra_atom);
   memory->sfree(fextra_atom);
   memory->destroy(extra_peratom);
   memory->destroy(extra_nlen);
   memory->destroy(extra_max);
   memory->sfree(requestor);
   xextra_atom = fextra_atom = NULL;
   extra_peratom = extra_nlen = NULL;
   extra_max = NULL;
   requestor = NULL;
 
   // virial_style:
   // 1 if computed explicitly by pair->compute via sum over pair interactions
   // 2 if computed implicitly by pair->virial_compute via sum over ghost atoms
 
   if (force->newton_pair) virial_style = 2;
   else virial_style = 1;
 
   // setup lists of computes for global and per-atom PE and pressure
 
   ev_setup();
 
   // detect if fix omp is present for clearing force arrays
 
   int ifix = modify->find_fix("package_omp");
   if (ifix >= 0) external_force_clear = 1;
 
-  // set flags for what arrays to clear in force_clear()
-  // need to clear additionals arrays if they exist
+  // set flags for arrays to clear in force_clear()
 
-  torqueflag = 0;
+  torqueflag = extraflag = 0;
   if (atom->torque_flag) torqueflag = 1;
-  erforceflag = 0;
-  if (atom->erforce_flag) erforceflag = 1;
-  e_flag = 0;
-  if (atom->e_flag) e_flag = 1;
-  rho_flag = 0;
-  if (atom->rho_flag) rho_flag = 1;
+  if (atom->avec->forceclearflag) extraflag = 1;
 
   // allow pair and Kspace compute() to be turned off via modify flags
 
   if (force->pair && force->pair->compute_flag) pair_compute_flag = 1;
   else pair_compute_flag = 0;
   if (force->kspace && force->kspace->compute_flag) kspace_compute_flag = 1;
   else kspace_compute_flag = 0;
 
   // orthogonal vs triclinic simulation box
 
   triclinic = domain->triclinic;
 
   // reset reneighboring criteria if necessary
 
   neigh_every = neighbor->every;
   neigh_delay = neighbor->delay;
   neigh_dist_check = neighbor->dist_check;
 
   if (neigh_every != 1 || neigh_delay != 0 || neigh_dist_check != 1) {
     if (comm->me == 0)
       error->warning(FLERR,
                      "Resetting reneighboring criteria during minimization");
   }
 
   neighbor->every = 1;
   neighbor->delay = 0;
   neighbor->dist_check = 1;
 
   niter = neval = 0;
 }
 
 /* ----------------------------------------------------------------------
    setup before run
 ------------------------------------------------------------------------- */
 
 void Min::setup()
 {
   if (comm->me == 0 && screen) fprintf(screen,"Setting up minimization ...\n");
 
   update->setupflag = 1;
 
   // setup extra global dof due to fixes
   // cannot be done in init() b/c update init() is before modify init()
 
   nextra_global = modify->min_dof();
   if (nextra_global) fextra = new double[nextra_global];
 
   // compute for potential energy
 
   int id = modify->find_compute("thermo_pe");
   if (id < 0) error->all(FLERR,"Minimization could not find thermo_pe compute");
   pe_compute = modify->compute[id];
 
   // style-specific setup does two tasks
   // setup extra global dof vectors
   // setup extra per-atom dof vectors due to requests from Pair classes
   // cannot be done in init() b/c update init() is before modify/pair init()
 
   setup_style();
 
   // ndoftotal = total dof for entire minimization problem
   // dof for atoms, extra per-atom, extra global
 
   bigint ndofme = 3*atom->nlocal;
   for (int m = 0; m < nextra_atom; m++)
     ndofme += extra_peratom[m]*atom->nlocal;
   MPI_Allreduce(&ndofme,&ndoftotal,1,MPI_LMP_BIGINT,MPI_SUM,world);
   ndoftotal += nextra_global;
 
   // setup domain, communication and neighboring
   // acquire ghosts
   // build neighbor lists
 
   atom->setup();
   modify->setup_pre_exchange();
   if (triclinic) domain->x2lamda(atom->nlocal);
   domain->pbc();
   domain->reset_box();
   comm->setup();
   if (neighbor->style) neighbor->setup_bins();
   comm->exchange();
   if (atom->sortfreq > 0) atom->sort();
   comm->borders();
   if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
   domain->image_check();
   domain->box_too_small_check();
   modify->setup_pre_neighbor();
   neighbor->build();
   neighbor->ncalls = 0;
 
   // remove these restriction eventually
 
   if (nextra_global && searchflag == 0)
     error->all(FLERR,
                "Cannot use a damped dynamics min style with fix box/relax");
   if (nextra_atom && searchflag == 0)
     error->all(FLERR,
                "Cannot use a damped dynamics min style with per-atom DOF");
 
   // atoms may have migrated in comm->exchange()
 
   reset_vectors();
 
   // compute all forces
 
   ev_set(update->ntimestep);
   force_clear();
   modify->setup_pre_force(vflag);
 
   if (pair_compute_flag) force->pair->compute(eflag,vflag);
   else if (force->pair) force->pair->compute_dummy(eflag,vflag);
 
   if (atom->molecular) {
     if (force->bond) force->bond->compute(eflag,vflag);
     if (force->angle) force->angle->compute(eflag,vflag);
     if (force->dihedral) force->dihedral->compute(eflag,vflag);
     if (force->improper) force->improper->compute(eflag,vflag);
   }
 
   if (force->kspace) {
     force->kspace->setup();
     if (kspace_compute_flag) force->kspace->compute(eflag,vflag);
     else force->kspace->compute_dummy(eflag,vflag);
   }
 
   if (force->newton) comm->reverse_comm();
 
   // update per-atom minimization variables stored by pair styles
 
   if (nextra_atom)
     for (int m = 0; m < nextra_atom; m++)
       requestor[m]->min_xf_get(m);
 
   modify->setup(vflag);
   output->setup();
   update->setupflag = 0;
 
   // stats for Finish to print
 
   ecurrent = pe_compute->compute_scalar();
   if (nextra_global) ecurrent += modify->min_energy(fextra);
   if (output->thermo->normflag) ecurrent /= atom->natoms;
 
   einitial = ecurrent;
   fnorm2_init = sqrt(fnorm_sqr());
   fnorminf_init = fnorm_inf();
 }
 
 /* ----------------------------------------------------------------------
    setup without output or one-time post-init setup
    flag = 0 = just force calculation
    flag = 1 = reneighbor and force calculation
 ------------------------------------------------------------------------- */
 
 void Min::setup_minimal(int flag)
 {
   update->setupflag = 1;
 
   // setup domain, communication and neighboring
   // acquire ghosts
   // build neighbor lists
 
   if (flag) {
     modify->setup_pre_exchange();
     if (triclinic) domain->x2lamda(atom->nlocal);
     domain->pbc();
     domain->reset_box();
     comm->setup();
     if (neighbor->style) neighbor->setup_bins();
     comm->exchange();
     comm->borders();
     if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
     domain->image_check();
     domain->box_too_small_check();
     modify->setup_pre_neighbor();
     neighbor->build();
     neighbor->ncalls = 0;
   }
 
   // atoms may have migrated in comm->exchange()
 
   reset_vectors();
 
   // compute all forces
 
   ev_set(update->ntimestep);
   force_clear();
   modify->setup_pre_force(vflag);
 
   if (pair_compute_flag) force->pair->compute(eflag,vflag);
   else if (force->pair) force->pair->compute_dummy(eflag,vflag);
 
   if (atom->molecular) {
     if (force->bond) force->bond->compute(eflag,vflag);
     if (force->angle) force->angle->compute(eflag,vflag);
     if (force->dihedral) force->dihedral->compute(eflag,vflag);
     if (force->improper) force->improper->compute(eflag,vflag);
   }
 
   if (force->kspace) {
     force->kspace->setup();
     if (kspace_compute_flag) force->kspace->compute(eflag,vflag);
     else force->kspace->compute_dummy(eflag,vflag);
   }
 
   if (force->newton) comm->reverse_comm();
 
   // update per-atom minimization variables stored by pair styles
 
   if (nextra_atom)
     for (int m = 0; m < nextra_atom; m++)
       requestor[m]->min_xf_get(m);
 
   modify->setup(vflag);
   update->setupflag = 0;
 
   // stats for Finish to print
 
   ecurrent = pe_compute->compute_scalar();
   if (nextra_global) ecurrent += modify->min_energy(fextra);
   if (output->thermo->normflag) ecurrent /= atom->natoms;
 
   einitial = ecurrent;
   fnorm2_init = sqrt(fnorm_sqr());
   fnorminf_init = fnorm_inf();
 }
 
 /* ----------------------------------------------------------------------
    perform minimization, calling iterate() for N steps
 ------------------------------------------------------------------------- */
 
 void Min::run(int n)
 {
   // minimizer iterations
 
   stop_condition = iterate(n);
   stopstr = stopstrings(stop_condition);
 
   // if early exit from iterate loop:
   // set update->nsteps to niter for Finish stats to print
   // set output->next values to this timestep
   // call energy_force() to insure vflag is set when forces computed
   // output->write does final output for thermo, dump, restart files
   // add ntimestep to all computes that store invocation times
   //   since are hardwiring call to thermo/dumps and computes may not be ready
 
   if (stop_condition) {
     update->nsteps = niter;
 
     if (update->restrict_output == 0) {
       for (int idump = 0; idump < output->ndump; idump++)
         output->next_dump[idump] = update->ntimestep;
       output->next_dump_any = update->ntimestep;
       if (output->restart_flag) {
         output->next_restart = update->ntimestep;
         if (output->restart_every_single)
           output->next_restart_single = update->ntimestep;
         if (output->restart_every_double)
           output->next_restart_double = update->ntimestep;
       }
     }
     output->next_thermo = update->ntimestep;
 
     modify->addstep_compute_all(update->ntimestep);
     ecurrent = energy_force(0);
     output->write(update->ntimestep);
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Min::cleanup()
 {
   // stats for Finish to print
 
   efinal = ecurrent;
   fnorm2_final = sqrt(fnorm_sqr());
   fnorminf_final = fnorm_inf();
 
   // reset reneighboring criteria
 
   neighbor->every = neigh_every;
   neighbor->delay = neigh_delay;
   neighbor->dist_check = neigh_dist_check;
 
   // delete fix at end of run, so its atom arrays won't persist
 
   modify->delete_fix("MINIMIZE");
   domain->box_too_small_check();
 }
 
 /* ----------------------------------------------------------------------
    evaluate potential energy and forces
    may migrate atoms due to reneighboring
    return new energy, which should include nextra_global dof
    return negative gradient stored in atom->f
    return negative gradient for nextra_global dof in fextra
 ------------------------------------------------------------------------- */
 
 double Min::energy_force(int resetflag)
 {
   // check for reneighboring
   // always communicate since minimizer moved atoms
 
   int nflag = neighbor->decide();
 
   if (nflag == 0) {
     timer->stamp();
     comm->forward_comm();
     timer->stamp(TIME_COMM);
   } else {
     if (modify->n_min_pre_exchange) modify->min_pre_exchange();
     if (triclinic) domain->x2lamda(atom->nlocal);
     domain->pbc();
     if (domain->box_change) {
       domain->reset_box();
       comm->setup();
       if (neighbor->style) neighbor->setup_bins();
     }
     timer->stamp();
     comm->exchange();
     if (atom->sortfreq > 0 &&
         update->ntimestep >= atom->nextsort) atom->sort();
     comm->borders();
     if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
     timer->stamp(TIME_COMM);
     neighbor->build();
     timer->stamp(TIME_NEIGHBOR);
   }
 
   ev_set(update->ntimestep);
   force_clear();
   if (modify->n_min_pre_force) modify->min_pre_force(vflag);
 
   timer->stamp();
 
   if (pair_compute_flag) {
     force->pair->compute(eflag,vflag);
     timer->stamp(TIME_PAIR);
   }
 
   if (atom->molecular) {
     if (force->bond) force->bond->compute(eflag,vflag);
     if (force->angle) force->angle->compute(eflag,vflag);
     if (force->dihedral) force->dihedral->compute(eflag,vflag);
     if (force->improper) force->improper->compute(eflag,vflag);
     timer->stamp(TIME_BOND);
   }
 
   if (kspace_compute_flag) {
     force->kspace->compute(eflag,vflag);
     timer->stamp(TIME_KSPACE);
   }
 
   if (force->newton) {
     comm->reverse_comm();
     timer->stamp(TIME_COMM);
   }
 
   // update per-atom minimization variables stored by pair styles
 
   if (nextra_atom)
     for (int m = 0; m < nextra_atom; m++)
       requestor[m]->min_xf_get(m);
 
   // fixes that affect minimization
 
   if (modify->n_min_post_force) modify->min_post_force(vflag);
 
   // compute potential energy of system
   // normalize if thermo PE does
 
   double energy = pe_compute->compute_scalar();
   if (nextra_global) energy += modify->min_energy(fextra);
   if (output->thermo->normflag) energy /= atom->natoms;
 
   // if reneighbored, atoms migrated
   // if resetflag = 1, update x0 of atoms crossing PBC
   // reset vectors used by lo-level minimizer
 
   if (nflag) {
     if (resetflag) fix_minimize->reset_coords();
     reset_vectors();
   }
 
   return energy;
 }
 
 /* ----------------------------------------------------------------------
    clear force on own & ghost atoms
-   setup and clear other arrays as needed
+   clear other arrays as needed
 ------------------------------------------------------------------------- */
 
 void Min::force_clear()
 {
   if (external_force_clear) return;
 
   // clear global force array
-  // nall includes ghosts only if either newton flag is set
+  // if either newton flag is set, also include ghosts
 
-  int nall;
-  if (force->newton) nall = atom->nlocal + atom->nghost;
-  else nall = atom->nlocal;
-
-  size_t nbytes = sizeof(double) * nall;
+  size_t nbytes = sizeof(double) * atom->nlocal;
+  if (force->newton) nbytes += sizeof(double) * atom->nghost;
 
   if (nbytes) {
-    memset(&(atom->f[0][0]),0,3*nbytes);
-    if (torqueflag)  memset(&(atom->torque[0][0]),0,3*nbytes);
-    if (erforceflag) memset(&(atom->erforce[0]),  0,  nbytes);
-    if (e_flag)      memset(&(atom->de[0]),       0,  nbytes);
-    if (rho_flag)    memset(&(atom->drho[0]),     0,  nbytes);
+    memset(&atom->f[0][0],0,3*nbytes);
+    if (torqueflag) memset(&atom->torque[0][0],0,3*nbytes);
+    if (extraflag) atom->avec->force_clear(0,nbytes);
   }
 }
 
 /* ----------------------------------------------------------------------
    pair style makes request to add a per-atom variables to minimization
    requestor stores callback to pair class to invoke during min
      to get current variable and forces on it and to update the variable
    return flag that pair can use if it registers multiple variables
 ------------------------------------------------------------------------- */
 
 int Min::request(Pair *pair, int peratom, double maxvalue)
 {
   int n = nextra_atom + 1;
   xextra_atom = (double **) memory->srealloc(xextra_atom,n*sizeof(double *),
                                              "min:xextra_atom");
   fextra_atom = (double **) memory->srealloc(fextra_atom,n*sizeof(double *),
                                              "min:fextra_atom");
   memory->grow(extra_peratom,n,"min:extra_peratom");
   memory->grow(extra_nlen,n,"min:extra_nlen");
   memory->grow(extra_max,n,"min:extra_max");
   requestor = (Pair **) memory->srealloc(requestor,n*sizeof(Pair *),
                                          "min:requestor");
 
   requestor[nextra_atom] = pair;
   extra_peratom[nextra_atom] = peratom;
   extra_max[nextra_atom] = maxvalue;
   nextra_atom++;
   return nextra_atom-1;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Min::modify_params(int narg, char **arg)
 {
   if (narg == 0) error->all(FLERR,"Illegal min_modify command");
 
   int iarg = 0;
   while (iarg < narg) {
     if (strcmp(arg[iarg],"dmax") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command");
       dmax = force->numeric(FLERR,arg[iarg+1]);
       iarg += 2;
     } else if (strcmp(arg[iarg],"line") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command");
       if (strcmp(arg[iarg+1],"backtrack") == 0) linestyle = 0;
       else if (strcmp(arg[iarg+1],"quadratic") == 0) linestyle = 1;
       else if (strcmp(arg[iarg+1],"forcezero") == 0) linestyle = 2;
       else error->all(FLERR,"Illegal min_modify command");
       iarg += 2;
     } else error->all(FLERR,"Illegal min_modify command");
   }
 }
 
 /* ----------------------------------------------------------------------
    setup lists of computes for global and per-atom PE and pressure
 ------------------------------------------------------------------------- */
 
 void Min::ev_setup()
 {
   delete [] elist_global;
   delete [] elist_atom;
   delete [] vlist_global;
   delete [] vlist_atom;
   elist_global = elist_atom = NULL;
   vlist_global = vlist_atom = NULL;
 
   nelist_global = nelist_atom = 0;
   nvlist_global = nvlist_atom = 0;
   for (int i = 0; i < modify->ncompute; i++) {
     if (modify->compute[i]->peflag) nelist_global++;
     if (modify->compute[i]->peatomflag) nelist_atom++;
     if (modify->compute[i]->pressflag) nvlist_global++;
     if (modify->compute[i]->pressatomflag) nvlist_atom++;
   }
 
   if (nelist_global) elist_global = new Compute*[nelist_global];
   if (nelist_atom) elist_atom = new Compute*[nelist_atom];
   if (nvlist_global) vlist_global = new Compute*[nvlist_global];
   if (nvlist_atom) vlist_atom = new Compute*[nvlist_atom];
 
   nelist_global = nelist_atom = 0;
   nvlist_global = nvlist_atom = 0;
   for (int i = 0; i < modify->ncompute; i++) {
     if (modify->compute[i]->peflag)
       elist_global[nelist_global++] = modify->compute[i];
     if (modify->compute[i]->peatomflag)
       elist_atom[nelist_atom++] = modify->compute[i];
     if (modify->compute[i]->pressflag)
       vlist_global[nvlist_global++] = modify->compute[i];
     if (modify->compute[i]->pressatomflag)
       vlist_atom[nvlist_atom++] = modify->compute[i];
   }
 }
 
 /* ----------------------------------------------------------------------
    set eflag,vflag for current iteration
    invoke matchstep() on all timestep-dependent computes to clear their arrays
    eflag/vflag based on computes that need info on this ntimestep
    always set eflag_global = 1, since need energy every iteration
    eflag = 0 = no energy computation
    eflag = 1 = global energy only
    eflag = 2 = per-atom energy only
    eflag = 3 = both global and per-atom energy
    vflag = 0 = no virial computation (pressure)
    vflag = 1 = global virial with pair portion via sum of pairwise interactions
    vflag = 2 = global virial with pair portion via F dot r including ghosts
    vflag = 4 = per-atom virial only
    vflag = 5 or 6 = both global and per-atom virial
 ------------------------------------------------------------------------- */
 
 void Min::ev_set(bigint ntimestep)
 {
   int i,flag;
 
   int eflag_global = 1;
   for (i = 0; i < nelist_global; i++)
     elist_global[i]->matchstep(ntimestep);
 
   flag = 0;
   int eflag_atom = 0;
   for (i = 0; i < nelist_atom; i++)
     if (elist_atom[i]->matchstep(ntimestep)) flag = 1;
   if (flag) eflag_atom = 2;
 
   if (eflag_global) update->eflag_global = update->ntimestep;
   if (eflag_atom) update->eflag_atom = update->ntimestep;
   eflag = eflag_global + eflag_atom;
 
   flag = 0;
   int vflag_global = 0;
   for (i = 0; i < nvlist_global; i++)
     if (vlist_global[i]->matchstep(ntimestep)) flag = 1;
   if (flag) vflag_global = virial_style;
 
   flag = 0;
   int vflag_atom = 0;
   for (i = 0; i < nvlist_atom; i++)
     if (vlist_atom[i]->matchstep(ntimestep)) flag = 1;
   if (flag) vflag_atom = 4;
 
   if (vflag_global) update->vflag_global = update->ntimestep;
   if (vflag_atom) update->vflag_atom = update->ntimestep;
   vflag = vflag_global + vflag_atom;
 }
 
 /* ----------------------------------------------------------------------
    compute and return ||force||_2^2
 ------------------------------------------------------------------------- */
 
 double Min::fnorm_sqr()
 {
   int i,n;
   double *fatom;
 
   double local_norm2_sqr = 0.0;
   for (i = 0; i < nvec; i++) local_norm2_sqr += fvec[i]*fvec[i];
   if (nextra_atom) {
     for (int m = 0; m < nextra_atom; m++) {
       fatom = fextra_atom[m];
       n = extra_nlen[m];
       for (i = 0; i < n; i++)
         local_norm2_sqr += fatom[i]*fatom[i];
     }
   }
 
   double norm2_sqr = 0.0;
   MPI_Allreduce(&local_norm2_sqr,&norm2_sqr,1,MPI_DOUBLE,MPI_SUM,world);
 
   if (nextra_global)
     for (i = 0; i < nextra_global; i++)
       norm2_sqr += fextra[i]*fextra[i];
 
   return norm2_sqr;
 }
 
 /* ----------------------------------------------------------------------
    compute and return ||force||_inf
 ------------------------------------------------------------------------- */
 
 double Min::fnorm_inf()
 {
   int i,n;
   double *fatom;
 
   double local_norm_inf = 0.0;
   for (i = 0; i < nvec; i++)
     local_norm_inf = MAX(fabs(fvec[i]),local_norm_inf);
   if (nextra_atom) {
     for (int m = 0; m < nextra_atom; m++) {
       fatom = fextra_atom[m];
       n = extra_nlen[m];
       for (i = 0; i < n; i++)
         local_norm_inf = MAX(fabs(fatom[i]),local_norm_inf);
     }
   }
 
   double norm_inf = 0.0;
   MPI_Allreduce(&local_norm_inf,&norm_inf,1,MPI_DOUBLE,MPI_MAX,world);
 
   if (nextra_global)
     for (i = 0; i < nextra_global; i++)
       norm_inf = MAX(fabs(fextra[i]),norm_inf);
 
   return norm_inf;
 }
 
 /* ----------------------------------------------------------------------
    possible stop conditions
 ------------------------------------------------------------------------- */
 
 char *Min::stopstrings(int n)
 {
   const char *strings[] = {"max iterations",
                            "max force evaluations",
                            "energy tolerance",
                            "force tolerance",
                            "search direction is not downhill",
                            "linesearch alpha is zero",
                            "forces are zero",
                            "quadratic factors are zero",
                            "trust region too small",
                            "HFTN minimizer error"};
   return (char *) strings[n];
 }
diff --git a/src/min.h b/src/min.h
index 2ba3e5cb2..fc78fd873 100644
--- a/src/min.h
+++ b/src/min.h
@@ -1,144 +1,143 @@
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
 
    Copyright (2003) Sandia Corporation.  Under the terms of Contract
    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
    certain rights in this software.  This software is distributed under
    the GNU General Public License.
 
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
 #ifndef LMP_MIN_H
 #define LMP_MIN_H
 
 #include "pointers.h"
 
 namespace LAMMPS_NS {
 
 class Min : protected Pointers {
  public:
   double einitial,efinal,eprevious;
   double fnorm2_init,fnorminf_init,fnorm2_final,fnorminf_final;
   double alpha_final;
   int niter,neval;
   int stop_condition;
   char *stopstr;
   int searchflag;     // 0 if damped dynamics, 1 if sub-cycles on local search
 
   Min(class LAMMPS *);
   virtual ~Min();
   virtual void init();
   void setup();
   void setup_minimal(int);
   void run(int);
   void cleanup();
   int request(class Pair *, int, double);
   virtual bigint memory_usage() {return 0;}
   void modify_params(int, char **);
   double fnorm_sqr();
   double fnorm_inf();
 
   virtual void init_style() {}
   virtual void setup_style() = 0;
   virtual void reset_vectors() = 0;
   virtual int iterate(int) = 0;
 
  protected:
   int eflag,vflag;            // flags for energy/virial computation
   int virial_style;           // compute virial explicitly or implicitly
   int external_force_clear;   // clear forces locally or externally
 
   double dmax;                // max dist to move any atom in one step
   int linestyle;              // 0 = backtrack, 1 = quadratic, 2 = forcezero
 
   int nelist_global,nelist_atom;    // # of PE,virial computes to check
   int nvlist_global,nvlist_atom;
   class Compute **elist_global;     // lists of PE,virial Computes
   class Compute **elist_atom;
   class Compute **vlist_global;
   class Compute **vlist_atom;
 
   int triclinic;              // 0 if domain is orthog, 1 if triclinic
   int pairflag;
-  int torqueflag,erforceflag;
-  int e_flag,rho_flag;
+  int torqueflag,extraflag;
 
   int pair_compute_flag;            // 0 if pair->compute is skipped
   int kspace_compute_flag;          // 0 if kspace->compute is skipped
 
   int narray;                       // # of arrays stored by fix_minimize
   class FixMinimize *fix_minimize;  // fix that stores auxiliary data
 
   class Compute *pe_compute;        // compute for potential energy
   double ecurrent;                  // current potential energy
 
   bigint ndoftotal;           // total dof for entire problem
 
   int nvec;                   // local atomic dof = length of xvec
   double *xvec;               // variables for atomic dof, as 1d vector
   double *fvec;               // force vector for atomic dof, as 1d vector
 
   int nextra_global;          // # of extra global dof due to fixes
   double *fextra;             // force vector for extra global dof
                               // xextra is stored by fix
 
   int nextra_atom;            // # of extra per-atom variables
   double **xextra_atom;       // ptr to the variable
   double **fextra_atom;       // ptr to the force on the variable
   int *extra_peratom;         // # of values in variable, e.g. 3 in x
   int *extra_nlen;            // total local length of variable, e.g 3*nlocal
   double *extra_max;          // max allowed change per iter for atom's var
   class Pair **requestor;     // Pair that stores/manipulates the variable
 
   int neigh_every,neigh_delay,neigh_dist_check;  // neighboring params
 
   double energy_force(int);
   void force_clear();
 
   double compute_force_norm_sqr();
   double compute_force_norm_inf();
 
   void ev_setup();
   void ev_set(bigint);
 
   char *stopstrings(int);
 };
 
 }
 
 #endif
 
 /* ERROR/WARNING messages:
 
 W: Resetting reneighboring criteria during minimization
 
 Minimization requires that neigh_modify settings be delay = 0, every =
 1, check = yes.  Since these settings were not in place, LAMMPS
 changed them and will restore them to their original values after the
 minimization.
 
 E: Minimization could not find thermo_pe compute
 
 This compute is created by the thermo command.  It must have been
 explicitly deleted by a uncompute command.
 
 E: Cannot use a damped dynamics min style with fix box/relax
 
 This is a current restriction in LAMMPS.  Use another minimizer
 style.
 
 E: Cannot use a damped dynamics min style with per-atom DOF
 
 This is a current restriction in LAMMPS.  Use another minimizer
 style.
 
 E: Illegal ... command
 
 Self-explanatory.  Check the input script syntax and compare to the
 documentation for the command.  You can use -echo screen as a
 command-line option when running LAMMPS to see the offending line.
 
 */
diff --git a/src/respa.cpp b/src/respa.cpp
index 720527d6f..b648b8bfc 100644
--- a/src/respa.cpp
+++ b/src/respa.cpp
@@ -1,706 +1,697 @@
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
 
    Copyright (2003) Sandia Corporation.  Under the terms of Contract
    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
    certain rights in this software.  This software is distributed under
    the GNU General Public License.
 
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
 /* ----------------------------------------------------------------------
    Contributing authors: Mark Stevens (SNL), Paul Crozier (SNL)
 ------------------------------------------------------------------------- */
 
 #include "stdlib.h"
 #include "string.h"
 #include "respa.h"
 #include "neighbor.h"
 #include "atom.h"
+#include "atom_vec.h"
 #include "domain.h"
 #include "comm.h"
 #include "force.h"
 #include "pair.h"
 #include "bond.h"
 #include "angle.h"
 #include "dihedral.h"
 #include "improper.h"
 #include "kspace.h"
 #include "output.h"
 #include "update.h"
 #include "modify.h"
 #include "compute.h"
 #include "fix_respa.h"
 #include "timer.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
 Respa::Respa(LAMMPS *lmp, int narg, char **arg) : Integrate(lmp, narg, arg)
 {
   if (narg < 1) error->all(FLERR,"Illegal run_style respa command");
 
   nlevels = force->inumeric(FLERR,arg[0]);
   if (nlevels < 1) error->all(FLERR,"Respa levels must be >= 1");
 
   if (narg < nlevels) error->all(FLERR,"Illegal run_style respa command");
   loop = new int[nlevels];
   for (int iarg = 1; iarg < nlevels; iarg++) {
     loop[iarg-1] = force->inumeric(FLERR,arg[iarg]);
     if (loop[iarg-1] <= 0) error->all(FLERR,"Illegal run_style respa command");
   }
   loop[nlevels-1] = 1;
 
   // set level at which each force is computed
   // argument settings override defaults
 
   level_bond = level_angle = level_dihedral = level_improper = -1;
   level_pair = level_kspace = -1;
   level_inner = level_middle = level_outer = -1;
 
   int iarg = nlevels;
   while (iarg < narg) {
     if (strcmp(arg[iarg],"bond") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal run_style respa command");
       level_bond = force->inumeric(FLERR,arg[iarg+1]) - 1;
       iarg += 2;
     } else if (strcmp(arg[iarg],"angle") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal run_style respa command");
       level_angle = force->inumeric(FLERR,arg[iarg+1]) - 1;
       iarg += 2;
     } else if (strcmp(arg[iarg],"dihedral") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal run_style respa command");
       level_dihedral = force->inumeric(FLERR,arg[iarg+1]) - 1;
       iarg += 2;
     } else if (strcmp(arg[iarg],"improper") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal run_style respa command");
       level_improper = force->inumeric(FLERR,arg[iarg+1]) - 1;
       iarg += 2;
     } else if (strcmp(arg[iarg],"pair") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal run_style respa command");
       level_pair = force->inumeric(FLERR,arg[iarg+1]) - 1;
       iarg += 2;
     } else if (strcmp(arg[iarg],"inner") == 0) {
       if (iarg+4 > narg) error->all(FLERR,"Illegal run_style respa command");
       level_inner = force->inumeric(FLERR,arg[iarg+1]) - 1;
       cutoff[0] = force->numeric(FLERR,arg[iarg+2]);
       cutoff[1] = force->numeric(FLERR,arg[iarg+3]);
       iarg += 4;
     } else if (strcmp(arg[iarg],"middle") == 0) {
       if (iarg+4 > narg) error->all(FLERR,"Illegal run_style respa command");
       level_middle = force->inumeric(FLERR,arg[iarg+1]) - 1;
       cutoff[2] = force->numeric(FLERR,arg[iarg+2]);
       cutoff[3] = force->numeric(FLERR,arg[iarg+3]);
       iarg += 4;
     } else if (strcmp(arg[iarg],"outer") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal run_style respa command");
       level_outer = force->inumeric(FLERR,arg[iarg+1]) - 1;
       iarg += 2;
     } else if (strcmp(arg[iarg],"kspace") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal run_style respa command");
       level_kspace = force->inumeric(FLERR,arg[iarg+1]) - 1;
       iarg += 2;
     } else error->all(FLERR,"Illegal run_style respa command");
   }
 
   // cannot specify both pair and inner/middle/outer
 
   if (level_pair >= 0 &&
       (level_inner >= 0 || level_middle >= 0 || level_outer >= 0))
     error->all(FLERR,"Cannot set both respa pair and inner/middle/outer");
 
   // if either inner and outer is specified, then both must be
 
   if ((level_inner >= 0 && level_outer == -1) ||
       (level_outer >= 0 && level_inner == -1))
     error->all(FLERR,"Must set both respa inner and outer");
 
   // middle cannot be set without inner/outer
 
   if (level_middle >= 0 && level_inner == -1)
     error->all(FLERR,"Cannot set respa middle without inner/outer");
 
   // set defaults if user did not specify level
   // bond to innermost level
   // angle same as bond, dihedral same as angle, improper same as dihedral
   // pair to outermost level if no inner/middle/outer
   // inner/middle/outer have no defaults
   // kspace same as pair or outer
 
   if (level_bond == -1) level_bond = 0;
   if (level_angle == -1) level_angle = level_bond;
   if (level_dihedral == -1) level_dihedral = level_angle;
   if (level_improper == -1) level_improper = level_dihedral;
   if (level_pair == -1 && level_inner == -1) level_pair = nlevels-1;
   if (level_kspace == -1 && level_pair >= 0) level_kspace = level_pair;
   if (level_kspace == -1 && level_pair == -1) level_kspace = level_outer;
 
   // print respa levels
 
   if (comm->me == 0) {
     if (screen) {
       fprintf(screen,"Respa levels:\n");
       for (int i = 0; i < nlevels; i++) {
         fprintf(screen,"  %d =",i+1);
         if (level_bond == i) fprintf(screen," bond");
         if (level_angle == i) fprintf(screen," angle");
         if (level_dihedral == i) fprintf(screen," dihedral");
         if (level_improper == i) fprintf(screen," improper");
         if (level_pair == i) fprintf(screen," pair");
         if (level_inner == i) fprintf(screen," pair-inner");
         if (level_middle == i) fprintf(screen," pair-middle");
         if (level_outer == i) fprintf(screen," pair-outer");
         if (level_kspace == i) fprintf(screen," kspace");
         fprintf(screen,"\n");
       }
     }
     if (logfile) {
       fprintf(logfile,"Respa levels:\n");
       for (int i = 0; i < nlevels; i++) {
         fprintf(logfile,"  %d =",i+1);
         if (level_bond == i) fprintf(logfile," bond");
         if (level_angle == i) fprintf(logfile," angle");
         if (level_dihedral == i) fprintf(logfile," dihedral");
         if (level_improper == i) fprintf(logfile," improper");
         if (level_pair == i) fprintf(logfile," pair");
         if (level_inner == i) fprintf(logfile," pair-inner");
         if (level_middle == i) fprintf(logfile," pair-middle");
         if (level_outer == i) fprintf(logfile," pair-outer");
         if (level_kspace == i) fprintf(logfile," kspace");
         fprintf(logfile,"\n");
       }
     }
   }
 
   // check that levels are in correct order
 
   if (level_angle < level_bond || level_dihedral < level_angle ||
       level_improper < level_dihedral)
     error->all(FLERR,"Invalid order of forces within respa levels");
   if (level_pair >= 0) {
     if (level_pair < level_improper || level_kspace < level_pair)
       error->all(FLERR,"Invalid order of forces within respa levels");
   }
   if (level_pair == -1 && level_middle == -1) {
     if (level_inner < level_improper || level_outer < level_inner ||
         level_kspace < level_outer)
       error->all(FLERR,"Invalid order of forces within respa levels");
   }
   if (level_pair == -1 && level_middle >= 0) {
     if (level_inner < level_improper || level_middle < level_inner ||
         level_outer < level_inner || level_kspace < level_outer)
       error->all(FLERR,"Invalid order of forces within respa levels");
   }
 
   // warn if any levels are devoid of forces
 
   int flag = 0;
   for (int i = 0; i < nlevels; i++)
     if (level_bond != i && level_angle != i && level_dihedral != i &&
         level_improper != i && level_pair != i && level_inner != i &&
         level_middle != i && level_outer != i && level_kspace != i) flag = 1;
   if (flag && comm->me == 0)
     error->warning(FLERR,"One or more respa levels compute no forces");
 
   // check cutoff consistency if inner/middle/outer are enabled
 
   if (level_inner >= 0 && cutoff[1] < cutoff[0])
     error->all(FLERR,"Respa inner cutoffs are invalid");
   if (level_middle >= 0 && (cutoff[3] < cutoff[2] || cutoff[2] < cutoff[1]))
     error->all(FLERR,"Respa middle cutoffs are invalid");
 
   // set outer pair of cutoffs to inner pair if middle is not enabled
 
   if (level_inner >= 0 && level_middle < 0) {
     cutoff[2] = cutoff[0];
     cutoff[3] = cutoff[1];
   }
 
   // allocate other needed arrays
 
   newton = new int[nlevels];
   step = new double[nlevels];
 }
 
 /* ---------------------------------------------------------------------- */
 
 Respa::~Respa()
 {
   delete [] loop;
   delete [] newton;
   delete [] step;
 }
 
 /* ----------------------------------------------------------------------
    initialization before run
 ------------------------------------------------------------------------- */
 
 void Respa::init()
 {
   Integrate::init();
 
   // warn if no fixes
 
   if (modify->nfix == 0 && comm->me == 0)
     error->warning(FLERR,"No fixes defined, atoms won't move");
 
   // create fix needed for storing atom-based respa level forces
   // will delete it at end of run
 
   char **fixarg = new char*[4];
   fixarg[0] = (char *) "RESPA";
   fixarg[1] = (char *) "all";
   fixarg[2] = (char *) "RESPA";
   fixarg[3] = new char[8];
   sprintf(fixarg[3],"%d",nlevels);
   modify->add_fix(4,fixarg);
   delete [] fixarg[3];
   delete [] fixarg;
   fix_respa = (FixRespa *) modify->fix[modify->nfix-1];
 
   // insure respa inner/middle/outer is using Pair class that supports it
 
   if (level_inner >= 0)
     if (force->pair && force->pair->respa_enable == 0)
       error->all(FLERR,"Pair style does not support rRESPA inner/middle/outer");
 
   // virial_style = 1 (explicit) since never computed implicitly like Verlet
 
   virial_style = 1;
 
   // setup lists of computes for global and per-atom PE and pressure
 
   ev_setup();
 
   // detect if fix omp is present and will clear force arrays
 
   int ifix = modify->find_fix("package_omp");
   if (ifix >= 0) external_force_clear = 1;
 
-  // set flags for what arrays to clear in force_clear()
-  // need to clear additionals arrays if they exist
+  // set flags for arrays to clear in force_clear()
 
-  torqueflag = 0;
+  torqueflag = extraflag = 0;
   if (atom->torque_flag) torqueflag = 1;
-  erforceflag = 0;
-  if (atom->erforce_flag) erforceflag = 1;
-  e_flag = 0;
-  if (atom->e_flag) e_flag = 1;
-  rho_flag = 0;
-  if (atom->rho_flag) rho_flag = 1;
+  if (atom->avec->forceclearflag) extraflag = 1;
 
   // step[] = timestep for each level
 
   step[nlevels-1] = update->dt;
   for (int ilevel = nlevels-2; ilevel >= 0; ilevel--)
     step[ilevel] = step[ilevel+1]/loop[ilevel];
 
   // set newton flag for each level
 
   for (int ilevel = 0; ilevel < nlevels; ilevel++) {
     newton[ilevel] = 0;
     if (force->newton_bond) {
       if (level_bond == ilevel || level_angle == ilevel ||
           level_dihedral == ilevel || level_improper == ilevel)
         newton[ilevel] = 1;
     }
     if (force->newton_pair) {
       if (level_pair == ilevel || level_inner == ilevel ||
           level_middle == ilevel || level_outer == ilevel)
         newton[ilevel] = 1;
     }
   }
 
   // orthogonal vs triclinic simulation box
 
   triclinic = domain->triclinic;
 }
 
 /* ----------------------------------------------------------------------
    setup before run
 ------------------------------------------------------------------------- */
 
 void Respa::setup()
 {
   if (comm->me == 0 && screen) fprintf(screen,"Setting up run ...\n");
 
   update->setupflag = 1;
 
   // setup domain, communication and neighboring
   // acquire ghosts
   // build neighbor lists
 
   atom->setup();
   modify->setup_pre_exchange();
   if (triclinic) domain->x2lamda(atom->nlocal);
   domain->pbc();
   domain->reset_box();
   comm->setup();
   if (neighbor->style) neighbor->setup_bins();
   comm->exchange();
   if (atom->sortfreq > 0) atom->sort();
   comm->borders();
   if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
   domain->image_check();
   domain->box_too_small_check();
   modify->setup_pre_neighbor();
   neighbor->build();
   neighbor->ncalls = 0;
 
   // compute all forces
 
   ev_set(update->ntimestep);
 
   for (int ilevel = 0; ilevel < nlevels; ilevel++) {
     force_clear(newton[ilevel]);
     modify->setup_pre_force_respa(vflag,ilevel);
     if (level_pair == ilevel && pair_compute_flag)
       force->pair->compute(eflag,vflag);
     if (level_inner == ilevel && pair_compute_flag)
       force->pair->compute_inner();
     if (level_middle == ilevel && pair_compute_flag)
       force->pair->compute_middle();
     if (level_outer == ilevel && pair_compute_flag)
       force->pair->compute_outer(eflag,vflag);
     if (level_bond == ilevel && force->bond)
       force->bond->compute(eflag,vflag);
     if (level_angle == ilevel && force->angle)
       force->angle->compute(eflag,vflag);
     if (level_dihedral == ilevel && force->dihedral)
       force->dihedral->compute(eflag,vflag);
     if (level_improper == ilevel && force->improper)
       force->improper->compute(eflag,vflag);
     if (level_kspace == ilevel && force->kspace) {
       force->kspace->setup();
       if (kspace_compute_flag) force->kspace->compute(eflag,vflag);
     }
     if (newton[ilevel]) comm->reverse_comm();
     copy_f_flevel(ilevel);
   }
 
   sum_flevel_f();
   modify->setup(vflag);
   output->setup();
   update->setupflag = 0;
 }
 
 /* ----------------------------------------------------------------------
    setup without output
    flag = 0 = just force calculation
    flag = 1 = reneighbor and force calculation
 ------------------------------------------------------------------------- */
 
 void Respa::setup_minimal(int flag)
 {
   update->setupflag = 1;
 
   // setup domain, communication and neighboring
   // acquire ghosts
   // build neighbor lists
 
   if (flag) {
     modify->setup_pre_exchange();
     if (triclinic) domain->x2lamda(atom->nlocal);
     domain->pbc();
     domain->reset_box();
     comm->setup();
     if (neighbor->style) neighbor->setup_bins();
     comm->exchange();
     comm->borders();
     if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
     domain->image_check();
     domain->box_too_small_check();
     modify->setup_pre_neighbor();
     neighbor->build();
     neighbor->ncalls = 0;
   }
 
   // compute all forces
 
   ev_set(update->ntimestep);
 
   for (int ilevel = 0; ilevel < nlevels; ilevel++) {
     force_clear(newton[ilevel]);
     modify->setup_pre_force_respa(vflag,ilevel);
     if (level_pair == ilevel && pair_compute_flag)
       force->pair->compute(eflag,vflag);
     if (level_inner == ilevel && pair_compute_flag)
       force->pair->compute_inner();
     if (level_middle == ilevel && pair_compute_flag)
       force->pair->compute_middle();
     if (level_outer == ilevel && pair_compute_flag)
       force->pair->compute_outer(eflag,vflag);
     if (level_bond == ilevel && force->bond)
       force->bond->compute(eflag,vflag);
     if (level_angle == ilevel && force->angle)
       force->angle->compute(eflag,vflag);
     if (level_dihedral == ilevel && force->dihedral)
       force->dihedral->compute(eflag,vflag);
     if (level_improper == ilevel && force->improper)
       force->improper->compute(eflag,vflag);
     if (level_kspace == ilevel && force->kspace) {
       force->kspace->setup();
       if (kspace_compute_flag) force->kspace->compute(eflag,vflag);
     }
     if (newton[ilevel]) comm->reverse_comm();
     copy_f_flevel(ilevel);
   }
 
   sum_flevel_f();
   modify->setup(vflag);
   update->setupflag = 0;
 }
 
 /* ----------------------------------------------------------------------
    run for N steps
 ------------------------------------------------------------------------- */
 
 void Respa::run(int n)
 {
   bigint ntimestep;
 
   for (int i = 0; i < n; i++) {
 
     ntimestep = ++update->ntimestep;
     ev_set(ntimestep);
 
     recurse(nlevels-1);
 
     // needed in case end_of_step() or output() use total force
 
     sum_flevel_f();
 
     if (modify->n_end_of_step) modify->end_of_step();
 
     if (ntimestep == output->next) {
       timer->stamp();
       output->write(update->ntimestep);
       timer->stamp(TIME_OUTPUT);
     }
   }
 }
 
 /* ----------------------------------------------------------------------
    delete rRESPA fix at end of run, so its atom arrays won't persist
 ------------------------------------------------------------------------- */
 
 void Respa::cleanup()
 {
   modify->post_run();
   modify->delete_fix("RESPA");
   domain->box_too_small_check();
   update->update_time();
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Respa::reset_dt()
 {
   step[nlevels-1] = update->dt;
   for (int ilevel = nlevels-2; ilevel >= 0; ilevel--)
     step[ilevel] = step[ilevel+1]/loop[ilevel];
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Respa::recurse(int ilevel)
 {
   copy_flevel_f(ilevel);
 
   for (int iloop = 0; iloop < loop[ilevel]; iloop++) {
 
     modify->initial_integrate_respa(vflag,ilevel,iloop);
     if (modify->n_post_integrate_respa)
       modify->post_integrate_respa(ilevel,iloop);
 
     // at outermost level, check on rebuilding neighbor list
     // at innermost level, communicate
     // at middle levels, do nothing
 
     if (ilevel == nlevels-1) {
       int nflag = neighbor->decide();
       if (nflag) {
         if (modify->n_pre_exchange) modify->pre_exchange();
         if (triclinic) domain->x2lamda(atom->nlocal);
         domain->pbc();
         if (domain->box_change) {
           domain->reset_box();
           comm->setup();
           if (neighbor->style) neighbor->setup_bins();
         }
         timer->stamp();
         comm->exchange();
         if (atom->sortfreq > 0 &&
             update->ntimestep >= atom->nextsort) atom->sort();
         comm->borders();
         if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
         timer->stamp(TIME_COMM);
         if (modify->n_pre_neighbor) modify->pre_neighbor();
         neighbor->build();
         timer->stamp(TIME_NEIGHBOR);
       } else if (ilevel == 0) {
         timer->stamp();
         comm->forward_comm();
         timer->stamp(TIME_COMM);
       }
 
     } else if (ilevel == 0) {
       timer->stamp();
       comm->forward_comm();
       timer->stamp(TIME_COMM);
     }
 
     // rRESPA recursion thru all levels
     // this used to be before neigh list build,
     // which prevented per-atom energy/stress being tallied correctly
     // b/c atoms migrated to new procs between short/long force calls
     // now they migrate at very start of rRESPA timestep, before all forces
 
     if (ilevel) recurse(ilevel-1);
 
     // force computations
     // important that ordering is same as Verlet
     // so that any order dependencies are the same
     // when potentials are invoked at same level
 
     force_clear(newton[ilevel]);
     if (modify->n_pre_force_respa)
       modify->pre_force_respa(vflag,ilevel,iloop);
 
     timer->stamp();
     if (level_pair == ilevel && pair_compute_flag) {
       force->pair->compute(eflag,vflag);
       timer->stamp(TIME_PAIR);
     }
     if (level_inner == ilevel && pair_compute_flag) {
       force->pair->compute_inner();
       timer->stamp(TIME_PAIR);
     }
     if (level_middle == ilevel && pair_compute_flag) {
       force->pair->compute_middle();
       timer->stamp(TIME_PAIR);
     }
     if (level_outer == ilevel && pair_compute_flag) {
       force->pair->compute_outer(eflag,vflag);
       timer->stamp(TIME_PAIR);
     }
     if (level_bond == ilevel && force->bond) {
       force->bond->compute(eflag,vflag);
       timer->stamp(TIME_BOND);
     }
     if (level_angle == ilevel && force->angle) {
       force->angle->compute(eflag,vflag);
       timer->stamp(TIME_BOND);
     }
     if (level_dihedral == ilevel && force->dihedral) {
       force->dihedral->compute(eflag,vflag);
       timer->stamp(TIME_BOND);
     }
     if (level_improper == ilevel && force->improper) {
       force->improper->compute(eflag,vflag);
       timer->stamp(TIME_BOND);
     }
     if (level_kspace == ilevel && kspace_compute_flag) {
       force->kspace->compute(eflag,vflag);
       timer->stamp(TIME_KSPACE);
     }
 
     if (newton[ilevel]) {
       comm->reverse_comm();
       timer->stamp(TIME_COMM);
     }
 
     if (modify->n_post_force_respa)
       modify->post_force_respa(vflag,ilevel,iloop);
     modify->final_integrate_respa(ilevel,iloop);
   }
 
   copy_f_flevel(ilevel);
 }
 
 /* ----------------------------------------------------------------------
    clear force on own & ghost atoms
+   clear other arrays as needed
 ------------------------------------------------------------------------- */
 
 void Respa::force_clear(int newtonflag)
 {
   if (external_force_clear) return;
 
   // clear global force array
-  // nall includes ghosts only if newton flag is set
+  // if either newton flag is set, also include ghosts
 
-  int nall;
-  if (newtonflag) nall = atom->nlocal + atom->nghost;
-  else nall = atom->nlocal;
+  size_t nbytes = sizeof(double) * atom->nlocal;
+  if (force->newton) nbytes += sizeof(double) * atom->nghost;
 
-  size_t nbytes = sizeof(double) * nall;
-
-  if (nbytes > 0 ) {
-    memset(&(atom->f[0][0]),0,3*nbytes);
-    if (torqueflag)  memset(&(atom->torque[0][0]),0,3*nbytes);
-    if (erforceflag) memset(&(atom->erforce[0]),  0,  nbytes);
-    if (e_flag)      memset(&(atom->de[0]),       0,  nbytes);
-    if (rho_flag)    memset(&(atom->drho[0]),     0,  nbytes);
+  if (nbytes) {
+    memset(&atom->f[0][0],0,3*nbytes);
+    if (torqueflag) memset(&atom->torque[0][0],0,3*nbytes);
+    if (extraflag) atom->avec->force_clear(0,nbytes);
   }
 }
 
 /* ----------------------------------------------------------------------
    copy force components from atom->f to FixRespa->f_level
 ------------------------------------------------------------------------- */
 
 void Respa::copy_f_flevel(int ilevel)
 {
   double ***f_level = fix_respa->f_level;
   double **f = atom->f;
   int n = atom->nlocal;
 
   for (int i = 0; i < n; i++) {
     f_level[i][ilevel][0] = f[i][0];
     f_level[i][ilevel][1] = f[i][1];
     f_level[i][ilevel][2] = f[i][2];
   }
 }
 
 /* ----------------------------------------------------------------------
    copy force components from FixRespa->f_level to atom->f
 ------------------------------------------------------------------------- */
 
 void Respa::copy_flevel_f(int ilevel)
 {
   double ***f_level = fix_respa->f_level;
   double **f = atom->f;
   int n = atom->nlocal;
 
   for (int i = 0; i < n; i++) {
     f[i][0] = f_level[i][ilevel][0];
     f[i][1] = f_level[i][ilevel][1];
     f[i][2] = f_level[i][ilevel][2];
   }
 }
 
 /* ----------------------------------------------------------------------
    sum all force components from FixRespa->f_level to create full atom->f
 ------------------------------------------------------------------------- */
 
 void Respa::sum_flevel_f()
 {
   copy_flevel_f(0);
 
   double ***f_level = fix_respa->f_level;
   double **f = atom->f;
   int n = atom->nlocal;
 
   for (int ilevel = 1; ilevel < nlevels; ilevel++) {
     for (int i = 0; i < n; i++) {
       f[i][0] += f_level[i][ilevel][0];
       f[i][1] += f_level[i][ilevel][1];
       f[i][2] += f_level[i][ilevel][2];
     }
   }
 }
diff --git a/src/respa.h b/src/respa.h
index 26a719867..f8dbb9d06 100644
--- a/src/respa.h
+++ b/src/respa.h
@@ -1,129 +1,128 @@
 /* -*- c++ -*- ----------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
 
    Copyright (2003) Sandia Corporation.  Under the terms of Contract
    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
    certain rights in this software.  This software is distributed under
    the GNU General Public License.
 
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
 #ifdef INTEGRATE_CLASS
 
 IntegrateStyle(respa,Respa)
 
 #else
 
 #ifndef LMP_RESPA_H
 #define LMP_RESPA_H
 
 #include "integrate.h"
 
 namespace LAMMPS_NS {
 
 class Respa : public Integrate {
  public:
                           // public so Fixes, Pairs, Neighbor can see them
   int nlevels;            // number of rRESPA levels
                           // 0 = innermost level, nlevels-1 = outermost level
   double *step;           // timestep at each level
   int *loop;              // sub-cycling factor at each level
   double cutoff[4];       // cutoff[0] and cutoff[1] = between inner and middle
                           // cutoff[2] and cutoff[3] = between middle and outer
                           // if no middle then 0,1 = 2,3
 
   int level_bond,level_angle,level_dihedral;   // level to compute forces at
   int level_improper,level_pair,level_kspace;
   int level_inner,level_middle,level_outer;
 
   Respa(class LAMMPS *, int, char **);
   virtual ~Respa();
   virtual void init();
   virtual void setup();
   virtual void setup_minimal(int);
   virtual void run(int);
   virtual void cleanup();
   virtual void reset_dt();
 
   void copy_f_flevel(int);
   void copy_flevel_f(int);
 
  protected:
   int triclinic;                    // 0 if domain is orthog, 1 if triclinic
-  int torqueflag,erforceflag;
-  int e_flag,rho_flag;
+  int torqueflag,extraflag;
 
   int *newton;                      // newton flag at each level
   class FixRespa *fix_respa;        // Fix to store the force level array
 
   virtual void recurse(int);
   void force_clear(int);
   void sum_flevel_f();
 };
 
 }
 
 #endif
 #endif
 
 /* ERROR/WARNING messages:
 
 E: Illegal ... command
 
 Self-explanatory.  Check the input script syntax and compare to the
 documentation for the command.  You can use -echo screen as a
 command-line option when running LAMMPS to see the offending line.
 
 E: Respa levels must be >= 1
 
 Self-explanatory.
 
 E: Cannot set both respa pair and inner/middle/outer
 
 In the rRESPA integrator, you must compute pairwise potentials either
 all together (pair), or in pieces (inner/middle/outer).  You can't do
 both.
 
 E: Must set both respa inner and outer
 
 Cannot use just the inner or outer option with respa without using the
 other.
 
 E: Cannot set respa middle without inner/outer
 
 In the rRESPA integrator, you must define both a inner and outer
 setting in order to use a middle setting.
 
 E: Invalid order of forces within respa levels
 
 For respa, ordering of force computations within respa levels must
 obey certain rules.  E.g. bonds cannot be compute less frequently than
 angles, pairwise forces cannot be computed less frequently than
 kspace, etc.
 
 W: One or more respa levels compute no forces
 
 This is computationally inefficient.
 
 E: Respa inner cutoffs are invalid
 
 The first cutoff must be <= the second cutoff.
 
 E: Respa middle cutoffs are invalid
 
 The first cutoff must be <= the second cutoff.
 
 W: No fixes defined, atoms won't move
 
 If you are not using a fix like nve, nvt, npt then atom velocities and
 coordinates will not be updated during timestepping.
 
 E: Pair style does not support rRESPA inner/middle/outer
 
 You are attempting to use rRESPA options with a pair style that
 does not support them.
 
 */
diff --git a/src/verlet.cpp b/src/verlet.cpp
index f6c6a3ccf..fa5c11e3d 100644
--- a/src/verlet.cpp
+++ b/src/verlet.cpp
@@ -1,427 +1,371 @@
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
 
    Copyright (2003) Sandia Corporation.  Under the terms of Contract
    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
    certain rights in this software.  This software is distributed under
    the GNU General Public License.
 
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
 #include "string.h"
 #include "verlet.h"
 #include "neighbor.h"
 #include "domain.h"
 #include "comm.h"
 #include "atom.h"
+#include "atom_vec.h"
 #include "force.h"
 #include "pair.h"
 #include "bond.h"
 #include "angle.h"
 #include "dihedral.h"
 #include "improper.h"
 #include "kspace.h"
 #include "output.h"
 #include "update.h"
 #include "modify.h"
 #include "compute.h"
 #include "fix.h"
 #include "timer.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
 Verlet::Verlet(LAMMPS *lmp, int narg, char **arg) :
   Integrate(lmp, narg, arg) {}
 
 /* ----------------------------------------------------------------------
    initialization before run
 ------------------------------------------------------------------------- */
 
 void Verlet::init()
 {
   Integrate::init();
 
   // warn if no fixes
 
   if (modify->nfix == 0 && comm->me == 0)
     error->warning(FLERR,"No fixes defined, atoms won't move");
 
   // virial_style:
   // 1 if computed explicitly by pair->compute via sum over pair interactions
   // 2 if computed implicitly by pair->virial_fdotr_compute via sum over ghosts
 
   if (force->newton_pair) virial_style = 2;
   else virial_style = 1;
 
   // setup lists of computes for global and per-atom PE and pressure
 
   ev_setup();
 
   // detect if fix omp is present for clearing force arrays
 
   int ifix = modify->find_fix("package_omp");
   if (ifix >= 0) external_force_clear = 1;
 
-  // set flags for what arrays to clear in force_clear()
-  // need to clear additionals arrays if they exist
+  // set flags for arrays to clear in force_clear()
 
-  torqueflag = 0;
+  torqueflag = extraflag = 0;
   if (atom->torque_flag) torqueflag = 1;
-  erforceflag = 0;
-  if (atom->erforce_flag) erforceflag = 1;
-  e_flag = 0;
-  if (atom->e_flag) e_flag = 1;
-  rho_flag = 0;
-  if (atom->rho_flag) rho_flag = 1;
+  if (atom->avec->forceclearflag) extraflag = 1;
 
   // orthogonal vs triclinic simulation box
 
   triclinic = domain->triclinic;
 }
 
 /* ----------------------------------------------------------------------
    setup before run
 ------------------------------------------------------------------------- */
 
 void Verlet::setup()
 {
   if (comm->me == 0 && screen) fprintf(screen,"Setting up run ...\n");
 
   update->setupflag = 1;
 
   // setup domain, communication and neighboring
   // acquire ghosts
   // build neighbor lists
 
   atom->setup();
   modify->setup_pre_exchange();
   if (triclinic) domain->x2lamda(atom->nlocal);
   domain->pbc();
   domain->reset_box();
   comm->setup();
   if (neighbor->style) neighbor->setup_bins();
   comm->exchange();
   if (atom->sortfreq > 0) atom->sort();
   comm->borders();
   if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
   domain->image_check();
   domain->box_too_small_check();
   modify->setup_pre_neighbor();
   neighbor->build();
   neighbor->ncalls = 0;
 
   // compute all forces
 
   ev_set(update->ntimestep);
   force_clear();
   modify->setup_pre_force(vflag);
 
   if (pair_compute_flag) force->pair->compute(eflag,vflag);
   else if (force->pair) force->pair->compute_dummy(eflag,vflag);
 
   if (atom->molecular) {
     if (force->bond) force->bond->compute(eflag,vflag);
     if (force->angle) force->angle->compute(eflag,vflag);
     if (force->dihedral) force->dihedral->compute(eflag,vflag);
     if (force->improper) force->improper->compute(eflag,vflag);
   }
 
   if (force->kspace) {
     force->kspace->setup();
     if (kspace_compute_flag) force->kspace->compute(eflag,vflag);
     else force->kspace->compute_dummy(eflag,vflag);
   }
 
   if (force->newton) comm->reverse_comm();
 
   modify->setup(vflag);
   output->setup();
   update->setupflag = 0;
 }
 
 /* ----------------------------------------------------------------------
    setup without output
    flag = 0 = just force calculation
    flag = 1 = reneighbor and force calculation
 ------------------------------------------------------------------------- */
 
 void Verlet::setup_minimal(int flag)
 {
   update->setupflag = 1;
 
   // setup domain, communication and neighboring
   // acquire ghosts
   // build neighbor lists
 
   if (flag) {
     modify->setup_pre_exchange();
     if (triclinic) domain->x2lamda(atom->nlocal);
     domain->pbc();
     domain->reset_box();
     comm->setup();
     if (neighbor->style) neighbor->setup_bins();
     comm->exchange();
     comm->borders();
     if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
     domain->image_check();
     domain->box_too_small_check();
     modify->setup_pre_neighbor();
     neighbor->build();
     neighbor->ncalls = 0;
   }
 
   // compute all forces
 
   ev_set(update->ntimestep);
   force_clear();
   modify->setup_pre_force(vflag);
 
   if (pair_compute_flag) force->pair->compute(eflag,vflag);
   else if (force->pair) force->pair->compute_dummy(eflag,vflag);
 
   if (atom->molecular) {
     if (force->bond) force->bond->compute(eflag,vflag);
     if (force->angle) force->angle->compute(eflag,vflag);
     if (force->dihedral) force->dihedral->compute(eflag,vflag);
     if (force->improper) force->improper->compute(eflag,vflag);
   }
 
   if (force->kspace) {
     force->kspace->setup();
     if (kspace_compute_flag) force->kspace->compute(eflag,vflag);
     else force->kspace->compute_dummy(eflag,vflag);
   }
 
   if (force->newton) comm->reverse_comm();
 
   modify->setup(vflag);
   update->setupflag = 0;
 }
 
 /* ----------------------------------------------------------------------
    run for N steps
 ------------------------------------------------------------------------- */
 
 void Verlet::run(int n)
 {
   bigint ntimestep;
   int nflag,sortflag;
 
   int n_post_integrate = modify->n_post_integrate;
   int n_pre_exchange = modify->n_pre_exchange;
   int n_pre_neighbor = modify->n_pre_neighbor;
   int n_pre_force = modify->n_pre_force;
   int n_post_force = modify->n_post_force;
   int n_end_of_step = modify->n_end_of_step;
 
   if (atom->sortfreq > 0) sortflag = 1;
   else sortflag = 0;
 
   for (int i = 0; i < n; i++) {
 
     ntimestep = ++update->ntimestep;
     ev_set(ntimestep);
 
     // initial time integration
 
     modify->initial_integrate(vflag);
     if (n_post_integrate) modify->post_integrate();
 
     // regular communication vs neighbor list rebuild
 
     nflag = neighbor->decide();
 
     if (nflag == 0) {
       timer->stamp();
       comm->forward_comm();
       timer->stamp(TIME_COMM);
     } else {
       if (n_pre_exchange) modify->pre_exchange();
       if (triclinic) domain->x2lamda(atom->nlocal);
       domain->pbc();
       if (domain->box_change) {
         domain->reset_box();
         comm->setup();
         if (neighbor->style) neighbor->setup_bins();
       }
       timer->stamp();
       comm->exchange();
       if (sortflag && ntimestep >= atom->nextsort) atom->sort();
       comm->borders();
       if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
       timer->stamp(TIME_COMM);
       if (n_pre_neighbor) modify->pre_neighbor();
       neighbor->build();
       timer->stamp(TIME_NEIGHBOR);
     }
 
     // force computations
     // important for pair to come before bonded contributions
     // since some bonded potentials tally pairwise energy/virial
     // and Pair:ev_tally() needs to be called before any tallying
 
     force_clear();
     if (n_pre_force) modify->pre_force(vflag);
 
     timer->stamp();
 
     if (pair_compute_flag) {
       force->pair->compute(eflag,vflag);
       timer->stamp(TIME_PAIR);
     }
 
     if (atom->molecular) {
       if (force->bond) force->bond->compute(eflag,vflag);
       if (force->angle) force->angle->compute(eflag,vflag);
       if (force->dihedral) force->dihedral->compute(eflag,vflag);
       if (force->improper) force->improper->compute(eflag,vflag);
       timer->stamp(TIME_BOND);
     }
 
     if (kspace_compute_flag) {
       force->kspace->compute(eflag,vflag);
       timer->stamp(TIME_KSPACE);
     }
 
     // reverse communication of forces
 
     if (force->newton) {
       comm->reverse_comm();
       timer->stamp(TIME_COMM);
     }
 
     // force modifications, final time integration, diagnostics
 
     if (n_post_force) modify->post_force(vflag);
     modify->final_integrate();
     if (n_end_of_step) modify->end_of_step();
 
     // all output
 
     if (ntimestep == output->next) {
       timer->stamp();
       output->write(ntimestep);
       timer->stamp(TIME_OUTPUT);
     }
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Verlet::cleanup()
 {
   modify->post_run();
   domain->box_too_small_check();
   update->update_time();
 }
 
 /* ----------------------------------------------------------------------
    clear force on own & ghost atoms
    clear other arrays as needed
 ------------------------------------------------------------------------- */
 
 void Verlet::force_clear()
 {
   int i;
+  size_t nbytes;
 
   if (external_force_clear) return;
 
   // clear force on all particles
   // if either newton flag is set, also include ghosts
   // when using threads always clear all forces.
 
-  if (neighbor->includegroup == 0) {
-    int nall;
-    if (force->newton) nall = atom->nlocal + atom->nghost;
-    else nall = atom->nlocal;
+  int nlocal = atom->nlocal;
 
-    size_t nbytes = sizeof(double) * nall;
+  if (neighbor->includegroup == 0) {
+    nbytes = sizeof(double) * nlocal;
+    if (force->newton) nbytes += sizeof(double) * atom->nghost;
 
     if (nbytes) {
-      memset(&(atom->f[0][0]),0,3*nbytes);
-      if (torqueflag)  memset(&(atom->torque[0][0]),0,3*nbytes);
-      if (erforceflag) memset(&(atom->erforce[0]),  0,  nbytes);
-      if (e_flag)      memset(&(atom->de[0]),       0,  nbytes);
-      if (rho_flag)    memset(&(atom->drho[0]),     0,  nbytes);
+      memset(&atom->f[0][0],0,3*nbytes);
+      if (torqueflag) memset(&atom->torque[0][0],0,3*nbytes);
+      if (extraflag) atom->avec->force_clear(0,nbytes);
     }
 
   // neighbor includegroup flag is set
   // clear force only on initial nfirst particles
   // if either newton flag is set, also include ghosts
 
   } else {
-    int nall = atom->nfirst;
-
-    double **f = atom->f;
-    for (i = 0; i < nall; i++) {
-      f[i][0] = 0.0;
-      f[i][1] = 0.0;
-      f[i][2] = 0.0;
-    }
-
-    if (torqueflag) {
-      double **torque = atom->torque;
-      for (i = 0; i < nall; i++) {
-        torque[i][0] = 0.0;
-        torque[i][1] = 0.0;
-        torque[i][2] = 0.0;
-      }
-    }
-
-    if (erforceflag) {
-      double *erforce = atom->erforce;
-      for (i = 0; i < nall; i++) erforce[i] = 0.0;
-    }
-
-    if (e_flag) {
-      double *de = atom->de;
-      for (i = 0; i < nall; i++) de[i] = 0.0;
-    }
+    nbytes = sizeof(double) * atom->nfirst;
 
-    if (rho_flag) {
-      double *drho = atom->drho;
-      for (i = 0; i < nall; i++) drho[i] = 0.0;
+    if (nbytes) {
+      memset(&atom->f[0][0],0,3*nbytes);
+      if (torqueflag) memset(&atom->torque[0][0],0,3*nbytes);
+      if (extraflag) atom->avec->force_clear(0,nbytes);
     }
 
     if (force->newton) {
-      nall = atom->nlocal + atom->nghost;
-
-      for (i = atom->nlocal; i < nall; i++) {
-        f[i][0] = 0.0;
-        f[i][1] = 0.0;
-        f[i][2] = 0.0;
-      }
-
-      if (torqueflag) {
-        double **torque = atom->torque;
-        for (i = atom->nlocal; i < nall; i++) {
-          torque[i][0] = 0.0;
-          torque[i][1] = 0.0;
-          torque[i][2] = 0.0;
-        }
-      }
-
-      if (erforceflag) {
-        double *erforce = atom->erforce;
-        for (i = atom->nlocal; i < nall; i++) erforce[i] = 0.0;
-      }
-
-      if (e_flag) {
-        double *de = atom->de;
-        for (i = 0; i < nall; i++) de[i] = 0.0;
-      }
+      nbytes = sizeof(double) * atom->nghost;
 
-      if (rho_flag) {
-        double *drho = atom->drho;
-        for (i = 0; i < nall; i++) drho[i] = 0.0;
+      if (nbytes) {
+        memset(&atom->f[nlocal][0],0,3*nbytes);
+        if (torqueflag) memset(&atom->torque[nlocal][0],0,3*nbytes);
+        if (extraflag) atom->avec->force_clear(nlocal,nbytes);
       }
     }
   }
 }
diff --git a/src/verlet.h b/src/verlet.h
index 970d9b1c1..07aaf012e 100644
--- a/src/verlet.h
+++ b/src/verlet.h
@@ -1,57 +1,56 @@
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
 
    Copyright (2003) Sandia Corporation.  Under the terms of Contract
    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
    certain rights in this software.  This software is distributed under
    the GNU General Public License.
 
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
 #ifdef INTEGRATE_CLASS
 
 IntegrateStyle(verlet,Verlet)
 
 #else
 
 #ifndef LMP_VERLET_H
 #define LMP_VERLET_H
 
 #include "integrate.h"
 
 namespace LAMMPS_NS {
 
 class Verlet : public Integrate {
  public:
   Verlet(class LAMMPS *, int, char **);
   virtual ~Verlet() {}
   virtual void init();
   virtual void setup();
   virtual void setup_minimal(int);
   virtual void run(int);
   void cleanup();
 
  protected:
   int triclinic;                    // 0 if domain is orthog, 1 if triclinic
-  int torqueflag,erforceflag;
-  int e_flag,rho_flag;
+  int torqueflag,extraflag;
 
   virtual void force_clear();
 };
 
 }
 
 #endif
 #endif
 
 /* ERROR/WARNING messages:
 
 W: No fixes defined, atoms won't move
 
 If you are not using a fix like nve, nvt, npt then atom velocities and
 coordinates will not be updated during timestepping.
 
 */