From e0a528a42da5f705bcda90e40fbb9c80dd9b6d7c Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Thu, 24 Apr 2025 13:20:00 -0600 Subject: [PATCH 1/2] =?UTF-8?q?rigid-shear:=20squash=20onto=20develop=20?= =?UTF-8?q?=E2=80=94=20by=20wrhamilt?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit (cherry picked from commit ac69b58f0b073c230bb444d8c5dfea1efb3f2d01) --- doc/src/dump_image.rst | 202 +++++++++-------- src/.gitignore | 2 + src/KOKKOS/atom_vec_sphere_kokkos.cpp | 5 +- src/RIGID/compute_rigid_atom.cpp | 305 ++++++++++++++++++++++++++ src/RIGID/compute_rigid_atom.h | 49 +++++ src/RIGID/fix_rigid.cpp | 14 +- src/RIGID/fix_rigid_small.cpp | 59 ++++- src/RIGID/fix_rigid_small.h | 6 + src/atom_vec.cpp | 4 +- src/domain.cpp | 33 ++- src/domain.h | 2 +- src/fix_ave_chunk.cpp | 4 +- src/image.cpp | 41 +++- src/image.h | 5 + 14 files changed, 608 insertions(+), 123 deletions(-) create mode 100644 src/RIGID/compute_rigid_atom.cpp create mode 100644 src/RIGID/compute_rigid_atom.h diff --git a/doc/src/dump_image.rst b/doc/src/dump_image.rst index f513630b93..03f74d8dea 100644 --- a/doc/src/dump_image.rst +++ b/doc/src/dump_image.rst @@ -828,109 +828,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. - -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 *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. + +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. - -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. - -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". - -For discrete color maps, each entry has a *lo* and *hi* value and a +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 *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 *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 *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. - -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. - -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. +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 *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*\ . 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 @@ -942,8 +967,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. ---------- diff --git a/src/.gitignore b/src/.gitignore index c22710678d..5a1f3cc573 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -703,6 +703,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 diff --git a/src/KOKKOS/atom_vec_sphere_kokkos.cpp b/src/KOKKOS/atom_vec_sphere_kokkos.cpp index 85912b2579..fa46e6af2a 100644 --- a/src/KOKKOS/atom_vec_sphere_kokkos.cpp +++ b/src/KOKKOS/atom_vec_sphere_kokkos.cpp @@ -1148,13 +1148,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); diff --git a/src/RIGID/compute_rigid_atom.cpp b/src/RIGID/compute_rigid_atom.cpp new file mode 100644 index 0000000000..2f5e0091b1 --- /dev/null +++ b/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 + +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(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; +} diff --git a/src/RIGID/compute_rigid_atom.h b/src/RIGID/compute_rigid_atom.h new file mode 100644 index 0000000000..e581ec8cdf --- /dev/null +++ b/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 diff --git a/src/RIGID/fix_rigid.cpp b/src/RIGID/fix_rigid.cpp index d18c733401..2af66d2bce 100644 --- a/src/RIGID/fix_rigid.cpp +++ b/src/RIGID/fix_rigid.cpp @@ -970,20 +970,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(); } diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp index 43b48779ec..bbf3fe6592 100644 --- a/src/RIGID/fix_rigid_small.cpp +++ b/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 "input.h" @@ -561,6 +563,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(fix))->remapflag == Domain::V_REMAP) + deform_vremap = 1; + } + // add gravity forces based on gravity vector from fix if (id_gravity) { @@ -789,7 +801,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 @@ -797,18 +810,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; @@ -968,9 +980,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; @@ -994,6 +1008,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() diff --git a/src/RIGID/fix_rigid_small.h b/src/RIGID/fix_rigid_small.h index c176912998..5befe75314 100644 --- a/src/RIGID/fix_rigid_small.h +++ b/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 **); @@ -153,6 +154,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 @@ -202,6 +206,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 *); diff --git a/src/atom_vec.cpp b/src/atom_vec.cpp index ef01447547..567d05af79 100644 --- a/src/atom_vec.cpp +++ b/src/atom_vec.cpp @@ -487,7 +487,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; @@ -948,7 +948,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; diff --git a/src/domain.cpp b/src/domain.cpp index b06529e56a..2cd089612c 100644 --- a/src/domain.cpp +++ b/src/domain.cpp @@ -1538,16 +1538,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]; @@ -1569,6 +1572,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--; @@ -1577,6 +1581,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++; @@ -1589,6 +1594,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--; @@ -1597,6 +1606,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++; @@ -1609,6 +1622,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--; @@ -1617,6 +1635,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++; @@ -1630,7 +1653,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 @@ -1677,7 +1700,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 diff --git a/src/domain.h b/src/domain.h index 545fcf7ad1..9202c899ac 100644 --- a/src/domain.h +++ b/src/domain.h @@ -140,7 +140,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 *); diff --git a/src/fix_ave_chunk.cpp b/src/fix_ave_chunk.cpp index ed556f85b7..5f922809b6 100644 --- a/src/fix_ave_chunk.cpp +++ b/src/fix_ave_chunk.cpp @@ -272,7 +272,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) { @@ -289,7 +289,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()); diff --git a/src/image.cpp b/src/image.cpp index 4c65d0b25a..4040504946 100644 --- a/src/image.cpp +++ b/src/image.cpp @@ -1782,7 +1782,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 ------------------------------------------------------------------------- */ @@ -1846,13 +1846,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; @@ -1864,7 +1862,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; @@ -1880,6 +1878,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; } @@ -1887,8 +1888,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) @@ -1927,7 +1931,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) { @@ -1949,6 +1953,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 ((hicurrent-locurrent) * mbinsizeinv); + int jbin = static_cast ((hicurrent-locurrent-epsbin) * mbinsizeinv); + if (jbin < ibin) rounddown_flag = 1; + else rounddown_flag = 0; } return 0; @@ -1961,7 +1982,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); @@ -1970,10 +1991,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) { diff --git a/src/image.h b/src/image.h index 71d671c041..8669593049 100644 --- a/src/image.h +++ b/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 From 2158bedbcf7ad80897cf60586919d99bedff7438 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 31 Aug 2025 23:04:34 -0400 Subject: [PATCH 2/2] add missing override (cherry picked from commit 9e5500b734b6e8d1dffd61059086cb7fe1461eba) (cherry picked from commit 809ea52d6a984f819d4924c483fc4ff0b6d629bc)