Fixed triclinic 1xNx1 domain decomposition

With the Verlet scheme, 1D triclinic domain decomposition along
dimension y produces incorrect bounding boxes for the non-bonded grid.
This led to a lot of missing non-bonded interactions, which quickly
crashes any simulation affected by this.
Fixes #1631.

Change-Id: I9bd1fc9d983be839e0c9a8e62d47f6cf17684a03

diff --git a/src/mdlib/domdec.c b/src/mdlib/domdec.c
index ac8627bb57e48a190706c46496df47c5410b3ede..175cae83a21f7472e8d6f05d5a96371bdb8385ec 100644 (file)
--- a/src/mdlib/domdec.c
+++ b/src/mdlib/domdec.c
@@ -8763,20 +8763,21 @@ static void set_zones_size(gmx_domdec_t *dd,
             {
                 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 triclinic x-y and rectangular y-z.
-                 * Shift y back, so it will later end up at 0.
-                 */
-                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++)
                 {
-                    corner[j] += corner[i]*box[i][j]/box[i][i];
+                    /* With 1D domain decomposition the cg's are not in
+                     * a triclinic box, but triclinic x-y and rectangular y/x-z.
+                     * So we should ignore the coupling for the non
+                     * domain-decomposed dimension of the pair x and y.
+                     */
+                    if (!(dd->ndim == 1 && ((dd->dim[0] == XX && j == YY) ||
+                                            (dd->dim[0] == YY && j == XX))))
+                    {
+                        corner[j] += corner[i]*box[i][j]/box[i][i];
+                    }
                 }
             }
             if (c == 0)