git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@5023 f3b2605a-c512-4ea7-a41b-209d697bcdaa

src/COLLOID/fix_wall_colloid.cpp

@@ -84,11 +84,13 @@ void FixWallColloid::precompute(int m)
 }
 
 /* ----------------------------------------------------------------------
-   interaction of all particles in group with all 6 walls (if defined)
+   interaction of all particles in group with a wall
+   m = index of wall coeffs
+   which = xlo,xhi,ylo,yhi,zlo,zhi
    error if any finite-size particle is touching or penetrating wall
 ------------------------------------------------------------------------- */
 
-void FixWallColloid::wall_particle(int m, double coord)
+void FixWallColloid::wall_particle(int m, int which, double coord)
 {
   double delta,delta2,rinv,r2inv,r4inv,r8inv,fwall;
   double r2,rinv2,r2inv2,r4inv2,r6inv2;
@@ -102,8 +104,8 @@ void FixWallColloid::wall_particle(int m, double coord)
   int *type = atom->type;
   int nlocal = atom->nlocal;
 
-  int dim = m/2;
+  int dim = which  / 2;
-  int side = m % 2;
+  int side = which % 2;
   if (side == 0) side = -1;
 
   int onflag = 0;

src/COLLOID/fix_wall_colloid.h

@@ -29,7 +29,7 @@ class FixWallColloid : public FixWall {
   FixWallColloid(class LAMMPS *, int, char **);
   void init();
   void precompute(int);
-  void wall_particle(int, double);
+  void wall_particle(int, int, double);
 
  private:
   double coeff1[6],coeff2[6],coeff3[6],coeff4[6],offset[6];

src/fix_gravity.cpp

@@ -62,8 +62,8 @@ FixGravity::FixGravity(LAMMPS *lmp, int narg, char **arg) :
     zdir = atof(arg[7]);
   } else error->all("Illegal fix gravity command");
 
-  double PI = 2.0 * asin(1.0);
+  double PI = 4.0*atan(1.0);
-  degree2rad = 2.0*PI / 360.0;
+  degree2rad = PI/180.0;
 
   if (style == CHUTE || style == SPHERICAL || style == GRADIENT) {
     if (domain->dimension == 3) {

src/fix_wall.cpp

@@ -16,16 +16,18 @@
 #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 "output.h"
 #include "respa.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
 
 enum{XLO,XHI,YLO,YHI,ZLO,ZHI};
+enum{NONE,EDGE,CONSTANT,VARIABLE};
 
 /* ---------------------------------------------------------------------- */
 
@@ -34,17 +36,13 @@ FixWall::FixWall(LAMMPS *lmp, int narg, char **arg) :
 {
   scalar_flag = 1;
   vector_flag = 1;
-  size_vector = 6;
   global_freq = 1;
   extscalar = 1;
   extvector = 1;
-  time_depend = 1;
 
   // parse args
 
-  for (int m = 0; m < 6; m++) wallflag[m] = 0;
+  nwall = 0;
-  velflag = 0;
-  wigflag = 0;
   int scaleflag = 1;
 
   int iarg = 3;
@@ -53,39 +51,42 @@ FixWall::FixWall(LAMMPS *lmp, int narg, char **arg) :
 	(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("Illegal fix wall command");
-      int m;
+
-      if (strcmp(arg[iarg],"xlo") == 0) m = XLO;
+      int newwall;
-      else if (strcmp(arg[iarg],"xhi") == 0) m = XHI;
+      if (strcmp(arg[iarg],"xlo") == 0) newwall = XLO;
-      else if (strcmp(arg[iarg],"ylo") == 0) m = YLO;
+      else if (strcmp(arg[iarg],"xhi") == 0) newwall = XHI;
-      else if (strcmp(arg[iarg],"yhi") == 0) m = YHI;
+      else if (strcmp(arg[iarg],"ylo") == 0) newwall = YLO;
-      else if (strcmp(arg[iarg],"zlo") == 0) m = ZLO;
+      else if (strcmp(arg[iarg],"yhi") == 0) newwall = YHI;
-      else if (strcmp(arg[iarg],"zhi") == 0) m = ZHI;
+      else if (strcmp(arg[iarg],"zlo") == 0) newwall = ZLO;
-      wallflag[m] = 1;
+      else if (strcmp(arg[iarg],"zhi") == 0) newwall = ZHI;
-      coord0[m] = atof(arg[iarg+1]);
+
-      epsilon[m] = atof(arg[iarg+2]);
+      for (int m = 0; m < nwall; m++)
-      sigma[m] = atof(arg[iarg+3]);
+	if (newwall == wallwhich[m])
-      cutoff[m] = atof(arg[iarg+4]);
+	  error->all("Wall defined twice in fix wall 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] = atof(arg[iarg+1]);
+      }
+
+      epsilon[nwall] = atof(arg[iarg+2]);
+      sigma[nwall] = atof(arg[iarg+3]);
+      cutoff[nwall] = atof(arg[iarg+4]);
+      nwall++;
       iarg += 5;
-    } else if (strcmp(arg[iarg],"vel") == 0) {
+
-      if (iarg+2 > narg) error->all("Illegal fix wall command");
-      velflag = 1;
-      double vtmp = atof(arg[iarg+1]);
-      for (int m = 0; m < 6; m++) vel[m] = vtmp;
-      iarg += 2;
-    } else if (strcmp(arg[iarg],"wiggle/sin") == 0) {
-      if (iarg+3 > narg) error->all("Illegal fix wall command");
-      wigflag = 1;
-      double atmp = atof(arg[iarg+1]);
-      period = atof(arg[iarg+2]);
-      for (int m = 0; m < 6; m++) amplitude[m] = atmp;
-      iarg += 3;
-    } else if (strcmp(arg[iarg],"wiggle/cos") == 0) {
-      if (iarg+3 > narg) error->all("Illegal fix wall command");
-      wigflag = 2;
-      double atmp = atof(arg[iarg+1]);
-      period = atof(arg[iarg+2]);
-      for (int m = 0; m < 6; m++) amplitude[m] = atmp;
-      iarg += 3;
     } else if (strcmp(arg[iarg],"units") == 0) {
       if (iarg+2 > narg) error->all("Illegal fix wall command");
       if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0;
@@ -95,66 +96,70 @@ FixWall::FixWall(LAMMPS *lmp, int narg, char **arg) :
     } else error->all("Illegal fix wall command");
   }
 
+  size_vector = nwall;
+
   // error check
 
-  int flag = 0;
+  if (nwall == 0) error->all("Illegal fix wall command");
-  for (int m = 0; m < 6; m++) if (wallflag[m]) flag = 1;
+  for (int m = 0; m < nwall; m++)
-  if (!flag) error->all("Illegal fix wall command");
+    if (cutoff[m] <= 0.0)
-
-  for (int m = 0; m < 6; m++)
-    if (wallflag[m] && cutoff[m] <= 0.0)
       error->all("Fix wall cutoff <= 0.0");
 
-  if (velflag && wigflag) 
+  for (int m = 0; m < nwall; m++) {
-    error->all("Cannot set both vel and wiggle in fix wall command");
+    if ((wallwhich[m] == XLO || wallwhich[m] == XHI) && domain->xperiodic)
-
+      error->all("Cannot use fix wall in periodic dimension");
-  if ((wallflag[XLO] || wallflag[XHI]) && domain->xperiodic)
+    if ((wallwhich[m] == YLO || wallwhich[m] == YHI) && domain->yperiodic)
-    error->all("Cannot use fix wall in periodic dimension");
+      error->all("Cannot use fix wall in periodic dimension");
-  if ((wallflag[YLO] || wallflag[YHI]) && domain->yperiodic)
+    if ((wallwhich[m] == ZLO || wallwhich[m] == ZHI) && domain->zperiodic)
-    error->all("Cannot use fix wall in periodic dimension");
+      error->all("Cannot use fix wall in periodic dimension");
-  if ((wallflag[ZLO] || wallflag[ZHI]) && domain->zperiodic)
-    error->all("Cannot use fix wall in periodic dimension");
-
-  if ((wallflag[ZLO] || wallflag[ZHI]) && domain->dimension == 2)
-    error->all("Cannot use fix wall zlo/zhi for a 2d simulation");
-
-  // setup scaling
-
-  if (scaleflag && domain->lattice == NULL)
-    error->all("Use of fix wall 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;
-
-  // apply scaling factors to coord0, vel, amplitude
-
-  double scale;
-  for (int m = 0; m < 6; m++) {
-    if (wallflag[m] == 0) continue;
-    if (m < 2) scale = xscale;
-    else if (m < 4) scale = yscale;
-    else scale = zscale;
-    coord0[m] *= scale;
-    if (velflag) vel[m] *= scale;
-    if (wigflag) amplitude[m] *= scale;
   }
 
-  // setup oscillations
+  for (int m = 0; m < nwall; m++)
+    if ((wallwhich[m] == ZLO || wallwhich[m] == ZHI) && domain->dimension == 2)
+      error->all("Cannot use fix wall zlo/zhi for a 2d simulation");
   
-  if (wigflag) {
+  // scale coord for CONSTANT walls
-    double PI = 4.0 * atan(1.0);
+
-    omega = 2.0*PI / period;
+  int flag = 0;
+  for (int m = 0; m < nwall; m++)
+    if (wallstyle[m] == CONSTANT) flag = 1;
+
+  if (flag) {
+    if (scaleflag && domain->lattice == NULL)
+      error->all("Use of fix wall 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;
+
+    double scale;
+    for (int m = 0; m < nwall; m++) {
+      if (wallwhich[m] < YLO) scale = xscale;
+      else if (wallwhich[m] < ZLO) scale = yscale;
+      else scale = zscale;
+      if (wallstyle[m] == CONSTANT) coord0[m] *= scale;
+    }
   }
 
+  // set time_depend if any wall positions are variable
+
+  for (int m = 0; m < nwall; m++)
+    if (wallstyle[m] == VARIABLE) time_depend = 1;
+
   eflag = 0;
-  for (int m = 0; m < 7; m++) ewall[m] = 0.0;
+  for (int m = 0; m <= nwall; m++) ewall[m] = 0.0;
+}
 
-  time_origin = update->ntimestep;
+/* ---------------------------------------------------------------------- */
+
+FixWall::~FixWall()
+{
+  for (int m = 0; m < nwall; m++)
+    if (wallstyle[m] == VARIABLE) delete [] varstr[m];
 }
 
 /* ---------------------------------------------------------------------- */
@@ -175,10 +180,18 @@ void FixWall::init()
 {
   dt = update->dt;
 
+  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("Variable name for fix wall does not exist");
+    if (!input->variable->equalstyle(varindex[m]))
+      error->all("Variable for fix wall is invalid style");
+  }
+
   // setup coefficients
 
-  for (int m = 0; m < 6; m++)
+  for (int m = 0; m < nwall; m++) precompute(m);
-    if (wallflag[m]) precompute(m);
 
   if (strcmp(update->integrate_style,"respa") == 0)
     nlevels_respa = ((Respa *) update->integrate)->nlevels;
@@ -209,29 +222,18 @@ void FixWall::min_setup(int vflag)
 void FixWall::post_force(int vflag)
 {
   eflag = 0;
-  for (int m = 0; m < 7; m++) ewall[m] = 0.0;
+  for (int m = 0; m <= nwall; m++) ewall[m] = 0.0;
 
-  // coord0 = initial position of wall
   // coord = current position of wall
+  // evaluate variable if necessary
 
-  double delta = (update->ntimestep - time_origin) * dt;
   double coord;
+  for (int m = 0; m < nwall; m++) {
+    if (wallstyle[m] == VARIABLE)
+      coord = input->variable->compute_equal(varindex[m]);
+    else coord = coord0[m];
 
-  for (int m = 0; m < 6; m++) {
+    wall_particle(m,wallwhich[m],coord);
-    if (wallflag[m] == 0) continue;
-
-    if (velflag) {
-      if (m % 2 == 0) coord = coord0[m] + delta*vel[m];
-      else coord = coord0[m] - delta*vel[m];
-    } else if (wigflag == 1) {
-      if (m % 2 == 0) coord = coord0[m] + amplitude[m]*sin(omega*delta);
-      else coord = coord0[m] - amplitude[m]*sin(omega*delta);
-    } else if (wigflag == 2) {
-      if (m % 2 == 0) coord = coord0[m] + amplitude[m]*(1.0-cos(omega*delta));
-      else coord = coord0[m] - amplitude[m]*(1.0-cos(omega*delta));
-    } else coord = coord0[m];
-
-    wall_particle(m,coord);
   }
 }
 
@@ -258,7 +260,7 @@ double FixWall::compute_scalar()
   // only sum across procs one time
 
   if (eflag == 0) {
-    MPI_Allreduce(ewall,ewall_all,7,MPI_DOUBLE,MPI_SUM,world);
+    MPI_Allreduce(ewall,ewall_all,nwall+1,MPI_DOUBLE,MPI_SUM,world);
     eflag = 1;
   }
   return ewall_all[0];
@@ -273,7 +275,7 @@ double FixWall::compute_vector(int n)
   // only sum across procs one time
 
   if (eflag == 0) {
-    MPI_Allreduce(ewall,ewall_all,7,MPI_DOUBLE,MPI_SUM,world);
+    MPI_Allreduce(ewall,ewall_all,nwall+1,MPI_DOUBLE,MPI_SUM,world);
     eflag = 1;
   }
   return ewall_all[n+1];

src/fix_wall.h

@@ -21,7 +21,7 @@ namespace LAMMPS_NS {
 class FixWall : public Fix {
  public:
   FixWall(class LAMMPS *, int, char **);
-  virtual ~FixWall() {}
+  virtual ~FixWall();
   int setmask();
   virtual void init();
   void setup(int);
@@ -33,14 +33,14 @@ class FixWall : public Fix {
   double compute_vector(int);
 
   virtual void precompute(int) = 0;
-  virtual void wall_particle(int, double) = 0;
+  virtual void wall_particle(int, int, double) = 0;
 
  protected:
-  int wallflag[6];
+  int nwall;
+  int wallwhich[6],wallstyle[6];
   double coord0[6],epsilon[6],sigma[6],cutoff[6];
-  int velflag,wigflag;
+  char *varstr[6];
-  double vel[6],amplitude[6];
+  int varindex[6];
-  double period,omega;
   int eflag;
   double ewall[7],ewall_all[7];
   int nlevels_respa;

src/fix_wall_harmonic.cpp

@@ -24,11 +24,13 @@ FixWallHarmonic::FixWallHarmonic(LAMMPS *lmp, int narg, char **arg) :
   FixWall(lmp, narg, arg) {}
 
 /* ----------------------------------------------------------------------
-   interaction of all particles in group with all 6 walls (if defined)
+   interaction of all particles in group with a wall
+   m = index of wall coeffs
+   which = xlo,xhi,ylo,yhi,zlo,zhi
    error if any particle is on or behind wall
 ------------------------------------------------------------------------- */
 
-void FixWallHarmonic::wall_particle(int m, double coord)
+void FixWallHarmonic::wall_particle(int m, int which, double coord)
 {
   double delta,dr,fwall;
 
@@ -37,8 +39,8 @@ void FixWallHarmonic::wall_particle(int m, double coord)
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
-  int dim = m/2;
+  int dim = which / 2;
-  int side = m % 2;
+  int side = which % 2;
   if (side == 0) side = -1;
 
   int onflag = 0;

src/fix_wall_harmonic.h

@@ -28,7 +28,7 @@ class FixWallHarmonic : public FixWall {
  public:
   FixWallHarmonic(class LAMMPS *, int, char **);
   void precompute(int) {}
-  void wall_particle(int, double);
+  void wall_particle(int, int, double);
 
  private:
   double offset[6];

src/fix_wall_lj126.cpp

@@ -38,11 +38,13 @@ void FixWallLJ126::precompute(int m)
 }
 
 /* ----------------------------------------------------------------------
-   interaction of all particles in group with all 6 walls (if defined)
+   interaction of all particles in group with a wall
+   m = index of wall coeffs
+   which = xlo,xhi,ylo,yhi,zlo,zhi
    error if any particle is on or behind wall
 ------------------------------------------------------------------------- */
 
-void FixWallLJ126::wall_particle(int m, double coord)
+void FixWallLJ126::wall_particle(int m, int which, double coord)
 {
   double delta,rinv,r2inv,r6inv,fwall;
 
@@ -51,8 +53,8 @@ void FixWallLJ126::wall_particle(int m, double coord)
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
-  int dim = m/2;
+  int dim = which / 2;
-  int side = m % 2;
+  int side = which % 2;
   if (side == 0) side = -1;
 
   int onflag = 0;

src/fix_wall_lj126.h

@@ -28,7 +28,7 @@ class FixWallLJ126 : public FixWall {
  public:
   FixWallLJ126(class LAMMPS *, int, char **);
   void precompute(int);
-  void wall_particle(int, double);
+  void wall_particle(int, int, double);
 
  private:
   double coeff1[6],coeff2[6],coeff3[6],coeff4[6],offset[6];

src/fix_wall_lj93.cpp

@@ -39,11 +39,13 @@ void FixWallLJ93::precompute(int m)
 }
 
 /* ----------------------------------------------------------------------
-   interaction of all particles in group with all 6 walls (if defined)
+   interaction of all particles in group with a wall
+   m = index of wall coeffs
+   which = xlo,xhi,ylo,yhi,zlo,zhi
    error if any particle is on or behind wall
 ------------------------------------------------------------------------- */
 
-void FixWallLJ93::wall_particle(int m, double coord)
+void FixWallLJ93::wall_particle(int m, int which, double coord)
 {
   double delta,rinv,r2inv,r4inv,r10inv,fwall;
 
@@ -52,8 +54,8 @@ void FixWallLJ93::wall_particle(int m, double coord)
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
-  int dim = m/2;
+  int dim = which / 2;
-  int side = m % 2;
+  int side = which % 2;
   if (side == 0) side = -1;
 
   int onflag = 0;

src/fix_wall_lj93.h

@@ -28,7 +28,7 @@ class FixWallLJ93 : public FixWall {
  public:
   FixWallLJ93(class LAMMPS *, int, char **);
   void precompute(int);
-  void wall_particle(int, double);
+  void wall_particle(int, int, double);
 
  private:
   double coeff1[6],coeff2[6],coeff3[6],coeff4[6],offset[6];

src/fix_wall_reflect.cpp

@@ -25,6 +25,7 @@
 
 using namespace LAMMPS_NS;
 
+enum{XLO,XHI,YLO,YHI,ZLO,ZHI};
 enum{NONE,EDGE,CONSTANT,VARIABLE};
 
 /* ---------------------------------------------------------------------- */
@@ -34,101 +35,48 @@ FixWallReflect::FixWallReflect(LAMMPS *lmp, int narg, char **arg) :
 {
   if (narg < 4) error->all("Illegal fix wall/reflect command");
 
-  xloflag = xhiflag = yloflag = yhiflag = zloflag = zhiflag = NONE;
+  // parse args
-  xlostr = xhistr = ylostr = yhistr = zlostr = zhistr = NULL;
+
+  nwall = 0;
   int scaleflag = 1;
 
   int iarg = 3;
   while (iarg < narg) {
-    if (strcmp(arg[iarg],"xlo") == 0) {
+    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("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++)
+	if (newwall == wallwhich[m])
+	  error->all("Wall defined twice in fix wall/reflect command");
+
+      wallwhich[nwall] = newwall;
       if (strcmp(arg[iarg+1],"EDGE") == 0) {
-	xloflag = EDGE;
+	wallstyle[nwall] = EDGE;
-	xlo = domain->boxlo[0];
+	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]) {
-	xloflag = VARIABLE;
+	wallstyle[nwall] = VARIABLE;
 	int n = strlen(&arg[iarg+1][2]) + 1;
-	xlostr = new char[n];
+	varstr[nwall] = new char[n];
-	strcpy(xlostr,&arg[iarg+1][2]);
+	strcpy(varstr[nwall],&arg[iarg+1][2]);
       } else {
-	xloflag = CONSTANT;
+	wallstyle[nwall] = CONSTANT;
-	xlo = atof(arg[iarg+1]);
+	coord0[nwall] = atof(arg[iarg+1]);
-      }
-      iarg += 2;
-    } else if (strcmp(arg[iarg],"xhi") == 0) {
-      if (iarg+2 > narg) error->all("Illegal fix wall/reflect command");
-      if (strcmp(arg[iarg+1],"EDGE") == 0) {
-	xhiflag = EDGE;
-	xhi = domain->boxhi[0];
-      } else if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
-	xhiflag = VARIABLE;
-	int n = strlen(&arg[iarg+1][2]) + 1;
-	xhistr = new char[n];
-	strcpy(xhistr,&arg[iarg+1][2]);
-      } else {
-	xhiflag = CONSTANT;
-	xhi = atof(arg[iarg+1]);
-      }
-      iarg += 2;
-    } else if (strcmp(arg[iarg],"ylo") == 0) {
-      if (iarg+2 > narg) error->all("Illegal fix wall/reflect command");
-      if (strcmp(arg[iarg+1],"EDGE") == 0) {
-	yloflag = EDGE;
-	ylo = domain->boxlo[1];
-      } else if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
-	yloflag = VARIABLE;
-	int n = strlen(&arg[iarg+1][2]) + 1;
-	ylostr = new char[n];
-	strcpy(ylostr,&arg[iarg+1][2]);
-      } else {
-	yloflag = CONSTANT;
-	ylo = atof(arg[iarg+1]);
-      }
-      iarg += 2;
-    } else if (strcmp(arg[iarg],"yhi") == 0) {
-      if (iarg+2 > narg) error->all("Illegal fix wall/reflect command");
-      if (strcmp(arg[iarg+1],"EDGE") == 0) {
-	yhiflag = EDGE;
-	yhi = domain->boxhi[1];
-      } else if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
-	yhiflag = VARIABLE;
-	int n = strlen(&arg[iarg+1][2]) + 1;
-	yhistr = new char[n];
-	strcpy(yhistr,&arg[iarg+1][2]);
-      } else {
-	yhiflag = CONSTANT;
-	yhi = atof(arg[iarg+1]);
-      }
-      iarg += 2;
-    } else if (strcmp(arg[iarg],"zlo") == 0) {
-      if (iarg+2 > narg) error->all("Illegal fix wall/reflect command");
-      if (strcmp(arg[iarg+1],"EDGE") == 0) {
-	zloflag = EDGE;
-	zlo = domain->boxlo[2];
-      } else if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
-	zloflag = VARIABLE;
-	int n = strlen(&arg[iarg+1][2]) + 1;
-	zlostr = new char[n];
-	strcpy(zlostr,&arg[iarg+1][2]);
-      } else {
-	zloflag = CONSTANT;
-	zlo = atof(arg[iarg+1]);
-      }
-      iarg += 2;
-    } else if (strcmp(arg[iarg],"zhi") == 0) {
-      if (iarg+2 > narg) error->all("Illegal fix wall/reflect command");
-      if (strcmp(arg[iarg+1],"EDGE") == 0) {
-	zhiflag = EDGE;
-	zhi = domain->boxhi[2];
-      } else if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
-	zhiflag = VARIABLE;
-	int n = strlen(&arg[iarg+1][2]) + 1;
-	zhistr = new char[n];
-	strcpy(zhistr,&arg[iarg+1][2]);
-      } else {
-	zhiflag = CONSTANT;
-	zhi = atof(arg[iarg+1]);
       }
+
+      nwall++;
       iarg += 2;
 
     } else if (strcmp(arg[iarg],"units") == 0) {
@@ -140,15 +88,32 @@ FixWallReflect::FixWallReflect(LAMMPS *lmp, int narg, char **arg) :
     } else error->all("Illegal fix wall/reflect command");
   }
 
-  // setup scaling if any wall positions were set as CONSTANT
+  // error check
-  // apply scaling factors to wall positions
 
-  if (xloflag == CONSTANT || xhiflag == CONSTANT ||
+  if (nwall == 0) error->all("Illegal fix wall command");
-      yloflag == CONSTANT || yhiflag == CONSTANT ||
-      zloflag == CONSTANT || zhiflag == CONSTANT) {
 
+  for (int m = 0; m < nwall; m++) {
+    if ((wallwhich[m] == XLO || wallwhich[m] == XHI) && domain->xperiodic)
+      error->all("Cannot use fix wall/reflect in periodic dimension");
+    if ((wallwhich[m] == YLO || wallwhich[m] == YHI) && domain->yperiodic)
+      error->all("Cannot use fix wall/reflect in periodic dimension");
+    if ((wallwhich[m] == ZLO || wallwhich[m] == ZHI) && domain->zperiodic)
+      error->all("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("Cannot use fix wall/reflect zlo/zhi for a 2d simulation");
+
+  // scale coord for CONSTANT walls
+
+  int flag = 0;
+  for (int m = 0; m < nwall; m++)
+    if (wallstyle[m] == CONSTANT) flag = 1;
+
+  if (flag) {
     if (scaleflag && domain->lattice == NULL)
-      error->all("Use of fix wall/reflect with undefined lattice");
+      error->all("Use of fix wall with undefined lattice");
 
     double xscale,yscale,zscale;
     if (scaleflag) {
@@ -158,37 +123,27 @@ FixWallReflect::FixWallReflect(LAMMPS *lmp, int narg, char **arg) :
     }
     else xscale = yscale = zscale = 1.0;
 
-    if (xloflag == CONSTANT) xlo *= xscale;
+    double scale;
-    if (xhiflag == CONSTANT) xhi *= xscale;
+    for (int m = 0; m < nwall; m++) {
-    if (yloflag == CONSTANT) ylo *= xscale;
+      if (wallwhich[m] < YLO) scale = xscale;
-    if (yhiflag == CONSTANT) yhi *= xscale;
+      else if (wallwhich[m] < ZLO) scale = yscale;
-    if (zloflag == CONSTANT) zlo *= xscale;
+      else scale = zscale;
-    if (zhiflag == CONSTANT) zhi *= xscale;
+      if (wallstyle[m] == CONSTANT) coord0[m] *= scale;
+    }
   }
 
-  // error checks
+  // set time_depend if any wall positions are variable
 
-  if (!xloflag && !xhiflag && !yloflag && !yhiflag && !zloflag && !zhiflag)
+  for (int m = 0; m < nwall; m++)
-    error->all("Illegal fix wall/reflect command");
+    if (wallstyle[m] == VARIABLE) time_depend = 1;
-
-  if ((xloflag || xhiflag) && domain->xperiodic)
-    error->all("Cannot use wall in periodic dimension");
-  if ((yloflag || yhiflag) && domain->yperiodic)
-    error->all("Cannot use wall in periodic dimension");
-  if ((zloflag || zhiflag) && domain->zperiodic)
-    error->all("Cannot use wall in periodic dimension");
 }
 
 /* ---------------------------------------------------------------------- */
 
 FixWallReflect::~FixWallReflect()
 {
-  delete [] xlostr;
+  for (int m = 0; m < nwall; m++)
-  delete [] xhistr;
+    if (wallstyle[m] == VARIABLE) delete [] varstr[m];
-  delete [] ylostr;
-  delete [] yhistr;
-  delete [] zlostr;
-  delete [] zhistr;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -205,46 +160,12 @@ int FixWallReflect::setmask()
 
 void FixWallReflect::init()
 {
-  if (xlostr) {
+  for (int m = 0; m < nwall; m++) {
-    xlovar = input->variable->find(xlostr);
+    if (wallstyle[m] != VARIABLE) continue;
-    if (xlovar < 0)
+    varindex[m] = input->variable->find(varstr[m]);
-      error->all("Variable name for fix wall/relect does not exist");
+    if (varindex[m] < 0)
-    if (!input->variable->equalstyle(xlovar))
+      error->all("Variable name for fix wall/reflect does not exist");
-      error->all("Variable for fix wall/reflect is invalid style");
+    if (!input->variable->equalstyle(varindex[m]))
-  }
-  if (xhistr) {
-    xhivar = input->variable->find(xhistr);
-    if (xhivar < 0)
-      error->all("Variable name for fix wall/relect does not exist");
-    if (!input->variable->equalstyle(xhivar))
-      error->all("Variable for fix wall/reflect is invalid style");
-  }
-  if (ylostr) {
-    ylovar = input->variable->find(ylostr);
-    if (ylovar < 0)
-      error->all("Variable name for fix wall/relect does not exist");
-    if (!input->variable->equalstyle(ylovar))
-      error->all("Variable for fix wall/reflect is invalid style");
-  }
-  if (yhistr) {
-    yhivar = input->variable->find(yhistr);
-    if (yhivar < 0)
-      error->all("Variable name for fix wall/relect does not exist");
-    if (!input->variable->equalstyle(yhivar))
-      error->all("Variable for fix wall/reflect is invalid style");
-  }
-  if (zlostr) {
-    zlovar = input->variable->find(zlostr);
-    if (zlovar < 0)
-      error->all("Variable name for fix wall/relect does not exist");
-    if (!input->variable->equalstyle(zlovar))
-      error->all("Variable for fix wall/reflect is invalid style");
-  }
-  if (zhistr) {
-    zhivar = input->variable->find(zhistr);
-    if (zhivar < 0)
-      error->all("Variable name for fix wall/relect does not exist");
-    if (!input->variable->equalstyle(zhivar))
       error->all("Variable for fix wall/reflect is invalid style");
   }
 
@@ -261,44 +182,35 @@ void FixWallReflect::init()
 
 void FixWallReflect::post_integrate()
 {
-  if (xloflag == VARIABLE) xlo = input->variable->compute_equal(xlovar);
+  int i,dim,side;
-  if (xhiflag == VARIABLE) xhi = input->variable->compute_equal(xhivar);
+  double coord;
-  if (yloflag == VARIABLE) ylo = input->variable->compute_equal(ylovar);
-  if (yhiflag == VARIABLE) yhi = input->variable->compute_equal(yhivar);
-  if (zloflag == VARIABLE) zlo = input->variable->compute_equal(zlovar);
-  if (zhiflag == VARIABLE) zhi = input->variable->compute_equal(zhivar);
 
   double **x = atom->x;
   double **v = atom->v;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
-  for (int i = 0; i < nlocal; i++) {
+  for (int m = 0; m < nwall; m++) {
-    if (mask[i] & groupbit) {
+    if (wallstyle[m] == VARIABLE)
-      if (xloflag && x[i][0] < xlo) {
+      coord = input->variable->compute_equal(varindex[m]);
-	x[i][0] = xlo + (xlo - x[i][0]);
+    else coord = coord0[m];
-	v[i][0] = -v[i][0];
+
+    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][0] > coord) {
+	    x[i][dim] = coord - (x[i][dim] - coord);
+	    v[i][dim] = -v[i][dim];
+	  }
+	}
       }
-      if (xhiflag && x[i][0] > xhi) {
-	x[i][0] = xhi - (x[i][0] - xhi);
-	v[i][0] = -v[i][0];
-      }
-      if (yloflag && x[i][1] < ylo) {
-	x[i][1] = ylo + (ylo - x[i][1]);
-	v[i][1] = -v[i][1];
-      }
-      if (yhiflag && x[i][1] > yhi) {
-	x[i][1] = yhi - (x[i][1] - yhi);
-	v[i][1] = -v[i][1];
-      }
-      if (zloflag && x[i][2] < zlo) {
-	x[i][2] = zlo + (zlo - x[i][2]);
-	v[i][2] = -v[i][2];
-      }
-      if (zhiflag && x[i][2] > zhi) {
-	x[i][2] = zhi - (x[i][2] - zhi);
-	v[i][2] = -v[i][2];
-      }
-    }
   }
 }

src/fix_wall_reflect.h

@@ -33,10 +33,11 @@ class FixWallReflect : public Fix {
   void post_integrate();
 
  private:
-  int xloflag,xhiflag,yloflag,yhiflag,zloflag,zhiflag;
+  int nwall;
-  char *xlostr,*xhistr,*ylostr,*yhistr,*zlostr,*zhistr;
+  int wallwhich[6],wallstyle[6];
-  int xlovar,xhivar,ylovar,yhivar,zlovar,zhivar;
+  double coord0[6];
-  double xlo,xhi,ylo,yhi,zlo,zhi;
+  char *varstr[6];
+  int varindex[6];
 };
 
 }

src/variable.cpp

@@ -43,11 +43,18 @@ using namespace LAMMPS_NS;
 
 enum{INDEX,LOOP,WORLD,UNIVERSE,ULOOP,STRING,EQUAL,ATOM};
 enum{ARG,OP};
+
+// customize by adding a math function
+
 enum{DONE,ADD,SUBTRACT,MULTIPLY,DIVIDE,CARAT,UNARY,
        EQ,NE,LT,LE,GT,GE,AND,OR,
        SQRT,EXP,LN,LOG,SIN,COS,TAN,ASIN,ACOS,ATAN,ATAN2,
-       RANDOM,NORMAL,CEIL,FLOOR,ROUND,
+       RANDOM,NORMAL,CEIL,FLOOR,ROUND,RAMP,STAGGER,LOGFREQ,
+       VLINEAR,SWIGGLE,CWIGGLE,
        VALUE,ATOMARRAY,TYPEARRAY,INTARRAY};
+
+// customize by adding a special function
+
 enum{SUM,XMIN,XMAX,AVE,TRAP};
 
 #define INVOKED_SCALAR 1
@@ -83,6 +90,8 @@ Variable::Variable(LAMMPS *lmp) : Pointers(lmp)
   precedence[MULTIPLY] = precedence[DIVIDE] = 6;
   precedence[CARAT] = 7;
   precedence[UNARY] = 8;
+
+  PI = 4.0*atan(1.0);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -595,6 +604,7 @@ void Variable::copy(int narg, char **from, char **to)
                       sqrt(x),exp(x),ln(x),log(x),
 		      sin(x),cos(x),tan(x),asin(x),atan2(y,x),...
      group function = count(group), mass(group), xcm(group,x), ...
+     special function = sum(x),min(x), ...
      atom value = x[i], y[i], vx[i], ...
      atom vector = x, y, vx, ...
      compute = c_ID, c_ID[i], c_ID[i][j]
@@ -686,10 +696,10 @@ double Variable::evaluate(char *str, Tree **tree)
 
     // ----------------
     // letter: c_ID, c_ID[], c_ID[][], f_ID, f_ID[], f_ID[][],
-    //         v_name, v_name[], exp(), xcm(,), x, x[], vol
+    //         v_name, v_name[], exp(), xcm(,), x, x[], PI, vol
     // ----------------
 
-    } else if (islower(onechar)) {
+    } else if (isalpha(onechar)) {
       if (expect == OP) error->all("Invalid syntax in variable formula");
       expect = OP;
 
@@ -1134,7 +1144,8 @@ double Variable::evaluate(char *str, Tree **tree)
 	delete [] id;
 
       // ----------------
-      // math/group/special function or atom value/vector or thermo keyword
+      // math/group/special function or atom value/vector or
+      // constant or thermo keyword
       // ----------------
 
       } else {
@@ -1183,6 +1194,20 @@ double Variable::evaluate(char *str, Tree **tree)
 
 	  atom_vector(word,tree,treestack,ntreestack);
 
+	// ----------------
+	// constant
+	// ----------------
+
+	} else if (is_constant(word)) {
+	  value1 = constant(word);
+	  if (tree) {
+	    Tree *newtree = new Tree();
+	    newtree->type = VALUE;
+	    newtree->value = value1;
+	    newtree->left = newtree->middle = newtree->right = NULL;
+	    treestack[ntreestack++] = newtree;
+	  } else argstack[nargstack++] = value1;
+
 	// ----------------
 	// thermo keyword
 	// ----------------
@@ -1350,7 +1375,9 @@ double Variable::evaluate(char *str, Tree **tree)
    process an evaulation tree
    customize by adding a math function:
      sqrt(),exp(),ln(),log(),sin(),cos(),tan(),asin(),acos(),atan(),
-     atan2(y,x),random(x,y,z),normal(x,y,z),ceil(),floor(),round()
+     atan2(y,x),random(x,y,z),normal(x,y,z),ceil(),floor(),round(),
+     ramp(x,y),stagger(x,y),logfreq(x,y,z),
+     vlinear(x,y),swiggle(x,y,z),cwiggle(x,y,z)
 ---------------------------------------------------------------------- */
 
 double Variable::eval_tree(Tree *tree, int i)
@@ -1368,12 +1395,12 @@ double Variable::eval_tree(Tree *tree, int i)
     return eval_tree(tree->left,i) * eval_tree(tree->right,i);
   if (tree->type == DIVIDE) {
     double denom = eval_tree(tree->right,i);
-    if (denom == 0.0) error->all("Divide by 0 in variable formula");
+    if (denom == 0.0) error->one("Divide by 0 in variable formula");
     return eval_tree(tree->left,i) / denom;
   }
   if (tree->type == CARAT) {
     double exponent = eval_tree(tree->right,i);
-    if (exponent == 0.0) error->all("Power by 0 in variable formula");
+    if (exponent == 0.0) error->one("Power by 0 in variable formula");
     return pow(eval_tree(tree->left,i),exponent);
   }
   if (tree->type == UNARY)
@@ -1416,7 +1443,7 @@ double Variable::eval_tree(Tree *tree, int i)
 
   if (tree->type == SQRT) {
     double arg = eval_tree(tree->left,i);
-    if (arg < 0.0) error->all("Sqrt of negative value in variable formula");
+    if (arg < 0.0) error->one("Sqrt of negative value in variable formula");
     return sqrt(arg);
   }
   if (tree->type == EXP)
@@ -1424,13 +1451,13 @@ double Variable::eval_tree(Tree *tree, int i)
   if (tree->type == LN) {
     double arg = eval_tree(tree->left,i);
     if (arg <= 0.0) 
-      error->all("Log of zero/negative value in variable formula");
+      error->one("Log of zero/negative value in variable formula");
     return log(arg);
   }
   if (tree->type == LOG) {
     double arg = eval_tree(tree->left,i);
     if (arg <= 0.0) 
-      error->all("Log of zero/negative value in variable formula");
+      error->one("Log of zero/negative value in variable formula");
     return log10(arg);
   }
 
@@ -1444,13 +1471,13 @@ double Variable::eval_tree(Tree *tree, int i)
   if (tree->type == ASIN) {
     double arg = eval_tree(tree->left,i);
     if (arg < -1.0 || arg > 1.0)
-      error->all("Arcsin of invalid value in variable formula");
+      error->one("Arcsin of invalid value in variable formula");
     return asin(arg);
   }
   if (tree->type == ACOS) {
     double arg = eval_tree(tree->left,i);
     if (arg < -1.0 || arg > 1.0)
-      error->all("Arccos of invalid value in variable formula");
+      error->one("Arccos of invalid value in variable formula");
     return acos(arg);
   }
   if (tree->type == ATAN)
@@ -1463,6 +1490,7 @@ double Variable::eval_tree(Tree *tree, int i)
     double upper = eval_tree(tree->middle,i);
     if (randomatom == NULL) {
       int seed = static_cast<int> (eval_tree(tree->right,i));
+      if (seed <= 0) error->one("Invalid math function in variable formula");
       randomatom = new RanMars(lmp,seed+me);
     }
     return randomatom->uniform()*(upper-lower)+lower;
@@ -1470,11 +1498,13 @@ double Variable::eval_tree(Tree *tree, int i)
   if (tree->type == NORMAL) {
     double mu = eval_tree(tree->left,i);
     double sigma = eval_tree(tree->middle,i);
+    if (sigma < 0.0) error->one("Invalid math function in variable formula");
     if (randomatom == NULL) {
       int seed = static_cast<int> (eval_tree(tree->right,i));
+      if (seed <= 0) error->one("Invalid math function in variable formula");
       randomatom = new RanMars(lmp,seed+me);
     }
-    return randomatom->gaussian()*sigma + mu;
+    return mu + sigma*randomatom->gaussian();
   }
 
   if (tree->type == CEIL)
@@ -1484,6 +1514,76 @@ double Variable::eval_tree(Tree *tree, int i)
   if (tree->type == ROUND)
     return MYROUND(eval_tree(tree->left,i));
 
+  if (tree->type == RAMP) {
+    double arg1 = eval_tree(tree->left,i);
+    double arg2 = eval_tree(tree->right,i);
+    double delta = update->ntimestep - update->beginstep;
+    delta /= update->endstep - update->beginstep;
+    double arg = arg1 + delta*(arg2-arg1);
+    return arg;
+  }
+
+  if (tree->type == STAGGER) {
+    int ivalue1 = static_cast<int> (eval_tree(tree->left,i));
+    int ivalue2 = static_cast<int> (eval_tree(tree->right,i));
+    if (ivalue1 <= 0 || ivalue2 <= 0 || ivalue1 <= ivalue2)
+      error->one("Invalid math function in variable formula");
+    int lower = update->ntimestep/ivalue1 * ivalue1;
+    int delta = update->ntimestep - lower;
+    double arg;
+    if (delta < ivalue2) arg = lower+ivalue2;
+    else arg = lower+ivalue1;
+    return arg;
+  }
+
+  if (tree->type == LOGFREQ) {
+    int ivalue1 = static_cast<int> (eval_tree(tree->left,i));
+    int ivalue2 = static_cast<int> (eval_tree(tree->middle,i));
+    int ivalue3 = static_cast<int> (eval_tree(tree->right,i));
+    if (ivalue1 <= 0 || ivalue2 <= 0 || ivalue3 <= 0 || ivalue2 >= ivalue3)
+      error->one("Invalid math function in variable formula");
+    double arg;
+    if (update->ntimestep < ivalue1) arg = ivalue1;
+    else {
+      int lower = ivalue1;
+      while (update->ntimestep >= ivalue3*lower) lower *= ivalue3;
+      int multiple = update->ntimestep/lower;
+      if (multiple < ivalue2) arg = (multiple+1)*lower;
+      else arg = lower*ivalue3;
+    }
+    return arg;
+  }
+
+  if (tree->type == VLINEAR) {
+    double arg1 = eval_tree(tree->left,i);
+    double arg2 = eval_tree(tree->right,i);
+    double delta = update->ntimestep - update->beginstep;
+    double arg = arg1 + arg2*delta*update->dt;
+    return arg;
+  }
+
+  if (tree->type == SWIGGLE) {
+    double arg1 = eval_tree(tree->left,i);
+    double arg2 = eval_tree(tree->middle,i);
+    double arg3 = eval_tree(tree->right,i);
+    if (arg3 == 0.0) error->one("Invalid math function in variable formula");
+    double delta = update->ntimestep - update->beginstep;
+    double omega = 2.0*PI/arg3;
+    double arg = arg1 + arg2*sin(omega*delta*update->dt);
+    return arg;
+  }
+
+  if (tree->type == CWIGGLE) {
+    double arg1 = eval_tree(tree->left,i);
+    double arg2 = eval_tree(tree->middle,i);
+    double arg3 = eval_tree(tree->right,i);
+    if (arg3 == 0.0) error->one("Invalid math function in variable formula");
+    double delta = update->ntimestep - update->beginstep;
+    double omega = 2.0*PI/arg3;
+    double arg = arg1 + arg2*(1.0-cos(omega*delta*update->dt));
+    return arg;
+  }
+
   return 0.0;
 }
 
@@ -1563,9 +1663,11 @@ int Variable::int_between_brackets(char *&ptr)
    word = math function
    contents = str between parentheses with one,two,three args
    return 0 if not a match, 1 if successfully processed
-   customize by adding a math function in 2 places:
+   customize by adding a math function:
      sqrt(),exp(),ln(),log(),sin(),cos(),tan(),asin(),acos(),atan(),
-     atan2(y,x),random(x,y,z),normal(x,y,z),ceil(),floor(),round()
+     atan2(y,x),random(x,y,z),normal(x,y,z),ceil(),floor(),round(),
+     ramp(x,y),stagger(x,y),logfreq(x,y,z),
+     vlinear(x,y),swiggle(x,y,z),cwiggle(x,y,z)
 ------------------------------------------------------------------------- */
 
 int Variable::math_function(char *word, char *contents, Tree **tree,
@@ -1581,7 +1683,10 @@ int Variable::math_function(char *word, char *contents, Tree **tree,
       strcmp(word,"acos") && strcmp(word,"atan") && 
       strcmp(word,"atan2") && strcmp(word,"random") && 
       strcmp(word,"normal") && strcmp(word,"ceil") && 
-      strcmp(word,"floor") && strcmp(word,"round"))
+      strcmp(word,"floor") && strcmp(word,"round") &&
+      strcmp(word,"ramp") && strcmp(word,"stagger") &&
+      strcmp(word,"logfreq") && strcmp(word,"vlinear") &&
+      strcmp(word,"swiggle") && strcmp(word,"cwiggle"))
     return 0;
 
   // parse contents for arg1,arg2,arg3 separated by commas
@@ -1728,6 +1833,7 @@ int Variable::math_function(char *word, char *contents, Tree **tree,
     else {
       if (randomequal == NULL) {
 	int seed = static_cast<int> (value3);
+	if (seed <= 0) error->all("Invalid math function in variable formula");
 	randomequal = new RanMars(lmp,seed);
       }
       argstack[nargstack++] = randomequal->uniform()*(value2-value1) + value1;
@@ -1736,11 +1842,14 @@ int Variable::math_function(char *word, char *contents, Tree **tree,
     if (narg != 3) error->all("Invalid math function in variable formula");
     if (tree) newtree->type = NORMAL;
     else {
+      if (value2 < 0.0) 
+	error->all("Invalid math function in variable formula");
       if (randomequal == NULL) {
 	int seed = static_cast<int> (value3);
+	if (seed <= 0) error->all("Invalid math function in variable formula");
 	randomequal = new RanMars(lmp,seed);
       }
-      argstack[nargstack++] = randomequal->gaussian()*value2 + value1;
+      argstack[nargstack++] = value1 + value2*randomequal->gaussian();
     }
 
   } else if (strcmp(word,"ceil") == 0) {
@@ -1757,6 +1866,94 @@ int Variable::math_function(char *word, char *contents, Tree **tree,
     if (narg != 1) error->all("Invalid math function in variable formula");
     if (tree) newtree->type = ROUND;
     else argstack[nargstack++] = MYROUND(value1);
+
+  } else if (strcmp(word,"ramp") == 0) {
+    if (narg != 2) error->all("Invalid math function in variable formula");
+    if (update->whichflag == 0)
+      error->all("Cannot use ramp in variable formula between runs");
+    if (tree) newtree->type = RAMP;
+    else {
+      double delta = update->ntimestep - update->beginstep;
+      delta /= update->endstep - update->beginstep;
+      double value = value1 + delta*(value2-value1);
+      argstack[nargstack++] = value;
+    }
+
+  } else if (strcmp(word,"stagger") == 0) {
+    if (narg != 2) error->all("Invalid math function in variable formula");
+    if (tree) newtree->type = STAGGER;
+    else {
+      int ivalue1 = static_cast<int> (value1);
+      int ivalue2 = static_cast<int> (value2);
+      if (ivalue1 <= 0 || ivalue2 <= 0 || ivalue1 <= ivalue2)
+	error->all("Invalid math function in variable formula");
+      int lower = update->ntimestep/ivalue1 * ivalue1;
+      int delta = update->ntimestep - lower;
+      double value;
+      if (delta < ivalue2) value = lower+ivalue2;
+      else value = lower+ivalue1;
+      argstack[nargstack++] = value;
+    }
+
+  } else if (strcmp(word,"logfreq") == 0) {
+    if (narg != 3) error->all("Invalid math function in variable formula");
+    if (tree) newtree->type = LOGFREQ;
+    else {
+      int ivalue1 = static_cast<int> (value1);
+      int ivalue2 = static_cast<int> (value2);
+      int ivalue3 = static_cast<int> (value3);
+      if (ivalue1 <= 0 || ivalue2 <= 0 || ivalue3 <= 0 || ivalue2 >= ivalue3)
+	error->all("Invalid math function in variable formula");
+      double value;
+      if (update->ntimestep < ivalue1) value = ivalue1;
+      else {
+	int lower = ivalue1;
+	while (update->ntimestep >= ivalue3*lower) lower *= ivalue3;
+	int multiple = update->ntimestep/lower;
+	if (multiple < ivalue2) value = (multiple+1)*lower;
+	else value = lower*ivalue3;
+      }
+      argstack[nargstack++] = value;
+    }
+
+  } else if (strcmp(word,"vlinear") == 0) {
+    if (narg != 2) error->all("Invalid math function in variable formula");
+    if (update->whichflag == 0)
+      error->all("Cannot use vlinear in variable formula between runs");
+    if (tree) newtree->type = VLINEAR;
+    else {
+      double delta = update->ntimestep - update->beginstep;
+      double value = value1 + value2*delta*update->dt;
+      argstack[nargstack++] = value;
+    }
+
+  } else if (strcmp(word,"swiggle") == 0) {
+    if (narg != 3) error->all("Invalid math function in variable formula");
+    if (update->whichflag == 0)
+      error->all("Cannot use swiggle in variable formula between runs");
+    if (tree) newtree->type = CWIGGLE;
+    else {
+      if (value3 == 0.0)
+	error->all("Invalid math function in variable formula");
+      double delta = update->ntimestep - update->beginstep;
+      double omega = 2.0*PI/value3;
+      double value = value1 + value2*sin(omega*delta*update->dt);
+      argstack[nargstack++] = value;
+    }
+
+  } else if (strcmp(word,"cwiggle") == 0) {
+    if (narg != 3) error->all("Invalid math function in variable formula");
+    if (update->whichflag == 0)
+      error->all("Cannot use cwiggle in variable formula between runs");
+    if (tree) newtree->type = CWIGGLE;
+    else {
+      if (value3 == 0.0)
+	error->all("Invalid math function in variable formula");
+      double delta = update->ntimestep - update->beginstep;
+      double omega = 2.0*PI/value3;
+      double value = value1 + value2*(1.0-cos(omega*delta*update->dt));
+      argstack[nargstack++] = value;
+    }
   }
 
   delete [] arg1;
@@ -1772,7 +1969,7 @@ int Variable::math_function(char *word, char *contents, Tree **tree,
    word = group function
    contents = str between parentheses with one,two,three args
    return 0 if not a match, 1 if successfully processed
-   customize by adding a group function in 2 places with optional region arg:
+   customize by adding a group function with optional region arg:
      count(group),mass(group),charge(group),
      xcm(group,dim),vcm(group,dim),fcm(group,dim),
      bound(group,xmin),gyration(group),ke(group),angmom(group),
@@ -2017,9 +2214,8 @@ int Variable::region_function(char *id)
    word = special function
    contents = str between parentheses with one,two,three args
    return 0 if not a match, 1 if successfully processed
-   customize by adding a special function in 2 places
+   customize by adding a special function:
-     ramp(x,y),stagger(x,y),logfreq(x,y,z),sum(x),min(x),max(x),
+     sum(x),min(x),max(x),ave(x),trap(x)
-     ave(x),trap(x)
 ------------------------------------------------------------------------- */
 
 int Variable::special_function(char *word, char *contents, Tree **tree,
@@ -2028,9 +2224,7 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
 {
   // word not a match to any special function
 
-  if (strcmp(word,"ramp") && strcmp(word,"stagger") && 
+  if (strcmp(word,"sum") && strcmp(word,"min") && strcmp(word,"max") &&
-      strcmp(word,"logfreq") && strcmp(word,"sum") &&
-      strcmp(word,"min") && strcmp(word,"max") &&
       strcmp(word,"ave") && strcmp(word,"trap"))
     return 0;
 
@@ -2066,164 +2260,123 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
 
   // match word to special function
 
-  double value;
+  int method;
+  if (strcmp(word,"sum") == 0) method = SUM;
+  else if (strcmp(word,"min") == 0) method = XMIN;
+  else if (strcmp(word,"max") == 0) method = XMAX;
+  else if (strcmp(word,"ave") == 0) method = AVE;
+  else if (strcmp(word,"trap") == 0) method = TRAP;
 
-  if (strcmp(word,"ramp") == 0) {
+  if (narg != 1) error->all("Invalid special function in variable formula");
-    if (narg != 2) error->all("Invalid special function in variable formula");
-    if (update->whichflag == 0)
-      error->all("Cannot use ramp in variable formula between runs");
-    double value1 = numeric(arg1);
-    double value2 = numeric(arg2);
-    double delta = update->ntimestep - update->beginstep;
-    delta /= update->endstep - update->beginstep;
-    value = value1 + delta*(value2-value1);
   
-  } else if (strcmp(word,"stagger") == 0) {
+  Compute *compute = NULL;
-    if (narg != 2) error->all("Invalid special function in variable formula");
+  Fix *fix = NULL;
-    int ivalue1 = inumeric(arg1);
+  int index,nvec,nstride;
-    int ivalue2 = inumeric(arg2);
-    if (ivalue1 <= 0 || ivalue2 <= 0 || ivalue1 <= ivalue2)
-      error->all("Invalid special function in variable formula");
-    int lower = update->ntimestep/ivalue1 * ivalue1;
-    int delta = update->ntimestep - lower;
-    if (delta < ivalue2) value = lower+ivalue2;
-    else value = lower+ivalue1;
 
-  } else if (strcmp(word,"logfreq") == 0) {
+  if (strstr(arg1,"c_") == arg1) {
-    if (narg != 3) error->all("Invalid special function in variable formula");
+    ptr1 = strchr(arg1,'[');
-    int ivalue1 = inumeric(arg1);
+    if (ptr1) {
-    int ivalue2 = inumeric(arg2);
+      ptr2 = ptr1;
-    int ivalue3 = inumeric(arg3);
+      index = int_between_brackets(ptr2);
-    if (ivalue1 <= 0 || ivalue2 <= 0 || ivalue3 <= 0 || ivalue2 >= ivalue3)
+      *ptr1 = '\0';
-      error->all("Invalid special function in variable formula");
+    } else index = 0;
-    if (update->ntimestep < ivalue1) value = ivalue1;
-    else {
-      int lower = ivalue1;
-      while (update->ntimestep >= ivalue3*lower) lower *= ivalue3;
-      int multiple = update->ntimestep/lower;
-      if (multiple < ivalue2) value = (multiple+1)*lower;
-      else value = lower*ivalue3;
-    }
 
-  } else {
+    int icompute = modify->find_compute(&arg1[2]);
-    if (narg != 1) error->all("Invalid special function in variable formula");
+    if (icompute < 0) error->all("Invalid compute ID in variable formula");
-
+    compute = modify->compute[icompute];
-    int method;
+    if (index == 0 && compute->vector_flag) {
-    if (strcmp(word,"sum") == 0) method = SUM;
+      if (update->whichflag == 0) {
-    else if (strcmp(word,"min") == 0) method = XMIN;
+	if (compute->invoked_vector != update->ntimestep)
-    else if (strcmp(word,"max") == 0) method = XMAX;
+	  error->all("Compute used in variable between runs is not current");
-    else if (strcmp(word,"ave") == 0) method = AVE;
+      } else if (!(compute->invoked_flag & INVOKED_VECTOR)) {
-    else if (strcmp(word,"trap") == 0) method = TRAP;
+	compute->compute_vector();
-
+	compute->invoked_flag |= INVOKED_VECTOR;
-    Compute *compute = NULL;
-    Fix *fix = NULL;
-    int index,nvec,nstride;
-
-    if (strstr(arg1,"c_") == arg1) {
-      ptr1 = strchr(arg1,'[');
-      if (ptr1) {
-	ptr2 = ptr1;
-	index = int_between_brackets(ptr2);
-	*ptr1 = '\0';
-      } else index = 0;
-
-      int icompute = modify->find_compute(&arg1[2]);
-      if (icompute < 0) error->all("Invalid compute ID in variable formula");
-      compute = modify->compute[icompute];
-      if (index == 0 && compute->vector_flag) {
-	if (update->whichflag == 0) {
-	  if (compute->invoked_vector != update->ntimestep)
-	    error->all("Compute used in variable between runs is not current");
-	} else if (!(compute->invoked_flag & INVOKED_VECTOR)) {
-	  compute->compute_vector();
-	  compute->invoked_flag |= INVOKED_VECTOR;
-	}
-	nvec = compute->size_vector;
-	nstride = 1;
-      } else if (index && compute->array_flag) {
-	if (index > compute->size_array_cols)
-	  error->all("Variable formula compute array "
-		     "is accessed out-of-range");
-	if (update->whichflag == 0) {
-	  if (compute->invoked_array != update->ntimestep)
-	    error->all("Compute used in variable between runs is not current");
-	} else if (!(compute->invoked_flag & INVOKED_ARRAY)) {
-	  compute->compute_array();
-	  compute->invoked_flag |= INVOKED_ARRAY;
-	}
-	nvec = compute->size_array_rows;
-	nstride = compute->size_array_cols;
-      } else error->all("Mismatched compute in variable formula");
-
-    } else if (strstr(arg1,"f_") == arg1) {
-      ptr1 = strchr(arg1,'[');
-      if (ptr1) {
-	ptr2 = ptr1;
-	index = int_between_brackets(ptr2);
-	*ptr1 = '\0';
-      } else index = 0;
-
-      int ifix = modify->find_fix(&arg1[2]);
-      if (ifix < 0) error->all("Invalid fix ID in variable formula");
-      fix = modify->fix[ifix];
-      if (index == 0 && fix->vector_flag) {
-	if (update->whichflag > 0 && update->ntimestep % fix->global_freq)
-	  error->all("Fix in variable not computed at compatible time");
-	nvec = fix->size_vector;
-	nstride = 1;
-      } else if (index && fix->array_flag) {
-	if (index > fix->size_array_cols)
-	  error->all("Variable formula fix array is accessed out-of-range");
-	if (update->whichflag > 0 && update->ntimestep % fix->global_freq)
-	  error->all("Fix in variable not computed at compatible time");
-	nvec = fix->size_array_rows;
-	nstride = fix->size_array_cols;
-      } else error->all("Mismatched fix in variable formula");
-
-    } else error->all("Invalid special function in variable formula");
-    
-    value = 0.0;
-    if (method == XMIN) value = BIG;
-    if (method == XMAX) value = -BIG;
-
-    if (compute) {
-      double *vec;
-      if (index) vec = &compute->array[0][index];
-      else vec = compute->vector; 
-
-      int j = 0;
-      for (int i = 0; i < nvec; i++) {
-	if (method == SUM) value += vec[j];
-	else if (method == XMIN) value = MIN(value,vec[j]);
-	else if (method == XMAX) value = MAX(value,vec[j]);
-	else if (method == AVE) value += vec[j];
-	else if (method == TRAP) {
-	  if (i > 0 && i < nvec-1) value += vec[j];
-	  else value += 0.5*vec[j];
-	}
-	j += nstride;
       }
-    }
+      nvec = compute->size_vector;
-
+      nstride = 1;
-    if (fix) {
+    } else if (index && compute->array_flag) {
-      double one;
+      if (index > compute->size_array_cols)
-      for (int i = 0; i < nvec; i++) {
+	error->all("Variable formula compute array "
-	if (index) one = fix->compute_array(i,index-1);
+		   "is accessed out-of-range");
-	else one = fix->compute_vector(i);
+      if (update->whichflag == 0) {
-	if (method == SUM) value += one;
+	if (compute->invoked_array != update->ntimestep)
-	else if (method == XMIN) value = MIN(value,one);
+	  error->all("Compute used in variable between runs is not current");
-	else if (method == XMAX) value = MAX(value,one);
+      } else if (!(compute->invoked_flag & INVOKED_ARRAY)) {
-	else if (method == AVE) value += one;
+	compute->compute_array();
-	else if (method == TRAP) {
+	compute->invoked_flag |= INVOKED_ARRAY;
-	  if (i > 1 && i < nvec) value += one;
-	  else value += 0.5*one;
-	}
       }
-    }
+      nvec = compute->size_array_rows;
+      nstride = compute->size_array_cols;
+    } else error->all("Mismatched compute in variable formula");
     
-    if (method == AVE) value /= nvec;
+  } else if (strstr(arg1,"f_") == arg1) {
+    ptr1 = strchr(arg1,'[');
+    if (ptr1) {
+      ptr2 = ptr1;
+      index = int_between_brackets(ptr2);
+      *ptr1 = '\0';
+    } else index = 0;
+    
+    int ifix = modify->find_fix(&arg1[2]);
+    if (ifix < 0) error->all("Invalid fix ID in variable formula");
+    fix = modify->fix[ifix];
+    if (index == 0 && fix->vector_flag) {
+      if (update->whichflag > 0 && update->ntimestep % fix->global_freq)
+	error->all("Fix in variable not computed at compatible time");
+      nvec = fix->size_vector;
+      nstride = 1;
+    } else if (index && fix->array_flag) {
+      if (index > fix->size_array_cols)
+	error->all("Variable formula fix array is accessed out-of-range");
+      if (update->whichflag > 0 && update->ntimestep % fix->global_freq)
+	error->all("Fix in variable not computed at compatible time");
+      nvec = fix->size_array_rows;
+      nstride = fix->size_array_cols;
+    } else error->all("Mismatched fix in variable formula");
+    
+  } else error->all("Invalid special function in variable formula");
+  
+  double value = 0.0;
+  if (method == XMIN) value = BIG;
+  if (method == XMAX) value = -BIG;
+  
+  if (compute) {
+    double *vec;
+    if (index) vec = &compute->array[0][index];
+    else vec = compute->vector; 
+    
+    int j = 0;
+    for (int i = 0; i < nvec; i++) {
+      if (method == SUM) value += vec[j];
+      else if (method == XMIN) value = MIN(value,vec[j]);
+      else if (method == XMAX) value = MAX(value,vec[j]);
+      else if (method == AVE) value += vec[j];
+      else if (method == TRAP) {
+	if (i > 0 && i < nvec-1) value += vec[j];
+	else value += 0.5*vec[j];
+      }
+      j += nstride;
+    }
   }
   
+  if (fix) {
+    double one;
+    for (int i = 0; i < nvec; i++) {
+      if (index) one = fix->compute_array(i,index-1);
+      else one = fix->compute_vector(i);
+      if (method == SUM) value += one;
+      else if (method == XMIN) value = MIN(value,one);
+      else if (method == XMAX) value = MAX(value,one);
+      else if (method == AVE) value += one;
+      else if (method == TRAP) {
+	if (i > 1 && i < nvec) value += one;
+	else value += 0.5*one;
+      }
+    }
+  }
+  
+  if (method == AVE) value /= nvec;
+    
   delete [] arg1;
   delete [] arg2;
   delete [] arg3;
@@ -2247,7 +2400,8 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
    flag = 1 -> vector is a per-atom compute or fix quantity with nstride
    id = positive global ID of atom, converted to local index
    push result onto tree or arg stack
-   customize by adding an atom vector: mass,type,x,y,z,vx,vy,vz,fx,fy,fz
+   customize by adding an atom vector:
+     mass,type,x,y,z,vx,vy,vz,fx,fy,fz
 ------------------------------------------------------------------------- */
 
 void Variable::peratom2global(int flag, char *word,
@@ -2300,7 +2454,8 @@ void Variable::peratom2global(int flag, char *word,
 /* ----------------------------------------------------------------------
    check if word matches an atom vector
    return 1 if yes, else 0
-   customize by adding an atom vector: mass,type,x,y,z,vx,vy,vz,fx,fy,fz
+   customize by adding an atom vector:
+     mass,type,x,y,z,vx,vy,vz,fx,fy,fz
 ------------------------------------------------------------------------- */
 
 int Variable::is_atom_vector(char *word)
@@ -2323,7 +2478,8 @@ int Variable::is_atom_vector(char *word)
    process an atom vector in formula
    push result onto tree
    word = atom vector
-   customize by adding an atom vector: mass,type,x,y,z,vx,vy,vz,fx,fy,fz
+   customize by adding an atom vector:
+     mass,type,x,y,z,vx,vy,vz,fx,fy,fz
 ------------------------------------------------------------------------- */
 
 void Variable::atom_vector(char *word, Tree **tree,
@@ -2362,6 +2518,29 @@ void Variable::atom_vector(char *word, Tree **tree,
   else if (strcmp(word,"fz") == 0) newtree->array = &atom->f[0][2];
 }
 
+/* ----------------------------------------------------------------------
+   check if word matches a constant
+   return 1 if yes, else 0
+   customize by adding a constant: PI
+------------------------------------------------------------------------- */
+
+int Variable::is_constant(char *word)
+{
+  if (strcmp(word,"PI") == 0) return 1;
+  return 0;
+}
+
+/* ----------------------------------------------------------------------
+   process a constant in formula
+   customize by adding a constant: PI
+------------------------------------------------------------------------- */
+
+double Variable::constant(char *word)
+{
+  if (strcmp(word,"PI") == 0) return PI;
+  return 0.0;
+}
+
 /* ----------------------------------------------------------------------
    read a floating point value from a string
    generate an error if not a legitimate floating point value

src/variable.h

@@ -43,6 +43,7 @@ class Variable : protected Pointers {
   int *which;              // next available value for each variable
   int *pad;                // 1 = pad loop/uloop variables with 0s, 0 = no pad
   char ***data;            // str value of each variable's values
+  double PI;
 
   class RanMars *randomequal;   // random number generator for equal-style vars
   class RanMars *randomatom;    // random number generator for atom-style vars
@@ -74,6 +75,8 @@ class Variable : protected Pointers {
 		      Tree **, Tree **, int &, double *, int &);
   int is_atom_vector(char *);
   void atom_vector(char *, Tree **, Tree **, int &);
+  int is_constant(char *);
+  double constant(char *);
   double numeric(char *);
   int inumeric(char *);
   void print_tree(Tree *, int);