initial implementation of minimizer support in fix shake/rattle

doc/src/fix_shake.rst

@@ -33,12 +33,14 @@ Syntax
        *m* value = one or more mass values
 
 * zero or more keyword/value pairs may be appended
-* keyword = *mol*
+* keyword = *mol* or *kbond*
 
   .. parsed-literal::
 
        *mol* value = template-ID
          template-ID = ID of molecule template specified in a separate :doc:`molecule <molecule>` command
+       *kbond* value = force constant
+         force constant = force constant used to apply a restraint force when used during minimization
 
 Examples
 """"""""
@@ -152,17 +154,23 @@ for.
 
 ----------
 
-The *mol* keyword should be used when other commands, such as :doc:`fix deposit <fix_deposit>` or :doc:`fix pour <fix_pour>`, add molecules
+The *mol* keyword should be used when other commands, such as :doc:`fix
+deposit <fix_deposit>` or :doc:`fix pour <fix_pour>`, add molecules
 on-the-fly during a simulation, and you wish to constrain the new
 molecules via SHAKE.  You specify a *template-ID* previously defined
 using the :doc:`molecule <molecule>` command, which reads a file that
 defines the molecule.  You must use the same *template-ID* that the
 command adding molecules uses.  The coordinates, atom types, special
-bond restrictions, and SHAKE info can be specified in the molecule
+bond restrictions, and SHAKE info can be specified in the molecule file.
-file.  See the :doc:`molecule <molecule>` command for details.  The only
+See the :doc:`molecule <molecule>` command for details.  The only
 settings required to be in this file (by this command) are the SHAKE
 info of atoms in the molecule.
 
+The *kbond* keyword allows to set the restraint force constant when
+fix shake or fix rattle are used during minimization. In that case
+the constraint algorithms are **not** applied and restraint
+forces are used instead to help maintaining the geometries.
+
 ----------
 
 .. include:: accel_styles.rst
@@ -205,8 +213,10 @@ setting for this fix is :doc:`fix_modify virial yes <fix_modify>`.
 No global or per-atom quantities are stored by these fixes for access
 by various :doc:`output commands <Howto_output>`.  No parameter of
 these fixes can be used with the *start/stop* keywords of the
-:doc:`run <run>` command.  These fixes are not invoked during
+:doc:`run <run>` command.
-:doc:`energy minimization <minimize>`.
+
+When used during minimization, the SHAKE or RATTLE algorithms are **not**
+applied.  Strong restraint forces are applied instead.
 
 Restrictions
 """"""""""""
@@ -232,13 +242,14 @@ make a linear molecule rigid.
 Related commands
 """"""""""""""""
 
-none
+`fix rigid <fix_rigid>`, `fix ehex <fix_ehex>`,
+`fix nve/manifold/rattle <fix_nve_manifold_rattle>`
 
 
 Default
 """""""
 
-none
+kbond = 1.0e6
 
 ----------
 

src/KOKKOS/fix_shake_kokkos.cpp

@@ -12,12 +12,8 @@
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
-#include <mpi.h>
-#include <cmath>
-#include <cstdlib>
-#include <cstring>
-#include <cstdio>
 #include "fix_shake_kokkos.h"
+
 #include "fix_rattle.h"
 #include "atom_kokkos.h"
 #include "atom_vec.h"
@@ -38,6 +34,9 @@
 #include "kokkos.h"
 #include "atom_masks.h"
 
+#include <cmath>
+#include <cstring>
+
 using namespace LAMMPS_NS;
 using namespace FixConst;
 using namespace MathConst;
@@ -292,7 +291,7 @@ void FixShakeKokkos<DeviceType>::pre_neighbor()
 
   if (h_error_flag() == 1) {
     error->one(FLERR,"Shake atoms missing on proc "
-                                 "{} at step {}",me,update->ntimestep);
+                                 "{} at step {}",comm->me,update->ntimestep);
   }
 }
 
@@ -341,7 +340,7 @@ void FixShakeKokkos<DeviceType>::post_force(int vflag)
   // communicate results if necessary
 
   unconstrained_update();
-  if (nprocs > 1) comm->forward_comm(this);
+  if (comm->nprocs > 1) comm->forward_comm(this);
   k_xshake.sync<DeviceType>();
 
   // virial setup
@@ -1702,7 +1701,7 @@ void FixShakeKokkos<DeviceType>::correct_coordinates(int vflag) {
 
   double **xtmp = xshake;
   xshake = x;
-  if (nprocs > 1) {
+  if (comm->nprocs > 1) {
     forward_comm_device = 0;
     comm->forward_comm(this);
     forward_comm_device = 1;

src/RIGID/fix_rattle.cpp

@@ -81,7 +81,7 @@ FixRattle::~FixRattle()
 {
   memory->destroy(vp);
 
-  if (RATTLE_DEBUG) {
+#if RATTLE_DEBUG
 
     // communicate maximum distance error
 
@@ -91,13 +91,11 @@ FixRattle::~FixRattle()
     MPI_Reduce(&derr_max, &global_derr_max, 1 , MPI_DOUBLE, MPI_MAX, 0, world);
     MPI_Reduce(&verr_max, &global_verr_max, 1 , MPI_DOUBLE, MPI_MAX, 0, world);
 
-    MPI_Comm_rank (world, &npid); // Find out process rank
+    if (comm->me == 0 && screen) {
-
-    if (npid == 0 && screen) {
       fprintf(screen, "RATTLE: Maximum overall relative position error ( (r_ij-d_ij)/d_ij ): %.10g\n", global_derr_max);
       fprintf(screen, "RATTLE: Maximum overall absolute velocity error (r_ij * v_ij): %.10g\n", global_verr_max);
     }
-  }
+#endif
 }
 
 /* ---------------------------------------------------------------------- */
@@ -110,7 +108,10 @@ int FixRattle::setmask()
   mask |= POST_FORCE_RESPA;
   mask |= FINAL_INTEGRATE;
   mask |= FINAL_INTEGRATE_RESPA;
-  if (RATTLE_DEBUG) mask |= END_OF_STEP;
+  mask |= MIN_POST_FORCE;
+#if RATTLE_DEBUG
+  mask |= END_OF_STEP;
+#endif
   return mask;
 }
 
@@ -156,7 +157,7 @@ void FixRattle::post_force(int vflag)
 
   // communicate the unconstrained velocities
 
-  if (nprocs > 1) {
+  if (comm->nprocs > 1) {
     comm_mode = VP;
     comm->forward_comm(this);
   }
@@ -188,7 +189,7 @@ void FixRattle::post_force_respa(int vflag, int ilevel, int /*iloop*/)
 
   // communicate the unconstrained velocities
 
-  if (nprocs > 1) {
+  if (comm->nprocs > 1) {
     comm_mode = VP;
     comm->forward_comm(this);
   }
@@ -718,7 +719,7 @@ void FixRattle::unpack_forward_comm(int n, int first, double *buf)
 
 void FixRattle::shake_end_of_step(int vflag) {
 
-  if (nprocs > 1) {
+  if (comm->nprocs > 1) {
     comm_mode = V;
     comm->forward_comm(this);
   }
@@ -738,7 +739,6 @@ void FixRattle::correct_coordinates(int vflag) {
   FixShake::correct_coordinates(vflag);
 }
 
-
 /* ----------------------------------------------------------------------
    Remove the velocity component along any bond.
 ------------------------------------------------------------------------- */
@@ -759,7 +759,7 @@ void FixRattle::correct_velocities() {
 
   // communicate the unconstrained velocities
 
-  if (nprocs > 1) {
+  if (comm->nprocs > 1) {
     comm_mode = VP;
     comm->forward_comm(this);
   }
@@ -769,14 +769,13 @@ void FixRattle::correct_velocities() {
   int m;
   for (int i = 0; i < nlist; i++) {
     m = list[i];
-    if      (shake_flag[m] == 2)        vrattle2(m);
+    if      (shake_flag[m] == 2)  vrattle2(m);
-    else if (shake_flag[m] == 3)        vrattle3(m);
+    else if (shake_flag[m] == 3)  vrattle3(m);
-    else if (shake_flag[m] == 4)        vrattle4(m);
+    else if (shake_flag[m] == 4)  vrattle4(m);
-    else                                vrattle3angle(m);
+    else                          vrattle3angle(m);
   }
 }
 
-
 /* ----------------------------------------------------------------------
    DEBUGGING methods
    The functions below allow you to check whether the
@@ -788,16 +787,16 @@ void FixRattle::correct_velocities() {
 
 void FixRattle::end_of_step()
 {
-  if (nprocs > 1) {
+  if (comm->nprocs > 1) {
-       comm_mode = V;
+    comm_mode = V;
-       comm->forward_comm(this);
+    comm->forward_comm(this);
   }
-  if (!check_constraints(v, RATTLE_TEST_POS, RATTLE_TEST_VEL) && RATTLE_RAISE_ERROR) {
+#if RATTLE_RAISE_ERROR
+  if (!check_constraints(v, RATTLE_TEST_POS, RATTLE_TEST_VEL))
     error->one(FLERR, "Rattle failed ");
-  }
+#endif
 }
 
-
 /* ---------------------------------------------------------------------- */
 
 bool FixRattle::check_constraints(double **v, bool checkr, bool checkv)
@@ -807,12 +806,11 @@ bool FixRattle::check_constraints(double **v, bool checkr, bool checkv)
   int i=0;
   while (i < nlist && ret) {
     m = list[i];
-    if      (shake_flag[m] == 2)     ret =   check2(v,m,checkr,checkv);
+    if      (shake_flag[m] == 2)     ret = check2(v,m,checkr,checkv);
-    else if (shake_flag[m] == 3)     ret =   check3(v,m,checkr,checkv);
+    else if (shake_flag[m] == 3)     ret = check3(v,m,checkr,checkv);
-    else if (shake_flag[m] == 4)     ret =   check4(v,m,checkr,checkv);
+    else if (shake_flag[m] == 4)     ret = check4(v,m,checkr,checkv);
-    else                             ret =   check3angle(v,m,checkr,checkv);
+    else                             ret = check3angle(v,m,checkr,checkv);
     i++;
-    if (!RATTLE_RAISE_ERROR)         ret = true;
   }
   return ret;
 }
@@ -834,14 +832,10 @@ bool FixRattle::check2(double **v, int m, bool checkr, bool checkv)
   MathExtra::sub3(v[i1],v[i0],v01);
 
   stat = !(checkr && (fabs(sqrt(MathExtra::dot3(r01,r01)) - bond1) > tol));
-  if (!stat)
+  if (!stat) error->one(FLERR,"Coordinate constraints are not satisfied up to desired tolerance ");
-     error->one(FLERR,"Coordinate constraints are not satisfied "
-                "up to desired tolerance ");
 
   stat = !(checkv && (fabs(MathExtra::dot3(r01,v01)) > tol));
-  if (!stat)
+  if (!stat) error->one(FLERR,"Velocity constraints are not satisfied up to desired tolerance ");
-     error->one(FLERR,"Velocity constraints are not satisfied "
-                "up to desired tolerance ");
   return stat;
 }
 
@@ -871,15 +865,11 @@ bool FixRattle::check3(double **v, int m, bool checkr, bool checkv)
 
   stat = !(checkr && (fabs(sqrt(MathExtra::dot3(r01,r01)) - bond1) > tol ||
                       fabs(sqrt(MathExtra::dot3(r02,r02))-bond2) > tol));
-  if (!stat)
+  if (!stat) error->one(FLERR,"Coordinate constraints are not satisfied up to desired tolerance ");
-     error->one(FLERR,"Coordinate constraints are not satisfied "
-                "up to desired tolerance ");
 
   stat = !(checkv && (fabs(MathExtra::dot3(r01,v01)) > tol ||
                       fabs(MathExtra::dot3(r02,v02)) > tol));
-  if (!stat)
+  if (!stat) error->one(FLERR,"Velocity constraints are not satisfied up to desired tolerance ");
-     error->one(FLERR,"Velocity constraints are not satisfied "
-                "up to desired tolerance ");
   return stat;
 }
 
@@ -914,16 +904,12 @@ bool FixRattle::check4(double **v, int m, bool checkr, bool checkv)
   stat = !(checkr && (fabs(sqrt(MathExtra::dot3(r01,r01)) - bond1) > tol ||
                       fabs(sqrt(MathExtra::dot3(r02,r02))-bond2) > tol ||
                       fabs(sqrt(MathExtra::dot3(r03,r03))-bond3) > tol));
-  if (!stat)
+  if (!stat) error->one(FLERR,"Coordinate constraints are not satisfied up to desired tolerance ");
-     error->one(FLERR,"Coordinate constraints are not satisfied "
-                "up to desired tolerance ");
 
   stat = !(checkv && (fabs(MathExtra::dot3(r01,v01)) > tol ||
                       fabs(MathExtra::dot3(r02,v02)) > tol ||
                       fabs(MathExtra::dot3(r03,v03)) > tol));
-  if (!stat)
+  if (!stat) error->one(FLERR,"Velocity constraints are not satisfied up to desired tolerance ");
-     error->one(FLERR,"Velocity constraints are not satisfied "
-                "up to desired tolerance ");
   return stat;
 }
 
@@ -954,25 +940,19 @@ bool FixRattle::check3angle(double **v, int m, bool checkr, bool checkv)
   MathExtra::sub3(v[i2],v[i0],v02);
   MathExtra::sub3(v[i2],v[i1],v12);
 
-
-
   double db1 = fabs(sqrt(MathExtra::dot3(r01,r01)) - bond1);
   double db2 = fabs(sqrt(MathExtra::dot3(r02,r02))-bond2);
   double db12 = fabs(sqrt(MathExtra::dot3(r12,r12))-bond12);
 
-
+  stat = !(checkr && (db1 > tol || db2 > tol || db12 > tol));
-  stat = !(checkr && (db1 > tol ||
-                      db2 > tol ||
-                      db12 > tol));
 
   if (derr_max < db1/bond1)    derr_max = db1/bond1;
   if (derr_max < db2/bond2)    derr_max = db2/bond2;
   if (derr_max < db12/bond12)  derr_max = db12/bond12;
 
-
+#if RATTLE_RAISE_ERROR
-  if (!stat && RATTLE_RAISE_ERROR)
+  if (!stat) error->one(FLERR,"Coordinate constraints are not satisfied up to desired tolerance ");
-     error->one(FLERR,"Coordinate constraints are not satisfied "
+#endif
-                "up to desired tolerance ");
 
   double dv1 = fabs(MathExtra::dot3(r01,v01));
   double dv2 = fabs(MathExtra::dot3(r02,v02));
@@ -982,16 +962,10 @@ bool FixRattle::check3angle(double **v, int m, bool checkr, bool checkv)
   if (verr_max < dv2)    verr_max = dv2;
   if (verr_max < dv12)   verr_max = dv12;
 
+  stat = !(checkv && (dv1 > tol || dv2 > tol || dv12> tol));
 
-  stat = !(checkv && (dv1 > tol ||
+#if RATTLE_RAISE_ERROR
-                      dv2 > tol ||
+  if (!stat) error->one(FLERR,"Velocity constraints are not satisfied up to desired tolerance!");
-                      dv12> tol));
+#endif
-
-
-  if (!stat && RATTLE_RAISE_ERROR)
-     error->one(FLERR,"Velocity constraints are not satisfied "
-                "up to desired tolerance!");
-
-
   return stat;
 }

src/RIGID/fix_shake.cpp

@@ -41,8 +41,8 @@ using namespace MathConst;
 
 #define RVOUS 1   // 0 for irregular, 1 for all2all
 
-#define BIG 1.0e20
+static constexpr double BIG = 1.0e20;
-#define MASSDELTA 0.1
+static constexpr double MASSDELTA = 0.1;
 
 /* ---------------------------------------------------------------------- */
 
@@ -57,21 +57,20 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
   a_count_all(nullptr), a_ave(nullptr), a_max(nullptr), a_min(nullptr), a_ave_all(nullptr),
   a_max_all(nullptr), a_min_all(nullptr), atommols(nullptr), onemols(nullptr)
 {
-  MPI_Comm_rank(world,&me);
+  energy_global_flag = energy_peratom_flag = 1;
-  MPI_Comm_size(world,&nprocs);
-
   virial_global_flag = virial_peratom_flag = 1;
-  thermo_virial = 1;
+  thermo_energy = thermo_virial = 1;
   create_attribute = 1;
   dof_flag = 1;
   stores_ids = 1;
   centroidstressflag = CENTROID_AVAIL;
+  next_output = -1;
 
   // error check
 
   molecular = atom->molecular;
   if (molecular == Atom::ATOMIC)
-    error->all(FLERR,"Cannot use fix shake with non-molecular system");
+    error->all(FLERR,"Cannot use fix {} with non-molecular system", style);
 
   // perform initial allocation of atom-based arrays
   // register with Atom class
@@ -92,8 +91,9 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
   comm_forward = 3;
 
   // parse SHAKE args
+  auto mystyle = fmt::format("fix {}",style);
 
-  if (narg < 8) error->all(FLERR,"Illegal fix shake command");
+  if (narg < 8) utils::missing_cmd_args(FLERR,mystyle, error);
 
   tolerance = utils::numeric(FLERR,arg[3],false,lmp);
   max_iter = utils::inumeric(FLERR,arg[4],false,lmp);
@@ -133,49 +133,55 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
     else if (mode == 'b') {
       int i = utils::inumeric(FLERR,arg[next],false,lmp);
       if (i < 1 || i > atom->nbondtypes)
-        error->all(FLERR,"Invalid bond type index for fix shake");
+        error->all(FLERR,"Invalid bond type index for {}", mystyle);
       bond_flag[i] = 1;
 
     } else if (mode == 'a') {
       int i = utils::inumeric(FLERR,arg[next],false,lmp);
       if (i < 1 || i > atom->nangletypes)
-        error->all(FLERR,"Invalid angle type index for fix shake");
+        error->all(FLERR,"Invalid angle type index for {}", mystyle);
       angle_flag[i] = 1;
 
     } else if (mode == 't') {
       int i = utils::inumeric(FLERR,arg[next],false,lmp);
       if (i < 1 || i > atom->ntypes)
-        error->all(FLERR,"Invalid atom type index for fix shake");
+        error->all(FLERR,"Invalid atom type index for {}", mystyle);
       type_flag[i] = 1;
 
     } else if (mode == 'm') {
       double massone = utils::numeric(FLERR,arg[next],false,lmp);
-      if (massone == 0.0) error->all(FLERR,"Invalid atom mass for fix shake");
+      if (massone == 0.0) error->all(FLERR,"Invalid atom mass for {}", mystyle);
       if (nmass == atom->ntypes)
-        error->all(FLERR,"Too many masses for fix shake");
+        error->all(FLERR,"Too many masses for {}", mystyle);
       mass_list[nmass++] = massone;
 
-    } else error->all(FLERR,"Illegal fix shake command");
+    } else error->all(FLERR,"Unknown {} command option: {}", mystyle, arg[next]);
     next++;
   }
 
   // parse optional args
 
   onemols = nullptr;
+  kbond = 1.0e6;
 
   int iarg = next;
   while (iarg < narg) {
-    if (strcmp(arg[next],"mol") == 0) {
+    if (strcmp(arg[iarg],"mol") == 0) {
-      if (iarg+2 > narg) error->all(FLERR,"Illegal fix shake command");
+      if (iarg+2 > narg) utils::missing_cmd_args(FLERR,mystyle+" mol",error);
       int imol = atom->find_molecule(arg[iarg+1]);
       if (imol == -1)
-        error->all(FLERR,"Molecule template ID for fix shake does not exist");
+        error->all(FLERR,"Molecule template ID {} for {} does not exist", mystyle, arg[iarg+1]);
-      if (atom->molecules[imol]->nset > 1 && comm->me == 0)
+      if ((atom->molecules[imol]->nset > 1) && (comm->me == 0))
-        error->warning(FLERR,"Molecule template for fix shake has multiple molecules");
+        error->warning(FLERR,"Molecule template for {} has multiple molecules", mystyle);
       onemols = &atom->molecules[imol];
       nmol = onemols[0]->nset;
       iarg += 2;
-    } else error->all(FLERR,"Illegal fix shake command");
+    } else if (strcmp(arg[iarg],"kbond") == 0) {
+      if (iarg+2 > narg) utils::missing_cmd_args(FLERR,mystyle+" kbond",error);
+      kbond = utils::numeric(FLERR, arg[iarg+1], false, lmp);
+      if (kbond < 0) error->all(FLERR,"Illegal {} kbond value {}. Must be >= 0.0", mystyle, kbond);
+      iarg += 2;
+    } else error->all(FLERR,"Unknown {} command option: {}", mystyle, arg[iarg]);
   }
 
   // error check for Molecule template
@@ -183,7 +189,7 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
   if (onemols) {
     for (int i = 0; i < nmol; i++)
       if (onemols[i]->shakeflag == 0)
-        error->all(FLERR,"Fix shake molecule template must have shake info");
+        error->all(FLERR,"Fix {} molecule template must have shake info", style);
   }
 
   // allocate bond and angle distance arrays, indexed from 1 to n
@@ -321,6 +327,7 @@ int FixShake::setmask()
   mask |= PRE_NEIGHBOR;
   mask |= POST_FORCE;
   mask |= POST_FORCE_RESPA;
+  mask |= MIN_POST_FORCE;
   return mask;
 }
 
@@ -335,28 +342,24 @@ void FixShake::init()
   double rsq,angle;
 
   // error if more than one shake fix
+  auto pattern = fmt::format("^{}",style);
 
-  int count = 0;
+  if (modify->get_fix_by_style(pattern).size() > 1)
-  for (i = 0; i < modify->nfix; i++)
+    error->all(FLERR,"More than one fix {} instance",style);
-    if (strcmp(modify->fix[i]->style,"shake") == 0) count++;
-  if (count > 1) error->all(FLERR,"More than one fix shake");
 
   // cannot use with minimization since SHAKE turns off bonds
   // that should contribute to potential energy
 
-  if (update->whichflag == 2)
+  if ((comm->me == 0) && (update->whichflag == 2))
-    error->all(FLERR,"Fix shake cannot be used with minimization");
+    error->warning(FLERR,"Using fix {} with minimization. Substituting constraints with "
+                   "restraint forces using k={:.4g}", style, kbond);
 
-  // error if npt,nph fix comes before shake fix
+  // error if a fix changing the box comes before shake fix
-
+  bool boxflag = false;
-  for (i = 0; i < modify->nfix; i++) {
+  for (auto ifix : modify->get_fix_list()) {
-    if (strcmp(modify->fix[i]->style,"npt") == 0) break;
+   if (boxflag && utils::strmatch(ifix->style,pattern))
-    if (strcmp(modify->fix[i]->style,"nph") == 0) break;
+     error->all(FLERR,"Fix {} must come before any box changing fix", style);
-  }
+    if (ifix->box_change) boxflag = true;
-  if (i < modify->nfix) {
-    for (int j = i; j < modify->nfix; j++)
-      if (strcmp(modify->fix[j]->style,"shake") == 0)
-        error->all(FLERR,"Shake fix must come before NPT/NPH fix");
   }
 
   // if rRESPA, find associated fix that must exist
@@ -379,7 +382,7 @@ void FixShake::init()
   // set equilibrium bond distances
 
   if (force->bond == nullptr)
-    error->all(FLERR,"Bond potential must be defined for SHAKE");
+    error->all(FLERR,"Bond style must be defined for fix {}",style);
   for (i = 1; i <= atom->nbondtypes; i++)
     bond_distance[i] = force->bond->equilibrium_distance(i);
 
@@ -390,7 +393,7 @@ void FixShake::init()
   for (i = 1; i <= atom->nangletypes; i++) {
     if (angle_flag[i] == 0) continue;
     if (force->angle == nullptr)
-      error->all(FLERR,"Angle potential must be defined for SHAKE");
+      error->all(FLERR,"Angle style must be defined for fix {}",style);
 
     // scan all atoms for a SHAKE angle cluster
     // extract bond types for the 2 bonds in the cluster
@@ -417,7 +420,7 @@ void FixShake::init()
     // error check for any bond types that are not the same
 
     MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_MAX,world);
-    if (flag_all) error->all(FLERR,"Shake angles have different bond types");
+    if (flag_all) error->all(FLERR,"Fix {} angles have different bond types", style);
 
     // insure all procs have bond types
 
@@ -494,6 +497,16 @@ void FixShake::setup(int vflag)
   shake_end_of_step(vflag);
 }
 
+/* ----------------------------------------------------------------------
+   during minimization fix SHAKE adds strong bond forces
+------------------------------------------------------------------------- */
+
+void FixShake::min_setup(int vflag)
+{
+  pre_neighbor();
+  min_post_force(vflag);
+}
+
 /* ----------------------------------------------------------------------
    build list of SHAKE clusters to constrain
    if one or more atoms in cluster are on this proc,
@@ -533,19 +546,17 @@ void FixShake::pre_neighbor()
         atom1 = atom->map(shake_atom[i][0]);
         atom2 = atom->map(shake_atom[i][1]);
         if (atom1 == -1 || atom2 == -1)
-          error->one(FLERR,"Shake atoms {} {} missing on proc "
+          error->one(FLERR,"Shake atoms {} {} missing on proc {} at step {}",shake_atom[i][0],
-                                       "{} at step {}",shake_atom[i][0],
+                     shake_atom[i][1],comm->me,update->ntimestep);
-                                       shake_atom[i][1],me,update->ntimestep);
         if (i <= atom1 && i <= atom2) list[nlist++] = i;
       } else if (shake_flag[i] % 2 == 1) {
         atom1 = atom->map(shake_atom[i][0]);
         atom2 = atom->map(shake_atom[i][1]);
         atom3 = atom->map(shake_atom[i][2]);
         if (atom1 == -1 || atom2 == -1 || atom3 == -1)
-          error->one(FLERR,"Shake atoms {} {} {} missing on proc "
+          error->one(FLERR,"Shake atoms {} {} {} missing on proc {} at step {}",shake_atom[i][0],
-                                       "{} at step {}",shake_atom[i][0],
                                        shake_atom[i][1],shake_atom[i][2],
-                                       me,update->ntimestep);
+                                       comm->me,update->ntimestep);
         if (i <= atom1 && i <= atom2 && i <= atom3) list[nlist++] = i;
       } else {
         atom1 = atom->map(shake_atom[i][0]);
@@ -553,10 +564,9 @@ void FixShake::pre_neighbor()
         atom3 = atom->map(shake_atom[i][2]);
         atom4 = atom->map(shake_atom[i][3]);
         if (atom1 == -1 || atom2 == -1 || atom3 == -1 || atom4 == -1)
-          error->one(FLERR,"Shake atoms {} {} {} {} missing on "
+          error->one(FLERR,"Shake atoms {} {} {} {} missing on proc {} at step {}",shake_atom[i][0],
-                                       "proc {} at step {}",shake_atom[i][0],
                                        shake_atom[i][1],shake_atom[i][2],
-                                       shake_atom[i][3],me,update->ntimestep);
+                                       shake_atom[i][3],comm->me,update->ntimestep);
         if (i <= atom1 && i <= atom2 && i <= atom3 && i <= atom4)
           list[nlist++] = i;
       }
@@ -575,7 +585,7 @@ void FixShake::post_force(int vflag)
   // communicate results if necessary
 
   unconstrained_update();
-  if (nprocs > 1) comm->forward_comm(this);
+  if (comm->nprocs > 1) comm->forward_comm(this);
 
   // virial setup
 
@@ -619,7 +629,7 @@ void FixShake::post_force_respa(int vflag, int ilevel, int iloop)
   // communicate results if necessary
 
   unconstrained_update_respa(ilevel);
-  if (nprocs > 1) comm->forward_comm(this);
+  if (comm->nprocs > 1) comm->forward_comm(this);
 
   // virial setup only needed on last iteration of innermost level
   //   and if pressure is requested
@@ -644,6 +654,48 @@ void FixShake::post_force_respa(int vflag, int ilevel, int iloop)
   vflag_post_force = vflag;
 }
 
+/* ----------------------------------------------------------------------
+   substitute shake constraints with very strong bonds
+------------------------------------------------------------------------- */
+
+void FixShake::min_post_force(int vflag)
+{
+  if (output_every) {
+    bigint ntimestep = update->ntimestep;
+    if (next_output == ntimestep) stats();
+
+    next_output = ntimestep + output_every;
+    if (ntimestep % output_every != 0)
+      next_output = (ntimestep/output_every)*output_every + output_every;
+  } else next_output = -1;
+
+  v_init(vflag);
+
+  x = atom->x;
+  f = atom->f;
+  nlocal = atom->nlocal;
+
+  // loop over clusters to add strong restraint forces
+
+  for (int i = 0; i < nlist; i++) {
+    int m = list[i];
+    if (shake_flag[m] == 2) {
+      bond_force(shake_atom[m][0], shake_atom[m][1], bond_distance[shake_type[m][0]]);
+    } else if (shake_flag[m] == 3) {
+      bond_force(shake_atom[m][0], shake_atom[m][1], bond_distance[shake_type[m][0]]);
+      bond_force(shake_atom[m][0], shake_atom[m][2], bond_distance[shake_type[m][1]]);
+    } else if (shake_flag[m] == 4) {
+      bond_force(shake_atom[m][0], shake_atom[m][1], bond_distance[shake_type[m][0]]);
+      bond_force(shake_atom[m][0], shake_atom[m][2], bond_distance[shake_type[m][1]]);
+      bond_force(shake_atom[m][0], shake_atom[m][3], bond_distance[shake_type[m][2]]);
+    } else {
+      bond_force(shake_atom[m][0], shake_atom[m][1], bond_distance[shake_type[m][0]]);
+      bond_force(shake_atom[m][0], shake_atom[m][2], bond_distance[shake_type[m][1]]);
+      bond_force(shake_atom[m][1], shake_atom[m][2], angle_distance[shake_type[m][2]]);
+    }
+  }
+}
+
 /* ----------------------------------------------------------------------
    count # of degrees-of-freedom removed by SHAKE for atoms in igroup
 ------------------------------------------------------------------------- */
@@ -690,10 +742,7 @@ void FixShake::find_clusters()
   tagint tagprev;
   double massone;
 
-  if ((me == 0) && screen) {
+  if (comm->me == 0) utils::logmesg(lmp, "Finding {} clusters ...\n",utils::uppercase(style));
-    if (!rattle) fputs("Finding SHAKE clusters ...\n",screen);
-    else fputs("Finding RATTLE clusters ...\n",screen);
-  }
   atommols = atom->avec->onemols;
 
   tagint *tag = atom->tag;
@@ -805,7 +854,7 @@ void FixShake::find_clusters()
     }
 
   MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_SUM,world);
-  if (flag_all) error->all(FLERR,"Did not find fix shake partner info");
+  if (flag_all) error->all(FLERR,"Did not find fix {} partner info", style);
 
   // -----------------------------------------------------
   // identify SHAKEable bonds
@@ -1015,7 +1064,7 @@ void FixShake::find_clusters()
   tmp = count4;
   MPI_Allreduce(&tmp,&count4,1,MPI_INT,MPI_SUM,world);
 
-  if (me == 0) {
+  if (comm->me == 0) {
     utils::logmesg(lmp,"{:>8} = # of size 2 clusters\n"
                    "{:>8} = # of size 3 clusters\n"
                    "{:>8} = # of size 4 clusters\n"
@@ -1043,8 +1092,8 @@ void FixShake::atom_owners()
   // one datum for each owned atom: datum = owning proc, atomID
 
   for (int i = 0; i < nlocal; i++) {
-    proclist[i] = tag[i] % nprocs;
+    proclist[i] = tag[i] % comm->nprocs;
-    idbuf[i].me = me;
+    idbuf[i].me = comm->me;
     idbuf[i].atomID = tag[i];
   }
 
@@ -1089,7 +1138,7 @@ void FixShake::partner_info(int *npartner, tagint **partner_tag,
   // set values in 4 partner arrays for all partner atoms I own
   // also setup input buf to rendezvous comm
   // input datums = pair of bonded atoms where I do not own partner
-  // owning proc for each datum = partner_tag % nprocs
+  // owning proc for each datum = partner_tag % comm->nprocs
   // datum: atomID = partner_tag (off-proc), partnerID = tag (on-proc)
   //        4 values for my owned atom
 
@@ -1127,7 +1176,7 @@ void FixShake::partner_info(int *npartner, tagint **partner_tag,
         }
 
       } else {
-        proclist[nsend] = partner_tag[i][j] % nprocs;
+        proclist[nsend] = partner_tag[i][j] % comm->nprocs;
         inbuf[nsend].atomID = partner_tag[i][j];
         inbuf[nsend].partnerID = tag[i];
         inbuf[nsend].mask = mask[i];
@@ -1217,7 +1266,7 @@ void FixShake::nshake_info(int *npartner, tagint **partner_tag,
   // set partner_nshake for all partner atoms I own
   // also setup input buf to rendezvous comm
   // input datums = pair of bonded atoms where I do not own partner
-  // owning proc for each datum = partner_tag % nprocs
+  // owning proc for each datum = partner_tag % comm->nprocs
   // datum: atomID = partner_tag (off-proc), partnerID = tag (on-proc)
   //        nshake value for my owned atom
 
@@ -1231,7 +1280,7 @@ void FixShake::nshake_info(int *npartner, tagint **partner_tag,
       if (m >= 0 && m < nlocal) {
         partner_nshake[i][j] = nshake[m];
       } else {
-        proclist[nsend] = partner_tag[i][j] % nprocs;
+        proclist[nsend] = partner_tag[i][j] % comm->nprocs;
         inbuf[nsend].atomID = partner_tag[i][j];
         inbuf[nsend].partnerID = tag[i];
         inbuf[nsend].nshake = nshake[i];
@@ -1295,7 +1344,7 @@ void FixShake::shake_info(int *npartner, tagint **partner_tag,
   // set 3 shake arrays for all partner atoms I own
   // also setup input buf to rendezvous comm
   // input datums = partner atom where I do not own partner
-  // owning proc for each datum = partner_tag % nprocs
+  // owning proc for each datum = partner_tag % comm->nprocs
   // datum: atomID = partner_tag (off-proc)
   //        values in 3 shake arrays
 
@@ -1317,7 +1366,7 @@ void FixShake::shake_info(int *npartner, tagint **partner_tag,
         shake_type[m][2] = shake_type[i][2];
 
       } else {
-        proclist[nsend] = partner_tag[i][j] % nprocs;
+        proclist[nsend] = partner_tag[i][j] % comm->nprocs;
         inbuf[nsend].atomID = partner_tag[i][j];
         inbuf[nsend].shake_flag = shake_flag[i];
         inbuf[nsend].shake_atom[0] = shake_atom[i][0];
@@ -2451,6 +2500,53 @@ void FixShake::shake3angle(int m)
   }
 }
 
+/* ----------------------------------------------------------------------
+   apply bond force for minimization
+------------------------------------------------------------------------- */
+
+void FixShake::bond_force(tagint id1, tagint id2, double length)
+{
+
+  int i1 = atom->map(id1);
+  int i2 = atom->map(id2);
+
+  if ((i1 < 0) || (i2 < 0)) return;
+
+  // distance vec between atoms, with PBC
+
+  double delx = x[i1][0] - x[i2][0];
+  double dely = x[i1][1] - x[i2][1];
+  double delz = x[i1][2] - x[i2][2];
+  domain->minimum_image(delx, dely, delz);
+
+  // compute and apply force
+
+  const double r = sqrt(delx * delx + dely * dely + delz * delz);
+  const double dr = r - length;
+  const double rk = kbond * dr;
+  const double fbond = (r > 0.0) ? -2.0 * rk / r : 0.0;
+  double v[6];
+  v[0] = 0.5 * delx * delx * fbond;
+  v[1] = 0.5 * dely * dely * fbond;
+  v[2] = 0.5 * delz * delz * fbond;
+  v[3] = 0.5 * delx * dely * fbond;
+  v[4] = 0.5 * delx * delz * fbond;
+  v[5] = 0.5 * dely * delz * fbond;
+
+  if (i1 < nlocal) {
+    f[i1][0] += delx * fbond;
+    f[i1][1] += dely * fbond;
+    f[i1][2] += delz * fbond;
+    if (evflag) v_tally(i1, v);
+  }
+  if (i2 < nlocal) {
+    f[i2][0] -= delx * fbond;
+    f[i2][1] -= dely * fbond;
+    f[i2][2] -= delz * fbond;
+    if (evflag) v_tally(i2, v);
+  }
+}
+
 /* ----------------------------------------------------------------------
    print-out bond & angle statistics
 ------------------------------------------------------------------------- */
@@ -2558,9 +2654,10 @@ void FixShake::stats()
 
   // print stats only for non-zero counts
 
-  if (me == 0) {
+  if (comm->me == 0) {
     const int width = log10((MAX(MAX(1,nb),na)))+2;
-    auto mesg = fmt::format("SHAKE stats (type/ave/delta/count) on step {}\n", update->ntimestep);
+    auto mesg = fmt::format("{} stats (type/ave/delta/count) on step {}\n",
+                            utils::uppercase(style), update->ntimestep);
     for (i = 1; i < nb; i++) {
       const auto bcnt = b_count_all[i];
       if (bcnt)
@@ -3109,7 +3206,7 @@ void FixShake::correct_coordinates(int vflag) {
 
   double **xtmp = xshake;
   xshake = x;
-  if (nprocs > 1) {
+  if (comm->nprocs > 1) {
     comm->forward_comm(this);
   }
   xshake = xtmp;

src/RIGID/fix_shake.h

@@ -34,9 +34,11 @@ class FixShake : public Fix {
   int setmask() override;
   void init() override;
   void setup(int) override;
+  void min_setup(int) override;
   void pre_neighbor() override;
   void post_force(int) override;
   void post_force_respa(int, int, int) override;
+  void min_post_force(int) override;
 
   double memory_usage() override;
   void grow_arrays(int) override;
@@ -61,12 +63,11 @@ class FixShake : public Fix {
  protected:
   int vflag_post_force;    // store the vflag of last post_force call
   int respa;               // 0 = vel. Verlet, 1 = respa
-  int me, nprocs;
+  int rattle;              // 0 = SHAKE, 1 = RATTLE
-  int rattle;            // 0 = SHAKE, 1 = RATTLE
+  double tolerance;        // SHAKE tolerance
-  double tolerance;      // SHAKE tolerance
+  int max_iter;            // max # of SHAKE iterations
-  int max_iter;          // max # of SHAKE iterations
+  int output_every;        // SHAKE stat output every so often
-  int output_every;      // SHAKE stat output every so often
+  bigint next_output;      // timestep for next output
-  bigint next_output;    // timestep for next output
 
   // settings from input command
   int *bond_flag, *angle_flag;    // bond/angle types to constrain
@@ -76,6 +77,7 @@ class FixShake : public Fix {
 
   int molecular;                             // copy of atom->molecular
   double *bond_distance, *angle_distance;    // constraint distances
+  double kbond;                              // force constant for restraint
 
   class FixRespa *fix_respa;    // rRESPA fix needed by SHAKE
   int nlevels_respa;            // copies of needed rRESPA variables
@@ -133,6 +135,7 @@ class FixShake : public Fix {
   void shake3(int);
   void shake4(int);
   void shake3angle(int);
+  void bond_force(tagint, tagint, double);
   void stats();
   int bondtype_findset(int, tagint, tagint, int);
   int angletype_findset(int, tagint, tagint, int);

src/fix_restrain.cpp

@@ -19,17 +19,18 @@
 
 #include "fix_restrain.h"
 
-#include <cmath>
-#include <cstring>
 #include "atom.h"
-#include "force.h"
-#include "update.h"
-#include "domain.h"
 #include "comm.h"
-#include "respa.h"
+#include "domain.h"
+#include "error.h"
+#include "force.h"
 #include "math_const.h"
 #include "memory.h"
-#include "error.h"
+#include "respa.h"
+#include "update.h"
+
+#include <cmath>
+#include <cstring>
 
 using namespace LAMMPS_NS;
 using namespace FixConst;
@@ -271,14 +272,12 @@ void FixRestrain::restrain_bond(int m)
   if (newton_bond) {
     if (i2 == -1 || i2 >= nlocal) return;
     if (i1 == -1)
-      error->one(FLERR,"Restrain atoms {} {} missing on "
+      error->one(FLERR,"Restrain atoms {} {} missing on proc {} at step {}", ids[m][0],ids[m][1],
-                                   "proc {} at step {}", ids[m][0],ids[m][1],
                                    comm->me,update->ntimestep);
   } else {
     if ((i1 == -1 || i1 >= nlocal) && (i2 == -1 || i2 >= nlocal)) return;
     if (i1 == -1 || i2 == -1)
-      error->one(FLERR,"Restrain atoms {} {} missing on "
+      error->one(FLERR,"Restrain atoms {} {} missing on proc {} at step {}", ids[m][0],ids[m][1],
-                                   "proc {} at step {}", ids[m][0],ids[m][1],
                                    comm->me,update->ntimestep);
   }