Fix harmless bug with combination of group and twin-range

The combination of
* group cut-off scheme,
* coulomb-modifier != none,
* rcoulomb > rlist (ie. twin-range), and
* rlistlong == -1

was being treated by the wrong clause in readir.c because the check
for "might be zero at cutoff" did not consider the modifier. The only
behavioural difference was issuing a warning that there will be no
buffering.

The modifier checks have been refactored to make it possible to check
the modifier.

Also removed some useless checks for EVDW_SWITCH

Change-Id: Ic64f3f48be17f6d87230bc4060a31e1fb4343973

diff --git a/src/gromacs/gmxana/gmx_tune_pme.c b/src/gromacs/gmxana/gmx_tune_pme.c
index 121e763e6385dcf175305e22bdbdca97ca3fa979..e63be167e18cbb6dff8c15422c00042d0a0f44ab 100644 (file)
--- a/src/gromacs/gmxana/gmx_tune_pme.c
+++ b/src/gromacs/gmxana/gmx_tune_pme.c
@@ -58,6 +58,7 @@
 #include "gmx_ana.h"
 #include "names.h"
 #include "perf_est.h"
+#include "inputrec.h"
 #include "gromacs/timing/walltime_accounting.h"
 
 
@@ -930,7 +931,7 @@ static void make_benchmark_tprs(
     fprintf(fp, "   Grid spacing x y z   : %f %f %f\n",
             box_size[XX]/ir->nkx, box_size[YY]/ir->nky, box_size[ZZ]/ir->nkz);
     fprintf(fp, "   Van der Waals type   : %s\n", EVDWTYPE(ir->vdwtype));
-    if (EVDW_SWITCHED(ir->vdwtype))
+    if (ir_vdw_switched(ir))
     {
         fprintf(fp, "   rvdw_switch          : %f nm\n", ir->rvdw_switch);
     }

diff --git a/src/gromacs/gmxlib/inputrec.c b/src/gromacs/gmxlib/inputrec.c
index 7514fad18e457743e7a4c193e6553dc61933df0d..74e186cc4011e0a0c96ada8bd581d6e67ba3425e 100644 (file)
--- a/src/gromacs/gmxlib/inputrec.c
+++ b/src/gromacs/gmxlib/inputrec.c
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2010, The GROMACS development team.
- * Copyright (c) 2012, by the GROMACS development team, led by
+ * Copyright (c) 2012,2014, 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.
@@ -195,3 +195,42 @@ int ir_optimal_nstpcouple(const t_inputrec *ir)
 
     return n;
 }
+
+gmx_bool ir_coulomb_switched(const t_inputrec *ir)
+{
+    return (ir->coulombtype == eelSWITCH ||
+            ir->coulombtype == eelSHIFT ||
+            ir->coulombtype == eelENCADSHIFT ||
+            ir->coulombtype == eelPMESWITCH ||
+            ir->coulombtype == eelPMEUSERSWITCH ||
+            ir->coulomb_modifier == eintmodPOTSWITCH);
+}
+
+gmx_bool ir_coulomb_is_zero_at_cutoff(const t_inputrec *ir)
+{
+    return (ir_coulomb_switched(ir) || ir->coulomb_modifier != eintmodNONE ||
+            ir->coulombtype == eelRF_ZERO);
+}
+
+gmx_bool ir_coulomb_might_be_zero_at_cutoff(const t_inputrec *ir)
+{
+    return (ir_coulomb_is_zero_at_cutoff(ir) || ir->coulombtype == eelUSER || ir->coulombtype == eelPMEUSER);
+}
+
+gmx_bool ir_vdw_switched(const t_inputrec *ir)
+{
+    return (ir->vdwtype == evdwSWITCH ||
+            ir->vdwtype == evdwSHIFT ||
+            ir->vdwtype == evdwENCADSHIFT ||
+            ir->vdw_modifier == eintmodPOTSWITCH);
+}
+
+gmx_bool ir_vdw_is_zero_at_cutoff(const t_inputrec *ir)
+{
+    return (ir_vdw_switched(ir) || ir->vdw_modifier != eintmodNONE);
+}
+
+gmx_bool ir_vdw_might_be_zero_at_cutoff(const t_inputrec *ir)
+{
+    return (ir_vdw_is_zero_at_cutoff(ir) || ir->vdwtype == evdwUSER);
+}

diff --git a/src/gromacs/gmxpreprocess/readir.c b/src/gromacs/gmxpreprocess/readir.c
index bc5b57e7f4864d5189b2a7d6d5990dee5bbd8d19..4fde014d5f7e50cbd7c40b080afd6b727f9a8d2d 100644 (file)
--- a/src/gromacs/gmxpreprocess/readir.c
+++ b/src/gromacs/gmxpreprocess/readir.c
@@ -274,8 +274,8 @@ void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
     {
         /* BASIC CUT-OFF STUFF */
         if (ir->rlist == 0 ||
-            !((EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) && ir->rcoulomb > ir->rlist) ||
-              (EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype)    && ir->rvdw     > ir->rlist)))
+            !((ir_coulomb_might_be_zero_at_cutoff(ir) && ir->rcoulomb > ir->rlist) ||
+              (ir_vdw_might_be_zero_at_cutoff(ir)     && ir->rvdw     > ir->rlist)))
         {
             /* No switched potential and/or no twin-range:
              * we can set the long-range cut-off to the maximum of the other cut-offs.
@@ -1076,9 +1076,9 @@ void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
      * means the interaction is zero outside rcoulomb, but it helps to
      * provide accurate energy conservation.
      */
-    if (EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype))
+    if (ir_coulomb_might_be_zero_at_cutoff(ir))
     {
-        if (EEL_SWITCHED(ir->coulombtype))
+        if (ir_coulomb_switched(ir))
         {
             sprintf(err_buf,
                     "With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the potential modifier options!",
@@ -1147,7 +1147,7 @@ void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
 
     if (EVDW_PME(ir->vdwtype))
     {
-        if (EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype))
+        if (ir_vdw_might_be_zero_at_cutoff(ir))
         {
             sprintf(err_buf, "With vdwtype = %s, rvdw must be <= rlist",
                     evdw_names[ir->vdwtype]);
@@ -1183,7 +1183,7 @@ void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
         CHECK(ir->wall_ewald_zfac < 2);
     }
 
-    if (EVDW_SWITCHED(ir->vdwtype))
+    if (ir_vdw_switched(ir))
     {
         sprintf(err_buf, "With vdwtype = %s rvdw-switch must be < rvdw. Or, better - use a potential modifier.",
                 evdw_names[ir->vdwtype]);
@@ -1209,14 +1209,13 @@ void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
                          "between neighborsearch steps");
         }
 
-        if (EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype)
-            && (ir->rlistlong <= ir->rcoulomb))
+        if (ir_coulomb_is_zero_at_cutoff(ir) && ir->rlistlong <= ir->rcoulomb)
         {
             sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rcoulomb.",
                     IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
             warning_note(wi, warn_buf);
         }
-        if (EVDW_SWITCHED(ir->vdwtype) && (ir->rlistlong <= ir->rvdw))
+        if (ir_vdw_switched(ir) && (ir->rlistlong <= ir->rvdw))
         {
             sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rvdw.",
                     IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
@@ -1251,13 +1250,13 @@ void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
     /* ENERGY CONSERVATION */
     if (ir_NVE(ir) && ir->cutoff_scheme == ecutsGROUP)
     {
-        if (!EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype) && ir->rvdw > 0 && ir->vdw_modifier == eintmodNONE)
+        if (!ir_vdw_might_be_zero_at_cutoff(ir) && ir->rvdw > 0 && ir->vdw_modifier == eintmodNONE)
         {
             sprintf(warn_buf, "You are using a cut-off for VdW interactions with NVE, for good energy conservation use vdwtype = %s (possibly with DispCorr)",
                     evdw_names[evdwSHIFT]);
             warning_note(wi, warn_buf);
         }
-        if (!EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) && ir->rcoulomb > 0 && ir->coulomb_modifier == eintmodNONE)
+        if (!ir_coulomb_might_be_zero_at_cutoff(ir) && ir->rcoulomb > 0 && ir->coulomb_modifier == eintmodNONE)
         {
             sprintf(warn_buf, "You are using a cut-off for electrostatics with NVE, for good energy conservation use coulombtype = %s or %s",
                     eel_names[eelPMESWITCH], eel_names[eelRF_ZERO]);
@@ -4207,8 +4206,7 @@ void check_chargegroup_radii(const gmx_mtop_t *mtop, const t_inputrec *ir,
              * since user defined interactions might purposely
              * not be zero at the cut-off.
              */
-            if ((EVDW_IS_ZERO_AT_CUTOFF(ir->vdwtype) ||
-                 ir->vdw_modifier != eintmodNONE) &&
+            if ((ir_vdw_is_zero_at_cutoff(ir) || ir->vdw_modifier != eintmodNONE) &&
                 rvdw1 + rvdw2 > ir->rlistlong - ir->rvdw)
             {
                 sprintf(warn_buf, "The sum of the two largest charge group "
@@ -4229,7 +4227,7 @@ void check_chargegroup_radii(const gmx_mtop_t *mtop, const t_inputrec *ir,
                     warning_note(wi, warn_buf);
                 }
             }
-            if ((EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype) ||
+            if ((ir_coulomb_is_zero_at_cutoff(ir) ||
                  ir->coulomb_modifier != eintmodNONE) &&
                 rcoul1 + rcoul2 > ir->rlistlong - ir->rcoulomb)
             {

diff --git a/src/gromacs/legacyheaders/inputrec.h b/src/gromacs/legacyheaders/inputrec.h
index f3d0e33874fadecfc3d2b9ed3f8737d53a0d85ba..9bebf01934974388e030ef120fca7c38dc84943c 100644 (file)
--- a/src/gromacs/legacyheaders/inputrec.h
+++ b/src/gromacs/legacyheaders/inputrec.h
@@ -57,6 +57,32 @@ int pcouple_min_integration_steps(int epc);
 
 int ir_optimal_nstpcouple(const t_inputrec *ir);
 
+/* Returns if the Coulomb force or potential is switched to zero */
+gmx_bool ir_coulomb_switched(const t_inputrec *ir);
+
+/* Returns if the Coulomb interactions are zero beyond the rcoulomb.
+ * Note: always returns TRUE for the Verlet cut-off scheme.
+ */
+gmx_bool ir_coulomb_is_zero_at_cutoff(const t_inputrec *ir);
+
+/* As ir_coulomb_is_zero_at_cutoff, but also returns TRUE for user tabulated
+ * interactions, since these might be zero beyond rcoulomb.
+ */
+gmx_bool ir_coulomb_might_be_zero_at_cutoff(const t_inputrec *ir);
+
+/* Returns if the Van der Waals force or potential is switched to zero */
+gmx_bool ir_vdw_switched(const t_inputrec *ir);
+
+/* Returns if the Van der Waals interactions are zero beyond the rvdw.
+ * Note: always returns TRUE for the Verlet cut-off scheme.
+ */
+gmx_bool ir_vdw_is_zero_at_cutoff(const t_inputrec *ir);
+
+/* As ir_vdw_is_zero_at_cutoff, but also returns TRUE for user tabulated
+ * interactions, since these might be zero beyond rvdw.
+ */
+gmx_bool ir_vdw_might_be_zero_at_cutoff(const t_inputrec *ir);
+
 #ifdef __cplusplus
 }
 #endif

diff --git a/src/gromacs/legacyheaders/types/enums.h b/src/gromacs/legacyheaders/types/enums.h
index 16ca1b92a18a40d600d29ef5cf8bb3c1f7af0e47..a77e714e365dc1efb813cfb56fb8927f5431f30f 100644 (file)
--- a/src/gromacs/legacyheaders/types/enums.h
+++ b/src/gromacs/legacyheaders/types/enums.h
@@ -114,15 +114,8 @@ enum {
 #define EEL_PME(e)  ((e) == eelPME || (e) == eelPMESWITCH || (e) == eelPMEUSER || (e) == eelPMEUSERSWITCH || (e) == eelP3M_AD)
 #define EEL_EWALD(e)  (EEL_PME(e) || (e) == eelEWALD)
 #define EEL_FULL(e) (EEL_PME(e) || (e) == eelPOISSON || (e) == eelEWALD)
-
-#define EEL_SWITCHED(e) ((e) == eelSWITCH || (e) == eelSHIFT || (e) == eelENCADSHIFT || (e) == eelPMESWITCH || (e) == eelPMEUSERSWITCH)
-
 #define EEL_USER(e) ((e) == eelUSER || (e) == eelPMEUSER || (e) == (eelPMEUSERSWITCH))
 
-#define EEL_IS_ZERO_AT_CUTOFF(e) (EEL_SWITCHED(e) || (e) == eelRF_ZERO)
-
-#define EEL_MIGHT_BE_ZERO_AT_CUTOFF(e) (EEL_IS_ZERO_AT_CUTOFF(e) || (e) == eelUSER || (e) == eelPMEUSER)
-
 enum {
     evdwCUT, evdwSWITCH, evdwSHIFT, evdwUSER, evdwENCADSHIFT,
     evdwPME, evdwNR
@@ -134,12 +127,6 @@ enum {
 
 #define EVDW_PME(e) ((e) == evdwPME)
 
-#define EVDW_SWITCHED(e) ((e) == evdwSWITCH || (e) == evdwSHIFT || (e) == evdwENCADSHIFT)
-
-#define EVDW_IS_ZERO_AT_CUTOFF(e) EVDW_SWITCHED(e)
-
-#define EVDW_MIGHT_BE_ZERO_AT_CUTOFF(e) (EVDW_IS_ZERO_AT_CUTOFF(e) || (e) == evdwUSER)
-
 enum {
     ensGRID, ensSIMPLE, ensNR
 };

diff --git a/src/gromacs/mdlib/ns.c b/src/gromacs/mdlib/ns.c
index f776001169b312ecdd338545bc1be4dafe995bdc..c54aef43b1b7d00518bf1555fd5a4f266d3e1cc6 100644 (file)
--- a/src/gromacs/mdlib/ns.c
+++ b/src/gromacs/mdlib/ns.c
@@ -2046,26 +2046,36 @@ static void get_cutoff2(t_forcerec *fr, gmx_bool bDoLongRange,
 
     if (bDoLongRange && fr->bTwinRange)
     {
-        /* The VdW and elec. LR cut-off's could be different,
+        /* With plain cut-off or RF we need to make the list exactly
+         * up to the cut-off and the cut-off's can be different,
          * so we can not simply set them to rlistlong.
+         * To keep this code compatible with (exotic) old cases,
+         * we also create lists up to rvdw/rcoulomb for PME and Ewald.
+         * The interaction check should correspond to:
+         * !ir_vdw/coulomb_might_be_zero_at_cutoff from inputrec.c.
          */
-        if (EVDW_MIGHT_BE_ZERO_AT_CUTOFF(fr->vdwtype) &&
-            fr->rvdw > fr->rlist)
+        if (((fr->vdwtype == evdwCUT || fr->vdwtype == evdwPME) &&
+             fr->vdw_modifier == eintmodNONE) ||
+            fr->rvdw <= fr->rlist)
         {
-            *rvdw2  = sqr(fr->rlistlong);
+            *rvdw2  = sqr(fr->rvdw);
         }
         else
         {
-            *rvdw2  = sqr(fr->rvdw);
+            *rvdw2  = sqr(fr->rlistlong);
         }
-        if (EEL_MIGHT_BE_ZERO_AT_CUTOFF(fr->eeltype) &&
-            fr->rcoulomb > fr->rlist)
+        if (((fr->eeltype == eelCUT ||
+              (EEL_RF(fr->eeltype) && fr->eeltype != eelRF_ZERO) ||
+              fr->eeltype == eelPME ||
+              fr->eeltype == eelEWALD) &&
+             fr->coulomb_modifier == eintmodNONE) ||
+            fr->rcoulomb <= fr->rlist)
         {
-            *rcoul2 = sqr(fr->rlistlong);
+            *rcoul2 = sqr(fr->rcoulomb);
         }
         else
         {
-            *rcoul2 = sqr(fr->rcoulomb);
+            *rcoul2 = sqr(fr->rlistlong);
         }
     }
     else

diff --git a/src/gromacs/mdlib/sim_util.c b/src/gromacs/mdlib/sim_util.c
index 6b6cf4cb4b1120609f5d3a5218d55accfe517d04..2aacb9b32c46d6ef65edd0d097927c41cdb526b0 100644 (file)
--- a/src/gromacs/mdlib/sim_util.c
+++ b/src/gromacs/mdlib/sim_util.c
@@ -2228,13 +2228,6 @@ void calc_enervirdiff(FILE *fplog, int eDispCorr, t_forcerec *fr)
         }
         else if (EVDW_PME(fr->vdwtype))
         {
-            if (EVDW_SWITCHED(fr->vdwtype) && fr->rvdw_switch == 0)
-            {
-                gmx_fatal(FARGS,
-                          "With dispersion correction rvdw-switch can not be zero "
-                          "for vdw-type = %s", evdw_names[fr->vdwtype]);
-            }
-
             scale  = fr->nblists[0].table_vdw.scale;
             vdwtab = fr->nblists[0].table_vdw.data;
 
@@ -2251,15 +2244,6 @@ void calc_enervirdiff(FILE *fplog, int eDispCorr, t_forcerec *fr)
              */
             fr->enershiftsix = pow(fr->ewaldcoeff_lj, 6) / 6.0;
 
-            /* Calculate C12 values as without PME. */
-            if (EVDW_SWITCHED(fr->vdwtype))
-            {
-                enersum = 0;
-                virsum  = 0;
-                integrate_table(vdwtab, scale, 4, ri0, ri1, &enersum, &virsum);
-                eners[1] -= enersum;
-                virs[1]  -= virsum;
-            }
             /* Add analytical corrections, C6 for the whole range, C12
              * from rvdw_switch to infinity.
              */

diff --git a/src/programs/mdrun/runner.c b/src/programs/mdrun/runner.c
index bdebcd9b139c62a83b2bfd3435b2b6ef970e0dd2..a6b838578a089d58893ed45a499b7f5359db669e 100644 (file)
--- a/src/programs/mdrun/runner.c
+++ b/src/programs/mdrun/runner.c
@@ -77,6 +77,7 @@
 #include "membed.h"
 #include "macros.h"
 #include "gmx_thread_affinity.h"
+#include "inputrec.h"
 
 #include "gromacs/fileio/tpxio.h"
 #include "gromacs/mdlib/nbnxn_search.h"
@@ -762,15 +763,15 @@ static void convert_to_verlet_scheme(FILE *fplog,
     {
         gmx_fatal(FARGS, "User non-bonded potentials are not (yet) supported with the Verlet scheme");
     }
-    else if (EVDW_SWITCHED(ir->vdwtype) || EEL_SWITCHED(ir->coulombtype))
+    else if (ir_vdw_switched(ir) || ir_coulomb_switched(ir))
     {
         md_print_warn(NULL, fplog, "Converting switched or shifted interactions to a shifted potential (without force shift), this will lead to slightly different interaction potentials");
 
-        if (EVDW_SWITCHED(ir->vdwtype))
+        if (ir_vdw_switched(ir))
         {
             ir->vdwtype = evdwCUT;
         }
-        if (EEL_SWITCHED(ir->coulombtype))
+        if (ir_coulomb_switched(ir))
         {
             if (EEL_FULL(ir->coulombtype))
             {