diff --git a/src/KSPACE/ewald.cpp b/src/KSPACE/ewald.cpp
index b6cff3d932..4e2782ff04 100644
--- a/src/KSPACE/ewald.cpp
+++ b/src/KSPACE/ewald.cpp
@@ -95,7 +95,6 @@ void Ewald::init()
 
   // extract short-range Coulombic cutoff from pair style
 
-  qqrd2e = force->qqrd2e;
   scale = 1.0;
 
   if (force->pair == NULL)
@@ -270,10 +269,12 @@ void Ewald::compute(int eflag, int vflag)
 
   // convert E-field to force
 
+  const double qscale = force->qqrd2e * scale;
+
   for (i = 0; i < nlocal; i++) {
-    f[i][0] += qqrd2e*scale * q[i]*ek[i][0];
-    f[i][1] += qqrd2e*scale * q[i]*ek[i][1];
-    f[i][2] += qqrd2e*scale * q[i]*ek[i][2];
+    f[i][0] += qscale * q[i]*ek[i][0];
+    f[i][1] += qscale * q[i]*ek[i][1];
+    f[i][2] += qscale * q[i]*ek[i][2];
   }
  
   // energy if requested
@@ -284,7 +285,7 @@ void Ewald::compute(int eflag, int vflag)
 			 sfacim_all[k]*sfacim_all[k]);
     energy -= g_ewald*qsqsum/MY_PIS +
       MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);
-    energy *= qqrd2e*scale;
+    energy *= qscale;
   }
 
   // virial if requested
@@ -295,7 +296,7 @@ void Ewald::compute(int eflag, int vflag)
       uk = ug[k] * (sfacrl_all[k]*sfacrl_all[k] + sfacim_all[k]*sfacim_all[k]);
       for (n = 0; n < 6; n++) virial[n] += uk*vg[k][n];
     }
-    for (n = 0; n < 6; n++) virial[n] *= qqrd2e*scale;
+    for (n = 0; n < 6; n++) virial[n] *= qscale;
   }
 
   if (slabflag) slabcorr(eflag);
@@ -817,16 +818,17 @@ void Ewald::slabcorr(int eflag)
 
   // compute corrections
   
-  double e_slabcorr = 2.0*MY_PI*dipole_all*dipole_all/volume;
+  const double e_slabcorr = 2.0*MY_PI*dipole_all*dipole_all/volume;
+  const double qscale = force->qqrd2e * scale;
   
-  if (eflag) energy += qqrd2e*scale * e_slabcorr;
+  if (eflag) energy += qscale * e_slabcorr;
 
   // add on force corrections
 
   double ffact = -4.0*MY_PI*dipole_all/volume; 
   double **f = atom->f;
 
-  for (int i = 0; i < nlocal; i++) f[i][2] += qqrd2e*scale * q[i]*ffact;
+  for (int i = 0; i < nlocal; i++) f[i][2] += qscale * q[i]*ffact;
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/KSPACE/ewald.h b/src/KSPACE/ewald.h
index cd22f987f7..e88b8989ac 100644
--- a/src/KSPACE/ewald.h
+++ b/src/KSPACE/ewald.h
@@ -36,7 +36,6 @@ class Ewald : public KSpace {
  protected:
   double precision;
   int kcount,kmax,kmax3d,kmax_created;
-  double qqrd2e;
   double gsqmx,qsum,qsqsum,volume;
   int nmax;
 
diff --git a/src/KSPACE/pppm.cpp b/src/KSPACE/pppm.cpp
index e2d017a6e5..4e810bbc23 100644
--- a/src/KSPACE/pppm.cpp
+++ b/src/KSPACE/pppm.cpp
@@ -139,7 +139,6 @@ void PPPM::init()
 
   // extract short-range Coulombic cutoff from pair style
 
-  qqrd2e = force->qqrd2e;
   scale = 1.0;
 
   if (force->pair == NULL)
@@ -714,6 +713,8 @@ void PPPM::compute(int eflag, int vflag)
 
   // sum energy across procs and add in volume-dependent term
 
+  const double qscale = force->qqrd2e * scale;
+
   if (eflag) {
     double energy_all;
     MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
@@ -722,7 +723,7 @@ void PPPM::compute(int eflag, int vflag)
     energy *= 0.5*volume;
     energy -= g_ewald*qsqsum/MY_PIS +
       MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);
-    energy *= qqrd2e*scale;
+    energy *= qscale;
   }
 
   // sum virial across procs
@@ -730,7 +731,7 @@ void PPPM::compute(int eflag, int vflag)
   if (vflag) {
     double virial_all[6];
     MPI_Allreduce(virial,virial_all,6,MPI_DOUBLE,MPI_SUM,world);
-    for (i = 0; i < 6; i++) virial[i] = 0.5*qqrd2e*scale*volume*virial_all[i];
+    for (i = 0; i < 6; i++) virial[i] = 0.5*qscale*volume*virial_all[i];
   }
 
   // 2d slab correction
@@ -1734,7 +1735,7 @@ void PPPM::fieldforce()
     }
 
     // convert E-field to force
-    const double qfactor = qqrd2e*scale*q[i];
+    const double qfactor = force->qqrd2e * scale * q[i];
     f[i][0] += qfactor*ekx;
     f[i][1] += qfactor*eky;
     f[i][2] += qfactor*ekz;
@@ -1887,16 +1888,17 @@ void PPPM::slabcorr(int eflag)
 
   // compute corrections
   
-  double e_slabcorr = 2.0*MY_PI*dipole_all*dipole_all/volume;
+  const double e_slabcorr = 2.0*MY_PI*dipole_all*dipole_all/volume;
+  const double qscale = force->qqrd2e * scale;
   
-  if (eflag) energy += qqrd2e*scale * e_slabcorr;
+  if (eflag) energy += qscale * e_slabcorr;
 
   // add on force corrections
 
   double ffact = -4.0*MY_PI*dipole_all/volume; 
   double **f = atom->f;
 
-  for (int i = 0; i < nlocal; i++) f[i][2] += qqrd2e*scale * q[i]*ffact;
+  for (int i = 0; i < nlocal; i++) f[i][2] += qscale * q[i]*ffact;
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/KSPACE/pppm.h b/src/KSPACE/pppm.h
index e2a58e8804..1ee2cfbb4f 100644
--- a/src/KSPACE/pppm.h
+++ b/src/KSPACE/pppm.h
@@ -51,7 +51,6 @@ class PPPM : public KSpace {
   int nfactors;
   int *factors;
   double qsum,qsqsum;
-  double qqrd2e;
   double cutoff;
   double volume;
   double delxinv,delyinv,delzinv,delvolinv;
diff --git a/src/KSPACE/pppm_cg.cpp b/src/KSPACE/pppm_cg.cpp
index daffb5ff2d..663d263a20 100644
--- a/src/KSPACE/pppm_cg.cpp
+++ b/src/KSPACE/pppm_cg.cpp
@@ -23,6 +23,7 @@
 #include "atom.h"
 #include "domain.h"
 #include "error.h"
+#include "force.h"
 #include "memory.h"
 #include "pppm_cg.h"
 
@@ -172,6 +173,8 @@ void PPPMCG::compute(int eflag, int vflag)
 
   // sum energy across procs and add in volume-dependent term
 
+  const double qscale = force->qqrd2e * scale;
+
   if (eflag) {
     double energy_all;
     MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
@@ -180,7 +183,7 @@ void PPPMCG::compute(int eflag, int vflag)
     energy *= 0.5*volume;
     energy -= g_ewald*qsqsum/MY_PIS +
       MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);
-    energy *= qqrd2e*scale;
+    energy *= qscale;
   }
 
   // sum virial across procs
@@ -188,7 +191,7 @@ void PPPMCG::compute(int eflag, int vflag)
   if (vflag) {
     double virial_all[6];
     MPI_Allreduce(virial,virial_all,6,MPI_DOUBLE,MPI_SUM,world);
-    for (i = 0; i < 6; i++) virial[i] = 0.5*qqrd2e*scale*volume*virial_all[i];
+    for (i = 0; i < 6; i++) virial[i] = 0.5*qscale*volume*virial_all[i];
   }
 
   // 2d slab correction
@@ -341,7 +344,7 @@ void PPPMCG::fieldforce()
     }
 
     // convert E-field to force
-    const double qfactor = qqrd2e*scale*q[i];
+    const double qfactor = force->qqrd2e * scale * q[i];
     f[i][0] += qfactor*ekx;
     f[i][1] += qfactor*eky;
     f[i][2] += qfactor*ekz;
@@ -375,13 +378,14 @@ void PPPMCG::slabcorr(int eflag)
 
   // compute corrections
   
-  double e_slabcorr = 2.0*MY_PI*dipole_all*dipole_all/volume;
+  const double e_slabcorr = 2.0*MY_PI*dipole_all*dipole_all/volume;
+  const double qscale = force->qqrd2e * scale;
   
-  if (eflag) energy += qqrd2e*scale * e_slabcorr;
+  if (eflag) energy += qscale * e_slabcorr;
 
   // add on force corrections
 
-  double ffact = -4.0*MY_PI*dipole_all/volume * qqrd2e * scale; 
+  const double ffact = -4.0*MY_PI*dipole_all/volume * qscale; 
   double **f = atom->f;
 
   for (int j = 0; j < num_charged; j++) {
diff --git a/src/KSPACE/pppm_tip4p.cpp b/src/KSPACE/pppm_tip4p.cpp
index a191fa742c..16a0f7a04c 100644
--- a/src/KSPACE/pppm_tip4p.cpp
+++ b/src/KSPACE/pppm_tip4p.cpp
@@ -19,6 +19,7 @@
 #include "pppm_tip4p.h"
 #include "atom.h"
 #include "domain.h"
+#include "force.h"
 #include "memory.h"
 #include "error.h"
 
@@ -205,7 +206,7 @@ void PPPMTIP4P::fieldforce()
     }
 
     // convert E-field to force
-    const double qfactor = qqrd2e*scale*q[i];
+    const double qfactor = force->qqrd2e * scale * q[i];
     if (type[i] != typeO) {
       f[i][0] += qfactor*ekx;
       f[i][1] += qfactor*eky;
diff --git a/src/USER-EWALDN/ewald_n.cpp b/src/USER-EWALDN/ewald_n.cpp
index ed05ca5044..f3d7fe08d3 100644
--- a/src/USER-EWALDN/ewald_n.cpp
+++ b/src/USER-EWALDN/ewald_n.cpp
@@ -87,7 +87,6 @@ void EwaldN::init()
       error->all(FLERR,"Incorrect boundaries with slab EwaldN");
   }
 
-  qqrd2e = force->qqrd2e;				// check pair_style
   scale = 1.0;
   //mumurd2e = force->mumurd2e;
   //dielectric = force->dielectric;
@@ -368,9 +367,11 @@ void EwaldN::init_self()
 
   memset(energy_self, 0, EWALD_NFUNCS*sizeof(double));	// self energy
   memset(virial_self, 0, EWALD_NFUNCS*sizeof(double));
+  const double qscale = force->qqrd2e * scale;
+
   if (function[0]) {					// 1/r
-    virial_self[0] = -0.5*M_PI*qqrd2e*scale/(g2*volume)*sum[0].x*sum[0].x;
-    energy_self[0] = sum[0].x2*qqrd2e*scale*g1/sqrt(M_PI)-virial_self[0];
+    virial_self[0] = -0.5*M_PI*qscale/(g2*volume)*sum[0].x*sum[0].x;
+    energy_self[0] = sum[0].x2*qscale*g1/sqrt(M_PI)-virial_self[0];
   }
   if (function[1]) {					// geometric 1/r^6
     virial_self[1] = M_PI*sqrt(M_PI)*g3/(6.0*volume)*sum[1].x*sum[1].x;
@@ -484,8 +485,9 @@ void EwaldN::compute_force()
   complex *cek, zc, zx, zxy;
   double *f = atom->f[0], *fn = f+3*atom->nlocal, *q = atom->q, *t = NULL;
   double *mu = atom->mu ? atom->mu[0] : NULL;
+  const double qscale = force->qqrd2e * scale;
   double *ke, c[EWALD_NFUNCS] = {
-    8.0*M_PI*qqrd2e*scale/volume, 2.0*M_PI*sqrt(M_PI)/(12.0*volume),
+    8.0*M_PI*qscale/volume, 2.0*M_PI*sqrt(M_PI)/(12.0*volume),
     2.0*M_PI*sqrt(M_PI)/(192.0*volume), 8.0*M_PI*mumurd2e/volume};
   double kt = 4.0*pow(g_ewald, 3.0)/3.0/sqrt(M_PI)/c[3];
   int i, kx, ky, lbytes = (2*nbox+1)*sizeof(cvector), *type = atom->type;
@@ -587,8 +589,9 @@ void EwaldN::compute_energy(int eflag)
   
   complex *cek = cek_global;
   double *ke = kenergy;
+  const double qscale = force->qqrd2e * scale;
   double c[EWALD_NFUNCS] = {
-    4.0*M_PI*qqrd2e*scale/volume, 2.0*M_PI*sqrt(M_PI)/(24.0*volume),
+    4.0*M_PI*qscale/volume, 2.0*M_PI*sqrt(M_PI)/(24.0*volume),
     2.0*M_PI*sqrt(M_PI)/(192.0*volume), 4.0*M_PI*mumurd2e/volume};
   double sum[EWALD_NFUNCS];
   int func[EWALD_NFUNCS];
@@ -625,8 +628,9 @@ void EwaldN::compute_virial(int vflag)
 
   complex *cek = cek_global;
   double *kv = kvirial;
+  const double qscale = force->qqrd2e * scale;
   double c[EWALD_NFUNCS] = {
-    4.0*M_PI*qqrd2e*scale/volume, 2.0*M_PI*sqrt(M_PI)/(24.0*volume),
+    4.0*M_PI*qscale/volume, 2.0*M_PI*sqrt(M_PI)/(24.0*volume),
     2.0*M_PI*sqrt(M_PI)/(192.0*volume), 4.0*M_PI*mumurd2e/volume};
   shape sum[EWALD_NFUNCS];
   int func[EWALD_NFUNCS];
@@ -685,14 +689,16 @@ void EwaldN::compute_slabcorr(int eflag)
   
   while ((x+=3)<xn) dipole += *x * *(q++);
   MPI_Allreduce(&dipole, &dipole_all, 1, MPI_DOUBLE, MPI_SUM, world);
+
+  const double qscale = force->qqrd2e * scale;
   
-  double ffact = -4.0*M_PI*qqrd2e*scale*dipole_all/volume;
+  double ffact = -4.0*M_PI*qscale*dipole_all/volume;
   double *f = atom->f[0]-1, *fn = f+3*atom->nlocal-3;
   
   q = atom->q;
   while ((f+=3)<fn) *f += ffact* *(q++);
 
   if (eflag) 						// energy correction
-    energy += qqrd2e*scale*(2.0*M_PI*dipole_all*dipole_all/volume);
+    energy += qscale*(2.0*M_PI*dipole_all*dipole_all/volume);
 }
 
diff --git a/src/USER-EWALDN/ewald_n.h b/src/USER-EWALDN/ewald_n.h
index 721b45ccee..2bec485323 100644
--- a/src/USER-EWALDN/ewald_n.h
+++ b/src/USER-EWALDN/ewald_n.h
@@ -57,7 +57,7 @@ class EwaldN : public KSpace {
   hvector *hvec;
   kvector *kvec;
 
-  double qqrd2e, mumurd2e, dielectric, *B, volume;
+  double mumurd2e, dielectric, *B, volume;
   struct Sum { double x, x2; } sum[EWALD_MAX_NSUMS];
   complex *cek_local, *cek_global;
  
diff --git a/src/USER-OMP/ewald_omp.cpp b/src/USER-OMP/ewald_omp.cpp
index 613b5a5e4f..4daf255c22 100644
--- a/src/USER-OMP/ewald_omp.cpp
+++ b/src/USER-OMP/ewald_omp.cpp
@@ -19,6 +19,7 @@
 #include "ewald_omp.h"
 #include "atom.h"
 #include "comm.h"
+#include "force.h"
 #include "memory.h"
 
 #include <math.h>
@@ -125,9 +126,10 @@ void EwaldOMP::compute(int eflag, int vflag)
     }
 
     // convert E-field to force
+    const double qscale = force->qqrd2e * scale;
 
     for (i = ifrom; i < ito; i++) {
-      const double fac = qqrd2e*scale*q[i];
+      const double fac = qscale*q[i];
       f[i][0] += fac*ek[i][0];
       f[i][1] += fac*ek[i][1];
       f[i][2] += fac*ek[i][2];
@@ -146,7 +148,7 @@ void EwaldOMP::compute(int eflag, int vflag)
 
 	eng_tmp -= g_ewald*qsqsum/MY_PIS +
 	  MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);
-	eng_tmp *= qqrd2e*scale;
+	eng_tmp *= qscale;
 	energy = eng_tmp;
       }
 
@@ -160,7 +162,7 @@ void EwaldOMP::compute(int eflag, int vflag)
 	  for (i = 0; i < 6; i++) v[i] += uk*vg[k][i];
 	}
 
-	for (i = 0; i < 6; i++) virial[i] = v[i] * qqrd2e*scale;
+	for (i = 0; i < 6; i++) virial[i] = v[i] * qscale;
       }
     }
 
@@ -180,7 +182,7 @@ void EwaldOMP::eik_dot_r()
   const int nthreads = comm->nthreads;
   
 #if defined(_OPENMP)
-#pragma omp parallel
+#pragma omp parallel default(none)
 #endif
   {
     int i,ifrom,ito,k,l,m,n,ic,tid;
diff --git a/src/USER-OMP/pppm_omp.cpp b/src/USER-OMP/pppm_omp.cpp
index bf4613feda..58d67fbad2 100644
--- a/src/USER-OMP/pppm_omp.cpp
+++ b/src/USER-OMP/pppm_omp.cpp
@@ -18,6 +18,7 @@
 #include "pppm_omp.h"
 #include "atom.h"
 #include "comm.h"
+#include "force.h"
 #include "memory.h"
 
 #include <string.h>
@@ -112,7 +113,7 @@ static void data_reduce_fft(FFT_SCALAR *dall, int nall, int nthreads, int ndim,
 void PPPMOMP::deallocate()
 {
   PPPM::deallocate();
-  for (int i; i < comm->nthreads; ++i) {
+  for (int i=0; i < comm->nthreads; ++i) {
     ThrData * thr = fix->get_thr(i);
     double ** rho1d_thr = static_cast<double **>(thr->get_rho1d());
     memory->destroy2d_offset(rho1d_thr,-order/2);
@@ -153,7 +154,7 @@ void PPPMOMP::make_rho()
   const int nthreads = comm->nthreads;
 
 #if defined(_OPENMP)
-#pragma omp parallel default(none)
+#pragma omp parallel default(none) shared(atom,fix,density_brick,part2grid,boxlo,delxinv,delyinv,delzinv,delvolinv,ngrid,nlower,nupper,shiftone,nxlo_out,nylo_out,nzhi_out,nzlo_out)
   {  
     // each thread works on a fixed chunk of atoms.
     const int tid = omp_get_thread_num();
@@ -233,6 +234,8 @@ void PPPMOMP::fieldforce()
 
   const double * const q = atom->q;
   const double * const * const x = atom->x;
+  const int nthreads = comm->nthreads;
+  const double qqrd2e = force->qqrd2e;
 
 #if defined(_OPENMP)
 #pragma omp parallel default(none)
@@ -240,7 +243,7 @@ void PPPMOMP::fieldforce()
     // each thread works on a fixed chunk of atoms.
     const int tid = omp_get_thread_num();
     const int inum = atom->nlocal;
-    const int idelta = 1 + inum/comm->nthreads;
+    const int idelta = 1 + inum/nthreads;
     const int ifrom = tid*idelta;
     const int ito = ((ifrom + idelta) > inum) ? inum : ifrom + idelta;
 #else