From 2c8a49bb2600d4ebcf54c8a991a72d2845d36d8a Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Sat, 29 Apr 2023 18:14:24 -0600 Subject: [PATCH 01/23] new compute count_type --- src/compute_count_type.cpp | 165 +++++++++++++++++++++++++++++++++++++ src/compute_count_type.h | 46 +++++++++++ src/lmptype.h | 5 +- 3 files changed, 215 insertions(+), 1 deletion(-) create mode 100644 src/compute_count_type.cpp create mode 100644 src/compute_count_type.h diff --git a/src/compute_count_type.cpp b/src/compute_count_type.cpp new file mode 100644 index 0000000000..707a4bb943 --- /dev/null +++ b/src/compute_count_type.cpp @@ -0,0 +1,165 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "compute_count_type.h" + +#include "atom.h" +#include "atom_vec.h" +#include "domain.h" +#include "error.h" +#include "force.h" +#include "group.h" +#include "update.h" + +using namespace LAMMPS_NS; + +enum{ATOM,BOND}; + +/* ---------------------------------------------------------------------- */ + +ComputeCountType::ComputeCountType(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) +{ + if (narg != 4) error->all(FLERR, "Illegal compute count/type command"); + + // process args + + if (strcmp(arg[3],"atom") == 0) mode = ATOM; + else if (strcmp(arg[3],"bond") == 0) mode = BOND; + else error->all(FLERR, "Invalid compute count/type keyword {}",arg[3]); + + if (mode == ATOM) { + vector_flag = 1; + size_vector = atom->ntypes; + extvector = 1; + } else if (mode == BOND) { + scalar_flag = vector_flag = 1; + size_vector = atom->nbondtypes; + extscalar = 1; + extvector = 1; + } + + if (mode == BOND && !atom->avec->bonds_allow) + error->all(FLERR,"Cannot use compute count/type bond command with no bonds allowed"); + + vector = new double[size_vector]; + + // work vectors + + count = new int[size_vector]; + bcount_me = new bigint[size_vector]; + bcount = new bigint[size_vector]; +} + +/* ---------------------------------------------------------------------- */ + +ComputeCountType::~ComputeCountType() +{ + delete[] vector; +} + +/* ---------------------------------------------------------------------- + only invoked for mode = BOND to count broken bonds + broken bonds have bond_type = 0 +---------------------------------------------------------------------- */ + +double ComputeCountType::compute_scalar() +{ + invoked_scalar = update->ntimestep; + + int *num_bond = atom->num_bond; + int **bond_type = atom->bond_type; + int nlocal = atom->nlocal; + + int nbond; + int count = 0; + + for (int i = 0; i < nlocal; i++) { + nbond = num_bond[i]; + for (int m = 0; m < nbond; m++) + if (bond_type[i][m] == 0) count++; + } + + // sum across procs as bigint, then convert to double + // correct for double counting if newton_bond off + + bigint bcount_me = count; + bigint bcount; + MPI_Allreduce(&bcount_me, &bcount, 1, MPI_LMP_BIGINT, MPI_SUM, world); + if (force->newton_bond == 0) bcount /= 2; + + if (bcount > MAXDOUBLEINT) + error->all(FLERR,"Compute count/type overflow"); + double scalar = bcount; + return scalar; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeCountType::compute_vector() +{ + invoked_vector = update->ntimestep; + + int n; + + // count atoms by type + + if (mode == ATOM) { + int *type = atom->type; + int nlocal = atom->nlocal; + int ntypes = atom->ntypes; + + for (int m = 0; m < ntypes; m++) count[m] = 0; + for (int i = 0; i < nlocal; i++) count[type[i]-1]++; + + n = ntypes; + } + + // count bonds by type + // skip type = 0 bonds, they are counted by compute_scalar + // bond types can be negative for SHAKE + + else if (mode == BOND) { + int *num_bond = atom->num_bond; + int **bond_type = atom->bond_type; + int nlocal = atom->nlocal; + int nbondtypes = atom->nbondtypes; + + int nbond,itype; + for (int m = 0; m < nbondtypes; m++) count[m] = 0; + + for (int i = 0; i < nlocal; i++) { + nbond = num_bond[i]; + for (int m = 0; m < nbond; m++) { + itype = bond_type[i][m]; + if (itype == 0) continue; + if (itype > 0) count[itype-1]++; + else count[-itype-1]++; + } + } + + n = nbondtypes; + } + + // sum across procs as bigint, then convert to double + // correct for double counting if newton_bond off + + for (int m = 0; m < n; m++) bcount_me[m] = count[m]; + MPI_Allreduce(&bcount_me, &bcount, 1, MPI_LMP_BIGINT, MPI_SUM, world); + if (force->newton_bond == 0) + for (int m = 0; m < n; m++) bcount[m] /= 2; + + for (int m = 0; m < n; m++) + if (bcount[m] > MAXDOUBLEINT) + error->all(FLERR,"Compute count/type overflow"); + for (int m = 0; m < n; m++) vector[m] *= bcount[m]; +} diff --git a/src/compute_count_type.h b/src/compute_count_type.h new file mode 100644 index 0000000000..0ed7a22467 --- /dev/null +++ b/src/compute_count_type.h @@ -0,0 +1,46 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef COMPUTE_CLASS +// clang-format off +ComputeStyle(count/type,ComputeCountType); +// clang-format on +#else + +#ifndef LMP_COMPUTE_COMPUTE_TYPE_H +#define LMP_COMPUTE_COMPUTE_TYPE_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeCountType : public Compute { + public: + ComputeCountType(class LAMMPS *, int, char **); + ~ComputeCountType() override; + void init() override {} + double compute_scalar() override; + void compute_vector() override; + + protected: + int mode; + + int *count; + bigint *bcount_me; + bigint *bcount; +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/lmptype.h b/src/lmptype.h index bf56d07d77..891d5bdf89 100644 --- a/src/lmptype.h +++ b/src/lmptype.h @@ -95,7 +95,8 @@ typedef int64_t bigint; #define MAXSMALLINT INT_MAX #define MAXTAGINT INT_MAX #define MAXBIGINT INT64_MAX - +#define MAXDOUBLEINT 9007199254740992 // 2^53 + #define MPI_LMP_TAGINT MPI_INT #define MPI_LMP_IMAGEINT MPI_INT #define MPI_LMP_BIGINT MPI_LL @@ -132,6 +133,7 @@ typedef int64_t bigint; #define MAXSMALLINT INT_MAX #define MAXTAGINT INT64_MAX #define MAXBIGINT INT64_MAX +#define MAXDOUBLEINT 9007199254740992 // 2^53 #define MPI_LMP_TAGINT MPI_LL #define MPI_LMP_IMAGEINT MPI_LL @@ -168,6 +170,7 @@ typedef int bigint; #define MAXSMALLINT INT_MAX #define MAXTAGINT INT_MAX #define MAXBIGINT INT_MAX +#define MAXDOUBLEINT INT_MAX #define MPI_LMP_TAGINT MPI_INT #define MPI_LMP_IMAGEINT MPI_INT From f644f7078cb782702df393dafce4248c36c84ff9 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Sat, 29 Apr 2023 18:50:15 -0600 Subject: [PATCH 02/23] doc page for new command --- doc/src/Commands_compute.rst | 1 + doc/src/bond_style.rst | 14 +++--- doc/src/compute.rst | 1 + doc/src/compute_count_type.rst | 89 ++++++++++++++++++++++++++++++++++ src/compute_count_type.cpp | 8 ++- 5 files changed, 105 insertions(+), 8 deletions(-) create mode 100644 doc/src/compute_count_type.rst diff --git a/doc/src/Commands_compute.rst b/doc/src/Commands_compute.rst index 2f9c8b03fb..755000c976 100644 --- a/doc/src/Commands_compute.rst +++ b/doc/src/Commands_compute.rst @@ -46,6 +46,7 @@ KOKKOS, o = OPENMP, t = OPT. * :doc:`com/chunk ` * :doc:`contact/atom ` * :doc:`coord/atom (k) ` + * :doc:`count/type ` * :doc:`damage/atom ` * :doc:`dihedral ` * :doc:`dihedral/local ` diff --git a/doc/src/bond_style.rst b/doc/src/bond_style.rst index b33d0a9e9a..95f463e695 100644 --- a/doc/src/bond_style.rst +++ b/doc/src/bond_style.rst @@ -32,13 +32,13 @@ Set the formula(s) LAMMPS uses to compute bond interactions between pairs of atoms. In LAMMPS, a bond differs from a pairwise interaction, which are set via the :doc:`pair_style ` command. Bonds are defined between specified pairs of atoms and -remain in force for the duration of the simulation (unless the bond -breaks which is possible in some bond potentials). The list of bonded -atoms is read in by a :doc:`read_data ` or -:doc:`read_restart ` command from a data or restart file. -By contrast, pair potentials are typically defined between all pairs -of atoms within a cutoff distance and the set of active interactions -changes over time. +remain in force for the duration of the simulation (unless new bonds +are created or existing bonds break, which is possible in some fixes +and bond potentials). The list of bonded atoms is read in by a +:doc:`read_data ` or :doc:`read_restart ` +command from a data or restart file. By contrast, pair potentials are +typically defined between all pairs of atoms within a cutoff distance +and the set of active interactions changes over time. Hybrid models where bonds are computed using different bond potentials can be setup using the *hybrid* bond style. diff --git a/doc/src/compute.rst b/doc/src/compute.rst index 880f60a8a6..950487cb72 100644 --- a/doc/src/compute.rst +++ b/doc/src/compute.rst @@ -200,6 +200,7 @@ The individual style names on the :doc:`Commands compute ` pag * :doc:`com/chunk ` - center of mass for each chunk * :doc:`contact/atom ` - contact count for each spherical particle * :doc:`coord/atom ` - coordination number for each atom +* :doc:`count/type ` - count of atoms or bonds by type * :doc:`damage/atom ` - Peridynamic damage for each atom * :doc:`dihedral ` - energy of each dihedral sub-style * :doc:`dihedral/local ` - angle of each dihedral diff --git a/doc/src/compute_count_type.rst b/doc/src/compute_count_type.rst new file mode 100644 index 0000000000..17fd9ed784 --- /dev/null +++ b/doc/src/compute_count_type.rst @@ -0,0 +1,89 @@ +.. index:: compute count/type + +compute count/type command +==================== + +Syntax +"""""" + +.. code-block:: LAMMPS + + compute ID group-ID count/type mode + +* ID, group-ID are documented in :doc:`compute ` command +* count/type = style name of this compute command +* mode = {atom} or {bond} + +Examples +"""""""" + +.. code-block:: LAMMPS + + compute 1 all count/type atom + compute 1 flowmols count/type bond + +Description +""""""""""" + +Define a computation that counts the current number of atoms by atom +type or the number of bonds by bond type. The latter can be useful in +reactive simulations where bonds are broken or created. + +Note that for this command, bonds are the topological kind enumerated +in a data file, initially read by the :doc:`read_data ` +command. They do not refer to bonds defined on-the-fly by bond-order +or reactive pair styles. + +These commands can create and break toplogical bonds: + +* :doc:`fix bond/react ` +* :doc:`fix bond/create ` +* :doc:`fix bond/break ` +* :doc:`bond_style quartic ` +* :doc:`BPM package ` bond styles + +If the {mode} setting is {atom} then the count of atoms for each atom +type is tallied. Only atoms in the specified group are counted. + +If the {mode} setting is {bond} then the count of bonds for each bond +type is tallied. Only bonds with both atoms in the specified group +are counted. + +For {mode} = {bond}, broken bonds with a bond type or zero are also +counted. Some commands flag broken bonds by setting their bond type +to zero. See the :doc:`Howto broken bonds ` doc +page for details. Note that the group setting is ignored for broken +bonds; all broken bonds in the system are counted. + +---------- + +Output info +""""""""""" + +This compute calculates a global vector of counts. If the mode is +{atom}, the vector length is the number of atom types. If the mode is +{bond}, the vector length is the number of bond types. + +If the mode is {bond} this compute also calculates a global scalar which +counts the number of broken bonds. + +These values can be used by any command that uses global scalar or +vector values from a compute as input. See the :doc:`Howto output +` page for an overview of LAMMPS output options. + +The scalar and vector values calculated by this compute are "extensive". + +Restrictions +"""""""""""" + +none + +Related commands +"""""""""""""""" + +none + +Default +""""""" + +none diff --git a/src/compute_count_type.cpp b/src/compute_count_type.cpp index 707a4bb943..c3014394d8 100644 --- a/src/compute_count_type.cpp +++ b/src/compute_count_type.cpp @@ -82,6 +82,8 @@ double ComputeCountType::compute_scalar() int nbond; int count = 0; + + // NOTE: respect group setting for (int i = 0; i < nlocal; i++) { nbond = num_bond[i]; @@ -113,6 +115,8 @@ void ComputeCountType::compute_vector() // count atoms by type + // NOTE: respect group setting + if (mode == ATOM) { int *type = atom->type; int nlocal = atom->nlocal; @@ -127,7 +131,9 @@ void ComputeCountType::compute_vector() // count bonds by type // skip type = 0 bonds, they are counted by compute_scalar // bond types can be negative for SHAKE - + + // NOTE: respect group setting + else if (mode == BOND) { int *num_bond = atom->num_bond; int **bond_type = atom->bond_type; From c282d8d5d0eb80c9050f8b141639ef1e31617e06 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Mon, 1 May 2023 09:01:03 -0600 Subject: [PATCH 03/23] add group support to new command --- doc/src/compute_count_type.rst | 11 ++++----- src/compute_count_type.cpp | 45 ++++++++++++++++++++++++---------- 2 files changed, 37 insertions(+), 19 deletions(-) diff --git a/doc/src/compute_count_type.rst b/doc/src/compute_count_type.rst index 17fd9ed784..ca56648297 100644 --- a/doc/src/compute_count_type.rst +++ b/doc/src/compute_count_type.rst @@ -29,7 +29,7 @@ Define a computation that counts the current number of atoms by atom type or the number of bonds by bond type. The latter can be useful in reactive simulations where bonds are broken or created. -Note that for this command, bonds are the topological kind enumerated +Note that for this command, bonds are the topological ones enumerated in a data file, initially read by the :doc:`read_data ` command. They do not refer to bonds defined on-the-fly by bond-order or reactive pair styles. @@ -49,9 +49,8 @@ If the {mode} setting is {bond} then the count of bonds for each bond type is tallied. Only bonds with both atoms in the specified group are counted. -For {mode} = {bond}, broken bonds with a bond type or zero are also -counted. Some commands flag broken bonds by setting their bond type -to zero. See the :doc:`Howto broken bonds ` doc +For {mode} = {bond}, broken bonds with a bond type of zero are also +counted. See the :doc:`Howto broken bonds ` doc page for details. Note that the group setting is ignored for broken bonds; all broken bonds in the system are counted. @@ -64,8 +63,8 @@ This compute calculates a global vector of counts. If the mode is {atom}, the vector length is the number of atom types. If the mode is {bond}, the vector length is the number of bond types. -If the mode is {bond} this compute also calculates a global scalar which -counts the number of broken bonds. +If the mode is {bond} this compute also calculates a global scalar +which is the number of broken bonds. These values can be used by any command that uses global scalar or vector values from a compute as input. See the :doc:`Howto output diff --git a/src/compute_count_type.cpp b/src/compute_count_type.cpp index c3014394d8..9de5f39cab 100644 --- a/src/compute_count_type.cpp +++ b/src/compute_count_type.cpp @@ -80,14 +80,16 @@ double ComputeCountType::compute_scalar() int **bond_type = atom->bond_type; int nlocal = atom->nlocal; - int nbond; + int m,nbond; int count = 0; - // NOTE: respect group setting + // count broken bonds with bond_type = 0 + // ignore group setting since 2 atoms in a broken bond + // can be arbitrarily far apart for (int i = 0; i < nlocal; i++) { nbond = num_bond[i]; - for (int m = 0; m < nbond; m++) + for (m = 0; m < nbond; m++) if (bond_type[i][m] == 0) count++; } @@ -114,45 +116,62 @@ void ComputeCountType::compute_vector() int n; // count atoms by type - - // NOTE: respect group setting + // atom must be in group to be counted if (mode == ATOM) { int *type = atom->type; + int *mask = atom->mask; int nlocal = atom->nlocal; int ntypes = atom->ntypes; for (int m = 0; m < ntypes; m++) count[m] = 0; - for (int i = 0; i < nlocal; i++) count[type[i]-1]++; + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + count[type[i]-1]++; n = ntypes; } // count bonds by type + // both atoms in bond must be in group to be counted // skip type = 0 bonds, they are counted by compute_scalar // bond types can be negative for SHAKE - // NOTE: respect group setting - else if (mode == BOND) { - int *num_bond = atom->num_bond; + int **bond_atom = atom->bond_atom; int **bond_type = atom->bond_type; + int *num_bond = atom->num_bond; + int *mask = atom->mask; int nlocal = atom->nlocal; int nbondtypes = atom->nbondtypes; - int nbond,itype; + int j,m,nbond,itype; + int flag = 0; for (int m = 0; m < nbondtypes; m++) count[m] = 0; for (int i = 0; i < nlocal; i++) { nbond = num_bond[i]; - for (int m = 0; m < nbond; m++) { + for (m = 0; m < nbond; m++) { itype = bond_type[i][m]; if (itype == 0) continue; - if (itype > 0) count[itype-1]++; - else count[-itype-1]++; + + j = atom->map(bond_atom[i][m]); + if (j < 0) { + flag = 1; + continue; + } + + if ((mask[i] & groupbit) && (mask[j] & groupbit)) { + if (itype > 0) count[itype-1]++; + else count[-itype-1]++; + } } } + int flagany; + MPI_Allreduce(&flag, &flagany, 1, MPI_INT, MPI_SUM, world); + if (flagany) error->all(FLERR,"Missing bond atom in compute count/type"); + n = nbondtypes; } From 2fe423cc70a104514b8f2c396e3b732a86c74d0d Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Mon, 1 May 2023 09:20:20 -0600 Subject: [PATCH 04/23] bug fix --- src/compute_count_type.cpp | 24 +++++++++++++++--------- 1 file changed, 15 insertions(+), 9 deletions(-) diff --git a/src/compute_count_type.cpp b/src/compute_count_type.cpp index 9de5f39cab..feadcb04ba 100644 --- a/src/compute_count_type.cpp +++ b/src/compute_count_type.cpp @@ -51,6 +51,8 @@ ComputeCountType::ComputeCountType(LAMMPS *lmp, int narg, char **arg) : Compute( if (mode == BOND && !atom->avec->bonds_allow) error->all(FLERR,"Cannot use compute count/type bond command with no bonds allowed"); + // output vector + vector = new double[size_vector]; // work vectors @@ -65,6 +67,10 @@ ComputeCountType::ComputeCountType(LAMMPS *lmp, int narg, char **arg) : Compute( ComputeCountType::~ComputeCountType() { delete[] vector; + + delete[] count; + delete[] bcount_me; + delete[] bcount; } /* ---------------------------------------------------------------------- @@ -103,7 +109,7 @@ double ComputeCountType::compute_scalar() if (bcount > MAXDOUBLEINT) error->all(FLERR,"Compute count/type overflow"); - double scalar = bcount; + scalar = bcount; return scalar; } @@ -113,7 +119,7 @@ void ComputeCountType::compute_vector() { invoked_vector = update->ntimestep; - int n; + int nvec; // count atoms by type // atom must be in group to be counted @@ -129,7 +135,7 @@ void ComputeCountType::compute_vector() if (mask[i] & groupbit) count[type[i]-1]++; - n = ntypes; + nvec = ntypes; } // count bonds by type @@ -172,19 +178,19 @@ void ComputeCountType::compute_vector() MPI_Allreduce(&flag, &flagany, 1, MPI_INT, MPI_SUM, world); if (flagany) error->all(FLERR,"Missing bond atom in compute count/type"); - n = nbondtypes; + nvec = nbondtypes; } // sum across procs as bigint, then convert to double // correct for double counting if newton_bond off - for (int m = 0; m < n; m++) bcount_me[m] = count[m]; - MPI_Allreduce(&bcount_me, &bcount, 1, MPI_LMP_BIGINT, MPI_SUM, world); + for (int m = 0; m < nvec; m++) bcount_me[m] = count[m]; + MPI_Allreduce(bcount_me, bcount, nvec, MPI_LMP_BIGINT, MPI_SUM, world); if (force->newton_bond == 0) - for (int m = 0; m < n; m++) bcount[m] /= 2; + for (int m = 0; m < nvec; m++) bcount[m] /= 2; - for (int m = 0; m < n; m++) + for (int m = 0; m < nvec; m++) if (bcount[m] > MAXDOUBLEINT) error->all(FLERR,"Compute count/type overflow"); - for (int m = 0; m < n; m++) vector[m] *= bcount[m]; + for (int m = 0; m < nvec; m++) vector[m] = bcount[m]; } From ef9ce62aa8f89b9ce067cda9a1b8fc6f6bb78cd4 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Mon, 1 May 2023 15:17:19 -0600 Subject: [PATCH 05/23] update broken bond doc page --- doc/src/Howto_broken_bonds.rst | 75 ++++++++++++++++++---------------- 1 file changed, 40 insertions(+), 35 deletions(-) diff --git a/doc/src/Howto_broken_bonds.rst b/doc/src/Howto_broken_bonds.rst index 88c8acb87f..4fd0bcc30f 100644 --- a/doc/src/Howto_broken_bonds.rst +++ b/doc/src/Howto_broken_bonds.rst @@ -1,48 +1,53 @@ Broken Bonds ============ -Typically, bond interactions persist for the duration of a simulation in -LAMMPS. However, there are some exceptions that allow for bonds to -break, including the :doc:`quartic bond style ` and the -bond styles in the :doc:`BPM package ` which contains the -:doc:`bpm/spring ` and :doc:`bpm/rotational -` bond styles. In these cases, a bond can be broken -if it is stretched beyond a user-defined threshold. LAMMPS accomplishes -this by setting the bond type to 0, such that the bond force is no -longer computed. +Typically, bond interactions persist for the duration of a simulation +in LAMMPS. However, some commands allow bonds to break, including the +following: -Users are normally able to weight the contribution of pair forces to atoms -that are bonded using the :doc:`special_bonds command `. -When bonds break, this is not always the case. For the quartic bond style, -pair forces are always turned off between bonded particles. LAMMPS does -this via a computational sleight-of-hand. It subtracts the pairwise -interaction as part of the bond computation. When the bond breaks, the -subtraction stops. For this to work, the pairwise interaction must always -be computed by the :doc:`pair_style ` command, whether the bond -is broken or not. This means that :doc:`special_bonds ` must -be set to 1,1,1. After the bond breaks, the pairwise interaction between the -two atoms is turned on, since they are no longer bonded. +* :doc:`bond_style quartic ` +* :doc:`fix bond/break ` +* :doc:`fix bond/react ` +* :doc:`BPM package ` bond styles -In the BPM package, one can either turn off all pair interactions between -bonded particles or leave them on, overlaying pair forces on top of bond -forces. To remove pair forces, the special bond list is dynamically -updated. More details can be found on the :doc:`Howto BPM ` -page. +A bond can break if it is stretched beyond a user-defined threshold or +more generally if other criteria are met. -Bonds can also be broken by fixes which change bond topology, including -:doc:`fix bond/break ` and -:doc:`fix bond/react `. These fixes will automatically -trigger a rebuild of the neighbor list and update special bond data structures -when bonds are broken. +For the quartic bond style, when a bond is broken its bond type is set +to 0 and pairwise forces between the two atoms in the broken bond are +"turned on". Angles, dihedrals, etc cannot be defined for the system +when :doc:`bond_style quartic ` is used. -Note that when bonds are dumped to a file via the :doc:`dump local ` command, bonds with type 0 are not included. The -:doc:`delete_bonds ` command can also be used to query the -status of broken bonds or permanently delete them, e.g.: +The :doc:`fix bond/break ` and :doc:`fix bond/react +` commands allow breaking of bonds within a molecular +topology with also defines angles, dihedrals, etc. These fixes will +update internal topology data structures when bonds are broken, so +that the appropriate angle, dihederal, etc interactions are also +turned off. They will also trigger a rebuild of the neighbor list +when this occurs, to turn on the appropriate pairwise forces. + +In the BPM package, one can either turn off all pair interactions +between bonded particles or leave them on, overlaying pair forces on +top of bond forces. To remove pair forces, the special bond list is +dynamically updated. More details can be found on the :doc:`Howto BPM +` page. + +Note that when bonds are dumped to a file via the :doc:`dump local +` command, bonds with type 0 are not included. + +The :doc:`delete_bonds ` command can also be used to +query the status of broken bonds (type = 0) or permanently delete +them, e.g.: .. code-block:: LAMMPS delete_bonds all stats delete_bonds all bond 0 remove -The compute :doc:`nbond/atom ` can also be used -to tally the current number of bonds per atom, excluding broken bonds. +The compute :doc:`count/type bond ` command +tallies the current number of bonds for each bond type. It also +tallies broken bonds with type = 0. + +The compute :doc:`nbond/atom ` command tallies the +current number of bonds each atom is part of, excluding broken bonds +with type = 0. From fae0ef6cdea8a531bb9a5c77bb8e32c8ab8284da Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Mon, 1 May 2023 15:21:56 -0600 Subject: [PATCH 06/23] update count/type doc page --- doc/src/compute_count_type.rst | 30 +++++++++++++++++++----------- 1 file changed, 19 insertions(+), 11 deletions(-) diff --git a/doc/src/compute_count_type.rst b/doc/src/compute_count_type.rst index ca56648297..9ab0b1009f 100644 --- a/doc/src/compute_count_type.rst +++ b/doc/src/compute_count_type.rst @@ -25,16 +25,22 @@ Examples Description """"""""""" -Define a computation that counts the current number of atoms by atom -type or the number of bonds by bond type. The latter can be useful in -reactive simulations where bonds are broken or created. +Define a computation that counts the current number of atoms for each +atom type or the number of bonds for each bond type. The former can +be useful if atoms are added to or deleted from the system in random +ways, e.g. via the :doc:`fix deposit `, :doc:`fix pour +`, or :doc:`fix evaporate ` commands. The +latter can be useful in reactive simulations where molecular bonds are +broken or created. -Note that for this command, bonds are the topological ones enumerated -in a data file, initially read by the :doc:`read_data ` -command. They do not refer to bonds defined on-the-fly by bond-order -or reactive pair styles. +Note that for this command, bonds are defined as the topological kind +enumerated in a data file, initially read by the :doc:`read_data +` command or defined by the :doc:`molecule ` +command. They do not refer to implicit bonds defined on-the-fly by +bond-order or reactive pair styles based on the current conformation +of small clusters of atoms. -These commands can create and break toplogical bonds: +These commands can create and/or break topological bonds: * :doc:`fix bond/react ` * :doc:`fix bond/create ` @@ -50,9 +56,11 @@ type is tallied. Only bonds with both atoms in the specified group are counted. For {mode} = {bond}, broken bonds with a bond type of zero are also -counted. See the :doc:`Howto broken bonds ` doc -page for details. Note that the group setting is ignored for broken -bonds; all broken bonds in the system are counted. +counted. Some of the commands listed above which break bonds, do this +by setting their types to zero. See the :doc:`Howto broken bonds +` doc page for details. Note that the group +setting is ignored for broken bonds; all broken bonds in the system +are counted. ---------- From 73b8bb8617cdddc2a2a812090d40bb26e2a9b539 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 1 May 2023 20:41:44 -0400 Subject: [PATCH 07/23] fix bug that breaks compilation with -DLAMMPS_BIGBIG --- src/compute_count_type.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/compute_count_type.cpp b/src/compute_count_type.cpp index feadcb04ba..4fd814cdff 100644 --- a/src/compute_count_type.cpp +++ b/src/compute_count_type.cpp @@ -144,7 +144,7 @@ void ComputeCountType::compute_vector() // bond types can be negative for SHAKE else if (mode == BOND) { - int **bond_atom = atom->bond_atom; + tagint **bond_atom = atom->bond_atom; int **bond_type = atom->bond_type; int *num_bond = atom->num_bond; int *mask = atom->mask; From b38544a9e86612e130690dcd56460c3d50802a0a Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 1 May 2023 20:43:05 -0400 Subject: [PATCH 08/23] apply clang-format --- src/compute_count_type.cpp | 50 ++++++++++++++++++++------------------ src/compute_count_type.h | 2 +- 2 files changed, 27 insertions(+), 25 deletions(-) diff --git a/src/compute_count_type.cpp b/src/compute_count_type.cpp index 4fd814cdff..786506280e 100644 --- a/src/compute_count_type.cpp +++ b/src/compute_count_type.cpp @@ -23,7 +23,7 @@ using namespace LAMMPS_NS; -enum{ATOM,BOND}; +enum { ATOM, BOND }; /* ---------------------------------------------------------------------- */ @@ -33,9 +33,12 @@ ComputeCountType::ComputeCountType(LAMMPS *lmp, int narg, char **arg) : Compute( // process args - if (strcmp(arg[3],"atom") == 0) mode = ATOM; - else if (strcmp(arg[3],"bond") == 0) mode = BOND; - else error->all(FLERR, "Invalid compute count/type keyword {}",arg[3]); + if (strcmp(arg[3], "atom") == 0) + mode = ATOM; + else if (strcmp(arg[3], "bond") == 0) + mode = BOND; + else + error->all(FLERR, "Invalid compute count/type keyword {}", arg[3]); if (mode == ATOM) { vector_flag = 1; @@ -49,10 +52,10 @@ ComputeCountType::ComputeCountType(LAMMPS *lmp, int narg, char **arg) : Compute( } if (mode == BOND && !atom->avec->bonds_allow) - error->all(FLERR,"Cannot use compute count/type bond command with no bonds allowed"); + error->all(FLERR, "Cannot use compute count/type bond command with no bonds allowed"); // output vector - + vector = new double[size_vector]; // work vectors @@ -86,13 +89,13 @@ double ComputeCountType::compute_scalar() int **bond_type = atom->bond_type; int nlocal = atom->nlocal; - int m,nbond; + int m, nbond; int count = 0; // count broken bonds with bond_type = 0 // ignore group setting since 2 atoms in a broken bond // can be arbitrarily far apart - + for (int i = 0; i < nlocal; i++) { nbond = num_bond[i]; for (m = 0; m < nbond; m++) @@ -101,14 +104,13 @@ double ComputeCountType::compute_scalar() // sum across procs as bigint, then convert to double // correct for double counting if newton_bond off - + bigint bcount_me = count; bigint bcount; MPI_Allreduce(&bcount_me, &bcount, 1, MPI_LMP_BIGINT, MPI_SUM, world); if (force->newton_bond == 0) bcount /= 2; - if (bcount > MAXDOUBLEINT) - error->all(FLERR,"Compute count/type overflow"); + if (bcount > MAXDOUBLEINT) error->all(FLERR, "Compute count/type overflow"); scalar = bcount; return scalar; } @@ -132,12 +134,11 @@ void ComputeCountType::compute_vector() for (int m = 0; m < ntypes; m++) count[m] = 0; for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) - count[type[i]-1]++; + if (mask[i] & groupbit) count[type[i] - 1]++; nvec = ntypes; } - + // count bonds by type // both atoms in bond must be in group to be counted // skip type = 0 bonds, they are counted by compute_scalar @@ -150,11 +151,11 @@ void ComputeCountType::compute_vector() int *mask = atom->mask; int nlocal = atom->nlocal; int nbondtypes = atom->nbondtypes; - - int j,m,nbond,itype; + + int j, m, nbond, itype; int flag = 0; for (int m = 0; m < nbondtypes; m++) count[m] = 0; - + for (int i = 0; i < nlocal; i++) { nbond = num_bond[i]; for (m = 0; m < nbond; m++) { @@ -166,17 +167,19 @@ void ComputeCountType::compute_vector() flag = 1; continue; } - + if ((mask[i] & groupbit) && (mask[j] & groupbit)) { - if (itype > 0) count[itype-1]++; - else count[-itype-1]++; + if (itype > 0) + count[itype - 1]++; + else + count[-itype - 1]++; } } } int flagany; MPI_Allreduce(&flag, &flagany, 1, MPI_INT, MPI_SUM, world); - if (flagany) error->all(FLERR,"Missing bond atom in compute count/type"); + if (flagany) error->all(FLERR, "Missing bond atom in compute count/type"); nvec = nbondtypes; } @@ -186,11 +189,10 @@ void ComputeCountType::compute_vector() for (int m = 0; m < nvec; m++) bcount_me[m] = count[m]; MPI_Allreduce(bcount_me, bcount, nvec, MPI_LMP_BIGINT, MPI_SUM, world); - if (force->newton_bond == 0) + if (force->newton_bond == 0) for (int m = 0; m < nvec; m++) bcount[m] /= 2; for (int m = 0; m < nvec; m++) - if (bcount[m] > MAXDOUBLEINT) - error->all(FLERR,"Compute count/type overflow"); + if (bcount[m] > MAXDOUBLEINT) error->all(FLERR, "Compute count/type overflow"); for (int m = 0; m < nvec; m++) vector[m] = bcount[m]; } diff --git a/src/compute_count_type.h b/src/compute_count_type.h index 0ed7a22467..27b1669b1d 100644 --- a/src/compute_count_type.h +++ b/src/compute_count_type.h @@ -34,7 +34,7 @@ class ComputeCountType : public Compute { protected: int mode; - + int *count; bigint *bcount_me; bigint *bcount; From ff29ef7d3160c9597b9fb637619982afcf8e0adc Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 1 May 2023 20:57:32 -0400 Subject: [PATCH 09/23] use more obvious representation of 2^53 --- src/lmptype.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/lmptype.h b/src/lmptype.h index 891d5bdf89..769f730d6a 100644 --- a/src/lmptype.h +++ b/src/lmptype.h @@ -95,8 +95,8 @@ typedef int64_t bigint; #define MAXSMALLINT INT_MAX #define MAXTAGINT INT_MAX #define MAXBIGINT INT64_MAX -#define MAXDOUBLEINT 9007199254740992 // 2^53 - +#define MAXDOUBLEINT (1L<<53) + #define MPI_LMP_TAGINT MPI_INT #define MPI_LMP_IMAGEINT MPI_INT #define MPI_LMP_BIGINT MPI_LL @@ -133,7 +133,7 @@ typedef int64_t bigint; #define MAXSMALLINT INT_MAX #define MAXTAGINT INT64_MAX #define MAXBIGINT INT64_MAX -#define MAXDOUBLEINT 9007199254740992 // 2^53 +#define MAXDOUBLEINT (1L<<53) #define MPI_LMP_TAGINT MPI_LL #define MPI_LMP_IMAGEINT MPI_LL From c8d5b9e4d072c46c1795e696e26fc09b48016510 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 1 May 2023 23:52:04 -0400 Subject: [PATCH 10/23] avoid uninitialized pointers --- src/compute_count_type.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/compute_count_type.cpp b/src/compute_count_type.cpp index 786506280e..0cd109a1eb 100644 --- a/src/compute_count_type.cpp +++ b/src/compute_count_type.cpp @@ -27,9 +27,10 @@ enum { ATOM, BOND }; /* ---------------------------------------------------------------------- */ -ComputeCountType::ComputeCountType(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) +ComputeCountType::ComputeCountType(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg), count(nullptr), bcount(nullptr), bcount_me(nullptr) { - if (narg != 4) error->all(FLERR, "Illegal compute count/type command"); + if (narg != 4) error->all(FLERR, "Incorrect number of args for compute count/type command"); // process args @@ -51,7 +52,7 @@ ComputeCountType::ComputeCountType(LAMMPS *lmp, int narg, char **arg) : Compute( extvector = 1; } - if (mode == BOND && !atom->avec->bonds_allow) + if ((mode == BOND) && !atom->avec->bonds_allow) error->all(FLERR, "Cannot use compute count/type bond command with no bonds allowed"); // output vector From 12135bac770ecc591b05a31ad437a72038553b14 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 2 May 2023 10:22:37 -0400 Subject: [PATCH 11/23] add unit tests for compute count/type --- unittest/commands/test_compute_global.cpp | 54 +++++++++++++++++++++++ 1 file changed, 54 insertions(+) diff --git a/unittest/commands/test_compute_global.cpp b/unittest/commands/test_compute_global.cpp index d7573d8d16..4efc3f1b91 100644 --- a/unittest/commands/test_compute_global.cpp +++ b/unittest/commands/test_compute_global.cpp @@ -294,6 +294,60 @@ TEST_F(ComputeGlobalTest, Reduction) EXPECT_DOUBLE_EQ(rep[2], 26); EXPECT_DOUBLE_EQ(rep[3], max[0]); } + +TEST_F(ComputeGlobalTest, Counts) +{ + if (lammps_get_natoms(lmp) == 0.0) GTEST_SKIP(); + + BEGIN_HIDE_OUTPUT(); + command("pair_style lj/cut 10.0"); + command("pair_coeff * * 0.01 3.0"); + command("bond_style harmonic"); + command("bond_coeff * 100.0 1.5"); + + command("variable t1 atom type==1"); + command("variable t2 atom type==2"); + command("variable t3 atom type==3"); + command("variable t4 atom type==4"); + command("variable t5 atom type==5"); + command("compute asum all reduce sum v_t1 v_t2 v_t3 v_t4 v_t5"); + command("compute acnt all count/type atom"); + command("compute bcnt all count/type bond"); + command("thermo_style custom c_asum[*] c_acnt[*]"); + command("run 0 post no"); + END_HIDE_OUTPUT(); + + auto asum = get_vector("asum"); + auto acnt = get_vector("acnt"); + auto bcnt = get_vector("bcnt"); + auto bbrk = get_scalar("bcnt"); + + EXPECT_DOUBLE_EQ(bbrk, 0.0); + + EXPECT_DOUBLE_EQ(asum[0], acnt[0]); + EXPECT_DOUBLE_EQ(asum[1], acnt[1]); + EXPECT_DOUBLE_EQ(asum[2], acnt[2]); + EXPECT_DOUBLE_EQ(asum[3], acnt[3]); + EXPECT_DOUBLE_EQ(asum[4], acnt[4]); + + EXPECT_DOUBLE_EQ(bcnt[0], 3.0); + EXPECT_DOUBLE_EQ(bcnt[1], 6.0); + EXPECT_DOUBLE_EQ(bcnt[2], 3.0); + EXPECT_DOUBLE_EQ(bcnt[3], 2.0); + EXPECT_DOUBLE_EQ(bcnt[4], 10.0); + + BEGIN_HIDE_OUTPUT(); + command("delete_bonds all bond 3"); + command("run 0 post no"); + END_HIDE_OUTPUT(); + + EXPECT_DOUBLE_EQ(bbrk, 0.0); // should be 3 + EXPECT_DOUBLE_EQ(bcnt[0], 3.0); + EXPECT_DOUBLE_EQ(bcnt[1], 6.0); + EXPECT_DOUBLE_EQ(bcnt[2], 3.0); // should be 0 + EXPECT_DOUBLE_EQ(bcnt[3], 2.0); + EXPECT_DOUBLE_EQ(bcnt[4], 10.0); +} } // namespace LAMMPS_NS int main(int argc, char **argv) From eacb420e211452392442b34b26add5693dd3c092 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 2 May 2023 10:34:36 -0400 Subject: [PATCH 12/23] whitespace --- doc/src/compute_count_type.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/compute_count_type.rst b/doc/src/compute_count_type.rst index 9ab0b1009f..8cbd6aec9b 100644 --- a/doc/src/compute_count_type.rst +++ b/doc/src/compute_count_type.rst @@ -13,7 +13,7 @@ Syntax * ID, group-ID are documented in :doc:`compute ` command * count/type = style name of this compute command * mode = {atom} or {bond} - + Examples """""""" From d2bb1b420da791f73b141fc1ad52027b17846d16 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 2 May 2023 11:49:13 -0400 Subject: [PATCH 13/23] simplify and cleanup --- unittest/commands/test_compute_global.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/unittest/commands/test_compute_global.cpp b/unittest/commands/test_compute_global.cpp index 4efc3f1b91..7e5b3d4483 100644 --- a/unittest/commands/test_compute_global.cpp +++ b/unittest/commands/test_compute_global.cpp @@ -300,10 +300,8 @@ TEST_F(ComputeGlobalTest, Counts) if (lammps_get_natoms(lmp) == 0.0) GTEST_SKIP(); BEGIN_HIDE_OUTPUT(); - command("pair_style lj/cut 10.0"); - command("pair_coeff * * 0.01 3.0"); - command("bond_style harmonic"); - command("bond_coeff * 100.0 1.5"); + command("pair_style zero 10.0"); + command("pair_coeff * *"); command("variable t1 atom type==1"); command("variable t2 atom type==2"); @@ -341,6 +339,8 @@ TEST_F(ComputeGlobalTest, Counts) command("run 0 post no"); END_HIDE_OUTPUT(); + bcnt = get_vector("bcnt"); + bbrk = get_scalar("bcnt"); EXPECT_DOUBLE_EQ(bbrk, 0.0); // should be 3 EXPECT_DOUBLE_EQ(bcnt[0], 3.0); EXPECT_DOUBLE_EQ(bcnt[1], 6.0); From aa4447413ae32c3c51db83cb75facb8eafcb0d12 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Tue, 2 May 2023 13:39:37 -0600 Subject: [PATCH 14/23] expand to angles, dihedrals, impropers --- doc/src/Howto_broken_bonds.rst | 35 ++-- doc/src/compute_count_type.rst | 68 +++++-- src/compute_count_type.cpp | 330 ++++++++++++++++++++++++++------- src/compute_count_type.h | 6 + 4 files changed, 335 insertions(+), 104 deletions(-) diff --git a/doc/src/Howto_broken_bonds.rst b/doc/src/Howto_broken_bonds.rst index 4fd0bcc30f..282d2b32a4 100644 --- a/doc/src/Howto_broken_bonds.rst +++ b/doc/src/Howto_broken_bonds.rst @@ -1,9 +1,9 @@ Broken Bonds ============ -Typically, bond interactions persist for the duration of a simulation -in LAMMPS. However, some commands allow bonds to break, including the -following: +Typically, molecular bond interactions persist for the duration of a +simulation in LAMMPS. However, some commands break bonds dynamically, +including the following: * :doc:`bond_style quartic ` * :doc:`fix bond/break ` @@ -14,39 +14,40 @@ A bond can break if it is stretched beyond a user-defined threshold or more generally if other criteria are met. For the quartic bond style, when a bond is broken its bond type is set -to 0 and pairwise forces between the two atoms in the broken bond are -"turned on". Angles, dihedrals, etc cannot be defined for the system -when :doc:`bond_style quartic ` is used. +to 0 to effectively break it and pairwise forces between the two atoms +in the broken bond are "turned on". Angles, dihedrals, etc cannot be +defined for a system when :doc:`bond_style quartic ` is +used. The :doc:`fix bond/break ` and :doc:`fix bond/react ` commands allow breaking of bonds within a molecular -topology with also defines angles, dihedrals, etc. These fixes will -update internal topology data structures when bonds are broken, so -that the appropriate angle, dihederal, etc interactions are also -turned off. They will also trigger a rebuild of the neighbor list +topology with may also define angles, dihedrals, etc. These commands +update internal topology data structures to remove broken bonds, as +well as the appropriate angle, dihederal, etc interactions which +include the bond. They also trigger a rebuild of the neighbor list when this occurs, to turn on the appropriate pairwise forces. In the BPM package, one can either turn off all pair interactions between bonded particles or leave them on, overlaying pair forces on top of bond forces. To remove pair forces, the special bond list is -dynamically updated. More details can be found on the :doc:`Howto BPM +dynamically updated. More details can be found on the :doc:`Howto BPM ` page. Note that when bonds are dumped to a file via the :doc:`dump local ` command, bonds with type 0 are not included. -The :doc:`delete_bonds ` command can also be used to -query the status of broken bonds (type = 0) or permanently delete -them, e.g.: +The :doc:`delete_bonds ` command can be used to query +the status of broken bonds with type = 0 or permanently delete them, +e.g.: .. code-block:: LAMMPS delete_bonds all stats delete_bonds all bond 0 remove -The compute :doc:`count/type bond ` command -tallies the current number of bonds for each bond type. It also -tallies broken bonds with type = 0. +The compute :doc:`count/type ` command tallies the +current number of bonds (or angles, etc) for each bond (angle, etc) +type. It also tallies broken bonds with type = 0. The compute :doc:`nbond/atom ` command tallies the current number of bonds each atom is part of, excluding broken bonds diff --git a/doc/src/compute_count_type.rst b/doc/src/compute_count_type.rst index 9ab0b1009f..e9f531ef64 100644 --- a/doc/src/compute_count_type.rst +++ b/doc/src/compute_count_type.rst @@ -12,7 +12,7 @@ Syntax * ID, group-ID are documented in :doc:`compute ` command * count/type = style name of this compute command -* mode = {atom} or {bond} +* mode = {atom} or {bond} or {angle} or {dihedral} or {improper} Examples """""""" @@ -26,28 +26,46 @@ Description """"""""""" Define a computation that counts the current number of atoms for each -atom type or the number of bonds for each bond type. The former can -be useful if atoms are added to or deleted from the system in random -ways, e.g. via the :doc:`fix deposit `, :doc:`fix pour -`, or :doc:`fix evaporate ` commands. The -latter can be useful in reactive simulations where molecular bonds are -broken or created. +atom type. Or the number of bonds (angles, dihedrals, impropers) for +each bond (angle, dihedral, improper) type. -Note that for this command, bonds are defined as the topological kind -enumerated in a data file, initially read by the :doc:`read_data +The former can be useful if atoms are added to or deleted from the +system in random ways, e.g. via the :doc:`fix deposit `, +:doc:`fix pour `, or :doc:`fix evaporate ` +commands. The latter can be useful in reactive simulations where +molecular bonds are broken or created, as well as angles, dihedrals, +impropers. + +Note that for this command, bonds (angles, etc) are the topological +kind enumerated in a data file, initially read by the :doc:`read_data ` command or defined by the :doc:`molecule ` command. They do not refer to implicit bonds defined on-the-fly by bond-order or reactive pair styles based on the current conformation of small clusters of atoms. -These commands can create and/or break topological bonds: +These commands can turn off topological bonds (angles, etc) by setting +their bond (angle, etc) types to negative values. This command +includes the turned-off bonds (angles, etc) in the count for each +type: +* :doc:`fix shake ` +* :doc:`delete_bonds ` + +These commands can create and/or break topological bonds (angles, +etc). In the case of breaking, they remove the bond (angle, etc) from +the system, so that they no longer exist (:doc:`bond_style quartic +` is the exception, see the discussion below). Thus +they are not included in the counts for each type: + +* :doc:`delete_bonds remove ` +* :doc:`bond_style quartic ` * :doc:`fix bond/react ` * :doc:`fix bond/create ` * :doc:`fix bond/break ` -* :doc:`bond_style quartic ` * :doc:`BPM package ` bond styles +---------- + If the {mode} setting is {atom} then the count of atoms for each atom type is tallied. Only atoms in the specified group are counted. @@ -56,11 +74,22 @@ type is tallied. Only bonds with both atoms in the specified group are counted. For {mode} = {bond}, broken bonds with a bond type of zero are also -counted. Some of the commands listed above which break bonds, do this -by setting their types to zero. See the :doc:`Howto broken bonds -` doc page for details. Note that the group -setting is ignored for broken bonds; all broken bonds in the system -are counted. +counted. The :doc:`bond_style quartic ` breaks a bond +by doing this. See the :doc:`Howto broken bonds ` +doc page for more details. Note that the group setting is ignored for +broken bonds; all broken bonds in the system are counted. + +If the {mode} setting is {angle} then the count of angles for each +angle type is tallied. Only angles with all 3 atoms in the specified +group are counted. + +If the {mode} setting is {dihedral} then the count of dihedrals for +each dihedral type is tallied. Only dihedrals with all 4 atoms in the +specified group are counted. + +If the {mode} setting is {improper} then the count of impropers for +each improper type is tallied. Only impropers with all 4 atoms in the +specified group are counted. ---------- @@ -68,11 +97,12 @@ Output info """"""""""" This compute calculates a global vector of counts. If the mode is -{atom}, the vector length is the number of atom types. If the mode is -{bond}, the vector length is the number of bond types. +{atom} or {bond} or {angle} or {dihedral} or {improper}, then the +vector length is the number of atom types or bond types or angle types +or dihedral types or improper types, respectively. If the mode is {bond} this compute also calculates a global scalar -which is the number of broken bonds. +which is the number of broken bonds with type = 0, as explained above. These values can be used by any command that uses global scalar or vector values from a compute as input. See the :doc:`Howto output diff --git a/src/compute_count_type.cpp b/src/compute_count_type.cpp index feadcb04ba..36aa21da70 100644 --- a/src/compute_count_type.cpp +++ b/src/compute_count_type.cpp @@ -14,7 +14,6 @@ #include "compute_count_type.h" #include "atom.h" -#include "atom_vec.h" #include "domain.h" #include "error.h" #include "force.h" @@ -23,7 +22,7 @@ using namespace LAMMPS_NS; -enum{ATOM,BOND}; +enum{ATOM,BOND,ANGLE,DIHEDRAL,IMPROPER}; /* ---------------------------------------------------------------------- */ @@ -35,8 +34,24 @@ ComputeCountType::ComputeCountType(LAMMPS *lmp, int narg, char **arg) : Compute( if (strcmp(arg[3],"atom") == 0) mode = ATOM; else if (strcmp(arg[3],"bond") == 0) mode = BOND; + else if (strcmp(arg[3],"angle") == 0) mode = ANGLE; + else if (strcmp(arg[3],"dihedral") == 0) mode = DIHEDRAL; + else if (strcmp(arg[3],"improper") == 0) mode = IMPROPER; else error->all(FLERR, "Invalid compute count/type keyword {}",arg[3]); + // error check + + if (mode == BOND && !atom->nbondtypes) + error->all(FLERR,"Compute count/type bond command with no bonds defined"); + if (mode == ANGLE && !atom->nangletypes) + error->all(FLERR,"Compute count/type bond command with no angles defined"); + if (mode == DIHEDRAL && !atom->ndihedraltypes) + error->all(FLERR,"Compute count/type dihedral command with no dihedrals defined"); + if (mode == IMPROPER && !atom->nimpropertypes) + error->all(FLERR,"Compute count/type improper command with no impropers defined"); + + // set vector lengths + if (mode == ATOM) { vector_flag = 1; size_vector = atom->ntypes; @@ -46,11 +61,20 @@ ComputeCountType::ComputeCountType(LAMMPS *lmp, int narg, char **arg) : Compute( size_vector = atom->nbondtypes; extscalar = 1; extvector = 1; + } else if (mode == ANGLE) { + vector_flag = 1; + size_vector = atom->nangletypes; + extvector = 1; + } else if (mode == DIHEDRAL) { + vector_flag = 1; + size_vector = atom->ndihedraltypes; + extvector = 1; + } else if (mode == IMPROPER) { + vector_flag = 1; + size_vector = atom->nimpropertypes; + extvector = 1; } - if (mode == BOND && !atom->avec->bonds_allow) - error->all(FLERR,"Cannot use compute count/type bond command with no bonds allowed"); - // output vector vector = new double[size_vector]; @@ -121,76 +145,246 @@ void ComputeCountType::compute_vector() int nvec; - // count atoms by type - // atom must be in group to be counted - - if (mode == ATOM) { - int *type = atom->type; - int *mask = atom->mask; - int nlocal = atom->nlocal; - int ntypes = atom->ntypes; - - for (int m = 0; m < ntypes; m++) count[m] = 0; - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) - count[type[i]-1]++; - - nvec = ntypes; - } + if (mode == ATOM) nvec = count_atoms(); + else if (mode == BOND) nvec = count_bonds(); + else if (mode == ANGLE) nvec = count_angles(); + else if (mode == DIHEDRAL) nvec = count_dihedrals(); + else if (mode == IMPROPER) nvec = count_impropers(); - // count bonds by type - // both atoms in bond must be in group to be counted - // skip type = 0 bonds, they are counted by compute_scalar - // bond types can be negative for SHAKE - - else if (mode == BOND) { - int **bond_atom = atom->bond_atom; - int **bond_type = atom->bond_type; - int *num_bond = atom->num_bond; - int *mask = atom->mask; - int nlocal = atom->nlocal; - int nbondtypes = atom->nbondtypes; - - int j,m,nbond,itype; - int flag = 0; - for (int m = 0; m < nbondtypes; m++) count[m] = 0; - - for (int i = 0; i < nlocal; i++) { - nbond = num_bond[i]; - for (m = 0; m < nbond; m++) { - itype = bond_type[i][m]; - if (itype == 0) continue; - - j = atom->map(bond_atom[i][m]); - if (j < 0) { - flag = 1; - continue; - } - - if ((mask[i] & groupbit) && (mask[j] & groupbit)) { - if (itype > 0) count[itype-1]++; - else count[-itype-1]++; - } - } - } - - int flagany; - MPI_Allreduce(&flag, &flagany, 1, MPI_INT, MPI_SUM, world); - if (flagany) error->all(FLERR,"Missing bond atom in compute count/type"); - - nvec = nbondtypes; - } - // sum across procs as bigint, then convert to double - // correct for double counting if newton_bond off + // correct for multiple counting if newton_bond off for (int m = 0; m < nvec; m++) bcount_me[m] = count[m]; MPI_Allreduce(bcount_me, bcount, nvec, MPI_LMP_BIGINT, MPI_SUM, world); - if (force->newton_bond == 0) - for (int m = 0; m < nvec; m++) bcount[m] /= 2; - + + if (force->newton_bond == 0) { + if (mode == BOND) + for (int m = 0; m < nvec; m++) bcount[m] /= 2; + else if (mode == ANGLE) + for (int m = 0; m < nvec; m++) bcount[m] /= 3; + if (mode == DIHEDRAL || mode == IMPROPER) + for (int m = 0; m < nvec; m++) bcount[m] /= 4; + } + for (int m = 0; m < nvec; m++) if (bcount[m] > MAXDOUBLEINT) error->all(FLERR,"Compute count/type overflow"); for (int m = 0; m < nvec; m++) vector[m] = bcount[m]; } + +/* ---------------------------------------------------------------------- + count atoms by type + atom must be in group to be counted +---------------------------------------------------------------------- */ + +int ComputeCountType::count_atoms() +{ + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + int ntypes = atom->ntypes; + + for (int m = 0; m < ntypes; m++) count[m] = 0; + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + count[type[i]-1]++; + + return ntypes; +} + +/* ---------------------------------------------------------------------- + count bonds by type + both atoms in bond must be in group to be counted + skip type = 0 bonds, they are counted by compute_scalar() + bond types can be negative, count them as if positive +---------------------------------------------------------------------- */ + +int ComputeCountType::count_bonds() +{ + int **bond_atom = atom->bond_atom; + int **bond_type = atom->bond_type; + int *num_bond = atom->num_bond; + int *mask = atom->mask; + int nlocal = atom->nlocal; + int nbondtypes = atom->nbondtypes; + + int j,m,nbond,itype; + int flag = 0; + for (int m = 0; m < nbondtypes; m++) count[m] = 0; + + for (int i = 0; i < nlocal; i++) { + nbond = num_bond[i]; + for (m = 0; m < nbond; m++) { + itype = bond_type[i][m]; + if (itype == 0) continue; + + j = atom->map(bond_atom[i][m]); + if (j < 0) { + flag = 1; + continue; + } + + if ((mask[i] & groupbit) && (mask[j] & groupbit)) { + if (itype > 0) count[itype-1]++; + else count[-itype-1]++; + } + } + } + + int flagany; + MPI_Allreduce(&flag, &flagany, 1, MPI_INT, MPI_SUM, world); + if (flagany) error->all(FLERR,"Missing bond atom in compute count/type"); + + return nbondtypes; +} + +/* ---------------------------------------------------------------------- + count angles by type + all 3 atoms in angle must be in group to be counted + angle types can be negative, count them as if positive +---------------------------------------------------------------------- */ + +int ComputeCountType::count_angles() +{ + int **angle_atom1 = atom->angle_atom1; + int **angle_atom2 = atom->angle_atom2; + int **angle_atom3 = atom->angle_atom3; + int **angle_type = atom->angle_type; + int *num_angle = atom->num_angle; + int *mask = atom->mask; + int nlocal = atom->nlocal; + int nangletypes = atom->nangletypes; + + int j1,j2,j3,m,nangle,itype; + int flag = 0; + for (int m = 0; m < nangletypes; m++) count[m] = 0; + + for (int i = 0; i < nlocal; i++) { + nangle = num_angle[i]; + for (m = 0; m < nangle; m++) { + itype = angle_type[i][m]; + + j1 = atom->map(angle_atom1[i][m]); + j2 = atom->map(angle_atom2[i][m]); + j3 = atom->map(angle_atom3[i][m]); + if (j1 < 0 || j2 < 0 || j3 < 0) { + flag = 1; + continue; + } + + if ((mask[j1] & groupbit) && (mask[j2] & groupbit) && + (mask[j3] & groupbit)) { + if (itype > 0) count[itype-1]++; + else count[-itype-1]++; + } + } + } + + int flagany; + MPI_Allreduce(&flag, &flagany, 1, MPI_INT, MPI_SUM, world); + if (flagany) error->all(FLERR,"Missing angle atom in compute count/type"); + + return nangletypes; +} + +/* ---------------------------------------------------------------------- + count dihedrals by type + all 4 atoms in dihedral must be in group to be counted + dihedral types can be negative, count them as if positive +---------------------------------------------------------------------- */ + +int ComputeCountType::count_dihedrals() +{ + int **dihedral_atom1 = atom->dihedral_atom1; + int **dihedral_atom2 = atom->dihedral_atom2; + int **dihedral_atom3 = atom->dihedral_atom3; + int **dihedral_atom4 = atom->dihedral_atom4; + int **dihedral_type = atom->dihedral_type; + int *num_dihedral = atom->num_dihedral; + int *mask = atom->mask; + int nlocal = atom->nlocal; + int ndihedraltypes = atom->ndihedraltypes; + + int j1,j2,j3,j4,m,ndihedral,itype; + int flag = 0; + for (int m = 0; m < ndihedraltypes; m++) count[m] = 0; + + for (int i = 0; i < nlocal; i++) { + ndihedral = num_dihedral[i]; + for (m = 0; m < ndihedral; m++) { + itype = dihedral_type[i][m]; + + j1 = atom->map(dihedral_atom1[i][m]); + j2 = atom->map(dihedral_atom2[i][m]); + j3 = atom->map(dihedral_atom3[i][m]); + j4 = atom->map(dihedral_atom4[i][m]); + if (j1 < 0 || j2 < 0 || j3 < 0 || j4 < 0) { + flag = 1; + continue; + } + + if ((mask[j1] & groupbit) && (mask[j2] & groupbit) && + (mask[j3] & groupbit) && (mask[j4] & groupbit)) { + if (itype > 0) count[itype-1]++; + else count[-itype-1]++; + } + } + } + + int flagany; + MPI_Allreduce(&flag, &flagany, 1, MPI_INT, MPI_SUM, world); + if (flagany) error->all(FLERR,"Missing dihedral atom in compute count/type"); + + return ndihedraltypes; +} + +/* ---------------------------------------------------------------------- + count impropers by type + all 4 atoms in improper must be in group to be counted + improper types can be negative, count them as if positive +---------------------------------------------------------------------- */ + +int ComputeCountType::count_impropers() +{ + int **improper_atom1 = atom->improper_atom1; + int **improper_atom2 = atom->improper_atom2; + int **improper_atom3 = atom->improper_atom3; + int **improper_atom4 = atom->improper_atom4; + int **improper_type = atom->improper_type; + int *num_improper = atom->num_improper; + int *mask = atom->mask; + int nlocal = atom->nlocal; + int nimpropertypes = atom->nimpropertypes; + + int j1,j2,j3,j4,m,nimproper,itype; + int flag = 0; + for (int m = 0; m < nimpropertypes; m++) count[m] = 0; + + for (int i = 0; i < nlocal; i++) { + nimproper = num_improper[i]; + for (m = 0; m < nimproper; m++) { + itype = improper_type[i][m]; + + j1 = atom->map(improper_atom1[i][m]); + j2 = atom->map(improper_atom2[i][m]); + j3 = atom->map(improper_atom3[i][m]); + j4 = atom->map(improper_atom4[i][m]); + if (j1 < 0 || j2 < 0 || j3 < 0 || j4 < 0) { + flag = 1; + continue; + } + + if ((mask[j1] & groupbit) && (mask[j2] & groupbit) && + (mask[j3] & groupbit) && (mask[j4] & groupbit)) { + if (itype > 0) count[itype-1]++; + else count[-itype-1]++; + } + } + } + + int flagany; + MPI_Allreduce(&flag, &flagany, 1, MPI_INT, MPI_SUM, world); + if (flagany) error->all(FLERR,"Missing improper atom in compute count/type"); + + return nimpropertypes; +} diff --git a/src/compute_count_type.h b/src/compute_count_type.h index 0ed7a22467..68a51cb42a 100644 --- a/src/compute_count_type.h +++ b/src/compute_count_type.h @@ -38,6 +38,12 @@ class ComputeCountType : public Compute { int *count; bigint *bcount_me; bigint *bcount; + + int count_atoms(); + int count_bonds(); + int count_angles(); + int count_dihedrals(); + int count_impropers(); }; } // namespace LAMMPS_NS From 8a9091595df80f2834cb7452e8d1cc0a8cd1a7b0 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 2 May 2023 15:56:45 -0400 Subject: [PATCH 15/23] Revert "use more obvious representation of 2^53" This reverts commit ff29ef7d3160c9597b9fb637619982afcf8e0adc as it conflicts with how MSVC interprets 1L (it would require 1LL which creates issues with other compilers). --- src/lmptype.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/lmptype.h b/src/lmptype.h index 769f730d6a..891d5bdf89 100644 --- a/src/lmptype.h +++ b/src/lmptype.h @@ -95,8 +95,8 @@ typedef int64_t bigint; #define MAXSMALLINT INT_MAX #define MAXTAGINT INT_MAX #define MAXBIGINT INT64_MAX -#define MAXDOUBLEINT (1L<<53) - +#define MAXDOUBLEINT 9007199254740992 // 2^53 + #define MPI_LMP_TAGINT MPI_INT #define MPI_LMP_IMAGEINT MPI_INT #define MPI_LMP_BIGINT MPI_LL @@ -133,7 +133,7 @@ typedef int64_t bigint; #define MAXSMALLINT INT_MAX #define MAXTAGINT INT64_MAX #define MAXBIGINT INT64_MAX -#define MAXDOUBLEINT (1L<<53) +#define MAXDOUBLEINT 9007199254740992 // 2^53 #define MPI_LMP_TAGINT MPI_LL #define MPI_LMP_IMAGEINT MPI_LL From 3119434932f05b1276ece2c7d17bfbd1076f44fd Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 2 May 2023 15:59:34 -0400 Subject: [PATCH 16/23] apply clang-format --- src/compute_count_type.cpp | 147 +++++++++++++++++++++---------------- 1 file changed, 82 insertions(+), 65 deletions(-) diff --git a/src/compute_count_type.cpp b/src/compute_count_type.cpp index ed0467c435..5c27dd39ca 100644 --- a/src/compute_count_type.cpp +++ b/src/compute_count_type.cpp @@ -22,7 +22,7 @@ using namespace LAMMPS_NS; -enum{ATOM, BOND, ANGLE, DIHEDRAL, IMPROPER}; +enum { ATOM, BOND, ANGLE, DIHEDRAL, IMPROPER }; /* ---------------------------------------------------------------------- */ @@ -33,23 +33,29 @@ ComputeCountType::ComputeCountType(LAMMPS *lmp, int narg, char **arg) : // process args - if (strcmp(arg[3],"atom") == 0) mode = ATOM; - else if (strcmp(arg[3],"bond") == 0) mode = BOND; - else if (strcmp(arg[3],"angle") == 0) mode = ANGLE; - else if (strcmp(arg[3],"dihedral") == 0) mode = DIHEDRAL; - else if (strcmp(arg[3],"improper") == 0) mode = IMPROPER; - else error->all(FLERR, "Invalid compute count/type keyword {}",arg[3]); + if (strcmp(arg[3], "atom") == 0) + mode = ATOM; + else if (strcmp(arg[3], "bond") == 0) + mode = BOND; + else if (strcmp(arg[3], "angle") == 0) + mode = ANGLE; + else if (strcmp(arg[3], "dihedral") == 0) + mode = DIHEDRAL; + else if (strcmp(arg[3], "improper") == 0) + mode = IMPROPER; + else + error->all(FLERR, "Invalid compute count/type keyword {}", arg[3]); // error check - + if (mode == BOND && !atom->nbondtypes) - error->all(FLERR,"Compute count/type bond command with no bonds defined"); + error->all(FLERR, "Compute count/type bond command with no bonds defined"); if (mode == ANGLE && !atom->nangletypes) - error->all(FLERR,"Compute count/type bond command with no angles defined"); + error->all(FLERR, "Compute count/type bond command with no angles defined"); if (mode == DIHEDRAL && !atom->ndihedraltypes) - error->all(FLERR,"Compute count/type dihedral command with no dihedrals defined"); + error->all(FLERR, "Compute count/type dihedral command with no dihedrals defined"); if (mode == IMPROPER && !atom->nimpropertypes) - error->all(FLERR,"Compute count/type improper command with no impropers defined"); + error->all(FLERR, "Compute count/type improper command with no impropers defined"); // set vector lengths @@ -145,18 +151,23 @@ void ComputeCountType::compute_vector() int nvec; - if (mode == ATOM) nvec = count_atoms(); - else if (mode == BOND) nvec = count_bonds(); - else if (mode == ANGLE) nvec = count_angles(); - else if (mode == DIHEDRAL) nvec = count_dihedrals(); - else if (mode == IMPROPER) nvec = count_impropers(); - + if (mode == ATOM) + nvec = count_atoms(); + else if (mode == BOND) + nvec = count_bonds(); + else if (mode == ANGLE) + nvec = count_angles(); + else if (mode == DIHEDRAL) + nvec = count_dihedrals(); + else if (mode == IMPROPER) + nvec = count_impropers(); + // sum across procs as bigint, then convert to double // correct for multiple counting if newton_bond off for (int m = 0; m < nvec; m++) bcount_me[m] = count[m]; MPI_Allreduce(bcount_me, bcount, nvec, MPI_LMP_BIGINT, MPI_SUM, world); - + if (force->newton_bond == 0) { if (mode == BOND) for (int m = 0; m < nvec; m++) bcount[m] /= 2; @@ -165,7 +176,7 @@ void ComputeCountType::compute_vector() if (mode == DIHEDRAL || mode == IMPROPER) for (int m = 0; m < nvec; m++) bcount[m] /= 4; } - + for (int m = 0; m < nvec; m++) if (bcount[m] > MAXDOUBLEINT) error->all(FLERR, "Compute count/type overflow"); for (int m = 0; m < nvec; m++) vector[m] = bcount[m]; @@ -182,11 +193,10 @@ int ComputeCountType::count_atoms() int *mask = atom->mask; int nlocal = atom->nlocal; int ntypes = atom->ntypes; - + for (int m = 0; m < ntypes; m++) count[m] = 0; for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) - count[type[i]-1]++; + if (mask[i] & groupbit) count[type[i] - 1]++; return ntypes; } @@ -206,34 +216,36 @@ int ComputeCountType::count_bonds() int *mask = atom->mask; int nlocal = atom->nlocal; int nbondtypes = atom->nbondtypes; - - int j,m,nbond,itype; + + int j, m, nbond, itype; int flag = 0; for (int m = 0; m < nbondtypes; m++) count[m] = 0; - + for (int i = 0; i < nlocal; i++) { nbond = num_bond[i]; for (m = 0; m < nbond; m++) { itype = bond_type[i][m]; if (itype == 0) continue; - + j = atom->map(bond_atom[i][m]); if (j < 0) { flag = 1; continue; } - + if ((mask[i] & groupbit) && (mask[j] & groupbit)) { - if (itype > 0) count[itype-1]++; - else count[-itype-1]++; + if (itype > 0) + count[itype - 1]++; + else + count[-itype - 1]++; } } } int flagany; MPI_Allreduce(&flag, &flagany, 1, MPI_INT, MPI_SUM, world); - if (flagany) error->all(FLERR,"Missing bond atom in compute count/type"); - + if (flagany) error->all(FLERR, "Missing bond atom in compute count/type"); + return nbondtypes; } @@ -253,16 +265,16 @@ int ComputeCountType::count_angles() int *mask = atom->mask; int nlocal = atom->nlocal; int nangletypes = atom->nangletypes; - - int j1,j2,j3,m,nangle,itype; + + int j1, j2, j3, m, nangle, itype; int flag = 0; for (int m = 0; m < nangletypes; m++) count[m] = 0; - + for (int i = 0; i < nlocal; i++) { nangle = num_angle[i]; for (m = 0; m < nangle; m++) { itype = angle_type[i][m]; - + j1 = atom->map(angle_atom1[i][m]); j2 = atom->map(angle_atom2[i][m]); j3 = atom->map(angle_atom3[i][m]); @@ -270,19 +282,20 @@ int ComputeCountType::count_angles() flag = 1; continue; } - - if ((mask[j1] & groupbit) && (mask[j2] & groupbit) && - (mask[j3] & groupbit)) { - if (itype > 0) count[itype-1]++; - else count[-itype-1]++; + + if ((mask[j1] & groupbit) && (mask[j2] & groupbit) && (mask[j3] & groupbit)) { + if (itype > 0) + count[itype - 1]++; + else + count[-itype - 1]++; } } } int flagany; MPI_Allreduce(&flag, &flagany, 1, MPI_INT, MPI_SUM, world); - if (flagany) error->all(FLERR,"Missing angle atom in compute count/type"); - + if (flagany) error->all(FLERR, "Missing angle atom in compute count/type"); + return nangletypes; } @@ -303,16 +316,16 @@ int ComputeCountType::count_dihedrals() int *mask = atom->mask; int nlocal = atom->nlocal; int ndihedraltypes = atom->ndihedraltypes; - - int j1,j2,j3,j4,m,ndihedral,itype; + + int j1, j2, j3, j4, m, ndihedral, itype; int flag = 0; for (int m = 0; m < ndihedraltypes; m++) count[m] = 0; - + for (int i = 0; i < nlocal; i++) { ndihedral = num_dihedral[i]; for (m = 0; m < ndihedral; m++) { itype = dihedral_type[i][m]; - + j1 = atom->map(dihedral_atom1[i][m]); j2 = atom->map(dihedral_atom2[i][m]); j3 = atom->map(dihedral_atom3[i][m]); @@ -321,19 +334,21 @@ int ComputeCountType::count_dihedrals() flag = 1; continue; } - - if ((mask[j1] & groupbit) && (mask[j2] & groupbit) && - (mask[j3] & groupbit) && (mask[j4] & groupbit)) { - if (itype > 0) count[itype-1]++; - else count[-itype-1]++; + + if ((mask[j1] & groupbit) && (mask[j2] & groupbit) && (mask[j3] & groupbit) && + (mask[j4] & groupbit)) { + if (itype > 0) + count[itype - 1]++; + else + count[-itype - 1]++; } } } int flagany; MPI_Allreduce(&flag, &flagany, 1, MPI_INT, MPI_SUM, world); - if (flagany) error->all(FLERR,"Missing dihedral atom in compute count/type"); - + if (flagany) error->all(FLERR, "Missing dihedral atom in compute count/type"); + return ndihedraltypes; } @@ -354,16 +369,16 @@ int ComputeCountType::count_impropers() int *mask = atom->mask; int nlocal = atom->nlocal; int nimpropertypes = atom->nimpropertypes; - - int j1,j2,j3,j4,m,nimproper,itype; + + int j1, j2, j3, j4, m, nimproper, itype; int flag = 0; for (int m = 0; m < nimpropertypes; m++) count[m] = 0; - + for (int i = 0; i < nlocal; i++) { nimproper = num_improper[i]; for (m = 0; m < nimproper; m++) { itype = improper_type[i][m]; - + j1 = atom->map(improper_atom1[i][m]); j2 = atom->map(improper_atom2[i][m]); j3 = atom->map(improper_atom3[i][m]); @@ -372,18 +387,20 @@ int ComputeCountType::count_impropers() flag = 1; continue; } - - if ((mask[j1] & groupbit) && (mask[j2] & groupbit) && - (mask[j3] & groupbit) && (mask[j4] & groupbit)) { - if (itype > 0) count[itype-1]++; - else count[-itype-1]++; + + if ((mask[j1] & groupbit) && (mask[j2] & groupbit) && (mask[j3] & groupbit) && + (mask[j4] & groupbit)) { + if (itype > 0) + count[itype - 1]++; + else + count[-itype - 1]++; } } } int flagany; MPI_Allreduce(&flag, &flagany, 1, MPI_INT, MPI_SUM, world); - if (flagany) error->all(FLERR,"Missing improper atom in compute count/type"); - + if (flagany) error->all(FLERR, "Missing improper atom in compute count/type"); + return nimpropertypes; } From e8a77c61ac72a9933f9d10d68026b5e0e9901214 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 2 May 2023 16:19:53 -0400 Subject: [PATCH 17/23] whitespace --- doc/src/compute_count_type.rst | 2 +- src/lmptype.h | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/src/compute_count_type.rst b/doc/src/compute_count_type.rst index 2e3098b5ea..bb276ea21b 100644 --- a/doc/src/compute_count_type.rst +++ b/doc/src/compute_count_type.rst @@ -50,7 +50,7 @@ type: * :doc:`fix shake ` * :doc:`delete_bonds ` - + These commands can create and/or break topological bonds (angles, etc). In the case of breaking, they remove the bond (angle, etc) from the system, so that they no longer exist (:doc:`bond_style quartic diff --git a/src/lmptype.h b/src/lmptype.h index 891d5bdf89..ecd6ee5761 100644 --- a/src/lmptype.h +++ b/src/lmptype.h @@ -96,7 +96,7 @@ typedef int64_t bigint; #define MAXTAGINT INT_MAX #define MAXBIGINT INT64_MAX #define MAXDOUBLEINT 9007199254740992 // 2^53 - + #define MPI_LMP_TAGINT MPI_INT #define MPI_LMP_IMAGEINT MPI_INT #define MPI_LMP_BIGINT MPI_LL From b6e211dd62147c4f35ee70396f2831dfc2613492 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 2 May 2023 16:20:20 -0400 Subject: [PATCH 18/23] consistently declare variables when used only --- src/compute_count_type.cpp | 64 +++++++++++++++++--------------------- 1 file changed, 29 insertions(+), 35 deletions(-) diff --git a/src/compute_count_type.cpp b/src/compute_count_type.cpp index 5c27dd39ca..47bbe36dd2 100644 --- a/src/compute_count_type.cpp +++ b/src/compute_count_type.cpp @@ -27,7 +27,7 @@ enum { ATOM, BOND, ANGLE, DIHEDRAL, IMPROPER }; /* ---------------------------------------------------------------------- */ ComputeCountType::ComputeCountType(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), count(nullptr), bcount(nullptr), bcount_me(nullptr) + Compute(lmp, narg, arg), count(nullptr), bcount_me(nullptr), bcount(nullptr) { if (narg != 4) error->all(FLERR, "Incorrect number of args for compute count/type command"); @@ -117,24 +117,22 @@ double ComputeCountType::compute_scalar() int **bond_type = atom->bond_type; int nlocal = atom->nlocal; - int m, nbond; - int count = 0; - // count broken bonds with bond_type = 0 // ignore group setting since 2 atoms in a broken bond // can be arbitrarily far apart + int count = 0; for (int i = 0; i < nlocal; i++) { - nbond = num_bond[i]; - for (m = 0; m < nbond; m++) + int nbond = num_bond[i]; + for (int m = 0; m < nbond; m++) if (bond_type[i][m] == 0) count++; } // sum across procs as bigint, then convert to double // correct for double counting if newton_bond off + bigint bcount = 0; bigint bcount_me = count; - bigint bcount; MPI_Allreduce(&bcount_me, &bcount, 1, MPI_LMP_BIGINT, MPI_SUM, world); if (force->newton_bond == 0) bcount /= 2; @@ -217,17 +215,16 @@ int ComputeCountType::count_bonds() int nlocal = atom->nlocal; int nbondtypes = atom->nbondtypes; - int j, m, nbond, itype; int flag = 0; for (int m = 0; m < nbondtypes; m++) count[m] = 0; for (int i = 0; i < nlocal; i++) { - nbond = num_bond[i]; - for (m = 0; m < nbond; m++) { - itype = bond_type[i][m]; + int nbond = num_bond[i]; + for (int m = 0; m < nbond; m++) { + int itype = bond_type[i][m]; if (itype == 0) continue; - j = atom->map(bond_atom[i][m]); + int j = atom->map(bond_atom[i][m]); if (j < 0) { flag = 1; continue; @@ -266,18 +263,17 @@ int ComputeCountType::count_angles() int nlocal = atom->nlocal; int nangletypes = atom->nangletypes; - int j1, j2, j3, m, nangle, itype; int flag = 0; for (int m = 0; m < nangletypes; m++) count[m] = 0; for (int i = 0; i < nlocal; i++) { - nangle = num_angle[i]; - for (m = 0; m < nangle; m++) { - itype = angle_type[i][m]; + int nangle = num_angle[i]; + for (int m = 0; m < nangle; m++) { + int itype = angle_type[i][m]; - j1 = atom->map(angle_atom1[i][m]); - j2 = atom->map(angle_atom2[i][m]); - j3 = atom->map(angle_atom3[i][m]); + int j1 = atom->map(angle_atom1[i][m]); + int j2 = atom->map(angle_atom2[i][m]); + int j3 = atom->map(angle_atom3[i][m]); if (j1 < 0 || j2 < 0 || j3 < 0) { flag = 1; continue; @@ -317,19 +313,18 @@ int ComputeCountType::count_dihedrals() int nlocal = atom->nlocal; int ndihedraltypes = atom->ndihedraltypes; - int j1, j2, j3, j4, m, ndihedral, itype; int flag = 0; for (int m = 0; m < ndihedraltypes; m++) count[m] = 0; for (int i = 0; i < nlocal; i++) { - ndihedral = num_dihedral[i]; - for (m = 0; m < ndihedral; m++) { - itype = dihedral_type[i][m]; + int ndihedral = num_dihedral[i]; + for (int m = 0; m < ndihedral; m++) { + int itype = dihedral_type[i][m]; - j1 = atom->map(dihedral_atom1[i][m]); - j2 = atom->map(dihedral_atom2[i][m]); - j3 = atom->map(dihedral_atom3[i][m]); - j4 = atom->map(dihedral_atom4[i][m]); + int j1 = atom->map(dihedral_atom1[i][m]); + int j2 = atom->map(dihedral_atom2[i][m]); + int j3 = atom->map(dihedral_atom3[i][m]); + int j4 = atom->map(dihedral_atom4[i][m]); if (j1 < 0 || j2 < 0 || j3 < 0 || j4 < 0) { flag = 1; continue; @@ -370,19 +365,18 @@ int ComputeCountType::count_impropers() int nlocal = atom->nlocal; int nimpropertypes = atom->nimpropertypes; - int j1, j2, j3, j4, m, nimproper, itype; int flag = 0; for (int m = 0; m < nimpropertypes; m++) count[m] = 0; for (int i = 0; i < nlocal; i++) { - nimproper = num_improper[i]; - for (m = 0; m < nimproper; m++) { - itype = improper_type[i][m]; + int nimproper = num_improper[i]; + for (int m = 0; m < nimproper; m++) { + int itype = improper_type[i][m]; - j1 = atom->map(improper_atom1[i][m]); - j2 = atom->map(improper_atom2[i][m]); - j3 = atom->map(improper_atom3[i][m]); - j4 = atom->map(improper_atom4[i][m]); + int j1 = atom->map(improper_atom1[i][m]); + int j2 = atom->map(improper_atom2[i][m]); + int j3 = atom->map(improper_atom3[i][m]); + int j4 = atom->map(improper_atom4[i][m]); if (j1 < 0 || j2 < 0 || j3 < 0 || j4 < 0) { flag = 1; continue; From 3f8cd4577cd161a7ac9a04781bcf2882336a46fb Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 2 May 2023 16:20:29 -0400 Subject: [PATCH 19/23] update unit test --- unittest/commands/test_compute_global.cpp | 66 ++++++++++++++++++----- 1 file changed, 52 insertions(+), 14 deletions(-) diff --git a/unittest/commands/test_compute_global.cpp b/unittest/commands/test_compute_global.cpp index 7e5b3d4483..6c365c8c2b 100644 --- a/unittest/commands/test_compute_global.cpp +++ b/unittest/commands/test_compute_global.cpp @@ -308,45 +308,83 @@ TEST_F(ComputeGlobalTest, Counts) command("variable t3 atom type==3"); command("variable t4 atom type==4"); command("variable t5 atom type==5"); - command("compute asum all reduce sum v_t1 v_t2 v_t3 v_t4 v_t5"); - command("compute acnt all count/type atom"); + command("compute tsum all reduce sum v_t1 v_t2 v_t3 v_t4 v_t5"); + command("compute tcnt all count/type atom"); command("compute bcnt all count/type bond"); - command("thermo_style custom c_asum[*] c_acnt[*]"); + command("compute acnt all count/type angle"); + command("compute dcnt all count/type dihedral"); + command("compute icnt all count/type improper"); + command("thermo_style custom c_tsum[*] c_tcnt[*] c_bcnt[*] c_acnt[*] c_dcnt[*] c_icnt[*]"); command("run 0 post no"); END_HIDE_OUTPUT(); - auto asum = get_vector("asum"); - auto acnt = get_vector("acnt"); + auto tsum = get_vector("tsum"); + auto tcnt = get_vector("tcnt"); auto bcnt = get_vector("bcnt"); auto bbrk = get_scalar("bcnt"); + auto acnt = get_vector("acnt"); + auto dcnt = get_vector("dcnt"); + auto icnt = get_vector("icnt"); + + EXPECT_DOUBLE_EQ(tsum[0], tcnt[0]); + EXPECT_DOUBLE_EQ(tsum[1], tcnt[1]); + EXPECT_DOUBLE_EQ(tsum[2], tcnt[2]); + EXPECT_DOUBLE_EQ(tsum[3], tcnt[3]); + EXPECT_DOUBLE_EQ(tsum[4], tcnt[4]); EXPECT_DOUBLE_EQ(bbrk, 0.0); - EXPECT_DOUBLE_EQ(asum[0], acnt[0]); - EXPECT_DOUBLE_EQ(asum[1], acnt[1]); - EXPECT_DOUBLE_EQ(asum[2], acnt[2]); - EXPECT_DOUBLE_EQ(asum[3], acnt[3]); - EXPECT_DOUBLE_EQ(asum[4], acnt[4]); - EXPECT_DOUBLE_EQ(bcnt[0], 3.0); EXPECT_DOUBLE_EQ(bcnt[1], 6.0); EXPECT_DOUBLE_EQ(bcnt[2], 3.0); EXPECT_DOUBLE_EQ(bcnt[3], 2.0); EXPECT_DOUBLE_EQ(bcnt[4], 10.0); + EXPECT_DOUBLE_EQ(acnt[0], 6.0); + EXPECT_DOUBLE_EQ(acnt[1], 10.0); + EXPECT_DOUBLE_EQ(acnt[2], 5.0); + EXPECT_DOUBLE_EQ(acnt[3], 9.0); + + EXPECT_DOUBLE_EQ(dcnt[0], 3.0); + EXPECT_DOUBLE_EQ(dcnt[1], 8.0); + EXPECT_DOUBLE_EQ(dcnt[2], 3.0); + EXPECT_DOUBLE_EQ(dcnt[3], 4.0); + EXPECT_DOUBLE_EQ(dcnt[4], 13.0); + + EXPECT_DOUBLE_EQ(icnt[0], 1.0); + EXPECT_DOUBLE_EQ(icnt[1], 1.0); + BEGIN_HIDE_OUTPUT(); - command("delete_bonds all bond 3"); + command("delete_bonds all bond 3 remove"); command("run 0 post no"); END_HIDE_OUTPUT(); bcnt = get_vector("bcnt"); bbrk = get_scalar("bcnt"); - EXPECT_DOUBLE_EQ(bbrk, 0.0); // should be 3 + acnt = get_vector("acnt"); + dcnt = get_vector("dcnt"); + icnt = get_vector("icnt"); + + EXPECT_DOUBLE_EQ(bbrk, 0.0); EXPECT_DOUBLE_EQ(bcnt[0], 3.0); EXPECT_DOUBLE_EQ(bcnt[1], 6.0); - EXPECT_DOUBLE_EQ(bcnt[2], 3.0); // should be 0 + EXPECT_DOUBLE_EQ(bcnt[2], 0.0); EXPECT_DOUBLE_EQ(bcnt[3], 2.0); EXPECT_DOUBLE_EQ(bcnt[4], 10.0); + + EXPECT_DOUBLE_EQ(acnt[0], 6.0); + EXPECT_DOUBLE_EQ(acnt[1], 10.0); + EXPECT_DOUBLE_EQ(acnt[2], 5.0); + EXPECT_DOUBLE_EQ(acnt[3], 9.0); + + EXPECT_DOUBLE_EQ(dcnt[0], 3.0); + EXPECT_DOUBLE_EQ(dcnt[1], 8.0); + EXPECT_DOUBLE_EQ(dcnt[2], 3.0); + EXPECT_DOUBLE_EQ(dcnt[3], 4.0); + EXPECT_DOUBLE_EQ(dcnt[4], 13.0); + + EXPECT_DOUBLE_EQ(icnt[0], 1.0); + EXPECT_DOUBLE_EQ(icnt[1], 1.0); } } // namespace LAMMPS_NS From 6d050374c34b18ea65b06fe06feed36834d61e35 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 2 May 2023 16:24:43 -0400 Subject: [PATCH 20/23] fix compilation with -DLAMMPS_BIGBIG --- src/compute_count_type.cpp | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/src/compute_count_type.cpp b/src/compute_count_type.cpp index 47bbe36dd2..84b31e8b85 100644 --- a/src/compute_count_type.cpp +++ b/src/compute_count_type.cpp @@ -208,7 +208,7 @@ int ComputeCountType::count_atoms() int ComputeCountType::count_bonds() { - int **bond_atom = atom->bond_atom; + tagint **bond_atom = atom->bond_atom; int **bond_type = atom->bond_type; int *num_bond = atom->num_bond; int *mask = atom->mask; @@ -254,9 +254,9 @@ int ComputeCountType::count_bonds() int ComputeCountType::count_angles() { - int **angle_atom1 = atom->angle_atom1; - int **angle_atom2 = atom->angle_atom2; - int **angle_atom3 = atom->angle_atom3; + tagint **angle_atom1 = atom->angle_atom1; + tagint **angle_atom2 = atom->angle_atom2; + tagint **angle_atom3 = atom->angle_atom3; int **angle_type = atom->angle_type; int *num_angle = atom->num_angle; int *mask = atom->mask; @@ -303,10 +303,10 @@ int ComputeCountType::count_angles() int ComputeCountType::count_dihedrals() { - int **dihedral_atom1 = atom->dihedral_atom1; - int **dihedral_atom2 = atom->dihedral_atom2; - int **dihedral_atom3 = atom->dihedral_atom3; - int **dihedral_atom4 = atom->dihedral_atom4; + tagint **dihedral_atom1 = atom->dihedral_atom1; + tagint **dihedral_atom2 = atom->dihedral_atom2; + tagint **dihedral_atom3 = atom->dihedral_atom3; + tagint **dihedral_atom4 = atom->dihedral_atom4; int **dihedral_type = atom->dihedral_type; int *num_dihedral = atom->num_dihedral; int *mask = atom->mask; @@ -355,10 +355,10 @@ int ComputeCountType::count_dihedrals() int ComputeCountType::count_impropers() { - int **improper_atom1 = atom->improper_atom1; - int **improper_atom2 = atom->improper_atom2; - int **improper_atom3 = atom->improper_atom3; - int **improper_atom4 = atom->improper_atom4; + tagint **improper_atom1 = atom->improper_atom1; + tagint **improper_atom2 = atom->improper_atom2; + tagint **improper_atom3 = atom->improper_atom3; + tagint **improper_atom4 = atom->improper_atom4; int **improper_type = atom->improper_type; int *num_improper = atom->num_improper; int *mask = atom->mask; From 36632e33566bd6a93f08539461d5051bd71a1097 Mon Sep 17 00:00:00 2001 From: Joel Thomas Clemmer Date: Tue, 2 May 2023 14:28:17 -0600 Subject: [PATCH 21/23] Clarifying how BPM package works --- doc/src/Howto_broken_bonds.rst | 14 ++++++++------ doc/src/compute_count_type.rst | 14 ++++++++------ 2 files changed, 16 insertions(+), 12 deletions(-) diff --git a/doc/src/Howto_broken_bonds.rst b/doc/src/Howto_broken_bonds.rst index 282d2b32a4..9fcd676c36 100644 --- a/doc/src/Howto_broken_bonds.rst +++ b/doc/src/Howto_broken_bonds.rst @@ -19,6 +19,14 @@ in the broken bond are "turned on". Angles, dihedrals, etc cannot be defined for a system when :doc:`bond_style quartic ` is used. +Similarly, bond styles in the BPM package are also incompatible with +angles, dihedrals, etc. and when a bond breaks its type is set to zero. +However, in the BPM package one can either turn off all pair interactions +between bonded particles or leave them on, overlaying pair forces on +top of bond forces. To remove pair forces, the special bond list is +dynamically updated. More details can be found on the :doc:`Howto BPM +` page. + The :doc:`fix bond/break ` and :doc:`fix bond/react ` commands allow breaking of bonds within a molecular topology with may also define angles, dihedrals, etc. These commands @@ -27,12 +35,6 @@ well as the appropriate angle, dihederal, etc interactions which include the bond. They also trigger a rebuild of the neighbor list when this occurs, to turn on the appropriate pairwise forces. -In the BPM package, one can either turn off all pair interactions -between bonded particles or leave them on, overlaying pair forces on -top of bond forces. To remove pair forces, the special bond list is -dynamically updated. More details can be found on the :doc:`Howto BPM -` page. - Note that when bonds are dumped to a file via the :doc:`dump local ` command, bonds with type 0 are not included. diff --git a/doc/src/compute_count_type.rst b/doc/src/compute_count_type.rst index bb276ea21b..b92fe2f2fb 100644 --- a/doc/src/compute_count_type.rst +++ b/doc/src/compute_count_type.rst @@ -54,8 +54,9 @@ type: These commands can create and/or break topological bonds (angles, etc). In the case of breaking, they remove the bond (angle, etc) from the system, so that they no longer exist (:doc:`bond_style quartic -` is the exception, see the discussion below). Thus -they are not included in the counts for each type: +` and :doc:`BPM bond styles ` are exceptions, +see the discussion below). Thus they are not included in the counts +for each type: * :doc:`delete_bonds remove ` * :doc:`bond_style quartic ` @@ -74,10 +75,11 @@ type is tallied. Only bonds with both atoms in the specified group are counted. For {mode} = {bond}, broken bonds with a bond type of zero are also -counted. The :doc:`bond_style quartic ` breaks a bond -by doing this. See the :doc:`Howto broken bonds ` -doc page for more details. Note that the group setting is ignored for -broken bonds; all broken bonds in the system are counted. +counted. The :doc:`bond_style quartic ` and :doc:`BPM +bond styles ` break bonds by doing this. See the :doc:` +Howto broken bonds ` doc page for more details. +Note that the group setting is ignored for broken bonds; all broken +bonds in the system are counted. If the {mode} setting is {angle} then the count of angles for each angle type is tallied. Only angles with all 3 atoms in the specified From 3e3ed89f33428e6324bca28f90b3bc966f7c804a Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 2 May 2023 17:30:23 -0400 Subject: [PATCH 22/23] avoid out-of-range access if an angle/dihedral/improper type is set to 0 --- src/compute_count_type.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/compute_count_type.cpp b/src/compute_count_type.cpp index 84b31e8b85..3cb2e6bc3a 100644 --- a/src/compute_count_type.cpp +++ b/src/compute_count_type.cpp @@ -282,7 +282,7 @@ int ComputeCountType::count_angles() if ((mask[j1] & groupbit) && (mask[j2] & groupbit) && (mask[j3] & groupbit)) { if (itype > 0) count[itype - 1]++; - else + else if (itype < 0) count[-itype - 1]++; } } @@ -334,7 +334,7 @@ int ComputeCountType::count_dihedrals() (mask[j4] & groupbit)) { if (itype > 0) count[itype - 1]++; - else + else if (itype < 0) count[-itype - 1]++; } } @@ -386,7 +386,7 @@ int ComputeCountType::count_impropers() (mask[j4] & groupbit)) { if (itype > 0) count[itype - 1]++; - else + else if (itype < 0) count[-itype - 1]++; } } From 7ab30aa4684bf2bbf843f1b930274f5d0aa500de Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 2 May 2023 17:32:57 -0400 Subject: [PATCH 23/23] add versionadded tag --- doc/src/compute_count_type.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/doc/src/compute_count_type.rst b/doc/src/compute_count_type.rst index b92fe2f2fb..73e9d50f73 100644 --- a/doc/src/compute_count_type.rst +++ b/doc/src/compute_count_type.rst @@ -25,6 +25,8 @@ Examples Description """"""""""" +.. versionadded:: TBD + Define a computation that counts the current number of atoms for each atom type. Or the number of bonds (angles, dihedrals, impropers) for each bond (angle, dihedral, improper) type.