new compute: reaxff/atom

[rbberger]

Summary

This new compute allows extracting the local bond information from ReaxFF.

Related Issue(s)

Author(s)

@rbberger

Licensing

By submitting this pull request, I agree, that my contribution will be included in LAMMPS and redistributed under either the GNU General Public License version 2 (GPL v2) or the GNU Lesser General Public License version 2.1 (LGPL v2.1).

Backward Compatibility

Implementation Notes

Post Submission Checklist

• The feature or features in this pull request is complete
• Licensing information is complete
• Corresponding author information is complete
• The source code follows the LAMMPS formatting guidelines
• Suitable new documentation files and/or updates to the existing docs are included
• The added/updated documentation is integrated and tested with the documentation build system
• The feature has been verified to work with the conventional build system
• The feature has been verified to work with the CMake based build system
• Suitable tests have been added to the unittest tree.
• A package specific README file has been included or updated
• One or more example input decks are included

Further Information, Files, and Links

[sjplimp]

@rbberger Thanks for working on this. A few comments ...

(1) This appears to be a per-atom compute, not a local compute ?
A local compute would output one line of data per bond, not one line per atom.
E.g. see compute bond/local for an example.

(2) I don't see any logic using the group arg of the compute
This could be chosen by the user to generate info on bonds for a subset of the system,
e.g. a subset of the atoms or a subset of bonds (pairs of atoms).
Again, see compute bond/local for an example.

(3) If you stick with a per-atom compute, it appears you are making size_local_rows
a dynamic quantity (during a simulation). I think this may cause problems with
users of the compute output (like a dump custom command) which check their indexing
args before a run starts, and will expect to be able to access MAXREAXBOND columns
from compute->array_atom (or array_local as written now). Which is probably another reason to
make this a local compute, not per-atom.

[sjplimp]

@rbberger I should also note that it is now possible for a single compute to generate both per-atom and local data. So if it is important to be able dump per-atom values (total_bond_order, nlp, numbonds) along with per-bond info (i,j,abo) then this
compute could do both. The user could access the data via 2 dump files, custom and local. There could be a flag in the compute for the user to specify which data they want (either or both). See compute voronoi/atom for an example. If this is what you want to do, then the compute should be named something like compute reaxff/atom or compute reaxff/bond/atom.

[akohlmey]

@rbberger I just merged PR #3896 which standardizes access to multiple types of computed data in fixes and computes from variables and similar. This is implementing some of the things that Steve mentioned but that were not available to the standard code at the time. So it may be advantageous to merge with upstream and check your code.

[rbberger]

@akohlmey thanks for the heads up. I was just running into this while trying to create the "Steve" version.

[rbberger]

@akohlmey if you look at the modified example you'll see I needed to hard code the columns for dump local, since the wildcard expansion prefers the peratom columns. It might be beneficial to allow the expansion function to set a preference.

[rbberger]

@sjplimp updated the PR to now do both peratom and local arrays in the compute. You were right that it was breaking dump local. I'm still missing the group logic, but since that was also missing from the fix, I don't see it as a priority. I want to first nail down this base version.

I've also created a Kokkos version that reuses the existing PackBuffer logic and then unpacks the buffer onto host arrays. @stanmoore1 can you take a look?

[stanmoore1]

Hi Richard, this generally looks really good. I have a few minor points in the comments.

[stanmoore1]

Pair hybrid will break this logic, if there is more than 1 instance of pair reaxff

[rbberger]

Please see if fd83ed4 addresses your concerns.

[stanmoore1]

I think this may break if using the suffix command--it looks like force->pair_match doesn't automatically append the suffix? You can check for reax instead, then use the kokkosable flag which is set by every pair style.

[rbberger]

Please see if fd83ed4 addresses your concerns.

[sjplimp]

@rbberger - thanks for the changes - a few more comments:

(1) As mentioned earlier, for consistency with compute voronoi/atom, the name should probably be compute reaxff/atom or compute reaxff/bonds/atom

(2) Similar to voronoi/atom, we should give users the option to (not) produce the (large) list of bonds. Voronoi/atom does this
with a "neighbors no/yes" keyword (it always produces the peratom info). Maybe a "bonds no/yes" keyword for this compute
with no = default.

(3) I think the per-atom output should only be 3 values: total_bond_order, nlp, numneigh. The rest of the values can be accessed with standard dump custom keywords. Again, this is similar to compute voronoi/atom, with numneigh last.

(4) @akohlmey I agree with Richard that utils::expand_args() should work
when called from dump_custom.cpp or dump_local.cpp to expand a per-atom or local array. Perhaps this
could be encoded somehow in the "mode" arg ? Or a new optional arg ?
Note that compute reduce also has a new option to operate on
peratom or local data and thus its expand_args() call should match what dump custom versus dump local does
to expand either case. This is the logic in utils::expand_args() that needs to be enriched and ditto for the fix section:

// check for global vector/array, peratom array, local array

    if (compute) {
      if (mode == 0 && compute->vector_flag) {
        nmax = compute->size_vector;
        expandflag = 1;
      } else if (mode == 1 && compute->array_flag) {
        nmax = compute->size_array_cols;
        expandflag = 1;
      } else if (compute->peratom_flag && compute->size_peratom_cols) {
        nmax = compute->size_peratom_cols;
        expandflag = 1;
      } else if (compute->local_flag && compute->size_local_cols) {
        nmax = compute->size_local_cols;
        expandflag = 1;
      }
    }

[sjplimp]

@rbberger Also, the use of numneigh[i] in the code was a bit confusing to me since that is neighbor list nomenclature in LAMMPS.
In your case it is the number of bonds for each atom (bondorder > threshold). So I suggest bondcount[i] or something similar.

[rbberger]

@sjplimp that's just duplicating existing reaxff code (see fix reaxff/bonds). But of course, we can definitely give it better names.

[rbberger]

@sjplimp implemented (1), (2) and (3) from your last comment. Still missing @stanmoore1's comments.

[sjplimp]

Thanks @rbberger - changes look good. Just this small item:
(1) In compute_local(), I think it should be bond[0] = tag[i], not i ? So that it is an atom ID, like bond[1] ?

Also, I'm still hoping you can add group logic (item (2) in my Oct 12 comment), before this is merged.
Should be pretty simple I think.

[rbberger]

@sjplimp (1) makes sense, I actually got a question of how to transform it to tag[i], so it makes more sense to directly return it here. I'll work on group logic and documentation next week.

[stanmoore1]

Looks good to me.

[stanmoore1]

I ran this through Kokkos regression testing, everything passed (though the modified example is not included, so the new code was not run).

[stanmoore1]

Python unit tests failed, seems unrelated to this PR:

The following tests FAILED:
	 60 - PythonNumpy (Failed)
	 62 - PythonPyLammps (Failed)

[akohlmey]

Python unit tests failed, seems unrelated to this PR:

The following tests FAILED:
	 60 - PythonNumpy (Failed)
	 62 - PythonPyLammps (Failed)

I suspect it may be related to PR #4011 in upstream, but the failures are somewhat random. I only see them on the Azure macOS instances that have python 3.12.