Fix recent bug with trilinic 1D DD

A recent bug-fix (c8d919a3) for triclinic 1D domain decompostion
introduced a bug for boxes with box[YY][XX]!=0.
Fixes #1656.
Refs #1631.

Change-Id: I06b9376212390b73e90a3ce9704dee2bad9693fb

diff --git a/src/mdlib/domdec.c b/src/mdlib/domdec.c
index 175cae83a21f7472e8d6f05d5a96371bdb8385ec..ca21e5da3d45bf504ce7cbadee421d0e60321f41 100644 (file)
--- a/src/mdlib/domdec.c
+++ b/src/mdlib/domdec.c
@@ -8763,21 +8763,26 @@ static void set_zones_size(gmx_domdec_t *dd,
             {
                 corner[ZZ] = zones->size[z].x1[ZZ];
             }
+            if (dd->ndim == 1 && dd->dim[0] < ZZ && ZZ < dd->npbcdim &&
+                box[ZZ][1 - dd->dim[0]] != 0)
+            {
+                /* With 1D domain decomposition the cg's are not in
+                 * the triclinic box, but triclinic x-y and rectangular y/x-z.
+                 * Shift the corner of the z-vector back to along the box
+                 * vector of dimension d, so it will later end up at 0 along d.
+                 * This can affect the location of this corner along dd->dim[0]
+                 * through the matrix operation below if box[d][dd->dim[0]]!=0.
+                 */
+                int d = 1 - dd->dim[0];
+
+                corner[d] -= corner[ZZ]*box[ZZ][d]/box[ZZ][ZZ];
+            }
             /* Apply the triclinic couplings */
             for (i = YY; i < ddbox->npbcdim; i++)
             {
                 for (j = XX; j < i; j++)
                 {
-                    /* 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];
-                    }
+                    corner[j] += corner[i]*box[i][j]/box[i][i];
                 }
             }
             if (c == 0)