Fix g_mindist with pbc=XY
authorBerk Hess <hess@kth.se>
Fri, 20 Jun 2014 07:01:33 +0000 (09:01 +0200)
committerGerrit Code Review <gerrit@gerrit.gromacs.org>
Fri, 20 Jun 2014 14:27:05 +0000 (16:27 +0200)
Fixes #1189.

Change-Id: I7c8242df6f2644bebc49c13c86412f8bc1d37428

src/tools/gmx_mindist.c

index 94718f20108ef3def64c60649419e61c2dbc43d8..2116cf270e768c1b43ca8da97dfedb94750dcbe2 100644 (file)
 #include "rmpbc.h"
 #include "xtcio.h"
 #include "gmx_ana.h"
+#include "names.h"
 
 
-static void periodic_dist(matrix box, rvec x[], int n, atom_id index[],
+static void periodic_dist(int ePBC,
+                          matrix box, rvec x[], int n, atom_id index[],
                           real *rmin, real *rmax, int *min_ind)
 {
-#define NSHIFT 26
-    int  sx, sy, sz, i, j, s;
+#define NSHIFT_MAX 26
+    int  nsz, nshift, sx, sy, sz, i, j, s;
     real sqr_box, r2min, r2max, r2;
-    rvec shift[NSHIFT], d0, d;
+    rvec shift[NSHIFT_MAX], d0, d;
 
-    sqr_box = sqr(min(norm(box[XX]), min(norm(box[YY]), norm(box[ZZ]))));
+    sqr_box = min(norm2(box[XX]), norm2(box[YY]));
+    if (ePBC == epbcXYZ)
+    {
+        sqr_box = min(sqr_box, norm2(box[ZZ]));
+        nsz = 1;
+    }
+    else if (ePBC == epbcXY)
+    {
+        nsz = 0;
+    }
+    else
+    {
+        gmx_fatal(FARGS, "pbc = %s is not supported by g_mindist",
+                  epbc_names[ePBC]);
+        nsz = 0; /* Keep compilers quiet */
+    }
 
-    s = 0;
-    for (sz = -1; sz <= 1; sz++)
+    nshift = 0;
+    for (sz = -nsz; sz <= nsz; sz++)
     {
         for (sy = -1; sy <= 1; sy++)
         {
@@ -81,9 +98,10 @@ static void periodic_dist(matrix box, rvec x[], int n, atom_id index[],
                 {
                     for (i = 0; i < DIM; i++)
                     {
-                        shift[s][i] = sx*box[XX][i]+sy*box[YY][i]+sz*box[ZZ][i];
+                        shift[nshift][i] =
+                            sx*box[XX][i] + sy*box[YY][i] + sz*box[ZZ][i];
                     }
-                    s++;
+                    nshift++;
                 }
             }
         }
@@ -102,7 +120,7 @@ static void periodic_dist(matrix box, rvec x[], int n, atom_id index[],
             {
                 r2max = r2;
             }
-            for (s = 0; s < NSHIFT; s++)
+            for (s = 0; s < nshift; s++)
             {
                 rvec_add(d0, shift[s], d);
                 r2 = norm2(d);
@@ -164,7 +182,7 @@ static void periodic_mindist_plot(const char *trxfn, const char *outfn,
             gmx_rmpbc(gpbc, natoms, box, x);
         }
 
-        periodic_dist(box, x, n, index, &rmin, &rmax, ind_min);
+        periodic_dist(ePBC, box, x, n, index, &rmin, &rmax, ind_min);
         if (rmin < rmint)
         {
             rmint    = rmin;