fixed segv in Verlet pair-search with trilinic DD

With the Verlet scheme and domain decomposition the bounding box
of the zone, used for the pair-search grid, could be calculated
incorrectly for certain triclinic boxes, causing array indices
to get out of bounds (usually causing a segv).
Fixes #1059

Change-Id: I5bb146baaec59271dd5e68b40c545ed29d061831

diff --git a/src/mdlib/domdec.c b/src/mdlib/domdec.c
index 74940478d76b655db9283cff43fd75b7b9864c03..2f19d21b3b39fbdd665e5720ac333fa245dcffc9 100644 (file)
--- a/src/mdlib/domdec.c
+++ b/src/mdlib/domdec.c
@@ -8448,39 +8448,74 @@ static void set_zones_size(gmx_domdec_t *dd,
 
     for(z=zone_start; z<zone_end; z++)
     {
-        for(i=0; i<DIM; i++)
-        {
-            zones->size[z].bb_x0[i] = zones->size[z].x0[i];
-            zones->size[z].bb_x1[i] = zones->size[z].x1[i];
+        /* Initialization only required to keep the compiler happy */
+        rvec corner_min={0,0,0},corner_max={0,0,0},corner;
+        int  nc,c;
 
-            for(j=i+1; j<ddbox->npbcdim; j++)
+        /* To determine the bounding box for a zone we need to find
+         * the extreme corners of 4, 2 or 1 corners.
+         */
+        nc = 1 << (ddbox->npbcdim - 1);
+
+        for(c=0; c<nc; c++)
+        {
+            /* Set up a zone corner at x=0, ignoring trilinic couplings */
+            corner[XX] = 0;
+            if ((c & 1) == 0)
+            {
+                corner[YY] = zones->size[z].x0[YY];
+            }
+            else
+            {
+                corner[YY] = zones->size[z].x1[YY];
+            }
+            if ((c & 2) == 0)
+            {
+                corner[ZZ] = zones->size[z].x0[ZZ];
+            }
+            else
+            {
+                corner[ZZ] = zones->size[z].x1[ZZ];
+            }
+            if (dd->ndim == 1 && box[ZZ][YY] != 0)
             {
                 /* With 1D domain decomposition the cg's are not in
-                 * the triclinic box, but trilinic x-y and rectangular y-z.
+                 * the triclinic box, but triclinic x-y and rectangular y-z.
+                 * Shift y back, so it will later end up at 0.
                  */
-                if (box[j][i] != 0 &&
-                    !(dd->ndim == 1 && i == YY && j == ZZ))
+                corner[YY] -= corner[ZZ]*box[ZZ][YY]/box[ZZ][ZZ];
+            }
+            /* Apply the triclinic couplings */
+            for(i=YY; i<ddbox->npbcdim; i++)
+            {
+                for(j=XX; j<i; j++)
                 {
-                    /* Correct for triclinic offset of the lower corner */
-                    add_tric = zones->size[z].x0[j]*box[j][i]/box[j][j];
-                    zones->size[z].bb_x0[i] += add_tric;
-                    zones->size[z].bb_x1[i] += add_tric;
-
-                    /* Correct for triclinic offset of the upper corner */
-                    size_j = zones->size[z].x1[j] - zones->size[z].x0[j];
-                    add_tric = size_j*box[j][i]/box[j][j];
-
-                    if (box[j][i] < 0)
-                    {
-                        zones->size[z].bb_x0[i] += add_tric;
-                    }
-                    else
-                    {
-                        zones->size[z].bb_x1[i] += add_tric;
-                    }
+                    corner[j] += corner[i]*box[i][j]/box[i][i];
                 }
             }
+            if (c == 0)
+            {
+                copy_rvec(corner,corner_min);
+                copy_rvec(corner,corner_max);
+            }
+            else
+            {
+                for(i=0; i<DIM; i++)
+                {
+                    corner_min[i] = min(corner_min[i],corner[i]);
+                    corner_max[i] = max(corner_max[i],corner[i]);
+                }
+            }
+        }
+        /* Copy the extreme cornes without offset along x */
+        for(i=0; i<DIM; i++)
+        {
+            zones->size[z].bb_x0[i] = corner_min[i];
+            zones->size[z].bb_x1[i] = corner_max[i];
         }
+        /* Add the offset along x */
+        zones->size[z].bb_x0[XX] += zones->size[z].x0[XX];
+        zones->size[z].bb_x1[XX] += zones->size[z].x1[XX];
     }
 
     if (zone_start == 0)

diff --git a/src/mdlib/nbnxn_search.c b/src/mdlib/nbnxn_search.c
index f674699423a1717b126f6cbc70964fdd4b8a34dd..493025b07cf1c7dfabeb5bf5a47fa637cf6645c5 100644 (file)
--- a/src/mdlib/nbnxn_search.c
+++ b/src/mdlib/nbnxn_search.c
@@ -1381,8 +1381,10 @@ static void calc_column_indices(nbnxn_grid_t *grid,
             if (cell[i] < 0 || cell[i] >= grid->ncx*grid->ncy)
             {
                 gmx_fatal(FARGS,
-                          "grid cell cx %d cy %d out of range (max %d %d)",
-                          cx,cy,grid->ncx,grid->ncy);
+                          "grid cell cx %d cy %d out of range (max %d %d)\n"
+                          "atom %f %f %f, grid->c0 %f %f",
+                          cx,cy,grid->ncx,grid->ncy,
+                          x[i][XX],x[i][YY],x[i][ZZ],grid->c0[XX],grid->c0[YY]);
             }
 #endif
         }