diff --git a/src/POEMS/fix_poems.cpp b/src/POEMS/fix_poems.cpp index b8fdb03447..4ed43bb55c 100644 --- a/src/POEMS/fix_poems.cpp +++ b/src/POEMS/fix_poems.cpp @@ -280,6 +280,10 @@ FixPOEMS::FixPOEMS(LAMMPS *lmp, int narg, char **arg) : fprintf(logfile,"%d clusters, %d bodies, %d joints, %d atoms\n", ncluster,nbody,njoint,nsum); } + + // compute per body forces and torques inside final_integrate() by default + + earlyflag = 0; } /* ---------------------------------------------------------------------- @@ -330,6 +334,25 @@ FixPOEMS::~FixPOEMS() /* ---------------------------------------------------------------------- */ +int FixPOEMS::modify_param(int narg, char **arg) +{ + if (strcmp(arg[0],"bodyforces") == 0) { + if (narg < 2) + error->all(FLERR,"Illegal fix_modify command"); + if (strcmp(arg[1],"early") == 0) + earlyflag = 1; + else if (strcmp(arg[1],"late") == 0) + earlyflag = 0; + else + error->all(FLERR,"Illegal fix_modify command"); + return 2; + } + else + return 0; +} + +/* ---------------------------------------------------------------------- */ + int FixPOEMS::setmask() { int mask = 0; @@ -352,9 +375,17 @@ void FixPOEMS::init() // warn if more than one POEMS fix int count = 0; - for (int i = 0; i < modify->nfix; i++) - if (strcmp(modify->fix[i]->style,"poems") == 0) count++; - if (count > 1 && comm->me == 0) error->warning(FLERR,"More than one fix poems"); + int myindex, myposition; + for (i = 0; i < modify->nfix; i++) + if (strcmp(modify->fix[i]->style,"poems") == 0) { + count++; + if (strcmp(modify->fix[i]->id,id) != 0) { + myindex = i; + myposition = count; + } + } + if (count > 1 && myposition == 1 && comm->me == 0) + error->warning(FLERR,"More than one fix poems"); // error if npt,nph fix comes before rigid fix @@ -368,21 +399,26 @@ void FixPOEMS::init() error->all(FLERR,"POEMS fix must come before NPT/NPH fix"); } - // warn if any non-rigid fix with post_force succeeds any fix rigid: - int first_rigid = 0; - while (!modify->fix[first_rigid]->rigid_flag) - first_rigid++; - count = 0; - for (i = first_rigid + 1; i < modify->nfix; i++) { - Fix *ifix = modify->fix[i]; - if ( (modify->fmask[i] & POST_FORCE) && (!ifix->rigid_flag) ) { - count++; - if (comm->me == 0) { - if (count == 1) - error->warning(FLERR,"The following fixes must preceed any rigid-body fixes"); - if (screen) fprintf(screen,"> fix %s %s\n",ifix->id,ifix->style); - if (logfile) fprintf(logfile,"> fix %s %s\n",ifix->id,ifix->style); - } + // warn if fix poems preceeds non-rigid fixes with post-force tasks + // when computing body forces and torques in post_force() instead of final_integrate() + + if (earlyflag) { + int has_post_force[modify->nfix - myindex]; + count = 0; + for (i = myindex + 1; i < modify->nfix; i++) + if ( (modify->fmask[i] & POST_FORCE) && (!modify->fix[i]->rigid_flag) ) + has_post_force[count++] = i; + if (count) { + FILE *p[2] = {screen, logfile}; + for (int j = 0; j < 2; j++) + if (p[j]) { + fprintf(p[j],"WARNING: fix %s %s",id,style); + fprintf(p[j]," will add up forces before they are handled by:\n"); + for (int k = 0; k < count; k++) { + Fix *fix = modify->fix[has_post_force[k]]; + fprintf(p[j]," => fix %s %s\n",fix->id,fix->style); + } + } } } @@ -747,7 +783,7 @@ void FixPOEMS::initial_integrate(int vflag) only count joint atoms in 1st body ------------------------------------------------------------------------- */ -void FixPOEMS::post_force(int vflag) +void FixPOEMS::compute_forces_and_torques() { int i,ibody; int xbox,ybox,zbox; @@ -798,6 +834,13 @@ void FixPOEMS::post_force(int vflag) } } +/* ---------------------------------------------------------------------- */ + +void FixPOEMS::post_force(int vflag) +{ + if (earlyflag) compute_forces_and_torques(); +} + /* ---------------------------------------------------------------------- update vcm,omega by last 1/2 step set v of body atoms accordingly @@ -805,6 +848,8 @@ void FixPOEMS::post_force(int vflag) void FixPOEMS::final_integrate() { + if (!earlyflag) compute_forces_and_torques(); + // perform POEMS integration poems->LobattoTwo(vcm,omega,torque,fcm); diff --git a/src/POEMS/fix_poems.h b/src/POEMS/fix_poems.h index 4d5e5c3198..eec146a8dd 100644 --- a/src/POEMS/fix_poems.h +++ b/src/POEMS/fix_poems.h @@ -28,6 +28,7 @@ class FixPOEMS : public Fix { public: FixPOEMS(class LAMMPS *, int narg, char **arg); ~FixPOEMS(); + int modify_param(int, char **); int setmask(); void init(); void setup(int); @@ -56,6 +57,8 @@ class FixPOEMS : public Fix { int nlevels_respa; double total_ke; + int earlyflag; // 1 if forces and torques are computed at post_force() + // atom assignment to rigid bodies // double count joint atoms as being in multiple bodies @@ -96,6 +99,7 @@ class FixPOEMS : public Fix { // internal class functions + void compute_forces_and_torques(); void readfile(char *); int readline(FILE *, char **, int *); void jointbuild(); diff --git a/src/RIGID/fix_rigid.cpp b/src/RIGID/fix_rigid.cpp index a801c2225f..47079b4188 100644 --- a/src/RIGID/fix_rigid.cpp +++ b/src/RIGID/fix_rigid.cpp @@ -639,6 +639,25 @@ FixRigid::~FixRigid() /* ---------------------------------------------------------------------- */ +int FixRigid::modify_param(int narg, char **arg) +{ + if (strcmp(arg[0],"bodyforces") == 0) { + if (narg < 2) + error->all(FLERR,"Illegal fix_modify command"); + if (strcmp(arg[1],"early") == 0) + earlyflag = 1; + else if (strcmp(arg[1],"late") == 0) + earlyflag = 0; + else + error->all(FLERR,"Illegal fix_modify command"); + return 2; + } + else + return 0; +} + +/* ---------------------------------------------------------------------- */ + int FixRigid::setmask() { int mask = 0; @@ -668,9 +687,17 @@ void FixRigid::init() // warn if more than one rigid fix int count = 0; + int myindex, myposition; for (i = 0; i < modify->nfix; i++) - if (strcmp(modify->fix[i]->style,"rigid") == 0) count++; - if (count > 1 && me == 0) error->warning(FLERR,"More than one fix rigid"); + if (strcmp(modify->fix[i]->style,"rigid") == 0) { + count++; + if (strcmp(modify->fix[i]->id,id) != 0) { + myindex = i; + myposition = count; + } + } + if (count > 1 && myposition == 1 && comm->me == 0) + error->warning(FLERR,"More than one fix rigid"); // error if npt,nph fix comes before rigid fix @@ -684,25 +711,26 @@ void FixRigid::init() error->all(FLERR,"Rigid fix must come before NPT/NPH fix"); } - // warn if any non-rigid fix with post_force succeeds any fix rigid instance + // warn if fix rigid preceeds non-rigid fixes with post-force tasks // when computing body forces and torques in post_force() instead of final_integrate() if (earlyflag) { - int first_rigid = 0; - while (!modify->fix[first_rigid]->rigid_flag) - first_rigid++; + int has_post_force[modify->nfix - myindex]; count = 0; - for (i = first_rigid + 1; i < modify->nfix; i++) { - Fix *ifix = modify->fix[i]; - if ( (modify->fmask[i] & POST_FORCE) && (!ifix->rigid_flag) ) { - count++; - if (comm->me == 0) { - if (count == 1) - error->warning(FLERR,"The following fixes must preceed any rigid-body fixes"); - if (screen) fprintf(screen,"> fix %s %s\n",ifix->id,ifix->style); - if (logfile) fprintf(logfile,"> fix %s %s\n",ifix->id,ifix->style); + for (i = myindex + 1; i < modify->nfix; i++) + if ( (modify->fmask[i] & POST_FORCE) && (!modify->fix[i]->rigid_flag) ) + has_post_force[count++] = i; + if (count) { + FILE *p[2] = {screen, logfile}; + for (int j = 0; j < 2; j++) + if (p[j]) { + fprintf(p[j],"WARNING: fix %s %s",id,style); + fprintf(p[j]," will add up forces before they are handled by:\n"); + for (int k = 0; k < count; k++) { + Fix *fix = modify->fix[has_post_force[k]]; + fprintf(p[j]," => fix %s %s\n",fix->id,fix->style); + } } - } } } @@ -1054,8 +1082,10 @@ void FixRigid::compute_forces_and_torques() void FixRigid::post_force(int vflag) { - if (langflag) apply_langevin_thermostat(); - if (earlyflag) compute_forces_and_torques(); + if (earlyflag) { + if (langflag) apply_langevin_thermostat(); + compute_forces_and_torques(); + } } /* ---------------------------------------------------------------------- */ @@ -1065,7 +1095,10 @@ void FixRigid::final_integrate() int ibody; double dtfm; - if (!earlyflag) compute_forces_and_torques(); + if (!earlyflag) { + if (langflag) apply_langevin_thermostat(); + compute_forces_and_torques(); + } // update vcm and angmom // fflag,tflag = 0 for some dimensions in 2d @@ -2701,18 +2734,3 @@ double FixRigid::compute_array(int i, int j) if (j == 13) return (imagebody[i] >> IMGBITS & IMGMASK) - IMGMAX; return (imagebody[i] >> IMG2BITS) - IMGMAX; } - -/* ---------------------------------------------------------------------- */ - -int FixRigid::modify_param(int narg, char **arg) -{ - if (strcmp(arg[0],"bodyforces") == 0) { - if (narg < 2) error->all(FLERR,"Illegal fix_modify command"); - if (strcmp(arg[1],"early") == 0) earlyflag=1; - else if (strcmp(arg[1],"late") == 0) earlyflag=0; - else error->all(FLERR,"Illegal fix_modify command"); - - return 2; - } - return 0; -} diff --git a/src/RIGID/fix_rigid.h b/src/RIGID/fix_rigid.h index 6c2f18b498..6e0f67faaf 100644 --- a/src/RIGID/fix_rigid.h +++ b/src/RIGID/fix_rigid.h @@ -28,6 +28,7 @@ class FixRigid : public Fix { public: FixRigid(class LAMMPS *, int, char **); virtual ~FixRigid(); + int modify_param(int, char **); virtual int setmask(); virtual void init(); virtual void setup(int); @@ -38,7 +39,6 @@ class FixRigid : public Fix { void final_integrate_respa(int, int); void write_restart_file(char *); virtual double compute_scalar(); - virtual int modify_param(int, char **); double memory_usage(); void grow_arrays(int); @@ -72,6 +72,7 @@ class FixRigid : public Fix { char *infile; // 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 and torques are computed at post_force() int dimension; // # of dimensions int nbody; // # of rigid bodies @@ -107,7 +108,6 @@ class FixRigid : public Fix { int orientflag; // 1 if particles store spatial orientation int dorientflag; // 1 if particles store dipole orientation int reinitflag; // 1 if re-initialize rigid bodies between runs - int earlyflag; // 1 if compute body forces and torques early, i.e. in post_force() imageint *xcmimage; // internal image flags for atoms in rigid bodies // set relative to in-box xcm of each body diff --git a/src/RIGID/fix_rigid_nh_small.cpp b/src/RIGID/fix_rigid_nh_small.cpp index 199c04bd9c..2b55c73ec6 100644 --- a/src/RIGID/fix_rigid_nh_small.cpp +++ b/src/RIGID/fix_rigid_nh_small.cpp @@ -1438,9 +1438,8 @@ int FixRigidNHSmall::modify_param(int narg, char **arg) if (pressure->pressflag == 0) error->all(FLERR,"Fix_modify pressure ID does not compute pressure"); return 2; - } - - return 0; + } else + return FixRigidSmall::modify_param(narg,arg); } /* ---------------------------------------------------------------------- */ diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp index 317d0e36b4..5fb7e7a062 100644 --- a/src/RIGID/fix_rigid_small.cpp +++ b/src/RIGID/fix_rigid_small.cpp @@ -481,6 +481,25 @@ FixRigidSmall::~FixRigidSmall() /* ---------------------------------------------------------------------- */ +int FixRigidSmall::modify_param(int narg, char **arg) +{ + if (strcmp(arg[0],"bodyforces") == 0) { + if (narg < 2) + error->all(FLERR,"Illegal fix_modify command"); + if (strcmp(arg[1],"early") == 0) + earlyflag = 1; + else if (strcmp(arg[1],"late") == 0) + earlyflag = 0; + else + error->all(FLERR,"Illegal fix_modify command"); + return 2; + } + else + return 0; +} + +/* ---------------------------------------------------------------------- */ + int FixRigidSmall::setmask() { int mask = 0; @@ -504,9 +523,17 @@ void FixRigidSmall::init() // warn if more than one rigid fix int count = 0; + int myindex, myposition; for (i = 0; i < modify->nfix; i++) - if (strcmp(modify->fix[i]->style,"rigid") == 0) count++; - if (count > 1 && me == 0) error->warning(FLERR,"More than one fix rigid"); + if (strcmp(modify->fix[i]->style,"rigid") == 0) { + count++; + if (strcmp(modify->fix[i]->id,id) != 0) { + myindex = i; + myposition = count; + } + } + if (count > 1 && myposition == 1 && comm->me == 0) + error->warning(FLERR,"More than one fix rigid"); // error if npt,nph fix comes before rigid fix @@ -520,27 +547,28 @@ void FixRigidSmall::init() error->all(FLERR,"Rigid fix must come before NPT/NPH fix"); } - // warn if any non-rigid fix with post_force succeeds any fix rigid instance + // warn if fix rigid preceeds non-rigid fixes with post-force tasks // when computing body forces and torques in post_force() instead of final_integrate() if (earlyflag) { - int first_rigid = 0; - while (!modify->fix[first_rigid]->rigid_flag) - first_rigid++; + int has_post_force[modify->nfix - myindex]; count = 0; - for (i = first_rigid + 1; i < modify->nfix; i++) { - Fix *ifix = modify->fix[i]; - if ( (modify->fmask[i] & POST_FORCE) && (!ifix->rigid_flag) ) { - count++; - if (comm->me == 0) { - if (count == 1) - error->warning(FLERR,"The following fixes must preceed any rigid-body fixes"); - if (screen) fprintf(screen,"> fix %s %s\n",ifix->id,ifix->style); - if (logfile) fprintf(logfile,"> fix %s %s\n",ifix->id,ifix->style); + for (i = myindex + 1; i < modify->nfix; i++) + if ( (modify->fmask[i] & POST_FORCE) && (!modify->fix[i]->rigid_flag) ) + has_post_force[count++] = i; + if (count) { + FILE *p[2] = {screen, logfile}; + for (int j = 0; j < 2; j++) + if (p[j]) { + fprintf(p[j],"WARNING: fix %s %s",id,style); + fprintf(p[j]," will add up forces before they are handled by:\n"); + for (int k = 0; k < count; k++) { + Fix *fix = modify->fix[has_post_force[k]]; + fprintf(p[j]," => fix %s %s\n",fix->id,fix->style); + } } - } - } - } + } + } // timestep info @@ -908,8 +936,10 @@ void FixRigidSmall::compute_forces_and_torques() void FixRigidSmall::post_force(int vflag) { - if (langflag) apply_langevin_thermostat(); - if (earlyflag) compute_forces_and_torques(); + if (earlyflag) { + if (langflag) apply_langevin_thermostat(); + compute_forces_and_torques(); + } } /* ---------------------------------------------------------------------- */ @@ -918,9 +948,12 @@ void FixRigidSmall::final_integrate() { double dtfm; - //check(3) + //check(3); - if (!earlyflag) compute_forces_and_torques(); + if (!earlyflag) { + if (langflag) apply_langevin_thermostat(); + compute_forces_and_torques(); + } // update vcm and angmom, recompute omega @@ -3510,21 +3543,6 @@ double FixRigidSmall::compute_scalar() return tall; } -/* ---------------------------------------------------------------------- */ - -int FixRigidSmall::modify_param(int narg, char **arg) -{ - if (strcmp(arg[0],"bodyforces") == 0) { - if (narg < 2) error->all(FLERR,"Illegal fix_modify command"); - if (strcmp(arg[1],"early") == 0) earlyflag=1; - else if (strcmp(arg[1],"late") == 0) earlyflag=0; - else error->all(FLERR,"Illegal fix_modify command"); - - return 2; - } - return 0; -} - /* ---------------------------------------------------------------------- memory usage of local atom-based arrays ------------------------------------------------------------------------- */ diff --git a/src/RIGID/fix_rigid_small.h b/src/RIGID/fix_rigid_small.h index 37fb27b7cf..71a0af729a 100644 --- a/src/RIGID/fix_rigid_small.h +++ b/src/RIGID/fix_rigid_small.h @@ -33,11 +33,12 @@ class FixRigidSmall : public Fix { public: FixRigidSmall(class LAMMPS *, int, char **); virtual ~FixRigidSmall(); + int modify_param(int, char **); virtual int setmask(); virtual void init(); virtual void setup(int); virtual void initial_integrate(int); - virtual void post_force(int); + void post_force(int); virtual void final_integrate(); void initial_integrate_respa(int, int, int); void final_integrate_respa(int, int); @@ -69,7 +70,6 @@ class FixRigidSmall : public Fix { double extract_ke(); double extract_erotational(); double compute_scalar(); - virtual int modify_param(int, char **); double memory_usage(); protected: @@ -81,6 +81,7 @@ class FixRigidSmall : public Fix { char *infile; // file to read rigid body attributes from int setupflag; // 1 if body properties are setup, else 0 + int earlyflag; // 1 if forces and torques are computed at post_force() int commflag; // various modes of forward/reverse comm int nbody; // total # of rigid bodies int nlinear; // total # of linear rigid bodies @@ -131,7 +132,6 @@ class FixRigidSmall : public Fix { int orientflag; // 1 if particles store spatial orientation int dorientflag; // 1 if particles store dipole orientation int reinitflag; // 1 if re-initialize rigid bodies between runs - int earlyflag; // 1 if compute body forces and torques early, i.e. in post_force() int POINT,SPHERE,ELLIPSOID,LINE,TRIANGLE,DIPOLE; // bitmasks for eflags int OMEGA,ANGMOM,TORQUE;