fix ID group-ID style bodystyle args keyword values ...
single args = none molecule args = none group args = N groupID1 groupID2 ... N = # of groups groupID1, groupID2, ... = list of N group IDs
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) temp values = Tstart Tstop Tdamp Tstart,Tstop = desired temperature at start/stop of run (temperature units) Tdamp = temperature damping parameter (time units) iso or aniso values = Pstart Pstop Pdamp Pstart,Pstop = scalar external pressure at start/end of run (pressure units) Pdamp = pressure damping parameter (time units) x or y or z values = Pstart Pstop Pdamp Pstart,Pstop = external stress tensor component at start/end of run (pressure units) Pdamp = stress damping parameter (time units) couple = none or xyz or xy or yz or xz tparam 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 pchain values = Pchain Pchain = length of the Nose/Hoover thermostat chain coupled with the barostat dilate value = dilate-group-ID dilate-group-ID = only dilate atoms in this group due to barostat volume changes force 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 torque 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 infile filename filename = file with per-body values of mass, center-of-mass, moments of inertia mol value = template-ID template-ID = ID of molecule template specified in a separate molecule command
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 fix 1 particles rigid/npt/small molecule temp 1.0 1.0 1.0 iso 0.5 0.5 1.0
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. 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.
Examples of large rigid bodies are a colloidal particle, or portions of a biomolecule such as a protein.
Example of small rigid bodies are patchy nanoparticles, such as those modeled in this paper 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 fix shake 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.
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.
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.
IMPORTANT NOTE: Not all of the bodystyle options and keyword/value options are available for both the rigid and rigid/small variants. See details below.
The rigid 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.
The rigid/small 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.
IMPORTANT NOTE: To use rigid/small 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 pair_style cutoff plus neighbor skin does not span this distance, then you should use the communicate cutoff command with a setting epsilon larger than the distance.
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 may accumulate to produce divergent trajectories.
IMPORTANT NOTE: You should not update the atoms in rigid bodies via other time-integration fixes (e.g. fix nve, fix nvt, fix npt), 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 fix nve or fix nvt should be used for the non-rigid particles.
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 velocity command), setting the force on them to 0.0 (via the fix setforce command), and integrating them as usual (e.g. via the fix nve command).
IMPORTANT NOTE: The aggregate properties of each rigid body are calculated at the start of each simulation run. These include its center of mass, moments of inertia, and net velocity and angular momentum. This means that before or between runs, per-atom properties can be changed, e.g. via the set or velocity command, which will affect the bodies. An exception is if the infile keyword is used, then all the body properties (except net velocity and angular momentum) are only calculated once so that values from the file are valid.
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.
IMPORTANT NOTE: With fix rigid/small, which requires bodystyle molecule, you can define a system that has no rigid bodies initially. This is useful when you are adding rigid bodies on-the-fly via commands such as fix deposit or fix pour.
For bodystyle single the entire fix group of atoms is treated as one rigid body. This option is only allowed for fix rigid and its sub-styles.
For bodystyle molecule, 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.
For bodystyle group, 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.
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 set 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.
The force and torque keywords discussed next are only allowed for fix rigid and its sub-styles.
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 force and torque 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.
An xflag, yflag, or zflag set to off means turn off the component of force of torque in that dimension. A setting of on 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 force and torque 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 force or torque 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.
IMPORTANT NOTE: 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 neigh_modify exclude and delete_bonds commands are used to do this. If the rigid bodies have strongly overalapping atoms, you may need to turn off these interactions to avoid numerical problems due to large equal/opposite intra-body forces swamping the contribution of small inter-body forces.
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.
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 atom_style sphere and ellipsoid and line and tri 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 frictional granular interactions) and have an orientation. These contributions are accounted for by these fixes.
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 neigh_modify exclude and delete_bonds commands are used to do this. For finite-size particles this also means the particles can be highly overlapped when creating the rigid body.
The rigid and rigid/small and rigid/nve styles perform constant NVE time integration. The only difference is that the rigid and rigid/small styles use an integration technique based on Richardson iterations. The rigid/nve style uses the methods described in the paper by Miller, which are thought to provide better energy conservation than an iterative approach.
The rigid/nvt and rigid/nvt/small styles performs constant NVT integration using a Nose/Hoover thermostat with chains as described originally in (Hoover) and (Martyna), which thermostats both the translational and rotational degrees of freedom of the rigid bodies. The rigid-body algorithm used by rigid/nvt is described in the paper by Kamberaj.
The rigid/npt and rigid/nph (and their /small counterparts) 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 rigid/nvt.
The barostat parameters are specified using one or more of the iso, aniso, x, y, z and couple 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 fix npt/nph
NOTE: Currently the rigid/npt and rigid/nph (and their /small counterparts) styles do not support triclinic (non-orthongonal) boxes.
The target pressures for each of the 6 components of the stress tensor can be specified independently via the x, y, z 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 Pstart to Pstop. If a target pressure is specified for a component, then the corresponding box dimension will change during a simulation. For example, if the y 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 fix deform command.
For all barostat keywords, the Pdamp parameter operates like the Tdamp 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 units command).
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 dilate is specified with a dilate-group-ID 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.
The couple 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, xz means the Pxx and Pzz components of the stress tensor are coupled. Xyz 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 Pstart, Pstop, Pdamp parameters for any coupled dimensions must be identical. Couple xyz can be used for a 2d simulation; the z dimension is simply ignored.
The iso and aniso keywords are simply shortcuts that are equivalent to specifying several other keywords together.
The keyword iso 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:
x Pstart Pstop Pdamp y Pstart Pstop Pdamp z Pstart Pstop Pdamp couple xyz
The keyword aniso means x, y, and z dimensions are controlled independently using the Pxx, Pyy, and Pzz 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:
x Pstart Pstop Pdamp y Pstart Pstop Pdamp z Pstart Pstop Pdamp couple none
The keyword/value option pairs are used in the following ways.
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 rigid bodies and how it can be monitored via the fix output is discussed below.
The langevin keyword applies a Langevin thermostat to the constant NVE time integration performed by either the rigid or rigid/small or rigid/nve styles. It cannot be used with the rigid/nvt style. The desired temperature at each timestep is a ramped value during the run from Tstart to Tstop. The Tdamp 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 units command). The random # seed must be a positive integer.
The way that Langevin thermostatting operates is explained on the fix langevin doc page. If you wish to simply viscously damp the rotational motion without thermostatting, you can set Tstart and Tstop to 0.0, which means only the viscous drag term in the Langevin thermostat will be applied. See the discussion on the fix viscous doc page for details.
IMPORTANT NOTE: When the langevin 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.
The temp and tparam keywords apply a Nose/Hoover thermostat to the NVT time integration performed by the rigid/nvt style. They cannot be used with the rigid or rigid/small or rigid/nve styles. The desired temperature at each timestep is a ramped value during the run from Tstart to Tstop. The Tdamp 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 units command).
Nose/Hoover chains are used in conjunction with this thermostat. The tparam keyword can optionally be used to change the chain settings used. Tchain is the number of thermostats in the Nose Hoover chain. This value, along with Tdamp 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 pchain specifies the number of thermostats in the chain thermostatting the barostat degrees of freedom.
IMPORTANT NOTE: There are alternate ways to thermostat a system of rigid bodies. You can use fix langevin 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.
The mol keyword can only be used with fix rigid/small. It should be used when other commands, such as fix deposit or fix pour, add rigid bodies on-the-fly during a simulation. You specify a template-ID previously defined using the molecule command, which reads a file that defines the molecule. You must use the same template-ID that the command adding rigid bodies uses. The coordinates, atom types, atom diameters, center-of-mass, and moments of inertia can be specified in the molecule file. See the molecule command for details. The only settings required to be in this file are the coordinates and types of atoms in the molecule.
The infile 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).
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:
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
The rigid body IDs are all positive integers. For the single bodystyle, only an ID of 1 can be used. For the group bodystyle, IDs from 1 to Ng can be used where Ng is the number of specified groups. For the molecule bodystyle, use the molecule ID for the atoms in a specific rigid body as the rigid body ID.
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.
IMPORTANT NOTE: If you use the infile 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 infile 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.
IMPORTANT NOTE: If you are using fix rigid/small and defining a system that has no rigid bodies initially, because they will be added on-the-fly by commands such as fix deposit or fix pour, you may still wish to use the infile keyword. This is so that restart files written during the simulation will output an auxiliary restart file as described above with information on the new rigid bodies. In this case the initial infile file should use N = 0.
If you use a temperature compute 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.
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).
IMPORTANT NOTE: You may wish to explicitly subtract additional degrees-of-freedom if you use the force and torque keywords to eliminate certain motions of one or more rigid bodies. LAMMPS does not do this automatically.
The rigid body contribution to the pressure of the system (virial) is also accounted for by this fix.
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 fix nve, fix nvt, fix npt, and fix nph.
If you wish to perform NVE dynamics (no thermostatting or barostatting), use fix rigid or fix rigid/nve to integrate the rigid bodies, and fix nve to integrate the non-rigid particles.
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 (fix nvt, fix langevin, fix temp/berendsen). You can also use fix rigid or fix rigid/nve for the rigid bodies and thermostat them using fix langevin on the group that contains all the particles in the rigid bodies. The net force added by fix langevin to each rigid body effectively thermostats its translational center-of-mass motion. Not sure how well it does at thermostatting its rotational motion.
If you with to perform NPT or NPH dynamics (barostatting), you cannot use both fix npt 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:
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.
You could even use the 2nd and 3rd options for a non-hybrid simulation consisting of only rigid bodies, assuming you give fix npt an empty group, though it's an odd thing to do. The barostatting fixes (fix npt and fix press/berensen) 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.
Styles with a cuda, gpu, intel, kk, omp, or opt 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 Section_accelerate of the manual. The accelerated styles take the same arguments and should produce the same results, except for round-off and precision issues.
These accelerated styles are part of the USER-CUDA, GPU, USER-INTEL, KOKKOS, USER-OMP and OPT packages, respectively. They are only enabled if LAMMPS was built with those packages. See the Making LAMMPS section for more info.
You can specify the accelerated styles explicitly in your input script by including their suffix, or you can use the -suffix command-line switch when you invoke LAMMPS, or you can use the suffix command in your input script.
See Section_accelerate of the manual for more instructions on how to use the accelerated styles effectively.
Restart, fix_modify, output, run start/stop, minimize info:
No information about the rigid and rigid/small and rigid/nve fixes are written to binary restart files. For style rigid/nvt the state of the Nose/Hoover thermostat is written to binary restart files. See the read_restart 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.
The fix_modify energy 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 thermodynamic output.
The fix_modify temp and press 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.
The rigid and rigid/small and rigid/nve fixes compute a global scalar which can be accessed by various output commands. 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 force and torque keywords are removed from this calculation, but only for the rigid and rigid/nve fixes.
The rigid/nvt, rigid/npt, and rigid/nph fixes compute a global scalar which can be accessed by various output commands. 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.
All of the rigid fixes except rigid/small compute a global array of values which can be accessed by various output commands. 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.
The center of mass (COM) for each body is similar to unwrapped coordinates written to a dump file. It will always be inside (or slightly outside) the simulation box. The image flags have the same meaning as image flags for atom positions (see the "dump" command). This means you can calculate the unwrapped COM by applying the image flags to the COM, the same as when unwrapped coordinates are written to a dump file.
The force and torque values in the array are not affected by the force and torque keywords in the fix rigid command; they reflect values before any changes are made by those keywords.
The ordering of the rigid bodies (by row in the array) is as follows. For the single keyword there is just one rigid body. For the molecule keyword, the bodies are ordered by ascending molecule ID. For the group keyword, the list of group IDs determines the ordering of bodies.
The array values calculated by these fixes are "intensive", meaning they are independent of the number of atoms in the simulation.
No parameter of these fixes can be used with the start/stop keywords of the run command. These fixes are not invoked during energy minimization.
These fixes are all part of the RIGID package. It is only enabled if LAMMPS was built with that package. See the Making LAMMPS section for more info.
Assigning a temperature via the velocity create command to a system with rigid bodies may not have the desired outcome for two reasons. First, the velocity command can be invoked before the rigid-body fix is invoked or initialized and the number of adjusted degrees of freedom (DOFs) is known. Thus it is not possible to compute the target temperature correctly. Second, the assigned velocities may be partially canceled when constraints are first enforced, leading to a different temperature than desired. A workaround for this is to perform a run 0 command, which insures all DOFs are accounted for properly, and then rescale the temperature to the desired value before performing a simulation. For example:
velocity all create 300.0 12345 run 0 # temperature may not be 300K velocity all scale 300.0 # now it should be
delete_bonds, neigh_modify exclude, fix shake
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.
(Hoover) Hoover, Phys Rev A, 31, 1695 (1985).
(Kamberaj) Kamberaj, Low, Neal, J Chem Phys, 122, 224114 (2005).
(Martyna) Martyna, Klein, Tuckerman, J Chem Phys, 97, 2635 (1992); Martyna, Tuckerman, Tobias, Klein, Mol Phys, 87, 1117.
(Miller) Miller, Eleftheriou, Pattnaik, Ndirango, and Newns, J Chem Phys, 116, 8649 (2002).
(Zhang) Zhang, Glotzer, Nanoletters, 4, 1407-1413 (2004).