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

src/domain.cpp

@@ -119,6 +119,10 @@ void Domain::init()
 	deform_groupbit = modify->fix[i]->groupbit;
       } else deform_remap = 0;
     }
+
+  // region inits
+
+  for (int i = 0; i < nregion; i++) regions[i]->init();
 }
 
 /* ----------------------------------------------------------------------

src/fix_wall_harmonic.cpp

@@ -0,0 +1,63 @@
+/* ----------------------------------------------------------------------
+   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 "fix_wall_harmonic.h"
+#include "atom.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+FixWallHarmonic::FixWallHarmonic(LAMMPS *lmp, int narg, char **arg) : 
+  FixWall(lmp, narg, arg) {}
+
+/* ----------------------------------------------------------------------
+   interaction of all particles in group with all 6 walls (if defined)
+   error if any particle is on or behind wall
+------------------------------------------------------------------------- */
+
+void FixWallHarmonic::wall_particle(int m, double coord)
+{
+  double delta,dr,fwall;
+
+  double **x = atom->x;
+  double **f = atom->f;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  int dim = m/2;
+  int side = m % 2;
+  if (side == 0) side = -1;
+
+  int onflag = 0;
+
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+      if (side < 0) delta = x[i][dim] - coord;
+      else delta = coord - x[i][dim];
+      if (delta >= cutoff[m]) continue;
+      if (delta <= 0.0) {
+	onflag = 1;
+	continue;
+      }
+      dr = cutoff[m]-delta;
+      fwall = side * 2.0*epsilon[m]*dr;
+      f[i][dim] -= fwall;
+      ewall[0] += epsilon[m]*dr*dr;
+      ewall[m+1] += fwall;
+    }
+
+  if (onflag) error->one("Particle on or inside fix wall surface");
+}

src/fix_wall_harmonic.h

@@ -0,0 +1,33 @@
+/* ----------------------------------------------------------------------
+   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_WALL_HARMONIC_H
+#define FIX_WALL_HARMONIC_H
+
+#include "fix_wall.h"
+
+namespace LAMMPS_NS {
+
+class FixWallHarmonic : public FixWall {
+ public:
+  FixWallHarmonic(class LAMMPS *, int, char **);
+  void precompute(int) {}
+  void wall_particle(int, double);
+
+ private:
+  double offset[6];
+};
+
+}
+
+#endif

src/fix_wall_region.cpp

@@ -204,7 +204,7 @@ void FixWallRegion::post_force(int vflag)
       if (style == COLLOID) tooclose = shape[type[i]][0];
       else tooclose = 0.0;
 
-      n = region->surface(x[i],cutoff);
+      n = region->surface(x[i][0],x[i][1],x[i][2],cutoff);
 
       for (m = 0; m < n; m++) {
 	if (region->contact[m].r <= tooclose) {

src/region.cpp

@@ -15,12 +15,15 @@
 #include "stdlib.h"
 #include "string.h"
 #include "region.h"
+#include "update.h"
 #include "domain.h"
 #include "lattice.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
 
+enum{NONE,LINEAR,WIGGLE,ROTATE,VARIABLE};
+
 /* ---------------------------------------------------------------------- */
 
 Region::Region(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
@@ -32,6 +35,8 @@ Region::Region(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
   n = strlen(arg[1]) + 1;
   style = new char[n];
   strcpy(style,arg[1]);
+
+  time_origin = update->ntimestep;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -42,6 +47,13 @@ Region::~Region()
   delete [] style;
 }
 
+/* ---------------------------------------------------------------------- */
+
+void Region::init()
+{
+  dt = update->dt;
+}
+
 /* ----------------------------------------------------------------------
    parse optional parameters at end of region input line
 ------------------------------------------------------------------------- */
@@ -54,6 +66,7 @@ void Region::options(int narg, char **arg)
 
   interior = 1;
   scaleflag = 1;
+  dynamic = NONE;
 
   int iarg = 0;
   while (iarg < narg) {
@@ -69,9 +82,41 @@ void Region::options(int narg, char **arg)
       else if (strcmp(arg[iarg+1],"out") == 0) interior = 0;
       else error->all("Illegal region command");
       iarg += 2;
+    } else if (strcmp(arg[iarg],"linear") == 0) {
+      if (iarg+4 > narg) error->all("Illegal region command");
+      vx = atof(arg[iarg+1]);
+      vy = atof(arg[iarg+2]);
+      vz = atof(arg[iarg+3]);
+      dynamic = LINEAR;
+      iarg += 4;
+    } else if (strcmp(arg[iarg],"wiggle") == 0) {
+      if (iarg+5 > narg) error->all("Illegal region command");
+      ax = atof(arg[iarg+1]);
+      ay = atof(arg[iarg+2]);
+      az = atof(arg[iarg+3]);
+      period = atof(arg[iarg+4]);
+      dynamic = WIGGLE;
+      iarg += 5;
+    } else if (strcmp(arg[iarg],"rotate") == 0) {
+      if (iarg+8 > narg) error->all("Illegal region command");
+      point[0] = atof(arg[iarg+1]);
+      point[1] = atof(arg[iarg+2]);
+      point[2] = atof(arg[iarg+3]);
+      axis[0] = atof(arg[iarg+4]);
+      axis[1] = atof(arg[iarg+5]);
+      axis[2] = atof(arg[iarg+6]);
+      period = atof(arg[iarg+7]);
+      dynamic = ROTATE;
+      iarg += 8;
     } else error->all("Illegal region command");
   }
 
+  // error check
+  
+  if (dynamic && 
+      (strcmp(style,"union") == 0 || strcmp(style,"intersect") == 0))
+    error->all("Region union or intersect cannot be dynamic");
+
   // setup scaling
 
   if (scaleflag && domain->lattice == NULL)
@@ -83,16 +128,128 @@ void Region::options(int narg, char **arg)
     zscale = domain->lattice->zlattice;
   }
   else xscale = yscale = zscale = 1.0;
+
+  if (dynamic == LINEAR) {
+    vx *= xscale;
+    vy *= yscale;
+    vz *= zscale;
+  } else if (dynamic == WIGGLE) {
+    ax *= xscale;
+    ay *= yscale;
+    az *= zscale;
+  } else if (dynamic == ROTATE) {
+    point[0] *= xscale;
+    point[1] *= yscale;
+    point[2] *= zscale;
+  }
+
+  if (dynamic == WIGGLE || dynamic == ROTATE) {
+    double PI = 4.0 * atan(1.0);
+    omega_rotate = 2.0*PI / period;
+  }
+
+  // runit = unit vector along rotation axis
+
+  if (dynamic == ROTATE) {
+    double len = sqrt(axis[0]*axis[0] + axis[1]*axis[1] + axis[2]*axis[2]);
+    if (len == 0.0) 
+      error->all("Region cannot have 0 length rotation vector");
+    runit[0] = axis[0]/len;
+    runit[1] = axis[1]/len;
+    runit[2] = axis[2]/len;
+  }
+}
+
+/* ----------------------------------------------------------------------
+   determine if point x,y,z is a match to region volume
+   XOR computes 0 if 2 args are the same, 1 if different
+   note that inside() returns 1 for points on surface of region
+   thus point on surface of exterior region will not match
+   if region is dynamic, apply inverse of region change to x
+------------------------------------------------------------------------- */
+
+int Region::match(double x, double y, double z)
+{
+  double a[3],b[3],c[3],d[3],disp[0];
+
+  if (dynamic) {
+    double delta = (update->ntimestep - time_origin) * dt;
+    if (dynamic == LINEAR) {
+      x -= vx*delta;
+      y -= vy*delta;
+      z -= vz*delta;
+    } else if (dynamic == WIGGLE) {
+      double arg = omega_rotate * delta;
+      double sine = sin(arg);
+      x -= ax*sine;
+      y -= ay*sine;
+      z -= az*sine;
+    } else if (dynamic == ROTATE) {
+      double angle = -omega_rotate*delta;
+      rotate(x,y,z,angle);
+    }
+  }
+
+  return !(inside(x,y,z) ^ interior);
 }
 
 /* ----------------------------------------------------------------------
    generate list of contact points for interior or exterior regions
+   if region is dynamic:
+     change x by inverse of region change
+     change contact point by region change
 ------------------------------------------------------------------------- */
 
-int Region::surface(double *x, double cutoff)
+int Region::surface(double x, double y, double z, double cutoff)
 {
-  if (interior) return surface_interior(x,cutoff);
+  int ncontact;
-  return surface_exterior(x,cutoff);
+  double xnear[3],xhold[3];
+
+  if (dynamic) {
+    double delta = (update->ntimestep - time_origin) * dt;
+    if (dynamic == LINEAR) {
+      x -= vx*delta;
+      y -= vy*delta;
+      z -= vz*delta;
+    } else if (dynamic == WIGGLE) {
+      double arg = omega_rotate * delta;
+      double sine = sin(arg);
+      x -= ax*sine;
+      y -= ay*sine;
+      z -= az*sine;
+    } else if (dynamic == ROTATE) {
+      xhold[0] = x;
+      xhold[1] = y;
+      xhold[2] = z;
+      double angle = -omega_rotate*delta;
+      rotate(x,y,z,angle);
+    }
+  }
+
+  xnear[0] = x;
+  xnear[1] = y;
+  xnear[2] = z;
+
+  if (interior) ncontact = surface_interior(xnear,cutoff);
+  else ncontact = surface_exterior(xnear,cutoff);
+
+  if (dynamic && ncontact) {
+    double delta = (update->ntimestep - time_origin) * dt;
+    if (dynamic == ROTATE) {
+      for (int i = 0; i < ncontact; i++) {
+	x -= contact[i].delx;
+	y -= contact[i].dely;
+	z -= contact[i].delz;
+	double angle = omega_rotate*delta;
+	rotate(x,y,z,angle);
+	contact[i].delx = xhold[0] - x;
+	contact[i].dely = xhold[1] - y;
+	contact[i].delz = xhold[2] - z;
+      }
+    }
+  }
+
+  return ncontact;
 }
 
 /* ----------------------------------------------------------------------
@@ -111,3 +268,47 @@ void Region::add_contact(int n, double *x, double xp, double yp, double zp)
   contact[n].dely = dely;
   contact[n].delz = delz;
 }
+
+/* ----------------------------------------------------------------------
+   rotate x,y,z by angle via right-hand rule around point and runit normal
+   sign of angle determines whether rotating forward/backward in time
+   return updated x,y,z
+   P = point = vector = point of rotation
+   R = vector = axis of rotation
+   w = omega of rotation (from period)
+   X0 = x,y,z = initial coord of atom
+   R0 = runit = 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,y,z = P + C + A cos(w*dt) + B sin(w*dt)
+------------------------------------------------------------------------- */
+
+void Region::rotate(double &x, double &y, double &z, double angle)
+{
+  double a[3],b[3],c[3],d[3],disp[0];
+
+  double sine = sin(angle);
+  double cosine = cos(angle);
+  double x0dotr = x*runit[0] + y*runit[1] + z*runit[2]; 
+  c[0] = x0dotr * runit[0];
+  c[1] = x0dotr * runit[1];
+  c[2] = x0dotr * runit[2];
+  d[0] = x - point[0];
+  d[1] = y - point[1];
+  d[2] = z - 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 = point[0] + c[0] + disp[0];
+  y = point[1] + c[1] + disp[1];
+  z = point[2] + c[2] + disp[2];
+}

src/region.h

@@ -28,6 +28,7 @@ class Region : protected Pointers {
   double extent_ylo,extent_yhi;
   double extent_zlo,extent_zhi;
   int bboxflag;                     // 1 if bounding box is computable
+  int dynamic;                      // 1 if region changes over time
 
   // contact = particle near region surface
 
@@ -40,14 +41,26 @@ class Region : protected Pointers {
  
   Region(class LAMMPS *, int, char **);
   virtual ~Region();
-  virtual int match(double, double, double) = 0;
+  void init();
-  int surface(double *, double);
+  int match(double, double, double);
+  int surface(double, double, double, double);
+
+  virtual int inside(double, double, double) = 0;
   virtual int surface_interior(double *, double) = 0;
   virtual int surface_exterior(double *, double) = 0;
 
  protected:
   void options(int, char **);
   void add_contact(int, double *, double, double, double);
+
+ private:
+  int time_origin;
+  double dt,period,omega_rotate;
+  double vx,vy,vz;
+  double ax,ay,az;
+  double point[3],axis[3],runit[3];
+
+  void rotate(double &, double &, double &, double);
 };
 
 }

src/region_block.cpp

@@ -29,7 +29,7 @@ using namespace LAMMPS_NS;
 RegBlock::RegBlock(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg)
 {
   options(narg-8,&arg[8]);
-
+  
   if (strcmp(arg[2],"INF") == 0 || strcmp(arg[2],"EDGE") == 0) {
     if (domain->box_exist == 0) 
       error->all("Cannot use region INF or EDGE when box does not exist");
@@ -113,13 +113,11 @@ RegBlock::~RegBlock()
    inside = 0 if x,y,z is outside and not on surface
 ------------------------------------------------------------------------- */
 
-int RegBlock::match(double x, double y, double z)
+int RegBlock::inside(double x, double y, double z)
 {
-  int inside = 0;
   if (x >= xlo && x <= xhi && y >= ylo && y <= yhi && z >= zlo && z <= zhi)
-    inside = 1;
+    return 1;
-
+  return 0;
-  return !(inside ^ interior);         // 1 if same, 0 if different
 }
 
 /* ----------------------------------------------------------------------

src/region_block.h

@@ -24,7 +24,7 @@ class RegBlock : public Region {
  public:
   RegBlock(class LAMMPS *, int, char **);
   ~RegBlock();
-  int match(double, double, double);
+  int inside(double, double, double);
   int surface_interior(double *, double);
   int surface_exterior(double *, double);
 

src/region_cone.cpp

@@ -162,7 +162,7 @@ RegCone::~RegCone()
    inside = 0 if x,y,z is outside and not on surface
 ------------------------------------------------------------------------- */
 
-int RegCone::match(double x, double y, double z)
+int RegCone::inside(double x, double y, double z)
 {
   double del1,del2,dist;
   double currentradius;
@@ -193,7 +193,7 @@ int RegCone::match(double x, double y, double z)
     else inside = 0;
   }
 
-  return !(inside ^ interior);         // 1 if same, 0 if different
+  return inside;
 }
 
 /* ----------------------------------------------------------------------

src/region_cone.h

@@ -22,7 +22,7 @@ class RegCone : public Region {
  public:
   RegCone(class LAMMPS *, int, char **);
   ~RegCone();
-  int match(double, double, double);
+  int inside(double, double, double);
   int surface_interior(double *, double);
   int surface_exterior(double *, double);
 

src/region_cylinder.cpp

@@ -150,7 +150,7 @@ RegCylinder::~RegCylinder()
    inside = 0 if x,y,z is outside and not on surface
 ------------------------------------------------------------------------- */
 
-int RegCylinder::match(double x, double y, double z)
+int RegCylinder::inside(double x, double y, double z)
 {
   double del1,del2,dist;
   int inside;
@@ -175,7 +175,7 @@ int RegCylinder::match(double x, double y, double z)
     else inside = 0;
   }
 
-  return !(inside ^ interior);         // 1 if same, 0 if different
+  return inside;
 }
 
 /* ----------------------------------------------------------------------

src/region_cylinder.h

@@ -24,7 +24,7 @@ class RegCylinder : public Region {
  public:
   RegCylinder(class LAMMPS *, int, char **);
   ~RegCylinder();
-  int match(double, double, double);
+  int inside(double, double, double);
   int surface_interior(double *, double);
   int surface_exterior(double *, double);
 

src/region_intersect.cpp

@@ -98,17 +98,15 @@ RegIntersect::~RegIntersect()
    else inside = 0
 ------------------------------------------------------------------------- */
 
-int RegIntersect::match(double x, double y, double z)
+int RegIntersect::inside(double x, double y, double z)
 {
   int ilist;
   Region **regions = domain->regions;
   for (ilist = 0; ilist < nregion; ilist++)
     if (!regions[list[ilist]]->match(x,y,z)) break;
 
-  int inside = 0;                      // inside if matched all regions
+  if (ilist == nregion) return 1;
-  if (ilist == nregion) inside = 1;
+  return 0;
-
-  return !(inside ^ interior);         // 1 if same, 0 if different
 }
 
 /* ----------------------------------------------------------------------
@@ -127,7 +125,7 @@ int RegIntersect::surface_interior(double *x, double cutoff)
 
   for (ilist = 0; ilist < nregion; ilist++) {
     iregion = list[ilist];
-    ncontacts = regions[iregion]->surface(x,cutoff);
+    ncontacts = regions[iregion]->surface(x[0],x[1],x[2],cutoff);
     for (m = 0; m < ncontacts; m++) {
       xs = x[0] - regions[iregion]->contact[m].delx;
       ys = x[1] - regions[iregion]->contact[m].dely;
@@ -172,7 +170,7 @@ int RegIntersect::surface_exterior(double *x, double cutoff)
 
   for (ilist = 0; ilist < nregion; ilist++) {
     iregion = list[ilist];
-    ncontacts = regions[iregion]->surface(x,cutoff);
+    ncontacts = regions[iregion]->surface(x[0],x[1],x[2],cutoff);
     for (m = 0; m < ncontacts; m++) {
       xs = x[0] - regions[iregion]->contact[m].delx;
       ys = x[1] - regions[iregion]->contact[m].dely;

src/region_intersect.h

@@ -22,7 +22,7 @@ class RegIntersect : public Region {
  public:
   RegIntersect(class LAMMPS *, int, char **);
   ~RegIntersect();
-  int match(double, double, double);
+  int inside(double, double, double);
   int surface_interior(double *, double);
   int surface_exterior(double *, double);
 

src/region_plane.cpp

@@ -0,0 +1,110 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "stdlib.h"
+#include "string.h"
+#include "region_plane.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+RegPlane::RegPlane(LAMMPS *lmp, int narg, char **arg) :
+  Region(lmp, narg, arg)
+{
+  options(narg-8,&arg[8]);
+
+  xp = xscale*atof(arg[2]);
+  yp = yscale*atof(arg[3]);
+  zp = zscale*atof(arg[4]);
+  normal[0] = xscale*atof(arg[5]);
+  normal[1] = yscale*atof(arg[6]);
+  normal[2] = zscale*atof(arg[7]);
+
+  // enforce unit normal
+
+  double rsq = normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2];
+  if (rsq == 0.0) error->all("Illegal region plane command");
+  normal[0] /= sqrt(rsq);
+  normal[1] /= sqrt(rsq);
+  normal[2] /= sqrt(rsq);
+
+  // plane has no bounding box
+
+  bboxflag = 0;
+
+  cmax = 1;
+  contact = new Contact[cmax];
+}
+
+/* ---------------------------------------------------------------------- */
+
+RegPlane::~RegPlane()
+{
+  delete [] contact;
+}
+
+/* ----------------------------------------------------------------------
+   inside = 1 if x,y,z is on normal side of plane or on plane
+   inside = 0 if x,y,z is on non-normal side of plane and not on plane
+   x,y,z is inside if (x-xp) dot normal >= 0
+------------------------------------------------------------------------- */
+
+int RegPlane::inside(double x, double y, double z)
+{
+  double dot = (x-xp)*normal[0] + (y-yp)*normal[1] + (z-zp)*normal[2];
+
+  if (dot >= 0.0) return 1;
+  return 0;
+}
+
+/* ----------------------------------------------------------------------
+   one contact if 0 <= x < cutoff from normal side of plane
+   no contact if on other side (possible if called from union/intersect)
+   delxyz = vector from nearest projected point on plane to x
+------------------------------------------------------------------------- */
+
+int RegPlane::surface_interior(double *x, double cutoff)
+{
+  double dot = (x[0]-xp)*normal[0] + (x[1]-yp)*normal[1] + (x[2]-zp)*normal[2];
+  if (dot < cutoff && dot >= 0.0) {
+    contact[0].r = dot;
+    contact[0].delx = dot*normal[0];
+    contact[0].dely = dot*normal[1];
+    contact[0].delz = dot*normal[2];
+    return 1;
+  }
+  return 0;
+}
+
+/* ----------------------------------------------------------------------
+   one contact if 0 <= x < cutoff from non-normal side of plane
+   no contact if on other side (possible if called from union/intersect)
+   delxyz = vector from nearest projected point on plane to x
+------------------------------------------------------------------------- */
+
+int RegPlane::surface_exterior(double *x, double cutoff)
+{
+  double dot = (x[0]-xp)*normal[0] + (x[1]-yp)*normal[1] + (x[2]-zp)*normal[2];
+  dot = -dot;
+  if (dot < cutoff && dot >= 0.0) {
+    contact[0].r = dot;
+    contact[0].delx = dot*normal[0];
+    contact[0].dely = dot*normal[1];
+    contact[0].delz = dot*normal[2];
+    return 1;
+  }
+  return 0;
+}

src/region_plane.h

@@ -0,0 +1,36 @@
+/* ----------------------------------------------------------------------
+   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 REGION_PLANE_H
+#define REGION_PLANE_H
+
+#include "region.h"
+
+namespace LAMMPS_NS {
+
+class RegPlane : public Region {
+ public:
+  RegPlane(class LAMMPS *, int, char **);
+  ~RegPlane();
+  int inside(double, double, double);
+  int surface_interior(double *, double);
+  int surface_exterior(double *, double);
+
+ private:
+  double xp,yp,zp;
+  double normal[3];
+};
+
+}
+
+#endif

src/region_prism.cpp

@@ -240,17 +240,15 @@ RegPrism::~RegPrism()
    xyz/lo = lower-left corner of prism
 ------------------------------------------------------------------------- */
 
-int RegPrism::match(double x, double y, double z)
+int RegPrism::inside(double x, double y, double z)
 {
   double a = hinv[0][0]*(x-xlo) + hinv[0][1]*(y-ylo) + hinv[0][2]*(z-zlo);
   double b = hinv[1][1]*(y-ylo) + hinv[1][2]*(z-zlo);
   double c = hinv[2][2]*(z-zlo);
 
-  int inside = 0;
   if (a >= 0.0 && a <= 1.0 && b >= 0.0 && b <= 1.0 && c >= 0.0 && c <= 1.0)
-    inside = 1;
+    return 1;
-
+  return 0;
-  return !(inside ^ interior);         // 1 if same, 0 if different
 }
 
 /* ----------------------------------------------------------------------

src/region_prism.h

@@ -24,7 +24,7 @@ class RegPrism : public Region {
  public:
   RegPrism(class LAMMPS *, int, char **);
   ~RegPrism();
-  int match(double, double, double);
+  int inside(double, double, double);
   int surface_interior(double *, double);
   int surface_exterior(double *, double);
 

src/region_sphere.cpp

@@ -63,21 +63,19 @@ RegSphere::~RegSphere()
    inside = 0 if x,y,z is outside and not on surface
 ------------------------------------------------------------------------- */
 
-int RegSphere::match(double x, double y, double z)
+int RegSphere::	inside(double x, double y, double z)
 {
   double delx = x - xc;
   double dely = y - yc;
   double delz = z - zc;
   double r = sqrt(delx*delx + dely*dely + delz*delz);
 
-  int inside = 0;
+  if (r <= radius) return 1;
-  if (r <= radius) inside = 1;
+  return 0;
-
-  return !(inside ^ interior);         // 1 if same, 0 if different
 }
 
 /* ----------------------------------------------------------------------
-   one contact if 0 <= x < cutoff away from inner surface of sphere
+   one contact if 0 <= x < cutoff from inner surface of sphere
    no contact if outside (possible if called from union/intersect)
    delxyz = vector from nearest point on sphere to x
    special case: no contact if x is at center of sphere
@@ -103,7 +101,7 @@ int RegSphere::surface_interior(double *x, double cutoff)
 }
 
 /* ----------------------------------------------------------------------
-   one contact if 0 <= x < cutoff away from outer surface of sphere
+   one contact if 0 <= x < cutoff from outer surface of sphere
    no contact if inside (possible if called from union/intersect)
    delxyz = vector from nearest point on sphere to x
 ------------------------------------------------------------------------- */

src/region_sphere.h

@@ -22,7 +22,7 @@ class RegSphere : public Region {
  public:
   RegSphere(class LAMMPS *, int, char **);
   ~RegSphere();
-  int match(double, double, double);
+  int inside(double, double, double);
   int surface_interior(double *, double);
   int surface_exterior(double *, double);
 

src/region_union.cpp

@@ -90,17 +90,15 @@ RegUnion::~RegUnion()
    else inside = 0
 ------------------------------------------------------------------------- */
 
-int RegUnion::match(double x, double y, double z)
+int RegUnion::inside(double x, double y, double z)
 {
   int ilist;
   Region **regions = domain->regions;
   for (ilist = 0; ilist < nregion; ilist++)
     if (regions[list[ilist]]->match(x,y,z)) break;
 
-  int inside = 1;                      // inside if matched any region
+  if (ilist == nregion) return 0;
-  if (ilist == nregion) inside = 0;
+  return 1;
-
-  return !(inside ^ interior);         // 1 if same, 0 if different
 }
 
 /* ----------------------------------------------------------------------
@@ -119,7 +117,7 @@ int RegUnion::surface_interior(double *x, double cutoff)
 
   for (ilist = 0; ilist < nregion; ilist++) {
     iregion = list[ilist];
-    ncontacts = regions[iregion]->surface(x,cutoff);
+    ncontacts = regions[iregion]->surface(x[0],x[1],x[2],cutoff);
     for (m = 0; m < ncontacts; m++) {
       xs = x[0] - regions[iregion]->contact[m].delx;
       ys = x[1] - regions[iregion]->contact[m].dely;
@@ -164,7 +162,7 @@ int RegUnion::surface_exterior(double *x, double cutoff)
 
   for (ilist = 0; ilist < nregion; ilist++) {
     iregion = list[ilist];
-    ncontacts = regions[iregion]->surface(x,cutoff);
+    ncontacts = regions[iregion]->surface(x[0],x[1],x[2],cutoff);
     for (m = 0; m < ncontacts; m++) {
       xs = x[0] - regions[iregion]->contact[m].delx;
       ys = x[1] - regions[iregion]->contact[m].dely;

src/region_union.h

@@ -22,7 +22,7 @@ class RegUnion : public Region {
  public:
   RegUnion(class LAMMPS *, int, char **);
   ~RegUnion();
-  int match(double, double, double);
+  int inside(double, double, double);
   int surface_interior(double *, double);
   int surface_exterior(double *, double);
 

src/update.cpp

@@ -18,6 +18,8 @@
 #include "force.h"
 #include "modify.h"
 #include "fix.h"
+#include "domain.h"
+#include "region.h"
 #include "compute.h"
 #include "output.h"
 #include "memory.h"
@@ -222,6 +224,7 @@ void Update::create_minimize(int narg, char **arg)
    reset timestep from input script
    do not allow dump files or a restart to be defined
    do not allow any timestep-dependent fixes to be defined
+   do not allow any dynamic regions to be defined
    reset eflag/vflag global so nothing will think eng/virial are current
    reset invoked flags of computes, so nothing will think they are current
    clear timestep list of computes that store future invocation times
@@ -241,6 +244,10 @@ void Update::reset_timestep(int narg, char **arg)
     if (modify->fix[i]->time_depend)
       error->all("Cannot reset timestep with a time-dependent fix defined");
 
+  for (int i = 0; i < domain->nregion; i++)
+    if (domain->regions[i]->dynamic)
+      error->all("Cannot reset timestep with a dynamic region defined");
+
   eflag_global = vflag_global = -1;
 
   for (int i = 0; i < modify->ncompute; i++) {