Merge pull request #516 from akohlmey/check-rigid-overlap

Implement check whether commands or styles try to change cached properties in rigid body integrators

src/change_box.cpp

@@ -316,6 +316,9 @@ void ChangeBox::command(int narg, char **arg)
 
     } else if (ops[m].style == REMAP) {
 
+      if (modify->check_rigid_group_overlap(groupbit))
+        error->warning(FLERR,"Attempting to remap atoms in rigid bodies");
+
       // convert atoms to lamda coords, using last box state
       // convert atoms back to box coords, using current box state
       // save current box state

src/delete_atoms.cpp

@@ -69,6 +69,15 @@ void DeleteAtoms::command(int narg, char **arg)
   else if (strcmp(arg[0],"porosity") == 0) delete_porosity(narg,arg);
   else error->all(FLERR,"Illegal delete_atoms command");
 
+  if (allflag) {
+    int igroup = group->find("all");
+    if ((igroup >= 0) && modify->check_rigid_group_overlap(group->bitmask[igroup]))
+      error->warning(FLERR,"Attempting to delete atoms in rigid bodies");
+  } else {
+    if (modify->check_rigid_list_overlap(dlist))
+      error->warning(FLERR,"Attempting to delete atoms in rigid bodies");
+  }
+
   // if allflag = 1, just reset atom->nlocal
   // else delete atoms one by one
 
@@ -89,16 +98,16 @@ void DeleteAtoms::command(int narg, char **arg)
     int i = 0;
     while (i < nlocal) {
       if (dlist[i]) {
-	avec->copy(nlocal-1,i,1);
-	dlist[i] = dlist[nlocal-1];
-	nlocal--;
+        avec->copy(nlocal-1,i,1);
+        dlist[i] = dlist[nlocal-1];
+        nlocal--;
       } else i++;
     }
-    
+
     atom->nlocal = nlocal;
     memory->destroy(dlist);
   }
-  
+
   // if non-molecular system and compress flag set,
   // reset atom tags to be contiguous
   // set all atom IDs to 0, call tag_extend()
@@ -201,7 +210,7 @@ void DeleteAtoms::delete_group(int narg, char **arg)
     allflag = 1;
     return;
   }
-  
+
   // allocate and initialize deletion list
 
   int nlocal = atom->nlocal;

src/displace_atoms.cpp

@@ -75,6 +75,9 @@ void DisplaceAtoms::command(int narg, char **arg)
   if (igroup == -1) error->all(FLERR,"Could not find displace_atoms group ID");
   groupbit = group->bitmask[igroup];
 
+  if (modify->check_rigid_group_overlap(groupbit))
+    error->warning(FLERR,"Attempting to displace atoms in rigid bodies");
+
   int style = -1;
   if (strcmp(arg[1],"move") == 0) style = MOVE;
   else if (strcmp(arg[1],"ramp") == 0) style = RAMP;

src/fix_heat.cpp

@@ -64,7 +64,7 @@ idregion(NULL), hstr(NULL), vheat(NULL), vscale(NULL)
   // optional args
 
   iregion = -1;
-  
+
   int iarg = 5;
   while (iarg < narg) {
     if (strcmp(arg[iarg],"region") == 0) {
@@ -126,6 +126,10 @@ void FixHeat::init()
     else error->all(FLERR,"Variable for fix heat is invalid style");
   }
 
+  // check for rigid bodies in region (done here for performance reasons)
+  if (modify->check_rigid_region_overlap(groupbit,domain->regions[iregion]))
+    error->warning(FLERR,"Cannot apply fix heat to atoms in rigid bodies");
+
   // cannot have 0 atoms in group
 
   if (group->count(igroup) == 0)

src/fix_temp_berendsen.cpp

@@ -128,6 +128,9 @@ void FixTempBerendsen::init()
     error->all(FLERR,"Temperature ID for fix temp/berendsen does not exist");
   temperature = modify->compute[icompute];
 
+  if (modify->check_rigid_group_overlap(groupbit))
+    error->warning(FLERR,"Cannot thermostat atoms in rigid bodies");
+
   if (temperature->tempbias) which = BIAS;
   else which = NOBIAS;
 }

src/fix_temp_csld.cpp

@@ -155,6 +155,9 @@ void FixTempCSLD::init()
     error->all(FLERR,"Temperature ID for fix temp/csld does not exist");
   temperature = modify->compute[icompute];
 
+  if (modify->check_rigid_group_overlap(groupbit))
+    error->warning(FLERR,"Cannot thermostat atoms in rigid bodies");
+
   if (temperature->tempbias) which = BIAS;
   else which = NOBIAS;
 }

src/modify.cpp

@@ -23,6 +23,7 @@
 #include "group.h"
 #include "update.h"
 #include "domain.h"
+#include "region.h"
 #include "input.h"
 #include "variable.h"
 #include "memory.h"
@@ -995,6 +996,99 @@ int Modify::check_package(const char *package_fix_name)
   return 1;
 }
 
+
+/* ----------------------------------------------------------------------
+   check if the group indicated by groupbit overlaps with any
+   currently existing rigid fixes. return 1 in this case otherwise 0
+------------------------------------------------------------------------- */
+
+int Modify::check_rigid_group_overlap(int groupbit)
+{
+  const int * const mask = atom->mask;
+  const int nlocal = atom->nlocal;
+  int dim;
+
+  int n = 0;
+  for (int ifix = 0; ifix < nfix; ifix++) {
+    if (strncmp("rigid",fix[ifix]->style,5) == 0) {
+      const int * const body = (const int *)fix[ifix]->extract("body",dim);
+      if ((body == NULL) || (dim != 1)) break;
+
+      for (int i=0; (i < nlocal) && (n == 0); ++i)
+        if ((mask[i] & groupbit) && (body[i] >= 0)) ++n;
+    }
+  }
+
+  int n_all = 0;
+  MPI_Allreduce(&n,&n_all,1,MPI_INT,MPI_SUM,world);
+
+  if (n_all > 0) return 1;
+  return 0;
+}
+
+/* ----------------------------------------------------------------------
+   check if the atoms in the group indicated by groupbit _and_ region
+   indicated by regionid overlap with any currently existing rigid fixes.
+   return 1 in this case, otherwise 0
+------------------------------------------------------------------------- */
+
+int Modify::check_rigid_region_overlap(int groupbit, Region *reg)
+{
+  const int * const mask = atom->mask;
+  const double * const * const x = atom->x;
+  const int nlocal = atom->nlocal;
+  int dim;
+
+  int n = 0;
+  reg->prematch();
+  for (int ifix = 0; ifix < nfix; ifix++) {
+    if (strncmp("rigid",fix[ifix]->style,5) == 0) {
+      const int * const body = (const int *)fix[ifix]->extract("body",dim);
+      if ((body == NULL) || (dim != 1)) break;
+
+      for (int i=0; (i < nlocal) && (n == 0); ++i)
+        if ((mask[i] & groupbit) && (body[i] >= 0)
+            && reg->match(x[i][0],x[i][1],x[i][2])) ++n;
+    }
+  }
+
+  int n_all = 0;
+  MPI_Allreduce(&n,&n_all,1,MPI_INT,MPI_SUM,world);
+
+  if (n_all > 0) return 1;
+  return 0;
+}
+
+/* ----------------------------------------------------------------------
+   check if the atoms in the selection list (length atom->nlocal,
+   content: 1 if atom is contained, 0 if not) overlap with currently
+   existing rigid fixes. return 1 in this case otherwise 0
+------------------------------------------------------------------------- */
+
+int Modify::check_rigid_list_overlap(int *select)
+{
+  const int * const mask = atom->mask;
+  const int nlocal = atom->nlocal;
+  int dim;
+
+  int n = 0;
+  for (int ifix = 0; ifix < nfix; ifix++) {
+    if (strncmp("rigid",fix[ifix]->style,5) == 0) {
+      const int * const body = (const int *)fix[ifix]->extract("body",dim);
+      if ((body == NULL) || (dim != 1)) break;
+
+      for (int i=0; (i < nlocal) && (n == 0); ++i)
+        if ((body[i] >= 0) && select[i]) ++n;
+    }
+  }
+
+  int n_all = 0;
+  MPI_Allreduce(&n,&n_all,1,MPI_INT,MPI_SUM,world);
+
+  if (n_all > 0) return 1;
+  return 0;
+}
+
 /* ----------------------------------------------------------------------
    add a new compute
 ------------------------------------------------------------------------- */

src/modify.h

@@ -98,6 +98,9 @@ class Modify : protected Pointers {
   int find_fix(const char *);
   int find_fix_by_style(const char *);
   int check_package(const char *);
+  int check_rigid_group_overlap(int);
+  int check_rigid_region_overlap(int, class Region *);
+  int check_rigid_list_overlap(int *);
 
   void add_compute(int, char **, int trysuffix=1);
   void modify_compute(int, char **);

src/set.cpp

@@ -585,6 +585,28 @@ void Set::set(int keyword)
     }
   }
 
+  // check if properties of atoms in rigid bodies are updated
+  // that are cached as per-body data.
+  switch (keyword) {
+  case X:
+  case Y:
+  case Z:
+  case MOLECULE:
+  case MASS:
+  case ANGMOM:
+  case SHAPE:
+  case DIAMETER:
+  case DENSITY:
+  case QUAT:
+  case IMAGE:
+    if (modify->check_rigid_list_overlap(select))
+      error->warning(FLERR,"Changing a property of atoms in rigid bodies "
+                     "that has no effect unless rigid bodies are rebuild");
+    break;
+  default: // assume no conflict for all other properties
+    break;
+  }
+
   // loop over selected atoms
 
   AtomVecEllipsoid *avec_ellipsoid =

src/velocity.cpp

@@ -68,6 +68,12 @@ void Velocity::command(int narg, char **arg)
   if (igroup == -1) error->all(FLERR,"Could not find velocity group ID");
   groupbit = group->bitmask[igroup];
 
+  // check if velocities of atoms in rigid bodies are updated
+
+  if (modify->check_rigid_group_overlap(groupbit))
+    error->warning(FLERR,"Changing velocities of atoms in rigid bodies. "
+                     "This has no effect unless rigid bodies are rebuild");
+
   // identify style
 
   if (strcmp(arg[1],"create") == 0) style = CREATE;