From a7b6dc7b5954077cd3da19a3c77798eca683739b Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 30 Apr 2022 19:03:28 -0400 Subject: [PATCH 01/10] initial implementation of minimizer support in fix shake/rattle --- doc/src/fix_shake.rst | 27 ++-- src/KOKKOS/fix_shake_kokkos.cpp | 15 +- src/RIGID/fix_rattle.cpp | 104 ++++++-------- src/RIGID/fix_shake.cpp | 233 ++++++++++++++++++++++---------- src/RIGID/fix_shake.h | 15 +- src/fix_restrain.cpp | 21 ++- 6 files changed, 249 insertions(+), 166 deletions(-) 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); } From 322bf1ef477e928a49796eacdfe686e897a3a4c1 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 9 Jun 2022 22:34:35 -0400 Subject: [PATCH 02/10] compute energy due to restraint forces during minimization. output stats. --- doc/src/fix_shake.rst | 31 ++++++++++++++--------- src/RIGID/fix_shake.cpp | 54 ++++++++++++++++++++++++++--------------- src/RIGID/fix_shake.h | 3 +++ 3 files changed, 57 insertions(+), 31 deletions(-) diff --git a/doc/src/fix_shake.rst b/doc/src/fix_shake.rst index d723f28fc0..137cb1aedb 100644 --- a/doc/src/fix_shake.rst +++ b/doc/src/fix_shake.rst @@ -58,7 +58,9 @@ Description Apply bond and angle constraints to specified bonds and angles in the simulation by either the SHAKE or RATTLE algorithms. This typically -enables a longer timestep. +enables a longer timestep. This SHAKE or RATTLE algorithms can *only* +be applied during molecular dynamics runs. When this fix is used during +a minimization, the constraints are replaced by strong harmonic restraints. **SHAKE vs RATTLE:** @@ -166,10 +168,13 @@ 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. +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. How well the geometries +are maintained and how quickly a minimization will converge depends on +the magnitude of the force constant (kbond). If it is chosen too large +the minimization may converge slowly. The default is 1.0e6*k_B. ---------- @@ -202,6 +207,9 @@ Restart, fix_modify, output, run start/stop, minimize info No information about these fixes is written to :doc:`binary restart files `. +When used during minimization, the SHAKE or RATTLE algorithms are **not** +applied. Strong restraint forces are applied instead. + The :doc:`fix_modify ` *virial* option is supported by these fixes to add the contribution due to the added forces on atoms to both the global pressure and per-atom stress of the system via the @@ -209,14 +217,15 @@ to both the global pressure and per-atom stress of the system via the stress/atom ` commands. The former can be accessed by :doc:`thermodynamic output `. The default setting for this fix is :doc:`fix_modify virial yes `. +During minimization, the virial contribution is *NOT* available. -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 +No global or per-atom quantities are stored by these fixes for access by +various :doc:`output commands ` during a run. During +minimization, this fix computes a global scalar which is the energy of +the restraint forces applied insteat of the constraints. No parameter +of these fixes can be used with the *start/stop* keywords of the :doc:`run ` command. -When used during minimization, the SHAKE or RATTLE algorithms are **not** -applied. Strong restraint forces are applied instead. Restrictions """""""""""" @@ -249,7 +258,7 @@ Related commands Default """"""" -kbond = 1.0e6 +kbond = 1.0e9*k_B ---------- diff --git a/src/RIGID/fix_shake.cpp b/src/RIGID/fix_shake.cpp index 41d8c1599c..99836483dd 100644 --- a/src/RIGID/fix_shake.cpp +++ b/src/RIGID/fix_shake.cpp @@ -52,16 +52,17 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) : loop_respa(nullptr), step_respa(nullptr), x(nullptr), v(nullptr), f(nullptr), ftmp(nullptr), vtmp(nullptr), mass(nullptr), rmass(nullptr), type(nullptr), shake_flag(nullptr), shake_atom(nullptr), shake_type(nullptr), xshake(nullptr), nshake(nullptr), list(nullptr), - b_count(nullptr), b_count_all(nullptr), b_atom(nullptr), b_atom_all(nullptr), b_ave(nullptr), b_max(nullptr), - b_min(nullptr), b_ave_all(nullptr), b_max_all(nullptr), b_min_all(nullptr), a_count(nullptr), - 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) + b_count(nullptr), b_count_all(nullptr), b_atom(nullptr), b_atom_all(nullptr), b_ave(nullptr), + b_max(nullptr), b_min(nullptr), b_ave_all(nullptr), b_max_all(nullptr), b_min_all(nullptr), + a_count(nullptr), 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) { energy_global_flag = energy_peratom_flag = 1; virial_global_flag = virial_peratom_flag = 1; thermo_energy = thermo_virial = 1; create_attribute = 1; dof_flag = 1; + scalar_flag = 1; stores_ids = 1; centroidstressflag = CENTROID_AVAIL; next_output = -1; @@ -162,7 +163,7 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) : // parse optional args onemols = nullptr; - kbond = 1.0e6; + kbond = 1.0e6*force->boltz; int iarg = next; while (iarg < narg) { @@ -328,6 +329,7 @@ int FixShake::setmask() mask |= POST_FORCE; mask |= POST_FORCE_RESPA; mask |= MIN_POST_FORCE; + mask |= POST_RUN; return mask; } @@ -351,8 +353,8 @@ void FixShake::init() // that should contribute to potential energy 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->warning(FLERR,"Using fix {} with minimization.\n Substituting constraints with " + "harmonic restraint forces using kbond={:.4g}", style, kbond); // error if a fix changing the box comes before shake fix bool boxflag = false; @@ -590,6 +592,7 @@ void FixShake::post_force(int vflag) // virial setup v_init(vflag); + ebond = 0.0; // loop over clusters to add constraint forces @@ -669,13 +672,12 @@ void FixShake::min_post_force(int vflag) 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; + ebond = 0.0; - // loop over clusters to add strong restraint forces + // loop over shake clusters to add restraint forces for (int i = 0; i < nlist; i++) { int m = list[i]; @@ -2506,7 +2508,6 @@ void FixShake::shake3angle(int m) void FixShake::bond_force(tagint id1, tagint id2, double length) { - int i1 = atom->map(id1); int i2 = atom->map(id2); @@ -2525,25 +2526,18 @@ void FixShake::bond_force(tagint id1, tagint id2, double length) 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); + ebond += 0.5*rk*dr; } if (i2 < nlocal) { f[i2][0] -= delx * fbond; f[i2][1] -= dely * fbond; f[i2][2] -= delz * fbond; - if (evflag) v_tally(i2, v); + ebond += 0.5*rk*dr; } } @@ -3102,6 +3096,26 @@ void *FixShake::extract(const char *str, int &dim) return nullptr; } +/* ---------------------------------------------------------------------- + energy due to restraint forces +------------------------------------------------------------------------- */ + +double FixShake::compute_scalar() +{ + double all; + MPI_Allreduce(&ebond, &all, 1, MPI_DOUBLE, MPI_SUM, world); + return all; +} + +/* ---------------------------------------------------------------------- + print shake stats at the end of a minimization +------------------------------------------------------------------------- */ +void FixShake::post_run() +{ + if ((update->whichflag == 2) && (output_every > 0)) stats(); +} + + /* ---------------------------------------------------------------------- add coordinate constraining forces this method is called at the end of a timestep diff --git a/src/RIGID/fix_shake.h b/src/RIGID/fix_shake.h index a0f743d7d2..d6f9a8f5d6 100644 --- a/src/RIGID/fix_shake.h +++ b/src/RIGID/fix_shake.h @@ -39,6 +39,7 @@ class FixShake : public Fix { void post_force(int) override; void post_force_respa(int, int, int) override; void min_post_force(int) override; + void post_run() override; double memory_usage() override; void grow_arrays(int) override; @@ -59,6 +60,7 @@ class FixShake : public Fix { int dof(int) override; void reset_dt() override; void *extract(const char *, int &) override; + double compute_scalar() override; protected: int vflag_post_force; // store the vflag of last post_force call @@ -78,6 +80,7 @@ class FixShake : public Fix { int molecular; // copy of atom->molecular double *bond_distance, *angle_distance; // constraint distances double kbond; // force constant for restraint + double ebond; // energy of bond restraints class FixRespa *fix_respa; // rRESPA fix needed by SHAKE int nlevels_respa; // copies of needed rRESPA variables From 1ee35bea6117c0df955326d4857a7a4475ff696e Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 10 Jun 2022 01:41:14 -0400 Subject: [PATCH 03/10] fix shake stats (again) --- src/RIGID/fix_shake.cpp | 60 ++++++++++++++++++++--------------------- src/RIGID/fix_shake.h | 3 +-- 2 files changed, 30 insertions(+), 33 deletions(-) diff --git a/src/RIGID/fix_shake.cpp b/src/RIGID/fix_shake.cpp index 99836483dd..c284d99011 100644 --- a/src/RIGID/fix_shake.cpp +++ b/src/RIGID/fix_shake.cpp @@ -52,10 +52,10 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) : loop_respa(nullptr), step_respa(nullptr), x(nullptr), v(nullptr), f(nullptr), ftmp(nullptr), vtmp(nullptr), mass(nullptr), rmass(nullptr), type(nullptr), shake_flag(nullptr), shake_atom(nullptr), shake_type(nullptr), xshake(nullptr), nshake(nullptr), list(nullptr), - b_count(nullptr), b_count_all(nullptr), b_atom(nullptr), b_atom_all(nullptr), b_ave(nullptr), - b_max(nullptr), b_min(nullptr), b_ave_all(nullptr), b_max_all(nullptr), b_min_all(nullptr), - a_count(nullptr), 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) + b_count(nullptr), b_count_all(nullptr), b_ave(nullptr), b_max(nullptr), b_min(nullptr), + b_ave_all(nullptr), b_max_all(nullptr), b_min_all(nullptr), a_count(nullptr), + 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) { energy_global_flag = energy_peratom_flag = 1; virial_global_flag = virial_peratom_flag = 1; @@ -204,8 +204,6 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) : int nb = atom->nbondtypes + 1; b_count = new int[nb]; b_count_all = new int[nb]; - b_atom = new int[nb]; - b_atom_all = new int[nb]; b_ave = new double[nb]; b_ave_all = new double[nb]; b_max = new double[nb]; @@ -298,8 +296,6 @@ FixShake::~FixShake() if (output_every) { delete[] b_count; delete[] b_count_all; - delete[] b_atom; - delete[] b_atom_all; delete[] b_ave; delete[] b_ave_all; delete[] b_max; @@ -2547,7 +2543,6 @@ void FixShake::bond_force(tagint id1, tagint id2, double length) void FixShake::stats() { - int i,j,m,n,iatom,jatom,katom; double delx,dely,delz; double r,r1,r2,r3,angle; @@ -2556,13 +2551,12 @@ void FixShake::stats() int nb = atom->nbondtypes + 1; int na = atom->nangletypes + 1; - for (i = 0; i < nb; i++) { + for (int i = 0; i < nb; i++) { b_count[i] = 0; b_ave[i] = b_max[i] = 0.0; b_min[i] = BIG; - b_atom[i] = -1; } - for (i = 0; i < na; i++) { + for (int i = 0; i < na; i++) { a_count[i] = 0; a_ave[i] = a_max[i] = 0.0; a_min[i] = BIG; @@ -2574,25 +2568,26 @@ void FixShake::stats() double **x = atom->x; int nlocal = atom->nlocal; - for (i = 0; i < nlocal; i++) { - if (shake_flag[i] == 0) continue; + for (int ii = 0; ii < nlist; ++ii) { + int i = list[ii]; + int n = shake_flag[i]; + if (n == 0) continue; // bond stats - n = shake_flag[i]; if (n == 1) n = 3; - iatom = atom->map(shake_atom[i][0]); - for (j = 1; j < n; j++) { - jatom = atom->map(shake_atom[i][j]); + int iatom = atom->map(shake_atom[i][0]); + for (int j = 1; j < n; j++) { + int jatom = atom->map(shake_atom[i][j]); + if (jatom >= nlocal) continue; delx = x[iatom][0] - x[jatom][0]; dely = x[iatom][1] - x[jatom][1]; delz = x[iatom][2] - x[jatom][2]; domain->minimum_image(delx,dely,delz); - r = sqrt(delx*delx + dely*dely + delz*delz); - m = shake_type[i][j-1]; + r = sqrt(delx*delx + dely*dely + delz*delz); + int m = shake_type[i][j-1]; b_count[m]++; - b_atom[m] = n; b_ave[m] += r; b_max[m] = MAX(b_max[m],r); b_min[m] = MIN(b_min[m],r); @@ -2601,9 +2596,13 @@ void FixShake::stats() // angle stats if (shake_flag[i] == 1) { - iatom = atom->map(shake_atom[i][0]); - jatom = atom->map(shake_atom[i][1]); - katom = atom->map(shake_atom[i][2]); + int iatom = atom->map(shake_atom[i][0]); + int jatom = atom->map(shake_atom[i][1]); + int katom = atom->map(shake_atom[i][2]); + int n = 0; + if (iatom < nlocal) ++n; + if (jatom < nlocal) ++n; + if (katom < nlocal) ++n; delx = x[iatom][0] - x[jatom][0]; dely = x[iatom][1] - x[jatom][1]; @@ -2625,9 +2624,9 @@ void FixShake::stats() angle = acos((r1*r1 + r2*r2 - r3*r3) / (2.0*r1*r2)); angle *= 180.0/MY_PI; - m = shake_type[i][2]; - a_count[m]++; - a_ave[m] += angle; + int m = shake_type[i][2]; + a_count[m] += n; + a_ave[m] += n*angle; a_max[m] = MAX(a_max[m],angle); a_min[m] = MIN(a_min[m],angle); } @@ -2636,7 +2635,6 @@ void FixShake::stats() // sum across all procs MPI_Allreduce(b_count,b_count_all,nb,MPI_INT,MPI_SUM,world); - MPI_Allreduce(b_atom,b_atom_all,nb,MPI_INT,MPI_MAX,world); MPI_Allreduce(b_ave,b_ave_all,nb,MPI_DOUBLE,MPI_SUM,world); MPI_Allreduce(b_max,b_max_all,nb,MPI_DOUBLE,MPI_MAX,world); MPI_Allreduce(b_min,b_min_all,nb,MPI_DOUBLE,MPI_MIN,world); @@ -2652,13 +2650,13 @@ void FixShake::stats() const int width = log10((MAX(MAX(1,nb),na)))+2; auto mesg = fmt::format("{} stats (type/ave/delta/count) on step {}\n", utils::uppercase(style), update->ntimestep); - for (i = 1; i < nb; i++) { + for (int i = 1; i < nb; i++) { const auto bcnt = b_count_all[i]; if (bcnt) mesg += fmt::format("Bond: {:>{}d} {:<9.6} {:<11.6} {:>8d}\n",i,width, - b_ave_all[i]/bcnt,b_max_all[i]-b_min_all[i],bcnt/b_atom_all[i]); + b_ave_all[i]/bcnt,b_max_all[i]-b_min_all[i],bcnt); } - for (i = 1; i < na; i++) { + for (int i = 1; i < na; i++) { const auto acnt = a_count_all[i]; if (acnt) mesg += fmt::format("Angle: {:>{}d} {:<9.6} {:<11.6} {:>8d}\n",i,width, diff --git a/src/RIGID/fix_shake.h b/src/RIGID/fix_shake.h index d6f9a8f5d6..820d68ebe7 100644 --- a/src/RIGID/fix_shake.h +++ b/src/RIGID/fix_shake.h @@ -113,8 +113,7 @@ class FixShake : public Fix { int nlist, maxlist; // size and max-size of list // stat quantities - int *b_count, *b_count_all, *b_atom, - *b_atom_all; // counts for each bond type, atoms in bond cluster + int *b_count, *b_count_all; // counts for each bond type, atoms in bond cluster double *b_ave, *b_max, *b_min; // ave/max/min dist for each bond type double *b_ave_all, *b_max_all, *b_min_all; // MPI summing arrays int *a_count, *a_count_all; // ditto for angle types From c4a76103667823d30ce4ce1b4d89ff90889270a0 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 10 Jun 2022 11:07:50 -0400 Subject: [PATCH 04/10] update docs and include suggestions --- doc/src/fix_shake.rst | 56 ++++++++++++--------- doc/utils/sphinx-config/false_positives.txt | 1 + 2 files changed, 33 insertions(+), 24 deletions(-) diff --git a/doc/src/fix_shake.rst b/doc/src/fix_shake.rst index 137cb1aedb..4ba2eb7f4d 100644 --- a/doc/src/fix_shake.rst +++ b/doc/src/fix_shake.rst @@ -168,13 +168,17 @@ 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. How well the geometries -are maintained and how quickly a minimization will converge depends on -the magnitude of the force constant (kbond). If it is chosen too large -the minimization may converge slowly. The default is 1.0e6*k_B. +The *kbond* keyword sets 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 +maintain the geometries similar to the constraints. How well the +geometries are maintained and how quickly a minimization converges, +depends on the force constant *kbond*: larger values will reduce the +deviation from the desired geometry, but can also lead to slower +convergence of the minimization or lead to instabilities depending on +the minimization algorithm requiring to reduce the value of +:doc:`timestep `. The default value for *kbond* depends on +the :doc:`units ` setting and is 1.0e6*k_B. ---------- @@ -190,7 +194,7 @@ LAMMPS closely follows (:ref:`Andersen (1983) `). .. note:: - The fix rattle command modifies forces and velocities and thus + The *fix rattle* command modifies forces and velocities and thus should be defined after all other integration fixes in your input script. If you define other fixes that modify velocities or forces after fix rattle operates, then fix rattle will not take them into @@ -207,24 +211,28 @@ Restart, fix_modify, output, run start/stop, minimize info No information about these fixes is written to :doc:`binary restart files `. -When used during minimization, the SHAKE or RATTLE algorithms are **not** -applied. Strong restraint forces are applied instead. +Fix *shake* and *rattle* behave differently during minimization and +during a molecular dynamics run. -The :doc:`fix_modify ` *virial* option is supported by -these fixes to add the contribution due to the added forces on atoms -to both the global pressure and per-atom stress of the system via the -:doc:`compute pressure ` and :doc:`compute -stress/atom ` commands. The former can be -accessed by :doc:`thermodynamic output `. The default -setting for this fix is :doc:`fix_modify virial yes `. -During minimization, the virial contribution is *NOT* available. +When used during minimization, the SHAKE or RATTLE algorithms are +**not** applied. The constraints are replaced by restraint forces +instead. The energy due to restraint forces is included in the global +potential energy, but virial contributions from them are not included in +the global pressure. The restraint energy is also accessible as a +global scalar property of the fix. -No global or per-atom quantities are stored by these fixes for access by -various :doc:`output commands ` during a run. During -minimization, this fix computes a global scalar which is the energy of -the restraint forces applied insteat of the constraints. No parameter -of these fixes can be used with the *start/stop* keywords of the -:doc:`run ` command. +During molecular dynamics runs, the fixes apply the requested +constraints. The :doc:`fix_modify ` *virial* option is in +this case supported by these fixes to add the contribution due to the +added constraint forces on atoms to both the global pressure and +per-atom stress of the system via the :doc:`compute pressure +` and :doc:`compute stress/atom ` +commands. The former can be accessed by :doc:`thermodynamic output +`. The default 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 ` +during a run. No parameter of these fixes can be used with the +*start/stop* keywords of the :doc:`run ` command. Restrictions diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index b42cff262a..ff1d17ee07 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -1609,6 +1609,7 @@ kb kB kbit kbits +kbond kcal kcl Kd From f05fcaf0d5f8b086a0deaf07839bd8c1686ba68e Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 17 Jun 2022 05:31:42 -0400 Subject: [PATCH 05/10] change energy tally during minimize --- src/RIGID/fix_shake.cpp | 45 ++++++++++++++++++++++++++-- src/RIGID/fix_shake.h | 3 ++ src/compute_pe_atom.cpp | 66 ++++++++++++++++++++++------------------- 3 files changed, 82 insertions(+), 32 deletions(-) diff --git a/src/RIGID/fix_shake.cpp b/src/RIGID/fix_shake.cpp index c284d99011..5b2cc7a33a 100644 --- a/src/RIGID/fix_shake.cpp +++ b/src/RIGID/fix_shake.cpp @@ -67,6 +67,11 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) : centroidstressflag = CENTROID_AVAIL; next_output = -1; + // to avoid uninitialized access + vflag_post_force = 0; + eflag_pre_reverse = 0; + ebond = 0.0; + // error check molecular = atom->molecular; @@ -324,6 +329,7 @@ int FixShake::setmask() mask |= PRE_NEIGHBOR; mask |= POST_FORCE; mask |= POST_FORCE_RESPA; + mask |= MIN_PRE_REVERSE; mask |= MIN_POST_FORCE; mask |= POST_RUN; return mask; @@ -505,6 +511,13 @@ void FixShake::min_setup(int vflag) min_post_force(vflag); } +/* --------------------------------------------------------------------- */ + +void FixShake::setup_pre_reverse(int eflag, int vflag) +{ + min_pre_reverse(eflag,vflag); +} + /* ---------------------------------------------------------------------- build list of SHAKE clusters to constrain if one or more atoms in cluster are on this proc, @@ -513,6 +526,7 @@ void FixShake::min_setup(int vflag) void FixShake::pre_neighbor() { + ebond = 0.0; int atom1,atom2,atom3,atom4; // local copies of atom quantities @@ -653,6 +667,15 @@ void FixShake::post_force_respa(int vflag, int ilevel, int iloop) vflag_post_force = vflag; } +/* ---------------------------------------------------------------------- + store eflag so it can be used in min_post_force +------------------------------------------------------------------------- */ + +void FixShake::min_pre_reverse(int eflag, int /*vflag*/) +{ + eflag_pre_reverse = eflag; +} + /* ---------------------------------------------------------------------- substitute shake constraints with very strong bonds ------------------------------------------------------------------------- */ @@ -668,6 +691,9 @@ void FixShake::min_post_force(int vflag) next_output = (ntimestep/output_every)*output_every + output_every; } else next_output = -1; + int eflag = eflag_pre_reverse; + ev_init(eflag, vflag); + x = atom->x; f = atom->f; nlocal = atom->nlocal; @@ -2522,18 +2548,33 @@ void FixShake::bond_force(tagint id1, tagint id2, double length) const double dr = r - length; const double rk = kbond * dr; const double fbond = (r > 0.0) ? -2.0 * rk / r : 0.0; + const double eb = rk*dr; + int list[2]; + int nlist = 0; if (i1 < nlocal) { f[i1][0] += delx * fbond; f[i1][1] += dely * fbond; f[i1][2] += delz * fbond; - ebond += 0.5*rk*dr; + list[nlist++] = i1; + ebond += 0.5*eb; } if (i2 < nlocal) { f[i2][0] -= delx * fbond; f[i2][1] -= dely * fbond; f[i2][2] -= delz * fbond; - ebond += 0.5*rk*dr; + list[nlist++] = i2; + ebond += 0.5*eb; + } + if (evflag) { + 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; + ev_tally(nlist, list, 2.0, eb, v); } } diff --git a/src/RIGID/fix_shake.h b/src/RIGID/fix_shake.h index 820d68ebe7..914b73ea34 100644 --- a/src/RIGID/fix_shake.h +++ b/src/RIGID/fix_shake.h @@ -34,10 +34,12 @@ class FixShake : public Fix { int setmask() override; void init() override; void setup(int) override; + void setup_pre_reverse(int, 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_pre_reverse(int, int) override; void min_post_force(int) override; void post_run() override; @@ -64,6 +66,7 @@ class FixShake : public Fix { protected: int vflag_post_force; // store the vflag of last post_force call + int eflag_pre_reverse; // store the eflag of last pre_reverse call int respa; // 0 = vel. Verlet, 1 = respa int rattle; // 0 = SHAKE, 1 = RATTLE double tolerance; // SHAKE tolerance diff --git a/src/compute_pe_atom.cpp b/src/compute_pe_atom.cpp index a627e133e5..b1b62ec16d 100644 --- a/src/compute_pe_atom.cpp +++ b/src/compute_pe_atom.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -13,30 +12,31 @@ ------------------------------------------------------------------------- */ #include "compute_pe_atom.h" -#include -#include "atom.h" -#include "update.h" -#include "comm.h" -#include "force.h" -#include "pair.h" -#include "bond.h" + #include "angle.h" +#include "atom.h" +#include "bond.h" +#include "comm.h" #include "dihedral.h" +#include "error.h" +#include "force.h" #include "improper.h" #include "kspace.h" -#include "modify.h" #include "memory.h" -#include "error.h" +#include "modify.h" +#include "pair.h" +#include "update.h" + +#include using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ ComputePEAtom::ComputePEAtom(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), - energy(nullptr) + Compute(lmp, narg, arg), energy(nullptr) { - if (narg < 3) error->all(FLERR,"Illegal compute pe/atom command"); + if (narg < 3) error->all(FLERR, "Illegal compute pe/atom command"); peratom_flag = 1; size_peratom_cols = 0; @@ -56,14 +56,22 @@ ComputePEAtom::ComputePEAtom(LAMMPS *lmp, int narg, char **arg) : fixflag = 0; int iarg = 3; while (iarg < narg) { - if (strcmp(arg[iarg],"pair") == 0) pairflag = 1; - else if (strcmp(arg[iarg],"bond") == 0) bondflag = 1; - else if (strcmp(arg[iarg],"angle") == 0) angleflag = 1; - else if (strcmp(arg[iarg],"dihedral") == 0) dihedralflag = 1; - else if (strcmp(arg[iarg],"improper") == 0) improperflag = 1; - else if (strcmp(arg[iarg],"kspace") == 0) kspaceflag = 1; - else if (strcmp(arg[iarg],"fix") == 0) fixflag = 1; - else error->all(FLERR,"Illegal compute pe/atom command"); + if (strcmp(arg[iarg], "pair") == 0) + pairflag = 1; + else if (strcmp(arg[iarg], "bond") == 0) + bondflag = 1; + else if (strcmp(arg[iarg], "angle") == 0) + angleflag = 1; + else if (strcmp(arg[iarg], "dihedral") == 0) + dihedralflag = 1; + else if (strcmp(arg[iarg], "improper") == 0) + improperflag = 1; + else if (strcmp(arg[iarg], "kspace") == 0) + kspaceflag = 1; + else if (strcmp(arg[iarg], "fix") == 0) + fixflag = 1; + else + error->all(FLERR, "Illegal compute pe/atom command"); iarg++; } } @@ -86,7 +94,7 @@ void ComputePEAtom::compute_peratom() invoked_peratom = update->ntimestep; if (update->eflag_atom != invoked_peratom) - error->all(FLERR,"Per-atom energy was not tallied on needed timestep"); + error->all(FLERR, "Per-atom energy was not tallied on needed timestep"); // grow local energy array if necessary // needs to be atom->nmax in length @@ -94,7 +102,7 @@ void ComputePEAtom::compute_peratom() if (atom->nmax > nmax) { memory->destroy(energy); nmax = atom->nmax; - memory->create(energy,nmax,"pe/atom:energy"); + memory->create(energy, nmax, "pe/atom:energy"); vector_atom = energy; } @@ -153,13 +161,11 @@ void ComputePEAtom::compute_peratom() // add in per-atom contributions from relevant fixes // always only for owned atoms, not ghost - if (fixflag && modify->n_energy_atom) - modify->energy_atom(nlocal,energy); + if (fixflag && modify->n_energy_atom) modify->energy_atom(nlocal, energy); // communicate ghost energy between neighbor procs - if (force->newton || (force->kspace && force->kspace->tip4pflag)) - comm->reverse_comm(this); + if (force->newton || (force->kspace && force->kspace->tip4pflag)) comm->reverse_comm(this); // zero energy of atoms not in group // only do this after comm since ghost contributions must be included @@ -174,7 +180,7 @@ void ComputePEAtom::compute_peratom() int ComputePEAtom::pack_reverse_comm(int n, int first, double *buf) { - int i,m,last; + int i, m, last; m = 0; last = first + n; @@ -186,7 +192,7 @@ int ComputePEAtom::pack_reverse_comm(int n, int first, double *buf) void ComputePEAtom::unpack_reverse_comm(int n, int *list, double *buf) { - int i,j,m; + int i, j, m; m = 0; for (i = 0; i < n; i++) { @@ -201,6 +207,6 @@ void ComputePEAtom::unpack_reverse_comm(int n, int *list, double *buf) double ComputePEAtom::memory_usage() { - double bytes = (double)nmax * sizeof(double); + double bytes = (double) nmax * sizeof(double); return bytes; } From 0ad45a02249d6c60e2c80a34330fb76f08f6f620 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 17 Jun 2022 05:53:34 -0400 Subject: [PATCH 06/10] correctly produce eatom (=0) for MD runs --- src/RIGID/fix_shake.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/RIGID/fix_shake.cpp b/src/RIGID/fix_shake.cpp index 5b2cc7a33a..6673472458 100644 --- a/src/RIGID/fix_shake.cpp +++ b/src/RIGID/fix_shake.cpp @@ -601,7 +601,8 @@ void FixShake::post_force(int vflag) // virial setup - v_init(vflag); + int eflag = eflag_pre_reverse; + ev_init(eflag, vflag); ebond = 0.0; // loop over clusters to add constraint forces From e66229dadb61f60c784c7f1863331a6afefcf4de Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 17 Jun 2022 06:19:57 -0400 Subject: [PATCH 07/10] update docs --- doc/src/fix_shake.rst | 37 ++++++++++++++++++++++--------------- 1 file changed, 22 insertions(+), 15 deletions(-) diff --git a/doc/src/fix_shake.rst b/doc/src/fix_shake.rst index 4ba2eb7f4d..b53f8602ec 100644 --- a/doc/src/fix_shake.rst +++ b/doc/src/fix_shake.rst @@ -216,23 +216,25 @@ during a molecular dynamics run. When used during minimization, the SHAKE or RATTLE algorithms are **not** applied. The constraints are replaced by restraint forces -instead. The energy due to restraint forces is included in the global -potential energy, but virial contributions from them are not included in -the global pressure. The restraint energy is also accessible as a -global scalar property of the fix. +instead. The energy and virial contributions due to the restraint +forces are tallied into global and per-atom accumulators. The total +restraint energy is also accessible as a global scalar property of the +fix. During molecular dynamics runs, the fixes apply the requested -constraints. The :doc:`fix_modify ` *virial* option is in -this case supported by these fixes to add the contribution due to the -added constraint forces on atoms to both the global pressure and -per-atom stress of the system via the :doc:`compute pressure -` and :doc:`compute stress/atom ` -commands. The former can be accessed by :doc:`thermodynamic output -`. The default 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 ` -during a run. No parameter of these fixes can be used with the -*start/stop* keywords of the :doc:`run ` command. +constraints. + +The :doc:`fix_modify ` *virial* option is supported by these +fixes to add the contribution due to the added constraint forces on +atoms to both the global pressure and per-atom stress of the system via +the :doc:`compute pressure ` and :doc:`compute +stress/atom ` commands. The former can be accessed +by :doc:`thermodynamic output `. The default 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 ` during an MD run. No parameter of +these fixes can be used with the *start/stop* keywords of the :doc:`run +` command. Restrictions @@ -256,6 +258,11 @@ degrees (e.g. linear CO2 molecule). This causes numeric difficulties. You can use :doc:`fix rigid or fix rigid/small ` instead to make a linear molecule rigid. +When used during minimization choosing a too large value of the *kbond* +can make minimization very inefficent and also cause stability problems +with some minimization algorithms. Sometimes those can be avoided by +reducing the :doc:`timestep `. + Related commands """""""""""""""" From 114b19f620b3d7904cfe9bb5cbf1a1621767bd40 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 23 Jun 2022 06:51:13 -0400 Subject: [PATCH 08/10] make certain that the fix energy is properly reset to zero --- src/KOKKOS/fix_shake_kokkos.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/KOKKOS/fix_shake_kokkos.cpp b/src/KOKKOS/fix_shake_kokkos.cpp index 5ed8cfccb3..40f641780a 100644 --- a/src/KOKKOS/fix_shake_kokkos.cpp +++ b/src/KOKKOS/fix_shake_kokkos.cpp @@ -192,6 +192,7 @@ void FixShakeKokkos::pre_neighbor() // local copies of atom quantities // used by SHAKE until next re-neighboring + ebond = 0.0; x = atom->x; v = atom->v; f = atom->f; @@ -302,6 +303,7 @@ void FixShakeKokkos::pre_neighbor() template void FixShakeKokkos::post_force(int vflag) { + ebond = 0.0; copymode = 1; d_x = atomKK->k_x.view(); From caf21d09b43c2820518e00e033168e767b10119b Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 5 Aug 2022 18:12:29 -0400 Subject: [PATCH 09/10] disallow using fix shake/rattle with KOKKOS during minimization (for now) --- src/KOKKOS/fix_shake_kokkos.cpp | 10 ++++++++++ src/KOKKOS/fix_shake_kokkos.h | 1 + 2 files changed, 11 insertions(+) diff --git a/src/KOKKOS/fix_shake_kokkos.cpp b/src/KOKKOS/fix_shake_kokkos.cpp index 40f641780a..d622da94ec 100644 --- a/src/KOKKOS/fix_shake_kokkos.cpp +++ b/src/KOKKOS/fix_shake_kokkos.cpp @@ -180,6 +180,16 @@ void FixShakeKokkos::init() } +/* ---------------------------------------------------------------------- + run setup for minimization. +------------------------------------------------------------------------- */ + +template +void FixShakeKokkos::min_setup(int /*vflag*/) +{ + error->all(FLERR, "Cannot yet use fix {} during minimization with Kokkos", style); +} + /* ---------------------------------------------------------------------- build list of SHAKE clusters to constrain if one or more atoms in cluster are on this proc, diff --git a/src/KOKKOS/fix_shake_kokkos.h b/src/KOKKOS/fix_shake_kokkos.h index 027b36b100..df5256422f 100644 --- a/src/KOKKOS/fix_shake_kokkos.h +++ b/src/KOKKOS/fix_shake_kokkos.h @@ -51,6 +51,7 @@ class FixShakeKokkos : public FixShake, public KokkosBase { FixShakeKokkos(class LAMMPS *, int, char **); ~FixShakeKokkos() override; void init() override; + void min_setup(int) override; void pre_neighbor() override; void post_force(int) override; From daa6b78c43af7f32291a7f7c9629ecb69e829b83 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 17 Aug 2022 15:22:59 -0400 Subject: [PATCH 10/10] update documentation --- doc/src/fix_shake.rst | 100 ++++++++++++++++++++++-------------------- 1 file changed, 53 insertions(+), 47 deletions(-) diff --git a/doc/src/fix_shake.rst b/doc/src/fix_shake.rst index b53f8602ec..f5742b8f9f 100644 --- a/doc/src/fix_shake.rst +++ b/doc/src/fix_shake.rst @@ -58,9 +58,10 @@ Description Apply bond and angle constraints to specified bonds and angles in the simulation by either the SHAKE or RATTLE algorithms. This typically -enables a longer timestep. This SHAKE or RATTLE algorithms can *only* -be applied during molecular dynamics runs. When this fix is used during -a minimization, the constraints are replaced by strong harmonic restraints. +enables a longer timestep. The SHAKE or RATTLE algorithms, however, can +*only* be applied during molecular dynamics runs. When this fix is used +during a minimization, the constraints are *approximated* by strong +harmonic restraints. **SHAKE vs RATTLE:** @@ -94,7 +95,7 @@ The SHAKE algorithm satisfies the first condition, i.e. the sites at time *n+1* will have the desired separations Dij immediately after the coordinates are integrated. If we also enforce the second condition, the velocity components along the bonds will vanish. RATTLE satisfies -both conditions. As implemented in LAMMPS, fix rattle uses fix shake +both conditions. As implemented in LAMMPS, *fix rattle* uses fix shake for satisfying the coordinate constraints. Therefore the settings and optional keywords are the same for both fixes, and all the information below about SHAKE is also relevant for RATTLE. @@ -143,16 +144,16 @@ for. .. note:: - This command works by using the current forces on atoms to - calculate an additional constraint force which when added will leave - the atoms in positions that satisfy the SHAKE constraints (e.g. bond - length) after the next time integration step. If you define fixes - (e.g. :doc:`fix efield `) that add additional force to the - atoms after fix shake operates, then this fix will not take them into - account and the time integration will typically not satisfy the SHAKE - constraints. The solution for this is to make sure that fix shake is - defined in your input script after any other fixes which add or change - forces (to atoms that fix shake operates on). + This command works by using the current forces on atoms to calculate + an additional constraint force which when added will leave the atoms + in positions that satisfy the SHAKE constraints (e.g. bond length) + after the next time integration step. If you define fixes + (e.g. :doc:`fix efield `) that add additional force to + the atoms after *fix shake* operates, then this fix will not take them + into account and the time integration will typically not satisfy the + SHAKE constraints. The solution for this is to make sure that fix + shake is defined in your input script after any other fixes which add + or change forces (to atoms that *fix shake* operates on). ---------- @@ -168,17 +169,21 @@ 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 sets the restraint force constant when fix shake or -fix rattle are used during minimization. In that case the constraint +The *kbond* keyword sets 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 maintain the geometries similar to the constraints. How well the geometries are maintained and how quickly a minimization converges, -depends on the force constant *kbond*: larger values will reduce the -deviation from the desired geometry, but can also lead to slower +depends largely on the force constant *kbond*: larger values will reduce +the deviation from the desired geometry, but can also lead to slower convergence of the minimization or lead to instabilities depending on the minimization algorithm requiring to reduce the value of -:doc:`timestep `. The default value for *kbond* depends on -the :doc:`units ` setting and is 1.0e6*k_B. +:doc:`timestep `. Even though the restraints will not +preserve the bond lengths and angles as closely as the constraints +during the MD, they are generally close enough so that the constraints +will be fulfilled to the desired accuracy within a few MD steps +following the minimization. The default value for *kbond* depends on the +:doc:`units ` setting and is 1.0e6*k_B. ---------- @@ -197,7 +202,7 @@ LAMMPS closely follows (:ref:`Andersen (1983) `). The *fix rattle* command modifies forces and velocities and thus should be defined after all other integration fixes in your input script. If you define other fixes that modify velocities or forces - after fix rattle operates, then fix rattle will not take them into + after *fix rattle* operates, then *fix rattle* will not take them into account and the overall time integration will typically not satisfy the RATTLE constraints. You can check whether the constraints work correctly by setting the value of RATTLE_DEBUG in src/fix_rattle.cpp @@ -211,30 +216,31 @@ Restart, fix_modify, output, run start/stop, minimize info No information about these fixes is written to :doc:`binary restart files `. -Fix *shake* and *rattle* behave differently during minimization and -during a molecular dynamics run. +Both fix *shake* and fix *rattle* behave differently during a minimization +in comparison to a molecular dynamics run: -When used during minimization, the SHAKE or RATTLE algorithms are -**not** applied. The constraints are replaced by restraint forces -instead. The energy and virial contributions due to the restraint -forces are tallied into global and per-atom accumulators. The total -restraint energy is also accessible as a global scalar property of the -fix. +- When used during a minimization, the SHAKE or RATTLE constraint + algorithms themselves are **not** applied. Instead the constraints + are replaced by harmonic restraint forces. The energy and virial + contributions due to the restraint forces are tallied into global and + per-atom accumulators. The total restraint energy is also accessible + as a global scalar property of the fix. -During molecular dynamics runs, the fixes apply the requested -constraints. +- During molecular dynamics runs, however, the fixes do apply the + requested SHAKE or RATTLE constraint algorithms. -The :doc:`fix_modify ` *virial* option is supported by these -fixes to add the contribution due to the added constraint forces on -atoms to both the global pressure and per-atom stress of the system via -the :doc:`compute pressure ` and :doc:`compute -stress/atom ` commands. The former can be accessed -by :doc:`thermodynamic output `. The default 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 ` during an MD run. No parameter of -these fixes can be used with the *start/stop* keywords of the :doc:`run -` command. + The :doc:`fix_modify ` *virial* option is supported by + these fixes to add the contribution due to the added constraint forces + on atoms to both the global pressure and per-atom stress of the system + via the :doc:`compute pressure ` and :doc:`compute + stress/atom ` commands. The former can be + accessed by :doc:`thermodynamic output `. + + The default 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 ` during + an MD run. No parameter of these fixes can be used with the + *start/stop* keywords of the :doc:`run ` command. Restrictions @@ -253,13 +259,13 @@ which can lead to poor energy conservation. You can test for this in your system by running a constant NVE simulation with a particular set of SHAKE parameters and monitoring the energy versus time. -SHAKE or RATTLE should not be used to constrain an angle at 180 -degrees (e.g. linear CO2 molecule). This causes numeric difficulties. -You can use :doc:`fix rigid or fix rigid/small ` instead to -make a linear molecule rigid. +SHAKE or RATTLE should not be used to constrain an angle at 180 degrees +(e.g. a linear CO2 molecule). This causes a divergence when solving the +constraint equations numerically. You can use :doc:`fix rigid or fix +rigid/small ` instead to make a linear molecule rigid. When used during minimization choosing a too large value of the *kbond* -can make minimization very inefficent and also cause stability problems +can make minimization very inefficient and also cause stability problems with some minimization algorithms. Sometimes those can be avoided by reducing the :doc:`timestep `.