diff --git a/doc/src/fix_restrain.txt b/doc/src/fix_restrain.txt index c8ec20daaa..1a7c7b6ba5 100644 --- a/doc/src/fix_restrain.txt +++ b/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 diff --git a/src/fix_restrain.cpp b/src/fix_restrain.cpp index 0f6c86ba5b..09cb8946e5 100644 --- a/src/fix_restrain.cpp +++ b/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; diff --git a/src/fix_restrain.h b/src/fix_restrain.h index a6f23f86c6..14699ed5af 100644 --- a/src/fix_restrain.h +++ b/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;