bodystyle = single or molecule or group
@@ -75,9 +77,11 @@
Examples:
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
@@ -110,6 +114,41 @@ the constituent atoms as point masses.
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
+summed and exchanged via local communication between neighboring
+processors when ghost atom info is accumlated. This means that the
+ghost atom cutoff be large enough to cover the distance between the
+atom that owns the body and every other atom in the body. If the
+pair_style cutoff is not of sufficient extent to
+insure this, then the communicate cutoff command
+can be used.
+
+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.
+
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
@@ -136,14 +175,22 @@ most one rigid body. Which atoms are in which bodies can be defined
via several options.
For bodystyle single the entire fix group of atoms is treated as one
-rigid body.
+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.
+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.
+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
@@ -157,6 +204,9 @@ the image flags for each atom consistently or that you have used the
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
@@ -192,8 +242,9 @@ to the motion. The neigh_modify exclude and
delete_bonds commands are used to do this.
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.
+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.
@@ -218,12 +269,12 @@ body.
-The rigid and rigid/nve styles perform constant NVE time
-integration. The only difference is that the rigid style uses 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 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 style performs constant NVT integration using a
Nose/Hoover thermostat with chains as described originally in
@@ -232,17 +283,18 @@ 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 styles performs constant NPT or NPH
+
The rigid/npt and rigid/nph 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.
+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
+simulation. The effects of these keywords are similar to those
+defined in fix npt/nph
NOTE: Currently the rigid/npt and rigid/nph styles do not support
triclinic (non-orthongonal) boxes.
@@ -327,25 +379,31 @@ 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/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
+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 the Langevin thermostatting operates is explained on
the fix langevin doc page.
+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/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
+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
@@ -371,8 +429,9 @@ temperature as well without use of the Langevin or Nose/Hoover options
associated with the fix rigid commands.
The infile keyword allows a file of rigid body attributes to be read
-in from a file, rather then letting LAMMPS compute them. There are 3
-such attributes: the total mass of the rigid body, its center-of-mass
+in from a file, rather then letting LAMMPS compute them. It can only
+be used with the fix rigid command and its variants. 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
@@ -382,14 +441,13 @@ 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
+
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:
+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
@@ -525,13 +583,13 @@ rigid/nvt.
Restart, fix_modify, output, run start/stop, minimize info:
-No information about the rigid 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.
+
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
@@ -543,7 +601,10 @@ 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/nve fixes compute a global scalar which can be
+
No global or per-atom quantities are stored by the rigid/small fix
+for access by various output commands.
+
+The rigid 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
@@ -555,24 +616,25 @@ 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.
-The rigid/nvt, rigid/npt and rigid/nph fix compute a global scalar
-which can be accessed by various output
+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 these fixes 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, 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 force and
-torque keywords in the fix rigid command; they reflect values before
-any changes are made by those keywords.
+
All of the rigid fixes (but not the rigid/small fix) 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,
+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 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
@@ -595,11 +657,6 @@ of the run command. These fixes are not invoked during
LAMMPS was built with that package. See the Making
LAMMPS section for more info.
-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.
-
Related commands:
delete_bonds, neigh_modify
diff --git a/doc/fix_rigid.txt b/doc/fix_rigid.txt
index 1001dbcbed..e03e6c15c5 100644
--- a/doc/fix_rigid.txt
+++ b/doc/fix_rigid.txt
@@ -11,13 +11,14 @@ fix rigid/nve command :h3
fix rigid/nvt command :h3
fix rigid/npt command :h3
fix rigid/nph command :h3
+fix rigid/small command :h3
[Syntax:]
fix ID group-ID style bodystyle args keyword values ... :pre
ID, group-ID are documented in "fix"_fix.html command :ulb,l
-style = {rigid} or {rigid/nve} or {rigid/nvt} or {rigid/npt} or {rigid/nph} :l
+style = {rigid} or {rigid/nve} or {rigid/nvt} or {rigid/npt} or {rigid/nph} or {rigid/small} :l
bodystyle = {single} or {molecule} or {group} :l
{single} args = none
{molecule} args = none
@@ -62,9 +63,11 @@ keyword = {langevin} or {temp} or {iso} or {aniso} or {x} or {y} or {z} or {coup
[Examples:]
fix 1 clump rigid single
+fix 1 clump rigid/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
@@ -97,6 +100,41 @@ 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
+summed and exchanged via local communication between neighboring
+processors when ghost atom info is accumlated. This means that the
+ghost atom cutoff be large enough to cover the distance between the
+atom that owns the body and every other atom in the body. If the
+"pair_style"_pair_style.html cutoff is not of sufficient extent to
+insure this, then the "communicate cutoff"_communicate.html command
+can be used.
+
+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.
+
IMPORTANT NOTE: You should not update the atoms in rigid bodies via
other time-integration fixes (e.g. "fix nve"_fix_nve.html, "fix
nvt"_fix_nvt.html, "fix npt"_fix_npt.html), or you will be integrating
@@ -123,14 +161,22 @@ most one rigid body. Which atoms are in which bodies can be defined
via several options.
For bodystyle {single} the entire fix group of atoms is treated as one
-rigid body.
+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.
+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.
+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
@@ -144,6 +190,9 @@ the image flags for each atom consistently or that you have used the
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
@@ -179,8 +228,9 @@ to the motion. The "neigh_modify exclude"_neigh_modify.html and
"delete_bonds"_delete_bonds.html commands are used to do this.
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.
+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.
:line
@@ -205,12 +255,12 @@ body.
:line
-The {rigid} and {rigid/nve} styles perform constant NVE time
-integration. The only difference is that the {rigid} style uses an
-integration technique based on Richardson iterations. The {rigid/nve}
-style uses the methods described in the paper by "Miller"_#Miller,
-which are thought to provide better energy conservation than an
-iterative approach.
+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"_#Miller, which are thought to provide better energy
+conservation than an iterative approach.
The {rigid/nvt} style performs constant NVT integration using a
Nose/Hoover thermostat with chains as described originally in
@@ -219,17 +269,18 @@ 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"_#Kamberaj.
-The {rigid/npt} and {rigid/nph} styles performs constant NPT or NPH
+The {rigid/npt} and {rigid/nph} 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}.
+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"_fix_nh.html
+simulation. The effects of these keywords are similar to those
+defined in "fix npt/nph"_fix_nh.html
NOTE: Currently the {rigid/npt} and {rigid/nph} styles do not support
triclinic (non-orthongonal) boxes.
@@ -314,25 +365,31 @@ 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/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
+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"_units.html command). The random # {seed} must be a positive
integer. The way the Langevin thermostatting operates is explained on
the "fix langevin"_fix_langevin.html doc page.
+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/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
+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"_units.html command).
Nose/Hoover chains are used in conjunction with this thermostat. The
@@ -358,8 +415,9 @@ temperature as well without use of the Langevin or Nose/Hoover options
associated with the fix rigid commands.
The {infile} keyword allows a file of rigid body attributes to be read
-in from a file, rather then letting LAMMPS compute them. There are 3
-such attributes: the total mass of the rigid body, its center-of-mass
+in from a file, rather then letting LAMMPS compute them. It can only
+be used with the fix rigid command and its variants. 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
@@ -369,14 +427,13 @@ 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
+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:
+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
@@ -512,13 +569,13 @@ rigid/nvt.
[Restart, fix_modify, output, run start/stop, minimize info:]
-No information about the {rigid} and {rigid/nve} fixes are written to
-"binary restart files"_restart.html. For style {rigid/nvt} the state
-of the Nose/Hoover thermostat is written to "binary restart
-files"_restart.html. See the "read_restart"_read_restart.html 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.
+No information about the {rigid} and {rigid/small} and {rigid/nve}
+fixes are written to "binary restart files"_restart.html. For style
+{rigid/nvt} the state of the Nose/Hoover thermostat is written to
+"binary restart files"_restart.html. See the
+"read_restart"_read_restart.html 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"_fix_modify.html {energy} option is supported by the
rigid/nvt fix to add the energy change induced by the thermostatting
@@ -530,7 +587,10 @@ 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/nve fixes compute a global scalar which can be
+No global or per-atom quantities are stored by the {rigid/small} fix
+for access by various "output commands"_Section_howto.html#howto_15.
+
+The {rigid} and {rigid/nve} fixes compute a global scalar which can be
accessed by various "output commands"_Section_howto.html#howto_15.
The scalar value calculated by these fixes is "intensive". The scalar
is the current temperature of the collection of rigid bodies. This is
@@ -542,24 +602,25 @@ 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.
-The rigid/nvt, rigid/npt and rigid/nph fix compute a global scalar
-which can be accessed by various "output
+The {rigid/nvt}, {rigid/npt}, and {rigid/nph} fixes compute a global
+scalar which can be accessed by various "output
commands"_Section_howto.html#howto_15. 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 these fixes compute a global array of values which can be
-accessed by various "output commands"_Section_howto.html#howto_15.
-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 {force} and
-{torque} keywords in the fix rigid command; they reflect values before
-any changes are made by those keywords.
+All of the {rigid} fixes (but not the {rigid/small} fix) compute a
+global array of values which can be accessed by various "output
+commands"_Section_howto.html#howto_15. 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 {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
@@ -582,11 +643,6 @@ 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_start.html#start_3 section for more info.
-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.
-
[Related commands:]
"delete_bonds"_delete_bonds.html, "neigh_modify"_neigh_modify.html
diff --git a/doc/kspace_modify.html b/doc/kspace_modify.html
index 3da9a66b4f..8138692d0b 100644
--- a/doc/kspace_modify.html
+++ b/doc/kspace_modify.html
@@ -41,6 +41,7 @@
nozforce turns off kspace forces in the z direction
compute value = yes or no
cutoff/adjust value = yes or no
+ fftbench value = yes or no
diff value = ad or ik = 2 or 4 FFTs for PPPM in smoothed or non-smoothed mode
@@ -191,6 +192,11 @@ grid that minimizes cost using an estimate given by
(Hardy). Note that this cost estimate is not exact, somewhat
experimental, and still may not yield the optimal parameters.
+The fftbench keyword applies only to PPPM. It is on by default. If
+this option is turned off, LAMMPS will not take the time at the end
+of a run to give FFT benchmark timings, and will finish a few seconds
+faster than it would if this option were on.
+
The diff keyword specifies the differentiation scheme used by the
PPPM method to compute forces on particles given electrostatic
potentials on the PPPM mesh. The ik approach is the default for
@@ -221,7 +227,7 @@ option. Support for those PPPM variants will be added later.
The option defaults are mesh = mesh/disp = 0 0 0, order = order/disp =
5 (PPPM), order = 8 (MSM), minorder = 2, overlap = yes, force = -1.0,
gewald = gewald/disp = 0.0, slab = 1.0, compute = yes, cutoff/adjust =
-yes (MSM), and diff = ik (PPPM).
+yes (MSM), fftbench = yes (PPPM), and diff = ik (PPPM).
diff --git a/doc/kspace_modify.txt b/doc/kspace_modify.txt
index 3a21672458..b232f5f22b 100644
--- a/doc/kspace_modify.txt
+++ b/doc/kspace_modify.txt
@@ -36,6 +36,7 @@ keyword = {mesh} or {order} or {order/disp} or {overlap} or {minorder} or {force
{nozforce} turns off kspace forces in the z direction
{compute} value = {yes} or {no}
{cutoff/adjust} value = {yes} or {no}
+ {fftbench} value = {yes} or {no}
{diff} value = {ad} or {ik} = 2 or 4 FFTs for PPPM in smoothed or non-smoothed mode :pre
:ule
@@ -185,6 +186,11 @@ grid that minimizes cost using an estimate given by
"(Hardy)"_#Hardy. Note that this cost estimate is not exact, somewhat
experimental, and still may not yield the optimal parameters.
+The {fftbench} keyword applies only to PPPM. It is on by default. If
+this option is turned off, LAMMPS will not take the time at the end
+of a run to give FFT benchmark timings, and will finish a few seconds
+faster than it would if this option were on.
+
The {diff} keyword specifies the differentiation scheme used by the
PPPM method to compute forces on particles given electrostatic
potentials on the PPPM mesh. The {ik} approach is the default for
@@ -215,7 +221,7 @@ option. Support for those PPPM variants will be added later.
The option defaults are mesh = mesh/disp = 0 0 0, order = order/disp =
5 (PPPM), order = 8 (MSM), minorder = 2, overlap = yes, force = -1.0,
gewald = gewald/disp = 0.0, slab = 1.0, compute = yes, cutoff/adjust =
-yes (MSM), and diff = ik (PPPM).
+yes (MSM), fftbench = yes (PPPM), and diff = ik (PPPM).
:line
diff --git a/src/ASPHERE/compute_temp_asphere.cpp b/src/ASPHERE/compute_temp_asphere.cpp
index b61aa0607d..0bad2d713f 100644
--- a/src/ASPHERE/compute_temp_asphere.cpp
+++ b/src/ASPHERE/compute_temp_asphere.cpp
@@ -56,14 +56,16 @@ ComputeTempAsphere::ComputeTempAsphere(LAMMPS *lmp, int narg, char **arg) :
int iarg = 3;
while (iarg < narg) {
if (strcmp(arg[iarg],"bias") == 0) {
- if (iarg+2 > narg) error->all(FLERR,"Illegal compute temp/asphere command");
+ if (iarg+2 > narg)
+ error->all(FLERR,"Illegal compute temp/asphere command");
tempbias = 1;
int n = strlen(arg[iarg+1]) + 1;
id_bias = new char[n];
strcpy(id_bias,arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"dof") == 0) {
- if (iarg+2 > narg) error->all(FLERR,"Illegal compute temp/asphere command");
+ if (iarg+2 > narg)
+ error->all(FLERR,"Illegal compute temp/asphere command");
if (strcmp(arg[iarg+1],"rotate") == 0) mode = ROTATE;
else if (strcmp(arg[iarg+1],"all") == 0) mode = ALL;
else error->all(FLERR,"Illegal compute temp/asphere command");
@@ -106,7 +108,8 @@ void ComputeTempAsphere::init()
if (tempbias) {
int i = modify->find_compute(id_bias);
- if (i < 0) error->all(FLERR,"Could not find compute ID for temperature bias");
+ if (i < 0)
+ error->all(FLERR,"Could not find compute ID for temperature bias");
tbias = modify->compute[i];
if (tbias->tempflag == 0)
error->all(FLERR,"Bias compute does not calculate temperature");
@@ -115,10 +118,16 @@ void ComputeTempAsphere::init()
if (tbias->igroup != igroup)
error->all(FLERR,"Bias compute group does not match compute group");
tbias->init();
+ tbias->setup();
if (strcmp(tbias->style,"temp/region") == 0) tempbias = 2;
else tempbias = 1;
}
+}
+/* ---------------------------------------------------------------------- */
+
+void ComputeTempAsphere::setup()
+{
fix_dof = 0;
for (int i = 0; i < modify->nfix; i++)
fix_dof += modify->fix[i]->dof(igroup);
diff --git a/src/ASPHERE/compute_temp_asphere.h b/src/ASPHERE/compute_temp_asphere.h
index 1e393568cf..b746743823 100644
--- a/src/ASPHERE/compute_temp_asphere.h
+++ b/src/ASPHERE/compute_temp_asphere.h
@@ -29,6 +29,7 @@ class ComputeTempAsphere : public Compute {
ComputeTempAsphere(class LAMMPS *, int, char **);
~ComputeTempAsphere();
void init();
+ void setup();
double compute_scalar();
void compute_vector();
diff --git a/src/RIGID/Install.sh b/src/RIGID/Install.sh
index 3979831f8e..89c951ce03 100644
--- a/src/RIGID/Install.sh
+++ b/src/RIGID/Install.sh
@@ -8,6 +8,7 @@ if (test $1 = 1) then
cp fix_rigid_npt.cpp ..
cp fix_rigid_nve.cpp ..
cp fix_rigid_nvt.cpp ..
+ cp fix_rigid_small.cpp ..
cp fix_rigid.h ..
cp fix_rigid_nh.h ..
@@ -15,6 +16,7 @@ if (test $1 = 1) then
cp fix_rigid_npt.h ..
cp fix_rigid_nve.h ..
cp fix_rigid_nvt.h ..
+ cp fix_rigid_small.h ..
elif (test $1 = 0) then
@@ -24,6 +26,7 @@ elif (test $1 = 0) then
rm -f ../fix_rigid_npt.cpp
rm -f ../fix_rigid_nve.cpp
rm -f ../fix_rigid_nvt.cpp
+ rm -f ../fix_rigid_small.cpp
rm -f ../fix_rigid.h
rm -f ../fix_rigid_nh.h
@@ -31,6 +34,7 @@ elif (test $1 = 0) then
rm -f ../fix_rigid_npt.h
rm -f ../fix_rigid_nve.h
rm -f ../fix_rigid_nvt.h
+ rm -f ../fix_rigid_small.h
fi
diff --git a/src/RIGID/fix_rigid.cpp b/src/RIGID/fix_rigid.cpp
index b6fdbf56e3..6a146818f3 100644
--- a/src/RIGID/fix_rigid.cpp
+++ b/src/RIGID/fix_rigid.cpp
@@ -1131,12 +1131,13 @@ void FixRigid::pre_neighbor()
}
/* ----------------------------------------------------------------------
- count # of degrees-of-freedom removed by rigid bodies for atoms in igroup
+ count # of DOF removed by rigid bodies for atoms in igroup
+ return total count of DOF
------------------------------------------------------------------------- */
-int FixRigid::dof(int igroup)
+int FixRigid::dof(int tgroup)
{
- int groupbit = group->bitmask[igroup];
+ int tgroupbit = group->bitmask[tgroup];
// nall = # of point particles in each rigid body
// mall = # of finite-size particles in each rigid body
@@ -1151,7 +1152,7 @@ int FixRigid::dof(int igroup)
ncount[ibody] = mcount[ibody] = 0;
for (int i = 0; i < nlocal; i++)
- if (body[i] >= 0 && mask[i] & groupbit) {
+ if (body[i] >= 0 && mask[i] & tgroupbit) {
if (extended && eflags[i]) mcount[body[i]]++;
else ncount[body[i]]++;
}
diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp
new file mode 100644
index 0000000000..261eeef816
--- /dev/null
+++ b/src/RIGID/fix_rigid_small.cpp
@@ -0,0 +1,2550 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "stdio.h"
+#include "stdlib.h"
+#include "string.h"
+#include "fix_rigid_small.h"
+#include "math_extra.h"
+#include "atom.h"
+#include "atom_vec_ellipsoid.h"
+#include "atom_vec_line.h"
+#include "atom_vec_tri.h"
+#include "domain.h"
+#include "update.h"
+#include "respa.h"
+#include "modify.h"
+#include "group.h"
+#include "comm.h"
+#include "force.h"
+#include "output.h"
+#include "random_mars.h"
+#include "math_const.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+using namespace MathConst;
+
+// allocate space for static class variable
+
+FixRigidSmall *FixRigidSmall::frsptr;
+
+#define MAXLINE 256
+#define CHUNK 1024
+#define ATTRIBUTE_PERBODY 11
+
+#define TOLERANCE 1.0e-6
+#define EPSILON 1.0e-7
+#define BIG 1.0e20
+
+#define SINERTIA 0.4 // moment of inertia prefactor for sphere
+#define EINERTIA 0.4 // moment of inertia prefactor for ellipsoid
+#define LINERTIA (1.0/12.0) // moment of inertia prefactor for line segment
+
+#define DELTA_BODY 10000
+
+enum{FULL_BODY,INITIAL,FINAL,FORCE_TORQUE,VCM_ANGMOM,XCM_MASS,ITENSOR,DOF};
+
+/* ---------------------------------------------------------------------- */
+
+FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
+ Fix(lmp, narg, arg)
+{
+ int i,ibody;
+
+ scalar_flag = 1;
+ extscalar = 0;
+ time_integrate = 1;
+ rigid_flag = 1;
+ virial_flag = 1;
+ create_attribute = 1;
+
+
+ MPI_Comm_rank(world,&me);
+ MPI_Comm_size(world,&nprocs);
+
+ // perform initial allocation of atom-based arrays
+ // register with Atom class
+
+ extended = orientflag = dorientflag = 0;
+ bodyown = NULL;
+ bodytag = NULL;
+ atom2body = NULL;
+ displace = NULL;
+ eflags = NULL;
+ orient = NULL;
+ dorient = NULL;
+ grow_arrays(atom->nmax);
+ atom->add_callback(0);
+
+ // parse args for rigid body specification
+
+ if (narg < 4) error->all(FLERR,"Illegal fix rigid/small command");
+ if (strcmp(arg[3],"molecule") != 0)
+ error->all(FLERR,"Illegal fix rigid/small command");
+
+ // parse optional args
+
+ int seed;
+ langflag = 0;
+
+ int iarg = 4;
+ while (iarg < narg) {
+ if (strcmp(arg[iarg],"langevin") == 0) {
+ if (iarg+5 > narg) error->all(FLERR,"Illegal fix rigid/small command");
+ if (strcmp(style,"rigid/small") != 0)
+ error->all(FLERR,"Illegal fix rigid/small command");
+ langflag = 1;
+ t_start = atof(arg[iarg+1]);
+ t_stop = atof(arg[iarg+2]);
+ t_period = atof(arg[iarg+3]);
+ seed = atoi(arg[iarg+4]);
+ if (t_period <= 0.0)
+ error->all(FLERR,"Fix rigid/small langevin period must be > 0.0");
+ if (seed <= 0) error->all(FLERR,"Illegal fix rigid/small command");
+ iarg += 5;
+ } else error->all(FLERR,"Illegal fix rigid/small command");
+ }
+
+ // create rigid bodies based on molecule ID
+ // sets bodytag for owned atoms
+ // body attributes are computed later by setup_bodies()
+
+ if (atom->molecule_flag == 0)
+ error->all(FLERR,"Fix rigid/small requires atom attribute molecule");
+
+ create_bodies();
+
+ // set nlocal_body and allocate bodies I own
+
+ int *tag = atom->tag;
+ int nlocal = atom->nlocal;
+
+ nlocal_body = nghost_body = 0;
+ for (i = 0; i < nlocal; i++)
+ if (bodytag[i] == tag[i]) nlocal_body++;
+
+ nmax_body = 0;
+ while (nmax_body < nlocal_body) nmax_body += DELTA_BODY;
+ body = (Body *) memory->smalloc(nmax_body*sizeof(Body),
+ "rigid/small:body");
+
+ // set bodyown for owned atoms
+
+ nlocal_body = 0;
+ for (i = 0; i < nlocal; i++)
+ if (bodytag[i] == tag[i]) {
+ body[nlocal_body].ilocal = i;
+ bodyown[i] = nlocal_body++;
+ } else bodyown[i] = -1;
+
+ // bodysize = sizeof(Body) in doubles
+
+ bodysize = sizeof(Body)/sizeof(double);
+ if (bodysize*sizeof(double) != sizeof(Body)) bodysize++;
+
+ // set max comm sizes needed by this fix
+
+ comm_forward = 1 + bodysize;
+ comm_reverse = 6;
+
+ // bitmasks for properties of extended particles
+
+ POINT = 1;
+ SPHERE = 2;
+ ELLIPSOID = 4;
+ LINE = 8;
+ TRIANGLE = 16;
+ DIPOLE = 32;
+ OMEGA = 64;
+ ANGMOM = 128;
+ TORQUE = 256;
+
+ MINUSPI = -MY_PI;
+ TWOPI = 2.0*MY_PI;
+
+ // atom style pointers to particles that store extra info
+
+ avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
+ avec_line = (AtomVecLine *) atom->style_match("line");
+ avec_tri = (AtomVecTri *) atom->style_match("tri");
+
+ // print statistics
+
+ int one = 0;
+ bigint atomone = 0;
+ for (int i = 0; i < nlocal; i++) {
+ if (bodyown[i] >= 0) one++;
+ if (bodytag[i] > 0) atomone++;
+ }
+ MPI_Allreduce(&one,&nbody,1,MPI_INT,MPI_SUM,world);
+ bigint atomall;
+ MPI_Allreduce(&atomone,&atomall,1,MPI_LMP_BIGINT,MPI_SUM,world);
+ if (nbody == 0) error->all(FLERR,"No rigid bodies defined");
+
+ if (me == 0) {
+ if (screen) fprintf(screen,
+ "%d rigid bodies with " BIGINT_FORMAT " atoms\n",
+ nbody,atomall);
+ if (logfile) fprintf(logfile,
+ "%d rigid bodies with " BIGINT_FORMAT " atoms\n",
+ nbody,atomall);
+ }
+
+ // initialize Marsaglia RNG with processor-unique seed
+
+ maxlang = 0;
+ langextra = NULL;
+ random = NULL;
+ if (langflag) random = new RanMars(lmp,seed + comm->me);
+
+ // firstflag = 1 triggers one-time initialization of rigid body attributes
+
+ firstflag = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixRigidSmall::~FixRigidSmall()
+{
+ // unregister callbacks to this fix from Atom class
+
+ atom->delete_callback(id,0);
+
+ // delete locally stored arrays
+
+ memory->sfree(body);
+
+ memory->destroy(bodyown);
+ memory->destroy(bodytag);
+ memory->destroy(atom2body);
+ memory->destroy(displace);
+ memory->destroy(eflags);
+ memory->destroy(orient);
+ memory->destroy(dorient);
+
+ delete random;
+ memory->destroy(langextra);
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixRigidSmall::setmask()
+{
+ int mask = 0;
+ mask |= INITIAL_INTEGRATE;
+ mask |= FINAL_INTEGRATE;
+ if (langflag) mask |= POST_FORCE;
+ mask |= PRE_NEIGHBOR;
+ mask |= INITIAL_INTEGRATE_RESPA;
+ mask |= FINAL_INTEGRATE_RESPA;
+ return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixRigidSmall::init()
+{
+ int i,ibody;
+
+ triclinic = domain->triclinic;
+
+ // warn if more than one rigid fix
+
+ int count = 0;
+ for (i = 0; i < modify->nfix; i++)
+ if (strcmp(modify->fix[i]->style,"rigid") == 0) count++;
+ if (count > 1 && me == 0) error->warning(FLERR,"More than one fix rigid");
+
+ // error if npt,nph fix comes before rigid fix
+
+ for (i = 0; i < modify->nfix; i++) {
+ if (strcmp(modify->fix[i]->style,"npt") == 0) break;
+ if (strcmp(modify->fix[i]->style,"nph") == 0) break;
+ }
+ if (i < modify->nfix) {
+ for (int j = i; j < modify->nfix; j++)
+ if (strcmp(modify->fix[j]->style,"rigid") == 0)
+ error->all(FLERR,"Rigid fix must come before NPT/NPH fix");
+ }
+
+ // timestep info
+
+ dtv = update->dt;
+ dtf = 0.5 * update->dt * force->ftm2v;
+ dtq = 0.5 * update->dt;
+
+ if (strstr(update->integrate_style,"respa"))
+ step_respa = ((Respa *) update->integrate)->step;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixRigidSmall::setup(int vflag)
+{
+ int i,n,ibody;
+ double massone,radone;
+
+ // one-time initialization of rigid body attributes via local comm
+ // extended flags, mass, COM, inertia tensor, displacement of each atom
+
+ if (firstflag) {
+ firstflag = 0;
+ setup_bodies();
+ }
+
+ //check(1);
+
+ // sum vcm, fcm, angmom, torque across all rigid bodies
+ // vcm = velocity of COM
+ // fcm = force on COM
+ // angmom = angular momentum around COM
+ // torque = torque around COM
+
+ double **x = atom->x;
+ double **v = atom->v;
+ double **f = atom->f;
+ double *rmass = atom->rmass;
+ double *mass = atom->mass;
+ int *type = atom->type;
+ int *image = atom->image;
+ int nlocal = atom->nlocal;
+
+ double *xcm,*vcm,*fcm,*acm,*tcm;
+ double dx,dy,dz;
+ double unwrap[3];
+
+ for (ibody = 0; ibody < nlocal_body+nghost_body; ibody++) {
+ vcm = body[ibody].vcm;
+ vcm[0] = vcm[1] = vcm[2] = 0.0;
+ fcm = body[ibody].fcm;
+ fcm[0] = fcm[1] = fcm[2] = 0.0;
+ acm = body[ibody].angmom;
+ acm[0] = acm[1] = acm[2] = 0.0;
+ tcm = body[ibody].torque;
+ tcm[0] = tcm[1] = tcm[2] = 0.0;
+ }
+
+ for (i = 0; i < nlocal; i++) {
+ if (atom2body[i] < 0) continue;
+ Body *b = &body[atom2body[i]];
+
+ if (rmass) massone = rmass[i];
+ else massone = mass[type[i]];
+
+ vcm = b->vcm;
+ vcm[0] += v[i][0] * massone;
+ vcm[1] += v[i][1] * massone;
+ vcm[2] += v[i][2] * massone;
+ fcm = b->fcm;
+ fcm[0] += f[i][0];
+ fcm[1] += f[i][1];
+ fcm[2] += f[i][2];
+
+ domain->unmap(x[i],image[i],unwrap);
+ xcm = b->xcm;
+ dx = unwrap[0] - xcm[0];
+ dy = unwrap[1] - xcm[1];
+ dz = unwrap[2] - xcm[2];
+
+ acm = b->angmom;
+ acm[0] += dy * massone*v[i][2] - dz * massone*v[i][1];
+ acm[1] += dz * massone*v[i][0] - dx * massone*v[i][2];
+ acm[2] += dx * massone*v[i][1] - dy * massone*v[i][0];
+ tcm = b->torque;
+ tcm[0] += dy * f[i][2] - dz * f[i][1];
+ tcm[1] += dz * f[i][0] - dx * f[i][2];
+ tcm[2] += dx * f[i][1] - dy * f[i][0];
+ }
+
+ // extended particles add their rotation/torque to angmom/torque of body
+
+ if (extended) {
+ AtomVecLine::Bonus *lbonus;
+ if (avec_line) lbonus = avec_line->bonus;
+ double **omega = atom->omega;
+ double **angmom = atom->angmom;
+ double **torque = atom->torque;
+ double *radius = atom->radius;
+ int *line = atom->line;
+
+ for (i = 0; i < nlocal; i++) {
+ if (atom2body[i] < 0) continue;
+ Body *b = &body[atom2body[i]];
+
+ if (eflags[i] & OMEGA) {
+ if (eflags[i] & SPHERE) {
+ radone = radius[i];
+ acm = b->angmom;
+ acm[0] += SINERTIA*rmass[i] * radone*radone * omega[i][0];
+ acm[1] += SINERTIA*rmass[i] * radone*radone * omega[i][1];
+ acm[2] += SINERTIA*rmass[i] * radone*radone * omega[i][2];
+ } else if (eflags[i] & LINE) {
+ radone = lbonus[line[i]].length;
+ b->angmom[2] += LINERTIA*rmass[i] * radone*radone * omega[i][2];
+ }
+ }
+ if (eflags[i] & ANGMOM) {
+ acm = b->angmom;
+ acm[0] += angmom[i][0];
+ acm[1] += angmom[i][1];
+ acm[2] += angmom[i][2];
+ }
+ if (eflags[i] & TORQUE) {
+ tcm = b->torque;
+ tcm[0] += torque[i][0];
+ tcm[1] += torque[i][1];
+ tcm[2] += torque[i][2];
+ }
+ }
+ }
+
+ // reverse communicate vcm, fcm, angmom, torque of all bodies
+
+ commflag = VCM_ANGMOM;
+ comm->reverse_comm_variable_fix(this);
+ commflag = FORCE_TORQUE;
+ comm->reverse_comm_variable_fix(this);
+
+ // normalize velocity of COM
+
+ for (ibody = 0; ibody < nlocal_body; ibody++) {
+ vcm = body[ibody].vcm;
+ vcm[0] /= body[ibody].mass;
+ vcm[1] /= body[ibody].mass;
+ vcm[2] /= body[ibody].mass;
+ }
+
+ // virial setup before call to set_v
+
+ if (vflag) v_setup(vflag);
+ else evflag = 0;
+
+ // compute and forward communicate updated info of all bodies
+
+ for (ibody = 0; ibody < nlocal_body; ibody++) {
+ Body *b = &body[ibody];
+ MathExtra::angmom_to_omega(b->angmom,b->ex_space,b->ey_space,
+ b->ez_space,b->inertia,b->omega);
+ }
+
+ commflag = FINAL;
+ comm->forward_comm_variable_fix(this);
+
+ // set velocity/rotation of atoms in rigid bodues
+
+ set_v();
+
+ // guesstimate virial as 2x the set_v contribution
+
+ if (vflag_global)
+ for (n = 0; n < 6; n++) virial[n] *= 2.0;
+ if (vflag_atom) {
+ for (i = 0; i < nlocal; i++)
+ for (n = 0; n < 6; n++)
+ vatom[i][n] *= 2.0;
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixRigidSmall::initial_integrate(int vflag)
+{
+ double dtfm;
+
+ //check(2);
+
+ for (int ibody = 0; ibody < nlocal_body; ibody++) {
+ Body *b = &body[ibody];
+
+ // update vcm by 1/2 step
+
+ dtfm = dtf / b->mass;
+ b->vcm[0] += dtfm * b->fcm[0];
+ b->vcm[1] += dtfm * b->fcm[1];
+ b->vcm[2] += dtfm * b->fcm[2];
+
+ // update xcm by full step
+
+ b->xcm[0] += dtv * b->vcm[0];
+ b->xcm[1] += dtv * b->vcm[1];
+ b->xcm[2] += dtv * b->vcm[2];
+
+ // update angular momentum by 1/2 step
+
+ b->angmom[0] += dtf * b->torque[0];
+ b->angmom[1] += dtf * b->torque[1];
+ b->angmom[2] += dtf * b->torque[2];
+
+ // compute omega at 1/2 step from angmom at 1/2 step and current q
+ // update quaternion a full step via Richardson iteration
+ // returns new normalized quaternion, also updated omega at 1/2 step
+ // update ex,ey,ez to reflect new quaternion
+
+ MathExtra::angmom_to_omega(b->angmom,b->ex_space,b->ey_space,
+ b->ez_space,b->inertia,b->omega);
+ MathExtra::richardson(b->quat,b->angmom,b->omega,b->inertia,dtq);
+ MathExtra::q_to_exyz(b->quat,b->ex_space,b->ey_space,b->ez_space);
+ }
+
+ // virial setup before call to set_xv
+
+ if (vflag) v_setup(vflag);
+ else evflag = 0;
+
+ // forward communicate updated info of all bodies
+
+ commflag = INITIAL;
+ comm->forward_comm_variable_fix(this);
+
+ // set coords/orient and velocity/rotation of atoms in rigid bodies
+
+ set_xv();
+}
+
+/* ----------------------------------------------------------------------
+ apply Langevin thermostat to all 6 DOF of rigid bodies I own
+ unlike fix langevin, this stores extra force in extra arrays,
+ which are added in when final_integrate() calculates a new fcm/torque
+------------------------------------------------------------------------- */
+
+void FixRigidSmall::post_force(int vflag)
+{
+ double gamma1,gamma2;
+
+ // grow langextra if needed
+
+ if (nlocal_body > maxlang) {
+ memory->destroy(langextra);
+ maxlang = nlocal_body + nghost_body;
+ memory->create(langextra,maxlang,6,"rigid/small:langextra");
+ }
+
+ double delta = update->ntimestep - update->beginstep;
+ delta /= update->endstep - update->beginstep;
+ double t_target = t_start + delta * (t_stop-t_start);
+ double tsqrt = sqrt(t_target);
+
+ double boltz = force->boltz;
+ double dt = update->dt;
+ double mvv2e = force->mvv2e;
+ double ftm2v = force->ftm2v;
+
+ double *vcm,*omega,*inertia;
+
+ for (int ibody = 0; ibody < nlocal_body; ibody++) {
+ vcm = body[ibody].vcm;
+ omega = body[ibody].omega;
+ inertia = body[ibody].inertia;
+
+ gamma1 = -body[ibody].mass / t_period / ftm2v;
+ gamma2 = sqrt(body[ibody].mass) * tsqrt *
+ sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
+ langextra[ibody][0] = gamma1*vcm[0] + gamma2*(random->uniform()-0.5);
+ langextra[ibody][1] = gamma1*vcm[1] + gamma2*(random->uniform()-0.5);
+ langextra[ibody][2] = gamma1*vcm[2] + gamma2*(random->uniform()-0.5);
+
+ gamma1 = -1.0 / t_period / ftm2v;
+ gamma2 = tsqrt * sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
+ langextra[ibody][3] = inertia[0]*gamma1*omega[0] +
+ sqrt(inertia[0])*gamma2*(random->uniform()-0.5);
+ langextra[ibody][4] = inertia[1]*gamma1*omega[1] +
+ sqrt(inertia[1])*gamma2*(random->uniform()-0.5);
+ langextra[ibody][5] = inertia[2]*gamma1*omega[2] +
+ sqrt(inertia[2])*gamma2*(random->uniform()-0.5);
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixRigidSmall::final_integrate()
+{
+ int i,ibody;
+ double dtfm,xy,xz,yz;
+
+ //check(3);
+
+ // sum over atoms to get force and torque on rigid body
+
+ tagint *image = atom->image;
+ double **x = atom->x;
+ double **f = atom->f;
+ int nlocal = atom->nlocal;
+
+ double dx,dy,dz;
+ double unwrap[3];
+ double *xcm,*fcm,*tcm;
+
+ for (ibody = 0; ibody < nlocal_body+nghost_body; ibody++) {
+ fcm = body[ibody].fcm;
+ fcm[0] = fcm[1] = fcm[2] = 0.0;
+ tcm = body[ibody].torque;
+ tcm[0] = tcm[1] = tcm[2] = 0.0;
+ }
+
+ for (i = 0; i < nlocal; i++) {
+ if (atom2body[i] < 0) continue;
+ Body *b = &body[atom2body[i]];
+
+ fcm = b->fcm;
+ fcm[0] += f[i][0];
+ fcm[1] += f[i][1];
+ fcm[2] += f[i][2];
+
+ domain->unmap(x[i],image[i],unwrap);
+ xcm = b->xcm;
+ dx = unwrap[0] - xcm[0];
+ dy = unwrap[1] - xcm[1];
+ dz = unwrap[2] - xcm[2];
+
+ tcm = b->torque;
+ tcm[0] += dy*f[i][2] - dz*f[i][1];
+ tcm[1] += dz*f[i][0] - dx*f[i][2];
+ tcm[2] += dx*f[i][1] - dy*f[i][0];
+ }
+
+ // extended particles add their torque to torque of body
+
+ if (extended) {
+ double **torque = atom->torque;
+
+ for (i = 0; i < nlocal; i++) {
+ if (atom2body[i] < 0) continue;
+
+ if (eflags[i] & TORQUE) {
+ tcm = body[atom2body[i]].torque;
+ tcm[0] += torque[i][0];
+ tcm[1] += torque[i][1];
+ tcm[2] += torque[i][2];
+ }
+ }
+ }
+
+ // reverse communicate fcm, torque of all bodies
+
+ commflag = FORCE_TORQUE;
+ comm->reverse_comm_variable_fix(this);
+
+ // include Langevin thermostat forces and torques
+
+ if (langflag) {
+ for (int ibody = 0; ibody < nlocal_body; ibody++) {
+ fcm = body[atom2body[i]].fcm;
+ fcm[0] += langextra[ibody][0];
+ fcm[1] += langextra[ibody][1];
+ fcm[2] += langextra[ibody][2];
+ tcm = body[atom2body[i]].torque;
+ tcm[0] += langextra[ibody][3];
+ tcm[1] += langextra[ibody][4];
+ tcm[2] += langextra[ibody][5];
+ }
+ }
+
+ // update vcm and angmom, recompute omega
+
+ for (int ibody = 0; ibody < nlocal_body; ibody++) {
+ Body *b = &body[ibody];
+
+ // update vcm by 1/2 step
+
+ dtfm = dtf / b->mass;
+ b->vcm[0] += dtfm * b->fcm[0];
+ b->vcm[1] += dtfm * b->fcm[1];
+ b->vcm[2] += dtfm * b->fcm[2];
+
+ // update angular momentum by 1/2 step
+
+ b->angmom[0] += dtf * b->torque[0];
+ b->angmom[1] += dtf * b->torque[1];
+ b->angmom[2] += dtf * b->torque[2];
+
+ MathExtra::angmom_to_omega(b->angmom,b->ex_space,b->ey_space,
+ b->ez_space,b->inertia,b->omega);
+ }
+
+ // forward communicate updated info of all bodies
+
+ commflag = FINAL;
+ comm->forward_comm_variable_fix(this);
+
+ // set velocity/rotation of atoms in rigid bodies
+ // virial is already setup from initial_integrate
+
+ set_v();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixRigidSmall::initial_integrate_respa(int vflag, int ilevel, int iloop)
+{
+ dtv = step_respa[ilevel];
+ dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
+ dtq = 0.5 * step_respa[ilevel];
+
+ if (ilevel == 0) initial_integrate(vflag);
+ else final_integrate();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixRigidSmall::final_integrate_respa(int ilevel, int iloop)
+{
+ dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
+ final_integrate();
+}
+
+/* ----------------------------------------------------------------------
+ atom to processor assignments have changed,
+ so acquire ghost bodies and setup atom2body
+ remap xcm of each rigid body back into periodic simulation box
+ done during pre_neighbor so will be after call to pbc()
+ and after fix_deform::pre_exchange() may have flipped box
+ use domain->remap() in case xcm is far away from box
+ due to 1st definition of rigid body or due to box flip
+ if don't do this, then atoms of a body which drifts far away
+ from a triclinic box will be remapped back into box
+ with huge displacements when the box tilt changes via set_x()
+ adjust image flag of body and image flags of all atoms in body
+------------------------------------------------------------------------- */
+
+void FixRigidSmall::pre_neighbor()
+{
+ // do nothing if bodies are not yet initialized
+ // i.e. pre_neighbor() is being called from Verlet::setup()
+
+ if (firstflag) return;
+
+ // remap xcm and image flags of each body as needed
+
+ tagint original,oldimage,newimage;
+
+ for (int ibody = 0; ibody < nlocal_body; ibody++) {
+ Body *b = &body[ibody];
+
+ original = b->image;
+ domain->remap(b->xcm,b->image);
+
+ if (original == b->image) b->remapflag[3] = 0;
+ else {
+ oldimage = original & IMGMASK;
+ newimage = b->image & IMGMASK;
+ b->remapflag[0] = newimage - oldimage;
+ oldimage = (original >> IMGBITS) & IMGMASK;
+ newimage = (b->image >> IMGBITS) & IMGMASK;
+ b->remapflag[1] = newimage - oldimage;
+ oldimage = original >> IMG2BITS;
+ newimage = b->image >> IMG2BITS;
+ b->remapflag[2] = newimage - oldimage;
+ b->remapflag[3] = 1;
+ }
+ }
+
+ // acquire ghost bodies via forward comm
+ // also gets new remapflags needed for adjusting atom image flags
+ // reset atom2body for owned atoms
+ // forward comm sets atom2body for ghost atoms
+
+ nghost_body = 0;
+ commflag = FULL_BODY;
+ comm->forward_comm_variable_fix(this);
+ reset_atom2body();
+ //check(4);
+
+ // adjust image flags of any atom in a rigid body whose xcm was remapped
+
+ tagint *image = atom->image;
+ int nlocal = atom->nlocal;
+
+ int ibody;
+ tagint idim,otherdims;
+
+ for (int i = 0; i < nlocal; i++) {
+ if (atom2body[i] < 0) continue;
+ Body *b = &body[atom2body[i]];
+
+ if (b->remapflag[3] == 0) continue;
+
+ if (b->remapflag[0]) {
+ idim = image[i] & IMGMASK;
+ otherdims = image[i] ^ idim;
+ idim -= b->remapflag[0];
+ idim &= IMGMASK;
+ image[i] = otherdims | idim;
+ }
+ if (b->remapflag[1]) {
+ idim = (image[i] >> IMGBITS) & IMGMASK;
+ otherdims = image[i] ^ (idim << IMGBITS);
+ idim -= b->remapflag[1];
+ idim &= IMGMASK;
+ image[i] = otherdims | (idim << IMGBITS);
+ }
+ if (b->remapflag[2]) {
+ idim = image[i] >> IMG2BITS;
+ otherdims = image[i] ^ (idim << IMG2BITS);
+ idim -= b->remapflag[2];
+ idim &= IMGMASK;
+ image[i] = otherdims | (idim << IMG2BITS);
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ count # of DOF removed by rigid bodies for atoms in igroup
+ return total count of DOF
+------------------------------------------------------------------------- */
+
+int FixRigidSmall::dof(int tgroup)
+{
+ int i,j;
+
+ int tgroupbit = group->bitmask[tgroup];
+
+ // if firstflag, just return 0
+ // is being called by first run init()
+ // local comm stencil is not setup and rigid body inertia is not calculated
+ // will be called again by temperature compute setup()
+
+ if (firstflag) return 0;
+
+ // counts = 3 values per rigid body I own
+ // 0 = # of point particles in rigid body and in temperature group
+ // 1 = # of finite-size particles in rigid body and in temperature group
+ // 2 = # of particles in rigid body, disregarding temperature group
+
+ memory->create(counts,nlocal_body+nghost_body,3,"rigid/small:counts");
+ for (int i = 0; i < nlocal_body+nghost_body; i++)
+ counts[i][0] = counts[i][1] = counts[i][2] = 0;
+
+ // tally counts from my owned atoms
+ // 0 = # of point particles in rigid body and in temperature group
+ // 1 = # of finite-size particles in rigid body and in temperature group
+ // 2 = # of particles in rigid body, disregarding temperature group
+
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+
+ for (i = 0; i < nlocal; i++) {
+ if (atom2body[i] < 0) continue;
+ j = atom2body[i];
+ counts[j][2]++;
+ if (mask[i] & tgroupbit) {
+ if (extended && eflags[i]) counts[j][1]++;
+ else counts[j][0]++;
+ }
+ }
+
+ commflag = DOF;
+ comm->reverse_comm_variable_fix(this);
+
+ // nall = count0 = # of point particles in each rigid body
+ // mall = count1 = # of finite-size particles in each rigid body
+ // warn if nall+mall != nrigid for any body included in temperature group
+
+ int flag = 0;
+ for (int ibody = 0; ibody < nlocal_body; ibody++) {
+ if (counts[ibody][0]+counts[ibody][1] > 0 &&
+ counts[ibody][0]+counts[ibody][1] != counts[ibody][2]) flag = 1;
+ }
+ int flagall;
+ MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_MAX,world);
+ if (flagall && me == 0)
+ error->warning(FLERR,"Computing temperature of portions of rigid bodies");
+
+ // remove appropriate DOFs for each rigid body wholly in temperature group
+ // N = # of point particles in body
+ // M = # of finite-size particles in body
+ // 3d body has 3N + 6M dof to start with
+ // 2d body has 2N + 3M dof to start with
+ // 3d point-particle body with all non-zero I should have 6 dof, remove 3N-6
+ // 3d point-particle body (linear) with a 0 I should have 5 dof, remove 3N-5
+ // 2d point-particle body should have 3 dof, remove 2N-3
+ // 3d body with any finite-size M should have 6 dof, remove (3N+6M) - 6
+ // 2d body with any finite-size M should have 3 dof, remove (2N+3M) - 3
+
+ double *inertia;
+
+ int n = 0;
+ if (domain->dimension == 3) {
+ for (int ibody = 0; ibody < nlocal_body; ibody++) {
+ if (counts[ibody][0]+counts[ibody][1] == counts[ibody][2]) {
+ n += 3*counts[ibody][0] + 6*counts[ibody][1] - 6;
+ inertia = body[ibody].inertia;
+ if (inertia[0] == 0.0 || inertia[1] == 0.0 || inertia[2] == 0.0) n++;
+ }
+ }
+ } else if (domain->dimension == 2) {
+ for (int ibody = 0; ibody < nlocal_body; ibody++)
+ if (counts[ibody][0]+counts[ibody][1] == counts[ibody][2])
+ n += 2*counts[ibody][0] + 3*counts[ibody][1] - 3;
+ }
+
+ memory->destroy(counts);
+
+ int nall;
+ MPI_Allreduce(&n,&nall,1,MPI_INT,MPI_SUM,world);
+ return nall;
+}
+
+/* ----------------------------------------------------------------------
+ adjust xcm of each rigid body due to box deformation
+ called by various fixes that change box size/shape
+ flag = 0/1 means map from box to lamda coords or vice versa
+------------------------------------------------------------------------- */
+
+void FixRigidSmall::deform(int flag)
+{
+ if (flag == 0)
+ for (int ibody = 0; ibody < nlocal_body; ibody++)
+ domain->x2lamda(body[ibody].xcm,body[ibody].xcm);
+ else
+ for (int ibody = 0; ibody < nlocal_body; ibody++)
+ domain->lamda2x(body[ibody].xcm,body[ibody].xcm);
+}
+
+/* ----------------------------------------------------------------------
+ set space-frame coords and velocity of each atom in each rigid body
+ set orientation and rotation of extended particles
+ x = Q displace + Xcm, mapped back to periodic box
+ v = Vcm + (W cross (x - Xcm))
+------------------------------------------------------------------------- */
+
+void FixRigidSmall::set_xv()
+{
+ int ibody,itype;
+ int xbox,ybox,zbox;
+ double x0,x1,x2,v0,v1,v2,fc0,fc1,fc2,massone;
+ double ione[3],exone[3],eyone[3],ezone[3],vr[6],p[3][3];
+
+ double xprd = domain->xprd;
+ double yprd = domain->yprd;
+ double zprd = domain->zprd;
+ double xy = domain->xy;
+ double xz = domain->xz;
+ double yz = domain->yz;
+
+ tagint *image = atom->image;
+ double **x = atom->x;
+ double **v = atom->v;
+ double **f = atom->f;
+ double *rmass = atom->rmass;
+ double *mass = atom->mass;
+ int *type = atom->type;
+ int nlocal = atom->nlocal;
+
+ // set x and v of each atom
+
+ for (int i = 0; i < nlocal; i++) {
+ if (atom2body[i] < 0) continue;
+ Body *b = &body[atom2body[i]];
+
+ xbox = (image[i] & IMGMASK) - IMGMAX;
+ ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
+ zbox = (image[i] >> IMG2BITS) - IMGMAX;
+
+ // save old positions and velocities for virial
+
+ if (evflag) {
+ if (triclinic == 0) {
+ x0 = x[i][0] + xbox*xprd;
+ x1 = x[i][1] + ybox*yprd;
+ x2 = x[i][2] + zbox*zprd;
+ } else {
+ x0 = x[i][0] + xbox*xprd + ybox*xy + zbox*xz;
+ x1 = x[i][1] + ybox*yprd + zbox*yz;
+ x2 = x[i][2] + zbox*zprd;
+ }
+ v0 = v[i][0];
+ v1 = v[i][1];
+ v2 = v[i][2];
+ }
+
+ // x = displacement from center-of-mass, based on body orientation
+ // v = vcm + omega around center-of-mass
+
+ MathExtra::matvec(b->ex_space,b->ey_space,b->ez_space,displace[i],x[i]);
+
+ v[i][0] = b->omega[1]*x[i][2] - b->omega[2]*x[i][1] + b->vcm[0];
+ v[i][1] = b->omega[2]*x[i][0] - b->omega[0]*x[i][2] + b->vcm[1];
+ v[i][2] = b->omega[0]*x[i][1] - b->omega[1]*x[i][0] + b->vcm[2];
+
+ // add center of mass to displacement
+ // map back into periodic box via xbox,ybox,zbox
+ // for triclinic, add in box tilt factors as well
+
+ if (triclinic == 0) {
+ x[i][0] += b->xcm[0] - xbox*xprd;
+ x[i][1] += b->xcm[1] - ybox*yprd;
+ x[i][2] += b->xcm[2] - zbox*zprd;
+ } else {
+ x[i][0] += b->xcm[0] - xbox*xprd - ybox*xy - zbox*xz;
+ x[i][1] += b->xcm[1] - ybox*yprd - zbox*yz;
+ x[i][2] += b->xcm[2] - zbox*zprd;
+ }
+
+ // virial = unwrapped coords dotted into body constraint force
+ // body constraint force = implied force due to v change minus f external
+ // assume f does not include forces internal to body
+ // 1/2 factor b/c final_integrate contributes other half
+ // assume per-atom contribution is due to constraint force on that atom
+
+ if (evflag) {
+ if (rmass) massone = rmass[i];
+ else massone = mass[type[i]];
+ fc0 = massone*(v[i][0] - v0)/dtf - f[i][0];
+ fc1 = massone*(v[i][1] - v1)/dtf - f[i][1];
+ fc2 = massone*(v[i][2] - v2)/dtf - f[i][2];
+
+ vr[0] = 0.5*x0*fc0;
+ vr[1] = 0.5*x1*fc1;
+ vr[2] = 0.5*x2*fc2;
+ vr[3] = 0.5*x0*fc1;
+ vr[4] = 0.5*x0*fc2;
+ vr[5] = 0.5*x1*fc2;
+
+ v_tally(1,&i,1.0,vr);
+ }
+ }
+
+ // set orientation, omega, angmom of each extended particle
+
+ if (extended) {
+ double theta_body,theta;
+ double *shape,*quatatom,*inertiaatom;
+
+ AtomVecEllipsoid::Bonus *ebonus;
+ if (avec_ellipsoid) ebonus = avec_ellipsoid->bonus;
+ AtomVecLine::Bonus *lbonus;
+ if (avec_line) lbonus = avec_line->bonus;
+ AtomVecTri::Bonus *tbonus;
+ if (avec_tri) tbonus = avec_tri->bonus;
+ double **omega = atom->omega;
+ double **angmom = atom->angmom;
+ double **mu = atom->mu;
+ int *ellipsoid = atom->ellipsoid;
+ int *line = atom->line;
+ int *tri = atom->tri;
+
+ for (int i = 0; i < nlocal; i++) {
+ if (atom2body[i] < 0) continue;
+ Body *b = &body[atom2body[i]];
+
+ if (eflags[i] & SPHERE) {
+ omega[i][0] = b->omega[0];
+ omega[i][1] = b->omega[1];
+ omega[i][2] = b->omega[2];
+ } else if (eflags[i] & ELLIPSOID) {
+ shape = ebonus[ellipsoid[i]].shape;
+ quatatom = ebonus[ellipsoid[i]].quat;
+ MathExtra::quatquat(b->quat,orient[i],quatatom);
+ MathExtra::qnormalize(quatatom);
+ ione[0] = EINERTIA*rmass[i] * (shape[1]*shape[1] + shape[2]*shape[2]);
+ ione[1] = EINERTIA*rmass[i] * (shape[0]*shape[0] + shape[2]*shape[2]);
+ ione[2] = EINERTIA*rmass[i] * (shape[0]*shape[0] + shape[1]*shape[1]);
+ MathExtra::q_to_exyz(quatatom,exone,eyone,ezone);
+ MathExtra::omega_to_angmom(b->omega,exone,eyone,ezone,ione,angmom[i]);
+ } else if (eflags[i] & LINE) {
+ if (b->quat[3] >= 0.0) theta_body = 2.0*acos(b->quat[0]);
+ else theta_body = -2.0*acos(b->quat[0]);
+ theta = orient[i][0] + theta_body;
+ while (theta <= MINUSPI) theta += TWOPI;
+ while (theta > MY_PI) theta -= TWOPI;
+ lbonus[line[i]].theta = theta;
+ omega[i][0] = b->omega[0];
+ omega[i][1] = b->omega[1];
+ omega[i][2] = b->omega[2];
+ } else if (eflags[i] & TRIANGLE) {
+ inertiaatom = tbonus[tri[i]].inertia;
+ quatatom = tbonus[tri[i]].quat;
+ MathExtra::quatquat(b->quat,orient[i],quatatom);
+ MathExtra::qnormalize(quatatom);
+ MathExtra::q_to_exyz(quatatom,exone,eyone,ezone);
+ MathExtra::omega_to_angmom(b->omega,exone,eyone,ezone,
+ inertiaatom,angmom[i]);
+ }
+ if (eflags[i] & DIPOLE) {
+ MathExtra::quat_to_mat(b->quat,p);
+ MathExtra::matvec(p,dorient[i],mu[i]);
+ MathExtra::snormalize3(mu[i][3],mu[i],mu[i]);
+ }
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ set space-frame velocity of each atom in a rigid body
+ set omega and angmom of extended particles
+ v = Vcm + (W cross (x - Xcm))
+------------------------------------------------------------------------- */
+
+void FixRigidSmall::set_v()
+{
+ int ibody,itype;
+ int xbox,ybox,zbox;
+ double dx,dy,dz;
+ double x0,x1,x2,v0,v1,v2,fc0,fc1,fc2,massone;
+ double ione[3],exone[3],eyone[3],ezone[3],delta[3],vr[6];
+
+ double xprd = domain->xprd;
+ double yprd = domain->yprd;
+ double zprd = domain->zprd;
+ double xy = domain->xy;
+ double xz = domain->xz;
+ double yz = domain->yz;
+
+ double **x = atom->x;
+ double **v = atom->v;
+ double **f = atom->f;
+ double *rmass = atom->rmass;
+ double *mass = atom->mass;
+ int *type = atom->type;
+ tagint *image = atom->image;
+ int nlocal = atom->nlocal;
+
+ // set v of each atom
+
+ for (int i = 0; i < nlocal; i++) {
+ if (atom2body[i] < 0) continue;
+ Body *b = &body[atom2body[i]];
+
+ MathExtra::matvec(b->ex_space,b->ey_space,b->ez_space,displace[i],delta);
+
+ // save old velocities for virial
+
+ if (evflag) {
+ v0 = v[i][0];
+ v1 = v[i][1];
+ v2 = v[i][2];
+ }
+
+ v[i][0] = b->omega[1]*delta[2] - b->omega[2]*delta[1] + b->vcm[0];
+ v[i][1] = b->omega[2]*delta[0] - b->omega[0]*delta[2] + b->vcm[1];
+ v[i][2] = b->omega[0]*delta[1] - b->omega[1]*delta[0] + b->vcm[2];
+
+ // virial = unwrapped coords dotted into body constraint force
+ // body constraint force = implied force due to v change minus f external
+ // assume f does not include forces internal to body
+ // 1/2 factor b/c initial_integrate contributes other half
+ // assume per-atom contribution is due to constraint force on that atom
+
+ if (evflag) {
+ if (rmass) massone = rmass[i];
+ else massone = mass[type[i]];
+ fc0 = massone*(v[i][0] - v0)/dtf - f[i][0];
+ fc1 = massone*(v[i][1] - v1)/dtf - f[i][1];
+ fc2 = massone*(v[i][2] - v2)/dtf - f[i][2];
+
+ xbox = (image[i] & IMGMASK) - IMGMAX;
+ ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
+ zbox = (image[i] >> IMG2BITS) - IMGMAX;
+
+ if (triclinic == 0) {
+ x0 = x[i][0] + xbox*xprd;
+ x1 = x[i][1] + ybox*yprd;
+ x2 = x[i][2] + zbox*zprd;
+ } else {
+ x0 = x[i][0] + xbox*xprd + ybox*xy + zbox*xz;
+ x1 = x[i][1] + ybox*yprd + zbox*yz;
+ x2 = x[i][2] + zbox*zprd;
+ }
+
+ vr[0] = 0.5*x0*fc0;
+ vr[1] = 0.5*x1*fc1;
+ vr[2] = 0.5*x2*fc2;
+ vr[3] = 0.5*x0*fc1;
+ vr[4] = 0.5*x0*fc2;
+ vr[5] = 0.5*x1*fc2;
+
+ v_tally(1,&i,1.0,vr);
+ }
+ }
+
+ // set omega, angmom of each extended particle
+
+ if (extended) {
+ double *shape,*quatatom,*inertiaatom;
+
+ AtomVecEllipsoid::Bonus *ebonus;
+ if (avec_ellipsoid) ebonus = avec_ellipsoid->bonus;
+ AtomVecTri::Bonus *tbonus;
+ if (avec_tri) tbonus = avec_tri->bonus;
+ double **omega = atom->omega;
+ double **angmom = atom->angmom;
+ int *ellipsoid = atom->ellipsoid;
+ int *tri = atom->tri;
+
+ for (int i = 0; i < nlocal; i++) {
+ if (atom2body[i] < 0) continue;
+ Body *b = &body[atom2body[i]];
+
+ if (eflags[i] & SPHERE) {
+ omega[i][0] = b->omega[0];
+ omega[i][1] = b->omega[1];
+ omega[i][2] = b->omega[2];
+ } else if (eflags[i] & ELLIPSOID) {
+ shape = ebonus[ellipsoid[i]].shape;
+ quatatom = ebonus[ellipsoid[i]].quat;
+ ione[0] = EINERTIA*rmass[i] * (shape[1]*shape[1] + shape[2]*shape[2]);
+ ione[1] = EINERTIA*rmass[i] * (shape[0]*shape[0] + shape[2]*shape[2]);
+ ione[2] = EINERTIA*rmass[i] * (shape[0]*shape[0] + shape[1]*shape[1]);
+ MathExtra::q_to_exyz(quatatom,exone,eyone,ezone);
+ MathExtra::omega_to_angmom(b->omega,exone,eyone,ezone,ione,
+ angmom[i]);
+ } else if (eflags[i] & LINE) {
+ omega[i][0] = b->omega[0];
+ omega[i][1] = b->omega[1];
+ omega[i][2] = b->omega[2];
+ } else if (eflags[i] & TRIANGLE) {
+ inertiaatom = tbonus[tri[i]].inertia;
+ quatatom = tbonus[tri[i]].quat;
+ MathExtra::q_to_exyz(quatatom,exone,eyone,ezone);
+ MathExtra::omega_to_angmom(b->omega,exone,eyone,ezone,
+ inertiaatom,angmom[i]);
+ }
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ one-time identification of which atoms are in which rigid bodies
+ set bodytag for all owned atoms
+------------------------------------------------------------------------- */
+
+void FixRigidSmall::create_bodies()
+{
+ int i,m,n;
+ double unwrap[3];
+
+ // error check on image flags of atoms in rigid bodies
+
+ int *image = atom->image;
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+
+ int *periodicity = domain->periodicity;
+ int xbox,ybox,zbox;
+
+ int flag = 0;
+ for (i = 0; i < nlocal; i++) {
+ if (!(mask[i] & groupbit)) continue;
+ xbox = (image[i] & IMGMASK) - IMGMAX;
+ ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
+ zbox = (image[i] >> IMG2BITS) - IMGMAX;
+ if ((xbox && !periodicity[0]) || (ybox && !periodicity[1]) ||
+ (zbox && !periodicity[2])) flag = 1;
+ }
+
+ int flagall;
+ MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
+ if (flagall) error->all(FLERR,"Fix rigid/small atom has non-zero image flag "
+ "in a non-periodic dimension");
+
+ // allocate buffer for passing messages around ring of procs
+ // percount = max number of values to put in buffer for each of ncount
+
+ int ncount = 0;
+ for (i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit) ncount++;
+
+ int percount = 5;
+ double *buf;
+ memory->create(buf,ncount*percount,"rigid/small:buf");
+
+ // create map hash for storing unique molecule IDs of my atoms
+ // key = molecule ID
+ // value = index into per-body data structure
+ // n = # of entries in hash
+
+ hash = new std::map();
+ hash->clear();
+
+ // setup hash
+ // key = body ID
+ // value = index into N-length data structure
+ // n = count of unique bodies my atoms are part of
+
+ int *molecule = atom->molecule;
+
+ n = 0;
+ for (i = 0; i < nlocal; i++) {
+ if (!(mask[i] & groupbit)) continue;
+ if (hash->find(molecule[i]) == hash->end()) {
+ hash->insert(std::pair (molecule[i],n));
+ n++;
+ }
+ }
+
+ // bbox = bounding box of each rigid body my atoms are part of
+
+ memory->create(bbox,n,6,"rigid/small:bbox");
+
+ for (i = 0; i < n; i++) {
+ bbox[i][0] = bbox[i][2] = bbox[i][4] = BIG;
+ bbox[i][1] = bbox[i][3] = bbox[i][5] = -BIG;
+ }
+
+ // pack my atoms into buffer as molecule ID, unwrapped coords
+
+ double **x = atom->x;
+
+ m = 0;
+ for (i = 0; i < nlocal; i++) {
+ if (!(mask[i] & groupbit)) continue;
+ domain->unmap(x[i],image[i],unwrap);
+ buf[m++] = molecule[i];
+ buf[m++] = unwrap[0];
+ buf[m++] = unwrap[1];
+ buf[m++] = unwrap[2];
+ }
+
+ // pass buffer around ring of procs
+ // func = update bbox with atom coords from every proc
+ // when done, have full bbox for every rigid body my atoms are part of
+
+ frsptr = this;
+ comm->ring(m,sizeof(double),buf,1,bbox_update,NULL);
+
+ // ctr = center pt of each rigid body my atoms are part of
+
+ memory->create(ctr,n,6,"rigid/small:bbox");
+
+ for (i = 0; i < n; i++) {
+ ctr[i][0] = 0.5 * (bbox[i][0] + bbox[i][1]);
+ ctr[i][1] = 0.5 * (bbox[i][2] + bbox[i][3]);
+ ctr[i][2] = 0.5 * (bbox[i][4] + bbox[i][5]);
+ }
+
+ // idclose = ID of atom in body closest to center pt (smaller ID if tied)
+ // rsq = distance squared from idclose to center pt
+
+ memory->create(idclose,n,"rigid/small:idclose");
+ memory->create(rsqclose,n,"rigid/small:rsqclose");
+
+ for (i = 0; i < n; i++) rsqclose[i] = BIG;
+
+ // pack my atoms into buffer as molecule ID, atom ID, unwrapped coords
+
+ int *tag = atom->tag;
+
+ m = 0;
+ for (i = 0; i < nlocal; i++) {
+ if (!(mask[i] & groupbit)) continue;
+ domain->unmap(x[i],image[i],unwrap);
+ buf[m++] = molecule[i];
+ buf[m++] = tag[i];
+ buf[m++] = unwrap[0];
+ buf[m++] = unwrap[1];
+ buf[m++] = unwrap[2];
+ }
+
+ // pass buffer around ring of procs
+ // func = idclose,rsqclose with atom IDs from every proc
+ // when done, have idclose for every rigid body my atoms are part of
+
+ frsptr = this;
+ comm->ring(m,sizeof(double),buf,2,nearest_update,NULL);
+
+ // set bodytag of all owned atoms, based on idclose
+
+ for (i = 0; i < nlocal; i++) {
+ bodytag[i] = 0;
+ if (!(mask[i] & groupbit)) continue;
+ m = hash->find(molecule[i])->second;
+ bodytag[i] = idclose[m];
+ }
+
+ // clean up
+
+ delete hash;
+ memory->destroy(buf);
+ memory->destroy(bbox);
+ memory->destroy(ctr);
+ memory->destroy(idclose);
+ memory->destroy(rsqclose);
+}
+
+/* ----------------------------------------------------------------------
+ process rigid body atoms from another proc
+ update bounding box for rigid bodies my atoms are part of
+------------------------------------------------------------------------- */
+
+void FixRigidSmall::bbox_update(int n, char *cbuf)
+{
+ std::map *hash = frsptr->hash;
+ double **bbox = frsptr->bbox;
+
+ double *buf = (double *) cbuf;
+ int ndatums = n/4;
+
+ int j,imol;
+ double *x;
+
+ int m = 0;
+ for (int i = 0; i < ndatums; i++, m += 4) {
+ imol = static_cast (buf[m]);
+ if (hash->find(imol) != hash->end()) {
+ j = hash->find(imol)->second;
+ x = &buf[m+1];
+ bbox[j][0] = MIN(bbox[j][0],x[0]);
+ bbox[j][1] = MAX(bbox[j][1],x[0]);
+ bbox[j][2] = MIN(bbox[j][2],x[1]);
+ bbox[j][3] = MAX(bbox[j][3],x[1]);
+ bbox[j][4] = MIN(bbox[j][4],x[2]);
+ bbox[j][5] = MAX(bbox[j][5],x[2]);
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ process rigid body atoms from another proc
+ update nearest atom to body center for rigid bodies my atoms are part of
+------------------------------------------------------------------------- */
+
+void FixRigidSmall::nearest_update(int n, char *cbuf)
+{
+ std::map *hash = frsptr->hash;
+ double **ctr = frsptr->ctr;
+ int *idclose = frsptr->idclose;
+ double *rsqclose = frsptr->rsqclose;
+
+ double *buf = (double *) cbuf;
+ int ndatums = n/5;
+
+ int j,imol,tag;
+ double delx,dely,delz,rsq;
+ double *x;
+
+ int m = 0;
+ for (int i = 0; i < ndatums; i++, m += 5) {
+ imol = static_cast (buf[m]);
+ if (hash->find(imol) != hash->end()) {
+ j = hash->find(imol)->second;
+ tag = static_cast (buf[m+1]);
+ x = &buf[m+2];
+ delx = x[0] - ctr[j][0];
+ dely = x[1] - ctr[j][1];
+ delz = x[2] - ctr[j][2];
+ rsq = delx*delx + dely*dely + delz*delz;
+ if (rsq <= rsqclose[j]) {
+ if (rsq == rsqclose[j] && tag > idclose[j]) continue;
+ idclose[j] = tag;
+ rsqclose[j] = rsq;
+ }
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ one-time initialization of rigid body attributes
+ extended flags, mass, center-of-mass
+ Cartesian and diagonalized inertia tensor
+------------------------------------------------------------------------- */
+
+void FixRigidSmall::setup_bodies()
+{
+ int i,itype,ibody;
+
+ // extended = 1 if any particle in a rigid body is finite size
+ // or has a dipole moment
+
+ extended = orientflag = dorientflag = 0;
+
+ AtomVecEllipsoid::Bonus *ebonus;
+ if (avec_ellipsoid) ebonus = avec_ellipsoid->bonus;
+ AtomVecLine::Bonus *lbonus;
+ if (avec_line) lbonus = avec_line->bonus;
+ AtomVecTri::Bonus *tbonus;
+ if (avec_tri) tbonus = avec_tri->bonus;
+ double **mu = atom->mu;
+ double *radius = atom->radius;
+ double *rmass = atom->rmass;
+ double *mass = atom->mass;
+ int *ellipsoid = atom->ellipsoid;
+ int *line = atom->line;
+ int *tri = atom->tri;
+ int *type = atom->type;
+ int nlocal = atom->nlocal;
+
+ if (atom->radius_flag || atom->ellipsoid_flag || atom->line_flag ||
+ atom->tri_flag || atom->mu_flag) {
+ int flag = 0;
+ for (i = 0; i < nlocal; i++) {
+ if (atom2body[i] < 0) continue;
+ if (radius && radius[i] > 0.0) flag = 1;
+ if (ellipsoid && ellipsoid[i] >= 0) flag = 1;
+ if (line && line[i] >= 0) flag = 1;
+ if (tri && tri[i] >= 0) flag = 1;
+ if (mu && mu[i][3] > 0.0) flag = 1;
+ }
+
+ MPI_Allreduce(&flag,&extended,1,MPI_INT,MPI_MAX,world);
+ }
+
+ // grow extended arrays and set extended flags for each particle
+ // orientflag = 4 if any particle stores ellipsoid or tri orientation
+ // orientflag = 1 if any particle stores line orientation
+ // dorientflag = 1 if any particle stores dipole orientation
+
+ if (extended) {
+ if (atom->ellipsoid_flag) orientflag = 4;
+ if (atom->line_flag) orientflag = 1;
+ if (atom->tri_flag) orientflag = 4;
+ if (atom->mu_flag) dorientflag = 1;
+ grow_arrays(atom->nmax);
+
+ for (i = 0; i < nlocal; i++) {
+ eflags[i] = 0;
+ if (atom2body[i] < 0) continue;
+
+ // set to POINT or SPHERE or ELLIPSOID or LINE
+
+ if (radius && radius[i] > 0.0) {
+ eflags[i] |= SPHERE;
+ eflags[i] |= OMEGA;
+ eflags[i] |= TORQUE;
+ } else if (ellipsoid && ellipsoid[i] >= 0) {
+ eflags[i] |= ELLIPSOID;
+ eflags[i] |= ANGMOM;
+ eflags[i] |= TORQUE;
+ } else if (line && line[i] >= 0) {
+ eflags[i] |= LINE;
+ eflags[i] |= OMEGA;
+ eflags[i] |= TORQUE;
+ } else if (tri && tri[i] >= 0) {
+ eflags[i] |= TRIANGLE;
+ eflags[i] |= ANGMOM;
+ eflags[i] |= TORQUE;
+ } else eflags[i] |= POINT;
+
+ // set DIPOLE if atom->mu and mu[3] > 0.0
+
+ if (atom->mu_flag && mu[i][3] > 0.0)
+ eflags[i] |= DIPOLE;
+ }
+ }
+
+ // acquire ghost bodies via forward comm
+ // set atom2body for ghost atoms via forward comm
+ // set atom2body for other owned atoms via reset_atom2body()
+
+ nghost_body = 0;
+ commflag = FULL_BODY;
+ comm->forward_comm_variable_fix(this);
+ reset_atom2body();
+
+ // compute mass & center-of-mass of each rigid body
+ // error if image flag is not 0 in a non-periodic dim
+
+ double **x = atom->x;
+ tagint *image = atom->image;
+
+ double *xcm;
+
+ for (ibody = 0; ibody < nlocal_body+nghost_body; ibody++) {
+ xcm = body[ibody].xcm;
+ xcm[0] = xcm[1] = xcm[2] = 0.0;
+ body[ibody].mass = 0.0;
+ }
+
+ double unwrap[3];
+ double massone;
+
+ for (i = 0; i < nlocal; i++) {
+ if (atom2body[i] < 0) continue;
+ Body *b = &body[atom2body[i]];
+
+ if (rmass) massone = rmass[i];
+ else massone = mass[type[i]];
+
+ domain->unmap(x[i],image[i],unwrap);
+ xcm = b->xcm;
+ xcm[0] += unwrap[0] * massone;
+ xcm[1] += unwrap[1] * massone;
+ xcm[2] += unwrap[2] * massone;
+ b->mass += massone;
+ }
+
+ // reverse communicate xcm, mass of all bodies
+
+ commflag = XCM_MASS;
+ comm->reverse_comm_variable_fix(this);
+
+ for (ibody = 0; ibody < nlocal_body; ibody++) {
+ xcm = body[ibody].xcm;
+ xcm[0] /= body[ibody].mass;
+ xcm[1] /= body[ibody].mass;
+ xcm[2] /= body[ibody].mass;
+ }
+
+ // set image flags for each rigid body to default values
+ // then remap the xcm of each body back into simulation box if needed
+
+ for (ibody = 0; ibody < nlocal_body; ibody++)
+ body[ibody].image = ((tagint) IMGMAX << IMG2BITS) |
+ ((tagint) IMGMAX << IMGBITS) | IMGMAX;
+
+ pre_neighbor();
+
+ // compute 6 moments of inertia of each body in Cartesian reference frame
+ // dx,dy,dz = coords relative to center-of-mass
+ // symmetric 3x3 inertia tensor stored in Voigt notation as 6-vector
+
+ double dx,dy,dz,rad;
+ double *inertia;
+
+ for (ibody = 0; ibody < nlocal_body+nghost_body; ibody++)
+ for (i = 0; i < 6; i++) body[ibody].itensor[i] = 0.0;
+
+ for (i = 0; i < nlocal; i++) {
+ if (atom2body[i] < 0) continue;
+ Body *b = &body[atom2body[i]];
+
+ domain->unmap(x[i],image[i],unwrap);
+ xcm = b->xcm;
+ dx = unwrap[0] - xcm[0];
+ dy = unwrap[1] - xcm[1];
+ dz = unwrap[2] - xcm[2];
+
+ if (rmass) massone = rmass[i];
+ else massone = mass[type[i]];
+
+ inertia = b->itensor;
+ inertia[0] += massone * (dy*dy + dz*dz);
+ inertia[1] += massone * (dx*dx + dz*dz);
+ inertia[2] += massone * (dx*dx + dy*dy);
+ inertia[3] -= massone * dy*dz;
+ inertia[4] -= massone * dx*dz;
+ inertia[5] -= massone * dx*dy;
+ }
+
+ // extended particles may contribute extra terms to moments of inertia
+
+ if (extended) {
+ double ivec[6];
+ double *shape,*quatatom,*inertiaatom;
+ double length,theta;
+
+ for (i = 0; i < nlocal; i++) {
+ if (atom2body[i] < 0) continue;
+ inertia = body[atom2body[i]].itensor;
+
+ if (rmass) massone = rmass[i];
+ else massone = mass[type[i]];
+
+ if (eflags[i] & SPHERE) {
+ inertia[0] += SINERTIA*massone * radius[i]*radius[i];
+ inertia[1] += SINERTIA*massone * radius[i]*radius[i];
+ inertia[2] += SINERTIA*massone * radius[i]*radius[i];
+ } else if (eflags[i] & ELLIPSOID) {
+ shape = ebonus[ellipsoid[i]].shape;
+ quatatom = ebonus[ellipsoid[i]].quat;
+ MathExtra::inertia_ellipsoid(shape,quatatom,massone,ivec);
+ inertia[0] += ivec[0];
+ inertia[1] += ivec[1];
+ inertia[2] += ivec[2];
+ inertia[3] += ivec[3];
+ inertia[4] += ivec[4];
+ inertia[5] += ivec[5];
+ } else if (eflags[i] & LINE) {
+ length = lbonus[line[i]].length;
+ theta = lbonus[line[i]].theta;
+ MathExtra::inertia_line(length,theta,massone,ivec);
+ inertia[0] += ivec[0];
+ inertia[1] += ivec[1];
+ inertia[2] += ivec[2];
+ inertia[3] += ivec[3];
+ inertia[4] += ivec[4];
+ inertia[5] += ivec[5];
+ } else if (eflags[i] & TRIANGLE) {
+ inertiaatom = tbonus[tri[i]].inertia;
+ quatatom = tbonus[tri[i]].quat;
+ MathExtra::inertia_triangle(inertiaatom,quatatom,massone,ivec);
+ inertia[0] += ivec[0];
+ inertia[1] += ivec[1];
+ inertia[2] += ivec[2];
+ inertia[3] += ivec[3];
+ inertia[4] += ivec[4];
+ inertia[5] += ivec[5];
+ }
+ }
+ }
+
+ for (ibody = 0; ibody < nlocal_body; ibody++) {
+ inertia = body[atom2body[i]].itensor;
+ }
+
+ // reverse communicate inertia tensor of all bodies
+
+ commflag = ITENSOR;
+ comm->reverse_comm_variable_fix(this);
+
+ // diagonalize inertia tensor for each body via Jacobi rotations
+ // inertia = 3 eigenvalues = principal moments of inertia
+ // evectors and exzy_space = 3 evectors = principal axes of rigid body
+
+ int ierror;
+ double cross[3];
+ double tensor[3][3],evectors[3][3];
+ double *itensor,*ex,*ey,*ez;
+
+ for (ibody = 0; ibody < nlocal_body; ibody++) {
+ itensor = body[ibody].itensor;
+
+ tensor[0][0] = itensor[0];
+ tensor[1][1] = itensor[1];
+ tensor[2][2] = itensor[2];
+ tensor[1][2] = tensor[2][1] = itensor[3];
+ tensor[0][2] = tensor[2][0] = itensor[4];
+ tensor[0][1] = tensor[1][0] = itensor[5];
+
+ inertia = body[ibody].inertia;
+ ierror = MathExtra::jacobi(tensor,inertia,evectors);
+ if (ierror) error->all(FLERR,
+ "Insufficient Jacobi rotations for rigid body");
+
+ ex = body[ibody].ex_space;
+ ex[0] = evectors[0][0];
+ ex[1] = evectors[1][0];
+ ex[2] = evectors[2][0];
+ ey = body[ibody].ey_space;
+ ey[0] = evectors[0][1];
+ ey[1] = evectors[1][1];
+ ey[2] = evectors[2][1];
+ ez = body[ibody].ez_space;
+ ez[0] = evectors[0][2];
+ ez[1] = evectors[1][2];
+ ez[2] = evectors[2][2];
+
+ // if any principal moment < scaled EPSILON, set to 0.0
+
+ double max;
+ max = MAX(inertia[0],inertia[1]);
+ max = MAX(max,inertia[2]);
+
+ if (inertia[0] < EPSILON*max) inertia[0] = 0.0;
+ if (inertia[1] < EPSILON*max) inertia[1] = 0.0;
+ if (inertia[2] < EPSILON*max) inertia[2] = 0.0;
+
+ // enforce 3 evectors as a right-handed coordinate system
+ // flip 3rd vector if needed
+
+ MathExtra::cross3(ex,ey,cross);
+ if (MathExtra::dot3(cross,ez) < 0.0) MathExtra::negate3(ez);
+
+ // create initial quaternion
+
+ MathExtra::exyz_to_q(ex,ey,ez,body[ibody].quat);
+ }
+
+ // forward communicate updated info of all bodies
+
+ commflag = INITIAL;
+ comm->forward_comm_variable_fix(this);
+
+ // displace = initial atom coords in basis of principal axes
+ // set displace = 0.0 for atoms not in any rigid body
+ // for extended particles, set their orientation wrt to rigid body
+
+ double qc[4],delta[3];
+ double *quatatom;
+ double theta_body;
+
+ for (i = 0; i < nlocal; i++) {
+ if (atom2body[i] < 0) {
+ displace[i][0] = displace[i][1] = displace[i][2] = 0.0;
+ continue;
+ }
+
+ Body *b = &body[atom2body[i]];
+
+ domain->unmap(x[i],image[i],unwrap);
+ xcm = b->xcm;
+ delta[0] = unwrap[0] - xcm[0];
+ delta[1] = unwrap[1] - xcm[1];
+ delta[2] = unwrap[2] - xcm[2];
+ MathExtra::transpose_matvec(b->ex_space,b->ey_space,b->ez_space,
+ delta,displace[i]);
+
+ if (extended) {
+ if (eflags[i] & ELLIPSOID) {
+ quatatom = ebonus[ellipsoid[i]].quat;
+ MathExtra::qconjugate(b->quat,qc);
+ MathExtra::quatquat(qc,quatatom,orient[i]);
+ MathExtra::qnormalize(orient[i]);
+ } else if (eflags[i] & LINE) {
+ if (b->quat[3] >= 0.0) theta_body = 2.0*acos(b->quat[0]);
+ else theta_body = -2.0*acos(b->quat[0]);
+ orient[i][0] = lbonus[line[i]].theta - theta_body;
+ while (orient[i][0] <= MINUSPI) orient[i][0] += TWOPI;
+ while (orient[i][0] > MY_PI) orient[i][0] -= TWOPI;
+ if (orientflag == 4) orient[i][1] = orient[i][2] = orient[i][3] = 0.0;
+ } else if (eflags[i] & TRIANGLE) {
+ quatatom = tbonus[tri[i]].quat;
+ MathExtra::qconjugate(b->quat,qc);
+ MathExtra::quatquat(qc,quatatom,orient[i]);
+ MathExtra::qnormalize(orient[i]);
+ } else if (orientflag == 4) {
+ orient[i][0] = orient[i][1] = orient[i][2] = orient[i][3] = 0.0;
+ } else if (orientflag == 1)
+ orient[i][0] = 0.0;
+
+ if (eflags[i] & DIPOLE) {
+ MathExtra::transpose_matvec(b->ex_space,b->ey_space,b->ez_space,
+ mu[i],dorient[i]);
+ MathExtra::snormalize3(mu[i][3],dorient[i],dorient[i]);
+ } else if (dorientflag)
+ dorient[i][0] = dorient[i][1] = dorient[i][2] = 0.0;
+ }
+ }
+
+ // test for valid principal moments & axes
+ // recompute moments of inertia around new axes
+ // 3 diagonal moments should equal principal moments
+ // 3 off-diagonal moments should be 0.0
+ // extended particles may contribute extra terms to moments of inertia
+
+ for (ibody = 0; ibody < nlocal_body+nghost_body; ibody++)
+ for (i = 0; i < 6; i++) body[ibody].itensor[i] = 0.0;
+
+ for (i = 0; i < nlocal; i++) {
+ if (atom2body[i] < 0) continue;
+ inertia = body[atom2body[i]].itensor;
+
+ if (rmass) massone = rmass[i];
+ else massone = mass[type[i]];
+
+ inertia[0] += massone *
+ (displace[i][1]*displace[i][1] + displace[i][2]*displace[i][2]);
+ inertia[1] += massone *
+ (displace[i][0]*displace[i][0] + displace[i][2]*displace[i][2]);
+ inertia[2] += massone *
+ (displace[i][0]*displace[i][0] + displace[i][1]*displace[i][1]);
+ inertia[3] -= massone * displace[i][1]*displace[i][2];
+ inertia[4] -= massone * displace[i][0]*displace[i][2];
+ inertia[5] -= massone * displace[i][0]*displace[i][1];
+ }
+
+ if (extended) {
+ double ivec[6];
+ double *shape,*inertiaatom;
+ double length;
+
+ for (i = 0; i < nlocal; i++) {
+ if (atom2body[i] < 0) continue;
+ inertia = body[atom2body[i]].itensor;
+
+ if (rmass) massone = rmass[i];
+ else massone = mass[type[i]];
+
+ if (eflags[i] & SPHERE) {
+ inertia[0] += SINERTIA*massone * radius[i]*radius[i];
+ inertia[1] += SINERTIA*massone * radius[i]*radius[i];
+ inertia[2] += SINERTIA*massone * radius[i]*radius[i];
+ } else if (eflags[i] & ELLIPSOID) {
+ shape = ebonus[ellipsoid[i]].shape;
+ MathExtra::inertia_ellipsoid(shape,orient[i],massone,ivec);
+ inertia[0] += ivec[0];
+ inertia[1] += ivec[1];
+ inertia[2] += ivec[2];
+ inertia[3] += ivec[3];
+ inertia[4] += ivec[4];
+ inertia[5] += ivec[5];
+ } else if (eflags[i] & LINE) {
+ length = lbonus[line[i]].length;
+ MathExtra::inertia_line(length,orient[i][0],massone,ivec);
+ inertia[0] += ivec[0];
+ inertia[1] += ivec[1];
+ inertia[2] += ivec[2];
+ inertia[3] += ivec[3];
+ inertia[4] += ivec[4];
+ inertia[5] += ivec[5];
+ } else if (eflags[i] & TRIANGLE) {
+ inertiaatom = tbonus[tri[i]].inertia;
+ MathExtra::inertia_triangle(inertiaatom,orient[i],massone,ivec);
+ inertia[0] += ivec[0];
+ inertia[1] += ivec[1];
+ inertia[2] += ivec[2];
+ inertia[3] += ivec[3];
+ inertia[4] += ivec[4];
+ inertia[5] += ivec[5];
+ }
+ }
+ }
+
+ // reverse communicate inertia tensor of all bodies
+
+ commflag = ITENSOR;
+ comm->reverse_comm_variable_fix(this);
+
+ // error check that re-computed momemts of inertia match diagonalized ones
+
+ double *inew;
+
+ double norm;
+ for (ibody = 0; ibody < nlocal_body; ibody++) {
+ inertia = body[ibody].inertia;
+ itensor = body[ibody].itensor;
+
+ if (inertia[0] == 0.0) {
+ if (fabs(itensor[0]) > TOLERANCE)
+ error->all(FLERR,"Fix rigid: Bad principal moments");
+ } else {
+ if (fabs((itensor[0]-inertia[0])/inertia[0]) >
+ TOLERANCE) error->all(FLERR,"Fix rigid: Bad principal moments");
+ }
+ if (inertia[1] == 0.0) {
+ if (fabs(itensor[1]) > TOLERANCE)
+ error->all(FLERR,"Fix rigid: Bad principal moments");
+ } else {
+ if (fabs((itensor[1]-inertia[1])/inertia[1]) >
+ TOLERANCE) error->all(FLERR,"Fix rigid: Bad principal moments");
+ }
+ if (inertia[2] == 0.0) {
+ if (fabs(itensor[2]) > TOLERANCE)
+ error->all(FLERR,"Fix rigid: Bad principal moments");
+ } else {
+ if (fabs((itensor[2]-inertia[2])/inertia[2]) >
+ TOLERANCE) error->all(FLERR,"Fix rigid: Bad principal moments");
+ }
+ norm = (inertia[0] + inertia[1] + inertia[2]) / 3.0;
+ if (fabs(itensor[3]/norm) > TOLERANCE ||
+ fabs(itensor[4]/norm) > TOLERANCE ||
+ fabs(itensor[5]/norm) > TOLERANCE)
+ error->all(FLERR,"Fix rigid: Bad principal moments");
+ }
+}
+
+/* ----------------------------------------------------------------------
+ memory usage of local atom-based arrays
+------------------------------------------------------------------------- */
+
+double FixRigidSmall::memory_usage()
+{
+ int nmax = atom->nmax;
+ double bytes = 2 * nmax * sizeof(int);
+ bytes += nmax*3 * sizeof(double);
+ bytes += maxvatom*6 * sizeof(double); // vatom
+ if (extended) {
+ bytes += nmax * sizeof(int);
+ if (orientflag) bytes = nmax*orientflag * sizeof(double);
+ if (dorientflag) bytes = nmax*3 * sizeof(double);
+ }
+ bytes += nmax_body * sizeof(Body);
+ return bytes;
+}
+
+/* ----------------------------------------------------------------------
+ allocate local atom-based arrays
+------------------------------------------------------------------------- */
+
+void FixRigidSmall::grow_arrays(int nmax)
+{
+ memory->grow(bodyown,nmax,"rigid/small:bodyown");
+ memory->grow(bodytag,nmax,"rigid/small:bodytag");
+ memory->grow(atom2body,nmax,"rigid/small:atom2body");
+ memory->grow(displace,nmax,3,"rigid/small:displace");
+ if (extended) {
+ memory->grow(eflags,nmax,"rigid/small:eflags");
+ if (orientflag) memory->grow(orient,nmax,orientflag,"rigid/small:orient");
+ if (dorientflag) memory->grow(dorient,nmax,3,"rigid/small:dorient");
+ }
+}
+
+/* ----------------------------------------------------------------------
+ copy values within local atom-based arrays
+------------------------------------------------------------------------- */
+
+void FixRigidSmall::copy_arrays(int i, int j, int delflag)
+{
+ bodytag[j] = bodytag[i];
+ displace[j][0] = displace[i][0];
+ displace[j][1] = displace[i][1];
+ displace[j][2] = displace[i][2];
+
+ if (extended) {
+ eflags[j] = eflags[i];
+ for (int k = 0; k < orientflag; k++)
+ orient[j][k] = orient[i][k];
+ if (dorientflag) {
+ dorient[j][0] = dorient[i][0];
+ dorient[j][1] = dorient[i][1];
+ dorient[j][2] = dorient[i][2];
+ }
+ }
+
+ // if deleting atom J via delflag and J owns a body, then delete it
+
+ if (delflag && bodyown[j] >= 0) {
+ bodyown[body[nlocal_body-1].ilocal] = bodyown[j];
+ memcpy(&body[bodyown[j]],&body[nlocal_body-1],sizeof(Body));
+ nlocal_body--;
+ }
+
+ // if atom I owns a body, reset I's body.ilocal to loc J
+ // do NOT do this if self-copy (I=J) since I's body is already deleted
+
+ if (bodyown[i] >= 0 && i != j) body[bodyown[i]].ilocal = j;
+ bodyown[j] = bodyown[i];
+}
+
+/* ----------------------------------------------------------------------
+ initialize one atom's array values, called when atom is created
+------------------------------------------------------------------------- */
+
+void FixRigidSmall::set_arrays(int i)
+{
+ bodyown[i] = -1;
+ bodytag[i] = 0;
+ atom2body[i] = -1;
+ displace[i][0] = 0.0;
+ displace[i][1] = 0.0;
+ displace[i][2] = 0.0;
+}
+
+/* ----------------------------------------------------------------------
+ pack values in local atom-based arrays for exchange with another proc
+------------------------------------------------------------------------- */
+
+int FixRigidSmall::pack_exchange(int i, double *buf)
+{
+ buf[0] = bodytag[i];
+ buf[1] = displace[i][0];
+ buf[2] = displace[i][1];
+ buf[3] = displace[i][2];
+
+ // extended attribute info
+
+ int m = 4;
+ if (extended) {
+ buf[m++] = eflags[i];
+ for (int j = 0; j < orientflag; j++)
+ buf[m++] = orient[i][j];
+ if (dorientflag) {
+ buf[m++] = dorient[i][0];
+ buf[m++] = dorient[i][1];
+ buf[m++] = dorient[i][2];
+ }
+ }
+
+ // atom not in a rigid body
+
+ if (!bodytag[i]) return m;
+
+ // atom does not own its rigid body
+
+ if (bodyown[i] < 0) {
+ buf[m++] = 0;
+ return m;
+ }
+
+ // body info for atom that owns a rigid body
+
+ buf[m++] = 1;
+ memcpy(&buf[m],&body[bodyown[i]],sizeof(Body));
+ m += bodysize;
+ return m;
+}
+
+/* ----------------------------------------------------------------------
+ unpack values in local atom-based arrays from exchange with another proc
+------------------------------------------------------------------------- */
+
+int FixRigidSmall::unpack_exchange(int nlocal, double *buf)
+{
+ bodytag[nlocal] = static_cast (buf[0]);
+ displace[nlocal][0] = buf[1];
+ displace[nlocal][1] = buf[2];
+ displace[nlocal][2] = buf[3];
+
+ // extended attribute info
+
+ int m = 4;
+ if (extended) {
+ eflags[nlocal] = static_cast (buf[m++]);
+ for (int j = 0; j < orientflag; j++)
+ orient[nlocal][j] = buf[m++];
+ if (dorientflag) {
+ dorient[nlocal][0] = buf[m++];
+ dorient[nlocal][1] = buf[m++];
+ dorient[nlocal][2] = buf[m++];
+ }
+ }
+
+ // atom not in a rigid body
+
+ if (!bodytag[nlocal]) return m;
+
+ // atom does not own its rigid body
+
+ bodyown[nlocal] = static_cast (buf[m++]);
+ if (bodyown[nlocal] == 0) {
+ bodyown[nlocal] = -1;
+ return m;
+ }
+
+ // body info for atom that owns a rigid body
+
+ if (nlocal_body == nmax_body) grow_body();
+ memcpy(&body[nlocal_body],&buf[m],sizeof(Body));
+ m += bodysize;
+ body[nlocal_body].ilocal = nlocal;
+ bodyown[nlocal] = nlocal_body++;
+
+ return m;
+}
+
+/* ----------------------------------------------------------------------
+ only pack body info if own or ghost atom owns the body
+ for FULL_BODY, send 0/1 flag with every atom
+------------------------------------------------------------------------- */
+
+int FixRigidSmall::pack_comm(int n, int *list, double *buf,
+ int pbc_flag, int *pbc)
+{
+ int i,j;
+ double *xcm,*vcm,*quat,*omega,*ex_space,*ey_space,*ez_space;
+
+ int nlocal = atom->nlocal;
+
+ int m = 0;
+
+ if (commflag == INITIAL) {
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ if (bodyown[j] < 0) continue;
+ xcm = body[bodyown[j]].xcm;
+ buf[m++] = xcm[0];
+ buf[m++] = xcm[1];
+ buf[m++] = xcm[2];
+ vcm = body[bodyown[j]].vcm;
+ buf[m++] = vcm[0];
+ buf[m++] = vcm[1];
+ buf[m++] = vcm[2];
+ quat = body[bodyown[j]].quat;
+ buf[m++] = quat[0];
+ buf[m++] = quat[1];
+ buf[m++] = quat[2];
+ buf[m++] = quat[3];
+ omega = body[bodyown[j]].omega;
+ buf[m++] = omega[0];
+ buf[m++] = omega[1];
+ buf[m++] = omega[2];
+ ex_space = body[bodyown[j]].ex_space;
+ buf[m++] = ex_space[0];
+ buf[m++] = ex_space[1];
+ buf[m++] = ex_space[2];
+ ey_space = body[bodyown[j]].ey_space;
+ buf[m++] = ey_space[0];
+ buf[m++] = ey_space[1];
+ buf[m++] = ey_space[2];
+ ez_space = body[bodyown[j]].ez_space;
+ buf[m++] = ez_space[0];
+ buf[m++] = ez_space[1];
+ buf[m++] = ez_space[2];
+ }
+
+ } else if (commflag == FINAL) {
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ if (bodyown[j] < 0) continue;
+ vcm = body[bodyown[j]].vcm;
+ buf[m++] = vcm[0];
+ buf[m++] = vcm[1];
+ buf[m++] = vcm[2];
+ omega = body[bodyown[j]].omega;
+ buf[m++] = omega[0];
+ buf[m++] = omega[1];
+ buf[m++] = omega[2];
+ }
+
+ } else if (commflag == FULL_BODY) {
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ if (bodyown[j] < 0) buf[m++] = 0;
+ else {
+ buf[m++] = 1;
+ memcpy(&buf[m],&body[bodyown[j]],sizeof(Body));
+ m += bodysize;
+ }
+ }
+ }
+
+ return m;
+}
+
+/* ----------------------------------------------------------------------
+ only ghost atoms are looped over
+ for FULL_BODY, store a new ghost body if this atom owns it
+ for other commflag values, only unpack body info if atom owns it
+------------------------------------------------------------------------- */
+
+void FixRigidSmall::unpack_comm(int n, int first, double *buf)
+{
+ int i,j,last,flag;
+ double *xcm,*vcm,*quat,*omega,*ex_space,*ey_space,*ez_space;
+
+ int m = 0;
+ last = first + n;
+
+ if (commflag == INITIAL) {
+ for (i = first; i < last; i++) {
+ if (bodyown[i] < 0) continue;
+ xcm = body[bodyown[i]].xcm;
+ xcm[0] = buf[m++];
+ xcm[1] = buf[m++];
+ xcm[2] = buf[m++];
+ vcm = body[bodyown[i]].vcm;
+ vcm[0] = buf[m++];
+ vcm[1] = buf[m++];
+ vcm[2] = buf[m++];
+ quat = body[bodyown[i]].quat;
+ quat[0] = buf[m++];
+ quat[1] = buf[m++];
+ quat[2] = buf[m++];
+ quat[3] = buf[m++];
+ omega = body[bodyown[i]].omega;
+ omega[0] = buf[m++];
+ omega[1] = buf[m++];
+ omega[2] = buf[m++];
+ ex_space = body[bodyown[i]].ex_space;
+ ex_space[0] = buf[m++];
+ ex_space[1] = buf[m++];
+ ex_space[2] = buf[m++];
+ ey_space = body[bodyown[i]].ey_space;
+ ey_space[0] = buf[m++];
+ ey_space[1] = buf[m++];
+ ey_space[2] = buf[m++];
+ ez_space = body[bodyown[i]].ez_space;
+ ez_space[0] = buf[m++];
+ ez_space[1] = buf[m++];
+ ez_space[2] = buf[m++];
+ }
+
+ } else if (commflag == FINAL) {
+ for (i = first; i < last; i++) {
+ if (bodyown[i] < 0) continue;
+ vcm = body[bodyown[i]].vcm;
+ vcm[0] = buf[m++];
+ vcm[1] = buf[m++];
+ vcm[2] = buf[m++];
+ omega = body[bodyown[i]].omega;
+ omega[0] = buf[m++];
+ omega[1] = buf[m++];
+ omega[2] = buf[m++];
+ }
+
+ } else if (commflag == FULL_BODY) {
+ for (i = first; i < last; i++) {
+ bodyown[i] = static_cast (buf[m++]);
+ if (bodyown[i] == 0) bodyown[i] = -1;
+ else {
+ j = nlocal_body + nghost_body;
+ if (j == nmax_body) grow_body();
+ memcpy(&body[j],&buf[m],sizeof(Body));
+ m += bodysize;
+ body[j].ilocal = i;
+ bodyown[i] = j;
+ nghost_body++;
+ }
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ only ghost atoms are looped over
+ only pack body info if atom owns it
+------------------------------------------------------------------------- */
+
+int FixRigidSmall::pack_reverse_comm(int n, int first, double *buf)
+{
+ int i,j,m,last;
+ double *fcm,*torque,*vcm,*angmom,*xcm,*itensor;
+
+ m = 0;
+ last = first + n;
+
+ if (commflag == FORCE_TORQUE) {
+ for (i = first; i < last; i++) {
+ if (bodyown[i] < 0) continue;
+ fcm = body[bodyown[i]].fcm;
+ buf[m++] = fcm[0];
+ buf[m++] = fcm[1];
+ buf[m++] = fcm[2];
+ torque = body[bodyown[i]].torque;
+ buf[m++] = torque[0];
+ buf[m++] = torque[1];
+ buf[m++] = torque[2];
+ }
+
+ } else if (commflag == VCM_ANGMOM) {
+ for (i = first; i < last; i++) {
+ if (bodyown[i] < 0) continue;
+ vcm = body[bodyown[i]].vcm;
+ buf[m++] = vcm[0];
+ buf[m++] = vcm[1];
+ buf[m++] = vcm[2];
+ angmom = body[bodyown[i]].angmom;
+ buf[m++] = angmom[0];
+ buf[m++] = angmom[1];
+ buf[m++] = angmom[2];
+ }
+
+ } else if (commflag == XCM_MASS) {
+ for (i = first; i < last; i++) {
+ if (bodyown[i] < 0) continue;
+ xcm = body[bodyown[i]].xcm;
+ buf[m++] = xcm[0];
+ buf[m++] = xcm[1];
+ buf[m++] = xcm[2];
+ buf[m++] = body[bodyown[i]].mass;
+ }
+
+ } else if (commflag == ITENSOR) {
+ for (i = first; i < last; i++) {
+ if (bodyown[i] < 0) continue;
+ itensor = body[bodyown[i]].itensor;
+ buf[m++] = itensor[0];
+ buf[m++] = itensor[1];
+ buf[m++] = itensor[2];
+ buf[m++] = itensor[3];
+ buf[m++] = itensor[4];
+ buf[m++] = itensor[5];
+ }
+
+ } else if (commflag == DOF) {
+ for (i = first; i < last; i++) {
+ if (bodyown[i] < 0) continue;
+ j = bodyown[i];
+ buf[m++] = counts[j][0];
+ buf[m++] = counts[j][1];
+ buf[m++] = counts[j][2];
+ }
+ }
+
+ return m;
+}
+
+/* ----------------------------------------------------------------------
+ only unpack body info if own or ghost atom owns the body
+------------------------------------------------------------------------- */
+
+void FixRigidSmall::unpack_reverse_comm(int n, int *list, double *buf)
+{
+ int i,j,k;
+ double *fcm,*torque,*vcm,*angmom,*xcm,*itensor;
+
+ int *tag = atom->tag;
+ int nlocal = atom->nlocal;
+
+ int m = 0;
+
+ if (commflag == FORCE_TORQUE) {
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ if (bodyown[j] < 0) continue;
+ fcm = body[bodyown[j]].fcm;
+ fcm[0] += buf[m++];
+ fcm[1] += buf[m++];
+ fcm[2] += buf[m++];
+ torque = body[bodyown[j]].torque;
+ torque[0] += buf[m++];
+ torque[1] += buf[m++];
+ torque[2] += buf[m++];
+ }
+
+ } else if (commflag == VCM_ANGMOM) {
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ if (bodyown[j] < 0) continue;
+ vcm = body[bodyown[j]].vcm;
+ vcm[0] += buf[m++];
+ vcm[1] += buf[m++];
+ vcm[2] += buf[m++];
+ angmom = body[bodyown[j]].angmom;
+ angmom[0] += buf[m++];
+ angmom[1] += buf[m++];
+ angmom[2] += buf[m++];
+ }
+
+ } else if (commflag == XCM_MASS) {
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ if (bodyown[j] < 0) continue;
+ xcm = body[bodyown[j]].xcm;
+ xcm[0] += buf[m++];
+ xcm[1] += buf[m++];
+ xcm[2] += buf[m++];
+ body[bodyown[j]].mass += buf[m++];
+ }
+
+ } else if (commflag == ITENSOR) {
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ if (bodyown[j] < 0) continue;
+ itensor = body[bodyown[j]].itensor;
+ itensor[0] += buf[m++];
+ itensor[1] += buf[m++];
+ itensor[2] += buf[m++];
+ itensor[3] += buf[m++];
+ itensor[4] += buf[m++];
+ itensor[5] += buf[m++];
+ }
+
+ } else if (commflag == DOF) {
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ if (bodyown[j] < 0) continue;
+ k = bodyown[j];
+ counts[k][0] += buf[m++];
+ counts[k][1] += buf[m++];
+ counts[k][2] += buf[m++];
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ grow body data structure
+------------------------------------------------------------------------- */
+
+void FixRigidSmall::grow_body()
+{
+ nmax_body += DELTA_BODY;
+ body = (Body *) memory->srealloc(body,nmax_body*sizeof(Body),
+ "rigid/small:body");
+}
+
+/* ----------------------------------------------------------------------
+ reset atom2body for all owned atoms
+ do this via bodyown of atom that owns the body the owned atom is in
+ atom2body values can point to original body or any image of the body
+------------------------------------------------------------------------- */
+
+void FixRigidSmall::reset_atom2body()
+{
+ int iowner;
+
+ // iowner = index of atom that owns the body that atom I is in
+
+ int nlocal = atom->nlocal;
+
+ for (int i = 0; i < nlocal; i++) {
+ atom2body[i] = -1;
+ if (bodytag[i]) {
+ iowner = atom->map(bodytag[i]);
+ if (iowner == -1) {
+ char str[128];
+ sprintf(str,
+ "Rigid body atoms %d %d missing on proc %d at step "
+ BIGINT_FORMAT,
+ atom->tag[i],bodytag[i],comm->me,update->ntimestep);
+ error->one(FLERR,str);
+
+ }
+ atom2body[i] = bodyown[iowner];
+ }
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixRigidSmall::reset_dt()
+{
+ dtv = update->dt;
+ dtf = 0.5 * update->dt * force->ftm2v;
+ dtq = 0.5 * update->dt;
+}
+
+/* ----------------------------------------------------------------------
+ debug method for sanity checking of atom/body data pointers
+------------------------------------------------------------------------- */
+
+void FixRigidSmall::check(int flag)
+{
+ /*
+ for (int i = 0; i < atom->nlocal; i++) {
+ if (bodyown[i] >= 0) {
+ if (bodytag[i] != atom->tag[i]) {
+ printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
+ error->one(FLERR,"BAD AAA");
+ }
+ if (bodyown[i] < 0 || bodyown[i] >= nlocal_body) {
+ printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
+ error->one(FLERR,"BAD BBB");
+ }
+ if (atom2body[i] != bodyown[i]) {
+ printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
+ error->one(FLERR,"BAD CCC");
+ }
+ if (body[bodyown[i]].ilocal != i) {
+ printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
+ error->one(FLERR,"BAD DDD");
+ }
+ }
+ }
+
+ for (int i = 0; i < atom->nlocal; i++) {
+ if (bodyown[i] < 0 && bodytag[i] > 0) {
+ if (atom2body[i] < 0 || atom2body[i] >= nlocal_body+nghost_body) {
+ printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
+ error->one(FLERR,"BAD EEE");
+ }
+ if (bodytag[i] != atom->tag[body[atom2body[i]].ilocal]) {
+ printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
+ error->one(FLERR,"BAD FFF");
+ }
+ }
+ }
+
+ for (int i = atom->nlocal; i < atom->nlocal + atom->nghost; i++) {
+ if (bodyown[i] >= 0) {
+ if (bodyown[i] < nlocal_body ||
+ bodyown[i] >= nlocal_body+nghost_body) {
+ printf("Values %d %d: %d %d %d\n",
+ i,atom->tag[i],bodyown[i],nlocal_body,nghost_body);
+ printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
+ error->one(FLERR,"BAD GGG");
+ }
+ if (body[bodyown[i]].ilocal != i) {
+ printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
+ error->one(FLERR,"BAD HHH");
+ }
+ }
+ }
+
+ for (int i = 0; i < nlocal_body; i++) {
+ if (body[i].ilocal < 0 || body[i].ilocal >= atom->nlocal) {
+ printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
+ error->one(FLERR,"BAD III");
+ }
+ if (bodytag[body[i].ilocal] != atom->tag[body[i].ilocal] ||
+ bodyown[body[i].ilocal] != i) {
+ printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
+ error->one(FLERR,"BAD JJJ");
+ }
+ }
+
+ for (int i = nlocal_body; i < nlocal_body + nghost_body; i++) {
+ if (body[i].ilocal < atom->nlocal ||
+ body[i].ilocal >= atom->nlocal + atom->nghost) {
+ printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
+ error->one(FLERR,"BAD KKK");
+ }
+ if (bodyown[body[i].ilocal] != i) {
+ printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
+ error->one(FLERR,"BAD LLL");
+ }
+ }
+ */
+}
diff --git a/src/RIGID/fix_rigid_small.h b/src/RIGID/fix_rigid_small.h
new file mode 100644
index 0000000000..7ddd49139d
--- /dev/null
+++ b/src/RIGID/fix_rigid_small.h
@@ -0,0 +1,252 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+
+FixStyle(rigid/small,FixRigidSmall)
+
+#else
+
+#ifndef LMP_FIX_RIGID_SMALL_H
+#define LMP_FIX_RIGID_SMALL_H
+
+#include "fix.h"
+
+// replace this later
+#include