Converted pbcutils to C++
authorRossen Apostolov <rossen@kth.se>
Wed, 3 Sep 2014 10:36:48 +0000 (12:36 +0200)
committerGerrit Code Review <gerrit@gerrit.gromacs.org>
Sat, 12 Sep 2015 03:52:26 +0000 (05:52 +0200)
Used static_cast for double/float comparisons; removed unused
variables.
Replaced call to assert by GMX_ASSERT.
Replaced trjconv warning for failed put_atoms_in_compact_unitcell()
by a fatal_error in that function.
Added documentation.

Change-Id: I1e8c45d629b33974e655798b7b94edc545715f92

src/gromacs/gmxana/gmx_trjconv.cpp
src/gromacs/pbcutil/CMakeLists.txt
src/gromacs/pbcutil/pbc.cpp [moved from src/gromacs/pbcutil/pbc.c with 92% similarity]
src/gromacs/pbcutil/pbc.h
src/programs/view/manager.cpp

index 5e02c927c4b331e5dfcfc8a3a3ca9caaefefba26..4546240fe1a85148f7c7c5198ade2dc344304747 100644 (file)
@@ -905,8 +905,6 @@ int gmx_trjconv(int argc, char *argv[])
     char             *outf_base = NULL;
     const char       *outf_ext  = NULL;
     char              top_title[256], title[256], filemode[5];
-    gmx_bool          bWarnCompact = FALSE;
-    const char       *warn;
     output_env_t      oenv;
 
     t_filenm          fnm[] = {
@@ -1669,13 +1667,8 @@ int gmx_trjconv(int argc, char *argv[])
                                     put_atoms_in_triclinic_unitcell(ecenter, fr.box, natoms, fr.x);
                                     break;
                                 case euCompact:
-                                    warn = put_atoms_in_compact_unitcell(ePBC, ecenter, fr.box,
-                                                                         natoms, fr.x);
-                                    if (warn && !bWarnCompact)
-                                    {
-                                        fprintf(stderr, "\n%s\n", warn);
-                                        bWarnCompact = TRUE;
-                                    }
+                                    put_atoms_in_compact_unitcell(ePBC, ecenter, fr.box,
+                                                                  natoms, fr.x);
                                     break;
                             }
                         }
index 793ba194906e5f4409da15673ca1f23c25a07bb5..b87026387633b7956d02d2adaf85c9211a0c825d 100644 (file)
@@ -1,7 +1,7 @@
 #
 # This file is part of the GROMACS molecular simulation package.
 #
-# Copyright (c) 2014, by the GROMACS development team, led by
+# Copyright (c) 2014,2015, by the GROMACS development team, led by
 # Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
 # and including many others, as listed in the AUTHORS file in the
 # top-level source directory and at http://www.gromacs.org.
@@ -32,7 +32,7 @@
 # To help us fund GROMACS development, we humbly ask that you cite
 # the research papers on the package. Check out http://www.gromacs.org.
 
-file(GLOB PBCUTIL_SOURCES *.cpp *.c)
+file(GLOB PBCUTIL_SOURCES *.cpp)
 set(LIBGROMACS_SOURCES ${LIBGROMACS_SOURCES} ${PBCUTIL_SOURCES} PARENT_SCOPE)
 
 gmx_install_headers(
similarity index 92%
rename from src/gromacs/pbcutil/pbc.c
rename to src/gromacs/pbcutil/pbc.cpp
index c5ac5a2492cae647547efd6a9f155f62809a094f..efe150471e4cb987461c1dec6cb2892143b3551d 100644 (file)
  * To help us fund GROMACS development, we humbly ask that you cite
  * the research papers on the package. Check out http://www.gromacs.org.
  */
+/*! \internal \file
+ * \brief
+ * Implements routines in pbc.h.
+ *
+ * Utility functions for handling periodic boundary conditions.
+ * Mainly used in analysis tools.
+ */
 #include "gmxpre.h"
 
 #include "pbc.h"
 
-#include <assert.h>
-#include <math.h>
+#include <cmath>
+
+#include <algorithm>
 
 #include "gromacs/legacyheaders/gmx_omp_nthreads.h"
 #include "gromacs/legacyheaders/macros.h"
@@ -52,6 +60,7 @@
 #include "gromacs/pbcutil/ishift.h"
 #include "gromacs/pbcutil/mshift.h"
 #include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/gmxassert.h"
 #include "gromacs/utility/smalloc.h"
 
 /* Skip 0 so we have more chance of detecting if we forgot to call set_pbc. */
@@ -63,8 +72,9 @@ enum {
     epbcdxNOPBC,         epbcdxUNSUPPORTED
 };
 
-/* Margin factor for error message and correction if the box is too skewed */
+//! Margin factor for error message
 #define BOX_MARGIN         1.0010
+//! Margin correction if the box is too skewed
 #define BOX_MARGIN_CORRECT 1.0005
 
 int ePBC2npbcdim(int ePBC)
@@ -107,8 +117,6 @@ void dump_pbc(FILE *fp, t_pbc *pbc)
     rvec_add(pbc->hbox_diag, pbc->mhbox_diag, sum_box);
     pr_rvecs(fp, 0, "sum of the above two", &sum_box, 1);
     fprintf(fp, "max_cutoff2 = %g\n", pbc->max_cutoff2);
-    fprintf(fp, "bLimitDistance = %s\n", EBOOL(pbc->bLimitDistance));
-    fprintf(fp, "limit_distance2 = %g\n", pbc->limit_distance2);
     fprintf(fp, "ntric_vec = %d\n", pbc->ntric_vec);
     if (pbc->ntric_vec > 0)
     {
@@ -139,10 +147,10 @@ const char *check_box(int ePBC, matrix box)
     {
         ptr = "The unit cell can not have off-diagonal x-components with screw pbc";
     }
-    else if (fabs(box[YY][XX]) > BOX_MARGIN*0.5*box[XX][XX] ||
+    else if (std::fabs(box[YY][XX]) > BOX_MARGIN*0.5*box[XX][XX] ||
              (ePBC != epbcXY &&
-              (fabs(box[ZZ][XX]) > BOX_MARGIN*0.5*box[XX][XX] ||
-               fabs(box[ZZ][YY]) > BOX_MARGIN*0.5*box[YY][YY])))
+              (std::fabs(box[ZZ][XX]) > BOX_MARGIN*0.5*box[XX][XX] ||
+               std::fabs(box[ZZ][YY]) > BOX_MARGIN*0.5*box[YY][YY])))
     {
         ptr = "Triclinic box is too skewed.";
     }
@@ -156,15 +164,16 @@ const char *check_box(int ePBC, matrix box)
 
 real max_cutoff2(int ePBC, matrix box)
 {
-    real min_hv2, min_ss;
+    real       min_hv2, min_ss;
+    const real oneFourth = 0.25;
 
     /* Physical limitation of the cut-off
      * by half the length of the shortest box vector.
      */
-    min_hv2 = min(0.25*norm2(box[XX]), 0.25*norm2(box[YY]));
+    min_hv2 = oneFourth * std::min(norm2(box[XX]), norm2(box[YY]));
     if (ePBC != epbcXY)
     {
-        min_hv2 = min(min_hv2, 0.25*norm2(box[ZZ]));
+        min_hv2 = std::min(min_hv2, oneFourth * norm2(box[ZZ]));
     }
 
     /* Limitation to the smallest diagonal element due to optimizations:
@@ -174,17 +183,17 @@ real max_cutoff2(int ePBC, matrix box)
      */
     if (ePBC == epbcXY)
     {
-        min_ss = min(box[XX][XX], box[YY][YY]);
+        min_ss = std::min(box[XX][XX], box[YY][YY]);
     }
     else
     {
-        min_ss = min(box[XX][XX], min(box[YY][YY]-fabs(box[ZZ][YY]), box[ZZ][ZZ]));
+        min_ss = std::min(box[XX][XX], std::min(box[YY][YY] - std::fabs(box[ZZ][YY]), box[ZZ][ZZ]));
     }
 
-    return min(min_hv2, min_ss*min_ss);
+    return std::min(min_hv2, min_ss*min_ss);
 }
 
-/* this one is mostly harmless... */
+//! Set to true if warning has been printed
 static gmx_bool bWarnedGuess = FALSE;
 
 int guess_ePBC(matrix box)
@@ -223,6 +232,7 @@ int guess_ePBC(matrix box)
     return ePBC;
 }
 
+//! Check if the box still obeys the restrictions, if not, correct it
 static int correct_box_elem(FILE *fplog, int step, tensor box, int v, int d)
 {
     int shift, maxshift = 10;
@@ -279,7 +289,6 @@ gmx_bool correct_box(FILE *fplog, int step, tensor box, t_graph *graph)
     int      zy, zx, yx, i;
     gmx_bool bCorrected;
 
-    /* check if the box still obeys the restrictions, if not, correct it */
     zy = correct_box_elem(fplog, step, box, ZZ, YY);
     zx = correct_box_elem(fplog, step, box, ZZ, XX);
     yx = correct_box_elem(fplog, step, box, YY, XX);
@@ -323,6 +332,7 @@ int ndof_com(t_inputrec *ir)
     return n;
 }
 
+//! Do the real arithmetic for filling the pbc struct
 static void low_set_pbc(t_pbc *pbc, int ePBC, ivec *dd_nc, matrix box)
 {
     int         order[5] = {0, -1, 1, -2, 2};
@@ -330,14 +340,13 @@ static void low_set_pbc(t_pbc *pbc, int ePBC, ivec *dd_nc, matrix box)
     ivec        bPBC;
     real        d2old, d2new, d2new_c;
     rvec        trial, pos;
-    gmx_bool    bXY, bUse;
+    gmx_bool    bUse;
     const char *ptr;
 
     pbc->ePBC      = ePBC;
     pbc->ndim_ePBC = ePBC2npbcdim(ePBC);
 
     copy_mat(box, pbc->box);
-    pbc->bLimitDistance = FALSE;
     pbc->max_cutoff2    = 0;
     pbc->dim            = -1;
     pbc->ntric_vec      = 0;
@@ -358,14 +367,12 @@ static void low_set_pbc(t_pbc *pbc, int ePBC, ivec *dd_nc, matrix box)
     {
         fprintf(stderr,   "Warning: %s\n", ptr);
         pr_rvecs(stderr, 0, "         Box", box, DIM);
-        fprintf(stderr,   "         Can not fix pbc.\n");
-        pbc->ePBCDX          = epbcdxUNSUPPORTED;
-        pbc->bLimitDistance  = TRUE;
-        pbc->limit_distance2 = 0;
+        fprintf(stderr,   "         Can not fix pbc.\n\n");
+        pbc->ePBCDX = epbcdxUNSUPPORTED;
     }
     else
     {
-        if (ePBC == epbcSCREW && dd_nc)
+        if (ePBC == epbcSCREW && NULL != dd_nc)
         {
             /* This combinated should never appear here */
             gmx_incons("low_set_pbc called with screw pbc and dd_nc != NULL");
@@ -398,7 +405,7 @@ static void low_set_pbc(t_pbc *pbc, int ePBC, ivec *dd_nc, matrix box)
                         pbc->dim = i;
                     }
                 }
-                assert(pbc->dim < DIM);
+                GMX_ASSERT(pbc->dim < DIM, "Dimension for PBC incorrect");
                 for (i = 0; i < pbc->dim; i++)
                 {
                     if (pbc->box[pbc->dim][i] != 0)
@@ -514,11 +521,11 @@ static void low_set_pbc(t_pbc *pbc, int ePBC, ivec *dd_nc, matrix box)
                                 {
                                     if (trial[d] < 0)
                                     {
-                                        pos[d] = min( pbc->hbox_diag[d], -trial[d]);
+                                        pos[d] = std::min( pbc->hbox_diag[d], -trial[d]);
                                     }
                                     else
                                     {
-                                        pos[d] = max(-pbc->hbox_diag[d], -trial[d]);
+                                        pos[d] = std::max(-pbc->hbox_diag[d], -trial[d]);
                                     }
                                 }
                                 d2old += sqr(pos[d]);
@@ -530,6 +537,12 @@ static void low_set_pbc(t_pbc *pbc, int ePBC, ivec *dd_nc, matrix box)
                                 {
                                     /* Check if there is a single shift vector
                                      * that decreases this distance even more.
+                                     * This code can probably be removed, as
+                                     * our unit-cell restrictions checked
+                                     * before in check_box should avoid this.
+                                     * But because we have historically had
+                                     * some issues with triclinic boxes, we'll
+                                     * keep this check for now.
                                      */
                                     jc = 0;
                                     kc = 0;
@@ -547,17 +560,12 @@ static void low_set_pbc(t_pbc *pbc, int ePBC, ivec *dd_nc, matrix box)
                                         d2new_c += sqr(pos[d] + trial[d]
                                                        - jc*box[YY][d] - kc*box[ZZ][d]);
                                     }
-                                    if (d2new_c > BOX_MARGIN*d2new)
-                                    {
-                                        /* Reject this shift vector, as there is no a priori limit
-                                         * to the number of shifts that decrease distances.
-                                         */
-                                        if (!pbc->bLimitDistance || d2new <  pbc->limit_distance2)
-                                        {
-                                            pbc->limit_distance2 = d2new;
-                                        }
-                                        pbc->bLimitDistance = TRUE;
-                                    }
+                                    /* This should never happen, but if it does
+                                     * there might be no limit on the number
+                                     * of triclinic shift vectors.
+                                     */
+                                    GMX_RELEASE_ASSERT(d2new_c <= BOX_MARGIN*d2new, "Invalid triclinic box, but passed check_box()");
+                                    gmx_incons("ai");
                                 }
                                 else
                                 {
@@ -595,17 +603,18 @@ static void low_set_pbc(t_pbc *pbc, int ePBC, ivec *dd_nc, matrix box)
                                             pbc->tric_shift[pbc->ntric_vec][YY] = j;
                                             pbc->tric_shift[pbc->ntric_vec][ZZ] = k;
                                             pbc->ntric_vec++;
+
+                                            if (debug)
+                                            {
+                                                fprintf(debug, "  tricvec %2d = %2d %2d %2d  %5.2f %5.2f  %5.2f %5.2f %5.2f  %5.2f %5.2f %5.2f\n",
+                                                        pbc->ntric_vec, i, j, k,
+                                                        sqrt(d2old), sqrt(d2new),
+                                                        trial[XX], trial[YY], trial[ZZ],
+                                                        pos[XX], pos[YY], pos[ZZ]);
+                                            }
                                         }
                                     }
                                 }
-                                if (debug)
-                                {
-                                    fprintf(debug, "  tricvec %2d = %2d %2d %2d  %5.2f %5.2f  %5.2f %5.2f %5.2f  %5.2f %5.2f %5.2f\n",
-                                            pbc->ntric_vec, i, j, k,
-                                            sqrt(d2old), sqrt(d2new),
-                                            trial[XX], trial[YY], trial[ZZ],
-                                            pos[XX], pos[YY], pos[ZZ]);
-                                }
                             }
                         }
                     }
@@ -1114,6 +1123,7 @@ int pbc_dx_aiuc(const t_pbc *pbc, const rvec x1, const rvec x2, rvec dx)
     return is;
 }
 
+//! Compute distance vector in double precision
 void pbc_dx_d(const t_pbc *pbc, const dvec x1, const dvec x2, dvec dx)
 {
     int      i, j;
@@ -1237,6 +1247,7 @@ void pbc_dx_d(const t_pbc *pbc, const dvec x1, const dvec x2, dvec dx)
     }
 }
 
+//! Compute the box image corresponding to input vectors
 gmx_bool image_rect(ivec xi, ivec xj, ivec box_size, real rlong2, int *shift, real *r2)
 {
     int     m, t;
@@ -1415,8 +1426,9 @@ void calc_triclinic_images(matrix box, rvec img[])
 
 void calc_compact_unitcell_vertices(int ecenter, matrix box, rvec vert[])
 {
-    rvec img[NTRICIMG], box_center;
-    int  n, i, j, tmp[4], d;
+    rvec       img[NTRICIMG], box_center;
+    int        n, i, j, tmp[4], d;
+    const real oneFourth = 0.25;
 
     calc_triclinic_images(box, img);
 
@@ -1500,7 +1512,7 @@ void calc_compact_unitcell_vertices(int ecenter, matrix box, rvec vert[])
     {
         for (d = 0; d < DIM; d++)
         {
-            vert[i][d] = vert[i][d]*0.25+box_center[d];
+            vert[i][d] = vert[i][d]*oneFourth+box_center[d];
         }
     }
 }
@@ -1516,27 +1528,21 @@ int *compact_unitcell_edges()
         8, 20, 10, 18, 12, 16, 14, 22
     };
     int              e, i, j;
-    gmx_bool         bFirst = TRUE;
 
     snew(edge, NCUCEDGE*2);
 
-    if (bFirst)
+    e = 0;
+    for (i = 0; i < 6; i++)
     {
-        e = 0;
-        for (i = 0; i < 6; i++)
-        {
-            for (j = 0; j < 4; j++)
-            {
-                edge[e++] = 4*i + j;
-                edge[e++] = 4*i + (j+1) % 4;
-            }
-        }
-        for (i = 0; i < 12*2; i++)
+        for (j = 0; j < 4; j++)
         {
-            edge[e++] = hexcon[i];
+            edge[e++] = 4*i + j;
+            edge[e++] = 4*i + (j+1) % 4;
         }
-
-        bFirst = FALSE;
+    }
+    for (i = 0; i < 12*2; i++)
+    {
+        edge[e++] = hexcon[i];
     }
 
     return edge;
@@ -1676,23 +1682,24 @@ void put_atoms_in_triclinic_unitcell(int ecenter, matrix box,
     }
 }
 
-const char *
-put_atoms_in_compact_unitcell(int ePBC, int ecenter, matrix box,
-                              int natoms, rvec x[])
+void put_atoms_in_compact_unitcell(int ePBC, int ecenter, matrix box,
+                                   int natoms, rvec x[])
 {
     t_pbc pbc;
     rvec  box_center, dx;
     int   i;
 
     set_pbc(&pbc, ePBC, box);
+
+    if (pbc.ePBCDX == epbcdxUNSUPPORTED)
+    {
+        gmx_fatal(FARGS, "Can not put atoms in compact unitcell with unsupported PBC");
+    }
+
     calc_box_center(ecenter, box, box_center);
     for (i = 0; i < natoms; i++)
     {
         pbc_dx(&pbc, x[i], box_center, dx);
         rvec_add(box_center, dx, x[i]);
     }
-
-    return pbc.bLimitDistance ?
-           "WARNING: Could not put all atoms in the compact unitcell\n"
-           : NULL;
 }
index 1fdd756bba2296902c468221aa67c2f2d9ffe2e2..aeefec137eac778e4739ae9ae5cbdb4e64cfda87 100644 (file)
@@ -55,20 +55,49 @@ extern "C" {
  */
 #define MAX_NTRICVEC 12
 
+/*! \brief Structure containing info on periodic boundary conditions */
 typedef struct t_pbc {
+    //! The PBC type
     int        ePBC;
+    //! Number of dimensions in which PBC is exerted
     int        ndim_ePBC;
+    /*! \brief Determines how to compute distance vectors.
+     *
+     *  Indicator of how to compute distance vectors, depending
+     *  on PBC type (depends on ePBC and dimensions with(out) DD)
+     *  and the box angles.
+     */
     int        ePBCDX;
+    /*! \brief Used for selecting which dimensions to use in PBC.
+     *
+     *  In case of 1-D PBC this indicates which dimension is used,
+     *  in case of 2-D PBC this indicates the opposite
+     */
     int        dim;
+    //! The simulation box
     matrix     box;
+    //! The lengths of the diagonal of the full box
     rvec       fbox_diag;
+    //! Halve of the above
     rvec       hbox_diag;
+    //! Negative of the above
     rvec       mhbox_diag;
+    //! Maximum allowed cutoff squared for the box and PBC used
     real       max_cutoff2;
-    gmx_bool   bLimitDistance;
-    real       limit_distance2;
+    /*! \brief Number of triclinic shift vectors.
+     *
+     *  Number of triclinic shift vectors depends on the skewedness
+     *  of the box, that is mostly on the angles. For triclinic boxes
+     *  we first take the closest image along each Cartesian dimension
+     *  independently. When the resulting distance^2 is larger than
+     *  max_cutoff2, up to ntric_vec triclinic shift vectors need to
+     *  be tried. Because of the restrictions imposed on the unit-cell
+     *  by GROMACS, ntric_vec <= MAX_NTRICVEC = 12.
+     */
     int        ntric_vec;
+    //! The triclinic shift vectors in grid cells. Internal use only.
     ivec       tric_shift[MAX_NTRICVEC];
+    //!  The triclinic shift vectors in length units
     rvec       tric_vec[MAX_NTRICVEC];
 } t_pbc;
 
@@ -87,168 +116,284 @@ enum {
 
 struct t_graph;
 
+/*! \brief Returns the number of dimensions that use pbc
+ *
+ * \param[in] ePBC The periodic boundary condition type
+ * \return the number of dimensions that use pbc, starting at X
+ */
 int ePBC2npbcdim(int ePBC);
-/* Returns the number of dimensions that use pbc, starting at X */
 
-int inputrec2nboundeddim(const t_inputrec *ir);
-/* Returns the number of dimensions in which
+/*! \brief Return the number of bounded directories
+ *
+ * \param[in] ir The input record with MD parameters
+ * \return the number of dimensions in which
  * the coordinates of the particles are bounded, starting at X.
  */
+int inputrec2nboundeddim(const t_inputrec *ir);
 
+/*! \brief Dump the contents of the pbc structure to the file
+ *
+ * \param[in] fp  The file pointer to write to
+ * \param[in] pbc The periodic boundary condition information structure
+ */
 void dump_pbc(FILE *fp, t_pbc *pbc);
-/* Dump the contents of the pbc structure to the file */
 
-const char *check_box(int ePBC, matrix box);
-/* Returns NULL if the box is supported by Gromacs.
- * Otherwise is returns a string with the problem.
+/*! \brief Check the box for consistency
+ *
+ * \param[in] ePBC The pbc identifier
+ * \param[in] box  The box matrix
+ * \return NULL if the box is supported by Gromacs.
+ * Otherwise returns a string with the problem.
  * When ePBC=-1, the type of pbc is guessed from the box matrix.
  */
+const char *check_box(int ePBC, matrix box);
 
-real max_cutoff2(int ePBC, matrix box);
-/* Returns the square of the maximum cut-off allowed for the box,
+/*! \brief Compute the maximum cutoff for the box
+
+ * Returns the square of the maximum cut-off allowed for the box,
  * taking into account that the grid neighborsearch code and pbc_dx
  * only check combinations of single box-vector shifts.
+ * \param[in] ePBC The pbc identifier
+ * \param[in] box  The box matrix
+ * \return the maximum cut-off.
  */
+real max_cutoff2(int ePBC, matrix box);
 
+/*! \brief Guess PBC typr
+ *
+ * Guesses the type of periodic boundary conditions using the box
+ * \param[in] box  The box matrix
+ * \return The pbc identifier
+ */
 int guess_ePBC(matrix box);
-/* Guesses the type of periodic boundary conditions using the box */
 
-gmx_bool correct_box(FILE *fplog, int step, tensor box, struct t_graph *graph);
-/* Checks for un-allowed box angles and corrects the box
+/*! \brief Corrects the box if necessary
+ *
+ * Checks for un-allowed box angles and corrects the box
  * and the integer shift vectors in the graph (if graph!=NULL) if necessary.
- * Returns TRUE when the box was corrected.
+ * \param[in] fplog File for debug output
+ * \param[in] step  The MD step number
+ * \param[in] box   The simulation cell
+ * \param[in] graph Information about molecular connectivity
+ * \return TRUE when the box was corrected.
  */
+gmx_bool correct_box(FILE *fplog, int step, tensor box, struct t_graph *graph);
 
+/*! \brief Returns the number of degrees of freedom in center of mass motion
+ *
+ * \param[in] ir the inputrec structure
+ * \return the number of degrees of freedom of the center of mass
+ */
 int ndof_com(t_inputrec *ir);
-/* Returns the number of degrees of freedom of the center of mass */
 
-void set_pbc(t_pbc *pbc, int ePBC, matrix box);
-/* Initiate the periodic boundary conditions.
+/*! \brief Initiate the periodic boundary condition algorithms.
+ *
  * pbc_dx will not use pbc and return the normal difference vector
  * when one or more of the diagonal elements of box are zero.
  * When ePBC=-1, the type of pbc is guessed from the box matrix.
+ * \param[inout] pbc The pbc information structure
+ * \param[in] ePBC The PBC identifier
+ * \param[in] box  The box tensor
  */
+void set_pbc(t_pbc *pbc, int ePBC, matrix box);
 
-t_pbc *set_pbc_dd(t_pbc *pbc, int ePBC,
-                  gmx_domdec_t *dd, gmx_bool bSingleDir, matrix box);
-/* As set_pbc, but additionally sets that correct distances can
+/*! \brief Initiate the periodic boundary condition algorithms.
+ *
+ * As set_pbc, but additionally sets that correct distances can
  * be obtained using (combinations of) single box-vector shifts.
  * Should be used with pbc_dx_aiuc.
  * If dd!=NULL pbc is not used for directions
  * with dd->nc[i]==1 with bSingleDir==TRUE or
  * with dd->nc[i]<=2 with bSingleDir==FALSE.
- * Returns pbc when pbc operations are required, NULL otherwise.
+ * \param[inout] pbc The pbc information structure
+ * \param[in] ePBC       The PBC identifier
+ * \param[in] dd         The domain decomposition struct
+ * \param[in] bSingleDir TRUE if DD communicates only in one direction along dimensions
+ * \param[in] box        The box tensor
+ * \return the pbc structure when pbc operations are required, NULL otherwise.
  */
+t_pbc *set_pbc_dd(t_pbc *pbc, int ePBC,
+                  gmx_domdec_t *dd, gmx_bool bSingleDir, matrix box);
 
-void pbc_dx(const t_pbc *pbc, const rvec x1, const rvec x2, rvec dx);
-/* Calculate the correct distance vector from x2 to x1 and put it in dx.
+/*! \brief Compute distance with PBC
+ *
+ * Calculate the correct distance vector from x2 to x1 and put it in dx.
  * set_pbc must be called before ever calling this routine.
  *
- * For triclinic boxes pbc_dx does not necessarily return the shortest
- * distance vector. If pbc->bLimitDistance=TRUE an atom pair with
- * distance vector dx with norm2(dx) > pbc->limit_distance2 could
- * have a shorter distance, but not shorter than sqrt(pbc->limit_distance2).
- * pbc->limit_distance2 is always larger than max_cutoff2(box).
- * For the standard rhombic dodecahedron and truncated octahedron
- * pbc->bLimitDistance=FALSE and thus all distances are correct.
+ * Note that for triclinic boxes that do not obey the GROMACS unit-cell
+ * restrictions, pbc_dx and pbc_dx_aiuc will not correct for PBC.
+ * \param[inout] pbc The pbc information structure
+ * \param[in]    x1  Coordinates for particle 1
+ * \param[in]    x2  Coordinates for particle 2
+ * \param[out]   dx  Distance vector
  */
+void pbc_dx(const t_pbc *pbc, const rvec x1, const rvec x2, rvec dx);
 
-int pbc_dx_aiuc(const t_pbc *pbc, const rvec x1, const rvec x2, rvec dx);
-/* Calculate the correct distance vector from x2 to x1 and put it in dx,
+/*! \brief Compute distance vector for simple PBC types
+ *
+ * Calculate the correct distance vector from x2 to x1 and put it in dx,
  * This function can only be used when all atoms are in the rectangular
  * or triclinic unit-cell.
- * Returns the ishift required to shift x1 at closest distance to x2;
+ * set_pbc_dd or set_pbc must be called before ever calling this routine.
+ * \param[inout] pbc The pbc information structure
+ * \param[in]    x1  Coordinates for particle 1
+ * \param[in]    x2  Coordinates for particle 2
+ * \param[out]   dx  Distance vector
+ * \return the ishift required to shift x1 at closest distance to x2;
  * i.e. if 0<=ishift<SHIFTS then x1 - x2 + shift_vec[ishift] = dx
  * (see calc_shifts below on how to obtain shift_vec)
- * set_pbc_dd or set_pbc must be called before ever calling this routine.
  */
-void pbc_dx_d(const t_pbc *pbc, const dvec x1, const dvec x2, dvec dx);
-/* As pbc_dx, but for double precision vectors.
+int pbc_dx_aiuc(const t_pbc *pbc, const rvec x1, const rvec x2, rvec dx);
+
+/*\brief Compute distance with PBC
+ *
+ * As pbc_dx, but for double precision vectors.
  * set_pbc must be called before ever calling this routine.
+ * \param[inout] pbc The pbc information structure
+ * \param[in]    x1  Coordinates for particle 1
+ * \param[in]    x2  Coordinates for particle 2
+ * \param[out]   dx  Distance vector
  */
+void pbc_dx_d(const t_pbc *pbc, const dvec x1, const dvec x2, dvec dx);
 
-gmx_bool image_rect(ivec xi, ivec xj, ivec box_size,
-                    real rlong2, int *shift, real *r2);
-/* Calculate the distance between xi and xj for a rectangular box.
- * When the distance is SMALLER than rlong2 return TRUE, return
- * the shift code in shift and the distance in r2. When the distance is
- * >= rlong2 return FALSE;
+/*! \brief Calculate the distance between xi and xj for a rectangular box.
+ *
  * It is assumed that rlong2 is scaled the same way as the ivecs xi and xj.
+ * \param[in]  xi     Box index
+ * \param[in]  xj     Box index
+ * \param[in]  box    The box of grid cells
+ * \param[in]  rlong2 Cutoff squared
+ * \param[out] shift  The shift code
+ * \param[out] r2     The distance (squared???)
+ * \return TRUE when the distance is SMALLER than rlong2
  */
+gmx_bool image_rect(ivec xi, ivec xj, imatrix box,
+                    real rlong2, int *shift, real *r2);
 
-gmx_bool image_tri(ivec xi, ivec xj, imatrix box,
-                   real rlong2, int *shift, real *r2);
-/* Calculate the distance between xi and xj for a triclinic box.
- * When the distance is SMALLER than rlong2 return TRUE, return
- * the shift code in shift and the distance in r2. When the distance is
- * >= rlong2 return FALSE;
+/*! \brief Calculate the distance between xi and xj for a triclinic box.
+ *
  * It is assumed that rlong2 is scaled the same way as the ivecs xi and xj.
+ * \param[in]  xi     Box index
+ * \param[in]  xj     Box index
+ * \param[in]  box    Matrix of box grid cells
+ * \param[in]  rlong2 Cutoff squared
+ * \param[out] shift  The shift code
+ * \param[out] r2     The distance (squared???)
+ * \return TRUE when the distance is SMALLER than rlong2
  */
+gmx_bool image_tri(ivec xi, ivec xj, imatrix box,
+                   real rlong2, int *shift, real *r2);
 
-gmx_bool image_cylindric(ivec xi, ivec xj, ivec box_size, real rlong2,
-                         int *shift, real *r2);
-/* Calculate the distance between xi and xj for a rectangular box
+/*! \brief Compute distance vector when using cylindrical cutoff
+ *
+ * Calculate the distance between xi and xj for a rectangular box
  * using a cylindric cutoff for long-range only.
- * When the distance is SMALLER than rlong2 (in X and Y dir.)
- * return TRUE, return
- * the shift code in shift and the distance in r2. When the distance is
- * >= rlong2 return FALSE;
  * It is assumed that rlong2 is scaled the same way as the ivecs xi and xj.
+ * \param[in]  xi       Box index
+ * \param[in]  xj       Box index
+ * \param[in]  box_size Number of box grid cells
+ * \param[in]  rlong2   Cutoff squared
+ * \param[out] shift    The shift code
+ * \param[out] r2       The distance (squared???)
+ * \return TRUE when the distance is SMALLER than rlong2 (in X and Y dir)
  */
+gmx_bool image_cylindric(ivec xi, ivec xj, ivec box_size, real rlong2,
+                         int *shift, real *r2);
 
-void calc_shifts(matrix box, rvec shift_vec[]);
-/* This routine calculates ths shift vectors necessary to use the
- * ns routine.
+/*! \brief Computes shift vectors
+ *
+ * This routine calculates ths shift vectors necessary to use the
+ * neighbor searching routine.
+ * \param[in]  box       The simulation box
+ * \param[out] shift_vec The shifting vectors
  */
+void calc_shifts(matrix box, rvec shift_vec[]);
 
-void calc_box_center(int ecenter, matrix box, rvec box_center);
-/* Calculates the center of the box.
+/*! \brief Calculates the center of the box.
+ *
  * See the description for the enum ecenter above.
+ * \param[in]  ecenter    Description of center type
+ * \param[in]  box        The simulation box
+ * \param[out] box_center The center of the box
  */
+void calc_box_center(int ecenter, matrix box, rvec box_center);
 
+/*! \brief Calculates the NTRICIMG box images
+ *
+ * \param[in]  box The simulation box
+ * \param[out] img The triclinic box images
+ */
 void calc_triclinic_images(matrix box, rvec img[]);
-/* Calculates the NTRICIMG box images */
 
+/*! \brief Calculates the NCUCVERT vertices of a compact unitcell
+ *
+ * \param[in]  ecenter The center type
+ * \param[in]  box     The simulation box
+ * \param[out] vert    The vertices
+ */
 void calc_compact_unitcell_vertices(int ecenter, matrix box,
                                     rvec vert[]);
-/* Calculates the NCUCVERT vertices of a compact unitcell */
 
-int *compact_unitcell_edges(void);
-/* Return an array of unitcell edges of length NCUCEDGE*2,
+/*! \brief Compute unitcell edges
+ *
+ * \return an array of unitcell edges of length NCUCEDGE*2,
  * this is an index in vert[], which is calculated by calc_unitcell_vertices.
  * The index consists of NCUCEDGE pairs of vertex indices.
  * The index does not change, so it needs to be retrieved only once.
  */
+int *compact_unitcell_edges(void);
 
-void put_atoms_in_box_omp(int ePBC, matrix box, int natoms, rvec x[]);
-/* This wrapper function around put_atoms_in_box() with the ugly manual
- * workload splitting is needed toavoid silently introducing multithreading
+/*! \brief Parallellizes put_atoms_box
+ *
+ * This wrapper function around put_atoms_in_box() with the ugly manual
+ * workload splitting is needed to avoid silently introducing multithreading
  * in tools.
- * */
+ * \param[in]    ePBC   The pbc type
+ * \param[in]    box    The simulation box
+ * \param[in]    natoms The number of atoms
+ * \param[inout] x      The coordinates of the atoms
+ */
+void put_atoms_in_box_omp(int ePBC, matrix box, int natoms, rvec x[]);
 
 
-void put_atoms_in_box(int ePBC, matrix box, int natoms, rvec x[]);
-/* These routines puts ONE or ALL atoms in the box, not caring
+/*! \brief Put atoms inside the simulations box
+ *
+ * These routines puts ONE or ALL atoms in the box, not caring
  * about charge groups!
  * Also works for triclinic cells.
+ * \param[in]    ePBC   The pbc type
+ * \param[in]    box    The simulation box
+ * \param[in]    natoms The number of atoms
+ * \param[inout] x      The coordinates of the atoms
  */
+void put_atoms_in_box(int ePBC, matrix box, int natoms, rvec x[]);
 
-void put_atoms_in_triclinic_unitcell(int ecenter, matrix box,
-                                     int natoms, rvec x[]);
-/* This puts ALL atoms in the triclinic unit cell, centered around the
+/*! \brief Put atoms inside triclinic box
+ *
+ * This puts ALL atoms in the triclinic unit cell, centered around the
  * box center as calculated by calc_box_center.
+ * \param[in]    ecenter The pbc center type
+ * \param[in]    box     The simulation box
+ * \param[in]    natoms  The number of atoms
+ * \param[inout] x       The coordinates of the atoms
  */
+void put_atoms_in_triclinic_unitcell(int ecenter, matrix box,
+                                     int natoms, rvec x[]);
 
-const char *put_atoms_in_compact_unitcell(int ePBC, int ecenter,
-                                          matrix box,
-                                          int natoms, rvec x[]);
-/* This puts ALL atoms at the closest distance for the center of the box
+/*! \brief Put atoms inside the unitcell
+ *
+ * This puts ALL atoms at the closest distance for the center of the box
  * as calculated by calc_box_center.
- * Will return NULL is everything went ok and a warning string if not
- * all atoms could be placed in the unitcell. This can happen for some
- * triclinic unitcells, see the comment at pbc_dx above.
  * When ePBC=-1, the type of pbc is guessed from the box matrix.
+ * \param[in]    ePBC    The pbc type
+ * \param[in]    ecenter The pbc center type
+ * \param[in]    box     The simulation box
+ * \param[in]    natoms  The number of atoms
+ * \param[inout] x       The coordinates of the atoms
  */
+void put_atoms_in_compact_unitcell(int ePBC, int ecenter,
+                                   matrix box,
+                                   int natoms, rvec x[]);
 
 #ifdef __cplusplus
 }
index cfa9ca3932502cfcfbb5a7798f18dfacb8ea28b3..a2538f883f1e2e6d030e909f362f1c5119ce9991 100644 (file)
@@ -350,9 +350,7 @@ static void reset_mols(t_block *mols, matrix box, rvec x[])
 static bool step_man(t_manager *man, int *nat)
 {
     static int      ncount = 0;
-    static bool     bWarn  = false;
     bool            bEof;
-    const char     *warn;
 
     if (!man->natom)
     {
@@ -369,13 +367,8 @@ static bool step_man(t_manager *man, int *nat)
                 put_atoms_in_triclinic_unitcell(ecenterDEF, man->box, man->natom, man->x);
                 break;
             case esbTrunc:
-                warn = put_atoms_in_compact_unitcell(man->molw->ePBC, ecenterDEF, man->box,
-                                                     man->natom, man->x);
-                if (warn && !bWarn)
-                {
-                    fprintf(stderr, "\n%s\n", warn);
-                    bWarn = true;
-                }
+                put_atoms_in_compact_unitcell(man->molw->ePBC, ecenterDEF, man->box,
+                                              man->natom, man->x);
                 break;
             case esbRect:
             case esbNone: