From bb2846578399375a305318923715fbf80b9cbf9c Mon Sep 17 00:00:00 2001 From: sjplimp Date: Tue, 31 Mar 2015 23:50:13 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@13351 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/RIGID/fix_rigid.cpp | 66 +++++++++------------- src/RIGID/fix_rigid.h | 2 +- src/RIGID/fix_rigid_small.cpp | 103 +++++++++++++++++----------------- src/RIGID/fix_rigid_small.h | 2 +- src/compute_temp.cpp | 4 +- 5 files changed, 85 insertions(+), 92 deletions(-) diff --git a/src/RIGID/fix_rigid.cpp b/src/RIGID/fix_rigid.cpp index 30cd1eb1f8..e21a90d27d 100644 --- a/src/RIGID/fix_rigid.cpp +++ b/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 - // this means users cannot change atom properties like mass between runs - // if they do, bodies will not reflect the changes + // wait to setup bodies until first init() using current atom properties - staticflag = 0; - if (infile) setup_bodies_static(); + setupflag = 0; // 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 this if not using an infile, else was called in constructor + // only do initialization once, b/c properties may not be re-computable + // 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 - called from init() so body/atom properties can be changed between runs - unless reading from infile, in which case called once from constructor + one-time initialization of static rigid body attributes 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; - } + 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 - // staticflag insures this is only done once, not on successive runs + // set rigid body image flags to default values // 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; - } + 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, - since is just setting vcm/angmom at start of run, - which can be estimated value + not done if body properites read from file, e.g. for overlapping particles ------------------------------------------------------------------------- */ void FixRigid::setup_bodies_dynamic() diff --git a/src/RIGID/fix_rigid.h b/src/RIGID/fix_rigid.h index 42ad2c4600..31a29bc80d 100644 --- a/src/RIGID/fix_rigid.h +++ b/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 diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp index e4cdaaca1c..01e1880517 100644 --- a/src/RIGID/fix_rigid_small.cpp +++ b/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 static initialization must be done once - do not do again if using an infile or mol keyword - this is b/c the infile and Molecule class define the COM, etc - so should not reset those values + only do initialization once, b/c properties may not be re-computable + especially if overlapping particles or bodies inserted from mol template + 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 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 (!staticflag || (!infile)) setup_bodies_static(); - //else pre_neighbor(); - - //setup_bodies_dynamic(); - if (!onemols) setup_bodies_dynamic(); + if (!setupflag && !infile) setup_bodies_dynamic(); + setupflag = 1; } /* ---------------------------------------------------------------------- @@ -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 - 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 + one-time initialization of rigid body attributes 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; - } + 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 - // staticflag insures this is only done once, not on successive runs + // set rigid body image flags to default values // 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; - } + 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, - since is just setting initial time=0 vcm and angmom of the body - which can be estimated value + not done if body properites read from file, e.g. for overlapping particles ------------------------------------------------------------------------- */ 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 = 1 to read 6 moments of inertia, store in array - flag inbody = 0 for local bodies whose info is read from file + which = 0 to read everthing except 6 moments of inertia + which = 1 to read just 6 moments of inertia, store in array + 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 (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 { diff --git a/src/RIGID/fix_rigid_small.h b/src/RIGID/fix_rigid_small.h index 5dacdbf9ac..d1305d813c 100644 --- a/src/RIGID/fix_rigid_small.h +++ b/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 diff --git a/src/compute_temp.cpp b/src/compute_temp.cpp index 5a472549ab..0e27aa4f5c 100644 --- a/src/compute_temp.cpp +++ b/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) - error->all(FLERR,"Temperature compute degrees of freedom < 0"); + //if (dof < 0.0 && natoms > 0.0) + // error->all(FLERR,"Temperature compute degrees of freedom < 0"); if (dof > 0.0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; }