Merge pull request #3761 from lammps/count-types
New compute count/type command
This commit is contained in:
Axel Kohlmeyer
@ -46,6 +46,7 @@ KOKKOS, o = OPENMP, t = OPT.
|
||||
* :doc:`com/chunk <compute_com_chunk>`
|
||||
* :doc:`contact/atom <compute_contact_atom>`
|
||||
* :doc:`coord/atom (k) <compute_coord_atom>`
|
||||
* :doc:`count/type <compute_count_type>`
|
||||
* :doc:`damage/atom <compute_damage_atom>`
|
||||
* :doc:`dihedral <compute_dihedral>`
|
||||
* :doc:`dihedral/local <compute_dihedral_local>`
|
||||
|
||||
@ -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 <bond_quartic>` and the
|
||||
bond styles in the :doc:`BPM package <Howto_bpm>` which contains the
|
||||
:doc:`bpm/spring <bond_bpm_spring>` and :doc:`bpm/rotational
|
||||
<bond_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 <special_bonds>`.
|
||||
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 <pair_style>` command, whether the bond
|
||||
is broken or not. This means that :doc:`special_bonds <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 <bond_quartic>`
|
||||
* :doc:`fix bond/break <fix_bond_break>`
|
||||
* :doc:`fix bond/react <fix_bond_react>`
|
||||
* :doc:`BPM package <Howto_bpm>` 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 <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 <fix_bond_break>` and
|
||||
:doc:`fix bond/react <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 <bond_quartic>` is
|
||||
used.
|
||||
|
||||
Note that when bonds are dumped to a file via the :doc:`dump local <dump>` command, bonds with type 0 are not included. The
|
||||
:doc:`delete_bonds <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
|
||||
<Howto_bpm>` page.
|
||||
|
||||
The :doc:`fix bond/break <fix_bond_break>` and :doc:`fix bond/react
|
||||
<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
|
||||
<dump>` command, bonds with type 0 are not included.
|
||||
|
||||
The :doc:`delete_bonds <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 <compute_nbond_atom>` can also be used
|
||||
to tally the current number of bonds per atom, excluding broken bonds.
|
||||
The compute :doc:`count/type <compute_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 <compute_nbond_atom>` command tallies the
|
||||
current number of bonds each atom is part of, excluding broken bonds
|
||||
with type = 0.
|
||||
|
||||
@ -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 <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 <read_data>` or
|
||||
:doc:`read_restart <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 <read_data>` or :doc:`read_restart <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.
|
||||
|
||||
@ -200,6 +200,7 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` pag
|
||||
* :doc:`com/chunk <compute_com_chunk>` - center of mass for each chunk
|
||||
* :doc:`contact/atom <compute_contact_atom>` - contact count for each spherical particle
|
||||
* :doc:`coord/atom <compute_coord_atom>` - coordination number for each atom
|
||||
* :doc:`count/type <compute_count_type>` - count of atoms or bonds by type
|
||||
* :doc:`damage/atom <compute_damage_atom>` - Peridynamic damage for each atom
|
||||
* :doc:`dihedral <compute_dihedral>` - energy of each dihedral sub-style
|
||||
* :doc:`dihedral/local <compute_dihedral_local>` - angle of each dihedral
|
||||
|
||||
130
doc/src/compute_count_type.rst
Normal file
130
doc/src/compute_count_type.rst
Normal file
@ -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 <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 <fix_deposit>`,
|
||||
:doc:`fix pour <fix_pour>`, or :doc:`fix evaporate <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
|
||||
<read_data>` command or defined by the :doc:`molecule <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 <fix_shake>`
|
||||
* :doc:`delete_bonds <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
|
||||
<bond_quartic>` and :doc:`BPM bond styles <Howto_bpm>` are exceptions,
|
||||
see the discussion below). Thus they are not included in the counts
|
||||
for each type:
|
||||
|
||||
* :doc:`delete_bonds remove <delete_bonds>`
|
||||
* :doc:`bond_style quartic <bond_quartic>`
|
||||
* :doc:`fix bond/react <fix_bond_react>`
|
||||
* :doc:`fix bond/create <fix_bond_create>`
|
||||
* :doc:`fix bond/break <fix_bond_break>`
|
||||
* :doc:`BPM package <Howto_bpm>` 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 <bond_quartic>` and :doc:`BPM
|
||||
bond styles <Howto_bpm>` break bonds by doing this. See the :doc:`
|
||||
Howto broken bonds <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
|
||||
<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
|
||||
400
src/compute_count_type.cpp
Normal file
400
src/compute_count_type.cpp
Normal file
@ -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;
|
||||
}
|
||||
52
src/compute_count_type.h
Normal file
52
src/compute_count_type.h
Normal file
@ -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
|
||||
@ -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
|
||||
|
||||
@ -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)
|
||||
|
||||
Reference in New Issue
Block a user