tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'fortran-fix-external' of github.com:hammondkd/lammps into fortran-fix-external

Browse Source
This commit is contained in:
Karl Hammond
2022-12-02 00:02:10 -06:00
parent b1664ce8ea 4c9ca8761c
commit 43dca96ca4
58 changed files with 1795 additions and 745 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
.github/CODEOWNERS vendored
View File

@ -37,6 +37,7 @@ src/MESONT/* @iafoss
src/ML-HDNNP/* @singraber
src/ML-IAP/* @athomps
src/ML-PACE/* @yury-lysogorskiy
src/ML-POD/* @exapde @rohskopf
src/MOFFF/* @hheenen
src/MOLFILE/* @akohlmey
src/NETCDF/* @pastewka

3
doc/src/Commands_all.rst
View File

@ -89,8 +89,7 @@ table above.
* :doc:`region <region>`
* :doc:`replicate <replicate>`
* :doc:`rerun <rerun>`
* :doc:`reset_atom_ids <reset_atom_ids>`
* :doc:`reset_mol_ids <reset_mol_ids>`
* :doc:`reset_atoms <reset_atoms>`
* :doc:`reset_timestep <reset_timestep>`
* :doc:`restart <restart>`
* :doc:`run <run>`

33
doc/src/Commands_removed.rst
View File

@ -26,10 +26,15 @@ Box command
The *box* command has been removed and the LAMMPS code changed so it won't
be needed. If present, LAMMPS will ignore the command and print a warning.
Reset_ids command
-----------------
Reset_ids, reset_atom_ids, reset_mol_ids commands
-------------------------------------------------
The reset_ids command has been renamed to :doc:`reset_atom_ids <reset_atom_ids>`.
.. deprecated:: TBD
The *reset_ids*, *reset_atom_ids*, and *reset_mol_ids* commands have
been folded into the :doc:`reset_atoms <reset_atoms>` command. If
present, LAMMPS will replace the commands accordingly and print a
warning.
MEAM package
------------
@ -39,18 +44,21 @@ The code in the :ref:`MEAM package <PKG-MEAM>` is a translation of the
Fortran code of MEAM into C++, which removes several restrictions
(e.g. there can be multiple instances in hybrid pair styles) and allows
for some optimizations leading to better performance. The pair style
:doc:`meam <pair_meam>` has the exact same syntax.
:doc:`meam <pair_meam>` has the exact same syntax. For a transition
period the C++ version of MEAM was called USER-MEAMC so it could
coexist with the Fortran version.
REAX package
------------
The REAX package has been removed since it was superseded by the
:ref:`REAXFF package <PKG-REAXFF>`. The REAXFF
package has been tested to yield equivalent results to the REAX package,
offers better performance, supports OpenMP multi-threading via OPENMP,
and GPU and threading parallelization through KOKKOS. The new pair styles
are not syntax compatible with the removed reax pair style, so input
files will have to be adapted.
:ref:`REAXFF package <PKG-REAXFF>`. The REAXFF package has been tested
to yield equivalent results to the REAX package, offers better
performance, supports OpenMP multi-threading via OPENMP, and GPU and
threading parallelization through KOKKOS. The new pair styles are not
syntax compatible with the removed reax pair style, so input files will
have to be adapted. The REAXFF package was originally called
USER-REAXC.
USER-CUDA package
-----------------
@ -69,5 +77,6 @@ restart2data tool
The functionality of the restart2data tool has been folded into the
LAMMPS executable directly instead of having a separate tool. A
combination of the commands :doc:`read_restart <read_restart>` and
:doc:`write_data <write_data>` can be used to the same effect. For added
convenience this conversion can also be triggered by :doc:`command line flags <Run_options>`
:doc:`write_data <write_data>` can be used to the same effect. For
added convenience this conversion can also be triggered by
:doc:`command line flags <Run_options>`

8
doc/src/Developer_code_design.rst
View File

@ -78,7 +78,7 @@ LAMMPS makes extensive use of the object oriented programming (OOP)
principles of *compositing* and *inheritance*. Classes like the
``LAMMPS`` class are a **composite** containing pointers to instances
of other classes like ``Atom``, ``Comm``, ``Force``, ``Neighbor``,
``Modify``, and so on. Each of these classes implement certain
``Modify``, and so on. Each of these classes implements certain
functionality by storing and manipulating data related to the
simulation and providing member functions that trigger certain
actions. Some of those classes like ``Force`` are themselves
@ -87,9 +87,9 @@ interactions. Similarly the ``Modify`` class contains a list of
``Fix`` and ``Compute`` classes. If the input commands that
correspond to these classes include the word *style*, then LAMMPS
stores only a single instance of that class. E.g. *atom_style*,
*comm_style*, *pair_style*, *bond_style*. It the input command does
not include the word *style*, there can be many instances of that
class defined. E.g. *region*, *fix*, *compute*, *dump*.
*comm_style*, *pair_style*, *bond_style*. If the input command does
**not** include the word *style*, then there may be many instances of
that class defined, for example *region*, *fix*, *compute*, *dump*.
**Inheritance** enables creation of *derived* classes that can share
common functionality in their base class while providing a consistent

6
doc/src/Developer_unittest.rst
View File

@ -226,9 +226,9 @@ The following test programs are currently available:
* - ``test_kim_commands.cpp``
- KimCommands
- Tests for several commands from the :ref:`KIM package <PKG-KIM>`
* - ``test_reset_ids.cpp``
- ResetIDs
- Tests to validate the :doc:`reset_atom_ids <reset_atom_ids>` and :doc:`reset_mol_ids <reset_mol_ids>` commands
* - ``test_reset_atoms.cpp``
- ResetAtoms
- Tests to validate the :doc:`reset_atoms <reset_atoms>` sub-commands
Tests for the C-style library interface

3
doc/src/commands_list.rst
View File

@ -89,8 +89,7 @@ Commands
region
replicate
rerun
reset_atom_ids
reset_mol_ids
reset_atoms
reset_timestep
restart
run

2
doc/src/compute_pair_local.rst
View File

@ -59,7 +59,7 @@ commands.
The value *dist* is the distance between the pair of atoms.
The values *dx*, *dy*, and *dz* are the :math:`(x,y,z)` components of the
*distance* between the pair of atoms. This value is always the
distance from the atom of lower to the one with the higher id.
distance from the atom of higher to the one with the lower atom ID.
The value *eng* is the interaction energy for the pair of atoms.

4
doc/src/delete_atoms.rst
View File

@ -135,7 +135,7 @@ number of atoms in the system. Note that this is not done for
molecular systems (see the :doc:`atom_style <atom_style>` command),
regardless of the *compress* setting, since it would foul up the bond
connectivity that has already been assigned. However, the
:doc:`reset_atom_ids <reset_atom_ids>` command can be used after this
:doc:`reset_atoms id <reset_atoms>` command can be used after this
command to accomplish the same thing.
Note that the re-assignment of IDs is not really a compression, where
@ -203,7 +203,7 @@ using molecule template files via the :doc:`molecule <molecule>` and
Related commands
""""""""""""""""
:doc:`create_atoms <create_atoms>`, :doc:`reset_atom_ids <reset_atom_ids>`
:doc:`create_atoms <create_atoms>`, :doc:`reset_atoms id <reset_atoms>`
Default
"""""""

12
doc/src/fix_bond_react.rst
View File

@ -177,12 +177,12 @@ due to the internal dynamic grouping performed by fix bond/react.
If the group-ID is an existing static group, react-group-IDs
should also be specified as this static group or a subset.
The *reset_mol_ids* keyword invokes the :doc:`reset_mol_ids <reset_mol_ids>`
command after a reaction occurs, to ensure that molecule IDs are
consistent with the new bond topology. The group-ID used for
:doc:`reset_mol_ids <reset_mol_ids>` is the group-ID for this fix.
Resetting molecule IDs is necessarily a global operation, so it can
be slow for very large systems.
The *reset_mol_ids* keyword invokes the :doc:`reset_atoms mol
<reset_atoms>` command after a reaction occurs, to ensure that
molecule IDs are consistent with the new bond topology. The group-ID
used for :doc:`reset_atoms mol <reset_atoms>` is the group-ID for this
fix. Resetting molecule IDs is necessarily a global operation, so it
can be slow for very large systems.
The following comments pertain to each *react* argument (in other
words, they can be customized for each reaction, or reaction step):

6
doc/src/fix_ipi.rst
View File

@ -90,6 +90,12 @@ coordinates are transferred. However, one could use this strategy to
define an external potential acting on the atoms that are moved by
i-PI.
Since the i-PI code uses atomic units internally, this fix needs to
convert LAMMPS data to and from its :doc:`specified units <units>`
accordingly when communicating with i-PI. This is not possible for
reduced units ("units lj") and thus *fix ipi* will stop with an error in
this case.
This fix is part of the MISC package. It is only enabled if
LAMMPS was built with that package. See the
:doc:`Build package <Build_package>` page for more info.

2
doc/src/fix_pimd.rst
View File

@ -156,6 +156,8 @@ This fix is part of the REPLICA package. It is only enabled if
LAMMPS was built with that package. See the
:doc:`Build package <Build_package>` page for more info.
Fix pid cannot be used with :doc:`lj units <units>`.
A PIMD simulation can be initialized with a single data file read via
the :doc:`read_data <read_data>` command. However, this means all
quasi-beads in a ring polymer will have identical positions and

10
doc/src/pair_sw.rst
View File

@ -176,9 +176,13 @@ are placeholders for atom types that will be used with other potentials.
.. note::
When the *threebody off* keyword is used, multiple pair_coeff commands may
be used to specific the pairs of atoms which don't require three-body term.
In these cases, the first 2 arguments are not required to be \* \*.
When the *threebody off* keyword is used, multiple pair_coeff
commands may be used to specific the pairs of atoms which don't
require three-body term. In these cases, the first 2 arguments are
not required to be \* \*, the potential parameter file is only read
by the first :doc:`pair_coeff command <pair_coeff>` and the element
to atom type mappings must be consistent across all *pair_coeff*
statements. If not LAMMPS will abort with an error.
Stillinger-Weber files in the *potentials* directory of the LAMMPS
distribution have a ".sw" suffix. Lines that are not blank or

94
doc/src/reset_atom_ids.rst
View File

@ -1,94 +0,0 @@
.. index:: reset_atom_ids
reset_atom_ids command
======================
Syntax
""""""
.. code-block:: LAMMPS
reset_atom_ids keyword values ...
* zero or more keyword/value pairs may be appended
* keyword = *sort*
.. parsed-literal::
*sort* value = *yes* or *no*
Examples
""""""""
.. code-block:: LAMMPS
reset_atom_ids
reset_atom_ids sort yes
Description
"""""""""""
Reset atom IDs for the system, including all the global IDs stored
for bond, angle, dihedral, improper topology data. This will
create a set of IDs that are numbered contiguously from 1 to N
for a N atoms system.
This can be useful to do after performing a "delete_atoms" command for
a molecular system. The delete_atoms compress yes option will not
perform this operation due to the existence of bond topology. It can
also be useful to do after any simulation which has lost atoms,
e.g. due to atoms moving outside a simulation box with fixed
boundaries (see the "boundary command"), or due to evaporation (see
the "fix evaporate" command).
If the *sort* keyword is used with a setting of *yes*, then the
assignment of new atom IDs will be the same no matter how many
processors LAMMPS is running on. This is done by first doing a
spatial sort of all the atoms into bins and sorting them within each
bin. Because the set of bins is independent of the number of
processors, this enables a consistent assignment of new IDs to each
atom.
This can be useful to do after using the "create_atoms" command and/or
"replicate" command. In general those commands do not guarantee
assignment of the same atom ID to the same physical atom when LAMMPS
is run on different numbers of processors. Enforcing consistent IDs
can be useful for debugging or comparing output from two different
runs.
Note that the spatial sort requires communication of atom IDs and
coordinates between processors in an all-to-all manner. This is done
efficiently in LAMMPS, but it is more expensive than how atom IDs are
reset without sorting.
Note that whether sorting or not, the resetting of IDs is not a
compression, where gaps in atom IDs are removed by decrementing atom
IDs that are larger. Instead the IDs for all atoms are erased, and
new IDs are assigned so that the atoms owned by an individual
processor have consecutive IDs, as the :doc:`create_atoms
<create_atoms>` command explains.
.. note::
If this command is used before a :doc:`pair style <pair_style>` is
defined, an error about bond topology atom IDs not being found may
result. This is because the cutoff distance for ghost atom
communication was not sufficient to find atoms in bonds, angles, etc
that are owned by other processors. The :doc:`comm_modify cutoff <comm_modify>` command can be used to correct this issue.
Or you can define a pair style before using this command. If you do
the former, you should unset the comm_modify cutoff after using
reset_atom_ids so that subsequent communication is not inefficient.
Restrictions
""""""""""""
none
Related commands
""""""""""""""""
:doc:`delete_atoms <delete_atoms>`
Default
"""""""
By default, *sort* is no.

283
doc/src/reset_atoms.rst Normal file
View File

@ -0,0 +1,283 @@
.. index:: reset_atoms
reset_atoms command
===================
Syntax
""""""
.. code-block:: LAMMPS
reset_atoms property arguments ...
* property = *id* or *image* or *mol*
* additional arguments depend on the property
.. code-block:: LAMMPS
reset_atoms id keyword value ...
* zero or more keyword/value pairs can be appended
* keyword = *sort*
.. parsed-literal::
*sort* value = *yes* or *no*
.. code-block:: LAMMPS
reset_atoms image group-ID
* group-ID = ID of group of atoms whose image flags will be reset
.. code-block:: LAMMPS
reset atoms mol group-ID keyword value ...
* group-ID = ID of group of atoms whose molecule IDs will be reset
* zero or more keyword/value pairs can be appended
* keyword = *compress* or *offset* or *single*
.. parsed-literal::
*compress* value = *yes* or *no*
*offset* value = *Noffset* >= -1
*single* value = *yes* or *no* to treat single atoms (no bonds) as molecules
Examples
""""""""
.. code-block:: LAMMPS
reset_atoms id
reset_atoms id sort yes
reset_atoms image all
reset_atoms image mobile
reset_atoms mol all
reset_atoms mol all offset 10 single yes
reset_atoms mol solvent compress yes offset 100
reset_atoms mol solvent compress no
Description
"""""""""""
.. versionadded:: TBD
The *reset_atoms* command resets the values of a specified atom
property. In contrast to the set command, it does this in a
collective manner which resets the values for many atoms in a
self-consistent way. This is often useful when the simulated system
has undergone significant modifications like adding or removing atoms
or molecules, joining data files, changing bonds, or large-scale
diffusion.
The new values can be thought of as a *reset*, similar to values atoms
would have if a new data file were being read or a new simulation
performed. Note that the set command also resets atom properties to
new values, but it treats each atom independently.
The *property* setting can be *id* or *image* or *mol*. For *id*, the
IDs of all the atoms are reset to contiguous values. For *image*, the
image flags of atoms in the specified *group-ID* are reset so that at
least one atom in each molecule is in the simulation box (image flag =
0). For *mol*, the molecule IDs of all atoms are reset to contiguous
values.
More details on these operations and their arguments or optional
keyword/value settings are given below.
----------
*Property id*
Reset atom IDs for the entire system, including all the global IDs
stored for bond, angle, dihedral, improper topology data. This will
create a set of IDs that are numbered contiguously from 1 to N for a N
atoms system.
This can be useful to do after performing a "delete_atoms" command for
a molecular system. The delete_atoms compress yes option will not
perform this operation due to the existence of bond topology. It can
also be useful to do after any simulation which has lost atoms,
e.g. due to atoms moving outside a simulation box with fixed
boundaries (see the "boundary command"), or due to evaporation (see
the "fix evaporate" command).
If the *sort* keyword is used with a setting of *yes*, then the
assignment of new atom IDs will be the same no matter how many
processors LAMMPS is running on. This is done by first doing a
spatial sort of all the atoms into bins and sorting them within each
bin. Because the set of bins is independent of the number of
processors, this enables a consistent assignment of new IDs to each
atom.
This can be useful to do after using the "create_atoms" command and/or
"replicate" command. In general those commands do not guarantee
assignment of the same atom ID to the same physical atom when LAMMPS
is run on different numbers of processors. Enforcing consistent IDs
can be useful for debugging or comparing output from two different
runs.
Note that the spatial sort requires communication of atom IDs and
coordinates between processors in an all-to-all manner. This is done
efficiently in LAMMPS, but it is more expensive than how atom IDs are
reset without sorting.
Note that whether sorting or not, the resetting of IDs is not a
compression, where gaps in atom IDs are removed by decrementing atom
IDs that are larger. Instead the IDs for all atoms are erased, and
new IDs are assigned so that the atoms owned by an individual
processor have consecutive IDs, as the :doc:`create_atoms
<create_atoms>` command explains.
.. note::
If this command is used before a :doc:`pair style <pair_style>` is
defined, an error about bond topology atom IDs not being found may
result. This is because the cutoff distance for ghost atom
communication was not sufficient to find atoms in bonds, angles, etc
that are owned by other processors. The :doc:`comm_modify cutoff
<comm_modify>` command can be used to correct this issue. Or you can
define a pair style before using this command. If you do the former,
you should unset the *comm_modify cutoff* after using *reset
atoms id* so that subsequent communication is not inefficient.
----------
*Property image*
Reset the image flags of atoms so that at least one atom in each
molecule has an image flag of 0. Molecular topology is respected so
that if the molecule straddles a periodic simulation box boundary, the
images flags of all atoms in the molecule will be consistent. This
avoids inconsistent image flags that could result from resetting all
image flags to zero with the :doc:`set <set>` command.
.. note::
If the system has no bonds, there is no reason to use this command,
since image flags for different atoms do not need to be
consistent. Use the :doc:`set <set>` command with its *image*
keyword instead.
Only image flags for atoms in the specified *group-ID* are reset; all
others remain unchanged. No check is made for whether the group
covers complete molecule fragments and thus whether the command will
result in inconsistent image flags.
Molecular fragments are identified by the algorithm used by the
:doc:`compute fragment/atom <compute_cluster_atom>` command. For each
fragment the average of the largest and the smallest image flag in
each direction across all atoms in the fragment is computed and
subtracted from the current image flag in the same direction.
This can be a useful operation to perform after running longer
equilibration runs of mobile systems where molecules would pass
through the system multiple times and thus produce non-zero image
flags.
.. note::
Same as explained for the :doc:`compute fragment/atom
<compute_cluster_atom>` command, molecules are identified using the
current bond topology. This will **not** account for bonds broken by
the :doc:`bond_style quartic <bond_quartic>` command, because this
bond style does not perform a full update of the bond topology data
structures within LAMMPS. In that case, using the :doc:`delete_bonds
all bond 0 remove <delete_bonds>` will permanently delete such
broken bonds and should thus be used first.
----------
*Property mol*
Reset molecule IDs for a specified group of atoms based on current
bond connectivity. This will typically create a new set of molecule
IDs for atoms in the group. Only molecule IDs for atoms in the
specified *group-ID* are reset; molecule IDs for atoms not in the
group are not changed.
For purposes of this operation, molecules are identified by the current
bond connectivity in the system, which may or may not be consistent with
the current molecule IDs. A molecule in this context is a set of atoms
connected to each other with explicit bonds. The specific algorithm
used is the one of :doc:`compute fragment/atom <compute_cluster_atom>`
Once the molecules are identified and a new molecule ID computed for
each, this command will update the current molecule ID for all atoms in
the group with the new molecule ID. Note that if the group excludes
atoms within molecules, one (physical) molecule may become two or more
(logical) molecules. For example if the group excludes atoms in the
middle of a linear chain, then each end of the chain is considered an
independent molecule and will be assigned a different molecule ID.
This can be a useful operation to perform after running reactive
molecular dynamics run with :doc:`fix bond/react <fix_bond_react>`,
:doc:`fix bond/create <fix_bond_create>`, or :doc:`fix bond/break
<fix_bond_break>`, all of which can change molecule topologies. It can
also be useful after molecules have been deleted with the
:doc:`delete_atoms <delete_atoms>` command or after a simulation which
has lost molecules, e.g. via the :doc:`fix evaporate <fix_evaporate>`
command.
The *compress* keyword determines how new molecule IDs are computed. If
the setting is *yes* (the default) and there are N molecules in the
group, the new molecule IDs will be a set of N contiguous values. See
the *offset* keyword for details on selecting the range of these values.
If the setting is *no*, the molecule ID of every atom in the molecule
will be set to the smallest atom ID of any atom in the molecule.
The *single* keyword determines whether single atoms (not bonded to
another atom) are treated as one-atom molecules or not, based on the
*yes* or *no* setting. If the setting is *no* (the default), their
molecule IDs are set to 0. This setting can be important if the new
molecule IDs will be used as input to other commands such as
:doc:`compute chunk/atom molecule <compute_chunk_atom>` or :doc:`fix
rigid molecule <fix_rigid>`.
The *offset* keyword is only used if the *compress* setting is *yes*.
Its default value is *Noffset* = -1. In that case, if the specified
group is *all*, then the new compressed molecule IDs will range from 1
to N. If the specified group is not *all* and the largest molecule ID
of atoms outside that group is M, then the new compressed molecule IDs will
range from M+1 to M+N, to avoid collision with existing molecule
IDs. If an *Noffset* >= 0 is specified, then the new compressed
molecule IDs will range from *Noffset*\ +1 to *Noffset*\ +N. If the group
is not *all* there may be collisions with the molecule IDs of other atoms.
.. note::
Same as explained for the :doc:`compute fragment/atom
<compute_cluster_atom>` command, molecules are identified using the
current bond topology. This will **not** account for bonds broken by
the :doc:`bond_style quartic <bond_quartic>` command, because this
bond style does not perform a full update of the bond topology data
structures within LAMMPS. In that case, using the :doc:`delete_bonds
all bond 0 remove <delete_bonds>` will permanently delete such broken
bonds and should thus be used first.
Restrictions
""""""""""""
The *image* property can only be used when the atom style supports bonds.
Related commands
""""""""""""""""
:doc:`compute fragment/atom <compute_cluster_atom>`
:doc:`fix bond/react <fix_bond_react>`,
:doc:`fix bond/create <fix_bond_create>`,
:doc:`fix bond/break <fix_bond_break>`,
:doc:`fix evaporate <fix_evaporate>`,
:doc:`delete_atoms <delete_atoms>`,
:doc:`delete_bonds <delete_bonds>`
Defaults
""""""""
For property *id*, the default keyword setting is sort = no.
For property *mol*, the default keyword settings are compress = yes,
single = no, and offset = -1.

116
doc/src/reset_mol_ids.rst
View File

@ -1,116 +0,0 @@
.. index:: reset_mol_ids
reset_mol_ids command
=====================
Syntax
""""""
.. parsed-literal::
reset_mol_ids group-ID keyword value ...
* group-ID = ID of group of atoms whose molecule IDs will be reset
* zero or more keyword/value pairs may be appended
* keyword = *compress* or *offset* or *single*
.. parsed-literal::
*compress* value = *yes* or *no*
*offset* value = *Noffset* >= -1
*single* value = *yes* or *no* to treat single atoms (no bonds) as molecules
Examples
""""""""
.. code-block:: LAMMPS
reset_mol_ids all
reset_mol_ids all offset 10 single yes
reset_mol_ids solvent compress yes offset 100
reset_mol_ids solvent compress no
Description
"""""""""""
Reset molecule IDs for a group of atoms based on current bond
connectivity. This will typically create a new set of molecule IDs
for atoms in the group. Only molecule IDs for atoms in the specified
group are reset; molecule IDs for atoms not in the group are not
changed.
For purposes of this operation, molecules are identified by the current
bond connectivity in the system, which may or may not be consistent with
the current molecule IDs. A molecule in this context is a set of atoms
connected to each other with explicit bonds. The specific algorithm
used is the one of :doc:`compute fragment/atom <compute_cluster_atom>`
Once the molecules are identified and a new molecule ID computed for
each, this command will update the current molecule ID for all atoms in
the group with the new molecule ID. Note that if the group excludes
atoms within molecules, one (physical) molecule may become two or more
(logical) molecules. For example if the group excludes atoms in the
middle of a linear chain, then each end of the chain is considered an
independent molecule and will be assigned a different molecule ID.
This can be a useful operation to perform after running reactive
molecular dynamics run with :doc:`fix bond/react <fix_bond_react>`,
:doc:`fix bond/create <fix_bond_create>`, or :doc:`fix bond/break
<fix_bond_break>`, all of which can change molecule topologies. It can
also be useful after molecules have been deleted with the
:doc:`delete_atoms <delete_atoms>` command or after a simulation which
has lost molecules, e.g. via the :doc:`fix evaporate <fix_evaporate>`
command.
The *compress* keyword determines how new molecule IDs are computed. If
the setting is *yes* (the default) and there are N molecules in the
group, the new molecule IDs will be a set of N contiguous values. See
the *offset* keyword for details on selecting the range of these values.
If the setting is *no*, the molecule ID of every atom in the molecule
will be set to the smallest atom ID of any atom in the molecule.
The *single* keyword determines whether single atoms (not bonded to
another atom) are treated as one-atom molecules or not, based on the
*yes* or *no* setting. If the setting is *no* (the default), their
molecule IDs are set to 0. This setting can be important if the new
molecule IDs will be used as input to other commands such as
:doc:`compute chunk/atom molecule <compute_chunk_atom>` or :doc:`fix
rigid molecule <fix_rigid>`.
The *offset* keyword is only used if the *compress* setting is *yes*.
Its default value is *Noffset* = -1. In that case, if the specified
group is *all*, then the new compressed molecule IDs will range from 1
to N. If the specified group is not *all* and the largest molecule ID
of atoms outside that group is M, then the new compressed molecule IDs will
range from M+1 to M+N, to avoid collision with existing molecule
IDs. If an *Noffset* >= 0 is specified, then the new compressed
molecule IDs will range from *Noffset*\ +1 to *Noffset*\ +N. If the group
is not *all* there may be collisions with the molecule IDs of other atoms.
.. note::
The same as explained for the :doc:`compute fragment/atom
<compute_cluster_atom>` command, molecules are identified using the
current bond topology. This will not account for bonds broken by
the :doc:`bond_style quartic <bond_quartic>` command because it
does not perform a full update of the bond topology data structures
within LAMMPS.
Restrictions
""""""""""""
none
Related commands
""""""""""""""""
:doc:`reset_atom_ids <reset_atom_ids>`, :doc:`fix bond/react <fix_bond_react>`,
:doc:`fix bond/create <fix_bond_create>`,
:doc:`fix bond/break <fix_bond_break>`,
:doc:`fix evaporate <fix_evaporate>`,
:doc:`delete_atoms <delete_atoms>`,
:doc:`compute fragment/atom <compute_cluster_atom>`
Default
"""""""
The default keyword settings are compress = yes, single = no, and
offset = -1.

16
lib/linalg/README
View File

@ -1,13 +1,16 @@
This directory has BLAS and LAPACK files needed by the ATC and
AWPMD packages, and possibly by other packages in the future.
This directory has generic BLAS and LAPACK source files needed by the
ATC, AWPMD, ELECTRODE, LATTE, and ML-POD packages (and possibly by other
packages) in the future that can be used instead of platform or vendor
optimized BLAS/LAPACK library.
Note that this is an *incomplete* subset of full BLAS/LAPACK.
The files correspond to LAPACK version 3.11.0.
You should only need to build and use the library in this directory if
you want to build LAMMPS with the ATC and/or AWPMD packages
AND you do not have any other suitable BLAS and LAPACK libraries
installed on your system. E.g. ATLAS, GOTO-BLAS, OpenBLAS, ACML, or
MKL.
you want to build LAMMPS with one of the listed packages AND you do not
have any other suitable BLAS and LAPACK libraries installed on your
system (like ATLAS, GOTO-BLAS, OpenBLAS, ACML, MKL and so on).
You can type "make lib-linalg" from the src directory to see help on
how to build this library via make commands, or you can do the same
@ -27,3 +30,4 @@ liblinalg.a the library LAMMPS will link against
You can then include this library and its path in the Makefile.lammps
file of any packages that need it. As an example, see the
lib/atc/Makefile.lammps.linalg file.

26
lib/linalg/dgelss.f
View File

@ -254,11 +254,11 @@
*
* Compute space needed for DGEQRF
CALL DGEQRF( M, N, A, LDA, DUM(1), DUM(1), -1, INFO )
LWORK_DGEQRF=DUM(1)
LWORK_DGEQRF = INT( DUM(1) )
* Compute space needed for DORMQR
CALL DORMQR( 'L', 'T', M, NRHS, N, A, LDA, DUM(1), B,
$ LDB, DUM(1), -1, INFO )
LWORK_DORMQR=DUM(1)
LWORK_DORMQR = INT( DUM(1) )
MM = N
MAXWRK = MAX( MAXWRK, N + LWORK_DGEQRF )
MAXWRK = MAX( MAXWRK, N + LWORK_DORMQR )
@ -273,15 +273,15 @@
* Compute space needed for DGEBRD
CALL DGEBRD( MM, N, A, LDA, S, DUM(1), DUM(1),
$ DUM(1), DUM(1), -1, INFO )
LWORK_DGEBRD=DUM(1)
LWORK_DGEBRD = INT( DUM(1) )
* Compute space needed for DORMBR
CALL DORMBR( 'Q', 'L', 'T', MM, NRHS, N, A, LDA, DUM(1),
$ B, LDB, DUM(1), -1, INFO )
LWORK_DORMBR=DUM(1)
LWORK_DORMBR = INT( DUM(1) )
* Compute space needed for DORGBR
CALL DORGBR( 'P', N, N, N, A, LDA, DUM(1),
$ DUM(1), -1, INFO )
LWORK_DORGBR=DUM(1)
LWORK_DORGBR = INT( DUM(1) )
* Compute total workspace needed
MAXWRK = MAX( MAXWRK, 3*N + LWORK_DGEBRD )
MAXWRK = MAX( MAXWRK, 3*N + LWORK_DORMBR )
@ -305,23 +305,23 @@
* Compute space needed for DGELQF
CALL DGELQF( M, N, A, LDA, DUM(1), DUM(1),
$ -1, INFO )
LWORK_DGELQF=DUM(1)
LWORK_DGELQF = INT( DUM(1) )
* Compute space needed for DGEBRD
CALL DGEBRD( M, M, A, LDA, S, DUM(1), DUM(1),
$ DUM(1), DUM(1), -1, INFO )
LWORK_DGEBRD=DUM(1)
LWORK_DGEBRD = INT( DUM(1) )
* Compute space needed for DORMBR
CALL DORMBR( 'Q', 'L', 'T', M, NRHS, N, A, LDA,
$ DUM(1), B, LDB, DUM(1), -1, INFO )
LWORK_DORMBR=DUM(1)
LWORK_DORMBR = INT( DUM(1) )
* Compute space needed for DORGBR
CALL DORGBR( 'P', M, M, M, A, LDA, DUM(1),
$ DUM(1), -1, INFO )
LWORK_DORGBR=DUM(1)
LWORK_DORGBR = INT( DUM(1) )
* Compute space needed for DORMLQ
CALL DORMLQ( 'L', 'T', N, NRHS, M, A, LDA, DUM(1),
$ B, LDB, DUM(1), -1, INFO )
LWORK_DORMLQ=DUM(1)
LWORK_DORMLQ = INT( DUM(1) )
* Compute total workspace needed
MAXWRK = M + LWORK_DGELQF
MAXWRK = MAX( MAXWRK, M*M + 4*M + LWORK_DGEBRD )
@ -341,15 +341,15 @@
* Compute space needed for DGEBRD
CALL DGEBRD( M, N, A, LDA, S, DUM(1), DUM(1),
$ DUM(1), DUM(1), -1, INFO )
LWORK_DGEBRD=DUM(1)
LWORK_DGEBRD = INT( DUM(1) )
* Compute space needed for DORMBR
CALL DORMBR( 'Q', 'L', 'T', M, NRHS, M, A, LDA,
$ DUM(1), B, LDB, DUM(1), -1, INFO )
LWORK_DORMBR=DUM(1)
LWORK_DORMBR = INT( DUM(1) )
* Compute space needed for DORGBR
CALL DORGBR( 'P', M, N, M, A, LDA, DUM(1),
$ DUM(1), -1, INFO )
LWORK_DORGBR=DUM(1)
LWORK_DORGBR = INT( DUM(1) )
MAXWRK = 3*M + LWORK_DGEBRD
MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORMBR )
MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORGBR )

7
lib/linalg/dlaed4.f
View File

@ -328,9 +328,12 @@
IF( C.LT.ZERO )
$ C = ABS( C )
IF( C.EQ.ZERO ) THEN
* ETA = B/A
* ETA = B/A
* ETA = RHO - TAU
ETA = DLTUB - TAU
* ETA = DLTUB - TAU
*
* Update proposed by Li, Ren-Cang:
ETA = -W / ( DPSI+DPHI )
ELSE IF( A.GE.ZERO ) THEN
ETA = ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
ELSE

2
lib/linalg/dlascl.f
View File

@ -272,6 +272,8 @@
ELSE
MUL = CTOC / CFROMC
DONE = .TRUE.
IF (MUL .EQ. ONE)
$ RETURN
END IF
END IF
*

69
lib/linalg/dlatrs.f
View File

@ -264,8 +264,8 @@
* .. External Functions ..
LOGICAL LSAME
INTEGER IDAMAX
DOUBLE PRECISION DASUM, DDOT, DLAMCH
EXTERNAL LSAME, IDAMAX, DASUM, DDOT, DLAMCH
DOUBLE PRECISION DASUM, DDOT, DLAMCH, DLANGE
EXTERNAL LSAME, IDAMAX, DASUM, DDOT, DLAMCH, DLANGE
* ..
* .. External Subroutines ..
EXTERNAL DAXPY, DSCAL, DTRSV, XERBLA
@ -304,6 +304,7 @@
*
* Quick return if possible
*
SCALE = ONE
IF( N.EQ.0 )
$ RETURN
*
@ -311,7 +312,6 @@
*
SMLNUM = DLAMCH( 'Safe minimum' ) / DLAMCH( 'Precision' )
BIGNUM = ONE / SMLNUM
SCALE = ONE
*
IF( LSAME( NORMIN, 'N' ) ) THEN
*
@ -343,8 +343,67 @@
IF( TMAX.LE.BIGNUM ) THEN
TSCAL = ONE
ELSE
TSCAL = ONE / ( SMLNUM*TMAX )
CALL DSCAL( N, TSCAL, CNORM, 1 )
*
* Avoid NaN generation if entries in CNORM exceed the
* overflow threshold
*
IF( TMAX.LE.DLAMCH('Overflow') ) THEN
* Case 1: All entries in CNORM are valid floating-point numbers
TSCAL = ONE / ( SMLNUM*TMAX )
CALL DSCAL( N, TSCAL, CNORM, 1 )
ELSE
* Case 2: At least one column norm of A cannot be represented
* as floating-point number. Find the offdiagonal entry A( I, J )
* with the largest absolute value. If this entry is not +/- Infinity,
* use this value as TSCAL.
TMAX = ZERO
IF( UPPER ) THEN
*
* A is upper triangular.
*
DO J = 2, N
TMAX = MAX( DLANGE( 'M', J-1, 1, A( 1, J ), 1, SUMJ ),
$ TMAX )
END DO
ELSE
*
* A is lower triangular.
*
DO J = 1, N - 1
TMAX = MAX( DLANGE( 'M', N-J, 1, A( J+1, J ), 1,
$ SUMJ ), TMAX )
END DO
END IF
*
IF( TMAX.LE.DLAMCH('Overflow') ) THEN
TSCAL = ONE / ( SMLNUM*TMAX )
DO J = 1, N
IF( CNORM( J ).LE.DLAMCH('Overflow') ) THEN
CNORM( J ) = CNORM( J )*TSCAL
ELSE
* Recompute the 1-norm without introducing Infinity
* in the summation
CNORM( J ) = ZERO
IF( UPPER ) THEN
DO I = 1, J - 1
CNORM( J ) = CNORM( J ) +
$ TSCAL * ABS( A( I, J ) )
END DO
ELSE
DO I = J + 1, N
CNORM( J ) = CNORM( J ) +
$ TSCAL * ABS( A( I, J ) )
END DO
END IF
END IF
END DO
ELSE
* At least one entry of A is not a valid floating-point entry.
* Rely on TRSV to propagate Inf and NaN.
CALL DTRSV( UPLO, TRANS, DIAG, N, A, LDA, X, 1 )
RETURN
END IF
END IF
END IF
*
* Compute a bound on the computed solution vector to see if the

2
lib/linalg/dorgbr.f
View File

@ -232,7 +232,7 @@
END IF
END IF
END IF
LWKOPT = WORK( 1 )
LWKOPT = INT( WORK( 1 ) )
LWKOPT = MAX (LWKOPT, MN)
END IF
*

190
lib/linalg/dposv.f Normal file
View File

@ -0,0 +1,190 @@
*> \brief <b> DPOSV computes the solution to system of linear equations A * X = B for PO matrices</b>
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DPOSV + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dposv.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dposv.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dposv.f">
*> [TXT]</a>
*> \endhtmlonly
*
* Definition:
* ===========
*
* SUBROUTINE DPOSV( UPLO, N, NRHS, A, LDA, B, LDB, INFO )
*
* .. Scalar Arguments ..
* CHARACTER UPLO
* INTEGER INFO, LDA, LDB, N, NRHS
* ..
* .. Array Arguments ..
* DOUBLE PRECISION A( LDA, * ), B( LDB, * )
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> DPOSV computes the solution to a real system of linear equations
*> A * X = B,
*> where A is an N-by-N symmetric positive definite matrix and X and B
*> are N-by-NRHS matrices.
*>
*> The Cholesky decomposition is used to factor A as
*> A = U**T* U, if UPLO = 'U', or
*> A = L * L**T, if UPLO = 'L',
*> where U is an upper triangular matrix and L is a lower triangular
*> matrix. The factored form of A is then used to solve the system of
*> equations A * X = B.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] UPLO
*> \verbatim
*> UPLO is CHARACTER*1
*> = 'U': Upper triangle of A is stored;
*> = 'L': Lower triangle of A is stored.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> The number of linear equations, i.e., the order of the
*> matrix A. N >= 0.
*> \endverbatim
*>
*> \param[in] NRHS
*> \verbatim
*> NRHS is INTEGER
*> The number of right hand sides, i.e., the number of columns
*> of the matrix B. NRHS >= 0.
*> \endverbatim
*>
*> \param[in,out] A
*> \verbatim
*> A is DOUBLE PRECISION array, dimension (LDA,N)
*> On entry, the symmetric matrix A. If UPLO = 'U', the leading
*> N-by-N upper triangular part of A contains the upper
*> triangular part of the matrix A, and the strictly lower
*> triangular part of A is not referenced. If UPLO = 'L', the
*> leading N-by-N lower triangular part of A contains the lower
*> triangular part of the matrix A, and the strictly upper
*> triangular part of A is not referenced.
*>
*> On exit, if INFO = 0, the factor U or L from the Cholesky
*> factorization A = U**T*U or A = L*L**T.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*> LDA is INTEGER
*> The leading dimension of the array A. LDA >= max(1,N).
*> \endverbatim
*>
*> \param[in,out] B
*> \verbatim
*> B is DOUBLE PRECISION array, dimension (LDB,NRHS)
*> On entry, the N-by-NRHS right hand side matrix B.
*> On exit, if INFO = 0, the N-by-NRHS solution matrix X.
*> \endverbatim
*>
*> \param[in] LDB
*> \verbatim
*> LDB is INTEGER
*> The leading dimension of the array B. LDB >= max(1,N).
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
*> = 0: successful exit
*> < 0: if INFO = -i, the i-th argument had an illegal value
*> > 0: if INFO = i, the leading minor of order i of A is not
*> positive definite, so the factorization could not be
*> completed, and the solution has not been computed.
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup doublePOsolve
*
* =====================================================================
SUBROUTINE DPOSV( UPLO, N, NRHS, A, LDA, B, LDB, INFO )
*
* -- LAPACK driver routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
CHARACTER UPLO
INTEGER INFO, LDA, LDB, N, NRHS
* ..
* .. Array Arguments ..
DOUBLE PRECISION A( LDA, * ), B( LDB, * )
* ..
*
* =====================================================================
*
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* ..
* .. External Subroutines ..
EXTERNAL DPOTRF, DPOTRS, XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
INFO = 0
IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
INFO = -1
ELSE IF( N.LT.0 ) THEN
INFO = -2
ELSE IF( NRHS.LT.0 ) THEN
INFO = -3
ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
INFO = -5
ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
INFO = -7
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'DPOSV ', -INFO )
RETURN
END IF
*
* Compute the Cholesky factorization A = U**T*U or A = L*L**T.
*
CALL DPOTRF( UPLO, N, A, LDA, INFO )
IF( INFO.EQ.0 ) THEN
*
* Solve the system A*X = B, overwriting B with X.
*
CALL DPOTRS( UPLO, N, NRHS, A, LDA, B, LDB, INFO )
*
END IF
RETURN
*
* End of DPOSV
*
END

201
lib/linalg/dpotrs.f Normal file
View File

@ -0,0 +1,201 @@
*> \brief \b DPOTRS
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DPOTRS + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dpotrs.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dpotrs.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dpotrs.f">
*> [TXT]</a>
*> \endhtmlonly
*
* Definition:
* ===========
*
* SUBROUTINE DPOTRS( UPLO, N, NRHS, A, LDA, B, LDB, INFO )
*
* .. Scalar Arguments ..
* CHARACTER UPLO
* INTEGER INFO, LDA, LDB, N, NRHS
* ..
* .. Array Arguments ..
* DOUBLE PRECISION A( LDA, * ), B( LDB, * )
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> DPOTRS solves a system of linear equations A*X = B with a symmetric
*> positive definite matrix A using the Cholesky factorization
*> A = U**T*U or A = L*L**T computed by DPOTRF.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] UPLO
*> \verbatim
*> UPLO is CHARACTER*1
*> = 'U': Upper triangle of A is stored;
*> = 'L': Lower triangle of A is stored.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> The order of the matrix A. N >= 0.
*> \endverbatim
*>
*> \param[in] NRHS
*> \verbatim
*> NRHS is INTEGER
*> The number of right hand sides, i.e., the number of columns
*> of the matrix B. NRHS >= 0.
*> \endverbatim
*>
*> \param[in] A
*> \verbatim
*> A is DOUBLE PRECISION array, dimension (LDA,N)
*> The triangular factor U or L from the Cholesky factorization
*> A = U**T*U or A = L*L**T, as computed by DPOTRF.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*> LDA is INTEGER
*> The leading dimension of the array A. LDA >= max(1,N).
*> \endverbatim
*>
*> \param[in,out] B
*> \verbatim
*> B is DOUBLE PRECISION array, dimension (LDB,NRHS)
*> On entry, the right hand side matrix B.
*> On exit, the solution matrix X.
*> \endverbatim
*>
*> \param[in] LDB
*> \verbatim
*> LDB is INTEGER
*> The leading dimension of the array B. LDB >= max(1,N).
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
*> = 0: successful exit
*> < 0: if INFO = -i, the i-th argument had an illegal value
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup doublePOcomputational
*
* =====================================================================
SUBROUTINE DPOTRS( UPLO, N, NRHS, A, LDA, B, LDB, INFO )
*
* -- LAPACK computational routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
CHARACTER UPLO
INTEGER INFO, LDA, LDB, N, NRHS
* ..
* .. Array Arguments ..
DOUBLE PRECISION A( LDA, * ), B( LDB, * )
* ..
*
* =====================================================================
*
* .. Parameters ..
DOUBLE PRECISION ONE
PARAMETER ( ONE = 1.0D+0 )
* ..
* .. Local Scalars ..
LOGICAL UPPER
* ..
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* ..
* .. External Subroutines ..
EXTERNAL DTRSM, XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
INFO = 0
UPPER = LSAME( UPLO, 'U' )
IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
INFO = -1
ELSE IF( N.LT.0 ) THEN
INFO = -2
ELSE IF( NRHS.LT.0 ) THEN
INFO = -3
ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
INFO = -5
ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
INFO = -7
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'DPOTRS', -INFO )
RETURN
END IF
*
* Quick return if possible
*
IF( N.EQ.0 .OR. NRHS.EQ.0 )
$ RETURN
*
IF( UPPER ) THEN
*
* Solve A*X = B where A = U**T *U.
*
* Solve U**T *X = B, overwriting B with X.
*
CALL DTRSM( 'Left', 'Upper', 'Transpose', 'Non-unit', N, NRHS,
$ ONE, A, LDA, B, LDB )
*
* Solve U*X = B, overwriting B with X.
*
CALL DTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', N,
$ NRHS, ONE, A, LDA, B, LDB )
ELSE
*
* Solve A*X = B where A = L*L**T.
*
* Solve L*X = B, overwriting B with X.
*
CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Non-unit', N,
$ NRHS, ONE, A, LDA, B, LDB )
*
* Solve L**T *X = B, overwriting B with X.
*
CALL DTRSM( 'Left', 'Lower', 'Transpose', 'Non-unit', N, NRHS,
$ ONE, A, LDA, B, LDB )
END IF
*
RETURN
*
* End of DPOTRS
*
END

5
lib/linalg/dscal.f
View File

@ -93,11 +93,14 @@
*
* .. Local Scalars ..
INTEGER I,M,MP1,NINCX
* .. Parameters ..
DOUBLE PRECISION ONE
PARAMETER (ONE=1.0D+0)
* ..
* .. Intrinsic Functions ..
INTRINSIC MOD
* ..
IF (N.LE.0 .OR. INCX.LE.0) RETURN
IF (N.LE.0 .OR. INCX.LE.0 .OR. DA.EQ.ONE) RETURN
IF (INCX.EQ.1) THEN
*
* code for increment equal to 1

2
lib/linalg/dsyevd.f
View File

@ -257,7 +257,7 @@
LWMIN = 2*N + 1
END IF
LOPT = MAX( LWMIN, 2*N +
$ ILAENV( 1, 'DSYTRD', UPLO, N, -1, -1, -1 ) )
$ N*ILAENV( 1, 'DSYTRD', UPLO, N, -1, -1, -1 ) )
LIOPT = LIWMIN
END IF
WORK( 1 ) = LOPT

4
lib/linalg/dsygvd.f
View File

@ -330,8 +330,8 @@
CALL DSYGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
CALL DSYEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, IWORK, LIWORK,
$ INFO )
LOPT = MAX( DBLE( LOPT ), DBLE( WORK( 1 ) ) )
LIOPT = MAX( DBLE( LIOPT ), DBLE( IWORK( 1 ) ) )
LOPT = INT( MAX( DBLE( LOPT ), DBLE( WORK( 1 ) ) ) )
LIOPT = INT( MAX( DBLE( LIOPT ), DBLE( IWORK( 1 ) ) ) )
*
IF( WANTZ .AND. INFO.EQ.0 ) THEN
*

2
lib/linalg/ieeeck.f
View File

@ -41,7 +41,7 @@
*> \param[in] ISPEC
*> \verbatim
*> ISPEC is INTEGER
*> Specifies whether to test just for inifinity arithmetic
*> Specifies whether to test just for infinity arithmetic
*> or whether to test for infinity and NaN arithmetic.
*> = 0: Verify infinity arithmetic only.
*> = 1: Verify infinity and NaN arithmetic.

15
lib/linalg/ilaenv.f
View File

@ -469,6 +469,15 @@
ELSE
NB = 64
END IF
ELSE IF( C3.EQ.'SYL' ) THEN
* The upper bound is to prevent overly aggressive scaling.
IF( SNAME ) THEN
NB = MIN( MAX( 48, INT( ( MIN( N1, N2 ) * 16 ) / 100) ),
$ 240 )
ELSE
NB = MIN( MAX( 24, INT( ( MIN( N1, N2 ) * 8 ) / 100) ),
$ 80 )
END IF
END IF
ELSE IF( C2.EQ.'LA' ) THEN
IF( C3.EQ.'UUM' ) THEN
@ -477,6 +486,12 @@
ELSE
NB = 64
END IF
ELSE IF( C3.EQ.'TRS' ) THEN
IF( SNAME ) THEN
NB = 32
ELSE
NB = 32
END IF
END IF
ELSE IF( SNAME .AND. C2.EQ.'ST' ) THEN
IF( C3.EQ.'EBZ' ) THEN

11
lib/linalg/zdscal.f
View File

@ -92,17 +92,20 @@
*
* .. Local Scalars ..
INTEGER I,NINCX
* .. Parameters ..
DOUBLE PRECISION ONE
PARAMETER (ONE=1.0D+0)
* ..
* .. Intrinsic Functions ..
INTRINSIC DCMPLX
INTRINSIC DBLE, DCMPLX, DIMAG
* ..
IF (N.LE.0 .OR. INCX.LE.0) RETURN
IF (N.LE.0 .OR. INCX.LE.0 .OR. DA.EQ.ONE) RETURN
IF (INCX.EQ.1) THEN
*
* code for increment equal to 1
*
DO I = 1,N
ZX(I) = DCMPLX(DA,0.0d0)*ZX(I)
ZX(I) = DCMPLX(DA*DBLE(ZX(I)),DA*DIMAG(ZX(I)))
END DO
ELSE
*
@ -110,7 +113,7 @@
*
NINCX = N*INCX
DO I = 1,NINCX,INCX
ZX(I) = DCMPLX(DA,0.0d0)*ZX(I)
ZX(I) = DCMPLX(DA*DBLE(ZX(I)),DA*DIMAG(ZX(I)))
END DO
END IF
RETURN

2
lib/linalg/zheevd.f
View File

@ -284,7 +284,7 @@
LIWMIN = 1
END IF
LOPT = MAX( LWMIN, N +
$ ILAENV( 1, 'ZHETRD', UPLO, N, -1, -1, -1 ) )
$ N*ILAENV( 1, 'ZHETRD', UPLO, N, -1, -1, -1 ) )
LROPT = LRWMIN
LIOPT = LIWMIN
END IF

2
lib/linalg/zlascl.f
View File

@ -272,6 +272,8 @@
ELSE
MUL = CTOC / CFROMC
DONE = .TRUE.
IF (MUL .EQ. ONE)
$ RETURN
END IF
END IF
*

6
lib/linalg/zscal.f
View File

@ -93,7 +93,11 @@
* .. Local Scalars ..
INTEGER I,NINCX
* ..
IF (N.LE.0 .OR. INCX.LE.0) RETURN
* .. Parameters ..
COMPLEX*16 ONE
PARAMETER (ONE= (1.0D+0,0.0D+0))
* ..
IF (N.LE.0 .OR. INCX.LE.0 .OR. ZA.EQ.ONE) RETURN
IF (INCX.EQ.1) THEN
*
* code for increment equal to 1

2
src/CG-SPICA/angle_spica.cpp
View File

@ -297,7 +297,7 @@ void AngleSPICA::init_style()
repflag = 0;
for (int i = 1; i <= atom->nangletypes; i++)
if (repscale[i] > 0.0) repflag = 1;
if (repscale && (repscale[i] > 0.0)) repflag = 1;
// set up pointers to access SPICA LJ parameters for 1-3 interactions

56
src/MANYBODY/pair_sw.cpp
View File

@ -46,6 +46,7 @@ PairSW::PairSW(LAMMPS *lmp) : Pair(lmp)
centroidstressflag = CENTROID_NOTAVAIL;
unit_convert_flag = utils::get_supported_conversions(utils::ENERGY);
skip_threebody_flag = false;
params_mapped = 0;
params = nullptr;
@ -225,12 +226,12 @@ void PairSW::compute(int eflag, int vflag)
void PairSW::allocate()
{
allocated = 1;
int n = atom->ntypes;
int np1 = atom->ntypes + 1;
memory->create(setflag,n+1,n+1,"pair:setflag");
memory->create(cutsq,n+1,n+1,"pair:cutsq");
memory->create(neighshort,maxshort,"pair:neighshort");
map = new int[n+1];
memory->create(setflag, np1, np1, "pair:setflag");
memory->create(cutsq, np1, np1, "pair:cutsq");
memory->create(neighshort, maxshort, "pair:neighshort");
map = new int[np1];
}
/* ----------------------------------------------------------------------
@ -262,12 +263,49 @@ void PairSW::coeff(int narg, char **arg)
{
if (!allocated) allocate();
map_element2type(narg-3,arg+3);
// read potential file and set up element maps only once
if (one_coeff || !params_mapped) {
// make certain that the setflag array is always fully initialized
// the sw/intel pair style depends on it
if (!one_coeff) {
for (int i = 0; i <= atom->ntypes; i++) {
for (int j = 0; j <= atom->ntypes; j++) {
setflag[i][j] = 0;
}
}
}
// read potential file and initialize potential parameters
map_element2type(narg-3, arg+3, (one_coeff != 0));
read_file(arg[2]);
setup_params();
// read potential file and initialize potential parameters
read_file(arg[2]);
setup_params();
params_mapped = 1;
}
if (!one_coeff) {
int ilo, ihi, jlo, jhi;
utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error);
utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
if (((map[i] >= 0) && (strcmp(arg[i+2], elements[map[i]]) != 0)) ||
((map[i] < 0) && (strcmp(arg[i+2], "NULL") != 0)))
error->all(FLERR, "Must use consistent type to element mappings with threebody off");
if (map[i] < 0) error->all(FLERR, "Must not set pair_coeff mapped to NULL element");
for (int j = MAX(jlo, i); j <= jhi; j++) {
if (((map[j] >= 0) && (strcmp(arg[j+2], elements[map[j]]) != 0)) ||
((map[j] < 0) && (strcmp(arg[j+2], "NULL") != 0)))
error->all(FLERR, "Must use consistent type to element mappings with threebody off");
if (map[j] < 0) error->all(FLERR, "Must not set pair_coeff mapped to NULL element");
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients");
}
}
/* ----------------------------------------------------------------------

1
src/MANYBODY/pair_sw.h
View File

@ -54,6 +54,7 @@ class PairSW : public Pair {
int maxshort; // size of short neighbor list array
int *neighshort; // short neighbor list array
int skip_threebody_flag; // whether to run threebody loop
int params_mapped; // whether parameters have been read and mapped to elements
void settings(int, char **) override;
virtual void allocate();

115
src/MISC/fix_ipi.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -47,12 +46,12 @@ using namespace FixConst;
// socket interface
#ifndef _WIN32
#include <unistd.h>
#include <sys/types.h>
#include <sys/socket.h>
#include <netinet/in.h>
#include <sys/un.h>
#include <netdb.h>
#include <netinet/in.h>
#include <sys/socket.h>
#include <sys/types.h>
#include <sys/un.h>
#include <unistd.h>
#else
#ifndef WIN32_LEAN_AND_MEAN
#define WIN32_LEAN_AND_MEAN
@ -65,8 +64,7 @@ using namespace FixConst;
/* Utility functions to simplify the interface with POSIX sockets */
static void open_socket(int &sockfd, int inet, int port, char* host,
Error *error)
static void open_socket(int &sockfd, int inet, int port, char *host, Error *error)
/* Opens a socket.
Args:
@ -83,9 +81,9 @@ static void open_socket(int &sockfd, int inet, int port, char* host,
int ai_err;
#ifdef _WIN32
error->one(FLERR,"i-PI socket implementation requires UNIX environment");
error->one(FLERR, "i-PI socket implementation requires UNIX environment");
#else
if (inet>0) { // creates an internet socket
if (inet > 0) { // creates an internet socket
// fetches information on the host
struct addrinfo hints, *res;
@ -96,40 +94,39 @@ static void open_socket(int &sockfd, int inet, int port, char* host,
hints.ai_flags = AI_PASSIVE;
ai_err = getaddrinfo(host, std::to_string(port).c_str(), &hints, &res);
if (ai_err!=0)
error->one(FLERR,"Error fetching host data. Wrong host name?");
if (ai_err != 0) error->one(FLERR, "Error fetching host data. Wrong host name?");
// creates socket
sockfd = socket(res->ai_family, res->ai_socktype, res->ai_protocol);
if (sockfd < 0)
error->one(FLERR,"Error opening socket");
if (sockfd < 0) error->one(FLERR, "Error opening socket");
// makes connection
if (connect(sockfd, res->ai_addr, res->ai_addrlen) < 0)
error->one(FLERR,"Error opening INET socket: wrong port or server unreachable");
error->one(FLERR, "Error opening INET socket: wrong port or server unreachable");
freeaddrinfo(res);
} else { // creates a unix socket
} else { // creates a unix socket
struct sockaddr_un serv_addr;
// fills up details of the socket address
memset(&serv_addr, 0, sizeof(serv_addr));
serv_addr.sun_family = AF_UNIX;
strcpy(serv_addr.sun_path, "/tmp/ipi_");
strcpy(serv_addr.sun_path+9, host);
strcpy(serv_addr.sun_path + 9, host);
// creates the socket
sockfd = socket(AF_UNIX, SOCK_STREAM, 0);
// connects
if (connect(sockfd, (struct sockaddr *) &serv_addr, sizeof(serv_addr)) < 0)
error->one(FLERR,"Error opening UNIX socket: server may not be running "
error->one(FLERR,
"Error opening UNIX socket: server may not be running "
"or the path to the socket unavailable");
}
#endif
}
static void writebuffer(int sockfd, const char *data, int len, Error* error)
static void writebuffer(int sockfd, const char *data, int len, Error *error)
/* Writes to a socket.
Args:
@ -140,14 +137,11 @@ static void writebuffer(int sockfd, const char *data, int len, Error* error)
{
int n;
n = write(sockfd,data,len);
if (n < 0)
error->one(FLERR,"Error writing to socket: broken connection");
n = write(sockfd, data, len);
if (n < 0) error->one(FLERR, "Error writing to socket: broken connection");
}
static void readbuffer(int sockfd, char *data, int len, Error* error)
static void readbuffer(int sockfd, char *data, int len, Error *error)
/* Reads from a socket.
Args:
@ -158,44 +152,52 @@ static void readbuffer(int sockfd, char *data, int len, Error* error)
{
int n, nr;
n = nr = read(sockfd,data,len);
n = nr = read(sockfd, data, len);
while (nr>0 && n<len) {
nr=read(sockfd,&data[n],len-n);
n+=nr;
while (nr > 0 && n < len) {
nr = read(sockfd, &data[n], len - n);
n += nr;
}
if (n == 0)
error->one(FLERR,"Error reading from socket: broken connection");
if (n == 0) error->one(FLERR, "Error reading from socket: broken connection");
}
/* ---------------------------------------------------------------------- */
FixIPI::FixIPI(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), irregular(nullptr)
FixIPI::FixIPI(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), irregular(nullptr)
{
/* format for fix:
* fix num group_id ipi host port [unix]
*/
if (strcmp(style,"ipi") != 0 && narg < 5)
error->all(FLERR,"Illegal fix ipi command");
if (narg < 5) utils::missing_cmd_args(FLERR, "fix ipi", error);
if (atom->tag_enable == 0)
error->all(FLERR,"Cannot use fix ipi without atom IDs");
if (atom->tag_enable == 0) error->all(FLERR, "Cannot use fix ipi without atom IDs");
if (atom->tag_consecutive() == 0) error->all(FLERR, "Fix ipi requires consecutive atom IDs");
if (strcmp(update->unit_style, "lj") == 0) error->all(FLERR, "Fix ipi does not support lj units");
if (atom->tag_consecutive() == 0)
error->all(FLERR,"Fix ipi requires consecutive atom IDs");
if (strcmp(arg[1],"all") != 0)
error->warning(FLERR,"Fix ipi always uses group all");
if ((strcmp(arg[1], "all") != 0) && (comm->me == 0))
error->warning(FLERR, "Not using group 'all' with fix ipi can result in undefined behavior");
host = strdup(arg[3]);
port = utils::inumeric(FLERR,arg[4],false,lmp);
port = utils::inumeric(FLERR, arg[4], false, lmp);
inet = ((narg > 5) && (strcmp(arg[5],"unix") == 0) ) ? 0 : 1;
master = (comm->me==0) ? 1 : 0;
// check if forces should be reinitialized and set flag
reset_flag = ((narg > 6 && (strcmp(arg[5],"reset") == 0 )) || ((narg > 5) && (strcmp(arg[5],"reset") == 0)) ) ? 1 : 0;
master = (comm->me == 0) ? 1 : 0;
inet = 1;
reset_flag = 0;
int iarg = 5;
while (iarg < narg) {
if (strcmp(arg[iarg], "unix") == 0) {
inet = 0;
++iarg;
} else if (strcmp(arg[iarg], "reset") == 0) {
reset_flag = 1;
++iarg;
} else {
error->all(FLERR, "Unknown fix ipi keyword: {}", arg[iarg]);
}
}
// sanity check
if (inet && ((port <= 1024) || (port > 65536)))
error->all(FLERR, "Invalid port for fix ipi: {}", port);
hasdata = bsize = 0;
@ -223,7 +225,6 @@ FixIPI::~FixIPI()
delete irregular;
}
/* ---------------------------------------------------------------------- */
int FixIPI::setmask()
@ -240,9 +241,11 @@ void FixIPI::init()
{
//only opens socket on master process
if (master) {
if (!socketflag) open_socket(ipisock, inet, port, host, error);
} else ipisock=0;
//! should check for success in socket opening -- but the current open_socket routine dies brutally if unsuccessful
if (!socketflag) open_socket(ipisock, inet, port, host, error);
} else
ipisock = 0;
// TODO: should check for success in socket opening,
// but the current open_socket routine dies brutally if unsuccessful
// tell lammps we have assigned a socket
socketflag = 1;
@ -252,11 +255,13 @@ void FixIPI::init()
kspace_flag = (force->kspace) ? 1 : 0;
// makes sure that neighbor lists are re-built at each step (cannot make assumptions when cycling over beads!)
// makes sure that neighbor lists are re-built at each step
// (cannot make assumptions when cycling over beads!)
neighbor->delay = 0;
neighbor->every = 1;
}
// clang-format off
void FixIPI::initial_integrate(int /*vflag*/)
{
/* This is called at the beginning of the integration loop,

4
src/REACTION/fix_bond_react.cpp
View File

@ -39,7 +39,7 @@ Contributing Author: Jacob Gissinger (jacob.r.gissinger@gmail.com)
#include "neighbor.h"
#include "pair.h"
#include "random_mars.h"
#include "reset_mol_ids.h"
#include "reset_atoms_mol.h"
#include "respa.h"
#include "update.h"
#include "variable.h"
@ -207,7 +207,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
if (reset_mol_ids_flag) {
delete reset_mol_ids;
reset_mol_ids = new ResetMolIDs(lmp);
reset_mol_ids = new ResetAtomsMol(lmp);
reset_mol_ids->create_computes(id,group->names[igroup]);
}

2
src/REACTION/fix_bond_react.h
View File

@ -110,7 +110,7 @@ class FixBondReact : public Fix {
class RanMars **random; // random number for 'prob' keyword
class RanMars **rrhandom; // random number for Arrhenius constraint
class NeighList *list;
class ResetMolIDs *reset_mol_ids; // class for resetting mol IDs
class ResetAtomsMol *reset_mol_ids; // class for resetting mol IDs
int *reacted_mol, *unreacted_mol;
int *limit_duration; // indicates how long to relax

3
src/REPLICA/fix_pimd.cpp
View File

@ -102,6 +102,9 @@ FixPIMD::FixPIMD(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
error->universe_all(FLERR, fmt::format("Unknown keyword {} for fix pimd", arg[i]));
}
if (strcmp(update->unit_style, "lj") == 0)
error->all(FLERR, "Fix pimd does not support lj units");
/* Initiation */
size_peratom_cols = 12 * nhc_nchain + 3;

184
src/REPLICA/neb.cpp
View File

@ -108,9 +108,9 @@ NEB::~NEB()
void NEB::command(int narg, char **arg)
{
if (domain->box_exist == 0)
error->all(FLERR,"NEB command before simulation box is defined");
error->universe_all(FLERR,"NEB command before simulation box is defined");
if (narg < 6) error->universe_all(FLERR,"Illegal NEB command");
if (narg < 6) error->universe_all(FLERR,"Illegal NEB command: missing argument(s)");
etol = utils::numeric(FLERR,arg[0],false,lmp);
ftol = utils::numeric(FLERR,arg[1],false,lmp);
@ -120,11 +120,14 @@ void NEB::command(int narg, char **arg)
// error checks
if (etol < 0.0) error->all(FLERR,"Illegal NEB command");
if (ftol < 0.0) error->all(FLERR,"Illegal NEB command");
if (nevery <= 0) error->universe_all(FLERR,"Illegal NEB command");
if (n1steps % nevery || n2steps % nevery)
error->universe_all(FLERR,"Illegal NEB command");
if (etol < 0.0) error->universe_all(FLERR, fmt::format("Illegal NEB energy tolerance: {}", etol));
if (ftol < 0.0) error->universe_all(FLERR, fmt::format("Illegal NEB force tolerance: {}", ftol));
if (nevery <= 0)
error->universe_all(FLERR,fmt::format("Illegal NEB command every parameter: {}", nevery));
if (n1steps % nevery)
error->all(FLERR, fmt::format("NEB N1 value {} incompatible with every {}", n1steps, nevery));
if (n2steps % nevery)
error->all(FLERR, fmt::format("NEB N2 value {} incompatible with every {}", n2steps, nevery));
// replica info
@ -136,26 +139,38 @@ void NEB::command(int narg, char **arg)
// error checks
if (nreplica == 1) error->all(FLERR,"Cannot use NEB with a single replica");
if (nreplica == 1) error->universe_all(FLERR,"Cannot use NEB with a single replica");
if (atom->map_style == Atom::MAP_NONE)
error->all(FLERR,"Cannot use NEB unless atom map exists");
error->universe_all(FLERR,"Cannot use NEB without an atom map");
// process file-style setting to setup initial configs for all replicas
if (strcmp(arg[5],"final") == 0) {
if (narg != 7 && narg !=8) error->universe_all(FLERR,"Illegal NEB command");
inpfile = arg[6];
readfile(inpfile,0);
} else if (strcmp(arg[5],"each") == 0) {
if (narg != 7 && narg !=8) error->universe_all(FLERR,"Illegal NEB command");
inpfile = arg[6];
readfile(inpfile,1);
} else if (strcmp(arg[5],"none") == 0) {
if (narg != 6 && narg !=7) error->universe_all(FLERR,"Illegal NEB command");
} else error->universe_all(FLERR,"Illegal NEB command");
int iarg = 5;
int filecmd = 0;
verbose=false;
if (strcmp(arg[narg-1],"verbose") == 0) verbose=true;
while (iarg < narg) {
if (strcmp(arg[iarg],"final") == 0) {
if (iarg + 2 > narg) error->universe_all(FLERR,"Illegal NEB command: missing arguments");
inpfile = arg[iarg+1];
readfile(inpfile,0);
filecmd = 1;
iarg += 2;
} else if (strcmp(arg[iarg],"each") == 0) {
if (iarg + 2 > narg) error->universe_all(FLERR,"Illegal NEB command: missing arguments");
inpfile = arg[iarg+1];
readfile(inpfile,1);
filecmd = 1;
iarg += 2;
} else if (strcmp(arg[iarg],"none") == 0) {
filecmd = 1;
++iarg;
} else if (strcmp(arg[iarg],"verbose") == 0) {
verbose=true;
++iarg;
} else error->universe_all(FLERR,fmt::format("Unknown NEB command keyword: {}", arg[iarg]));
}
if (!filecmd) error->universe_all(FLERR, "NEB is missing 'final', 'each', or 'none' keyword");
// run the NEB calculation
run();
@ -176,7 +191,7 @@ void NEB::run()
auto fixes = modify->get_fix_by_style("^neb$");
if (fixes.size() != 1)
error->all(FLERR,"NEB requires use of exactly one fix neb instance");
error->universe_all(FLERR,"NEB requires use of exactly one fix neb instance");
fneb = dynamic_cast<FixNEB *>(fixes[0]);
if (verbose) numall =7;
@ -194,7 +209,7 @@ void NEB::run()
lmp->init();
if (update->minimize->searchflag)
error->all(FLERR,"NEB requires damped dynamics minimizer");
error->universe_all(FLERR,"NEB requires a damped dynamics minimizer");
// setup regular NEB minimization
FILE *uscreen = universe->uscreen;
@ -207,38 +222,43 @@ void NEB::run()
update->endstep = update->laststep = update->firststep + n1steps;
update->nsteps = n1steps;
update->max_eval = n1steps;
if (update->laststep < 0)
error->all(FLERR,"Too many timesteps for NEB");
if (update->laststep < 0) error->universe_all(FLERR,"Too many timesteps for NEB");
update->minimize->setup();
if (me_universe == 0) {
if (uscreen) {
fmt::print(uscreen," Step {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} ",
"MaxReplicaForce", "MaxAtomForce", "GradV0","GradV1","GradVc","EBF", "EBR", "RDT");
for (int i = 1; i <= nreplica; ++i)
fmt::print(uscreen, "{:^14} {:^14} ", "RD"+std::to_string(i), "PE"+std::to_string(i));
if (verbose) {
fprintf(uscreen,"Step MaxReplicaForce MaxAtomForce "
"GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... "
"RDN PEN pathangle1 angletangrad1 anglegrad1 gradV1 "
"ReplicaForce1 MaxAtomForce1 pathangle2 angletangrad2 "
"... ReplicaForceN MaxAtomForceN\n");
} else {
fprintf(uscreen,"Step MaxReplicaForce MaxAtomForce "
"GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... "
"RDN PEN\n");
for (int i = 1; i <= nreplica; ++i) {
auto idx = std::to_string(i);
fmt::print(uscreen, "{:^12}{:^12}{:^12} {:^12} {:^12}{:^12} ",
"pathangle"+idx, "angletangrad"+idx, "anglegrad"+idx, "gradV"+idx,
"RepForce"+idx, "MaxAtomForce"+idx);
}
}
fprintf(uscreen,"\n");
}
if (ulogfile) {
fmt::print(ulogfile," Step {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} ",
"MaxReplicaForce", "MaxAtomForce", "GradV0","GradV1","GradVc","EBF", "EBR", "RDT");
for (int i = 1; i <= nreplica; ++i)
fmt::print(ulogfile, "{:^14} {:^14} ", "RD"+std::to_string(i), "PE"+std::to_string(i));
if (verbose) {
fprintf(ulogfile,"Step MaxReplicaForce MaxAtomForce "
"GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... "
"RDN PEN pathangle1 angletangrad1 anglegrad1 gradV1 "
"ReplicaForce1 MaxAtomForce1 pathangle2 angletangrad2 "
"... ReplicaForceN MaxAtomForceN\n");
} else {
fprintf(ulogfile,"Step MaxReplicaForce MaxAtomForce "
"GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... "
"RDN PEN\n");
for (int i = 1; i <= nreplica; ++i) {
auto idx = std::to_string(i);
fmt::print(ulogfile, "{:^12}{:^12}{:^12} {:^12} {:^12}{:^12} ",
"pathangle"+idx, "angletangrad"+idx, "anglegrad"+idx, "gradV"+idx,
"RepForce"+idx, "MaxAtomForce"+idx);
}
}
fprintf(ulogfile,"\n");
}
}
print_status();
@ -292,8 +312,7 @@ void NEB::run()
update->endstep = update->laststep = update->firststep + n2steps;
update->nsteps = n2steps;
update->max_eval = n2steps;
if (update->laststep < 0)
error->all(FLERR,"Too many timesteps");
if (update->laststep < 0) error->universe_all(FLERR,"Too many timesteps");
update->minimize->init();
fneb->rclimber = top;
@ -301,34 +320,37 @@ void NEB::run()
if (me_universe == 0) {
if (uscreen) {
fmt::print(uscreen," Step {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} ",
"MaxReplicaForce", "MaxAtomForce", "GradV0","GradV1","GradVc","EBF", "EBR", "RDT");
for (int i = 1; i <= nreplica; ++i)
fmt::print(uscreen, "{:^14} {:^14} ", "RD"+std::to_string(i), "PE"+std::to_string(i));
if (verbose) {
fprintf(uscreen,"Step MaxReplicaForce MaxAtomForce "
"GradV0 GradV1 GradVc EBF EBR RDT "
"RD1 PE1 RD2 PE2 ... RDN PEN "
"pathangle1 angletangrad1 anglegrad1 gradV1 "
"ReplicaForce1 MaxAtomForce1 pathangle2 angletangrad2 "
"... ReplicaForceN MaxAtomForceN\n");
} else {
fprintf(uscreen,"Step MaxReplicaForce MaxAtomForce "
"GradV0 GradV1 GradVc "
"EBF EBR RDT "
"RD1 PE1 RD2 PE2 ... RDN PEN\n");
for (int i = 1; i <= nreplica; ++i) {
auto idx = std::to_string(i);
fmt::print(uscreen, "{:^12}{:^12}{:^12} {:^12} {:^12}{:^12} ",
"pathangle"+idx, "angletangrad"+idx, "anglegrad"+idx, "gradV"+idx,
"RepForce"+idx, "MaxAtomForce"+idx);
}
}
fprintf(uscreen,"\n");
}
if (ulogfile) {
fmt::print(ulogfile," Step {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} ",
"MaxReplicaForce", "MaxAtomForce", "GradV0","GradV1","GradVc","EBF", "EBR", "RDT");
for (int i = 1; i <= nreplica; ++i)
fmt::print(ulogfile, "{:^14} {:^14} ", "RD"+std::to_string(i), "PE"+std::to_string(i));
if (verbose) {
fprintf(ulogfile,"Step MaxReplicaForce MaxAtomForce "
"GradV0 GradV1 GradVc EBF EBR RDT "
"RD1 PE1 RD2 PE2 ... RDN PEN "
"pathangle1 angletangrad1 anglegrad1 gradV1 "
"ReplicaForce1 MaxAtomForce1 pathangle2 angletangrad2 "
"... ReplicaForceN MaxAtomForceN\n");
} else {
fprintf(ulogfile,"Step MaxReplicaForce MaxAtomForce "
"GradV0 GradV1 GradVc "
"EBF EBR RDT "
"RD1 PE1 RD2 PE2 ... RDN PEN\n");
for (int i = 1; i <= nreplica; ++i) {
auto idx = std::to_string(i);
fmt::print(ulogfile, "{:^12}{:^12}{:^12} {:^12} {:^12}{:^12} ",
"pathangle"+idx, "angletangrad"+idx, "anglegrad"+idx, "gradV"+idx,
"RepForce"+idx, "MaxAtomForce"+idx);
}
}
fprintf(ulogfile,"\n");
}
}
print_status();
@ -420,7 +442,7 @@ void NEB::readfile(char *file, int flag)
}
MPI_Bcast(&nlines,1,MPI_INT,0,world);
if (nlines < 0)
error->all(FLERR,"Incorrectly formatted NEB file");
error->universe_all(FLERR,"Incorrectly formatted NEB file");
}
auto buffer = new char[CHUNK*MAXLINE];
@ -542,11 +564,11 @@ void NEB::open(char *file)
if (platform::has_compress_extension(file)) {
compressed = 1;
fp = platform::compressed_read(file);
if (!fp) error->one(FLERR,"Cannot open compressed file");
if (!fp) error->one(FLERR,"Cannot open compressed file {}: {}", file, utils::getsyserror());
} else fp = fopen(file,"r");
if (fp == nullptr)
error->one(FLERR,"Cannot open file {}: {}",file,utils::getsyserror());
error->one(FLERR,"Cannot open file {}: {}", file, utils::getsyserror());
}
/* ----------------------------------------------------------------------
@ -626,19 +648,19 @@ void NEB::print_status()
if (me_universe == 0) {
constexpr double todeg=180.0/MY_PI;
std::string mesg = fmt::format("{} {:12.8g} {:12.8g} ",update->ntimestep,fmaxreplica,fmaxatom);
mesg += fmt::format("{:12.8g} {:12.8g} {:12.8g} ",gradvnorm0,gradvnorm1,gradvnormc);
mesg += fmt::format("{:12.8g} {:12.8g} {:12.8g} ",ebf,ebr,endpt);
for (int i = 0; i < nreplica; i++) mesg += fmt::format("{:12.8g} {:12.8g} ",rdist[i],all[i][0]);
std::string mesg = fmt::format("{:10} {:<14.8g} {:<14.8g} ",update->ntimestep,fmaxreplica,fmaxatom);
mesg += fmt::format("{:<14.8g} {:<14.8g} {:<14.8g} ",gradvnorm0,gradvnorm1,gradvnormc);
mesg += fmt::format("{:<14.8g} {:<14.8g} {:<14.8g} ",ebf,ebr,endpt);
for (int i = 0; i < nreplica; i++) mesg += fmt::format("{:<14.8g} {:<14.8g} ",rdist[i],all[i][0]);
if (verbose) {
mesg += fmt::format("{:12.5g} {:12.5g} {:12.5g} {:12.5g} {:12.5g} {:12.5g}",
mesg += fmt::format("{:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g}",
NAN,180-acos(all[0][5])*todeg,180-acos(all[0][6])*todeg,
all[0][3],freplica[0],fmaxatomInRepl[0]);
for (int i = 1; i < nreplica-1; i++)
mesg += fmt::format("{:12.5g} {:12.5g} {:12.5g} {:12.5g} {:12.5g} {:12.5g}",
mesg += fmt::format("{:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g}",
180-acos(all[i][4])*todeg,180-acos(all[i][5])*todeg,
180-acos(all[i][6])*todeg,all[i][3],freplica[i],fmaxatomInRepl[i]);
mesg += fmt::format("{:12.5g} {:12.5g} {:12.5g} {:12.5g} {:12.5g} {:12.5g}",
mesg += fmt::format("{:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g}",
NAN,180-acos(all[nreplica-1][5])*todeg,NAN,all[nreplica-1][3],
freplica[nreplica-1],fmaxatomInRepl[nreplica-1]);
}
@ -650,4 +672,8 @@ void NEB::print_status()
fflush(universe->ulogfile);
}
}
if (verbose) {
delete[] freplica;
delete[] fmaxatomInRepl;
}
}

10
src/angle_hybrid.cpp
View File

@ -306,6 +306,16 @@ void AngleHybrid::coeff(int narg, char **arg)
void AngleHybrid::init_style()
{
// error if sub-style is not used
int used;
for (int istyle = 0; istyle < nstyles; ++istyle) {
used = 0;
for (int itype = 1; itype <= atom->nangletypes; ++itype)
if (map[itype] == istyle) used = 1;
if (used == 0) error->all(FLERR, "Angle hybrid sub-style {} is not used", keywords[istyle]);
}
for (int m = 0; m < nstyles; m++)
if (styles[m]) styles[m]->init_style();
}

10
src/bond_hybrid.cpp
View File

@ -334,6 +334,16 @@ void BondHybrid::coeff(int narg, char **arg)
void BondHybrid::init_style()
{
// error if sub-style is not used
int used;
for (int istyle = 0; istyle < nstyles; ++istyle) {
used = 0;
for (int itype = 1; itype <= atom->nbondtypes; ++itype)
if (map[itype] == istyle) used = 1;
if (used == 0) error->all(FLERR, "Bond hybrid sub-style {} is not used", keywords[istyle]);
}
for (int m = 0; m < nstyles; m++)
if (styles[m]) styles[m]->init_style();

23
src/deprecated.cpp
View File

@ -36,16 +36,25 @@ void Deprecated::command(int narg, char **arg)
if (lmp->comm->me == 0)
utils::logmesg(lmp, "\nThe 'box' command has been removed and will be ignored\n\n");
return;
} else if (cmd == "reset_ids") {
if (lmp->comm->me == 0)
utils::logmesg(lmp, "\n'reset_ids' has been renamed to 'reset_atom_ids'\n\n");
} else if (utils::strmatch(cmd, "^kim_")) {
if (lmp->comm->me == 0)
utils::logmesg(lmp,
"\nWARNING: 'kim_<command>' has been renamed to 'kim <command>'. "
"Please update your input.\n\n");
std::string newcmd("kim");
newcmd += " " + cmd.substr(4);
if (lmp->comm->me == 0)
utils::logmesg(lmp, "\nWARNING: '{}' has been renamed to '{}'. Please update your input.\n\n",
cmd, newcmd);
for (int i = 0; i < narg; ++i) {
newcmd.append(1, ' ');
newcmd.append(arg[i]);
}
input->one(newcmd);
return;
} else if (utils::strmatch(cmd, "^reset_")) {
std::string newcmd("reset_atoms");
if ((cmd == "reset_ids") || (cmd == "reset_atom_ids")) newcmd += " id";
if (cmd == "reset_mol_ids") newcmd += " mol";
if (lmp->comm->me == 0)
utils::logmesg(lmp, "\nWARNING: '{}' has been renamed to '{}'. Please update your input.\n\n",
cmd, newcmd);
for (int i = 0; i < narg; ++i) {
newcmd.append(1, ' ');
newcmd.append(arg[i]);

4
src/deprecated.h
View File

@ -15,12 +15,14 @@
// clang-format off
CommandStyle(DEPRECATED,Deprecated);
CommandStyle(box,Deprecated);
CommandStyle(reset_ids,Deprecated);
CommandStyle(kim_init,Deprecated);
CommandStyle(kim_interactions,Deprecated);
CommandStyle(kim_param,Deprecated);
CommandStyle(kim_property,Deprecated);
CommandStyle(kim_query,Deprecated);
CommandStyle(reset_ids,Deprecated);
CommandStyle(reset_atom_ids,Deprecated);
CommandStyle(reset_mol_ids,Deprecated);
CommandStyle(message,Deprecated);
CommandStyle(server,Deprecated);
// clang-format on

10
src/dihedral_hybrid.cpp
View File

@ -313,6 +313,16 @@ void DihedralHybrid::coeff(int narg, char **arg)
void DihedralHybrid::init_style()
{
// error if sub-style is not used
int used;
for (int istyle = 0; istyle < nstyles; ++istyle) {
used = 0;
for (int itype = 1; itype <= atom->ndihedraltypes; ++itype)
if (map[itype] == istyle) used = 1;
if (used == 0) error->all(FLERR, "Dihedral hybrid sub-style {} is not used", keywords[istyle]);
}
for (int m = 0; m < nstyles; m++)
if (styles[m]) styles[m]->init_style();
}

10
src/improper_hybrid.cpp
View File

@ -305,6 +305,16 @@ void ImproperHybrid::coeff(int narg, char **arg)
void ImproperHybrid::init_style()
{
// error if sub-style is not used
int used;
for (int istyle = 0; istyle < nstyles; ++istyle) {
used = 0;
for (int itype = 1; itype <= atom->nimpropertypes; ++itype)
if (map[itype] == istyle) used = 1;
if (used == 0) error->all(FLERR, "Improper hybrid sub-style {} is not used", keywords[istyle]);
}
for (int m = 0; m < nstyles; m++)
if (styles[m]) styles[m]->init_style();
}

166
src/input.cpp
View File

@ -749,76 +749,77 @@ int Input::execute_command()
{
int flag = 1;
if (!strcmp(command,"clear")) clear();
else if (!strcmp(command,"echo")) echo();
else if (!strcmp(command,"if")) ifthenelse();
else if (!strcmp(command,"include")) include();
else if (!strcmp(command,"jump")) jump();
else if (!strcmp(command,"label")) label();
else if (!strcmp(command,"log")) log();
else if (!strcmp(command,"next")) next_command();
else if (!strcmp(command,"partition")) partition();
else if (!strcmp(command,"print")) print();
else if (!strcmp(command,"python")) python();
else if (!strcmp(command,"quit")) quit();
else if (!strcmp(command,"shell")) shell();
else if (!strcmp(command,"variable")) variable_command();
std::string mycmd = command;
if (mycmd == "clear") clear();
else if (mycmd == "echo") echo();
else if (mycmd == "if") ifthenelse();
else if (mycmd == "include") include();
else if (mycmd == "jump") jump();
else if (mycmd == "label") label();
else if (mycmd == "log") log();
else if (mycmd == "next") next_command();
else if (mycmd == "partition") partition();
else if (mycmd == "print") print();
else if (mycmd == "python") python();
else if (mycmd == "quit") quit();
else if (mycmd == "shell") shell();
else if (mycmd == "variable") variable_command();
else if (!strcmp(command,"angle_coeff")) angle_coeff();
else if (!strcmp(command,"angle_style")) angle_style();
else if (!strcmp(command,"atom_modify")) atom_modify();
else if (!strcmp(command,"atom_style")) atom_style();
else if (!strcmp(command,"bond_coeff")) bond_coeff();
else if (!strcmp(command,"bond_style")) bond_style();
else if (!strcmp(command,"bond_write")) bond_write();
else if (!strcmp(command,"boundary")) boundary();
else if (!strcmp(command,"comm_modify")) comm_modify();
else if (!strcmp(command,"comm_style")) comm_style();
else if (!strcmp(command,"compute")) compute();
else if (!strcmp(command,"compute_modify")) compute_modify();
else if (!strcmp(command,"dielectric")) dielectric();
else if (!strcmp(command,"dihedral_coeff")) dihedral_coeff();
else if (!strcmp(command,"dihedral_style")) dihedral_style();
else if (!strcmp(command,"dimension")) dimension();
else if (!strcmp(command,"dump")) dump();
else if (!strcmp(command,"dump_modify")) dump_modify();
else if (!strcmp(command,"fix")) fix();
else if (!strcmp(command,"fix_modify")) fix_modify();
else if (!strcmp(command,"group")) group_command();
else if (!strcmp(command,"improper_coeff")) improper_coeff();
else if (!strcmp(command,"improper_style")) improper_style();
else if (!strcmp(command,"kspace_modify")) kspace_modify();
else if (!strcmp(command,"kspace_style")) kspace_style();
else if (!strcmp(command,"labelmap")) labelmap();
else if (!strcmp(command,"lattice")) lattice();
else if (!strcmp(command,"mass")) mass();
else if (!strcmp(command,"min_modify")) min_modify();
else if (!strcmp(command,"min_style")) min_style();
else if (!strcmp(command,"molecule")) molecule();
else if (!strcmp(command,"neigh_modify")) neigh_modify();
else if (!strcmp(command,"neighbor")) neighbor_command();
else if (!strcmp(command,"newton")) newton();
else if (!strcmp(command,"package")) package();
else if (!strcmp(command,"pair_coeff")) pair_coeff();
else if (!strcmp(command,"pair_modify")) pair_modify();
else if (!strcmp(command,"pair_style")) pair_style();
else if (!strcmp(command,"pair_write")) pair_write();
else if (!strcmp(command,"processors")) processors();
else if (!strcmp(command,"region")) region();
else if (!strcmp(command,"reset_timestep")) reset_timestep();
else if (!strcmp(command,"restart")) restart();
else if (!strcmp(command,"run_style")) run_style();
else if (!strcmp(command,"special_bonds")) special_bonds();
else if (!strcmp(command,"suffix")) suffix();
else if (!strcmp(command,"thermo")) thermo();
else if (!strcmp(command,"thermo_modify")) thermo_modify();
else if (!strcmp(command,"thermo_style")) thermo_style();
else if (!strcmp(command,"timestep")) timestep();
else if (!strcmp(command,"timer")) timer_command();
else if (!strcmp(command,"uncompute")) uncompute();
else if (!strcmp(command,"undump")) undump();
else if (!strcmp(command,"unfix")) unfix();
else if (!strcmp(command,"units")) units();
else if (mycmd == "angle_coeff") angle_coeff();
else if (mycmd == "angle_style") angle_style();
else if (mycmd == "atom_modify") atom_modify();
else if (mycmd == "atom_style") atom_style();
else if (mycmd == "bond_coeff") bond_coeff();
else if (mycmd == "bond_style") bond_style();
else if (mycmd == "bond_write") bond_write();
else if (mycmd == "boundary") boundary();
else if (mycmd == "comm_modify") comm_modify();
else if (mycmd == "comm_style") comm_style();
else if (mycmd == "compute") compute();
else if (mycmd == "compute_modify") compute_modify();
else if (mycmd == "dielectric") dielectric();
else if (mycmd == "dihedral_coeff") dihedral_coeff();
else if (mycmd == "dihedral_style") dihedral_style();
else if (mycmd == "dimension") dimension();
else if (mycmd == "dump") dump();
else if (mycmd == "dump_modify") dump_modify();
else if (mycmd == "fix") fix();
else if (mycmd == "fix_modify") fix_modify();
else if (mycmd == "group") group_command();
else if (mycmd == "improper_coeff") improper_coeff();
else if (mycmd == "improper_style") improper_style();
else if (mycmd == "kspace_modify") kspace_modify();
else if (mycmd == "kspace_style") kspace_style();
else if (mycmd == "labelmap") labelmap();
else if (mycmd == "lattice") lattice();
else if (mycmd == "mass") mass();
else if (mycmd == "min_modify") min_modify();
else if (mycmd == "min_style") min_style();
else if (mycmd == "molecule") molecule();
else if (mycmd == "neigh_modify") neigh_modify();
else if (mycmd == "neighbor") neighbor_command();
else if (mycmd == "newton") newton();
else if (mycmd == "package") package();
else if (mycmd == "pair_coeff") pair_coeff();
else if (mycmd == "pair_modify") pair_modify();
else if (mycmd == "pair_style") pair_style();
else if (mycmd == "pair_write") pair_write();
else if (mycmd == "processors") processors();
else if (mycmd == "region") region();
else if (mycmd == "reset_timestep") reset_timestep();
else if (mycmd == "restart") restart();
else if (mycmd == "run_style") run_style();
else if (mycmd == "special_bonds") special_bonds();
else if (mycmd == "suffix") suffix();
else if (mycmd == "thermo") thermo();
else if (mycmd == "thermo_modify") thermo_modify();
else if (mycmd == "thermo_style") thermo_style();
else if (mycmd == "timestep") timestep();
else if (mycmd == "timer") timer_command();
else if (mycmd == "uncompute") uncompute();
else if (mycmd == "undump") undump();
else if (mycmd == "unfix") unfix();
else if (mycmd == "units") units();
else flag = 0;
@ -826,10 +827,15 @@ int Input::execute_command()
if (flag) return 0;
// process "meta-commands", i.e. commands that may have sub-commands
// they return 1 if there was a match and 0 if not
if (mycmd == "reset_atoms") flag = meta(mycmd);
if (flag) return 0;
// invoke commands added via style_command.h
// try suffixed version first
std::string mycmd = command;
if (lmp->suffix_enable && lmp->non_pair_suffix()) {
mycmd = command + std::string("/") + lmp->non_pair_suffix();
if (command_map->find(mycmd) == command_map->end()) {
@ -1718,7 +1724,7 @@ void Input::pair_coeff()
if (force->pair == nullptr) error->all(FLERR,"Pair_coeff command without a pair style");
if (narg < 2) utils::missing_cmd_args(FLERR,"pair_coeff", error);
if (force->pair->one_coeff && ((strcmp(arg[0],"*") != 0) || (strcmp(arg[1],"*") != 0)))
error->all(FLERR,"Pair_coeff must start with * * for this pair style");
error->all(FLERR,"Pair_coeff must start with * * for pair style {}", force->pair_style);
char *newarg0 = utils::expand_type(FLERR, arg[0], Atom::ATOM, lmp);
if (newarg0) arg[0] = newarg0;
@ -1978,3 +1984,21 @@ void Input::units()
error->all(FLERR,"Units command after simulation box is defined");
update->set_units(arg[0]);
}
/* ---------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
function for meta commands
------------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */
int Input::meta(const std::string &prefix)
{
auto mycmd = fmt::format("{}_{}", utils::uppercase(prefix), utils::uppercase(arg[0]));
if (command_map->find(mycmd) != command_map->end()) {
CommandCreator &command_creator = (*command_map)[mycmd];
Command *cmd = command_creator(lmp);
cmd->command(narg-1,arg+1);
delete cmd;
return 1;
} else return 0;
}

4
src/input.h
View File

@ -70,6 +70,8 @@ class Input : protected Pointers {
void reallocate(char *&, int &, int); // reallocate a char string
int execute_command(); // execute a single command
int meta(const std::string &); // process meta-commands
void clear(); // input script commands
void echo();
void ifthenelse();
@ -142,7 +144,5 @@ class Input : protected Pointers {
void unfix();
void units();
};
} // namespace LAMMPS_NS
#endif

2
src/pair.cpp
View File

@ -821,7 +821,7 @@ void Pair::map_element2type(int narg, char **arg, bool update_setflag)
// elements = list of element names
if (narg != ntypes)
error->all(FLERR,"Incorrect args for pair coefficients");
error->all(FLERR,"Number of element to type mappings does not match number of atom types");
if (elements) {
for (i = 0; i < nelements; i++) delete[] elements[i];

19
src/pair_hybrid.cpp
View File

@ -271,10 +271,9 @@ void PairHybrid::allocate()
void PairHybrid::settings(int narg, char **arg)
{
if (narg < 1) error->all(FLERR,"Illegal pair_style command");
if (narg < 1) utils::missing_cmd_args(FLERR, "pair_style hybrid", error);
if (lmp->kokkos && !utils::strmatch(force->pair_style,"^hybrid.*/kk$"))
error->all(FLERR,"Must use pair_style {}/kk with Kokkos",
force->pair_style);
error->all(FLERR,"Must use pair_style {}/kk with Kokkos", force->pair_style);
// delete old lists, since cannot just change settings
@ -326,9 +325,9 @@ void PairHybrid::settings(int narg, char **arg)
nstyles = 0;
while (iarg < narg) {
if (utils::strmatch(arg[iarg],"^hybrid"))
error->all(FLERR,"Pair style hybrid cannot have hybrid as an argument");
error->all(FLERR,"Pair style hybrid cannot have hybrid as a sub-style");
if (strcmp(arg[iarg],"none") == 0)
error->all(FLERR,"Pair style hybrid cannot have none as an argument");
error->all(FLERR,"Pair style hybrid cannot have none as a sub-style");
styles[nstyles] = force->new_pair(arg[iarg],1,dummy);
keywords[nstyles] = force->store_style(arg[iarg],0);
@ -478,7 +477,7 @@ void PairHybrid::init_svector()
void PairHybrid::coeff(int narg, char **arg)
{
if (narg < 3) error->all(FLERR,"Incorrect args for pair coefficients");
if (narg < 3) utils::missing_cmd_args(FLERR,"pair_coeff", error);
if (!allocated) allocate();
int ilo,ihi,jlo,jhi;
@ -497,7 +496,7 @@ void PairHybrid::coeff(int narg, char **arg)
if (strcmp(arg[2],keywords[m]) == 0) {
if (multiple[m]) {
multflag = 1;
if (narg < 4) error->all(FLERR,"Incorrect args for pair coefficients");
if (narg < 4) utils::missing_cmd_args(FLERR, "pair_coeff", error);
if (multiple[m] == utils::inumeric(FLERR,arg[3],false,lmp)) break;
else continue;
} else break;
@ -507,7 +506,7 @@ void PairHybrid::coeff(int narg, char **arg)
int none = 0;
if (m == nstyles) {
if (strcmp(arg[2],"none") == 0) none = 1;
else error->all(FLERR,"Pair coeff for hybrid has invalid style: {}",arg[2]);
else error->all(FLERR,"Pair coeff for hybrid has invalid style: {}", arg[2]);
}
// move 1st/2nd args to 2nd/3rd args
@ -527,7 +526,7 @@ void PairHybrid::coeff(int narg, char **arg)
if (!none && styles[m]->one_coeff) {
if ((strcmp(arg[0],"*") != 0) || (strcmp(arg[1],"*") != 0))
error->all(FLERR,"Incorrect args for pair coefficients");
error->all(FLERR,"Pair_coeff must start with * * for sub-style {}", keywords[m]);
for (int i = 1; i <= atom->ntypes; i++)
for (int j = i; j <= atom->ntypes; j++)
if (nmap[i][j] && map[i][j][0] == m) {
@ -578,7 +577,7 @@ void PairHybrid::init_style()
for (jtype = itype; jtype <= ntypes; jtype++)
for (m = 0; m < nmap[itype][jtype]; m++)
if (map[itype][jtype][m] == istyle) used = 1;
if (used == 0) error->all(FLERR,"Pair hybrid sub-style is not used");
if (used == 0) error->all(FLERR,"Pair hybrid sub-style {} is not used", keywords[istyle]);
}
// The GPU library uses global data for each pair style, so the

250
src/reset_atom_ids.cpp → src/reset_atoms_id.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -12,7 +11,7 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "reset_atom_ids.h"
#include "reset_atoms_id.h"
#include "atom.h"
#include "atom_vec.h"
@ -32,7 +31,7 @@ using namespace LAMMPS_NS;
#if defined(LMP_QSORT)
// allocate space for static class variable
// prototype for non-class function
ResetIDs::AtomRvous *ResetIDs::sortrvous;
ResetAtomsID::AtomRvous *ResetAtomsID::sortrvous;
static int compare_coords(const void *, const void *);
#else
// prototype for non-class function
@ -44,22 +43,21 @@ static int compare_coords(const int, const int, void *);
/* ---------------------------------------------------------------------- */
ResetIDs::ResetIDs(LAMMPS *lmp) : Command(lmp) {}
ResetAtomsID::ResetAtomsID(LAMMPS *lmp) : Command(lmp) {}
/* ---------------------------------------------------------------------- */
void ResetIDs::command(int narg, char **arg)
void ResetAtomsID::command(int narg, char **arg)
{
if (domain->box_exist == 0)
error->all(FLERR,"Reset_ids command before simulation box is defined");
if (atom->tag_enable == 0)
error->all(FLERR,"Cannot use reset_atom_ids unless atoms have IDs");
error->all(FLERR, "Reset_atoms id command before simulation box is defined");
if (atom->tag_enable == 0) error->all(FLERR, "Cannot use reset_atoms id unless atoms have IDs");
for (int i = 0; i < modify->nfix; i++)
if (modify->fix[i]->stores_ids)
error->all(FLERR,"Cannot use reset_atom_ids when a fix exists that stores atom IDs");
error->all(FLERR, "Cannot use reset_atoms id when a fix exists that stores atom IDs");
if (comm->me == 0) utils::logmesg(lmp,"Resetting atom IDs ...\n");
if (comm->me == 0) utils::logmesg(lmp, "Resetting atom IDs ...\n");
// process args
@ -67,11 +65,12 @@ void ResetIDs::command(int narg, char **arg)
int iarg = 0;
while (iarg < narg) {
if (strcmp(arg[iarg],"sort") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal reset_atom_ids command");
sortflag = utils::logical(FLERR,arg[iarg+1],false,lmp);
if (strcmp(arg[iarg], "sort") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "reset_atoms id", error);
sortflag = utils::logical(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else error->all(FLERR,"Illegal reset_atom_ids command");
} else
error->all(FLERR, "Unknown reset_atoms id keyword: {}", arg[iarg]);
}
// create an atom map if one doesn't exist already
@ -99,7 +98,7 @@ void ResetIDs::command(int narg, char **arg)
comm->setup();
comm->exchange();
comm->borders();
if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
if (domain->triclinic) domain->lamda2x(atom->nlocal + atom->nghost);
// oldIDs = copy of current owned IDs
@ -108,7 +107,7 @@ void ResetIDs::command(int narg, char **arg)
int nall = nlocal + atom->nghost;
tagint *oldIDs;
memory->create(oldIDs,nlocal,"reset_atom_ids:oldIDs");
memory->create(oldIDs, nlocal, "reset_atom_ids:oldIDs");
for (int i = 0; i < nlocal; i++) {
oldIDs[i] = tag[i];
@ -119,22 +118,24 @@ void ResetIDs::command(int narg, char **arg)
// if sortflag = no: use fast tag_extend()
// if sortflag = yes: use slower full spatial sort plus rendezvous comm
if (sortflag == 0) atom->tag_extend();
else sort();
if (sortflag == 0)
atom->tag_extend();
else
sort();
// newIDs = copy of new IDs
// restore old IDs, consistent with existing atom map
// forward_comm_array acquires new IDs for ghost atoms
double **newIDs;
memory->create(newIDs,nall,1,"reset_atom_ids:newIDs");
memory->create(newIDs, nall, 1, "reset_atom_ids:newIDs");
for (int i = 0; i < nlocal; i++) {
newIDs[i][0] = ubuf(tag[i]).d;
tag[i] = oldIDs[i];
}
comm->forward_comm_array(1,newIDs);
comm->forward_comm_array(1, newIDs);
// loop over bonds, angles, etc and reset IDs in stored topology arrays
// only necessary for molecular = Atom::MOLECULAR, not molecular = Atom::TEMPLATE
@ -143,7 +144,7 @@ void ResetIDs::command(int narg, char **arg)
int badcount = 0;
if (atom->molecular == Atom::MOLECULAR) {
int j,m;
int j, m;
tagint oldID;
if (atom->avec->bonds_allow) {
@ -153,8 +154,10 @@ void ResetIDs::command(int narg, char **arg)
for (j = 0; j < num_bond[i]; j++) {
oldID = bond_atom[i][j];
m = atom->map(oldID);
if (m >= 0) bond_atom[i][j] = (tagint) ubuf(newIDs[m][0]).i;
else badcount++;
if (m >= 0)
bond_atom[i][j] = (tagint) ubuf(newIDs[m][0]).i;
else
badcount++;
}
}
}
@ -168,18 +171,24 @@ void ResetIDs::command(int narg, char **arg)
for (j = 0; j < num_angle[i]; j++) {
oldID = angle_atom1[i][j];
m = atom->map(oldID);
if (m >= 0) angle_atom1[i][j] = (tagint) ubuf(newIDs[m][0]).i;
else badcount++;
if (m >= 0)
angle_atom1[i][j] = (tagint) ubuf(newIDs[m][0]).i;
else
badcount++;
oldID = angle_atom2[i][j];
m = atom->map(oldID);
if (m >= 0) angle_atom2[i][j] = (tagint) ubuf(newIDs[m][0]).i;
else badcount++;
if (m >= 0)
angle_atom2[i][j] = (tagint) ubuf(newIDs[m][0]).i;
else
badcount++;
oldID = angle_atom3[i][j];
m = atom->map(oldID);
if (m >= 0) angle_atom3[i][j] = (tagint) ubuf(newIDs[m][0]).i;
else badcount++;
if (m >= 0)
angle_atom3[i][j] = (tagint) ubuf(newIDs[m][0]).i;
else
badcount++;
}
}
}
@ -194,23 +203,31 @@ void ResetIDs::command(int narg, char **arg)
for (j = 0; j < num_dihedral[i]; j++) {
oldID = dihedral_atom1[i][j];
m = atom->map(oldID);
if (m >= 0) dihedral_atom1[i][j] = (tagint) ubuf(newIDs[m][0]).i;
else badcount++;
if (m >= 0)
dihedral_atom1[i][j] = (tagint) ubuf(newIDs[m][0]).i;
else
badcount++;
oldID = dihedral_atom2[i][j];
m = atom->map(oldID);
if (m >= 0) dihedral_atom2[i][j] = (tagint) ubuf(newIDs[m][0]).i;
else badcount++;
if (m >= 0)
dihedral_atom2[i][j] = (tagint) ubuf(newIDs[m][0]).i;
else
badcount++;
oldID = dihedral_atom3[i][j];
m = atom->map(oldID);
if (m >= 0) dihedral_atom3[i][j] = (tagint) ubuf(newIDs[m][0]).i;
else badcount++;
if (m >= 0)
dihedral_atom3[i][j] = (tagint) ubuf(newIDs[m][0]).i;
else
badcount++;
oldID = dihedral_atom4[i][j];
m = atom->map(oldID);
if (m >= 0) dihedral_atom4[i][j] = (tagint) ubuf(newIDs[m][0]).i;
else badcount++;
if (m >= 0)
dihedral_atom4[i][j] = (tagint) ubuf(newIDs[m][0]).i;
else
badcount++;
}
}
}
@ -225,23 +242,31 @@ void ResetIDs::command(int narg, char **arg)
for (j = 0; j < num_improper[i]; j++) {
oldID = improper_atom1[i][j];
m = atom->map(oldID);
if (m >= 0) improper_atom1[i][j] = (tagint) ubuf(newIDs[m][0]).i;
else badcount++;
if (m >= 0)
improper_atom1[i][j] = (tagint) ubuf(newIDs[m][0]).i;
else
badcount++;
oldID = improper_atom2[i][j];
m = atom->map(oldID);
if (m >= 0) improper_atom2[i][j] = (tagint) ubuf(newIDs[m][0]).i;
else badcount++;
if (m >= 0)
improper_atom2[i][j] = (tagint) ubuf(newIDs[m][0]).i;
else
badcount++;
oldID = improper_atom3[i][j];
m = atom->map(oldID);
if (m >= 0) improper_atom3[i][j] = (tagint) ubuf(newIDs[m][0]).i;
else badcount++;
if (m >= 0)
improper_atom3[i][j] = (tagint) ubuf(newIDs[m][0]).i;
else
badcount++;
oldID = improper_atom4[i][j];
m = atom->map(oldID);
if (m >= 0) improper_atom4[i][j] = (tagint) ubuf(newIDs[m][0]).i;
else badcount++;
if (m >= 0)
improper_atom4[i][j] = (tagint) ubuf(newIDs[m][0]).i;
else
badcount++;
}
}
}
@ -250,10 +275,10 @@ void ResetIDs::command(int narg, char **arg)
// error check
int all;
MPI_Allreduce(&badcount,&all,1,MPI_INT,MPI_SUM,world);
MPI_Allreduce(&badcount, &all, 1, MPI_INT, MPI_SUM, world);
if (all)
error->all(FLERR,"Reset_ids missing {} bond topology atom IDs - "
"use comm_modify cutoff",all);
error->all(FLERR,
"reset_atoms id is missing {} bond topology atom IDs - use comm_modify cutoff", all);
// reset IDs and atom map for owned atoms
@ -287,9 +312,9 @@ void ResetIDs::command(int narg, char **arg)
spatial sort of atoms followed by rendezvous comm to reset atom IDs
------------------------------------------------------------------------- */
void ResetIDs::sort()
void ResetAtomsID::sort()
{
double mylo[3],myhi[3],bboxlo[3],bboxhi[3];
double mylo[3], myhi[3], bboxlo[3], bboxhi[3];
int me = comm->me;
int nprocs = comm->nprocs;
@ -306,12 +331,12 @@ void ResetIDs::sort()
myhi[0] = myhi[1] = myhi[2] = -BIG;
for (int i = 0; i < nlocal; i++) {
mylo[0] = MIN(mylo[0],x[i][0]);
mylo[1] = MIN(mylo[1],x[i][1]);
mylo[2] = MIN(mylo[2],x[i][2]);
myhi[0] = MAX(myhi[0],x[i][0]);
myhi[1] = MAX(myhi[1],x[i][1]);
myhi[2] = MAX(myhi[2],x[i][2]);
mylo[0] = MIN(mylo[0], x[i][0]);
mylo[1] = MIN(mylo[1], x[i][1]);
mylo[2] = MIN(mylo[2], x[i][2]);
myhi[0] = MAX(myhi[0], x[i][0]);
myhi[1] = MAX(myhi[1], x[i][1]);
myhi[2] = MAX(myhi[2], x[i][2]);
}
if (dim == 2) mylo[2] = myhi[2] = 0.0;
@ -325,15 +350,15 @@ void ResetIDs::sort()
}
}
MPI_Allreduce(mylo,bboxlo,3,MPI_DOUBLE,MPI_MIN,world);
MPI_Allreduce(myhi,bboxhi,3,MPI_DOUBLE,MPI_MAX,world);
MPI_Allreduce(mylo, bboxlo, 3, MPI_DOUBLE, MPI_MIN, world);
MPI_Allreduce(myhi, bboxhi, 3, MPI_DOUBLE, MPI_MAX, world);
bboxlo[0] -= 0.0001 * (bboxhi[0]-bboxlo[0]);
bboxlo[1] -= 0.0001 * (bboxhi[1]-bboxlo[1]);
bboxlo[2] -= 0.0001 * (bboxhi[2]-bboxlo[2]);
bboxhi[0] += 0.0001 * (bboxhi[0]-bboxlo[0]);
bboxhi[1] += 0.0001 * (bboxhi[1]-bboxlo[1]);
bboxhi[2] += 0.0001 * (bboxhi[2]-bboxlo[2]);
bboxlo[0] -= 0.0001 * (bboxhi[0] - bboxlo[0]);
bboxlo[1] -= 0.0001 * (bboxhi[1] - bboxlo[1]);
bboxlo[2] -= 0.0001 * (bboxhi[2] - bboxlo[2]);
bboxhi[0] += 0.0001 * (bboxhi[0] - bboxlo[0]);
bboxhi[1] += 0.0001 * (bboxhi[1] - bboxlo[1]);
bboxhi[2] += 0.0001 * (bboxhi[2] - bboxlo[2]);
// nbin_estimate = estimate of total number of bins, each with PERBIN atoms
// binsize = edge length of a cubic bin
@ -342,19 +367,23 @@ void ResetIDs::sort()
bigint nbin_estimate = atom->natoms / PERBIN + 1;
double vol;
if (dim == 2) vol = (bboxhi[0]-bboxlo[0]) * (bboxhi[1]-bboxlo[1]);
else vol = (bboxhi[0]-bboxlo[0]) * (bboxhi[1]-bboxlo[1]) * (bboxhi[2]-bboxlo[2]);
double binsize = pow(vol/nbin_estimate,1.0/dim);
if (dim == 2)
vol = (bboxhi[0] - bboxlo[0]) * (bboxhi[1] - bboxlo[1]);
else
vol = (bboxhi[0] - bboxlo[0]) * (bboxhi[1] - bboxlo[1]) * (bboxhi[2] - bboxlo[2]);
double binsize = pow(vol / nbin_estimate, 1.0 / dim);
int nbinx = static_cast<int> ((bboxhi[0]-bboxlo[0]) / binsize) + 1;
int nbiny = static_cast<int> ((bboxhi[1]-bboxlo[1]) / binsize) + 1;
int nbinz = static_cast<int> ((bboxhi[2]-bboxlo[2]) / binsize) + 1;
int nbinx = static_cast<int>((bboxhi[0] - bboxlo[0]) / binsize) + 1;
int nbiny = static_cast<int>((bboxhi[1] - bboxlo[1]) / binsize) + 1;
int nbinz = static_cast<int>((bboxhi[2] - bboxlo[2]) / binsize) + 1;
double invx = 1.0 / (bboxhi[0]-bboxlo[0]);
double invy = 1.0 / (bboxhi[1]-bboxlo[1]);
double invx = 1.0 / (bboxhi[0] - bboxlo[0]);
double invy = 1.0 / (bboxhi[1] - bboxlo[1]);
double invz;
if (dim == 2) invz = 0.0;
else invz = 1.0 / (bboxhi[2]-bboxlo[2]);
if (dim == 2)
invz = 0.0;
else
invz = 1.0 / (bboxhi[2] - bboxlo[2]);
// nbins = actual total number of bins
// bins are split evenly across procs
@ -366,18 +395,18 @@ void ResetIDs::sort()
// binlo to binhi-1 = bin indices this proc will own in Rvous decomp
// bins are numbered from 0 to Nbins-1
bigint nbins = (bigint) nbinx*nbiny*nbinz;
bigint nbins = (bigint) nbinx * nbiny * nbinz;
bigint nlo = nbins / nprocs;
bigint nhi = nlo + 1;
bigint nplo = nprocs - (nbins % nprocs);
bigint nbinlo = nplo*nlo;
bigint nbinlo = nplo * nlo;
if (me < nplo) {
binlo = me * nlo;
binhi = (me+1) * nlo;
binhi = (me + 1) * nlo;
} else {
binlo = nbinlo + (me-nplo) * nhi;
binhi = nbinlo + (me+1-nplo) * nhi;
binlo = nbinlo + (me - nplo) * nhi;
binhi = nbinlo + (me + 1 - nplo) * nhi;
}
// fill atombuf with info on my atoms
@ -385,20 +414,23 @@ void ResetIDs::sort()
// proclist = proc that owns ibin
int *proclist;
memory->create(proclist,nlocal,"special:proclist");
auto atombuf = (AtomRvous *) memory->smalloc((bigint) nlocal*sizeof(AtomRvous),"resetIDs:idbuf");
memory->create(proclist, nlocal, "special:proclist");
auto atombuf =
(AtomRvous *) memory->smalloc((bigint) nlocal * sizeof(AtomRvous), "resetIDs:idbuf");
int ibinx,ibiny,ibinz,iproc;
int ibinx, ibiny, ibinz, iproc;
bigint ibin;
for (int i = 0; i < nlocal; i++) {
ibinx = static_cast<int> ((x[i][0]-bboxlo[0])*invx * nbinx);
ibiny = static_cast<int> ((x[i][1]-bboxlo[1])*invy * nbiny);
ibinz = static_cast<int> ((x[i][2]-bboxlo[2])*invz * nbinz);
ibin = (bigint) ibinz*nbiny*nbinx + (bigint) ibiny*nbinx + ibinx;
ibinx = static_cast<int>((x[i][0] - bboxlo[0]) * invx * nbinx);
ibiny = static_cast<int>((x[i][1] - bboxlo[1]) * invy * nbiny);
ibinz = static_cast<int>((x[i][2] - bboxlo[2]) * invz * nbinz);
ibin = (bigint) ibinz * nbiny * nbinx + (bigint) ibiny * nbinx + ibinx;
if (ibin < nbinlo) iproc = ibin / nlo;
else iproc = nplo + (ibin-nbinlo) / nhi;
if (ibin < nbinlo)
iproc = ibin / nlo;
else
iproc = nplo + (ibin - nbinlo) / nhi;
proclist[i] = iproc;
atombuf[i].ibin = ibin;
@ -412,8 +444,8 @@ void ResetIDs::sort()
// perform rendezvous operation, send atombuf to other procs
char *buf;
int nreturn = comm->rendezvous(1,nlocal,(char *) atombuf,sizeof(AtomRvous),0,proclist,
sort_bins,0,buf,sizeof(IDRvous),(void *) this);
int nreturn = comm->rendezvous(1, nlocal, (char *) atombuf, sizeof(AtomRvous), 0, proclist,
sort_bins, 0, buf, sizeof(IDRvous), (void *) this);
auto outbuf = (IDRvous *) buf;
memory->destroy(proclist);
@ -437,11 +469,11 @@ void ResetIDs::sort()
outbuf = list of N IDRvous datums, sent back to sending proc
------------------------------------------------------------------------- */
int ResetIDs::sort_bins(int n, char *inbuf, int &flag, int *&proclist, char *&outbuf, void *ptr)
int ResetAtomsID::sort_bins(int n, char *inbuf, int &flag, int *&proclist, char *&outbuf, void *ptr)
{
int i,ibin,index;
int i, ibin, index;
auto rptr = (ResetIDs *) ptr;
auto rptr = (ResetAtomsID *) ptr;
Memory *memory = rptr->memory;
Error *error = rptr->error;
MPI_Comm world = rptr->world;
@ -451,13 +483,13 @@ int ResetIDs::sort_bins(int n, char *inbuf, int &flag, int *&proclist, char *&ou
// nbins = my subset of bins from binlo to binhi-1
// initialize linked lists of my Rvous atoms in each bin
int *next,*head,*last,*count;
int *next, *head, *last, *count;
int nbins = binhi - binlo;
memory->create(next,n,"resetIDs:next");
memory->create(head,nbins,"resetIDs:head");
memory->create(last,nbins,"resetIDs:last");
memory->create(count,nbins,"resetIDs:count");
memory->create(next, n, "resetIDs:next");
memory->create(head, nbins, "resetIDs:head");
memory->create(last, nbins, "resetIDs:last");
memory->create(count, nbins, "resetIDs:count");
for (i = 0; i < n; i++) next[i] = -1;
@ -472,7 +504,7 @@ int ResetIDs::sort_bins(int n, char *inbuf, int &flag, int *&proclist, char *&ou
if (in[i].ibin < binlo || in[i].ibin >= binhi) {
//printf("BAD me %d i %d %d ibin %d binlohi %d %d\n",
// rptr->comm->me,i,n,in[i].ibin,binlo,binhi);
error->one(FLERR,"Bad spatial bin assignment in reset_atom_ids sort");
error->one(FLERR, "Bad spatial bin assignment in reset_atoms id sort");
}
ibin = in[i].ibin - binlo;
if (head[ibin] < 0) head[ibin] = i;
@ -482,11 +514,10 @@ int ResetIDs::sort_bins(int n, char *inbuf, int &flag, int *&proclist, char *&ou
}
int maxcount = 0;
for (ibin = 0; ibin < nbins; ibin++)
maxcount = MAX(maxcount,count[ibin]);
for (ibin = 0; ibin < nbins; ibin++) maxcount = MAX(maxcount, count[ibin]);
int *order;
memory->create(order,maxcount,"resetIDs:order");
memory->create(order, maxcount, "resetIDs:order");
// sort atoms in each bin by their spatial position
// recreate linked list for each bin based on sorted order
@ -498,12 +529,11 @@ int ResetIDs::sort_bins(int n, char *inbuf, int &flag, int *&proclist, char *&ou
index = next[index];
}
#if defined(LMP_QSORT)
sortrvous = in;
qsort(order,count[ibin],sizeof(int),compare_coords);
qsort(order, count[ibin], sizeof(int), compare_coords);
#else
utils::merge_sort(order,count[ibin],(void *) in,compare_coords);
utils::merge_sort(order, count[ibin], (void *) in, compare_coords);
#endif
head[ibin] = last[ibin] = -1;
@ -518,15 +548,15 @@ int ResetIDs::sort_bins(int n, char *inbuf, int &flag, int *&proclist, char *&ou
tagint ntag = n;
tagint nprev;
MPI_Scan(&ntag,&nprev,1,MPI_LMP_TAGINT,MPI_SUM,world);
MPI_Scan(&ntag, &nprev, 1, MPI_LMP_TAGINT, MPI_SUM, world);
nprev -= n;
// loop over bins and sorted atoms in each bin, reset ids to be consecutive
// pass outbuf of IDRvous datums back to comm->rendezvous
int nout = n;
memory->create(proclist,nout,"resetIDs:proclist");
auto out = (IDRvous *) memory->smalloc(nout*sizeof(IDRvous),"resetIDs:out");
memory->create(proclist, nout, "resetIDs:proclist");
auto out = (IDRvous *) memory->smalloc(nout * sizeof(IDRvous), "resetIDs:out");
tagint one = nprev + 1;
for (ibin = 0; ibin < nbins; ibin++) {
@ -567,7 +597,7 @@ int compare_coords(const void *iptr, const void *jptr)
{
int i = *((int *) iptr);
int j = *((int *) jptr);
ResetIDs::AtomRvous *rvous = ResetIDs::sortrvous;
ResetAtomsID::AtomRvous *rvous = ResetAtomsID::sortrvous;
double *xi = rvous[i].x;
double *xj = rvous[j].x;
if (xi[0] < xj[0]) return -1;
@ -588,7 +618,7 @@ int compare_coords(const void *iptr, const void *jptr)
int compare_coords(const int i, const int j, void *ptr)
{
auto rvous = (ResetIDs::AtomRvous *) ptr;
auto rvous = (ResetAtomsID::AtomRvous *) ptr;
double *xi = rvous[i].x;
double *xj = rvous[j].x;
if (xi[0] < xj[0]) return -1;

12
src/reset_atom_ids.h → src/reset_atoms_id.h
View File

@ -13,18 +13,18 @@
#ifdef COMMAND_CLASS
// clang-format off
CommandStyle(reset_atom_ids,ResetIDs);
CommandStyle(RESET_ATOMS_ID,ResetAtomsID);
// clang-format on
#else
#ifndef LMP_RESET_IDS_H
#define LMP_RESET_IDS_H
#ifndef LMP_RESET_ATOMS_ID_H
#define LMP_RESET_ATOMS_ID_H
#include "command.h"
namespace LAMMPS_NS {
class ResetIDs : public Command {
class ResetAtomsID : public Command {
public:
struct AtomRvous {
bigint ibin;
@ -42,7 +42,7 @@ class ResetIDs : public Command {
static AtomRvous *sortrvous;
#endif
ResetIDs(class LAMMPS *);
ResetAtomsID(class LAMMPS *);
void command(int, char **) override;
private:
@ -54,8 +54,6 @@ class ResetIDs : public Command {
void sort();
};
} // namespace LAMMPS_NS
#endif
#endif

141
src/reset_atoms_image.cpp Normal file
View File

@ -0,0 +1,141 @@
/* ----------------------------------------------------------------------
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 "reset_atoms_image.h"
#include "atom.h"
#include "atom_vec.h"
#include "comm.h"
#include "compute.h"
#include "domain.h"
#include "error.h"
#include "group.h"
#include "input.h"
#include "modify.h"
#include "variable.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ResetAtomsImage::ResetAtomsImage(LAMMPS *lmp) : Command(lmp) {}
/* ----------------------------------------------------------------------
called as reset_atoms image command in input script
------------------------------------------------------------------------- */
void ResetAtomsImage::command(int narg, char **arg)
{
if (domain->box_exist == 0)
error->all(FLERR, "Reset_atoms image command before simulation box is defined");
if (atom->tag_enable == 0)
error->all(FLERR, "Cannot use reset_atoms image unless atoms have IDs");
if (atom->avec->bonds_allow == 0)
error->all(FLERR, "Cannot use reset_atoms image used when bonds are not allowed");
// process args
if (narg < 1) utils::missing_cmd_args(FLERR, "reset_atoms image", error);
if (narg > 1) error->all(FLERR, "Unknown reset_atoms image keyword: {}", arg[1]);
int igroup = group->find(arg[0]);
if (igroup < 0) error->all(FLERR, "Could not find reset_atoms image group {}", arg[0]);
int groupbit = group->bitmask[igroup];
if (comm->me == 0) utils::logmesg(lmp, "Resetting image flags ...\n");
// record wall time for resetting molecule IDs
MPI_Barrier(world);
double time1 = platform::walltime();
// initialize system since comm->borders() will be invoked
lmp->init();
// setup domain, communication
// exchange will clear map, borders will reset
// this is the map needed to lookup current global IDs for bond topology
if (domain->triclinic) domain->x2lamda(atom->nlocal);
domain->pbc();
domain->reset_box();
comm->setup();
comm->exchange();
comm->borders();
if (domain->triclinic) domain->lamda2x(atom->nlocal + atom->nghost);
// create computes and variables
auto frags = modify->add_compute("frags_r_i_f all fragment/atom single yes");
auto chunk = modify->add_compute("chunk_r_i_f all chunk/atom c_frags_r_i_f compress yes");
auto flags = modify->add_compute("flags_r_i_f all property/atom ix iy iz");
input->variable->set("ix_r_i_f atom c_flags_r_i_f[1]");
input->variable->set("iy_r_i_f atom c_flags_r_i_f[2]");
input->variable->set("iz_r_i_f atom c_flags_r_i_f[3]");
auto ifmin = modify->add_compute("ifmin_r_i_f all reduce/chunk chunk_r_i_f min "
"v_ix_r_i_f v_iy_r_i_f v_iz_r_i_f");
auto ifmax = modify->add_compute("ifmax_r_i_f all reduce/chunk chunk_r_i_f max "
"v_ix_r_i_f v_iy_r_i_f v_iz_r_i_f");
auto cdist = modify->add_compute("cdist_r_i_f all chunk/spread/atom chunk_r_i_f "
"c_ifmax_r_i_f[*] c_ifmin_r_i_f[*]");
// trigger computes
frags->compute_peratom();
chunk->compute_peratom();
flags->compute_peratom();
ifmin->compute_array();
ifmax->compute_array();
cdist->compute_peratom();
// reset image flags for atoms in group
const int *const mask = atom->mask;
const int nlocal = atom->nlocal;
imageint *image = atom->image;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
int ix = (image[i] & IMGMASK) - IMGMAX;
int iy = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
int iz = (image[i] >> IMG2BITS) - IMGMAX;
ix -= (int) floor(0.5 * (cdist->array_atom[i][0] + cdist->array_atom[i][3]));
iy -= (int) floor(0.5 * (cdist->array_atom[i][1] + cdist->array_atom[i][4]));
iz -= (int) floor(0.5 * (cdist->array_atom[i][2] + cdist->array_atom[i][5]));
image[i] = ((imageint) (ix + IMGMAX) & IMGMASK) |
(((imageint) (iy + IMGMAX) & IMGMASK) << IMGBITS) |
(((imageint) (iz + IMGMAX) & IMGMASK) << IMG2BITS);
}
}
// cleanup
modify->delete_compute("cdist_r_i_f");
modify->delete_compute("ifmax_r_i_f");
modify->delete_compute("ifmin_r_i_f");
modify->delete_compute("flags_r_i_f");
modify->delete_compute("chunk_r_i_f");
modify->delete_compute("frags_r_i_f");
input->variable->set("ix_r_i_f delete");
input->variable->set("iy_r_i_f delete");
input->variable->set("iz_r_i_f delete");
// total time
MPI_Barrier(world);
if (comm->me == 0)
utils::logmesg(lmp, " reset_atoms image CPU = {:.3f} seconds\n", platform::walltime() - time1);
}

34
src/reset_atoms_image.h Normal file
View File

@ -0,0 +1,34 @@
/* -*- 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 COMMAND_CLASS
// clang-format off
CommandStyle(RESET_ATOMS_IMAGE,ResetAtomsImage);
// clang-format on
#else
#ifndef LMP_RESET_ATOMS_IMAGE_H
#define LMP_RESET_ATOMS_IMAGE_H
#include "command.h"
namespace LAMMPS_NS {
class ResetAtomsImage : public Command {
public:
ResetAtomsImage(class LAMMPS *);
void command(int, char **) override;
};
} // namespace LAMMPS_NS
#endif
#endif

104
src/reset_mol_ids.cpp → src/reset_atoms_mol.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -16,7 +15,7 @@
Contributing author: Jacob Gissinger (jacob.gissinger@colorado.edu)
------------------------------------------------------------------------- */
#include "reset_mol_ids.h"
#include "reset_atoms_mol.h"
#include "atom.h"
#include "comm.h"
@ -33,10 +32,8 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ResetMolIDs::ResetMolIDs(LAMMPS *lmp) : Command(lmp) {
cfa = nullptr;
cca = nullptr;
ResetAtomsMol::ResetAtomsMol(LAMMPS *lmp) : Command(lmp), cfa(nullptr), cca(nullptr)
{
// default settings
compressflag = 1;
@ -49,49 +46,49 @@ ResetMolIDs::ResetMolIDs(LAMMPS *lmp) : Command(lmp) {
/* ---------------------------------------------------------------------- */
ResetMolIDs::~ResetMolIDs()
ResetAtomsMol::~ResetAtomsMol()
{
if (!idfrag.empty()) modify->delete_compute(idfrag);
if (compressflag && !idchunk.empty()) modify->delete_compute(idchunk);
}
/* ----------------------------------------------------------------------
called as reset_mol_ids command in input script
called as reset_atoms mol command in input script
------------------------------------------------------------------------- */
void ResetMolIDs::command(int narg, char **arg)
void ResetAtomsMol::command(int narg, char **arg)
{
if (domain->box_exist == 0)
error->all(FLERR,"Reset_mol_ids command before simulation box is defined");
if (atom->tag_enable == 0)
error->all(FLERR,"Cannot use reset_mol_ids unless atoms have IDs");
error->all(FLERR, "Reset_atoms mol command before simulation box is defined");
if (atom->tag_enable == 0) error->all(FLERR, "Cannot use reset_atoms mol unless atoms have IDs");
if (atom->molecular != Atom::MOLECULAR)
error->all(FLERR,"Can only use reset_mol_ids on molecular systems");
error->all(FLERR, "Can only use reset_atoms mol on molecular systems");
// process args
if (narg < 1) error->all(FLERR,"Illegal reset_mol_ids command");
if (narg < 1) utils::missing_cmd_args(FLERR, "reset_atoms mol", error);
char *groupid = arg[0];
int iarg = 1;
while (iarg < narg) {
if (strcmp(arg[iarg],"compress") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal reset_mol_ids command");
compressflag = utils::logical(FLERR,arg[iarg+1],false,lmp);
if (strcmp(arg[iarg], "compress") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "reset_atoms mol compress", error);
compressflag = utils::logical(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"single") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal reset_mol_ids command");
singleflag = utils::logical(FLERR,arg[iarg+1],false,lmp);
} else if (strcmp(arg[iarg], "single") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "reset_atoms mol single", error);
singleflag = utils::logical(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"offset") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal reset_mol_ids command");
offset = utils::tnumeric(FLERR,arg[iarg+1],true,lmp);
if (offset < -1) error->all(FLERR,"Illegal reset_mol_ids command");
} else if (strcmp(arg[iarg], "offset") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "reset_atoms mol offset", error);
offset = utils::tnumeric(FLERR, arg[iarg + 1], true, lmp);
if (offset < -1) error->all(FLERR, "Illegal reset_atoms mol offset: {}", offset);
iarg += 2;
} else error->all(FLERR,"Illegal reset_mol_ids command");
} else
error->all(FLERR, "Unknown reset_atoms mol keyword: {}", arg[iarg]);
}
if (comm->me == 0) utils::logmesg(lmp,"Resetting molecule IDs ...\n");
if (comm->me == 0) utils::logmesg(lmp, "Resetting molecule IDs ...\n");
// record wall time for resetting molecule IDs
@ -112,11 +109,11 @@ void ResetMolIDs::command(int narg, char **arg)
comm->setup();
comm->exchange();
comm->borders();
if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
if (domain->triclinic) domain->lamda2x(atom->nlocal + atom->nghost);
// create computes
create_computes((char *) "COMMAND",groupid);
create_computes((char *) "COMMAND", groupid);
// reset molecule IDs
@ -128,44 +125,43 @@ void ResetMolIDs::command(int narg, char **arg)
if (comm->me == 0) {
if (nchunk < 0)
utils::logmesg(lmp," number of new molecule IDs = unknown\n");
utils::logmesg(lmp, " number of new molecule IDs = unknown\n");
else
utils::logmesg(lmp," number of new molecule IDs = {}\n",nchunk);
utils::logmesg(lmp," reset_mol_ids CPU = {:.3f} seconds\n",
platform::walltime()-time1);
utils::logmesg(lmp, " number of new molecule IDs = {}\n", nchunk);
utils::logmesg(lmp, " reset_atoms mol CPU = {:.3f} seconds\n", platform::walltime() - time1);
}
}
/* ----------------------------------------------------------------------
create computes used by reset_mol_ids
create computes used by reset_atoms_mol
------------------------------------------------------------------------- */
void ResetMolIDs::create_computes(char *fixid, char *groupid)
void ResetAtomsMol::create_computes(char *fixid, char *groupid)
{
int igroup = group->find(groupid);
if (igroup == -1) error->all(FLERR,"Could not find reset_mol_ids group ID");
if (igroup < 0) error->all(FLERR, "Could not find reset_atoms mol group {}", groupid);
groupbit = group->bitmask[igroup];
// create instances of compute fragment/atom, compute reduce (if needed),
// and compute chunk/atom. all use the group-ID for this command.
// 'fixid' allows for creating independent instances of the computes
idfrag = fmt::format("{}_reset_mol_ids_FRAGMENT_ATOM",fixid);
idfrag = fmt::format("{}_reset_atoms_mol_FRAGMENT_ATOM", fixid);
auto use_single = singleflag ? "yes" : "no";
cfa = dynamic_cast<ComputeFragmentAtom *>(
modify->add_compute(fmt::format("{} {} fragment/atom single {}",idfrag,groupid,use_single)));
cfa = dynamic_cast<ComputeFragmentAtom *>(modify->add_compute(
fmt::format("{} {} fragment/atom single {}", idfrag, groupid, use_single)));
idchunk = fmt::format("{}_reset_mol_ids_CHUNK_ATOM",fixid);
idchunk = fmt::format("{}_reset_atoms_mol_CHUNK_ATOM", fixid);
if (compressflag)
cca = dynamic_cast<ComputeChunkAtom *>(
modify->add_compute(fmt::format("{} {} chunk/atom molecule compress yes",idchunk,groupid)));
cca = dynamic_cast<ComputeChunkAtom *>(modify->add_compute(
fmt::format("{} {} chunk/atom molecule compress yes", idchunk, groupid)));
}
/* ----------------------------------------------------------------------
called from command() and directly from fixes that update molecules
------------------------------------------------------------------------- */
void ResetMolIDs::reset()
void ResetAtomsMol::reset()
{
// invoke peratom method of compute fragment/atom
// walks bond connectivity and assigns each atom a fragment ID
@ -181,8 +177,7 @@ void ResetMolIDs::reset()
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
molecule[i] = static_cast<tagint> (fragIDs[i]);
if (mask[i] & groupbit) molecule[i] = static_cast<tagint>(fragIDs[i]);
// if compressflag = 0, done
// set nchunk = -1 since cannot easily determine # of new molecule IDs
@ -208,7 +203,7 @@ void ResetMolIDs::reset()
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
if (fragIDs[i] == 0.0) mysingle = 1;
MPI_Allreduce(&mysingle,&singleexist,1,MPI_INT,MPI_MAX,world);
MPI_Allreduce(&mysingle, &singleexist, 1, MPI_INT, MPI_MAX, world);
if (singleexist) nchunk--;
}
@ -220,10 +215,10 @@ void ResetMolIDs::reset()
if (groupbit != 1) {
tagint mymol = 0;
for (int i = 0; i < nlocal; i++)
if (!(mask[i] & groupbit))
mymol = MAX(mymol,molecule[i]);
MPI_Allreduce(&mymol,&offset,1,MPI_LMP_TAGINT,MPI_MAX,world);
} else offset = 0;
if (!(mask[i] & groupbit)) mymol = MAX(mymol, molecule[i]);
MPI_Allreduce(&mymol, &offset, 1, MPI_LMP_TAGINT, MPI_MAX, world);
} else
offset = 0;
}
// reset molecule ID for all atoms in group
@ -231,11 +226,14 @@ void ResetMolIDs::reset()
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
auto newid = static_cast<tagint>(chunkIDs[i]);
auto newid = static_cast<tagint>(chunkIDs[i]);
if (singleexist) {
if (newid == 1) newid = 0;
else newid += offset - 1;
} else newid += offset;
if (newid == 1)
newid = 0;
else
newid += offset - 1;
} else
newid += offset;
molecule[i] = newid;
}
}

13
src/reset_mol_ids.h → src/reset_atoms_mol.h
View File

@ -13,21 +13,21 @@
#ifdef COMMAND_CLASS
// clang-format off
CommandStyle(reset_mol_ids,ResetMolIDs);
CommandStyle(RESET_ATOMS_MOL,ResetAtomsMol);
// clang-format on
#else
#ifndef LMP_RESET_MOL_IDS_H
#define LMP_RESET_MOL_IDS_H
#ifndef LMP_RESET_ATOMS_MOL_H
#define LMP_RESET_ATOMS_MOL_H
#include "command.h"
namespace LAMMPS_NS {
class ResetMolIDs : public Command {
class ResetAtomsMol : public Command {
public:
ResetMolIDs(class LAMMPS *);
~ResetMolIDs() override;
ResetAtomsMol(class LAMMPS *);
~ResetAtomsMol() override;
void command(int, char **) override;
void create_computes(char *, char *);
void reset();
@ -43,7 +43,6 @@ class ResetMolIDs : public Command {
class ComputeFragmentAtom *cfa;
class ComputeChunkAtom *cca;
};
} // namespace LAMMPS_NS
#endif

8
unittest/commands/CMakeLists.txt
View File

@ -53,10 +53,10 @@ endif()
target_link_libraries(test_kim_commands PRIVATE lammps GTest::GMock)
add_test(NAME KimCommands COMMAND test_kim_commands)
add_executable(test_reset_ids test_reset_ids.cpp)
target_compile_definitions(test_reset_ids PRIVATE -DTEST_INPUT_FOLDER=${CMAKE_CURRENT_SOURCE_DIR})
target_link_libraries(test_reset_ids PRIVATE lammps GTest::GMock)
add_test(NAME ResetIDs COMMAND test_reset_ids)
add_executable(test_reset_atoms test_reset_atoms.cpp)
target_compile_definitions(test_reset_atoms PRIVATE -DTEST_INPUT_FOLDER=${CMAKE_CURRENT_SOURCE_DIR})
target_link_libraries(test_reset_atoms PRIVATE lammps GTest::GMock)
add_test(NAME ResetAtoms COMMAND test_reset_atoms)
if(PKG_MOLECULE)
add_executable(test_compute_global test_compute_global.cpp)

208
unittest/commands/test_reset_ids.cpp → unittest/commands/test_reset_atoms.cpp
View File

@ -17,6 +17,7 @@
#include "info.h"
#include "input.h"
#include "lammps.h"
#include "library.h"
#include "output.h"
#include "update.h"
#include "utils.h"
@ -36,11 +37,11 @@ namespace LAMMPS_NS {
#define STRINGIFY(val) XSTR(val)
#define XSTR(val) #val
class ResetIDsTest : public LAMMPSTest {
class ResetAtomsIDTest : public LAMMPSTest {
protected:
void SetUp() override
{
testbinary = "ResetIDsTest";
testbinary = "ResetAtomsIDTest";
LAMMPSTest::SetUp();
if (info->has_style("atom", "full")) {
BEGIN_HIDE_OUTPUT();
@ -51,7 +52,7 @@ protected:
}
};
TEST_F(ResetIDsTest, MolIDAll)
TEST_F(ResetAtomsIDTest, MolIDAll)
{
if (lmp->atom->natoms == 0) GTEST_SKIP();
@ -89,7 +90,7 @@ TEST_F(ResetIDsTest, MolIDAll)
// the original data file has two different molecule IDs
// for two residues of the same molecule/fragment.
BEGIN_HIDE_OUTPUT();
command("reset_mol_ids all");
command("reset_atoms mol all");
END_HIDE_OUTPUT();
ASSERT_EQ(molid[GETIDX(1)], 1);
@ -123,7 +124,7 @@ TEST_F(ResetIDsTest, MolIDAll)
ASSERT_EQ(molid[GETIDX(29)], 5);
}
TEST_F(ResetIDsTest, DeletePlusAtomID)
TEST_F(ResetAtomsIDTest, DeletePlusAtomID)
{
if (lmp->atom->natoms == 0) GTEST_SKIP();
@ -134,7 +135,7 @@ TEST_F(ResetIDsTest, DeletePlusAtomID)
command("group allwater molecule 3:6");
command("group twowater molecule 4:6:2");
command("delete_atoms group twowater compress no bond yes");
command("reset_mol_ids all");
command("reset_atoms mol all");
END_HIDE_OUTPUT();
ASSERT_EQ(lmp->atom->natoms, 23);
ASSERT_EQ(lmp->atom->map_tag_max, 26);
@ -193,7 +194,7 @@ TEST_F(ResetIDsTest, DeletePlusAtomID)
ASSERT_GE(GETIDX(26), 0);
BEGIN_HIDE_OUTPUT();
command("reset_atom_ids");
command("reset_atoms id");
END_HIDE_OUTPUT();
ASSERT_EQ(lmp->atom->map_tag_max, 23);
@ -201,7 +202,7 @@ TEST_F(ResetIDsTest, DeletePlusAtomID)
ASSERT_GE(GETIDX(i), 0);
}
TEST_F(ResetIDsTest, PartialOffset)
TEST_F(ResetAtomsIDTest, PartialOffset)
{
if (lmp->atom->natoms == 0) GTEST_SKIP();
@ -211,7 +212,7 @@ TEST_F(ResetIDsTest, PartialOffset)
BEGIN_HIDE_OUTPUT();
command("group allwater molecule 3:6");
command("group nowater subtract all allwater");
command("reset_mol_ids allwater offset 4");
command("reset_atoms mol allwater offset 4");
END_HIDE_OUTPUT();
ASSERT_EQ(lmp->atom->natoms, 29);
ASSERT_EQ(lmp->atom->map_tag_max, 29);
@ -247,7 +248,7 @@ TEST_F(ResetIDsTest, PartialOffset)
ASSERT_EQ(molid[GETIDX(29)], 8);
BEGIN_HIDE_OUTPUT();
command("reset_mol_ids nowater offset 0");
command("reset_atoms mol nowater offset 0");
END_HIDE_OUTPUT();
ASSERT_EQ(molid[GETIDX(1)], 1);
@ -281,7 +282,7 @@ TEST_F(ResetIDsTest, PartialOffset)
ASSERT_EQ(molid[GETIDX(29)], 8);
}
TEST_F(ResetIDsTest, DeleteAdd)
TEST_F(ResetAtomsIDTest, DeleteAdd)
{
if (lmp->atom->natoms == 0) GTEST_SKIP();
@ -293,7 +294,7 @@ TEST_F(ResetIDsTest, DeleteAdd)
command("group twowater molecule 4:6:2");
command("group nowater subtract all allwater");
command("delete_atoms group twowater compress no bond yes mol yes");
command("reset_mol_ids allwater");
command("reset_atoms mol allwater");
END_HIDE_OUTPUT();
ASSERT_EQ(lmp->atom->natoms, 23);
ASSERT_EQ(lmp->atom->map_tag_max, 26);
@ -352,7 +353,7 @@ TEST_F(ResetIDsTest, DeleteAdd)
ASSERT_GE(GETIDX(26), 0);
BEGIN_HIDE_OUTPUT();
command("reset_atom_ids sort yes");
command("reset_atoms id sort yes");
END_HIDE_OUTPUT();
ASSERT_EQ(lmp->atom->map_tag_max, 23);
@ -360,7 +361,7 @@ TEST_F(ResetIDsTest, DeleteAdd)
ASSERT_GE(GETIDX(i), 0);
BEGIN_HIDE_OUTPUT();
command("reset_mol_ids nowater offset 1");
command("reset_atoms mol nowater offset 1");
END_HIDE_OUTPUT();
ASSERT_EQ(molid[GETIDX(1)], 2);
@ -392,7 +393,7 @@ TEST_F(ResetIDsTest, DeleteAdd)
command("create_atoms 2 single 1.0 0.0 0.0");
command("create_atoms 3 single 2.0 0.0 0.0");
command("create_atoms 4 single 3.0 0.0 0.0");
command("reset_mol_ids all single yes");
command("reset_atoms mol all single yes");
END_HIDE_OUTPUT();
ASSERT_EQ(lmp->atom->natoms, 27);
ASSERT_EQ(lmp->atom->map_tag_max, 27);
@ -406,7 +407,7 @@ TEST_F(ResetIDsTest, DeleteAdd)
ASSERT_EQ(molid[GETIDX(27)], 7);
BEGIN_HIDE_OUTPUT();
command("reset_mol_ids all single no");
command("reset_atoms mol all single no");
END_HIDE_OUTPUT();
ASSERT_EQ(molid[GETIDX(21)], 3);
@ -418,7 +419,7 @@ TEST_F(ResetIDsTest, DeleteAdd)
ASSERT_EQ(molid[GETIDX(27)], 0);
BEGIN_HIDE_OUTPUT();
command("reset_mol_ids all compress no single yes");
command("reset_atoms mol all compress no single yes");
END_HIDE_OUTPUT();
ASSERT_EQ(molid[GETIDX(21)], 21);
@ -430,7 +431,7 @@ TEST_F(ResetIDsTest, DeleteAdd)
ASSERT_EQ(molid[GETIDX(27)], 27);
}
TEST_F(ResetIDsTest, TopologyData)
TEST_F(ResetAtomsIDTest, TopologyData)
{
if (lmp->atom->natoms == 0) GTEST_SKIP();
@ -525,7 +526,7 @@ TEST_F(ResetIDsTest, TopologyData)
ASSERT_EQ(angle_atom3[GETIDX(24)][0], 26);
BEGIN_HIDE_OUTPUT();
command("reset_atom_ids sort yes");
command("reset_atoms id sort yes");
END_HIDE_OUTPUT();
num_bond = lmp->atom->num_bond;
@ -618,63 +619,186 @@ TEST_F(ResetIDsTest, TopologyData)
ASSERT_EQ(angle_atom3[GETIDX(22)][0], 23);
}
TEST_F(ResetIDsTest, DeathTests)
TEST_F(ResetAtomsIDTest, DeathTests)
{
if (lmp->atom->natoms == 0) GTEST_SKIP();
TEST_FAILURE(".*ERROR: Illegal reset_mol_ids command.*", command("reset_mol_ids"););
TEST_FAILURE(".*ERROR: Illegal reset_mol_ids command.*",
command("reset_mol_ids all offset 1 1"););
TEST_FAILURE(".*ERROR: Illegal reset_mol_ids command.*",
command("reset_mol_ids all offset -2"););
TEST_FAILURE(".*ERROR on proc 0: Expected integer.*", command("reset_mol_ids all offset xxx"););
TEST_FAILURE(".*ERROR: Illegal reset_atoms mol command.*", command("reset_atoms mol"););
TEST_FAILURE(".*ERROR: Unknown reset_atoms mol keyword: 1.*",
command("reset_atoms mol all offset 1 1"););
TEST_FAILURE(".*ERROR: Illegal reset_atoms mol offset: -2.*",
command("reset_atoms mol all offset -2"););
TEST_FAILURE(".*ERROR on proc 0: Expected integer.*",
command("reset_mol_ids all compress yes single no offset xxx"););
TEST_FAILURE(".*ERROR: Illegal reset_mol_ids command.*", command("reset_mol_ids all offset"););
TEST_FAILURE(".*ERROR: Illegal reset_mol_ids command.*",
command("reset_mol_ids all compress"););
command("reset_atoms mol all offset xxx"););
TEST_FAILURE(".*ERROR on proc 0: Expected integer.*",
command("reset_atoms mol all compress yes single no offset xxx"););
TEST_FAILURE(".*ERROR: Illegal reset_atoms mol offset command: missing argument.*",
command("reset_atoms mol all offset"););
TEST_FAILURE(".*ERROR: Illegal reset_atoms mol compress command: missing argument.*",
command("reset_atoms mol all compress"););
TEST_FAILURE(".*ERROR: Expected boolean parameter instead of 'xxx'.*",
command("reset_mol_ids all compress xxx"););
TEST_FAILURE(".*ERROR: Illegal reset_mol_ids command.*", command("reset_mol_ids all single"););
command("reset_atoms mol all compress xxx"););
TEST_FAILURE(".*ERROR: Illegal reset_atoms mol single command: missing argument.*",
command("reset_atoms mol all single"););
TEST_FAILURE(".*ERROR: Expected boolean parameter instead of 'xxx'.*",
command("reset_mol_ids all single xxx"););
command("reset_atoms mol all single xxx"););
TEST_FAILURE(".*ERROR: Illegal reset_atoms image command: missing argument.*",
command("reset_atoms image"););
TEST_FAILURE(".*ERROR: Unknown reset_atoms image keyword: xxx.*",
command("reset_atoms image all xxx"););
TEST_FAILURE(".*ERROR: Could not find reset_atoms image group xxx.*",
command("reset_atoms image xxx"););
}
class ResetMolIDsTest : public LAMMPSTest {
class ResetAtomsMolTest : public LAMMPSTest {
protected:
void SetUp() override
{
testbinary = "ResetIDsTest";
testbinary = "ResetAtomsMolTest";
LAMMPSTest::SetUp();
}
};
TEST_F(ResetMolIDsTest, FailBeforeBox)
TEST_F(ResetAtomsMolTest, FailBeforeBox)
{
TEST_FAILURE(".*ERROR: Reset_mol_ids command before simulation box is.*",
command("reset_mol_ids all"););
TEST_FAILURE(".*ERROR: Reset_atoms id command before simulation box is.*",
command("reset_atoms id"););
TEST_FAILURE(".*ERROR: Reset_atoms mol command before simulation box is.*",
command("reset_atoms mol all"););
TEST_FAILURE(".*ERROR: Reset_atoms image command before simulation box is.*",
command("reset_atoms image all"););
}
TEST_F(ResetMolIDsTest, FailMissingId)
TEST_F(ResetAtomsMolTest, FailMissingId)
{
BEGIN_HIDE_OUTPUT();
command("atom_modify id no");
command("region box block 0 1 0 1 0 1");
command("create_box 1 box");
END_HIDE_OUTPUT();
TEST_FAILURE(".*ERROR: Cannot use reset_mol_ids unless.*", command("reset_mol_ids all"););
TEST_FAILURE(".*ERROR: Cannot use reset_atoms mol unless.*", command("reset_atoms mol all"););
TEST_FAILURE(".*ERROR: Cannot use reset_atoms image unless.*",
command("reset_atoms image all"););
}
TEST_F(ResetMolIDsTest, FailOnlyMolecular)
TEST_F(ResetAtomsMolTest, FailOnlyMolecular)
{
BEGIN_HIDE_OUTPUT();
command("clear");
command("region box block 0 1 0 1 0 1");
command("create_box 1 box");
END_HIDE_OUTPUT();
TEST_FAILURE(".*ERROR: Can only use reset_mol_ids.*", command("reset_mol_ids all"););
TEST_FAILURE(".*ERROR: Can only use reset_atoms mol.*", command("reset_atoms mol all"););
}
class ResetAtomsImageTest : public LAMMPSTest {
protected:
void SetUp() override
{
testbinary = "ResetAtomsImageTest";
LAMMPSTest::SetUp();
if (info->has_style("atom", "full")) {
BEGIN_HIDE_OUTPUT();
command("variable input_dir index \"" STRINGIFY(TEST_INPUT_FOLDER) "\"");
command("include \"${input_dir}/in.fourmol\"");
command("create_atoms 1 single 1.0 1.0 1.0");
command("create_atoms 1 single 2.0 1.0 1.0");
command("create_atoms 1 single 1.0 2.0 1.0");
command("set atom 1*7 image 1 1 -1");
command("set atom 8*9 image 2 -1 0");
command("set atom 8*9 image 2 -1 0");
command("set atom 10 image 2 0 0");
command("set atom 11*12 image 2 0 -1");
command("set atom 13*15 image 1 0 0");
command("set atom 16*17 image 0 1 1");
command("set atom 18*19 image 0 1 0");
command("set atom 20 image 0 2 0");
command("set atom 21*23 image 20 -1 0");
command("set atom 27*28 image 1 0 0");
command("set atom 29 image 1 2 0");
command("set atom 31 image -2 0 1");
command("set atom 32 image 0 20 0");
END_HIDE_OUTPUT();
}
}
};
TEST_F(ResetAtomsImageTest, ResetAtomsImage)
{
if (lmp->atom->natoms == 0) GTEST_SKIP();
EXPECT_EQ(lmp->atom->image[GETIDX(1)], lammps_encode_image_flags(1, 1, -1));
EXPECT_EQ(lmp->atom->image[GETIDX(2)], lammps_encode_image_flags(1, 1, -1));
EXPECT_EQ(lmp->atom->image[GETIDX(3)], lammps_encode_image_flags(1, 1, -1));
EXPECT_EQ(lmp->atom->image[GETIDX(4)], lammps_encode_image_flags(1, 1, -1));
EXPECT_EQ(lmp->atom->image[GETIDX(5)], lammps_encode_image_flags(1, 1, -1));
EXPECT_EQ(lmp->atom->image[GETIDX(6)], lammps_encode_image_flags(1, 1, -1));
EXPECT_EQ(lmp->atom->image[GETIDX(7)], lammps_encode_image_flags(1, 1, -1));
EXPECT_EQ(lmp->atom->image[GETIDX(8)], lammps_encode_image_flags(2, -1, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(9)], lammps_encode_image_flags(2, -1, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(10)], lammps_encode_image_flags(2, 0, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(11)], lammps_encode_image_flags(2, 0, -1));
EXPECT_EQ(lmp->atom->image[GETIDX(12)], lammps_encode_image_flags(2, 0, -1));
EXPECT_EQ(lmp->atom->image[GETIDX(13)], lammps_encode_image_flags(1, 0, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(14)], lammps_encode_image_flags(1, 0, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(15)], lammps_encode_image_flags(1, 0, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(16)], lammps_encode_image_flags(0, 1, 1));
EXPECT_EQ(lmp->atom->image[GETIDX(17)], lammps_encode_image_flags(0, 1, 1));
EXPECT_EQ(lmp->atom->image[GETIDX(18)], lammps_encode_image_flags(0, 1, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(19)], lammps_encode_image_flags(0, 1, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(20)], lammps_encode_image_flags(0, 2, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(21)], lammps_encode_image_flags(20, -1, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(22)], lammps_encode_image_flags(20, -1, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(23)], lammps_encode_image_flags(20, -1, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(24)], lammps_encode_image_flags(0, 0, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(25)], lammps_encode_image_flags(0, 0, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(26)], lammps_encode_image_flags(0, 0, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(27)], lammps_encode_image_flags(1, 0, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(28)], lammps_encode_image_flags(1, 0, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(29)], lammps_encode_image_flags(1, 2, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(30)], lammps_encode_image_flags(0, 0, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(31)], lammps_encode_image_flags(-2, 0, 1));
EXPECT_EQ(lmp->atom->image[GETIDX(32)], lammps_encode_image_flags(0, 20, 0));
BEGIN_HIDE_OUTPUT();
command("group subset id 5:32");
command("reset_atoms image subset");
END_HIDE_OUTPUT();
EXPECT_EQ(lmp->atom->image[GETIDX(1)], lammps_encode_image_flags(1, 1, -1));
EXPECT_EQ(lmp->atom->image[GETIDX(2)], lammps_encode_image_flags(1, 1, -1));
EXPECT_EQ(lmp->atom->image[GETIDX(3)], lammps_encode_image_flags(1, 1, -1));
EXPECT_EQ(lmp->atom->image[GETIDX(4)], lammps_encode_image_flags(1, 1, -1));
EXPECT_EQ(lmp->atom->image[GETIDX(5)], lammps_encode_image_flags(0, 1, -1));
EXPECT_EQ(lmp->atom->image[GETIDX(6)], lammps_encode_image_flags(0, 1, -1));
EXPECT_EQ(lmp->atom->image[GETIDX(7)], lammps_encode_image_flags(0, 1, -1));
EXPECT_EQ(lmp->atom->image[GETIDX(8)], lammps_encode_image_flags(1, -1, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(9)], lammps_encode_image_flags(1, -1, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(10)], lammps_encode_image_flags(1, 0, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(11)], lammps_encode_image_flags(1, 0, -1));
EXPECT_EQ(lmp->atom->image[GETIDX(12)], lammps_encode_image_flags(1, 0, -1));
EXPECT_EQ(lmp->atom->image[GETIDX(13)], lammps_encode_image_flags(0, 0, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(14)], lammps_encode_image_flags(0, 0, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(15)], lammps_encode_image_flags(0, 0, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(16)], lammps_encode_image_flags(-1, 1, 1));
EXPECT_EQ(lmp->atom->image[GETIDX(17)], lammps_encode_image_flags(-1, 1, 1));
EXPECT_EQ(lmp->atom->image[GETIDX(18)], lammps_encode_image_flags(0, 0, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(19)], lammps_encode_image_flags(0, 0, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(20)], lammps_encode_image_flags(0, 1, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(21)], lammps_encode_image_flags(0, 0, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(22)], lammps_encode_image_flags(0, 0, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(23)], lammps_encode_image_flags(0, 0, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(24)], lammps_encode_image_flags(0, 0, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(25)], lammps_encode_image_flags(0, 0, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(26)], lammps_encode_image_flags(0, 0, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(27)], lammps_encode_image_flags(0, -1, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(28)], lammps_encode_image_flags(0, -1, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(29)], lammps_encode_image_flags(0, 1, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(30)], lammps_encode_image_flags(0, 0, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(30)], lammps_encode_image_flags(0, 0, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(31)], lammps_encode_image_flags(0, 0, 0));
EXPECT_EQ(lmp->atom->image[GETIDX(32)], lammps_encode_image_flags(0, 0, 0));
}
} // namespace LAMMPS_NS
int main(int argc, char **argv)