diff --git a/doc/src/Commands_compute.rst b/doc/src/Commands_compute.rst index 7c73583a4f..b53d9d6820 100644 --- a/doc/src/Commands_compute.rst +++ b/doc/src/Commands_compute.rst @@ -178,6 +178,7 @@ KOKKOS, o = OPENMP, t = OPT. * :doc:`ti ` * :doc:`torque/chunk ` * :doc:`vacf ` + * :doc:`vacf/chunk ` * :doc:`vcm/chunk ` * :doc:`viscosity/cos ` * :doc:`voronoi/atom ` diff --git a/doc/src/compute.rst b/doc/src/compute.rst index 9a8a1734fb..323f3ed3e5 100644 --- a/doc/src/compute.rst +++ b/doc/src/compute.rst @@ -356,6 +356,7 @@ The individual style names on the :doc:`Commands compute ` pag * :doc:`ti ` - thermodynamic integration free energy values * :doc:`torque/chunk ` - torque applied on each chunk * :doc:`vacf ` - velocity auto-correlation function of group of atoms +* :doc:`vacf/chunk ` - velocity auto-correlation for the center of mass velocities of chunks of atoms * :doc:`vcm/chunk ` - velocity of center-of-mass for each chunk * :doc:`viscosity/cos ` - velocity profile under cosine-shaped acceleration * :doc:`voronoi/atom ` - Voronoi volume and neighbors for each atom diff --git a/doc/src/compute_msd.rst b/doc/src/compute_msd.rst index bc16a3de6f..8f2726e8a3 100644 --- a/doc/src/compute_msd.rst +++ b/doc/src/compute_msd.rst @@ -116,7 +116,9 @@ Compute *msd* cannot be used with a dynamic group. Related commands """""""""""""""" -:doc:`compute msd/nongauss `, :doc:`compute displace_atom `, :doc:`fix store/state `, :doc:`compute msd/chunk ` +:doc:`compute msd/nongauss `, +:doc:`compute displace_atom `, :doc:`fix store/state `, +:doc:`compute msd/chunk ` Default """"""" diff --git a/doc/src/compute_msd_chunk.rst b/doc/src/compute_msd_chunk.rst index db6e1e6fc2..863be4db35 100644 --- a/doc/src/compute_msd_chunk.rst +++ b/doc/src/compute_msd_chunk.rst @@ -131,7 +131,7 @@ Restrictions Related commands """""""""""""""" -:doc:`compute msd ` +:doc:`compute msd `, :doc:`compute vacf/chunk ` Default """"""" diff --git a/doc/src/compute_vacf.rst b/doc/src/compute_vacf.rst index 704e597e18..72f675360c 100644 --- a/doc/src/compute_vacf.rst +++ b/doc/src/compute_vacf.rst @@ -76,7 +76,7 @@ Restrictions Related commands """""""""""""""" -:doc:`compute msd ` +:doc:`compute msd `, :doc:`compute vacf/chunk ` Default """"""" diff --git a/doc/src/compute_vacf_chunk.rst b/doc/src/compute_vacf_chunk.rst new file mode 100644 index 0000000000..214abe3845 --- /dev/null +++ b/doc/src/compute_vacf_chunk.rst @@ -0,0 +1,122 @@ +.. index:: compute vacf/chunk + +compute vacf/chunk command +========================== + +Syntax +"""""" + +.. code-block:: LAMMPS + + compute ID group-ID vacf/chunk chunkID + +* ID, group-ID are documented in :doc:`compute ` command +* vacf/chunk = style name of this compute command +* chunkID = ID of :doc:`compute chunk/atom ` command + +Examples +"""""""" + +.. code-block:: LAMMPS + + compute 1 all vacf/chunk molchunk + +Description +""""""""""" + +Define a computation that calculates the velocity auto-correlation +function (VACF) for multiple chunks of atoms. + +In LAMMPS, chunks are collections of atoms defined by a :doc:`compute +chunk/atom ` command, which assigns each atom to a +single chunk (or no chunk). The ID for this command is specified as +chunkID. For example, a single chunk could be the atoms in a molecule +or atoms in a spatial bin. See the :doc:`compute chunk/atom +` and :doc:`Howto chunk ` doc pages for +details of how chunks can be defined and examples of how they can be +used to measure properties of a system. + +Four quantities are calculated by this compute for each chunk. The +first 3 quantities are the product of the initial center of mass +velocity (VCM) for each chunk in *x*, *y*, and *z* direction with the +current center of mass velocity in the same direction. The fourth +component is the total VACF, i.e. the sum of the three components. + +Note that only atoms in the specified group contribute to the +calculation. The :doc:`compute chunk/atom ` command +defines its own group; atoms will have a chunk ID = 0 if they are not in +that group, signifying they are not assigned to a chunk, and will thus +also not contribute to this calculation. You can specify the "all" +group for this command if you simply want to include atoms with non-zero +chunk IDs. + +The integral of the VACF versus time is proportional to the diffusion +coefficient of the diffusing chunks. + +.. note:: + + The number of chunks *Nchunk* calculated by the + :doc:`compute chunk/atom ` command must remain constant + each time this compute is invoked, so that the dot product for each chunk + from its original position can be computed consistently. If *Nchunk* + does not remain constant, an error will be generated. If needed, you + can enforce a constant *Nchunk* by using the *nchunk once* or *ids once* + options when specifying the :doc:`compute chunk/atom ` + command. + +.. note:: + + This compute stores the original center-of-mass velocities of each + chunk. When a VACF is calculated on a later timestep, it is assumed + that the same atoms are assigned to the same chunk ID. However + LAMMPS has no simple way to ensure this is the case, though you can + use the *ids once* option when specifying the :doc:`compute + chunk/atom ` command. Note that if this is not + the case, the VACF calculation does not have a sensible meaning. + +.. note:: + + If you want the quantities calculated by this compute to be + continuous when running from a :doc:`restart file `, then + you should use the same ID for this compute, as in the original run. + This is so that the fix this compute creates to store per-chunk + quantities will also have the same ID, and thus be initialized + correctly with chunk reference positions from the restart file. + +The simplest way to output the results of the compute vacf/chunk +calculation to a file is to use the :doc:`fix ave/time ` +command, for example: + +.. code-block:: LAMMPS + + compute cc1 all chunk/atom molecule + compute myChunk all vacf/chunk cc1 + fix 1 all ave/time 100 1 100 c_myChunk[*] file tmp.out mode vector + +Output info +""""""""""" + +This compute calculates a global array where the number of rows = the +number of chunks *Nchunk* as calculated by the specified :doc:`compute +chunk/atom ` command. The number of columns = 4 for +the *x*, *y*, *z*, component and the total VACF. These values can be +accessed by any command that uses global array values from a compute as +input. See the :doc:`Howto output ` page for an overview +of LAMMPS output options. + +The array values are "intensive". The array values will be in +distance\ :math:`^2` divided by time\ :math:`^2` :doc:`units `. + +Restrictions +"""""""""""" + none + +Related commands +"""""""""""""""" + +:doc:`compute vacf `, :doc:`compute msd/chunk ` + +Default +""""""" + +none diff --git a/src/compute_msd_chunk.cpp b/src/compute_msd_chunk.cpp index d738b835fa..98d08f1fcc 100644 --- a/src/compute_msd_chunk.cpp +++ b/src/compute_msd_chunk.cpp @@ -21,6 +21,7 @@ #include "group.h" #include "memory.h" #include "modify.h" +#include "update.h" using namespace LAMMPS_NS; @@ -30,7 +31,7 @@ ComputeMSDChunk::ComputeMSDChunk(LAMMPS *lmp, int narg, char **arg) : ComputeChunk(lmp, narg, arg), id_fix(nullptr), fix(nullptr), massproc(nullptr), masstotal(nullptr), com(nullptr), comall(nullptr), msd(nullptr) { - if (narg != 4) error->all(FLERR, "Illegal compute msd/chunk command"); + if (narg != 4) error->all(FLERR, "Incorrect number of arguments for compute msd/chunk"); msdnchunk = 0; array_flag = 1; @@ -114,6 +115,8 @@ void ComputeMSDChunk::setup() void ComputeMSDChunk::compute_array() { + invoked_array = update->ntimestep; + int index; double massone; double unwrap[3]; @@ -127,7 +130,7 @@ void ComputeMSDChunk::compute_array() if (firstflag) msdnchunk = nchunk; else if (msdnchunk != nchunk) - error->all(FLERR, "Compute msd/chunk nchunk is not static"); + error->all(FLERR, Error::NOLASTLINE, "Compute msd/chunk nchunk is not static"); // zero local per-chunk values diff --git a/src/compute_vacf_chunk.cpp b/src/compute_vacf_chunk.cpp new file mode 100644 index 0000000000..7e058b4e85 --- /dev/null +++ b/src/compute_vacf_chunk.cpp @@ -0,0 +1,225 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + 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_vacf_chunk.h" + +#include "atom.h" +#include "compute_chunk_atom.h" +#include "error.h" +#include "fix_store_global.h" +#include "group.h" +#include "memory.h" +#include "modify.h" +#include "update.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +ComputeVACFChunk::ComputeVACFChunk(LAMMPS *lmp, int narg, char **arg) : + ComputeChunk(lmp, narg, arg), id_fix(nullptr), fix(nullptr), massproc(nullptr), + masstotal(nullptr), vcm(nullptr), vcmall(nullptr), vacf(nullptr) +{ + if (narg != 4) error->all(FLERR, "Incorrect number of arguments for compute vacf/chunk"); + + vacfnchunk = 0; + array_flag = 1; + size_array_cols = 4; + size_array_rows = 0; + size_array_rows_variable = 1; + extarray = 0; + + ComputeVACFChunk::init(); + + // create a new fix STORE style for reference velocities + // id = compute-ID + COMPUTE_STORE, fix group = compute group + // do not know size of array at this point, just allocate 1x1 array + // fix creation must be done now so that a restart run can + // potentially re-populate the fix array (and change it to correct size) + // otherwise size reset and init will be done in setup() + + id_fix = utils::strdup(std::string(id) + "_COMPUTE_STORE"); + fix = dynamic_cast( + modify->add_fix(fmt::format("{} {} STORE/GLOBAL 1 1", id_fix, group->names[igroup]))); +} + +/* ---------------------------------------------------------------------- */ + +ComputeVACFChunk::~ComputeVACFChunk() +{ + // check nfix in case all fixes have already been deleted + + if (modify->nfix) modify->delete_fix(id_fix); + + delete[] id_fix; + memory->destroy(massproc); + memory->destroy(masstotal); + memory->destroy(vcm); + memory->destroy(vcmall); + memory->destroy(vacf); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeVACFChunk::init() +{ + ComputeChunk::init(); + + // set fix which stores reference atom coords + // if firstflag, will be created in setup() + + if (!firstflag) { + fix = dynamic_cast(modify->get_fix_by_id(id_fix)); + if (!fix) error->all(FLERR, "Could not find compute vacf/chunk fix with ID {}", id_fix); + } +} + +/* ---------------------------------------------------------------------- + compute initial VCM for each chunk + only once on timestep compute is defined, when firstflag = 1 +------------------------------------------------------------------------- */ + +void ComputeVACFChunk::setup() +{ + if (!firstflag) return; + compute_array(); + firstflag = 0; + + // if fix->astore is already correct size, restart file set it up + // otherwise reset its size now and initialize to current VCM + + if (fix->nrow == nchunk && fix->ncol == 3) return; + fix->reset_global(nchunk, 3); + + double **vcminit = fix->astore; + for (int i = 0; i < nchunk; i++) { + vcminit[i][0] = vcmall[i][0]; + vcminit[i][1] = vcmall[i][1]; + vcminit[i][2] = vcmall[i][2]; + vacf[i][0] = vacf[i][1] = vacf[i][2] = vacf[i][3] = 1.0; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeVACFChunk::compute_array() +{ + invoked_array = update->ntimestep; + + int index; + double massone; + + ComputeChunk::compute_array(); + int *ichunk = cchunk->ichunk; + + // first time call, allocate per-chunk arrays + // thereafter, require nchunk remain the same + + if (firstflag) + vacfnchunk = nchunk; + else if (vacfnchunk != nchunk) + error->all(FLERR, Error::NOLASTLINE, "Compute vacf/chunk nchunk is not static"); + + // zero local per-chunk values + + for (int i = 0; i < nchunk; i++) { + massproc[i] = 0.0; + vcm[i][0] = vcm[i][1] = vcm[i][2] = 0.0; + } + + // compute current VCM for each chunk + + double **v = atom->v; + int *mask = atom->mask; + int *type = atom->type; + double *mass = atom->mass; + double *rmass = atom->rmass; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + index = ichunk[i] - 1; + if (index < 0) continue; + if (rmass) + massone = rmass[i]; + else + massone = mass[type[i]]; + massproc[index] += massone; + vcm[index][0] += v[i][0] * massone; + vcm[index][1] += v[i][1] * massone; + vcm[index][2] += v[i][2] * massone; + } + + MPI_Allreduce(massproc, masstotal, nchunk, MPI_DOUBLE, MPI_SUM, world); + MPI_Allreduce(&vcm[0][0], &vcmall[0][0], 3 * nchunk, MPI_DOUBLE, MPI_SUM, world); + + for (int i = 0; i < nchunk; i++) { + if (masstotal[i] > 0.0) { + vcmall[i][0] /= masstotal[i]; + vcmall[i][1] /= masstotal[i]; + vcmall[i][2] /= masstotal[i]; + } + } + + // VACF is dot product between current and initial VCM + // vcminit is initilialized by setup() when firstflag is set + + if (firstflag) return; + + double vxsq, vysq, vzsq; + double **vcminit = fix->astore; + + for (int i = 0; i < nchunk; i++) { + vxsq = vcmall[i][0] * vcminit[i][0]; + vysq = vcmall[i][1] * vcminit[i][1]; + vzsq = vcmall[i][2] * vcminit[i][2]; + vacf[i][0] = vxsq; + vacf[i][1] = vysq; + vacf[i][2] = vzsq; + vacf[i][3] = vxsq + vysq + vzsq; + } +} + +/* ---------------------------------------------------------------------- + one-time allocate of per-chunk arrays +------------------------------------------------------------------------- */ + +void ComputeVACFChunk::allocate() +{ + ComputeChunk::allocate(); + memory->destroy(massproc); + memory->destroy(masstotal); + memory->destroy(vcm); + memory->destroy(vcmall); + memory->destroy(vacf); + + memory->create(massproc, nchunk, "vacf/chunk:massproc"); + memory->create(masstotal, nchunk, "vacf/chunk:masstotal"); + memory->create(vcm, nchunk, 3, "vacf/chunk:vcm"); + memory->create(vcmall, nchunk, 3, "vacf/chunk:vcmall"); + memory->create(vacf, nchunk, 4, "vacf/chunk:vacf"); + array = vacf; +} + +/* ---------------------------------------------------------------------- + memory usage of local data +------------------------------------------------------------------------- */ + +double ComputeVACFChunk::memory_usage() +{ + double bytes = ComputeChunk::memory_usage(); + bytes += (bigint) nchunk * 2 * sizeof(double); + bytes += (double) nchunk * 2 * 3 * sizeof(double); + bytes += (double) nchunk * 4 * sizeof(double); + return bytes; +} diff --git a/src/compute_vacf_chunk.h b/src/compute_vacf_chunk.h new file mode 100644 index 0000000000..1501960422 --- /dev/null +++ b/src/compute_vacf_chunk.h @@ -0,0 +1,52 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + 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(vacf/chunk,ComputeVACFChunk); +// clang-format on +#else + +#ifndef LMP_COMPUTE_VACF_CHUNK_H +#define LMP_COMPUTE_VACF_CHUNK_H + +#include "compute_chunk.h" + +namespace LAMMPS_NS { + +class ComputeVACFChunk : public ComputeChunk { + public: + ComputeVACFChunk(class LAMMPS *, int, char **); + ~ComputeVACFChunk() override; + void init() override; + void setup() override; + + void compute_array() override; + + double memory_usage() override; + + protected: + bigint vacfnchunk; + + char *id_fix; + class FixStoreGlobal *fix; + + double *massproc, *masstotal; + double **vcm, **vcmall; + double **vacf; + + void allocate() override; +}; +} // namespace LAMMPS_NS +#endif +#endif