diff --git a/doc/src/Howto_viz.rst b/doc/src/Howto_viz.rst index 80ec18c8020..517fc025484 100644 --- a/doc/src/Howto_viz.rst +++ b/doc/src/Howto_viz.rst @@ -995,7 +995,7 @@ input file: group ogroup type 1 group hgroup type 2 - compute hb all hbond/local 3.3 30.0 ogroup ogroup hgroup + compute hb all hbond/local 3.5 30.0 ogroup ogroup hgroup dump viz all image 100 water-*.png element type size 600 600 zoom 1.331 view 70 20 & shiny 0.2 ssao yes 348276 0.6 fsaa yes box yes 0.025 & @@ -1017,25 +1017,19 @@ and the surrounding water molecules in both directions. .. code-block:: LAMMPS - # define ellipsoid region around peptide for hbond analysis and visualization - group peptide type <= 12 - variable comx equal xcm(peptide,x) - variable comy equal xcm(peptide,y) - variable comz equal xcm(peptide,z) - region shell ellipsoid v_comx v_comy v_comz 7.0 8.0 16.0 - group viz dynamic all region shell include molecule - - # define groups of donor and acceptor atoms for peptide and water - group pdonor type 5 - group wdonor type 13 - group pacceptor type 3 5 9 12 - group wacceptor type 13 - group hydrogen type 10 14 + # select atoms for visualization: peptide and water molecules within 3.5 angstrom + group viz dynamic peptide within 3.5 include molecule every 100 + + # define groups of donor, acceptor, and hydrogen atoms for peptide and water + group pdonor type 5 9 # peptide donors : nitrogens and phenol oxygen + group woxygen type 13 # water oxygens are donor and acceptor + group pacceptor type 3 5 9 12 # peptide acceptors: oxygens, nitrogens, and sulfur + group hydrogen type 4 10 14 # hydrogens bonded to oxygens and nitrogens # peptide-water hydrogen bonds where the peptide is the donor - compute hb1 all hbond/local 3.3 30.0 pdonor wacceptor hydrogen + compute hb1 all hbond/local 3.5 30.0 pdonor woxygen hydrogen # peptide-water hydrogen bonds where the peptide is the acceptor - compute hb2 all hbond/local 3.7 30.0 wdonor pacceptor hydrogen + compute hb2 all hbond/local 3.7 30.0 woxygen pacceptor hydrogen # create donor/acceptor hydrogen bond info text fix label all graphics/labels 100 text "Hydrogen bonds donated: $(c_hb1:%02.0f)" 207 72 0.0 & diff --git a/doc/src/group.rst b/doc/src/group.rst index 02d43c5e7e1..12ff3141053 100644 --- a/doc/src/group.rst +++ b/doc/src/group.rst @@ -38,13 +38,15 @@ Syntax *intersect* args = two or more group IDs *dynamic* args = parent-ID keyword value ... one or more keyword/value pairs may be appended - keyword = *region* or *var* or *property* or *every* + keyword = *region* or *var* or *property* or *every* or *exclude* or *include* or *within* *region* value = region-ID *var* value = name of variable *property* value = name of custom integer or floating point vector *every* value = N = update group every this many timesteps + *exclude* value = group-ID = exclude atoms in this group from dynamic group *include* args = molecule - molecule = add atoms to group with same molecule ID as atoms already in group + molecule = add atoms to the dynamic group with same molecule ID as atoms already selected for the dynamic group + *within* value = cutoff = add atoms to the dynamic group that are within the given cutoff distance of the atoms already selected for the dynamic group *static* = no args Examples @@ -67,6 +69,7 @@ Examples group boundary intersect upper flow group boundary delete group mine dynamic all region myRegion every 100 + group solvation dynamic peptide within 3.5 include molecule every 100 exclude peptide Description """"""""""" @@ -227,10 +230,12 @@ added to the specified group. The *dynamic* style flags an existing or new group as dynamic. This means atoms will be (re)assigned to the group periodically as a simulation runs. This is in contrast to static groups where atoms are -permanently assigned to the group. The way the assignment occurs is -as follows. Only atoms in the group specified as the parent group via -the parent-ID are assigned to the dynamic group before the following -conditions are applied. +permanently assigned to the group. The way the assignment occurs is as +follows. Only atoms in the group specified as the parent group via the +parent-ID are initially assigned to the dynamic group before the +following selection conditions are applied. The *include* and *within* +keywords can then trigger including additional atoms which may be +outside the parent group. If the *region* keyword is used, atoms not in the specified region are removed from the dynamic group. @@ -243,43 +248,69 @@ If the *property* keyword is used, the name refers to a custom integer or floating point per-atom vector defined via the :doc:`fix property/atom ` command. This means the values in the vector can be read as part of a data file with the :doc:`read_data -` command or specified with the :doc:`set ` command. -Or accessed and changed via the :doc:`library interface to LAMMPS +` command or specified with the :doc:`set ` command. Or +accessed and changed via the :doc:`library interface to LAMMPS `, or by styles you add to LAMMPS (pair, fix, compute, -etc) which access the custom vector and modify its values. Which -means the values can be modified between or during simulations. Atoms -whose values in the custom vector are zero are removed from the -dynamic group. Note that the name of the custom per-atom vector is -specified just as *name*, not as *i_name* or *d_name* as it is for -other commands that use different kinds of custom atom vectors or -arrays as arguments. +etc) which access the custom vector and modify its values. Which means +the values can be modified between or during simulations. Atoms whose +values in the custom vector are zero are removed from the dynamic group. +Note that the name of the custom per-atom vector is specified just as +*name*, not as *i_name* or *d_name* as it is for other commands that use +different kinds of custom atom vectors or arrays as arguments. The assignment of atoms to a dynamic group is done at the beginning of -each run and on every timestep that is a multiple of *N*\ , which is -the argument for the *every* keyword (:math:`N = 1` is the default). For an +each run and on every timestep that is a multiple of *N*\ , which is the +argument for the *every* keyword (:math:`N = 1` is the default). For an energy minimization, via the :doc:`minimize ` command, an -assignment is made at the beginning of the minimization, but not -during the iterations of the minimizer. - -The optional *include* keyword with its arg *molecule* adds atoms to a -dynamic group that have the same molecule ID as atoms already in the -group. The molecule ID = 0 is ignored in this operation, since it is -assumed to flag isolated atoms that are not part of molecules. An -example of where this operation is useful is if the dynamic group is -defined using a *region*. If molecules straddle the region boundary, -then atoms outside the region that are part of molecules with atoms -inside the region will be added to the group (see above). - -The point in the timestep at which atoms are assigned to a dynamic -group is after interatomic forces have been computed, but before any +assignment is made at the beginning of the minimization, but not during +the iterations of the minimizer. + +.. versionadded:: TBD + +The optional *exclude* keyword has a group ID as argument and after all +other atom selections for the dynamic group have been performed, it +removes selected atoms that are in the specified group from the dynamic +group. This can for example be used to first select additional solvent +atoms up to a given distance from a group using the *within* keyword +(see below), but then exclude the non-solvent atoms. + +.. versionadded:: 11Feb2026 + +The optional *include* keyword with its argument *molecule* adds atoms +to a dynamic group that have the same molecule ID as atoms already in +the group. Atoms with the molecule ID = 0 are ignored in this +operation, since that ID is generally assumed to flag isolated atoms +that are not part of molecules. An example of where this operation is +useful is if the dynamic group is defined using a *region*. If +molecules straddle the region boundary, then atoms outside the region +that are part of molecules with atoms inside the region will be added to +the group (see above). + +.. versionadded:: TBD + +The optional *within* keyword with its *cutoff* argument adds atoms to a +dynamic group that are within the given cutoff distance from any of the +atoms selected for the dynamic group. Since this requires processing a +neighbor list and computing distances for potentially many pairs of +atoms plus and additional forward communication step, this can be an +expensive operation that should be used in combination with the *every* +option and a value of (e.g. 10 or larger) to reduce the overhead for +updating the selection. Outside of flow simulations and similar, atoms +will mostly oscillate around their positions and diffuse on a much +larger timescale, so infrequent group updates will not be a problem. +Also when selecting atoms for output to a :doc:`dump file ` it is +not necessary to update the selection more often than the dump interval. + +The point in the timestep at which atoms are assigned to a dynamic group +is **after** interatomic forces have been computed, but **before** any fixes which alter forces or otherwise update the system have been invoked. This means that atom positions have been updated, neighbor lists and ghost atoms are current, and both intermolecular and -intramolecular forces have been calculated based on the new -coordinates. Thus the region criterion, if applied, should be -accurate. Also, any computes invoked by an atom-style variable should -use updated information for that timestep (e.g., potential energy/atom -or coordination number/atom). Similarly, fixes or computes which are +intramolecular forces have been calculated based on the new coordinates. +Thus the region criterion, if applied, should be accurate. Also, any +computes invoked by an atom-style variable should use updated +information for that timestep (e.g., potential energy/atom or +coordination number/atom). Similarly, fixes or computes which are invoked after that point in the timestep, should operate on the new group of atoms. diff --git a/examples/GRAPHICS/in.peptide-hbonds b/examples/GRAPHICS/in.peptide-hbonds index 228f54e671c..fbaf19c752d 100644 --- a/examples/GRAPHICS/in.peptide-hbonds +++ b/examples/GRAPHICS/in.peptide-hbonds @@ -1,4 +1,9 @@ -# Solvated 5-mer peptide +# Visualize 5-mer peptide with solvation shell and hydrogen bonds + +# output dump image every this many steps +variable vizsteps index 100 +# run for this many steps +variable runsteps index 20000 units real atom_style full @@ -23,40 +28,53 @@ thermo 1000 fix 1 all nvt temp 275.0 275.0 100.0 tchain 1 fix 2 all shake 0.0001 10 1000 b 4 6 8 10 12 14 18 a 31 -# define region around peptide for hbond analysis and visualization group peptide type <= 12 + +# compute center of mass of peptide and center of box and displace atoms so they will match variable comx equal xcm(peptide,x) variable comy equal xcm(peptide,y) variable comz equal xcm(peptide,z) -region water ellipsoid v_comx v_comy v_comz 7.0 8.0 16.0 -group viz dynamic all region water include molecule +variable midx equal 0.5*(xhi+xlo) +variable midy equal 0.5*(yhi+ylo) +variable midz equal 0.5*(zhi+zlo) +variable xshift equal -v_comx+v_midx +variable yshift equal -v_comy+v_midy +variable zshift equal -v_comz+v_midz +displace_atoms all move v_xshift v_yshift v_zshift + +# select atoms for visualization: peptide and water molecules within 3.5 angstrom +group viz dynamic peptide within 3.5 include molecule every ${vizsteps} # define groups of donor and acceptor atoms for peptide and water -group pdonor type 5 9 -group wdonor type 13 -group pacceptor type 3 5 9 12 -group wacceptor type 13 -group hydrogen type 10 14 +group pdonor type 5 9 # peptide donors : nitrogens and phenol oxygen +group woxygen type 13 # water oxygens are donor and acceptor +group pacceptor type 3 5 9 12 # peptide acceptors: oxygens, nitrogens, and sulfur +group hydrogen type 4 10 14 # hydrogens bonded to oxygens and nitrogens # peptide-water hydrogen bonds where the peptide is the donor -compute hb1 all hbond/local 3.3 30.0 pdonor wacceptor hydrogen +compute hb1 all hbond/local 3.5 30.0 pdonor woxygen hydrogen # peptide-water hydrogen bonds where the peptide is the acceptor -compute hb2 all hbond/local 3.7 30.0 wdonor pacceptor hydrogen +compute hb2 all hbond/local 3.7 30.0 woxygen pacceptor hydrogen # create donor/acceptor hydrogen bonds info text -fix label all graphics/labels 100 text " Hydrogen bonds donated: $(c_hb1:%02.0f)" 207 72 0.0 & - size 24 backcolor darkgray & - text "Hydrogen bonds accepted: $(c_hb2:%02.0f)" 210 30 0.0 & - size 24 backcolor darkgray +fix label all graphics/labels ${vizsteps} text " Hydrogen bonds donated: $(c_hb1:%02.0f)" 207 72 0.0 & + size 24 backcolor darkgray & + text "Hydrogen bonds accepted: $(c_hb2:%02.0f)" 210 30 0.0 & + size 24 backcolor darkgray # create colored arrows to go with text labels -fix obj all graphics/objects 100 arrow 5 80.0 61.0 39 80.0 67.0 39 0.3 0.2 & - arrow 6 80.0 61.0 37.2 80.0 67.0 37.2 0.3 0.2 +fix obj all graphics/objects ${vizsteps} arrow 5 80.0 57.1 27.5 80.0 65.0 27.5 0.35 0.2 & + arrow 6 80.0 57.1 25.7 80.0 65.0 25.7 0.35 0.2 +# add progress bar +variable prog equal step/${runsteps} +fix graphics all graphics/objects ${vizsteps} progbar 10 13 z $(xhi) $(ylo-0.5) $(zlo+0.5*(zhi-zlo)) $(lz-3.0) 0.5 v_prog 10 # combine the graphics into visualization -dump viz viz image 100 hbonds-*.png element element size 600 600 zoom 2.1 view 80 0 bond atom 0.3 & - fsaa yes ssao yes 12384 0.6 shiny 0.1 box no 0.1 center s 0.5 0.52 0.52 & - compute hb1 const -0.4 0.3 compute hb2 const -0.4 0.3 fix label const 1 0 fix obj type 0.0 0.0 +dump viz viz image ${vizsteps} hbonds-*.png element element size 800 800 zoom 1.7 view 90 0 & + fsaa yes ssao yes 12384 0.6 shiny 0.1 box no 0.1 bond atom 0.3 & + center s 0.5 0.5 0.4 axes no 0.5 0.1 & + compute hb1 const -0.4 0.4 compute hb2 const -0.4 0.4 & + fix label const 1 0 fix obj type 0.0 0.0 fix graphics element 1 0.1 dump_modify viz pad 5 boxcolor white backcolor darkgray backcolor2 silver & element C C O H N C C C O H H S O H ccolor hb1 cyan ccolor hb2 magenta -run 20000 +run ${runsteps} diff --git a/src/fix_group.cpp b/src/fix_group.cpp index 3c72d4ae071..8074054151e 100644 --- a/src/fix_group.cpp +++ b/src/fix_group.cpp @@ -21,6 +21,9 @@ #include "input.h" #include "memory.h" #include "modify.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "neighbor.h" #include "region.h" #include "respa.h" #include "update.h" @@ -34,7 +37,8 @@ using namespace FixConst; /* ---------------------------------------------------------------------- */ FixGroup::FixGroup(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), idregion(nullptr), idvar(nullptr), idprop(nullptr), region(nullptr) + Fix(lmp, narg, arg), idregion(nullptr), idvar(nullptr), idprop(nullptr), idexclude(nullptr), + region(nullptr), list(nullptr) { // gbit = bitmask of dynamic group // group ID is last part of fix ID @@ -51,14 +55,25 @@ FixGroup::FixGroup(LAMMPS *lmp, int narg, char **arg) : varflag = 0; propflag = 0; moleculeflag = 0; + withinflag = 0; + excludeflag = 0; + cutoff = 0.0; nevery = 1; + excludebit = 0; + + int ioffset = 0; + if (lmp->input->arg) { + for (int i = 0; i < lmp->input->narg; ++i) + if (lmp->input->arg[i] == arg[3]) ioffset = i; + } int iarg = 3; while (iarg < narg) { if (strcmp(arg[iarg], "region") == 0) { if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "group dynamic region", error); if (!domain->get_region_by_id(arg[iarg + 1])) - error->all(FLERR, "Region {} for dynamic group {} does not exist", arg[iarg + 1], dgroupid); + error->all(FLERR, ioffset + iarg + 1, "Region {} for dynamic group {} does not exist", + arg[iarg + 1], dgroupid); regionflag = 1; delete[] idregion; idregion = utils::strdup(arg[iarg + 1]); @@ -67,8 +82,8 @@ FixGroup::FixGroup(LAMMPS *lmp, int narg, char **arg) : } else if (strcmp(arg[iarg], "var") == 0) { if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "group dynamic var", error); if (input->variable->find(arg[iarg + 1]) < 0) - error->all(FLERR, "Variable '{}' for dynamic group {} does not exist", arg[iarg + 1], - dgroupid); + error->all(FLERR, ioffset + iarg + 1, "Variable '{}' for dynamic group {} does not exist", + arg[iarg + 1], dgroupid); varflag = 1; delete[] idvar; idvar = utils::strdup(arg[iarg + 1]); @@ -79,8 +94,9 @@ FixGroup::FixGroup(LAMMPS *lmp, int narg, char **arg) : int flag, cols; iprop = atom->find_custom(arg[iarg + 1], flag, cols); if (iprop < 0 || cols) - error->all(FLERR, "Custom per-atom vector {} for dynamic group {} does not exist", - arg[iarg + 1], dgroupid); + error->all(FLERR, ioffset + iarg + 1, + "Custom per-atom vector {} for dynamic group {} does not exist", arg[iarg + 1], + dgroupid); propflag = 1; delete[] idprop; idprop = utils::strdup(arg[iarg + 1]); @@ -90,21 +106,37 @@ FixGroup::FixGroup(LAMMPS *lmp, int narg, char **arg) : if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "group dynamic every", error); nevery = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); if (nevery <= 0) - error->all(FLERR, "Illegal every value {} for dynamic group {}", nevery, dgroupid); + error->all(FLERR, ioffset + iarg + 1, "Illegal every value {} for dynamic group {}", nevery, + dgroupid); iarg += 2; } else if (strcmp(arg[iarg], "include") == 0) { if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "group dynamic include", error); if (strcmp(arg[iarg + 1], "molecule") == 0) { moleculeflag = 1; if (!atom->molecule_flag) - error->all(FLERR, + error->all(FLERR, ioffset + iarg + 1, "Dynamic Group include molecule setting requires atom attribute molecule"); iarg += 2; } else { - error->all(FLERR, "Unknown include setting {} in dynamic group command", arg[iarg + 1]); + error->all(FLERR, ioffset + iarg + 1, "Unknown include setting {} in dynamic group command", + arg[iarg + 1]); } + } else if (strcmp(arg[iarg], "within") == 0) { + if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "group dynamic within", error); + withinflag = 1; + cutoff = utils::numeric(FLERR, arg[iarg + 1], false, lmp); + if (cutoff <= 0.0) + error->all(FLERR, ioffset + iarg + 1, "Illegal within cutoff value {} for dynamic group {}", + cutoff, dgroupid); + iarg += 2; + } else if (strcmp(arg[iarg], "exclude") == 0) { + if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "group dynamic exclude", error); + excludeflag = 1; + delete[] idexclude; + idexclude = utils::strdup(arg[iarg + 1]); + iarg += 2; } else - error->all(FLERR, "Unknown keyword {} in dynamic group command", arg[iarg]); + error->all(FLERR, ioffset + iarg, "Unknown keyword {} in dynamic group command", arg[iarg]); } } @@ -115,6 +147,7 @@ FixGroup::~FixGroup() delete[] idregion; delete[] idvar; delete[] idprop; + delete[] idexclude; } /* ---------------------------------------------------------------------- */ @@ -127,6 +160,13 @@ int FixGroup::setmask() /* ---------------------------------------------------------------------- */ +void FixGroup::init_list(int /*id*/, NeighList *ptr) +{ + list = ptr; +} + +/* ---------------------------------------------------------------------- */ + void FixGroup::init() { std::string dyngroup = group->names[igroup]; @@ -162,6 +202,18 @@ void FixGroup::init() error->all(FLERR, "Custom per-atom property vector {} for dynamic group {} does not exist", idprop, dyngroup); } + + if (withinflag) { + NeighRequest *req = nullptr; + if (nevery == 1) { + req = neighbor->add_request(this, NeighConst::REQ_FULL); + } else { + req = neighbor->add_request(this, NeighConst::REQ_FULL | NeighConst::REQ_OCCASIONAL); + } + req->set_cutoff(cutoff); + } + + if (excludeflag) excludebit = group->get_bitmask_by_id(FLERR, idexclude, "group dynamic exclude"); } /* ---------------------------------------------------------------------- @@ -221,6 +273,10 @@ void FixGroup::set_group() if (regionflag) region->prematch(); + // re-build occasional neighbor list for "within" processing + + if (withinflag && (nevery != 1)) neighbor->build_one(list); + // set mask for each atom // only in group if in parent group, in region, variable is non-zero @@ -246,13 +302,93 @@ void FixGroup::set_group() mask[i] &= gbitinverse; } - if (varflag) memory->destroy(var); - // ensure ghost atom masks are also updated comm->forward_comm(this); - if (moleculeflag) group->add_molecules(0,gbit); + // select additional atoms that are within cutoff distance of already selected atoms + + if (withinflag) { + int i, j, ii, jj, inum, jnum; + const int *ilist, *jlist; + double dx, dy, dz, rsq; + + inum = list->inum; + ilist = list->ilist; + const int *const numneigh = list->numneigh; + const int *const *const firstneigh = list->firstneigh; + + const double cutsq = cutoff * cutoff; + for (ii = 0; ii < inum; ++ii) { + i = ilist[ii]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + for (jj = 0; jj < jnum; ++jj) { + j = jlist[jj] & NEIGHMASK; + + // check pairs where only one of the two atoms is in the parent group and passes the + // selection criteria from above. we cannot use mask[] & gbit for this for local atoms, + // since we add atoms to the group and we don't want to add neighbors of these added atoms + if ((mask[i] & groupbit) && !(mask[j] & groupbit)) { + inflag = 1; + if (regionflag && !region->match(x[i][0], x[i][1], x[i][2])) inflag = 0; + if (varflag && var[i] == 0.0) inflag = 0; + if (propflag) { + if (!proptype && ivector[i] == 0) inflag = 0; + if (proptype && dvector[i] == 0.0) inflag = 0; + } + if (inflag && (j < nlocal)) { + dx = x[j][0] - x[i][0]; + dy = x[j][1] - x[i][1]; + dz = x[j][2] - x[i][2]; + rsq = dx * dx + dy * dy + dz * dz; + if (rsq <= cutsq) mask[j] |= gbit; + } + } else if (!(mask[i] & groupbit) && (mask[j] & groupbit)) { + inflag = 1; + // for ghost atoms, we do not have access to all data, but we can use gbit instead + // since it has already been forward communicated and we will only select local atoms + if (j >= nlocal) { + if (!(mask[j] & gbit)) inflag = 0; + } else { + if (regionflag && !region->match(x[j][0], x[j][1], x[j][2])) inflag = 0; + if (varflag && var[j] == 0.0) inflag = 0; + if (propflag) { + if (!proptype && ivector[j] == 0) inflag = 0; + if (proptype && dvector[j] == 0.0) inflag = 0; + } + } + if (inflag) { + dx = x[j][0] - x[i][0]; + dy = x[j][1] - x[i][1]; + dz = x[j][2] - x[i][2]; + rsq = dx * dx + dy * dy + dz * dz; + if (rsq <= cutsq) mask[i] |= gbit; + } + } + } + } + + // we need a second forward communication, since we could only update the masks of local atoms + comm->forward_comm(this); + } + + // no longer needed + + if (varflag) memory->destroy(var); + + // add atoms that have the same molecule ID as selected atoms + if (moleculeflag) group->add_molecules(0, gbit); + + // exclude selected atoms in the excluded group + if (excludeflag) { + int nall = nlocal + atom->nghost; + for (int i = 0; i < nall; ++i) { + if (mask[i] & gbit) { + if (mask[i] & excludebit) mask[i] &= gbitinverse; + } + } + } } /* ---------------------------------------------------------------------- */ diff --git a/src/fix_group.h b/src/fix_group.h index 7e77cc18946..3b327317366 100644 --- a/src/fix_group.h +++ b/src/fix_group.h @@ -29,6 +29,7 @@ class FixGroup : public Fix { FixGroup(class LAMMPS *, int, char **); ~FixGroup() override; int setmask() override; + void init_list(int, class NeighList *) override; void init() override; void setup(int) override; void post_force(int) override; @@ -39,10 +40,13 @@ class FixGroup : public Fix { private: int gbit, gbitinverse; - int regionflag, varflag, propflag, proptype, moleculeflag; + int regionflag, varflag, propflag, proptype, moleculeflag, withinflag, excludeflag; int ivar, iprop; - char *idregion, *idvar, *idprop; + int excludebit; + double cutoff; + char *idregion, *idvar, *idprop, *idexclude; class Region *region; + class NeighList *list; int nlevels_respa; diff --git a/unittest/commands/test_groups.cpp b/unittest/commands/test_groups.cpp index 4694569dda7..148f143c23f 100644 --- a/unittest/commands/test_groups.cpp +++ b/unittest/commands/test_groups.cpp @@ -291,7 +291,7 @@ TEST_F(GroupTest, Bitmap) group->get_bitmask_by_id(FLERR, "five", "unittest 5");); } -TEST_F(GroupTest, Dynamic) +TEST_F(GroupTest, DynamicAtomic) { atomic_system(); @@ -339,6 +339,27 @@ TEST_F(GroupTest, Dynamic) END_HIDE_OUTPUT(); ASSERT_EQ(group->ngroup, 3); + BEGIN_HIDE_OUTPUT(); + command("region chunk block -1.1 1.1 -1.1 1.1 -1.1 0.1"); + command("group chunk dynamic all region chunk every 1"); + END_HIDE_OUTPUT(); + ASSERT_EQ(group->count(group->find("all")), 64); + ASSERT_EQ(group->count(group->find("chunk")), 0); + ASSERT_EQ(group->ngroup, 4); + BEGIN_HIDE_OUTPUT(); + command("run 10 post no"); + END_HIDE_OUTPUT(); + ASSERT_EQ(group->count(group->find("chunk")), 4); + BEGIN_HIDE_OUTPUT(); + command("group chunk delete"); + command("group chunk region chunk"); + command("group within dynamic chunk region chunk every 1 within 2.0"); + command("comm_modify cutoff 4.1"); + command("run 10 post no"); + END_HIDE_OUTPUT(); + ASSERT_EQ(group->count(group->find("chunk")), 4); + ASSERT_EQ(group->count(group->find("within")), 52); + TEST_FAILURE(".*ERROR: Group dynamic cannot reference itself.*", command("group half dynamic half region top");); TEST_FAILURE(".*ERROR: Group dynamic parent group dummy does not exist.*", @@ -350,6 +371,35 @@ TEST_F(GroupTest, Dynamic) command("group ramp variable grow");); } +TEST_F(GroupTest, DynamicMolecular) +{ + molecular_system(); + + BEGIN_HIDE_OUTPUT(); + command("region chunk block -1.1 1.1 -1.1 1.1 -1.1 0.1"); + command("group chunk dynamic all region chunk every 10 include molecule"); + END_HIDE_OUTPUT(); + EXPECT_EQ(group->count(group->find("all")), 64); + EXPECT_EQ(group->count(group->find("chunk")), 0); + ASSERT_EQ(group->ngroup, 2); + BEGIN_HIDE_OUTPUT(); + command("run 10 post no"); + END_HIDE_OUTPUT(); + EXPECT_EQ(group->count(group->find("chunk")), 8); + + BEGIN_HIDE_OUTPUT(); + command("group chunk delete"); + command("group chunk region chunk"); + command("group within dynamic chunk within 2.0 include molecule"); + command("group exclude dynamic chunk within 2.0 exclude chunk"); + command("comm_modify cutoff 4.1"); + command("run 10 post no"); + END_HIDE_OUTPUT(); + EXPECT_EQ(group->count(group->find("chunk")), 4); + EXPECT_EQ(group->count(group->find("within")), 59); + EXPECT_EQ(group->count(group->find("exclude")), 48); +} + static constexpr double EPSILON = 1.0e-13; TEST_F(GroupTest, VariableFunctions)