From a1011b606e0e07399705234b7fa04ee43b16a36a Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Mon, 6 Jul 2020 10:50:34 -0600 Subject: [PATCH 01/29] new reset_mol_ids command --- doc/src/Commands_all.rst | 1 + doc/src/commands_list.rst | 1 + doc/src/compute_cluster_atom.rst | 25 +++-- doc/src/reset_mol_ids.rst | 53 ++++++++++ src/compute_fragment_atom.cpp | 176 ++++++++++++++++++------------- src/compute_fragment_atom.h | 3 +- src/reset_mol_ids.cpp | 131 +++++++++++++++++++++++ src/reset_mol_ids.h | 60 +++++++++++ 8 files changed, 365 insertions(+), 85 deletions(-) create mode 100644 doc/src/reset_mol_ids.rst create mode 100644 src/reset_mol_ids.cpp create mode 100644 src/reset_mol_ids.h diff --git a/doc/src/Commands_all.rst b/doc/src/Commands_all.rst index 458238cca6..55ed872948 100644 --- a/doc/src/Commands_all.rst +++ b/doc/src/Commands_all.rst @@ -109,6 +109,7 @@ An alphabetic list of all general LAMMPS commands. * :doc:`replicate ` * :doc:`rerun ` * :doc:`reset_ids ` + * :doc:`reset_mol_ids ` * :doc:`reset_timestep ` * :doc:`restart ` * :doc:`run ` diff --git a/doc/src/commands_list.rst b/doc/src/commands_list.rst index 349d388923..49cf7ea588 100644 --- a/doc/src/commands_list.rst +++ b/doc/src/commands_list.rst @@ -89,6 +89,7 @@ Commands replicate rerun reset_ids + reset_mol_ids reset_timestep restart run diff --git a/doc/src/compute_cluster_atom.rst b/doc/src/compute_cluster_atom.rst index d9742a4a4b..502d2858d0 100644 --- a/doc/src/compute_cluster_atom.rst +++ b/doc/src/compute_cluster_atom.rst @@ -29,7 +29,6 @@ Examples compute 1 all cluster/atom 3.5 compute 1 all fragment/atom - compute 1 all aggregate/atom 3.5 Description @@ -43,13 +42,15 @@ cutoff distance from one or more other atoms in the cluster. If an atom has no neighbors within the cutoff distance, then it is a 1-atom cluster. -A fragment is similarly defined as a set of atoms, each of -which has an explicit bond (i.e. defined via a :doc:`data file `, -the :doc:`create_bonds ` command, or through fixes like -:doc:`fix bond/create `, :doc:`fix bond/swap `, -or :doc:`fix bond/break `). The cluster ID or fragment ID -of every atom in the cluster will be set to the smallest atom ID of any atom -in the cluster or fragment, respectively. +A fragment is similarly defined as a set of atoms, each of which has a +bond to another atom in the fragment, as defined initially via the +:doc:`data file ` or :doc:`create_bonds ` +commands, or by fixes which dynamically create or break bonds like +:doc:`fix bond/react `, :doc:`fix bond/create +`, :doc:`fix bond/swap `, or :doc:`fix +bond/break `. The cluster ID or fragment ID of every +atom in the cluster will be set to the smallest atom ID of any atom in +the cluster or fragment, respectively. An aggregate is defined by combining the rules for clusters and fragments, i.e. a set of atoms, where each of it is within the cutoff @@ -89,6 +90,14 @@ multiple compute/dump commands, each of a *cluster/atom* or :doc:`special_bonds ` command that includes all pairs in the neighbor list. +.. note:: + + For the compute fragment/atom style, each fragment is identified + using the current bond topology within each fragment. This will + not account for bonds broken by the :doc:`bond_style quartic + ` command because it does not perform a full update + of the bond topology data structures within LAMMPS. + **Output info:** This compute calculates a per-atom vector, which can be accessed by diff --git a/doc/src/reset_mol_ids.rst b/doc/src/reset_mol_ids.rst new file mode 100644 index 0000000000..085f912307 --- /dev/null +++ b/doc/src/reset_mol_ids.rst @@ -0,0 +1,53 @@ +.. index:: reset_mol_ids + +reset_mol_ids command +================= + +Syntax +"""""" + +Syntax +"""""" + +.. parsed-literal:: + + reset_mol_ids group-ID + +* group-ID = ID of group of atoms whose molecule IDs will be reset + +Examples +"""""""" + +.. code-block:: LAMMPS + + reset_mol_ids all + reset_mol_ids solvent + +Description +""""""""""" + +Reset molecule IDs for a group of atoms. This will create a new set +of molecule IDs that are numbered contiguously from 1 to N, if there +are N different molecule IDs used by the group. Only molecule IDs for +atoms in the specified group are reset. + +This can be useful to invoke after performing a reactive molecular +dynamics run with :doc:`fix bond/react `, :doc:`fix +bond/create `, or :doc:`fix bond/delete +`. It can also be useful after any simulation which +has lost molecules, e.g. via the :doc:`fix evaporate ` +command. + +Restrictions +"""""""""""" +none + +Related commands +"""""""""""""""" + +:doc:`reset_ids `, :doc:`fix bond/react `, +:doc:`fix bond/create `, +:doc:`fix bond/delete `, +:doc:`fix evaporate ` + +**Default:** none diff --git a/src/compute_fragment_atom.cpp b/src/compute_fragment_atom.cpp index f9b68dd217..6239ea9101 100644 --- a/src/compute_fragment_atom.cpp +++ b/src/compute_fragment_atom.cpp @@ -23,14 +23,15 @@ #include "update.h" #include "modify.h" #include "force.h" +#include "group.h" #include "comm.h" #include "memory.h" #include "error.h" -#include "group.h" - using namespace LAMMPS_NS; +#define BIG 1.0e20 + /* ---------------------------------------------------------------------- */ ComputeFragmentAtom::ComputeFragmentAtom(LAMMPS *lmp, int narg, char **arg) : @@ -45,15 +46,20 @@ ComputeFragmentAtom::ComputeFragmentAtom(LAMMPS *lmp, int narg, char **arg) : peratom_flag = 1; size_peratom_cols = 0; comm_forward = 1; - comm_reverse = 1; nmax = 0; + stack = NULL; + clist = NULL; + markflag = NULL; } /* ---------------------------------------------------------------------- */ ComputeFragmentAtom::~ComputeFragmentAtom() { + memory->destroy(stack); + memory->destroy(clist); + memory->destroy(markflag); memory->destroy(fragmentID); } @@ -63,8 +69,8 @@ void ComputeFragmentAtom::init() { if (atom->tag_enable == 0) error->all(FLERR,"Cannot use compute fragment/atom unless atoms have IDs"); - if (force->bond == NULL) - error->all(FLERR,"Compute fragment/atom requires a bond style to be defined"); + if (!atom->molecular) + error->all(FLERR,"Compute fragment/atom requires a molecular system"); int count = 0; for (int i = 0; i < modify->ncompute; i++) @@ -77,15 +83,24 @@ void ComputeFragmentAtom::init() void ComputeFragmentAtom::compute_peratom() { - int i,j,k; + int i,j,k,m,n; + int nstack,ncluster,done,alldone; + double newID,cID; + tagint *list; invoked_peratom = update->ntimestep; - // grow fragmentID array if necessary + // grow work and fragmentID vectors if necessary if (atom->nmax > nmax) { + memory->destroy(stack); + memory->destroy(clist); + memory->destroy(markflag); memory->destroy(fragmentID); nmax = atom->nmax; + memory->create(stack,nmax,"fragment/atom:stack"); + memory->create(clist,nmax,"fragment/atom:clist"); + memory->create(markflag,nmax,"fragment/atom:markflag"); memory->create(fragmentID,nmax,"fragment/atom:fragmentID"); vector_atom = fragmentID; } @@ -97,63 +112,104 @@ void ComputeFragmentAtom::compute_peratom() comm->forward_comm_compute(this); } - // each atom starts in its own fragment, + // owned + ghost atoms start with fragmentID = atomID + // atoms not in group have fragmentID = 0 - int nlocal = atom->nlocal; tagint *tag = atom->tag; int *mask = atom->mask; - int *num_bond = atom->num_bond; - int **bond_type = atom->bond_type; - tagint **bond_atom = atom->bond_atom; - - for (i = 0; i < nlocal + atom->nghost; i++) + tagint **special = atom->special; + int **nspecial = atom->nspecial; + int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; + + for (i = 0; i < nall; i++) if (mask[i] & groupbit) fragmentID[i] = tag[i]; else fragmentID[i] = 0; - // loop until no more changes on any proc: + // loop until no ghost atom fragment ID is changed // acquire fragmentIDs of ghost atoms - // loop over my atoms, and check atoms bound to it - // if both atoms are in fragment, assign lowest fragmentID to both - // iterate until no changes in my atoms - // then check if any proc made changes + // loop over clusters of atoms, which include ghost atoms + // set fragmentIDs for each cluster to min framentID in the clusters + // iterate until no changes to ghost atom fragmentIDs commflag = 1; - int change,done,anychange; - + int iteration = 0; + while (1) { + iteration++; + comm->forward_comm_compute(this); + done = 1; - // reverse communication when bonds are not stored on every processor + // set markflag = 0 for all owned atoms, for new iteration - if (force->newton_bond) - comm->reverse_comm_compute(this); + for (i = 0; i < nlocal; i++) markflag[i] = 0; - change = 0; - while (1) { - done = 1; - for (i = 0; i < nlocal; i++) { - if (!(mask[i] & groupbit)) continue; + // loop over all owned atoms + // each unmarked atom starts a cluster search + + for (i = 0; i < nlocal; i++) { - for (j = 0; j < num_bond[i]; j++) { - if (bond_type[i][j] == 0) continue; - k = atom->map(bond_atom[i][j]); - if (k < 0) continue; - if (!(mask[k] & groupbit)) continue; - if (fragmentID[i] == fragmentID[k]) continue; + // skip atom I if not in group or already marked + + if (!(mask[i] & groupbit)) continue; + if (markflag[i]) continue; - fragmentID[i] = fragmentID[k] = MIN(fragmentID[i],fragmentID[k]); - done = 0; - } + // find one cluster of bond-connected atoms + // ncluster = # of owned and ghost atoms in cluster + // clist = vector of local indices of the ncluster atoms + // stack is used to walk the bond topology + + ncluster = nstack = 0; + stack[nstack++] = i; + + while (nstack) { + j = stack[--nstack]; + clist[ncluster++] = j; + markflag[j] = 1; + + n = nspecial[j][0]; + list = special[j]; + for (m = 0; m < n; m++) { + k = atom->map(list[m]); + + // skip bond neighbor K if not in group or already marked + + if (k < 0) continue; + if (!(mask[k] & groupbit)) continue; + if (k < nlocal && markflag[k]) continue; + + // owned bond neighbors are added to stack for further walking + // ghost bond neighbors are added directly w/out use of stack + + if (k < nlocal) stack[nstack++] = k; + else clist[ncluster++] = k; + } } - if (!done) change = 1; - if (done) break; - } + // newID = minimum fragment ID in cluster list, including ghost atoms + + newID = BIG; + for (m = 0; m < ncluster; m++) { + cID = fragmentID[clist[m]]; + newID = MIN(newID,cID); + } + + // set fragmentID = newID for all atoms in cluster, including ghost atoms + // not done with iterations if change the fragmentID of a ghost atom + + for (m = 0; m < ncluster; m++) { + j = clist[m]; + if (j >= nlocal && fragmentID[j] != newID) done = 0; + fragmentID[j] = newID; + } + } + // stop if all procs are done - MPI_Allreduce(&change,&anychange,1,MPI_INT,MPI_MAX,world); - if (!anychange) break; + MPI_Allreduce(&done,&alldone,1,MPI_INT,MPI_MIN,world); + if (alldone) break; } } @@ -203,43 +259,13 @@ void ComputeFragmentAtom::unpack_forward_comm(int n, int first, double *buf) } } -/* ---------------------------------------------------------------------- */ - -int ComputeFragmentAtom::pack_reverse_comm(int n, int first, double *buf) -{ - int i,m,last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - buf[m++] = fragmentID[i]; - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -void ComputeFragmentAtom::unpack_reverse_comm(int n, int *list, double *buf) -{ - int i,j,m; - - m = 0; - for (i = 0; i < n; i++) { - j = list[i]; - double x = buf[m++]; - - // only overwrite local IDs with values lower than current ones - - fragmentID[j] = MIN(x,fragmentID[j]); - } -} - /* ---------------------------------------------------------------------- - memory usage of local atom-based array + memory usage of local atom-based arrays ------------------------------------------------------------------------- */ double ComputeFragmentAtom::memory_usage() { double bytes = nmax * sizeof(double); + bytes += 3*nmax * sizeof(int); return bytes; } diff --git a/src/compute_fragment_atom.h b/src/compute_fragment_atom.h index 6031d43841..1422a71fb4 100644 --- a/src/compute_fragment_atom.h +++ b/src/compute_fragment_atom.h @@ -32,12 +32,11 @@ class ComputeFragmentAtom : public Compute { void compute_peratom(); int pack_forward_comm(int, int *, double *, int, int *); void unpack_forward_comm(int, int, double *); - int pack_reverse_comm(int, int, double *); - void unpack_reverse_comm(int, int *, double *); double memory_usage(); private: int nmax,commflag; + int *stack,*clist,*markflag; double *fragmentID; }; diff --git a/src/reset_mol_ids.cpp b/src/reset_mol_ids.cpp new file mode 100644 index 0000000000..0b1b705279 --- /dev/null +++ b/src/reset_mol_ids.cpp @@ -0,0 +1,131 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Jacob Gissinger (jacob.gissinger@colorado.edu) +------------------------------------------------------------------------- */ + +#include "reset_mol_ids.h" +#include +#include +#include "atom.h" +#include "domain.h" +#include "comm.h" +#include "group.h" +#include "modify.h" +#include "compute_fragment_atom.h" +#include "compute_chunk_atom.h" +#include "utils.h" +#include "error.h" +#include "fmt/format.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +ResetMolIDs::ResetMolIDs(LAMMPS *lmp) : Pointers(lmp) {} + +/* ---------------------------------------------------------------------- */ + +void ResetMolIDs::command(int narg, char **arg) +{ + if (domain->box_exist == 0) + error->all(FLERR,"Reset_mol_ids command before simulation box is defined"); + if (atom->tag_enable == 0) + error->all(FLERR,"Cannot use reset_mol_ids unless atoms have IDs"); + if (atom->molecular != 1) + error->all(FLERR,"Can only use reset_mol_ids on molecular systems"); + + if (narg != 1) error->all(FLERR,"Illegal reset_mol_ids command"); + int igroup = group->find(arg[0]); + if (igroup == -1) error->all(FLERR,"Could not find reset_mol_ids group ID"); + + if (comm->me == 0) { + if (screen) fprintf(screen,"Resetting molecule IDs ...\n"); + if (logfile) fprintf(logfile,"Resetting molecule IDs ...\n"); + } + + // record wall time for resetting molecule IDs + + MPI_Barrier(world); + double time1 = MPI_Wtime(); + + // create instances of compute fragment/atom and compute chunk/atom + + std::string idfrag = "reset_mol_ids_FRAGMENT_ATOM"; + std::string fragcmd = idfrag + fmt::format(" {} fragment/atom",arg[0]); + modify->add_compute(fragcmd); + + std::string idchunk = "reset_mol_ids_CHUNK_ATOM"; + std::string chunkcmd = idchunk + fmt::format(" all chunk/atom molecule compress yes"); + modify->add_compute(chunkcmd); + + // initialize system since comm->borders() will be invoked + + lmp->init(); + + // setup domain, communication + // exchange will clear map, borders will reset + // this is the map needed to lookup current global IDs for bond topology + + if (domain->triclinic) domain->x2lamda(atom->nlocal); + domain->pbc(); + domain->reset_box(); + comm->setup(); + comm->exchange(); + comm->borders(); + if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost); + + // invoke peratom method of compute fragment/atom + // walks bond connectivity and assigns each atom a fragment ID + + int icompute = modify->find_compute(idfrag); + ComputeFragmentAtom *cfa = (ComputeFragmentAtom *) modify->compute[icompute]; + cfa->compute_peratom(); + double *ids = cfa->vector_atom; + + // copy fragmentID to molecule ID + + tagint *molecule = atom->molecule; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + molecule[i] = static_cast (ids[i]); + + // invoke peratom method of compute chunk/atom + // compresses new molecule IDs to be contiguous 1 to Nmol + + icompute = modify->find_compute(idchunk); + ComputeChunkAtom *cca = (ComputeChunkAtom *) modify->compute[icompute]; + cca->compute_peratom(); + ids = cca->vector_atom; + int nchunk = cca->nchunk; + + // copy chunkID to molecule ID + + for (int i = 0; i < nlocal; i++) + molecule[i] = static_cast (ids[i]); + + // clean up + + modify->delete_compute(idchunk); + modify->delete_compute(idfrag); + + // total time + + MPI_Barrier(world); + + if (comm->me == 0) + utils::logmesg(lmp,fmt::format(" reset_mol_ids CPU = {:.3f} seconds\n", + MPI_Wtime()-time1)); +} diff --git a/src/reset_mol_ids.h b/src/reset_mol_ids.h new file mode 100644 index 0000000000..83c0b5d80f --- /dev/null +++ b/src/reset_mol_ids.h @@ -0,0 +1,60 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef COMMAND_CLASS + +CommandStyle(reset_mol_ids,ResetMolIDs) + +#else + +#ifndef LMP_RESET_MOL_IDS_H +#define LMP_RESET_MOL_IDS_H + +#include "pointers.h" + +namespace LAMMPS_NS { + +class ResetMolIDs : protected Pointers { + public: + ResetMolIDs(class LAMMPS *); + void command(int, char **); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Reset_mol_ids command before simulation box is defined + +UNDOCUMENTED + +E: Can only use reset_mol_ids on molecular systems + +UNDOCUMENTED + +E: Illegal ... command + +UNDOCUMENTED + +E: Cannot use reset_mol_ids unless molecules have IDs + +UNDOCUMENTED + +E: Reset_mol_ids missing %d bond topology atom IDs - use comm_modify cutoff + +UNDOCUMENTED + +*/ From 836570ec26838061bf02aba36304bec5b7a624cc Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 6 Jul 2020 14:12:37 -0400 Subject: [PATCH 02/29] update docs --- doc/src/reset_mol_ids.rst | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/doc/src/reset_mol_ids.rst b/doc/src/reset_mol_ids.rst index 085f912307..8d2d6eeab5 100644 --- a/doc/src/reset_mol_ids.rst +++ b/doc/src/reset_mol_ids.rst @@ -1,7 +1,7 @@ .. index:: reset_mol_ids reset_mol_ids command -================= +===================== Syntax """""" @@ -33,10 +33,11 @@ atoms in the specified group are reset. This can be useful to invoke after performing a reactive molecular dynamics run with :doc:`fix bond/react `, :doc:`fix -bond/create `, or :doc:`fix bond/delete -`. It can also be useful after any simulation which -has lost molecules, e.g. via the :doc:`fix evaporate ` -command. +bond/create `, or :doc:`fix bond/break +`. It can also be useful after molecules have been +deleted with :doc:`delete_atoms ` or after a simulation +which has lost molecules, e.g. via the :doc:`fix evaporate +` command. Restrictions """""""""""" @@ -47,7 +48,8 @@ Related commands :doc:`reset_ids `, :doc:`fix bond/react `, :doc:`fix bond/create `, -:doc:`fix bond/delete `, -:doc:`fix evaporate ` +:doc:`fix bond/break `, +:doc:`fix evaporate `, +:doc:`delete_atoms ` **Default:** none From d37e943e8d5ed2eb9a9ed502b73168a2f807d661 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 6 Jul 2020 14:12:50 -0400 Subject: [PATCH 03/29] refactor tester --- unittest/formats/test_potential_file_reader.cpp | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/unittest/formats/test_potential_file_reader.cpp b/unittest/formats/test_potential_file_reader.cpp index 0533d547c2..8a6db0197c 100644 --- a/unittest/formats/test_potential_file_reader.cpp +++ b/unittest/formats/test_potential_file_reader.cpp @@ -39,6 +39,13 @@ using namespace LAMMPS_NS; using ::testing::MatchesRegex; using utils::split_words; +#define TEST_FAILURE(...) \ + if (Info::has_exceptions()) { \ + ASSERT_ANY_THROW({__VA_ARGS__}); \ + } else { \ + ASSERT_DEATH({__VA_ARGS__}, ""); \ + } + // whether to print verbose output (i.e. not capturing LAMMPS screen output). bool verbose = false; @@ -108,15 +115,9 @@ TEST_F(PotentialFileReaderTest, Sw_noconv) if (!verbose) ::testing::internal::CaptureStdout(); lmp->input->one("units real"); if (!verbose) ::testing::internal::GetCapturedStdout(); + ::testing::internal::CaptureStdout(); - if (Info::has_exceptions()) { - ASSERT_ANY_THROW( - PotentialFileReader reader(lmp, "Si.sw", "Stillinger-Weber", utils::REAL2METAL);); - } else { - ASSERT_DEATH( - { PotentialFileReader reader(lmp, "Si.sw", "Stillinger-Weber", utils::REAL2METAL); }, - ""); - } + TEST_FAILURE(PotentialFileReader reader(lmp, "Si.sw", "Stillinger-Weber", utils::REAL2METAL);); std::string mesg = ::testing::internal::GetCapturedStdout(); ASSERT_THAT(mesg, MatchesRegex(".*ERROR on proc.*potential.*requires metal units but real.*")); } From 89f0116eab777938260f91fea56ee8c67b642742 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 6 Jul 2020 21:12:55 -0400 Subject: [PATCH 04/29] fix communication data conversion bug corrupting bond list --- src/reset_ids.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/reset_ids.cpp b/src/reset_ids.cpp index 671ba7940f..6a7a8e8b2a 100644 --- a/src/reset_ids.cpp +++ b/src/reset_ids.cpp @@ -155,7 +155,7 @@ void ResetIDs::command(int narg, char **arg) for (j = 0; j < num_bond[i]; j++) { oldID = bond_atom[i][j]; m = atom->map(oldID); - if (m >= 0) bond_atom[i][j] = static_cast (newIDs[m][0]); + if (m >= 0) bond_atom[i][j] = (tagint) ubuf(newIDs[m][0]).i; else badcount++; } } From 2351f99bef18f6ce7ddb7f635c2839d9b21c6a86 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 6 Jul 2020 21:26:11 -0400 Subject: [PATCH 05/29] add (incomplete) tester for reset_mol_ids (with a little use of reset_ids, too) --- unittest/commands/CMakeLists.txt | 5 + unittest/commands/data.fourmol | 1 + unittest/commands/in.fourmol | 1 + unittest/commands/test_reset_mol_ids.cpp | 250 +++++++++++++++++++++++ 4 files changed, 257 insertions(+) create mode 120000 unittest/commands/data.fourmol create mode 120000 unittest/commands/in.fourmol create mode 100644 unittest/commands/test_reset_mol_ids.cpp diff --git a/unittest/commands/CMakeLists.txt b/unittest/commands/CMakeLists.txt index 626301d404..d0793b707f 100644 --- a/unittest/commands/CMakeLists.txt +++ b/unittest/commands/CMakeLists.txt @@ -2,3 +2,8 @@ add_executable(test_simple_commands test_simple_commands.cpp) target_link_libraries(test_simple_commands PRIVATE lammps GTest::GMock GTest::GTest) add_test(NAME SimpleCommands COMMAND test_simple_commands WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) + +add_executable(test_reset_mol_ids test_reset_mol_ids.cpp) +target_compile_definitions(test_reset_mol_ids PRIVATE -DTEST_INPUT_FOLDER=${CMAKE_CURRENT_SOURCE_DIR}) +target_link_libraries(test_reset_mol_ids PRIVATE lammps GTest::GMock GTest::GTest) +add_test(NAME ResetMolIDs COMMAND test_reset_mol_ids WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) diff --git a/unittest/commands/data.fourmol b/unittest/commands/data.fourmol new file mode 120000 index 0000000000..ebbe325d24 --- /dev/null +++ b/unittest/commands/data.fourmol @@ -0,0 +1 @@ +../force-styles/tests/data.fourmol \ No newline at end of file diff --git a/unittest/commands/in.fourmol b/unittest/commands/in.fourmol new file mode 120000 index 0000000000..cc47b5adb9 --- /dev/null +++ b/unittest/commands/in.fourmol @@ -0,0 +1 @@ +../force-styles/tests/in.fourmol \ No newline at end of file diff --git a/unittest/commands/test_reset_mol_ids.cpp b/unittest/commands/test_reset_mol_ids.cpp new file mode 100644 index 0000000000..bc38c6404c --- /dev/null +++ b/unittest/commands/test_reset_mol_ids.cpp @@ -0,0 +1,250 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "atom.h" +#include "fmt/format.h" +#include "info.h" +#include "input.h" +#include "lammps.h" +#include "output.h" +#include "update.h" +#include "utils.h" +#include "gmock/gmock.h" +#include "gtest/gtest.h" + +#include +#include + +// whether to print verbose output (i.e. not capturing LAMMPS screen output). +bool verbose = false; + +using LAMMPS_NS::utils::split_words; + +namespace LAMMPS_NS { +using ::testing::MatchesRegex; + +#define GETIDX(i) lmp->atom->map(i) + +#define TEST_FAILURE(...) \ + if (Info::has_exceptions()) { \ + ASSERT_ANY_THROW({__VA_ARGS__}); \ + } else { \ + ASSERT_DEATH({__VA_ARGS__}, ""); \ + } + +#define STRINGIFY(val) XSTR(val) +#define XSTR(val) #val + +class ResetMolIDsTest : public ::testing::Test { +protected: + LAMMPS *lmp; + + void SetUp() override + { + const char *args[] = {"ResetMolIDsTest", "-log", "none", "-nocite", "-echo", "screen"}; + char **argv = (char **)args; + int argc = sizeof(args) / sizeof(char *); + if (!verbose) ::testing::internal::CaptureStdout(); + lmp = new LAMMPS(argc, argv, MPI_COMM_WORLD); + Info *info = new Info(lmp); + if (info->has_style("atom","full")) { + lmp->input->one("variable input_dir index " STRINGIFY(TEST_INPUT_FOLDER)); + lmp->input->one("include ${input_dir}/in.fourmol"); + } + if (!verbose) ::testing::internal::GetCapturedStdout(); + } + + void TearDown() override + { + if (!verbose) ::testing::internal::CaptureStdout(); + delete lmp; + if (!verbose) ::testing::internal::GetCapturedStdout(); + } +}; + +TEST_F(ResetMolIDsTest, Plain) +{ + if (lmp->atom->natoms == 0) GTEST_SKIP(); + + auto molid = lmp->atom->molecule; + ASSERT_EQ(molid[GETIDX(1)], 1); + ASSERT_EQ(molid[GETIDX(2)], 1); + ASSERT_EQ(molid[GETIDX(3)], 1); + ASSERT_EQ(molid[GETIDX(4)], 1); + ASSERT_EQ(molid[GETIDX(5)], 1); + ASSERT_EQ(molid[GETIDX(6)], 1); + ASSERT_EQ(molid[GETIDX(7)], 1); + ASSERT_EQ(molid[GETIDX(8)], 2); + ASSERT_EQ(molid[GETIDX(9)], 2); + ASSERT_EQ(molid[GETIDX(10)], 2); + ASSERT_EQ(molid[GETIDX(11)], 2); + ASSERT_EQ(molid[GETIDX(12)], 2); + ASSERT_EQ(molid[GETIDX(13)], 2); + ASSERT_EQ(molid[GETIDX(14)], 2); + ASSERT_EQ(molid[GETIDX(15)], 2); + ASSERT_EQ(molid[GETIDX(16)], 2); + ASSERT_EQ(molid[GETIDX(17)], 2); + ASSERT_EQ(molid[GETIDX(18)], 3); + ASSERT_EQ(molid[GETIDX(19)], 3); + ASSERT_EQ(molid[GETIDX(20)], 3); + ASSERT_EQ(molid[GETIDX(21)], 4); + ASSERT_EQ(molid[GETIDX(22)], 4); + ASSERT_EQ(molid[GETIDX(23)], 4); + ASSERT_EQ(molid[GETIDX(24)], 5); + ASSERT_EQ(molid[GETIDX(25)], 5); + ASSERT_EQ(molid[GETIDX(26)], 5); + ASSERT_EQ(molid[GETIDX(27)], 6); + ASSERT_EQ(molid[GETIDX(28)], 6); + ASSERT_EQ(molid[GETIDX(29)], 6); + + // the original data file has two different molecule IDs + // for two residues of the same molecule/fragment. + if (!verbose) ::testing::internal::CaptureStdout(); + lmp->input->one("reset_mol_ids all"); + if (!verbose) ::testing::internal::GetCapturedStdout(); + + ASSERT_EQ(molid[GETIDX(1)], 1); + ASSERT_EQ(molid[GETIDX(2)], 1); + ASSERT_EQ(molid[GETIDX(3)], 1); + ASSERT_EQ(molid[GETIDX(4)], 1); + ASSERT_EQ(molid[GETIDX(5)], 1); + ASSERT_EQ(molid[GETIDX(6)], 1); + ASSERT_EQ(molid[GETIDX(7)], 1); + ASSERT_EQ(molid[GETIDX(8)], 1); + ASSERT_EQ(molid[GETIDX(9)], 1); + ASSERT_EQ(molid[GETIDX(10)], 1); + ASSERT_EQ(molid[GETIDX(11)], 1); + ASSERT_EQ(molid[GETIDX(12)], 1); + ASSERT_EQ(molid[GETIDX(13)], 1); + ASSERT_EQ(molid[GETIDX(14)], 1); + ASSERT_EQ(molid[GETIDX(15)], 1); + ASSERT_EQ(molid[GETIDX(16)], 1); + ASSERT_EQ(molid[GETIDX(17)], 1); + ASSERT_EQ(molid[GETIDX(18)], 2); + ASSERT_EQ(molid[GETIDX(19)], 2); + ASSERT_EQ(molid[GETIDX(20)], 2); + ASSERT_EQ(molid[GETIDX(21)], 3); + ASSERT_EQ(molid[GETIDX(22)], 3); + ASSERT_EQ(molid[GETIDX(23)], 3); + ASSERT_EQ(molid[GETIDX(24)], 4); + ASSERT_EQ(molid[GETIDX(25)], 4); + ASSERT_EQ(molid[GETIDX(26)], 4); + ASSERT_EQ(molid[GETIDX(27)], 5); + ASSERT_EQ(molid[GETIDX(28)], 5); + ASSERT_EQ(molid[GETIDX(29)], 5); +} + +TEST_F(ResetMolIDsTest, DeletePlusAtomID) +{ + if (lmp->atom->natoms == 0) GTEST_SKIP(); + + auto molid = lmp->atom->molecule; + + // delete two water molecules + if (!verbose) ::testing::internal::CaptureStdout(); + lmp->input->one("group allwater molecule 3:6"); + lmp->input->one("group twowater molecule 4:6:2"); + lmp->input->one("delete_atoms group twowater compress no bond yes"); + lmp->input->one("reset_mol_ids all"); + if (!verbose) ::testing::internal::GetCapturedStdout(); + ASSERT_EQ(lmp->atom->natoms, 23); + ASSERT_EQ(lmp->atom->map_tag_max, 26); + + EXPECT_EQ(molid[GETIDX(1)], 1); + EXPECT_EQ(molid[GETIDX(2)], 1); + EXPECT_EQ(molid[GETIDX(3)], 1); + EXPECT_EQ(molid[GETIDX(4)], 1); + EXPECT_EQ(molid[GETIDX(5)], 1); + EXPECT_EQ(molid[GETIDX(6)], 1); + EXPECT_EQ(molid[GETIDX(7)], 1); + EXPECT_EQ(molid[GETIDX(8)], 1); + EXPECT_EQ(molid[GETIDX(9)], 1); + EXPECT_EQ(molid[GETIDX(10)], 1); + EXPECT_EQ(molid[GETIDX(11)], 1); + EXPECT_EQ(molid[GETIDX(12)], 1); + EXPECT_EQ(molid[GETIDX(13)], 1); + EXPECT_EQ(molid[GETIDX(14)], 1); + EXPECT_EQ(molid[GETIDX(15)], 1); + EXPECT_EQ(molid[GETIDX(16)], 1); + EXPECT_EQ(molid[GETIDX(17)], 1); + EXPECT_EQ(molid[GETIDX(18)], 2); + EXPECT_EQ(molid[GETIDX(19)], 2); + EXPECT_EQ(molid[GETIDX(20)], 2); + EXPECT_EQ(molid[GETIDX(24)], 3); + EXPECT_EQ(molid[GETIDX(25)], 3); + EXPECT_EQ(molid[GETIDX(26)], 3); + + // now also check and reset the atom ids + + EXPECT_GE(GETIDX(1), 0); + EXPECT_GE(GETIDX(2), 0); + EXPECT_GE(GETIDX(3), 0); + EXPECT_GE(GETIDX(4), 0); + EXPECT_GE(GETIDX(5), 0); + EXPECT_GE(GETIDX(6), 0); + EXPECT_GE(GETIDX(7), 0); + EXPECT_GE(GETIDX(8), 0); + EXPECT_GE(GETIDX(9), 0); + EXPECT_GE(GETIDX(10), 0); + EXPECT_GE(GETIDX(11), 0); + EXPECT_GE(GETIDX(12), 0); + EXPECT_GE(GETIDX(13), 0); + EXPECT_GE(GETIDX(14), 0); + EXPECT_GE(GETIDX(15), 0); + EXPECT_GE(GETIDX(16), 0); + EXPECT_GE(GETIDX(17), 0); + EXPECT_GE(GETIDX(18), 0); + EXPECT_GE(GETIDX(19), 0); + EXPECT_GE(GETIDX(20), 0); + EXPECT_EQ(GETIDX(21), -1); + EXPECT_EQ(GETIDX(22), -1); + EXPECT_EQ(GETIDX(23), -1); + EXPECT_GE(GETIDX(24), 0); + EXPECT_GE(GETIDX(25), 0); + EXPECT_GE(GETIDX(26), 0); + + if (!verbose) ::testing::internal::CaptureStdout(); + lmp->input->one("reset_ids"); + if (!verbose) ::testing::internal::GetCapturedStdout(); + + ASSERT_EQ(lmp->atom->map_tag_max, 23); + for (int i = 1; i <= 23; ++i) + EXPECT_GE(GETIDX(i), 0); +} + +// TEST_FAILURE(lmp->input->one("reset_mol_ids all");); +// ASSERT_THAT(mesg, MatchesRegex(".*ERROR: Unknown command.*")); + +} // namespace LAMMPS_NS + +int main(int argc, char **argv) +{ + MPI_Init(&argc, &argv); + ::testing::InitGoogleMock(&argc, argv); + + // handle arguments passed via environment variable + if (const char *var = getenv("TEST_ARGS")) { + std::vector env = split_words(var); + for (auto arg : env) { + if (arg == "-v") { + verbose = true; + } + } + } + + if ((argc > 1) && (strcmp(argv[1], "-v") == 0)) verbose = true; + + int rv = RUN_ALL_TESTS(); + MPI_Finalize(); + return rv; +} From 960addcc2cc9e385187c5091331eb3a761d396b3 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 6 Jul 2020 21:26:49 -0400 Subject: [PATCH 06/29] simplify code string/fmtlib code some more --- src/reset_mol_ids.cpp | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/src/reset_mol_ids.cpp b/src/reset_mol_ids.cpp index 0b1b705279..81ecf12a42 100644 --- a/src/reset_mol_ids.cpp +++ b/src/reset_mol_ids.cpp @@ -50,10 +50,8 @@ void ResetMolIDs::command(int narg, char **arg) int igroup = group->find(arg[0]); if (igroup == -1) error->all(FLERR,"Could not find reset_mol_ids group ID"); - if (comm->me == 0) { - if (screen) fprintf(screen,"Resetting molecule IDs ...\n"); - if (logfile) fprintf(logfile,"Resetting molecule IDs ...\n"); - } + if (comm->me == 0) + utils::logmesg(lmp,"Resetting molecule IDs ...\n"); // record wall time for resetting molecule IDs @@ -62,13 +60,12 @@ void ResetMolIDs::command(int narg, char **arg) // create instances of compute fragment/atom and compute chunk/atom - std::string idfrag = "reset_mol_ids_FRAGMENT_ATOM"; - std::string fragcmd = idfrag + fmt::format(" {} fragment/atom",arg[0]); - modify->add_compute(fragcmd); + const std::string idfrag = "reset_mol_ids_FRAGMENT_ATOM"; + modify->add_compute(fmt::format("{} {} fragment/atom",idfrag, arg[0])); - std::string idchunk = "reset_mol_ids_CHUNK_ATOM"; - std::string chunkcmd = idchunk + fmt::format(" all chunk/atom molecule compress yes"); - modify->add_compute(chunkcmd); + const std::string idchunk = "reset_mol_ids_CHUNK_ATOM"; + modify->add_compute(fmt::format("{} all chunk/atom molecule compress yes", + idchunk)); // initialize system since comm->borders() will be invoked From 53d20c9ebc9903642aa973369108880ed08210ae Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 6 Jul 2020 21:44:39 -0400 Subject: [PATCH 07/29] add missing symlink --- examples/threebody/CdTeZnSeHgS0.sw | 1 + 1 file changed, 1 insertion(+) create mode 120000 examples/threebody/CdTeZnSeHgS0.sw diff --git a/examples/threebody/CdTeZnSeHgS0.sw b/examples/threebody/CdTeZnSeHgS0.sw new file mode 120000 index 0000000000..405221bd72 --- /dev/null +++ b/examples/threebody/CdTeZnSeHgS0.sw @@ -0,0 +1 @@ +../../potentials/CdTeZnSeHgS0.sw \ No newline at end of file From 0b1443ed23032d940a18eb152dfc24bb418c3d41 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 7 Jul 2020 09:51:02 -0400 Subject: [PATCH 08/29] add prototypes for exception handling functions --- python/lammps.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/python/lammps.py b/python/lammps.py index ffcdbac456..feda10802e 100644 --- a/python/lammps.py +++ b/python/lammps.py @@ -204,6 +204,12 @@ class lammps(object): self.lib.lammps_neighlist_element_neighbors.argtypes = [c_void_p, c_int, c_int, POINTER(c_int), POINTER(c_int), POINTER(POINTER(c_int))] self.lib.lammps_neighlist_element_neighbors.restype = None + self.lib.lammps_has_error.argtypes = [v_void_p] + self.lib.lammps_has_error.restype = c_bool + + self.lib.lammps_get_last_error_message.argtypes = [c_void_p, c_char_p, c_int] + self.lib.lammps_get_last_error_message.restype = c_int + # if no ptr provided, create an instance of LAMMPS # don't know how to pass an MPI communicator from PyPar # but we can pass an MPI communicator from mpi4py v2.0.0 and later From bb9ab025c1c0965ab0195062660d3ec8011a6b09 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 7 Jul 2020 09:51:39 -0400 Subject: [PATCH 09/29] avoid exception when having a command fail due to an empty string --- src/input.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/input.cpp b/src/input.cpp index 27f8db7012..4a0f880505 100644 --- a/src/input.cpp +++ b/src/input.cpp @@ -227,7 +227,7 @@ void Input::file() // execute the command - if (execute_command()) + if (execute_command() && line) error->all(FLERR,fmt::format("Unknown command: {}",line)); } } From de7f02e48b79f398d3de08971648c15386e5f433 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 7 Jul 2020 09:54:16 -0400 Subject: [PATCH 10/29] fix typo --- python/lammps.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/lammps.py b/python/lammps.py index feda10802e..8bd1fc693b 100644 --- a/python/lammps.py +++ b/python/lammps.py @@ -204,7 +204,7 @@ class lammps(object): self.lib.lammps_neighlist_element_neighbors.argtypes = [c_void_p, c_int, c_int, POINTER(c_int), POINTER(c_int), POINTER(POINTER(c_int))] self.lib.lammps_neighlist_element_neighbors.restype = None - self.lib.lammps_has_error.argtypes = [v_void_p] + self.lib.lammps_has_error.argtypes = [c_void_p] self.lib.lammps_has_error.restype = c_bool self.lib.lammps_get_last_error_message.argtypes = [c_void_p, c_char_p, c_int] From fcc6ed3a588e8dc5d9da1721d473099e3e01cc03 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Tue, 7 Jul 2020 08:37:48 -0600 Subject: [PATCH 11/29] add offset option to reset_mol_ids command --- doc/src/reset_mol_ids.rst | 36 +++++++++++++++++-------- src/reset_mol_ids.cpp | 55 +++++++++++++++++++++++++++------------ 2 files changed, 64 insertions(+), 27 deletions(-) diff --git a/doc/src/reset_mol_ids.rst b/doc/src/reset_mol_ids.rst index 8d2d6eeab5..0bfcccb219 100644 --- a/doc/src/reset_mol_ids.rst +++ b/doc/src/reset_mol_ids.rst @@ -11,33 +11,47 @@ Syntax .. parsed-literal:: - reset_mol_ids group-ID + reset_mol_ids group-ID keyword value ... * group-ID = ID of group of atoms whose molecule IDs will be reset +* zero or more keyword/value pairs may be appended +* keyword = *offset* + .. parsed-literal:: + + *offset* value = *Noffset* + Examples """""""" .. code-block:: LAMMPS reset_mol_ids all - reset_mol_ids solvent + reset_mol_ids solvent offset 1000 Description """"""""""" Reset molecule IDs for a group of atoms. This will create a new set of molecule IDs that are numbered contiguously from 1 to N, if there -are N different molecule IDs used by the group. Only molecule IDs for -atoms in the specified group are reset. +are N different molecules in the group. Only molecule IDs for atoms +in the specified group are reset. -This can be useful to invoke after performing a reactive molecular -dynamics run with :doc:`fix bond/react `, :doc:`fix -bond/create `, or :doc:`fix bond/break -`. It can also be useful after molecules have been -deleted with :doc:`delete_atoms ` or after a simulation -which has lost molecules, e.g. via the :doc:`fix evaporate -` command. +For purposes of this operation, molecules are identified by the bond +topology of the system, not by the current molecule IDs. A molecule +is a set of atoms, each is which is bonded to one or more atoms in the +set. If an atom is not bonded to any other atom, it becomes its own +1-atom molecule. Once new molecules are identified, this command will +overwrite the current molecule ID for each atom with a new ID. + +This can be a useful operation to perform after running reactive +molecular dynamics run with :doc:`fix bond/react `, +:doc:`fix bond/create `, or :doc:`fix bond/break +`, all of which can change molecule topologies. It can +also be useful after molecules have been deleted with the +:doc:`delete_atoms ` command or after a simulation which +has lost molecules, e.g. via the :doc:`fix evaporate ` +command. Restrictions """""""""""" diff --git a/src/reset_mol_ids.cpp b/src/reset_mol_ids.cpp index 81ecf12a42..0b51dd08c8 100644 --- a/src/reset_mol_ids.cpp +++ b/src/reset_mol_ids.cpp @@ -46,12 +46,26 @@ void ResetMolIDs::command(int narg, char **arg) if (atom->molecular != 1) error->all(FLERR,"Can only use reset_mol_ids on molecular systems"); - if (narg != 1) error->all(FLERR,"Illegal reset_mol_ids command"); + // process args + + if (narg < 1) error->all(FLERR,"Illegal reset_mol_ids command"); int igroup = group->find(arg[0]); if (igroup == -1) error->all(FLERR,"Could not find reset_mol_ids group ID"); + int groupbit = group->bitmask[igroup]; - if (comm->me == 0) - utils::logmesg(lmp,"Resetting molecule IDs ...\n"); + bigint offset = 0; + + int iarg = 1; + while (iarg < narg) { + if (strcmp(arg[iarg],"offset") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal reset_mol_ids command"); + offset = utils::tnumeric(FLERR,arg[iarg+1],1,lmp); + if (offset < 0) error->all(FLERR,"Illegal reset_mol_ids command"); + iarg += 2; + } else error->all(FLERR,"Illegal reset_mol_ids command"); + } + + if (comm->me == 0) utils::logmesg(lmp,"Resetting molecule IDs ...\n"); // record wall time for resetting molecule IDs @@ -59,13 +73,14 @@ void ResetMolIDs::command(int narg, char **arg) double time1 = MPI_Wtime(); // create instances of compute fragment/atom and compute chunk/atom - + // both use the group-ID for this command + const std::string idfrag = "reset_mol_ids_FRAGMENT_ATOM"; - modify->add_compute(fmt::format("{} {} fragment/atom",idfrag, arg[0])); + modify->add_compute(fmt::format("{} {} fragment/atom",idfrag,arg[0])); const std::string idchunk = "reset_mol_ids_CHUNK_ATOM"; - modify->add_compute(fmt::format("{} all chunk/atom molecule compress yes", - idchunk)); + modify->add_compute(fmt::format("{} {} chunk/atom molecule compress yes", + idchunk,arg[0])); // initialize system since comm->borders() will be invoked @@ -92,27 +107,33 @@ void ResetMolIDs::command(int narg, char **arg) double *ids = cfa->vector_atom; // copy fragmentID to molecule ID - + // only for atoms in the group + tagint *molecule = atom->molecule; + int *mask = atom->mask; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) - molecule[i] = static_cast (ids[i]); + if (mask[i] & groupbit) + molecule[i] = static_cast (ids[i]); // invoke peratom method of compute chunk/atom - // compresses new molecule IDs to be contiguous 1 to Nmol - + // compress new molecule IDs to be contiguous 1 to Nmol + // NOTE: use of compute chunk/atom limits # of molecules to 32-bit int + icompute = modify->find_compute(idchunk); ComputeChunkAtom *cca = (ComputeChunkAtom *) modify->compute[icompute]; cca->compute_peratom(); ids = cca->vector_atom; int nchunk = cca->nchunk; - // copy chunkID to molecule ID - + // copy chunkID to molecule ID + offset + // only for atoms in the group + for (int i = 0; i < nlocal; i++) - molecule[i] = static_cast (ids[i]); - + if (mask[i] & groupbit) + molecule[i] = static_cast (ids[i]) + offset; + // clean up modify->delete_compute(idchunk); @@ -122,7 +143,9 @@ void ResetMolIDs::command(int narg, char **arg) MPI_Barrier(world); - if (comm->me == 0) + if (comm->me == 0) { + utils::logmesg(lmp,fmt::format(" number of new molecule IDs = {}\n",nchunk)); utils::logmesg(lmp,fmt::format(" reset_mol_ids CPU = {:.3f} seconds\n", MPI_Wtime()-time1)); + } } From 7870a89133f4e750f80af46cff1d2c469adb00d1 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Tue, 7 Jul 2020 08:48:24 -0600 Subject: [PATCH 12/29] added note to doc page --- doc/src/reset_mol_ids.rst | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/doc/src/reset_mol_ids.rst b/doc/src/reset_mol_ids.rst index 0bfcccb219..3c56080002 100644 --- a/doc/src/reset_mol_ids.rst +++ b/doc/src/reset_mol_ids.rst @@ -53,6 +53,15 @@ also be useful after molecules have been deleted with the has lost molecules, e.g. via the :doc:`fix evaporate ` command. +.. note:: + + The same as explained for the :doc:`compute fragment/atom + ` command, molecules are identified using the + current bond topology within each fragment. This will not account + for bonds broken by the :doc:`bond_style quartic ` + command because it does not perform a full update of the bond + topology data structures within LAMMPS. + Restrictions """""""""""" none From 0944eda3911ee93817a86a05831d7ca1f4a9d215 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Tue, 7 Jul 2020 09:07:48 -0600 Subject: [PATCH 13/29] final details for offset keyword added to reset_mol_ids --- doc/src/reset_mol_ids.rst | 18 +++++++++++++++++- src/reset_mol_ids.cpp | 15 +++++++++++++-- 2 files changed, 30 insertions(+), 3 deletions(-) diff --git a/doc/src/reset_mol_ids.rst b/doc/src/reset_mol_ids.rst index 3c56080002..2bb4557e60 100644 --- a/doc/src/reset_mol_ids.rst +++ b/doc/src/reset_mol_ids.rst @@ -53,6 +53,19 @@ also be useful after molecules have been deleted with the has lost molecules, e.g. via the :doc:`fix evaporate ` command. +The *offset* keyword can be used to change the range of new molecule +IDs assigned. If *Noffset* = 0 (the default) and the specified group +is *all*, then new molecule IDs will be from 1 to N. If *Noffset* = 0 +and the group is not all, then new molecule IDs will be from M+1 to +M+N, where M is the largest molecule ID for any atom not in the group. +This insures new molecule IDs do not collide with existing molecule +IDs that are not changed by this command. + +If *Noffset* is set to a value greater than 0, then new molecule IDs +will be from Noffset+1 to Noffset+N. If the group is not all, it is +up to you to insure the choice of *Noffset* does not produce +collisions between existing and new molecule IDs. + .. note:: The same as explained for the :doc:`compute fragment/atom @@ -75,4 +88,7 @@ Related commands :doc:`fix evaporate `, :doc:`delete_atoms ` -**Default:** none +**Default:** + +The default keyword setting is offset = 0. + diff --git a/src/reset_mol_ids.cpp b/src/reset_mol_ids.cpp index 0b51dd08c8..2f9591a594 100644 --- a/src/reset_mol_ids.cpp +++ b/src/reset_mol_ids.cpp @@ -53,7 +53,7 @@ void ResetMolIDs::command(int narg, char **arg) if (igroup == -1) error->all(FLERR,"Could not find reset_mol_ids group ID"); int groupbit = group->bitmask[igroup]; - bigint offset = 0; + tagint offset = 0; int iarg = 1; while (iarg < narg) { @@ -119,13 +119,24 @@ void ResetMolIDs::command(int narg, char **arg) // invoke peratom method of compute chunk/atom // compress new molecule IDs to be contiguous 1 to Nmol - // NOTE: use of compute chunk/atom limits # of molecules to 32-bit int + // NOTE: use of compute chunk/atom limits # of molecules to a 32-bit int icompute = modify->find_compute(idchunk); ComputeChunkAtom *cca = (ComputeChunkAtom *) modify->compute[icompute]; cca->compute_peratom(); ids = cca->vector_atom; int nchunk = cca->nchunk; + + // if offset = 0 and group != all, + // reset offset to largest molID of non-group atoms + + if (offset == 0 && groupbit != 1) { + tagint mymol = 0; + for (int i = 0; i < nlocal; i++) + if (!(mask[i] & groupbit)) + mymol = MAX(mymol,molecule[i]); + MPI_Allreduce(&mymol,&offset,1,MPI_LMP_TAGINT,MPI_MAX,world); + } // copy chunkID to molecule ID + offset // only for atoms in the group From 94e9b3bc82e6e565316e17ca91afc9eaf1ec7bce Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 7 Jul 2020 13:25:24 -0400 Subject: [PATCH 14/29] rename reset_ids command to reset_atom_ids --- doc/src/Commands_all.rst | 2 +- doc/src/commands_list.rst | 2 +- doc/src/delete_atoms.rst | 4 ++-- doc/src/{reset_ids.rst => reset_atom_ids.rst} | 14 +++++++------- doc/src/reset_mol_ids.rst | 2 +- src/deprecated.cpp | 17 +++++++++-------- src/deprecated.h | 1 + src/{reset_ids.cpp => reset_atom_ids.cpp} | 18 +++++++++--------- src/{reset_ids.h => reset_atom_ids.h} | 4 ++-- unittest/commands/CMakeLists.txt | 8 ++++---- ...st_reset_mol_ids.cpp => test_reset_ids.cpp} | 10 +++++----- unittest/formats/test_atom_styles.cpp | 6 +++--- 12 files changed, 45 insertions(+), 43 deletions(-) rename doc/src/{reset_ids.rst => reset_atom_ids.rst} (92%) rename src/{reset_ids.cpp => reset_atom_ids.cpp} (96%) rename src/{reset_ids.h => reset_atom_ids.h} (94%) rename unittest/commands/{test_reset_mol_ids.cpp => test_reset_ids.cpp} (96%) diff --git a/doc/src/Commands_all.rst b/doc/src/Commands_all.rst index 068545564a..a38bd5f0db 100644 --- a/doc/src/Commands_all.rst +++ b/doc/src/Commands_all.rst @@ -101,7 +101,7 @@ An alphabetic list of all general LAMMPS commands. * :doc:`region ` * :doc:`replicate ` * :doc:`rerun ` - * :doc:`reset_ids ` + * :doc:`reset_atom_ids ` * :doc:`reset_mol_ids ` * :doc:`reset_timestep ` * :doc:`restart ` diff --git a/doc/src/commands_list.rst b/doc/src/commands_list.rst index 49cf7ea588..2ec20ac220 100644 --- a/doc/src/commands_list.rst +++ b/doc/src/commands_list.rst @@ -88,7 +88,7 @@ Commands region replicate rerun - reset_ids + reset_atom_ids reset_mol_ids reset_timestep restart diff --git a/doc/src/delete_atoms.rst b/doc/src/delete_atoms.rst index d695603ded..127518535f 100644 --- a/doc/src/delete_atoms.rst +++ b/doc/src/delete_atoms.rst @@ -89,7 +89,7 @@ number of atoms in the system. Note that this is not done for molecular systems (see the :doc:`atom_style ` command), regardless of the *compress* setting, since it would foul up the bond connectivity that has already been assigned. However, the -:doc:`reset_ids ` command can be used after this command to +:doc:`reset_atom_ids ` command can be used after this command to accomplish the same thing. Note that the re-assignment of IDs is not really a compression, where @@ -157,7 +157,7 @@ using molecule template files via the :doc:`molecule ` and Related commands """""""""""""""" -:doc:`create_atoms `, :doc:`reset_ids ` +:doc:`create_atoms `, :doc:`reset_atom_ids ` Default """"""" diff --git a/doc/src/reset_ids.rst b/doc/src/reset_atom_ids.rst similarity index 92% rename from doc/src/reset_ids.rst rename to doc/src/reset_atom_ids.rst index bc7ad31927..e83d65d546 100644 --- a/doc/src/reset_ids.rst +++ b/doc/src/reset_atom_ids.rst @@ -1,14 +1,14 @@ -.. index:: reset_ids +.. index:: reset_atom_ids -reset_ids command -================= +reset_atom_ids command +====================== Syntax """""" .. code-block:: LAMMPS - reset_ids keyword values ... + reset_atom_ids keyword values ... * zero or more keyword/value pairs may be appended * keyword = *sort* @@ -22,8 +22,8 @@ Examples .. code-block:: LAMMPS - reset_ids - reset_ids sort yes + reset_atom_ids + reset_atom_ids sort yes Description """"""""""" @@ -77,7 +77,7 @@ processor have consecutive IDs, as the :doc:`create_atoms that are owned by other processors. The :doc:`comm_modify cutoff ` command can be used to correct this issue. Or you can define a pair style before using this command. If you do the former, you should unset the comm_modify cutoff after using - reset_ids so that subsequent communication is not inefficient. + reset_atom_ids so that subsequent communication is not inefficient. Restrictions """""""""""" diff --git a/doc/src/reset_mol_ids.rst b/doc/src/reset_mol_ids.rst index 8d2d6eeab5..9f2cc8db8a 100644 --- a/doc/src/reset_mol_ids.rst +++ b/doc/src/reset_mol_ids.rst @@ -46,7 +46,7 @@ none Related commands """""""""""""""" -:doc:`reset_ids `, :doc:`fix bond/react `, +:doc:`reset_atom_ids `, :doc:`fix bond/react `, :doc:`fix bond/create `, :doc:`fix bond/break `, :doc:`fix evaporate `, diff --git a/src/deprecated.cpp b/src/deprecated.cpp index 86af54fbfd..8cd3fc8cb1 100644 --- a/src/deprecated.cpp +++ b/src/deprecated.cpp @@ -16,19 +16,17 @@ ------------------------------------------------------------------------- */ #include "deprecated.h" -#include +#include #include "comm.h" #include "error.h" #include "input.h" +#include "utils.h" using namespace LAMMPS_NS; -static void writemsg(LAMMPS *lmp, const char *msg, int abend=1) +static void writemsg(LAMMPS *lmp, const std::string &msg, int abend=1) { - if (lmp->comm->me == 0) { - if (lmp->screen) fputs(msg,lmp->screen); - if (lmp->logfile) fputs(msg,lmp->logfile); - } + if (lmp->comm->me == 0) utils::logmesg(lmp,msg); if (abend) lmp->error->all(FLERR,"This command is no longer available"); } @@ -37,8 +35,11 @@ static void writemsg(LAMMPS *lmp, const char *msg, int abend=1) void Deprecated::command(int /* narg */, char ** /* arg */) { - if (strcmp(input->command,"DEPRECATED") == 0) { - writemsg(lmp,"\nCommand 'DEPRECATED' is a dummy command\n\n",0); + const std::string cmd = input->command; + if (cmd == "DEPRECATED") { + writemsg(lmp,"\nCommand 'DEPRECATED' is a dummy command\n\n",0); + } else if (cmd == "reset_ids") { + writemsg(lmp,"\n'reset_ids' has been renamed to 'reset_atom_ids'\n\n"); } } diff --git a/src/deprecated.h b/src/deprecated.h index 97a77591c5..4a7fbc44cf 100644 --- a/src/deprecated.h +++ b/src/deprecated.h @@ -14,6 +14,7 @@ #ifdef COMMAND_CLASS CommandStyle(DEPRECATED,Deprecated) +CommandStyle(reset_ids,Deprecated) #else diff --git a/src/reset_ids.cpp b/src/reset_atom_ids.cpp similarity index 96% rename from src/reset_ids.cpp rename to src/reset_atom_ids.cpp index 6a7a8e8b2a..7e493ca28e 100644 --- a/src/reset_ids.cpp +++ b/src/reset_atom_ids.cpp @@ -11,7 +11,7 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#include "reset_ids.h" +#include "reset_atom_ids.h" #include #include #include "atom.h" @@ -53,11 +53,11 @@ void ResetIDs::command(int narg, char **arg) if (domain->box_exist == 0) error->all(FLERR,"Reset_ids command before simulation box is defined"); if (atom->tag_enable == 0) - error->all(FLERR,"Cannot use reset_ids unless atoms have IDs"); + error->all(FLERR,"Cannot use reset_atom_ids unless atoms have IDs"); for (int i = 0; i < modify->nfix; i++) if (modify->fix[i]->stores_ids) - error->all(FLERR,"Cannot use reset_ids when a fix exists that stores atom IDs"); + error->all(FLERR,"Cannot use reset_atom_ids when a fix exists that stores atom IDs"); if (comm->me == 0) utils::logmesg(lmp,"Resetting atom IDs ...\n"); @@ -68,12 +68,12 @@ void ResetIDs::command(int narg, char **arg) int iarg = 0; while (iarg < narg) { if (strcmp(arg[iarg],"sort") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal reset_ids command"); + if (iarg+2 > narg) error->all(FLERR,"Illegal reset_atom_ids command"); if (strcmp(arg[iarg+1],"yes") == 0) sortflag = 1; else if (strcmp(arg[iarg+1],"no") == 0) sortflag = 0; - else error->all(FLERR,"Illegal reset_ids command"); + else error->all(FLERR,"Illegal reset_atom_ids command"); iarg += 2; - } else error->all(FLERR,"Illegal reset_ids command"); + } else error->all(FLERR,"Illegal reset_atom_ids command"); } // create an atom map if one doesn't exist already @@ -110,7 +110,7 @@ void ResetIDs::command(int narg, char **arg) int nall = nlocal + atom->nghost; tagint *oldIDs; - memory->create(oldIDs,nlocal,"reset_ids:oldIDs"); + memory->create(oldIDs,nlocal,"reset_atom_ids:oldIDs"); for (int i = 0; i < nlocal; i++) { oldIDs[i] = tag[i]; @@ -129,7 +129,7 @@ void ResetIDs::command(int narg, char **arg) // forward_comm_array acquires new IDs for ghost atoms double **newIDs; - memory->create(newIDs,nall,1,"reset_ids:newIDs"); + memory->create(newIDs,nall,1,"reset_atom_ids:newIDs"); for (int i = 0; i < nlocal; i++) { newIDs[i][0] = ubuf(tag[i]).d; @@ -468,7 +468,7 @@ int ResetIDs::sort_bins(int n, char *inbuf, if (in[i].ibin < binlo || in[i].ibin >= binhi) { //printf("BAD me %d i %d %d ibin %d binlohi %d %d\n", // rptr->comm->me,i,n,in[i].ibin,binlo,binhi); - error->one(FLERR,"Bad spatial bin assignment in reset_ids sort"); + error->one(FLERR,"Bad spatial bin assignment in reset_atom_ids sort"); } ibin = in[i].ibin - binlo; if (head[ibin] < 0) head[ibin] = i; diff --git a/src/reset_ids.h b/src/reset_atom_ids.h similarity index 94% rename from src/reset_ids.h rename to src/reset_atom_ids.h index 902fe51519..7c5c53e2ba 100644 --- a/src/reset_ids.h +++ b/src/reset_atom_ids.h @@ -13,7 +13,7 @@ #ifdef COMMAND_CLASS -CommandStyle(reset_ids,ResetIDs) +CommandStyle(reset_atom_ids,ResetIDs) #else @@ -70,7 +70,7 @@ E: Illegal ... command UNDOCUMENTED -E: Cannot use reset_ids unless atoms have IDs +E: Cannot use reset_atom_ids unless atoms have IDs UNDOCUMENTED diff --git a/unittest/commands/CMakeLists.txt b/unittest/commands/CMakeLists.txt index d0793b707f..275e4b973f 100644 --- a/unittest/commands/CMakeLists.txt +++ b/unittest/commands/CMakeLists.txt @@ -3,7 +3,7 @@ add_executable(test_simple_commands test_simple_commands.cpp) target_link_libraries(test_simple_commands PRIVATE lammps GTest::GMock GTest::GTest) add_test(NAME SimpleCommands COMMAND test_simple_commands WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) -add_executable(test_reset_mol_ids test_reset_mol_ids.cpp) -target_compile_definitions(test_reset_mol_ids PRIVATE -DTEST_INPUT_FOLDER=${CMAKE_CURRENT_SOURCE_DIR}) -target_link_libraries(test_reset_mol_ids PRIVATE lammps GTest::GMock GTest::GTest) -add_test(NAME ResetMolIDs COMMAND test_reset_mol_ids WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) +add_executable(test_reset_ids test_reset_ids.cpp) +target_compile_definitions(test_reset_ids PRIVATE -DTEST_INPUT_FOLDER=${CMAKE_CURRENT_SOURCE_DIR}) +target_link_libraries(test_reset_ids PRIVATE lammps GTest::GMock GTest::GTest) +add_test(NAME ResetIDs COMMAND test_reset_ids WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) diff --git a/unittest/commands/test_reset_mol_ids.cpp b/unittest/commands/test_reset_ids.cpp similarity index 96% rename from unittest/commands/test_reset_mol_ids.cpp rename to unittest/commands/test_reset_ids.cpp index bc38c6404c..198d7d3af2 100644 --- a/unittest/commands/test_reset_mol_ids.cpp +++ b/unittest/commands/test_reset_ids.cpp @@ -45,13 +45,13 @@ using ::testing::MatchesRegex; #define STRINGIFY(val) XSTR(val) #define XSTR(val) #val -class ResetMolIDsTest : public ::testing::Test { +class ResetIDsTest : public ::testing::Test { protected: LAMMPS *lmp; void SetUp() override { - const char *args[] = {"ResetMolIDsTest", "-log", "none", "-nocite", "-echo", "screen"}; + const char *args[] = {"ResetIDsTest", "-log", "none", "-nocite", "-echo", "screen"}; char **argv = (char **)args; int argc = sizeof(args) / sizeof(char *); if (!verbose) ::testing::internal::CaptureStdout(); @@ -72,7 +72,7 @@ protected: } }; -TEST_F(ResetMolIDsTest, Plain) +TEST_F(ResetIDsTest, MolIDAll) { if (lmp->atom->natoms == 0) GTEST_SKIP(); @@ -144,7 +144,7 @@ TEST_F(ResetMolIDsTest, Plain) ASSERT_EQ(molid[GETIDX(29)], 5); } -TEST_F(ResetMolIDsTest, DeletePlusAtomID) +TEST_F(ResetIDsTest, DeletePlusAtomID) { if (lmp->atom->natoms == 0) GTEST_SKIP(); @@ -214,7 +214,7 @@ TEST_F(ResetMolIDsTest, DeletePlusAtomID) EXPECT_GE(GETIDX(26), 0); if (!verbose) ::testing::internal::CaptureStdout(); - lmp->input->one("reset_ids"); + lmp->input->one("reset_atom_ids"); if (!verbose) ::testing::internal::GetCapturedStdout(); ASSERT_EQ(lmp->atom->map_tag_max, 23); diff --git a/unittest/formats/test_atom_styles.cpp b/unittest/formats/test_atom_styles.cpp index a405a240aa..491504d804 100644 --- a/unittest/formats/test_atom_styles.cpp +++ b/unittest/formats/test_atom_styles.cpp @@ -415,7 +415,7 @@ TEST_F(AtomStyleTest, atomic) ASSERT_EQ(lmp->atom->map_user, 1); ASSERT_EQ(lmp->atom->map_tag_max, 3); if (!verbose) ::testing::internal::CaptureStdout(); - lmp->input->one("reset_ids"); + lmp->input->one("reset_atom_ids"); if (!verbose) ::testing::internal::GetCapturedStdout(); ASSERT_EQ(lmp->atom->map_tag_max, 2); ASSERT_EQ(lmp->atom->tag_consecutive(), 1); @@ -802,7 +802,7 @@ TEST_F(AtomStyleTest, charge) ASSERT_EQ(lmp->atom->mass_setflag[2], 1); if (!verbose) ::testing::internal::CaptureStdout(); - lmp->input->one("reset_ids"); + lmp->input->one("reset_atom_ids"); lmp->input->one("replicate 2 2 2"); if (!verbose) ::testing::internal::GetCapturedStdout(); ASSERT_EQ(lmp->atom->map_tag_max, 16); @@ -1129,7 +1129,7 @@ TEST_F(AtomStyleTest, sphere) EXPECT_THAT(std::string(lmp->atom->atom_style), Eq("atomic")); lmp->input->one("read_restart test_atom_styles.restart"); lmp->input->one("replicate 1 1 2"); - lmp->input->one("reset_ids"); + lmp->input->one("reset_atom_ids"); if (!verbose) ::testing::internal::GetCapturedStdout(); ASSERT_THAT(std::string(lmp->atom->atom_style), Eq("sphere")); ASSERT_NE(lmp->atom->avec, nullptr); From 37d56a6bf6708dd0f456c2dbe0217b422daaa70c Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 7 Jul 2020 13:40:15 -0400 Subject: [PATCH 15/29] small cleanup in "deprecated" styles --- src/angle_deprecated.cpp | 3 +-- src/bond_deprecated.cpp | 3 +-- src/compute_deprecated.cpp | 3 +-- src/deprecated.cpp | 15 ++++++--------- src/dihedral_deprecated.cpp | 3 +-- src/dump_deprecated.cpp | 3 +-- src/fix_deprecated.cpp | 3 +-- src/improper_deprecated.cpp | 3 +-- src/kspace_deprecated.cpp | 3 +-- src/pair_deprecated.cpp | 3 +-- src/region_deprecated.cpp | 2 +- 11 files changed, 16 insertions(+), 28 deletions(-) diff --git a/src/angle_deprecated.cpp b/src/angle_deprecated.cpp index 968445f3f1..fb3f90ec28 100644 --- a/src/angle_deprecated.cpp +++ b/src/angle_deprecated.cpp @@ -44,8 +44,7 @@ void AngleDeprecated::settings(int, char **) utils::logmesg(lmp,"\nAngle style 'DEPRECATED' is a dummy style\n\n"); return; } - - lmp->error->all(FLERR,"This angle style is no longer available"); + error->all(FLERR,"This angle style is no longer available"); } diff --git a/src/bond_deprecated.cpp b/src/bond_deprecated.cpp index 7bd65b795d..0d91f73ea7 100644 --- a/src/bond_deprecated.cpp +++ b/src/bond_deprecated.cpp @@ -44,8 +44,7 @@ void BondDeprecated::settings(int, char **) utils::logmesg(lmp,"\nBond style 'DEPRECATED' is a dummy style\n\n"); return; } - - lmp->error->all(FLERR,"This bond style is no longer available"); + error->all(FLERR,"This bond style is no longer available"); } diff --git a/src/compute_deprecated.cpp b/src/compute_deprecated.cpp index 356e6b372a..2419e5bc8b 100644 --- a/src/compute_deprecated.cpp +++ b/src/compute_deprecated.cpp @@ -31,6 +31,5 @@ ComputeDeprecated::ComputeDeprecated(LAMMPS *lmp, int narg, char **arg) : utils::logmesg(lmp,"\nCompute style 'DEPRECATED' is a dummy style\n\n"); return; } - - lmp->error->all(FLERR,"This compute style is no longer available"); + error->all(FLERR,"This compute style is no longer available"); } diff --git a/src/deprecated.cpp b/src/deprecated.cpp index 8cd3fc8cb1..dc0ab400d4 100644 --- a/src/deprecated.cpp +++ b/src/deprecated.cpp @@ -24,13 +24,6 @@ using namespace LAMMPS_NS; -static void writemsg(LAMMPS *lmp, const std::string &msg, int abend=1) -{ - if (lmp->comm->me == 0) utils::logmesg(lmp,msg); - if (abend) - lmp->error->all(FLERR,"This command is no longer available"); -} - /* ---------------------------------------------------------------------- */ void Deprecated::command(int /* narg */, char ** /* arg */) @@ -38,8 +31,12 @@ void Deprecated::command(int /* narg */, char ** /* arg */) const std::string cmd = input->command; if (cmd == "DEPRECATED") { - writemsg(lmp,"\nCommand 'DEPRECATED' is a dummy command\n\n",0); + if (lmp->comm->me == 0) + utils::logmesg(lmp,"\nCommand 'DEPRECATED' is a dummy command\n\n"); + return; } else if (cmd == "reset_ids") { - writemsg(lmp,"\n'reset_ids' has been renamed to 'reset_atom_ids'\n\n"); + if (lmp->comm->me == 0) + utils::logmesg(lmp,"\n'reset_ids' has been renamed to 'reset_atom_ids'\n\n"); } + error->all(FLERR,"This command is no longer available"); } diff --git a/src/dihedral_deprecated.cpp b/src/dihedral_deprecated.cpp index a44dcae9d9..1709f7d3d7 100644 --- a/src/dihedral_deprecated.cpp +++ b/src/dihedral_deprecated.cpp @@ -45,6 +45,5 @@ void DihedralDeprecated::settings(int, char **) utils::logmesg(lmp,"\nDihedral style 'DEPRECATED' is a dummy style\n\n"); return; } - - lmp->error->all(FLERR,"This dihedral style is no longer available"); + error->all(FLERR,"This dihedral style is no longer available"); } diff --git a/src/dump_deprecated.cpp b/src/dump_deprecated.cpp index 3dae66248d..822bed7832 100644 --- a/src/dump_deprecated.cpp +++ b/src/dump_deprecated.cpp @@ -31,6 +31,5 @@ DumpDeprecated::DumpDeprecated(LAMMPS *lmp, int narg, char **arg) : utils::logmesg(lmp,"\nDump style 'DEPRECATED' is a dummy style\n\n"); return; } - - lmp->error->all(FLERR,"This dump style is no longer available"); + error->all(FLERR,"This dump style is no longer available"); } diff --git a/src/fix_deprecated.cpp b/src/fix_deprecated.cpp index c6313bd872..d3ceb69ea0 100644 --- a/src/fix_deprecated.cpp +++ b/src/fix_deprecated.cpp @@ -43,6 +43,5 @@ FixDeprecated::FixDeprecated(LAMMPS *lmp, int narg, char **arg) : "compute chunk/atom:\n dim, origin, delta, region, " "bound, discard, units\n\n"); } - - lmp->error->all(FLERR,"This fix style is no longer available"); + error->all(FLERR,"This fix style is no longer available"); } diff --git a/src/improper_deprecated.cpp b/src/improper_deprecated.cpp index ea089b5da6..0c1a45603e 100644 --- a/src/improper_deprecated.cpp +++ b/src/improper_deprecated.cpp @@ -45,8 +45,7 @@ void ImproperDeprecated::settings(int, char **) utils::logmesg(lmp,"\nImproper style 'DEPRECATED' is a dummy style\n\n"); return; } - - lmp->error->all(FLERR,"This improper style is no longer available"); + error->all(FLERR,"This improper style is no longer available"); } diff --git a/src/kspace_deprecated.cpp b/src/kspace_deprecated.cpp index 076bd8a186..5dc2b50e2b 100644 --- a/src/kspace_deprecated.cpp +++ b/src/kspace_deprecated.cpp @@ -35,8 +35,7 @@ void KSpaceDeprecated::settings(int, char **) utils::logmesg(lmp,"\nKSpace style 'DEPRECATED' is a dummy style\n\n"); return; } - - lmp->error->all(FLERR,"This kspace style is no longer available"); + error->all(FLERR,"This kspace style is no longer available"); } diff --git a/src/pair_deprecated.cpp b/src/pair_deprecated.cpp index d48638a456..75cb7efcdf 100644 --- a/src/pair_deprecated.cpp +++ b/src/pair_deprecated.cpp @@ -50,6 +50,5 @@ void PairDeprecated::settings(int, char **) utils::logmesg(lmp,"\nPair style 'reax' has been removed from LAMMPS " "after the 12 December 2018 version\n\n"); } - - lmp->error->all(FLERR,"This pair style is no longer available"); + error->all(FLERR,"This pair style is no longer available"); } diff --git a/src/region_deprecated.cpp b/src/region_deprecated.cpp index d0f1fb8e1f..afbe9e4189 100644 --- a/src/region_deprecated.cpp +++ b/src/region_deprecated.cpp @@ -31,5 +31,5 @@ RegionDeprecated::RegionDeprecated(LAMMPS *lmp, int narg, char **arg) : utils::logmesg(lmp,"\nRegion style 'DEPRECATED' is a dummy style\n\n"); return; } - lmp->error->all(FLERR,"This region style is no longer available"); + error->all(FLERR,"This region style is no longer available"); } From 12f62583f9ddba2b331e6c4b3f09c8898be19ab1 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 7 Jul 2020 13:56:50 -0400 Subject: [PATCH 16/29] whitespace cleanup --- src/reset_mol_ids.cpp | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/reset_mol_ids.cpp b/src/reset_mol_ids.cpp index 2f9591a594..4b26016125 100644 --- a/src/reset_mol_ids.cpp +++ b/src/reset_mol_ids.cpp @@ -47,14 +47,14 @@ void ResetMolIDs::command(int narg, char **arg) error->all(FLERR,"Can only use reset_mol_ids on molecular systems"); // process args - + if (narg < 1) error->all(FLERR,"Illegal reset_mol_ids command"); int igroup = group->find(arg[0]); if (igroup == -1) error->all(FLERR,"Could not find reset_mol_ids group ID"); int groupbit = group->bitmask[igroup]; tagint offset = 0; - + int iarg = 1; while (iarg < narg) { if (strcmp(arg[iarg],"offset") == 0) { @@ -74,7 +74,7 @@ void ResetMolIDs::command(int narg, char **arg) // create instances of compute fragment/atom and compute chunk/atom // both use the group-ID for this command - + const std::string idfrag = "reset_mol_ids_FRAGMENT_ATOM"; modify->add_compute(fmt::format("{} {} fragment/atom",idfrag,arg[0])); @@ -108,11 +108,11 @@ void ResetMolIDs::command(int narg, char **arg) // copy fragmentID to molecule ID // only for atoms in the group - + tagint *molecule = atom->molecule; int *mask = atom->mask; int nlocal = atom->nlocal; - + for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) molecule[i] = static_cast (ids[i]); @@ -120,7 +120,7 @@ void ResetMolIDs::command(int narg, char **arg) // invoke peratom method of compute chunk/atom // compress new molecule IDs to be contiguous 1 to Nmol // NOTE: use of compute chunk/atom limits # of molecules to a 32-bit int - + icompute = modify->find_compute(idchunk); ComputeChunkAtom *cca = (ComputeChunkAtom *) modify->compute[icompute]; cca->compute_peratom(); @@ -134,17 +134,17 @@ void ResetMolIDs::command(int narg, char **arg) tagint mymol = 0; for (int i = 0; i < nlocal; i++) if (!(mask[i] & groupbit)) - mymol = MAX(mymol,molecule[i]); + mymol = MAX(mymol,molecule[i]); MPI_Allreduce(&mymol,&offset,1,MPI_LMP_TAGINT,MPI_MAX,world); } - + // copy chunkID to molecule ID + offset // only for atoms in the group - + for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) molecule[i] = static_cast (ids[i]) + offset; - + // clean up modify->delete_compute(idchunk); @@ -153,7 +153,7 @@ void ResetMolIDs::command(int narg, char **arg) // total time MPI_Barrier(world); - + if (comm->me == 0) { utils::logmesg(lmp,fmt::format(" number of new molecule IDs = {}\n",nchunk)); utils::logmesg(lmp,fmt::format(" reset_mol_ids CPU = {:.3f} seconds\n", From d3853af4bedba45bf1308704bad29e528399bcfb Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 7 Jul 2020 17:13:43 -0400 Subject: [PATCH 17/29] add singlezero keyword to compute fragment/atom to give all single non-bonded atoms an ID of 0 --- doc/src/compute_cluster_atom.rst | 9 ++++- src/compute_fragment_atom.cpp | 60 ++++++++++++++++++++------------ src/compute_fragment_atom.h | 2 +- 3 files changed, 46 insertions(+), 25 deletions(-) diff --git a/doc/src/compute_cluster_atom.rst b/doc/src/compute_cluster_atom.rst index 502d2858d0..a4017171c9 100644 --- a/doc/src/compute_cluster_atom.rst +++ b/doc/src/compute_cluster_atom.rst @@ -15,12 +15,13 @@ Syntax .. parsed-literal:: compute ID group-ID cluster/atom cutoff - compute ID group-ID fragment/atom + compute ID group-ID fragment/atom [singlezero] compute ID group-ID aggregate/atom cutoff * ID, group-ID are documented in :doc:`compute ` command * *cluster/atom* or *fragment/atom* or *aggregate/atom* = style name of this compute command * cutoff = distance within which to label atoms as part of same cluster (distance units) +* *singlezero* = keyword to trigger assigning an ID of 0 to fragments with a single atom (optional) Examples """""""" @@ -52,6 +53,12 @@ bond/break `. The cluster ID or fragment ID of every atom in the cluster will be set to the smallest atom ID of any atom in the cluster or fragment, respectively. +The *singlezero* keyword turns on a special treatment for fragments, +where all fragments within the compute group that contain only a single +atom will have a fragment ID of 0. This can be useful in cases where +the fragment IDs are used as input for other commands in LAMMPS that +treat such single atoms different. + An aggregate is defined by combining the rules for clusters and fragments, i.e. a set of atoms, where each of it is within the cutoff distance from one or more atoms within a fragment that is part of diff --git a/src/compute_fragment_atom.cpp b/src/compute_fragment_atom.cpp index 6239ea9101..e71593ce3d 100644 --- a/src/compute_fragment_atom.cpp +++ b/src/compute_fragment_atom.cpp @@ -38,7 +38,12 @@ ComputeFragmentAtom::ComputeFragmentAtom(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), fragmentID(NULL) { - if (narg != 3) error->all(FLERR,"Illegal compute fragment/atom command"); + singleflag = 0; + + if ((narg < 3) || (narg > 4)) + error->all(FLERR,"Illegal compute fragment/atom command"); + if ((narg == 4) && (strcmp(arg[3],"singlezero") == 0)) + singleflag = 1; if (atom->avec->bonds_allow == 0) error->all(FLERR,"Compute fragment/atom used when bonds are not allowed"); @@ -114,6 +119,7 @@ void ComputeFragmentAtom::compute_peratom() // owned + ghost atoms start with fragmentID = atomID // atoms not in group have fragmentID = 0 + // if singleflag is set atoms without bonds have fragmentID 0 as well. tagint *tag = atom->tag; int *mask = atom->mask; @@ -122,9 +128,10 @@ void ComputeFragmentAtom::compute_peratom() int nlocal = atom->nlocal; int nall = nlocal + atom->nghost; - for (i = 0; i < nall; i++) + for (i = 0; i < nall; i++) { if (mask[i] & groupbit) fragmentID[i] = tag[i]; else fragmentID[i] = 0; + } // loop until no ghost atom fragment ID is changed // acquire fragmentIDs of ghost atoms @@ -152,9 +159,16 @@ void ComputeFragmentAtom::compute_peratom() for (i = 0; i < nlocal; i++) { // skip atom I if not in group or already marked + // also skip and set fragment ID to zero if singleflag is set + // and the atom is an isolated atom without bonds if (!(mask[i] & groupbit)) continue; if (markflag[i]) continue; + if (singleflag && (nspecial[i][0] == 0)) { + fragmentID[i] = 0; + markflag[i] = 1; + continue; + } // find one cluster of bond-connected atoms // ncluster = # of owned and ghost atoms in cluster @@ -165,44 +179,44 @@ void ComputeFragmentAtom::compute_peratom() stack[nstack++] = i; while (nstack) { - j = stack[--nstack]; - clist[ncluster++] = j; - markflag[j] = 1; + j = stack[--nstack]; + clist[ncluster++] = j; + markflag[j] = 1; - n = nspecial[j][0]; - list = special[j]; + n = nspecial[j][0]; + list = special[j]; for (m = 0; m < n; m++) { k = atom->map(list[m]); - // skip bond neighbor K if not in group or already marked - - if (k < 0) continue; - if (!(mask[k] & groupbit)) continue; - if (k < nlocal && markflag[k]) continue; + // skip bond neighbor K if not in group or already marked - // owned bond neighbors are added to stack for further walking - // ghost bond neighbors are added directly w/out use of stack - - if (k < nlocal) stack[nstack++] = k; - else clist[ncluster++] = k; - } + if (k < 0) continue; + if (!(mask[k] & groupbit)) continue; + if (k < nlocal && markflag[k]) continue; + + // owned bond neighbors are added to stack for further walking + // ghost bond neighbors are added directly w/out use of stack + + if (k < nlocal) stack[nstack++] = k; + else clist[ncluster++] = k; + } } // newID = minimum fragment ID in cluster list, including ghost atoms newID = BIG; for (m = 0; m < ncluster; m++) { - cID = fragmentID[clist[m]]; - newID = MIN(newID,cID); + cID = fragmentID[clist[m]]; + newID = MIN(newID,cID); } // set fragmentID = newID for all atoms in cluster, including ghost atoms // not done with iterations if change the fragmentID of a ghost atom for (m = 0; m < ncluster; m++) { - j = clist[m]; - if (j >= nlocal && fragmentID[j] != newID) done = 0; - fragmentID[j] = newID; + j = clist[m]; + if (j >= nlocal && fragmentID[j] != newID) done = 0; + fragmentID[j] = newID; } } diff --git a/src/compute_fragment_atom.h b/src/compute_fragment_atom.h index 1422a71fb4..e5cc2a5b6f 100644 --- a/src/compute_fragment_atom.h +++ b/src/compute_fragment_atom.h @@ -35,7 +35,7 @@ class ComputeFragmentAtom : public Compute { double memory_usage(); private: - int nmax,commflag; + int nmax,commflag,singleflag; int *stack,*clist,*markflag; double *fragmentID; }; From fd95fc98c5c8f11518bee9b47009cf83f41648b5 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 7 Jul 2020 17:14:29 -0400 Subject: [PATCH 18/29] add support for auto offset and singlezero option --- doc/src/reset_mol_ids.rst | 28 ++++++---- doc/utils/sphinx-config/false_positives.txt | 1 + src/reset_mol_ids.cpp | 59 +++++++++++++++++---- 3 files changed, 68 insertions(+), 20 deletions(-) diff --git a/doc/src/reset_mol_ids.rst b/doc/src/reset_mol_ids.rst index c6d1cfab03..b56596cc8a 100644 --- a/doc/src/reset_mol_ids.rst +++ b/doc/src/reset_mol_ids.rst @@ -15,19 +15,22 @@ Syntax * group-ID = ID of group of atoms whose molecule IDs will be reset * zero or more keyword/value pairs may be appended -* keyword = *offset* +* keyword = *offset* or *singlezero* .. parsed-literal:: - *offset* value = *Noffset* - + *offset* value = *Noffset* or auto + *singlezero* = assign molecule ID 0 to atoms without bonds + Examples """""""" .. code-block:: LAMMPS reset_mol_ids all + reset_mol_ids all offset 10 singlezero reset_mol_ids solvent offset 1000 + reset_mol_ids solvent offset auto Description """"""""""" @@ -39,7 +42,7 @@ in the specified group are reset. For purposes of this operation, molecules are identified by the bond topology of the system, not by the current molecule IDs. A molecule -is a set of atoms, each is which is bonded to one or more atoms in the +is a set of atoms, each of which is bonded to one or more atoms in the set. If an atom is not bonded to any other atom, it becomes its own 1-atom molecule. Once new molecules are identified, this command will overwrite the current molecule ID for each atom with a new ID. @@ -54,18 +57,25 @@ has lost molecules, e.g. via the :doc:`fix evaporate ` command. The *offset* keyword can be used to change the range of new molecule -IDs assigned. If *Noffset* = 0 (the default) and the specified group -is *all*, then new molecule IDs will be from 1 to N. If *Noffset* = 0 -and the group is not all, then new molecule IDs will be from M+1 to -M+N, where M is the largest molecule ID for any atom not in the group. +IDs assigned. If *Noffset* = auto (the default) and the specified group +is *all*, then new molecule IDs will be from 1 to N. If *Noffset* = auto +and the group is **not** all, then new molecule IDs will be from M+1 to +M+N, where M is the largest molecule ID for any atom **not** in the group. This insures new molecule IDs do not collide with existing molecule IDs that are not changed by this command. -If *Noffset* is set to a value greater than 0, then new molecule IDs +If *Noffset* is set to a number instead of 'auto', then new molecule IDs will be from Noffset+1 to Noffset+N. If the group is not all, it is up to you to insure the choice of *Noffset* does not produce collisions between existing and new molecule IDs. +The *singlezero* keyword turns on a special treatment for single atoms +without bonds. Instead of assigning each atom a different molecule ID +those atoms will all be given the molecule ID 0. This can be useful +when the molecule ID is used to group atoms where atoms with a molecule +ID of 0 will be considered as individual atoms; for example when using +:doc:`fix rigid ` with the *molecule* option. + .. note:: The same as explained for the :doc:`compute fragment/atom diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index 682969714c..be37f37c16 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -2058,6 +2058,7 @@ nodeless nodeset nodesets Noehring +Noffset noforce Noid nolib diff --git a/src/reset_mol_ids.cpp b/src/reset_mol_ids.cpp index 4b26016125..69c2939282 100644 --- a/src/reset_mol_ids.cpp +++ b/src/reset_mol_ids.cpp @@ -25,6 +25,7 @@ #include "modify.h" #include "compute_fragment_atom.h" #include "compute_chunk_atom.h" +#include "compute_reduce.h" #include "utils.h" #include "error.h" #include "fmt/format.h" @@ -53,15 +54,23 @@ void ResetMolIDs::command(int narg, char **arg) if (igroup == -1) error->all(FLERR,"Could not find reset_mol_ids group ID"); int groupbit = group->bitmask[igroup]; - tagint offset = 0; + tagint offset = -1; + int singleflag = 0; int iarg = 1; while (iarg < narg) { if (strcmp(arg[iarg],"offset") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal reset_mol_ids command"); - offset = utils::tnumeric(FLERR,arg[iarg+1],1,lmp); - if (offset < 0) error->all(FLERR,"Illegal reset_mol_ids command"); + if (strcmp(arg[iarg+1],"auto") == 0) { + offset = -1; + } else { + offset = utils::tnumeric(FLERR,arg[iarg+1],1,lmp); + if (offset < 0) error->all(FLERR,"Illegal reset_mol_ids command"); + } iarg += 2; + } else if (strcmp(arg[iarg],"singlezero") == 0) { + singleflag = 1; + ++iarg; } else error->all(FLERR,"Illegal reset_mol_ids command"); } @@ -72,11 +81,18 @@ void ResetMolIDs::command(int narg, char **arg) MPI_Barrier(world); double time1 = MPI_Wtime(); - // create instances of compute fragment/atom and compute chunk/atom - // both use the group-ID for this command + // create instances of compute fragment/atom, compute reduce (if needed), + // and compute chunk/atom. all use the group-ID for this command const std::string idfrag = "reset_mol_ids_FRAGMENT_ATOM"; - modify->add_compute(fmt::format("{} {} fragment/atom",idfrag,arg[0])); + if (singleflag) + modify->add_compute(fmt::format("{} {} fragment/atom singlezero",idfrag,arg[0])); + else + modify->add_compute(fmt::format("{} {} fragment/atom",idfrag,arg[0])); + + const std::string idmin = "reset_mol_ids_COMPUTE_MINFRAG"; + if (singleflag) + modify->add_compute(fmt::format("{} {} reduce min c_{}",idmin,arg[0],idfrag)); const std::string idchunk = "reset_mol_ids_CHUNK_ATOM"; modify->add_compute(fmt::format("{} {} chunk/atom molecule compress yes", @@ -106,6 +122,16 @@ void ResetMolIDs::command(int narg, char **arg) cfa->compute_peratom(); double *ids = cfa->vector_atom; + // when the lowest fragment ID is 0, in case the singlezero option is used + // we need to remove adjust the chunk ID, so that the resuling molecule ID is correct. + + int adjid = 0; + if (singleflag) { + icompute = modify->find_compute(idmin); + ComputeReduce *crd = (ComputeReduce *) modify->compute[icompute]; + if (crd->compute_scalar() == 0.0) adjid = 1; + } + // copy fragmentID to molecule ID // only for atoms in the group @@ -127,28 +153,39 @@ void ResetMolIDs::command(int narg, char **arg) ids = cca->vector_atom; int nchunk = cca->nchunk; - // if offset = 0 and group != all, + // if offset = -1 (i.e. "auto") and group != all, // reset offset to largest molID of non-group atoms + // otherwise set to 0. - if (offset == 0 && groupbit != 1) { + if (offset == -1) { + if (groupbit != 1) { tagint mymol = 0; for (int i = 0; i < nlocal; i++) if (!(mask[i] & groupbit)) mymol = MAX(mymol,molecule[i]); MPI_Allreduce(&mymol,&offset,1,MPI_LMP_TAGINT,MPI_MAX,world); + } else offset = 0; } // copy chunkID to molecule ID + offset // only for atoms in the group - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) - molecule[i] = static_cast (ids[i]) + offset; + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + tagint newid = static_cast(ids[i]); + if (singleflag && adjid) { + if (newid == 1) newid = 0; + else newid += offset - 1; + molecule[i] = newid; + } else molecule[i] = newid + offset; + } + } // clean up modify->delete_compute(idchunk); modify->delete_compute(idfrag); + if (singleflag) modify->delete_compute(idmin); // total time From 416467a154bb83db7fbcf4b35f29419e7f7481f1 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 7 Jul 2020 17:14:54 -0400 Subject: [PATCH 19/29] update reset_mol_ids tests for the latest changes --- unittest/commands/test_reset_ids.cpp | 349 ++++++++++++++++++++++----- 1 file changed, 295 insertions(+), 54 deletions(-) diff --git a/unittest/commands/test_reset_ids.cpp b/unittest/commands/test_reset_ids.cpp index 198d7d3af2..8929c35dab 100644 --- a/unittest/commands/test_reset_ids.cpp +++ b/unittest/commands/test_reset_ids.cpp @@ -55,9 +55,9 @@ protected: char **argv = (char **)args; int argc = sizeof(args) / sizeof(char *); if (!verbose) ::testing::internal::CaptureStdout(); - lmp = new LAMMPS(argc, argv, MPI_COMM_WORLD); + lmp = new LAMMPS(argc, argv, MPI_COMM_WORLD); Info *info = new Info(lmp); - if (info->has_style("atom","full")) { + if (info->has_style("atom", "full")) { lmp->input->one("variable input_dir index " STRINGIFY(TEST_INPUT_FOLDER)); lmp->input->one("include ${input_dir}/in.fourmol"); } @@ -160,58 +160,58 @@ TEST_F(ResetIDsTest, DeletePlusAtomID) ASSERT_EQ(lmp->atom->natoms, 23); ASSERT_EQ(lmp->atom->map_tag_max, 26); - EXPECT_EQ(molid[GETIDX(1)], 1); - EXPECT_EQ(molid[GETIDX(2)], 1); - EXPECT_EQ(molid[GETIDX(3)], 1); - EXPECT_EQ(molid[GETIDX(4)], 1); - EXPECT_EQ(molid[GETIDX(5)], 1); - EXPECT_EQ(molid[GETIDX(6)], 1); - EXPECT_EQ(molid[GETIDX(7)], 1); - EXPECT_EQ(molid[GETIDX(8)], 1); - EXPECT_EQ(molid[GETIDX(9)], 1); - EXPECT_EQ(molid[GETIDX(10)], 1); - EXPECT_EQ(molid[GETIDX(11)], 1); - EXPECT_EQ(molid[GETIDX(12)], 1); - EXPECT_EQ(molid[GETIDX(13)], 1); - EXPECT_EQ(molid[GETIDX(14)], 1); - EXPECT_EQ(molid[GETIDX(15)], 1); - EXPECT_EQ(molid[GETIDX(16)], 1); - EXPECT_EQ(molid[GETIDX(17)], 1); - EXPECT_EQ(molid[GETIDX(18)], 2); - EXPECT_EQ(molid[GETIDX(19)], 2); - EXPECT_EQ(molid[GETIDX(20)], 2); - EXPECT_EQ(molid[GETIDX(24)], 3); - EXPECT_EQ(molid[GETIDX(25)], 3); - EXPECT_EQ(molid[GETIDX(26)], 3); + ASSERT_EQ(molid[GETIDX(1)], 1); + ASSERT_EQ(molid[GETIDX(2)], 1); + ASSERT_EQ(molid[GETIDX(3)], 1); + ASSERT_EQ(molid[GETIDX(4)], 1); + ASSERT_EQ(molid[GETIDX(5)], 1); + ASSERT_EQ(molid[GETIDX(6)], 1); + ASSERT_EQ(molid[GETIDX(7)], 1); + ASSERT_EQ(molid[GETIDX(8)], 1); + ASSERT_EQ(molid[GETIDX(9)], 1); + ASSERT_EQ(molid[GETIDX(10)], 1); + ASSERT_EQ(molid[GETIDX(11)], 1); + ASSERT_EQ(molid[GETIDX(12)], 1); + ASSERT_EQ(molid[GETIDX(13)], 1); + ASSERT_EQ(molid[GETIDX(14)], 1); + ASSERT_EQ(molid[GETIDX(15)], 1); + ASSERT_EQ(molid[GETIDX(16)], 1); + ASSERT_EQ(molid[GETIDX(17)], 1); + ASSERT_EQ(molid[GETIDX(18)], 2); + ASSERT_EQ(molid[GETIDX(19)], 2); + ASSERT_EQ(molid[GETIDX(20)], 2); + ASSERT_EQ(molid[GETIDX(24)], 3); + ASSERT_EQ(molid[GETIDX(25)], 3); + ASSERT_EQ(molid[GETIDX(26)], 3); // now also check and reset the atom ids - EXPECT_GE(GETIDX(1), 0); - EXPECT_GE(GETIDX(2), 0); - EXPECT_GE(GETIDX(3), 0); - EXPECT_GE(GETIDX(4), 0); - EXPECT_GE(GETIDX(5), 0); - EXPECT_GE(GETIDX(6), 0); - EXPECT_GE(GETIDX(7), 0); - EXPECT_GE(GETIDX(8), 0); - EXPECT_GE(GETIDX(9), 0); - EXPECT_GE(GETIDX(10), 0); - EXPECT_GE(GETIDX(11), 0); - EXPECT_GE(GETIDX(12), 0); - EXPECT_GE(GETIDX(13), 0); - EXPECT_GE(GETIDX(14), 0); - EXPECT_GE(GETIDX(15), 0); - EXPECT_GE(GETIDX(16), 0); - EXPECT_GE(GETIDX(17), 0); - EXPECT_GE(GETIDX(18), 0); - EXPECT_GE(GETIDX(19), 0); - EXPECT_GE(GETIDX(20), 0); - EXPECT_EQ(GETIDX(21), -1); - EXPECT_EQ(GETIDX(22), -1); - EXPECT_EQ(GETIDX(23), -1); - EXPECT_GE(GETIDX(24), 0); - EXPECT_GE(GETIDX(25), 0); - EXPECT_GE(GETIDX(26), 0); + ASSERT_GE(GETIDX(1), 0); + ASSERT_GE(GETIDX(2), 0); + ASSERT_GE(GETIDX(3), 0); + ASSERT_GE(GETIDX(4), 0); + ASSERT_GE(GETIDX(5), 0); + ASSERT_GE(GETIDX(6), 0); + ASSERT_GE(GETIDX(7), 0); + ASSERT_GE(GETIDX(8), 0); + ASSERT_GE(GETIDX(9), 0); + ASSERT_GE(GETIDX(10), 0); + ASSERT_GE(GETIDX(11), 0); + ASSERT_GE(GETIDX(12), 0); + ASSERT_GE(GETIDX(13), 0); + ASSERT_GE(GETIDX(14), 0); + ASSERT_GE(GETIDX(15), 0); + ASSERT_GE(GETIDX(16), 0); + ASSERT_GE(GETIDX(17), 0); + ASSERT_GE(GETIDX(18), 0); + ASSERT_GE(GETIDX(19), 0); + ASSERT_GE(GETIDX(20), 0); + ASSERT_EQ(GETIDX(21), -1); + ASSERT_EQ(GETIDX(22), -1); + ASSERT_EQ(GETIDX(23), -1); + ASSERT_GE(GETIDX(24), 0); + ASSERT_GE(GETIDX(25), 0); + ASSERT_GE(GETIDX(26), 0); if (!verbose) ::testing::internal::CaptureStdout(); lmp->input->one("reset_atom_ids"); @@ -219,12 +219,253 @@ TEST_F(ResetIDsTest, DeletePlusAtomID) ASSERT_EQ(lmp->atom->map_tag_max, 23); for (int i = 1; i <= 23; ++i) - EXPECT_GE(GETIDX(i), 0); + ASSERT_GE(GETIDX(i), 0); } -// TEST_FAILURE(lmp->input->one("reset_mol_ids all");); -// ASSERT_THAT(mesg, MatchesRegex(".*ERROR: Unknown command.*")); +TEST_F(ResetIDsTest, PartialOffset) +{ + if (lmp->atom->natoms == 0) GTEST_SKIP(); + auto molid = lmp->atom->molecule; + + // delete two water molecules + if (!verbose) ::testing::internal::CaptureStdout(); + lmp->input->one("group allwater molecule 3:6"); + lmp->input->one("group nowater subtract all allwater"); + lmp->input->one("reset_mol_ids allwater offset 4"); + if (!verbose) ::testing::internal::GetCapturedStdout(); + ASSERT_EQ(lmp->atom->natoms, 29); + ASSERT_EQ(lmp->atom->map_tag_max, 29); + + ASSERT_EQ(molid[GETIDX(1)], 1); + ASSERT_EQ(molid[GETIDX(2)], 1); + ASSERT_EQ(molid[GETIDX(3)], 1); + ASSERT_EQ(molid[GETIDX(4)], 1); + ASSERT_EQ(molid[GETIDX(5)], 1); + ASSERT_EQ(molid[GETIDX(6)], 1); + ASSERT_EQ(molid[GETIDX(7)], 1); + ASSERT_EQ(molid[GETIDX(8)], 2); + ASSERT_EQ(molid[GETIDX(9)], 2); + ASSERT_EQ(molid[GETIDX(10)], 2); + ASSERT_EQ(molid[GETIDX(11)], 2); + ASSERT_EQ(molid[GETIDX(12)], 2); + ASSERT_EQ(molid[GETIDX(13)], 2); + ASSERT_EQ(molid[GETIDX(14)], 2); + ASSERT_EQ(molid[GETIDX(15)], 2); + ASSERT_EQ(molid[GETIDX(16)], 2); + ASSERT_EQ(molid[GETIDX(17)], 2); + ASSERT_EQ(molid[GETIDX(18)], 5); + ASSERT_EQ(molid[GETIDX(19)], 5); + ASSERT_EQ(molid[GETIDX(20)], 5); + ASSERT_EQ(molid[GETIDX(21)], 6); + ASSERT_EQ(molid[GETIDX(22)], 6); + ASSERT_EQ(molid[GETIDX(23)], 6); + ASSERT_EQ(molid[GETIDX(24)], 7); + ASSERT_EQ(molid[GETIDX(25)], 7); + ASSERT_EQ(molid[GETIDX(26)], 7); + ASSERT_EQ(molid[GETIDX(27)], 8); + ASSERT_EQ(molid[GETIDX(28)], 8); + ASSERT_EQ(molid[GETIDX(29)], 8); + + if (!verbose) ::testing::internal::CaptureStdout(); + lmp->input->one("reset_mol_ids nowater offset 0"); + if (!verbose) ::testing::internal::GetCapturedStdout(); + + ASSERT_EQ(molid[GETIDX(1)], 1); + ASSERT_EQ(molid[GETIDX(2)], 1); + ASSERT_EQ(molid[GETIDX(3)], 1); + ASSERT_EQ(molid[GETIDX(4)], 1); + ASSERT_EQ(molid[GETIDX(5)], 1); + ASSERT_EQ(molid[GETIDX(6)], 1); + ASSERT_EQ(molid[GETIDX(7)], 1); + ASSERT_EQ(molid[GETIDX(8)], 1); + ASSERT_EQ(molid[GETIDX(9)], 1); + ASSERT_EQ(molid[GETIDX(10)], 1); + ASSERT_EQ(molid[GETIDX(11)], 1); + ASSERT_EQ(molid[GETIDX(12)], 1); + ASSERT_EQ(molid[GETIDX(13)], 1); + ASSERT_EQ(molid[GETIDX(14)], 1); + ASSERT_EQ(molid[GETIDX(15)], 1); + ASSERT_EQ(molid[GETIDX(16)], 1); + ASSERT_EQ(molid[GETIDX(17)], 1); + ASSERT_EQ(molid[GETIDX(18)], 5); + ASSERT_EQ(molid[GETIDX(19)], 5); + ASSERT_EQ(molid[GETIDX(20)], 5); + ASSERT_EQ(molid[GETIDX(21)], 6); + ASSERT_EQ(molid[GETIDX(22)], 6); + ASSERT_EQ(molid[GETIDX(23)], 6); + ASSERT_EQ(molid[GETIDX(24)], 7); + ASSERT_EQ(molid[GETIDX(25)], 7); + ASSERT_EQ(molid[GETIDX(26)], 7); + ASSERT_EQ(molid[GETIDX(27)], 8); + ASSERT_EQ(molid[GETIDX(28)], 8); + ASSERT_EQ(molid[GETIDX(29)], 8); +} + +TEST_F(ResetIDsTest, DeleteAdd) +{ + if (lmp->atom->natoms == 0) GTEST_SKIP(); + + auto molid = lmp->atom->molecule; + + // delete two water molecules + if (!verbose) ::testing::internal::CaptureStdout(); + lmp->input->one("group allwater molecule 3:6"); + lmp->input->one("group twowater molecule 4:6:2"); + lmp->input->one("group nowater subtract all allwater"); + lmp->input->one("delete_atoms group twowater compress no bond yes mol yes"); + lmp->input->one("reset_mol_ids allwater offset auto"); + if (!verbose) ::testing::internal::GetCapturedStdout(); + ASSERT_EQ(lmp->atom->natoms, 23); + ASSERT_EQ(lmp->atom->map_tag_max, 26); + + ASSERT_EQ(molid[GETIDX(1)], 1); + ASSERT_EQ(molid[GETIDX(2)], 1); + ASSERT_EQ(molid[GETIDX(3)], 1); + ASSERT_EQ(molid[GETIDX(4)], 1); + ASSERT_EQ(molid[GETIDX(5)], 1); + ASSERT_EQ(molid[GETIDX(6)], 1); + ASSERT_EQ(molid[GETIDX(7)], 1); + ASSERT_EQ(molid[GETIDX(8)], 2); + ASSERT_EQ(molid[GETIDX(9)], 2); + ASSERT_EQ(molid[GETIDX(10)], 2); + ASSERT_EQ(molid[GETIDX(11)], 2); + ASSERT_EQ(molid[GETIDX(12)], 2); + ASSERT_EQ(molid[GETIDX(13)], 2); + ASSERT_EQ(molid[GETIDX(14)], 2); + ASSERT_EQ(molid[GETIDX(15)], 2); + ASSERT_EQ(molid[GETIDX(16)], 2); + ASSERT_EQ(molid[GETIDX(17)], 2); + ASSERT_EQ(molid[GETIDX(18)], 3); + ASSERT_EQ(molid[GETIDX(19)], 3); + ASSERT_EQ(molid[GETIDX(20)], 3); + ASSERT_EQ(molid[GETIDX(24)], 4); + ASSERT_EQ(molid[GETIDX(25)], 4); + ASSERT_EQ(molid[GETIDX(26)], 4); + + // now also check and reset the atom ids + + ASSERT_GE(GETIDX(1), 0); + ASSERT_GE(GETIDX(2), 0); + ASSERT_GE(GETIDX(3), 0); + ASSERT_GE(GETIDX(4), 0); + ASSERT_GE(GETIDX(5), 0); + ASSERT_GE(GETIDX(6), 0); + ASSERT_GE(GETIDX(7), 0); + ASSERT_GE(GETIDX(8), 0); + ASSERT_GE(GETIDX(9), 0); + ASSERT_GE(GETIDX(10), 0); + ASSERT_GE(GETIDX(11), 0); + ASSERT_GE(GETIDX(12), 0); + ASSERT_GE(GETIDX(13), 0); + ASSERT_GE(GETIDX(14), 0); + ASSERT_GE(GETIDX(15), 0); + ASSERT_GE(GETIDX(16), 0); + ASSERT_GE(GETIDX(17), 0); + ASSERT_GE(GETIDX(18), 0); + ASSERT_GE(GETIDX(19), 0); + ASSERT_GE(GETIDX(20), 0); + ASSERT_EQ(GETIDX(21), -1); + ASSERT_EQ(GETIDX(22), -1); + ASSERT_EQ(GETIDX(23), -1); + ASSERT_GE(GETIDX(24), 0); + ASSERT_GE(GETIDX(25), 0); + ASSERT_GE(GETIDX(26), 0); + + if (!verbose) ::testing::internal::CaptureStdout(); + lmp->input->one("reset_atom_ids sort yes"); + if (!verbose) ::testing::internal::GetCapturedStdout(); + + ASSERT_EQ(lmp->atom->map_tag_max, 23); + for (int i = 1; i <= 23; ++i) + ASSERT_GE(GETIDX(i), 0); + + if (!verbose) ::testing::internal::CaptureStdout(); + lmp->input->one("reset_mol_ids nowater offset 1"); + if (!verbose) ::testing::internal::GetCapturedStdout(); + + EXPECT_EQ(molid[GETIDX(1)], 2); + EXPECT_EQ(molid[GETIDX(2)], 2); + EXPECT_EQ(molid[GETIDX(3)], 2); + EXPECT_EQ(molid[GETIDX(4)], 2); + EXPECT_EQ(molid[GETIDX(5)], 2); + EXPECT_EQ(molid[GETIDX(6)], 2); + EXPECT_EQ(molid[GETIDX(7)], 2); + EXPECT_EQ(molid[GETIDX(8)], 2); + EXPECT_EQ(molid[GETIDX(9)], 2); + EXPECT_EQ(molid[GETIDX(10)], 2); + EXPECT_EQ(molid[GETIDX(11)], 2); + EXPECT_EQ(molid[GETIDX(12)], 2); + EXPECT_EQ(molid[GETIDX(13)], 3); + EXPECT_EQ(molid[GETIDX(14)], 3); + EXPECT_EQ(molid[GETIDX(15)], 3); + EXPECT_EQ(molid[GETIDX(16)], 2); + EXPECT_EQ(molid[GETIDX(17)], 2); + EXPECT_EQ(molid[GETIDX(18)], 2); + EXPECT_EQ(molid[GETIDX(19)], 2); + EXPECT_EQ(molid[GETIDX(20)], 2); + EXPECT_EQ(molid[GETIDX(21)], 4); + EXPECT_EQ(molid[GETIDX(22)], 4); + EXPECT_EQ(molid[GETIDX(23)], 4); + + if (!verbose) ::testing::internal::CaptureStdout(); + lmp->input->one("create_atoms 1 single 0.0 0.0 0.0"); + lmp->input->one("create_atoms 2 single 1.0 0.0 0.0"); + lmp->input->one("create_atoms 3 single 2.0 0.0 0.0"); + lmp->input->one("create_atoms 4 single 3.0 0.0 0.0"); + lmp->input->one("reset_mol_ids all"); + if (!verbose) ::testing::internal::GetCapturedStdout(); + EXPECT_EQ(lmp->atom->natoms, 27); + EXPECT_EQ(lmp->atom->map_tag_max, 27); + + EXPECT_EQ(molid[GETIDX(21)], 3); + EXPECT_EQ(molid[GETIDX(22)], 3); + EXPECT_EQ(molid[GETIDX(23)], 3); + EXPECT_EQ(molid[GETIDX(24)], 4); + EXPECT_EQ(molid[GETIDX(25)], 5); + EXPECT_EQ(molid[GETIDX(26)], 6); + EXPECT_EQ(molid[GETIDX(27)], 7); + + if (!verbose) ::testing::internal::CaptureStdout(); + lmp->input->one("reset_mol_ids all singlezero"); + if (!verbose) ::testing::internal::GetCapturedStdout(); + + EXPECT_EQ(molid[GETIDX(21)], 3); + EXPECT_EQ(molid[GETIDX(22)], 3); + EXPECT_EQ(molid[GETIDX(23)], 3); + EXPECT_EQ(molid[GETIDX(24)], 0); + EXPECT_EQ(molid[GETIDX(25)], 0); + EXPECT_EQ(molid[GETIDX(26)], 0); + EXPECT_EQ(molid[GETIDX(27)], 0); +} + +TEST_F(ResetIDsTest, DeathTests) +{ + ::testing::internal::CaptureStdout(); + TEST_FAILURE(lmp->input->one("reset_mol_ids");); + auto mesg = ::testing::internal::GetCapturedStdout(); + ASSERT_THAT(mesg, MatchesRegex(".*ERROR: Illegal reset_mol_ids command.*")); + + ::testing::internal::CaptureStdout(); + TEST_FAILURE(lmp->input->one("reset_mol_ids all offset 1 1");); + mesg = ::testing::internal::GetCapturedStdout(); + ASSERT_THAT(mesg, MatchesRegex(".*ERROR: Illegal reset_mol_ids command.*")); + + ::testing::internal::CaptureStdout(); + TEST_FAILURE(lmp->input->one("reset_mol_ids all offset -2");); + mesg = ::testing::internal::GetCapturedStdout(); + ASSERT_THAT(mesg, MatchesRegex(".*ERROR: Illegal reset_mol_ids command.*")); + + ::testing::internal::CaptureStdout(); + TEST_FAILURE(lmp->input->one("reset_mol_ids all offset xxx");); + mesg = ::testing::internal::GetCapturedStdout(); + ASSERT_THAT(mesg, MatchesRegex(".*ERROR on proc 0: Expected integer.*")); + + ::testing::internal::CaptureStdout(); + TEST_FAILURE(lmp->input->one("reset_mol_ids all offset");); + mesg = ::testing::internal::GetCapturedStdout(); + ASSERT_THAT(mesg, MatchesRegex(".*ERROR: Illegal reset_mol_ids command.*")); +} } // namespace LAMMPS_NS int main(int argc, char **argv) From 9e83279887ad791d15747eecb7f285ee581bc4a8 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Thu, 9 Jul 2020 15:40:24 -0600 Subject: [PATCH 20/29] minor adjustments to new reset_mol_ids command --- doc/src/compute_cluster_atom.rst | 72 +++++++------- doc/src/reset_mol_ids.rst | 84 +++++++++------- src/compute_fragment_atom.cpp | 36 ++++--- src/reset_mol_ids.cpp | 159 +++++++++++++++++-------------- 4 files changed, 196 insertions(+), 155 deletions(-) diff --git a/doc/src/compute_cluster_atom.rst b/doc/src/compute_cluster_atom.rst index a4017171c9..3772b0dd0b 100644 --- a/doc/src/compute_cluster_atom.rst +++ b/doc/src/compute_cluster_atom.rst @@ -15,13 +15,18 @@ Syntax .. parsed-literal:: compute ID group-ID cluster/atom cutoff - compute ID group-ID fragment/atom [singlezero] + compute ID group-ID fragment/atom keyword value ... compute ID group-ID aggregate/atom cutoff * ID, group-ID are documented in :doc:`compute ` command * *cluster/atom* or *fragment/atom* or *aggregate/atom* = style name of this compute command * cutoff = distance within which to label atoms as part of same cluster (distance units) -* *singlezero* = keyword to trigger assigning an ID of 0 to fragments with a single atom (optional) +* zero or more keyword/value pairs may be appended to *fragment/atom* +* keyword = *single* + + .. parsed-literal:: + + *single* value = *yes* or *no* to treat single atoms (no bonds) as fragments Examples """""""" @@ -30,13 +35,16 @@ Examples compute 1 all cluster/atom 3.5 compute 1 all fragment/atom + compute 1 all fragment/atom single no compute 1 all aggregate/atom 3.5 Description """"""""""" -Define a computation that assigns each atom a cluster, fragment, -or aggregate ID. +Define a computation that assigns each atom a cluster, fragment, or +aggregate ID. Only atoms in the compute group are clustered and +assigned cluster IDs. Atoms not in the compute group are assigned an +ID = 0. A cluster is defined as a set of atoms, each of which is within the cutoff distance from one or more other atoms in the cluster. If an @@ -44,20 +52,19 @@ atom has no neighbors within the cutoff distance, then it is a 1-atom cluster. A fragment is similarly defined as a set of atoms, each of which has a -bond to another atom in the fragment, as defined initially via the -:doc:`data file ` or :doc:`create_bonds ` -commands, or by fixes which dynamically create or break bonds like -:doc:`fix bond/react `, :doc:`fix bond/create -`, :doc:`fix bond/swap `, or :doc:`fix -bond/break `. The cluster ID or fragment ID of every -atom in the cluster will be set to the smallest atom ID of any atom in -the cluster or fragment, respectively. +bond to another atom in the fragment. Bonds can be defined initially +via the :doc:`data file ` or :doc:`create_bonds +` commands, or dynamically by fixes which create or +break bonds like :doc:`fix bond/react `, :doc:`fix +bond/create `, :doc:`fix bond/swap `, +or :doc:`fix bond/break `. The cluster ID or fragment +ID of every atom in the cluster will be set to the smallest atom ID of +any atom in the cluster or fragment, respectively. -The *singlezero* keyword turns on a special treatment for fragments, -where all fragments within the compute group that contain only a single -atom will have a fragment ID of 0. This can be useful in cases where -the fragment IDs are used as input for other commands in LAMMPS that -treat such single atoms different. +For the *fragment/atom* style, the *single* keyword determines whether +single atoms (not bonded to another atom) are treated as one-atom +fragments or not, based on the *yes* or *no* setting. If the setting +is *no* (the default), their fragment IDs are set to 0. An aggregate is defined by combining the rules for clusters and fragments, i.e. a set of atoms, where each of it is within the cutoff @@ -65,19 +72,11 @@ distance from one or more atoms within a fragment that is part of the same cluster. This measure can be used to track molecular assemblies like micelles. -Only atoms in the compute group are clustered and assigned cluster -IDs. Atoms not in the compute group are assigned a cluster ID = 0. -For fragments, only bonds where **both** atoms of the bond are included -in the compute group are assigned to fragments, so that only fragments -are detected where **all** atoms are in the compute group. Thus atoms -may be included in the compute group, yes still have a fragment ID of 0. - -For computes *cluster/atom* and *aggregate/atom* the neighbor list needed -to compute this quantity is constructed each time the calculation is -performed (i.e. each time a snapshot of atoms is dumped). Thus it can be -inefficient to compute/dump this quantity too frequently or to have -multiple compute/dump commands, each of a *cluster/atom* or -*aggregate/atom* style. +For computes *cluster/atom* and *aggregate/atom* a neighbor list +needed to compute cluster IDs is constructed each time the compute is +invoked. Thus it can be inefficient to compute/dump this quantity too +frequently or to have multiple *cluster/atom* or *aggregate/atom* +style computes. .. note:: @@ -100,10 +99,10 @@ multiple compute/dump commands, each of a *cluster/atom* or .. note:: For the compute fragment/atom style, each fragment is identified - using the current bond topology within each fragment. This will - not account for bonds broken by the :doc:`bond_style quartic - ` command because it does not perform a full update - of the bond topology data structures within LAMMPS. + using the current bond topology. This will not account for bonds + broken by the :doc:`bond_style quartic ` command + because it does not perform a full update of the bond topology data + structures within LAMMPS. **Output info:** @@ -123,4 +122,7 @@ Related commands :doc:`compute coord/atom ` -**Default:** none +**Default:** + +The default for fragment/atom is single no. + diff --git a/doc/src/reset_mol_ids.rst b/doc/src/reset_mol_ids.rst index b56596cc8a..bf5782fd39 100644 --- a/doc/src/reset_mol_ids.rst +++ b/doc/src/reset_mol_ids.rst @@ -15,12 +15,13 @@ Syntax * group-ID = ID of group of atoms whose molecule IDs will be reset * zero or more keyword/value pairs may be appended -* keyword = *offset* or *singlezero* +* keyword = *compress* or *offset* or *single* .. parsed-literal:: - *offset* value = *Noffset* or auto - *singlezero* = assign molecule ID 0 to atoms without bonds + *compress* value = *yes* or *no* + *offset* value = *Noffset* >= -1 + *single* value = *yes* or *no* to treat single atoms (no bonds) as molecules Examples """""""" @@ -35,17 +36,22 @@ Examples Description """"""""""" -Reset molecule IDs for a group of atoms. This will create a new set -of molecule IDs that are numbered contiguously from 1 to N, if there -are N different molecules in the group. Only molecule IDs for atoms -in the specified group are reset. +Reset molecule IDs for a group of atoms based on current bond +connectivity. This will typically create a new set of molecule IDs +for atoms in the group. Only molecule IDs for atoms in the specified +group are reset; molecule IDs for atoms not in the group are not +changed. -For purposes of this operation, molecules are identified by the bond -topology of the system, not by the current molecule IDs. A molecule -is a set of atoms, each of which is bonded to one or more atoms in the -set. If an atom is not bonded to any other atom, it becomes its own -1-atom molecule. Once new molecules are identified, this command will -overwrite the current molecule ID for each atom with a new ID. +For purposes of this operation, molecules are identified by the +current bond connectivity in the system, which may or may not be +consistent with current molecule IDs. A molecule is a set of atoms, +each of which is bonded to one or more atoms in the set. Once new +molecules are identified and a molecule ID assigned to each one, this +command will update the current molecule ID for each atom in the group +with a (potentially) new ID. Note that if the group is set so as to +exclude atoms within molecules, one molecule may become several. For +example if the group excludes atoms in the midddle of a linear chain, +then each end of the chain becomes an independent molecules. This can be a useful operation to perform after running reactive molecular dynamics run with :doc:`fix bond/react `, @@ -56,34 +62,40 @@ also be useful after molecules have been deleted with the has lost molecules, e.g. via the :doc:`fix evaporate ` command. -The *offset* keyword can be used to change the range of new molecule -IDs assigned. If *Noffset* = auto (the default) and the specified group -is *all*, then new molecule IDs will be from 1 to N. If *Noffset* = auto -and the group is **not** all, then new molecule IDs will be from M+1 to -M+N, where M is the largest molecule ID for any atom **not** in the group. -This insures new molecule IDs do not collide with existing molecule -IDs that are not changed by this command. +The *compress* keyword determines how new molecule IDs are assigned. +If the setting is *no* (the default), the molecule ID of every atom in +the molecule will be set to the smallest atom ID of any atom in the +molecule. If the setting is *yes*, and there are N molecules in the +group, the new molecule IDs will be a set of N contiguous values. See +the *offset* keyword for more details. -If *Noffset* is set to a number instead of 'auto', then new molecule IDs -will be from Noffset+1 to Noffset+N. If the group is not all, it is -up to you to insure the choice of *Noffset* does not produce -collisions between existing and new molecule IDs. +The *single* keyword determines whether single atoms (not bonded to +another atom) are treated as one-atom molecules or not, based on the +*yes* or *no* setting. If the setting is *no* (the default), their +molecule IDs are set to 0. This setting can be important if the new +molecule IDs will be used as input to other commands such as +:doc:`compute chunk/atom molecule ` or :doc:`fix +rigid molecule `. -The *singlezero* keyword turns on a special treatment for single atoms -without bonds. Instead of assigning each atom a different molecule ID -those atoms will all be given the molecule ID 0. This can be useful -when the molecule ID is used to group atoms where atoms with a molecule -ID of 0 will be considered as individual atoms; for example when using -:doc:`fix rigid ` with the *molecule* option. +The *offset* keyword is only used if the *compress* setting is *yes*. +Its default value is *Noffset* = -1. In that case, if the specified +group is *all*, then the new compressed molecule IDs will range from 1 +to N. If the specified group is not *all* and the largest molecule ID +in the non-group atoms is M, then the new compressed molecule IDs will +range from M+1 to M+N, so as to not collide with existing molecule +IDs. If an *Noffset* >= 0 is specified, then the new compressed +molecule IDs will range from *Noffset*+1 to *Noffset*+N. If the group +is not *all* it is up to you to insure there are no collisions with +the molecule IDs of non-group atoms. .. note:: The same as explained for the :doc:`compute fragment/atom ` command, molecules are identified using the - current bond topology within each fragment. This will not account - for bonds broken by the :doc:`bond_style quartic ` - command because it does not perform a full update of the bond - topology data structures within LAMMPS. + current bond topology. This will not account for bonds broken by + the :doc:`bond_style quartic ` command because it + does not perform a full update of the bond topology data structures + within LAMMPS. Restrictions """""""""""" @@ -100,5 +112,5 @@ Related commands **Default:** -The default keyword setting is offset = 0. - +The default keyword settings are compress = no, single = no, and +offset = -1. diff --git a/src/compute_fragment_atom.cpp b/src/compute_fragment_atom.cpp index e71593ce3d..aad427a7c8 100644 --- a/src/compute_fragment_atom.cpp +++ b/src/compute_fragment_atom.cpp @@ -38,13 +38,6 @@ ComputeFragmentAtom::ComputeFragmentAtom(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), fragmentID(NULL) { - singleflag = 0; - - if ((narg < 3) || (narg > 4)) - error->all(FLERR,"Illegal compute fragment/atom command"); - if ((narg == 4) && (strcmp(arg[3],"singlezero") == 0)) - singleflag = 1; - if (atom->avec->bonds_allow == 0) error->all(FLERR,"Compute fragment/atom used when bonds are not allowed"); @@ -52,6 +45,21 @@ ComputeFragmentAtom::ComputeFragmentAtom(LAMMPS *lmp, int narg, char **arg) : size_peratom_cols = 0; comm_forward = 1; + // process optional args + + singleflag = 0; + + int iarg = 3; + while (iarg < narg) { + if (strcmp(arg[iarg],"single") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal compute fragment/atom command"); + if (strcmp(arg[iarg+1],"yes") == 0) singleflag = 1; + else if (strcmp(arg[iarg+1],"no") == 0) singleflag = 0; + else error->all(FLERR,"Illegal compute fragment/atom command"); + iarg += 2; + } else error->all(FLERR,"Illegal compute fragment/atom command"); + } + nmax = 0; stack = NULL; clist = NULL; @@ -119,7 +127,6 @@ void ComputeFragmentAtom::compute_peratom() // owned + ghost atoms start with fragmentID = atomID // atoms not in group have fragmentID = 0 - // if singleflag is set atoms without bonds have fragmentID 0 as well. tagint *tag = atom->tag; int *mask = atom->mask; @@ -137,6 +144,7 @@ void ComputeFragmentAtom::compute_peratom() // acquire fragmentIDs of ghost atoms // loop over clusters of atoms, which include ghost atoms // set fragmentIDs for each cluster to min framentID in the clusters + // if singleflag = 0 atoms without bonds are assigned fragmentID = 0 // iterate until no changes to ghost atom fragmentIDs commflag = 1; @@ -159,17 +167,15 @@ void ComputeFragmentAtom::compute_peratom() for (i = 0; i < nlocal; i++) { // skip atom I if not in group or already marked - // also skip and set fragment ID to zero if singleflag is set - // and the atom is an isolated atom without bonds + // if singleflag = 0 and atom has no bond partners, fragID = 0 and done if (!(mask[i] & groupbit)) continue; if (markflag[i]) continue; - if (singleflag && (nspecial[i][0] == 0)) { - fragmentID[i] = 0; - markflag[i] = 1; - continue; + if (!singleflag && (nspecial[i][0] == 0)) { + fragmentID[i] = 0.0; + continue; } - + // find one cluster of bond-connected atoms // ncluster = # of owned and ghost atoms in cluster // clist = vector of local indices of the ncluster atoms diff --git a/src/reset_mol_ids.cpp b/src/reset_mol_ids.cpp index 69c2939282..3ab28d173a 100644 --- a/src/reset_mol_ids.cpp +++ b/src/reset_mol_ids.cpp @@ -25,7 +25,6 @@ #include "modify.h" #include "compute_fragment_atom.h" #include "compute_chunk_atom.h" -#include "compute_reduce.h" #include "utils.h" #include "error.h" #include "fmt/format.h" @@ -54,23 +53,29 @@ void ResetMolIDs::command(int narg, char **arg) if (igroup == -1) error->all(FLERR,"Could not find reset_mol_ids group ID"); int groupbit = group->bitmask[igroup]; - tagint offset = -1; + int compressflag = 0; int singleflag = 0; + tagint offset = -1; int iarg = 1; while (iarg < narg) { - if (strcmp(arg[iarg],"offset") == 0) { + if (strcmp(arg[iarg],"compress") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal reset_mol_ids command"); - if (strcmp(arg[iarg+1],"auto") == 0) { - offset = -1; - } else { - offset = utils::tnumeric(FLERR,arg[iarg+1],1,lmp); - if (offset < 0) error->all(FLERR,"Illegal reset_mol_ids command"); - } + if (strcmp(arg[iarg+1],"yes") == 0) compressflag = 1; + else if (strcmp(arg[iarg+1],"no") == 0) compressflag = 0; + else error->all(FLERR,"Illegal reset_mol_ids command"); + iarg += 2; + } else if (strcmp(arg[iarg],"single") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal reset_mol_ids command"); + if (strcmp(arg[iarg+1],"yes") == 0) singleflag = 1; + else if (strcmp(arg[iarg+1],"no") == 0) singleflag = 0; + else error->all(FLERR,"Illegal reset_mol_ids command"); + iarg += 2; + } else if (strcmp(arg[iarg],"offset") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal reset_mol_ids command"); + offset = utils::tnumeric(FLERR,arg[iarg+1],1,lmp); + if (offset < -1) error->all(FLERR,"Illegal reset_mol_ids command"); iarg += 2; - } else if (strcmp(arg[iarg],"singlezero") == 0) { - singleflag = 1; - ++iarg; } else error->all(FLERR,"Illegal reset_mol_ids command"); } @@ -86,17 +91,14 @@ void ResetMolIDs::command(int narg, char **arg) const std::string idfrag = "reset_mol_ids_FRAGMENT_ATOM"; if (singleflag) - modify->add_compute(fmt::format("{} {} fragment/atom singlezero",idfrag,arg[0])); + modify->add_compute(fmt::format("{} {} fragment/atom single yes",idfrag,arg[0])); else - modify->add_compute(fmt::format("{} {} fragment/atom",idfrag,arg[0])); - - const std::string idmin = "reset_mol_ids_COMPUTE_MINFRAG"; - if (singleflag) - modify->add_compute(fmt::format("{} {} reduce min c_{}",idmin,arg[0],idfrag)); + modify->add_compute(fmt::format("{} {} fragment/atom single no",idfrag,arg[0])); const std::string idchunk = "reset_mol_ids_CHUNK_ATOM"; - modify->add_compute(fmt::format("{} {} chunk/atom molecule compress yes", - idchunk,arg[0])); + if (compressflag) + modify->add_compute(fmt::format("{} {} chunk/atom molecule compress yes", + idchunk,arg[0])); // initialize system since comm->borders() will be invoked @@ -116,24 +118,17 @@ void ResetMolIDs::command(int narg, char **arg) // invoke peratom method of compute fragment/atom // walks bond connectivity and assigns each atom a fragment ID - + // if singleflag = 0, atoms w/out bonds will be assigned fragID = 0 + int icompute = modify->find_compute(idfrag); ComputeFragmentAtom *cfa = (ComputeFragmentAtom *) modify->compute[icompute]; cfa->compute_peratom(); - double *ids = cfa->vector_atom; + double *fragIDs = cfa->vector_atom; - // when the lowest fragment ID is 0, in case the singlezero option is used - // we need to remove adjust the chunk ID, so that the resuling molecule ID is correct. - - int adjid = 0; - if (singleflag) { - icompute = modify->find_compute(idmin); - ComputeReduce *crd = (ComputeReduce *) modify->compute[icompute]; - if (crd->compute_scalar() == 0.0) adjid = 1; - } - - // copy fragmentID to molecule ID - // only for atoms in the group + for (int i = 0; i < atom->nlocal; i++) + printf("FragID %d %g\n",atom->tag[i],fragIDs[i]); + + // copy fragID to molecule ID for atoms in group tagint *molecule = atom->molecule; int *mask = atom->mask; @@ -141,58 +136,84 @@ void ResetMolIDs::command(int narg, char **arg) for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) - molecule[i] = static_cast (ids[i]); + molecule[i] = static_cast (fragIDs[i]); - // invoke peratom method of compute chunk/atom - // compress new molecule IDs to be contiguous 1 to Nmol - // NOTE: use of compute chunk/atom limits # of molecules to a 32-bit int + // if compressflag = 0, done + // set nchunk = -1 since cannot easily determine # of new molecule IDs - icompute = modify->find_compute(idchunk); - ComputeChunkAtom *cca = (ComputeChunkAtom *) modify->compute[icompute]; - cca->compute_peratom(); - ids = cca->vector_atom; - int nchunk = cca->nchunk; + int nchunk = -1; + + // if compressflag = 1, invoke peratom method of compute chunk/atom + // will compress new molecule IDs to be contiguous 1 to Nmol + // NOTE: use of compute chunk/atom limits Nmol to a 32-bit int - // if offset = -1 (i.e. "auto") and group != all, - // reset offset to largest molID of non-group atoms - // otherwise set to 0. + if (compressflag) { + icompute = modify->find_compute(idchunk); + ComputeChunkAtom *cca = (ComputeChunkAtom *) modify->compute[icompute]; + cca->compute_peratom(); + double *chunkIDs = cca->vector_atom; + nchunk = cca->nchunk; + + for (int i = 0; i < atom->nlocal; i++) + printf("ChunkID %d %g\n",atom->tag[i],chunkIDs[i]); - if (offset == -1) { - if (groupbit != 1) { - tagint mymol = 0; - for (int i = 0; i < nlocal; i++) - if (!(mask[i] & groupbit)) - mymol = MAX(mymol,molecule[i]); - MPI_Allreduce(&mymol,&offset,1,MPI_LMP_TAGINT,MPI_MAX,world); - } else offset = 0; - } + // if singleflag = 0, check if any single (no-bond) atoms exist + // if so, they have fragID = 0, and compression just set their chunkID = 1 + // singleexist = 0/1 if no/yes single atoms exist with chunkID = 1 + + int singleexist = 0; + if (!singleflag) { + int mysingle = 0; + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + if (fragIDs[i] == 0.0) mysingle = 1; + MPI_Allreduce(&mysingle,&singleexist,1,MPI_INT,MPI_MAX,world); + if (singleexist) nchunk--; + } - // copy chunkID to molecule ID + offset - // only for atoms in the group + // if offset < 0 (default), reset it + // if group = all, offset = 0 + // else offset = largest molID of non-group atoms + + if (offset < 0) { + if (groupbit != 1) { + tagint mymol = 0; + for (int i = 0; i < nlocal; i++) + if (!(mask[i] & groupbit)) + mymol = MAX(mymol,molecule[i]); + MPI_Allreduce(&mymol,&offset,1,MPI_LMP_TAGINT,MPI_MAX,world); + } else offset = 0; + } - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - tagint newid = static_cast(ids[i]); - if (singleflag && adjid) { - if (newid == 1) newid = 0; - else newid += offset - 1; - molecule[i] = newid; - } else molecule[i] = newid + offset; + // reset molecule ID for all atoms in group + // newID = chunkID + offset + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + tagint newid = static_cast(chunkIDs[i]); + if (singleexist) { + if (newid == 1) newid = 0; + else newid += offset - 1; + } else newid += offset; + molecule[i] = newid; + } } } - + // clean up - modify->delete_compute(idchunk); modify->delete_compute(idfrag); - if (singleflag) modify->delete_compute(idmin); + if (compressflag) modify->delete_compute(idchunk); // total time MPI_Barrier(world); if (comm->me == 0) { - utils::logmesg(lmp,fmt::format(" number of new molecule IDs = {}\n",nchunk)); + if (nchunk < 0) + utils::logmesg(lmp,fmt::format(" number of new molecule IDs = unknown\n")); + else + utils::logmesg(lmp,fmt::format(" number of new molecule IDs = {}\n",nchunk)); utils::logmesg(lmp,fmt::format(" reset_mol_ids CPU = {:.3f} seconds\n", MPI_Wtime()-time1)); } From 4a1acffde989f37b50b010c75589980bc67d98c4 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 9 Jul 2020 19:36:58 -0400 Subject: [PATCH 21/29] whitespace fixes --- src/compute_fragment_atom.cpp | 22 +++++++++--------- src/reset_atom_ids.cpp | 10 ++++---- src/reset_mol_ids.cpp | 44 +++++++++++++++++------------------ 3 files changed, 38 insertions(+), 38 deletions(-) diff --git a/src/compute_fragment_atom.cpp b/src/compute_fragment_atom.cpp index aad427a7c8..664e2fa9ec 100644 --- a/src/compute_fragment_atom.cpp +++ b/src/compute_fragment_atom.cpp @@ -48,7 +48,7 @@ ComputeFragmentAtom::ComputeFragmentAtom(LAMMPS *lmp, int narg, char **arg) : // process optional args singleflag = 0; - + int iarg = 3; while (iarg < narg) { if (strcmp(arg[iarg],"single") == 0) { @@ -134,7 +134,7 @@ void ComputeFragmentAtom::compute_peratom() int **nspecial = atom->nspecial; int nlocal = atom->nlocal; int nall = nlocal + atom->nghost; - + for (i = 0; i < nall; i++) { if (mask[i] & groupbit) fragmentID[i] = tag[i]; else fragmentID[i] = 0; @@ -150,7 +150,7 @@ void ComputeFragmentAtom::compute_peratom() commflag = 1; int iteration = 0; - + while (1) { iteration++; @@ -163,24 +163,24 @@ void ComputeFragmentAtom::compute_peratom() // loop over all owned atoms // each unmarked atom starts a cluster search - + for (i = 0; i < nlocal; i++) { // skip atom I if not in group or already marked // if singleflag = 0 and atom has no bond partners, fragID = 0 and done - + if (!(mask[i] & groupbit)) continue; if (markflag[i]) continue; if (!singleflag && (nspecial[i][0] == 0)) { - fragmentID[i] = 0.0; - continue; + fragmentID[i] = 0.0; + continue; } - + // find one cluster of bond-connected atoms // ncluster = # of owned and ghost atoms in cluster // clist = vector of local indices of the ncluster atoms // stack is used to walk the bond topology - + ncluster = nstack = 0; stack[nstack++] = i; @@ -218,14 +218,14 @@ void ComputeFragmentAtom::compute_peratom() // set fragmentID = newID for all atoms in cluster, including ghost atoms // not done with iterations if change the fragmentID of a ghost atom - + for (m = 0; m < ncluster; m++) { j = clist[m]; if (j >= nlocal && fragmentID[j] != newID) done = 0; fragmentID[j] = newID; } } - + // stop if all procs are done MPI_Allreduce(&done,&alldone,1,MPI_INT,MPI_MIN,world); diff --git a/src/reset_atom_ids.cpp b/src/reset_atom_ids.cpp index 7e493ca28e..5ab7400a52 100644 --- a/src/reset_atom_ids.cpp +++ b/src/reset_atom_ids.cpp @@ -406,8 +406,8 @@ void ResetIDs::sort() char *buf; int nreturn = comm->rendezvous(1,nlocal,(char *) atombuf,sizeof(AtomRvous), - 0,proclist, - sort_bins,0,buf,sizeof(IDRvous),(void *) this); + 0,proclist, + sort_bins,0,buf,sizeof(IDRvous),(void *) this); IDRvous *outbuf = (IDRvous *) buf; memory->destroy(proclist); @@ -432,8 +432,8 @@ void ResetIDs::sort() ------------------------------------------------------------------------- */ int ResetIDs::sort_bins(int n, char *inbuf, - int &flag, int *&proclist, char *&outbuf, - void *ptr) + int &flag, int *&proclist, char *&outbuf, + void *ptr) { int i,ibin,index; @@ -467,7 +467,7 @@ int ResetIDs::sort_bins(int n, char *inbuf, for (i = 0; i < n; i++) { if (in[i].ibin < binlo || in[i].ibin >= binhi) { //printf("BAD me %d i %d %d ibin %d binlohi %d %d\n", - // rptr->comm->me,i,n,in[i].ibin,binlo,binhi); + // rptr->comm->me,i,n,in[i].ibin,binlo,binhi); error->one(FLERR,"Bad spatial bin assignment in reset_atom_ids sort"); } ibin = in[i].ibin - binlo; diff --git a/src/reset_mol_ids.cpp b/src/reset_mol_ids.cpp index 3ab28d173a..4f648734bd 100644 --- a/src/reset_mol_ids.cpp +++ b/src/reset_mol_ids.cpp @@ -98,7 +98,7 @@ void ResetMolIDs::command(int narg, char **arg) const std::string idchunk = "reset_mol_ids_CHUNK_ATOM"; if (compressflag) modify->add_compute(fmt::format("{} {} chunk/atom molecule compress yes", - idchunk,arg[0])); + idchunk,arg[0])); // initialize system since comm->borders() will be invoked @@ -119,7 +119,7 @@ void ResetMolIDs::command(int narg, char **arg) // invoke peratom method of compute fragment/atom // walks bond connectivity and assigns each atom a fragment ID // if singleflag = 0, atoms w/out bonds will be assigned fragID = 0 - + int icompute = modify->find_compute(idfrag); ComputeFragmentAtom *cfa = (ComputeFragmentAtom *) modify->compute[icompute]; cfa->compute_peratom(); @@ -127,7 +127,7 @@ void ResetMolIDs::command(int narg, char **arg) for (int i = 0; i < atom->nlocal; i++) printf("FragID %d %g\n",atom->tag[i],fragIDs[i]); - + // copy fragID to molecule ID for atoms in group tagint *molecule = atom->molecule; @@ -142,7 +142,7 @@ void ResetMolIDs::command(int narg, char **arg) // set nchunk = -1 since cannot easily determine # of new molecule IDs int nchunk = -1; - + // if compressflag = 1, invoke peratom method of compute chunk/atom // will compress new molecule IDs to be contiguous 1 to Nmol // NOTE: use of compute chunk/atom limits Nmol to a 32-bit int @@ -153,20 +153,20 @@ void ResetMolIDs::command(int narg, char **arg) cca->compute_peratom(); double *chunkIDs = cca->vector_atom; nchunk = cca->nchunk; - + for (int i = 0; i < atom->nlocal; i++) printf("ChunkID %d %g\n",atom->tag[i],chunkIDs[i]); // if singleflag = 0, check if any single (no-bond) atoms exist // if so, they have fragID = 0, and compression just set their chunkID = 1 // singleexist = 0/1 if no/yes single atoms exist with chunkID = 1 - + int singleexist = 0; if (!singleflag) { int mysingle = 0; for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) - if (fragIDs[i] == 0.0) mysingle = 1; + if (mask[i] & groupbit) + if (fragIDs[i] == 0.0) mysingle = 1; MPI_Allreduce(&mysingle,&singleexist,1,MPI_INT,MPI_MAX,world); if (singleexist) nchunk--; } @@ -174,32 +174,32 @@ void ResetMolIDs::command(int narg, char **arg) // if offset < 0 (default), reset it // if group = all, offset = 0 // else offset = largest molID of non-group atoms - + if (offset < 0) { if (groupbit != 1) { - tagint mymol = 0; - for (int i = 0; i < nlocal; i++) - if (!(mask[i] & groupbit)) - mymol = MAX(mymol,molecule[i]); - MPI_Allreduce(&mymol,&offset,1,MPI_LMP_TAGINT,MPI_MAX,world); + tagint mymol = 0; + for (int i = 0; i < nlocal; i++) + if (!(mask[i] & groupbit)) + mymol = MAX(mymol,molecule[i]); + MPI_Allreduce(&mymol,&offset,1,MPI_LMP_TAGINT,MPI_MAX,world); } else offset = 0; } // reset molecule ID for all atoms in group // newID = chunkID + offset - + for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - tagint newid = static_cast(chunkIDs[i]); - if (singleexist) { - if (newid == 1) newid = 0; - else newid += offset - 1; - } else newid += offset; - molecule[i] = newid; + tagint newid = static_cast(chunkIDs[i]); + if (singleexist) { + if (newid == 1) newid = 0; + else newid += offset - 1; + } else newid += offset; + molecule[i] = newid; } } } - + // clean up modify->delete_compute(idfrag); From e0e24799c2ccb4b55e7b5aa01a49276ece2633ad Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 9 Jul 2020 19:52:49 -0400 Subject: [PATCH 22/29] simplify formulations a little bit. update example command lines --- doc/src/reset_mol_ids.rst | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/doc/src/reset_mol_ids.rst b/doc/src/reset_mol_ids.rst index bf5782fd39..31fe582648 100644 --- a/doc/src/reset_mol_ids.rst +++ b/doc/src/reset_mol_ids.rst @@ -29,7 +29,7 @@ Examples .. code-block:: LAMMPS reset_mol_ids all - reset_mol_ids all offset 10 singlezero + reset_mol_ids all offset 10 single yes compress yes reset_mol_ids solvent offset 1000 reset_mol_ids solvent offset auto @@ -48,10 +48,11 @@ consistent with current molecule IDs. A molecule is a set of atoms, each of which is bonded to one or more atoms in the set. Once new molecules are identified and a molecule ID assigned to each one, this command will update the current molecule ID for each atom in the group -with a (potentially) new ID. Note that if the group is set so as to -exclude atoms within molecules, one molecule may become several. For -example if the group excludes atoms in the midddle of a linear chain, -then each end of the chain becomes an independent molecules. +with a (potentially) new ID. Note that if the group excludes atoms +within molecules, one molecule may become two or more. For +example if the group excludes atoms in the middle of a linear chain, +then each end of the chain is considered an independent molecule +and will be assigned a different molecule ID. This can be a useful operation to perform after running reactive molecular dynamics run with :doc:`fix bond/react `, @@ -67,7 +68,7 @@ If the setting is *no* (the default), the molecule ID of every atom in the molecule will be set to the smallest atom ID of any atom in the molecule. If the setting is *yes*, and there are N molecules in the group, the new molecule IDs will be a set of N contiguous values. See -the *offset* keyword for more details. +the *offset* keyword for details on the selecting the range of these values. The *single* keyword determines whether single atoms (not bonded to another atom) are treated as one-atom molecules or not, based on the @@ -81,12 +82,11 @@ The *offset* keyword is only used if the *compress* setting is *yes*. Its default value is *Noffset* = -1. In that case, if the specified group is *all*, then the new compressed molecule IDs will range from 1 to N. If the specified group is not *all* and the largest molecule ID -in the non-group atoms is M, then the new compressed molecule IDs will -range from M+1 to M+N, so as to not collide with existing molecule +of atoms outside that group is M, then the new compressed molecule IDs will +range from M+1 to M+N, to avoid collision with existing molecule IDs. If an *Noffset* >= 0 is specified, then the new compressed -molecule IDs will range from *Noffset*+1 to *Noffset*+N. If the group -is not *all* it is up to you to insure there are no collisions with -the molecule IDs of non-group atoms. +molecule IDs will range from *Noffset*\ +1 to *Noffset*\ +N. If the group +is not *all* there may be collisions with the molecule IDs of other atoms. .. note:: From 9c97ca11fee1e8d27deccd431b98e2df5d8e394c Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 9 Jul 2020 20:46:31 -0400 Subject: [PATCH 23/29] adjust regex for removed styles to correctly handle command styles --- doc/utils/check-styles.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/utils/check-styles.py b/doc/utils/check-styles.py index a5bf31801c..55933cc6c4 100755 --- a/doc/utils/check-styles.py +++ b/doc/utils/check-styles.py @@ -60,7 +60,7 @@ kokkos = re.compile("(.+)/kk$") kokkos_skip = re.compile("(.+)/kk/(host|device)$") omp = re.compile("(.+)/omp$") opt = re.compile("(.+)/opt$") -removed = re.compile("(.+)Deprecated$") +removed = re.compile("(.*)Deprecated$") def register_style(list,style,info): if style in list.keys(): From 9ec77585ea95f6b0434e53c2d666a691801ed50f Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 9 Jul 2020 20:46:45 -0400 Subject: [PATCH 24/29] make compress yes the default --- doc/src/reset_mol_ids.rst | 50 +++++++++++++++++++-------------------- src/reset_mol_ids.cpp | 2 +- 2 files changed, 26 insertions(+), 26 deletions(-) diff --git a/doc/src/reset_mol_ids.rst b/doc/src/reset_mol_ids.rst index 31fe582648..0d6063b3ef 100644 --- a/doc/src/reset_mol_ids.rst +++ b/doc/src/reset_mol_ids.rst @@ -6,9 +6,6 @@ reset_mol_ids command Syntax """""" -Syntax -"""""" - .. parsed-literal:: reset_mol_ids group-ID keyword value ... @@ -29,9 +26,9 @@ Examples .. code-block:: LAMMPS reset_mol_ids all - reset_mol_ids all offset 10 single yes compress yes - reset_mol_ids solvent offset 1000 - reset_mol_ids solvent offset auto + reset_mol_ids all offset 10 single yes + reset_mol_ids solvent compress yes offset 100 + reset_mol_ids solvent compress no Description """"""""""" @@ -42,17 +39,18 @@ for atoms in the group. Only molecule IDs for atoms in the specified group are reset; molecule IDs for atoms not in the group are not changed. -For purposes of this operation, molecules are identified by the -current bond connectivity in the system, which may or may not be -consistent with current molecule IDs. A molecule is a set of atoms, -each of which is bonded to one or more atoms in the set. Once new -molecules are identified and a molecule ID assigned to each one, this -command will update the current molecule ID for each atom in the group -with a (potentially) new ID. Note that if the group excludes atoms -within molecules, one molecule may become two or more. For -example if the group excludes atoms in the middle of a linear chain, -then each end of the chain is considered an independent molecule -and will be assigned a different molecule ID. +For purposes of this operation, molecules are identified by the current +bond connectivity in the system, which may or may not be consistent with +the current molecule IDs. A molecule in this context is a set of atoms +connected to each other with explicit bonds. The specific algorithm +used is the one of :doc:`compute fragment/atom ` +Once the molecules are identified and a new molecule ID computed for +each, this command will update the current molecule ID for all atoms in +the group with the new molecule ID. Note that if the group excludes +atoms within molecules, one (physical) molecule may become two or more +(logical) molecules. For example if the group excludes atoms in the +middle of a linear chain, then each end of the chain is considered an +independent molecule and will be assigned a different molecule ID. This can be a useful operation to perform after running reactive molecular dynamics run with :doc:`fix bond/react `, @@ -63,12 +61,12 @@ also be useful after molecules have been deleted with the has lost molecules, e.g. via the :doc:`fix evaporate ` command. -The *compress* keyword determines how new molecule IDs are assigned. -If the setting is *no* (the default), the molecule ID of every atom in -the molecule will be set to the smallest atom ID of any atom in the -molecule. If the setting is *yes*, and there are N molecules in the +The *compress* keyword determines how new molecule IDs are computed. If +the setting is *yes* (the default) and there are N molecules in the group, the new molecule IDs will be a set of N contiguous values. See -the *offset* keyword for details on the selecting the range of these values. +the *offset* keyword for details on selecting the range of these values. +If the setting is *no*, the molecule ID of every atom in the molecule +will be set to the smallest atom ID of any atom in the molecule. The *single* keyword determines whether single atoms (not bonded to another atom) are treated as one-atom molecules or not, based on the @@ -108,9 +106,11 @@ Related commands :doc:`fix bond/create `, :doc:`fix bond/break `, :doc:`fix evaporate `, -:doc:`delete_atoms ` +:doc:`delete_atoms `, +:doc:`compute fragment/atom ` -**Default:** +Default +""""""" -The default keyword settings are compress = no, single = no, and +The default keyword settings are compress = yes, single = no, and offset = -1. diff --git a/src/reset_mol_ids.cpp b/src/reset_mol_ids.cpp index 4f648734bd..6df76501e4 100644 --- a/src/reset_mol_ids.cpp +++ b/src/reset_mol_ids.cpp @@ -53,7 +53,7 @@ void ResetMolIDs::command(int narg, char **arg) if (igroup == -1) error->all(FLERR,"Could not find reset_mol_ids group ID"); int groupbit = group->bitmask[igroup]; - int compressflag = 0; + int compressflag = 1; int singleflag = 0; tagint offset = -1; From 49780480a80d0ebdf062ac982ed2d983eb01beb5 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 9 Jul 2020 20:52:57 -0400 Subject: [PATCH 25/29] count total number of styles including aliases, suffixes, and undocumented --- doc/utils/check-styles.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/doc/utils/check-styles.py b/doc/utils/check-styles.py index 55933cc6c4..40e6109f5b 100755 --- a/doc/utils/check-styles.py +++ b/doc/utils/check-styles.py @@ -52,6 +52,7 @@ minimize = {} pair = {} reader = {} region = {} +total = 0 upper = re.compile("[A-Z]+") gpu = re.compile("(.+)/gpu$") @@ -119,6 +120,7 @@ for h in headers: # skip over internal styles w/o explicit documentation style = m[1] + total += 1 if upper.match(style): continue @@ -197,12 +199,13 @@ print("""Parsed style names w/o suffixes from C++ tree in %s: Fix styles: %3d Improper styles: %3d Integrate styles: %3d Kspace styles: %3d Minimize styles: %3d Pair styles: %3d - Reader styles: %3d Region styles: %3d""" \ + Reader styles: %3d Region styles: %3d +---------------------------------------------------- +Total number of styles (including suffixes): %d""" \ % (src, len(angle), len(atom), len(body), len(bond), \ len(command), len(compute), len(dihedral), len(dump), \ len(fix), len(improper), len(integrate), len(kspace), \ - len(minimize), len(pair), len(reader), len(region))) - + len(minimize), len(pair), len(reader), len(region), total)) counter = 0 From bade009b6c3084c5c48440d5a4c7791602a7d9e1 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 9 Jul 2020 20:58:39 -0400 Subject: [PATCH 26/29] remove debug code --- src/reset_mol_ids.cpp | 6 ------ 1 file changed, 6 deletions(-) diff --git a/src/reset_mol_ids.cpp b/src/reset_mol_ids.cpp index 6df76501e4..e3cc877548 100644 --- a/src/reset_mol_ids.cpp +++ b/src/reset_mol_ids.cpp @@ -125,9 +125,6 @@ void ResetMolIDs::command(int narg, char **arg) cfa->compute_peratom(); double *fragIDs = cfa->vector_atom; - for (int i = 0; i < atom->nlocal; i++) - printf("FragID %d %g\n",atom->tag[i],fragIDs[i]); - // copy fragID to molecule ID for atoms in group tagint *molecule = atom->molecule; @@ -154,9 +151,6 @@ void ResetMolIDs::command(int narg, char **arg) double *chunkIDs = cca->vector_atom; nchunk = cca->nchunk; - for (int i = 0; i < atom->nlocal; i++) - printf("ChunkID %d %g\n",atom->tag[i],chunkIDs[i]); - // if singleflag = 0, check if any single (no-bond) atoms exist // if so, they have fragID = 0, and compression just set their chunkID = 1 // singleexist = 0/1 if no/yes single atoms exist with chunkID = 1 From d15264a668b3ac325672668177cff0b1dc0edf40 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 9 Jul 2020 21:18:58 -0400 Subject: [PATCH 27/29] update and expand tester code for reset_mol_ids changes --- unittest/commands/test_reset_ids.cpp | 117 +++++++++++++++++---------- 1 file changed, 75 insertions(+), 42 deletions(-) diff --git a/unittest/commands/test_reset_ids.cpp b/unittest/commands/test_reset_ids.cpp index 8929c35dab..81d4b01ab5 100644 --- a/unittest/commands/test_reset_ids.cpp +++ b/unittest/commands/test_reset_ids.cpp @@ -314,7 +314,7 @@ TEST_F(ResetIDsTest, DeleteAdd) lmp->input->one("group twowater molecule 4:6:2"); lmp->input->one("group nowater subtract all allwater"); lmp->input->one("delete_atoms group twowater compress no bond yes mol yes"); - lmp->input->one("reset_mol_ids allwater offset auto"); + lmp->input->one("reset_mol_ids allwater"); if (!verbose) ::testing::internal::GetCapturedStdout(); ASSERT_EQ(lmp->atom->natoms, 23); ASSERT_EQ(lmp->atom->map_tag_max, 26); @@ -384,59 +384,71 @@ TEST_F(ResetIDsTest, DeleteAdd) lmp->input->one("reset_mol_ids nowater offset 1"); if (!verbose) ::testing::internal::GetCapturedStdout(); - EXPECT_EQ(molid[GETIDX(1)], 2); - EXPECT_EQ(molid[GETIDX(2)], 2); - EXPECT_EQ(molid[GETIDX(3)], 2); - EXPECT_EQ(molid[GETIDX(4)], 2); - EXPECT_EQ(molid[GETIDX(5)], 2); - EXPECT_EQ(molid[GETIDX(6)], 2); - EXPECT_EQ(molid[GETIDX(7)], 2); - EXPECT_EQ(molid[GETIDX(8)], 2); - EXPECT_EQ(molid[GETIDX(9)], 2); - EXPECT_EQ(molid[GETIDX(10)], 2); - EXPECT_EQ(molid[GETIDX(11)], 2); - EXPECT_EQ(molid[GETIDX(12)], 2); - EXPECT_EQ(molid[GETIDX(13)], 3); - EXPECT_EQ(molid[GETIDX(14)], 3); - EXPECT_EQ(molid[GETIDX(15)], 3); - EXPECT_EQ(molid[GETIDX(16)], 2); - EXPECT_EQ(molid[GETIDX(17)], 2); - EXPECT_EQ(molid[GETIDX(18)], 2); - EXPECT_EQ(molid[GETIDX(19)], 2); - EXPECT_EQ(molid[GETIDX(20)], 2); - EXPECT_EQ(molid[GETIDX(21)], 4); - EXPECT_EQ(molid[GETIDX(22)], 4); - EXPECT_EQ(molid[GETIDX(23)], 4); + ASSERT_EQ(molid[GETIDX(1)], 2); + ASSERT_EQ(molid[GETIDX(2)], 2); + ASSERT_EQ(molid[GETIDX(3)], 2); + ASSERT_EQ(molid[GETIDX(4)], 2); + ASSERT_EQ(molid[GETIDX(5)], 2); + ASSERT_EQ(molid[GETIDX(6)], 2); + ASSERT_EQ(molid[GETIDX(7)], 2); + ASSERT_EQ(molid[GETIDX(8)], 2); + ASSERT_EQ(molid[GETIDX(9)], 2); + ASSERT_EQ(molid[GETIDX(10)], 2); + ASSERT_EQ(molid[GETIDX(11)], 2); + ASSERT_EQ(molid[GETIDX(12)], 2); + ASSERT_EQ(molid[GETIDX(13)], 3); + ASSERT_EQ(molid[GETIDX(14)], 3); + ASSERT_EQ(molid[GETIDX(15)], 3); + ASSERT_EQ(molid[GETIDX(16)], 2); + ASSERT_EQ(molid[GETIDX(17)], 2); + ASSERT_EQ(molid[GETIDX(18)], 2); + ASSERT_EQ(molid[GETIDX(19)], 2); + ASSERT_EQ(molid[GETIDX(20)], 2); + ASSERT_EQ(molid[GETIDX(21)], 4); + ASSERT_EQ(molid[GETIDX(22)], 4); + ASSERT_EQ(molid[GETIDX(23)], 4); if (!verbose) ::testing::internal::CaptureStdout(); lmp->input->one("create_atoms 1 single 0.0 0.0 0.0"); lmp->input->one("create_atoms 2 single 1.0 0.0 0.0"); lmp->input->one("create_atoms 3 single 2.0 0.0 0.0"); lmp->input->one("create_atoms 4 single 3.0 0.0 0.0"); - lmp->input->one("reset_mol_ids all"); + lmp->input->one("reset_mol_ids all single yes"); if (!verbose) ::testing::internal::GetCapturedStdout(); - EXPECT_EQ(lmp->atom->natoms, 27); - EXPECT_EQ(lmp->atom->map_tag_max, 27); + ASSERT_EQ(lmp->atom->natoms, 27); + ASSERT_EQ(lmp->atom->map_tag_max, 27); - EXPECT_EQ(molid[GETIDX(21)], 3); - EXPECT_EQ(molid[GETIDX(22)], 3); - EXPECT_EQ(molid[GETIDX(23)], 3); - EXPECT_EQ(molid[GETIDX(24)], 4); - EXPECT_EQ(molid[GETIDX(25)], 5); - EXPECT_EQ(molid[GETIDX(26)], 6); - EXPECT_EQ(molid[GETIDX(27)], 7); + ASSERT_EQ(molid[GETIDX(21)], 3); + ASSERT_EQ(molid[GETIDX(22)], 3); + ASSERT_EQ(molid[GETIDX(23)], 3); + ASSERT_EQ(molid[GETIDX(24)], 4); + ASSERT_EQ(molid[GETIDX(25)], 5); + ASSERT_EQ(molid[GETIDX(26)], 6); + ASSERT_EQ(molid[GETIDX(27)], 7); if (!verbose) ::testing::internal::CaptureStdout(); - lmp->input->one("reset_mol_ids all singlezero"); + lmp->input->one("reset_mol_ids all single no"); if (!verbose) ::testing::internal::GetCapturedStdout(); - EXPECT_EQ(molid[GETIDX(21)], 3); - EXPECT_EQ(molid[GETIDX(22)], 3); - EXPECT_EQ(molid[GETIDX(23)], 3); - EXPECT_EQ(molid[GETIDX(24)], 0); - EXPECT_EQ(molid[GETIDX(25)], 0); - EXPECT_EQ(molid[GETIDX(26)], 0); - EXPECT_EQ(molid[GETIDX(27)], 0); + ASSERT_EQ(molid[GETIDX(21)], 3); + ASSERT_EQ(molid[GETIDX(22)], 3); + ASSERT_EQ(molid[GETIDX(23)], 3); + ASSERT_EQ(molid[GETIDX(24)], 0); + ASSERT_EQ(molid[GETIDX(25)], 0); + ASSERT_EQ(molid[GETIDX(26)], 0); + ASSERT_EQ(molid[GETIDX(27)], 0); + + if (!verbose) ::testing::internal::CaptureStdout(); + lmp->input->one("reset_mol_ids all compress no single yes"); + if (!verbose) ::testing::internal::GetCapturedStdout(); + + ASSERT_EQ(molid[GETIDX(21)], 21); + ASSERT_EQ(molid[GETIDX(22)], 21); + ASSERT_EQ(molid[GETIDX(23)], 21); + ASSERT_EQ(molid[GETIDX(24)], 24); + ASSERT_EQ(molid[GETIDX(25)], 25); + ASSERT_EQ(molid[GETIDX(26)], 26); + ASSERT_EQ(molid[GETIDX(27)], 27); } TEST_F(ResetIDsTest, DeathTests) @@ -461,10 +473,31 @@ TEST_F(ResetIDsTest, DeathTests) mesg = ::testing::internal::GetCapturedStdout(); ASSERT_THAT(mesg, MatchesRegex(".*ERROR on proc 0: Expected integer.*")); + ::testing::internal::CaptureStdout(); + TEST_FAILURE(lmp->input->one("reset_mol_ids all compress yes single no offset xxx");); + mesg = ::testing::internal::GetCapturedStdout(); + ASSERT_THAT(mesg, MatchesRegex(".*ERROR on proc 0: Expected integer.*")); + ::testing::internal::CaptureStdout(); TEST_FAILURE(lmp->input->one("reset_mol_ids all offset");); mesg = ::testing::internal::GetCapturedStdout(); ASSERT_THAT(mesg, MatchesRegex(".*ERROR: Illegal reset_mol_ids command.*")); + ::testing::internal::CaptureStdout(); + TEST_FAILURE(lmp->input->one("reset_mol_ids all compress");); + mesg = ::testing::internal::GetCapturedStdout(); + ASSERT_THAT(mesg, MatchesRegex(".*ERROR: Illegal reset_mol_ids command.*")); + ::testing::internal::CaptureStdout(); + TEST_FAILURE(lmp->input->one("reset_mol_ids all compress xxx");); + mesg = ::testing::internal::GetCapturedStdout(); + ASSERT_THAT(mesg, MatchesRegex(".*ERROR: Illegal reset_mol_ids command.*")); + ::testing::internal::CaptureStdout(); + TEST_FAILURE(lmp->input->one("reset_mol_ids all single");); + mesg = ::testing::internal::GetCapturedStdout(); + ASSERT_THAT(mesg, MatchesRegex(".*ERROR: Illegal reset_mol_ids command.*")); + ::testing::internal::CaptureStdout(); + TEST_FAILURE(lmp->input->one("reset_mol_ids all single xxx");); + mesg = ::testing::internal::GetCapturedStdout(); + ASSERT_THAT(mesg, MatchesRegex(".*ERROR: Illegal reset_mol_ids command.*")); } } // namespace LAMMPS_NS From f0af7c686acb58508d6f8f70d794781493301ed5 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 9 Jul 2020 21:39:55 -0400 Subject: [PATCH 28/29] more death tests to reach 100% coverage --- unittest/commands/test_reset_ids.cpp | 38 ++++++++++++++++++++++++++++ 1 file changed, 38 insertions(+) diff --git a/unittest/commands/test_reset_ids.cpp b/unittest/commands/test_reset_ids.cpp index 81d4b01ab5..f6570dd4d5 100644 --- a/unittest/commands/test_reset_ids.cpp +++ b/unittest/commands/test_reset_ids.cpp @@ -61,6 +61,7 @@ protected: lmp->input->one("variable input_dir index " STRINGIFY(TEST_INPUT_FOLDER)); lmp->input->one("include ${input_dir}/in.fourmol"); } + delete info; if (!verbose) ::testing::internal::GetCapturedStdout(); } @@ -499,6 +500,43 @@ TEST_F(ResetIDsTest, DeathTests) mesg = ::testing::internal::GetCapturedStdout(); ASSERT_THAT(mesg, MatchesRegex(".*ERROR: Illegal reset_mol_ids command.*")); } + +TEST(ResetMolIds, CMDFail) +{ + LAMMPS *lmp; + const char *args[] = {"ResetIDsTest", "-log", "none", "-nocite", "-echo", "screen"}; + char **argv = (char **)args; + int argc = sizeof(args) / sizeof(char *); + if (!verbose) ::testing::internal::CaptureStdout(); + lmp = new LAMMPS(argc, argv, MPI_COMM_WORLD); + if (!verbose) ::testing::internal::GetCapturedStdout(); + + ::testing::internal::CaptureStdout(); + TEST_FAILURE(lmp->input->one("reset_mol_ids all");); + auto mesg = ::testing::internal::GetCapturedStdout(); + ASSERT_THAT(mesg, MatchesRegex(".*ERROR: Reset_mol_ids command before.*")); + + ::testing::internal::CaptureStdout(); + lmp->input->one("atom_modify id no"); + lmp->input->one("region box block 0 1 0 1 0 1"); + lmp->input->one("create_box 1 box"); + TEST_FAILURE(lmp->input->one("reset_mol_ids all");); + mesg = ::testing::internal::GetCapturedStdout(); + ASSERT_THAT(mesg, MatchesRegex(".*ERROR: Cannot use reset_mol_ids unl.*")); + + ::testing::internal::CaptureStdout(); + lmp->input->one("clear"); + lmp->input->one("region box block 0 1 0 1 0 1"); + lmp->input->one("create_box 1 box"); + TEST_FAILURE(lmp->input->one("reset_mol_ids all");); + mesg = ::testing::internal::GetCapturedStdout(); + ASSERT_THAT(mesg, MatchesRegex(".*ERROR: Can only use reset_mol_ids.*")); + + if (!verbose) ::testing::internal::CaptureStdout(); + delete lmp; + if (!verbose) ::testing::internal::GetCapturedStdout(); +}; + } // namespace LAMMPS_NS int main(int argc, char **argv) From 0c89b517a5b6052de2bd804315676461df05c748 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 10 Jul 2020 08:25:37 -0400 Subject: [PATCH 29/29] avoid segfaults if fewer than 10 atoms or bounding box length is zero --- src/reset_atom_ids.cpp | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/src/reset_atom_ids.cpp b/src/reset_atom_ids.cpp index 5ab7400a52..8c7d1cafc9 100644 --- a/src/reset_atom_ids.cpp +++ b/src/reset_atom_ids.cpp @@ -318,6 +318,15 @@ void ResetIDs::sort() if (dim == 2) mylo[2] = myhi[2] = 0.0; + // must ensure that bounding box volume is > 0.0 + + for (int i = 0; i < 3; ++i) { + if (mylo[i] == myhi[i]) { + mylo[i] -= 0.5; + myhi[i] += 0.5; + } + } + MPI_Allreduce(mylo,bboxlo,3,MPI_DOUBLE,MPI_MIN,world); MPI_Allreduce(myhi,bboxhi,3,MPI_DOUBLE,MPI_MAX,world); @@ -332,7 +341,8 @@ void ResetIDs::sort() // binsize = edge length of a cubic bin // nbin xyz = bin count in each dimension - bigint nbin_estimate = atom->natoms / PERBIN; + bigint nbin_estimate = atom->natoms / PERBIN + 1; + double vol; if (dim == 2) vol = (bboxhi[0]-bboxlo[0]) * (bboxhi[1]-bboxlo[1]); else vol = (bboxhi[0]-bboxlo[0]) * (bboxhi[1]-bboxlo[1]) * (bboxhi[2]-bboxlo[2]);