diff --git a/src/fix_wall.cpp b/src/fix_wall.cpp
index 83d34c306..503b87f4a 100644
--- a/src/fix_wall.cpp
+++ b/src/fix_wall.cpp
@@ -1,378 +1,378 @@
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
 
    Copyright (2003) Sandia Corporation.  Under the terms of Contract
    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
    certain rights in this software.  This software is distributed under
    the GNU General Public License.
 
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
 #include <math.h>
 #include <stdlib.h>
 #include <string.h>
 #include "fix_wall.h"
 #include "atom.h"
 #include "input.h"
 #include "variable.h"
 #include "domain.h"
 #include "lattice.h"
 #include "update.h"
 #include "modify.h"
 #include "respa.h"
 #include "error.h"
 #include "force.h"
 
 using namespace LAMMPS_NS;
 using namespace FixConst;
 
 enum{XLO=0,XHI=1,YLO=2,YHI=3,ZLO=4,ZHI=5};
 enum{NONE=0,EDGE,CONSTANT,VARIABLE};
 
 /* ---------------------------------------------------------------------- */
 
 FixWall::FixWall(LAMMPS *lmp, int narg, char **arg) :
-  Fix(lmp, narg, arg)
+  Fix(lmp, narg, arg),
+  nwall(0)
 {
   scalar_flag = 1;
   vector_flag = 1;
   global_freq = 1;
   extscalar = 1;
   extvector = 1;
   respa_level_support = 1;
   ilevel_respa = 0;
 
   // parse args
 
-  nwall = 0;
   int scaleflag = 1;
   fldflag = 0;
   int pbcflag = 0;
 
   for (int i = 0; i < 6; i++) xstr[i] = estr[i] = sstr[i] = NULL;
 
   int iarg = 3;
   while (iarg < narg) {
     if ((strcmp(arg[iarg],"xlo") == 0) || (strcmp(arg[iarg],"xhi") == 0) ||
         (strcmp(arg[iarg],"ylo") == 0) || (strcmp(arg[iarg],"yhi") == 0) ||
         (strcmp(arg[iarg],"zlo") == 0) || (strcmp(arg[iarg],"zhi") == 0)) {
       if (iarg+5 > narg) error->all(FLERR,"Illegal fix wall command");
 
       int newwall;
       if (strcmp(arg[iarg],"xlo") == 0) newwall = XLO;
       else if (strcmp(arg[iarg],"xhi") == 0) newwall = XHI;
       else if (strcmp(arg[iarg],"ylo") == 0) newwall = YLO;
       else if (strcmp(arg[iarg],"yhi") == 0) newwall = YHI;
       else if (strcmp(arg[iarg],"zlo") == 0) newwall = ZLO;
       else if (strcmp(arg[iarg],"zhi") == 0) newwall = ZHI;
 
       for (int m = 0; (m < nwall) && (m < 6); m++)
         if (newwall == wallwhich[m])
           error->all(FLERR,"Wall defined twice in fix wall command");
 
       wallwhich[nwall] = newwall;
       if (strcmp(arg[iarg+1],"EDGE") == 0) {
         xstyle[nwall] = EDGE;
         int dim = wallwhich[nwall] / 2;
         int side = wallwhich[nwall] % 2;
         if (side == 0) coord0[nwall] = domain->boxlo[dim];
         else coord0[nwall] = domain->boxhi[dim];
       } else if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
         xstyle[nwall] = VARIABLE;
         int n = strlen(&arg[iarg+1][2]) + 1;
         xstr[nwall] = new char[n];
         strcpy(xstr[nwall],&arg[iarg+1][2]);
       } else {
         xstyle[nwall] = CONSTANT;
         coord0[nwall] = force->numeric(FLERR,arg[iarg+1]);
       }
 
       if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) {
         int n = strlen(&arg[iarg+2][2]) + 1;
         estr[nwall] = new char[n];
         strcpy(estr[nwall],&arg[iarg+2][2]);
         estyle[nwall] = VARIABLE;
       } else {
         epsilon[nwall] = force->numeric(FLERR,arg[iarg+2]);
         estyle[nwall] = CONSTANT;
       }
 
       if (strstr(arg[iarg+3],"v_") == arg[iarg+3]) {
         int n = strlen(&arg[iarg+3][2]) + 1;
         sstr[nwall] = new char[n];
         strcpy(sstr[nwall],&arg[iarg+3][2]);
         sstyle[nwall] = VARIABLE;
       } else {
         sigma[nwall] = force->numeric(FLERR,arg[iarg+3]);
         sstyle[nwall] = CONSTANT;
       }
 
       cutoff[nwall] = force->numeric(FLERR,arg[iarg+4]);
       nwall++;
       iarg += 5;
 
     } else if (strcmp(arg[iarg],"units") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix wall command");
       if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0;
       else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1;
       else error->all(FLERR,"Illegal fix wall command");
       iarg += 2;
     } else if (strcmp(arg[iarg],"fld") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix wall command");
       if (strcmp(arg[iarg+1],"no") == 0) fldflag = 0;
       else if (strcmp(arg[iarg+1],"yes") == 0) fldflag = 1;
       else error->all(FLERR,"Illegal fix wall command");
       iarg += 2;
     } else if (strcmp(arg[iarg],"pbc") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix wall command");
       if (strcmp(arg[iarg+1],"yes") == 0) pbcflag = 1;
       else if (strcmp(arg[iarg+1],"no") == 0) pbcflag = 0;
       else error->all(FLERR,"Illegal fix wall command");
       iarg += 2;
     } else error->all(FLERR,"Illegal fix wall command");
   }
 
   size_vector = nwall;
 
   // error checks
 
   if (nwall == 0) error->all(FLERR,"Illegal fix wall command");
   for (int m = 0; m < nwall; m++)
     if (cutoff[m] <= 0.0)
       error->all(FLERR,"Fix wall cutoff <= 0.0");
 
   for (int m = 0; m < nwall; m++)
     if ((wallwhich[m] == ZLO || wallwhich[m] == ZHI) && domain->dimension == 2)
       error->all(FLERR,"Cannot use fix wall zlo/zhi for a 2d simulation");
 
   if (!pbcflag) {
     for (int m = 0; m < nwall; m++) {
       if ((wallwhich[m] == XLO || wallwhich[m] == XHI) && domain->xperiodic)
         error->all(FLERR,"Cannot use fix wall in periodic dimension");
       if ((wallwhich[m] == YLO || wallwhich[m] == YHI) && domain->yperiodic)
         error->all(FLERR,"Cannot use fix wall in periodic dimension");
       if ((wallwhich[m] == ZLO || wallwhich[m] == ZHI) && domain->zperiodic)
         error->all(FLERR,"Cannot use fix wall in periodic dimension");
     }
   }
 
   // scale factors for wall position for CONSTANT and VARIABLE walls
 
   int flag = 0;
   for (int m = 0; m < nwall; m++)
     if (xstyle[m] != EDGE) flag = 1;
 
   if (flag) {
     if (scaleflag) {
       xscale = domain->lattice->xlattice;
       yscale = domain->lattice->ylattice;
       zscale = domain->lattice->zlattice;
     }
     else xscale = yscale = zscale = 1.0;
 
     for (int m = 0; m < nwall; m++) {
       if (xstyle[m] != CONSTANT) continue;
       if (wallwhich[m] < YLO) coord0[m] *= xscale;
       else if (wallwhich[m] < ZLO) coord0[m] *= yscale;
       else coord0[m] *= zscale;
     }
   }
 
   // set xflag if any wall positions are variable
   // set varflag if any wall positions or parameters are variable
   // set wstyle to VARIABLE if either epsilon or sigma is a variable
 
   varflag = xflag = 0;
   for (int m = 0; m < nwall; m++) {
     if (xstyle[m] == VARIABLE) xflag = 1;
     if (xflag || estyle[m] == VARIABLE || sstyle[m] == VARIABLE) varflag = 1;
     if (estyle[m] == VARIABLE || sstyle[m] == VARIABLE) wstyle[m] = VARIABLE;
     else wstyle[m] = CONSTANT;
   }
 
   eflag = 0;
   for (int m = 0; m <= nwall; m++) ewall[m] = 0.0;
 }
 
 /* ---------------------------------------------------------------------- */
 
 FixWall::~FixWall()
 {
   for (int m = 0; m < nwall; m++) {
     delete [] xstr[m];
     delete [] estr[m];
     delete [] sstr[m];
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 int FixWall::setmask()
 {
   int mask = 0;
 
   // FLD implicit needs to invoke wall forces before pair style
 
   if (fldflag) mask |= PRE_FORCE;
   else mask |= POST_FORCE;
 
   mask |= THERMO_ENERGY;
   mask |= POST_FORCE_RESPA;
   mask |= MIN_POST_FORCE;
   return mask;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixWall::init()
 {
   for (int m = 0; m < nwall; m++) {
     if (xstyle[m] == VARIABLE) {
       xindex[m] = input->variable->find(xstr[m]);
       if (xindex[m] < 0)
         error->all(FLERR,"Variable name for fix wall does not exist");
       if (!input->variable->equalstyle(xindex[m]))
         error->all(FLERR,"Variable for fix wall is invalid style");
     }
     if (estyle[m] == VARIABLE) {
       eindex[m] = input->variable->find(estr[m]);
       if (eindex[m] < 0)
         error->all(FLERR,"Variable name for fix wall does not exist");
       if (!input->variable->equalstyle(eindex[m]))
         error->all(FLERR,"Variable for fix wall is invalid style");
     }
     if (sstyle[m] == VARIABLE) {
       sindex[m] = input->variable->find(sstr[m]);
       if (sindex[m] < 0)
         error->all(FLERR,"Variable name for fix wall does not exist");
       if (!input->variable->equalstyle(sindex[m]))
         error->all(FLERR,"Variable for fix wall is invalid style");
     }
   }
 
   // setup coefficients
 
   for (int m = 0; m < nwall; m++) precompute(m);
 
   if (strstr(update->integrate_style,"respa")) {
     ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
     if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixWall::setup(int vflag)
 {
   if (strstr(update->integrate_style,"verlet")) {
     if (!fldflag) post_force(vflag);
   } else {
     ((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
     post_force_respa(vflag,ilevel_respa,0);
     ((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixWall::min_setup(int vflag)
 {
   post_force(vflag);
 }
 
 /* ----------------------------------------------------------------------
    only called if fldflag set, in place of post_force
 ------------------------------------------------------------------------- */
 
 void FixWall::pre_force(int vflag)
 {
   post_force(vflag);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixWall::post_force(int vflag)
 {
   eflag = 0;
   for (int m = 0; m <= nwall; m++) ewall[m] = 0.0;
 
   // coord = current position of wall
   // evaluate variables if necessary, wrap with clear/add
   // for epsilon/sigma variables need to re-invoke precompute()
 
   if (varflag) modify->clearstep_compute();
 
   double coord;
   for (int m = 0; m < nwall; m++) {
     if (xstyle[m] == VARIABLE) {
       coord = input->variable->compute_equal(xindex[m]);
       if (wallwhich[m] < YLO) coord *= xscale;
       else if (wallwhich[m] < ZLO) coord *= yscale;
       else coord *= zscale;
     } else coord = coord0[m];
     if (wstyle[m] == VARIABLE) {
       if (estyle[m] == VARIABLE) {
         epsilon[m] = input->variable->compute_equal(eindex[m]);
         if (epsilon[m] < 0.0)
           error->all(FLERR,"Variable evaluation in fix wall gave bad value");
       }
       if (sstyle[m] == VARIABLE) {
         sigma[m] = input->variable->compute_equal(sindex[m]);
         if (sigma[m] < 0.0)
           error->all(FLERR,"Variable evaluation in fix wall gave bad value");
       }
       precompute(m);
     }
 
     wall_particle(m,wallwhich[m],coord);
   }
 
   if (varflag) modify->addstep_compute(update->ntimestep + 1);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixWall::post_force_respa(int vflag, int ilevel, int iloop)
 {
   if (ilevel == ilevel_respa) post_force(vflag);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixWall::min_post_force(int vflag)
 {
   post_force(vflag);
 }
 
 /* ----------------------------------------------------------------------
    energy of wall interaction
 ------------------------------------------------------------------------- */
 
 double FixWall::compute_scalar()
 {
   // only sum across procs one time
 
   if (eflag == 0) {
     MPI_Allreduce(ewall,ewall_all,nwall+1,MPI_DOUBLE,MPI_SUM,world);
     eflag = 1;
   }
   return ewall_all[0];
 }
 
 /* ----------------------------------------------------------------------
    components of force on wall
 ------------------------------------------------------------------------- */
 
 double FixWall::compute_vector(int n)
 {
   // only sum across procs one time
 
   if (eflag == 0) {
     MPI_Allreduce(ewall,ewall_all,nwall+1,MPI_DOUBLE,MPI_SUM,world);
     eflag = 1;
   }
   return ewall_all[n+1];
 }
diff --git a/src/fix_wall_reflect.cpp b/src/fix_wall_reflect.cpp
index b5c70037b..dba7cf703 100644
--- a/src/fix_wall_reflect.cpp
+++ b/src/fix_wall_reflect.cpp
@@ -1,228 +1,229 @@
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
 
    Copyright (2003) Sandia Corporation.  Under the terms of Contract
    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
    certain rights in this software.  This software is distributed under
    the GNU General Public License.
 
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
 #include <stdlib.h>
 #include <string.h>
 #include "fix_wall_reflect.h"
 #include "atom.h"
 #include "comm.h"
 #include "update.h"
 #include "modify.h"
 #include "domain.h"
 #include "lattice.h"
 #include "input.h"
 #include "variable.h"
 #include "error.h"
 #include "force.h"
 
 using namespace LAMMPS_NS;
 using namespace FixConst;
 
 enum{XLO=0,XHI=1,YLO=2,YHI=3,ZLO=4,ZHI=5};
 enum{NONE=0,EDGE,CONSTANT,VARIABLE};
 
 /* ---------------------------------------------------------------------- */
 
 FixWallReflect::FixWallReflect(LAMMPS *lmp, int narg, char **arg) :
-  Fix(lmp, narg, arg)
+  Fix(lmp, narg, arg),
+  nwall(0)
 {
   if (narg < 4) error->all(FLERR,"Illegal fix wall/reflect command");
 
   // parse args
 
   nwall = 0;
   int scaleflag = 1;
 
   int iarg = 3;
   while (iarg < narg) {
     if ((strcmp(arg[iarg],"xlo") == 0) || (strcmp(arg[iarg],"xhi") == 0) ||
         (strcmp(arg[iarg],"ylo") == 0) || (strcmp(arg[iarg],"yhi") == 0) ||
         (strcmp(arg[iarg],"zlo") == 0) || (strcmp(arg[iarg],"zhi") == 0)) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix wall/reflect command");
 
       int newwall;
       if (strcmp(arg[iarg],"xlo") == 0) newwall = XLO;
       else if (strcmp(arg[iarg],"xhi") == 0) newwall = XHI;
       else if (strcmp(arg[iarg],"ylo") == 0) newwall = YLO;
       else if (strcmp(arg[iarg],"yhi") == 0) newwall = YHI;
       else if (strcmp(arg[iarg],"zlo") == 0) newwall = ZLO;
       else if (strcmp(arg[iarg],"zhi") == 0) newwall = ZHI;
 
       for (int m = 0; (m < nwall) && (m < 6); m++)
         if (newwall == wallwhich[m])
           error->all(FLERR,"Wall defined twice in fix wall/reflect command");
 
       wallwhich[nwall] = newwall;
       if (strcmp(arg[iarg+1],"EDGE") == 0) {
         wallstyle[nwall] = EDGE;
         int dim = wallwhich[nwall] / 2;
         int side = wallwhich[nwall] % 2;
         if (side == 0) coord0[nwall] = domain->boxlo[dim];
         else coord0[nwall] = domain->boxhi[dim];
       } else if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
         wallstyle[nwall] = VARIABLE;
         int n = strlen(&arg[iarg+1][2]) + 1;
         varstr[nwall] = new char[n];
         strcpy(varstr[nwall],&arg[iarg+1][2]);
       } else {
         wallstyle[nwall] = CONSTANT;
         coord0[nwall] = force->numeric(FLERR,arg[iarg+1]);
       }
 
       nwall++;
       iarg += 2;
 
     } else if (strcmp(arg[iarg],"units") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal wall/reflect command");
       if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0;
       else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1;
       else error->all(FLERR,"Illegal fix wall/reflect command");
       iarg += 2;
     } else error->all(FLERR,"Illegal fix wall/reflect command");
   }
 
   // error check
 
   if (nwall == 0) error->all(FLERR,"Illegal fix wall command");
 
   for (int m = 0; m < nwall; m++) {
     if ((wallwhich[m] == XLO || wallwhich[m] == XHI) && domain->xperiodic)
       error->all(FLERR,"Cannot use fix wall/reflect in periodic dimension");
     if ((wallwhich[m] == YLO || wallwhich[m] == YHI) && domain->yperiodic)
       error->all(FLERR,"Cannot use fix wall/reflect in periodic dimension");
     if ((wallwhich[m] == ZLO || wallwhich[m] == ZHI) && domain->zperiodic)
       error->all(FLERR,"Cannot use fix wall/reflect in periodic dimension");
   }
 
   for (int m = 0; m < nwall; m++)
     if ((wallwhich[m] == ZLO || wallwhich[m] == ZHI) && domain->dimension == 2)
       error->all(FLERR,
                  "Cannot use fix wall/reflect zlo/zhi for a 2d simulation");
 
   // scale factors for CONSTANT and VARIABLE walls
 
   int flag = 0;
   for (int m = 0; m < nwall; m++)
     if (wallstyle[m] != EDGE) flag = 1;
 
   if (flag) {
     if (scaleflag) {
       xscale = domain->lattice->xlattice;
       yscale = domain->lattice->ylattice;
       zscale = domain->lattice->zlattice;
     }
     else xscale = yscale = zscale = 1.0;
 
     for (int m = 0; m < nwall; m++) {
       if (wallstyle[m] != CONSTANT) continue;
       if (wallwhich[m] < YLO) coord0[m] *= xscale;
       else if (wallwhich[m] < ZLO) coord0[m] *= yscale;
       else coord0[m] *= zscale;
     }
   }
 
   // set varflag if any wall positions are variable
 
   varflag = 0;
   for (int m = 0; m < nwall; m++)
     if (wallstyle[m] == VARIABLE) varflag = 1;
 }
 
 /* ---------------------------------------------------------------------- */
 
 FixWallReflect::~FixWallReflect()
 {
   if (copymode) return;
 
   for (int m = 0; m < nwall; m++)
     if (wallstyle[m] == VARIABLE) delete [] varstr[m];
 }
 
 /* ---------------------------------------------------------------------- */
 
 int FixWallReflect::setmask()
 {
   int mask = 0;
   mask |= POST_INTEGRATE;
   mask |= POST_INTEGRATE_RESPA;
   return mask;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixWallReflect::init()
 {
   for (int m = 0; m < nwall; m++) {
     if (wallstyle[m] != VARIABLE) continue;
     varindex[m] = input->variable->find(varstr[m]);
     if (varindex[m] < 0)
       error->all(FLERR,"Variable name for fix wall/reflect does not exist");
     if (!input->variable->equalstyle(varindex[m]))
       error->all(FLERR,"Variable for fix wall/reflect is invalid style");
   }
 
   int nrigid = 0;
   for (int i = 0; i < modify->nfix; i++)
     if (modify->fix[i]->rigid_flag) nrigid++;
 
   if (nrigid && comm->me == 0)
     error->warning(FLERR,"Should not allow rigid bodies to bounce off "
                    "relecting walls");
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixWallReflect::post_integrate()
 {
   int i,dim,side;
   double coord;
 
   // coord = current position of wall
   // evaluate variable if necessary, wrap with clear/add
 
   double **x = atom->x;
   double **v = atom->v;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
   if (varflag) modify->clearstep_compute();
 
   for (int m = 0; m < nwall; m++) {
     if (wallstyle[m] == VARIABLE) {
       coord = input->variable->compute_equal(varindex[m]);
       if (wallwhich[m] < YLO) coord *= xscale;
       else if (wallwhich[m] < ZLO) coord *= yscale;
       else coord *= zscale;
     } else coord = coord0[m];
 
     dim = wallwhich[m] / 2;
     side = wallwhich[m] % 2;
 
     for (i = 0; i < nlocal; i++)
       if (mask[i] & groupbit) {
         if (side == 0) {
           if (x[i][dim] < coord) {
             x[i][dim] = coord + (coord - x[i][dim]);
             v[i][dim] = -v[i][dim];
           }
         } else {
           if (x[i][dim] > coord) {
             x[i][dim] = coord - (x[i][dim] - coord);
             v[i][dim] = -v[i][dim];
           }
         }
       }
   }
 
   if (varflag) modify->addstep_compute(update->ntimestep + 1);
 }
diff --git a/src/fix_wall_region.cpp b/src/fix_wall_region.cpp
index 824f92eb8..f543c1217 100644
--- a/src/fix_wall_region.cpp
+++ b/src/fix_wall_region.cpp
@@ -1,367 +1,368 @@
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
 
    Copyright (2003) Sandia Corporation.  Under the terms of Contract
    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
    certain rights in this software.  This software is distributed under
    the GNU General Public License.
 
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
 #include <math.h>
 #include <stdlib.h>
 #include <string.h>
 #include "fix_wall_region.h"
 #include "atom.h"
 #include "atom_vec.h"
 #include "domain.h"
 #include "region.h"
 #include "force.h"
 #include "lattice.h"
 #include "update.h"
 #include "output.h"
 #include "respa.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
 using namespace FixConst;
 
 enum{LJ93,LJ126,COLLOID,HARMONIC};
 
 /* ---------------------------------------------------------------------- */
 
 FixWallRegion::FixWallRegion(LAMMPS *lmp, int narg, char **arg) :
-  Fix(lmp, narg, arg)
+  Fix(lmp, narg, arg),
+  idregion(NULL)
 {
   if (narg != 8) error->all(FLERR,"Illegal fix wall/region command");
 
   scalar_flag = 1;
   vector_flag = 1;
   size_vector = 3;
   global_freq = 1;
   extscalar = 1;
   extvector = 1;
   respa_level_support = 1;
   ilevel_respa = 0;
 
   // parse args
 
   iregion = domain->find_region(arg[3]);
   if (iregion == -1)
     error->all(FLERR,"Region ID for fix wall/region does not exist");
   int n = strlen(arg[3]) + 1;
   idregion = new char[n];
   strcpy(idregion,arg[3]);
 
   if (strcmp(arg[4],"lj93") == 0) style = LJ93;
   else if (strcmp(arg[4],"lj126") == 0) style = LJ126;
   else if (strcmp(arg[4],"colloid") == 0) style = COLLOID;
   else if (strcmp(arg[4],"harmonic") == 0) style = HARMONIC;
   else error->all(FLERR,"Illegal fix wall/region command");
 
   epsilon = force->numeric(FLERR,arg[5]);
   sigma = force->numeric(FLERR,arg[6]);
   cutoff = force->numeric(FLERR,arg[7]);
 
   if (cutoff <= 0.0) error->all(FLERR,"Fix wall/region cutoff <= 0.0");
 
   eflag = 0;
   ewall[0] = ewall[1] = ewall[2] = ewall[3] = 0.0;
 }
 
 /* ---------------------------------------------------------------------- */
 
 FixWallRegion::~FixWallRegion()
 {
   delete [] idregion;
 }
 
 /* ---------------------------------------------------------------------- */
 
 int FixWallRegion::setmask()
 {
   int mask = 0;
   mask |= POST_FORCE;
   mask |= THERMO_ENERGY;
   mask |= POST_FORCE_RESPA;
   mask |= MIN_POST_FORCE;
   return mask;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixWallRegion::init()
 {
   // set index and check validity of region
 
   iregion = domain->find_region(idregion);
   if (iregion == -1)
     error->all(FLERR,"Region ID for fix wall/region does not exist");
 
   // error checks for style COLLOID
   // insure all particles in group are extended particles
 
   if (style == COLLOID) {
     if (!atom->sphere_flag)
       error->all(FLERR,"Fix wall/region colloid requires atom style sphere");
 
     double *radius = atom->radius;
     int *mask = atom->mask;
     int nlocal = atom->nlocal;
 
     int flag = 0;
     for (int i = 0; i < nlocal; i++)
       if (mask[i] & groupbit)
         if (radius[i] == 0.0) flag = 1;
 
     int flagall;
     MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
     if (flagall)
       error->all(FLERR,"Fix wall/region colloid requires extended particles");
   }
 
   // setup coefficients for each style
 
   if (style == LJ93) {
     coeff1 = 6.0/5.0 * epsilon * pow(sigma,9.0);
     coeff2 = 3.0 * epsilon * pow(sigma,3.0);
     coeff3 = 2.0/15.0 * epsilon * pow(sigma,9.0);
     coeff4 = epsilon * pow(sigma,3.0);
     double rinv = 1.0/cutoff;
     double r2inv = rinv*rinv;
     double r4inv = r2inv*r2inv;
     offset = coeff3*r4inv*r4inv*rinv - coeff4*r2inv*rinv;
   } else if (style == LJ126) {
     coeff1 = 48.0 * epsilon * pow(sigma,12.0);
     coeff2 = 24.0 * epsilon * pow(sigma,6.0);
     coeff3 = 4.0 * epsilon * pow(sigma,12.0);
     coeff4 = 4.0 * epsilon * pow(sigma,6.0);
     double r2inv = 1.0/(cutoff*cutoff);
     double r6inv = r2inv*r2inv*r2inv;
     offset = r6inv*(coeff3*r6inv - coeff4);
   } else if (style == COLLOID) {
     coeff1 = -4.0/315.0 * epsilon * pow(sigma,6.0);
     coeff2 = -2.0/3.0 * epsilon;
     coeff3 = epsilon * pow(sigma,6.0)/7560.0;
     coeff4 = epsilon/6.0;
     double rinv = 1.0/cutoff;
     double r2inv = rinv*rinv;
     double r4inv = r2inv*r2inv;
     offset = coeff3*r4inv*r4inv*rinv - coeff4*r2inv*rinv;
   }
 
   if (strstr(update->integrate_style,"respa")) {
     ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
     if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixWallRegion::setup(int vflag)
 {
   if (strstr(update->integrate_style,"verlet"))
     post_force(vflag);
   else {
     ((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
     post_force_respa(vflag,ilevel_respa,0);
     ((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixWallRegion::min_setup(int vflag)
 {
   post_force(vflag);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixWallRegion::post_force(int vflag)
 {
   int i,m,n;
   double rinv,fx,fy,fz,tooclose;
 
   double **x = atom->x;
   double **f = atom->f;
   double *radius = atom->radius;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
   Region *region = domain->regions[iregion];
   region->prematch();
 
   int onflag = 0;
 
   // region->match() insures particle is in region or on surface, else error
   // if returned contact dist r = 0, is on surface, also an error
   // in COLLOID case, r <= radius is an error
   // initilize ewall after region->prematch(),
   //   so a dynamic region can access last timestep values
 
   eflag = 0;
   ewall[0] = ewall[1] = ewall[2] = ewall[3] = 0.0;
 
   for (i = 0; i < nlocal; i++)
     if (mask[i] & groupbit) {
       if (!region->match(x[i][0],x[i][1],x[i][2])) {
         onflag = 1;
         continue;
       }
       if (style == COLLOID) tooclose = radius[i];
       else tooclose = 0.0;
 
       n = region->surface(x[i][0],x[i][1],x[i][2],cutoff);
 
       for (m = 0; m < n; m++) {
         if (region->contact[m].r <= tooclose) {
           onflag = 1;
           continue;
         } else rinv = 1.0/region->contact[m].r;
 
         if (style == LJ93) lj93(region->contact[m].r);
         else if (style == LJ126) lj126(region->contact[m].r);
         else if (style == COLLOID) colloid(region->contact[m].r,radius[i]);
         else harmonic(region->contact[m].r);
 
         ewall[0] += eng;
         fx = fwall * region->contact[m].delx * rinv;
         fy = fwall * region->contact[m].dely * rinv;
         fz = fwall * region->contact[m].delz * rinv;
         f[i][0] += fx;
         f[i][1] += fy;
         f[i][2] += fz;
         ewall[1] -= fx;
         ewall[2] -= fy;
         ewall[3] -= fz;
       }
     }
 
   if (onflag) error->one(FLERR,"Particle outside surface of region "
                          "used in fix wall/region");
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixWallRegion::post_force_respa(int vflag, int ilevel, int iloop)
 {
   if (ilevel == ilevel_respa) post_force(vflag);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixWallRegion::min_post_force(int vflag)
 {
   post_force(vflag);
 }
 
 /* ----------------------------------------------------------------------
    energy of wall interaction
 ------------------------------------------------------------------------- */
 
 double FixWallRegion::compute_scalar()
 {
   // only sum across procs one time
 
   if (eflag == 0) {
     MPI_Allreduce(ewall,ewall_all,4,MPI_DOUBLE,MPI_SUM,world);
     eflag = 1;
   }
   return ewall_all[0];
 }
 
 /* ----------------------------------------------------------------------
    components of force on wall
 ------------------------------------------------------------------------- */
 
 double FixWallRegion::compute_vector(int n)
 {
   // only sum across procs one time
 
   if (eflag == 0) {
     MPI_Allreduce(ewall,ewall_all,4,MPI_DOUBLE,MPI_SUM,world);
     eflag = 1;
   }
   return ewall_all[n+1];
 }
 
 /* ----------------------------------------------------------------------
    LJ 9/3 interaction for particle with wall
    compute eng and fwall = magnitude of wall force
 ------------------------------------------------------------------------- */
 
 void FixWallRegion::lj93(double r)
 {
   double rinv = 1.0/r;
   double r2inv = rinv*rinv;
   double r4inv = r2inv*r2inv;
   double r10inv = r4inv*r4inv*r2inv;
   fwall = coeff1*r10inv - coeff2*r4inv;
   eng = coeff3*r4inv*r4inv*rinv - coeff4*r2inv*rinv - offset;
 }
 
 /* ----------------------------------------------------------------------
    LJ 12/6 interaction for particle with wall
    compute eng and fwall = magnitude of wall force
 ------------------------------------------------------------------------- */
 
 void FixWallRegion::lj126(double r)
 {
   double rinv = 1.0/r;
   double r2inv = rinv*rinv;
   double r6inv = r2inv*r2inv*r2inv;
   fwall = r6inv*(coeff1*r6inv - coeff2) * rinv;
   eng = r6inv*(coeff3*r6inv - coeff4) - offset;
 }
 
 /* ----------------------------------------------------------------------
    colloid interaction for finite-size particle of rad with wall
    compute eng and fwall = magnitude of wall force
 ------------------------------------------------------------------------- */
 
 void FixWallRegion::colloid(double r, double rad)
 {
   double new_coeff2 = coeff2*rad*rad*rad;
   double diam = 2.0*rad;
 
   double rad2 = rad*rad;
   double rad4 = rad2*rad2;
   double rad8 = rad4*rad4;
   double delta2 = rad2 - r*r;
   double rinv = 1.0/delta2;
   double r2inv = rinv*rinv;
   double r4inv = r2inv*r2inv;
   double r8inv = r4inv*r4inv;
   fwall = coeff1*(rad8*rad + 27.0*rad4*rad2*rad*pow(r,2.0)
                   + 63.0*rad4*rad*pow(r,4.0)
                   + 21.0*rad2*rad*pow(r,6.0))*r8inv - new_coeff2*r2inv;
 
   double r2 = 0.5*diam - r;
   double rinv2 = 1.0/r2;
   double r2inv2 = rinv2*rinv2;
   double r4inv2 = r2inv2*r2inv2;
   double r3 = r + 0.5*diam;
   double rinv3 = 1.0/r3;
   double r2inv3 = rinv3*rinv3;
   double r4inv3 = r2inv3*r2inv3;
   eng = coeff3*((-3.5*diam+r)*r4inv2*r2inv2*rinv2
                 + (3.5*diam+r)*r4inv3*r2inv3*rinv3) -
     coeff4*((-diam*r+r2*r3*(log(-r2)-log(r3)))*
             (-rinv2)*rinv3) - offset;
 }
 
 /* ----------------------------------------------------------------------
    harmonic interaction for particle with wall
    compute eng and fwall = magnitude of wall force
 ------------------------------------------------------------------------- */
 
 void FixWallRegion::harmonic(double r)
 {
   double dr = cutoff - r;
   fwall = 2.0*epsilon*dr;
   eng = epsilon*dr*dr;
 }