Fix issues in compute ave/sphere/atom

doc/src/compute_ave_sphere_atom.rst

@@ -35,7 +35,7 @@ Examples
 Description
 """""""""""
 
-Define a computation that calculates the local density and temperature
+Define a computation that calculates the local mass density and temperature
 for each atom and neighbors inside a spherical cutoff.
 
 The optional keyword *cutoff* defines the distance cutoff
@@ -58,7 +58,7 @@ too frequently.
    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
+   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
@@ -77,7 +77,7 @@ too frequently.
 Output info
 """""""""""
 
-This compute calculates a per-atom array with two columns: density and temperature.
+This compute calculates a per-atom array with two columns: mass 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

src/EXTRA-COMPUTE/compute_ave_sphere_atom.cpp

@@ -15,6 +15,7 @@
 
 #include "atom.h"
 #include "comm.h"
+#include "domain.h"
 #include "error.h"
 #include "force.h"
 #include "math_const.h"
@@ -98,7 +99,10 @@ void ComputeAveSphereAtom::init()
   }
 
   cutsq = cutoff * cutoff;
-  sphere_vol = 4.0 / 3.0 * MY_PI * cutsq * cutoff;
+  if (domain->dimension == 3)
+    volume = 4.0 / 3.0 * MY_PI * cutsq * cutoff;
+  else
+    volume = MY_PI * cutsq;
 
   // need an occasional full neighbor list
 
@@ -152,12 +156,25 @@ void ComputeAveSphereAtom::compute_peratom()
 
   double **x = atom->x;
   double **v = atom->v;
+  double *mass = atom->mass;
+  double *rmass = atom->rmass;
+  int *type = atom->type;
   int *mask = atom->mask;
+  double massone_i,massone_j;
+  double totalmass = 0.0;
+
+  double adof = domain->dimension;
+  double mvv2e = force->mvv2e;
+  double mv2d = force->mv2d;
+  double boltz = force->boltz;
 
   for (ii = 0; ii < inum; ii++) {
     i = ilist[ii];
 
     if (mask[i] & groupbit) {
+      if (rmass) massone_i = rmass[i];
+      else massone_i = mass[type[i]];
+
       xtmp = x[i][0];
       ytmp = x[i][1];
       ztmp = x[i][2];
@@ -174,6 +191,8 @@ void ComputeAveSphereAtom::compute_peratom()
       for (jj = 0; jj < jnum; jj++) {
         j = jlist[jj];
         j &= NEIGHMASK;
+        if (rmass) massone_j = rmass[i];
+        else massone_j = mass[type[i]];
 
         delx = xtmp - x[j][0];
         dely = ytmp - x[j][1];
@@ -194,10 +213,11 @@ void ComputeAveSphereAtom::compute_peratom()
       // i atom contribution
 
       count = 1;
+      totalmass = massone_i;
       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];
+      double ke_sum = massone_i * (vnet[0] * vnet[0] + vnet[1] * vnet[1] + vnet[2] * vnet[2]);
 
       for (jj = 0; jj < jnum; jj++) {
         j = jlist[jj];
@@ -209,14 +229,15 @@ void ComputeAveSphereAtom::compute_peratom()
         rsq = delx * delx + dely * dely + delz * delz;
         if (rsq < cutsq) {
           count++;
+          totalmass += massone_j;
           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];
+          ke_sum += massone_j * (vnet[0] * vnet[0] + vnet[1] * vnet[1] + vnet[2] * vnet[2]);
         }
       }
-      double density = count / sphere_vol;
+      double density = mv2d * totalmass / volume;
-      double temp = ke_sum / 3.0 / count;
+      double temp = mvv2e * ke_sum / (adof * count * boltz);
       result[i][0] = density;
       result[i][1] = temp;
     }

src/EXTRA-COMPUTE/compute_ave_sphere_atom.h

@@ -37,7 +37,7 @@ class ComputeAveSphereAtom : public Compute {
 
  protected:
   int nmax;
-  double cutoff, cutsq, sphere_vol;
+  double cutoff, cutsq, volume;
   class NeighList *list;
 
   double **result;

src/KOKKOS/compute_ave_sphere_atom_kokkos.cpp

@@ -16,6 +16,7 @@
 #include "atom_kokkos.h"
 #include "atom_masks.h"
 #include "comm.h"
+#include "domain.h"
 #include "error.h"
 #include "force.h"
 #include "memory_kokkos.h"
@@ -105,11 +106,19 @@ void ComputeAveSphereAtomKokkos<DeviceType>::compute_peratom()
   // 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);
+  atomKK->sync(execution_space,X_MASK|V_MASK|RMASS_MASK|TYPE_MASK|MASK_MASK);
   x = atomKK->k_x.view<DeviceType>();
   v = atomKK->k_v.view<DeviceType>();
+  rmass = atomKK->k_rmass.view<DeviceType>();
+  mass = atomKK->k_mass.view<DeviceType>();
+  type = atomKK->k_type.view<DeviceType>();
   mask = atomKK->k_mask.view<DeviceType>();
 
+  adof = domain->dimension;
+  mvv2e = force->mvv2e;
+  mv2d = force->mv2d;
+  boltz = force->boltz;
+
   Kokkos::deep_copy(d_result,0.0);
 
   copymode = 1;
@@ -125,8 +134,13 @@ template<class DeviceType>
 KOKKOS_INLINE_FUNCTION
 void ComputeAveSphereAtomKokkos<DeviceType>::operator()(TagComputeAveSphereAtom, const int &ii) const
 {
+  double massone_i,massone_j,totalmass;
+
   const int i = d_ilist[ii];
   if (mask[i] & groupbit) {
+    if (rmass.data()) massone_i = rmass[i];
+    else massone_i = mass[type[i]];
+
     const X_FLOAT xtmp = x(i,0);
     const X_FLOAT ytmp = x(i,1);
     const X_FLOAT ztmp = x(i,2);
@@ -164,15 +178,18 @@ void ComputeAveSphereAtomKokkos<DeviceType>::operator()(TagComputeAveSphereAtom,
     // i atom contribution
 
     count = 1;
+    totalmass = massone_i;
     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];
+    double ke_sum = massone_i * (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;
+      if (rmass.data()) massone_j = rmass[i];
+      else massone_j = mass[type[i]];
 
       const F_FLOAT delx = x(j,0) - xtmp;
       const F_FLOAT dely = x(j,1) - ytmp;
@@ -180,14 +197,15 @@ void ComputeAveSphereAtomKokkos<DeviceType>::operator()(TagComputeAveSphereAtom,
       const F_FLOAT rsq = delx*delx + dely*dely + delz*delz;
       if (rsq < cutsq) {
         count++;
+        totalmass += massone_j;
         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];
+        ke_sum += massone_j * (vnet[0]*vnet[0] + vnet[1]*vnet[1] + vnet[2]*vnet[2]);
       }
     }
-    double density = count/sphere_vol;
+    double density = mv2d*totalmass/volume;
-    double temp = ke_sum/3.0/count;
+    double temp = mvv2e*ke_sum/(adof*count*boltz);
     d_result(i,0) = density;
     d_result(i,1) = temp;
   }

src/KOKKOS/compute_ave_sphere_atom_kokkos.h

@@ -46,13 +46,18 @@ template <class DeviceType> class ComputeAveSphereAtomKokkos : public ComputeAve
   void operator()(TagComputeAveSphereAtom, const int &) const;
 
  private:
-  typename AT::t_x_array_randomread x;
+  double adof,mvv2e,mv2d,boltz;
-  typename AT::t_v_array_randomread v;
+
+  typename AT::t_x_array x;
+  typename AT::t_v_array v;
+  typename ArrayTypes<DeviceType>::t_float_1d rmass;
+  typename ArrayTypes<DeviceType>::t_float_1d mass;
+  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_randomread d_ilist;
+  typename AT::t_int_1d d_ilist;
-  typename AT::t_int_1d_randomread d_numneigh;
+  typename AT::t_int_1d d_numneigh;
 
   DAT::tdual_float_2d k_result;
   typename AT::t_float_2d d_result;