Move old neighborsearch code to C++
authorErik Lindahl <erik@kth.se>
Tue, 7 Jul 2015 20:02:24 +0000 (22:02 +0200)
committerGerrit Code Review <gerrit@gerrit.gromacs.org>
Wed, 8 Jul 2015 07:14:59 +0000 (09:14 +0200)
This will eventually be removed in favor of the new
verlet code, but for now it will make life simpler to
have everything moved to C++.

Change-Id: I3242e5f345769a42ffb99d1953bc1bef39324ff5

src/gromacs/mdlib/ns.cpp [moved from src/gromacs/mdlib/ns.c with 97% similarity]
src/gromacs/mdlib/nsgrid.cpp [moved from src/gromacs/mdlib/nsgrid.c with 97% similarity]

similarity index 97%
rename from src/gromacs/mdlib/ns.c
rename to src/gromacs/mdlib/ns.cpp
index d6e9314fb16a70eeca45b33d9fb3774d8463a34c..e3f37fc45a6b089a0ae6715d3b82b81f0bc98023 100644 (file)
 #include <stdlib.h>
 #include <string.h>
 
+#include <cmath>
+
+#include <algorithm>
+
 #include "gromacs/domdec/domdec.h"
 #include "gromacs/legacyheaders/force.h"
 #include "gromacs/legacyheaders/macros.h"
@@ -93,7 +97,7 @@ static gmx_bool NOTEXCL_(t_excl e[], atom_id i, atom_id j)
 static int
 round_up_to_simd_width(int length, int simd_width)
 {
-    int offset, newlength;
+    int offset;
 
     offset = (simd_width > 0) ? length % simd_width : 0;
 
@@ -132,7 +136,7 @@ static void init_nblist(FILE *log, t_nblist *nl_sr, t_nblist *nl_lr,
 {
     t_nblist *nl;
     int       homenr;
-    int       i, nn;
+    int       i;
 
     for (i = 0; (i < 2); i++)
     {
@@ -199,7 +203,6 @@ void init_neighbor_list(FILE *log, t_forcerec *fr, int homenr)
      */
     int        maxsr, maxsr_wat, maxlr, maxlr_wat;
     int        ielec, ivdw, ielecmod, ivdwmod, type;
-    int        solvent;
     int        igeometry_def, igeometry_w, igeometry_ww;
     int        i;
     gmx_bool   bElecAndVdwSwitchDiffers;
@@ -217,11 +220,11 @@ void init_neighbor_list(FILE *log, t_forcerec *fr, int homenr)
      * all the nlist arrays many times in a row.
      * The numbers seem very accurate, but they are uncritical.
      */
-    maxsr_wat = min(fr->nWatMol, (homenr+2)/3);
+    maxsr_wat = std::min(fr->nWatMol, (homenr+2)/3);
     if (fr->bTwinRange)
     {
         maxlr     = 50;
-        maxlr_wat = min(maxsr_wat, maxlr);
+        maxlr_wat = std::min(maxsr_wat, maxlr);
     }
     else
     {
@@ -363,9 +366,7 @@ static void reset_neighbor_lists(t_forcerec *fr, gmx_bool bResetSR, gmx_bool bRe
 
 static gmx_inline void new_i_nblist(t_nblist *nlist, atom_id i_atom, int shift, int gid)
 {
-    int    i, k, nri, nshift;
-
-    nri = nlist->nri;
+    int    nri = nlist->nri;
 
     /* Check whether we have to increase the i counter */
     if ((nri == -1) ||
@@ -594,14 +595,14 @@ put_in_list_at(gmx_bool              bHaveVdW[],
     t_nblist  *   vdwc_ww    = NULL;
     t_nblist  *   coul_ww    = NULL;
 
-    int           i, j, jcg, igid, gid, nbl_ind, ind_ij;
+    int           i, j, jcg, igid, gid, nbl_ind;
     atom_id       jj, jj0, jj1, i_atom;
-    int           i0, nicg, len;
+    int           i0, nicg;
 
     int          *cginfo;
     int          *type, *typeB;
     real         *charge, *chargeB;
-    real          qi, qiB, qq, rlj;
+    real          qi, qiB;
     gmx_bool      bFreeEnergy, bFree, bFreeJ, bNotEx, *bPert;
     gmx_bool      bDoVdW_i, bDoCoul_i, bDoCoul_i_sol;
     int           iwater, jwater;
@@ -1094,7 +1095,7 @@ put_in_list_adress(gmx_bool              bHaveVdW[],
                    gmx_bool              bLR,
                    gmx_bool              bDoVdW,
                    gmx_bool              bDoCoul,
-                   int                   solvent_opt)
+                   int       gmx_unused  solvent_opt)
 {
     /* The a[] index has been removed,
      * to put it back in i_atom should be a[i0] and jj should be a[jj].
@@ -1105,33 +1106,26 @@ put_in_list_adress(gmx_bool              bHaveVdW[],
     t_nblist  *   vdwc_adress  = NULL;
     t_nblist  *   vdw_adress   = NULL;
     t_nblist  *   coul_adress  = NULL;
-    t_nblist  *   vdwc_ww      = NULL;
-    t_nblist  *   coul_ww      = NULL;
 
     int           i, j, jcg, igid, gid, nbl_ind, nbl_ind_adress;
     atom_id       jj, jj0, jj1, i_atom;
-    int           i0, nicg, len;
+    int           i0, nicg;
 
     int          *cginfo;
-    int          *type, *typeB;
-    real         *charge, *chargeB;
+    int          *type;
+    real         *charge;
     real         *wf;
-    real          qi, qiB, qq, rlj;
-    gmx_bool      bFreeEnergy, bFree, bFreeJ, bNotEx, *bPert;
-    gmx_bool      bDoVdW_i, bDoCoul_i, bDoCoul_i_sol;
+    real          qi;
+    gmx_bool      bNotEx;
+    gmx_bool      bDoVdW_i, bDoCoul_i;
     gmx_bool      b_hybrid;
-    gmx_bool      j_all_atom;
-    int           iwater, jwater;
     t_nblist     *nlist, *nlist_adress;
     gmx_bool      bEnergyGroupCG;
 
     /* Copy some pointers */
     cginfo  = fr->cginfo;
     charge  = md->chargeA;
-    chargeB = md->chargeB;
     type    = md->typeA;
-    typeB   = md->typeB;
-    bPert   = md->bPerturbed;
     wf      = md->wf;
 
     /* Get atom range */
@@ -1141,8 +1135,6 @@ put_in_list_adress(gmx_bool              bHaveVdW[],
     /* Get the i charge group info */
     igid   = GET_CGINFO_GID(cginfo[icg]);
 
-    iwater = (solvent_opt != esolNO) ? GET_CGINFO_SOLOPT(cginfo[icg]) : esolNO;
-
     if (md->nPerturbed)
     {
         gmx_fatal(FARGS, "AdResS does not support free energy pertubation\n");
@@ -1585,18 +1577,18 @@ static real calc_image_tric(rvec xi, rvec xj, matrix box,
     /* Perform NINT operation, using trunc operation, therefore
      * we first add 2.5 then subtract 2 again
      */
-    tz  = dz*b_inv[ZZ] + h25;
+    tz  = static_cast<int>(dz*b_inv[ZZ] + h25);
     tz -= 2;
     dz -= tz*box[ZZ][ZZ];
     dy -= tz*box[ZZ][YY];
     dx -= tz*box[ZZ][XX];
 
-    ty  = dy*b_inv[YY] + h25;
+    ty  = static_cast<int>(dy*b_inv[YY] + h25);
     ty -= 2;
     dy -= ty*box[YY][YY];
     dx -= ty*box[YY][XX];
 
-    tx  = dx*b_inv[XX]+h25;
+    tx  = static_cast<int>(dx*b_inv[XX]+h25);
     tx -= 2;
     dx -= tx*box[XX][XX];
 
@@ -1625,9 +1617,9 @@ static real calc_image_rect(rvec xi, rvec xj, rvec box_size,
     /* Perform NINT operation, using trunc operation, therefore
      * we first add 1.5 then subtract 1 again
      */
-    tx = dx*b_inv[XX] + h15;
-    ty = dy*b_inv[YY] + h15;
-    tz = dz*b_inv[ZZ] + h15;
+    tx = static_cast<int>(dx*b_inv[XX] + h15);
+    ty = static_cast<int>(dy*b_inv[YY] + h15);
+    tz = static_cast<int>(dz*b_inv[ZZ] + h15);
     tx--;
     ty--;
     tz--;
@@ -1645,7 +1637,7 @@ static real calc_image_rect(rvec xi, rvec xj, rvec box_size,
     return r2;
 }
 
-static void add_simple(t_ns_buf *nsbuf, int nrj, atom_id cg_j,
+static void add_simple(t_ns_buf * nsbuf, int nrj, atom_id cg_j,
                        gmx_bool bHaveVdW[], int ngid, t_mdatoms *md,
                        int icg, int jgid, t_block *cgs, t_excl bexcl[],
                        int shift, t_forcerec *fr, put_in_list_t *put_in_list)
@@ -1673,7 +1665,6 @@ static void ns_inner_tric(rvec x[], int icg, int *i_egp_flags,
     int       j, nrj, jgid;
     int      *cginfo = fr->cginfo;
     atom_id   cg_j, *cgindex;
-    t_ns_buf *nsbuf;
 
     cgindex = cgs->index;
     shift   = CENTRAL;
@@ -1706,7 +1697,6 @@ static void ns_inner_rect(rvec x[], int icg, int *i_egp_flags,
     int       j, nrj, jgid;
     int      *cginfo = fr->cginfo;
     atom_id   cg_j, *cgindex;
-    t_ns_buf *nsbuf;
 
     cgindex = cgs->index;
     if (bBox)
@@ -1760,7 +1750,7 @@ static int ns_simple_core(t_forcerec *fr,
 {
     int          naaj, k;
     real         rlist2;
-    int          nsearch, icg, jcg, igid, i0, nri, nn;
+    int          nsearch, icg, igid, nn;
     int         *cginfo;
     t_ns_buf    *nsbuf;
     /* atom_id  *i_atoms; */
@@ -2044,8 +2034,8 @@ static void get_cutoff2(t_forcerec *fr, gmx_bool bDoLongRange,
         *rvdw2  = *rs2;
         *rcoul2 = *rs2;
     }
-    *rm2 = min(*rvdw2, *rcoul2);
-    *rl2 = max(*rvdw2, *rcoul2);
+    *rm2 = std::min(*rvdw2, *rcoul2);
+    *rl2 = std::max(*rvdw2, *rcoul2);
 }
 
 static void init_nsgrid_lists(t_forcerec *fr, int ngid, gmx_ns_t *ns)
@@ -2107,15 +2097,15 @@ static int nsgrid_core(t_commrec *cr, t_forcerec *fr,
 #endif
     int           dx0, dx1, dy0, dy1, dz0, dz1;
     int           Nx, Ny, Nz, shift = -1, j, nrj, nns, nn = -1;
-    real          gridx, gridy, gridz, grid_x, grid_y, grid_z;
+    real          gridx, gridy, gridz, grid_x, grid_y;
     real         *dcx2, *dcy2, *dcz2;
     int           zgi, ygi, xgi;
-    int           cg0, cg1, icg = -1, cgsnr, i0, igid, nri, naaj, max_jcg;
+    int           cg0, cg1, icg = -1, cgsnr, i0, igid, naaj, max_jcg;
     int           jcg0, jcg1, jjcg, cgj0, jgid;
     int          *grida, *gridnra, *gridind;
     gmx_bool      rvdw_lt_rcoul, rcoul_lt_rvdw;
-    rvec          xi, *cgcm, grid_offset;
-    real          r2, rs2, rvdw2, rcoul2, rm2, rl2, XI, YI, ZI, dcx, dcy, dcz, tmp1, tmp2;
+    rvec         *cgcm, grid_offset;
+    real          r2, rs2, rvdw2, rcoul2, rm2, rl2, XI, YI, ZI, tmp1, tmp2;
     int          *i_egp_flags;
     gmx_bool      bDomDec, bTriclinicX, bTriclinicY;
     ivec          ncpddc;
@@ -2167,7 +2157,6 @@ static int nsgrid_core(t_commrec *cr, t_forcerec *fr,
     gridz      = grid->cell_size[ZZ];
     grid_x     = 1/gridx;
     grid_y     = 1/gridy;
-    grid_z     = 1/gridz;
     copy_rvec(grid->cell_offset, grid_offset);
     copy_ivec(grid->ncpddc, ncpddc);
     dcx2       = grid->dcx2;
@@ -2213,7 +2202,7 @@ static int nsgrid_core(t_commrec *cr, t_forcerec *fr,
         else
         {
             if (d == XX &&
-                box[XX][XX] - fabs(box[YY][XX]) - fabs(box[ZZ][XX]) < sqrt(rl2))
+                box[XX][XX] - fabs(box[YY][XX]) - fabs(box[ZZ][XX]) < std::sqrt(rl2))
             {
                 shp[d] = 2;
             }
@@ -2512,7 +2501,6 @@ static int nsgrid_core(t_commrec *cr, t_forcerec *fr,
                 }
             }
         }
-        /* setexcl(nri,i_atoms,&top->atoms.excl,FALSE,bexcl); */
         setexcl(cgs->index[icg], cgs->index[icg+1], &top->excls, FALSE, bexcl);
     }
     /* No need to perform any left-over force calculations anymore (as we used to do here)
@@ -2548,7 +2536,6 @@ void init_ns(FILE *fplog, const t_commrec *cr,
 {
     int  mt, icg, nr_in_cg, maxcg, i, j, jcg, ngid, ncg;
     t_block *cgs;
-    char *ptr;
 
     /* Compute largest charge groups size (# atoms) */
     nr_in_cg = 1;
@@ -2557,7 +2544,7 @@ void init_ns(FILE *fplog, const t_commrec *cr,
         cgs = &mtop->moltype[mt].cgs;
         for (icg = 0; (icg < cgs->nr); icg++)
         {
-            nr_in_cg = max(nr_in_cg, (int)(cgs->index[icg+1]-cgs->index[icg]));
+            nr_in_cg = std::max(nr_in_cg, (int)(cgs->index[icg+1]-cgs->index[icg]));
         }
     }
 
@@ -2669,13 +2656,11 @@ int search_neighbours(FILE *log, t_forcerec *fr,
 {
     t_block  *cgs = &(top->cgs);
     rvec     box_size, grid_x0, grid_x1;
-    int      i, j, m, ngid;
+    int      m, ngid;
     real     min_size, grid_dens;
     int      nsearch;
     gmx_bool     bGrid;
-    char     *ptr;
-    gmx_bool     *i_egp_flags;
-    int      cg_start, cg_end, start, end;
+    int      start, end;
     gmx_ns_t *ns;
     t_grid   *grid;
     gmx_domdec_zones_t *dd_zones;
@@ -2700,7 +2685,7 @@ int search_neighbours(FILE *log, t_forcerec *fr,
         }
         if (!bGrid)
         {
-            min_size = min(box_size[XX], min(box_size[YY], box_size[ZZ]));
+            min_size = std::min(box_size[XX], std::min(box_size[YY], box_size[ZZ]));
             if (2*fr->rlistlong >= min_size)
             {
                 gmx_fatal(FARGS, "One of the box diagonal elements has become smaller than twice the cut-off length.");
similarity index 97%
rename from src/gromacs/mdlib/nsgrid.c
rename to src/gromacs/mdlib/nsgrid.cpp
index 75e689e8f4112e029502bd3310de05d301667817..b92d5a5c765bf73755c01e1355f10872341d303c 100644 (file)
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2013,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.
 #include <stdio.h>
 #include <stdlib.h>
 
+#include <cmath>
+
+#include <algorithm>
+
 #include "gromacs/domdec/domdec.h"
 #include "gromacs/fileio/pdbio.h"
 #include "gromacs/legacyheaders/macros.h"
@@ -90,7 +94,7 @@ static void calc_x_av_stddev(int n, rvec *x, rvec av, rvec stddev)
     for (d = 0; d < DIM; d++)
     {
         av[d]     = s1[d];
-        stddev[d] = sqrt(s2[d] - s1[d]*s1[d]);
+        stddev[d] = std::sqrt(s2[d] - s1[d]*s1[d]);
     }
 }
 
@@ -202,7 +206,6 @@ static void set_grid_sizes(matrix box, rvec izones_x0, rvec izones_x1, real rlis
 {
     int      i, j;
     gmx_bool bDD, bDDRect;
-    rvec     av, stddev;
     rvec     izones_size;
     real     inv_r_ideal, size, add_tric, radd;
 
@@ -218,7 +221,7 @@ static void set_grid_sizes(matrix box, rvec izones_x0, rvec izones_x1, real rlis
     }
 
     /* Use the ideal number of cg's per cell to set the ideal cell size */
-    inv_r_ideal = pow(grid_density/grid->ncg_ideal, 1.0/3.0);
+    inv_r_ideal = std::pow((real)(grid_density/grid->ncg_ideal), (real)(1.0/3.0));
     if (rlist > 0 && inv_r_ideal*rlist < 1)
     {
         inv_r_ideal = 1/rlist;
@@ -273,8 +276,11 @@ static void set_grid_sizes(matrix box, rvec izones_x0, rvec izones_x1, real rlis
 
             /* Check if the cell boundary in this direction is
              * perpendicular to the Cartesian axis.
+             * Since grid->npbcdim isan integer that in principle can take
+             * any value, we help the compiler avoid warnings and potentially
+             * optimize by ensuring that j < DIM here.
              */
-            for (j = i+1; j < grid->npbcdim; j++)
+            for (j = i+1; j < grid->npbcdim && j < DIM; j++)
             {
                 if (box[j][i] != 0)
                 {
@@ -357,7 +363,6 @@ static void set_grid_sizes(matrix box, rvec izones_x0, rvec izones_x1, real rlis
 
 t_grid *init_grid(FILE *fplog, t_forcerec *fr)
 {
-    int     d, m;
     char   *ptr;
     t_grid *grid;
 
@@ -475,7 +480,6 @@ void grid_first(FILE *fplog, t_grid *grid,
                 real rlistlong, real grid_density)
 {
     int    i, m;
-    ivec   cx;
 
     set_grid_sizes(box, izones_x0, izones_x1, rlistlong, dd, ddbox, grid, grid_density);
 
@@ -495,7 +499,7 @@ void grid_first(FILE *fplog, t_grid *grid,
         }
     }
 
-    m = max(grid->n[XX], max(grid->n[YY], grid->n[ZZ]));
+    m = std::max(grid->n[XX], std::max(grid->n[YY], grid->n[ZZ]));
     if (m > grid->dc_nalloc)
     {
         /* Allocate with double the initial size for box scaling */
@@ -642,7 +646,7 @@ void fill_grid(gmx_domdec_zones_t *dd_zones,
                int cg0, int cg1, rvec cg_cm[])
 {
     int       *cell_index;
-    int        nrx, nry, nrz;
+    int        nry, nrz;
     rvec       n_box, offset;
     int        zone, ccg0, ccg1, cg, d, not_used;
     ivec       shift0, useall, b0, b1, ind;
@@ -661,7 +665,6 @@ void fill_grid(gmx_domdec_zones_t *dd_zones,
     cell_index = grid->cell_index;
 
     /* Initiate cell borders */
-    nrx = grid->n[XX];
     nry = grid->n[YY];
     nrz = grid->n[ZZ];
     for (d = 0; d < DIM; d++)
@@ -689,7 +692,7 @@ void fill_grid(gmx_domdec_zones_t *dd_zones,
         {
             for (d = 0; d < DIM; d++)
             {
-                ind[d] = (cg_cm[cg][d] - offset[d])*n_box[d];
+                ind[d] = static_cast<int>((cg_cm[cg][d] - offset[d])*n_box[d]);
                 /* With pbc we should be done here.
                  * Without pbc cg's outside the grid
                  * should be assigned to the closest grid cell.
@@ -764,7 +767,7 @@ void fill_grid(gmx_domdec_zones_t *dd_zones,
                 bUse = TRUE;
                 for (d = 0; d < DIM; d++)
                 {
-                    ind[d] = (cg_cm[cg][d] - offset[d])*n_box[d];
+                    ind[d] = static_cast<int>((cg_cm[cg][d] - offset[d])*n_box[d]);
                     /* Here we have to correct for rounding problems,
                      * as this cg_cm to cell index operation is not necessarily
                      * binary identical to the operation for the DD zone assignment