Added flat-bottom cylindrical restraints in x and y.

The code previously only allowed the long axis of the cylinder
to be aligned along z.  This patch allows the cylinder to be
aligned along any axis.

Old topologies and .tpr files will still work for all geometries with
g values from 2 to 5, as 2 is now just a synonym for a z-directional
cylinder.  The new values are g=6 for x, g=7 for y, and g=8 for z.
The manual has been updated to reflect these changes.

Introduced a new static function that computes forces and energy for the
flat-bottom restraint.

Refs #1577

Change-Id: I40e0bd35225170e13a18ee54f6798470e697f29d

diff --git a/docs/manual/forcefield.tex b/docs/manual/forcefield.tex
index f20986fc617ff7547f5e8cb5dad0fd3a0ff52bb2..70829ad13e07568f4187669fcc171ec40c4bbb55 100644 (file)
--- a/docs/manual/forcefield.tex
+++ b/docs/manual/forcefield.tex
@@ -1087,9 +1087,12 @@ radius. The force acts towards the center of the sphere. The following distance
 \beq
   d_g(\ve{r}_i;\ve{R}_i) = |\ve{r}_i-\ve{R}_i|
 \eeq
-{\bfseries Cylinder} ($g=2$): The particle is kept in a cylinder of given radius
-parallel to the $z$-axis. The force from the flat-bottomed potential acts
-towards the axis of the cylinder. The $z$-component of the force is zero.
+{\bfseries Cylinder} ($g=6,7,8$): The particle is kept in a cylinder of given radius
+parallel to the $x$ ($g=6$), $y$ ($g=7$), or $z$-axis ($g=8$). For backwards compatibility, setting
+$g=2$ is mapped to $g=8$ in the code so that old {\tt .tpr} files and topologies work.  
+The force from the flat-bottomed potential acts towards the axis of the cylinder. 
+The component of the force parallel to the cylinder axis is zero.
+For a cylinder aligned along the $z$-axis:
 \beq
  d_g(\ve{r}_i;\ve{R}_i) = \sqrt{ (x_i-X_i)^2 + (y_i - Y_i)^2 }
 \eeq

diff --git a/src/gromacs/bonded/bonded.cpp b/src/gromacs/bonded/bonded.cpp
index a8757886d410cf35c224baf0894e719fa840b51f..0ad3139599b1f8aa49333ea3e932247626d4e381 100644 (file)
--- a/src/gromacs/bonded/bonded.cpp
+++ b/src/gromacs/bonded/bonded.cpp
@@ -2255,6 +2255,44 @@ static void posres_dx(const rvec x, const rvec pos0A, const rvec pos0B,
     }
 }
 
+/*! \brief Computes forces and potential for flat-bottom cylindrical restraints.
+ *         Returns the flat-bottom potential. */
+static real do_fbposres_cylinder(int fbdim, rvec fm, rvec dx, real rfb, real kk, gmx_bool bInvert)
+{
+    int     d;
+    real    dr, dr2, invdr, v, rfb2;
+
+    dr2  = 0.0;
+    rfb2 = sqr(rfb);
+    v    = 0.0;
+
+    for (d = 0; d < DIM; d++)
+    {
+        if (d != fbdim)
+        {
+            dr2 += sqr(dx[d]);
+        }
+    }
+
+    if  (dr2 > 0.0 &&
+         ( (dr2 > rfb2 && bInvert == FALSE ) || (dr2 < rfb2 && bInvert == TRUE ) )
+         )
+    {
+        dr     = sqrt(dr2);
+        invdr  = 1./dr;
+        v      = 0.5*kk*sqr(dr - rfb);
+        for (d = 0; d < DIM; d++)
+        {
+            if (d != fbdim)
+            {
+                fm[d] = -kk*(dr-rfb)*dx[d]*invdr; /* Force pointing to the center */
+            }
+        }
+    }
+
+    return v;
+}
+
 /*! \brief Adds forces of flat-bottomed positions restraints to f[]
  *         and fixes vir_diag. Returns the flat-bottomed potential. */
 real fbposres(int nbonds,
@@ -2267,7 +2305,7 @@ real fbposres(int nbonds,
     int              i, ai, m, d, type, npbcdim = 0, fbdim;
     const t_iparams *pr;
     real             vtot, kk, v;
-    real             dr, dr2, rfb, rfb2, fact, invdr;
+    real             dr, dr2, rfb, rfb2, fact;
     rvec             com_sc, rdist, dx, dpdl, fm;
     gmx_bool         bInvert;
 
@@ -2329,19 +2367,22 @@ real fbposres(int nbonds,
                     svmul(fact, dx, fm);
                 }
                 break;
+            case efbposresCYLINDERX:
+                /* cylindrical flat-bottom posres in y-z plane. fm[XX] = 0. */
+                fbdim = XX;
+                v     = do_fbposres_cylinder(fbdim, fm, dx, rfb, kk, bInvert);
+                break;
+            case efbposresCYLINDERY:
+                /* cylindrical flat-bottom posres in x-z plane. fm[YY] = 0. */
+                fbdim = YY;
+                v     = do_fbposres_cylinder(fbdim, fm, dx, rfb, kk, bInvert);
+                break;
             case efbposresCYLINDER:
-                /* cylidrical flat-bottom posres in x-y plane. fm[ZZ] = 0. */
-                dr2 = sqr(dx[XX])+sqr(dx[YY]);
-                if  (dr2 > 0.0 &&
-                     ( (dr2 > rfb2 && bInvert == FALSE ) || (dr2 < rfb2 && bInvert == TRUE ) )
-                     )
-                {
-                    dr     = sqrt(dr2);
-                    invdr  = 1./dr;
-                    v      = 0.5*kk*sqr(dr - rfb);
-                    fm[XX] = -kk*(dr-rfb)*dx[XX]*invdr; /* Force pointing to the center */
-                    fm[YY] = -kk*(dr-rfb)*dx[YY]*invdr;
-                }
+            /* equivalent to efbposresCYLINDERZ for backwards compatibility */
+            case efbposresCYLINDERZ:
+                /* cylindrical flat-bottom posres in x-y plane. fm[ZZ] = 0. */
+                fbdim = ZZ;
+                v     = do_fbposres_cylinder(fbdim, fm, dx, rfb, kk, bInvert);
                 break;
             case efbposresX: /* fbdim=XX */
             case efbposresY: /* fbdim=YY */

diff --git a/src/gromacs/gmxpreprocess/readir.c b/src/gromacs/gmxpreprocess/readir.c
index b86bbef806e7f0a6317fe8900035573c4612b2d3..0322b5b28c490ee4807eb4a3d05f38fd8234d28e 100644 (file)
--- a/src/gromacs/gmxpreprocess/readir.c
+++ b/src/gromacs/gmxpreprocess/readir.c
@@ -3823,7 +3823,15 @@ static gmx_bool absolute_reference(t_inputrec *ir, gmx_mtop_t *sys,
                         case efbposresSPHERE:
                             AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1;
                             break;
+                        case efbposresCYLINDERX:
+                            AbsRef[YY] = AbsRef[ZZ] = 1;
+                            break;
+                        case efbposresCYLINDERY:
+                            AbsRef[XX] = AbsRef[ZZ] = 1;
+                            break;
                         case efbposresCYLINDER:
+                        /* efbposres is a synonym for efbposresCYLINDERZ for backwards compatibility */
+                        case efbposresCYLINDERZ:
                             AbsRef[XX] = AbsRef[YY] = 1;
                             break;
                         case efbposresX: /* d=XX */

diff --git a/src/gromacs/legacyheaders/types/enums.h b/src/gromacs/legacyheaders/types/enums.h
index 4b1161a4fa7d10e1b68a355778a64e9fd8c2ec65..8756469d3dd5a19057e7b2acb2ca4e33b30f2216 100644 (file)
--- a/src/gromacs/legacyheaders/types/enums.h
+++ b/src/gromacs/legacyheaders/types/enums.h
@@ -346,7 +346,7 @@ enum {
 /* flat-bottom posres geometries */
 enum {
     efbposresZERO, efbposresSPHERE, efbposresCYLINDER, efbposresX, efbposresY, efbposresZ,
-    efbposresNR
+    efbposresCYLINDERX, efbposresCYLINDERY, efbposresCYLINDERZ, efbposresNR
 };
 
 enum {