diff --git a/src/DIELECTRIC/atom_vec_dielectric.cpp b/src/DIELECTRIC/atom_vec_dielectric.cpp
index ca89258bb1..3b25ad4e17 100644
--- a/src/DIELECTRIC/atom_vec_dielectric.cpp
+++ b/src/DIELECTRIC/atom_vec_dielectric.cpp
@@ -60,28 +60,28 @@ AtomVecDielectric::AtomVecDielectric(LAMMPS *_lmp) : AtomVec(_lmp)
     "angle_atom1", "angle_atom2", "angle_atom3", "num_dihedral", "dihedral_type", "dihedral_atom1",
     "dihedral_atom2", "dihedral_atom3", "dihedral_atom4", "num_improper", "improper_type",
     "improper_atom1", "improper_atom2", "improper_atom3", "improper_atom4", "nspecial", "special",
-    "mu", "area", "ed", "em", "epsilon", "curvature", "q_unscaled"};
+    "mu", "area", "ed", "em", "epsilon", "curvature", "q_scaled"};
   fields_copy = {"q", "molecule", "num_bond", "bond_type", "bond_atom", "num_angle", "angle_type",
     "angle_atom1", "angle_atom2", "angle_atom3", "num_dihedral", "dihedral_type", "dihedral_atom1",
     "dihedral_atom2", "dihedral_atom3", "dihedral_atom4", "num_improper", "improper_type",
     "improper_atom1", "improper_atom2", "improper_atom3", "improper_atom4", "nspecial", "special",
-    "mu", "area", "ed", "em", "epsilon", "curvature", "q_unscaled"};
-  fields_comm = {"q", "mu", "area", "ed", "em", "epsilon", "curvature", "q_unscaled"};
-  fields_border = {"q", "molecule", "mu", "area", "ed", "em", "epsilon", "curvature", "q_unscaled"};
+    "mu", "area", "ed", "em", "epsilon", "curvature", "q_scaled"};
+  fields_comm = {"q", "mu", "area", "ed", "em", "epsilon", "curvature", "q_scaled"};
+  fields_border = {"q", "molecule", "mu", "area", "ed", "em", "epsilon", "curvature", "q_scaled"};
   fields_border_vel = {"q", "molecule", "mu", "area", "ed", "em", "epsilon", "curvature",
-    "q_unscaled"};
+    "q_scaled"};
   fields_exchange = {"q", "molecule", "num_bond", "bond_type", "bond_atom", "num_angle",
     "angle_type", "angle_atom1", "angle_atom2", "angle_atom3", "num_dihedral", "dihedral_type",
     "dihedral_atom1", "dihedral_atom2", "dihedral_atom3", "dihedral_atom4", "num_improper",
     "improper_type", "improper_atom1", "improper_atom2", "improper_atom3", "improper_atom4",
-    "nspecial", "special", "mu", "area", "ed", "em", "epsilon", "curvature", "q_unscaled"};
+    "nspecial", "special", "mu", "area", "ed", "em", "epsilon", "curvature", "q_scaled"};
   fields_restart = {"q", "molecule", "num_bond", "bond_type", "bond_atom", "num_angle",
     "angle_type", "angle_atom1", "angle_atom2", "angle_atom3", "num_dihedral", "dihedral_type",
     "dihedral_atom1", "dihedral_atom2", "dihedral_atom3", "dihedral_atom4", "num_improper",
     "improper_type", "improper_atom1", "improper_atom2", "improper_atom3", "improper_atom4",
-    "mu", "area", "ed", "em", "epsilon", "curvature", "q_unscaled"};
+    "mu", "area", "ed", "em", "epsilon", "curvature", "q_scaled"};
   fields_create = {"q", "molecule", "num_bond", "num_angle", "num_dihedral", "num_improper",
-    "nspecial", "mu", "area", "ed", "em", "epsilon", "curvature", "q_unscaled"};
+    "nspecial", "mu", "area", "ed", "em", "epsilon", "curvature", "q_scaled"};
   fields_data_atom = {"id", "molecule", "type", "q", "x", "mu3", "area", "ed", "em", "epsilon",
     "curvature"};
   fields_data_vel = {"id", "v"};
@@ -151,7 +151,7 @@ void AtomVecDielectric::grow_pointers()
   em = atom->em;
   epsilon = atom->epsilon;
   curvature = atom->curvature;
-  q_unscaled = atom->q_unscaled;
+  q_scaled = atom->q_scaled;
 }
 
 /* ----------------------------------------------------------------------
@@ -181,31 +181,12 @@ void AtomVecDielectric::data_atom_post(int ilocal)
   nspecial[ilocal][2] = 0;
 
   double *q = atom->q;
-  q_unscaled[ilocal] = q[ilocal];
-  q[ilocal] /= epsilon[ilocal];
+  q_scaled[ilocal] = q[ilocal]/epsilon[ilocal];
 
   double *mu_one = mu[ilocal];
   mu_one[3] = sqrt(mu_one[0] * mu_one[0] + mu_one[1] * mu_one[1] + mu_one[2] * mu_one[2]);
 }
 
-/* ----------------------------------------------------------------------
-   restore original data for writing the data file
-------------------------------------------------------------------------- */
-
-void AtomVecDielectric::pack_data_pre(int ilocal)
-{
-  atom->q[ilocal] = q_unscaled[ilocal];
-}
-
-/* ----------------------------------------------------------------------
-   undo restore and get back to post read data state
-------------------------------------------------------------------------- */
-
-void AtomVecDielectric::pack_data_post(int ilocal)
-{
-  atom->q[ilocal] /= epsilon[ilocal];
-}
-
 /* ----------------------------------------------------------------------
    initialize other atom quantities after AtomVec::unpack_restart()
 ------------------------------------------------------------------------- */
@@ -229,7 +210,7 @@ int AtomVecDielectric::property_atom(const std::string &name)
   if (name == "em") return 2;
   if (name == "epsilon") return 3;
   if (name == "curvature") return 4;
-  if (name == "q_unscaled") return 5;
+  if (name == "q_scaled") return 5;
   return -1;
 }
 
@@ -287,7 +268,7 @@ void AtomVecDielectric::pack_property_atom(int index, double *buf, int nvalues,
   } else if (index == 5) {
     for (int i = 0; i < nlocal; i++) {
       if (mask[i] & groupbit)
-        buf[n] = q_unscaled[i];
+        buf[n] = q_scaled[i];
       else
         buf[n] = 0.0;
       n += nvalues;
diff --git a/src/DIELECTRIC/atom_vec_dielectric.h b/src/DIELECTRIC/atom_vec_dielectric.h
index 34fe11b249..28bf7abb33 100644
--- a/src/DIELECTRIC/atom_vec_dielectric.h
+++ b/src/DIELECTRIC/atom_vec_dielectric.h
@@ -38,8 +38,6 @@ class AtomVecDielectric : virtual public AtomVec {
   void unpack_restart_init(int) override;
   int property_atom(const std::string &) override;
   void pack_property_atom(int, double *, int, int) override;
-  void pack_data_pre(int) override;
-  void pack_data_post(int) override;
 
  protected:
   int *num_bond, *num_angle, *num_dihedral, *num_improper;
@@ -49,7 +47,7 @@ class AtomVecDielectric : virtual public AtomVec {
   int bond_per_atom, angle_per_atom, dihedral_per_atom, improper_per_atom;
 
   double **mu;
-  double *area, *ed, *em, *epsilon, *curvature, *q_unscaled;
+  double *area, *ed, *em, *epsilon, *curvature, *q_scaled;
 };
 
 }    // namespace LAMMPS_NS
diff --git a/src/DIELECTRIC/compute_efield_atom.cpp b/src/DIELECTRIC/compute_efield_atom.cpp
index a6f2ff6d2c..5e847a2a28 100644
--- a/src/DIELECTRIC/compute_efield_atom.cpp
+++ b/src/DIELECTRIC/compute_efield_atom.cpp
@@ -27,6 +27,7 @@
 #include "pair_coul_cut_dielectric.h"
 #include "pair_coul_long_dielectric.h"
 #include "pair_lj_cut_coul_cut_dielectric.h"
+#include "pair_lj_cut_coul_debye_dielectric.h"
 #include "pair_lj_cut_coul_long_dielectric.h"
 #include "pair_lj_cut_coul_msm_dielectric.h"
 #include "pppm_dielectric.h"
@@ -68,8 +69,6 @@ ComputeEfieldAtom::ComputeEfieldAtom(LAMMPS *_lmp, int narg, char **arg) :
   }
 
   nmax = 0;
-
-  comm_reverse = 1;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -100,6 +99,10 @@ void ComputeEfieldAtom::setup()
   else if (strcmp(force->pair_style, "lj/cut/coul/cut/dielectric") == 0)
     efield_pair = (dynamic_cast<PairLJCutCoulCutDielectric *>(force->pair))->efield;
   else if (strcmp(force->pair_style, "lj/cut/coul/cut/dielectric/omp") == 0)
+    efield_pair = (dynamic_cast<PairLJCutCoulDebyeDielectric *>(force->pair))->efield;
+  else if (strcmp(force->pair_style, "lj/cut/coul/debye/dielectric") == 0)
+    efield_pair = (dynamic_cast<PairLJCutCoulDebyeDielectric *>(force->pair))->efield;
+  else if (strcmp(force->pair_style, "lj/cut/coul/debye/dielectric/omp") == 0)
     efield_pair = (dynamic_cast<PairLJCutCoulCutDielectric *>(force->pair))->efield;
   else if (strcmp(force->pair_style, "coul/long/dielectric") == 0)
     efield_pair = (dynamic_cast<PairCoulLongDielectric *>(force->pair))->efield;
diff --git a/src/DIELECTRIC/fix_polarize_bem_gmres.cpp b/src/DIELECTRIC/fix_polarize_bem_gmres.cpp
index 6f76bdee6b..e31337d0fd 100644
--- a/src/DIELECTRIC/fix_polarize_bem_gmres.cpp
+++ b/src/DIELECTRIC/fix_polarize_bem_gmres.cpp
@@ -41,6 +41,7 @@
 #include "comm.h"
 #include "error.h"
 #include "force.h"
+#include "group.h"
 #include "kspace.h"
 #include "math_const.h"
 #include "memory.h"
@@ -234,7 +235,7 @@ void FixPolarizeBEMGMRES::init()
   }
 
   if (comm->me == 0)
-    utils::logmesg(lmp, "GMRES solver for {} induced charges using maximum {} q-vectors\n",
+    utils::logmesg(lmp, "BEM/GMRES solver for {} induced charges using maximum {} q-vectors\n",
                    num_induced_charges, mr);
 }
 
@@ -265,6 +266,8 @@ void FixPolarizeBEMGMRES::setup(int /*vflag*/)
   else
     error->all(FLERR, "Pair style not compatible with fix polarize/bem/gmres");
 
+  // check if kspace is used for force computation
+
   if (force->kspace) {
     kspaceflag = 1;
     if (strcmp(force->kspace_style, "pppm/dielectric") == 0)
@@ -340,32 +343,35 @@ void FixPolarizeBEMGMRES::pre_force(int)
 
 void FixPolarizeBEMGMRES::compute_induced_charges()
 {
+  double *q_scaled = atom->q_scaled;
   double *q = atom->q;
-  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;
-  int eflag = 0;
+  int eflag = 1;
   int vflag = 0;
 
   // compute the right hand side (vector b) of Eq. (40) according to Eq. (42)
   // keep the scaled real charges intact here to compute efield for the right hand side (b)
   //   and backup all the charges
-  // for induced charges q_real stores the free surface charge
+  // for induced charges q stores the free surface charge
   // set the induced charges to be zero to compute the right hand side (b)
   // the current value can be accessed via induced_charges[induced_charge_idx[i]]
 
   for (int i = 0; i < nlocal; i++) {
-    q_backup[i] = q[i];
-    if (induced_charge_idx[i] >= 0) q[i] = 0;
+    q_backup[i] = q_scaled[i];
+    if (induced_charge_idx[i] >= 0) q_scaled[i] = 0;
   }
 
+  // communicate q_scaled between the neighboring procs
+
   comm->forward_comm(this);
 
-  // note here q[i] are the bound charges including area
+  // note here q_scaled[i] are the bound charges including area
   // so that kspace solver can be used directly with the charge values
   // for the moment, require that newton off and full neighbor list for pair
   // Note that in the definition of the electrical fields in Equations (41) and (53)
@@ -396,7 +402,7 @@ void FixPolarizeBEMGMRES::compute_induced_charges()
       Ez += efield_kspace[i][2];
     }
     double ndotE = epsilon0e2q * (Ex * norm[i][0] + Ey * norm[i][1] + Ez * norm[i][2]) / epsilon[i];
-    double sigma_f = q_real[i] / area[i];
+    double sigma_f = q[i] / area[i];
     buffer[idx] = (1 - em[i]) * sigma_f - ed[i] * ndotE / (4 * MY_PI);
   }
 
@@ -429,15 +435,41 @@ void FixPolarizeBEMGMRES::compute_induced_charges()
   for (int i = 0; i < nlocal; i++) {
     if (induced_charge_idx[i] >= 0) {
       int idx = induced_charge_idx[i];
-      q[i] = induced_charges[idx] * area[i] + q_real[i];
+      q_scaled[i] = induced_charges[idx] * area[i] + q[i];
     } else {
-      q[i] = q_backup[i];
+      q_scaled[i] = q_backup[i];
     }
   }
 
+  // communicate q_scaled to neighboring procs
+
   comm->forward_comm(this);
 
+  // compute the total induced charges of the interface particles
+  // for interface particles: set the charge to be the sum of unscaled (free) charges and induced charges
+
+  double tmp = 0;
+  for (int i = 0; i < nlocal; i++) {
+    if (!(mask[i] & groupbit)) continue;
+
+    double q_bound = q_scaled[i] - q[i];
+    tmp += q_bound;
+    q[i] = q_scaled[i];
+  }
+
   if (first) first = 0;
+
+  // ensure sum of all induced charges being zero
+
+  int ncount = group->count(igroup);
+  double sum = 0;
+  MPI_Allreduce(&tmp, &sum, 1, MPI_DOUBLE, MPI_SUM, world);
+  double qboundave = sum/(double)ncount;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (!(mask[i] & groupbit)) continue;
+    q[i] -=  qboundave;
+  }
 }
 
 /* ---------------------------------------------------------------------- */
@@ -609,14 +641,14 @@ void FixPolarizeBEMGMRES::gmres_solve(double *x, double *r)
 void FixPolarizeBEMGMRES::apply_operator(double *w, double *Aw, int /*n*/)
 {
   int i;
-  double *q = atom->q;
+  double *q_scaled = atom->q_scaled;
   double **norm = atom->mu;
   double *area = atom->area;
   double *ed = atom->ed;
   double *em = atom->em;
   double *epsilon = atom->epsilon;
   int nlocal = atom->nlocal;
-  int eflag = 0;
+  int eflag = 1;
   int vflag = 0;
 
   // set the induced charges to be w
@@ -626,10 +658,10 @@ void FixPolarizeBEMGMRES::apply_operator(double *w, double *Aw, int /*n*/)
 
   for (i = 0; i < nlocal; i++) {
     if (induced_charge_idx[i] < 0) {
-      q[i] = 0;
+      q_scaled[i] = 0;
     } else {
       int idx = induced_charge_idx[i];
-      q[i] = w[idx] * area[i];
+      q_scaled[i] = w[idx] * area[i];
     }
   }
 
@@ -677,15 +709,15 @@ void FixPolarizeBEMGMRES::apply_operator(double *w, double *Aw, int /*n*/)
 void FixPolarizeBEMGMRES::update_residual(double *w, double *r, int /*n*/)
 {
   int i;
+  double *q_scaled = atom->q_scaled;
   double *q = atom->q;
-  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 nlocal = atom->nlocal;
-  int eflag = 0;
+  int eflag = 1;
   int vflag = 0;
 
   // compute the Coulombic forces and electrical field E
@@ -695,13 +727,15 @@ void FixPolarizeBEMGMRES::update_residual(double *w, double *r, int /*n*/)
 
   for (i = 0; i < nlocal; i++) {
     if (induced_charge_idx[i] < 0) {
-      q[i] = q_backup[i];
+      q_scaled[i] = q_backup[i];
     } else {
       int idx = induced_charge_idx[i];
-      q[i] = w[idx] * area[i] + q_real[i];
+      q_scaled[i] = w[idx] * area[i] + q[i];
     }
   }
 
+  // communicate q_scaled between neighboring procs
+
   comm->forward_comm(this);
 
   force_clear();
@@ -733,7 +767,7 @@ void FixPolarizeBEMGMRES::update_residual(double *w, double *r, int /*n*/)
     }
     double ndotE = epsilon0e2q * (Ex * norm[i][0] + Ey * norm[i][1] + Ez * norm[i][2]) /
         epsilon[i] / (4 * MY_PI);
-    double sigma_f = q_real[i] / area[i];
+    double sigma_f = q[i] / area[i];
     buffer[idx] = (1 - em[i]) * sigma_f - em[i] * w[idx] - ed[i] * ndotE;
   }
 
@@ -895,7 +929,7 @@ int FixPolarizeBEMGMRES::pack_forward_comm(int n, int *list, double *buf, int /*
                                            int * /*pbc*/)
 {
   int m;
-  for (m = 0; m < n; m++) buf[m] = atom->q[list[m]];
+  for (m = 0; m < n; m++) buf[m] = atom->q_scaled[list[m]];
   return n;
 }
 
@@ -904,7 +938,7 @@ int FixPolarizeBEMGMRES::pack_forward_comm(int n, int *list, double *buf, int /*
 void FixPolarizeBEMGMRES::unpack_forward_comm(int n, int first, double *buf)
 {
   int i, m;
-  for (m = 0, i = first; m < n; m++, i++) atom->q[i] = buf[m];
+  for (m = 0, i = first; m < n; m++, i++) atom->q_scaled[i] = buf[m];
 }
 
 /* ----------------------------------------------------------------------
@@ -952,7 +986,7 @@ void FixPolarizeBEMGMRES::set_dielectric_params(double ediff, double emean, doub
   double *area = atom->area;
   double *ed = atom->ed;
   double *em = atom->em;
-  double *q_unscaled = atom->q_unscaled;
+  double *q = atom->q;
   double *epsilon = atom->epsilon;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
@@ -963,7 +997,7 @@ void FixPolarizeBEMGMRES::set_dielectric_params(double ediff, double emean, doub
       em[i] = emean;
       if (areai > 0) area[i] = areai;
       if (epsiloni > 0) epsilon[i] = epsiloni;
-      if (set_charge) q_unscaled[i] = qvalue;
+      if (set_charge) q[i] = qvalue;
     }
   }
 }
diff --git a/src/DIELECTRIC/fix_polarize_bem_icc.cpp b/src/DIELECTRIC/fix_polarize_bem_icc.cpp
index 2f9339d483..0d503d87a7 100644
--- a/src/DIELECTRIC/fix_polarize_bem_icc.cpp
+++ b/src/DIELECTRIC/fix_polarize_bem_icc.cpp
@@ -104,7 +104,11 @@ int FixPolarizeBEMICC::setmask()
 void FixPolarizeBEMICC::init()
 {
   int ncount = group->count(igroup);
-  if (comm->me == 0) utils::logmesg(lmp, "BEM/ICC solver for {} induced charges\n", ncount);
+  if (comm->me == 0) {
+    utils::logmesg(lmp, "BEM/ICC solver for {} induced charges\n", ncount);
+    utils::logmesg(lmp, " using pair style {}\n", force->pair_style);
+    if (force->kspace) utils::logmesg(lmp, " using kspace style {}\n", force->kspace_style);
+  }
 
   // initialize random induced charges with zero sum
 
@@ -241,8 +245,8 @@ void FixPolarizeBEMICC::pre_force(int)
 
 void FixPolarizeBEMICC::compute_induced_charges()
 {
+  double *q_scaled = atom->q_scaled;
   double *q = atom->q;
-  double *q_real = atom->q_unscaled;
   double **norm = atom->mu;
   double *area = atom->area;
   double *ed = atom->ed;
@@ -261,7 +265,8 @@ void FixPolarizeBEMICC::compute_induced_charges()
   //   q_real are read from the data file
   // Note that the electrical fields here are due to the rescaled real charges,
   //   and also multiplied by epsilon[i]
-  // Let's choose that epsilon[i] = em[i] for the interface particles
+  // for the interface particles assume that epsilon[i] = em[i]
+  // NOTE: 13Dec2022 pair and kspace with dielectric suffix operate on q_scaled
 
   force_clear();
   force->pair->compute(eflag, vflag);
@@ -283,15 +288,17 @@ void FixPolarizeBEMICC::compute_induced_charges()
     // divide (Ex,Ey,Ez) by epsilon[i] here
     double ndotE = epsilon0e2q * (Ex * norm[i][0] + Ey * norm[i][1] + Ez * norm[i][2]) /
         epsilon[i] / (2 * MY_PI);
-    double q_free = q_real[i];
+    double q_free = q[i];
     double q_bound = 0;
     q_bound = (1.0 / em[i] - 1) * q_free - (ed[i] / (2 * em[i])) * ndotE * area[i];
-    q[i] = q_free + q_bound;
+    q_scaled[i] = q_free + q_bound;
   }
 
+  // communicate q_scaled between neighboring procs
+
   comm->forward_comm(this);
 
-  // iterate
+  // iterate until convergence
 
   for (itr = 0; itr < itr_max; itr++) {
 
@@ -304,8 +311,8 @@ void FixPolarizeBEMICC::compute_induced_charges()
     for (int i = 0; i < nlocal; i++) {
       if (!(mask[i] & groupbit)) continue;
 
-      double q_free = q_real[i];
-      double qtmp = q[i] - q_free;
+      double q_free = q[i];
+      double qtmp = q_scaled[i] - q_free;
       double Ex = efield_pair[i][0];
       double Ey = efield_pair[i][1];
       double Ez = efield_pair[i][2];
@@ -321,10 +328,10 @@ void FixPolarizeBEMICC::compute_induced_charges()
 
       double ndotE = epsilon0e2q * (Ex * norm[i][0] + Ey * norm[i][1] + Ez * norm[i][2]) /
           (4 * MY_PI) / epsilon[i];
-      double q_bound = q[i] - q_free;
+      double q_bound = q_scaled[i] - q_free;
       q_bound = (1 - omega) * q_bound +
           omega * ((1.0 / em[i] - 1) * q_free - (ed[i] / em[i]) * ndotE * area[i]);
-      q[i] = q_free + q_bound;
+      q_scaled[i] = q_free + q_bound;
 
       // Eq. (11) in Tyagi et al., with f from Eq. (6)
       // NOTE: Tyagi et al. defined the normal vector n_i pointing
@@ -344,6 +351,8 @@ void FixPolarizeBEMICC::compute_induced_charges()
       if (tol < r) tol = r;
     }
 
+    // communicate q_scaled for efield compute in the next iteration
+
     comm->forward_comm(this);
 
     MPI_Allreduce(&tol, &rho, 1, MPI_DOUBLE, MPI_MAX, world);
@@ -351,6 +360,30 @@ void FixPolarizeBEMICC::compute_induced_charges()
   }
 
   iterations = itr;
+
+  // compute the total induced charges of the interface particles
+  // for interface particles: set the charge to be the sum of unscaled (free) charges and induced charges
+
+  double tmp = 0;
+  for (int i = 0; i < nlocal; i++) {
+    if (!(mask[i] & groupbit)) continue;
+
+    double q_bound = q_scaled[i] - q[i];
+    tmp += q_bound;
+    q[i] = q_scaled[i];
+  }
+
+  // ensure sum of all induced charges being zero
+
+  int ncount = group->count(igroup);
+  double sum = 0;
+  MPI_Allreduce(&tmp, &sum, 1, MPI_DOUBLE, MPI_SUM, world);
+  double qboundave = sum/(double)ncount;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (!(mask[i] & groupbit)) continue;
+    q[i] -=  qboundave;
+  }
 }
 
 /* ---------------------------------------------------------------------- */
@@ -422,7 +455,7 @@ int FixPolarizeBEMICC::pack_forward_comm(int n, int *list, double *buf, int /*pb
                                          int * /*pbc*/)
 {
   int m;
-  for (m = 0; m < n; m++) buf[m] = atom->q[list[m]];
+  for (m = 0; m < n; m++) buf[m] = atom->q_scaled[list[m]];
   return n;
 }
 
@@ -431,9 +464,10 @@ int FixPolarizeBEMICC::pack_forward_comm(int n, int *list, double *buf, int /*pb
 void FixPolarizeBEMICC::unpack_forward_comm(int n, int first, double *buf)
 {
   int i, m;
-  for (m = 0, i = first; m < n; m++, i++) atom->q[i] = buf[m];
+  for (m = 0, i = first; m < n; m++, i++) atom->q_scaled[i] = buf[m];
 }
 
+
 /* ----------------------------------------------------------------------
    set dielectric params for the atoms in the group
 ------------------------------------------------------------------------- */
@@ -444,7 +478,7 @@ void FixPolarizeBEMICC::set_dielectric_params(double ediff, double emean, double
   double *area = atom->area;
   double *ed = atom->ed;
   double *em = atom->em;
-  double *q_unscaled = atom->q_unscaled;
+  double *q = atom->q;
   double *epsilon = atom->epsilon;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
@@ -455,7 +489,7 @@ void FixPolarizeBEMICC::set_dielectric_params(double ediff, double emean, double
       em[i] = emean;
       if (areai > 0) area[i] = areai;
       if (epsiloni > 0) epsilon[i] = epsiloni;
-      if (set_charge) q_unscaled[i] = qvalue;
+      if (set_charge) q[i] = qvalue;
     }
   }
 }
diff --git a/src/DIELECTRIC/fix_polarize_functional.cpp b/src/DIELECTRIC/fix_polarize_functional.cpp
index fa9629e54a..d751bcd6ce 100644
--- a/src/DIELECTRIC/fix_polarize_functional.cpp
+++ b/src/DIELECTRIC/fix_polarize_functional.cpp
@@ -59,6 +59,8 @@ using namespace MathExtra;
 using namespace MathConst;
 using namespace MathSpecial;
 
+//#define _POLARIZE_DEBUG
+
 enum { REAL2SCALED = 0, SCALED2REAL = 1 };
 
 static constexpr double EPSILON = 1.0e-6;
@@ -263,7 +265,7 @@ void FixPolarizeFunctional::init()
   // need a full neighbor list w/ Newton off and ghost neighbors
   // built whenever re-neighboring occurs
 
-  neighbor->add_request(this, NeighConst::REQ_FULL | NeighConst::REQ_OCCASIONAL);
+  neighbor->add_request(this, NeighConst::REQ_FULL);
 
   if (force->kspace)
     g_ewald = force->kspace->g_ewald;
@@ -356,7 +358,7 @@ void FixPolarizeFunctional::pre_force(int)
 
 void FixPolarizeFunctional::update_induced_charges()
 {
-  // convert all ions from scaled charges (q) to real q by multiplying with epsilon
+  // convert all ion charges from scaled (q_scaled) to real q by multiplying with epsilon
 
   charge_rescaled(SCALED2REAL);
 
@@ -375,13 +377,27 @@ void FixPolarizeFunctional::update_induced_charges()
 
   // assign charges to the particles in the group
 
+  double *q_scaled = atom->q_scaled;
   double *q = atom->q;
+  double *epsilon = atom->epsilon;
   int nlocal = atom->nlocal;
+  double tmp = 0;
 
   for (int i = 0; i < nlocal; i++) {
     if (induced_charge_idx[i] < 0) continue;
     int idx = induced_charge_idx[i];
     q[i] = -induced_charges[idx] / (4 * MY_PI);
+    q_scaled[i] = q[i] / epsilon[i];
+    tmp += q_scaled[i];
+  }
+
+  double sum = 0;
+  MPI_Allreduce(&tmp, &sum, 1, MPI_DOUBLE, MPI_SUM, world);
+  double qboundave = sum/(double)num_induced_charges;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (induced_charge_idx[i] < 0) continue;
+    q_scaled[i] -=  qboundave;
   }
 
   // revert to scaled charges to calculate forces
@@ -396,19 +412,21 @@ void FixPolarizeFunctional::update_induced_charges()
 
 void FixPolarizeFunctional::charge_rescaled(int scaled2real)
 {
+  double *q_scaled = atom->q_scaled;
   double *q = atom->q;
-  double *q_real = atom->q_unscaled;
   double *epsilon = atom->epsilon;
   int nlocal = atom->nlocal;
 
-  if (scaled2real) {
+  if (scaled2real == SCALED2REAL) {
     for (int i = 0; i < nlocal; i++)
-      if (induced_charge_idx[i] < 0) q[i] = q_real[i];
+      if (induced_charge_idx[i] < 0) q_scaled[i] = q[i];
   } else {
     for (int i = 0; i < nlocal; i++)
-      if (induced_charge_idx[i] < 0) q[i] = q_real[i] / epsilon[i];
+      if (induced_charge_idx[i] < 0) q_scaled[i] = q[i] / epsilon[i];
   }
 
+  // communicate q_scaled to neighboring procs
+
   comm->forward_comm(this);
 }
 
@@ -530,7 +548,7 @@ int FixPolarizeFunctional::pack_forward_comm(int n, int *list, double *buf, int
                                              int * /*pbc*/)
 {
   int m;
-  for (m = 0; m < n; m++) buf[m] = atom->q[list[m]];
+  for (m = 0; m < n; m++) buf[m] = atom->q_scaled[list[m]];
   return n;
 }
 
@@ -539,7 +557,7 @@ int FixPolarizeFunctional::pack_forward_comm(int n, int *list, double *buf, int
 void FixPolarizeFunctional::unpack_forward_comm(int n, int first, double *buf)
 {
   int i, m;
-  for (m = 0, i = first; m < n; m++, i++) atom->q[i] = buf[m];
+  for (m = 0, i = first; m < n; m++, i++) atom->q_scaled[i] = buf[m];
 }
 
 /* ----------------------------------------------------------------------
@@ -818,7 +836,7 @@ void FixPolarizeFunctional::calculate_qiRqw_cutoff()
   int *mask = atom->mask;
   tagint *tag = atom->tag;
   double **x = atom->x;
-  double *q = atom->q_unscaled;
+  double *q = atom->q_scaled;
   double *epsilon = atom->epsilon;
   double *area = atom->area;
   double **norm = atom->mu;
@@ -1027,7 +1045,7 @@ void FixPolarizeFunctional::set_dielectric_params(double ediff, double emean, do
   double *area = atom->area;
   double *ed = atom->ed;
   double *em = atom->em;
-  double *q_unscaled = atom->q_unscaled;
+  double *q = atom->q;
   double *epsilon = atom->epsilon;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
@@ -1038,7 +1056,7 @@ void FixPolarizeFunctional::set_dielectric_params(double ediff, double emean, do
       em[i] = emean;
       if (areai > 0) area[i] = areai;
       if (epsiloni > 0) epsilon[i] = epsiloni;
-      if (set_charge) q_unscaled[i] = qvalue;
+      if (set_charge) q[i] = qvalue;
     }
   }
 }
diff --git a/src/DIELECTRIC/msm_dielectric.cpp b/src/DIELECTRIC/msm_dielectric.cpp
index 5e74b36937..fd09e40acc 100644
--- a/src/DIELECTRIC/msm_dielectric.cpp
+++ b/src/DIELECTRIC/msm_dielectric.cpp
@@ -245,7 +245,7 @@ void MSMDielectric::compute(int eflag, int vflag)
   // energy includes self-energy correction
 
   if (evflag_atom) {
-    double *q = atom->q;
+    double *q = atom->q_scaled;
     int nlocal = atom->nlocal;
 
     if (eflag_atom) {
@@ -287,7 +287,7 @@ void MSMDielectric::fieldforce()
   // (mx,my,mz) = global coords of moving stencil pt
   // ek = 3 components of E-field on particle
 
-  double *q = atom->q;
+  double *q = atom->q_scaled;
   double **x = atom->x;
   double **f = atom->f;
   double *eps = atom->epsilon;
diff --git a/src/DIELECTRIC/pair_coul_cut_dielectric.cpp b/src/DIELECTRIC/pair_coul_cut_dielectric.cpp
index cc9ac0e986..c44552176d 100644
--- a/src/DIELECTRIC/pair_coul_cut_dielectric.cpp
+++ b/src/DIELECTRIC/pair_coul_cut_dielectric.cpp
@@ -69,7 +69,7 @@ void PairCoulCutDielectric::compute(int eflag, int vflag)
 
   double **x = atom->x;
   double **f = atom->f;
-  double *q = atom->q;
+  double *q = atom->q_scaled;
   double *eps = atom->epsilon;
   double **norm = atom->mu;
   double *curvature = atom->curvature;
@@ -165,10 +165,11 @@ double PairCoulCutDielectric::single(int i, int j, int /*itype*/, int /*jtype*/,
                                      double factor_coul, double /*factor_lj*/, double &fforce)
 {
   double r2inv, phicoul, ei, ej;
+  double *q = atom->q_scaled;
   double *eps = atom->epsilon;
 
   r2inv = 1.0 / rsq;
-  fforce = force->qqrd2e * atom->q[i] * atom->q[j] * sqrt(r2inv) * eps[i];
+  fforce = force->qqrd2e * q[i] * q[j] * sqrt(r2inv) * eps[i];
 
   double eng = 0.0;
   if (eps[i] == 1)
@@ -179,7 +180,7 @@ double PairCoulCutDielectric::single(int i, int j, int /*itype*/, int /*jtype*/,
     ej = 0;
   else
     ej = eps[j];
-  phicoul = force->qqrd2e * atom->q[i] * atom->q[j] * sqrt(r2inv);
+  phicoul = force->qqrd2e * q[i] * q[j] * sqrt(r2inv);
   phicoul *= 0.5 * (ei + ej);
   eng += factor_coul * phicoul;
 
diff --git a/src/DIELECTRIC/pair_coul_long_dielectric.cpp b/src/DIELECTRIC/pair_coul_long_dielectric.cpp
index f84adac488..a30bd6e45e 100644
--- a/src/DIELECTRIC/pair_coul_long_dielectric.cpp
+++ b/src/DIELECTRIC/pair_coul_long_dielectric.cpp
@@ -73,7 +73,7 @@ void PairCoulLongDielectric::compute(int eflag, int vflag)
 
   double **x = atom->x;
   double **f = atom->f;
-  double *q = atom->q;
+  double *q = atom->q_scaled;
   double *eps = atom->epsilon;
   double **norm = atom->mu;
   double *curvature = atom->curvature;
diff --git a/src/DIELECTRIC/pair_lj_cut_coul_cut_dielectric.cpp b/src/DIELECTRIC/pair_lj_cut_coul_cut_dielectric.cpp
index 3bcd1e39ff..3d8793457b 100644
--- a/src/DIELECTRIC/pair_lj_cut_coul_cut_dielectric.cpp
+++ b/src/DIELECTRIC/pair_lj_cut_coul_cut_dielectric.cpp
@@ -74,7 +74,7 @@ void PairLJCutCoulCutDielectric::compute(int eflag, int vflag)
 
   double **x = atom->x;
   double **f = atom->f;
-  double *q = atom->q;
+  double *q = atom->q_scaled;
   double *eps = atom->epsilon;
   double **norm = atom->mu;
   double *curvature = atom->curvature;
@@ -195,11 +195,12 @@ double PairLJCutCoulCutDielectric::single(int i, int j, int itype, int jtype, do
                                           double factor_coul, double factor_lj, double &fforce)
 {
   double r2inv, r6inv, forcecoul, forcelj, phicoul, ei, ej, philj;
+  double *q = atom->q_scaled;
   double *eps = atom->epsilon;
 
   r2inv = 1.0 / rsq;
   if (rsq < cut_coulsq[itype][jtype])
-    forcecoul = force->qqrd2e * atom->q[i] * atom->q[j] * sqrt(r2inv) * eps[i];
+    forcecoul = force->qqrd2e * q[i] * q[j] * sqrt(r2inv) * eps[i];
   else
     forcecoul = 0.0;
   if (rsq < cut_ljsq[itype][jtype]) {
@@ -219,7 +220,7 @@ double PairLJCutCoulCutDielectric::single(int i, int j, int itype, int jtype, do
   else
     ej = eps[j];
   if (rsq < cut_coulsq[itype][jtype]) {
-    phicoul = force->qqrd2e * atom->q[i] * atom->q[j] * sqrt(r2inv);
+    phicoul = force->qqrd2e * q[i] * q[j] * sqrt(r2inv);
     phicoul *= 0.5 * (ei + ej);
     eng += factor_coul * phicoul;
   }
diff --git a/src/DIELECTRIC/pair_lj_cut_coul_debye_dielectric.cpp b/src/DIELECTRIC/pair_lj_cut_coul_debye_dielectric.cpp
index 17a534f814..edc5ede368 100644
--- a/src/DIELECTRIC/pair_lj_cut_coul_debye_dielectric.cpp
+++ b/src/DIELECTRIC/pair_lj_cut_coul_debye_dielectric.cpp
@@ -75,7 +75,7 @@ void PairLJCutCoulDebyeDielectric::compute(int eflag, int vflag)
 
   double **x = atom->x;
   double **f = atom->f;
-  double *q = atom->q;
+  double *q = atom->q_scaled;
   double *eps = atom->epsilon;
   double **norm = atom->mu;
   double *curvature = atom->curvature;
@@ -197,6 +197,7 @@ double PairLJCutCoulDebyeDielectric::single(int i, int j, int itype, int jtype,
 {
   double r2inv, r6inv, forcecoul, forcelj, phicoul, ei, ej, philj;
   double r, rinv, screening;
+  double *q = atom->q_scaled;
   double *eps = atom->epsilon;
 
   r2inv = 1.0 / rsq;
@@ -204,7 +205,7 @@ double PairLJCutCoulDebyeDielectric::single(int i, int j, int itype, int jtype,
     r = sqrt(rsq);
     rinv = 1.0 / r;
     screening = exp(-kappa * r);
-    forcecoul = force->qqrd2e * atom->q[i] * atom->q[j] * screening * (kappa + rinv) * eps[i];
+    forcecoul = force->qqrd2e * q[i] * q[j] * screening * (kappa + rinv) * eps[i];
   } else
     forcecoul = 0.0;
   if (rsq < cut_ljsq[itype][jtype]) {
@@ -224,7 +225,7 @@ double PairLJCutCoulDebyeDielectric::single(int i, int j, int itype, int jtype,
   else
     ej = eps[j];
   if (rsq < cut_coulsq[itype][jtype]) {
-    phicoul = force->qqrd2e * atom->q[i] * atom->q[j] * rinv * screening;
+    phicoul = force->qqrd2e * q[i] * q[j] * rinv * screening;
     phicoul *= 0.5 * (ei + ej);
     eng += factor_coul * phicoul;
   }
diff --git a/src/DIELECTRIC/pair_lj_cut_coul_long_dielectric.cpp b/src/DIELECTRIC/pair_lj_cut_coul_long_dielectric.cpp
index 7a4d071d89..ab64d2bd10 100644
--- a/src/DIELECTRIC/pair_lj_cut_coul_long_dielectric.cpp
+++ b/src/DIELECTRIC/pair_lj_cut_coul_long_dielectric.cpp
@@ -81,7 +81,7 @@ void PairLJCutCoulLongDielectric::compute(int eflag, int vflag)
 
   double **x = atom->x;
   double **f = atom->f;
-  double *q = atom->q;
+  double *q = atom->q_scaled;
   double *eps = atom->epsilon;
   double **norm = atom->mu;
   double *curvature = atom->curvature;
@@ -252,6 +252,7 @@ double PairLJCutCoulLongDielectric::single(int i, int j, int itype, int jtype, d
   double r2inv, r6inv, r, grij, expm2, t, erfc, ei, ej, prefactor;
   double fraction, table, forcecoul, forcelj, phicoul, philj;
   int itable;
+  double *q = atom->q_scaled;
   double *eps = atom->epsilon;
 
   r2inv = 1.0 / rsq;
@@ -262,7 +263,7 @@ double PairLJCutCoulLongDielectric::single(int i, int j, int itype, int jtype, d
       expm2 = exp(-grij * grij);
       t = 1.0 / (1.0 + EWALD_P * grij);
       erfc = t * (A1 + t * (A2 + t * (A3 + t * (A4 + t * A5)))) * expm2;
-      prefactor = force->qqrd2e * atom->q[i] * atom->q[j] / r;
+      prefactor = force->qqrd2e * q[i] * q[j] / r;
       forcecoul = prefactor * (erfc + EWALD_F * grij * expm2);
       if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * prefactor;
     } else {
@@ -272,10 +273,10 @@ double PairLJCutCoulLongDielectric::single(int i, int j, int itype, int jtype, d
       itable >>= ncoulshiftbits;
       fraction = (rsq_lookup_single.f - rtable[itable]) * drtable[itable];
       table = ftable[itable] + fraction * dftable[itable];
-      forcecoul = atom->q[i] * atom->q[j] * table;
+      forcecoul = q[i] * q[j] * table;
       if (factor_coul < 1.0) {
         table = ctable[itable] + fraction * dctable[itable];
-        prefactor = atom->q[i] * atom->q[j] * table;
+        prefactor = q[i] * q[j] * table;
         forcecoul -= (1.0 - factor_coul) * prefactor;
       }
     }
@@ -301,12 +302,11 @@ double PairLJCutCoulLongDielectric::single(int i, int j, int itype, int jtype, d
     ej = eps[j];
   if (rsq < cut_coulsq) {
     if (!ncoultablebits || rsq <= tabinnersq)
-      phicoul = prefactor * (ei + ej) * erfc;
+      phicoul = prefactor * 0.5 * (ei + ej) * erfc;
     else {
       table = etable[itable] + fraction * detable[itable];
-      phicoul = atom->q[i] * atom->q[j] * (ei + ej) * table;
+      phicoul = q[i] * q[j] * 0.5 * (ei + ej) * table;
     }
-    phicoul *= 0.5;
     if (factor_coul < 1.0) phicoul -= (1.0 - factor_coul) * prefactor;
     eng += phicoul;
   }
diff --git a/src/DIELECTRIC/pair_lj_cut_coul_msm_dielectric.cpp b/src/DIELECTRIC/pair_lj_cut_coul_msm_dielectric.cpp
index ffcfc8c76b..4bec063e78 100644
--- a/src/DIELECTRIC/pair_lj_cut_coul_msm_dielectric.cpp
+++ b/src/DIELECTRIC/pair_lj_cut_coul_msm_dielectric.cpp
@@ -101,7 +101,7 @@ void PairLJCutCoulMSMDielectric::compute(int eflag, int vflag)
 
   double **x = atom->x;
   double **f = atom->f;
-  double *q = atom->q;
+  double *q = atom->q_scaled;
   double *eps = atom->epsilon;
   double **norm = atom->mu;
   double *curvature = atom->curvature;
@@ -269,15 +269,17 @@ void PairLJCutCoulMSMDielectric::compute(int eflag, int vflag)
 double PairLJCutCoulMSMDielectric::single(int i, int j, int itype, int jtype, double rsq,
                                           double factor_coul, double factor_lj, double &fforce)
 {
-  double r2inv, r6inv, r, egamma, fgamma, prefactor;
+  double r2inv, r6inv, r, egamma, fgamma, ei, ej, prefactor;
   double fraction, table, forcecoul, forcelj, phicoul, philj;
   int itable;
+  double *q = atom->q_scaled;
+  double *eps = atom->epsilon;
 
   r2inv = 1.0 / rsq;
   if (rsq < cut_coulsq) {
     if (!ncoultablebits || rsq <= tabinnersq) {
       r = sqrt(rsq);
-      prefactor = force->qqrd2e * atom->q[i] * atom->q[j] / r;
+      prefactor = force->qqrd2e * q[i] * q[j] / r;
       egamma = 1.0 - (r / cut_coul) * force->kspace->gamma(r / cut_coul);
       fgamma = 1.0 + (rsq / cut_coulsq) * force->kspace->dgamma(r / cut_coul);
       forcecoul = prefactor * fgamma;
@@ -289,10 +291,10 @@ double PairLJCutCoulMSMDielectric::single(int i, int j, int itype, int jtype, do
       itable >>= ncoulshiftbits;
       fraction = (rsq_lookup_single.f - rtable[itable]) * drtable[itable];
       table = ftable[itable] + fraction * dftable[itable];
-      forcecoul = atom->q[i] * atom->q[j] * table;
+      forcecoul = q[i] * q[j] * table;
       if (factor_coul < 1.0) {
         table = ctable[itable] + fraction * dctable[itable];
-        prefactor = atom->q[i] * atom->q[j] * table;
+        prefactor = q[i] * q[j] * table;
         forcecoul -= (1.0 - factor_coul) * prefactor;
       }
     }
@@ -308,12 +310,20 @@ double PairLJCutCoulMSMDielectric::single(int i, int j, int itype, int jtype, do
   fforce = (forcecoul + factor_lj * forcelj) * r2inv;
 
   double eng = 0.0;
+  if (eps[i] == 1)
+    ei = 0;
+  else
+    ei = eps[i];
+  if (eps[j] == 1)
+    ej = 0;
+  else
+    ej = eps[j];
   if (rsq < cut_coulsq) {
     if (!ncoultablebits || rsq <= tabinnersq)
-      phicoul = prefactor * egamma;
+      phicoul = prefactor * 0.5 * (ei + ej) * egamma;
     else {
       table = etable[itable] + fraction * detable[itable];
-      phicoul = atom->q[i] * atom->q[j] * table;
+      phicoul = q[i] * q[j] * 0.5 * (ei + ej) * table;
     }
     if (factor_coul < 1.0) phicoul -= (1.0 - factor_coul) * prefactor;
     eng += phicoul;
diff --git a/src/DIELECTRIC/pair_lj_long_coul_long_dielectric.cpp b/src/DIELECTRIC/pair_lj_long_coul_long_dielectric.cpp
index 069567822a..c1696170cc 100644
--- a/src/DIELECTRIC/pair_lj_long_coul_long_dielectric.cpp
+++ b/src/DIELECTRIC/pair_lj_long_coul_long_dielectric.cpp
@@ -92,7 +92,7 @@ void PairLJLongCoulLongDielectric::compute(int eflag, int vflag)
 
   double **x = atom->x, *x0 = x[0];
   double **f = atom->f;
-  double *q = atom->q;
+  double *q = atom->q_scaled;
   double *eps = avec->epsilon;
   double **norm = avec->mu;
   double *curvature = avec->curvature;
@@ -292,7 +292,8 @@ double PairLJLongCoulLongDielectric::single(int i, int j, int itype, int jtype,
                                             double factor_coul, double factor_lj, double &fforce)
 {
   double r2inv, r6inv, force_coul, force_lj, ei, ej;
-  double g2 = g_ewald_6 * g_ewald_6, g6 = g2 * g2 * g2, g8 = g6 * g2, *q = atom->q;
+  double g2 = g_ewald_6 * g_ewald_6, g6 = g2 * g2 * g2, g8 = g6 * g2;
+  double *q = atom->q_scaled;
   double *eps = avec->epsilon;
 
   double eng = 0.0;
diff --git a/src/DIELECTRIC/pppm_dielectric.cpp b/src/DIELECTRIC/pppm_dielectric.cpp
index e11ed107ad..ca87964e76 100644
--- a/src/DIELECTRIC/pppm_dielectric.cpp
+++ b/src/DIELECTRIC/pppm_dielectric.cpp
@@ -20,8 +20,11 @@
 
 #include "atom.h"
 #include "atom_vec_dielectric.h"
+#include "comm.h"
 #include "domain.h"
 #include "error.h"
+#include "fft3d_wrap.h"
+#include "force.h"
 #include "grid3d.h"
 #include "math_const.h"
 #include "math_special.h"
@@ -35,8 +38,8 @@ using namespace MathSpecial;
 
 #define SMALL 0.00001
 
-enum{REVERSE_RHO};
-enum{FORWARD_IK,FORWARD_AD,FORWARD_IK_PERATOM,FORWARD_AD_PERATOM};
+enum {REVERSE_RHO};
+enum {FORWARD_IK,FORWARD_AD,FORWARD_IK_PERATOM,FORWARD_AD_PERATOM};
 
 #ifdef FFT_SINGLE
 #define ZEROF 0.0f
@@ -55,6 +58,7 @@ PPPMDielectric::PPPMDielectric(LAMMPS *_lmp) : PPPM(_lmp)
   efield = nullptr;
   phi = nullptr;
   potflag = 0;
+  use_qscaled = true;
 
   // no warnings about non-neutral systems from qsum_qsq()
   warn_nonneutral = 2;
@@ -93,6 +97,10 @@ void PPPMDielectric::compute(int eflag, int vflag)
     natoms_original = atom->natoms;
   }
 
+  // recompute the average epsilon of all the atoms
+
+  compute_ave_epsilon();
+
   // return if there are no charges or dipoles
 
   if (qsqsum == 0.0) return;
@@ -136,6 +144,7 @@ void PPPMDielectric::compute(int eflag, int vflag)
   // return gradients (electric fields) in 3d brick decomposition
   // also performs per-atom calculations via poisson_peratom()
 
+  double energy_before_poisson = energy;
   poisson();
 
   // all procs communicate E-field values
@@ -168,18 +177,47 @@ void PPPMDielectric::compute(int eflag, int vflag)
   if (evflag_atom) fieldforce_peratom();
 
   // sum global energy across procs and add in volume-dependent term
+  // NOTE: electrostatic energy is not linearly dependent on charge density (unlike forces)
+  //       recall that we are using atom->q_scaled for make_rho()
+  //       need to switch to atom->q to compute energy (elong)
+  //       also need to use average epsilon (assuming that global dielectric = 1, not set in the input script)
 
   const double qscale = qqrd2e * scale;
 
   if (eflag_global) {
+
+    energy = energy_before_poisson;
+
+    // switch to unscaled charges to find charge density
+
+    use_qscaled = false;
+
+    // redo the charge density
+
+    make_rho();
+
+    // communicate for charge density
+
+    gc->reverse_comm(Grid3d::KSPACE,this,REVERSE_RHO,1,sizeof(FFT_SCALAR),
+                     gc_buf1,gc_buf2,MPI_FFT_SCALAR);
+
+    brick2fft();
+
+    // compute electrostatic energy with the unscaled charges and average epsilon
+
+    poisson();
+
     double energy_all;
     MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
     energy = energy_all;
-
     energy *= 0.5*volume;
-    energy -= g_ewald*qsqsum/MY_PIS +
+      energy -= g_ewald*qsqsum/MY_PIS +
       MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);
-    energy *= qscale;
+    energy *= (qscale/epsilon_ave);
+
+    // revert to qscaled charges (for force in the next time step)
+
+    use_qscaled = true;
   }
 
   // sum global virial across procs
@@ -225,6 +263,130 @@ void PPPMDielectric::compute(int eflag, int vflag)
   if (triclinic) domain->lamda2x(atom->nlocal);
 }
 
+/* ----------------------------------------------------------------------
+   compute the average dielectric constant of all the atoms
+   NOTE: for dielectric use cases
+------------------------------------------------------------------------- */
+
+void PPPMDielectric::compute_ave_epsilon()
+{
+  const double * const epsilon = atom->epsilon;
+  const int nlocal = atom->nlocal;
+  double epsilon_local(0.0);
+
+#if defined(_OPENMP)
+#pragma omp parallel for default(shared) reduction(+:epsilon_local)
+#endif
+  for (int i = 0; i < nlocal; i++) {
+    epsilon_local += epsilon[i];
+  }
+
+  MPI_Allreduce(&epsilon_local,&epsilon_ave,1,MPI_DOUBLE,MPI_SUM,world);
+  epsilon_ave /= (double)atom->natoms;
+}
+
+/* ----------------------------------------------------------------------
+   compute qsum,qsqsum,q2 and give error/warning if not charge neutral
+   called initially, when particle count changes, when charges are changed
+------------------------------------------------------------------------- */
+
+void PPPMDielectric::qsum_qsq(int warning_flag)
+{
+  const double * const q = atom->q;
+  const double * const epsilon = atom->epsilon;
+  const int nlocal = atom->nlocal;
+  double qsum_local(0.0), qsqsum_local(0.0), qsqsume_local(0.0);
+  double qsqsume;
+
+#if defined(_OPENMP)
+#pragma omp parallel for default(shared) reduction(+:qsum_local,qsqsum_local)
+#endif
+  for (int i = 0; i < nlocal; i++) {
+    qsum_local += q[i];
+    qsqsum_local += q[i]*q[i];
+    qsqsume_local += q[i]*q[i]/epsilon[i];
+  }
+
+  MPI_Allreduce(&qsum_local,&qsum,1,MPI_DOUBLE,MPI_SUM,world);
+  MPI_Allreduce(&qsqsum_local,&qsqsum,1,MPI_DOUBLE,MPI_SUM,world);
+  MPI_Allreduce(&qsqsume_local,&qsqsume,1,MPI_DOUBLE,MPI_SUM,world);
+
+  if ((qsqsum == 0.0) && (comm->me == 0) && warn_nocharge && warning_flag) {
+    error->warning(FLERR,"Using kspace solver on system with no charge");
+    warn_nocharge = 0;
+  }
+
+  // q2 is used to compute the mesh spacing, here using qsqsume to match with regular pppm
+  q2 = qsqsume * force->qqrd2e; //q2 = qsqsum * force->qqrd2e;
+
+  // not yet sure of the correction needed for non-neutral systems
+  // so issue warning or error
+
+  if (fabs(qsum) > SMALL) {
+    std::string message = fmt::format("System is not charge neutral, net "
+                                      "charge = {:.8}",qsum);
+    if (!warn_nonneutral) error->all(FLERR,message);
+    if (warn_nonneutral == 1 && comm->me == 0) error->warning(FLERR,message);
+    warn_nonneutral = 2;
+  }
+}
+
+/* ----------------------------------------------------------------------
+   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
+   NOTE: compute charge density using q_scaled if use_qscaled==true
+                        else using unscaled charge values
+------------------------------------------------------------------------- */
+
+void PPPMDielectric::make_rho()
+{
+  int l,m,n,nx,ny,nz,mx,my,mz;
+  FFT_SCALAR dx,dy,dz,x0,y0,z0;
+
+  // clear 3d density array
+
+  memset(&(density_brick[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 *q = atom->q_scaled;
+  if (use_qscaled == false) q = atom->q;
+  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 * q[i];
+    for (n = nlower; n <= nupper; n++) {
+      mz = n+nz;
+      y0 = z0*rho1d[2][n];
+      for (m = nlower; m <= nupper; m++) {
+        my = m+ny;
+        x0 = y0*rho1d[1][m];
+        for (l = nlower; l <= nupper; l++) {
+          mx = l+nx;
+          density_brick[mz][my][mx] += x0*rho1d[0][l];
+        }
+      }
+    }
+  }
+}
+
 /* ----------------------------------------------------------------------
    interpolate from grid to get electric field & force on my particles for ik
 ------------------------------------------------------------------------- */
@@ -241,7 +403,7 @@ void PPPMDielectric::fieldforce_ik()
   // (mx,my,mz) = global coords of moving stencil pt
   // ek = 3 components of E-field on particle
 
-  double *q = atom->q;
+  double *q = atom->q_scaled;
   double **x = atom->x;
   double **f = atom->f;
   double *eps = atom->epsilon;
@@ -321,7 +483,7 @@ void PPPMDielectric::fieldforce_ad()
   // (mx,my,mz) = global coords of moving stencil pt
   // ek = 3 components of E-field on particle
 
-  double *q = atom->q;
+  double *q = atom->q_scaled;
   double **x = atom->x;
   double **f = atom->f;
   double *eps = atom->epsilon;
@@ -362,8 +524,8 @@ void PPPMDielectric::fieldforce_ad()
 
     // convert E-field to force and substract self forces
 
-    const double qfactor = qqrd2e * scale;
-    double qtmp = eps[i]*q[i];
+    const double qfactor = qqrd2e * scale * eps[i];
+    double qtmp = q[i];
 
     s1 = x[i][0]*hx_inv;
     s2 = x[i][1]*hy_inv;
@@ -406,10 +568,9 @@ void PPPMDielectric::slabcorr()
 {
   // compute local contribution to global dipole moment
 
-  double *q = atom->q;
+  double *q = atom->q_scaled;
   double **x = atom->x;
   double *eps = atom->epsilon;
-
   double zprd_slab = domain->zprd*slab_volfactor;
   int nlocal = atom->nlocal;
 
@@ -426,7 +587,6 @@ void PPPMDielectric::slabcorr()
 
   double dipole_r2 = 0.0;
   if (eflag_atom || fabs(qsum) > SMALL) {
-
     for (int i = 0; i < nlocal; i++)
       dipole_r2 += q[i]*x[i][2]*x[i][2];
 
diff --git a/src/DIELECTRIC/pppm_dielectric.h b/src/DIELECTRIC/pppm_dielectric.h
index 69d7fd66a0..c266488561 100644
--- a/src/DIELECTRIC/pppm_dielectric.h
+++ b/src/DIELECTRIC/pppm_dielectric.h
@@ -36,11 +36,16 @@ class PPPMDielectric : public PPPM {
 
  protected:
   void slabcorr() override;
-
+  void make_rho() override;
   void fieldforce_ik() override;
   void fieldforce_ad() override;
+  void qsum_qsq(int warning_flag = 1) override;
 
   class AtomVecDielectric *avec;
+  bool use_qscaled;
+
+  void compute_ave_epsilon();
+  double epsilon_ave;
 };
 
 }    // namespace LAMMPS_NS
diff --git a/src/DIELECTRIC/pppm_disp_dielectric.cpp b/src/DIELECTRIC/pppm_disp_dielectric.cpp
index e8d4d621df..e8a954616b 100644
--- a/src/DIELECTRIC/pppm_disp_dielectric.cpp
+++ b/src/DIELECTRIC/pppm_disp_dielectric.cpp
@@ -20,6 +20,7 @@
 
 #include "atom.h"
 #include "atom_vec_dielectric.h"
+#include "comm.h"
 #include "domain.h"
 #include "error.h"
 #include "force.h"
@@ -63,6 +64,7 @@ PPPMDispDielectric::PPPMDispDielectric(LAMMPS *_lmp) : PPPMDisp(_lmp)
   group_group_enable = 0;
 
   mu_flag = 0;
+  use_qscaled = true;
 
   // no warnings about non-neutral systems from qsum_qsq()
   warn_nonneutral = 2;
@@ -415,11 +417,56 @@ void PPPMDispDielectric::compute(int eflag, int vflag)
     natoms_original = atom->natoms;
   }
 
+  // recompute the average epsilon of all the atoms
+
+  compute_ave_epsilon();
+
   // sum energy across procs and add in volume-dependent term
 
   const double qscale = force->qqrd2e * scale;
 
   if (eflag_global) {
+
+    if (function[0]) {
+
+      // switch to unscaled charges to find charge density
+
+      use_qscaled = false;
+
+      make_rho_c();
+
+      gc->reverse_comm(Grid3d::KSPACE,this,REVERSE_RHO,1,sizeof(FFT_SCALAR),
+                     gc_buf1,gc_buf2,MPI_FFT_SCALAR);
+
+      brick2fft(nxlo_in,nylo_in,nzlo_in,nxhi_in,nyhi_in,nzhi_in,
+                density_brick,density_fft,work1,remap);
+
+      // compute electrostatic energy with the unscaled charges and average epsilon
+
+      if (differentiation_flag == 1) {
+        poisson_ad(work1,work2,density_fft,fft1,fft2,
+                  nx_pppm,ny_pppm,nz_pppm,nfft,
+                  nxlo_fft,nylo_fft,nzlo_fft,nxhi_fft,nyhi_fft,nzhi_fft,
+                  nxlo_in,nylo_in,nzlo_in,nxhi_in,nyhi_in,nzhi_in,
+                  energy_1,greensfn,
+                  virial_1,vg,vg2,
+                  u_brick,v0_brick,v1_brick,v2_brick,v3_brick,v4_brick,v5_brick);
+
+      } else {
+        poisson_ik(work1,work2,density_fft,fft1,fft2,
+                  nx_pppm,ny_pppm,nz_pppm,nfft,
+                  nxlo_fft,nylo_fft,nzlo_fft,nxhi_fft,nyhi_fft,nzhi_fft,
+                  nxlo_in,nylo_in,nzlo_in,nxhi_in,nyhi_in,nzhi_in,
+                  energy_1,greensfn,
+                  fkx,fky,fkz,fkx2,fky2,fkz2,
+                  vdx_brick,vdy_brick,vdz_brick,virial_1,vg,vg2,
+                  u_brick,v0_brick,v1_brick,v2_brick,v3_brick,v4_brick,v5_brick);
+
+        gc->forward_comm(Grid3d::KSPACE,this,FORWARD_IK,3,sizeof(FFT_SCALAR),
+                        gc_buf1,gc_buf2,MPI_FFT_SCALAR);
+      }
+    }
+
     double energy_all;
     MPI_Allreduce(&energy_1,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
     energy_1 = energy_all;
@@ -434,6 +481,10 @@ void PPPMDispDielectric::compute(int eflag, int vflag)
     energy_6 += - MY_PI*MY_PIS/(6*volume)*pow(g_ewald_6,3)*csumij +
       1.0/12.0*pow(g_ewald_6,6)*csum;
     energy_1 *= qscale;
+
+    // revert to qscaled charges (for force in the next time step)
+
+    use_qscaled = true;
   }
 
   // sum virial across procs
@@ -494,6 +545,127 @@ void PPPMDispDielectric::compute(int eflag, int vflag)
   if (triclinic) domain->lamda2x(atom->nlocal);
 }
 
+/* ----------------------------------------------------------------------
+   compute the average dielectric constant of all the atoms
+   NOTE: for dielectric use cases
+------------------------------------------------------------------------- */
+
+void PPPMDispDielectric::compute_ave_epsilon()
+{
+  const double * const epsilon = atom->epsilon;
+  const int nlocal = atom->nlocal;
+  double epsilon_local(0.0);
+
+#if defined(_OPENMP)
+#pragma omp parallel for default(shared) reduction(+:epsilon_local)
+#endif
+  for (int i = 0; i < nlocal; i++) {
+    epsilon_local += epsilon[i];
+  }
+
+  MPI_Allreduce(&epsilon_local,&epsilon_ave,1,MPI_DOUBLE,MPI_SUM,world);
+  epsilon_ave /= (double)atom->natoms;
+}
+
+/* ----------------------------------------------------------------------
+   compute qsum,qsqsum,q2 and give error/warning if not charge neutral
+   called initially, when particle count changes, when charges are changed
+------------------------------------------------------------------------- */
+
+void PPPMDispDielectric::qsum_qsq(int warning_flag)
+{
+  const double * const q = atom->q;
+  const double * const epsilon = atom->epsilon;
+  const int nlocal = atom->nlocal;
+  double qsum_local(0.0), qsqsum_local(0.0), qsqsume_local(0.0);
+  double qsqsume;
+
+#if defined(_OPENMP)
+#pragma omp parallel for default(shared) reduction(+:qsum_local,qsqsum_local)
+#endif
+  for (int i = 0; i < nlocal; i++) {
+    qsum_local += q[i];
+    qsqsum_local += q[i]*q[i];
+    qsqsume_local += q[i]*q[i]/epsilon[i];
+  }
+
+  MPI_Allreduce(&qsum_local,&qsum,1,MPI_DOUBLE,MPI_SUM,world);
+  MPI_Allreduce(&qsqsum_local,&qsqsum,1,MPI_DOUBLE,MPI_SUM,world);
+  MPI_Allreduce(&qsqsume_local,&qsqsume,1,MPI_DOUBLE,MPI_SUM,world);
+
+  if ((qsqsum == 0.0) && (comm->me == 0) && warn_nocharge && warning_flag) {
+    error->warning(FLERR,"Using kspace solver on system with no charge");
+    warn_nocharge = 0;
+  }
+
+  // q2 is used to compute the mesh spacing, here using qsqsume to match with regular pppm
+  q2 = qsqsume * force->qqrd2e; //q2 = qsqsum * force->qqrd2e;
+
+  // not yet sure of the correction needed for non-neutral systems
+  // so issue warning or error
+
+  if (fabs(qsum) > SMALL) {
+    std::string message = fmt::format("System is not charge neutral, net "
+                                      "charge = {:.8}",qsum);
+    if (!warn_nonneutral) error->all(FLERR,message);
+    if (warn_nonneutral == 1 && comm->me == 0) error->warning(FLERR,message);
+    warn_nonneutral = 2;
+  }
+}
+
+/* ----------------------------------------------------------------------
+   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 PPPMDispDielectric::make_rho_c()
+{
+  int l,m,n,nx,ny,nz,mx,my,mz;
+  FFT_SCALAR dx,dy,dz,x0,y0,z0;
+
+  // clear 3d density array
+
+  memset(&(density_brick[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 *q = atom->q_scaled;
+  if (use_qscaled == false) q = atom->q;
+  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, order, rho_coeff, rho1d);
+
+    z0 = delvolinv * q[i];
+    for (n = nlower; n <= nupper; n++) {
+      mz = n+nz;
+      y0 = z0*rho1d[2][n];
+      for (m = nlower; m <= nupper; m++) {
+        my = m+ny;
+        x0 = y0*rho1d[1][m];
+        for (l = nlower; l <= nupper; l++) {
+          mx = l+nx;
+          density_brick[mz][my][mx] += x0*rho1d[0][l];
+        }
+      }
+    }
+  }
+}
+
 /* ----------------------------------------------------------------------
    interpolate from grid to get electric field & force on my particles
    for ik scheme
diff --git a/src/DIELECTRIC/pppm_disp_dielectric.h b/src/DIELECTRIC/pppm_disp_dielectric.h
index 043834ca7b..7412762818 100644
--- a/src/DIELECTRIC/pppm_disp_dielectric.h
+++ b/src/DIELECTRIC/pppm_disp_dielectric.h
@@ -37,11 +37,17 @@ class PPPMDispDielectric : public PPPMDisp {
   int potflag;    // 1/0 if per-atom electrostatic potential phi is needed
 
  protected:
+  void make_rho_c() override;
   void fieldforce_c_ik() override;
   void fieldforce_c_ad() override;
   void fieldforce_c_peratom() override;
+  void qsum_qsq(int warning_flag = 1) override;
 
   class AtomVecDielectric *avec;
+  bool use_qscaled;
+
+  void compute_ave_epsilon();
+  double epsilon_ave;
   int mu_flag;
 };
 
diff --git a/src/OPENMP/pair_lj_cut_coul_cut_dielectric_omp.cpp b/src/OPENMP/pair_lj_cut_coul_cut_dielectric_omp.cpp
index d312914ef1..6546b60830 100644
--- a/src/OPENMP/pair_lj_cut_coul_cut_dielectric_omp.cpp
+++ b/src/OPENMP/pair_lj_cut_coul_cut_dielectric_omp.cpp
@@ -101,7 +101,7 @@ void PairLJCutCoulCutDielectricOMP::eval(int iifrom, int iito, ThrData *const th
 
   const auto *_noalias const x = (dbl3_t *) atom->x[0];
   auto *_noalias const f = (dbl3_t *) thr->get_f()[0];
-  const double *_noalias const q = atom->q;
+  const double *_noalias const q = atom->q_scaled;
   const double *_noalias const eps = atom->epsilon;
   const auto *_noalias const norm = (dbl3_t *) atom->mu[0];
   const double *_noalias const curvature = atom->curvature;
diff --git a/src/OPENMP/pair_lj_cut_coul_debye_dielectric_omp.cpp b/src/OPENMP/pair_lj_cut_coul_debye_dielectric_omp.cpp
index ef754decfc..79f76c8343 100644
--- a/src/OPENMP/pair_lj_cut_coul_debye_dielectric_omp.cpp
+++ b/src/OPENMP/pair_lj_cut_coul_debye_dielectric_omp.cpp
@@ -101,7 +101,7 @@ void PairLJCutCoulDebyeDielectricOMP::eval(int iifrom, int iito, ThrData *const
 
   const dbl3_t *_noalias const x = (dbl3_t *) atom->x[0];
   dbl3_t *_noalias const f = (dbl3_t *) thr->get_f()[0];
-  const double *_noalias const q = atom->q;
+  const double *_noalias const q = atom->q_scaled;
   const double *_noalias const eps = atom->epsilon;
   const dbl3_t *_noalias const norm = (dbl3_t *) atom->mu[0];
   const double *_noalias const curvature = atom->curvature;
diff --git a/src/OPENMP/pair_lj_cut_coul_long_dielectric_omp.cpp b/src/OPENMP/pair_lj_cut_coul_long_dielectric_omp.cpp
index c311a76ef7..c0f4693c58 100644
--- a/src/OPENMP/pair_lj_cut_coul_long_dielectric_omp.cpp
+++ b/src/OPENMP/pair_lj_cut_coul_long_dielectric_omp.cpp
@@ -103,7 +103,7 @@ void PairLJCutCoulLongDielectricOMP::eval(int iifrom, int iito, ThrData *const t
 
   const auto *_noalias const x = (dbl3_t *) atom->x[0];
   auto *_noalias const f = (dbl3_t *) thr->get_f()[0];
-  const double *_noalias const q = atom->q;
+  const double *_noalias const q = atom->q_scaled;
   const double *_noalias const eps = atom->epsilon;
   const auto *_noalias const norm = (dbl3_t *) atom->mu[0];
   const double *_noalias const curvature = atom->curvature;
diff --git a/src/atom.cpp b/src/atom.cpp
index 3828df6be9..947023b496 100644
--- a/src/atom.cpp
+++ b/src/atom.cpp
@@ -214,7 +214,7 @@ Atom::Atom(LAMMPS *_lmp) : Pointers(_lmp)
 
   // DIELECTRIC package
 
-  area = ed = em = epsilon = curvature = q_unscaled = nullptr;
+  area = ed = em = epsilon = curvature = q_scaled = nullptr;
 
   // end of customization section
   // --------------------------------------------------------------------
@@ -563,7 +563,7 @@ void Atom::peratom_create()
   add_peratom("em",&em,DOUBLE,0);
   add_peratom("epsilon",&epsilon,DOUBLE,0);
   add_peratom("curvature",&curvature,DOUBLE,0);
-  add_peratom("q_unscaled",&q_unscaled,DOUBLE,0);
+  add_peratom("q_scaled",&q_scaled,DOUBLE,0);
 
   // end of customization section
   // --------------------------------------------------------------------
@@ -2951,7 +2951,7 @@ void *Atom::extract(const char *name)
   if (strcmp(name,"em") == 0) return (void *) em;
   if (strcmp(name,"epsilon") == 0) return (void *) epsilon;
   if (strcmp(name,"curvature") == 0) return (void *) curvature;
-  if (strcmp(name,"q_unscaled") == 0) return (void *) q_unscaled;
+  if (strcmp(name,"q_scaled") == 0) return (void *) q_scaled;
 
   // end of customization section
   // --------------------------------------------------------------------
diff --git a/src/atom.h b/src/atom.h
index e8b1330651..2e663c628a 100644
--- a/src/atom.h
+++ b/src/atom.h
@@ -173,7 +173,7 @@ class Atom : protected Pointers {
 
   // DIELECTRIC package
 
-  double *area, *ed, *em, *epsilon, *curvature, *q_unscaled;
+  double *area, *ed, *em, *epsilon, *curvature, *q_scaled;
 
   // end of customization section
   // --------------------------------------------------------------------
diff --git a/src/kspace.h b/src/kspace.h
index 6cc8f89138..61ab15c1d9 100644
--- a/src/kspace.h
+++ b/src/kspace.h
@@ -116,7 +116,7 @@ class KSpace : protected Pointers {
 
   // public so can be called by commands that change charge
 
-  void qsum_qsq(int warning_flag = 1);
+  virtual void qsum_qsq(int warning_flag = 1);
 
   // general child-class methods
 
diff --git a/src/set.cpp b/src/set.cpp
index 67ae59cd30..3d2bf0ed1a 100644
--- a/src/set.cpp
+++ b/src/set.cpp
@@ -855,8 +855,11 @@ void Set::set(int keyword)
     else if (keyword == VX) atom->v[i][0] = dvalue;
     else if (keyword == VY) atom->v[i][1] = dvalue;
     else if (keyword == VZ) atom->v[i][2] = dvalue;
-    else if (keyword == CHARGE) atom->q[i] = dvalue;
-    else if (keyword == MASS) {
+    else if (keyword == CHARGE) {
+      atom->q[i] = dvalue;
+      // ensure that scaled charges are consistent the new charge value
+      if (atom->epsilon) atom->q_scaled[i] = dvalue / atom->epsilon[i];
+    } else if (keyword == MASS) {
       if (dvalue <= 0.0) error->one(FLERR,"Invalid mass in set command");
       atom->rmass[i] = dvalue;
     }
@@ -1087,14 +1090,11 @@ void Set::set(int keyword)
     else if (keyword == EPSILON) {
       if (dvalue >= 0.0) {
 
-        // compute the unscaled charge value (i.e. atom->q_unscaled)
         // assign the new local dielectric constant
-        // update both the scaled and unscaled charge values
+        // update both the scaled charge value
 
-        double q_unscaled = atom->q[i] * atom->epsilon[i];
         atom->epsilon[i] = dvalue;
-        atom->q[i] = q_unscaled / dvalue;
-        atom->q_unscaled[i] = q_unscaled;
+        atom->q_scaled[i] = atom->q[i] / dvalue;
       }
     }
 
diff --git a/unittest/force-styles/tests/kspace-pppm_dielectric.yaml b/unittest/force-styles/tests/kspace-pppm_dielectric.yaml
index 2ad84b3478..a14356301c 100644
--- a/unittest/force-styles/tests/kspace-pppm_dielectric.yaml
+++ b/unittest/force-styles/tests/kspace-pppm_dielectric.yaml
@@ -1,6 +1,6 @@
 ---
-lammps_version: 2 Jun 2022
-date_generated: Sun Jun  5 11:47:58 2022
+lammps_version: 22 Dec 2022
+date_generated: Sun Jan  8 14:06:48 2023
 epsilon: 1.5e-13
 skip_tests:
 prerequisites: ! |
@@ -23,67 +23,67 @@ init_coul: 0
 init_stress: ! |2-
    0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
 init_forces: ! |2
-    1 -9.5704063239170824e-03  1.3440401839656547e-03  3.9938351166047105e-03
-    2  4.2546747552535428e-03 -4.6394704924901569e-03 -1.8639383905026156e-03
-    3 -6.2689436339003596e-04 -1.1923832294754119e-04  3.1024334072044251e-04
-    4  2.8904653542133046e-03  4.0048286752677591e-04 -1.1881079811383200e-03
-    5  2.8067481142527392e-03  1.1249538111229643e-03 -7.8992757739952202e-04
-    6  1.0523465137115526e-02  5.7444835008325897e-03 -9.4171156316220847e-03
-    7 -6.0442808532858950e-03 -5.8528007629445105e-03  5.2900162100588104e-03
-    8 -2.9095138717290920e-03 -9.0050388447901180e-03  4.4014558257614057e-03
-    9  3.2886707142652893e-03  4.4648703522243427e-03  3.8291771586486379e-04
-   10 -9.0909991198697099e-04  1.7647045143455965e-03 -2.4537591861546526e-06
-   11 -1.5002813696751375e-03  2.4776755287400877e-03 -2.6041426800205674e-04
-   12  8.2255495247906246e-03 -6.8126223463819621e-03 -8.3464377417941820e-04
-   13 -2.8154328633420100e-03  1.8626508137035600e-03  9.2570078297929057e-04
-   14 -2.9991011482475407e-03  2.2193552769033294e-03  4.5675765704394564e-04
-   15 -2.4422674405748931e-03  1.3779216258678446e-03  1.5087422119394080e-04
-   16 -6.3065226817821248e-03  8.0994785421873144e-03  7.7948169511488970e-03
-   17  2.1712932472423164e-03 -7.7092752926742779e-03 -9.9843183915653676e-03
-   18  1.2250633518865794e-02  2.3931580719148072e-02 -2.3032783639349450e-02
-   19 -4.5739325868170509e-03 -1.2083425237677518e-02  1.3127920656307673e-02
-   20 -6.3509675186534543e-03 -1.1554317452880160e-02  1.2513366661911501e-02
-   21  7.2370604334770303e-03  9.4214770437382137e-03 -2.0582480710702272e-02
-   22 -5.0302962829462793e-03 -2.3954421805390562e-03  9.6853611331667918e-03
-   23 -3.1397831881129165e-03 -4.7636081827106625e-03  9.7394906276265063e-03
-   24  9.1592606135343238e-04  2.7705153422388124e-02 -5.8648725733196610e-03
-   25  1.9389266275156363e-03 -1.0450682262511200e-02  2.6060479368781131e-03
-   26 -3.5468346395038284e-03 -1.5670009179897540e-02  1.9825521897942600e-03
-   27 -1.3499137089964794e-02  2.6091966390167483e-02 -1.5098953238706988e-02
-   28  9.4647573055174836e-03 -1.4728281444592208e-02  8.8519568594155693e-03
-   29  6.2965813400663821e-03 -1.2246582589825026e-02  6.7066960491971583e-03
+    1 -9.0240864973911948e-03  1.4038916148678216e-03  3.7037121250649771e-03
+    2  3.7877030534016216e-03 -4.8550977429375076e-03 -2.3689918060965398e-03
+    3 -5.9655389539660028e-04 -1.6265556490413675e-04  3.4628862920435083e-04
+    4  2.8246800647794197e-03  5.0706139861645446e-04 -1.3613405543800795e-03
+    5  2.7706763415996510e-03  1.3192664085953656e-03 -6.4693146379363187e-04
+    6  9.8521150188877558e-03  7.2316371685612465e-03 -1.1869654289057200e-02
+    7 -5.9825556046998350e-03 -6.9265836429670815e-03  6.8668494607014635e-03
+    8 -2.5144593251176608e-03 -1.0702207637180206e-02  5.9414687019552171e-03
+    9  3.1921755971977732e-03  5.5790682925965689e-03  9.1596521903762108e-04
+   10 -8.9583980735863833e-04  1.9183664785442267e-03 -2.6841558552602918e-04
+   11 -1.4675402166109490e-03  2.6118239195873335e-03 -6.9553971036785136e-04
+   12  7.9459309545569581e-03 -7.4250715157140254e-03  6.5918066145004848e-04
+   13 -2.7128154347354000e-03  2.0238602255809855e-03  4.5300484074527176e-04
+   14 -2.9991301386514883e-03  2.3792297052533519e-03  1.8849132343493573e-04
+   15 -2.3885112958716326e-03  1.4922929171770360e-03 -2.7211922381807559e-04
+   16 -5.9894972992190739e-03  7.5170483687749361e-03  9.2783833083170258e-03
+   17  2.3467271943640022e-03 -7.1634029709802790e-03 -1.3724376775234111e-02
+   18  1.2628594191704279e-02  2.6944878208688711e-02 -2.4146105092988153e-02
+   19 -4.4694266391438467e-03 -1.3503908711633225e-02  1.3447394876877700e-02
+   20 -6.8292702443397109e-03 -1.2242896845309152e-02  1.2749828615434848e-02
+   21  8.9532168947175951e-03  9.4647709143936692e-03 -2.0318893339010807e-02
+   22 -5.1021986421169729e-03 -2.1596132999930984e-03  1.0204115119223008e-02
+   23 -4.9434674418902015e-03 -5.0964957389006815e-03  9.7116877087425283e-03
+   24  9.8240113311758438e-04  3.0357874946112112e-02 -4.8489323929702798e-03
+   25  2.3014421568304448e-03 -1.2251539432392653e-02  3.9375681971606583e-03
+   26 -3.8219153726559581e-03 -1.6978221833764549e-02  1.3052944135472525e-03
+   27 -1.4886455790134088e-02  2.8641111159034106e-02 -1.6294597887912415e-02
+   28  9.8847687498804567e-03 -1.5935389770801115e-02  9.4004886798317450e-03
+   29  7.1532922942957209e-03 -1.3989097018906178e-02  7.7061762404265401e-03
 run_vdwl: 0
 run_coul: 0
 run_stress: ! |2-
    0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
 run_forces: ! |2
-    1 -9.5512367296693348e-03  1.3481986228832394e-03  4.0360496920819845e-03
-    2  4.2310095077217249e-03 -4.6520267411891266e-03 -1.8873891848711444e-03
-    3 -6.2657635097199582e-04 -1.1896339414013587e-04  3.1228053952145037e-04
-    4  2.8929119769660746e-03  3.9833099588961228e-04 -1.1948966254162630e-03
-    5  2.8036607264009267e-03  1.1250824895483655e-03 -7.9977001904481899e-04
-    6  1.0517775415299555e-02  5.7363868487737019e-03 -9.4753365739491131e-03
-    7 -6.0484479858843492e-03 -5.8590676715711935e-03  5.3263656833568585e-03
-    8 -2.8902731245966352e-03 -9.0037174743433958e-03  4.4633941645762475e-03
-    9  3.2719440551636165e-03  4.4604792993012321e-03  3.4264400361289614e-04
-   10 -9.1214291344286910e-04  1.7671332579526759e-03 -1.0718060181645637e-05
-   11 -1.5036983613157409e-03  2.4849579312452809e-03 -2.6767344531034297e-04
-   12  8.2370000539625232e-03 -6.8149518655650299e-03 -7.9828139230475636e-04
-   13 -2.8188931760921637e-03  1.8650744557518759e-03  9.1426138691727598e-04
-   14 -3.0020345361480650e-03  2.2216374831535396e-03  4.4819872834903844e-04
-   15 -2.4434489580208332e-03  1.3746131141371646e-03  1.3692983090348620e-04
-   16 -6.3290560550783867e-03  8.1123385965123744e-03  7.7558668570284657e-03
-   17  2.1893044938489075e-03 -7.6981064388344927e-03 -9.9464802743007962e-03
-   18  1.2326575373962534e-02  2.4035949736671783e-02 -2.2966596407517550e-02
-   19 -4.5945299269544680e-03 -1.2120131566157018e-02  1.3105619820168259e-02
-   20 -6.3956048954198911e-03 -1.1610478112584980e-02  1.2479235193738019e-02
-   21  7.2449337028217256e-03  9.2775525645935504e-03 -2.0499726162547745e-02
-   22 -5.0276728424607717e-03 -2.3257295140813555e-03  9.6505434927086856e-03
-   23 -3.1414718347009523e-03 -4.7050318098388480e-03  9.7005400386281204e-03
-   24  9.4171169493705606e-04  2.7631000472134227e-02 -5.8324765686082400e-03
-   25  1.9197188048042241e-03 -1.0419799501124184e-02  2.5868533348110123e-03
-   26 -3.5519466566892540e-03 -1.5629585791677170e-02  1.9678442785116389e-03
-   27 -1.3518109991418836e-02  2.6086722910789121e-02 -1.5022292876012720e-02
-   28  9.4744709425719827e-03 -1.4720145605462569e-02  8.8095347267501904e-03
-   29  6.3041275904037030e-03 -1.2247723292768253e-02  6.6654758184014956e-03
+    1 -9.0041469902975083e-03  1.4079128645891239e-03  3.7454874558270081e-03
+    2  3.7645787277636967e-03 -4.8683493328976549e-03 -2.3920956093199269e-03
+    3 -5.9624682956190922e-04 -1.6235758995387017e-04  3.4824512095250714e-04
+    4  2.8273064132848390e-03  5.0498009019496828e-04 -1.3680472151764869e-03
+    5  2.7676323391099414e-03  1.3190769273285616e-03 -6.5648622300019488e-04
+    6  9.8455723855407316e-03  7.2241664557800610e-03 -1.1926898901957301e-02
+    7 -5.9860293302087903e-03 -6.9349881691303499e-03  6.9034801543381387e-03
+    8 -2.4951707858607963e-03 -1.0699604634547560e-02  6.0017607963033636e-03
+    9  3.1758296804222314e-03  5.5730049066296042e-03  8.7764129419339075e-04
+   10 -8.9871735397559530e-04  1.9204608360439955e-03 -2.7650828365437897e-04
+   11 -1.4708156364411449e-03  2.6184932735207446e-03 -7.0255510826413572e-04
+   12  7.9563087247412623e-03 -7.4261878100576606e-03  6.9343641340987520e-04
+   13 -2.7163219002503949e-03  2.0255295949472103e-03  4.4263040419780290e-04
+   14 -3.0018627327104740e-03  2.3815439416672098e-03  1.8058576207462599e-04
+   15 -2.3893711712619039e-03  1.4900511475556947e-03 -2.8561791170206334e-04
+   16 -6.0101351141619999e-03  7.5306307779994938e-03  9.2409981605329223e-03
+   17  2.3632302529971224e-03 -7.1550898973440896e-03 -1.3685304978041121e-02
+   18  1.2708227256223214e-02  2.7049952948538600e-02 -2.4071103911405968e-02
+   19 -4.4911945713642055e-03 -1.3539135003090926e-02  1.3425202892529920e-02
+   20 -6.8770712710612223e-03 -1.2299906456646187e-02  1.2713597968374448e-02
+   21  8.9604122906157130e-03  9.3126036471350218e-03 -2.0232925593705698e-02
+   22 -5.0980838764751795e-03 -2.0851853158015952e-03  1.0166764658842283e-02
+   23 -4.9465919390184829e-03 -5.0342268804687137e-03  9.6703603172035180e-03
+   24  1.0076873262936329e-03  3.0277862968197698e-02 -4.8175616854022434e-03
+   25  2.2805335344714185e-03 -1.2216062086673396e-02  3.9138629522124067e-03
+   26 -3.8256342727503762e-03 -1.6935631184734987e-02  1.2877457401954645e-03
+   27 -1.4905842380482990e-02  2.8637174134366720e-02 -1.6214467912845618e-02
+   28  9.8946081070436093e-03 -1.5927143980741659e-02  9.3566889655639814e-03
+   29  7.1613091173755462e-03 -1.3989576172406041e-02  7.6610842777234383e-03
 ...
diff --git a/unittest/force-styles/tests/kspace-pppm_disp_dielectric.yaml b/unittest/force-styles/tests/kspace-pppm_disp_dielectric.yaml
index 351583833e..f9c27ae6b1 100644
--- a/unittest/force-styles/tests/kspace-pppm_disp_dielectric.yaml
+++ b/unittest/force-styles/tests/kspace-pppm_disp_dielectric.yaml
@@ -1,7 +1,7 @@
 ---
-lammps_version: 2 Jun 2022
+lammps_version: 22 Dec 2022
 tags: slow
-date_generated: Sun Jun  5 11:47:58 2022
+date_generated: Sun Jan  8 14:06:49 2023
 epsilon: 2.5e-13
 skip_tests: intel
 prerequisites: ! |
@@ -34,67 +34,67 @@ init_coul: 0
 init_stress: ! |2-
    0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
 init_forces: ! |2
-    1  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-    2  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-    3  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-    4  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-    5  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-    6  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-    7  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-    8  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-    9  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   10  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   11  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   12  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   13  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   14  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   15  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   16  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   17  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   18  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   19  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   20  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   21  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   22  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   23  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   24  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   25  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   26  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   27  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   28  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   29  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+    1 -9.1192078690927932e-03  1.4495617840522479e-03  3.7798782771053465e-03
+    2  3.7720899254057526e-03 -4.8637874581151279e-03 -2.3371541480691365e-03
+    3 -6.0032537658623244e-04 -1.6174100878128198e-04  3.4786454607935922e-04
+    4  2.8398487096903918e-03  5.0004390045606500e-04 -1.3558332192793528e-03
+    5  2.7943938157179654e-03  1.3113410511708477e-03 -6.6215289856192978e-04
+    6  9.8362075577179461e-03  7.2739615272235623e-03 -1.1714152557391479e-02
+    7 -5.9559839224051260e-03 -6.9918320119199745e-03  6.8227727728859596e-03
+    8 -2.4399896047000913e-03 -1.0754611680064429e-02  5.8948703359515253e-03
+    9  3.1597958400271190e-03  5.5744919679502365e-03  8.5623544196936118e-04
+   10 -8.9988895267770891e-04  1.9288816927157165e-03 -2.4382039547014042e-04
+   11 -1.4730127425939767e-03  2.6302448858035342e-03 -6.7933682321038571e-04
+   12  7.9689158236380796e-03 -7.4262782929543298e-03  5.8484918329561363e-04
+   13 -2.7165428920156833e-03  2.0206464997513823e-03  4.7064076015990104e-04
+   14 -2.9948781398817007e-03  2.3772997029477038e-03  1.7855804285807838e-04
+   15 -2.4027123761803755e-03  1.4895389891037350e-03 -2.4169692099536508e-04
+   16 -5.9752060401817844e-03  7.5458977707643935e-03  9.2254634359853933e-03
+   17  2.3367828791080640e-03 -7.1804354861216687e-03 -1.3662338955720804e-02
+   18  1.2723013957950361e-02  2.6938230938849794e-02 -2.4150886335690219e-02
+   19 -4.5207624388967646e-03 -1.3481361111124537e-02  1.3395118968446813e-02
+   20 -6.8542055845031563e-03 -1.2254240542987327e-02  1.2708592528169912e-02
+   21  9.0822804375539287e-03  9.4741623354682569e-03 -2.0208748390807822e-02
+   22 -5.1607002273109147e-03 -2.1554386441792364e-03  1.0137593221094414e-02
+   23 -5.0294017225099618e-03 -5.1075253489796983e-03  9.6593012011709153e-03
+   24  1.1608971654866597e-03  3.0346220531046218e-02 -4.8624863038438425e-03
+   25  2.2197610206980654e-03 -1.2254371301472900e-02  3.9404364020363142e-03
+   26 -3.8910412828591224e-03 -1.6968379488344357e-02  1.3010339808255957e-03
+   27 -1.5014639958520915e-02  2.8778190147924231e-02 -1.6304337765055573e-02
+   28  9.9536295161445250e-03 -1.5982701867618014e-02  9.4090685576554995e-03
+   29  7.2008824817774543e-03 -1.4056009482565078e-02  7.7106670584060595e-03
 run_vdwl: 0
 run_coul: 0
 run_stress: ! |2-
    0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
 run_forces: ! |2
-    1  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-    2  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-    3  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-    4  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-    5  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-    6  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-    7  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-    8  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-    9  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   10  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   11  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   12  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   13  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   14  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   15  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   16  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   17  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   18  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   19  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   20  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   21  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   22  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   23  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   24  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   25  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   26  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   27  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   28  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
-   29  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
+    1 -9.0991065994639846e-03  1.4535750492104861e-03  3.8217957245221270e-03
+    2  3.7489325156386984e-03 -4.8769311623033090e-03 -2.3604975864010302e-03
+    3 -6.0001553006534965e-04 -1.6144187873488151e-04  3.4982440017513830e-04
+    4  2.8424883809995676e-03  4.9795877852022675e-04 -1.3625138311279131e-03
+    5  2.7913227296444891e-03  1.3111481937246542e-03 -6.7172814410274863e-04
+    6  9.8296998153309476e-03  7.2663180983279533e-03 -1.1771445504611623e-02
+    7 -5.9594475967707497e-03 -7.0000828437108455e-03  6.8592323866196077e-03
+    8 -2.4207569182026468e-03 -1.0751853186430514e-02  5.9552517741290291e-03
+    9  3.1435054964006007e-03  5.5683293372686578e-03  8.1787001068474329e-04
+   10 -9.0276632005639463e-04  1.9309601881084811e-03 -2.5191931983487042e-04
+   11 -1.4762793586003183e-03  2.6368857130867926e-03 -6.8637341481350104e-04
+   12  7.9792853590517701e-03 -7.4273234202696155e-03  6.1909389673923422e-04
+   13 -2.7200609199609506e-03  2.0222931599477288e-03  4.6025076840837842e-04
+   14 -2.9975714842595757e-03  2.3795677147271693e-03  1.7070935034618187e-04
+   15 -2.4035803890593327e-03  1.4872986431882908e-03 -2.5521007706313665e-04
+   16 -5.9958239437963873e-03  7.5593948903379454e-03  9.1879038237176434e-03
+   17  2.3532867821412715e-03 -7.1720854721499779e-03 -1.3623031991652995e-02
+   18  1.2802728849711620e-02  2.7043356451850818e-02 -2.4076055049592836e-02
+   19 -4.5427582713879946e-03 -1.3516544861853890e-02  1.3373151256499361e-02
+   20 -6.9019406524318157e-03 -1.2311361735980344e-02  1.2672469255604397e-02
+   21  9.0894838284619490e-03  9.3219704340619085e-03 -2.0122780448886900e-02
+   22 -5.1565421781962365e-03 -2.0809952855138998e-03  1.0100250978387133e-02
+   23 -5.0325884857908955e-03 -5.0452106679046201e-03  9.6179958286591281e-03
+   24  1.1857944847112256e-03  3.0266030255059431e-02 -4.8310785851609058e-03
+   25  2.1989851701805347e-03 -1.2218856194106672e-02  3.9166848091745530e-03
+   26 -3.8946116887646418e-03 -1.6925708327337474e-02  1.2834628694644208e-03
+   27 -1.5033949860253316e-02  2.8774217256229383e-02 -1.6224064740787608e-02
+   28  9.9633496453770745e-03 -1.5974432317747287e-02  9.3651772703528761e-03
+   29  7.2089371394108469e-03 -1.4056476809606566e-02  7.6655742905521115e-03
 ...