add ability to set multiplicity for dihedrals with fix restrain

this introduces an optional keyword, so backward compatibility is preserved.
this closes #884

doc/src/fix_restrain.txt

@@ -24,10 +24,12 @@ keyword = {bond} or {angle} or {dihedral} :l
     atom1,atom2,atom3 = IDs of 3 atoms in angle, atom2 = middle atom
     Kstart,Kstop = restraint coefficients at start/end of run (energy units)
     theta0 = equilibrium angle theta (degrees)
-  {dihedral} args = atom1 atom2 atom3 atom4 Kstart Kstop phi0
+  {dihedral} args = atom1 atom2 atom3 atom4 Kstart Kstop phi0 keyword/value
     atom1,atom2,atom3,atom4 = IDs of 4 atoms in dihedral in linear order
     Kstart,Kstop = restraint coefficients at start/end of run (energy units)
-    phi0 = equilibrium dihedral angle phi (degrees) :pre
+    phi0 = equilibrium dihedral angle phi (degrees)
+    keyword/value = optional keyword value pairs. supported keyword/value pairs:
+      {mult} n = dihedral multiplicity n (integer >= 0, default = 1) :pre
 :ule
 
 [Examples:]
@@ -155,11 +157,13 @@ associated with the restraint is
 with the following coefficients:
 
 K (energy)
-n = 1
+n (multiplicity, >= 0)
 d (degrees) = phi0 + 180 :ul
 
-K and phi0 are specified with the fix.  Note that the value of n is
-hard-wired to 1.  Also note that the energy will be a minimum when the
+K and phi0 are specified with the fix.  Note that the value of the
+dihedral multiplicity {n} is set by default to 1. You can use the
+optional {mult} keyword to set it to a different positive integer.
+Also note that the energy will be a minimum when the
 current dihedral angle phi is equal to phi0.
 
 :line

src/fix_restrain.cpp

@@ -45,7 +45,7 @@ enum{BOND,ANGLE,DIHEDRAL};
 
 FixRestrain::FixRestrain(LAMMPS *lmp, int narg, char **arg) :
   Fix(lmp, narg, arg),
-  rstyle(NULL), ids(NULL), kstart(NULL), kstop(NULL), target(NULL),
+  rstyle(NULL), mult(NULL), ids(NULL), kstart(NULL), kstop(NULL), target(NULL),
   cos_target(NULL), sin_target(NULL)
 {
   if (narg < 4) error->all(FLERR,"Illegal fix restrain command");
@@ -65,6 +65,7 @@ FixRestrain::FixRestrain(LAMMPS *lmp, int narg, char **arg) :
     if (nrestrain == maxrestrain) {
       maxrestrain += DELTA;
       memory->grow(rstyle,maxrestrain,"restrain:rstyle");
+      memory->grow(mult,maxrestrain,"restrain:mult");
       memory->grow(ids,maxrestrain,4,"restrain:ids");
       memory->grow(kstart,maxrestrain,"restrain:kstart");
       memory->grow(kstop,maxrestrain,"restrain:kstop");
@@ -96,6 +97,7 @@ FixRestrain::FixRestrain(LAMMPS *lmp, int narg, char **arg) :
     } else if (strcmp(arg[iarg],"dihedral") == 0) {
       if (iarg+8 > narg) error->all(FLERR,"Illegal fix restrain command");
       rstyle[nrestrain] = DIHEDRAL;
+      mult[nrestrain]   = 1;
       ids[nrestrain][0] = force->inumeric(FLERR,arg[iarg+1]);
       ids[nrestrain][1] = force->inumeric(FLERR,arg[iarg+2]);
       ids[nrestrain][2] = force->inumeric(FLERR,arg[iarg+3]);
@@ -107,6 +109,13 @@ FixRestrain::FixRestrain(LAMMPS *lmp, int narg, char **arg) :
       cos_target[nrestrain] = cos(target[nrestrain]);
       sin_target[nrestrain] = sin(target[nrestrain]);
       iarg += 8;
+      if ((iarg < narg) && (strcmp("mult",arg[iarg]) == 0)) {
+        if (iarg+1 > narg) error->all(FLERR,"Illegal fix restrain command");
+        mult[nrestrain] = force->inumeric(FLERR,arg[iarg+1]);
+        if (mult[nrestrain] < 0)
+          error->all(FLERR,"Illegal fix restrain command");
+        iarg += 2;
+      }
     } else error->all(FLERR,"Illegal fix restrain command");
 
     nrestrain++;
@@ -123,6 +132,7 @@ FixRestrain::FixRestrain(LAMMPS *lmp, int narg, char **arg) :
 FixRestrain::~FixRestrain()
 {
   memory->destroy(rstyle);
+  memory->destroy(mult);
   memory->destroy(ids);
   memory->destroy(kstart);
   memory->destroy(kstop);
@@ -410,7 +420,7 @@ void FixRestrain::restrain_angle(int m)
 
 void FixRestrain::restrain_dihedral(int m)
 {
-  int i1,i2,i3,i4,i,mult;
+  int i1,i2,i3,i4,i;
   double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm;
   double f1[3],f2[3],f3[3],f4[3];
   double ax,ay,az,bx,by,bz,rasq,rbsq,rgsq,rg,rginv,ra2inv,rb2inv,rabinv;
@@ -534,11 +544,10 @@ void FixRestrain::restrain_dihedral(int m)
   if (c > 1.0) c = 1.0;
   if (c < -1.0) c = -1.0;
 
-  mult = 1;  // multiplicity
   p = 1.0;
   df1 = 0.0;
 
-  for (i = 0; i < mult; i++) {
+  for (i = 0; i < mult[m]; i++) {
     ddf1 = p*c - df1*s;
     df1 = p*s + df1*c;
     p = ddf1;
@@ -546,7 +555,7 @@ void FixRestrain::restrain_dihedral(int m)
 
   p = p*cos_target[m] + df1*sin_target[m];
   df1 = df1*cos_target[m] - ddf1*sin_target[m];
-  df1 *= -mult;
+  df1 *= -mult[m];
   p += 1.0;
 
   energy += k * p;

src/fix_restrain.h

@@ -41,6 +41,7 @@ class FixRestrain : public Fix {
   int ilevel_respa;
   int nrestrain,maxrestrain;
   int *rstyle;
+  int *mult;
   int **ids;
   double *kstart,*kstop,*target;
   double *cos_target,*sin_target;