diff --git a/src/GRANULAR/fix_wall_gran.cpp b/src/GRANULAR/fix_wall_gran.cpp
index 52b801ad1b..fadddc79d8 100644
--- a/src/GRANULAR/fix_wall_gran.cpp
+++ b/src/GRANULAR/fix_wall_gran.cpp
@@ -46,8 +46,6 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
 {
   if (narg < 10) error->all("Illegal fix wall/gran command");
 
-  time_depend = 1;
-
   if (!atom->radius_flag || !atom->omega_flag || !atom->torque_flag)
     error->all("Fix wall/gran requires atom attributes radius, omega, torque");
 
@@ -166,7 +164,6 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
   if (wiggle) {
     double PI = 4.0 * atan(1.0);
     omega = 2.0*PI / period;
-    time_origin = update->ntimestep;
   }
 
   // perform initial allocation of atom-based arrays
@@ -254,7 +251,7 @@ void FixWallGran::post_force(int vflag)
   double whi = hi;
   vwall[0] = vwall[1] = vwall[2] = 0.0;
   if (wiggle) {
-    double arg = omega * (update->ntimestep - time_origin) * dt;
+    double arg = omega * (update->ntimestep - update->beginstep) * dt;
     wlo = lo + amplitude - amplitude*cos(arg);
     whi = hi + amplitude - amplitude*cos(arg);
     vwall[axis] = amplitude*omega*sin(arg);
diff --git a/src/PERI/compute_damage_atom.cpp b/src/PERI/compute_damage_atom.cpp
index 8df2ec603a..b8bbc935d8 100644
--- a/src/PERI/compute_damage_atom.cpp
+++ b/src/PERI/compute_damage_atom.cpp
@@ -81,7 +81,7 @@ void ComputeDamageAtom::compute_peratom()
     memory->sfree(damage);
     nmax = atom->nmax;
     damage = (double *) memory->smalloc(nmax*sizeof(double),
-					"compute/damage/atom:damage");
+					"damage/atom:damage");
     scalar_atom = damage;
   }
 
diff --git a/src/compute_centro_atom.cpp b/src/compute_centro_atom.cpp
index c534b91541..e23ff74de9 100644
--- a/src/compute_centro_atom.cpp
+++ b/src/compute_centro_atom.cpp
@@ -100,7 +100,7 @@ void ComputeCentroAtom::compute_peratom()
     memory->sfree(centro);
     nmax = atom->nmax;
     centro = (double *) 
-      memory->smalloc(nmax*sizeof(double),"compute/centro:centro");
+      memory->smalloc(nmax*sizeof(double),"centro/atom:centro");
     scalar_atom = centro;
   }
 
@@ -137,9 +137,9 @@ void ComputeCentroAtom::compute_peratom()
 	memory->sfree(nearest);
 	maxneigh = jnum;
 	distsq = (double *) memory->smalloc(maxneigh*sizeof(double),
-					    "compute/centro:distsq");
+					    "centro/atom:distsq");
 	nearest = (int *) memory->smalloc(maxneigh*sizeof(int),
-					  "compute/centro:nearest");
+					  "centro/atom:nearest");
       }
 
       // loop over list of all neighbors within force cutoff
diff --git a/src/compute_coord_atom.cpp b/src/compute_coord_atom.cpp
index 657fe3f985..c74f6fb6c2 100644
--- a/src/compute_coord_atom.cpp
+++ b/src/compute_coord_atom.cpp
@@ -101,7 +101,7 @@ void ComputeCoordAtom::compute_peratom()
     memory->sfree(coordination);
     nmax = atom->nmax;
     coordination = (double *) 
-      memory->smalloc(nmax*sizeof(double),"compute/coord/atom:coordination");
+      memory->smalloc(nmax*sizeof(double),"coord/atom:coordination");
     scalar_atom = coordination;
   }
 
diff --git a/src/compute_displace_atom.cpp b/src/compute_displace_atom.cpp
index 5174e80843..ac7dcb54e6 100644
--- a/src/compute_displace_atom.cpp
+++ b/src/compute_displace_atom.cpp
@@ -82,7 +82,7 @@ void ComputeDisplaceAtom::compute_peratom()
     memory->destroy_2d_double_array(displace);
     nmax = atom->nmax;
     displace =
-      memory->create_2d_double_array(nmax,4,"compute/displace/atom:displace");
+      memory->create_2d_double_array(nmax,4,"displace/atom:displace");
     vector_atom = displace;
   }
 
diff --git a/src/compute_ke_atom.cpp b/src/compute_ke_atom.cpp
index 432f7d5298..766c45479f 100644
--- a/src/compute_ke_atom.cpp
+++ b/src/compute_ke_atom.cpp
@@ -66,7 +66,7 @@ void ComputeKEAtom::compute_peratom()
   if (atom->nlocal > nmax) {
     memory->sfree(ke);
     nmax = atom->nmax;
-    ke = (double *) memory->smalloc(nmax*sizeof(double),"compute/ke/atom:ke");
+    ke = (double *) memory->smalloc(nmax*sizeof(double),"ke/atom:ke");
     scalar_atom = ke;
   }
 
diff --git a/src/compute_pe_atom.cpp b/src/compute_pe_atom.cpp
index 87042ec430..980c869283 100755
--- a/src/compute_pe_atom.cpp
+++ b/src/compute_pe_atom.cpp
@@ -86,7 +86,7 @@ void ComputePEAtom::compute_peratom()
     memory->sfree(energy);
     nmax = atom->nmax;
     energy = (double *) 
-      memory->smalloc(nmax*sizeof(double),"compute/pe/atom:energy");
+      memory->smalloc(nmax*sizeof(double),"pe/atom:energy");
     scalar_atom = energy;
   }
 
diff --git a/src/compute_reduce.cpp b/src/compute_reduce.cpp
index 143fbbf09c..2b8f3f4224 100644
--- a/src/compute_reduce.cpp
+++ b/src/compute_reduce.cpp
@@ -336,7 +336,7 @@ double ComputeReduce::compute_one(int m)
       maxatom = atom->nmax;
       memory->sfree(varatom);
       varatom =	(double *) 
-	memory->smalloc(maxatom*sizeof(double),"compute/reduce:varatom");
+	memory->smalloc(maxatom*sizeof(double),"reduce:varatom");
     }
 
     input->variable->compute_atom(n,igroup,varatom,1,0);
diff --git a/src/compute_reduce_region.cpp b/src/compute_reduce_region.cpp
index 8c80d75392..61c0da7e13 100644
--- a/src/compute_reduce_region.cpp
+++ b/src/compute_reduce_region.cpp
@@ -125,7 +125,7 @@ double ComputeReduceRegion::compute_one(int m)
       maxatom = atom->nmax;
       memory->sfree(varatom);
       varatom =	(double *) 
-	memory->smalloc(maxatom*sizeof(double),"compute/reduce:varatom");
+	memory->smalloc(maxatom*sizeof(double),"reduce/region:varatom");
     }
 
     input->variable->compute_atom(n,igroup,varatom,1,0);
diff --git a/src/compute_stress_atom.cpp b/src/compute_stress_atom.cpp
index 63e5dc8b2c..48337ed252 100644
--- a/src/compute_stress_atom.cpp
+++ b/src/compute_stress_atom.cpp
@@ -100,7 +100,7 @@ void ComputeStressAtom::compute_peratom()
     memory->destroy_2d_double_array(stress);
     nmax = atom->nmax;
     stress =
-      memory->create_2d_double_array(nmax,6,"compute/stress/atom:stress");
+      memory->create_2d_double_array(nmax,6,"stress/atom:stress");
     vector_atom = stress;
   }
 
diff --git a/src/compute_temp_deform.cpp b/src/compute_temp_deform.cpp
index 5017dc527b..3475ecc441 100644
--- a/src/compute_temp_deform.cpp
+++ b/src/compute_temp_deform.cpp
@@ -227,7 +227,7 @@ void ComputeTempDeform::remove_bias_all()
     memory->destroy_2d_double_array(vbiasall);
     maxbias = atom->nmax;
     vbiasall = memory->create_2d_double_array(maxbias,3,
-					      "compute/temp:vbiasall");
+					      "temp/deform:vbiasall");
   }
 
   double lamda[3];
diff --git a/src/compute_temp_partial.cpp b/src/compute_temp_partial.cpp
index adecda74b0..03b6f948f6 100644
--- a/src/compute_temp_partial.cpp
+++ b/src/compute_temp_partial.cpp
@@ -190,7 +190,7 @@ void ComputeTempPartial::remove_bias_all()
     memory->destroy_2d_double_array(vbiasall);
     maxbias = atom->nmax;
     vbiasall = memory->create_2d_double_array(maxbias,3,
-					      "compute/temp:vbiasall");
+					      "temp/partial:vbiasall");
   }
 
   if (!xflag) {
diff --git a/src/compute_temp_profile.cpp b/src/compute_temp_profile.cpp
index 140e906e91..7e98a8145f 100644
--- a/src/compute_temp_profile.cpp
+++ b/src/compute_temp_profile.cpp
@@ -286,7 +286,7 @@ void ComputeTempProfile::remove_bias_all()
     memory->destroy_2d_double_array(vbiasall);
     maxbias = atom->nmax;
     vbiasall = memory->create_2d_double_array(maxbias,3,
-					      "compute/temp:vbiasall");
+					      "temp/profile:vbiasall");
   }
 
   int ibin;
diff --git a/src/compute_temp_ramp.cpp b/src/compute_temp_ramp.cpp
index 897b8b8744..b9e6cdff9a 100644
--- a/src/compute_temp_ramp.cpp
+++ b/src/compute_temp_ramp.cpp
@@ -254,7 +254,7 @@ void ComputeTempRamp::remove_bias_all()
     memory->destroy_2d_double_array(vbiasall);
     maxbias = atom->nmax;
     vbiasall = memory->create_2d_double_array(maxbias,3,
-					      "compute/temp:vbiasall");
+					      "temp/ramp:vbiasall");
   }
 
   double fraction;
diff --git a/src/compute_temp_region.cpp b/src/compute_temp_region.cpp
index 3cdaa0606f..03fabcda0f 100644
--- a/src/compute_temp_region.cpp
+++ b/src/compute_temp_region.cpp
@@ -182,7 +182,7 @@ void ComputeTempRegion::remove_bias_all()
     memory->destroy_2d_double_array(vbiasall);
     maxbias = atom->nmax;
     vbiasall = memory->create_2d_double_array(maxbias,3,
-					      "compute/temp:vbiasall");
+					      "temp/region:vbiasall");
   }
   
   Region *region = domain->regions[iregion];
diff --git a/src/domain.cpp b/src/domain.cpp
index 2ffa88d88b..1180dfb6c6 100644
--- a/src/domain.cpp
+++ b/src/domain.cpp
@@ -78,6 +78,7 @@ Domain::Domain(LAMMPS *lmp) : Pointers(lmp)
   h_ratelo[0] = h_ratelo[1] = h_ratelo[2] = 0.0;
   
   prd_lamda[0] = prd_lamda[1] = prd_lamda[2] = 1.0;
+  prd_half_lamda[0] = prd_half_lamda[1] = prd_half_lamda[2] = 0.5;
   boxlo_lamda[0] = boxlo_lamda[1] = boxlo_lamda[2] = 0.0;
   boxhi_lamda[0] = boxhi_lamda[1] = boxhi_lamda[2] = 1.0;
 
@@ -186,9 +187,9 @@ void Domain::set_global_box()
   h_inv[1] = 1.0/h[1];
   h_inv[2] = 1.0/h[2];
 
-  xprd_half = 0.5*xprd;
-  yprd_half = 0.5*yprd;
-  zprd_half = 0.5*zprd;
+  prd_half[0] = xprd_half = 0.5*xprd;
+  prd_half[1] = yprd_half = 0.5*yprd;
+  prd_half[2] = zprd_half = 0.5*zprd;
 
   if (triclinic) {
     h[3] = yz;
@@ -680,7 +681,7 @@ void Domain::remap(double *x, int &image)
    remap the point into the periodic box no matter how far away
    resulting coord must satisfy lo <= coord < hi
    MAX is important since coord - prd < lo can happen when coord = hi
-   for triclinic, point is converted to lamda coords (0-1) before doing remap
+   for triclinic, point is converted to lamda coords (0-1) before remap
 ------------------------------------------------------------------------- */
 
 void Domain::remap(double *x)
@@ -722,6 +723,83 @@ void Domain::remap(double *x)
   if (triclinic) lamda2x(coord,x);
 }
 
+/* ----------------------------------------------------------------------
+   remap xnew to be within half box length of xold
+   do it directly, not iteratively, in case is far away
+   for triclinic, both points are converted to lamda coords (0-1) before remap
+------------------------------------------------------------------------- */
+
+void Domain::remap_near(double *xnew, double *xold)
+{
+  int n;
+  double *coordnew,*coordold,*period,*half;
+  double lamdanew[3],lamdaold[3];
+
+  if (triclinic == 0) {
+    period = prd;
+    half = prd_half;
+    coordnew = xnew;
+    coordold = xold;
+  } else {
+    period = prd_lamda;
+    half = prd_half_lamda;
+    x2lamda(xnew,lamdanew);
+    coordnew = lamdanew;
+    x2lamda(xold,lamdaold);
+    coordold = lamdaold;
+  }
+
+  // iterative form
+  // if (xperiodic) {
+  //   while (xnew[0]-xold[0] > half[0]) xnew[0] -= period[0];
+  //   while (xold[0]-xnew[0] > half[0]) xnew[0] += period[0];
+  // }
+
+  if (xperiodic) {
+    if (coordnew[0]-coordold[0] > period[0]) {
+      n = static_cast<int> ((coordnew[0]-coordold[0])/period[0]);
+      xnew[0] -= n*period[0];
+    }
+    while (xnew[0]-xold[0] > half[0]) xnew[0] -= period[0];
+    if (xold[0]-xnew[0] > period[0]) {
+      n = static_cast<int> ((xold[0]-xnew[0])/period[0]);
+      xnew[0] += n*period[0];
+    }
+    while (xold[0]-xnew[0] > half[0]) xnew[0] += period[0];
+  }
+
+  if (yperiodic) {
+    if (coordnew[1]-coordold[1] > period[1]) {
+      n = static_cast<int> ((coordnew[1]-coordold[1])/period[1]);
+      xnew[1] -= n*period[1];
+    }
+    while (xnew[1]-xold[1] > half[1]) xnew[1] -= period[1];
+    if (xold[1]-xnew[1] > period[1]) {
+      n = static_cast<int> ((xold[1]-xnew[1])/period[1]);
+      xnew[1] += n*period[1];
+    }
+    while (xold[1]-xnew[1] > half[1]) xnew[1] += period[1];
+  }
+
+  if (zperiodic) {
+    if (coordnew[2]-coordold[2] > period[2]) {
+      n = static_cast<int> ((coordnew[2]-coordold[2])/period[2]);
+      xnew[2] -= n*period[2];
+    }
+    while (xnew[2]-xold[2] > half[2]) xnew[2] -= period[2];
+    if (xold[2]-xnew[2] > period[2]) {
+      n = static_cast<int> ((xold[2]-xnew[2])/period[2]);
+      xnew[2] += n*period[2];
+    }
+    while (xold[2]-xnew[2] > half[2]) xnew[2] += period[2];
+  }
+
+  if (triclinic) {
+    lamda2x(coordnew,xnew);
+    lamda2x(coordold,xold);
+  }
+}
+
 /* ----------------------------------------------------------------------
    unmap the point via image flags
    x overwritten with result, don't reset image flag
diff --git a/src/domain.h b/src/domain.h
index 5732670d66..fff4f37837 100644
--- a/src/domain.h
+++ b/src/domain.h
@@ -40,15 +40,18 @@ class Domain : protected Pointers {
   double xprd,yprd,zprd;                 // global box dimensions
   double xprd_half,yprd_half,zprd_half;  // half dimensions
   double prd[3];                         // array form of dimensions
+  double prd_half[3];                    // array form of half dimensions
 
                                          // triclinic box
-                                         // xprd,half,prd = same (as if untilt)
+                                         // xprd,xprd_half,prd,prd_half =
+                                         // same as if untilted
   double prd_lamda[3];                   // lamda box = (1,1,1)
+  double prd_half_lamda[3];              // lamda half box = (0.5,0.5,0.5)
 
   double boxlo[3],boxhi[3];              // orthogonal box global bounds
 
                                          // triclinic box
-                                         // boxlo/hi = same (as if untilted)
+                                         // boxlo/hi = same as if untilted
   double boxlo_lamda[3],boxhi_lamda[3];  // lamda box = (0,1)
   double boxlo_bound[3],boxhi_bound[3];  // bounding box of tilted domain
 
@@ -91,6 +94,7 @@ class Domain : protected Pointers {
   void pbc();
   void remap(double *, int &);
   void remap(double *);
+  void remap_near(double *, double *);
   void unmap(double *, int);
   void unmap(double *, int, double *);
   void minimum_image(double &, double &, double &);
diff --git a/src/fix_gravity.cpp b/src/fix_gravity.cpp
index 87dea960ba..c8e7bc7c0a 100644
--- a/src/fix_gravity.cpp
+++ b/src/fix_gravity.cpp
@@ -33,8 +33,6 @@ FixGravity::FixGravity(LAMMPS *lmp, int narg, char **arg) :
 {
   if (narg < 5) error->all("Illegal fix gravity command");
 
-  time_depend = 1;
-
   magnitude = atof(arg[3]);
 
   if (strcmp(arg[4],"chute") == 0) {
@@ -88,8 +86,6 @@ FixGravity::FixGravity(LAMMPS *lmp, int narg, char **arg) :
       zgrav = 0.0;
     }
   }
-
-  time_initial = update->ntimestep;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -138,15 +134,15 @@ void FixGravity::post_force(int vflag)
   if (style == GRADIENT) {
     if (domain->dimension == 3) {
       double phi_current = degree2rad * 
-	(phi + (update->ntimestep-time_initial)*dt*phigrad*360.0);
+	(phi + (update->ntimestep-update->beginstep)*dt*phigrad*360.0);
       double theta_current = degree2rad * 
-	(theta + (update->ntimestep-time_initial)*dt*thetagrad*360.0);
+	(theta + (update->ntimestep-update->beginstep)*dt*thetagrad*360.0);
       xgrav = sin(theta_current) * cos(phi_current);
       ygrav = sin(theta_current) * sin(phi_current);
       zgrav = cos(theta_current);
     } else {
       double theta_current = degree2rad * 
-	(theta + (update->ntimestep-time_initial)*dt*thetagrad*360.0);
+	(theta + (update->ntimestep-update->beginstep)*dt*thetagrad*360.0);
       xgrav = sin(theta_current);
       ygrav = cos(theta_current);
     }
diff --git a/src/fix_gravity.h b/src/fix_gravity.h
index 0933d27ae8..3d6dc42754 100644
--- a/src/fix_gravity.h
+++ b/src/fix_gravity.h
@@ -30,7 +30,7 @@ class FixGravity : public Fix {
   void post_force_respa(int, int, int);
 
  private:
-  int style,time_initial;
+  int style;
   double magnitude,dt;
   double phi,theta,phigrad,thetagrad;
   double xdir,ydir,zdir;
diff --git a/src/fix_move.cpp b/src/fix_move.cpp
new file mode 100644
index 0000000000..cd0e7f6889
--- /dev/null
+++ b/src/fix_move.cpp
@@ -0,0 +1,890 @@
+/* ----------------------------------------------------------------------
+   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 "string.h"
+#include "stdlib.h"
+#include "math.h"
+#include "fix_move.h"
+#include "atom.h"
+#include "group.h"
+#include "update.h"
+#include "force.h"
+#include "domain.h"
+#include "lattice.h"
+#include "comm.h"
+#include "respa.h"
+#include "input.h"
+#include "variable.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+enum{LINEAR,WIGGLE,ROTATE,VARIABLE};
+enum{EQUAL,ATOM};
+
+/* ---------------------------------------------------------------------- */
+
+FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) :
+  Fix(lmp, narg, arg)
+{
+  if (narg < 4) error->all("Illegal fix move command");
+
+  restart_peratom = 1;
+  peratom_flag = 1;
+  size_peratom = 3;
+  peratom_freq = 1;
+  time_integrate = 1;
+
+  // parse args
+
+  int iarg;
+  xvarstr = yvarstr = zvarstr = NULL;
+  vxvarstr = vyvarstr = vzvarstr = NULL;
+
+  if (strcmp(arg[3],"linear") == 0) {
+    if (narg < 7) error->all("Illegal fix move command");
+    iarg = 7;
+    mstyle = LINEAR;
+    if (strcmp(arg[4],"NULL") == 0) vxflag = 0;
+    else {
+      vxflag = 1;
+      vx = atof(arg[4]);
+    }
+    if (strcmp(arg[5],"NULL") == 0) vyflag = 0;
+    else {
+      vyflag = 1;
+      vy = atof(arg[5]);
+    }
+    if (strcmp(arg[6],"NULL") == 0) vzflag = 0;
+    else {
+      vzflag = 1;
+      vz = atof(arg[6]);
+    }
+
+  } else if (strcmp(arg[3],"wiggle") == 0) {
+    if (narg < 8) error->all("Illegal fix move command");
+    iarg = 8;
+    mstyle = WIGGLE;
+    if (strcmp(arg[4],"NULL") == 0) axflag = 0;
+    else {
+      axflag = 1;
+      ax = atof(arg[4]);
+    }
+    if (strcmp(arg[5],"NULL") == 0) ayflag = 0;
+    else {
+      ayflag = 1;
+      ay = atof(arg[5]);
+    }
+    if (strcmp(arg[6],"NULL") == 0) azflag = 0;
+    else {
+      azflag = 1;
+      az = atof(arg[6]);
+    }
+    period = atof(arg[7]);
+
+  } else if (strcmp(arg[3],"rotate") == 0) {
+    if (narg < 11) error->all("Illegal fix move command");
+    iarg = 11;
+    mstyle = ROTATE;
+    point[0] = atof(arg[4]);
+    point[1] = atof(arg[5]);
+    point[2] = atof(arg[6]);
+    axis[0] = atof(arg[7]);
+    axis[1] = atof(arg[8]);
+    axis[2] = atof(arg[9]);
+    period = atof(arg[10]);
+
+  } else if (strcmp(arg[3],"variable") == 0) {
+    if (narg < 10) error->all("Illegal fix move command");
+    iarg = 10;
+    mstyle = VARIABLE;
+    if (strcmp(arg[4],"NULL") == 0) xvarstr = NULL;
+    else if (strstr(arg[4],"v_") == arg[4]) {
+      int n = strlen(&arg[4][2]) + 1;
+      xvarstr = new char[n];
+      strcpy(xvarstr,&arg[4][2]);
+    } else error->all("Illegal fix move command");
+    if (strcmp(arg[5],"NULL") == 0) yvarstr = NULL;
+    else if (strstr(arg[5],"v_") == arg[5]) {
+      int n = strlen(&arg[5][2]) + 1;
+      yvarstr = new char[n];
+      strcpy(yvarstr,&arg[5][2]);
+    } else error->all("Illegal fix move command");
+    if (strcmp(arg[6],"NULL") == 0) zvarstr = NULL;
+    else if (strstr(arg[6],"v_") == arg[6]) {
+      int n = strlen(&arg[6][2]) + 1;
+      zvarstr = new char[n];
+      strcpy(zvarstr,&arg[6][2]);
+    } else error->all("Illegal fix move command");
+    if (strcmp(arg[7],"NULL") == 0) vxvarstr = NULL;
+    else if (strstr(arg[7],"v_") == arg[7]) {
+      int n = strlen(&arg[7][2]) + 1;
+      vxvarstr = new char[n];
+      strcpy(vxvarstr,&arg[7][2]);
+    } else error->all("Illegal fix move command");
+    if (strcmp(arg[8],"NULL") == 0) vyvarstr = NULL;
+    else if (strstr(arg[8],"v_") == arg[8]) {
+      int n = strlen(&arg[8][2]) + 1;
+      vyvarstr = new char[n];
+      strcpy(vyvarstr,&arg[8][2]);
+    } else error->all("Illegal fix move command");
+    if (strcmp(arg[9],"NULL") == 0) vzvarstr = NULL;
+    else if (strstr(arg[9],"v_") == arg[9]) {
+      int n = strlen(&arg[9][2]) + 1;
+      vzvarstr = new char[n];
+      strcpy(vzvarstr,&arg[9][2]);
+    } else error->all("Illegal fix move command");
+    
+  } else error->all("Illegal fix move command");
+
+  // optional args
+
+  int scaleflag = 1;
+
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"units") == 0) {
+      if (iarg+2 > narg) error->all("Illegal fix move command");
+      if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0;
+      else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1;
+      else error->all("Illegal fix move command");
+      iarg += 2;
+    } else error->all("Illegal fix move command");
+  }
+
+  // error checks and warnings
+
+  if (domain->dimension == 2) {
+    if (mstyle == LINEAR && vzflag && vz != 0.0)
+      error->all("Fix move cannot set linear z motion for 2d problem");
+    if (mstyle == WIGGLE && azflag && az != 0.0)
+      error->all("Fix move cannot set wiggle z motion for 2d problem");
+    if (mstyle == ROTATE && (axis[0] != 0.0 || axis[1] != 0.0))
+      error->all("Fix move cannot rotate aroung non z-axis for 2d problem");
+    if (mstyle == VARIABLE && (zvarstr || vzvarstr))
+      error->all("Fix move cannot define z or vz variable for 2d problem");
+  }
+
+  if (atom->angmom_flag && comm->me == 0)
+    error->warning("Fix move does not update angular momentum");
+  if (atom->quat_flag && comm->me == 0)
+    error->warning("Fix move does not update quaternions");
+
+  // setup scaling and apply scaling factors to velocity & amplitude
+
+  if ((mstyle == LINEAR || mstyle == WIGGLE || mstyle == ROTATE) && 
+      scaleflag) {
+    if (domain->lattice == NULL)
+      error->all("Use of fix move with undefined lattice");
+
+    double xscale,yscale,zscale;
+    if (scaleflag) {
+      xscale = domain->lattice->xlattice;
+      yscale = domain->lattice->ylattice;
+      zscale = domain->lattice->zlattice;
+    }
+    else xscale = yscale = zscale = 1.0;
+
+    if (mstyle == LINEAR) {
+      if (vxflag) vx *= xscale;
+      if (vyflag) vy *= yscale;
+      if (vzflag) vz *= zscale;
+    } else if (mstyle == WIGGLE) {
+      if (axflag) ax *= xscale;
+      if (ayflag) ay *= yscale;
+      if (azflag) az *= zscale;
+    } else if (mstyle == ROTATE) {
+      point[0] *= xscale;
+      point[1] *= yscale;
+      point[2] *= zscale;
+    }
+  }
+
+  // set omega_rotate from period
+
+  if (mstyle == WIGGLE || mstyle == ROTATE) {
+    double PI = 4.0 * atan(1.0);
+    omega_rotate = 2.0*PI / period;
+  }
+
+  // runit = unit vector along rotation axis
+
+  if (mstyle == ROTATE) {
+    double len = sqrt(axis[0]*axis[0] + axis[1]*axis[1] + axis[2]*axis[2]);
+    if (len == 0.0) 
+      error->all("Fix move cannot have 0 length rotation vector");
+    runit[0] = axis[0]/len;
+    runit[1] = axis[1]/len;
+    runit[2] = axis[2]/len;
+  }
+
+  // set omega_flag if particles store omega
+
+  omega_flag = atom->omega_flag;
+ 
+  // perform initial allocation of atom-based array
+  // register with Atom class
+
+  xoriginal = NULL;
+  grow_arrays(atom->nmax);
+  atom->add_callback(0);
+  atom->add_callback(1);
+
+  maxatom = 0;
+  displace = velocity = NULL;
+
+  // xoriginal = initial unwrapped positions of atoms
+
+  double **x = atom->x;
+  int *image = atom->image;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) domain->unmap(x[i],image[i],xoriginal[i]);
+    else xoriginal[i][0] = xoriginal[i][1] = xoriginal[i][2] = 0.0;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixMove::~FixMove()
+{
+  // unregister callbacks to this fix from Atom class
+ 
+  atom->delete_callback(id,0);
+  atom->delete_callback(id,1);
+
+  // delete locally stored arrays
+
+  memory->destroy_2d_double_array(xoriginal);
+  memory->destroy_2d_double_array(displace);
+  memory->destroy_2d_double_array(velocity);
+
+  delete [] xvarstr;
+  delete [] yvarstr;
+  delete [] zvarstr;
+  delete [] vxvarstr;
+  delete [] vyvarstr;
+  delete [] vzvarstr;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixMove::setmask()
+{
+  int mask = 0;
+  mask |= INITIAL_INTEGRATE;
+  mask |= INITIAL_INTEGRATE_RESPA;
+  mask |= FINAL_INTEGRATE;
+  mask |= FINAL_INTEGRATE_RESPA;
+  return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixMove::init()
+{
+  dt = update->dt;
+  dtv = update->dt;
+  dtf = 0.5 * update->dt * force->ftm2v;
+
+  // set indices and style of all variables
+
+  if (mstyle == VARIABLE) {
+    if (xvarstr) {
+      xvar = input->variable->find(xvarstr);
+      if (xvar < 0) error->all("Variable name for fix move does not exist");
+      if (input->variable->equalstyle(xvar)) xvarstyle = EQUAL;
+      else if (input->variable->atomstyle(xvar)) xvarstyle = ATOM;
+      else error->all("Variable for fix move is invalid style");
+    }
+    if (yvarstr) {
+      yvar = input->variable->find(yvarstr);
+      if (yvar < 0) error->all("Variable name for fix move does not exist");
+      if (input->variable->equalstyle(yvar)) yvarstyle = EQUAL;
+      else if (input->variable->atomstyle(yvar)) yvarstyle = ATOM;
+      else error->all("Variable for fix move is invalid style");
+    }
+    if (zvarstr) {
+      zvar = input->variable->find(zvarstr);
+      if (zvar < 0) error->all("Variable name for fix move does not exist");
+      if (input->variable->equalstyle(zvar)) zvarstyle = EQUAL;
+      else if (input->variable->atomstyle(zvar)) zvarstyle = ATOM;
+      else error->all("Variable for fix move is invalid style");
+    }
+    if (vxvarstr) {
+      vxvar = input->variable->find(vxvarstr);
+      if (vxvar < 0) error->all("Variable name for fix move does not exist");
+      if (input->variable->equalstyle(vxvar)) vxvarstyle = EQUAL;
+      else if (input->variable->atomstyle(vxvar)) vxvarstyle = ATOM;
+      else error->all("Variable for fix move is invalid style");
+    }
+    if (vyvarstr) {
+      vyvar = input->variable->find(vyvarstr);
+      if (vyvar < 0) error->all("Variable name for fix move does not exist");
+      if (input->variable->equalstyle(vyvar)) vyvarstyle = EQUAL;
+      else if (input->variable->atomstyle(vyvar)) vyvarstyle = ATOM;
+      else error->all("Variable for fix move is invalid style");
+    }
+    if (vzvarstr) {
+      vzvar = input->variable->find(vzvarstr);
+      if (vzvar < 0) error->all("Variable name for fix move does not exist");
+      if (input->variable->equalstyle(vzvar)) vzvarstyle = EQUAL;
+      else if (input->variable->atomstyle(vzvar)) vzvarstyle = ATOM;
+      else error->all("Variable for fix move is invalid style");
+    }
+
+    displaceflag = velocityflag = 0;
+    if (xvarstr && xvarstyle == ATOM) displaceflag = 1;
+    if (yvarstr && yvarstyle == ATOM) displaceflag = 1;
+    if (zvarstr && zvarstyle == ATOM) displaceflag = 1;
+    if (vxvarstr && vxvarstyle == ATOM) velocityflag = 1;
+    if (vyvarstr && vyvarstyle == ATOM) velocityflag = 1;
+    if (vzvarstr && vzvarstyle == ATOM) velocityflag = 1;
+  }
+
+  if (strcmp(update->integrate_style,"respa") == 0)
+    nlevels_respa = ((Respa *) update->integrate)->nlevels;
+}
+
+/* ----------------------------------------------------------------------
+   set x,v of particles
+------------------------------------------------------------------------- */
+
+void FixMove::initial_integrate(int vflag)
+{
+  double dtfm;
+  double xold[3],a[3],b[3],c[3],d[3],disp[3];
+  double x0dotr,dx,dy,dz;
+
+  double delta = (update->ntimestep - update->beginstep) * dt;
+
+  double **x = atom->x;
+  double **v = atom->v;
+  double **f = atom->f;
+  double **omega = atom->omega;
+  double *rmass = atom->rmass;
+  double *mass = atom->mass;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  // for linear: X = X0 + V*dt
+
+  if (mstyle == LINEAR) {
+    for (int i = 0; i < nlocal; i++) {
+      if (mask[i] & groupbit) {
+	xold[0] = x[i][0];
+	xold[1] = x[i][1];
+	xold[2] = x[i][2];
+
+	if (vxflag) {
+	  v[i][0] = vx;
+	  x[i][0] = xoriginal[i][0] + vx*delta;
+	} else if (rmass) {
+	  dtfm = dtf / rmass[i];
+	  v[i][0] += dtfm * f[i][0];
+	  x[i][0] += dtv * v[i][0];
+	} else {
+	  dtfm = dtf / mass[type[i]];
+	  v[i][0] += dtfm * f[i][0];
+	  x[i][0] += dtv * v[i][0];
+	}
+
+	if (vyflag) {
+	  v[i][1] = vy;
+	  x[i][1] = xoriginal[i][1] + vy*delta;
+	} else if (rmass) {
+	  dtfm = dtf / rmass[i];
+	  v[i][1] += dtfm * f[i][1];
+	  x[i][1] += dtv * v[i][1];
+	} else {
+	  dtfm = dtf / mass[type[i]];
+	  v[i][1] += dtfm * f[i][1];
+	  x[i][1] += dtv * v[i][1];
+	}
+
+	if (vzflag) {
+	  v[i][2] = vz;
+	  x[i][2] = xoriginal[i][2] + vz*delta;
+	} else if (rmass) {
+	  dtfm = dtf / rmass[i];
+	  v[i][2] += dtfm * f[i][2];
+	  x[i][2] += dtv * v[i][2];
+	} else {
+	  dtfm = dtf / mass[type[i]];
+	  v[i][2] += dtfm * f[i][2];
+	  x[i][2] += dtv * v[i][2];
+	}
+
+	domain->remap_near(x[i],xold);
+      }
+    }
+
+  // for wiggle: X = X0 + A sin(w*dt)
+
+  } else if (mstyle == WIGGLE) {
+    double arg = omega_rotate * delta;
+    double sine = sin(arg);
+    double cosine = cos(arg);
+
+    for (int i = 0; i < nlocal; i++) {
+      if (mask[i] & groupbit) {
+	xold[0] = x[i][0];
+	xold[1] = x[i][1];
+	xold[2] = x[i][2];
+
+	if (axflag) {
+	  v[i][0] = ax*omega_rotate*cosine;
+	  x[i][0] = xoriginal[i][0] + ax*sine;
+	} else if (rmass) {
+	  dtfm = dtf / rmass[i];
+	  v[i][0] += dtfm * f[i][0];
+	  x[i][0] += dtv * v[i][0];
+	} else {
+	  dtfm = dtf / mass[type[i]];
+	  v[i][0] += dtfm * f[i][0];
+	  x[i][0] += dtv * v[i][0];
+	}
+
+	if (ayflag) {
+	  v[i][1] = ay*omega_rotate*cosine;
+	  x[i][1] = xoriginal[i][1] + ay*sine;
+	} else if (rmass) {
+	  dtfm = dtf / rmass[i];
+	  v[i][1] += dtfm * f[i][1];
+	  x[i][1] += dtv * v[i][1];
+	} else {
+	  dtfm = dtf / mass[type[i]];
+	  v[i][1] += dtfm * f[i][1];
+	  x[i][1] += dtv * v[i][1];
+	}
+
+	if (azflag) {
+	  v[i][2] = az*omega_rotate*cosine;
+	  x[i][2] = xoriginal[i][2] + az*sine;
+	} else if (rmass) {
+	  dtfm = dtf / rmass[i];
+	  v[i][2] += dtfm * f[i][2];
+	  x[i][2] += dtv * v[i][2];
+	} else {
+	  dtfm = dtf / mass[type[i]];
+	  v[i][2] += dtfm * f[i][2];
+	  x[i][2] += dtv * v[i][2];
+	}
+
+	domain->remap_near(x[i],xold);
+      }
+    }
+
+  // for rotate:
+  // P = vector = point of rotation
+  // R = vector = axis of rotation
+  // w = omega of rotation (from period)
+  // X0 = initial coord of atom
+  // R0 = unit vector for R
+  // C = (X0 dot R0) R0 = projection of atom coord onto R
+  // D = X0 - P = vector from P to X0
+  // A = D - C = vector from R line to X0
+  // B = R0 cross A = vector perp to A in plane of rotation
+  // A,B define plane of circular rotation around R line
+  // X = P + C + A cos(w*dt) + B sin(w*dt)
+  // V = w R0 cross (A cos(w*dt) + B sin(w*dt))
+
+  } else if (mstyle == ROTATE) {
+    double arg = omega_rotate * delta;
+    double sine = sin(arg);
+    double cosine = cos(arg);
+
+    for (int i = 0; i < nlocal; i++) {
+      if (mask[i] & groupbit) {
+	xold[0] = x[i][0];
+	xold[1] = x[i][1];
+	xold[2] = x[i][2];
+
+	x0dotr = xoriginal[i][0]*runit[0] + 
+	  xoriginal[i][1]*runit[1] + xoriginal[i][2]*runit[2]; 
+	c[0] = x0dotr * runit[0];
+	c[1] = x0dotr * runit[1];
+	c[2] = x0dotr * runit[2];
+	d[0] = xoriginal[i][0] - point[0];
+	d[1] = xoriginal[i][1] - point[1];
+	d[2] = xoriginal[i][2] - point[2];
+	a[0] = d[0] - c[0];
+	a[1] = d[1] - c[1];
+	a[2] = d[2] - c[2];
+	b[0] = runit[1]*a[2] - runit[2]*a[1];
+	b[1] = runit[2]*a[0] - runit[0]*a[2];
+	b[2] = runit[0]*a[1] - runit[1]*a[0];
+	disp[0] = a[0]*cosine  + b[0]*sine;
+	disp[1] = a[1]*cosine  + b[1]*sine;
+	disp[2] = a[2]*cosine  + b[2]*sine;
+
+	x[i][0] = point[0] + c[0] + disp[0];
+	x[i][1] = point[1] + c[1] + disp[1];
+	x[i][2] = point[2] + c[2] + disp[2];
+	v[i][0] = omega_rotate * (runit[1]*disp[2] - runit[2]*disp[1]);
+	v[i][1] = omega_rotate * (runit[2]*disp[0] - runit[0]*disp[2]);
+	v[i][2] = omega_rotate * (runit[0]*disp[1] - runit[1]*disp[0]);
+	if (omega_flag) {
+	  omega[i][0] = omega_rotate*runit[0];
+	  omega[i][1] = omega_rotate*runit[1];
+	  omega[i][2] = omega_rotate*runit[2];
+	}
+
+	domain->remap_near(x[i],xold);
+      }
+    }
+
+  // for variable: compute x,v from variables
+
+  } else if (mstyle == VARIABLE) {
+
+    // reallocate displace and velocity arrays as necessary
+
+    if ((displaceflag || velocityflag) && nlocal > maxatom) {
+      maxatom = atom->nmax;
+      if (displaceflag) {
+	memory->destroy_2d_double_array(displace);
+	displace = memory->create_2d_double_array(maxatom,3,"move:displace");
+      }
+      if (velocityflag) {
+	memory->destroy_2d_double_array(velocity);
+	velocity = memory->create_2d_double_array(maxatom,3,"move:velocity");
+      }
+    }
+
+    // pre-compute variable values
+
+    if (xvarstr) {
+      if (xvarstyle == EQUAL) dx = input->variable->compute_equal(xvar);
+      else input->variable->compute_atom(xvar,igroup,&displace[0][0],3,0);
+    }
+    if (yvarstr) {
+      if (yvarstyle == EQUAL) dy = input->variable->compute_equal(yvar);
+      else input->variable->compute_atom(yvar,igroup,&displace[0][1],3,0);
+    }
+    if (zvarstr) {
+      if (zvarstyle == EQUAL) dz = input->variable->compute_equal(zvar);
+      else input->variable->compute_atom(zvar,igroup,&displace[0][2],3,0);
+    }
+    if (vxvarstr) {
+      if (vxvarstyle == EQUAL) vx = input->variable->compute_equal(vxvar);
+      else input->variable->compute_atom(vxvar,igroup,&velocity[0][0],3,0);
+    }
+    if (vyvarstr) {
+      if (vyvarstyle == EQUAL) vy = input->variable->compute_equal(vyvar);
+      else input->variable->compute_atom(vyvar,igroup,&velocity[0][1],3,0);
+    }
+    if (vzvarstr) {
+      if (vzvarstyle == EQUAL) vz = input->variable->compute_equal(vzvar);
+      else input->variable->compute_atom(vzvar,igroup,&velocity[0][2],3,0);
+    }
+
+    // update x,v
+
+    for (int i = 0; i < nlocal; i++) {
+      if (mask[i] & groupbit) {
+	xold[0] = x[i][0];
+	xold[1] = x[i][1];
+	xold[2] = x[i][2];
+
+	if (xvarstr && vxvarstr) {
+	  if (vxvarstyle == EQUAL) v[i][0] = vx;
+	  else v[i][0] = velocity[i][0];
+	  if (xvarstyle == EQUAL) x[i][0] = xoriginal[i][0] + dx;
+	  else x[i][0] = xoriginal[i][0] + displace[i][0];
+	} else if (xvarstr) {
+	  if (xvarstyle == EQUAL) x[i][0] = xoriginal[i][0] + dx;
+	  else x[i][0] = xoriginal[i][0] + displace[i][0];
+	} else if (vxvarstr) {
+	  if (vxvarstyle == EQUAL) v[i][0] = vx;
+	  else v[i][0] = velocity[i][0];
+	  if (rmass) {
+	    dtfm = dtf / rmass[i];
+	    x[i][0] += dtv * v[i][0];
+	  } else {
+	    dtfm = dtf / mass[type[i]];
+	    x[i][0] += dtv * v[i][0];
+	  }
+	} else {
+	  if (rmass) {
+	    dtfm = dtf / rmass[i];
+	    v[i][0] += dtfm * f[i][0];
+	    x[i][0] += dtv * v[i][0];
+	  } else {
+	    dtfm = dtf / mass[type[i]];
+	    v[i][0] += dtfm * f[i][0];
+	    x[i][0] += dtv * v[i][0];
+	  }
+	}
+
+	if (yvarstr && vyvarstr) {
+	  if (vyvarstyle == EQUAL) v[i][1] = vy;
+	  else v[i][1] = velocity[i][1];
+	  if (yvarstyle == EQUAL) x[i][1] = xoriginal[i][1] + dy;
+	  else x[i][1] = xoriginal[i][1] + displace[i][1];
+	} else if (yvarstr) {
+	  if (yvarstyle == EQUAL) x[i][1] = xoriginal[i][1] + dy;
+	  else x[i][1] = xoriginal[i][1] + displace[i][1];
+	} else if (vyvarstr) {
+	  if (vyvarstyle == EQUAL) v[i][1] = vy;
+	  else v[i][1] = velocity[i][1];
+	  if (rmass) {
+	    dtfm = dtf / rmass[i];
+	    x[i][1] += dtv * v[i][1];
+	  } else {
+	    dtfm = dtf / mass[type[i]];
+	    x[i][1] += dtv * v[i][1];
+	  }
+	} else {
+	  if (rmass) {
+	    dtfm = dtf / rmass[i];
+	    v[i][1] += dtfm * f[i][1];
+	    x[i][1] += dtv * v[i][1];
+	  } else {
+	    dtfm = dtf / mass[type[i]];
+	    v[i][1] += dtfm * f[i][1];
+	    x[i][1] += dtv * v[i][1];
+	  }
+	}
+
+	if (zvarstr && vzvarstr) {
+	  if (vzvarstyle == EQUAL) v[i][2] = vz;
+	  else v[i][2] = velocity[i][2];
+	  if (zvarstyle == EQUAL) x[i][2] = xoriginal[i][2] + dz;
+	  else x[i][2] = xoriginal[i][2] + displace[i][2];
+	} else if (zvarstr) {
+	  if (zvarstyle == EQUAL) x[i][2] = xoriginal[i][2] + dz;
+	  else x[i][2] = xoriginal[i][2] + displace[i][2];
+	} else if (vzvarstr) {
+	  if (vzvarstyle == EQUAL) v[i][2] = vz;
+	  else v[i][2] = velocity[i][2];
+	  if (rmass) {
+	    dtfm = dtf / rmass[i];
+	    x[i][2] += dtv * v[i][2];
+	  } else {
+	    dtfm = dtf / mass[type[i]];
+	    x[i][2] += dtv * v[i][2];
+	  }
+	} else {
+	  if (rmass) {
+	    dtfm = dtf / rmass[i];
+	    v[i][2] += dtfm * f[i][2];
+	    x[i][2] += dtv * v[i][2];
+	  } else {
+	    dtfm = dtf / mass[type[i]];
+	    v[i][2] += dtfm * f[i][2];
+	    x[i][2] += dtv * v[i][2];
+	  }
+	}
+
+	domain->remap_near(x[i],xold);
+      }
+    }
+  }
+}
+
+/* ----------------------------------------------------------------------
+   final NVE of particles with NULL components
+------------------------------------------------------------------------- */
+
+void FixMove::final_integrate()
+{
+  double dtfm;
+
+  int xflag = 1;
+  if (mstyle == LINEAR && vxflag) xflag = 0;
+  else if (mstyle == WIGGLE && axflag) xflag = 0;
+  else if (mstyle == ROTATE) xflag = 0;
+  else if (mstyle == VARIABLE && (xvarstr || vxvarstr)) xflag = 0;
+
+  int yflag = 1;
+  if (mstyle == LINEAR && vyflag) yflag = 0;
+  else if (mstyle == WIGGLE && ayflag) yflag = 0;
+  else if (mstyle == ROTATE) yflag = 0;
+  else if (mstyle == VARIABLE && (yvarstr || vyvarstr)) yflag = 0;
+
+  int zflag = 1;
+  if (mstyle == LINEAR && vzflag) zflag = 0;
+  else if (mstyle == WIGGLE && azflag) zflag = 0;
+  else if (mstyle == ROTATE) zflag = 0;
+  else if (mstyle == VARIABLE && (zvarstr || vzvarstr)) zflag = 0;
+
+  double **x = atom->x;
+  double **v = atom->v;
+  double **f = atom->f;
+  double **omega = atom->omega;
+  double *rmass = atom->rmass;
+  double *mass = atom->mass;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) {
+      if (xflag)
+	if (rmass) {
+	  dtfm = dtf / rmass[i];
+	  v[i][0] += dtfm * f[i][0];
+	} else {
+	  dtfm = dtf / mass[type[i]];
+	  v[i][0] += dtfm * f[i][0];
+	}
+      
+      if (yflag)
+	if (rmass) {
+	  dtfm = dtf / rmass[i];
+	  v[i][1] += dtfm * f[i][1];
+	} else {
+	  dtfm = dtf / mass[type[i]];
+	  v[i][1] += dtfm * f[i][1];
+	}
+      
+      if (zflag)
+	if (rmass) {
+	  dtfm = dtf / rmass[i];
+	  v[i][2] += dtfm * f[i][2];
+	} else {
+	  dtfm = dtf / mass[type[i]];
+	  v[i][2] += dtfm * f[i][2];
+	}
+    }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixMove::initial_integrate_respa(int vflag, int ilevel, int flag)
+{
+  if (flag) return;             // only used by NPT,NPH
+
+  // outermost level - update v and x
+  // all other levels - nothing
+
+  if (ilevel == nlevels_respa-1) initial_integrate(vflag);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixMove::final_integrate_respa(int ilevel)
+{
+  if (ilevel == nlevels_respa-1) final_integrate();
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of local atom-based array
+------------------------------------------------------------------------- */
+
+double FixMove::memory_usage()
+{
+  double bytes = atom->nmax*3 * sizeof(double);
+  return bytes;
+}
+
+/* ----------------------------------------------------------------------
+   allocate atom-based array
+------------------------------------------------------------------------- */
+
+void FixMove::grow_arrays(int nmax)
+{
+  xoriginal =
+    memory->grow_2d_double_array(xoriginal,nmax,3,"move:xoriginal");
+  vector_atom = xoriginal;
+}
+
+/* ----------------------------------------------------------------------
+   copy values within local atom-based array
+------------------------------------------------------------------------- */
+
+void FixMove::copy_arrays(int i, int j)
+{
+  xoriginal[j][0] = xoriginal[i][0];
+  xoriginal[j][1] = xoriginal[i][1];
+  xoriginal[j][2] = xoriginal[i][2];
+}
+
+/* ----------------------------------------------------------------------
+   pack values in local atom-based array for exchange with another proc
+------------------------------------------------------------------------- */
+
+int FixMove::pack_exchange(int i, double *buf)
+{
+  buf[0] = xoriginal[i][0];
+  buf[1] = xoriginal[i][1];
+  buf[2] = xoriginal[i][2];
+  return 3;
+}
+
+/* ----------------------------------------------------------------------
+   unpack values in local atom-based array from exchange with another proc
+------------------------------------------------------------------------- */
+
+int FixMove::unpack_exchange(int nlocal, double *buf)
+{
+  xoriginal[nlocal][0] = buf[0];
+  xoriginal[nlocal][1] = buf[1];
+  xoriginal[nlocal][2] = buf[2];
+  return 3;
+}
+
+/* ----------------------------------------------------------------------
+   pack values in local atom-based arrays for restart file
+------------------------------------------------------------------------- */
+
+int FixMove::pack_restart(int i, double *buf)
+{
+  buf[0] = 4;
+  buf[1] = xoriginal[i][0];
+  buf[2] = xoriginal[i][1];
+  buf[3] = xoriginal[i][2];
+  return 4;
+}
+
+/* ----------------------------------------------------------------------
+   unpack values from atom->extra array to restart the fix
+------------------------------------------------------------------------- */
+
+void FixMove::unpack_restart(int nlocal, int nth)
+{
+  double **extra = atom->extra;
+
+  // skip to Nth set of extra values
+
+  int m = 0;
+  for (int i = 0; i < nth; i++) m += static_cast<int> (extra[nlocal][m]);
+  m++;
+
+  xoriginal[nlocal][0] = extra[nlocal][m++];
+  xoriginal[nlocal][1] = extra[nlocal][m++];
+  xoriginal[nlocal][2] = extra[nlocal][m++];
+}
+
+/* ----------------------------------------------------------------------
+   maxsize of any atom's restart data
+------------------------------------------------------------------------- */
+
+int FixMove::maxsize_restart()
+{
+  return 4;
+}
+
+/* ----------------------------------------------------------------------
+   size of atom nlocal's restart data
+------------------------------------------------------------------------- */
+
+int FixMove::size_restart(int nlocal)
+{
+  return 4;
+}
diff --git a/src/fix_move.h b/src/fix_move.h
new file mode 100644
index 0000000000..86adbc2731
--- /dev/null
+++ b/src/fix_move.h
@@ -0,0 +1,62 @@
+/* ----------------------------------------------------------------------
+   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.
+------------------------------------------------------------------------- */
+
+#ifndef FIX_MOVE_H
+#define FIX_MOVE_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixMove : public Fix {
+ public:
+  FixMove(class LAMMPS *, int, char **);
+  ~FixMove();
+  int setmask();
+  void init();
+  void initial_integrate(int);
+  void final_integrate();
+  void initial_integrate_respa(int, int, int);
+  void final_integrate_respa(int);
+
+  double memory_usage();
+  void grow_arrays(int);
+  void copy_arrays(int, int);
+  int pack_exchange(int, double *);
+  int unpack_exchange(int, double *);
+  int pack_restart(int, double *);
+  void unpack_restart(int, int);
+  int maxsize_restart();
+  int size_restart(int);
+
+ private:
+  char *xvarstr,*yvarstr,*zvarstr,*vxvarstr,*vyvarstr,*vzvarstr;
+  int mstyle;
+  int vxflag,vyflag,vzflag,axflag,ayflag,azflag;
+  double vx,vy,vz,ax,ay,az;
+  double period,omega_rotate;
+  double point[3],axis[3],runit[3];
+  double dt,dtv,dtf;
+  int xvar,yvar,zvar,vxvar,vyvar,vzvar;
+  int xvarstyle,yvarstyle,zvarstyle,vxvarstyle,vyvarstyle,vzvarstyle;
+  int omega_flag,nlevels_respa;
+
+  double **xoriginal;         // original coords of atoms
+  int displaceflag,velocityflag;
+  int maxatom;
+  double **displace,**velocity;
+};
+
+}
+
+#endif
diff --git a/src/fix_wiggle.cpp b/src/fix_wiggle.cpp
deleted file mode 100644
index 4eaf97616e..0000000000
--- a/src/fix_wiggle.cpp
+++ /dev/null
@@ -1,186 +0,0 @@
-/* ----------------------------------------------------------------------
-   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 "string.h"
-#include "stdlib.h"
-#include "fix_wiggle.h"
-#include "atom.h"
-#include "update.h"
-#include "respa.h"
-#include "memory.h"
-#include "error.h"
-
-using namespace LAMMPS_NS;
-
-/* ---------------------------------------------------------------------- */
-
-FixWiggle::FixWiggle(LAMMPS *lmp, int narg, char **arg) :
-  Fix(lmp, narg, arg)
-{
-  if (narg != 6) error->all("Illegal fix wiggle command");
-
-  time_depend = 1;
-
-  // parse command-line args
-
-  if (strcmp(arg[3],"x") == 0) axis = 0;
-  else if (strcmp(arg[3],"y") == 0) axis = 1;
-  else if (strcmp(arg[3],"z") == 0) axis = 2;
-  else error->all("Illegal fix wiggle command");
-  
-  amplitude = atof(arg[4]);
-  period = atof(arg[5]);
-
-  // perform initial allocation of atom-based array
-  // register with Atom class
-  
-  original = NULL;
-  grow_arrays(atom->nmax);
-  atom->add_callback(0);
-
-  // oscillation parameters
-
-  double PI = 4.0 * atan(1.0);
-  omega = 2.0*PI / period;
-  time_origin = update->ntimestep;
-
-  // original = initial atom coord in wiggled direction
-  
-  double **x = atom->x;
-  int *mask = atom->mask;
-  int nlocal = atom->nlocal;
-
-  for (int i = 0; i < nlocal; i++)
-    if (mask[i] & groupbit) original[i] = x[i][axis];
-    else original[i] = 0.0;
-}
-
-/* ---------------------------------------------------------------------- */
-
-FixWiggle::~FixWiggle()
-{
-  // unregister callbacks to this fix from Atom class
-
-  atom->delete_callback(id,0);
-
-  // delete locally stored array
-
-  memory->sfree(original);
-}
-
-/* ---------------------------------------------------------------------- */
-
-int FixWiggle::setmask()
-{
-  int mask = 0;
-  mask |= POST_FORCE;
-  mask |= POST_FORCE_RESPA;
-  return mask;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixWiggle::init()
-{
-  dt = update->dt;
-
-  if (strcmp(update->integrate_style,"respa") == 0)
-    nlevels_respa = ((Respa *) update->integrate)->nlevels;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixWiggle::post_force(int vflag)
-{
-  double arg = omega * (update->ntimestep - time_origin) * dt;
-  double cosine = cos(arg);
-  double sine = sin(arg);
-
-  double **x = atom->x;
-  double **v = atom->v;
-  double **f = atom->f;
-  int *mask = atom->mask;
-  int nlocal = atom->nlocal;
-
-  for (int i = 0; i < nlocal; i++) {
-    if (mask[i] & groupbit) {
-      x[i][axis] = original[i] + amplitude - amplitude*cosine;
-      v[i][axis] = amplitude*omega*sine;
-      f[i][axis] = 0.0;
-    }
-  }
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixWiggle::post_force_respa(int vflag, int ilevel, int iloop)
-{
-  if (ilevel == nlevels_respa-1) post_force(vflag);
-}
-
-/* ----------------------------------------------------------------------
-   memory usage of local atom-based array 
-------------------------------------------------------------------------- */
-
-double FixWiggle::memory_usage()
-{
-  double bytes = atom->nmax * sizeof(double);
-  return bytes;
-}
-
-/* ----------------------------------------------------------------------
-   allocate local atom-based array 
-------------------------------------------------------------------------- */
-
-void FixWiggle::grow_arrays(int nmax)
-{
-  original = (double *)
-    memory->srealloc(original,nmax*sizeof(double),"wiggle:original");
-}
-
-/* ----------------------------------------------------------------------
-   copy values within local atom-based array 
-------------------------------------------------------------------------- */
-
-void FixWiggle::copy_arrays(int i, int j)
-{
-  original[j] = original[i];
-}
-
-/* ----------------------------------------------------------------------
-   pack values in local atom-based array for exchange with another proc 
-------------------------------------------------------------------------- */
-
-int FixWiggle::pack_exchange(int i, double *buf)
-{
-  buf[0] = original[i];
-  return 1;
-}
-
-/* ----------------------------------------------------------------------
-   unpack values in local atom-based array from exchange with another proc 
-------------------------------------------------------------------------- */
-
-int FixWiggle::unpack_exchange(int nlocal, double *buf)
-{
-  original[nlocal] = buf[0];
-  return 1;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixWiggle::reset_dt()
-{
-  dt = update->dt;
-}
diff --git a/src/fix_wiggle.h b/src/fix_wiggle.h
deleted file mode 100644
index 620d2e78f5..0000000000
--- a/src/fix_wiggle.h
+++ /dev/null
@@ -1,47 +0,0 @@
-/* ----------------------------------------------------------------------
-   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.
-------------------------------------------------------------------------- */
-
-#ifndef FIX_WIGGLE_H
-#define FIX_WIGGLE_H
-
-#include "fix.h"
-
-namespace LAMMPS_NS {
-
-class FixWiggle : public Fix {
- public:
-  FixWiggle(class LAMMPS *, int, char **);
-  ~FixWiggle();
-  int setmask();
-  void init();
-  void post_force(int);
-  void post_force_respa(int, int, int);
-
-  double memory_usage();
-  void grow_arrays(int);
-  void copy_arrays(int, int);
-  int pack_exchange(int, double *);
-  int unpack_exchange(int, double *);
-  void reset_dt();
-
- private:
-  double dt;
-  double *original;
-  double amplitude,period,omega;
-  int axis,time_origin;
-  int nlevels_respa;
-};
-
-}
-
-#endif
diff --git a/src/style.h b/src/style.h
index d628ed144f..01a676fb72 100644
--- a/src/style.h
+++ b/src/style.h
@@ -174,6 +174,7 @@ DumpStyle(xyz,DumpXYZ)
 #include "fix_minimize.h"
 #include "fix_msd.h"
 #include "fix_momentum.h"
+#include "fix_move.h"
 #include "fix_nph.h"
 #include "fix_npt.h"
 #include "fix_npt_sphere.h"
@@ -208,7 +209,6 @@ DumpStyle(xyz,DumpXYZ)
 #include "fix_wall_lj126.h"
 #include "fix_wall_lj93.h"
 #include "fix_wall_reflect.h"
-#include "fix_wiggle.h"
 #endif
 
 #ifdef FixClass
@@ -235,6 +235,7 @@ FixStyle(langevin,FixLangevin)
 FixStyle(lineforce,FixLineForce)
 FixStyle(MINIMIZE,FixMinimize)
 FixStyle(momentum,FixMomentum)
+FixStyle(move,FixMove)
 FixStyle(msd,FixMSD)
 FixStyle(nph,FixNPH)
 FixStyle(npt,FixNPT)
@@ -270,7 +271,6 @@ FixStyle(viscous,FixViscous)
 FixStyle(wall/lj126,FixWallLJ126)
 FixStyle(wall/lj93,FixWallLJ93)
 FixStyle(wall/reflect,FixWallReflect)
-FixStyle(wiggle,FixWiggle)
 #endif
 
 #ifdef ImproperInclude
diff --git a/src/variable.cpp b/src/variable.cpp
index 73761a66f8..1aa2bd4ba2 100644
--- a/src/variable.cpp
+++ b/src/variable.cpp
@@ -398,7 +398,8 @@ double Variable::compute_equal(int ivar)
 
 /* ----------------------------------------------------------------------
    compute result of atom-style variable evaluation
-   stride used since result may not be contiguous memory locs
+   only computed for atoms in igroup, else result is 0.0
+   answers are placed every stride locations into result
    if sumflag, add variable values to existing result
 ------------------------------------------------------------------------- */