Merge branch 'develop' into time-dumps2

doc/src/Commands_compute.rst

@@ -28,6 +28,7 @@ KOKKOS, o = OPENMP, t = OPT.
    * :doc:`angle <compute_angle>`
    * :doc:`angle/local <compute_angle_local>`
    * :doc:`angmom/chunk <compute_angmom_chunk>`
+   * :doc:`ave/sphere/atom (k) <compute_ave_sphere_atom>`
    * :doc:`basal/atom <compute_basal_atom>`
    * :doc:`body/local <compute_body_local>`
    * :doc:`bond <compute_bond>`

doc/src/Errors_messages.rst

@@ -7772,9 +7772,6 @@ keyword to allow for additional bonds to be formed
    The system size must fit in a 32-bit integer to use this dump
    style.
 
-*Too many atoms to dump sort*
-   Cannot sort when running with more than 2\^31 atoms.
-
 *Too many elements extracted from MEAM library.*
    Increase 'maxelt' in meam.h and recompile.
 

doc/src/compute.rst

@@ -174,6 +174,7 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` pag
 * :doc:`angle <compute_angle>` - energy of each angle sub-style
 * :doc:`angle/local <compute_angle_local>` - theta and energy of each angle
 * :doc:`angmom/chunk <compute_angmom_chunk>` - angular momentum for each chunk
+* :doc:`ave/sphere/atom <compute_ave_sphere_atom>` - compute local density and temperature around each atom
 * :doc:`basal/atom <compute_basal_atom>` - calculates the hexagonal close-packed "c" lattice vector of each atom
 * :doc:`body/local <compute_body_local>` - attributes of body sub-particles
 * :doc:`bond <compute_bond>` - energy of each bond sub-style

doc/src/compute_ave_sphere_atom.rst

@@ -0,0 +1,101 @@
+.. index:: compute ave/sphere/atom
+.. index:: compute ave/sphere/atom/kk
+
+compute ave/sphere/atom command
+================================
+
+Accelerator Variants: *ave/sphere/atom/kk*
+
+Syntax
+""""""
+
+.. parsed-literal::
+
+   compute ID group-ID ave/sphere/atom keyword values ...
+
+* ID, group-ID are documented in :doc:`compute <compute>` command
+* ave/sphere/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 ave/sphere/atom
+
+   compute 1 all ave/sphere/atom cutoff 5.0
+   comm_modify cutoff 5.0
+
+Description
+"""""""""""
+
+Define a computation that calculates the local density and temperature
+for each atom and neighbors inside a spherical cutoff.
+
+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 fix 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.
+
+----------
+
+
+.. include:: accel_styles.rst
+
+
+----------
+
+Output info
+"""""""""""
+
+This compute calculates a per-atom array with two columns: density and temperature.
+
+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/src/dump_modify.rst

@@ -662,7 +662,9 @@ The dump *local* style cannot be sorted by atom ID, since there are
 typically multiple lines of output per atom.  Some dump styles, such
 as *dcd* and *xtc*, require sorting by atom ID to format the output
 file correctly.  If multiple processors are writing the dump file, via
-the "%" wildcard in the dump filename, then sorting cannot be
+the "%" wildcard in the dump filename and the *nfile* or *fileper*
+keywords are set to non-default values (i.e. the number of dump file
+pieces is not equal to the number of procs), then sorting cannot be
 performed.
 
 .. note::

src/.gitignore

@@ -429,6 +429,8 @@
 /commgrid.h
 /compute_ackland_atom.cpp
 /compute_ackland_atom.h
+/compute_ave_sphere_atom.cpp
+/compute_ave_sphere_atom.h
 /compute_basal_atom.cpp
 /compute_basal_atom.h
 /compute_body_local.cpp

src/Depend.sh

@@ -77,6 +77,10 @@ if (test $1 = "DPD-BASIC") then
   depend INTEL
 fi
 
+if (test $1 = "EXTRA-COMPUTE") then
+  depend KOKKOS
+fi
+
 if (test $1 = "EXTRA-MOLECULE") then
   depend GPU
   depend OPENMP

src/EXTRA-COMPUTE/compute_ave_sphere_atom.cpp

@@ -0,0 +1,278 @@
+/* ----------------------------------------------------------------------
+   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_ave_sphere_atom.h"
+
+#include "atom.h"
+#include "comm.h"
+#include "error.h"
+#include "force.h"
+#include "group.h"
+#include "memory.h"
+#include "modify.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "neighbor.h"
+#include "pair.h"
+#include "update.h"
+#include "math_const.h"
+
+#include <cmath>
+#include <cstring>
+
+using namespace LAMMPS_NS;
+using namespace MathConst;
+
+
+/* ---------------------------------------------------------------------- */
+
+ComputeAveSphereAtom::ComputeAveSphereAtom(LAMMPS *lmp, int narg, char **arg) :
+  Compute(lmp, narg, arg),
+  result(nullptr)
+{
+  if (narg < 3 || narg > 5) error->all(FLERR,"Illegal compute ave/sphere/atom command");
+
+  // process optional args
+
+  cutoff = 0.0;
+
+  int iarg = 3;
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"cutoff") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal compute ave/sphere/atom command");
+      cutoff = utils::numeric(FLERR,arg[iarg+1],false,lmp);
+      if (cutoff <= 0.0) error->all(FLERR,"Illegal compute ave/sphere/atom command");
+      iarg += 2;
+    } else error->all(FLERR,"Illegal compute ave/sphere/atom command");
+  }
+
+  peratom_flag = 1;
+  size_peratom_cols = 2;
+  comm_forward = 3;
+
+  nmax = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeAveSphereAtom::~ComputeAveSphereAtom()
+{
+  if (copymode) return;
+
+  memory->destroy(result);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeAveSphereAtom::init()
+{
+  if (!force->pair && cutoff == 0.0)
+    error->all(FLERR,"Compute ave/sphere/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 ave/sphere/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;
+  sphere_vol = 4.0/3.0*MY_PI*cutsq*cutoff;
+
+  // need an occasional full neighbor list
+
+  int irequest = neighbor->request(this,instance_me);
+  neighbor->requests[irequest]->pair = 0;
+  neighbor->requests[irequest]->compute = 1;
+  neighbor->requests[irequest]->half = 0;
+  neighbor->requests[irequest]->full = 1;
+  neighbor->requests[irequest]->occasional = 1;
+  if (cutflag) {
+    neighbor->requests[irequest]->cut = 1;
+    neighbor->requests[irequest]->cutoff = cutoff;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeAveSphereAtom::init_list(int /*id*/, NeighList *ptr)
+{
+  list = ptr;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeAveSphereAtom::compute_peratom()
+{
+  int i,j,ii,jj,inum,jnum;
+  double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
+  int *ilist,*jlist,*numneigh,**firstneigh;
+  int count;
+  double vsum[3],vavg[3],vnet[3];
+
+  invoked_peratom = update->ntimestep;
+
+  // grow result array if necessary
+
+  if (atom->nmax > nmax) {
+    memory->destroy(result);
+    nmax = atom->nmax;
+    memory->create(result,nmax,2,"ave/sphere/atom:result");
+    array_atom = result;
+  }
+
+  // need velocities of ghost atoms
+
+  comm->forward_comm_compute(this);
+
+  // 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;
+  double **v = atom->v;
+  int *mask = atom->mask;
+
+  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;
+      vsum[0] = v[i][0];
+      vsum[1] = v[i][1];
+      vsum[2] = v[i][2];
+
+      for (jj = 0; jj < jnum; jj++) {
+        j = jlist[jj];
+        j &= NEIGHMASK;
+
+        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++;
+          vsum[0] += v[j][0];
+          vsum[1] += v[j][1];
+          vsum[2] += v[j][2];
+        }
+      }
+
+      vavg[0] = vsum[0]/count;
+      vavg[1] = vsum[1]/count;
+      vavg[2] = vsum[2]/count;
+
+      // i atom contribution
+
+      count = 1;
+      vnet[0] = v[i][0] - vavg[0];
+      vnet[1] = v[i][1] - vavg[1];
+      vnet[2] = v[i][2] - vavg[2];
+      double ke_sum = vnet[0]*vnet[0] + vnet[1]*vnet[1] + vnet[2]*vnet[2];
+
+      for (jj = 0; jj < jnum; jj++) {
+        j = jlist[jj];
+        j &= NEIGHMASK;
+
+        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++;
+          vnet[0] = v[j][0] - vavg[0];
+          vnet[1] = v[j][1] - vavg[1];
+          vnet[2] = v[j][2] - vavg[2];
+          ke_sum += vnet[0]*vnet[0] + vnet[1]*vnet[1] + vnet[2]*vnet[2];
+        }
+      }
+      double density = count/sphere_vol;
+      double temp = ke_sum/3.0/count;
+      result[i][0] = density;
+      result[i][1] = temp;
+    }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+int ComputeAveSphereAtom::pack_forward_comm(int n, int *list, double *buf,
+                                        int /*pbc_flag*/, int * /*pbc*/)
+{
+  double **v = atom->v;
+
+  int i,m=0;
+  for (i = 0; i < n; ++i) {
+    buf[m++] = v[list[i]][0];
+    buf[m++] = v[list[i]][1];
+    buf[m++] = v[list[i]][2];
+  }
+
+  return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeAveSphereAtom::unpack_forward_comm(int n, int first, double *buf)
+{
+  double **v = atom->v;
+
+  int i,last,m=0;
+  last = first + n;
+  for (i = first; i < last; ++i) {
+    v[i][0] = buf[m++];
+    v[i][1] = buf[m++];
+    v[i][2] = buf[m++];
+  }
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of local atom-based array
+------------------------------------------------------------------------- */
+
+double ComputeAveSphereAtom::memory_usage()
+{
+  double bytes = (double)2*nmax * sizeof(double);
+  return bytes;
+}

src/EXTRA-COMPUTE/compute_ave_sphere_atom.h

@@ -0,0 +1,67 @@
+/* -*- 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
+
+ComputeStyle(ave/sphere/atom,ComputeAveSphereAtom)
+
+#else
+
+#ifndef LMP_COMPUTE_AVE_SPHERE_ATOM_H
+#define LMP_COMPUTE_AVE_SPHERE_ATOM_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeAveSphereAtom : public Compute {
+ public:
+  ComputeAveSphereAtom(class LAMMPS *, int, char **);
+  virtual ~ComputeAveSphereAtom();
+  virtual void init();
+  void init_list(int, class NeighList *);
+  virtual void compute_peratom();
+  int pack_forward_comm(int, int *, double *, int, int *);
+  void unpack_forward_comm(int, int, double *);
+  double memory_usage();
+
+ protected:
+  int nmax;
+  double cutoff,cutsq,sphere_vol;
+  class NeighList *list;
+
+  double **result;
+};
+
+}
+
+#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 ave/sphere/atom requires a cutoff be specified or a pair style be defined
+
+Self-explanatory.
+
+E: Compute ave/sphere/atom cutoff exceeds ghost atom range - use comm_modify cutoff command
+
+Self-explanatory.
+
+*/

src/H5MD/dump_h5md.cpp

@@ -181,6 +181,7 @@ DumpH5MD::DumpH5MD(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, arg)
   // allocate global array for atom coords
 
   bigint n = group->count(igroup);
+  if ((bigint) domain->dimension*n > MAXSMALLINT) error->all(FLERR,"Too many atoms for dump h5md");
   natoms = static_cast<int> (n);
 
   if (every_position>=0)

src/H5MD/dump_h5md.h

@@ -90,6 +90,11 @@ E: Dump h5md requires sorting by atom ID
 
 Use the dump_modify sort command to enable this.
 
+E: Too many atoms for dump h5md
+
+The system size must fit in a 32-bit integer to use this dump
+style.
+
 E: Cannot use variable every setting for dump xtc
 
 The format of this file requires snapshots at regular intervals.

src/KOKKOS/Install.sh

@@ -88,6 +88,8 @@ action comm_kokkos.cpp
 action comm_kokkos.h
 action comm_tiled_kokkos.cpp
 action comm_tiled_kokkos.h
+action compute_ave_sphere_atom_kokkos.cpp compute_ave_sphere_atom.cpp
+action compute_ave_sphere_atom_kokkos.h compute_ave_sphere_atom.h
 action compute_coord_atom_kokkos.cpp
 action compute_coord_atom_kokkos.h
 action compute_orientorder_atom_kokkos.cpp

src/KOKKOS/compute_ave_sphere_atom_kokkos.cpp

@@ -0,0 +1,209 @@
+/* ----------------------------------------------------------------------
+   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_ave_sphere_atom_kokkos.h"
+
+#include "atom_kokkos.h"
+#include "atom_masks.h"
+#include "comm.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 "math_const.h"
+
+#include <cmath>
+#include <cstring>
+
+using namespace LAMMPS_NS;
+using namespace MathConst;
+
+
+/* ---------------------------------------------------------------------- */
+
+template<class DeviceType>
+ComputeAveSphereAtomKokkos<DeviceType>::ComputeAveSphereAtomKokkos(LAMMPS *lmp, int narg, char **arg) :
+  ComputeAveSphereAtom(lmp, narg, arg)
+{
+  kokkosable = 1;
+  atomKK = (AtomKokkos *) atom;
+  execution_space = ExecutionSpaceFromDevice<DeviceType>::space;
+  datamask_read = EMPTY_MASK;
+  datamask_modify = EMPTY_MASK;
+}
+
+/* ---------------------------------------------------------------------- */
+
+template<class DeviceType>
+ComputeAveSphereAtomKokkos<DeviceType>::~ComputeAveSphereAtomKokkos()
+{
+  if (copymode) return;
+
+  memoryKK->destroy_kokkos(k_result,result);
+}
+
+/* ---------------------------------------------------------------------- */
+
+template<class DeviceType>
+void ComputeAveSphereAtomKokkos<DeviceType>::init()
+{
+  ComputeAveSphereAtom::init();
+
+  // need an occasional full neighbor list
+
+  // irequest = neigh request made by parent class
+
+  int irequest = neighbor->nrequest - 1;
+
+  neighbor->requests[irequest]->
+    kokkos_host = std::is_same<DeviceType,LMPHostType>::value &&
+    !std::is_same<DeviceType,LMPDeviceType>::value;
+  neighbor->requests[irequest]->
+    kokkos_device = std::is_same<DeviceType,LMPDeviceType>::value;
+}
+
+/* ---------------------------------------------------------------------- */
+
+template<class DeviceType>
+void ComputeAveSphereAtomKokkos<DeviceType>::compute_peratom()
+{
+  invoked_peratom = update->ntimestep;
+
+  // grow result array if necessary
+
+  if (atom->nmax > nmax) {
+    memoryKK->destroy_kokkos(k_result,result);
+    nmax = atom->nmax;
+    memoryKK->create_kokkos(k_result,result,nmax,2,"ave/sphere/atom:result");
+    d_result = k_result.view<DeviceType>();
+    array_atom = result;
+  }
+
+  // need velocities of ghost atoms
+
+  atomKK->sync(Host,V_MASK);
+  comm->forward_comm_compute(this);
+  atomKK->modified(Host,V_MASK);
+
+  // 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|V_MASK|TYPE_MASK|MASK_MASK);
+  x = atomKK->k_x.view<DeviceType>();
+  v = atomKK->k_v.view<DeviceType>();
+  mask = atomKK->k_mask.view<DeviceType>();
+
+  Kokkos::deep_copy(d_result,0.0);
+
+  copymode = 1;
+  typename Kokkos::RangePolicy<DeviceType, TagComputeAveSphereAtom> policy(0,inum);
+  Kokkos::parallel_for("ComputeAveSphereAtom",policy,*this);
+  copymode = 0;
+
+  k_result.modify<DeviceType>();
+  k_result.sync_host();
+}
+
+template<class DeviceType>
+KOKKOS_INLINE_FUNCTION
+void ComputeAveSphereAtomKokkos<DeviceType>::operator()(TagComputeAveSphereAtom, 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;
+    double vsum[3];
+    vsum[0] = v(i,0);
+    vsum[1] = v(i,1);
+    vsum[2] = v(i,2);
+
+    for (int jj = 0; jj < jnum; jj++) {
+      int j = d_neighbors(i,jj);
+      j &= NEIGHMASK;
+
+      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++;
+        vsum[0] += v(j,0);
+        vsum[1] += v(j,1);
+        vsum[2] += v(j,2);
+      }
+    }
+
+    double vavg[3];
+    vavg[0] = vsum[0]/count;
+    vavg[1] = vsum[1]/count;
+    vavg[2] = vsum[2]/count;
+
+    // i atom contribution
+
+    count = 1;
+    double vnet[3];
+    vnet[0] = v(i,0) - vavg[0];
+    vnet[1] = v(i,1) - vavg[1];
+    vnet[2] = v(i,2) - vavg[2];
+    double ke_sum = vnet[0]*vnet[0] + vnet[1]*vnet[1] + vnet[2]*vnet[2];
+
+    for (int jj = 0; jj < jnum; jj++) {
+      int j = d_neighbors(i,jj);
+      j &= NEIGHMASK;
+
+      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++;
+        vnet[0] = v(j,0) - vavg[0];
+        vnet[1] = v(j,1) - vavg[1];
+        vnet[2] = v(j,2) - vavg[2];
+        ke_sum += vnet[0]*vnet[0] + vnet[1]*vnet[1] + vnet[2]*vnet[2];
+      }
+    }
+    double density = count/sphere_vol;
+    double temp = ke_sum/3.0/count;
+    d_result(i,0) = density;
+    d_result(i,1) = temp;
+  }
+}
+
+namespace LAMMPS_NS {
+template class ComputeAveSphereAtomKokkos<LMPDeviceType>;
+#ifdef LMP_KOKKOS_GPU
+template class ComputeAveSphereAtomKokkos<LMPHostType>;
+#endif
+}

src/KOKKOS/compute_ave_sphere_atom_kokkos.h

@@ -0,0 +1,66 @@
+/* -*- 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
+
+ComputeStyle(ave/sphere/atom/kk,ComputeAveSphereAtomKokkos<LMPDeviceType>)
+ComputeStyle(ave/sphere/atom/kk/device,ComputeAveSphereAtomKokkos<LMPDeviceType>)
+ComputeStyle(ave/sphere/atom/kk/host,ComputeAveSphereAtomKokkos<LMPHostType>)
+
+#else
+
+#ifndef LMP_COMPUTE_AVE_SPHERE_ATOM_KOKKOS_H
+#define LMP_COMPUTE_AVE_SPHERE_ATOM_KOKKOS_H
+
+#include "compute_ave_sphere_atom.h"
+#include "kokkos_type.h"
+
+namespace LAMMPS_NS {
+
+struct TagComputeAveSphereAtom{};
+
+template<class DeviceType>
+class ComputeAveSphereAtomKokkos : public ComputeAveSphereAtom {
+ public:
+  typedef DeviceType device_type;
+  typedef ArrayTypes<DeviceType> AT;
+
+  ComputeAveSphereAtomKokkos(class LAMMPS *, int, char **);
+  virtual ~ComputeAveSphereAtomKokkos();
+  void init();
+  void compute_peratom();
+
+  KOKKOS_INLINE_FUNCTION
+  void operator()(TagComputeAveSphereAtom, const int&) const;
+
+ private:
+  typename AT::t_x_array_randomread x;
+  typename AT::t_v_array_randomread v;
+  typename ArrayTypes<DeviceType>::t_int_1d mask;
+
+  typename AT::t_neighbors_2d d_neighbors;
+  typename AT::t_int_1d_randomread d_ilist;
+  typename AT::t_int_1d_randomread d_numneigh;
+
+  DAT::tdual_float_2d k_result;
+  typename AT::t_float_2d d_result;
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+*/

src/dump.cpp

@@ -228,7 +228,7 @@ void Dump::init()
   if (sort_flag) {
     if (multiproc > 1)
       error->all(FLERR,
-                 "Cannot dump sort when multiple dump files are written");
+                 "Cannot dump sort when 'nfile' or 'fileper' keywords are set to non-default values");
     if (sortcol == 0 && atom->tag_enable == 0)
       error->all(FLERR,"Cannot dump sort on atom IDs with no atom IDs defined");
     if (sortcol && sortcol > size_one)
@@ -237,8 +237,6 @@ void Dump::init()
       irregular = new Irregular(lmp);
 
     bigint size = group->count(igroup);
-    if (size > MAXSMALLINT) error->all(FLERR,"Too many atoms to dump sort");
-    int isize = static_cast<int> (size);
 
     // set reorderflag = 1 if can simply reorder local atoms rather than sort
     // criteria: sorting by ID, atom IDs are consecutive from 1 to Natoms
@@ -268,7 +266,7 @@ void Dump::init()
       MPI_Allreduce(&min,&minall,1,MPI_LMP_TAGINT,MPI_MIN,world);
       MPI_Allreduce(&max,&maxall,1,MPI_LMP_TAGINT,MPI_MAX,world);
 
-      if (maxall-minall+1 == isize) {
+      if (maxall-minall+1 == size) {
         reorderflag = 1;
         double range = maxall-minall + EPSILON;
         idlo = static_cast<tagint> (range*me/nprocs + minall);
@@ -284,7 +282,7 @@ void Dump::init()
         else if (me+1 != hi) idhi++;
 
         nme_reorder = idhi-idlo;
-        ntotal_reorder = isize;
+        ntotal_reorder = size;
       }
     }
   }
@@ -369,16 +367,6 @@ void Dump::write()
   if (multiproc != nprocs) MPI_Allreduce(&nme,&nmax,1,MPI_INT,MPI_MAX,world);
   else nmax = nme;
 
-  // write timestep header
-  // for multiproc,
-  //   nheader = # of lines in this file via Allreduce on clustercomm
-
-  bigint nheader = ntotal;
-  if (multiproc)
-    MPI_Allreduce(&bnme,&nheader,1,MPI_LMP_BIGINT,MPI_SUM,clustercomm);
-
-  if (filewriter && write_header_flag) write_header(nheader);
-
   // insure buf is sized for packing and communicating
   // use nmax to insure filewriter proc can receive info from others
   // limit nmax*size_one to int since used as arg in MPI calls
@@ -431,6 +419,19 @@ void Dump::write()
   else pack(nullptr);
   if (sort_flag) sort();
 
+  // write timestep header
+  // for multiproc,
+  //   nheader = # of lines in this file via Allreduce on clustercomm
+  //   must come after sort, which can change nme
+
+  bigint nheader = ntotal;
+  if (multiproc) {
+    bnme = nme;
+    MPI_Allreduce(&bnme,&nheader,1,MPI_LMP_BIGINT,MPI_SUM,clustercomm);
+  }
+
+  if (filewriter && write_header_flag) write_header(nheader);
+
   // if buffering, convert doubles into strings
   // insure sbuf is sized for communicating
   // cannot buffer if output is to binary file

src/dump.h

@@ -116,7 +116,7 @@ class Dump : protected Pointers {
 
   bigint ntotal;         // total # of per-atom lines in snapshot
   int reorderflag;       // 1 if OK to reorder instead of sort
-  int ntotal_reorder;    // # of atoms that must be in snapshot
+  bigint ntotal_reorder;    // # of atoms that must be in snapshot
   int nme_reorder;       // # of atoms I must own in snapshot
   tagint idlo;           // lowest ID I own when reordering
 
@@ -173,10 +173,9 @@ E: Dump file MPI-IO output not allowed with % in filename
 This is because a % signifies one file per processor and MPI-IO
 creates one large file for all processors.
 
-E: Cannot dump sort when multiple dump files are written
+E: Cannot dump sort when 'nfile' or 'fileper' keywords are set to non-default values
 
-In this mode, each processor dumps its atoms to a file, so
-no sorting is allowed.
+Can only dump sort when the number of dump file pieces using % in filename equals the number of processors
 
 E: Cannot dump sort on atom IDs with no atom IDs defined
 
@@ -186,10 +185,6 @@ E: Dump sort column is invalid
 
 Self-explanatory.
 
-E: Too many atoms to dump sort
-
-Cannot sort when running with more than 2^31 atoms.
-
 E: Dump could not find refresh compute ID
 
 UNDOCUMENTED