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