From f6523851b590019e67b39e812fdb6620ce49682e Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Wed, 20 Apr 2022 08:58:53 -0600 Subject: [PATCH 01/15] Explcitly state that fix nve uses velocity-Verlet integrator --- doc/src/fix_nve.rst | 8 ++++++-- doc/src/run_style.rst | 3 ++- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/doc/src/fix_nve.rst b/doc/src/fix_nve.rst index da2184c99f..09eea1299d 100644 --- a/doc/src/fix_nve.rst +++ b/doc/src/fix_nve.rst @@ -35,9 +35,13 @@ consistent with the microcanonical ensemble (NVE) provided there are (full) periodic boundary conditions and no other "manipulations" of the system (e.g. fixes that modify forces or velocities). +By default, this fix invokes the velocity form of the +Störmer-Verlet time integration algorithm (velocity-Verlet). Other schemes +can be invoked using the :doc:`run_style ` command. + ---------- -.. include:: accel_styles.rst +.. Include:: accel_styles.rst ---------- @@ -57,7 +61,7 @@ Restrictions Related commands """""""""""""""" -:doc:`fix nvt `, :doc:`fix npt ` +:doc:`fix nvt `, :doc:`fix npt `, :doc:`run_style ` Default """"""" diff --git a/doc/src/run_style.rst b/doc/src/run_style.rst index fd63c82b90..8becbec671 100644 --- a/doc/src/run_style.rst +++ b/doc/src/run_style.rst @@ -67,7 +67,8 @@ Description Choose the style of time integrator used for molecular dynamics simulations performed by LAMMPS. -The *verlet* style is a standard velocity-Verlet integrator. +The *verlet* style is the velocity form of the +Störmer-Verlet time integration algorithm (velocity-Verlet) ---------- From 6dcafd693fa56985e161d5af448f05c308864fc6 Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Thu, 21 Apr 2022 10:13:05 -0600 Subject: [PATCH 02/15] Update fix_nve.rst --- doc/src/fix_nve.rst | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/doc/src/fix_nve.rst b/doc/src/fix_nve.rst index 09eea1299d..a0e186df29 100644 --- a/doc/src/fix_nve.rst +++ b/doc/src/fix_nve.rst @@ -35,9 +35,9 @@ consistent with the microcanonical ensemble (NVE) provided there are (full) periodic boundary conditions and no other "manipulations" of the system (e.g. fixes that modify forces or velocities). -By default, this fix invokes the velocity form of the -Störmer-Verlet time integration algorithm (velocity-Verlet). Other schemes -can be invoked using the :doc:`run_style ` command. +This fix invokes the velocity form of the +Störmer-Verlet time integration algorithm (velocity-Verlet). Other +time integration options can be invoked using the :doc:`run_style ` command. ---------- From 8cb47c850439102993fdad05f094d6ce0e1f0e58 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Thu, 21 Apr 2022 15:02:47 -0600 Subject: [PATCH 03/15] change update of dynamic groups to post_force location in timestep --- doc/src/group.rst | 16 +++++--- src/compute_cluster_atom.cpp | 61 ++++------------------------- src/compute_coord_atom.cpp | 22 ++++++----- src/fix_group.cpp | 76 ++++++++++++++++++++---------------- src/fix_group.h | 6 ++- src/modify.cpp | 38 ++++++++++++++++++ src/modify.h | 5 +++ src/update.cpp | 1 - src/update.h | 1 - 9 files changed, 119 insertions(+), 107 deletions(-) diff --git a/doc/src/group.rst b/doc/src/group.rst index 1720ecfe1a..9d93e949f9 100644 --- a/doc/src/group.rst +++ b/doc/src/group.rst @@ -258,11 +258,17 @@ assignment is made at the beginning of the minimization, but not during the iterations of the minimizer. The point in the timestep at which atoms are assigned to a dynamic -group is after the initial stage of velocity Verlet time integration -has been performed, and before neighbor lists or forces are computed. -This is the point in the timestep where atom positions have just -changed due to the time integration, so the region criterion should be -accurate, if applied. +group is after interatomic forces have been computed, but before any +fixes which alter forces or otherwise update the system have been +invoked. This means that atom positions have been updated, neighbor +lists and ghost atoms are current, and both intermolecular and +intramolecular forces have been calculated based on the new +coordinates. Thus the region criterion, if applied, should be +accurate. Also, any computes invoked by an atom-style variable should +use updated information for that timestep, e.g. potential energy/atom +or coordination number/atom. Similarly, fixes or computes which are +invoked after that point in the timestep, should operate on the new +group of atoms. .. note:: diff --git a/src/compute_cluster_atom.cpp b/src/compute_cluster_atom.cpp index 07a0c700ee..bb4a7122e7 100644 --- a/src/compute_cluster_atom.cpp +++ b/src/compute_cluster_atom.cpp @@ -30,8 +30,6 @@ using namespace LAMMPS_NS; -enum { CLUSTER, MASK, COORDS }; - /* ---------------------------------------------------------------------- */ ComputeClusterAtom::ComputeClusterAtom(LAMMPS *lmp, int narg, char **arg) : @@ -44,7 +42,7 @@ ComputeClusterAtom::ComputeClusterAtom(LAMMPS *lmp, int narg, char **arg) : peratom_flag = 1; size_peratom_cols = 0; - comm_forward = 3; + comm_forward = 1; nmax = 0; } @@ -117,22 +115,6 @@ void ComputeClusterAtom::compute_peratom() numneigh = list->numneigh; firstneigh = list->firstneigh; - // if update->post_integrate set: - // a dynamic group in FixGroup is invoking a variable with this compute - // thus ghost atom coords need to be up-to-date after initial_integrate() - - if (update->post_integrate) { - commflag = COORDS; - comm->forward_comm(this); - } - - // if group is dynamic, insure ghost atom masks are current - - if (group->dynamic[igroup]) { - commflag = MASK; - comm->forward_comm(this); - } - // every atom starts in its own cluster, with clusterID = atomID tagint *tag = atom->tag; @@ -153,7 +135,6 @@ void ComputeClusterAtom::compute_peratom() // iterate until no changes in my atoms // then check if any proc made changes - commflag = CLUSTER; double **x = atom->x; int change, done, anychange; @@ -203,31 +184,15 @@ void ComputeClusterAtom::compute_peratom() /* ---------------------------------------------------------------------- */ -int ComputeClusterAtom::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, - int * /*pbc*/) +int ComputeClusterAtom::pack_forward_comm(int n, int *list, double *buf, + int /*pbc_flag*/, int * /*pbc*/) { int i, j, m; m = 0; - if (commflag == CLUSTER) { - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = clusterID[j]; - } - } else if (commflag == MASK) { - int *mask = atom->mask; - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = ubuf(mask[j]).d; - } - } else if (commflag == COORDS) { - double **x = atom->x; - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0]; - buf[m++] = x[j][1]; - buf[m++] = x[j][2]; - } + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = clusterID[j]; } return m; @@ -241,19 +206,7 @@ void ComputeClusterAtom::unpack_forward_comm(int n, int first, double *buf) m = 0; last = first + n; - if (commflag == CLUSTER) { - for (i = first; i < last; i++) clusterID[i] = buf[m++]; - } else if (commflag == MASK) { - int *mask = atom->mask; - for (i = first; i < last; i++) mask[i] = (int) ubuf(buf[m++]).i; - } else if (commflag == COORDS) { - double **x = atom->x; - for (i = first; i < last; i++) { - x[i][0] = buf[m++]; - x[i][1] = buf[m++]; - x[i][2] = buf[m++]; - } - } + for (i = first; i < last; i++) clusterID[i] = buf[m++]; } /* ---------------------------------------------------------------------- diff --git a/src/compute_coord_atom.cpp b/src/compute_coord_atom.cpp index e7cd73df19..48646b3f38 100644 --- a/src/compute_coord_atom.cpp +++ b/src/compute_coord_atom.cpp @@ -259,14 +259,16 @@ void ComputeCoordAtom::compute_peratom() j = jlist[jj]; j &= NEIGHMASK; - jtype = type[j]; - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx * delx + dely * dely + delz * delz; - if (rsq < cutsq) { - for (m = 0; m < ncol; m++) - if (jtype >= typelo[m] && jtype <= typehi[m]) count[m] += 1.0; + if (mask[j] & jgroupbit) { + jtype = type[j]; + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx * delx + dely * dely + delz * delz; + if (rsq < cutsq) { + for (m = 0; m < ncol; m++) + if (jtype >= typelo[m] && jtype <= typehi[m]) count[m] += 1.0; + } } } } @@ -309,8 +311,8 @@ void ComputeCoordAtom::compute_peratom() /* ---------------------------------------------------------------------- */ -int ComputeCoordAtom::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, - int * /*pbc*/) +int ComputeCoordAtom::pack_forward_comm(int n, int *list, double *buf, + int /*pbc_flag*/, int * /*pbc*/) { int i, m = 0, j; for (i = 0; i < n; ++i) { diff --git a/src/fix_group.cpp b/src/fix_group.cpp index 1291721f0a..73143b2f4f 100644 --- a/src/fix_group.cpp +++ b/src/fix_group.cpp @@ -44,6 +44,8 @@ idregion(nullptr), idvar(nullptr), idprop(nullptr) gbit = group->bitmask[group->find(dgroupid)]; gbitinverse = group->inversemask[group->find(dgroupid)]; + comm_forward = 1; + // process optional args regionflag = 0; @@ -106,8 +108,6 @@ FixGroup::~FixGroup() int FixGroup::setmask() { int mask = 0; - mask |= POST_INTEGRATE; - mask |= POST_INTEGRATE_RESPA; return mask; } @@ -147,29 +147,6 @@ void FixGroup::init() if (iprop < 0 || cols) error->all(FLERR,"Group dynamic command custom property vector does not exist"); } - - // warn if any FixGroup is not at tail end of all post_integrate fixes - - Fix **fix = modify->fix; - int *fmask = modify->fmask; - int nfix = modify->nfix; - - int n = 0; - for (int i = 0; i < nfix; i++) if (POST_INTEGRATE & fmask[i]) n++; - int warn = 0; - for (int i = 0; i < nfix; i++) { - if (POST_INTEGRATE & fmask[i]) { - for (int j = i+1; j < nfix; j++) { - if (POST_INTEGRATE & fmask[j]) { - if (strstr(fix[j]->id,"GROUP_") != fix[j]->id) warn = 1; - } - } - } - } - - if (warn && comm->me == 0) - error->warning(FLERR,"One or more dynamic groups may not be " - "updated at correct point in timestep"); } /* ---------------------------------------------------------------------- @@ -183,7 +160,7 @@ void FixGroup::setup(int /*vflag*/) /* ---------------------------------------------------------------------- */ -void FixGroup::post_integrate() +void FixGroup::post_force(int /*vflag*/) { // only assign atoms to group on steps that are multiples of nevery @@ -192,9 +169,9 @@ void FixGroup::post_integrate() /* ---------------------------------------------------------------------- */ -void FixGroup::post_integrate_respa(int ilevel, int /*iloop*/) +void FixGroup::post_force_respa(int vflag, int ilevel, int /*iloop*/) { - if (ilevel == nlevels_respa-1) post_integrate(); + if (ilevel == nlevels_respa-1) post_force(vflag); } /* ---------------------------------------------------------------------- */ @@ -204,7 +181,6 @@ void FixGroup::set_group() int nlocal = atom->nlocal; // invoke atom-style variable if defined - // set post_integrate flag to 1, then unset after // this is for any compute to check if it needs to // operate differently due to invocation this early in timestep // e.g. perform ghost comm update due to atoms having just moved @@ -214,12 +190,10 @@ void FixGroup::set_group() double *dvector = nullptr; if (varflag) { - update->post_integrate = 1; modify->clearstep_compute(); - memory->create(var,nlocal,"fix/group:varvalue"); + memory->create(var,nlocal,"fix/group:var"); input->variable->compute_atom(ivar,igroup,var,1,0); modify->addstep_compute(update->ntimestep + nevery); - update->post_integrate = 0; } // set ptr to custom atom vector @@ -233,8 +207,6 @@ void FixGroup::set_group() // set mask for each atom // only in group if in parent group, in region, variable is non-zero - // if compute, fix, etc needs updated masks of ghost atoms, - // it must do forward_comm() to update them double **x = atom->x; int *mask = atom->mask; @@ -256,6 +228,42 @@ void FixGroup::set_group() } if (varflag) memory->destroy(var); + + // insure ghost atom masks are also updated + + comm->forward_comm(this); +} + +/* ---------------------------------------------------------------------- */ + +int FixGroup::pack_forward_comm(int n, int *list, double *buf, + int /*pbc_flag*/, int * /*pbc*/) +{ + int i, j, m; + + int *mask = atom->mask; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = ubuf(mask[j]).d; + } + + return m; +} + +/* ---------------------------------------------------------------------- */ + +void FixGroup::unpack_forward_comm(int n, int first, double *buf) +{ + int i, m, last; + + m = 0; + last = first + n; + + int *mask = atom->mask; + + for (i = first; i < last; i++) mask[i] = (int) ubuf(buf[m++]).i; } /* ---------------------------------------------------------------------- */ diff --git a/src/fix_group.h b/src/fix_group.h index ffa2afcfc3..e4c4c9bd32 100644 --- a/src/fix_group.h +++ b/src/fix_group.h @@ -31,8 +31,10 @@ class FixGroup : public Fix { int setmask() override; void init() override; void setup(int) override; - void post_integrate() override; - void post_integrate_respa(int, int) override; + void post_force(int) override; + void post_force_respa(int, int, int) override; + int pack_forward_comm(int, int *, double *, int, int *) override; + void unpack_forward_comm(int, int, double *) override; void *extract(const char *, int &) override; private: diff --git a/src/modify.cpp b/src/modify.cpp index 7554079e2a..ec74b03554 100644 --- a/src/modify.cpp +++ b/src/modify.cpp @@ -58,6 +58,8 @@ Modify::Modify(LAMMPS *lmp) : Pointers(lmp) n_pre_force_respa = n_post_force_respa = n_final_integrate_respa = 0; n_min_pre_exchange = n_min_pre_force = n_min_pre_reverse = 0; n_min_post_force = n_min_energy = 0; + + n_post_force_group = 0; n_timeflag = -1; fix = nullptr; @@ -76,6 +78,8 @@ Modify::Modify(LAMMPS *lmp) : Pointers(lmp) end_of_step_every = nullptr; + list_post_force_group = nullptr; + list_timeflag = nullptr; nfix_restart_global = 0; @@ -156,6 +160,7 @@ Modify::~Modify() delete[] list_min_pre_reverse; delete[] list_min_post_force; delete[] list_min_energy; + delete[] list_post_force_group; delete[] end_of_step_every; delete[] list_timeflag; @@ -226,6 +231,7 @@ void Modify::init() list_init_energy_couple(n_energy_couple, list_energy_couple); list_init_energy_global(n_energy_global, list_energy_global); list_init_energy_atom(n_energy_atom, list_energy_atom); + list_init_post_force_group(n_post_force_group, list_post_force_group); list_init(INITIAL_INTEGRATE_RESPA, n_initial_integrate_respa, list_initial_integrate_respa); list_init(POST_INTEGRATE_RESPA, n_post_integrate_respa, list_post_integrate_respa); @@ -441,10 +447,16 @@ void Modify::pre_reverse(int eflag, int vflag) /* ---------------------------------------------------------------------- post_force call, only for relevant fixes + first call any instances of fix GROUP if they exist ------------------------------------------------------------------------- */ void Modify::post_force(int vflag) { + if (n_post_force_group) { + for (int i = 0; i < n_post_force_group; i++) + fix[list_post_force_group[i]]->post_force(vflag); + } + for (int i = 0; i < n_post_force; i++) fix[list_post_force[i]]->post_force(vflag); } @@ -585,10 +597,16 @@ void Modify::pre_force_respa(int vflag, int ilevel, int iloop) /* ---------------------------------------------------------------------- rRESPA post_force call, only for relevant fixes + first call any instances of fix GROUP if they exist ------------------------------------------------------------------------- */ void Modify::post_force_respa(int vflag, int ilevel, int iloop) { + if (n_post_force_group) { + for (int i = 0; i < n_post_force_group; i++) + fix[list_post_force_group[i]]->post_force_respa(vflag, ilevel, iloop); + } + for (int i = 0; i < n_post_force_respa; i++) fix[list_post_force_respa[i]]->post_force_respa(vflag, ilevel, iloop); } @@ -1716,6 +1734,26 @@ void Modify::list_init_energy_atom(int &n, int *&list) if (fix[i]->energy_peratom_flag && fix[i]->thermo_energy) list[n++] = i; } +/* ---------------------------------------------------------------------- + create list of fix indices for fix GROUP + are invoked first in post_force() or post_force_respa() +------------------------------------------------------------------------- */ + +void Modify::list_init_post_force_group(int &n, int *&list) +{ + delete[] list; + + n = 0; + for (int i = 0; i < nfix; i++) + if (strcmp(fix[i]->style,"GROUP") == 0) n++; + list = new int[n]; + + n = 0; + for (int i = 0; i < nfix; i++) + if (strcmp(fix[i]->style,"GROUP") == 0) + list[n++] = i; +} + /* ---------------------------------------------------------------------- create list of compute indices for computes which store invocation times ------------------------------------------------------------------------- */ diff --git a/src/modify.h b/src/modify.h index df8e52752f..da5d00a799 100644 --- a/src/modify.h +++ b/src/modify.h @@ -163,6 +163,9 @@ class Modify : protected Pointers { int *end_of_step_every; + int n_post_force_group; // list of fix GROUPs for post_force invocation + int *list_post_force_group; + int n_timeflag; // list of computes that store time invocation int *list_timeflag; @@ -187,6 +190,8 @@ class Modify : protected Pointers { void list_init_energy_couple(int &, int *&); void list_init_energy_global(int &, int *&); void list_init_energy_atom(int &, int *&); + void list_init_post_force_group(int &, int *&); + void list_init_post_force_respa_group(int &, int *&); void list_init_dofflag(int &, int *&); void list_init_compute(); diff --git a/src/update.cpp b/src/update.cpp index 3870f6d38e..297053e39d 100644 --- a/src/update.cpp +++ b/src/update.cpp @@ -60,7 +60,6 @@ Update::Update(LAMMPS *lmp) : Pointers(lmp) beginstep = endstep = 0; restrict_output = 0; setupflag = 0; - post_integrate = 0; multireplica = 0; eflag_global = vflag_global = -1; diff --git a/src/update.h b/src/update.h index afcafb87cb..c3e79b72a2 100644 --- a/src/update.h +++ b/src/update.h @@ -35,7 +35,6 @@ class Update : protected Pointers { int max_eval; // max force evaluations for minimizer int restrict_output; // 1 if output should not write dump/restart int setupflag; // set when setup() is computing forces - int post_integrate; // 1 if now at post_integrate() in timestep int multireplica; // 1 if min across replicas, else 0 int dt_default; // 1 if dt is at default value, else 0 From 0ece11c491043c79aee34804dc0f02489c8e10d5 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Thu, 21 Apr 2022 15:19:04 -0600 Subject: [PATCH 04/15] add a code comment --- src/fix_group.cpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/fix_group.cpp b/src/fix_group.cpp index 73143b2f4f..112776251a 100644 --- a/src/fix_group.cpp +++ b/src/fix_group.cpp @@ -181,9 +181,11 @@ void FixGroup::set_group() int nlocal = atom->nlocal; // invoke atom-style variable if defined - // this is for any compute to check if it needs to - // operate differently due to invocation this early in timestep - // e.g. perform ghost comm update due to atoms having just moved + // NOTE: after variable invocation could reset invoked computes to not-invoked + // this would avoid an issue where other post-force fixes + // change the compute result since it will not be re-invoked at end-of-step, + // e.g. if compute pe/atom includes pe contributions from fixes + double *var = nullptr; int *ivector = nullptr; From 3f3c481554ee190fe352a4566bb228445b55ed7f Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Thu, 21 Apr 2022 17:00:11 -0600 Subject: [PATCH 05/15] add support to fix adapt for angle coeffs --- doc/src/fix_adapt.rst | 81 +++++++++++++++------ src/EXTRA-MOLECULE/bond_fene_nm.cpp | 2 +- src/EXTRA-MOLECULE/bond_gaussian.cpp | 4 +- src/MOLECULE/angle_harmonic.cpp | 12 ++++ src/MOLECULE/angle_harmonic.h | 1 + src/MOLECULE/bond_fene.cpp | 2 +- src/MOLECULE/bond_gromos.cpp | 7 +- src/MOLECULE/bond_harmonic.cpp | 10 ++- src/angle.cpp | 12 ++++ src/angle.h | 5 ++ src/bond.cpp | 7 +- src/bond.h | 6 +- src/fix_adapt.cpp | 102 ++++++++++++++++++++++++--- src/fix_adapt.h | 8 ++- 14 files changed, 202 insertions(+), 57 deletions(-) diff --git a/doc/src/fix_adapt.rst b/doc/src/fix_adapt.rst index b7151334db..0f68219358 100644 --- a/doc/src/fix_adapt.rst +++ b/doc/src/fix_adapt.rst @@ -14,7 +14,7 @@ Syntax * adapt = style name of this fix command * N = adapt simulation settings every this many timesteps * one or more attribute/arg pairs may be appended -* attribute = *pair* or *bond* or *kspace* or *atom* +* attribute = *pair* or *bond* or *angle* or *kspace* or *atom* .. parsed-literal:: @@ -28,11 +28,16 @@ Syntax bparam = parameter to adapt over time I = type bond to set parameter for v_name = variable with name that calculates value of bparam + *angle* args = astyle aparam I v_name + astyle = angle style name, e.g. harmonic + aparam = parameter to adapt over time + I = type angle to set parameter for + v_name = variable with name that calculates value of aparam *kspace* arg = v_name v_name = variable with name that calculates scale factor on K-space terms - *atom* args = aparam v_name - aparam = parameter to adapt over time - v_name = variable with name that calculates value of aparam + *atom* args = atomparam v_name + atomparam = parameter to adapt over time + v_name = variable with name that calculates value of atomparam * zero or more keyword/value pairs may be appended * keyword = *scale* or *reset* or *mass* @@ -283,30 +288,60 @@ operates. The only difference is that now a bond coefficient for a given bond type is adapted. A wild-card asterisk can be used in place of or in conjunction with -the bond type argument to set the coefficients for multiple bond types. -This takes the form "\*" or "\*n" or "n\*" or "m\*n". If N = the number of -atom types, then an asterisk with no numeric values means all types -from 1 to N. A leading asterisk means all types from 1 to n (inclusive). -A trailing asterisk means all types from n to N (inclusive). A middle -asterisk means all types from m to n (inclusive). +the bond type argument to set the coefficients for multiple bond +types. This takes the form "\*" or "\*n" or "n\*" or "m\*n". If N = +the number of bond types, then an asterisk with no numeric values +means all types from 1 to N. A leading asterisk means all types from +1 to n (inclusive). A trailing asterisk means all types from n to N +(inclusive). A middle asterisk means all types from m to n +(inclusive). Currently *bond* does not support bond_style hybrid nor bond_style hybrid/overlay as bond styles. The only bonds that currently are working with fix_adapt are -+------------------------------------+-------+------------+ -| :doc:`class2 ` | r0 | type bonds | -+------------------------------------+-------+------------+ -| :doc:`fene ` | k, r0 | type bonds | -+------------------------------------+-------+------------+ -| :doc:`gromos ` | k, r0 | type bonds | -+------------------------------------+-------+------------+ -| :doc:`harmonic ` | k,r0 | type bonds | -+------------------------------------+-------+------------+ -| :doc:`morse ` | r0 | type bonds | -+------------------------------------+-------+------------+ -| :doc:`nonlinear ` | r0 | type bonds | -+------------------------------------+-------+------------+ ++------------------------------------+-------+-----------------+ +| :doc:`class2 ` | r0 | type bonds | ++------------------------------------+-------+-----------------+ +| :doc:`fene ` | k, r0 | type bonds | ++------------------------------------+-------+-----------------+ +| :doc:`fene/nm ` | k, r0 | type bonds | ++------------------------------------+-------+-----------------+ +| :doc:`gromos ` | k, r0 | type bonds | ++------------------------------------+-------+-----------------+ +| :doc:`harmonic ` | k,r0 | type bonds | ++------------------------------------+-------+-----------------+ +| :doc:`morse ` | r0 | type bonds | ++------------------------------------+-------+-----------------+ +| :doc:`nonlinear ` | epsilon,r0 | type bonds | ++------------------------------------+-------+-----------------+ + +---------- + +The *angle* keyword uses the specified variable to change the value of +an angle coefficient over time, very similar to how the *pair* keyword +operates. The only difference is that now an angle coefficient for a +given angle type is adapted. + +A wild-card asterisk can be used in place of or in conjunction with +the angle type argument to set the coefficients for multiple angle +types. This takes the form "\*" or "\*n" or "n\*" or "m\*n". If N = +the number of angle types, then an asterisk with no numeric values +means all types from 1 to N. A leading asterisk means all types from +1 to n (inclusive). A trailing asterisk means all types from n to N +(inclusive). A middle asterisk means all types from m to n +(inclusive). + +Currently *angle* does not support angle_style hybrid nor angle_style +hybrid/overlay as angle styles. The only angles that currently are +working with fix_adapt are + ++------------------------------------+-------+-----------------+ +| :doc:`harmonic ` | k,theta0 | type angles | ++------------------------------------+-------+-----------------+ + +Note that internally, theta0 is stored in radians, so the variable +this fix uses to reset theta0 needs to generate values in radians. ---------- diff --git a/src/EXTRA-MOLECULE/bond_fene_nm.cpp b/src/EXTRA-MOLECULE/bond_fene_nm.cpp index 839217df78..bd9d16b95d 100644 --- a/src/EXTRA-MOLECULE/bond_fene_nm.cpp +++ b/src/EXTRA-MOLECULE/bond_fene_nm.cpp @@ -273,7 +273,7 @@ double BondFENENM::single(int type, double rsq, int /*i*/, int /*j*/, double &ff void *BondFENENM::extract(const char *str, int &dim) { dim = 1; - if (strcmp(str, "kappa") == 0) return (void *) k; + if (strcmp(str, "k") == 0) return (void *) k; if (strcmp(str, "r0") == 0) return (void *) r0; return nullptr; } diff --git a/src/EXTRA-MOLECULE/bond_gaussian.cpp b/src/EXTRA-MOLECULE/bond_gaussian.cpp index c2ab00dfde..0ad125aeca 100644 --- a/src/EXTRA-MOLECULE/bond_gaussian.cpp +++ b/src/EXTRA-MOLECULE/bond_gaussian.cpp @@ -34,9 +34,7 @@ using namespace MathConst; BondGaussian::BondGaussian(LAMMPS *lmp) : Bond(lmp), nterms(nullptr), bond_temperature(nullptr), alpha(nullptr), width(nullptr), r0(nullptr) -{ - reinitflag = 1; -} +{} /* ---------------------------------------------------------------------- */ diff --git a/src/MOLECULE/angle_harmonic.cpp b/src/MOLECULE/angle_harmonic.cpp index 261770d9a3..aa24fa27b4 100644 --- a/src/MOLECULE/angle_harmonic.cpp +++ b/src/MOLECULE/angle_harmonic.cpp @@ -264,3 +264,15 @@ double AngleHarmonic::single(int type, int i1, int i2, int i3) double tk = k[type] * dtheta; return tk * dtheta; } + +/* ---------------------------------------------------------------------- + return ptr to internal members upon request +------------------------------------------------------------------------ */ + +void *AngleHarmonic::extract(const char *str, int &dim) +{ + dim = 1; + if (strcmp(str, "k") == 0) return (void *) k; + if (strcmp(str, "theta0") == 0) return (void *) theta0; + return nullptr; +} diff --git a/src/MOLECULE/angle_harmonic.h b/src/MOLECULE/angle_harmonic.h index 718ac4bd0a..2643ea3a17 100644 --- a/src/MOLECULE/angle_harmonic.h +++ b/src/MOLECULE/angle_harmonic.h @@ -35,6 +35,7 @@ class AngleHarmonic : public Angle { void read_restart(FILE *) override; void write_data(FILE *) override; double single(int, int, int, int) override; + void *extract(const char *, int &) override; protected: double *k, *theta0; diff --git a/src/MOLECULE/bond_fene.cpp b/src/MOLECULE/bond_fene.cpp index 8daf6e862c..2950d48ca6 100644 --- a/src/MOLECULE/bond_fene.cpp +++ b/src/MOLECULE/bond_fene.cpp @@ -265,7 +265,7 @@ double BondFENE::single(int type, double rsq, int /*i*/, int /*j*/, double &ffor void *BondFENE::extract(const char *str, int &dim) { dim = 1; - if (strcmp(str, "kappa") == 0) return (void *) k; + if (strcmp(str, "k") == 0) return (void *) k; if (strcmp(str, "r0") == 0) return (void *) r0; return nullptr; } diff --git a/src/MOLECULE/bond_gromos.cpp b/src/MOLECULE/bond_gromos.cpp index adf3f91044..4572d2c900 100644 --- a/src/MOLECULE/bond_gromos.cpp +++ b/src/MOLECULE/bond_gromos.cpp @@ -30,10 +30,7 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -BondGromos::BondGromos(LAMMPS *_lmp) : Bond(_lmp) -{ - reinitflag = 1; -} +BondGromos::BondGromos(LAMMPS *_lmp) : Bond(_lmp) {} /* ---------------------------------------------------------------------- */ @@ -200,7 +197,7 @@ double BondGromos::single(int type, double rsq, int /*i*/, int /*j*/, double &ff void *BondGromos::extract(const char *str, int &dim) { dim = 1; - if (strcmp(str, "kappa") == 0) return (void *) k; + if (strcmp(str, "k") == 0) return (void *) k; if (strcmp(str, "r0") == 0) return (void *) r0; return nullptr; } diff --git a/src/MOLECULE/bond_harmonic.cpp b/src/MOLECULE/bond_harmonic.cpp index cff766aa3b..687a871f17 100644 --- a/src/MOLECULE/bond_harmonic.cpp +++ b/src/MOLECULE/bond_harmonic.cpp @@ -27,10 +27,7 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -BondHarmonic::BondHarmonic(LAMMPS *_lmp) : Bond(_lmp) -{ - reinitflag = 1; -} +BondHarmonic::BondHarmonic(LAMMPS *_lmp) : Bond(_lmp) {} /* ---------------------------------------------------------------------- */ @@ -201,12 +198,13 @@ double BondHarmonic::single(int type, double rsq, int /*i*/, int /*j*/, double & } /* ---------------------------------------------------------------------- - Return ptr to internal members upon request. + return ptr to internal members upon request ------------------------------------------------------------------------ */ + void *BondHarmonic::extract(const char *str, int &dim) { dim = 1; - if (strcmp(str, "kappa") == 0) return (void *) k; + if (strcmp(str, "k") == 0) return (void *) k; if (strcmp(str, "r0") == 0) return (void *) r0; return nullptr; } diff --git a/src/angle.cpp b/src/angle.cpp index 52d92b72b2..caa86d691e 100644 --- a/src/angle.cpp +++ b/src/angle.cpp @@ -353,3 +353,15 @@ double Angle::memory_usage() bytes += (double) comm->nthreads * maxcvatom * 9 * sizeof(double); return bytes; } + +/* ----------------------------------------------------------------------- + reset all type-based angle params via init() +-------------------------------------------------------------------------- */ + +void Angle::reinit() +{ + if (!reinitflag) + error->all(FLERR, "Fix adapt interface to this angle style not supported"); + + init(); +} diff --git a/src/angle.h b/src/angle.h index 12443fa4f3..8e200ce37b 100644 --- a/src/angle.h +++ b/src/angle.h @@ -36,6 +36,9 @@ class Angle : protected Pointers { // CENTROID_AVAIL = different and implemented // CENTROID_NOTAVAIL = different, not yet implemented + int reinitflag; // 0 if not compatible with fix adapt + // extract() method may still need to be added + // KOKKOS host/device flag and data masks ExecutionSpace execution_space; @@ -57,6 +60,8 @@ class Angle : protected Pointers { virtual void write_data(FILE *) {} virtual double single(int, int, int, int) = 0; virtual double memory_usage(); + virtual void *extract(const char *, int &) { return nullptr; } + void reinit(); protected: int suffix_flag; // suffix compatibility flag diff --git a/src/bond.cpp b/src/bond.cpp index 6ae78bc429..7f7b6aa3a3 100644 --- a/src/bond.cpp +++ b/src/bond.cpp @@ -43,6 +43,7 @@ Bond::Bond(LAMMPS *_lmp) : Pointers(_lmp) energy = 0.0; virial[0] = virial[1] = virial[2] = virial[3] = virial[4] = virial[5] = 0.0; writedata = 1; + reinitflag = 1; comm_forward = comm_reverse = comm_reverse_off = 0; @@ -336,11 +337,13 @@ double Bond::memory_usage() } /* ----------------------------------------------------------------------- - Reset all type-based bond params via init. + reset all type-based bond params via init() -------------------------------------------------------------------------- */ + void Bond::reinit() { - if (!reinitflag) error->all(FLERR, "Fix adapt interface to this bond style not supported"); + if (!reinitflag) + error->all(FLERR, "Fix adapt interface to this bond style not supported"); init(); } diff --git a/src/bond.h b/src/bond.h index 47c8687f50..7b6a70dc7d 100644 --- a/src/bond.h +++ b/src/bond.h @@ -37,7 +37,8 @@ class Bond : protected Pointers { int comm_reverse; // size of reverse communication (0 if none) int comm_reverse_off; // size of reverse comm even if newton off - int reinitflag; // 1 if compatible with fix adapt and alike + int reinitflag; // 0 if not compatible with fix adapt + // extract() method may still need to be added // KOKKOS host/device flag and data masks @@ -61,7 +62,8 @@ class Bond : protected Pointers { virtual double single(int, double, int, int, double &) = 0; virtual double memory_usage(); virtual void *extract(const char *, int &) { return nullptr; } - virtual void reinit(); + void reinit(); + virtual int pack_forward_comm(int, int *, double *, int, int *) {return 0;} virtual void unpack_forward_comm(int, int, double *) {} virtual int pack_reverse_comm(int, int, double *) {return 0;} diff --git a/src/fix_adapt.cpp b/src/fix_adapt.cpp index 2632ccf597..0d9aa00b7a 100644 --- a/src/fix_adapt.cpp +++ b/src/fix_adapt.cpp @@ -14,6 +14,7 @@ #include "fix_adapt.h" +#include "angle.h" #include "atom.h" #include "bond.h" #include "domain.h" @@ -38,7 +39,7 @@ using namespace LAMMPS_NS; using namespace FixConst; using namespace MathConst; -enum{PAIR,KSPACE,ATOM,BOND}; +enum{PAIR,KSPACE,ATOM,BOND,ANGLE}; enum{DIAMETER,CHARGE}; /* ---------------------------------------------------------------------- */ @@ -75,6 +76,10 @@ nadapt(0), id_fix_diam(nullptr), id_fix_chg(nullptr), adapt(nullptr) if (iarg+5 > narg) error->all(FLERR,"Illegal fix adapt command"); nadapt++; iarg += 5; + } else if (strcmp(arg[iarg],"angle") == 0) { + if (iarg+5 > narg) error->all(FLERR,"Illegal fix adapt command"); + nadapt++; + iarg += 5; } else break; } @@ -119,6 +124,20 @@ nadapt(0), id_fix_diam(nullptr), id_fix_chg(nullptr), adapt(nullptr) nadapt++; iarg += 5; + } else if (strcmp(arg[iarg],"angle") == 0) { + if (iarg+5 > narg) error->all(FLERR, "Illegal fix adapt command"); + adapt[nadapt].which = ANGLE; + adapt[nadapt].angle = nullptr; + adapt[nadapt].astyle = utils::strdup(arg[iarg+1]); + adapt[nadapt].aparam = utils::strdup(arg[iarg+2]); + utils::bounds(FLERR,arg[iarg+3],1,atom->nangletypes, + adapt[nadapt].ilo,adapt[nadapt].ihi,error); + if (utils::strmatch(arg[iarg+4],"^v_")) { + adapt[nadapt].var = utils::strdup(arg[iarg+4]+2); + } else error->all(FLERR,"Illegal fix adapt command"); + nadapt++; + iarg += 5; + } else if (strcmp(arg[iarg],"kspace") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt command"); adapt[nadapt].which = KSPACE; @@ -133,12 +152,12 @@ nadapt(0), id_fix_diam(nullptr), id_fix_chg(nullptr), adapt(nullptr) adapt[nadapt].which = ATOM; if (strcmp(arg[iarg+1],"diameter") == 0 || strcmp(arg[iarg+1],"diameter/disc") == 0) { - adapt[nadapt].aparam = DIAMETER; + adapt[nadapt].atomparam = DIAMETER; diamflag = 1; discflag = 0; if (strcmp(arg[iarg+1],"diameter/disc") == 0) discflag = 1; } else if (strcmp(arg[iarg+1],"charge") == 0) { - adapt[nadapt].aparam = CHARGE; + adapt[nadapt].atomparam = CHARGE; chgflag = 1; } else error->all(FLERR,"Illegal fix adapt command"); if (utils::strmatch(arg[iarg+2],"^v_")) { @@ -191,6 +210,13 @@ nadapt(0), id_fix_diam(nullptr), id_fix_chg(nullptr), adapt(nullptr) for (int m = 0; m < nadapt; ++m) if (adapt[m].which == BOND) memory->create(adapt[m].vector_orig,n+1,"adapt:vector_orig"); + + // allocate angle style arrays: + + n = atom->nbondtypes; + for (int m = 0; m < nadapt; ++m) + if (adapt[m].which == ANGLE) + memory->create(adapt[m].vector_orig,n+1,"adapt:vector_orig"); } /* ---------------------------------------------------------------------- */ @@ -207,6 +233,10 @@ FixAdapt::~FixAdapt() delete [] adapt[m].bstyle; delete [] adapt[m].bparam; memory->destroy(adapt[m].vector_orig); + } else if (adapt[m].which == ANGLE) { + delete [] adapt[m].astyle; + delete [] adapt[m].aparam; + memory->destroy(adapt[m].vector_orig); } } delete [] adapt; @@ -357,6 +387,7 @@ void FixAdapt::init() } delete [] pstyle; + } else if (ad->which == BOND) { ad->bond = nullptr; anybond = 1; @@ -383,13 +414,39 @@ void FixAdapt::init() delete [] bstyle; + } else if (ad->which == ANGLE) { + ad->angle = nullptr; + anyangle = 1; + + char *astyle = utils::strdup(ad->astyle); + if (lmp->suffix_enable) + ad->angle = force->angle_match(fmt::format("{}/{}",astyle,lmp->suffix)); + + if (ad->angle == nullptr) ad->angle = force->angle_match(astyle); + if (ad->angle == nullptr ) + error->all(FLERR,"Fix adapt angle style does not exist"); + + void *ptr = ad->angle->extract(ad->aparam,ad->bdim); + + if (ptr == nullptr) + error->all(FLERR,"Fix adapt angle style param not supported"); + + // for angle styles, use a vector + + if (ad->adim == 1) ad->vector = (double *) ptr; + + if (utils::strmatch(force->angle_style,"^hybrid")) + error->all(FLERR,"Fix adapt does not support angle_style hybrid"); + + delete [] astyle; + } else if (ad->which == KSPACE) { if (force->kspace == nullptr) error->all(FLERR,"Fix adapt kspace style does not exist"); kspace_scale = (double *) force->kspace->extract("scale"); } else if (ad->which == ATOM) { - if (ad->aparam == DIAMETER) { + if (ad->atomparam == DIAMETER) { if (!atom->radius_flag) error->all(FLERR,"Fix adapt requires atom attribute diameter"); if (!atom->rmass_flag) @@ -398,7 +455,7 @@ void FixAdapt::init() error->all(FLERR,"Fix adapt requires 2d simulation"); if (!restart_reset) previous_diam_scale = 1.0; } - if (ad->aparam == CHARGE) { + if (ad->atomparam == CHARGE) { if (!atom->q_flag) error->all(FLERR,"Fix adapt requires atom attribute charge"); if (!restart_reset) previous_chg_scale = 1.0; @@ -408,7 +465,7 @@ void FixAdapt::init() if (restart_reset) restart_reset = 0; - // make copy of original pair/bond array values + // make copy of original pair/bond/angle array values for (int m = 0; m < nadapt; m++) { Adapt *ad = &adapt[m]; @@ -422,6 +479,10 @@ void FixAdapt::init() } else if (ad->which == BOND && ad->bdim == 1) { for (i = ad->ilo; i <= ad->ihi; ++i ) ad->vector_orig[i] = ad->vector[i]; + + } else if (ad->which == ANGLE && ad->adim == 1) { + for (i = ad->ilo; i <= ad->ihi; ++i ) + ad->vector_orig[i] = ad->vector[i]; } } @@ -483,7 +544,7 @@ void FixAdapt::post_run() } /* ---------------------------------------------------------------------- - change pair,kspace,atom parameters based on variable evaluation + change pair,bond,angle,kspace,atom parameters based on variable evaluation ------------------------------------------------------------------------- */ void FixAdapt::change_settings() @@ -527,6 +588,18 @@ void FixAdapt::change_settings() ad->vector[i] = value; } + // set angle type array values: + + } else if (ad->which == ANGLE) { + if (ad->adim == 1) { + if (scaleflag) + for (i = ad->ilo; i <= ad->ihi; ++i ) + ad->vector[i] = value*ad->vector_orig[i]; + else + for (i = ad->ilo; i <= ad->ihi; ++i ) + ad->vector[i] = value; + } + // set kspace scale factor } else if (ad->which == KSPACE) { @@ -540,7 +613,7 @@ void FixAdapt::change_settings() // also reset rmass to new value assuming density remains constant // for scaleflag, previous_diam_scale is the scale factor on previous step - if (ad->aparam == DIAMETER) { + if (ad->atomparam == DIAMETER) { double scale; double *radius = atom->radius; double *rmass = atom->rmass; @@ -567,7 +640,7 @@ void FixAdapt::change_settings() // reset charge to new value, for both owned and ghost atoms // for scaleflag, previous_chg_scale is the scale factor on previous step - } else if (ad->aparam == CHARGE) { + } else if (ad->atomparam == CHARGE) { double scale; double *q = atom->q; int *mask = atom->mask; @@ -591,7 +664,7 @@ void FixAdapt::change_settings() modify->addstep_compute(update->ntimestep + nevery); // re-initialize pair styles if any PAIR settings were changed - // ditto for bond styles if any BOND settings were changed + // ditto for bond/angle styles if any BOND/ANGLE settings were changed // this resets other coeffs that may depend on changed values, // and also offset and tail corrections // we must call force->pair->reinit() instead of the individual @@ -601,6 +674,7 @@ void FixAdapt::change_settings() if (anypair) force->pair->reinit(); if (anybond) force->bond->reinit(); + if (anyangle) force->angle->reinit(); // reset KSpace charges if charges have changed @@ -624,7 +698,13 @@ void FixAdapt::restore_settings() } } else if (ad->which == BOND) { - if (ad->pdim == 1) { + if (ad->bdim == 1) { + for (int i = ad->ilo; i <= ad->ihi; i++) + ad->vector[i] = ad->vector_orig[i]; + } + + } else if (ad->which == ANGLE) { + if (ad->adim == 1) { for (int i = ad->ilo; i <= ad->ihi; i++) ad->vector[i] = ad->vector_orig[i]; } diff --git a/src/fix_adapt.h b/src/fix_adapt.h index 939f46f8a2..121ef06ece 100644 --- a/src/fix_adapt.h +++ b/src/fix_adapt.h @@ -45,7 +45,7 @@ class FixAdapt : public Fix { private: int nadapt, resetflag, scaleflag, massflag; - int anypair, anybond; + int anypair, anybond, anyangle; int nlevels_respa; char *id_fix_diam, *id_fix_chg; class FixStore *fix_diam, *fix_chg; @@ -57,14 +57,16 @@ class FixAdapt : public Fix { char *var; char *pstyle, *pparam; char *bstyle, *bparam; + char *astyle, *aparam; int ilo, ihi, jlo, jhi; - int pdim, bdim; + int pdim, bdim, adim; double *scalar, scalar_orig; double *vector, *vector_orig; double **array, **array_orig; - int aparam; + int atomparam; class Pair *pair; class Bond *bond; + class Angle *angle; }; Adapt *adapt; From 13664a01851627801b2754989a9990391dab249a Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 22 Apr 2022 05:25:50 -0400 Subject: [PATCH 06/15] add utility to print url with pointer to error message --- doc/src/Developer_utils.rst | 3 +++ doc/src/Errors.rst | 1 + doc/src/Errors_details.rst | 27 +++++++++++++++++++++++++++ src/read_data.cpp | 4 ++-- src/utils.cpp | 13 ++++++++++--- src/utils.h | 14 +++++++++++++- 6 files changed, 56 insertions(+), 6 deletions(-) create mode 100644 doc/src/Errors_details.rst diff --git a/doc/src/Developer_utils.rst b/doc/src/Developer_utils.rst index 9c6ef67945..aceecd7a46 100644 --- a/doc/src/Developer_utils.rst +++ b/doc/src/Developer_utils.rst @@ -211,6 +211,9 @@ Convenience functions .. doxygenfunction:: logmesg(LAMMPS *lmp, const std::string &mesg) :project: progguide +.. doxygenfunction:: errorurl + :project: progguide + .. doxygenfunction:: flush_buffers(LAMMPS *lmp) :project: progguide diff --git a/doc/src/Errors.rst b/doc/src/Errors.rst index 5975c22c41..3e8ebc7f8e 100644 --- a/doc/src/Errors.rst +++ b/doc/src/Errors.rst @@ -11,6 +11,7 @@ them. :maxdepth: 1 Errors_common + Errors_details Errors_bugs Errors_debug Errors_messages diff --git a/doc/src/Errors_details.rst b/doc/src/Errors_details.rst new file mode 100644 index 0000000000..70fe0e90ff --- /dev/null +++ b/doc/src/Errors_details.rst @@ -0,0 +1,27 @@ +Detailed discussion of errors and warnings +========================================== + +Many errors or warnings are self-explanatory and thus straightforward to +resolve. However, there are also cases, where there is no single cause +and explanation, where LAMMPS can only detect symptoms of an error but +not the exact cause, or where the explanation needs to be more detailed than +what can be fit into a message printed by the program. The following are +discussions of such cases. + +.. _err0001: + +Unknown identifier in data file +------------------------------- + +This error happens when LAMMPS encounters a line of text in an unexpected format +while reading a data file. This is most commonly cause by inconsistent header and +section data. The header section informs LAMMPS how many entries or lines are expected in the +various sections (like Atoms, Masses, Pair Coeffs, *etc.*\ ) of the data file. +If there is a mismatch, LAMMPS will either keep reading beyond the end of a section +or stop reading before the section has ended. + +Such a mismatch can happen unexpectedly when the first line of the data +is *not* a comment as required by the format. That would result in +LAMMPS expecting, for instance, 0 atoms because the "atoms" header line +is treated as a comment. + diff --git a/src/read_data.cpp b/src/read_data.cpp index 5505b5a5a0..b77394bcb3 100644 --- a/src/read_data.cpp +++ b/src/read_data.cpp @@ -744,9 +744,9 @@ void ReadData::command(int narg, char **arg) break; } if (i == nfix) - error->all(FLERR,"Unknown identifier in data file: {}",keyword); + error->all(FLERR,"Unknown identifier in data file: {}{}", keyword, utils::errorurl(1)); - } else error->all(FLERR,"Unknown identifier in data file: {}",keyword); + } else error->all(FLERR,"Unknown identifier in data file: {}{}", keyword, utils::errorurl(1)); parse_keyword(0); } diff --git a/src/utils.cpp b/src/utils.cpp index a8e1b14104..72944c1838 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -22,8 +22,8 @@ #include "memory.h" #include "modify.h" #include "text_file_reader.h" -#include "update.h" #include "universe.h" +#include "update.h" #include #include @@ -137,11 +137,18 @@ void utils::fmtargs_logmesg(LAMMPS *lmp, fmt::string_view format, fmt::format_ar } } +std::string utils::errorurl(int errorcode) +{ + return fmt::format( + "\nFor more information please go to https://docs.lammps.org/Errors_details.html#err{:04d}", + errorcode); +} + void utils::flush_buffers(LAMMPS *lmp) { if (lmp->screen) fflush(lmp->screen); if (lmp->logfile) fflush(lmp->logfile); - if (lmp->universe->uscreen) fflush(lmp->universe->uscreen); + if (lmp->universe->uscreen) fflush(lmp->universe->uscreen); if (lmp->universe->ulogfile) fflush(lmp->universe->ulogfile); } @@ -805,7 +812,7 @@ std::string utils::star_subst(const std::string &name, bigint step, int pad) auto star = name.find('*'); if (star == std::string::npos) return name; - return fmt::format("{}{:0{}}{}",name.substr(0,star),step,pad,name.substr(star+1)); + return fmt::format("{}{:0{}}{}", name.substr(0, star), step, pad, name.substr(star + 1)); } /* ---------------------------------------------------------------------- diff --git a/src/utils.h b/src/utils.h index 86f31508e7..4dbb2fb8bb 100644 --- a/src/utils.h +++ b/src/utils.h @@ -21,7 +21,7 @@ #include -#include // IWYU pragma: export +#include // IWYU pragma: export namespace LAMMPS_NS { @@ -74,6 +74,18 @@ namespace utils { void logmesg(LAMMPS *lmp, const std::string &mesg); + /*! Return text redirecting the user to a specific paragraph in the manual + * + * The LAMMPS manual contains detailed detailed explanations for errors and + * warnings where a simple error message may not be sufficient. These can + * be reached through URLs with a numeric code. This function creates the + * corresponding text to be included into the error message that redirects + * the user to that URL. + * + * \param errorcode number pointing to a paragraph in the manual */ + + std::string errorurl(int errorcode); + /*! Flush output buffers * * This function calls fflush() on screen and logfile FILE pointers From b4e2e2ec34a359079cae7bfe2d3d429cce7ad679 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Fri, 22 Apr 2022 09:12:12 -0600 Subject: [PATCH 07/15] add support for angle cosine --- doc/src/fix_adapt.rst | 16 +++++++++------- src/MOLECULE/angle_cosine.cpp | 11 +++++++++++ src/MOLECULE/angle_cosine.h | 1 + 3 files changed, 21 insertions(+), 7 deletions(-) diff --git a/doc/src/fix_adapt.rst b/doc/src/fix_adapt.rst index 0f68219358..1276adf444 100644 --- a/doc/src/fix_adapt.rst +++ b/doc/src/fix_adapt.rst @@ -297,17 +297,17 @@ means all types from 1 to N. A leading asterisk means all types from (inclusive). Currently *bond* does not support bond_style hybrid nor bond_style -hybrid/overlay as bond styles. The only bonds that currently are -working with fix_adapt are +hybrid/overlay as bond styles. The bond styles that currently work +with fix_adapt are +------------------------------------+-------+-----------------+ | :doc:`class2 ` | r0 | type bonds | +------------------------------------+-------+-----------------+ -| :doc:`fene ` | k, r0 | type bonds | +| :doc:`fene ` | k,r0 | type bonds | +------------------------------------+-------+-----------------+ -| :doc:`fene/nm ` | k, r0 | type bonds | +| :doc:`fene/nm ` | k,r0 | type bonds | +------------------------------------+-------+-----------------+ -| :doc:`gromos ` | k, r0 | type bonds | +| :doc:`gromos ` | k,r0 | type bonds | +------------------------------------+-------+-----------------+ | :doc:`harmonic ` | k,r0 | type bonds | +------------------------------------+-------+-----------------+ @@ -333,12 +333,14 @@ means all types from 1 to N. A leading asterisk means all types from (inclusive). Currently *angle* does not support angle_style hybrid nor angle_style -hybrid/overlay as angle styles. The only angles that currently are -working with fix_adapt are +hybrid/overlay as angle styles. The angle styles that currently work +with fix_adapt are +------------------------------------+-------+-----------------+ | :doc:`harmonic ` | k,theta0 | type angles | +------------------------------------+-------+-----------------+ +| :doc:`cosine ` | k | type angles | ++------------------------------------+-------+-----------------+ Note that internally, theta0 is stored in radians, so the variable this fix uses to reset theta0 needs to generate values in radians. diff --git a/src/MOLECULE/angle_cosine.cpp b/src/MOLECULE/angle_cosine.cpp index d32dc4559d..260ebc3948 100644 --- a/src/MOLECULE/angle_cosine.cpp +++ b/src/MOLECULE/angle_cosine.cpp @@ -234,3 +234,14 @@ double AngleCosine::single(int type, int i1, int i2, int i3) return k[type] * (1.0 + c); } + +/* ---------------------------------------------------------------------- + return ptr to internal members upon request +------------------------------------------------------------------------ */ + +void *AngleCosine::extract(const char *str, int &dim) +{ + dim = 1; + if (strcmp(str, "k") == 0) return (void *) k; + return nullptr; +} diff --git a/src/MOLECULE/angle_cosine.h b/src/MOLECULE/angle_cosine.h index 19b6222c87..e249e7a44c 100644 --- a/src/MOLECULE/angle_cosine.h +++ b/src/MOLECULE/angle_cosine.h @@ -35,6 +35,7 @@ class AngleCosine : public Angle { void read_restart(FILE *) override; void write_data(FILE *) override; double single(int, int, int, int) override; + void *extract(const char *, int &) override; protected: double *k; From c0d0c84f7d6735263a4560327f96e0e9a61e627b Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 22 Apr 2022 12:08:23 -0400 Subject: [PATCH 08/15] update unit test files to implementation changes --- unittest/force-styles/tests/angle-cosine.yaml | 3 ++- unittest/force-styles/tests/angle-harmonic.yaml | 4 +++- unittest/force-styles/tests/bond-fene.yaml | 2 +- unittest/force-styles/tests/bond-fene_nm.yaml | 2 +- unittest/force-styles/tests/bond-gromos.yaml | 2 +- unittest/force-styles/tests/bond-harmonic.yaml | 2 +- 6 files changed, 9 insertions(+), 6 deletions(-) diff --git a/unittest/force-styles/tests/angle-cosine.yaml b/unittest/force-styles/tests/angle-cosine.yaml index 5a59fcce86..43629712d4 100644 --- a/unittest/force-styles/tests/angle-cosine.yaml +++ b/unittest/force-styles/tests/angle-cosine.yaml @@ -16,7 +16,8 @@ angle_coeff: ! | 3 50.0 4 100.0 equilibrium: 4 3.141592653589793 3.141592653589793 3.141592653589793 3.141592653589793 -extract: ! "" +extract: ! | + k 1 natoms: 29 init_energy: 1347.8670856939623 init_stress: ! |2- diff --git a/unittest/force-styles/tests/angle-harmonic.yaml b/unittest/force-styles/tests/angle-harmonic.yaml index dee700aa2c..d2164af8c5 100644 --- a/unittest/force-styles/tests/angle-harmonic.yaml +++ b/unittest/force-styles/tests/angle-harmonic.yaml @@ -16,7 +16,9 @@ angle_coeff: ! | 3 50.0 120.0 4 100.0 108.5 equilibrium: 4 1.9216075064457567 1.9373154697137058 2.0943951023931953 1.8936822384138476 -extract: ! "" +extract: ! | + k 1 + theta0 1 natoms: 29 init_energy: 41.53081789649104 init_stress: ! |2- diff --git a/unittest/force-styles/tests/bond-fene.yaml b/unittest/force-styles/tests/bond-fene.yaml index 88c02574cb..e5077eda0e 100644 --- a/unittest/force-styles/tests/bond-fene.yaml +++ b/unittest/force-styles/tests/bond-fene.yaml @@ -18,7 +18,7 @@ bond_coeff: ! | 5 450 2 0.018 1 equilibrium: 5 1.455 1.067 1.261 1.164 0.97 extract: ! | - kappa 1 + k 1 r0 1 natoms: 29 init_energy: 7104.900486467235 diff --git a/unittest/force-styles/tests/bond-fene_nm.yaml b/unittest/force-styles/tests/bond-fene_nm.yaml index c6be31a1c3..892d26c7aa 100644 --- a/unittest/force-styles/tests/bond-fene_nm.yaml +++ b/unittest/force-styles/tests/bond-fene_nm.yaml @@ -18,7 +18,7 @@ bond_coeff: ! | 5 450 2 0.018 1 12 6 equilibrium: 5 1.455 1.067 1.261 1.164 0.97 extract: ! | - kappa 1 + k 1 r0 1 natoms: 29 init_energy: 7104.538647187164 diff --git a/unittest/force-styles/tests/bond-gromos.yaml b/unittest/force-styles/tests/bond-gromos.yaml index a2f7e7ef3e..18abc99a3c 100644 --- a/unittest/force-styles/tests/bond-gromos.yaml +++ b/unittest/force-styles/tests/bond-gromos.yaml @@ -18,7 +18,7 @@ bond_coeff: ! | 5 450.0 1.0 equilibrium: 5 1.5 1.1 1.3 1.2 1 extract: ! | - kappa 1 + k 1 r0 1 natoms: 29 init_energy: 33.70930417641326 diff --git a/unittest/force-styles/tests/bond-harmonic.yaml b/unittest/force-styles/tests/bond-harmonic.yaml index 9ee14c07b9..bf686558d7 100644 --- a/unittest/force-styles/tests/bond-harmonic.yaml +++ b/unittest/force-styles/tests/bond-harmonic.yaml @@ -18,7 +18,7 @@ bond_coeff: ! | 5 450.0 1.0 equilibrium: 5 1.5 1.1 1.3 1.2 1 extract: ! | - kappa 1 + k 1 r0 1 natoms: 29 init_energy: 4.789374024601252 From c054edda6b2b4b317180ddb96d1485498095ac3e Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 22 Apr 2022 13:22:01 -0400 Subject: [PATCH 09/15] allow larger error margin for pressure computes --- unittest/commands/test_compute_global.cpp | 28 +++++++++++------------ 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/unittest/commands/test_compute_global.cpp b/unittest/commands/test_compute_global.cpp index 1c9be99ba4..e21acdbca0 100644 --- a/unittest/commands/test_compute_global.cpp +++ b/unittest/commands/test_compute_global.cpp @@ -102,24 +102,24 @@ TEST_F(ComputeGlobalTest, Energy) EXPECT_DOUBLE_EQ(get_scalar("pe1"), 24155.155261642241); EXPECT_DOUBLE_EQ(get_scalar("pe2"), 361.37528652881286); EXPECT_DOUBLE_EQ(get_scalar("pe3"), 0.0); - EXPECT_DOUBLE_EQ(get_scalar("pr1"), 1956948.4735454607); - EXPECT_DOUBLE_EQ(get_scalar("pr2"), 1956916.7725807722); + EXPECT_NEAR(get_scalar("pr1"), 1956948.4735454607, 0.0000000005); + EXPECT_NEAR(get_scalar("pr2"), 1956916.7725807722, 0.0000000005); EXPECT_DOUBLE_EQ(get_scalar("pr3"), 0.0); auto pr1 = get_vector("pr1"); auto pr2 = get_vector("pr2"); auto pr3 = get_vector("pr3"); - EXPECT_DOUBLE_EQ(pr1[0], 2150600.9207200543); - EXPECT_DOUBLE_EQ(pr1[1], 1466949.7512112649); - EXPECT_DOUBLE_EQ(pr1[2], 2253294.7487050635); - EXPECT_DOUBLE_EQ(pr1[3], 856643.16926486336); - EXPECT_DOUBLE_EQ(pr1[4], 692710.86929464422); - EXPECT_DOUBLE_EQ(pr1[5], -44403.909298603547); - EXPECT_DOUBLE_EQ(pr2[0], 2150575.6989334146); - EXPECT_DOUBLE_EQ(pr2[1], 1466911.3911461537); - EXPECT_DOUBLE_EQ(pr2[2], 2253263.2276627473); - EXPECT_DOUBLE_EQ(pr2[3], 856632.34707690508); - EXPECT_DOUBLE_EQ(pr2[4], 692712.89222328411); - EXPECT_DOUBLE_EQ(pr2[5], -44399.277068014424); + EXPECT_NEAR(pr1[0], 2150600.9207200543, 0.0000000005); + EXPECT_NEAR(pr1[1], 1466949.7512112649, 0.0000000005); + EXPECT_NEAR(pr1[2], 2253294.7487050635, 0.0000000005); + EXPECT_NEAR(pr1[3], 856643.16926486336, 0.0000000005); + EXPECT_NEAR(pr1[4], 692710.86929464422, 0.0000000005); + EXPECT_NEAR(pr1[5], -44403.909298603547, 0.0000000005); + EXPECT_NEAR(pr2[0], 2150575.6989334146, 0.0000000005); + EXPECT_NEAR(pr2[1], 1466911.3911461537, 0.0000000005); + EXPECT_NEAR(pr2[2], 2253263.2276627473, 0.0000000005); + EXPECT_NEAR(pr2[3], 856632.34707690508, 0.0000000005); + EXPECT_NEAR(pr2[4], 692712.89222328411, 0.0000000005); + EXPECT_NEAR(pr2[5], -44399.277068014424, 0.0000000005); EXPECT_DOUBLE_EQ(pr3[0], 0.0); EXPECT_DOUBLE_EQ(pr3[1], 0.0); EXPECT_DOUBLE_EQ(pr3[2], 0.0); From 1568974e8e3216d3dcef85e2c22cd5f4ddc5e828 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 22 Apr 2022 13:39:05 -0400 Subject: [PATCH 10/15] whitespace --- src/angle.cpp | 2 +- src/bond.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/angle.cpp b/src/angle.cpp index caa86d691e..8132146ba4 100644 --- a/src/angle.cpp +++ b/src/angle.cpp @@ -360,7 +360,7 @@ double Angle::memory_usage() void Angle::reinit() { - if (!reinitflag) + if (!reinitflag) error->all(FLERR, "Fix adapt interface to this angle style not supported"); init(); diff --git a/src/bond.cpp b/src/bond.cpp index 7f7b6aa3a3..a208b5a7e4 100644 --- a/src/bond.cpp +++ b/src/bond.cpp @@ -342,7 +342,7 @@ double Bond::memory_usage() void Bond::reinit() { - if (!reinitflag) + if (!reinitflag) error->all(FLERR, "Fix adapt interface to this bond style not supported"); init(); From fec5538d3cbbb1c00d987389a8565c8fb39ffd97 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 22 Apr 2022 13:52:15 -0400 Subject: [PATCH 11/15] fix initialization bugs --- src/fix_adapt.cpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/fix_adapt.cpp b/src/fix_adapt.cpp index 0d9aa00b7a..edaa5c9866 100644 --- a/src/fix_adapt.cpp +++ b/src/fix_adapt.cpp @@ -44,8 +44,9 @@ enum{DIAMETER,CHARGE}; /* ---------------------------------------------------------------------- */ -FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), -nadapt(0), id_fix_diam(nullptr), id_fix_chg(nullptr), adapt(nullptr) +FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg), nadapt(0), anypair(0), anybond(0), anyangle(0), + id_fix_diam(nullptr), id_fix_chg(nullptr), adapt(nullptr) { if (narg < 5) error->all(FLERR,"Illegal fix adapt command"); nevery = utils::inumeric(FLERR,arg[3],false,lmp); @@ -329,6 +330,7 @@ void FixAdapt::init() anypair = 0; anybond = 0; + anyangle = 0; for (int m = 0; m < nadapt; m++) { Adapt *ad = &adapt[m]; @@ -748,6 +750,7 @@ void FixAdapt::restore_settings() if (anypair) force->pair->reinit(); if (anybond) force->bond->reinit(); + if (anyangle) force->angle->reinit(); if (chgflag && force->kspace) force->kspace->qsum_qsq(); } From 01e75309026d34fb649e09f68f0512a9cc893ebc Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 22 Apr 2022 17:00:55 -0400 Subject: [PATCH 12/15] shorten URL message text --- src/utils.cpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/utils.cpp b/src/utils.cpp index 72944c1838..0d0bc91227 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -139,9 +139,7 @@ void utils::fmtargs_logmesg(LAMMPS *lmp, fmt::string_view format, fmt::format_ar std::string utils::errorurl(int errorcode) { - return fmt::format( - "\nFor more information please go to https://docs.lammps.org/Errors_details.html#err{:04d}", - errorcode); + return fmt::format("\nFor more information see https://docs.lammps.org/err{:04d}", errorcode); } void utils::flush_buffers(LAMMPS *lmp) From a90c632ae23c5b400662e501237494d92797845b Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Fri, 22 Apr 2022 17:10:56 -0600 Subject: [PATCH 13/15] fix bug when no other post_force fixes are defined --- src/fix_group.cpp | 1 - src/modify.cpp | 31 ++++++++++++++++++++----------- src/modify.h | 14 ++++++++------ src/respa.cpp | 5 +++-- src/verlet.cpp | 4 ++-- 5 files changed, 33 insertions(+), 22 deletions(-) diff --git a/src/fix_group.cpp b/src/fix_group.cpp index 112776251a..76b5bdbb4f 100644 --- a/src/fix_group.cpp +++ b/src/fix_group.cpp @@ -186,7 +186,6 @@ void FixGroup::set_group() // change the compute result since it will not be re-invoked at end-of-step, // e.g. if compute pe/atom includes pe contributions from fixes - double *var = nullptr; int *ivector = nullptr; double *dvector = nullptr; diff --git a/src/modify.cpp b/src/modify.cpp index ec74b03554..9a359e1971 100644 --- a/src/modify.cpp +++ b/src/modify.cpp @@ -51,22 +51,22 @@ Modify::Modify(LAMMPS *lmp) : Pointers(lmp) nfix = maxfix = 0; n_initial_integrate = n_post_integrate = 0; n_pre_exchange = n_pre_neighbor = n_post_neighbor = 0; - n_pre_force = n_pre_reverse = n_post_force = 0; + n_pre_force = n_pre_reverse = n_post_force_any = 0; n_final_integrate = n_end_of_step = 0; n_energy_couple = n_energy_global = n_energy_atom = 0; n_initial_integrate_respa = n_post_integrate_respa = 0; - n_pre_force_respa = n_post_force_respa = n_final_integrate_respa = 0; + n_pre_force_respa = n_post_force_respa_any = n_final_integrate_respa = 0; n_min_pre_exchange = n_min_pre_force = n_min_pre_reverse = 0; n_min_post_force = n_min_energy = 0; - n_post_force_group = 0; n_timeflag = -1; fix = nullptr; fmask = nullptr; list_initial_integrate = list_post_integrate = nullptr; list_pre_exchange = list_pre_neighbor = list_post_neighbor = nullptr; - list_pre_force = list_pre_reverse = list_post_force = nullptr; + list_pre_force = list_pre_reverse = nullptr; + list_post_force = list_post_force_group = nullptr; list_final_integrate = list_end_of_step = nullptr; list_energy_couple = list_energy_global = list_energy_atom = nullptr; list_initial_integrate_respa = list_post_integrate_respa = nullptr; @@ -78,8 +78,6 @@ Modify::Modify(LAMMPS *lmp) : Pointers(lmp) end_of_step_every = nullptr; - list_post_force_group = nullptr; - list_timeflag = nullptr; nfix_restart_global = 0; @@ -143,6 +141,7 @@ Modify::~Modify() delete[] list_pre_force; delete[] list_pre_reverse; delete[] list_post_force; + delete[] list_post_force_group; delete[] list_final_integrate; delete[] list_end_of_step; delete[] list_energy_couple; @@ -160,7 +159,6 @@ Modify::~Modify() delete[] list_min_pre_reverse; delete[] list_min_post_force; delete[] list_min_energy; - delete[] list_post_force_group; delete[] end_of_step_every; delete[] list_timeflag; @@ -226,12 +224,12 @@ void Modify::init() list_init(PRE_FORCE, n_pre_force, list_pre_force); list_init(PRE_REVERSE, n_pre_reverse, list_pre_reverse); list_init(POST_FORCE, n_post_force, list_post_force); + list_init_post_force_group(n_post_force_group, list_post_force_group); list_init(FINAL_INTEGRATE, n_final_integrate, list_final_integrate); list_init_end_of_step(END_OF_STEP, n_end_of_step, list_end_of_step); list_init_energy_couple(n_energy_couple, list_energy_couple); list_init_energy_global(n_energy_global, list_energy_global); list_init_energy_atom(n_energy_atom, list_energy_atom); - list_init_post_force_group(n_post_force_group, list_post_force_group); list_init(INITIAL_INTEGRATE_RESPA, n_initial_integrate_respa, list_initial_integrate_respa); list_init(POST_INTEGRATE_RESPA, n_post_integrate_respa, list_post_integrate_respa); @@ -247,6 +245,11 @@ void Modify::init() list_init(MIN_POST_FORCE, n_min_post_force, list_min_post_force); list_init(MIN_ENERGY, n_min_energy, list_min_energy); + // two post_force_any counters used by integrators add in post_force_group + + n_post_force_any = n_post_force + n_post_force_group; + n_post_force_respa_any = n_post_force_respa + n_post_force_group; + // create list of computes that store invocation times list_init_compute(); @@ -448,6 +451,7 @@ void Modify::pre_reverse(int eflag, int vflag) /* ---------------------------------------------------------------------- post_force call, only for relevant fixes first call any instances of fix GROUP if they exist + they are not in n_post_force count ------------------------------------------------------------------------- */ void Modify::post_force(int vflag) @@ -457,7 +461,10 @@ void Modify::post_force(int vflag) fix[list_post_force_group[i]]->post_force(vflag); } - for (int i = 0; i < n_post_force; i++) fix[list_post_force[i]]->post_force(vflag); + if (n_post_force) { + for (int i = 0; i < n_post_force; i++) + fix[list_post_force[i]]->post_force(vflag); + } } /* ---------------------------------------------------------------------- @@ -607,8 +614,10 @@ void Modify::post_force_respa(int vflag, int ilevel, int iloop) fix[list_post_force_group[i]]->post_force_respa(vflag, ilevel, iloop); } - for (int i = 0; i < n_post_force_respa; i++) - fix[list_post_force_respa[i]]->post_force_respa(vflag, ilevel, iloop); + if (n_post_force_respa) { + for (int i = 0; i < n_post_force_respa; i++) + fix[list_post_force_respa[i]]->post_force_respa(vflag, ilevel, iloop); + } } /* ---------------------------------------------------------------------- diff --git a/src/modify.h b/src/modify.h index da5d00a799..5e05cfce6a 100644 --- a/src/modify.h +++ b/src/modify.h @@ -33,11 +33,11 @@ class Modify : protected Pointers { public: int n_initial_integrate, n_post_integrate, n_pre_exchange; int n_pre_neighbor, n_post_neighbor; - int n_pre_force, n_pre_reverse, n_post_force; + int n_pre_force, n_pre_reverse, n_post_force_any; int n_final_integrate, n_end_of_step; int n_energy_couple, n_energy_global, n_energy_atom; int n_initial_integrate_respa, n_post_integrate_respa; - int n_pre_force_respa, n_post_force_respa, n_final_integrate_respa; + int n_pre_force_respa, n_post_force_respa_any, n_final_integrate_respa; int n_min_pre_exchange, n_min_pre_neighbor, n_min_post_neighbor; int n_min_pre_force, n_min_pre_reverse, n_min_post_force, n_min_energy; @@ -147,11 +147,16 @@ class Modify : protected Pointers { double memory_usage(); protected: + // internal fix counts + + int n_post_force, n_post_force_group, n_post_force_respa; + // lists of fixes to apply at different stages of timestep int *list_initial_integrate, *list_post_integrate; int *list_pre_exchange, *list_pre_neighbor, *list_post_neighbor; - int *list_pre_force, *list_pre_reverse, *list_post_force; + int *list_pre_force, *list_pre_reverse; + int *list_post_force, *list_post_force_group; int *list_final_integrate, *list_end_of_step; int *list_energy_couple, *list_energy_global, *list_energy_atom; int *list_initial_integrate_respa, *list_post_integrate_respa; @@ -163,9 +168,6 @@ class Modify : protected Pointers { int *end_of_step_every; - int n_post_force_group; // list of fix GROUPs for post_force invocation - int *list_post_force_group; - int n_timeflag; // list of computes that store time invocation int *list_timeflag; diff --git a/src/respa.cpp b/src/respa.cpp index 97356eac47..50d5e3250e 100644 --- a/src/respa.cpp +++ b/src/respa.cpp @@ -372,7 +372,7 @@ void Respa::setup(int flag) mesg += fmt::format(" {}:{}", ilevel + 1, step[ilevel]); mesg += "\n r-RESPA fixes :"; - for (int l = 0; l < modify->n_post_force_respa; ++l) { + for (int l = 0; l < modify->n_post_force_respa_any; ++l) { Fix *f = modify->get_fix_by_index(modify->list_post_force_respa[l]); if (f->respa_level >= 0) mesg += fmt::format(" {}:{}[{}]", MIN(f->respa_level + 1, nlevels), f->style, f->id); @@ -704,7 +704,8 @@ void Respa::recurse(int ilevel) timer->stamp(Timer::COMM); } timer->stamp(); - if (modify->n_post_force_respa) modify->post_force_respa(vflag, ilevel, iloop); + if (modify->n_post_force_respa_any) + modify->post_force_respa(vflag, ilevel, iloop); modify->final_integrate_respa(ilevel, iloop); timer->stamp(Timer::MODIFY); } diff --git a/src/verlet.cpp b/src/verlet.cpp index 342dc3d951..15e1ab0800 100644 --- a/src/verlet.cpp +++ b/src/verlet.cpp @@ -237,7 +237,7 @@ void Verlet::run(int n) int n_post_neighbor = modify->n_post_neighbor; int n_pre_force = modify->n_pre_force; int n_pre_reverse = modify->n_pre_reverse; - int n_post_force = modify->n_post_force; + int n_post_force_any = modify->n_post_force_any; int n_end_of_step = modify->n_end_of_step; if (atom->sortfreq > 0) sortflag = 1; @@ -344,7 +344,7 @@ void Verlet::run(int n) // force modifications, final time integration, diagnostics - if (n_post_force) modify->post_force(vflag); + if (n_post_force_any) modify->post_force(vflag); modify->final_integrate(); if (n_end_of_step) modify->end_of_step(); timer->stamp(Timer::MODIFY); From 3a5ab301b2e5a2e6ba9ef17d522542151116249b Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 22 Apr 2022 21:46:24 -0400 Subject: [PATCH 14/15] propagate changes from previous commit to other packages that were missing them --- src/INTEL/verlet_lrt_intel.cpp | 2 +- src/KOKKOS/dynamical_matrix_kokkos.cpp | 2 +- src/KOKKOS/third_order_kokkos.cpp | 2 +- src/KOKKOS/verlet_kokkos.cpp | 2 +- src/MC/fix_atom_swap.cpp | 2 +- src/MC/fix_charge_regulation.cpp | 2 +- src/MC/fix_gcmc.cpp | 2 +- src/MC/fix_mol_swap.cpp | 2 +- src/OPENMP/respa_omp.cpp | 4 ++-- src/PHONON/dynamical_matrix.cpp | 2 +- src/PHONON/third_order.cpp | 2 +- src/REPLICA/verlet_split.cpp | 2 +- 12 files changed, 13 insertions(+), 13 deletions(-) diff --git a/src/INTEL/verlet_lrt_intel.cpp b/src/INTEL/verlet_lrt_intel.cpp index 0d8eb2e468..4102510c66 100644 --- a/src/INTEL/verlet_lrt_intel.cpp +++ b/src/INTEL/verlet_lrt_intel.cpp @@ -215,7 +215,7 @@ void VerletLRTIntel::run(int n) int n_pre_neighbor = modify->n_pre_neighbor; int n_pre_force = modify->n_pre_force; int n_pre_reverse = modify->n_pre_reverse; - int n_post_force = modify->n_post_force; + int n_post_force = modify->n_post_force_any; int n_end_of_step = modify->n_end_of_step; if (atom->sortfreq > 0) sortflag = 1; diff --git a/src/KOKKOS/dynamical_matrix_kokkos.cpp b/src/KOKKOS/dynamical_matrix_kokkos.cpp index dc180a1743..7b85a3766d 100644 --- a/src/KOKKOS/dynamical_matrix_kokkos.cpp +++ b/src/KOKKOS/dynamical_matrix_kokkos.cpp @@ -159,7 +159,7 @@ void DynamicalMatrixKokkos::update_force() { int n_pre_force = modify->n_pre_force; int n_pre_reverse = modify->n_pre_reverse; - int n_post_force = modify->n_post_force; + int n_post_force = modify->n_post_force_any; lmp->kokkos->auto_sync = 0; diff --git a/src/KOKKOS/third_order_kokkos.cpp b/src/KOKKOS/third_order_kokkos.cpp index 2aeb9152a1..56c9295e70 100644 --- a/src/KOKKOS/third_order_kokkos.cpp +++ b/src/KOKKOS/third_order_kokkos.cpp @@ -160,7 +160,7 @@ void ThirdOrderKokkos::update_force() { int n_pre_force = modify->n_pre_force; int n_pre_reverse = modify->n_pre_reverse; - int n_post_force = modify->n_post_force; + int n_post_force = modify->n_post_force_any; lmp->kokkos->auto_sync = 0; diff --git a/src/KOKKOS/verlet_kokkos.cpp b/src/KOKKOS/verlet_kokkos.cpp index 2d142edd83..1816415b47 100644 --- a/src/KOKKOS/verlet_kokkos.cpp +++ b/src/KOKKOS/verlet_kokkos.cpp @@ -271,7 +271,7 @@ void VerletKokkos::run(int n) int n_post_neighbor = modify->n_post_neighbor; int n_pre_force = modify->n_pre_force; int n_pre_reverse = modify->n_pre_reverse; - int n_post_force = modify->n_post_force; + int n_post_force = modify->n_post_force_any; int n_end_of_step = modify->n_end_of_step; lmp->kokkos->auto_sync = 0; diff --git a/src/MC/fix_atom_swap.cpp b/src/MC/fix_atom_swap.cpp index 690e628fe3..8405919db8 100644 --- a/src/MC/fix_atom_swap.cpp +++ b/src/MC/fix_atom_swap.cpp @@ -532,7 +532,7 @@ double FixAtomSwap::energy_full() if (force->kspace) force->kspace->compute(eflag,vflag); - if (modify->n_post_force) modify->post_force(vflag); + if (modify->n_post_force_any) modify->post_force(vflag); update->eflag_global = update->ntimestep; double total_energy = c_pe->compute_scalar(); diff --git a/src/MC/fix_charge_regulation.cpp b/src/MC/fix_charge_regulation.cpp index a5a2eaff82..81f722bb6c 100644 --- a/src/MC/fix_charge_regulation.cpp +++ b/src/MC/fix_charge_regulation.cpp @@ -1129,7 +1129,7 @@ double FixChargeRegulation::energy_full() { if (force->kspace) force->kspace->compute(eflag, vflag); if (modify->n_pre_reverse) modify->pre_reverse(eflag,vflag); - if (modify->n_post_force) modify->post_force(vflag); + if (modify->n_post_force_any) modify->post_force(vflag); update->eflag_global = update->ntimestep; double total_energy = c_pe->compute_scalar(); diff --git a/src/MC/fix_gcmc.cpp b/src/MC/fix_gcmc.cpp index 6937b90202..a526ba4756 100644 --- a/src/MC/fix_gcmc.cpp +++ b/src/MC/fix_gcmc.cpp @@ -2316,7 +2316,7 @@ double FixGCMC::energy_full() // but Modify::pre_reverse() is needed for INTEL if (modify->n_pre_reverse) modify->pre_reverse(eflag,vflag); - if (modify->n_post_force) modify->post_force(vflag); + if (modify->n_post_force_any) modify->post_force(vflag); // NOTE: all fixes with energy_global_flag set and which // operate at pre_force() or post_force() diff --git a/src/MC/fix_mol_swap.cpp b/src/MC/fix_mol_swap.cpp index c93f0093a8..c23ceb7d36 100644 --- a/src/MC/fix_mol_swap.cpp +++ b/src/MC/fix_mol_swap.cpp @@ -390,7 +390,7 @@ double FixMolSwap::energy_full() if (force->kspace) force->kspace->compute(eflag,vflag); - if (modify->n_post_force) modify->post_force(vflag); + if (modify->n_post_force_any) modify->post_force(vflag); update->eflag_global = update->ntimestep; double total_energy = c_pe->compute_scalar(); diff --git a/src/OPENMP/respa_omp.cpp b/src/OPENMP/respa_omp.cpp index c6500185ad..93eb3c41ef 100644 --- a/src/OPENMP/respa_omp.cpp +++ b/src/OPENMP/respa_omp.cpp @@ -77,7 +77,7 @@ void RespaOMP::setup(int flag) mesg += fmt::format(" {}:{}", ilevel + 1, step[ilevel]); mesg += "\n r-RESPA fixes :"; - for (int l = 0; l < modify->n_post_force_respa; ++l) { + for (int l = 0; l < modify->n_post_force_respa_any; ++l) { Fix *f = modify->get_fix_by_index(modify->list_post_force_respa[l]); if (f->respa_level >= 0) mesg += fmt::format(" {}:{}[{}]", MIN(f->respa_level + 1, nlevels), f->style, f->id); @@ -420,7 +420,7 @@ void RespaOMP::recurse(int ilevel) timer->stamp(Timer::COMM); } timer->stamp(); - if (modify->n_post_force_respa) + if (modify->n_post_force_respa_any) modify->post_force_respa(vflag,ilevel,iloop); modify->final_integrate_respa(ilevel,iloop); timer->stamp(Timer::MODIFY); diff --git a/src/PHONON/dynamical_matrix.cpp b/src/PHONON/dynamical_matrix.cpp index 3723a7a467..8f667f51f2 100644 --- a/src/PHONON/dynamical_matrix.cpp +++ b/src/PHONON/dynamical_matrix.cpp @@ -431,7 +431,7 @@ void DynamicalMatrix::update_force() force_clear(); int n_pre_force = modify->n_pre_force; int n_pre_reverse = modify->n_pre_reverse; - int n_post_force = modify->n_post_force; + int n_post_force = modify->n_post_force_any; if (n_pre_force) { modify->pre_force(vflag); diff --git a/src/PHONON/third_order.cpp b/src/PHONON/third_order.cpp index 76d90c65bd..bbd6818b95 100644 --- a/src/PHONON/third_order.cpp +++ b/src/PHONON/third_order.cpp @@ -487,7 +487,7 @@ void ThirdOrder::update_force() neighbor->ago = 0; if (modify->get_fix_by_id("package_intel")) neighbor->decide(); force_clear(); - int n_post_force = modify->n_post_force; + int n_post_force = modify->n_post_force_any; int n_pre_force = modify->n_pre_force; int n_pre_reverse = modify->n_pre_reverse; diff --git a/src/REPLICA/verlet_split.cpp b/src/REPLICA/verlet_split.cpp index 95cd54119c..661b86129c 100644 --- a/src/REPLICA/verlet_split.cpp +++ b/src/REPLICA/verlet_split.cpp @@ -300,7 +300,7 @@ void VerletSplit::run(int n) int n_pre_neighbor = modify->n_pre_neighbor; int n_pre_force = modify->n_pre_force; int n_pre_reverse = modify->n_pre_reverse; - int n_post_force = modify->n_post_force; + int n_post_force = modify->n_post_force_any; int n_end_of_step = modify->n_end_of_step; if (atom->sortfreq > 0) sortflag = 1; From 1319cb2cf5b12d36b7e2f599dda5d26d7f3f89e9 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 22 Apr 2022 22:08:13 -0400 Subject: [PATCH 15/15] fix typo --- doc/src/fix_nve.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/fix_nve.rst b/doc/src/fix_nve.rst index a0e186df29..457ee62133 100644 --- a/doc/src/fix_nve.rst +++ b/doc/src/fix_nve.rst @@ -41,7 +41,7 @@ time integration options can be invoked using the :doc:`run_style ` c ---------- -.. Include:: accel_styles.rst +.. include:: accel_styles.rst ----------