diff --git a/doc/src/compute_ave_sphere_atom.rst b/doc/src/compute_ave_sphere_atom.rst index db04682865..0cf631a941 100644 --- a/doc/src/compute_ave_sphere_atom.rst +++ b/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 ` 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 ` 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 ` doc diff --git a/src/EXTRA-COMPUTE/compute_ave_sphere_atom.cpp b/src/EXTRA-COMPUTE/compute_ave_sphere_atom.cpp index 07803aca20..b63c5c2b07 100644 --- a/src/EXTRA-COMPUTE/compute_ave_sphere_atom.cpp +++ b/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 temp = ke_sum / 3.0 / count; + double density = mv2d * totalmass / volume; + double temp = mvv2e * ke_sum / (adof * count * boltz); result[i][0] = density; result[i][1] = temp; } diff --git a/src/EXTRA-COMPUTE/compute_ave_sphere_atom.h b/src/EXTRA-COMPUTE/compute_ave_sphere_atom.h index ffed09bae5..76350997f9 100644 --- a/src/EXTRA-COMPUTE/compute_ave_sphere_atom.h +++ b/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; diff --git a/src/KOKKOS/compute_ave_sphere_atom_kokkos.cpp b/src/KOKKOS/compute_ave_sphere_atom_kokkos.cpp index d2cb6682a7..df70a27738 100644 --- a/src/KOKKOS/compute_ave_sphere_atom_kokkos.cpp +++ b/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::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(); v = atomKK->k_v.view(); + rmass = atomKK->k_rmass.view(); + mass = atomKK->k_mass.view(); + type = atomKK->k_type.view(); mask = atomKK->k_mask.view(); + 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 KOKKOS_INLINE_FUNCTION void ComputeAveSphereAtomKokkos::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::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::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 temp = ke_sum/3.0/count; + double density = mv2d*totalmass/volume; + double temp = mvv2e*ke_sum/(adof*count*boltz); d_result(i,0) = density; d_result(i,1) = temp; } diff --git a/src/KOKKOS/compute_ave_sphere_atom_kokkos.h b/src/KOKKOS/compute_ave_sphere_atom_kokkos.h index 75b5ca3aba..1ddf943e7d 100644 --- a/src/KOKKOS/compute_ave_sphere_atom_kokkos.h +++ b/src/KOKKOS/compute_ave_sphere_atom_kokkos.h @@ -46,13 +46,18 @@ template class ComputeAveSphereAtomKokkos : public ComputeAve void operator()(TagComputeAveSphereAtom, const int &) const; private: - typename AT::t_x_array_randomread x; - typename AT::t_v_array_randomread v; + double adof,mvv2e,mv2d,boltz; + + typename AT::t_x_array x; + typename AT::t_v_array v; + typename ArrayTypes::t_float_1d rmass; + typename ArrayTypes::t_float_1d mass; + typename ArrayTypes::t_int_1d type; typename ArrayTypes::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; + 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;