diff --git a/doc/src/Commands_compute.rst b/doc/src/Commands_compute.rst index 2c67d44f24..91910ccdaa 100644 --- a/doc/src/Commands_compute.rst +++ b/doc/src/Commands_compute.rst @@ -28,6 +28,7 @@ KOKKOS, o = OPENMP, t = OPT. * :doc:`angle ` * :doc:`angle/local ` * :doc:`angmom/chunk ` + * :doc:`ave/sphere/atom (k) ` * :doc:`basal/atom ` * :doc:`body/local ` * :doc:`bond ` diff --git a/doc/src/Errors_messages.rst b/doc/src/Errors_messages.rst index 3a593b5a3f..c06f4c86e3 100644 --- a/doc/src/Errors_messages.rst +++ b/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. diff --git a/doc/src/compute.rst b/doc/src/compute.rst index 71d3bada76..9edd7e8474 100644 --- a/doc/src/compute.rst +++ b/doc/src/compute.rst @@ -174,6 +174,7 @@ The individual style names on the :doc:`Commands compute ` pag * :doc:`angle ` - energy of each angle sub-style * :doc:`angle/local ` - theta and energy of each angle * :doc:`angmom/chunk ` - angular momentum for each chunk +* :doc:`ave/sphere/atom ` - compute local density and temperature around each atom * :doc:`basal/atom ` - calculates the hexagonal close-packed "c" lattice vector of each atom * :doc:`body/local ` - attributes of body sub-particles * :doc:`bond ` - energy of each bond sub-style diff --git a/doc/src/compute_ave_sphere_atom.rst b/doc/src/compute_ave_sphere_atom.rst new file mode 100644 index 0000000000..db04682865 --- /dev/null +++ b/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 ` 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 ` command can remove pairwise + 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 + 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 + snapshots in the dump file. The rerun script can use a + :doc:`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 ` 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 ` page for more info. + +Related commands +"""""""""""""""" + +:doc:`comm_modify ` + +Default +""""""" + +The option defaults are *cutoff* = pair style cutoff + diff --git a/doc/src/dump_modify.rst b/doc/src/dump_modify.rst index 0094fceb7b..be75153f6f 100644 --- a/doc/src/dump_modify.rst +++ b/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:: diff --git a/src/.gitignore b/src/.gitignore index 19bafa4b52..695c9a19af 100644 --- a/src/.gitignore +++ b/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 diff --git a/src/Depend.sh b/src/Depend.sh index af88f24bb4..a8e17e0546 100755 --- a/src/Depend.sh +++ b/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 diff --git a/src/EXTRA-COMPUTE/compute_ave_sphere_atom.cpp b/src/EXTRA-COMPUTE/compute_ave_sphere_atom.cpp new file mode 100644 index 0000000000..14a4c364a1 --- /dev/null +++ b/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 +#include + +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; +} diff --git a/src/EXTRA-COMPUTE/compute_ave_sphere_atom.h b/src/EXTRA-COMPUTE/compute_ave_sphere_atom.h new file mode 100644 index 0000000000..9b5e38750b --- /dev/null +++ b/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. + +*/ diff --git a/src/H5MD/dump_h5md.cpp b/src/H5MD/dump_h5md.cpp index bc9c98caa0..a59de9773c 100644 --- a/src/H5MD/dump_h5md.cpp +++ b/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 (n); if (every_position>=0) diff --git a/src/H5MD/dump_h5md.h b/src/H5MD/dump_h5md.h index 28b74649c6..5e3f3f8279 100644 --- a/src/H5MD/dump_h5md.h +++ b/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. diff --git a/src/KOKKOS/Install.sh b/src/KOKKOS/Install.sh index 04bf84ed31..45fa0654a9 100755 --- a/src/KOKKOS/Install.sh +++ b/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 diff --git a/src/KOKKOS/compute_ave_sphere_atom_kokkos.cpp b/src/KOKKOS/compute_ave_sphere_atom_kokkos.cpp new file mode 100644 index 0000000000..3f83c24fb6 --- /dev/null +++ b/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 +#include + +using namespace LAMMPS_NS; +using namespace MathConst; + + +/* ---------------------------------------------------------------------- */ + +template +ComputeAveSphereAtomKokkos::ComputeAveSphereAtomKokkos(LAMMPS *lmp, int narg, char **arg) : + ComputeAveSphereAtom(lmp, narg, arg) +{ + kokkosable = 1; + atomKK = (AtomKokkos *) atom; + execution_space = ExecutionSpaceFromDevice::space; + datamask_read = EMPTY_MASK; + datamask_modify = EMPTY_MASK; +} + +/* ---------------------------------------------------------------------- */ + +template +ComputeAveSphereAtomKokkos::~ComputeAveSphereAtomKokkos() +{ + if (copymode) return; + + memoryKK->destroy_kokkos(k_result,result); +} + +/* ---------------------------------------------------------------------- */ + +template +void ComputeAveSphereAtomKokkos::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::value && + !std::is_same::value; + neighbor->requests[irequest]-> + kokkos_device = std::is_same::value; +} + +/* ---------------------------------------------------------------------- */ + +template +void ComputeAveSphereAtomKokkos::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(); + 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* k_list = static_cast*>(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(); + v = atomKK->k_v.view(); + mask = atomKK->k_mask.view(); + + Kokkos::deep_copy(d_result,0.0); + + copymode = 1; + typename Kokkos::RangePolicy policy(0,inum); + Kokkos::parallel_for("ComputeAveSphereAtom",policy,*this); + copymode = 0; + + k_result.modify(); + k_result.sync_host(); +} + +template +KOKKOS_INLINE_FUNCTION +void ComputeAveSphereAtomKokkos::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; +#ifdef LMP_KOKKOS_GPU +template class ComputeAveSphereAtomKokkos; +#endif +} diff --git a/src/KOKKOS/compute_ave_sphere_atom_kokkos.h b/src/KOKKOS/compute_ave_sphere_atom_kokkos.h new file mode 100644 index 0000000000..42607e5239 --- /dev/null +++ b/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) +ComputeStyle(ave/sphere/atom/kk/device,ComputeAveSphereAtomKokkos) +ComputeStyle(ave/sphere/atom/kk/host,ComputeAveSphereAtomKokkos) + +#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 ComputeAveSphereAtomKokkos : public ComputeAveSphereAtom { + public: + typedef DeviceType device_type; + typedef ArrayTypes 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::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: + +*/ diff --git a/src/dump.cpp b/src/dump.cpp index 0f338e99b5..e02116afab 100644 --- a/src/dump.cpp +++ b/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 (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 (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 diff --git a/src/dump.h b/src/dump.h index caeef827c7..35da154d7c 100644 --- a/src/dump.h +++ b/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