diff --git a/doc/src/fix_shake.rst b/doc/src/fix_shake.rst index f0c847cb5e..d723f28fc0 100644 --- a/doc/src/fix_shake.rst +++ b/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 ` 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 ` or :doc:`fix pour `, add molecules +The *mol* keyword should be used when other commands, such as :doc:`fix +deposit ` or :doc:`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 ` 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 -file. See the :doc:`molecule ` command for details. The only +bond restrictions, and SHAKE info can be specified in the molecule file. +See the :doc:`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 `. No global or per-atom quantities are stored by these fixes for access by various :doc:`output commands `. No parameter of these fixes can be used with the *start/stop* keywords of the -:doc:`run ` command. These fixes are not invoked during -:doc:`energy minimization `. +:doc:`run ` command. + +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 ehex `, +`fix nve/manifold/rattle ` Default """"""" -none +kbond = 1.0e6 ---------- diff --git a/src/KOKKOS/fix_shake_kokkos.cpp b/src/KOKKOS/fix_shake_kokkos.cpp index a24c47c1e2..5ed8cfccb3 100644 --- a/src/KOKKOS/fix_shake_kokkos.cpp +++ b/src/KOKKOS/fix_shake_kokkos.cpp @@ -12,12 +12,8 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#include -#include -#include -#include -#include #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 +#include + using namespace LAMMPS_NS; using namespace FixConst; using namespace MathConst; @@ -292,7 +291,7 @@ void FixShakeKokkos::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::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(); // virial setup @@ -1702,7 +1701,7 @@ void FixShakeKokkos::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; diff --git a/src/RIGID/fix_rattle.cpp b/src/RIGID/fix_rattle.cpp index afe08415c3..1d6336712b 100644 --- a/src/RIGID/fix_rattle.cpp +++ b/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 (npid == 0 && screen) { + if (comm->me == 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); - else if (shake_flag[m] == 3) vrattle3(m); - else if (shake_flag[m] == 4) vrattle4(m); - else vrattle3angle(m); + if (shake_flag[m] == 2) vrattle2(m); + else if (shake_flag[m] == 3) vrattle3(m); + else if (shake_flag[m] == 4) vrattle4(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) { - comm_mode = V; - comm->forward_comm(this); + if (comm->nprocs > 1) { + comm_mode = V; + 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); - 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 ret = check3angle(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] == 4) ret = check4(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) - error->one(FLERR,"Coordinate constraints are not satisfied " - "up to desired tolerance "); + if (!stat) error->one(FLERR,"Coordinate constraints are not satisfied up to desired tolerance "); stat = !(checkv && (fabs(MathExtra::dot3(r01,v01)) > tol)); - if (!stat) - error->one(FLERR,"Velocity constraints are not satisfied " - "up to desired tolerance "); + if (!stat) 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) - error->one(FLERR,"Coordinate constraints are not satisfied " - "up to desired tolerance "); + if (!stat) 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) - error->one(FLERR,"Velocity constraints are not satisfied " - "up to desired tolerance "); + if (!stat) 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) - error->one(FLERR,"Coordinate constraints are not satisfied " - "up to desired tolerance "); + if (!stat) 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) - error->one(FLERR,"Velocity constraints are not satisfied " - "up to desired tolerance "); + if (!stat) 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 (!stat && RATTLE_RAISE_ERROR) - error->one(FLERR,"Coordinate constraints are not satisfied " - "up to desired tolerance "); +#if RATTLE_RAISE_ERROR + if (!stat) error->one(FLERR,"Coordinate constraints are not satisfied up to desired tolerance "); +#endif 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 || - dv2 > tol || - dv12> tol)); - - - if (!stat && RATTLE_RAISE_ERROR) - error->one(FLERR,"Velocity constraints are not satisfied " - "up to desired tolerance!"); - - +#if RATTLE_RAISE_ERROR + if (!stat) error->one(FLERR,"Velocity constraints are not satisfied up to desired tolerance!"); +#endif return stat; } diff --git a/src/RIGID/fix_shake.cpp b/src/RIGID/fix_shake.cpp index d74f72fb69..41d8c1599c 100644 --- a/src/RIGID/fix_shake.cpp +++ b/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 -#define MASSDELTA 0.1 +static constexpr double BIG = 1.0e20; +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); - MPI_Comm_size(world,&nprocs); - + energy_global_flag = energy_peratom_flag = 1; 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 (iarg+2 > narg) error->all(FLERR,"Illegal fix shake command"); + if (strcmp(arg[iarg],"mol") == 0) { + 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"); - if (atom->molecules[imol]->nset > 1 && comm->me == 0) - error->warning(FLERR,"Molecule template for fix shake has multiple molecules"); + error->all(FLERR,"Molecule template ID {} for {} does not exist", mystyle, arg[iarg+1]); + if ((atom->molecules[imol]->nset > 1) && (comm->me == 0)) + 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; - for (i = 0; i < modify->nfix; i++) - if (strcmp(modify->fix[i]->style,"shake") == 0) count++; - if (count > 1) error->all(FLERR,"More than one fix shake"); + if (modify->get_fix_by_style(pattern).size() > 1) + error->all(FLERR,"More than one fix {} instance",style); // cannot use with minimization since SHAKE turns off bonds // that should contribute to potential energy - if (update->whichflag == 2) - error->all(FLERR,"Fix shake cannot be used with minimization"); + if ((comm->me == 0) && (update->whichflag == 2)) + 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 - - for (i = 0; i < modify->nfix; i++) { - if (strcmp(modify->fix[i]->style,"npt") == 0) break; - if (strcmp(modify->fix[i]->style,"nph") == 0) break; - } - 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"); + // error if a fix changing the box comes before shake fix + bool boxflag = false; + for (auto ifix : modify->get_fix_list()) { + if (boxflag && utils::strmatch(ifix->style,pattern)) + error->all(FLERR,"Fix {} must come before any box changing fix", style); + if (ifix->box_change) boxflag = true; } // 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 " - "{} at step {}",shake_atom[i][0], - shake_atom[i][1],me,update->ntimestep); + error->one(FLERR,"Shake atoms {} {} missing on proc {} at step {}",shake_atom[i][0], + shake_atom[i][1],comm->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 " - "{} at step {}",shake_atom[i][0], + error->one(FLERR,"Shake atoms {} {} {} missing on proc {} 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 " - "proc {} at step {}",shake_atom[i][0], + error->one(FLERR,"Shake atoms {} {} {} {} missing on 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 (!rattle) fputs("Finding SHAKE clusters ...\n",screen); - else fputs("Finding RATTLE clusters ...\n",screen); - } + if (comm->me == 0) utils::logmesg(lmp, "Finding {} clusters ...\n",utils::uppercase(style)); 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; - idbuf[i].me = me; + proclist[i] = tag[i] % comm->nprocs; + 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; diff --git a/src/RIGID/fix_shake.h b/src/RIGID/fix_shake.h index 12e24fb350..a0f743d7d2 100644 --- a/src/RIGID/fix_shake.h +++ b/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 - double tolerance; // SHAKE tolerance - int max_iter; // max # of SHAKE iterations - int output_every; // SHAKE stat output every so often - bigint next_output; // timestep for next output + int rattle; // 0 = SHAKE, 1 = RATTLE + double tolerance; // SHAKE tolerance + int max_iter; // max # of SHAKE iterations + int output_every; // SHAKE stat output every so often + 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); diff --git a/src/fix_restrain.cpp b/src/fix_restrain.cpp index 8b97715ff6..e5440830f5 100644 --- a/src/fix_restrain.cpp +++ b/src/fix_restrain.cpp @@ -19,17 +19,18 @@ #include "fix_restrain.h" -#include -#include #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 +#include 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 " - "proc {} at step {}", ids[m][0],ids[m][1], + error->one(FLERR,"Restrain atoms {} {} missing on 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 " - "proc {} at step {}", ids[m][0],ids[m][1], + error->one(FLERR,"Restrain atoms {} {} missing on proc {} at step {}", ids[m][0],ids[m][1], comm->me,update->ntimestep); }