From 0f00bb1ca6c08c3d5fcc14857f2a629fd9c8292d Mon Sep 17 00:00:00 2001 From: david-castillo Date: Mon, 23 Mar 2020 16:48:08 +0100 Subject: [PATCH 1/6] Added new argument r0stop to fix restrain bond Added new fix restrain lowerbound harmonic --- doc/src/fix_restrain.rst | 33 +++++++++- src/fix_restrain.cpp | 133 ++++++++++++++++++++++++++++++++++++--- src/fix_restrain.h | 4 +- 3 files changed, 160 insertions(+), 10 deletions(-) diff --git a/doc/src/fix_restrain.rst b/doc/src/fix_restrain.rst index 3965cee04d..9b9beaa5ac 100644 --- a/doc/src/fix_restrain.rst +++ b/doc/src/fix_restrain.rst @@ -17,10 +17,18 @@ Syntax .. parsed-literal:: - *bond* args = atom1 atom2 Kstart Kstop r0 + *bond* args = atom1 atom2 Kstart Kstop r0start r0stop atom1,atom2 = IDs of 2 atoms in bond Kstart,Kstop = restraint coefficients at start/end of run (energy units) - r0 = equilibrium bond distance (distance units) + r0start = equilibrium bond distance at start of run (distance units) + r0stop = equilibrium bond distance at end of run (distance units). If not + specified it's assumed to be equal to r0start + *lbond* args = atom1 atom2 Kstart Kstop r0start r0stop + atom1,atom2 = IDs of 2 atoms in bond + Kstart,Kstop = restraint coefficients at start/end of run (energy units) + r0start = equilibrium bond distance at start of run (distance units) + r0stop = equilibrium bond distance at end of run (distance units). If not + specified it's assumed to be equal to r0start *angle* args = atom1 atom2 atom3 Kstart Kstop theta0 atom1,atom2,atom3 = IDs of 3 atoms in angle, atom2 = middle atom Kstart,Kstop = restraint coefficients at start/end of run (energy units) @@ -38,6 +46,7 @@ Examples .. code-block:: LAMMPS fix holdem all restrain bond 45 48 2000.0 2000.0 2.75 + fix holdem all restrain lbond 45 48 2000.0 2000.0 2.75 fix holdem all restrain dihedral 1 2 3 4 2000.0 2000.0 120.0 fix holdem all restrain bond 45 48 2000.0 2000.0 2.75 dihedral 1 2 3 4 2000.0 2000.0 120.0 fix texas_holdem all restrain dihedral 1 2 3 4 0.0 2000.0 120.0 dihedral 1 2 3 5 0.0 2000.0 -120.0 dihedral 1 2 3 6 0.0 2000.0 0.0 @@ -141,6 +150,26 @@ is included in :math:`K`. ---------- +The *lbond* keyword applies a lowerbound bond restraint to the specified atoms +using the same functional form used by the :doc:`bond_style harmonic ` command if the distance between +the atoms is smaller than the equilibrium bond distance and 0 otherwise. The potential associated with +the restraint is + +.. math:: + + E = K (r - r_0)^2 if r < r_0 + E = 0 if r >= r_0 + +with the following coefficients: + +* :math:`K` (energy/distance\^2) +* :math:`r_0` (distance) + +:math:`K` and :math:`r_0` are specified with the fix. Note that the usual 1/2 factor +is included in :math:`K`. + +---------- + The *angle* keyword applies an angle restraint to the specified atoms using the same functional form used by the :doc:`angle_style harmonic ` command. The potential associated with the restraint is diff --git a/src/fix_restrain.cpp b/src/fix_restrain.cpp index 5a0d4fb662..5f65d7069b 100644 --- a/src/fix_restrain.cpp +++ b/src/fix_restrain.cpp @@ -34,7 +34,7 @@ using namespace LAMMPS_NS; using namespace FixConst; using namespace MathConst; -enum{BOND,ANGLE,DIHEDRAL}; +enum{BOND,LBOUND,ANGLE,DIHEDRAL}; #define TOLERANCE 0.05 #define SMALL 0.001 @@ -45,7 +45,7 @@ enum{BOND,ANGLE,DIHEDRAL}; FixRestrain::FixRestrain(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), rstyle(NULL), mult(NULL), ids(NULL), kstart(NULL), kstop(NULL), target(NULL), - cos_target(NULL), sin_target(NULL) + deqstart(NULL), deqstop(NULL), cos_target(NULL), sin_target(NULL) { if (narg < 4) error->all(FLERR,"Illegal fix restrain command"); @@ -72,6 +72,8 @@ FixRestrain::FixRestrain(LAMMPS *lmp, int narg, char **arg) : memory->grow(kstart,maxrestrain,"restrain:kstart"); memory->grow(kstop,maxrestrain,"restrain:kstop"); memory->grow(target,maxrestrain,"restrain:target"); + memory->grow(deqstart,maxrestrain,"restrain:deqstart"); + memory->grow(deqstop,maxrestrain,"restrain:deqstop"); memory->grow(cos_target,maxrestrain,"restrain:cos_target"); memory->grow(sin_target,maxrestrain,"restrain:sin_target"); } @@ -83,8 +85,29 @@ FixRestrain::FixRestrain(LAMMPS *lmp, int narg, char **arg) : ids[nrestrain][1] = force->inumeric(FLERR,arg[iarg+2]); kstart[nrestrain] = force->numeric(FLERR,arg[iarg+3]); kstop[nrestrain] = force->numeric(FLERR,arg[iarg+4]); - target[nrestrain] = force->numeric(FLERR,arg[iarg+5]); - iarg += 6; + deqstart[nrestrain] = force->numeric(FLERR,arg[iarg+5]); + if (iarg+6 == narg) { + deqstop[nrestrain] = force->numeric(FLERR,arg[iarg+5]); + iarg += 6; + } else { + deqstop[nrestrain] = force->numeric(FLERR,arg[iarg+6]); + iarg += 7; + } + } else if (strcmp(arg[iarg],"lbound") == 0) { + if (iarg+6 > narg) error->all(FLERR,"Illegal fix restrain command"); + rstyle[nrestrain] = LBOUND; + ids[nrestrain][0] = force->inumeric(FLERR,arg[iarg+1]); + ids[nrestrain][1] = force->inumeric(FLERR,arg[iarg+2]); + kstart[nrestrain] = force->numeric(FLERR,arg[iarg+3]); + kstop[nrestrain] = force->numeric(FLERR,arg[iarg+4]); + deqstart[nrestrain] = force->numeric(FLERR,arg[iarg+5]); + if (iarg+6 == narg) { + deqstop[nrestrain] = force->numeric(FLERR,arg[iarg+5]); + iarg += 6; + } else { + deqstop[nrestrain] = force->numeric(FLERR,arg[iarg+6]); + iarg += 7; + } } else if (strcmp(arg[iarg],"angle") == 0) { if (iarg+7 > narg) error->all(FLERR,"Illegal fix restrain command"); rstyle[nrestrain] = ANGLE; @@ -139,6 +162,8 @@ FixRestrain::~FixRestrain() memory->destroy(kstart); memory->destroy(kstop); memory->destroy(target); + memory->destroy(deqstart); + memory->destroy(deqstop); memory->destroy(cos_target); memory->destroy(sin_target); } @@ -192,11 +217,13 @@ void FixRestrain::post_force(int /*vflag*/) energy = 0.0; ebond = 0.0; + elbound = 0.0; eangle = 0.0; edihed = 0.0; for (int m = 0; m < nrestrain; m++) if (rstyle[m] == BOND) restrain_bond(m); + else if (rstyle[m] == LBOUND) restrain_lbound(m); else if (rstyle[m] == ANGLE) restrain_angle(m); else if (rstyle[m] == DIHEDRAL) restrain_dihedral(m); } @@ -233,6 +260,7 @@ void FixRestrain::restrain_bond(int m) double delta = update->ntimestep - update->beginstep; if (delta != 0.0) delta /= update->endstep - update->beginstep; double k = kstart[m] + delta * (kstop[m] - kstart[m]); + double deq = deqstart[m] + delta * (deqstop[m] - deqstart[m]); i1 = atom->map(ids[m][0]); i2 = atom->map(ids[m][1]); @@ -269,7 +297,7 @@ void FixRestrain::restrain_bond(int m) rsq = delx*delx + dely*dely + delz*delz; r = sqrt(rsq); - dr = r - target[m]; + dr = r - deq; rk = k * dr; // force & energy @@ -277,7 +305,7 @@ void FixRestrain::restrain_bond(int m) if (r > 0.0) fbond = -2.0*rk/r; else fbond = 0.0; - ebond += rk*dr; + ebond += rk*dr; energy += rk*dr; // apply force to each of 2 atoms @@ -295,6 +323,94 @@ void FixRestrain::restrain_bond(int m) } } +/* ---------------------------------------------------------------------- + apply harmonic lower-bound bond restraints +---------------------------------------------------------------------- */ + +void FixRestrain::restrain_lbound(int m) +{ + int i1,i2; + double delx,dely,delz,fbond; + double rsq,r,dr,rk; + + double **x = atom->x; + double **f = atom->f; + int nlocal = atom->nlocal; + int newton_bond = force->newton_bond; + + double delta = update->ntimestep - update->beginstep; + if (delta != 0.0) delta /= update->endstep - update->beginstep; + double k = kstart[m] + delta * (kstop[m] - kstart[m]); + double deq = deqstart[m] + delta * (deqstop[m] - deqstart[m]); + + i1 = atom->map(ids[m][0]); + i2 = atom->map(ids[m][1]); + + // newton_bond on: only processor owning i2 computes restraint + // newton_bond off: only processors owning either of i1,i2 computes restraint + + if (newton_bond) { + if (i2 == -1 || i2 >= nlocal) return; + if (i1 == -1) { + char str[128]; + sprintf(str, + "Restrain atoms %d %d missing on proc %d at step " BIGINT_FORMAT, + ids[m][0],ids[m][1], + comm->me,update->ntimestep); + error->one(FLERR,str); + } + } else { + if ((i1 == -1 || i1 >= nlocal) && (i2 == -1 || i2 >= nlocal)) return; + if (i1 == -1 || i2 == -1) { + char str[128]; + sprintf(str, + "Restrain atoms %d %d missing on proc %d at step " BIGINT_FORMAT, + ids[m][0],ids[m][1], + comm->me,update->ntimestep); + error->one(FLERR,str); + } + } + + delx = x[i1][0] - x[i2][0]; + dely = x[i1][1] - x[i2][1]; + delz = x[i1][2] - x[i2][2]; + domain->minimum_image(delx,dely,delz); + + rsq = delx*delx + dely*dely + delz*delz; + r = sqrt(rsq); + dr = r - deq; + rk = k * dr; + + // force & energy + + if (dr < 0) { + if (r > 0.0) fbond = -2.0*rk/r; + else fbond = 0.0; + + elbound += rk*dr; + energy += rk*dr; + } else { + fbond = 0.0; + + elbound += 0.0; + energy += 0.0; + } + + // apply force to each of 2 atoms + + if (newton_bond || i1 < nlocal) { + f[i1][0] += delx*fbond; + f[i1][1] += dely*fbond; + f[i1][2] += delz*fbond; + } + + if (newton_bond || i2 < nlocal) { + f[i2][0] -= delx*fbond; + f[i2][1] -= dely*fbond; + f[i2][2] -= delz*fbond; + } +} + /* ---------------------------------------------------------------------- apply harmonic angle restraints ---------------------------------------------------------------------- */ @@ -655,9 +771,12 @@ double FixRestrain::compute_vector(int n) MPI_Allreduce(&ebond,&ebond_all,1,MPI_DOUBLE,MPI_SUM,world); return ebond_all; } else if (n == 1) { + MPI_Allreduce(&elbound,&elbound_all,1,MPI_DOUBLE,MPI_SUM,world); + return elbound_all; + } else if (n == 3) { MPI_Allreduce(&eangle,&eangle_all,1,MPI_DOUBLE,MPI_SUM,world); return eangle_all; - } else if (n == 2) { + } else if (n == 4) { MPI_Allreduce(&edihed,&edihed_all,1,MPI_DOUBLE,MPI_SUM,world); return edihed_all; } else { diff --git a/src/fix_restrain.h b/src/fix_restrain.h index 4572905d46..3497b7daf5 100644 --- a/src/fix_restrain.h +++ b/src/fix_restrain.h @@ -44,14 +44,16 @@ class FixRestrain : public Fix { int *rstyle; int *mult; int **ids; - double *kstart,*kstop,*target; + double *kstart,*kstop,*deqstart,*deqstop,*target; double *cos_target,*sin_target; double energy,energy_all; double ebond,ebond_all; + double elbound,elbound_all; double eangle,eangle_all; double edihed,edihed_all; void restrain_bond(int); + void restrain_lbound(int); void restrain_angle(int); void restrain_dihedral(int); }; From 6209479797e29df62dd17476127e30a64d73d08f Mon Sep 17 00:00:00 2001 From: david-castillo Date: Mon, 23 Mar 2020 17:57:58 +0100 Subject: [PATCH 2/6] Removed some tabs Changed mathjax for multiline formula --- doc/src/fix_restrain.rst | 9 ++++++--- src/fix_restrain.cpp | 16 ++++++++-------- 2 files changed, 14 insertions(+), 11 deletions(-) diff --git a/doc/src/fix_restrain.rst b/doc/src/fix_restrain.rst index 9b9beaa5ac..110b26680d 100644 --- a/doc/src/fix_restrain.rst +++ b/doc/src/fix_restrain.rst @@ -150,16 +150,19 @@ is included in :math:`K`. ---------- -The *lbond* keyword applies a lowerbound bond restraint to the specified atoms +The *lbond* keyword applies a lower bound bond restraint to the specified atoms using the same functional form used by the :doc:`bond_style harmonic ` command if the distance between the atoms is smaller than the equilibrium bond distance and 0 otherwise. The potential associated with the restraint is .. math:: - E = K (r - r_0)^2 if r < r_0 - E = 0 if r >= r_0 + E = K (r - r_0)^2 ,if\;r < r_0 +.. math:: + + E = 0 \qquad\quad\quad ,if\;r \ge r_0 + with the following coefficients: * :math:`K` (energy/distance\^2) diff --git a/src/fix_restrain.cpp b/src/fix_restrain.cpp index 5f65d7069b..a6793a36be 100644 --- a/src/fix_restrain.cpp +++ b/src/fix_restrain.cpp @@ -87,11 +87,11 @@ FixRestrain::FixRestrain(LAMMPS *lmp, int narg, char **arg) : kstop[nrestrain] = force->numeric(FLERR,arg[iarg+4]); deqstart[nrestrain] = force->numeric(FLERR,arg[iarg+5]); if (iarg+6 == narg) { - deqstop[nrestrain] = force->numeric(FLERR,arg[iarg+5]); - iarg += 6; + deqstop[nrestrain] = force->numeric(FLERR,arg[iarg+5]); + iarg += 6; } else { - deqstop[nrestrain] = force->numeric(FLERR,arg[iarg+6]); - iarg += 7; + deqstop[nrestrain] = force->numeric(FLERR,arg[iarg+6]); + iarg += 7; } } else if (strcmp(arg[iarg],"lbound") == 0) { if (iarg+6 > narg) error->all(FLERR,"Illegal fix restrain command"); @@ -102,11 +102,11 @@ FixRestrain::FixRestrain(LAMMPS *lmp, int narg, char **arg) : kstop[nrestrain] = force->numeric(FLERR,arg[iarg+4]); deqstart[nrestrain] = force->numeric(FLERR,arg[iarg+5]); if (iarg+6 == narg) { - deqstop[nrestrain] = force->numeric(FLERR,arg[iarg+5]); - iarg += 6; + deqstop[nrestrain] = force->numeric(FLERR,arg[iarg+5]); + iarg += 6; } else { - deqstop[nrestrain] = force->numeric(FLERR,arg[iarg+6]); - iarg += 7; + deqstop[nrestrain] = force->numeric(FLERR,arg[iarg+6]); + iarg += 7; } } else if (strcmp(arg[iarg],"angle") == 0) { if (iarg+7 > narg) error->all(FLERR,"Illegal fix restrain command"); From e94576ecd8f66c8cd693c032ebb6d4e80af267ff Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 15 Apr 2020 11:21:25 -0400 Subject: [PATCH 3/6] update documentation as requested --- doc/src/fix_restrain.rst | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/doc/src/fix_restrain.rst b/doc/src/fix_restrain.rst index 110b26680d..44f5e7585d 100644 --- a/doc/src/fix_restrain.rst +++ b/doc/src/fix_restrain.rst @@ -17,12 +17,12 @@ Syntax .. parsed-literal:: - *bond* args = atom1 atom2 Kstart Kstop r0start r0stop + *bond* args = atom1 atom2 Kstart Kstop r0start (r0stop) atom1,atom2 = IDs of 2 atoms in bond Kstart,Kstop = restraint coefficients at start/end of run (energy units) r0start = equilibrium bond distance at start of run (distance units) - r0stop = equilibrium bond distance at end of run (distance units). If not - specified it's assumed to be equal to r0start + r0stop = equilibrium bond distance at end of run (optional) (distance units). If not + specified it is assumed to be equal to r0start *lbond* args = atom1 atom2 Kstart Kstop r0start r0stop atom1,atom2 = IDs of 2 atoms in bond Kstart,Kstop = restraint coefficients at start/end of run (energy units) From 0796a7032e4856e276b91923f8bc4422333b66e9 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 15 Apr 2020 11:22:48 -0400 Subject: [PATCH 4/6] same doc change for "lbond" option --- doc/src/fix_restrain.rst | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/doc/src/fix_restrain.rst b/doc/src/fix_restrain.rst index 44f5e7585d..692336bec6 100644 --- a/doc/src/fix_restrain.rst +++ b/doc/src/fix_restrain.rst @@ -23,12 +23,12 @@ Syntax r0start = equilibrium bond distance at start of run (distance units) r0stop = equilibrium bond distance at end of run (optional) (distance units). If not specified it is assumed to be equal to r0start - *lbond* args = atom1 atom2 Kstart Kstop r0start r0stop + *lbond* args = atom1 atom2 Kstart Kstop r0start (r0stop) atom1,atom2 = IDs of 2 atoms in bond Kstart,Kstop = restraint coefficients at start/end of run (energy units) r0start = equilibrium bond distance at start of run (distance units) - r0stop = equilibrium bond distance at end of run (distance units). If not - specified it's assumed to be equal to r0start + r0stop = equilibrium bond distance at end of run (optional) (distance units). If not + specified it is assumed to be equal to r0start *angle* args = atom1 atom2 atom3 Kstart Kstop theta0 atom1,atom2,atom3 = IDs of 3 atoms in angle, atom2 = middle atom Kstart,Kstop = restraint coefficients at start/end of run (energy units) From adf9602ae3b6f81c4fbf87696a318e128772cd2f Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 15 Apr 2020 14:26:37 -0400 Subject: [PATCH 5/6] whitespace fix --- src/KOKKOS/min_kokkos.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/KOKKOS/min_kokkos.cpp b/src/KOKKOS/min_kokkos.cpp index 30b362dc51..e5f22411f6 100644 --- a/src/KOKKOS/min_kokkos.cpp +++ b/src/KOKKOS/min_kokkos.cpp @@ -636,7 +636,7 @@ double MinKokkos::fnorm_inf() double MinKokkos::fnorm_max() { - + double local_norm_max = 0.0; { // local variables for lambda capture From c5f6b5c6873795d742d637cf6cd6eabfe8b5fe05 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 15 Apr 2020 14:29:37 -0400 Subject: [PATCH 6/6] one more whitespace fix --- src/fix_restrain.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/fix_restrain.cpp b/src/fix_restrain.cpp index a6793a36be..f1032c1bc1 100644 --- a/src/fix_restrain.cpp +++ b/src/fix_restrain.cpp @@ -382,11 +382,11 @@ void FixRestrain::restrain_lbound(int m) rk = k * dr; // force & energy - + if (dr < 0) { if (r > 0.0) fbond = -2.0*rk/r; else fbond = 0.0; - + elbound += rk*dr; energy += rk*dr; } else {