diff --git a/src/KSPACE/ewald.cpp b/src/KSPACE/ewald.cpp
index a684ce80a5..1acf04e220 100644
--- a/src/KSPACE/ewald.cpp
+++ b/src/KSPACE/ewald.cpp
@@ -106,8 +106,7 @@ void Ewald::init()
 
   // extract short-range Coulombic cutoff from pair style
 
-  scale = 1.0;
-
+  triclinic = domain->triclinic;
   pair_check();
 
   int itmp;
@@ -116,25 +115,12 @@ void Ewald::init()
     error->all(FLERR,"KSpace style is incompatible with Pair style");
   double cutoff = *p_cutoff;
 
-  qsum = qsqsum = 0.0;
-  for (int i = 0; i < atom->nlocal; i++) {
-    qsum += atom->q[i];
-    qsqsum += atom->q[i]*atom->q[i];
-  }
+  // compute qsum & qsqsum and warn if not charge-neutral
 
-  double tmp;
-  MPI_Allreduce(&qsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
-  qsum = tmp;
-  MPI_Allreduce(&qsqsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
-  qsqsum = tmp;
-
-  if (qsqsum == 0.0)
-    error->all(FLERR,"Cannot use kspace solver on system with no charge");
-  if (fabs(qsum) > SMALL && comm->me == 0) {
-    char str[128];
-    sprintf(str,"System is not charge neutral, net charge = %g",qsum);
-    error->warning(FLERR,str);
-  }
+  scale = 1.0;
+  qqrd2e = force->qqrd2e;
+  qsum_qsq(0);
+  natoms_original = atom->natoms;
 
   // set accuracy (force units) from accuracy_relative or accuracy_absolute
 
@@ -143,11 +129,8 @@ void Ewald::init()
 
   // setup K-space resolution
 
-  q2 = qsqsum * force->qqrd2e;
   bigint natoms = atom->natoms;
 
-  triclinic = domain->triclinic;
-
   // use xprd,yprd,zprd even if triclinic so grid size is the same
   // adjust z dimension for 2d slab Ewald
   // 3d Ewald just uses zprd since slab_volfactor = 1.0
@@ -440,7 +423,7 @@ void Ewald::compute(int eflag, int vflag)
 
   // convert E-field to force
 
-  const double qscale = force->qqrd2e * scale;
+  const double qscale = qqrd2e * scale;
 
   for (i = 0; i < nlocal; i++) {
     f[i][0] += qscale * q[i]*ek[i][0];
@@ -448,12 +431,19 @@ void Ewald::compute(int eflag, int vflag)
     if (slabflag != 2) f[i][2] += qscale * q[i]*ek[i][2];
   }
 
-  // global energy
+  // sum global energy across Kspace vevs and add in volume-dependent term
+  // reset qsum and qsqsum if atom count has changed
 
   if (eflag_global) {
     for (k = 0; k < kcount; k++)
       energy += ug[k] * (sfacrl_all[k]*sfacrl_all[k] +
                          sfacim_all[k]*sfacim_all[k]);
+
+    if (atom->natoms != natoms_original) {
+      qsum_qsq(0);
+      natoms_original = atom->natoms;
+    }
+
     energy -= g_ewald*qsqsum/MY_PIS +
       MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);
     energy *= qscale;
@@ -1206,7 +1196,7 @@ void Ewald::slabcorr()
 
   const double e_slabcorr = MY_2PI*(dipole_all*dipole_all -
     qsum*dipole_r2 - qsum*qsum*zprd*zprd/12.0)/volume;
-  const double qscale = force->qqrd2e * scale;
+  const double qscale = qqrd2e * scale;
 
   if (eflag_global) energy += qscale * e_slabcorr;
 
@@ -1337,7 +1327,7 @@ void Ewald::compute_group_group(int groupbit_A, int groupbit_B, int AA_flag)
   MPI_Allreduce(sfacrl_B,sfacrl_B_all,kcount,MPI_DOUBLE,MPI_SUM,world);
   MPI_Allreduce(sfacim_B,sfacim_B_all,kcount,MPI_DOUBLE,MPI_SUM,world);
 
-  const double qscale = force->qqrd2e * scale;
+  const double qscale = qqrd2e * scale;
   double partial_group;
 
   // total group A <--> group B energy
@@ -1437,7 +1427,7 @@ void Ewald::slabcorr_groups(int groupbit_A, int groupbit_B, int AA_flag)
 
   // compute corrections
 
-  const double qscale = force->qqrd2e * scale;
+  const double qscale = qqrd2e * scale;
   const double efact = qscale * MY_2PI/volume;
 
   e2group += efact * (dipole_A*dipole_B - 0.5*(qsum_A*dipole_r2_B +
diff --git a/src/KSPACE/msm.cpp b/src/KSPACE/msm.cpp
index 809de53d7a..523916519c 100644
--- a/src/KSPACE/msm.cpp
+++ b/src/KSPACE/msm.cpp
@@ -166,46 +166,23 @@ void MSM::init()
     error->all(FLERR,"Cannot (yet) use single precision with MSM "
                "(remove -DFFT_SINGLE from Makefile and recompile)");
 
-  pair_check();
-
   // extract short-range Coulombic cutoff from pair style
 
+  triclinic = domain->triclinic;
+  pair_check();
+
   int itmp;
   double *p_cutoff = (double *) force->pair->extract("cut_coul",itmp);
   if (p_cutoff == NULL)
     error->all(FLERR,"KSpace style is incompatible with Pair style");
   cutoff = *p_cutoff;
 
-  // compute qsum & qsqsum and give error if not charge-neutral
+  // compute qsum & qsqsum and error if not charge-neutral
 
-  qsum = qsqsum = 0.0;
-  for (int i = 0; i < atom->nlocal; i++) {
-    qsum += atom->q[i];
-    qsqsum += atom->q[i]*atom->q[i];
-  }
-
-  qqrd2e = force->qqrd2e;
   scale = 1.0;
-
-  double tmp;
-  MPI_Allreduce(&qsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
-  qsum = tmp;
-  MPI_Allreduce(&qsqsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
-  qsqsum = tmp;
-  q2 = qsqsum * force->qqrd2e;
-
-  if (qsqsum == 0.0)
-    error->all(FLERR,"Cannot use kspace solver on system with no charge");
-
-  // not yet sure of the correction needed for non-neutral systems
-
-  if (fabs(qsum) > SMALL) {
-    char str[128];
-    sprintf(str,"System is not charge neutral, net charge = %g",qsum);
-    error->all(FLERR,str);
-  }
-
-  triclinic = domain->triclinic;
+  qqrd2e = force->qqrd2e;
+  qsum_qsq(1);
+  natoms_original = atom->natoms;
 
   // set accuracy (force units) from accuracy_relative or accuracy_absolute
 
@@ -586,16 +563,22 @@ void MSM::compute(int eflag, int vflag)
 
   if (evflag_atom) fieldforce_peratom();
 
-  // total long-range energy
+  // sum global energy across procs and add in self-energy term
+  // reset qsum and qsqsum if atom count has changed
 
-  const double qscale = force->qqrd2e * scale;
+  const double qscale = qqrd2e * scale;
 
   if (eflag_global) {
     double energy_all;
     MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
     energy = energy_all;
 
-    double e_self = qsqsum*gamma(0.0)/cutoff;  // Self-energy term
+    if (atom->natoms != natoms_original) {
+      qsum_qsq(0);
+      natoms_original = atom->natoms;
+    }
+
+    double e_self = qsqsum*gamma(0.0)/cutoff;
     energy -= e_self;
     energy *= 0.5*qscale;
   }
@@ -2796,7 +2779,7 @@ void MSM::fieldforce()
 
     // convert E-field to force
 
-    const double qfactor = force->qqrd2e*scale*q[i];
+    const double qfactor = qqrd2e*scale*q[i];
     f[i][0] += qfactor*ekx;
     f[i][1] += qfactor*eky;
     f[i][2] += qfactor*ekz;
diff --git a/src/KSPACE/msm.h b/src/KSPACE/msm.h
index 3d333e9073..82be217f5b 100644
--- a/src/KSPACE/msm.h
+++ b/src/KSPACE/msm.h
@@ -40,7 +40,6 @@ class MSM : public KSpace {
   double precision;
   int nfactors;
   int *factors;
-  double qsum,qsqsum,q2;
   double qqrd2e;
   double cutoff;
   double volume;
diff --git a/src/KSPACE/pppm.cpp b/src/KSPACE/pppm.cpp
index ac217a6f1d..35a78408ea 100644
--- a/src/KSPACE/pppm.cpp
+++ b/src/KSPACE/pppm.cpp
@@ -48,8 +48,8 @@ using namespace MathSpecial;
 
 #define MAXORDER 7
 #define OFFSET 16384
-#define SMALL 0.00001
 #define LARGE 10000.0
+#define SMALL 0.00001
 #define EPS_HOC 1.0e-7
 
 enum{REVERSE_RHO};
@@ -203,8 +203,6 @@ void PPPM::init()
   // extract short-range Coulombic cutoff from pair style
 
   triclinic = domain->triclinic;
-  scale = 1.0;
-
   pair_check();
 
   int itmp = 0;
@@ -250,26 +248,10 @@ void PPPM::init()
 
   // compute qsum & qsqsum and warn if not charge-neutral
 
-  qsum = qsqsum = 0.0;
-  for (int i = 0; i < atom->nlocal; i++) {
-    qsum += atom->q[i];
-    qsqsum += atom->q[i]*atom->q[i];
-  }
-
-  double tmp;
-  MPI_Allreduce(&qsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
-  qsum = tmp;
-  MPI_Allreduce(&qsqsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
-  qsqsum = tmp;
-  q2 = qsqsum * force->qqrd2e;
-
-  if (qsqsum == 0.0)
-    error->all(FLERR,"Cannot use kspace solver on system with no charge");
-  if (fabs(qsum) > SMALL && me == 0) {
-    char str[128];
-    sprintf(str,"System is not charge neutral, net charge = %g",qsum);
-    error->warning(FLERR,str);
-  }
+  scale = 1.0;
+  qqrd2e = force->qqrd2e;
+  qsum_qsq(0);
+  natoms_original = atom->natoms;
 
   // set accuracy (force units) from accuracy_relative or accuracy_absolute
 
@@ -678,14 +660,20 @@ void PPPM::compute(int eflag, int vflag)
   if (evflag_atom) fieldforce_peratom();
 
   // sum global energy across procs and add in volume-dependent term
+  // reset qsum and qsqsum if atom count has changed
 
-  const double qscale = force->qqrd2e * scale;
+  const double qscale = qqrd2e * scale;
 
   if (eflag_global) {
     double energy_all;
     MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
     energy = energy_all;
 
+    if (atom->natoms != natoms_original) {
+      qsum_qsq(0);
+      natoms_original = atom->natoms;
+    }
+
     energy *= 0.5*volume;
     energy -= g_ewald*qsqsum/MY_PIS +
       MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);
@@ -2460,7 +2448,7 @@ void PPPM::fieldforce_ik()
 
     // convert E-field to force
 
-    const double qfactor = force->qqrd2e * scale * q[i];
+    const double qfactor = qqrd2e * scale * q[i];
     f[i][0] += qfactor*ekx;
     f[i][1] += qfactor*eky;
     if (slabflag != 2) f[i][2] += qfactor*ekz;
@@ -2531,7 +2519,7 @@ void PPPM::fieldforce_ad()
 
     // convert E-field to force and substract self forces
 
-    const double qfactor = force->qqrd2e * scale;
+    const double qfactor = qqrd2e * scale;
 
     s1 = x[i][0]*hx_inv;
     s2 = x[i][1]*hy_inv;
@@ -2954,7 +2942,7 @@ void PPPM::slabcorr()
 
   const double e_slabcorr = MY_2PI*(dipole_all*dipole_all -
     qsum*dipole_r2 - qsum*qsum*zprd*zprd/12.0)/volume;
-  const double qscale = force->qqrd2e * scale;
+  const double qscale = qqrd2e * scale;
 
   if (eflag_global) energy += qscale * e_slabcorr;
 
@@ -3141,7 +3129,7 @@ void PPPM::compute_group_group(int groupbit_A, int groupbit_B, int AA_flag)
 
   poisson_groups(AA_flag);
 
-  const double qscale = force->qqrd2e * scale;
+  const double qscale = qqrd2e * scale;
 
   // total group A <--> group B energy
   // self and boundary correction terms are in compute_group_group.cpp
@@ -3487,7 +3475,7 @@ void PPPM::slabcorr_groups(int groupbit_A, int groupbit_B, int AA_flag)
 
   // compute corrections
 
-  const double qscale = force->qqrd2e * scale;
+  const double qscale = qqrd2e * scale;
   const double efact = qscale * MY_2PI/volume;
 
   e2group += efact * (dipole_A*dipole_B - 0.5*(qsum_A*dipole_r2_B +
diff --git a/src/KSPACE/pppm.h b/src/KSPACE/pppm.h
index bf6cb9bc89..a061c656be 100644
--- a/src/KSPACE/pppm.h
+++ b/src/KSPACE/pppm.h
@@ -53,7 +53,6 @@ class PPPM : public KSpace {
   int me,nprocs;
   int nfactors;
   int *factors;
-  double qsum,qsqsum,q2;
   double cutoff;
   double volume;
   double delxinv,delyinv,delzinv,delvolinv;
diff --git a/src/KSPACE/pppm_cg.cpp b/src/KSPACE/pppm_cg.cpp
index ca8aa2c0c3..53da74f97f 100644
--- a/src/KSPACE/pppm_cg.cpp
+++ b/src/KSPACE/pppm_cg.cpp
@@ -209,14 +209,20 @@ void PPPMCG::compute(int eflag, int vflag)
   if (evflag_atom) fieldforce_peratom();
 
   // sum global energy across procs and add in volume-dependent term
+  // reset qsum and qsqsum if atom count has changed
 
-  const double qscale = force->qqrd2e * scale;
+  const double qscale = qqrd2e * scale;
 
   if (eflag_global) {
     double energy_all;
     MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
     energy = energy_all;
 
+    if (atom->natoms != natoms_original) {
+      qsum_qsq(0);
+      natoms_original = atom->natoms;
+    }
+
     energy *= 0.5*volume;
     energy -= g_ewald*qsqsum/MY_PIS +
       MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);
@@ -406,7 +412,7 @@ void PPPMCG::fieldforce_ik()
 
     // convert E-field to force
 
-    const double qfactor = force->qqrd2e * scale * q[i];
+    const double qfactor = qqrd2e * scale * q[i];
     f[i][0] += qfactor*ekx;
     f[i][1] += qfactor*eky;
     if (slabflag != 2) f[i][2] += qfactor*ekz;
@@ -479,7 +485,7 @@ void PPPMCG::fieldforce_ad()
 
     // convert E-field to force and substract self forces
 
-    const double qfactor = force->qqrd2e * scale;
+    const double qfactor = qqrd2e * scale;
 
     s1 = x[i][0]*hx_inv;
     s2 = x[i][1]*hy_inv;
@@ -618,7 +624,7 @@ void PPPMCG::slabcorr()
 
   const double e_slabcorr = MY_2PI*(dipole_all*dipole_all -
     qsum*dipole_r2 - qsum*qsum*zprd*zprd/12.0)/volume;
-  const double qscale = force->qqrd2e * scale;
+  const double qscale = qqrd2e * scale;
 
   if (eflag_global) energy += qscale * e_slabcorr;
 
diff --git a/src/KSPACE/pppm_disp.cpp b/src/KSPACE/pppm_disp.cpp
index 0750ec02e9..4439502d9f 100755
--- a/src/KSPACE/pppm_disp.cpp
+++ b/src/KSPACE/pppm_disp.cpp
@@ -233,14 +233,9 @@ void PPPMDisp::init()
   deallocate();
   deallocate_peratom(); 
 
-  // set scale
-
-  scale = 1.0;
-
-  triclinic = domain->triclinic;
-
   // check whether cutoff and pair style are set
 
+  triclinic = domain->triclinic;
   pair_check();
 
   int tmp;
@@ -269,7 +264,8 @@ void PPPMDisp::init()
 	case 1:
 	  k = 0; break;
 	case 6:
-	  if ((ewald_mix==GEOMETRIC || ewald_mix==SIXTHPOWER|| mixflag == 1) && mixflag!= 2) { k = 1; break; }
+	  if ((ewald_mix==GEOMETRIC || ewald_mix==SIXTHPOWER || 
+               mixflag == 1) && mixflag!= 2) { k = 1; break; }
 	  else if (ewald_mix==ARITHMETIC && mixflag!=2) { k = 2; break; }
 	  else if (mixflag == 2) { k = 3; break; }
 	default:
@@ -282,41 +278,20 @@ void PPPMDisp::init()
  
 
   // warn, if function[0] is not set but charge attribute is set!
+
   if (!function[0] && atom->q_flag && me == 0) {
     char str[128];
     sprintf(str, "Charges are set, but coulombic solver is not used");
     error->warning(FLERR, str);
   }
 
-  // compute qsum & qsqsum, if function[0] is set, print error if no charges are set or warn if not charge-neutral  
+  // compute qsum & qsqsum, if function[0] is set, warn if not charge-neutral
+
+  scale = 1.0;
+  qqrd2e = force->qqrd2e;
+  natoms_original = atom->natoms;
  
-  if (function[0]) {
-    if (!atom->q_flag) 
-      error->all(FLERR,"Kspace style with selected options "
-                 "requires atom attribute q");
- 
-    qsum = qsqsum = 0.0;
-    for (int i = 0; i < atom->nlocal; i++) {
-      qsum += atom->q[i];
-      qsqsum += atom->q[i]*atom->q[i];
-
-    }
-
-    double tmp;
-    MPI_Allreduce(&qsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
-    qsum = tmp;
-    MPI_Allreduce(&qsqsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
-    qsqsum = tmp;
-
-    if (qsqsum == 0.0)
-      error->all(FLERR,"Cannot use kspace solver with selected options "
-                 "on system with no charge");
-    if (fabs(qsum) > SMALL && me == 0) {
-      char str[128];
-      sprintf(str,"System is not charge neutral, net charge = %g",qsum);
-      error->warning(FLERR,str);
-    }
-  }
+  if (function[0]) qsum_qsq(0);
 
   // if kspace is TIP4P, extract TIP4P params from pair style
   // bond/angle are not yet init(), so insure equilibrium request is valid
@@ -351,8 +326,8 @@ void PPPMDisp::init()
     alpha = qdist / (cos(0.5*theta) * blen);
   }
 
-
   // initialize the pair style to get the coefficients
+
   neighrequest_flag = 0;
   pair->init();
   neighrequest_flag = 1;
@@ -1144,6 +1119,7 @@ void PPPMDisp::compute(int eflag, int vflag)
   }
 
   // sum energy across procs and add in volume-dependent term
+  // reset qsum and qsqsum if atom count has changed
 
   const double qscale = force->qqrd2e * scale;
   if (eflag_global) {
@@ -1156,6 +1132,11 @@ void PPPMDisp::compute(int eflag, int vflag)
     energy_1 *= 0.5*volume;
     energy_6 *= 0.5*volume;
     
+    if (atom->natoms != natoms_original) {
+      qsum_qsq(0);
+      natoms_original = atom->natoms;
+    }
+
     energy_1 -= g_ewald*qsqsum/MY_PIS +
       MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);
     energy_6 += - MY_PI*MY_PIS/(6*volume)*pow(g_ewald_6,3)*csumij +
diff --git a/src/KSPACE/pppm_disp.h b/src/KSPACE/pppm_disp.h
index 26c56d1a0b..e1ca2747a2 100755
--- a/src/KSPACE/pppm_disp.h
+++ b/src/KSPACE/pppm_disp.h
@@ -62,7 +62,6 @@ Variables needed for calculating the 1/r and 1/r^6 potential
   int me,nprocs;
   int nfactors;
   int *factors;
-  double qsum,qsqsum;
   double csumij;
   double csum;
   double *csumi;  //needed as correction term for per atom calculations!
diff --git a/src/KSPACE/pppm_disp_tip4p.cpp b/src/KSPACE/pppm_disp_tip4p.cpp
index d3c2484b00..6967f795ab 100755
--- a/src/KSPACE/pppm_disp_tip4p.cpp
+++ b/src/KSPACE/pppm_disp_tip4p.cpp
@@ -97,7 +97,6 @@ void PPPMDispTIP4P::particle_map_c(double delx, double dely, double delz,
     p2g[i][1] = ny;
     p2g[i][2] = nz;
 
-
     // check that entire stencil around nx,ny,nz will fit in my 3d brick
 
     if (nx+nlow < nxlo || nx+nup > nxhi ||
diff --git a/src/KSPACE/pppm_old.cpp b/src/KSPACE/pppm_old.cpp
index c2959823ff..c2d1cad053 100644
--- a/src/KSPACE/pppm_old.cpp
+++ b/src/KSPACE/pppm_old.cpp
@@ -154,8 +154,6 @@ void PPPMOld::init()
 
   // extract short-range Coulombic cutoff from pair style
 
-  scale = 1.0;
-
   pair_check();
 
   int itmp=0;
@@ -198,25 +196,10 @@ void PPPMOld::init()
 
   // compute qsum & qsqsum and warn if not charge-neutral
 
-  qsum = qsqsum = 0.0;
-  for (int i = 0; i < atom->nlocal; i++) {
-    qsum += atom->q[i];
-    qsqsum += atom->q[i]*atom->q[i];
-  }
-
-  double tmp;
-  MPI_Allreduce(&qsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
-  qsum = tmp;
-  MPI_Allreduce(&qsqsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
-  qsqsum = tmp;
-
-  if (qsqsum == 0.0)
-    error->all(FLERR,"Cannot use kspace solver on system with no charge");
-  if (fabs(qsum) > SMALL && me == 0) {
-    char str[128];
-    sprintf(str,"System is not charge neutral, net charge = %g",qsum);
-    error->warning(FLERR,str);
-  }
+  scale = 1.0;
+  qqrd2e = force->qqrd2e;
+  qsum_qsq(0);
+  natoms_original = atom->natoms;
 
   // set accuracy (force units) from accuracy_relative or accuracy_absolute
 
@@ -742,14 +725,20 @@ void PPPMOld::compute(int eflag, int vflag)
   if (evflag_atom) fieldforce_peratom();
 
   // sum global energy across procs and add in volume-dependent term
+  // reset qsum and qsqsum if atom count has changed
 
-  const double qscale = force->qqrd2e * scale;
+  const double qscale = qqrd2e * scale;
 
   if (eflag_global) {
     double energy_all;
     MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
     energy = energy_all;
 
+    if (atom->natoms != natoms_original) {
+      qsum_qsq(0);
+      natoms_original = atom->natoms;
+    }
+
     energy *= 0.5*volume;
     energy -= g_ewald*qsqsum/MY_PIS +
       MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);
@@ -971,8 +960,6 @@ void PPPMOld::set_grid()
   acons[7][5] = 1755948832039.0 / 36229939200000.0;
   acons[7][6] = 4887769399.0 / 37838389248.0;
 
-  double q2 = qsqsum * force->qqrd2e;
-
   // use xprd,yprd,zprd even if triclinic so grid size is the same
   // adjust z dimension for 2d slab PPPM
   // 3d PPPM just uses zprd since slab_volfactor = 1.0
@@ -2232,7 +2219,7 @@ void PPPMOld::fieldforce()
 
     // convert E-field to force
 
-    const double qfactor = force->qqrd2e * scale * q[i];
+    const double qfactor = qqrd2e * scale * q[i];
     f[i][0] += qfactor*ekx;
     f[i][1] += qfactor*eky;
     if (slabflag != 2) f[i][2] += qfactor*ekz;
@@ -2471,7 +2458,7 @@ void PPPMOld::slabcorr()
 
   const double e_slabcorr = MY_2PI*(dipole_all*dipole_all -
     qsum*dipole_r2 - qsum*qsum*zprd*zprd/12.0)/volume;
-  const double qscale = force->qqrd2e * scale;
+  const double qscale = qqrd2e * scale;
 
   if (eflag_global) energy += qscale * e_slabcorr;
 
@@ -2645,7 +2632,7 @@ void PPPMOld::compute_group_group(int groupbit_A, int groupbit_B, int BA_flag)
 
   poisson_groups(BA_flag);
 
-  const double qscale = force->qqrd2e * scale;
+  const double qscale = qqrd2e * scale;
 
   // total group A <--> group B energy
   // self and boundary correction terms are in compute_group_group.cpp
diff --git a/src/KSPACE/pppm_old.h b/src/KSPACE/pppm_old.h
index 758dd79a49..f8382008f9 100644
--- a/src/KSPACE/pppm_old.h
+++ b/src/KSPACE/pppm_old.h
@@ -52,7 +52,6 @@ class PPPMOld : public KSpace {
   int me,nprocs;
   int nfactors;
   int *factors;
-  double qsum,qsqsum;
   double cutoff;
   double volume;
   double delxinv,delyinv,delzinv,delvolinv;
diff --git a/src/KSPACE/pppm_stagger.cpp b/src/KSPACE/pppm_stagger.cpp
index c8612f090f..23c55b8cee 100755
--- a/src/KSPACE/pppm_stagger.cpp
+++ b/src/KSPACE/pppm_stagger.cpp
@@ -203,14 +203,20 @@ void PPPMStagger::compute(int eflag, int vflag)
   }
 
   // sum global energy across procs and add in volume-dependent term
+  // reset qsum and qsqsum if atom count has changed
 
-  const double qscale = force->qqrd2e * scale;
+  const double qscale = qqrd2e * scale;
 
   if (eflag_global) {
     double energy_all;
     MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
     energy = energy_all;
 
+    if (atom->natoms != natoms_original) {
+      qsum_qsq(0);
+      natoms_original = atom->natoms;
+    }
+
     energy *= 0.5*volume/float(nstagger);
     energy -= g_ewald*qsqsum/MY_PIS +
       MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);
@@ -813,7 +819,7 @@ void PPPMStagger::fieldforce_ik()
 
     // convert E-field to force
 
-    const double qfactor = force->qqrd2e * scale * q[i] / float(nstagger);
+    const double qfactor = qqrd2e * scale * q[i] / float(nstagger);
     f[i][0] += qfactor*ekx;
     f[i][1] += qfactor*eky;
     if (slabflag != 2) f[i][2] += qfactor*ekz;
@@ -884,7 +890,7 @@ void PPPMStagger::fieldforce_ad()
 
     // convert E-field to force and substract self forces
 
-    const double qfactor = force->qqrd2e * scale / float(nstagger);
+    const double qfactor = qqrd2e * scale / float(nstagger);
 
     s1 = x[i][0]*hx_inv + stagger;
     s2 = x[i][1]*hy_inv + stagger;
@@ -899,7 +905,6 @@ void PPPMStagger::fieldforce_ad()
     sf *= 2*q[i]*q[i];
     f[i][1] += qfactor*(eky*q[i] - sf);
 
-
     sf = sf_coeff[4]*sin(2*MY_PI*s3);
     sf += sf_coeff[5]*sin(4*MY_PI*s3);
     sf *= 2*q[i]*q[i];
diff --git a/src/KSPACE/pppm_stagger.h b/src/KSPACE/pppm_stagger.h
index 2b2b5ae1c0..b51f83819b 100755
--- a/src/KSPACE/pppm_stagger.h
+++ b/src/KSPACE/pppm_stagger.h
@@ -52,7 +52,8 @@ class PPPMStagger : public PPPM {
 
 
   inline double gf_denom2(const double &x, const double &y,
-                         const double &z) const {
+                         const double &z) const
+  {
     double sx,sy,sz;
     double x2 = x*x;
     double y2 = y*y;
diff --git a/src/KSPACE/pppm_tip4p.cpp b/src/KSPACE/pppm_tip4p.cpp
index 528a3c1848..decdb1bad4 100644
--- a/src/KSPACE/pppm_tip4p.cpp
+++ b/src/KSPACE/pppm_tip4p.cpp
@@ -223,7 +223,7 @@ void PPPMTIP4P::fieldforce_ik()
 
     // convert E-field to force
 
-    const double qfactor = force->qqrd2e * scale * q[i];
+    const double qfactor = qqrd2e * scale * q[i];
     if (type[i] != typeO) {
       f[i][0] += qfactor*ekx;
       f[i][1] += qfactor*eky;
@@ -326,7 +326,8 @@ void PPPMTIP4P::fieldforce_ad()
     ekz *= hz_inv;
     
     // convert E-field to force and substract self forces
-    const double qfactor = force->qqrd2e * scale;
+
+    const double qfactor = qqrd2e * scale;
 
     s1 = x[i][0]*hx_inv;
     s2 = x[i][1]*hy_inv;