diff --git a/src/GPU/pppm_gpu.cpp b/src/GPU/pppm_gpu.cpp
index ef79325845..75773ba318 100644
--- a/src/GPU/pppm_gpu.cpp
+++ b/src/GPU/pppm_gpu.cpp
@@ -290,13 +290,11 @@ void PPPMGPU::compute(int eflag, int vflag)
 
   if (evflag_atom) fieldforce_peratom();
 
-  // update qsum and qsqsum, if needed
+  // update qsum and qsqsum, if atom count has changed and energy needed
 
-  if (eflag_global || eflag_atom) {
-    if (qsum_update_flag || (atom->natoms != natoms_original)) {
-      qsum_qsq(0);
-      natoms_original = atom->natoms;
-    }
+  if ((eflag_global || eflag_atom) && atom->natoms != natoms_original) {
+    qsum_qsq(0);
+    natoms_original = atom->natoms;
   }
 
   // sum energy across procs and add in volume-dependent term
diff --git a/src/KSPACE/ewald.cpp b/src/KSPACE/ewald.cpp
index 115dc95f94..107dbfbaff 100644
--- a/src/KSPACE/ewald.cpp
+++ b/src/KSPACE/ewald.cpp
@@ -119,7 +119,7 @@ void Ewald::init()
 
   scale = 1.0;
   qqrd2e = force->qqrd2e;
-  qsum_qsq(0);
+  qsum_qsq();
   natoms_original = atom->natoms;
 
   // set accuracy (force units) from accuracy_relative or accuracy_absolute
@@ -431,13 +431,11 @@ void Ewald::compute(int eflag, int vflag)
     if (slabflag != 2) f[i][2] += qscale * q[i]*ek[i][2];
   }
 
-  // update qsum and qsqsum, if needed
+  // update qsum and qsqsum, if atom count has changed and energy needed
 
-  if (eflag_global || eflag_atom) {
-    if (qsum_update_flag || (atom->natoms != natoms_original)) {
-      qsum_qsq(0);
-      natoms_original = atom->natoms;
-    }
+  if ((eflag_global || eflag_atom) && atom->natoms != natoms_original) {
+    qsum_qsq();
+    natoms_original = atom->natoms;
   }
 
   // sum global energy across Kspace vevs and add in volume-dependent term
diff --git a/src/KSPACE/ewald_disp.cpp b/src/KSPACE/ewald_disp.cpp
index f7c2275aa3..c107fc85c7 100644
--- a/src/KSPACE/ewald_disp.cpp
+++ b/src/KSPACE/ewald_disp.cpp
@@ -140,13 +140,7 @@ void EwaldDisp::init()
   pair->init();  // so B is defined
   init_coeffs();
   init_coeff_sums();
-
-  double qsum, qsqsum, bsbsum;
-  qsum = qsqsum = bsbsum = 0.0;
-  if (function[0]) {
-    qsum = sum[0].x;
-    qsqsum = sum[0].x2;
-  }
+  if (function[0]) qsum_qsq();
 
   // turn off coulombic if no charge
 
@@ -156,6 +150,7 @@ void EwaldDisp::init()
     nsums -= 1;
   }
 
+  double bsbsum = 0.0;
   if (function[1]) bsbsum = sum[1].x2;
   if (function[2]) bsbsum = sum[2].x2;
 
@@ -246,9 +241,6 @@ void EwaldDisp::setup()
     error->all(FLERR,"KSpace accuracy too low");
   }
 
-  if (qsum_update_flag)
-    error->all(FLERR,"Kspace style ewald/disp does not support dynamic charges");
-
   bigint natoms = atom->natoms;
   double err;
   int kxmax = 1;
@@ -297,7 +289,8 @@ void EwaldDisp::setup()
    compute RMS accuracy for a dimension
 ------------------------------------------------------------------------- */
 
-double EwaldDisp::rms(int km, double prd, bigint natoms, double q2, double b2, double M2)
+double EwaldDisp::rms(int km, double prd, bigint natoms, 
+                      double q2, double b2, double M2)
 {
   double value = 0.0;
 
@@ -520,11 +513,18 @@ void EwaldDisp::init_coeff_sums()
   Sum sum_local[EWALD_MAX_NSUMS];
 
   memset(sum_local, 0, EWALD_MAX_NSUMS*sizeof(Sum));
-  if (function[0]) {                                        // 1/r
-    double *q = atom->q, *qn = q+atom->nlocal;
-    for (double *i=q; i<qn; ++i) {
-      sum_local[0].x += i[0]; sum_local[0].x2 += i[0]*i[0]; }
-  }
+
+  // now perform qsum and qsq via parent qsum_qsq()
+
+  sum_local[0].x = 0.0;
+  sum_local[0].x2 = 0.0;
+
+  //if (function[0]) {                                        // 1/r
+  //  double *q = atom->q, *qn = q+atom->nlocal;
+  //  for (double *i=q; i<qn; ++i) {
+  //    sum_local[0].x += i[0]; sum_local[0].x2 += i[0]*i[0]; }
+  //}
+
   if (function[1]) {                                        // geometric 1/r^6
     int *type = atom->type, *ntype = type+atom->nlocal;
     for (int *i=type; i<ntype; ++i) {
@@ -557,8 +557,8 @@ void EwaldDisp::init_self()
   memset(virial_self, 0, EWALD_NFUNCS*sizeof(double));
 
   if (function[0]) {                                        // 1/r
-    virial_self[0] = -0.5*MY_PI*qscale/(g2*volume)*sum[0].x*sum[0].x;
-    energy_self[0] = sum[0].x2*qscale*g1/MY_PIS-virial_self[0];
+    virial_self[0] = -0.5*MY_PI*qscale/(g2*volume)*qsum*qsum;
+    energy_self[0] = qsqsum*qscale*g1/MY_PIS-virial_self[0];
   }
   if (function[1]) {                                        // geometric 1/r^6
     virial_self[1] = MY_PI*MY_PIS*g3/(6.0*volume)*sum[1].x*sum[1].x;
@@ -597,7 +597,7 @@ void EwaldDisp::init_self_peratom()
     double *qi = atom->q, *qn = qi + nlocal;
     for (; qi < qn; qi++, vi += EWALD_NFUNCS, ei += EWALD_NFUNCS) {
       double q = *qi;
-      *vi = cv*q*sum[0].x;
+      *vi = cv*q*qsum;
       *ei = ce*q*q-vi[0];
     }
   }
@@ -676,6 +676,14 @@ void EwaldDisp::compute(int eflag, int vflag)
   compute_ek();
   compute_force();
   //compute_surface(); // assume conducting metal (tinfoil) boundary conditions
+
+  // update qsum and qsqsum, if atom count has changed and energy needed
+
+  if ((eflag_global || eflag_atom) && atom->natoms != natoms_original) {
+    qsum_qsq();
+    natoms_original = atom->natoms;
+  }
+
   compute_energy();
   compute_energy_peratom();
   compute_virial();
@@ -1196,7 +1204,6 @@ void EwaldDisp::compute_virial_dipole()
     for (int n = 0; n < 6; n++)
       virial[n] += sum[n];
   }
-
 }
 
 void EwaldDisp::compute_virial_peratom()
@@ -1325,7 +1332,6 @@ void EwaldDisp::compute_virial_peratom()
   }
 }
 
-
 /* ----------------------------------------------------------------------
    Slab-geometry correction term to dampen inter-slab interactions between
    periodically repeating slabs.  Yields good approximation to 2D Ewald if
@@ -1343,9 +1349,6 @@ void EwaldDisp::compute_slabcorr()
   double zprd = domain->zprd;
   int nlocal = atom->nlocal;
 
-  double qsum = 0.0;
-  if (function[0]) qsum = sum[0].x;
-
   double dipole = 0.0;
   for (int i = 0; i < nlocal; i++) dipole += q[i]*x[i][2];
 
@@ -1401,7 +1404,8 @@ void EwaldDisp::compute_slabcorr()
   double ffact = qscale * (-4.0*MY_PI/volume);
   double **f = atom->f;
 
-  for (int i = 0; i < nlocal; i++) f[i][2] += ffact * q[i]*(dipole_all - qsum*x[i][2]);
+  for (int i = 0; i < nlocal; i++) 
+    f[i][2] += ffact * q[i]*(dipole_all - qsum*x[i][2]);
 
   // add on torque corrections
 
@@ -1416,8 +1420,8 @@ void EwaldDisp::compute_slabcorr()
 }
 
 /* ----------------------------------------------------------------------
-  Newton solver used to find g_ewald for LJ systems
- ------------------------------------------------------------------------- */
+   Newton solver used to find g_ewald for LJ systems
+------------------------------------------------------------------------- */
 
 double EwaldDisp::NewtonSolve(double x, double Rc, 
                               bigint natoms, double vol, double b2)
diff --git a/src/KSPACE/msm.cpp b/src/KSPACE/msm.cpp
index 1f97100285..d3c72306be 100644
--- a/src/KSPACE/msm.cpp
+++ b/src/KSPACE/msm.cpp
@@ -100,6 +100,7 @@ MSM::MSM(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg)
 
   peratom_allocate_flag = 0;
   scalar_pressure_flag = 1;
+  warn_nonneutral = 0;
 
   order = 10;
 }
@@ -183,7 +184,7 @@ void MSM::init()
 
   scale = 1.0;
   qqrd2e = force->qqrd2e;
-  qsum_qsq(1);
+  qsum_qsq();
   natoms_original = atom->natoms;
 
   // set accuracy (force units) from accuracy_relative or accuracy_absolute
@@ -507,7 +508,6 @@ void MSM::compute(int eflag, int vflag)
     restriction(n);
   }
 
-
   // compute direct interation for top grid level for nonperiodic
   //   and for second from top grid level for periodic
 
@@ -565,13 +565,11 @@ void MSM::compute(int eflag, int vflag)
 
   if (evflag_atom) fieldforce_peratom();
 
-  // update qsum and qsqsum, if needed
+  // update qsum and qsqsum, if atom count has changed and energy needed
 
-  if (eflag_global || eflag_atom) {
-    if (qsum_update_flag || (atom->natoms != natoms_original)) {
-      qsum_qsq(0);
-      natoms_original = atom->natoms;
-    }
+  if ((eflag_global || eflag_atom) && atom->natoms != natoms_original) {
+    qsum_qsq();
+    natoms_original = atom->natoms;
   }
 
   // sum global energy across procs and add in self-energy term
diff --git a/src/KSPACE/msm_cg.cpp b/src/KSPACE/msm_cg.cpp
index 7c1c0fc448..30dcbeb3f6 100644
--- a/src/KSPACE/msm_cg.cpp
+++ b/src/KSPACE/msm_cg.cpp
@@ -247,13 +247,11 @@ void MSMCG::compute(int eflag, int vflag)
 
   if (evflag_atom) fieldforce_peratom();
 
-  // update qsum and qsqsum, if needed
+  // update qsum and qsqsum, if atom count has changed and energy needed
 
-  if (eflag_global || eflag_atom) {
-    if (qsum_update_flag || (atom->natoms != natoms_original)) {
-      qsum_qsq(0);
-      natoms_original = atom->natoms;
-    }
+  if ((eflag_global || eflag_atom) && atom->natoms != natoms_original) {
+    qsum_qsq();
+    natoms_original = atom->natoms;
   }
 
   // sum global energy across procs and add in self-energy term
@@ -301,7 +299,6 @@ void MSMCG::compute(int eflag, int vflag)
       }
     }
   }
-
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/KSPACE/pppm.cpp b/src/KSPACE/pppm.cpp
index 7a941d841c..814ea1e63e 100644
--- a/src/KSPACE/pppm.cpp
+++ b/src/KSPACE/pppm.cpp
@@ -257,7 +257,7 @@ void PPPM::init()
 
   scale = 1.0;
   qqrd2e = force->qqrd2e;
-  qsum_qsq(0);
+  qsum_qsq();
   natoms_original = atom->natoms;
 
   // set accuracy (force units) from accuracy_relative or accuracy_absolute
@@ -666,13 +666,11 @@ void PPPM::compute(int eflag, int vflag)
 
   if (evflag_atom) fieldforce_peratom();
 
-  // update qsum and qsqsum, if needed
+  // update qsum and qsqsum, if atom count has changed and energy needed
 
-  if (eflag_global || eflag_atom) {
-    if (qsum_update_flag || (atom->natoms != natoms_original)) {
-      qsum_qsq(0);
-      natoms_original = atom->natoms;
-    }
+  if ((eflag_global || eflag_atom) && atom->natoms != natoms_original) {
+    qsum_qsq();
+    natoms_original = atom->natoms;
   }
 
   // sum global energy across procs and add in volume-dependent term
diff --git a/src/KSPACE/pppm_cg.cpp b/src/KSPACE/pppm_cg.cpp
index 79cf7afe83..21dd09cae1 100644
--- a/src/KSPACE/pppm_cg.cpp
+++ b/src/KSPACE/pppm_cg.cpp
@@ -207,13 +207,11 @@ void PPPMCG::compute(int eflag, int vflag)
 
   if (evflag_atom) fieldforce_peratom();
 
-  // update qsum and qsqsum, if needed
+  // update qsum and qsqsum, if atom count has changed and energy needed
 
-  if (eflag_global || eflag_atom) {
-    if (qsum_update_flag || (atom->natoms != natoms_original)) {
-      qsum_qsq(0);
-      natoms_original = atom->natoms;
-    }
+  if ((eflag_global || eflag_atom) && atom->natoms != natoms_original) {
+    qsum_qsq();
+    natoms_original = atom->natoms;
   }
 
   // sum global energy across procs and add in volume-dependent term
diff --git a/src/KSPACE/pppm_disp.cpp b/src/KSPACE/pppm_disp.cpp
index d177cf7d94..d2dbe52f0c 100755
--- a/src/KSPACE/pppm_disp.cpp
+++ b/src/KSPACE/pppm_disp.cpp
@@ -293,7 +293,7 @@ void PPPMDisp::init()
   qqrd2e = force->qqrd2e;
   natoms_original = atom->natoms;
  
-  if (function[0]) qsum_qsq(0);
+  if (function[0]) qsum_qsq();
 
   // if kspace is TIP4P, extract TIP4P params from pair style
   // bond/angle are not yet init(), so insure equilibrium request is valid
@@ -1114,19 +1114,17 @@ void PPPMDisp::compute(int eflag, int vflag)
 
       fieldforce_none_ik();
 
-
-      if (evflag_atom) cg_peratom_6->forward_comm(this, FORWARD_IK_PERATOM_NONE);
+      if (evflag_atom) 
+        cg_peratom_6->forward_comm(this, FORWARD_IK_PERATOM_NONE);
     }
     if (evflag_atom) fieldforce_none_peratom();
   }
 
-  // update qsum and qsqsum, if needed
+  // update qsum and qsqsum, if atom count has changed and energy needed
 
-  if (eflag_global || eflag_atom) {
-    if (qsum_update_flag || (atom->natoms != natoms_original)) {
-      qsum_qsq(0);
-      natoms_original = atom->natoms;
-    }
+  if ((eflag_global || eflag_atom) && atom->natoms != natoms_original) {
+    qsum_qsq();
+    natoms_original = atom->natoms;
   }
 
   // sum energy across procs and add in volume-dependent term
diff --git a/src/KSPACE/pppm_stagger.cpp b/src/KSPACE/pppm_stagger.cpp
index 8ba10efd19..0e317fd64f 100755
--- a/src/KSPACE/pppm_stagger.cpp
+++ b/src/KSPACE/pppm_stagger.cpp
@@ -201,13 +201,11 @@ void PPPMStagger::compute(int eflag, int vflag)
     stagger += 1.0/float(nstagger);
   }
 
-  // update qsum and qsqsum, if needed
+  // update qsum and qsqsum, if atom count has changed and energy needed
 
-  if (eflag_global || eflag_atom) {
-    if (qsum_update_flag || (atom->natoms != natoms_original)) {
-      qsum_qsq(0);
-      natoms_original = atom->natoms;
-    }
+  if ((eflag_global || eflag_atom) && atom->natoms != natoms_original) {
+    qsum_qsq();
+    natoms_original = atom->natoms;
   }
 
   // sum global energy across procs and add in volume-dependent term
diff --git a/src/QEQ/fix_qeq.cpp b/src/QEQ/fix_qeq.cpp
index 2e735ef10c..999b9f371e 100644
--- a/src/QEQ/fix_qeq.cpp
+++ b/src/QEQ/fix_qeq.cpp
@@ -107,7 +107,6 @@ FixQEq::FixQEq(LAMMPS *lmp, int narg, char **arg) :
       s_hist[i][j] = t_hist[i][j] = 0;
 
   read_file(arg[7]);
-
 }
 
 /* ---------------------------------------------------------------------- */
@@ -262,7 +261,6 @@ void FixQEq::reallocate_matrix()
 void FixQEq::init_list(int id, NeighList *ptr)
 {
   list = ptr;
-  if (force->kspace) force->kspace->qsum_update_flag = 1;
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/QEQ/fix_qeq_dynamic.cpp b/src/QEQ/fix_qeq_dynamic.cpp
index 6099560245..f52f9e0c54 100644
--- a/src/QEQ/fix_qeq_dynamic.cpp
+++ b/src/QEQ/fix_qeq_dynamic.cpp
@@ -64,7 +64,8 @@ FixQEqDynamic::FixQEqDynamic(LAMMPS *lmp, int narg, char **arg) :
 
 void FixQEqDynamic::init()
 {
-  if (!atom->q_flag) error->all(FLERR,"Fix qeq/dynamic requires atom attribute q");
+  if (!atom->q_flag)
+    error->all(FLERR,"Fix qeq/dynamic requires atom attribute q");
 
   ngroup = group->count(igroup);
   if (ngroup == 0) error->all(FLERR,"Fix qeq/dynamic group has no atoms");
@@ -168,8 +169,7 @@ void FixQEqDynamic::pre_force(int vflag)
     }
   }
 
-  if (force->kspace) force->kspace->setup();
-
+  if (force->kspace) force->kspace->qsum_qsq();
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/QEQ/fix_qeq_point.cpp b/src/QEQ/fix_qeq_point.cpp
index 28499b5f8d..37194ba25a 100644
--- a/src/QEQ/fix_qeq_point.cpp
+++ b/src/QEQ/fix_qeq_point.cpp
@@ -45,7 +45,8 @@ FixQEqPoint::FixQEqPoint(LAMMPS *lmp, int narg, char **arg) :
 
 void FixQEqPoint::init()
 {
-  if (!atom->q_flag) error->all(FLERR,"Fix qeq/point requires atom attribute q");
+  if (!atom->q_flag) 
+    error->all(FLERR,"Fix qeq/point requires atom attribute q");
 
   ngroup = group->count(igroup);
   if (ngroup == 0) error->all(FLERR,"Fix qeq/point group has no atoms");
@@ -83,8 +84,7 @@ void FixQEqPoint::pre_force(int vflag)
   matvecs += CG(b_t, t); 	// CG on t - parallel
   calculate_Q();
 
-  if (force->kspace) force->kspace->setup();
-
+  if (force->kspace) force->kspace->qsum_qsq();
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/QEQ/fix_qeq_shielded.cpp b/src/QEQ/fix_qeq_shielded.cpp
index d4b061575c..0b6bfed0ef 100644
--- a/src/QEQ/fix_qeq_shielded.cpp
+++ b/src/QEQ/fix_qeq_shielded.cpp
@@ -45,7 +45,8 @@ FixQEqShielded::FixQEqShielded(LAMMPS *lmp, int narg, char **arg) :
 
 void FixQEqShielded::init()
 {
-  if (!atom->q_flag) error->all(FLERR,"Fix qeq/shielded requires atom attribute q");
+  if (!atom->q_flag) 
+    error->all(FLERR,"Fix qeq/shielded requires atom attribute q");
 
   ngroup = group->count(igroup);
   if (ngroup == 0) error->all(FLERR,"Fix qeq/shielded group has no atoms");
@@ -63,7 +64,8 @@ void FixQEqShielded::init()
 
   int i;
   for (i = 1; i <= ntypes; i++) {
-    if (gamma[i] == 0.0) error->all(FLERR,"Invalid param file for fix qeq/shielded");
+    if (gamma[i] == 0.0) 
+      error->all(FLERR,"Invalid param file for fix qeq/shielded");
   }
 
   if (strstr(update->integrate_style,"respa"))
@@ -126,8 +128,7 @@ void FixQEqShielded::pre_force(int vflag)
   matvecs += CG(b_t, t); 	// CG on t - parallel
   calculate_Q();
 
-  if (force->kspace) force->kspace->setup();
-
+  if (force->kspace) force->kspace->qsum_qsq();
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/QEQ/fix_qeq_slater.cpp b/src/QEQ/fix_qeq_slater.cpp
index dbae687597..acf3cc9f65 100644
--- a/src/QEQ/fix_qeq_slater.cpp
+++ b/src/QEQ/fix_qeq_slater.cpp
@@ -62,7 +62,8 @@ FixQEqSlater::FixQEqSlater(LAMMPS *lmp, int narg, char **arg) :
 
 void FixQEqSlater::init()
 {
-  if (!atom->q_flag) error->all(FLERR,"Fix qeq/slater requires atom attribute q");
+  if (!atom->q_flag) 
+    error->all(FLERR,"Fix qeq/slater requires atom attribute q");
 
   ngroup = group->count(igroup);
   if (ngroup == 0) error->all(FLERR,"Fix qeq/slater group has no atoms");
@@ -75,7 +76,8 @@ void FixQEqSlater::init()
 
   int ntypes = atom->ntypes;
   for (int i = 1; i <= ntypes; i++) {
-    if (zeta[i] == 0.0) error->all(FLERR,"Invalid param file for fix qeq/slater");
+    if (zeta[i] == 0.0) 
+      error->all(FLERR,"Invalid param file for fix qeq/slater");
   }
 
   if (strstr(update->integrate_style,"respa"))
@@ -101,8 +103,7 @@ void FixQEqSlater::pre_force(int vflag)
   matvecs += CG(b_t, t); 	// CG on t - parallel
   calculate_Q();
 
-  if (force->kspace) force->kspace->setup();
-
+  if (force->kspace) force->kspace->qsum_qsq();
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/USER-CUDA/pppm_cuda.cpp b/src/USER-CUDA/pppm_cuda.cpp
index e65c88fc2f..bb531f751e 100644
--- a/src/USER-CUDA/pppm_cuda.cpp
+++ b/src/USER-CUDA/pppm_cuda.cpp
@@ -691,7 +691,7 @@ void PPPMCuda::compute(int eflag, int vflag)
   // extend size of per-atom arrays if necessary
   // and force update of device data, if charges can change
 
-  if ((cu_atom->update_nmax)||(old_nmax==0)||qsum_update_flag) {
+  if ((cu_atom->update_nmax)||(old_nmax==0)||qdynamic) {
     memory->destroy(part2grid);
     nmax = atom->nmax;
     memory->create(part2grid,nmax,3,"pppm:part2grid");
diff --git a/src/USER-CUDA/pppm_old.cpp b/src/USER-CUDA/pppm_old.cpp
index be0a67040a..083e1d0db9 100755
--- a/src/USER-CUDA/pppm_old.cpp
+++ b/src/USER-CUDA/pppm_old.cpp
@@ -722,13 +722,11 @@ void PPPMOld::compute(int eflag, int vflag)
 
   if (evflag_atom) fieldforce_peratom();
 
-  // update qsum and qsqsum, if needed
+  // update qsum and qsqsum, if atom count has changed and energy needed
 
-  if (eflag_global || eflag_atom) {
-    if (qsum_update_flag || (atom->natoms != natoms_original)) {
-      qsum_qsq(0);
-      natoms_original = atom->natoms;
-    }
+  if ((eflag_global || eflag_atom) && atom->natoms != natoms_original) {
+    qsum_qsq(0);
+    natoms_original = atom->natoms;
   }
 
   // sum global energy across procs and add in volume-dependent term
diff --git a/src/USER-FEP/compute_fep.cpp b/src/USER-FEP/compute_fep.cpp
index 3bd5c186bb..06103f874b 100644
--- a/src/USER-FEP/compute_fep.cpp
+++ b/src/USER-FEP/compute_fep.cpp
@@ -412,15 +412,9 @@ void ComputeFEP::perturb_params()
 
   if (pairflag) force->pair->reinit();
 
-  // when perturbing charge and using kspace, 
-  // need to recompute additional params in kspace->setup()
-  // first backup the state of kspace->qsum_update_flag
+  // reset KSpace charges if charges have changed
 
-  if (chgflag && force->kspace) {
-    sys_qsum_update_flag = force->kspace->qsum_update_flag;
-    force->kspace->qsum_update_flag = 1;
-    force->kspace->setup();
-  }
+  if (chgflag && force->kspace) force->kspace->qsum_qsq();
 }
 
 
@@ -461,11 +455,10 @@ void ComputeFEP::restore_params()
   }
 
   if (pairflag) force->pair->reinit();
-  if (chgflag && force->kspace) {
-    force->kspace->setup();
-    // restore kspace->qsum_update_flag to original state
-    force->kspace->qsum_update_flag = sys_qsum_update_flag; 
-  }
+
+  // reset KSpace charges if charges have changed
+
+  if (chgflag && force->kspace) force->kspace->qsum_qsq();
 }
 
 
diff --git a/src/USER-FEP/fix_adapt_fep.cpp b/src/USER-FEP/fix_adapt_fep.cpp
index ffd7a6d99a..6d2efa44dc 100644
--- a/src/USER-FEP/fix_adapt_fep.cpp
+++ b/src/USER-FEP/fix_adapt_fep.cpp
@@ -172,11 +172,6 @@ FixAdaptFEP::FixAdaptFEP(LAMMPS *lmp, int narg, char **arg) :
     if (adapt[m].which == PAIR)
       memory->create(adapt[m].array_orig,n+1,n+1,"adapt:array_orig");
 
-  // when adapting charge and using kspace, 
-  // need to recompute additional params in kspace->setup()
-
-  if (chgflag && force->kspace) force->kspace->qsum_update_flag = 1;
-
   id_fix_diam = id_fix_chg = NULL;
 }
 
@@ -194,8 +189,6 @@ FixAdaptFEP::~FixAdaptFEP()
   }
   delete [] adapt;
 
-  if (chgflag && force->kspace) force->kspace->qsum_update_flag = 0;
-
   // check nfix in case all fixes have already been deleted
 
   if (id_fix_diam && modify->nfix) modify->delete_fix(id_fix_diam);
@@ -533,10 +526,9 @@ void FixAdaptFEP::change_settings()
 
   if (anypair) force->pair->reinit();
 
-  // re-setup KSpace if it exists and adapting charges
-  // since charges have changed
+  // reset KSpace charges if charges have changed
 
-  if (chgflag && force->kspace) force->kspace->setup();
+  if (chgflag && force->kspace) force->kspace->qsum_qsq();
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/USER-OMP/ewald_omp.cpp b/src/USER-OMP/ewald_omp.cpp
index 9d06f321a9..fb3470bd1e 100644
--- a/src/USER-OMP/ewald_omp.cpp
+++ b/src/USER-OMP/ewald_omp.cpp
@@ -86,14 +86,12 @@ void EwaldOMP::compute(int eflag, int vflag)
   MPI_Allreduce(sfacrl,sfacrl_all,kcount,MPI_DOUBLE,MPI_SUM,world);
   MPI_Allreduce(sfacim,sfacim_all,kcount,MPI_DOUBLE,MPI_SUM,world);
 
-  // update qsum and qsqsum, if needed
+  // update qsum and qsqsum, if atom count has changed and energy needed
   // (n.b. needs to be done outside of the multi-threaded region)
 
-  if (eflag_global || eflag_atom) {
-    if (qsum_update_flag || (atom->natoms != natoms_original)) {
-      qsum_qsq(0);
-      natoms_original = atom->natoms;
-    }
+  if ((eflag_global || eflag_atom) && atom->natoms != natoms_original) {
+    qsum_qsq(0);
+    natoms_original = atom->natoms;
   }
 
   // K-space portion of electric field
diff --git a/src/USER-OMP/msm_cg_omp.cpp b/src/USER-OMP/msm_cg_omp.cpp
index 340f736ff6..cb91dc8cc9 100644
--- a/src/USER-OMP/msm_cg_omp.cpp
+++ b/src/USER-OMP/msm_cg_omp.cpp
@@ -249,13 +249,11 @@ void MSMCGOMP::compute(int eflag, int vflag)
 
   if (evflag_atom) fieldforce_peratom();
 
-  // update qsum and qsqsum, if needed
+  // update qsum and qsqsum, if atom count has changed and energy needed
 
-  if (eflag_global || eflag_atom) {
-    if (qsum_update_flag || (atom->natoms != natoms_original)) {
-      qsum_qsq(0);
-      natoms_original = atom->natoms;
-    }
+  if ((eflag_global || eflag_atom) && atom->natoms != natoms_original) {
+    qsum_qsq(0);
+    natoms_original = atom->natoms;
   }
 
   // sum global energy across procs and add in self-energy term
diff --git a/src/fix_adapt.cpp b/src/fix_adapt.cpp
index f12f644de8..471485857e 100644
--- a/src/fix_adapt.cpp
+++ b/src/fix_adapt.cpp
@@ -175,8 +175,6 @@ FixAdapt::~FixAdapt()
   }
   delete [] adapt;
 
-  if (chgflag && force->kspace) force->kspace->qsum_update_flag = 0;
-
   // check nfix in case all fixes have already been deleted
 
   if (id_fix_diam && modify->nfix) modify->delete_fix(id_fix_diam);
@@ -353,11 +351,6 @@ void FixAdapt::init()
     }
   }
 
-  // when adapting charge and using kspace, 
-  // need to recompute additional params in kspace->setup()
-
-  if (chgflag && force->kspace) force->kspace->qsum_update_flag = 1;
-
   // fixes that store initial per-atom values
   
   if (id_fix_diam) {
@@ -504,10 +497,9 @@ void FixAdapt::change_settings()
 
   if (anypair) force->pair->reinit();
 
-  // re-setup KSpace if it exists and adapting charges
-  // since charges have changed
+  // reset KSpace charges if charges have changed
 
-  if (chgflag && force->kspace) force->kspace->setup();
+  if (chgflag && force->kspace) force->kspace->qsum_qsq();
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/kspace.cpp b/src/kspace.cpp
index 5942a5d0c1..4202b18f94 100644
--- a/src/kspace.cpp
+++ b/src/kspace.cpp
@@ -68,8 +68,7 @@ KSpace::KSpace(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
   suffix_flag = Suffix::NONE;
   adjust_cutoff_flag = 1;
   scalar_pressure_flag = 0;
-  qsum_update_flag = 0;
-  warn_neutral = 1;
+  warn_nonneutral = 1;
 
   accuracy_absolute = -1.0;
   accuracy_real_6 = -1.0;
@@ -259,10 +258,10 @@ void KSpace::ev_setup(int eflag, int vflag)
 
 /* ----------------------------------------------------------------------
    compute qsum,qsqsum,q2 and give error/warning if not charge neutral
-   only called initially or when particle count changes
+   called initially, when particle count changes, when charges are changed
 ------------------------------------------------------------------------- */
 
-void KSpace::qsum_qsq(int flag)
+void KSpace::qsum_qsq()
 {
   const double * const q = atom->q;
   const int nlocal = atom->nlocal;
@@ -285,15 +284,14 @@ void KSpace::qsum_qsq(int flag)
   q2 = qsqsum * force->qqrd2e;
 
   // not yet sure of the correction needed for non-neutral systems
+  // so issue warning or error
 
   if (fabs(qsum) > SMALL) {
     char str[128];
     sprintf(str,"System is not charge neutral, net charge = %g",qsum);
-    if (warn_neutral && (comm->me == 0)) {
-      if (flag) error->all(FLERR,str);
-      else error->warning(FLERR,str);
-    }
-    warn_neutral = 0;
+    if (!warn_nonneutral) error->all(FLERR,str);
+    if (warn_nonneutral == 1 && comm->me == 0) error->warning(FLERR,str);
+    warn_nonneutral = 2;
   }
 }
 
diff --git a/src/kspace.h b/src/kspace.h
index 09b8bc759d..32dd474cc6 100644
--- a/src/kspace.h
+++ b/src/kspace.h
@@ -51,9 +51,12 @@ class KSpace : protected Pointers {
                                  // for LJ coefficients
   int slabflag;
   int scalar_pressure_flag;      // 1 if using MSM fast scalar pressure
-  int qsum_update_flag;          // 1 if setup() needs to call qsum_qsq()
   double slab_volfactor;
 
+  int warn_nonneutral;           // 0 = error if non-neutral system
+                                 // 1 = warn once if non-neutral system
+                                 // 2 = warn, but already warned
+
   int order,order_6,order_allocated;
   double accuracy;                  // accuracy of KSpace solver (force units)
   double accuracy_absolute;         // user-specifed accuracy in force units
@@ -101,6 +104,10 @@ class KSpace : protected Pointers {
   void lamda2xvector(double *, double *);
   void kspacebbox(double, double *);
 
+  // public so can be called by commands that change charge
+
+  void qsum_qsq();
+
   // general child-class methods
 
   virtual void init() = 0;
@@ -168,7 +175,6 @@ class KSpace : protected Pointers {
   int minorder,overlap_allowed;
   int adjust_cutoff_flag;
   int suffix_flag;                  // suffix compatibility flag
-  int warn_neutral;                 // warn about non-neutral system if 1
   bigint natoms_original;
   double scale,qqrd2e;
   double qsum,qsqsum,q2;
@@ -182,7 +188,6 @@ class KSpace : protected Pointers {
   int kewaldflag;                   // 1 if kspace range set for Ewald sum
   int kx_ewald,ky_ewald,kz_ewald;   // kspace settings for Ewald sum
 
-  void qsum_qsq(int);
   void pair_check();
   void ev_setup(int, int);
   double estimate_table_accuracy(double, double);