diff --git a/src/QEQ/fix_qeq.cpp b/src/QEQ/fix_qeq.cpp
index c2aab73c358d85cdce4d0c75a5b4f22222ae485d..2e735ef10cb926bdb8d95b31749789ab52be3eb3 100644
--- a/src/QEQ/fix_qeq.cpp
+++ b/src/QEQ/fix_qeq.cpp
@@ -269,6 +269,9 @@ void FixQEq::init_list(int id, NeighList *ptr)
 
 void FixQEq::setup_pre_force(int vflag)
 {
+  if (force->newton_pair == 0)
+    error->all(FLERR,"QEQ with 'newton pair off' not supported");
+
   neighbor->build_one(list);
 
   deallocate_storage();
@@ -405,11 +408,9 @@ void FixQEq::sparse_matvec( sparse_matrix *A, double *x, double *b )
 {
   int i, j, itr_j;
   int nn, NN;
-  int *ilist;
 
   nn = atom->nlocal;
   NN = atom->nlocal + atom->nghost;
-  ilist = list->ilist;
 
   for( i = 0; i < nn; ++i ) {
     if (atom->mask[i] & groupbit)
@@ -696,7 +697,7 @@ void FixQEq::vector_add( double* dest, double c, double* v, int k )
 
 void FixQEq::read_file(char *file)
 {
-  int i,itype,ntypes;
+  int itype,ntypes;
   int params_per_line = 6;
   char **words = new char*[params_per_line+1];
 
diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp
index 1e44b345089add6154516ebcbb09077c676d7091..b014f44ba29d4942a156fa7b6ea3b9afa28311fb 100644
--- a/src/RIGID/fix_rigid_small.cpp
+++ b/src/RIGID/fix_rigid_small.cpp
@@ -1734,6 +1734,7 @@ void FixRigidSmall::setup_bodies_static()
   // dorientflag = 1 if any particle stores dipole orientation
 
   if (extended) {
+    grow_arrays(atom->nmax);
     if (atom->ellipsoid_flag) orientflag = 4;
     if (atom->line_flag) orientflag = 1;
     if (atom->tri_flag) orientflag = 4;
diff --git a/src/USER-FEP/compute_fep.cpp b/src/USER-FEP/compute_fep.cpp
index 847502cabc17dcb99678a79415668aec7056a96f..3bd5c186bbcb74e040248d1ee4740e4d4c1fa8ba 100644
--- a/src/USER-FEP/compute_fep.cpp
+++ b/src/USER-FEP/compute_fep.cpp
@@ -28,6 +28,8 @@
 #include "pair_hybrid.h"
 #include "kspace.h"
 #include "input.h"
+#include "fix.h"
+#include "modify.h"
 #include "variable.h"
 #include "timer.h"
 #include "memory.h"
@@ -39,9 +41,6 @@ using namespace LAMMPS_NS;
 enum{PAIR,ATOM};
 enum{CHARGE};
 
-#undef FEP_DEBUG
-#undef FEP_MAXDEBUG
-
 /* ---------------------------------------------------------------------- */
 
 ComputeFEP::ComputeFEP(LAMMPS *lmp, int narg, char **arg) :
@@ -66,11 +65,13 @@ ComputeFEP::ComputeFEP(LAMMPS *lmp, int narg, char **arg) :
   int iarg = 4;
   while (iarg < narg) {
     if (strcmp(arg[iarg],"pair") == 0) {
-      if (iarg+6 > narg) error->all(FLERR,"Illegal pair attribute in compute fep");
+      if (iarg+6 > narg) error->all(FLERR,
+                                    "Illegal pair attribute in compute fep");
       npert++;
       iarg += 6;
     } else if (strcmp(arg[iarg],"atom") == 0) {
-      if (iarg+4 > narg) error->all(FLERR,"Illegal atom attribute in compute fep");
+      if (iarg+4 > narg) error->all(FLERR,
+                                    "Illegal atom attribute in compute fep");
       npert++;
       iarg += 4;
     } else break;
@@ -159,6 +160,7 @@ ComputeFEP::ComputeFEP(LAMMPS *lmp, int narg, char **arg) :
 
   allocate_storage();
 
+  fixgpu = NULL;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -244,6 +246,11 @@ void ComputeFEP::init()
                  "compute tail corrections");
   }
 
+  // detect if package gpu is present
+
+  int ifixgpu = modify->find_fix("package_gpu");
+  if (ifixgpu >= 0) fixgpu = modify->fix[ifixgpu];
+
   if (comm->me == 0) {
     if (screen) {
       fprintf(screen, "FEP settings ...\n");
@@ -282,6 +289,9 @@ void ComputeFEP::compute_vector()
 {
   double pe0,pe1;
 
+  eflag = 1;
+  vflag = 0;
+
   invoked_vector = update->ntimestep;
 
   if (atom->nmax > nmax) {  // reallocate working arrays if necessary
@@ -294,26 +304,36 @@ void ComputeFEP::compute_vector()
 
   timer->stamp();
   if (force->pair && force->pair->compute_flag) {
-    force->pair->compute(1,0);
+    force->pair->compute(eflag,vflag);
     timer->stamp(TIME_PAIR);
   }
-  if (force->kspace && force->kspace->compute_flag) {
-    force->kspace->compute(1,0);
+  if (chgflag && force->kspace && force->kspace->compute_flag) {
+    force->kspace->compute(eflag,vflag);
     timer->stamp(TIME_KSPACE);
   }
+
+  // accumulate force/energy/virial from /gpu pair styles
+  if (fixgpu) fixgpu->post_force(vflag);
+
   pe0 = compute_epair();
 
   perturb_params();
 
   timer->stamp();
   if (force->pair && force->pair->compute_flag) {
-    force->pair->compute(1,0);
+    force->pair->compute(eflag,vflag);
     timer->stamp(TIME_PAIR);
   }
-  if (force->kspace && force->kspace->compute_flag) {
-    force->kspace->compute(1,0);
+  if (chgflag && force->kspace && force->kspace->compute_flag) {
+    force->kspace->compute(eflag,vflag);
     timer->stamp(TIME_KSPACE);
   }
+
+  // accumulate force/energy/virial from /gpu pair styles
+  // this is required as to empty the answer queue, 
+  // otherwise the force compute on the GPU in the next step would be incorrect
+  if (fixgpu) fixgpu->post_force(vflag);
+    
   pe1 = compute_epair();
 
   restore_qfev();   // restore charge, force, energy, virial array values
@@ -324,12 +344,6 @@ void ComputeFEP::compute_vector()
   vector[2] = domain->xprd * domain->yprd * domain->zprd;
   if (volumeflag)
     vector[1] *= vector[2];
-
-#ifdef FEP_DEBUG
-  if (comm->me == 0 && screen)
-    fprintf(screen, "FEP U0 = %f  U1 = %f  DU = %f  exp(-DU/kT) = %f\n",
-            pe0,pe1,vector[0],vector[1]);
-#endif
 }
 
 
@@ -351,7 +365,7 @@ double ComputeFEP::compute_epair()
     eng_pair += force->pair->etail / volume;
   }
 
-  if (force->kspace) eng_pair += force->kspace->energy;
+  if (chgflag && force->kspace) eng_pair += force->kspace->energy;
 
   return eng_pair;
 }
@@ -374,18 +388,6 @@ void ComputeFEP::perturb_params()
       for (i = pert->ilo; i <= pert->ihi; i++)
         for (j = MAX(pert->jlo,i); j <= pert->jhi; j++)
           pert->array[i][j] = pert->array_orig[i][j] + delta;
-      
-#ifdef FEP_MAXDEBUG
-      if (comm->me == 0 && screen) {
-        fprintf(screen, "###FEP change %s %s, delta = %f\n",
-                pert->pstyle, pert->pparam, delta);
-        fprintf(screen, "###FEP  I  J   old_param new_param\n");
-        for (i = pert->ilo; i <= pert->ihi; i++)
-          for (j = MAX(pert->jlo,i); j <= pert->jhi; j++)
-            fprintf(screen, "###FEP %2d %2d %9.5f %9.5f\n", i, j, 
-                    pert->array_orig[i][j], pert->array[i][j]);
-      }
-#endif
 
     } else if (pert->which == ATOM) {
 
@@ -400,18 +402,6 @@ void ComputeFEP::perturb_params()
             if (mask[i] & groupbit)
               q[i] += delta; 
 
-#ifdef FEP_MAXDEBUG
-        if (comm->me == 0 && screen) {
-          fprintf(screen, "###FEP change charge, delta = %f\n", delta);
-          fprintf(screen, "###FEP  atom  I   old_q     new_q\n");
-          for (i = 0; i < atom->nlocal; i++)
-            if (atype[i] >= pert->ilo && atype[i] <= pert->ihi)
-              if (mask[i] & groupbit)
-                fprintf(screen, "###FEP %5d %2d %9.5f %9.5f\n", i, atype[i],
-                        q_orig[i], q[i]);
-        }
-#endif
-
       }
     }
   }
@@ -421,6 +411,16 @@ void ComputeFEP::perturb_params()
   // and also offset and tail corrections
 
   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
+
+  if (chgflag && force->kspace) {
+    sys_qsum_update_flag = force->kspace->qsum_update_flag;
+    force->kspace->qsum_update_flag = 1;
+    force->kspace->setup();
+  }
 }
 
 
@@ -457,40 +457,15 @@ void ComputeFEP::restore_params()
       for (i = pert->ilo; i <= pert->ihi; i++)
         for (j = MAX(pert->jlo,i); j <= pert->jhi; j++)
           pert->array[i][j] = pert->array_orig[i][j];
-
-#ifdef FEP_MAXDEBUG
-      if (comm->me == 0 && screen) {
-        fprintf(screen, "###FEP restore %s %s\n", pert->pstyle, pert->pparam);
-        fprintf(screen, "###FEP  I  J   param\n");
-        for (i = pert->ilo; i <= pert->ihi; i++)
-          for (j = MAX(pert->jlo,i); j <= pert->jhi; j++)
-            fprintf(screen, "###FEP %2d %2d %9.5f\n", i, j, pert->array[i][j]);
-      }
-#endif
     }
-#ifdef FEP_MAXDEBUG
-    if (pert->which == ATOM) {
-      if (comm->me == 0 && screen) {
-        int *atype = atom->type;
-        double *q = atom->q; 
-        int *mask = atom->mask;
-        int natom = atom->nlocal;
-        fprintf(screen, "###FEP restore charge\n");
-        fprintf(screen, "###FEP  atom  I   q\n");
-        for (i = 0; i < natom; i++)
-          if (atype[i] >= pert->ilo && atype[i] <= pert->ihi)
-            if (mask[i] & groupbit)
-              fprintf(screen, "###FEP %5d %2d %9.5f\n", i, atype[i], q[i]);
-      }
-    }
-#endif
   }
 
-  // re-initialize pair styles if any PAIR settings were changed
-  // this resets other coeffs that may depend on changed values,
-  // and also offset and tail corrections
-
   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; 
+  }
 }
 
 
@@ -501,14 +476,15 @@ void ComputeFEP::restore_params()
 void ComputeFEP::allocate_storage()
 {
   nmax = atom->nmax;
-  if (chgflag)
-    memory->create(q_orig,nmax,"fep:q_orig");
   memory->create(f_orig,nmax,3,"fep:f_orig");
   memory->create(peatom_orig,nmax,"fep:peatom_orig");
   memory->create(pvatom_orig,nmax,6,"fep:pvatom_orig");
-  if (force->kspace) {
-    memory->create(keatom_orig,nmax,"fep:keatom_orig");
-    memory->create(kvatom_orig,nmax,6,"fep:kvatom_orig");
+  if (chgflag) {
+    memory->create(q_orig,nmax,"fep:q_orig");
+    if (force->kspace) {
+      memory->create(keatom_orig,nmax,"fep:keatom_orig");
+      memory->create(kvatom_orig,nmax,6,"fep:kvatom_orig");
+    }
   }
 }
 
@@ -516,14 +492,15 @@ void ComputeFEP::allocate_storage()
 
 void ComputeFEP::deallocate_storage()
 {
-  if (chgflag)
-    memory->destroy(q_orig);
   memory->destroy(f_orig);
   memory->destroy(peatom_orig);
   memory->destroy(pvatom_orig);
-  if (force->kspace) {
-    memory->destroy(keatom_orig);
-    memory->destroy(kvatom_orig);
+  if (chgflag) {
+    memory->destroy(q_orig);
+    if (force->kspace) {
+      memory->destroy(keatom_orig);
+      memory->destroy(kvatom_orig);
+    }
   }
 }
 
@@ -541,12 +518,6 @@ void ComputeFEP::backup_qfev()
   if (force->newton || force->kspace->tip4pflag)
       natom += atom->nghost;
 
-  if (chgflag) {
-    double *q = atom->q; 
-    for (i = 0; i < nall; i++)
-      q_orig[i] = q[i];
-  }
-
   double **f = atom->f;
   for (i = 0; i < natom; i++) {
     f_orig[i][0] = f[i][0]; 
@@ -581,29 +552,35 @@ void ComputeFEP::backup_qfev()
     }
   }
 
-  if (force->kspace) {
-    energy_orig = force->kspace->energy;
-    kvirial_orig[0] = force->kspace->virial[0];
-    kvirial_orig[1] = force->kspace->virial[1];
-    kvirial_orig[2] = force->kspace->virial[2];
-    kvirial_orig[3] = force->kspace->virial[3];
-    kvirial_orig[4] = force->kspace->virial[4];
-    kvirial_orig[5] = force->kspace->virial[5];
-    
-    if (update->eflag_atom) {
-      double *keatom = force->kspace->eatom;
-      for (i = 0; i < natom; i++)
-        keatom_orig[i] = keatom[i];
-    }
-    if (update->vflag_atom) {
-      double **kvatom = force->kspace->vatom;
-      for (i = 0; i < natom; i++) {
-        kvatom_orig[i][0] = kvatom[i][0];
-        kvatom_orig[i][1] = kvatom[i][1];
-        kvatom_orig[i][2] = kvatom[i][2];
-        kvatom_orig[i][3] = kvatom[i][3];
-        kvatom_orig[i][4] = kvatom[i][4];
-        kvatom_orig[i][5] = kvatom[i][5];
+  if (chgflag) {
+    double *q = atom->q; 
+    for (i = 0; i < nall; i++)
+      q_orig[i] = q[i];
+
+    if (force->kspace) {
+      energy_orig = force->kspace->energy;
+      kvirial_orig[0] = force->kspace->virial[0];
+      kvirial_orig[1] = force->kspace->virial[1];
+      kvirial_orig[2] = force->kspace->virial[2];
+      kvirial_orig[3] = force->kspace->virial[3];
+      kvirial_orig[4] = force->kspace->virial[4];
+      kvirial_orig[5] = force->kspace->virial[5];
+      
+      if (update->eflag_atom) {
+        double *keatom = force->kspace->eatom;
+        for (i = 0; i < natom; i++)
+          keatom_orig[i] = keatom[i];
+      }
+      if (update->vflag_atom) {
+        double **kvatom = force->kspace->vatom;
+        for (i = 0; i < natom; i++) {
+          kvatom_orig[i][0] = kvatom[i][0];
+          kvatom_orig[i][1] = kvatom[i][1];
+          kvatom_orig[i][2] = kvatom[i][2];
+          kvatom_orig[i][3] = kvatom[i][3];
+          kvatom_orig[i][4] = kvatom[i][4];
+          kvatom_orig[i][5] = kvatom[i][5];
+        }
       }
     }
   }
@@ -620,12 +597,6 @@ void ComputeFEP::restore_qfev()
   if (force->newton || force->kspace->tip4pflag)
       natom += atom->nghost;
 
-  if (chgflag) {
-    double *q = atom->q; 
-    for (i = 0; i < nall; i++)
-      q[i] = q_orig[i];
-  }
-
   double **f = atom->f;
   for (i = 0; i < natom; i++) {
     f[i][0] = f_orig[i][0]; 
@@ -660,29 +631,35 @@ void ComputeFEP::restore_qfev()
     }
   }
 
-  if (force->kspace) {
-    force->kspace->energy = energy_orig;
-    force->kspace->virial[0] = kvirial_orig[0];
-    force->kspace->virial[1] = kvirial_orig[1];
-    force->kspace->virial[2] = kvirial_orig[2];
-    force->kspace->virial[3] = kvirial_orig[3];
-    force->kspace->virial[4] = kvirial_orig[4];
-    force->kspace->virial[5] = kvirial_orig[5];
+  if (chgflag) {
+    double *q = atom->q; 
+    for (i = 0; i < nall; i++)
+      q[i] = q_orig[i];
+
+    if (force->kspace) {
+      force->kspace->energy = energy_orig;
+      force->kspace->virial[0] = kvirial_orig[0];
+      force->kspace->virial[1] = kvirial_orig[1];
+      force->kspace->virial[2] = kvirial_orig[2];
+      force->kspace->virial[3] = kvirial_orig[3];
+      force->kspace->virial[4] = kvirial_orig[4];
+      force->kspace->virial[5] = kvirial_orig[5];
     
-    if (update->eflag_atom) {
-      double *keatom = force->kspace->eatom;
-      for (i = 0; i < natom; i++)
-        keatom[i] = keatom_orig[i];
-    }
-    if (update->vflag_atom) {
-      double **kvatom = force->kspace->vatom;
-      for (i = 0; i < natom; i++) {
-        kvatom[i][0] = kvatom_orig[i][0];
-        kvatom[i][1] = kvatom_orig[i][1];
-        kvatom[i][2] = kvatom_orig[i][2];
-        kvatom[i][3] = kvatom_orig[i][3];
-        kvatom[i][4] = kvatom_orig[i][4];
-        kvatom[i][5] = kvatom_orig[i][5];
+      if (update->eflag_atom) {
+        double *keatom = force->kspace->eatom;
+        for (i = 0; i < natom; i++)
+          keatom[i] = keatom_orig[i];
+      }
+      if (update->vflag_atom) {
+        double **kvatom = force->kspace->vatom;
+        for (i = 0; i < natom; i++) {
+          kvatom[i][0] = kvatom_orig[i][0];
+          kvatom[i][1] = kvatom_orig[i][1];
+          kvatom[i][2] = kvatom_orig[i][2];
+          kvatom[i][3] = kvatom_orig[i][3];
+          kvatom[i][4] = kvatom_orig[i][4];
+          kvatom[i][5] = kvatom_orig[i][5];
+        }
       }
     }
   }
diff --git a/src/USER-FEP/compute_fep.h b/src/USER-FEP/compute_fep.h
index 7cc96034e7b9f89a3386d3cc5f1c919bf0c64704..fdbd4ed6e3d58ed6f2ddddcad907370a57a98f62 100644
--- a/src/USER-FEP/compute_fep.h
+++ b/src/USER-FEP/compute_fep.h
@@ -41,6 +41,8 @@ class ComputeFEP : public Compute {
   int chgflag;
   int tailflag, volumeflag;
   int fepinitflag;
+  int eflag, vflag;
+  int sys_qsum_update_flag;
   double temp_fep;
 
   int nmax;
@@ -53,6 +55,8 @@ class ComputeFEP : public Compute {
   double kvirial_orig[6];
   double *keatom_orig,**kvatom_orig;
 
+  class Fix *fixgpu;
+
   struct Perturb {
     int which,ivar;
     char *var;
diff --git a/src/USER-FEP/fix_adapt_fep.cpp b/src/USER-FEP/fix_adapt_fep.cpp
index 8cda1b69169767914a20e43da2fe0560781d295d..793b7a3793bb7b2351acdf4cc5d883d8c2b907d5 100644
--- a/src/USER-FEP/fix_adapt_fep.cpp
+++ b/src/USER-FEP/fix_adapt_fep.cpp
@@ -171,6 +171,11 @@ 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;
 }
 
@@ -188,7 +193,7 @@ FixAdaptFEP::~FixAdaptFEP()
   }
   delete [] adapt;
 
-  if (force->kspace) force->kspace->qsum_update_flag = 0;
+  if (chgflag && force->kspace) force->kspace->qsum_update_flag = 0;
 
   // check nfix in case all fixes have already been deleted
 
@@ -291,11 +296,6 @@ void FixAdaptFEP::init()
       if (adapt[i].which == ATOM) 
         error->all(FLERR,"Cannot use dynamic group with fix adapt/fep atom");
 
-  // when using kspace, we need to recompute 
-  // some additional parameters in kspace->setup()
-
-  if (force->kspace) force->kspace->qsum_update_flag = 1;
-  
   // setup and error checks
 
   anypair = 0;
@@ -455,6 +455,8 @@ void FixAdaptFEP::change_settings()
     } else if (ad->which == KSPACE) {
       *kspace_scale = value;
 
+    // set per atom values, also make changes for ghost atoms
+
     } else if (ad->which == ATOM) {
 
       // reset radius from diameter
@@ -510,13 +512,14 @@ void FixAdaptFEP::change_settings()
 
   if (anypair) force->pair->reinit();
 
-  // re-setup KSpace if using it, since charges may have changed
+  // re-setup KSpace if it exists and adapting charges
+  // since charges have changed
 
-  if (force->kspace) force->kspace->setup();
+  if (chgflag && force->kspace) force->kspace->setup();
 }
 
 /* ----------------------------------------------------------------------
-   restore pair,kspace,atom parameters to original values
+   restore pair,kspace.atom parameters to original values
 ------------------------------------------------------------------------- */
 
 void FixAdaptFEP::restore_settings()
@@ -565,5 +568,5 @@ void FixAdaptFEP::restore_settings()
   }
 
   if (anypair) force->pair->reinit();
-  if (force->kspace) force->kspace->setup();
+  if (chgflag && force->kspace) force->kspace->setup();
 }
diff --git a/src/VORONOI/compute_voronoi_atom.cpp b/src/VORONOI/compute_voronoi_atom.cpp
index f312476836d63cced8d37751d17b87a4ae25e598..ac3bcbda8d3d562c48dc207d8e195e9e5e072fad 100644
--- a/src/VORONOI/compute_voronoi_atom.cpp
+++ b/src/VORONOI/compute_voronoi_atom.cpp
@@ -15,6 +15,7 @@
    Contributing author: Daniel Schwen
 ------------------------------------------------------------------------- */
 
+#include "mpi.h"
 #include "math.h"
 #include "string.h"
 #include "stdlib.h"