Merge branch 'rigid-gravity' of github.com:lammps/lammps into rigid-gravity

src/RIGID/fix_rigid.cpp

@@ -52,8 +52,8 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
   torque(NULL), quat(NULL), imagebody(NULL), fflag(NULL),
   tflag(NULL), langextra(NULL), sum(NULL), all(NULL),
   remapflag(NULL), xcmimage(NULL), eflags(NULL), orient(NULL),
-  dorient(NULL), id_dilate(NULL), random(NULL), avec_ellipsoid(NULL),
+  dorient(NULL), id_dilate(NULL), id_gravity(NULL), random(NULL), 
-  avec_line(NULL), avec_tri(NULL)
+  avec_ellipsoid(NULL), avec_line(NULL), avec_tri(NULL)
 {
   int i,ibody;
 
@@ -124,14 +124,17 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
     if (custom_flag) {
       if (narg < 5) error->all(FLERR,"Illegal fix rigid command");
 
-      // determine whether atom-style variable or atom property is used.
+      // determine whether atom-style variable or atom property is used
+
       if (strstr(arg[4],"i_") == arg[4]) {
         int is_double=0;
         int custom_index = atom->find_custom(arg[4]+2,is_double);
         if (custom_index == -1)
-          error->all(FLERR,"Fix rigid custom requires previously defined property/atom");
+          error->all(FLERR,"Fix rigid custom requires "
+                     "previously defined property/atom");
         else if (is_double)
-          error->all(FLERR,"Fix rigid custom requires integer-valued property/atom");
+          error->all(FLERR,"Fix rigid custom requires "
+                     "integer-valued property/atom");
         int minval = INT_MAX;
         int *value = atom->ivector[custom_index];
         for (i = 0; i < nlocal; i++)
@@ -163,12 +166,15 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
           if (mask[i] & groupbit)
             molecule[i] = (tagint)((tagint)value[i] - minval + 1);
         delete[] value;
+
       } else error->all(FLERR,"Unsupported fix rigid custom property");
+
     } else {
       if (atom->molecule_flag == 0)
         error->all(FLERR,"Fix rigid molecule requires atom attribute molecule");
       molecule = atom->molecule;
     }
+
     iarg = 4 + custom_flag;
 
     tagint maxmol_tag = -1;
@@ -297,6 +303,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
   }
 
   // number of linear rigid bodies is counted later
+
   nlinear = 0;
 
   // parse optional args
@@ -308,12 +315,14 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
   tstat_flag = 0;
   pstat_flag = 0;
   allremap = 1;
-  id_dilate = NULL;
   t_chain = 10;
   t_iter = 1;
   t_order = 3;
   p_chain = 10;
+
   inpfile = NULL;
+  id_gravity = NULL;
+  id_dilate = NULL;
 
   pcouple = NONE;
   pstyle = ANISO;
@@ -543,12 +552,20 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
       iarg += 2;
 
     } else if (strcmp(arg[iarg],"reinit") == 0) {
-      if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid/small command");
+      if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid command");
       if (strcmp("yes",arg[iarg+1]) == 0) reinitflag = 1;
       else if  (strcmp("no",arg[iarg+1]) == 0) reinitflag = 0;
       else error->all(FLERR,"Illegal fix rigid command");
       iarg += 2;
 
+    } else if (strcmp(arg[iarg],"gravity") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid command");
+      delete [] id_gravity;
+      int n = strlen(arg[iarg+1]) + 1;
+      id_gravity = new char[n];
+      strcpy(id_gravity,arg[iarg+1]);
+      iarg += 2;
+
     } else error->all(FLERR,"Illegal fix rigid command");
   }
 
@@ -621,6 +638,9 @@ FixRigid::~FixRigid()
 
   delete random;
   delete [] inpfile;
+  delete [] id_dilate;
+  delete [] id_gravity;
+
   memory->destroy(mol2body);
   memory->destroy(body2mol);
 
@@ -687,7 +707,7 @@ void FixRigid::init()
   avec_tri = (AtomVecTri *) atom->style_match("tri");
 
   // warn if more than one rigid fix
-  // if earlyflag, warn if any post-force fixes come after POEMS fix
+  // if earlyflag, warn if any post-force fixes come after a rigid fix
 
   int count = 0;
   for (i = 0; i < modify->nfix; i++)
@@ -701,7 +721,8 @@ void FixRigid::init()
       if (rflag && (modify->fmask[i] & POST_FORCE) &&
           !modify->fix[i]->rigid_flag) {
         char str[128];
-        snprintf(str,128,"Fix %s alters forces after fix rigid",modify->fix[i]->id);
+        snprintf(str,128,"Fix %s alters forces after fix rigid",
+                 modify->fix[i]->id);
         error->warning(FLERR,str);
       }
     }
@@ -713,12 +734,24 @@ void FixRigid::init()
     if (strcmp(modify->fix[i]->style,"npt") == 0) break;
     if (strcmp(modify->fix[i]->style,"nph") == 0) break;
   }
+
   if (i < modify->nfix) {
     for (int j = i; j < modify->nfix; j++)
       if (strcmp(modify->fix[j]->style,"rigid") == 0)
         error->all(FLERR,"Rigid fix must come before NPT/NPH fix");
   }
 
+  // add gravity forces based on gravity vector from fix
+
+  if (id_gravity) {
+    int ifix = modify->find_fix(id_gravity);
+    if (ifix < 0) error->all(FLERR,"Fix rigid cannot find fix gravity ID");
+    if (strcmp(modify->fix[ifix]->style,"gravity") != 0) 
+      error->all(FLERR,"Fix rigid gravity fix is invalid");
+    int tmp;
+    gvec = (double *) modify->fix[ifix]->extract("gvec",tmp);
+  }
+
   // timestep info
 
   dtv = update->dt;
@@ -1061,8 +1094,17 @@ void FixRigid::compute_forces_and_torques()
     torque[ibody][1] = all[ibody][4] + langextra[ibody][4];
     torque[ibody][2] = all[ibody][5] + langextra[ibody][5];
   }
-}
 
+  // add gravity force to COM of each body
+
+  if (id_gravity) {
+    for (ibody = 0; ibody < nbody; ibody++) {
+      fcm[ibody][0] += gvec[0]*masstotal[ibody];
+      fcm[ibody][1] += gvec[1]*masstotal[ibody];
+      fcm[ibody][2] += gvec[2]*masstotal[ibody];
+    }
+  }
+}
 
 /* ---------------------------------------------------------------------- */
 

src/RIGID/fix_rigid.h

@@ -67,6 +67,7 @@ class FixRigid : public Fix {
   int triclinic;
 
   char *inpfile;            // file to read rigid body attributes from
+
   int rstyle;               // SINGLE,MOLECULE,GROUP
   int setupflag;            // 1 if body properties are setup, else 0
   int earlyflag;            // 1 if forces/torques computed at post_force()
@@ -131,6 +132,9 @@ class FixRigid : public Fix {
   int dilate_group_bit;      // mask for dilation group
   char *id_dilate;           // group name to dilate
 
+  char *id_gravity;         // ID of fix gravity command to add gravity forces
+  double *gvec;             // ptr to gravity vector inside the fix
+
   class RanMars *random;
   class AtomVecEllipsoid *avec_ellipsoid;
   class AtomVecLine *avec_line;

src/RIGID/fix_rigid_small.cpp

@@ -57,7 +57,7 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
   xcmimage(NULL), displace(NULL), eflags(NULL), orient(NULL), dorient(NULL),
   avec_ellipsoid(NULL), avec_line(NULL), avec_tri(NULL), counts(NULL),
   itensor(NULL), mass_body(NULL), langextra(NULL), random(NULL),
-  id_dilate(NULL), onemols(NULL)
+  id_dilate(NULL), id_gravity(NULL), onemols(NULL)
 {
   int i;
 
@@ -107,7 +107,8 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
       bodyID = new tagint[nlocal];
       customflag = 1;
 
-      // determine whether atom-style variable or atom property is used.
+      // determine whether atom-style variable or atom property is used
+
       if (strstr(arg[4],"i_") == arg[4]) {
         int is_double=0;
         int custom_index = atom->find_custom(arg[4]+2,is_double);
@@ -356,6 +357,13 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
       p_chain = force->numeric(FLERR,arg[iarg+1]);
       iarg += 2;
 
+    } else if (strcmp(arg[iarg],"gravity") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid/small command");
+      delete [] id_gravity;
+      int n = strlen(arg[iarg+1]) + 1;
+      id_gravity = new char[n];
+      strcpy(id_gravity,arg[iarg+1]);
+      iarg += 2;
 
     } else error->all(FLERR,"Illegal fix rigid/small command");
   }
@@ -515,6 +523,8 @@ FixRigidSmall::~FixRigidSmall()
 
   delete random;
   delete [] inpfile;
+  delete [] id_dilate;
+  delete [] id_gravity;
 
   memory->destroy(langextra);
   memory->destroy(mass_body);
@@ -543,6 +553,7 @@ void FixRigidSmall::init()
   triclinic = domain->triclinic;
 
   // warn if more than one rigid fix
+  // if earlyflag, warn if any post-force fixes come after a rigid fix
 
   int count = 0;
   for (i = 0; i < modify->nfix; i++)
@@ -575,6 +586,17 @@ void FixRigidSmall::init()
         error->all(FLERR,"Rigid fix must come before NPT/NPH fix");
   }
 
+  // add gravity forces based on gravity vector from fix
+
+  if (id_gravity) {
+    int ifix = modify->find_fix(id_gravity);
+    if (ifix < 0) error->all(FLERR,"Fix rigid/small cannot find fix gravity ID");
+    if (strcmp(modify->fix[ifix]->style,"gravity") != 0) 
+      error->all(FLERR,"Fix rigid/small gravity fix is invalid");
+    int tmp;
+    gvec = (double *) modify->fix[ifix]->extract("gvec",tmp);
+  }
+
   // timestep info
 
   dtv = update->dt;
@@ -954,8 +976,20 @@ void FixRigidSmall::compute_forces_and_torques()
       tcm[2] += langextra[ibody][5];
     }
   }
-}
 
+  // add gravity force to COM of each body
+
+  if (id_gravity) {
+    double mass;
+    for (ibody = 0; ibody < nlocal_body; ibody++) {
+      mass = body[ibody].mass;
+      fcm = body[ibody].fcm;
+      fcm[0] += gvec[0]*mass;
+      fcm[1] += gvec[1]*mass;
+      fcm[2] += gvec[2]*mass;
+    }
+  }
+}
 
 /* ---------------------------------------------------------------------- */
 

src/RIGID/fix_rigid_small.h

@@ -164,6 +164,9 @@ class FixRigidSmall : public Fix {
   int dilate_group_bit;      // mask for dilation group
   char *id_dilate;           // group name to dilate
 
+  char *id_gravity;         // ID of fix gravity command to add gravity forces
+  double *gvec;             // ptr to gravity vector inside the fix
+
   double p_current[3],p_target[3];
 
   // molecules added on-the-fly as rigid bodies

src/fix_gravity.cpp

@@ -37,7 +37,8 @@ enum{CONSTANT,EQUAL};
 
 FixGravity::FixGravity(LAMMPS *lmp, int narg, char **arg) :
   Fix(lmp, narg, arg),
-  mstr(NULL), vstr(NULL), pstr(NULL), tstr(NULL), xstr(NULL), ystr(NULL), zstr(NULL)
+  mstr(NULL), vstr(NULL), pstr(NULL), tstr(NULL), 
+  xstr(NULL), ystr(NULL), zstr(NULL)
 {
   if (narg < 5) error->all(FLERR,"Illegal fix gravity command");
 
@@ -61,8 +62,10 @@ FixGravity::FixGravity(LAMMPS *lmp, int narg, char **arg) :
     mstyle = CONSTANT;
   }
 
+  int iarg;
+
   if (strcmp(arg[4],"chute") == 0) {
-    if (narg != 6) error->all(FLERR,"Illegal fix gravity command");
+    if (narg < 6) error->all(FLERR,"Illegal fix gravity command");
     style = CHUTE;
     if (strstr(arg[5],"v_") == arg[5]) {
       int n = strlen(&arg[5][2]) + 1;
@@ -73,9 +76,10 @@ FixGravity::FixGravity(LAMMPS *lmp, int narg, char **arg) :
       vert = force->numeric(FLERR,arg[5]);
       vstyle = CONSTANT;
     }
+    iarg = 6;
 
   } else if (strcmp(arg[4],"spherical") == 0) {
-    if (narg != 7) error->all(FLERR,"Illegal fix gravity command");
+    if (narg < 7) error->all(FLERR,"Illegal fix gravity command");
     style = SPHERICAL;
     if (strstr(arg[5],"v_") == arg[5]) {
       int n = strlen(&arg[5][2]) + 1;
@@ -95,9 +99,10 @@ FixGravity::FixGravity(LAMMPS *lmp, int narg, char **arg) :
       theta = force->numeric(FLERR,arg[6]);
       tstyle = CONSTANT;
     }
+    iarg = 7;
 
   } else if (strcmp(arg[4],"vector") == 0) {
-    if (narg != 8) error->all(FLERR,"Illegal fix gravity command");
+    if (narg < 8) error->all(FLERR,"Illegal fix gravity command");
     style = VECTOR;
     if (strstr(arg[5],"v_") == arg[5]) {
       int n = strlen(&arg[5][2]) + 1;
@@ -126,9 +131,23 @@ FixGravity::FixGravity(LAMMPS *lmp, int narg, char **arg) :
       zdir = force->numeric(FLERR,arg[7]);
       zstyle = CONSTANT;
     }
+    iarg = 8;
 
   } else error->all(FLERR,"Illegal fix gravity command");
 
+  // optional keywords
+
+  disable = 0;
+
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"disable") == 0) {
+      disable = 1;
+      iarg++;
+    } else error->all(FLERR,"Illegal fix gravity command");
+  }
+
+  // initializations
+
   degree2rad = MY_PI/180.0;
   time_origin = update->ntimestep;
 
@@ -266,6 +285,12 @@ void FixGravity::post_force(int /*vflag*/)
     set_acceleration();
   }
 
+  // just exit if application of force is disabled
+
+  if (disable) return;
+
+  // apply gravity force to each particle
+
   double **x = atom->x;
   double **f = atom->f;
   double *rmass = atom->rmass;
@@ -338,9 +363,9 @@ void FixGravity::set_acceleration()
     }
   }
 
-  xacc = magnitude*xgrav;
+  gvec[0] = xacc = magnitude*xgrav;
-  yacc = magnitude*ygrav;
+  gvec[1] = yacc = magnitude*ygrav;
-  zacc = magnitude*zgrav;
+  gvec[2] = zacc = magnitude*zgrav;
 }
 
 /* ----------------------------------------------------------------------
@@ -357,3 +382,16 @@ double FixGravity::compute_scalar()
   }
   return egrav_all;
 }
+
+/* ----------------------------------------------------------------------
+   extract current gravity direction vector
+------------------------------------------------------------------------- */
+
+void *FixGravity::extract(const char *name, int &dim)
+{
+  if (strcmp(name,"gvec") == 0) {
+    dim = 1;
+    return (void *) gvec;
+  }
+  return NULL;
+}

src/fix_gravity.h

@@ -36,9 +36,10 @@ class FixGravity : public Fix {
   virtual void post_force(int);
   virtual void post_force_respa(int, int, int);
   double compute_scalar();
+  void *extract(const char *, int &);
 
  protected:
-  int style;
+  int style,disable;
   double magnitude;
   double vert,phi,theta;
   double xdir,ydir,zdir;
@@ -46,6 +47,8 @@ class FixGravity : public Fix {
   double degree2rad;
   int ilevel_respa;
   int time_origin;
+  double gvec[3];
+
   int eflag;
   double egrav,egrav_all;