diff --git a/doc/src/fix_rigid.txt b/doc/src/fix_rigid.txt
index 03edf61ed8fd997370513f209bf4fc417bace630..dbadd3fa634a173b15485554a7f7a79bd3a79e91 100644
--- a/doc/src/fix_rigid.txt
+++ b/doc/src/fix_rigid.txt
@@ -31,11 +31,12 @@ bodystyle = {single} or {molecule} or {group} :l
groupID1, groupID2, ... = list of N group IDs :pre
zero or more keyword/value pairs may be appended :l
-keyword = {langevin} or {temp} or {iso} or {aniso} or {x} or {y} or {z} or {couple} or {tparam} or {pchain} or {dilate} or {force} or {torque} or {infile} :l
+keyword = {langevin} or {reinit} or {temp} or {iso} or {aniso} or {x} or {y} or {z} or {couple} or {tparam} or {pchain} or {dilate} or {force} or {torque} or {infile} :l
{langevin} values = Tstart Tstop Tperiod seed
Tstart,Tstop = desired temperature at start/stop of run (temperature units)
Tdamp = temperature damping parameter (time units)
seed = random number seed to use for white noise (positive integer)
+ {reinit} = {yes} or {no}
{temp} values = Tstart Tstop Tdamp
Tstart,Tstop = desired temperature at start/stop of run (temperature units)
Tdamp = temperature damping parameter (time units)
@@ -68,10 +69,10 @@ keyword = {langevin} or {temp} or {iso} or {aniso} or {x} or {y} or {z} or {coup
[Examples:]
-fix 1 clump rigid single
+fix 1 clump rigid single reinit yes
fix 1 clump rigid/small molecule
fix 1 clump rigid single force 1 off off on langevin 1.0 1.0 1.0 428984
-fix 1 polychains rigid/nvt molecule temp 1.0 1.0 5.0
+fix 1 polychains rigid/nvt molecule temp 1.0 1.0 5.0 reinit no
fix 1 polychains rigid molecule force 1*5 off off off force 6*10 off off on
fix 1 polychains rigid/small molecule langevin 1.0 1.0 1.0 428984
fix 2 fluid rigid group 3 clump1 clump2 clump3 torque * off off off
@@ -87,7 +88,12 @@ means that each timestep the total force and torque on each rigid body
is computed as the sum of the forces and torques on its constituent
particles. The coordinates, velocities, and orientations of the atoms
in each body are then updated so that the body moves and rotates as a
-single entity.
+single entity. This is implemented by creating internal data structures
+for each rigid body and performing time integration on these data
+structures. Positions, velocities, and orientations of the constituent
+particles are regenerated from the rigid body data structures in every
+time step. This restricts which operations and fixes can be applied to
+rigid bodies. See below for a detailed discussion.
Examples of large rigid bodies are a colloidal particle, or portions
of a biomolecule such as a protein.
@@ -148,8 +154,9 @@ differences may accumulate to produce divergent trajectories.
NOTE: You should not update the atoms in rigid bodies via other
time-integration fixes (e.g. "fix nve"_fix_nve.html, "fix
-nvt"_fix_nh.html, "fix npt"_fix_nh.html), or you will be integrating
-their motion more than once each timestep. When performing a hybrid
+nvt"_fix_nh.html, "fix npt"_fix_nh.html, "fix move"_fix_move.html),
+or you will have conflicting updates to positions and velocities
+resulting in unphysical behavior in most cases. When performing a hybrid
simulation with some atoms in rigid bodies, and some not, a separate
time integration fix like "fix nve"_fix_nve.html or "fix
nvt"_fix_nh.html should be used for the non-rigid particles.
@@ -165,23 +172,29 @@ setting the force on them to 0.0 (via the "fix
setforce"_fix_setforce.html command), and integrating them as usual
(e.g. via the "fix nve"_fix_nve.html command).
-NOTE: The aggregate properties of each rigid body are calculated one
-time at the start of the first simulation run after these fixes are
-specified. The properties include the position and velocity of the
-center-of-mass of the body, its moments of inertia, and its angular
-momentum. This is done using the properties of the constituent atoms
-of the body at that point in time (or see the {infile} keyword
-option). Thereafter, changing properties of individual atoms in the
-body will have no effect on a rigid body's dynamics, unless they
-affect the "pair_style"_pair_style.html interactions that individual
-particles are part of. For example, you might think you could
-displace the atoms in a body or add a large velocity to each atom in a
-body to make it move in a desired direction before a 2nd run is
+IMPORTANT NOTE: The aggregate properties of each rigid body are
+calculated at the start of a simulation run and are maintained in
+internal data structures. The properties include the position and
+velocity of the center-of-mass of the body, its moments of inertia, and
+its angular momentum. This is done using the properties of the
+constituent atoms of the body at that point in time (or see the {infile}
+keyword option). Thereafter, changing these properties of individual
+atoms in the body will have no effect on a rigid body's dynamics, unless
+they effect any computation of per-atom forces or torques. If the
+keyword {reinit} is set to {yes} (the default), the rigid body data
+structures will be recreated at the beginning of each {run} command;
+if the keyword {reinit} is set to {no}, the rigid body data structures
+will be built only at the very first {run} command and maintained for
+as long as the rigid fix is defined. For example, you might think you
+could displace the atoms in a body or add a large velocity to each atom
+in a body to make it move in a desired direction before a 2nd run is
performed, using the "set"_set.html or
"displace_atoms"_displace_atoms.html or "velocity"_velocity.html
-command. But these commands will not affect the internal attributes
-of the body, and the position and velocity of individual atoms in the
-body will be reset when time integration starts.
+commands. But these commands will not affect the internal attributes
+of the body unless {reinit} is set to {yes}. With {reinit} set to {no}
+(or using the {infile} option, which implies {reinit} {no}) the position
+and velocity of individual atoms in the body will be reset when time
+integration starts again.
:line
@@ -401,6 +414,14 @@ couple none :pre
The keyword/value option pairs are used in the following ways.
+The {reinit} keyword determines, whether the rigid body properties
+are reinitialized between run commands. With the option {yes} (the
+default) this is done, with the option {no} this is not done. Turning
+off the reinitialization can be helpful to protect rigid bodies against
+unphysical manipulations between runs or when properties cannot be
+easily recomputed (e.g. when read from a file). When using the {infile}
+keyword, the {reinit} option is automatically set to {no}.
+
The {langevin} and {temp} and {tparam} keywords perform thermostatting
of the rigid bodies, altering both their translational and rotational
degrees of freedom. What is meant by "temperature" of a collection of
@@ -778,7 +799,7 @@ exclude, "fix shake"_fix_shake.html
The option defaults are force * on on on and torque * on on on,
meaning all rigid bodies are acted on by center-of-mass force and
-torque. Also Tchain = Pchain = 10, Titer = 1, Torder = 3.
+torque. Also Tchain = Pchain = 10, Titer = 1, Torder = 3, reinit = yes.
:line
diff --git a/src/RIGID/fix_rigid.cpp b/src/RIGID/fix_rigid.cpp
index e9c06dee5c42b27806fb9133bfc38abf5bcc6fca..83022cbcba9de4e98965b44786c9d20c48f47896 100644
--- a/src/RIGID/fix_rigid.cpp
+++ b/src/RIGID/fix_rigid.cpp
@@ -267,6 +267,8 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
int seed;
langflag = 0;
+ reinitflag = 1;
+
tstat_flag = 0;
pstat_flag = 0;
allremap = 1;
@@ -501,6 +503,14 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
infile = new char[n];
strcpy(infile,arg[iarg+1]);
restart_file = 1;
+ reinitflag = 0;
+ iarg += 2;
+
+ } else if (strcmp(arg[iarg],"reinit") == 0) {
+ if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid/small 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 error->all(FLERR,"Illegal fix rigid command");
@@ -679,15 +689,15 @@ void FixRigid::init()
if (strstr(update->integrate_style,"respa"))
step_respa = ((Respa *) update->integrate)->step;
- // setup rigid bodies, using current atom info
- // 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 rigid bodies, using current atom info. if reinitflag is not set,
+ // do the initialization only 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 (!setupflag) {
+ if (reinitflag || !setupflag) {
setup_bodies_static();
if (!infile) setup_bodies_dynamic();
setupflag = 1;
diff --git a/src/RIGID/fix_rigid.h b/src/RIGID/fix_rigid.h
index a6d1f65e1c089aa9c76125bc46679259105c64c1..12439d42cfbc452022595dd9f6dfd27285096ea0 100644
--- a/src/RIGID/fix_rigid.h
+++ b/src/RIGID/fix_rigid.h
@@ -104,6 +104,7 @@ class FixRigid : public Fix {
int extended; // 1 if any particles have extended attributes
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
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_small.cpp b/src/RIGID/fix_rigid_small.cpp
index 2d8e736a1eef1ebf811c555259695ace712b8038..f40b52e1626c937ea1079f3f552fe95b58636570 100644
--- a/src/RIGID/fix_rigid_small.cpp
+++ b/src/RIGID/fix_rigid_small.cpp
@@ -138,6 +138,7 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
langflag = 0;
infile = NULL;
onemols = NULL;
+ reinitflag = 1;
tstat_flag = 0;
pstat_flag = 0;
@@ -173,6 +174,7 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
error->all(FLERR,"Fix rigid/small langevin period must be > 0.0");
if (seed <= 0) error->all(FLERR,"Illegal fix rigid/small command");
iarg += 5;
+
} else if (strcmp(arg[iarg],"infile") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid/small command");
delete [] infile;
@@ -180,7 +182,16 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
infile = new char[n];
strcpy(infile,arg[iarg+1]);
restart_file = 1;
+ reinitflag = 0;
+ iarg += 2;
+
+ } else if (strcmp(arg[iarg],"reinit") == 0) {
+ if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid/small 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/small command");
iarg += 2;
+
} else if (strcmp(arg[iarg],"mol") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid/small command");
int imol = atom->find_molecule(arg[iarg+1]);
@@ -520,14 +531,15 @@ void FixRigidSmall::init()
}
/* ----------------------------------------------------------------------
- setup static/dynamic properties of rigid bodies, using current atom info
- 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
+ setup static/dynamic properties of rigid bodies, using current atom info.
+ if reinitflag is not set, do the initialization only once, b/c properties
+ may not be re-computable especially if overlapping particles or bodies
+ are 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
setup_bodies_static() invokes pre_neighbor itself
@@ -535,9 +547,13 @@ void FixRigidSmall::init()
void FixRigidSmall::setup_pre_neighbor()
{
- if (!setupflag) setup_bodies_static();
+ if (reinitflag || !setupflag)
+ setup_bodies_static();
else pre_neighbor();
- if (!setupflag && !infile) setup_bodies_dynamic();
+
+ if ((reinitflag || !setupflag) && !infile)
+ setup_bodies_dynamic();
+
setupflag = 1;
}
diff --git a/src/RIGID/fix_rigid_small.h b/src/RIGID/fix_rigid_small.h
index 9c89bab885cf2365a24351b56e3eafc41510d853..72e32719c0071d77f023304a1d89b26c08bd5773 100644
--- a/src/RIGID/fix_rigid_small.h
+++ b/src/RIGID/fix_rigid_small.h
@@ -131,6 +131,7 @@ class FixRigidSmall : public Fix {
int extended; // 1 if any particles have extended attributes
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 POINT,SPHERE,ELLIPSOID,LINE,TRIANGLE,DIPOLE; // bitmasks for eflags
int OMEGA,ANGMOM,TORQUE;