From cb466e44535ac45f2ee8341dcae39207b5813abe Mon Sep 17 00:00:00 2001 From: Justin Lemkul Date: Thu, 14 Aug 2014 19:04:34 -0400 Subject: [PATCH] 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 --- docs/manual/forcefield.tex | 9 ++-- src/gromacs/bonded/bonded.cpp | 67 ++++++++++++++++++++----- src/gromacs/gmxpreprocess/readir.c | 8 +++ src/gromacs/legacyheaders/types/enums.h | 2 +- 4 files changed, 69 insertions(+), 17 deletions(-) diff --git a/docs/manual/forcefield.tex b/docs/manual/forcefield.tex index f20986fc61..70829ad13e 100644 --- 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 a8757886d4..0ad3139599 100644 --- 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 b86bbef806..0322b5b28c 100644 --- 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 4b1161a4fa..8756469d3d 100644 --- 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 { -- 2.22.0