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

src/RIGID/fix_rigid.cpp

@@ -548,12 +548,9 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
   MINUSPI = -MY_PI;
   TWOPI = 2.0*MY_PI;
 
-  // if infile set, only setup rigid bodies once, using info from file
+  // wait to setup bodies until first init() using current atom properties
-  // this means users cannot change atom properties like mass between runs
-  // if they do, bodies will not reflect the changes
 
-  staticflag = 0;
+  setupflag = 0;
-  if (infile) setup_bodies_static();
 
   // print statistics
 
@@ -670,10 +667,17 @@ void FixRigid::init()
     step_respa = ((Respa *) update->integrate)->step;
 
   // setup rigid bodies, using current atom info
-  // allows resetting of atom properties like mass between runs
+  // only do initialization once, b/c properties may not be re-computable
-  // only do this if not using an infile, else was called in constructor
+  //   especially if overlapping particles
+  // do not do dynamic init if read body properties from infile
+  //   this is b/c the infile defines the static and dynamic properties
+  //   and may not be computable if contain overlapping particles
+  //   setup_bodies_static() reads infile itself
 
-  if (!infile) setup_bodies_static();
+  if (!setupflag) {
+    setup_bodies_static();
+    if (!infile) setup_bodies_dynamic();
+  }
 
   // temperature scale factor
 
@@ -709,7 +713,6 @@ void FixRigid::setup(int vflag)
   // so angmom_to_omega() and set_v() below will set per-atom vels correctly
   // re-calling it every run allows reset of body/atom velocities between runs
 
-  setup_bodies_dynamic();
 
   // fcm = force on center-of-mass of each rigid body
 
@@ -1085,7 +1088,7 @@ int FixRigid::dof(int tgroup)
 {
   // cannot count DOF correctly unless setup_bodies_static() has been called
 
-  if (!staticflag) {
+  if (!setupflag) {
     if (comm->me == 0) 
       error->warning(FLERR,"Cannot count rigid body degrees-of-freedom "
                      "before bodies are initialized");
@@ -1494,12 +1497,11 @@ void FixRigid::set_v()
 }
 
 /* ----------------------------------------------------------------------
-   initialization of static rigid body attributes
+   one-time initialization of static rigid body attributes
-   called from init() so body/atom properties can be changed between runs
-   unless reading from infile, in which case called once from constructor
    sets extended flags, masstotal, center-of-mass
    sets Cartesian and diagonalized inertia tensor
-   sets body image flags, but only on first call
+   sets body image flags
+   may read some properties from infile
 ------------------------------------------------------------------------- */
 
 void FixRigid::setup_bodies_static()
@@ -1585,14 +1587,12 @@ void FixRigid::setup_bodies_static()
     }
   }
 
-  // first-time setting of body xcmimage flags = true image flags
+  // set body xcmimage flags = true image flags
 
-  if (!staticflag) {
   imageint *image = atom->image;
   for (i = 0; i < nlocal; i++)
     if (body[i] >= 0) xcmimage[i] = image[i];
     else xcmimage[i] = 0;
-  }
 
   // compute masstotal & center-of-mass of each rigid body
   // error if image flag is not 0 in a non-periodic dim
@@ -1662,16 +1662,13 @@ void FixRigid::setup_bodies_static()
     readfile(0,masstotal,xcm,inbody);
   }
 
-  // one-time set of rigid body image flags to default values
+  // set rigid body image flags to default values
-  //   staticflag insures this is only done once, not on successive runs
   // then remap the xcm of each body back into simulation box
-  //   and reset body xcmimage flags via pre_neighbor()
+  //   and reset body and atom xcmimage flags via pre_neighbor()
 
-  if (!staticflag) {
   for (ibody = 0; ibody < nbody; ibody++)
     imagebody[ibody] = ((imageint) IMGMAX << IMG2BITS) | 
       ((imageint) IMGMAX << IMGBITS) | IMGMAX;
-  }
 
   pre_neighbor();
 
@@ -2004,19 +2001,12 @@ void FixRigid::setup_bodies_static()
   }
 
   if (infile) memory->destroy(inbody);
-
-  // static properties have now been initialized once
-  // used to prevent re-initialization which would re-read infile
-
-  staticflag = 1;
 }
 
 /* ----------------------------------------------------------------------
-   initialization of dynamic rigid body attributes
+   one-time initialization of dynamic rigid body attributes
    set vcm and angmom, computed explicitly from constituent particles
-   OK if wrong for overlapping particles,
+   not done if body properites read from file, e.g. for overlapping particles
-     since is just setting vcm/angmom at start of run,
-     which can be estimated value
 ------------------------------------------------------------------------- */
 
 void FixRigid::setup_bodies_dynamic()

src/RIGID/fix_rigid.h

@@ -68,7 +68,7 @@ class FixRigid : public Fix {
 
   char *infile;             // file to read rigid body attributes from
   int rstyle;               // SINGLE,MOLECULE,GROUP
-  int staticflag;           // 1 if static body properties are setup, else 0
+  int setupflag;            // 1 if body properties are setup, else 0
 
   int dimension;            // # of dimensions
   int nbody;                // # of rigid bodies

src/RIGID/fix_rigid_small.cpp

@@ -45,9 +45,9 @@ using namespace MathConst;
 
 FixRigidSmall *FixRigidSmall::frsptr;
 
-#define MAXLINE 256
+#define MAXLINE 1024
 #define CHUNK 1024
-#define ATTRIBUTE_PERBODY 11
+#define ATTRIBUTE_PERBODY 17
 
 #define TOLERANCE 1.0e-6
 #define EPSILON 1.0e-7
@@ -178,6 +178,7 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
                    "fix rigid/small does not exist");
       onemols = &atom->molecules[imol];
       nmol = onemols[0]->nset;
+      restart_file = 1;
       iarg += 2;
 
     } else if (strcmp(arg[iarg],"temp") == 0) {
@@ -423,7 +424,9 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
   mass_body = NULL;
   nmax_mass = 0;
 
-  staticflag = 0;
+  // wait to setup bodies until comm stencils are defined
+
+  setupflag = 0;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -507,11 +510,12 @@ void FixRigidSmall::init()
 
 /* ----------------------------------------------------------------------
    setup static/dynamic properties of rigid bodies, using current atom info
-   allows resetting of atom properties like mass between runs
+   only do initialization once, b/c properties may not be re-computable
-   only do static initialization must be done once
+     especially if overlapping particles or bodies inserted from mol template
-     do not do again if using an infile or mol keyword
+   do not do dynamic init if read body properties from infile
-     this is b/c the infile and Molecule class define the COM, etc
+     this is b/c the infile defines the static and dynamic properties
-       so should not reset those values
+     and may not be computable if contain overlapping particles
+     setup_bodies_static() reads infile itself
    cannot do this until now, b/c requires comm->setup() to have setup stencil
    invoke pre_neighbor() to insure body xcmimage flags are reset
      needed if Verlet::setup::pbc() has remapped/migrated atoms for 2nd run
@@ -520,14 +524,10 @@ void FixRigidSmall::init()
 
 void FixRigidSmall::setup_pre_neighbor()
 {
-  if (!staticflag || (!infile && !onemols)) setup_bodies_static();
+  if (!setupflag) setup_bodies_static();
   else pre_neighbor();
-
+  if (!setupflag && !infile) setup_bodies_dynamic();
-  //if (!staticflag || (!infile)) setup_bodies_static();
+  setupflag = 1;
-  //else pre_neighbor();
-
-  //setup_bodies_dynamic();
-  if (!onemols) setup_bodies_dynamic();
 }
 
 /* ----------------------------------------------------------------------
@@ -961,7 +961,7 @@ int FixRigidSmall::dof(int tgroup)
 
   // cannot count DOF correctly unless setup_bodies_static() has been called
 
-  if (!staticflag) {
+  if (!setupflag) {
     if (comm->me == 0) 
       error->warning(FLERR,"Cannot count rigid body degrees-of-freedom "
                      "before bodies are fully initialized");
@@ -1659,13 +1659,11 @@ void FixRigidSmall::ring_farthest(int n, char *cbuf)
 }
 
 /* ----------------------------------------------------------------------
-   initialization of rigid body attributes
+   one-time initialization of rigid body attributes
-   called at setup, so body/atom properties can be changed between runs
-   unless reading from infile or defined by Molecule class,
-     in which case only called once
    sets extended flags, masstotal, center-of-mass
    sets Cartesian and diagonalized inertia tensor
-   sets body image flags, but only on first call
+   sets body image flags
+   may read some properties from infile
 ------------------------------------------------------------------------- */
 
 void FixRigidSmall::setup_bodies_static()
@@ -1764,14 +1762,12 @@ void FixRigidSmall::setup_bodies_static()
     }
   }
 
-  // first-time setting of body xcmimage flags = true image flags
+  // set body xcmimage flags = true image flags
 
-  if (!staticflag) {
   imageint *image = atom->image;
   for (i = 0; i < nlocal; i++)
     if (bodytag[i] >= 0) xcmimage[i] = image[i];
     else xcmimage[i] = 0;
-  }
 
   // acquire ghost bodies via forward comm
   // set atom2body for ghost atoms via forward comm
@@ -1834,16 +1830,13 @@ void FixRigidSmall::setup_bodies_static()
     readfile(0,NULL,inbody);
   }
 
-  // one-time set of rigid body image flags to default values
+  // set rigid body image flags to default values
-  //   staticflag insures this is only done once, not on successive runs
   // then remap the xcm of each body back into simulation box
-  //   and reset body xcmimage flags via pre_neighbor()
+  //   and reset body and atom xcmimage flags via pre_neighbor()
 
-  if (!staticflag) {
   for (ibody = 0; ibody < nlocal_body; ibody++)
     body[ibody].image = ((imageint) IMGMAX << IMG2BITS) | 
       ((imageint) IMGMAX << IMGBITS) | IMGMAX;
-  }
 
   pre_neighbor();
 
@@ -2177,19 +2170,12 @@ void FixRigidSmall::setup_bodies_static()
 
   memory->destroy(itensor);
   if (infile) memory->destroy(inbody);
-
-  // static properties have now been initialized once
-  // used to prevent re-initialization which would re-read infile
-
-  staticflag = 1;
 }
 
 /* ----------------------------------------------------------------------
    one-time initialization of dynamic rigid body attributes
    vcm and angmom, computed explicitly from constituent particles
-   even if wrong for overlapping particles, is OK,
+   not done if body properites read from file, e.g. for overlapping particles
-     since is just setting initial time=0 vcm and angmom of the body
-     which can be estimated value
 ------------------------------------------------------------------------- */
 
 void FixRigidSmall::setup_bodies_dynamic()
@@ -2295,9 +2281,9 @@ void FixRigidSmall::setup_bodies_dynamic()
 
 /* ----------------------------------------------------------------------
    read per rigid body info from user-provided file
-   which = 0 to read total mass and center-of-mass
+   which = 0 to read everthing except 6 moments of inertia
-   which = 1 to read 6 moments of inertia, store in array
+   which = 1 to read just 6 moments of inertia, store in array
-   flag inbody = 0 for local bodies whose info is read from file
+   flag inbody = 0 for local bodies this proc initializes from file
    nlines = # of lines of rigid body info, 0 is OK
    one line = rigid-ID mass xcm ycm zcm ixx iyy izz ixy ixz iyz
    and rigid-ID = mol-ID for fix rigid/small
@@ -2390,6 +2376,12 @@ void FixRigidSmall::readfile(int which, double **array, int *inbody)
         body[m].xcm[0] = atof(values[2]);
         body[m].xcm[1] = atof(values[3]);
         body[m].xcm[2] = atof(values[4]);
+        body[m].vcm[0] = atof(values[11]);
+        body[m].vcm[1] = atof(values[12]);
+        body[m].vcm[2] = atof(values[13]);
+        body[m].angmom[0] = atof(values[14]);
+        body[m].angmom[1] = atof(values[15]);
+        body[m].angmom[2] = atof(values[16]);
       } else {
         array[m][0] = atof(values[5]);
         array[m][1] = atof(values[6]);
@@ -2424,7 +2416,7 @@ void FixRigidSmall::write_restart_file(char *file)
 
   // do not write file if bodies have not yet been intialized
 
-  if (!staticflag) return;
+  if (!setupflag) return;
 
   // proc 0 opens file and writes header
 
@@ -2446,8 +2438,9 @@ void FixRigidSmall::write_restart_file(char *file)
 
   // communication buffer for all my rigid body info
   // max_size = largest buffer needed by any proc
+  // ncol = # of values per line in output file
 
-  int ncol = 11;
+  int ncol = ATTRIBUTE_PERBODY;
   int sendrow = nlocal_body;
   int maxrow;
   MPI_Allreduce(&sendrow,&maxrow,1,MPI_INT,MPI_MAX,world);
@@ -2479,6 +2472,13 @@ void FixRigidSmall::write_restart_file(char *file)
     buf[i][8] = ispace[0][1];
     buf[i][9] = ispace[0][2];
     buf[i][10] = ispace[1][2];
+    buf[i][10] = ispace[1][2];
+    buf[i][11] = body[i].vcm[0];
+    buf[i][12] = body[i].vcm[1];
+    buf[i][13] = body[i].vcm[2];
+    buf[i][14] = body[i].angmom[0];
+    buf[i][15] = body[i].angmom[1];
+    buf[i][16] = body[i].angmom[2];
   }
 
   // write one chunk of rigid body info per proc to file
@@ -2501,10 +2501,13 @@ void FixRigidSmall::write_restart_file(char *file)
 
       for (int i = 0; i < recvrow; i++)
         fprintf(fp,"%d %-1.16e %-1.16e %-1.16e %-1.16e "
+                "%-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e "
                 "%-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e\n",
                 static_cast<int> (buf[i][0]),buf[i][1],
                 buf[i][2],buf[i][3],buf[i][4],
-                buf[i][5],buf[i][6],buf[i][7],buf[i][8],buf[i][9],buf[i][10]);
+                buf[i][5],buf[i][6],buf[i][7],buf[i][8],buf[i][9],buf[i][10],
+                buf[i][11],buf[i][12],buf[i][13],
+                buf[i][14],buf[i][15],buf[i][16]);
     }
     
   } else {

src/RIGID/fix_rigid_small.h

@@ -78,7 +78,7 @@ class FixRigidSmall : public Fix {
   double MINUSPI,TWOPI;
 
   char *infile;             // file to read rigid body attributes from
-  int staticflag;           // 1 if static body properties are setup, else 0
+  int setupflag;            // 1 if body properties are setup, else 0
   int commflag;             // various modes of forward/reverse comm
   int nbody;                // total # of rigid bodies
   tagint maxmol;            // max mol-ID

src/compute_temp.cpp

@@ -64,8 +64,8 @@ void ComputeTemp::dof_compute()
   double natoms = group->count(igroup);
   dof = domain->dimension * natoms;
   dof -= extra_dof + fix_dof;
-  if (dof < 0.0 && natoms > 0.0) 
+  //if (dof < 0.0 && natoms > 0.0) 
-    error->all(FLERR,"Temperature compute degrees of freedom < 0");
+  //  error->all(FLERR,"Temperature compute degrees of freedom < 0");
   if (dof > 0.0) tfactor = force->mvv2e / (dof * force->boltz);
   else tfactor = 0.0;
 }