lammps/doc/fix_rigid.html

403 lines
19 KiB
HTML

<HTML>
<CENTER><A HREF = "http://lammps.sandia.gov">LAMMPS WWW Site</A> - <A HREF = "Manual.html">LAMMPS Documentation</A> - <A HREF = "Section_commands.html#comm">LAMMPS Commands</A>
</CENTER>
<HR>
<H3>fix rigid command
</H3>
<H3>fix rigid/nve command
</H3>
<H3>fix rigid/nvt command
</H3>
<P><B>Syntax:</B>
</P>
<PRE>fix ID group-ID style bodystyle args keyword values ...
</PRE>
<UL><LI>ID, group-ID are documented in <A HREF = "fix.html">fix</A> command
<LI>style = <I>rigid</I> or <I>rigid/nve</I> or <I>rigid/nvt</I>
<LI>bodystyle = <I>single</I> or <I>molecule</I> or <I>group</I>
<PRE> <I>single</I> args = none
<I>molecule</I> args = none
<I>group</I> args = N groupID1 groupID2 ...
N = # of groups
groupID1, groupID2, ... = list of N group IDs
</PRE>
<LI>zero or more keyword/value pairs may be appended
<LI>keyword = <I>langevin</I> or <I>temp</I> or <I>tparam</I> or <I>force</I> or <I>torque</I>
<PRE> <I>langevin</I> 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)
<I>temp</I> values = Tstart Tstop Tdamp
Tstart,Tstop = desired temperature at start/stop of run (temperature units)
Tdamp = temperature damping parameter (time units)
<I>tparam</I> values = Tchain Titer Torder
Tchain = length of Nose/Hoover thermostat chain
Titer = number of thermostat iterations performed
Torder = 3 or 5 = Yoshida-Suzuki integration parameters
<I>force</I> values = M xflag yflag zflag
M = which rigid body from 1-Nbody (see asterisk form below)
xflag,yflag,zflag = off/on if component of center-of-mass force is active
<I>torque</I> values = M xflag yflag zflag
M = which rigid body from 1-Nbody (see asterisk form below)
xflag,yflag,zflag = off/on if component of center-of-mass torque is active
</PRE>
</UL>
<P><B>Examples:</B>
</P>
<PRE>fix 1 clump rigid single
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 molecule force 1*5 off off off force 6*10 off off on
fix 2 fluid rigid group 3 clump1 clump2 clump3 torque * off off off
</PRE>
<P><B>Description:</B>
</P>
<P>Treat one or more sets of atoms as independent rigid bodies. This
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 and the coordinates, velocities, and orientations of the
atoms in each body are updated so that the body moves and rotates as a
single entity.
</P>
<P>Examples of large rigid bodies are a large colloidal particle, or
portions of a large biomolecule such as a protein.
</P>
<P>Example of small rigid bodies are patchy nanoparticles, such as those
modeled in <A HREF = "#Zhang">this paper</A> by Sharon Glotzer's group, clumps of
granular particles, lipid molecules consiting of one or more point
dipoles connected to other spheroids or ellipsoids, and coarse-grain
models of nano or colloidal particles consisting of a small number of
constituent particles. Note that the <A HREF = "fix_shake.html">fix shake</A>
command can also be used to rigidify small molecules of 2, 3, or 4
atoms, e.g. water molecules. That fix treats the constituent atoms as
point masses.
</P>
<P>These fixes also update the positions and velocities of the atoms in
each rigid body via time integration. The <I>rigid</I> and <I>rigid/nve</I>
styles do this via constant NVE integration. The only difference is
that the <I>rigid</I> style uses an integration technique based on
Richardson iterations. The <I>rigid/nve</I> style uses the methods
described in the paper by <A HREF = "#Miller">Miller</A>, which are thought to
provide better energy conservation than an iterative approach.
</P>
<P>The <I>rigid/nvt</I> style performs constant NVT integration using a
Nose/Hoover thermostat with chains as described originally in
<A HREF = "#Hoover">(Hoover)</A> and <A HREF = "#Martyna">(Martyna)</A>, which thermostats both
the translational and rotational degrees of freedom of the rigid
bodies. The rigid-body algorithm used by <I>rigid/nvt</I> is described in
the paper by <A HREF = "#Kamberaj">Kamberaj</A>.
</P>
<P>IMPORTANT NOTE: You should not update the atoms in rigid bodies via
other time-integration fixes (e.g. nve, nvt, npt), or you will be
integrating their motion more than once each timestep.
</P>
<P>IMPORTANT NOTE: These fixes are overkill if you simply want to hold a
collection of atoms stationary or have them move with a constant
velocity. A simpler way to hold atoms stationary is to not include
those atoms in your time integration fix. E.g. use "fix 1 mobile nve"
instead of "fix 1 all nve", where "mobile" is the group of atoms that
you want to move. You can move atoms with a constant velocity by
assigning them an initial velocity (via the <A HREF = "velocity.html">velocity</A>
command), setting the force on them to 0.0 (via the <A HREF = "fix_setforce.html">fix
setforce</A> command), and integrating them as usual
(e.g. via the <A HREF = "fix_nve.html">fix nve</A> command).
</P>
<HR>
<P>The constituent particles within a rigid body can be point particles
(the default in LAMMPS) or finite-size particles, such as spheres and
ellipsoids. See the <A HREF = "atom_style.html">atom_style sphere and ellipsoid</A>
commands for more details on these kinds of particles. Finite-size
particles contribute differently to the moment of inertia of a rigid
body than do point particles. Finite-size particles can also
experience torque (e.g. due to <A HREF = "pair_gran.html">frictional granular
interactions</A>) and have an orientation. These
contributions are accounted for by these fixes.
</P>
<P>Forces between particles within a body do not contribute to the
external force or torque on the body. Thus for computational
efficiency, you may wish to turn off pairwise and bond interactions
between particles within each rigid body. The <A HREF = "neigh_modify.html">neigh_modify
exclude</A> and <A HREF = "delete_bonds.html">delete_bonds</A>
commands are used to do this. For finite-size particles this also
means the particles can be highly overlapped when creating the rigid
body.
</P>
<HR>
<P>Each body must have two or more atoms. An atom can belong to at most
one rigid body. Which atoms are in which bodies can be defined via
several options.
</P>
<P>For bodystyle <I>single</I> the entire fix group of atoms is treated as one
rigid body.
</P>
<P>For bodystyle <I>molecule</I>, each set of atoms in the fix group with a
different molecule ID is treated as a rigid body.
</P>
<P>For bodystyle <I>group</I>, each of the listed groups is treated as a
separate rigid body. Only atoms that are also in the fix group are
included in each rigid body.
</P>
<P>IMPORTANT NOTE: To compute the initial center-of-mass position and
other properties of each rigid body, the image flags for each atom in
the body are used to "unwrap" the atom coordinates. Thus you must
insure that these image flags are consistent so that the unwrapping
creates a valid rigid body (one where the atoms are close together),
particularly if the atoms in a single rigid body straddle a periodic
boundary. This means the input data file or restart file must define
the image flags for each atom consistently or that you have used the
<A HREF = "set.html">set</A> command to specify them correctly. If a dimension is
non-periodic then the image flag of each atom must be 0 in that
dimension, else an error is generated.
</P>
<P>By default, each rigid body is acted on by other atoms which induce an
external force and torque on its center of mass, causing it to
translate and rotate. Components of the external center-of-mass force
and torque can be turned off by the <I>force</I> and <I>torque</I> keywords.
This may be useful if you wish a body to rotate but not translate, or
vice versa, or if you wish it to rotate or translate continuously
unaffected by interactions with other particles. Note that if you
expect a rigid body not to move or rotate by using these keywords, you
must insure its initial center-of-mass translational or angular
velocity is 0.0. Otherwise the initial translational or angular
momentum the body has will persist.
</P>
<P>An xflag, yflag, or zflag set to <I>off</I> means turn off the component of
force of torque in that dimension. A setting of <I>on</I> means turn on
the component, which is the default. Which rigid body(s) the settings
apply to is determined by the first argument of the <I>force</I> and
<I>torque</I> keywords. It can be an integer M from 1 to Nbody, where
Nbody is the number of rigid bodies defined. A wild-card asterisk can
be used in place of, or in conjunction with, the M argument to set the
flags for multiple rigid bodies. This takes the form "*" or "*n" or
"n*" or "m*n". If N = the number of rigid bodies, then an asterisk
with no numeric values means all bodies from 1 to N. A leading
asterisk means all bodies from 1 to n (inclusive). A trailing
asterisk means all bodies from n to N (inclusive). A middle asterisk
means all types from m to n (inclusive). Note that you can use the
<I>force</I> or <I>torque</I> keywords as many times as you like. If a
particular rigid body has its component flags set multiple times, the
settings from the final keyword are used.
</P>
<P>For computational efficiency, you may wish to turn off pairwise and
bond interactions within each rigid body, as they no longer contribute
to the motion. The <A HREF = "neigh_modify.html">neigh_modify exclude</A> and
<A HREF = "delete_bonds.html">delete_bonds</A> commands are used to do this.
</P>
<P>For computational efficiency, you should typically define one fix
rigid which includes all the desired rigid bodies. LAMMPS will allow
multiple rigid fixes to be defined, but it is more expensive.
</P>
<HR>
<P>The keyword/value option pairs are used in the following ways.
</P>
<P>The <I>langevin</I> and <I>temp</I> and <I>tparam</I> 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
rigid bodies and how it can be monitored via the fix output is
discussed below.
</P>
<P>The <I>langevin</I> keyword applies a Langevin thermostat to the constant
NVE time integration performed by either the <I>rigid</I> or <I>rigid/nve</I>
styles. It cannot be used with the <I>rigid/nvt</I> style. The desired
temperature at each timestep is a ramped value during the run from
<I>Tstart</I> to <I>Tstop</I>. The <I>Tdamp</I> parameter is specified in time units
and determines how rapidly the temperature is relaxed. For example, a
value of 100.0 means to relax the temperature in a timespan of
(roughly) 100 time units (tau or fmsec or psec - see the
<A HREF = "units.html">units</A> command). The random # <I>seed</I> must be a positive
integer. The way the Langevin thermostatting operates is explained on
the <A HREF = "fix_langevin.html">fix langevin</A> doc page.
</P>
<P>The <I>temp</I> and <I>tparam</I> keywords apply a Nose/Hoover thermostat to the
NVT time integration performed by the <I>rigid/nvt</I> style. They cannot
be used with the <I>rigid</I> or <I>rigid/nve</I> styles. The desired
temperature at each timestep is a ramped value during the run from
<I>Tstart</I> to <I>Tstop</I>. The <I>Tdamp</I> parameter is specified in time units
and determines how rapidly the temperature is relaxed. For example, a
value of 100.0 means to relax the temperature in a timespan of
(roughly) 100 time units (tau or fmsec or psec - see the
<A HREF = "units.html">units</A> command).
</P>
<P>Nose/Hoover chains are used in conjunction with this thermostat. The
<I>tparam</I> keyword can optionally be used to change the chain settings
used. <I>Tchain</I> is the number of thermostats in the Nose Hoover chain.
This value, along with <I>Tdamp</I> can be varied to dampen undesirable
oscillations in temperature that can occur in a simulation. As a rule
of thumb, increasing the chain length should lead to smaller
oscillations.
</P>
<P>IMPORTANT NOTE: There are alternate ways to thermostat a system of
rigid bodies. You can use <A HREF = "fix_langevin.html">fix langevin</A> to treat
the individual particles in the rigid bodies as effectively immersed
in an implicit solvent, e.g. a Brownian dynamics model. For hybrid
systems with both rigid bodies and solvent particles, you can
thermostat only the solvent particles that surround one or more rigid
bodies by appropriate choice of groups in the compute and fix commands
for temperature and thermostatting. The solvent interactions with the
rigid bodies should then effectively thermostat the rigid body
temperature as well without use of the Langevin or Nose/Hoover options
associated with the fix rigid commands.
</P>
<HR>
<P>The keyword/value option pairs are used in the following ways.
</P>
<P>If you use a <A HREF = "compute.html">temperature compute</A> with a group that
includes particles in rigid bodies, the degrees-of-freedom removed by
each rigid body are accounted for in the temperature (and pressure)
computation, but only if the temperature group includes all the
particles in a particular rigid body.
</P>
<P>A 3d rigid body has 6 degrees of freedom (3 translational, 3
rotational), except for a collection of point particles lying on a
straight line, which has only 5, e.g a dimer. A 2d rigid body has 3
degrees of freedom (2 translational, 1 rotational).
</P>
<P>IMPORTANT NOTE: You may wish to explicitly subtract additional
degrees-of-freedom if you use the <I>force</I> and <I>torque</I> keywords to
eliminate certain motions of one or more rigid bodies. LAMMPS does
not do this automatically.
</P>
<P>The rigid body contribution to the pressure of the system (virial) is
also accounted for by this fix.
</P>
<P>IMPORTANT NOTE: The periodic image flags of atoms in rigid bodies are
altered so that the rigid body can be reconstructed correctly when it
straddles periodic boundaries. The atom image flags are not
incremented/decremented as they would be for non-rigid atoms as the
rigid body crosses periodic boundaries. This means you cannot
interpret them as you normally would. For example, the image flag
values written to a <A HREF = "dump.html">dump file</A> will be different than they
would be if the atoms were not in a rigid body. Likewise the <A HREF = "compute_msd.html">compute
msd</A> will not compute the expected mean-squared
displacement for such atoms if the body moves across periodic
boundaries. It also means that if you have bonds between a pair of
rigid bodies and the bond straddles a periodic boundary, you cannot
use the <A HREF = "replicate.html">replicate</A> command to increase the system
size. Note that this fix does define image flags for each rigid body,
which are incremented when the rigid body crosses a periodic boundary
in the usual way. These image flags have the same meaning as atom
images (see the "dump" command) and can be accessed and output as
described below.
</P>
<HR>
<P><B>Restart, fix_modify, output, run start/stop, minimize info:</B>
</P>
<P>No information about the <I>rigid</I> and <I>rigid/nve</I> fixes are written to
<A HREF = "restart.html">binary restart files</A>. For style <I>rigid/nvt</I> the state
of the Nose/Hoover thermostat is written to <A HREF = "restart.html">binary restart
files</A>. See the <A HREF = "read_restart.html">read_restart</A> command
for info on how to re-specify a fix in an input script that reads a
restart file, so that the operation of the fix continues in an
uninterrupted fashion.
</P>
<P>The <A HREF = "fix_modify.html">fix_modify</A> <I>energy</I> option is supported by the
rigid/nvt fix to add the energy change induced by the thermostatting
to the system's potential energy as part of <A HREF = "thermo_style.html">thermodynamic
output</A>.
</P>
<P>The rigid and rigid/nve fixes computes a global scalar which can be
accessed by various <A HREF = "Section_howto.html#howto_15">output commands</A>.
The scalar value calculated by these fixes is "intensive". The scalar
is the current temperature of the collection of rigid bodies. This is
averaged over all rigid bodies and their translational and rotational
degrees of freedom. The translational energy of a rigid body is 1/2 m
v^2, where m = total mass of the body and v = the velocity of its
center of mass. The rotational energy of a rigid body is 1/2 I w^2,
where I = the moment of inertia tensor of the body and w = its angular
velocity. Degrees of freedom constrained by the <I>force</I> and <I>torque</I>
keywords are removed from this calculation.
</P>
<P>The rigid/nvt fix computes a global scalar which can be accessed by
various <A HREF = "Section_howto.html#howto_15">output commands</A>. The scalar
value calculated by the rigid/nvt fix is "extensive". The scalar is
the cumulative energy change due to the thermostatting the fix
performs.
</P>
<P>All of these fixes compute a global array of values which can be
accessed by various <A HREF = "Section_howto.html#howto_15">output commands</A>.
The number of rows in the array is equal to the number of rigid
bodies. The number of columns is 15. Thus for each rigid body, 15
values are stored: the xyz coords of the center of mass (COM), the xyz
components of the COM velocity, the xyz components of the force acting
on the COM, the xyz components of the torque acting on the COM, and
the xyz image flags of the COM, which have the same meaning as image
flags for atom positions (see the "dump" command). The force and
torque values in the array are not affected by the <I>force</I> and
<I>torque</I> keywords in the fix rigid command; they reflect values before
any changes are made by those keywords.
</P>
<P>The ordering of the rigid bodies (by row in the array) is as follows.
For the <I>single</I> keyword there is just one rigid body. For the
<I>molecule</I> keyword, the bodies are ordered by ascending molecule ID.
For the <I>group</I> keyword, the list of group IDs determines the ordering
of bodies.
</P>
<P>The array values calculated by these fixes are "intensive", meaning
they are independent of the number of atoms in the simulation.
</P>
<P>No parameter of these fixes can be used with the <I>start/stop</I> keywords
of the <A HREF = "run.html">run</A> command. These fixes are not invoked during
<A HREF = "minimize.html">energy minimization</A>.
</P>
<P><B>Restrictions:</B>
</P>
<P>These fixes performs an MPI_Allreduce each timestep that is
proportional in length to the number of rigid bodies. Hence they will
not scale well in parallel if large numbers of rigid bodies are
simulated.
</P>
<P><B>Related commands:</B>
</P>
<P><A HREF = "delete_bonds.html">delete_bonds</A>, <A HREF = "neigh_modify.html">neigh_modify</A>
exclude
</P>
<P><B>Default:</B>
</P>
<P>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 = 10, Titer = 1, Torder = 3.
</P>
<HR>
<A NAME = "Hoover"></A>
<P><B>(Hoover)</B> Hoover, Phys Rev A, 31, 1695 (1985).
</P>
<A NAME = "Kamberaj"></A>
<P><B>(Kamberaj)</B> Kamberaj, Low, Neal, J Chem Phys, 122, 224114 (2005).
</P>
<A NAME = "Martyna"></A>
<P><B>(Martyna)</B> Martyna, Klein, Tuckerman, J Chem Phys, 97, 2635 (1992);
Martyna, Tuckerman, Tobias, Klein, Mol Phys, 87, 1117.
</P>
<A NAME = "Miller"></A>
<P><B>(Miller)</B> Miller, Eleftheriou, Pattnaik, Ndirango, and Newns,
J Chem Phys, 116, 8649 (2002).
</P>
<A NAME = "Zhang"></A>
<P><B>(Zhang)</B> Zhang, Glotzer, Nanoletters, 4, 1407-1413 (2004).
</P>
</HTML>