diff --git a/doc/src/fix_modify.txt b/doc/src/fix_modify.txt index be3a1669d217550e576031061b3c685111a65718..308bba1ac31fdf9a7cbd45ec1da50122cd46bcb9 100644 --- a/doc/src/fix_modify.txt +++ b/doc/src/fix_modify.txt @@ -23,7 +23,7 @@ keyword = {temp} or {press} or {energy} or {virial} or {respa} or {dynamic/dof} {dynamic/dof} value = {yes} or {no} yes/no = do or do not recompute the number of degrees of freedom (DOF) contributing to the temperature {bodyforces} value = {early} or {late} - early/late = compute per-rigid-body forces and torques at post_force (early) or at final_integrate (late) :pre + early/late = compute rigid-body forces/torques early or late in the timestep :pre :ule [Examples:] @@ -86,9 +86,8 @@ if you want virial contribution of the fix to be part of the relaxation criteria, although this seems unlikely. NOTE: This option is only supported by fixes that explicitly say -so. For some of these (e.g. the -"fix shake"_fix_shake.html command) the default setting is -{virial yes}, for others it is {virial no}. +so. For some of these (e.g. the "fix shake"_fix_shake.html command) +the default setting is {virial yes}, for others it is {virial no}. For fixes that set or modify forces, it may be possible to select at which "r-RESPA"_run_style.html level the fix operates via the {respa} @@ -122,22 +121,27 @@ compute to calculate temperature. See the "compute_modify dynamic/dof"_compute_modify.html command for a similar way to insure correct temperature normalization for those thermostats. -The {bodyforces} keyword determines whether the resultant forces and -torques acting on rigid bodies are computed at the post-force stage of -a time step (i.e. right after per-atom forces have been computed and -communicated among processors) or at the final-integrate stage (i.e. -after all other fixes have finished their post-force tasks). This option -applies for "fix rigid"_fix_rigid.html and "fix rigid/small"_fix_rigid.html, -along with their nve, nvt, npt, and nph versions. It also applies for -"fix poems"_fix_poems.html. Few fix styles actually do post-force -tasks. Some of them cause modifications in the computed per-atoms forces -(e.g. "fix addforce"_fix_addforce.html, "fix setforce"_fix_setforce.html, -"fix spring"_fix_spring.html, "fix shake"_fix_shake.html, and -"fix rattle"_fix_shake.html). These tasks are executed sequentially for -each fix, following the order of their definitions in the input script. -Therefore, once the {bodyforces} keyword is set as {early} for a given -rigid-style fix, per-atom force modifications done by other fixes defined -afterwards will have no effect on the per-body forces/torques it computes. +The {bodyforces} keyword determines whether the forces and torques +acting on rigid bodies are computed {early} at the post-force stage of +each timestep (right after per-atom forces have been computed and +communicated among processors), or {late} at the final-integrate stage +of each timestep (after any other fixes have finished their post-force +tasks). Only the rigid-body integration fixes use this option, which +includes "fix rigid"_fix_rigid.html and "fix +rigid/small"_fix_rigid.html, and their variants, and also "fix +poems"_fix_poems.html. + +The default is {late}. If there are other fixes that add forces to +individual atoms, then the rigid-body constraints will include these +forces when time-integrating the rigid bodies. If {early} is +specified, then new fixes can be written that use or modify the +per-body force and torque, before time-integration of the rigid bodies +occurs. Note however this has the side effect, that fixes such as +"fix addforce"_fix_addforce.html, "fix setforce"_fix_setforce.html, +"fix spring"_fix_spring.html, which add forces to individual atoms +will have no effect on the motion of the rigid bodies if they are +specified in the input script after the fix rigid command. LAMMPS +will give a warning if that is the case. [Restrictions:] none @@ -149,4 +153,5 @@ pressure"_compute_pressure.html, "thermo_style"_thermo_style.html [Default:] The option defaults are temp = ID defined by fix, press = ID defined -by fix, energy = no, virial = varies by fix style, respa = 0, bodyforces = late. \ No newline at end of file +by fix, energy = no, virial = different for each fix style, respa = 0, +bodyforce = late. diff --git a/doc/src/fix_poems.txt b/doc/src/fix_poems.txt index f949acc43357f824e4536de65afcc240410776d3..03abc058b8b0d6e982f6291e1fecbb887cf2ad9a 100644 --- a/doc/src/fix_poems.txt +++ b/doc/src/fix_poems.txt @@ -105,17 +105,19 @@ off, and there is only a single fix poems defined. [Restart, fix_modify, output, run start/stop, minimize info:] -The "fix_modify"_fix_modify.html {bodyforces} option is supported by -this fix style to set whether per-body forces and torques are -computed early or late in a time step, i.e., at the post-force stage -or at the final-integrate stage, respectively. - No information about this fix is written to "binary restart -files"_restart.html. No global or per-atom quantities are stored -by this fix for access by various "output -commands"_Section_howto.html#howto_15. No parameter of this fix can -be used with the {start/stop} keywords of the "run"_run.html command. -This fix is not invoked during "energy minimization"_minimize.html. +files"_restart.html. + +The "fix_modify"_fix_modify.html {bodyforces} option is supported by +this fix style to set whether per-body forces and torques are computed +early or late in a timestep, i.e. at the post-force stage or at the +final-integrate stage, respectively. + +No global or per-atom quantities are stored by this fix for access by +various "output commands"_Section_howto.html#howto_15. No parameter +of this fix can be used with the {start/stop} keywords of the +"run"_run.html command. This fix is not invoked during "energy +minimization"_minimize.html. [Restrictions:] diff --git a/doc/src/fix_rigid.txt b/doc/src/fix_rigid.txt index 9995a77ae24d598846b55bd74a521f82d0c69363..ec7ed4f2b1e16374b0dd318a75cf31b3618711ce 100644 --- a/doc/src/fix_rigid.txt +++ b/doc/src/fix_rigid.txt @@ -223,10 +223,10 @@ via several options. NOTE: With the {rigid/small} styles, which require that {bodystyle} be specified as {molecule} or {custom}, you can define a system that has -no rigid bodies initially. This is useful when you are using the {mol} -keyword in conjunction with another fix that is adding rigid bodies -on-the-fly as molecules, such as "fix deposit"_fix_deposit.html or -"fix pour"_fix_pour.html. +no rigid bodies initially. This is useful when you are using the +{mol} keyword in conjunction with another fix that is adding rigid +bodies on-the-fly as molecules, such as "fix deposit"_fix_deposit.html +or "fix pour"_fix_pour.html. For bodystyle {single} the entire fix group of atoms is treated as one rigid body. This option is only allowed for the {rigid} styles. @@ -744,8 +744,8 @@ instantaneous temperature. The "fix_modify"_fix_modify.html {bodyforces} option is supported by all rigid styles to set whether per-body forces and torques are -computed early or late in a time step, i.e., at the post-force stage -or at the final-integrate stage, respectively. +computed early or late in a timestep, i.e. at the post-force stage or +at the final-integrate stage or the timestep, respectively. The 2 NVE rigid fixes compute a global scalar which can be accessed by various "output commands"_Section_howto.html#howto_15. The scalar diff --git a/src/RIGID/fix_rigid.cpp b/src/RIGID/fix_rigid.cpp index 5d9ed88e311fb93b04a7d09daae9999561db51ec..73a8e487cee29c3400173a14878dd8566a73e222 100644 --- a/src/RIGID/fix_rigid.cpp +++ b/src/RIGID/fix_rigid.cpp @@ -174,7 +174,8 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) : MPI_Allreduce(&vmin,&minval,1,MPI_INT,MPI_MIN,world); molecule = new tagint[nlocal]; for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) molecule[i] = (tagint)((tagint)value[i] - minval + 1); + if (mask[i] & groupbit) + molecule[i] = (tagint)((tagint)value[i] - minval + 1); delete[] value; } else error->all(FLERR,"Unsupported fix rigid custom property"); } else { @@ -624,7 +625,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) : setupflag = 0; - // compute per body forces and torques inside final_integrate() by default + // compute per body forces and torques at final_integrate() by default earlyflag = 0; @@ -688,29 +689,12 @@ 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; mask |= INITIAL_INTEGRATE; mask |= FINAL_INTEGRATE; - mask |= POST_FORCE; + if (langflag || earlyflag) mask |= POST_FORCE; mask |= PRE_NEIGHBOR; mask |= INITIAL_INTEGRATE_RESPA; mask |= FINAL_INTEGRATE_RESPA; @@ -732,19 +716,25 @@ 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 int count = 0; - int myindex, myposition; for (i = 0; i < modify->nfix; i++) - if (strcmp(modify->fix[i]->style,"rigid") == 0) { - count++; - if (strcmp(modify->fix[i]->id,this->id) == 0) { - myindex = i; - myposition = count; + if (modify->fix[i]->rigid_flag) count++; + if (count > 1 && me == 0) error->warning(FLERR,"More than one fix rigid"); + + if (earlyflag) { + int rflag = 0; + for (i = 0; i < modify->nfix; i++) { + if (modify->fix[i]->rigid_flag) rflag = 1; + if (rflag && (modify->fmask[i] & POST_FORCE) && + !modify->fix[i]->rigid_flag) { + char str[128]; + sprintf(str,"Fix %d alters forces after fix rigid",modify->fix[i]->id); + error->warning(FLERR,str); } } - 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 @@ -758,29 +748,6 @@ void FixRigid::init() error->all(FLERR,"Rigid fix must come before NPT/NPH fix"); } - // 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 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); - } - } - } - } - // timestep info dtv = update->dt; @@ -993,7 +960,7 @@ void FixRigid::initial_integrate(int vflag) apply Langevin thermostat to all 6 DOF of rigid bodies computed by proc 0, broadcast to other procs unlike fix langevin, this stores extra force in extra arrays, - which are added in when one calculates a new fcm/torque + which are added in when final_integrate() calculates a new fcm/torque ------------------------------------------------------------------------- */ void FixRigid::apply_langevin_thermostat() @@ -1063,6 +1030,7 @@ void FixRigid::enforce2d() void FixRigid::compute_forces_and_torques() { int i,ibody; + double dtfm; // sum over atoms to get force and torque on rigid body @@ -1125,6 +1093,7 @@ void FixRigid::compute_forces_and_torques() } } + /* ---------------------------------------------------------------------- */ void FixRigid::post_force(int vflag) @@ -2660,6 +2629,21 @@ void FixRigid::zero_rotation() set_v(); } +/* ---------------------------------------------------------------------- */ + +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; +} + /* ---------------------------------------------------------------------- return temperature of collection of rigid bodies non-active DOF are removed by fflag/tflag and in tfactor diff --git a/src/RIGID/fix_rigid.h b/src/RIGID/fix_rigid.h index 0e1936a816e849fbb0514bf9eb19f967eabdea14..507e4c755356686ae3b7a5f3f5c352e49932ae0f 100644 --- a/src/RIGID/fix_rigid.h +++ b/src/RIGID/fix_rigid.h @@ -28,7 +28,6 @@ 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); @@ -55,6 +54,7 @@ class FixRigid : public Fix { void reset_dt(); void zero_momentum(); void zero_rotation(); + virtual int modify_param(int, char **); virtual void *extract(const char*, int &); double extract_ke(); double extract_erotational(); @@ -70,7 +70,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 earlyflag; // 1 if forces/torques computed at post_force() int dimension; // # of dimensions int nbody; // # of rigid bodies @@ -146,7 +146,7 @@ class FixRigid : public Fix { void setup_bodies_static(); void setup_bodies_dynamic(); void apply_langevin_thermostat(); - virtual void compute_forces_and_torques(); + void compute_forces_and_torques(); void readfile(int, double *, double **, double **, double **, imageint *, int *); }; diff --git a/src/RIGID/fix_rigid_nh.cpp b/src/RIGID/fix_rigid_nh.cpp index 672f6db83273d9cf0629abf41db2053126d717d5..96c44d15b5b0e1861803203e4d20e0d0e052e891 100644 --- a/src/RIGID/fix_rigid_nh.cpp +++ b/src/RIGID/fix_rigid_nh.cpp @@ -591,9 +591,9 @@ void FixRigidNH::initial_integrate(int vflag) void FixRigidNH::final_integrate() { - int ibody; + int i,ibody; double tmp,scale_t[3],scale_r; - double dtfm; + double dtfm,xy,xz,yz; double mbody[3],tbody[3],fquat[4]; double dtf2 = dtf * 2.0; @@ -1249,7 +1249,9 @@ int FixRigidNH::modify_param(int narg, char **arg) if (pressure->pressflag == 0) error->all(FLERR,"Fix_modify pressure ID does not compute pressure"); return 2; - } else return FixRigid::modify_param(narg,arg); + } + + return FixRigid::modify_param(narg,arg); } /* ---------------------------------------------------------------------- */ diff --git a/src/RIGID/fix_rigid_nh_small.cpp b/src/RIGID/fix_rigid_nh_small.cpp index c97c1e6d1a557ae59fc6738544daeac1cd1f71b9..135a1fb4bdb0e8ee37b8e958bb31e5b6117a1a33 100644 --- a/src/RIGID/fix_rigid_nh_small.cpp +++ b/src/RIGID/fix_rigid_nh_small.cpp @@ -618,7 +618,7 @@ void FixRigidNHSmall::initial_integrate(int vflag) void FixRigidNHSmall::final_integrate() { - int ibody; + int i,ibody; double tmp,scale_t[3],scale_r; double dtfm; double mbody[3],tbody[3],fquat[4]; @@ -1367,7 +1367,9 @@ 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; - } else return FixRigidSmall::modify_param(narg,arg); + } + + return FixRigidSmall::modify_param(narg,arg); } /* ---------------------------------------------------------------------- */ diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp index a079f30967ac5ccb85734939da6b7295144d67f8..5533098318e30370a88403dbac1a6c9e6d76b2c6 100644 --- a/src/RIGID/fix_rigid_small.cpp +++ b/src/RIGID/fix_rigid_small.cpp @@ -535,29 +535,12 @@ 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; mask |= INITIAL_INTEGRATE; mask |= FINAL_INTEGRATE; - mask |= POST_FORCE; + if (langflag || earlyflag) mask |= POST_FORCE; mask |= PRE_NEIGHBOR; mask |= INITIAL_INTEGRATE_RESPA; mask |= FINAL_INTEGRATE_RESPA; @@ -575,17 +558,22 @@ 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 (strcmp(modify->fix[i]->id,this->id) == 0) { - myindex = i; - myposition = count; + if (modify->fix[i]->rigid_flag) count++; + if (count > 1 && me == 0) error->warning(FLERR,"More than one fix rigid"); + + if (earlyflag) { + int rflag = 0; + for (i = 0; i < modify->nfix; i++) { + if (modify->fix[i]->rigid_flag) rflag = 1; + if (rflag && (modify->fmask[i] & POST_FORCE) && + !modify->fix[i]->rigid_flag) { + char str[128]; + sprintf(str,"Fix %d alters forces after fix rigid",modify->fix[i]->id); + error->warning(FLERR,str); } } - 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 @@ -599,29 +587,6 @@ void FixRigidSmall::init() error->all(FLERR,"Rigid fix must come before NPT/NPH fix"); } - // 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 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); - } - } - } - } - // timestep info dtv = update->dt; @@ -822,10 +787,10 @@ void FixRigidSmall::initial_integrate(int vflag) /* ---------------------------------------------------------------------- apply Langevin thermostat to all 6 DOF of rigid bodies I own unlike fix langevin, this stores extra force in extra arrays, - which are added in when one calculates a new fcm/torque + which are added in when a new fcm/torque are calculated ------------------------------------------------------------------------- */ -void FixRigidSmall::apply_langevin_thermostat() +void FixRigidSmall::apply_langevin_thermostat(int vflag) { double gamma1,gamma2; @@ -984,13 +949,6 @@ void FixRigidSmall::compute_forces_and_torques() } } -/* ---------------------------------------------------------------------- */ - -void FixRigidSmall::post_force(int vflag) -{ - if (langflag) apply_langevin_thermostat(); - if (earlyflag) compute_forces_and_torques(); -} /* ---------------------------------------------------------------------- */ @@ -3454,6 +3412,21 @@ void FixRigidSmall::zero_rotation() /* ---------------------------------------------------------------------- */ +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; +} + +/* ---------------------------------------------------------------------- */ + void *FixRigidSmall::extract(const char *str, int &dim) { if (strcmp(str,"body") == 0) { diff --git a/src/RIGID/fix_rigid_small.h b/src/RIGID/fix_rigid_small.h index 6c4a002059a990005d5cc777903bc54c65ad8fe0..3f6826f9bb140bf131d1ccb2a4c260fee6fd01e5 100644 --- a/src/RIGID/fix_rigid_small.h +++ b/src/RIGID/fix_rigid_small.h @@ -33,7 +33,6 @@ 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); @@ -64,6 +63,7 @@ class FixRigidSmall : public Fix { void reset_dt(); void zero_momentum(); void zero_rotation(); + int modify_param(int, char **); void *extract(const char*, int &); double extract_ke(); double extract_erotational(); @@ -79,7 +79,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 earlyflag; // 1 if forces/torques are computed at post_force() int commflag; // various modes of forward/reverse comm int customflag; // 1 if custom property/variable define bodies int nbody; // total # of rigid bodies @@ -194,7 +194,7 @@ class FixRigidSmall : public Fix { void setup_bodies_static(); void setup_bodies_dynamic(); void apply_langevin_thermostat(); - virtual void compute_forces_and_torques(); + void compute_forces_and_torques(); void readfile(int, double **, int *); void grow_body(); void reset_atom2body(); diff --git a/src/USER-OMP/fix_rigid_omp.h b/src/USER-OMP/fix_rigid_omp.h index 372486aa859b199a027e61240b32e6be2e64ba1c..b666a173a87c0a949ba8e8aba33f81da493f345a 100644 --- a/src/USER-OMP/fix_rigid_omp.h +++ b/src/USER-OMP/fix_rigid_omp.h @@ -45,4 +45,3 @@ class FixRigidOMP : public FixRigid { #endif #endif - diff --git a/src/USER-OMP/fix_rigid_small_omp.cpp b/src/USER-OMP/fix_rigid_small_omp.cpp index acf4b7919f22cca685106cc243cc1468cb537caa..a0495b3b933a4fa5a852b079dfd5936cb8b10088 100644 --- a/src/USER-OMP/fix_rigid_small_omp.cpp +++ b/src/USER-OMP/fix_rigid_small_omp.cpp @@ -618,3 +618,4 @@ void FixRigidSmallOMP::set_v_thr() } } } +