whitespace

doc/src/dump_image.rst

@@ -783,109 +783,134 @@ x-component of velocity if the atom-attribute "vx" was specified.
 The basic idea of a color map is that the atom-attribute will be
 within a range of values, and that range is associated with a series
 of colors (e.g. red, blue, green).  An atom's specific value (vx =
--3.2) can then mapped to the series of colors (e.g. halfway between
-red and blue), and a specific color is determined via an interpolation
-procedure.
+-3.2) can then be mapped to a specific color.  The details of the
+mapping procedure depend on the *style* of the color map, as explained
+below.
 
-There are many possible options for the color map, enabled by the
-*amap* keyword.  Here are the details.
+There are many possible options for defining a color map, enabled by
+the *amap* keyword.  Here are the details.
 
 The *lo* and *hi* settings determine the range of values allowed for
 the atom attribute.  If numeric values are used for *lo* and/or *hi*,
-then values that are lower/higher than that value are set to the
-value.  I.e. the range is static.  If *lo* is specified as *min* or
-*hi* as *max* then the range is dynamic, and the lower and/or
-upper bound will be calculated each time an image is drawn, based
-on the set of atoms being visualized.
+then individual atom values which are lower/higher than lo/hi are set
+to lo/hi for purposes of determining the atom's color.  I.e. the range
+is static.  If *lo* is specified as *min* or *hi* as *max* then the
+range is dynamic.  The lower and/or upper bound will be calculated
+each time an image is drawn, based on the current atom values of all
+the atoms being visualized.
 
 The *style* setting is two letters, such as "ca".  The first letter is
 either "c" for continuous, "d" for discrete, or "s" for sequential.
 The second letter is either "a" for absolute, or "f" for fractional.
 
-A continuous color map is one in which the color changes continuously
-from value to value within the range.  A discrete color map is one in
-which discrete colors are assigned to sub-ranges of values within the
-range.  A sequential color map is one in which discrete colors are
-assigned to a sequence of sub-ranges of values covering the entire
-range.
+A *continuous* color map is one in which the color of an atom changes
+continuously as its attribute value increases within the range.
+Colors are assigned to specific values within the range; an atom with
+an attribue value between two adjacent specific values is assigned a
+color interpolated between the two adjacent colors.
 
-An absolute color map is one in which the values to which colors are
-assigned are specified explicitly as values within the range.  A
-fractional color map is one in which the values to which colors are
-assigned are specified as a fractional portion of the range.  For
-example if the range is from -10.0 to 10.0, and the color red is to be
-assigned to atoms with a value of 5.0, then for an absolute color map
-the number 5.0 would be used.  But for a fractional map, the number
-0.75 would be used since 5.0 is 3/4 of the way from -10.0 to 10.0.
+A *discrete* color map is one in which discrete colors are assigned to
+sub-ranges of values within the overall range.  Each sub-range can be
+of variable width and overlap with other sub-ranges.  An atom with an
+attribute value is mapped to one of the sub-ranges and assigned that
+color.
+
+A *sequential* color map is similar to a discrete color map except that
+all sub-ranges are of equal width and discrete colors are assigned to
+each sub-range in a round-robin fashion until the overall range is
+covered from *lo* to *hi*.
+
+An *absolute* color map is one in which the numeric settings
+associated with assigned colors are specified explicitly as values
+within the range.
+
+A *fractional* color map is one in which the numeric settings
+associated with assigned colors are specified as a fractional position
+within the range.
+
+For a continuous color map, the numeric settings are the specific
+values each color is assigned to.  For a discrete color map, the
+numeric settings are the bounds of each sub-range.  For a sequential
+color map, the numeric settings is the width of all the sub-ranges.
+For example if the overall range is from -10.0 to 10.0, and the color
+red is to be assigned to atoms with an attribute value of 5.0, then
+for an absolute color map the numeric value of 5.0 should map to red.
+But for a fractional map, the numeric value of 0.75 should map to red,
+since 5.0 is 3/4 of the way from -10.0 to 10.0.
 
 The *delta* setting must be specified for all styles, but is only used
-for the sequential style; otherwise the value is ignored.  It
-specifies the bin size to use within the range for assigning
-consecutive colors to.  For example, if the range is from :math:`-10.0` to
-:math:`10.0` and a *delta* of :math:`1.0` is used, then 20 colors will be
-assigned to the range.  The first will be from
-:math:`-10.0 \le \text{color1} < -9.0`, then second from
+for the *sequential* style; otherwise the setting is ignored.  It
+specifies the bin size of the sub-ranges of values described above,
+For example, if the overall range is from :math:`-10.0` to
+:math:`10.0` and a *delta* of :math:`1.0` is used, then 20 colors will
+be assigned to a series of sub-ranges.  The first sub-range will be
+from :math:`-10.0 \le \text{color1} < -9.0`, the second from
 :math:`-9.0 \le color2 < -8.0`, etc.
 
-The *N* setting is how many entries follow.  The format of the entries
-depends on whether the color map style is continuous, discrete or
-sequential.  In all cases the *color* setting can be any of the 140
-pre-defined colors (see below) or a color name defined by the
-dump_modify color option.
+The *N* setting is how many color entries follow.  The format of each
+color entry depends on whether the color map style is continuous,
+discrete, or sequential.  For each entry, the specified *color* can be
+any of the 140 pre-defined colors (see below) or a color name defined
+by the dump_modify color option.
 
-For continuous color maps, each entry has a *value* and a *color*\ .
-The *value* is either a number within the range of values or *min* or
-*max*\ .  The *value* of the first entry must be *min* and the *value*
-of the last entry must be *max*\ .  Any entries in between must have
-increasing values.  Note that numeric values can be specified either
-as absolute numbers or as fractions (0.0 to 1.0) of the range,
-depending on the "a" or "f" in the style setting for the color map.
+For *continuous* color maps, each entry has a *value* and a *color*\ .
+The *value* is either a number within the *lo/hi* range of values or
+*min* or *max*\ .  The *value* for the first entry must be *min* and
+the *value* for the last entry must be *max*\ .  In-between entries
+must have increasing numeric values.  There must be 2 or more entries.
+Note that numeric values are specified either as absolute numbers or
+as fractions (0.0 to 1.0) of the range, depending on the "a" or "f" in
+the style setting for the color map.
 
-Here is how the entries are used to determine the color of an individual
-atom, given the value :math:`X` of its atom attribute.  :math:`X` will
-fall between 2 of the entry values.  The color of the atom is linearly
-interpolated (in each of the RGB values) between the 2 colors associated
-with those entries.  For example, if :math:`X = -5.0` and the two
-surrounding entries are "red" at :math:`-10.0` and "blue" at
-:math:`0.0`, then the atom's color will be halfway between "red" and
-"blue", which happens to be "purple".
+Here is how the *N* entries are used to determine the color of an
+individual atom, based on the value :math:`X` of its atom attribute.
+:math:`X` will fall between 2 of the entry values.  The color of the
+atom is linearly interpolated (in each of the RGB values) between the
+2 colors associated with those entries.  For example, if :math:`X =
+-5.0` and the two surrounding entries are "red" at :math:`-10.0` and
+"blue" at :math:`0.0`, then the atom's color will be halfway between
+"red" and "blue", which in this case is "purple".
 
-For discrete color maps, each entry has a *lo* and *hi* value and a
+For *discrete* color maps, each entry has a *lo* and *hi* value and a
 *color*\ .  The *lo* and *hi* settings are either numbers within the
-range of values or *lo* can be *min* or *hi* can be *max*\ .  The *lo*
-and *hi* settings of the last entry must be *min* and *max*\ .  Other
+range of values or *min* (for *lo) or *max* (for *hi*).  The *lo* and
+*hi* settings of the last entry must be *min* and *max*\ .  Other
 entries can have any *lo* and *hi* values and the sub-ranges of
-different values can overlap.  Note that numeric *lo* and *hi* values
-can be specified either as absolute numbers or as fractions (0.0 to 1.0)
-of the range, depending on the "a" or "f" in the style setting for the
-color map.
+different entries can overlap.  There must be one or more entries.
+Note that numeric *lo* and *hi* values are specified either as
+absolute numbers or as fractions (0.0 to 1.0) of the range, depending
+on the "a" or "f" in the style setting for the color map.  The lo/hi
+values in each sub-range must satisfy lo < hi.
 
-Here is how the entries are used to determine the color of an individual
-atom, given the value X of its atom attribute.  The entries are scanned
-from first to last.  The first time that *lo* <= X <= *hi*, X is
-assigned the color associated with that entry.  You can think of the
-last entry as assigning a default color (since it will always be matched
-by X), and the earlier entries as colors that override the default.
-Also note that no interpolation of a color RGB is done.  All atoms will
-be drawn with one of the colors in the list of entries.
+Here is how the *N* entries are used to determine the color of an
+individual atom, based on the value :math:`X` of its atom attribute.
+The entries are scanned from first to last.  The first time that *lo*
+<= X <= *hi*, X is assigned the color associated with that entry.
+This means can the last entry can be thought of as a default color
+(since it will always be matched by X); the earlier entries override
+the default.  Note that for a *discrete* map, no interpolation of a
+color RGB values is done.  All atoms will be drawn with one of the
+colors in the list of entries.
 
-For sequential color maps, each entry has only a *color*\ .  Here is how
-the entries are used to determine the color of an individual atom,
-given the value X of its atom attribute.  The range is partitioned
-into N bins of width *binsize*\ .  Thus X will fall in a specific bin
-from 1 to N, say the Mth bin.  If it falls on a boundary between 2
-bins, it is considered to be in the higher of the 2 bins.  Each bin is
-assigned a color from the E entries.  If E < N, then the colors are
-repeated.  For example if 2 entries with colors red and green are
-specified, then the odd numbered bins will be red and the even bins
-green.  The color of the atom is the color of its bin.  Note that the
-sequential color map is really a shorthand way of defining a discrete
-color map without having to specify where all the bin boundaries are.
+For *sequential* color maps, each entry has only a *color*\ .  There
+must be 1 or more entries.  Here is how the *N* entries are used to
+determine the color of an individual atom, given the value X of its
+atom attribute.  The range is overlaid with M bins of width *delta*\ ,
+the last of which may extend beyond the *hi* boundary of the range.
+Thus X will fall in a specific bin from 1 to M.  If it falls on a
+boundary between 2 bins, it is considered to be in the higher of the 2
+bins (except in the case of 2 bins whose boundary is the *hi*
+boundary, it is considered to be in the lower of the 2 bins).  Each of
+the M bins is assigned a color from the *N* entries.  If M < *N*, then
+the colors are repeated in a round-robin fashion.  For example if 2
+entries with colors red and green are specified, then the odd numbered
+bins will be red and the even bins green.  An atom's color is the
+color of its bin.
 
 Here is an example of using a sequential color map to color all the
-atoms in individual molecules with a different color.  See the
-examples/pour/in.pour.2d.molecule input script for an example of how
-this is used.
+atoms in individual molecules with a different color based on its
+molecule ID.  See the examples/pour/in.pour.2d.molecule input script
+for an example of how this is used.
 
 .. code-block:: LAMMPS
 
@@ -897,8 +922,9 @@ this is used.
                    zoom 1.6 adiam 1.5
    dump_modify     1 pad 5 amap 0 10 sa 1 10 ${colors}
 
-In this case, 10 colors are defined, and molecule IDs are
-mapped to one of the colors, even if there are 1000s of molecules.
+In this case, the sequential color map has 10 color entries, and
+molecule IDs are mapped to one of the colors, even if there are 1000s
+of molecules.
 
 ----------
 

src/.gitignore

@@ -677,6 +677,8 @@
 /compute_rattlers_atom.h
 /compute_reaxff_atom.cpp
 /compute_reaxff_atom.h
+/compute_rigid_atom.cpp
+/compute_rigid_atom.h
 /compute_rigid_local.cpp
 /compute_rigid_local.h
 /compute_slcsa_atom.cpp

src/KOKKOS/atom_vec_sphere_kokkos.cpp

@@ -1141,13 +1141,12 @@ struct AtomVecSphereKokkos_PackBorderVel {
     _buf(i,6) = _radius(j);
     _buf(i,7) = _rmass(j);
     if (DEFORM_VREMAP) {
-      if (_mask(i) & _deform_groupbit) {
+      if (_mask(j) & _deform_groupbit) {
         _buf(i,8) = _v(j,0) + _dvx;
         _buf(i,9) = _v(j,1) + _dvy;
         _buf(i,10) = _v(j,2) + _dvz;
       }
-    }
-    else {
+    } else {
       _buf(i,8) = _v(j,0);
       _buf(i,9) = _v(j,1);
       _buf(i,10) = _v(j,2);

src/RIGID/compute_rigid_atom.cpp

@@ -0,0 +1,305 @@
+// clang-format off
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   LAMMPS development team: developers@lammps.org
+
+   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 "compute_rigid_atom.h"
+
+#include "atom.h"
+#include "domain.h"
+#include "error.h"
+#include "modify.h"
+#include "fix_rigid_small.h"
+#include "memory.h"
+#include "update.h"
+
+#include <cstring>
+
+using namespace LAMMPS_NS;
+
+static constexpr int DELTA = 10000;
+
+enum{ID,MOL,MASS,X,Y,Z,XU,YU,ZU,VX,VY,VZ,FX,FY,FZ,IX,IY,IZ,
+     TQX,TQY,TQZ,OMEGAX,OMEGAY,OMEGAZ,ANGMOMX,ANGMOMY,ANGMOMZ,
+     QUATW,QUATI,QUATJ,QUATK,INERTIAX,INERTIAY,INERTIAZ};
+
+/* ---------------------------------------------------------------------- */
+
+ComputeRigidAtom::ComputeRigidAtom(LAMMPS *lmp, int narg, char **arg) :
+  Compute(lmp, narg, arg),
+  rstyle(nullptr), idrigid(nullptr), fixrigid(nullptr)
+{
+  if (narg < 5) error->all(FLERR,"Illegal compute rigid/local command");
+
+  peratom_flag = 1;
+
+  idrigid = utils::strdup(arg[3]);
+
+  nvalues = narg - 4;
+  rstyle = new int[nvalues];
+
+  nvalues = 0;
+  for (int iarg = 4; iarg < narg; iarg++) {
+    if (strcmp(arg[iarg],"id") == 0) rstyle[nvalues++] = ID;
+    else if (strcmp(arg[iarg],"mol") == 0) rstyle[nvalues++] = MOL;
+    else if (strcmp(arg[iarg],"mass") == 0) rstyle[nvalues++] = MASS;
+    else if (strcmp(arg[iarg],"x") == 0) rstyle[nvalues++] = X;
+    else if (strcmp(arg[iarg],"y") == 0) rstyle[nvalues++] = Y;
+    else if (strcmp(arg[iarg],"z") == 0) rstyle[nvalues++] = Z;
+    else if (strcmp(arg[iarg],"xu") == 0) rstyle[nvalues++] = XU;
+    else if (strcmp(arg[iarg],"yu") == 0) rstyle[nvalues++] = YU;
+    else if (strcmp(arg[iarg],"zu") == 0) rstyle[nvalues++] = ZU;
+    else if (strcmp(arg[iarg],"vx") == 0) rstyle[nvalues++] = VX;
+    else if (strcmp(arg[iarg],"vy") == 0) rstyle[nvalues++] = VY;
+    else if (strcmp(arg[iarg],"vz") == 0) rstyle[nvalues++] = VZ;
+    else if (strcmp(arg[iarg],"fx") == 0) rstyle[nvalues++] = FX;
+    else if (strcmp(arg[iarg],"fy") == 0) rstyle[nvalues++] = FY;
+    else if (strcmp(arg[iarg],"fz") == 0) rstyle[nvalues++] = FZ;
+    else if (strcmp(arg[iarg],"ix") == 0) rstyle[nvalues++] = IX;
+    else if (strcmp(arg[iarg],"iy") == 0) rstyle[nvalues++] = IY;
+    else if (strcmp(arg[iarg],"iz") == 0) rstyle[nvalues++] = IZ;
+    else if (strcmp(arg[iarg],"tqx") == 0) rstyle[nvalues++] = TQX;
+    else if (strcmp(arg[iarg],"tqy") == 0) rstyle[nvalues++] = TQY;
+    else if (strcmp(arg[iarg],"tqz") == 0) rstyle[nvalues++] = TQZ;
+    else if (strcmp(arg[iarg],"omegax") == 0) rstyle[nvalues++] = OMEGAX;
+    else if (strcmp(arg[iarg],"omegay") == 0) rstyle[nvalues++] = OMEGAY;
+    else if (strcmp(arg[iarg],"omegaz") == 0) rstyle[nvalues++] = OMEGAZ;
+    else if (strcmp(arg[iarg],"angmomx") == 0) rstyle[nvalues++] = ANGMOMX;
+    else if (strcmp(arg[iarg],"angmomy") == 0) rstyle[nvalues++] = ANGMOMY;
+    else if (strcmp(arg[iarg],"angmomz") == 0) rstyle[nvalues++] = ANGMOMZ;
+    else if (strcmp(arg[iarg],"quatw") == 0) rstyle[nvalues++] = QUATW;
+    else if (strcmp(arg[iarg],"quati") == 0) rstyle[nvalues++] = QUATI;
+    else if (strcmp(arg[iarg],"quatj") == 0) rstyle[nvalues++] = QUATJ;
+    else if (strcmp(arg[iarg],"quatk") == 0) rstyle[nvalues++] = QUATK;
+    else if (strcmp(arg[iarg],"inertiax") == 0) rstyle[nvalues++] = INERTIAX;
+    else if (strcmp(arg[iarg],"inertiay") == 0) rstyle[nvalues++] = INERTIAY;
+    else if (strcmp(arg[iarg],"inertiaz") == 0) rstyle[nvalues++] = INERTIAZ;
+    else error->all(FLERR,"Invalid keyword in compute rigid/local command");
+  }
+
+  if (nvalues == 1) size_peratom_cols = 0;
+  else size_peratom_cols = nvalues;
+
+  maxatom = 0;
+  vector_atom = nullptr;
+  array_atom = nullptr;
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeRigidAtom::~ComputeRigidAtom()
+{
+  memory->destroy(vector_atom);
+  memory->destroy(array_atom);
+  delete[] idrigid;
+  delete[] rstyle;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeRigidAtom::init()
+{
+  // set fixrigid
+
+  auto ifix = modify->get_fix_by_id(idrigid);
+  if (!ifix) error->all(FLERR,"FixRigidSmall ID {} for compute rigid/atom does not exist", idrigid);
+  fixrigid = dynamic_cast<FixRigidSmall *>(ifix);
+  if (!fixrigid)
+    error->all(FLERR,"Fix ID {} for compute rigid/atom does not point to fix rigid/small", idrigid);
+
+  // do initial memory allocation so that memory_usage() is correct
+
+  if (atom->nmax > maxatom) {
+    maxatom = atom->nmax;
+    if (nvalues == 1) {
+      memory->destroy(vector_atom);
+      memory->create(vector_atom, maxatom, "rigid/atom:vector_atom");
+    } else {
+      memory->destroy(array_atom);
+      memory->create(array_atom, maxatom, nvalues, "rigid/atom:array_atom");
+    }
+  }
+}
+
+/* ----------------------------------------------------------------------
+   loop over atoms, store rigid body info for the body it is in
+   set output to zero if atom not in group or not in a body
+------------------------------------------------------------------------- */
+
+void ComputeRigidAtom::compute_peratom()
+{
+  invoked_peratom = update->ntimestep;
+
+  int i,m,n,ibody;
+  double *ptr;
+  FixRigidSmall::Body *body;
+
+  double xprd = domain->xprd;
+  double yprd = domain->yprd;
+  double zprd = domain->zprd;
+
+  tagint *tag = atom->tag;
+  tagint *molecule = atom->molecule;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  m = 0;
+  for (i = 0; i < nlocal; i++) {
+    if (!(mask[i] & groupbit)) {
+      zero(i);
+      continue;
+    }
+
+    ibody = fixrigid->atom2body[i];
+    if (ibody < 0) {
+      zero(i);
+      continue;
+    }
+
+    body = &fixrigid->body[ibody];
+
+    if (nvalues == 1) ptr = &vector_atom[i];
+    else ptr = array_atom[i];
+
+    for (n = 0; n < nvalues; n++) {
+      switch (rstyle[n]) {
+      case ID:
+        ptr[n] = tag[body->ilocal];
+        break;
+      case MOL:
+        ptr[n] = molecule[body->ilocal];
+        break;
+      case MASS:
+        ptr[n] = body->mass;
+        break;
+      case X:
+        ptr[n] = body->xcm[0];
+        break;
+      case Y:
+        ptr[n] = body->xcm[1];
+        break;
+      case Z:
+        ptr[n] = body->xcm[2];
+        break;
+      case XU:
+        ptr[n] = body->xcm[0] +
+          ((body->image & IMGMASK) - IMGMAX) * xprd;
+        break;
+      case YU:
+        ptr[n] = body->xcm[1] +
+          ((body->image >> IMGBITS & IMGMASK) - IMGMAX) * yprd;
+        break;
+      case ZU:
+        ptr[n] = body->xcm[2] +
+          ((body->image >> IMG2BITS) - IMGMAX) * zprd;
+        break;
+      case VX:
+        ptr[n] = body->vcm[0];
+        break;
+      case VY:
+        ptr[n] = body->vcm[1];
+        break;
+      case VZ:
+        ptr[n] = body->vcm[2];
+        break;
+      case FX:
+        ptr[n] = body->fcm[0];
+        break;
+      case FY:
+        ptr[n] = body->fcm[1];
+        break;
+      case FZ:
+        ptr[n] = body->fcm[2];
+        break;
+      case IX:
+        ptr[n] = (body->image & IMGMASK) - IMGMAX;
+        break;
+      case IY:
+        ptr[n] = (body->image >> IMGBITS & IMGMASK) - IMGMAX;
+        break;
+      case IZ:
+        ptr[n] = (body->image >> IMG2BITS) - IMGMAX;
+        break;
+      case TQX:
+        ptr[n] = body->torque[0];
+        break;
+      case TQY:
+        ptr[n] = body->torque[1];
+        break;
+      case TQZ:
+        ptr[n] = body->torque[2];
+        break;
+      case OMEGAX:
+        ptr[n] = body->omega[0];
+        break;
+      case OMEGAY:
+        ptr[n] = body->omega[1];
+        break;
+      case OMEGAZ:
+        ptr[n] = body->omega[2];
+        break;
+      case ANGMOMX:
+        ptr[n] = body->angmom[0];
+        break;
+      case ANGMOMY:
+        ptr[n] = body->angmom[1];
+        break;
+      case ANGMOMZ:
+        ptr[n] = body->angmom[2];
+        break;
+      case QUATW:
+        ptr[n] = body->quat[0];
+        break;
+      case QUATI:
+        ptr[n] = body->quat[1];
+        break;
+      case QUATJ:
+        ptr[n] = body->quat[2];
+        break;
+      case QUATK:
+        ptr[n] = body->quat[3];
+        break;
+      case INERTIAX:
+        ptr[n] = body->inertia[0];
+        break;
+      case INERTIAY:
+        ptr[n] = body->inertia[1];
+        break;
+      case INERTIAZ:
+        ptr[n] = body->inertia[2];
+        break;
+      }
+    }
+  }
+}
+
+/* ----------------------------------------------------------------------
+   zero the peratom output for atom I
+------------------------------------------------------------------------- */
+
+void ComputeRigidAtom::zero(int i)
+{
+  if (nvalues == 1) vector_atom[i] = 0.0;
+  else {
+    for (int n = 0; n < nvalues; n++) array_atom[i][n] = 0.0;
+  }
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of peratom data
+------------------------------------------------------------------------- */
+
+double ComputeRigidAtom::memory_usage()
+{
+  double bytes = (double)maxatom*nvalues * sizeof(double);
+  return bytes;
+}

src/RIGID/compute_rigid_atom.h

@@ -0,0 +1,49 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   LAMMPS development team: developers@lammps.org
+
+   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 COMPUTE_CLASS
+// clang-format off
+ComputeStyle(rigid/atom,ComputeRigidAtom);
+// clang-format on
+#else
+
+#ifndef LMP_COMPUTE_RIGID_ATOM_H
+#define LMP_COMPUTE_RIGID_ATOM_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeRigidAtom : public Compute {
+ public:
+  ComputeRigidAtom(class LAMMPS *, int, char **);
+  ~ComputeRigidAtom() override;
+  void init() override;
+  void compute_peratom() override;
+  double memory_usage() override;
+
+ private:
+  int nvalues;
+  int *rstyle;
+
+  int maxatom;
+  char *idrigid;
+  class FixRigidSmall *fixrigid;
+
+  void zero(int);
+};
+
+}    // namespace LAMMPS_NS
+
+#endif
+#endif

src/RIGID/compute_rigid_local.cpp

@@ -89,7 +89,7 @@ ComputeRigidLocal::ComputeRigidLocal(LAMMPS *lmp, int narg, char **arg) :
   if (nvalues == 1) size_local_cols = 0;
   else size_local_cols = nvalues;
 
-ncount = nmax = 0;
+  ncount = nmax = 0;
   vlocal = nullptr;
   alocal = nullptr;
 }

src/RIGID/fix_rigid.cpp

@@ -969,20 +969,24 @@ void FixRigid::initial_integrate(int vflag)
    use domain->remap() in case xcm is far away from box
      due to first-time definition of rigid body in setup_bodies_static()
      or due to box flip
+   also remap vcm if xcm crosses periodic shearing boundary
    also adjust imagebody = rigid body image flags, due to xcm remap
    also reset body xcmimage flags of all atoms in bodies
    xcmimage flags are relative to xcm so that body can be unwrapped
    if don't do this, would need xcm to move with true image flags
      then a body could end up very far away from box
-     set_xv() will then compute huge displacements every step to
-       reset coords of all body atoms to be back inside the box,
-       ditto for triclinic box flip, which causes numeric problems
+   set_xv() would then compute huge displacements every step to
+     reset coords of all body atoms to be back inside the box,
+     ditto for triclinic box flip which could cause numeric problems
 ------------------------------------------------------------------------- */
 
 void FixRigid::pre_neighbor()
 {
-  for (int ibody = 0; ibody < nbody; ibody++)
-    domain->remap(xcm[ibody],imagebody[ibody]);
+  for (int ibody = 0; ibody < nbody; ibody++) {
+    //domain->remap(xcm[ibody],imagebody[ibody]);
+    domain->remap(xcm[ibody],imagebody[ibody],vcm[ibody]);
+  }
+
   image_shift();
 }
 

src/RIGID/fix_rigid_small.cpp

@@ -21,6 +21,8 @@
 #include "comm.h"
 #include "domain.h"
 #include "error.h"
+#include "fix.h"
+#include "fix_deform.h"
 #include "force.h"
 #include "group.h"
 #include "hashlittle.h"
@@ -559,6 +561,16 @@ void FixRigidSmall::init()
     if (ifix->box_change) boxflag = true;
   }
 
+  // check for fix deform with V_REMAP set
+
+  deform_vremap = 0;
+  const auto &fixes = modify->get_fix_list();
+  for (const auto &fix : fixes)
+    if (utils::strmatch(fix->style,"^deform")) {
+      if ((dynamic_cast<FixDeform *>(fix))->remapflag == Domain::V_REMAP)
+        deform_vremap = 1;
+    }
+
   // add gravity forces based on gravity vector from fix
 
   if (id_gravity) {
@@ -787,7 +799,8 @@ void FixRigidSmall::initial_integrate(int vflag)
      due to first-time definition of rigid body in setup_bodies_static()
      or due to box flip
    also adjust imagebody = rigid body image flags, due to xcm remap
-   then communicate bodies so other procs will know of changes to body xcm
+   also remap vcm if xcm crosses periodic shearing boundary
+   then communicate bodies so other procs will know of changes to body xcm/vcm
    then adjust xcmimage flags of all atoms in bodies via image_shift()
      for two effects
      (1) change in true image flags due to pbc() call during exchange
@@ -795,18 +808,17 @@ void FixRigidSmall::initial_integrate(int vflag)
    xcmimage flags are always -1,0,-1 so that body can be unwrapped
      around in-box xcm and stay close to simulation box
    if just inferred unwrapped from atom image flags,
-     then a body could end up very far away
-     when unwrapped by true image flags
-   then set_xv() will compute huge displacements every step to reset coords of
-     all the body atoms to be back inside the box, ditto for triclinic box flip
-     note: so just want to avoid that numeric problem?
+     then an unwrapped body could end up very far away from box
+   set_xv() would then compute huge displacements every step to
+     reset coords of all body atoms to be back inside the box,
+     ditto for triclinic box flip which could cause numeric problems
 ------------------------------------------------------------------------- */
 
 void FixRigidSmall::pre_neighbor()
 {
   for (int ibody = 0; ibody < nlocal_body; ibody++) {
     Body *b = &body[ibody];
-    domain->remap(b->xcm,b->image);
+    domain->remap(b->xcm,b->image,b->vcm);
   }
 
   nghost_body = 0;
@@ -966,9 +978,11 @@ void FixRigidSmall::apply_langevin_thermostat()
     gamma1 = -body[ibody].mass / t_period / ftm2v;
     gamma2 = sqrt(body[ibody].mass) * tsqrt *
       sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
+    if (deform_vremap) remove_bias(ibody,vcm);
     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);
+    if (deform_vremap) restore_bias(ibody,vcm);
 
     gamma1 = -1.0 / t_period / ftm2v;
     gamma2 = tsqrt * sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
@@ -992,6 +1006,37 @@ void FixRigidSmall::apply_langevin_thermostat()
   }
 }
 
+/* ----------------------------------------------------------------------
+   remove velocity bias from VCM of Body ibody to leave thermal VCM
+------------------------------------------------------------------------- */
+
+void FixRigidSmall::remove_bias(int ibody, double *vcm)
+{
+  double lamda[3];
+  double *h_rate = domain->h_rate;
+  double *h_ratelo = domain->h_ratelo;
+
+  domain->x2lamda(body[ibody].xcm, lamda);
+  vbias[0] = h_rate[0] * lamda[0] + h_rate[5] * lamda[1] + h_rate[4] * lamda[2] + h_ratelo[0];
+  vbias[1] = h_rate[1] * lamda[1] + h_rate[3] * lamda[2] + h_ratelo[1];
+  vbias[2] = h_rate[2] * lamda[2] + h_ratelo[2];
+  vcm[0] -= vbias[0];
+  vcm[1] -= vbias[1];
+  vcm[2] -= vbias[2];
+}
+
+/* ----------------------------------------------------------------------
+   add back velocity bias to VCM of Body ibody removed by remove_bias()
+   assume remove_bias() was previously called
+------------------------------------------------------------------------- */
+
+void FixRigidSmall::restore_bias(int /*i*/, double *vcm)
+{
+  vcm[0] += vbias[0];
+  vcm[1] += vbias[1];
+  vcm[2] += vbias[2];
+}
+
 /* ---------------------------------------------------------------------- */
 
 void FixRigidSmall::compute_forces_and_torques()

src/RIGID/fix_rigid_small.h

@@ -26,6 +26,7 @@ namespace LAMMPS_NS {
 
 class FixRigidSmall : public Fix {
   friend class ComputeRigidLocal;
+  friend class ComputeRigidAtom;
 
  public:
   FixRigidSmall(class LAMMPS *, int, char **);
@@ -151,6 +152,9 @@ class FixRigidSmall : public Fix {
   double **langextra;                  // Langevin thermostat forces and torques
   int maxlang;                         // max size of langextra
   class RanMars *random;               // RNG
+  int deform_vremap;                   // 1 if fix deform with V_REMAP exists
+                                       //   if so, add/sub bias around Langevin
+  double vbias[3];                     // store deformation bias for one body
 
   int tstat_flag, pstat_flag;    // 0/1 = no/yes thermostat/barostat
 
@@ -200,6 +204,8 @@ class FixRigidSmall : public Fix {
   void setup_bodies_static();
   void setup_bodies_dynamic();
   void apply_langevin_thermostat();
+  void remove_bias(int, double *);
+  void restore_bias(int, double *);
   virtual void compute_forces_and_torques();
   void enforce2d();
   void readfile(int, double **, int *);

src/atom_vec.cpp

@@ -485,7 +485,7 @@ int AtomVec::pack_comm_vel(int n, int *list, double *buf, int pbc_flag, int *pbc
         buf[m++] = x[j][0] + dx;
         buf[m++] = x[j][1] + dy;
         buf[m++] = x[j][2] + dz;
-        if (mask[i] & deform_groupbit) {
+        if (mask[j] & deform_groupbit) {
           buf[m++] = v[j][0] + dvx;
           buf[m++] = v[j][1] + dvy;
           buf[m++] = v[j][2] + dvz;
@@ -946,7 +946,7 @@ int AtomVec::pack_border_vel(int n, int *list, double *buf, int pbc_flag, int *p
         buf[m++] = ubuf(tag[j]).d;
         buf[m++] = ubuf(type[j]).d;
         buf[m++] = ubuf(mask[j]).d;
-        if (mask[i] & deform_groupbit) {
+        if (mask[j] & deform_groupbit) {
           buf[m++] = v[j][0] + dvx;
           buf[m++] = v[j][1] + dvy;
           buf[m++] = v[j][2] + dvz;

src/domain.cpp

@@ -1527,16 +1527,19 @@ void Domain::closest_image(const double * const xi, const double * const xj, dou
 }
 
 /* ----------------------------------------------------------------------
-   remap the point into the periodic box no matter how far away
+   remap the point X into the periodic box no matter how far away
    adjust 3 image flags encoded in image accordingly
    resulting coord must satisfy lo <= coord < hi
    MAX is important since coord - prd < lo can happen when coord = hi
    for triclinic, point is converted to lamda coords (0-1) before doing remap
-   image = 10 bits for each dimension
+   image = 10 or 20 bits for each dimension
    increment/decrement in wrap-around fashion
+   if V is specified (default = NULL) and deform_vremap set by fix deform:
+     also remap v via h_rate calclated by fix deform
+     currently only used by fix rigid commands to remap body VCM
 ------------------------------------------------------------------------- */
 
-void Domain::remap(double *x, imageint &image)
+void Domain::remap(double *x, imageint &image, double *v)
 {
   double *lo,*hi,*period,*coord;
   double lamda[3];
@@ -1558,6 +1561,7 @@ void Domain::remap(double *x, imageint &image)
   if (xperiodic) {
     while (coord[0] < lo[0]) {
       coord[0] += period[0];
+      if (deform_vremap && v) v[0] += h_rate[0];
       idim = image & IMGMASK;
       otherdims = image ^ idim;
       idim--;
@@ -1566,6 +1570,7 @@ void Domain::remap(double *x, imageint &image)
     }
     while (coord[0] >= hi[0]) {
       coord[0] -= period[0];
+      if (deform_vremap && v) v[0] -= h_rate[0];
       idim = image & IMGMASK;
       otherdims = image ^ idim;
       idim++;
@@ -1578,6 +1583,10 @@ void Domain::remap(double *x, imageint &image)
   if (yperiodic) {
     while (coord[1] < lo[1]) {
       coord[1] += period[1];
+      if (deform_vremap && v) {
+        v[0] += h_rate[5];
+        v[1] += h_rate[1];
+      }
       idim = (image >> IMGBITS) & IMGMASK;
       otherdims = image ^ (idim << IMGBITS);
       idim--;
@@ -1586,6 +1595,10 @@ void Domain::remap(double *x, imageint &image)
     }
     while (coord[1] >= hi[1]) {
       coord[1] -= period[1];
+      if (deform_vremap && v) {
+        v[0] -= h_rate[5];
+        v[1] -= h_rate[1];
+      }
       idim = (image >> IMGBITS) & IMGMASK;
       otherdims = image ^ (idim << IMGBITS);
       idim++;
@@ -1598,6 +1611,11 @@ void Domain::remap(double *x, imageint &image)
   if (zperiodic) {
     while (coord[2] < lo[2]) {
       coord[2] += period[2];
+      if (deform_vremap && v) {
+        v[0] += h_rate[4];
+        v[1] += h_rate[3];
+        v[2] += h_rate[2];
+      }
       idim = image >> IMG2BITS;
       otherdims = image ^ (idim << IMG2BITS);
       idim--;
@@ -1606,6 +1624,11 @@ void Domain::remap(double *x, imageint &image)
     }
     while (coord[2] >= hi[2]) {
       coord[2] -= period[2];
+      if (deform_vremap && v) {
+        v[0] -= h_rate[4];
+        v[1] -= h_rate[3];
+        v[2] -= h_rate[2];
+      }
       idim = image >> IMG2BITS;
       otherdims = image ^ (idim << IMG2BITS);
       idim++;
@@ -1619,7 +1642,7 @@ void Domain::remap(double *x, imageint &image)
 }
 
 /* ----------------------------------------------------------------------
-   remap the point into the periodic box no matter how far away
+   remap the point X into the periodic box no matter how far away
    no image flag calculation
    resulting coord must satisfy lo <= coord < hi
    MAX is important since coord - prd < lo can happen when coord = hi
@@ -1666,7 +1689,7 @@ void Domain::remap(double *x)
 }
 
 /* ----------------------------------------------------------------------
-   remap all points into the periodic box no matter how far away
+   remap all atom coords into the periodic box no matter how far away
    adjust 3 image flags encoded in image accordingly
    resulting coord must satisfy lo <= coord < hi
    MAX is important since coord - prd < lo can happen when coord = hi

src/domain.h

@@ -134,7 +134,7 @@ class Domain : protected Pointers {
   int closest_image(int, int);
   int closest_image(const double *const, int);
   void closest_image(const double *const, const double *const, double *const);
-  void remap(double *, imageint &);
+  void remap(double *, imageint &, double *v = nullptr);
   void remap(double *);
   void remap_all();
   void remap_near(double *, double *);

src/fix_ave_chunk.cpp

@@ -271,7 +271,7 @@ FixAveChunk::FixAveChunk(LAMMPS *lmp, int narg, char **arg) :
         error->all(FLERR, val.iarg, "Fix ave/chunk compute {} does not calculate a per-atom array",
                    val.id);
       if (val.argindex && (val.argindex > val.val.c->size_peratom_cols))
-        error->all(FLERR, val.iarg, "Fix ave/chunk compute {} vector is accessed out-of-range{}",
+        error->all(FLERR, val.iarg, "Fix ave/chunk compute {} array is accessed out-of-range{}",
                    val.id, utils::errorurl(20));
 
     } else if (val.which == ArgInfo::FIX) {
@@ -288,7 +288,7 @@ FixAveChunk::FixAveChunk(LAMMPS *lmp, int narg, char **arg) :
         error->all(FLERR, val.iarg, "Fix ave/chunk fix {} does not calculate a per-atom array",
                    val.id);
       if (val.argindex && val.argindex > val.val.f->size_peratom_cols)
-        error->all(FLERR, val.iarg, "Fix ave/chunk fix {} vector is accessed out-of-range{}",
+        error->all(FLERR, val.iarg, "Fix ave/chunk fix {} array is accessed out-of-range{}",
                    val.id, utils::errorurl(20));
     } else if (val.which == ArgInfo::VARIABLE) {
       val.val.v = input->variable->find(val.id.c_str());

src/image.cpp

@@ -1785,7 +1785,7 @@ ColorMap::~ColorMap()
 }
 
 /* ----------------------------------------------------------------------
-   redefine color map
+   redefine a color map
    args = lo hi style delta N entry1 entry2 ... entryN as defined by caller
    return 1 if any error in args, else return 0
 ------------------------------------------------------------------------- */
@@ -1849,13 +1849,11 @@ int ColorMap::reset(int narg, char **arg)
         mentry[i].lo = NUMERIC;
         mentry[i].lvalue = utils::numeric(FLERR,arg[n],false,lmp);
       } else if (strcmp(arg[n],"min") == 0) mentry[i].lo = MINVALUE;
-      else if (strcmp(arg[n],"max") == 0) mentry[i].lo = MAXVALUE;
       else return 1;
       if (!islower(arg[n+1][0])) {
         mentry[i].hi = NUMERIC;
         mentry[i].hvalue = utils::numeric(FLERR,arg[n+1],false,lmp);
-      } else if (strcmp(arg[n+1],"min") == 0) mentry[i].hi = MINVALUE;
-      else if (strcmp(arg[n+1],"max") == 0) mentry[i].hi = MAXVALUE;
+      } else if (strcmp(arg[n+1],"max") == 0) mentry[i].hi = MAXVALUE;
       else return 1;
       mentry[i].color = image->color2rgb(arg[n+2]);
       n += 3;
@@ -1867,7 +1865,7 @@ int ColorMap::reset(int narg, char **arg)
       // current code is just 1st nentry values of ALL or USER
       // need to comment out error check in DumpImage::modify_param()
       //   for amap check on (narg < n) to get it to work
-      // need to add extra logic here to check not accessing undefined colors
+      // need to add extra logic here to ensure not accessinging undefined colors
       if (i == 0) {
         if (n+1 > narg) return 1;
         if (strcmp(arg[n],"ALL") == 0) expandflag = 1;
@@ -1883,6 +1881,9 @@ int ColorMap::reset(int narg, char **arg)
       }
       n += 1;
     }
+
+    // error return if color string was not recognized
+
     if (mentry[i].color == nullptr) return 1;
   }
 
@@ -1890,8 +1891,11 @@ int ColorMap::reset(int narg, char **arg)
     if (nentry < 2) return 1;
     if (mentry[0].single != MINVALUE || mentry[nentry-1].single != MAXVALUE)
       return 1;
-    for (int i = 2; i < nentry-1; i++)
+    for (int i = 1; i < nentry-1; i++)
+      if (mentry[i].single != NUMERIC) return 1;
+    for (int i = 2; i < nentry-1; i++) {
       if (mentry[i].svalue <= mentry[i-1].svalue) return 1;
+    }
   } else if (mstyle == DISCRETE) {
     if (nentry < 1) return 1;
     if (mentry[nentry-1].lo != MINVALUE || mentry[nentry-1].hi != MAXVALUE)
@@ -1930,7 +1934,7 @@ int ColorMap::minmax(double mindynamic, double maxdynamic)
     if (mrange == ABSOLUTE) mentry[nentry-1].svalue = hicurrent;
     else mentry[nentry-1].svalue = 1.0;
 
-    // error in ABSOLUTE mode if new lo/hi current cause
+    // error in ABSOLUTE mode if current lo/hi values cause
     // first/last entry to become lo > hi with adjacent entry
 
     if (mrange == ABSOLUTE) {
@@ -1952,6 +1956,23 @@ int ColorMap::minmax(double mindynamic, double maxdynamic)
         else mentry[i].hvalue = 1.0;
       }
     }
+
+  // set rounddown_flag if bin boundary is essentially at hicurrent
+  // this is to prevent all atoms with values >= hicurrent being assigned
+  //   the color for a bin starting at highcurrent,
+  //   more sensible to assign them color for bin ending at hicurrent
+  // rounddown_flag is used in value2color()
+
+  } else if (mstyle == SEQUENTIAL) {
+    // NOTE: need to formalize and pre-set these 2 epsilon constants
+    double epsbin;
+    if (mrange == ABSOLUTE) epsbin = (hicurrent-locurrent) * 1.0e-12;
+    else epsbin = 1.0e-9;
+
+    int ibin = static_cast<int> ((hicurrent-locurrent) * mbinsizeinv);
+    int jbin = static_cast<int> ((hicurrent-locurrent-epsbin) * mbinsizeinv);
+    if (jbin < ibin) rounddown_flag = 1;
+    else rounddown_flag = 0;
   }
 
   return 0;
@@ -1964,7 +1985,7 @@ int ColorMap::minmax(double mindynamic, double maxdynamic)
 
 double *ColorMap::value2color(double value)
 {
-  double lo;//,hi;
+  double lo,hi;
 
   value = MAX(value,locurrent);
   value = MIN(value,hicurrent);
@@ -1973,10 +1994,10 @@ double *ColorMap::value2color(double value)
     if (locurrent == hicurrent) value = 0.0;
     else value = (value-locurrent) / (hicurrent-locurrent);
     lo = 0.0;
-    //hi = 1.0;
+    hi = 1.0;
   } else {
     lo = locurrent;
-    //hi = hicurrent;
+    hi = hicurrent;
   }
 
   if (mstyle == CONTINUOUS) {
@@ -1998,6 +2019,7 @@ double *ColorMap::value2color(double value)
         return mentry[i].color;
   } else {
     int ibin = static_cast<int> ((value-lo) * mbinsizeinv);
+    if (value == hi && rounddown_flag) ibin--;
     return mentry[ibin%nentry].color;
   }
 

src/image.h

@@ -165,8 +165,13 @@ class ColorMap : protected Pointers {
   class Image *image;              // caller with color2rgb() method
   int mstyle, mrange;              // 2-letter style/range of color map
   int mlo, mhi;                    // bounds = NUMERIC or MINVALUE or MAXVALUE
+  int rounddown_flag;              // for sequential color map,
+                                   // ensure value at hi end of range is not
+                                   //   assigned color of bin starting at hi,
+                                   //   but rather the color of bin ending at hi
   double mlovalue, mhivalue;       // user bounds if NUMERIC
   double locurrent, hicurrent;     // current bounds for this snapshot
+                                   //   from user bounds or dynamic bounds
   double mbinsize, mbinsizeinv;    // bin size for sequential color map
   double interpolate[3];           // local storage for returned RGB color