jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
lammps/doc/fix_rigid.html

747 lines
38 KiB
HTML
Raw Blame History

<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>
<H3>fix rigid/npt command
</H3>
<H3>fix rigid/nph command
</H3>
<H3>fix rigid/small 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> or <I>rigid/npt</I> or <I>rigid/nph</I> or <I>rigid/small</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>iso</I> or <I>aniso</I> or <I>x</I> or <I>y</I> or <I>z</I> or <I>couple</I> or <I>tparam</I> or <I>pchain</I> or <I>dilate</I> or <I>force</I> or <I>torque</I> or <I>infile</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>iso</I> or <I>aniso</I> values = Pstart Pstop Pdamp
Pstart,Pstop = scalar external pressure at start/end of run (pressure units)
Pdamp = pressure damping parameter (time units)
<I>x</I> or <I>y</I> or <I>z</I> values = Pstart Pstop Pdamp
Pstart,Pstop = external stress tensor component at start/end of run (pressure units)
Pdamp = stress damping parameter (time units)
<I>couple</I> = <I>none</I> or <I>xyz</I> or <I>xy</I> or <I>yz</I> or <I>xz</I>
<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>pchain</I> values = Pchain
Pchain = length of the Nose/Hoover thermostat chain coupled with the barostat
<I>dilate</I> value = dilate-group-ID
dilate-group-ID = only dilate atoms in this group due to barostat volume changes
<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
<I>infile</I> filename
filename = file with per-body values of mass, center-of-mass, moments of inertia
</PRE>
</UL>
<P><B>Examples:</B>
</P>
<PRE>fix 1 clump rigid single
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 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
fix 1 rods rigid/npt molecule temp 300.0 300.0 100.0 iso 0.5 0.5 10.0
fix 1 particles rigid/npt molecule temp 1.0 1.0 5.0 x 0.5 0.5 1.0 z 0.5 0.5 1.0 couple xz
fix 1 water rigid/nph molecule iso 0.5 0.5 1.0
</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, irregular
particles built from line segments (2d) or triangles (3d), 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, in the NVE, NVT, NPT, or NPH
ensemble, as described below.
</P>
<P>There are two main variants of this fix, fix rigid and fix
rigid/small. The NVE/NVT/NPT/NHT versions belong to one of the two
variants, as their style names indicate.
</P>
<P>IMPORTANT NOTE: Not all of the bodystyle options and keyword/value
options are available for both the <I>rigid</I> and <I>rigid/small</I> variants.
See details below.
</P>
<P>The <I>rigid</I> variant is typically the best choice for a system with a
small number of large rigid bodies, each of which can extend across
the domain of many processors. It operates by creating a single
global list of rigid bodies, which all processors contribute to.
MPI_Allreduce operations are performed each timestep to sum the
contributions from each processor to the force and torque on all the
bodies. This operation will not scale well in parallel if large
numbers of rigid bodies are simulated.
</P>
<P>The <I>rigid/small</I> variant is typically best for a system with a large
number of small rigid bodies. Each body is assigned to the atom
closest to the geometrical center of the body. The fix operates using
local lists of rigid bodies owned by each processor and information is
exchanged and summed via local communication between neighboring
processors when ghost atom info is accumlated.
</P>
<P>IMPORTANT NOTE: To use <I>rigid/small</I> the ghost atom cutoff must be
large enough to span the distance between the atom that owns the body
and every other atom in the body. This distance value is printed out
when the rigid bodies are defined. If the
<A HREF = "pair_style.html">pair_style</A> cutoff plus neighbor skin does not span
this distance, then you should use the <A HREF = "communicate.html">communicate
cutoff</A> command with a setting epsilon larger than
the distance.
</P>
<P>Which of the two variants is faster for a particular problem is hard
to predict. The best way to decide is to perform a short test run.
Both variants should give identical numerical answers for short runs.
Long runs should give statistically similar results, but round-off
differences will accumulate to produce divergent trajectories.
</P>
<P>IMPORTANT NOTE: You should not update the atoms in rigid bodies via
other time-integration fixes (e.g. <A HREF = "fix_nve.html">fix nve</A>, <A HREF = "fix_nvt.html">fix
nvt</A>, <A HREF = "fix_npt.html">fix npt</A>), or you will be integrating
their motion more than once each timestep. When performing a hybrid
simulation with some atoms in rigid bodies, and some not, a separate
time integration fix like <A HREF = "fix_nve.html">fix nve</A> or <A HREF = "fix_nh.html">fix
nvt</A> should be used for the non-rigid particles.
</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>Each rigid 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. This option is only allowed for fix rigid and its
sub-styles.
</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. This option is
allowed for fix rigid and fix rigid/small, and their sub-styles. Note
that atoms with a molecule ID = 0 will be treated as a single rigid
body. For a system with atomic solvent (typically this is atoms with
molecule ID = 0) surrounding rigid bodies, this may not be what you
want. Thus you should be careful to use a fix group that only
includes atoms you want to be part of rigid bodies.
</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. This option is only allowed for fix
rigid and its sub-styles.
</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>The <I>force</I> and <I>torque</I> keywords discussed next are only allowed for
fix rigid and its sub-styles.
</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 or fix rigid/small command 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 constituent particles within a rigid body can be point particles
(the default in LAMMPS) or finite-size particles, such as spheres or
ellipsoids or line segments or triangles. See the <A HREF = "atom_style.html">atom_style sphere
and ellipsoid and line and tri</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>The <I>rigid</I> and <I>rigid/small</I> and <I>rigid/nve</I> styles perform constant
NVE time integration. The only difference is that the <I>rigid</I> and
<I>rigid/small</I> styles use 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>The <I>rigid/npt</I> and <I>rigid/nph</I> styles perform constant NPT or NPH
integration using a Nose/Hoover barostat with chains. For the NPT
case, the same Nose/Hoover thermostat is also used as with
<I>rigid/nvt</I>.
</P>
<P>The barostat parameters are specified using one or more of the <I>iso</I>,
<I>aniso</I>, <I>x</I>, <I>y</I>, <I>z</I> and <I>couple</I> keywords. These keywords give you
the ability to specify 3 diagonal components of the external stress
tensor, and to couple these components together so that the dimensions
they represent are varied together during a constant-pressure
simulation. The effects of these keywords are similar to those
defined in <A HREF = "fix_nh.html">fix npt/nph</A>
</P>
<P>NOTE: Currently the <I>rigid/npt</I> and <I>rigid/nph</I> styles do not support
triclinic (non-orthongonal) boxes.
</P>
<P>The target pressures for each of the 6 components of the stress tensor
can be specified independently via the <I>x</I>, <I>y</I>, <I>z</I> keywords, which
correspond to the 3 simulation box dimensions. For each component,
the external pressure or tensor component at each timestep is a ramped
value during the run from <I>Pstart</I> to <I>Pstop</I>. If a target pressure is
specified for a component, then the corresponding box dimension will
change during a simulation. For example, if the <I>y</I> keyword is used,
the y-box length will change. A box dimension will not change if that
component is not specified, although you have the option to change
that dimension via the <A HREF = "fix_deform.html">fix deform</A> command.
</P>
<P>For all barostat keywords, the <I>Pdamp</I> parameter operates like the
<I>Tdamp</I> parameter, determining the time scale on which pressure is
relaxed. For example, a value of 10.0 means to relax the pressure in
a timespan of (roughly) 10 time units (e.g. tau or fmsec or psec - see
the <A HREF = "units.html">units</A> command).
</P>
<P>Regardless of what atoms are in the fix group (the only atoms which
are time integrated), a global pressure or stress tensor is computed
for all atoms. Similarly, when the size of the simulation box is
changed, all atoms are re-scaled to new positions, unless the keyword
<I>dilate</I> is specified with a <I>dilate-group-ID</I> for a group that
represents a subset of the atoms. This can be useful, for example, to
leave the coordinates of atoms in a solid substrate unchanged and
controlling the pressure of a surrounding fluid. Another example is a
system consisting of rigid bodies and point particles where the
barostat is only coupled with the rigid bodies. This option should be
used with care, since it can be unphysical to dilate some atoms and
not others, because it can introduce large, instantaneous
displacements between a pair of atoms (one dilated, one not) that are
far from the dilation origin.
</P>
<P>The <I>couple</I> keyword allows two or three of the diagonal components of
the pressure tensor to be "coupled" together. The value specified
with the keyword determines which are coupled. For example, <I>xz</I>
means the <I>Pxx</I> and <I>Pzz</I> components of the stress tensor are coupled.
<I>Xyz</I> means all 3 diagonal components are coupled. Coupling means two
things: the instantaneous stress will be computed as an average of the
corresponding diagonal components, and the coupled box dimensions will
be changed together in lockstep, meaning coupled dimensions will be
dilated or contracted by the same percentage every timestep. The
<I>Pstart</I>, <I>Pstop</I>, <I>Pdamp</I> parameters for any coupled dimensions must
be identical. <I>Couple xyz</I> can be used for a 2d simulation; the <I>z</I>
dimension is simply ignored.
</P>
<P>The <I>iso</I> and <I>aniso</I> keywords are simply shortcuts that are
equivalent to specifying several other keywords together.
</P>
<P>The keyword <I>iso</I> means couple all 3 diagonal components together when
pressure is computed (hydrostatic pressure), and dilate/contract the
dimensions together. Using "iso Pstart Pstop Pdamp" is the same as
specifying these 4 keywords:
</P>
<PRE>x Pstart Pstop Pdamp
y Pstart Pstop Pdamp
z Pstart Pstop Pdamp
couple xyz
</PRE>
<P>The keyword <I>aniso</I> means <I>x</I>, <I>y</I>, and <I>z</I> dimensions are controlled
independently using the <I>Pxx</I>, <I>Pyy</I>, and <I>Pzz</I> components of the
stress tensor as the driving forces, and the specified scalar external
pressure. Using "aniso Pstart Pstop Pdamp" is the same as specifying
these 4 keywords:
</P>
<PRE>x Pstart Pstop Pdamp
y Pstart Pstop Pdamp
z Pstart Pstop Pdamp
couple none
</PRE>
<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/small</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.
</P>
<P>The way that Langevin thermostatting operates is explained on the <A HREF = "fix_langevin.html">fix
langevin</A> doc page. If you wish to simply viscously
damp the rotational motion without thermostatting, you can set
<I>Tstart</I> and <I>Tstop</I> to 0.0, which means only the viscous drag term in
the Langevin thermostat will be applied. See the discussion on the
<A HREF = "doc/fix_viscous.html">fix viscous</A> doc page for details.
</P>
<P>IMPORTANT NOTE: When the <I>langevin</I> keyword is used with fix rigid
versus fix rigid/small, different dynamics will result for parallel
runs. This is because of the way random numbers are used in the two
cases. The dynamics for the two cases should be statistically
similar, but will not be identical, even for a single timestep.
</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/small</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. The keyword <I>pchain</I> specifies the number of
thermostats in the chain thermostatting the barostat degrees of
freedom.
</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>
<P>The <I>infile</I> keyword allows a file of rigid body attributes to be read
in from a file, rather then having LAMMPS compute them. There are 3
such attributes: the total mass of the rigid body, its center-of-mass
position, and its 6 moments of inertia. For rigid bodies consisting
of point particles or non-overlapping finite-size particles, LAMMPS
can compute these values accurately. However, for rigid bodies
consisting of finite-size particles which overlap each other, LAMMPS
will ignore the overlaps when computing these 3 attributes. The
amount of error this induces depends on the amount of overlap. To
avoid this issue, the values can be pre-computed (e.g. using Monte
Carlo integration).
</P>
<P>The format of the file is as follows. Note that the file does not
have to list attributes for every rigid body integrated by fix rigid.
Only bodies which the file specifies will have their computed
attributes overridden. The file can contain initial blank lines or
comment lines starting with "#" which are ignored. The first
non-blank, non-comment line should list N = the number of lines to
follow. The N successive lines contain the following information:
</P>
<PRE>ID1 masstotal xcm ycm zcm ixx iyy izz ixy ixz iyz
ID2 masstotal xcm ycm zcm ixx iyy izz ixy ixz iyz
...
IDN masstotal xcm ycm zcm ixx iyy izz ixy ixz iyz
</PRE>
<P>The rigid body IDs are all positive integers. For the <I>single</I>
bodystyle, only an ID of 1 can be used. For the <I>group</I> bodystyle,
IDs from 1 to Ng can be used where Ng is the number of specified
groups. For the <I>molecule</I> bodystyle, use the molecule ID for the
atoms in a specific rigid body as the rigid body ID.
</P>
<P>The masstotal and center-of-mass coordinates (xcm,ycm,zcm) are
self-explanatory. The center-of-mass should be consistent with what
is calculated for the position of the rigid body with all its atoms
unwrapped by their respective image flags. If this produces a
center-of-mass that is outside the simulation box, LAMMPS wraps it
back into the box. The 6 moments of inertia (ixx,iyy,izz,ixy,ixz,iyz)
should be the values consistent with the current orientation of the
rigid body around its center of mass. The values are with respect to
the simulation box XYZ axes, not with respect to the prinicpal axes of
the rigid body itself. LAMMPS performs the latter calculation
internally.
</P>
<P>IMPORTANT NOTE: If you use the <I>infile</I> keyword and write restart
files during a simulation, then each time a restart file is written,
the fix also write an auxiliary restart file with the name
rfile.rigid, where "rfile" is the name of the restart file,
e.g. tmp.restart.10000 and tmp.restart.10000.rigid. This auxiliary
file is in the same format described above and contains info on the
current center-of-mass and 6 moments of inertia. Thus it can be used
in a new input script that restarts the run and re-specifies a rigid
fix using an <I>infile</I> keyword and the appropriate filename. Note that
the auxiliary file will contain one line for every rigid body, even if
the original file only listed a subset of the rigid bodies.
</P>
<HR>
<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. Specifically, they are set so
that the center-of-mass (COM) of the rigid body always remains inside
the simulation box.
</P>
<P>This means that if you output per-atom image flags you cannot
interpret them as you normally would. I.e. 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.
</P>
<P>Here are details on how, you can post-process a dump file to calculate
a diffusion coefficient for rigid bodies, using the altered per-atom
image flags written to a dump file. The image flags for atoms in the
same rigid body can be used to unwrap the body and calculate its
center-of-mass (COM). As mentioned above, this COM will always be
inside the simulation box. Thus it will "jump" from one side of the
box to the other when the COM crosses a periodic boundary. If you
keep track of the jumps, you can effectively "unwrap" the COM and use
that value to track the displacement of each rigid body, and thus the
mean-squared displacement (MSD) of an ensemble of bodies, and thus a
diffusion coefficient.
</P>
<P>Note that fix rigid does define image flags for each rigid body, which
are incremented when the center-of-mass of 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>If your simlulation is a hybrid model with a mixture of rigid bodies
and non-rigid particles (e.g. solvent) there are several ways these
rigid fixes can be used in tandem with <A HREF = "fix_nve.html">fix nve</A>, <A HREF = "fix_nh.html">fix
nvt</A>, <A HREF = "fix_nh.html">fix npt</A>, and <A HREF = "fix_nh.html">fix nph</A>.
</P>
<P>If you wish to perform NVE dynamics (no thermostatting or
barostatting), use fix rigid or fix rigid/nve to integrate the rigid
bodies, and <A HREF = "fix_nve.html">fix nve</A> to integrate the non-rigid
particles.
</P>
<P>If you wish to perform NVT dynamics (thermostatting, but no
barostatting), you can use fix rigid/nvt for the rigid bodies, and any
thermostatting fix for the non-rigid particles (<A HREF = "fix_nh.html">fix nvt</A>,
<A HREF = "fix_langevin.html">fix langevin</A>, <A HREF = "fix_temp_berendsen.html">fix
temp/berendsen</A>). You can also use fix rigid
or fix rigid/nve for the rigid bodies and thermostat them using <A HREF = "fix_langevin.html">fix
langevin</A> on the group that contains all the
particles in the rigid bodies. The net force added by <A HREF = "fix_langevin.html">fix
langevin</A> to each rigid body effectively thermostats
its translational center-of-mass motion. Not sure how well it does at
thermostatting its rotational motion.
</P>
<P>If you with to perform NPT or NPH dynamics (barostatting), you cannot
use both <A HREF = "fix_nh.html">fix npt</A> and fix rigid/npt (or the nph
variants). This is because there can only be one fix which monitors
the global pressure and changes the simulation box dimensions. So you
have 3 choices:
</P>
<UL><LI>Use fix rigid/npt for the rigid bodies. Use the <I>dilate</I> all option
so that it will dilate the positions of the non-rigid particles as
well. Use <A HREF = "fix_nh.html">fix nvt</A> (or any other thermostat) for the
non-rigid particles.
<LI>Use <A HREF = "fix_nh.html">fix npt</A> for the group of non-rigid particles. Use
the <I>dilate</I> all option so that it will dilate the center-of-mass
positions of the rigid bodies as well. Use fix rigid/nvt for the
rigid bodies.
<LI>Use <A HREF = "fix_press_berendsen.html">fix press/berendsen</A> to compute the
pressure and change the box dimensions. Use fix rigid/nvt for the
rigid bodies. Use <A HREF = "fix_nh.thml">fix nvt</A> (or any other thermostat) for
the non-rigid particles.
</UL>
<P>In all case, the rigid bodies and non-rigid particles both contribute
to the global pressure and the box is scaled the same by any of the
barostatting fixes.
</P>
<P>You could even use the 2nd and 3rd options for a non-hybrid simulation
consisting of only rigid bodies, assuming you give <A HREF = "fix_nh.html">fix
npt</A> an empty group, though it's an odd thing to do. The
barostatting fixes (<A HREF = "fix_nh.html">fix npt</A> and <A HREF = "fix_press_berendsen.html">fix
press/berensen</A>) will monitor the pressure
and change the box dimensions, but not time integrate any particles.
The integration of the rigid bodies will be performed by fix
rigid/nvt.
</P>
<HR>
<P>Styles with a <I>cuda</I>, <I>gpu</I>, <I>omp</I>, or <I>opt</I> suffix are functionally
the same as the corresponding style without the suffix. They have
been optimized to run faster, depending on your available hardware, as
discussed in <A HREF = "Section_accelerate.html">Section_accelerate</A> of the
manual. The accelerated styles take the same arguments and should
produce the same results, except for round-off and precision issues.
</P>
<P>These accelerated styles are part of the USER-CUDA, GPU, USER-OMP and OPT
packages, respectively. They are only enabled if LAMMPS was built with
those packages. See the <A HREF = "Section_start.html#start_3">Making LAMMPS</A>
section for more info.
</P>
<P>You can specify the accelerated styles explicitly in your input script
by including their suffix, or you can use the <A HREF = "Section_start.html#start_7">-suffix command-line
switch</A> when you invoke LAMMPS, or you can
use the <A HREF = "suffix.html">suffix</A> command in your input script.
</P>
<P>See <A HREF = "Section_accelerate.html">Section_accelerate</A> of the manual for
more instructions on how to use the accelerated styles effectively.
</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/small</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 <A HREF = "fix_modify.html">fix_modify</A> <I>temp</I> and <I>press</I> options are
supported by the rigid/npt and rigid/nph fixes to change the computes used
to calculate the instantaneous pressure tensor. Note that the rigid/nvt fix
does not use any external compute to compute instantaneous temperature.
</P>
<P>The <I>rigid</I> and <I>rigid/small</I> and <I>rigid/nve</I> fixes compute 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, but only for the <I>rigid</I> and
<I>rigid/nve</I> fixes.
</P>
<P>The <I>rigid/nvt</I>, <I>rigid/npt</I>, and <I>rigid/nph</I> fixes compute 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 "extensive". The scalar is the cumulative energy
change due to the thermostatting and barostatting the fix performs.
</P>
<P>All of the <I>rigid</I> fixes except <I>rigid/small</I> 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>
<HR>
<P><B>Restrictions:</B>
</P>
<P>These fixes are all part of the RIGID package. It is only enabled if
LAMMPS was built with that package. See the <A HREF = "Section_start.html#start_3">Making
LAMMPS</A> section for more info.
</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 = Pchain = 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>
Reference in New Issue View Git Blame Copy Permalink