From c9e300f76f4ab31bb1a5fa1d8aa03bfcb3c71a56 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 3 Mar 2023 16:18:31 -0500 Subject: [PATCH 01/14] implementation of bond style harmonic/restrain /w docs --- doc/src/Commands_bond.rst | 1 + doc/src/bond_harmonic_restrain.rst | 79 ++++++ doc/src/bond_style.rst | 3 +- doc/utils/sphinx-config/false_positives.txt | 1 + src/.gitignore | 2 + src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp | 244 ++++++++++++++++++ src/EXTRA-MOLECULE/bond_harmonic_restrain.h | 52 ++++ 7 files changed, 381 insertions(+), 1 deletion(-) create mode 100644 doc/src/bond_harmonic_restrain.rst create mode 100644 src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp create mode 100644 src/EXTRA-MOLECULE/bond_harmonic_restrain.h diff --git a/doc/src/Commands_bond.rst b/doc/src/Commands_bond.rst index cfc896aa0e..aaf706b5df 100644 --- a/doc/src/Commands_bond.rst +++ b/doc/src/Commands_bond.rst @@ -42,6 +42,7 @@ OPT. * :doc:`gaussian ` * :doc:`gromos (o) ` * :doc:`harmonic (iko) ` + * :doc:`harmonic/restrain ` * :doc:`harmonic/shift (o) ` * :doc:`harmonic/shift/cut (o) ` * :doc:`lepton (o) ` diff --git a/doc/src/bond_harmonic_restrain.rst b/doc/src/bond_harmonic_restrain.rst new file mode 100644 index 0000000000..a9918d8b5d --- /dev/null +++ b/doc/src/bond_harmonic_restrain.rst @@ -0,0 +1,79 @@ +.. index:: bond_style harmonic/restrain + +bond_style harmonic/restrain command +==================================== + +Syntax +"""""" + +.. code-block:: LAMMPS + + bond_style harmonic/restrain + +Examples +"""""""" + +.. code-block:: LAMMPS + + bond_style harmonic + bond_coeff 5 80.0 + +Description +""""""""""" + +The *harmonic/restrain* bond style uses the potential + +.. math:: + + E = K (r - r_{t=0})^2 + +where :math:`r_{t=0}` is the distance between the bond atoms at the +beginning of the first :doc:`run ` or :doc:`minimize ` +command after the bond style has been defined (*t=0*). Note that the +usual 1/2 factor is included in :math:`K`. This will effectively +restrain bonds to their initial length, whatever that is. + +The following coefficient must be defined for each bond type via the +:doc:`bond_coeff ` command as in the example above, or in +the data file or restart files read by the :doc:`read_data ` +or :doc:`read_restart ` commands + +* :math:`K` (energy/distance\^2) + +Restart info +"""""""""""" + +This bond style supports the :doc:`write_restart ` and +:doc:`read_restart ` commands. The state of the initial +bond lengths is stored with the restart files and read back. + +Restrictions +"""""""""""" + +This bond style can only be used if LAMMPS was built with the +EXTRA-MOLECULE package. See the :doc:`Build package ` +page for more info. + +This bond style does not write its coefficients to a data file since the +:math:`r_{t=0}` values are for individual bonds and not bond types. +Since it uses :doc:`fix property/atom ` internally, +you should use the "nofix" argument to the :doc:`write_data command +` to avoid writing a section to the data file can cannot be +read back in by the bond style. + +This bond style cannot be used with :doc:`fix shake or fix rattle +`, with :doc:`fix filter/corotate `, or +any :doc:`tip4p pair style ` since there is no specific +equilibrium distance for a given bond type. + + +Related commands +"""""""""""""""" + +:doc:`bond_coeff `, :doc:`delete_bonds `, +:doc:`bond_harmonic ` + +Default +""""""" + +none diff --git a/doc/src/bond_style.rst b/doc/src/bond_style.rst index 23b89d00a2..b33d0a9e9a 100644 --- a/doc/src/bond_style.rst +++ b/doc/src/bond_style.rst @@ -10,7 +10,7 @@ Syntax bond_style style args -* style = *none* or *zero* or *hybrid* or *bpm/rotational* or *bpm/spring* or *class2* or *fene* or *fene/expand* or *fene/nm* or *gaussian* or *gromos* or *harmonic* or *harmonic/shift* or *harmonic/shift/cut* or *lepton* or *morse* or *nonlinear* or *oxdna/fene* or *oxdena2/fene* or *oxrna2/fene* or *quartic* or *special* or *table* +* style = *none* or *zero* or *hybrid* or *bpm/rotational* or *bpm/spring* or *class2* or *fene* or *fene/expand* or *fene/nm* or *gaussian* or *gromos* or *harmonic* or *harmonic/restrain* *harmonic/shift* or *harmonic/shift/cut* or *lepton* or *morse* or *nonlinear* or *oxdna/fene* or *oxdena2/fene* or *oxrna2/fene* or *quartic* or *special* or *table* * args = none for any style except *hybrid* @@ -93,6 +93,7 @@ accelerated styles exist. * :doc:`gaussian ` - multicentered Gaussian-based bond potential * :doc:`gromos ` - GROMOS force field bond * :doc:`harmonic ` - harmonic bond +* :doc:`harmonic/restrain ` - harmonic bond to restrain to original bond distance * :doc:`harmonic/shift ` - shifted harmonic bond * :doc:`harmonic/shift/cut ` - shifted harmonic bond with a cutoff * :doc:`lepton ` - bond potential from evaluating a string diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index d03300439f..ff8f0d455b 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -2463,6 +2463,7 @@ nodeless nodeset nodesets Noehring +nofix Noffset noforce noguess diff --git a/src/.gitignore b/src/.gitignore index df02ab3397..aa57927669 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -492,6 +492,8 @@ /bond_gromos.h /bond_harmonic.cpp /bond_harmonic.h +/bond_harmonic_restrain.cpp +/bond_harmonic_restrain.h /bond_harmonic_shift.cpp /bond_harmonic_shift.h /bond_harmonic_shift_cut.cpp diff --git a/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp b/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp new file mode 100644 index 0000000000..413d00226d --- /dev/null +++ b/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp @@ -0,0 +1,244 @@ +/* ---------------------------------------------------------------------- + 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 "bond_harmonic_restrain.h" + +#include "atom.h" +#include "comm.h" +#include "error.h" +#include "fix_property_atom.h" +#include "force.h" +#include "memory.h" +#include "modify.h" +#include "neighbor.h" + +#include +#include + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +BondHarmonicRestrain::BondHarmonicRestrain(LAMMPS *_lmp) : Bond(_lmp), initial(nullptr) +{ + writedata = 0; + natoms = -1; +} + +/* ---------------------------------------------------------------------- */ + +BondHarmonicRestrain::~BondHarmonicRestrain() +{ + modify->delete_fix("BOND_RESTRAIN"); + if (allocated) { + memory->destroy(setflag); + memory->destroy(k); + } +} + +/* ---------------------------------------------------------------------- */ + +void BondHarmonicRestrain::compute(int eflag, int vflag) +{ + int i1, i2, n, type; + double delx, dely, delz, ebond, fbond; + double rsq, r, r0, dr, rk; + + ebond = 0.0; + ev_init(eflag, vflag); + + double **x = atom->x; + double **x0 = (double **) atom->extract("d2_x0"); + double **f = atom->f; + int **bondlist = neighbor->bondlist; + int nbondlist = neighbor->nbondlist; + int nlocal = atom->nlocal; + int newton_bond = force->newton_bond; + + for (n = 0; n < nbondlist; n++) { + i1 = bondlist[n][0]; + i2 = bondlist[n][1]; + type = bondlist[n][2]; + + delx = x0[i1][0] - x0[i2][0]; + dely = x0[i1][1] - x0[i2][1]; + delz = x0[i1][2] - x0[i2][2]; + rsq = delx * delx + dely * dely + delz * delz; + r0 = sqrt(rsq); + + delx = x[i1][0] - x[i2][0]; + dely = x[i1][1] - x[i2][1]; + delz = x[i1][2] - x[i2][2]; + + rsq = delx * delx + dely * dely + delz * delz; + r = sqrt(rsq); + dr = r - r0; + rk = k[type] * dr; + + // force & energy + + if (r > 0.0) + fbond = -2.0 * rk / r; + else + fbond = 0.0; + + if (eflag) ebond = rk * dr; + + // apply force to each of 2 atoms + + if (newton_bond || i1 < nlocal) { + f[i1][0] += delx * fbond; + f[i1][1] += dely * fbond; + f[i1][2] += delz * fbond; + } + + if (newton_bond || i2 < nlocal) { + f[i2][0] -= delx * fbond; + f[i2][1] -= dely * fbond; + f[i2][2] -= delz * fbond; + } + + if (evflag) ev_tally(i1, i2, nlocal, newton_bond, ebond, fbond, delx, dely, delz); + } +} + +/* ---------------------------------------------------------------------- */ + +void BondHarmonicRestrain::allocate() +{ + allocated = 1; + const int np1 = atom->nbondtypes + 1; + + memory->create(k, np1, "bond:k"); + + memory->create(setflag, np1, "bond:setflag"); + for (int i = 1; i < np1; i++) setflag[i] = 0; +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more types +------------------------------------------------------------------------- */ + +void BondHarmonicRestrain::coeff(int narg, char **arg) +{ + if (narg != 2) error->all(FLERR, "Incorrect args for bond coefficients"); + if (!allocated) allocate(); + + int ilo, ihi; + utils::bounds(FLERR, arg[0], 1, atom->nbondtypes, ilo, ihi, error); + + double k_one = utils::numeric(FLERR, arg[1], false, lmp); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + k[i] = k_one; + setflag[i] = 1; + count++; + } + + if (count == 0) error->all(FLERR, "Incorrect args for bond coefficients"); +} + +/* ---------------------------------------------------------------------- + initialize custom data storage +------------------------------------------------------------------------- */ + +void BondHarmonicRestrain::init_style() +{ + // create internal fix to store initial positions + initial = dynamic_cast( + modify->add_fix("BOND_RESTRAIN all property/atom d2_x0 3 ghost yes")); + if (!initial) error->all(FLERR, "Failure to create internal per-atom storage"); + + // store initial positions + if (natoms < 0) { + natoms = atom->natoms; + double *const *const x0 = (double **) atom->extract("d2_x0"); + const double *const *const x = atom->x; + for (int i = 0; i < atom->nlocal; ++i) + for (int j = 0; j < 3; ++j) x0[i][j] = x[i][j]; + } + + // must not add atoms + if (natoms < atom->natoms) + error->all(FLERR, "Bond style harmonic/restrain does not support adding atoms"); +} + +/* ---------------------------------------------------------------------- + proc 0 writes out coeffs to restart file +------------------------------------------------------------------------- */ + +void BondHarmonicRestrain::write_restart(FILE *fp) +{ + fwrite(&natoms, sizeof(bigint), 1, fp); + fwrite(&k[1], sizeof(double), atom->nbondtypes, fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads coeffs from restart file, bcasts them +------------------------------------------------------------------------- */ + +void BondHarmonicRestrain::read_restart(FILE *fp) +{ + allocate(); + + if (comm->me == 0) { + utils::sfread(FLERR, &natoms, sizeof(bigint), 1, fp, nullptr, error); + utils::sfread(FLERR, &k[1], sizeof(double), atom->nbondtypes, fp, nullptr, error); + } + MPI_Bcast(&natoms, 1, MPI_LMP_BIGINT, 0, world); + MPI_Bcast(&k[1], atom->nbondtypes, MPI_DOUBLE, 0, world); + + for (int i = 1; i <= atom->nbondtypes; i++) setflag[i] = 1; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to data file +------------------------------------------------------------------------- */ + +void BondHarmonicRestrain::write_data(FILE *fp) +{ + for (int i = 1; i <= atom->nbondtypes; i++) fprintf(fp, "%d %g\n", i, k[i]); +} + +/* ---------------------------------------------------------------------- */ + +double BondHarmonicRestrain::single(int type, double rsq, int i, int j, double &fforce) +{ + double **x0 = (double **) atom->extract("d2_x0"); + + double delx = x0[i][0] - x0[j][0]; + double dely = x0[i][1] - x0[j][1]; + double delz = x0[i][2] - x0[j][2]; + double r0 = sqrt(delx * delx + dely * dely + delz * delz); + + double r = sqrt(rsq); + double dr = r - r0; + double rk = k[type] * dr; + fforce = 0; + if (r > 0.0) fforce = -2.0 * rk / r; + return rk * dr; +} + +/* ---------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + return ptr to internal members upon request +------------------------------------------------------------------------ */ + +void *BondHarmonicRestrain::extract(const char *str, int &dim) +{ + dim = 1; + if (strcmp(str, "k") == 0) return (void *) k; + return nullptr; +} diff --git a/src/EXTRA-MOLECULE/bond_harmonic_restrain.h b/src/EXTRA-MOLECULE/bond_harmonic_restrain.h new file mode 100644 index 0000000000..e17c21e125 --- /dev/null +++ b/src/EXTRA-MOLECULE/bond_harmonic_restrain.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 BOND_CLASS +// clang-format off +BondStyle(harmonic/restrain,BondHarmonicRestrain); +// clang-format on +#else + +#ifndef LMP_BOND_HARMONIC_RESTRAIN_H +#define LMP_BOND_HARMONIC_RESTRAIN_H + +#include "bond.h" + +namespace LAMMPS_NS { + +class BondHarmonicRestrain : public Bond { + public: + BondHarmonicRestrain(class LAMMPS *); + ~BondHarmonicRestrain() override; + void compute(int, int) override; + void coeff(int, char **) override; + void init_style() override; + double equilibrium_distance(int) override { return -1.0; }; // return invalid value since not applicable + void write_restart(FILE *) override; + void read_restart(FILE *) override; + void write_data(FILE *) override; + double single(int, double, int, int, double &) override; + void *extract(const char *, int &) override; + + protected: + double *k; + bigint natoms; + class FixPropertyAtom *initial; + + virtual void allocate(); +}; + +} // namespace LAMMPS_NS + +#endif +#endif From 77ae215047b6b5ec65239f42aac8e9028347cfb2 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 3 Mar 2023 17:16:12 -0500 Subject: [PATCH 02/14] add versionadded tag --- doc/src/bond_harmonic_restrain.rst | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/doc/src/bond_harmonic_restrain.rst b/doc/src/bond_harmonic_restrain.rst index a9918d8b5d..5ffbfb424c 100644 --- a/doc/src/bond_harmonic_restrain.rst +++ b/doc/src/bond_harmonic_restrain.rst @@ -21,6 +21,8 @@ Examples Description """"""""""" +.. versionadded:: TBD + The *harmonic/restrain* bond style uses the potential .. math:: @@ -46,7 +48,7 @@ Restart info This bond style supports the :doc:`write_restart ` and :doc:`read_restart ` commands. The state of the initial bond lengths is stored with the restart files and read back. - + Restrictions """""""""""" From d32411f61a37f671afdc610316982a51f66e9c3d Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 7 Mar 2023 16:37:14 -0500 Subject: [PATCH 03/14] recover restart ability and avoid memory leak --- src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp | 27 ++++++++++++++----- 1 file changed, 21 insertions(+), 6 deletions(-) diff --git a/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp b/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp index 413d00226d..76e7347ada 100644 --- a/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp +++ b/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp @@ -15,6 +15,7 @@ #include "atom.h" #include "comm.h" +#include "domain.h" #include "error.h" #include "fix_property_atom.h" #include "force.h" @@ -73,6 +74,7 @@ void BondHarmonicRestrain::compute(int eflag, int vflag) delx = x0[i1][0] - x0[i2][0]; dely = x0[i1][1] - x0[i2][1]; delz = x0[i1][2] - x0[i2][2]; + domain->minimum_image(delx, dely, delz); rsq = delx * delx + dely * dely + delz * delz; r0 = sqrt(rsq); @@ -155,19 +157,32 @@ void BondHarmonicRestrain::coeff(int narg, char **arg) void BondHarmonicRestrain::init_style() { - // create internal fix to store initial positions - initial = dynamic_cast( - modify->add_fix("BOND_RESTRAIN all property/atom d2_x0 3 ghost yes")); - if (!initial) error->all(FLERR, "Failure to create internal per-atom storage"); - // store initial positions + if (natoms < 0) { + + // create internal fix to store initial positions + initial = dynamic_cast( + modify->add_fix("BOND_RESTRAIN all property/atom d2_x0 3 ghost yes")); + if (!initial) error->all(FLERR, "Failure to create internal per-atom storage"); + natoms = atom->natoms; double *const *const x0 = (double **) atom->extract("d2_x0"); const double *const *const x = atom->x; for (int i = 0; i < atom->nlocal; ++i) for (int j = 0; j < 3; ++j) x0[i][j] = x[i][j]; - } + + } else { + + // after a restart, natoms is set but initial is a null pointer. + // we add the fix, but do not initialize it. It will pull the data from the restart. + + if (!initial) { + initial = dynamic_cast( + modify->add_fix("BOND_RESTRAIN all property/atom d2_x0 3 ghost yes")); + if (!initial) error->all(FLERR, "Failure to create internal per-atom storage"); + } + } // must not add atoms if (natoms < atom->natoms) From 010b030b567d9c057b92f227a3ad6f38190fe777 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 7 Mar 2023 19:38:34 -0500 Subject: [PATCH 04/14] update docs and mention alternatives and describe differences --- doc/src/bond_harmonic_restrain.rst | 27 ++++++++++++++++++--------- 1 file changed, 18 insertions(+), 9 deletions(-) diff --git a/doc/src/bond_harmonic_restrain.rst b/doc/src/bond_harmonic_restrain.rst index 5ffbfb424c..d4c34ae622 100644 --- a/doc/src/bond_harmonic_restrain.rst +++ b/doc/src/bond_harmonic_restrain.rst @@ -29,11 +29,13 @@ The *harmonic/restrain* bond style uses the potential E = K (r - r_{t=0})^2 -where :math:`r_{t=0}` is the distance between the bond atoms at the +where :math:`r_{t=0}` is the distance between the bonded atoms at the beginning of the first :doc:`run ` or :doc:`minimize ` command after the bond style has been defined (*t=0*). Note that the usual 1/2 factor is included in :math:`K`. This will effectively -restrain bonds to their initial length, whatever that is. +restrain bonds to their initial length, whatever that is. This is where +this bond style differs from :doc:`bond style harmonic ` +where the bond length is set through the per bond type coefficients. The following coefficient must be defined for each bond type via the :doc:`bond_coeff ` command as in the example above, or in @@ -42,12 +44,19 @@ or :doc:`read_restart ` commands * :math:`K` (energy/distance\^2) +This bond style differs from other options to add harmonic restraints +like :doc:`fix restrain ` or :doc:`pair style list +` or :doc:`fix colvars ` in that it requires a +bond topology, and thus the defined bonds will trigger exclusion of +special neighbors from the neighbor list according to the +:doc:`special_bonds ` settings. + Restart info """""""""""" This bond style supports the :doc:`write_restart ` and :doc:`read_restart ` commands. The state of the initial -bond lengths is stored with the restart files and read back. +bond lengths is stored with restart files and read back. Restrictions """""""""""" @@ -56,24 +65,24 @@ This bond style can only be used if LAMMPS was built with the EXTRA-MOLECULE package. See the :doc:`Build package ` page for more info. -This bond style does not write its coefficients to a data file since the +This bond style does NOT write its coefficients to a data file since the :math:`r_{t=0}` values are for individual bonds and not bond types. + Since it uses :doc:`fix property/atom ` internally, you should use the "nofix" argument to the :doc:`write_data command -` to avoid writing a section to the data file can cannot be -read back in by the bond style. +` to avoid writing a section to the data file that cannot be +read back by the bond style. This bond style cannot be used with :doc:`fix shake or fix rattle `, with :doc:`fix filter/corotate `, or any :doc:`tip4p pair style ` since there is no specific equilibrium distance for a given bond type. - Related commands """""""""""""""" -:doc:`bond_coeff `, :doc:`delete_bonds `, -:doc:`bond_harmonic ` +:doc:`bond_coeff `, :doc:`bond_harmonic `, +:doc:`fix restrain `, :doc:`pair style list ` Default """"""" From fd2cda66bbb22ccfcb2db522a6e4866834074e77 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 7 Mar 2023 20:18:41 -0500 Subject: [PATCH 05/14] handle the case when the bond style is recreated after a restart --- src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp b/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp index 76e7347ada..32e1226347 100644 --- a/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp +++ b/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp @@ -40,7 +40,7 @@ BondHarmonicRestrain::BondHarmonicRestrain(LAMMPS *_lmp) : Bond(_lmp), initial(n BondHarmonicRestrain::~BondHarmonicRestrain() { - modify->delete_fix("BOND_RESTRAIN"); + if (initial) modify->delete_fix("BOND_RESTRAIN"); if (allocated) { memory->destroy(setflag); memory->destroy(k); From 9a8b5ebae0e1d93136fb1b833a250252ddf917af Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 7 Mar 2023 20:36:28 -0500 Subject: [PATCH 06/14] use more unusual name for per-atom property --- src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp b/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp index 32e1226347..2054884e3e 100644 --- a/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp +++ b/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp @@ -40,7 +40,7 @@ BondHarmonicRestrain::BondHarmonicRestrain(LAMMPS *_lmp) : Bond(_lmp), initial(n BondHarmonicRestrain::~BondHarmonicRestrain() { - if (initial) modify->delete_fix("BOND_RESTRAIN"); + if (initial) modify->delete_fix("BOND_RESTRAIN_X0"); if (allocated) { memory->destroy(setflag); memory->destroy(k); @@ -59,7 +59,7 @@ void BondHarmonicRestrain::compute(int eflag, int vflag) ev_init(eflag, vflag); double **x = atom->x; - double **x0 = (double **) atom->extract("d2_x0"); + double **x0 = (double **) atom->extract("d2_BOND_RESTRAIN_X0"); double **f = atom->f; int **bondlist = neighbor->bondlist; int nbondlist = neighbor->nbondlist; @@ -163,11 +163,11 @@ void BondHarmonicRestrain::init_style() // create internal fix to store initial positions initial = dynamic_cast( - modify->add_fix("BOND_RESTRAIN all property/atom d2_x0 3 ghost yes")); + modify->add_fix("BOND_RESTRAIN_X0 all property/atom d2_BOND_RESTRAIN_X0 3 ghost yes")); if (!initial) error->all(FLERR, "Failure to create internal per-atom storage"); natoms = atom->natoms; - double *const *const x0 = (double **) atom->extract("d2_x0"); + double *const *const x0 = (double **) atom->extract("d2_BOND_RESTRAIN_X0"); const double *const *const x = atom->x; for (int i = 0; i < atom->nlocal; ++i) for (int j = 0; j < 3; ++j) x0[i][j] = x[i][j]; @@ -179,7 +179,7 @@ void BondHarmonicRestrain::init_style() if (!initial) { initial = dynamic_cast( - modify->add_fix("BOND_RESTRAIN all property/atom d2_x0 3 ghost yes")); + modify->add_fix("BOND_RESTRAIN_X0 all property/atom d2_BOND_RESTRAIN_X0 3 ghost yes")); if (!initial) error->all(FLERR, "Failure to create internal per-atom storage"); } } @@ -230,11 +230,12 @@ void BondHarmonicRestrain::write_data(FILE *fp) double BondHarmonicRestrain::single(int type, double rsq, int i, int j, double &fforce) { - double **x0 = (double **) atom->extract("d2_x0"); + double **x0 = (double **) atom->extract("d2_BOND_RESTRAIN_X0"); double delx = x0[i][0] - x0[j][0]; double dely = x0[i][1] - x0[j][1]; double delz = x0[i][2] - x0[j][2]; + domain->minimum_image(delx, dely, delz); double r0 = sqrt(delx * delx + dely * dely + delz * delz); double r = sqrt(rsq); From fc8a048662107dc9185331264e59edb1cd07caa2 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 7 Mar 2023 21:22:57 -0500 Subject: [PATCH 07/14] add unit test --- src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp | 4 +- unittest/force-styles/test_bond_style.cpp | 2 +- .../tests/bond-harmonic_restrain.yaml | 90 +++++++++++++++++++ 3 files changed, 93 insertions(+), 3 deletions(-) create mode 100644 unittest/force-styles/tests/bond-harmonic_restrain.yaml diff --git a/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp b/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp index 2054884e3e..a15642f90b 100644 --- a/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp +++ b/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp @@ -163,7 +163,7 @@ void BondHarmonicRestrain::init_style() // create internal fix to store initial positions initial = dynamic_cast( - modify->add_fix("BOND_RESTRAIN_X0 all property/atom d2_BOND_RESTRAIN_X0 3 ghost yes")); + modify->add_fix("BOND_RESTRAIN_X0 all property/atom d2_BOND_RESTRAIN_X0 3 ghost yes")); if (!initial) error->all(FLERR, "Failure to create internal per-atom storage"); natoms = atom->natoms; @@ -182,7 +182,7 @@ void BondHarmonicRestrain::init_style() modify->add_fix("BOND_RESTRAIN_X0 all property/atom d2_BOND_RESTRAIN_X0 3 ghost yes")); if (!initial) error->all(FLERR, "Failure to create internal per-atom storage"); } - } + } // must not add atoms if (natoms < atom->natoms) diff --git a/unittest/force-styles/test_bond_style.cpp b/unittest/force-styles/test_bond_style.cpp index d056fdc876..c723541366 100644 --- a/unittest/force-styles/test_bond_style.cpp +++ b/unittest/force-styles/test_bond_style.cpp @@ -123,7 +123,7 @@ LAMMPS *init_lammps(int argc, char **argv, const TestConfig &cfg, const bool new command("run 0 post no"); command("write_restart " + cfg.basename + ".restart"); - command("write_data " + cfg.basename + ".data"); + command("write_data " + cfg.basename + ".data nofix"); command("write_coeff " + cfg.basename + "-coeffs.in"); return lmp; diff --git a/unittest/force-styles/tests/bond-harmonic_restrain.yaml b/unittest/force-styles/tests/bond-harmonic_restrain.yaml new file mode 100644 index 0000000000..1fd04b818e --- /dev/null +++ b/unittest/force-styles/tests/bond-harmonic_restrain.yaml @@ -0,0 +1,90 @@ +--- +lammps_version: 8 Feb 2023 +tags: generated +date_generated: Tue Mar 7 21:07:27 2023 +epsilon: 2.5e-13 +skip_tests: extract +prerequisites: ! | + atom full + bond harmonic/restrain +pre_commands: ! "" +post_commands: ! "" +input_file: in.fourmol +bond_style: harmonic/restrain +bond_coeff: ! | + 1 250.0 + 2 300.0 + 3 350.0 + 4 650.0 + 5 450.0 +equilibrium: 5 -1 -1 -1 -1 -1 +extract: ! | + k 1 +natoms: 29 +init_energy: 0 +init_stress: ! |2- + 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 +init_forces: ! |2 + 1 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 2 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 3 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 4 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 5 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 6 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 7 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 8 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 9 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 10 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 11 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 12 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 13 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 14 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 15 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 16 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 17 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 18 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 19 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 20 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 21 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 22 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 23 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 24 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 25 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 26 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 27 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 28 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 + 29 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 +run_energy: 0.009384424442608962 +run_stress: ! |2- + 1.3334287184173488e+00 3.1110677999939773e-01 -1.4064657849540130e+00 2.6546279385354110e-01 1.5492956756466219e+00 7.2358876924785698e-01 +run_forces: ! |2 + 1 1.1485941046856593e-01 1.4088049450172524e-01 -8.5852076971833030e-02 + 2 -9.3217481807658653e-02 -7.6867386712704044e-02 1.0932532178570745e-01 + 3 -2.3741413942894266e-01 -8.5667734110971339e-02 -8.2143831089792743e-02 + 4 2.3324188168788909e-01 -6.5094431822773302e-02 1.6591701889629748e-01 + 5 2.6793200848186677e-02 4.1227376782380774e-02 -1.3685685314886964e-01 + 6 2.4862628338772436e-01 -3.6125357500552457e-01 -7.8230634578348013e-01 + 7 -3.8706943762359766e-02 1.9888379369193138e-01 9.8933719851963053e-01 + 8 -3.5267248957722913e-01 4.2911503092501313e-01 1.0190469073820968e-01 + 9 -5.1632048293607458e-02 -5.5343430872859360e-02 -2.2821606693951479e-01 + 10 1.2860877826855494e-01 -7.5971599169706849e-01 -5.5343112593256227e-01 + 11 1.6956002318038377e-01 4.2001247662003161e-01 6.8896503378985918e-01 + 12 3.4231534762621078e-02 -2.4857235824638585e-01 -6.3964642589377518e-01 + 13 -2.7815178906130594e-01 1.0946454990993748e-01 5.1669158882660616e-03 + 14 2.2738751895410908e-01 -4.5437525741145839e-02 5.8956676893113813e-01 + 15 -4.6207378210972273e-03 1.6438094307388113e-01 5.8917445017986604e-02 + 16 -1.1399994473799732e-02 1.1329499720204761e-01 -5.0720152745025260e-01 + 17 -1.1549300733203431e-01 8.0692771502484495e-02 3.0655385964298520e-01 + 18 1.9145703373728828e+00 1.8373130373081787e+00 -3.9519344330792983e-01 + 19 -4.7317441908503255e-01 -3.0353033418925196e-01 -5.7175303978447201e-01 + 20 -1.4413959182878502e+00 -1.5337827031189268e+00 9.6694648309240183e-01 + 21 -3.3244533656237202e-01 -3.0309080808086836e-01 6.5775553694406208e-01 + 22 -2.9549353211149859e-01 -7.2150425050573716e-02 -4.6194669575592789e-01 + 23 6.2793886867387061e-01 3.7524123313144209e-01 -1.9580884118813421e-01 + 24 8.2698432709098157e-01 -1.4448474780063469e+00 1.2133188519739595e+00 + 25 -1.5417378702379472e+00 2.4470962384652978e-01 -1.2964818258023694e+00 + 26 7.1475354314696560e-01 1.2001378541598171e+00 8.3162973828409884e-02 + 27 2.1885007758012770e-01 -5.7423924185981945e-01 2.7745249319135684e-01 + 28 -5.2036258696458004e-01 2.2330647386191435e-01 -3.5136945845348061e-01 + 29 3.0151250938445234e-01 3.5093276799790507e-01 7.3916965262123782e-02 +... From 4220944ace092cb3cf02236745a84ce2e776bf8f Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 8 Mar 2023 00:25:08 -0500 Subject: [PATCH 08/14] cosmetic --- unittest/force-styles/test_angle_style.cpp | 4 ++-- unittest/force-styles/tests/bond-harmonic_restrain.yaml | 1 - 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/unittest/force-styles/test_angle_style.cpp b/unittest/force-styles/test_angle_style.cpp index d9f541df96..e9f4a3f7fc 100644 --- a/unittest/force-styles/test_angle_style.cpp +++ b/unittest/force-styles/test_angle_style.cpp @@ -691,8 +691,8 @@ TEST(AngleStyle, extract) } auto angle = lmp->force->angle; - void *ptr = nullptr; - int dim = 0; + void *ptr = nullptr; + int dim = 0; for (auto extract : test_config.extract) { ptr = angle->extract(extract.first.c_str(), dim); EXPECT_NE(ptr, nullptr); diff --git a/unittest/force-styles/tests/bond-harmonic_restrain.yaml b/unittest/force-styles/tests/bond-harmonic_restrain.yaml index 1fd04b818e..07546775ab 100644 --- a/unittest/force-styles/tests/bond-harmonic_restrain.yaml +++ b/unittest/force-styles/tests/bond-harmonic_restrain.yaml @@ -1,6 +1,5 @@ --- lammps_version: 8 Feb 2023 -tags: generated date_generated: Tue Mar 7 21:07:27 2023 epsilon: 2.5e-13 skip_tests: extract From b6756c03194fbe26781f46bedc448a5d618825ca Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 9 Mar 2023 13:00:59 -0500 Subject: [PATCH 09/14] update docs about using data files --- doc/src/bond_harmonic_restrain.rst | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/doc/src/bond_harmonic_restrain.rst b/doc/src/bond_harmonic_restrain.rst index d4c34ae622..c48e68917c 100644 --- a/doc/src/bond_harmonic_restrain.rst +++ b/doc/src/bond_harmonic_restrain.rst @@ -65,13 +65,18 @@ This bond style can only be used if LAMMPS was built with the EXTRA-MOLECULE package. See the :doc:`Build package ` page for more info. -This bond style does NOT write its coefficients to a data file since the -:math:`r_{t=0}` values are for individual bonds and not bond types. - -Since it uses :doc:`fix property/atom ` internally, -you should use the "nofix" argument to the :doc:`write_data command -` to avoid writing a section to the data file that cannot be -read back by the bond style. +This bond style uses :doc:`fix property/atom ` +internally to be able to determine the original bond lengths. This fix +will write its state to :doc:`data files `, but it is +currently not possible to read that stored data back with +:doc:`read_data ` for such an internally defined fix. You +are thus advised to use the "nofix" argument to the :doc:`write_data +command ` in order to avoid writing that section to the data +file that would cause an error when reading the data file. When +continuing from such a "nofix" data file, the reference bond lengths +will not be the original lengths but those of the current geometry. +Thus restarting is currently *only* possible with :doc:`binary restart +files `. This bond style cannot be used with :doc:`fix shake or fix rattle `, with :doc:`fix filter/corotate `, or From 3815c0ef76b100bdee330a8ed64d2de76e879dae Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Mon, 13 Mar 2023 17:00:54 -0600 Subject: [PATCH 10/14] modify fix STORE/PERATOM and callers to allow ghost comm --- src/AMOEBA/pair_amoeba.cpp | 6 +- src/CORESHELL/compute_temp_cs.cpp | 2 +- src/EXTRA-COMPUTE/compute_hma.cpp | 2 +- src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp | 18 +-- src/EXTRA-MOLECULE/bond_harmonic_restrain.h | 2 +- src/FEP/fix_adapt_fep.cpp | 4 +- src/REPLICA/tad.cpp | 2 +- src/VTK/dump_vtk.cpp | 2 +- src/balance.cpp | 2 +- src/compute_chunk_atom.cpp | 2 +- src/compute_displace_atom.cpp | 2 +- src/compute_msd.cpp | 2 +- src/compute_vacf.cpp | 2 +- src/dump_custom.cpp | 2 +- src/fix_adapt.cpp | 4 +- src/fix_store_peratom.cpp | 145 +++++++++++++----- src/fix_store_peratom.h | 8 +- src/variable.cpp | 2 +- 18 files changed, 140 insertions(+), 69 deletions(-) diff --git a/src/AMOEBA/pair_amoeba.cpp b/src/AMOEBA/pair_amoeba.cpp index 0812fe43f0..319bbd8925 100644 --- a/src/AMOEBA/pair_amoeba.cpp +++ b/src/AMOEBA/pair_amoeba.cpp @@ -861,7 +861,7 @@ void PairAmoeba::init_style() Fix *myfix; if (first_flag) { id_pole = utils::strdup("AMOEBA_pole"); - myfix = modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 13",id_pole,group->names[0])); + myfix = modify->add_fix(fmt::format("{} {} STORE/PERATOM 13 0 0 1",id_pole,group->names[0])); fixpole = dynamic_cast(myfix); } @@ -873,12 +873,12 @@ void PairAmoeba::init_style() if (first_flag && use_pred) { id_udalt = utils::strdup("AMOEBA_udalt"); - myfix = modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 {} 3", + myfix = modify->add_fix(fmt::format("{} {} STORE/PERATOM {} 3 0 1", id_udalt, group->names[0], maxualt)); fixudalt = dynamic_cast(myfix); id_upalt = utils::strdup("AMOEBA_upalt"); - myfix = modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 {} 3", + myfix = modify->add_fix(fmt::format("{} {} STORE/PERATOM {} 3 0 1", id_upalt, group->names[0], maxualt)); fixupalt = dynamic_cast(myfix); } diff --git a/src/CORESHELL/compute_temp_cs.cpp b/src/CORESHELL/compute_temp_cs.cpp index c28d52e5b2..b0ec39e542 100644 --- a/src/CORESHELL/compute_temp_cs.cpp +++ b/src/CORESHELL/compute_temp_cs.cpp @@ -68,7 +68,7 @@ ComputeTempCS::ComputeTempCS(LAMMPS *lmp, int narg, char **arg) : id_fix = utils::strdup(id + std::string("_COMPUTE_STORE")); fix = dynamic_cast( - modify->add_fix(fmt::format("{} {} STORE/PERATOM 0 1", id_fix, group->names[igroup]))); + modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 0 0 0", id_fix, group->names[igroup]))); // set fix store values = 0 for now // fill them in via setup() once Comm::borders() has been called diff --git a/src/EXTRA-COMPUTE/compute_hma.cpp b/src/EXTRA-COMPUTE/compute_hma.cpp index 0f7279a5f5..806228faef 100644 --- a/src/EXTRA-COMPUTE/compute_hma.cpp +++ b/src/EXTRA-COMPUTE/compute_hma.cpp @@ -91,7 +91,7 @@ ComputeHMA::ComputeHMA(LAMMPS *lmp, int narg, char **arg) : id_fix = utils::strdup(std::string(id)+"_COMPUTE_STORE"); fix = dynamic_cast( - modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 3", id_fix, group->names[igroup]))); + modify->add_fix(fmt::format("{} {} STORE/PERATOM 3 0 0 1", id_fix, group->names[igroup]))); // calculate xu,yu,zu for fix store array // skip if reset from restart file diff --git a/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp b/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp index a15642f90b..59bab65713 100644 --- a/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp +++ b/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp @@ -17,7 +17,7 @@ #include "comm.h" #include "domain.h" #include "error.h" -#include "fix_property_atom.h" +#include "fix_store_peratom.h" #include "force.h" #include "memory.h" #include "modify.h" @@ -59,7 +59,7 @@ void BondHarmonicRestrain::compute(int eflag, int vflag) ev_init(eflag, vflag); double **x = atom->x; - double **x0 = (double **) atom->extract("d2_BOND_RESTRAIN_X0"); + double **x0 = initial->astore; double **f = atom->f; int **bondlist = neighbor->bondlist; int nbondlist = neighbor->nbondlist; @@ -162,12 +162,12 @@ void BondHarmonicRestrain::init_style() if (natoms < 0) { // create internal fix to store initial positions - initial = dynamic_cast( - modify->add_fix("BOND_RESTRAIN_X0 all property/atom d2_BOND_RESTRAIN_X0 3 ghost yes")); + initial = dynamic_cast( + modify->add_fix("BOND_RESTRAIN_X0 all STORE/PERATOM 3 0 1 1")); if (!initial) error->all(FLERR, "Failure to create internal per-atom storage"); natoms = atom->natoms; - double *const *const x0 = (double **) atom->extract("d2_BOND_RESTRAIN_X0"); + double **x0 = initial->astore; const double *const *const x = atom->x; for (int i = 0; i < atom->nlocal; ++i) for (int j = 0; j < 3; ++j) x0[i][j] = x[i][j]; @@ -176,10 +176,10 @@ void BondHarmonicRestrain::init_style() // after a restart, natoms is set but initial is a null pointer. // we add the fix, but do not initialize it. It will pull the data from the restart. - + if (!initial) { - initial = dynamic_cast( - modify->add_fix("BOND_RESTRAIN_X0 all property/atom d2_BOND_RESTRAIN_X0 3 ghost yes")); + initial = dynamic_cast( + modify->add_fix("BOND_RESTRAIN_X0 all STORE/PERATOM 3 0 1 1")); if (!initial) error->all(FLERR, "Failure to create internal per-atom storage"); } } @@ -230,7 +230,7 @@ void BondHarmonicRestrain::write_data(FILE *fp) double BondHarmonicRestrain::single(int type, double rsq, int i, int j, double &fforce) { - double **x0 = (double **) atom->extract("d2_BOND_RESTRAIN_X0"); + double **x0 = initial->astore; double delx = x0[i][0] - x0[j][0]; double dely = x0[i][1] - x0[j][1]; diff --git a/src/EXTRA-MOLECULE/bond_harmonic_restrain.h b/src/EXTRA-MOLECULE/bond_harmonic_restrain.h index e17c21e125..61d01ab9b2 100644 --- a/src/EXTRA-MOLECULE/bond_harmonic_restrain.h +++ b/src/EXTRA-MOLECULE/bond_harmonic_restrain.h @@ -41,7 +41,7 @@ class BondHarmonicRestrain : public Bond { protected: double *k; bigint natoms; - class FixPropertyAtom *initial; + class FixStorePeratom *initial; virtual void allocate(); }; diff --git a/src/FEP/fix_adapt_fep.cpp b/src/FEP/fix_adapt_fep.cpp index 2d2409d27f..3f3d26c4d9 100644 --- a/src/FEP/fix_adapt_fep.cpp +++ b/src/FEP/fix_adapt_fep.cpp @@ -212,7 +212,7 @@ void FixAdaptFEP::post_constructor() if (diamflag) { id_fix_diam = utils::strdup(id + std::string("_FIX_STORE_DIAM")); fix_diam = dynamic_cast( - modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 1", id_fix_diam,group->names[igroup]))); + modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 0 0 1", id_fix_diam,group->names[igroup]))); if (fix_diam->restart_reset) fix_diam->restart_reset = 0; else { double *vec = fix_diam->vstore; @@ -230,7 +230,7 @@ void FixAdaptFEP::post_constructor() if (chgflag) { id_fix_chg = utils::strdup(id + std::string("_FIX_STORE_CHG")); fix_chg = dynamic_cast( - modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 1",id_fix_chg,group->names[igroup]))); + modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 0 0 1",id_fix_chg,group->names[igroup]))); if (fix_chg->restart_reset) fix_chg->restart_reset = 0; else { double *vec = fix_chg->vstore; diff --git a/src/REPLICA/tad.cpp b/src/REPLICA/tad.cpp index d89cbe80ba..e98cd39233 100644 --- a/src/REPLICA/tad.cpp +++ b/src/REPLICA/tad.cpp @@ -135,7 +135,7 @@ void TAD::command(int narg, char **arg) // create FixStorePeratom object to store revert state - fix_revert = dynamic_cast(modify->add_fix("tad_revert all STORE/PERATOM 0 7")); + fix_revert = dynamic_cast(modify->add_fix("tad_revert all STORE/PERATOM 7 0 0 0")); // create Finish for timing output diff --git a/src/VTK/dump_vtk.cpp b/src/VTK/dump_vtk.cpp index a188a1f060..c4dd68bf3c 100644 --- a/src/VTK/dump_vtk.cpp +++ b/src/VTK/dump_vtk.cpp @@ -2365,7 +2365,7 @@ int DumpVTK::modify_param(int narg, char **arg) std::string threshid = fmt::format("{}{}_DUMP_STORE",id,nthreshlast); thresh_fixID[nthreshlast] = utils::strdup(threshid); - threshid += fmt::format(" {} STORE/PERATOM 1 1", group->names[igroup]); + threshid += fmt::format(" {} STORE/PERATOM 1 0 0 1", group->names[igroup]); thresh_fix[nthreshlast] = dynamic_cast(modify->add_fix(threshid)); thresh_last[nthreshlast] = nthreshlast; diff --git a/src/balance.cpp b/src/balance.cpp index 91b9d70378..8d1ba460db 100644 --- a/src/balance.cpp +++ b/src/balance.cpp @@ -515,7 +515,7 @@ void Balance::weight_storage(char *prefix) fixstore = dynamic_cast(modify->get_fix_by_id(cmd)); if (!fixstore) - fixstore = dynamic_cast(modify->add_fix(cmd + " all STORE/PERATOM 0 1")); + fixstore = dynamic_cast(modify->add_fix(cmd + " all STORE/PERATOM 1 0 0 0")); // do not carry weights with atoms during normal atom migration diff --git a/src/compute_chunk_atom.cpp b/src/compute_chunk_atom.cpp index f9723497b1..b59449989b 100644 --- a/src/compute_chunk_atom.cpp +++ b/src/compute_chunk_atom.cpp @@ -572,7 +572,7 @@ void ComputeChunkAtom::init() if ((idsflag == ONCE || lockcount) && !fixstore) { id_fix = utils::strdup(id + std::string("_COMPUTE_STORE")); fixstore = dynamic_cast( - modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 1", id_fix, group->names[igroup]))); + modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 0 0 1", id_fix, group->names[igroup]))); } if ((idsflag != ONCE && !lockcount) && fixstore) { diff --git a/src/compute_displace_atom.cpp b/src/compute_displace_atom.cpp index 0b1ffff687..6f352fa111 100644 --- a/src/compute_displace_atom.cpp +++ b/src/compute_displace_atom.cpp @@ -75,7 +75,7 @@ ComputeDisplaceAtom::ComputeDisplaceAtom(LAMMPS *lmp, int narg, char **arg) : id_fix = utils::strdup(std::string(id) + "_COMPUTE_STORE"); fix = dynamic_cast( - modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 3", id_fix, group->names[igroup]))); + modify->add_fix(fmt::format("{} {} STORE/PERATOM 3 0 0 1", id_fix, group->names[igroup]))); // calculate xu,yu,zu for fix store array // skip if reset from restart file diff --git a/src/compute_msd.cpp b/src/compute_msd.cpp index 8b7043e20d..ecbb540a48 100644 --- a/src/compute_msd.cpp +++ b/src/compute_msd.cpp @@ -64,7 +64,7 @@ ComputeMSD::ComputeMSD(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, a id_fix = utils::strdup(id + std::string("_COMPUTE_STORE")); fix = dynamic_cast( - modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 3", id_fix, group->names[igroup]))); + modify->add_fix(fmt::format("{} {} STORE/PERATOM 3 0 0 1", id_fix, group->names[igroup]))); // calculate xu,yu,zu for fix store array // skip if reset from restart file diff --git a/src/compute_vacf.cpp b/src/compute_vacf.cpp index 4e209f8612..9d6018d05b 100644 --- a/src/compute_vacf.cpp +++ b/src/compute_vacf.cpp @@ -40,7 +40,7 @@ ComputeVACF::ComputeVACF(LAMMPS *lmp, int narg, char **arg) : id_fix = utils::strdup(id + std::string("_COMPUTE_STORE")); fix = dynamic_cast( - modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 3", id_fix, group->names[igroup]))); + modify->add_fix(fmt::format("{} {} STORE/PERATOM 3 0 0 1", id_fix, group->names[igroup]))); // store current velocities in fix store array // skip if reset from restart file diff --git a/src/dump_custom.cpp b/src/dump_custom.cpp index 8d64e532a4..922fd4f6fc 100644 --- a/src/dump_custom.cpp +++ b/src/dump_custom.cpp @@ -2013,7 +2013,7 @@ int DumpCustom::modify_param(int narg, char **arg) std::string threshid = fmt::format("{}{}_DUMP_STORE",id,nthreshlast); thresh_fixID[nthreshlast] = utils::strdup(threshid); - threshid += fmt::format(" {} STORE/PERATOM 1 1", group->names[igroup]); + threshid += fmt::format(" {} STORE/PERATOM 1 0 0 1", group->names[igroup]); thresh_fix[nthreshlast] = dynamic_cast(modify->add_fix(threshid)); thresh_last[nthreshlast] = nthreshlast; diff --git a/src/fix_adapt.cpp b/src/fix_adapt.cpp index e0eaeb864e..98cf955485 100644 --- a/src/fix_adapt.cpp +++ b/src/fix_adapt.cpp @@ -279,7 +279,7 @@ void FixAdapt::post_constructor() if (diamflag && atom->radius_flag) { id_fix_diam = utils::strdup(id + std::string("_FIX_STORE_DIAM")); fix_diam = dynamic_cast( - modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 1", id_fix_diam,group->names[igroup]))); + modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 0 0 1", id_fix_diam,group->names[igroup]))); if (fix_diam->restart_reset) fix_diam->restart_reset = 0; else { double *vec = fix_diam->vstore; @@ -297,7 +297,7 @@ void FixAdapt::post_constructor() if (chgflag && atom->q_flag) { id_fix_chg = utils::strdup(id + std::string("_FIX_STORE_CHG")); fix_chg = dynamic_cast( - modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 1",id_fix_chg,group->names[igroup]))); + modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 0 0 1",id_fix_chg,group->names[igroup]))); if (fix_chg->restart_reset) fix_chg->restart_reset = 0; else { double *vec = fix_chg->vstore; diff --git a/src/fix_store_peratom.cpp b/src/fix_store_peratom.cpp index 6f5637e3f6..f9f3b5540f 100644 --- a/src/fix_store_peratom.cpp +++ b/src/fix_store_peratom.cpp @@ -22,50 +22,57 @@ using namespace LAMMPS_NS; using namespace FixConst; +// INTERNAL fix for storing/communicating per-atom values +// syntax: id group style n1 n2 gflag rflag +// N1 = 1, N2 = 0 is per-atom vector, single value per atom +// N1 > 1, N2 = 0 is per-atom array, N1 values per atom +// N1 > 0, N2 > 0 is per-atom tensor, N1xN2 array per atom +// gflag = 0/1, no/yes communicate per-atom values with ghost atoms +// rflag = 0/1, no/yes store per-atom values in restart file + /* ---------------------------------------------------------------------- */ FixStorePeratom::FixStorePeratom(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), vstore(nullptr), astore(nullptr) { - if (narg != 5 && narg != 6) + if (narg != 8) error->all(FLERR, "Illegal fix STORE/PERATOM command: number of args"); - // syntax: id group style 0/1 n1 n2 (n3), last arg optional - // 0/1 flag = not-store or store peratom values in restart file - // N2 = 1 and no n3 is vector, N2 > 1 and no n3 is array, N3 = tensor - // nvalues = # of peratom values, N = 1 is vector, N > 1 is array - disable = 0; - vecflag = arrayflag = tensorflag = 0; - restart_peratom = utils::inumeric(FLERR, arg[3], false, lmp); - n2 = utils::inumeric(FLERR, arg[4], false, lmp); - if (narg == 6) - n3 = utils::inumeric(FLERR, arg[5], false, lmp); - else - n3 = 1; - if (restart_peratom < 0 || restart_peratom > 1) - error->all(FLERR, "Illegal fix STORE/PERATOM restart flag: {}", restart_peratom); - if (n2 <= 0 || n3 <= 0) - error->all(FLERR, "Illegal fix STORE/PERATOM dimension args: must be >0: {} {}", n2, n3); - if (n2 == 1 && narg == 5) - vecflag = 1; - else if (narg == 5) - arrayflag = 1; - else - tensorflag = 1; - nvalues = n2 * n3; + n1 = utils::inumeric(FLERR, arg[4], false, lmp); + n2 = utils::inumeric(FLERR, arg[5], false, lmp); + ghostflag = utils::inumeric(FLERR, arg[6], false, lmp); + restartflag = utils::inumeric(FLERR, arg[7], false, lmp); + + vecflag = arrayflag = tensorflag = 0; + if (n1 == 1 && n2 == 0) vecflag = 1; + else if (n1 > 1 && n2 == 0) arrayflag = 1; + else if (n1 > 0 && n2 > 0) tensorflag = 1; + else error->all(FLERR, "Illegal fix STORE/PERATOM dimension args: {} {}", n1, n2); + + if (restartflag < 0 || restartflag > 1) + error->all(FLERR, "Illegal fix STORE/PERATOM restart flag: {}", restartflag); + if (ghostflag < 0 || ghostflag > 1) + error->all(FLERR, "Illegal fix STORE/PERATOM ghost flag: {}", ghostflag); + + if (vecflag || arrayflag) nvalues = n1; + else nvalues = n1 * n2; nbytes = nvalues * sizeof(double); + if (ghostflag) comm_border = nvalues; + maxexchange = nvalues; + + // allocate data structs and register with Atom class + vstore = nullptr; astore = nullptr; tstore = nullptr; - // allocate data structs and register with Atom class - FixStorePeratom::grow_arrays(atom->nmax); atom->add_callback(Atom::GROW); - if (restart_peratom) atom->add_callback(Atom::RESTART); + if (restartflag) atom->add_callback(Atom::RESTART); + if (ghostflag) atom->add_callback(Atom::BORDER); // zero the storage @@ -74,13 +81,12 @@ FixStorePeratom::FixStorePeratom(LAMMPS *lmp, int narg, char **arg) : for (int i = 0; i < nlocal; i++) vstore[i] = 0.0; } else if (arrayflag) { for (int i = 0; i < nlocal; i++) - for (int j = 0; j < n2; j++) astore[i][j] = 0.0; + for (int j = 0; j < n1; j++) astore[i][j] = 0.0; } else if (tensorflag) { for (int i = 0; i < nlocal; i++) - for (int j = 0; j < n2; j++) - for (int k = 0; k < n3; k++) tstore[i][j][k] = 0.0; + for (int j = 0; j < n1; j++) + for (int k = 0; k < n2; k++) tstore[i][j][k] = 0.0; } - maxexchange = nvalues; } /* ---------------------------------------------------------------------- */ @@ -90,7 +96,8 @@ FixStorePeratom::~FixStorePeratom() // unregister callbacks to this fix from Atom class atom->delete_callback(id, Atom::GROW); - if (restart_peratom) atom->delete_callback(id, Atom::RESTART); + if (restartflag) atom->delete_callback(id, Atom::RESTART); + if (ghostflag) atom->delete_callback(id, Atom::BORDER); memory->destroy(vstore); memory->destroy(astore); @@ -114,9 +121,9 @@ void FixStorePeratom::grow_arrays(int nmax) if (vecflag) memory->grow(vstore, nmax, "store:vstore"); else if (arrayflag) - memory->grow(astore, nmax, n2, "store:astore"); + memory->grow(astore, nmax, n1, "store:astore"); else if (tensorflag) - memory->grow(tstore, nmax, n2, n3, "store:tstore"); + memory->grow(tstore, nmax, n1, n2, "store:tstore"); } /* ---------------------------------------------------------------------- @@ -136,6 +143,62 @@ void FixStorePeratom::copy_arrays(int i, int j, int /*delflag*/) } } +/* ---------------------------------------------------------------------- + pack values for border communication at re-neighboring +------------------------------------------------------------------------- */ + +int FixStorePeratom::pack_border(int n, int *list, double *buf) +{ + int i, j, k, ncol; + + int m = 0; + if (vecflag) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = vstore[j]; + } + } else if (arrayflag) { + for (i = 0; i < n; i++) { + j = list[i]; + for (k = 0; k < nvalues; k++) + buf[m++] = astore[j][k]; + } + } else if (tensorflag) { + for (i = 0; i < n; i++) { + j = list[i]; + memcpy(&buf[m], &tstore[j][0][0], nbytes); + m += nvalues; + } + } + + return m; +} + +/* ---------------------------------------------------------------------- + unpack values for border communication at re-neighboring +------------------------------------------------------------------------- */ + +int FixStorePeratom::unpack_border(int n, int first, double *buf) +{ + int i, k, last, ncol; + + int m = 0; + last = first + n; + if (vecflag) { + for (i = first; i < last; i++) + vstore[i] = buf[m++]; + } else if (arrayflag) { + for (i = first; i < last; i++) + for (k = 0; k < nvalues; k++) + astore[i][k] = buf[m++]; + } else if (tensorflag) { + for (i = first; i < last; i++) { + memcpy(&tstore[i][0][0], &buf[m], nbytes); + m += nvalues; + } + } +} + /* ---------------------------------------------------------------------- pack values in local atom-based array for exchange with another proc ------------------------------------------------------------------------- */ @@ -147,7 +210,8 @@ int FixStorePeratom::pack_exchange(int i, double *buf) if (vecflag) { buf[0] = vstore[i]; } else if (arrayflag) { - for (int m = 0; m < nvalues; m++) buf[m] = astore[i][m]; + for (int m = 0; m < nvalues; m++) + buf[m] = astore[i][m]; } else if (tensorflag) { memcpy(buf, &tstore[i][0][0], nbytes); } @@ -166,7 +230,8 @@ int FixStorePeratom::unpack_exchange(int nlocal, double *buf) if (vecflag) { vstore[nlocal] = buf[0]; } else if (arrayflag) { - for (int m = 0; m < nvalues; m++) astore[nlocal][m] = buf[m]; + for (int m = 0; m < nvalues; m++) + astore[nlocal][m] = buf[m]; } else if (tensorflag) { memcpy(&tstore[nlocal][0][0], buf, nbytes); } @@ -191,7 +256,8 @@ int FixStorePeratom::pack_restart(int i, double *buf) if (vecflag) { buf[1] = vstore[i]; } else if (arrayflag) { - for (int m = 0; m < nvalues; m++) buf[m + 1] = astore[i][m]; + for (int m = 0; m < nvalues; m++) + buf[m + 1] = astore[i][m]; } else if (tensorflag) { memcpy(&buf[1], &tstore[i][0][0], nbytes); } @@ -219,7 +285,8 @@ void FixStorePeratom::unpack_restart(int nlocal, int nth) if (vecflag) { vstore[nlocal] = extra[nlocal][m]; } else if (arrayflag) { - for (int i = 0; i < nvalues; i++) astore[nlocal][i] = extra[nlocal][m++]; + for (int i = 0; i < nvalues; i++) + astore[nlocal][i] = extra[nlocal][m++]; } else if (tensorflag) { memcpy(&tstore[nlocal][0][0], &extra[nlocal][m], nbytes); } @@ -251,5 +318,5 @@ int FixStorePeratom::size_restart(int /*nlocal*/) double FixStorePeratom::memory_usage() { - return (double) atom->nmax * n2 * n3 * sizeof(double); + return (double) atom->nmax * nvalues * sizeof(double); } diff --git a/src/fix_store_peratom.h b/src/fix_store_peratom.h index d3efe4142b..dcf28be265 100644 --- a/src/fix_store_peratom.h +++ b/src/fix_store_peratom.h @@ -37,6 +37,8 @@ class FixStorePeratom : public Fix { void grow_arrays(int) override; void copy_arrays(int, int, int) override; + int pack_border(int, int *, double *) override; + int unpack_border(int, int, double *) override; int pack_exchange(int, double *) override; int unpack_exchange(int, double *) override; int pack_restart(int, double *) override; @@ -50,8 +52,10 @@ class FixStorePeratom : public Fix { int vecflag; // 1 if ncol=1 int arrayflag; // 1 if a 2d array (vector per atom) int tensorflag; // 1 if a 3d array (array per atom) - - int n2, n3; // size of 3d dims of per-atom data struct + int ghostflag; // 0/1 to communicate values with ghost atoms + int restartflag; // 0/1 to store values in restart files + + int n1, n2; // size of 3d dims of per-atom data struct int nvalues; // number of per-atom values int nbytes; // number of per-atom bytes }; diff --git a/src/variable.cpp b/src/variable.cpp index 2ac9bd0364..56a0604a7e 100644 --- a/src/variable.cpp +++ b/src/variable.cpp @@ -5028,7 +5028,7 @@ VarReader::VarReader(LAMMPS *lmp, char *name, char *file, int flag) : id_fix = utils::strdup(std::string(name) + "_VARIABLE_STORE"); fixstore = dynamic_cast( - modify->add_fix(std::string(id_fix) + " all STORE/PERATOM 0 1")); + modify->add_fix(std::string(id_fix) + " all STORE/PERATOM 1 0 0 0")); buffer = new char[CHUNK*MAXLINE]; } } From c369c0252f66a75b4d0036054ce6046515e477b1 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 13 Mar 2023 19:43:08 -0400 Subject: [PATCH 11/14] fix issues with changes in fix STORE/PERATOM --- src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp | 2 +- src/fix_store_peratom.cpp | 66 +++++++++---------- 2 files changed, 32 insertions(+), 36 deletions(-) diff --git a/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp b/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp index 59bab65713..b71e60c286 100644 --- a/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp +++ b/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp @@ -179,7 +179,7 @@ void BondHarmonicRestrain::init_style() if (!initial) { initial = dynamic_cast( - modify->add_fix("BOND_RESTRAIN_X0 all STORE/PERATOM 3 0 1 1")); + modify->add_fix("BOND_RESTRAIN_X0 all STORE/PERATOM 3 0 1 1")); if (!initial) error->all(FLERR, "Failure to create internal per-atom storage"); } } diff --git a/src/fix_store_peratom.cpp b/src/fix_store_peratom.cpp index f9f3b5540f..2a8429e42f 100644 --- a/src/fix_store_peratom.cpp +++ b/src/fix_store_peratom.cpp @@ -23,7 +23,7 @@ using namespace LAMMPS_NS; using namespace FixConst; // INTERNAL fix for storing/communicating per-atom values -// syntax: id group style n1 n2 gflag rflag +// syntax: id group style n1 n2 gflag rflag // N1 = 1, N2 = 0 is per-atom vector, single value per atom // N1 > 1, N2 = 0 is per-atom array, N1 values per atom // N1 > 0, N2 > 0 is per-atom tensor, N1xN2 array per atom @@ -35,34 +35,36 @@ using namespace FixConst; FixStorePeratom::FixStorePeratom(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), vstore(nullptr), astore(nullptr) { - if (narg != 8) - error->all(FLERR, "Illegal fix STORE/PERATOM command: number of args"); + if (narg != 7) error->all(FLERR, "Illegal fix STORE/PERATOM command: number of args"); disable = 0; - n1 = utils::inumeric(FLERR, arg[4], false, lmp); - n2 = utils::inumeric(FLERR, arg[5], false, lmp); - ghostflag = utils::inumeric(FLERR, arg[6], false, lmp); - restartflag = utils::inumeric(FLERR, arg[7], false, lmp); + n1 = utils::inumeric(FLERR, arg[3], false, lmp); + n2 = utils::inumeric(FLERR, arg[4], false, lmp); + ghostflag = utils::logical(FLERR, arg[5], false, lmp); + restartflag = utils::logical(FLERR, arg[6], false, lmp); vecflag = arrayflag = tensorflag = 0; - if (n1 == 1 && n2 == 0) vecflag = 1; - else if (n1 > 1 && n2 == 0) arrayflag = 1; - else if (n1 > 0 && n2 > 0) tensorflag = 1; - else error->all(FLERR, "Illegal fix STORE/PERATOM dimension args: {} {}", n1, n2); - - if (restartflag < 0 || restartflag > 1) - error->all(FLERR, "Illegal fix STORE/PERATOM restart flag: {}", restartflag); - if (ghostflag < 0 || ghostflag > 1) - error->all(FLERR, "Illegal fix STORE/PERATOM ghost flag: {}", ghostflag); + if (n1 == 1 && n2 == 0) + vecflag = 1; + else if (n1 > 1 && n2 == 0) + arrayflag = 1; + else if (n1 > 0 && n2 > 0) + tensorflag = 1; + else + error->all(FLERR, "Illegal fix STORE/PERATOM dimension args: {} {}", n1, n2); - if (vecflag || arrayflag) nvalues = n1; - else nvalues = n1 * n2; + if (vecflag || arrayflag) + nvalues = n1; + else + nvalues = n1 * n2; nbytes = nvalues * sizeof(double); if (ghostflag) comm_border = nvalues; maxexchange = nvalues; - + + if (restartflag) restart_peratom = 1; + // allocate data structs and register with Atom class vstore = nullptr; @@ -149,7 +151,7 @@ void FixStorePeratom::copy_arrays(int i, int j, int /*delflag*/) int FixStorePeratom::pack_border(int n, int *list, double *buf) { - int i, j, k, ncol; + int i, j, k; int m = 0; if (vecflag) { @@ -160,8 +162,7 @@ int FixStorePeratom::pack_border(int n, int *list, double *buf) } else if (arrayflag) { for (i = 0; i < n; i++) { j = list[i]; - for (k = 0; k < nvalues; k++) - buf[m++] = astore[j][k]; + for (k = 0; k < nvalues; k++) buf[m++] = astore[j][k]; } } else if (tensorflag) { for (i = 0; i < n; i++) { @@ -180,23 +181,22 @@ int FixStorePeratom::pack_border(int n, int *list, double *buf) int FixStorePeratom::unpack_border(int n, int first, double *buf) { - int i, k, last, ncol; + int i, k, last; int m = 0; last = first + n; if (vecflag) { - for (i = first; i < last; i++) - vstore[i] = buf[m++]; + for (i = first; i < last; i++) vstore[i] = buf[m++]; } else if (arrayflag) { for (i = first; i < last; i++) - for (k = 0; k < nvalues; k++) - astore[i][k] = buf[m++]; + for (k = 0; k < nvalues; k++) astore[i][k] = buf[m++]; } else if (tensorflag) { for (i = first; i < last; i++) { memcpy(&tstore[i][0][0], &buf[m], nbytes); m += nvalues; } } + return m; } /* ---------------------------------------------------------------------- @@ -210,8 +210,7 @@ int FixStorePeratom::pack_exchange(int i, double *buf) if (vecflag) { buf[0] = vstore[i]; } else if (arrayflag) { - for (int m = 0; m < nvalues; m++) - buf[m] = astore[i][m]; + for (int m = 0; m < nvalues; m++) buf[m] = astore[i][m]; } else if (tensorflag) { memcpy(buf, &tstore[i][0][0], nbytes); } @@ -230,8 +229,7 @@ int FixStorePeratom::unpack_exchange(int nlocal, double *buf) if (vecflag) { vstore[nlocal] = buf[0]; } else if (arrayflag) { - for (int m = 0; m < nvalues; m++) - astore[nlocal][m] = buf[m]; + for (int m = 0; m < nvalues; m++) astore[nlocal][m] = buf[m]; } else if (tensorflag) { memcpy(&tstore[nlocal][0][0], buf, nbytes); } @@ -256,8 +254,7 @@ int FixStorePeratom::pack_restart(int i, double *buf) if (vecflag) { buf[1] = vstore[i]; } else if (arrayflag) { - for (int m = 0; m < nvalues; m++) - buf[m + 1] = astore[i][m]; + for (int m = 0; m < nvalues; m++) buf[m + 1] = astore[i][m]; } else if (tensorflag) { memcpy(&buf[1], &tstore[i][0][0], nbytes); } @@ -285,8 +282,7 @@ void FixStorePeratom::unpack_restart(int nlocal, int nth) if (vecflag) { vstore[nlocal] = extra[nlocal][m]; } else if (arrayflag) { - for (int i = 0; i < nvalues; i++) - astore[nlocal][i] = extra[nlocal][m++]; + for (int i = 0; i < nvalues; i++) astore[nlocal][i] = extra[nlocal][m++]; } else if (tensorflag) { memcpy(&tstore[nlocal][0][0], &extra[nlocal][m], nbytes); } From ef67f790a2e17b6c53edc821be0ded31a065cdaa Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 13 Mar 2023 19:43:44 -0400 Subject: [PATCH 12/14] whitespace --- src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp | 2 +- src/fix_store_peratom.h | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp b/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp index b71e60c286..7f3a7519b6 100644 --- a/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp +++ b/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp @@ -176,7 +176,7 @@ void BondHarmonicRestrain::init_style() // after a restart, natoms is set but initial is a null pointer. // we add the fix, but do not initialize it. It will pull the data from the restart. - + if (!initial) { initial = dynamic_cast( modify->add_fix("BOND_RESTRAIN_X0 all STORE/PERATOM 3 0 1 1")); diff --git a/src/fix_store_peratom.h b/src/fix_store_peratom.h index dcf28be265..9434ffe799 100644 --- a/src/fix_store_peratom.h +++ b/src/fix_store_peratom.h @@ -54,7 +54,7 @@ class FixStorePeratom : public Fix { int tensorflag; // 1 if a 3d array (array per atom) int ghostflag; // 0/1 to communicate values with ghost atoms int restartflag; // 0/1 to store values in restart files - + int n1, n2; // size of 3d dims of per-atom data struct int nvalues; // number of per-atom values int nbytes; // number of per-atom bytes From 1031110d930826e6d752439fa7e428715249cbff Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 13 Mar 2023 19:55:52 -0400 Subject: [PATCH 13/14] update docs --- doc/src/bond_harmonic_restrain.rst | 19 +++++++------------ 1 file changed, 7 insertions(+), 12 deletions(-) diff --git a/doc/src/bond_harmonic_restrain.rst b/doc/src/bond_harmonic_restrain.rst index c48e68917c..c9707f5546 100644 --- a/doc/src/bond_harmonic_restrain.rst +++ b/doc/src/bond_harmonic_restrain.rst @@ -65,18 +65,13 @@ This bond style can only be used if LAMMPS was built with the EXTRA-MOLECULE package. See the :doc:`Build package ` page for more info. -This bond style uses :doc:`fix property/atom ` -internally to be able to determine the original bond lengths. This fix -will write its state to :doc:`data files `, but it is -currently not possible to read that stored data back with -:doc:`read_data ` for such an internally defined fix. You -are thus advised to use the "nofix" argument to the :doc:`write_data -command ` in order to avoid writing that section to the data -file that would cause an error when reading the data file. When -continuing from such a "nofix" data file, the reference bond lengths -will not be the original lengths but those of the current geometry. -Thus restarting is currently *only* possible with :doc:`binary restart -files `. +This bond style maintains internal data to determine the original bond +lengths :math:`r_{t=0}`. This information will be written to +:doc:`binary restart files ` but **not** to :doc:`data +files `. Thus, continuing a simulation is *only* possible +with :doc:`read_restart `. When using the :doc:`read_data +command `, the reference bond lengths :math:`r_{t=0}` will be +re-initialized from the current geometry. This bond style cannot be used with :doc:`fix shake or fix rattle `, with :doc:`fix filter/corotate `, or From 17f39d9d2cf8e5d85ba3a9e3fe9901d0b1bae501 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 13 Mar 2023 22:28:32 -0400 Subject: [PATCH 14/14] rename fix STORE/PERATOM to STORE/ATOM --- doc/src/Developer_notes.rst | 25 ++++++++++++ doc/src/Developer_updating.rst | 29 ++++++++++++++ src/AMOEBA/amoeba_induce.cpp | 2 +- src/AMOEBA/amoeba_utils.cpp | 2 +- src/AMOEBA/pair_amoeba.cpp | 26 ++++++------- src/AMOEBA/pair_amoeba.h | 6 +-- src/CORESHELL/compute_temp_cs.cpp | 6 +-- src/CORESHELL/compute_temp_cs.h | 2 +- src/EXTRA-COMPUTE/compute_hma.cpp | 8 ++-- src/EXTRA-COMPUTE/compute_hma.h | 2 +- src/EXTRA-COMPUTE/compute_msd_nongauss.cpp | 2 +- src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp | 10 ++--- src/EXTRA-MOLECULE/bond_harmonic_restrain.h | 2 +- src/FEP/fix_adapt_fep.cpp | 14 +++---- src/FEP/fix_adapt_fep.h | 2 +- src/GPU/pair_amoeba_gpu.cpp | 2 +- src/GPU/pair_hippo_gpu.cpp | 2 +- src/REPLICA/tad.cpp | 6 +-- src/REPLICA/tad.h | 18 ++++----- src/VTK/dump_vtk.cpp | 10 ++--- src/balance.cpp | 6 +-- src/balance.h | 10 ++--- src/compute_chunk_atom.cpp | 6 +-- src/compute_chunk_atom.h | 2 +- src/compute_displace_atom.cpp | 8 ++-- src/compute_displace_atom.h | 2 +- src/compute_msd.cpp | 8 ++-- src/compute_msd.h | 2 +- src/compute_vacf.cpp | 8 ++-- src/compute_vacf.h | 2 +- src/dump_custom.cpp | 10 ++--- src/dump_custom.h | 26 ++++++------- src/fix_adapt.cpp | 14 +++---- src/fix_adapt.h | 2 +- src/fix_balance.cpp | 2 +- ...x_store_peratom.cpp => fix_store_atom.cpp} | 38 +++++++++---------- src/{fix_store_peratom.h => fix_store_atom.h} | 12 +++--- src/fix_store_global.cpp | 2 +- src/variable.cpp | 6 +-- src/variable.h | 2 +- 40 files changed, 199 insertions(+), 145 deletions(-) rename src/{fix_store_peratom.cpp => fix_store_atom.cpp} (88%) rename src/{fix_store_peratom.h => fix_store_atom.h} (89%) diff --git a/doc/src/Developer_notes.rst b/doc/src/Developer_notes.rst index 92121cca15..07b289b583 100644 --- a/doc/src/Developer_notes.rst +++ b/doc/src/Developer_notes.rst @@ -11,6 +11,7 @@ Available topics are: - `Reading and parsing of text and text files`_ - `Requesting and accessing neighbor lists`_ +- `Choosing between a custom atom style, fix property/atom, and fix STORE/ATOM`_ - `Fix contributions to instantaneous energy, virial, and cumulative energy`_ - `KSpace PPPM FFT grids`_ @@ -216,6 +217,30 @@ command: neighbor->add_request(this, "delete_atoms", NeighConst::REQ_FULL); +Choosing between a custom atom style, fix property/atom, and fix STORE/ATOM +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +There are multiple ways to manage per-atom data within LAMMPS. Often +the per-atom storage is only used locally and managed by the class that +uses it. If the data has to persist between multiple time steps and +migrate with atoms when they move from sub-domain to sub-domain or +across periodic boundaries, then using a custom atom style, or :doc:`fix +property/atom `, or the internal fix STORE/ATOM are +possible options. + +- Using the atom style is usually the most programming effort and mostly + needed when the per-atom data is an integral part of the model like a + per-atom charge or diameter and thus should be part of the Atoms + section of a :doc:`data file `. + +- Fix property/atom is useful if the data is optional or should be + entered by the user, or accessed as a (named) custom property. In this + case the fix should be entered as part of the input (and not + internally) which allows to enter and store its content with data files. + +- Fix STORE/ATOM should be used when the data should be accessed internally + only and thus the fix can be created internally. + Fix contributions to instantaneous energy, virial, and cumulative energy ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ diff --git a/doc/src/Developer_updating.rst b/doc/src/Developer_updating.rst index 6089969c38..28b3712ae0 100644 --- a/doc/src/Developer_updating.rst +++ b/doc/src/Developer_updating.rst @@ -24,6 +24,7 @@ Available topics in mostly chronological order are: - `Use of "override" instead of "virtual"`_ - `Simplified and more compact neighbor list requests`_ - `Split of fix STORE into fix STORE/GLOBAL and fix STORE/PERATOM`_ +- `Rename of fix STORE/PERATOM to fix STORE/ATOM and change of arguments`_ - `Use Output::get_dump_by_id() instead of Output::find_dump()`_ - `Refactored grid communication using Grid3d/Grid2d classes instead of GridComm`_ @@ -385,6 +386,34 @@ New: This change is **required** or else the code will not compile. +Rename of fix STORE/PERATOM to fix STORE/ATOM and change of arguments +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +.. versionchanged:: TBD + +The available functionality of the internal fix to store per-atom +properties was expanded to enable storing data with ghost atoms and to +support binary restart files. With those changes, the fix was renamed +to fix STORE/ATOM and the number and order of (required) arguments has +changed. + +Old syntax: ``ID group-ID STORE/PERATOM rflag n1 n2 [n3]`` + +- *rflag* = 0/1, *no*/*yes* store per-atom values in restart file +- :math:`n1 = 1, n2 = 1, \mathrm{no}\;n3 \to` per-atom vector, single value per atom +- :math:`n1 = 1, n2 > 1, \mathrm{no}\;n3 \to` per-atom array, *n2* values per atom +- :math:`n1 = 1, n2 > 0, n3 > 0 \to` per-atom tensor, *n2* x *n3* values per atom + +New syntax: ``ID group-ID STORE/ATOM n1 n2 gflag rflag`` + +- :math:`n1 = 1, n2 = 0 \to` per-atom vector, single value per atom +- :math:`n1 > 1, n2 = 0 \to` per-atom array, *n1* values per atom +- :math:`n1 > 0, n2 > 0 \to` per-atom tensor, *n1* x *n2* values per atom +- *gflag* = 0/1, *no*/*yes* communicate per-atom values with ghost atoms +- *rflag* = 0/1, *no*/*yes* store per-atom values in restart file + +Since this is an internal fix, there is no user visible change. + Use Output::get_dump_by_id() instead of Output::find_dump() ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ diff --git a/src/AMOEBA/amoeba_induce.cpp b/src/AMOEBA/amoeba_induce.cpp index ecc20a198c..6017b775ca 100644 --- a/src/AMOEBA/amoeba_induce.cpp +++ b/src/AMOEBA/amoeba_induce.cpp @@ -19,7 +19,7 @@ #include "comm.h" #include "domain.h" #include "error.h" -#include "fix_store_peratom.h" +#include "fix_store_atom.h" #include "math_const.h" #include "math_special.h" #include "my_page.h" diff --git a/src/AMOEBA/amoeba_utils.cpp b/src/AMOEBA/amoeba_utils.cpp index 5a4057930c..8fb839b693 100644 --- a/src/AMOEBA/amoeba_utils.cpp +++ b/src/AMOEBA/amoeba_utils.cpp @@ -17,7 +17,7 @@ #include "atom.h" #include "domain.h" #include "error.h" -#include "fix_store_peratom.h" +#include "fix_store_atom.h" #include "neigh_list.h" #include diff --git a/src/AMOEBA/pair_amoeba.cpp b/src/AMOEBA/pair_amoeba.cpp index 319bbd8925..72efa76523 100644 --- a/src/AMOEBA/pair_amoeba.cpp +++ b/src/AMOEBA/pair_amoeba.cpp @@ -20,7 +20,7 @@ #include "domain.h" #include "error.h" #include "fix.h" -#include "fix_store_peratom.h" +#include "fix_store_atom.h" #include "force.h" #include "group.h" #include "math_special.h" @@ -861,8 +861,8 @@ void PairAmoeba::init_style() Fix *myfix; if (first_flag) { id_pole = utils::strdup("AMOEBA_pole"); - myfix = modify->add_fix(fmt::format("{} {} STORE/PERATOM 13 0 0 1",id_pole,group->names[0])); - fixpole = dynamic_cast(myfix); + myfix = modify->add_fix(fmt::format("{} {} STORE/ATOM 13 0 0 1",id_pole,group->names[0])); + fixpole = dynamic_cast(myfix); } // creation of per-atom storage @@ -873,14 +873,14 @@ void PairAmoeba::init_style() if (first_flag && use_pred) { id_udalt = utils::strdup("AMOEBA_udalt"); - myfix = modify->add_fix(fmt::format("{} {} STORE/PERATOM {} 3 0 1", + myfix = modify->add_fix(fmt::format("{} {} STORE/ATOM {} 3 0 1", id_udalt, group->names[0], maxualt)); - fixudalt = dynamic_cast(myfix); + fixudalt = dynamic_cast(myfix); id_upalt = utils::strdup("AMOEBA_upalt"); - myfix = modify->add_fix(fmt::format("{} {} STORE/PERATOM {} 3 0 1", + myfix = modify->add_fix(fmt::format("{} {} STORE/ATOM {} 3 0 1", id_upalt, group->names[0], maxualt)); - fixupalt = dynamic_cast(myfix); + fixupalt = dynamic_cast(myfix); } // create pages for storing pairwise data: @@ -995,21 +995,21 @@ void PairAmoeba::init_style() if (id_pole) { myfix = modify->get_fix_by_id(id_pole); if (!myfix) - error->all(FLERR,"Could not find internal pair amoeba fix STORE/PERATOM id {}", id_pole); - fixpole = dynamic_cast(myfix); + error->all(FLERR,"Could not find internal pair amoeba fix STORE/ATOM id {}", id_pole); + fixpole = dynamic_cast(myfix); } if (id_udalt) { myfix = modify->get_fix_by_id(id_udalt); if (!myfix) - error->all(FLERR,"Could not find internal pair amoeba fix STORE/PERATOM id {}", id_udalt); - fixudalt = dynamic_cast(myfix); + error->all(FLERR,"Could not find internal pair amoeba fix STORE/ATOM id {}", id_udalt); + fixudalt = dynamic_cast(myfix); myfix = modify->get_fix_by_id(id_upalt); if (!myfix) - error->all(FLERR,"Could not find internal pair amoeba fix STORE/PERATOM id {}", id_upalt); - fixupalt = dynamic_cast(myfix); + error->all(FLERR,"Could not find internal pair amoeba fix STORE/ATOM id {}", id_upalt); + fixupalt = dynamic_cast(myfix); } // assign hydrogen neighbors (redID) to each owned atom diff --git a/src/AMOEBA/pair_amoeba.h b/src/AMOEBA/pair_amoeba.h index cdeee6c95f..1f3a4b799a 100644 --- a/src/AMOEBA/pair_amoeba.h +++ b/src/AMOEBA/pair_amoeba.h @@ -166,9 +166,9 @@ class PairAmoeba : public Pair { int *amgroup; // AMOEBA polarization group, 1 to Ngroup char *id_pole, *id_udalt, *id_upalt; - class FixStorePeratom *fixpole; // stores pole = multipole components - class FixStorePeratom *fixudalt; // stores udalt = induced dipole history - class FixStorePeratom *fixupalt; // stores upalt = induced dipole history + class FixStoreAtom *fixpole; // stores pole = multipole components + class FixStoreAtom *fixudalt; // stores udalt = induced dipole history + class FixStoreAtom *fixupalt; // stores upalt = induced dipole history // static per-type properties defined in force-field file diff --git a/src/CORESHELL/compute_temp_cs.cpp b/src/CORESHELL/compute_temp_cs.cpp index b0ec39e542..b66db30bc5 100644 --- a/src/CORESHELL/compute_temp_cs.cpp +++ b/src/CORESHELL/compute_temp_cs.cpp @@ -24,7 +24,7 @@ #include "comm.h" #include "domain.h" #include "error.h" -#include "fix_store_peratom.h" +#include "fix_store_atom.h" #include "force.h" #include "group.h" #include "memory.h" @@ -67,8 +67,8 @@ ComputeTempCS::ComputeTempCS(LAMMPS *lmp, int narg, char **arg) : // id = compute-ID + COMPUTE_STORE, fix group = compute group id_fix = utils::strdup(id + std::string("_COMPUTE_STORE")); - fix = dynamic_cast( - modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 0 0 0", id_fix, group->names[igroup]))); + fix = dynamic_cast( + modify->add_fix(fmt::format("{} {} STORE/ATOM 1 0 0 0", id_fix, group->names[igroup]))); // set fix store values = 0 for now // fill them in via setup() once Comm::borders() has been called diff --git a/src/CORESHELL/compute_temp_cs.h b/src/CORESHELL/compute_temp_cs.h index 3fdf0f3711..1e7c537a83 100644 --- a/src/CORESHELL/compute_temp_cs.h +++ b/src/CORESHELL/compute_temp_cs.h @@ -54,7 +54,7 @@ class ComputeTempCS : public Compute { double **vint; char *id_fix; - class FixStorePeratom *fix; + class FixStoreAtom *fix; void dof_compute(); void vcm_pairs(); diff --git a/src/EXTRA-COMPUTE/compute_hma.cpp b/src/EXTRA-COMPUTE/compute_hma.cpp index 806228faef..4b949f6592 100644 --- a/src/EXTRA-COMPUTE/compute_hma.cpp +++ b/src/EXTRA-COMPUTE/compute_hma.cpp @@ -54,7 +54,7 @@ https://doi.org/10.1103/PhysRevE.92.043303 #include "domain.h" #include "error.h" #include "fix.h" -#include "fix_store_peratom.h" +#include "fix_store_atom.h" #include "force.h" #include "group.h" #include "improper.h" @@ -90,8 +90,8 @@ ComputeHMA::ComputeHMA(LAMMPS *lmp, int narg, char **arg) : // our new fix's group = same as compute group id_fix = utils::strdup(std::string(id)+"_COMPUTE_STORE"); - fix = dynamic_cast( - modify->add_fix(fmt::format("{} {} STORE/PERATOM 3 0 0 1", id_fix, group->names[igroup]))); + fix = dynamic_cast( + modify->add_fix(fmt::format("{} {} STORE/ATOM 3 0 0 1", id_fix, group->names[igroup]))); // calculate xu,yu,zu for fix store array // skip if reset from restart file @@ -196,7 +196,7 @@ void ComputeHMA::setup() // set fix which stores original atom coords - fix = dynamic_cast(modify->get_fix_by_id(id_fix)); + fix = dynamic_cast(modify->get_fix_by_id(id_fix)); if (!fix) error->all(FLERR,"Could not find hma per-atom store fix ID {}", id_fix); } diff --git a/src/EXTRA-COMPUTE/compute_hma.h b/src/EXTRA-COMPUTE/compute_hma.h index d8b8e55b96..8b5ea1a8dc 100644 --- a/src/EXTRA-COMPUTE/compute_hma.h +++ b/src/EXTRA-COMPUTE/compute_hma.h @@ -43,7 +43,7 @@ class ComputeHMA : public Compute { char *id_fix; char *id_temp; double finaltemp; - class FixStorePeratom *fix; + class FixStoreAtom *fix; double boltz, nktv2p, inv_volume; double deltaPcap; double virial_compute(int); diff --git a/src/EXTRA-COMPUTE/compute_msd_nongauss.cpp b/src/EXTRA-COMPUTE/compute_msd_nongauss.cpp index f75ab69a01..42af0f0b6e 100644 --- a/src/EXTRA-COMPUTE/compute_msd_nongauss.cpp +++ b/src/EXTRA-COMPUTE/compute_msd_nongauss.cpp @@ -19,7 +19,7 @@ #include "atom.h" #include "domain.h" -#include "fix_store_peratom.h" +#include "fix_store_atom.h" #include "group.h" #include "update.h" diff --git a/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp b/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp index 7f3a7519b6..87b6e5a80e 100644 --- a/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp +++ b/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp @@ -17,7 +17,7 @@ #include "comm.h" #include "domain.h" #include "error.h" -#include "fix_store_peratom.h" +#include "fix_store_atom.h" #include "force.h" #include "memory.h" #include "modify.h" @@ -162,8 +162,8 @@ void BondHarmonicRestrain::init_style() if (natoms < 0) { // create internal fix to store initial positions - initial = dynamic_cast( - modify->add_fix("BOND_RESTRAIN_X0 all STORE/PERATOM 3 0 1 1")); + initial = dynamic_cast( + modify->add_fix("BOND_RESTRAIN_X0 all STORE/ATOM 3 0 1 1")); if (!initial) error->all(FLERR, "Failure to create internal per-atom storage"); natoms = atom->natoms; @@ -178,8 +178,8 @@ void BondHarmonicRestrain::init_style() // we add the fix, but do not initialize it. It will pull the data from the restart. if (!initial) { - initial = dynamic_cast( - modify->add_fix("BOND_RESTRAIN_X0 all STORE/PERATOM 3 0 1 1")); + initial = dynamic_cast( + modify->add_fix("BOND_RESTRAIN_X0 all STORE/ATOM 3 0 1 1")); if (!initial) error->all(FLERR, "Failure to create internal per-atom storage"); } } diff --git a/src/EXTRA-MOLECULE/bond_harmonic_restrain.h b/src/EXTRA-MOLECULE/bond_harmonic_restrain.h index 61d01ab9b2..66ee1300c1 100644 --- a/src/EXTRA-MOLECULE/bond_harmonic_restrain.h +++ b/src/EXTRA-MOLECULE/bond_harmonic_restrain.h @@ -41,7 +41,7 @@ class BondHarmonicRestrain : public Bond { protected: double *k; bigint natoms; - class FixStorePeratom *initial; + class FixStoreAtom *initial; virtual void allocate(); }; diff --git a/src/FEP/fix_adapt_fep.cpp b/src/FEP/fix_adapt_fep.cpp index 3f3d26c4d9..6e38c90ef2 100644 --- a/src/FEP/fix_adapt_fep.cpp +++ b/src/FEP/fix_adapt_fep.cpp @@ -20,7 +20,7 @@ #include "atom.h" #include "error.h" -#include "fix_store_peratom.h" +#include "fix_store_atom.h" #include "force.h" #include "group.h" #include "input.h" @@ -211,8 +211,8 @@ void FixAdaptFEP::post_constructor() if (diamflag) { id_fix_diam = utils::strdup(id + std::string("_FIX_STORE_DIAM")); - fix_diam = dynamic_cast( - modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 0 0 1", id_fix_diam,group->names[igroup]))); + fix_diam = dynamic_cast( + modify->add_fix(fmt::format("{} {} STORE/ATOM 1 0 0 1", id_fix_diam,group->names[igroup]))); if (fix_diam->restart_reset) fix_diam->restart_reset = 0; else { double *vec = fix_diam->vstore; @@ -229,8 +229,8 @@ void FixAdaptFEP::post_constructor() if (chgflag) { id_fix_chg = utils::strdup(id + std::string("_FIX_STORE_CHG")); - fix_chg = dynamic_cast( - modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 0 0 1",id_fix_chg,group->names[igroup]))); + fix_chg = dynamic_cast( + modify->add_fix(fmt::format("{} {} STORE/ATOM 1 0 0 1",id_fix_chg,group->names[igroup]))); if (fix_chg->restart_reset) fix_chg->restart_reset = 0; else { double *vec = fix_chg->vstore; @@ -333,11 +333,11 @@ void FixAdaptFEP::init() // fixes that store initial per-atom values if (id_fix_diam) { - fix_diam = dynamic_cast(modify->get_fix_by_id(id_fix_diam)); + fix_diam = dynamic_cast(modify->get_fix_by_id(id_fix_diam)); if (!fix_diam) error->all(FLERR,"Could not find fix adapt/fep storage fix ID {}", id_fix_diam); } if (id_fix_chg) { - fix_chg = dynamic_cast(modify->get_fix_by_id(id_fix_chg)); + fix_chg = dynamic_cast(modify->get_fix_by_id(id_fix_chg)); if (!fix_chg) error->all(FLERR,"Could not find fix adapt/fep storage fix ID {}", id_fix_chg); } diff --git a/src/FEP/fix_adapt_fep.h b/src/FEP/fix_adapt_fep.h index a93349a780..370a34a0f3 100644 --- a/src/FEP/fix_adapt_fep.h +++ b/src/FEP/fix_adapt_fep.h @@ -46,7 +46,7 @@ class FixAdaptFEP : public Fix { int anypair; int nlevels_respa; char *id_fix_diam, *id_fix_chg; - class FixStorePeratom *fix_diam, *fix_chg; + class FixStoreAtom *fix_diam, *fix_chg; struct Adapt { int which, ivar; diff --git a/src/GPU/pair_amoeba_gpu.cpp b/src/GPU/pair_amoeba_gpu.cpp index fd423486fd..c513d32081 100644 --- a/src/GPU/pair_amoeba_gpu.cpp +++ b/src/GPU/pair_amoeba_gpu.cpp @@ -23,7 +23,7 @@ #include "comm.h" #include "domain.h" #include "error.h" -#include "fix_store_peratom.h" +#include "fix_store_atom.h" #include "force.h" #include "gpu_extra.h" #include "info.h" diff --git a/src/GPU/pair_hippo_gpu.cpp b/src/GPU/pair_hippo_gpu.cpp index 9d286d5db7..aa10ad6c5f 100644 --- a/src/GPU/pair_hippo_gpu.cpp +++ b/src/GPU/pair_hippo_gpu.cpp @@ -23,7 +23,7 @@ #include "comm.h" #include "domain.h" #include "error.h" -#include "fix_store_peratom.h" +#include "fix_store_atom.h" #include "force.h" #include "gpu_extra.h" #include "info.h" diff --git a/src/REPLICA/tad.cpp b/src/REPLICA/tad.cpp index e98cd39233..b9b04a63d5 100644 --- a/src/REPLICA/tad.cpp +++ b/src/REPLICA/tad.cpp @@ -24,7 +24,7 @@ #include "error.h" #include "finish.h" #include "fix_event_tad.h" -#include "fix_store_peratom.h" +#include "fix_store_atom.h" #include "force.h" #include "integrate.h" #include "memory.h" @@ -133,9 +133,9 @@ void TAD::command(int narg, char **arg) fix_event = dynamic_cast(modify->add_fix("tad_event all EVENT/TAD")); - // create FixStorePeratom object to store revert state + // create FixStoreAtom object to store revert state - fix_revert = dynamic_cast(modify->add_fix("tad_revert all STORE/PERATOM 7 0 0 0")); + fix_revert = dynamic_cast(modify->add_fix("tad_revert all STORE/ATOM 7 0 0 0")); // create Finish for timing output diff --git a/src/REPLICA/tad.h b/src/REPLICA/tad.h index 8f51ca8c91..88379794cc 100644 --- a/src/REPLICA/tad.h +++ b/src/REPLICA/tad.h @@ -48,15 +48,15 @@ class TAD : public Command { double time_dynamics, time_quench, time_neb, time_comm, time_output; double time_start; - class NEB *neb; // NEB object - class Fix *fix_neb; // FixNEB object - class Compute *compute_event; // compute to detect event - class FixEventTAD *fix_event; // current event/state - class FixStorePeratom *fix_revert; // revert state - FixEventTAD **fix_event_list; // list of possible events - int n_event_list; // number of events - int nmax_event_list; // allocated events - int nmin_event_list; // minimum allocation + class NEB *neb; // NEB object + class Fix *fix_neb; // FixNEB object + class Compute *compute_event; // compute to detect event + class FixEventTAD *fix_event; // current event/state + class FixStoreAtom *fix_revert; // revert state + FixEventTAD **fix_event_list; // list of possible events + int n_event_list; // number of events + int nmax_event_list; // allocated events + int nmin_event_list; // minimum allocation char *neb_logfilename; // filename for ulogfile_neb FILE *uscreen_neb; // neb universe screen output diff --git a/src/VTK/dump_vtk.cpp b/src/VTK/dump_vtk.cpp index c4dd68bf3c..75033684ae 100644 --- a/src/VTK/dump_vtk.cpp +++ b/src/VTK/dump_vtk.cpp @@ -31,7 +31,7 @@ #include "domain.h" #include "error.h" #include "fix.h" -#include "fix_store_peratom.h" +#include "fix_store_atom.h" #include "force.h" #include "group.h" #include "input.h" @@ -2357,16 +2357,16 @@ int DumpVTK::modify_param(int narg, char **arg) thresh_value[nthresh] = utils::numeric(FLERR,arg[3],false,lmp); thresh_last[nthresh] = -1; } else { - thresh_fix = (FixStorePeratom **) - memory->srealloc(thresh_fix,(nthreshlast+1)*sizeof(FixStorePeratom *),"dump:thresh_fix"); + thresh_fix = (FixStoreAtom **) + memory->srealloc(thresh_fix,(nthreshlast+1)*sizeof(FixStoreAtom *),"dump:thresh_fix"); thresh_fixID = (char **) memory->srealloc(thresh_fixID,(nthreshlast+1)*sizeof(char *),"dump:thresh_fixID"); memory->grow(thresh_first,(nthreshlast+1),"dump:thresh_first"); std::string threshid = fmt::format("{}{}_DUMP_STORE",id,nthreshlast); thresh_fixID[nthreshlast] = utils::strdup(threshid); - threshid += fmt::format(" {} STORE/PERATOM 1 0 0 1", group->names[igroup]); - thresh_fix[nthreshlast] = dynamic_cast(modify->add_fix(threshid)); + threshid += fmt::format(" {} STORE/ATOM 1 0 0 1", group->names[igroup]); + thresh_fix[nthreshlast] = dynamic_cast(modify->add_fix(threshid)); thresh_last[nthreshlast] = nthreshlast; thresh_first[nthreshlast] = 1; diff --git a/src/balance.cpp b/src/balance.cpp index 8d1ba460db..3bd083e2b9 100644 --- a/src/balance.cpp +++ b/src/balance.cpp @@ -26,7 +26,7 @@ #include "neighbor.h" #include "comm.h" #include "domain.h" -#include "fix_store_peratom.h" +#include "fix_store_atom.h" #include "force.h" #include "imbalance.h" #include "imbalance_group.h" @@ -513,9 +513,9 @@ void Balance::weight_storage(char *prefix) if (prefix) cmd = prefix; cmd += "IMBALANCE_WEIGHTS"; - fixstore = dynamic_cast(modify->get_fix_by_id(cmd)); + fixstore = dynamic_cast(modify->get_fix_by_id(cmd)); if (!fixstore) - fixstore = dynamic_cast(modify->add_fix(cmd + " all STORE/PERATOM 1 0 0 0")); + fixstore = dynamic_cast(modify->add_fix(cmd + " all STORE/ATOM 1 0 0 0")); // do not carry weights with atoms during normal atom migration diff --git a/src/balance.h b/src/balance.h index 3a6e8bf457..669d21e2b5 100644 --- a/src/balance.h +++ b/src/balance.h @@ -27,11 +27,11 @@ namespace LAMMPS_NS { class Balance : public Command { public: class RCB *rcb; - class FixStorePeratom *fixstore; // per-atom weights stored in FixStorePeratom - int wtflag; // 1 if particle weighting is used - int varflag; // 1 if weight style var(iable) is used - int sortflag; // 1 if sorting of comm messages is done - int outflag; // 1 for output of balance results to file + class FixStoreAtom *fixstore; // per-atom weights stored in FixStorePeratom + int wtflag; // 1 if particle weighting is used + int varflag; // 1 if weight style var(iable) is used + int sortflag; // 1 if sorting of comm messages is done + int outflag; // 1 for output of balance results to file Balance(class LAMMPS *); ~Balance() override; diff --git a/src/compute_chunk_atom.cpp b/src/compute_chunk_atom.cpp index b59449989b..f5c68c8f25 100644 --- a/src/compute_chunk_atom.cpp +++ b/src/compute_chunk_atom.cpp @@ -21,7 +21,7 @@ #include "domain.h" #include "error.h" #include "fix.h" -#include "fix_store_peratom.h" +#include "fix_store_atom.h" #include "group.h" #include "input.h" #include "lattice.h" @@ -571,8 +571,8 @@ void ComputeChunkAtom::init() if ((idsflag == ONCE || lockcount) && !fixstore) { id_fix = utils::strdup(id + std::string("_COMPUTE_STORE")); - fixstore = dynamic_cast( - modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 0 0 1", id_fix, group->names[igroup]))); + fixstore = dynamic_cast( + modify->add_fix(fmt::format("{} {} STORE/ATOM 1 0 0 1", id_fix, group->names[igroup]))); } if ((idsflag != ONCE && !lockcount) && fixstore) { diff --git a/src/compute_chunk_atom.h b/src/compute_chunk_atom.h index 6850003b3e..943e36c985 100644 --- a/src/compute_chunk_atom.h +++ b/src/compute_chunk_atom.h @@ -93,7 +93,7 @@ class ComputeChunkAtom : public Compute { double *varatom; char *id_fix; - class FixStorePeratom *fixstore; + class FixStoreAtom *fixstore; class Fix *lockfix; // ptr to FixAveChunk that is locking out setups // null pointer if no lock currently in place diff --git a/src/compute_displace_atom.cpp b/src/compute_displace_atom.cpp index 6f352fa111..9f37ff2a78 100644 --- a/src/compute_displace_atom.cpp +++ b/src/compute_displace_atom.cpp @@ -17,7 +17,7 @@ #include "atom.h" #include "domain.h" #include "error.h" -#include "fix_store_peratom.h" +#include "fix_store_atom.h" #include "group.h" #include "input.h" #include "memory.h" @@ -74,8 +74,8 @@ ComputeDisplaceAtom::ComputeDisplaceAtom(LAMMPS *lmp, int narg, char **arg) : // id = compute-ID + COMPUTE_STORE, fix group = compute group id_fix = utils::strdup(std::string(id) + "_COMPUTE_STORE"); - fix = dynamic_cast( - modify->add_fix(fmt::format("{} {} STORE/PERATOM 3 0 0 1", id_fix, group->names[igroup]))); + fix = dynamic_cast( + modify->add_fix(fmt::format("{} {} STORE/ATOM 3 0 0 1", id_fix, group->names[igroup]))); // calculate xu,yu,zu for fix store array // skip if reset from restart file @@ -120,7 +120,7 @@ void ComputeDisplaceAtom::init() { // set fix which stores original atom coords - fix = dynamic_cast(modify->get_fix_by_id(id_fix)); + fix = dynamic_cast(modify->get_fix_by_id(id_fix)); if (!fix) error->all(FLERR,"Could not find compute displace/atom fix with ID {}", id_fix); if (refreshflag) { diff --git a/src/compute_displace_atom.h b/src/compute_displace_atom.h index dc1be7434f..006928bd42 100644 --- a/src/compute_displace_atom.h +++ b/src/compute_displace_atom.h @@ -38,7 +38,7 @@ class ComputeDisplaceAtom : public Compute { int nmax; double **displace; char *id_fix; - class FixStorePeratom *fix; + class FixStoreAtom *fix; int refreshflag, ivar, nvmax; // refresh option is enabled char *rvar; // for incremental dumps diff --git a/src/compute_msd.cpp b/src/compute_msd.cpp index ecbb540a48..e73dbd3d53 100644 --- a/src/compute_msd.cpp +++ b/src/compute_msd.cpp @@ -16,7 +16,7 @@ #include "atom.h" #include "domain.h" #include "error.h" -#include "fix_store_peratom.h" +#include "fix_store_atom.h" #include "group.h" #include "modify.h" #include "update.h" @@ -63,8 +63,8 @@ ComputeMSD::ComputeMSD(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, a // id = compute-ID + COMPUTE_STORE, fix group = compute group id_fix = utils::strdup(id + std::string("_COMPUTE_STORE")); - fix = dynamic_cast( - modify->add_fix(fmt::format("{} {} STORE/PERATOM 3 0 0 1", id_fix, group->names[igroup]))); + fix = dynamic_cast( + modify->add_fix(fmt::format("{} {} STORE/ATOM 3 0 0 1", id_fix, group->names[igroup]))); // calculate xu,yu,zu for fix store array // skip if reset from restart file @@ -127,7 +127,7 @@ void ComputeMSD::init() { // set fix which stores reference atom coords - fix = dynamic_cast(modify->get_fix_by_id(id_fix)); + fix = dynamic_cast(modify->get_fix_by_id(id_fix)); if (!fix) error->all(FLERR, "Could not find compute msd fix with ID {}", id_fix); // nmsd = # of atoms in group diff --git a/src/compute_msd.h b/src/compute_msd.h index a273e9cf58..ffb0c18c8c 100644 --- a/src/compute_msd.h +++ b/src/compute_msd.h @@ -39,7 +39,7 @@ class ComputeMSD : public Compute { bigint nmsd; double masstotal; char *id_fix; - class FixStorePeratom *fix; + class FixStoreAtom *fix; }; } // namespace LAMMPS_NS diff --git a/src/compute_vacf.cpp b/src/compute_vacf.cpp index 9d6018d05b..acf1405092 100644 --- a/src/compute_vacf.cpp +++ b/src/compute_vacf.cpp @@ -18,7 +18,7 @@ #include "update.h" #include "group.h" #include "modify.h" -#include "fix_store_peratom.h" +#include "fix_store_atom.h" #include "error.h" using namespace LAMMPS_NS; @@ -39,8 +39,8 @@ ComputeVACF::ComputeVACF(LAMMPS *lmp, int narg, char **arg) : // id = compute-ID + COMPUTE_STORE, fix group = compute group id_fix = utils::strdup(id + std::string("_COMPUTE_STORE")); - fix = dynamic_cast( - modify->add_fix(fmt::format("{} {} STORE/PERATOM 3 0 0 1", id_fix, group->names[igroup]))); + fix = dynamic_cast( + modify->add_fix(fmt::format("{} {} STORE/ATOM 3 0 0 1", id_fix, group->names[igroup]))); // store current velocities in fix store array // skip if reset from restart file @@ -84,7 +84,7 @@ void ComputeVACF::init() { // set fix which stores original atom velocities - fix = dynamic_cast(modify->get_fix_by_id(id_fix)); + fix = dynamic_cast(modify->get_fix_by_id(id_fix)); if (!fix) error->all(FLERR,"Could not find compute vacf fix ID {}", id_fix); // nvacf = # of atoms in group diff --git a/src/compute_vacf.h b/src/compute_vacf.h index 42b3ed4adc..8e5e57ab65 100644 --- a/src/compute_vacf.h +++ b/src/compute_vacf.h @@ -35,7 +35,7 @@ class ComputeVACF : public Compute { protected: bigint nvacf; char *id_fix; - class FixStorePeratom *fix; + class FixStoreAtom *fix; }; } // namespace LAMMPS_NS diff --git a/src/dump_custom.cpp b/src/dump_custom.cpp index 922fd4f6fc..8667c6af12 100644 --- a/src/dump_custom.cpp +++ b/src/dump_custom.cpp @@ -20,7 +20,7 @@ #include "domain.h" #include "error.h" #include "fix.h" -#include "fix_store_peratom.h" +#include "fix_store_atom.h" #include "group.h" #include "input.h" #include "memory.h" @@ -2005,16 +2005,16 @@ int DumpCustom::modify_param(int narg, char **arg) thresh_value[nthresh] = utils::numeric(FLERR,arg[3],false,lmp); thresh_last[nthresh] = -1; } else { - thresh_fix = (FixStorePeratom **) - memory->srealloc(thresh_fix,(nthreshlast+1)*sizeof(FixStorePeratom *),"dump:thresh_fix"); + thresh_fix = (FixStoreAtom **) + memory->srealloc(thresh_fix,(nthreshlast+1)*sizeof(FixStoreAtom *),"dump:thresh_fix"); thresh_fixID = (char **) memory->srealloc(thresh_fixID,(nthreshlast+1)*sizeof(char *),"dump:thresh_fixID"); memory->grow(thresh_first,(nthreshlast+1),"dump:thresh_first"); std::string threshid = fmt::format("{}{}_DUMP_STORE",id,nthreshlast); thresh_fixID[nthreshlast] = utils::strdup(threshid); - threshid += fmt::format(" {} STORE/PERATOM 1 0 0 1", group->names[igroup]); - thresh_fix[nthreshlast] = dynamic_cast(modify->add_fix(threshid)); + threshid += fmt::format(" {} STORE/ATOM 1 0 0 1", group->names[igroup]); + thresh_fix[nthreshlast] = dynamic_cast(modify->add_fix(threshid)); thresh_last[nthreshlast] = nthreshlast; thresh_first[nthreshlast] = 1; diff --git a/src/dump_custom.h b/src/dump_custom.h index da49ffdc22..b600bd60b8 100644 --- a/src/dump_custom.h +++ b/src/dump_custom.h @@ -37,19 +37,19 @@ class DumpCustom : public Dump { int nevery; // dump frequency for output char *idregion; // region ID, nullptr if no region - int nthresh; // # of defined thresholds - int nthreshlast; // # of defined thresholds with value = LAST - // - int *thresh_array; // array to threshold on for each nthresh - int *thresh_op; // threshold operation for each nthresh - double *thresh_value; // threshold value for each nthresh - int *thresh_last; // for threshold value = LAST, - // index into thresh_fix - // -1 if not LAST, value is numeric - // - class FixStorePeratom **thresh_fix; // stores values for each threshold LAST - char **thresh_fixID; // IDs of thresh_fixes - int *thresh_first; // 1 the first time a FixStore values accessed + int nthresh; // # of defined thresholds + int nthreshlast; // # of defined thresholds with value = LAST + // + int *thresh_array; // array to threshold on for each nthresh + int *thresh_op; // threshold operation for each nthresh + double *thresh_value; // threshold value for each nthresh + int *thresh_last; // for threshold value = LAST, + // index into thresh_fix + // -1 if not LAST, value is numeric + // + class FixStoreAtom **thresh_fix; // stores values for each threshold LAST + char **thresh_fixID; // IDs of thresh_fixes + int *thresh_first; // 1 the first time a FixStore values accessed int expand; // flag for whether field args were expanded char **earg; // field names with wildcard expansion diff --git a/src/fix_adapt.cpp b/src/fix_adapt.cpp index 98cf955485..8dd97250a3 100644 --- a/src/fix_adapt.cpp +++ b/src/fix_adapt.cpp @@ -19,7 +19,7 @@ #include "bond.h" #include "domain.h" #include "error.h" -#include "fix_store_peratom.h" +#include "fix_store_atom.h" #include "force.h" #include "group.h" #include "input.h" @@ -278,8 +278,8 @@ void FixAdapt::post_constructor() if (diamflag && atom->radius_flag) { id_fix_diam = utils::strdup(id + std::string("_FIX_STORE_DIAM")); - fix_diam = dynamic_cast( - modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 0 0 1", id_fix_diam,group->names[igroup]))); + fix_diam = dynamic_cast( + modify->add_fix(fmt::format("{} {} STORE/ATOM 1 0 0 1", id_fix_diam,group->names[igroup]))); if (fix_diam->restart_reset) fix_diam->restart_reset = 0; else { double *vec = fix_diam->vstore; @@ -296,8 +296,8 @@ void FixAdapt::post_constructor() if (chgflag && atom->q_flag) { id_fix_chg = utils::strdup(id + std::string("_FIX_STORE_CHG")); - fix_chg = dynamic_cast( - modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 0 0 1",id_fix_chg,group->names[igroup]))); + fix_chg = dynamic_cast( + modify->add_fix(fmt::format("{} {} STORE/ATOM 1 0 0 1",id_fix_chg,group->names[igroup]))); if (fix_chg->restart_reset) fix_chg->restart_reset = 0; else { double *vec = fix_chg->vstore; @@ -494,11 +494,11 @@ void FixAdapt::init() // fixes that store initial per-atom values if (id_fix_diam) { - fix_diam = dynamic_cast(modify->get_fix_by_id(id_fix_diam)); + fix_diam = dynamic_cast(modify->get_fix_by_id(id_fix_diam)); if (!fix_diam) error->all(FLERR,"Could not find fix adapt storage fix ID {}", id_fix_diam); } if (id_fix_chg) { - fix_chg = dynamic_cast(modify->get_fix_by_id(id_fix_chg)); + fix_chg = dynamic_cast(modify->get_fix_by_id(id_fix_chg)); if (!fix_chg) error->all(FLERR,"Could not find fix adapt storage fix ID {}", id_fix_chg); } diff --git a/src/fix_adapt.h b/src/fix_adapt.h index 82dc739c67..f8477f7259 100644 --- a/src/fix_adapt.h +++ b/src/fix_adapt.h @@ -48,7 +48,7 @@ class FixAdapt : public Fix { int anypair, anybond, anyangle; int nlevels_respa; char *id_fix_diam, *id_fix_chg; - class FixStorePeratom *fix_diam, *fix_chg; + class FixStoreAtom *fix_diam, *fix_chg; double previous_diam_scale, previous_chg_scale; int discflag; diff --git a/src/fix_balance.cpp b/src/fix_balance.cpp index 06274eea03..844ffe474e 100644 --- a/src/fix_balance.cpp +++ b/src/fix_balance.cpp @@ -19,7 +19,7 @@ #include "comm.h" #include "domain.h" #include "error.h" -#include "fix_store_peratom.h" +#include "fix_store_atom.h" #include "force.h" #include "irregular.h" #include "kspace.h" diff --git a/src/fix_store_peratom.cpp b/src/fix_store_atom.cpp similarity index 88% rename from src/fix_store_peratom.cpp rename to src/fix_store_atom.cpp index 2a8429e42f..461959d264 100644 --- a/src/fix_store_peratom.cpp +++ b/src/fix_store_atom.cpp @@ -11,7 +11,7 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#include "fix_store_peratom.h" +#include "fix_store_atom.h" #include "atom.h" #include "error.h" @@ -32,10 +32,10 @@ using namespace FixConst; /* ---------------------------------------------------------------------- */ -FixStorePeratom::FixStorePeratom(LAMMPS *lmp, int narg, char **arg) : +FixStoreAtom::FixStoreAtom(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), vstore(nullptr), astore(nullptr) { - if (narg != 7) error->all(FLERR, "Illegal fix STORE/PERATOM command: number of args"); + if (narg != 7) error->all(FLERR, "Illegal fix STORE/ATOM command: number of args"); disable = 0; @@ -52,7 +52,7 @@ FixStorePeratom::FixStorePeratom(LAMMPS *lmp, int narg, char **arg) : else if (n1 > 0 && n2 > 0) tensorflag = 1; else - error->all(FLERR, "Illegal fix STORE/PERATOM dimension args: {} {}", n1, n2); + error->all(FLERR, "Illegal fix STORE/ATOM dimension args: {} {}", n1, n2); if (vecflag || arrayflag) nvalues = n1; @@ -71,7 +71,7 @@ FixStorePeratom::FixStorePeratom(LAMMPS *lmp, int narg, char **arg) : astore = nullptr; tstore = nullptr; - FixStorePeratom::grow_arrays(atom->nmax); + FixStoreAtom::grow_arrays(atom->nmax); atom->add_callback(Atom::GROW); if (restartflag) atom->add_callback(Atom::RESTART); if (ghostflag) atom->add_callback(Atom::BORDER); @@ -93,7 +93,7 @@ FixStorePeratom::FixStorePeratom(LAMMPS *lmp, int narg, char **arg) : /* ---------------------------------------------------------------------- */ -FixStorePeratom::~FixStorePeratom() +FixStoreAtom::~FixStoreAtom() { // unregister callbacks to this fix from Atom class @@ -108,7 +108,7 @@ FixStorePeratom::~FixStorePeratom() /* ---------------------------------------------------------------------- */ -int FixStorePeratom::setmask() +int FixStoreAtom::setmask() { int mask = 0; return mask; @@ -118,7 +118,7 @@ int FixStorePeratom::setmask() allocate atom-based array ------------------------------------------------------------------------- */ -void FixStorePeratom::grow_arrays(int nmax) +void FixStoreAtom::grow_arrays(int nmax) { if (vecflag) memory->grow(vstore, nmax, "store:vstore"); @@ -132,7 +132,7 @@ void FixStorePeratom::grow_arrays(int nmax) copy values within local atom-based array ------------------------------------------------------------------------- */ -void FixStorePeratom::copy_arrays(int i, int j, int /*delflag*/) +void FixStoreAtom::copy_arrays(int i, int j, int /*delflag*/) { if (disable) return; @@ -149,7 +149,7 @@ void FixStorePeratom::copy_arrays(int i, int j, int /*delflag*/) pack values for border communication at re-neighboring ------------------------------------------------------------------------- */ -int FixStorePeratom::pack_border(int n, int *list, double *buf) +int FixStoreAtom::pack_border(int n, int *list, double *buf) { int i, j, k; @@ -179,7 +179,7 @@ int FixStorePeratom::pack_border(int n, int *list, double *buf) unpack values for border communication at re-neighboring ------------------------------------------------------------------------- */ -int FixStorePeratom::unpack_border(int n, int first, double *buf) +int FixStoreAtom::unpack_border(int n, int first, double *buf) { int i, k, last; @@ -203,7 +203,7 @@ int FixStorePeratom::unpack_border(int n, int first, double *buf) pack values in local atom-based array for exchange with another proc ------------------------------------------------------------------------- */ -int FixStorePeratom::pack_exchange(int i, double *buf) +int FixStoreAtom::pack_exchange(int i, double *buf) { if (disable) return 0; @@ -222,7 +222,7 @@ int FixStorePeratom::pack_exchange(int i, double *buf) unpack values in local atom-based array from exchange with another proc ------------------------------------------------------------------------- */ -int FixStorePeratom::unpack_exchange(int nlocal, double *buf) +int FixStoreAtom::unpack_exchange(int nlocal, double *buf) { if (disable) return 0; @@ -241,7 +241,7 @@ int FixStorePeratom::unpack_exchange(int nlocal, double *buf) pack values in local atom-based arrays for restart file ------------------------------------------------------------------------- */ -int FixStorePeratom::pack_restart(int i, double *buf) +int FixStoreAtom::pack_restart(int i, double *buf) { if (disable) { buf[0] = 0; @@ -266,7 +266,7 @@ int FixStorePeratom::pack_restart(int i, double *buf) unpack values from atom->extra array to restart the fix ------------------------------------------------------------------------- */ -void FixStorePeratom::unpack_restart(int nlocal, int nth) +void FixStoreAtom::unpack_restart(int nlocal, int nth) { if (disable) return; @@ -292,7 +292,7 @@ void FixStorePeratom::unpack_restart(int nlocal, int nth) maxsize of any atom's restart data ------------------------------------------------------------------------- */ -int FixStorePeratom::maxsize_restart() +int FixStoreAtom::maxsize_restart() { if (disable) return 1; return nvalues + 1; @@ -302,17 +302,17 @@ int FixStorePeratom::maxsize_restart() size of atom nlocal's restart data ------------------------------------------------------------------------- */ -int FixStorePeratom::size_restart(int /*nlocal*/) +int FixStoreAtom::size_restart(int /*nlocal*/) { if (disable) return 1; return nvalues + 1; } /* ---------------------------------------------------------------------- - memory usage of global or peratom atom-based array + memory usage of per-atom atom-based array ------------------------------------------------------------------------- */ -double FixStorePeratom::memory_usage() +double FixStoreAtom::memory_usage() { return (double) atom->nmax * nvalues * sizeof(double); } diff --git a/src/fix_store_peratom.h b/src/fix_store_atom.h similarity index 89% rename from src/fix_store_peratom.h rename to src/fix_store_atom.h index 9434ffe799..e13b8aa37d 100644 --- a/src/fix_store_peratom.h +++ b/src/fix_store_atom.h @@ -13,26 +13,26 @@ #ifdef FIX_CLASS // clang-format off -FixStyle(STORE/PERATOM,FixStorePeratom); +FixStyle(STORE/ATOM,FixStoreAtom); // clang-format on #else -#ifndef LMP_FIX_STORE_PERATOM_H -#define LMP_FIX_STORE_PERATOM_H +#ifndef LMP_FIX_STORE_ATOM_H +#define LMP_FIX_STORE_ATOM_H #include "fix.h" namespace LAMMPS_NS { -class FixStorePeratom : public Fix { +class FixStoreAtom : public Fix { public: double *vstore; // vector storage double **astore; // array storage double ***tstore; // tensor (3d array) storage int disable; // 1 if operations (except grow) are currently disabled - FixStorePeratom(class LAMMPS *, int, char **); - ~FixStorePeratom() override; + FixStoreAtom(class LAMMPS *, int, char **); + ~FixStoreAtom() override; int setmask() override; void grow_arrays(int) override; diff --git a/src/fix_store_global.cpp b/src/fix_store_global.cpp index cebf4f7690..028d35810e 100644 --- a/src/fix_store_global.cpp +++ b/src/fix_store_global.cpp @@ -182,7 +182,7 @@ void FixStoreGlobal::restart(char *buf) } /* ---------------------------------------------------------------------- - memory usage of global or peratom atom-based array + memory usage of global data ------------------------------------------------------------------------- */ double FixStoreGlobal::memory_usage() diff --git a/src/variable.cpp b/src/variable.cpp index 56a0604a7e..19aed73734 100644 --- a/src/variable.cpp +++ b/src/variable.cpp @@ -20,7 +20,7 @@ #include "domain.h" #include "error.h" #include "fix.h" -#include "fix_store_peratom.h" +#include "fix_store_atom.h" #include "group.h" #include "info.h" #include "input.h" @@ -5027,8 +5027,8 @@ VarReader::VarReader(LAMMPS *lmp, char *name, char *file, int flag) : error->all(FLERR,"Cannot use atomfile-style variable unless an atom map exists"); id_fix = utils::strdup(std::string(name) + "_VARIABLE_STORE"); - fixstore = dynamic_cast( - modify->add_fix(std::string(id_fix) + " all STORE/PERATOM 1 0 0 0")); + fixstore = dynamic_cast( + modify->add_fix(std::string(id_fix) + " all STORE/ATOM 1 0 0 0")); buffer = new char[CHUNK*MAXLINE]; } } diff --git a/src/variable.h b/src/variable.h index 6326c5cc78..825a207d78 100644 --- a/src/variable.h +++ b/src/variable.h @@ -152,7 +152,7 @@ class Variable : protected Pointers { class VarReader : protected Pointers { public: - class FixStorePeratom *fixstore; + class FixStoreAtom *fixstore; char *id_fix; VarReader(class LAMMPS *, char *, char *, int);