diff --git a/src/GPU/pair_lj_cut_gpu.cpp b/src/GPU/pair_lj_cut_gpu.cpp
index 29b14c033f..ea4e9532fd 100644
--- a/src/GPU/pair_lj_cut_gpu.cpp
+++ b/src/GPU/pair_lj_cut_gpu.cpp
@@ -135,27 +135,22 @@ void PairLJCutGPU::init_style()
 {
   cut_respa = nullptr;
 
-  //if (force->newton_pair)
-//    error->all(FLERR,"Cannot use newton pair with lj/cut/gpu pair style");
+  if (force->newton_pair)
+    error->all(FLERR,"Cannot use newton pair with lj/cut/gpu pair style");
 
   // Repeat cutsq calculation because done after call to init_style
   double maxcut = -1.0;
   double cut;
   for (int i = 1; i <= atom->ntypes; i++) {
     for (int j = i; j <= atom->ntypes; j++) {
-      cut = init_one(i,j);
       if (setflag[i][j] != 0 || (setflag[i][i] != 0 && setflag[j][j] != 0)) {
-      
         cut = init_one(i,j);
-        //printf("lj/cut/gpu: i = %d; j = %d: setflag = %d cut = %f\n", i, j, setflag[i][j], cut);
         cut *= cut;
         if (cut > maxcut)
           maxcut = cut;
         cutsq[i][j] = cutsq[j][i] = cut;
-      } else {
-        //printf("lj/cut/gpu: i = %d; j = %d: setflag = %d cut = %f\n", i, j, setflag[i][j], cut);
+      } else
         cutsq[i][j] = cutsq[j][i] = 0.0;
-      }
     }
   }
   double cell_size = sqrt(maxcut) + neighbor->skin;
diff --git a/src/USER-DIELECTRIC/atom_vec_dielectric.cpp b/src/USER-DIELECTRIC/atom_vec_dielectric.cpp
index 07ae1abd06..1c2a99eb19 100644
--- a/src/USER-DIELECTRIC/atom_vec_dielectric.cpp
+++ b/src/USER-DIELECTRIC/atom_vec_dielectric.cpp
@@ -26,1245 +26,83 @@ using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-AtomVecDielectric::AtomVecDielectric(LAMMPS *lmp) : AtomVecFull(lmp),
-  mu(NULL), area(NULL), ed(NULL), em(NULL), q_real(NULL), epsilon(NULL),
-  curvature(NULL)
+AtomVecDielectric::AtomVecDielectric(LAMMPS *lmp) : AtomVec(lmp)
 {
-  comm_x_only = 0;        // 0 to exchange other data in addition to x in pack comm
-  size_forward = 13;      // # items in func: pack_comm
-  size_reverse = 3;
-  size_border = 18;       // # items in func: pack_border
-  size_velocity = 3;
-  size_data_atom = 15;    // # items from read data file, func: data_atom.
-  size_data_vel = 4;
-  xcol_data = 5;          // column in data files where coordinates start
-                          // the first four columns are tag, mol, type, q
+  molecular = Atom::ATOMIC;
+  mass_type = PER_TYPE;
+  
+  atom->q_flag = atom->mu_flag = 1;
 
-  atom->mu_flag = 1;
-}
+  // strings with peratom variables to include in each AtomVec method
+  // strings cannot contain fields in corresponding AtomVec default strings
+  // order of fields in a string does not matter
+  // except: fields_data_atom & fields_data_vel must match data file
 
-/* ---------------------------------------------------------------------- */
+  fields_grow = (char *) "q mu3 area ed em epsilon curvature";
+  fields_copy = (char *) "q mu3 area ed em epsilon curvature";
+  fields_comm = (char *) "q mu3 area ed em epsilon curvature";
+  fields_comm_vel = (char *) "q mu3 area ed em epsilon curvature";
+  fields_reverse = (char *) "";
+  fields_border = (char *) "q mu3 area ed em epsilon curvature";
+  fields_border_vel = (char *) "q mu3 area ed em epsilon curvature";
+  fields_exchange = (char *) "q mu3 area ed em epsilon curvature";
+  fields_restart = (char * ) "q mu3 area ed em epsilon curvature";
+  fields_create = (char *) "q mu3 area ed em epsilon curvature";
+  fields_data_atom = (char *) "id molecule type q x mu3 area ed em epsilon curvature";
+  fields_data_vel = (char *) "id v";
 
-AtomVecDielectric::~AtomVecDielectric()
-{
-  memory->destroy(area);
-  memory->destroy(ed);
-  memory->destroy(em);
-  memory->destroy(epsilon);
-  memory->destroy(q_real);
-  memory->destroy(curvature);
+  setup_fields();
 }
 
 /* ----------------------------------------------------------------------
-   grow atom arrays
-   n = 0 grows arrays by a chunk
-   n > 0 allocates arrays to size n
+   set local copies of all grow ptrs used by this class, except defaults
+   needed in replicate when 2 atom classes exist and it calls pack_restart()
 ------------------------------------------------------------------------- */
 
-void AtomVecDielectric::grow(int n)
+void AtomVecDielectric::grow_pointers()
 {
-  AtomVecFull::grow(n);
-
-  mu = memory->grow(atom->mu,nmax,4,"atom:mu");
-
-  // USER-DIELECTRIC specifics
-
-  area = memory->grow(area,nmax,"atom:area");
-  ed = memory->grow(ed,nmax,"atom:ed");
-  em = memory->grow(em,nmax,"atom:em");
-  epsilon = memory->grow(epsilon,nmax,"atom:epsilon");
-  q_real = memory->grow(q_real,nmax,"atom:q_real");
-  curvature = memory->grow(curvature,nmax,"atom:curvature");
-}
-
-/* ----------------------------------------------------------------------
-   reset local array ptrs
-------------------------------------------------------------------------- */
-
-void AtomVecDielectric::grow_reset()
-{
-  AtomVecFull::grow_reset();
-
   mu = atom->mu;
+  area = atom->area;
+  ed = atom->ed;
+  em = atom->em;
+  epsilon = atom->epsilon;
+  curvature = atom->curvature;
+  q_unscaled = atom->q_unscaled;
 }
 
 /* ----------------------------------------------------------------------
-   copy atom I info to atom J
+   initialize non-zero atom quantities
 ------------------------------------------------------------------------- */
 
-void AtomVecDielectric::copy(int i, int j, int delflag)
+void AtomVecDielectric::create_atom_post(int ilocal)
 {
-  AtomVecFull::copy(i, j, delflag);
-
-  // USER-DIELECTRIC specifics
-
-  mu[j][0] = mu[i][0];
-  mu[j][1] = mu[i][1];
-  mu[j][2] = mu[i][2];
-  mu[j][3] = mu[i][3];
-
-  area[j] = area[i];
-  ed[j] = ed[i];
-  em[j] = em[i];
-  epsilon[j] = epsilon[i];
-  q_real[j] = q_real[i];
-  curvature[j] = curvature[i];
-}
-
-/* ---------------------------------------------------------------------- */
-
-int AtomVecDielectric::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++] = q[j];
-      buf[m++] = mu[j][0];
-      buf[m++] = mu[j][1];
-      buf[m++] = mu[j][2];
-      buf[m++] = area[j];
-      buf[m++] = ed[j];
-      buf[m++] = em[j];
-      buf[m++] = epsilon[j];
-      buf[m++] = q_real[j];
-      buf[m++] = curvature[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++] = q[j];
-      buf[m++] = mu[j][0];
-      buf[m++] = mu[j][1];
-      buf[m++] = mu[j][2];
-      buf[m++] = area[j];
-      buf[m++] = ed[j];
-      buf[m++] = em[j];
-      buf[m++] = epsilon[j];
-      buf[m++] = q_real[j];
-      buf[m++] = curvature[j];
-    }
-  }
-  return m;
-}
-
-/* ---------------------------------------------------------------------- */
-
-int AtomVecDielectric::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++] = q[j];
-      buf[m++] = mu[j][0];
-      buf[m++] = mu[j][1];
-      buf[m++] = mu[j][2];
-      buf[m++] = area[j];
-      buf[m++] = ed[j];
-      buf[m++] = em[j];
-      buf[m++] = epsilon[j];
-      buf[m++] = q_real[j];
-      buf[m++] = curvature[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++] = q[j];
-        buf[m++] = mu[j][0];
-        buf[m++] = mu[j][1];
-        buf[m++] = mu[j][2];
-        buf[m++] = area[j];
-        buf[m++] = ed[j];
-        buf[m++] = em[j];
-        buf[m++] = epsilon[j];
-        buf[m++] = q_real[j];
-        buf[m++] = curvature[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++] = q[j];
-        buf[m++] = mu[j][0];
-        buf[m++] = mu[j][1];
-        buf[m++] = mu[j][2];
-        buf[m++] = area[j];
-        buf[m++] = ed[j];
-        buf[m++] = em[j];
-        buf[m++] = epsilon[j];
-        buf[m++] = q_real[j];
-        buf[m++] = curvature[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 AtomVecDielectric::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++] = q[j];
-    buf[m++] = mu[j][0];
-    buf[m++] = mu[j][1];
-    buf[m++] = mu[j][2];
-    buf[m++] = area[j];
-    buf[m++] = ed[j];
-    buf[m++] = em[j];
-    buf[m++] = epsilon[j];
-    buf[m++] = q_real[j];
-    buf[m++] = curvature[j];
-  }
-  return m;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void AtomVecDielectric::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++];
-    q[i] = buf[m++];
-    mu[i][0] = buf[m++];
-    mu[i][1] = buf[m++];
-    mu[i][2] = buf[m++];
-    area[i] = buf[m++];
-    ed[i] = buf[m++];
-    em[i] = buf[m++];
-    epsilon[i] = buf[m++];
-    q_real[i] = buf[m++];
-    curvature[i] = buf[m++];
-  }
-}
-
-/* ---------------------------------------------------------------------- */
-
-void AtomVecDielectric::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++];
-    q[i] = buf[m++];
-    mu[i][0] = buf[m++];
-    mu[i][1] = buf[m++];
-    mu[i][2] = buf[m++];
-    area[i] = buf[m++];
-    ed[i] = buf[m++];
-    em[i] = buf[m++];
-    epsilon[i] = buf[m++];
-    q_real[i] = buf[m++];
-    curvature[i] = buf[m++];
-    v[i][0] = buf[m++];
-    v[i][1] = buf[m++];
-    v[i][2] = buf[m++];
-  }
-}
-
-/* ---------------------------------------------------------------------- */
-
-int AtomVecDielectric::unpack_comm_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++];
-    mu[i][0] = buf[m++];
-    mu[i][1] = buf[m++];
-    mu[i][2] = buf[m++];
-    area[i] = buf[m++];
-    ed[i] = buf[m++];
-    em[i] = buf[m++];
-    epsilon[i] = buf[m++];
-    q_real[i] = buf[m++];
-    curvature[i] = buf[m++];
-  }
-  return m;
-}
-
-/* ---------------------------------------------------------------------- */
-
-int AtomVecDielectric::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];
-  }
-  return m;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void AtomVecDielectric::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++];
-  }
-}
-
-/* ---------------------------------------------------------------------- */
-
-int AtomVecDielectric::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++] = ubuf(molecule[j]).d;
-      buf[m++] = q[j];
-      buf[m++] = mu[j][0];
-      buf[m++] = mu[j][1];
-      buf[m++] = mu[j][2];
-      buf[m++] = mu[j][3];
-      buf[m++] = area[j];
-      buf[m++] = ed[j];
-      buf[m++] = em[j];
-      buf[m++] = epsilon[j];
-      buf[m++] = q_real[j];
-      buf[m++] = curvature[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++] = ubuf(molecule[j]).d;
-      buf[m++] = q[j];
-      buf[m++] = mu[j][0];
-      buf[m++] = mu[j][1];
-      buf[m++] = mu[j][2];
-      buf[m++] = mu[j][3];
-      buf[m++] = area[j];
-      buf[m++] = ed[j];
-      buf[m++] = em[j];
-      buf[m++] = epsilon[j];
-      buf[m++] = q_real[j];
-      buf[m++] = curvature[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 AtomVecDielectric::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++] = ubuf(molecule[j]).d;
-      buf[m++] = q[j];
-      buf[m++] = mu[j][0];
-      buf[m++] = mu[j][1];
-      buf[m++] = mu[j][2];
-      buf[m++] = mu[j][3];
-      buf[m++] = area[j];
-      buf[m++] = ed[j];
-      buf[m++] = em[j];
-      buf[m++] = epsilon[j];
-      buf[m++] = q_real[j];
-      buf[m++] = curvature[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 (!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++] = ubuf(molecule[j]).d;
-        buf[m++] = q[j];
-        buf[m++] = mu[j][0];
-        buf[m++] = mu[j][1];
-        buf[m++] = mu[j][2];
-        buf[m++] = mu[j][3];
-        buf[m++] = area[j];
-        buf[m++] = ed[j];
-        buf[m++] = em[j];
-        buf[m++] = epsilon[j];
-        buf[m++] = q_real[j];
-        buf[m++] = curvature[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++] = ubuf(molecule[j]).d;
-        buf[m++] = q[j];
-        buf[m++] = mu[j][0];
-        buf[m++] = mu[j][1];
-        buf[m++] = mu[j][2];
-        buf[m++] = mu[j][3];
-        buf[m++] = area[j];
-        buf[m++] = ed[j];
-        buf[m++] = em[j];
-        buf[m++] = epsilon[j];
-        buf[m++] = curvature[j];
-        buf[m++] = q_real[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 AtomVecDielectric::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++] = ubuf(molecule[j]).d;
-    buf[m++] = q[j];
-    buf[m++] = mu[j][0];
-    buf[m++] = mu[j][1];
-    buf[m++] = mu[j][2];
-    buf[m++] = mu[j][3];
-    buf[m++] = area[j];
-    buf[m++] = ed[j];
-    buf[m++] = em[j];
-    buf[m++] = epsilon[j];
-    buf[m++] = q_real[j];
-    buf[m++] = curvature[j];
-  }
-  return m;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void AtomVecDielectric::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;
-    molecule[i] = (tagint) ubuf(buf[m++]).i;
-    q[i] = buf[m++];
-    mu[i][0] = buf[m++];
-    mu[i][1] = buf[m++];
-    mu[i][2] = buf[m++];
-    mu[i][3] = buf[m++];
-    area[i] = buf[m++];
-    ed[i] = buf[m++];
-    em[i] = buf[m++];
-    epsilon[i] = buf[m++];
-    q_real[i] = buf[m++];
-    curvature[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 AtomVecDielectric::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;
-    molecule[i] = (tagint) ubuf(buf[m++]).i;
-    q[i] = buf[m++];
-    mu[i][0] = buf[m++];
-    mu[i][1] = buf[m++];
-    mu[i][2] = buf[m++];
-    mu[i][3] = buf[m++];
-    area[i] = buf[m++];
-    ed[i] = buf[m++];
-    em[i] = buf[m++];
-    epsilon[i] = buf[m++];
-    q_real[i] = buf[m++];
-    curvature[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 AtomVecDielectric::unpack_border_hybrid(int n, int first, double *buf)
-{
-  int i,m,last;
-
-  m = 0;
-  last = first + n;
-  for (i = first; i < last; i++) {
-    molecule[i] = (tagint) ubuf(buf[m++]).i;
-    q[i] = buf[m++];
-    mu[i][0] = buf[m++];
-    mu[i][1] = buf[m++];
-    mu[i][2] = buf[m++];
-    mu[i][3] = buf[m++];
-    area[i] = buf[m++];
-    ed[i] = buf[m++];
-    em[i] = buf[m++];
-    epsilon[i] = buf[m++];
-    q_real[i] = buf[m++];
-    curvature[i] = buf[m++];
-  }
-  return m;
+  area[ilocal] = 1.0;
+  em[ilocal] = 1.0;
+  epsilon[ilocal] = 1.0;
 }
 
 /* ----------------------------------------------------------------------
-   pack all atom quantities for shipping to another proc
-   xyz must be 1st 3 values, so that comm::exchange can test on them
+   modify what AtomVec::data_atom() just unpacked
+   or initialize other atom quantities
 ------------------------------------------------------------------------- */
 
-int AtomVecDielectric::pack_exchange(int i, double *buf)
+void AtomVecDielectric::data_atom_post(int ilocal)
 {
-  int k;
-
-  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++] = ubuf(molecule[i]).d;
-
-  buf[m++] = q[i];
-  buf[m++] = mu[i][0];
-  buf[m++] = mu[i][1];
-  buf[m++] = mu[i][2];
-  buf[m++] = mu[i][3];
-  buf[m++] = area[i];
-  buf[m++] = ed[i];
-  buf[m++] = em[i];
-  buf[m++] = epsilon[i];
-  buf[m++] = q_real[i];
-  buf[m++] = curvature[i];
-
-  buf[m++] = ubuf(num_bond[i]).d;
-  for (k = 0; k < num_bond[i]; k++) {
-    buf[m++] = ubuf(bond_type[i][k]).d;
-    buf[m++] = ubuf(bond_atom[i][k]).d;
-  }
-
-  buf[m++] = ubuf(num_angle[i]).d;
-  for (k = 0; k < num_angle[i]; k++) {
-    buf[m++] = ubuf(angle_type[i][k]).d;
-    buf[m++] = ubuf(angle_atom1[i][k]).d;
-    buf[m++] = ubuf(angle_atom2[i][k]).d;
-    buf[m++] = ubuf(angle_atom3[i][k]).d;
-  }
-
-  buf[m++] = ubuf(num_dihedral[i]).d;
-  for (k = 0; k < num_dihedral[i]; k++) {
-    buf[m++] = ubuf(dihedral_type[i][k]).d;
-    buf[m++] = ubuf(dihedral_atom1[i][k]).d;
-    buf[m++] = ubuf(dihedral_atom2[i][k]).d;
-    buf[m++] = ubuf(dihedral_atom3[i][k]).d;
-    buf[m++] = ubuf(dihedral_atom4[i][k]).d;
-  }
-
-  buf[m++] = ubuf(num_improper[i]).d;
-  for (k = 0; k < num_improper[i]; k++) {
-    buf[m++] = ubuf(improper_type[i][k]).d;
-    buf[m++] = ubuf(improper_atom1[i][k]).d;
-    buf[m++] = ubuf(improper_atom2[i][k]).d;
-    buf[m++] = ubuf(improper_atom3[i][k]).d;
-    buf[m++] = ubuf(improper_atom4[i][k]).d;
-  }
-
-  buf[m++] = ubuf(nspecial[i][0]).d;
-  buf[m++] = ubuf(nspecial[i][1]).d;
-  buf[m++] = ubuf(nspecial[i][2]).d;
-  for (k = 0; k < nspecial[i][2]; k++) buf[m++] = ubuf(special[i][k]).d;
-
-  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 AtomVecDielectric::unpack_exchange(double *buf)
-{
-  int k;
-
-  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;
-
-  molecule[nlocal] = (tagint) ubuf(buf[m++]).i;
-
-  q[nlocal] = buf[m++];
-  mu[nlocal][0] = buf[m++];
-  mu[nlocal][1] = buf[m++];
-  mu[nlocal][2] = buf[m++];
-  mu[nlocal][3] = buf[m++];
-  area[nlocal] = buf[m++];
-  ed[nlocal] = buf[m++];
-  em[nlocal] = buf[m++];
-  epsilon[nlocal] = buf[m++];
-  q_real[nlocal] = buf[m++];
-  curvature[nlocal] = buf[m++];
-
-  num_bond[nlocal] = (int) ubuf(buf[m++]).i;
-  for (k = 0; k < num_bond[nlocal]; k++) {
-    bond_type[nlocal][k] = (int) ubuf(buf[m++]).i;
-    bond_atom[nlocal][k] = (tagint) ubuf(buf[m++]).i;
-  }
-
-  num_angle[nlocal] = (int) ubuf(buf[m++]).i;
-  for (k = 0; k < num_angle[nlocal]; k++) {
-    angle_type[nlocal][k] = (int) ubuf(buf[m++]).i;
-    angle_atom1[nlocal][k] = (tagint) ubuf(buf[m++]).i;
-    angle_atom2[nlocal][k] = (tagint) ubuf(buf[m++]).i;
-    angle_atom3[nlocal][k] = (tagint) ubuf(buf[m++]).i;
-  }
-
-  num_dihedral[nlocal] = (int) ubuf(buf[m++]).i;
-  for (k = 0; k < num_dihedral[nlocal]; k++) {
-    dihedral_type[nlocal][k] = (int) ubuf(buf[m++]).i;
-    dihedral_atom1[nlocal][k] = (tagint) ubuf(buf[m++]).i;
-    dihedral_atom2[nlocal][k] = (tagint) ubuf(buf[m++]).i;
-    dihedral_atom3[nlocal][k] = (tagint) ubuf(buf[m++]).i;
-    dihedral_atom4[nlocal][k] = (tagint) ubuf(buf[m++]).i;
-  }
-
-  num_improper[nlocal] = (int) ubuf(buf[m++]).i;
-  for (k = 0; k < num_improper[nlocal]; k++) {
-    improper_type[nlocal][k] = (int) ubuf(buf[m++]).i;
-    improper_atom1[nlocal][k] = (tagint) ubuf(buf[m++]).i;
-    improper_atom2[nlocal][k] = (tagint) ubuf(buf[m++]).i;
-    improper_atom3[nlocal][k] = (tagint) ubuf(buf[m++]).i;
-    improper_atom4[nlocal][k] = (tagint) ubuf(buf[m++]).i;
-  }
-
-  nspecial[nlocal][0] = (int) ubuf(buf[m++]).i;
-  nspecial[nlocal][1] = (int) ubuf(buf[m++]).i;
-  nspecial[nlocal][2] = (int) ubuf(buf[m++]).i;
-  for (k = 0; k < nspecial[nlocal][2]; k++)
-    special[nlocal][k] = (tagint) ubuf(buf[m++]).i;
-
-  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;
+  double* q = atom->q;
+  q_unscaled[ilocal] = q[ilocal];
+  q[ilocal] /= epsilon[ilocal];
 }
 
 /* ----------------------------------------------------------------------
-   size of restart data for all atoms owned by this proc
-   include extra data stored by fixes
+   unmodify values packed by AtomVec::pack_data()
 ------------------------------------------------------------------------- */
 
-int AtomVecDielectric::size_restart()
+void AtomVecDielectric::pack_data_post(int ilocal)
 {
-  int i;
-
-  int nlocal = atom->nlocal;
-  int n = 0;
-  // n is calculated from pack_restart: similar to full (17)
-  //   17 = 3 for x[], 4 for (tag, type, mask, image), 3 for v[], 1 for molecule ID, 1 for q,
-  //        4 for num bonds/angles/dihedrals/impropers and 1 for the first element (buf[0])
-  // plus 10 = 4 for mu (norm) and 6 (area, ed, em, epsilon, q_real, curvature) 
-  // totaling 27
-  for (i = 0; i < nlocal; i++)
-    n += 27 + 2*num_bond[i] + 4*num_angle[i] +
-      5*num_dihedral[i] + 5*num_improper[i];
-
-  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;
+  double* q = atom->q;
+  q_unscaled[ilocal] = q[ilocal]/epsilon[ilocal];
 }
 
-/* ----------------------------------------------------------------------
-   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 AtomVecDielectric::pack_restart(int i, double *buf)
-{
-  int k;
 
-  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++] = ubuf(molecule[i]).d;
 
-  buf[m++] = q[i];
-  buf[m++] = mu[i][0];
-  buf[m++] = mu[i][1];
-  buf[m++] = mu[i][2];
-  buf[m++] = mu[i][3];
-  buf[m++] = area[i];
-  buf[m++] = ed[i];
-  buf[m++] = em[i];
-  buf[m++] = epsilon[i];
-  buf[m++] = q_real[i];
-  buf[m++] = curvature[i];
-
-  buf[m++] = ubuf(num_bond[i]).d;
-  for (k = 0; k < num_bond[i]; k++) {
-    buf[m++] = ubuf(MAX(bond_type[i][k],-bond_type[i][k])).d;
-    buf[m++] = ubuf(bond_atom[i][k]).d;
-  }
-
-  buf[m++] = ubuf(num_angle[i]).d;
-  for (k = 0; k < num_angle[i]; k++) {
-    buf[m++] = ubuf(MAX(angle_type[i][k],-angle_type[i][k])).d;
-    buf[m++] = ubuf(angle_atom1[i][k]).d;
-    buf[m++] = ubuf(angle_atom2[i][k]).d;
-    buf[m++] = ubuf(angle_atom3[i][k]).d;
-  }
-
-  buf[m++] = ubuf(num_dihedral[i]).d;
-  for (k = 0; k < num_dihedral[i]; k++) {
-    buf[m++] = ubuf(MAX(dihedral_type[i][k],-dihedral_type[i][k])).d;
-    buf[m++] = ubuf(dihedral_atom1[i][k]).d;
-    buf[m++] = ubuf(dihedral_atom2[i][k]).d;
-    buf[m++] = ubuf(dihedral_atom3[i][k]).d;
-    buf[m++] = ubuf(dihedral_atom4[i][k]).d;
-  }
-
-  buf[m++] = ubuf(num_improper[i]).d;
-  for (k = 0; k < num_improper[i]; k++) {
-    buf[m++] = ubuf(MAX(improper_type[i][k],-improper_type[i][k])).d;
-    buf[m++] = ubuf(improper_atom1[i][k]).d;
-    buf[m++] = ubuf(improper_atom2[i][k]).d;
-    buf[m++] = ubuf(improper_atom3[i][k]).d;
-    buf[m++] = ubuf(improper_atom4[i][k]).d;
-  }
-
-  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 AtomVecDielectric::unpack_restart(double *buf)
-{
-  int k;
-
-  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++];
-
-  molecule[nlocal] = (tagint) ubuf(buf[m++]).i;
-
-  q[nlocal] = buf[m++];
-  mu[nlocal][0] = buf[m++];
-  mu[nlocal][1] = buf[m++];
-  mu[nlocal][2] = buf[m++];
-  mu[nlocal][3] = buf[m++];
-  area[nlocal] = buf[m++];
-  ed[nlocal] = buf[m++];
-  em[nlocal] = buf[m++];
-  epsilon[nlocal] = buf[m++];
-  q_real[nlocal] = buf[m++];
-  curvature[nlocal] = buf[m++];
-
-  num_bond[nlocal] = (int) ubuf(buf[m++]).i;
-  for (k = 0; k < num_bond[nlocal]; k++) {
-    bond_type[nlocal][k] = (int) ubuf(buf[m++]).i;
-    bond_atom[nlocal][k] = (tagint) ubuf(buf[m++]).i;
-  }
-
-  num_angle[nlocal] = (int) ubuf(buf[m++]).i;
-  for (k = 0; k < num_angle[nlocal]; k++) {
-    angle_type[nlocal][k] = (int) ubuf(buf[m++]).i;
-    angle_atom1[nlocal][k] = (tagint) ubuf(buf[m++]).i;
-    angle_atom2[nlocal][k] = (tagint) ubuf(buf[m++]).i;
-    angle_atom3[nlocal][k] = (tagint) ubuf(buf[m++]).i;
-  }
-
-  num_dihedral[nlocal] = (int) ubuf(buf[m++]).i;
-  for (k = 0; k < num_dihedral[nlocal]; k++) {
-    dihedral_type[nlocal][k] = (int) ubuf(buf[m++]).i;
-    dihedral_atom1[nlocal][k] = (tagint) ubuf(buf[m++]).i;
-    dihedral_atom2[nlocal][k] = (tagint) ubuf(buf[m++]).i;
-    dihedral_atom3[nlocal][k] = (tagint) ubuf(buf[m++]).i;
-    dihedral_atom4[nlocal][k] = (tagint) ubuf(buf[m++]).i;
-  }
-
-  num_improper[nlocal] = (int) ubuf(buf[m++]).i;
-  for (k = 0; k < num_improper[nlocal]; k++) {
-    improper_type[nlocal][k] = (int) ubuf(buf[m++]).i;
-    improper_atom1[nlocal][k] = (tagint) ubuf(buf[m++]).i;
-    improper_atom2[nlocal][k] = (tagint) ubuf(buf[m++]).i;
-    improper_atom3[nlocal][k] = (tagint) ubuf(buf[m++]).i;
-    improper_atom4[nlocal][k] = (tagint) ubuf(buf[m++]).i;
-  }
-
-  nspecial[nlocal][0] = nspecial[nlocal][1] = nspecial[nlocal][2] = 0;
-
-  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 AtomVecDielectric::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;
-  mu[nlocal][0] = 0.0;
-  mu[nlocal][1] = 0.0;
-  mu[nlocal][2] = 1.0;
-  mu[nlocal][3] = 1.0;
-  area[nlocal] = 1.0;
-  ed[nlocal] = 0.0;
-  em[nlocal] = 1.0;
-  epsilon[nlocal] = 1.0;
-  q_real[nlocal] = 0.0;
-  curvature[nlocal] = 0.0;
-
-  molecule[nlocal] = 0;
-  num_bond[nlocal] = 0;
-  num_angle[nlocal] = 0;
-  num_dihedral[nlocal] = 0;
-  num_improper[nlocal] = 0;
-  nspecial[nlocal][0] = nspecial[nlocal][1] = nspecial[nlocal][2] = 0;
-
-  atom->nlocal++;
-}
-
-/* ----------------------------------------------------------------------
-   unpack one line from Atoms section of data file
-   initialize other atom quantities
-------------------------------------------------------------------------- */
-
-void AtomVecDielectric::data_atom(double *coord, imageint imagetmp, char **values)
-{
-  int nlocal = atom->nlocal;
-  if (nlocal == nmax) grow(0);
-
-  tag[nlocal] = ATOTAGINT(values[0]);
-  molecule[nlocal] = ATOTAGINT(values[1]);
-  type[nlocal] = atoi(values[2]);
-  if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes)
-    error->one(FLERR,"Invalid atom type in Atoms section of data file");
-
-  q_real[nlocal] = atof(values[3]);
-
-  x[nlocal][0] = coord[0];
-  x[nlocal][1] = coord[1];
-  x[nlocal][2] = coord[2];
-
-  mu[nlocal][0] = atof(values[7]);
-  mu[nlocal][1] = atof(values[8]);
-  mu[nlocal][2] = atof(values[9]);
-  mu[nlocal][3] = 1;
-
-  area[nlocal] = atof(values[10]);
-  ed[nlocal] = atof(values[11]);
-  em[nlocal] = atof(values[12]);
-  epsilon[nlocal] = atof(values[13]);
-  curvature[nlocal] = atof(values[14]);
-  q[nlocal] = q_real[nlocal] / epsilon[nlocal];
-
-  image[nlocal] = imagetmp;
-
-  mask[nlocal] = 1;
-  v[nlocal][0] = 0.0;
-  v[nlocal][1] = 0.0;
-  v[nlocal][2] = 0.0;
-  num_bond[nlocal] = 0;
-  num_angle[nlocal] = 0;
-  num_dihedral[nlocal] = 0;
-  num_improper[nlocal] = 0;
-
-  atom->nlocal++;
-}
-
-/* ----------------------------------------------------------------------
-   unpack hybrid quantities from one line in Atoms section of data file
-   initialize other atom quantities for this sub-style
-------------------------------------------------------------------------- */
-
-int AtomVecDielectric::data_atom_hybrid(int nlocal, char **values)
-{
-  molecule[nlocal] = ATOTAGINT(values[0]);
-  q_real[nlocal] = atof(values[1]);
-  mu[nlocal][0] = atof(values[2]);
-  mu[nlocal][1] = atof(values[3]);
-  mu[nlocal][2] = atof(values[4]);
-  mu[nlocal][3] = 1;
-  area[nlocal] = atof(values[5]);
-  ed[nlocal] = atof(values[6]);
-  em[nlocal] = atof(values[7]);
-  epsilon[nlocal] = atof(values[8]);
-  curvature[nlocal] = atof(values[9]);
-  q[nlocal] = q_real[nlocal] / epsilon[nlocal];
-
-  num_bond[nlocal] = 0;
-  num_angle[nlocal] = 0;
-  num_dihedral[nlocal] = 0;
-  num_improper[nlocal] = 0;
-
-  return 10;
-}
-
-/* ----------------------------------------------------------------------
-   pack atom info for data file including 3 image flags
-------------------------------------------------------------------------- */
-
-void AtomVecDielectric::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(molecule[i]).d;
-    buf[i][2] = ubuf(type[i]).d;
-    buf[i][3] = q[i];
-    buf[i][4] = x[i][0];
-    buf[i][5] = x[i][1];
-    buf[i][6] = x[i][2];
-    buf[i][7] = mu[i][0];
-    buf[i][8] = mu[i][1];
-    buf[i][9] = mu[i][2];
-    buf[i][10] = area[i];
-    buf[i][11] = ed[i];
-    buf[i][12] = em[i];
-    buf[i][13] = epsilon[i];
-    buf[i][14] = curvature[i];
-    buf[i][15] = ubuf((image[i] & IMGMASK) - IMGMAX).d;
-    buf[i][16] = ubuf((image[i] >> IMGBITS & IMGMASK) - IMGMAX).d;
-    buf[i][17] = ubuf((image[i] >> IMG2BITS) - IMGMAX).d;
-  }
-}
-
-/* ----------------------------------------------------------------------
-   pack hybrid atom info for data file
-------------------------------------------------------------------------- */
-
-int AtomVecDielectric::pack_data_hybrid(int i, double *buf)
-{
-  buf[0] = ubuf(molecule[i]).d;
-  buf[1] = q[i];
-  buf[2] = mu[i][0];
-  buf[3] = mu[i][1];
-  buf[4] = mu[i][2];
-  buf[5] = area[i];
-  buf[6] = ed[i];
-  buf[7] = em[i];
-  buf[8] = epsilon[i];
-  buf[9] = q_real[i];
-  buf[10] = curvature[i];
-  return 11;
-}
-
-/* ----------------------------------------------------------------------
-   write atom info to data file including 3 image flags
-------------------------------------------------------------------------- */
-
-void AtomVecDielectric::write_data(FILE *fp, int n, double **buf)
-{
-  for (int i = 0; i < n; i++)
-    fprintf(fp,TAGINT_FORMAT " " TAGINT_FORMAT
-//            " %d %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e "
-//            "%-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %d %d %d\n",
-            " %d %g %g %g %g %g %g "
-            "%g %g %g %g %g %g %d %d %d\n",
-            (tagint) ubuf(buf[i][0]).i,(tagint) ubuf(buf[i][1]).i, (int) ubuf(buf[i][2]).i,
-            buf[i][3],buf[i][4],buf[i][5], buf[i][6],buf[i][7],buf[i][8],
-            buf[i][9],buf[i][10],buf[i][11], buf[i][12], buf[i][13], buf[i][14],
-            (int) ubuf(buf[i][15]).i,(int) ubuf(buf[i][16]).i,
-            (int) ubuf(buf[i][17]).i);
-}
-
-/* ----------------------------------------------------------------------
-   write hybrid atom info to data file
-------------------------------------------------------------------------- */
-
-int AtomVecDielectric::write_data_hybrid(FILE *fp, double *buf)
-{
-  fprintf(fp, TAGINT_FORMAT
-          " %-1.16e %-1.16e %-1.16e %-1.16e %-1.16 %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e",
-    buf[0],buf[1],buf[2],buf[3],buf[4],buf[5],buf[6],buf[7],buf[8], buf[9], buf[10]);
-  return 11;
-}
-
-/* ----------------------------------------------------------------------
-   return # of bytes of allocated memory
-------------------------------------------------------------------------- */
-
-bigint AtomVecDielectric::memory_usage()
-{
-  bigint bytes = AtomVecFull::memory_usage();
-
-  if (atom->memcheck("mu")) bytes += memory->usage(mu,nmax,4);
-
-  if (atom->memcheck("area")) bytes += memory->usage(area,nmax);
-  if (atom->memcheck("ed")) bytes += memory->usage(ed,nmax);
-  if (atom->memcheck("em")) bytes += memory->usage(em,nmax);
-  if (atom->memcheck("epsilon")) bytes += memory->usage(epsilon,nmax);
-  if (atom->memcheck("q_real")) bytes += memory->usage(q_real,nmax);
-  if (atom->memcheck("curvature")) bytes += memory->usage(curvature,nmax);
-
-  return bytes;
-}
diff --git a/src/USER-DIELECTRIC/atom_vec_dielectric.h b/src/USER-DIELECTRIC/atom_vec_dielectric.h
index bb932daeb0..797e41f112 100644
--- a/src/USER-DIELECTRIC/atom_vec_dielectric.h
+++ b/src/USER-DIELECTRIC/atom_vec_dielectric.h
@@ -20,54 +20,24 @@ AtomStyle(dielectric,AtomVecDielectric)
 #ifndef LMP_ATOM_VEC_DIELECTRIC_H
 #define LMP_ATOM_VEC_DIELECTRIC_H
 
-#include "atom_vec_full.h"
+#include "atom_vec.h"
 
 namespace LAMMPS_NS {
 
-class AtomVecDielectric : public AtomVecFull {
+class AtomVecDielectric : public AtomVec {
  public:
   AtomVecDielectric(class LAMMPS *);
-  ~AtomVecDielectric();
-  void grow(int);
-  void grow_reset();
 
-  void copy(int, int, int);
-  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 *);
-  void unpack_reverse(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 pack_data(double **);
-  int pack_data_hybrid(int, double *);
-  void write_data(FILE *, int, double **);
-  int write_data_hybrid(FILE *, double *);
-  bigint memory_usage();
+  void grow_pointers();
+  void create_atom_post(int);
+  void data_atom_post(int);
+  void pack_data_post(int);
 
- public:
-  double **mu;         // normal vector at the patch
-  double *area;        // patch area
-  double *em;          // mean dielectric constant at the patch
-  double *ed;          // difference in dielectric constants at the patch
-  double *epsilon;     // dielectric at the patch and real charges
-  double *curvature;   // curvature at the patch
-  double *q_real;      // unscaled charge: q_real = value read in from data file
+  double **mu;
+  double *area,*ed,*em,*epsilon,*curvature,*q_unscaled;
+
+ private:
+  
 };
 
 }
diff --git a/src/USER-DIELECTRIC/fix_polarize_bem_icc.cpp b/src/USER-DIELECTRIC/fix_polarize_bem_icc.cpp
index d606765a78..786ba5426a 100644
--- a/src/USER-DIELECTRIC/fix_polarize_bem_icc.cpp
+++ b/src/USER-DIELECTRIC/fix_polarize_bem_icc.cpp
@@ -21,21 +21,21 @@
      Barros, Sinkovits, Luijten, J. Chem. Phys 2014, 140, 064903
 ------------------------------------------------------------------------- */
 
-#include <cstring>
-#include <cstdlib>
 #include "fix_polarize_bem_icc.h"
-#include "atom_vec_dielectric.h"
-#include "update.h"
+
 #include "atom.h"
+#include "atom_vec_dielectric.h"
 #include "comm.h"
 #include "domain.h"
-#include "neighbor.h"
+#include "error.h"
 #include "force.h"
 #include "group.h"
 #include "kspace.h"
+#include "math_const.h"
 #include "memory.h"
 #include "modify.h"
-#include "math_const.h"
+#include "msm_dielectric.h"
+#include "neighbor.h"
 #include "pair_coul_cut_dielectric.h"
 #include "pair_coul_long_dielectric.h"
 #include "pair_lj_cut_coul_cut_dielectric.h"
@@ -45,7 +45,10 @@
 #include "msm_dielectric.h"
 #include "random_park.h"
 #include "timer.h"
-#include "error.h"
+#include "update.h"
+
+#include <cmath>
+#include <cstring>
 
 using namespace LAMMPS_NS;
 using namespace FixConst;
@@ -65,9 +68,9 @@ FixPolarizeBEMICC::FixPolarizeBEMICC(LAMMPS *lmp, int narg, char **arg) :
 
   // parse required arguments
 
-  nevery = force->inumeric(FLERR,arg[3]);
+  nevery = utils::inumeric(FLERR,arg[3],false,lmp);
   if (nevery < 0) error->all(FLERR,"Illegal fix polarize/bem/icc command");
-  double tol =  force->numeric(FLERR,arg[4]);
+  double tol =  utils::numeric(FLERR,arg[4],false,lmp);
   tol_abs = tol_rel = tol;
 
   itr_max = 20;
@@ -214,12 +217,12 @@ void FixPolarizeBEMICC::pre_force(int)
 void FixPolarizeBEMICC::compute_induced_charges()
 {
   double *q = atom->q;
-  double *q_real = avec->q_real;
-  double **norm = avec->mu;
-  double *area = avec->area;
-  double *ed = avec->ed;
-  double *em = avec->em;
-  double *epsilon = avec->epsilon;
+  double *q_real = atom->q_unscaled;
+  double **norm = atom->mu;
+  double *area = atom->area;
+  double *ed = atom->ed;
+  double *em = atom->em;
+  double *epsilon = atom->epsilon;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
   double epsilon0 = force->dielectric; 
@@ -354,11 +357,11 @@ int FixPolarizeBEMICC::modify_param(int narg, char **arg)
   while (iarg < narg) {
     if (strcmp(arg[iarg],"itr_max") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix_modify command");
-      itr_max = force->numeric(FLERR,arg[iarg+1]);
+      itr_max = utils::numeric(FLERR,arg[iarg+1],false,lmp);
       iarg += 2;
     } else if (strcmp(arg[iarg],"omega") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix_modify command");
-      omega = force->numeric(FLERR,arg[iarg+1]);
+      omega = utils::numeric(FLERR,arg[iarg+1],false,lmp);
       iarg += 2;
     } else if (strcmp(arg[iarg],"kspace") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix_modify command");
@@ -371,12 +374,12 @@ int FixPolarizeBEMICC::modify_param(int narg, char **arg)
       double epsiloni=-1, areai=-1;
       double qreali=0;
       int set_charge=0;
-      double ediff = force->numeric(FLERR,arg[iarg+1]);
-      double emean = force->numeric(FLERR,arg[iarg+2]);
-      if (strcmp(arg[iarg+3],"NULL") != 0) epsiloni = force->numeric(FLERR,arg[iarg+3]);
-      if (strcmp(arg[iarg+4],"NULL") != 0) areai = force->numeric(FLERR,arg[iarg+4]);
+      double ediff = utils::numeric(FLERR,arg[iarg+1],false,lmp);
+      double emean = utils::numeric(FLERR,arg[iarg+2],false,lmp);
+      if (strcmp(arg[iarg+3],"NULL") != 0) epsiloni = utils::numeric(FLERR,arg[iarg+3],false,lmp);
+      if (strcmp(arg[iarg+4],"NULL") != 0) areai = utils::numeric(FLERR,arg[iarg+4],false,lmp);
       if (strcmp(arg[iarg+5],"NULL") != 0) {
-        qreali = force->numeric(FLERR,arg[iarg+5]);
+        qreali = utils::numeric(FLERR,arg[iarg+5],false,lmp);
         set_charge = 1;
       }
       set_dielectric_params(ediff, emean, epsiloni, areai, set_charge, qreali);
@@ -384,8 +387,8 @@ int FixPolarizeBEMICC::modify_param(int narg, char **arg)
       iarg += 6;
     } else if (strcmp(arg[iarg],"rand") == 0) {
       if (iarg+3 > narg) error->all(FLERR,"Illegal fix_modify command");
-      ave_charge = force->numeric(FLERR,arg[iarg+1]);
-      seed_charge = force->numeric(FLERR,arg[iarg+2]);
+      ave_charge = utils::numeric(FLERR,arg[iarg+1],false,lmp);
+      seed_charge = utils::numeric(FLERR,arg[iarg+2],false,lmp);
       randomized = 1;
       iarg += 3;
     } else error->all(FLERR,"Illegal fix_modify command");
@@ -419,11 +422,11 @@ void FixPolarizeBEMICC::unpack_forward_comm(int n, int first, double *buf)
 void FixPolarizeBEMICC::set_dielectric_params(double ediff, double emean,
    double epsiloni, double areai, int set_charge, double qreali)
 {
-  double *area = avec->area;
-  double *ed = avec->ed;
-  double *em = avec->em;
-  double *q_real = avec->q_real;
-  double *epsilon = avec->epsilon;
+  double *area = atom->area;
+  double *ed = atom->ed;
+  double *em = atom->em;
+  double *q_real = atom->q_unscaled;
+  double *epsilon = atom->epsilon;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
   
diff --git a/src/USER-DIELECTRIC/msm_dielectric.cpp b/src/USER-DIELECTRIC/msm_dielectric.cpp
index 89308e5be6..b6290405c1 100644
--- a/src/USER-DIELECTRIC/msm_dielectric.cpp
+++ b/src/USER-DIELECTRIC/msm_dielectric.cpp
@@ -15,24 +15,22 @@
    Contributing authors: Trung Nguyen (Northwestern)
 ------------------------------------------------------------------------- */
 
-#include <mpi.h>
-#include <string.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>
 #include "msm_dielectric.h"
+
 #include "atom.h"
 #include "atom_vec_dielectric.h"
 #include "comm.h"
-#include "gridcomm.h"
-#include "neighbor.h"
-#include "force.h"
-#include "pair.h"
 #include "domain.h"
-#include "memory.h"
 #include "error.h"
-
+#include "force.h"
+#include "gridcomm.h"
 #include "math_const.h"
+#include "memory.h"
+#include "neighbor.h"
+#include "pair.h"
+
+#include <cstring>
+#include <cmath>
 
 using namespace LAMMPS_NS;
 using namespace MathConst;
@@ -86,7 +84,7 @@ void MSMDielectric::compute(int eflag, int vflag)
   if (scalar_pressure_flag && vflag_either) {
     if (vflag_atom)
       error->all(FLERR,"Must use 'kspace_modify pressure/scalar no' to obtain "
-        "per-atom virial with kspace_style MSMDielectric");
+        "per-atom virial with kspace_style msm/dielectric");
 
     // must switch on global energy computation if not already on
 
@@ -109,16 +107,7 @@ void MSMDielectric::compute(int eflag, int vflag)
 
   // invoke allocate_peratom() if needed for first time
 
-  if (vflag_atom && !peratom_allocate_flag) {
-    allocate_peratom();
-    cg_peratom_all->ghost_notify();
-    cg_peratom_all->setup();
-    for (int n=0; n<levels; n++) {
-      if (!active_flag[n]) continue;
-      cg_peratom[n]->ghost_notify();
-      cg_peratom[n]->setup();
-    }
-  }
+  if (vflag_atom && !peratom_allocate_flag) allocate_peratom();
 
   // convert atoms from box to lamda coords
 
@@ -146,7 +135,8 @@ void MSMDielectric::compute(int eflag, int vflag)
   // to fully sum contribution in their 3d grid
 
   current_level = 0;
-  cg_all->reverse_comm(this,REVERSE_RHO);
+  gcall->reverse_comm_kspace(this,1,sizeof(double),REVERSE_RHO,
+                             gcall_buf1,gcall_buf2,MPI_DOUBLE);
 
   // forward communicate charge density values to fill ghost grid points
   // compute direct sum interaction and then restrict to coarser grid
@@ -154,23 +144,30 @@ void MSMDielectric::compute(int eflag, int vflag)
   for (int n=0; n<=levels-2; n++) {
     if (!active_flag[n]) continue;
     current_level = n;
-    cg[n]->forward_comm(this,FORWARD_RHO);
-
+    gc[n]->forward_comm_kspace(this,1,sizeof(double),FORWARD_RHO,
+                               gc_buf1[n],gc_buf2[n],MPI_DOUBLE);
     direct(n);
     restriction(n);
   }
 
-  // compute direct interation for top grid level for nonperiodic
+  // compute direct interation for top grid level for non-periodic
   //   and for second from top grid level for periodic
 
   if (active_flag[levels-1]) {
     if (domain->nonperiodic) {
       current_level = levels-1;
-      cg[levels-1]->forward_comm(this,FORWARD_RHO);
+      gc[levels-1]->
+        forward_comm_kspace(this,1,sizeof(double),FORWARD_RHO,
+                            gc_buf1[levels-1],gc_buf2[levels-1],MPI_DOUBLE);
       direct_top(levels-1);
-      cg[levels-1]->reverse_comm(this,REVERSE_AD);
+      gc[levels-1]->
+        reverse_comm_kspace(this,1,sizeof(double),REVERSE_AD,
+                            gc_buf1[levels-1],gc_buf2[levels-1],MPI_DOUBLE);
       if (vflag_atom)
-        cg_peratom[levels-1]->reverse_comm(this,REVERSE_AD_PERATOM);
+        gc[levels-1]->
+          reverse_comm_kspace(this,6,sizeof(double),REVERSE_AD_PERATOM,
+                              gc_buf1[levels-1],gc_buf2[levels-1],MPI_DOUBLE);
+
     } else {
       // Here using MPI_Allreduce is cheaper than using commgrid
       grid_swap_forward(levels-1,qgrid[levels-1]);
@@ -178,7 +175,9 @@ void MSMDielectric::compute(int eflag, int vflag)
       grid_swap_reverse(levels-1,egrid[levels-1]);
       current_level = levels-1;
       if (vflag_atom)
-        cg_peratom[levels-1]->reverse_comm(this,REVERSE_AD_PERATOM);
+        gc[levels-1]->
+          reverse_comm_kspace(this,6,sizeof(double),REVERSE_AD_PERATOM,
+                              gc_buf1[levels-1],gc_buf2[levels-1],MPI_DOUBLE);
     }
   }
 
@@ -190,24 +189,28 @@ void MSMDielectric::compute(int eflag, int vflag)
     prolongation(n);
 
     current_level = n;
-    cg[n]->reverse_comm(this,REVERSE_AD);
+    gc[n]->reverse_comm_kspace(this,1,sizeof(double),REVERSE_AD,
+                               gc_buf1[n],gc_buf2[n],MPI_DOUBLE);
 
     // extra per-atom virial communication
 
     if (vflag_atom)
-      cg_peratom[n]->reverse_comm(this,REVERSE_AD_PERATOM);
+      gc[n]->reverse_comm_kspace(this,6,sizeof(double),REVERSE_AD_PERATOM,
+                                 gc_buf1[n],gc_buf2[n],MPI_DOUBLE);
   }
 
   // all procs communicate E-field values
   // to fill ghost cells surrounding their 3d bricks
 
   current_level = 0;
-  cg_all->forward_comm(this,FORWARD_AD);
+  gcall->forward_comm_kspace(this,1,sizeof(double),FORWARD_AD,
+                             gcall_buf1,gcall_buf2,MPI_DOUBLE);
 
   // extra per-atom energy/virial communication
 
   if (vflag_atom)
-    cg_peratom_all->forward_comm(this,FORWARD_AD_PERATOM);
+    gcall->forward_comm_kspace(this,6,sizeof(double),FORWARD_AD_PERATOM,
+                               gcall_buf1,gcall_buf2,MPI_DOUBLE);
 
   // calculate the force on my particles (interpolation)
 
@@ -266,8 +269,7 @@ void MSMDielectric::compute(int eflag, int vflag)
 
   // convert atoms back from lamda to box coords
 
-  if (triclinic)
-    domain->lamda2x(atom->nlocal);
+  if (triclinic) domain->lamda2x(atom->nlocal);
 }
 
 /* ----------------------------------------------------------------------
@@ -276,8 +278,6 @@ void MSMDielectric::compute(int eflag, int vflag)
 
 void MSMDielectric::fieldforce()
 {
-  //fprintf(screen,"MSM interpolation\n\n");
-
   double ***egridn = egrid[0];
 
   int i,l,m,n,nx,ny,nz,mx,my,mz;
@@ -296,7 +296,7 @@ void MSMDielectric::fieldforce()
   double *q = atom->q;
   double **x = atom->x;
   double **f = atom->f;
-  double *eps = avec->epsilon;
+  double *eps = atom->epsilon;
   int nlocal = atom->nlocal;
 
   for (i = 0; i < nlocal; i++) {
diff --git a/src/USER-DIELECTRIC/pair_coul_cut_dielectric.cpp b/src/USER-DIELECTRIC/pair_coul_cut_dielectric.cpp
index 3d76d475c1..2cb411e68f 100644
--- a/src/USER-DIELECTRIC/pair_coul_cut_dielectric.cpp
+++ b/src/USER-DIELECTRIC/pair_coul_cut_dielectric.cpp
@@ -74,7 +74,7 @@ void PairCoulCutDielectric::compute(int eflag, int vflag)
   double **x = atom->x;
   double **f = atom->f;
   double *q = atom->q;
-  double *q_real = avec->q_real;
+  double *q_real = avec->q_unscaled;
   double* eps = avec->epsilon;
   double** norm = avec->mu;
   double* curvature = avec->curvature;
@@ -128,7 +128,7 @@ void PairCoulCutDielectric::compute(int eflag, int vflag)
       if (rsq < cutsq[itype][jtype] && rsq > EPSILON) {
         r2inv = 1.0/rsq;
         rinv = sqrt(r2inv);
-        efield_i = qqrd2e * scale[itype][jtype] * q[j]*rinv;
+        efield_i = scale[itype][jtype] * q[j]*rinv;
         forcecoul = qtmp*efield_i;
 
         fpair_i = factor_coul*etmp*forcecoul*r2inv;
@@ -175,3 +175,27 @@ void PairCoulCutDielectric::init_style()
   neighbor->requests[irequest]->full = 1;
 }
 
+/* ---------------------------------------------------------------------- */
+
+double PairCoulCutDielectric::single(int i, int j, int itype, int jtype,
+                                double rsq,
+                                double factor_coul, double factor_lj,
+                                double &fforce)
+{
+  double r2inv,forcecoul,phicoul,ei,ej;
+  double* eps = avec->epsilon;
+
+  r2inv = 1.0/rsq;
+  forcecoul = force->qqrd2e * atom->q[i]*atom->q[j]*sqrt(r2inv)*eps[i];
+  
+  double eng = 0.0;
+  if (eps[i] == 1) ei = 0;
+  else ei = eps[i];
+  if (eps[j] == 1) ej = 0;
+  else ej = eps[j];
+  phicoul = force->qqrd2e * atom->q[i]*atom->q[j]*sqrt(r2inv);
+  phicoul *= 0.5*(ei+ej);
+  eng += factor_coul*phicoul;
+
+  return eng;
+}
\ No newline at end of file
diff --git a/src/USER-DIELECTRIC/pair_coul_cut_dielectric.h b/src/USER-DIELECTRIC/pair_coul_cut_dielectric.h
index fee444536a..2abc6b906e 100644
--- a/src/USER-DIELECTRIC/pair_coul_cut_dielectric.h
+++ b/src/USER-DIELECTRIC/pair_coul_cut_dielectric.h
@@ -29,6 +29,7 @@ class PairCoulCutDielectric : public PairCoulCut {
   PairCoulCutDielectric(class LAMMPS *);
   virtual ~PairCoulCutDielectric();
   virtual void compute(int, int);
+  virtual double single(int, int, int, int, double, double, double, double &);
   void init_style();
 
   double** efield;
diff --git a/src/USER-DIELECTRIC/pair_coul_long_dielectric.cpp b/src/USER-DIELECTRIC/pair_coul_long_dielectric.cpp
index 6a87e39fe0..6b3a38e00f 100644
--- a/src/USER-DIELECTRIC/pair_coul_long_dielectric.cpp
+++ b/src/USER-DIELECTRIC/pair_coul_long_dielectric.cpp
@@ -15,25 +15,22 @@
    Contributing author: Trung Nguyen (Northwestern)
 ------------------------------------------------------------------------- */
 
-#include <math.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
 #include "pair_coul_long_dielectric.h"
+
 #include "atom.h"
 #include "atom_vec_dielectric.h"
 #include "comm.h"
+#include "error.h"
 #include "force.h"
 #include "kspace.h"
 #include "neighbor.h"
 #include "neigh_list.h"
 #include "neigh_request.h"
-#include "update.h"
-#include "integrate.h"
-#include "respa.h"
 #include "math_const.h"
 #include "memory.h"
-#include "error.h"
+
+#include <cmath>
+#include <cstring>
 
 using namespace LAMMPS_NS;
 using namespace MathConst;
@@ -87,10 +84,10 @@ void PairCoulLongDielectric::compute(int eflag, int vflag)
   double **x = atom->x;
   double **f = atom->f;
   double *q = atom->q;
-  double** norm = avec->mu;
-  double* curvature = avec->curvature;
-  double* area = avec->area;
-  double *eps = avec->epsilon;
+  double** norm = atom->mu;
+  double* curvature = atom->curvature;
+  double* area = atom->area;
+  double *eps = atom->epsilon;
   int *type = atom->type;
   int nlocal = atom->nlocal;
   double *special_coul = force->special_coul;
diff --git a/src/USER-DIELECTRIC/pair_lj_cut_coul_cut_dielectric.cpp b/src/USER-DIELECTRIC/pair_lj_cut_coul_cut_dielectric.cpp
index 51be72acaa..56bdd579d9 100644
--- a/src/USER-DIELECTRIC/pair_lj_cut_coul_cut_dielectric.cpp
+++ b/src/USER-DIELECTRIC/pair_lj_cut_coul_cut_dielectric.cpp
@@ -15,21 +15,21 @@
    Contributing author: Trung Nguyen (Northwestern)
 ------------------------------------------------------------------------- */
 
-#include <math.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
 #include "pair_lj_cut_coul_cut_dielectric.h"
+
 #include "atom.h"
 #include "atom_vec_dielectric.h"
 #include "comm.h"
+#include "error.h"
 #include "force.h"
 #include "neighbor.h"
 #include "neigh_list.h"
 #include "neigh_request.h"
 #include "math_const.h"
 #include "memory.h"
-#include "error.h"
+
+#include <cmath>
+#include <cstring>
 
 using namespace LAMMPS_NS;
 using namespace MathConst;
@@ -78,11 +78,11 @@ void PairLJCutCoulCutDielectric::compute(int eflag, int vflag)
   double **x = atom->x;
   double **f = atom->f;
   double *q = atom->q;
-  double *q_real = avec->q_real;
-  double* eps = avec->epsilon;
-  double** norm = avec->mu;
-  double* curvature = avec->curvature;
-  double* area = avec->area;
+  double *q_real = atom->q_unscaled;
+  double* eps = atom->epsilon;
+  double** norm = atom->mu;
+  double* curvature = atom->curvature;
+  double* area = atom->area;
   int *type = atom->type;
   int nlocal = atom->nlocal;
   double *special_coul = force->special_coul;
diff --git a/src/USER-DIELECTRIC/pair_lj_cut_coul_long_dielectric.cpp b/src/USER-DIELECTRIC/pair_lj_cut_coul_long_dielectric.cpp
index 5af71e3d33..613defa2f6 100644
--- a/src/USER-DIELECTRIC/pair_lj_cut_coul_long_dielectric.cpp
+++ b/src/USER-DIELECTRIC/pair_lj_cut_coul_long_dielectric.cpp
@@ -15,25 +15,25 @@
    Contributing author: Trung Nguyen (Northwestern)
 ------------------------------------------------------------------------- */
 
-#include <math.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
 #include "pair_lj_cut_coul_long_dielectric.h"
+
 #include "atom.h"
 #include "atom_vec_dielectric.h"
 #include "comm.h"
+#include "error.h"
 #include "force.h"
-#include "kspace.h"
-#include "update.h"
 #include "integrate.h"
-#include "respa.h"
+#include "kspace.h"
 #include "neighbor.h"
 #include "neigh_list.h"
 #include "neigh_request.h"
 #include "math_const.h"
 #include "memory.h"
-#include "error.h"
+#include "respa.h"
+#include "update.h"
+
+#include <cmath>
+#include <cstring>
 
 using namespace LAMMPS_NS;
 using namespace MathConst;
@@ -95,10 +95,10 @@ void PairLJCutCoulLongDielectric::compute(int eflag, int vflag)
   double **x = atom->x;
   double **f = atom->f;
   double *q = atom->q;
-  double *eps = avec->epsilon;
-  double** norm = avec->mu;
-  double* curvature = avec->curvature;
-  double* area = avec->area;
+  double *eps = atom->epsilon;
+  double** norm = atom->mu;
+  double* curvature = atom->curvature;
+  double* area = atom->area;
   int *type = atom->type;
   int nlocal = atom->nlocal;
   double *special_coul = force->special_coul;
diff --git a/src/USER-DIELECTRIC/pair_lj_cut_coul_msm_dielectric.cpp b/src/USER-DIELECTRIC/pair_lj_cut_coul_msm_dielectric.cpp
index 3aa55a1bce..c87c46ade6 100644
--- a/src/USER-DIELECTRIC/pair_lj_cut_coul_msm_dielectric.cpp
+++ b/src/USER-DIELECTRIC/pair_lj_cut_coul_msm_dielectric.cpp
@@ -15,25 +15,25 @@
    Contributing author: Trung Nguyen (Northwestern)
 ------------------------------------------------------------------------- */
 
-#include <math.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
 #include "pair_lj_cut_coul_msm_dielectric.h"
+
 #include "atom.h"
 #include "atom_vec_dielectric.h"
 #include "comm.h"
+#include "error.h"
 #include "force.h"
 #include "kspace.h"
-#include "update.h"
 #include "integrate.h"
-#include "respa.h"
 #include "neighbor.h"
 #include "neigh_list.h"
 #include "neigh_request.h"
 #include "math_const.h"
 #include "memory.h"
-#include "error.h"
+#include "respa.h"
+#include "update.h"
+
+#include <cmath>
+#include <cstring>
 
 using namespace LAMMPS_NS;
 using namespace MathConst;
@@ -108,10 +108,10 @@ void PairLJCutCoulMSMDielectric::compute(int eflag, int vflag)
   double **x = atom->x;
   double **f = atom->f;
   double *q = atom->q;
-  double *eps = avec->epsilon;
-  double** norm = avec->mu;
-  double* curvature = avec->curvature;
-  double* area = avec->area;
+  double *eps = atom->epsilon;
+  double** norm = atom->mu;
+  double* curvature = atom->curvature;
+  double* area = atom->area;
   int *type = atom->type;
   int nlocal = atom->nlocal;
   double *special_coul = force->special_coul;
diff --git a/src/USER-DIELECTRIC/pppm_dielectric.cpp b/src/USER-DIELECTRIC/pppm_dielectric.cpp
index 084e0ab689..62779d79c4 100644
--- a/src/USER-DIELECTRIC/pppm_dielectric.cpp
+++ b/src/USER-DIELECTRIC/pppm_dielectric.cpp
@@ -13,30 +13,29 @@
 
 /* ----------------------------------------------------------------------
    Contributing authors: Trung Nguyen (Northwestern)
-     point-dipoles by Stan Moore (SNL)
 ------------------------------------------------------------------------- */
 
-#include <mpi.h>
-#include <cstring>
-#include <cstdio>
-#include <cstdlib>
-#include <cmath>
 #include "pppm_dielectric.h"
+
+#include "angle.h"
 #include "atom.h"
 #include "atom_vec_dielectric.h"
+#include "bond.h"
 #include "comm.h"
-#include "gridcomm.h"
-#include "neighbor.h"
-#include "force.h"
-#include "pair.h"
 #include "domain.h"
-#include "fft3d_wrap.h"
-#include "remap_wrap.h"
-#include "memory.h"
 #include "error.h"
-
+#include "fft3d_wrap.h"
+#include "force.h"
+#include "gridcomm.h"
 #include "math_const.h"
 #include "math_special.h"
+#include "memory.h"
+#include "neighbor.h"
+#include "pair.h"
+#include "remap_wrap.h"
+
+#include <cmath>
+#include <cstring>
 
 using namespace LAMMPS_NS;
 using namespace MathConst;
@@ -44,6 +43,9 @@ using namespace MathSpecial;
 
 #define SMALL 0.00001
 
+enum{REVERSE_RHO};
+enum{FORWARD_IK,FORWARD_AD,FORWARD_IK_PERATOM,FORWARD_AD_PERATOM};
+
 #ifdef FFT_SINGLE
 #define ZEROF 0.0f
 #define ONEF  1.0f
@@ -52,23 +54,12 @@ using namespace MathSpecial;
 #define ONEF  1.0
 #endif
 
-enum{REVERSE_RHO,REVERSE_MU};
-enum{FORWARD_IK,FORWARD_MU,FORWARD_AD,FORWARD_IK_PERATOM,FORWARD_MU_PERATOM,FORWARD_AD_PERATOM};
-
 /* ---------------------------------------------------------------------- */
 
-PPPMDielectric::PPPMDielectric(LAMMPS *lmp) : PPPM(lmp), 
-  cg_mu(NULL), densityx_brick_dipole(NULL), densityy_brick_dipole(NULL), densityz_brick_dipole(NULL),
-  u_brick_dipole(NULL), ux_brick_dipole(NULL), uy_brick_dipole(NULL), uz_brick_dipole(NULL),
-  vdxx_brick_dipole(NULL), vdxy_brick_dipole(NULL),  vdyy_brick_dipole(NULL), vdxz_brick_dipole(NULL),
-  vdyz_brick_dipole(NULL), vdzz_brick_dipole(NULL), work3(NULL), work4(NULL),
-  densityx_fft_dipole(NULL), densityy_fft_dipole(NULL), densityz_fft_dipole(NULL)
+PPPMDielectric::PPPMDielectric(LAMMPS *lmp) : PPPM(lmp)
 {
-  dipoleflag = 0; // turned off for now, until dipole works
   group_group_enable = 0;
 
-  mu_flag = 0;
-
   efield = NULL;
   phi = NULL;
   potflag = 0;
@@ -85,44 +76,6 @@ PPPMDielectric::~PPPMDielectric()
   memory->destroy(phi);
 }
 
-/* ----------------------------------------------------------------------
-   called once before run
-------------------------------------------------------------------------- */
-
-void PPPMDielectric::init()
-{
-  PPPM::init();
-
-  if (mu_flag) musum_musq();
-
-  if (mu_flag) {
-    cg_mu->ghost_notify();
-    cg_mu->setup();
-  }
-}
-
-/* ----------------------------------------------------------------------
-   reset local grid arrays and communication stencils
-   called by fix balance b/c it changed sizes of processor sub-domains
-------------------------------------------------------------------------- */
-
-void PPPMDielectric::setup_grid()
-{
-  PPPM::setup_grid();
-
-  if (mu_flag) {
-    cg_mu->ghost_notify();
-    if (overlap_allowed == 0 && cg_mu->ghost_overlap())
-      error->all(FLERR,"PPPMDielectric grid stencil extends "
-                 "beyond nearest neighbor processor");
-    cg_mu->setup();
-  }
-
-  // pre-compute volume-dependent coeffs
-
-  PPPM::setup();
-}
-
 /* ----------------------------------------------------------------------
    compute the PPPMDielectric long-range force, energy, virial
 ------------------------------------------------------------------------- */
@@ -134,29 +87,20 @@ void PPPMDielectric::compute(int eflag, int vflag)
   // set energy/virial flags
   // invoke allocate_peratom() if needed for first time
 
-  if (eflag || vflag) ev_setup(eflag,vflag);
-  else evflag = evflag_atom = eflag_global = vflag_global =
-         eflag_atom = vflag_atom = 0;
+  ev_init(eflag,vflag);
 
-  if (potflag) evflag_atom = 1;
-
-  if (evflag_atom && !peratom_allocate_flag) {
-    allocate_peratom();
-    cg_peratom->ghost_notify();
-    cg_peratom->setup();
-  }
+  if (evflag_atom && !peratom_allocate_flag) allocate_peratom();
 
   // if atom count has changed, update qsum and qsqsum
 
   if (atom->natoms != natoms_original) {
     qsum_qsq();
-    musum_musq();
     natoms_original = atom->natoms;
   }
 
   // return if there are no charges or dipoles
 
-  //if (qsqsum == 0.0 && musqsum == 0.0) return;
+  if (qsqsum == 0.0) return;
 
   // convert atoms from box to lamda coords
 
@@ -184,21 +128,14 @@ void PPPMDielectric::compute(int eflag, int vflag)
   particle_map();
   make_rho();
 
-  if (mu_flag)
-    make_rho_dipole();
-
   // all procs communicate density values from their ghost cells
   //   to fully sum contribution in their 3d bricks
   // remap from 3d decomposition to FFT decomposition
 
-  cg->reverse_comm(this,REVERSE_RHO);
+  gc->reverse_comm_kspace(this,1,sizeof(FFT_SCALAR),REVERSE_RHO,
+                          gc_buf1,gc_buf2,MPI_FFT_SCALAR);
   brick2fft();
 
-  if (mu_flag) {
-    cg_mu->reverse_comm(this,REVERSE_MU);
-  	brick2fft_dipole();
-  }
-
   // compute potential gradient on my FFT grid and
   //   portion of e_long on this proc's FFT grid
   // return gradients (electric fields) in 3d brick decomposition
@@ -206,34 +143,31 @@ void PPPMDielectric::compute(int eflag, int vflag)
 
   poisson();
 
-  if (mu_flag)
-    poisson_ik_dipole();
-
   // all procs communicate E-field values
   // to fill ghost cells surrounding their 3d bricks
 
-  if (differentiation_flag == 1) cg->forward_comm(this,FORWARD_AD);
-  else cg->forward_comm(this,FORWARD_IK);
-
-  if (mu_flag)
-    cg_mu->forward_comm(this,FORWARD_MU);
+  if (differentiation_flag == 1)
+    gc->forward_comm_kspace(this,1,sizeof(FFT_SCALAR),FORWARD_AD,
+                            gc_buf1,gc_buf2,MPI_FFT_SCALAR);
+  else
+    gc->forward_comm_kspace(this,3,sizeof(FFT_SCALAR),FORWARD_IK,
+                            gc_buf1,gc_buf2,MPI_FFT_SCALAR);
 
   // extra per-atom energy/virial communication
 
   if (evflag_atom) {
     if (differentiation_flag == 1 && vflag_atom)
-      cg_peratom->forward_comm(this,FORWARD_AD_PERATOM);
+      gc->forward_comm_kspace(this,6,sizeof(FFT_SCALAR),FORWARD_AD_PERATOM,
+                              gc_buf1,gc_buf2,MPI_FFT_SCALAR);
     else if (differentiation_flag == 0)
-      cg_peratom->forward_comm(this,FORWARD_IK_PERATOM);
+      gc->forward_comm_kspace(this,7,sizeof(FFT_SCALAR),FORWARD_IK_PERATOM,
+                              gc_buf1,gc_buf2,MPI_FFT_SCALAR);
   }
 
   // calculate the force on my particles
 
   fieldforce();
 
-  if (mu_flag)
-    fieldforce_ik_dipole();
-
   // extra per-atom energy/virial communication
 
   if (evflag_atom) fieldforce_peratom();
@@ -241,7 +175,6 @@ void PPPMDielectric::compute(int eflag, int vflag)
   // sum global energy across procs and add in volume-dependent term
 
   const double qscale = qqrd2e * scale;
-  const double g3 = g_ewald*g_ewald*g_ewald;
 
   if (eflag_global) {
     double energy_all;
@@ -251,8 +184,6 @@ void PPPMDielectric::compute(int eflag, int vflag)
     energy *= 0.5*volume;
     energy -= g_ewald*qsqsum/MY_PIS +
       MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);
-    if (mu_flag)
-      energy -= musqsum*qqrd2e*2.0*g3/3.0/MY_PIS;
     energy *= qscale;
   }
 
@@ -270,7 +201,6 @@ void PPPMDielectric::compute(int eflag, int vflag)
 
   if (evflag_atom) {
     double *q = atom->q;
-    double **mu = atom->mu;
     int nlocal = atom->nlocal;
     int ntotal = nlocal;
     if (tip4pflag) ntotal += atom->nghost;
@@ -280,10 +210,6 @@ void PPPMDielectric::compute(int eflag, int vflag)
         eatom[i] *= 0.5;
         eatom[i] -= g_ewald*q[i]*q[i]/MY_PIS + MY_PI2*q[i]*qsum /
           (g_ewald*g_ewald*volume);
-
-      if (mu_flag)
-        eatom[i] -= (mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] + mu[i][2]*mu[i][2])*qqrd2e*2.0*g3/3.0/MY_PIS;
-
         eatom[i] *= qscale;
       }
       for (i = nlocal; i < ntotal; i++) eatom[i] *= 0.5*qscale;
@@ -304,491 +230,6 @@ void PPPMDielectric::compute(int eflag, int vflag)
   if (triclinic) domain->lamda2x(atom->nlocal);
 }
 
-/* ----------------------------------------------------------------------
-   allocate memory that depends on # of K-vectors and order
-------------------------------------------------------------------------- */
-
-void PPPMDielectric::allocate()
-{
-  PPPM::allocate();
-
-  if (mu_flag) {
-    memory->create3d_offset(densityx_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out,
-                            nxlo_out,nxhi_out,"pppm:densityx_brick_dipole");
-    memory->create3d_offset(densityy_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out,
-                            nxlo_out,nxhi_out,"pppm:densityy_brick_dipole");
-    memory->create3d_offset(densityz_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out,
-                            nxlo_out,nxhi_out,"pppm:densityz_brick_dipole");
-
-    memory->create(densityx_fft_dipole,nfft_both,"pppm:densityy_fft_dipole");
-    memory->create(densityy_fft_dipole,nfft_both,"pppm:densityy_fft_dipole");
-    memory->create(densityz_fft_dipole,nfft_both,"pppm:densityz_fft_dipole");
-
-    memory->create(work3,2*nfft_both,"pppm:work3");
-    memory->create(work4,2*nfft_both,"pppm:work4");
-
-    memory->create3d_offset(u_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out,
-                            nxlo_out,nxhi_out,"pppm:u_brick_dipole");
-
-    memory->create3d_offset(ux_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out,
-                            nxlo_out,nxhi_out,"pppm:ux_brick_dipole");
-    memory->create3d_offset(uy_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out,
-                            nxlo_out,nxhi_out,"pppm:uy_brick_dipole");
-    memory->create3d_offset(uz_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out,
-                            nxlo_out,nxhi_out,"pppm:uz_brick_dipole");
-
-    memory->create3d_offset(vdxx_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out,
-                            nxlo_out,nxhi_out,"pppm:vdxx_brick_dipole");
-    memory->create3d_offset(vdxy_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out,
-                            nxlo_out,nxhi_out,"pppm:vdxy_brick_dipole");
-    memory->create3d_offset(vdyy_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out,
-                            nxlo_out,nxhi_out,"pppm:vdyy_brick_dipole");
-    memory->create3d_offset(vdxz_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out,
-                            nxlo_out,nxhi_out,"pppm:vdxz_brick_dipole");
-    memory->create3d_offset(vdyz_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out,
-                            nxlo_out,nxhi_out,"pppm:vdyz_brick_dipole");
-    memory->create3d_offset(vdzz_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out,
-                            nxlo_out,nxhi_out,"pppm:vdzz_brick_dipole");
-
-    int (*procneigh)[2] = comm->procneigh;
-    cg_mu = new GridComm(lmp,world,9,3,
-                         nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in,
-                         nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out,
-                         procneigh[0][0],procneigh[0][1],procneigh[1][0],
-                         procneigh[1][1],procneigh[2][0],procneigh[2][1]);
-  }
-}
-
-/* ----------------------------------------------------------------------
-   deallocate memory that depends on # of K-vectors and order
-------------------------------------------------------------------------- */
-
-void PPPMDielectric::deallocate()
-{
-  PPPM::deallocate();
-
-  if (mu_flag) {
-    memory->destroy3d_offset(densityx_brick_dipole,nzlo_out,nylo_out,nxlo_out);
-    memory->destroy3d_offset(densityy_brick_dipole,nzlo_out,nylo_out,nxlo_out);
-    memory->destroy3d_offset(densityz_brick_dipole,nzlo_out,nylo_out,nxlo_out);
-
-    memory->destroy3d_offset(u_brick_dipole,nzlo_out,nylo_out,nxlo_out);
-
-    memory->destroy3d_offset(ux_brick_dipole,nzlo_out,nylo_out,nxlo_out);
-    memory->destroy3d_offset(uy_brick_dipole,nzlo_out,nylo_out,nxlo_out);
-    memory->destroy3d_offset(uz_brick_dipole,nzlo_out,nylo_out,nxlo_out);
-
-    memory->destroy3d_offset(vdxx_brick_dipole,nzlo_out,nylo_out,nxlo_out);
-    memory->destroy3d_offset(vdxy_brick_dipole,nzlo_out,nylo_out,nxlo_out);
-    memory->destroy3d_offset(vdyy_brick_dipole,nzlo_out,nylo_out,nxlo_out);
-    memory->destroy3d_offset(vdxz_brick_dipole,nzlo_out,nylo_out,nxlo_out);
-    memory->destroy3d_offset(vdyz_brick_dipole,nzlo_out,nylo_out,nxlo_out);
-    memory->destroy3d_offset(vdzz_brick_dipole,nzlo_out,nylo_out,nxlo_out);
-
-    memory->destroy(densityx_fft_dipole);
-    memory->destroy(densityy_fft_dipole);
-    memory->destroy(densityz_fft_dipole);
-
-    memory->destroy(work3);
-    memory->destroy(work4);
-
-    delete cg_mu;
-  }
-}
-
-/* ----------------------------------------------------------------------
-   create discretized "density" on section of global grid due to my particles
-   density(x,y,z) = charge "density" at grid points of my 3d brick
-   (nxlo:nxhi,nylo:nyhi,nzlo:nzhi) is extent of my brick (including ghosts)
-   in global grid
-------------------------------------------------------------------------- */
-
-void PPPMDielectric::make_rho_dipole()
-{
-  int l,m,n,nx,ny,nz,mx,my,mz;
-  FFT_SCALAR dx,dy,dz;
-  FFT_SCALAR x0,y0,z0;
-  FFT_SCALAR x1,y1,z1;
-  FFT_SCALAR x2,y2,z2;
-
-  // clear 3d density array
-
-  memset(&(densityx_brick_dipole[nzlo_out][nylo_out][nxlo_out]),0,
-         ngrid*sizeof(FFT_SCALAR));
-  memset(&(densityy_brick_dipole[nzlo_out][nylo_out][nxlo_out]),0,
-         ngrid*sizeof(FFT_SCALAR));
-  memset(&(densityz_brick_dipole[nzlo_out][nylo_out][nxlo_out]),0,
-         ngrid*sizeof(FFT_SCALAR));
-
-  // loop over my charges, add their contribution to nearby grid points
-  // (nx,ny,nz) = global coords of grid pt to "lower left" of charge
-  // (dx,dy,dz) = distance to "lower left" grid pt
-  // (mx,my,mz) = global coords of moving stencil pt
-
-  double **mu = atom->mu;
-  double **x = atom->x;
-  int nlocal = atom->nlocal;
-
-  for (int i = 0; i < nlocal; i++) {
-
-    nx = part2grid[i][0];
-    ny = part2grid[i][1];
-    nz = part2grid[i][2];
-    dx = nx+shiftone - (x[i][0]-boxlo[0])*delxinv;
-    dy = ny+shiftone - (x[i][1]-boxlo[1])*delyinv;
-    dz = nz+shiftone - (x[i][2]-boxlo[2])*delzinv;
-
-    compute_rho1d(dx,dy,dz);
-
-    z0 = delvolinv * mu[i][0];
-    z1 = delvolinv * mu[i][1];
-    z2 = delvolinv * mu[i][2];
-    for (n = nlower; n <= nupper; n++) {
-      mz = n+nz;
-      y0 = z0*rho1d[2][n];
-      y1 = z1*rho1d[2][n];
-      y2 = z2*rho1d[2][n];
-      for (m = nlower; m <= nupper; m++) {
-        my = m+ny;
-        x0 = y0*rho1d[1][m];
-        x1 = y1*rho1d[1][m];
-        x2 = y2*rho1d[1][m];
-        for (l = nlower; l <= nupper; l++) {
-          mx = l+nx;
-          densityx_brick_dipole[mz][my][mx] += x0*rho1d[0][l];
-          densityy_brick_dipole[mz][my][mx] += x1*rho1d[0][l];
-          densityz_brick_dipole[mz][my][mx] += x2*rho1d[0][l];
-        }
-      }
-    }
-  }
-}
-
-/* ----------------------------------------------------------------------
-   remap density from 3d brick decomposition to FFT decomposition
-------------------------------------------------------------------------- */
-
-void PPPMDielectric::brick2fft_dipole()
-{
-  int n,ix,iy,iz;
-
-  // copy grabs inner portion of density from 3d brick
-  // remap could be done as pre-stage of FFT,
-  //   but this works optimally on only double values, not complex values
-
-  n = 0;
-  for (iz = nzlo_in; iz <= nzhi_in; iz++)
-    for (iy = nylo_in; iy <= nyhi_in; iy++)
-      for (ix = nxlo_in; ix <= nxhi_in; ix++) {
-        densityx_fft_dipole[n] = densityx_brick_dipole[iz][iy][ix];
-        densityy_fft_dipole[n] = densityy_brick_dipole[iz][iy][ix];
-        densityz_fft_dipole[n] = densityz_brick_dipole[iz][iy][ix];
-        n++;
-      }
-
-  remap->perform(densityx_fft_dipole,densityx_fft_dipole,work1);
-  remap->perform(densityy_fft_dipole,densityy_fft_dipole,work1);
-  remap->perform(densityz_fft_dipole,densityz_fft_dipole,work1);
-}
-
-/* ----------------------------------------------------------------------
-   FFT-based Poisson solver for ik
-------------------------------------------------------------------------- */
-
-void PPPMDielectric::poisson_ik_dipole()
-{
-  int i,j,k,n,ii;
-  double eng;
-  double wreal,wimg;
-
-  // transform dipole density (r -> k)
-
-  n = 0;
-  for (i = 0; i < nfft; i++) {
-    work1[n] = densityx_fft_dipole[i];
-    work1[n+1] = ZEROF;
-    work2[n] = densityy_fft_dipole[i];
-    work2[n+1] = ZEROF;
-    work3[n] = densityz_fft_dipole[i];
-    work3[n+1] = ZEROF;
-    n += 2;
-  }
-
-  fft1->compute(work1,work1,1);
-  fft1->compute(work2,work2,1);
-  fft1->compute(work3,work3,1);
-
-  // global energy and virial contribution
-
-  double scaleinv = 1.0/(nx_pppm*ny_pppm*nz_pppm);
-  double s2 = scaleinv*scaleinv;
-
-  if (eflag_global || vflag_global) {
-    if (vflag_global) {
-      n = 0;
-      ii = 0;
-      for (k = nzlo_fft; k <= nzhi_fft; k++)
-        for (j = nylo_fft; j <= nyhi_fft; j++)
-          for (i = nxlo_fft; i <= nxhi_fft; i++) {
-            wreal = (work1[n]*fkx[i] + work2[n]*fky[j] + work3[n]*fkz[k]);
-            wimg = (work1[n+1]*fkx[i] + work2[n+1]*fky[j] + work3[n+1]*fkz[k]);
-            eng = s2 * greensfn[ii] * (wreal*wreal + wimg*wimg);
-            //eng = s2 * greensfn[ii] * (wreal+wimg)*(wreal+wimg);
-            for (int jj = 0; jj < 6; jj++) virial[jj] += eng*vg[ii][jj];
-            if (eflag_global) energy += eng;
-            ii++;
-            n += 2;
-          }
-    } else {
-      n = 0;
-      ii = 0;
-      for (k = nzlo_fft; k <= nzhi_fft; k++)
-        for (j = nylo_fft; j <= nyhi_fft; j++)
-          for (i = nxlo_fft; i <= nxhi_fft; i++) {
-            wreal = (work1[n]*fkx[i] + work2[n]*fky[j] + work3[n]*fkz[k]);
-            wimg = (work1[n+1]*fkx[i] + work2[n+1]*fky[j] + work3[n+1]*fkz[k]);
-            energy +=
-            //s2 * greensfn[ii] * (wreal+wimg)*(wreal+wimg);
-            s2 * greensfn[ii] * (wreal*wreal + wimg*wimg);
-	          ii++;
-            n += 2;
-          }
-    }
-  }
-
-  // scale by 1/total-grid-pts to get rho(k)
-  // multiply by Green's function to get V(k)
-
-  n = 0;
-  for (i = 0; i < nfft; i++) {
-    work1[n]   *= scaleinv * greensfn[i];
-    work1[n+1] *= scaleinv * greensfn[i];
-    work2[n]   *= scaleinv * greensfn[i];
-    work2[n+1] *= scaleinv * greensfn[i];
-    work3[n]   *= scaleinv * greensfn[i];
-    work3[n+1] *= scaleinv * greensfn[i];
-    n += 2;
-  }
-
-  // triclinic system
-
-  /*if (triclinic) {
-    poisson_ik_triclinic();
-    return;
-  }*/
-
-  // compute electric potential
-  // FFT leaves data in 3d brick decomposition
-
-  n = 0;
-  for (k = nzlo_fft; k <= nzhi_fft; k++)
-    for (j = nylo_fft; j <= nyhi_fft; j++)
-      for (i = nxlo_fft; i <= nxhi_fft; i++) {
-        work4[n] = work1[n]*fkx[i] + work2[n]*fky[j] + work3[n]*fkz[k];
-        work4[n+1] = work1[n+1]*fkx[i] + work2[n+1]*fky[j] + work3[n+1]*fkz[k];
-        n += 2;
-      }
-
-  fft2->compute(work4,work4,-1);
-
-  n = 0;
-  for (k = nzlo_in; k <= nzhi_in; k++)
-    for (j = nylo_in; j <= nyhi_in; j++)
-      for (i = nxlo_in; i <= nxhi_in; i++) {
-        u_brick_dipole[k][j][i] = work4[n];
-        n += 2;
-      }
-
-  // Ex
-
-  n = 0;
-  for (k = nzlo_fft; k <= nzhi_fft; k++)
-    for (j = nylo_fft; j <= nyhi_fft; j++)
-      for (i = nxlo_fft; i <= nxhi_fft; i++) {
-        work4[n] = fkx[i]*(work1[n]*fkx[i] + work2[n]*fky[j] + work3[n]*fkz[k]);
-        work4[n+1] = fkx[i]*(work1[n+1]*fkx[i] + work2[n+1]*fky[j] + work3[n+1]*fkz[k]);
-        n += 2;
-      }
-
-  fft2->compute(work4,work4,-1);
-
-  n = 0;
-  for (k = nzlo_in; k <= nzhi_in; k++)
-    for (j = nylo_in; j <= nyhi_in; j++)
-      for (i = nxlo_in; i <= nxhi_in; i++) {
-        ux_brick_dipole[k][j][i] = work4[n];
-        n += 2;
-      }
-
-  // Ey
-
-  n = 0;
-  for (k = nzlo_fft; k <= nzhi_fft; k++)
-    for (j = nylo_fft; j <= nyhi_fft; j++)
-      for (i = nxlo_fft; i <= nxhi_fft; i++) {
-        work4[n] = fky[j]*(work1[n]*fkx[i] + work2[n]*fky[j] + work3[n]*fkz[k]);
-        work4[n+1] = fky[j]*(work1[n+1]*fkx[i] + work2[n+1]*fky[j] + work3[n+1]*fkz[k]);
-        n += 2;
-      }
-
-  fft2->compute(work4,work4,-1);
-
-  n = 0;
-  for (k = nzlo_in; k <= nzhi_in; k++)
-    for (j = nylo_in; j <= nyhi_in; j++)
-      for (i = nxlo_in; i <= nxhi_in; i++) {
-        uy_brick_dipole[k][j][i] = work4[n];
-        n += 2;
-      }
-
-  // Ez
-
-  n = 0;
-  for (k = nzlo_fft; k <= nzhi_fft; k++)
-    for (j = nylo_fft; j <= nyhi_fft; j++)
-      for (i = nxlo_fft; i <= nxhi_fft; i++) {
-        work4[n] = fkz[k]*(work1[n]*fkx[i] + work2[n]*fky[j] + work3[n]*fkz[k]);
-        work4[n+1] = fkz[k]*(work1[n+1]*fkx[i] + work2[n+1]*fky[j] + work3[n+1]*fkz[k]);
-        n += 2;
-      }
-
-  fft2->compute(work4,work4,-1);
-
-  n = 0;
-  for (k = nzlo_in; k <= nzhi_in; k++)
-    for (j = nylo_in; j <= nyhi_in; j++)
-      for (i = nxlo_in; i <= nxhi_in; i++) {
-        uz_brick_dipole[k][j][i] = work4[n];
-        n += 2;
-      }
-
-  // Vxx
-
-  n = 0;
-  for (k = nzlo_fft; k <= nzhi_fft; k++)
-    for (j = nylo_fft; j <= nyhi_fft; j++)
-      for (i = nxlo_fft; i <= nxhi_fft; i++) {
-        work4[n] = fkx[i]*fkx[i]*(work1[n+1]*fkx[i] + work2[n+1]*fky[j] + work3[n+1]*fkz[k]);
-        work4[n+1] = -fkx[i]*fkx[i]*(work1[n]*fkx[i] + work2[n]*fky[j] + work3[n]*fkz[k]);
-        n += 2;
-      }
-
-  fft2->compute(work4,work4,-1);
-
-  n = 0;
-  for (k = nzlo_in; k <= nzhi_in; k++)
-    for (j = nylo_in; j <= nyhi_in; j++)
-      for (i = nxlo_in; i <= nxhi_in; i++) {
-        vdxx_brick_dipole[k][j][i] = work4[n];
-        n += 2;
-      }
-
-  // Vyy
-
-  n = 0;
-  for (k = nzlo_fft; k <= nzhi_fft; k++)
-    for (j = nylo_fft; j <= nyhi_fft; j++)
-      for (i = nxlo_fft; i <= nxhi_fft; i++) {
-        work4[n] = fky[j]*fky[j]*(work1[n+1]*fkx[i] + work2[n+1]*fky[j] + work3[n+1]*fkz[k]);
-        work4[n+1] = -fky[j]*fky[j]*(work1[n]*fkx[i] + work2[n]*fky[j] + work3[n]*fkz[k]);
-        n += 2;
-      }
-
-  fft2->compute(work4,work4,-1);
-
-  n = 0;
-  for (k = nzlo_in; k <= nzhi_in; k++)
-    for (j = nylo_in; j <= nyhi_in; j++)
-      for (i = nxlo_in; i <= nxhi_in; i++) {
-        vdyy_brick_dipole[k][j][i] = work4[n];
-        n += 2;
-      }
-
-  // Vzz
-
-  n = 0;
-  for (k = nzlo_fft; k <= nzhi_fft; k++)
-    for (j = nylo_fft; j <= nyhi_fft; j++)
-      for (i = nxlo_fft; i <= nxhi_fft; i++) {
-        work4[n] = fkz[k]*fkz[k]*(work1[n+1]*fkx[i] + work2[n+1]*fky[j] + work3[n+1]*fkz[k]);
-        work4[n+1] = -fkz[k]*fkz[k]*(work1[n]*fkx[i] + work2[n]*fky[j] + work3[n]*fkz[k]);
-        n += 2;
-      }
-
-  fft2->compute(work4,work4,-1);
-
-  n = 0;
-  for (k = nzlo_in; k <= nzhi_in; k++)
-    for (j = nylo_in; j <= nyhi_in; j++)
-      for (i = nxlo_in; i <= nxhi_in; i++) {
-        vdzz_brick_dipole[k][j][i] = work4[n];
-        n += 2;
-      }
-
-  // Vxy
-
-  n = 0;
-  for (k = nzlo_fft; k <= nzhi_fft; k++)
-    for (j = nylo_fft; j <= nyhi_fft; j++)
-      for (i = nxlo_fft; i <= nxhi_fft; i++) {
-        work4[n] = fkx[i]*fky[j]*(work1[n+1]*fkx[i] + work2[n+1]*fky[j] + work3[n+1]*fkz[k]);
-        work4[n+1] = -fkx[i]*fky[j]*(work1[n]*fkx[i] + work2[n]*fky[j] + work3[n]*fkz[k]);
-        n += 2;
-      }
-
-  fft2->compute(work4,work4,-1);
-
-  n = 0;
-  for (k = nzlo_in; k <= nzhi_in; k++)
-    for (j = nylo_in; j <= nyhi_in; j++)
-      for (i = nxlo_in; i <= nxhi_in; i++) {
-        vdxy_brick_dipole[k][j][i] = work4[n];
-        n += 2;
-      }
-
-  // Vxz
-
-  n = 0;
-  for (k = nzlo_fft; k <= nzhi_fft; k++)
-    for (j = nylo_fft; j <= nyhi_fft; j++)
-      for (i = nxlo_fft; i <= nxhi_fft; i++) {
-        work4[n] = fkx[i]*fkz[k]*(work1[n+1]*fkx[i] + work2[n+1]*fky[j] + work3[n+1]*fkz[k]);
-        work4[n+1] = -fkx[i]*fkz[k]*(work1[n]*fkx[i] + work2[n]*fky[j] + work3[n]*fkz[k]);
-        n += 2;
-      }
-
-  fft2->compute(work4,work4,-1);
-
-  n = 0;
-  for (k = nzlo_in; k <= nzhi_in; k++)
-    for (j = nylo_in; j <= nyhi_in; j++)
-      for (i = nxlo_in; i <= nxhi_in; i++) {
-        vdxz_brick_dipole[k][j][i] = work4[n];
-        n += 2;
-      }
-
-  // Vyz
-
-  n = 0;
-  for (k = nzlo_fft; k <= nzhi_fft; k++)
-    for (j = nylo_fft; j <= nyhi_fft; j++)
-      for (i = nxlo_fft; i <= nxhi_fft; i++) {
-        work4[n] = fky[j]*fkz[k]*(work1[n+1]*fkx[i] + work2[n+1]*fky[j] + work3[n+1]*fkz[k]);
-        work4[n+1] = -fky[j]*fkz[k]*(work1[n]*fkx[i] + work2[n]*fky[j] + work3[n]*fkz[k]);
-        n += 2;
-      }
-
-  fft2->compute(work4,work4,-1);
-
-  n = 0;
-  for (k = nzlo_in; k <= nzhi_in; k++)
-    for (j = nylo_in; j <= nyhi_in; j++)
-      for (i = nxlo_in; i <= nxhi_in; i++) {
-        vdyz_brick_dipole[k][j][i] = work4[n];
-        n += 2;
-      }
-}
-
 /* ----------------------------------------------------------------------
    interpolate from grid to get electric field & force on my particles for ik
 ------------------------------------------------------------------------- */
@@ -808,7 +249,7 @@ void PPPMDielectric::fieldforce_ik()
   double *q = atom->q;
   double **x = atom->x;
   double **f = atom->f;
-  double *eps = avec->epsilon;
+  double *eps = atom->epsilon;
   int nlocal = atom->nlocal;
 
   for (i = 0; i < nlocal; i++) {
@@ -849,91 +290,13 @@ void PPPMDielectric::fieldforce_ik()
     efield[i][1] = efactor*eky;
     efield[i][2] = efactor*ekz;
 
-    const double qfactor = qqrd2e * scale * q[i] * eps[i];
+    const double qfactor = qqrd2e * efactor * q[i];
     f[i][0] += qfactor*ekx;
     f[i][1] += qfactor*eky;
     if (slabflag != 2) f[i][2] += qfactor*ekz;
   }
 }
 
-
-/* ----------------------------------------------------------------------
-   interpolate from grid to get electric field & force on my particles for ik
-------------------------------------------------------------------------- */
-
-void PPPMDielectric::fieldforce_ik_dipole()
-{
-  int i,l,m,n,nx,ny,nz,mx,my,mz;
-  FFT_SCALAR dx,dy,dz,u;
-  FFT_SCALAR x0,y0,z0;
-  FFT_SCALAR ex,ey,ez;
-  FFT_SCALAR vxx,vyy,vzz,vxy,vxz,vyz;
-
-  // loop over my charges, interpolate electric field from nearby grid points
-  // (nx,ny,nz) = global coords of grid pt to "lower left" of charge
-  // (dx,dy,dz) = distance to "lower left" grid pt
-  // (mx,my,mz) = global coords of moving stencil pt
-
-
-  double **mu = atom->mu;
-  double **x = atom->x;
-  double **f = atom->f;
-  double **t = atom->torque;
-
-  int nlocal = atom->nlocal;
-
-  for (i = 0; i < nlocal; i++) {
-    nx = part2grid[i][0];
-    ny = part2grid[i][1];
-    nz = part2grid[i][2];
-    dx = nx+shiftone - (x[i][0]-boxlo[0])*delxinv;
-    dy = ny+shiftone - (x[i][1]-boxlo[1])*delyinv;
-    dz = nz+shiftone - (x[i][2]-boxlo[2])*delzinv;
-
-    compute_rho1d(dx,dy,dz);
-
-    u = ex = ey = ez = ZEROF;
-    vxx = vyy = vzz = vxy = vxz = vyz = ZEROF;
-    for (n = nlower; n <= nupper; n++) {
-      mz = n+nz;
-      z0 = rho1d[2][n];
-      for (m = nlower; m <= nupper; m++) {
-        my = m+ny;
-        y0 = z0*rho1d[1][m];
-        for (l = nlower; l <= nupper; l++) {
-          mx = l+nx;
-          x0 = y0*rho1d[0][l];
-          u += x0*u_brick_dipole[mz][my][mx];
-          ex -= x0*ux_brick_dipole[mz][my][mx];
-          ey -= x0*uy_brick_dipole[mz][my][mx];
-          ez -= x0*uz_brick_dipole[mz][my][mx];
-          vxx -= x0*vdxx_brick_dipole[mz][my][mx];
-          vyy -= x0*vdyy_brick_dipole[mz][my][mx];
-          vzz -= x0*vdzz_brick_dipole[mz][my][mx];
-          vxy -= x0*vdxy_brick_dipole[mz][my][mx];
-          vxz -= x0*vdxz_brick_dipole[mz][my][mx];
-          vyz -= x0*vdyz_brick_dipole[mz][my][mx];
-        }
-      }
-    }
-
-    // electrical potential due to dipoles
-
-    if (potflag) phi[i] = u;
-
-    // convert E-field to torque
-
-    const double mufactor = qqrd2e * scale;
-    f[i][0] += mufactor*(vxx*mu[i][0] + vxy*mu[i][1] + vxz*mu[i][2]);
-    f[i][1] += mufactor*(vxy*mu[i][0] + vyy*mu[i][1] + vyz*mu[i][2]);
-    f[i][2] += mufactor*(vxz*mu[i][0] + vyz*mu[i][1] + vzz*mu[i][2]);
-
-    t[i][0] += mufactor*(mu[i][1]*ez - mu[i][2]*ey);
-    t[i][1] += mufactor*(mu[i][2]*ex - mu[i][0]*ez);
-    t[i][2] += mufactor*(mu[i][0]*ey - mu[i][1]*ex);
-  }
-}
-
 /* ----------------------------------------------------------------------
    interpolate from grid to get electric field & force on my particles for ad
 ------------------------------------------------------------------------- */
@@ -966,7 +329,7 @@ void PPPMDielectric::fieldforce_ad()
   double *q = atom->q;
   double **x = atom->x;
   double **f = atom->f;
-  double *eps = avec->epsilon;
+  double *eps = atom->epsilon;
   int nlocal = atom->nlocal;
 
   for (i = 0; i < nlocal; i++) {
@@ -1036,210 +399,6 @@ void PPPMDielectric::fieldforce_ad()
   }
 }
 
-/* ----------------------------------------------------------------------
-   pack own values to buf to send to another proc
-------------------------------------------------------------------------- */
-
-void PPPMDielectric::pack_forward(int flag, FFT_SCALAR *buf, int nlist, int *list)
-{
-  int n = 0;
-
-  if (flag == FORWARD_IK) {
-    FFT_SCALAR *xsrc = &vdx_brick[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *ysrc = &vdy_brick[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *zsrc = &vdz_brick[nzlo_out][nylo_out][nxlo_out];
-    for (int i = 0; i < nlist; i++) {
-      buf[n++] = xsrc[list[i]];
-      buf[n++] = ysrc[list[i]];
-      buf[n++] = zsrc[list[i]];
-    }
-  } else if (flag == FORWARD_MU) {
-    FFT_SCALAR *src_ux = &ux_brick_dipole[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *src_uy = &uy_brick_dipole[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *src_uz = &uz_brick_dipole[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *src_vxx = &vdxx_brick_dipole[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *src_vyy = &vdyy_brick_dipole[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *src_vzz = &vdzz_brick_dipole[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *src_vxy = &vdxy_brick_dipole[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *src_vxz = &vdxz_brick_dipole[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *src_vyz = &vdyz_brick_dipole[nzlo_out][nylo_out][nxlo_out];
-    for (int i = 0; i < nlist; i++) {
-      buf[n++] = src_ux[list[i]];
-      buf[n++] = src_uy[list[i]];
-      buf[n++] = src_uz[list[i]];
-      buf[n++] = src_vxx[list[i]];
-      buf[n++] = src_vyy[list[i]];
-      buf[n++] = src_vzz[list[i]];
-      buf[n++] = src_vxy[list[i]];
-      buf[n++] = src_vxz[list[i]];
-      buf[n++] = src_vyz[list[i]];
-    }
-  } else if (flag == FORWARD_AD) {
-    FFT_SCALAR *src = &u_brick[nzlo_out][nylo_out][nxlo_out];
-    for (int i = 0; i < nlist; i++)
-      buf[i] = src[list[i]];
-  } else if (flag == FORWARD_IK_PERATOM) {
-    FFT_SCALAR *esrc = &u_brick[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *v0src = &v0_brick[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *v1src = &v1_brick[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *v2src = &v2_brick[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *v3src = &v3_brick[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *v4src = &v4_brick[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *v5src = &v5_brick[nzlo_out][nylo_out][nxlo_out];
-    for (int i = 0; i < nlist; i++) {
-      if (eflag_atom) buf[n++] = esrc[list[i]];
-      if (vflag_atom) {
-        buf[n++] = v0src[list[i]];
-        buf[n++] = v1src[list[i]];
-        buf[n++] = v2src[list[i]];
-        buf[n++] = v3src[list[i]];
-        buf[n++] = v4src[list[i]];
-        buf[n++] = v5src[list[i]];
-      }
-    }
-  } else if (flag == FORWARD_AD_PERATOM) {
-    FFT_SCALAR *v0src = &v0_brick[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *v1src = &v1_brick[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *v2src = &v2_brick[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *v3src = &v3_brick[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *v4src = &v4_brick[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *v5src = &v5_brick[nzlo_out][nylo_out][nxlo_out];
-    for (int i = 0; i < nlist; i++) {
-      buf[n++] = v0src[list[i]];
-      buf[n++] = v1src[list[i]];
-      buf[n++] = v2src[list[i]];
-      buf[n++] = v3src[list[i]];
-      buf[n++] = v4src[list[i]];
-      buf[n++] = v5src[list[i]];
-    }
-  }
-}
-
-/* ----------------------------------------------------------------------
-   unpack another proc's own values from buf and set own ghost values
-------------------------------------------------------------------------- */
-
-void PPPMDielectric::unpack_forward(int flag, FFT_SCALAR *buf, int nlist, int *list)
-{
-  int n = 0;
-
-  if (flag == FORWARD_IK) {
-    FFT_SCALAR *xdest = &vdx_brick[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *ydest = &vdy_brick[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *zdest = &vdz_brick[nzlo_out][nylo_out][nxlo_out];
-    for (int i = 0; i < nlist; i++) {
-      xdest[list[i]] = buf[n++];
-      ydest[list[i]] = buf[n++];
-      zdest[list[i]] = buf[n++];
-    }
-  } else if (flag == FORWARD_MU) {
-    FFT_SCALAR *dest_ux = &ux_brick_dipole[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *dest_uy = &uy_brick_dipole[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *dest_uz = &uz_brick_dipole[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *dest_vxx = &vdxx_brick_dipole[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *dest_vyy = &vdyy_brick_dipole[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *dest_vzz = &vdzz_brick_dipole[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *dest_vxy = &vdxy_brick_dipole[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *dest_vxz = &vdxz_brick_dipole[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *dest_vyz = &vdyz_brick_dipole[nzlo_out][nylo_out][nxlo_out];
-    for (int i = 0; i < nlist; i++) {
-      dest_ux[list[i]] = buf[n++];
-      dest_uy[list[i]] = buf[n++];
-      dest_uz[list[i]] = buf[n++];
-      dest_vxx[list[i]] = buf[n++];
-      dest_vyy[list[i]] = buf[n++];
-      dest_vzz[list[i]] = buf[n++];
-      dest_vxy[list[i]] = buf[n++];
-      dest_vxz[list[i]] = buf[n++];
-      dest_vyz[list[i]] = buf[n++];
-    }
-  } else if (flag == FORWARD_AD) {
-    FFT_SCALAR *dest = &u_brick[nzlo_out][nylo_out][nxlo_out];
-    for (int i = 0; i < nlist; i++)
-      dest[list[i]] = buf[i];
-  } else if (flag == FORWARD_IK_PERATOM) {
-    FFT_SCALAR *esrc = &u_brick[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *v0src = &v0_brick[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *v1src = &v1_brick[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *v2src = &v2_brick[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *v3src = &v3_brick[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *v4src = &v4_brick[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *v5src = &v5_brick[nzlo_out][nylo_out][nxlo_out];
-    for (int i = 0; i < nlist; i++) {
-      if (eflag_atom) esrc[list[i]] = buf[n++];
-      if (vflag_atom) {
-        v0src[list[i]] = buf[n++];
-        v1src[list[i]] = buf[n++];
-        v2src[list[i]] = buf[n++];
-        v3src[list[i]] = buf[n++];
-        v4src[list[i]] = buf[n++];
-        v5src[list[i]] = buf[n++];
-      }
-    }
-  } else if (flag == FORWARD_AD_PERATOM) {
-    FFT_SCALAR *v0src = &v0_brick[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *v1src = &v1_brick[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *v2src = &v2_brick[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *v3src = &v3_brick[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *v4src = &v4_brick[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *v5src = &v5_brick[nzlo_out][nylo_out][nxlo_out];
-    for (int i = 0; i < nlist; i++) {
-      v0src[list[i]] = buf[n++];
-      v1src[list[i]] = buf[n++];
-      v2src[list[i]] = buf[n++];
-      v3src[list[i]] = buf[n++];
-      v4src[list[i]] = buf[n++];
-      v5src[list[i]] = buf[n++];
-    }
-  }
-}
-
-/* ----------------------------------------------------------------------
-   pack ghost values into buf to send to another proc
-------------------------------------------------------------------------- */
-
-void PPPMDielectric::pack_reverse(int flag, FFT_SCALAR *buf, int nlist, int *list)
-{
-  int n = 0;
-  if (flag == REVERSE_RHO) {
-    FFT_SCALAR *src = &density_brick[nzlo_out][nylo_out][nxlo_out];
-    for (int i = 0; i < nlist; i++)
-      buf[n++] = src[list[i]];
-  } else if (flag == REVERSE_MU) {
-    FFT_SCALAR *src_mu0 = &densityx_brick_dipole[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *src_mu1 = &densityy_brick_dipole[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *src_mu2 = &densityz_brick_dipole[nzlo_out][nylo_out][nxlo_out];
-    for (int i = 0; i < nlist; i++) {
-      buf[n++] = src_mu0[list[i]];
-      buf[n++] = src_mu1[list[i]];
-      buf[n++] = src_mu2[list[i]];
-    }
-  }
-}
-
-/* ----------------------------------------------------------------------
-   unpack another proc's ghost values from buf and add to own values
-------------------------------------------------------------------------- */
-
-void PPPMDielectric::unpack_reverse(int flag, FFT_SCALAR *buf, int nlist, int *list)
-{
-  int n = 0;
-  if (flag == REVERSE_RHO) {
-    FFT_SCALAR *dest = &density_brick[nzlo_out][nylo_out][nxlo_out];
-    for (int i = 0; i < nlist; i++)
-      dest[list[i]] += buf[n++];
-  } else if (flag == REVERSE_MU) {
-    FFT_SCALAR *dest_mu0 = &densityx_brick_dipole[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *dest_mu1 = &densityy_brick_dipole[nzlo_out][nylo_out][nxlo_out];
-    FFT_SCALAR *dest_mu2 = &densityz_brick_dipole[nzlo_out][nylo_out][nxlo_out];
-    for (int i = 0; i < nlist; i++) {
-      dest_mu0[list[i]] += buf[n++];
-      dest_mu1[list[i]] += buf[n++];
-      dest_mu2[list[i]] += buf[n++];
-    }
-  }
-}
-
 /* ----------------------------------------------------------------------
    Slab-geometry correction term to dampen inter-slab interactions between
    periodically repeating slabs.  Yields good approximation to 2D Ewald if
@@ -1254,18 +413,14 @@ void PPPMDielectric::slabcorr()
 
   double *q = atom->q;
   double **x = atom->x;
-  double *eps = avec->epsilon;
+  double *eps = atom->epsilon;
+  
   double zprd = domain->zprd;
   int nlocal = atom->nlocal;
 
   double dipole = 0.0;
   for (int i = 0; i < nlocal; i++) dipole += q[i]*x[i][2];
 
-  if (mu_flag) {
-    double **mu = atom->mu;
-    for (int i = 0; i < nlocal; i++) dipole += mu[i][2];
-  }
-
   // sum local contributions to get global dipole moment
 
   double dipole_all;
@@ -1277,10 +432,6 @@ void PPPMDielectric::slabcorr()
   double dipole_r2 = 0.0;
   if (eflag_atom || fabs(qsum) > SMALL) {
 
-    if (mu_flag)
-      error->all(FLERR,"Cannot (yet) use kspace slab correction with "
-        "long-range dipoles and non-neutral systems or per-atom energy");
-
     for (int i = 0; i < nlocal; i++)
       dipole_r2 += q[i]*x[i][2]*x[i][2];
 
@@ -1317,17 +468,6 @@ void PPPMDielectric::slabcorr()
     f[i][2] += ffact * eps[i]*q[i]*(dipole_all - qsum*x[i][2]);
     efield[i][2] += ffact * eps[i]*(dipole_all - qsum*x[i][2]);
   }
-
-  // add on torque corrections
-
-  if (mu_flag && atom->torque) {
-    double **mu = atom->mu;
-    double **torque = atom->torque;
-    for (int i = 0; i < nlocal; i++) {
-      torque[i][0] += ffact * dipole_all * mu[i][1];
-      torque[i][1] += -ffact * dipole_all * mu[i][0];
-    }
-  }
 }
 
 /* ----------------------------------------------------------------------
@@ -1338,7 +478,7 @@ void PPPMDielectric::slabcorr()
 void PPPMDielectric::qsum_qsq()
 {
   const double * const q = atom->q;
-  const double * const eps = avec->epsilon;
+  const double * const eps = atom->epsilon;
   const int nlocal = atom->nlocal;
   double qsum_local(0.0), qsqsum_local(0.0);
 
@@ -1356,84 +496,3 @@ void PPPMDielectric::qsum_qsq()
 
   q2 = qsqsum * force->qqrd2e;
 }
-
-/* ----------------------------------------------------------------------
-   memory usage of local arrays
-------------------------------------------------------------------------- */
-
-double PPPMDielectric::memory_usage()
-{
-  double bytes = nmax*3 * sizeof(double);
-  int nbrick = (nxhi_out-nxlo_out+1) * (nyhi_out-nylo_out+1) *
-    (nzhi_out-nzlo_out+1);
-  if (differentiation_flag == 1) {
-    bytes += 2 * nbrick * sizeof(FFT_SCALAR);
-  } else {
-    bytes += 4 * nbrick * sizeof(FFT_SCALAR);
-  }
-  if (triclinic) bytes += 3 * nfft_both * sizeof(double);
-  bytes += 6 * nfft_both * sizeof(double);
-  bytes += nfft_both * sizeof(double);
-  bytes += nfft_both*5 * sizeof(FFT_SCALAR);
-
-  if (mu_flag) {
-    bytes += 3 * nbrick * sizeof(FFT_SCALAR);
-    //work?
-  }
-
-  if (peratom_allocate_flag)
-    bytes += 6 * nbrick * sizeof(FFT_SCALAR);
-
-  if (group_allocate_flag) {
-    bytes += 2 * nbrick * sizeof(FFT_SCALAR);
-    bytes += 2 * nfft_both * sizeof(FFT_SCALAR);;
-  }
-
-  if (cg) bytes += cg->memory_usage();
-
-  if (mu_flag)
-    bytes += cg_mu->memory_usage();
-
-  return bytes;
-}
-
-
-/* ----------------------------------------------------------------------
-   compute qsum,qsqsum,q2 and give error/warning if not charge neutral
-   called initially, when particle count changes, when charges are changed
-------------------------------------------------------------------------- */
-
-void PPPMDielectric::musum_musq()
-{
-  double** mu = atom->mu;
-  const int nlocal = atom->nlocal;
-  double musum_local(0.0), musqsum_local(0.0);
-
-  for (int i = 0; i < nlocal; i++) {
-    musum_local += mu[i][0] + mu[i][1] + mu[i][2];
-    musqsum_local += mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] + mu[i][2]*mu[i][2];
-  }
-
-  MPI_Allreduce(&musum_local,&musum,1,MPI_DOUBLE,MPI_SUM,world);
-  MPI_Allreduce(&musqsum_local,&musqsum,1,MPI_DOUBLE,MPI_SUM,world);
-
-  /*
-  if ((qsqsum == 0.0) && (comm->me == 0) && warn_nocharge) {
-    error->warning(FLERR,"Using kspace solver on system with no charge");
-    warn_nocharge = 0;
-  }*/
-
-  mu2 = musqsum * force->qqrd2e;
-
-  // not yet sure of the correction needed for non-neutral systems
-  // so issue warning or error
-  /*
-  if (fabs(qsum) > SMALL) {
-    char str[128];
-    sprintf(str,"System is not charge neutral, net charge = %g",qsum);
-    if (!warn_nonneutral) error->all(FLERR,str);
-    if (warn_nonneutral == 1 && comm->me == 0) error->warning(FLERR,str);
-    warn_nonneutral = 2;
-  }*/
-}
-
diff --git a/src/USER-DIELECTRIC/pppm_dielectric.h b/src/USER-DIELECTRIC/pppm_dielectric.h
index e30d03cd9b..b7803011f8 100644
--- a/src/USER-DIELECTRIC/pppm_dielectric.h
+++ b/src/USER-DIELECTRIC/pppm_dielectric.h
@@ -28,10 +28,7 @@ class PPPMDielectric : public PPPM {
  public:
   PPPMDielectric(class LAMMPS *);
   virtual ~PPPMDielectric();
-  virtual void init();
-  void setup_grid();
   virtual void compute(int, int);
-  virtual double memory_usage();
 
   double** efield;
   double* phi;
@@ -40,40 +37,11 @@ class PPPMDielectric : public PPPM {
   void qsum_qsq();
 
  protected:
-  FFT_SCALAR ***densityx_brick_dipole,***densityy_brick_dipole,***densityz_brick_dipole;
-  FFT_SCALAR ***vdxx_brick_dipole,***vdyy_brick_dipole,***vdzz_brick_dipole;
-  FFT_SCALAR ***vdxy_brick_dipole,***vdxz_brick_dipole,***vdyz_brick_dipole;
-  FFT_SCALAR ***u_brick_dipole;
-  FFT_SCALAR ***ux_brick_dipole,***uy_brick_dipole,***uz_brick_dipole;
-  FFT_SCALAR *densityx_fft_dipole,*densityy_fft_dipole,*densityz_fft_dipole;
-  FFT_SCALAR *work3,*work4;
-
-  class GridComm *cg_mu;
-
-  virtual void allocate();
-  virtual void deallocate();
   void slabcorr();
 
   void fieldforce_ik();
   void fieldforce_ad();
 
-  // grid communication
-
-  virtual void pack_forward(int, FFT_SCALAR *, int, int *);
-  virtual void unpack_forward(int, FFT_SCALAR *, int, int *);
-  virtual void pack_reverse(int, FFT_SCALAR *, int, int *);
-  virtual void unpack_reverse(int, FFT_SCALAR *, int, int *);
-
-  // dipole
-
-  int mu_flag;
-  double musqsum,musum,mu2;
-  void make_rho_dipole();
-  void brick2fft_dipole();
-  void poisson_ik_dipole();
-  void fieldforce_ik_dipole();
-  void musum_musq();
-
   class AtomVecDielectric* avec;
 
 };
diff --git a/src/atom.cpp b/src/atom.cpp
index ecb82993ce..63e90692e0 100644
--- a/src/atom.cpp
+++ b/src/atom.cpp
@@ -197,6 +197,10 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
   rho = drho = esph = desph = cv = nullptr;
   vest = nullptr;
 
+  // USER-DIELECTRIC package
+
+  area = ed = em = epsilon = curvature = q_unscaled = nullptr;
+
   // end of customization section
   // --------------------------------------------------------------------
 
@@ -509,6 +513,15 @@ void Atom::peratom_create()
   add_peratom("eff_plastic_strain_rate",&eff_plastic_strain_rate,DOUBLE,0);
   add_peratom("damage",&damage,DOUBLE,0);
 
+  // USER-DIELECTRIC package
+
+  add_peratom("area",&area,DOUBLE,0);
+  add_peratom("ed",&ed,DOUBLE,0);
+  add_peratom("em",&em,DOUBLE,0);
+  add_peratom("epsilon",&epsilon,DOUBLE,0);
+  add_peratom("curvature",&curvature,DOUBLE,0);
+  add_peratom("q_unscaled",&curvature,DOUBLE,0);
+
   // end of customization section
   // --------------------------------------------------------------------
 }
diff --git a/src/atom.h b/src/atom.h
index bc69d3b27a..bd1594480c 100644
--- a/src/atom.h
+++ b/src/atom.h
@@ -157,6 +157,10 @@ class Atom : protected Pointers {
   double *rho,*drho,*esph,*desph,*cv;
   double **vest;
 
+  // USER-DIELECTRIC package
+
+  double *area,*ed,*em,*epsilon,*curvature,*q_unscaled;
+
   // end of customization section
   // --------------------------------------------------------------------