Merge pull request #3871 from megmcca/local-composition

Per-atom compute to calculate local composition of atom types

doc/src/Commands_compute.rst

@@ -91,6 +91,7 @@ KOKKOS, o = OPENMP, t = OPT.
    * :doc:`ke/atom/eff <compute_ke_atom_eff>`
    * :doc:`ke/eff <compute_ke_eff>`
    * :doc:`ke/rigid <compute_ke_rigid>`
+   * :doc:`local/comp/atom (k) <compute_local_comp_atom>`
    * :doc:`mliap <compute_mliap>`
    * :doc:`momentum <compute_momentum>`
    * :doc:`msd <compute_msd>`

doc/src/compute.rst

@@ -245,6 +245,7 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` pag
 * :doc:`ke/atom/eff <compute_ke_atom_eff>` - per-atom translational and radial kinetic energy in the electron force field model
 * :doc:`ke/eff <compute_ke_eff>` - kinetic energy of a group of nuclei and electrons in the electron force field model
 * :doc:`ke/rigid <compute_ke_rigid>` - translational kinetic energy of rigid bodies
+* :doc:`local/comp/atom <compute_local_comp_atom>` - local composition for each atom
 * :doc:`mliap <compute_mliap>` - gradients of energy and forces with respect to model parameters and related quantities for training machine learning interatomic potentials
 * :doc:`momentum <compute_momentum>` - translational momentum
 * :doc:`msd <compute_msd>` - mean-squared displacement of group of atoms

doc/src/compute_local_comp_atom.rst

@@ -0,0 +1,118 @@
+.. index:: compute local/comp/atom
+.. index:: compute local/comp/atom/kk
+
+compute local/comp/atom command
+===============================
+
+Accelerator Variants: *local/comp/atom/kk*
+
+Syntax
+""""""
+
+.. code-block:: LAMMPS
+
+   compute ID group-ID local/comp/atom keyword values ...
+
+* ID, group-ID are documented in :doc:`compute <compute>` command
+* local/comp/atom = style name of this compute command
+* one or more keyword/value pairs may be appended
+
+  .. parsed-literal::
+
+     keyword = *cutoff*
+       *cutoff* value = distance cutoff
+
+Examples
+""""""""
+
+.. code-block:: LAMMPS
+
+   compute 1 all local/comp/atom
+
+   compute 1 all local/comp/atom cutoff 9.0
+   comm_modify cutoff 9.0
+
+
+Description
+"""""""""""
+
+.. versionadded:: TBD
+
+Define a computation that calculates a local composition vector for each
+atom. For a central atom with :math:`M` neighbors within the neighbor cutoff sphere,
+composition is defined as the number of atoms of a given type
+(including the central atom) divided by (:math:`M+1`).  For a given central atom,
+the sum of all compositions equals one.
+
+.. note::
+
+   This compute uses the number of atom types, not chemical species, assigned in
+   :doc:`pair_coeff <pair_coeff>` command.  If an interatomic potential has two
+   species (i.e., Cu and Ni) assigned to four different atom types in
+   :doc:`pair_coeff <pair_coeff>` (i.e., 'Cu Cu Ni Ni'), the compute will
+   output four fractional values.  In those cases, the user may desire an extra
+   calculation step to consolidate per-type fractions into per-species fractions.
+   This calculation can be conducted within LAMMPS using another compute such as
+   :doc:`compute reduce <compute_reduce>`, an atom-style :doc:`variable`, or as a
+   post-processing step.
+
+----------
+
+The optional keyword *cutoff* defines the distance cutoff used when
+searching for neighbors. The default value is the cutoff specified by
+the pair style. If no pair style is defined, then a cutoff must be
+defined using this keyword. If the specified cutoff is larger than
+that of the pair_style plus neighbor skin (or no pair style is
+defined), the *comm_modify cutoff* option must also be set to match
+that of the *cutoff* keyword.
+
+The neighbor list needed to compute this quantity is constructed each
+time the calculation is performed (i.e. each time a snapshot of atoms
+is dumped).  Thus it can be inefficient to compute/dump this quantity
+too frequently.
+
+.. note::
+
+   If you have a bonded system, then the settings of
+   :doc:`special_bonds <special_bonds>` command can remove pairwise
+   interactions between atoms in the same bond, angle, or dihedral.
+   This is the default setting for the :doc:`special_bonds
+   <special_bonds>` command, and means those pairwise interactions do
+   not appear in the neighbor list.  Because this compute uses the
+   neighbor list, it also means those pairs will not be included in
+   the order parameter.  This difficulty can be circumvented by
+   writing a dump file, and using the :doc:`rerun <rerun>` command to
+   compute the order parameter for snapshots in the dump file.  The
+   rerun script can use a :doc:`special_bonds <special_bonds>` command
+   that includes all pairs in the neighbor list.
+
+----------
+
+Output info
+"""""""""""
+
+This compute calculates a per-atom array with :math:`1 + N` columns, where :math:`N`
+is the number of atom types. The first column is a count of the number of atoms
+used to calculate composition (including the central atom), and each subsequent
+column indicates the fraction of that atom type within the cutoff sphere.
+
+These values can be accessed by any command that uses per-atom values
+from a compute as input.  See the :doc:`Howto output <Howto_output>`
+doc page for an overview of LAMMPS output options.
+
+Restrictions
+""""""""""""
+
+This compute is part of the EXTRA-COMPUTE package.  It is only enabled
+if LAMMPS was built with that package.  See the :doc:`Build package
+<Build_package>` page for more info.
+
+Related commands
+""""""""""""""""
+
+:doc:`comm_modify <comm_modify>`
+
+Default
+"""""""
+
+The option defaults are *cutoff* = pair style cutoff.

doc/utils/sphinx-config/false_positives.txt

@@ -3621,6 +3621,7 @@ Tk
 Tkin
 tloop
 tlsph
+tm
 tmax
 Tmax
 tmd

src/.gitignore

@@ -580,6 +580,8 @@
 /compute_ke_eff.h
 /compute_ke_rigid.cpp
 /compute_ke_rigid.h
+/compute_local_comp_atom.cpp
+/compute_local_comp_atom.h
 /compute_meso_e_atom.cpp
 /compute_meso_e_atom.h
 /compute_meso_rho_atom.cpp

src/EXTRA-COMPUTE/compute_local_comp_atom.cpp

@@ -0,0 +1,214 @@
+/* ----------------------------------------------------------------------
+   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.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing author: Megan McCarthy (SNL)
+------------------------------------------------------------------------- */
+
+#include "compute_local_comp_atom.h"
+
+#include "atom.h"
+#include "comm.h"
+#include "domain.h"
+#include "error.h"
+#include "force.h"
+#include "math_const.h"
+#include "memory.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "neighbor.h"
+#include "pair.h"
+#include "update.h"
+
+#include <cstring>
+
+using namespace LAMMPS_NS;
+using namespace MathConst;
+
+/* ---------------------------------------------------------------------- */
+
+ComputeLocalCompAtom::ComputeLocalCompAtom(LAMMPS *lmp, int narg, char **arg) :
+    Compute(lmp, narg, arg), result(nullptr)
+{
+  if (narg < 3 || narg > 5) error->all(FLERR, "Illegal compute local/comp/atom command");
+
+  cutoff = 0.0;
+
+  int iarg = 3;
+  while (iarg < narg) {
+    if (strcmp(arg[iarg], "cutoff") == 0) {
+      if (iarg + 2 > narg) error->all(FLERR, "Illegal compute local/comp/atom command");
+      cutoff = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
+      if (cutoff <= 0.0) error->all(FLERR, "Illegal compute local/comp/atom command");
+      iarg += 2;
+    } else
+      error->all(FLERR, "Illegal compute local/comp/atom command");
+  }
+
+  peratom_flag = 1;
+
+  ntypes = atom->ntypes;
+  size_peratom_cols = 1 + ntypes;
+
+  nmax = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeLocalCompAtom::~ComputeLocalCompAtom()
+{
+  if (copymode) return;
+
+  memory->destroy(result);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeLocalCompAtom::init()
+{
+  if (!force->pair && cutoff == 0.0)
+    error->all(FLERR,
+               "Compute local/comp/atom requires a cutoff be specified "
+               "or a pair style be defined");
+
+  double skin = neighbor->skin;
+  if (cutoff != 0.0) {
+    double cutghost;    // as computed by Neighbor and Comm
+    if (force->pair)
+      cutghost = MAX(force->pair->cutforce + skin, comm->cutghostuser);
+    else
+      cutghost = comm->cutghostuser;
+
+    if (cutoff > cutghost)
+      error->all(FLERR,
+                 "Compute local/comp/atom cutoff exceeds ghost atom range - "
+                 "use comm_modify cutoff command");
+  }
+
+  int cutflag = 1;
+  if (force->pair) {
+    if (cutoff == 0.0) { cutoff = force->pair->cutforce; }
+    if (cutoff <= force->pair->cutforce + skin) cutflag = 0;
+  }
+
+  cutsq = cutoff * cutoff;
+
+  // need an occasional full neighbor list
+
+  auto req = neighbor->add_request(this, NeighConst::REQ_FULL | NeighConst::REQ_OCCASIONAL);
+  if (cutflag) req->set_cutoff(cutoff);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeLocalCompAtom::init_list(int /*id*/, NeighList *ptr)
+{
+  list = ptr;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeLocalCompAtom::compute_peratom()
+{
+  int i, j, ii, jj, inum, jnum;
+  double xtmp, ytmp, ztmp, delx, dely, delz, rsq;
+  int *ilist, *jlist, *numneigh, **firstneigh;
+  int count, itype, jtype;
+
+  invoked_peratom = update->ntimestep;
+
+  // grow result array if necessary
+
+  if (atom->nmax > nmax) {
+    memory->destroy(result);
+    nmax = atom->nmax;
+    memory->create(result, nmax, size_peratom_cols, "local/comp/atom:result");
+    array_atom = result;
+  }
+
+  // invoke full neighbor list (will copy or build if necessary)
+
+  neighbor->build_one(list);
+
+  inum = list->inum;
+  ilist = list->ilist;
+  numneigh = list->numneigh;
+  firstneigh = list->firstneigh;
+
+  // compute properties for each atom in group
+  // use full neighbor list to count atoms less than cutoff
+
+  double **x = atom->x;
+  int *type = atom->type;
+  int *mask = atom->mask;
+
+  // get per-atom local compositions
+
+  for (ii = 0; ii < inum; ii++) {
+
+    i = ilist[ii];
+
+    if (mask[i] & groupbit) {
+
+      xtmp = x[i][0];
+      ytmp = x[i][1];
+      ztmp = x[i][2];
+      jlist = firstneigh[i];
+      jnum = numneigh[i];
+
+      // i atom contribution
+
+      count = 1;
+
+      itype = type[i];
+      result[i][itype]++;
+
+      for (jj = 0; jj < jnum; jj++) {
+        j = jlist[jj];
+        j &= NEIGHMASK;
+
+        jtype = type[j];
+
+        delx = xtmp - x[j][0];
+        dely = ytmp - x[j][1];
+        delz = ztmp - x[j][2];
+        rsq = delx * delx + dely * dely + delz * delz;
+        if (rsq < cutsq) {
+          count++;
+          result[i][jtype]++;
+        }
+      }
+
+      // total count of atoms found in sampled radius range
+
+      result[i][0] = count;
+
+      // local comp fractions per element
+
+      double lfac = 1.0 / count;
+      for (int n = 1;  n < size_peratom_cols; n++)
+        result[i][n+1] *= lfac;
+
+    }
+  }
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of local atom-based array
+------------------------------------------------------------------------- */
+
+double ComputeLocalCompAtom::memory_usage()
+{
+  double bytes = (double) 2 * nmax * sizeof(double);
+  return bytes;
+}

src/EXTRA-COMPUTE/compute_local_comp_atom.h

@@ -0,0 +1,47 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   LAMMPS development team: developers@lammps.org
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef COMPUTE_CLASS
+// clang-format off
+ComputeStyle(local/comp/atom,ComputeLocalCompAtom);
+// clang-format on
+#else
+
+#ifndef LMP_COMPUTE_LOCAL_COMP_ATOM_H
+#define LMP_COMPUTE_LOCAL_COMP_ATOM_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeLocalCompAtom : public Compute {
+ public:
+  ComputeLocalCompAtom(class LAMMPS *, int, char **);
+  ~ComputeLocalCompAtom() override;
+  void init() override;
+  void init_list(int, class NeighList *) override;
+  void compute_peratom() override;
+  double memory_usage() override;
+
+ protected:
+  int nmax, ntypes;
+  double cutoff;            // global cutoff distance
+  double cutsq;             // cutoff**2
+  class NeighList *list;    // neighbor list
+  double **result;          // peratom array of local compositions
+};
+
+}    // namespace LAMMPS_NS
+
+#endif
+#endif

src/KOKKOS/Install.sh

@@ -96,6 +96,8 @@ action compute_coord_atom_kokkos.cpp
 action compute_coord_atom_kokkos.h
 action compute_erotate_sphere_kokkos.cpp
 action compute_erotate_sphere_kokkos.h
+action compute_local_comp_atom_kokkos.cpp compute_local_comp_atom.cpp
+action compute_local_comp_atom_kokkos.h compute_local_comp_atom.h
 action compute_orientorder_atom_kokkos.cpp
 action compute_orientorder_atom_kokkos.h
 action compute_temp_deform_kokkos.cpp

src/KOKKOS/compute_local_comp_atom_kokkos.cpp

@@ -0,0 +1,180 @@
+/* ----------------------------------------------------------------------
+   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 Gdirectoryeneral Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing author: Megan McCarthy (SNL)
+------------------------------------------------------------------------- */
+
+#include "compute_local_comp_atom_kokkos.h"
+
+#include "atom_kokkos.h"
+#include "atom_masks.h"
+#include "comm.h"
+#include "domain.h"
+#include "error.h"
+#include "force.h"
+#include "memory_kokkos.h"
+#include "modify.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "neighbor_kokkos.h"
+#include "pair.h"
+#include "update.h"
+
+#include <cmath>
+#include <cstring>
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+template<class DeviceType>
+ComputeLocalCompAtomKokkos<DeviceType>::ComputeLocalCompAtomKokkos(LAMMPS *lmp, int narg, char **arg) :
+  ComputeLocalCompAtom(lmp, narg, arg)
+{
+  kokkosable = 1;
+  atomKK = (AtomKokkos *) atom;
+  execution_space = ExecutionSpaceFromDevice<DeviceType>::space;
+  datamask_read = EMPTY_MASK;
+  datamask_modify = EMPTY_MASK;
+}
+
+/* ---------------------------------------------------------------------- */
+
+template<class DeviceType>
+ComputeLocalCompAtomKokkos<DeviceType>::~ComputeLocalCompAtomKokkos()
+{
+  if (copymode) return;
+
+  memoryKK->destroy_kokkos(k_result,result);
+}
+
+/* ---------------------------------------------------------------------- */
+
+template<class DeviceType>
+void ComputeLocalCompAtomKokkos<DeviceType>::init()
+{
+  ComputeLocalCompAtom::init();
+
+  // adjust neighbor list request for KOKKOS
+
+  auto request = neighbor->find_request(this);
+  request->set_kokkos_host(std::is_same<DeviceType,LMPHostType>::value &&
+                           !std::is_same<DeviceType,LMPDeviceType>::value);
+  request->set_kokkos_device(std::is_same<DeviceType,LMPDeviceType>::value);
+}
+
+/* ---------------------------------------------------------------------- */
+
+template<class DeviceType>
+void ComputeLocalCompAtomKokkos<DeviceType>::compute_peratom()
+{
+  invoked_peratom = update->ntimestep;
+
+  // grow result array if necessary
+
+  int size_peratom_cols = 1 + atom->ntypes;
+  if (atom->nmax > nmax) {
+    memoryKK->destroy_kokkos(k_result,result);
+    nmax = atom->nmax;
+    memoryKK->create_kokkos(k_result,result,nmax,size_peratom_cols,"local/comp/atom:result");
+    d_result = k_result.view<DeviceType>();
+    array_atom = result;
+  }
+
+  // invoke full neighbor list (will copy or build if necessary)
+
+  neighbor->build_one(list);
+  int inum = list->inum;
+
+  NeighListKokkos<DeviceType>* k_list = static_cast<NeighListKokkos<DeviceType>*>(list);
+  d_numneigh = k_list->d_numneigh;
+  d_neighbors = k_list->d_neighbors;
+  d_ilist = k_list->d_ilist;
+
+  // compute properties for each atom in group
+  // use full neighbor list to count atoms less than cutoff
+
+  atomKK->sync(execution_space,X_MASK|TYPE_MASK|MASK_MASK);
+  x = atomKK->k_x.view<DeviceType>();
+  type = atomKK->k_type.view<DeviceType>();
+  mask = atomKK->k_mask.view<DeviceType>();
+  ntypes = atom->ntypes;
+
+  Kokkos::deep_copy(d_result,0.0);
+
+  copymode = 1;
+  typename Kokkos::RangePolicy<DeviceType, TagComputeLocalCompAtom> policy(0,inum);
+  Kokkos::parallel_for("ComputeLocalComp",policy,*this);
+  copymode = 0;
+
+  k_result.modify<DeviceType>();
+  k_result.sync_host();
+}
+
+template<class DeviceType>
+KOKKOS_INLINE_FUNCTION
+void ComputeLocalCompAtomKokkos<DeviceType>::operator()(TagComputeLocalCompAtom, const int &ii) const
+{
+  const int i = d_ilist[ii];
+
+  if (mask[i] & groupbit) {
+
+    const X_FLOAT xtmp = x(i,0);
+    const X_FLOAT ytmp = x(i,1);
+    const X_FLOAT ztmp = x(i,2);
+    const int jnum = d_numneigh[i];
+
+    // i atom contribution
+
+    int count = 1.0;
+    int itype = type[i];
+    d_result(i,itype)++;
+
+    for (int jj = 0; jj < jnum; jj++) {
+
+      int j = d_neighbors(i,jj);
+      j &= NEIGHMASK;
+
+      int jtype = type[j];
+
+      const F_FLOAT delx = x(j,0) - xtmp;
+      const F_FLOAT dely = x(j,1) - ytmp;
+      const F_FLOAT delz = x(j,2) - ztmp;
+      const F_FLOAT rsq = delx*delx + dely*dely + delz*delz;
+      if (rsq < cutsq) {
+        count++;
+        d_result(i,jtype) += 1.0;
+      }
+    }
+
+    // total count of atoms found in sampled radius range
+
+    d_result(i,0) = count;
+
+    // local comp fractions per atom type
+
+    double lfac = 1.0 / count;
+
+    for (int n = 1; n < size_peratom_cols; n++) {
+      d_result(i,n) *= lfac;
+    }
+  }
+}
+
+namespace LAMMPS_NS {
+template class ComputeLocalCompAtomKokkos<LMPDeviceType>;
+#ifdef LMP_KOKKOS_GPU
+template class ComputeLocalCompAtomKokkos<LMPHostType>;
+#endif
+}

src/KOKKOS/compute_local_comp_atom_kokkos.h

@@ -0,0 +1,64 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   LAMMPS development team: developers@lammps.org
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef COMPUTE_CLASS
+// clang-format off
+ComputeStyle(local/comp/atom/kk,ComputeLocalCompAtomKokkos<LMPDeviceType>);
+ComputeStyle(local/comp/atom/kk/device,ComputeLocalCompAtomKokkos<LMPDeviceType>);
+ComputeStyle(local/comp/atom/kk/host,ComputeLocalCompAtomKokkos<LMPHostType>);
+// clang-format on
+
+#else
+
+#ifndef LMP_COMPUTE_LOCAL_COMP_ATOM_KOKKOS_H
+#define LMP_COMPUTE_LOCAL_COMP_ATOM_KOKKOS_H
+
+#include "compute_local_comp_atom.h"
+#include "kokkos_type.h"
+
+namespace LAMMPS_NS {
+
+// clang-format off
+struct TagComputeLocalCompAtom {};
+// clang-format on
+
+template <class DeviceType> class ComputeLocalCompAtomKokkos : public ComputeLocalCompAtom {
+ public:
+  typedef DeviceType device_type;
+  typedef ArrayTypes<DeviceType> AT;
+
+  ComputeLocalCompAtomKokkos(class LAMMPS *, int, char **);
+  ~ComputeLocalCompAtomKokkos() override;
+  void init() override;
+  void compute_peratom() override;
+
+  KOKKOS_INLINE_FUNCTION
+  void operator()(TagComputeLocalCompAtom, const int &) const;
+
+ private:
+
+  typename AT::t_x_array x;
+  typename ArrayTypes<DeviceType>::t_int_1d type;
+  typename ArrayTypes<DeviceType>::t_int_1d mask;
+
+  typename AT::t_neighbors_2d d_neighbors;
+  typename AT::t_int_1d d_ilist;
+  typename AT::t_int_1d d_numneigh;
+  DAT::tdual_float_2d k_result;
+  typename AT::t_float_2d d_result;
+};
+
+}    // namespace LAMMPS_NS
+
+#endif
+#endif

unittest/formats/CMakeLists.txt

@@ -41,6 +41,8 @@ set_tests_properties(TextFileReader PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${
 add_executable(test_file_operations test_file_operations.cpp)
 target_link_libraries(test_file_operations PRIVATE lammps GTest::GMock)
 add_test(NAME FileOperations COMMAND test_file_operations)
+# try to mitigate possible OpenMPI bug
+set_tests_properties(TextFileReader PROPERTIES ENVIRONMENT "OMPI_MCA_sharedfp=\"^sm\"")
 
 add_executable(test_dump_atom test_dump_atom.cpp)
 target_link_libraries(test_dump_atom PRIVATE lammps GTest::GMock)