Merge branch 'develop' into safer-data-file-parsing

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:: missing_cmd_args
    :project: progguide
 

doc/src/Errors.rst

@@ -11,6 +11,7 @@ them.
    :maxdepth: 1
 
    Errors_common
+   Errors_details
    Errors_bugs
    Errors_debug
    Errors_messages

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.
+

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
+       *atom* args = atomparam v_name
-         aparam = parameter to adapt over time
+         atomparam = parameter to adapt over time
-         v_name = variable with name that calculates value of aparam
+         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,62 @@ 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.
+the bond type argument to set the coefficients for multiple bond
-This takes the form "\*" or "\*n" or "n\*" or "m\*n".  If N = the number of
+types.  This takes the form "\*" or "\*n" or "n\*" or "m\*n".  If N =
-atom types, then an asterisk with no numeric values means all types
+the number of bond types, then an asterisk with no numeric values
-from 1 to N.  A leading asterisk means all types from 1 to n (inclusive).
+means all types from 1 to N.  A leading asterisk means all types from
-A trailing asterisk means all types from n to N (inclusive).  A middle
+1 to n (inclusive).  A trailing asterisk means all types from n to N
-asterisk means all types from m to n (inclusive).
+(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
+hybrid/overlay as bond styles. The bond styles that currently work
-working with fix_adapt are
+with fix_adapt are
 
-+------------------------------------+-------+------------+
++------------------------------------+-------+-----------------+
-| :doc:`class2 <bond_class2>`        | r0    | type bonds |
+| :doc:`class2 <bond_class2>`        | r0         | type bonds |
-+------------------------------------+-------+------------+
++------------------------------------+-------+-----------------+
-| :doc:`fene <bond_fene>`            | k, r0 | type bonds |
+| :doc:`fene <bond_fene>`            | k,r0       | type bonds |
-+------------------------------------+-------+------------+
++------------------------------------+-------+-----------------+
-| :doc:`gromos <bond_gromos>`        | k, r0 | type bonds |
+| :doc:`fene/nm <bond_fene_nm>`      | k,r0       | type bonds |
-+------------------------------------+-------+------------+
++------------------------------------+-------+-----------------+
-| :doc:`harmonic <bond_harmonic>`    | k,r0  | type bonds |
+| :doc:`gromos <bond_gromos>`        | k,r0       | type bonds |
-+------------------------------------+-------+------------+
++------------------------------------+-------+-----------------+
-| :doc:`morse <bond_morse>`          | r0    | type bonds |
+| :doc:`harmonic <bond_harmonic>`    | k,r0       | type bonds |
-+------------------------------------+-------+------------+
++------------------------------------+-------+-----------------+
-| :doc:`nonlinear <bond_nonlinear>`  | r0    | type bonds |
+| :doc:`morse <bond_morse>`          | r0         | type bonds |
-+------------------------------------+-------+------------+
++------------------------------------+-------+-----------------+
+| :doc:`nonlinear <bond_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 angle styles that currently work
+with fix_adapt are
+
++------------------------------------+-------+-----------------+
+| :doc:`harmonic <angle_harmonic>`    | k,theta0 | type angles |
++------------------------------------+-------+-----------------+
+| :doc:`cosine <angle_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.
 
 ----------
 

doc/src/fix_nve.rst

@@ -35,6 +35,10 @@ 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).
 
+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 <run_style>` command.
+
 ----------
 
 .. include:: accel_styles.rst
@@ -57,7 +61,7 @@ Restrictions
 Related commands
 """"""""""""""""
 
-:doc:`fix nvt <fix_nh>`, :doc:`fix npt <fix_nh>`
+:doc:`fix nvt <fix_nh>`, :doc:`fix npt <fix_nh>`, :doc:`run_style <run_style>`
 
 Default
 """""""

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
+group is after interatomic forces have been computed, but before any
-has been performed, and before neighbor lists or forces are computed.
+fixes which alter forces or otherwise update the system have been
-This is the point in the timestep where atom positions have just
+invoked.  This means that atom positions have been updated, neighbor
-changed due to the time integration, so the region criterion should be
+lists and ghost atoms are current, and both intermolecular and
-accurate, if applied.
+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::
 

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)
 
 ----------
 

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;
 }

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;
-}
 
 /* ---------------------------------------------------------------------- */
 

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;

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;
 

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;
 

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;

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();

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();

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()

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();

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;
+}

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;

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;
+}

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;

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;
 }

src/MOLECULE/bond_gromos.cpp

@@ -30,10 +30,7 @@ using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-BondGromos::BondGromos(LAMMPS *_lmp) : Bond(_lmp)
+BondGromos::BondGromos(LAMMPS *_lmp) : Bond(_lmp) {}
-{
-  reinitflag = 1;
-}
 
 /* ---------------------------------------------------------------------- */
 
@@ -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;
 }

src/MOLECULE/bond_harmonic.cpp

@@ -27,10 +27,7 @@ using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-BondHarmonic::BondHarmonic(LAMMPS *_lmp) : Bond(_lmp)
+BondHarmonic::BondHarmonic(LAMMPS *_lmp) : Bond(_lmp) {}
-{
-  reinitflag = 1;
-}
 
 /* ---------------------------------------------------------------------- */
 
@@ -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;
 }

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);

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);

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;
 

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;

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();
+}

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

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();
 }

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;}

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 ComputeClusterAtom::pack_forward_comm(int n, int *list, double *buf, 
-                                          int * /*pbc*/)
+                                          int /*pbc_flag*/, int * /*pbc*/)
 {
   int i, j, m;
 
   m = 0;
-  if (commflag == CLUSTER) {
+  for (i = 0; i < n; i++) {
-    for (i = 0; i < n; i++) {
+    j = list[i];
-      j = list[i];
+    buf[m++] = clusterID[j];
-      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];
-    }
   }
 
   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++];
-    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++];
-    }
-  }
 }
 
 /* ----------------------------------------------------------------------

src/compute_coord_atom.cpp

@@ -259,14 +259,16 @@ void ComputeCoordAtom::compute_peratom()
             j = jlist[jj];
             j &= NEIGHMASK;
 
-            jtype = type[j];
+            if (mask[j] & jgroupbit) {
-            delx = xtmp - x[j][0];
+              jtype = type[j];
-            dely = ytmp - x[j][1];
+              delx = xtmp - x[j][0];
-            delz = ztmp - x[j][2];
+              dely = ytmp - x[j][1];
-            rsq = delx * delx + dely * dely + delz * delz;
+              delz = ztmp - x[j][2];
-            if (rsq < cutsq) {
+              rsq = delx * delx + dely * dely + delz * delz;
-              for (m = 0; m < ncol; m++)
+              if (rsq < cutsq) {
-                if (jtype >= typelo[m] && jtype <= typehi[m]) count[m] += 1.0;
+                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 ComputeCoordAtom::pack_forward_comm(int n, int *list, double *buf, 
-                                        int * /*pbc*/)
+                                        int /*pbc_flag*/, int * /*pbc*/)
 {
   int i, m = 0, j;
   for (i = 0; i < n; ++i) {

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,13 +39,14 @@ using namespace LAMMPS_NS;
 using namespace FixConst;
 using namespace MathConst;
 
-enum{PAIR,KSPACE,ATOM,BOND};
+enum{PAIR,KSPACE,ATOM,BOND,ANGLE};
 enum{DIAMETER,CHARGE};
 
 /* ---------------------------------------------------------------------- */
 
-FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg),
+FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) :
-nadapt(0), id_fix_diam(nullptr), id_fix_chg(nullptr), adapt(nullptr)
+  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);
@@ -75,6 +77,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 +125,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 +153,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 +211,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 +234,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;
@@ -299,6 +330,7 @@ void FixAdapt::init()
 
   anypair = 0;
   anybond = 0;
+  anyangle = 0;
 
   for (int m = 0; m < nadapt; m++) {
     Adapt *ad = &adapt[m];
@@ -357,6 +389,7 @@ void FixAdapt::init()
       }
 
       delete [] pstyle;
+
     } else if (ad->which == BOND) {
       ad->bond = nullptr;
       anybond = 1;
@@ -383,13 +416,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 +457,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 +467,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 +481,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 +546,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 +590,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 +615,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 +642,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 +666,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 +676,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 +700,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];
       }
@@ -668,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();
 }
 

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;

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,22 +181,20 @@ void FixGroup::set_group()
   int nlocal = atom->nlocal;
 
   // invoke atom-style variable if defined
-  // set post_integrate flag to 1, then unset after
+  // NOTE: after variable invocation could reset invoked computes to not-invoked
-  // this is for any compute to check if it needs to
+  //   this would avoid an issue where other post-force fixes 
-  //   operate differently due to invocation this early in timestep
+  //   change the compute result since it will not be re-invoked at end-of-step,
-  // e.g. perform ghost comm update due to atoms having just moved
+  //   e.g. if compute pe/atom includes pe contributions from fixes
 
   double *var = nullptr;
   int *ivector = nullptr;
   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 +208,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 +229,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;
 }
 
 /* ---------------------------------------------------------------------- */

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_force(int) override;
-  void post_integrate_respa(int, 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:

src/modify.cpp

@@ -51,20 +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_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;
@@ -139,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;
@@ -221,6 +224,7 @@ 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);
@@ -241,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();
@@ -441,11 +450,21 @@ 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)
 {
-  for (int i = 0; i < n_post_force; i++) fix[list_post_force[i]]->post_force(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);
+  }
+
+  if (n_post_force) {
+    for (int i = 0; i < n_post_force; i++) 
+      fix[list_post_force[i]]->post_force(vflag);
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -585,12 +604,20 @@ 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)
 {
-  for (int i = 0; i < n_post_force_respa; i++)
+  if (n_post_force_group) {
-    fix[list_post_force_respa[i]]->post_force_respa(vflag, ilevel, iloop);
+    for (int i = 0; i < n_post_force_group; i++) 
+      fix[list_post_force_group[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);
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -1716,6 +1743,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
 ------------------------------------------------------------------------- */

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;
@@ -187,6 +192,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();
 

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);
     }

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);
   }

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;

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
 

src/utils.cpp

@@ -143,6 +143,11 @@ 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 see https://docs.lammps.org/err{:04d}", errorcode);
+}
+
 void utils::flush_buffers(LAMMPS *lmp)
 {
   if (lmp->screen) fflush(lmp->screen);

src/utils.h

@@ -86,6 +86,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

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);

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_NEAR(get_scalar("pr1"), 1956948.4735454607, 0.0000000005);
-    EXPECT_DOUBLE_EQ(get_scalar("pr2"), 1956916.7725807722);
+    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_NEAR(pr1[0], 2150600.9207200543, 0.0000000005);
-    EXPECT_DOUBLE_EQ(pr1[1], 1466949.7512112649);
+    EXPECT_NEAR(pr1[1], 1466949.7512112649, 0.0000000005);
-    EXPECT_DOUBLE_EQ(pr1[2], 2253294.7487050635);
+    EXPECT_NEAR(pr1[2], 2253294.7487050635, 0.0000000005);
-    EXPECT_DOUBLE_EQ(pr1[3], 856643.16926486336);
+    EXPECT_NEAR(pr1[3], 856643.16926486336, 0.0000000005);
-    EXPECT_DOUBLE_EQ(pr1[4], 692710.86929464422);
+    EXPECT_NEAR(pr1[4], 692710.86929464422, 0.0000000005);
-    EXPECT_DOUBLE_EQ(pr1[5], -44403.909298603547);
+    EXPECT_NEAR(pr1[5], -44403.909298603547, 0.0000000005);
-    EXPECT_DOUBLE_EQ(pr2[0], 2150575.6989334146);
+    EXPECT_NEAR(pr2[0], 2150575.6989334146, 0.0000000005);
-    EXPECT_DOUBLE_EQ(pr2[1], 1466911.3911461537);
+    EXPECT_NEAR(pr2[1], 1466911.3911461537, 0.0000000005);
-    EXPECT_DOUBLE_EQ(pr2[2], 2253263.2276627473);
+    EXPECT_NEAR(pr2[2], 2253263.2276627473, 0.0000000005);
-    EXPECT_DOUBLE_EQ(pr2[3], 856632.34707690508);
+    EXPECT_NEAR(pr2[3], 856632.34707690508, 0.0000000005);
-    EXPECT_DOUBLE_EQ(pr2[4], 692712.89222328411);
+    EXPECT_NEAR(pr2[4], 692712.89222328411, 0.0000000005);
-    EXPECT_DOUBLE_EQ(pr2[5], -44399.277068014424);
+    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);

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-

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-

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

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

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

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