Fix PME order and grid size checks

The PME grid size restrictions were not checked in grompp.
Now both grompp and mdrun check the grid size. Both checks are done
with the same minimum grid size independent of DD and OpenMP settings.
Renamed calc_grid to calcFftGrid and let it take the minimum grid size
into account.
Also corrected an assert that checked for pme_order>=4 iso 3.

Change-Id: I57b300f4a413ea2942fa671be67839be7ae16c39

diff --git a/src/gromacs/ewald/pme-internal.h b/src/gromacs/ewald/pme-internal.h
index 74273b17e4f20815932ade96ab57cd4128f260c1..4cf2a6a87352c1cf1820b4df36852cd9f5ef8943 100644 (file)
--- a/src/gromacs/ewald/pme-internal.h
+++ b/src/gromacs/ewald/pme-internal.h
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2016, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -341,7 +341,6 @@ typedef struct gmx_pme_t {
 void gmx_pme_check_restrictions(int pme_order,
                                 int nkx, int nky, int nkz,
                                 int nnodes_major,
-                                int nnodes_minor,
                                 gmx_bool bUseThreads,
                                 gmx_bool bFatal,
                                 gmx_bool *bValidSettings);

diff --git a/src/gromacs/ewald/pme-load-balancing.cpp b/src/gromacs/ewald/pme-load-balancing.cpp
index a559a454fe0cc25bad557eeacab7d781ae7cac2e..40f12fac3feb8e499aee377e335467e32c404dba 100644 (file)
--- a/src/gromacs/ewald/pme-load-balancing.cpp
+++ b/src/gromacs/ewald/pme-load-balancing.cpp
@@ -56,6 +56,7 @@
 #include "gromacs/domdec/domdec.h"
 #include "gromacs/domdec/domdec_network.h"
 #include "gromacs/domdec/domdec_struct.h"
+#include "gromacs/ewald/pme.h"
 #include "gromacs/fft/calcgrid.h"
 #include "gromacs/gmxlib/network.h"
 #include "gromacs/math/functions.h"
@@ -320,11 +321,12 @@ static gmx_bool pme_loadbal_increase_cutoff(pme_load_balancing_t *pme_lb,
 
         fac *= 1.01;
         clear_ivec(set->grid);
-        sp = calc_grid(nullptr, pme_lb->box_start,
-                       fac*pme_lb->setup[pme_lb->cur].spacing,
-                       &set->grid[XX],
-                       &set->grid[YY],
-                       &set->grid[ZZ]);
+        sp = calcFftGrid(nullptr, pme_lb->box_start,
+                         fac*pme_lb->setup[pme_lb->cur].spacing,
+                         minimalPmeGridSize(pme_order),
+                         &set->grid[XX],
+                         &set->grid[YY],
+                         &set->grid[ZZ]);
 
         /* As here we can't easily check if one of the PME ranks
          * uses threading, we do a conservative grid check.
@@ -333,7 +335,7 @@ static gmx_bool pme_loadbal_increase_cutoff(pme_load_balancing_t *pme_lb,
          */
         gmx_pme_check_restrictions(pme_order,
                                    set->grid[XX], set->grid[YY], set->grid[ZZ],
-                                   npmeranks_x, npmeranks_y,
+                                   npmeranks_x,
                                    TRUE,
                                    FALSE,
                                    &grid_ok);

diff --git a/src/gromacs/ewald/pme-spread.cpp b/src/gromacs/ewald/pme-spread.cpp
index 477a24937869d47172199d321c7be6acbe2d6e70..cb35e529a32863324776147d9e07c1dbea6ea2d8 100644 (file)
--- a/src/gromacs/ewald/pme-spread.cpp
+++ b/src/gromacs/ewald/pme-spread.cpp
@@ -275,7 +275,7 @@ static void make_bsplines(splinevec theta, splinevec dtheta, int order,
         if (bDoSplines || coefficient[ii] != 0.0)
         {
             xptr = fractx[ii];
-            assert(order >= 4 && order <= PME_ORDER_MAX);
+            assert(order >= 3 && order <= PME_ORDER_MAX);
             switch (order)
             {
                 case 4:  CALC_SPLINE(4);     break;

diff --git a/src/gromacs/ewald/pme.cpp b/src/gromacs/ewald/pme.cpp
index b66892b8f55c0b855f491f74d26eccf83397b390..b2f33aec543f694892450c1ad9655b3f581df8a3 100644 (file)
--- a/src/gromacs/ewald/pme.cpp
+++ b/src/gromacs/ewald/pme.cpp
@@ -419,10 +419,26 @@ destroy_overlap_comm(const pme_overlap_t *ol)
     sfree(ol->recvbuf);
 }
 
+int minimalPmeGridSize(int pmeOrder)
+{
+    /* The actual grid size limitations are:
+     *   serial:        >= pme_order
+     *   DD, no OpenMP: >= 2*(pme_order - 1)
+     *   DD, OpenMP:    >= pme_order + 1
+     * But we use the maximum for simplicity since in practice there is not
+     * much performance difference between pme_order and 2*(pme_order -1).
+     */
+    int minimalSize = 2*(pmeOrder - 1);
+
+    GMX_RELEASE_ASSERT(pmeOrder >= 3, "pmeOrder has to be >= 3");
+    GMX_RELEASE_ASSERT(minimalSize >= pmeOrder + 1, "The grid size should be >= pmeOrder + 1");
+
+    return minimalSize;
+}
+
 void gmx_pme_check_restrictions(int pme_order,
                                 int nkx, int nky, int nkz,
                                 int nnodes_major,
-                                int nnodes_minor,
                                 gmx_bool bUseThreads,
                                 gmx_bool bFatal,
                                 gmx_bool *bValidSettings)
@@ -438,17 +454,18 @@ void gmx_pme_check_restrictions(int pme_order,
                   pme_order, PME_ORDER_MAX);
     }
 
-    if (nkx <= pme_order*(nnodes_major > 1 ? 2 : 1) ||
-        nky <= pme_order*(nnodes_minor > 1 ? 2 : 1) ||
-        nkz <= pme_order)
+    const int minGridSize = minimalPmeGridSize(pme_order);
+    if (nkx < minGridSize ||
+        nky < minGridSize ||
+        nkz < minGridSize)
     {
         if (!bFatal)
         {
             *bValidSettings = FALSE;
             return;
         }
-        gmx_fatal(FARGS, "The PME grid sizes need to be larger than pme_order (%d) and for dimensions with domain decomposition larger than 2*pme_order",
-                  pme_order);
+        gmx_fatal(FARGS, "The PME grid sizes need to be >= 2*(pme_order-1) (%d)",
+                  minGridSize);
     }
 
     /* Check for a limitation of the (current) sum_fftgrid_dd code.
@@ -644,7 +661,6 @@ int gmx_pme_init(struct gmx_pme_t **pmedata,
     gmx_pme_check_restrictions(pme->pme_order,
                                pme->nkx, pme->nky, pme->nkz,
                                pme->nnodes_major,
-                               pme->nnodes_minor,
                                pme->bUseThreads,
                                TRUE,
                                nullptr);

diff --git a/src/gromacs/ewald/pme.h b/src/gromacs/ewald/pme.h
index 2712cc7e56bdbdc6fc8f381897c14642e6b55220..d13c5be1e37b33a75b2cb74fe0606b337f746f8f 100644 (file)
--- a/src/gromacs/ewald/pme.h
+++ b/src/gromacs/ewald/pme.h
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2016, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -66,6 +66,9 @@ enum {
     GMX_SUM_GRID_FORWARD, GMX_SUM_GRID_BACKWARD
 };
 
+/*! \brief Return the smallest allowed PME grid size for \p pmeOrder */
+int minimalPmeGridSize(int pmeOrder);
+
 /*! \brief Initialize \p pmedata
  *
  * Return value 0 indicates all well, non zero is an error code.

diff --git a/src/gromacs/fft/calcgrid.cpp b/src/gromacs/fft/calcgrid.cpp
index 5e826ada8957323999c84c090abdfaecce254ea4..4850ee2590beff7b0ab16a9b9837db3a13935fc7 100644 (file)
--- a/src/gromacs/fft/calcgrid.cpp
+++ b/src/gromacs/fft/calcgrid.cpp
@@ -59,8 +59,9 @@ const int grid_init[g_initNR] = { 6, 8, 10, 12, 14, 16, 20, 24, 25, 28, 32, 36,
 #define g_baseNR 14
 const int grid_base[g_baseNR] = { 45, 48, 50, 52, 54, 56, 60, 64, 70, 72, 75, 80, 81, 84 };
 
-real calc_grid(FILE *fp, const matrix box, real gr_sp,
-               int *nx, int *ny, int *nz)
+real calcFftGrid(FILE *fp,
+                 const matrix box, real gridSpacing, int minGridPointsPerDim,
+                 int *nx, int *ny, int *nz)
 {
     int  d, n[DIM];
     int  i;
@@ -69,9 +70,9 @@ real calc_grid(FILE *fp, const matrix box, real gr_sp,
     rvec spacing;
     real max_spacing;
 
-    if ((*nx <= 0 || *ny <= 0 || *nz <= 0) && gr_sp <= 0)
+    if ((*nx <= 0 || *ny <= 0 || *nz <= 0) && gridSpacing <= 0)
     {
-        gmx_fatal(FARGS, "invalid fourier grid spacing: %g", gr_sp);
+        gmx_fatal(FARGS, "invalid fourier grid spacing: %g", gridSpacing);
     }
 
     if (grid_base[g_baseNR-1] % 4 != 0)
@@ -119,7 +120,8 @@ real calc_grid(FILE *fp, const matrix box, real gr_sp,
     {
         if (n[d] <= 0)
         {
-            nmin = static_cast<int>(box_size[d]/gr_sp + 0.999);
+            nmin = static_cast<int>(box_size[d]/gridSpacing + 0.999);
+            nmin = std::max(nmin, minGridPointsPerDim);
 
             i = g_initNR - 1;
             if (grid_init[i] >= nmin)

diff --git a/src/gromacs/fft/calcgrid.h b/src/gromacs/fft/calcgrid.h
index 86eb76987240a3bdcd21edbe6f69e0fb780e9450..4da63261320d3196b30d4f7e62cf88270e057049 100644 (file)
--- a/src/gromacs/fft/calcgrid.h
+++ b/src/gromacs/fft/calcgrid.h
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2010,2014,2015,2016, by the GROMACS development team, led by
+ * Copyright (c) 2010,2014,2015,2016,2017, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -42,11 +42,12 @@
 #include "gromacs/math/vectypes.h"
 #include "gromacs/utility/real.h"
 
-real calc_grid(FILE *fp,
-               const matrix box, real gr_sp,
-               int *nx, int *ny, int *nz);
+real calcFftGrid(FILE *fp,
+                 const matrix box, real gridSpacing, int minGridPointsPerDim,
+                 int *nx, int *ny, int *nz);
 /* Sets the number of grid points for each zero n* to the first reasonable
- * number which gives a spacing equal to or smaller than gr_sp.
+ * number which gives a spacing equal to or smaller than gridSpacing
+ * and is >= minGridPointsPerDim.
  * Returns the maximum grid spacing.
  */
 

diff --git a/src/gromacs/gmxana/gmx_pme_error.cpp b/src/gromacs/gmxana/gmx_pme_error.cpp
index 46373cff09fe2490bf824f39fa9077af8d6a72d9..cfbfed333c45ebc43c18875a641025ed9c4d1a25 100644 (file)
--- a/src/gromacs/gmxana/gmx_pme_error.cpp
+++ b/src/gromacs/gmxana/gmx_pme_error.cpp
@@ -41,6 +41,7 @@
 #include <algorithm>
 
 #include "gromacs/commandline/pargs.h"
+#include "gromacs/ewald/pme.h"
 #include "gromacs/fft/calcgrid.h"
 #include "gromacs/fileio/checkpoint.h"
 #include "gromacs/fileio/tpxio.h"
@@ -1170,7 +1171,8 @@ int gmx_pme_error(int argc, char *argv[])
         info.nkx[0] = 0;
         info.nky[0] = 0;
         info.nkz[0] = 0;
-        calc_grid(stdout, state.box, info.fourier_sp[0], &(info.nkx[0]), &(info.nky[0]), &(info.nkz[0]));
+        calcFftGrid(stdout, state.box, info.fourier_sp[0], minimalPmeGridSize(info.pme_order[0]),
+                    &(info.nkx[0]), &(info.nky[0]), &(info.nkz[0]));
         if ( (ir->nkx != info.nkx[0]) || (ir->nky != info.nky[0]) || (ir->nkz != info.nkz[0]) )
         {
             gmx_fatal(FARGS, "Wrong fourierspacing %f nm, input file grid = %d x %d x %d, computed grid = %d x %d x %d",

diff --git a/src/gromacs/gmxana/gmx_tune_pme.cpp b/src/gromacs/gmxana/gmx_tune_pme.cpp
index ca8eae0e00e50f2c38f769230384b36d455814a3..f2087a4069f29e5e7d5fb237e96afbf584dea5ab 100644 (file)
--- a/src/gromacs/gmxana/gmx_tune_pme.cpp
+++ b/src/gromacs/gmxana/gmx_tune_pme.cpp
@@ -48,6 +48,7 @@
 #endif
 
 #include "gromacs/commandline/pargs.h"
+#include "gromacs/ewald/pme.h"
 #include "gromacs/fft/calcgrid.h"
 #include "gromacs/fileio/checkpoint.h"
 #include "gromacs/fileio/tpxio.h"
@@ -1078,7 +1079,8 @@ static void make_benchmark_tprs(
 
             /* Scale the Fourier grid spacing */
             ir->nkx = ir->nky = ir->nkz = 0;
-            calc_grid(nullptr, state.box, fourierspacing*fac, &ir->nkx, &ir->nky, &ir->nkz);
+            calcFftGrid(nullptr, state.box, fourierspacing*fac, minimalPmeGridSize(ir->pme_order),
+                        &ir->nkx, &ir->nky, &ir->nkz);
 
             /* Adjust other radii since various conditions need to be fulfilled */
             if (eelPME == ir->coulombtype)

diff --git a/src/gromacs/gmxpreprocess/grompp.cpp b/src/gromacs/gmxpreprocess/grompp.cpp
index e088d43be7435c84575bce79be680a5cd2683a41..023e3808ea2693a6164737ae82c49398f87f2619 100644 (file)
--- a/src/gromacs/gmxpreprocess/grompp.cpp
+++ b/src/gromacs/gmxpreprocess/grompp.cpp
@@ -49,6 +49,7 @@
 #include <sys/types.h>
 
 #include "gromacs/commandline/pargs.h"
+#include "gromacs/ewald/pme.h"
 #include "gromacs/fft/calcgrid.h"
 #include "gromacs/fileio/confio.h"
 #include "gromacs/fileio/enxio.h"
@@ -2076,8 +2077,15 @@ int gmx_grompp(int argc, char *argv[])
             set_warning_line(wi, mdparin, -1);
             warning_error(wi, "Some of the Fourier grid sizes are set, but all of them need to be set.");
         }
-        calc_grid(stdout, box, ir->fourier_spacing,
-                  &(ir->nkx), &(ir->nky), &(ir->nkz));
+        const int minGridSize = minimalPmeGridSize(ir->pme_order);
+        calcFftGrid(stdout, box, ir->fourier_spacing, minGridSize,
+                    &(ir->nkx), &(ir->nky), &(ir->nkz));
+        if (ir->nkx < minGridSize ||
+            ir->nky < minGridSize ||
+            ir->nkz < minGridSize)
+        {
+            warning_error(wi, "The PME grid size should be >= 2*(pme-order - 1); either manually increase the grid size or decrease pme-order");
+        }
     }
 
     /* MRS: eventually figure out better logic for initializing the fep