Merge branch 'master' into package-reorganization-step1

# Conflicts:
#	doc/src/Packages_details.rst

doc/src/Commands_compute.rst

@@ -60,6 +60,7 @@ KOKKOS, o = OPENMP, t = OPT.
    * :doc:`erotate/sphere <compute_erotate_sphere>`
    * :doc:`erotate/sphere/atom <compute_erotate_sphere_atom>`
    * :doc:`event/displace <compute_event_displace>`
+   * :doc:`fabric <compute_fabric>`
    * :doc:`fep <compute_fep>`
    * :doc:`force/tally <compute_tally>`
    * :doc:`fragment/atom <compute_cluster_atom>`

doc/src/Howto_granular.rst

@@ -17,6 +17,11 @@ This compute
 * :doc:`compute erotate/sphere <compute_erotate_sphere>`
 
 calculates rotational kinetic energy which can be :doc:`output with thermodynamic info <Howto_output>`.
+The compute
+
+* :doc:`compute fabric <compute_fabric>`
+
+calculates various versions of the fabric tensor for granular and non-granular pair styles.
 
 Use one of these 4 pair potentials, which compute forces and torques
 between interacting pairs of particles:

doc/src/Library_create.rst

@@ -9,6 +9,7 @@ This section documents the following functions:
 - :cpp:func:`lammps_close`
 - :cpp:func:`lammps_mpi_init`
 - :cpp:func:`lammps_mpi_finalize`
+- :cpp:func:`lammps_kokkos_finalize`
 
 --------------------
 

doc/src/Python_create.rst

@@ -134,7 +134,10 @@ compiled with.
 The :py:func:`lmp.close() <lammps.lammps.close()>` call is
 optional since the LAMMPS class instance will also be deleted
 automatically during the :py:class:`lammps <lammps.lammps>` class
-destructor.
+destructor.  Instead of :py:func:`lmp.close() <lammps.lammps.close()>`
+it is also possible to call :py:func:`lmp.finalize() <lammps.lammps.finalize()>`;
+this will destruct the LAMMPS instance, but also finalized and release
+the MPI and/or Kokkos environment if enabled and active.
 
 Note that you can create multiple LAMMPS objects in your Python
 script, and coordinate and run multiple simulations, e.g.

doc/src/compute.rst

@@ -206,6 +206,7 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` doc
 * :doc:`erotate/sphere <compute_erotate_sphere>` - rotational energy of spherical particles
 * :doc:`erotate/sphere/atom <compute_erotate_sphere_atom>` - rotational energy for each spherical particle
 * :doc:`event/displace <compute_event_displace>` - detect event on atom displacement
+* :doc:`fabric <compute_fabric>` - calculates fabric tensors from pair interactions
 * :doc:`fep <compute_fep>` -
 * :doc:`force/tally <compute_tally>` -
 * :doc:`fragment/atom <compute_cluster_atom>` - fragment ID for each atom

doc/src/compute_fabric.rst

@@ -0,0 +1,185 @@
+.. index:: compute fabric
+
+compute fabric command
+======================
+
+Syntax
+""""""
+
+.. parsed-literal::
+
+   compute ID group-ID fabric cutoff attribute1 attribute2 ... keyword values ...
+
+* ID, group-ID are documented in :doc:`compute <compute>` command
+* fabric = style name of this compute command
+* cutoff = *type* or *radius*
+
+  .. parsed-literal::
+
+       *type* = cutoffs determined based on atom types
+       *radius* = cutoffs determined based on atom diameters (atom style sphere)
+
+* one or more attributes may be appended
+
+  .. parsed-literal::
+
+       *contact* = contact tensor
+       *branch* = branch tensor
+       *force/normal* = normal force tensor
+       *force/tangential* = tangential force tensor
+
+* zero or more keyword/value pairs may be appended
+* keyword = *type/include*
+
+  .. parsed-literal::
+
+       *type/include* value = arg1 arg2
+         arg = separate lists of types (see below)
+
+Examples
+""""""""
+
+.. code-block:: LAMMPS
+
+   compute 1 all fabric type contact force/normal type/include 1,2 3*4
+   compute 1 all fabric radius force/normal force/tangential
+
+Description
+"""""""""""
+
+Define a compute that calculates various fabric tensors for pairwise
+interaction :ref:`(Ouadfel) <Ouadfel>`. Fabric tensors are commonly used
+to quantify the anisotropy or orientation of granular contacts but can also
+be used to characterize the direction of pairwise interactions in general
+systems. The *type* and *radius* settings are used to select whether interactions
+cutoffs are determined by atom types or by the sum of atomic radii (atom
+style sphere), respectively. Calling this compute is roughly the cost of a
+pair style invocation as it involves a loop over the neighbor list. If the
+normal or tangential force tensors are requested, it will be more expensive
+than a pair style invocation as it will also recalculate all pair forces.
+
+Four fabric tensors are available: the contact, branch, normal force, or
+tangential force tensor. The contact tensor is calculated as
+
+.. math::
+
+   C_{ab}  =  \frac{15}{2} (\phi_{ab} - \mathrm{tr}(\phi) \delta_{ab})
+
+where :math:`a` and :math:`b` are the :math:`x`, :math:`y`, :math:`z`
+directions, :math:`\delta_{ab}` is the Kronecker delta function, and
+the tensor :math:`\phi` is defined as
+
+.. math::
+
+   \phi_{ab}  =  \sum_{n = 1}^{N_p} \frac{r_{a} r_{b}}{r^2}
+
+where :math:`n` loops over the :math:`N_p` pair interactions in the simulation,
+:math:`r_{a}` is the :math:`a` component of the radial vector between the
+two pairwise interacting particles, and :math:`r` is the magnitude of the radial vector.
+
+The branch tensor is calculated as
+
+.. math::
+
+   B_{ab}  =  \frac{15}{6 \mathrm{tr}(D)} (D_{ab} - \mathrm{tr}(D) \delta_{ab})
+
+where the tensor :math:`D` is defined as
+
+.. math::
+
+   D_{ab}  =  \sum_{n = 1}^{N_p}
+                \frac{1}{N_c (r^2 + C_{cd} r_c r_d)}
+                \frac{r_{a} r_{b}}{r}
+
+where :math:`N_c` is the total number of contacts in the system and the subscripts
+:math:`c` and :math:`d` indices are summed according to Einstein notation.
+
+The normal force fabric tensor is calculated as
+
+.. math::
+
+   F^n_{ab}  =  \frac{15}{6 \mathrm{tr}(N)} (N_{ab} - \mathrm{tr}(N) \delta_{ab})
+
+where the tensor :math:`N` is defined as
+
+.. math::
+
+   N_{ab}  =  \sum_{n = 1}^{N_p}
+                \frac{1}{N_c (r^2 + C_{cd} r_c r_d)}
+                \frac{r_{a} r_{b}}{r^2} f_n
+
+and :math:`f_n` is the magnitude of the normal, central-body force between the two atoms.
+
+Finally, the tangential force fabric tensor is only defined for pair styles that
+apply tangential forces to particles, namely granular pair styles. It is calculated
+as
+
+.. math::
+
+   F^t_{ab}  =  \frac{15}{9 \mathrm{tr}(N)} (T_{ab} - \mathrm{tr}(T) \delta_{ab})
+
+where the tensor :math:`T` is defined as
+
+.. math::
+
+   T_{ab}  =  \sum_{n = 1}^{N_p}
+                \frac{1}{N_c (r^2 + C_{cd} r_c r_d)}
+                \frac{r_{a} r_{b}}{r^2} f_t
+
+and :math:`f_t` is the magnitude of the tangential force between the two atoms.
+
+The *type/include* keyword filters interactions based on the types of the two atoms.
+Interactions between two atoms are only included in calculations if the atom types
+are in the two lists. Each list consists of a series of type
+ranges separated by commas. The range can be specified as a
+single numeric value, or a wildcard asterisk can be used to specify a range
+of values.  This takes the form "\*" or "\*n" or "n\*" or "m\*n".  For
+example, if M = the number of atom types, then an asterisk with no numeric
+values means all types from 1 to M.  A leading asterisk means all types
+from 1 to n (inclusive).  A trailing asterisk means all types from n to M
+(inclusive).  A middle asterisk means all types from m to n (inclusive).
+Multiple *type/include* keywords may be added.
+
+Output info
+"""""""""""
+
+This compute calculates a local vector of doubles and a scalar. The vector stores the
+unique components of the first requested tensor in the order xx, yy, zz, xy, xz, yz
+followed by the same components for all subsequent tensors. The length of the vector
+is therefore six times the number of requested tensors. The scalar output is the
+number of pairwise interactions included in the calculation of the fabric tensor.
+
+Restrictions
+""""""""""""
+
+This fix is part of the GRANULAR package.  It is only enabled if LAMMPS
+was built with that package.  See the :doc:`Build package <Build_package>`
+doc page for more info.
+
+Currently, compute *fabric* does not support pair styles
+with many-body interactions.  It also does not
+support models with long-range Coulombic or dispersion forces,
+i.e. the kspace_style command in LAMMPS.  It also does not support the
+following fixes which add rigid-body constraints: :doc:`fix shake
+<fix_shake>`, :doc:`fix rattle <fix_shake>`, :doc:`fix rigid
+<fix_rigid>`, :doc:`fix rigid/small <fix_rigid>`. It does not support
+granular pair styles that extend beyond the contact of atomic radii
+(e.g. JKR and DMT).
+
+Related commands
+""""""""""""""""
+
+none
+
+Default
+"""""""
+
+none
+
+----------
+
+.. _Ouadfel:
+
+**(Ouadfel)** Ouadfel and Rothenburg
+"Stress-force-fabric relationship for assemblies of ellipsoids",
+Mechanics of Materials (2001). (`link to paper <https://doi.org/10.1016/S0167-6636(00)00057-0>`_)

doc/utils/sphinx-config/false_positives.txt

@@ -2680,6 +2680,7 @@ qoverride
 qqr
 qqrd
 qtb
+Quadfel
 quadratically
 quadrupolar
 Quant
@@ -2861,6 +2862,7 @@ Rossky
 rosybrown
 rotationally
 Rotenberg
+Rothenburg
 Rovigatti
 royalblue
 rozero

fortran/lammps.f90

@@ -76,17 +76,15 @@ MODULE LIBLAMMPS
         TYPE(c_ptr), VALUE :: handle
       END SUBROUTINE lammps_close
 
-      SUBROUTINE lammps_mpi_init(handle) BIND(C, name='lammps_mpi_init')
+      SUBROUTINE lammps_mpi_init() BIND(C, name='lammps_mpi_init')
-        IMPORT :: c_ptr
-        TYPE(c_ptr), VALUE :: handle
       END SUBROUTINE lammps_mpi_init
 
-      SUBROUTINE lammps_mpi_finalize(handle) &
+      SUBROUTINE lammps_mpi_finalize() BIND(C, name='lammps_mpi_finalize')
-          BIND(C, name='lammps_mpi_finalize')
-        IMPORT :: c_ptr
-        TYPE(c_ptr), VALUE :: handle
       END SUBROUTINE lammps_mpi_finalize
 
+      SUBROUTINE lammps_kokkos_finalize() BIND(C, name='lammps_kokkos_finalize')
+      END SUBROUTINE lammps_kokkos_finalize
+
       SUBROUTINE lammps_file(handle,filename) BIND(C, name='lammps_file')
         IMPORT :: c_ptr
         TYPE(c_ptr), VALUE :: handle
@@ -188,7 +186,8 @@ CONTAINS
 
     IF (PRESENT(finalize)) THEN
         IF (finalize) THEN
-            CALL lammps_mpi_finalize(self%handle)
+            CALL lammps_kokkos_finalize()
+            CALL lammps_mpi_finalize()
         END IF
     END IF
   END SUBROUTINE lmp_close

python/lammps/core.py

@@ -460,10 +460,16 @@ class lammps(object):
   # -------------------------------------------------------------------------
 
   def finalize(self):
-    """Shut down the MPI communication through the library interface by calling :cpp:func:`lammps_finalize`.
+    """Shut down the MPI communication and Kokkos environment (if active) through the
+       library interface by  calling :cpp:func:`lammps_mpi_finalize` and
+       :cpp:func:`lammps_kokkos_finalize`.
+
+       You cannot create or use any LAMMPS instances after this function is called
+       unless LAMMPS was compiled without MPI and without Kokkos support.
     """
     self.close()
-    self.lib.lammps_finalize()
+    self.lib.lammps_kokkos_finalize()
+    self.lib.lammps_mpi_finalize()
 
   # -------------------------------------------------------------------------
 

src/.gitignore

@@ -386,6 +386,8 @@
 /compute_erotate_rigid.h
 /compute_event_displace.cpp
 /compute_event_displace.h
+/compute_fabric.cpp
+/compute_fabric.h
 /compute_fep.cpp
 /compute_fep.h
 /compute_force_tally.cpp

src/GRANULAR/compute_fabric.cpp

@@ -0,0 +1,595 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "compute_fabric.h"
+
+#include "atom.h"
+#include "error.h"
+#include "force.h"
+#include "memory.h"
+#include "modify.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "neighbor.h"
+#include "pair.h"
+#include "tokenizer.h"
+#include "update.h"
+
+#include <cmath>
+#include <cstring>
+
+using namespace LAMMPS_NS;
+
+enum { OTHER, GRANULAR };
+enum { TYPE, RADIUS };
+enum { CN, BR, FN, FT };
+
+/* ---------------------------------------------------------------------- */
+
+ComputeFabric::ComputeFabric(LAMMPS *lmp, int narg, char **arg) :
+    Compute(lmp, narg, arg), tensor_style(NULL)
+{
+  if (narg < 3) error->all(FLERR, "Illegal compute fabric command");
+
+  if (strcmp(arg[3], "type") == 0)
+    cutstyle = TYPE;
+  else if (strcmp(arg[3], "radius") == 0)
+    cutstyle = RADIUS;
+  else
+    error->all(FLERR, "Illegal compute fabric command");
+
+  if (cutstyle == RADIUS && !atom->radius_flag)
+    error->all(FLERR, "Compute fabric radius style requires atom attribute radius");
+
+  // If optional arguments included, this will be oversized
+  ntensors = narg - 4;
+  tensor_style = new int[ntensors];
+
+  cn_flag = 0;
+  br_flag = 0;
+  fn_flag = 0;
+  ft_flag = 0;
+  type_filter = nullptr;
+
+  ntensors = 0;
+  int iarg = 4;
+  while (iarg < narg) {
+    if (strcmp(arg[iarg], "contact") == 0) {
+      cn_flag = 1;
+      tensor_style[ntensors++] = CN;
+    } else if (strcmp(arg[iarg], "branch") == 0) {
+      br_flag = 1;
+      tensor_style[ntensors++] = BR;
+    } else if (strcmp(arg[iarg], "force/normal") == 0) {
+      fn_flag = 1;
+      tensor_style[ntensors++] = FN;
+    } else if (strcmp(arg[iarg], "force/tangential") == 0) {
+      ft_flag = 1;
+      tensor_style[ntensors++] = FT;
+    } else if (strcmp(arg[iarg], "type/include") == 0) {
+      if (iarg + 1 >= narg) error->all(FLERR, "Invalid keyword in compute fabric command");
+      int ntypes = atom->ntypes;
+
+      int i, j, itype, jtype, in, jn, infield, jnfield;
+      int inlo, inhi, jnlo, jnhi;
+      char *istr, *jstr;
+      if (!type_filter) {
+        memory->create(type_filter, ntypes + 1, ntypes + 1, "compute/fabric:type_filter");
+
+        for (i = 0; i <= ntypes; i++) {
+          for (j = 0; j <= ntypes; j++) { type_filter[i][j] = 0; }
+        }
+      }
+
+      in = strlen(arg[iarg + 1]) + 1;
+      istr = new char[in];
+      strcpy(istr, arg[iarg + 1]);
+      std::vector<std::string> iwords = Tokenizer(istr, ",").as_vector();
+      infield = iwords.size();
+
+      jn = strlen(arg[iarg + 2]) + 1;
+      jstr = new char[jn];
+      strcpy(jstr, arg[iarg + 2]);
+      std::vector<std::string> jwords = Tokenizer(jstr, ",").as_vector();
+      jnfield = jwords.size();
+
+      for (i = 0; i < infield; i++) {
+        const char *ifield = iwords[i].c_str();
+        utils::bounds(FLERR, ifield, 1, ntypes, inlo, inhi, error);
+
+        for (j = 0; j < jnfield; j++) {
+          const char *jfield = jwords[j].c_str();
+          utils::bounds(FLERR, jfield, 1, ntypes, jnlo, jnhi, error);
+
+          for (itype = inlo; itype <= inhi; itype++) {
+            for (jtype = jnlo; jtype <= jnhi; jtype++) {
+              type_filter[itype][jtype] = 1;
+              type_filter[jtype][itype] = 1;
+            }
+          }
+        }
+      }
+
+      delete[] istr;
+      delete[] jstr;
+
+      iarg += 2;
+    } else
+      error->all(FLERR, "Illegal compute fabric command");
+    iarg++;
+  }
+
+  vector_flag = 1;
+  size_vector = ntensors * 6;
+  extvector = 0;
+
+  scalar_flag = 1;
+  extscalar = 1;
+
+  vector = new double[size_vector];
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeFabric::~ComputeFabric()
+{
+  delete[] vector;
+  delete[] tensor_style;
+  memory->destroy(type_filter);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeFabric::init()
+{
+  if (force->pair == NULL) error->all(FLERR, "No pair style is defined for compute fabric");
+  if (force->pair->single_enable == 0 && (fn_flag || ft_flag))
+    error->all(FLERR, "Pair style does not support compute fabric normal or tangential force");
+
+  // Find if granular or gran
+  pstyle = OTHER;
+  if (force->pair_match("^granular", 0) || force->pair_match("^gran/", 0)) pstyle = GRANULAR;
+
+  if (pstyle != GRANULAR && ft_flag)
+    error->all(FLERR, "Pair style does not calculate tangential forces for compute fabric");
+
+  if (force->pair->beyond_contact)
+    error->all(FLERR, "Compute fabric does not support pair styles that extend beyond contact");
+
+  // need an occasional half neighbor list
+  // set size to same value as request made by force->pair
+  // this should enable it to always be a copy list (e.g. for granular pstyle)
+
+  int irequest = neighbor->request(this, instance_me);
+  neighbor->requests[irequest]->pair = 0;
+  neighbor->requests[irequest]->compute = 1;
+  neighbor->requests[irequest]->occasional = 1;
+  NeighRequest *pairrequest = neighbor->find_request((void *) force->pair);
+  if (pairrequest) neighbor->requests[irequest]->size = pairrequest->size;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeFabric::init_list(int /*id*/, NeighList *ptr)
+{
+  list = ptr;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeFabric::compute_vector()
+{
+  invoked_vector = update->ntimestep;
+
+  int i, j, ii, jj, inum, jnum, itype, jtype;
+  tagint itag, jtag;
+  double xtmp, ytmp, ztmp, delx, dely, delz;
+  double r, rinv, rsq, radsum, eng, fpair;
+
+  double nx, ny, nz;
+  double ncinv, denom, fn, ft, prefactor;
+  double br_tensor[6], ft_tensor[6], fn_tensor[6];
+  double trace_phi, trace_D, trace_Xfn, trace_Xft;
+  double phi_ij[6] = {0.0};
+  double Ac_ij[6] = {0.0};
+  double D_ij[6] = {0.0};
+  double Xfn_ij[6] = {0.0};
+  double Xft_ij[6] = {0.0};
+  double temp_dbl[6];
+
+  int *ilist, *jlist, *numneigh, **firstneigh;
+
+  double **x = atom->x;
+  double *radius = atom->radius;
+  tagint *tag = atom->tag;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+  int newton_pair = force->newton_pair;
+
+  // invoke half neighbor list (will copy or build if necessary)
+  neighbor->build_one(list);
+
+  inum = list->inum;
+  ilist = list->ilist;
+  numneigh = list->numneigh;
+  firstneigh = list->firstneigh;
+
+  Pair *pair = force->pair;
+  double **cutsq = force->pair->cutsq;
+
+  // invoke compute_scalar() to update the number of contacts, if needed
+  nc = compute_scalar();
+
+  // If no contacts, everything will be zero
+  if (nc == 0) {
+    for (i = 0; i < size_vector; i++) vector[i] = 0.0;
+    return;
+  }
+  ncinv = 1.0 / nc;
+
+  // First loop through and calculate contact fabric tensor
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    if (!(mask[i] & groupbit)) continue;
+
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+    itag = tag[i];
+    itype = type[i];
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+
+    for (jj = 0; jj < jnum; jj++) {
+      j = jlist[jj];
+      j &= NEIGHMASK;
+
+      if (!(mask[j] & groupbit)) continue;
+
+      // itag = jtag is possible for long cutoffs that include images of self
+
+      if (newton_pair == 0 && j >= nlocal) {
+        jtag = tag[j];
+        if (itag > jtag) {
+          if ((itag + jtag) % 2 == 0) continue;
+        } else if (itag < jtag) {
+          if ((itag + jtag) % 2 == 1) continue;
+        } else {
+          if (x[j][2] < ztmp) continue;
+          if (x[j][2] == ztmp) {
+            if (x[j][1] < ytmp) continue;
+            if (x[j][1] == ytmp && x[j][0] < xtmp) continue;
+          }
+        }
+      }
+
+      jtype = type[j];
+
+      if (type_filter)
+        if (type_filter[itype][jtype] == 0) continue;
+
+      delx = xtmp - x[j][0];
+      dely = ytmp - x[j][1];
+      delz = ztmp - x[j][2];
+      rsq = delx * delx + dely * dely + delz * delz;
+
+      if (cutstyle == TYPE) {
+        if (rsq >= cutsq[itype][jtype]) continue;
+      } else {
+        radsum = radius[i] + radius[j];
+        if (rsq >= radsum * radsum) continue;
+      }
+
+      r = sqrt(rsq);
+      rinv = 1.0 / r;
+      nx = delx * rinv;
+      ny = dely * rinv;
+      nz = delz * rinv;
+
+      phi_ij[0] += nx * nx;
+      phi_ij[1] += ny * ny;
+      phi_ij[2] += nz * nz;
+      phi_ij[3] += nx * ny;
+      phi_ij[4] += nx * nz;
+      phi_ij[5] += ny * nz;
+    }
+  }
+
+  //Sum phi across processors
+  MPI_Allreduce(phi_ij, temp_dbl, 6, MPI_DOUBLE, MPI_SUM, world);
+  for (i = 0; i < 6; i++) phi_ij[i] = temp_dbl[i] * ncinv;
+
+  trace_phi = (1.0 / 3.0) * (phi_ij[0] + phi_ij[1] + phi_ij[2]);
+
+  Ac_ij[0] = (15.0 / 2.0) * (phi_ij[0] - trace_phi);
+  Ac_ij[1] = (15.0 / 2.0) * (phi_ij[1] - trace_phi);
+  Ac_ij[2] = (15.0 / 2.0) * (phi_ij[2] - trace_phi);
+  Ac_ij[3] = (15.0 / 2.0) * (phi_ij[3]);
+  Ac_ij[4] = (15.0 / 2.0) * (phi_ij[4]);
+  Ac_ij[5] = (15.0 / 2.0) * (phi_ij[5]);
+
+  // If needed, loop through and calculate other fabric tensors
+  if (br_flag || fn_flag || ft_flag) {
+
+    for (ii = 0; ii < inum; ii++) {
+      i = ilist[ii];
+      if (!(mask[i] & groupbit)) continue;
+
+      xtmp = x[i][0];
+      ytmp = x[i][1];
+      ztmp = x[i][2];
+      itag = tag[i];
+      itype = type[i];
+      jlist = firstneigh[i];
+      jnum = numneigh[i];
+
+      for (jj = 0; jj < jnum; jj++) {
+        j = jlist[jj];
+        j &= NEIGHMASK;
+
+        if (!(mask[j] & groupbit)) continue;
+
+        // itag = jtag is possible for long cutoffs that include images of self
+
+        if (newton_pair == 0 && j >= nlocal) {
+          jtag = tag[j];
+          if (itag > jtag) {
+            if ((itag + jtag) % 2 == 0) continue;
+          } else if (itag < jtag) {
+            if ((itag + jtag) % 2 == 1) continue;
+          } else {
+            if (x[j][2] < ztmp) continue;
+            if (x[j][2] == ztmp) {
+              if (x[j][1] < ytmp) continue;
+              if (x[j][1] == ytmp && x[j][0] < xtmp) continue;
+            }
+          }
+        }
+
+        jtype = type[j];
+
+        if (type_filter)
+          if (type_filter[itype][jtype] == 0) continue;
+
+        delx = xtmp - x[j][0];
+        dely = ytmp - x[j][1];
+        delz = ztmp - x[j][2];
+        rsq = delx * delx + dely * dely + delz * delz;
+
+        if (cutstyle == TYPE) {
+          if (rsq >= cutsq[itype][jtype]) continue;
+        } else {
+          radsum = radius[i] + radius[j];
+          if (rsq >= radsum * radsum) continue;
+        }
+
+        if (fn_flag || ft_flag) eng = pair->single(i, j, itype, jtype, rsq, 1.0, 1.0, fpair);
+
+        r = sqrt(rsq);
+        rinv = 1.0 / r;
+        nx = delx * rinv;
+        ny = dely * rinv;
+        nz = delz * rinv;
+
+        denom = 1 + Ac_ij[0] * nx * nx + Ac_ij[1] * ny * ny + Ac_ij[2] * nz * nz;
+        denom += 2 * Ac_ij[3] * nx * ny + 2 * Ac_ij[4] * nx * nz + 2 * Ac_ij[5] * ny * nz;
+        prefactor = ncinv / denom;
+
+        if (br_flag) {
+          D_ij[0] += prefactor * nx * nx * r;
+          D_ij[1] += prefactor * ny * ny * r;
+          D_ij[2] += prefactor * nz * nz * r;
+          D_ij[3] += prefactor * nx * ny * r;
+          D_ij[4] += prefactor * nx * nz * r;
+          D_ij[5] += prefactor * ny * nz * r;
+        }
+
+        if (fn_flag || ft_flag) {
+          fn = r * fpair;
+
+          Xfn_ij[0] += prefactor * nx * nx * fn;
+          Xfn_ij[1] += prefactor * ny * ny * fn;
+          Xfn_ij[2] += prefactor * nz * nz * fn;
+          Xfn_ij[3] += prefactor * nx * ny * fn;
+          Xfn_ij[4] += prefactor * nx * nz * fn;
+          Xfn_ij[5] += prefactor * ny * nz * fn;
+
+          if (ft_flag) {
+            ft = force->pair->svector[3];
+
+            Xft_ij[0] += prefactor * nx * nx * ft;
+            Xft_ij[1] += prefactor * ny * ny * ft;
+            Xft_ij[2] += prefactor * nz * nz * ft;
+            Xft_ij[3] += prefactor * nx * ny * ft;
+            Xft_ij[4] += prefactor * nx * nz * ft;
+            Xft_ij[5] += prefactor * ny * nz * ft;
+          }
+        }
+      }
+    }
+  }
+
+  // Output results
+
+  if (cn_flag) {
+    for (i = 0; i < ntensors; i++) {
+      if (tensor_style[i] == CN) {
+        for (j = 0; j < 6; j++) vector[6 * i + j] = Ac_ij[j];
+      }
+    }
+  }
+
+  if (br_flag) {
+    MPI_Allreduce(D_ij, temp_dbl, 6, MPI_DOUBLE, MPI_SUM, world);
+    for (i = 0; i < 6; i++) D_ij[i] = temp_dbl[i];
+
+    trace_D = (1.0 / 3.0) * (D_ij[0] + D_ij[1] + D_ij[2]);
+
+    br_tensor[0] = (15.0 / (6.0 * trace_D)) * (D_ij[0] - trace_D);
+    br_tensor[1] = (15.0 / (6.0 * trace_D)) * (D_ij[1] - trace_D);
+    br_tensor[2] = (15.0 / (6.0 * trace_D)) * (D_ij[2] - trace_D);
+    br_tensor[3] = (15.0 / (6.0 * trace_D)) * (D_ij[3]);
+    br_tensor[4] = (15.0 / (6.0 * trace_D)) * (D_ij[4]);
+    br_tensor[5] = (15.0 / (6.0 * trace_D)) * (D_ij[5]);
+
+    for (i = 0; i < ntensors; i++) {
+      if (tensor_style[i] == BR) {
+        for (j = 0; j < 6; j++) vector[6 * i + j] = br_tensor[j];
+      }
+    }
+  }
+
+  if (fn_flag || ft_flag) {
+    MPI_Allreduce(Xfn_ij, temp_dbl, 6, MPI_DOUBLE, MPI_SUM, world);
+    for (i = 0; i < 6; i++) Xfn_ij[i] = temp_dbl[i];
+
+    trace_Xfn = (1.0 / 3.0) * (Xfn_ij[0] + Xfn_ij[1] + Xfn_ij[2]);
+  }
+
+  if (fn_flag) {
+
+    fn_tensor[0] = (15.0 / (6.0 * trace_Xfn)) * (Xfn_ij[0] - trace_Xfn);
+    fn_tensor[1] = (15.0 / (6.0 * trace_Xfn)) * (Xfn_ij[1] - trace_Xfn);
+    fn_tensor[2] = (15.0 / (6.0 * trace_Xfn)) * (Xfn_ij[2] - trace_Xfn);
+    fn_tensor[3] = (15.0 / (6.0 * trace_Xfn)) * (Xfn_ij[3]);
+    fn_tensor[4] = (15.0 / (6.0 * trace_Xfn)) * (Xfn_ij[4]);
+    fn_tensor[5] = (15.0 / (6.0 * trace_Xfn)) * (Xfn_ij[5]);
+
+    for (i = 0; i < ntensors; i++) {
+      if (tensor_style[i] == FN) {
+        for (j = 0; j < 6; j++) vector[6 * i + j] = fn_tensor[j];
+      }
+    }
+  }
+
+  if (ft_flag) {
+    MPI_Allreduce(Xft_ij, temp_dbl, 6, MPI_DOUBLE, MPI_SUM, world);
+    for (i = 0; i < 6; i++) Xft_ij[i] = temp_dbl[i];
+
+    trace_Xft = (1.0 / 3.0) * (Xft_ij[0] + Xft_ij[1] + Xft_ij[2]);
+
+    ft_tensor[0] = (15.0 / (9.0 * trace_Xfn)) * (Xft_ij[0] - trace_Xft);
+    ft_tensor[1] = (15.0 / (9.0 * trace_Xfn)) * (Xft_ij[1] - trace_Xft);
+    ft_tensor[2] = (15.0 / (9.0 * trace_Xfn)) * (Xft_ij[2] - trace_Xft);
+    ft_tensor[3] = (15.0 / (9.0 * trace_Xfn)) * (Xft_ij[3]);
+    ft_tensor[4] = (15.0 / (9.0 * trace_Xfn)) * (Xft_ij[4]);
+    ft_tensor[5] = (15.0 / (9.0 * trace_Xfn)) * (Xft_ij[5]);
+
+    for (i = 0; i < ntensors; i++) {
+      if (tensor_style[i] == FT) {
+        for (j = 0; j < 6; j++) vector[6 * i + j] = ft_tensor[j];
+      }
+    }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+double ComputeFabric::compute_scalar()
+{
+  // Skip if already calculated on this timestep
+  if (invoked_scalar == update->ntimestep) return nc;
+
+  invoked_scalar = update->ntimestep;
+
+  int i, j, ii, jj, inum, jnum, itype, jtype;
+  tagint itag, jtag;
+  double xtmp, ytmp, ztmp, delx, dely, delz;
+  double rsq, radsum, temp_dbl;
+
+  int *ilist, *jlist, *numneigh, **firstneigh;
+
+  double **x = atom->x;
+  double *radius = atom->radius;
+  tagint *tag = atom->tag;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+  int newton_pair = force->newton_pair;
+
+  // invoke half neighbor list (will copy or build if necessary)
+  neighbor->build_one(list);
+
+  inum = list->inum;
+  ilist = list->ilist;
+  numneigh = list->numneigh;
+  firstneigh = list->firstneigh;
+
+  double **cutsq = force->pair->cutsq;
+
+  // First loop through and calculate contact fabric tensor
+  nc = 0;
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    if (!(mask[i] & groupbit)) continue;
+
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+    itag = tag[i];
+    itype = type[i];
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+
+    for (jj = 0; jj < jnum; jj++) {
+      j = jlist[jj];
+      j &= NEIGHMASK;
+
+      if (!(mask[j] & groupbit)) continue;
+
+      // itag = jtag is possible for long cutoffs that include images of self
+
+      if (newton_pair == 0 && j >= nlocal) {
+        jtag = tag[j];
+        if (itag > jtag) {
+          if ((itag + jtag) % 2 == 0) continue;
+        } else if (itag < jtag) {
+          if ((itag + jtag) % 2 == 1) continue;
+        } else {
+          if (x[j][2] < ztmp) continue;
+          if (x[j][2] == ztmp) {
+            if (x[j][1] < ytmp) continue;
+            if (x[j][1] == ytmp && x[j][0] < xtmp) continue;
+          }
+        }
+      }
+
+      jtype = type[j];
+
+      if (type_filter)
+        if (type_filter[itype][jtype] == 0) continue;
+
+      delx = xtmp - x[j][0];
+      dely = ytmp - x[j][1];
+      delz = ztmp - x[j][2];
+      rsq = delx * delx + dely * dely + delz * delz;
+
+      if (cutstyle == TYPE) {
+        if (rsq >= cutsq[itype][jtype]) continue;
+      } else {
+        radsum = radius[i] + radius[j];
+        if (rsq >= radsum * radsum) continue;
+      }
+
+      nc += 1.0;
+    }
+  }
+  //Count total contacts across processors
+  MPI_Allreduce(&nc, &temp_dbl, 1, MPI_DOUBLE, MPI_SUM, world);
+  nc = temp_dbl;
+
+  scalar = nc;
+  return nc;
+}
+

src/GRANULAR/compute_fabric.h

@@ -0,0 +1,79 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   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(fabric,ComputeFabric);
+// clang-format on
+#else
+
+#ifndef LMP_COMPUTE_FABRIC_H
+#define LMP_COMPUTE_FABRIC_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeFabric : public Compute {
+ public:
+  ComputeFabric(class LAMMPS *, int, char **);
+  ~ComputeFabric();
+  void init();
+  void init_list(int, class NeighList *);
+  void compute_vector();
+  double compute_scalar();
+
+ private:
+  int ntensors, pstyle, cutstyle;
+  double nc;
+  int *tensor_style;
+  int **type_filter;
+  class NeighList *list;
+
+  int cn_flag, br_flag, fn_flag, ft_flag;
+};
+
+}    // namespace LAMMPS_NS
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Compute fabric radius style requires atom attribute radius
+
+Self-explanatory.
+
+E: No pair style is defined for compute fabric
+
+Self-explanatory.
+
+E: Pair style does not support compute fabric normal or tangential force
+
+Pair style must be single enabled to calculate the normal or tangential force tensors
+
+E: Pair style does not calculate tangential forces for compute fabric
+
+The tangential force tensor can only be calculated for granular pair styles with tangential forces
+
+E: Compute fabric does not support pair styles that extend beyond contact
+
+Granular pair styles that extend beyond contact such as JKR or DMT are not supported
+
+*/

src/KOKKOS/kokkos.cpp

@@ -69,6 +69,10 @@ GPU_AWARE_UNKNOWN
 
 using namespace LAMMPS_NS;
 
+Kokkos::InitArguments KokkosLMP::args{-1, -1, -1, false};
+int KokkosLMP::is_finalized = 0;
+int KokkosLMP::init_ngpus = 0;
+
 /* ---------------------------------------------------------------------- */
 
 KokkosLMP::KokkosLMP(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
@@ -155,6 +159,10 @@ KokkosLMP::KokkosLMP(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
     } else if (strcmp(arg[iarg],"t") == 0 ||
                strcmp(arg[iarg],"threads") == 0) {
       nthreads = atoi(arg[iarg+1]);
+
+      if (nthreads <= 0)
+        error->all(FLERR,"Invalid number of threads requested for Kokkos: must be 1 or greater");
+
       iarg += 2;
 
     } else if (strcmp(arg[iarg],"n") == 0 ||
@@ -165,13 +173,27 @@ KokkosLMP::KokkosLMP(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
     } else error->all(FLERR,"Invalid Kokkos command-line args");
   }
 
-  // initialize Kokkos
+  // Initialize Kokkos. However, we cannot change any
+  // Kokkos library parameters after the first initalization
 
-  if (me == 0) {
+  if (args.num_threads != -1) {
-    if (screen) fprintf(screen,"  will use up to %d GPU(s) per node\n",ngpus);
+    if (args.num_threads != nthreads || args.num_numa != numa || args.device_id != device)
-    if (logfile) fprintf(logfile,"  will use up to %d GPU(s) per node\n",ngpus);
+      if (me == 0)
+        error->warning(FLERR,"Kokkos package already initalized, cannot reinitialize with different parameters");
+    nthreads = args.num_threads;
+    numa = args.num_numa;
+    device = args.device_id;
+    ngpus = init_ngpus;
+  } else {
+    args.num_threads = nthreads;
+    args.num_numa = numa;
+    args.device_id = device;
+    init_ngpus = ngpus;
   }
 
+  if (me == 0)
+    utils::logmesg(lmp, "  will use up to {} GPU(s) per node\n",ngpus);
+
 #ifdef LMP_KOKKOS_GPU
   if (ngpus <= 0)
     error->all(FLERR,"Kokkos has been compiled for CUDA, HIP, or SYCL but no GPUs are requested");
@@ -184,12 +206,7 @@ KokkosLMP::KokkosLMP(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
                          "than the OpenMP backend");
 #endif
 
-  Kokkos::InitArguments args;
+  KokkosLMP::initialize(args,error);
-  args.num_threads = nthreads;
-  args.num_numa = numa;
-  args.device_id = device;
-
-  Kokkos::initialize(args);
 
   // default settings for package kokkos command
 
@@ -299,9 +316,27 @@ KokkosLMP::KokkosLMP(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
 
 KokkosLMP::~KokkosLMP()
 {
-  // finalize Kokkos
 
+}
+
+/* ---------------------------------------------------------------------- */
+
+void KokkosLMP::initialize(Kokkos::InitArguments args, Error *error)
+{
+  if (!Kokkos::is_initialized()) {
+    if (is_finalized)
+      error->all(FLERR,"Kokkos package already finalized, cannot re-initialize");
+    Kokkos::initialize(args);
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void KokkosLMP::finalize()
+{
+  if (Kokkos::is_initialized() && !is_finalized)
     Kokkos::finalize();
+  is_finalized = 1;
 }
 
 /* ----------------------------------------------------------------------

src/KOKKOS/kokkos.h

@@ -49,8 +49,14 @@ class KokkosLMP : protected Pointers {
   int newtonflag;
   double binsize;
 
+  static int is_finalized;
+  static Kokkos::InitArguments args;
+  static int init_ngpus;
+
   KokkosLMP(class LAMMPS *, int, char **);
   ~KokkosLMP();
+  static void initialize(Kokkos::InitArguments, Error *);
+  static void finalize();
   void accelerator(int, char **);
   int neigh_count(int);
 
@@ -84,13 +90,21 @@ because MPI library not recognized
 
 The local MPI rank was not found in one of four supported environment variables.
 
+E: Invalid number of threads requested for Kokkos: must be 1 or greater
+
+Self-explanatory.
+
 E: GPUs are requested but Kokkos has not been compiled for CUDA
 
 Recompile Kokkos with CUDA support to use GPUs.
 
-E: Kokkos has been compiled for CUDA but no GPUs are requested
+E: Kokkos has been compiled for CUDA, HIP, or SYCL but no GPUs are requested
 
-One or more GPUs must be used when Kokkos is compiled for CUDA.
+One or more GPUs must be used when Kokkos is compiled for CUDA/HIP/SYCL.
+
+W: Kokkos package already initalized, cannot reinitialize with different parameters
+
+Self-explanatory.
 
 E: Illegal ... command
 

src/SPIN/fix_nve_spin.cpp

@@ -176,6 +176,7 @@ void FixNVESpin::init()
 
   // loop 1: obtain # of Pairs, and # of Pair/Spin styles
 
+  npairspin = 0;
   PairHybrid *hybrid = (PairHybrid *)force->pair_match("^hybrid",0);
   if (force->pair_match("^spin",0,0)) {        // only one Pair/Spin style
     pair = force->pair_match("^spin",0,0);
@@ -232,6 +233,7 @@ void FixNVESpin::init()
   // loop 1: obtain # of fix precession/spin styles
 
   int iforce;
+  nprecspin = 0;
   for (iforce = 0; iforce < modify->nfix; iforce++) {
     if (utils::strmatch(modify->fix[iforce]->style,"^precession/spin")) {
       nprecspin++;
@@ -264,6 +266,7 @@ void FixNVESpin::init()
 
   // loop 1: obtain # of fix langevin/spin styles
 
+  nlangspin = 0;
   for (iforce = 0; iforce < modify->nfix; iforce++) {
     if (utils::strmatch(modify->fix[iforce]->style,"^langevin/spin")) {
       nlangspin++;

src/accelerator_kokkos.h

@@ -56,16 +56,12 @@ class KokkosLMP {
 
   KokkosLMP(class LAMMPS *, int, char **) { kokkos_exists = 0; }
   ~KokkosLMP() {}
+  static void finalize() {}
   void accelerator(int, char **) {}
   int neigh_list_kokkos(int) { return 0; }
   int neigh_count(int) { return 0; }
 };
 
-class Kokkos {
- public:
-  static void finalize() {}
-};
-
 class AtomKokkos : public Atom {
  public:
   tagint **k_special;

src/error.cpp

@@ -81,7 +81,7 @@ void Error::universe_all(const std::string &file, int line, const std::string &s
 
   throw LAMMPSException(mesg);
 #else
-  if (lmp->kokkos) Kokkos::finalize();
+  KokkosLMP::finalize();
   MPI_Finalize();
   exit(1);
 #endif
@@ -107,6 +107,7 @@ void Error::universe_one(const std::string &file, int line, const std::string &s
 
   throw LAMMPSAbortException(mesg, universe->uworld);
 #else
+  KokkosLMP::finalize();
   MPI_Abort(universe->uworld,1);
   exit(1); // to trick "smart" compilers into believing this does not return
 #endif
@@ -173,8 +174,8 @@ void Error::all(const std::string &file, int line, const std::string &str)
   if (screen && screen != stdout) fclose(screen);
   if (logfile) fclose(logfile);
 
+  KokkosLMP::finalize();
   if (universe->nworlds > 1) MPI_Abort(universe->uworld,1);
-  if (lmp->kokkos) Kokkos::finalize();
   MPI_Finalize();
   exit(1);
 #endif
@@ -213,6 +214,7 @@ void Error::one(const std::string &file, int line, const std::string &str)
 #else
   if (screen) fflush(screen);
   if (logfile) fflush(logfile);
+  KokkosLMP::finalize();
   MPI_Abort(world,1);
   exit(1); // to trick "smart" compilers into believing this does not return
 #endif
@@ -315,7 +317,7 @@ void Error::done(int status)
   if (screen && screen != stdout) fclose(screen);
   if (logfile) fclose(logfile);
 
-  if (lmp->kokkos) Kokkos::finalize();
+  KokkosLMP::finalize();
   MPI_Finalize();
   exit(status);
 }

src/library.cpp

@@ -19,6 +19,7 @@
 #include "library.h"
 #include <mpi.h>
 
+#include "accelerator_kokkos.h"
 #include "atom.h"
 #include "atom_vec.h"
 #include "comm.h"
@@ -333,8 +334,8 @@ The MPI standard requires that any MPI application calls
 do any MPI calls, MPI is still initialized internally to avoid errors
 accessing any MPI functions.  This function should then be called right
 before exiting the program to wait until all (parallel) tasks are
-completed and then MPI is cleanly shut down.  After this function no
+completed and then MPI is cleanly shut down.  After calling this
-more MPI calls may be made.
+function no more MPI calls may be made.
 
 .. versionadded:: 18Sep2020
 
@@ -353,6 +354,28 @@ void lammps_mpi_finalize()
   }
 }
 
+/* ---------------------------------------------------------------------- */
+
+/** Shut down the Kokkos library environment.
+ *
+\verbatim embed:rst
+
+The Kokkos library may only be initialized once during the execution of
+a process.  This is done automatically the first time Kokkos
+functionality is used.  This requires that the Kokkos environment
+must be explicitly shut down after any LAMMPS instance using it is
+closed (to release associated resources).
+After calling this function no Kokkos functionality may be used.
+
+.. versionadded:: TBD
+
+\endverbatim */
+
+void lammps_kokkos_finalize()
+{
+  KokkosLMP::finalize();
+}
+
 // ----------------------------------------------------------------------
 // Library functions to process commands
 // ----------------------------------------------------------------------

src/library.h

@@ -94,6 +94,7 @@ void lammps_close(void *handle);
 
 void lammps_mpi_init();
 void lammps_mpi_finalize();
+void lammps_kokkos_finalize();
 
 /* ----------------------------------------------------------------------
  * Library functions to process commands

src/main.cpp

@@ -14,6 +14,7 @@
 #include "lammps.h"
 
 #include "input.h"
+#include "accelerator_kokkos.h"
 #if defined(LAMMPS_EXCEPTIONS)
 #include "exceptions.h"
 #endif
@@ -77,13 +78,16 @@ int main(int argc, char **argv)
     lammps->input->file();
     delete lammps;
   } catch (LAMMPSAbortException &ae) {
+    KokkosLMP::finalize();
     MPI_Abort(ae.universe, 1);
   } catch (LAMMPSException &e) {
+    KokkosLMP::finalize();
     MPI_Barrier(lammps_comm);
     MPI_Finalize();
     exit(1);
   } catch (fmt::format_error &fe) {
     fprintf(stderr, "fmt::format_error: %s\n", fe.what());
+    KokkosLMP::finalize();
     MPI_Abort(MPI_COMM_WORLD, 1);
     exit(1);
   }
@@ -94,10 +98,12 @@ int main(int argc, char **argv)
     delete lammps;
   } catch (fmt::format_error &fe) {
     fprintf(stderr, "fmt::format_error: %s\n", fe.what());
+    KokkosLMP::finalize();
     MPI_Abort(MPI_COMM_WORLD, 1);
     exit(1);
   }
 #endif
+  KokkosLMP::finalize();
   MPI_Barrier(lammps_comm);
   MPI_Finalize();
 }

tools/swig/lammps.i

@@ -63,6 +63,7 @@ extern void  *lammps_open_fortran(int argc, char **argv, int f_comm);
 extern void   lammps_close(void *handle);
 extern void   lammps_mpi_init();
 extern void   lammps_mpi_finalize();
+extern void   lammps_kokkos_finalize();
 extern void   lammps_file(void *handle, const char *file);
 extern char  *lammps_command(void *handle, const char *cmd);
 extern void   lammps_commands_list(void *handle, int ncmd, const char **cmds);
@@ -185,6 +186,7 @@ extern void  *lammps_open_fortran(int argc, char **argv, int f_comm);
 extern void   lammps_close(void *handle);
 extern void   lammps_mpi_init();
 extern void   lammps_mpi_finalize();
+extern void   lammps_kokkos_finalize();
 extern void   lammps_file(void *handle, const char *file);
 extern char  *lammps_command(void *handle, const char *cmd);
 extern void   lammps_commands_list(void *handle, int ncmd, const char **cmds);

unittest/python/python-capabilities.py

@@ -165,6 +165,13 @@ class PythonCapabilities(unittest.TestCase):
             if self.cmake_cache['GPU_PREC'].lower() == 'single':
                  self.assertIn('single',settings['GPU']['precision'])
 
+        if self.cmake_cache['PKG_KOKKOS']:
+            if self.cmake_cache['Kokkos_ENABLE_OPENMP']:
+                self.assertIn('openmp',settings['KOKKOS']['api'])
+            if self.cmake_cache['Kokkos_ENABLE_SERIAL']:
+                self.assertIn('serial',settings['KOKKOS']['api'])
+            self.assertIn('double',settings['KOKKOS']['precision'])
+
     def test_gpu_device(self):
 
         info = self.lmp.get_gpu_device_info()