Sort all includes in src/gromacs

[alexxy/gromacs.git] / src / gromacs / bonded / bonded.cpp

diff --git a/src/gromacs/bonded/bonded.cpp b/src/gromacs/bonded/bonded.cpp
index a8757886d410cf35c224baf0894e719fa840b51f..481c34ae8a75ca9e1a305112dbea8582e8b0e5fa 100644 (file)
--- a/src/gromacs/bonded/bonded.cpp
+++ b/src/gromacs/bonded/bonded.cpp
@@ -41,22 +41,23 @@
 #include "config.h"
 
 #include <assert.h>
+
 #include <cmath>
 
 #include <algorithm>
 
-#include "gromacs/math/units.h"
-#include "gromacs/math/vec.h"
-#include "gromacs/math/utilities.h"
-#include "gromacs/legacyheaders/txtdump.h"
-#include "gromacs/legacyheaders/ns.h"
-#include "gromacs/legacyheaders/macros.h"
-#include "gromacs/legacyheaders/names.h"
+#include "gromacs/bonded/restcbt.h"
 #include "gromacs/legacyheaders/disre.h"
-#include "gromacs/legacyheaders/orires.h"
 #include "gromacs/legacyheaders/force.h"
+#include "gromacs/legacyheaders/macros.h"
+#include "gromacs/legacyheaders/names.h"
 #include "gromacs/legacyheaders/nonbonded.h"
-
+#include "gromacs/legacyheaders/ns.h"
+#include "gromacs/legacyheaders/orires.h"
+#include "gromacs/legacyheaders/txtdump.h"
+#include "gromacs/math/units.h"
+#include "gromacs/math/utilities.h"
+#include "gromacs/math/vec.h"
 #include "gromacs/pbcutil/ishift.h"
 #include "gromacs/pbcutil/mshift.h"
 #include "gromacs/pbcutil/pbc.h"
@@ -66,8 +67,6 @@
 #include "gromacs/utility/fatalerror.h"
 #include "gromacs/utility/smalloc.h"
 
-#include "restcbt.h"
-
 /* Find a better place for this? */
 const int cmap_coeff_matrix[] = {
     1, 0, -3,  2, 0, 0,  0,  0, -3,  0,  9, -6,  2,  0, -6,  4,
@@ -2255,6 +2254,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 +2304,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 +2366,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 */