From: Mark Abraham <mark.j.abraham@gmail.com>
Date: Thu, 4 Sep 2014 20:21:35 +0000 (+0200)
Subject: Convert domdec code to C++
X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?p=alexxy%2Fgromacs.git;a=commitdiff_plain;h=fdbd674832f07e7388761a004321edd1089e8d39

Convert domdec code to C++

Used std:: versions of at least some calls to min, max, sqrt and pow,
eliminated unused variables, used a named floating-point constant,
used a static_cast to make clear that we intend to truncate to integer
in one place, made some iteration variables local to their block,
protected some uses of fplog and fshift.

Eliminated some redundant #ifdef GMX_MPI, changed the range of others
to declare variables only where used.

Converted some neighbour-search constants to be static const real,
rather than macros, with the side effect of not having to deal
with overloaded std::sqrt for constants.

Introduced one work-around (duplicating the code of add_ifunc) for
dealing with the way that nral is a complicated function of
ftype. Static analysis was not being consistent in the way it was
choosing code paths with possible values of nral in different
functions. Co-locating two loops over nral resolves the issue.

Added assertions that help static analysis understand that settles are
kinds of constraints and that sane things happen elsewhere, so that
code paths that look like they might be called when only
inter-charge-group settles exist are actually well formed.

Added yet more assertions to reassure the analyzer that we handle
memory allocation with charge-group linking correctly.

Clarified a PBC-related variable name.

Change-Id: I5d1cf06c866098d183d7fe5aa023a292ca3d109b
---

diff --git a/src/gromacs/legacyheaders/nsgrid.h b/src/gromacs/legacyheaders/nsgrid.h
index 23c2fca243..c00df0078f 100644
--- a/src/gromacs/legacyheaders/nsgrid.h
+++ b/src/gromacs/legacyheaders/nsgrid.h
@@ -41,16 +41,22 @@
 extern "C" {
 #endif
 
-#define   GRID_STDDEV_FAC  sqrt(3)
-#define NSGRID_STDDEV_FAC  2.0
-/*
- * GRID_STDDEV_FAC * stddev is used to estimate the interaction density.
- * sqrt(3) gives a uniform load for a rectangular block of cg's.
- * For a sphere it is not a bad approximation for 4x1x1 up to 4x2x2.
+/*! \brief Used when estimating the interaction density.
+ *
+ * GRID_STDDEV_FAC * stddev estimates the interaction density. The
+ * value sqrt(3) == 1.73205080757 gives a uniform load for a
+ * rectangular 3D block of charge groups. For a sphere, it is not a
+ * bad approximation for 4x1x1 up to 4x2x2.
  *
- * The extent of the neighborsearch grid is a bit larger than sqrt(3)
+ * \todo It would be nicer to use sqrt(3) here, when all code that
+ * #includes this file is in C++, which will let us cope with the
+ * std::sqrt<T> on Windows. */
+static const real GRID_STDDEV_FAC = 1.73205080757;
+
+/*! \brief The extent of the neighborsearch grid is a bit larger than sqrt(3)
  * to account for less dense regions at the edges of the system.
  */
+static const real NSGRID_STDDEV_FAC = 2.0;
 
 #define NSGRID_SIGNAL_MOVED_FAC  4
 /* A cell index of NSGRID_SIGNAL_MOVED_FAC*ncells signals
diff --git a/src/gromacs/mdlib/domdec.c b/src/gromacs/mdlib/domdec.cpp
similarity index 97%
rename from src/gromacs/mdlib/domdec.c
rename to src/gromacs/mdlib/domdec.cpp
index bad4092e3e..6c182cc12d 100644
--- a/src/gromacs/mdlib/domdec.c
+++ b/src/gromacs/mdlib/domdec.cpp
@@ -39,6 +39,8 @@
 
 #include "config.h"
 
+#include <algorithm>
+
 #include <assert.h>
 #include <math.h>
 #include <stdio.h>
@@ -770,7 +772,7 @@ void dd_move_f(gmx_domdec_t *dd, rvec f[], rvec *fshift)
     rvec                  *buf, *sbuf;
     ivec                   vis;
     int                    is;
-    gmx_bool               bPBC, bScrew;
+    gmx_bool               bShiftForcesNeedPbc, bScrew;
 
     comm = dd->comm;
 
@@ -778,16 +780,17 @@ void dd_move_f(gmx_domdec_t *dd, rvec f[], rvec *fshift)
 
     buf = comm->vbuf.v;
 
-    n       = 0;
     nzone   = comm->zones.n/2;
     nat_tot = dd->nat_tot;
     for (d = dd->ndim-1; d >= 0; d--)
     {
-        bPBC   = (dd->ci[dd->dim[d]] == 0);
-        bScrew = (bPBC && dd->bScrewPBC && dd->dim[d] == XX);
+        /* Only forces in domains near the PBC boundaries need to
+           consider PBC in the treatment of fshift */
+        bShiftForcesNeedPbc   = (dd->ci[dd->dim[d]] == 0);
+        bScrew = (bShiftForcesNeedPbc && dd->bScrewPBC && dd->dim[d] == XX);
         if (fshift == NULL && !bScrew)
         {
-            bPBC = FALSE;
+            bShiftForcesNeedPbc = FALSE;
         }
         /* Determine which shift vector we need */
         clear_ivec(vis);
@@ -823,7 +826,7 @@ void dd_move_f(gmx_domdec_t *dd, rvec f[], rvec *fshift)
             index = ind->index;
             /* Add the received forces */
             n = 0;
-            if (!bPBC)
+            if (!bShiftForcesNeedPbc)
             {
                 for (i = 0; i < ind->nsend[nzone]; i++)
                 {
@@ -838,6 +841,9 @@ void dd_move_f(gmx_domdec_t *dd, rvec f[], rvec *fshift)
             }
             else if (!bScrew)
             {
+                /* fshift should always be defined if this function is
+                 * called when bShiftForcesNeedPbc is true */
+                assert(NULL != fshift);
                 for (i = 0; i < ind->nsend[nzone]; i++)
                 {
                     at0 = cgindex[index[i]];
@@ -958,7 +964,6 @@ void dd_atom_sum_real(gmx_domdec_t *dd, real v[])
 
     buf = &comm->vbuf.v[0][0];
 
-    n       = 0;
     nzone   = comm->zones.n/2;
     nat_tot = dd->nat_tot;
     for (d = dd->ndim-1; d >= 0; d--)
@@ -1064,7 +1069,7 @@ static void dd_sendrecv_ddzone(const gmx_domdec_t *dd,
 static void dd_move_cellx(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
                           rvec cell_ns_x0, rvec cell_ns_x1)
 {
-    int                d, d1, dim, dim1, pos, buf_size, i, j, k, p, npulse, npulse_min;
+    int                d, d1, dim, pos, buf_size, i, j, p, npulse, npulse_min;
     gmx_ddzone_t      *zp;
     gmx_ddzone_t       buf_s[DDZONECOMM_MAXZONE];
     gmx_ddzone_t       buf_r[DDZONECOMM_MAXZONE];
@@ -1136,7 +1141,7 @@ static void dd_move_cellx(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
         if (bPBC)
         {
             /* Take the minimum to avoid double communication */
-            npulse_min = min(npulse, dd->nc[dim]-1-npulse);
+            npulse_min = std::min(npulse, dd->nc[dim]-1-npulse);
         }
         else
         {
@@ -1160,9 +1165,9 @@ static void dd_move_cellx(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
             {
                 for (d1 = d; d1 < dd->ndim-1; d1++)
                 {
-                    extr_s[d1][0] = max(extr_s[d1][0], extr_r[d1][0]);
-                    extr_s[d1][1] = min(extr_s[d1][1], extr_r[d1][1]);
-                    extr_s[d1][2] = min(extr_s[d1][2], extr_r[d1][2]);
+                    extr_s[d1][0] = std::max(extr_s[d1][0], extr_r[d1][0]);
+                    extr_s[d1][1] = std::min(extr_s[d1][1], extr_r[d1][1]);
+                    extr_s[d1][2] = std::min(extr_s[d1][2], extr_r[d1][2]);
                 }
             }
         }
@@ -1223,9 +1228,9 @@ static void dd_move_cellx(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
                 {
                     if (bUse)
                     {
-                        buf_e[i].min0 = min(buf_e[i].min0, buf_r[i].min0);
-                        buf_e[i].max1 = max(buf_e[i].max1, buf_r[i].max1);
-                        buf_e[i].min1 = min(buf_e[i].min1, buf_r[i].min1);
+                        buf_e[i].min0 = std::min(buf_e[i].min0, buf_r[i].min0);
+                        buf_e[i].max1 = std::max(buf_e[i].max1, buf_r[i].max1);
+                        buf_e[i].min1 = std::min(buf_e[i].min1, buf_r[i].min1);
                     }
 
                     if (dd->ndim == 3 && d == 0 && i == buf_size - 1)
@@ -1238,8 +1243,8 @@ static void dd_move_cellx(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
                     }
                     if (bUse && dh[d1] >= 0)
                     {
-                        buf_e[i].mch0 = max(buf_e[i].mch0, buf_r[i].mch0-dh[d1]);
-                        buf_e[i].mch1 = max(buf_e[i].mch1, buf_r[i].mch1-dh[d1]);
+                        buf_e[i].mch0 = std::max(buf_e[i].mch0, buf_r[i].mch0-dh[d1]);
+                        buf_e[i].mch1 = std::max(buf_e[i].mch1, buf_r[i].mch1-dh[d1]);
                     }
                 }
                 /* Copy the received buffer to the send buffer,
@@ -1255,9 +1260,9 @@ static void dd_move_cellx(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
 
                 for (d1 = d; d1 < dd->ndim-1; d1++)
                 {
-                    extr_s[d1][1] = min(extr_s[d1][1], buf_e[pos].min0);
-                    extr_s[d1][0] = max(extr_s[d1][0], buf_e[pos].max1);
-                    extr_s[d1][2] = min(extr_s[d1][2], buf_e[pos].min1);
+                    extr_s[d1][1] = std::min(extr_s[d1][1], buf_e[pos].min0);
+                    extr_s[d1][0] = std::max(extr_s[d1][0], buf_e[pos].max1);
+                    extr_s[d1][2] = std::min(extr_s[d1][2], buf_e[pos].min1);
                     pos++;
                 }
 
@@ -1287,8 +1292,8 @@ static void dd_move_cellx(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
             {
                 print_ddzone(debug, 1, i, 0, &comm->zone_d1[i]);
             }
-            cell_ns_x0[dim] = min(cell_ns_x0[dim], comm->zone_d1[i].min0);
-            cell_ns_x1[dim] = max(cell_ns_x1[dim], comm->zone_d1[i].max1);
+            cell_ns_x0[dim] = std::min(cell_ns_x0[dim], comm->zone_d1[i].min0);
+            cell_ns_x1[dim] = std::max(cell_ns_x1[dim], comm->zone_d1[i].max1);
         }
     }
     if (dd->ndim >= 3)
@@ -1302,8 +1307,8 @@ static void dd_move_cellx(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
                 {
                     print_ddzone(debug, 2, i, j, &comm->zone_d2[i][j]);
                 }
-                cell_ns_x0[dim] = min(cell_ns_x0[dim], comm->zone_d2[i][j].min0);
-                cell_ns_x1[dim] = max(cell_ns_x1[dim], comm->zone_d2[i][j].max1);
+                cell_ns_x0[dim] = std::min(cell_ns_x0[dim], comm->zone_d2[i][j].min0);
+                cell_ns_x1[dim] = std::max(cell_ns_x1[dim], comm->zone_d2[i][j].max1);
             }
         }
     }
@@ -1529,10 +1534,6 @@ static void dd_collect_vec_gatherv(gmx_domdec_t *dd,
 void dd_collect_vec(gmx_domdec_t *dd,
                     t_state *state_local, rvec *lv, rvec *v)
 {
-    gmx_domdec_master_t *ma;
-    int                  n, i, c, a, nalloc = 0;
-    rvec                *buf = NULL;
-
     dd_collect_cg(dd, state_local);
 
     if (dd->nnodes <= GMX_DD_NNODES_SENDRECV)
@@ -1754,7 +1755,7 @@ static void dd_distribute_vec_scatterv(gmx_domdec_t *dd, t_block *cgs,
 {
     gmx_domdec_master_t *ma;
     int                 *scounts = NULL, *disps = NULL;
-    int                  n, i, c, a, nalloc = 0;
+    int                  n, i, c, a;
     rvec                *buf = NULL;
 
     if (DDMASTER(dd))
@@ -2091,15 +2092,15 @@ real dd_cutoff_mbody(gmx_domdec_t *dd)
             r = comm->cellsize_min[dd->dim[0]];
             for (di = 1; di < dd->ndim; di++)
             {
-                r = min(r, comm->cellsize_min[dd->dim[di]]);
+                r = std::min(r, comm->cellsize_min[dd->dim[di]]);
             }
             if (comm->bBondComm)
             {
-                r = max(r, comm->cutoff_mbody);
+                r = std::max(r, comm->cutoff_mbody);
             }
             else
             {
-                r = min(r, comm->cutoff);
+                r = std::min(r, comm->cutoff);
             }
         }
     }
@@ -2113,7 +2114,7 @@ real dd_cutoff_twobody(gmx_domdec_t *dd)
 
     r_mb = dd_cutoff_mbody(dd);
 
-    return max(dd->comm->cutoff, r_mb);
+    return std::max(dd->comm->cutoff, r_mb);
 }
 
 
@@ -2175,7 +2176,7 @@ static int *dd_pmenodes(t_commrec *cr)
 static int gmx_ddcoord2pmeindex(t_commrec *cr, int x, int y, int z)
 {
     gmx_domdec_t *dd;
-    ivec          coords, coords_pme, nc;
+    ivec          coords;
     int           slab;
 
     dd = cr->dd;
@@ -2244,7 +2245,6 @@ static int dd_simnode2pmenode(t_commrec *cr, int sim_nodeid)
 {
     gmx_domdec_t      *dd;
     gmx_domdec_comm_t *comm;
-    ivec               coord, coord_pme;
     int                i;
     int                pmenode = -1;
 
@@ -2255,6 +2255,7 @@ static int dd_simnode2pmenode(t_commrec *cr, int sim_nodeid)
     if (comm->bCartesianPP_PME)
     {
 #ifdef GMX_MPI
+        ivec coord, coord_pme;
         MPI_Cart_coords(cr->mpi_comm_mysim, sim_nodeid, DIM, coord);
         if (coord[comm->cartpmedim] < dd->nc[comm->cartpmedim])
         {
@@ -2392,7 +2393,7 @@ void get_pme_ddnodes(t_commrec *cr, int pmenodeid,
 static gmx_bool receive_vir_ener(t_commrec *cr)
 {
     gmx_domdec_comm_t *comm;
-    int                pmenode, coords[DIM], rank;
+    int                pmenode;
     gmx_bool           bReceive;
 
     bReceive = TRUE;
@@ -2403,10 +2404,12 @@ static gmx_bool receive_vir_ener(t_commrec *cr)
         {
             pmenode = dd_simnode2pmenode(cr, cr->sim_nodeid);
 #ifdef GMX_MPI
+            ivec coords;
             MPI_Cart_coords(cr->mpi_comm_mysim, cr->sim_nodeid, DIM, coords);
             coords[comm->cartpmedim]++;
             if (coords[comm->cartpmedim] < cr->dd->nc[comm->cartpmedim])
             {
+                int rank;
                 MPI_Cart_rank(cr->mpi_comm_mysim, coords, &rank);
                 if (dd_simnode2pmenode(cr, rank) == pmenode)
                 {
@@ -2514,12 +2517,8 @@ static void make_dd_indices(gmx_domdec_t *dd,
 {
     int          nzone, zone, zone1, cg0, cg1, cg1_p1, cg, cg_gl, a, a_gl;
     int         *zone2cg, *zone_ncg1, *index_gl, *gatindex;
-    gmx_ga2la_t *ga2la;
-    char        *bLocalCG;
     gmx_bool     bCGs;
 
-    bLocalCG = dd->comm->bLocalCG;
-
     if (dd->nat_tot > dd->gatindex_nalloc)
     {
         dd->gatindex_nalloc = over_alloc_dd(dd->nat_tot);
@@ -2584,7 +2583,7 @@ static void make_dd_indices(gmx_domdec_t *dd,
 static int check_bLocalCG(gmx_domdec_t *dd, int ncg_sys, const char *bLocalCG,
                           const char *where)
 {
-    int ncg, i, ngl, nerr;
+    int i, ngl, nerr;
 
     nerr = 0;
     if (bLocalCG == NULL)
@@ -2744,12 +2743,12 @@ static real cellsize_min_dlb(gmx_domdec_comm_t *comm, int dim_ind, int dim)
         /* The cut-off might have changed, e.g. by PME load balacning,
          * from the value used to set comm->cellsize_min, so check it.
          */
-        cellsize_min = max(cellsize_min, comm->cutoff/comm->cd[dim_ind].np_dlb);
+        cellsize_min = std::max(cellsize_min, comm->cutoff/comm->cd[dim_ind].np_dlb);
 
         if (comm->bPMELoadBalDLBLimits)
         {
             /* Check for the cut-off limit set by the PME load balancing */
-            cellsize_min = max(cellsize_min, comm->PMELoadBal_max_cutoff/comm->cd[dim_ind].np_dlb);
+            cellsize_min = std::max(cellsize_min, comm->PMELoadBal_max_cutoff/comm->cd[dim_ind].np_dlb);
         }
     }
 
@@ -2772,10 +2771,10 @@ static real grid_jump_limit(gmx_domdec_comm_t *comm, real cutoff,
     {
         if (comm->bPMELoadBalDLBLimits)
         {
-            cutoff = max(cutoff, comm->PMELoadBal_max_cutoff);
+            cutoff = std::max(cutoff, comm->PMELoadBal_max_cutoff);
         }
-        grid_jump_limit = max(grid_jump_limit,
-                              cutoff/comm->cd[dim_ind].np);
+        grid_jump_limit = std::max(grid_jump_limit,
+                                   cutoff/comm->cd[dim_ind].np);
     }
 
     return grid_jump_limit;
@@ -2962,8 +2961,8 @@ static void init_ddpme(gmx_domdec_t *dd, gmx_ddpme_t *ddpme, int dimind)
             {
                 slab = pmeindex % ddpme->nslab;
             }
-            ddpme->pp_min[slab] = min(ddpme->pp_min[slab], xyz[dimind]);
-            ddpme->pp_max[slab] = max(ddpme->pp_max[slab], xyz[dimind]);
+            ddpme->pp_min[slab] = std::min(ddpme->pp_min[slab], xyz[dimind]);
+            ddpme->pp_max[slab] = std::max(ddpme->pp_max[slab], xyz[dimind]);
         }
     }
 
@@ -3164,7 +3163,7 @@ static void set_dd_cell_sizes_slb(gmx_domdec_t *dd, gmx_ddbox_t *ddbox,
                 {
                     npulse[d]++;
                 }
-                cellsize_min[d] = min(cellsize_min[d], cellsize);
+                cellsize_min[d] = std::min(cellsize_min[d], cellsize);
             }
             if (setmode == setcellsizeslbLOCAL)
             {
@@ -3434,7 +3433,7 @@ static void set_dd_cell_sizes_dlb_root(gmx_domdec_t *dd,
                                        gmx_bool bUniform, gmx_int64_t step)
 {
     gmx_domdec_comm_t *comm;
-    int                ncd, d1, i, j, pos;
+    int                ncd, d1, i, pos;
     real              *cell_size;
     real               load_aver, load_i, imbalance, change, change_max, sc;
     real               cellsize_limit_f, dist_min_f, dist_min_f_hard, space;
@@ -3477,7 +3476,7 @@ static void set_dd_cell_sizes_dlb_root(gmx_domdec_t *dd,
             imbalance = (load_i - load_aver)/(load_aver > 0 ? load_aver : 1);
             /* Determine the change of the cell size using underrelaxation */
             change     = -relax*imbalance;
-            change_max = max(change_max, max(change, -change));
+            change_max = std::max(change_max, std::max(change, -change));
         }
         /* Limit the amount of scaling.
          * We need to use the same rescaling for all cells in one row,
@@ -3620,7 +3619,7 @@ static void distribute_dd_cell_sizes_dlb(gmx_domdec_t *dd,
                                          gmx_ddbox_t *ddbox)
 {
     gmx_domdec_comm_t *comm;
-    int                d1, dim1, pos;
+    int                d1, pos;
 
     comm = dd->comm;
 
@@ -3908,13 +3907,13 @@ static void distribute_cg(FILE *fplog, gmx_int64_t step,
 {
     gmx_domdec_master_t *ma;
     int                **tmp_ind = NULL, *tmp_nalloc = NULL;
-    int                  i, icg, j, k, k0, k1, d, npbcdim;
+    int                  i, icg, j, k, k0, k1, d;
     matrix               tcm;
-    rvec                 box_size, cg_cm;
+    rvec                 cg_cm;
     ivec                 ind;
     real                 nrcg, inv_ncg, pos_d;
     atom_id             *cgindex;
-    gmx_bool             bUnbounded, bScrew;
+    gmx_bool             bScrew;
 
     ma = dd->ma;
 
@@ -4122,9 +4121,9 @@ static void get_cg_distribution(FILE *fplog, gmx_int64_t step, gmx_domdec_t *dd,
     }
 
     dd_scatterv(dd,
-                DDMASTER(dd) ? ma->ibuf : NULL,
-                DDMASTER(dd) ? ma->ibuf+dd->nnodes : NULL,
-                DDMASTER(dd) ? ma->cg : NULL,
+                bMaster ? ma->ibuf : NULL,
+                bMaster ? ma->ibuf+dd->nnodes : NULL,
+                bMaster ? ma->cg : NULL,
                 dd->ncg_home*sizeof(int), dd->index_gl);
 
     /* Determine the home charge group sizes */
@@ -4455,8 +4454,8 @@ static void calc_cg_move(FILE *fplog, gmx_int64_t step,
                          int *move)
 {
     int      npbcdim;
-    int      c, i, cg, k, k0, k1, d, dim, dim2, dir, d2, d3, d4, cell_d;
-    int      mc, cdd, nrcg, ncg_recv, nat_recv, nvs, nvr, nvec, vec;
+    int      cg, k, k0, k1, d, dim, d2;
+    int      mc, nrcg;
     int      flag;
     gmx_bool bScrew;
     ivec     dev;
@@ -4630,17 +4629,15 @@ static void dd_redistribute_cg(FILE *fplog, gmx_int64_t step,
     int               *move;
     int                npbcdim;
     int                ncg[DIM*2], nat[DIM*2];
-    int                c, i, cg, k, k0, k1, d, dim, dim2, dir, d2, d3, d4, cell_d;
-    int                mc, cdd, nrcg, ncg_recv, nat_recv, nvs, nvr, nvec, vec;
+    int                c, i, cg, k, d, dim, dim2, dir, d2, d3;
+    int                mc, cdd, nrcg, ncg_recv, nvs, nvr, nvec, vec;
     int                sbuf[2], rbuf[2];
     int                home_pos_cg, home_pos_at, buf_pos;
     int                flag;
     gmx_bool           bV = FALSE, bSDX = FALSE, bCGP = FALSE;
-    gmx_bool           bScrew;
-    ivec               dev;
-    real               inv_ncg, pos_d;
+    real               pos_d;
     matrix             tcm;
-    rvec              *cg_cm = NULL, cell_x0, cell_x1, limitd, limit0, limit1, cm_new;
+    rvec              *cg_cm = NULL, cell_x0, cell_x1, limitd, limit0, limit1;
     atom_id           *cgindex;
     cginfo_mb_t       *cginfo_mb;
     gmx_domdec_comm_t *comm;
@@ -4897,7 +4894,6 @@ static void dd_redistribute_cg(FILE *fplog, gmx_int64_t step,
     {
         dim      = dd->dim[d];
         ncg_recv = 0;
-        nat_recv = 0;
         nvr      = 0;
         for (dir = 0; dir < (dd->nc[dim] == 2 ? 1 : 2); dir++)
         {
@@ -4932,7 +4928,6 @@ static void dd_redistribute_cg(FILE *fplog, gmx_int64_t step,
                              comm->cgcm_state[cdd], nvs,
                              comm->vbuf.v+nvr, i);
             ncg_recv += rbuf[0];
-            nat_recv += rbuf[1];
             nvr      += i;
         }
 
@@ -5236,7 +5231,7 @@ static void get_load_distribution(gmx_domdec_t *dd, gmx_wallcycle_t wcycle)
     gmx_domdec_comm_t *comm;
     gmx_domdec_load_t *load;
     gmx_domdec_root_t *root = NULL;
-    int                d, dim, cid, i, pos;
+    int                d, dim, i, pos;
     float              cell_frac = 0, sbuf[DD_NLOAD_MAX];
     gmx_bool           bSepPME;
 
@@ -5332,7 +5327,7 @@ static void get_load_distribution(gmx_domdec_t *dd, gmx_wallcycle_t wcycle)
                 for (i = 0; i < dd->nc[dim]; i++)
                 {
                     load->sum += load->load[pos++];
-                    load->max  = max(load->max, load->load[pos]);
+                    load->max  = std::max(load->max, load->load[pos]);
                     pos++;
                     if (dd->bGridJump)
                     {
@@ -5341,14 +5336,14 @@ static void get_load_distribution(gmx_domdec_t *dd, gmx_wallcycle_t wcycle)
                             /* This direction could not be load balanced properly,
                              * therefore we need to use the maximum iso the average load.
                              */
-                            load->sum_m = max(load->sum_m, load->load[pos]);
+                            load->sum_m = std::max(load->sum_m, load->load[pos]);
                         }
                         else
                         {
                             load->sum_m += load->load[pos];
                         }
                         pos++;
-                        load->cvol_min = min(load->cvol_min, load->load[pos]);
+                        load->cvol_min = std::min(load->cvol_min, load->load[pos]);
                         pos++;
                         if (d < dd->ndim-1)
                         {
@@ -5362,9 +5357,9 @@ static void get_load_distribution(gmx_domdec_t *dd, gmx_wallcycle_t wcycle)
                     }
                     if (bSepPME)
                     {
-                        load->mdf = max(load->mdf, load->load[pos]);
+                        load->mdf = std::max(load->mdf, load->load[pos]);
                         pos++;
-                        load->pme = max(load->pme, load->load[pos]);
+                        load->pme = std::max(load->pme, load->load[pos]);
                         pos++;
                     }
                 }
@@ -5758,7 +5753,6 @@ static void make_load_communicators(gmx_domdec_t gmx_unused *dd)
 
 void setup_dd_grid(FILE *fplog, gmx_domdec_t *dd)
 {
-    gmx_bool                bZYX;
     int                     d, dim, i, j, m;
     ivec                    tmp, s;
     int                     nzone, nzonep;
@@ -5917,17 +5911,16 @@ void setup_dd_grid(FILE *fplog, gmx_domdec_t *dd)
 static void make_pp_communicator(FILE *fplog, t_commrec *cr, int gmx_unused reorder)
 {
     gmx_domdec_t      *dd;
+    dd   = cr->dd;
+
+#ifdef GMX_MPI
     gmx_domdec_comm_t *comm;
-    int                i, rank, *buf;
+    int                rank, *buf;
     ivec               periods;
-#ifdef GMX_MPI
     MPI_Comm           comm_cart;
-#endif
 
-    dd   = cr->dd;
     comm = dd->comm;
 
-#ifdef GMX_MPI
     if (comm->bCartesianPP)
     {
         /* Set up cartesian communication for the particle-particle part */
@@ -5937,7 +5930,7 @@ static void make_pp_communicator(FILE *fplog, t_commrec *cr, int gmx_unused reor
                     dd->nc[XX], dd->nc[YY], dd->nc[ZZ]);
         }
 
-        for (i = 0; i < DIM; i++)
+        for (int i = 0; i < DIM; i++)
         {
             periods[i] = TRUE;
         }
@@ -6001,7 +5994,7 @@ static void make_pp_communicator(FILE *fplog, t_commrec *cr, int gmx_unused reor
         /* Determine the master coordinates and rank.
          * The DD master should be the same node as the master of this sim.
          */
-        for (i = 0; i < dd->nnodes; i++)
+        for (int i = 0; i < dd->nnodes; i++)
         {
             if (comm->ddindex2simnodeid[i] == 0)
             {
@@ -6039,30 +6032,27 @@ static void make_pp_communicator(FILE *fplog, t_commrec *cr, int gmx_unused reor
     }
 }
 
-static void receive_ddindex2simnodeid(t_commrec *cr)
+static void receive_ddindex2simnodeid(t_commrec gmx_unused *cr)
 {
+#ifdef GMX_MPI
     gmx_domdec_t      *dd;
-
     gmx_domdec_comm_t *comm;
-    int               *buf;
 
     dd   = cr->dd;
     comm = dd->comm;
 
-#ifdef GMX_MPI
     if (!comm->bCartesianPP_PME && comm->bCartesianPP)
     {
+        int *buf;
         snew(comm->ddindex2simnodeid, dd->nnodes);
         snew(buf, dd->nnodes);
         if (cr->duty & DUTY_PP)
         {
             buf[dd_index(dd->nc, dd->ci)] = cr->sim_nodeid;
         }
-#ifdef GMX_MPI
         /* Communicate the ddindex to simulation nodeid index */
         MPI_Allreduce(buf, comm->ddindex2simnodeid, dd->nnodes, MPI_INT, MPI_SUM,
                       cr->mpi_comm_mysim);
-#endif
         sfree(buf);
     }
 #endif
@@ -6104,9 +6094,8 @@ static void split_communicator(FILE *fplog, t_commrec *cr, int gmx_unused dd_nod
 {
     gmx_domdec_t      *dd;
     gmx_domdec_comm_t *comm;
-    int                i, rank;
+    int                i;
     gmx_bool           bDiv[DIM];
-    ivec               periods;
 #ifdef GMX_MPI
     MPI_Comm           comm_cart;
 #endif
@@ -6154,6 +6143,9 @@ static void split_communicator(FILE *fplog, t_commrec *cr, int gmx_unused dd_nod
 #ifdef GMX_MPI
     if (comm->bCartesianPP_PME)
     {
+        int  rank;
+        ivec periods;
+
         if (fplog)
         {
             fprintf(fplog, "Will use a Cartesian communicator for PP <-> PME: %d x %d x %d\n", comm->ntot[XX], comm->ntot[YY], comm->ntot[ZZ]);
@@ -6165,7 +6157,6 @@ static void split_communicator(FILE *fplog, t_commrec *cr, int gmx_unused dd_nod
         }
         MPI_Cart_create(cr->mpi_comm_mysim, DIM, comm->ntot, periods, reorder,
                         &comm_cart);
-
         MPI_Comm_rank(comm_cart, &rank);
         if (MASTERNODE(cr) && rank != 0)
         {
@@ -6349,7 +6340,7 @@ static real *get_slb_frac(FILE *fplog, const char *dir, int nc, const char *size
         for (i = 0; i < nc; i++)
         {
             dbl = 0;
-            sscanf(size_string, "%lf%n", &dbl, &n);
+            sscanf(size_string, "%20lf%n", &dbl, &n);
             if (dbl == 0)
             {
                 gmx_fatal(FARGS, "Incorrect or not enough DD cell size entries for direction %s: '%s'", dir, size_string);
@@ -6412,7 +6403,7 @@ static int dd_getenv(FILE *fplog, const char *env_var, int def)
     val = getenv(env_var);
     if (val)
     {
-        if (sscanf(val, "%d", &nst) <= 0)
+        if (sscanf(val, "%20d", &nst) <= 0)
         {
             nst = 1;
         }
@@ -6473,7 +6464,7 @@ static real average_cellsize_min(gmx_domdec_t *dd, gmx_ddbox_t *ddbox)
     {
         d = dd->dim[di];
         /* Check using the initial average cell size */
-        r = min(r, ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
+        r = std::min(r, ddbox->box_size[d]*ddbox->skew_fac[d]/dd->nc[d]);
     }
 
     return r;
@@ -6483,7 +6474,6 @@ static int check_dlb_support(FILE *fplog, t_commrec *cr,
                              const char *dlb_opt, gmx_bool bRecordLoad,
                              unsigned long Flags, t_inputrec *ir)
 {
-    gmx_domdec_t *dd;
     int           eDLB = -1;
     char          buf[STRLEN];
 
@@ -6625,10 +6615,10 @@ gmx_domdec_t *init_domain_decomposition(FILE *fplog, t_commrec *cr,
     gmx_domdec_t      *dd;
     gmx_domdec_comm_t *comm;
     int                recload;
-    int                d, i, j;
     real               r_2b, r_mb, r_bonded = -1, r_bonded_limit = -1, limit, acs;
     gmx_bool           bC;
     char               buf[STRLEN];
+    const real         tenPercentMargin = 1.1;
 
     if (fplog)
     {
@@ -6760,7 +6750,7 @@ gmx_domdec_t *init_domain_decomposition(FILE *fplog, t_commrec *cr,
             }
             else
             {
-                comm->cutoff = max(comm->cutoff, comm->cutoff_mbody);
+                comm->cutoff = std::max(comm->cutoff, comm->cutoff_mbody);
             }
             r_bonded_limit = comm->cutoff_mbody;
         }
@@ -6786,16 +6776,16 @@ gmx_domdec_t *init_domain_decomposition(FILE *fplog, t_commrec *cr,
              */
             if (Flags & MD_DDBONDCOMM)
             {
-                if (max(r_2b, r_mb) > comm->cutoff)
+                if (std::max(r_2b, r_mb) > comm->cutoff)
                 {
-                    r_bonded        = max(r_2b, r_mb);
-                    r_bonded_limit  = 1.1*r_bonded;
+                    r_bonded        = std::max(r_2b, r_mb);
+                    r_bonded_limit  = tenPercentMargin*r_bonded;
                     comm->bBondComm = TRUE;
                 }
                 else
                 {
                     r_bonded       = r_mb;
-                    r_bonded_limit = min(1.1*r_bonded, comm->cutoff);
+                    r_bonded_limit = std::min(tenPercentMargin*r_bonded, comm->cutoff);
                 }
                 /* We determine cutoff_mbody later */
             }
@@ -6804,12 +6794,12 @@ gmx_domdec_t *init_domain_decomposition(FILE *fplog, t_commrec *cr,
                 /* No special bonded communication,
                  * simply increase the DD cut-off.
                  */
-                r_bonded_limit     = 1.1*max(r_2b, r_mb);
+                r_bonded_limit     = tenPercentMargin*std::max(r_2b, r_mb);
                 comm->cutoff_mbody = r_bonded_limit;
-                comm->cutoff       = max(comm->cutoff, comm->cutoff_mbody);
+                comm->cutoff       = std::max(comm->cutoff, comm->cutoff_mbody);
             }
         }
-        comm->cellsize_limit = max(comm->cellsize_limit, r_bonded_limit);
+        comm->cellsize_limit = std::max(comm->cellsize_limit, r_bonded_limit);
         if (fplog)
         {
             fprintf(fplog,
@@ -6843,7 +6833,7 @@ gmx_domdec_t *init_domain_decomposition(FILE *fplog, t_commrec *cr,
                 "User supplied maximum distance required for P-LINCS: %.3f nm\n",
                 rconstr);
     }
-    comm->cellsize_limit = max(comm->cellsize_limit, rconstr);
+    comm->cellsize_limit = std::max(comm->cellsize_limit, rconstr);
 
     comm->cgs_gl = gmx_mtop_global_cgs(mtop);
 
@@ -6999,12 +6989,12 @@ gmx_domdec_t *init_domain_decomposition(FILE *fplog, t_commrec *cr,
             if (comm->eDLB != edlbNO)
             {
                 /* Check if this does not limit the scaling */
-                comm->cutoff_mbody = min(comm->cutoff_mbody, dlb_scale*acs);
+                comm->cutoff_mbody = std::min(comm->cutoff_mbody, dlb_scale*acs);
             }
             if (!comm->bBondComm)
             {
                 /* Without bBondComm do not go beyond the n.b. cut-off */
-                comm->cutoff_mbody = min(comm->cutoff_mbody, comm->cutoff);
+                comm->cutoff_mbody = std::min(comm->cutoff_mbody, comm->cutoff);
                 if (comm->cellsize_limit >= comm->cutoff)
                 {
                     /* We don't loose a lot of efficieny
@@ -7016,7 +7006,7 @@ gmx_domdec_t *init_domain_decomposition(FILE *fplog, t_commrec *cr,
                 }
             }
             /* Check if we did not end up below our original limit */
-            comm->cutoff_mbody = max(comm->cutoff_mbody, r_bonded_limit);
+            comm->cutoff_mbody = std::max(comm->cutoff_mbody, r_bonded_limit);
 
             if (comm->cutoff_mbody > comm->cellsize_limit)
             {
@@ -7079,7 +7069,7 @@ static void turn_on_dlb(FILE *fplog, t_commrec *cr, gmx_int64_t step)
     cellsize_min = comm->cellsize_min[dd->dim[0]];
     for (d = 1; d < dd->ndim; d++)
     {
-        cellsize_min = min(cellsize_min, comm->cellsize_min[dd->dim[d]]);
+        cellsize_min = std::min(cellsize_min, comm->cellsize_min[dd->dim[d]]);
     }
 
     if (cellsize_min < comm->cellsize_limit*1.05)
@@ -7144,8 +7134,6 @@ void dd_init_bondeds(FILE *fplog,
                      t_inputrec *ir, gmx_bool bBCheck, cginfo_mb_t *cginfo_mb)
 {
     gmx_domdec_comm_t *comm;
-    gmx_bool           bBondComm;
-    int                d;
 
     dd_make_reverse_top(fplog, dd, mtop, vsite, ir, bBCheck);
 
@@ -7256,7 +7244,7 @@ static void print_dd_settings(FILE *fplog, gmx_domdec_t *dd,
             limit = dd->comm->cellsize_min[XX];
             for (d = 1; d < DIM; d++)
             {
-                limit = min(limit, dd->comm->cellsize_min[d]);
+                limit = std::min(limit, dd->comm->cellsize_min[d]);
             }
         }
 
@@ -7264,10 +7252,10 @@ static void print_dd_settings(FILE *fplog, gmx_domdec_t *dd,
         {
             fprintf(fplog, "%40s  %-7s %6.3f nm\n",
                     "two-body bonded interactions", "(-rdd)",
-                    max(comm->cutoff, comm->cutoff_mbody));
+                    std::max(comm->cutoff, comm->cutoff_mbody));
             fprintf(fplog, "%40s  %-7s %6.3f nm\n",
                     "multi-body bonded interactions", "(-rdd)",
-                    (comm->bBondComm || dd->bGridJump) ? comm->cutoff_mbody : min(comm->cutoff, limit));
+                    (comm->bBondComm || dd->bGridJump) ? comm->cutoff_mbody : std::min(comm->cutoff, limit));
         }
         if (dd->vsite_comm)
         {
@@ -7302,7 +7290,7 @@ static void set_cell_limits_dlb(gmx_domdec_t      *dd,
 
     /* Determine the maximum number of comm. pulses in one dimension */
 
-    comm->cellsize_limit = max(comm->cellsize_limit, comm->cutoff_mbody);
+    comm->cellsize_limit = std::max(comm->cellsize_limit, comm->cutoff_mbody);
 
     /* Determine the maximum required number of grid pulses */
     if (comm->cellsize_limit >= comm->cutoff)
@@ -7323,7 +7311,7 @@ static void set_cell_limits_dlb(gmx_domdec_t      *dd,
     else
     {
         /* There is no cell size limit */
-        npulse = max(dd->nc[XX]-1, max(dd->nc[YY]-1, dd->nc[ZZ]-1));
+        npulse = std::max(dd->nc[XX]-1, std::max(dd->nc[YY]-1, dd->nc[ZZ]-1));
     }
 
     if (!bNoCutOff && npulse > 1)
@@ -7335,9 +7323,9 @@ static void set_cell_limits_dlb(gmx_domdec_t      *dd,
             dim      = dd->dim[d];
             npulse_d = (int)(1 + dd->nc[dim]*comm->cutoff
                              /(ddbox->box_size[dim]*ddbox->skew_fac[dim]*dlb_scale));
-            npulse_d_max = max(npulse_d_max, npulse_d);
+            npulse_d_max = std::max(npulse_d_max, npulse_d);
         }
-        npulse = min(npulse, npulse_d_max);
+        npulse = std::min(npulse, npulse_d_max);
     }
 
     /* This env var can override npulse */
@@ -7351,10 +7339,10 @@ static void set_cell_limits_dlb(gmx_domdec_t      *dd,
     comm->bVacDLBNoLimit = (ir->ePBC == epbcNONE);
     for (d = 0; d < dd->ndim; d++)
     {
-        comm->cd[d].np_dlb    = min(npulse, dd->nc[dd->dim[d]]-1);
+        comm->cd[d].np_dlb    = std::min(npulse, dd->nc[dd->dim[d]]-1);
         comm->cd[d].np_nalloc = comm->cd[d].np_dlb;
         snew(comm->cd[d].ind, comm->cd[d].np_nalloc);
-        comm->maxpulse = max(comm->maxpulse, comm->cd[d].np_dlb);
+        comm->maxpulse = std::max(comm->maxpulse, comm->cd[d].np_dlb);
         if (comm->cd[d].np_dlb < dd->nc[dd->dim[d]]-1)
         {
             comm->bVacDLBNoLimit = FALSE;
@@ -7364,10 +7352,10 @@ static void set_cell_limits_dlb(gmx_domdec_t      *dd,
     /* cellsize_limit is set for LINCS in init_domain_decomposition */
     if (!comm->bVacDLBNoLimit)
     {
-        comm->cellsize_limit = max(comm->cellsize_limit,
-                                   comm->cutoff/comm->maxpulse);
+        comm->cellsize_limit = std::max(comm->cellsize_limit,
+                                        comm->cutoff/comm->maxpulse);
     }
-    comm->cellsize_limit = max(comm->cellsize_limit, comm->cutoff_mbody);
+    comm->cellsize_limit = std::max(comm->cellsize_limit, comm->cutoff_mbody);
     /* Set the minimum cell size for each DD dimension */
     for (d = 0; d < dd->ndim; d++)
     {
@@ -7384,7 +7372,7 @@ static void set_cell_limits_dlb(gmx_domdec_t      *dd,
     }
     if (comm->cutoff_mbody <= 0)
     {
-        comm->cutoff_mbody = min(comm->cutoff, comm->cellsize_limit);
+        comm->cutoff_mbody = std::min(comm->cutoff, comm->cellsize_limit);
     }
     if (comm->bDynLoadBal)
     {
@@ -7475,7 +7463,7 @@ void set_dd_parameters(FILE *fplog, gmx_domdec_t *dd, real dlb_scale,
     }
     natoms_tot = comm->cgs_gl.index[comm->cgs_gl.nr];
 
-    dd->ga2la = ga2la_init(natoms_tot, vol_frac*natoms_tot);
+    dd->ga2la = ga2la_init(natoms_tot, static_cast<int>(vol_frac*natoms_tot));
 }
 
 static gmx_bool test_dd_cutoff(t_commrec *cr,
@@ -7739,11 +7727,11 @@ set_dd_corners(const gmx_domdec_t *dd,
         c->c[1][1] = comm->cell_x0[dim1];
         if (dd->bGridJump)
         {
-            c->c[1][1] = max(comm->cell_x0[dim1], comm->zone_d1[1].mch0);
+            c->c[1][1] = std::max(comm->cell_x0[dim1], comm->zone_d1[1].mch0);
             if (bDistMB)
             {
                 /* For the multi-body distance we need the maximum */
-                c->bc[1] = max(comm->cell_x0[dim1], comm->zone_d1[1].p1_0);
+                c->bc[1] = std::max(comm->cell_x0[dim1], comm->zone_d1[1].p1_0);
             }
         }
         /* Set the upper-right corner for rounding */
@@ -7766,8 +7754,8 @@ set_dd_corners(const gmx_domdec_t *dd,
                         if (j >= 4)
                         {
                             c->c[2][j-4] =
-                                max(c->c[2][j-4],
-                                    comm->zone_d2[zones->shift[i][dim0]][zones->shift[i][dim1]].mch0);
+                                std::max(c->c[2][j-4],
+                                         comm->zone_d2[zones->shift[i][dim0]][zones->shift[i][dim1]].mch0);
                         }
                     }
                 }
@@ -7779,7 +7767,7 @@ set_dd_corners(const gmx_domdec_t *dd,
                     {
                         for (j = 0; j < 2; j++)
                         {
-                            c->bc[2] = max(c->bc[2], comm->zone_d2[i][j].p1_0);
+                            c->bc[2] = std::max(c->bc[2], comm->zone_d2[i][j].p1_0);
                         }
                     }
                 }
@@ -7793,11 +7781,11 @@ set_dd_corners(const gmx_domdec_t *dd,
             c->cr1[3] = comm->cell_x1[dim1];
             if (dd->bGridJump)
             {
-                c->cr1[0] = max(comm->cell_x1[dim1], comm->zone_d1[1].mch1);
+                c->cr1[0] = std::max(comm->cell_x1[dim1], comm->zone_d1[1].mch1);
                 if (bDistMB)
                 {
                     /* For the multi-body distance we need the maximum */
-                    c->bcr1 = max(comm->cell_x1[dim1], comm->zone_d1[1].p1_1);
+                    c->bcr1 = std::max(comm->cell_x1[dim1], comm->zone_d1[1].p1_1);
                 }
             }
         }
@@ -8077,7 +8065,7 @@ static void setup_dd_communication(gmx_domdec_t *dd,
 {
     int                    dim_ind, dim, dim0, dim1, dim2, dimd, p, nat_tot;
     int                    nzone, nzone_send, zone, zonei, cg0, cg1;
-    int                    c, i, j, cg, cg_gl, nrcg;
+    int                    c, i, cg, cg_gl, nrcg;
     int                   *zone_cg_range, pos_cg, *index_gl, *cgindex, *recv_i;
     gmx_domdec_comm_t     *comm;
     gmx_domdec_zones_t    *zones;
@@ -8085,7 +8073,7 @@ static void setup_dd_communication(gmx_domdec_t *dd,
     gmx_domdec_ind_t      *ind;
     cginfo_mb_t           *cginfo_mb;
     gmx_bool               bBondComm, bDist2B, bDistMB, bDistBonded;
-    real                   r_mb, r_comm2, r_scomm2, r_bcomm2, r_0, r_1, r2inc, inv_ncg;
+    real                   r_comm2, r_bcomm2;
     dd_corners_t           corners;
     ivec                   tric_dist;
     rvec                  *cg_cm, *normal, *v_d, *v_0 = NULL, *v_1 = NULL, *recv_vr;
@@ -8116,8 +8104,6 @@ static void setup_dd_communication(gmx_domdec_t *dd,
 
     for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
     {
-        dim = dd->dim[dim_ind];
-
         /* Check if we need to use triclinic distances */
         tric_dist[dim_ind] = 0;
         for (i = 0; i <= dim_ind; i++)
@@ -8409,7 +8395,7 @@ static void setup_dd_communication(gmx_domdec_t *dd,
                     /* The rvec buffer is also required for atom buffers
                      * of size nrecv in dd_move_x and dd_move_f.
                      */
-                    i = max(cd->ind[0].nrecv[nzone+1], ind->nrecv[nzone+1]);
+                    i = std::max(cd->ind[0].nrecv[nzone+1], ind->nrecv[nzone+1]);
                     vec_rvec_check_alloc(&comm->vbuf2, i);
                 }
             }
@@ -8554,10 +8540,9 @@ static void set_zones_size(gmx_domdec_t *dd,
     gmx_domdec_comm_t  *comm;
     gmx_domdec_zones_t *zones;
     gmx_bool            bDistMB;
-    int                 z, zi, zj0, zj1, d, dim;
+    int                 z, zi, d, dim;
     real                rcs, rcmbs;
     int                 i, j;
-    real                size_j, add_tric;
     real                vol;
 
     comm = dd->comm;
@@ -8655,8 +8640,8 @@ static void set_zones_size(gmx_domdec_t *dd,
                                  * With multiple pulses this will lead
                                  * to a larger zone then strictly necessary.
                                  */
-                                zones->size[z].x1[dim] = max(zones->size[z].x1[dim],
-                                                             zones->size[zi].x1[dim]+rcmbs);
+                                zones->size[z].x1[dim] = std::max(zones->size[z].x1[dim],
+                                                                  zones->size[zi].x1[dim]+rcmbs);
                             }
                         }
                     }
@@ -8675,8 +8660,8 @@ static void set_zones_size(gmx_domdec_t *dd,
                 {
                     if (zones->shift[z][dim] > 0)
                     {
-                        zones->size[z].x1[dim] = max(zones->size[z].x1[dim],
-                                                     zones->size[zi].x1[dim]+rcs);
+                        zones->size[z].x1[dim] = std::max(zones->size[z].x1[dim],
+                                                          zones->size[zi].x1[dim]+rcs);
                     }
                 }
             }
@@ -8740,8 +8725,8 @@ static void set_zones_size(gmx_domdec_t *dd,
             {
                 for (i = 0; i < DIM; i++)
                 {
-                    corner_min[i] = min(corner_min[i], corner[i]);
-                    corner_max[i] = max(corner_max[i], corner[i]);
+                    corner_min[i] = std::min(corner_min[i], corner[i]);
+                    corner_max[i] = std::max(corner_max[i], corner[i]);
                 }
             }
         }
@@ -8911,8 +8896,7 @@ static int dd_sort_order(gmx_domdec_t *dd, t_forcerec *fr, int ncg_home_old)
 {
     gmx_domdec_sort_t *sort;
     gmx_cgsort_t      *cgsort, *sort_i;
-    int                ncg_new, nsort2, nsort_new, i, *a, moved, *ibuf;
-    int                sort_last, sort_skip;
+    int                ncg_new, nsort2, nsort_new, i, *a, moved;
 
     sort = dd->comm->sort;
 
@@ -9023,7 +9007,7 @@ static void dd_sort_state(gmx_domdec_t *dd, rvec *cgcm, t_forcerec *fr, t_state
                           int ncg_home_old)
 {
     gmx_domdec_sort_t *sort;
-    gmx_cgsort_t      *cgsort, *sort_i;
+    gmx_cgsort_t      *cgsort;
     int               *cgindex;
     int                ncg_new, i, *ibuf, cgsize;
     rvec              *vbuf;
@@ -9280,7 +9264,7 @@ void dd_partition_system(FILE                *fplog,
     t_block           *cgs_gl;
     gmx_int64_t        step_pcoupl;
     rvec               cell_ns_x0, cell_ns_x1;
-    int                i, j, n, ncgindex_set, ncg_home_old = -1, ncg_moved, nat_f_novirsum;
+    int                i, n, ncgindex_set, ncg_home_old = -1, ncg_moved, nat_f_novirsum;
     gmx_bool           bBoxChanged, bNStGlobalComm, bDoDLB, bCheckDLB, bTurnOnDLB, bLogLoad;
     gmx_bool           bRedist, bSortCG, bResortAll;
     ivec               ncells_old = {0, 0, 0}, ncells_new = {0, 0, 0}, np;
@@ -9348,7 +9332,7 @@ void dd_partition_system(FILE                *fplog,
              * and every 100 partitionings,
              * so the extra communication cost is negligible.
              */
-            n         = max(100, nstglobalcomm);
+            n         = std::max(100, nstglobalcomm);
             bCheckDLB = (comm->n_load_collect == 0 ||
                          comm->n_load_have % n == n-1);
         }
diff --git a/src/gromacs/mdlib/domdec_box.c b/src/gromacs/mdlib/domdec_box.cpp
similarity index 100%
rename from src/gromacs/mdlib/domdec_box.c
rename to src/gromacs/mdlib/domdec_box.cpp
diff --git a/src/gromacs/mdlib/domdec_con.c b/src/gromacs/mdlib/domdec_con.cpp
similarity index 97%
rename from src/gromacs/mdlib/domdec_con.c
rename to src/gromacs/mdlib/domdec_con.cpp
index e538a6a504..22cf39f781 100644
--- a/src/gromacs/mdlib/domdec_con.c
+++ b/src/gromacs/mdlib/domdec_con.cpp
@@ -35,6 +35,7 @@
 
 #include "gmxpre.h"
 
+#include <algorithm>
 #include <assert.h>
 
 #include "gromacs/legacyheaders/constr.h"
@@ -49,6 +50,7 @@
 #include "gromacs/pbcutil/ishift.h"
 #include "gromacs/topology/mtop_util.h"
 #include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/gmxassert.h"
 #include "gromacs/utility/smalloc.h"
 
 typedef struct {
@@ -467,8 +469,6 @@ void dd_clear_local_constraint_indices(gmx_domdec_t *dd)
 
 void dd_clear_local_vsite_indices(gmx_domdec_t *dd)
 {
-    int i;
-
     if (dd->vsite_comm)
     {
         gmx_hash_clear_and_optimize(dd->ga2la_vsite);
@@ -487,7 +487,7 @@ static int setup_specat_communication(gmx_domdec_t             *dd,
     int               nsend[2], nlast, nsend_zero[2] = {0, 0}, *nsend_ptr;
     int               d, dim, ndir, dir, nr, ns, i, nrecv_local, n0, start, indr, ind, buf[2];
     int               nat_tot_specat, nat_tot_prev, nalloc_old;
-    gmx_bool          bPBC, bFirst;
+    gmx_bool          bPBC;
     gmx_specatsend_t *spas;
 
     if (debug)
@@ -556,7 +556,6 @@ static int setup_specat_communication(gmx_domdec_t             *dd,
     nrecv_local    = 0;
     for (d = 0; d < dd->ndim; d++)
     {
-        bFirst = (d == 0);
         /* Pulse the grid forward and backward */
         if (dd->dim[d] >= dd->npbcdim || dd->nc[dd->dim[d]] > 2)
         {
@@ -1073,6 +1072,18 @@ int dd_make_local_constraints(gmx_domdec_t *dd, int at_start,
     int at_end, i, j;
     t_iatom                  *iap;
 
+    // This code should not be called unless this condition is true,
+    // because that's the only time init_domdec_constraints is
+    // called...
+    GMX_RELEASE_ASSERT(dd->bInterCGcons || dd->bInterCGsettles, "dd_make_local_constraints called when there are no local constraints");
+    // ... and init_domdec_constraints always sets
+    // dd->constraint_comm...
+    GMX_RELEASE_ASSERT(dd->constraint_comm, "Invalid use of dd_make_local_constraints before construction of constraint_comm");
+    // ... which static analysis needs to be reassured about, because
+    // otherwise, when dd->bInterCGsettles is
+    // true. dd->constraint_comm is unilaterally dereferenced before
+    // the call to atoms_to_settles.
+
     dc = dd->constraints;
 
     ilc_local = &il_local[F_CONSTR];
@@ -1088,6 +1099,7 @@ int dd_make_local_constraints(gmx_domdec_t *dd, int at_start,
     }
     else
     {
+        // Currently unreachable
         at2con_mt = NULL;
         ireq      = NULL;
     }
@@ -1243,6 +1255,7 @@ int dd_make_local_constraints(gmx_domdec_t *dd, int at_start,
     }
     else
     {
+        // Currently unreachable
         at_end = at_start;
     }
 
@@ -1254,7 +1267,7 @@ int dd_make_local_vsites(gmx_domdec_t *dd, int at_start, t_ilist *lil)
     gmx_domdec_specat_comm_t *spac;
     ind_req_t                *ireq;
     gmx_hash_t                ga2la_specat;
-    int  ftype, nral, i, j, gat, a;
+    int  ftype, nral, i, j, a;
     t_ilist                  *lilf;
     t_iatom                  *iatoms;
     int  at_end;
@@ -1347,7 +1360,7 @@ void init_domdec_constraints(gmx_domdec_t *dd,
 {
     gmx_domdec_constraints_t *dc;
     gmx_molblock_t           *molb;
-    int mb, ncon, c, a;
+    int mb, ncon, c;
 
     if (debug)
     {
@@ -1383,8 +1396,8 @@ void init_domdec_constraints(gmx_domdec_t *dd,
     /* Use a hash table for the global to local index.
      * The number of keys is a rough estimate, it will be optimized later.
      */
-    dc->ga2la = gmx_hash_init(min(mtop->natoms/20,
-                                  mtop->natoms/(2*dd->nnodes)));
+    dc->ga2la = gmx_hash_init(std::min(mtop->natoms/20,
+                                       mtop->natoms/(2*dd->nnodes)));
 
     dc->nthread = gmx_omp_nthreads_get(emntDomdec);
     snew(dc->ils, dc->nthread);
@@ -1394,9 +1407,6 @@ void init_domdec_constraints(gmx_domdec_t *dd,
 
 void init_domdec_vsites(gmx_domdec_t *dd, int n_intercg_vsite)
 {
-    int i;
-    gmx_domdec_constraints_t *dc;
-
     if (debug)
     {
         fprintf(debug, "Begin init_domdec_vsites\n");
@@ -1405,8 +1415,8 @@ void init_domdec_vsites(gmx_domdec_t *dd, int n_intercg_vsite)
     /* Use a hash table for the global to local index.
      * The number of keys is a rough estimate, it will be optimized later.
      */
-    dd->ga2la_vsite = gmx_hash_init(min(n_intercg_vsite/20,
-                                        n_intercg_vsite/(2*dd->nnodes)));
+    dd->ga2la_vsite = gmx_hash_init(std::min(n_intercg_vsite/20,
+                                             n_intercg_vsite/(2*dd->nnodes)));
 
     dd->vsite_comm = specat_comm_init(1);
 }
diff --git a/src/gromacs/mdlib/domdec_network.c b/src/gromacs/mdlib/domdec_network.cpp
similarity index 100%
rename from src/gromacs/mdlib/domdec_network.c
rename to src/gromacs/mdlib/domdec_network.cpp
diff --git a/src/gromacs/mdlib/domdec_setup.c b/src/gromacs/mdlib/domdec_setup.cpp
similarity index 98%
rename from src/gromacs/mdlib/domdec_setup.c
rename to src/gromacs/mdlib/domdec_setup.cpp
index 0903c32d14..1d03b26bea 100644
--- a/src/gromacs/mdlib/domdec_setup.c
+++ b/src/gromacs/mdlib/domdec_setup.cpp
@@ -36,6 +36,7 @@
 #include "gmxpre.h"
 
 #include <assert.h>
+#include <cmath>
 #include <stdio.h>
 
 #include "gromacs/legacyheaders/domdec.h"
@@ -125,8 +126,8 @@ static gmx_bool fits_pp_pme_perf(int nnodes, int npme, float ratio)
     sfree(div);
     sfree(mdiv);
 
-    npp_root3  = (int)(pow(nnodes-npme, 1.0/3.0) + 0.5);
-    npme_root2 = (int)(sqrt(npme) + 0.5);
+    npp_root3  = static_cast<int>(std::pow(nnodes-npme, 1.0/3.0) + 0.5);
+    npme_root2 = static_cast<int>(std::sqrt(static_cast<double>(npme)) + 0.5);
 
     /* The check below gives a reasonable division:
      * factor 5 allowed at 5 or more PP nodes,
@@ -156,8 +157,7 @@ static int guess_npme(FILE *fplog, gmx_mtop_t *mtop, t_inputrec *ir, matrix box,
                       int nnodes)
 {
     float      ratio;
-    int        npme, nkx, nky;
-    t_inputrec ir_try;
+    int        npme;
 
     ratio = pme_load_estimate(mtop, ir, box);
 
@@ -247,7 +247,7 @@ static int div_up(int n, int f)
 
 real comm_box_frac(ivec dd_nc, real cutoff, gmx_ddbox_t *ddbox)
 {
-    int  i, j, k, npp;
+    int  i, j, k;
     rvec bt, nw;
     real comm_vol;
 
@@ -257,13 +257,11 @@ real comm_box_frac(ivec dd_nc, real cutoff, gmx_ddbox_t *ddbox)
         nw[i] = dd_nc[i]*cutoff/bt[i];
     }
 
-    npp      = 1;
     comm_vol = 0;
     for (i = 0; i < DIM; i++)
     {
         if (dd_nc[i] > 1)
         {
-            npp      *= dd_nc[i];
             comm_vol += nw[i];
             for (j = i+1; j < DIM; j++)
             {
@@ -311,7 +309,7 @@ static float comm_cost_est(real limit, real cutoff,
                            int npme_tot, ivec nc)
 {
     ivec  npme = {1, 1, 1};
-    int   i, j, k, nk, overlap;
+    int   i, j, nk, overlap;
     rvec  bt;
     float comm_vol, comm_vol_xf, comm_pme, cost_pbcdx;
     /* This is the cost of a pbc_dx call relative to the cost
@@ -506,7 +504,7 @@ static void assign_factors(gmx_domdec_t *dd,
                            float pbcdxr, int npme,
                            int ndiv, int *div, int *mdiv, ivec ir_try, ivec opt)
 {
-    int   x, y, z, i;
+    int   x, y, i;
     float ce;
 
     if (ndiv == 0)
diff --git a/src/gromacs/mdlib/domdec_top.c b/src/gromacs/mdlib/domdec_top.cpp
similarity index 97%
rename from src/gromacs/mdlib/domdec_top.c
rename to src/gromacs/mdlib/domdec_top.cpp
index cf45ba3ff4..029904b6dc 100644
--- a/src/gromacs/mdlib/domdec_top.c
+++ b/src/gromacs/mdlib/domdec_top.cpp
@@ -55,6 +55,7 @@
 #include "gromacs/topology/topsort.h"
 #include "gromacs/utility/cstringutil.h"
 #include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/gmxassert.h"
 #include "gromacs/utility/smalloc.h"
 
 /* for dd_init_local_state */
@@ -147,7 +148,7 @@ static void print_missing_interactions_mb(FILE *fplog, t_commrec *cr,
                                           t_idef *idef)
 {
     int      nril_mol, *assigned, *gatindex;
-    int      ftype, ftype_j, nral, i, j_mol, j, k, a0, a0_mol, mol, a, a_gl;
+    int      ftype, ftype_j, nral, i, j_mol, j, a0, a0_mol, mol, a;
     int      nprint;
     t_ilist *il;
     t_iatom *ia;
@@ -333,7 +334,10 @@ void dd_print_missing_interactions(FILE *fplog, t_commrec *cr, int local_count,
 
     if (DDMASTER(dd))
     {
-        fprintf(fplog, "\nA list of missing interactions:\n");
+        if (fplog)
+        {
+            fprintf(fplog, "\nA list of missing interactions:\n");
+        }
         fprintf(stderr, "\nA list of missing interactions:\n");
         rest_global = dd->nbonded_global;
         rest_local  = local_count;
@@ -349,7 +353,6 @@ void dd_print_missing_interactions(FILE *fplog, t_commrec *cr, int local_count,
                 || (dd->reverse_top->bConstr && ftype == F_CONSTR)
                 || (dd->reverse_top->bSettle && ftype == F_SETTLE))
             {
-                nral = NRAL(ftype);
                 n    = gmx_mtop_ftype_count(err_top_global, ftype);
                 if (ftype == F_CONSTR)
                 {
@@ -360,7 +363,10 @@ void dd_print_missing_interactions(FILE *fplog, t_commrec *cr, int local_count,
                 {
                     sprintf(buf, "%20s of %6d missing %6d",
                             interaction_function[ftype].longname, n, -ndiff);
-                    fprintf(fplog, "%s\n", buf);
+                    if (fplog)
+                    {
+                        fprintf(fplog, "%s\n", buf);
+                    }
                     fprintf(stderr, "%s\n", buf);
                 }
                 rest_global -= n;
@@ -373,7 +379,10 @@ void dd_print_missing_interactions(FILE *fplog, t_commrec *cr, int local_count,
         {
             sprintf(buf, "%20s of %6d missing %6d", "exclusions",
                     rest_global, -ndiff);
-            fprintf(fplog, "%s\n", buf);
+            if (fplog)
+            {
+                fprintf(fplog, "%s\n", buf);
+            }
             fprintf(stderr, "%s\n", buf);
         }
     }
@@ -398,9 +407,6 @@ void dd_print_missing_interactions(FILE *fplog, t_commrec *cr, int local_count,
 static void global_atomnr_to_moltype_ind(gmx_reverse_top_t *rt, int i_gl,
                                          int *mb, int *mt, int *mol, int *i_mol)
 {
-    int molb;
-
-
     gmx_molblock_ind_t *mbi   = rt->mbi;
     int                 start = 0;
     int                 end   =  rt->nmolblock; /* exclusive */
@@ -434,7 +440,7 @@ static void global_atomnr_to_moltype_ind(gmx_reverse_top_t *rt, int i_gl,
 
 static int count_excls(t_block *cgs, t_blocka *excls, int *n_intercg_excl)
 {
-    int n, n_inter, cg, at0, at1, at, excl, atj;
+    int n, cg, at0, at1, at, excl, atj;
 
     n               = 0;
     *n_intercg_excl = 0;
@@ -488,7 +494,6 @@ static int low_make_reverse_ilist(t_ilist *il_mt, t_atom *atom,
             bVSite = (interaction_function[ftype].flags & IF_VSITE);
             nral   = NRAL(ftype);
             il     = &il_mt[ftype];
-            ia     = il->iatoms;
             for (i = 0; i < il->nr; i += 1+nral)
             {
                 ia = il->iatoms + i;
@@ -706,7 +711,7 @@ void dd_make_reverse_top(FILE *fplog,
                          gmx_vsite_t *vsite,
                          t_inputrec *ir, gmx_bool bBCheck)
 {
-    int             mb, n_recursive_vsite, nexcl, nexcl_icg, a;
+    int             mb, nexcl, nexcl_icg;
     gmx_molblock_t *molb;
     gmx_moltype_t  *molt;
 
@@ -785,6 +790,60 @@ void dd_make_reverse_top(FILE *fplog,
     }
 }
 
+/* This routine is very similar to add_ifunc, but vsites interactions
+ * have more work to do than other kinds of interactions, and the
+ * complex way nral (and thus vector contents) depends on ftype
+ * confuses static analysis tools unless we fuse the vsite
+ * atom-indexing organization code with the ifunc-adding code, so that
+ * they can see that nral is the same value. */
+static gmx_inline void
+add_ifunc_for_vsites(t_iatom *tiatoms, gmx_ga2la_t ga2la,
+                     int nral, gmx_bool bHomeA,
+                     int a, int a_gl, int a_mol,
+                     const t_iatom *iatoms,
+                     t_ilist *il)
+{
+    t_iatom *liatoms;
+
+    if (il->nr+1+nral > il->nalloc)
+    {
+        il->nalloc = over_alloc_large(il->nr+1+nral);
+        srenew(il->iatoms, il->nalloc);
+    }
+    liatoms = il->iatoms + il->nr;
+    il->nr += 1 + nral;
+
+    /* Copy the type */
+    tiatoms[0] = iatoms[0];
+
+    if (bHomeA)
+    {
+        /* We know the local index of the first atom */
+        tiatoms[1] = a;
+    }
+    else
+    {
+        /* Convert later in make_local_vsites */
+        tiatoms[1] = -a_gl - 1;
+    }
+
+    for (int k = 2; k < 1+nral; k++)
+    {
+        int ak_gl = a_gl + iatoms[k] - a_mol;
+        if (!ga2la_get_home(ga2la, ak_gl, &tiatoms[k]))
+        {
+            /* Copy the global index, convert later in make_local_vsites */
+            tiatoms[k] = -(ak_gl + 1);
+        }
+        // Note that ga2la_get_home always sets the third parameter if
+        // it returns TRUE
+    }
+    for (int k = 0; k < 1+nral; k++)
+    {
+        liatoms[k] = tiatoms[k];
+    }
+}
+
 static gmx_inline void add_ifunc(int nral, t_iatom *tiatoms, t_ilist *il)
 {
     t_iatom *liatoms;
@@ -891,36 +950,13 @@ static void add_vsite(gmx_ga2la_t ga2la, int *index, int *rtil,
                       t_iatom *iatoms,
                       t_idef *idef, int **vsite_pbc, int *vsite_pbc_nalloc)
 {
-    int     k, ak_gl, vsi, pbc_a_mol;
-    t_iatom tiatoms[1+MAXATOMLIST], *iatoms_r;
+    int     k, vsi, pbc_a_mol;
+    t_iatom tiatoms[1+MAXATOMLIST];
     int     j, ftype_r, nral_r;
 
-    /* Copy the type */
-    tiatoms[0] = iatoms[0];
-
-    if (bHomeA)
-    {
-        /* We know the local index of the first atom */
-        tiatoms[1] = a;
-    }
-    else
-    {
-        /* Convert later in make_local_vsites */
-        tiatoms[1] = -a_gl - 1;
-    }
-
-    for (k = 2; k < 1+nral; k++)
-    {
-        ak_gl = a_gl + iatoms[k] - a_mol;
-        if (!ga2la_get_home(ga2la, ak_gl, &tiatoms[k]))
-        {
-            /* Copy the global index, convert later in make_local_vsites */
-            tiatoms[k] = -(ak_gl + 1);
-        }
-    }
-
     /* Add this interaction to the local topology */
-    add_ifunc(nral, tiatoms, &idef->il[ftype]);
+    add_ifunc_for_vsites(tiatoms, ga2la, nral, bHomeA, a, a_gl, a_mol, iatoms, &idef->il[ftype]);
+
     if (vsite_pbc)
     {
         vsi = idef->il[ftype].nr/(1+nral) - 1;
@@ -1417,11 +1453,10 @@ static int make_exclusions_zone(gmx_domdec_t *dd, gmx_domdec_zones_t *zones,
                                 int iz,
                                 int cg_start, int cg_end)
 {
-    int             nizone, n, count, jla0, jla1, jla;
+    int             n, count, jla0, jla1, jla;
     int             cg, la0, la1, la, a_gl, mb, mt, mol, a_mol, j, aj_mol;
     const t_blocka *excls;
     gmx_ga2la_t     ga2la;
-    int             a_loc;
     int             cell;
 
     ga2la = dd->ga2la;
@@ -1672,7 +1707,6 @@ static int make_local_bondeds_excls(gmx_domdec_t *dd,
         {
             int       cg0t, cg1t;
             t_idef   *idef_t;
-            int       ftype;
             int     **vsite_pbc;
             int      *vsite_pbc_nalloc;
             t_blocka *excl_t;
@@ -1799,7 +1833,7 @@ void dd_make_local_top(gmx_domdec_t *dd, gmx_domdec_zones_t *zones,
                        gmx_vsite_t *vsite,
                        gmx_mtop_t *mtop, gmx_localtop_t *ltop)
 {
-    gmx_bool bUniqueExcl, bRCheckMB, bRCheck2B, bRCheckExcl;
+    gmx_bool bRCheckMB, bRCheck2B;
     real     rc = -1;
     ivec     rcheck;
     int      d, nexcl;
@@ -1814,7 +1848,6 @@ void dd_make_local_top(gmx_domdec_t *dd, gmx_domdec_zones_t *zones,
 
     bRCheckMB   = FALSE;
     bRCheck2B   = FALSE;
-    bRCheckExcl = FALSE;
 
     if (dd->reverse_top->bMultiCGmols)
     {
@@ -1853,17 +1886,6 @@ void dd_make_local_top(gmx_domdec_t *dd, gmx_domdec_zones_t *zones,
                         d, cellsize_min[d], d, rcheck[d], bRCheck2B);
             }
         }
-        if (dd->reverse_top->bExclRequired)
-        {
-            bRCheckExcl = bRCheck2B;
-        }
-        else
-        {
-            /* If we don't have forces on exclusions,
-             * we don't care about exclusions being assigned mulitple times.
-             */
-            bRCheckExcl = FALSE;
-        }
         if (bRCheckMB || bRCheck2B)
         {
             make_la2lc(dd);
@@ -1964,12 +1986,13 @@ void dd_init_local_state(gmx_domdec_t *dd,
 
 static void check_link(t_blocka *link, int cg_gl, int cg_gl_j)
 {
-    int      k, aj;
+    int      k;
     gmx_bool bFound;
 
     bFound = FALSE;
     for (k = link->index[cg_gl]; k < link->index[cg_gl+1]; k++)
     {
+        GMX_RELEASE_ASSERT(link->a, "Inconsistent NULL pointer while making charge-group links");
         if (link->a[k] == cg_gl_j)
         {
             bFound = TRUE;
@@ -1977,6 +2000,8 @@ static void check_link(t_blocka *link, int cg_gl, int cg_gl_j)
     }
     if (!bFound)
     {
+        GMX_RELEASE_ASSERT(link->a || link->index[cg_gl+1]+1 > link->nalloc_a,
+                           "Inconsistent allocation of link");
         /* Add this charge group link */
         if (link->index[cg_gl+1]+1 > link->nalloc_a)
         {
@@ -2286,7 +2311,7 @@ void dd_bonded_cg_distance(FILE *fplog,
                            real *r_2b, real *r_mb)
 {
     gmx_bool        bExclRequired;
-    int             mb, cg_offset, at_offset, *at2cg, mol;
+    int             mb, at_offset, *at2cg, mol;
     t_graph         graph;
     gmx_vsite_t    *vsite;
     gmx_molblock_t *molb;
@@ -2302,7 +2327,6 @@ void dd_bonded_cg_distance(FILE *fplog,
 
     *r_2b     = 0;
     *r_mb     = 0;
-    cg_offset = 0;
     at_offset = 0;
     for (mb = 0; mb < mtop->nmolblock; mb++)
     {
@@ -2310,7 +2334,6 @@ void dd_bonded_cg_distance(FILE *fplog,
         molt = &mtop->moltype[molb->type];
         if (molt->cgs.nr == 1 || molb->nmol == 0)
         {
-            cg_offset += molb->nmol*molt->cgs.nr;
             at_offset += molb->nmol*molt->atoms.nr;
         }
         else
@@ -2348,7 +2371,6 @@ void dd_bonded_cg_distance(FILE *fplog,
                     a_mb_2 = at_offset + amol_mb_2;
                 }
 
-                cg_offset += molt->cgs.nr;
                 at_offset += molt->atoms.nr;
             }
             sfree(cg_cm);