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/Howto_broken_bonds.rst b/doc/src/Howto_broken_bonds.rst index 88c8acb87f..9fcd676c36 100644 --- a/doc/src/Howto_broken_bonds.rst +++ b/doc/src/Howto_broken_bonds.rst @@ -1,48 +1,56 @@ 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, molecular bond interactions persist for the duration of a +simulation in LAMMPS. However, some commands break bonds dynamically, +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 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. -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.: +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 +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. + +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 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:`nbond/atom ` can also be used -to tally the current number of bonds per atom, excluding broken bonds. +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 +with type = 0. 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..73e9d50f73 --- /dev/null +++ b/doc/src/compute_count_type.rst @@ -0,0 +1,130 @@ +.. 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} or {angle} or {dihedral} or {improper} + +Examples +"""""""" + +.. code-block:: LAMMPS + + compute 1 all count/type atom + compute 1 flowmols count/type bond + +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. + +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 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 +` 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 ` +* :doc:`fix bond/react ` +* :doc:`fix bond/create ` +* :doc:`fix bond/break ` +* :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 of zero are also +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 +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. + +---------- + +Output info +""""""""""" + +This compute calculates a global vector of counts. If the mode is +{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 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 +` 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 new file mode 100644 index 0000000000..3cb2e6bc3a --- /dev/null +++ b/src/compute_count_type.cpp @@ -0,0 +1,400 @@ +/* ---------------------------------------------------------------------- + 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 "domain.h" +#include "error.h" +#include "force.h" +#include "group.h" +#include "update.h" + +using namespace LAMMPS_NS; + +enum { ATOM, BOND, ANGLE, DIHEDRAL, IMPROPER }; + +/* ---------------------------------------------------------------------- */ + +ComputeCountType::ComputeCountType(LAMMPS *lmp, int narg, char **arg) : + 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"); + + // 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]); + + // 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; + extvector = 1; + } else if (mode == BOND) { + scalar_flag = vector_flag = 1; + 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; + } + + // output vector + + 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; + + delete[] count; + delete[] bcount_me; + delete[] bcount; +} + +/* ---------------------------------------------------------------------- + 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; + + // 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++) { + 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; + 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"); + scalar = bcount; + return scalar; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeCountType::compute_vector() +{ + invoked_vector = update->ntimestep; + + 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(); + + // 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; + 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() +{ + tagint **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 flag = 0; + for (int m = 0; m < nbondtypes; m++) count[m] = 0; + + for (int i = 0; i < nlocal; i++) { + int nbond = num_bond[i]; + for (int m = 0; m < nbond; m++) { + int itype = bond_type[i][m]; + if (itype == 0) continue; + + int 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() +{ + 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; + int nlocal = atom->nlocal; + int nangletypes = atom->nangletypes; + + int flag = 0; + for (int m = 0; m < nangletypes; m++) count[m] = 0; + + for (int i = 0; i < nlocal; i++) { + int nangle = num_angle[i]; + for (int m = 0; m < nangle; m++) { + int itype = angle_type[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; + } + + if ((mask[j1] & groupbit) && (mask[j2] & groupbit) && (mask[j3] & groupbit)) { + if (itype > 0) + count[itype - 1]++; + else if (itype < 0) + 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() +{ + 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; + int nlocal = atom->nlocal; + int ndihedraltypes = atom->ndihedraltypes; + + int flag = 0; + for (int m = 0; m < ndihedraltypes; m++) count[m] = 0; + + for (int i = 0; i < nlocal; i++) { + int ndihedral = num_dihedral[i]; + for (int m = 0; m < ndihedral; m++) { + int itype = dihedral_type[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; + } + + if ((mask[j1] & groupbit) && (mask[j2] & groupbit) && (mask[j3] & groupbit) && + (mask[j4] & groupbit)) { + if (itype > 0) + count[itype - 1]++; + else if (itype < 0) + 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() +{ + 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; + int nlocal = atom->nlocal; + int nimpropertypes = atom->nimpropertypes; + + int flag = 0; + for (int m = 0; m < nimpropertypes; m++) count[m] = 0; + + for (int i = 0; i < nlocal; i++) { + int nimproper = num_improper[i]; + for (int m = 0; m < nimproper; m++) { + int itype = improper_type[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; + } + + if ((mask[j1] & groupbit) && (mask[j2] & groupbit) && (mask[j3] & groupbit) && + (mask[j4] & groupbit)) { + if (itype > 0) + count[itype - 1]++; + else if (itype < 0) + 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 new file mode 100644 index 0000000000..b766edfae8 --- /dev/null +++ b/src/compute_count_type.h @@ -0,0 +1,52 @@ +/* -*- 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; + + int count_atoms(); + int count_bonds(); + int count_angles(); + int count_dihedrals(); + int count_impropers(); +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/lmptype.h b/src/lmptype.h index bf56d07d77..ecd6ee5761 100644 --- a/src/lmptype.h +++ b/src/lmptype.h @@ -95,6 +95,7 @@ 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 @@ -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 diff --git a/unittest/commands/test_compute_global.cpp b/unittest/commands/test_compute_global.cpp index d7573d8d16..6c365c8c2b 100644 --- a/unittest/commands/test_compute_global.cpp +++ b/unittest/commands/test_compute_global.cpp @@ -294,6 +294,98 @@ 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 zero 10.0"); + command("pair_coeff * *"); + + 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 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("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 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(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 remove"); + command("run 0 post no"); + END_HIDE_OUTPUT(); + + bcnt = get_vector("bcnt"); + bbrk = get_scalar("bcnt"); + 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], 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 int main(int argc, char **argv)