Merge pull request #517 from akohlmey/select-rigid-reinit-option

Add `reinit` keyword to rigid body fixes
This commit is contained in:
sjplimp
2017-06-15 08:46:29 -06:00
committed by GitHub
5 changed files with 88 additions and 39 deletions

View File

@ -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

View File

@ -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;

View File

@ -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

View File

@ -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;
}

View File

@ -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;