Added new argument r0stop to fix restrain bond

Added new fix restrain lowerbound harmonic

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 <bond_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 <angle_harmonic>` command.  The potential associated with
 the restraint is

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]);
+      deqstart[nrestrain] = force->numeric(FLERR,arg[iarg+5]);
-      iarg += 6;
+      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 {

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