diff --git a/src/USER-CG-CMM/Package.sh b/src/USER-CG-CMM/Package.sh
index db51eed742..9b6101979d 100644
--- a/src/USER-CG-CMM/Package.sh
+++ b/src/USER-CG-CMM/Package.sh
@@ -3,10 +3,10 @@
 # do not copy molecular and kspace files if corresponding versions do not exist
 
 for file in *.cpp *.h; do
-  if (test $file = angle_cg_cmm.cpp -a ! -e ../pair_angle_harmonic.cpp) then
+  if (test $file = angle_cg_cmm.cpp -a ! -e ../angle_harmonic.cpp) then
     continue
   fi
-  if (test $file = angle_cg_cmm.h -a ! -e ../pair_angle_harmonic.h) then
+  if (test $file = angle_cg_cmm.h -a ! -e ../angle_harmonic.h) then
     continue
   fi
   if (test $file = pair_cg_cmm_coul_long.cpp -a ! -e ../pair_lj_cut_coul_long.cpp) then
diff --git a/src/USER-CG-CMM/pair_cg_cmm_coul_long.cpp b/src/USER-CG-CMM/pair_cg_cmm_coul_long.cpp
index fa7769927a..89eda915f1 100644
--- a/src/USER-CG-CMM/pair_cg_cmm_coul_long.cpp
+++ b/src/USER-CG-CMM/pair_cg_cmm_coul_long.cpp
@@ -444,7 +444,7 @@ double PairCGCMMCoulLong::single(int i, int j, int itype, int jtype, double rsq,
 
 /* ---------------------------------------------------------------------- */
 
-void *PairCGCMMCoulLong::extract(char *str, int &dim)
+void *PairCGCMMCoulLong::extract(const char *str, int &dim)
 {
   dim = 0;
   if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul_global;
diff --git a/src/USER-CG-CMM/pair_cg_cmm_coul_long.h b/src/USER-CG-CMM/pair_cg_cmm_coul_long.h
index 7bc15a853c..ce1fa044bc 100644
--- a/src/USER-CG-CMM/pair_cg_cmm_coul_long.h
+++ b/src/USER-CG-CMM/pair_cg_cmm_coul_long.h
@@ -1,4 +1,4 @@
-/* ----------------------------------------------------------------------
+/* -*- c++ -*- ----------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
@@ -43,7 +43,7 @@ namespace LAMMPS_NS {
     double memory_usage();
 
     double single(int, int, int, int, double, double, double, double &);
-    void *extract(char *str, int &);
+    void *extract(const char *str, int &);
 
     protected:
     void allocate();
diff --git a/src/USER-CG-CMM/pair_lj_sdk.cpp b/src/USER-CG-CMM/pair_lj_sdk.cpp
index a31b9eac48..9dac84d313 100644
--- a/src/USER-CG-CMM/pair_lj_sdk.cpp
+++ b/src/USER-CG-CMM/pair_lj_sdk.cpp
@@ -453,7 +453,7 @@ double PairLJSDK::single(int, int, int itype, int jtype, double rsq,
 
 /* ---------------------------------------------------------------------- */
 
-void *PairLJSDK::extract(char *str, int &dim)
+void *PairLJSDK::extract(const char *str, int &dim)
 {
   dim = 2;
   if (strcmp(str,"epsilon") == 0) return (void *) epsilon;
diff --git a/src/USER-CG-CMM/pair_lj_sdk.h b/src/USER-CG-CMM/pair_lj_sdk.h
index 65aea005c3..9953b4aede 100644
--- a/src/USER-CG-CMM/pair_lj_sdk.h
+++ b/src/USER-CG-CMM/pair_lj_sdk.h
@@ -43,7 +43,7 @@ class PairLJSDK : public Pair {
   void write_restart_settings(FILE *);
   void read_restart_settings(FILE *);
   double single(int, int, int, int, double, double, double, double &);
-  void *extract(char *, int &);
+  void *extract(const char *, int &);
   virtual double memory_usage();
 
  protected:
diff --git a/src/USER-CG-CMM/pair_lj_sdk_coul_long.cpp b/src/USER-CG-CMM/pair_lj_sdk_coul_long.cpp
index 24a0718b70..90e102e993 100644
--- a/src/USER-CG-CMM/pair_lj_sdk_coul_long.cpp
+++ b/src/USER-CG-CMM/pair_lj_sdk_coul_long.cpp
@@ -726,7 +726,7 @@ double PairLJSDKCoulLong::single(int i, int j, int itype, int jtype,
 
 /* ---------------------------------------------------------------------- */
 
-void *PairLJSDKCoulLong::extract(char *str, int &dim)
+void *PairLJSDKCoulLong::extract(const char *str, int &dim)
 {
   dim = 2;
   if (strcmp(str,"epsilon") == 0) return (void *) epsilon;
diff --git a/src/USER-CG-CMM/pair_lj_sdk_coul_long.h b/src/USER-CG-CMM/pair_lj_sdk_coul_long.h
index f4ed4b4d66..61043ad2b0 100644
--- a/src/USER-CG-CMM/pair_lj_sdk_coul_long.h
+++ b/src/USER-CG-CMM/pair_lj_sdk_coul_long.h
@@ -43,7 +43,7 @@ class PairLJSDKCoulLong : public Pair {
   virtual void write_restart_settings(FILE *);
   virtual void read_restart_settings(FILE *);
   virtual double single(int, int, int, int, double, double, double, double &);
-  void *extract(char *, int &);
+  void *extract(const char *, int &);
   virtual double memory_usage();
 
  protected:
diff --git a/src/USER-EWALDN/pair_buck_coul.cpp b/src/USER-EWALDN/pair_buck_coul.cpp
index 74e9e56476..5e3ebdb1f5 100644
--- a/src/USER-EWALDN/pair_buck_coul.cpp
+++ b/src/USER-EWALDN/pair_buck_coul.cpp
@@ -168,7 +168,7 @@ void PairBuckCoul::allocate()
    extract protected data from object
 ------------------------------------------------------------------------- */
 
-void *PairBuckCoul::extract(char *id, int &dim)
+void *PairBuckCoul::extract(const char *id, int &dim)
 {
   char *ids[] = {
     "B", "ewald_order", "ewald_cut", "ewald_mix", "cut_coul", NULL};
diff --git a/src/USER-EWALDN/pair_buck_coul.h b/src/USER-EWALDN/pair_buck_coul.h
index 9c2856f42f..5b9d406cbc 100644
--- a/src/USER-EWALDN/pair_buck_coul.h
+++ b/src/USER-EWALDN/pair_buck_coul.h
@@ -1,4 +1,4 @@
-/* ----------------------------------------------------------------------
+/* -*- c++ -*- ----------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
@@ -43,7 +43,7 @@ class PairBuckCoul : public Pair {
   void write_restart_settings(FILE *);
   void read_restart_settings(FILE *);
   double single(int, int, int, int, double, double, double, double &);
-  void *extract(char *, int &);
+  void *extract(const char *, int &);
 
   void compute_inner();
   void compute_middle();
diff --git a/src/USER-EWALDN/pair_lj_coul.cpp b/src/USER-EWALDN/pair_lj_coul.cpp
index 4e1692680f..b2bbd144b0 100644
--- a/src/USER-EWALDN/pair_lj_coul.cpp
+++ b/src/USER-EWALDN/pair_lj_coul.cpp
@@ -164,7 +164,7 @@ void PairLJCoul::allocate()
    extract protected data from object
 ------------------------------------------------------------------------- */
 
-void *PairLJCoul::extract(char *id, int &dim)
+void *PairLJCoul::extract(const char *id, int &dim)
 {
   char *ids[] = {
     "B", "sigma", "epsilon", "ewald_order", "ewald_cut", "ewald_mix",
diff --git a/src/USER-EWALDN/pair_lj_coul.h b/src/USER-EWALDN/pair_lj_coul.h
index cc74f35eee..4f17095770 100644
--- a/src/USER-EWALDN/pair_lj_coul.h
+++ b/src/USER-EWALDN/pair_lj_coul.h
@@ -1,4 +1,4 @@
-/* ----------------------------------------------------------------------
+/* -*- c++ -*- ----------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
@@ -42,7 +42,7 @@ class PairLJCoul : public Pair {
   void write_restart_settings(FILE *);
   void read_restart_settings(FILE *);
   double single(int, int, int, int, double, double, double, double &);
-  void *extract(char *, int &);
+  void *extract(const char *, int &);
 
   void compute_inner();
   void compute_middle();
diff --git a/src/USER-OMP/pair_born_coul_long_omp.cpp b/src/USER-OMP/pair_born_coul_long_omp.cpp
index cf409f3cfc..5a51bdd878 100644
--- a/src/USER-OMP/pair_born_coul_long_omp.cpp
+++ b/src/USER-OMP/pair_born_coul_long_omp.cpp
@@ -171,8 +171,8 @@ void PairBornCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
 	    evdwl = a[itype][jtype]*rexp - c[itype][jtype]*r6inv 
 	      + d[itype][jtype]*r6inv*r2inv - offset[itype][jtype];
 	    evdwl *= factor_lj;
-	  }
-	} else evdwl = 0.0;
+	  } else evdwl = 0.0;
+	}
 
 	if (EVFLAG) ev_tally_thr(this, i,j,nlocal,NEWTON_PAIR,
 				 evdwl,ecoul,fpair,delx,dely,delz,thr);
diff --git a/src/USER-OMP/pair_buck_coul_cut_omp.cpp b/src/USER-OMP/pair_buck_coul_cut_omp.cpp
index 235f1c4f2c..0369852f47 100644
--- a/src/USER-OMP/pair_buck_coul_cut_omp.cpp
+++ b/src/USER-OMP/pair_buck_coul_cut_omp.cpp
@@ -153,8 +153,8 @@ void PairBuckCoulCutOMP::eval(int iifrom, int iito, ThrData * const thr)
 	    evdwl = a[itype][jtype]*rexp - c[itype][jtype]*r6inv -
 	      offset[itype][jtype];
 	    evdwl *= factor_lj;
-	  }
-	} else evdwl = 0.0;
+	  } else evdwl = 0.0;
+	}
 
 	if (EVFLAG) ev_tally_thr(this, i,j,nlocal,NEWTON_PAIR,
 				 evdwl,ecoul,fpair,delx,dely,delz,thr);
diff --git a/src/USER-OMP/pair_buck_coul_long_omp.cpp b/src/USER-OMP/pair_buck_coul_long_omp.cpp
index 083b9acc6e..54a6f32803 100644
--- a/src/USER-OMP/pair_buck_coul_long_omp.cpp
+++ b/src/USER-OMP/pair_buck_coul_long_omp.cpp
@@ -171,8 +171,8 @@ void PairBuckCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
 	    evdwl = a[itype][jtype]*rexp - c[itype][jtype]*r6inv -
 	      offset[itype][jtype];
 	    evdwl *= factor_lj;
-	  }
-	} else evdwl = 0.0;
+	  } else evdwl = 0.0;
+	}
 
 	if (EVFLAG) ev_tally_thr(this, i,j,nlocal,NEWTON_PAIR,
 				 evdwl,ecoul,fpair,delx,dely,delz,thr);
diff --git a/src/USER-OMP/pair_dipole_cut_omp.cpp b/src/USER-OMP/pair_dipole_cut_omp.cpp
index 85079dd718..3ec3445396 100644
--- a/src/USER-OMP/pair_dipole_cut_omp.cpp
+++ b/src/USER-OMP/pair_dipole_cut_omp.cpp
@@ -107,6 +107,7 @@ void PairDipoleCutOMP::eval(int iifrom, int iito, ThrData * const thr)
     xtmp = x[i][0];
     ytmp = x[i][1];
     ztmp = x[i][2];
+    qtmp = q[i];
     itype = type[i];
     jlist = firstneigh[i];
     jnum = numneigh[i];
diff --git a/src/USER-OMP/pair_dipole_sf_omp.cpp b/src/USER-OMP/pair_dipole_sf_omp.cpp
index b920ff5c83..316a4e75ac 100644
--- a/src/USER-OMP/pair_dipole_sf_omp.cpp
+++ b/src/USER-OMP/pair_dipole_sf_omp.cpp
@@ -111,6 +111,7 @@ void PairDipoleSFOMP::eval(int iifrom, int iito, ThrData * const thr)
     xtmp = x[i][0];
     ytmp = x[i][1];
     ztmp = x[i][2];
+    qtmp = q[i];
     itype = type[i];
     jlist = firstneigh[i];
     jnum = numneigh[i];
diff --git a/src/USER-OMP/pair_lj_charmm_coul_pppm_omp.cpp b/src/USER-OMP/pair_lj_charmm_coul_pppm_omp.cpp
index 3b12668549..ecad175c18 100644
--- a/src/USER-OMP/pair_lj_charmm_coul_pppm_omp.cpp
+++ b/src/USER-OMP/pair_lj_charmm_coul_pppm_omp.cpp
@@ -85,7 +85,6 @@ void PairLJCharmmCoulPPPMOMP::compute(int eflag, int vflag)
 
     // thread id 0 runs pppm, the rest the pair style
     if (tid < nproxy) {
-      if (update->setupflag) kspace->setup_proxy();
       kspace->compute_proxy(eflag,vflag);
     } else {
       if (evflag) {
diff --git a/src/USER-OMP/pair_lj_cut_coul_pppm_omp.cpp b/src/USER-OMP/pair_lj_cut_coul_pppm_omp.cpp
index c048db9d4c..d72b9c0b0a 100644
--- a/src/USER-OMP/pair_lj_cut_coul_pppm_omp.cpp
+++ b/src/USER-OMP/pair_lj_cut_coul_pppm_omp.cpp
@@ -85,7 +85,6 @@ void PairLJCutCoulPPPMOMP::compute(int eflag, int vflag)
 
     // thread id 0 runs pppm, the rest the pair style
     if (tid < nproxy) {
-      if (update->setupflag) kspace->setup_proxy();
       kspace->compute_proxy(eflag,vflag);
     } else {
       if (evflag) {
diff --git a/src/USER-OMP/pair_lj_cut_coul_pppm_tip4p_omp.cpp b/src/USER-OMP/pair_lj_cut_coul_pppm_tip4p_omp.cpp
index add0e295e2..a4aa811075 100644
--- a/src/USER-OMP/pair_lj_cut_coul_pppm_tip4p_omp.cpp
+++ b/src/USER-OMP/pair_lj_cut_coul_pppm_tip4p_omp.cpp
@@ -133,7 +133,6 @@ void PairLJCutCoulPPPMTIP4POMP::compute(int eflag, int vflag)
 
     // thread id 0 runs pppm, the rest the pair style
     if (tid < nproxy) {
-      if (update->setupflag) kspace->setup_proxy();
       kspace->compute_proxy(eflag,vflag);
     } else {
       if (evflag) {
diff --git a/src/USER-OMP/pppm_omp.cpp b/src/USER-OMP/pppm_omp.cpp
index 50c0b5422b..70a66f27f2 100644
--- a/src/USER-OMP/pppm_omp.cpp
+++ b/src/USER-OMP/pppm_omp.cpp
@@ -18,12 +18,16 @@
 #include "pppm_omp.h"
 #include "atom.h"
 #include "comm.h"
+#include "domain.h"
 #include "force.h"
 #include "memory.h"
+#include "math_const.h"
 
 #include <string.h>
+#include <math.h>
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #ifdef FFT_SINGLE
 #define ZEROF 0.0f
@@ -33,6 +37,8 @@ using namespace LAMMPS_NS;
 #define ONEF  1.0
 #endif
 
+#define EPS_HOC 1.0e-7
+
 /* ---------------------------------------------------------------------- */
 
 PPPMOMP::PPPMOMP(LAMMPS *lmp, int narg, char **arg) :
@@ -40,7 +46,6 @@ PPPMOMP::PPPMOMP(LAMMPS *lmp, int narg, char **arg) :
 {
 }
 
-
 /* ----------------------------------------------------------------------
    allocate memory that depends on # of K-vectors and order 
 ------------------------------------------------------------------------- */
@@ -82,7 +87,7 @@ void PPPMOMP::allocate()
 // multi-dimensional arrays like force stored in this order
 // x1,y1,z1,x2,y2,z2,...
 // we need to post a barrier to wait until all threads are done
-// with writing to the array .
+// writing to the array.
 static void data_reduce_fft(FFT_SCALAR *dall, int nall, int nthreads, int ndim, int tid)
 {
 #if defined(_OPENMP)
@@ -95,10 +100,13 @@ static void data_reduce_fft(FFT_SCALAR *dall, int nall, int nthreads, int ndim,
     const int ifrom = tid*idelta;
     const int ito   = ((ifrom + idelta) > nvals) ? nvals : (ifrom + idelta);
 
-    for (int m = ifrom; m < ito; ++m) {
-      for (int n = 1; n < nthreads; ++n) {
-	dall[m] += dall[n*nvals + m];
-	dall[n*nvals + m] = 0.0;
+    // this if protects against having more threads than atoms
+    if (ifrom < nall) { 
+      for (int m = ifrom; m < ito; ++m) {
+	for (int n = 1; n < nthreads; ++n) {
+	  dall[m] += dall[n*nvals + m];
+	  dall[n*nvals + m] = 0.0;
+	}
       }
     }
   }
@@ -122,6 +130,185 @@ void PPPMOMP::deallocate()
   }
 }
 
+/* ----------------------------------------------------------------------
+   adjust PPPM coeffs, called initially and whenever volume has changed 
+------------------------------------------------------------------------- */
+
+void PPPMOMP::setup()
+{
+  int i,j,k,n;
+  double *prd;
+
+  // volume-dependent factors
+  // adjust z dimension for 2d slab PPPM
+  // z dimension for 3d PPPM is zprd since slab_volfactor = 1.0
+
+  if (triclinic == 0) prd = domain->prd;
+  else prd = domain->prd_lamda;
+
+  const double xprd = prd[0];
+  const double yprd = prd[1];
+  const double zprd = prd[2];
+  const double zprd_slab = zprd*slab_volfactor;
+  volume = xprd * yprd * zprd_slab;
+    
+  delxinv = nx_pppm/xprd;
+  delyinv = ny_pppm/yprd;
+  delzinv = nz_pppm/zprd_slab;
+
+  delvolinv = delxinv*delyinv*delzinv;
+
+  const double unitkx = (2.0*MY_PI/xprd);
+  const double unitky = (2.0*MY_PI/yprd);
+  const double unitkz = (2.0*MY_PI/zprd_slab);
+
+  // fkx,fky,fkz for my FFT grid pts
+
+  double per;
+
+  for (i = nxlo_fft; i <= nxhi_fft; i++) {
+    per = i - nx_pppm*(2*i/nx_pppm);
+    fkx[i] = unitkx*per;
+  }
+
+  for (i = nylo_fft; i <= nyhi_fft; i++) {
+    per = i - ny_pppm*(2*i/ny_pppm);
+    fky[i] = unitky*per;
+  }
+
+  for (i = nzlo_fft; i <= nzhi_fft; i++) {
+    per = i - nz_pppm*(2*i/nz_pppm);
+    fkz[i] = unitkz*per;
+  }
+
+  // virial coefficients
+
+  double sqk,vterm;
+
+  n = 0;
+  for (k = nzlo_fft; k <= nzhi_fft; k++) {
+    for (j = nylo_fft; j <= nyhi_fft; j++) {
+      for (i = nxlo_fft; i <= nxhi_fft; i++) {
+	sqk = fkx[i]*fkx[i] + fky[j]*fky[j] + fkz[k]*fkz[k];
+	if (sqk == 0.0) {
+	  vg[n][0] = 0.0;
+	  vg[n][1] = 0.0;
+	  vg[n][2] = 0.0;
+	  vg[n][3] = 0.0;
+	  vg[n][4] = 0.0;
+	  vg[n][5] = 0.0;
+	} else {
+	  vterm = -2.0 * (1.0/sqk + 0.25/(g_ewald*g_ewald));
+	  vg[n][0] = 1.0 + vterm*fkx[i]*fkx[i];
+	  vg[n][1] = 1.0 + vterm*fky[j]*fky[j];
+	  vg[n][2] = 1.0 + vterm*fkz[k]*fkz[k];
+	  vg[n][3] = vterm*fkx[i]*fky[j];
+	  vg[n][4] = vterm*fkx[i]*fkz[k];
+	  vg[n][5] = vterm*fky[j]*fkz[k];
+	}
+	n++;
+      }
+    }
+  }
+
+  // modified (Hockney-Eastwood) Coulomb Green's function
+
+  const int nbx = static_cast<int> ((g_ewald*xprd/(MY_PI*nx_pppm)) * 
+				    pow(-log(EPS_HOC),0.25));
+  const int nby = static_cast<int> ((g_ewald*yprd/(MY_PI*ny_pppm)) * 
+				    pow(-log(EPS_HOC),0.25));
+  const int nbz = static_cast<int> ((g_ewald*zprd_slab/(MY_PI*nz_pppm)) * 
+				    pow(-log(EPS_HOC),0.25));
+
+  const double form = 1.0;
+
+#if defined(_OPENMP)
+#pragma omp parallel default(none)
+#endif
+  {
+    int tid,nn,nnfrom,nnto,nx,ny,nz,k,l,m;
+    double snx,sny,snz,sqk;
+    double argx,argy,argz,wx,wy,wz,sx,sy,sz,qx,qy,qz;
+    double sum1,dot1,dot2;
+    double numerator,denominator;
+    const double gew2 = -4.0*g_ewald*g_ewald;
+    
+    const int nnx = nxhi_fft-nxlo_fft+1;
+    const int nny = nyhi_fft-nylo_fft+1;
+    
+    loop_setup_thr(nnfrom, nnto, tid, nfft, comm->nthreads);
+    
+    for (m = nzlo_fft; m <= nzhi_fft; m++) {
+      
+      const double fkzm = fkz[m];
+      snz = sin(0.5*fkzm*zprd_slab/nz_pppm);
+      snz *= snz;
+
+      for (l = nylo_fft; l <= nyhi_fft; l++) {
+	const double fkyl = fky[l];
+	sny = sin(0.5*fkyl*yprd/ny_pppm);
+	sny *= sny;
+
+	for (k = nxlo_fft; k <= nxhi_fft; k++) {
+
+	  /* only compute the part designated to this thread */
+	  nn = k-nxlo_fft + nnx*(l-nylo_fft + nny*(m-nzlo_fft));
+	  if ((nn < nnfrom) || (nn >=nnto)) continue;
+
+	  const double fkxk = fkx[k];
+	  snx = sin(0.5*fkxk*xprd/nx_pppm);
+	  snx *= snx;
+      
+	  sqk = fkxk*fkxk + fkyl*fkyl + fkzm*fkzm;
+
+	  if (sqk != 0.0) {
+	    numerator = form*MY_4PI/sqk;
+	    denominator = gf_denom(snx,sny,snz);  
+	    sum1 = 0.0;
+	    for (nx = -nbx; nx <= nbx; nx++) {
+	      qx = fkxk + unitkx*nx_pppm*nx;
+	      sx = exp(qx*qx/gew2);
+	      wx = 1.0;
+	      argx = 0.5*qx*xprd/nx_pppm;
+	      if (argx != 0.0) wx = pow(sin(argx)/argx,order);
+	      wx *=wx;
+
+	      for (ny = -nby; ny <= nby; ny++) {
+		qy = fkyl + unitky*ny_pppm*ny;
+		sy = exp(qy*qy/gew2);
+		wy = 1.0;
+		argy = 0.5*qy*yprd/ny_pppm;
+		if (argy != 0.0) wy = pow(sin(argy)/argy,order);
+		wy *= wy;
+
+		for (nz = -nbz; nz <= nbz; nz++) {
+		  qz = fkzm + unitkz*nz_pppm*nz;
+		  sz = exp(qz*qz/gew2);
+		  wz = 1.0;
+		  argz = 0.5*qz*zprd_slab/nz_pppm;
+		  if (argz != 0.0) wz = pow(sin(argz)/argz,order);
+		  wz *= wz;
+
+		  dot1 = fkxk*qx + fkyl*qy + fkzm*qz;
+		  dot2 = qx*qx+qy*qy+qz*qz;
+		  sum1 += (dot1/dot2) * sx*sy*sz * wx*wy*wz;
+		}
+	      }
+	    }
+	    greensfn[nn] = numerator*sum1/denominator;
+	  } else greensfn[nn] = 0.0;
+	}
+      }
+    }
+  }
+}
+
+/* ----------------------------------------------------------------------
+   run the regular toplevel compute method from plain PPPPM 
+   which will have individual methods replaced by our threaded
+   versions and then call the obligatory force reduction.
+------------------------------------------------------------------------- */
+
 void PPPMOMP::compute(int eflag, int vflag)
 {
 
@@ -153,20 +340,21 @@ void PPPMOMP::make_rho()
   const double * const q = atom->q;
   const double * const * const x = atom->x;
   const int nthreads = comm->nthreads;
+  const int nlocal = atom->nlocal;
 
 #if defined(_OPENMP)
 #pragma omp parallel default(none)
   {  
     // each thread works on a fixed chunk of atoms.
     const int tid = omp_get_thread_num();
-    const int inum = atom->nlocal;
+    const int inum = nlocal;
     const int idelta = 1 + inum/nthreads;
     const int ifrom = tid*idelta;
     const int ito = ((ifrom + idelta) > inum) ? inum : ifrom + idelta;
 #else
     const int tid = 0;
     const int ifrom = 0;
-    const int ito = atom->nlocal;
+    const int ito = nlocal;
 #endif
 
     int l,m,n,nx,ny,nz,mx,my,mz;
@@ -184,28 +372,31 @@ void PPPMOMP::make_rho()
     // (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
+    
+    // this if protects against having more threads than local atoms
+    if (ifrom < nlocal) { 
+      for (int i = ifrom; i < ito; i++) {
 
-    for (int i = ifrom; i < ito; 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;
 
-      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_thr(r1d,dx,dy,dz);
 
-      compute_rho1d_thr(r1d,dx,dy,dz);
-
-      z0 = delvolinv * q[i];
-      for (n = nlower; n <= nupper; n++) {
-	mz = n+nz;
-	y0 = z0*r1d[2][n];
-	for (m = nlower; m <= nupper; m++) {
-	  my = m+ny;
-	  x0 = y0*r1d[1][m];
-	  for (l = nlower; l <= nupper; l++) {
-	    mx = l+nx;
-	    db[mz][my][mx] += x0*r1d[0][l];
+	z0 = delvolinv * q[i];
+	for (n = nlower; n <= nupper; n++) {
+	  mz = n+nz;
+	  y0 = z0*r1d[2][n];
+	  for (m = nlower; m <= nupper; m++) {
+	    my = m+ny;
+	    x0 = y0*r1d[1][m];
+	    for (l = nlower; l <= nupper; l++) {
+	      mx = l+nx;
+	      db[mz][my][mx] += x0*r1d[0][l];
+	    }
 	  }
 	}
       }
@@ -213,7 +404,6 @@ void PPPMOMP::make_rho()
 #if defined(_OPENMP)
     // reduce 3d density array
     if (nthreads > 1) {
-      sync_threads();
       data_reduce_fft(&(density_brick[nzlo_out][nylo_out][nxlo_out]),ngrid,nthreads,1,tid);
     }
   }
@@ -235,6 +425,7 @@ void PPPMOMP::fieldforce()
   const double * const q = atom->q;
   const double * const * const x = atom->x;
   const int nthreads = comm->nthreads;
+  const int nlocal = atom->nlocal;
   const double qqrd2e = force->qqrd2e;
 
 #if defined(_OPENMP)
@@ -242,13 +433,13 @@ 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 inum = nlocal;
     const int idelta = 1 + inum/nthreads;
     const int ifrom = tid*idelta;
     const int ito = ((ifrom + idelta) > inum) ? inum : ifrom + idelta;
 #else
     const int ifrom = 0;
-    const int ito = atom->nlocal;
+    const int ito = nlocal;
     const int tid = 0;
 #endif
     ThrData *thr = fix->get_thr(tid);
@@ -259,39 +450,42 @@ void PPPMOMP::fieldforce()
     FFT_SCALAR dx,dy,dz,x0,y0,z0;
     FFT_SCALAR ekx,eky,ekz;
 
-    for (int i = ifrom; i < ito; i++) {
+    // this if protects against having more threads than local atoms
+    if (ifrom < nlocal) { 
+      for (int i = ifrom; i < ito; 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;
+	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_thr(r1d,dx,dy,dz);
+	compute_rho1d_thr(r1d,dx,dy,dz);
 
-      ekx = eky = ekz = ZEROF;
-      for (n = nlower; n <= nupper; n++) {
-	mz = n+nz;
-	z0 = r1d[2][n];
-	for (m = nlower; m <= nupper; m++) {
-	  my = m+ny;
-	  y0 = z0*r1d[1][m];
-	  for (l = nlower; l <= nupper; l++) {
-	    mx = l+nx;
-	    x0 = y0*r1d[0][l];
-	    ekx -= x0*vdx_brick[mz][my][mx];
-	    eky -= x0*vdy_brick[mz][my][mx];
-	    ekz -= x0*vdz_brick[mz][my][mx];
+	ekx = eky = ekz = ZEROF;
+	for (n = nlower; n <= nupper; n++) {
+	  mz = n+nz;
+	  z0 = r1d[2][n];
+	  for (m = nlower; m <= nupper; m++) {
+	    my = m+ny;
+	    y0 = z0*r1d[1][m];
+	    for (l = nlower; l <= nupper; l++) {
+	      mx = l+nx;
+	      x0 = y0*r1d[0][l];
+	      ekx -= x0*vdx_brick[mz][my][mx];
+	      eky -= x0*vdy_brick[mz][my][mx];
+	      ekz -= x0*vdz_brick[mz][my][mx];
+	    }
 	  }
 	}
-      }
 
-      // convert E-field to force
-      const double qfactor = qqrd2e*scale*q[i];
-      f[i][0] += qfactor*ekx;
-      f[i][1] += qfactor*eky;
-      f[i][2] += qfactor*ekz;
+	// convert E-field to force
+	const double qfactor = qqrd2e*scale*q[i];
+	f[i][0] += qfactor*ekx;
+	f[i][1] += qfactor*eky;
+	f[i][2] += qfactor*ekz;
+      }
     }
 #if defined(_OPENMP)
   }
diff --git a/src/USER-OMP/pppm_omp.h b/src/USER-OMP/pppm_omp.h
index 0701cb3396..8abcd64aa6 100644
--- a/src/USER-OMP/pppm_omp.h
+++ b/src/USER-OMP/pppm_omp.h
@@ -1,4 +1,4 @@
-/* ----------------------------------------------------------------------
+/* -*- c++ -*- ----------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
@@ -29,13 +29,15 @@ namespace LAMMPS_NS {
  public:
   PPPMOMP(class LAMMPS *, int, char **);
   virtual ~PPPMOMP () {};
+  virtual void setup();
+  virtual void compute(int, int);
 
  protected:
   virtual void allocate();
   virtual void deallocate();
-  virtual void compute(int, int);
   virtual void fieldforce();
   virtual void make_rho();
+
   void compute_rho1d_thr(FFT_SCALAR * const * const, const FFT_SCALAR &,
 			 const FFT_SCALAR &, const FFT_SCALAR &);
 //  void compute_rho_coeff();
diff --git a/src/USER-OMP/pppm_proxy.cpp b/src/USER-OMP/pppm_proxy.cpp
index 1b8a9c5453..4e8112060f 100644
--- a/src/USER-OMP/pppm_proxy.cpp
+++ b/src/USER-OMP/pppm_proxy.cpp
@@ -24,13 +24,16 @@ using namespace LAMMPS_NS;
 /* ---------------------------------------------------------------------- */
 
 PPPMProxy::PPPMProxy(LAMMPS *lmp, int narg, char **arg) :
-  PPPM(lmp, narg, arg), ThrOMP(lmp, THR_KSPACE|THR_PROXY) {}
+  PPPM(lmp, narg, arg), ThrOMP(lmp, THR_KSPACE|THR_PROXY) { need_setup=1; }
 
 /* ---------------------------------------------------------------------- */
 
 void PPPMProxy::setup_proxy()
 {
-  PPPM::setup();
+  if (need_setup)
+    PPPM::setup();
+
+  need_setup = 0;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -56,6 +59,7 @@ void PPPMProxy::compute(int eflag, int vflag)
 
 void PPPMProxy::compute_proxy(int eflag, int vflag)
 {
+  setup_proxy();
   PPPM::compute(eflag,vflag);
 }
 
diff --git a/src/USER-OMP/pppm_proxy.h b/src/USER-OMP/pppm_proxy.h
index e7387569d8..505efe2e81 100644
--- a/src/USER-OMP/pppm_proxy.h
+++ b/src/USER-OMP/pppm_proxy.h
@@ -1,4 +1,4 @@
-/* ----------------------------------------------------------------------
+/* -*- c++ -*- ----------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
@@ -33,9 +33,12 @@ namespace LAMMPS_NS {
   virtual void compute(int, int);
   virtual void compute_proxy(int eflag, int vflag);
 
-  // nothing to do in setup
-  virtual void setup() {};
+  // setup is delayed until the compute proxy is called
+  virtual void setup() { need_setup=1; };
   virtual void setup_proxy();
+
+ protected:
+  int need_setup;
 };
 
 }
diff --git a/src/USER-OMP/pppm_tip4p_proxy.cpp b/src/USER-OMP/pppm_tip4p_proxy.cpp
index fa8d382c27..367b51e48d 100644
--- a/src/USER-OMP/pppm_tip4p_proxy.cpp
+++ b/src/USER-OMP/pppm_tip4p_proxy.cpp
@@ -24,13 +24,16 @@ using namespace LAMMPS_NS;
 /* ---------------------------------------------------------------------- */
 
 PPPMTIP4PProxy::PPPMTIP4PProxy(LAMMPS *lmp, int narg, char **arg) :
-  PPPMTIP4P(lmp, narg, arg), ThrOMP(lmp, THR_KSPACE|THR_PROXY) {}
+  PPPMTIP4P(lmp, narg, arg), ThrOMP(lmp, THR_KSPACE|THR_PROXY) { need_setup=1; }
 
 /* ---------------------------------------------------------------------- */
 
 void PPPMTIP4PProxy::setup_proxy()
 {
-  PPPMTIP4P::setup();
+  if (need_setup)
+    PPPMTIP4P::setup();
+
+  need_setup=0;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -56,6 +59,7 @@ void PPPMTIP4PProxy::compute(int eflag, int vflag)
 
 void PPPMTIP4PProxy::compute_proxy(int eflag, int vflag)
 {
+  setup_proxy();
   PPPMTIP4P::compute(eflag,vflag);
 }
 
diff --git a/src/USER-OMP/pppm_tip4p_proxy.h b/src/USER-OMP/pppm_tip4p_proxy.h
index b69cab2ad6..fee06800ff 100644
--- a/src/USER-OMP/pppm_tip4p_proxy.h
+++ b/src/USER-OMP/pppm_tip4p_proxy.h
@@ -1,4 +1,4 @@
-/* ----------------------------------------------------------------------
+/* -*- c++ -*- ----------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
@@ -33,9 +33,12 @@ namespace LAMMPS_NS {
   virtual void compute(int, int);
   virtual void compute_proxy(int eflag, int vflag);
 
-  // nothing to do in setup
-  virtual void setup() {};
+  // setup is delayed until the compute_proxy is called
+  virtual void setup() { need_setup=1; };
   virtual void setup_proxy();
+
+ protected:
+  int need_setup;
 };
 
 }
diff --git a/src/USER-OMP/thr_data.cpp b/src/USER-OMP/thr_data.cpp
index 069a5fc756..01fed3bc5f 100644
--- a/src/USER-OMP/thr_data.cpp
+++ b/src/USER-OMP/thr_data.cpp
@@ -215,10 +215,13 @@ void LAMMPS_NS::data_reduce_thr(double *dall, int nall, int nthreads, int ndim,
     const int ifrom = tid*idelta;
     const int ito   = ((ifrom + idelta) > nvals) ? nvals : (ifrom + idelta);
 
-    for (int m = ifrom; m < ito; ++m) {
-      for (int n = 1; n < nthreads; ++n) {
-	dall[m] += dall[n*nvals + m];
-	dall[n*nvals + m] = 0.0;
+    // this if protects against having more threads than atoms
+    if (ifrom < nvals) { 
+      for (int m = ifrom; m < ito; ++m) {
+	for (int n = 1; n < nthreads; ++n) {
+	  dall[m] += dall[n*nvals + m];
+	  dall[n*nvals + m] = 0.0;
+	}
       }
     }
   }
diff --git a/src/USER-REAXC/pair_reax_c.cpp b/src/USER-REAXC/pair_reax_c.cpp
index 0c30b765d9..0d54910465 100644
--- a/src/USER-REAXC/pair_reax_c.cpp
+++ b/src/USER-REAXC/pair_reax_c.cpp
@@ -745,7 +745,7 @@ void PairReaxC::read_reax_forces()
 
 /* ---------------------------------------------------------------------- */
 
-void *PairReaxC::extract(char *str, int &dim)
+void *PairReaxC::extract(const char *str, int &dim)
 {
   dim = 1;
   if (strcmp(str,"chi") == 0 && chi) {
diff --git a/src/fix_deposit.cpp b/src/fix_deposit.cpp
index 5f975dbee2..8bc15c34ba 100644
--- a/src/fix_deposit.cpp
+++ b/src/fix_deposit.cpp
@@ -61,6 +61,7 @@ FixDeposit::FixDeposit(LAMMPS *lmp, int narg, char **arg) :
   rateflag = 0;
   vxlo = vxhi = vylo = vyhi = vzlo = vzhi = 0.0;
   scaleflag = 1;
+  targetflag = 0;
 
   // read options from end of input line
 
@@ -125,6 +126,9 @@ FixDeposit::FixDeposit(LAMMPS *lmp, int narg, char **arg) :
   vyhi *= yscale;
   vzlo *= zscale;
   vzhi *= zscale;
+  tx *= xscale;
+  ty *= yscale;
+  tz *= zscale;
 
   // random number generator, same for all procs
 
@@ -284,6 +288,21 @@ void FixDeposit::pre_exchange()
     double vytmp = vylo + random->uniform() * (vyhi-vylo);
     double vztmp = vzlo + random->uniform() * (vzhi-vzlo);
 
+    // if we have a sputter target change velocity vector accordingly
+    if (targetflag) {
+      double vel = sqrt(vxtmp*vxtmp + vytmp*vytmp + vztmp*vztmp);
+      delx = tx - coord[0];
+      dely = ty - coord[1];
+      delz = tz - coord[2];
+      double rsq = delx*delx + dely*dely + delz*delz;
+      if (rsq > 0.0) {
+	double rinv = sqrt(1.0/rsq);
+	vxtmp = delx*rinv*vel;
+	vytmp = dely*rinv*vel;
+	vztmp = delz*rinv*vel;
+      }
+    }
+
     // check if new atom is in my sub-box or above it if I'm highest proc
     // if so, add to my list via create_atom()
     // initialize info about the atoms
@@ -419,6 +438,13 @@ void FixDeposit::options(int narg, char **arg)
       else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1;
       else error->all(FLERR,"Illegal fix deposit command");
       iarg += 2;
+    } else if (strcmp(arg[iarg],"target") == 0) {
+      if (iarg+4 > narg) error->all(FLERR,"Illegal fix deposit command");
+      tx = atof(arg[iarg+1]);
+      ty = atof(arg[iarg+2]);
+      tz = atof(arg[iarg+3]);
+      targetflag = 1;
+      iarg += 4;
     } else error->all(FLERR,"Illegal fix deposit command");
   }
 }
diff --git a/src/fix_deposit.h b/src/fix_deposit.h
index b844e54b3d..0b1f4b8028 100644
--- a/src/fix_deposit.h
+++ b/src/fix_deposit.h
@@ -1,4 +1,4 @@
-/* ----------------------------------------------------------------------
+/* -*- c++ -*- ----------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
@@ -37,11 +37,12 @@ class FixDeposit : public Fix {
 
  private:
   int ninsert,ntype,nfreq,seed;
-  int iregion,globalflag,localflag,maxattempt,rateflag,scaleflag;
+  int iregion,globalflag,localflag,maxattempt,rateflag,scaleflag,targetflag;
   char *idregion;
   double lo,hi,deltasq,nearsq,rate;
   double vxlo,vxhi,vylo,vyhi,vzlo,vzhi;
   double xlo,xhi,ylo,yhi,zlo,zhi;
+  double tx,ty,tz;
   int nfirst,ninserted;
   class RanPark *random;
 
diff --git a/src/fix_spring_self.cpp b/src/fix_spring_self.cpp
index 6fb3472fb6..4be64ed271 100644
--- a/src/fix_spring_self.cpp
+++ b/src/fix_spring_self.cpp
@@ -32,7 +32,8 @@ using namespace LAMMPS_NS;
 FixSpringSelf::FixSpringSelf(LAMMPS *lmp, int narg, char **arg) :
   Fix(lmp, narg, arg)
 {
-  if (narg != 4) error->all(FLERR,"Illegal fix spring/self command");
+  if ((narg < 4) || (narg > 5))
+    error->all(FLERR,"Illegal fix spring/self command");
 
   restart_peratom = 1;
   scalar_flag = 1;
@@ -42,6 +43,25 @@ FixSpringSelf::FixSpringSelf(LAMMPS *lmp, int narg, char **arg) :
   k = atof(arg[3]);
   if (k <= 0.0) error->all(FLERR,"Illegal fix spring/self command");
 
+  xflag = yflag = zflag = 1;
+  if (narg == 5) {
+    if (strcmp(arg[4],"xyz") == 0) {
+      ; /* default */
+    } else if (strcmp(arg[4],"xy") == 0) {
+      zflag = 0;
+    } else if (strcmp(arg[4],"xz") == 0) {
+      yflag = 0;
+    } else if (strcmp(arg[4],"yz") == 0) {
+      xflag = 0;
+    } else if (strcmp(arg[4],"x") == 0) {
+      yflag = zflag = 0;
+    } else if (strcmp(arg[4],"y") == 0) {
+      xflag = zflag = 0;
+    } else if (strcmp(arg[4],"z") == 0) {
+      xflag = yflag = 0;
+    } else error->all(FLERR,"Illegal fix spring/self command");
+  }
+
   // perform initial allocation of atom-based array
   // register with Atom class
 
@@ -155,6 +175,9 @@ void FixSpringSelf::post_force(int vflag)
       dx = x[i][0] + xbox*xprd - xoriginal[i][0];
       dy = x[i][1] + ybox*yprd - xoriginal[i][1];
       dz = x[i][2] + zbox*zprd - xoriginal[i][2];
+      if (!xflag) dx = 0.0;
+      if (!yflag) dy = 0.0;
+      if (!zflag) dz = 0.0;
       f[i][0] -= k*dx;
       f[i][1] -= k*dy;
       f[i][2] -= k*dz;
diff --git a/src/fix_spring_self.h b/src/fix_spring_self.h
index 5f65aa5263..e545b8fbf5 100644
--- a/src/fix_spring_self.h
+++ b/src/fix_spring_self.h
@@ -1,4 +1,4 @@
-/* ----------------------------------------------------------------------
+/* -*- c++ -*- ----------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
@@ -50,6 +50,7 @@ class FixSpringSelf : public Fix {
  private:
   double k,espring;
   double **xoriginal;         // original coords of atoms
+  int xflag, yflag, zflag;
   int nlevels_respa;
 };
 
diff --git a/src/lmptype.h b/src/lmptype.h
index a88e6d498a..f5739c3934 100644
--- a/src/lmptype.h
+++ b/src/lmptype.h
@@ -1,4 +1,4 @@
-/* ----------------------------------------------------------------------
+/* -*- c++ -*- ----------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
@@ -27,8 +27,13 @@
 #ifndef LMP_LMPTYPE_H
 #define LMP_LMPTYPE_H
 
+#ifndef __STDC_LIMIT_MACROS
 #define __STDC_LIMIT_MACROS
+#endif
+
+#ifndef __STDC_FORMAT_MACROS
 #define __STDC_FORMAT_MACROS
+#endif
 
 #include "limits.h"
 #include "stdint.h"
diff --git a/src/lmpwindows.h b/src/lmpwindows.h
index a3324f1f9d..a1cfd0c75e 100644
--- a/src/lmpwindows.h
+++ b/src/lmpwindows.h
@@ -1,9 +1,13 @@
 #include <iso646.h>
+#if !defined(__MINGW32_VERSION)
 #include "erf.h"
+#endif
 #include "direct.h"
 #include "math.h"
 // LAMMPS uses usleep with 100 ms arguments, no microsecond precision needed
+#if !defined(__MINGW32_VERSION)
 #include "sleep.h"
+#endif
 
 // some symbols have different names in Windows
 
@@ -20,6 +24,10 @@ inline double pow(int i, int j){
   return pow((double)i,(double) j);
 }
 
+inline double pow(double i, int j){
+  return pow(i,(double) j);
+}
+
 inline double sqrt(int i){
   return sqrt((double) i);
 }
@@ -36,6 +44,6 @@ inline double trunc(double x) {
 # define S_IRWXU 0
 # define S_IRGRP 0
 # define S_IXGRP 0
-inline int mkdir(const char *path, int flags){
+inline int mkdir(const char *path, int){
   return _mkdir(path);
 }
diff --git a/src/math_const.h b/src/math_const.h
index 5cbc2aee28..2d1cb3819f 100644
--- a/src/math_const.h
+++ b/src/math_const.h
@@ -1,4 +1,4 @@
-/* ----------------------------------------------------------------------
+/* -*- c++ -*- ----------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
@@ -21,6 +21,7 @@ namespace MathConst {
   static const double MY_PI  = 3.14159265358979323846; // pi
   static const double MY_2PI = 6.28318530717958647692; // 2pi
   static const double MY_3PI = 9.42477796076937971538; // 3pi
+  static const double MY_4PI = 12.56637061435917295384; // 4pi
   static const double MY_PI2 = 1.57079632679489661923; // pi/2
   static const double MY_PI4 = 0.78539816339744830962; // pi/4
   static const double MY_PIS = 1.77245385090551602729; // sqrt(pi)
diff --git a/src/memory.h b/src/memory.h
index 1e66207666..681af216ab 100644
--- a/src/memory.h
+++ b/src/memory.h
@@ -1,4 +1,4 @@
-/* ----------------------------------------------------------------------
+/* -*- c++ -*- ----------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
@@ -14,6 +14,7 @@
 #ifndef LMP_MEMORY_H
 #define LMP_MEMORY_H
 
+#include "lmptype.h"
 #include "pointers.h"
 
 namespace LAMMPS_NS {
diff --git a/src/pair_born_coul_wolf.cpp b/src/pair_born_coul_wolf.cpp
index 647c6bdaf6..0b33eb7549 100644
--- a/src/pair_born_coul_wolf.cpp
+++ b/src/pair_born_coul_wolf.cpp
@@ -69,7 +69,7 @@ void PairBornCoulWolf::compute(int eflag, int vflag)
   double prefactor;
   double r,rexp;
   int *ilist,*jlist,*numneigh,**firstneigh;
-  double erfcc,erfcd,v_sh,dvdrr,e_self,qisq;
+  double erfcc,erfcd,v_sh,dvdrr,e_self,e_shift,f_shift,qisq;
 
   evdwl = ecoul = 0.0;
   if (eflag || vflag) ev_setup(eflag,vflag);
@@ -116,13 +116,9 @@ void PairBornCoulWolf::compute(int eflag, int vflag)
 
     for (jj = 0; jj < jnum; jj++) {
       j = jlist[jj];
-
-      if (j < nall) factor_coul = factor_lj = 1.0;
-      else {
-	factor_coul = special_coul[j/nall];
-	factor_lj = special_lj[j/nall];
-	j %= nall;
-      }
+      factor_lj = special_lj[sbmask(j)];
+      factor_coul = special_coul[sbmask(j)];
+      j &= NEIGHMASK;
 
       delx = xtmp - x[j][0];
       dely = ytmp - x[j][1];
diff --git a/src/pair_born_coul_wolf.h b/src/pair_born_coul_wolf.h
index 56493c8d24..0792682a2a 100644
--- a/src/pair_born_coul_wolf.h
+++ b/src/pair_born_coul_wolf.h
@@ -1,4 +1,4 @@
-/* ----------------------------------------------------------------------
+/* -*- c++ -*- ----------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
@@ -27,8 +27,8 @@ namespace LAMMPS_NS {
 class PairBornCoulWolf : public Pair {
  public:
   PairBornCoulWolf(class LAMMPS *);
-  ~PairBornCoulWolf();
-  void compute(int, int);
+  virtual ~PairBornCoulWolf();
+  virtual void compute(int, int);
   void settings(int, char **);
   void coeff(int, char **);
   void init_style();
@@ -39,13 +39,12 @@ class PairBornCoulWolf : public Pair {
   void read_restart_settings(FILE *);
   double single(int, int, int, int, double, double, double, double &);
   
- private:
-  double cut_lj_global;
+ protected:
+  double cut_lj_global,alf;
   double **cut_lj,**cut_ljsq;
   double cut_coul,cut_coulsq;
   double **a,**rho,**sigma,**c,**d;
   double **rhoinv,**born1,**born2,**born3,**offset;
-  double alf,e_shift,f_shift;
 
   void allocate();
 };
diff --git a/src/pair_coul_wolf.cpp b/src/pair_coul_wolf.cpp
index d44a245426..a1382913c9 100644
--- a/src/pair_coul_wolf.cpp
+++ b/src/pair_coul_wolf.cpp
@@ -56,7 +56,7 @@ void PairCoulWolf::compute(int eflag, int vflag)
   double prefactor;
   double r,rexp;
   int *ilist,*jlist,*numneigh,**firstneigh;
-  double erfcc,erfcd,v_sh,dvdrr,e_self,qisq;
+  double erfcc,erfcd,v_sh,dvdrr,e_self,e_shift,f_shift,qisq;
 
   ecoul = 0.0;
   if (eflag || vflag) ev_setup(eflag,vflag);
@@ -102,12 +102,8 @@ void PairCoulWolf::compute(int eflag, int vflag)
 
     for (jj = 0; jj < jnum; jj++) {
       j = jlist[jj];
-
-      if (j < nall) factor_coul  = 1.0;
-      else {
-	factor_coul = special_coul[j/nall];
-	j %= nall;
-      }
+      factor_coul = special_coul[sbmask(j)];
+      j &= NEIGHMASK;
 
       delx = xtmp - x[j][0];
       dely = ytmp - x[j][1];
diff --git a/src/pair_coul_wolf.h b/src/pair_coul_wolf.h
index a7ed784b83..f4b6f37759 100644
--- a/src/pair_coul_wolf.h
+++ b/src/pair_coul_wolf.h
@@ -1,4 +1,4 @@
-/* ----------------------------------------------------------------------
+/* -*- c++ -*- ----------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
@@ -27,8 +27,8 @@ namespace LAMMPS_NS {
 class PairCoulWolf : public Pair {
  public:
   PairCoulWolf(class LAMMPS *);
-  ~PairCoulWolf();
-  void compute(int, int);
+  virtual ~PairCoulWolf();
+  virtual void compute(int, int);
   void settings(int, char **);
   void coeff(int, char **);
   void init_style();
@@ -39,9 +39,8 @@ class PairCoulWolf : public Pair {
   void read_restart_settings(FILE *);
   double single(int, int, int, int, double, double, double, double &);
 
- private:
-  double cut_coul,cut_coulsq;
-  double alf,e_shift,f_shift; 
+ protected:
+  double cut_coul,cut_coulsq,alf; 
 
   void allocate();
 };
diff --git a/src/region_cone.cpp b/src/region_cone.cpp
index bfa4e79716..adb409118e 100644
--- a/src/region_cone.cpp
+++ b/src/region_cone.cpp
@@ -317,7 +317,7 @@ int RegCone::surface_interior(double *x, double cutoff)
     r = sqrt(del1*del1 + del2*del2);
     currentradius = radiuslo + (x[2]-lo)*(radiushi-radiuslo)/(hi-lo);
 
-    // x is exterior to cone
+    // z is exterior to cone
     
     if (r > currentradius || x[2] < lo || x[2] > hi) return 0;